-
Notifications
You must be signed in to change notification settings - Fork 122
/
UnwrapSNS.cpp
333 lines (286 loc) · 11.9 KB
/
UnwrapSNS.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
#include "MantidAlgorithms/UnwrapSNS.h"
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/RawCountValidator.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidDataObjects/EventList.h"
#include "MantidGeometry/Instrument.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/UnitFactory.h"
#include <limits>
namespace Mantid {
namespace Algorithms {
DECLARE_ALGORITHM(UnwrapSNS)
using namespace Kernel;
using namespace API;
using DataObjects::EventList;
using DataObjects::EventWorkspace;
using Kernel::Exception::InstrumentDefinitionError;
using Kernel::Exception::NotFoundError;
using std::size_t;
/// Default constructor
UnwrapSNS::UnwrapSNS()
: m_conversionConstant(0.), m_inputWS(), m_inputEvWS(), m_LRef(0.),
m_Tmin(0.), m_Tmax(0.), m_frameWidth(0.), m_numberOfSpectra(0),
m_XSize(0), m_progress(nullptr) {}
/// Destructor
UnwrapSNS::~UnwrapSNS() {
if (m_progress)
delete m_progress;
m_progress = nullptr;
}
/// Algorithm's name for identification overriding a virtual method
const std::string UnwrapSNS::name() const { return "UnwrapSNS"; }
/// Algorithm's version for identification overriding a virtual method
int UnwrapSNS::version() const { return 1; }
/// Algorithm's category for identification overriding a virtual method
const std::string UnwrapSNS::category() const {
return "CorrectionFunctions\\InstrumentCorrections";
}
/// Initialisation method
void UnwrapSNS::init() {
auto wsValidator = boost::make_shared<CompositeValidator>();
wsValidator->add<WorkspaceUnitValidator>("TOF");
wsValidator->add<HistogramValidator>();
wsValidator->add<RawCountValidator>();
wsValidator->add<InstrumentValidator>();
declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
"InputWorkspace", "", Direction::Input, wsValidator),
"Contains numbers counts against time of flight (TOF).");
declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
"OutputWorkspace", "", Direction::Output),
"This workspace will be in the units of time of flight. (See "
"http://www.mantidproject.org/Units)");
auto validator = boost::make_shared<BoundedValidator<double>>();
validator->setLower(0.01);
declareProperty("LRef", 0.0, validator,
"A distance at which it is possible to deduce if a particle "
"is from the current or a past frame based on its arrival "
"time. This time criterion can be set with the property "
"below e.g. correct when arrival time < Tmin.");
validator->setLower(0.01);
declareProperty("Tmin", Mantid::EMPTY_DBL(), validator,
"With LRef this defines the maximum speed expected for "
"particles. For each count or time bin the mean particle "
"speed is calculated and if this is greater than LRef/Tmin "
"its TOF is corrected.");
validator->setLower(0.01);
declareProperty("Tmax", Mantid::EMPTY_DBL(), validator,
"The maximum time of flight of the data used for the width "
"of the frame. If not set the maximum time of flight of the "
"data is used.");
// Calculate and set the constant factor for the conversion to wavelength
const double TOFisinMicroseconds = 1e6;
const double toAngstroms = 1e10;
m_conversionConstant = (PhysicalConstants::h * toAngstroms) /
(PhysicalConstants::NeutronMass * TOFisinMicroseconds);
}
/** Executes the algorithm
* @throw std::runtime_error if the workspace is invalid or a child algorithm
*fails
* @throw Kernel::Exception::InstrumentDefinitionError if detector, source or
*sample positions cannot be calculated
*
*/
void UnwrapSNS::exec() {
// Get the input workspace
m_inputWS = getProperty("InputWorkspace");
// Get the "reference" flightpath (currently passed in as a property)
m_LRef = getProperty("LRef");
m_XSize = static_cast<int>(m_inputWS->x(0).size());
m_numberOfSpectra = static_cast<int>(m_inputWS->getNumberHistograms());
g_log.debug() << "Number of spectra in input workspace: " << m_numberOfSpectra
<< "\n";
// go off and do the event version if appropriate
m_inputEvWS = boost::dynamic_pointer_cast<const EventWorkspace>(m_inputWS);
if ((m_inputEvWS != nullptr)) // && ! this->getProperty("ForceHist")) // TODO
// remove ForceHist option
{
this->execEvent();
return;
}
this->getTofRangeData(false);
// set up the progress bar
m_progress = new Progress(this, 0.0, 1.0, m_numberOfSpectra);
MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
if (outputWS != m_inputWS) {
outputWS = WorkspaceFactory::Instance().create(m_inputWS, m_numberOfSpectra,
m_XSize, m_XSize - 1);
setProperty("OutputWorkspace", outputWS);
}
// without the primary flight path the algorithm cannot work
const auto &spectrumInfo = m_inputWS->spectrumInfo();
const double L1 = spectrumInfo.l1();
PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWS, *outputWS))
for (int workspaceIndex = 0; workspaceIndex < m_numberOfSpectra;
workspaceIndex++) {
PARALLEL_START_INTERUPT_REGION
if (!spectrumInfo.hasDetectors(workspaceIndex)) {
// If the detector flightpath is missing, zero the data
g_log.debug() << "Detector information for workspace index "
<< workspaceIndex << " is not available.\n";
outputWS->setSharedX(workspaceIndex, m_inputWS->sharedX(workspaceIndex));
outputWS->mutableY(workspaceIndex) = 0.0;
outputWS->mutableE(workspaceIndex) = 0.0;
} else {
const double Ld = L1 + spectrumInfo.l2(workspaceIndex);
// fix the x-axis
std::vector<double> timeBins;
size_t pivot = this->unwrapX(m_inputWS->x(workspaceIndex), timeBins, Ld);
outputWS->setBinEdges(workspaceIndex, std::move(timeBins));
pivot++; // one-off difference between x and y
// fix the counts using the pivot point
auto &yIn = m_inputWS->y(workspaceIndex);
auto &yOut = outputWS->mutableY(workspaceIndex);
auto lengthFirstPartY = std::distance(yIn.begin() + pivot, yIn.end());
std::copy(yIn.begin() + pivot, yIn.end(), yOut.begin());
std::copy(yIn.begin(), yIn.begin() + pivot,
yOut.begin() + lengthFirstPartY);
// fix the uncertainties using the pivot point
auto &eIn = m_inputWS->e(workspaceIndex);
auto &eOut = outputWS->mutableE(workspaceIndex);
auto lengthFirstPartE = std::distance(eIn.begin() + pivot, eIn.end());
std::copy(eIn.begin() + pivot, eIn.end(), eOut.begin());
std::copy(eIn.begin(), eIn.begin() + pivot,
eOut.begin() + lengthFirstPartE);
}
m_progress->report();
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
m_inputWS.reset();
this->runMaskDetectors();
}
void UnwrapSNS::execEvent() {
// set up the output workspace
MatrixWorkspace_sptr matrixOutW = this->getProperty("OutputWorkspace");
if (matrixOutW != m_inputWS) {
matrixOutW = m_inputWS->clone();
setProperty("OutputWorkspace", matrixOutW);
}
auto outW = boost::dynamic_pointer_cast<EventWorkspace>(matrixOutW);
// set up the progress bar
m_progress = new Progress(this, 0.0, 1.0, m_numberOfSpectra * 2);
// algorithm assumes the data is sorted so it can jump out early
outW->sortAll(Mantid::DataObjects::TOF_SORT, m_progress);
this->getTofRangeData(true);
// without the primary flight path the algorithm cannot work
const auto &spectrumInfo = m_inputWS->spectrumInfo();
const double L1 = spectrumInfo.l1();
// do the actual work
for (int workspaceIndex = 0; workspaceIndex < m_numberOfSpectra;
workspaceIndex++) {
std::size_t numEvents = outW->getSpectrum(workspaceIndex).getNumberEvents();
double Ld = -1.0;
if (spectrumInfo.hasDetectors(workspaceIndex))
Ld = L1 + spectrumInfo.l2(workspaceIndex);
std::vector<double> time_bins;
if (outW->x(0).size() > 2) {
this->unwrapX(m_inputWS->x(workspaceIndex), time_bins, Ld);
outW->setBinEdges(workspaceIndex, std::move(time_bins));
} else {
outW->setSharedX(workspaceIndex, m_inputWS->sharedX(workspaceIndex));
}
if (numEvents > 0) {
std::vector<double> times(numEvents);
outW->getSpectrum(workspaceIndex).getTofs(times);
double filterVal = m_Tmin * Ld / m_LRef;
for (size_t j = 0; j < numEvents; j++) {
if (times[j] < filterVal)
times[j] += m_frameWidth;
else
break; // stop filtering
}
outW->getSpectrum(workspaceIndex).setTofs(times);
}
m_progress->report();
}
outW->clearMRU();
this->runMaskDetectors();
}
int UnwrapSNS::unwrapX(const Mantid::HistogramData::HistogramX &datain,
std::vector<double> &dataout, const double &Ld) {
std::vector<double> tempX_L; // lower half - to be frame wrapped
tempX_L.reserve(m_XSize);
tempX_L.clear();
std::vector<double> tempX_U; // upper half - to not be frame wrapped
tempX_U.reserve(m_XSize);
tempX_U.clear();
double filterVal = m_Tmin * Ld / m_LRef;
dataout.clear();
int specialBin = 0;
for (int bin = 0; bin < m_XSize; ++bin) {
// This is the time-of-flight value under consideration in the current
// iteration of the loop
const double tof = datain[bin];
if (tof < filterVal) {
tempX_L.push_back(tof + m_frameWidth);
// Record the bins that fall in this range for copying over the data &
// errors
if (specialBin < bin)
specialBin = bin;
} else {
tempX_U.push_back(tof);
}
} // loop over X values
// now put it back into the vector supplied
dataout.clear();
dataout.insert(dataout.begin(), tempX_U.begin(), tempX_U.end());
dataout.insert(dataout.end(), tempX_L.begin(), tempX_L.end());
assert(datain.size() == dataout.size());
return specialBin;
}
void UnwrapSNS::getTofRangeData(const bool isEvent) {
// get the Tmin/Tmax properties
m_Tmin = this->getProperty("Tmin");
m_Tmax = this->getProperty("Tmax");
// if either the values are not specified by properties, find them from the
// data
double empty = Mantid::EMPTY_DBL();
if ((m_Tmin == empty) || (m_Tmax == empty)) {
// get data min/max values
double dataTmin;
double dataTmax;
if (isEvent) {
m_inputEvWS->sortAll(DataObjects::TOF_SORT, nullptr);
m_inputEvWS->getEventXMinMax(dataTmin, dataTmax);
} else {
m_inputWS->getXMinMax(dataTmin, dataTmax);
}
// fix the unspecified values
if (m_Tmin == empty) {
m_Tmin = dataTmin;
}
if (m_Tmax == empty) {
m_Tmax = dataTmax;
}
}
// check the frame width
m_frameWidth = m_Tmax - m_Tmin;
g_log.information() << "Frame range in microseconds is: " << m_Tmin << " - "
<< m_Tmax << "\n";
if (m_Tmin < 0.)
throw std::runtime_error("Cannot have Tmin less than zero");
if (m_Tmin > m_Tmax)
throw std::runtime_error("Have case of Tmin > Tmax");
g_log.information() << "Wavelength cuttoff is : "
<< (m_conversionConstant * m_Tmin / m_LRef)
<< "Angstrom, Frame width is: " << m_frameWidth
<< "microseconds\n";
}
void UnwrapSNS::runMaskDetectors() {
IAlgorithm_sptr alg = createChildAlgorithm("MaskDetectors");
alg->setProperty<MatrixWorkspace_sptr>("Workspace",
this->getProperty("OutputWorkspace"));
alg->setProperty<MatrixWorkspace_sptr>("MaskedWorkspace",
this->getProperty("InputWorkspace"));
if (!alg->execute())
throw std::runtime_error(
"MaskDetectors Child Algorithm has not executed successfully");
}
} // namespace Algorithm
} // namespace Mantid