Mantid
Loading...
Searching...
No Matches
LoadEMU.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2020 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
7
10#include "MantidAPI/Axis.h"
14#include "MantidAPI/Run.h"
26
27#include <boost/filesystem.hpp>
28#include <boost/math/special_functions/round.hpp>
29#include <boost/math/tools/minima.hpp>
30
31#include <Poco/AutoPtr.h>
32#include <Poco/TemporaryFile.h>
33#include <Poco/Util/PropertyFileConfiguration.h>
34
35#include <algorithm>
36#include <cmath>
37#include <cstdio>
38#include <fstream>
39#include <utility>
40
43namespace {
44
45// number of physical detectors
46constexpr size_t HORIZONTAL_TUBES = 16;
47constexpr size_t VERTICAL_TUBES = 35;
48constexpr size_t DETECTOR_TUBES = HORIZONTAL_TUBES + VERTICAL_TUBES;
49// analysed and direct detectors
50constexpr size_t HISTO_BINS_X = DETECTOR_TUBES * 2;
51constexpr size_t HISTO_BINS_Y = 1024;
52constexpr size_t HISTO_BINS_Y_DENUMERATOR = 16;
53constexpr size_t PIXELS_PER_TUBE = HISTO_BINS_Y / HISTO_BINS_Y_DENUMERATOR;
54
55constexpr size_t BM_HISTOGRAMS = HISTO_BINS_X * PIXELS_PER_TUBE;
56constexpr size_t HISTOGRAMS = BM_HISTOGRAMS + PIXELS_PER_TUBE;
57constexpr size_t BEAM_MONITOR_BINS = 100;
58constexpr size_t PSEUDO_BM_TUBE = 55;
59// the half window for the running average matches the plateau width for the peak
60constexpr size_t BM_HALF_WINDOW = 5;
61
62// File loading progress boundaries
63constexpr size_t Progress_LoadBinFile = 48;
64constexpr size_t Progress_ReserveMemory = 4;
67// Algorithm parameter names
68constexpr char FilenameStr[] = "Filename";
69constexpr char MaskStr[] = "Mask";
70constexpr char SelectDetectorTubesStr[] = "SelectDetectorTubes";
71constexpr char SelectDatasetStr[] = "SelectDataset";
72constexpr char OverrideDopplerFreqStr[] = "OverrideDopplerFrequency";
73constexpr char OverrideDopplerPhaseStr[] = "OverrideDopplerPhase";
74constexpr char FilterByTimeStartStr[] = "FilterByTimeStart";
75constexpr char FilterByTimeStopStr[] = "FilterByTimeStop";
76constexpr char RawDopplerTimeStr[] = "LoadAsRawDopplerTime";
77constexpr char IncludePseudoBMStr[] = "IncludeBeamMonitor";
78constexpr char CalibrateDopplerPhaseStr[] = "CalibrateDopplerPhase";
79constexpr char PathToBinaryStr[] = "BinaryEventPath";
80
81// Common pairing of limits
82using TimeLimits = std::pair<double, double>;
83
84template <typename Type>
85void AddSinglePointTimeSeriesProperty(API::LogManager &logManager, const std::string &time, const std::string &name,
86 const Type value) {
87 // create time series property and add single value
88 auto p = new Kernel::TimeSeriesProperty<Type>(name);
89 p->addValue(time, value);
90
91 // add to log manager
92 logManager.addProperty(p);
93}
94
95// Utility functions for loading values with defaults
96// Single value properties only support int, double, string and bool
97template <typename Type>
98Type GetNeXusValue(const NeXus::NXEntry &entry, const std::string &path, const Type &defval, int32_t index) {
99 try {
100 NeXus::NXDataSetTyped<Type> dataSet = entry.openNXDataSet<Type>(path);
101 dataSet.load();
102
103 return dataSet()[index];
104 } catch (std::runtime_error &) {
105 return defval;
106 }
107}
108// string and double are special cases
109template <>
110double GetNeXusValue<double>(const NeXus::NXEntry &entry, const std::string &path, const double &defval,
111 int32_t index) {
112 try {
113 NeXus::NXDataSetTyped<float> dataSet = entry.openNXDataSet<float>(path);
114 dataSet.load();
115
116 return dataSet()[index];
117 } catch (std::runtime_error &) {
118 return defval;
119 }
120}
121template <>
122std::string GetNeXusValue<std::string>(const NeXus::NXEntry &entry, const std::string &path, const std::string &defval,
123 int32_t /*unused*/) {
124
125 try {
126 NeXus::NXChar dataSet = entry.openNXChar(path);
127 dataSet.load();
128
129 return std::string(dataSet(), dataSet.dim0());
130 } catch (std::runtime_error &) {
131 return defval;
132 }
133}
134
135template <typename T>
136void MapNeXusToProperty(NeXus::NXEntry &entry, const std::string &path, const T &defval, API::LogManager &logManager,
137 const std::string &name, const T &factor, int32_t index) {
138
139 T value = GetNeXusValue<T>(entry, path, defval, index);
140 logManager.addProperty<T>(name, value * factor);
141}
142
143// sting is a special case
144template <>
145void MapNeXusToProperty<std::string>(NeXus::NXEntry &entry, const std::string &path, const std::string &defval,
146 API::LogManager &logManager, const std::string &name,
147 const std::string & /*unused*/, int32_t index) {
148
149 std::string value = GetNeXusValue<std::string>(entry, path, defval, index);
150 logManager.addProperty<std::string>(name, value);
151}
152
153template <typename T>
154void MapNeXusToSeries(NeXus::NXEntry &entry, const std::string &path, const T &defval, API::LogManager &logManager,
155 const std::string &time, const std::string &name, const T &factor, int32_t index) {
156
157 auto value = GetNeXusValue<T>(entry, path, defval, index);
158 AddSinglePointTimeSeriesProperty<T>(logManager, time, name, value * factor);
159}
160
161// map the comma separated range of indexes to the vector via a lambda function
162// throws an exception if it is outside the vector range
163//
164template <typename T, typename F> void mapRangeToIndex(const std::string &line, std::vector<T> &result, const F &fn) {
165
166 std::stringstream ss(line);
167 std::string item;
168 size_t index = 0;
169 while (std::getline(ss, item, ',')) {
170 auto k = item.find('-');
171
172 size_t p0, p1;
173 if (k != std::string::npos) {
174 p0 = boost::lexical_cast<size_t>(item.substr(0, k));
175 p1 = boost::lexical_cast<size_t>(item.substr(k + 1, item.size() - k - 1));
176 } else {
177 p0 = boost::lexical_cast<size_t>(item);
178 p1 = p0;
179 }
180
181 if (p1 < result.size() && p0 <= p1) {
182 while (p0 <= p1) {
183 result[p0++] = fn(index);
184 index++;
185 }
186 } else if (p0 < result.size() && p1 < p0) {
187 do {
188 result[p0] = fn(index);
189 index++;
190 } while (p1 < p0--);
191 } else
192 throw std::invalid_argument("invalid range specification");
193 }
194}
195
196// secant invert fucntion
197//
198template <typename F> double invert(double y, const F &f, double x0 = 0.0, const double eps = 1e-16) {
199 // secant method
200 double e0 = f(x0) - y;
201
202 double x1 = x0 + eps;
203 double e1 = f(x1) - y;
204 int loop = 16;
205 while (fabs(e0) > eps && loop-- > 0) {
206 double x = (x1 * e0 - x0 * e1) / (e0 - e1);
207
208 x1 = x0;
209 e1 = e0;
210
211 x0 = x;
212 e0 = f(x0) - y;
213 }
214
215 return x0;
216}
217
218// Need to convert two different methods for direct and analysed data
219// provide two distinct methods that may be called but the distances between
220// the detectors is not constant so it needs to store the distance for each
221// Note that observation time and TOF are in usec
222//
223using TofData = std::tuple<double, double>;
224
225class ConvertTOF {
226 const double m_w;
227 const double m_phi;
228 const double m_L0;
229 const double m_v2;
230 const double m_A;
231 const std::vector<double> &m_L2;
232
233 inline double L1(double t) const { return m_L0 + m_A * sin(m_w * t + m_phi); }
234
235 inline double v1(double t) const { return m_v2 - m_A * m_w * cos(m_w * t + m_phi); }
236
237public:
238 ConvertTOF(double Amp, double freq, double phase, double L1, double v2, const std::vector<double> &L2)
239 : m_w(2 * M_PI * freq), m_phi(M_PI * phase / 180.0), m_L0(L1), m_v2(v2), m_A(Amp), m_L2(L2) {}
240
241 TofData directTOF(size_t detID, double tobs) const {
242
243 // observation time and tof are in usec
244 auto tn = [=](double t) { return t + (L1(t) + m_L2[detID]) / v1(t); };
245
246 double tsec = tobs * 1.0e-6;
247 double t0 = tsec - (m_L0 + m_L2[detID]) / m_v2;
248 double tinv = invert(tsec, tn, t0);
249 double tof = (m_L0 + m_L2[detID]) / v1(tinv);
250
251 return TofData(tinv * 1.0e6, tof * 1.0e6);
252 }
253
254 TofData analysedTOF(size_t detID, double tobs) const {
255 // observation time and tof are in usec
256 auto tn = [=](double t) { return t + L1(t) / v1(t) + m_L2[detID] / m_v2; };
257
258 double tsec = tobs * 1.0e-6;
259 double t0 = tsec - (m_L0 + m_L2[detID]) / m_v2;
260 double t = invert(tsec, tn, t0);
261 double tof = m_L0 / v1(t) + m_L2[detID] / m_v2;
262
263 return TofData(t * 1.0e6, tof * 1.0e6);
264 }
265};
266
267// calculate mean of a subset of the vector
268double maskedMean(const std::vector<double> &vec, const std::vector<bool> &mask) {
269 if (vec.size() == 0 || vec.size() != mask.size())
270 throw std::runtime_error("masked mean of empty or mismatched vectors");
271 double sum = 0.0;
272 size_t count = 0;
273 for (size_t i = 0; i != vec.size(); i++) {
274 if (!mask[i])
275 continue;
276 sum += vec[i];
277 count++;
278 }
279 if (count == 0)
280 throw std::runtime_error("mean of empty vector");
281 return sum / static_cast<double>(count);
282}
283
284// calculate stdev for a subset of the vector
285double maskedStdev(const std::vector<double> &vec, const std::vector<bool> &mask) {
286
287 auto avg = maskedMean(vec, mask);
288 size_t count = 0;
289 double sum = 0.0;
290 for (size_t i = 0; i != vec.size(); i++) {
291 if (!mask[i])
292 continue;
293 sum += (vec[i] - avg) * (vec[i] - avg);
294 count++;
295 }
296 return std::sqrt(sum / static_cast<double>(count));
297}
298
299// Calculates a running average for a given half window size assuming the data is wrapped.
300// As the data are integer values, it is possible to maintain an exact sum without any
301// loss of accuracy where the next filtered value is calculated by adding the leading
302// point and subtracting the trailing point from the sum. As the data is wrapped the edge
303// points are handled by the modulus of the array index.
304//
305// |...........[.....x.....].............|
306// 0 |totalWindow| |N
307//
308std::vector<double> runningAverage(const std::vector<size_t> &data, size_t halfWindow) {
309 const auto N = data.size();
310 const size_t totalWindow = 2 * halfWindow + 1;
311
312 // Step 1 is to calculate the sum for the first point where startIndex is the
313 // first point in the total window. The start index wraps from the end of the data.
314 size_t sum{0};
315 size_t startIndex = N - halfWindow;
316 for (size_t i = 0; i < totalWindow; i++) {
317 size_t ix = (startIndex + i) % N;
318 sum += data[ix];
319 }
320
321 // Save the average for the current point, drop the first point from the total window
322 // sum, add the next point to the sum and shift the start index for the window.
323 std::vector<double> filtered(N, 0);
324 for (size_t i = 0; i < N; i++) {
325 // save previous and then move to next
326 filtered[i] = static_cast<double>(sum) / static_cast<double>(totalWindow);
327
328 sum -= data[startIndex];
329 sum += data[(startIndex + totalWindow) % N];
330 startIndex = (startIndex + 1) % N;
331 }
332 return filtered;
333}
334
335// Simple reader that is compatible with the ASNTO event file loader
336class FileLoader {
337 std::ifstream _ifs;
338 size_t _size;
339
340public:
341 explicit FileLoader(const char *filename) : _ifs(filename, std::ios::binary | std::ios::in) {
342 if (!_ifs.is_open() || _ifs.fail())
343 throw std::runtime_error("unable to open file");
344
345 _ifs.seekg(0, _ifs.end);
346 _size = _ifs.tellg();
347 _ifs.seekg(0, _ifs.beg);
348 }
349
350 bool read(char *s, std::streamsize n) { return static_cast<bool>(_ifs.read(s, n)); }
351
352 size_t size() { return _size; }
353
354 size_t position() { return _ifs.tellg(); }
355
356 size_t selected_position() { return _ifs.tellg(); }
357};
358
359} // anonymous namespace
360
361namespace EMU {
362
363// Implement emu specific handlers for the more general EMU events. The main
364// differences with the current impementation of handlers:
365// - tof needs to be derived from the observed time because of the doppler drive
366// - emu includes direct and indirect virtual detectors that account for
367// alternate paths
368// - the loader returns the doppler and auxillary time along with the
369// absolute time rather than just the primary observed time that is
370// equivalent to tof
371//
372// In the future the ANSTO helper and event file loader will be generalized to
373// handle the instruments consistently.
374
376protected:
377 // fields
378 const std::vector<bool> &m_roi;
379 const std::vector<size_t> &m_mapIndex;
380 const size_t m_stride;
381 const double m_framePeriod;
382 const double m_gatePeriod;
383
384 // number of frames
385 size_t m_frames;
387
388 // time boundaries
389 const TimeLimits m_timeBoundary; // seconds
390 const TimeLimits m_directTaux; // microsec
391 const TimeLimits m_analysedTaux; // microsec
392
393 const bool m_includeBM;
394
395 virtual void addEventImpl(size_t id, size_t x, size_t y, double tof) = 0;
396 virtual void addPseudoBMEventImpl(size_t id, double tobs) = 0;
397
398public:
399 EventProcessor(const std::vector<bool> &roi, const std::vector<size_t> &mapIndex, const size_t stride,
400 const double framePeriod, const double gatePeriod, TimeLimits timeBoundary, TimeLimits directLimits,
401 TimeLimits analysedLimits, bool includeBM)
402 : m_roi(roi), m_mapIndex(mapIndex), m_stride(stride), m_framePeriod(framePeriod), m_gatePeriod(gatePeriod),
403 m_frames(0), m_framesValid(0), m_timeBoundary(std::move(timeBoundary)), m_directTaux(std::move(directLimits)),
404 m_analysedTaux(std::move(analysedLimits)), m_includeBM(includeBM) {}
405
406 void newFrame() {
407 m_frames++;
408 if (validFrame())
410 }
411
412 inline bool validFrame() const {
413 double frameTime = m_framePeriod * static_cast<double>(m_frames) * 1.0e-6;
414 return (frameTime >= m_timeBoundary.first && frameTime <= m_timeBoundary.second);
415 }
416
417 double duration() const {
418 // test length in seconds
419 return m_framePeriod * static_cast<double>(m_frames) * 1.0e-6;
420 }
421
422 inline int64_t frameStart() const {
423 // returns time in nanoseconds from start of test
424 auto start = m_framePeriod * static_cast<double>(m_frames);
425 return static_cast<int64_t>(start * 1.0e3);
426 }
427
428 void addEvent(size_t x, size_t p, double tdop, double taux) {
429
430 // check if in time boundaries
431 if (!validFrame())
432 return;
433
434 // group pixels
435 auto y = static_cast<size_t>(p / HISTO_BINS_Y_DENUMERATOR);
436
437 // check for beam monitor tube
438 if (x == PSEUDO_BM_TUBE && y < m_stride) {
439 size_t id = BM_HISTOGRAMS + y;
440 double ptaux = fmod(taux, m_gatePeriod);
441 if (ptaux < 0)
442 ptaux = ptaux + m_gatePeriod;
443 addPseudoBMEventImpl(id, ptaux);
444 return;
445 }
446
447 // determine detector id and check limits
448 if (x >= DETECTOR_TUBES || y >= m_stride)
449 return;
450
451 // map the raw detector index to the physical model
452 size_t xid = m_mapIndex[x];
453
454 // take the modulus of the taux time to account for the
455 // longer background chopper rate
456 double ptaux = fmod(taux, m_gatePeriod);
457 if (ptaux >= m_directTaux.first && ptaux <= m_directTaux.second)
458 xid = xid + DETECTOR_TUBES;
459 else if (!(ptaux >= m_analysedTaux.first && ptaux <= m_analysedTaux.second))
460 return;
461
462 size_t id = m_stride * xid + y;
463 if (id >= m_roi.size())
464 return;
465
466 // check if neutron is in region of interest
467 if (!m_roi[id])
468 return;
469
470 // finally pass to specific handler
471 addEventImpl(id, xid, y, tdop);
472 }
473};
474
476protected:
477 // fields
478 std::vector<size_t> &m_eventCounts;
479
480 void addEventImpl(size_t id, size_t /*x*/, size_t /*y*/, double /*tof*/) override { m_eventCounts[id]++; }
481 void addPseudoBMEventImpl(size_t id, double) override {
482 if (m_includeBM) {
483 m_eventCounts[id]++;
484 }
485 }
486
487public:
488 // construction
489 EventCounter(const std::vector<bool> &roi, const std::vector<size_t> &mapIndex, const size_t stride,
490 const double framePeriod, const double gatePeriod, const TimeLimits &timeBoundary,
491 const TimeLimits &directLimits, const TimeLimits &analysedLimits, std::vector<size_t> &eventCounts,
492 bool includeBM)
493 : EventProcessor(roi, mapIndex, stride, framePeriod, gatePeriod, timeBoundary, directLimits, analysedLimits,
494 includeBM),
495 m_eventCounts(eventCounts) {}
496
497 size_t numFrames() const { return m_framesValid; }
498};
499
501protected:
502 // fields
503 std::vector<EventVector_pt> &m_eventVectors;
504 const ConvertTOF &m_convertTOF;
505 double m_tofMin;
506 double m_tofMax;
507 int64_t m_startTime;
509 const double m_binSize;
510 std::vector<size_t> m_bmCounts;
511
512 void addEventImpl(size_t id, size_t x, size_t /*y*/, double tobs) override {
513
514 // get the absolute time for the start of the frame
515 auto offset = m_startTime + frameStart();
516
517 // convert observation time to tof and set the pulse time
518 // relative to the start of the doppler cycle
519 double tof = tobs;
520
521 if (m_saveAsTOF) {
522 double pulse;
523 if (x < DETECTOR_TUBES)
524 std::tie(pulse, tof) = m_convertTOF.analysedTOF(id, tobs);
525 else
526 std::tie(pulse, tof) = m_convertTOF.directTOF(id, tobs);
527 offset += static_cast<int64_t>(pulse * 1e3);
528 }
529
530 if (m_tofMin > tof)
531 m_tofMin = tof;
532 if (m_tofMax < tof)
533 m_tofMax = tof;
534
535 auto ev = Types::Event::TofEvent(tof, Types::Core::DateAndTime(offset));
536 m_eventVectors[id]->emplace_back(ev);
537 }
538
539 void addPseudoBMEventImpl(size_t id, double tobs) override {
540 // get the absolute time for the start of the frame
541 if (m_includeBM) {
542 auto offset = m_startTime + frameStart();
543 auto ev = Types::Event::TofEvent(tobs, Types::Core::DateAndTime(offset));
544 m_eventVectors[id]->emplace_back(ev);
545 }
546
547 // add to the binned counts
548 auto index = static_cast<size_t>(tobs / m_binSize);
549 m_bmCounts[index] += 1;
550 }
551
552public:
553 EventAssigner(const std::vector<bool> &roi, const std::vector<size_t> &mapIndex, const size_t stride,
554 const double framePeriod, const double gatePeriod, const TimeLimits &timeBoundary,
555 const TimeLimits &directLimits, const TimeLimits &analysedLimits, ConvertTOF &convert,
556 std::vector<EventVector_pt> &eventVectors, int64_t startTime, bool saveAsTOF, bool includeBM)
557 : EventProcessor(roi, mapIndex, stride, framePeriod, gatePeriod, timeBoundary, directLimits, analysedLimits,
558 includeBM),
559 m_eventVectors(eventVectors), m_convertTOF(convert), m_tofMin(std::numeric_limits<double>::max()),
560 m_tofMax(std::numeric_limits<double>::min()), m_startTime(startTime), m_saveAsTOF(saveAsTOF),
561 m_binSize(gatePeriod / BEAM_MONITOR_BINS), m_bmCounts(BEAM_MONITOR_BINS, 0) {}
562
563 double tofMin() const { return m_tofMin <= m_tofMax ? m_tofMin : 0.0; }
564 double tofMax() const { return m_tofMin <= m_tofMax ? m_tofMax : 0.0; }
565 const std::vector<size_t> &beamMonitorCounts() const { return m_bmCounts; }
566 double binSize() const { return m_binSize; }
567 size_t numBins() const { return BEAM_MONITOR_BINS; }
568 size_t bmCounts() const { return std::accumulate(m_bmCounts.begin(), m_bmCounts.end(), (size_t)0); }
569};
570
571template <typename EP>
572void loadEvents(API::Progress &prog, const char *progMsg, const std::string &eventFile, EP &eventProcessor) {
573
574 using namespace ANSTO;
575
576 prog.doReport(progMsg);
577
578 FileLoader loader(eventFile.c_str());
579
580 // for progress notifications
581 ANSTO::ProgressTracker progTracker(prog, progMsg, loader.size(), Progress_LoadBinFile);
582
583 ReadEventFile(loader, eventProcessor, progTracker, 100, false);
584}
585
586} // namespace EMU
587
590template <typename FD> void LoadEMU<FD>::init(bool hdfLoader) {
591
592 // Specify file extensions which can be associated with a specific file.
593 std::vector<std::string> exts;
594
595 // Declare the Filename algorithm property. Mandatory. Sets the path to the
596 // file to load.
597 exts.clear();
598 if (hdfLoader)
599 exts.emplace_back(".hdf");
600 else
601 exts.emplace_back(".tar");
602 Base::declareProperty(std::make_unique<API::FileProperty>(FilenameStr, "", API::FileProperty::Load, exts),
603 "The input filename of the stored data");
604
605 if (hdfLoader) {
606 Base::declareProperty(PathToBinaryStr, "",
607 "Relative or absolute path to the compressed binary\n"
608 "event file linked to the HDF file, eg /storage/data/");
609 }
610
611 // mask
612 exts.clear();
613 exts.emplace_back(".xml");
614 Base::declareProperty(std::make_unique<API::FileProperty>(MaskStr, "", API::FileProperty::OptionalLoad, exts),
615 "The input filename of the mask data");
616
617 Base::declareProperty(SelectDetectorTubesStr, "",
618 "Comma separated range of detectors tubes to be loaded,\n"
619 " eg 16,19-45,47");
620
621 Base::declareProperty(
622 std::make_unique<API::WorkspaceProperty<API::IEventWorkspace>>("OutputWorkspace", "", Kernel::Direction::Output));
623
624 if (hdfLoader) {
625 Base::declareProperty(SelectDatasetStr, 0, "Select the index for the dataset to be loaded.");
626 }
627
628 Base::declareProperty(OverrideDopplerFreqStr, EMPTY_DBL(), "Override the Doppler frequency, in Hertz.");
629
630 Base::declareProperty(OverrideDopplerPhaseStr, EMPTY_DBL(), "Override the Doppler phase, in degrees.");
631
632 Base::declareProperty(CalibrateDopplerPhaseStr, false,
633 "Calibrate the Doppler phase prior to TOF conversion,\n"
634 "ignored if imported as Doppler time or phase entered");
635
636 Base::declareProperty(RawDopplerTimeStr, false,
637 "Import file as observed time relative the Doppler\n"
638 "drive, in microsecs.");
639
640 Base::declareProperty(IncludePseudoBMStr, false, "Include the individual beam monitor events as spectra.");
641
642 Base::declareProperty(FilterByTimeStartStr, 0.0,
643 "Only include events after the provided start time, in "
644 "seconds (relative to the start of the run).");
645
646 Base::declareProperty(FilterByTimeStopStr, EMPTY_DBL(),
647 "Only include events before the provided stop time, in "
648 "seconds (relative to the start of the run).");
649
650 std::string grpOptional = "Filters";
651 Base::setPropertyGroup(FilterByTimeStartStr, grpOptional);
652 Base::setPropertyGroup(FilterByTimeStopStr, grpOptional);
653}
654
656template <typename FD> void LoadEMU<FD>::createWorkspace(const std::string &title) {
657
658 // Create the workspace
659 m_localWorkspace = std::make_shared<DataObjects::EventWorkspace>();
660 m_localWorkspace->initialize(HISTOGRAMS, 2, 1);
661
662 // set the units
663 m_localWorkspace->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("TOF");
664 m_localWorkspace->setYUnit("Counts");
665
666 // set title
667 m_localWorkspace->setTitle(title);
668}
669
678template <typename FD> void LoadEMU<FD>::exec(const std::string &hdfFile, const std::string &eventFile) {
679
680 namespace fs = boost::filesystem;
681
682 // Create workspace
683 // ----------------
684 fs::path p = hdfFile;
685 for (; !p.extension().empty();)
686 p = p.stem();
687 std::string title = p.generic_string();
688 createWorkspace(title);
689 API::LogManager &logManager = m_localWorkspace->mutableRun();
690 API::Progress prog(this, 0.0, 1.0, Progress_Total);
691
692 // Load instrument and workspace properties
693 // ----------------------------------------
694 logManager.addProperty(SelectDatasetStr, m_datasetIndex);
695 loadParameters(hdfFile, logManager);
696 prog.doReport("creating instrument");
697 loadInstrument();
698
699 // Get the region of interest and filters and save to log
700 //
701 std::string maskfile = Base::getPropertyValue(MaskStr);
702 std::string seltubes = Base::getPropertyValue(SelectDetectorTubesStr);
703 logManager.addProperty(SelectDetectorTubesStr, seltubes);
704 logManager.addProperty(MaskStr, maskfile);
705
706 std::vector<bool> roi = createRoiVector(seltubes, maskfile);
707 double timeMaxBoundary = Base::getProperty(FilterByTimeStopStr);
708 if (Base::isEmpty(timeMaxBoundary))
709 timeMaxBoundary = std::numeric_limits<double>::infinity();
710 TimeLimits timeBoundary(Base::getProperty(FilterByTimeStartStr), timeMaxBoundary);
711
712 // lambda to simplify loading instrument parameters
713 auto instr = m_localWorkspace->getInstrument();
714 auto iparam = [&instr](const std::string &tag) { return instr->getNumberParameter(tag)[0]; };
715
716 // Update the neutronic positions for the indirect detectors
717 // noting that the vertical and horizontal are a continuous
718 // sequence starting from 0
719 //
720 double sampleAnalyser = iparam("SampleAnalyser");
721 auto endID = static_cast<detid_t>(DETECTOR_TUBES * PIXELS_PER_TUBE);
722 for (detid_t detID = 0; detID < endID; ++detID)
723 updateNeutronicPostions(detID, sampleAnalyser);
724
725 // get the detector map from raw input to a physical detector
726 //
727 std::string dmapStr = instr->getParameterAsString("DetectorMap");
728 std::vector<size_t> detMapIndex = std::vector<size_t>(DETECTOR_TUBES, 0);
729 mapRangeToIndex(dmapStr, detMapIndex, [](size_t n) { return n; });
730
731 // Collect the L2 distances, Doppler characteristics and
732 // initiate the TOF converter
733 //
734 loadDetectorL2Values();
735 loadDopplerParameters(logManager);
736 double v2 = iparam("AnalysedV2"); // analysed velocity in metres per sec
737 double framePeriod = 1.0e6 / m_dopplerFreq; // period and max direct as microsec
738 double sourceSample = iparam("SourceSample");
739 ConvertTOF convertTOF(m_dopplerAmpl * m_dopplerRun, m_dopplerFreq, m_dopplerPhase, sourceSample, v2, m_detectorL2);
740
741 // Load the events file
742 // --------------------
743 // First count the number of events to reserve memory and then assign the
744 // events to the detectors
745
746 // load events
747 size_t numberHistograms = m_localWorkspace->getNumberHistograms();
748 std::vector<EventVector_pt> eventVectors(numberHistograms, nullptr);
749 std::vector<size_t> eventCounts(numberHistograms, 0);
750
751 // Discriminating between direct and analysed is based on the auxillary time
752 // and is determined by the graphite chopper frequency and v2 which are stable
753 // so these limits are kept in the instrument parameter file. Convert from
754 // milsec to microsec.
755 TimeLimits directLimits(1000.0 * iparam("DirectTauxMin"), 1000.0 * iparam("DirectTauxMax"));
756 TimeLimits analysedLimits(1000.0 * iparam("AnalysedTauxMin"), 1000.0 * iparam("AnalysedTauxMax"));
757
758 // fabs because the value can be negative
759 double gatePeriod = 1.0e6 / fabs(logManager.getTimeSeriesProperty<double>("GraphiteChopperFrequency")->firstValue());
760
761 // count total events per pixel and reserve necessary memory
762 bool includeBM = Base::getProperty(IncludePseudoBMStr);
763 EMU::EventCounter eventCounter(roi, detMapIndex, PIXELS_PER_TUBE, framePeriod, gatePeriod, timeBoundary, directLimits,
764 analysedLimits, eventCounts, includeBM);
765 EMU::loadEvents(prog, "loading neutron counts", eventFile, eventCounter);
766 ANSTO::ProgressTracker progTracker(prog, "creating neutron event lists", numberHistograms, Progress_ReserveMemory);
767 prepareEventStorage(progTracker, eventCounts, eventVectors);
768
769 // now perform the actual event collection and TOF convert if necessary
770 // if a phase calibration is required then load it as raw doppler time
771 // perform the calibration and then convert to TOF
772 Types::Core::DateAndTime startTime(m_startRun);
773 auto start_nanosec = startTime.totalNanoseconds();
774 bool saveAsTOF = !Base::getProperty(RawDopplerTimeStr);
775 bool loadAsTOF = !m_calibrateDoppler && saveAsTOF;
776 EMU::EventAssigner eventAssigner(roi, detMapIndex, PIXELS_PER_TUBE, framePeriod, gatePeriod, timeBoundary,
777 directLimits, analysedLimits, convertTOF, eventVectors, start_nanosec, loadAsTOF,
778 includeBM);
779 EMU::loadEvents(prog, "loading neutron events (TOF)", eventFile, eventAssigner);
780
781 // determine the minimum and maximum beam rate per sec
782 auto filteredBM = runningAverage(eventAssigner.beamMonitorCounts(), BM_HALF_WINDOW);
783 auto res = std::minmax_element(filteredBM.begin(), filteredBM.end());
784 auto ratePerSec = static_cast<double>(eventAssigner.numBins()) / eventCounter.duration();
785 auto minBM = *res.first * ratePerSec;
786 auto maxBM = *res.second * ratePerSec;
787 AddSinglePointTimeSeriesProperty<double>(logManager, m_startRun, "BeamMonitorBkgRate", minBM);
788 AddSinglePointTimeSeriesProperty<double>(logManager, m_startRun, "BeamMonitorRate", maxBM);
789 AddSinglePointTimeSeriesProperty<int>(logManager, m_startRun, "MonitorCounts",
790 static_cast<int>(eventAssigner.bmCounts()));
791
792 // perform a calibration and then TOF conversion if necessary
793 // and update the tof limits
794 auto minTOF = eventAssigner.tofMin();
795 auto maxTOF = eventAssigner.tofMax();
796 if (m_calibrateDoppler) {
797 calibrateDopplerPhase(eventCounts, eventVectors);
798 if (saveAsTOF) {
799 dopplerTimeToTOF(eventVectors, minTOF, maxTOF);
800 }
801 }
802 AddSinglePointTimeSeriesProperty<double>(logManager, m_startRun, "DopplerPhase", m_dopplerPhase);
803
804 // just to make sure the bins hold it all and setup the detector masks
805 m_localWorkspace->setAllX(HistogramData::BinEdges{std::max(0.0, floor(minTOF)), maxTOF + 1});
806 setupDetectorMasks(roi);
807
808 // set log values
809 auto frame_count = eventCounter.numFrames();
810 AddSinglePointTimeSeriesProperty<int>(logManager, m_startRun, "frame_count", static_cast<int>(frame_count));
811
812 // add the scan period in secs to the log
813 auto scan_period = static_cast<double>(frame_count + 1) / m_dopplerFreq;
814 AddSinglePointTimeSeriesProperty<double>(logManager, m_startRun, "ScanPeriod", scan_period);
815
816 std::string filename = Base::getPropertyValue(FilenameStr);
817 logManager.addProperty("filename", filename);
818
819 Types::Core::time_duration duration =
820 boost::posix_time::microseconds(static_cast<boost::int64_t>(eventCounter.duration() * 1.0e6));
821 Types::Core::DateAndTime endTime(startTime + duration);
822 logManager.addProperty("start_time", startTime.toISO8601String());
823 logManager.addProperty("end_time", endTime.toISO8601String());
824 logManager.addProperty<double>("dur", eventCounter.duration());
825
826 // Finally add the time-series parameter explicitly
827 loadEnvironParameters(hdfFile, logManager);
828
829 Base::setProperty("OutputWorkspace", m_localWorkspace);
830}
831
833template <typename FD> void LoadEMU<FD>::setupDetectorMasks(const std::vector<bool> &roi) {
834
835 // count total number of masked bins
836 size_t maskedBins = 0;
837 for (size_t i = 0; i != roi.size(); i++)
838 if (!roi[i])
839 maskedBins++;
840
841 if (maskedBins > 0) {
842 // create list of masked bins
843 std::vector<size_t> maskIndexList(maskedBins);
844 size_t maskIndex = 0;
845
846 for (size_t i = 0; i != roi.size(); i++)
847 if (!roi[i])
848 maskIndexList[maskIndex++] = i;
849
850 auto maskingAlg = Base::createChildAlgorithm("MaskDetectors");
851 maskingAlg->setProperty("Workspace", m_localWorkspace);
852 maskingAlg->setProperty("WorkspaceIndexList", maskIndexList);
853 maskingAlg->executeAsChildAlg();
854 }
855}
856
859template <typename FD>
860void LoadEMU<FD>::prepareEventStorage(ANSTO::ProgressTracker &progTracker, const std::vector<size_t> &eventCounts,
861 std::vector<EventVector_pt> &eventVectors) {
862
863 size_t numberHistograms = eventCounts.size();
864 for (size_t i = 0; i != numberHistograms; ++i) {
865 DataObjects::EventList &eventList = m_localWorkspace->getSpectrum(i);
866
868 eventList.reserve(eventCounts[i]);
869
870 eventList.setDetectorID(static_cast<detid_t>(i));
871 eventList.setSpectrumNo(static_cast<detid_t>(i));
872
873 DataObjects::getEventsFrom(eventList, eventVectors[i]);
874
875 progTracker.update(i);
876 }
877 progTracker.complete();
878}
879
881template <typename FD> void LoadEMU<FD>::loadDopplerParameters(API::LogManager &logm) {
882
883 auto instr = m_localWorkspace->getInstrument();
884
885 // use nominal frequency based on amp and velocity unless overridden
886 m_dopplerAmpl = logm.getTimeSeriesProperty<double>("DopplerAmplitude")->firstValue();
887 m_dopplerRun = logm.getTimeSeriesProperty<int32_t>("DopplerRun")->firstValue();
888 m_dopplerFreq = Base::getProperty(OverrideDopplerFreqStr);
889 if (Base::isEmpty(m_dopplerFreq)) {
890 auto doppVel = logm.getTimeSeriesProperty<double>("DopplerVelocity")->firstValue();
891 m_dopplerFreq = 0.5 * doppVel / (M_PI * m_dopplerAmpl);
892 }
893 AddSinglePointTimeSeriesProperty<double>(logm, m_startRun, "DopplerFrequency", m_dopplerFreq);
894
895 m_dopplerPhase = Base::getProperty(OverrideDopplerPhaseStr);
896 m_calibrateDoppler = Base::getProperty(CalibrateDopplerPhaseStr) && Base::isEmpty(m_dopplerPhase);
897 if (Base::isEmpty(m_dopplerPhase)) {
898 // sinusoidal motion crossing a threshold with a delay
899 double doppThreshold = instr->getNumberParameter("DopplerReferenceThreshold")[0];
900 double doppDelay = instr->getNumberParameter("DopplerReferenceDelay")[0];
901 m_dopplerPhase = 180.0 - asin(0.001 * doppThreshold / m_dopplerAmpl) * 180.0 / M_PI + doppDelay * m_dopplerFreq;
902 }
903
904 // problem adding 'bool' to log
905 int32_t calPhase = (m_calibrateDoppler ? 1 : 0);
906 logm.addProperty("CalibratePhase", calPhase);
907}
908
911template <typename FD>
912void LoadEMU<FD>::calibrateDopplerPhase(const std::vector<size_t> &eventCounts,
913 const std::vector<EventVector_pt> &eventVectors) {
914
915 // get the doppler parameters
916 auto instr = m_localWorkspace->getInstrument();
917 double v2 = instr->getNumberParameter("AnalysedV2")[0];
918 double l1 = instr->getNumberParameter("SourceSample")[0];
919
920 // get the number of analysed events and initial doppler time
921 auto startID = static_cast<size_t>(HORIZONTAL_TUBES * PIXELS_PER_TUBE);
922 auto endID = static_cast<size_t>(DETECTOR_TUBES * PIXELS_PER_TUBE);
923 size_t numEvents = 0;
924 for (size_t i = startID; i < endID; i++)
925 numEvents += eventCounts[i];
926 if (numEvents == 0)
927 throw std::runtime_error("no analysed events for phase calibration");
928 std::vector<double> nVel(numEvents);
929 std::vector<size_t> nMap(numEvents);
930 std::vector<bool> nCnd(numEvents);
931 constexpr size_t NHIST = 100;
932 std::vector<int> histogram(NHIST + 1, 0);
933
934 // define the cost function to optimize phase
935 auto costFn = [&, this](double phase) {
936 ConvertTOF convTOF(m_dopplerAmpl * m_dopplerRun, m_dopplerFreq, phase, l1, v2, m_detectorL2);
937
938 // convert each analysed event to source velocity
939 size_t ix = 0;
940 double tof, pulse;
941 for (size_t i = startID; i < endID; i++) {
942 for (auto const &x : *eventVectors[i]) {
943 std::tie(pulse, tof) = convTOF.analysedTOF(i, x.tof());
944 auto tof1 = 1e-6 * tof - m_detectorL2[i] / v2;
945 nVel[ix++] = l1 / tof1;
946 }
947 }
948
949 // now histogram the data and create the map from velocity to hist
950 auto ixlim = std::minmax_element(nVel.begin(), nVel.end());
951 auto vmin = nVel[ixlim.first - nVel.begin()];
952 auto vmax = nVel[ixlim.second - nVel.begin()];
953 int maxHist = 0;
954 std::fill(histogram.begin(), histogram.end(), 0);
955 auto delta = (vmax - vmin) / NHIST;
956 for (size_t i = 0; i < numEvents; i++) {
957 auto v = nVel[i];
958 auto j = static_cast<size_t>(std::floor((v - vmin) / delta));
959 histogram[j]++;
960 if (histogram[j] > maxHist)
961 maxHist = histogram[j];
962 nMap[i] = j;
963 }
964
965 // determine the points above the 25% threshold
966 auto minLevel = static_cast<int>(maxHist / 4);
967 for (size_t i = 0; i < numEvents; i++) {
968 nCnd[i] = (histogram[nMap[i]] >= minLevel ? true : false);
969 }
970
971 // calculate the standard deviation for the points above the threshold
972 auto cost = maskedStdev(nVel, nCnd);
973 return cost;
974 };
975
976 // call the optimizer and update the doppler phase value
977 // limit the optimization to 30 iterations
978 int bits = std::numeric_limits<double>::digits;
979 boost::uintmax_t itn = 30;
980 using boost::math::tools::brent_find_minima;
981 auto minPhase = m_dopplerPhase - 5.0;
982 auto maxPhase = m_dopplerPhase + 5.0;
983 auto r = brent_find_minima(costFn, minPhase, maxPhase, bits, itn);
984 m_dopplerPhase = r.first;
985}
986
989template <typename FD>
990void LoadEMU<FD>::dopplerTimeToTOF(std::vector<EventVector_pt> &eventVectors, double &minTOF, double &maxTOF) {
991
992 // get the doppler parameters and initialise TOD converter
993 auto instr = m_localWorkspace->getInstrument();
994 double v2 = instr->getNumberParameter("AnalysedV2")[0];
995 double l1 = instr->getNumberParameter("SourceSample")[0];
996 ConvertTOF convTOF(m_dopplerAmpl * m_dopplerRun, m_dopplerFreq, m_dopplerPhase, l1, v2, m_detectorL2);
997
998 // run through all the events noting that analysed event are in
999 // the bottom half of the detector ids
1000 auto start = true;
1001 auto directID = static_cast<size_t>(DETECTOR_TUBES * PIXELS_PER_TUBE);
1002 for (size_t id = 0; id < eventVectors.size(); id++) {
1003 for (auto &x : *eventVectors[id]) {
1004 double tof, pulse;
1005 if (id < directID)
1006 std::tie(pulse, tof) = convTOF.analysedTOF(id, x.tof());
1007 else
1008 std::tie(pulse, tof) = convTOF.directTOF(id, x.tof());
1009
1010 // update the pulse time and tof
1011 int64_t pulseTime = x.pulseTime().totalNanoseconds();
1012 pulseTime += static_cast<int64_t>(pulse * 1000);
1013 x = Types::Event::TofEvent(tof, Types::Core::DateAndTime(pulseTime));
1014
1015 if (start) {
1016 minTOF = maxTOF = x.tof();
1017 start = false;
1018 } else {
1019 minTOF = std::min(minTOF, x.tof());
1020 maxTOF = std::max(maxTOF, x.tof());
1021 }
1022 }
1023 }
1024}
1025
1027template <typename FD> void LoadEMU<FD>::loadDetectorL2Values() {
1028
1029 m_detectorL2 = std::vector<double>(HISTOGRAMS);
1030 const auto &detectorInfo = m_localWorkspace->detectorInfo();
1031 auto detIDs = detectorInfo.detectorIDs();
1032 for (const auto detID : detIDs) {
1033 auto ix = detectorInfo.indexOf(detID);
1034 double l2 = detectorInfo.l2(ix);
1035 m_detectorL2[detID] = l2;
1036 }
1037}
1038
1041template <typename FD> void LoadEMU<FD>::updateNeutronicPostions(detid_t detID, double sampleAnalyser) {
1042
1043 Geometry::Instrument_const_sptr instrument = m_localWorkspace->getInstrument();
1044 auto &compInfo = m_localWorkspace->mutableComponentInfo();
1045
1046 try {
1047 auto component = instrument->getDetector(detID);
1048 double rho, theta, phi;
1049 V3D position = component->getPos();
1050 position.getSpherical(rho, theta, phi);
1051
1052 double scale = -(2 * sampleAnalyser + rho) / rho;
1053 position *= scale;
1054
1055 const auto componentIndex = compInfo.indexOf(component->getComponentID());
1056 compInfo.setPosition(componentIndex, position);
1057 } catch (const std::runtime_error &) {
1058 // just continue with the remainder
1059 }
1060}
1061
1064template <typename FD>
1065std::vector<bool> LoadEMU<FD>::createRoiVector(const std::string &selected, const std::string &maskfile) {
1066
1067 std::vector<bool> result(HISTOGRAMS, true);
1068
1069 // turn off pixels linked to missing tubes
1070 if (!selected.empty()) {
1071 std::vector<bool> tubes(HISTO_BINS_X, false);
1072 mapRangeToIndex(selected, tubes, [](size_t) { return true; });
1073 for (size_t i = 0; i < HISTO_BINS_X; i++) {
1074 if (tubes[i] == false) {
1075 for (size_t j = 0; j < PIXELS_PER_TUBE; j++) {
1076 result[i * PIXELS_PER_TUBE + j] = false;
1077 }
1078 }
1079 }
1080 }
1081
1082 if (maskfile.length() == 0)
1083 return result;
1084
1085 std::ifstream input(maskfile.c_str());
1086 if (!input.good())
1087 throw std::invalid_argument("invalid mask file");
1088
1089 std::string line;
1090 while (std::getline(input, line)) {
1091 auto i0 = line.find("<detids>");
1092 auto iN = line.find("</detids>");
1093
1094 if ((i0 != std::string::npos) && (iN != std::string::npos) && (i0 < iN)) {
1095 line = line.substr(i0 + 8, iN - i0 - 8); // 8 = len("<detids>")
1096 mapRangeToIndex(line, result, [](size_t) { return false; });
1097 }
1098 }
1099
1100 return result;
1101}
1102
1104template <typename FD> void LoadEMU<FD>::loadParameters(const std::string &hdfFile, API::LogManager &logm) {
1105
1106 NeXus::NXRoot root(hdfFile);
1107 NeXus::NXEntry entry = root.openFirstEntry();
1108
1109 MapNeXusToProperty<std::string>(entry, "sample/name", "unknown", logm, "SampleName", "", 0);
1110 MapNeXusToProperty<std::string>(entry, "sample/description", "unknown", logm, "SampleDescription", "", 0);
1111
1112 // if dataset index > 0 need to add an offset to the start time
1113 Types::Core::DateAndTime startTime(GetNeXusValue<std::string>(entry, "start_time", "2000-01-01T00:00:00", 0));
1114 if (m_datasetIndex > 0) {
1115 auto baseTime = GetNeXusValue<int32_t>(entry, "instrument/detector/start_time", 0, 0);
1116 auto nthTime = GetNeXusValue<int32_t>(entry, "instrument/detector/start_time", 0, m_datasetIndex);
1117
1118 Types::Core::time_duration duration =
1119 boost::posix_time::microseconds(static_cast<boost::int64_t>((nthTime - baseTime) * 1.0e6));
1120 Types::Core::DateAndTime startDataset(startTime + duration);
1121 m_startRun = startDataset.toISO8601String();
1122 } else {
1123 m_startRun = startTime.toISO8601String();
1124 }
1125
1126 MapNeXusToSeries<double>(entry, "instrument/doppler/ctrl/amplitude", 75.0, logm, m_startRun, "DopplerAmplitude",
1127 0.001, m_datasetIndex);
1128 MapNeXusToSeries<double>(entry, "instrument/doppler/ctrl/velocity", 4.7, logm, m_startRun, "DopplerVelocity", 1,
1129 m_datasetIndex);
1130 MapNeXusToSeries<int>(entry, "instrument/doppler/ctrl/run_cmd", 1, logm, m_startRun, "DopplerRun", 1, m_datasetIndex);
1131
1132 MapNeXusToSeries<double>(entry, "instrument/chpr/background/actspeed", 1272.8, logm, m_startRun,
1133 "BackgroundChopperFrequency", 1.0 / 60, 0);
1134 MapNeXusToSeries<double>(entry, "instrument/chpr/graphite/actspeed", 2545.6, logm, m_startRun,
1135 "GraphiteChopperFrequency", 1.0 / 60, 0);
1136 // hz tube gap or equivalent to be added later - reverts to default
1137 MapNeXusToSeries<double>(entry, "instrument/hztubegap", 0.02, logm, m_startRun, "horizontal_tubes_gap", 1.0, 0);
1138 // add the reactor power to the log
1139 MapNeXusToSeries<double>(entry, "instrument/source/power", 20.0, logm, m_startRun, "ReactorPower", 1.0,
1140 m_datasetIndex);
1141 // fix for source position when loading IDF
1142 MapNeXusToProperty<double>(entry, "instrument/doppler/tosource", 2.035, logm, "SourceSample", 1.0, 0);
1143}
1144
1147template <typename FD> void LoadEMU<FD>::loadEnvironParameters(const std::string &hdfFile, API::LogManager &logm) {
1148
1149 NeXus::NXRoot root(hdfFile);
1150 NeXus::NXEntry entry = root.openFirstEntry();
1151 auto time_str = logm.getPropertyValueAsType<std::string>("end_time");
1152
1153 // load the environment variables for the dataset loaded
1154 std::vector<std::string> tags = {"P01PS03", "P01PSP03", "T01S00", "T01S05", "T01S06", "T01SP00", "T01SP06",
1155 "T02S00", "T02S04", "T02S05", "T02S06", "T02SP00", "T02SP06", "T3S1",
1156 "T3S2", "T3S3", "T3S4", "T3SP1", "T3SP2", "T3SP3", "T3SP4"};
1157
1158 for (const auto &tag : tags) {
1159 MapNeXusToSeries<double>(entry, "data/" + tag, 0.0, logm, time_str, "env_" + tag, 1.0, m_datasetIndex);
1160 }
1161}
1162
1164template <typename FD> void LoadEMU<FD>::loadInstrument() {
1165
1166 // loads the IDF and parameter file
1167 auto loadInstrumentAlg = Base::createChildAlgorithm("LoadInstrument");
1168 loadInstrumentAlg->setProperty("Workspace", m_localWorkspace);
1169 loadInstrumentAlg->setPropertyValue("InstrumentName", "EMUau");
1170 loadInstrumentAlg->setProperty("RewriteSpectraMap", Mantid::Kernel::OptionalBool(false));
1171 loadInstrumentAlg->executeAsChildAlg();
1172}
1173
1174//--------- explicit template instantiation -----------
1175
1176// Instantiate base class for LoadEMU's
1177template class LoadEMU<Kernel::FileDescriptor>;
1179
1180// -------- EMU Hdf loader -----------------------
1181
1183int LoadEMUHdf::version() const { return 1; }
1184
1186const std::vector<std::string> LoadEMUHdf::seeAlso() const { return {"Load", "LoadQKK"}; }
1188const std::string LoadEMUHdf::category() const { return "DataHandling\\ANSTO"; }
1189
1191const std::string LoadEMUHdf::name() const { return "LoadEMUHdf"; }
1192
1194const std::string LoadEMUHdf::summary() const { return "Loads an EMU Hdf and linked event file into a workspace."; }
1195
1199 if (descriptor.extension() != ".hdf")
1200 return 0;
1201
1202 if (descriptor.pathExists("/entry1/site_name") && descriptor.pathExists("/entry1/instrument/doppler/ctrl/velocity") &&
1203 descriptor.pathExists("/entry1/instrument/doppler/ctrl/amplitude") &&
1204 descriptor.pathExists("/entry1/instrument/detector/daq_dirname") &&
1205 descriptor.pathExists("/entry1/instrument/detector/dataset_number") &&
1206 descriptor.pathExists("/entry1/data/hmm_total_t_ds0") && descriptor.pathExists("/entry1/data/hmm_total_t_ds1") &&
1207 descriptor.pathExists("/entry1/data/hmm_total_xt_ds0") &&
1208 descriptor.pathExists("/entry1/data/hmm_total_xt_ds1")) {
1209 return 80;
1210 } else {
1211 return 0;
1212 }
1213}
1214
1218
1221// exec() function that works with the two files.
1223
1224 namespace fs = boost::filesystem;
1225
1226 // Open the hdf file and find the dirname and dataset number
1227 std::string hdfFile = Base::getPropertyValue(FilenameStr);
1228 std::string evtPath = Base::getPropertyValue(PathToBinaryStr);
1229 if (evtPath.empty())
1230 evtPath = "./";
1231
1232 // if relative ./ or ../ then append to the directory for the hdf file
1233 if (evtPath.rfind("./") == 0 || evtPath.rfind("../") == 0) {
1234 fs::path hp = hdfFile;
1235 evtPath = fs::canonical(evtPath, hp.parent_path()).generic_string();
1236 }
1237
1238 // dataset index to be loaded
1239 m_datasetIndex = Base::getProperty(SelectDatasetStr);
1240
1241 // if path provided build the file path
1242 if (fs::is_directory(evtPath)) {
1243 NeXus::NXRoot root(hdfFile);
1244 NeXus::NXEntry entry = root.openFirstEntry();
1245 auto eventDir = GetNeXusValue<std::string>(entry, "instrument/detector/daq_dirname", "./", 0);
1246 auto dataset = GetNeXusValue<int32_t>(entry, "instrument/detector/dataset_number", 0, m_datasetIndex);
1247
1248 // build the path to the event file as if a relative or absolute
1249 // path is passed:
1250 // 'relpath/[daq_dirname]/DATASET_[n]/EOS.bin' or the
1251 char buffer[255] = {};
1252 snprintf(buffer, sizeof(buffer), "%s/DATASET_%d/EOS.bin", eventDir.c_str(), dataset);
1253 fs::path path = evtPath;
1254 path /= buffer;
1255 path = fs::absolute(path);
1256 evtPath = path.generic_string();
1257 }
1258
1259 // finally check that the event file exists
1260 if (!fs::is_regular_file(evtPath)) {
1261 std::string msg = "Check path, cannot open binary event file: " + evtPath;
1262 throw std::runtime_error(msg);
1263 }
1264
1266}
1267
1268// -------- EMU Tar loader -----------------------
1269
1271int LoadEMUTar::version() const { return 1; }
1272
1274const std::vector<std::string> LoadEMUTar::seeAlso() const { return {"Load", "LoadQKK"}; }
1276const std::string LoadEMUTar::category() const { return "DataHandling\\ANSTO"; }
1277
1279const std::string LoadEMUTar::name() const { return "LoadEMU"; }
1280
1282const std::string LoadEMUTar::summary() const {
1283 return "Loads an EMU tar file, containing the Hdf and event file, into a "
1284 "workspace.";
1285}
1286
1290 if (descriptor.extension() != ".tar")
1291 return 0;
1292
1293 ANSTO::Tar::File file(descriptor.filename());
1294 if (!file.good())
1295 return 0;
1296
1297 size_t hdfFiles = 0;
1298 size_t binFiles = 0;
1299 const std::vector<std::string> &subFiles = file.files();
1300 for (const auto &subFile : subFiles) {
1301 auto len = subFile.length();
1302 if ((len > 4) && (subFile.find_first_of("\\/", 0, 2) == std::string::npos)) {
1303 if ((subFile.rfind(".hdf") == len - 4) && (subFile.compare(0, 3, "EMU") == 0))
1304 hdfFiles++;
1305 else if (subFile.rfind(".bin") == len - 4)
1306 binFiles++;
1307 }
1308 }
1309
1310 return (hdfFiles == 1) && (binFiles == 1) ? 50 : 0;
1311}
1312
1316
1321
1322 // Opens the tar and extracts the hdf and event data to temporary files
1323 std::string filename = Base::getPropertyValue(FilenameStr);
1324 ANSTO::Tar::File tarFile(filename);
1325 if (!tarFile.good())
1326 throw std::invalid_argument("invalid EMU tar file");
1327
1328 // dataset selection not supported in tar version - order is not guaranteed
1329 m_datasetIndex = 0;
1330
1331 // lambda functions to find the first file of extension and to extract
1332 // the file
1333 const std::vector<std::string> &files = tarFile.files();
1334 auto selectFile = [&](const std::string &ext) {
1335 auto itf = std::find_if(files.cbegin(), files.cend(),
1336 [&ext](const std::string &file) { return file.rfind(ext) == file.length() - 4; });
1337 if (itf == files.end())
1338 throw std::runtime_error("missing tar file data");
1339 else
1340 tarFile.select(itf->c_str());
1341 };
1342 auto extractFile = [&](Poco::TemporaryFile &tfile) {
1343 std::shared_ptr<FILE> handle(fopen(tfile.path().c_str(), "wb"), fclose);
1344 if (handle) {
1345 // copy content
1346 char buffer[4096];
1347 size_t bytesRead;
1348 while (0 != (bytesRead = tarFile.read(buffer, sizeof(buffer))))
1349 fwrite(buffer, bytesRead, 1, handle.get());
1350 handle.reset();
1351 }
1352 };
1353
1354 // extract hdf file into tmp file
1355 selectFile(".hdf");
1356 Poco::TemporaryFile hdfFile;
1357 extractFile(hdfFile);
1358
1359 // extract the event file
1360 selectFile(".bin");
1361 Poco::TemporaryFile eventFile;
1362 extractFile(eventFile);
1363
1364 // call the common loader
1365 LoadEMU<Kernel::FileDescriptor>::exec(hdfFile.path(), eventFile.path());
1366}
1367
1368// register the algorithms into the AlgorithmFactory
1371
1372} // namespace Mantid::DataHandling
double value
The value of the point.
Definition: FitMW.cpp:51
double position
Definition: GetAllEi.cpp:154
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
const double m_L0
Definition: LoadEMU.cpp:228
size_t _size
Definition: LoadEMU.cpp:338
const double m_A
Definition: LoadEMU.cpp:230
std::ifstream _ifs
Definition: LoadEMU.cpp:337
const double m_phi
Definition: LoadEMU.cpp:227
const double m_w
Definition: LoadEMU.cpp:226
const double m_v2
Definition: LoadEMU.cpp:229
const std::vector< double > & m_L2
Definition: LoadEMU.cpp:231
#define fabs(x)
Definition: Matrix.cpp:22
int count
counter
Definition: Matrix.cpp:37
#define DECLARE_NEXUS_FILELOADER_ALGORITHM(classname)
DECLARE_NEXUS_FILELOADER_ALGORITHM should be used in place of the standard DECLARE_ALGORITHM macro wh...
#define DECLARE_FILELOADER_ALGORITHM(classname)
DECLARE_FILELOADER_ALGORITHM should be used in place of the standard DECLARE_ALGORITHM macro when wri...
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
virtual void init()=0
Virtual method - must be overridden by concrete algorithm.
virtual void exec()=0
Virtual method - must be overridden by concrete algorithm.
@ OptionalLoad
to specify a file to read but the file doesn't have to exist
Definition: FileProperty.h:53
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
void setDetectorID(const detid_t detID)
Clear the list of detector IDs, then add one.
Definition: ISpectrum.cpp:84
void setSpectrumNo(specnum_t num)
Sets the the spectrum number of this spectrum.
Definition: ISpectrum.cpp:127
This class contains the information about the log entries.
Definition: LogManager.h:44
void addProperty(Kernel::Property *prop, bool overwrite=false)
Add data to the object in the form of a property.
Definition: LogManager.h:79
HeldType getPropertyValueAsType(const std::string &name) const
Get the value of a property as the given TYPE.
Definition: LogManager.cpp:332
Kernel::TimeSeriesProperty< T > * getTimeSeriesProperty(const std::string &name) const
Returns a property as a time series property.
Definition: LogManager.cpp:308
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
void doReport(const std::string &msg="") override
Actually do the reporting, without changing the loop counter.
Definition: Progress.cpp:70
A property class for workspaces.
helper class to keep track of progress
size_t read(void *dst, size_t size)
const std::vector< std::string > & files() const
std::vector< EventVector_pt > & m_eventVectors
Definition: LoadEMU.cpp:503
const std::vector< size_t > & beamMonitorCounts() const
Definition: LoadEMU.cpp:565
EventAssigner(const std::vector< bool > &roi, const std::vector< size_t > &mapIndex, const size_t stride, const double framePeriod, const double gatePeriod, const TimeLimits &timeBoundary, const TimeLimits &directLimits, const TimeLimits &analysedLimits, ConvertTOF &convert, std::vector< EventVector_pt > &eventVectors, int64_t startTime, bool saveAsTOF, bool includeBM)
Definition: LoadEMU.cpp:553
void addEventImpl(size_t id, size_t x, size_t, double tobs) override
Definition: LoadEMU.cpp:512
void addPseudoBMEventImpl(size_t id, double tobs) override
Definition: LoadEMU.cpp:539
EventCounter(const std::vector< bool > &roi, const std::vector< size_t > &mapIndex, const size_t stride, const double framePeriod, const double gatePeriod, const TimeLimits &timeBoundary, const TimeLimits &directLimits, const TimeLimits &analysedLimits, std::vector< size_t > &eventCounts, bool includeBM)
Definition: LoadEMU.cpp:489
void addEventImpl(size_t id, size_t, size_t, double) override
Definition: LoadEMU.cpp:480
void addPseudoBMEventImpl(size_t id, double) override
Definition: LoadEMU.cpp:481
std::vector< size_t > & m_eventCounts
Definition: LoadEMU.cpp:478
const std::vector< bool > & m_roi
Definition: LoadEMU.cpp:378
void addEvent(size_t x, size_t p, double tdop, double taux)
Definition: LoadEMU.cpp:428
const std::vector< size_t > & m_mapIndex
Definition: LoadEMU.cpp:379
virtual void addEventImpl(size_t id, size_t x, size_t y, double tof)=0
EventProcessor(const std::vector< bool > &roi, const std::vector< size_t > &mapIndex, const size_t stride, const double framePeriod, const double gatePeriod, TimeLimits timeBoundary, TimeLimits directLimits, TimeLimits analysedLimits, bool includeBM)
Definition: LoadEMU.cpp:399
virtual void addPseudoBMEventImpl(size_t id, double tobs)=0
LoadEMUHdf : Loads an ANSTO EMU Hdf and linked event file into a workspace.
Definition: LoadEMU.h:157
const std::vector< std::string > seeAlso() const override
Similar algorithms.
Definition: LoadEMU.cpp:1186
void exec() override
Execute the algorithm.
Definition: LoadEMU.cpp:1222
int version() const override
Algorithm's version for identification.
Definition: LoadEMU.cpp:1183
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
Definition: LoadEMU.cpp:1194
const std::string name() const override
Algorithms name for identification.
Definition: LoadEMU.cpp:1191
const std::string category() const override
Algorithm's category for identification.
Definition: LoadEMU.cpp:1188
void init() override
Initialise the algorithm and declare the properties for the nexus descriptor.
Definition: LoadEMU.cpp:1217
int confidence(Kernel::NexusDescriptor &descriptor) const override
Return the confidence as an integer value that this algorithm can load the file descriptor.
Definition: LoadEMU.cpp:1198
LoadEMUTar : Loads a merged ANSTO EMU Hdf and event file into a workspace.
Definition: LoadEMU.h:120
int version() const override
Algorithm's version for identification.
Definition: LoadEMU.cpp:1271
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
Definition: LoadEMU.cpp:1282
const std::string name() const override
Algorithms name for identification.
Definition: LoadEMU.cpp:1279
int confidence(Kernel::FileDescriptor &descriptor) const override
Return the confidence as an integer value that this algorithm can load the file descriptor.
Definition: LoadEMU.cpp:1289
const std::string category() const override
Algorithm's category for identification.
Definition: LoadEMU.cpp:1276
void exec() override
Execute the algorithm.
Definition: LoadEMU.cpp:1320
const std::vector< std::string > seeAlso() const override
Similar algorithms.
Definition: LoadEMU.cpp:1274
void init() override
Initialise the algorithm and declare the standard properties for the general file descriptor.
Definition: LoadEMU.cpp:1315
void prepareEventStorage(ANSTO::ProgressTracker &prog, const std::vector< size_t > &eventCounts, std::vector< EventVector_pt > &eventVectors)
Allocate space for the event storage in eventVectors after the eventCounts have been determined.
Definition: LoadEMU.cpp:860
void calibrateDopplerPhase(const std::vector< size_t > &eventCounts, const std::vector< EventVector_pt > &eventVectors)
Calibrate the doppler phase based on the analysed events using the eventCounts and eventVectors.
Definition: LoadEMU.cpp:912
std::vector< bool > createRoiVector(const std::string &seltubes, const std::string &maskfile)
Region of interest is defined by the selected detectors and the maskfile.
Definition: LoadEMU.cpp:1065
void loadInstrument()
Load the instrument definition.
Definition: LoadEMU.cpp:1164
void loadParameters(const std::string &hdfFile, API::LogManager &logm)
Load parameters from input hdfFile and save to the log manager, logm.
Definition: LoadEMU.cpp:1104
void dopplerTimeToTOF(std::vector< EventVector_pt > &eventVectors, double &minTOF, double &maxTOF)
Convert the doppler time to TOF for all the events in eventVectors and time of flight range as minTOF...
Definition: LoadEMU.cpp:990
void loadDetectorL2Values()
Recovers the L2 neutronic distance for each detector.
Definition: LoadEMU.cpp:1027
void setupDetectorMasks(const std::vector< bool > &roi)
Set up the detector masks to the region of interest roi.
Definition: LoadEMU.cpp:833
void updateNeutronicPostions(detid_t detID, double sampleAnalyser)
Update the neutronic position for the detID using the distance from the source and the sample to anal...
Definition: LoadEMU.cpp:1041
void loadDopplerParameters(API::LogManager &logm)
Get the Doppler parameters and record to the log manager, logm.
Definition: LoadEMU.cpp:881
void createWorkspace(const std::string &title)
Creates an event workspace and sets the title.
Definition: LoadEMU.cpp:656
void loadEnvironParameters(const std::string &hdfFile, API::LogManager &logm)
Load the environment variables from the hdfFile and save as time series to the log manager,...
Definition: LoadEMU.cpp:1147
A class for holding :
Definition: EventList.h:56
void setSortOrder(const EventSortType order) const
Manually set the event list sort order value.
Definition: EventList.cpp:943
void reserve(size_t num) override
Reserve a certain number of entries in event list of the specified eventType.
Definition: EventList.cpp:895
Defines a wrapper around an open file.
const std::string & filename() const
Access the filename.
const std::string & extension() const
Access the file extension.
Defines a wrapper around a file whose internal structure can be accessed using the NeXus API.
const std::string & extension() const
Access the file extension.
bool pathExists(const std::string &path) const
Query if a path exists.
OptionalBool : Tri-state bool.
Definition: OptionalBool.h:25
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
A specialised Property class for holding a series of time-value pairs.
Class for 3D vectors.
Definition: V3D.h:34
NXDataSetTyped< T > openNXDataSet(const std::string &name) const
Templated method for creating datasets.
Definition: NexusClasses.h:536
Templated class implementation of NXDataSet.
Definition: NexusClasses.h:203
void load(const int blocksize=1, int i=-1, int j=-1, int k=-1, int l=-1) override
Implementation of the virtual NXDataSet::load(...) method.
Definition: NexusClasses.h:289
Implements NXentry Nexus class.
Definition: NexusClasses.h:898
Implements NXroot Nexus class.
Definition: NexusClasses.h:922
NXEntry openFirstEntry()
Open the first NXentry in the file.
std::shared_ptr< T > createWorkspace(InitArgs... args)
void loadEvents(API::Progress &prog, const char *progMsg, const std::string &eventFile, EP &eventProcessor)
Definition: LoadEMU.cpp:572
static char const *const FilenameStr
Definition: LoadBBY.cpp:54
static const size_t HISTO_BINS_X
Definition: LoadBBY.cpp:47
static const size_t HISTO_BINS_Y
Definition: LoadBBY.cpp:48
static const size_t Progress_LoadBinFile
Definition: LoadBBY.cpp:50
void AddSinglePointTimeSeriesProperty(API::LogManager &logManager, const std::string &time, const std::string &name, const TYPE value)
Definition: LoadBBY.cpp:66
std::size_t numEvents(::NeXus::File &file, bool &hasTotalCounts, bool &oldNeXusFileNames, const std::string &prefix, const NexusHDF5Descriptor &descriptor)
Get the number of events in the currently opened group.
static char const *const MaskStr
Definition: LoadBBY.cpp:55
static char const *const FilterByTimeStartStr
Definition: LoadBBY.cpp:60
static const size_t Progress_Total
Definition: LoadBBY.cpp:52
static char const *const FilterByTimeStopStr
Definition: LoadBBY.cpp:61
static const size_t Progress_ReserveMemory
Definition: LoadBBY.cpp:51
DLLExport void getEventsFrom(EventList &el, std::vector< Types::Event::TofEvent > *&events)
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
NXDataSetTyped< char > NXChar
The char dataset type.
Definition: NexusClasses.h:469
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
STL namespace.
@ Output
An output workspace.
Definition: Property.h:54