LoadNexusProcessed.cpp

Go to the documentation of this file.

// 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 "MantidDataHandling/ISISRunLogs.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/LeanElasticPeaksWorkspace.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/FacilityInfo.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/StringTokenizer.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidNexus/NexusClasses.h"
#include "MantidNexus/NexusFileIO.h"
 
#include <boost/algorithm/string/trim.hpp>
#include <boost/regex.hpp>
#include <nexus/NeXusException.hpp>
 
#include <map>
#include <memory>
#include <string>
#include <vector>
 
namespace Mantid::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 = std::vector<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 spectraNumbers;
  // Index of the detector in the workspace.
  IntArray detectorIndex;
  // Number of detectors associated with each spectra
  IntArray detectorCount;
  // Detector list contains a list of all of the detector numbers
  IntArray detectorList;
};
 
// Helper typdef.
using SpectraInfo_optional = boost::optional<SpectraInfo>;
 
SpectraInfo extractMappingInfo(const 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.vecBuffer();
    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.vecBuffer();
  } 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.vecBuffer();
  // 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.vecBuffer();
 
  return spectraInfo;
}
 
bool isMultiPeriodFile(int nWorkspaceEntries, const Workspace_sptr &sampleWS, Logger &log) {
  bool isMultiPeriod = false;
  if (ExperimentInfo_sptr expInfo = std::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
 
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() {}
 
LoadNexusProcessed::~LoadNexusProcessed() = default;
 
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);
}
 
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 = std::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.");
}
 
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;
}
 
//-------------------------------------------------------------------------------------------------
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.
      applyLogFiltering(tempWS);
      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 = std::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 = std::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);
        }
 
        applyLogFiltering(local_workspace);
        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", std::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
 
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;
}
 
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;
    }
  }
}
 
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;
}
 
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;
}
 
//-------------------------------------------------------------------------------------------------
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 = std::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
  std::vector<int64_t> pulsetimes;
  if (wksp_cls.isValid("pulsetime")) {
    NXDataSetTyped<int64_t> pulsetime = wksp_cls.openNXDataSet<int64_t>("pulsetime");
    pulsetime.load();
    pulsetimes = pulsetime.vecBuffer();
  }
 
  std::vector<double> tofs;
  if (wksp_cls.isValid("tof")) {
    NXDouble tof = wksp_cls.openNXDouble("tof");
    tof.load();
    tofs = tof.vecBuffer();
  }
 
  std::vector<float> error_squareds;
  if (wksp_cls.isValid("error_squared")) {
    NXFloat error_squared = wksp_cls.openNXFloat("error_squared");
    error_squared.load();
    error_squareds = error_squared.vecBuffer();
  }
 
  std::vector<float> weights;
  if (wksp_cls.isValid("weight")) {
    NXFloat weight = wksp_cls.openNXFloat("weight");
    weight.load();
    weights = weight.vecBuffer();
  }
 
  // What type of event lists?
  EventType type = TOF;
  if (!tofs.empty() && !pulsetimes.empty() && !weights.empty() && !error_squareds.empty())
    type = WEIGHTED;
  else if ((!tofs.empty() && !weights.empty() && !error_squareds.empty()))
    type = WEIGHTED_NOTIME;
  else if (!pulsetimes.empty() && !tofs.empty())
    type = TOF;
  else
    throw std::runtime_error("Could not figure out the type of event list!");
 
  // indices of events
  std::vector<int64_t> indices = indices_data.vecBuffer();
  // 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_INTERRUPT_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_INTERRUPT_REGION
  }
  PARALLEL_CHECK_INTERRUPT_REGION
 
  return ws;
}
 
//-------------------------------------------------------------------------------------------------
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));
    }
  }
}
 
//-------------------------------------------------------------------------------------------------
API::Workspace_sptr LoadNexusProcessed::loadTableEntry(const 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 std::static_pointer_cast<API::Workspace>(workspace);
}
 
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);
    }
  }
}
 
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));
    }
  }
}
 
API::Workspace_sptr LoadNexusProcessed::loadLeanElasticPeaksEntry(NXEntry &entry) {
  g_log.notice("Load as LeanElasticPeaks");
 
  // API::IPeaksWorkspace_sptr workspace;
  API::ITableWorkspace_sptr tWorkspace;
  // PeaksWorkspace_sptr workspace;
  tWorkspace = Mantid::API::WorkspaceFactory::Instance().createTable("LeanElasticPeaksWorkspace");
 
  IPeaksWorkspace_sptr peakWS = std::dynamic_pointer_cast<LeanElasticPeaksWorkspace>(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;
    }
 
    // store column names
    columnNames.emplace_back(str);
 
    // determine number of peaks
    // here we assume that a peaks_table has always one column of doubles
 
    if (info.type == NX_FLOAT64) {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      std::string columnTitle = nxDouble.attributes("name");
      if (!columnTitle.empty() && numberPeaks == 0) {
        numberPeaks = nxDouble.dim0();
      }
    }
 
    columnNumber++;
 
  } while (true);
 
  // Get information from all but data group
  std::string parameterStr;
  // Hop to the right point /mantid_workspace_1
  try {
    m_nexusFile->openPath(entry.path()); // This is
  } catch (std::runtime_error &re) {
    throw std::runtime_error("Error while opening a path in a Peaks entry in a "
                             "Nexus processed file. "
                             "This path is wrong: " +
                             entry.path() + ". Lower level error description: " + re.what());
  }
  try {
    // This loads logs, sample, and instrument.
    peakWS->loadExperimentInfoNexus(getPropertyValue("Filename"), m_nexusFile.get(), parameterStr);
    // Populate the instrument parameters in this workspace
    peakWS->readParameterMap(parameterStr);
  } catch (std::exception &e) {
    g_log.information("Error loading Instrument section of nxs file");
    g_log.information(e.what());
  }
 
  // Coordinates - Older versions did not have the separate field but used a log
  // value
  const std::string peaksWSName = "peaks_workspace";
  try {
    m_nexusFile->openGroup(peaksWSName, "NXentry");
  } catch (std::runtime_error &re) {
    throw std::runtime_error("Error while opening a peaks workspace in a Nexus processed file. "
                             "Cannot open gropu " +
                             peaksWSName + ". Lower level error description: " + re.what());
  }
  try {
    uint32_t loadCoord(0);
    m_nexusFile->readData("coordinate_system", loadCoord);
    peakWS->setCoordinateSystem(static_cast<Kernel::SpecialCoordinateSystem>(loadCoord));
  } catch (::NeXus::Exception &) {
    // Check for a log value
    auto logs = peakWS->logs();
    if (logs->hasProperty("CoordinateSystem")) {
      auto *prop = dynamic_cast<PropertyWithValue<int> *>(logs->getProperty("CoordinateSystem"));
      if (prop) {
        int value((*prop)());
        peakWS->setCoordinateSystem(static_cast<Kernel::SpecialCoordinateSystem>(value));
      }
    }
  }
 
  std::string m_QConvention = "Inelastic";
  try {
    m_nexusFile->getAttr("QConvention", m_QConvention);
  } catch (std::exception &) {
  }
 
  // peaks_workspace
  m_nexusFile->closeGroup();
 
  // Change convention of loaded file to that in Preferen
  double qSign = 1.0;
  std::string convention = ConfigService::Instance().getString("Q.convention");
  if (convention != m_QConvention)
    qSign = -1.0;
 
  for (int r = 0; r < numberPeaks; r++) {
    // Create individual LeanElasticPeak
    const auto goniometer = peakWS->run().getGoniometer();
    LeanElasticPeak peak;
    peak.setGoniometerMatrix(goniometer.getR());
    peak.setRunNumber(peakWS->getRunNumber());
    peakWS->addPeak(peak);
  }
 
  for (const auto &str : columnNames) {
    if (str == "column_1") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = qSign * nxDouble[r];
        peakWS->getPeak(r).setH(val);
      }
    } else if (str == "column_2") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = qSign * nxDouble[r];
        peakWS->getPeak(r).setK(val);
      }
    } else if (str == "column_3") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = qSign * nxDouble[r];
        peakWS->getPeak(r).setL(val);
      }
    } else if (str == "column_4") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = nxDouble[r];
        peakWS->getPeak(r).setIntensity(val);
      }
    } else if (str == "column_5") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = nxDouble[r];
        peakWS->getPeak(r).setSigmaIntensity(val);
      }
    } else if (str == "column_6") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = nxDouble[r];
        peakWS->getPeak(r).setBinCount(val);
      }
    } else if (str == "column_7") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = nxDouble[r];
        peakWS->getPeak(r).setWavelength(val);
      }
    } else if (str == "column_10") {
      NXInt nxInt = nx_tw.openNXInt(str);
      nxInt.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        int ival = nxInt[r];
        if (ival != -1)
          peakWS->getPeak(r).setRunNumber(ival);
      }
    } else if (str == "column_11") {
      NXInt nxInt = nx_tw.openNXInt(str);
      nxInt.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        int ival = nxInt[r];
        peakWS->getPeak(r).setPeakNumber(ival);
      }
    } else if (str == "column_12") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = nxDouble[r];
        peakWS->getPeak(r).setAbsorptionWeightedPathLength(val);
      }
    } else if (str == "column_13") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
      Kernel::Matrix<double> gm(3, 3, false);
      int k = 0;
      for (int r = 0; r < numberPeaks; r++) {
        for (int j = 0; j < 9; j++) {
          double val = nxDouble[k];
          k++;
          gm[j % 3][j / 3] = val;
        }
        peakWS->getPeak(r).setGoniometerMatrix(gm);
      }
    } else if (str == "column_15") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
      V3D qlab;
      for (int r = 0; r < numberPeaks; ++r) {
        qlab = V3D(nxDouble[r * 3], nxDouble[r * 3 + 1], nxDouble[r * 3 + 2]);
        peakWS->getPeak(r).setQLabFrame(qlab, 0.0);
      }
    }
 
    // After all columns read set IntHKL
    for (int r = 0; r < numberPeaks; r++) {
      V3D intHKL = V3D(peakWS->getPeak(r).getH(), peakWS->getPeak(r).getK(), peakWS->getPeak(r).getL());
      peakWS->getPeak(r).setIntHKL(intHKL);
    }
  }
 
  return std::static_pointer_cast<API::Workspace>(peakWS);
}
 
//-------------------------------------------------------------------------------------------------
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 = std::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;
    }
 
    // store column names
    columnNames.emplace_back(str);
 
    // determine number of peaks
    // here we assume that a peaks_table has always one column of doubles
 
    if (info.type == NX_FLOAT64) {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      std::string columnTitle = nxDouble.attributes("name");
      if (!columnTitle.empty() && numberPeaks == 0) {
        numberPeaks = nxDouble.dim0();
      }
    }
 
    columnNumber++;
 
  } while (true);
 
  // Get information from all but data group
  std::string parameterStr;
  // Hop to the right point /mantid_workspace_1
  try {
    m_nexusFile->openPath(entry.path()); // This is
  } catch (std::runtime_error &re) {
    throw std::runtime_error("Error while opening a path in a Peaks entry in a "
                             "Nexus processed file. "
                             "This path is wrong: " +
                             entry.path() + ". Lower level error description: " + re.what());
  }
  try {
    // This loads logs, sample, and instrument.
    peakWS->loadExperimentInfoNexus(getPropertyValue("Filename"), m_nexusFile.get(), parameterStr);
    // Populate the instrument parameters in this workspace
    peakWS->readParameterMap(parameterStr);
  } catch (std::exception &e) {
    g_log.information("Error loading Instrument section of nxs file");
    g_log.information(e.what());
  }
 
  // Coordinates - Older versions did not have the separate field but used a log
  // value
  const std::string peaksWSName = "peaks_workspace";
  try {
    m_nexusFile->openGroup(peaksWSName, "NXentry");
  } catch (std::runtime_error &re) {
    throw std::runtime_error("Error while opening a peaks workspace in a Nexus processed file. "
                             "Cannot open gropu " +
                             peaksWSName + ". Lower level error description: " + re.what());
  }
  try {
    uint32_t loadCoord(0);
    m_nexusFile->readData("coordinate_system", loadCoord);
    peakWS->setCoordinateSystem(static_cast<Kernel::SpecialCoordinateSystem>(loadCoord));
  } catch (::NeXus::Exception &) {
    // Check for a log value
    auto logs = peakWS->logs();
    if (logs->hasProperty("CoordinateSystem")) {
      auto *prop = dynamic_cast<PropertyWithValue<int> *>(logs->getProperty("CoordinateSystem"));
      if (prop) {
        int value((*prop)());
        peakWS->setCoordinateSystem(static_cast<Kernel::SpecialCoordinateSystem>(value));
      }
    }
  }
 
  std::string m_QConvention = "Inelastic";
  try {
    m_nexusFile->getAttr("QConvention", m_QConvention);
  } catch (std::exception &) {
  }
 
  // peaks_workspace
  m_nexusFile->closeGroup();
 
  // Change convention of loaded file to that in Preferen
  double qSign = 1.0;
  std::string convention = ConfigService::Instance().getString("Q.convention");
  if (convention != m_QConvention)
    qSign = -1.0;
 
  for (int r = 0; r < numberPeaks; r++) {
    // Warning! Do not use anything other than the default constructor here
    // It is currently important (10/05/17) that the DetID (set in the loop
    // below this one) is set before QLabFrame as this causes Peak to ray trace
    // to find the location of the detector, which significantly increases
    // loading times.
    const auto goniometer = peakWS->run().getGoniometer();
    Peak peak;
    peak.setInstrument(peakWS->getInstrument());
    peak.setGoniometerMatrix(goniometer.getR());
    peak.setRunNumber(peakWS->getRunNumber());
    peakWS->addPeak(std::move(peak));
  }
 
  for (const auto &str : columnNames) {
    if (str == "column_1") {
      NXInt nxInt = nx_tw.openNXInt(str);
      nxInt.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        int ival = nxInt[r];
        if (ival != -1)
          peakWS->getPeak(r).setDetectorID(ival);
      }
    } else if (str == "column_2") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = qSign * nxDouble[r];
        peakWS->getPeak(r).setH(val);
      }
    } else if (str == "column_3") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = qSign * nxDouble[r];
        peakWS->getPeak(r).setK(val);
      }
    } else if (str == "column_4") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = qSign * nxDouble[r];
        peakWS->getPeak(r).setL(val);
      }
    } else if (str == "column_5") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = nxDouble[r];
        peakWS->getPeak(r).setIntensity(val);
      }
    } else if (str == "column_6") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = nxDouble[r];
        peakWS->getPeak(r).setSigmaIntensity(val);
      }
    } else if (str == "column_7") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = nxDouble[r];
        peakWS->getPeak(r).setBinCount(val);
      }
    } else if (str == "column_10") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = nxDouble[r];
        peakWS->getPeak(r).setWavelength(val);
      }
    } else if (str == "column_14") {
      NXInt nxInt = nx_tw.openNXInt(str);
      nxInt.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        int ival = nxInt[r];
        if (ival != -1)
          peakWS->getPeak(r).setRunNumber(ival);
      }
    } else if (str == "column_17") {
      NXInt nxInt = nx_tw.openNXInt(str);
      nxInt.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        int ival = nxInt[r];
        peakWS->getPeak(r).setPeakNumber(ival);
      }
    } else if (str == "column_18") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
 
      for (int r = 0; r < numberPeaks; r++) {
        double val = nxDouble[r];
        peakWS->getPeak(r).setAbsorptionWeightedPathLength(val);
      }
    } else if (str == "column_15") {
      NXDouble nxDouble = nx_tw.openNXDouble(str);
      nxDouble.load();
      Kernel::Matrix<double> gm(3, 3, false);
      int k = 0;
      for (int r = 0; r < numberPeaks; r++) {
        for (int j = 0; j < 9; j++) {
          double val = nxDouble[k];
          k++;
          gm[j % 3][j / 3] = val;
        }
        peakWS->getPeak(r).setGoniometerMatrix(gm);
      }
    } else if (str == "column_16") {
      // Read shape information
      using namespace Mantid::DataObjects;
 
      PeakShapeFactory_sptr peakFactoryEllipsoid = std::make_shared<PeakShapeEllipsoidFactory>();
      PeakShapeFactory_sptr peakFactorySphere = std::make_shared<PeakShapeSphericalFactory>();
      PeakShapeFactory_sptr peakFactoryNone = std::make_shared<PeakNoShapeFactory>();
 
      peakFactoryEllipsoid->setSuccessor(peakFactorySphere);
      peakFactorySphere->setSuccessor(peakFactoryNone);
 
      NXInfo info = nx_tw.getDataSetInfo(str);
      NXChar data = nx_tw.openNXChar(str);
 
      const int maxShapeJSONLength = info.dims[1];
      data.load();
      for (int i = 0; i < numberPeaks; ++i) {
 
        // iR = peak row number
        auto startPoint = data() + (maxShapeJSONLength * i);
        std::string shapeJSON(startPoint, startPoint + maxShapeJSONLength);
        boost::trim_right(shapeJSON);
 
        // Make the shape
        Mantid::Geometry::PeakShape *peakShape = peakFactoryEllipsoid->create(shapeJSON);
 
        // Set the shape
        peakWS->getPeak(i).setPeakShape(peakShape);
      }
    }
    // After all columns read set IntHKL
    for (int r = 0; r < numberPeaks; r++) {
      V3D intHKL = V3D(peakWS->getPeak(r).getH(), peakWS->getPeak(r).getK(), peakWS->getPeak(r).getL());
      peakWS->getPeak(r).setIntHKL(intHKL);
    }
  }
 
  return std::static_pointer_cast<API::Workspace>(peakWS);
}
 
//-------------------------------------------------------------------------------------------------
API::MatrixWorkspace_sptr LoadNexusProcessed::loadNonEventEntry(NXData &wksp_cls, NXDouble &xbins,
                                                                const double &progressStart,
                                                                const double &progressRange, const NXEntry &mtd_entry,
                                                                const int xlength, std::string &workspaceType) {
  // Filter the list of spectra to process, applying min/max/list options
  NXDataSetTyped<double> data = wksp_cls.openDoubleData();
  int nchannels = data.dim1();
  size_t nspectra = data.dim0();
  // process optional spectrum parameters, if set
  checkOptionalProperties(nspectra);
  // Actual number of spectra in output workspace (if only a range was going to be loaded)
  size_t total_specs = calculateWorkspaceSize(nspectra);
 
  if (nchannels == 1 && nspectra == 1) {
    // if there is only one value of channels and nspectra, it may be a WorkspaceSingleValue
    // check for instrument
    bool hasInstrument = mtd_entry.containsGroup("instrument");
    if (hasInstrument) {
      std::string inst_name = mtd_entry.getString("instrument/name");
      boost::algorithm::trim(inst_name);
      if (inst_name == "")
        hasInstrument = false;
    } else {
      // data saved with SaveNexusESS will have the instrument in a directory named after it
      // check for special types of instrument: "basic_rect" and "unspecified_instrument":
      if (mtd_entry.containsGroup("basic_rect") || mtd_entry.containsGroup("unspecified_instrument")) {
        hasInstrument = true;
      } else {
        // check for other possible instruments
        for (auto facility : ConfigService::Instance().getFacilities()) {
          for (auto instrumentName : facility->instruments()) {
            if (instrumentName.name() != "" && mtd_entry.containsGroup(instrumentName.name())) {
              hasInstrument = true;
              break;
            }
          }
        }
      }
    }
    // check for metadata
    bool hasMetadata = mtd_entry.containsGroup("logs");
    if (hasMetadata) {
      // if there is more than one log (called "goniometer") then it's not a single-valued ws
      const auto nLogs = mtd_entry.openNXGroup("logs").groups().size();
      if (nLogs <= 1) { // only "goniometer" group is present, thus it's a single-valued ws
        hasMetadata = false;
      }
    }
    // a workspace with no instrument and no metadata, and only one entry is a single-valued ws
    if (!hasInstrument && !hasMetadata)
      workspaceType = "WorkspaceSingleValue";
  }
  bool hasFracArea = false;
  if (wksp_cls.isValid("frac_area")) {
    // frac_area entry is the signal for a RebinnedOutput workspace
    hasFracArea = true;
    workspaceType = "RebinnedOutput";
  }
 
  API::MatrixWorkspace_sptr local_workspace = std::dynamic_pointer_cast<API::MatrixWorkspace>(
      WorkspaceFactory::Instance().create(workspaceType, total_specs, xlength, nchannels));
  try {
    local_workspace->setTitle(mtd_entry.getString("title"));
  } catch (std::runtime_error &) {
    g_log.debug() << "No title was found in the input file, " << getPropertyValue("Filename") << '\n';
  }
 
  // Set the YUnit label
  local_workspace->setYUnit(data.attributes("units"));
  std::string unitLabel = data.attributes("unit_label");
  if (unitLabel.empty())
    unitLabel = data.attributes("units");
  local_workspace->setYUnitLabel(unitLabel);
 
  readBinMasking(wksp_cls, local_workspace);
  NXDataSetTyped<double> errors = wksp_cls.openNXDouble("errors");
  NXDataSetTyped<double> fracarea = errors;
  if (hasFracArea) {
    fracarea = wksp_cls.openNXDouble("frac_area");
 
    // Set the fractional area attributes, default values consistent with
    // previous assumptions: finalized = true, sqrdErrs = false
    auto rbWS = std::dynamic_pointer_cast<RebinnedOutput>(local_workspace);
    auto finalizedValue = fracarea.attributes("finalized");
    auto finalized = (finalizedValue.empty() ? true : finalizedValue == "1");
    rbWS->setFinalized(finalized);
    auto sqrdErrsValue = fracarea.attributes("sqrd_errors");
    auto sqrdErrs = (sqrdErrsValue.empty() ? false : sqrdErrsValue == "1");
    rbWS->setSqrdErrors(sqrdErrs);
  }
 
  // Check for x errors; as with fracArea we set it to xbins
  // although in this case it would never be used.
  auto hasXErrors = wksp_cls.isValid("xerrors");
  auto xErrors = hasXErrors ? wksp_cls.openNXDouble("xerrors") : errors;
  if (hasXErrors) {
    if (xErrors.dim1() == nchannels + 1)
      g_log.warning() << "Legacy X uncertainty found in input file, i.e., "
                         "delta-Q for each BIN EDGE. Uncertainties will be "
                         "re-interpreted as delta-Q of the BIN CENTRE and the "
                         "last value will be dropped.\n";
  }
 
  int blocksize = 8;
  // const int fullblocks = nspectra / blocksize;
  // size of the workspace
  // have to cast down to int as later functions require ints
  int fullblocks = static_cast<int>(total_specs) / blocksize;
  int read_stop = (fullblocks * blocksize);
  const double progressBegin = progressStart + 0.25 * progressRange;
  const double progressScaler = 0.75 * progressRange;
  int hist_index = 0;
  int wsIndex = 0;
  if (m_shared_bins) {
    // if spectrum min,max,list properties are set
    if (m_interval || m_list) {
      // if spectrum max,min properties are set read the data as a
      // block(multiple of 8) and
      // then read the remaining data as finalblock
      if (m_interval) {
        // specs at the min-max interval
        int interval_specs = m_spec_max - m_spec_min;
        fullblocks = (interval_specs) / blocksize;
        read_stop = (fullblocks * blocksize) + m_spec_min - 1;
 
        if (interval_specs < blocksize) {
          blocksize = static_cast<int>(total_specs);
          read_stop = m_spec_max - 1;
        }
        hist_index = m_spec_min - 1;
 
        for (; hist_index < read_stop;) {
          progress(progressBegin + progressScaler * static_cast<double>(hist_index) / static_cast<double>(read_stop),
                   "Reading workspace data...");
          loadBlock(data, errors, fracarea, hasFracArea, xErrors, hasXErrors, blocksize, nchannels, hist_index, wsIndex,
                    local_workspace);
        }
        size_t finalblock = m_spec_max - 1 - read_stop;
        if (finalblock > 0) {
          loadBlock(data, errors, fracarea, hasFracArea, xErrors, hasXErrors, static_cast<int>(finalblock), nchannels,
                    hist_index, wsIndex, local_workspace);
        }
      }
      // if spectrum list property is set read each spectrum separately by
      // setting blocksize=1
      if (m_list) {
        for (const auto itr : m_spec_list) {
          int specIndex = itr - 1;
          progress(progressBegin +
                       progressScaler * static_cast<double>(specIndex) / static_cast<double>(m_spec_list.size()),
                   "Reading workspace data...");
          loadBlock(data, errors, fracarea, hasFracArea, xErrors, hasXErrors, 1, nchannels, specIndex, wsIndex,
                    local_workspace);
        }
      }
    } else {
      for (; hist_index < read_stop;) {
        progress(progressBegin + progressScaler * static_cast<double>(hist_index) / static_cast<double>(read_stop),
                 "Reading workspace data...");
        loadBlock(data, errors, fracarea, hasFracArea, xErrors, hasXErrors, blocksize, nchannels, hist_index, wsIndex,
                  local_workspace);
      }
      size_t finalblock = total_specs - read_stop;
      if (finalblock > 0) {
        loadBlock(data, errors, fracarea, hasFracArea, xErrors, hasXErrors, static_cast<int>(finalblock), nchannels,
                  hist_index, wsIndex, local_workspace);
      }
    }
 
  } else {
    if (m_interval || m_list) {
      if (m_interval) {
        int interval_specs = m_spec_max - m_spec_min;
        fullblocks = (interval_specs) / blocksize;
        read_stop = (fullblocks * blocksize) + m_spec_min - 1;
 
        if (interval_specs < blocksize) {
          blocksize = interval_specs;
          read_stop = m_spec_max - 1;
        }
        hist_index = m_spec_min - 1;
 
        for (; hist_index < read_stop;) {
          progress(progressBegin + progressScaler * static_cast<double>(hist_index) / static_cast<double>(read_stop),
                   "Reading workspace data...");
          loadBlock(data, errors, fracarea, hasFracArea, xErrors, hasXErrors, xbins, blocksize, nchannels, hist_index,
                    wsIndex, local_workspace);
        }
        size_t finalblock = m_spec_max - 1 - read_stop;
        if (finalblock > 0) {
          loadBlock(data, errors, fracarea, hasFracArea, xErrors, hasXErrors, xbins, static_cast<int>(finalblock),
                    nchannels, hist_index, wsIndex, local_workspace);
        }
      }
      //
      if (m_list) {
        for (const auto itr : m_spec_list) {
          int specIndex = itr - 1;
          progress(progressBegin + progressScaler * static_cast<double>(specIndex) / static_cast<double>(read_stop),
                   "Reading workspace data...");
          loadBlock(data, errors, fracarea, hasFracArea, xErrors, hasXErrors, xbins, 1, nchannels, specIndex, wsIndex,
                    local_workspace);
        }
      }
    } else {
      for (; hist_index < read_stop;) {
        progress(progressBegin + progressScaler * static_cast<double>(hist_index) / static_cast<double>(read_stop),
                 "Reading workspace data...");
        loadBlock(data, errors, fracarea, hasFracArea, xErrors, hasXErrors, xbins, blocksize, nchannels, hist_index,
                  wsIndex, local_workspace);
      }
      size_t finalblock = total_specs - read_stop;
      if (finalblock > 0) {
        loadBlock(data, errors, fracarea, hasFracArea, xErrors, hasXErrors, xbins, static_cast<int>(finalblock),
                  nchannels, hist_index, wsIndex, local_workspace);
      }
    }
  }
  return local_workspace;
}
 
//-------------------------------------------------------------------------------------------------
API::Workspace_sptr LoadNexusProcessed::loadEntry(NXRoot &root, const std::string &entry_name,
                                                  const double &progressStart, const double &progressRange) {
  progress(progressStart, "Opening entry " + entry_name + "...");
 
  NXEntry mtd_entry = root.openEntry(entry_name);
 
  if (mtd_entry.containsGroup("table_workspace")) {
    return loadTableEntry(mtd_entry);
  }
 
  if (mtd_entry.containsGroup("peaks_workspace")) {
    try {
      // try standard PeakWorkspace first
      return loadPeaksEntry(mtd_entry);
    } catch (std::exception &) {
      return loadLeanElasticPeaksEntry(mtd_entry);
    }
  }
 
  // Determine workspace type and name of group containing workspace
  // characteristics
  bool isEvent = false;
  std::string workspaceType = "Workspace2D";
  std::string group_name = "workspace";
  if (mtd_entry.containsGroup("event_workspace")) {
    isEvent = true;
    group_name = "event_workspace";
  } else if (mtd_entry.containsGroup("offsets_workspace")) {
    workspaceType = "OffsetsWorkspace";
    group_name = "offsets_workspace";
  } else if (mtd_entry.containsGroup("mask_workspace")) {
    workspaceType = "MaskWorkspace";
    group_name = "mask_workspace";
  }
 
  // Get workspace characteristics
  NXData wksp_cls = mtd_entry.openNXData(group_name);
 
  // Axis information
  // "X" axis
 
  NXDouble xbins = wksp_cls.openNXDouble("axis1");
  xbins.load();
  std::string unit1 = xbins.attributes("units");
  // Non-uniform x bins get saved as a 2D 'axis1' dataset
  int xlength(-1);
  if (xbins.rank() == 2) {
    xlength = xbins.dim1();
    m_shared_bins = false;
  } else if (xbins.rank() == 1) {
    xlength = xbins.dim0();
    m_shared_bins = true;
    xbins.load();
    m_xbins = HistogramData::HistogramX(xbins(), xbins() + xlength);
  } else {
    throw std::runtime_error("Unknown axis1 dimension encountered.");
  }
 
  // MatrixWorkspace axis 1
  NXDouble axis2 = wksp_cls.openNXDouble("axis2");
  std::string unit2 = axis2.attributes("units");
 
  // --- Load workspace (as event_workspace or workspace2d) ---
  API::MatrixWorkspace_sptr local_workspace;
  if (isEvent) {
    local_workspace = loadEventEntry(wksp_cls, xbins, progressStart, progressRange);
  } else {
    local_workspace =
        loadNonEventEntry(wksp_cls, xbins, progressStart, progressRange, mtd_entry, xlength, workspaceType);
  }
  size_t nspectra = local_workspace->getNumberHistograms();
 
  // Units
  bool verticalHistogram(false);
  try {
    local_workspace->getAxis(0)->unit() = UnitFactory::Instance().create(unit1);
    if (unit1 == "Label") {
      auto label = std::dynamic_pointer_cast<Mantid::Kernel::Units::Label>(local_workspace->getAxis(0)->unit());
      auto ax = wksp_cls.openNXDouble("axis1");
      label->setLabel(ax.attributes("caption"), ax.attributes("label"));
    }
 
    // If this doesn't throw then it is a numeric access so grab the data so
    // we
    // can set it later
    axis2.load();
    if (static_cast<size_t>(axis2.size()) == nspectra + 1)
      verticalHistogram = true;
    m_axis1vals = MantidVec(axis2(), axis2() + axis2.dim0());
  } catch (std::runtime_error &) {
    g_log.information() << "Axis 0 set to unitless quantity \"" << unit1 << "\"\n";
  }
 
  // Setting a unit onto a TextAxis makes no sense.
  if (unit2 == "TextAxis") {
    auto newAxis = std::make_unique<Mantid::API::TextAxis>(nspectra);
    local_workspace->replaceAxis(1, std::move(newAxis));
  } else if (unit2 != "spectraNumber") {
    try {
      auto newAxis = (verticalHistogram) ? std::make_unique<API::BinEdgeAxis>(nspectra + 1)
                                         : std::make_unique<API::NumericAxis>(nspectra);
      auto newAxisRaw = newAxis.get();
      local_workspace->replaceAxis(1, std::move(newAxis));
      newAxisRaw->unit() = UnitFactory::Instance().create(unit2);
      if (unit2 == "Label") {
        auto label = std::dynamic_pointer_cast<Mantid::Kernel::Units::Label>(newAxisRaw->unit());
        auto ax = wksp_cls.openNXDouble("axis2");
        label->setLabel(ax.attributes("caption"), ax.attributes("label"));
      }
    } catch (std::runtime_error &) {
      g_log.information() << "Axis 1 set to unitless quantity \"" << unit2 << "\"\n";
    }
  }
 
  // Are we a distribution
  std::string dist = xbins.attributes("distribution");
  if (dist == "1") {
    local_workspace->setDistribution(true);
  } else {
    local_workspace->setDistribution(false);
  }
 
  // Get information from all but data group
  std::string parameterStr;
 
  progress(progressStart + 0.05 * progressRange, "Reading the sample details...");
 
  // Hop to the right point
  m_nexusFile->openPath(mtd_entry.path());
  try {
    // This loads logs, sample, and instrument.
    local_workspace->loadExperimentInfoNexus(getPropertyValue("Filename"), m_nexusFile.get(),
                                             parameterStr); // REQUIRED PER PERIOD
 
    // Parameter map parsing only if instrument loaded OK.
    progress(progressStart + 0.11 * progressRange, "Reading the parameter maps...");
    local_workspace->readParameterMap(parameterStr);
  } catch (std::exception &e) {
    // TODO. For workspaces saved via SaveNexusESS, these warnings are not
    // relevant. Unfortunately we need to close all file handles before we can
    // attempt loading the new way see loadNexusGeometry function . A better
    // solution should be found
    g_log.warning("Error loading Instrument section of nxs file");
    g_log.warning(e.what());
    g_log.warning("Try running LoadInstrument Algorithm on the Workspace to "
                  "update the geometry");
  }
 
  readSpectraToDetectorMapping(mtd_entry, *local_workspace);
 
  if (!local_workspace->getAxis(1)->isSpectra()) { // If not a spectra axis, load the axis data into
                                                   // the workspace. (MW 25/11/10)
    loadNonSpectraAxis(local_workspace, wksp_cls);
  }
 
  progress(progressStart + 0.15 * progressRange, "Reading the workspace history...");
  m_nexusFile->openPath(mtd_entry.path());
  try {
    bool load_history = getProperty("LoadHistory");
    if (load_history)
      local_workspace->history().loadNexus(m_nexusFile.get());
  } catch (std::out_of_range &) {
    g_log.warning() << "Error in the workspaces algorithm list, its processing "
                       "history is incomplete\n";
  }
 
  progress(progressStart + 0.2 * progressRange, "Reading the workspace history...");
 
  try {
    if (local_workspace->getTitle().empty())
      local_workspace->setTitle(mtd_entry.getString("title"));
  } catch (std::runtime_error &) {
    g_log.debug() << "No title was found in the input file, " << getPropertyValue("Filename") << '\n';
  }
 
  return std::static_pointer_cast<API::Workspace>(local_workspace);
}
 
//-------------------------------------------------------------------------------------------------
void LoadNexusProcessed::readInstrumentGroup(NXEntry &mtd_entry, API::MatrixWorkspace &local_workspace) {
  // Get spectrum information for the current entry.
 
  SpectraInfo spectraInfo = extractMappingInfo(mtd_entry, this->g_log);
 
  // Now build the spectra list
  int index = 0;
  bool haveSpectraAxis = local_workspace.getAxis(1)->isSpectra();
 
  for (int i = 1; i <= spectraInfo.nSpectra; ++i) {
    int spectrum;
    // prefer the spectra number from the instrument section
    // over anything else. If not there then use a spectra axis
    // number if we have one, else make one up as nothing was
    // written to the file. We should always set it so that
    // CompareWorkspaces gives the expected answer on a Save/Load
    // round trip.
    if (spectraInfo.hasSpectra) {
      spectrum = spectraInfo.spectraNumbers[i - 1];
    } else if (haveSpectraAxis && !m_axis1vals.empty()) {
      spectrum = static_cast<specnum_t>(m_axis1vals[i - 1]);
    } else {
      spectrum = i + 1;
    }
 
    if ((i >= m_spec_min && i < m_spec_max) ||
        (m_list && find(m_spec_list.begin(), m_spec_list.end(), i) != m_spec_list.end())) {
      auto &spec = local_workspace.getSpectrum(index);
      spec.setSpectrumNo(spectrum);
      ++index;
 
      if (!spectraInfo.detectorIndex.empty()) {
        const int start = spectraInfo.detectorIndex[i - 1];
        const int end = start + spectraInfo.detectorCount[i - 1];
        spec.setDetectorIDs(
            std::set<detid_t>(spectraInfo.detectorList.data() + start, spectraInfo.detectorList.data() + end));
      }
    }
  }
}
 
std::map<std::string, std::string> LoadNexusProcessed::validateInputs() {
  using namespace std;
  map<string, string> errorList;
 
  int specMin = getProperty("SpectrumMin");
  int specMax = getProperty("SpectrumMax");
 
  // Check our range is not reversed
  if (specMax < specMin) {
    errorList["SpectrumMin"] = "SpectrumMin must be smaller than SpectrumMax";
    errorList["SpectrumMax"] = "SpectrumMax must be larger than SpectrumMin";
  }
 
  // Finished testing return any errors
  return errorList;
}
 
//-------------------------------------------------------------------------------------------------
void LoadNexusProcessed::loadNonSpectraAxis(const API::MatrixWorkspace_sptr &local_workspace, NXData &data) {
  Axis *axis = local_workspace->getAxis(1);
 
  if (axis->isNumeric()) {
    NXDouble axisData = data.openNXDouble("axis2");
    axisData.load();
    for (int i = 0; i < static_cast<int>(axis->length()); i++) {
      axis->setValue(i, axisData[i]);
    }
  } else if (axis->isText()) {
    NXChar axisData = data.openNXChar("axis2");
    std::string axisLabels;
    try {
      axisData.load();
      axisLabels = std::string(axisData(), axisData.dim0());
    } catch (std::runtime_error &) {
      axisLabels = "";
    }
    // Use boost::tokenizer to split up the input
    Mantid::Kernel::StringTokenizer tokenizer(axisLabels, "\n", Mantid::Kernel::StringTokenizer::TOK_IGNORE_EMPTY);
    // We must cast the axis object to TextAxis so we may use ->setLabel
    auto *textAxis = static_cast<TextAxis *>(axis);
    int i = 0;
    for (const auto &tokIter : tokenizer) {
      textAxis->setLabel(i, tokIter);
      i++;
    }
  }
}
 
bool UDlesserExecCount(const NXClassInfo &elem1, const NXClassInfo &elem2) {
  std::string::size_type index1, index2;
  std::string num1, num2;
  // find the number after "_" in algorithm name ( eg:MantidAlogorthm_1)
  index1 = elem1.nxname.find('_');
  if (index1 != std::string::npos) {
    num1 = elem1.nxname.substr(index1 + 1, elem1.nxname.length() - index1);
  }
  index2 = elem2.nxname.find('_');
  if (index2 != std::string::npos) {
    num2 = elem2.nxname.substr(index2 + 1, elem2.nxname.length() - index2);
  }
  std::stringstream is1, is2;
  is1 << num1;
  is2 << num2;
 
  int execNum1 = -1;
  int execNum2 = -1;
  is1 >> execNum1;
  is2 >> execNum2;
 
  return execNum1 < execNum2;
}
 
//-------------------------------------------------------------------------------------------------
void LoadNexusProcessed::getWordsInString(const std::string &words3, std::string &w1, std::string &w2,
                                          std::string &w3) {
  Mantid::Kernel::StringTokenizer data(words3, " ", Mantid::Kernel::StringTokenizer::TOK_TRIM);
  if (data.count() != 3) {
    g_log.warning() << "Algorithm list line " + words3 + " is not of the correct format\n";
    throw std::out_of_range(words3);
  }
 
  w1 = data[0];
  w2 = data[1];
  w3 = data[2];
}
 
//-------------------------------------------------------------------------------------------------
void LoadNexusProcessed::getWordsInString(const std::string &words4, std::string &w1, std::string &w2, std::string &w3,
                                          std::string &w4) {
  Mantid::Kernel::StringTokenizer data(words4, " ", Mantid::Kernel::StringTokenizer::TOK_TRIM);
  if (data.count() != 4) {
    g_log.warning() << "Algorithm list line " + words4 + " is not of the correct format\n";
    throw std::out_of_range(words4);
  }
 
  w1 = data[0];
  w2 = data[1];
  w3 = data[2];
  w4 = data[3];
}
 
//-------------------------------------------------------------------------------------------------
void LoadNexusProcessed::readBinMasking(const NXData &wksp_cls, const API::MatrixWorkspace_sptr &local_workspace) {
  if (wksp_cls.getDataSetInfo("masked_spectra").stat == NX_ERROR) {
    return;
  }
  NXInt spec = wksp_cls.openNXInt("masked_spectra");
  spec.load();
  NXSize bins = wksp_cls.openNXSize("masked_bins");
  bins.load();
  NXDouble weights = wksp_cls.openNXDouble("mask_weights");
  weights.load();
  const int n = spec.dim0();
  const int n1 = n - 1;
  for (int i = 0; i < n; ++i) {
    int si = spec(i, 0);
    int j0 = spec(i, 1);
    int j1 = i < n1 ? spec(i + 1, 1) : bins.dim0();
    for (int j = j0; j < j1; ++j) {
      local_workspace->flagMasked(si, bins[j], weights[j]);
    }
  }
}
 
void LoadNexusProcessed::loadBlock(NXDataSetTyped<double> &data, NXDataSetTyped<double> &errors,
                                   NXDataSetTyped<double> &farea, bool hasFArea, NXDouble &xErrors, bool hasXErrors,
                                   int blocksize, int nchannels, int &hist,
                                   const API::MatrixWorkspace_sptr &local_workspace) {
  data.load(blocksize, hist);
  errors.load(blocksize, hist);
  double *data_start = data();
  double *data_end = data_start + nchannels;
  double *err_start = errors();
  double *err_end = err_start + nchannels;
  double *farea_start = nullptr;
  double *farea_end = nullptr;
  double *xErrors_start = nullptr;
  double *xErrors_end = nullptr;
  size_t dx_increment = nchannels;
  // NexusFileIO stores Dx data for all spectra (sharing not preserved) so dim0
  // is the histograms, dim1 is Dx length. For old files this is nchannels+1,
  // otherwise nchannels. See #16298.
  // WARNING: We are dropping the last Dx value for old files!
  size_t dx_input_increment = xErrors.dim1();
  RebinnedOutput_sptr rb_workspace;
  if (hasFArea) {
    farea.load(blocksize, hist);
    farea_start = farea();
    farea_end = farea_start + nchannels;
    rb_workspace = std::dynamic_pointer_cast<RebinnedOutput>(local_workspace);
  }
  if (hasXErrors) {
    xErrors.load(blocksize, hist);
    xErrors_start = xErrors();
    xErrors_end = xErrors_start + dx_increment;
  }
 
  int final(hist + blocksize);
  while (hist < final) {
    auto &Y = local_workspace->mutableY(hist);
    Y.assign(data_start, data_end);
    data_start += nchannels;
    data_end += nchannels;
    auto &E = local_workspace->mutableE(hist);
    E.assign(err_start, err_end);
    err_start += nchannels;
    err_end += nchannels;
    if (hasFArea) {
      MantidVec &F = rb_workspace->dataF(hist);
      F.assign(farea_start, farea_end);
      farea_start += nchannels;
      farea_end += nchannels;
    }
    if (hasXErrors) {
      local_workspace->setSharedDx(hist, Kernel::make_cow<HistogramData::HistogramDx>(xErrors_start, xErrors_end));
      xErrors_start += dx_input_increment;
      xErrors_end += dx_input_increment;
    }
 
    local_workspace->setSharedX(hist, m_xbins.cowData());
    ++hist;
  }
}
 
void LoadNexusProcessed::loadBlock(NXDataSetTyped<double> &data, NXDataSetTyped<double> &errors,
                                   NXDataSetTyped<double> &farea, bool hasFArea, NXDouble &xErrors, bool hasXErrors,
                                   int blocksize, int nchannels, int &hist, int &wsIndex,
                                   const API::MatrixWorkspace_sptr &local_workspace) {
  data.load(blocksize, hist);
  errors.load(blocksize, hist);
  double *data_start = data();
  double *data_end = data_start + nchannels;
  double *err_start = errors();
  double *err_end = err_start + nchannels;
  double *farea_start = nullptr;
  double *farea_end = nullptr;
  double *xErrors_start = nullptr;
  double *xErrors_end = nullptr;
  size_t dx_increment = nchannels;
  // NexusFileIO stores Dx data for all spectra (sharing not preserved) so dim0
  // is the histograms, dim1 is Dx length. For old files this is nchannels+1,
  // otherwise nchannels. See #16298.
  // WARNING: We are dropping the last Dx value for old files!
  size_t dx_input_increment = xErrors.dim1();
  RebinnedOutput_sptr rb_workspace;
  if (hasFArea) {
    farea.load(blocksize, hist);
    farea_start = farea();
    farea_end = farea_start + nchannels;
    rb_workspace = std::dynamic_pointer_cast<RebinnedOutput>(local_workspace);
  }
  if (hasXErrors) {
    xErrors.load(blocksize, hist);
    xErrors_start = xErrors();
    xErrors_end = xErrors_start + dx_increment;
  }
 
  int final(hist + blocksize);
  while (hist < final) {
    auto &Y = local_workspace->mutableY(wsIndex);
    Y.assign(data_start, data_end);
    data_start += nchannels;
    data_end += nchannels;
    auto &E = local_workspace->mutableE(wsIndex);
    E.assign(err_start, err_end);
    err_start += nchannels;
    err_end += nchannels;
    if (hasFArea) {
      MantidVec &F = rb_workspace->dataF(wsIndex);
      F.assign(farea_start, farea_end);
      farea_start += nchannels;
      farea_end += nchannels;
    }
    if (hasXErrors) {
      local_workspace->setSharedDx(wsIndex, Kernel::make_cow<HistogramData::HistogramDx>(xErrors_start, xErrors_end));
      xErrors_start += dx_input_increment;
      xErrors_end += dx_input_increment;
    }
    local_workspace->setSharedX(wsIndex, m_xbins.cowData());
    ++hist;
    ++wsIndex;
  }
}
 
void LoadNexusProcessed::loadBlock(NXDataSetTyped<double> &data, NXDataSetTyped<double> &errors,
                                   NXDataSetTyped<double> &farea, bool hasFArea, NXDouble &xErrors, bool hasXErrors,
                                   NXDouble &xbins, int blocksize, int nchannels, int &hist, int &wsIndex,
                                   const API::MatrixWorkspace_sptr &local_workspace) {
  data.load(blocksize, hist);
  double *data_start = data();
  double *data_end = data_start + nchannels;
  errors.load(blocksize, hist);
  double *err_start = errors();
  double *err_end = err_start + nchannels;
  double *farea_start = nullptr;
  double *farea_end = nullptr;
  double *xErrors_start = nullptr;
  double *xErrors_end = nullptr;
  size_t dx_increment = nchannels;
  // NexusFileIO stores Dx data for all spectra (sharing not preserved) so dim0
  // is the histograms, dim1 is Dx length. For old files this is nchannels+1,
  // otherwise nchannels. See #16298.
  // WARNING: We are dropping the last Dx value for old files!
  size_t dx_input_increment = xErrors.dim1();
  RebinnedOutput_sptr rb_workspace;
  if (hasFArea) {
    farea.load(blocksize, hist);
    farea_start = farea();
    farea_end = farea_start + nchannels;
    rb_workspace = std::dynamic_pointer_cast<RebinnedOutput>(local_workspace);
  }
  xbins.load(blocksize, hist);
  const int nxbins(xbins.dim1());
  double *xbin_start = xbins();
  double *xbin_end = xbin_start + nxbins;
  int final(hist + blocksize);
 
  if (hasXErrors) {
    xErrors.load(blocksize, hist);
    xErrors_start = xErrors();
    xErrors_end = xErrors_start + dx_increment;
  }
 
  while (hist < final) {
    auto &Y = local_workspace->mutableY(wsIndex);
    Y.assign(data_start, data_end);
    data_start += nchannels;
    data_end += nchannels;
    auto &E = local_workspace->mutableE(wsIndex);
    E.assign(err_start, err_end);
    err_start += nchannels;
    err_end += nchannels;
    if (hasFArea) {
      MantidVec &F = rb_workspace->dataF(wsIndex);
      F.assign(farea_start, farea_end);
      farea_start += nchannels;
      farea_end += nchannels;
    }
    if (hasXErrors) {
      local_workspace->setSharedDx(wsIndex, Kernel::make_cow<HistogramData::HistogramDx>(xErrors_start, xErrors_end));
      xErrors_start += dx_input_increment;
      xErrors_end += dx_input_increment;
    }
    auto &X = local_workspace->mutableX(wsIndex);
    X.assign(xbin_start, xbin_end);
    xbin_start += nxbins;
    xbin_end += nxbins;
    ++hist;
    ++wsIndex;
  }
}
 
void LoadNexusProcessed::checkOptionalProperties(const std::size_t numberofspectra) {
  // read in the settings passed to the algorithm
  m_spec_list = getProperty("SpectrumList");
  m_spec_max = getProperty("SpectrumMax");
  m_spec_min = getProperty("SpectrumMin");
  // Are we using a list of spectra or all the spectra in a range?
  m_list = !m_spec_list.empty();
  m_interval = (m_spec_max != Mantid::EMPTY_INT()) || (m_spec_min != 1);
  if (m_spec_max == Mantid::EMPTY_INT())
    m_spec_max = 1;
 
  // Check validity of spectra list property, if set
  if (m_list) {
    const int minlist = *min_element(m_spec_list.begin(), m_spec_list.end());
    const int maxlist = *max_element(m_spec_list.begin(), m_spec_list.end());
    // Need to check before casting
    if (maxlist < 0) {
      g_log.error("Invalid list of spectra");
      throw std::invalid_argument("Spectra max is less than 0");
    }
 
    if (maxlist > static_cast<int>(numberofspectra) || minlist == 0) {
      g_log.error("Invalid list of spectra");
      throw std::invalid_argument("Inconsistent properties defined");
    }
  }
 
  // Check validity of spectra range, if set
  if (m_interval) {
    m_interval = true;
    m_spec_min = getProperty("SpectrumMin");
    if (m_spec_min != 1 && m_spec_max == 1) {
      m_spec_max = static_cast<int>(numberofspectra);
    }
    if (m_spec_max < m_spec_min || m_spec_max > static_cast<int>(numberofspectra)) {
      g_log.error("Invalid Spectrum min/max properties");
      throw std::invalid_argument("Inconsistent properties defined");
    }
  }
}
 
size_t LoadNexusProcessed::calculateWorkspaceSize(const std::size_t numberofspectra, bool gen_filtered_list) {
  // Calculate the size of a workspace, given its number of spectra to read
  size_t total_specs;
  if (m_interval || m_list) {
    if (m_interval) {
      if (m_spec_min != 1 && m_spec_max == 1) {
        m_spec_max = static_cast<int>(numberofspectra);
      }
      total_specs = m_spec_max - m_spec_min + 1;
      m_spec_max += 1;
 
      if (gen_filtered_list) {
        m_filtered_spec_idxs.resize(total_specs);
        size_t j = 0;
        for (int si = m_spec_min; si < m_spec_max; si++, j++)
          m_filtered_spec_idxs[j] = si;
      }
    } else {
      total_specs = 0;
    }
 
    if (m_list) {
      if (m_interval) {
        for (auto it = m_spec_list.begin(); it != m_spec_list.end();)
          if (*it >= m_spec_min && *it < m_spec_max) {
            it = m_spec_list.erase(it);
          } else
            ++it;
      }
      if (m_spec_list.empty())
        m_list = false;
      total_specs += m_spec_list.size();
 
      if (gen_filtered_list) {
        // range list + spare indices from list
        // example: min: 2, max: 8, list: 3,4,5,10,12;
        //          result: 2,3,...,7,8,10,12
        m_filtered_spec_idxs.insert(m_filtered_spec_idxs.end(), m_spec_list.begin(), m_spec_list.end());
      }
    }
  } else {
    total_specs = numberofspectra;
    m_spec_min = 1;
    m_spec_max = static_cast<int>(numberofspectra) + 1;
 
    if (gen_filtered_list) {
      m_filtered_spec_idxs.resize(total_specs, 0);
      for (int j = 0; j < static_cast<int>(total_specs); j++)
        m_filtered_spec_idxs[j] = m_spec_min + j;
    }
  }
  return total_specs;
}
 
void LoadNexusProcessed::applyLogFiltering(const Mantid::API::Workspace_sptr &local_workspace) {
  auto mWorkspace = std::dynamic_pointer_cast<MatrixWorkspace>(local_workspace);
  if (mWorkspace) {
    auto run = mWorkspace->run();
    // check for presence of filterable logs that suggest this is ISIS data
    if (run.hasProperty(LogParser::statusLogName()) || run.hasProperty(LogParser::periodsLogName())) {
      ISISRunLogs::applyLogFiltering(mWorkspace->mutableRun());
    }
  }
}
 
} // namespace Mantid::DataHandling