-
Notifications
You must be signed in to change notification settings - Fork 122
/
NormaliseByDetector.cpp
202 lines (173 loc) · 8.34 KB
/
NormaliseByDetector.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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
// 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 "MantidAlgorithms/NormaliseByDetector.h"
#include "MantidAPI/FunctionDomain1D.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/FunctionValues.h"
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/IFunction.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidGeometry/Instrument/Component.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
#include "MantidGeometry/Instrument/FitParameter.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidGeometry/muParser_Silent.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/UnitFactory.h"
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
namespace Mantid::Algorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(NormaliseByDetector)
//----------------------------------------------------------------------------------------------
/** Constructor
*/
NormaliseByDetector::NormaliseByDetector(bool parallelExecution) : m_parallelExecution(parallelExecution) {}
//----------------------------------------------------------------------------------------------
/// Algorithm's name for identification. @see Algorithm::name
const std::string NormaliseByDetector::name() const { return "NormaliseByDetector"; }
/// Algorithm's version for identification. @see Algorithm::version
int NormaliseByDetector::version() const { return 1; }
/// Algorithm's category for identification. @see Algorithm::category
const std::string NormaliseByDetector::category() const { return "CorrectionFunctions\\NormalisationCorrections"; }
//----------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void NormaliseByDetector::init() {
auto compositeValidator = std::make_shared<CompositeValidator>();
compositeValidator->add(std::make_shared<API::WorkspaceUnitValidator>("Wavelength"));
compositeValidator->add(std::make_shared<API::HistogramValidator>());
declareProperty(
std::make_unique<WorkspaceProperty<MatrixWorkspace>>("InputWorkspace", "", Direction::Input, compositeValidator),
"An input workspace in wavelength");
declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
"An output workspace.");
}
const Geometry::FitParameter NormaliseByDetector::tryParseFunctionParameter(const Geometry::Parameter_sptr ¶meter,
const Geometry::IDetector &det) {
if (parameter == nullptr) {
std::stringstream stream;
stream << det.getName()
<< " and all of it's parent components, have no "
"fitting type parameters. This algorithm "
"cannot be run without fitting parameters. See "
"wiki help for details on setup.";
this->g_log.warning(stream.str());
throw std::runtime_error(stream.str());
}
return parameter->value<Geometry::FitParameter>();
}
/**
Process each histogram of the input workspace, extracting the detector/component
and looking up the efficiency function.
Efficiency functions are then executed against the X data of the input workspace
to generate new Y and E outputs for the denominatorWS.
@param wsIndex: The index of the histogram in the input workspace to process.
@param denominatorWS : Workspace that will become the denominator in the
normalisation routine.
@param inWS: Workspace input. Contains instrument to use as well as X data to
use.
@param prog: progress reporting object.
*/
void NormaliseByDetector::processHistogram(size_t wsIndex, const MatrixWorkspace_const_sptr &inWS,
const MatrixWorkspace_sptr &denominatorWS, Progress &prog) {
const auto ¶mMap = inWS->constInstrumentParameters();
const auto &spectrumInfo = inWS->spectrumInfo();
const auto &det = spectrumInfo.detector(wsIndex);
const std::string type = "fitting";
Geometry::Parameter_sptr foundParam = paramMap.getRecursiveByType(&det, type);
const Geometry::FitParameter &foundFittingParam = tryParseFunctionParameter(foundParam, det);
const std::string &fitFunctionName = foundFittingParam.getFunction();
IFunction_sptr function = FunctionFactory::Instance().createFunction(fitFunctionName);
using ParamNames = std::vector<std::string>;
ParamNames allParamNames = function->getParameterNames();
// Lookup each parameter name.
for (auto &name : allParamNames) {
Geometry::Parameter_sptr param = paramMap.getRecursive(&det, name, type);
const Geometry::FitParameter &fitParam = tryParseFunctionParameter(param, det);
if (fitParam.getFormula().empty()) {
throw std::runtime_error("A Forumla has not been provided for a fit function");
} else {
const std::string &resultUnitStr = fitParam.getResultUnit();
if (!resultUnitStr.empty() && resultUnitStr != "Wavelength") {
throw std::runtime_error("Units for function parameters must be in Wavelength");
}
}
mu::Parser p;
p.SetExpr(fitParam.getFormula());
double paramValue = p.Eval();
// Set the function coeffiecents.
function->setParameter(fitParam.getName(), paramValue);
}
auto wavelengths = inWS->points(wsIndex);
FunctionDomain1DVector domain(wavelengths.rawData());
FunctionValues values(domain);
function->function(domain, values);
auto &Y = denominatorWS->mutableY(wsIndex);
for (size_t i = 0; i < domain.size(); ++i) {
Y[i] = values[i];
}
denominatorWS->mutableE(wsIndex) = 0.0;
prog.report();
}
/**
Controlling function. Processes the histograms either in parallel or
sequentially.
@param inWS: Workspace input. Contains instrument to use as well as X data to
use.
*/
MatrixWorkspace_sptr NormaliseByDetector::processHistograms(const MatrixWorkspace_sptr &inWS) {
const size_t nHistograms = inWS->getNumberHistograms();
const auto progress_items = static_cast<size_t>(double(nHistograms) * 1.2);
Progress prog(this, 0.0, 1.0, progress_items);
// Clone the input workspace to create a template for the denominator
// workspace.
auto cloneAlg = createChildAlgorithm("CloneWorkspace", 0.0, 0.1, true);
cloneAlg->setProperty("InputWorkspace", inWS);
cloneAlg->setPropertyValue("OutputWorkspace", "temp");
cloneAlg->executeAsChildAlg();
Workspace_sptr temp = cloneAlg->getProperty("OutputWorkspace");
MatrixWorkspace_sptr denominatorWS = std::dynamic_pointer_cast<MatrixWorkspace>(temp);
// Choose between parallel execution and sequential execution then, process
// histograms accordingly.
if (m_parallelExecution) {
PARALLEL_FOR_IF(Kernel::threadSafe(*inWS, *denominatorWS))
for (int wsIndex = 0; wsIndex < static_cast<int>(nHistograms); ++wsIndex) {
PARALLEL_START_INTERRUPT_REGION
this->processHistogram(wsIndex, inWS, denominatorWS, prog);
PARALLEL_END_INTERRUPT_REGION
}
PARALLEL_CHECK_INTERRUPT_REGION
} else {
for (size_t wsIndex = 0; wsIndex < nHistograms; ++wsIndex) {
this->processHistogram(wsIndex, inWS, denominatorWS, prog);
}
}
return denominatorWS;
}
//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void NormaliseByDetector::exec() {
MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");
// Do the work of extracting functions and applying them to each bin on each
// histogram. The denominator workspace is mutable.
MatrixWorkspace_sptr denominatorWS = processHistograms(inWS);
// Perform the normalisation.
auto divideAlg = createChildAlgorithm("Divide", 0.9, 1.0, true);
divideAlg->setRethrows(true);
divideAlg->setProperty("LHSWorkspace", inWS);
divideAlg->setProperty("RHSWorkspace", denominatorWS);
divideAlg->executeAsChildAlg();
MatrixWorkspace_sptr outputWS = divideAlg->getProperty("OutputWorkspace");
setProperty("OutputWorkspace", outputWS);
}
} // namespace Mantid::Algorithms