24#include <boost/filesystem.hpp>
33using namespace Kernel;
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;
53constexpr char MaskStr[] =
"Mask";
54constexpr char SelectDetectorTubesStr[] =
"SelectDetectorTubes";
55constexpr char SelectDatasetStr[] =
"SelectDataset";
58constexpr char PathToBinaryStr[] =
"BinaryEventPath";
59constexpr char TOFBiasStr[] =
"TimeOfFlightBias";
60constexpr char CalibrateTOFStr[] =
"CalibrateTOFBias";
61constexpr char LambdaOnTwoStr[] =
"LambdaOnTwoMode";
64using TimeLimits = std::pair<double, double>;
66template <
typename Type>
70 auto p =
new Kernel::TimeSeriesProperty<Type>(name);
71 p->addValue(time,
value);
74 logManager.addProperty(p);
79template <
typename Type>
80Type GetNeXusValue(
const NeXus::NXEntry &entry,
const std::string &path,
const Type &defval, int32_t
index) {
82 NeXus::NXDataSetTyped<Type> dataSet = entry.openNXDataSet<Type>(path);
85 return dataSet()[
index];
86 }
catch (std::runtime_error &) {
93double GetNeXusValue<double>(
const NeXus::NXEntry &entry,
const std::string &path,
const double &defval,
96 NeXus::NXDataSetTyped<float> dataSet = entry.openNXDataSet<
float>(path);
99 return dataSet()[
index];
100 }
catch (std::runtime_error &) {
106std::string GetNeXusValue<std::string>(
const NeXus::NXEntry &entry,
const std::string &path,
const std::string &defval,
113 return std::string(dataSet(), dataSet.dim0());
114 }
catch (std::runtime_error &) {
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) {
123 T
value = GetNeXusValue<T>(entry, path, defval,
index);
124 logManager.addProperty<T>(name,
value * factor);
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 & , int32_t
index) {
133 std::string
value = GetNeXusValue<std::string>(entry, path, defval,
index);
134 logManager.addProperty<std::string>(name,
value);
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,
142 auto value = GetNeXusValue<T>(entry, path, defval,
index);
143 AddSinglePointTimeSeriesProperty<T>(logManager, time, name,
value * factor);
149template <
typename T,
typename F>
void mapRangeToIndex(
const std::string &line, std::vector<T> &result,
const F &fn) {
151 std::stringstream ss(line);
154 while (std::getline(ss, item,
',')) {
155 auto const k = item.find(
'-');
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));
162 p0 = boost::lexical_cast<size_t>(item);
166 if (p1 < result.size() && p0 <= p1) {
168 result[p0++] = fn(
index);
171 }
else if (p0 < result.size() && p1 < p0) {
173 result[p0] = fn(
index);
177 throw std::invalid_argument(
"invalid range specification");
187 explicit FileLoader(
const char *filename) :
_ifs(filename,
std::ios::binary |
std::ios::in) {
189 throw std::runtime_error(
"unable to open file");
196 bool read(
char *s, std::streamsize
n) {
return static_cast<bool>(
_ifs.read(s,
n)); }
198 size_t size()
const {
return _size; }
202 size_t selected_position() {
return _ifs.tellg(); }
223 m_M = (
static_cast<double>(N) / (maxVal - minVal));
229 inline double ival(
double val)
const {
return m_M * val +
m_B; }
231 inline double xval(
double ix)
const {
return (ix -
m_B) /
m_M; }
233 inline void add(
double val) {
234 auto ix =
static_cast<size_t>(std::floor(
ival(val)));
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)
299 return static_cast<int64_t
>(start * 1.0e3);
308 void addEvent(
size_t x,
size_t p,
double tof,
double ) {
315 auto y =
static_cast<size_t>(p / HISTO_BINS_Y_DENUMERATOR);
324 size_t id = xid < DETECTOR_TUBES ? PIXELS_PER_TUBE * xid +
y : DETECTOR_SPECTRA + xid;
325 if (
id >=
m_roi.size())
354 const std::vector<double> &
m_L2;
360 double deltaT = 1.0e6 * (
m_L1 +
m_L2[id]) /
m_V0 - tobs;
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)
380 for (
size_t i = 0; i < hvec.size(); i++) {
381 if (hvec[i] >= minLevel) {
382 auto ix =
static_cast<double>(i);
388 return (
count > 0 ? sum /
static_cast<double>(
count) : 0.0);
418 auto ev = Types::Event::TofEvent(tof, Types::Core::DateAndTime(offset));
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)
434template <
typename EP>
437 using namespace ANSTO;
441 FileLoader loader(eventFile.c_str());
446 ReadEventFile(loader, eventProcessor, progTracker, 100,
false);
455 std::vector<std::string> exts;
460 exts.emplace_back(
".hdf");
462 "The input filename of the stored data");
465 "Relative or absolute path to the compressed binary\n"
466 "event file linked to the HDF file, eg /storage/data/");
470 exts.emplace_back(
".xml");
472 "The input filename of the mask data");
475 "Comma separated range of detectors tubes to be loaded,\n"
481 declareProperty(SelectDatasetStr, 0,
"Select the index for the dataset to be loaded.");
483 declareProperty(TOFBiasStr, 0.0,
"Time of flight correction in micro-sec.");
485 declareProperty(CalibrateTOFStr,
false,
"Calibrate the TOF correction from the elastic pulse.");
487 declareProperty(LambdaOnTwoStr,
false,
"Instrument is operating in Lambda on Two mode.");
490 "Only include events after the provided start time, in "
491 "seconds (relative to the start of the run).");
494 "Only include events before the provided stop time, in "
495 "seconds (relative to the start of the run).");
497 std::string grpOptional =
"Filters";
525void LoadPLN::exec(
const std::string &hdfFile,
const std::string &eventFile) {
527 namespace fs = boost::filesystem;
531 fs::path p = hdfFile;
532 for (; !p.extension().empty();)
534 std::string title = p.generic_string();
542 prog.
doReport(
"creating instrument");
548 logManager.
addProperty(SelectDetectorTubesStr, seltubes);
554 timeMaxBoundary = std::numeric_limits<double>::infinity();
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; });
566 std::vector<EventVector_pt> eventVectors(numberHistograms,
nullptr);
567 std::vector<size_t> eventCounts(numberHistograms, 0);
571 double framePeriod = 1.0e6 / masterRpm;
575 double gatePeriod = (std::round(masterRpm / slaveRpm) == 1.0 ? 0.5 * framePeriod : framePeriod);
576 AddSinglePointTimeSeriesProperty<double>(logManager,
m_startRun,
"GatePeriod", gatePeriod);
579 size_t hdfCounts =
static_cast<size_t>(logManager.
getTimeSeriesProperty<
int>(
"TotalCounts")->firstValue());
581 double sourceSample =
fabs(instr->getSource()->getPos().Z());
584 double sampleTime = 1.0e6 * sourceSample / velocity;
585 PLN::EventCounter eventCounter(roi, detMapIndex, framePeriod, gatePeriod, timeBoundary, eventCounts, sourceSample,
587 PLN::loadEvents(prog,
"loading neutron counts", eventFile, eventCounter);
601 Types::Core::DateAndTime startTime(
m_startRun);
602 auto const start_nanosec = startTime.totalNanoseconds();
603 bool const calibrateTOF =
getProperty(CalibrateTOFStr);
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);
616 auto minTOF = eventAssigner.
tofMin();
617 auto maxTOF = eventAssigner.
tofMax();
620 m_localWorkspace->setAllX(HistogramData::BinEdges{std::max(0.0, floor(minTOF)), maxTOF + 1});
624 auto frame_count =
static_cast<int>(eventCounter.
numFrames());
625 AddSinglePointTimeSeriesProperty<int>(logManager,
m_startRun,
"frame_count", frame_count);
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());
648 auto detIDs = detectorInfo.detectorIDs();
649 for (
const auto detID : detIDs) {
650 auto ix = detectorInfo.indexOf(detID);
651 double l2 = detectorInfo.l2(ix);
660 size_t maskedBins = 0;
661 for (
size_t i = 0; i != roi.size(); i++)
665 if (maskedBins > 0) {
667 std::vector<size_t> maskIndexList(maskedBins);
668 size_t maskIndex = 0;
670 for (
size_t i = 0; i != roi.size(); i++)
672 maskIndexList[maskIndex++] = i;
676 maskingAlg->setProperty(
"WorkspaceIndexList", maskIndexList);
677 maskingAlg->executeAsChildAlg();
684 std::vector<EventVector_pt> &eventVectors) {
686 size_t numberHistograms = eventCounts.size();
687 for (
size_t i = 0; i != numberHistograms; ++i) {
691 eventList.
reserve(eventCounts[i]);
707 std::vector<bool> result(HISTOGRAMS,
true);
710 if (!selected.empty()) {
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;
720 for (
size_t i = 0; i < MONITORS; i++) {
721 result[DETECTOR_SPECTRA + i] = tubes[DETECTOR_TUBES + i];
725 if (maskfile.length() == 0)
728 std::ifstream input(maskfile.c_str());
730 throw std::invalid_argument(
"invalid mask file");
733 while (std::getline(input, line)) {
734 auto i0 = line.find(
"<detids>");
735 auto iN = line.find(
"</detids>");
737 if ((i0 != std::string::npos) && (iN != std::string::npos) && (i0 < iN)) {
738 line = line.substr(i0 + 8, iN - i0 - 8);
739 mapRangeToIndex(line, result, [](
size_t) {
return false; });
752 MapNeXusToProperty<std::string>(entry,
"sample/name",
"unknown", logm,
"SampleName",
"", 0);
753 MapNeXusToProperty<std::string>(entry,
"sample/description",
"unknown", logm,
"SampleDescription",
"", 0);
756 Types::Core::DateAndTime startTime(GetNeXusValue<std::string>(entry,
"start_time",
"2000-01-01T00:00:00", 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);
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);
773 bool const lambdaOnTwoMode =
getProperty(LambdaOnTwoStr);
774 double lambdaFactor = (lambdaOnTwoMode ? 0.5 : 1.0);
775 logm.
addProperty(
"LambdaOnTwoMode", (lambdaOnTwoMode ? 1 : 0));
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",
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,
800 std::vector<std::string> tags = {
"P01PS05",
"P01PSP05",
"T01S00",
"T01S06",
"T01S07",
"T01S08",
"T01SP00",
"T01SP06",
801 "T01SP07",
"T01SP08",
"T2S1",
"T2S2",
"T2S3",
"T2S4",
"T2SP1",
"T2SP2"};
803 for (
const auto &tag : tags) {
804 MapNeXusToSeries<double>(entry,
"data/" + tag, 0.0, logm, time_str,
"env_" + tag, 1.0,
m_datasetIndex);
814 loadInstrumentAlg->setPropertyValue(
"InstrumentName",
"PELICAN");
816 loadInstrumentAlg->executeAsChildAlg();
831const std::string
LoadPLN::summary()
const {
return "Loads a PLN Hdf and linked event file into a workspace."; }
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")) {
857 namespace fs = boost::filesystem;
866 if (evtPath.rfind(
"./") == 0 || evtPath.rfind(
"../") == 0) {
867 fs::path hp = hdfFile;
868 evtPath = fs::canonical(evtPath, hp.parent_path()).generic_string();
877 if (fs::is_directory(evtPath)) {
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);
883 std::string message(
"Negative dataset index recorded in HDF, reset to zero!");
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)) {
900 fs::path hp = hdfFile;
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();
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);
917 exec(hdfFile, evtPath);
double value
The value of the point.
std::map< DeltaEMode::Type, std::string > index
#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.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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.
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
@ Load
allowed here which will be passed to the algorithm
void setDetectorID(const detid_t detID)
Clear the list of detector IDs, then add one.
void setSpectrumNo(specnum_t num)
Sets the the spectrum number of this spectrum.
This class contains the information about the log entries.
void addProperty(Kernel::Property *prop, bool overwrite=false)
Add data to the object in the form of a property.
HeldType getPropertyValueAsType(const std::string &name) const
Get the value of a property as the given TYPE.
Kernel::TimeSeriesProperty< T > * getTimeSeriesProperty(const std::string &name) const
Returns a property as a time series property.
Helper class for reporting progress from algorithms.
void doReport(const std::string &msg="") override
Actually do the reporting, without changing the loop counter.
A property class for workspaces.
helper class to keep track of progress
void update(int64_t position)
LoadPLN : Loads an ANSTO PLN Hdf and linked event file into a workspace.
std::vector< double > m_detectorL2
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
const std::vector< std::string > seeAlso() const override
Similar algorithms.
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.
int version() const override
Algorithm's version for identification.
void setupDetectorMasks(const std::vector< bool > &roi)
Set up the detector masks to the region of interest roi.
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,...
DataObjects::EventWorkspace_sptr m_localWorkspace
void init() override
Initialise the algorithm and declare the properties for the nexus descriptor.
void createWorkspace(const std::string &title)
Creates an event workspace and sets the title.
void loadParameters(const std::string &hdfFile, API::LogManager &logm)
Load parameters from input hdfFile and save to the log manager, logm.
void loadDetectorL2Values()
Recovers the L2 neutronic distance for each detector.
void loadInstrument()
Load the instrument definition.
std::vector< bool > createRoiVector(const std::string &seltubes, const std::string &maskfile)
Region of interest is defined by the selected detectors and the maskfile.
const std::string name() const override
Algorithms name for identification.
void exec() override
Execute the algorithm.
int confidence(Kernel::NexusDescriptor &descriptor) const override
Return the confidence as an integer value that this algorithm can load the file descriptor.
const std::string category() const override
Algorithm's category for identification.
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)
std::vector< EventVector_pt > & m_eventVectors
void addEventImpl(size_t id, size_t, size_t, double tobs) override
std::vector< size_t > & m_eventCounts
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)
void addEventImpl(size_t id, size_t, size_t, double tobs) override
const std::vector< double > & m_L2
const std::vector< bool > & m_roi
void addEvent(size_t x, size_t p, double tof, double)
const std::vector< size_t > & m_mapIndex
size_t availableEvents() const
const double m_framePeriod
const TimeLimits m_timeBoundary
virtual void addEventImpl(size_t id, size_t x, size_t y, double tof)=0
int64_t frameStart() const
const double m_gatePeriod
size_t processedEvents() const
EventProcessor(const std::vector< bool > &roi, const std::vector< size_t > &mapIndex, const double framePeriod, const double gatePeriod, const TimeLimits &timeBoundary, size_t maxEvents)
double ival(double val) const
std::vector< size_t > m_hist
double xval(double ix) const
SimpleHist(size_t N, double minVal, double maxVal)
const std::vector< size_t > & histogram() const
void setSortOrder(const EventSortType order) const
Manually set the event list sort order value.
void reserve(size_t num) override
Reserve a certain number of entries in event list of the specified eventType.
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.
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.
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.
Implements NXroot Nexus class.
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)
static char const *const FilenameStr
static const size_t HISTO_BINS_X
static const size_t Progress_LoadBinFile
void AddSinglePointTimeSeriesProperty(API::LogManager &logManager, const std::string &time, const std::string &name, const TYPE value)
static char const *const MaskStr
static char const *const FilterByTimeStartStr
static const size_t Progress_Total
static char const *const FilterByTimeStopStr
static const size_t Progress_ReserveMemory
DLLExport void getEventsFrom(EventList &el, std::vector< Types::Event::TofEvent > *&events)
NXDataSetTyped< char > NXChar
The char dataset type.
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.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Output
An output workspace.