forked from npshub/mantid
-
Notifications
You must be signed in to change notification settings - Fork 0
/
RRFMuon.cpp
122 lines (103 loc) · 4.46 KB
/
RRFMuon.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
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
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidMuon/RRFMuon.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidKernel/Unit.h"
namespace Mantid::Algorithms {
using namespace Kernel;
// Register the class into the algorithm factory
DECLARE_ALGORITHM(RRFMuon)
/** Initialisation method. Declares properties to be used in algorithm.
*
*/
void RRFMuon::init() {
declareProperty(
std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>("InputWorkspace", "", Direction::Input),
"Name of the input workspace containing the spectra in the lab frame");
declareProperty(
std::make_unique<API::WorkspaceProperty<API::MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
"Name of the output workspace containing the spectra in the RRF");
declareProperty(std::make_unique<PropertyWithValue<double>>("Frequency", 0, Direction::Input),
"Frequency of the oscillations");
std::vector<std::string> unitOptions{"MHz", "Gauss", "Mrad/s"};
declareProperty("FrequencyUnits", "MHz", std::make_shared<StringListValidator>(unitOptions), "The frequency units");
declareProperty(std::make_unique<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
double freq = getProperty("Frequency");
// Get units
std::string units = getProperty("FrequencyUnits");
// Convert frequency to input workspace X units
double factor = unitConversionFactor(inputWs->getAxis(0)->unit()->label().ascii(), units);
// Get phase
double phase = getProperty("Phase");
// Get number of histograms
const size_t nHisto = inputWs->getNumberHistograms();
if (nHisto != 2) {
throw std::runtime_error("Invalid number of spectra in input workspace");
}
// Set number of data points
const size_t nData = inputWs->blocksize();
// Compute the RRF polarization
const double twoPiFreq = 2. * M_PI * freq * factor;
const auto &time = inputWs->x(0); // X axis: time
const auto &labRe = inputWs->y(0); // Lab-frame polarization (real part)
const auto &labIm = inputWs->y(1); // Lab-frame polarization (imaginary part)
MantidVec rrfRe(nData),
rrfIm(nData); // Rotating Reference frame (RRF) polarizations
for (size_t t = 0; t < nData; t++) {
rrfRe[t] = labRe[t] * cos(twoPiFreq * time[t] + phase) + labIm[t] * sin(twoPiFreq * time[t] + phase);
rrfIm[t] = -labRe[t] * sin(twoPiFreq * time[t] + phase) + labIm[t] * cos(twoPiFreq * time[t] + phase);
}
// Create output workspace to put results into
API::MatrixWorkspace_sptr outputWs = std::dynamic_pointer_cast<API::MatrixWorkspace>(
API::WorkspaceFactory::Instance().create("Workspace2D", nHisto, nData + 1, nData));
outputWs->getAxis(0)->unit() = inputWs->getAxis(0)->unit();
// Put results into output workspace
// Real RRF polarization
outputWs->setSharedX(0, inputWs->sharedX(0));
outputWs->mutableY(0) = rrfRe;
// Imaginary RRF polarization
outputWs->setSharedX(1, inputWs->sharedX(1));
outputWs->mutableY(1) = rrfIm;
// Set output workspace
setProperty("OutputWorkspace", outputWs);
}
/** Gets factor to convert frequency units to input workspace units, typically
* in microseconds
* @param uin :: [input] input workspace units
* @param uuser :: [input] units selected by user
*/
double RRFMuon::unitConversionFactor(const std::string &uin, const std::string &uuser) {
if ((uin == "microsecond")) {
if (uuser == "MHz") {
return 1.0;
} else if (uuser == "Gauss") {
// Factor = 2 * PI * MU
// where MU is the muon gyromagnetic ratio:
// MU = 135.538817 MHz/T, 1T = 10000 Gauss
return 2.0 * M_PI * 135.538817 * 0.0001;
} else if (uuser == "Mrad/s") {
// Factor = 2 * PI
return 2.0 * M_PI;
} else {
throw std::runtime_error("Could not find units");
}
} else {
throw std::runtime_error("X units must be in microseconds");
}
}
} // namespace Mantid::Algorithms