26#include <boost/algorithm/string.hpp>
27#include <boost/algorithm/string/trim.hpp>
28#include <boost/math/special_functions/round.hpp>
30#include <Poco/AutoPtr.h>
31#include <Poco/DOM/AutoPtr.h>
32#include <Poco/DOM/DOMParser.h>
33#include <Poco/DOM/Document.h>
34#include <Poco/DOM/Element.h>
35#include <Poco/DOM/NodeFilter.h>
36#include <Poco/DOM/NodeIterator.h>
37#include <Poco/DOM/NodeList.h>
38#include <Poco/TemporaryFile.h>
39#include <Poco/Util/PropertyFileConfiguration.h>
63using ANSTO::EventVector_pt;
65template <typename
TYPE>
70 p->addValue(time,
value);
73 logManager.addProperty(p);
92 const std::vector<std::string> &subFiles = file.
files();
93 for (
const auto &subFile : subFiles) {
94 auto len = subFile.length();
95 if ((len > 4) && (subFile.find_first_of(
"\\/", 0, 2) == std::string::npos)) {
96 if ((subFile.rfind(
".hdf") == len - 4) && (subFile.compare(0, 3,
"BBY") == 0))
98 else if (subFile.rfind(
".bin") == len - 4)
103 return (hdfFiles == 1) && (binFiles == 1) ? 50 : 0;
112 std::vector<std::string> exts;
117 exts.emplace_back(
".tar");
119 "The input filename of the stored data");
123 exts.emplace_back(
".xml");
125 "The input filename of the mask data");
133 "Optional: To exclude events that do not fall within a range "
134 "of times-of-flight. "
135 "This is the minimum accepted value in microseconds. Keep "
136 "blank to load all events.");
141 "Optional: To exclude events that do not fall within a range "
142 "of times-of-flight. "
143 "This is the maximum accepted value in microseconds. Keep "
144 "blank to load all events.");
149 "Optional: To only include events after the provided start time, in "
150 "seconds (relative to the start of the run).");
155 "Optional: To only include events before the provided stop time, in "
156 "seconds (relative to the start of the run).");
158 std::string grpOptional =
"Filters";
177 throw std::invalid_argument(
"invalid BBY file");
189 tofMaxBoundary = std::numeric_limits<double>::infinity();
191 timeMaxBoundary = std::numeric_limits<double>::infinity();
194 prog.
doReport(
"creating instrument");
205 std::map<std::string, double> logParams;
206 std::map<std::string, std::string> logStrings;
207 std::map<std::string, std::string> allParams;
208 createInstrument(tarFile, instrumentInfo, logParams, logStrings, allParams);
211 if (instrumentInfo.
is_tof)
216 eventWS->setYUnit(
"Counts");
219 const std::vector<std::string> &subFiles = tarFile.
files();
220 for (
const auto &subFile : subFiles)
221 if (subFile.compare(0, 3,
"BBY") == 0) {
222 std::string title = subFile;
224 if (title.rfind(
".hdf") == title.length() - 4)
225 title.resize(title.length() - 4);
227 if (title.rfind(
".nx") == title.length() - 3)
228 title.resize(title.length() - 3);
230 eventWS->setTitle(title);
235 size_t numberHistograms = eventWS->getNumberHistograms();
237 std::vector<EventVector_pt> eventVectors(numberHistograms,
nullptr);
238 std::vector<size_t> eventCounts(numberHistograms, 0);
246 double period = periodSlave;
247 double shift = -1.0 / 6.0 * periodMaster - periodSlave * phaseSlave / 360.0;
250 Types::Core::DateAndTime startTime(instrumentInfo.
start_time);
251 auto startInNanosec = startTime.totalNanoseconds();
255 timeMinBoundary, timeMaxBoundary, eventCounts);
257 loadEvents(prog,
"loading neutron counts", tarFile, eventCounter);
262 for (
size_t i = 0; i != numberHistograms; ++i) {
266 eventList.
reserve(eventCounts[i]);
277 if (instrumentInfo.
is_tof) {
279 timeMinBoundary, timeMaxBoundary, eventVectors);
281 loadEvents(prog,
"loading neutron events (TOF)", tarFile, eventAssigner);
284 startInNanosec, tofMinBoundary, tofMaxBoundary, timeMinBoundary,
285 timeMaxBoundary, eventVectors);
287 loadEvents(prog,
"loading neutron events (Wavelength)", tarFile, eventAssigner);
290 auto getParam = [&allParams](
const std::string &tag,
double defValue) {
292 return std::stod(allParams[tag]);
293 }
catch (
const std::invalid_argument &) {
297 if (instrumentInfo.
is_tof) {
299 eventWS->setAllX(HistogramData::BinEdges{std::max(0.0, floor(eventCounter.
tofMin())), eventCounter.
tofMax() + 1});
301 double lof = getParam(
"wavelength_extn_lo", 0.95);
302 double hif = getParam(
"wavelength_extn_hi", 1.05);
303 eventWS->setAllX(HistogramData::BinEdges{instrumentInfo.
wavelength * lof, instrumentInfo.
wavelength * hif});
307 size_t maskedBins = 0;
308 for (
size_t i = 0; i != roi.size(); i++)
312 if (maskedBins > 0) {
314 std::vector<size_t> maskIndexList(maskedBins);
315 size_t maskIndex = 0;
317 for (
size_t i = 0; i != roi.size(); i++)
319 maskIndexList[maskIndex++] = i;
322 maskingAlg->setProperty(
"Workspace", eventWS);
323 maskingAlg->setProperty(
"WorkspaceIndexList", maskIndexList);
324 maskingAlg->executeAsChildAlg();
330 auto frame_count =
static_cast<int>(eventCounter.
numFrames());
334 logManager.
addProperty(
"frame_count", frame_count);
339 logManager.
addProperty(
"bm_counts",
static_cast<double>(frame_count) * period /
342 Types::Core::time_duration duration =
343 boost::posix_time::microseconds(
static_cast<boost::int64_t
>(
static_cast<double>(frame_count) * period));
345 Types::Core::DateAndTime start_time(instrumentInfo.
start_time);
346 Types::Core::DateAndTime end_time(start_time + duration);
348 logManager.
addProperty(
"start_time", start_time.toISO8601String());
349 logManager.
addProperty(
"run_start", start_time.toISO8601String());
350 logManager.
addProperty(
"end_time", end_time.toISO8601String());
353 std::string time_str = start_time.toISO8601String();
361 for (
auto &
x : logStrings) {
364 for (
auto &
x : logParams) {
369 loadInstrumentAlg->setProperty(
"Workspace", eventWS);
370 loadInstrumentAlg->setPropertyValue(
"InstrumentName",
"BILBY");
372 loadInstrumentAlg->executeAsChildAlg();
381 if (maskfile.length() == 0)
384 std::ifstream input(maskfile.c_str());
386 throw std::invalid_argument(
"invalid mask file");
389 while (std::getline(input, line)) {
390 auto i0 = line.find(
"<detids>");
391 auto iN = line.find(
"</detids>");
393 if ((i0 != std::string::npos) && (iN != std::string::npos) && (i0 < iN)) {
394 line = line.substr(i0 + 8, iN - i0 - 8);
395 std::stringstream ss(line);
398 while (std::getline(ss, item,
',')) {
399 auto k = item.find(
'-');
402 if (k != std::string::npos) {
403 p0 = boost::lexical_cast<size_t>(item.substr(0, k));
404 p1 = boost::lexical_cast<size_t>(item.substr(k + 1, item.size() - k - 1));
409 p0 = boost::lexical_cast<size_t>(item);
413 if (p0 < result.size()) {
414 if (p1 >= result.size())
415 p1 = result.size() - 1;
418 result[p0++] =
false;
429 std::map<std::string, std::string> &logStrings,
430 std::map<std::string, std::string> &allParams) {
435 std::string parameterFilename = idfDirectory +
"BILBY_Parameters.xml";
438 Poco::AutoPtr<Poco::XML::Document> pDoc;
440 pDoc = pParser.parse(parameterFilename);
444 NodeIterator it(pDoc, Poco::XML::NodeFilter::SHOW_ELEMENT);
445 Node *pNode = it.nextNode();
447 if (pNode->nodeName() ==
"parameter") {
448 auto pElem =
dynamic_cast<Element *
>(pNode);
449 std::string
name = pElem->getAttribute(
"name");
450 Poco::AutoPtr<NodeList> nodeList = pElem->childNodes();
451 for (
unsigned long i = 0; i < nodeList->length(); i++) {
452 auto cNode = nodeList->item(i);
453 if (cNode->nodeName() ==
"value") {
454 auto cElem =
dynamic_cast<Poco::XML::Element *
>(cNode);
455 std::string
value = cElem->getAttribute(
"val");
460 pNode = it.nextNode();
463 auto isNumeric = [](
const std::string &tag) {
465 auto stag = boost::algorithm::trim_copy(tag);
467 auto value = std::stod(stag, &sz);
468 return sz > 0 && stag.size() == sz && std::isfinite(
value);
469 }
catch (
const std::invalid_argument &) {
474 std::string tmpString;
475 float tmpFloat = 0.0f;
476 for (
auto &
x : allParams) {
477 if (
x.first.find(
"log_") == 0) {
478 auto logTag = boost::algorithm::trim_copy(
x.first.substr(4));
479 auto line =
x.second;
482 std::vector<std::string> details;
483 boost::split(details, line, boost::is_any_of(
","));
484 if (details.size() < 3) {
485 g_log.
warning() <<
"Invalid format for BILBY parameter " <<
x.first << std::endl;
488 auto hdfTag = boost::algorithm::trim_copy(details[0]);
492 auto updateOk =
false;
493 if (!hdfTag.empty()) {
494 if (isNumeric(details[1])) {
496 auto factor = std::stod(details[1]);
497 logParams[logTag] = factor * tmpFloat;
501 logStrings[logTag] = tmpString;
508 auto defValue = boost::algorithm::trim_copy(details[2]);
509 if (!defValue.empty()) {
510 if (isNumeric(defValue))
511 logParams[logTag] = std::stod(defValue);
513 logStrings[logTag] = defValue;
515 g_log.
warning() <<
"Cannot find hdf parameter " << hdfTag <<
", using default.\n";
518 }
catch (
const std::invalid_argument &) {
519 g_log.
warning() <<
"Invalid format for BILBY parameter " <<
x.first << std::endl;
523 }
catch (std::exception &ex) {
524 g_log.
warning() <<
"Failed to load instrument with error: " << ex.what()
525 <<
". The current facility may not be fully "
532 std::map<std::string, double> &logParams, std::map<std::string, std::string> &logStrings,
533 std::map<std::string, std::string> &allParams) {
535 const double toMeters = 1.0 / 1000;
539 instrumentInfo.
start_time =
"2000-01-01T00:00:00";
546 instrumentInfo.
is_tof =
true;
554 const std::vector<std::string> &files = tarFile.
files();
555 auto file_it = std::find_if(files.cbegin(), files.cend(),
556 [](
const std::string &file) { return file.rfind(
".hdf") == file.length() - 4; });
557 if (file_it != files.end()) {
558 tarFile.
select(file_it->c_str());
560 Poco::TemporaryFile hdfFile;
561 std::shared_ptr<FILE> handle(fopen(hdfFile.path().c_str(),
"wb"), fclose);
566 while (0 != (bytesRead = tarFile.
read(buffer,
sizeof(buffer))))
567 fwrite(buffer, bytesRead, 1, handle.get());
573 float tmp_float = 0.0f;
574 int32_t tmp_int32 = 0;
580 instrumentInfo.
att_pos = boost::math::iround(tmp_float);
589 if (
loadNXDataSet(entry,
"instrument/master1_chopper_id", tmp_int32))
591 if (
loadNXDataSet(entry,
"instrument/master2_chopper_id", tmp_int32))
594 if (
loadNXString(entry,
"instrument/detector/frame_source", tmp_str))
595 instrumentInfo.
is_tof = tmp_str ==
"EXTERNAL";
596 if (
loadNXDataSet(entry,
"instrument/nvs067/lambda", tmp_float))
599 if (
loadNXDataSet(entry,
"instrument/master_chopper_freq", tmp_float) && (tmp_float > 0.0f))
601 if (
loadNXDataSet(entry,
"instrument/t0_chopper_freq", tmp_float) && (tmp_float > 0.0f))
603 if (
loadNXDataSet(entry,
"instrument/t0_chopper_phase", tmp_float))
604 instrumentInfo.
phase_slave = tmp_float < 999.0 ? tmp_float : 0.0;
610 auto findLtof = logParams.find(
"Ltof_det_value");
611 if (findLtof != logParams.end()) {
612 logParams[
"L1_chopper_value"] = logParams[
"Ltof_det_value"] - logParams[
"L2_det_value"];
614 logParams[
"L1_chopper_value"] = 18.4726;
615 g_log.
warning() <<
"Cannot recover parameter 'L1_chopper_value'"
616 <<
", using default.\n";
622 file_it = std::find(files.cbegin(), files.cend(),
"History.log");
623 if (file_it != files.cend()) {
624 tarFile.
select(file_it->c_str());
625 std::string logContent;
627 tarFile.
read(&logContent[0], logContent.size());
628 std::istringstream data(logContent);
629 Poco::AutoPtr<Poco::Util::PropertyFileConfiguration> conf(
new Poco::Util::PropertyFileConfiguration(data));
631 if (conf->hasProperty(
"Bm1Counts"))
632 instrumentInfo.
bm_counts = conf->getInt(
"Bm1Counts");
633 if (conf->hasProperty(
"AttPos"))
634 instrumentInfo.
att_pos = boost::math::iround(conf->getDouble(
"AttPos"));
636 if (conf->hasProperty(
"SampleName"))
637 instrumentInfo.
sample_name = conf->getString(
"SampleName");
639 if (conf->hasProperty(
"MasterChopperFreq")) {
640 auto tmp = conf->getDouble(
"MasterChopperFreq");
644 if (conf->hasProperty(
"T0ChopperFreq")) {
645 auto tmp = conf->getDouble(
"T0ChopperFreq");
649 if (conf->hasProperty(
"T0ChopperPhase")) {
650 auto tmp = conf->getDouble(
"T0ChopperPhase");
654 if (conf->hasProperty(
"FrameSource"))
655 instrumentInfo.
is_tof = conf->getString(
"FrameSource") ==
"EXTERNAL";
656 if (conf->hasProperty(
"Wavelength"))
657 instrumentInfo.
wavelength = conf->getDouble(
"Wavelength");
659 if (conf->hasProperty(
"SampleAperture"))
660 logParams[
"sample_aperture"] = conf->getDouble(
"SampleAperture");
661 if (conf->hasProperty(
"SourceAperture"))
662 logParams[
"source_aperture"] = conf->getDouble(
"SourceAperture");
663 if (conf->hasProperty(
"L1"))
664 logParams[
"L1_source_value"] = conf->getDouble(
"L1") * toMeters;
665 if (conf->hasProperty(
"LTofDet"))
666 logParams[
"L1_chopper_value"] = conf->getDouble(
"LTofDet") * toMeters - logParams[
"L2_det_value"];
667 if (conf->hasProperty(
"L2Det"))
668 logParams[
"L2_det_value"] = conf->getDouble(
"L2Det") * toMeters;
670 if (conf->hasProperty(
"L2CurtainL"))
671 logParams[
"L2_curtainl_value"] = conf->getDouble(
"L2CurtainL") * toMeters;
672 if (conf->hasProperty(
"L2CurtainR"))
673 logParams[
"L2_curtainr_value"] = conf->getDouble(
"L2CurtainR") * toMeters;
674 if (conf->hasProperty(
"L2CurtainU"))
675 logParams[
"L2_curtainu_value"] = conf->getDouble(
"L2CurtainU") * toMeters;
676 if (conf->hasProperty(
"L2CurtainD"))
677 logParams[
"L2_curtaind_value"] = conf->getDouble(
"L2CurtainD") * toMeters;
679 if (conf->hasProperty(
"CurtainL"))
680 logParams[
"D_curtainl_value"] = conf->getDouble(
"CurtainL") * toMeters;
681 if (conf->hasProperty(
"CurtainR"))
682 logParams[
"D_curtainr_value"] = conf->getDouble(
"CurtainR") * toMeters;
683 if (conf->hasProperty(
"CurtainU"))
684 logParams[
"D_curtainu_value"] = conf->getDouble(
"CurtainU") * toMeters;
685 if (conf->hasProperty(
"CurtainD"))
686 logParams[
"D_curtaind_value"] = conf->getDouble(
"CurtainD") * toMeters;
698 }
catch (std::runtime_error &) {
707 value = std::string(dataSet(), dataSet.
dim0());
709 }
catch (std::runtime_error &) {
715template <
class EventProcessor>
717 EventProcessor &eventProcessor) {
721 int64_t fileSize = 0;
722 const std::vector<std::string> &files = tarFile.
files();
723 const auto found = std::find_if(files.cbegin(), files.cend(),
724 [](
const auto &file) { return file.rfind(
".bin") == file.length() - 4; });
725 if (found != files.cend()) {
726 tarFile.
select(found->c_str());
741 if ((fileSize == 0) || !tarFile.
skip(128))
745 int invalidEvents = 0;
747 while ((c =
static_cast<uint32_t
>(tarFile.
read_byte())) !=
static_cast<uint32_t
>(-1)) {
749 bool event_ended =
false;
756 x |= (c & 0x01) << 8;
761 event_ended = (c & 0xC0) != 0xC0;
765 y |= (c & 0x01) << 7;
766 dt = (c & 0xFE) >> 1;
774 event_ended = (c & 0xC0) != 0xC0;
779 dt |= (c & 0xFF) << (5 + 6 * (state - 3));
784 if (event_ended || (state == 8)) {
787 if ((
x == 0) && (
y == 0) && (dt == 0xFFFFFFFF)) {
789 eventProcessor.newFrame();
793 tof +=
static_cast<int>(dt) * 0.1;
797 tof +=
static_cast<int>(dt) * 0.1;
799 eventProcessor.addEvent(
x,
y, tof);
805 if (invalidEvents > 0) {
806 g_log.
warning() <<
"BILBY loader dropped " << invalidEvents <<
" invalid event(s)" << std::endl;
double value
The value of the point.
#define DECLARE_FILELOADER_ALGORITHM(classname)
DECLARE_FILELOADER_ALGORITHM should be used in place of the standard DECLARE_ALGORITHM macro when wri...
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.
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)
int64_t selected_position() const
size_t read(void *dst, size_t size)
bool select(const char *file)
const std::vector< std::string > & files() const
bool skip(uint64_t offset)
int64_t selected_size() const
static bool loadNXDataSet(const NeXus::NXEntry &entry, const std::string &path, T &value)
static std::vector< bool > createRoiVector(const std::string &maskfile)
void init() override
Initialise the algorithm.
void createInstrument(ANSTO::Tar::File &tarFile, InstrumentInfo &instrumentInfo, std::map< std::string, double > &logParams, std::map< std::string, std::string > &logStrings, std::map< std::string, std::string > &allParams)
const std::string name() const override
function to return a name of the algorithm, must be overridden in all algorithms
void loadEvents(API::Progress &prog, const char *progMsg, ANSTO::Tar::File &tarFile, EventProcessor &eventProcessor)
int confidence(Kernel::FileDescriptor &descriptor) const override
Return the confidence value that this algorithm can load the file.
void loadInstrumentParameters(const NeXus::NXEntry &entry, std::map< std::string, double > &logParams, std::map< std::string, std::string > &logStrings, std::map< std::string, std::string > &allParams)
void exec() override
Execute the algorithm.
bool loadNXString(const NeXus::NXEntry &entry, const std::string &path, std::string &value)
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.
Records the filename and the description of failure.
Defines a wrapper around an open file.
const std::string & filename() const
Access the filename.
const std::string & extension() const
Access the file extension.
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 warning(const std::string &msg)
Logs at warning level.
OptionalBool : Tri-state bool.
The concrete, templated class for properties.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
A specialised Property class for holding a series of time-value pairs.
NXDataSetTyped< T > openNXDataSet(const std::string &name) const
Templated method for creating datasets.
NXChar openNXChar(const std::string &name) const
Creates and opens a char dataset.
Templated class implementation of NXDataSet.
void load(const int blocksize=1, int i=-1, int j=-1, int k=-1, int l=-1) override
Implementation of the virtual NXDataSet::load(...) method.
int dim0() const
Returns the number of elements along the first dimension.
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
static char const *const FilenameStr
static char const *const FilterByTofMaxStr
static const size_t HISTO_BINS_X
static const size_t HISTO_BINS_Y
static const size_t Progress_LoadBinFile
static char const *const FilterByTofMinStr
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)
std::shared_ptr< EventWorkspace > EventWorkspace_sptr
shared pointer to the EventWorkspace class
int32_t detid_t
Typedef for a detector ID.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
int32_t master2_chopper_id
int32_t master1_chopper_id
std::string sample_description
@ Input
An input workspace.
@ Output
An output workspace.