diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/RemoveBackground.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/RemoveBackground.h new file mode 100644 index 000000000000..03139d781566 --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/RemoveBackground.h @@ -0,0 +1,122 @@ +#ifndef MANTID_ALGORITHMS_REBIN_H_ +#define MANTID_ALGORITHMS_REBIN_H_ + +//---------------------------------------------------------------------- +// Includes +//---------------------------------------------------------------------- +#include "MantidAPI/Algorithm.h" +#include "MantidAlgorithms/BackgroundHelper.h" + +namespace Mantid +{ + namespace Algorithms + { + /** Performs removal of constant (and possibly non-constant after simple modification) background calculated in TOF units + from a matrix workspace, expressed in units, different from TOF. + + @date 26/10/2014 + + + Copyright © 2008-9 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory + + This file is part of Mantid. + + Mantid is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. + + Mantid is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + + File change history is stored at: + Code Documentation is available at: + */ + + class DLLExport RemoveBackground : public API::Algorithm + { + public: + /// Default constructor + RemoveBackground() : API::Algorithm(),m_BackgroundHelper() {}; + /// Destructor + virtual ~RemoveBackground(){}; + /// Algorithm's name for identification overriding a virtual method + virtual const std::string name() const { return "RemoveBackground";} + ///Summary of algorithms purpose + virtual const std::string summary() const {return "Removes background (constant for now) calculated in TOF units " + "from a matrix workspace, expressed in units, different from TOF";} + + /// Algorithm's version for identification overriding a virtual method + virtual int version() const { return 1;} + /// Algorithm's category for identification overriding a virtual method + virtual const std::string category() const { return "Transforms\\RemoveBackground";} + + protected: + + // Overridden Algorithm methods + void init(); + virtual void exec(); + + + 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); + // class responsible for background removal + BackgroundHelper m_BackgroundHelper; + protected: // for testing + /**Class actually performing background removal from a workspace spectra */ + class BackgroundHelper + { + public: + BackgroundHelper(); + ~BackgroundHelper(); + void initialize(const API::MatrixWorkspace_const_sptr &bkgWS,const API::MatrixWorkspace_sptr &sourceWS,int emode,int nTreads=1); + + void removeBackground(int hist,const MantidVec &XValues,MantidVec &y_data,MantidVec &e_data,int tread_num=0)const; + + private: + //vector of pointers to the units conversion class for the working workspace; + std::vector m_WSUnit; + + // shared pointer to the workspace containing background + API::MatrixWorkspace_const_sptr m_bgWs; + // shared pointer to the workspace where background should be removed + API::MatrixWorkspace_const_sptr m_wkWS; + + // if the background workspace is single value workspace + bool m_singleValueBackground; + // average number of counts at background for first spectra of a background workspace + double m_NBg; + // time interval for measuring the background + double m_dtBg; + // Squared error of the background for first spectra of a background workspace + //double m_ErrSq; + // energy conversion mode + int m_Emode; + // source-sample distance + double m_L1; + // incident for direct or analysis for indirect energy for units conversion + double m_Efix; + // shared pointer to the sample + Geometry::IComponent_const_sptr m_Sample; + + // get Ei attached to direct or indirect instrument workspace + double getEi(const API::MatrixWorkspace_const_sptr &inputWS)const; + // the procedure user to delete existing unit converter pointers + void deleteUnitsConverters(); + }; + + + }; + + + + } // namespace Algorithms +} // namespace Mantid + +#endif /*MANTID_ALGORITHMS_REBIN_H_*/ diff --git a/Code/Mantid/Framework/Algorithms/src/RemoveBackground.cpp b/Code/Mantid/Framework/Algorithms/src/RemoveBackground.cpp new file mode 100644 index 000000000000..f5b0e909fb68 --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/src/RemoveBackground.cpp @@ -0,0 +1,553 @@ +//---------------------------------------------------------------------- +// Includes +//---------------------------------------------------------------------- +#include "MantidAlgorithms/RemoveBackground.h" + +#include "MantidKernel/ArrayProperty.h" +#include "MantidKernel/VectorHelper.h" +#include "MantidKernel/ListValidator.h" +#include "MantidKernel/RebinParamsValidator.h" +#include "MantidKernel/VisibleWhenProperty.h" + +#include "MantidAPI/WorkspaceValidators.h" +#include "MantidDataObjects/EventWorkspace.h" +#include "MantidDataObjects/EventList.h" + + +namespace Mantid +{ + namespace Algorithms + { + + // Register the class into the algorithm factory + DECLARE_ALGORITHM(RemoveBackground) + + using namespace Kernel; + using namespace API; + + //--------------------------------------------------------------------------------------------- + // Public methods + //--------------------------------------------------------------------------------------------- + + /** Initialisation method. Declares properties to be used in algorithm. + * + */ + void RemoveBackground::init() + { + auto sourceValidator = boost::make_shared(); + sourceValidator->add(); + sourceValidator->add(); + declareProperty( + new WorkspaceProperty<>("InputWorkspace", "", Direction::Input,sourceValidator), + "Workspace containing the input data"); + declareProperty( + new WorkspaceProperty<>("OutputWorkspace","",Direction::Output), + "The name to give the output workspace"); + + auto vsValidator = boost::make_shared(); + vsValidator->add("TOF"); + vsValidator->add(); + declareProperty(new WorkspaceProperty<>("FlatBkgWorkspace","",Direction::Input,vsValidator), + "An optional histogram workspace in the units of TOF defining 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 the background time region. " + "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 mainly for removing background while rebinning an event workspace in the units different from TOF."); + + std::vector dE_modes = Kernel::DeltaEMode().availableTypes(); + declareProperty("EMode",dE_modes[Kernel::DeltaEMode::Direct],boost::make_shared(dE_modes), + "If FlatBkgWorkspace is present, this property used to define the units for conversion from the units of the InputWorkspace to TOF" ,Direction::Input); + } + + + /** Executes the rebin algorithm + * + * @throw runtime_error Thrown if the bin range does not intersect the range of the input workspace + */ + void RemoveBackground::exec() + { + // Get the input workspace + MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace"); + MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace"); + + // Removing background in-place + bool inPlace = (inputWS == outputWS); + + + // workspace independent determination of length + const int histnumber = static_cast(inputWS->getNumberHistograms()); + + //-- Flat Background removal ---------------------------- + int eMode; // in convert units emode is still integer + API::MatrixWorkspace_const_sptr bkgWS = checkRemoveBackgroundParameters(inputWS,eMode,PreserveEvents); + const bool remove_background = bool(bkgWS); + if(remove_background) + { + int nThreads = omp_get_max_threads(); + m_BackgroundHelper.initialize(bkgWS,inputWS,eMode,nThreads); + } + //------------------------------------------------------- + + bool fullBinsOnly = getProperty("FullBinsOnly"); + + MantidVecPtr XValues_new; + // create new output X axis + const int ntcnew = VectorHelper::createAxisFromRebinParams(rbParams, XValues_new.access(), + true, fullBinsOnly); + + //--------------------------------------------------------------------------------- + //Now, determine if the input workspace is actually an EventWorkspace + EventWorkspace_const_sptr eventInputWS = boost::dynamic_pointer_cast(inputWS); + + if (eventInputWS != NULL) + { + //------- EventWorkspace as input ------------------------------------- + EventWorkspace_sptr eventOutputWS = boost::dynamic_pointer_cast(outputWS); + + if (inPlace && PreserveEvents) + { + // -------------Rebin in-place, preserving events ---------------------------------------------- + // This only sets the X axis. Actual rebinning will be done upon data access. + eventOutputWS->setAllX(XValues_new); + this->setProperty("OutputWorkspace", boost::dynamic_pointer_cast(eventOutputWS)); + } + else if (!inPlace && PreserveEvents) + { + // -------- NOT in-place, but you want to keep events for some reason. ---------------------- + // Must copy the event workspace to a new EventWorkspace (and bin that). + + //Make a brand new EventWorkspace + eventOutputWS = boost::dynamic_pointer_cast( + 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. + eventOutputWS->copyDataFrom( (*eventInputWS) ); + + // This only sets the X axis. Actual rebinning will be done upon data access. + eventOutputWS->setAllX(XValues_new); + + //Cast to the matrixOutputWS and save it + this->setProperty("OutputWorkspace", boost::dynamic_pointer_cast(eventOutputWS)); + } + else + { + //--------- Different output, OR you're inplace but not preserving Events --- create a Workspace2D ------- + g_log.information() << "Creating a Workspace2D from the EventWorkspace " << eventInputWS->getName() << ".\n"; + + //Create a Workspace2D + // This creates a new Workspace2D through a torturous route using the WorkspaceFactory. + // The Workspace2D is created with an EMPTY CONSTRUCTOR + outputWS = WorkspaceFactory::Instance().create("Workspace2D",histnumber,ntcnew,ntcnew-1); + WorkspaceFactory::Instance().initializeFromParent(inputWS, outputWS, true); + + //Initialize progress reporting. + Progress prog(this,0.0,1.0, histnumber); + + //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); + if(remove_background) + { + int id = omp_get_thread_num(); + m_BackgroundHelper.removeBackground(i,*XValues_new,y_data,e_data,id); + } + + + //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 + + else + + { //------- Workspace2D or other MatrixWorkspace --------------------------- + + if ( ! isHist ) + { + g_log.information() << "Rebin: Converting Data to Histogram." << std::endl; + Mantid::API::Algorithm_sptr ChildAlg = createChildAlgorithm("ConvertToHistogram"); + ChildAlg->initialize(); + ChildAlg->setProperty("InputWorkspace", inputWS); + ChildAlg->execute(); + inputWS = ChildAlg->getProperty("OutputWorkspace"); + } + + // This will be the output workspace (exact type may vary) + API::MatrixWorkspace_sptr outputWS; + + // make output Workspace the same type is the input, but with new length of signal array + outputWS = API::WorkspaceFactory::Instance().create(inputWS,histnumber,ntcnew,ntcnew-1); + + + // Copy over the 'vertical' axis + if (inputWS->axes() > 1) outputWS->replaceAxis( 1, inputWS->getAxis(1)->clone(outputWS.get()) ); + + 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); + if(remove_background) + { + int id = omp_get_thread_num(); + m_BackgroundHelper.removeBackground(hist,*XValues_new,YValues_new,YErrors_new,id); + } + } 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); + + 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(); + } + + if ( ! isHist ) + { + 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"); + } + + if(remove_background) + { + const std::list &failedBkgRemoalList= m_BackgroundHelper.getFailingSpectrsList(); + if(!failedBkgRemoalList.empty()) + { + size_t nFailed = failedBkgRemoalList.size(); + if(nFailed == static_cast(histnumber)) + { + g_log.warning()<<" has not been able to remove any background while rebinning workspace "<getName()<< + "\n possible reasons: wrong instrument or units conversion mode\n"; + } + else + { + g_log.debug()<<" has not removed background from "<getInstrument(); + if(pInstrument) + { + if(!pInstrument->getSample()) + { + 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(); + } + + } + 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(); + } + + + const std::string emodeStr = getProperty("EMode"); + eMode = static_cast(Kernel::DeltaEMode().fromString(emodeStr)); + + return bkgWksp; + } +//------------------------------------------------------------------------------------------------------------------------------- + /// Constructor + BackgroundHelper::BackgroundHelper(): + m_WSUnit(), + m_bgWs(),m_wkWS(), + m_singleValueBackground(false), + m_NBg(0),m_dtBg(1), //m_ErrSq(0), + m_Emode(0), + m_L1(0),m_Efix(0), + m_Sample(), + FailingSpectraList() + {}; + BackgroundHelper::~BackgroundHelper() + { + this->deleteUnitsConverters(); + } + + /** The method deletes all units converter allocated*/ + void BackgroundHelper::deleteUnitsConverters() + { + for(size_t i=0;igetAxis(0)->unit()->unitID(); + if(bgUnits!="TOF") + throw std::invalid_argument(" Background Workspace: "+bkgWS->getName()+" should be in the units of TOF"); + + if(!(bkgWS->getNumberHistograms() == 1 || sourceWS->getNumberHistograms()==bkgWS->getNumberHistograms())) + throw std::invalid_argument(" Background Workspace: "+bkgWS->getName()+" should have the same number of spectra as source workspace or be a single histogram workspace"); + + auto WSUnit = sourceWS->getAxis(0)->unit(); + if(!WSUnit) + throw std::invalid_argument(" Source Workspace: "+sourceWS->getName()+" should have units"); + + Geometry::IComponent_const_sptr source = sourceWS->getInstrument()->getSource(); + m_Sample = sourceWS->getInstrument()->getSample(); + if ((!source) || (!m_Sample)) + throw std::invalid_argument("Instrument on Source workspace:"+sourceWS->getName()+"is not sufficiently defined: failed to get source and/or sample"); + m_L1 = source->getDistance(*m_Sample); + + // just in case. + this->deleteUnitsConverters(); + // allocate the array of units converters to avoid units reallocation within a loop + m_WSUnit.assign(nThreads,NULL); + for(int i=0;iclone(); + } + + m_singleValueBackground = false; + if(bkgWS->getNumberHistograms()==0) + m_singleValueBackground = true; + const MantidVec& dataX = bkgWS->dataX(0); + const MantidVec& dataY = bkgWS->dataY(0); + //const MantidVec& dataE = bkgWS->dataE(0); + m_NBg = dataY[0]; + m_dtBg = dataX[1]-dataX[0]; + //m_ErrSq = dataE[0]*dataE[0]; // needs further clarification + + + m_Efix = this->getEi(sourceWS); + } + /**Method removes background from vectors which represent a histogram data for a single spectra + * @param nHist -- number (workspaceID) of the spectra in the workspace, where background going to be removed + * @param XValues -- the spectra x-values (presumably not in TOF units) + * @param y_data -- the spectra signal + * @param e_data -- the spectra errors + * @param threadNum -- number of thread doing conversion (by default 0, single thread) + */ + void BackgroundHelper::removeBackground(int nHist,const MantidVec &XValues,MantidVec &y_data,MantidVec &e_data,int threadNum)const + { + + double dtBg,IBg; + if(m_singleValueBackground) + { + dtBg = m_dtBg; + // ErrSq = m_ErrSq; + IBg = m_NBg; + } + else + { + const MantidVec& dataX = m_bgWs->dataX(nHist); + const MantidVec& dataY = m_bgWs->dataY(nHist); + //const MantidVec& dataE = m_bgWs->dataX(nHist); + dtBg = (dataX[1]-dataX[0]); + IBg = dataY[0]; + //ErrSq= dataE[0]*dataE[0]; // Needs further clarification + } + + try + { + auto detector = m_wkWS->getDetector(nHist); + // + double twoTheta = m_wkWS->detectorTwoTheta(detector); + double L2 = detector->getDistance(*m_Sample); + double delta(std::numeric_limits::quiet_NaN()); + // use thread-specific unit conversion class to avoid multithreading issues + Kernel::Unit * unitConv = m_WSUnit[threadNum]; + unitConv->initialize(m_L1, L2,twoTheta, m_Emode, m_Efix,delta); + double tof1 = unitConv->singleToTOF(XValues[0]); + for(size_t i=0;isingleToTOF(XValues[i+1]); + double Jack = std::fabs((tof2-tof1)/dtBg); + double normBkgrnd = IBg*Jack; + tof1=tof2; + y_data[i] -=normBkgrnd; + //e_data[i] =std::sqrt((ErrSq*Jack*Jack+e_data[i]*e_data[i])/2.); // needs further clarification -- Gaussian error summation? + //--> assume error for background is sqrt(signal): + e_data[i] =std::sqrt((normBkgrnd+e_data[i]*e_data[i])/2.); // needs further clarification -- Gaussian error summation? + } + + } + catch(...) + { + // no background removal for this spectra as it does not have a detector or other reason; How should it be reported? + FailingSpectraList.push_front(nHist); + } + + } + /** Method returns the efixed or Ei value stored in properties of the input workspace. + * Indirect instruments can have eFxed and Direct instruments can have Ei defined as the properties of the workspace. + * + * This method provide guess for efixed for all other kind of instruments. Correct indirect instrument will overwrite + * this value while wrongly defined or different types of instruments will provide the value of "Ei" property (log value) + * or undefined if "Ei" property is not found. + * + */ + double BackgroundHelper::getEi(const API::MatrixWorkspace_const_sptr &inputWS)const + { + double Efi = std::numeric_limits::quiet_NaN(); + + + // is Ei on workspace properties? (it can be defined for some reason if detectors do not have one, and then it would exist as Ei) + bool EiFound(false); + try + { + Efi = inputWS->run().getPropertyValueAsType("Ei"); + EiFound = true; + } + catch(Kernel::Exception::NotFoundError &) + {} + // try to get Efixed as property on a workspace, obtained for indirect instrument + //bool eFixedFound(false); + if (!EiFound) + { + try + { + Efi =inputWS->run().getPropertyValueAsType("eFixed"); + //eFixedFound = true; + } + catch(Kernel::Exception::NotFoundError &) + {} + } + + //if (!(EiFound||eFixedFound)) + // g_log.debug()<<" Ei/eFixed requested but have not been found\n"; + + return Efi; + } + + + + } // namespace Algorithm +} // namespace Mantid diff --git a/Code/Mantid/Framework/Algorithms/test/BackgroundHelperTest.h b/Code/Mantid/Framework/Algorithms/test/RemoveBackgroundTest.h similarity index 100% rename from Code/Mantid/Framework/Algorithms/test/BackgroundHelperTest.h rename to Code/Mantid/Framework/Algorithms/test/RemoveBackgroundTest.h diff --git a/Code/Mantid/docs/source/algorithms/RemoveBackground-v1.rst b/Code/Mantid/docs/source/algorithms/RemoveBackground-v1.rst new file mode 100644 index 000000000000..952817a48b66 --- /dev/null +++ b/Code/Mantid/docs/source/algorithms/RemoveBackground-v1.rst @@ -0,0 +1,411 @@ +.. algorithm:: + +.. summary:: + +.. alias:: + +.. properties:: + +Description +----------- + +The algorithm rebins data with new bin boundaries. The 'params' property +defines new boundaries in intervals :math:`x_i-x_{i+1}\,`. Positive +:math:`\Delta x_i\,` make constant width bins, whilst negative ones +create logarithmic binning using the formula +:math:`x(j+1)=x(j)(1+|\Delta x_i|)\,` + +This algorithms is useful both in data reduction, but also in remapping +`ragged workspaces `__ to a regular set of bin +boundaries. + +Unless the FullBinsOnly option is enabled, the bin immediately before +the specified boundaries :math:`x_2`, :math:`x_3`, ... :math:`x_i` is +likely to have a different width from its neighbours because there can +be no gaps between bins. Rebin ensures that any of these space filling +bins cannot be less than 25% or more than 125% of the width that was +specified. + + +.. _rebin-example-strings: + +Example Rebin param strings +########################### + +* **"-0.0001"**: from min(TOF) to max(TOF) among all events in Logarithmic bins of 0.0001 +* **"*0,100,20000"**: from 0 rebin in constant size bins of 100 up to 20,000 +* **"2,-0.035,10"**: from 10 rebin in Logarithmic bins of 0.035 up to 10 +* **"0,100,10000,200,20000"**: from 0 rebin in steps of 100 to 10,000 then steps of 200 to 20,000 + +For EventWorkspaces +################### + +If the input is an `EventWorkspace `__ and the "Preserve +Events" property is True, the rebinning is performed in place, and only +the X axes of the workspace are set. The actual Y histogram data will +only be requested as needed, for example, when plotting or displaying +the data. + +If "Preserve Events" is false., then the output workspace will be +created as a :ref:`Workspace2D `, with fixed histogram bins, +and all Y data will be computed immediately. All event-specific data is +lost at that point. + +For Data-Point Workspaces +######################### + +If the input workspace contains data points, rather than histograms, +then Rebin will automatically use the +:ref:`ConvertToHistogram ` and +:ref:`ConvertToHistogram ` algorithms before and after +the rebinning has taken place. + +FullBinsOnly option +################### + +If FullBinsOnly option is enabled, each range will only contain bins of +the size equal to the step specified. In other words, the will be no +space filling bins which are bigger or smaller than the other ones. + +This, however, means that specified bin boundaries might get amended in +the process of binning. For example, if rebin *Param* string is +specified as "0, 2, 4.5, 3, 11" and FullBinsOnly is enabled, the +following will happen: + +- From 0 rebin in bins of size 2 **up to 4**. 4.5 is ignored, because + otherwise we would need to create a filling bin of size 0.5. +- **From 4** rebin in bins of size 3 **up to 10**. + +Hence the actual *Param* string used is "0, 2, 4, 3, 10". + + +Remove Background during rebinning options +########################################## + +These options allow you to remove flat background, defined by the a single bin +histogram workspace with X-axis in the units of TOF from a workspace in any +units with known conversion into TOF. The background removal occurs during rebinning + +These options are especially useful during reduction +of event workspaces in multirep mode, where different event regions are associated with +different incident energies and rebinned into appropriate energy range together with +background removal on-the-fly. + +The algorithm used during background removal is equivalent to the proof-of concept one, +presented below, except intermediate workspaces are not created and the background removal calculations +performed during rebinning. +Errors of the background workspace are currently ignored and their value +is calculated as the square root of correspondent background signal:: + + from mantid.simpleapi import * + from mantid import config + import numpy as np + import sys + import os + + + maps_dir = '/home/user/InstrumentFiles/let/' + data_dir ='/home/user/results' + ref_data_dir = '/home/user/SystemTests/AnalysisTests/ReferenceResults' + config.setDataSearchDirs('{0};{1};{2}'.format(data_dir,maps_dir,ref_data_dir)) + config['defaultsave.directory'] = data_dir # folder to save resulting spe/nxspe files. Defaults are in + + # the name of a workspace containing background + filename = 'LET0007438' + groupedFilename = filename+'rings'; + # + Ei= 25 + e_min = -20 + e_max = 20 + dE = 0.1 + bgRange = [15000,18000] + + if not("Tgrid" in mtd): + + if not(groupedFilename in mtd): + Load(Filename=filename+'.nxs', OutputWorkspace=filename, LoadMonitors=True) + GroupDetectors(InputWorkspace=filename, OutputWorkspace=groupedFilename , MapFile='LET_one2one_123.map', Behaviour='Average') + + wsParent = mtd[groupedFilename]; + + nHist = wsParent.getNumberHistograms(); + print "Parent workspace contains {0:10} histograms".format(nHist) + # Get the energy binning correspondent to the binning produced by rebin function (not to re-implement the same function) + ws1s = ExtractSingleSpectrum(wsParent,0); + ws1s = ConvertUnits(ws1s,'DeltaE','Direct',Ei); + ws1s = Rebin(ws1s,Params=[e_min,dE,e_max]); + e_bins = ws1s.dataX(0); + nBins =e_bins.size; + + x=[e_bins[i] for i in xrange(0,nBins)] + y=[0 for xx in xrange(0,len(x)-1)]*nHist + x = x*nHist + DeleteWorkspace(ws1s); + + eGrid = CreateWorkspace(DataX=x,DataY=y,UnitX='DeltaE',NSpec=nHist,VerticalAxisUnit='SpectraNumber',ParentWorkspace=wsParent) + + Tgrid=ConvertUnits(eGrid,'TOF',Emode='Direct',EFixed=Ei) + + else: + Tgrid = mtd['Tgrid']; + eGrid = mtd['eGrid']; + nHist = Tgrid.getNumberHistograms(); + nBins = Tgrid.dataX(0).size; + + if not('Bg' in mtd): + Bg=Rebin(InputWorkspace=groupedFilename, Params=[bgRange[0],bgRange[1]-bgRange[0],bgRange[1]],PreserveEvents=False) + #Bg=CalculateFlatBackground(InputWorkspace=groupedFilename, StartX=bgRange[0], EndX=bgRange[1], Mode='Mean', OutputMode='Return Background', SkipMonitors=True) + else: + Bg = mtd['Bg'] + + # Assign constant background to the Time grid workspace, minding different time bin width + for nspec in xrange(0,nHist): + bg = Bg.dataY(nspec) + if bg[0]>0: + bgT = Bg.dataX(nspec) + TimeScale = Tgrid.dataX(nspec); + # jacobian for the unit conversion + Jac = (TimeScale[1:nBins]-TimeScale[0:nBins-1])*(bg[0]/(bgT[1]-bgT[0])); + error = np.sqrt(Jac); + eGrid.setY(nspec, Jac) + eGrid.setE(nspec, error) + else: # signal and error for background is 0 anyway. + pass + #print " bg at spectra {0} equal to : {1}".format(nspec,bg[0]) + + + background = eGrid; + resultEt = ConvertUnits(groupedFilename,'DeltaE',Emode='Direct',EFixed=Ei) + result = Rebin(InputWorkspace=resultEt, Params=[e_min,dE,e_max],PreserveEvents=False) + fr = result-background; + # + sourceSum = SumSpectra(result,0,nHist); + bckgrdSum = SumSpectra(background ,0,nHist); + removedBkgSum = SumSpectra(fr ,0,nHist); + +The results of executing this script on workspace contained measured background and the results of the background removal are +presented on the following picture: + +.. image:: /images/BgRemoval.png + +Blue line on this image represents the results, obtained using Rebin with background removal. The results produced using +the script below and shifted by one to show that there is another result plotted on the image, as both results +are identical:: + + from mantid.simpleapi import * + from mantid import config + import numpy as np + import sys + import os + + + maps_dir = '/home/user/InstrumentFiles/let/' + data_dir ='/home/user/results' + ref_data_dir = '/home/user/SystemTests/AnalysisTests/ReferenceResults' + config.setDataSearchDirs('{0};{1};{2}'.format(data_dir,maps_dir,ref_data_dir)) + config['defaultsave.directory'] = data_dir # folder to save resulting spe/nxspe files. Defaults are in + + # the name of a workspace containing background + filename = 'LET0007438' + groupedFilename = filename+'rings'; + # + Ei= 25 + e_min = -20 + e_max = 20 + dE = 0.1 + bgRange = [15000,18000] + + + if not(groupedFilename in mtd): + Load(Filename=filename+'.nxs', OutputWorkspace=filename, LoadMonitors=True) + GroupDetectors(InputWorkspace=filename, OutputWorkspace=groupedFilename , MapFile='LET_one2one_123.map', Behaviour='Average') + + + if not('Bg' in mtd): + Bg=Rebin(InputWorkspace=groupedFilename, Params=[bgRange[0],bgRange[1]-bgRange[0],bgRange[1]],PreserveEvents=False) + else: + Bg = mtd['Bg'] + + if 'resultEtransf' in mtd: + resultEtransf = mtd['resultEtransf'] + else: + resultEtransf = ConvertUnits(groupedFilename,'DeltaE',Emode='Direct',EFixed=Ei) + + noBgWorkspace= Rebin(InputWorkspace=resultEtransf, Params=[e_min,dE,e_max],PreserveEvents=False,FlatBkgWorkspace='Bg',EMode='Direct') + nHist = Bg.getNumberHistograms() + removedBkgSum = SumSpectra(noBgWorkspace ,0,nHist-1); + +.. _rebin-usage: + +Usage +----- + +**Example - simple rebin of a histogram workspace:** + +.. testcode:: ExHistSimple + + # create histogram workspace + dataX = [0,1,2,3,4,5,6,7,8,9] # or use dataX=range(0,10) + dataY = [1,1,1,1,1,1,1,1,1] # or use dataY=[1]*9 + ws = CreateWorkspace(dataX, dataY) + + # rebin from min to max with size bin = 2 + ws = Rebin(ws, 2) + + print "The rebinned X values are: " + str(ws.readX(0)) + print "The rebinned Y values are: " + str(ws.readY(0)) + +Output: + +.. testoutput:: ExHistSimple + + The rebinned X values are: [ 0. 2. 4. 6. 8. 9.] + The rebinned Y values are: [ 2. 2. 2. 2. 1.] + +**Example - logarithmic rebinning:** + +.. testcode:: ExHistLog + + # create histogram workspace + dataX = [1,2,3,4,5,6,7,8,9,10] # or use dataX=range(1,11) + dataY = [1,2,3,4,5,6,7,8,9] # or use dataY=range(1,10) + ws = CreateWorkspace(dataX, dataY) + + # rebin from min to max with logarithmic bins of 0.5 + ws = Rebin(ws, -0.5) + + print "The 2nd and 3rd rebinned X values are: " + str(ws.readX(0)[1:3]) + +Output: + +.. testoutput:: ExHistLog + + The 2nd and 3rd rebinned X values are: [ 1.5 2.25] + +**Example - custom two regions rebinning:** + +.. testcode:: ExHistCustom + + # create histogram workspace + dataX = [0,1,2,3,4,5,6,7,8,9] # or use dataX=range(0,10) + dataY = [0,1,2,3,4,5,6,7,8] # or use dataY=range(0,9) + ws = CreateWorkspace(dataX, dataY) + + # rebin from 0 to 3 in steps of 2 and from 3 to 9 in steps of 3 + ws = Rebin(ws, "1,2,3,3,9") + + print "The rebinned X values are: " + str(ws.readX(0)) + +Output: + +.. testoutput:: ExHistCustom + + The rebinned X values are: [ 1. 3. 6. 9.] + +**Example - use option FullBinsOnly:** + +.. testcode:: ExHistFullBinsOnly + + # create histogram workspace + dataX = [0,1,2,3,4,5,6,7,8,9] # or use dataX=range(0,10) + dataY = [1,1,1,1,1,1,1,1,1] # or use dataY=[1]*9 + ws = CreateWorkspace(dataX, dataY) + + # rebin from min to max with size bin = 2 + ws = Rebin(ws, 2, FullBinsOnly=True) + + print "The rebinned X values are: " + str(ws.readX(0)) + print "The rebinned Y values are: " + str(ws.readY(0)) + +Output: + +.. testoutput:: ExHistFullBinsOnly + + The rebinned X values are: [ 0. 2. 4. 6. 8.] + The rebinned Y values are: [ 2. 2. 2. 2.] + +**Example - use option PreserveEvents:** + +.. testcode:: ExEventRebin + + # create some event workspace + ws = CreateSampleWorkspace(WorkspaceType="Event") + + print "What type is the workspace before 1st rebin: " + str(type(ws)) + # rebin from min to max with size bin = 2 preserving event workspace (default behaviour) + ws = Rebin(ws, 2) + print "What type is the workspace after 1st rebin: " + str(type(ws)) + ws = Rebin(ws, 2, PreserveEvents=False) + print "What type is the workspace after 2nd rebin: " + str(type(ws)) + # note you can also check the type of a workspace using: print isinstance(ws, IEventWorkspace) + +Output: + +.. testoutput:: ExEventRebin + + What type is the workspace before 1st rebin: + What type is the workspace after 1st rebin: + What type is the workspace after 2nd rebin: + +**Example - Background removal during rebinning** + +.. testcode:: ExRebinWithBkgRemoval + + # Create sample workspace with events + Test=CreateSampleWorkspace(WorkspaceType='Event', Function='Flat background') + # Add sample log necessary for unit conversion + AddSampleLog(Test,'Ei',LogText='25.',LogType='Number'); + + # Calculate background + Bg = Rebin(Test,Params='15000,5000,20000',PreserveEvents=False); + + + # Convert event's units + Test_BgDE=ConvertUnits(Test,Target='DeltaE',EMode='Direct'); + + # Calculate histograms for event workspace in energy binning + Sample = Rebin(Test_BgDE,Params='-20,2,20',PreserveEvents=False); + # Calculate histograms for event workspace in energy binning and background removed + Result = Rebin(Test_BgDE,Params='-20,2,20',PreserveEvents=False,FlatBkgWorkspace='Bg',EMode='Direct'); + + # Get access to the results + XS = Sample.dataX(0); + XR = Result .dataX(0); + + YS = Sample.dataY(0); + YR = Result .dataY(0); + + ES = Sample.dataE(0); + ER = Result .dataE(0); + + # print first spectra, Note invalid error calculations + print "| x sampl | x result | S sample | S no bg | Err samp | Err no_bg|" + for i in xrange(0,20): + print "|{0:10}|{1:10}|{2:10.4f}|{3:10.3f}|{4:10.3f}|{5:10.3f}|".format(XS[i],XR[i],YS[i],YR[i],ES[i],ER[i]); + +.. testoutput:: ExRebinWithBkgRemoval + + | x sampl | x result | S sample | S no bg | Err samp | Err no_bg| + | -20.0| -20.0| 1.0000| -0.959| 1.000| 1.216| + | -18.0| -18.0| 2.0000| -0.101| 1.414| 1.432| + | -16.0| -16.0| 3.0000| 0.740| 1.732| 1.622| + | -14.0| -14.0| 1.0000| -1.441| 1.000| 1.312| + | -12.0| -12.0| 5.0000| 2.353| 2.236| 1.955| + | -10.0| -10.0| 2.0000| -0.885| 1.414| 1.563| + | -8.0| -8.0| 5.0000| 1.841| 2.236| 2.020| + | -6.0| -6.0| 2.0000| -1.481| 1.414| 1.655| + | -4.0| -4.0| 4.0000| 0.139| 2.000| 1.983| + | -2.0| -2.0| 3.0000| -1.315| 1.732| 1.912| + | 0.0| 0.0| 6.0000| 1.133| 2.449| 2.331| + | 2.0| 2.0| 7.0000| 1.454| 2.646| 2.505| + | 4.0| 4.0| 5.0000| -1.400| 2.236| 2.388| + | 6.0| 6.0| 7.0000| -0.499| 2.646| 2.692| + | 8.0| 8.0| 9.0000| 0.047| 3.000| 2.996| + | 10.0| 10.0| 11.0000| 0.054| 3.317| 3.313| + | 12.0| 12.0| 16.0000| 2.190| 4.000| 3.861| + | 14.0| 14.0| 16.0000| -2.188| 4.000| 4.135| + | 16.0| 16.0| 26.0000| 0.490| 5.099| 5.075| + | 18.0| 18.0| 39.0000| -0.581| 6.245| 6.268| + + +.. categories::