diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Rebin.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Rebin.h index 57386d01a4a2..87b5924515bf 100644 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Rebin.h +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Rebin.h @@ -83,6 +83,11 @@ class DLLExport Rebin : public API::Algorithm void propagateMasks(API::MatrixWorkspace_const_sptr inputW, API::MatrixWorkspace_sptr outputW, int hist); + +private: + // method to check if removing background is requested and possible + API::MatrixWorkspace_const_sptr checkRemoveBackgroundParameters(const API::MatrixWorkspace_sptr &inputWS,int &eMode, bool PreserveEvents); + }; } // namespace Algorithms diff --git a/Code/Mantid/Framework/Algorithms/src/Rebin.cpp b/Code/Mantid/Framework/Algorithms/src/Rebin.cpp index cae7feb175d2..ba3d4315f27f 100644 --- a/Code/Mantid/Framework/Algorithms/src/Rebin.cpp +++ b/Code/Mantid/Framework/Algorithms/src/Rebin.cpp @@ -21,7 +21,7 @@ namespace Mantid // Register the class into the algorithm factory DECLARE_ALGORITHM(Rebin) - + using namespace Kernel; using namespace API; using DataObjects::EventList; @@ -34,15 +34,15 @@ namespace Mantid //--------------------------------------------------------------------------------------------- /** - * Return the rebin parameters from a user input - * @param inParams Input vector from user - * @param inputWS Input workspace from user - * @param logger A reference to a logger - * @returns A new vector containing the rebin parameters - */ + * Return the rebin parameters from a user input + * @param inParams Input vector from user + * @param inputWS Input workspace from user + * @param logger A reference to a logger + * @returns A new vector containing the rebin parameters + */ std::vector Rebin::rebinParamsFromInput(const std::vector & inParams, - const API::MatrixWorkspace & inputWS, - Kernel::Logger & logger) + const API::MatrixWorkspace & inputWS, + Kernel::Logger & logger) { std::vector rbParams; // The validator only passes parameters with size 1, or 3xn. No need to check again here @@ -92,9 +92,9 @@ namespace Mantid "non-event Workspace, respectively. Negative width values indicate logarithmic binning. "); declareProperty("PreserveEvents", true,"Keep the output workspace as an EventWorkspace, " - "if the input has events. If the input and output EventWorkspace " - "names are the same, only the X bins are set, which is very quick. If false, " - "then the workspace gets converted to a Workspace2D histogram."); + "if the input has events. If the input and output EventWorkspace " + "names are the same, only the X bins are set, which is very quick. If false, " + "then the workspace gets converted to a Workspace2D histogram."); declareProperty("FullBinsOnly", false, "Omit the final bin if it's width is smaller than the step size"); //--------------------------------------------------------------- @@ -104,17 +104,18 @@ namespace Mantid vsValidator->add("TOF"); vsValidator->add(); declareProperty(new WorkspaceProperty<>("FlatBkgWorkspace","",Direction::Input,API::PropertyMode::Optional,vsValidator), - "An optional histogram workspace in the units of TOF, containing the same number of spectra as the \"InputWorkspace\" " - "and single Y value per each spectra, representing flat background in this time bin. " + "An optional histogram workspace in the units of TOF defined background for removal during rebinning." + "The workspace has to have single value or contain the same number of spectra as the \"InputWorkspace\" and single Y value per each spectra, + "representing flat background in this time bin. " "If such workspace is present, the value of the flat background provided by this workspace is removed " - "from each spectra of the rebinned workspace. This works for histogram and event workspace when events are not retained - "but actually useful for event workspace in the units different from TOF."); + "from each spectra of the rebinned workspace. This works for histogram and event workspace when events are not retained " + "but actually useful mainly for removing background from event workspace in the units different from TOF."); std::vector dE_modes = Kernel::DeltaEMode().availableTypes(); declareProperty("dEAnalysisMode",dE_modes[Kernel::DeltaEMode::Direct],boost::make_shared(dE_modes), "If FlatBkgWorkspace, this property is used to define the units conversion from TOF to the units of the InputWorkspace",Direction::Input); setPropertySettings("dEAnalysisMode", - new Kernel::VisibleWhenProperty("FlatBkgWorkspace", IS_NOT_EQUAL_TO, "")); + new Kernel::VisibleWhenProperty("FlatBkgWorkspace", IS_NOT_EQUAL_TO, "")); } @@ -143,16 +144,17 @@ namespace Mantid // workspace independent determination of length const int histnumber = static_cast(inputWS->getNumberHistograms()); - + int eMode; // in convert units emode is still integer - const bool remove_background = checkRemoveBaclgroundParameters(inputWS,eMode); + API::MatrixWorkspace_const_sptr bkgWS = checkRemoveBackgroundParameters(inputWS,eMode,PreserveEvents); + const bool remove_background = bool(bkgWS); bool fullBinsOnly = getProperty("FullBinsOnly"); MantidVecPtr XValues_new; // create new output X axis const int ntcnew = VectorHelper::createAxisFromRebinParams(rbParams, XValues_new.access(), - true, fullBinsOnly); + true, fullBinsOnly); //--------------------------------------------------------------------------------- //Now, determine if the input workspace is actually an EventWorkspace @@ -177,7 +179,7 @@ namespace Mantid //Make a brand new EventWorkspace eventOutputWS = boost::dynamic_pointer_cast( - API::WorkspaceFactory::Instance().create("EventWorkspace", inputWS->getNumberHistograms(), 2, 1)); + API::WorkspaceFactory::Instance().create("EventWorkspace", inputWS->getNumberHistograms(), 2, 1)); //Copy geometry over. API::WorkspaceFactory::Instance().initializeFromParent(inputWS, eventOutputWS, false); //You need to copy over the data as well. @@ -205,44 +207,44 @@ namespace Mantid //Go through all the histograms and set the data PARALLEL_FOR3(inputWS, eventInputWS, outputWS) - for (int i=0; i < histnumber; ++i) - { - PARALLEL_START_INTERUPT_REGION - - //Set the X axis for each output histogram - outputWS->setX(i, XValues_new); - - //Get a const event list reference. eventInputWS->dataY() doesn't work. - const EventList& el = eventInputWS->getEventList(i); - MantidVec y_data, e_data; - // The EventList takes care of histogramming. - el.generateHistogram(*XValues_new, y_data, e_data); - - //Copy the data over. - outputWS->dataY(i).assign(y_data.begin(), y_data.end()); - outputWS->dataE(i).assign(e_data.begin(), e_data.end()); - - //Report progress - prog.report(name()); - PARALLEL_END_INTERUPT_REGION - } - PARALLEL_CHECK_INTERUPT_REGION - - //Copy all the axes - for (int i=1; iaxes(); i++) - { - outputWS->replaceAxis( i, inputWS->getAxis(i)->clone(outputWS.get()) ); - outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit(); - } - - //Copy the units over too. - for (int i=0; i < outputWS->axes(); ++i) - outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit(); - outputWS->setYUnit(eventInputWS->YUnit()); - outputWS->setYUnitLabel(eventInputWS->YUnitLabel()); - - // Assign it to the output workspace property - setProperty("OutputWorkspace", outputWS); + for (int i=0; i < histnumber; ++i) + { + PARALLEL_START_INTERUPT_REGION + + //Set the X axis for each output histogram + outputWS->setX(i, XValues_new); + + //Get a const event list reference. eventInputWS->dataY() doesn't work. + const EventList& el = eventInputWS->getEventList(i); + MantidVec y_data, e_data; + // The EventList takes care of histogramming. + el.generateHistogram(*XValues_new, y_data, e_data); + + //Copy the data over. + outputWS->dataY(i).assign(y_data.begin(), y_data.end()); + outputWS->dataE(i).assign(e_data.begin(), e_data.end()); + + //Report progress + prog.report(name()); + PARALLEL_END_INTERUPT_REGION + } + PARALLEL_CHECK_INTERUPT_REGION + + //Copy all the axes + for (int i=1; iaxes(); i++) + { + outputWS->replaceAxis( i, inputWS->getAxis(i)->clone(outputWS.get()) ); + outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit(); + } + + //Copy the units over too. + for (int i=0; i < outputWS->axes(); ++i) + outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit(); + outputWS->setYUnit(eventInputWS->YUnit()); + outputWS->setYUnitLabel(eventInputWS->YUnitLabel()); + + // Assign it to the output workspace property + setProperty("OutputWorkspace", outputWS); } } // END ---- EventWorkspace @@ -273,92 +275,141 @@ namespace Mantid Progress prog(this,0.0,1.0,histnumber); PARALLEL_FOR2(inputWS,outputWS) - for (int hist=0; hist < histnumber;++hist) - { - PARALLEL_START_INTERUPT_REGION - // get const references to input Workspace arrays (no copying) - const MantidVec& XValues = inputWS->readX(hist); - const MantidVec& YValues = inputWS->readY(hist); - const MantidVec& YErrors = inputWS->readE(hist); - - //get references to output workspace data (no copying) - MantidVec& YValues_new=outputWS->dataY(hist); - MantidVec& YErrors_new=outputWS->dataE(hist); - - // output data arrays are implicitly filled by function - try { - VectorHelper::rebin(XValues,YValues,YErrors,*XValues_new,YValues_new,YErrors_new, dist); - } catch (std::exception& ex) + for (int hist=0; hist < histnumber;++hist) { - g_log.error() << "Error in rebin function: " << ex.what() << std::endl; - throw; - } + PARALLEL_START_INTERUPT_REGION + // get const references to input Workspace arrays (no copying) + const MantidVec& XValues = inputWS->readX(hist); + const MantidVec& YValues = inputWS->readY(hist); + const MantidVec& YErrors = inputWS->readE(hist); + + //get references to output workspace data (no copying) + MantidVec& YValues_new=outputWS->dataY(hist); + MantidVec& YErrors_new=outputWS->dataE(hist); + + // output data arrays are implicitly filled by function + try { + VectorHelper::rebin(XValues,YValues,YErrors,*XValues_new,YValues_new,YErrors_new, dist); + } catch (std::exception& ex) + { + g_log.error() << "Error in rebin function: " << ex.what() << std::endl; + throw; + } + + // Populate the output workspace X values + outputWS->setX(hist,XValues_new); - // Populate the output workspace X values - outputWS->setX(hist,XValues_new); + prog.report(name()); + PARALLEL_END_INTERUPT_REGION + } + PARALLEL_CHECK_INTERUPT_REGION + outputWS->isDistribution(dist); - prog.report(name()); - PARALLEL_END_INTERUPT_REGION - } - PARALLEL_CHECK_INTERUPT_REGION - outputWS->isDistribution(dist); + // Now propagate any masking correctly to the output workspace + // More efficient to have this in a separate loop because + // MatrixWorkspace::maskBins blocks multi-threading + for (int i=0; i < histnumber; ++i) + { + if ( inputWS->hasMaskedBins(i) ) // Does the current spectrum have any masked bins? + { + this->propagateMasks(inputWS,outputWS,i); + } + } + //Copy the units over too. + for (int i=0; i < outputWS->axes(); ++i) + { + outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit(); + } - // Now propagate any masking correctly to the output workspace - // More efficient to have this in a separate loop because - // MatrixWorkspace::maskBins blocks multi-threading - for (int i=0; i < histnumber; ++i) - { - if ( inputWS->hasMaskedBins(i) ) // Does the current spectrum have any masked bins? + if ( ! isHist ) { - this->propagateMasks(inputWS,outputWS,i); + g_log.information() << "Rebin: Converting Data back to Data Points." << std::endl; + Mantid::API::Algorithm_sptr ChildAlg = createChildAlgorithm("ConvertToPointData"); + ChildAlg->initialize(); + ChildAlg->setProperty("InputWorkspace", outputWS); + ChildAlg->execute(); + outputWS = ChildAlg->getProperty("OutputWorkspace"); } - } - //Copy the units over too. - for (int i=0; i < outputWS->axes(); ++i) - { - outputWS->getAxis(i)->unit() = inputWS->getAxis(i)->unit(); - } - if ( ! isHist ) + // Assign it to the output workspace property + setProperty("OutputWorkspace",outputWS); + + + } // END ---- Workspace2D + + return; + } + /*** method checks if removing background is requested and possible + and if yes sets up class parameter referring the workspace with flat background + + @param inputWS the input workspace where removing background is expected + @param eMode the energy conversion mode, defined/requested for the instrument + + @returns background workspace if removal is expected + */ + API::MatrixWorkspace_const_sptr Rebin::checkRemoveBackgroundParameters(const API::MatrixWorkspace_sptr &inputWS,int &eMode, bool PreserveEvents) + { + API::MatrixWorkspace_const_sptr bkgWksp = getProperty("FlatBkgWorkspace"); + + if (!bkgWksp) + return bkgWksp; + + // can not remove background while preserving events + if(PreserveEvents) + return API::MatrixWorkspace_const_sptr(); + + // source workspace has to have full instrument defined to perform background removal using this procedure. + auto pInstrument = inputWS->getInstrument(); + if(pInstrument) + { + if(!pInstrument->getSample()) { - g_log.information() << "Rebin: Converting Data back to Data Points." << std::endl; - Mantid::API::Algorithm_sptr ChildAlg = createChildAlgorithm("ConvertToPointData"); - ChildAlg->initialize(); - ChildAlg->setProperty("InputWorkspace", outputWS); - ChildAlg->execute(); - outputWS = ChildAlg->getProperty("OutputWorkspace"); + g_log.warning()<<" Workspace: "<< inputWS->getName()<< " to rebin does not have properly defined instrument. Background removal will not be performed"; + return API::MatrixWorkspace_const_sptr(); } - // Assign it to the output workspace property - setProperty("OutputWorkspace",outputWS); + } + else + { + g_log.warning()<<" Workspace: "<< inputWS->getName()<< " to rebin does not have instrument attached. Background removal will not be performed"; + return API::MatrixWorkspace_const_sptr(); + } + + if(!(bkgWksp->getNumberHistograms() == 1 || inputWS->getNumberHistograms()==bkgWksp->getNumberHistograms())) + { + g_log.warning()<<" Background Workspace: "<< bkgWksp->getName()<< " should have the same number of spectra as source workspace or be a single histogram workspace\n"; + return API::MatrixWorkspace_const_sptr(); + } - } // END ---- Workspace2D + const std::string emodeStr = getProperty("dEAnalysisMode"); + eMode = static_cast(Kernel::DeltaEMode().fromString(emodeStr)); - return; + return bkgWksp; } -// -// /** Continue execution for EventWorkspace scenario */ -// void Rebin::execEvent() -// { -// // retrieve the properties -// std::vector rb_params=getProperty("Params"); -// -// } + + // + // /** Continue execution for EventWorkspace scenario */ + // void Rebin::execEvent() + // { + // // retrieve the properties + // std::vector rb_params=getProperty("Params"); + // + // } /** Takes the masks in the input workspace and apportions the weights into the new bins that overlap - * with a masked bin. These bins are then masked with the calculated weight. - * - * @param inputWS :: The input workspace - * @param outputWS :: The output workspace - * @param hist :: The index of the current histogram - */ + * with a masked bin. These bins are then masked with the calculated weight. + * + * @param inputWS :: The input workspace + * @param outputWS :: The output workspace + * @param hist :: The index of the current histogram + */ void Rebin::propagateMasks(API::MatrixWorkspace_const_sptr inputWS, API::MatrixWorkspace_sptr outputWS, int hist) { // Not too happy with the efficiency of this way of doing it, but it's a lot simpler to use the // existing rebin algorithm to distribute the weights than to re-implement it for this - + MantidVec masked_bins,weights; // Get a reference to the list of masked bins for this spectrum const MatrixWorkspace::MaskList& mask = inputWS->maskedBins(hist); @@ -380,7 +431,7 @@ namespace Mantid weights.push_back((*it).second); masked_bins.push_back(XValues[(*it).first + 1]); } - + // Create a zero vector for the errors because we don't care about them here const MantidVec zeroes(weights.size(),0.0); // Create a vector to hold the redistributed weights @@ -388,13 +439,13 @@ namespace Mantid MantidVec newWeights(XValues_new.size()-1),zeroes2(XValues_new.size()-1); // Use rebin function to redistribute the weights. Note that distribution flag is set VectorHelper::rebin(masked_bins,weights,zeroes,XValues_new,newWeights,zeroes2,true); - + // Now process the output vector and fill the new masking list for (size_t index = 0; index < newWeights.size(); ++index) { if ( newWeights[index] > 0.0 ) outputWS->flagMasked(hist, index, newWeights[index]); } } - + } // namespace Algorithm } // namespace Mantid