UnitsConversionHelper.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 "MantidMDAlgorithms/UnitsConversionHelper.h"
#include "MantidAPI/NumericAxis.h"
#include "MantidKernel/Strings.h"
#include "MantidKernel/UnitFactory.h"
#include <cmath>
 
using Mantid::Kernel::UnitParams;
 
namespace Mantid::MDAlgorithms {
 
CnvrtToMD::ConvertUnits UnitsConversionHelper::analyzeUnitsConversion(const std::string &UnitsFrom,
                                                                      const std::string &UnitsTo, bool forceViaTOF) {
  // if units are equal, no conversion is necessary;
  if (UnitsFrom == UnitsTo)
    return CnvrtToMD::ConvertNo;
 
  // get all known units:
  std::vector<std::string> AllKnownUnits = Kernel::UnitFactory::Instance().getKeys();
 
  // check if unit conversion is possible at all:
  if (Kernel::Strings::isMember(AllKnownUnits, UnitsFrom) < 0)
    throw(std::invalid_argument(" Can not initate conversion from unknown unit: " + UnitsFrom));
 
  if (Kernel::Strings::isMember(AllKnownUnits, UnitsFrom) < 0)
    throw(std::invalid_argument(" Can not initiate conversion to unknown unit: " + UnitsTo));
 
  // is a quick conversion available?
  m_SourceWSUnit = Kernel::UnitFactory::Instance().create(UnitsFrom);
  if (m_SourceWSUnit->quickConversion(UnitsTo, m_Factor, m_Power) && !forceViaTOF) {
    return CnvrtToMD::ConvertFast;
  } else {
    // are the input units TOF?
    if (UnitsFrom == "TOF") {
      return CnvrtToMD::ConvertFromTOF;
    } else { // convert using TOF
      m_TargetUnit = Kernel::UnitFactory::Instance().create(UnitsTo);
      return CnvrtToMD::ConvertByTOF;
    }
  }
}
bool UnitsConversionHelper::isUnitConverted() const { return m_UnitCnvrsn != CnvrtToMD::ConvertNo; }
void UnitsConversionHelper::initialize(const MDWSDescription &targetWSDescr, const std::string &unitsTo,
                                       bool forceViaTOF) {
  // obtain input workspace units
  API::MatrixWorkspace_const_sptr inWS2D = targetWSDescr.getInWS();
  if (!inWS2D)
    throw(std::runtime_error("UnitsConversionHelper::initialize Should not be "
                             "able to call this function when workpsace is "
                             "undefined"));
 
  API::NumericAxis *pAxis = dynamic_cast<API::NumericAxis *>(inWS2D->getAxis(0));
  if (!pAxis)
    throw(std::invalid_argument("Cannot retrieve numeric X axis from the input workspace: " + inWS2D->getName()));
 
  std::string unitsFrom = inWS2D->getAxis(0)->unit()->unitID();
 
  // get detectors positions and other data needed for units conversion:
  if (!(targetWSDescr.m_PreprDetTable))
    throw std::runtime_error("MDWSDescription does not have a detector table");
 
  auto Emode = static_cast<int>(targetWSDescr.getEMode());
 
  this->initialize(unitsFrom, unitsTo, targetWSDescr.m_PreprDetTable, Emode, forceViaTOF);
}
// the helper function which used in the code below to simplify check if the
// variable is in range
inline bool inRange(const std::pair<double, double> &range, const double &val) {
  return val >= range.first && val <= range.second;
}
// the helper function which used in the code below to simplify check if the
// variable is in range
inline bool inRange(const double &xMin, const double &xMax, const double &val) { return val >= xMin && val <= xMax; }
 
std::pair<double, double> UnitsConversionHelper::getConversionRange(double x1, double x2) const {
  std::pair<double, double> range;
  range.first = std::min(x1, x2);
  range.second = std::max(x1, x2);
 
  switch (m_UnitCnvrsn) {
  case (CnvrtToMD::ConvertNo): {
    return range;
  }
  case (CnvrtToMD::ConvertFast): {
    auto trRange = m_TargetUnit->conversionRange();
    double u1 = this->convertUnits(x1);
    double u2 = this->convertUnits(x2);
 
    if (!inRange(trRange, u1) || !inRange(trRange, u2)) // hopefully it is a rare event
    {
      double uMin = std::min(u1, u2);
      double uMax = std::max(u1, u2);
      if (inRange(uMin, uMax, trRange.first)) {
        double t1 = m_TargetUnit->singleToTOF(trRange.first);
        range.first = m_TargetUnit->singleFromTOF(t1);
      }
      if (inRange(uMin, uMax, trRange.second)) {
        double t2 = m_TargetUnit->singleToTOF(trRange.second);
        range.second = m_SourceWSUnit->singleFromTOF(t2);
      }
    }
    return range;
  }
  case (CnvrtToMD::ConvertFromTOF): {
    double tMin = m_TargetUnit->conversionTOFMin();
    double tMax = m_TargetUnit->conversionTOFMax();
 
    if (inRange(tMin, tMax, x1) && inRange(tMin, tMax, x2)) {
      return range;
    } else {
      if (inRange(range, tMin))
        range.first = tMin;
      if (inRange(range, tMax))
        range.second = tMax;
    }
    return range;
  }
  case (CnvrtToMD::ConvertByTOF): {
    auto source_range = m_SourceWSUnit->conversionRange();
    if (!inRange(source_range, range.first) || !inRange(source_range, range.second)) {
      if (inRange(range, source_range.first))
        range.first = source_range.first;
      if (inRange(range, source_range.second))
        range.second = source_range.second;
    }
    double tof1 = m_SourceWSUnit->singleToTOF(range.first);
    double tof2 = m_SourceWSUnit->singleToTOF(range.second);
    if (std::isnan(tof1) || std::isnan(tof2)) {
      if (range.first < source_range.first)
        range.first = source_range.first;
      if (range.second > source_range.second)
        range.second = source_range.second;
      tof1 = m_SourceWSUnit->singleToTOF(range.first);
      tof2 = m_SourceWSUnit->singleToTOF(range.second);
    }
 
    double tMin = m_TargetUnit->conversionTOFMin();
    double tMax = m_TargetUnit->conversionTOFMax();
    if (inRange(tMin, tMax, tof1) && inRange(tMin, tMax, tof2)) {
      return range;
    } else {
      double u1 = range.first;
      double u2 = range.second;
      if (inRange(tof1, tof2, tMin))
        u1 = m_SourceWSUnit->singleFromTOF(tMin);
      if (inRange(tof1, tof2, tMax))
        u2 = m_SourceWSUnit->singleFromTOF(tMax);
      range.first = std::min(u1, u2);
      range.second = std::max(u1, u2);
    }
    return range;
  }
  default:
    throw std::runtime_error("updateConversion: unknown type of conversion requested");
  }
}
 
void UnitsConversionHelper::initialize(const std::string &unitsFrom, const std::string &unitsTo,
                                       const DataObjects::TableWorkspace_const_sptr &DetWS, int Emode,
                                       bool forceViaTOF) {
  m_Emode = Emode;
 
  if (!DetWS)
    throw std::runtime_error("UnitsConversionHelper::initialize called with "
                             "empty preprocessed detectors table");
 
  // Check how the source units relate to the units requested and create source
  // units
  m_UnitCnvrsn = analyzeUnitsConversion(unitsFrom, unitsTo, forceViaTOF);
 
  // create target units class
  m_TargetUnit = Kernel::UnitFactory::Instance().create(unitsTo);
  if (!m_TargetUnit)
    throw(std::runtime_error(" Cannot retrieve target unit from the units factory"));
 
  // get access to all values used by unit conversion.
  m_pTwoThetas = &(DetWS->getColVector<double>("TwoTheta"));
  m_pL2s = &(DetWS->getColVector<double>("L2"));
 
  m_L1 = DetWS->getLogs()->getPropertyValueAsType<double>("L1");
 
  // get efix
  m_Efix = DetWS->getLogs()->getPropertyValueAsType<double>("Ei");
  m_pEfixedArray = nullptr;
  if (m_Emode == static_cast<int>(Kernel::DeltaEMode::Indirect))
    m_pEfixedArray = DetWS->getColDataArray<float>("eFixed");
 
  m_pDIFAs = &(DetWS->getColVector<double>("DIFA"));
  m_pDIFCs = &(DetWS->getColVector<double>("DIFC"));
  m_pTZEROs = &(DetWS->getColVector<double>("TZERO"));
 
  // set up conversion to working state -- in some tests it can be used straight
  // from the beginning.
  m_TwoTheta = (*m_pTwoThetas)[0];
  m_L2 = (*m_pL2s)[0];
  double Efix = m_Efix;
  if (m_pEfixedArray)
    Efix = static_cast<double>(*(m_pEfixedArray + 0));
  m_DIFA = (*m_pDIFAs)[0];
  m_DIFC = (*m_pDIFCs)[0];
  m_TZERO = (*m_pTZEROs)[0];
 
  m_TargetUnit->initialize(m_L1, m_Emode,
                           {{UnitParams::l2, m_L2},
                            {UnitParams::twoTheta, m_TwoTheta},
                            {UnitParams::efixed, Efix},
                            {UnitParams::difa, m_DIFA},
                            {UnitParams::difc, m_DIFC},
                            {UnitParams::tzero, m_TZERO}});
  if (m_SourceWSUnit) {
    m_SourceWSUnit->initialize(m_L1, m_Emode,
                               {{UnitParams::l2, m_L2},
                                {UnitParams::twoTheta, m_TwoTheta},
                                {UnitParams::efixed, Efix},
                                {UnitParams::difa, m_DIFA},
                                {UnitParams::difc, m_DIFC},
                                {UnitParams::tzero, m_TZERO}});
  }
}
void UnitsConversionHelper::updateConversion(size_t i) {
  switch (m_UnitCnvrsn) {
  case (CnvrtToMD::ConvertNo):
    return;
  case (CnvrtToMD::ConvertFast):
    return;
  case (CnvrtToMD::ConvertFromTOF): {
    m_TwoTheta = (*m_pTwoThetas)[i];
    m_L2 = (*m_pL2s)[i];
    double Efix = m_Efix;
    if (m_pEfixedArray)
      Efix = static_cast<double>(*(m_pEfixedArray + i));
    m_DIFA = (*m_pDIFAs)[i];
    m_DIFC = (*m_pDIFCs)[i];
    m_TZERO = (*m_pTZEROs)[i];
 
    m_TargetUnit->initialize(m_L1, m_Emode,
                             {{UnitParams::l2, m_L2},
                              {UnitParams::twoTheta, m_TwoTheta},
                              {UnitParams::efixed, Efix},
                              {UnitParams::difa, m_DIFA},
                              {UnitParams::difc, m_DIFC},
                              {UnitParams::tzero, m_TZERO}});
    return;
  }
  case (CnvrtToMD::ConvertByTOF): {
    m_TwoTheta = (*m_pTwoThetas)[i];
    m_L2 = (*m_pL2s)[i];
    double Efix = m_Efix;
    if (m_pEfixedArray)
      Efix = static_cast<double>(*(m_pEfixedArray + i));
    m_DIFA = (*m_pDIFAs)[i];
    m_DIFC = (*m_pDIFCs)[i];
    m_TZERO = (*m_pTZEROs)[i];
 
    Kernel::UnitParametersMap pmap = {{UnitParams::l2, m_L2},     {UnitParams::twoTheta, m_TwoTheta},
                                      {UnitParams::efixed, Efix}, {UnitParams::difa, m_DIFA},
                                      {UnitParams::difc, m_DIFC}, {UnitParams::tzero, m_TZERO}};
    m_TargetUnit->initialize(m_L1, m_Emode, pmap);
    m_SourceWSUnit->initialize(m_L1, m_Emode, pmap);
    return;
  }
  default:
    throw std::runtime_error("updateConversion: unknown type of conversion requested");
  }
} // namespace MDAlgorithms
double UnitsConversionHelper::convertUnits(double val) const {
  switch (m_UnitCnvrsn) {
  case (CnvrtToMD::ConvertNo): {
    return val;
  }
  case (CnvrtToMD::ConvertFast): {
    return m_Factor * std::pow(val, m_Power);
  }
  case (CnvrtToMD::ConvertFromTOF): {
    return m_TargetUnit->singleFromTOF(val);
  }
  case (CnvrtToMD::ConvertByTOF): {
    double tof = m_SourceWSUnit->singleToTOF(val);
    return m_TargetUnit->singleFromTOF(tof);
  }
  default:
    throw std::runtime_error("updateConversion: unknown type of conversion requested");
  }
}
// copy constructor;
UnitsConversionHelper::UnitsConversionHelper(const UnitsConversionHelper &another) {
  m_UnitCnvrsn = another.m_UnitCnvrsn;
  m_Factor = another.m_Factor;
  m_Power = another.m_Power;
 
  m_Emode = another.m_Emode;
  m_L1 = another.m_L1;
  m_Efix = another.m_Efix;
  m_TwoTheta = another.m_TwoTheta;
  m_L2 = another.m_L2;
  m_DIFA = another.m_DIFA;
  m_DIFC = another.m_DIFC;
  m_TZERO = another.m_TZERO;
  m_pTwoThetas = another.m_pTwoThetas;
  m_pL2s = another.m_pL2s;
  m_pEfixedArray = another.m_pEfixedArray;
  m_pDIFAs = another.m_pDIFAs;
  m_pDIFCs = another.m_pDIFCs;
  m_pTZEROs = another.m_pTZEROs;
 
  if (another.m_SourceWSUnit)
    m_SourceWSUnit = Kernel::Unit_sptr(another.m_SourceWSUnit->clone());
  if (another.m_TargetUnit)
    m_TargetUnit = Kernel::Unit_sptr(another.m_TargetUnit->clone());
}
 
UnitsConversionHelper::UnitsConversionHelper()
    : m_UnitCnvrsn(CnvrtToMD::ConvertNo), m_Factor(1), m_Power(1), m_Emode(-1), // undefined
      m_L1(1), m_Efix(1), m_TwoTheta(0), m_L2(1), m_DIFA(0.), m_DIFC(0.), m_TZERO(0.), m_pTwoThetas(nullptr),
      m_pL2s(nullptr), m_pEfixedArray(nullptr), m_pDIFAs(nullptr), m_pDIFCs(nullptr), m_pTZEROs(nullptr) {}
 
} // namespace Mantid::MDAlgorithms