Mantid
Loading...
Searching...
No Matches
CalculateTransmissionBeamSpreader.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
7//----------------------------------------------------------------------
8// Includes
9//----------------------------------------------------------------------
19#include "MantidHistogramData/Histogram.h"
23
24namespace Mantid::Algorithms {
25
26// Register the algorithm into the AlgorithmFactory
27DECLARE_ALGORITHM(CalculateTransmissionBeamSpreader)
28
29using namespace Kernel;
30using namespace API;
31using namespace DataObjects;
32using namespace HistogramData;
33using std::size_t;
34
36 auto wsValidator = std::make_shared<CompositeValidator>();
37 wsValidator->add<WorkspaceUnitValidator>("Wavelength");
38 wsValidator->add<CommonBinsValidator>();
39 wsValidator->add<HistogramValidator>();
40
42 std::make_unique<WorkspaceProperty<>>("SampleSpreaderRunWorkspace", "", Direction::Input, wsValidator),
43 "The workspace containing the sample beam-spreader run");
45 std::make_unique<WorkspaceProperty<>>("DirectSpreaderRunWorkspace", "", Direction::Input, wsValidator),
46 "The workspace containing the direct beam-spreader run");
47 declareProperty(std::make_unique<WorkspaceProperty<>>("SampleScatterRunWorkspace", "", Direction::Input, wsValidator),
48 "The workspace containing the sample scattering run");
49 declareProperty(std::make_unique<WorkspaceProperty<>>("DirectScatterRunWorkspace", "", Direction::Input, wsValidator),
50 "The workspace containing the direct beam scattering run");
51 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
52 "The fitted transmission correction");
53
54 auto zeroOrMore = std::make_shared<BoundedValidator<int>>();
55 zeroOrMore->setLower(0);
56 // The defaults here are the correct detector numbers for LOQ
57 declareProperty("IncidentBeamMonitor", 2, zeroOrMore, "The UDET of the incident beam monitor");
58
59 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
60 mustBePositive->setLower(0.0);
61
62 declareProperty("SpreaderTransmissionValue", 1.0, mustBePositive, "Transmission coefficient of the beam spreader");
63 declareProperty("SpreaderTransmissionError", 0.0, mustBePositive,
64 "Uncertainty on the transmission coefficient of the beam spreader");
65
66 declareProperty("MinWavelength", 2.2, mustBePositive, "The minimum wavelength for the fit");
67 declareProperty("MaxWavelength", 10.0, mustBePositive, "The maximum wavelength for the fit");
68
69 std::vector<std::string> options(2);
70 options[0] = "Linear";
71 options[1] = "Log";
72 declareProperty("FitMethod", "Log", std::make_shared<StringListValidator>(options),
73 "Whether to fit directly to the transmission curve (Linear) "
74 "or to the log of it (Log)");
75
76 declareProperty("OutputUnfittedData", false);
77}
78
80 MatrixWorkspace_sptr sample_spreaderWS = getProperty("SampleSpreaderRunWorkspace");
81 MatrixWorkspace_sptr direct_spreaderWS = getProperty("DirectSpreaderRunWorkspace");
82 MatrixWorkspace_sptr sample_scatterWS = getProperty("SampleScatterRunWorkspace");
83 MatrixWorkspace_sptr direct_scatterWS = getProperty("DirectScatterRunWorkspace");
84
85 // Check that the two input workspaces are from the same instrument
86 if (sample_spreaderWS->getInstrument()->getName() != direct_spreaderWS->getInstrument()->getName() ||
87 sample_spreaderWS->getInstrument()->getName() != sample_scatterWS->getInstrument()->getName() ||
88 sample_spreaderWS->getInstrument()->getName() != direct_scatterWS->getInstrument()->getName()) {
89 g_log.error("The input workspaces do not come from the same instrument");
90 throw std::invalid_argument("The input workspaces do not come from the same instrument");
91 }
92 // Check that the two inputs have matching binning
93 if (!WorkspaceHelpers::matchingBins(*sample_spreaderWS, *direct_spreaderWS) ||
94 !WorkspaceHelpers::matchingBins(*sample_spreaderWS, *sample_scatterWS) ||
95 !WorkspaceHelpers::matchingBins(*sample_spreaderWS, *direct_scatterWS)) {
96 g_log.error("Input workspaces do not have matching binning");
97 throw std::invalid_argument("Input workspaces do not have matching binning");
98 }
99
100 // Extract the required spectra into separate workspaces
101 // The static_cast should not be necessary but it is required to avoid a
102 // "internal compiler error: segmentation fault" when compiling with gcc
103 // and std=c++1z
104 std::vector<detid_t> udets{static_cast<detid_t>(getProperty("IncidentBeamMonitor"))};
105
106 // Convert UDETs to workspace indices
107 // Get monitors (assume that the detector mapping is the same for all data
108 // sets)
109 std::vector<size_t> indices = sample_scatterWS->getIndicesFromDetectorIDs(udets);
110 if (indices.size() != 1) {
111 g_log.error() << "Could not find the incident monitor spectra\n";
112 throw std::invalid_argument("Could not find the incident monitor spectra\n");
113 }
114
115 MatrixWorkspace_sptr sample_scatter_mon = this->extractSpectrum(sample_scatterWS, indices[0]);
116 MatrixWorkspace_sptr direct_scatter_mon = this->extractSpectrum(direct_scatterWS, indices[0]);
117 MatrixWorkspace_sptr sample_spreader_mon = this->extractSpectrum(sample_spreaderWS, indices[0]);
118 MatrixWorkspace_sptr direct_spreader_mon = this->extractSpectrum(direct_spreaderWS, indices[0]);
119
120 // Sum the whole detector for each of the four data sets
121 MatrixWorkspace_sptr sample_scatter_sum;
122 MatrixWorkspace_sptr direct_scatter_sum;
123 MatrixWorkspace_sptr sample_spreader_sum;
124 MatrixWorkspace_sptr direct_spreader_sum;
125
126 // Note: Replaced PARALLEL_SECTION with this OMP for loop, due to occasional
127 // unexplained segfault.
128 std::vector<MatrixWorkspace_sptr> in_ws{sample_scatterWS, direct_scatterWS, sample_spreaderWS, direct_spreaderWS};
129
130 std::vector<MatrixWorkspace_sptr> out_ws(4);
131
132 PARALLEL_FOR_IF(true)
133 for (int i = 0; i < 4; i++) {
134 out_ws[i] = this->sumSpectra(in_ws[i]);
135 }
136 sample_scatter_sum = out_ws[0];
137 direct_scatter_sum = out_ws[1];
138 sample_spreader_sum = out_ws[2];
139 direct_spreader_sum = out_ws[3];
140
141 // Beam spreader transmission
142 MatrixWorkspace_sptr spreader_trans = create<WorkspaceSingleValue>(1, Points(1));
143 spreader_trans->setYUnit("");
144 spreader_trans->setDistribution(true);
145 spreader_trans->mutableX(0)[0] = 0.0;
146 spreader_trans->mutableY(0)[0] = getProperty("SpreaderTransmissionValue");
147 spreader_trans->mutableE(0)[0] = getProperty("SpreaderTransmissionError");
148
149 // The main calculation
150 MatrixWorkspace_sptr numerator =
151 sample_spreader_sum / sample_spreader_mon - spreader_trans * sample_scatter_sum / sample_scatter_mon;
152
153 MatrixWorkspace_sptr denominator =
154 direct_spreader_sum / direct_spreader_mon - spreader_trans * direct_scatter_sum / direct_scatter_mon;
155
156 MatrixWorkspace_sptr transmission = numerator / denominator;
157
158 // Output this data if requested
159 const bool outputRaw = getProperty("OutputUnfittedData");
160 if (outputRaw) {
161 std::string outputWSName = getPropertyValue("OutputWorkspace");
162 outputWSName += "_unfitted";
163 declareProperty(std::make_unique<WorkspaceProperty<>>("UnfittedData", outputWSName, Direction::Output));
164 setProperty("UnfittedData", transmission);
165 }
166
167 // Check that there are more than a single bin in the transmission
168 // workspace. Skip the fit it there isn't.
169 if (transmission->y(0).size() == 1) {
170 setProperty("OutputWorkspace", transmission);
171 } else {
173 const std::string fitMethod = getProperty("FitMethod");
174 logFit = (fitMethod == "Log");
175 if (logFit) {
176 g_log.debug("Fitting to the logarithm of the transmission");
177 // Take a copy of this workspace for the fitting
178 MatrixWorkspace_sptr logTransmission = this->extractSpectrum(transmission, 0);
179
180 // Take the log of each datapoint for fitting. Preserve errors
181 // percentage-wise.
182 auto &Y = logTransmission->mutableY(0);
183 auto &E = logTransmission->mutableE(0);
184 Progress progress(this, 0.4, 0.6, Y.size());
185 for (size_t i = 0; i < Y.size(); ++i) {
186 E[i] = std::abs(E[i] / Y[i]);
187 Y[i] = std::log10(Y[i]);
188 progress.report("Calculate Transmission");
189 }
190
191 // Now fit this to a straight line
192 fit = this->fitToData(logTransmission);
193 } // logFit true
194 else {
195 g_log.debug("Fitting directly to the data (i.e. linearly)");
196 fit = this->fitToData(transmission);
197 }
198
199 setProperty("OutputWorkspace", fit);
200 }
201}
202
208 Algorithm_sptr childAlg = createChildAlgorithm("SumSpectra");
209 childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", WS);
210 childAlg->setProperty<bool>("IncludeMonitors", false);
211 childAlg->executeAsChildAlg();
212 return childAlg->getProperty("OutputWorkspace");
213}
214
222 const size_t index) {
223 // Check that given spectra are monitors
224 if (!WS->spectrumInfo().isMonitor(index)) {
225 g_log.information("The Incident Beam Monitor UDET provided is not marked as a monitor");
226 }
227
228 Algorithm_sptr childAlg = createChildAlgorithm("ExtractSingleSpectrum", 0.0, 0.4);
229 childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", WS);
230 childAlg->setProperty<int>("WorkspaceIndex", static_cast<int>(index));
231 childAlg->executeAsChildAlg();
232 return childAlg->getProperty("OutputWorkspace");
233}
234
241 g_log.information("Fitting the experimental transmission curve");
242 Algorithm_sptr childAlg = createChildAlgorithm("Linear", 0.6, 1.0);
243 childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", WS);
244 const double lambdaMin = getProperty("MinWavelength");
245 const double lambdaMax = getProperty("MaxWavelength");
246 childAlg->setProperty<double>("StartX", lambdaMin);
247 childAlg->setProperty<double>("EndX", lambdaMax);
248 childAlg->executeAsChildAlg();
249
250 std::string fitStatus = childAlg->getProperty("FitStatus");
251 if (fitStatus != "success") {
252 g_log.error("Unable to successfully fit the data: " + fitStatus);
253 throw std::runtime_error("Unable to successfully fit the data");
254 }
255
256 // Only get to here if successful
257 MatrixWorkspace_sptr result = childAlg->getProperty("OutputWorkspace");
258
259 if (logFit) {
260 // Need to transform back to 'unlogged'
261 double b = childAlg->getProperty("FitIntercept");
262 double m = childAlg->getProperty("FitSlope");
263 b = std::pow(10, b);
264 m = std::pow(10, m);
265
266 auto X = result->points(0);
267 auto &Y = result->mutableY(0);
268 auto &E = result->mutableE(0);
269 for (size_t i = 0; i < Y.size(); ++i) {
270 Y[i] = b * (std::pow(m, X[i]));
271 E[i] = std::abs(E[i] * Y[i]);
272 }
273 }
274
275 return result;
276}
277
278} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
A validator which provides a TENTATIVE check that a workspace contains common bins in each spectrum.
A validator which checks that a workspace contains histogram data (the default) or point data as requ...
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
API::MatrixWorkspace_sptr fitToData(const API::MatrixWorkspace_sptr &WS)
Call the Linear fitting algorithm as a child algorithm.
API::MatrixWorkspace_sptr extractSpectrum(const API::MatrixWorkspace_sptr &WS, const size_t index)
Pull out a single spectrum from a 2D workspace.
bool logFit
If true, will take log of transmission curve before fitting.
API::MatrixWorkspace_sptr sumSpectra(const API::MatrixWorkspace_sptr &WS)
Sum the total detector, excluding masked pixels and monitors.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
std::shared_ptr< Algorithm > Algorithm_sptr
Typedef for a shared pointer to an Algorithm.
Definition: Algorithm.h:61
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
static bool matchingBins(const MatrixWorkspace &ws1, const MatrixWorkspace &ws2, const bool firstOnly=false)
Checks whether the bins (X values) of two workspace are the same.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54