Mantid
Loading...
Searching...
No Matches
LoadPLN.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
9#include "MantidAPI/Axis.h"
13#include "MantidAPI/Run.h"
23
24#include <boost/filesystem.hpp>
25
26#include <algorithm>
27#include <cmath>
28#include <fstream>
29#include <utility>
30
31namespace Mantid::DataHandling {
32
33using namespace Kernel;
34
35namespace {
36
37// number of physical detectors
38constexpr size_t MONITORS = 8;
39constexpr size_t DETECTOR_TUBES = 200;
40constexpr size_t HISTO_BINS_X = DETECTOR_TUBES + MONITORS;
41constexpr size_t HISTO_BINS_Y_DENUMERATOR = 16;
42constexpr size_t PIXELS_PER_TUBE = 1024 / HISTO_BINS_Y_DENUMERATOR;
43constexpr size_t DETECTOR_SPECTRA = DETECTOR_TUBES * PIXELS_PER_TUBE;
44constexpr size_t HISTOGRAMS = DETECTOR_SPECTRA + MONITORS;
45
46// File loading progress boundaries
47constexpr size_t Progress_LoadBinFile = 48;
48constexpr size_t Progress_ReserveMemory = 4;
50
51// Algorithm parameter names
52constexpr char FilenameStr[] = "Filename";
53constexpr char MaskStr[] = "Mask";
54constexpr char SelectDetectorTubesStr[] = "SelectDetectorTubes";
55constexpr char SelectDatasetStr[] = "SelectDataset";
56constexpr char FilterByTimeStartStr[] = "FilterByTimeStart";
57constexpr char FilterByTimeStopStr[] = "FilterByTimeStop";
58constexpr char PathToBinaryStr[] = "BinaryEventPath";
59constexpr char TOFBiasStr[] = "TimeOfFlightBias";
60constexpr char CalibrateTOFStr[] = "CalibrateTOFBias";
61constexpr char LambdaOnTwoStr[] = "LambdaOnTwoMode";
62
63// Common pairing of limits
64using TimeLimits = std::pair<double, double>;
65
66template <typename Type>
67void AddSinglePointTimeSeriesProperty(API::LogManager &logManager, const std::string &time, const std::string &name,
68 const Type value) {
69 // create time series property and add single value
70 auto p = new Kernel::TimeSeriesProperty<Type>(name);
71 p->addValue(time, value);
72
73 // add to log manager
74 logManager.addProperty(p);
75}
76
77// Utility functions for loading values with defaults
78// Single value properties only support int, double, string and bool
79template <typename Type>
80Type GetNeXusValue(const NeXus::NXEntry &entry, const std::string &path, const Type &defval, int32_t index) {
81 try {
82 NeXus::NXDataSetTyped<Type> dataSet = entry.openNXDataSet<Type>(path);
83 dataSet.load();
84
85 return dataSet()[index];
86 } catch (std::runtime_error &) {
87 return defval;
88 }
89}
90
91// string and double are special cases
92template <>
93double GetNeXusValue<double>(const NeXus::NXEntry &entry, const std::string &path, const double &defval,
94 int32_t index) {
95 try {
96 NeXus::NXDataSetTyped<float> dataSet = entry.openNXDataSet<float>(path);
97 dataSet.load();
98
99 return dataSet()[index];
100 } catch (std::runtime_error &) {
101 return defval;
102 }
103}
104
105template <>
106std::string GetNeXusValue<std::string>(const NeXus::NXEntry &entry, const std::string &path, const std::string &defval,
107 int32_t /*unused*/) {
108
109 try {
110 NeXus::NXChar dataSet = entry.openNXChar(path);
111 dataSet.load();
112
113 return std::string(dataSet(), dataSet.dim0());
114 } catch (std::runtime_error &) {
115 return defval;
116 }
117}
118
119template <typename T>
120void MapNeXusToProperty(const NeXus::NXEntry &entry, const std::string &path, const T &defval,
121 API::LogManager &logManager, const std::string &name, const T &factor, int32_t index) {
122
123 T value = GetNeXusValue<T>(entry, path, defval, index);
124 logManager.addProperty<T>(name, value * factor);
125}
126
127// string is a special case
128template <>
129void MapNeXusToProperty<std::string>(const NeXus::NXEntry &entry, const std::string &path, const std::string &defval,
130 API::LogManager &logManager, const std::string &name,
131 const std::string & /*unused*/, int32_t index) {
132
133 std::string value = GetNeXusValue<std::string>(entry, path, defval, index);
134 logManager.addProperty<std::string>(name, value);
135}
136
137template <typename T>
138void MapNeXusToSeries(const NeXus::NXEntry &entry, const std::string &path, const T &defval,
139 API::LogManager &logManager, const std::string &time, const std::string &name, const T &factor,
140 int32_t index) {
141
142 auto value = GetNeXusValue<T>(entry, path, defval, index);
143 AddSinglePointTimeSeriesProperty<T>(logManager, time, name, value * factor);
144}
145
146// map the comma separated range of indexes to the vector via a lambda function
147// throws an exception if it is outside the vector range
148//
149template <typename T, typename F> void mapRangeToIndex(const std::string &line, std::vector<T> &result, const F &fn) {
150
151 std::stringstream ss(line);
152 std::string item;
153 size_t index = 0;
154 while (std::getline(ss, item, ',')) {
155 auto const k = item.find('-');
156
157 size_t p0, p1;
158 if (k != std::string::npos) {
159 p0 = boost::lexical_cast<size_t>(item.substr(0, k));
160 p1 = boost::lexical_cast<size_t>(item.substr(k + 1, item.size() - k - 1));
161 } else {
162 p0 = boost::lexical_cast<size_t>(item);
163 p1 = p0;
164 }
165
166 if (p1 < result.size() && p0 <= p1) {
167 while (p0 <= p1) {
168 result[p0++] = fn(index);
169 index++;
170 }
171 } else if (p0 < result.size() && p1 < p0) {
172 do {
173 result[p0] = fn(index);
174 index++;
175 } while (p1 < p0--);
176 } else
177 throw std::invalid_argument("invalid range specification");
178 }
179}
180
181// Simple reader that is compatible with the ASNTO event file loader
182class FileLoader {
183 std::ifstream _ifs;
184 size_t _size;
185
186public:
187 explicit FileLoader(const char *filename) : _ifs(filename, std::ios::binary | std::ios::in) {
188 if (!_ifs.is_open() || _ifs.fail())
189 throw std::runtime_error("unable to open file");
190
191 _ifs.seekg(0, _ifs.end);
192 _size = _ifs.tellg();
193 _ifs.seekg(0, _ifs.beg);
194 }
195
196 bool read(char *s, std::streamsize n) { return static_cast<bool>(_ifs.read(s, n)); }
197
198 size_t size() const { return _size; }
199
200 size_t position() { return _ifs.tellg(); }
201
202 size_t selected_position() { return _ifs.tellg(); }
203};
204
205} // anonymous namespace
206
207namespace PLN {
208
209//
210// In the future the ANSTO helper and event file loader will be generalized to
211// handle the instruments consistently.
212
213// Simple 1D histogram class
215 std::vector<size_t> m_hist;
216 double m_M;
217 double m_B;
218 size_t m_peak;
219 size_t m_count;
220
221public:
222 SimpleHist(size_t N, double minVal, double maxVal) : m_hist(N, 0) {
223 m_M = (static_cast<double>(N) / (maxVal - minVal));
224 m_B = -m_M * minVal;
225 m_peak = 0;
226 m_count = 0;
227 }
228
229 inline double ival(double val) const { return m_M * val + m_B; }
230
231 inline double xval(double ix) const { return (ix - m_B) / m_M; }
232
233 inline void add(double val) {
234 auto ix = static_cast<size_t>(std::floor(ival(val)));
235 if (ix < m_hist.size()) {
236 m_hist[ix]++;
237 m_count++;
238 if (m_hist[ix] > m_peak)
239 m_peak = m_hist[ix];
240 }
241 }
242
243 const std::vector<size_t> &histogram() const { return m_hist; }
244
245 inline size_t peak() const { return m_peak; }
246 inline size_t count() const { return m_count; }
247};
248
250protected:
251 // fields
252 const std::vector<bool> &m_roi;
253 const std::vector<size_t> &m_mapIndex;
254 const double m_framePeriod;
255 const double m_gatePeriod;
256
257 // number of frames
258 size_t m_frames;
260
261 // number of events
265
266 // time boundaries
267 const TimeLimits m_timeBoundary; // seconds
268
269 virtual void addEventImpl(size_t id, size_t x, size_t y, double tof) = 0;
270
271public:
272 EventProcessor(const std::vector<bool> &roi, const std::vector<size_t> &mapIndex, const double framePeriod,
273 const double gatePeriod, const TimeLimits &timeBoundary, size_t maxEvents)
274 : m_roi(roi), m_mapIndex(mapIndex), m_framePeriod(framePeriod), m_gatePeriod(gatePeriod), m_frames(0),
276 m_timeBoundary(timeBoundary) {}
277
278 void newFrame() {
280 m_frames++;
281 if (validFrame())
283 }
284 }
285
286 inline bool validFrame() const {
287 double frameTime = m_framePeriod * static_cast<double>(m_frames) * 1.0e-6;
288 return (frameTime >= m_timeBoundary.first && frameTime <= m_timeBoundary.second);
289 }
290
291 double duration() const {
292 // length test in seconds
293 return m_framePeriod * static_cast<double>(m_frames) * 1.0e-6;
294 }
295
296 inline int64_t frameStart() const {
297 // returns time in nanoseconds from start of test
298 auto start = m_framePeriod * static_cast<double>(m_frames);
299 return static_cast<int64_t>(start * 1.0e3);
300 }
301
302 size_t numFrames() const { return m_framesValid; }
303
305
306 size_t processedEvents() const { return m_processedEvents; }
307
308 void addEvent(size_t x, size_t p, double tof, double /*taux*/) {
309
310 // check if in time boundaries
311 if (!validFrame())
312 return;
313
314 // group pixels
315 auto y = static_cast<size_t>(p / HISTO_BINS_Y_DENUMERATOR);
316
317 // determine detector id and check limits
318 if (x >= HISTO_BINS_X || y >= PIXELS_PER_TUBE)
319 return;
320
321 // map the raw detector index to the physical model
322 size_t xid = m_mapIndex[x];
323
324 size_t id = xid < DETECTOR_TUBES ? PIXELS_PER_TUBE * xid + y : DETECTOR_SPECTRA + xid;
325 if (id >= m_roi.size())
326 return;
327
328 // check if neutron is in region of interest
329 if (!m_roi[id])
330 return;
331
332 // finally pass to specific handler
334 // take the modulus of the tof time to account for the
335 // longer background chopper rate
336 double mtof = tof < 0.0 ? fmod(tof + m_gatePeriod, m_gatePeriod) : fmod(tof, m_gatePeriod);
337
338 addEventImpl(id, xid, y, mtof);
340 } else {
342 }
343 }
344};
345
346// The class determines the number of counts linked to the detectors and the
347// tof correction.
349protected:
350 // fields
351 std::vector<size_t> &m_eventCounts;
352 double m_L1;
353 double m_V0;
354 const std::vector<double> &m_L2;
356
357 void addEventImpl(size_t id, size_t /*x*/, size_t /*y*/, double tobs) override {
358 m_eventCounts[id]++;
359 // the maximum occurs at the elastic peak
360 double deltaT = 1.0e6 * (m_L1 + m_L2[id]) / m_V0 - tobs;
361 m_histogram.add(deltaT);
362 }
363
364public:
365 // construction
366 EventCounter(const std::vector<bool> &roi, const std::vector<size_t> &mapIndex, const double framePeriod,
367 const double gatePeriod, const TimeLimits &timeBoundary, std::vector<size_t> &eventCounts,
368 const double L1, const double V0, const std::vector<double> &vecL2, size_t maxEvents)
369 : EventProcessor(roi, mapIndex, framePeriod, gatePeriod, timeBoundary, maxEvents), m_eventCounts(eventCounts),
370 m_L1(L1), m_V0(V0), m_L2(vecL2), m_histogram(5000, -2500.0, 2500.0) {}
371
372 // clips the histogram above 25% and takes the mean of the values
373 double tofCorrection() {
374
375 // determine the points above the 25% threshold
376 auto minLevel = static_cast<size_t>(m_histogram.peak() / 4);
377 auto hvec = m_histogram.histogram();
378 double sum = 0.0;
379 size_t count = 0;
380 for (size_t i = 0; i < hvec.size(); i++) {
381 if (hvec[i] >= minLevel) {
382 auto ix = static_cast<double>(i);
383 sum += static_cast<double>(hvec[i]) * m_histogram.xval(ix + 0.5);
384 count += hvec[i];
385 }
386 }
387
388 return (count > 0 ? sum / static_cast<double>(count) : 0.0);
389 }
390};
391
393protected:
394 // fields
395 std::vector<EventVector_pt> &m_eventVectors;
396 double m_tofMin;
397 double m_tofMax;
398 int64_t m_startTime;
401
402 void addEventImpl(size_t id, size_t /*x*/, size_t /*y*/, double tobs) override {
403
404 // get the absolute time for the start of the frame
405 auto const offset = m_startTime + frameStart();
406
407 // adjust the the tof to account for the correction and allocate events
408 // that occur before the sample time as slow events from the previous pulse
409 double tof = tobs + m_tofCorrection - m_sampleTime;
410 if (tof < 0.0)
411 tof = fmod(tof + m_gatePeriod, m_gatePeriod);
412 tof += m_sampleTime;
413 if (m_tofMin > tof)
414 m_tofMin = tof;
415 if (m_tofMax < tof)
416 m_tofMax = tof;
417
418 auto ev = Types::Event::TofEvent(tof, Types::Core::DateAndTime(offset));
419 m_eventVectors[id]->emplace_back(ev);
420 }
421
422public:
423 EventAssigner(const std::vector<bool> &roi, const std::vector<size_t> &mapIndex, const double framePeriod,
424 const double gatePeriod, const TimeLimits &timeBoundary, std::vector<EventVector_pt> &eventVectors,
425 int64_t startTime, double tofCorrection, double sampleTime, size_t maxEvents)
426 : EventProcessor(roi, mapIndex, framePeriod, gatePeriod, timeBoundary, maxEvents), m_eventVectors(eventVectors),
427 m_tofMin(std::numeric_limits<double>::max()), m_tofMax(std::numeric_limits<double>::min()),
428 m_startTime(startTime), m_tofCorrection(tofCorrection), m_sampleTime(sampleTime) {}
429
430 double tofMin() const { return m_tofMin <= m_tofMax ? m_tofMin : 0.0; }
431 double tofMax() const { return m_tofMin <= m_tofMax ? m_tofMax : 0.0; }
432};
433
434template <typename EP>
435void loadEvents(API::Progress &prog, const char *progMsg, const std::string &eventFile, EP &eventProcessor) {
436
437 using namespace ANSTO;
438
439 prog.doReport(progMsg);
440
441 FileLoader loader(eventFile.c_str());
442
443 // for progress notifications
444 ANSTO::ProgressTracker progTracker(prog, progMsg, loader.size(), Progress_LoadBinFile);
445
446 ReadEventFile(loader, eventProcessor, progTracker, 100, false);
447}
448} // namespace PLN
449
453
454 // Specify file extensions which can be associated with a specific file.
455 std::vector<std::string> exts;
456
457 // Declare the Filename algorithm property. Mandatory. Sets the path to the
458 // file to load.
459 exts.clear();
460 exts.emplace_back(".hdf");
461 declareProperty(std::make_unique<API::FileProperty>(FilenameStr, "", API::FileProperty::Load, exts),
462 "The input filename of the stored data");
463
464 declareProperty(PathToBinaryStr, "./", std::make_shared<Kernel::MandatoryValidator<std::string>>(),
465 "Relative or absolute path to the compressed binary\n"
466 "event file linked to the HDF file, eg /storage/data/");
467
468 // mask
469 exts.clear();
470 exts.emplace_back(".xml");
471 declareProperty(std::make_unique<API::FileProperty>(MaskStr, "", API::FileProperty::OptionalLoad, exts),
472 "The input filename of the mask data");
473
474 declareProperty(SelectDetectorTubesStr, "",
475 "Comma separated range of detectors tubes to be loaded,\n"
476 " eg 16,19-45,47");
477
479 std::make_unique<API::WorkspaceProperty<API::IEventWorkspace>>("OutputWorkspace", "", Kernel::Direction::Output));
480
481 declareProperty(SelectDatasetStr, 0, "Select the index for the dataset to be loaded.");
482
483 declareProperty(TOFBiasStr, 0.0, "Time of flight correction in micro-sec.");
484
485 declareProperty(CalibrateTOFStr, false, "Calibrate the TOF correction from the elastic pulse.");
486
487 declareProperty(LambdaOnTwoStr, false, "Instrument is operating in Lambda on Two mode.");
488
490 "Only include events after the provided start time, in "
491 "seconds (relative to the start of the run).");
492
494 "Only include events before the provided stop time, in "
495 "seconds (relative to the start of the run).");
496
497 std::string grpOptional = "Filters";
500}
501
503void LoadPLN::createWorkspace(const std::string &title) {
504
505 // Create the workspace
506 m_localWorkspace = std::make_shared<DataObjects::EventWorkspace>();
507 m_localWorkspace->initialize(HISTOGRAMS, 2, 1);
508
509 // set the units
510 m_localWorkspace->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("TOF");
511 m_localWorkspace->setYUnit("Counts");
512
513 // set title
514 m_localWorkspace->setTitle(title);
515}
516
524
525void LoadPLN::exec(const std::string &hdfFile, const std::string &eventFile) {
526
527 namespace fs = boost::filesystem;
528
529 // Create workspace
530 // ----------------
531 fs::path p = hdfFile;
532 for (; !p.extension().empty();)
533 p = p.stem();
534 std::string title = p.generic_string();
535 createWorkspace(title);
536 API::LogManager &logManager = m_localWorkspace->mutableRun();
537 API::Progress prog(this, 0.0, 1.0, Progress_Total);
538
539 // Load instrument and workspace properties
540 logManager.addProperty(SelectDatasetStr, m_datasetIndex);
541 loadParameters(hdfFile, logManager);
542 prog.doReport("creating instrument");
544
545 // Get the region of interest and filters and save to log
546 std::string const maskfile = getPropertyValue(MaskStr);
547 std::string const seltubes = getPropertyValue(SelectDetectorTubesStr);
548 logManager.addProperty(SelectDetectorTubesStr, seltubes);
549 logManager.addProperty(MaskStr, maskfile);
550
551 std::vector<bool> roi = createRoiVector(seltubes, maskfile);
552 double timeMaxBoundary = getProperty(FilterByTimeStopStr);
553 if (isEmpty(timeMaxBoundary))
554 timeMaxBoundary = std::numeric_limits<double>::infinity();
555 TimeLimits timeBoundary(getProperty(FilterByTimeStartStr), timeMaxBoundary);
556
557 // get the detector map from raw input to a physical detector
558 auto instr = m_localWorkspace->getInstrument();
559 std::string dmapStr = instr->getParameterAsString("DetectorMap");
560 std::vector<size_t> detMapIndex = std::vector<size_t>(HISTO_BINS_X, 0);
561 mapRangeToIndex(dmapStr, detMapIndex, [](size_t n) { return n; });
562
563 // Load the events file. First count the number of events to reserve
564 // memory and then assign the events to the detectors
565 size_t numberHistograms = m_localWorkspace->getNumberHistograms();
566 std::vector<EventVector_pt> eventVectors(numberHistograms, nullptr);
567 std::vector<size_t> eventCounts(numberHistograms, 0);
568
569 double masterRpm = fabs(logManager.getTimeSeriesProperty<double>("FermiChopperFreq")->firstValue());
570 double slaveRpm = fabs(logManager.getTimeSeriesProperty<double>("OverlapChopperFreq")->firstValue());
571 double framePeriod = 1.0e6 / masterRpm;
572
573 // if fermi chopper freq equals the overlap freq then the gate period is
574 // half the frame period
575 double gatePeriod = (std::round(masterRpm / slaveRpm) == 1.0 ? 0.5 * framePeriod : framePeriod);
576 AddSinglePointTimeSeriesProperty<double>(logManager, m_startRun, "GatePeriod", gatePeriod);
577
578 // count total events per pixel and reserve necessary memory
579 size_t hdfCounts = static_cast<size_t>(logManager.getTimeSeriesProperty<int>("TotalCounts")->firstValue());
581 double sourceSample = fabs(instr->getSource()->getPos().Z());
582 double wavelength = logManager.getTimeSeriesProperty<double>("Wavelength")->firstValue();
583 double velocity = PhysicalConstants::h / (PhysicalConstants::NeutronMass * wavelength * 1e-10);
584 double sampleTime = 1.0e6 * sourceSample / velocity;
585 PLN::EventCounter eventCounter(roi, detMapIndex, framePeriod, gatePeriod, timeBoundary, eventCounts, sourceSample,
586 velocity, m_detectorL2, hdfCounts);
587 PLN::loadEvents(prog, "loading neutron counts", eventFile, eventCounter);
588 ANSTO::ProgressTracker progTracker(prog, "creating neutron event lists", numberHistograms, Progress_ReserveMemory);
589 prepareEventStorage(progTracker, eventCounts, eventVectors);
590
591 // log a message if the number of events in the event file does not match
592 // the total counts in the hdf
593 if (hdfCounts != eventCounter.availableEvents()) {
594 g_log.error("HDF and event counts differ: " + std::to_string(hdfCounts) + ", " +
595 std::to_string(eventCounter.availableEvents()));
596 }
597
598 // now perform the actual event collection and TOF convert if necessary
599 // if a phase calibration is required then load it as raw doppler time
600 // perform the calibration and then convert to TOF
601 Types::Core::DateAndTime startTime(m_startRun);
602 auto const start_nanosec = startTime.totalNanoseconds();
603 bool const calibrateTOF = getProperty(CalibrateTOFStr);
604 double tofCorrection = getProperty(TOFBiasStr);
605 if (calibrateTOF) {
606 tofCorrection = eventCounter.tofCorrection();
607 }
608 logManager.addProperty("CalibrateTOF", (calibrateTOF ? 1 : 0));
609 AddSinglePointTimeSeriesProperty<double>(logManager, m_startRun, "TOFCorrection", tofCorrection);
610 PLN::EventAssigner eventAssigner(roi, detMapIndex, framePeriod, gatePeriod, timeBoundary, eventVectors, start_nanosec,
611 tofCorrection, sampleTime, hdfCounts);
612 PLN::loadEvents(prog, "loading neutron events (TOF)", eventFile, eventAssigner);
613
614 // perform a calibration and then TOF conversion if necessary
615 // and update the tof limits
616 auto minTOF = eventAssigner.tofMin();
617 auto maxTOF = eventAssigner.tofMax();
618
619 // just to make sure the bins hold it all and setup the detector masks
620 m_localWorkspace->setAllX(HistogramData::BinEdges{std::max(0.0, floor(minTOF)), maxTOF + 1});
622
623 // set log values
624 auto frame_count = static_cast<int>(eventCounter.numFrames());
625 AddSinglePointTimeSeriesProperty<int>(logManager, m_startRun, "frame_count", frame_count);
626
627 std::string filename = getPropertyValue(FilenameStr);
628 logManager.addProperty("filename", filename);
629
630 Types::Core::time_duration duration =
631 boost::posix_time::microseconds(static_cast<boost::int64_t>(eventCounter.duration() * 1.0e6));
632 Types::Core::DateAndTime endTime(startTime + duration);
633 logManager.addProperty("start_time", startTime.toISO8601String());
634 logManager.addProperty("end_time", endTime.toISO8601String());
635 logManager.addProperty<double>("dur", eventCounter.duration());
636
637 // Finally add the time-series evironment parameters explicitly
638 loadEnvironParameters(hdfFile, logManager);
639
640 setProperty("OutputWorkspace", m_localWorkspace);
641}
642
645
646 m_detectorL2.resize(HISTOGRAMS, 0.0);
647 const auto &detectorInfo = m_localWorkspace->detectorInfo();
648 auto detIDs = detectorInfo.detectorIDs();
649 for (const auto detID : detIDs) {
650 auto ix = detectorInfo.indexOf(detID);
651 double l2 = detectorInfo.l2(ix);
652 m_detectorL2[detID] = l2;
653 }
654}
655
657void LoadPLN::setupDetectorMasks(const std::vector<bool> &roi) {
658
659 // count total number of masked bins
660 size_t maskedBins = 0;
661 for (size_t i = 0; i != roi.size(); i++)
662 if (!roi[i])
663 maskedBins++;
664
665 if (maskedBins > 0) {
666 // create list of masked bins
667 std::vector<size_t> maskIndexList(maskedBins);
668 size_t maskIndex = 0;
669
670 for (size_t i = 0; i != roi.size(); i++)
671 if (!roi[i])
672 maskIndexList[maskIndex++] = i;
673
674 auto maskingAlg = createChildAlgorithm("MaskDetectors");
675 maskingAlg->setProperty("Workspace", m_localWorkspace);
676 maskingAlg->setProperty("WorkspaceIndexList", maskIndexList);
677 maskingAlg->executeAsChildAlg();
678 }
679}
680
683void LoadPLN::prepareEventStorage(ANSTO::ProgressTracker &progTracker, std::vector<size_t> &eventCounts,
684 std::vector<EventVector_pt> &eventVectors) {
685
686 size_t numberHistograms = eventCounts.size();
687 for (size_t i = 0; i != numberHistograms; ++i) {
688 DataObjects::EventList &eventList = m_localWorkspace->getSpectrum(i);
689
691 eventList.reserve(eventCounts[i]);
692
693 eventList.setDetectorID(static_cast<detid_t>(i));
694 eventList.setSpectrumNo(static_cast<detid_t>(i));
695
696 DataObjects::getEventsFrom(eventList, eventVectors[i]);
697
698 progTracker.update(i);
699 }
700 progTracker.complete();
701}
702
705std::vector<bool> LoadPLN::createRoiVector(const std::string &selected, const std::string &maskfile) {
706
707 std::vector<bool> result(HISTOGRAMS, true);
708
709 // turn off pixels linked to missing tubes
710 if (!selected.empty()) {
711 std::vector<bool> tubes(HISTO_BINS_X, false);
712 mapRangeToIndex(selected, tubes, [](size_t) { return true; });
713 for (size_t i = 0; i < DETECTOR_TUBES; i++) {
714 if (tubes[i] == false) {
715 for (size_t j = 0; j < PIXELS_PER_TUBE; j++) {
716 result[i * PIXELS_PER_TUBE + j] = false;
717 }
718 }
719 }
720 for (size_t i = 0; i < MONITORS; i++) {
721 result[DETECTOR_SPECTRA + i] = tubes[DETECTOR_TUBES + i];
722 }
723 }
724
725 if (maskfile.length() == 0)
726 return result;
727
728 std::ifstream input(maskfile.c_str());
729 if (!input.good())
730 throw std::invalid_argument("invalid mask file");
731
732 std::string line;
733 while (std::getline(input, line)) {
734 auto i0 = line.find("<detids>");
735 auto iN = line.find("</detids>");
736
737 if ((i0 != std::string::npos) && (iN != std::string::npos) && (i0 < iN)) {
738 line = line.substr(i0 + 8, iN - i0 - 8); // 8 = len("<detids>")
739 mapRangeToIndex(line, result, [](size_t) { return false; });
740 }
741 }
742
743 return result;
744}
745
747void LoadPLN::loadParameters(const std::string &hdfFile, API::LogManager &logm) {
748
749 NeXus::NXRoot root(hdfFile);
750 NeXus::NXEntry entry = root.openFirstEntry();
751
752 MapNeXusToProperty<std::string>(entry, "sample/name", "unknown", logm, "SampleName", "", 0);
753 MapNeXusToProperty<std::string>(entry, "sample/description", "unknown", logm, "SampleDescription", "", 0);
754
755 // if dataset index > 0 need to add an offset to the start time
756 Types::Core::DateAndTime startTime(GetNeXusValue<std::string>(entry, "start_time", "2000-01-01T00:00:00", 0));
757 if (m_datasetIndex > 0) {
758 auto baseTime = GetNeXusValue<int32_t>(entry, "instrument/detector/start_time", 0, 0);
759 auto nthTime = GetNeXusValue<int32_t>(entry, "instrument/detector/start_time", 0, m_datasetIndex);
760
761 Types::Core::time_duration duration =
762 boost::posix_time::microseconds((static_cast<int64_t>(nthTime) - static_cast<int64_t>(baseTime)) * 1'000'000);
763 Types::Core::DateAndTime startDataset(startTime + duration);
764 m_startRun = startDataset.toISO8601String();
765 } else {
766 m_startRun = startTime.toISO8601String();
767 }
768
769 // Add support for instrument running in lambda on two mode.
770 // Added as UI option as the available instrument parameters
771 // cannot be reliably interpreted to predict the mode (as
772 // advised by the instrument scientist).
773 bool const lambdaOnTwoMode = getProperty(LambdaOnTwoStr);
774 double lambdaFactor = (lambdaOnTwoMode ? 0.5 : 1.0);
775 logm.addProperty("LambdaOnTwoMode", (lambdaOnTwoMode ? 1 : 0));
776
777 MapNeXusToSeries<double>(entry, "instrument/fermi_chopper/mchs", 0.0, logm, m_startRun, "FermiChopperFreq", 1.0 / 60,
779 MapNeXusToSeries<double>(entry, "instrument/fermi_chopper/schs", 0.0, logm, m_startRun, "OverlapChopperFreq",
780 1.0 / 60, m_datasetIndex);
781 MapNeXusToSeries<double>(entry, "instrument/crystal/wavelength", 0.0, logm, m_startRun, "Wavelength", lambdaFactor,
783 MapNeXusToSeries<double>(entry, "instrument/detector/stth", 0.0, logm, m_startRun, "DetectorTankAngle", 1.0,
785 MapNeXusToSeries<int32_t>(entry, "monitor/bm1_counts", 0, logm, m_startRun, "MonitorCounts", 1, m_datasetIndex);
786 MapNeXusToSeries<int32_t>(entry, "data/total_counts", 0, logm, m_startRun, "TotalCounts", 1, m_datasetIndex);
787 MapNeXusToSeries<double>(entry, "data/tofw", 5.0, logm, m_startRun, "ChannelWidth", 1, m_datasetIndex);
788 MapNeXusToSeries<double>(entry, "sample/mscor", 0.0, logm, m_startRun, "SampleRotation", 1, m_datasetIndex);
789}
790
793void LoadPLN::loadEnvironParameters(const std::string &hdfFile, API::LogManager &logm) {
794
795 NeXus::NXRoot root(hdfFile);
796 NeXus::NXEntry entry = root.openFirstEntry();
797 auto time_str = logm.getPropertyValueAsType<std::string>("end_time");
798
799 // load the environment variables for the dataset loaded
800 std::vector<std::string> tags = {"P01PS05", "P01PSP05", "T01S00", "T01S06", "T01S07", "T01S08", "T01SP00", "T01SP06",
801 "T01SP07", "T01SP08", "T2S1", "T2S2", "T2S3", "T2S4", "T2SP1", "T2SP2"};
802
803 for (const auto &tag : tags) {
804 MapNeXusToSeries<double>(entry, "data/" + tag, 0.0, logm, time_str, "env_" + tag, 1.0, m_datasetIndex);
805 }
806}
807
810
811 // loads the IDF and parameter file
812 auto loadInstrumentAlg = createChildAlgorithm("LoadInstrument");
813 loadInstrumentAlg->setProperty("Workspace", m_localWorkspace);
814 loadInstrumentAlg->setPropertyValue("InstrumentName", "PELICAN");
815 loadInstrumentAlg->setProperty("RewriteSpectraMap", Mantid::Kernel::OptionalBool(false));
816 loadInstrumentAlg->executeAsChildAlg();
817}
818
820int LoadPLN::version() const { return 1; }
821
823const std::vector<std::string> LoadPLN::seeAlso() const { return {"Load", "LoadEMU"}; }
825const std::string LoadPLN::category() const { return "DataHandling\\ANSTO"; }
826
828const std::string LoadPLN::name() const { return "LoadPLN"; }
829
831const std::string LoadPLN::summary() const { return "Loads a PLN Hdf and linked event file into a workspace."; }
832
836 if (descriptor.extension() != ".hdf")
837 return 0;
838
839 if (descriptor.pathExists("/entry1/site_name") && descriptor.pathExists("/entry1/instrument/fermi_chopper") &&
840 descriptor.pathExists("/entry1/instrument/aperture/sh1") &&
841 descriptor.pathExists("/entry1/instrument/ag1010/MEAS/Temperature") &&
842 descriptor.pathExists("/entry1/instrument/detector/daq_dirname") &&
843 descriptor.pathExists("/entry1/instrument/detector/dataset_number") &&
844 descriptor.pathExists("/entry1/data/hmm") && descriptor.pathExists("/entry1/data/time_of_flight") &&
845 descriptor.pathExists("/entry1/data/total_counts")) {
846 return 80;
847 } else {
848 return 0;
849 }
850}
851
854// exec() function that works with the two files.
856
857 namespace fs = boost::filesystem;
858
859 // Open the hdf file and find the dirname and dataset number
860 std::string hdfFile = getPropertyValue(FilenameStr);
861 std::string evtPath = getPropertyValue(PathToBinaryStr);
862 if (evtPath.empty())
863 evtPath = "./";
864
865 // if relative ./ or ../ then append to the directory for the hdf file
866 if (evtPath.rfind("./") == 0 || evtPath.rfind("../") == 0) {
867 fs::path hp = hdfFile;
868 evtPath = fs::canonical(evtPath, hp.parent_path()).generic_string();
869 }
870
871 // dataset index to be loaded
872 m_datasetIndex = getProperty(SelectDatasetStr);
873
874 // if path provided build the file path from the directory name and dataset
875 // number from the hdf file, however if this is not a valid path then try
876 // the basename with a '.bin' extension
877 if (fs::is_directory(evtPath)) {
878 NeXus::NXRoot root(hdfFile);
879 NeXus::NXEntry entry = root.openFirstEntry();
880 auto eventDir = GetNeXusValue<std::string>(entry, "instrument/detector/daq_dirname", "./", 0);
881 auto dataset = GetNeXusValue<int32_t>(entry, "instrument/detector/dataset_number", 0, m_datasetIndex);
882 if (dataset < 0) {
883 std::string message("Negative dataset index recorded in HDF, reset to zero!");
884 g_log.error(message);
885 dataset = 0;
886 }
887
888 // build the path to the event file using the standard storage convention at ansto:
889 // 'relpath/[daq_dirname]/DATASET_[n]/EOS.bin'
890 // but if the file is missing, try relpath/{source}.bin
891 std::stringstream buffer;
892 buffer << eventDir.c_str() << "/DATASET_" << dataset << "/EOS.bin";
893 fs::path filePath = evtPath;
894 filePath /= buffer.str();
895 filePath = fs::absolute(filePath);
896 std::string nomPath = filePath.generic_string();
897 if (fs::is_regular_file(nomPath)) {
898 evtPath = nomPath;
899 } else {
900 fs::path hp = hdfFile;
901 buffer.str("");
902 buffer.clear();
903 buffer << hp.stem().generic_string().c_str() << ".bin";
904 fs::path path = evtPath;
905 path /= buffer.str();
906 path = fs::absolute(path);
907 evtPath = path.generic_string();
908 }
909 }
910
911 // finally check that the event file exists
912 if (!fs::is_regular_file(evtPath)) {
913 std::string msg = "Check path, cannot open binary event file: " + evtPath;
914 throw std::runtime_error(msg);
915 }
916
917 exec(hdfFile, evtPath);
918}
919
920// register the algorithms into the AlgorithmFactory
922
923} // 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
size_t _size
Definition: LoadEMU.cpp:338
std::ifstream _ifs
Definition: LoadEMU.cpp:337
#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...
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
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 std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
static bool isEmpty(const NumT toCheck)
checks that the value was not set by users, uses the value in empty double/int.
@ 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
LoadPLN : Loads an ANSTO PLN Hdf and linked event file into a workspace.
Definition: LoadPLN.h:59
std::vector< double > m_detectorL2
Definition: LoadPLN.h:102
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
Definition: LoadPLN.cpp:831
const std::vector< std::string > seeAlso() const override
Similar algorithms.
Definition: LoadPLN.cpp:823
void prepareEventStorage(ANSTO::ProgressTracker &prog, 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: LoadPLN.cpp:683
int version() const override
Algorithm's version for identification.
Definition: LoadPLN.cpp:820
void setupDetectorMasks(const std::vector< bool > &roi)
Set up the detector masks to the region of interest roi.
Definition: LoadPLN.cpp:657
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: LoadPLN.cpp:793
DataObjects::EventWorkspace_sptr m_localWorkspace
Definition: LoadPLN.h:99
void init() override
Initialise the algorithm and declare the properties for the nexus descriptor.
Definition: LoadPLN.cpp:452
void createWorkspace(const std::string &title)
Creates an event workspace and sets the title.
Definition: LoadPLN.cpp:503
void loadParameters(const std::string &hdfFile, API::LogManager &logm)
Load parameters from input hdfFile and save to the log manager, logm.
Definition: LoadPLN.cpp:747
void loadDetectorL2Values()
Recovers the L2 neutronic distance for each detector.
Definition: LoadPLN.cpp:644
void loadInstrument()
Load the instrument definition.
Definition: LoadPLN.cpp:809
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: LoadPLN.cpp:705
const std::string name() const override
Algorithms name for identification.
Definition: LoadPLN.cpp:828
void exec() override
Execute the algorithm.
Definition: LoadPLN.cpp:855
int confidence(Kernel::NexusDescriptor &descriptor) const override
Return the confidence as an integer value that this algorithm can load the file descriptor.
Definition: LoadPLN.cpp:835
const std::string category() const override
Algorithm's category for identification.
Definition: LoadPLN.cpp:825
EventAssigner(const std::vector< bool > &roi, const std::vector< size_t > &mapIndex, const double framePeriod, const double gatePeriod, const TimeLimits &timeBoundary, std::vector< EventVector_pt > &eventVectors, int64_t startTime, double tofCorrection, double sampleTime, size_t maxEvents)
Definition: LoadPLN.cpp:423
std::vector< EventVector_pt > & m_eventVectors
Definition: LoadPLN.cpp:395
void addEventImpl(size_t id, size_t, size_t, double tobs) override
Definition: LoadPLN.cpp:402
std::vector< size_t > & m_eventCounts
Definition: LoadPLN.cpp:351
EventCounter(const std::vector< bool > &roi, const std::vector< size_t > &mapIndex, const double framePeriod, const double gatePeriod, const TimeLimits &timeBoundary, std::vector< size_t > &eventCounts, const double L1, const double V0, const std::vector< double > &vecL2, size_t maxEvents)
Definition: LoadPLN.cpp:366
void addEventImpl(size_t id, size_t, size_t, double tobs) override
Definition: LoadPLN.cpp:357
const std::vector< double > & m_L2
Definition: LoadPLN.cpp:354
const std::vector< bool > & m_roi
Definition: LoadPLN.cpp:252
void addEvent(size_t x, size_t p, double tof, double)
Definition: LoadPLN.cpp:308
const std::vector< size_t > & m_mapIndex
Definition: LoadPLN.cpp:253
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 double framePeriod, const double gatePeriod, const TimeLimits &timeBoundary, size_t maxEvents)
Definition: LoadPLN.cpp:272
double ival(double val) const
Definition: LoadPLN.cpp:229
std::vector< size_t > m_hist
Definition: LoadPLN.cpp:215
double xval(double ix) const
Definition: LoadPLN.cpp:231
SimpleHist(size_t N, double minVal, double maxVal)
Definition: LoadPLN.cpp:222
const std::vector< size_t > & histogram() const
Definition: LoadPLN.cpp:243
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
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertyGroup(const std::string &name, const std::string &group)
Set the group for a given property.
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
Validator to check that a property is not left empty.
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...
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.
Kernel::Logger g_log("ExperimentInfo")
static logger object
void loadEvents(API::Progress &prog, const char *progMsg, const std::string &eventFile, EP &eventProcessor)
Definition: LoadPLN.cpp:435
static char const *const FilenameStr
Definition: LoadBBY.cpp:54
static const size_t HISTO_BINS_X
Definition: LoadBBY.cpp:47
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
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)
NXDataSetTyped< char > NXChar
The char dataset type.
Definition: NexusClasses.h:469
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double h
Planck constant in J*s.
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.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Output
An output workspace.
Definition: Property.h:54