22#include <boost/algorithm/string.hpp>
24#include <nexus/NeXusFile.hpp>
25#include <nexus/NeXusException.hpp>
29using namespace Kernel;
31using namespace DataObjects;
52 const std::vector<std::string> exts{
".h5",
".nxs"};
54 "The name of the Nexus file to load");
57 "An output workspace.");
60 "When this property is set to false errors are set equal to data values, "
61 "and when set to true all errors are set equal to one. This property "
65 "When true the algorithm only outputs the sum of all event data into "
66 "one eventworkspace EventData + _ + name of the OutputWorkspace. "
67 "If false eventworkspaces are also returned for each individual "
68 "McStas components storing event data");
77 g_log.
debug() <<
"Opening file " << filename <<
'\n';
79 ::NeXus::File nxFile(filename);
80 auto entries = nxFile.getEntries();
83 auto entry = entries.begin();
84 std::string
name = entry->first;
85 std::string type = entry->second;
88 nxFile.openGroup(
name, type);
89 nxFile.openGroup(
"data",
"NXdetector");
91 auto dataEntries = nxFile.getEntries();
93 std::map<std::string, std::string> eventEntries;
94 std::map<std::string, std::string> histogramEntries;
97 for (
auto &dataEntry : dataEntries) {
98 std::string dataName = dataEntry.first;
99 std::string dataType = dataEntry.second;
100 if (dataName ==
"content_nxs" || dataType !=
"NXdata")
103 g_log.
debug() <<
"Opening " << dataName <<
" " << dataType <<
'\n';
106 nxFile.openGroup(dataName, dataType);
114 auto nxdataEntries = nxFile.getEntries();
116 for (
auto &nxdataEntry : nxdataEntries) {
117 if (nxdataEntry.second ==
"NXparameters")
119 nxFile.openData(nxdataEntry.first);
120 if (nxFile.hasAttr(
"long_name")) {
121 std::string nameAttrValue;
122 nxFile.getAttr(
"long_name", nameAttrValue);
124 if (nameAttrValue.find(
"Neutron_ID") != std::string::npos) {
125 eventEntries[dataEntry.first] = dataEntry.second;
127 histogramEntries[dataEntry.first] = dataEntry.second;
135 std::vector<std::string> scatteringWSNames;
136 std::vector<std::string> histoWSNames;
137 if (!eventEntries.empty()) {
143 scatteringWSNames.insert(scatteringWSNames.end(), histoWSNames.begin(), histoWSNames.end());
159 groupAlgorithm->setChild(
true);
160 groupAlgorithm->setLogging(
false);
161 groupAlgorithm->initialize();
162 groupAlgorithm->setProperty(
"InputWorkspaces", workspaces);
163 groupAlgorithm->setProperty(
"OutputWorkspace",
"__grouped");
164 groupAlgorithm->execute();
165 return groupAlgorithm->getProperty(
"OutputWorkspace");
175 ::NeXus::File &nxFile) {
178 std::vector<std::string> scatteringWSNames;
181 auto entries = nxFile.getEntries();
182 const bool errorBarsSetTo1 =
getProperty(
"ErrorBarsSetTo1");
198 const double progressFractionInitial = 0.1;
199 Progress progInitial(
this, 0.0, progressFractionInitial, reports);
201 std::string instrumentXML;
202 progInitial.
report(
"Loading instrument");
204 nxFile.openGroup(
"instrument",
"NXinstrument");
205 nxFile.openGroup(
"instrument_xml",
"NXnote");
206 nxFile.readData(
"data", instrumentXML);
210 g_log.
warning() <<
"\nCould not find the instrument description in the Nexus file:" << filename
211 <<
" Ignore eventdata from the Nexus file\n";
212 return scatteringWSNames;
217 std::string instrumentName =
"McStas";
227 instrument = parser.
parseXML(
nullptr);
232 g_log.
warning() <<
"When trying to read the instrument description in the Nexus file: " << filename
233 <<
" the following error is reported: " << e.
what() <<
" Ignore eventdata from the Nexus file\n";
234 return scatteringWSNames;
237 g_log.
warning() <<
"Could not parse instrument description in the Nexus file: " << filename
238 <<
" Ignore eventdata from the Nexus file\n";
239 return scatteringWSNames;
243 nxFile.openGroup(
"data",
"NXdetector");
246 progInitial.
report(
"Set up EventWorkspace");
250 eventWS->initialize(instrument->getNumberDetectors(), 1, 1);
253 eventWS->setYUnit(
"Counts");
255 eventWS->setInstrument(instrument);
258 std::vector<detid_t> detIDs = instrument->getDetectorIDs();
260 for (
size_t i = 0; i < instrument->getNumberDetectors(); i++) {
261 eventWS->getSpectrum(i).addDetectorID(detIDs[i]);
263 eventWS->getSpectrum(i).setSpectrumNo(detIDs[i]);
266 eventWS->rebuildSpectraMapping(
true);
268 bool isAnyNeutrons =
false;
270 double shortestTOF(0.0);
271 double longestTOF(0.0);
274 const size_t numEventEntries = eventEntries.size();
275 std::string nameOfGroupWS =
getProperty(
"OutputWorkspace");
276 const auto eventDataTotalName =
"EventData_" + nameOfGroupWS;
277 std::vector<std::pair<EventWorkspace_sptr, std::string>> allEventWS = {{eventWS, eventDataTotalName}};
279 const bool onlySummedEventWorkspace =
getProperty(
"OutputOnlySummedEventWorkspace");
280 if (!onlySummedEventWorkspace && numEventEntries > 1) {
281 for (
const auto &eventEntry : eventEntries) {
282 const std::string &dataName = eventEntry.first;
285 const auto ws_name = dataName +
"_" + nameOfGroupWS;
286 allEventWS.emplace_back(eventWS->clone(), ws_name);
290 Progress progEntries(
this, progressFractionInitial, 1.0, numEventEntries * 2);
293 auto eventWSIndex = 1u;
295 for (
const auto &eventEntry : eventEntries) {
296 const std::string &dataName = eventEntry.first;
297 const std::string &dataType = eventEntry.second;
300 nxFile.openGroup(dataName, dataType);
301 std::vector<double> data;
302 nxFile.openData(
"events");
303 progEntries.
report(
"read event data from nexus");
319 ::NeXus::Info id_info = nxFile.getInfo();
320 if (id_info.dims.size() != 2) {
321 g_log.
error() <<
"Event data in McStas nexus file not loaded. Expected "
322 "event data block to be two dimensional\n";
323 return scatteringWSNames;
326 int64_t nNeutrons = id_info.dims[0];
327 int64_t numberOfDataColumn = id_info.dims[1];
328 if (nNeutrons && numberOfDataColumn != 6) {
329 g_log.
error() <<
"Event data in McStas nexus file expecting 6 columns\n";
330 return scatteringWSNames;
333 if (!isAnyNeutrons && nNeutrons > 0)
334 isAnyNeutrons =
true;
336 std::vector<int64_t> start(2);
337 std::vector<int64_t> step(2);
341 int64_t nNeutronsInBlock = 1000000;
342 int64_t nOfFullBlocks = nNeutrons / nNeutronsInBlock;
343 int64_t nRemainingNeutrons = nNeutrons - nOfFullBlocks * nNeutronsInBlock;
345 for (int64_t iBlock = 0; iBlock < nOfFullBlocks + 1; iBlock++) {
346 if (iBlock == nOfFullBlocks) {
348 start[0] = nOfFullBlocks * nNeutronsInBlock;
350 step[0] = nRemainingNeutrons;
351 step[1] = numberOfDataColumn;
354 start[0] = iBlock * nNeutronsInBlock;
356 step[0] = nNeutronsInBlock;
357 step[1] = numberOfDataColumn;
359 const int64_t nNeutronsForthisBlock = step[0];
360 data.resize(nNeutronsForthisBlock * numberOfDataColumn);
363 if (id_info.type == ::NeXus::FLOAT64) {
364 nxFile.getSlab(&data[0], start, step);
366 g_log.
warning() <<
"Entry event field is not FLOAT64! It will be skipped.\n";
371 const detid2index_map detIDtoWSindex_map = allEventWS[0].first->getDetectorIDToWorkspaceIndexMap(
true);
373 progEntries.
report(
"read event data into workspace");
374 for (int64_t in = 0; in < nNeutronsForthisBlock; in++) {
375 const auto detectorID =
static_cast<int>(data[4 + numberOfDataColumn * in]);
376 const double detector_time = data[5 + numberOfDataColumn * in] * 1.0e6;
377 if (in == 0 && iBlock == 0) {
378 shortestTOF = detector_time;
379 longestTOF = detector_time;
381 if (detector_time < shortestTOF)
382 shortestTOF = detector_time;
383 if (detector_time > longestTOF)
384 longestTOF = detector_time;
387 const size_t workspaceIndex = detIDtoWSindex_map.find(detectorID)->second;
389 int64_t pulse_time = 0;
391 if (errorBarsSetTo1) {
392 weightedEvent =
WeightedEvent(detector_time, pulse_time, data[numberOfDataColumn * in], 1.0);
394 weightedEvent =
WeightedEvent(detector_time, pulse_time, data[numberOfDataColumn * in],
395 data[numberOfDataColumn * in] * data[numberOfDataColumn * in]);
397 allEventWS[0].first->getSpectrum(workspaceIndex) += weightedEvent;
398 if (!onlySummedEventWorkspace && numEventEntries > 1) {
399 allEventWS[eventWSIndex].first->getSpectrum(workspaceIndex) += weightedEvent;
415 auto axis = HistogramData::BinEdges{shortestTOF - 1, longestTOF + 1};
419 for (
const auto &wsAndName : allEventWS) {
420 const auto ws = wsAndName.first;
423 scatteringWSNames.emplace_back(wsAndName.second);
425 return scatteringWSNames;
435 ::NeXus::File &nxFile) {
437 std::string nameAttrValueYLABEL;
438 std::vector<std::string> histoWSNames;
440 for (
const auto &histogramEntry : histogramEntries) {
441 const std::string &dataName = histogramEntry.first;
442 const std::string &dataType = histogramEntry.second;
445 nxFile.openGroup(dataName, dataType);
448 std::string nameAttrValueTITLE;
449 nxFile.getAttr(
"filename", nameAttrValueTITLE);
451 if (nxFile.hasAttr(
"ylabel")) {
452 nxFile.getAttr(
"ylabel", nameAttrValueYLABEL);
456 auto nxdataEntries = nxFile.getEntries();
457 std::string axis1Name, axis2Name;
458 for (
auto &nxdataEntry : nxdataEntries) {
459 if (nxdataEntry.second ==
"NXparameters")
461 if (nxdataEntry.first ==
"ncount")
463 nxFile.openData(nxdataEntry.first);
465 if (nxFile.hasAttr(
"axis")) {
467 nxFile.getAttr(
"axis", axisNo);
469 axis1Name = nxdataEntry.first;
470 else if (axisNo == 2)
471 axis2Name = nxdataEntry.first;
473 throw std::invalid_argument(
"Unknown axis number");
478 std::vector<double> axis1Values;
479 std::vector<double> axis2Values;
480 nxFile.readData<
double>(axis1Name, axis1Values);
481 if (axis2Name.length() == 0) {
482 axis2Name = nameAttrValueYLABEL;
483 axis2Values.emplace_back(0.0);
485 nxFile.readData<
double>(axis2Name, axis2Values);
488 const size_t axis1Length = axis1Values.size();
489 const size_t axis2Length = axis2Values.size();
490 g_log.
debug() <<
"Axis lengths=" << axis1Length <<
" " << axis2Length <<
'\n';
493 std::vector<double> data;
494 nxFile.readData<
double>(
"data", data);
497 std::vector<double> errors;
499 nxFile.readData<
double>(
"errors", errors);
500 }
catch (::NeXus::Exception &) {
501 g_log.
information() <<
"Field " << dataName <<
" contains no error information.\n";
508 Axis *axis1 = ws->getAxis(0);
509 axis1->
title() = axis1Name;
511 auto lblUnit = std::make_shared<Units::Label>();
512 lblUnit->setLabel(axis1Name,
"");
513 axis1->
unit() = lblUnit;
515 auto axis2 = std::make_unique<NumericAxis>(axis2Length);
516 auto axis2Raw = axis2.get();
517 axis2->title() = axis2Name;
519 lblUnit = std::make_shared<Units::Label>();
520 lblUnit->setLabel(axis2Name,
"");
521 axis2->unit() = lblUnit;
523 ws->setYUnit(axis2Name);
524 ws->replaceAxis(1, std::move(axis2));
526 for (
size_t wsIndex = 0; wsIndex < axis2Length; ++wsIndex) {
527 auto &dataX = ws->mutableX(wsIndex);
528 auto &dataY = ws->mutableY(wsIndex);
529 auto &dataE = ws->mutableE(wsIndex);
531 for (
size_t j = 0; j < axis1Length; ++j) {
534 const size_t fileDataIndex = j * axis2Length + wsIndex;
536 dataX[j] = axis1Values[j];
537 dataY[j] = data[fileDataIndex];
539 dataE[j] = errors[fileDataIndex];
541 axis2Raw->setValue(wsIndex, axis2Values[wsIndex]);
545 ws->setTitle(nameAttrValueTITLE);
548 std::replace(nameAttrValueTITLE.begin(), nameAttrValueTITLE.end(),
' ',
'_');
552 const std::string outputWS =
getProperty(
"OutputWorkspace");
553 const std::string nameUserSee = nameAttrValueTITLE +
"_" + outputWS;
556 histoWSNames.emplace_back(ws->getName());
570 using namespace ::
NeXus;
574 if (descriptor.
pathExists(
"/entry1/simulation/name")) {
577 ::NeXus::File file = ::NeXus::File(descriptor.
filename());
579 file.openGroup(
"simulation",
"NXnote");
582 file.readData(
"name",
value);
583 if (boost::iequals(
value,
"mccode"))
587 }
catch (::NeXus::Exception &) {
double value
The value of the point.
#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.
Class to represent the axis of a workspace.
const std::string & title() const
Returns the user-defined title for this axis.
const std::shared_ptr< Kernel::Unit > & unit() const
The unit for this axis.
@ Load
allowed here which will be passed to the algorithm
Helper class for reporting progress from algorithms.
A property class for workspaces.
LoadMcStas : TODO: DESCRIPTION.
void init() override
Initialize the algorithm's properties.
int confidence(Kernel::NexusDescriptor &descriptor) const override
Returns a confidence value that this algorithm can load a file.
API::WorkspaceGroup_sptr groupWorkspaces(const std::vector< std::string > &workspaces) const
Group workspaces.
void exec() override
Execute the algorithm.
const std::string name() const override
function to return a name of the algorithm, must be overridden in all algorithms
const std::string category() const override
function to return a category of the algorithm.
int version() const override
function to return a version of the algorithm, must be overridden in all algorithms
std::vector< std::string > readEventData(const std::map< std::string, std::string > &eventEntries, ::NeXus::File &nxFile)
Read Event Data.
std::vector< std::string > readHistogramData(const std::map< std::string, std::string > &histogramEntries, ::NeXus::File &nxFile)
Read histogram data.
This class is intended to fulfill the design specified in <https://github.com/mantidproject/documents...
Info about a single neutron detection event, including a weight and error value:
Creates an instrument data from a XML instrument description file.
std::shared_ptr< Instrument > parseXML(Kernel::ProgressBase *progressReporter)
Parse XML contents.
std::string getMangledName()
Handle used in the singleton constructor for instrument file should append the value file sha-1 check...
Exception for errors associated with the instrument definition.
const char * what() const noexcept override
Writes out the range and limits.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
void error(const std::string &msg)
Logs at error level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
Defines a wrapper around a file whose internal structure can be accessed using the NeXus API.
const std::pair< std::string, std::string > & firstEntryNameType() const
Returns the name & type of the first entry in the file.
bool pathExists(const std::string &path) const
Query if a path exists.
const std::string & filename() const
Access the filename.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< WorkspaceGroup > WorkspaceGroup_sptr
shared pointer to Mantid::API::WorkspaceGroup
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< EventWorkspace > EventWorkspace_sptr
shared pointer to the EventWorkspace class
std::shared_ptr< Instrument > Instrument_sptr
Shared pointer to an instrument object.
std::unordered_map< detid_t, size_t > detid2index_map
Map with key = detector ID, value = workspace index.
@ Output
An output workspace.