From ba10bc8873c0da1c3fc2037e7ff7e75971917b35 Mon Sep 17 00:00:00 2001
From: Celine Durniak <c.durniak@fz-juelich.de>
Date: Mon, 20 Jul 2015 14:48:39 +0200
Subject: [PATCH] source files, unit test and documentation to normalize by
monitor counts in TOFTOF
---
.../Framework/Algorithms/CMakeLists.txt | 3 +
.../MonitorEfficiencyCorUser.h | 56 ++++++
.../src/MonitorEfficiencyCorUser.cpp | 173 ++++++++++++++++++
.../test/MonitorEfficiencyCorUserTest.h | 139 ++++++++++++++
.../MonitorEfficiencyCorUser-v1.rst | 59 ++++++
5 files changed, 430 insertions(+)
create mode 100644 Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/MonitorEfficiencyCorUser.h
create mode 100644 Code/Mantid/Framework/Algorithms/src/MonitorEfficiencyCorUser.cpp
create mode 100644 Code/Mantid/Framework/Algorithms/test/MonitorEfficiencyCorUserTest.h
create mode 100644 Code/Mantid/docs/source/algorithms/MonitorEfficiencyCorUser-v1.rst
diff --git a/Code/Mantid/Framework/Algorithms/CMakeLists.txt b/Code/Mantid/Framework/Algorithms/CMakeLists.txt
index 88c34d71dd2..d06562cce9f 100644
--- a/Code/Mantid/Framework/Algorithms/CMakeLists.txt
+++ b/Code/Mantid/Framework/Algorithms/CMakeLists.txt
@@ -154,6 +154,7 @@ set ( SRC_FILES
src/Minus.cpp
src/ModeratorTzero.cpp
src/ModeratorTzeroLinear.cpp
+ src/MonitorEfficiencyCorUser.cpp
src/MonteCarloAbsorption.cpp
src/MultipleScatteringCylinderAbsorption.cpp
src/Multiply.cpp
@@ -421,6 +422,7 @@ set ( INC_FILES
inc/MantidAlgorithms/Minus.h
inc/MantidAlgorithms/ModeratorTzero.h
inc/MantidAlgorithms/ModeratorTzeroLinear.h
+ inc/MantidAlgorithms/MonitorEfficiencyCorUser.h
inc/MantidAlgorithms/MonteCarloAbsorption.h
inc/MantidAlgorithms/MultipleScatteringCylinderAbsorption.h
inc/MantidAlgorithms/Multiply.h
@@ -688,6 +690,7 @@ set ( TEST_FILES
MinusTest.h
ModeratorTzeroLinearTest.h
ModeratorTzeroTest.h
+ MonitorEfficiencyCorUserTest.h
MonteCarloAbsorptionTest.h
MultipleScatteringCylinderAbsorptionTest.h
MultiplyRangeTest.h
diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/MonitorEfficiencyCorUser.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/MonitorEfficiencyCorUser.h
new file mode 100644
index 00000000000..3fe3c1b2245
--- /dev/null
+++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/MonitorEfficiencyCorUser.h
@@ -0,0 +1,56 @@
+#ifndef MANTID_ALGORITHMS_MONITOREFFICIENCYCORUSER_H_
+#define MANTID_ALGORITHMS_MONITOREFFICIENCYCORUSER_H_
+
+#include "MantidKernel/System.h"
+#include "MantidAPI/Algorithm.h"
+
+namespace Mantid
+{
+ namespace Algorithms
+ {
+
+class DLLExport MonitorEfficiencyCorUser : public API::Algorithm
+{
+public:
+ /// (Empty) Constructor
+ 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.";}
+
+ /// Algorithm's version
+ virtual int version() const { return 1; }
+ /// Algorithm's category for identification
+ virtual const std::string category() const { return "MLZ;CorrectionFunctions\\NormalisationCorrections"; }
+
+private:
+ /// Initialisation code
+ void init();
+ ///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);
+
+ /// The user selected (input) workspace
+ API::MatrixWorkspace_const_sptr m_inputWS;
+ /// The output workspace, maybe the same as the input one
+ API::MatrixWorkspace_sptr m_outputWS;
+ /// stores the incident energy of the neutrons
+ double m_Ei;
+ /// stores the total count of neutrons from the monitor
+ int m_monitorCounts;
+
+};
+
+} // namespace Algorithms
+} // namespace Mantid
+
+
+#endif /*MONITOREFFICIENCYCORUSER_H_*/
diff --git a/Code/Mantid/Framework/Algorithms/src/MonitorEfficiencyCorUser.cpp b/Code/Mantid/Framework/Algorithms/src/MonitorEfficiencyCorUser.cpp
new file mode 100644
index 00000000000..8aa2d99836f
--- /dev/null
+++ b/Code/Mantid/Framework/Algorithms/src/MonitorEfficiencyCorUser.cpp
@@ -0,0 +1,173 @@
+#include "MantidAlgorithms/MonitorEfficiencyCorUser.h"
+#include "MantidAPI/WorkspaceValidators.h"
+#include "MantidKernel/BoundedValidator.h"
+#include "MantidKernel/CompositeValidator.h"
+#include "MantidGeometry/muParser_Silent.h"
+#include "MantidKernel/MultiThreaded.h"
+#include <ctime>
+
+namespace Mantid {
+namespace Algorithms {
+
+using namespace Kernel;
+using namespace API;
+using namespace Geometry;
+
+// Register the algorithm into the AlgorithmFactory
+DECLARE_ALGORITHM(MonitorEfficiencyCorUser)
+
+
+//----------------------------------------------------------------------------------------------
+/** Constructor
+ */
+MonitorEfficiencyCorUser::MonitorEfficiencyCorUser() {
+}
+
+//----------------------------------------------------------------------------------------------
+/** Destructor
+ */
+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::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
+
+ setProperty("OutputWorkspace", m_outputWS);
+}
+
+
+/**
+ * Apply the monitor efficiency to a single spectrum
+ * @param numberOfChannels Number of channels in a spectra (nbins - 1)
+ * @param yIn spectrum counts
+ * @param eIn spectrum errors
+ * @param effVec efficiency values (to be divided by the counts)
+ * @param yOut corrected spectrum counts
+ * @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;
+ }
+
+}
+
+
+/**
+ * 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());
+ }
+}
+
+
+/**
+ * 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!");
+ }
+}
+
+} // namespace Algorithms
+} // namespace Mantid
diff --git a/Code/Mantid/Framework/Algorithms/test/MonitorEfficiencyCorUserTest.h b/Code/Mantid/Framework/Algorithms/test/MonitorEfficiencyCorUserTest.h
new file mode 100644
index 00000000000..5eb7527dc66
--- /dev/null
+++ b/Code/Mantid/Framework/Algorithms/test/MonitorEfficiencyCorUserTest.h
@@ -0,0 +1,139 @@
+#ifndef MANTID_ALGORITHMS_MONITORTOREFFICIENCYCORUSERTEST_H_
+#define MANTID_ALGORITHMS_MONITOREFFICIENCYCORUSERTEST_H_
+
+#include <cxxtest/TestSuite.h>
+
+#include "MantidAlgorithms/MonitorEfficiencyCorUser.h"
+#include "MantidDataObjects/Workspace2D.h"
+#include "MantidTestHelpers/WorkspaceCreationHelper.h"
+#include "MantidTestHelpers/ComponentCreationHelper.h"
+
+using Mantid::Algorithms::MonitorEfficiencyCorUser;
+
+using namespace Mantid;
+using namespace Mantid::API;
+using namespace Kernel;
+
+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);
+ }
+
+private:
+ double m_Ei;
+ int m_monitor_counts;
+ const std::string m_inWSName, m_outWSName;
+ MonitorEfficiencyCorUser alg;
+
+ void createInputWorkSpace() {
+ int numHist = 1;
+ int numBins = 20;
+
+ 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. };
+
+ std::vector<double> vec (arr, arr + sizeof(arr) / sizeof(arr[0]) );
+
+ 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->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_ */
diff --git a/Code/Mantid/docs/source/algorithms/MonitorEfficiencyCorUser-v1.rst b/Code/Mantid/docs/source/algorithms/MonitorEfficiencyCorUser-v1.rst
new file mode 100644
index 00000000000..2205cca40d1
--- /dev/null
+++ b/Code/Mantid/docs/source/algorithms/MonitorEfficiencyCorUser-v1.rst
@@ -0,0 +1,59 @@
+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+Normalizes the counts by monitor counts with additional efficiency correction.
+
+Corrects the monitor efficiency according to the formula contained in TOFTOF's Parameters file.
+
+The values of the monitor counts and incident energy are stored in the SampleLogs of the input workspace.
+
+To date this algorithm only supports the TOFTOF instrument.
+
+
+Restrictions
+###################################
+
+A formula named "formula\_eff" must be defined in the instrument
+parameters file.
+
+Usage
+-----
+
+**Example**
+
+.. testcode:: ExMonitorEfficiencyCorUser
+
+ import numpy as np
+ #Load data
+
+ wsSample = LoadMLZ(Filename='TOFTOFTestdata.nxs')
+
+ wsNorm = MonitorEfficiencyCorUser(wsSample)
+ # Check calculation of normalization coefficient
+ wsCheck = Divide(wsSample,wsNorm)
+ print 'Coefficient of proportionality between Input and Output of MonitorEfficiencyCorUser algorithm:'
+ print wsCheck.readY(102)[1]
+ 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)
+
+
+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
+
+.. categories::
+
+.. sourcelink::
--
GitLab