VesuvioCalculateMS.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 "MantidCurveFitting/Algorithms/VesuvioCalculateMS.h"
// Use helpers for storing detector/resolution parameters
#include "MantidCurveFitting/Algorithms/ConvertToYSpace.h"
#include "MantidCurveFitting/Functions/VesuvioResolution.h"
 
#include "MantidAPI/Axis.h"
#include "MantidAPI/Sample.h"
#include "MantidAPI/SampleShapeValidator.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
 
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/ComponentInfo.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidGeometry/Objects/Track.h"
 
#include "MantidKernel/ArrayLengthValidator.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/VectorHelper.h"
 
#include <memory>
 
namespace Mantid::CurveFitting::Algorithms {
using namespace API;
using namespace Kernel;
using namespace CurveFitting;
using namespace CurveFitting::Functions;
using Geometry::Track;
 
namespace {
const size_t MAX_SCATTER_PT_TRIES = 500;
const double MASS_TO_MEV = 0.5 * PhysicalConstants::NeutronMass / PhysicalConstants::meV;
} // end anonymous namespace
 
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(VesuvioCalculateMS)
 
 
VesuvioCalculateMS::VesuvioCalculateMS()
    : Algorithm(), m_acrossIdx(0), m_upIdx(1), m_beamIdx(3), m_beamDir(), m_srcR2(0.0), m_halfSampleHeight(0.0),
      m_halfSampleWidth(0.0), m_halfSampleThick(0.0), m_sampleShape(nullptr), m_sampleProps(nullptr), m_detHeight(-1.0),
      m_detWidth(-1.0), m_detThick(-1.0), m_tmin(-1.0), m_tmax(-1.0), m_delt(-1.0), m_foilRes(-1.0), m_nscatters(0),
      m_nruns(0), m_nevents(0), m_progress(nullptr), m_inputWS() {}
 
void VesuvioCalculateMS::init() {
  // Inputs
  auto inputWSValidator = std::make_shared<CompositeValidator>();
  inputWSValidator->add<WorkspaceUnitValidator>("TOF");
  inputWSValidator->add<SampleShapeValidator>();
  declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, inputWSValidator),
                  "Input workspace to be corrected, in units of TOF.");
 
  // -- Sample --
  auto positiveInt = std::make_shared<Kernel::BoundedValidator<int>>();
  positiveInt->setLower(1);
  declareProperty("NoOfMasses", -1, positiveInt, "The number of masses contained within the sample");
 
  auto positiveNonZero = std::make_shared<BoundedValidator<double>>();
  positiveNonZero->setLower(0.0);
  positiveNonZero->setLowerExclusive(true);
  declareProperty("SampleDensity", -1.0, positiveNonZero, "The density of the sample in gm/cm^3");
 
  auto nonEmptyArray = std::make_shared<ArrayLengthValidator<double>>();
  nonEmptyArray->setLengthMin(3);
  declareProperty(std::make_unique<ArrayProperty<double>>("AtomicProperties", nonEmptyArray),
                  "Atomic properties of masses within the sample. "
                  "The expected format is 3 consecutive values per mass: "
                  "mass(amu), cross-section (barns) & s.d of Compton profile.");
  setPropertyGroup("NoOfMasses", "Sample");
  setPropertyGroup("SampleDensity", "Sample");
  setPropertyGroup("AtomicProperties", "Sample");
 
  // -- Beam --
  declareProperty("BeamRadius", -1.0, positiveNonZero, "Radius, in cm, of beam");
 
  // -- Algorithm --
  declareProperty("Seed", 123456789, positiveInt, "Seed the random number generator with this value");
  declareProperty("NumScatters", 3, positiveInt, "Number of scattering orders to calculate");
  declareProperty("NumRuns", 10, positiveInt, "Number of simulated runs per spectrum");
  declareProperty("NumEventsPerRun", 50000, positiveInt, "Number of events per run");
  setPropertyGroup("Seed", "Algorithm");
  setPropertyGroup("NumScatters", "Algorithm");
  setPropertyGroup("NumRuns", "Algorithm");
  setPropertyGroup("NumEventsPerRun", "Algorithm");
 
  // Outputs
  declareProperty(std::make_unique<WorkspaceProperty<>>("TotalScatteringWS", "", Direction::Output),
                  "Workspace to store the calculated total scattering counts");
  declareProperty(std::make_unique<WorkspaceProperty<>>("MultipleScatteringWS", "", Direction::Output),
                  "Workspace to store the calculated multiple scattering counts summed for "
                  "all orders");
}
 
void VesuvioCalculateMS::exec() {
  m_inputWS = getProperty("InputWorkspace");
  cacheInputs();
 
  // Create new workspaces
  MatrixWorkspace_sptr totalsc = WorkspaceFactory::Instance().create(m_inputWS);
  MatrixWorkspace_sptr multsc = WorkspaceFactory::Instance().create(m_inputWS);
 
  // Setup progress
  const size_t nhist = m_inputWS->getNumberHistograms();
  m_progress = std::make_unique<Progress>(this, 0.0, 1.0, nhist * m_nruns * 2);
  const auto &spectrumInfo = m_inputWS->spectrumInfo();
 
  PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWS))
  for (int64_t i = 0; i < static_cast<int64_t>(nhist); ++i) {
    PARALLEL_START_INTERRUPT_REGION
 
    // set common X-values
    totalsc->setSharedX(i, m_inputWS->sharedX(i));
 
    multsc->setSharedX(i, m_inputWS->sharedX(i));
 
    // Final detector position
    if (!spectrumInfo.hasDetectors(i)) {
      std::ostringstream os;
      os << "No valid detector object found for spectrum at workspace index '" << i << "'. No correction calculated.";
      g_log.information(os.str());
      continue;
    }
 
    // Initialize random number generator
    const int seed = getProperty("Seed");
    CurveFitting::MSVesuvioHelper::RandomVariateGenerator rng(seed + int(i));
 
    // the output spectrum objects have references to where the data will be
    // stored
    calculateMS(rng, i, totalsc->getSpectrum(i), multsc->getSpectrum(i));
    PARALLEL_END_INTERRUPT_REGION
  }
  PARALLEL_CHECK_INTERRUPT_REGION
 
  setProperty("TotalScatteringWS", totalsc);
  setProperty("MultipleScatteringWS", multsc);
}
 
void VesuvioCalculateMS::cacheInputs() {
  // Algorithm
  int nscatters = getProperty("NumScatters");
  m_nscatters = static_cast<size_t>(nscatters);
  int nruns = getProperty("NumRuns");
  m_nruns = static_cast<size_t>(nruns);
  int nevents = getProperty("NumEventsPerRun");
  m_nevents = static_cast<size_t>(nevents);
 
  // -- Geometry --
  const auto instrument = m_inputWS->getInstrument();
  m_beamDir = instrument->getSample()->getPos() - instrument->getSource()->getPos();
  m_beamDir.normalize();
 
  const auto rframe = instrument->getReferenceFrame();
  m_acrossIdx = rframe->pointingHorizontal();
  m_upIdx = rframe->pointingUp();
  m_beamIdx = rframe->pointingAlongBeam();
 
  m_srcR2 = getProperty("BeamRadius");
  // Convert to metres
  m_srcR2 /= 100.0;
 
  // Sample shape
  m_sampleShape = &(m_inputWS->sample().getShape());
  // We know the shape is valid from the property validator
  // Use the bounding box as an approximation to determine the extents
  // as this will match in both height and width for a cuboid & cylinder
  // sample shape
  Geometry::BoundingBox bounds = m_sampleShape->getBoundingBox();
  V3D boxWidth = bounds.width();
  // Use half-width/height for easier calculation later
  m_halfSampleWidth = 0.5 * boxWidth[m_acrossIdx];
  m_halfSampleHeight = 0.5 * boxWidth[m_upIdx];
  m_halfSampleThick = 0.5 * boxWidth[m_beamIdx];
 
  // -- Workspace --
  const auto &inX = m_inputWS->x(0);
  m_tmin = inX.front() * 1e-06;
  m_tmax = inX.back() * 1e-06;
  m_delt = (inX[1] - inX.front());
 
  // -- Sample --
  int nmasses = getProperty("NoOfMasses");
  std::vector<double> sampleInfo = getProperty("AtomicProperties");
  const auto nInputAtomProps = static_cast<int>(sampleInfo.size());
  const int nExptdAtomProp(3);
  if (nInputAtomProps != nExptdAtomProp * nmasses) {
    std::ostringstream os;
    os << "Inconsistent AtomicProperties list defined. Expected " << nExptdAtomProp * nmasses
       << " values, however, only " << sampleInfo.size() << " have been given.";
    throw std::invalid_argument(os.str());
  }
  const int natoms = nInputAtomProps / 3;
  m_sampleProps = std::make_unique<SampleComptonProperties>(natoms);
  m_sampleProps->density = getProperty("SampleDensity");
 
  double totalMass(0.0); // total mass in grams
  m_sampleProps->totalxsec = 0.0;
  for (int i = 0; i < natoms; ++i) {
    auto &comptonAtom = m_sampleProps->atoms[i];
    comptonAtom.mass = sampleInfo[nExptdAtomProp * i];
    totalMass += comptonAtom.mass * PhysicalConstants::AtomicMassUnit * 1000;
 
    const double xsec = sampleInfo[nExptdAtomProp * i + 1];
    comptonAtom.sclength = sqrt(xsec / (4.0 * M_PI));
    const double factor = 1.0 + (PhysicalConstants::NeutronMassAMU / comptonAtom.mass);
    m_sampleProps->totalxsec += (xsec / (factor * factor));
 
    comptonAtom.profile = sampleInfo[nExptdAtomProp * i + 2];
  }
  const double numberDensity = m_sampleProps->density * 1e6 / totalMass; // formula units/m^3
  m_sampleProps->mu = numberDensity * m_sampleProps->totalxsec * 1e-28;
 
  // -- Detector geometry -- choose first detector that is not a monitor
  const auto &spectrumInfo = m_inputWS->spectrumInfo();
  int64_t index = -1;
  for (size_t i = 0; i < m_inputWS->getNumberHistograms(); ++i) {
    if (!spectrumInfo.hasDetectors(i))
      continue;
    if (!spectrumInfo.isMonitor(i)) {
      index = i;
      break;
    }
  }
  // Bounding box in detector frame
  if (index < 0) {
    throw std::runtime_error("Failed to get detector");
  }
  // If is a detector group then take shape of first pixel
  // All detectors in same bansk should be same shape anyway
  const size_t detIndex = spectrumInfo.spectrumDefinition(index)[0].first;
  const auto &componentInfo = m_inputWS->componentInfo();
  if (!componentInfo.hasValidShape(detIndex)) {
    throw std::invalid_argument("Detector pixel has no defined shape!");
  }
  Geometry::BoundingBox detBounds = componentInfo.shape(detIndex).getBoundingBox();
  V3D detBoxWidth = detBounds.width();
  m_detWidth = detBoxWidth[m_acrossIdx];
  m_detHeight = detBoxWidth[m_upIdx];
  m_detThick = detBoxWidth[m_beamIdx];
 
  // Foil resolution
  size_t foilIndex;
  try {
    foilIndex = componentInfo.indexOfAny("foil-pos0");
  } catch (const std::invalid_argument &) {
    throw std::runtime_error("Workspace has no gold foil component defined.");
  }
  const auto foilRes = componentInfo.getNumberParameter(foilIndex, "hwhm_lorentz", false);
  if (foilRes.empty()) {
    throw std::runtime_error("Foil component has no hwhm_lorentz parameter defined.");
  }
  m_foilRes = foilRes.front();
}
 
void VesuvioCalculateMS::calculateMS(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const size_t wsIndex,
                                     API::ISpectrum &totalsc, API::ISpectrum &multsc) const {
  // Detector information
  DetectorParams detpar = ConvertToYSpace::getDetectorParameters(m_inputWS, wsIndex);
  detpar.t0 *= 1e6; // t0 in microseconds here
  Functions::ResolutionParams respar = Functions::VesuvioResolution::getResolutionParameters(m_inputWS, wsIndex);
 
  // Final counts averaged over all simulations
  CurveFitting::MSVesuvioHelper::SimulationAggregator accumulator(m_nruns);
  for (size_t i = 0; i < m_nruns; ++i) {
    m_progress->report("MS calculation: idx=" + std::to_string(wsIndex) + ", run=" + std::to_string(i));
 
    simulate(rng, detpar, respar, accumulator.newSimulation(m_nscatters, m_inputWS->blocksize()));
 
    m_progress->report("MS calculation: idx=" + std::to_string(wsIndex) + ", run=" + std::to_string(i));
  }
 
  // Average over all runs and assign to output workspaces
  CurveFitting::MSVesuvioHelper::SimulationWithErrors avgCounts = accumulator.average();
  avgCounts.normalise();
 
  assignToOutput(avgCounts, totalsc, multsc);
}
 
void VesuvioCalculateMS::simulate(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng,
                                  const DetectorParams &detpar, const ResolutionParams &respar,
                                  CurveFitting::MSVesuvioHelper::Simulation &simulCounts) const {
  for (size_t i = 0; i < m_nevents; ++i) {
    calculateCounts(rng, detpar, respar, simulCounts);
  }
}
 
void VesuvioCalculateMS::assignToOutput(const CurveFitting::MSVesuvioHelper::SimulationWithErrors &avgCounts,
                                        API::ISpectrum &totalsc, API::ISpectrum &multsc) const {
  // Sum up all multiple scatter events
  auto &msscatY = multsc.mutableY();
  auto &msscatE = multsc.mutableE();
  for (size_t i = 1; i < m_nscatters; ++i) //(i >= 1 for multiple scatters)
  {
    const auto &counts = avgCounts.sim.counts[i];
 
    msscatY += counts;
 
    const auto &scerrors = avgCounts.errors[i];
    // sum errors in quadrature
    std::transform(scerrors.begin(), scerrors.end(), msscatE.begin(), msscatE.begin(),
                   VectorHelper::SumGaussError<double>());
  }
  // for total scattering add on single-scatter events
  auto &totalscY = totalsc.mutableY();
  auto &totalscE = totalsc.mutableE();
  const auto &counts0 = avgCounts.sim.counts.front();
  std::transform(counts0.begin(), counts0.end(), msscatY.begin(), totalscY.begin(), std::plus<double>());
  const auto &errors0 = avgCounts.errors.front();
  std::transform(errors0.begin(), errors0.end(), msscatE.begin(), totalscE.begin(),
                 VectorHelper::SumGaussError<double>());
}
 
double VesuvioCalculateMS::calculateCounts(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng,
                                           const DetectorParams &detpar, const ResolutionParams &respar,
                                           CurveFitting::MSVesuvioHelper::Simulation &simulation) const {
  double weightSum(0.0);
 
  // moderator coord in lab frame
  V3D srcPos = generateSrcPos(rng, detpar.l1);
  if (fabs(srcPos[m_acrossIdx]) > m_halfSampleWidth || fabs(srcPos[m_upIdx]) > m_halfSampleHeight) {
    return 0.0; // misses sample
  }
 
  // track various variables during calculation
  std::vector<double> weights(m_nscatters, 1.0), // start at 1.0
      tofs(m_nscatters,
           0.0), // tof accumulates for each piece of the calculation
      en1(m_nscatters, 0.0);
 
  const double vel2 = sqrt(detpar.efixed / MASS_TO_MEV);
  const double t2 = detpar.l2 / vel2;
  en1[0] = generateE0(rng, detpar.l1, t2, weights[0]);
  tofs[0] = generateTOF(rng, en1[0], respar.dtof,
                        respar.dl1); // correction for resolution in l1
 
  // Neutron path
  std::vector<V3D> scatterPts(m_nscatters), // track origin of each scatter
      neutronDirs(m_nscatters);             // neutron directions
  V3D startPos(srcPos);
  neutronDirs[0] = m_beamDir;
 
  generateScatter(rng, startPos, neutronDirs[0], weights[0], scatterPts[0]);
  double distFromStart = startPos.distance(scatterPts[0]);
  // Compute TOF for first scatter event
  const double vel0 = sqrt(en1[0] / MASS_TO_MEV);
  tofs[0] += (distFromStart * 1e6 / vel0);
 
  // multiple scatter events within sample, i.e not including zeroth
  for (size_t i = 1; i < m_nscatters; ++i) {
    weights[i] = weights[i - 1];
    tofs[i] = tofs[i - 1];
 
    // Generate a new direction of travel
    const V3D &prevSc = scatterPts[i - 1];
    V3D &curSc = scatterPts[i];
    const V3D &oldDir = neutronDirs[i - 1];
    V3D &newDir = neutronDirs[i];
    size_t ntries(0);
    do {
      const double randth = acos(2.0 * rng.flat() - 1.0);
      const double randphi = 2.0 * M_PI * rng.flat();
      newDir.azimuth_polar_SNS(1.0, randphi, randth);
 
      // Update weight
      const double wgt = weights[i];
      if (generateScatter(rng, prevSc, newDir, weights[i], curSc))
        break;
      else {
        weights[i] = wgt; // put it back to what it was
        ++ntries;
      }
    } while (ntries < MAX_SCATTER_PT_TRIES);
    if (ntries == MAX_SCATTER_PT_TRIES) {
      throw std::runtime_error("Cannot generate valid trajectory from within "
                               "the sample that intersects the sample. Does it "
                               "have a valid shape?");
    }
 
    const double scang = newDir.angle(oldDir);
    auto e1range = calculateE1Range(scang, en1[i - 1]);
 
    en1[i] = e1range.first + rng.flat() * (e1range.second - e1range.first);
 
    const double d2sig = partialDiffXSec(en1[i - 1], en1[i], scang);
    double weight = d2sig * 4.0 * M_PI * (e1range.second - e1range.first) / m_sampleProps->totalxsec;
    // accumulate total weight
    weightSum += weight;
    weights[i] *= weight; // account for this scatter on top of previous
 
    // Increment time of flight...
    const double veli = sqrt(en1[i] / MASS_TO_MEV);
    tofs[i] += (curSc.distance(prevSc) * 1e6 / veli);
  }
 
  // force all orders in to current detector
  const auto &inX = m_inputWS->x(0);
  for (size_t i = 0; i < m_nscatters; ++i) {
    double scang(0.0), distToExit(0.0);
 
    V3D detPos = generateDetectorPos(rng, detpar.pos, en1[i], scatterPts[i], neutronDirs[i], scang, distToExit);
 
    // Weight by probability neutron leaves sample
    double &curWgt = weights[i];
    curWgt *= exp(-m_sampleProps->mu * distToExit);
    // Weight by cross-section for the final energy
    const double efinal = generateE1(rng, detpar.theta, detpar.efixed, m_foilRes);
    curWgt *= partialDiffXSec(en1[i], efinal, scang) / m_sampleProps->totalxsec;
    // final TOF
    const double veli = sqrt(efinal / MASS_TO_MEV);
    tofs[i] += detpar.t0 + (scatterPts[i].distance(detPos) * 1e6) / veli;
    // "Bin" weight into appropriate place
    std::vector<double> &counts = simulation.counts[i];
    const double finalTOF = tofs[i];
 
    for (size_t it = 0; it < inX.size(); ++it) {
      if (inX[it] - 0.5 * m_delt < finalTOF && finalTOF < inX[it] + 0.5 * m_delt) {
        counts[it] += curWgt;
        break;
      }
    }
  }
 
  return weightSum;
}
 
V3D VesuvioCalculateMS::generateSrcPos(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng,
                                       const double l1) const {
  double radius(-1.0), widthPos(0.0), heightPos(0.0);
  do {
    widthPos = -m_srcR2 + 2.0 * m_srcR2 * rng.flat();
    heightPos = -m_srcR2 + 2.0 * m_srcR2 * rng.flat();
    using std::sqrt;
    radius = sqrt(widthPos * widthPos + heightPos * heightPos);
  } while (radius > m_srcR2);
  // assign to output
  V3D srcPos;
  srcPos[m_acrossIdx] = widthPos;
  srcPos[m_upIdx] = heightPos;
  srcPos[m_beamIdx] = -l1;
  return srcPos;
}
 
double VesuvioCalculateMS::generateE0(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const double l1,
                                      const double t2, double &weight) const {
  const double tof = m_tmin + (m_tmax - m_tmin) * rng.flat();
  const double t1 = (tof - t2);
  const double vel0 = l1 / t1;
  const double en0 = MASS_TO_MEV * vel0 * vel0;
 
  weight = 2.0 * en0 / t1 / pow(en0, 0.9);
  weight *= 1e-4; // Reduce weight to ~1
 
  return en0;
}
 
double VesuvioCalculateMS::generateTOF(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const double en0,
                                       const double dtof, const double dl1) const {
  const double vel1 = sqrt(en0 / MASS_TO_MEV);
  const double dt1 = (dl1 / vel1) * 1e6;
  const double xmin(0.0), xmax(15.0 * dt1);
  double dx = 0.5 * (xmax - xmin);
  // Generate a random y position in th distribution
  const double yv = rng.flat();
 
  double xt(xmin);
  double tof = rng.gaussian(0.0, dtof);
  while (true) {
    xt += dx;
    // Y=1-(0.5*X**2/T0**2+X/T0+1)*EXP(-X/T0)
    double y = 1.0 - (0.5 * xt * xt / (dt1 * dt1) + xt / dt1 + 1) * exp(-xt / dt1);
    if (fabs(y - yv) < 1e-4) {
      tof += xt - 3 * dt1;
      break;
    }
    if (y > yv) {
      dx = -fabs(0.5 * dx);
    } else {
      dx = fabs(0.5 * dx);
    }
  }
  return tof;
}
 
bool VesuvioCalculateMS::generateScatter(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng,
                                         const Kernel::V3D &startPos, const Kernel::V3D &direc, double &weight,
                                         V3D &scatterPt) const {
  Track scatterTrack(startPos, direc);
  if (m_sampleShape->interceptSurface(scatterTrack) != 1) {
    return false;
  }
  // Find distance inside object and compute probability of scattering
  const auto &link = scatterTrack.cbegin();
  double totalObjectDist = link->distInsideObject;
  const double scatterProb = 1.0 - exp(-m_sampleProps->mu * totalObjectDist);
  // Select a random point on the track that is the actual scatter point
  // from the scattering probability distribution
  const double dist = -log(1.0 - rng.flat() * scatterProb) / m_sampleProps->mu;
  const double fraction = dist / totalObjectDist;
  // Scatter point is then entry point + fraction of width in each direction
  scatterPt = link->entryPoint;
  V3D edgeDistances = (link->exitPoint - link->entryPoint);
  scatterPt += edgeDistances * fraction;
  // Update weight
  weight *= scatterProb;
  return true;
}
 
std::pair<double, double> VesuvioCalculateMS::calculateE1Range(const double theta, const double en0) const {
  const double k0 = sqrt(en0 / PhysicalConstants::E_mev_toNeutronWavenumberSq);
  const double sth(sin(theta)), cth(cos(theta));
 
  double e1min(1e10), e1max(-1e10); // large so that anything else is smaller
  const auto &atoms = m_sampleProps->atoms;
  for (const auto &atom : atoms) {
    const double mass = atom.mass;
 
    const double fraction = (cth + sqrt(mass * mass - sth * sth)) / (1.0 + mass);
    const double k1 = fraction * k0;
    const double en1 = PhysicalConstants::E_mev_toNeutronWavenumberSq * k1 * k1;
    const double qr = sqrt(k0 * k0 + k1 * k1 - 2.0 * k0 * k1 * cth);
    const double wr = en0 - en1;
    const double width = PhysicalConstants::E_mev_toNeutronWavenumberSq * atom.profile * qr / mass;
    const double e1a = en0 - wr - 10.0 * width;
    const double e1b = en0 - wr + 10.0 * width;
    if (e1a < e1min)
      e1min = e1a;
    if (e1b > e1max)
      e1max = e1b;
  }
  if (e1min < 0.0)
    e1min = 0.0;
  return std::make_pair(e1min, e1max);
}
 
double VesuvioCalculateMS::partialDiffXSec(const double en0, const double en1, const double theta) const {
  const double rt2pi = sqrt(2.0 * M_PI);
 
  const double k0 = sqrt(en0 / PhysicalConstants::E_mev_toNeutronWavenumberSq);
  const double k1 = sqrt(en1 / PhysicalConstants::E_mev_toNeutronWavenumberSq);
  const double q = sqrt(k0 * k0 + k1 * k1 - 2.0 * k0 * k1 * cos(theta));
  const double w = en0 - en1;
 
  double pdcs(0.0);
  const auto &atoms = m_sampleProps->atoms;
  if (q > 0.0) // avoid continuous checking in loop
  {
    for (const auto &atom : atoms) {
      const double jstddev = atom.profile;
      const double mass = atom.mass;
      const double y = mass * w / (4.18036 * q) - 0.5 * q;
      const double jy = exp(-0.5 * y * y / (jstddev * jstddev)) / (jstddev * rt2pi);
      const double sqw = mass * jy / (4.18036 * q);
 
      const double sclength = atom.sclength;
      pdcs += sclength * sclength * (k1 / k0) * sqw;
    }
  } else {
    for (const auto &atom : atoms) {
      const double sclength = atom.sclength;
      pdcs += sclength * sclength;
    }
  }
 
  return pdcs;
}
 
V3D VesuvioCalculateMS::generateDetectorPos(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng,
                                            const V3D &nominalPos, const double energy, const V3D &scatterPt,
                                            const V3D &direcBeforeSc, double &scang, double &distToExit) const {
  // Inverse attenuation length (m-1) for vesuvio det.
  const double mu = 7430.0 / sqrt(energy);
  // Probability of detection in path thickness.
  const double ps = 1.0 - exp(-mu * m_detThick);
  V3D detPos;
  scang = 0.0;
  distToExit = 0.0;
  size_t ntries(0);
  do {
    // Beam direction by moving to front of "box"define by detector dimensions
    // and then
    // computing expected distance travelled based on probability
    detPos[m_beamIdx] = (nominalPos[m_beamIdx] - 0.5 * m_detThick) - (log(1.0 - rng.flat() * ps) / mu);
    // perturb away from nominal position
    detPos[m_acrossIdx] = nominalPos[m_acrossIdx] + (rng.flat() - 0.5) * m_detWidth;
    detPos[m_upIdx] = nominalPos[m_upIdx] + (rng.flat() - 0.5) * m_detHeight;
 
    // Distance to exit the sample for this order
    const V3D scToDet = normalize(detPos - scatterPt);
    Geometry::Track scatterToDet(scatterPt, scToDet);
    if (m_sampleShape->interceptSurface(scatterToDet) > 0) {
      scang = direcBeforeSc.angle(scToDet);
      const auto &link = scatterToDet.cbegin();
      distToExit = link->distInsideObject;
      break;
    }
    // if point is very close surface then there may be no valid intercept so
    // try again
    ++ntries;
  } while (ntries < MAX_SCATTER_PT_TRIES);
  if (ntries == MAX_SCATTER_PT_TRIES) {
    // Assume it is very close to the surface so that the distance travelled
    // would
    // be a neglible contribution
    distToExit = 0.0;
  }
  return detPos;
}
 
double VesuvioCalculateMS::generateE1(CurveFitting::MSVesuvioHelper::RandomVariateGenerator &rng, const double angle,
                                      const double e1nom, const double e1res) const {
  if (e1res == 0.0)
    return e1nom;
 
  const double randv = rng.flat();
  if (e1nom < 5000.0) {
    if (angle > 90.0)
      return CurveFitting::MSVesuvioHelper::finalEnergyAuDD(randv);
    else
      return CurveFitting::MSVesuvioHelper::finalEnergyAuYap(randv);
  } else {
    return CurveFitting::MSVesuvioHelper::finalEnergyUranium(randv);
  }
}
 
} // namespace Mantid::CurveFitting::Algorithms