LoadNexusProcessed.cpp

// 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 & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidDataHandling/LoadNexusProcessed.h"
#include "MantidAPI/AlgorithmFactory.h"
#include "MantidAPI/AlgorithmManager.h"
#include "MantidAPI/BinEdgeAxis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/TextAxis.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidAPI/WorkspaceHistory.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/Peak.h"
#include "MantidDataObjects/PeakNoShapeFactory.h"
#include "MantidDataObjects/PeakShapeEllipsoidFactory.h"
#include "MantidDataObjects/PeakShapeSphericalFactory.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidDataObjects/RebinnedOutput.h"
#include "MantidGeometry/Instrument/Goniometer.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/DateAndTime.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/StringTokenizer.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidNexus/NexusClasses.h"
#include "MantidNexus/NexusFileIO.h"

#include <boost/regex.hpp>
#include <boost/shared_array.hpp>
#include <boost/shared_ptr.hpp>

#include <nexus/NeXusException.hpp>

#include <map>
#include <string>
#include <vector>

namespace Mantid {
namespace DataHandling {

// Register the algorithm into the algorithm factory
DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadNexusProcessed)

using namespace Mantid::NeXus;
using namespace DataObjects;
using namespace Kernel;
using namespace API;
using Geometry::Instrument_const_sptr;
using Mantid::Types::Event::TofEvent;
using Types::Core::DateAndTime;

namespace {

// Helper typedef
using IntArray_shared = boost::shared_array<int>;

// Struct to contain spectrum information.
struct SpectraInfo {
  // Number of spectra
  int nSpectra{0};
  // Do we have any spectra
  bool hasSpectra{false};
  // Contains spectrum numbers for each workspace index
  IntArray_shared spectraNumbers;
  // Index of the detector in the workspace.
  IntArray_shared detectorIndex;
  // Number of detectors associated with each spectra
  IntArray_shared detectorCount;
  // Detector list contains a list of all of the detector numbers
  IntArray_shared detectorList;
};

// Helper typdef.
using SpectraInfo_optional = boost::optional<SpectraInfo>;

/**
 * Extract ALL the detector, spectrum number and workspace index mapping
 * information.
 * @param mtd_entry
 * @param logger
 * @return
 */
SpectraInfo extractMappingInfo(NXEntry &mtd_entry, Logger &logger) {
  SpectraInfo spectraInfo;
  // Instrument information

  if (!mtd_entry.containsGroup("instrument")) {
    logger.information() << "No NXinstrument group called `instrument` under "
                            "NXEntry. The workspace will not "
                            "contain any detector information.\n";
    return spectraInfo;
  }
  NXInstrument inst = mtd_entry.openNXInstrument("instrument");
  if (!inst.containsGroup("detector")) {
    logger.information() << "Detector block not found. The workspace will not "
                            "contain any detector information.\n";
    return spectraInfo;
  }

  // Populate the spectra-detector map
  NXDetector detgroup = inst.openNXDetector("detector");

  // Spectra block - Contains spectrum numbers for each workspace index
  // This might not exist so wrap and check. If it doesn't exist create a
  // default mapping
  try {
    NXInt spectra_block = detgroup.openNXInt("spectra");
    spectra_block.load();
    spectraInfo.spectraNumbers = spectra_block.sharedBuffer();
    spectraInfo.nSpectra = spectra_block.dim0();
    spectraInfo.hasSpectra = true;
  } catch (std::runtime_error &) {
    spectraInfo.hasSpectra = false;
  }

  // Read necessary arrays from the file
  // Detector list contains a list of all of the detector numbers. If it not
  // present then we can't update the spectra
  // map
  try {
    NXInt detlist_group = detgroup.openNXInt("detector_list");
    detlist_group.load();
    spectraInfo.detectorList = detlist_group.sharedBuffer();
  } catch (std::runtime_error &) {
    logger.information() << "detector_list block not found. The workspace will "
                            "not contain any detector information.\n";
    return spectraInfo;
  }

  // Detector count contains the number of detectors associated with each
  // spectra
  NXInt det_count = detgroup.openNXInt("detector_count");
  det_count.load();
  spectraInfo.detectorCount = det_count.sharedBuffer();
  // Detector index - contains the index of the detector in the workspace
  NXInt det_index = detgroup.openNXInt("detector_index");
  det_index.load();
  spectraInfo.nSpectra = det_index.dim0();
  spectraInfo.detectorIndex = det_index.sharedBuffer();

  return spectraInfo;
}

/**
 * Is this file from a well-formed multiperiod group workspace.
 * @param nWorkspaceEntries : Number of entries in the group workspace
 * @param sampleWS : Sample workspace to inspect the logs of
 * @param log : Information logger object
 * @return True only if multiperiod.
 */
bool isMultiPeriodFile(int nWorkspaceEntries, const Workspace_sptr &sampleWS,
                       Logger &log) {
  bool isMultiPeriod = false;
  if (ExperimentInfo_sptr expInfo =
          boost::dynamic_pointer_cast<ExperimentInfo>(sampleWS)) {
    const std::string nPeriodsLogEntryName = "nperiods";
    const Run &run = expInfo->run();
    if (run.hasProperty(nPeriodsLogEntryName)) {
      const auto nPeriods =
          run.getPropertyValueAsType<int>(nPeriodsLogEntryName);
      if (nPeriods == nWorkspaceEntries) {
        isMultiPeriod = true;
        log.information("Loading as MultiPeriod group workspace.");
      }
    }
  }
  return isMultiPeriod;
}

} // namespace

/// Default constructor
LoadNexusProcessed::LoadNexusProcessed()
    : m_shared_bins(false), m_xbins(0), m_axis1vals(), m_list(false),
      m_interval(false), m_spec_min(0), m_spec_max(Mantid::EMPTY_INT()),
      m_spec_list(), m_filtered_spec_idxs(), m_nexusFile() {}

/// Destructor defined here so that NeXus::File can be forward declared
/// in header
LoadNexusProcessed::~LoadNexusProcessed() {}

/**
 * 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 LoadNexusProcessed::confidence(Kernel::NexusDescriptor &descriptor) const {
  if (descriptor.pathExists("/mantid_workspace_1"))
    return 80;
  else
    return 0;
}

void LoadNexusProcessed::readSpectraToDetectorMapping(
    NXEntry &mtd_entry, Mantid::API::MatrixWorkspace &ws) {
  readInstrumentGroup(mtd_entry, ws);
}

/** Initialisation method.
 *
 */
void LoadNexusProcessed::init() {
  // Declare required input parameters for algorithm
  const std::vector<std::string> exts{".nxs", ".nx5", ".xml"};
  declareProperty(
      std::make_unique<FileProperty>("Filename", "", FileProperty::Load, exts),
      "The name of the Nexus file to read, as a full or relative path.");
  declareProperty(std::make_unique<WorkspaceProperty<Workspace>>(
                      "OutputWorkspace", "", Direction::Output),
                  "The name of the workspace to be created as the output of "
                  "the algorithm.  A workspace of this name will be created "
                  "and stored in the Analysis Data Service. For multiperiod "
                  "files, one workspace may be generated for each period. "
                  "Currently only one workspace can be saved at a time so "
                  "multiperiod Mantid files are not generated.");

  // optional
  auto mustBePositive = boost::make_shared<BoundedValidator<int>>();
  mustBePositive->setLower(0);

  // Use a static cast as MSVC sometimes gets confused and casts as int64
  declareProperty("SpectrumMin", static_cast<int>(1), mustBePositive,
                  "Number of first spectrum to read.");
  declareProperty("SpectrumMax", static_cast<int>(Mantid::EMPTY_INT()),
                  mustBePositive, "Number of last spectrum to read.");
  declareProperty(std::make_unique<ArrayProperty<int>>("SpectrumList"),
                  "List of spectrum numbers to read.");
  declareProperty("EntryNumber", static_cast<int>(0), mustBePositive,
                  "0 indicates that every entry is loaded, into a separate "
                  "workspace within a group. "
                  "A positive number identifies one entry to be loaded, into "
                  "one workspace");
  declareProperty("LoadHistory", true,
                  "If true, the workspace history will be loaded");
  declareProperty(
      std::make_unique<PropertyWithValue<bool>>("FastMultiPeriod", true,
                                                Direction::Input),
      "For multiperiod workspaces. Copy instrument, parameter and x-data "
      "rather than loading it directly for each workspace. Y, E and log "
      "information is always loaded.");
}

/**
 * Loading specifically for mulitperiod group workspaces
 * @param root : NXRoot ref
 * @param entryName : Entry name to load.
 * @param tempMatrixWorkspace : Template workspace to base the next multiperiod
 * entry off.
 * @param nWorkspaceEntries : N entries in the file
 * @param p : index + 1 being processed.
 * @return Next multiperiod group workspace
 */
Workspace_sptr LoadNexusProcessed::doAccelleratedMultiPeriodLoading(
    NXRoot &root, const std::string &entryName,
    MatrixWorkspace_sptr &tempMatrixWorkspace, const size_t nWorkspaceEntries,
    const size_t p) {

  MatrixWorkspace_sptr periodWorkspace =
      WorkspaceFactory::Instance().create(tempMatrixWorkspace);

  const size_t nHistograms = periodWorkspace->getNumberHistograms();
  for (size_t i = 0; i < nHistograms; ++i) {
    periodWorkspace->setSharedX(i, tempMatrixWorkspace->sharedX(i));
  }

  // We avoid using `openEntry` or similar here because they're just wrappers
  // around `open`. `open` is slow for large multiperiod datasets, because it
  // does a search upon the entire HDF5 tree. `openLocal` is *much* quicker, as
  // it only searches the current group. It does, however, require that the
  // parent group is currently open.
  // Words of Warning: While the openLocal construct is an optimization,
  // it is very dangerous. Forgetting to close an entry of an NXEntry in a
  // completely unrelated part of the code can result in us opening the
  // wrong NXEntry here!
  NXEntry mtdEntry(root, entryName);
  mtdEntry.openLocal();

  NXData wsEntry(mtdEntry, "workspace");
  if (!wsEntry.openLocal()) {
    std::stringstream buffer;
    buffer << "Group entry " << p - 1
           << " is not a workspace 2D. Retry with "
              "FastMultiPeriod option set off.\n";
    throw std::runtime_error(buffer.str());
  }

  if (wsEntry.isValid("frac_area")) {
    std::stringstream buffer;
    buffer << "Group entry " << p - 1
           << " has fractional area present. Try "
              "reloading with FastMultiPeriod set "
              "off.\n";
    throw std::runtime_error(buffer.str());
  }

  NXDataSetTyped<double> data(wsEntry, "values");
  data.openLocal();
  NXDataSetTyped<double> errors(wsEntry, "errors");
  errors.openLocal();

  const int nChannels = data.dim1();

  int blockSize = 8; // Read block size. Set to 8 for efficiency. i.e. read
                     // 8 histograms at a time.
  int nFullBlocks =
      static_cast<int>(nHistograms) /
      blockSize; // Truncated number of full blocks to read. Remainder removed
  const int readOptimumStop = (nFullBlocks * blockSize);
  const int readStop = m_spec_max - 1;
  const int finalBlockSize = readStop - readOptimumStop;

  int wsIndex = 0;
  int histIndex = m_spec_min - 1;

  for (; histIndex < readStop;) {
    if (histIndex >= readOptimumStop) {
      blockSize = finalBlockSize;
    }

    data.load(blockSize, histIndex);
    errors.load(blockSize, histIndex);

    double *dataStart = data();
    double *dataEnd = dataStart + nChannels;

    double *errorStart = errors();
    double *errorEnd = errorStart + nChannels;

    int final(histIndex + blockSize);
    while (histIndex < final) {
      auto &Y = periodWorkspace->mutableY(wsIndex);
      Y.assign(dataStart, dataEnd);
      dataStart += nChannels;
      dataEnd += nChannels;
      auto &E = periodWorkspace->mutableE(wsIndex);
      E.assign(errorStart, errorEnd);
      errorStart += nChannels;
      errorEnd += nChannels;

      ++wsIndex;
      ++histIndex;
    }
  }

  // We always start one layer too deep
  // go from /workspace_{n}/{something} -> /workspace_{n}
  m_nexusFile->closeGroup();

  // Now move to the correct period group
  // /workspace_{n} -> /workspace_{n+1}
  m_nexusFile->closeGroup();
  m_nexusFile->openGroup(entryName, "NXentry");

  try {
    // This loads logs, sample, and instrument.
    periodWorkspace->loadSampleAndLogInfoNexus(m_nexusFile.get());
  } catch (std::exception &e) {
    g_log.information("Error loading Instrument section of nxs file");
    g_log.information(e.what());
  }

  // We make sure to close the current entries. Failing to do this can cause
  // strange off-by-one errors when loading the spectra.
  wsEntry.close();
  mtdEntry.close();

  const double fractionComplete = double(p - 1) / double(nWorkspaceEntries);
  progress(fractionComplete, "Loading multiperiod entry");
  return periodWorkspace;
}

//-------------------------------------------------------------------------------------------------
/** Executes the algorithm. Reading in the file and creating and populating
 *  the output workspace
 *
 *  @throw runtime_error Thrown if algorithm cannot execute
 */
void LoadNexusProcessed::exec() {

  API::Workspace_sptr tempWS;
  int nWorkspaceEntries = 0;
  // Start scoped block
  {
    progress(0, "Opening file...");

    // Throws an approriate exception if there is a problem with file access
    const std::string filename = getPropertyValue("Filename");
    NXRoot root(filename);

    // "Open" the same file but with the C++ interface
    m_nexusFile = std::make_unique<::NeXus::File>(root.m_fileID);

    // Find out how many first level entries there are
    // Cast down to int as another property later on is an int
    nWorkspaceEntries = static_cast<int>((root.groups().size()));

    // Check for an entry number property
    int entrynumber = getProperty("EntryNumber");
    Property const *const entryNumberProperty =
        this->getProperty("EntryNumber");
    bool bDefaultEntryNumber = entryNumberProperty->isDefault();

    if (!bDefaultEntryNumber && entrynumber > nWorkspaceEntries) {
      g_log.error() << "Invalid entry number specified. File only contains "
                    << nWorkspaceEntries << " entries.\n";
      throw std::invalid_argument("Invalid entry number specified.");
    }

    const std::string basename = "mantid_workspace_";

    std::ostringstream os;
    if (bDefaultEntryNumber) {
      // Set the entry number to 1 if not provided.
      entrynumber = 1;
    }
    os << basename << entrynumber;
    const std::string targetEntryName = os.str();

    // Take the first real workspace obtainable. We need it even if loading
    // groups.
    tempWS = loadEntry(root, targetEntryName, 0, 1);

    if (nWorkspaceEntries == 1 || !bDefaultEntryNumber) {
      // We have what we need.
      setProperty("OutputWorkspace", tempWS);
    } else {
      // We already know that this is a group workspace. Is it a true
      // multiperiod workspace.
      const bool bFastMultiPeriod = this->getProperty("FastMultiPeriod");
      const bool bIsMultiPeriod =
          isMultiPeriodFile(nWorkspaceEntries, tempWS, g_log);
      Property *specListProp = this->getProperty("SpectrumList");
      m_list = !specListProp->isDefault();

      // Load all first level entries
      auto wksp_group = boost::make_shared<WorkspaceGroup>();
      // This forms the name of the group
      std::string base_name = getPropertyValue("OutputWorkspace");
      // First member of group should be the group itself, for some reason!

      // load names of each of the workspaces. Note that if we have duplicate
      // names then we don't select them
      auto names =
          extractWorkspaceNames(root, static_cast<size_t>(nWorkspaceEntries));

      // remove existing workspace and replace with the one being loaded
      bool wsExists = AnalysisDataService::Instance().doesExist(base_name);
      if (wsExists) {
        Algorithm_sptr alg =
            AlgorithmManager::Instance().createUnmanaged("DeleteWorkspace");
        alg->initialize();
        alg->setChild(true);
        alg->setProperty("Workspace", base_name);
        alg->execute();
      }

      base_name += "_";
      const std::string prop_name = "OutputWorkspace_";

      MatrixWorkspace_sptr tempMatrixWorkspace =
          boost::dynamic_pointer_cast<Workspace2D>(tempWS);
      bool bAccelleratedMultiPeriodLoading = false;
      if (tempMatrixWorkspace) {
        // We only accelerate for simple scenarios for now. Spectrum lists are
        // too complicated to bother with.
        bAccelleratedMultiPeriodLoading =
            bIsMultiPeriod && bFastMultiPeriod && !m_list;
        // Strip out any loaded logs. That way we don't pay for copying that
        // information around.
        tempMatrixWorkspace->mutableRun().clearLogs();
      }

      if (bAccelleratedMultiPeriodLoading) {
        g_log.information("Accelerated multiperiod loading");
      } else {
        g_log.information("Individual group loading");
      }

      for (int p = 1; p <= nWorkspaceEntries; ++p) {
        const auto indexStr = std::to_string(p);

        // decide what the workspace should be called
        std::string wsName = buildWorkspaceName(names[p], base_name, p);

        Workspace_sptr local_workspace;

        /*
        For multiperiod workspaces we can accelerate the loading by making
        resonable assumptions about the differences between the workspaces
        Only Y, E and log data entries should vary. Therefore we can clone our
        temp workspace, and overwrite those things we are interested in.
        */
        if (bAccelleratedMultiPeriodLoading) {
          local_workspace = doAccelleratedMultiPeriodLoading(
              root, basename + indexStr, tempMatrixWorkspace, nWorkspaceEntries,
              p);
        } else // Fall-back for generic loading
        {
          const auto nWorkspaceEntries_d =
              static_cast<double>(nWorkspaceEntries);
          local_workspace =
              loadEntry(root, basename + indexStr,
                        static_cast<double>(p - 1) / nWorkspaceEntries_d,
                        1. / nWorkspaceEntries_d);
        }

        declareProperty(std::make_unique<WorkspaceProperty<API::Workspace>>(
            prop_name + indexStr, wsName, Direction::Output));

        wksp_group->addWorkspace(local_workspace);
        setProperty(prop_name + indexStr, local_workspace);
      }

      // The group is the root property value
      setProperty("OutputWorkspace",
                  boost::static_pointer_cast<Workspace>(wksp_group));
    }

    root.close();
  } // All file resources should be scoped to here. All previous file handles
  // must be cleared to release locks
  loadNexusGeometry(*tempWS, nWorkspaceEntries, g_log,
                    std::string(getProperty("Filename")));

  m_axis1vals.clear();
} // namespace DataHandling

/**
 * Decides what to call a child of a group workspace.
 *
 * This function builds the workspace name based on either a workspace name
 * which was stored in the file or the base name.
 *
 * @param name :: The name loaded from the file (possibly the empty string if
 *none was loaded)
 * @param baseName :: The name group workspace
 * @param wsIndex :: The current index of this workspace
 *
 * @return The name of the workspace
 */
std::string LoadNexusProcessed::buildWorkspaceName(const std::string &name,
                                                   const std::string &baseName,
                                                   size_t wsIndex) {
  std::string wsName;
  std::string index = std::to_string(wsIndex);

  if (!name.empty()) {
    wsName = name;
  } else {
    // we have a common stem so rename accordingly
    boost::smatch results;
    const boost::regex exp(".*_(\\d+$)");
    // if we have a common name stem then name is <OutputWorkspaceName>_n
    if (boost::regex_search(name, results, exp)) {
      wsName = baseName + std::string(results[1].first, results[1].second);
    } else {
      // if the name property wasn't defined just use <OutputWorkspaceName>_n
      wsName = baseName + index;
    }
  }

  correctForWorkspaceNameClash(wsName);

  return wsName;
}

/**
 * Append an index to the name if it already exists in the AnalysisDataService
 *
 * @param wsName :: Name to call the workspace
 */
void LoadNexusProcessed::correctForWorkspaceNameClash(std::string &wsName) {
  bool noClash(false);

  for (int i = 0; !noClash; ++i) {
    std::string wsIndex; // dont use an index if there is no other
                         // workspace
    if (i > 0) {
      wsIndex = "_" + std::to_string(i);
    }

    bool wsExists = AnalysisDataService::Instance().doesExist(wsName + wsIndex);
    if (!wsExists) {
      wsName += wsIndex;
      noClash = true;
    }
  }
}

/**
 * Extract the workspace names from the file (if any are stored)
 *
 * @param root :: the root for the NeXus document
 * @param nWorkspaceEntries :: the number of workspace entries
 */
std::vector<std::string>
LoadNexusProcessed::extractWorkspaceNames(NXRoot &root,
                                          size_t nWorkspaceEntries) {
  std::vector<std::string> names(nWorkspaceEntries + 1);
  for (size_t p = 1; p <= nWorkspaceEntries; ++p) {
    auto period = std::to_string(p);
    names[p] = loadWorkspaceName(root, "mantid_workspace_" + period);
  }

  // Check that there are no duplicates in the workspace name
  // This can cause severe problems
  auto it = std::unique(names.begin(), names.end());
  if (it != names.end()) {
    auto size = names.size();
    names.clear();
    names.resize(size);
  }
  return names;
}

/**
 * Load the workspace name, if the attribute exists
 *
 * @param root :: Root of NeXus file
 * @param entry_name :: Entry in NeXus file to look at
 * @return The workspace name. If none found an empty string is returned.
 */
std::string
LoadNexusProcessed::loadWorkspaceName(NXRoot &root,
                                      const std::string &entry_name) {
  NXEntry mtd_entry = root.openEntry(entry_name);
  std::string workspaceName = std::string();
  try {
    workspaceName = mtd_entry.getString("workspace_name");
  } catch (std::runtime_error &) {
  }
  mtd_entry.close();
  return workspaceName;
}

//-------------------------------------------------------------------------------------------------
/**
 * Load an event_workspace field
 *
 * @param wksp_cls  Nexus data for "event_workspace"
 * @param xbins bins on the "X" axis
 * @param progressStart algorithm progress (from 0)
 * @param progressRange  progress made after loading an entry
 *
 * @return event_workspace object with data
 */
API::MatrixWorkspace_sptr
LoadNexusProcessed::loadEventEntry(NXData &wksp_cls, NXDouble &xbins,
                                   const double &progressStart,
                                   const double &progressRange) {
  NXDataSetTyped<int64_t> indices_data =
      wksp_cls.openNXDataSet<int64_t>("indices");
  indices_data.load();
  size_t numspec = indices_data.dim0() - 1;

  // process optional spectrum parameters, if set
  checkOptionalProperties(numspec);
  // Actual number of spectra in output workspace (if only a user-specified
  // range and/or list was going to be loaded)
  numspec = calculateWorkspaceSize(numspec, true);

  int num_xbins = xbins.dim0();
  if (xbins.rank() == 2) {
    num_xbins = xbins.dim1();
  }
  if (num_xbins < 2)
    num_xbins = 2;
  EventWorkspace_sptr ws = boost::dynamic_pointer_cast<EventWorkspace>(
      WorkspaceFactory::Instance().create("EventWorkspace", numspec, num_xbins,
                                          num_xbins - 1));

  // Set the YUnit label
  ws->setYUnit(indices_data.attributes("units"));
  std::string unitLabel = indices_data.attributes("unit_label");
  if (unitLabel.empty())
    unitLabel = indices_data.attributes("units");
  ws->setYUnitLabel(unitLabel);

  // Handle optional fields.
  // TODO: Handle inconsistent sizes
  boost::shared_array<int64_t> pulsetimes;
  if (wksp_cls.isValid("pulsetime")) {
    NXDataSetTyped<int64_t> pulsetime =
        wksp_cls.openNXDataSet<int64_t>("pulsetime");
    pulsetime.load();
    pulsetimes = pulsetime.sharedBuffer();
  }

  boost::shared_array<double> tofs;
  if (wksp_cls.isValid("tof")) {
    NXDouble tof = wksp_cls.openNXDouble("tof");
    tof.load();
    tofs = tof.sharedBuffer();
  }

  boost::shared_array<float> error_squareds;
  if (wksp_cls.isValid("error_squared")) {
    NXFloat error_squared = wksp_cls.openNXFloat("error_squared");
    error_squared.load();
    error_squareds = error_squared.sharedBuffer();
  }

  boost::shared_array<float> weights;
  if (wksp_cls.isValid("weight")) {
    NXFloat weight = wksp_cls.openNXFloat("weight");
    weight.load();
    weights = weight.sharedBuffer();
  }

  // What type of event lists?
  EventType type = TOF;
  if (tofs && pulsetimes && weights && error_squareds)
    type = WEIGHTED;
  else if ((tofs && weights && error_squareds))
    type = WEIGHTED_NOTIME;
  else if (pulsetimes && tofs)
    type = TOF;
  else
    throw std::runtime_error("Could not figure out the type of event list!");

  // indices of events
  boost::shared_array<int64_t> indices = indices_data.sharedBuffer();
  // Create all the event lists
  auto max = static_cast<int64_t>(m_filtered_spec_idxs.size());
  Progress progress(this, progressStart, progressStart + progressRange, max);
  PARALLEL_FOR_NO_WSP_CHECK()
  for (int64_t j = 0; j < max; ++j) {
    PARALLEL_START_INTERUPT_REGION
    size_t wi = m_filtered_spec_idxs[j] - 1;
    int64_t index_start = indices[wi];
    int64_t index_end = indices[wi + 1];
    if (index_end >= index_start) {
      EventList &el = ws->getSpectrum(j);
      el.switchTo(type);

      // Allocate all the required memory
      el.reserve(index_end - index_start);
      el.clearDetectorIDs();

      for (int64_t i = index_start; i < index_end; i++)
        switch (type) {
        case TOF:
          el.addEventQuickly(TofEvent(tofs[i], DateAndTime(pulsetimes[i])));
          break;
        case WEIGHTED:
          el.addEventQuickly(WeightedEvent(tofs[i], DateAndTime(pulsetimes[i]),
                                           weights[i], error_squareds[i]));
          break;
        case WEIGHTED_NOTIME:
          el.addEventQuickly(
              WeightedEventNoTime(tofs[i], weights[i], error_squareds[i]));
          break;
        }

      // Set the X axis
      if (this->m_shared_bins)
        el.setHistogram(this->m_xbins);
      else {
        MantidVec x(xbins.dim1());

        for (int i = 0; i < xbins.dim1(); i++)
          x[i] = xbins(static_cast<int>(wi), i);
        // Workspace and el was just created, so we can just set a new histogram
        // We can move x as it is not longer used after this point
        el.setHistogram(HistogramData::BinEdges(std::move(x)));
      }
    }
    progress.report();
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION

  return ws;
}

//-------------------------------------------------------------------------------------------------
/**
 * Load a numeric column to the TableWorkspace.
 * @param tableData   :: Table data to load from
 * @param dataSetName :: Name of the data set to use to get column data
 * @param tableWs     :: Workspace to add column to
 * @param columnType  :: Name of the column type to create
 */
template <typename ColumnType, typename NexusType>
void LoadNexusProcessed::loadNumericColumn(
    const Mantid::NeXus::NXData &tableData, const std::string &dataSetName,
    const API::ITableWorkspace_sptr &tableWs, const std::string &columnType) {
  NXDataSetTyped<NexusType> data =
      tableData.openNXDataSet<NexusType>(dataSetName);
  std::string columnTitle = data.attributes("name");
  if (!columnTitle.empty()) {
    data.load();
    auto length = static_cast<size_t>(data.dim0());
    auto rowCount = tableWs->rowCount();
    // check that the row count is OK
    if (rowCount == 0) {
      tableWs->setRowCount(length);
    } else if (rowCount != length) {
      throw std::runtime_error("Columns have different sizes.");
    }
    // copy the data
    auto column = tableWs->addColumn(columnType, columnTitle);
    for (size_t i = 0; i < length; i++) {
      column->cell<ColumnType>(i) = static_cast<ColumnType>(*(data() + i));
    }
  }
}

//-------------------------------------------------------------------------------------------------
/**
 * Load a table
 */
API::Workspace_sptr LoadNexusProcessed::loadTableEntry(NXEntry &entry) {
  API::ITableWorkspace_sptr workspace;
  workspace =
      Mantid::API::WorkspaceFactory::Instance().createTable("TableWorkspace");

  NXData nx_tw = entry.openNXData("table_workspace");

  int columnNumber = 1;
  do {
    std::string dataSetName = "column_" + std::to_string(columnNumber);

    NXInfo info = nx_tw.getDataSetInfo(dataSetName);
    if (info.stat == NX_ERROR) {
      // Assume we done last column of table
      break;
    }

    if (info.rank == 1) {
      if (info.type == NX_FLOAT64) {
        loadNumericColumn<double, double>(nx_tw, dataSetName, workspace,
                                          "double");
      } else if (info.type == NX_INT32) {
        loadNumericColumn<int, int32_t>(nx_tw, dataSetName, workspace, "int");
      } else if (info.type == NX_UINT32) {
        loadNumericColumn<uint32_t, uint32_t>(nx_tw, dataSetName, workspace,
                                              "uint");
      } else if (info.type == NX_INT64) {
        loadNumericColumn<int64_t, int64_t>(nx_tw, dataSetName, workspace,
                                            "long64");
      } else if (info.type == NX_UINT64) {
        loadNumericColumn<size_t, uint64_t>(nx_tw, dataSetName, workspace,
                                            "size_t");
      } else if (info.type == NX_FLOAT32) {
        loadNumericColumn<float, float>(nx_tw, dataSetName, workspace, "float");
      } else if (info.type == NX_UINT8) {
        loadNumericColumn<bool, bool>(nx_tw, dataSetName, workspace, "bool");
      } else {
        throw std::logic_error("Column with Nexus data type " +
                               std::to_string(info.type) +
                               " cannot be loaded.");
      }
    } else if (info.rank == 2) {
      if (info.type == NX_CHAR) {
        NXChar data = nx_tw.openNXChar(dataSetName);
        std::string columnTitle = data.attributes("name");
        if (!columnTitle.empty()) {
          workspace->addColumn("str", columnTitle);
          int nRows = info.dims[0];
          workspace->setRowCount(nRows);

          const int maxStr = info.dims[1];
          data.load();
          for (int iR = 0; iR < nRows; ++iR) {
            auto &cellContents =
                workspace->cell<std::string>(iR, columnNumber - 1);
            auto startPoint = data() + maxStr * iR;
            cellContents.assign(startPoint, startPoint + maxStr);
            boost::trim_right(cellContents);
          }
        }
      } else if (info.type == NX_INT32) {
        loadVectorColumn<int>(nx_tw, dataSetName, workspace, "vector_int");
      } else if (info.type == NX_FLOAT64) {
        auto data = nx_tw.openNXDouble(dataSetName);
        if (data.attributes("interpret_as") == "V3D") {
          loadV3DColumn(data, workspace);
        } else {
          loadVectorColumn<double>(nx_tw, dataSetName, workspace,
                                   "vector_double");
        }
      }
    }

    columnNumber++;

  } while (true);

  return boost::static_pointer_cast<API::Workspace>(workspace);
}

/**
 * Loads a vector column to the TableWorkspace.
 * @param tableData   :: Table data to load from
 * @param dataSetName :: Name of the data set to use to get column data
 * @param tableWs     :: Workspace to add column to
 * @param columnType  :: Name of the column type to create
 */
template <typename Type>
void LoadNexusProcessed::loadVectorColumn(const NXData &tableData,
                                          const std::string &dataSetName,
                                          const ITableWorkspace_sptr &tableWs,
                                          const std::string &columnType) {
  NXDataSetTyped<Type> data = tableData.openNXDataSet<Type>(dataSetName);
  std::string columnTitle = data.attributes("name");
  if (!columnTitle.empty()) {
    tableWs->addColumn(columnType, columnTitle);

    NXInfo info = tableData.getDataSetInfo(dataSetName);
    const size_t rowCount = info.dims[0];
    const size_t blockSize = info.dims[1];

    // This might've been done already, but doing it twice should't do any harm
    tableWs->setRowCount(rowCount);

    data.load();

    for (size_t i = 0; i < rowCount; ++i) {
      auto &cell =
          tableWs->cell<std::vector<Type>>(i, tableWs->columnCount() - 1);

      Type *from = data() + blockSize * i;

      cell.assign(from, from + blockSize);

      std::ostringstream rowSizeAttrName;
      rowSizeAttrName << "row_size_" << i;

      // This is ugly, but I can only get attribute as a string using the API
      std::istringstream rowSizeStr(data.attributes(rowSizeAttrName.str()));

      int rowSize;
      rowSizeStr >> rowSize;

      cell.resize(rowSize);
    }
  }
}

/**
 * Loads a V3D column to the TableWorkspace.
 * @param data   :: Table data to load from
 * @param tableWs     :: Workspace to add column to
 */
void LoadNexusProcessed::loadV3DColumn(
    Mantid::NeXus::NXDouble &data, const API::ITableWorkspace_sptr &tableWs) {
  std::string columnTitle = data.attributes("name");
  if (!columnTitle.empty()) {
    ColumnVector<V3D> col = tableWs->addColumn("V3D", columnTitle);

    const int rowCount = data.dim0();

    // This might've been done already, but doing it twice should't do any harm
    tableWs->setRowCount(rowCount);

    data.load();

    for (int i = 0; i < rowCount; ++i) {
      auto &cell = col[i];
      cell(data(i, 0), data(i, 1), data(i, 2));
    }
  }
}

//-------------------------------------------------------------------------------------------------
/**
 * Load peaks
 */
API::Workspace_sptr LoadNexusProcessed::loadPeaksEntry(NXEntry &entry) {
  // API::IPeaksWorkspace_sptr workspace;
  API::ITableWorkspace_sptr tWorkspace;
  // PeaksWorkspace_sptr workspace;
  tWorkspace =
      Mantid::API::WorkspaceFactory::Instance().createTable("PeaksWorkspace");

  PeaksWorkspace_sptr peakWS =
      boost::dynamic_pointer_cast<PeaksWorkspace>(tWorkspace);

  NXData nx_tw = entry.openNXData("peaks_workspace");

  int columnNumber = 1;
  int numberPeaks = 0;
  std::vector<std::string> columnNames;
  do {
    std::string str = "column_" + std::to_string(columnNumber);

    NXInfo info = nx_tw.getDataSetInfo(str);
    if (info.stat == NX_ERROR) {
      // Assume we done last column of table
      break;