-
Notifications
You must be signed in to change notification settings - Fork 122
/
GetDetectorOffsets.cpp
313 lines (267 loc) · 12.1 KB
/
GetDetectorOffsets.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
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
// 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/GetDetectorOffsets.h"
#include "MantidAPI/CompositeFunction.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/IPeakFunction.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/TextAxis.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidAlgorithms/PeakParameterHelper.h"
#include "MantidDataObjects/MaskWorkspace.h"
#include "MantidDataObjects/OffsetsWorkspace.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/ListValidator.h"
namespace Mantid::Algorithms {
// Register the class into the algorithm factory
DECLARE_ALGORITHM(GetDetectorOffsets)
using namespace Kernel;
using namespace Algorithms::PeakParameterHelper;
using namespace API;
using std::size_t;
using namespace DataObjects;
//-----------------------------------------------------------------------------------------
/** Initialisation method. Declares properties to be used in algorithm.
*/
void GetDetectorOffsets::init() {
declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
"A 2D workspace with X values of d-spacing");
auto mustBePositive = std::make_shared<BoundedValidator<double>>();
mustBePositive->setLower(0);
declareProperty("Step", 0.001, mustBePositive, "Step size used to bin d-spacing data");
declareProperty("DReference", 2.0, mustBePositive, "Center of reference peak in d-space");
declareProperty("XMin", 0.0, "Minimum of CrossCorrelation data to search for peak, usually negative");
declareProperty("XMax", 0.0, "Maximum of CrossCorrelation data to search for peak, usually positive");
declareProperty(std::make_unique<FileProperty>("GroupingFileName", "", FileProperty::OptionalSave, ".cal"),
"Optional: The name of the output CalFile to save the "
"generated OffsetsWorkspace.");
declareProperty(std::make_unique<WorkspaceProperty<OffsetsWorkspace>>("OutputWorkspace", "", Direction::Output),
"An output workspace containing the offsets.");
declareProperty(std::make_unique<WorkspaceProperty<>>("MaskWorkspace", "Mask", Direction::Output),
"An output workspace containing the mask.");
// Only keep peaks
declareProperty("PeakFunction", "Gaussian",
std::make_shared<StringListValidator>(FunctionFactory::Instance().getFunctionNames<IPeakFunction>()),
"The function type for fitting the peaks.");
declareProperty(std::make_unique<Kernel::PropertyWithValue<bool>>("EstimateFWHM", false),
"Whether to esimate FWHM of peak function when estimating fit parameters");
declareProperty("MaxOffset", 1.0, "Maximum absolute value of offsets; default is 1");
/* Signed mode calculates offset in number of bins */
std::vector<std::string> modes{"Relative", "Absolute", "Signed"};
declareProperty("OffsetMode", "Relative", std::make_shared<StringListValidator>(modes),
"Whether to calculate a relative, absolute, or signed offset");
declareProperty("DIdeal", 2.0, mustBePositive,
"The known peak centre value from the NIST standard "
"information, this is only used in Absolute OffsetMode.");
}
std::map<std::string, std::string> GetDetectorOffsets::validateInputs() {
std::map<std::string, std::string> result;
MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
if (!inputWS) {
result["InputWorkspace"] = "The InputWorkspace must be a MatrixWorkspace.";
return result;
}
const auto unit = inputWS->getAxis(0)->unit()->caption();
const auto unitErrorMsg =
"GetDetectorOffsets only supports input workspaces with units 'Bins of Shift' or 'd-Spacing', your unit was : " +
unit;
if (unit != "Bins of Shift" && unit != "d-Spacing") {
result["InputWorkspace"] = unitErrorMsg;
}
return result;
}
//-----------------------------------------------------------------------------------------
/** Executes the algorithm
*
* @throw Exception::FileError If the grouping file cannot be opened or read
*successfully
*/
void GetDetectorOffsets::exec() {
inputW = getProperty("InputWorkspace");
m_Xmin = getProperty("XMin");
m_Xmax = getProperty("XMax");
m_maxOffset = getProperty("MaxOffset");
if (m_Xmin >= m_Xmax)
throw std::runtime_error("Must specify m_Xmin<m_Xmax");
m_dreference = getProperty("DReference");
m_step = getProperty("Step");
m_estimateFWHM = getProperty("EstimateFWHM");
std::string mode_str = getProperty("OffsetMode");
if (mode_str == "Absolute") {
mode = offset_mode::absolute_offset;
}
else if (mode_str == "Relative") {
mode = offset_mode::relative_offset;
}
else if (mode_str == "Signed") {
mode = offset_mode::signed_offset;
}
m_dideal = getProperty("DIdeal");
int64_t nspec = inputW->getNumberHistograms();
// Create the output OffsetsWorkspace
auto outputW = std::make_shared<OffsetsWorkspace>(inputW->getInstrument());
// Create the output MaskWorkspace
auto maskWS = std::make_shared<MaskWorkspace>(inputW->getInstrument());
// To get the workspace index from the detector ID
const detid2index_map pixel_to_wi = maskWS->getDetectorIDToWorkspaceIndexMap(true);
// Fit all the spectra with a gaussian
Progress prog(this, 0.0, 1.0, nspec);
auto &spectrumInfo = maskWS->mutableSpectrumInfo();
PARALLEL_FOR_IF(Kernel::threadSafe(*inputW))
for (int64_t wi = 0; wi < nspec; ++wi) {
PARALLEL_START_INTERRUPT_REGION
// Fit the peak
double offset = fitSpectra(wi);
double mask = 0.0;
if (std::abs(offset) > m_maxOffset) {
offset = 0.0;
mask = 1.0;
}
// Get the list of detectors in this pixel
const auto &dets = inputW->getSpectrum(wi).getDetectorIDs();
// Most of the exec time is in FitSpectra, so this critical block should not
// be a problem.
PARALLEL_CRITICAL(GetDetectorOffsets_setValue) {
// Use the same offset for all detectors from this pixel
for (const auto &det : dets) {
outputW->setValue(det, offset);
const auto mapEntry = pixel_to_wi.find(det);
if (mapEntry == pixel_to_wi.end())
continue;
const size_t workspaceIndex = mapEntry->second;
if (mask == 1.) {
// Being masked
maskWS->getSpectrum(workspaceIndex).clearData();
spectrumInfo.setMasked(workspaceIndex, true);
maskWS->mutableY(workspaceIndex)[0] = mask;
} else {
// Using the detector
maskWS->mutableY(workspaceIndex)[0] = mask;
}
}
}
prog.report();
PARALLEL_END_INTERRUPT_REGION
}
PARALLEL_CHECK_INTERRUPT_REGION
// Return the output
setProperty("OutputWorkspace", outputW);
setProperty("MaskWorkspace", maskWS);
// Also save to .cal file, if requested
std::string filename = getProperty("GroupingFileName");
if (!filename.empty()) {
progress(0.9, "Saving .cal file");
auto childAlg = createChildAlgorithm("SaveCalFile");
childAlg->setProperty("OffsetsWorkspace", outputW);
childAlg->setProperty("MaskWorkspace", maskWS);
childAlg->setPropertyValue("Filename", filename);
childAlg->executeAsChildAlg();
}
}
//-----------------------------------------------------------------------------------------
/** Calls Gaussian1D as a child algorithm to fit the offset peak in a spectrum
*
* @param s :: The Workspace Index to fit
* @return The calculated offset value
*/
double GetDetectorOffsets::fitSpectra(const int64_t s) {
// Find point of peak centre
const auto &yValues = inputW->y(s);
auto it = std::max_element(yValues.cbegin(), yValues.cend());
// Set the default peak height and location
double peakHeight = *it;
const double peakLoc = inputW->x(s)[it - yValues.begin()];
// Return if peak of Cross Correlation is nan (Happens when spectra is zero)
// Pixel with large offset will be masked
if (std::isnan(peakHeight))
return (1000.);
IFunction_sptr fun_ptr = createFunction(peakHeight, peakLoc);
// Try to observe the peak height and location
const auto &histogram = inputW->histogram(s);
const auto &vector_x = histogram.points();
const auto start_index = findXIndex(vector_x, m_Xmin);
const auto stop_index = findXIndex(vector_x, m_Xmax, start_index);
// observe parameters if we found a peak range, otherwise use defaults
if (start_index != stop_index) {
// create a background function
auto bkgdFunction = std::dynamic_pointer_cast<IBackgroundFunction>(fun_ptr->getFunction(0));
auto peakFunction = std::dynamic_pointer_cast<IPeakFunction>(fun_ptr->getFunction(1));
int result = estimatePeakParameters(histogram, std::pair<size_t, size_t>(start_index, stop_index), peakFunction,
bkgdFunction, m_estimateFWHM, EstimatePeakWidth::Observation, EMPTY_DBL(), 0.0);
if (result != PeakFitResult::GOOD) {
g_log.debug() << "ws index: " << s
<< " bad result for observing peak parameters, using default peak height and loc\n";
}
} else {
g_log.notice() << "ws index: " << s
<< " range size is zero in estimatePeakParameters, using default peak height and loc\n";
}
IAlgorithm_sptr fit_alg;
try {
// set the ChildAlgorithm no to log as this will be run once per spectra
fit_alg = createChildAlgorithm("Fit", -1, -1, false);
} catch (Exception::NotFoundError &) {
g_log.error("Can't locate Fit algorithm");
throw;
}
fit_alg->setProperty("Function", fun_ptr);
fit_alg->setProperty("InputWorkspace", inputW);
fit_alg->setProperty<int>("WorkspaceIndex",
static_cast<int>(s)); // TODO what is the right thing to do here?
fit_alg->setProperty("StartX", m_Xmin);
fit_alg->setProperty("EndX", m_Xmax);
fit_alg->setProperty("MaxIterations", 100);
fit_alg->executeAsChildAlg();
std::string fitStatus = fit_alg->getProperty("OutputStatus");
// Pixel with large offset will be masked
if (fitStatus != "success")
return (1000.);
// std::vector<double> params = fit_alg->getProperty("Parameters");
API::IFunction_sptr function = fit_alg->getProperty("Function");
double offset = function->getParameter(3); // params[3]; // f1.PeakCentre
if (mode == offset_mode::signed_offset) {
// factor := factor * (1+offset) for d-spacemap conversion so factor cannot be
// negative
offset *= -1;
}
/* offset relative to the reference */
else if (mode == offset_mode::relative_offset) {
// factor := factor * (1+offset) for d-spacemap conversion so factor cannot be
// negative
offset = -1. * offset * m_step / (m_dreference + offset * m_step);
}
/* Offset relative to the ideal */
else if (mode == offset_mode::absolute_offset) {
offset = -1. * offset * m_step / (m_dreference + offset * m_step);
// translated from(DIdeal - FittedPeakCentre)/(FittedPeakCentre)
// given by Matt Tucker in ticket #10642
offset += (m_dideal - m_dreference) / m_dreference;
}
return offset;
}
/**
* Create a function string from the given parameters and the algorithm inputs
* @param peakHeight :: The height of the peak
* @param peakLoc :: The location of the peak
*/
IFunction_sptr GetDetectorOffsets::createFunction(const double peakHeight, const double peakLoc) {
const FunctionFactoryImpl &creator = FunctionFactory::Instance();
auto background = creator.createFunction("LinearBackground");
auto peak = std::dynamic_pointer_cast<IPeakFunction>(creator.createFunction(getProperty("PeakFunction")));
peak->setHeight(peakHeight);
peak->setCentre(peakLoc);
const double sigma(10.0);
peak->setFwhm(2.0 * std::sqrt(2.0 * M_LN2) * sigma);
auto fitFunc = new CompositeFunction(); // Takes ownership of the functions
fitFunc->addFunction(background);
fitFunc->addFunction(peak);
return std::shared_ptr<IFunction>(fitFunc);
}
} // namespace Mantid::Algorithms