From 4fd3de6e4e551a36c8c088952fece559a101149e Mon Sep 17 00:00:00 2001
From: Geish Miladinovic <geish.miladinovic@ansto.gov.au>
Date: Thu, 6 Dec 2018 14:53:39 +1100
Subject: [PATCH] Merged emu workflow support (no ticket)
---
Framework/API/inc/MantidAPI/IEventList.h | 2 +
.../MantidDataHandling/LoadANSTOEventFile.h | 6 +-
.../inc/MantidDataHandling/LoadEMU.h | 43 +-
Framework/DataHandling/src/LoadEMU.cpp | 837 ++++++++++++++----
Framework/DataHandling/test/LoadEMUauTest.h | 7 +-
.../inc/MantidDataObjects/EventList.h | 5 +
Framework/DataObjects/src/EventList.cpp | 86 ++
Framework/DataObjects/test/EventListTest.h | 35 +
.../mantid/api/src/Exports/IEventList.cpp | 14 +-
.../api/src/Exports/IEventWorkspace.cpp | 6 +
10 files changed, 842 insertions(+), 199 deletions(-)
diff --git a/Framework/API/inc/MantidAPI/IEventList.h b/Framework/API/inc/MantidAPI/IEventList.h
index 5e699f6b510..24d9473b5f6 100644
--- a/Framework/API/inc/MantidAPI/IEventList.h
+++ b/Framework/API/inc/MantidAPI/IEventList.h
@@ -77,6 +77,8 @@ public:
virtual void addPulsetime(const double seconds) = 0;
/// Mask a given TOF range
virtual void maskTof(const double tofMin, const double tofMax) = 0;
+ /// Mask the events by the condition vector
+ virtual void maskCondition(const std::vector<bool> & mask) = 0;
/// Return the list of TOF values
virtual std::vector<double> getTofs() const = 0;
/// Return the list of TOF values
diff --git a/Framework/DataHandling/inc/MantidDataHandling/LoadANSTOEventFile.h b/Framework/DataHandling/inc/MantidDataHandling/LoadANSTOEventFile.h
index 549a68d77d8..006e80d7f75 100644
--- a/Framework/DataHandling/inc/MantidDataHandling/LoadANSTOEventFile.h
+++ b/Framework/DataHandling/inc/MantidDataHandling/LoadANSTOEventFile.h
@@ -91,11 +91,11 @@ void ReadEventFile(IReader &loader, IEventHandler &handler, IProgress &progress,
int32_t def_clock_scale, bool use_tx_chopper) {
// read file headers (base header then packed-format header)
EventFileHeader_Base hdr_base;
- if (!loader.read(reinterpret_cast<int8_t *>(&hdr_base), sizeof(hdr_base)))
+ if (!loader.read(reinterpret_cast<char *>(&hdr_base), sizeof(hdr_base)))
throw std::runtime_error("unable to load EventFileHeader-Base");
EventFileHeader_Packed hdr_packed;
- if (!loader.read(reinterpret_cast<int8_t *>(&hdr_packed), sizeof(hdr_packed)))
+ if (!loader.read(reinterpret_cast<char *>(&hdr_packed), sizeof(hdr_packed)))
throw std::runtime_error("unable to load EventFileHeader-Packed");
if (hdr_base.magic_number != EVENTFILEHEADER_BASE_MAGIC_NUMBER)
@@ -185,7 +185,7 @@ void ReadEventFile(IReader &loader, IEventHandler &handler, IProgress &progress,
// read next byte
uint8_t ch;
- if (!loader.read(reinterpret_cast<int8_t *>(&ch), 1))
+ if (!loader.read(reinterpret_cast<char *>(&ch), 1))
break;
int32_t nbits_ch_used = 0; // no bits used initially, 8 to go
diff --git a/Framework/DataHandling/inc/MantidDataHandling/LoadEMU.h b/Framework/DataHandling/inc/MantidDataHandling/LoadEMU.h
index 829ac8d3fd8..89b311cbc96 100644
--- a/Framework/DataHandling/inc/MantidDataHandling/LoadEMU.h
+++ b/Framework/DataHandling/inc/MantidDataHandling/LoadEMU.h
@@ -31,12 +31,15 @@ Required Properties:
Optional Properties:
<UL>
+<LI> PathToBinaryEventFile - Rel or abs path to event file linked to hdf</LI>
<LI> Mask - The input filename of the mask data</LI>
<LI> SelectDetectorTubes - Range of detector tubes to be loaded</LI>
+<LI> OverrideDopplerFrequency - Override the Doppler frequency (Hz)</LI>
<LI> OverrideDopplerPhase - Override the Doppler phase (degrees)</LI>
+<LI> CalibrateDopplerPhase - Calibrate the Doppler phase prior to TOF</LI>
+<LI> LoadAsRawDopplerTime - Save event time relative the Doppler</LI>
<LI> FilterByTimeStart - Only include events after the start time</LI>
<LI> FilterByTimeStop - Only include events before the stop time</LI>
-<LI> LoadAsRawDopplerTime - Save event time relative the Doppler</LI>
</UL>
@author Geish Miladinovic (ANSTO)
@@ -63,7 +66,8 @@ File change history is stored at: <https://github.com/mantidproject/mantid>.
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
-class DLLExport LoadEMU : public API::IFileLoader<Kernel::FileDescriptor> {
+template <typename FD>
+class DLLExport LoadEMU : public API::IFileLoader<FD> {
public:
// description
@@ -71,34 +75,36 @@ public:
const std::vector<std::string> seeAlso() const override {
return {"Load", "LoadQKK"};
}
- const std::string name() const override { return "LoadEMU"; }
+ const std::string name() const override;
const std::string category() const override { return "DataHandling\\ANSTO"; }
- const std::string summary() const override {
- return "Loads an EMU data file into a workspace.";
- }
+ const std::string summary() const override;
// returns a confidence value that this algorithm can load a specified file
- int confidence(Kernel::FileDescriptor &descriptor) const override;
+ int confidence(FD &descriptor) const override;
private:
// initialisation
void init() override;
+ void init(bool hdfLoader);
// execution
void exec() override;
+ void exec(const std::string &hdfFile, const std::string &eventFile);
// region of intreset
static std::vector<bool> createRoiVector(const std::string &seltubes,
const std::string &maskfile);
// load parameters from input file
- void loadParameters(ANSTO::Tar::File &tarFile, API::LogManager &logm);
+ void loadParameters(const std::string &hdfFile, API::LogManager &logm);
+ void loadEnvironParameters(const std::string &hdfFile,
+ API::LogManager &logm);
// load the instrument definition and instrument parameters
void loadInstrument();
// create workspace
- void createWorkspace(ANSTO::Tar::File &tarFile);
+ void createWorkspace(const std::string &title);
// dynamically update the neutronic position
void loadDetectorL2Values();
@@ -107,6 +113,12 @@ private:
// load and log the Doppler parameters
void loadDopplerParameters(API::LogManager &logm);
+ // calibrate doppler phase
+ void calibrateDopplerPhase(std::vector<size_t> &eventCounts,
+ std::vector<EventVector_pt> &eventVectors);
+ void dopplerTimeToTOF(std::vector<EventVector_pt> &eventVectors,
+ double& minTOF, double& maxTOF);
+
// prepare event storage
void prepareEventStorage(ANSTO::ProgressTracker &prog,
std::vector<size_t> &eventCounts,
@@ -116,21 +128,28 @@ private:
void setupDetectorMasks(std::vector<bool> &roi);
// binary file access
- template <class EventProcessor>
+ template <class Processor>
static void loadEvents(API::Progress &prog, const char *progMsg,
- ANSTO::Tar::File &tarFile,
- EventProcessor &eventProcessor);
+ const std::string &eventFile,
+ Processor &eventProcessor);
// shared member variables
DataObjects::EventWorkspace_sptr m_localWorkspace;
std::vector<double> m_detectorL2;
+ int32_t m_datasetIndex;
+ std::string m_startRun;
// Doppler characteristics
double m_dopplerAmpl;
double m_dopplerFreq;
double m_dopplerPhase;
+ int32_t m_dopplerRun;
+ bool m_calibrateDoppler;
};
+using LoadEMUHdf = LoadEMU<Kernel::NexusDescriptor>;
+using LoadEMUTar = LoadEMU<Kernel::FileDescriptor>;
+
} // namespace DataHandling
} // namespace Mantid
diff --git a/Framework/DataHandling/src/LoadEMU.cpp b/Framework/DataHandling/src/LoadEMU.cpp
index 3fe20129738..8086397332b 100644
--- a/Framework/DataHandling/src/LoadEMU.cpp
+++ b/Framework/DataHandling/src/LoadEMU.cpp
@@ -19,12 +19,18 @@
#include "MantidNexus/NexusClasses.h"
#include <boost/math/special_functions/round.hpp>
+#include <boost/math/tools/minima.hpp>
+#include <boost/filesystem.hpp>
#include <Poco/AutoPtr.h>
#include <Poco/TemporaryFile.h>
#include <Poco/Util/PropertyFileConfiguration.h>
-#include <math.h>
+#include <cstdio>
+#include <cmath>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
namespace Mantid {
namespace DataHandling {
@@ -36,9 +42,9 @@ namespace {
// number of physical detectors
constexpr size_t HORIZONTAL_TUBES = 16;
constexpr size_t VERTICAL_TUBES = 35;
-constexpr size_t DETECTORS = HORIZONTAL_TUBES + VERTICAL_TUBES;
+constexpr size_t DETECTOR_TUBES = HORIZONTAL_TUBES + VERTICAL_TUBES;
// analysed and direct detectors
-constexpr size_t HISTO_BINS_X = DETECTORS * 2;
+constexpr size_t HISTO_BINS_X = DETECTOR_TUBES * 2;
constexpr size_t HISTO_BINS_Y = 1024;
constexpr size_t HISTO_BINS_Y_DENUMERATOR = 16;
constexpr size_t PIXELS_PER_TUBE = HISTO_BINS_Y / HISTO_BINS_Y_DENUMERATOR;
@@ -55,24 +61,52 @@ constexpr size_t Progress_Total =
constexpr char FilenameStr[] = "Filename";
constexpr char MaskStr[] = "Mask";
constexpr char SelectDetectorTubesStr[] = "SelectDetectorTubes";
+constexpr char SelectDatasetStr[] = "SelectDataset";
+constexpr char OverrideDopplerFreqStr[] = "OverrideDopplerFrequency";
constexpr char OverrideDopplerPhaseStr[] = "OverrideDopplerPhase";
constexpr char FilterByTimeStartStr[] = "FilterByTimeStart";
constexpr char FilterByTimeStopStr[] = "FilterByTimeStop";
constexpr char RawDopplerTimeStr[] = "LoadAsRawDopplerTime";
+constexpr char CalibrateDopplerPhaseStr[] = "CalibrateDopplerPhase";
+constexpr char PathToBinaryStr[] = "BinaryEventPath";
// Common pairing of limits
using TimeLimits = std::pair<double, double>;
+template <typename Type>
+void AddSinglePointTimeSeriesProperty(API::LogManager &logManager,
+ const std::string &time,
+ const std::string &name,
+ const Type value) {
+ // create time series property and add single value
+ auto p = new Kernel::TimeSeriesProperty<Type>(name);
+ p->addValue(time, value);
+
+ // add to log manager
+ logManager.addProperty(p);
+}
+
+template <typename T>
+void MapNeXusToSeries(NeXus::NXEntry &entry, const std::string &path,
+ const T &defval, API::LogManager &logManager,
+ const std::string &time, const std::string &name,
+ const T &factor, int32_t index) {
+
+ T value = GetNeXusValue<T>(entry, path, defval, index);
+ AddSinglePointTimeSeriesProperty<T>(logManager, time,
+ name, value * factor);
+}
+
// Utility functions for loading values with defaults
// Single value properties only support int, double, string and bool
template <typename Type>
Type GetNeXusValue(NeXus::NXEntry &entry, const std::string &path,
- const Type &defval) {
+ const Type &defval, int32_t index) {
try {
NeXus::NXDataSetTyped<Type> dataSet = entry.openNXDataSet<Type>(path);
dataSet.load();
- return *dataSet();
+ return dataSet()[index];
} catch (std::runtime_error &) {
return defval;
}
@@ -80,12 +114,12 @@ Type GetNeXusValue(NeXus::NXEntry &entry, const std::string &path,
// string and double are special cases
template <>
double GetNeXusValue<double>(NeXus::NXEntry &entry, const std::string &path,
- const double &defval) {
+ const double &defval, int32_t index) {
try {
NeXus::NXDataSetTyped<float> dataSet = entry.openNXDataSet<float>(path);
dataSet.load();
- return *dataSet();
+ return dataSet()[index];
} catch (std::runtime_error &) {
return defval;
}
@@ -93,7 +127,7 @@ double GetNeXusValue<double>(NeXus::NXEntry &entry, const std::string &path,
template <>
std::string GetNeXusValue<std::string>(NeXus::NXEntry &entry,
const std::string &path,
- const std::string &defval) {
+ const std::string &defval, int32_t) {
try {
NeXus::NXChar dataSet = entry.openNXChar(path);
@@ -108,9 +142,9 @@ std::string GetNeXusValue<std::string>(NeXus::NXEntry &entry,
template <typename T>
void MapNeXusToProperty(NeXus::NXEntry &entry, const std::string &path,
const T &defval, API::LogManager &logManager,
- const std::string &name, const T &factor) {
+ const std::string &name, const T &factor, int32_t index) {
- T value = GetNeXusValue<T>(entry, path, defval);
+ T value = GetNeXusValue<T>(entry, path, defval, index);
logManager.addProperty<T>(name, value * factor);
}
@@ -118,9 +152,10 @@ void MapNeXusToProperty(NeXus::NXEntry &entry, const std::string &path,
template <>
void MapNeXusToProperty<std::string>(
NeXus::NXEntry &entry, const std::string &path, const std::string &defval,
- API::LogManager &logManager, const std::string &name, const std::string &) {
+ API::LogManager &logManager, const std::string &name, const std::string &,
+ int32_t index) {
- std::string value = GetNeXusValue<std::string>(entry, path, defval);
+ std::string value = GetNeXusValue<std::string>(entry, path, defval, index);
logManager.addProperty<std::string>(name, value);
}
@@ -189,13 +224,15 @@ double invert(double y, const F &f, double x0 = 0.0, const double eps = 1e-16) {
// the detectors is not constant so it needs to store the distance for each
// Note that observation time and TOF are in usec
//
+using TofData = std::tuple<double, double>;
+
class ConvertTOF {
const double m_w;
const double m_phi;
const double m_L0;
const double m_v2;
const double m_A;
- std::vector<double> m_L2;
+ const std::vector<double>& m_L2;
inline double L1(double t) const { return m_L0 + m_A * sin(m_w * t + m_phi); }
@@ -209,7 +246,7 @@ public:
: 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) {}
- double directTOF(size_t detID, double tobs) const {
+ TofData directTOF(size_t detID, double tobs) const {
// observation time and tof are in usec
auto tn = [=](double t) { return t + (L1(t) + m_L2[detID]) / v1(t); };
@@ -217,12 +254,12 @@ public:
double tsec = tobs * 1.0e-6;
double t0 = tsec - (m_L0 + m_L2[detID]) / m_v2;
double tinv = invert(tsec, tn, t0);
- double tof = m_L0 / v1(tinv) + m_L2[detID] / m_v2;
+ double tof = (m_L0 + m_L2[detID]) / v1(tinv);
- return tof * 1.0e6;
+ return TofData(tinv * 1.0e6, tof * 1.0e6);
}
- double analysedTOF(size_t detID, double tobs) const {
+ TofData analysedTOF(size_t detID, double tobs) const {
// observation time and tof are in usec
auto tn = [=](double t) { return t + L1(t) / v1(t) + m_L2[detID] / m_v2; };
@@ -231,10 +268,81 @@ public:
double t = invert(tsec, tn, t0);
double tof = m_L0 / v1(t) + m_L2[detID] / m_v2;
- return tof * 1.0e6;
+ return TofData(t * 1.0e6, tof * 1.0e6);
}
};
+// calculate mean of a subset of the vector
+double maskedMean(std::vector<double>& vec, std::vector<bool>& mask) {
+ if (vec.size() == 0 || vec.size() != mask.size())
+ throw std::runtime_error("masked mean of empty or mismatched vectors");
+ double sum = 0.0;
+ size_t count = 0;
+ for (size_t i = 0; i != vec.size(); i++) {
+ if (!mask[i])
+ continue;
+ sum += vec[i];
+ count++;
+ }
+ if (count == 0)
+ throw std::runtime_error("mean of empty vector");
+ return sum / count;
+}
+
+// calculate stdev fro a subset of the vector
+double maskedStdev(std::vector<double>& vec, std::vector<bool>& mask) {
+
+ auto avg = maskedMean(vec, mask);
+ size_t count = 0;
+ double sum = 0.0;
+ for (size_t i = 0; i != vec.size(); i++) {
+ if (!mask[i])
+ continue;
+ sum += (vec[i] - avg) * (vec[i] - avg);
+ count++;
+ }
+ return std::sqrt(sum / count);
+}
+
+// Simple reader that is compatible with the ASNTO event file loader
+class FileLoader
+{
+ std::ifstream _ifs;
+ size_t _size;
+
+public:
+ FileLoader(const char *filename)
+ : _ifs(filename, std::ios::binary | std::ios::in)
+ {
+ if (!_ifs.is_open() || _ifs.fail())
+ throw std::runtime_error("unable to open file");
+
+ _ifs.seekg(0, _ifs.end);
+ _size = _ifs.tellg();
+ _ifs.seekg(0, _ifs.beg);
+ }
+
+ bool read(char *s, std::streamsize n)
+ {
+ return static_cast<bool>(_ifs.read(s, n));
+ }
+
+ size_t size()
+ {
+ return _size;
+ }
+
+ size_t position()
+ {
+ return _ifs.tellg();
+ }
+
+ size_t selected_position()
+ {
+ return _ifs.tellg();
+ }
+};
+
} // anonymous namespace
namespace EMU {
@@ -298,6 +406,13 @@ public:
// length test in seconds
return m_framePeriod * static_cast<double>(m_frames) * 1.0e-6;
}
+
+ inline int64_t frameStart() const {
+ // returns time in nanoseconds from start of test
+ auto start = m_framePeriod * static_cast<double>(m_frames);
+ return static_cast<int64_t>(start * 1.0e3);
+ }
+
void addEvent(size_t x, size_t p, double tdop, double taux) {
// check if in time boundaries
@@ -308,7 +423,7 @@ public:
auto y = static_cast<size_t>(p / HISTO_BINS_Y_DENUMERATOR);
// determine detector id and check limits
- if (x >= DETECTORS || y >= m_stride)
+ if (x >= DETECTOR_TUBES || y >= m_stride)
return;
// map the raw detector index to the physical model
@@ -318,7 +433,7 @@ public:
// longer background chopper rate
double ptaux = fmod(taux, m_gatePeriod);
if (ptaux >= m_directTaux.first && ptaux <= m_directTaux.second)
- xid = xid + DETECTORS;
+ xid = xid + DETECTOR_TUBES;
else if (!(ptaux >= m_analysedTaux.first && ptaux <= m_analysedTaux.second))
return;
@@ -366,22 +481,33 @@ protected:
const ConvertTOF &m_convertTOF;
double m_tofMin;
double m_tofMax;
+ int64_t m_startTime;
bool m_saveAsTOF;
void addEventImpl(size_t id, size_t x, size_t, double tobs) override {
- // convert observation time to tof
+ // get the absolute time for the start of the frame
+ auto offset = m_startTime + frameStart();
+
+ // convert observation time to tof and set the pulse time
+ // relative to the start of the doppler cycle
double tof = tobs;
- if (m_saveAsTOF)
- tof = x < DETECTORS ? m_convertTOF.analysedTOF(id, tobs)
- : m_convertTOF.directTOF(id, tobs);
+ double pulse;
+ if (m_saveAsTOF) {
+ if (x < DETECTOR_TUBES)
+ std::tie(pulse, tof) = m_convertTOF.analysedTOF(id, tobs);
+ else
+ std::tie(pulse, tof) = m_convertTOF.directTOF(id, tobs);
+ offset += static_cast<int64_t>(pulse * 1e3);
+ }
if (m_tofMin > tof)
m_tofMin = tof;
if (m_tofMax < tof)
m_tofMax = tof;
- m_eventVectors[id]->push_back(tof);
+ auto ev = Types::Event::TofEvent(tof, Types::Core::DateAndTime(offset));
+ m_eventVectors[id]->push_back(ev);
}
public:
@@ -390,12 +516,14 @@ public:
const double framePeriod, const double gatePeriod,
const TimeLimits &timeBoundary, const TimeLimits &directLimits,
const TimeLimits &analysedLimits, ConvertTOF &convert,
- std::vector<EventVector_pt> &eventVectors, bool saveAsTOF)
+ std::vector<EventVector_pt> &eventVectors, int64_t startTime,
+ bool saveAsTOF)
: EventProcessor(roi, mapIndex, stride, framePeriod, gatePeriod,
timeBoundary, directLimits, analysedLimits),
m_eventVectors(eventVectors), m_convertTOF(convert),
m_tofMin(std::numeric_limits<double>::max()),
- m_tofMax(std::numeric_limits<double>::min()), m_saveAsTOF(saveAsTOF) {}
+ m_tofMax(std::numeric_limits<double>::min()),
+ m_startTime(startTime), m_saveAsTOF(saveAsTOF) {}
double tofMin() const { return m_tofMin <= m_tofMax ? m_tofMin : 0.0; }
double tofMax() const { return m_tofMin <= m_tofMax ? m_tofMax : 0.0; }
@@ -403,7 +531,24 @@ public:
} // namespace EMU
// register the algorithm into the AlgorithmFactory
-DECLARE_FILELOADER_ALGORITHM(LoadEMU)
+DECLARE_FILELOADER_ALGORITHM(LoadEMUTar)
+DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadEMUHdf)
+
+const std::string LoadEMU<Kernel::FileDescriptor>::name() const {
+ return "LoadEMUTar";
+}
+
+const std::string LoadEMU<Kernel::FileDescriptor>::summary() const {
+ return "Loads a merged EMU Hdf and event file into a workspace.";
+}
+
+const std::string LoadEMU<Kernel::NexusDescriptor>::name() const {
+ return "LoadEMUHdf";
+}
+
+const std::string LoadEMU<Kernel::NexusDescriptor>::summary() const {
+ return "Loads an EMU Hdf and linked event file into a workspace.";
+}
/**
* Return the confidence value that this algorithm can load the file
@@ -411,7 +556,8 @@ DECLARE_FILELOADER_ALGORITHM(LoadEMU)
* @returns An integer specifying the confidence level. 0 indicates it will not
* be used
*/
-int LoadEMU::confidence(Kernel::FileDescriptor &descriptor) const {
+template <>
+int LoadEMU<Kernel::FileDescriptor>::confidence(Kernel::FileDescriptor &descriptor) const {
if (descriptor.extension() != ".tar")
return 0;
@@ -437,23 +583,77 @@ int LoadEMU::confidence(Kernel::FileDescriptor &descriptor) const {
return (hdfFiles == 1) && (binFiles == 1) ? 50 : 0;
}
+/**
+* Return the confidence value that this algorithm can load the file
+* @param descriptor A descriptor for the file
+* @returns An integer specifying the confidence level. 0 indicates it will not
+* be used
+*/
+template <>
+int LoadEMU<Kernel::NexusDescriptor>::confidence(Kernel::NexusDescriptor &descriptor) const {
+ if (descriptor.extension() != ".hdf")
+ return 0;
+
+ if (descriptor.pathExists("/entry1/site_name") &&
+ descriptor.pathExists("/entry1/instrument/doppler/ctrl/velocity") &&
+ descriptor.pathExists("/entry1/instrument/doppler/ctrl/amplitude") &&
+ descriptor.pathExists("/entry1/instrument/detector/daq_dirname") &&
+ descriptor.pathExists("/entry1/instrument/detector/dataset_number") &&
+ descriptor.pathExists("/entry1/data/hmm_total_t_ds0") &&
+ descriptor.pathExists("/entry1/data/hmm_total_t_ds1") &&
+ descriptor.pathExists("/entry1/data/hmm_total_xt_ds0") &&
+ descriptor.pathExists("/entry1/data/hmm_total_xt_ds1")) {
+ return 80;
+ }
+ else {
+ return 0;
+ }
+}
+
/**
* Initialise the algorithm. Declare properties which can be set before
- * execution (input) and
- * read from after the execution (output).
+ * execution (input) and read from after the execution (output).
+ */
+template <>
+void LoadEMU<Kernel::FileDescriptor>::init() {
+ init(false);
+}
+
+/**
+* Initialise the algorithm. Declare properties which can be set before
+* execution (input) and read from after the execution (output).
+*/
+template <>
+void LoadEMU<Kernel::NexusDescriptor>::init() {
+ init(true);
+}
+
+/**
+ * Declares the properties for the two loader variants. Adds the path option
+ * path to the binary file if it is the hdf loader.
*/
-void LoadEMU::init() {
+template <typename FD>
+void LoadEMU<FD>::init(bool hdfLoader) {
// Specify file extensions which can be associated with a specific file.
std::vector<std::string> exts;
// Declare the Filename algorithm property. Mandatory. Sets the path to the
// file to load.
exts.clear();
- exts.emplace_back(".tar");
+ if (hdfLoader)
+ exts.emplace_back(".hdf");
+ else
+ exts.emplace_back(".tar");
declareProperty(Kernel::make_unique<API::FileProperty>(
FilenameStr, "", API::FileProperty::Load, exts),
"The input filename of the stored data");
+ if (hdfLoader) {
+ declareProperty(PathToBinaryStr, "",
+ "Relative or absolute path to the compressed binary\n"
+ "event file linked to the HDF file, eg /storage/data/");
+ }
+
// mask
exts.clear();
exts.emplace_back(".xml");
@@ -463,17 +663,29 @@ void LoadEMU::init() {
declareProperty(SelectDetectorTubesStr, "",
"Comma separated range of detectors tubes to be loaded,\n"
- " e.g. 16,19-45,47");
+ " eg 16,19-45,47");
declareProperty(
Kernel::make_unique<API::WorkspaceProperty<API::IEventWorkspace>>(
"OutputWorkspace", "", Kernel::Direction::Output));
+ if (hdfLoader) {
+ declareProperty(SelectDatasetStr, 0,
+ "Select the index for the dataset to be loaded.");
+ }
+
+ declareProperty(OverrideDopplerFreqStr, EMPTY_DBL(),
+ "Override the Doppler frequency, in Hertz.");
+
declareProperty(OverrideDopplerPhaseStr, EMPTY_DBL(),
"Override the Doppler phase, in degrees.");
+ declareProperty(CalibrateDopplerPhaseStr, false,
+ "Calibrate the Doppler phase prior to TOF conversion,\n"
+ "ignored if imported as Doppler time or phase entered");
+
declareProperty(RawDopplerTimeStr, false,
- "Import file as observed time relative the Doppler "
+ "Import file as observed time relative the Doppler\n"
"drive, in microsecs.");
declareProperty(FilterByTimeStartStr, 0.0,
@@ -492,7 +704,8 @@ void LoadEMU::init() {
/**
* Creates an event workspace and sets the title.
*/
-void LoadEMU::createWorkspace(ANSTO::Tar::File &tarFile) {
+template <typename FD>
+void LoadEMU<FD>::createWorkspace(const std::string &title) {
// Create the workspace
m_localWorkspace = boost::make_shared<DataObjects::EventWorkspace>();
@@ -504,53 +717,162 @@ void LoadEMU::createWorkspace(ANSTO::Tar::File &tarFile) {
m_localWorkspace->setYUnit("Counts");
// set title
- const std::vector<std::string> &subFiles = tarFile.files();
- for (const auto &subFile : subFiles)
- if (subFile.compare(0, 3, "EMU") == 0) {
- std::string title = subFile;
+ m_localWorkspace->setTitle(title);
+}
- if (title.rfind(".hdf") == title.length() - 4)
- title.resize(title.length() - 4);
+/**
+ * Execute the algorithm. Extracts the hdf and event file from the tar
+ * and invokes the invokes the common exec() function that works with
+ * the two files.
+ */
+template <>
+void LoadEMU<Kernel::FileDescriptor>::exec() {
+ /**
+ * Opens the tar and extracts the hdf and event data to temporary files
+ */
+ std::string filename = getPropertyValue(FilenameStr);
+ ANSTO::Tar::File tarFile(filename);
+ if (!tarFile.good())
+ throw std::invalid_argument("invalid EMU tar file");
- if (title.rfind(".nx") == title.length() - 3)
- title.resize(title.length() - 3);
+ // dataset selection not supported in tar version - order is not guaranteed
+ m_datasetIndex = 0;
- m_localWorkspace->setTitle(title);
- break;
- }
+ // lambda functions to find the first file of extension and to extract
+ // the file
+ const std::vector<std::string> &files = tarFile.files();
+ auto selectFile = [&](const std::string& ext) {
+ auto itf =
+ std::find_if(files.cbegin(), files.cend(), [&ext](const std::string &file) {
+ return file.rfind(ext) == file.length() - 4;
+ });
+ if (itf == files.end())
+ throw std::runtime_error("missing tar file data");
+ else
+ tarFile.select(itf->c_str());
+ };
+ auto extractFile = [&](Poco::TemporaryFile& tfile) {
+
+ boost::shared_ptr<FILE> handle(fopen(tfile.path().c_str(), "wb"), fclose);
+ if (handle) {
+ // copy content
+ char buffer[4096];
+ size_t bytesRead;
+ while (0 != (bytesRead = tarFile.read(buffer, sizeof(buffer))))
+ fwrite(buffer, bytesRead, 1, handle.get());
+ handle.reset();
+ }
+ };
+
+ // extract hdf file into tmp file
+ selectFile(".hdf");
+ Poco::TemporaryFile hdfFile;
+ extractFile(hdfFile);
+
+ // extract the event file
+ selectFile(".bin");
+ Poco::TemporaryFile eventFile;
+ extractFile(eventFile);
+
+ // call the common loader
+ exec(hdfFile.path(), eventFile.path());
}
/**
- * Execute the algorithm. The steps involved are:
- * Get the instrument properties and load options
- * Create the workspace
- * Load the instrument from the IDF
- * Reposition the relevant neutronic values for model based on the parameters
- * Load the data values and convert to TOF
- * Setting up the masks
- */
-void LoadEMU::exec() {
+* Execute the algorithm. Establishes the filepath to the event file
+* from the HDF link and the path provided and invokes the common
+* exec() function that works with the two files.
+*/
+template <>
+void LoadEMU<Kernel::NexusDescriptor>::exec() {
+
+ namespace fs = boost::filesystem;
+
+ // Open the hdf file and find the dirname and dataset number
+ std::string hdfFile = getPropertyValue(FilenameStr);
+ std::string evtPath = getPropertyValue(PathToBinaryStr);
+ if (evtPath.empty())
+ evtPath = "./";
+
+ // if relative ./ or ../ then append to the directory for the hdf file
+ if (evtPath.rfind("./") == 0 || evtPath.rfind("../") == 0) {
+ fs::path hp = hdfFile;
+ evtPath = fs::canonical(evtPath, hp.parent_path()).generic_string();
+ }
+
+ // dataset index to be loaded
+ m_datasetIndex = getProperty(SelectDatasetStr);
+
+ // if path provided build the file path
+ if (fs::is_directory(evtPath))
+ {
+ NeXus::NXRoot root(hdfFile);
+ NeXus::NXEntry entry = root.openFirstEntry();
+ auto eventDir = GetNeXusValue<std::string>(entry,
+ "instrument/detector/daq_dirname", "./", 0);
+ auto dataset = GetNeXusValue<int32_t>(entry,
+ "instrument/detector/dataset_number", 0,
+ m_datasetIndex);
+
+ // build the path to the event file as if a relative or absolute
+ // path is passed:
+ // 'relpath/[daq_dirname]/DATASET_[n]/EOS.bin' or the
+ char buffer[255] = {};
+ snprintf(buffer, sizeof(buffer), "%s/DATASET_%d/EOS.bin",
+ eventDir.c_str(), dataset);
+ fs::path path = evtPath;
+ path /= buffer;
+ path = fs::absolute(path);
+ evtPath = path.generic_string();
+ }
+
+ // finally check that the event file exists
+ if (!fs::is_regular_file(evtPath)) {
+ std::string msg = "Check path, cannot open binary event file: " + evtPath;
+ throw std::runtime_error(msg);
+ }
+
+ exec(hdfFile, evtPath);
+}
+
+/**
+* Execute the algorithm. The steps involved are:
+* Create the workspace
+* Get the instrument properties and load options
+* Load the instrument from the IDF
+* Reposition the relevant neutronic values for model based on the parameters
+* Load the data values and convert to TOF
+* Setting up the masks
+*/
+template <typename FD>
+void LoadEMU<FD>::exec(const std::string& hdfFile, const std::string& eventFile) {
+
+ namespace fs = boost::filesystem;
// Create workspace
// ----------------
- std::string filename = getPropertyValue(FilenameStr);
- ANSTO::Tar::File tarFile(filename);
- if (!tarFile.good())
- throw std::invalid_argument("invalid EMU file");
- createWorkspace(tarFile);
+ fs::path p = hdfFile;
+ for (; !p.extension().empty();)
+ p = p.stem();
+ std::string title = p.generic_string();
+ createWorkspace(title);
API::LogManager &logManager = m_localWorkspace->mutableRun();
API::Progress prog(this, 0.0, 1.0, Progress_Total);
// Load instrument and workspace properties
// ----------------------------------------
- loadParameters(tarFile, logManager);
+ logManager.addProperty(SelectDatasetStr, m_datasetIndex);
+ loadParameters(hdfFile, logManager);
prog.doReport("creating instrument");
loadInstrument();
- // Get the region of interest and filters
+ // Get the region of interest and filters and save to log
//
std::string maskfile = getPropertyValue(MaskStr);
std::string seltubes = getPropertyValue(SelectDetectorTubesStr);
+ logManager.addProperty(SelectDetectorTubesStr, seltubes);
+ logManager.addProperty(MaskStr, maskfile);
+
std::vector<bool> roi = createRoiVector(seltubes, maskfile);
double timeMaxBoundary = getProperty(FilterByTimeStopStr);
if (isEmpty(timeMaxBoundary))
@@ -568,14 +890,14 @@ void LoadEMU::exec() {
// sequence starting from 0
//
double sampleAnalyser = iparam("SampleAnalyser");
- auto endID = static_cast<detid_t>(DETECTORS * PIXELS_PER_TUBE);
+ auto endID = static_cast<detid_t>(DETECTOR_TUBES * PIXELS_PER_TUBE);
for (detid_t detID = 0; detID < endID; detID++)
updateNeutronicPostions(detID, sampleAnalyser);
// get the detector map from raw input to a physical detector
//
std::string dmapStr = instr->getParameterAsString("DetectorMap");
- std::vector<size_t> detMapIndex = std::vector<size_t>(DETECTORS, 0);
+ std::vector<size_t> detMapIndex = std::vector<size_t>(DETECTOR_TUBES, 0);
mapRangeToIndex(dmapStr, detMapIndex, [](size_t n) { return n; });
// Collect the L2 distances, Doppler characteristics and
@@ -587,12 +909,8 @@ void LoadEMU::exec() {
double framePeriod =
1.0e6 / m_dopplerFreq; // period and max direct as microsec
double sourceSample = iparam("SourceSample");
- ConvertTOF convertTOF(m_dopplerAmpl, m_dopplerFreq, m_dopplerPhase,
- sourceSample, v2, m_detectorL2);
-
- Types::Core::DateAndTime start_time(
- logManager.getPropertyValueAsType<std::string>("StartTime"));
- std::string time_str = start_time.toISO8601String();
+ ConvertTOF convertTOF(m_dopplerAmpl * m_dopplerRun, m_dopplerFreq,
+ m_dopplerPhase, sourceSample, v2, m_detectorL2);
// Load the events file
// --------------------
@@ -614,48 +932,73 @@ void LoadEMU::exec() {
1000.0 * iparam("AnalysedTauxMax"));
// fabs because the value can be negative
- double gatePeriod = 1.0e6 / fabs(logManager.getPropertyValueAsType<double>(
- "GraphiteChopperFrequency"));
+ double gatePeriod = 1.0e6 / fabs(logManager.getTimeSeriesProperty<double>(
+ "GraphiteChopperFrequency")->firstValue());
// count total events per pixel and reserve necessary memory
EMU::EventCounter eventCounter(roi, detMapIndex, PIXELS_PER_TUBE, framePeriod,
gatePeriod, timeBoundary, directLimits,
analysedLimits, eventCounts);
- loadEvents(prog, "loading neutron counts", tarFile, eventCounter);
+ loadEvents(prog, "loading neutron counts", eventFile, eventCounter);
ANSTO::ProgressTracker progTracker(prog, "creating neutron event lists",
numberHistograms, Progress_ReserveMemory);
prepareEventStorage(progTracker, eventCounts, eventVectors);
- // now perfrom the actual event collection and TOF convert if necessary
- bool saveAsTOF = !getProperty(RawDopplerTimeStr);
+ // now perform the actual event collection and TOF convert if necessary
+ // if a phase calibration is required then load it as raw doppler time
+ // perform the calibration and then convert to TOF
+ Types::Core::DateAndTime startTime(m_startRun);
+ auto start_nanosec = startTime.totalNanoseconds();
+ bool saveAsTOF = !getProperty(RawDopplerTimeStr);
+ bool loadAsTOF = !m_calibrateDoppler && saveAsTOF;
EMU::EventAssigner eventAssigner(
roi, detMapIndex, PIXELS_PER_TUBE, framePeriod, gatePeriod, timeBoundary,
- directLimits, analysedLimits, convertTOF, eventVectors, saveAsTOF);
- loadEvents(prog, "loading neutron events (TOF)", tarFile, eventAssigner);
+ directLimits, analysedLimits, convertTOF, eventVectors,
+ start_nanosec, loadAsTOF);
+ loadEvents(prog, "loading neutron events (TOF)", eventFile, eventAssigner);
+
+ // perform a calibration and then TOF conversion if necessary
+ // and update the tof limits
+ auto minTOF = eventAssigner.tofMin();
+ auto maxTOF = eventAssigner.tofMax();
+ if (m_calibrateDoppler) {
+ calibrateDopplerPhase(eventCounts, eventVectors);
+ if (saveAsTOF) {
+ dopplerTimeToTOF(eventVectors, minTOF, maxTOF);
+ }
+ }
+ AddSinglePointTimeSeriesProperty<double>(logManager, m_startRun,
+ "DopplerPhase", m_dopplerPhase);
// just to make sure the bins hold it all and setup the detector masks
m_localWorkspace->setAllX(
- HistogramData::BinEdges{std::max(0.0, floor(eventAssigner.tofMin())),
- eventAssigner.tofMax() + 1});
+ HistogramData::BinEdges{std::max(0.0, floor(minTOF)), maxTOF + 1});
setupDetectorMasks(roi);
// set log values
auto frame_count = static_cast<int>(eventCounter.numFrames());
+ AddSinglePointTimeSeriesProperty<int>(logManager, m_startRun,
+ "frame_count", frame_count);
+ std::string filename = getPropertyValue(FilenameStr);
logManager.addProperty("filename", filename);
- logManager.addProperty("frame_count", frame_count);
Types::Core::time_duration duration = boost::posix_time::microseconds(
static_cast<boost::int64_t>(eventCounter.duration() * 1.0e6));
- Types::Core::DateAndTime end_time(start_time + duration);
- logManager.addProperty("start_time", start_time.toISO8601String());
- logManager.addProperty("end_time", end_time.toISO8601String());
+ Types::Core::DateAndTime endTime(startTime + duration);
+ logManager.addProperty("start_time", startTime.toISO8601String());
+ logManager.addProperty("end_time", endTime.toISO8601String());
+ logManager.addProperty<double>("dur", eventCounter.duration());
+
+ // Finally add the time-series parameter explicitly
+ loadEnvironParameters(hdfFile, logManager);
setProperty("OutputWorkspace", m_localWorkspace);
}
// set up the detector masks
-void LoadEMU::setupDetectorMasks(std::vector<bool> &roi) {
+template <typename FD>
+void LoadEMU<FD>::setupDetectorMasks(std::vector<bool> &roi) {
// count total number of masked bins
size_t maskedBins = 0;
@@ -680,7 +1023,8 @@ void LoadEMU::setupDetectorMasks(std::vector<bool> &roi) {
}
// prepare the event storage
-void LoadEMU::prepareEventStorage(ANSTO::ProgressTracker &progTracker,
+template <typename FD>
+void LoadEMU<FD>::prepareEventStorage(ANSTO::ProgressTracker &progTracker,
std::vector<size_t> &eventCounts,
std::vector<EventVector_pt> &eventVectors) {
@@ -702,12 +1046,23 @@ void LoadEMU::prepareEventStorage(ANSTO::ProgressTracker &progTracker,
}
// get and log the Doppler parameters
-void LoadEMU::loadDopplerParameters(API::LogManager &logm) {
+template <typename FD>
+void LoadEMU<FD>::loadDopplerParameters(API::LogManager &logm) {
auto instr = m_localWorkspace->getInstrument();
- m_dopplerFreq = logm.getPropertyValueAsType<double>("DopplerFrequency");
- m_dopplerAmpl = logm.getPropertyValueAsType<double>("DopplerAmplitude");
+
+ // use nominal frequency based on amp and velocity unless overridden
+ m_dopplerAmpl = logm.getTimeSeriesProperty<double>("DopplerAmplitude")->firstValue();
+ m_dopplerRun = logm.getTimeSeriesProperty<int32_t>("DopplerRun")->firstValue();
+ m_dopplerFreq = getProperty(OverrideDopplerFreqStr);
+ if (isEmpty(m_dopplerFreq)) {
+ auto doppVel = logm.getTimeSeriesProperty<double>("DopplerVelocity")->firstValue();
+ m_dopplerFreq = 0.5 * doppVel / (M_PI * m_dopplerAmpl);
+ }
+ AddSinglePointTimeSeriesProperty<double>(logm, m_startRun, "DopplerFrequency", m_dopplerFreq);
+
m_dopplerPhase = getProperty(OverrideDopplerPhaseStr);
+ m_calibrateDoppler = getProperty(CalibrateDopplerPhaseStr) && isEmpty(m_dopplerPhase);
if (isEmpty(m_dopplerPhase)) {
// sinusoidal motion crossing a threshold with a delay
double doppThreshold =
@@ -717,11 +1072,135 @@ void LoadEMU::loadDopplerParameters(API::LogManager &logm) {
180.0 - asin(0.001 * doppThreshold / m_dopplerAmpl) * 180.0 / M_PI +
doppDelay * m_dopplerFreq;
}
- logm.addProperty<double>("DopplerPhase", m_dopplerPhase);
+
+ // problems adding 'bool' to log
+ int32_t calPhase = 1 ? m_calibrateDoppler : 0;
+ logm.addProperty("CalibratePhase", calPhase);
+}
+
+// calibrate doppler phase based on teh analysed events
+template <typename FD>
+void LoadEMU<FD>::calibrateDopplerPhase(std::vector<size_t> &eventCounts,
+ std::vector<EventVector_pt> &eventVectors) {
+
+ // get the doppler parameters
+ auto instr = m_localWorkspace->getInstrument();
+ double v2 = instr->getNumberParameter("AnalysedV2")[0];
+ double l1 = instr->getNumberParameter("SourceSample")[0];
+
+ // get the number of analysed events and initial doppler time
+ auto startID = static_cast<size_t>(HORIZONTAL_TUBES * PIXELS_PER_TUBE);
+ auto endID = static_cast<size_t>(DETECTOR_TUBES * PIXELS_PER_TUBE);
+ size_t numEvents = 0;
+ for (size_t i = startID; i < endID; i++)
+ numEvents += eventCounts[i];
+ if (numEvents == 0)
+ throw std::runtime_error("no analysed events for phase calibration");
+ std::vector<double> nVel(numEvents);
+ std::vector<size_t> nMap(numEvents);
+ std::vector<bool> nCnd(numEvents);
+ constexpr size_t NHIST = 100;
+ std::vector<int> histogram(NHIST+1, 0);
+
+ // define the cost function to optimize phase
+ auto costFn = [&, this](double phase) {
+
+ ConvertTOF convTOF(m_dopplerAmpl * m_dopplerRun, m_dopplerFreq, phase,
+ l1, v2, m_detectorL2);
+
+ // convert each analysed event to source velocity
+ size_t ix = 0;
+ double tof, pulse;
+ for (size_t i = startID; i < endID; i++) {
+ for (auto const& x : *eventVectors[i]) {
+ std::tie(pulse, tof) = convTOF.analysedTOF(i, x.tof());
+ auto tof1 = 1e-6 * tof - m_detectorL2[i] / v2;
+ nVel[ix++] = l1 / tof1;
+ }
+ }
+
+ // now histogram the data and create the map from velocity to hist
+ auto ixlim = std::minmax_element(nVel.begin(), nVel.end());
+ auto vmin = nVel[ixlim.first - nVel.begin()];
+ auto vmax = nVel[ixlim.second - nVel.begin()];
+ int maxHist = 0;
+ std::fill(histogram.begin(), histogram.end(), 0);
+ auto delta = (vmax - vmin) / NHIST;
+ for (int i = 0; i < numEvents; i++) {
+ auto v = nVel[i];
+ auto ix = static_cast<size_t>(std::floor((v - vmin) / delta));
+ histogram[ix]++;
+ if (histogram[ix] > maxHist)
+ maxHist = histogram[ix];
+ nMap[i] = ix;
+ }
+
+ // determine the points above the 25% threshold
+ int minLevel = static_cast<int>(maxHist / 4);
+ for (size_t i = 0; i < numEvents; i++) {
+ nCnd[i] = (histogram[nMap[i]] >= minLevel ? true : false);
+ }
+
+ // calculate the standard deviation for the points above the threshold
+ auto cost = maskedStdev(nVel, nCnd);
+ return cost;
+ };
+
+ // call the optimizer and update the doppler phase value
+ // limit the optimization to 30 iterations
+ int bits = std::numeric_limits<double>::digits;
+ boost::uintmax_t itn = 30;
+ using boost::math::tools::brent_find_minima;
+ auto minPhase = m_dopplerPhase - 5.0;
+ auto maxPhase = m_dopplerPhase + 5.0;
+ auto r = brent_find_minima(costFn, minPhase, maxPhase, bits, itn);
+ m_dopplerPhase = r.first;
+}
+
+// convert the doppler time to TOF for all the events in the workspace
+template <typename FD>
+void LoadEMU<FD>::dopplerTimeToTOF(std::vector<EventVector_pt> &eventVectors,
+ double& minTOF, double& maxTOF) {
+
+ // get the doppler parameters and initialise TOD converter
+ auto instr = m_localWorkspace->getInstrument();
+ double v2 = instr->getNumberParameter("AnalysedV2")[0];
+ double l1 = instr->getNumberParameter("SourceSample")[0];
+ ConvertTOF convTOF(m_dopplerAmpl * m_dopplerRun, m_dopplerFreq,
+ m_dopplerPhase, l1, v2, m_detectorL2);
+
+ // run through all the events noting that analysed event are in
+ // the bottom half of the detector ids
+ auto start = true;
+ auto directID = static_cast<size_t>(DETECTOR_TUBES * PIXELS_PER_TUBE);
+ for (size_t id = 0; id < eventVectors.size(); id++) {
+ for (auto &x : *eventVectors[id]) {
+ double tof, pulse;
+ if (id < directID)
+ std::tie(pulse, tof) = convTOF.analysedTOF(id, x.tof());
+ else
+ std::tie(pulse, tof) = convTOF.directTOF(id, x.tof());
+
+ // update the pulse time and tof
+ int64_t pulseTime = x.pulseTime().totalNanoseconds();
+ pulseTime += static_cast<int64_t>(pulse * 1000);
+ x = Types::Event::TofEvent(tof, Types::Core::DateAndTime(pulseTime));
+
+ if (start) {
+ minTOF = maxTOF = x.tof();
+ start = false;
+ }
+ else {
+ minTOF = std::min(minTOF, x.tof());
+ maxTOF = std::max(maxTOF, x.tof());
+ }
+ }
+ }
}
// Recovers the L2 neutronic distance for each detector.
-void LoadEMU::loadDetectorL2Values() {
+template <typename FD>
+void LoadEMU<FD>::loadDetectorL2Values() {
m_detectorL2 = std::vector<double>(HISTOGRAMS);
const auto &detectorInfo = m_localWorkspace->detectorInfo();
@@ -734,14 +1213,8 @@ void LoadEMU::loadDetectorL2Values() {
}
// update the neutronic positins
-void LoadEMU::updateNeutronicPostions(detid_t detID, double sampleAnalyser) {
-
- // get the instrument
-
- // get the list of indirect horizontal detectors
-
- // for each detector get the current position and scale
- // the detectors
+template <typename FD>
+void LoadEMU<FD>::updateNeutronicPostions(detid_t detID, double sampleAnalyser) {
Geometry::Instrument_const_sptr instrument =
m_localWorkspace->getInstrument();
@@ -765,7 +1238,8 @@ void LoadEMU::updateNeutronicPostions(detid_t detID, double sampleAnalyser) {
// region of interest is defined by the selected detectors and the mask
//
-std::vector<bool> LoadEMU::createRoiVector(const std::string &selected,
+template <typename FD>
+std::vector<bool> LoadEMU<FD>::createRoiVector(const std::string &selected,
const std::string &maskfile) {
std::vector<bool> result(HISTOGRAMS, true);
@@ -804,83 +1278,88 @@ std::vector<bool> LoadEMU::createRoiVector(const std::string &selected,
return result;
}
-// load parameters from input file
-void LoadEMU::loadParameters(ANSTO::Tar::File &tarFile, API::LogManager &logm) {
+// load parameters from input hdf file
+template <typename FD>
+void LoadEMU<FD>::loadParameters(const std::string& hdfFile, API::LogManager &logm) {
- // extract log and hdf file
- const std::vector<std::string> &files = tarFile.files();
- auto file_it =
- std::find_if(files.cbegin(), files.cend(), [](const std::string &file) {
- return file.rfind(".hdf") == file.length() - 4;
- });
- if (file_it != files.end()) {
- tarFile.select(file_it->c_str());
- // extract hdf file into tmp file
- Poco::TemporaryFile hdfFile;
- boost::shared_ptr<FILE> handle(fopen(hdfFile.path().c_str(), "wb"), fclose);
- if (handle) {
- // copy content
- char buffer[4096];
- size_t bytesRead;
- while (0 != (bytesRead = tarFile.read(buffer, sizeof(buffer))))
- fwrite(buffer, bytesRead, 1, handle.get());
- handle.reset();
-
- NeXus::NXRoot root(hdfFile.path());
+ NeXus::NXRoot root(hdfFile);
NeXus::NXEntry entry = root.openFirstEntry();
MapNeXusToProperty<std::string>(entry, "sample/name", "unknown", logm,
- "SampleName", "");
+ "SampleName", "", 0);
MapNeXusToProperty<std::string>(entry, "sample/description", "unknown",
- logm, "SampleDescription", "");
- MapNeXusToProperty<std::string>(
- entry, "start_time", "2000-01-01T00:00:00", logm, "StartTime", "");
-
- MapNeXusToProperty<double>(entry, "instrument/doppler/ctrl/amplitude",
- 75.0, logm, "DopplerAmplitude", 0.001);
- double speedToFreq =
- 0.5 / (M_PI * m_localWorkspace->run().getPropertyValueAsType<double>(
- "DopplerAmplitude"));
- MapNeXusToProperty<double>(entry, "instrument/doppler/ctrl/velocity", 4.7,
- logm, "DopplerFrequency", speedToFreq);
-
- MapNeXusToProperty<double>(entry, "instrument/chpr/background/actspeed",
- 1272.8, logm, "BackgroundChopperFrequency",
- 1.0 / 60);
- MapNeXusToProperty<double>(entry, "instrument/chpr/graphite/actspeed",
- 2545.6, logm, "GraphiteChopperFrequency",
- 1.0 / 60);
- MapNeXusToProperty<double>(entry, "instrument/hztubegap", 0.02, logm,
- "HorizontalTubesGap", 1.0);
- MapNeXusToProperty<int32_t>(entry, "monitor/bm1_counts", 0, logm,
- "MonitorCounts", 1);
-
- // temp fix for source position when loading IDF
+ logm, "SampleDescription", "", 0);
+
+ // if dataset index > 0 need to add an offset to the start time
+ Types::Core::DateAndTime startTime(GetNeXusValue<std::string>(
+ entry, "start_time",
+ "2000-01-01T00:00:00", 0));
+ if (m_datasetIndex > 0) {
+ auto baseTime = GetNeXusValue<int32_t>(entry,
+ "instrument/detector/start_time",
+ 0, 0);
+ auto nthTime = GetNeXusValue<int32_t>(entry,
+ "instrument/detector/start_time",
+ 0, m_datasetIndex);
+
+ Types::Core::time_duration duration = boost::posix_time::microseconds(
+ static_cast<boost::int64_t>((nthTime - baseTime) * 1.0e6));
+ Types::Core::DateAndTime startDataset(startTime + duration);
+ m_startRun = startDataset.toISO8601String();
+ }
+ else {
+ m_startRun = startTime.toISO8601String();
+ }
+
+ MapNeXusToSeries<double>(entry, "instrument/doppler/ctrl/amplitude",
+ 75.0, logm, m_startRun, "DopplerAmplitude",
+ 0.001, m_datasetIndex);
+ MapNeXusToSeries<double>(entry, "instrument/doppler/ctrl/velocity", 4.7,
+ logm, m_startRun, "DopplerVelocity", 1, m_datasetIndex);
+ MapNeXusToSeries<int>(entry, "instrument/doppler/ctrl/run_cmd",
+ 1, logm, m_startRun, "DopplerRun", 1, m_datasetIndex);
+
+ MapNeXusToSeries<double>(entry, "instrument/chpr/background/actspeed",
+ 1272.8, logm, m_startRun, "BackgroundChopperFrequency",
+ 1.0 / 60, 0);
+ MapNeXusToSeries<double>(entry, "instrument/chpr/graphite/actspeed",
+ 2545.6, logm, m_startRun, "GraphiteChopperFrequency",
+ 1.0 / 60, 0);
+ // hz tube gap or equivalent to be added later - reverts to default
+ MapNeXusToSeries<double>(entry, "instrument/hztubegap", 0.02, logm, m_startRun,
+ "horizontal_tubes_gap", 1.0, 0);
+ MapNeXusToSeries<int32_t>(entry, "monitor/bm1_counts", 0, logm, m_startRun,
+ "MonitorCounts", 1, m_datasetIndex);
+
+ // fix for source position when loading IDF
MapNeXusToProperty<double>(entry, "instrument/doppler/tosource", 2.035,
- logm, "SourceSample", 1.0);
- }
- }
+ logm, "SourceSample", 1.0, 0);
+}
- // patching
- file_it = std::find(files.cbegin(), files.cend(), "History.log");
- if (file_it != files.cend()) {
- tarFile.select(file_it->c_str());
- std::string logContent;
- logContent.resize(tarFile.selected_size());
- tarFile.read(&logContent[0], logContent.size());
- std::istringstream data(logContent);
- Poco::AutoPtr<Poco::Util::PropertyFileConfiguration> conf(
- new Poco::Util::PropertyFileConfiguration(data));
-
- if (conf->hasProperty("SampleName")) {
- std::string name = conf->getString("SampleName");
- logm.addProperty("SampleName", name);
- }
- }
+// the environment variable needs to be loaded at the end
+// as the
+template <typename FD>
+void LoadEMU<FD>::loadEnvironParameters(const std::string &hdfFile, API::LogManager &logm) {
+
+ NeXus::NXRoot root(hdfFile);
+ NeXus::NXEntry entry = root.openFirstEntry();
+ auto time_str = logm.getPropertyValueAsType<std::string>("end_time");
+
+ // load the environment variables for teh dataset loaded
+ std::vector<std::string> tags = {
+ "P01PS03", "P01PSP03", "T01S00", "T01S05", "T01S06", "T01SP00",
+ "T01SP06", "T02S00", "T02S04", "T02S05", "T02S06",
+ "T02SP00", "T02SP06" };
+
+ for (const auto &tag : tags) {
+ MapNeXusToSeries<double>(entry, "data/" + tag, 0.0, logm, time_str,
+ "env_" + tag, 1.0, m_datasetIndex);
+ }
}
// load the instrument definition and instrument parameters
-void LoadEMU::loadInstrument() {
+template <typename FD>
+void LoadEMU<FD>::loadInstrument() {
// loads the IDF and parameter file
API::IAlgorithm_sptr loadInstrumentAlg =
@@ -893,30 +1372,28 @@ void LoadEMU::loadInstrument() {
}
// read counts/events from binary file
+template <typename FD>
template <class EventProcessor>
-void LoadEMU::loadEvents(API::Progress &prog, const char *progMsg,
- ANSTO::Tar::File &tarFile,
+void LoadEMU<FD>::loadEvents(API::Progress &prog, const char *progMsg,
+ const std::string& eventFile,
EventProcessor &eventProcessor) {
using namespace ANSTO;
prog.doReport(progMsg);
- // select bin file
- int64_t fileSize = 0;
- const std::vector<std::string> &files = tarFile.files();
- for (const auto &file : files)
- if (file.rfind(".bin") == file.length() - 4) {
- tarFile.select(file.c_str());
- fileSize = tarFile.selected_size();
- break;
- }
+ FileLoader loader(eventFile.c_str());
// for progress notifications
- ANSTO::ProgressTracker progTracker(prog, progMsg, fileSize,
+ ANSTO::ProgressTracker progTracker(prog, progMsg, loader.size(),
Progress_LoadBinFile);
- ReadEventFile(tarFile, eventProcessor, progTracker, 100, false);
+ ReadEventFile(loader, eventProcessor, progTracker, 100, false);
}
+
+// finally instantiate the two loaders
+template class LoadEMU<Kernel::FileDescriptor>;
+template class LoadEMU<Kernel::NexusDescriptor>;
+
} // namespace DataHandling
} // namespace Mantid
\ No newline at end of file
diff --git a/Framework/DataHandling/test/LoadEMUauTest.h b/Framework/DataHandling/test/LoadEMUauTest.h
index 4c9dffc9726..ceb47bd38b0 100644
--- a/Framework/DataHandling/test/LoadEMUauTest.h
+++ b/Framework/DataHandling/test/LoadEMUauTest.h
@@ -23,14 +23,14 @@ using namespace Mantid::DataObjects;
class LoadEMUauTest : public CxxTest::TestSuite {
public:
void test_load_emu_algorithm_init() {
- LoadEMU algToBeTested;
+ LoadEMUTar algToBeTested;
TS_ASSERT_THROWS_NOTHING(algToBeTested.initialize());
TS_ASSERT(algToBeTested.isInitialized());
}
void test_load_emu_algorithm() {
- LoadEMU algToBeTested;
+ LoadEMUTar algToBeTested;
if (!algToBeTested.isInitialized())
algToBeTested.initialize();
@@ -72,7 +72,8 @@ public:
// test some data properties
auto logpm = [&run](std::string tag) {
- return run.getPropertyValueAsType<double>(tag);
+ return dynamic_cast<TimeSeriesProperty<double> *>(
+ run.getProperty(tag))->firstValue();
};
TS_ASSERT_DELTA(logpm("DopplerFrequency"), 9.974, 1.0e-3);
TS_ASSERT_DELTA(logpm("DopplerAmplitude"), 0.075, 1.0e-4);
diff --git a/Framework/DataObjects/inc/MantidDataObjects/EventList.h b/Framework/DataObjects/inc/MantidDataObjects/EventList.h
index 22e967eaa25..e6efdb71d92 100644
--- a/Framework/DataObjects/inc/MantidDataObjects/EventList.h
+++ b/Framework/DataObjects/inc/MantidDataObjects/EventList.h
@@ -244,6 +244,7 @@ public:
void addPulsetime(const double seconds) override;
void maskTof(const double tofMin, const double tofMax) override;
+ void maskCondition(const std::vector<bool> & mask) override;
void getTofs(std::vector<double> &tofs) const override;
double getTofMin() const override;
@@ -456,6 +457,10 @@ private:
static std::size_t maskTofHelper(std::vector<T> &events, const double tofMin,
const double tofMax);
template <class T>
+ static std::size_t maskConditionHelper(std::vector<T> &events,
+ const std::vector<bool> &mask);
+
+ template <class T>
static void getTofsHelper(const std::vector<T> &events,
std::vector<double> &tofs);
template <class T>
diff --git a/Framework/DataObjects/src/EventList.cpp b/Framework/DataObjects/src/EventList.cpp
index 4cef5659ab2..cb1ad8aff2b 100644
--- a/Framework/DataObjects/src/EventList.cpp
+++ b/Framework/DataObjects/src/EventList.cpp
@@ -2695,6 +2695,92 @@ void EventList::maskTof(const double tofMin, const double tofMax) {
this->clear(false);
}
+// --------------------------------------------------------------------------
+/** Mask out events by the condition vector.
+* Events are removed from the list.
+* @param events :: reference to a vector of events to change.
+* @param mask :: condition vector
+* @returns The number of events deleted.
+*/
+template <class T>
+std::size_t EventList::maskConditionHelper(std::vector<T> &events,
+ const std::vector<bool> &mask) {
+
+ // runs through the two synchronized vectors and delete elements
+ // for condition false
+ auto itm = std::find(mask.begin(), mask.end(), false);
+ auto first = events.begin() + (itm - mask.begin());
+
+ if (itm != mask.end()) {
+ for (auto ite = first; ++ite != events.end() && ++itm != mask.end();) {
+ if (*itm != false) {
+ *first++ = std::move(*ite);
+ }
+ }
+ }
+
+ auto n = events.end() - first;
+ if (n != 0)
+ events.erase(first, events.end());
+
+ return n;
+
+ /*
+ size_t count = 0;
+ auto it_e = events.begin();
+ auto it_m = mask.begin();
+ while (it_e != events.end() && it_m != mask.end()) {
+ if (*it_m) {
+ ++it_e;
+ }
+ else {
+ it_e = events.erase(it_e);
+ count++;
+ }
+ ++it_m;
+ }
+ return count;
+ */
+}
+
+// --------------------------------------------------------------------------
+/**
+* Mask out events by the condition vector.
+* Events are removed from the list.
+* @param mask :: condition vector
+*/
+void EventList::maskCondition(const std::vector<bool> & mask) {
+
+ // mask size must match the number of events
+ if (this->getNumberEvents() != mask.size())
+ throw std::runtime_error("EventList::maskTof: tofMax must be > tofMin");
+
+ // don't do anything with an emply list
+ if (this->getNumberEvents() == 0)
+ return;
+
+ // Convert the list
+ size_t numOrig = 0;
+ size_t numDel = 0;
+ switch (eventType) {
+ case TOF:
+ numOrig = this->events.size();
+ numDel = this->maskConditionHelper(this->events, mask);
+ break;
+ case WEIGHTED:
+ numOrig = this->weightedEvents.size();
+ numDel = this->maskConditionHelper(this->weightedEvents, mask);
+ break;
+ case WEIGHTED_NOTIME:
+ numOrig = this->weightedEventsNoTime.size();
+ numDel = this->maskConditionHelper(this->weightedEventsNoTime, mask);
+ break;
+ }
+
+ if (numDel >= numOrig)
+ this->clear(false);
+}
+
// --------------------------------------------------------------------------
/** Get the m_tof member of all events in a list
*
diff --git a/Framework/DataObjects/test/EventListTest.h b/Framework/DataObjects/test/EventListTest.h
index d29c5ac56fd..ffc6ac593c1 100644
--- a/Framework/DataObjects/test/EventListTest.h
+++ b/Framework/DataObjects/test/EventListTest.h
@@ -1314,6 +1314,41 @@ public:
}
}
+ //-----------------------------------------------------------------------------------------------
+ void test_maskCondition_allTypes() {
+ // Go through each possible EventType as the input
+ for (int this_type = 0; this_type < 3; this_type++) {
+ this->fake_uniform_data();
+ el.switchTo(static_cast<EventType>(this_type));
+
+ // tof steps of 5 microseconds, starting at 100 ns, up to 20 msec
+ // How many events did we make?
+ TS_ASSERT_EQUALS(el.getNumberEvents(), 2 * MAX_TOF / BIN_DELTA);
+
+ // Mask out 5-10 milliseconds
+ auto nlen = el.getNumberEvents();
+ std::vector<bool> mask(nlen, true);
+
+ // first check no removal
+ el.maskCondition(mask);
+ TS_ASSERT_EQUALS(el.getNumberEvents(), 2 * MAX_TOF / BIN_DELTA);
+
+ double min = MAX_TOF * 0.25;
+ double max = MAX_TOF * 0.5;
+ for (size_t i = 0; i < nlen; i++) {
+ if ((el.getEvent(i).tof() >= min) && (el.getEvent(i).tof() <= max)) {
+ mask[i] = false;
+ }
+ }
+ el.maskCondition(mask);
+ for (std::size_t i = 0; i < el.getNumberEvents(); i++) {
+ // No tofs in that range
+ TS_ASSERT((el.getEvent(i).tof() < min) || (el.getEvent(i).tof() > max));
+ }
+ TS_ASSERT_EQUALS(el.getNumberEvents(), 0.75 * 2 * MAX_TOF / BIN_DELTA);
+ }
+ }
+
//-----------------------------------------------------------------------------------------------
void test_getTofs_and_setTofs() {
// Go through each possible EventType as the input
diff --git a/Framework/PythonInterface/mantid/api/src/Exports/IEventList.cpp b/Framework/PythonInterface/mantid/api/src/Exports/IEventList.cpp
index 006e403de80..92581c82e27 100644
--- a/Framework/PythonInterface/mantid/api/src/Exports/IEventList.cpp
+++ b/Framework/PythonInterface/mantid/api/src/Exports/IEventList.cpp
@@ -7,6 +7,7 @@
#include "MantidAPI/IEventList.h"
#include "MantidPythonInterface/kernel/GetPointer.h"
#include "MantidPythonInterface/kernel/Policies/VectorToNumpy.h"
+#include "MantidPythonInterface/kernel/Converters/NDArrayToVector.h"
#include <boost/python/class.hpp>
#include <boost/python/enum.hpp>
#include <boost/python/register_ptr_to_python.hpp>
@@ -18,11 +19,16 @@ using Mantid::API::TOF;
using Mantid::API::WEIGHTED;
using Mantid::API::WEIGHTED_NOTIME;
-namespace Policies = Mantid::PythonInterface::Policies;
+using namespace Mantid::PythonInterface;
using namespace boost::python;
GET_POINTER_SPECIALIZATION(IEventList)
+namespace {
+void maskCondition(IEventList &self, const NDArray &data) {
+ self.maskCondition(Converters::NDArrayToVector<bool>(data)());
+}
+}
/// return_value_policy for copied numpy array
using return_clone_numpy = return_value_policy<Policies::VectorToNumpy>;
@@ -67,6 +73,8 @@ void export_IEventList() {
.def("maskTof", &IEventList::maskTof, args("self", "tofMin", "tofMax"),
"Mask out events that have a tof between tofMin and tofMax "
"(inclusively)")
+ .def("maskCondition", &maskCondition, args("self", "mask"),
+ "Mask out events by the condition vector")
.def("getTofs",
(std::vector<double>(IEventList::*)(void) const) &
IEventList::getTofs,
@@ -84,6 +92,10 @@ void export_IEventList() {
"Get a vector of the weights of the events")
.def("getPulseTimes", &IEventList::getPulseTimes, args("self"),
"Get a vector of the pulse times of the events")
+ .def("getPulseTimeMax", &IEventList::getPulseTimeMax, args("self"),
+ "The maximum pulse time for the list of the events.")
+ .def("getPulseTimeMin", &IEventList::getPulseTimeMin, args("self"),
+ "The minimum pulse time for the list of the events.")
.def("getTofMin", &IEventList::getTofMin, args("self"),
"The minimum tof value for the list of the events.")
.def("getTofMax", &IEventList::getTofMax, args("self"),
diff --git a/Framework/PythonInterface/mantid/api/src/Exports/IEventWorkspace.cpp b/Framework/PythonInterface/mantid/api/src/Exports/IEventWorkspace.cpp
index 291268be925..94a7ded037a 100644
--- a/Framework/PythonInterface/mantid/api/src/Exports/IEventWorkspace.cpp
+++ b/Framework/PythonInterface/mantid/api/src/Exports/IEventWorkspace.cpp
@@ -45,6 +45,12 @@ void export_IEventWorkspace() {
.def("getTofMax", &IEventWorkspace::getTofMax, args("self"),
"Returns the maximum TOF value (in microseconds) held by the "
":class:`~mantid.api.Workspace`")
+ .def("getPulseTimeMin", &IEventWorkspace::getPulseTimeMin, args("self"),
+ "Returns the minimum pulse time held by the "
+ ":class:`~mantid.api.Workspace`")
+ .def("getPulseTimeMax", &IEventWorkspace::getPulseTimeMax, args("self"),
+ "Returns the maximum pulse time held by the "
+ ":class:`~mantid.api.Workspace`")
.def("getEventList", &deprecatedGetEventList,
return_internal_reference<>(), args("self", "workspace_index"),
"Return the :class:`~mantid.api.IEventList` managing the events at "
--
GitLab