Skip to content

Commit

Permalink
refs #10384 Mainly finished helper and first unit test working
Browse files Browse the repository at this point in the history
  • Loading branch information
abuts committed Oct 28, 2014
1 parent 15e1f5b commit 466f9ad
Show file tree
Hide file tree
Showing 4 changed files with 102 additions and 16 deletions.
Expand Up @@ -7,12 +7,13 @@

#include "MantidAPI/DllConfig.h"
#include "MantidAPI/MatrixWorkspace.h"
#include <forward_list>

namespace Mantid
{
namespace Algorithms
{
/** Class helping to remove constant (and possibly non-constant after simple modification) background calculated in TOF units
/** 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
Expand Down Expand Up @@ -44,7 +45,10 @@ namespace Algorithms
BackgroundHelper();
void initialize(const API::MatrixWorkspace_const_sptr &bkgWS,const API::MatrixWorkspace_sptr &sourceWS,int emode);

void removeBackground(int hist,const MantidVec &XValues,MantidVec &y_data,MantidVec &e_data);
void removeBackground(int hist,const MantidVec &XValues,MantidVec &y_data,MantidVec &e_data)const;

//returns the list of the failing detectors
std::forward_list<int> & getFailingSpectrsList()const{return FailingSpectraList;}
private:
// pointer to the units conversion class for the working workspace;
Kernel::Unit_sptr m_WSUnit;
Expand All @@ -56,16 +60,24 @@ namespace Algorithms

// if the background workspace is single value workspace
bool m_singleValueBackground;
// the intensity of the background for a given spectra of a background workspace
// the intensity of the background for first spectra of a background workspace
double m_Jack05;
// 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;


// list of the spectra numbers for which detectors retrieval has been unsuccessful
mutable std::forward_list<int> FailingSpectraList;
};

}
Expand Down
72 changes: 62 additions & 10 deletions Code/Mantid/Framework/Algorithms/src/BackgroundHelper.cpp
Expand Up @@ -5,12 +5,16 @@ namespace Mantid
{
namespace Algorithms
{

/// Constructor
BackgroundHelper::BackgroundHelper():
m_WSUnit(),m_bgWs(),m_wkWS(),
m_singleValueBackground(false),
m_Jack05(0),m_ErrSq(0),
m_Emode(0),
m_Efix(0),
m_Jack05(0)
m_L1(0),m_Efix(0),
m_Sample(),
FailingSpectraList()
{};

/** Initialization method:
Expand All @@ -22,6 +26,7 @@ namespace Mantid
{
m_bgWs = bkgWS;
m_wkWS = sourceWS;
FailingSpectraList.clear();

std::string bgUnits = bkgWS->getAxis(0)->unit()->unitID();
if(bgUnits!="TOF")
Expand All @@ -44,30 +49,36 @@ namespace Mantid
m_singleValueBackground = false;
if(bkgWS->getNumberHistograms()==0)
m_singleValueBackground = true;
const MantidVec& dataY = bkgWS->dataY(0);
const MantidVec& dataX = bkgWS->dataX(0);
const MantidVec& dataY = bkgWS->dataY(0);
const MantidVec& dataE = bkgWS->dataE(0);
m_Jack05 = dataY[0]/(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
*/
void BackgroundHelper::removeBackground(int nHist,const MantidVec &XValues,MantidVec &y_data,MantidVec &e_data)
void BackgroundHelper::removeBackground(int nHist,const MantidVec &XValues,MantidVec &y_data,MantidVec &e_data)const
{
double Jack05;
double Jack05,ErrSq;
if(m_singleValueBackground)
{
Jack05= m_Jack05;
ErrSq = m_ErrSq;
}
else
{
const MantidVec& dataY = m_bgWs->dataY(nHist);
const MantidVec& dataX = m_bgWs->dataX(nHist);
const MantidVec& dataY = m_bgWs->dataY(nHist);
const MantidVec& dataE = m_bgWs->dataX(nHist);
Jack05 = dataY[0]/(dataX[1]-dataX[0]);
ErrSq = dataE[0]*dataE[0]; // Needs further clarification
}

try
Expand All @@ -84,18 +95,59 @@ namespace Mantid
double tof2=m_WSUnit->convertSingleToTOF(XValues[i+1], m_L1, L2,twoTheta, m_Emode, m_Efix, delta);
double normBkgrnd = (tof1-tof2)*Jack05;
tof1=tof2;
double normErr = std::sqrt(normBkgrnd);
y_data[i-1] -=normBkgrnd;
e_data[i-1] +=normErr;
y_data[i] -=normBkgrnd;
e_data[i] =std::sqrt((ErrSq+e_data[i]*e_data[i])/2.); // needs further clarification -- Gaussian error summation
}

}
catch(...)
{
// no background removal for this detector; How should it be reported?
// no background removal for this spectra as it does not have a detector; 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<double>::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<double>("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<double>("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
24 changes: 23 additions & 1 deletion Code/Mantid/Framework/Algorithms/test/BackgroundHelperTest.h
Expand Up @@ -73,13 +73,34 @@ class BackgroundHelperTest : public CxxTest::TestSuite

unitsConv.execute();

SourceWS = API::AnalysisDataService::Instance().retrieveWS<API::MatrixWorkspace>("sampleWSdE");


}

~BackgroundHelperTest()
{
BgWS.reset();
SourceWS.reset();
}
void testWrongInit()
{
Algorithms::BackgroundHelper bgRemoval;
// create workspace with units of energy transfer
auto bkgWS = WorkspaceCreationHelper::createProcessedInelasticWS(std::vector<double>(1,1.), std::vector<double>(1,20.), std::vector<double>(1,10.));
TSM_ASSERT_THROWS("Should throw if background workspace is not in TOF units",bgRemoval.initialize(bkgWS,SourceWS,0),std::invalid_argument);


bkgWS = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(2, 15);
TSM_ASSERT_THROWS("Should throw if background is not 1 or equal to source",bgRemoval.initialize(bkgWS,SourceWS,0),std::invalid_argument);

auto sourceWS = WorkspaceCreationHelper::Create2DWorkspace(5,10);
TSM_ASSERT_THROWS("Should throw if source workspace does not have units",bgRemoval.initialize(BgWS,sourceWS,0),std::invalid_argument);

sourceWS ->getAxis(0)->setUnit("TOF");
TSM_ASSERT_THROWS("Should throw if source workspace does not have proper instrument",bgRemoval.initialize(BgWS,sourceWS,0),std::invalid_argument);
}


void testBackgroundHelper()
{
Expand All @@ -93,7 +114,8 @@ class BackgroundHelperTest : public CxxTest::TestSuite
}

private:
API::MatrixWorkspace_sptr BgWS;
API::MatrixWorkspace_sptr BgWS;
API::MatrixWorkspace_sptr SourceWS;

};

Expand Down
4 changes: 2 additions & 2 deletions Code/Mantid/Framework/Kernel/inc/MantidKernel/Unit.h
Expand Up @@ -284,9 +284,9 @@ class MANTID_KERNEL_DLL TOF : public Unit
virtual double singleToTOF(const double x) const;
virtual double singleFromTOF(const double tof) const;
virtual Unit * clone() const;
///@return -DBL_MAX as ToF convetanble to TOF for in any time range
///@return -DBL_MAX as ToF convertible to TOF for in any time range
virtual double conversionTOFMin()const;
///@return DBL_MAX as ToF convetanble to TOF for in any time range
///@return DBL_MAX as ToF convertible to TOF for in any time range
virtual double conversionTOFMax()const;
};

Expand Down

0 comments on commit 466f9ad

Please sign in to comment.