diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/RRFMuon.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/RRFMuon.h new file mode 100644 index 000000000000..f05d549d894e --- /dev/null +++ b/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 +#include + +namespace Mantid +{ + namespace Algorithms + { + /**Algorithm for calculating Muon decay envelope. + + @author Raquel Alvarez, ISIS, RAL + @date 5/12/2014 + + Copyright © 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 . + + File change history is stored at: + Code Documentation is available at: + */ + 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 decay envelop from squashograms and input frequency.";} + + /// 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(); + /// Compute the Hilbert transform + void hilbert (API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspace_sptr outputWs); + + /// Frequency of the oscillations + double m_freq; + /// 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_*/ \ No newline at end of file diff --git a/Code/Mantid/Framework/Algorithms/src/RRFMuon.cpp b/Code/Mantid/Framework/Algorithms/src/RRFMuon.cpp new file mode 100644 index 000000000000..8da6f0ced83e --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/src/RRFMuon.cpp @@ -0,0 +1,159 @@ +//---------------------------------------------------------------------- +// Includes +//---------------------------------------------------------------------- +#include "MantidAlgorithms/RRFMuon.h" +#include + +#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("InputWorkspace", "", Direction::Input), + "Name of the input workspace containing the spectra"); + + declareProperty(new API::WorkspaceProperty("OutputWorkspace", "", Direction::Output), + "Name of the output workspace containing envelope of the decay" ); + + declareProperty(new PropertyWithValue("Frequency", 0, Direction::Input), + "Frequency of the oscillations"); + +} + +/** Executes the algorithm + * + */ +void RRFMuon::exec() +{ + // Get input workspace + API::MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace"); + + // Get number of histograms + m_nHisto = inputWs->getNumberHistograms(); + // Get number of data points + m_nData = inputWs->blocksize(); + + //if (m_nHisto!=2) + //{ + // throw std::runtime_error("RRFMuon: InputWorkspace must contain two spectra"); + //} + + + // Create output workspace with two spectra + API::MatrixWorkspace_sptr outputWs = boost::dynamic_pointer_cast( + API::WorkspaceFactory::Instance().create("Workspace2D", m_nHisto, m_nData+1, m_nData)); + outputWs->getAxis(0)->unit() = inputWs->getAxis(0)->unit(); + + hilbert(inputWs,outputWs); + + // Get input frequency + m_freq = getProperty("Frequency"); + + + // Set output workspace + setProperty("OutputWorkspace", outputWs); + +} + +void RRFMuon::hilbert (API::MatrixWorkspace_sptr inputWs, API::MatrixWorkspace_sptr outputWs) +{ + const int ySize = static_cast(inputWs->blocksize()); + + // Alloc memory + gsl_fft_complex_wavetable * wavetable = gsl_fft_complex_wavetable_alloc(ySize); + gsl_fft_complex_workspace * workspace = gsl_fft_complex_workspace_alloc(ySize); + boost::shared_array data(new double[2*ySize]); + + + // Data initialisation + for(int i=0; idataY(0)[i]; // even indexes filled with the real part + data[2*i+1] = 0; // odd indexes filled with the imaginary part + } + + // Forward FFT + gsl_fft_complex_forward (data.get(), 1, ySize, wavetable, workspace); + + // Compute Hilbert transform + boost::shared_array hilbertData(new double[2*ySize]); + for(int i=0; i ySize/2 ) + { + hilbertData[2*i ] = -data[2*i+1]; + hilbertData[2*i+1] = data[2*i]; + } + else + { + hilbertData[2*i ] = 0; + hilbertData[2*i+1] = 0; + } + } + gsl_fft_complex_inverse (hilbertData.get(), 1, ySize, wavetable, workspace); + + + // Write results into output workspace + outputWs->dataX(0) = inputWs->dataX(0); + auto oldData = inputWs->dataY(0); + for (int i=0; idataY(0)[i] = hilbertData[2*i]; // uncommet to print hilbert transform + outputWs->dataY(0)[i] = sqrt( oldData[i]*oldData[i] + hilbertData[2*i]*hilbertData[2*i]); + } + + // Free memory + gsl_fft_complex_wavetable_free (wavetable); + gsl_fft_complex_workspace_free (workspace); +} + + +//void hilbert(double *source, double *target, double *filt, int npt, int lfilt) +//{ +// int i=1; +// int l=1; +// double yt =0.0; +// for (l=1; l<=npt-lfilt+1; l++) +// { +// yt = 0.0; +// for (i=1; i<=lfilt; i++) +// yt = yt+source[l+i-1]*filt[lfilt+1-i]; +// target[l] = yt; +// } +// /* shifting points */ +// for (i=1; i<=npt-lfilt; i++) +// { +// target[i] = 0.5*(target[i]+target[i+1]); +// } +// for (i=npt-lfilt; i>=1; i--) +// target[i+lfilt/2]=target[i]; +// /* adding zeros */ +// for (i=1; i<=lfilt/2; i++) +// { +// target[i] = 0.0; +// target[npt+1-i] = 0.0; +// } +//} + +} +} \ No newline at end of file diff --git a/Code/Mantid/Framework/Algorithms/test/RRFMuonTest.h b/Code/Mantid/Framework/Algorithms/test/RRFMuonTest.h new file mode 100644 index 000000000000..c0bfb87782ad --- /dev/null +++ b/Code/Mantid/Framework/Algorithms/test/RRFMuonTest.h @@ -0,0 +1,37 @@ +#ifndef MANTID_ALGORITHMS_RRFMUON_H_ +#define MANTID_ALGORITHMS_RRFMUON_H_ + +#include +#include "MantidAlgorithms/RRFMuon.h" +#include +#include + +using namespace Mantid::Algorithms; +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 testInit() + { + TS_ASSERT_THROWS_NOTHING( rrfMuon.initialize() ); + TS_ASSERT( rrfMuon.isInitialized() ) + } + + +private: + RRFMuon rrfMuon; +}; + +#endif /* MANTID_ALGORITHMS_RRFMUON_H_ */ \ No newline at end of file diff --git a/Code/Mantid/docs/source/algorithms/RRFMuon-v1.rst b/Code/Mantid/docs/source/algorithms/RRFMuon-v1.rst new file mode 100644 index 000000000000..b932b35bffdb --- /dev/null +++ b/Code/Mantid/docs/source/algorithms/RRFMuon-v1.rst @@ -0,0 +1,15 @@ +.. algorithm:: + +.. summary:: + +.. alias:: + +.. properties:: + +Description +----------- + +Assuming that the *InputWorkspace* contains two squashograms, the algorithm returns a workspace +containing the envelope of the muon decay, for a given frequency. + +.. categories::