LoadHelper.cpp

/*
 * Helper file to gather common routines to the Loaders
 * */

#include "MantidDataHandling/LoadHelper.h"

#include "MantidAPI/MatrixWorkspace.h"
#include "MantidGeometry/Instrument/ComponentHelper.h"

#include <nexus/napi.h>

#include <boost/algorithm/string/predicate.hpp> //assert(boost::algorithm::ends_with("mystring", "ing"));

namespace Mantid {
namespace DataHandling {

namespace {
/// static logger
Kernel::Logger g_log("LoadHelper");
}

using namespace Kernel;
using namespace API;

/**
 * Finds the path for the instrument name in the nexus file
 * Usually of the form: entry0/\<NXinstrument class\>/name
 */
std::string
LoadHelper::findInstrumentNexusPath(const Mantid::NeXus::NXEntry &firstEntry) {
  std::string insNamePath;
  std::vector<Mantid::NeXus::NXClassInfo> v = firstEntry.groups();
  for (auto it = v.begin(); it < v.end(); it++) {
    if (it->nxclass == "NXinstrument") {
      insNamePath = it->nxname;
      break;
    }
  }
  return insNamePath;
}

std::string
LoadHelper::getStringFromNexusPath(const Mantid::NeXus::NXEntry &firstEntry,
                                   const std::string &nexusPath) {
  return firstEntry.getString(nexusPath);
}

double
LoadHelper::getDoubleFromNexusPath(const Mantid::NeXus::NXEntry &firstEntry,
                                   const std::string &nexusPath) {
  return firstEntry.getFloat(nexusPath);
}

/**
 * Gets the time binning from a Nexus float array
 * Adds an extra bin at the end
 */
std::vector<double> LoadHelper::getTimeBinningFromNexusPath(
    const Mantid::NeXus::NXEntry &firstEntry, const std::string &nexusPath) {

  Mantid::NeXus::NXFloat timeBinningNexus = firstEntry.openNXFloat(nexusPath);
  timeBinningNexus.load();

  size_t numberOfBins =
      static_cast<size_t>(timeBinningNexus.dim0()) + 1; // boundaries

  float *timeBinning_p = &timeBinningNexus[0];
  std::vector<double> timeBinning(numberOfBins);
  timeBinning.assign(timeBinning_p, timeBinning_p + numberOfBins);
  // calculate the extra bin at the end
  timeBinning[numberOfBins - 1] =
      timeBinning[numberOfBins - 2] + timeBinning[1] - timeBinning[0];

  return timeBinning;
}
/**
 * Calculate Neutron Energy from wavelength: \f$ E = h^2 / 2m\lambda ^2 \f$
 *  @param wavelength :: wavelength in \f$ \mbox{\AA} \f$
 *  @return tof in seconds
 */
double LoadHelper::calculateEnergy(double wavelength) {
  double e =
      (PhysicalConstants::h * PhysicalConstants::h) /
      (2 * PhysicalConstants::NeutronMass * wavelength * wavelength * 1e-20) /
      PhysicalConstants::meV;
  return e;
}

/**
 * Calculate TOF from distance
 *  @param distance :: distance in meters
 *  @param wavelength :: wavelength to calculate TOF from
 *  @return tof in seconds
 */
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;
}

double LoadHelper::getL1(const API::MatrixWorkspace_sptr &workspace) {
  Geometry::Instrument_const_sptr instrument = workspace->getInstrument();
  Geometry::IComponent_const_sptr sample = instrument->getSample();
  double l1 = instrument->getSource()->getDistance(*sample);
  return l1;
}

double LoadHelper::getL2(const API::MatrixWorkspace_sptr &workspace,
                         int detId) {
  // Get a pointer to the instrument contained in the workspace
  Geometry::Instrument_const_sptr instrument = workspace->getInstrument();
  // Get the distance between the source and the sample (assume in metres)
  Geometry::IComponent_const_sptr sample = instrument->getSample();
  // Get the sample-detector distance for this detector (in metres)
  double l2 =
      workspace->getDetector(detId)->getPos().distance(sample->getPos());
  return l2;
}

/*
 * Get instrument property as double
 * @s - input property name
 *
 */
double
LoadHelper::getInstrumentProperty(const API::MatrixWorkspace_sptr &workspace,
                                  std::string s) {
  std::vector<std::string> prop =
      workspace->getInstrument()->getStringParameter(s);
  if (prop.empty()) {
    g_log.debug("Property <" + s + "> doesn't exist!");
    return EMPTY_DBL();
  } else {
    g_log.debug() << "Property <" + s + "> = " << prop[0] << '\n';
    return boost::lexical_cast<double>(prop[0]);
  }
}

/**
 * Add properties from a nexus file to
 * the workspace run.
 * API entry for recursive routine below
 *
 *
 * @param nxfileID    :: Nexus file handle to be parsed, just after an
 *NXopengroup
 * @param runDetails  :: where to add properties
 *
 */
void LoadHelper::addNexusFieldsToWsRun(NXhandle nxfileID,
                                       API::Run &runDetails) {
  std::string emptyStr; // needed for first call
  int datatype;
  char nxname[NX_MAXNAMELEN], nxclass[NX_MAXNAMELEN];

  // As a workaround against some "not so good" old ILL nexus files
  // (ILLIN5_Vana_095893.nxs for example)
  // we begin the parse on the first entry (entry0). This allow to avoid the
  // bogus entries that follows.

  NXstatus getnextentry_status =
      NXgetnextentry(nxfileID, nxname, nxclass, &datatype);
  if (getnextentry_status == NX_OK) {
    NXstatus opengroup_status;
    if ((opengroup_status = NXopengroup(nxfileID, nxname, nxclass)) == NX_OK) {
      if (std::string(nxname) == "entry0") {
        recurseAndAddNexusFieldsToWsRun(nxfileID, runDetails, emptyStr,
                                        emptyStr, 1 /* level */);
      } else {
        g_log.debug() << "Unexpected group name in nexus file : " << nxname
                      << '\n';
      }
      NXclosegroup(nxfileID);
    }
  }
}

/**
 * Recursively add properties from a nexus file to
 * the workspace run.
 *
 * @param nxfileID    :: Nexus file handle to be parsed, just after an
 *NXopengroup
 * @param runDetails  :: where to add properties
 * @param parent_name :: nexus caller name
 * @param parent_class :: nexus caller class
 * @param level       :: current level in nexus tree
 *
 */
void LoadHelper::recurseAndAddNexusFieldsToWsRun(NXhandle nxfileID,
                                                 API::Run &runDetails,
                                                 std::string &parent_name,
                                                 std::string &parent_class,
                                                 int level) {

  std::string indent_str(level * 2, ' '); // Two space by indent level

  // Link ?

  // Attributes ?

  // Classes
  NXstatus getnextentry_status; ///< return status
  int datatype; ///< NX data type if a dataset, e.g. NX_CHAR, NX_FLOAT32, see
  /// napi.h
  char nxname[NX_MAXNAMELEN], nxclass[NX_MAXNAMELEN];
  nxname[0] = '0';
  nxclass[0] = '0';

  bool has_entry = true; // follows getnextentry_status
  while (has_entry) {
    getnextentry_status = NXgetnextentry(nxfileID, nxname, nxclass, &datatype);

    if (getnextentry_status == NX_OK) {
      g_log.debug() << indent_str << parent_name << "." << nxname << " ; "
                    << nxclass << '\n';

      NXstatus opengroup_status;
      NXstatus opendata_status;

      if ((opengroup_status = NXopengroup(nxfileID, nxname, nxclass)) ==
          NX_OK) {

        // Go down to one level
        std::string p_nxname(
            nxname); // current names can be useful for next level
        std::string p_nxclass(nxclass);

        recurseAndAddNexusFieldsToWsRun(nxfileID, runDetails, p_nxname,
                                        p_nxclass, level + 1);

        NXclosegroup(nxfileID);
      } // if(NXopengroup
      else if ((opendata_status = NXopendata(nxfileID, nxname)) == NX_OK) {
        // dump_attributes(nxfileID, indent_str);
        g_log.debug() << indent_str << nxname << " opened.\n";

        if (parent_class == "NXData" || parent_class == "NXMonitor" ||
            std::string(nxname) == "data") {
          g_log.debug() << indent_str << "skipping " << parent_class << " ("
                        << nxname << ")\n";
          /* nothing */
        } else { // create a property
          int rank;
          int dims[4];
          int type;
          dims[0] = dims[1] = dims[2] = dims[3] = 0;

          std::string property_name =
              (parent_name.empty() ? nxname : parent_name + "." + nxname);

          g_log.debug() << indent_str << "considering property "
                        << property_name << '\n';

          // Get the value
          NXgetinfo(nxfileID, &rank, dims, &type);

          // Note, we choose to only build properties on small float arrays
          // filter logic is below
          bool build_small_float_array = false; // default

          if ((type == NX_FLOAT32) || (type == NX_FLOAT64)) {
            if ((rank == 1) && (dims[0] <= 9)) {
              build_small_float_array = true;
            } else {
              g_log.debug() << indent_str
                            << "ignored multi dimension float data on "
                            << property_name << '\n';
            }
          } else if (type != NX_CHAR) {
            if ((rank != 1) || (dims[0] != 1) || (dims[1] != 1) ||
                (dims[2] != 1) || (dims[3] != 1)) {
              g_log.debug() << indent_str << "ignored multi dimension data on "
                            << property_name << '\n';
            }
          }

          void *dataBuffer;
          NXmalloc(&dataBuffer, rank, dims, type);

          if (NXgetdata(nxfileID, dataBuffer) != NX_OK) {
            NXfree(&dataBuffer);
            throw std::runtime_error("Cannot read data from NeXus file");
          }

          if (type == NX_CHAR) {
            std::string property_value(
                reinterpret_cast<const char *>(dataBuffer));
            if (boost::algorithm::ends_with(property_name, "_time")) {
              // That's a time value! Convert to Mantid standard
              property_value = dateTimeInIsoFormat(property_value);
            }
            runDetails.addProperty(property_name, property_value);

          } else if ((type == NX_FLOAT32) || (type == NX_FLOAT64) ||
                     (type == NX_INT16) || (type == NX_INT32) ||
                     (type == NX_UINT16)) {

            // Look for "units"
            NXstatus units_status;
            char units_sbuf[NX_MAXNAMELEN];
            int units_len = NX_MAXNAMELEN;
            int units_type = NX_CHAR;

            units_status = NXgetattr(nxfileID, "units", units_sbuf, &units_len,
                                     &units_type);
            if (units_status != NX_ERROR) {
              g_log.debug() << indent_str << "[ " << property_name
                            << " has unit " << units_sbuf << " ]\n";
            }

            if ((type == NX_FLOAT32) || (type == NX_FLOAT64)) {
              // Mantid numerical properties are double only.
              double property_double_value = 0.0;

              // Simple case, one value
              if (dims[0] == 1) {
                if (type == NX_FLOAT32) {
                  property_double_value =
                      *(reinterpret_cast<float *>(dataBuffer));
                } else if (type == NX_FLOAT64) {
                  property_double_value =
                      *(reinterpret_cast<double *>(dataBuffer));
                }
                if (units_status != NX_ERROR)
                  runDetails.addProperty(property_name, property_double_value,
                                         std::string(units_sbuf));
                else
                  runDetails.addProperty(property_name, property_double_value);
              } else if (build_small_float_array) {
                // An array, converted to "name_index", with index < 10 (see
                // test above)
                for (int dim_index = 0; dim_index < dims[0]; dim_index++) {
                  if (type == NX_FLOAT32) {
                    property_double_value =
                        (reinterpret_cast<float *>(dataBuffer))[dim_index];
                  } else if (type == NX_FLOAT64) {
                    property_double_value =
                        (reinterpret_cast<double *>(dataBuffer))[dim_index];
                  }
                  std::string indexed_property_name = property_name +
                                                      std::string("_") +
                                                      std::to_string(dim_index);
                  if (units_status != NX_ERROR)
                    runDetails.addProperty(indexed_property_name,
                                           property_double_value,
                                           std::string(units_sbuf));
                  else
                    runDetails.addProperty(indexed_property_name,
                                           property_double_value);
                }
              }

            } else {
              // int case
              int property_int_value = 0;
              if (type == NX_INT16) {
                property_int_value =
                    *(reinterpret_cast<short int *>(dataBuffer));
              } else if (type == NX_INT32) {
                property_int_value = *(reinterpret_cast<int *>(dataBuffer));
              } else if (type == NX_UINT16) {
                property_int_value =
                    *(reinterpret_cast<short unsigned int *>(dataBuffer));
              }

              if (units_status != NX_ERROR)
                runDetails.addProperty(property_name, property_int_value,
                                       std::string(units_sbuf));
              else
                runDetails.addProperty(property_name, property_int_value);

            } // if (type==...

          } else {
            g_log.debug() << indent_str << "unexpected data on "
                          << property_name << '\n';
          } // test on nxdata type

          NXfree(&dataBuffer);
          dataBuffer = nullptr;

        } // if (parent_class == "NXData" || parent_class == "NXMonitor") else

        NXclosedata(nxfileID);
      } else {
        g_log.debug() << indent_str << "unexpected status (" << opendata_status
                      << ") on " << nxname << '\n';
      }

    } else if (getnextentry_status == NX_EOD) {
      g_log.debug() << indent_str << "End of Dir\n";
      has_entry = false; // end of loop
    } else {
      g_log.debug() << indent_str << "unexpected status ("
                    << getnextentry_status << ")\n";
      has_entry = false; // end of loop
    }

  } // while

} // recurseAndAddNexusFieldsToWsRun

/**
 * Show attributes attached to the current Nexus entry
 *
 * @param nxfileID The Nexus entry
 * @param indentStr Indent spaces do display nexus entries as a tree
 *
 */
void LoadHelper::dumpNexusAttributes(NXhandle nxfileID,
                                     std::string &indentStr) {
  // Attributes
  NXname pName;
  int iLength, iType;
  int rank;
  int dims[4];
  int nbuff = 127;
  boost::shared_array<char> buff(new char[nbuff + 1]);

  while (NXgetnextattra(nxfileID, pName, &rank, dims, &iType) != NX_EOD) {
    g_log.debug() << indentStr << '@' << pName << " = ";
    if (rank > 1) { // mantid only supports single value attributes
      throw std::runtime_error(
          "Encountered attribute with multi-dimensional array value");
    }
    iLength = dims[0]; // to clarify things
    if (iType != NX_CHAR && iLength != 1) {
      throw std::runtime_error("Encountered attribute with array value");
    }

    switch (iType) {
    case NX_CHAR: {
      if (iLength > nbuff + 1) {
        nbuff = iLength;
        buff.reset(new char[nbuff + 1]);
      }
      int nz = iLength + 1;
      NXgetattr(nxfileID, pName, buff.get(), &nz, &iType);
      g_log.debug() << indentStr << buff.get() << '\n';
      break;
    }
    case NX_INT16: {
      short int value;
      NXgetattr(nxfileID, pName, &value, &iLength, &iType);
      g_log.debug() << indentStr << value << '\n';
      break;
    }
    case NX_INT32: {
      int value;
      NXgetattr(nxfileID, pName, &value, &iLength, &iType);
      g_log.debug() << indentStr << value << '\n';
      break;
    }
    case NX_UINT16: {
      short unsigned int value;
      NXgetattr(nxfileID, pName, &value, &iLength, &iType);
      g_log.debug() << indentStr << value << '\n';
      break;
    }
    } // switch
  }   // while
}
/**
 * Parses the date as formatted at the ILL:
 * 29-Jun-12 11:27:26
 * and converts it to the ISO format used in Mantid:
 * ISO8601 format string: "yyyy-mm-ddThh:mm:ss[Z+-]tz:tz"
 *
 *  @param dateToParse :: date as string
 *  @return date as required in Mantid
 */
std::string LoadHelper::dateTimeInIsoFormat(std::string dateToParse) {
  namespace bt = boost::posix_time;
  // 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(API::MatrixWorkspace_sptr ws,
                               const std::string &componentName,
                               const V3D &newPos) {

  try {

    Geometry::Instrument_const_sptr instrument = ws->getInstrument();
    Geometry::IComponent_const_sptr component =
        instrument->getComponentByName(componentName);

    // g_log.debug() << tube->getName() << " : t = " << theta << " ==> t = " <<
    // newTheta << "\n";
    Geometry::ParameterMap &pmap = ws->instrumentParameters();
    Geometry::ComponentHelper::moveComponent(
        *component, pmap, newPos, Geometry::ComponentHelper::Absolute);

  } catch (Mantid::Kernel::Exception::NotFoundError &) {
    throw std::runtime_error("Error when trying to move the " + componentName +
                             " : NotFoundError");
  } catch (std::runtime_error &) {
    throw std::runtime_error("Error when trying to move the " + componentName +
                             " : runtime_error");
  }
}

void LoadHelper::rotateComponent(API::MatrixWorkspace_sptr ws,
                                 const std::string &componentName,
                                 const Kernel::Quat &rot) {

  try {

    Geometry::Instrument_const_sptr instrument = ws->getInstrument();
    Geometry::IComponent_const_sptr component =
        instrument->getComponentByName(componentName);

    // g_log.debug() << tube->getName() << " : t = " << theta << " ==> t = " <<
    // newTheta << "\n";
    Geometry::ParameterMap &pmap = ws->instrumentParameters();
    Geometry::ComponentHelper::rotateComponent(
        *component, pmap, rot, Geometry::ComponentHelper::Absolute);

  } catch (Mantid::Kernel::Exception::NotFoundError &) {
    throw std::runtime_error("Error when trying to move the " + componentName +
                             " : NotFoundError");
  } catch (std::runtime_error &) {
    throw std::runtime_error("Error when trying to move the " + componentName +
                             " : runtime_error");
  }
}

V3D LoadHelper::getComponentPosition(API::MatrixWorkspace_sptr ws,
                                     const std::string &componentName) {
  try {
    Geometry::Instrument_const_sptr instrument = ws->getInstrument();
    Geometry::IComponent_const_sptr component =
        instrument->getComponentByName(componentName);
    V3D pos = component->getPos();
    return pos;
  } catch (Mantid::Kernel::Exception::NotFoundError &) {
    throw std::runtime_error("Error when trying to move the " + componentName +
                             " : NotFoundError");
  }
}

template <typename T>
T LoadHelper::getPropertyFromRun(API::MatrixWorkspace_const_sptr inputWS,
                                 const std::string &propertyName) {
  if (inputWS->run().hasProperty(propertyName)) {
    Kernel::Property *prop = inputWS->run().getProperty(propertyName);
    return boost::lexical_cast<T>(prop->value());
  } else {
    std::string mesg =
        "No '" + propertyName + "' property found in the input workspace....";
    throw std::runtime_error(mesg);
  }
}

} // namespace DataHandling
} // namespace Mantid