LoadEventNexus.h

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2010 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#pragma once
 
#include "MantidAPI/IFileLoader.h"
#include "MantidAPI/InstrumentFileFinder.h"
#include "MantidAPI/NexusFileLoader.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidDataHandling/BankPulseTimes.h"
#include "MantidDataHandling/DllConfig.h"
#include "MantidDataHandling/EventWorkspaceCollection.h"
#include "MantidDataHandling/LoadGeometry.h"
#include "MantidDataObjects/EventList.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/Events.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/NexusHDF5Descriptor.h"
#include "MantidKernel/OptionalBool.h"
#include "MantidKernel/TimeSeriesProperty.h"
 
#include <Poco/Path.h>
#include <boost/lexical_cast.hpp>
#include <boost/scoped_array.hpp>
 
// clang-format off
#include <nexus/NeXusFile.hpp>
#include <nexus/NeXusException.hpp>
// clang-format on
 
#include <algorithm>
#include <cstdint>
#include <functional>
#include <memory>
#include <mutex>
#include <numeric>
#include <random>
#include <string>
 
namespace Mantid {
namespace DataHandling {
 
class InvalidLogPeriods : public std::invalid_argument {
public:
  InvalidLogPeriods(const std::string &msg) : std::invalid_argument(msg) {}
};
 
bool exists(::NeXus::File &file, const std::string &name);
 
bool exists(const std::map<std::string, std::string> &entries, const std::string &name);
 
class MANTID_DATAHANDLING_DLL LoadEventNexus : public API::NexusFileLoader {
 
public:
  LoadEventNexus();
 
  const std::string name() const override { return "LoadEventNexus"; };
 
  const std::string summary() const override {
    return "Loads an Event NeXus file and stores as an "
           "EventWorkspace. Optionally, you can filter out events falling "
           "outside a range of times-of-flight and/or a time interval.";
  }
 
  int version() const override { return 1; };
  const std::vector<std::string> seeAlso() const override { return {"LoadISISNexus", "LoadEventAndCompress"}; }
 
  const std::string category() const override { return "DataHandling\\Nexus"; }
 
  int confidence(Kernel::NexusHDF5Descriptor &descriptor) const override;
 
  template <typename T>
  static std::shared_ptr<BankPulseTimes>
  runLoadNexusLogs(const std::string &nexusfilename, T localWorkspace, Algorithm &alg, bool returnpulsetimes,
                   int &nPeriods, std::unique_ptr<const Kernel::TimeSeriesProperty<int>> &periodLog);
 
  template <typename T>
  static std::shared_ptr<BankPulseTimes>
  runLoadNexusLogs(const std::string &nexusfilename, T localWorkspace, Algorithm &alg, bool returnpulsetimes,
                   int &nPeriods, std::unique_ptr<const Kernel::TimeSeriesProperty<int>> &periodLog,
                   const std::vector<std::string> &allow_list, const std::vector<std::string> &block_list);
 
  static void checkForCorruptedPeriods(std::unique_ptr<Kernel::TimeSeriesProperty<int>> tempPeriodLog,
                                       std::unique_ptr<const Kernel::TimeSeriesProperty<int>> &periodLog,
                                       const int &nPeriods, const std::string &nexusfilename);
 
  template <typename T>
  static void loadEntryMetadata(const std::string &nexusfilename, T WS, const std::string &entry_name,
                                const Kernel::NexusHDF5Descriptor &descriptor);
 
  template <typename T>
  static bool loadInstrument(const std::string &nexusfilename, T localWorkspace, const std::string &top_entry_name,
                             Algorithm *alg, const Kernel::NexusHDF5Descriptor *descriptor = nullptr);
 
  template <typename T>
  static bool runLoadIDFFromNexus(const std::string &nexusfilename, T localWorkspace, const std::string &top_entry_name,
                                  Algorithm *alg);
 
  template <typename T>
  static bool runLoadInstrument(const std::string &nexusfilename, T localWorkspace, const std::string &top_entry_name,
                                Algorithm *alg, const Kernel::NexusHDF5Descriptor *descriptor = nullptr);
 
  static void loadSampleDataISIScompatibility(::NeXus::File &file, EventWorkspaceCollection &WS);
 
  static std::string readInstrumentFromISIS_VMSCompat(::NeXus::File &hFile);
 
public:
  std::string m_filename;
 
  std::shared_ptr<EventWorkspaceCollection> m_ws;
 
  double filter_tof_min;
  double filter_tof_max;
 
  int32_t m_specMin;
  int32_t m_specMax;
 
  Mantid::Types::Core::DateAndTime filter_time_start;
  Mantid::Types::Core::DateAndTime filter_time_stop;
 
  std::mutex m_tofMutex;
 
  double longest_tof;
  double shortest_tof;
  size_t bad_tofs;
  size_t discarded_events;
 
  double compressTolerance;
 
  std::shared_ptr<BankPulseTimes> m_allBanksPulseTimes;
 
  std::string m_top_entry_name;
  std::unique_ptr<::NeXus::File> m_file;
 
protected:
  Parallel::ExecutionMode
  getParallelExecutionMode(const std::map<std::string, Parallel::StorageMode> &storageModes) const override;
 
private:
  enum class LoaderType;
 
  void init() override;
 
  void execLoader() override;
 
  LoadEventNexus::LoaderType defineLoaderType(const bool haveWeights, const bool oldNeXusFileNames,
                                              const std::string &classType) const;
 
  DataObjects::EventWorkspace_sptr createEmptyEventWorkspace();
 
  void loadEvents(API::Progress *const prog, const bool monitors);
  void createSpectraMapping(const std::string &nxsfile, const bool monitorsOnly,
                            const std::vector<std::string> &bankNames = std::vector<std::string>());
  void deleteBanks(const EventWorkspaceCollection_sptr &workspace, const std::vector<std::string> &bankNames);
  void runLoadMonitors();
  void setTimeFilters(const bool monitors);
  template <typename T>
  T filterEventsByTime(T workspace, Mantid::Types::Core::DateAndTime &startTime,
                       Mantid::Types::Core::DateAndTime &stopTime);
 
  std::unique_ptr<std::pair<std::vector<int32_t>, std::vector<int32_t>>>
  loadISISVMSSpectraMapping(const std::string &entry_name);
 
  template <typename T> void filterDuringPause(T workspace);
 
  void setTopEntryName();
 
  void safeOpenFile(const std::string &fname);
 
  bool m_instrument_loaded_correctly;
 
  bool loadlogs;
  bool event_id_is_spec;
};
 
//-----------------------------------------------------------------------------
//               ISIS event corrections
//-----------------------------------------------------------------------------
 
template <typename T>
void makeTimeOfFlightDataFuzzy(::NeXus::File &file, T localWorkspace, const std::string &binsName, size_t start_wi = 0,
                               size_t end_wi = 0) {
  const std::string EVENT_TIME_SHIFT_TAG("event_time_offset_shift");
  // first check if the data is already randomized
  const auto entries = file.getEntries();
  if (entries.find(EVENT_TIME_SHIFT_TAG) != entries.end()) {
    std::string event_shift_type;
    file.readData(EVENT_TIME_SHIFT_TAG, event_shift_type);
    if (event_shift_type == "random") {
      return;
    }
  }
 
  // if the data is not randomized randomize it uniformly within each bin
  file.openData(binsName);
  // time of flights of events
  std::vector<float> tofsFile;
  file.getData(tofsFile);
  file.closeData();
 
  // todo: try to find if tof can be reduced to just 3 numbers: start, end and
  // dt
  if (end_wi <= start_wi) {
    end_wi = localWorkspace->getNumberHistograms();
  }
 
  // random number generator
  std::mt19937 rng;
 
  // loop over spectra
  for (size_t wi = start_wi; wi < end_wi; ++wi) {
    DataObjects::EventList &event_list = dynamic_cast<DataObjects::EventList &>(localWorkspace->getSpectrum(wi));
    if (event_list.empty())
      continue;
    // sort the events
    event_list.sortTof();
    auto tofsEventList = event_list.getTofs();
 
    size_t n = tofsFile.size();
    // iterate over the events and time bins
    auto ev = tofsEventList.begin();
    auto ev_end = tofsEventList.end();
    for (size_t i = 1; i < n; ++i) {
      double right = double(tofsFile[i]);
      // find the right boundary for the current event
      if ((ev != ev_end) && (right < *ev)) {
        continue;
      }
      // count events which have the same right boundary
      size_t m = 0;
      while ((ev != ev_end) && (*ev < right)) {
        ++ev;
        ++m; // count events in the i-th bin
      }
 
      if (m > 0) { // m events in this bin
        double left = double(tofsFile[i - 1]);
        // spread the events uniformly inside the bin
        std::uniform_real_distribution<double> flat(left, right);
        std::vector<double> random_numbers(m);
        std::generate(random_numbers.begin(), random_numbers.end(), [&flat, &rng]() { return flat(rng); });
        std::sort(random_numbers.begin(), random_numbers.end());
        auto it = random_numbers.begin();
        for (auto ev1 = ev - m; ev1 != ev; ++ev1, ++it) {
          *ev1 = *it;
        }
      }
 
    } // for i
    event_list.setTofs(tofsEventList);
 
    event_list.sortTof();
  } // for wi
}
 
template <typename T>
void adjustTimeOfFlightISISLegacy(::NeXus::File &file, T localWorkspace, const std::string &entry_name,
                                  const std::string &classType,
                                  const Kernel::NexusHDF5Descriptor *descriptor = nullptr) {
  bool done = false;
  // Go to the root, and then top entry
  file.openPath("/");
  file.openGroup(entry_name, "NXentry");
 
  // NexusHDF5Descriptor
  if (descriptor != nullptr) {
    // not an ISIS file
    if (!descriptor->isEntry("/" + entry_name + "/detector_1_events")) {
      return;
    }
  }
 
  using string_map_t = std::map<std::string, std::string>;
  string_map_t entries = file.getEntries();
 
  if (entries.find("detector_1_events") == entries.end()) { // not an ISIS file
    return;
  }
 
  // try if monitors have their own bins
  if (classType == "NXmonitor") {
    std::vector<std::string> bankNames;
    for (string_map_t::const_iterator it = entries.begin(); it != entries.end(); ++it) {
      std::string entryName(it->first);
      std::string entry_class(it->second);
      if (entry_class == classType) {
        bankNames.emplace_back(entryName);
      }
    }
    for (size_t i = 0; i < bankNames.size(); ++i) {
      const std::string &mon = bankNames[i];
      file.openGroup(mon, classType);
      entries = file.getEntries();
      if (entries.find("event_time_bins") == entries.end()) {
        // bins = entries.find("time_of_flight"); // I think time_of_flight
        // doesn't work here
        // if (bins == entries.end())
        //{
        done = false;
        file.closeGroup();
        break; // done == false => use bins from the detectors
               //}
      }
      done = true;
      makeTimeOfFlightDataFuzzy(file, localWorkspace, "event_time_bins", i, i + 1);
      file.closeGroup();
    }
  }
 
  if (!done) {
    // first check detector_1_events
    file.openGroup("detector_1_events", "NXevent_data");
    entries = file.getEntries();
    for (string_map_t::const_iterator it = entries.begin(); it != entries.end(); ++it) {
      if (it->first == "time_of_flight" || it->first == "event_time_bins") {
        makeTimeOfFlightDataFuzzy(file, localWorkspace, it->first);
        done = true;
      }
    }
    file.closeGroup(); // detector_1_events
 
    if (!done) { // if time_of_flight was not found try
                 // instrument/dae/time_channels_#
      file.openGroup("instrument", "NXinstrument");
      file.openGroup("dae", "IXdae");
      entries = file.getEntries();
      size_t time_channels_number = 0;
      for (string_map_t::const_iterator it = entries.begin(); it != entries.end(); ++it) {
        // check if there are groups with names "time_channels_#" and select the
        // one with the highest number
        if (it->first.size() > 14 && it->first.substr(0, 14) == "time_channels_") {
          size_t n = boost::lexical_cast<size_t>(it->first.substr(14));
          if (n > time_channels_number) {
            time_channels_number = n;
          }
        }
      }
      if (time_channels_number > 0) // the numbers start with 1
      {
        file.openGroup("time_channels_" + std::to_string(time_channels_number), "IXtime_channels");
        entries = file.getEntries();
        for (string_map_t::const_iterator it = entries.begin(); it != entries.end(); ++it) {
          if (it->first == "time_of_flight" || it->first == "event_time_bins") {
            makeTimeOfFlightDataFuzzy(file, localWorkspace, it->first);
          }
        }
        file.closeGroup();
      }
      file.closeGroup(); // dae
      file.closeGroup(); // instrument
    }
  }
 
  // close top entry (or entry given in entry_name)
  file.closeGroup();
}
 
//-----------------------------------------------------------------------------
template <typename T>
bool LoadEventNexus::runLoadInstrument(const std::string &nexusfilename, T localWorkspace,
                                       const std::string &top_entry_name, Algorithm *alg,
                                       const Kernel::NexusHDF5Descriptor *descriptor) {
  std::string instrument;
  std::string instFilename;
 
  const bool isNexus = (descriptor == nullptr) ? LoadGeometry::isNexus(nexusfilename)
                                               : LoadGeometry::isNexus(nexusfilename, descriptor->getAllEntries());
 
  // Check if the geometry can be loaded directly from the Nexus file
  if (isNexus) {
    instFilename = nexusfilename;
  } else {
    // Get the instrument name
    ::NeXus::File nxfile(nexusfilename);
    // Start with the base entry
    nxfile.openGroup(top_entry_name, "NXentry");
    // Open the instrument
    nxfile.openGroup("instrument", "NXinstrument");
    try {
      nxfile.openData("name");
      instrument = nxfile.getStrData();
      alg->getLogger().debug() << "Instrument name read from NeXus file is " << instrument << '\n';
    } catch (::NeXus::Exception &) {
      // Try to fall back to isis compatibility options
      nxfile.closeGroup();
      instrument = readInstrumentFromISIS_VMSCompat(nxfile);
      if (instrument.empty()) {
        // Get the instrument name from the file instead
        size_t n = nexusfilename.rfind('/');
        if (n != std::string::npos) {
          std::string temp = nexusfilename.substr(n + 1, nexusfilename.size() - n - 1);
          n = temp.find('_');
          if (n != std::string::npos && n > 0) {
            instrument = temp.substr(0, n);
          }
        }
      }
    }
    if (instrument == "POWGEN3") // hack for powgen b/c of bad long name
      instrument = "POWGEN";
    if (instrument == "NOM") // hack for nomad
      instrument = "NOMAD";
 
    if (instrument.empty())
      throw std::runtime_error("Could not find the instrument name in the NXS "
                               "file or using the filename. Cannot load "
                               "instrument!");
 
    // Now let's close the file as we don't need it anymore to load the
    // instrument.
    nxfile.close();
  }
 
  if (instFilename.empty()) {
    try {
      instFilename =
          API::InstrumentFileFinder::getInstrumentFilename(instrument, localWorkspace->getWorkspaceStartDate());
    } catch (Kernel::Exception::NotFoundError &) {
      if (instFilename.empty()) {
        Poco::Path directory(Kernel::ConfigService::Instance().getInstrumentDirectory());
        Poco::Path file(instrument + "_Definition.xml");
        Poco::Path fullPath(directory, file);
        instFilename = fullPath.toString();
      }
    }
  }
 
  // do the actual work
  auto loadInst = alg->createChildAlgorithm("LoadInstrument");
 
  // Now execute the Child Algorithm. Catch and log any error, but don't stop.
  bool executionSuccessful(true);
  try {
    loadInst->setPropertyValue("Filename", instFilename);
    loadInst->setPropertyValue("InstrumentName", instrument);
    loadInst->setProperty<Mantid::API::MatrixWorkspace_sptr>("Workspace", localWorkspace);
    loadInst->setProperty("RewriteSpectraMap", Mantid::Kernel::OptionalBool(false));
    loadInst->execute();
 
    // Populate the instrument parameters in this workspace - this works around
    // a bug
    localWorkspace->populateInstrumentParameters();
  } catch (std::invalid_argument &e) {
    alg->getLogger().information() << "Invalid argument to LoadInstrument Child Algorithm : " << e.what() << '\n';
    executionSuccessful = false;
  } catch (std::runtime_error &e) {
    alg->getLogger().information("Unable to successfully run LoadInstrument Child Algorithm");
    alg->getLogger().information(e.what());
    executionSuccessful = false;
  }
 
  // If loading instrument definition file fails
  if (!executionSuccessful) {
    alg->getLogger().error() << "Error loading Instrument definition file\n";
    return false;
  }
 
  // Ticket #2049: Cleanup all loadinstrument members to a single instance
  // If requested update the instrument to positions in the data file
  const auto &pmap = localWorkspace->constInstrumentParameters();
  if (!pmap.contains(localWorkspace->getInstrument()->getComponentID(), "det-pos-source"))
    return executionSuccessful;
 
  std::shared_ptr<Geometry::Parameter> updateDets =
      pmap.get(localWorkspace->getInstrument()->getComponentID(), "det-pos-source");
  std::string value = updateDets->value<std::string>();
  if (value.substr(0, 8) == "datafile") {
    auto updateInst = alg->createChildAlgorithm("UpdateInstrumentFromFile");
    updateInst->setProperty<Mantid::API::MatrixWorkspace_sptr>("Workspace", localWorkspace);
    updateInst->setPropertyValue("Filename", nexusfilename);
    if (value == "datafile-ignore-phi") {
      updateInst->setProperty("IgnorePhi", true);
      alg->getLogger().information("Detector positions in IDF updated with "
                                   "positions in the data file except for the "
                                   "phi values");
    } else {
      alg->getLogger().information("Detector positions in IDF updated with positions in the data file");
    }
    // We want this to throw if it fails to warn the user that the information
    // is not correct.
    updateInst->execute();
  }
 
  return executionSuccessful;
}
 
//-----------------------------------------------------------------------------
template <typename T>
void LoadEventNexus::loadEntryMetadata(const std::string &nexusfilename, T WS, const std::string &entry_name,
                                       const Kernel::NexusHDF5Descriptor &descriptor) {
  // Open the file
  ::NeXus::File file(nexusfilename);
  file.openGroup(entry_name, "NXentry");
 
  // get the title
  if (descriptor.isEntry("/" + entry_name + "/title", "SDS")) {
    file.openData("title");
    if (file.getInfo().type == ::NeXus::CHAR) {
      std::string title = file.getStrData();
      if (!title.empty())
        WS->setTitle(title);
    }
    file.closeData();
  }
 
  // get the notes
  if (descriptor.isEntry("/" + entry_name + "/notes", "SDS")) {
    file.openData("notes");
    if (file.getInfo().type == ::NeXus::CHAR) {
      std::string notes = file.getStrData();
      if (!notes.empty())
        WS->mutableRun().addProperty("file_notes", notes, true);
    }
    file.closeData();
  }
 
  // Get the run number
  if (descriptor.isEntry("/" + entry_name + "/run_number", "SDS")) {
    file.openData("run_number");
    std::string run;
    if (file.getInfo().type == ::NeXus::CHAR) {
      run = file.getStrData();
    } else if (file.isDataInt()) {
      // inside ISIS the run_number type is int32
      std::vector<int> value;
      file.getData(value);
      if (!value.empty())
        run = std::to_string(value[0]);
    }
    if (!run.empty()) {
      WS->mutableRun().addProperty("run_number", run, true);
    }
    file.closeData();
  }
 
  // get the experiment identifier
  if (descriptor.isEntry("/" + entry_name + "/experiment_identifier", "SDS")) {
    file.openData("experiment_identifier");
    std::string expId;
    if (file.getInfo().type == ::NeXus::CHAR) {
      expId = file.getStrData();
    }
    if (!expId.empty()) {
      WS->mutableRun().addProperty("experiment_identifier", expId, true);
    }
    file.closeData();
  }
 
  // get the sample name - nested try/catch to leave the handle in an
  // appropriate state
  if (descriptor.isEntry("/" + entry_name + "/sample", "NXsample")) {
    file.openGroup("sample", "NXsample");
    try {
      if (descriptor.isEntry("/" + entry_name + "/sample/name", "SDS")) {
        file.openData("name");
        const auto info = file.getInfo();
        std::string name;
        if (info.type == ::NeXus::CHAR) {
          if (info.dims.size() == 1) {
            name = file.getStrData();
          } else { // something special for 2-d array
            const int64_t total_length = std::accumulate(info.dims.begin(), info.dims.end(), static_cast<int64_t>(1),
                                                         std::multiplies<int64_t>());
            boost::scoped_array<char> val_array(new char[total_length]);
            file.getData(val_array.get());
            name = std::string(val_array.get(), total_length);
          }
        }
        file.closeData();
        if (!name.empty()) {
          WS->mutableSample().setName(name);
        }
      }
    } catch (::NeXus::Exception &) {
      // let it drop on floor if an exception occurs while reading sample
    }
    file.closeGroup();
  }
 
  // get the duration
  if (descriptor.isEntry("/" + entry_name + "/duration", "SDS")) {
    file.openData("duration");
    std::vector<double> duration;
    file.getDataCoerce(duration);
    if (duration.size() == 1) {
      // get the units
      // clang-format off
    std::vector< ::NeXus::AttrInfo> infos = file.getAttrInfos();
    std::string units;
    for (std::vector< ::NeXus::AttrInfo>::const_iterator it = infos.begin();
         it != infos.end(); ++it) {
      if (it->name == "units") {
        units = file.getStrAttr(*it);
        break;
      }
    }
      // clang-format on
 
      // set the property
      WS->mutableRun().addProperty("duration", duration[0], units, true);
    }
    file.closeData();
  }
 
  // close the file
  file.close();
}
 
//-----------------------------------------------------------------------------
template <typename T>
bool LoadEventNexus::loadInstrument(const std::string &nexusfilename, T localWorkspace,
                                    const std::string &top_entry_name, Algorithm *alg,
                                    const Kernel::NexusHDF5Descriptor *descriptor) {
 
  bool loadNexusInstrumentXML = true;
  if (alg->existsProperty("LoadNexusInstrumentXML"))
    loadNexusInstrumentXML = alg->getProperty("LoadNexusInstrumentXML");
 
  bool foundInstrument = false;
  if (loadNexusInstrumentXML)
    foundInstrument = runLoadIDFFromNexus<T>(nexusfilename, localWorkspace, top_entry_name, alg);
  if (!foundInstrument)
    foundInstrument = runLoadInstrument<T>(nexusfilename, localWorkspace, top_entry_name, alg, descriptor);
  return foundInstrument;
}
 
//-----------------------------------------------------------------------------
template <typename T>
bool LoadEventNexus::runLoadIDFFromNexus(const std::string &nexusfilename, T localWorkspace,
                                         const std::string &top_entry_name, Algorithm *alg) {
  // Test if IDF exists in file, move on quickly if not
  try {
    ::NeXus::File nxsfile(nexusfilename);
    nxsfile.openPath(top_entry_name + "/instrument/instrument_xml");
  } catch (::NeXus::Exception &) {
    alg->getLogger().information("No instrument XML definition found in " + nexusfilename + " at " + top_entry_name +
                                 "/instrument");
    return false;
  }
 
  auto loadInst = alg->createChildAlgorithm("LoadIDFFromNexus");
 
  // Now execute the Child Algorithm. Catch and log any error, but don't stop.
  try {
    loadInst->setPropertyValue("Filename", nexusfilename);
    loadInst->setProperty<Mantid::API::MatrixWorkspace_sptr>("Workspace", localWorkspace);
    loadInst->setPropertyValue("InstrumentParentPath", top_entry_name);
    loadInst->execute();
  } catch (std::invalid_argument &) {
    alg->getLogger().error("Invalid argument to LoadIDFFromNexus Child Algorithm ");
  } catch (std::runtime_error &) {
    alg->getLogger().debug("No instrument definition found by LoadIDFFromNexus in " + nexusfilename + " at " +
                           top_entry_name + "/instrument");
  }
 
  if (!loadInst->isExecuted())
    alg->getLogger().information("No IDF loaded from Nexus file.");
  return loadInst->isExecuted();
}
} // namespace DataHandling
} // namespace Mantid