27#include <boost/algorithm/string.hpp>
28#include <boost/algorithm/string/trim.hpp>
29#include <boost/math/special_functions/round.hpp>
31#include <Poco/AutoPtr.h>
32#include <Poco/DOM/AutoPtr.h>
33#include <Poco/DOM/DOMParser.h>
34#include <Poco/DOM/Document.h>
35#include <Poco/DOM/Element.h>
36#include <Poco/DOM/NodeFilter.h>
37#include <Poco/DOM/NodeIterator.h>
38#include <Poco/DOM/NodeList.h>
39#include <Poco/TemporaryFile.h>
40#include <Poco/Util/PropertyFileConfiguration.h>
66template <typename
TYPE>
71 p->addValue(time,
value);
74 logManager.addProperty(p);
93 const std::vector<std::string> &subFiles = file.
files();
94 for (
const auto &subFile : subFiles) {
95 auto len = subFile.length();
96 if ((len > 4) && (subFile.find_first_of(
"\\/", 0, 2) == std::string::npos)) {
97 if ((subFile.rfind(
".hdf") == len - 4) && (subFile.compare(0, 3,
"BBY") == 0))
99 else if (subFile.rfind(
".bin") == len - 4)
104 return (hdfFiles == 1) && (binFiles == 1) ? 50 : 0;
113 std::vector<std::string> exts;
118 exts.emplace_back(
".tar");
120 "The input filename of the stored data");
124 exts.emplace_back(
".xml");
126 "The input filename of the mask data");
134 "Optional: To exclude events that do not fall within a range "
135 "of times-of-flight. "
136 "This is the minimum accepted value in microseconds. Keep "
137 "blank to load all events.");
142 "Optional: To exclude events that do not fall within a range "
143 "of times-of-flight. "
144 "This is the maximum accepted value in microseconds. Keep "
145 "blank to load all events.");
150 "Optional: To only include events after the provided start time, in "
151 "seconds (relative to the start of the run).");
156 "Optional: To only include events before the provided stop time, in "
157 "seconds (relative to the start of the run).");
159 std::string grpOptional =
"Filters";
171 if (API::AnalysisDataService::Instance().doesExist(outName))
172 API::AnalysisDataService::Instance().remove(outName);
178 throw std::invalid_argument(
"invalid BBY file");
190 tofMaxBoundary = std::numeric_limits<double>::infinity();
192 timeMaxBoundary = std::numeric_limits<double>::infinity();
195 prog.
doReport(
"creating instrument");
206 std::map<std::string, double> logParams;
207 std::map<std::string, std::string> logStrings;
208 std::map<std::string, std::string> allParams;
209 createInstrument(tarFile, instrumentInfo, logParams, logStrings, allParams);
212 if (instrumentInfo.
is_tof)
213 eventWS->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create(
"TOF");
215 eventWS->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create(
"Wavelength");
217 eventWS->setYUnit(
"Counts");
220 const std::vector<std::string> &subFiles = tarFile.
files();
221 const auto it = std::find_if(subFiles.cbegin(), subFiles.cend(),
222 [](
const auto &subFile) { return subFile.compare(0, 3,
"BBY") == 0; });
223 if (it != subFiles.cend()) {
224 std::string title = *it;
226 if (title.rfind(
".hdf") == title.length() - 4)
227 title.resize(title.length() - 4);
229 if (title.rfind(
".nx") == title.length() - 3)
230 title.resize(title.length() - 3);
232 eventWS->setTitle(title);
236 size_t numberHistograms = eventWS->getNumberHistograms();
238 std::vector<EventVector_pt> eventVectors(numberHistograms,
nullptr);
239 std::vector<size_t> eventCounts(numberHistograms, 0);
247 double period = periodSlave;
248 double shift = -1.0 / 6.0 * periodMaster - periodSlave * phaseSlave / 360.0;
251 Types::Core::DateAndTime startTime(instrumentInfo.
start_time);
252 auto startInNanosec = startTime.totalNanoseconds();
256 timeMinBoundary, timeMaxBoundary, eventCounts);
258 loadEvents(prog,
"loading neutron counts", tarFile, eventCounter);
263 for (
size_t i = 0; i != numberHistograms; ++i) {
267 eventList.
reserve(eventCounts[i]);
278 if (instrumentInfo.
is_tof) {
280 timeMinBoundary, timeMaxBoundary, eventVectors);
282 loadEvents(prog,
"loading neutron events (TOF)", tarFile, eventAssigner);
285 startInNanosec, tofMinBoundary, tofMaxBoundary, timeMinBoundary,
286 timeMaxBoundary, eventVectors);
288 loadEvents(prog,
"loading neutron events (Wavelength)", tarFile, eventAssigner);
291 auto getParam = [&allParams](
const std::string &tag,
double defValue) {
293 return std::stod(allParams[tag]);
294 }
catch (
const std::invalid_argument &) {
298 if (instrumentInfo.
is_tof) {
300 eventWS->setAllX(HistogramData::BinEdges{std::max(0.0, floor(eventCounter.
tofMin())), eventCounter.
tofMax() + 1});
302 double lof = getParam(
"wavelength_extn_lo", 0.95);
303 double hif = getParam(
"wavelength_extn_hi", 1.05);
304 eventWS->setAllX(HistogramData::BinEdges{instrumentInfo.
wavelength * lof, instrumentInfo.
wavelength * hif});
308 size_t maskedBins = std::count_if(roi.cbegin(), roi.cend(), [](
bool v) { return !v; });
310 if (maskedBins > 0) {
312 std::vector<size_t> maskIndexList(maskedBins);
313 size_t maskIndex = 0;
315 for (
size_t i = 0; i != roi.size(); i++)
317 maskIndexList[maskIndex++] = i;
320 maskingAlg->setProperty(
"Workspace", eventWS);
321 maskingAlg->setProperty(
"WorkspaceIndexList", maskIndexList);
322 maskingAlg->executeAsChildAlg();
328 auto frame_count =
static_cast<int>(eventCounter.
numFrames());
332 logManager.
addProperty(
"frame_count", frame_count);
337 logManager.
addProperty(
"bm_counts",
static_cast<double>(frame_count) * period /
340 Types::Core::time_duration duration =
341 boost::posix_time::microseconds(
static_cast<boost::int64_t
>(
static_cast<double>(frame_count) * period));
343 Types::Core::DateAndTime start_time(instrumentInfo.
start_time);
344 Types::Core::DateAndTime end_time(start_time + duration);
346 logManager.
addProperty(
"start_time", start_time.toISO8601String());
347 logManager.
addProperty(
"run_start", start_time.toISO8601String());
348 logManager.
addProperty(
"end_time", end_time.toISO8601String());
351 std::string time_str = start_time.toISO8601String();
359 for (
auto &
x : logStrings) {
362 for (
auto &
x : logParams) {
367 loadInstrumentAlg->setProperty(
"Workspace", eventWS);
368 loadInstrumentAlg->setPropertyValue(
"InstrumentName",
"BILBY");
370 loadInstrumentAlg->executeAsChildAlg();
379 if (maskfile.length() == 0)
382 std::ifstream input(maskfile.c_str());
384 throw std::invalid_argument(
"invalid mask file");
387 while (std::getline(input, line)) {
388 auto i0 = line.find(
"<detids>");
389 auto iN = line.find(
"</detids>");
391 if ((i0 != std::string::npos) && (iN != std::string::npos) && (i0 < iN)) {
392 line = line.substr(i0 + 8, iN - i0 - 8);
393 std::stringstream ss(line);
396 while (std::getline(ss, item,
',')) {
397 auto k = item.find(
'-');
400 if (k != std::string::npos) {
401 p0 = boost::lexical_cast<size_t>(item.substr(0, k));
402 p1 = boost::lexical_cast<size_t>(item.substr(k + 1, item.size() - k - 1));
407 p0 = boost::lexical_cast<size_t>(item);
411 if (p0 < result.size()) {
412 if (p1 >= result.size())
413 p1 = result.size() - 1;
416 result[p0++] =
false;
427 std::map<std::string, std::string> &logStrings,
428 std::map<std::string, std::string> &allParams) {
430 std::string idfDirectory = Mantid::Kernel::ConfigService::Instance().getString(
"instrumentDefinition.directory");
433 std::string parameterFilename = idfDirectory +
"BILBY_Parameters.xml";
436 Poco::AutoPtr<Poco::XML::Document> pDoc;
438 pDoc = pParser.parse(parameterFilename);
442 NodeIterator it(pDoc, Poco::XML::NodeFilter::SHOW_ELEMENT);
443 Node *pNode = it.nextNode();
445 if (pNode->nodeName() ==
"parameter") {
446 auto pElem =
dynamic_cast<Element *
>(pNode);
447 std::string paramName = pElem->getAttribute(
"name");
448 Poco::AutoPtr<NodeList> nodeList = pElem->childNodes();
449 for (
unsigned long i = 0; i < nodeList->length(); i++) {
450 auto cNode = nodeList->item(i);
451 if (cNode->nodeName() ==
"value") {
452 auto cElem =
dynamic_cast<Poco::XML::Element *
>(cNode);
453 std::string
value = cElem->getAttribute(
"val");
454 allParams[paramName] = std::move(
value);
458 pNode = it.nextNode();
461 auto isNumeric = [](
const std::string &tag) {
463 auto stag = boost::algorithm::trim_copy(tag);
465 auto value = std::stod(stag, &sz);
466 return sz > 0 && stag.size() == sz && std::isfinite(
value);
467 }
catch (
const std::invalid_argument &) {
472 std::string tmpString;
473 float tmpFloat = 0.0f;
474 for (
auto &
x : allParams) {
475 if (
x.first.find(
"log_") == 0) {
476 auto logTag = boost::algorithm::trim_copy(
x.first.substr(4));
477 auto line =
x.second;
480 std::vector<std::string> details;
481 boost::split(details, line, boost::is_any_of(
","));
482 if (details.size() < 3) {
483 g_log.
warning() <<
"Invalid format for BILBY parameter " <<
x.first << std::endl;
486 auto hdfTag = boost::algorithm::trim_copy(details[0]);
490 auto updateOk =
false;
491 if (!hdfTag.empty()) {
492 if (isNumeric(details[1])) {
494 auto factor = std::stod(details[1]);
495 logParams[logTag] = factor * tmpFloat;
499 logStrings[logTag] = tmpString;
506 auto defValue = boost::algorithm::trim_copy(details[2]);
507 if (!defValue.empty()) {
508 if (isNumeric(defValue))
509 logParams[logTag] = std::stod(defValue);
511 logStrings[logTag] = std::move(defValue);
513 g_log.
warning() <<
"Cannot find hdf parameter " << hdfTag <<
", using default.\n";
516 }
catch (
const std::invalid_argument &) {
517 g_log.
warning() <<
"Invalid format for BILBY parameter " <<
x.first << std::endl;
521 }
catch (std::exception &ex) {
522 g_log.
warning() <<
"Failed to load instrument with error: " << ex.what()
523 <<
". The current facility may not be fully "
530 std::map<std::string, double> &logParams, std::map<std::string, std::string> &logStrings,
531 std::map<std::string, std::string> &allParams) {
533 const double toMeters = 1.0 / 1000;
537 instrumentInfo.
start_time =
"2000-01-01T00:00:00";
544 instrumentInfo.
is_tof =
true;
552 const std::vector<std::string> &files = tarFile.
files();
553 auto file_it = std::find_if(files.cbegin(), files.cend(),
554 [](
const std::string &file) { return file.rfind(
".hdf") == file.length() - 4; });
555 if (file_it != files.end()) {
556 tarFile.
select(file_it->c_str());
558 Poco::TemporaryFile hdfFile;
559 std::shared_ptr<FILE> handle(fopen(hdfFile.path().c_str(),
"wb"), fclose);
564 while (0 != (bytesRead = tarFile.
read(buffer,
sizeof(buffer))))
565 fwrite(buffer, bytesRead, 1, handle.get());
571 float tmp_float = 0.0f;
572 int32_t tmp_int32 = 0;
578 instrumentInfo.
att_pos = boost::math::iround(tmp_float);
587 if (
loadNXDataSet(entry,
"instrument/master1_chopper_id", tmp_int32))
589 if (
loadNXDataSet(entry,
"instrument/master2_chopper_id", tmp_int32))
592 if (
loadNXString(entry,
"instrument/detector/frame_source", tmp_str))
593 instrumentInfo.
is_tof = tmp_str ==
"EXTERNAL";
594 if (
loadNXDataSet(entry,
"instrument/nvs067/lambda", tmp_float))
597 if (
loadNXDataSet(entry,
"instrument/master_chopper_freq", tmp_float) && (tmp_float > 0.0f))
599 if (
loadNXDataSet(entry,
"instrument/t0_chopper_freq", tmp_float) && (tmp_float > 0.0f))
601 if (
loadNXDataSet(entry,
"instrument/t0_chopper_phase", tmp_float))
602 instrumentInfo.
phase_slave = tmp_float < 999.0 ? tmp_float : 0.0;
608 auto findLtof = logParams.find(
"Ltof_det_value");
609 if (findLtof != logParams.end()) {
610 logParams[
"L1_chopper_value"] = logParams[
"Ltof_det_value"] - logParams[
"L2_det_value"];
612 logParams[
"L1_chopper_value"] = 18.4726;
613 g_log.
warning() <<
"Cannot recover parameter 'L1_chopper_value'"
614 <<
", using default.\n";
620 file_it = std::find(files.cbegin(), files.cend(),
"History.log");
621 if (file_it != files.cend()) {
622 tarFile.
select(file_it->c_str());
623 std::string logContent;
625 tarFile.
read(&logContent[0], logContent.size());
626 std::istringstream data(logContent);
627 Poco::AutoPtr<Poco::Util::PropertyFileConfiguration> conf(
new Poco::Util::PropertyFileConfiguration(data));
629 if (conf->hasProperty(
"Bm1Counts"))
630 instrumentInfo.
bm_counts = conf->getInt(
"Bm1Counts");
631 if (conf->hasProperty(
"AttPos"))
632 instrumentInfo.
att_pos = boost::math::iround(conf->getDouble(
"AttPos"));
634 if (conf->hasProperty(
"SampleName"))
635 instrumentInfo.
sample_name = conf->getString(
"SampleName");
637 if (conf->hasProperty(
"MasterChopperFreq")) {
638 auto tmp = conf->getDouble(
"MasterChopperFreq");
642 if (conf->hasProperty(
"T0ChopperFreq")) {
643 auto tmp = conf->getDouble(
"T0ChopperFreq");
647 if (conf->hasProperty(
"T0ChopperPhase")) {
648 auto tmp = conf->getDouble(
"T0ChopperPhase");
652 if (conf->hasProperty(
"FrameSource"))
653 instrumentInfo.
is_tof = conf->getString(
"FrameSource") ==
"EXTERNAL";
654 if (conf->hasProperty(
"Wavelength"))
655 instrumentInfo.
wavelength = conf->getDouble(
"Wavelength");
657 if (conf->hasProperty(
"SampleAperture"))
658 logParams[
"sample_aperture"] = conf->getDouble(
"SampleAperture");
659 if (conf->hasProperty(
"SourceAperture"))
660 logParams[
"source_aperture"] = conf->getDouble(
"SourceAperture");
661 if (conf->hasProperty(
"L1"))
662 logParams[
"L1_source_value"] = conf->getDouble(
"L1") * toMeters;
663 if (conf->hasProperty(
"LTofDet"))
664 logParams[
"L1_chopper_value"] = conf->getDouble(
"LTofDet") * toMeters - logParams[
"L2_det_value"];
665 if (conf->hasProperty(
"L2Det"))
666 logParams[
"L2_det_value"] = conf->getDouble(
"L2Det") * toMeters;
668 if (conf->hasProperty(
"L2CurtainL"))
669 logParams[
"L2_curtainl_value"] = conf->getDouble(
"L2CurtainL") * toMeters;
670 if (conf->hasProperty(
"L2CurtainR"))
671 logParams[
"L2_curtainr_value"] = conf->getDouble(
"L2CurtainR") * toMeters;
672 if (conf->hasProperty(
"L2CurtainU"))
673 logParams[
"L2_curtainu_value"] = conf->getDouble(
"L2CurtainU") * toMeters;
674 if (conf->hasProperty(
"L2CurtainD"))
675 logParams[
"L2_curtaind_value"] = conf->getDouble(
"L2CurtainD") * toMeters;
677 if (conf->hasProperty(
"CurtainL"))
678 logParams[
"D_curtainl_value"] = conf->getDouble(
"CurtainL") * toMeters;
679 if (conf->hasProperty(
"CurtainR"))
680 logParams[
"D_curtainr_value"] = conf->getDouble(
"CurtainR") * toMeters;
681 if (conf->hasProperty(
"CurtainU"))
682 logParams[
"D_curtainu_value"] = conf->getDouble(
"CurtainU") * toMeters;
683 if (conf->hasProperty(
"CurtainD"))
684 logParams[
"D_curtaind_value"] = conf->getDouble(
"CurtainD") * toMeters;
696 }
catch (std::runtime_error &) {
704 }
catch (
const std::runtime_error &) {
710template <
class EventProcessor>
716 int64_t fileSize = 0;
717 const std::vector<std::string> &files = tarFile.
files();
718 const auto found = std::find_if(files.cbegin(), files.cend(),
719 [](
const auto &file) { return file.rfind(
".bin") == file.length() - 4; });
720 if (found != files.cend()) {
721 tarFile.
select(found->c_str());
736 if ((fileSize == 0) || !tarFile.
skip(128))
740 int invalidEvents = 0;
742 while ((c =
static_cast<uint32_t
>(tarFile.
read_byte())) !=
static_cast<uint32_t
>(-1)) {
744 bool event_ended =
false;
751 x |= (c & 0x01) << 8;
756 event_ended = (c & 0xC0) != 0xC0;
760 y |= (c & 0x01) << 7;
761 dt = (c & 0xFE) >> 1;
769 event_ended = (c & 0xC0) != 0xC0;
774 dt |= (c & 0xFF) << (5 + 6 * (state - 3));
779 if (event_ended || (state == 8)) {
782 if ((
x == 0) && (
y == 0) && (dt == 0xFFFFFFFF)) {
788 tof +=
static_cast<int>(dt) * 0.1;
792 tof +=
static_cast<int>(dt) * 0.1;
800 if (invalidEvents > 0) {
801 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 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.
void addEvent(size_t x, size_t y, double tof)
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
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)
bool loadNXString(const Nexus::NXEntry &entry, const std::string &address, std::string &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)
static bool loadNXDataSet(const Nexus::NXEntry &entry, const std::string &address, T &value)
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 exec() override
Execute the algorithm.
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.
A specialised Property class for holding a series of time-value pairs.
std::string getString(const std::string &name) const
Returns a string.
NXDataSetTyped< T > openNXDataSet(const std::string &name) const
Templated method for creating datasets.
Templated class implementation of NXDataSet.
void load()
Read all of the datablock in.
Implements NXentry Nexus class.
Implements NXroot Nexus class.
NXEntry openFirstEntry()
Open the first NXentry in the file.
std::vector< Types::Event::TofEvent > * EventVector_pt
pointer to the vector of events
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.