source files, unit test and documentation to normalize by monitor counts in TOFTOF

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

Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/MonitorEfficiencyCorUser.h

+#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_*/

Code/Mantid/Framework/Algorithms/src/MonitorEfficiencyCorUser.cpp

+#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

Code/Mantid/Framework/Algorithms/test/MonitorEfficiencyCorUserTest.h

+#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_ */

Code/Mantid/docs/source/algorithms/MonitorEfficiencyCorUser-v1.rst

+.. 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::