MDTransfQ3D.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/MDTransfQ3D.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/RegistrationHelper.h"
 
namespace Mantid::MDAlgorithms {
// register the class, whith conversion factory under Q3D name
DECLARE_MD_TRANSFID(MDTransfQ3D, Q3D)
 
 
unsigned int MDTransfQ3D::getNMatrixDimensions(Kernel::DeltaEMode::Type mode,
                                               API::MatrixWorkspace_const_sptr inWS) const {
  UNUSED_ARG(inWS);
  switch (mode) {
  case (Kernel::DeltaEMode::Direct):
    return 4;
  case (Kernel::DeltaEMode::Indirect):
    return 4;
  case (Kernel::DeltaEMode::Elastic):
    return 3;
  default:
    throw(std::invalid_argument("Unknow or unsupported energy conversion mode"));
  }
}
 
bool MDTransfQ3D::calcMatrixCoord(const double &deltaEOrK0, std::vector<coord_t> &Coord, double &s, double &err) const {
  if (m_Emode == Kernel::DeltaEMode::Elastic) {
    return calcMatrixCoord3DElastic(deltaEOrK0, Coord, s, err);
  } else {
    return calcMatrixCoord3DInelastic(deltaEOrK0, Coord);
  }
}
 
bool MDTransfQ3D::calcMatrixCoord3DInelastic(const double &deltaE, std::vector<coord_t> &Coord) const {
  Coord[3] = static_cast<coord_t>(deltaE);
  if (Coord[3] < m_DimMin[3] || Coord[3] >= m_DimMax[3])
    return false;
 
  // 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 (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 {
    qx = -m_ex * m_kFixed;
    qy = -m_ey * m_kFixed;
    const double kInitial = sqrt((m_eFixed + deltaE) / PhysicalConstants::E_mev_toNeutronWavenumberSq);
    qz = kInitial - m_ez * m_kFixed;
  }
 
  if (convention == "Crystallography") {
    qx = -qx;
    qy = -qy;
    qz = -qz;
  }
 
  Coord[0] = static_cast<coord_t>(m_RotMat[0] * qx + m_RotMat[1] * qy + m_RotMat[2] * qz);
 
  if (Coord[0] < m_DimMin[0] || Coord[0] >= m_DimMax[0])
    return false;
  Coord[1] = static_cast<coord_t>(m_RotMat[3] * qx + m_RotMat[4] * qy + m_RotMat[5] * qz);
  if (Coord[1] < m_DimMin[1] || Coord[1] >= m_DimMax[1])
    return false;
  Coord[2] = static_cast<coord_t>(m_RotMat[6] * qx + m_RotMat[7] * qy + m_RotMat[8] * qz);
  if (Coord[2] < m_DimMin[2] || Coord[2] >= m_DimMax[2])
    return false;
 
  if (std::sqrt(Coord[0] * Coord[0] + Coord[1] * Coord[1] + Coord[2] * Coord[2]) < m_AbsMin)
    return false;
 
  return true;
}
 
bool MDTransfQ3D::calcMatrixCoord3DElastic(const double &k0, std::vector<coord_t> &Coord, double &signal,
                                           double &errSq) const {
 
  double qx = -m_ex * k0;
  double qy = -m_ey * k0;
  double qz = (1 - m_ez) * k0;
  if (convention == "Crystallography") {
    qx = -qx;
    qy = -qy;
    qz = -qz;
  }
 
  // Dimension limits have to be converted to coord_t, otherwise floating point
  // error will cause valid events to be discarded.
  Coord[0] = static_cast<coord_t>(m_RotMat[0] * qx + m_RotMat[1] * qy + m_RotMat[2] * qz);
  if (Coord[0] < static_cast<coord_t>(m_DimMin[0]) || Coord[0] >= static_cast<coord_t>(m_DimMax[0]))
    return false;
 
  Coord[1] = static_cast<coord_t>(m_RotMat[3] * qx + m_RotMat[4] * qy + m_RotMat[5] * qz);
  if (Coord[1] < static_cast<coord_t>(m_DimMin[1]) || Coord[1] >= static_cast<coord_t>(m_DimMax[1]))
    return false;
 
  Coord[2] = static_cast<coord_t>(m_RotMat[6] * qx + m_RotMat[7] * qy + m_RotMat[8] * qz);
  if (Coord[2] < static_cast<coord_t>(m_DimMin[2]) || Coord[2] >= static_cast<coord_t>(m_DimMax[2]))
    return false;
 
  if (std::sqrt(Coord[0] * Coord[0] + Coord[1] * Coord[1] + Coord[2] * Coord[2]) < m_AbsMin)
    return false;
 
  /*Apply Lorentz corrections if necessary */
  if (m_isLorentzCorrected) {
    double kdash = k0 / (2 * M_PI);
    double correct = m_SinThetaSq * kdash * kdash * kdash * kdash;
    signal *= correct;
    errSq *= (correct * correct);
  }
 
  return true;
}
 
std::vector<double> MDTransfQ3D::getExtremumPoints(const double xMin, const double xMax, size_t det_num) const {
  UNUSED_ARG(det_num);
 
  std::vector<double> rez(2);
  rez[0] = xMin;
  rez[1] = xMax;
 
  return rez;
}
 
bool MDTransfQ3D::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 Lorentz-corrected, retrieve the sin(Theta)^2 for the detector;
  if (m_isLorentzCorrected)
    m_SinThetaSq = *(m_SinThetaSqArray + i);
  // 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 masks are defined and detector masked -- no further calculations
  if (m_pDetMasks) {
    if (*(m_pDetMasks + i) > 0)
      return false;
  }
 
  return true;
}
 
void MDTransfQ3D::initialize(const MDWSDescription &ConvParams) {
  m_pEfixedArray = nullptr;
  m_pDetMasks = nullptr;
  convention = Kernel::ConfigService::Instance().getString("Q.convention");
  //********** Generic part of initialization, common for elastic and inelastic
  // modes:
  // get transformation matrix (needed for CrystalAsPoder mode)
  m_RotMat = ConvParams.getTransfMatrix();
 
  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 preprocessed 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);
  // get additional coordinates which are
  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
    m_eFixed = ConvParams.m_PreprDetTable->getLogs()->getPropertyValueAsType<double>("Ei");
    // 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::runtime_error("MDTransfQ3D::initialize::Unknown or "
                               "unsupported energy conversion mode"));
    // check if we need to calculate Lorentz corrections and if we do, prepare
    // values for their precalculation:
    m_isLorentzCorrected = ConvParams.isLorentsCorrections();
    if (m_isLorentzCorrected) {
      auto &TwoTheta = ConvParams.m_PreprDetTable->getColVector<double>("TwoTheta");
      SinThetaSq.resize(TwoTheta.size());
      for (size_t i = 0; i < TwoTheta.size(); i++) {
        double sth = sin(0.5 * TwoTheta[i]);
        SinThetaSq[i] = sth * sth;
      }
      m_SinThetaSqArray = &SinThetaSq[0];
      if (!m_SinThetaSqArray)
        throw(std::runtime_error("MDTransfQ3D::initialize::Uninitilized "
                                 "Sin(Theta)^2 array for calculating Lorentz "
                                 "corrections"));
    }
  }
  // use detectors masks untill signals are masked by 0 instead of NaN
  m_pDetMasks = ConvParams.m_PreprDetTable->getColDataArray<int>("detMask");
  m_AbsMin = ConvParams.absMin();
}
std::vector<std::string> MDTransfQ3D::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(3);
    break;
  }
  case (Kernel::DeltaEMode::Direct):
  case (Kernel::DeltaEMode::Indirect): {
    default_dim_ID.resize(4);
    default_dim_ID[3] = "DeltaE";
    break;
  }
  default:
    throw(std::invalid_argument("MDTransfQ3D::getDefaultDimID::Unknown energy conversion mode"));
  }
  default_dim_ID[0] = "Q1";
  default_dim_ID[1] = "Q2";
  default_dim_ID[2] = "Q3";
 
  return default_dim_ID;
}
 
std::vector<std::string> MDTransfQ3D::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)?
  std::string kUnits("MomentumTransfer");
  if (dEmode == Kernel::DeltaEMode::Elastic)
    kUnits = "Momentum";
 
  UnitID[0] = kUnits;
  UnitID[1] = kUnits;
  UnitID[2] = kUnits;
  return UnitID;
}
 
MDTransfQ3D::MDTransfQ3D()
    : m_isLorentzCorrected(false), m_SinThetaSqArray(nullptr), SinThetaSq(), m_SinThetaSq(0.), m_AbsMin(0.) {}
 
} // namespace Mantid::MDAlgorithms