-
Notifications
You must be signed in to change notification settings - Fork 122
/
MuonProcess.cpp
428 lines (379 loc) · 17.1 KB
/
MuonProcess.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
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
// 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 "MantidWorkflowAlgorithms/MuonProcess.h"
#include "MantidAPI/FileProperty.h"
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/MandatoryValidator.h"
#include "MantidWorkflowAlgorithms/MuonGroupAsymmetryCalculator.h"
#include "MantidWorkflowAlgorithms/MuonGroupCountsCalculator.h"
#include "MantidWorkflowAlgorithms/MuonPairAsymmetryCalculator.h"
// free functions
namespace {
/**
* Tests if argument is equal to one.
* @param i :: [input] integer to test
* @returns True if i == 1, else false.
*/
bool isOne(int i) { return (i == 1); }
} // namespace
namespace Mantid {
namespace WorkflowAlgorithms {
using namespace Kernel;
using namespace API;
using namespace DataObjects;
using API::WorkspaceGroup_sptr;
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(MuonProcess)
//----------------------------------------------------------------------------------------------
/// Algorithm's name for identification. @see Algorithm::name
const std::string MuonProcess::name() const { return "MuonProcess"; }
/// Algorithm's version for identification. @see Algorithm::version
int MuonProcess::version() const { return 1; }
/// Algorithm's category for identification. @see Algorithm::category
const std::string MuonProcess::category() const { return "Workflow\\Muon"; }
//----------------------------------------------------------------------------------------------
//----------------------------------------------------------------------------------------------
/*
* Initialize the algorithm's properties.
*/
void MuonProcess::init() {
declareProperty(
std::make_unique<WorkspaceProperty<Workspace>>("InputWorkspace", "", Direction::Input, PropertyMode::Mandatory),
"Input workspace loaded from file (e.g. by LoadMuonNexus)");
std::vector<std::string> allowedModes{"CorrectAndGroup", "Analyse", "Combined"};
auto modeVal = std::make_shared<CompositeValidator>();
modeVal->add(std::make_shared<StringListValidator>(allowedModes));
modeVal->add(std::make_shared<MandatoryValidator<std::string>>());
declareProperty("Mode", "Combined", modeVal,
"Mode to run in. CorrectAndGroup applies dead time "
"correction and grouping; Analyse changes bin offset, "
"crops, rebins and calculates asymmetry; Combined does all "
"of the above.");
declareProperty(std::make_unique<ArrayProperty<int>>("SummedPeriodSet", Direction::Input),
"Comma-separated list of periods to be summed");
declareProperty(std::make_unique<ArrayProperty<int>>("SubtractedPeriodSet", Direction::Input),
"Comma-separated list of periods to be subtracted from the "
"SummedPeriodSet");
declareProperty("ApplyDeadTimeCorrection", false,
"Whether dead time correction should be applied to loaded workspace");
declareProperty(std::make_unique<WorkspaceProperty<TableWorkspace>>("DeadTimeTable", "", Direction::Input,
PropertyMode::Optional),
"Table with dead time information, e.g. from LoadMuonNexus."
"Must be specified if ApplyDeadTimeCorrection is set true.");
declareProperty(std::make_unique<WorkspaceProperty<TableWorkspace>>("DetectorGroupingTable", "", Direction::Input,
PropertyMode::Optional),
"Table with detector grouping information, e.g. from LoadMuonNexus.");
declareProperty("TimeZero", EMPTY_DBL(), "Value used for Time Zero correction");
declareProperty("LoadedTimeZero", EMPTY_DBL(), std::make_shared<MandatoryValidator<double>>(),
"Time Zero value loaded from file, e.g. from LoadMuonNexus.");
declareProperty(std::make_unique<ArrayProperty<double>>("RebinParams"),
"Params used for rebinning. If empty - rebinning is not done.");
declareProperty("Xmin", EMPTY_DBL(), "Minimal X value to include");
declareProperty("Xmax", EMPTY_DBL(), "Maximal X value to include");
std::vector<std::string> allowedTypes{"PairAsymmetry", "GroupAsymmetry", "GroupCounts"};
declareProperty("OutputType", "PairAsymmetry", std::make_shared<StringListValidator>(allowedTypes),
"What kind of workspace required for analysis.");
declareProperty("PairFirstIndex", EMPTY_INT(), "Workspace index of the first pair group");
declareProperty("PairSecondIndex", EMPTY_INT(), "Workspace index of the second pair group");
declareProperty("Alpha", 1.0, "Alpha value of the pair");
declareProperty("GroupIndex", EMPTY_INT(), "Workspace index of the group");
declareProperty(std::make_unique<WorkspaceProperty<Workspace>>("OutputWorkspace", "", Direction::Output),
"An output workspace.");
declareProperty("CropWorkspace", true,
"Determines if the input workspace "
"should be cropped at Xmax, Xmin is "
"still aplied.");
declareProperty("WorkspaceName", "", "The name of the input workspace");
}
//----------------------------------------------------------------------------------------------
/**
* Execute the algorithm.
*/
void MuonProcess::exec() {
Progress progress(this, 0.0, 1.0, 5);
// Supplied input workspace
Workspace_sptr inputWS = getProperty("InputWorkspace");
std::vector<int> summedPeriods = getProperty("SummedPeriodSet");
std::vector<int> subtractedPeriods = getProperty("SubtractedPeriodSet");
auto allPeriodsWS = std::make_shared<WorkspaceGroup>();
progress.report();
// Deal with single-period workspace
if (auto ws = std::dynamic_pointer_cast<MatrixWorkspace>(inputWS)) {
allPeriodsWS->addWorkspace(ws);
}
// Deal with multi-period workspace
else if (auto group = std::dynamic_pointer_cast<WorkspaceGroup>(inputWS)) {
allPeriodsWS = group;
}
progress.report();
// Check mode
const std::string mode = getProperty("Mode");
// Dead time correction and grouping
if (mode != "Analyse") {
bool applyDtc = getProperty("ApplyDeadTimeCorrection");
// Deal with dead time correction (if required)
if (applyDtc) {
TableWorkspace_sptr deadTimes = getProperty("DeadTimeTable");
allPeriodsWS = applyDTC(allPeriodsWS, deadTimes);
}
progress.report();
TableWorkspace_sptr grouping;
grouping = getProperty("DetectorGroupingTable");
progress.report();
// Deal with grouping
allPeriodsWS = groupWorkspaces(allPeriodsWS, grouping);
} else {
progress.report();
progress.report();
}
// If not analysing, the present WS will be the output
Workspace_sptr outWS = allPeriodsWS;
if (mode != "CorrectAndGroup") {
// Correct bin values
double loadedTimeZero = getProperty("LoadedTimeZero");
allPeriodsWS = correctWorkspaces(allPeriodsWS, loadedTimeZero);
// Perform appropriate calculation
std::string outputType = getProperty("OutputType");
int groupIndex = getProperty("GroupIndex");
std::unique_ptr<IMuonAsymmetryCalculator> asymCalc;
if (outputType == "GroupCounts") {
asymCalc =
std::make_unique<MuonGroupCountsCalculator>(allPeriodsWS, summedPeriods, subtractedPeriods, groupIndex);
} else if (outputType == "GroupAsymmetry") {
asymCalc = std::make_unique<MuonGroupAsymmetryCalculator>(allPeriodsWS, summedPeriods, subtractedPeriods,
groupIndex, getProperty("Xmin"), getProperty("Xmax"),
getProperty("WorkspaceName"));
} else if (outputType == "PairAsymmetry") {
int first = getProperty("PairFirstIndex");
int second = getProperty("PairSecondIndex");
double alpha = getProperty("Alpha");
asymCalc = std::make_unique<MuonPairAsymmetryCalculator>(allPeriodsWS, summedPeriods, subtractedPeriods, first,
second, alpha);
}
progress.report();
outWS = asymCalc->calculate();
}
setProperty("OutputWorkspace", outWS);
}
/**
* Groups specified workspace group according to specified
* DetectorGroupingTable.
* @param wsGroup :: WorkspaceGroup to group
* @param grouping :: Detector grouping table to use
* @return Grouped workspaces
*/
WorkspaceGroup_sptr MuonProcess::groupWorkspaces(const WorkspaceGroup_sptr &wsGroup,
const TableWorkspace_sptr &grouping) {
WorkspaceGroup_sptr outWS = std::make_shared<WorkspaceGroup>();
for (int i = 0; i < wsGroup->getNumberOfEntries(); i++) {
auto ws = std::dynamic_pointer_cast<MatrixWorkspace>(wsGroup->getItem(i));
if (ws) {
MatrixWorkspace_sptr result;
auto group = createChildAlgorithm("MuonGroupDetectors");
group->setProperty("InputWorkspace", ws);
group->setProperty("DetectorGroupingTable", grouping);
group->execute();
result = group->getProperty("OutputWorkspace");
outWS->addWorkspace(result);
}
}
return outWS;
}
/**
* Applies dead time correction to the workspace group.
* @param wsGroup :: Workspace group to apply correction to
* @param dt :: Dead time table to use
* @return Corrected workspace group
*/
WorkspaceGroup_sptr MuonProcess::applyDTC(const WorkspaceGroup_sptr &wsGroup, const TableWorkspace_sptr &dt) {
WorkspaceGroup_sptr outWS = std::make_shared<WorkspaceGroup>();
for (int i = 0; i < wsGroup->getNumberOfEntries(); i++) {
auto ws = std::dynamic_pointer_cast<MatrixWorkspace>(wsGroup->getItem(i));
if (ws) {
MatrixWorkspace_sptr result;
auto dtc = createChildAlgorithm("ApplyDeadTimeCorr");
dtc->setProperty("InputWorkspace", ws);
dtc->setProperty("DeadTimeTable", dt);
dtc->execute();
result = dtc->getProperty("OutputWorkspace");
if (result) {
outWS->addWorkspace(result);
} else {
throw std::runtime_error("ApplyDeadTimeCorr failed to apply dead time "
"correction in MuonProcess");
}
}
}
return outWS;
}
/**
* Applies offset, crops and rebin the workspaces in the group according to
* specified params.
* @param wsGroup :: Workspaces to correct
* @param loadedTimeZero :: Time zero of the data, so we can calculate the
* offset
* @return Corrected workspaces
*/
WorkspaceGroup_sptr MuonProcess::correctWorkspaces(const WorkspaceGroup_sptr &wsGroup, double loadedTimeZero) {
WorkspaceGroup_sptr outWS = std::make_shared<WorkspaceGroup>();
for (int i = 0; i < wsGroup->getNumberOfEntries(); i++) {
auto ws = std::dynamic_pointer_cast<MatrixWorkspace>(wsGroup->getItem(i));
if (ws) {
MatrixWorkspace_sptr result;
result = correctWorkspace(ws, loadedTimeZero);
outWS->addWorkspace(result);
}
}
return outWS;
}
/**
* Applies offset, crops and rebin the workspace according to specified params.
* @param ws :: Workspace to correct
* @param loadedTimeZero :: Time zero of the data, so we can calculate the
* offset
* @return Corrected workspace
*/
MatrixWorkspace_sptr MuonProcess::correctWorkspace(MatrixWorkspace_sptr ws, double loadedTimeZero) {
// Offset workspace, if need to
double timeZero = getProperty("TimeZero");
if (timeZero != EMPTY_DBL()) {
double offset = loadedTimeZero - timeZero;
auto changeOffset = createChildAlgorithm("ChangeBinOffset");
changeOffset->setProperty("InputWorkspace", ws);
changeOffset->setProperty("Offset", offset);
changeOffset->execute();
ws = changeOffset->getProperty("OutputWorkspace");
}
// Crop workspace, if need to
double Xmin = getProperty("Xmin");
double Xmax = getProperty("Xmax");
if (Xmin != EMPTY_DBL() || Xmax != EMPTY_DBL()) {
auto crop = createChildAlgorithm("CropWorkspace");
crop->setProperty("InputWorkspace", ws);
if (Xmin != EMPTY_DBL())
crop->setProperty("Xmin", Xmin);
bool toCrop = getProperty("CropWorkspace");
if (toCrop) {
if (Xmax != EMPTY_DBL())
crop->setProperty("Xmax", Xmax);
}
crop->execute();
ws = crop->getProperty("OutputWorkspace");
}
// Rebin workspace if need to
std::vector<double> rebinParams = getProperty("RebinParams");
if (!rebinParams.empty()) {
auto rebin = createChildAlgorithm("Rebin");
rebin->setProperty("InputWorkspace", ws);
rebin->setProperty("Params", rebinParams);
rebin->setProperty("FullBinsOnly", true);
rebin->execute();
ws = rebin->getProperty("OutputWorkspace");
}
return ws;
}
/**
* Performs validation of inputs to the algorithm.
* - Input workspace must be single-period or a workspace group
* - Single-period input must only use period 1
* - Supplied period numbers must all be valid (between 1 and total number of
* periods)
* - If analysis will take place, SummedPeriodSet is mandatory
* - If ApplyDeadTimeCorrection is true, DeadTimeTable is mandatory
* @returns Map of parameter names to errors
*/
std::map<std::string, std::string> MuonProcess::validateInputs() {
std::map<std::string, std::string> errors;
// Supplied input workspace and sets of periods
const std::string propInputWS("InputWorkspace"), propSummedPeriodSet("SummedPeriodSet"),
propSubtractedPeriodSet("SubtractedPeriodSet");
Workspace_sptr inputWS = getProperty(propInputWS);
std::vector<int> summedPeriods = getProperty(propSummedPeriodSet);
std::vector<int> subtractedPeriods = getProperty(propSubtractedPeriodSet);
// If single-period data, test the sets of periods specified
if (auto ws = std::dynamic_pointer_cast<MatrixWorkspace>(inputWS)) {
if (std::find_if_not(summedPeriods.begin(), summedPeriods.end(), isOne) != summedPeriods.end()) {
errors[propSummedPeriodSet] = "Single period data but set of periods to "
"sum contains invalid values.";
}
if (!subtractedPeriods.empty()) {
errors[propSubtractedPeriodSet] = "Single period data but second set of periods specified";
}
} else {
// If not a MatrixWorkspace, must be a multi-period WorkspaceGroup
auto group = std::dynamic_pointer_cast<WorkspaceGroup>(inputWS);
if (group == nullptr) {
errors[propInputWS] = "Input workspace is of invalid type";
} else {
auto numPeriods = group->getNumberOfEntries();
if (numPeriods < 1) {
errors[propInputWS] = "Input workspace contains no periods";
}
// check summed period numbers
std::vector<int> invalidPeriods;
std::copy_if(summedPeriods.cbegin(), summedPeriods.cend(), std::back_inserter(invalidPeriods),
[numPeriods](auto period) { return period < 1 || period > numPeriods; });
if (!invalidPeriods.empty()) {
errors[propSummedPeriodSet] = buildErrorString(invalidPeriods);
invalidPeriods.clear();
}
// check subtracted period numbers
std::copy_if(subtractedPeriods.cbegin(), subtractedPeriods.cend(), std::back_inserter(invalidPeriods),
[numPeriods](auto period) { return period < 1 || period > numPeriods; });
if (!invalidPeriods.empty()) {
errors[propSubtractedPeriodSet] = buildErrorString(invalidPeriods);
}
}
}
// Some parameters can be mandatory or not depending on the mode
const std::string propMode("Mode"), propApplyDTC("ApplyDeadTimeCorrection"), propDeadTime("DeadTimeTable"),
propDetGroup("DetectorGroupingTable");
const std::string mode = getProperty(propMode);
// If analysis will take place, SummedPeriodSet is mandatory
if (mode != "CorrectAndGroup") {
if (summedPeriods.empty()) {
errors[propSummedPeriodSet] = "Cannot analyse: list of periods to sum was empty";
}
}
// If correcting/grouping will take place, DetectorGroupingTable is mandatory
if (mode != "Analyse") {
TableWorkspace_sptr grouping = getProperty(propDetGroup);
if (grouping == nullptr) {
errors[propDetGroup] = "No detector grouping table supplied";
}
}
// If dead time correction is to be applied, must supply dead times
bool applyDtc = getProperty(propApplyDTC);
if (applyDtc) {
TableWorkspace_sptr deadTimes = getProperty(propDeadTime);
if (!deadTimes) {
errors[propDeadTime] = "Cannot apply dead time correction as no dead times were supplied";
}
}
return errors;
}
/**
* Builds an error message from the supplied parameters.
* @param invalidPeriods :: [input] Vector containing invalid periods
* @returns An error message
*/
std::string MuonProcess::buildErrorString(const std::vector<int> &invalidPeriods) const {
std::stringstream message;
message << "Invalid periods specified: ";
for (auto it = invalidPeriods.begin(); it != invalidPeriods.end(); it++) {
message << *it;
if (it != invalidPeriods.end() - 1) {
message << ", ";
}
}
return message.str();
}
} // namespace WorkflowAlgorithms
} // namespace Mantid