#include "MantidMDAlgorithms/MDNormDirectSC.h"
#include "MantidAPI/CommonBinsValidator.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/Run.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceHistory.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/MDEventWorkspace.h"
#include "MantidDataObjects/MDHistoWorkspace.h"
#include "MantidGeometry/Instrument.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/Strings.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/VectorHelper.h"
namespace Mantid {
namespace MDAlgorithms {
using Mantid::API::WorkspaceProperty;
using Mantid::Kernel::Direction;
using namespace Mantid::DataObjects;
using namespace Mantid::API;
using namespace Mantid::Kernel;
namespace {
// function to compare two intersections (h,k,l,Momentum) by Momentum
bool compareMomentum(const std::array<double, 4> &v1,
const std::array<double, 4> &v2) {
return (v1[3] < v2[3]);
}
} // namespace
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(MDNormDirectSC)
/**
* Constructor
*/
MDNormDirectSC::MDNormDirectSC()
: m_normWS(), m_inputWS(), m_hmin(0.0f), m_hmax(0.0f), m_kmin(0.0f),
m_kmax(0.0f), m_lmin(0.0f), m_lmax(0.0f), m_dEmin(0.f), m_dEmax(0.f),
m_Ei(0.), m_ki(0.), m_kfmin(0.), m_kfmax(0.), m_hIntegrated(true),
m_kIntegrated(true), m_lIntegrated(true), m_dEIntegrated(false),
m_rubw(3, 3), m_hIdx(-1), m_kIdx(-1), m_lIdx(-1), m_eIdx(-1), m_hX(),
m_kX(), m_lX(), m_eX(), m_samplePos(), m_beamDir() {}
/// Algorithm's version for identification. @see Algorithm::version
int MDNormDirectSC::version() const { return 1; }
/// Algorithm's category for identification. @see Algorithm::category
const std::string MDNormDirectSC::category() const {
return "MDAlgorithms\\Normalisation";
}
/// Algorithm's summary for use in the GUI and help. @see Algorithm::summary
const std::string MDNormDirectSC::summary() const {
return "Calculate normalization for an MDEvent workspace for single crystal "
"direct geometry inelastic measurement.";
}
/// Algorithm's name for use in the GUI and help. @see Algorithm::name
const std::string MDNormDirectSC::name() const { return "MDNormDirectSC"; }
/**
* Initialize the algorithm's properties.
*/
void MDNormDirectSC::init() {
declareProperty(make_unique<WorkspaceProperty<IMDEventWorkspace>>(
"InputWorkspace", "", Direction::Input),
"An input MDWorkspace.");
std::string dimChars = getDimensionChars();
// --------------- Axis-aligned properties
// ---------------------------------------
for (size_t i = 0; i < dimChars.size(); i++) {
std::string dim(" ");
dim[0] = dimChars[i];
std::string propName = "AlignedDim" + dim;
declareProperty(
Kernel::make_unique<PropertyWithValue<std::string>>(propName, "",
Direction::Input),
"Binning parameters for the " + Strings::toString(i) +
"th dimension.\n"
"Enter it as a comma-separated list of values with the format: "
"'name,minimum,maximum,number_of_bins'. Leave blank for NONE.");
}
auto solidAngleValidator = boost::make_shared<CompositeValidator>();
solidAngleValidator->add<InstrumentValidator>();
solidAngleValidator->add<CommonBinsValidator>();
declareProperty(
make_unique<WorkspaceProperty<>>("SolidAngleWorkspace", "",
Direction::Input, PropertyMode::Optional,
solidAngleValidator),
"An input workspace containing integrated vanadium (a measure of the "
"solid angle).");
declareProperty(make_unique<PropertyWithValue<bool>>("SkipSafetyCheck", false,
Direction::Input),
"If set to true, the algorithm does "
"not check history if the workspace was modified since the"
"ConvertToMD algorithm was run, and assume that the direct "
"geometry inelastic mode is used.");
declareProperty(make_unique<WorkspaceProperty<IMDHistoWorkspace>>(
"TemporaryNormalizationWorkspace", "", Direction::Input,
PropertyMode::Optional),
"An input MDHistoWorkspace used to accumulate normalization "
"from multiple MDEventWorkspaces. If unspecified a blank "
"MDHistoWorkspace will be created.");
declareProperty(make_unique<WorkspaceProperty<IMDHistoWorkspace>>(
"TemporaryDataWorkspace", "", Direction::Input,
PropertyMode::Optional),
"An input MDHistoWorkspace used to accumulate data from "
"multiple MDEventWorkspaces. If unspecified a blank "
"MDHistoWorkspace will be created.");
declareProperty(make_unique<WorkspaceProperty<Workspace>>(
"OutputWorkspace", "", Direction::Output),
"A name for the output data MDHistoWorkspace.");
declareProperty(make_unique<WorkspaceProperty<Workspace>>(
"OutputNormalizationWorkspace", "", Direction::Output),
"A name for the output normalization MDHistoWorkspace.");
}
//----------------------------------------------------------------------------------------------
/**
* Execute the algorithm.
*/
void MDNormDirectSC::exec() {
cacheInputs();
auto outputWS = binInputWS();
convention = Kernel::ConfigService::Instance().getString("Q.convention");
outputWS->setDisplayNormalization(Mantid::API::NoNormalization);
setProperty<Workspace_sptr>("OutputWorkspace", outputWS);
createNormalizationWS(*outputWS);
m_normWS->setDisplayNormalization(Mantid::API::NoNormalization);
setProperty("OutputNormalizationWorkspace", m_normWS);
m_numExptInfos = outputWS->getNumExperimentInfo();
// loop over all experiment infos
for (uint16_t expInfoIndex = 0; expInfoIndex < m_numExptInfos;
expInfoIndex++) {
// Check for other dimensions if we could measure anything in the original
// data
bool skipNormalization = false;
const std::vector<coord_t> otherValues =
getValuesFromOtherDimensions(skipNormalization, expInfoIndex);
const auto affineTrans =
findIntergratedDimensions(otherValues, skipNormalization);
cacheDimensionXValues();
if (!skipNormalization) {
calculateNormalization(otherValues, affineTrans, expInfoIndex);
} else {
g_log.warning("Binning limits are outside the limits of the MDWorkspace. "
"Not applying normalization.");
}
// if more than one experiment info, keep accumulating
m_accumulate = true;
}
// Set the display normalization based on the input workspace
outputWS->setDisplayNormalization(m_inputWS->displayNormalizationHisto());
}
/**
* Set up starting values for cached variables
*/
void MDNormDirectSC::cacheInputs() {
m_inputWS = getProperty("InputWorkspace");
bool skipCheck = getProperty("SkipSafetyCheck");
if (!skipCheck && (inputEnergyMode() != "Direct")) {
throw std::invalid_argument("Invalid energy transfer mode. Algorithm only "
"supports direct geometry spectrometers.");
}
// Min/max dimension values
const auto hdim(m_inputWS->getDimension(0)), kdim(m_inputWS->getDimension(1)),
ldim(m_inputWS->getDimension(2)), edim(m_inputWS->getDimension(3));
m_hmin = hdim->getMinimum();
m_kmin = kdim->getMinimum();
m_lmin = ldim->getMinimum();
m_dEmin = edim->getMinimum();
m_hmax = hdim->getMaximum();
m_kmax = kdim->getMaximum();
m_lmax = ldim->getMaximum();
m_dEmax = edim->getMaximum();
const auto &exptInfoZero = *(m_inputWS->getExperimentInfo(0));
auto source = exptInfoZero.getInstrument()->getSource();
auto sample = exptInfoZero.getInstrument()->getSample();
if (source == nullptr || sample == nullptr) {
throw Kernel::Exception::InstrumentDefinitionError(
"Instrument not sufficiently defined: failed to get source and/or "
"sample");
}
m_samplePos = sample->getPos();
m_beamDir = m_samplePos - source->getPos();
m_beamDir.normalize();
double originaldEmin = exptInfoZero.run().getBinBoundaries().front();
double originaldEmax = exptInfoZero.run().getBinBoundaries().back();
if (exptInfoZero.run().hasProperty("Ei")) {
Kernel::Property *eiprop = exptInfoZero.run().getProperty("Ei");
m_Ei = boost::lexical_cast<double>(eiprop->value());
if (m_Ei <= 0) {
throw std::invalid_argument("Ei stored in the workspace is not positive");
}
} else {
throw std::invalid_argument("Could not find Ei value in the workspace.");
}
double eps = 1e-7;
if (m_Ei - originaldEmin < eps) {
originaldEmin = m_Ei - eps;
}
if (m_Ei - originaldEmax < eps) {
originaldEmax = m_Ei - 1e-7;
}
if (originaldEmin == originaldEmax) {
throw std::runtime_error("The limits of the original workspace used in "
"ConvertToMD are incorrect");
}
const double energyToK = 8.0 * M_PI * M_PI * PhysicalConstants::NeutronMass *
PhysicalConstants::meV * 1e-20 /
(PhysicalConstants::h * PhysicalConstants::h);
m_ki = std::sqrt(energyToK * m_Ei);
m_kfmin = std::sqrt(energyToK * (m_Ei - originaldEmin));
m_kfmax = std::sqrt(energyToK * (m_Ei - originaldEmax));
}
/**
* Currently looks for the ConvertToMD algorithm in the history
* @return A string donating the energy transfer mode of the input workspace
*/
std::string MDNormDirectSC::inputEnergyMode() const {
const auto &hist = m_inputWS->getHistory();
const size_t nalgs = hist.size();
const auto &lastAlgHist = hist.getAlgorithmHistory(nalgs - 1);
const auto &penultimateAlgHist = hist.getAlgorithmHistory(nalgs - 2);
std::string emode;
if (lastAlgHist->name() == "ConvertToMD") {
emode = lastAlgHist->getPropertyValue("dEAnalysisMode");
} else if ((lastAlgHist->name() == "Load" ||
lastAlgHist->name() == "LoadMD") &&
penultimateAlgHist->name() == "ConvertToMD") {
// get dEAnalysisMode
emode = penultimateAlgHist->getPropertyValue("dEAnalysisMode");
} else {
throw std::invalid_argument("The last algorithm in the history of the "
"input workspace is not ConvertToMD");
}
return emode;
}
/**
* Runs the BinMD algorithm on the input to provide the output workspace
* All slicing algorithm properties are passed along
* @return MDHistoWorkspace as a result of the binning
*/
MDHistoWorkspace_sptr MDNormDirectSC::binInputWS() {
const auto &props = getProperties();
IAlgorithm_sptr binMD = createChildAlgorithm("BinMD", 0.0, 0.3);
binMD->setPropertyValue("AxisAligned", "1");
for (auto prop : props) {
const auto &propName = prop->name();
if (propName != "SolidAngleWorkspace" &&
propName != "TemporaryNormalizationWorkspace" &&
propName != "OutputNormalizationWorkspace" &&
propName != "SkipSafetyCheck") {
binMD->setPropertyValue(propName, prop->value());
}
}
binMD->executeAsChildAlg();
Workspace_sptr outputWS = binMD->getProperty("OutputWorkspace");
return boost::dynamic_pointer_cast<MDHistoWorkspace>(outputWS);
}
/**
* Create & cached the normalization workspace
* @param dataWS The binned workspace that will be used for the data
*/
void MDNormDirectSC::createNormalizationWS(const MDHistoWorkspace &dataWS) {
// Copy the MDHisto workspace, and change signals and errors to 0.
boost::shared_ptr<IMDHistoWorkspace> tmp =
this->getProperty("TemporaryNormalizationWorkspace");
m_normWS = boost::dynamic_pointer_cast<MDHistoWorkspace>(tmp);
if (!m_normWS) {
m_normWS = dataWS.clone();
m_normWS->setTo(0., 0., 0.);
} else {
m_accumulate = true;
}
}
/**
* Retrieve logged values from non-HKL dimensions
* @param skipNormalization [InOut] Updated to false if any values are outside
* range measured by input workspace
* @param expInfoIndex current experiment info index
* @return A vector of values from other dimensions to be include in normalized
* MD position calculation
*/
std::vector<coord_t>
MDNormDirectSC::getValuesFromOtherDimensions(bool &skipNormalization,
uint16_t expInfoIndex) const {
const auto ¤tRun = m_inputWS->getExperimentInfo(expInfoIndex)->run();
std::vector<coord_t> otherDimValues;
for (size_t i = 4; i < m_inputWS->getNumDims(); i++) {
const auto dimension = m_inputWS->getDimension(i);
float dimMin = static_cast<float>(dimension->getMinimum());
float dimMax = static_cast<float>(dimension->getMaximum());
auto *dimProp = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(
currentRun.getProperty(dimension->getName()));
if (dimProp) {
coord_t value = static_cast<coord_t>(dimProp->firstValue());
otherDimValues.push_back(value);
// in the original MD data no time was spent measuring between dimMin and
// dimMax
if (value < dimMin || value > dimMax) {
skipNormalization = true;
}
}
}
return otherDimValues;
}
/**
* Checks the normalization workspace against the indices of the original
* dimensions.
* If not found, the corresponding dimension is integrated
* @param otherDimValues Values from non-HKL dimensions
* @param skipNormalization [InOut] Sets the flag true if normalization values
* are outside of original inputs
* @return Affine trasform matrix
*/
Kernel::Matrix<coord_t> MDNormDirectSC::findIntergratedDimensions(
const std::vector<coord_t> &otherDimValues, bool &skipNormalization) {
// Get indices of the original dimensions in the output workspace,
// and if not found, the corresponding dimension is integrated
Kernel::Matrix<coord_t> affineMat =
m_normWS->getTransformFromOriginal(0)->makeAffineMatrix();
const size_t nrm1 = affineMat.numRows() - 1;
const size_t ncm1 = affineMat.numCols() - 1;
for (size_t row = 0; row < nrm1; row++) // affine matrix, ignore last row
{
const auto dimen = m_normWS->getDimension(row);
const auto dimMin(dimen->getMinimum()), dimMax(dimen->getMaximum());
if (affineMat[row][0] == 1.) {
m_hIntegrated = false;
m_hIdx = row;
m_hmin = std::max(m_hmin, dimMin);
m_hmax = std::min(m_hmax, dimMax);
if (m_hmin > dimMax || m_hmax < dimMin) {
skipNormalization = true;
}
}
if (affineMat[row][1] == 1.) {
m_kIntegrated = false;
m_kIdx = row;
m_kmin = std::max(m_kmin, dimMin);
m_kmax = std::min(m_kmax, dimMax);
if (m_kmin > dimMax || m_kmax < dimMin) {
skipNormalization = true;
}
}
if (affineMat[row][2] == 1.) {
m_lIntegrated = false;
m_lIdx = row;
m_lmin = std::max(m_lmin, dimMin);
m_lmax = std::min(m_lmax, dimMax);
if (m_lmin > dimMax || m_lmax < dimMin) {
skipNormalization = true;
}
}
if (affineMat[row][3] == 1.) {
m_dEIntegrated = false;
m_eIdx = row;
m_dEmin = std::max(m_dEmin, dimMin);
m_dEmax = std::min(m_dEmax, dimMax);
if (m_dEmin > dimMax || m_dEmax < dimMin) {
skipNormalization = true;
}
}
for (size_t col = 4; col < ncm1; col++) // affine matrix, ignore last column
{
if (affineMat[row][col] == 1.) {
double val = otherDimValues.at(col - 3);
if (val > dimMax || val < dimMin) {
skipNormalization = true;
}
}
}
}
return affineMat;
}
/**
* Stores the X values from each H,K,L,E dimension as member variables
* Energy dimension is transformed to final wavevector.
*/
void MDNormDirectSC::cacheDimensionXValues() {
constexpr double energyToK = 8.0 * M_PI * M_PI *
PhysicalConstants::NeutronMass *
PhysicalConstants::meV * 1e-20 /
(PhysicalConstants::h * PhysicalConstants::h);
if (!m_hIntegrated) {
auto &hDim = *m_normWS->getDimension(m_hIdx);
m_hX.resize(hDim.getNBoundaries());
for (size_t i = 0; i < m_hX.size(); ++i) {
m_hX[i] = hDim.getX(i);
}
}
if (!m_kIntegrated) {
auto &kDim = *m_normWS->getDimension(m_kIdx);
m_kX.resize(kDim.getNBoundaries());
for (size_t i = 0; i < m_kX.size(); ++i) {
m_kX[i] = kDim.getX(i);
}
}
if (!m_lIntegrated) {
auto &lDim = *m_normWS->getDimension(m_lIdx);
m_lX.resize(lDim.getNBoundaries());
for (size_t i = 0; i < m_lX.size(); ++i) {
m_lX[i] = lDim.getX(i);
}
}
if (!m_dEIntegrated) {
// NOTE: store k final instead
auto &eDim = *m_normWS->getDimension(m_eIdx);
m_eX.resize(eDim.getNBoundaries());
for (size_t i = 0; i < m_eX.size(); ++i) {
double temp = m_Ei - eDim.getX(i);
temp = std::max(temp, 0.);
m_eX[i] = std::sqrt(energyToK * temp);
}
}
}
/**
* Computed the normalization for the input workspace. Results are stored in
* m_normWS
* @param otherValues non HKLE dimensions
* @param affineTrans affine matrix
* @param expInfoIndex current experiment info index
*/
void MDNormDirectSC::calculateNormalization(
const std::vector<coord_t> &otherValues,
const Kernel::Matrix<coord_t> &affineTrans, uint16_t expInfoIndex) {
constexpr double energyToK = 8.0 * M_PI * M_PI *
PhysicalConstants::NeutronMass *
PhysicalConstants::meV * 1e-20 /
(PhysicalConstants::h * PhysicalConstants::h);
const auto ¤tExptInfo = *(m_inputWS->getExperimentInfo(expInfoIndex));
using VectorDoubleProperty = Kernel::PropertyWithValue<std::vector<double>>;
auto *rubwLog = dynamic_cast<VectorDoubleProperty *>(
currentExptInfo.getLog("RUBW_MATRIX"));
if (!rubwLog) {
throw std::runtime_error(
"Wokspace does not contain a log entry for the RUBW matrix."
"Cannot continue.");
} else {
Kernel::DblMatrix rubwValue(
(*rubwLog)()); // includes the 2*pi factor but not goniometer for now :)
m_rubw = currentExptInfo.run().getGoniometerMatrix() * rubwValue;
m_rubw.Invert();
}
const double protonCharge = currentExptInfo.run().getProtonCharge();
const auto &spectrumInfo = currentExptInfo.spectrumInfo();
// Mapping
const int64_t ndets = static_cast<int64_t>(spectrumInfo.size());
bool haveSA = false;
API::MatrixWorkspace_const_sptr solidAngleWS =
getProperty("SolidAngleWorkspace");
detid2index_map solidAngDetToIdx;
if (solidAngleWS != nullptr) {
haveSA = true;
solidAngDetToIdx = solidAngleWS->getDetectorIDToWorkspaceIndexMap();
}
const size_t vmdDims = 4;
std::vector<std::atomic<signal_t>> signalArray(m_normWS->getNPoints());
std::vector<std::array<double, 4>> intersections;
std::vector<coord_t> pos, posNew;
double progStep = 0.7 / m_numExptInfos;
auto prog =
make_unique<API::Progress>(this, 0.3 + progStep * expInfoIndex,
0.3 + progStep * (expInfoIndex + 1.), ndets);
// cppcheck-suppress syntaxError
PRAGMA_OMP(parallel for private(intersections, pos, posNew))
for (int64_t i = 0; i < ndets; i++) {
PARALLEL_START_INTERUPT_REGION
if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i) ||
spectrumInfo.isMasked(i)) {
continue;
}
const auto &detector = spectrumInfo.detector(i);
double theta = detector.getTwoTheta(m_samplePos, m_beamDir);
double phi = detector.getPhi();
// If the detector is a group, this should be the ID of the first detector
const auto detID = detector.getID();
// Intersections
this->calculateIntersections(intersections, theta, phi);
if (intersections.empty())
continue;
// Get solid angle for this contribution
double solid = protonCharge;
if (haveSA) {
solid =
solidAngleWS->y(solidAngDetToIdx.find(detID)->second)[0] * protonCharge;
}
// Compute final position in HKL
// pre-allocate for efficiency and copy non-hkl dim values into place
pos.resize(vmdDims + otherValues.size() + 1);
std::copy(otherValues.begin(), otherValues.end(), pos.begin() + vmdDims);
pos.push_back(1.);
auto intersectionsBegin = intersections.begin();
for (auto it = intersectionsBegin + 1; it != intersections.end(); ++it) {
const auto &curIntSec = *it;
const auto &prevIntSec = *(it - 1);
// the full vector isn't used so compute only what is necessary
double delta =
(curIntSec[3] * curIntSec[3] - prevIntSec[3] * prevIntSec[3]) /
energyToK;
if (delta < 1e-10)
continue; // Assume zero contribution if difference is small
// Average between two intersections for final position
std::transform(curIntSec.data(), curIntSec.data() + vmdDims,
prevIntSec.data(), pos.begin(),
[](const double rhs, const double lhs) {
return static_cast<coord_t>(0.5 * (rhs + lhs));
});
// transform kf to energy transfer
pos[3] = static_cast<coord_t>(m_Ei - pos[3] * pos[3] / energyToK);
affineTrans.multiplyPoint(pos, posNew);
size_t linIndex = m_normWS->getLinearIndexAtCoord(posNew.data());
if (linIndex == size_t(-1))
continue;
// signal = integral between two consecutive intersections *solid angle
// *PC
double signal = solid * delta;
Mantid::Kernel::AtomicOp(signalArray[linIndex], signal,
std::plus<signal_t>());
}
prog->report();
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
if (m_accumulate) {
std::transform(
signalArray.cbegin(), signalArray.cend(), m_normWS->getSignalArray(),
m_normWS->getSignalArray(),
[](const std::atomic<signal_t> &a, const signal_t &b) { return a + b; });
} else {
std::copy(signalArray.cbegin(), signalArray.cend(),
m_normWS->getSignalArray());
}
}
/**
* Calculate the points of intersection for the given detector with cuboid
* surrounding the
* detector position in HKL
* @param intersections A list of intersections in HKL space
* @param theta Polar angle with detector
* @param phi Azimuthal angle with detector
*/
void MDNormDirectSC::calculateIntersections(
std::vector<std::array<double, 4>> &intersections, const double theta,
const double phi) {
V3D qout(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta)),
qin(0., 0., m_ki);
qout = m_rubw * qout;
qin = m_rubw * qin;
if (convention == "Crystallography") {
qout *= -1;
qin *= -1;
}
double hStart = qin.X() - qout.X() * m_kfmin,
hEnd = qin.X() - qout.X() * m_kfmax;
double kStart = qin.Y() - qout.Y() * m_kfmin,
kEnd = qin.Y() - qout.Y() * m_kfmax;
double lStart = qin.Z() - qout.Z() * m_kfmin,
lEnd = qin.Z() - qout.Z() * m_kfmax;
double eps = 1e-10;
auto hNBins = m_hX.size();
auto kNBins = m_kX.size();
auto lNBins = m_lX.size();
auto eNBins = m_eX.size();
intersections.clear();
intersections.reserve(hNBins + kNBins + lNBins + eNBins +
8); // 8 is 3*(min,max for each Q component)+kfmin+kfmax
// calculate intersections with planes perpendicular to h
if (fabs(hStart - hEnd) > eps) {
double fmom = (m_kfmax - m_kfmin) / (hEnd - hStart);
double fk = (kEnd - kStart) / (hEnd - hStart);
double fl = (lEnd - lStart) / (hEnd - hStart);
if (!m_hIntegrated) {
for (size_t i = 0; i < hNBins; i++) {
double hi = m_hX[i];
if ((hi >= m_hmin) && (hi <= m_hmax) &&
((hStart - hi) * (hEnd - hi) < 0)) {
// if hi is between hStart and hEnd, then ki and li will be between
// kStart, kEnd and lStart, lEnd and momi will be between m_kfmin and
// m_kfmax
double ki = fk * (hi - hStart) + kStart;
double li = fl * (hi - hStart) + lStart;
if ((ki >= m_kmin) && (ki <= m_kmax) && (li >= m_lmin) &&
(li <= m_lmax)) {
double momi = fmom * (hi - hStart) + m_kfmin;
intersections.push_back({{hi, ki, li, momi}});
}
}
}
}
double momhMin = fmom * (m_hmin - hStart) + m_kfmin;
if ((momhMin - m_kfmin) * (momhMin - m_kfmax) < 0) // m_kfmin>m_kfmax
{
// khmin and lhmin
double khmin = fk * (m_hmin - hStart) + kStart;
double lhmin = fl * (m_hmin - hStart) + lStart;
if ((khmin >= m_kmin) && (khmin <= m_kmax) && (lhmin >= m_lmin) &&
(lhmin <= m_lmax)) {
intersections.push_back({{m_hmin, khmin, lhmin, momhMin}});
}
}
double momhMax = fmom * (m_hmax - hStart) + m_kfmin;
if ((momhMax - m_kfmin) * (momhMax - m_kfmax) <= 0) {
// khmax and lhmax
double khmax = fk * (m_hmax - hStart) + kStart;
double lhmax = fl * (m_hmax - hStart) + lStart;
if ((khmax >= m_kmin) && (khmax <= m_kmax) && (lhmax >= m_lmin) &&
(lhmax <= m_lmax)) {
intersections.push_back({{m_hmax, khmax, lhmax, momhMax}});
}
}
}
// calculate intersections with planes perpendicular to k
if (fabs(kStart - kEnd) > eps) {
double fmom = (m_kfmax - m_kfmin) / (kEnd - kStart);
double fh = (hEnd - hStart) / (kEnd - kStart);
double fl = (lEnd - lStart) / (kEnd - kStart);
if (!m_kIntegrated) {
for (size_t i = 0; i < kNBins; i++) {
double ki = m_kX[i];
if ((ki >= m_kmin) && (ki <= m_kmax) &&
((kStart - ki) * (kEnd - ki) < 0)) {
// if ki is between kStart and kEnd, then hi and li will be between
// hStart, hEnd and lStart, lEnd and momi will be between m_kfmin and
// m_kfmax
double hi = fh * (ki - kStart) + hStart;
double li = fl * (ki - kStart) + lStart;
if ((hi >= m_hmin) && (hi <= m_hmax) && (li >= m_lmin) &&
(li <= m_lmax)) {
double momi = fmom * (ki - kStart) + m_kfmin;
intersections.push_back({{hi, ki, li, momi}});
}
}
}
}
double momkMin = fmom * (m_kmin - kStart) + m_kfmin;
if ((momkMin - m_kfmin) * (momkMin - m_kfmax) < 0) {
// hkmin and lkmin
double hkmin = fh * (m_kmin - kStart) + hStart;
double lkmin = fl * (m_kmin - kStart) + lStart;
if ((hkmin >= m_hmin) && (hkmin <= m_hmax) && (lkmin >= m_lmin) &&
(lkmin <= m_lmax)) {
intersections.push_back({{hkmin, m_kmin, lkmin, momkMin}});
}
}
double momkMax = fmom * (m_kmax - kStart) + m_kfmin;
if ((momkMax - m_kfmin) * (momkMax - m_kfmax) <= 0) {
// hkmax and lkmax
double hkmax = fh * (m_kmax - kStart) + hStart;
double lkmax = fl * (m_kmax - kStart) + lStart;
if ((hkmax >= m_hmin) && (hkmax <= m_hmax) && (lkmax >= m_lmin) &&
(lkmax <= m_lmax)) {
intersections.push_back({{hkmax, m_kmax, lkmax, momkMax}});
}
}
}
// calculate intersections with planes perpendicular to l
if (fabs(lStart - lEnd) > eps) {
double fmom = (m_kfmax - m_kfmin) / (lEnd - lStart);
double fh = (hEnd - hStart) / (lEnd - lStart);
double fk = (kEnd - kStart) / (lEnd - lStart);
if (!m_lIntegrated) {
for (size_t i = 0; i < lNBins; i++) {
double li = m_lX[i];
if ((li >= m_lmin) && (li <= m_lmax) &&
((lStart - li) * (lEnd - li) < 0)) {
double hi = fh * (li - lStart) + hStart;
double ki = fk * (li - lStart) + kStart;
if ((hi >= m_hmin) && (hi <= m_hmax) && (ki >= m_kmin) &&
(ki <= m_kmax)) {
double momi = fmom * (li - lStart) + m_kfmin;
intersections.push_back({{hi, ki, li, momi}});
}
}
}
}
double momlMin = fmom * (m_lmin - lStart) + m_kfmin;
if ((momlMin - m_kfmin) * (momlMin - m_kfmax) <= 0) {
// hlmin and klmin
double hlmin = fh * (m_lmin - lStart) + hStart;
double klmin = fk * (m_lmin - lStart) + kStart;
if ((hlmin >= m_hmin) && (hlmin <= m_hmax) && (klmin >= m_kmin) &&
(klmin <= m_kmax)) {
intersections.push_back({{hlmin, klmin, m_lmin, momlMin}});
}
}
double momlMax = fmom * (m_lmax - lStart) + m_kfmin;
if ((momlMax - m_kfmin) * (momlMax - m_kfmax) < 0) {
// hlmax and klmax
double hlmax = fh * (m_lmax - lStart) + hStart;
double klmax = fk * (m_lmax - lStart) + kStart;
if ((hlmax >= m_hmin) && (hlmax <= m_hmax) && (klmax >= m_kmin) &&
(klmax <= m_kmax)) {
intersections.push_back({{hlmax, klmax, m_lmax, momlMax}});
}
}
}
// intersections with dE
if (!m_dEIntegrated) {
for (size_t i = 0; i < eNBins; i++) {
double kfi = m_eX[i];
if ((kfi - m_kfmin) * (kfi - m_kfmax) <= 0) {
double h = qin.X() - qout.X() * kfi;
double k = qin.Y() - qout.Y() * kfi;
double l = qin.Z() - qout.Z() * kfi;
if ((h >= m_hmin) && (h <= m_hmax) && (k >= m_kmin) && (k <= m_kmax) &&
(l >= m_lmin) && (l <= m_lmax)) {
intersections.push_back({{h, k, l, kfi}});
}
}
}
}
// endpoints
if ((hStart >= m_hmin) && (hStart <= m_hmax) && (kStart >= m_kmin) &&
(kStart <= m_kmax) && (lStart >= m_lmin) && (lStart <= m_lmax)) {
intersections.push_back({{hStart, kStart, lStart, m_kfmin}});
}
if ((hEnd >= m_hmin) && (hEnd <= m_hmax) && (kEnd >= m_kmin) &&
(kEnd <= m_kmax) && (lEnd >= m_lmin) && (lEnd <= m_lmax)) {
intersections.push_back({{hEnd, kEnd, lEnd, m_kfmax}});
}
// sort intersections by final momentum
std::stable_sort(intersections.begin(), intersections.end(), compareMomentum);
}
} // namespace MDAlgorithms
} // namespace Mantid