MDTransfModQ.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/MDTransfModQ.h"
#include "Eigen/Core"
#include "Eigen/Dense"
#include "MantidKernel/RegistrationHelper.h"
#include "MantidMDAlgorithms/DisplayNormalizationSetter.h"
 
namespace Mantid::MDAlgorithms {
// register the class, whith conversion factory under ModQ name
// clang-format off
DECLARE_MD_TRANSFID(MDTransfModQ, |Q|)
// clang-format on
 
 
const std::string MDTransfModQ::inputUnitID(Kernel::DeltaEMode::Type dEmode,
                                            API::MatrixWorkspace_const_sptr inWS) const {
  UNUSED_ARG(inWS);
  switch (dEmode) {
  case (Kernel::DeltaEMode::Elastic):
    return "Momentum";
  case (Kernel::DeltaEMode::Direct):
    return "DeltaE";
  case (Kernel::DeltaEMode::Indirect):
    return "DeltaE";
  default:
    throw(std::invalid_argument(" MDTransfModQ::inputUnitID: this class "
                                "supports only conversion in Elastic and "
                                "Inelastic energy transfer modes"));
  }
}
 
unsigned int MDTransfModQ::getNMatrixDimensions(Kernel::DeltaEMode::Type mode,
                                                API::MatrixWorkspace_const_sptr inWS) const {
  UNUSED_ARG(inWS);
  switch (mode) {
  case (Kernel::DeltaEMode::Direct):
    return 2;
  case (Kernel::DeltaEMode::Indirect):
    return 2;
  case (Kernel::DeltaEMode::Elastic):
    return 1;
  default:
    throw(std::invalid_argument("Unknown or unsupported energy conversion mode"));
  }
}
 
bool MDTransfModQ::calcMatrixCoord(const double &deltaEOrK0, std::vector<coord_t> &Coord, double &signal,
                                   double &ErrSq) const {
  UNUSED_ARG(signal);
  UNUSED_ARG(ErrSq);
  if (m_Emode == Kernel::DeltaEMode::Elastic) {
    return calcMatrixCoordElastic(deltaEOrK0, Coord);
  } else {
    return calcMatrixCoordInelastic(deltaEOrK0, Coord);
  }
}
 
bool MDTransfModQ::calcGenericVariables(std::vector<coord_t> &Coord, size_t nd) {
  // sanity check. If fails, something went fundamentally wrong
  if (m_NMatrixDim + m_AddDimCoordinates.size() != nd) {
    std::string ERR = "Number of matrix dimensions: " + std::to_string(m_NMatrixDim) +
                      " plus number of additional dimensions: " + std::to_string(m_AddDimCoordinates.size()) +
                      " not equal to number of workspace dimensions: " + std::to_string(nd);
    throw(std::invalid_argument(ERR));
  }
 
  // in Elastic case, 1  coordinate (|Q|) came from workspace
  // in inelastic 2 coordinates (|Q| dE) came from workspace. All other are
  // defined by properties.
  // m_NMatrixDim is either 1 in elastic case or 2 in inelastic
  size_t ic(0);
  for (size_t i = m_NMatrixDim; i < nd; i++) {
    if (m_AddDimCoordinates[ic] < m_DimMin[i] || m_AddDimCoordinates[ic] >= m_DimMax[i])
      return false;
    Coord[i] = m_AddDimCoordinates[ic];
    ic++;
  }
  return true;
}
 
bool MDTransfModQ::calcYDepCoordinates(std::vector<coord_t> &Coord, size_t i) {
  UNUSED_ARG(Coord);
  m_ex = (m_DetDirecton + i)->X();
  m_ey = (m_DetDirecton + i)->Y();
  m_ez = (m_DetDirecton + i)->Z();
  // if input energy changes on each detector (efixed, indirect mode only), then
  // set up its value
  if (m_pEfixedArray) {
    m_eFixed = double(*(m_pEfixedArray + i));
    m_kFixed = sqrt(m_eFixed / PhysicalConstants::E_mev_toNeutronWavenumberSq);
  }
  // if spectra masked, this spectra should be excluded
  if (m_pDetMasks) {
    if (*(m_pDetMasks + i) > 0)
      return false;
  }
  return true;
}
 
bool MDTransfModQ::calcMatrixCoordInelastic(const double deltaE, std::vector<coord_t> &Coord) const {
  if (deltaE < m_DimMin[1] || deltaE >= m_DimMax[1])
    return false;
  Coord[1] = static_cast<coord_t>(deltaE);
 
  // x,y,z refer to internal coordinate system where Z is the beam direction
  double qx{0.0}, qy{0.0}, qz{0.0};
  if (this->m_Emode == Kernel::DeltaEMode::Direct) {
    const double kFinal = sqrt((m_eFixed - deltaE) / PhysicalConstants::E_mev_toNeutronWavenumberSq);
    qx = -m_ex * kFinal;
    qy = -m_ey * kFinal;
    qz = m_kFixed - m_ez * kFinal;
  } else {
    const double kInitial = sqrt((m_eFixed + deltaE) / PhysicalConstants::E_mev_toNeutronWavenumberSq);
    qx = -m_ex * m_kFixed;
    qy = -m_ey * m_kFixed;
    qz = kInitial - m_ez * m_kFixed;
  }
  return applyCoordTransf(qx, qy, qz, Coord);
}
bool MDTransfModQ::calcMatrixCoordElastic(const double k0, std::vector<coord_t> &Coord) const {
  double qx = -m_ex * k0;
  double qy = -m_ey * k0;
  double qz = (1 - m_ez) * k0;
  return applyCoordTransf(qx, qy, qz, Coord);
}
 
std::vector<double> MDTransfModQ::getExtremumPoints(const double eMin, const double eMax, size_t det_num) const {
  std::vector<double> rez(2);
  switch (m_Emode) {
  case (Kernel::DeltaEMode::Elastic): {
    rez[0] = eMin;
    rez[1] = eMax;
    return rez;
  }
  case (Kernel::DeltaEMode::Direct):
  case (Kernel::DeltaEMode::Indirect): {
    double ei = m_eFixed;
    if (m_pEfixedArray)
      ei = double(*(m_pEfixedArray + det_num));
 
    double ez = (m_DetDirecton + det_num)->Z();
    double eps_extr = ei * (1 - ez * ez);
    if (eps_extr > eMin && eps_extr < eMax) {
      rez.resize(3);
      rez[0] = eMin;
      rez[1] = eps_extr;
      rez[2] = eMax;
    } else {
      rez[0] = eMin;
      rez[1] = eMax;
    }
    return rez;
  }
  default: {
    throw std::invalid_argument("Undefined or unsupported energy conversion mode ");
  }
  }
  return rez;
}
 
void MDTransfModQ::initialize(const MDWSDescription &ConvParams) {
  //********** Generic part of initialization, common for elastic and inelastic
  // modes:
  // get transformation matrix (needed for CrystalAsPoder mode)
  m_RotMat = ConvParams.getTransfMatrix();
  m_pEfixedArray = nullptr;
  m_invertRot = false;
  if (!ConvParams.m_PreprDetTable)
    throw(std::runtime_error("The detectors have not been preprocessed but "
                             "they have to before running initialize"));
 
  // get pointer to the positions of the detectors
  std::vector<Kernel::V3D> const &DetDir = ConvParams.m_PreprDetTable->getColVector<Kernel::V3D>("DetDirections");
  m_DetDirecton = &DetDir[0]; //
 
  // get min and max values defined by the algorithm.
  ConvParams.getMinMax(m_DimMin, m_DimMax);
  // m_DimMin/max here are momentums and they are verified on momentum squared
  // base
  if (m_DimMin[0] < 0)
    m_DimMin[0] = 0;
  if (m_DimMax[0] < 0)
    m_DimMax[0] = 0;
 
  // m_DimMin here is a momentum and it is verified on momentum squared base
  m_DimMin[0] *= m_DimMin[0];
  m_DimMax[0] *= m_DimMax[0];
  if (std::fabs(m_DimMin[0] - m_DimMax[0]) < FLT_EPSILON || m_DimMax[0] < m_DimMin[0]) {
    std::string ERR =
        "ModQ coordinate transformation: Min Q^2 value: " + boost::lexical_cast<std::string>(m_DimMin[0]) +
        " is more or equal then Max Q^2 value: " + boost::lexical_cast<std::string>(m_DimMax[0]);
    throw(std::invalid_argument(ERR));
  }
  m_AddDimCoordinates = ConvParams.getAddCoord();
 
  //************   specific part of the initialization, dependent on emode:
  m_Emode = ConvParams.getEMode();
  m_NMatrixDim = getNMatrixDimensions(m_Emode);
  if (m_Emode == Kernel::DeltaEMode::Direct || m_Emode == Kernel::DeltaEMode::Indirect) {
    // energy needed in inelastic case
    volatile auto Ei = ConvParams.m_PreprDetTable->getLogs()->getPropertyValueAsType<double>("Ei");
    m_eFixed = Ei;
    if (Ei != m_eFixed) // Ei is NaN, try Efixed, but the value should be
                        // overridden later
    {
      try {
        m_eFixed = ConvParams.m_PreprDetTable->getLogs()->getPropertyValueAsType<double>("eFixed");
      } catch (...) {
      }
    }
 
    // the wave vector of incident neutrons;
    m_kFixed = sqrt(m_eFixed / PhysicalConstants::E_mev_toNeutronWavenumberSq);
    m_pEfixedArray = nullptr;
    if (m_Emode == static_cast<int>(Kernel::DeltaEMode::Indirect))
      m_pEfixedArray = ConvParams.m_PreprDetTable->getColDataArray<float>("eFixed");
  } else if (m_Emode != Kernel::DeltaEMode::Elastic)
    throw(std::invalid_argument("MDTransfModQ::initialize::Unknown energy conversion mode"));
 
  m_pDetMasks = ConvParams.m_PreprDetTable->getColDataArray<int>("detMask");
}
std::vector<std::string> MDTransfModQ::getDefaultDimID(Kernel::DeltaEMode::Type dEmode,
                                                       API::MatrixWorkspace_const_sptr inWS) const {
  UNUSED_ARG(inWS);
  std::vector<std::string> default_dim_ID;
  switch (dEmode) {
  case (Kernel::DeltaEMode::Elastic): {
    default_dim_ID.resize(1);
    break;
  }
  case (Kernel::DeltaEMode::Direct):
  case (Kernel::DeltaEMode::Indirect): {
    default_dim_ID.resize(2);
    default_dim_ID[1] = "DeltaE";
    break;
  }
  default:
    throw(std::invalid_argument("MDTransfModQ::getDefaultDimID::Unknown energy conversion mode"));
  }
  default_dim_ID[0] = "|Q|";
 
  return default_dim_ID;
}
 
std::vector<std::string> MDTransfModQ::outputUnitID(Kernel::DeltaEMode::Type dEmode,
                                                    API::MatrixWorkspace_const_sptr inWS) const {
  UNUSED_ARG(inWS);
  std::vector<std::string> UnitID = this->getDefaultDimID(dEmode, inWS);
  // TODO: is it really momentum transfer, as MomentumTransfer units are seems
  // bound to elastic mode only (at least accorting to Units description on
  // Wiki)?
  if (dEmode == Kernel::DeltaEMode::Elastic) {
    UnitID[0] = "Momentum";
  } else {
    UnitID[0] = "MomentumTransfer";
  }
  return UnitID;
}
 
MDTransfModQ::MDTransfModQ()
    : m_ex(0), m_ey(0), m_ez(1), m_DetDirecton(nullptr), //,m_NMatrixDim(-1)
      m_NMatrixDim(0),                                   // uninitialized
      m_Emode(Kernel::DeltaEMode::Undefined),            // uninitialized
      m_kFixed(1.), m_eFixed(1.), m_pEfixedArray(nullptr), m_pDetMasks(nullptr), m_invertRot(false) {}
 
std::pair<coord_t, coord_t> MDTransfModQ::getDimBounds(size_t dim) const {
  return std::make_pair(m_DimMin[dim], m_DimMax[dim]);
}
 
std::vector<std::string> MDTransfModQ::getEmodes() const { return Kernel::DeltaEMode::availableTypes(); }
 
void MDTransfModQ::setDisplayNormalization(Mantid::API::IMDWorkspace_sptr mdWorkspace,
                                           Mantid::API::MatrixWorkspace_sptr underlyingWorkspace) const {
  DisplayNormalizationSetter setter;
  auto isQ = true;
  setter(mdWorkspace, underlyingWorkspace, isQ, m_Emode);
}
 
bool MDTransfModQ::applyCoordTransf(double qx, double qy, double qz, std::vector<coord_t> &Coord) const {
  std::array<coord_t, 3> Q{};
  if (m_invertRot) {
    calcMatrixCoordLinSys(qx, qy, qz, Q);
  } else {
    // transformation matrix has to be here for "Crystal AS Powder conversion
    // mode, further specialization possible if "powder" mode defined"
    Q[0] = static_cast<coord_t>(m_RotMat[0] * qx + m_RotMat[1] * qy + m_RotMat[2] * qz);
    Q[1] = static_cast<coord_t>(m_RotMat[3] * qx + m_RotMat[4] * qy + m_RotMat[5] * qz);
    Q[2] = static_cast<coord_t>(m_RotMat[6] * qx + m_RotMat[7] * qy + m_RotMat[8] * qz);
  }
 
  const auto Qsq = Q[0] * Q[0] + Q[1] * Q[1] + Q[2] * Q[2];
  if (Qsq < static_cast<coord_t>(m_DimMin[0]) || Qsq >= static_cast<coord_t>(m_DimMax[0])) {
    return false;
  }
  Coord[0] = static_cast<coord_t>(std::sqrt(Qsq));
 
  return true;
}
 
void MDTransfModQ::calcMatrixCoordLinSys(double qx, double qy, double qz, std::array<coord_t, 3> &Coord) const {
  // For some computations, e.g. continuous rotation in ConvToMDEventsWS, the rotation matrix
  // has to be recomputed multiple times, so it is not inverted prior to calculating the coordinates.
  // Deferring to a linear system solution here makes it slightly more efficient and stable.
  Eigen::Map<const Eigen::Matrix<double, 3, 3, Eigen::RowMajor>> map_rm(m_RotMat.data());
  const Eigen::Vector3d qs(qx, qy, qz);
  const Eigen::PartialPivLU<Eigen::Matrix<double, 3, 3, Eigen::RowMajor>> lu(map_rm);
  Eigen::Vector3d coords = lu.solve(qs);
 
  Coord[0] = static_cast<coord_t>(coords[0]);
  Coord[1] = static_cast<coord_t>(coords[1]);
  Coord[2] = static_cast<coord_t>(coords[2]);
}
 
} // namespace Mantid::MDAlgorithms