Newer
Older
Peterson, Peter
committed
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include <cmath>
#include <vector>
#include "MantidAPI/Workspace.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidAlgorithms/MuonRemoveExpDecay.h"
namespace Mantid
{
namespace Algorithms
{
using namespace Kernel;
using API::Progress;
// Register the class into the algorithm factory
DECLARE_ALGORITHM(MuonRemoveExpDecay)
Janik Zikovsky
committed
/// Sets documentation strings for this algorithm
void MuonRemoveExpDecay::initDocs()
{
this->setWikiSummary("This algorithm removes the exponential decay from a muon workspace. ");
this->setOptionalMessage("This algorithm removes the exponential decay from a muon workspace.");
}
Peterson, Peter
committed
/** Initialisation method. Declares properties to be used in algorithm.
*
*/
void MuonRemoveExpDecay::init()
{
declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>("InputWorkspace", "",
Direction::Input), "Name of the input 2D workspace");
declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>("OutputWorkspace", "",
Direction::Output), "The name of the workspace to be created as the output of the algorithm");
std::vector<int> empty;
declareProperty(new Kernel::ArrayProperty<int>("Spectra", empty), "Spectra to remove the exponential decay from");
Peterson, Peter
committed
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
}
/** Executes the algorithm
*
*/
void MuonRemoveExpDecay::exec()
{
std::vector<int> Spectra = getProperty("Spectra");
//Get original workspace
API::MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
int numSpectra = inputWS->size() / inputWS->blocksize();
//Create output workspace with same dimensions as input
API::MatrixWorkspace_sptr outputWS = API::WorkspaceFactory::Instance().create(inputWS);
if (Spectra.size() == 0)
{
Progress prog(this, 0.0, 1.0, numSpectra);
//Do all the spectra
PARALLEL_FOR2(inputWS,outputWS)
for (int i = 0; i < numSpectra; ++i)
{
PARALLEL_START_INTERUPT_REGION
removeDecay(inputWS->readX(i), inputWS->readY(i), outputWS->dataY(i));
outputWS->dataX(i) = inputWS->readX(i);
//Need to do something about the errors?
prog.report();
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
}
else
{
Progress prog(this, 0.0, 1.0, numSpectra + Spectra.size());
if (getPropertyValue("InputWorkspace") != getPropertyValue("OutputWorkspace"))
{
//Copy all the X,Y and E data
PARALLEL_FOR2(inputWS,outputWS)
for (int i = 0; i < numSpectra; ++i)
{
PARALLEL_START_INTERUPT_REGION
outputWS->dataX(i) = inputWS->readX(i);
outputWS->dataY(i) = inputWS->readY(i);
outputWS->dataE(i) = inputWS->readE(i);
prog.report();
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
}
//Do the specified spectra only
int specLength = Spectra.size();
PARALLEL_FOR2(inputWS,outputWS)
for (int i = 0; i < specLength; ++i)
{
PARALLEL_START_INTERUPT_REGION
if (Spectra[i] > numSpectra)
{
g_log.error("Spectra size greater than the number of spectra!");
throw std::invalid_argument("Spectra size greater than the number of spectra!");
}
removeDecay(inputWS->readX(Spectra[i]), inputWS->readY(Spectra[i]), outputWS->dataY(Spectra[i]));
outputWS->dataX(Spectra[i]) = inputWS->readX(Spectra[i]);
//Need to do something about the errors?
prog.report();
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
}
setProperty("OutputWorkspace", outputWS);
}
/** This method corrects the data for one spectra.
* The muon lifetime is in microseconds not seconds, i.e. 2.2 rather than 0.0000022.
* This is because the data is in microseconds.
Janik Zikovsky
committed
* @param inX :: The X vector
* @param inY :: The input data vector
* @param outY :: The output data vector
Peterson, Peter
committed
*/
void MuonRemoveExpDecay::removeDecay(const MantidVec& inX, const MantidVec& inY,
MantidVec& outY)
{
//Do the removal
for (size_t i = 0; i < inY.size(); ++i)
{
outY[i] = inY[i] * exp(inX[i] / (Mantid::PhysicalConstants::MuonLifetime * 1000000.0));
}
}
} // namespace Algorithm
} // namespace Mantid