LoadSassena.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/LoadSassena.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/IAlgorithm.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceGroup.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/Unit.h"
#include "MantidKernel/UnitFactory.h"
 
#include <hdf5_hl.h>
 
#include <utility>
 
namespace Mantid::DataHandling {
 
DECLARE_NEXUS_FILELOADER_ALGORITHM(LoadSassena)
 
 
int LoadSassena::confidence(Kernel::NexusDescriptor &descriptor) const {
  if (descriptor.hasRootAttr("sassena_version") || descriptor.pathExists("/qvectors")) {
    return 99;
  }
  return 0;
}
 
void LoadSassena::registerWorkspace(const API::WorkspaceGroup_sptr &gws, const std::string &wsName,
                                    const DataObjects::Workspace2D_sptr &ws, const std::string &description) {
  UNUSED_ARG(description);
  API::AnalysisDataService::Instance().add(wsName, ws);
  gws->addWorkspace(ws);
}
 
herr_t LoadSassena::dataSetInfo(const hid_t &h5file, const std::string &setName, hsize_t *dims) const {
  H5T_class_t class_id;
  size_t type_size;
  herr_t errorcode = H5LTget_dataset_info(h5file, setName.c_str(), dims, &class_id, &type_size);
  if (errorcode < 0) {
    g_log.error("Unable to read " + setName + " dataset info");
  }
  return errorcode;
}
 
herr_t LoadSassena::dataSetDouble(const hid_t &h5file, const std::string &setName, std::vector<double> &buf) {
  herr_t errorcode = H5LTread_dataset_double(h5file, setName.c_str(), buf.data());
  if (errorcode < 0) {
    this->g_log.error("Cannot read " + setName + " dataset");
  }
  return errorcode;
}
 
/* Helper object and function to sort modulus of Q-vectors
 */
using mypair = std::pair<double, int>;
bool compare(const mypair &left, const mypair &right) { return left.first < right.first; }
HistogramData::Points LoadSassena::loadQvectors(const hid_t &h5file, const API::WorkspaceGroup_sptr &gws,
                                                std::vector<int> &sorting_indexes) {
 
  // store the modulus of the vector
  std::vector<double> qvmod;
 
  const std::string gwsName = this->getPropertyValue("OutputWorkspace");
  const std::string setName("qvectors");
 
  hsize_t dims[3];
  if (dataSetInfo(h5file, setName, dims) < 0) {
    throw Kernel::Exception::FileError("Unable to read " + setName + " dataset info:", m_filename);
  }
  auto nq = static_cast<int>(dims[0]); // number of q-vectors
  std::vector<double> buf(nq * 3);
 
  herr_t errorcode = this->dataSetDouble(h5file, "qvectors", buf);
  if (errorcode < 0) {
    this->g_log.error("LoadSassena::loadQvectors cannot proceed");
    qvmod.resize(0);
    return HistogramData::Points(std::move(qvmod));
  }
 
  qvmod.reserve(nq);
  for (auto curr = buf.cbegin(); curr != buf.cend(); curr += 3) {
    qvmod.emplace_back(sqrt(curr[0] * curr[0] + curr[1] * curr[1] + curr[2] * curr[2]));
  }
 
  if (getProperty("SortByQVectors")) {
    std::vector<mypair> qvmodpair;
    qvmodpair.reserve(nq);
    for (int iq = 0; iq < nq; iq++)
      qvmodpair.emplace_back(qvmod[iq], iq);
    std::sort(qvmodpair.begin(), qvmodpair.end(), compare);
    for (int iq = 0; iq < nq; iq++)
      sorting_indexes.emplace_back(qvmodpair[iq].second);
    std::sort(qvmod.begin(), qvmod.end());
  } else
    for (int iq = 0; iq < nq; iq++)
      sorting_indexes.emplace_back(iq);
 
  DataObjects::Workspace2D_sptr ws = std::dynamic_pointer_cast<DataObjects::Workspace2D>(
      API::WorkspaceFactory::Instance().create("Workspace2D", nq, 3, 3));
  std::string wsName = gwsName + std::string("_") + setName;
  ws->setTitle(wsName);
 
  for (int iq = 0; iq < nq; iq++) {
    auto &Y = ws->mutableY(iq);
    auto curr = std::next(buf.cbegin(), 3 * sorting_indexes[iq]);
    Y.assign(curr, std::next(curr, 3));
  }
 
  ws->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("MomentumTransfer"); // Set the Units
 
  this->registerWorkspace(gws, wsName, ws, "X-axis: origin of Q-vectors; Y-axis: tip of Q-vectors");
  return HistogramData::Points(std::move(qvmod));
}
 
void LoadSassena::loadFQ(const hid_t &h5file, const API::WorkspaceGroup_sptr &gws, const std::string &setName,
                         const HistogramData::Points &qvmod, const std::vector<int> &sorting_indexes) {
 
  auto nq = static_cast<int>(qvmod.size()); // number of q-vectors
  std::vector<double> buf(nq * 2);
  herr_t errorcode = this->dataSetDouble(h5file, setName, buf);
  if (errorcode < 0) {
    this->g_log.error("LoadSassena::loadFQ cannot proceed");
    return;
  }
 
  const std::string gwsName = this->getPropertyValue("OutputWorkspace");
 
  DataObjects::Workspace2D_sptr ws = std::dynamic_pointer_cast<DataObjects::Workspace2D>(
      API::WorkspaceFactory::Instance().create("Workspace2D", 2, nq, nq));
  const std::string wsName = gwsName + std::string("_") + setName;
  ws->setTitle(wsName);
 
  auto &re = ws->mutableY(0); // store the real part
  ws->setPoints(0, qvmod);    // X-axis values are the modulus of the q vector
  auto &im = ws->mutableY(1); // store the imaginary part
  ws->setPoints(1, qvmod);
 
  for (int iq = 0; iq < nq; ++iq) {
    auto curr = std::next(buf.cbegin(), 2 * sorting_indexes[iq]);
    re[iq] = *curr;
    im[iq] = *std::next(curr);
  }
 
  // Set the Units
  ws->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("MomentumTransfer");
 
  this->registerWorkspace(gws, wsName, ws, "X-axis: Q-vector modulus; Y-axis: intermediate structure factor");
}
 
void LoadSassena::loadFQT(const hid_t &h5file, const API::WorkspaceGroup_sptr &gws, const std::string &setName,
                          const HistogramData::Points &qvmod, const std::vector<int> &sorting_indexes) {
 
  hsize_t dims[3];
  if (dataSetInfo(h5file, setName, dims) < 0) {
    this->g_log.error("Unable to read " + setName + " dataset info");
    this->g_log.error("LoadSassena::loadFQT cannot proceed");
    return;
  }
  auto nnt = static_cast<int>(dims[1]); // number of non-negative time points
  int nt = 2 * nnt - 1;                 // number of time points
 
  auto nq = static_cast<int>(qvmod.size()); // number of q-vectors
  std::vector<double> buf(nq * nnt * 2);
  herr_t errorcode = this->dataSetDouble(h5file, setName, buf);
  if (errorcode < 0) {
    this->g_log.error("LoadSassena::loadFQT cannot proceed");
    return;
  }
 
  const std::string gwsName = this->getPropertyValue("OutputWorkspace");
  const double dt = getProperty("TimeUnit"); // time unit increment, in picoseconds;
  DataObjects::Workspace2D_sptr wsRe = std::dynamic_pointer_cast<DataObjects::Workspace2D>(
      API::WorkspaceFactory::Instance().create("Workspace2D", nq, nt, nt));
  const std::string wsReName = gwsName + std::string("_") + setName + std::string(".Re");
  wsRe->setTitle(wsReName);
 
  DataObjects::Workspace2D_sptr wsIm = std::dynamic_pointer_cast<DataObjects::Workspace2D>(
      API::WorkspaceFactory::Instance().create("Workspace2D", nq, nt, nt));
  const std::string wsImName = gwsName + std::string("_") + setName + std::string(".Im");
  wsIm->setTitle(wsImName);
 
  int origin = nnt - 1;
  for (int iq = 0; iq < nq; iq++) {
    auto &reX = wsRe->mutableX(iq);
    auto &imX = wsIm->mutableX(iq);
    auto &reY = wsRe->mutableY(iq);
    auto &imY = wsIm->mutableY(iq);
    const int index = sorting_indexes[iq];
    auto curr = std::next(buf.cbegin(), index * nnt * 2);
    for (int it = 0; it < nnt; it++) {
      reX[origin + it] = it * dt;  // time point for the real part
      reY[origin + it] = *curr;    // real part of the intermediate structure factor
      reX[origin - it] = -it * dt; // symmetric negative time
      reY[origin - it] = *curr;    // symmetric value for the negative time
      ++curr;
      imX[origin + it] = it * dt;
      imY[origin + it] = *curr;
      imX[origin - it] = -it * dt;
      imY[origin - it] = -(*curr); // antisymmetric value for negative times
      ++curr;
    }
  }
 
  // Set the Time unit for the X-axis
  wsRe->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("Label");
  auto unitPtr = std::dynamic_pointer_cast<Kernel::Units::Label>(wsRe->getAxis(0)->unit());
  unitPtr->setLabel("Time", "picoseconds");
 
  wsIm->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create("Label");
  unitPtr = std::dynamic_pointer_cast<Kernel::Units::Label>(wsIm->getAxis(0)->unit());
  unitPtr->setLabel("Time", "picoseconds");
 
  // Create a numeric axis to replace the default vertical one
  auto verticalAxisRe = std::make_unique<API::NumericAxis>(nq);
  auto verticalAxisIm = std::make_unique<API::NumericAxis>(nq);
  auto verticalAxisReRaw = verticalAxisRe.get();
  auto verticalAxisImRaw = verticalAxisIm.get();
  wsRe->replaceAxis(1, std::move(verticalAxisRe));
  wsIm->replaceAxis(1, std::move(verticalAxisIm));
 
  // Now set the axis values
  for (int i = 0; i < nq; ++i) {
    verticalAxisReRaw->setValue(i, qvmod[i]);
    verticalAxisImRaw->setValue(i, qvmod[i]);
  }
 
  // Set the axis units
  verticalAxisReRaw->unit() = Kernel::UnitFactory::Instance().create("MomentumTransfer");
  verticalAxisReRaw->title() = "|Q|";
  verticalAxisImRaw->unit() = Kernel::UnitFactory::Instance().create("MomentumTransfer");
  verticalAxisImRaw->title() = "|Q|";
 
  // Set the X axis title (for conversion to MD)
  wsRe->getAxis(0)->title() = "Energy transfer";
  wsIm->getAxis(0)->title() = "Energy transfer";
 
  // Register the workspaces
  registerWorkspace(gws, wsReName, wsRe, "X-axis: time; Y-axis: real part of intermediate structure factor");
  registerWorkspace(gws, wsImName, wsIm, "X-axis: time; Y-axis: imaginary part of intermediate structure factor");
}
 
void LoadSassena::init() {
  // Declare the Filename algorithm property. Mandatory. Sets the path to the
  // file to load.
  const std::vector<std::string> exts{".h5", ".hd5"};
  declareProperty(std::make_unique<API::FileProperty>("Filename", "", API::FileProperty::Load, exts), "A Sassena file");
  // Declare the OutputWorkspace property
  declareProperty(
      std::make_unique<API::WorkspaceProperty<API::Workspace>>("OutputWorkspace", "", Kernel::Direction::Output),
      "The name of the group workspace to be created.");
  declareProperty(std::make_unique<Kernel::PropertyWithValue<double>>("TimeUnit", 1.0, Kernel::Direction::Input),
                  "The Time unit in between data points, in picoseconds. "
                  "Default is 1.0 picosec.");
  declareProperty(std::make_unique<Kernel::PropertyWithValue<bool>>("SortByQVectors", true, Kernel::Direction::Input),
                  "Sort structure factors by increasing momentum transfer?");
}
 
void LoadSassena::exec() {
  // auto
  // gws=std::dynamic_pointer_cast<API::WorkspaceGroup>(getProperty("OutputWorkspace"));
  // API::WorkspaceGroup_sptr gws=getProperty("OutputWorkspace");
  API::Workspace_sptr ows = getProperty("OutputWorkspace");
 
  API::WorkspaceGroup_sptr gws = std::dynamic_pointer_cast<API::WorkspaceGroup>(ows);
  if (gws && API::AnalysisDataService::Instance().doesExist(gws->getName())) {
    // gws->deepRemoveAll(); // remove workspace members
    API::AnalysisDataService::Instance().deepRemoveGroup(gws->getName());
  } else {
    gws = std::make_shared<API::WorkspaceGroup>();
    setProperty("OutputWorkspace", std::dynamic_pointer_cast<API::Workspace>(gws));
  }
 
  // populate m_validSets
  int nvalidSets = 4;
  const char *validSets[] = {"fq", "fq0", "fq2", "fqt"};
  for (int iSet = 0; iSet < nvalidSets; iSet++)
    this->m_validSets.emplace_back(validSets[iSet]);
 
  // open the HDF5 file for reading
  m_filename = this->getPropertyValue("Filename");
  hid_t h5file = H5Fopen(m_filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
  if (h5file < 0) {
    this->g_log.error("Cannot open " + m_filename);
    throw Kernel::Exception::FileError("Unable to open:", m_filename);
  }
 
  // find out the sassena version used
  char cversion[16];
  if (H5LTget_attribute_string(h5file, "/", "sassena_version", cversion) < 0) {
    this->g_log.error("Unable to read Sassena version");
  }
  // const std::string version(cversion);
  // determine which loader protocol to use based on the version
  // to be done at a later time, maybe implement a Version class
 
  // Block to read the Q-vectors
  std::vector<int> sorting_indexes;
  const auto qvmod = this->loadQvectors(h5file, gws, sorting_indexes);
  if (qvmod.empty()) {
    this->g_log.error("No Q-vectors read. Unable to proceed");
    H5Fclose(h5file);
    return;
  }
 
  // iterate over the valid sets
  std::string setName;
  for (std::vector<std::string>::const_iterator it = this->m_validSets.begin(); it != this->m_validSets.end(); ++it) {
    setName = *it;
    if (H5Lexists(h5file, setName.c_str(), H5P_DEFAULT)) {
      if (setName == "fq" || setName == "fq0" || setName == "fq2")
        this->loadFQ(h5file, gws, setName, qvmod, sorting_indexes);
      else if (setName == "fqt")
        this->loadFQT(h5file, gws, setName, qvmod, sorting_indexes);
    } else
      this->g_log.information("Dataset " + setName + " not present in file");
  } // end of iterate over the valid sets
 
  H5Fclose(h5file);
} // end of LoadSassena::exec()
 
} // namespace Mantid::DataHandling