Detailed balance correction algorithm. Updated physical constants to calculate...

Detailed balance correction algorithm. Updated physical constants to calculate energy to wavevector constant. Updated physical constants with energy to temperature conversion constant. Refs #3659

Code/Mantid/Framework/Algorithms/CMakeLists.txt

@@ -4,6 +4,7 @@ set ( SRC_FILES
 	src/AddSampleLog.cpp
 	src/AlignDetectors.cpp
 	src/AnyShapeAbsorption.cpp
+	src/ApplyDetailedBalance.cpp
 	src/ApplyTransmissionCorrection.cpp
 	src/BinaryOperation.cpp
 	src/BlendSq.cpp
@@ -167,6 +168,7 @@ set ( INC_FILES
 	inc/MantidAlgorithms/AddSampleLog.h
 	inc/MantidAlgorithms/AlignDetectors.h
 	inc/MantidAlgorithms/AnyShapeAbsorption.h
+	inc/MantidAlgorithms/ApplyDetailedBalance.h
 	inc/MantidAlgorithms/ApplyTransmissionCorrection.h
 	inc/MantidAlgorithms/BinaryOperation.h
 	inc/MantidAlgorithms/BlendSq.h
@@ -329,6 +331,7 @@ set ( TEST_FILES
 	test/AddSampleLogTest.h
 	test/AlignDetectorsTest.h
 	test/AnyShapeAbsorptionTest.h
+	test/ApplyDetailedBalanceTest.h
 	test/ApplyTransmissionCorrectionTest.h
 	test/BinaryOperationTest.h
 	test/BlendSqTest.h

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

+#ifndef MANTID_ALGORITHMS_APPLYDETAILEDBALANCE_H_
+#define MANTID_ALGORITHMS_APPLYDETAILEDBALANCE_H_
+    
+#include "MantidKernel/System.h"
+#include "MantidAPI/Algorithm.h" 
+
+namespace Mantid
+{
+namespace Algorithms
+{
+
+  /** ApplyDetailedBalance : The algorithm calculates the imaginary part of the dynamic susceptibility chi''=Pi*(1-exp(-E/T))*S
+    
+    @author Andrei Savici, ORNL
+    @date 2011-09-01
+
+    Copyright © 2011 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 <http://www.gnu.org/licenses/>.
+
+    File change history is stored at: <https://svn.mantidproject.org/mantid/trunk/Code/Mantid>
+    Code Documentation is available at: <http://doxygen.mantidproject.org>
+  */
+  class DLLExport ApplyDetailedBalance  : public API::Algorithm
+  {
+  public:
+    ApplyDetailedBalance();
+    ~ApplyDetailedBalance();
+    
+    /// Algorithm's name for identification 
+    virtual const std::string name() const { return "ApplyDetailedBalance";};
+    /// Algorithm's version for identification 
+    virtual int version() const { return 1;};
+    /// Algorithm's category for identification
+    virtual const std::string category() const { return "General";}
+    
+  private:
+    /// Sets documentation strings for this algorithm
+    virtual void initDocs();
+    /// Initialise the properties
+    void init();
+    /// Run the algorithm
+    void exec();
+
+
+  };
+
+
+} // namespace Algorithms
+} // namespace Mantid
+
+#endif  /* MANTID_ALGORITHMS_APPLYDETAILEDBALANCE_H_ */

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

+#include "MantidAlgorithms/ApplyDetailedBalance.h"
+#include "MantidKernel/System.h"
+#include "MantidKernel/TimeSeriesProperty.h"
+#include "MantidKernel/PropertyWithValue.h"
+#include "boost/lexical_cast.hpp"
+#include <iostream>
+#include <cmath>
+#include "MantidAPI/WorkspaceValidators.h"
+#include "MantidKernel/PhysicalConstants.h"
+
+using std::string;
+using namespace Mantid::Kernel;
+using namespace Mantid::API;
+
+namespace Mantid
+{
+namespace Algorithms
+{
+
+  // Register the algorithm into the AlgorithmFactory
+  DECLARE_ALGORITHM(ApplyDetailedBalance)
+  
+
+
+  //----------------------------------------------------------------------------------------------
+  /** Constructor
+   */
+  ApplyDetailedBalance::ApplyDetailedBalance()
+  {
+    // TODO Auto-generated constructor stub
+  }
+    
+  //----------------------------------------------------------------------------------------------
+  /** Destructor
+   */
+  ApplyDetailedBalance::~ApplyDetailedBalance()
+  {
+    // TODO Auto-generated destructor stub
+  }
+  
+
+  //----------------------------------------------------------------------------------------------
+  /// Sets documentation strings for this algorithm
+  void ApplyDetailedBalance::initDocs()
+  {
+    this->setWikiSummary("TODO: Enter a quick description of your algorithm.");
+    this->setOptionalMessage("TODO: Enter a quick description of your algorithm.");
+    this->setWikiDescription("TODO: Enter the text to be placed in the Description section of the wiki page.");
+  }
+
+  //----------------------------------------------------------------------------------------------
+  /** Initialize the algorithm's properties.
+   */
+  void ApplyDetailedBalance::init()
+  {
+    CompositeWorkspaceValidator<> *wsValidator = new CompositeWorkspaceValidator<>;
+    wsValidator->add(new WorkspaceUnitValidator<>("DeltaE"));
+    declareProperty(new WorkspaceProperty<>("InputWorkspace","",Direction::Input,wsValidator), "An input workspace.");
+    declareProperty(new WorkspaceProperty<>("OutputWorkspace","",Direction::Output), "An output workspace.");
+    declareProperty(new PropertyWithValue<string>("Temperature","",Direction::Input),"SampleLog variable name that contains the temperature or a number");
+  }
+
+  //----------------------------------------------------------------------------------------------
+  /** Execute the algorithm.
+   */
+  void ApplyDetailedBalance::exec()
+  {
+     MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
+     MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
+     // If input and output workspaces are not the same, create a new workspace for the output
+     if (outputWS != inputWS)
+     {
+       outputWS = API::WorkspaceFactory::Instance().create(inputWS);
+     }
+
+     std::string Tstring=getProperty("Temperature");
+     double Temp;
+     if (inputWS->run().hasProperty(Tstring))
+       Temp=(dynamic_cast<Kernel::TimeSeriesProperty<double> *>(inputWS->run().getProperty(Tstring)))->getStatistics().mean;
+     else
+       Temp=boost::lexical_cast<double>(Tstring);
+     double oneOverT=1.0/(Temp*PhysicalConstants::meVtoKelvin);
+     // Run the exponential correction algorithm explicitly to enable progress reporting
+     IAlgorithm_sptr expcor = createSubAlgorithm("OneMinusExponentialCor",0.0,1.0);
+     expcor->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputWS);
+     expcor->setProperty<MatrixWorkspace_sptr>("OutputWorkspace", outputWS);
+     expcor->setProperty<double>("C1", M_PI);
+     expcor->setProperty<double>("C", oneOverT);
+     expcor->setPropertyValue("Operation","Multiply");
+     expcor->executeAsSubAlg();
+     // Get back the result
+     outputWS = expcor->getProperty("OutputWorkspace");
+
+     setProperty("OutputWorkspace",outputWS);
+  }
+
+
+
+} // namespace Mantid
+} // namespace Algorithms
+

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

+#ifndef MANTID_ALGORITHMS_APPLYDETAILEDBALANCETEST_H_
+#define MANTID_ALGORITHMS_APPLYDETAILEDBALANCETEST_H_
+
+#include <cxxtest/TestSuite.h>
+#include "MantidKernel/Timer.h"
+#include "MantidKernel/System.h"
+#include <iostream>
+#include <iomanip>
+#include "MantidAPI/AlgorithmManager.h"
+#include "MantidTestHelpers/WorkspaceCreationHelper.h"
+#include "MantidDataObjects/Workspace2D.h"
+#include "MantidAlgorithms/ApplyDetailedBalance.h"
+#include "MantidKernel/UnitFactory.h"
+
+using namespace Mantid;
+using namespace Mantid::Algorithms;
+using namespace Mantid::API;
+using namespace Mantid::DataObjects;
+using namespace Mantid::Kernel;
+//using namespace Mantid::DataHandling;
+
+class ApplyDetailedBalanceTest : public CxxTest::TestSuite
+{
+public:
+  // This pair of boilerplate methods prevent the suite being created statically
+  // This means the constructor isn't called when running other tests
+  static ApplyDetailedBalanceTest *createSuite() { return new ApplyDetailedBalanceTest(); }
+  static void destroySuite( ApplyDetailedBalanceTest *suite ) { delete suite; }
+  ApplyDetailedBalanceTest():inputWSname("testInput"),outputWSname("testOutput")
+  {
+  }
+
+  void test_Init()
+  {
+    TS_ASSERT_THROWS_NOTHING( alg.initialize() )
+    TS_ASSERT( alg.isInitialized() )
+  }
+
+
+  void test_exec()
+  {
+    createWorkspace2D(true);
+    TS_ASSERT_THROWS_NOTHING( alg.initialize() )
+    TS_ASSERT( alg.isInitialized() )
+    TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("InputWorkspace",inputWSname) );
+    TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace",outputWSname) );
+    TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("Temperature", "300.") );
+    TS_ASSERT_THROWS_NOTHING( alg.execute() );
+    TS_ASSERT( alg.isExecuted() );
+    /*
+    // Retrieve the workspace from data service. TODO: Change to your desired type
+    Workspace_sptr ws;
+    TS_ASSERT_THROWS_NOTHING( ws = boost::dynamic_pointer_cast<Workspace>(AnalysisDataService::Instance().retrieve(outWSName)) );
+    TS_ASSERT(ws);
+    if (!ws) return;
+    
+    // TODO: Check the results
+    
+    // Remove workspace from the data service.
+    AnalysisDataService::Instance().remove(outWSName);*/
+  }
+  
+  void test_Something()
+  {
+  }
+
+private:
+  ApplyDetailedBalance alg;
+  std::string inputWSname;
+  std::string outputWSname;
+
+  void createWorkspace2D(bool isHistogram)
+  {
+    const int nspecs(1);
+    const int nbins(5);
+    double h=0;
+
+    if(isHistogram) h=0.5;
+
+    Workspace2D_sptr ws2D(new Workspace2D);
+    ws2D->initialize(nspecs,nbins+1,nbins);
+    ws2D->getAxis(0)->unit() = UnitFactory::Instance().create("DeltaE");
+
+    Mantid::MantidVecPtr xv,yv,ev;
+    if (isHistogram)
+      {xv.access().resize(nbins + 1, 0.0);}
+    else
+      {xv.access().resize(nbins , 0.0);}
+    yv.access().resize(nbins, 0.0);
+    ev.access().resize(nbins, 0.0);
+    for (int i = 0; i < nbins; ++i)
+    {
+      xv.access()[i] = static_cast<double>((i-2.-h)*5.);
+      yv.access()[i] = 1.0+i;
+      ev.access()[i] = std::sqrt(1.0+i);
+    }
+    if (isHistogram)
+      {xv.access()[nbins] = static_cast<double>((nbins-2.5)*5.);}
+
+
+    for (int i=0; i< nspecs; i++)
+    {
+      ws2D->setX(i,xv);
+      ws2D->setData(i,yv,ev);
+      ws2D->getSpectrum(i)->setSpectrumNo(i);
+    }
+    ws2D->generateSpectraMap();
+
+    AnalysisDataService::Instance().add(inputWSname, ws2D);
+
+  }
+
+};
+
+
+#endif /* MANTID_ALGORITHMS_APPLYDETAILEDBALANCETEST_H_ */
+

Code/Mantid/Framework/Kernel/inc/MantidKernel/PhysicalConstants.h

@@ -58,8 +58,12 @@ namespace PhysicalConstants
   
   /** 1 meV in wavenumber (cm<SUP>-1</SUP>). Taken from <http://physics.nist.gov/cuu/Constants> on 02/04/2008. */
   static const double meVtoWavenumber = 8.06554465;
-  /**  K-neutron[m^-10] = sqrt(E_neutron[meV]/E_transtormation[mEv]); [G.L. Squires, 1978] may be 2.078 is more accurate?*/
-  static const double E_mev_toNeutronWavenumberSq=2.072; //[meV*Angstrom^2]
+
+  /** 1 meV in Kelvin. Taken from <http://physics.nist.gov/cuu/Constants> on 09/09/2011. */
+  static const double meVtoKelvin = 11.604519;
+
+  /**  K-neutron[m^-10] = sqrt(E_neutron[meV]/E_transtormation[mEv]);*/
+  static const double E_mev_toNeutronWavenumberSq=1.0e20*h_bar*h_bar/(2*NeutronMass*meV); //[meV*Angstrom^2]
 
   /**  Muon lifetime. Taken from MuLan experiment on 08/12/2008. */
   static const double MuonLifetime = 2.197019e-6;