-
Notifications
You must be signed in to change notification settings - Fork 122
/
MuonLoad.cpp
422 lines (334 loc) · 14.9 KB
/
MuonLoad.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
/*WIKI*
The algorithm replicates the sequence of actions undertaken by MuonAnalysis in order to produce a Muon workspace ready for fitting.
Specifically:
# Load the specified filename
# Apply dead time correction
# Group the workspace
# Offset, crop and rebin the workspace
# If the loaded data is multi-period - apply the specified operation to specified periods to get a single data set.
# Use [[MuonCalculateAsymmetry]] to get the resulting workspace.
*WIKI*/
#include "MantidWorkflowAlgorithms/MuonLoad.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/ListValidator.h"
#include "MantidAPI/FileProperty.h"
#include "MantidDataObjects/TableWorkspace.h"
namespace Mantid
{
namespace WorkflowAlgorithms
{
using namespace Kernel;
using namespace API;
using namespace DataObjects;
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(MuonLoad)
//----------------------------------------------------------------------------------------------
/**
* Constructor
*/
MuonLoad::MuonLoad()
{
}
//----------------------------------------------------------------------------------------------
/**
* Destructor
*/
MuonLoad::~MuonLoad()
{
}
//----------------------------------------------------------------------------------------------
/// Algorithm's name for identification. @see Algorithm::name
const std::string MuonLoad::name() const { return "MuonLoad";};
/// Algorithm's version for identification. @see Algorithm::version
int MuonLoad::version() const { return 1;};
/// Algorithm's category for identification. @see Algorithm::category
const std::string MuonLoad::category() const { return "Workflow\\Muon";}
//----------------------------------------------------------------------------------------------
/// Sets documentation strings for this algorithm
void MuonLoad::initDocs()
{
this->setWikiSummary("Loads Muon workspace ready for analysis.");
this->setOptionalMessage("Loads Muon workspace ready for analysis.");
}
//----------------------------------------------------------------------------------------------
/*
* Initialize the algorithm's properties.
*/
void MuonLoad::init()
{
declareProperty(new FileProperty("Filename", "", FileProperty::Load, ".nxs"),
"The name of the Nexus file to load" );
declareProperty("FirstPeriod", 0, "Group index of the first period workspace to use");
declareProperty("SecondPeriod", EMPTY_INT(), "Group index of the first period workspace to use");
std::vector<std::string> allowedOperations;
allowedOperations.push_back("+");
allowedOperations.push_back("-");
declareProperty("PeriodOperation","+", boost::make_shared<StringListValidator>(allowedOperations),
"If two periods specified, what operation to apply to workspaces to get a final one.");
declareProperty("ApplyDeadTimeCorrection", false,
"Whether dead time correction should be applied to loaded workspace");
declareProperty(new WorkspaceProperty<TableWorkspace>("CustomDeadTimeTable", "",Direction::Input,
PropertyMode::Optional),
"Table with dead time information. See LoadMuonNexus for format expected."
"If not specified -- algorithm tries to use dead times stored in the data file.");
declareProperty(new WorkspaceProperty<TableWorkspace>("DetectorGroupingTable","",Direction::Input,
PropertyMode::Optional),
"Table with detector grouping information. See LoadMuonNexus for format expected. "
"If not specified -- algorithm tries to get grouping information from the data file.");
declareProperty("TimeZero", EMPTY_DBL(), "Value used for Time Zero correction.");
declareProperty(new 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;
allowedTypes.push_back("PairAsymmetry");
allowedTypes.push_back("GroupAsymmetry");
allowedTypes.push_back("GroupCounts");
declareProperty("OutputType", "PairAsymmetry", boost::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(new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace","",Direction::Output),
"An output workspace.");
}
//----------------------------------------------------------------------------------------------
/**
* Execute the algorithm.
*/
void MuonLoad::exec()
{
const std::string filename = getProperty("Filename");
// Whether Dead Time Correction should be applied
bool applyDtc = getProperty("ApplyDeadTimeCorrection");
// If DetectorGropingTable not specified - use auto-grouping
bool autoGroup = ! static_cast<TableWorkspace_sptr>( getProperty("DetectorGroupingTable") );
// Load the file
IAlgorithm_sptr load = createChildAlgorithm("LoadMuonNexus");
load->setProperty("Filename", filename);
if ( applyDtc ) // Load dead times as well, if needed
load->setProperty("DeadTimeTable", "__NotUsed");
if ( autoGroup ) // Load grouping as well, if needed
load->setProperty("DetectorGroupingTable", "__NotUsed");
load->execute();
Workspace_sptr loadedWS = load->getProperty("OutputWorkspace");
MatrixWorkspace_sptr firstPeriodWS, secondPeriodWS;
// Deal with single-period workspace
if ( auto ws = boost::dynamic_pointer_cast<MatrixWorkspace>(loadedWS) )
{
if ( static_cast<int>( getProperty("FirstPeriod") ) != 0 )
throw std::invalid_argument("Single period data but first period is not 0.");
if ( static_cast<int>( getProperty("SecondPeriod") ) != EMPTY_INT() )
throw std::invalid_argument("Single period data but second period specified");
firstPeriodWS = ws;
}
// Deal with multi-period workspace
else if ( auto group = boost::dynamic_pointer_cast<WorkspaceGroup>(loadedWS) )
{
firstPeriodWS = getFirstPeriodWS(group);
secondPeriodWS = getSecondPeriodWS(group);
}
// Unexpected workspace type
else
{
throw std::runtime_error("Loaded workspace is of invalid type");
}
// Deal with dead time correction (if required)
if ( applyDtc )
{
TableWorkspace_sptr deadTimes = getProperty("CustomDeadTimeTable");
if ( ! deadTimes )
{
// If no dead times specified - try to use ones from the file
Workspace_sptr loadedDeadTimes = load->getProperty("DeadTimeTable");
if ( auto table = boost::dynamic_pointer_cast<TableWorkspace>(loadedDeadTimes) )
{
deadTimes = table;
}
else if ( auto group = boost::dynamic_pointer_cast<WorkspaceGroup>(loadedDeadTimes) )
{
// XXX: using first table only for now. Can use the one for appropriate period if necessary.
deadTimes = boost::dynamic_pointer_cast<TableWorkspace>( group->getItem(0) );
}
if ( ! deadTimes )
throw std::runtime_error("File doesn't contain any dead times");
}
firstPeriodWS = applyDTC(firstPeriodWS, deadTimes);
if ( secondPeriodWS )
secondPeriodWS = applyDTC(secondPeriodWS, deadTimes);
}
TableWorkspace_sptr grouping;
if ( autoGroup )
{
Workspace_sptr loadedGrouping = load->getProperty("DetectorGroupingTable");
if ( auto table = boost::dynamic_pointer_cast<TableWorkspace>(loadedGrouping) )
{
grouping = table;
}
else if ( auto group = boost::dynamic_pointer_cast<WorkspaceGroup>(loadedGrouping) )
{
// XXX: using first table only for now. Can use the one for appropriate period if necessary.
grouping = boost::dynamic_pointer_cast<TableWorkspace>( group->getItem(0) );
}
if ( ! grouping )
throw std::runtime_error("File doesn't contain grouping information");
}
else
{
grouping = getProperty("DetectorGroupingTable");
}
// Deal with grouping
firstPeriodWS = groupWorkspace(firstPeriodWS, grouping);
if ( secondPeriodWS )
secondPeriodWS = groupWorkspace(secondPeriodWS, grouping);
// Correct bin values
double loadedTimeZero = load->getProperty("TimeZero");
firstPeriodWS = correctWorkspace(firstPeriodWS, loadedTimeZero);
if ( secondPeriodWS )
secondPeriodWS = correctWorkspace(secondPeriodWS, loadedTimeZero);
IAlgorithm_sptr calcAssym = createChildAlgorithm("MuonCalculateAsymmetry");
// Set first period workspace
calcAssym->setProperty("FirstPeriodWorkspace", firstPeriodWS);
// Set second period workspace, if have one
if ( secondPeriodWS )
calcAssym->setProperty("SecondPeriodWorkspace", secondPeriodWS);
// Copy similar properties over
calcAssym->setProperty("PeriodOperation",
static_cast<std::string>( getProperty("PeriodOperation") ) );
calcAssym->setProperty("OutputType",
static_cast<std::string>( getProperty("OutputType") ) );
calcAssym->setProperty("PairFirstIndex",
static_cast<int>( getProperty("PairFirstIndex") ) );
calcAssym->setProperty("PairSecondIndex",
static_cast<int>( getProperty("PairSecondIndex") ) );
calcAssym->setProperty("Alpha",
static_cast<double>( getProperty("Alpha") ) );
calcAssym->setProperty("GroupIndex",
static_cast<int>( getProperty("GroupIndex") ) );
calcAssym->execute();
MatrixWorkspace_sptr outWS = calcAssym->getProperty("OutputWorkspace");
setProperty("OutputWorkspace", outWS);
}
/**
* Returns a workspace for the first period as specified using FirstPeriod property.
* @param group :: Loaded group of workspaces to use
* @return Workspace for the period
*/
MatrixWorkspace_sptr MuonLoad::getFirstPeriodWS(WorkspaceGroup_sptr group)
{
int firstPeriod = getProperty("FirstPeriod");
MatrixWorkspace_sptr resultWS;
if ( firstPeriod < 0 || firstPeriod >= static_cast<int>( group->size() ) )
throw std::invalid_argument("Workspace doesn't contain specified first period");
resultWS = boost::dynamic_pointer_cast<MatrixWorkspace>( group->getItem(firstPeriod) );
if ( ! resultWS )
throw std::invalid_argument("First period workspace is not a MatrixWorkspace");
return resultWS;
}
/**
* Returns a workspace for the second period as specified using SecondPeriod property.
* @param group :: Loaded group of workspaces to use
* @return Workspace for the period
*/
MatrixWorkspace_sptr MuonLoad::getSecondPeriodWS(WorkspaceGroup_sptr group)
{
int secondPeriod = getProperty("SecondPeriod");
MatrixWorkspace_sptr resultWS;
if ( secondPeriod != EMPTY_INT() )
{
if ( secondPeriod < 0 || secondPeriod >= static_cast<int>( group->size() ) )
throw std::invalid_argument("Workspace doesn't contain specified second period");
resultWS = boost::dynamic_pointer_cast<MatrixWorkspace>( group->getItem(secondPeriod) );
if ( ! resultWS )
throw std::invalid_argument("Second period workspace is not a MatrixWorkspace");
}
return resultWS;
}
/**
* Groups specified workspace according to specified DetectorGroupingTable.
* @param ws :: Workspace to group
* @param grouping :: Detector grouping table to use
* @return Grouped workspace
*/
MatrixWorkspace_sptr MuonLoad::groupWorkspace(MatrixWorkspace_sptr ws, TableWorkspace_sptr grouping)
{
IAlgorithm_sptr group = createChildAlgorithm("MuonGroupDetectors");
group->setProperty("InputWorkspace", ws);
group->setProperty("DetectorGroupingTable", grouping);
group->execute();
return group->getProperty("OutputWorkspace");
}
/**
* Applies dead time correction to the workspace.
* @param ws :: Workspace to apply correction
* @param dt :: Dead time table to use
* @return Corrected workspace
*/
MatrixWorkspace_sptr MuonLoad::applyDTC(MatrixWorkspace_sptr ws, TableWorkspace_sptr dt)
{
IAlgorithm_sptr dtc = createChildAlgorithm("ApplyDeadTimeCorr");
dtc->setProperty("InputWorkspace", ws);
dtc->setProperty("DeadTimeTable", dt);
dtc->execute();
return dtc->getProperty("OutputWorkspace");
}
/**
* 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 MuonLoad::correctWorkspace(MatrixWorkspace_sptr ws, double loadedTimeZero)
{
// Offset workspace, if need to
double timeZero = getProperty("TimeZero");
if ( timeZero != EMPTY_DBL() )
{
double offset = loadedTimeZero - timeZero;
IAlgorithm_sptr 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() )
{
IAlgorithm_sptr crop = createChildAlgorithm("CropWorkspace");
crop->setProperty("InputWorkspace", ws);
if ( Xmin != EMPTY_DBL() )
crop->setProperty("Xmin", Xmin);
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() )
{
IAlgorithm_sptr rebin = createChildAlgorithm("Rebin");
rebin->setProperty("InputWorkspace", ws);
rebin->setProperty("Params", rebinParams);
rebin->execute();
ws = rebin->getProperty("OutputWorkspace");
// TODO: following should be removed when FullBinsOnly function is added to Rebin algorithm
// Muon group don't want last bin if shorter than previous bins
double binSize = ws->dataX(0)[1] - ws->dataX(0)[0];
size_t numBins = ws->dataX(0).size();
if ( (ws->dataX(0)[numBins-1] - ws->dataX(0)[numBins-2]) != binSize )
{
IAlgorithm_sptr crop2 = createChildAlgorithm("CropWorkspace");
crop2->setProperty("InputWorkspace", ws);
crop2->setProperty("Xmax", ws->dataX(0)[numBins-2]);
crop2->execute();
ws = crop2->getProperty("OutputWorkspace");
}
}
return ws;
}
} // namespace WorkflowAlgorithms
} // namespace Mantid