Re #10499 Updating code according to RRF paper

Code/Mantid/Framework/Algorithms/CMakeLists.txt

@@ -203,6 +203,7 @@ set ( SRC_FILES
 	src/ResetNegatives.cpp
 	src/ResizeRectangularDetector.cpp
 	src/RingProfile.cpp
+	src/RRFMuon.cpp
 	src/SANSDirectBeamScaling.cpp
 	src/SassenaFFT.cpp
 	src/SaveGSASInstrumentFile.cpp
@@ -456,6 +457,7 @@ set ( INC_FILES
 	inc/MantidAlgorithms/ResetNegatives.h
 	inc/MantidAlgorithms/ResizeRectangularDetector.h
 	inc/MantidAlgorithms/RingProfile.h
+	inc/MantidAlgorithms/RRFMuon.h
 	inc/MantidAlgorithms/SANSDirectBeamScaling.h
 	inc/MantidAlgorithms/SassenaFFT.h
 	inc/MantidAlgorithms/SaveGSASInstrumentFile.h
@@ -703,6 +705,7 @@ set ( TEST_FILES
 	ResetNegativesTest.h
 	ResizeRectangularDetectorTest.h
 	RingProfileTest.h
+	RRFMuonTest.h
 	SassenaFFTTest.h
 	SaveGSASInstrumentFileTest.h
 	ScaleTest.h

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

@@ -0,0 +1,77 @@
+#ifndef MANTID_ALGORITHM_RRFMUON_H_
+#define MANTID_ALGORITHM_RRFMUON_H_
+
+//----------------------------------------------------------------------
+// Includes
+//----------------------------------------------------------------------
+#include "MantidAPI/Algorithm.h"
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_fft_complex.h>
+
+namespace Mantid
+{
+  namespace Algorithms
+  {
+    /**Algorithm for calculating Muon decay envelope.
+
+    @author Raquel Alvarez, ISIS, RAL
+    @date 5/12/2014
+
+    Copyright &copy; 2014-12 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://github.com/mantidproject/mantid>
+    Code Documentation is available at: <http://doxygen.mantidproject.org>
+    */
+    class DLLExport RRFMuon : public API::Algorithm
+    {
+    public:
+      /// Default constructor
+      RRFMuon() {};
+      /// Destructor
+      virtual ~RRFMuon() {};
+      /// Algorithm's name for identification overriding a virtual method
+      virtual const std::string name() const { return "RRFMuon";}
+      ///Summary of algorithm's purpose
+      virtual const std::string summary() const {return "Calculate Muon polarization in the rotating reference frame (RRF).";}
+
+      /// Algorithm's version for identification overriding a virtual method
+      virtual int version() const { return 1;}
+      /// Algorithm's category for identification overriding a virtual method
+      virtual const std::string category() const { return "Muon";}
+
+    private:
+
+      /// Initialise the properties
+      void init();
+      /// Run the algorithm
+      void exec();
+
+      /// Frequency of the oscillations
+      double m_freq;
+      /// Phase accounting for any misalignment of the counters
+      double m_phase;
+      /// Number of data points per histogram
+      size_t m_nData;
+      /// Number of histograms
+      size_t m_nHisto;
+    };
+
+  } // namespace Algorithms
+} // namespace Mantid
+
+#endif /*MANTID_ALGORITHM_RRFMUON_H_*/

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

@@ -0,0 +1,96 @@
+//----------------------------------------------------------------------
+// Includes
+//----------------------------------------------------------------------
+#include "MantidAlgorithms/RRFMuon.h"
+#include "MantidKernel/ArrayProperty.h"
+#include <boost/shared_array.hpp>
+
+#define REAL(z,i) ((z)[2*(i)])
+#define IMAG(z,i) ((z)[2*(i)+1])
+
+namespace Mantid
+{
+namespace Algorithms
+{
+
+using namespace Kernel;
+using API::Progress;
+
+// Register the class into the algorithm factory
+DECLARE_ALGORITHM(RRFMuon)
+
+/** Initialisation method. Declares properties to be used in algorithm.
+ *
+ */
+void RRFMuon::init()
+{
+
+  declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>("InputWorkspace", "", Direction::Input), 
+    "Name of the input workspace containing the spectra in the lab frame");
+
+  declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>("OutputWorkspace", "", Direction::Output), 
+    "Name of the output workspace containing the spectra in the RRF" );
+
+  declareProperty(new PropertyWithValue<double>("Frequency", 0, Direction::Input), 
+    "Frequency of the oscillations");
+
+  declareProperty(new PropertyWithValue<double>("Phase", 0, Direction::Input), 
+    "Phase accounting for any misalignment of the counters");
+
+}
+
+/** Executes the algorithm
+ *
+ */
+void RRFMuon::exec()
+{
+	// Get input workspace containing polarization in the lab-frame
+  API::MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace");
+  // Get frequency
+  m_freq = getProperty("Frequency");
+  // Get phase
+  m_phase = getProperty("Phase");
+  // Get number of histograms
+  m_nHisto = inputWs->getNumberHistograms();
+  if ( m_nHisto != 2 )
+  {
+    throw std::runtime_error("Invalid number of spectra in input workspace");
+  }
+  // Set number of data points
+  m_nData = inputWs->blocksize();
+
+  // Compute the RRF polarization
+  const double twoPiFreq = 2. * M_PI * m_freq;
+  MantidVec time  = inputWs->readX(0); // X axis: time
+  MantidVec labRe = inputWs->readY(0); // Lab-frame polarization (real part)
+  MantidVec labIm = inputWs->readY(1); // Lab-frame polarization (imaginary part)
+  MantidVec rrfRe(m_nData), rrfIm(m_nData); // Rotating Reference frame (RRF) polarizations
+  for (size_t t=0; t<m_nData; t++)
+  {
+    rrfRe [t] =  labRe [t] * cos(twoPiFreq * time[t] + m_phase) + labIm [t] * sin(twoPiFreq * time[t] + m_phase);
+    rrfIm [t] = -labRe [t] * sin(twoPiFreq * time[t] + m_phase) + labIm [t] * cos(twoPiFreq * time[t] + m_phase);
+  }
+
+  // Create output workspace to put results into
+  API::MatrixWorkspace_sptr outputWs = boost::dynamic_pointer_cast<API::MatrixWorkspace>(
+    API::WorkspaceFactory::Instance().create("Workspace2D", m_nHisto, m_nData+1, m_nData));
+  outputWs->getAxis(0)->unit() = inputWs->getAxis(0)->unit();
+
+  // Put results into output workspace
+  // Real RRF polarization
+  outputWs->getSpectrum(0)->setSpectrumNo(1);
+  outputWs->dataX(0) = inputWs->readX(0);
+  outputWs->dataY(0) = rrfRe;
+  // Imaginary RRF polarization
+  outputWs->getSpectrum(1)->setSpectrumNo(2);
+  outputWs->dataX(1) = inputWs->readX(1);
+  outputWs->dataY(1) = rrfIm;
+
+  // Set output workspace
+  setProperty("OutputWorkspace", outputWs);
+
+}
+
+
+}
+}

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

@@ -0,0 +1,127 @@
+#ifndef MANTID_ALGORITHMS_RRFMUON_H_
+#define MANTID_ALGORITHMS_RRFMUON_H_
+
+#include <cxxtest/TestSuite.h>
+#include "MantidAlgorithms/RRFMuon.h"
+#include "MantidDataObjects/Workspace2D.h"
+#include "MantidAPI/Workspace.h"
+#include <Poco/File.h>
+#include <stdexcept>
+
+using namespace Mantid::Algorithms;
+using namespace Mantid::DataObjects;
+using namespace Mantid::API;
+
+class RRFMuonTest : public CxxTest::TestSuite
+{
+public:
+
+  void testName()
+  {
+    TS_ASSERT_EQUALS( rrfMuon.name(), "RRFMuon" );
+  }
+
+  void testCategory()
+  {
+    TS_ASSERT_EQUALS( rrfMuon.category(), "Muon" )
+  }
+
+  void testRRFMuonZeroFrequency()
+  {
+    // Test of the algorithm at zero frequency
+    // At zero frequency input and output workspaces should contain the same X, Y data
+
+    // Create input workspace with three spectra
+    MatrixWorkspace_sptr ws = createDummyWorkspace();
+
+    // Initialise
+    TS_ASSERT_THROWS_NOTHING( rrfMuon.initialize() );
+    TS_ASSERT( rrfMuon.isInitialized() );
+    // Set Values
+    TS_ASSERT_THROWS_NOTHING( rrfMuon.setProperty("InputWorkspace", ws) );
+    TS_ASSERT_THROWS_NOTHING( rrfMuon.setProperty("OutputWorkspace", "outputWs") );
+    TS_ASSERT_THROWS_NOTHING( rrfMuon.setProperty("Frequency", "0") );
+    TS_ASSERT_THROWS_NOTHING( rrfMuon.setProperty("Phase", "0") );
+    // Execute
+    TS_ASSERT_THROWS_NOTHING(rrfMuon.execute());
+    TS_ASSERT(rrfMuon.isExecuted());
+    // Get result
+    MatrixWorkspace_const_sptr ows =
+      boost::dynamic_pointer_cast<MatrixWorkspace>(AnalysisDataService::Instance().retrieve("outputWs"));
+    TS_ASSERT(ows);
+
+    // Checks
+    // X values
+    TS_ASSERT_EQUALS( ws->readX(0), ows->readX(0) );
+    TS_ASSERT_EQUALS( ws->readX(1), ows->readX(1) );
+    // Y values
+    TS_ASSERT_EQUALS( ws->readY(0), ows->readY(0) );
+    TS_ASSERT_EQUALS( ws->readY(1), ows->readY(1) );
+  }
+
+    void testRRFMuonNonZeroFrequency()
+  {
+    // Test of the algorithm at non-zero frequency
+
+    // Create input workspace with three spectra
+    MatrixWorkspace_sptr ws = createDummyWorkspace();
+
+    // Initialise
+    TS_ASSERT_THROWS_NOTHING( rrfMuon.initialize() );
+    TS_ASSERT( rrfMuon.isInitialized() );
+    // Set Values
+    TS_ASSERT_THROWS_NOTHING( rrfMuon.setProperty("InputWorkspace", ws) );
+    TS_ASSERT_THROWS_NOTHING( rrfMuon.setProperty("OutputWorkspace", "outputWs") );
+    TS_ASSERT_THROWS_NOTHING( rrfMuon.setProperty("Frequency", "197") );
+    TS_ASSERT_THROWS_NOTHING( rrfMuon.setProperty("Phase", "0") );
+    // Execute
+    TS_ASSERT_THROWS_NOTHING(rrfMuon.execute());
+    TS_ASSERT(rrfMuon.isExecuted());
+    // Get result
+    MatrixWorkspace_const_sptr ows =
+      boost::dynamic_pointer_cast<MatrixWorkspace>(AnalysisDataService::Instance().retrieve("outputWs"));
+    TS_ASSERT(ows);
+
+    // Checks
+    // X values
+    TS_ASSERT_EQUALS( ws->readX(0), ows->readX(0) );
+    TS_ASSERT_EQUALS( ws->readX(1), ows->readX(1) );
+    // Y values
+    // The input frequency is close to the precession frequency, so:
+    // The real part of the RRF polarization should be close to 1 for all X values
+    // The imaginary part should be close to 0 for all X values
+    TS_ASSERT_DELTA( ows->readY(0)[ 0], 1, 0.001 );
+    TS_ASSERT_DELTA( ows->readY(0)[49], 1, 0.001 );
+    TS_ASSERT_DELTA( ows->readY(0)[99], 1, 0.001 );
+    TS_ASSERT_DELTA( ows->readY(1)[ 0], 0, 0.001 );
+    TS_ASSERT_DELTA( ows->readY(1)[49], 0, 0.001 );
+    TS_ASSERT_DELTA( ows->readY(1)[99], 0, 0.001 );
+  }
+
+private:
+  RRFMuon rrfMuon;
+
+  MatrixWorkspace_sptr createDummyWorkspace()
+  {
+    int nBins = 100;
+    double pi = 3.14159;
+    MatrixWorkspace_sptr ws = WorkspaceFactory::Instance().create("Workspace2D", 2, nBins+1, nBins);
+
+    for (int i=0; i<nBins; i++)
+    {
+      double x = (i-50.)/nBins;
+      ws->dataX(0)[i] = x;
+      ws->dataY(0)[i] = cos(-2*pi*3*x);
+      ws->dataX(1)[i] = x;
+      ws->dataY(1)[i] = sin(-2*pi*3*x);
+    }
+
+    ws->dataX(0)[nBins] = nBins;
+    ws->dataX(1)[nBins] = nBins;
+
+    return ws;
+  }
+
+};
+
+#endif /* MANTID_ALGORITHMS_RRFMUON_H_ */

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

@@ -0,0 +1,50 @@
+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+Assuming *InputWorkspace* contains the real and imaginary parts of the asymmetry in the lab-fram, the algorithm 
+returns the muon polarization in the Rotating Reference Frame (RRF), given the frequency of the oscillations and the phase of the 
+detectors, as given by the formula:
+
+.. math:: P_R^{RRF} (\nu_0, \phi, t)= + P_R\left(0,t\right) \cos\left(2\pi\nu_0 t + \phi\right) + P_I\left(0,t\right) \sin\left(2\pi\nu_0 t + \phi\right)
+.. math:: P_I^{RRF} (\nu_0, \phi, t)= - P_R\left(0,t\right) \sin\left(2\pi\nu_0 t + \phi\right) + P_I\left(0,t\right) \cos\left(2\pi\nu_0 t + \phi\right)
+
+where :math:`P_R\left(0,t\right)` and :math:`P_I\left(0,t\right)` are the real and imaginary part of the asymmetry in the lab-frame, 
+:math:`\nu_0` is the input frequency, and :math:`\phi` the input phase.
+
+Usage
+-----
+
+**Example - Computing polarization in RRF**
+
+.. testcode:: ExRRF
+
+   import math
+   # Create an input workspace with two spectra
+   datax = [(i-50)/100. for i in range(1,100)]
+   datay1 = [ math.cos(2*3.14159*3*(50-i)/100.) for i in range(1,99) ]
+   datay2 = [ math.sin(2*3.14159*3*(50-i)/100.) for i in range(1,99) ]
+   datay = datay1 + datay2
+   input = CreateWorkspace(dataX=datax, dataY=datay,Nspec=2)
+   # Compute polarization in RRF
+   output = RRFMuon(input,196,0)
+   # Print some values
+   print output.readY(0)[49]
+   print output.readY(1)[49]
+
+Output:
+
+.. testoutput:: ExRRF
+
+   1.0
+   0.0
+
+
+.. categories::