From 0cf37b67875ea1c8413fd9a56732688212785ee5 Mon Sep 17 00:00:00 2001 From: Alex Buts Date: Thu, 30 Oct 2014 19:50:50 +0000 Subject: [PATCH] refs #10384 Separated BackgroundHelper from Rebin and moved it to RemoveBackground. --- .../inc/MantidAlgorithms/BackgroundHelper.h | 90 ---- .../Algorithms/inc/MantidAlgorithms/Rebin.h | 2 - .../inc/MantidAlgorithms/RemoveBackground.h | 94 ++--- .../Algorithms/src/BackgroundHelper.cpp | 193 --------- .../Algorithms/src/RemoveBackground.cpp | 399 +++++------------- 5 files changed, 143 insertions(+), 635 deletions(-) delete mode 100644 Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/BackgroundHelper.h delete mode 100644 Code/Mantid/Framework/Algorithms/src/BackgroundHelper.cpp diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/BackgroundHelper.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/BackgroundHelper.h deleted file mode 100644 index cb1b9ea22b93..000000000000 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/BackgroundHelper.h +++ /dev/null @@ -1,90 +0,0 @@ -#ifndef MANTID_ALGORITHM_BACKGROUNDHELPER_H_ -#define MANTID_ALGORITHM_BACKGROUNDHELPER_H_ - - -#include "MantidKernel/Unit.h" -#include "MantidKernel/cow_ptr.h" - -#include "MantidAPI/DllConfig.h" -#include "MantidAPI/MatrixWorkspace.h" -#include - -namespace Mantid -{ -namespace Algorithms -{ -/** Class provides 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 © 2014 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 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; - - //returns the list of the failing detectors - std::list & getFailingSpectrsList()const{return FailingSpectraList;} - 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; - - // list of the spectra numbers for which detectors retrieval has been unsuccessful - mutable std::list FailingSpectraList; - - // 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(); - }; - -} -} -#endif \ No newline at end of file diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Rebin.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Rebin.h index da28874ca486..3f88aea104ba 100644 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Rebin.h +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/Rebin.h @@ -5,8 +5,6 @@ // Includes //---------------------------------------------------------------------- #include "MantidAPI/Algorithm.h" -#include "MantidAlgorithms/BackgroundHelper.h" - namespace Mantid { namespace Algorithms diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/RemoveBackground.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/RemoveBackground.h index 03139d781566..a484bde4560c 100644 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/RemoveBackground.h +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/RemoveBackground.h @@ -5,8 +5,6 @@ // Includes //---------------------------------------------------------------------- #include "MantidAPI/Algorithm.h" -#include "MantidAlgorithms/BackgroundHelper.h" - namespace Mantid { namespace Algorithms @@ -37,7 +35,52 @@ namespace Mantid File change history is stored at: Code Documentation is available at: */ - + /**Class actually performing background removal from a workspace spectra */ + class DLLExport BackgroundHelper + { + public: + BackgroundHelper(); + ~BackgroundHelper(); + void initialize(const API::MatrixWorkspace_const_sptr &bkgWS,const API::MatrixWorkspace_sptr &sourceWS, + int emode,int nTreads=1,bool inPlace=true); + void removeBackground(int hist,const MantidVec &XValues,MantidVec &y_data,MantidVec &e_data, + int tread_num=0,Kernel::Logger *pLog=NULL)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; + + // perform background removal in-place + bool m_inPlace; + + // 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(); + }; + class DLLExport RemoveBackground : public API::Algorithm { public: @@ -62,54 +105,9 @@ namespace Mantid 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(); - }; }; diff --git a/Code/Mantid/Framework/Algorithms/src/BackgroundHelper.cpp b/Code/Mantid/Framework/Algorithms/src/BackgroundHelper.cpp deleted file mode 100644 index 4519ede8ff88..000000000000 --- a/Code/Mantid/Framework/Algorithms/src/BackgroundHelper.cpp +++ /dev/null @@ -1,193 +0,0 @@ -#include "MantidAlgorithms/BackgroundHelper.h" - - -namespace Mantid -{ - namespace Algorithms - { - - /// 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; - } - - - } // end API -} // end Mantid diff --git a/Code/Mantid/Framework/Algorithms/src/RemoveBackground.cpp b/Code/Mantid/Framework/Algorithms/src/RemoveBackground.cpp index f5b0e909fb68..90cd63d4f760 100644 --- a/Code/Mantid/Framework/Algorithms/src/RemoveBackground.cpp +++ b/Code/Mantid/Framework/Algorithms/src/RemoveBackground.cpp @@ -24,7 +24,7 @@ namespace Mantid using namespace Kernel; using namespace API; - + //--------------------------------------------------------------------------------------------- // Public methods //--------------------------------------------------------------------------------------------- @@ -67,323 +67,103 @@ namespace Mantid */ 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); - + API::MatrixWorkspace_const_sptr bkgWksp = getProperty("FlatBkgWorkspace"); - // 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) + // source workspace has to have full instrument defined to perform background removal using this procedure. + auto pInstrument = inputWS->getInstrument(); + if(pInstrument) { - int nThreads = omp_get_max_threads(); - m_BackgroundHelper.initialize(bkgWS,inputWS,eMode,nThreads); + if(!pInstrument->getSample()) + throw std::invalid_argument(" Workspace: "+inputWS->getName()+" does not have properly defined instrument. Can not remove background"); } - //------------------------------------------------------- - - 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) + if(!(bkgWksp->getNumberHistograms() == 1 || inputWS->getNumberHistograms()==bkgWksp->getNumberHistograms())) { - //------- 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 + throw std::invalid_argument(" Background Workspace: "+bkgWksp->getName()+" should have the same number of spectra as source workspace or be a single histogram workspace"); + } - else + // + int eMode; // in convert units emode is still integer + const std::string emodeStr = getProperty("EMode"); + eMode = static_cast(Kernel::DeltaEMode().fromString(emodeStr)); - { //------- Workspace2D or other MatrixWorkspace --------------------------- + // Removing background in-place + bool inPlace = (inputWS == outputWS); + // + int nThreads = omp_get_max_threads(); + m_BackgroundHelper.initialize(bkgWksp,inputWS,eMode,nThreads,inPlace); - 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); + // workspace independent determination of length + const int histnumber = static_cast(inputWS->getNumberHistograms()); - // Copy over the 'vertical' axis - if (inputWS->axes() > 1) outputWS->replaceAxis( 1, inputWS->getAxis(1)->clone(outputWS.get()) ); + if(!inPlace) + { + // make the copy of the output Workspace from the input. Also copies X-vectors and axis + outputWS = API::WorkspaceFactory::Instance().create(inputWS); + } - 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"); - } + 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 = outputWS->readX(hist); + MantidVec& YValues = outputWS->dataY(hist); + MantidVec& YErrors = outputWS->dataE(hist); + + // output data arrays are implicitly filled by function + int id = omp_get_thread_num(); + m_BackgroundHelper.removeBackground(hist,XValues,YValues,YErrors,id,&g_log); + + prog.report(name()); + PARALLEL_END_INTERUPT_REGION + } + PARALLEL_CHECK_INTERUPT_REGION - 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; - } -//------------------------------------------------------------------------------------------------------------------------------- + //------------------------------------------------------------------------------------------------------------------------------- + //---------------- BACKGROUND HELPER --------------------------------------------------------- + //------------------------------------------------------------------------------------------------------------------------------- /// Constructor BackgroundHelper::BackgroundHelper(): m_WSUnit(), m_bgWs(),m_wkWS(), + m_inPlace(true), m_singleValueBackground(false), m_NBg(0),m_dtBg(1), //m_ErrSq(0), m_Emode(0), m_L1(0),m_Efix(0), - m_Sample(), - FailingSpectraList() + m_Sample() {}; + /// Desctructoe BackgroundHelper::~BackgroundHelper() { this->deleteUnitsConverters(); } - /** The method deletes all units converter allocated*/ - void BackgroundHelper::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"); @@ -455,14 +235,14 @@ namespace Mantid * @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 + void BackgroundHelper::removeBackground(int nHist,const MantidVec &XValues,MantidVec &y_data,MantidVec &e_data,int threadNum,Kernel::Logger *pLog)const { double dtBg,IBg; if(m_singleValueBackground) { dtBg = m_dtBg; - // ErrSq = m_ErrSq; + // ErrSq = m_ErrSq; IBg = m_NBg; } else @@ -474,7 +254,7 @@ namespace Mantid IBg = dataY[0]; //ErrSq= dataE[0]*dataE[0]; // Needs further clarification } - + try { auto detector = m_wkWS->getDetector(nHist); @@ -482,9 +262,15 @@ namespace Mantid double twoTheta = m_wkWS->detectorTwoTheta(detector); double L2 = detector->getDistance(*m_Sample); double delta(std::numeric_limits::quiet_NaN()); + // get access to source workspace in case if target is different from source + const MantidVec& YValues = m_wkWS->readY(nHist); + const MantidVec& YErrors = m_wkWS->readE(nHist); + // 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;i 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? + if(m_inPlace) + { + 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? + } + else + { + y_data[i] =YValues[i]-normBkgrnd; + e_data[i] =std::sqrt((normBkgrnd+YErrors[i]*YErrors[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); + // no background removal for this spectra as it does not have a detector or other reason + if(pLog) + pLog->debug()<<" Can not remove background for the spectra with number (id)"<::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; + Efi = inputWS->run().getPropertyValueAsType("Ei"); + EiFound = true; } catch(Kernel::Exception::NotFoundError &) {} @@ -534,8 +330,8 @@ namespace Mantid { try { - Efi =inputWS->run().getPropertyValueAsType("eFixed"); - //eFixedFound = true; + Efi =inputWS->run().getPropertyValueAsType("eFixed"); + //eFixedFound = true; } catch(Kernel::Exception::NotFoundError &) {} @@ -548,6 +344,5 @@ namespace Mantid } - } // namespace Algorithm } // namespace Mantid