From c774e4c9def3252796c74963b6bfa9e06d14d877 Mon Sep 17 00:00:00 2001 From: Wenduo Zhou Date: Fri, 17 Jan 2014 11:31:54 -0500 Subject: [PATCH] Enabled to split events using splitter matrix workspace. Refs #8685. --- .../inc/MantidAlgorithms/FilterEvents.h | 30 +- .../Framework/Algorithms/src/FilterEvents.cpp | 288 ++++++++++++++---- 2 files changed, 255 insertions(+), 63 deletions(-) diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FilterEvents.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FilterEvents.h index 25782f566da3..9b3b50fa1253 100644 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FilterEvents.h +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/FilterEvents.h @@ -61,18 +61,31 @@ namespace Algorithms // Implement abstract Algorithm methods void exec(); + /// Process user input properties + void processProperties(); + void processSplittersWorkspace(); - void createOutputWorkspaces(std::string outputwsnamebase); + /// + void processMatrixSplitterWorkspace(); + + void createOutputWorkspaces(); void importDetectorTOFCalibration(); void filterEventsBySplitters(double progressamount); + /// + void filterEventsByVectorSplitters(double progressamount); + DataObjects::EventWorkspace_sptr m_eventWS; - DataObjects::SplittersWorkspace_sptr mSplittersWorkspace; + DataObjects::SplittersWorkspace_sptr m_splittersWorkspace; + API::MatrixWorkspace_const_sptr m_matrixSplitterWS; DataObjects::TableWorkspace_sptr m_detCorrectWorkspace; + /// Flag to use matrix splitters or table splitters + bool m_useTableSplitters; + std::set m_workGroupIndexes; Kernel::TimeSplitterType m_splitters; std::map m_outputWS; @@ -83,7 +96,7 @@ namespace Algorithms bool mFilterByPulseTime; - DataObjects::TableWorkspace_sptr mInformationWS; + DataObjects::TableWorkspace_sptr m_informationWS; bool mWithInfo; double mProgress; @@ -99,6 +112,17 @@ namespace Algorithms /// Flag to generate TOF correction bool m_genTOFCorrection; + /// Base of output workspace's name + std::string m_outputWSNameBase; + + /// + bool m_toGroupWS; + + /// + std::vector m_vecSplitterTime; + /// + std::vector m_vecSplitterGroup; + }; diff --git a/Code/Mantid/Framework/Algorithms/src/FilterEvents.cpp b/Code/Mantid/Framework/Algorithms/src/FilterEvents.cpp index 542405cb33f2..f0c50359094e 100644 --- a/Code/Mantid/Framework/Algorithms/src/FilterEvents.cpp +++ b/Code/Mantid/Framework/Algorithms/src/FilterEvents.cpp @@ -83,17 +83,17 @@ namespace Algorithms void FilterEvents::init() { declareProperty( - new API::WorkspaceProperty("InputWorkspace","",Direction::Input), + new API::WorkspaceProperty("InputWorkspace","",Direction::Input), "An input event workspace" ); declareProperty("OutputWorkspaceBaseName", "OutputWorkspace", "The base name to use for the output workspace" ); - declareProperty(new WorkspaceProperty("InformationWorkspace", "", Direction::Input, PropertyMode::Optional), + declareProperty(new WorkspaceProperty("InformationWorkspace", "", Direction::Input, PropertyMode::Optional), "Optional output for the information of each splitter workspace index."); declareProperty( - new API::WorkspaceProperty("SplitterWorkspace", "", Direction::Input), + new API::WorkspaceProperty("SplitterWorkspace", "", Direction::Input), "An input SpilltersWorskpace for filtering"); auto tablewsprop = new WorkspaceProperty("DetectorTOFCorrectionWorkspace", "", Direction::Input, PropertyMode::Optional); @@ -121,58 +121,20 @@ namespace Algorithms void FilterEvents::exec() { // Process algorithm properties - m_eventWS = this->getProperty("InputWorkspace"); - if (!m_eventWS) - { - stringstream errss; - errss << "Inputworkspace is not event workspace. "; - g_log.error(errss.str()); - throw std::invalid_argument(errss.str()); - } - - mSplittersWorkspace = this->getProperty("SplitterWorkspace"); - mInformationWS = this->getProperty("InformationWorkspace"); - - std::string outputwsnamebase = this->getProperty("OutputWorkspaceBaseName"); - m_detCorrectWorkspace = getProperty("DetectorTOFCorrectionWorkspace"); - mFilterByPulseTime = this->getProperty("FilterByPulseTime"); - - bool togroupws = this->getProperty("GroupWorkspaces"); - - // Do correction or not? - m_doTOFCorrection = true; - m_genTOFCorrection = getProperty("GenerateTOFCorrection"); - if (m_detCorrectWorkspace) - { - // User specify detector TOF correction, then no need to generate TOF correction - m_genTOFCorrection = false; - } - else if (m_genTOFCorrection) - { - // If no detector TOF correction workspace is specified but specified to go generate TOF - m_doTOFCorrection = true; - } - else - { - // No correction is needed - m_doTOFCorrection = false; - } - - // Informatin workspace is specified? - if (!mInformationWS) - mWithInfo = false; - else - mWithInfo = true; + processProperties(); // Parse splitters mProgress = 0.0; progress(mProgress, "Processing SplittersWorkspace."); - processSplittersWorkspace(); + if (m_useTableSplitters) + processSplittersWorkspace(); + else + processMatrixSplitterWorkspace(); // Create output workspaces mProgress = 0.1; progress(mProgress, "Create Output Workspaces."); - createOutputWorkspaces(outputwsnamebase); + createOutputWorkspaces(); // Optionall import corrections mProgress = 0.20; @@ -184,19 +146,22 @@ namespace Algorithms mProgress = 0.30; progress(mProgress, "Filter Events."); double progressamount; - if (togroupws) + if (m_toGroupWS) progressamount = 0.6; else progressamount = 0.7; - filterEventsBySplitters(progressamount); + if (m_useTableSplitters) + filterEventsBySplitters(progressamount); + else + filterEventsByVectorSplitters(progressamount); // Optional to group detector - if (togroupws) + if (m_toGroupWS) { mProgress = 0.9; progress(mProgress, "Group workspaces"); - std::string groupname = outputwsnamebase; + std::string groupname = m_outputWSNameBase; API::IAlgorithm_sptr groupws = createChildAlgorithm("GroupWorkspaces", 0.99, 1.00, true); // groupws->initialize(); groupws->setAlwaysStoreInADS(true); @@ -215,6 +180,73 @@ namespace Algorithms return; } + + void FilterEvents::processProperties() + { + m_eventWS = this->getProperty("InputWorkspace"); + if (!m_eventWS) + { + stringstream errss; + errss << "Inputworkspace is not event workspace. "; + g_log.error(errss.str()); + throw std::invalid_argument(errss.str()); + } + + // Process splitting workspace (table or data) + API::Workspace_sptr tempws = this->getProperty("SplitterWorkspace"); + + m_splittersWorkspace = boost::dynamic_pointer_cast(tempws); + if (m_splittersWorkspace) + { + m_useTableSplitters = true; + } + else + { + m_matrixSplitterWS = boost::dynamic_pointer_cast(tempws); + if (m_matrixSplitterWS) + { + m_useTableSplitters = false; + } + else + { + throw runtime_error("Invalid type of input workspace, neither SplittersWorkspace nor MatrixWorkspace."); + } + } + + m_informationWS = this->getProperty("InformationWorkspace"); + + m_outputWSNameBase = this->getPropertyValue("OutputWorkspaceBaseName"); + m_detCorrectWorkspace = getProperty("DetectorTOFCorrectionWorkspace"); + mFilterByPulseTime = this->getProperty("FilterByPulseTime"); + + m_toGroupWS = this->getProperty("GroupWorkspaces"); + + // Do correction or not? + m_doTOFCorrection = true; + m_genTOFCorrection = getProperty("GenerateTOFCorrection"); + if (m_detCorrectWorkspace) + { + // User specify detector TOF correction, then no need to generate TOF correction + m_genTOFCorrection = false; + } + else if (m_genTOFCorrection) + { + // If no detector TOF correction workspace is specified but specified to go generate TOF + m_doTOFCorrection = true; + } + else + { + // No correction is needed + m_doTOFCorrection = false; + } + + // Informatin workspace is specified? + if (!m_informationWS) + mWithInfo = false; + else + mWithInfo = true; + } + //---------------------------------------------------------------------------------------------- /** Convert SplitterWorkspace object to TimeSplitterType (sorted vector) * and create a map for all workspace group number @@ -222,14 +254,14 @@ namespace Algorithms void FilterEvents::processSplittersWorkspace() { // 1. Init data structure - size_t numsplitters = mSplittersWorkspace->getNumberSplitters(); + size_t numsplitters = m_splittersWorkspace->getNumberSplitters(); m_splitters.reserve(numsplitters); // 2. Insert all splitters bool inorder = true; for (size_t i = 0; i < numsplitters; i ++) { - m_splitters.push_back(mSplittersWorkspace->getSplitter(i)); + m_splitters.push_back(m_splittersWorkspace->getSplitter(i)); m_workGroupIndexes.insert(m_splitters.back().index()); if (inorder && i > 0 && m_splitters[i] < m_splitters[i-1]) @@ -250,10 +282,10 @@ namespace Algorithms // 5. Add information if (mWithInfo) { - if (m_workGroupIndexes.size() > mInformationWS->rowCount()+1) + if (m_workGroupIndexes.size() > m_informationWS->rowCount()+1) { g_log.warning() << "Input Splitters Workspace has different entries (" << m_workGroupIndexes.size() -1 - << ") than input information workspaces (" << mInformationWS->rowCount() << "). " + << ") than input information workspaces (" << m_informationWS->rowCount() << "). " << " Information may not be accurate. " << std::endl; } } @@ -261,18 +293,49 @@ namespace Algorithms return; } + + void FilterEvents::processMatrixSplitterWorkspace() + { + // Check input workspace validity + const MantidVec& vecX = m_matrixSplitterWS->readX(0); + const MantidVec& vecY = m_matrixSplitterWS->readY(0); + size_t sizex = vecX.size(); + size_t sizey = vecY.size(); + if (sizex - sizey != 1) + throw runtime_error("Size must be N and N-1."); + + // Assign vectors for time comparison + m_vecSplitterTime.assign(vecX.size(), 0); + m_vecSplitterGroup.assign(vecY.size(), -1); + + // Transform vector + for (size_t i = 0; i < sizex; ++i) + { + m_vecSplitterTime[i] = static_cast(vecX[i]); + } + for (size_t i = 0; i < sizey; ++i) + { + m_vecSplitterGroup[i] = static_cast(vecY[i]); + m_workGroupIndexes.insert(m_vecSplitterGroup[i]); + } + + return; + } + + //---------------------------------------------------------------------------------------------- /** Create a list of EventWorkspace for output */ - void FilterEvents::createOutputWorkspaces(std::string outputwsnamebase) + void FilterEvents::createOutputWorkspaces() { + // Convert information workspace to map std::map infomap; if (mWithInfo) { - for (size_t ir = 0; ir < mInformationWS->rowCount(); ++ ir) + for (size_t ir = 0; ir < m_informationWS->rowCount(); ++ ir) { - API::TableRow row = mInformationWS->getRow(ir); + API::TableRow row = m_informationWS->getRow(ir); int& indexws = row.Int(0); std::string& info = row.String(1); infomap.insert(std::make_pair(indexws, info)); @@ -310,11 +373,11 @@ namespace Algorithms std::stringstream wsname; if (wsgroup >= 0) { - wsname << outputwsnamebase << "_" << (wsgroup+delta_wsindex); + wsname << m_outputWSNameBase << "_" << (wsgroup+delta_wsindex); } else { - wsname << outputwsnamebase << "_unfiltered"; + wsname << m_outputWSNameBase << "_unfiltered"; if (from1) add2output = false; } @@ -591,6 +654,111 @@ namespace Algorithms } + //---------------------------------------------------------------------------------------------- + /** Split events by splitters represented by vector + */ + void FilterEvents::filterEventsByVectorSplitters(double progressamount) + { + size_t numberOfSpectra = m_eventWS->getNumberHistograms(); + // FIXME : consider to use vector to index workspace and event list + std::map::iterator wsiter; + + // Loop over the histograms (detector spectra) to do split from 1 event list to N event list + g_log.debug() << "Number of spectra in input/source EventWorkspace = " << numberOfSpectra << ".\n"; + + // FIXME Make it parallel + // PARALLEL_FOR_NO_WSP_CHECK() + for (int64_t iws = 0; iws < int64_t(numberOfSpectra); ++iws) + { + // FIXME Make it parallel + // PARALLEL_START_INTERUPT_REGION + + // Get the output event lists (should be empty) to be a map + map outputs; + for (wsiter = m_outputWS.begin(); wsiter != m_outputWS.end(); ++ wsiter) + { + int index = wsiter->first; + DataObjects::EventList* output_el = wsiter->second->getEventListPtr(iws); + outputs.insert(std::make_pair(index, output_el)); + } + + // Get a holder on input workspace's event list of this spectrum + const DataObjects::EventList& input_el = m_eventWS->getEventList(iws); + + // Perform the filtering (using the splitting function and just one output) + if (mFilterByPulseTime) + { + throw runtime_error("It is not a good practice to split fast event by pulse time. "); + // input_el.splitByPulseTimeMatrixSplitter(m_vecSplitterTime, m_vecSplitterGroup, outputs); + } + else if (m_doTOFCorrection) + { + input_el.splitByFullTimeMatrixSplitter(m_vecSplitterTime, m_vecSplitterGroup, outputs, + m_detTofOffsets[iws], m_doTOFCorrection); + } + else + { + input_el.splitByFullTimeMatrixSplitter(m_vecSplitterTime, m_vecSplitterGroup, outputs, 1.0, m_doTOFCorrection); + } + + mProgress = 0.3+(progressamount-0.2)*static_cast(iws)/static_cast(numberOfSpectra); + progress(mProgress, "Filtering events"); + + // FIXME - Turn on parallel + // PARALLEL_END_INTERUPT_REGION + } // END FOR i = 0 + // PARALLEL_CHECK_INTERUPT_REGION + // FIXME - Turn on parallel + + + // Finish (1) adding events and splitting the sample logs in each target workspace. + progress(0.1+progressamount, "Splitting logs"); + + std::vector lognames; + this->getTimeSeriesLogNames(lognames); + g_log.debug() << "[FilterEvents D1214]: Number of TimeSeries Logs = " << lognames.size() + << " to " << m_outputWS.size() << " outptu workspaces. \n"; + + double numws = static_cast(m_outputWS.size()); + double outwsindex = 0.; + for (wsiter = m_outputWS.begin(); wsiter != m_outputWS.end(); ++wsiter) + { + int wsindex = wsiter->first; + DataObjects::EventWorkspace_sptr opws = wsiter->second; + + // Generate a list of splitters for current output workspace + Kernel::TimeSplitterType splitters; + generateSplitters(wsindex, splitters); + + g_log.debug() << "[FilterEvents D1215]: Output orkspace Index " << wsindex + << ": Name = " << opws->name() << "; Number of splitters = " << splitters.size() << ".\n"; + + // Skip output workspace has ZERO splitters + if (splitters.size() == 0) + { + g_log.warning() << "[FilterEvents] Workspace " << opws->name() << " Indexed @ " << wsindex << + " won't have logs splitted due to zero splitter size. " << ".\n"; + continue; + } + + // Split log + size_t numlogs = lognames.size(); + for (size_t ilog = 0; ilog < numlogs; ++ilog) + { + this->splitLog(opws, lognames[ilog], splitters); + } + opws->mutableRun().integrateProtonCharge(); + + progress(0.1+progressamount+outwsindex/numws*0.2, "Splitting logs"); + outwsindex += 1.; + } + + return; + + return; + } + + //---------------------------------------------------------------------------------------------- /** Generate splitters for specified workspace index as a subset of m_splitters */