Skip to content
Snippets Groups Projects
LoadEventNexus.cpp 60.3 KiB
Newer Older
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//     NScD Oak Ridge National Laboratory, European Spallation Source
//     & Institut Laue - Langevin
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidDataHandling/LoadEventNexus.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/Sample.h"
#include "MantidDataHandling/DefaultEventLoader.h"
#include "MantidDataHandling/EventWorkspaceCollection.h"
#include "MantidDataHandling/LoadEventNexusIndexSetup.h"
#include "MantidDataHandling/ParallelEventLoader.h"
#include "MantidGeometry/Instrument/Goniometer.h"
#include "MantidGeometry/Instrument/RectangularDetector.h"
#include "MantidIndexing/IndexInfo.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/DateAndTimeHelpers.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/Timer.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/VisibleWhenProperty.h"

#include <H5Cpp.h>
#include <boost/shared_array.hpp>
#include <boost/shared_ptr.hpp>
using Mantid::Types::Core::DateAndTime;
using std::map;
using std::string;
using std::vector;
using namespace ::NeXus;

namespace Mantid {
namespace DataHandling {
DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadEventNexus)
using namespace Kernel;
using namespace DateAndTimeHelpers;
using namespace Geometry;
using namespace API;
using namespace DataObjects;
using Types::Core::DateAndTime;
namespace {
// detnotes the end of iteration for NeXus::getNextEntry
const std::string NULL_STR("NULL");
} // namespace

/**
 * Based on the current group in the file, does the named sub-entry exist?
 * @param file : File handle. This is not modified, but cannot be const
 * @param name : sub entry name to look for
 * @return true only if it exists
 */
bool exists(::NeXus::File &file, const std::string &name) {
  auto entries = file.getEntries();
  return entries.find(name) != entries.end();
}

//----------------------------------------------------------------------------------------------
/** Empty default constructor
LoadEventNexus::LoadEventNexus()
    : filter_tof_min(0), filter_tof_max(0), m_specMin(0), m_specMax(0),
      longest_tof(0), shortest_tof(0), bad_tofs(0), discarded_events(0),
      compressTolerance(0), m_instrument_loaded_correctly(false),
      loadlogs(false), event_id_is_spec(false) {}
//----------------------------------------------------------------------------------------------
/**
 * Return the confidence with with this algorithm can load the file
 * @param descriptor A descriptor for the file
 * @returns An integer specifying the confidence level. 0 indicates it will not
 * be used
 */
int LoadEventNexus::confidence(Kernel::NexusDescriptor &descriptor) const {
  int confidence(0);
  if (descriptor.classTypeExists("NXevent_data")) {
    if (descriptor.pathOfTypeExists("/entry", "NXentry") ||
        descriptor.pathOfTypeExists("/raw_data_1", "NXentry")) {
      confidence = 80;
    }
  }
  return confidence;
}
//----------------------------------------------------------------------------------------------
/** Initialisation method.
void LoadEventNexus::init() {
  const std::vector<std::string> exts{".nxs.h5", ".nxs", "_event.nxs"};
  this->declareProperty(
Sam Jenkins's avatar
Sam Jenkins committed
      std::make_unique<FileProperty>("Filename", "", FileProperty::Load, exts),
      "The name of the Event NeXus file to read, including its full or "
      "relative path. "
      "The file name is typically of the form INST_####_event.nxs (N.B. case "
      "sensitive if running on Linux).");

  this->declareProperty(
      std::make_unique<WorkspaceProperty<Workspace>>("OutputWorkspace", "",
Sam Jenkins's avatar
Sam Jenkins committed
                                                     Direction::Output),
      "The name of the output EventWorkspace or WorkspaceGroup in which to "
      "load the EventNexus file.");
  declareProperty(
      std::make_unique<PropertyWithValue<string>>("NXentryName", "",
Sam Jenkins's avatar
Sam Jenkins committed
                                                  Direction::Input),
      "Optional: Name of the NXentry to load if it's not the default.");

  declareProperty(std::make_unique<PropertyWithValue<double>>(
                      "FilterByTofMin", EMPTY_DBL(), Direction::Input),
                  "Optional: To exclude events that do not fall within a range "
                  "of times-of-flight. "
                  "This is the minimum accepted value in microseconds. Keep "
                  "blank to load all events.");

  declareProperty(std::make_unique<PropertyWithValue<double>>(
                      "FilterByTofMax", EMPTY_DBL(), Direction::Input),
                  "Optional: To exclude events that do not fall within a range "
                  "of times-of-flight. "
                  "This is the maximum accepted value in microseconds. Keep "
                  "blank to load all events.");

  declareProperty(std::make_unique<PropertyWithValue<double>>(
                      "FilterByTimeStart", EMPTY_DBL(), Direction::Input),
                  "Optional: To only include events after the provided start "
                  "time, in seconds (relative to the start of the run).");

  declareProperty(std::make_unique<PropertyWithValue<double>>(
                      "FilterByTimeStop", EMPTY_DBL(), Direction::Input),
                  "Optional: To only include events before the provided stop "
                  "time, in seconds (relative to the start of the run).");

  std::string grp1 = "Filter Events";
  setPropertyGroup("FilterByTofMin", grp1);
  setPropertyGroup("FilterByTofMax", grp1);
  setPropertyGroup("FilterByTimeStart", grp1);
  setPropertyGroup("FilterByTimeStop", grp1);

      std::make_unique<ArrayProperty<string>>("BankName", Direction::Input),
      "Optional: To only include events from one bank. Any bank "
      "whose name does not match the given string will have no "
      "events.");
Sam Jenkins's avatar
Sam Jenkins committed
  declareProperty(std::make_unique<PropertyWithValue<bool>>(
                      "SingleBankPixelsOnly", true, Direction::Input),
                  "Optional: Only applies if you specified a single bank to "
                  "load with BankName. "
                  "Only pixels in the specified bank will be created if true; "
                  "all of the instrument's pixels will be created otherwise.");
Sam Jenkins's avatar
Sam Jenkins committed
  setPropertySettings(
      "SingleBankPixelsOnly",
      std::make_unique<VisibleWhenProperty>("BankName", IS_NOT_DEFAULT));

  std::string grp2 = "Loading a Single Bank";
  setPropertyGroup("BankName", grp2);
  setPropertyGroup("SingleBankPixelsOnly", grp2);

  declareProperty(
Sam Jenkins's avatar
Sam Jenkins committed
      std::make_unique<PropertyWithValue<bool>>("Precount", true,
                                                Direction::Input),
      "Pre-count the number of events in each pixel before allocating memory "
      "(optional, default True). "
      "This can significantly reduce memory use and memory fragmentation; it "
      "may also speed up loading.");

  declareProperty(std::make_unique<PropertyWithValue<double>>(
                      "CompressTolerance", -1.0, Direction::Input),
                  "Run CompressEvents while loading (optional, leave blank or "
                  "negative to not do). "
                  "This specified the tolerance to use (in microseconds) when "
                  "compressing.");

  auto mustBePositive = boost::make_shared<BoundedValidator<int>>();
  mustBePositive->setLower(1);
  declareProperty("ChunkNumber", EMPTY_INT(), mustBePositive,
                  "If loading the file by sections ('chunks'), this is the "
                  "section number of this execution of the algorithm.");
  declareProperty("TotalChunks", EMPTY_INT(), mustBePositive,
                  "If loading the file by sections ('chunks'), this is the "
                  "total number of sections.");
  // TotalChunks is only meaningful if ChunkNumber is set
  // Would be nice to be able to restrict ChunkNumber to be <= TotalChunks at
  // validation
  setPropertySettings("TotalChunks", std::make_unique<VisibleWhenProperty>(
                                         "ChunkNumber", IS_NOT_DEFAULT));

  std::string grp3 = "Reduce Memory Use";
  setPropertyGroup("Precount", grp3);
  setPropertyGroup("CompressTolerance", grp3);
  setPropertyGroup("ChunkNumber", grp3);
  setPropertyGroup("TotalChunks", grp3);

Sam Jenkins's avatar
Sam Jenkins committed
  declareProperty(std::make_unique<PropertyWithValue<bool>>(
                      "LoadMonitors", false, Direction::Input),
                  "Load the monitors from the file (optional, default False).");
  std::vector<std::string> options{"", "Events", "Histogram"};
  declareProperty("MonitorsLoadOnly", "",
                  boost::make_shared<Kernel::StringListValidator>(options),
                  "If multiple repesentations exist, which one to load. "
                  "Default is to load the one that is present.");
  declareProperty(std::make_unique<PropertyWithValue<double>>(
                      "FilterMonByTofMin", EMPTY_DBL(), Direction::Input),
                  "Optional: To exclude events from monitors that do not fall "
                  "within a range of times-of-flight. "
                  "This is the minimum accepted value in microseconds.");

  declareProperty(std::make_unique<PropertyWithValue<double>>(
                      "FilterMonByTofMax", EMPTY_DBL(), Direction::Input),
                  "Optional: To exclude events from monitors that do not fall "
                  "within a range of times-of-flight. "
                  "This is the maximum accepted value in microseconds.");

  declareProperty(std::make_unique<PropertyWithValue<double>>(
                      "FilterMonByTimeStart", EMPTY_DBL(), Direction::Input),
                  "Optional: To only include events from monitors after the "
                  "provided start time, in seconds (relative to the start of "
                  "the run).");

  declareProperty(std::make_unique<PropertyWithValue<double>>(
                      "FilterMonByTimeStop", EMPTY_DBL(), Direction::Input),
                  "Optional: To only include events from monitors before the "
                  "provided stop time, in seconds (relative to the start of "
                  "the run).");

  setPropertySettings(
      "MonitorsLoadOnly",
      std::make_unique<VisibleWhenProperty>("LoadMonitors", IS_EQUAL_TO, "1"));
  auto asEventsIsOn = [] {
    std::unique_ptr<IPropertySettings> prop =
        std::make_unique<VisibleWhenProperty>("LoadMonitors", IS_EQUAL_TO, "1");
    return prop;
  };
  setPropertySettings("FilterMonByTofMin", asEventsIsOn());
  setPropertySettings("FilterMonByTofMax", asEventsIsOn());
  setPropertySettings("FilterMonByTimeStart", asEventsIsOn());
  setPropertySettings("FilterMonByTimeStop", asEventsIsOn());

  std::string grp4 = "Monitors";
  setPropertyGroup("LoadMonitors", grp4);
  setPropertyGroup("MonitorsLoadOnly", grp4);
  setPropertyGroup("FilterMonByTofMin", grp4);
  setPropertyGroup("FilterMonByTofMax", grp4);
  setPropertyGroup("FilterMonByTimeStart", grp4);
  setPropertyGroup("FilterMonByTimeStop", grp4);

Hahn, Steven's avatar
Hahn, Steven committed
  declareProperty("SpectrumMin", EMPTY_INT(), mustBePositive,
                  "The number of the first spectrum to read.");
Hahn, Steven's avatar
Hahn, Steven committed
  declareProperty("SpectrumMax", EMPTY_INT(), mustBePositive,
                  "The number of the last spectrum to read.");
  declareProperty(std::make_unique<ArrayProperty<int32_t>>("SpectrumList"),
                  "A comma-separated list of individual spectra to read.");

  declareProperty(
      std::make_unique<PropertyWithValue<bool>>("MetaDataOnly", false,
Sam Jenkins's avatar
Sam Jenkins committed
                                                Direction::Input),
      "If true, only the meta data and sample logs will be loaded.");

Sam Jenkins's avatar
Sam Jenkins committed
  declareProperty(std::make_unique<PropertyWithValue<bool>>("LoadLogs", true,
                                                            Direction::Input),
                  "Load the Sample/DAS logs from the file (default True).");
  std::vector<std::string> loadType{"Default"};
  loadType.emplace_back("Multiprocess (experimental)");
  loadType.emplace_back("MPI");
  auto loadTypeValidator = boost::make_shared<StringListValidator>(loadType);
  declareProperty("LoadType", "Default", loadTypeValidator,
                  "Set type of loader. 2 options {Default, Multiproceess},"
                  "'Multiprocess' should work faster for big files and it is "
                  "experimental, available only in Linux");
Sam Jenkins's avatar
Sam Jenkins committed
  declareProperty(std::make_unique<PropertyWithValue<bool>>(
                      "LoadNexusInstrumentXML", true, Direction::Input),
                  "Reads the embedded Instrument XML from the NeXus file "
                  "(optional, default True). ");
//----------------------------------------------------------------------------------------------
/** set the name of the top level NXentry m_top_entry_name
void LoadEventNexus::setTopEntryName() {
  std::string nxentryProperty = getProperty("NXentryName");
  if (!nxentryProperty.empty()) {
    m_top_entry_name = nxentryProperty;
    return;
  }

  try {
    while (true) {
      const auto entry = m_file->getNextEntry();
      if (entry.second == "NXentry") {
        if ((entry.first == "entry") || (entry.first == "raw_data_1")) {
          m_top_entry_name = entry.first;
          break;
        }
      } else if (entry.first == NULL_STR && entry.second == NULL_STR) {
        g_log.error()
            << "Unable to determine name of top level NXentry - assuming "
               "\"entry\".\n";
        m_top_entry_name = "entry";
    }
  } catch (const std::exception &) {
Hahn, Steven's avatar
Hahn, Steven committed
    g_log.error() << "Unable to determine name of top level NXentry - assuming "
                     "\"entry\".\n";
    m_top_entry_name = "entry";
  }
}
template <typename T> void LoadEventNexus::filterDuringPause(T workspace) {
  try {
    if ((!ConfigService::Instance().hasProperty(
            "loadeventnexus.keeppausedevents")) &&
        (m_ws->run().getLogData("pause")->size() > 1)) {
      g_log.notice("Filtering out events when the run was marked as paused. "
                   "Set the loadeventnexus.keeppausedevents configuration "
                   "property to override this.");

      auto filter = createChildAlgorithm("FilterByLogValue");
      filter->setProperty("InputWorkspace", workspace);
      filter->setProperty("OutputWorkspace", workspace);
      filter->setProperty("LogName", "pause");
      // The log value is set to 1 when the run is paused, 0 otherwise.
      filter->setProperty("MinimumValue", 0.0);
      filter->setProperty("MaximumValue", 0.0);
      filter->setProperty("LogBoundary", "Left");
      filter->execute();
    }
  } catch (Exception::NotFoundError &) {
    // No "pause" log, just carry on
  }
}

template <>
void LoadEventNexus::filterDuringPause<EventWorkspaceCollection_sptr>(
    EventWorkspaceCollection_sptr workspace) {
  // We provide a function pointer to the filter method of the object
  using std::placeholders::_1;
  auto func = std::bind(
      &LoadEventNexus::filterDuringPause<MatrixWorkspace_sptr>, this, _1);
  workspace->applyFilter(func);
}

//------------------------------------------------------------------------------------------------
/** Executes the algorithm. Reading in the file and creating and populating
 *  the output workspace
 */
void LoadEventNexus::exec() {
  // Retrieve the filename from the properties
  m_filename = getPropertyValue("Filename");

  compressTolerance = getProperty("CompressTolerance");

  loadlogs = getProperty("LoadLogs");

  // Check to see if the monitors need to be loaded later
  bool load_monitors = this->getProperty("LoadMonitors");

  // this must make absolutely sure that m_file is a valid (and open)
  // NeXus::File object
  safeOpenFile(m_filename);

  setTopEntryName();

  // Initialize progress reporting.
  int reports = 3;
  if (load_monitors)
    reports++;
  Progress prog(this, 0.0, 0.3, reports);

  // Load the detector events
  m_ws = boost::make_shared<EventWorkspaceCollection>(); // Algorithm currently
                                                         // relies on an
  // object-level workspace ptr
  loadEvents(&prog, false); // Do not load monitor blocks

  if (discarded_events > 0) {
    g_log.information() << discarded_events
                        << " events were encountered coming from pixels which "
                           "are not in the Instrument Definition File."
                           "These events were discarded.\n";
  }

  // If the run was paused at any point, filter out those events (SNS only, I
  // think)
  filterDuringPause(m_ws->getSingleHeldWorkspace());
  m_ws->mutableRun().addProperty("Filename", m_filename);
  // Save output
  this->setProperty("OutputWorkspace", m_ws->combinedWorkspace());
Peterson, Peter's avatar
Peterson, Peter committed

  // close the file since LoadNexusMonitors will take care of its own file
  // handle
  m_file->close();

  // Load the monitors with child algorithm 'LoadNexusMonitors'
  if (load_monitors) {
    prog.report("Loading monitors");
Peterson, Peter's avatar
Peterson, Peter committed
    this->runLoadMonitors();
std::pair<DateAndTime, DateAndTime>
firstLastPulseTimes(::NeXus::File &file, Kernel::Logger &logger) {
  file.openData("event_time_zero");
  if (!file.hasAttr("offset"))
    throw std::runtime_error("No ISO8601 offset attribute provided");

  const auto heldTimeZeroType = file.getInfo().type;
  // Nexus only requires event_time_zero to be a NXNumber, we support two
  // possibilities for held type

  std::string isooffset; // ISO8601 offset
  file.getAttr("offset", isooffset);
  DateAndTime offset(isooffset);
  std::string units; // time units
  if (file.hasAttr("units"))
    file.getAttr("units", units);
  // TODO. Logic here is similar to BankPulseTimes (ctor) should be consolidated
  if (heldTimeZeroType == ::NeXus::UINT64) {
      logger.warning(
          "event_time_zero is uint64_t, but units not in ns. Found to be: " +
    std::vector<uint64_t> nanoseconds;
    file.readData("event_time_zero", nanoseconds);
    if (nanoseconds.empty())
      throw std::runtime_error(
          "No event time zeros. Cannot establish run start or end");
    auto absoluteFirst = DateAndTime(int64_t(0), int64_t(nanoseconds.front())) +
                         offset.totalNanoseconds();
    auto absoluteLast = DateAndTime(int64_t(0), int64_t(nanoseconds.back())) +
                        offset.totalNanoseconds();
    return std::make_pair(absoluteFirst, absoluteLast);
  } else if (heldTimeZeroType == ::NeXus::FLOAT64) {
    if (units != "second")
      logger.warning("event_time_zero is double_t, but units not in seconds. "
                     "Found to be: " +
    std::vector<double> seconds;
    file.readData("event_time_zero", seconds);
    if (seconds.empty())
      throw std::runtime_error(
          "No event time zeros. Cannot establish run start or end");
    auto absoluteFirst =
        DateAndTime(seconds.front(), double(0)) + offset.totalNanoseconds();
    auto absoluteLast =
        DateAndTime(seconds.back(), double(0)) + offset.totalNanoseconds();
    return std::make_pair(absoluteFirst, absoluteLast);
  }

  else {
    throw std::runtime_error("Unrecognised type for event_time_zero");
  }
} // namespace DataHandling
 * Get the number of events in the currently opened group.
 *
 * @param file The handle to the nexus file opened to the group to look at.
 * @param hasTotalCounts Whether to try looking at the total_counts field.
 * This variable will be changed if the field is not there.
 * @param oldNeXusFileNames Whether to try using old names. This variable will
 * be changed if it is determined that old names are being used.
 *
 * @return The number of events.
 */
std::size_t numEvents(::NeXus::File &file, bool &hasTotalCounts,
                      bool &oldNeXusFileNames) {
  // try getting the value of total_counts
  if (hasTotalCounts) {
    hasTotalCounts = false;
    if (exists(file, "total_counts")) {
      try {
        file.openData("total_counts");
        auto info = file.getInfo();
        file.closeData();
        if (info.type == NeXus::UINT64) {
          uint64_t eventCount;
          file.readData("total_counts", eventCount);
          hasTotalCounts = true;
          return eventCount;
        }
      } catch (::NeXus::Exception &) {
      }
    }
  }

  // just get the length of the event pixel ids
  try {
    if (oldNeXusFileNames)
      file.openData("event_pixel_id");
    else
      file.openData("event_id");
  } catch (::NeXus::Exception &) {
    // Older files (before Nov 5, 2010) used this field.
    try {
      file.openData("event_pixel_id");
      oldNeXusFileNames = true;
    } catch (::NeXus::Exception &) {
      // Some groups have neither indicating there are not events here
      return 0;
    }
  }
  size_t numEvents = static_cast<std::size_t>(file.getInfo().dims[0]);
  file.closeData();
  return numEvents;
}
/** Load the instrument from the nexus file
 *
 * @param nexusfilename :: The name of the nexus file being loaded
 * @param localWorkspace :: Templated workspace in which to put the instrument
 *geometry
 * @param alg :: Handle of the algorithm
 * @param returnpulsetimes :: flag to return shared pointer for
 *BankPulseTimes, otherwise NULL.
 * @param nPeriods : Number of periods (write to)
 * @param periodLog : Period logs DateAndTime to int map.
 *
 * @return Pulse times given in the DAS logs
 */
template <typename T>
boost::shared_ptr<BankPulseTimes> LoadEventNexus::runLoadNexusLogs(
    const std::string &nexusfilename, T localWorkspace, API::Algorithm &alg,
    bool returnpulsetimes, int &nPeriods,
Tom Titcombe's avatar
Tom Titcombe committed
    std::unique_ptr<const TimeSeriesProperty<int>> &periodLog) {
  // --------------------- Load DAS Logs -----------------
  // The pulse times will be empty if not specified in the DAS logs.
  // BankPulseTimes * out = NULL;
  boost::shared_ptr<BankPulseTimes> out;
  API::IAlgorithm_sptr loadLogs = alg.createChildAlgorithm("LoadNexusLogs");

  // Now execute the Child Algorithm. Catch and log any error, but don't stop.
  try {
    alg.getLogger().information() << "Loading logs from NeXus file..."
                                  << "\n";
    loadLogs->setPropertyValue("Filename", nexusfilename);
    loadLogs->setProperty<API::MatrixWorkspace_sptr>("Workspace",
                                                     localWorkspace);
    try {
      loadLogs->setPropertyValue("NXentryName",
                                 alg.getPropertyValue("NXentryName"));
    } catch (...) {
    const Run &run = localWorkspace->run();
    // Get the number of periods
    if (run.hasProperty("nperiods")) {
      nPeriods = run.getPropertyValueAsType<int>("nperiods");
    }
    // Get the period log. Map of DateAndTime to Period int values.
    if (run.hasProperty("period_log")) {
      auto *temp = run.getProperty("period_log");
Tom Titcombe's avatar
Tom Titcombe committed
      // Check for corrupted period logs
      std::unique_ptr<TimeSeriesProperty<int>> tempPeriodLog(
          dynamic_cast<TimeSeriesProperty<int> *>(temp->clone()));
      checkForCorruptedPeriods(std::move(tempPeriodLog), periodLog, nPeriods,
                               nexusfilename);
    }

    // If successful, we can try to load the pulse times
    std::vector<Types::Core::DateAndTime> temp;
    if (localWorkspace->run().hasProperty("proton_charge")) {
      auto *log = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(
          localWorkspace->mutableRun().getProperty("proton_charge"));
      if (log)
        temp = log->timesAsVector();
    }
    if (returnpulsetimes)
      out = boost::make_shared<BankPulseTimes>(temp);

    // Use the first pulse as the run_start time.
    if (!temp.empty()) {
      if (temp[0] < Types::Core::DateAndTime("1991-01-01T00:00:00"))
        alg.getLogger().warning() << "Found entries in the proton_charge "
                                     "sample log with invalid pulse time!\n";

      Types::Core::DateAndTime run_start = localWorkspace->getFirstPulseTime();
      // add the start of the run as a ISO8601 date/time string. The start =
      // first non-zero time.
      // (this is used in LoadInstrument to find the right instrument file to
      // use).
      localWorkspace->mutableRun().addProperty(
          "run_start", run_start.toISO8601String(), true);
    } else {
      alg.getLogger().warning() << "Empty proton_charge sample log. You will "
                                   "not be able to filter by time.\n";
    }
    /// Attempt to make a gonoimeter from the logs
    try {
      Geometry::Goniometer gm;
      gm.makeUniversalGoniometer();
      localWorkspace->mutableRun().setGoniometer(gm, true);
    } catch (std::runtime_error &) {
    }
  } catch (const InvalidLogPeriods &) {
    // Rethrow so LoadEventNexus fails.
    // If we don't, Mantid will crash.
  } catch (...) {
    alg.getLogger().error() << "Error while loading Logs from SNS Nexus. Some "
                               "sample logs may be missing."
                            << "\n";
    return out;
  }
  return out;
}

Tom Titcombe's avatar
Tom Titcombe committed
/** Check for corrupted period logs
 * If data is historical (1 periods, period is labelled 0) then change period
 * labels to 1 If number of periods does not match expected number of periods
 * then throw an error
 * @param tempPeriodLog :: a temporary local copy of period logs, which will
 * change
 * @param periodLog :: unique pointer which will point to period logs once they
 * have been changed
 * @param nPeriods :: the value in the nperiods log of the run. Number of
 * expected periods
 * @param nexusfilename :: the filename of the run to load
Tom Titcombe's avatar
Tom Titcombe committed
 */
void LoadEventNexus::checkForCorruptedPeriods(
    std::unique_ptr<TimeSeriesProperty<int>> tempPeriodLog,
    std::unique_ptr<const TimeSeriesProperty<int>> &periodLog,
    const int &nPeriods, const std::string &nexusfilename) {
  const auto valuesAsVector = tempPeriodLog->valuesAsVector();
  const auto nPeriodsInLog =
      *std::max_element(valuesAsVector.begin(), valuesAsVector.end());

  // Check for historic files
  if (nPeriodsInLog == 0 && nPeriods == 1) {
    // "modernize" the local copy here by making period_log
    // a vector of 1s
    const std::vector<int> newValues(tempPeriodLog->realSize(), 1);
    const auto times = tempPeriodLog->timesAsVector();
    periodLog.reset(
        new const TimeSeriesProperty<int>("period_log", times, newValues));
  } else if (nPeriodsInLog != nPeriods) {
    // Sanity check here that period_log only contains period numbers up to
    // nperiods. These values can be different due to instrument noise, and
    // cause undescriptive crashes if not caught.
    // We throw here to make it clear
    // that the file is corrupted and must be manually assessed.
    const auto msg = "File " + nexusfilename +
                     " has been corrupted. The log framelog/period_log/value "
                     "contains " +
                     std::to_string(nPeriodsInLog) +
                     " periods, but periods/number contains " +
                     std::to_string(nPeriods) +
                     ". This file should be manually inspected and corrected.";
    throw InvalidLogPeriods(msg);
  } else {
    // periodLog should point to a copy of the period logs
    periodLog = std::make_unique<const TimeSeriesProperty<int>>(*tempPeriodLog);
    tempPeriodLog.reset();
  }
}

/** Load the instrument from the nexus file
 *
 * @param nexusfilename :: The name of the nexus file being loaded
 * @param localWorkspace :: EventWorkspaceCollection in which to put the
 *instrument
 *geometry
 * @param alg :: Handle of the algorithm
 * @param returnpulsetimes :: flag to return shared pointer for
 *BankPulseTimes, otherwise NULL.
 * @param nPeriods : Number of periods (write to)
 * @param periodLog : Period logs DateAndTime to int map.
 *
 * @return Pulse times given in the DAS logs
 */
boost::shared_ptr<BankPulseTimes>
LoadEventNexus::runLoadNexusLogs<EventWorkspaceCollection_sptr>(
    const std::string &nexusfilename,
    EventWorkspaceCollection_sptr localWorkspace, API::Algorithm &alg,
    bool returnpulsetimes, int &nPeriods,
Tom Titcombe's avatar
Tom Titcombe committed
    std::unique_ptr<const TimeSeriesProperty<int>> &periodLog) {
  auto ws = localWorkspace->getSingleHeldWorkspace();
  auto ret = runLoadNexusLogs<MatrixWorkspace_sptr>(
      nexusfilename, ws, alg, returnpulsetimes, nPeriods, periodLog);
enum class LoadEventNexus::LoaderType { MPI, MULTIPROCESS, DEFAULT };
//-----------------------------------------------------------------------------
/**
 * Load events from the file.
 * @param prog :: A pointer to the progress reporting object
 * @param monitors :: If true the events from the monitors are loaded and not
 *the main banks
 *
 * This also loads the instrument, but only if it has not been set in the
 *workspace
 * being used as input (m_ws data member). Same applies to the logs.
 */
void LoadEventNexus::loadEvents(API::Progress *const prog,
                                const bool monitors) {
  bool metaDataOnly = getProperty("MetaDataOnly");

  // Get the time filters
  setTimeFilters(monitors);

  // The run_start will be loaded from the pulse times.
  DateAndTime run_start(0, 0);
  bool takeTimesFromEvents = false;
  // Initialize the counter of bad TOFs
  bad_tofs = 0;
Tom Titcombe's avatar
Tom Titcombe committed
  auto periodLog =
      std::make_unique<const TimeSeriesProperty<int>>("period_log");
  if (loadlogs) {
    prog->doReport("Loading DAS logs");
    m_allBanksPulseTimes = runLoadNexusLogs<EventWorkspaceCollection_sptr>(
Tom Titcombe's avatar
Tom Titcombe committed
        m_filename, m_ws, *this, true, nPeriods, periodLog);
      run_start = m_ws->getFirstPulseTime();
    } catch (Kernel::Exception::NotFoundError &) {
      /*
        This is added to (a) support legacy behaviour of continuing to take
        times from the proto_charge log, but (b) allowing a fall back of
        getting run start and end from actual pulse times within the
        NXevent_data group. Note that the latter is better Nexus compliant.
      */
      takeTimesFromEvents = true;
  } else {
    g_log.information() << "Skipping the loading of sample logs!\n"
                        << "Reading the start time directly from /"
                        << m_top_entry_name << "/start_time\n";
    // start_time is read and set
    m_file->openPath("/");
    m_file->openGroup(m_top_entry_name, "NXentry");
    std::string tmp;
    m_file->readData("start_time", tmp);
    m_file->closeGroup();
    run_start = createFromSanitizedISO8601(tmp);
    m_ws->mutableRun().addProperty("run_start", run_start.toISO8601String(),
                                   true);
  // set more properties on the workspace
  try {
    // this is a static method that is why it is passing the
    // file object and the file path
    loadEntryMetadata<EventWorkspaceCollection_sptr>(m_filename, m_ws,
                                                     m_top_entry_name);
  } catch (std::runtime_error &e) {
    // Missing metadata is not a fatal error. Log and go on with your life
    g_log.error() << "Error loading metadata: " << e.what() << '\n';
  }
  m_ws->setNPeriods(
      static_cast<size_t>(nPeriods),
      periodLog); // This is how many workspaces we are going to make.
  // Make sure you have a non-NULL m_allBanksPulseTimes
  if (m_allBanksPulseTimes == nullptr) {
    std::vector<DateAndTime> temp;
    m_allBanksPulseTimes = boost::make_shared<BankPulseTimes>(temp);
  }

  if (!m_ws->getInstrument() || !m_instrument_loaded_correctly) {
    // Load the instrument (if not loaded before)
    // Note that closing an re-opening the file is needed here for loading
    // instruments directly from the nexus file containing the event data.
    // This may not be needed in the future if both LoadEventNexus and
    // LoadInstrument are made to use the same Nexus/HDF5 library
    m_file->close();
        loadInstrument(m_filename, m_ws, m_top_entry_name, this);
    if (!m_instrument_loaded_correctly)
      throw std::runtime_error("Instrument was not initialized correctly! "
                               "Loading cannot continue.");
    // reopen file
    safeOpenFile(m_filename);

  // top level file information
  // Start with the base entry
  m_file->openGroup(m_top_entry_name, "NXentry");

  // Now we want to go through all the bankN_event entries
  vector<string> bankNames;
  vector<std::size_t> bankNumEvents;
  std::string classType = monitors ? "NXmonitor" : "NXevent_data";
  ::NeXus::Info info;
  bool oldNeXusFileNames(false);
  bool hasTotalCounts(true);
  auto firstPulseT = DateAndTime::maximum();
  while (true) { // should be broken when entry name is set
    const auto entry = m_file->getNextEntry();
    const std::string entry_name(entry.first);
    const std::string entry_class(entry.second);
    if (entry_name == NULL_STR && entry_class == NULL_STR)
      break;
    if (entry_class == classType) {
      // open the group
      m_file->openGroup(entry_name, classType);
      if (takeTimesFromEvents) {
        /* If we are here, we are loading logs, but have failed to establish
         * the run_start from the proton_charge log. We are going to get this
         * from our event_time_zero instead
        auto localFirstLast = firstLastPulseTimes(*m_file, this->g_log);
        firstPulseT = std::min(firstPulseT, localFirstLast.first);
      }
      // get the number of events
      std::size_t num = numEvents(*m_file, hasTotalCounts, oldNeXusFileNames);
      bankNames.emplace_back(entry_name);
      bankNumEvents.emplace_back(num);

      // Look for weights in simulated file
      haveWeights = exists(*m_file, "event_weight");
  if (takeTimesFromEvents)
    run_start = firstPulseT;
  loadSampleDataISIScompatibility(*m_file, *m_ws);
  // Close the 'top entry' group (raw_data_1 for NexusProcessed, etc.)
  m_file->closeGroup();

  // Delete the output workspace name if it existed
  std::string outName = getPropertyValue("OutputWorkspace");
  if (AnalysisDataService::Instance().doesExist(outName))
    AnalysisDataService::Instance().remove(outName);

  // --------------------------- Time filtering
  // ------------------------------------
  double filter_time_start_sec, filter_time_stop_sec;
  filter_time_start_sec = getProperty("FilterByTimeStart");
  filter_time_stop_sec = getProperty("FilterByTimeStop");

  // Default to ALL pulse times
  bool is_time_filtered = false;
  filter_time_start = Types::Core::DateAndTime::minimum();
  filter_time_stop = Types::Core::DateAndTime::maximum();

  if (m_allBanksPulseTimes->numPulses > 0) {
    // If not specified, use the limits of doubles. Otherwise, convert from
    // seconds to absolute PulseTime
    if (filter_time_start_sec != EMPTY_DBL()) {
      filter_time_start = run_start + filter_time_start_sec;
      is_time_filtered = true;
    if (filter_time_stop_sec != EMPTY_DBL()) {
      filter_time_stop = run_start + filter_time_stop_sec;
      is_time_filtered = true;
    // Silly values?
    if (filter_time_stop < filter_time_start) {
      std::string msg = "Your ";
      if (monitors)
        msg += "monitor ";
      msg += "filter for time's Stop value is smaller than the Start value.";
      throw std::invalid_argument(msg);
  }

  if (is_time_filtered) {
    // Now filter out the run, using the DateAndTime type.
    m_ws->mutableRun().filterByTime(filter_time_start, filter_time_stop);
  }

  if (metaDataOnly) {
    // Now, create a default X-vector for histogramming, with just 2 bins.
    auto axis = HistogramData::BinEdges{
        1, static_cast<double>(std::numeric_limits<uint32_t>::max()) * 0.1 - 1};
    // Set the binning axis using this.
    m_ws->setAllX(axis);
    createSpectraMapping(m_filename, monitors, std::vector<std::string>());
  }

  // --------- Loading only one bank ? ----------------------------------
  std::vector<std::string> someBanks = getProperty("BankName");
  const bool SingleBankPixelsOnly = getProperty("SingleBankPixelsOnly");
  if ((!someBanks.empty()) && (!monitors)) {
    std::vector<std::string> eventedBanks;
    eventedBanks.reserve(someBanks.size());
    for (const auto &bank : someBanks) {
      eventedBanks.emplace_back(bank + "_events");
    }
    // check that all of the requested banks are in the file
    const auto invalidBank =
        std::find_if(eventedBanks.cbegin(), eventedBanks.cend(),
                     [&bankNames](const auto &someBank) {
                       return std::none_of(bankNames.cbegin(), bankNames.cend(),
                                           [&someBank](const auto &name) {
                                             return name == someBank;
                                           });
                     });
    if (invalidBank != eventedBanks.cend()) {
      throw std::invalid_argument("No entry named '" + *invalidBank +
                                  "' was found in the .NXS file.");
    // change the number of banks to load
    bankNames.assign(eventedBanks.cbegin(), eventedBanks.cend());
    // TODO this equally weights the banks
    bankNumEvents.assign(someBanks.size(), 1);

    if (!SingleBankPixelsOnly)
      someBanks.clear(); // Marker to load all pixels
  }

  prog->report("Initializing all pixels");

  // Remove unused banks if parameter is set
  if (m_ws->getInstrument()->hasParameter("remove-unused-banks")) {
    std::vector<double> instrumentUnused =
        m_ws->getInstrument()->getNumberParameter("remove-unused-banks", true);
    if (!instrumentUnused.empty()) {
      const auto unused = static_cast<int>(instrumentUnused.front());
      if (unused == 1)
        deleteBanks(m_ws, bankNames);
  }
  //----------------- Pad Empty Pixels -------------------------------
  createSpectraMapping(m_filename, monitors, someBanks);

  // Set all (empty) event lists as sorted by pulse time. That way, calling
  // SortEvents will not try to sort these empty lists.
  for (size_t i = 0; i < m_ws->getNumberHistograms(); i++)
    m_ws->getSpectrum(i).setSortOrder(DataObjects::PULSETIME_SORT);

  // Count the limits to time of flight
  shortest_tof =
      static_cast<double>(std::numeric_limits<uint32_t>::max()) * 0.1;
  longest_tof = 0.;
  bool loaded{false};
  auto loaderType = defineLoaderType(haveWeights, oldNeXusFileNames, classType);
  if (loaderType != LoaderType::DEFAULT) {
    auto ws = m_ws->getSingleHeldWorkspace();
    m_file->close();
    if (loaderType == LoaderType::MPI) {
      try {
        ParallelEventLoader::loadMPI(*ws, m_filename, m_top_entry_name,
                                     bankNames, event_id_is_spec);
        g_log.information() << "Used MPI ParallelEventLoader.\n";
        loaded = true;
        shortest_tof = 0.0;
        longest_tof = 1e10;
      } catch (const std::runtime_error &) {
        g_log.warning()
            << "MPI event loader failed, falling back to default loader.\n";
      }
    } else {

      struct ExceptionOutput {
        static void out(decltype(g_log) &log, const std::exception &e,
                        int level = 0) {
          log.warning() << std::string(level, ' ') << "exception: " << e.what()
                        << '\n';
          try {
            std::rethrow_if_nested(e);