LoadHelper.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 +
/*
 * Helper file to gather common routines to the Loaders
 * */
 
#include "MantidDataHandling/LoadHelper.h"
#include "MantidAPI/AlgorithmManager.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/ComponentInfo.h"
#include "MantidKernel/OptionalBool.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidNexus/NexusClasses.h"
#include "MantidNexus/NexusException.h"
#include "MantidNexus/NexusFile.h"
#include <strstream>
 
namespace Mantid {
 
namespace {
Kernel::Logger g_log("LoadHelper");
} // namespace
 
namespace DataHandling {
using namespace Kernel;
using namespace API;
 
std::string LoadHelper::findInstrumentNexusAddress(const Mantid::Nexus::NXEntry &firstEntry) {
  std::string result("");
  std::vector<Mantid::Nexus::NXClassInfo> v = firstEntry.groups();
  const auto it = std::find_if(v.cbegin(), v.cend(), [](const auto &group) { return group.nxclass == "NXinstrument"; });
  if (it != v.cend())
    result = it->nxname;
 
  return result;
}
 
std::string LoadHelper::getStringFromNexusAddress(const Mantid::Nexus::NXEntry &firstEntry,
                                                  const std::string &nexusAddress) {
  return firstEntry.getString(nexusAddress);
}
 
double LoadHelper::getDoubleFromNexusAddress(const Mantid::Nexus::NXEntry &firstEntry,
                                             const std::string &nexusAddress) {
  return firstEntry.getFloat(nexusAddress);
}
 
std::vector<double> LoadHelper::getTimeBinningFromNexusAddress(const Mantid::Nexus::NXEntry &firstEntry,
                                                               const std::string &nexusAddress) {
 
  Mantid::Nexus::NXFloat timeBinningNexus = firstEntry.openNXFloat(nexusAddress);
  timeBinningNexus.load();
 
  size_t numberOfBins = static_cast<size_t>(timeBinningNexus.dim0()) + 1; // boundaries
 
  float *timeBinning_p = &timeBinningNexus[0];
  std::vector<double> timeBinning(numberOfBins);
  std::copy(timeBinning_p, timeBinning_p + timeBinningNexus.dim0(), std::begin(timeBinning));
  // calculate the extra bin at the end
  timeBinning[numberOfBins - 1] = timeBinning[numberOfBins - 2] + timeBinning[1] - timeBinning[0];
 
  return timeBinning;
}
double LoadHelper::calculateEnergy(double wavelength) {
  double e = (PhysicalConstants::h * PhysicalConstants::h) /
             (2 * PhysicalConstants::NeutronMass * wavelength * wavelength * 1e-20) / PhysicalConstants::meV;
  return e;
}
 
double LoadHelper::calculateTOF(double distance, double wavelength) {
  if (wavelength <= 0) {
    throw std::runtime_error("Wavelenght is <= 0");
  }
 
  double velocity = PhysicalConstants::h / (PhysicalConstants::NeutronMass * wavelength * 1e-10); // m/s
 
  return distance / velocity;
}
 
/*
 * Get instrument property as double
 * @s - input property name
 *
 */
double LoadHelper::getInstrumentProperty(const API::MatrixWorkspace_sptr &workspace, const std::string &s) {
  std::vector<std::string> prop = workspace->getInstrument()->getStringParameter(s);
  if (prop.empty()) {
    return EMPTY_DBL();
  } else {
    return boost::lexical_cast<double>(prop[0]);
  }
}
 
void LoadHelper::addNexusFieldsToWsRun(Nexus::File &filehandle, API::Run &runDetails, const std::string &entryName,
                                       bool useFullAddress) {
  // As a workaround against some "not so good" old ILL nexus files (ILLIN5_Vana_095893.nxs for example) by default we
  // begin the parse on the first entry (entry0), or from a chosen entryName. This allow to avoid the bogus entries that
  // follows
  std::string entryNameActual(entryName);
  if (entryName.empty()) {
    entryNameActual = filehandle.getTopLevelEntryName();
  }
 
  // open the group and parse down
  if (!entryNameActual.empty()) {
    constexpr int LEVEL_RECURSE{1};
    filehandle.openGroup(entryNameActual, "NXentry");
    const std::string EMPTY_STR;
    recurseAndAddNexusFieldsToWsRun(filehandle, runDetails, EMPTY_STR, EMPTY_STR, LEVEL_RECURSE, useFullAddress);
    filehandle.closeGroup();
  }
}
 
namespace {
template <typename NumericType>
void addNumericProperty(Nexus::File &filehandle, const Nexus::Info &nxinfo, const std::string &property_name,
                        API::Run &runDetails) {
  if (!runDetails.hasProperty(property_name)) {
    g_log.warning() << "Property " << property_name << " was set twice. Please check the Nexus file and your inputs.";
  }
  // this assumes all values are rank==1
 
  // Look for "units"
  std::string units; // empty string
  if (filehandle.hasAttr("units"))
    filehandle.getAttr("units", units);
 
  // read the field
  std::vector<NumericType> data_vec;
  data_vec.reserve(nxinfo.dims.front());
  filehandle.getDataCoerce(data_vec);
 
  // create the property on the run objects
  const auto dim1size = data_vec.size();
  if (dim1size == 1) {
    if (units.empty())
      runDetails.addProperty(property_name, data_vec.front());
    else
      runDetails.addProperty(property_name, data_vec.front(), units);
  } else {
    // create a bunch of little properties
    for (size_t index = 0; index < dim1size; ++index) {
      const auto property_name_indexed = property_name + "_" + std::to_string(index);
      if (units.empty())
        runDetails.addProperty(property_name_indexed, data_vec[index]);
      else
        runDetails.addProperty(property_name_indexed, data_vec[index], units);
    }
  }
}
} // namespace
 
void LoadHelper::recurseAndAddNexusFieldsToWsRun(Nexus::File &filehandle, API::Run &runDetails,
                                                 const std::string &parent_name, const std::string &parent_class,
                                                 int level, bool useFullAddress) {
  const std::string SDS("SDS"); // denotes data field
 
  for (const auto &entry : filehandle.getEntries()) {
    const auto nxname = entry.first;
    const auto nxclass = entry.second;
    if (nxclass == SDS) {
      if ((parent_class != "NXData") && parent_class != "NXMonitor" && (nxname != "data")) { // create a property
        filehandle.openData(nxname);
        const auto nxinfo = filehandle.getInfo(); // get information about the open data
        const auto rank = nxinfo.dims.size();
 
        // Note, we choose to only build properties on small float arrays filter logic is below
        bool read_property = (rank == 1);
        switch (nxinfo.type) {
        case (NXnumtype::CHAR):
          break; // everything is fine
        case (NXnumtype::FLOAT32):
        case (NXnumtype::FLOAT64):
          // some 1d float arrays may be loaded
          read_property = (nxinfo.dims[0] <= 9);
          break;
        case (NXnumtype::INT16):
        case (NXnumtype::INT32):
        case (NXnumtype::UINT16):
          // only read scalar values for everything else - e.g. integers
          read_property = (nxinfo.dims.front() == 1);
          break;
        default:
          read_property = false;
        }
 
        if (read_property) {
          // generate a name for the property
          const std::string property_name = (parent_name.empty() ? nxname : parent_name + "." + nxname);
 
          switch (nxinfo.type) {
          case (NXnumtype::CHAR): {
            std::string property_value = filehandle.getStrData();
            if (property_name.ends_with("_time")) {
              // That's a time value! Convert to Mantid standard
              property_value = dateTimeInIsoFormat(property_value);
              if (runDetails.hasProperty(property_name))
                runDetails.getProperty(property_name)->setValue(property_value);
              else
                runDetails.addProperty(property_name, property_value);
            } else {
              if (!runDetails.hasProperty(property_name)) {
                runDetails.addProperty(property_name, property_value);
              } else {
                g_log.warning() << "Property " << property_name
                                << " was set twice. Please check the Nexus file and your inputs.\n";
              }
            }
            break;
          }
          case (NXnumtype::FLOAT32):
          case (NXnumtype::FLOAT64):
            addNumericProperty<double>(filehandle, nxinfo, property_name, runDetails);
            break;
          case (NXnumtype::INT16):
          case (NXnumtype::INT32):
          case (NXnumtype::UINT16):
            addNumericProperty<int>(filehandle, nxinfo, property_name, runDetails);
            break;
          default:
            std::stringstream msg;
            msg << "Encountered unknown type: " << nxinfo.type;
            throw std::runtime_error(msg.str());
            break;
          }
        }
        filehandle.closeData();
      }
    } else {
      if ((nxclass != "ILL_data_scan_vars") && (nxclass != "NXill_data_scan_vars")) {
        // open the group and recurse down, if the group is known to nexus
        filehandle.openGroup(nxname, nxclass);
 
        const std::string property_name = (parent_name.empty() ? nxname : parent_name + "." + nxname);
        const std::string p_nxname =
            useFullAddress ? property_name : nxname; // current names can be useful for next level
        const std::string p_nxclass(nxclass);
 
        recurseAndAddNexusFieldsToWsRun(filehandle, runDetails, p_nxname, p_nxclass, level + 1, useFullAddress);
 
        filehandle.closeGroup();
      }
    }
  }
}
 
std::string LoadHelper::dateTimeInIsoFormat(const std::string &dateToParse) {
  namespace bt = boost::posix_time;
  try {
    // Try to convert to a boost date as an iso string. If it works, we are already in ISO and no conversion is
    // requiered.
    bt::from_iso_extended_string(dateToParse);
    return dateToParse;
 
  } catch (...) {
  }
  // parsing format
  const std::locale format = std::locale(std::locale::classic(), new bt::time_input_facet("%d-%b-%y %H:%M:%S"));
 
  bt::ptime pt;
  std::istringstream is(dateToParse);
  is.imbue(format);
  is >> pt;
 
  if (pt != bt::ptime()) {
    // Converts to ISO
    std::string s = bt::to_iso_extended_string(pt);
    return s;
  } else {
    return "";
  }
}
 
void LoadHelper::moveComponent(const API::MatrixWorkspace_sptr &ws, const std::string &componentName,
                               const V3D &newPos) {
  Geometry::IComponent_const_sptr component = ws->getInstrument()->getComponentByName(componentName);
  if (!component) {
    throw std::invalid_argument("Instrument component " + componentName + " not found");
  }
  auto &componentInfo = ws->mutableComponentInfo();
  const auto componentIndex = componentInfo.indexOf(component->getComponentID());
  componentInfo.setPosition(componentIndex, newPos);
}
 
void LoadHelper::rotateComponent(const API::MatrixWorkspace_sptr &ws, const std::string &componentName,
                                 const Kernel::Quat &rot) {
  Geometry::IComponent_const_sptr component = ws->getInstrument()->getComponentByName(componentName);
  if (!component) {
    throw std::invalid_argument("Instrument component " + componentName + " not found");
  }
  auto &componentInfo = ws->mutableComponentInfo();
  const auto componentIndex = componentInfo.indexOf(component->getComponentID());
  componentInfo.setRotation(componentIndex, rot);
}
 
V3D LoadHelper::getComponentPosition(const API::MatrixWorkspace_sptr &ws, const std::string &componentName) {
  Geometry::IComponent_const_sptr component = ws->getInstrument()->getComponentByName(componentName);
  if (!component) {
    throw std::invalid_argument("Instrument component " + componentName + " not found");
  }
  V3D pos = component->getPos();
  return pos;
}
 
void LoadHelper::loadEmptyInstrument(const API::MatrixWorkspace_sptr &ws, const std::string &instrumentName,
                                     const std::string &instrumentPath) {
  auto loadInst = AlgorithmManager::Instance().create("LoadInstrument");
  loadInst->initialize();
  loadInst->setChild(true);
  loadInst->setPropertyValue("InstrumentName", instrumentName);
  if (!instrumentPath.empty())
    loadInst->setPropertyValue("Filename", instrumentPath);
  loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", ws);
  loadInst->setProperty("RewriteSpectraMap", OptionalBool(true));
  loadInst->execute();
}
 
void LoadHelper::fillStaticWorkspace(const API::MatrixWorkspace_sptr &ws, const Mantid::Nexus::NXInt &data,
                                     const std::vector<double> &xAxis, int64_t initialSpectrum, bool pointData,
                                     const std::vector<detid_t> &detectorIDs,
                                     const std::set<detid_t> &acceptedDetectorIDs,
                                     const std::tuple<short, short, short> &axisOrder) {
 
  const bool customDetectorIDs = detectorIDs.size() != 0;
  const bool excludeDetectorIDs = acceptedDetectorIDs.size() != 0;
 
  std::array dims = {data.dim0(), data.dim1(), data.dim2()};
  const auto nTubes = dims[std::get<0>(axisOrder)];
  const auto nPixels = dims[std::get<1>(axisOrder)];
  const auto nChannels = dims[std::get<2>(axisOrder)];
 
  int loadOrder[3] = {0, 1, 2};
  loadingOrder(axisOrder, loadOrder);
 
  HistogramData::Points histoPoints;
  HistogramData::BinEdges binEdges;
  if (pointData)
    histoPoints = HistogramData::Points(xAxis);
  else
    binEdges = HistogramData::BinEdges(xAxis);
  int nSkipped = 0;
 
#pragma omp parallel for if (!excludeDetectorIDs && Kernel::threadSafe(*ws))
  for (specnum_t tube_no = 0; tube_no < static_cast<specnum_t>(nTubes); ++tube_no) {
    for (specnum_t pixel_no = 0; pixel_no < static_cast<specnum_t>(nPixels); ++pixel_no) {
      specnum_t currentSpectrum =
          static_cast<specnum_t>(initialSpectrum) + tube_no * static_cast<specnum_t>(nPixels) + pixel_no;
      if (excludeDetectorIDs != 0 && std::find(acceptedDetectorIDs.cbegin(), acceptedDetectorIDs.cend(),
                                               currentSpectrum) == acceptedDetectorIDs.end()) {
        nSkipped++;
        continue;
      }
      currentSpectrum -= nSkipped;
 
      std::vector<int> spectrum(nChannels);
      for (specnum_t channel_no = 0; channel_no < static_cast<specnum_t>(nChannels); ++channel_no) {
        const int dataIndices[3] = {tube_no, pixel_no, channel_no};
        spectrum[channel_no] = data(dataIndices[loadOrder[0]], dataIndices[loadOrder[1]], dataIndices[loadOrder[2]]);
      }
      const HistogramData::Counts counts(spectrum.begin(), spectrum.end());
      const HistogramData::CountVariances countVariances(spectrum.begin(), spectrum.end());
      if (pointData) {
        ws->setCounts(currentSpectrum, counts);
        ws->setCountVariances(currentSpectrum, countVariances);
        ws->setPoints(currentSpectrum, histoPoints);
      } else {
        ws->setHistogram(currentSpectrum, binEdges, counts);
      }
      const specnum_t detectorID = customDetectorIDs ? detectorIDs[currentSpectrum] : currentSpectrum;
      ws->getSpectrum(currentSpectrum).setSpectrumNo(detectorID);
    }
  }
}
 
void LoadHelper::loadingOrder(const std::tuple<short, short, short> &dataOrder, int *dataIndices) {
  if (std::get<0>(dataOrder) != 0)
    dataIndices[0] = std::get<1>(dataOrder) == 0 ? 1 : 2;
  if (std::get<1>(dataOrder) != 1)
    dataIndices[1] = std::get<0>(dataOrder) == 1 ? 0 : 2;
  if (std::get<2>(dataOrder) != 2)
    dataIndices[2] = std::get<1>(dataOrder) == 2 ? 1 : 0;
}
 
void LoadHelper::fillMovingWorkspace(const API::MatrixWorkspace_sptr &ws, const Mantid::Nexus::NXInt &data,
                                     const std::vector<double> &xAxis, int64_t initialSpectrum,
                                     const std::set<detid_t> &acceptedDetectorIDs,
                                     const std::vector<detid_t> &customDetectorIDs,
                                     const std::tuple<short, short, short> &axisOrder) {
 
  const auto useCustomSpectraMap = customDetectorIDs.size() != 0;
  const auto useAcceptedDetectorIDs = acceptedDetectorIDs.size() != 0;
 
  std::array<Nexus::dimsize_t, 3U> dims = {data.dim0(), data.dim1(), data.dim2()};
  const auto nTubes = dims[std::get<0>(axisOrder)];
  const auto nPixels = dims[std::get<1>(axisOrder)];
  const auto nScans = dims[std::get<2>(axisOrder)];
 
  int nSkipped = 0;
#pragma omp parallel for if (Kernel::threadSafe(*ws))
  for (specnum_t tube_no = 0; tube_no < static_cast<specnum_t>(nTubes); ++tube_no) {
    for (specnum_t pixel_no = 0; pixel_no < static_cast<specnum_t>(nPixels); ++pixel_no) {
      specnum_t currentDetector =
          static_cast<specnum_t>(initialSpectrum) + tube_no * static_cast<specnum_t>(nPixels) + pixel_no;
      if (useAcceptedDetectorIDs && std::find(acceptedDetectorIDs.cbegin(), acceptedDetectorIDs.cend(),
                                              currentDetector) == acceptedDetectorIDs.end()) {
        nSkipped++;
        continue;
      }
      currentDetector -= nSkipped;
      int64_t currentSpectrum;
      if (useCustomSpectraMap)
        currentSpectrum = customDetectorIDs[currentDetector - initialSpectrum];
      else
        currentSpectrum = currentDetector;
      currentSpectrum *= nScans;
      for (specnum_t channel_no = 0; channel_no < static_cast<specnum_t>(nScans); ++channel_no) {
        auto spectrumValue = data(channel_no, tube_no, pixel_no);
        ws->mutableY(currentSpectrum) = spectrumValue;
        ws->mutableE(currentSpectrum) = sqrt(spectrumValue);
        ws->mutableX(currentSpectrum) = xAxis;
        currentSpectrum++;
      }
    }
  }
}
 
void LoadHelper::replaceZeroErrors(const API::MatrixWorkspace_sptr &ws, double zeroCountsError) {
  for (size_t spectrum_no = 0; spectrum_no < ws->getNumberHistograms(); ++spectrum_no) {
    auto &errorAxis = ws->mutableE(spectrum_no);
    std::transform(errorAxis.begin(), errorAxis.end(), errorAxis.begin(),
                   [zeroCountsError](const auto &error) { return error == 0 ? zeroCountsError : error; });
  }
}
 
Nexus::NXInt LoadHelper::getIntDataset(const Nexus::NXEntry &entry, const std::string &groupName) {
  auto dataGroup = entry.openNXData(groupName);
  return dataGroup.openIntData();
}
 
Nexus::NXDouble LoadHelper::getDoubleDataset(const Nexus::NXEntry &entry, const std::string &groupName) {
  auto dataGroup = entry.openNXData(groupName);
  return dataGroup.openDoubleData();
}
 
} // namespace DataHandling
} // namespace Mantid