LoadNexusProcessed.cpp

//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/AlgorithmFactory.h"
#include "MantidAPI/AlgorithmManager.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/BinEdgeAxis.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/TextAxis.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidDataHandling/LoadNexusProcessed.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/RebinnedOutput.h"
#include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidDataObjects/PeakNoShapeFactory.h"
#include "MantidDataObjects/PeakShapeSphericalFactory.h"
#include "MantidDataObjects/PeakShapeEllipsoidFactory.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/DateAndTime.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidNexus/NexusClasses.h"
#include "MantidNexus/NexusFileIO.h"

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

#include <MantidKernel/StringTokenizer.h>

#include <nexus/NeXusException.hpp>

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;

namespace {

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

// Struct to contain spectrum information.
struct SpectraInfo {
  // Number of spectra
  int nSpectra;
  // Do we have any spectra
  bool hasSpectra;
  // 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;

  SpectraInfo() : nSpectra(0), hasSpectra(false) {}

  SpectraInfo(int _nSpectra, bool _hasSpectra, IntArray_shared _spectraNumbers,
              IntArray_shared _detectorIndex, IntArray_shared _detectorCount,
              IntArray_shared _detectorList)
      : nSpectra(_nSpectra), hasSpectra(_hasSpectra),
        spectraNumbers(_spectraNumbers), detectorIndex(_detectorIndex),
        detectorCount(_detectorCount), detectorList(_detectorList) {}

  SpectraInfo(const SpectraInfo &other)
      : nSpectra(other.nSpectra), hasSpectra(other.hasSpectra),
        spectraNumbers(other.spectraNumbers),
        detectorIndex(other.detectorIndex), detectorCount(other.detectorCount),
        detectorList(other.detectorList) {}

  /*
  SpectraInfo& operator=(const SpectraInfo& other)
  {
    if (&other != this)
    {
      nSpectra = other.nSpectra;
      hasSpectra = other.hasSpectra;
      spectraNumbers = other.spectraNumbers;
      detectorIndex = other.detectorIndex;
      detectorCount = other.detectorCount;
      detectorList = other.detectorList;
    }
    return *this;
  }
  */
};

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

/**
 * Extract ALL the detector, spectrum number and workspace index mapping
 * information.
 * @param mtd_entry
 * @param logger
 * @return
 */
SpectraInfo extractMappingInfo(NXEntry &mtd_entry, Logger &logger) {
  // Instrument information
  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");

  // 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
  boost::shared_array<int> detectorList;
  try {
    NXInt detlist_group = detgroup.openNXInt("detector_list");
    detlist_group.load();
    detectorList = detlist_group.sharedBuffer();
  } catch (std::runtime_error &) {
    logger.information() << "detector_list block not found. The workspace will "
                            "not contain any detector information."
                         << std::endl;
    return SpectraInfo();
  }

  // Detector count contains the number of detectors associated with each
  // spectra
  NXInt det_count = detgroup.openNXInt("detector_count");
  det_count.load();
  boost::shared_array<int> 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();
  int nspectra = det_index.dim0();
  boost::shared_array<int> detectorIndex = det_index.sharedBuffer();

  // 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
  bool have_spectra(true);
  boost::shared_array<int> spectra;
  try {
    NXInt spectra_block = detgroup.openNXInt("spectra");
    spectra_block.load();
    spectra = spectra_block.sharedBuffer();
  } catch (std::runtime_error &) {
    have_spectra = false;
  }
  return SpectraInfo(nspectra, have_spectra, spectra, detectorIndex,
                     detectorCount, detectorList);
}

/**
 * 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, 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 int nPeriods =
          run.getPropertyValueAsType<int>(nPeriodsLogEntryName);
      if (nPeriods == nWorkspaceEntries) {
        isMultiPeriod = true;
        log.information("Loading as MultiPeriod group workspace.");
      }
    }
  }
  return isMultiPeriod;
}
}

/// Default constructor
LoadNexusProcessed::LoadNexusProcessed()
    : m_shared_bins(false), m_xbins(), 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_cppFile(nullptr) {}

/// Delete NexusFileIO in destructor
LoadNexusProcessed::~LoadNexusProcessed() { delete m_cppFile; }

/**
 * 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;
}

/** Initialisation method.
 *
 */
void LoadNexusProcessed::init() {
  // Declare required input parameters for algorithm
  declareProperty(
      new FileProperty("Filename", "", FileProperty::Load,
                       {".nxs", ".nx5", ".xml"}),
      "The name of the Nexus file to read, as a full or relative path.");
  declareProperty(new 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<int64_t>>();
  mustBePositive->setLower(0);

  declareProperty("SpectrumMin", static_cast<int64_t>(1), mustBePositive,
                  "Number of first spectrum to read.");
  declareProperty("SpectrumMax", static_cast<int64_t>(Mantid::EMPTY_INT()),
                  mustBePositive, "Number of last spectrum to read.");
  declareProperty(new ArrayProperty<int64_t>("SpectrumList"),
                  "List of spectrum numbers to read.");
  declareProperty("EntryNumber", static_cast<int64_t>(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 worskspace");
  declareProperty("LoadHistory", true,
                  "If true, the workspace history will be loaded");
  declareProperty(
      new 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 int64_t nWorkspaceEntries,
    const int64_t p) {

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

  const size_t nHistograms = periodWorkspace->getNumberHistograms();
  for (size_t i = 0; i < nHistograms; ++i) {
    periodWorkspace->setX(i, tempMatrixWorkspace->refX(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.
  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."
        << std::endl;
    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." << std::endl;
    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();

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

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

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

    data.load(static_cast<int>(blockSize), static_cast<int>(histIndex));
    errors.load(static_cast<int>(blockSize), static_cast<int>(histIndex));

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

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

    int64_t final(histIndex + blockSize);
    while (histIndex < final) {
      MantidVec &Y = periodWorkspace->dataY(wsIndex);
      Y.assign(dataStart, dataEnd);
      dataStart += nChannels;
      dataEnd += nChannels;
      MantidVec &E = periodWorkspace->dataE(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_cppFile->closeGroup();

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

  try {
    // This loads logs, sample, and instrument.
    periodWorkspace->loadSampleAndLogInfoNexus(m_cppFile);
  } 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() {
  progress(0, "Opening file...");

  // Throws an approriate exception if there is a problem with file access
  NXRoot root(getPropertyValue("Filename"));

  // "Open" the same file but with the C++ interface
  m_cppFile = new ::NeXus::File(root.m_fileID);

  // Find out how many first level entries there are
  int64_t nWorkspaceEntries = static_cast<int64_t>(root.groups().size());

  // Check for an entry number property
  int64_t entrynumber = static_cast<int64_t>(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.
  API::Workspace_sptr 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(static_cast<int>(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 and check for a common stem
    std::vector<std::string> names(nWorkspaceEntries + 1);
    bool commonStem = bIsMultiPeriod || checkForCommonNameStem(root, names);

    // 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;
      tempMatrixWorkspace->mutableRun().clearLogs(); // Strip out any loaded
                                                     // logs. That way we don't
                                                     // pay for copying that
                                                     // information around.
    }

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

    for (int64_t p = 1; p <= nWorkspaceEntries; ++p) {
      std::ostringstream os;
      os << p;

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

      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 + os.str(), tempMatrixWorkspace, nWorkspaceEntries,
            p);
      } else // Fall-back for generic loading
      {
        const double nWorkspaceEntries_d =
            static_cast<double>(nWorkspaceEntries);
        local_workspace =
            loadEntry(root, basename + os.str(),
                      static_cast<double>(p - 1) / nWorkspaceEntries_d,
                      1. / nWorkspaceEntries_d);
      }

      declareProperty(new WorkspaceProperty<API::Workspace>(
          prop_name + os.str(), wsName, Direction::Output));

      wksp_group->addWorkspace(local_workspace);
      setProperty(prop_name + os.str(), local_workspace);
    }

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

  m_axis1vals.clear();
}

/**
 * Decides what to call a child of a group workspace.
 *
 * This function uses information about if the child workspace has a common stem
 * and checks if the file contained a workspace name to decide what it should be
 *called
 *
 * @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
 * @param commonStem :: Whether the workspaces share a common name stem
 *
 * @return The name of the workspace
 */
std::string LoadNexusProcessed::buildWorkspaceName(const std::string &name,
                                                   const std::string &baseName,
                                                   int64_t wsIndex,
                                                   bool commonStem) {
  std::string wsName;
  std::string index = boost::lexical_cast<std::string>(wsIndex);

  // if we don't have a common stem then use name tag
  if (!commonStem) {
    if (!name.empty()) {
      // use name loaded from file there's no common stem
      wsName = name;
    } else {
      // if the name property wasn't defined just use <OutputWorkspaceName>_n
      wsName = baseName + index;
    }
  } 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 {
      // use default name if we couldn't match for some reason
      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 = "_" + boost::lexical_cast<std::string>(i);
    }

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

/**
 * Check if the workspace name contains a common stem and load the workspace
 *names
 *
 * @param root :: the root for the NeXus document
 * @param names :: vector to store the names to be loaded.
 * @return Whether there was a common stem.
 */
bool LoadNexusProcessed::checkForCommonNameStem(
    NXRoot &root, std::vector<std::string> &names) {
  bool success(true);
  int64_t nWorkspaceEntries = static_cast<int64_t>(root.groups().size());
  for (int64_t p = 1; p <= nWorkspaceEntries; ++p) {
    std::ostringstream os;
    os << p;

    names[p] = loadWorkspaceName(root, "mantid_workspace_" + os.str());

    boost::smatch results;
    const boost::regex exp(".*_\\d+$");

    // check if the workspace name has an index on the end
    if (!boost::regex_match(names[p], results, exp)) {
      success = false;
    }
  }

  return success;
}

/**
 * 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);
  try {
    return mtd_entry.getString("workspace_name");
  } catch (std::runtime_error &) {
    return std::string();
  }
}

//-------------------------------------------------------------------------------------------------
/**
 * 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 (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
  PARALLEL_FOR_NO_WSP_CHECK()
  for (int64_t j = 0; j < static_cast<int64_t>(m_filtered_spec_idxs.size());
       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->getEventList(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.setX(this->m_xbins);
      else {
        MantidVec x;
        x.resize(xbins.dim0());
        for (int i = 0; i < xbins.dim0(); i++)
          x[i] = xbins(static_cast<int>(wi), i);
        el.setX(x);
      }
    }

    progress(progressStart +
             progressRange * (1.0 / static_cast<double>(numspec)));
    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_" + boost::lexical_cast<std::string>(columnNumber);

    NXInfo info = nx_tw.getDataSetInfo(dataSetName.c_str());
    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 " +
                               boost::lexical_cast<std::string>(info.type) +
                               " cannot be loaded.");
      }
    } else if (info.rank == 2) {
      if (info.type == NX_CHAR) {
        NXChar data = nx_tw.openNXChar(dataSetName.c_str());
        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.c_str());
        if (data.attributes("interpret_as") == "V3D") {
          loadV3DColumn(data, workspace);
        } else {
          loadVectorColumn<double>(nx_tw, dataSetName, workspace,
                                   "vector_double");
        }
      }
    }

    columnNumber++;

  } while (1);

  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.c_str());
  std::string columnTitle = data.attributes("name");
  if (!columnTitle.empty()) {
    tableWs->addColumn(columnType, columnTitle);

    NXInfo info = tableData.getDataSetInfo(dataSetName.c_str());
    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_" + boost::lexical_cast<std::string>(columnNumber);

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