Q1D2.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 <utility>
 
#include "MantidAPI/Axis.h"
#include "MantidAPI/CommonBinsValidator.h"
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/ISpectrum.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidAlgorithms/GravitySANSHelper.h"
#include "MantidAlgorithms/Q1D2.h"
#include "MantidAlgorithms/Qhelper.h"
#include "MantidDataObjects/Histogram1D.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidDataObjects/WorkspaceCreation.h"
#include "MantidGeometry/Instrument.h"
#include "MantidGeometry/Instrument/DetectorInfo.h"
#include "MantidIndexing/IndexInfo.h"
#include "MantidKernel/ArrayProperty.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/CompositeValidator.h"
#include "MantidKernel/RebinParamsValidator.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/VectorHelper.h"
#include "MantidTypes/SpectrumDefinition.h"
 
namespace Mantid::Algorithms {
 
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(Q1D2)
 
using namespace Kernel;
using namespace API;
using namespace Geometry;
using namespace DataObjects;
 
Q1D2::Q1D2() : API::Algorithm(), m_dataWS(), m_doSolidAngle(false) {}
 
void Q1D2::init() {
  auto dataVal = std::make_shared<CompositeValidator>();
  dataVal->add<WorkspaceUnitValidator>("Wavelength");
  dataVal->add<HistogramValidator>();
  dataVal->add<InstrumentValidator>();
  dataVal->add<CommonBinsValidator>();
  declareProperty(std::make_unique<WorkspaceProperty<>>("DetBankWorkspace", "", Direction::Input, dataVal),
                  "Particle counts as a function of wavelength");
  declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
                  "Name of the workspace that will contain the result of the calculation");
  declareProperty(std::make_unique<ArrayProperty<double>>("OutputBinning", std::make_shared<RebinParamsValidator>()),
                  "A comma separated list of first bin boundary, width, last bin boundary. "
                  "Optionally\n"
                  "this can be followed by a comma and more widths and last boundary "
                  "pairs.\n"
                  "Negative width values indicate logarithmic binning.");
  declareProperty(std::make_unique<WorkspaceProperty<>>("PixelAdj", "", Direction::Input, PropertyMode::Optional),
                  "Scaling to apply to each spectrum. Must have\n"
                  "the same number of spectra as the DetBankWorkspace");
  auto wavVal = std::make_shared<CompositeValidator>();
  wavVal->add<WorkspaceUnitValidator>("Wavelength");
  wavVal->add<HistogramValidator>();
  declareProperty(
      std::make_unique<WorkspaceProperty<>>("WavelengthAdj", "", Direction::Input, PropertyMode::Optional, wavVal),
      "Scaling to apply to each bin.\n"
      "Must have the same number of bins as the DetBankWorkspace");
  declareProperty(
      std::make_unique<WorkspaceProperty<>>("WavePixelAdj", "", Direction::Input, PropertyMode::Optional, dataVal),
      "Scaling that depends on both pixel and wavelength together.\n"
      "Must have the same number of bins and spectra as the DetBankWorkspace.");
  declareProperty("AccountForGravity", false, "Whether to correct for the effects of gravity");
  declareProperty("SolidAngleWeighting", true, "If true, pixels will be weighted by their solid angle.",
                  Direction::Input);
  auto mustBePositive = std::make_shared<BoundedValidator<double>>();
  mustBePositive->setLower(0.0);
  auto mustBeIntOver2 = std::make_shared<BoundedValidator<int>>();
  mustBeIntOver2->setLower(3);
 
  declareProperty("RadiusCut", 0.0, mustBePositive,
                  "To increase resolution some wavelengths are excluded within "
                  "this distance from the beam center (mm). Note that RadiusCut\n"
                  " and WaveCut both need to be larger than 0 to affect \n"
                  "the effective cutoff. See the algorithm description for\n"
                  " a detailed explanation of the cutoff.");
 
  declareProperty("WaveCut", 0.0, mustBePositive,
                  "To increase resolution by starting to remove some wavelengths below this"
                  " threshold (angstrom).  Note that WaveCut\n"
                  " and RadiusCut both need to be larger than 0 to affect \n"
                  "on the effective cutoff. See the algorithm description for\n"
                  " a detailed explanation of the cutoff.");
 
  declareProperty("OutputParts", false,
                  "Set to true to output two additional workspaces which will "
                  "have the names OutputWorkspace_sumOfCounts "
                  "OutputWorkspace_sumOfNormFactors. The division of "
                  "_sumOfCounts and _sumOfNormFactors equals the workspace"
                  " returned by the property OutputWorkspace "
                  "(default is false).");
  declareProperty("ExtraLength", 0.0, mustBePositive, "Additional length for gravity correction.");
  declareProperty(
      "SolidAngleNumberOfCylinderSlices", 10, mustBeIntOver2,
      "The number of angular slices used when triangulating a cylinder in order to calculate the solid "
      "angle of a tube detector. The default is 10 to preserve legacy behaviour, but increased accuracy has "
      "been observed when using values of 11+");
  declareProperty(
      std::make_unique<WorkspaceProperty<>>("QResolution", "", Direction::Input, PropertyMode::Optional, dataVal),
      "Workspace to calculate the Q resolution.\n");
}
void Q1D2::exec() {
  m_dataWS = getProperty("DetBankWorkspace");
  MatrixWorkspace_const_sptr waveAdj = getProperty("WavelengthAdj");
  MatrixWorkspace_const_sptr pixelAdj = getProperty("PixelAdj");
  MatrixWorkspace_const_sptr wavePixelAdj = getProperty("WavePixelAdj");
  MatrixWorkspace_const_sptr qResolution = getProperty("QResolution");
 
  const bool doGravity = getProperty("AccountForGravity");
  m_doSolidAngle = getProperty("SolidAngleWeighting");
 
  // throws if we don't have common binning or another incompatibility
  Qhelper helper;
  helper.examineInput(m_dataWS, waveAdj, pixelAdj, qResolution);
  // FIXME: how to examine the wavePixelAdj?
  g_log.debug() << "All input workspaces were found to be valid\n";
  // normalization as a function of wavelength (i.e. centers of x-value bins)
  double const *const binNorms = waveAdj ? &(waveAdj->y(0)[0]) : nullptr;
  // error on the wavelength normalization
  double const *const binNormEs = waveAdj ? &(waveAdj->e(0)[0]) : nullptr;
 
  // define the (large number of) data objects that are going to be used in all
  // iterations of the loop below
 
  // Flag to decide if Q Resolution is to be used
  auto useQResolution = static_cast<bool>(qResolution);
 
  // this will become the output workspace from this algorithm
  MatrixWorkspace_sptr outputWS = setUpOutputWorkspace(getProperty("OutputBinning"));
 
  auto &QOut = outputWS->x(0);
  auto &YOut = outputWS->mutableY(0);
  auto &EOutTo2 = outputWS->mutableE(0);
  // normalisation that is applied to counts in each Q bin
  HistogramData::HistogramY normSum(YOut.size(), 0.0);
  // the error on the normalisation
  HistogramData::HistogramE normError2(EOutTo2.size(), 0.0);
  // the averaged Q resolution.
  HistogramData::HistogramDx qResolutionOut(YOut.size(), 0.0);
 
  const auto numSpec = static_cast<int>(m_dataWS->getNumberHistograms());
  Progress progress(this, 0.05, 1.0, numSpec + 1);
 
  const auto &spectrumInfo = m_dataWS->spectrumInfo();
  PARALLEL_FOR_IF(Kernel::threadSafe(*m_dataWS, *outputWS, pixelAdj.get()))
  for (int i = 0; i < numSpec; ++i) {
    PARALLEL_START_INTERRUPT_REGION
    if (!spectrumInfo.hasDetectors(i)) {
      g_log.warning() << "Workspace index " << i << " (SpectrumIndex = " << m_dataWS->getSpectrum(i).getSpectrumNo()
                      << ") has no detector assigned to it - discarding\n";
      continue;
    }
    // Skip if we have a monitor or if the detector is masked.
    if (spectrumInfo.isMonitor(i) || spectrumInfo.isMasked(i))
      continue;
 
    // get the bins that are included inside the RadiusCut/WaveCutcut off, those
    // to calculate for
    // const size_t wavStart = waveLengthCutOff(i);
    const size_t wavStart =
        helper.waveLengthCutOff(m_dataWS, spectrumInfo, getProperty("RadiusCut"), getProperty("WaveCut"), i);
    if (wavStart >= m_dataWS->y(i).size()) {
      // all the spectra in this detector are out of range
      continue;
    }
 
    const size_t numWavbins = m_dataWS->y(i).size() - wavStart;
    // make just one call to new to reduce CPU overhead on each thread, access
    // to these
    // three "arrays" is via iterators
    HistogramData::HistogramY _noDirectUseStorage_(3 * numWavbins);
    // normalization term
    auto norms = _noDirectUseStorage_.begin();
    // the error on these weights, it contributes to the error calculation on
    // the output workspace
    auto normETo2s = norms + numWavbins;
    // the Q values calculated from input wavelength workspace
    auto QIn = normETo2s + numWavbins;
 
    // the weighting for this input spectrum that is added to the normalization
    calculateNormalization(wavStart, i, pixelAdj, wavePixelAdj, binNorms, binNormEs, norms, normETo2s);
 
    // now read the data from the input workspace, calculate Q for each bin
    convertWavetoQ(spectrumInfo, i, doGravity, wavStart, QIn, getProperty("ExtraLength"));
 
    // Pointers to the counts data and it's error
    auto YIn = m_dataWS->y(i).cbegin() + wavStart;
    auto EIn = m_dataWS->e(i).cbegin() + wavStart;
 
    // Pointers to the QResolution data. Note that the xdata was initially the
    // same, hence
    // the same indexing applies to the y values of m_dataWS and qResolution
    // If we want to use it set it to the correct value, else to YIN, although
    // that does not matter, as
    // we won't use it
    auto QResIn = useQResolution ? (qResolution->y(i).cbegin() + wavStart) : YIn;
 
    // when finding the output Q bin remember that the input Q bins (from the
    // convert to wavelength) start high and reduce
    auto loc = QOut.cend();
    // sum the Q contributions from each individual spectrum into the output
    // array
    const auto end = m_dataWS->y(i).cend();
    for (; YIn != end; ++YIn, ++EIn, ++QIn, ++norms, ++normETo2s) {
      // find the output bin that each input y-value will fall into, remembering
      // there is one more bin boundary than bins
      getQBinPlus1(QOut, *QIn, loc);
      // ignore counts that are out of the output range
      if ((loc != QOut.begin()) && (loc != QOut.end())) {
        // the actual Q-bin to add something to
        const size_t bin = loc - QOut.begin() - 1;
        PARALLEL_CRITICAL(q1d_counts_sum) {
          YOut[bin] += *YIn;
          normSum[bin] += *norms;
          // these are the errors squared which will be summed and square rooted
          // at the end
          EOutTo2[bin] += (*EIn) * (*EIn);
          normError2[bin] += *normETo2s;
          if (useQResolution) {
            auto QBin = (QOut[bin + 1] - QOut[bin]);
            // Here we need to take into account the Bin width and the count
            // weigthing. The
            // formula should be YIN* sqrt(QResIn^2 + (QBin/sqrt(12))^2)
            qResolutionOut[bin] += (*YIn) * std::sqrt((*QResIn) * (*QResIn) + QBin * QBin / 12.0);
          }
        }
      }
 
      // Increment the QResolution iterator
      if (useQResolution) {
        ++QResIn;
      }
    }
 
    PARALLEL_CRITICAL(q1d_spectra_map) {
      progress.report("Computing I(Q)");
      // Add up the detector IDs in the output spectrum at workspace index 0
      const auto &inSpec = m_dataWS->getSpectrum(i);
      auto &outSpec = outputWS->getSpectrum(0);
      outSpec.addDetectorIDs(inSpec.getDetectorIDs());
    }
 
    PARALLEL_END_INTERRUPT_REGION
  }
  PARALLEL_CHECK_INTERRUPT_REGION
 
  if (useQResolution) {
    // The number of Q (x)_ values is N, while the number of DeltaQ values is
    // N-1,
    // Richard Heenan suggested to duplicate the last entry of DeltaQ
    auto countsIterator = YOut.cbegin();
    auto qResolutionIterator = qResolutionOut.begin();
    for (; countsIterator != YOut.end(); ++countsIterator, ++qResolutionIterator) {
      // Divide by the counts of the Qbin, if the counts are 0, the the
      // qresolution will also be 0
      if ((*countsIterator) > 0.0) {
        *qResolutionIterator = (*qResolutionIterator) / (*countsIterator);
      }
    }
    outputWS->setPointStandardDeviations(0, std::move(qResolutionOut));
  }
 
  bool doOutputParts = getProperty("OutputParts");
  if (doOutputParts) {
    MatrixWorkspace_sptr ws_sumOfCounts = WorkspaceFactory::Instance().create(outputWS);
    ws_sumOfCounts->setSharedX(0, outputWS->sharedX(0));
    // Copy now as YOut is modified in normalize
    ws_sumOfCounts->mutableY(0) = YOut;
    ws_sumOfCounts->setSharedDx(0, outputWS->sharedDx(0));
    ws_sumOfCounts->setFrequencyVariances(0, outputWS->e(0));
 
    MatrixWorkspace_sptr ws_sumOfNormFactors = WorkspaceFactory::Instance().create(outputWS);
    ws_sumOfNormFactors->setSharedX(0, outputWS->sharedX(0));
    ws_sumOfNormFactors->mutableY(0) = normSum;
    ws_sumOfNormFactors->setSharedDx(0, outputWS->sharedDx(0));
    ws_sumOfNormFactors->setFrequencyVariances(0, normError2);
 
    helper.outputParts(this, ws_sumOfCounts, ws_sumOfNormFactors);
  }
 
  progress.report("Normalizing I(Q)");
  // finally divide the number of counts in each output Q bin by its weighting
  normalize(normSum, normError2, YOut, EOutTo2);
 
  setProperty("OutputWorkspace", outputWS);
}
 
API::MatrixWorkspace_sptr Q1D2::setUpOutputWorkspace(const std::vector<double> &binParams) const {
  // Calculate the output binning
  std::vector<double> xAxisTmp;
  VectorHelper::createAxisFromRebinParams(binParams, xAxisTmp);
  HistogramData::BinEdges XOut(std::move(xAxisTmp));
 
  // Create output workspace. On all but rank 0 this is a temporary workspace.
  Indexing::IndexInfo indexInfo(1);
  indexInfo.setSpectrumDefinitions(std::vector<SpectrumDefinition>(1));
  auto outputWS = create<MatrixWorkspace>(*m_dataWS, indexInfo, XOut);
  outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("MomentumTransfer");
  outputWS->setYUnitLabel("1/cm");
 
  outputWS->setDistribution(true);
 
  outputWS->getSpectrum(0).clearDetectorIDs();
  outputWS->getSpectrum(0).setSpectrumNo(1);
 
  return outputWS;
}
 
void Q1D2::calculateNormalization(const size_t wavStart, const size_t wsIndex,
                                  const API::MatrixWorkspace_const_sptr &pixelAdj,
                                  const API::MatrixWorkspace_const_sptr &wavePixelAdj, double const *const binNorms,
                                  double const *const binNormEs, HistogramData::HistogramY::iterator norm,
                                  HistogramData::HistogramY::iterator normETo2) const {
  double detectorAdj, detAdjErr;
  pixelWeight(pixelAdj, wsIndex, detectorAdj, detAdjErr);
  // use that the normalization array ends at the start of the error array
  for (auto n = norm, e = normETo2; n != normETo2; ++n, ++e) {
    *n = detectorAdj;
    *e = detAdjErr * detAdjErr;
  }
 
  if (binNorms && binNormEs) {
    if (wavePixelAdj)
      // pass the iterator for the wave pixel Adj dependent
      addWaveAdj(binNorms + wavStart, binNormEs + wavStart, norm, normETo2, wavePixelAdj->y(wsIndex).begin() + wavStart,
                 wavePixelAdj->e(wsIndex).begin() + wavStart);
    else
      addWaveAdj(binNorms + wavStart, binNormEs + wavStart, norm, normETo2);
  }
  normToMask(wavStart, wsIndex, norm, normETo2);
}
 
void Q1D2::pixelWeight(const API::MatrixWorkspace_const_sptr &pixelAdj, const size_t wsIndex, double &weight,
                       double &error) const {
  const auto &detectorInfo = m_dataWS->detectorInfo();
  const V3D samplePos = detectorInfo.samplePosition();
 
  if (m_doSolidAngle) {
    const int numberOfCylinderSlices = getProperty("SolidAngleNumberOfCylinderSlices");
    weight = 0.0;
    for (const auto detID : m_dataWS->getSpectrum(wsIndex).getDetectorIDs()) {
      const auto index = detectorInfo.indexOf(detID);
      if (!detectorInfo.isMasked(index))
        weight += detectorInfo.detector(index).solidAngle({samplePos, numberOfCylinderSlices});
    }
  } else
    weight = 1.0;
 
  if (weight < 1e-200) {
    throw std::logic_error("Invalid (zero or negative) solid angle for one detector");
  }
  // this input multiplies up the adjustment if it exists
  if (pixelAdj) {
    weight *= pixelAdj->readY(wsIndex)[0];
    error = weight * pixelAdj->readE(wsIndex)[0];
  } else {
    error = 0.0;
  }
}
 
void Q1D2::addWaveAdj(const double *c, const double *Dc, HistogramData::HistogramY::iterator bInOut,
                      HistogramData::HistogramY::iterator e2InOut) const {
  // normalize by the wavelength dependent correction, keeping the percentage
  // errors the same
  // the error when a = b*c, the formula for Da, the error on a, in terms of Db,
  // etc. is (Da/a)^2 = (Db/b)^2 + (Dc/c)^2
  //(Da)^2 = ((Db*a/b)^2 + (Dc*a/c)^2) = (Db*c)^2 + (Dc*b)^2
  // the variable names relate to those above as: existing values (b=bInOut)
  // multiplied by the additional errors (Dc=binNormEs), existing errors
  // (Db=sqrt(e2InOut)) times new factor (c=binNorms)
 
  // use the fact that error array follows straight after the normalization
  // array
  const auto end = e2InOut;
  for (; bInOut != end; ++e2InOut, ++c, ++Dc, ++bInOut) {
    // first the error
    *e2InOut = ((*e2InOut) * (*c) * (*c)) + ((*Dc) * (*Dc) * (*bInOut) * (*bInOut));
    // now the actual calculation a = b*c
    *bInOut = (*bInOut) * (*c);
  }
}
void Q1D2::addWaveAdj(const double *c, const double *Dc, HistogramData::HistogramY::iterator bInOut,
                      HistogramData::HistogramY::iterator e2InOut,
                      HistogramData::HistogramY::const_iterator wavePixelAdjData,
                      HistogramData::HistogramE::const_iterator wavePixelAdjError) const {
  // normalize by the wavelength dependent correction, keeping the percentage
  // errors the same
  // the error when a = b*c*e, the formula for Da, the error on a, in terms of
  // Db, etc. is
  // (Da/a)^2 = (Db/b)^2 + (Dc/c)^2 + (De/e)^2
  //(Da)^2 = ((Db*a/b)^2 + (Dc*a/c)^2) + (De * a/e)^2
  // But: a/b = c*e; a/c = b*e; a/e = b*c;
  // So:
  // (Da)^2 = (c*e*Db)^2 + (b*e*Dc)^2 + (b*c*De)^2
  //
  // Consider:
  // Da = Error (e2InOut)
  // Db^2 = PixelDependentError (e2InOut)
  // b  = PixelDependentValue (bInOut)
  // c  = WaveDependentValue (c)
  // Dc = WaveDependentError (Dc)
  // e  = PixelWaveDependentValue (wavePixelAdjData)
  // De = PiexlWaveDependentError (wavePixelAdjError)
 
  // use the fact that error array follows straight after the normalization
  // array
  const auto end = e2InOut;
  for (; bInOut != end; ++e2InOut, ++c, ++Dc, ++bInOut, ++wavePixelAdjData, ++wavePixelAdjError) {
    // first the error
    *e2InOut = ((*e2InOut) * (*c) * (*c) * (*wavePixelAdjData) * (*wavePixelAdjData)) +
               ((*Dc) * (*Dc) * (*bInOut) * (*bInOut) * (*wavePixelAdjData) * (*wavePixelAdjData)) +
               ((*wavePixelAdjError) * (*wavePixelAdjError) * (*c) * (*c) * (*bInOut) * (*bInOut));
    // now the actual calculation a = b*c*e : Pixel * Wave * PixelWave
    *bInOut = (*bInOut) * (*c) * (*wavePixelAdjData);
  }
}
 
void Q1D2::normToMask(const size_t offSet, const size_t wsIndex, const HistogramData::HistogramY::iterator theNorms,
                      const HistogramData::HistogramY::iterator errorSquared) const {
  // if any bins are masked it is normally a small proportion
  if (m_dataWS->hasMaskedBins(wsIndex)) {
    // Get a reference to the list of masked bins
    const MatrixWorkspace::MaskList &mask = m_dataWS->maskedBins(wsIndex);
    // Now iterate over the list, adjusting the weights for the affected bins
    MatrixWorkspace::MaskList::const_iterator it;
    for (it = mask.begin(); it != mask.end(); ++it) {
      size_t outBin = it->first;
      if (outBin < offSet) {
        // this masked bin wasn't in the range being delt with anyway
        continue;
      }
      outBin -= offSet;
      // The weight for this masked bin is 1 - the degree to which this bin is
      // masked
      const double factor = 1.0 - (it->second);
      *(theNorms + outBin) *= factor;
      *(errorSquared + outBin) *= factor * factor;
    }
  }
}
 
void Q1D2::convertWavetoQ(const SpectrumInfo &spectrumInfo, const size_t wsInd, const bool doGravity,
                          const size_t offset, HistogramData::HistogramY::iterator Qs, const double extraLength) const {
  static const double FOUR_PI = 4.0 * M_PI;
 
  // wavelengths (lamda) to be converted to Q
  auto waves = m_dataWS->x(wsInd).cbegin() + offset;
  // going from bin boundaries to bin centered x-values the size goes down one
  const auto end = m_dataWS->x(wsInd).end() - 1;
  if (doGravity) {
    GravitySANSHelper grav(spectrumInfo, wsInd, extraLength);
    for (; waves != end; ++Qs, ++waves) {
      // the HistogramValidator at the start should ensure that we have one more
      // bin on the input wavelengths
      const double lambda = 0.5 * (*(waves + 1) + (*waves));
      // as the fall under gravity is wavelength dependent sin theta is now
      // different for each bin with each detector
      const double sinTheta = grav.calcSinTheta(lambda);
      // Now we're ready to go to Q
      *Qs = FOUR_PI * sinTheta / lambda;
    }
  } else {
    // Calculate the Q values for the current spectrum, using Q =
    // 4*pi*sin(theta)/lambda
    const double factor = 2.0 * FOUR_PI * sin(spectrumInfo.twoTheta(wsInd) * 0.5);
    for (; waves != end; ++Qs, ++waves) {
      // the HistogramValidator at the start should ensure that we have one more
      // bin on the input wavelengths
      *Qs = factor / (*(waves + 1) + (*waves));
    }
  }
}
void Q1D2::getQBinPlus1(const HistogramData::HistogramX &OutQs, const double QToFind,
                        HistogramData::HistogramX::const_iterator &loc) const {
  if (loc != OutQs.end()) {
    while (loc != OutQs.begin()) {
      if ((QToFind >= *(loc - 1)) && (QToFind < *loc)) {
        return;
      }
      --loc;
    }
    if (QToFind < *loc) {
      // QToFind is outside the array leave loc == OutQs.begin()
      return;
    }
  } else // loc == OutQs.end()
  {
    if (OutQs.empty() || QToFind > *(loc - 1)) {
      // outside the array leave loc == OutQs.end()
      return;
    }
  }
 
  // we are lost, normally the order of the Q values means we only get here on
  // the first iteration. It's slow
  loc = std::lower_bound(OutQs.begin(), OutQs.end(), QToFind);
}
 
void Q1D2::normalize(const HistogramData::HistogramY &normSum, const HistogramData::HistogramE &normError2,
                     HistogramData::HistogramY &counts, HistogramData::HistogramE &errors) const {
  for (size_t k = 0; k < counts.size(); ++k) {
    // the normalisation is a = b/c where b = counts c =normalistion term
    const double c = normSum[k];
    const double a = counts[k] /= c;
    // when a = b/c, the formula for Da, the error on a, in terms of Db, etc. is
    // (Da/a)^2 = (Db/b)^2 + (Dc/c)^2
    //(Da)^2 = ((Db/b)^2 + (Dc/c)^2)*(b^2/c^2) = ((Db/c)^2 + (b*Dc/c^2)^2) =
    //(Db^2 + (b*Dc/c)^2)/c^2 = (Db^2 + (Dc*a)^2)/c^2
    // this will work as long as c>0, but then the above formula above can't
    // deal with 0 either
    const double aOverc = a / c;
    errors[k] = std::sqrt(errors[k] / (c * c) + normError2[k] * aOverc * aOverc);
  }
}
} // namespace Mantid::Algorithms