Re #10498 Adding cpp h and CMakeList files

Code/Mantid/Framework/Algorithms/CMakeLists.txt

@@ -14,6 +14,7 @@ set ( SRC_FILES
 	src/AppendSpectra.cpp
 	src/ApplyCalibration.cpp
 	src/ApplyDeadTimeCorr.cpp
+	src/PhaseQuadMuon.cpp
 	src/ApplyDetailedBalance.cpp
 	src/ApplyTransmissionCorrection.cpp
 	src/AsymmetryCalc.cpp
@@ -262,6 +263,7 @@ set ( INC_FILES
 	inc/MantidAlgorithms/AppendSpectra.h
 	inc/MantidAlgorithms/ApplyCalibration.h
 	inc/MantidAlgorithms/ApplyDeadTimeCorr.h
+	inc/MantidAlgorithms/PhaseQuadMuon.h
 	inc/MantidAlgorithms/ApplyDetailedBalance.h
 	inc/MantidAlgorithms/ApplyTransmissionCorrection.h
 	inc/MantidAlgorithms/AsymmetryCalc.h

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

@@ -0,0 +1,107 @@
+#ifndef MANTID_ALGORITHM_PHASEQUADMUON_H_
+#define MANTID_ALGORITHM_PHASEQUADMUON_H_
+
+//----------------------------------------------------------------------
+// Includes
+//----------------------------------------------------------------------
+#include "MantidAPI/Algorithm.h"
+
+namespace Mantid
+{
+  namespace Algorithms
+  {
+    /**Algorithm for calculating Muon spectra.
+
+    @author Raquel Alvarez, ISIS, RAL
+    @date 1/12/2011
+
+    Copyright © 2014-11 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 PhaseQuadMuon : public API::Algorithm
+    {
+    public:
+      /// Default constructor
+      PhaseQuadMuon() : API::Algorithm(), m_muLife(2.19703E-06) {};
+      /// Destructor
+      virtual ~PhaseQuadMuon() {};
+      /// Algorithm's name for identification overriding a virtual method
+      virtual const std::string name() const { return "PhaseQuad";}
+      ///Summary of algorithms purpose
+      virtual const std::string summary() const {return "Calculate Muon spectra from PhaseTable.";}
+
+      /// 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:
+
+      class HistData { // TODO: do I need all the members?
+      public:
+        double n0;    // 
+        double alpha; // Efective asymmetry
+        double phi;   // Phase
+        double lag;   // Time-shift
+        double dead;  // 
+        double deadm; //
+        double chisq; // Of silver fit
+      };
+
+      /// Initialise the properties
+      void init();
+      /// Run the algorithm
+      void exec();
+      /// Load the PhaseTable
+      void loadPhaseTable(const std::string& filename);
+      /// TODO
+      void removeExponentialDecay (API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspace_sptr outputWs);
+      /// TODO
+      void squash(API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspace_sptr outputWs);
+      /// Number of histograms (spectra)
+      size_t m_nHist;
+      /// Number of datapoints in each histogram
+      size_t m_nData;
+      /// TODO
+      double m_res;
+      /// Mean of time-shifts
+      double m_meanLag;
+      /// Good muons from here on (bin no)
+      size_t m_tValid;
+      /// Pulse definitely finished by here (bin no)
+      size_t m_tNoPulse;
+      /// Double-pulse flag
+      bool m_isDouble;
+      /// Muon decay curve
+      double m_tau;
+      /// Freq (of silver run)
+      double m_w;
+      /// Mu life
+      double m_muLife;
+      /// TODO
+      double m_bigNumber;
+      /// Vector of histograms data
+      std::vector<HistData> m_histData;
+    };
+
+  } // namespace Algorithms
+} // namespace Mantid
+
+#endif /*MANTID_ALGORITHM_PHASEQUAD_H_*/

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

@@ -0,0 +1,188 @@
+//----------------------------------------------------------------------
+// Includes
+//----------------------------------------------------------------------
+#include "MantidKernel/ArrayProperty.h"
+#include "MantidAlgorithms/PhaseQuadMuon.h"
+#include "MantidAPI/FileProperty.h"
+#include "MantidAPI/IFileLoader.h"
+
+namespace Mantid
+{
+namespace Algorithms
+{
+
+using namespace Kernel;
+using API::Progress;
+
+// Register the class into the algorithm factory
+DECLARE_ALGORITHM(PhaseQuadMuon)
+
+/** Initialisation method. Declares properties to be used in algorithm.
+ *
+ */
+void PhaseQuadMuon::init()
+{
+
+  declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>("InputWorkspace", "", Direction::Input), 
+    "Name of the input workspace");
+
+  declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>("OutputWorkspace","", Direction::Output), 
+    "Name of the output workspace" );
+
+  declareProperty(new API::FileProperty("PhaseTable", "", API::FileProperty::Load, ".INF"),
+    "The name of the list of phases for each spectrum");
+}
+
+/** Executes the algorithm
+ *
+ */
+void PhaseQuadMuon::exec()
+{
+	// TODO
+
+	// Get input workspace
+  API::MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace");
+
+  // Get phase table
+  std::string filename = getPropertyValue("PhaseTable");
+  loadPhaseTable ( filename );
+
+  m_nData = inputWs->getSpectrum(0)->readY().size();
+  m_nHist = inputWs->getNumberHistograms();
+
+  // Create output workspace with two spectra
+  API::MatrixWorkspace_sptr outputWs = getProperty("OutputWorkspace");
+  outputWs = boost::dynamic_pointer_cast<API::MatrixWorkspace>(
+    API::WorkspaceFactory::Instance().create("Workspace2D", inputWs->getNumberHistograms(), inputWs->getSpectrum(0)->readX().size(), inputWs->getSpectrum(0)->readY().size())); // TODO change number histograms
+  outputWs->getAxis(0)->unit() = inputWs->getAxis(0)->unit();
+
+  removeExponentialDecay(inputWs,outputWs);
+  //squash (inputWs, outputWs);
+
+  setProperty("OutputWorkspace", outputWs);
+}
+
+//----------------------------------------------------------------------------------------------
+/** Load PhaseTable file to a vector of HistData.
+* @param filename :: phase table .inf filename
+*/
+void PhaseQuadMuon::loadPhaseTable(const std::string& filename )
+{
+
+  std::ifstream input(filename.c_str(), std::ios_base::in);
+
+  if (input.is_open())
+  {
+
+    if ( input.eof() )
+    {
+      throw Exception::FileError("PhaseQuad: File is empty", filename);
+    }
+    else
+    {
+      std::string line;
+
+      // Header of .INF file is as follows:
+      //
+      // Comment on the output file
+      // Top row of numbers are:
+      // #histos, typ. first good bin#, typ. bin# when pulse over, mean lag.
+      // Tabulated numbers are, per histogram:
+      // det ok, asymmetry, phase, lag, deadtime_c, deadtime_m.
+      //
+      std::getline( input, line ); // Skip first line in header
+      std::getline( input, line ); // Skip second line
+      std::getline( input, line ); // ...
+      std::getline( input, line );
+      std::getline( input, line );
+
+      // Read first useful line
+      input >> m_nHist >> m_tValid >> m_tNoPulse >> m_meanLag;
+
+      int cont=0;
+      while( !input.eof() )
+      {
+        // Read histogram data
+        HistData tempData;
+        input >> tempData.n0 >> tempData.alpha >> tempData.phi >> tempData.lag >> tempData.dead >> tempData.deadm;
+        m_histData.push_back (tempData);
+        cont++;
+      }
+
+      if ( cont!= m_nHist )
+      {
+        throw Exception::FileError("PhaseQuad: File was not in expected format", filename); // TODO: specify that number of lines is wrong
+      }
+
+      if ( cont<3 )
+      {
+        throw std::runtime_error("PhaseQuad: Found less than 4 histograms");
+      }
+
+    }
+  }
+  else
+  {
+    // Throw exception if file cannot be opened
+    throw std::runtime_error("PhaseQuad: Unable to open PhaseTable");
+  }
+
+}
+
+
+//----------------------------------------------------------------------------------------------
+/** Remove exponential decay from input histograms
+* @param inputWs :: input workspace containing the spectra
+* @param outputWs :: output workspace to hold temporary results
+*/
+void PhaseQuadMuon::removeExponentialDecay (API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspace_sptr outputWs)
+{
+
+  for (size_t h=0; h<inputWs->getNumberHistograms(); h++)
+  {
+    auto specIn = inputWs->getSpectrum(h);
+    auto specOut = outputWs->getSpectrum(h);
+
+    for (int i=0; i<m_nData; i++)
+    {
+        double usey = specIn->readY()[i];
+        double oops = ( usey<= 0 ) || ( specIn->readE()[i] > m_bigNumber );
+        specOut->dataY()[i] = oops ? 0 : log(usey);
+        specOut->dataE()[i] = oops ? m_bigNumber : specIn->readE()[i] / usey;
+        std::cout << 
+    }
+    double s=0;
+    double sx=0;
+    double sy=0;
+
+    for (int i=0; i<m_nData; i++)
+    {
+      double sigma = specOut->dataE()[i] * specOut->dataE()[i];
+      s+= 1./sigma;
+      sx+= specOut->dataX()[i]/sigma;
+      sy+= specOut->dataY()[i]/sigma;
+    }
+    double N0 = (sy+sx/m_muLife/1E9)/s;
+//    std::cout << specIn->readY()[40] << " " << sy << " " << sx << " " << s << " " << N0 << std::endl;
+  }
+
+}
+
+
+void PhaseQuadMuon::squash(API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspace_sptr outputWs)
+{
+
+  if ( m_nHist != inputWs->getNumberHistograms() )
+  {
+    throw std::runtime_error("InputWorkspace and PhaseTable do not have the same number of spectra");
+  }
+  else
+  {
+    // TODO
+  }
+
+}
+
+
+}
+}