Skip to content
Snippets Groups Projects
Commit 7cde9908 authored by celinedurniak's avatar celinedurniak
Browse files

Modified algorithms and documentation for MonitorEfficiencyCorUser

parent 57e62192
No related merge requests found
......@@ -4,39 +4,42 @@
#include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h"
namespace Mantid
{
namespace Algorithms
{
namespace Mantid {
namespace Algorithms {
class DLLExport MonitorEfficiencyCorUser : public API::Algorithm
{
class DLLExport MonitorEfficiencyCorUser : public API::Algorithm {
public:
/// (Empty) Constructor
MonitorEfficiencyCorUser();// : Mantid::API::Algorithm() {}
MonitorEfficiencyCorUser(); // : Mantid::API::Algorithm() {}
/// Virtual destructor
virtual ~MonitorEfficiencyCorUser();
/// Algorithm's name
virtual const std::string name() const { return "MonitorEfficiencyCorUser"; }
///Summary of algorithms purpose
virtual const std::string summary() const {return "This algorithm normalises the counts by the monitor counts with additional efficiency correction according to the formula set in the instrument parameters file.";}
/// Summary of algorithms purpose
virtual const std::string summary() const {
return "This algorithm normalises the counts by the monitor counts with "
"additional efficiency correction according to the formula set in "
"the instrument parameters file.";
}
/// Algorithm's version
virtual int version() const { return 1; }
/// Algorithm's category for identification
virtual const std::string category() const { return "CorrectionFunctions\\NormalisationCorrections"; }
virtual const std::string category() const {
return "CorrectionFunctions\\NormalisationCorrections";
}
private:
/// Initialisation code
void init();
///Execution code
/// Execution code
void exec();
double calculateFormulaValue(const std::string&,double);
std::string getValFromInstrumentDef(const std::string&);
void applyMonEfficiency(const size_t numberOfChannels,
const MantidVec& yIn, const MantidVec& eIn, const double coeff,
MantidVec& yOut, MantidVec& eOut);
double calculateFormulaValue(const std::string &, double);
std::string getValFromInstrumentDef(const std::string &);
void applyMonEfficiency(const size_t numberOfChannels, const MantidVec &yIn,
const MantidVec &eIn, const double coeff,
MantidVec &yOut, MantidVec &eOut);
/// The user selected (input) workspace
API::MatrixWorkspace_const_sptr m_inputWS;
......@@ -46,11 +49,9 @@ private:
double m_Ei;
/// stores the total count of neutrons from the monitor
int m_monitorCounts;
};
} // namespace Algorithms
} // namespace Mantid
#endif /*MONITOREFFICIENCYCORUSER_H_*/
......@@ -16,89 +16,87 @@ using namespace Geometry;
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(MonitorEfficiencyCorUser)
//----------------------------------------------------------------------------------------------
/** Constructor
*/
MonitorEfficiencyCorUser::MonitorEfficiencyCorUser() {
}
MonitorEfficiencyCorUser::MonitorEfficiencyCorUser() {}
//----------------------------------------------------------------------------------------------
/** Destructor
*/
MonitorEfficiencyCorUser::~MonitorEfficiencyCorUser() {
}
MonitorEfficiencyCorUser::~MonitorEfficiencyCorUser() {}
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void MonitorEfficiencyCorUser::init()
{
declareProperty(
new WorkspaceProperty<>("InputWorkspace", "", Direction::Input,boost::make_shared<InstrumentValidator>()
), "The workspace to correct for monitor efficiency");
declareProperty(
new WorkspaceProperty<>("OutputWorkspace", "", Direction::Output),
"The name of the workspace in which to store the result.");
void MonitorEfficiencyCorUser::init() {
declareProperty(
new WorkspaceProperty<>("InputWorkspace", "", Direction::Input,
boost::make_shared<InstrumentValidator>()),
"The workspace to correct for monitor efficiency");
declareProperty(
new WorkspaceProperty<>("OutputWorkspace", "", Direction::Output),
"The name of the workspace in which to store the result.");
}
void MonitorEfficiencyCorUser::exec()
{
m_inputWS = this->getProperty("InputWorkspace");
m_outputWS = this->getProperty("OutputWorkspace");
if ( m_inputWS->getInstrument()->getName() != "TOFTOF" )
{ std::string message("The input workspace does not come from TOFTOF");
g_log.error(message);
throw std::invalid_argument(message);
}
// If input and output workspaces are not the same, create a new workspace for the output
if (m_outputWS != this->m_inputWS) {
m_outputWS = API::WorkspaceFactory::Instance().create(m_inputWS);
}
double val;
Strings::convert(m_inputWS->run().getProperty("Ei")->value(),val);
m_Ei = val;
Strings::convert(m_inputWS->run().getProperty("monitor_counts")->value(),m_monitorCounts);
// get Efficiency formula from the IDF - Parameters file
const std::string effFormula = getValFromInstrumentDef("formula_mon_eff");
// Calculate Efficiency for E = Ei
const double eff0 = m_monitorCounts*calculateFormulaValue(effFormula, m_Ei);
const size_t numberOfChannels = this->m_inputWS->blocksize();
// Calculate the number of spectra in this workspace
const int numberOfSpectra = static_cast<int>(this->m_inputWS->getNumberHistograms());
API::Progress prog(this, 0.0, 1.0, numberOfSpectra);
int64_t numberOfSpectra_i = static_cast<int64_t>(numberOfSpectra); // cast to make openmp happy
// Loop over the histograms (detector spectra)
PARALLEL_FOR2(m_outputWS,m_inputWS)
for (int64_t i = 0; i < numberOfSpectra_i; ++i) {
PARALLEL_START_INTERUPT_REGION
//MantidVec& xOut = m_outputWS->dataX(i);
MantidVec& yOut = m_outputWS->dataY(i);
MantidVec& eOut = m_outputWS->dataE(i);
const MantidVec& yIn = m_inputWS->readY(i);
const MantidVec& eIn = m_inputWS->readE(i);
m_outputWS->setX(i, m_inputWS->refX(i));
applyMonEfficiency(numberOfChannels, yIn, eIn, eff0, yOut, eOut);
prog.report("Detector Efficiency correction...");
void MonitorEfficiencyCorUser::exec() {
m_inputWS = this->getProperty("InputWorkspace");
m_outputWS = this->getProperty("OutputWorkspace");
if (m_inputWS->getInstrument()->getName() != "TOFTOF") {
std::string message("The input workspace does not come from TOFTOF");
g_log.error(message);
throw std::invalid_argument(message);
}
// If input and output workspaces are not the same, create a new workspace for
// the output
if (m_outputWS != this->m_inputWS) {
m_outputWS = API::WorkspaceFactory::Instance().create(m_inputWS);
}
double val;
Strings::convert(m_inputWS->run().getProperty("Ei")->value(), val);
m_Ei = val;
Strings::convert(m_inputWS->run().getProperty("monitor_counts")->value(),
m_monitorCounts);
// get Efficiency formula from the IDF - Parameters file
const std::string effFormula = getValFromInstrumentDef("formula_mon_eff");
// Calculate Efficiency for E = Ei
const double eff0 = m_monitorCounts * calculateFormulaValue(effFormula, m_Ei);
const size_t numberOfChannels = this->m_inputWS->blocksize();
// Calculate the number of spectra in this workspace
const int numberOfSpectra =
static_cast<int>(this->m_inputWS->getNumberHistograms());
API::Progress prog(this, 0.0, 1.0, numberOfSpectra);
int64_t numberOfSpectra_i =
static_cast<int64_t>(numberOfSpectra); // cast to make openmp happy
// Loop over the histograms (detector spectra)
PARALLEL_FOR2(m_outputWS, m_inputWS)
for (int64_t i = 0; i < numberOfSpectra_i; ++i) {
PARALLEL_START_INTERUPT_REGION
// MantidVec& xOut = m_outputWS->dataX(i);
MantidVec &yOut = m_outputWS->dataY(i);
MantidVec &eOut = m_outputWS->dataE(i);
const MantidVec &yIn = m_inputWS->readY(i);
const MantidVec &eIn = m_inputWS->readE(i);
m_outputWS->setX(i, m_inputWS->refX(i));
applyMonEfficiency(numberOfChannels, yIn, eIn, eff0, yOut, eOut);
prog.report("Detector Efficiency correction...");
PARALLEL_END_INTERUPT_REGION
} //end for i
PARALLEL_CHECK_INTERUPT_REGION
} // end for i
PARALLEL_CHECK_INTERUPT_REGION
setProperty("OutputWorkspace", m_outputWS);
setProperty("OutputWorkspace", m_outputWS);
}
/**
* Apply the monitor efficiency to a single spectrum
* @param numberOfChannels Number of channels in a spectra (nbins - 1)
......@@ -109,64 +107,61 @@ PARALLEL_CHECK_INTERUPT_REGION
* @param eOut corrected spectrum errors
*/
void MonitorEfficiencyCorUser::applyMonEfficiency(
const size_t numberOfChannels, const MantidVec& yIn,
const MantidVec& eIn, const double effVec, MantidVec& yOut,
MantidVec& eOut) {
for (unsigned int j = 0; j < numberOfChannels; ++j) {
yOut[j] = yIn[j] / effVec;
eOut[j] = eIn[j] / effVec;
}
const size_t numberOfChannels, const MantidVec &yIn, const MantidVec &eIn,
const double effVec, MantidVec &yOut, MantidVec &eOut) {
for (unsigned int j = 0; j < numberOfChannels; ++j) {
yOut[j] = yIn[j] / effVec;
eOut[j] = eIn[j] / effVec;
}
}
/**
* Calculate the value of a formula
* @param formula :: Formula to correct
* @param energy :: value of the energy to use in the formula
* @return calculated corrected value of the monitor count
*/
double MonitorEfficiencyCorUser::calculateFormulaValue(
const std::string &formula, double energy) {
try {
mu::Parser p;
p.DefineVar("e", &energy);
p.SetExpr(formula);
double eff = p.Eval();
g_log.debug() << "Formula: " << formula << " with: " << energy << "evaluated to: "<< eff << std::endl;
return eff;
} catch (mu::Parser::exception_type &e) {
throw Kernel::Exception::InstrumentDefinitionError(
"Error calculating formula from string. Muparser error message is: "
+ e.GetMsg());
}
double
MonitorEfficiencyCorUser::calculateFormulaValue(const std::string &formula,
double energy) {
try {
mu::Parser p;
p.DefineVar("e", &energy);
p.SetExpr(formula);
double eff = p.Eval();
g_log.debug() << "Formula: " << formula << " with: " << energy
<< "evaluated to: " << eff << std::endl;
return eff;
} catch (mu::Parser::exception_type &e) {
throw Kernel::Exception::InstrumentDefinitionError(
"Error calculating formula from string. Muparser error message is: " +
e.GetMsg());
}
}
/**
* Returns the value associated to a parameter name in the IDF
* @param parameterName :: parameter name in the IDF
* @return the value associated to the parameter name
*/
std::string MonitorEfficiencyCorUser::getValFromInstrumentDef(
const std::string& parameterName) {
const ParameterMap& pmap = m_inputWS->constInstrumentParameters();
Instrument_const_sptr instrument = m_inputWS->getInstrument();
Parameter_sptr par = pmap.getRecursive(instrument->getChild(0).get(),
parameterName);
if (par) {
std::string ret = par->asString();
g_log.debug() << "Parsed parameter " << parameterName << ": " << ret
<< "\n";
return ret;
} else {
throw Kernel::Exception::InstrumentDefinitionError(
"There is no <" + parameterName
+ "> in the instrument definition!");
}
const std::string &parameterName) {
const ParameterMap &pmap = m_inputWS->constInstrumentParameters();
Instrument_const_sptr instrument = m_inputWS->getInstrument();
Parameter_sptr par =
pmap.getRecursive(instrument->getChild(0).get(), parameterName);
if (par) {
std::string ret = par->asString();
g_log.debug() << "Parsed parameter " << parameterName << ": " << ret
<< "\n";
return ret;
} else {
throw Kernel::Exception::InstrumentDefinitionError(
"There is no <" + parameterName + "> in the instrument definition!");
}
}
} // namespace Algorithms
......
......@@ -14,126 +14,119 @@ using namespace Mantid;
using namespace Mantid::API;
using namespace Kernel;
class MonitorEfficiencyCorUserTest: public CxxTest::TestSuite {
class MonitorEfficiencyCorUserTest : public CxxTest::TestSuite {
public:
static MonitorEfficiencyCorUserTest *createSuite() {
return new MonitorEfficiencyCorUserTest();
}
static void destroySuite(MonitorEfficiencyCorUserTest *suite) {
delete suite;
}
// constructor
MonitorEfficiencyCorUserTest() :
m_inWSName("input_workspace"), m_outWSName("output_workspace") {
m_Ei = 3.27;
m_monitor_counts = 1000;
createInputWorkSpace();
}
// test init
void test_Init() {
TS_ASSERT_THROWS_NOTHING(alg.initialize())
TS_ASSERT(alg.isInitialized())
}
void testName()
{
TS_ASSERT_EQUALS( alg.name(), "MonitorEfficiencyCorUser" );
}
void testVersion()
{
TS_ASSERT_EQUALS( alg.version(), 1 );
}
void test_exec() {
if (!alg.isInitialized()) alg.initialize();
TS_ASSERT_THROWS_NOTHING(
alg.setPropertyValue("InputWorkspace", m_inWSName));
TS_ASSERT_THROWS_NOTHING(
alg.setPropertyValue("OutputWorkspace", m_outWSName));
TS_ASSERT_THROWS_NOTHING(alg.execute(); );
TS_ASSERT(alg.isExecuted());
// Retrieve the output workspace from data service.
MatrixWorkspace_sptr outWS;
TS_ASSERT_THROWS_NOTHING(
outWS = AnalysisDataService::Instance().retrieveWS<
MatrixWorkspace>(m_outWSName));
TS_ASSERT(outWS);
if (!outWS)
return;
// Retrieve the input workspace from data service.
MatrixWorkspace_sptr inWS;
TS_ASSERT_THROWS_NOTHING(
inWS = AnalysisDataService::Instance().retrieveWS<
MatrixWorkspace>(m_inWSName));
TS_ASSERT(inWS);
if (!inWS)
return;
// Test equality or proportionality between input and output workspaces
const size_t xsize = outWS->blocksize();
double proportionality_coeff = m_monitor_counts*sqrt(m_Ei/25.3);
for(size_t i = 0; i < outWS->getNumberHistograms(); ++i)
{
for(size_t j = 0; j < xsize; ++j)
{ // Same x-values
TS_ASSERT_DELTA(outWS->readX(i)[j], inWS->readX(i)[j], 1e-12);
// Output Y-values proportional to input
TS_ASSERT_DELTA(proportionality_coeff*outWS->readY(i)[j], inWS->readY(i)[j], 1e-12);
// Output Err-values proportional to input
TS_ASSERT_DELTA(proportionality_coeff*outWS->readE(i)[j], inWS->readE(i)[j], 1e-12);
}
}
// Remove workspace from the data service.
AnalysisDataService::Instance().remove(m_outWSName);
AnalysisDataService::Instance().remove(m_inWSName);
}
static MonitorEfficiencyCorUserTest *createSuite() {
return new MonitorEfficiencyCorUserTest();
}
static void destroySuite(MonitorEfficiencyCorUserTest *suite) {
delete suite;
}
// constructor
MonitorEfficiencyCorUserTest()
: m_inWSName("input_workspace"), m_outWSName("output_workspace") {
m_Ei = 3.27;
m_monitor_counts = 1000;
createInputWorkSpace();
}
// test init
void test_Init() {
TS_ASSERT_THROWS_NOTHING(alg.initialize())
TS_ASSERT(alg.isInitialized())
}
void testName() { TS_ASSERT_EQUALS(alg.name(), "MonitorEfficiencyCorUser"); }
void testVersion() { TS_ASSERT_EQUALS(alg.version(), 1); }
void test_exec() {
if (!alg.isInitialized())
alg.initialize();
TS_ASSERT_THROWS_NOTHING(
alg.setPropertyValue("InputWorkspace", m_inWSName));
TS_ASSERT_THROWS_NOTHING(
alg.setPropertyValue("OutputWorkspace", m_outWSName));
TS_ASSERT_THROWS_NOTHING(alg.execute(););
TS_ASSERT(alg.isExecuted());
// Retrieve the output workspace from data service.
MatrixWorkspace_sptr outWS;
TS_ASSERT_THROWS_NOTHING(
outWS = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(
m_outWSName));
TS_ASSERT(outWS);
if (!outWS)
return;
// Retrieve the input workspace from data service.
MatrixWorkspace_sptr inWS;
TS_ASSERT_THROWS_NOTHING(
inWS = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(
m_inWSName));
TS_ASSERT(inWS);
if (!inWS)
return;
// Test equality or proportionality between input and output workspaces
const size_t xsize = outWS->blocksize();
double proportionality_coeff = m_monitor_counts * sqrt(m_Ei / 25.3);
for (size_t i = 0; i < outWS->getNumberHistograms(); ++i) {
for (size_t j = 0; j < xsize; ++j) { // Same x-values
TS_ASSERT_DELTA(outWS->readX(i)[j], inWS->readX(i)[j], 1e-12);
// Output Y-values proportional to input
TS_ASSERT_DELTA(proportionality_coeff * outWS->readY(i)[j],
inWS->readY(i)[j], 1e-12);
// Output Err-values proportional to input
TS_ASSERT_DELTA(proportionality_coeff * outWS->readE(i)[j],
inWS->readE(i)[j], 1e-12);
}
}
// Remove workspace from the data service.
AnalysisDataService::Instance().remove(m_outWSName);
AnalysisDataService::Instance().remove(m_inWSName);
}
private:
double m_Ei;
int m_monitor_counts;
const std::string m_inWSName, m_outWSName;
MonitorEfficiencyCorUser alg;
double m_Ei;
int m_monitor_counts;
const std::string m_inWSName, m_outWSName;
MonitorEfficiencyCorUser alg;
void createInputWorkSpace() {
int numHist = 1;
int numBins = 20;
void createInputWorkSpace() {
int numHist = 1;
int numBins = 20;
DataObjects::Workspace2D_sptr dataws = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(
numHist, numBins, false, false, true, "TOFTOF");
DataObjects::Workspace2D_sptr dataws =
WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(
numHist, numBins, false, false, true, "TOFTOF");
static const double arr[] = {-10. , -9.25, -8.5 , -7.75, -7. , -6.25, -5.5 , -4.75,
-4. , -3.25, -2.5 , -1.75, -1. , -0.25, 0.5 , 1.25,
2. , 2.75, 3.5 , 4.25, 5. };
static const double arr[] = {
-10., -9.25, -8.5, -7.75, -7., -6.25, -5.5, -4.75, -4., -3.25, -2.5,
-1.75, -1., -0.25, 0.5, 1.25, 2., 2.75, 3.5, 4.25, 5.};
std::vector<double> vec (arr, arr + sizeof(arr) / sizeof(arr[0]) );
std::vector<double> vec(arr, arr + sizeof(arr) / sizeof(arr[0]));
for (size_t wi = 0; wi < dataws->getNumberHistograms(); wi++) {
dataws->setX(wi, vec);
}
for (size_t wi = 0; wi < dataws->getNumberHistograms(); wi++) {
dataws->setX(wi, vec);
}
dataws->getAxis(0)->setUnit("TOF");
dataws->mutableRun().addProperty("Ei",
boost::lexical_cast<std::string>(m_Ei));
dataws->mutableRun().addProperty(
"monitor_counts", boost::lexical_cast<std::string>(m_monitor_counts));
dataws->getAxis(0)->setUnit("TOF");
dataws->mutableRun().addProperty("Ei",
boost::lexical_cast<std::string>(m_Ei));
dataws->mutableRun().addProperty("monitor_counts",
boost::lexical_cast<std::string>(m_monitor_counts));
dataws->instrumentParameters().addString(
dataws->getInstrument()->getChild(0).get(), "formula_mon_eff", "sqrt(e/25.3)"); //TOFTOF
API::AnalysisDataService::Instance().addOrReplace(m_inWSName, dataws);
}
dataws->instrumentParameters().addString(
dataws->getInstrument()->getChild(0).get(), "formula_mon_eff",
"sqrt(e/25.3)"); // TOFTOF
API::AnalysisDataService::Instance().addOrReplace(m_inWSName, dataws);
}
};
#endif /* MANTID_ALGORITHMS_MONITOREFFICIENCYCORUSERTEST_H_ */
......@@ -9,20 +9,24 @@
Description
-----------
Normalizes the counts by monitor counts with additional efficiency correction.
This algorithm normalises the neutron counts by monitor counts with an additional efficiency correction.
Corrects the monitor efficiency according to the formula contained in TOFTOF's Parameters file.
To date this algorithm only supports the TOFTOF instrument.
The values of the monitor counts and incident energy are stored in the SampleLogs of the input workspace.
The monitor counts is the total count and it is stored in the SampleLogs of the input workspace.
To date this algorithm only supports the TOFTOF instrument.
This count is corrected taking into account the monitor efficiency. The formula used for the correction is stored in the Parameters file and requires the incident energy (Ei), which is stored in the SampleLogs of the input workspace.
The corrected value of the monitor counts is used to normalise the input workspace.
Restrictions
###################################
A formula named "formula\_eff" must be defined in the instrument
parameters file.
parameters file. It is defined as "monitor counts * sqrt(Ei/25.3)"
The incident energy Ei and the monitor counts are read in the SampleLogs of the input workspace.
Usage
-----
......@@ -32,27 +36,32 @@ Usage
.. testcode:: ExMonitorEfficiencyCorUser
import numpy as np
#Load data
wsSample = LoadMLZ(Filename='TOFTOFTestdata.nxs')
wsNorm = MonitorEfficiencyCorUser(wsSample)
# Check calculation of normalization coefficient
wsNorm = MonitorEfficiencyCorUser(wsSample )
# Input and output workspaces have the same structure
print 'Number of histograms of the input and output workspaces:'
print wsSample.getNumberHistograms(), wsNorm.getNumberHistograms()
print 'Number of time channels of the input an output workspaces:'
print wsSample.blocksize(), wsNorm.blocksize()
# Check calculation of normalisation coefficient between input and output workspaces
wsCheck = Divide(wsSample,wsNorm)
print 'Coefficient of proportionality between Input and Output of MonitorEfficiencyCorUser algorithm:'
print wsCheck.readY(102)[1]
print "Coefficient of proportionality between Input and Output of MonitorEfficiencyCorUser algorithm: %5.3f" % wsCheck.readY(102)[1]
# Read the values of the incident energy and of the monitor counts from the SampleLogs of wsSample
monitor_counts = float(mtd['wsSample'].getRun().getLogData('monitor_counts').value)
Ei = float(mtd['wsSample'].getRun().getLogData('Ei').value)
print 'Coefficient from theoretical formula = monitor_counts * sqrt(Ei/25.3) ', monitor_counts*np.sqrt(Ei/25.3)
print "Coefficient from theoretical formula = monitor_counts * sqrt(Ei/25.3): %5.3f" % (monitor_counts*np.sqrt(Ei/25.3))
Output:
.. testoutput:: ExMonitorEfficiencyCorUser
Coefficient of proportionality between Input and Output of MonitorEfficiencyCorUser algorithm:
41038.4323645
Coefficient from theoretical formula = monitor_counts * sqrt(Ei/25.3) 41038.4323645
Number of histograms of the input and output workspaces:
1006 1006
Number of time channels of the input an output workspaces:
1024 1024
Coefficient of proportionality between Input and Output of MonitorEfficiencyCorUser algorithm: 41038.432
Coefficient from theoretical formula = monitor_counts * sqrt(Ei/25.3): 41038.432
.. categories::
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment