Q1D2.cpp

//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/Q1D2.h"
#include "MantidAlgorithms/Qhelper.h"
#include "MantidAPI/CommonBinsValidator.h"
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/ISpectrum.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
#include "MantidDataObjects/Histogram1D.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"

namespace Mantid {
namespace Algorithms {

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(Q1D2)

using namespace Kernel;
using namespace API;
using namespace Geometry;

Q1D2::Q1D2() : API::Algorithm(), m_dataWS(), m_doSolidAngle(false) {}

void Q1D2::init() {
  auto dataVal = boost::make_shared<CompositeValidator>();
  dataVal->add<WorkspaceUnitValidator>("Wavelength");
  dataVal->add<HistogramValidator>();
  dataVal->add<InstrumentValidator>();
  dataVal->add<CommonBinsValidator>();
  declareProperty(new WorkspaceProperty<>("DetBankWorkspace", "",
                                          Direction::Input, dataVal),
                  "Particle counts as a function of wavelength");
  declareProperty(
      new WorkspaceProperty<>("OutputWorkspace", "", Direction::Output),
      "Name of the workspace that will contain the result of the calculation");
  declareProperty(
      new ArrayProperty<double>("OutputBinning",
                                boost::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(new 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 = boost::make_shared<CompositeValidator>();
  wavVal->add<WorkspaceUnitValidator>("Wavelength");
  wavVal->add<HistogramValidator>();
  declareProperty(new 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(
      new 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 = boost::make_shared<BoundedValidator<double>>();
  mustBePositive->setLower(0.0);
  declareProperty("RadiusCut", 0.0, mustBePositive,
                  "To increase resolution some wavelengths are excluded within "
                  "this distance from\n"
                  "the beam center (mm)");
  declareProperty(
      "WaveCut", 0.0, mustBePositive,
      "To increase resolution by starting to remove some wavelengths below this"
      "freshold (angstrom)");
  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(new WorkspaceProperty<>("QResolution", "", Direction::Input,
                                          PropertyMode::Optional, dataVal),
                  "Workspace to calculate the Q resolution.\n");
}
/**
  @ throw invalid_argument if the workspaces are not mututially compatible
*/
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->readY(0)[0]) : NULL;
  // error on the wavelength normalization
  double const *const binNormEs = waveAdj ? &(waveAdj->readE(0)[0]) : NULL;

  // 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 = qResolution ? true : false;

  // this will become the output workspace from this algorithm
  MatrixWorkspace_sptr outputWS =
      setUpOutputWorkspace(getProperty("OutputBinning"));

  const MantidVec &QOut = outputWS->readX(0);
  MantidVec &YOut = outputWS->dataY(0);
  MantidVec &EOutTo2 = outputWS->dataE(0);
  // normalisation that is applied to counts in each Q bin
  MantidVec normSum(YOut.size(), 0.0);
  // the error on the normalisation
  MantidVec normError2(YOut.size(), 0.0);

  // the averaged Q resolution. We need the a named dummy variable although it
  // won't be
  // used since we only want to create a reference to DX if it is really
  // required. Referencing
  // DX sets a flag which might not be desirable.
  MantidVec dummy;
  MantidVec &qResolutionOut =
      useQResolution ? outputWS->dataDx(0) : outputWS->dataY(0);

  const int numSpec = static_cast<int>(m_dataWS->getNumberHistograms());
  Progress progress(this, 0.05, 1.0, numSpec + 1);

  PARALLEL_FOR3(m_dataWS, outputWS, pixelAdj)
  for (int i = 0; i < numSpec; ++i) {
    PARALLEL_START_INTERUPT_REGION
    // Get the pixel relating to this spectrum
    IDetector_const_sptr det;
    try {
      det = m_dataWS->getDetector(i);
    } catch (Exception::NotFoundError &) {
      g_log.warning() << "Workspace index " << i << " (SpectrumIndex = "
                      << m_dataWS->getSpectrum(i)->getSpectrumNo()
                      << ") has no detector assigned to it - discarding"
                      << std::endl;
      // Catch if no detector. Next line tests whether this happened - test
      // placed
      // outside here because Mac Intel compiler doesn't like 'continue' in a
      // catch
      // in an openmp block.
    }
    // If no detector found or if detector is masked shouldn't be included skip
    // onto the next spectrum
    if (!det || det->isMonitor() || det->isMasked()) {
      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, getProperty("RadiusCut"), getProperty("WaveCut"), i);
    if (wavStart >= m_dataWS->readY(i).size()) {
      // all the spectra in this detector are out of range
      continue;
    }

    const size_t numWavbins = m_dataWS->readY(i).size() - wavStart;
    // make just one call to new to reduce CPU overhead on each thread, access
    // to these
    // three "arrays" is via iterators
    MantidVec _noDirectUseStorage_(3 * numWavbins);
    // normalization term
    MantidVec::iterator norms = _noDirectUseStorage_.begin();
    // the error on these weights, it contributes to the error calculation on
    // the output workspace
    MantidVec::iterator normETo2s = norms + numWavbins;
    // the Q values calculated from input wavelength workspace
    MantidVec::iterator 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(i, doGravity, wavStart, QIn, getProperty("ExtraLength"));

    // Pointers to the counts data and it's error
    MantidVec::const_iterator YIn = m_dataWS->readY(i).begin() + wavStart;
    MantidVec::const_iterator EIn = m_dataWS->readE(i).begin() + 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
    MantidVec::const_iterator QResIn =
        useQResolution ? (qResolution->readY(i).begin() + wavStart) : YIn;

    // when finding the output Q bin remember that the input Q bins (from the
    // convert to wavelength) start high and reduce
    MantidVec::const_iterator loc = QOut.end();
    // sum the Q contributions from each individual spectrum into the output
    // array
    const MantidVec::const_iterator end = m_dataWS->readY(i).end();
    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) {
            qResolutionOut[bin] += (*YIn) * (*QResIn);
          }
        }
      }

      // 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 ISpectrum *inSpec = m_dataWS->getSpectrum(i);
      ISpectrum *outSpec = outputWS->getSpectrum(0);
      outSpec->addDetectorIDs(inSpec->getDetectorIDs());
    }

    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_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
    Mantid::MantidVec::const_iterator countsIterator = YOut.begin();
    Mantid::MantidVec::iterator 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);
      }
    }
    // Now duplicate write the second to last element into the last element of
    // the deltaQ vector
    if (qResolutionOut.size() > 1) {
      qResolutionOut.rbegin()[0] = qResolutionOut.rbegin()[1];
    }
  }

  bool doOutputParts = getProperty("OutputParts");
  if (doOutputParts) {
    MatrixWorkspace_sptr ws_sumOfCounts =
        WorkspaceFactory::Instance().create(outputWS);
    ws_sumOfCounts->dataX(0) = outputWS->dataX(0);
    ws_sumOfCounts->dataY(0) = outputWS->dataY(0);
    if (useQResolution) {
      ws_sumOfCounts->dataDx(0) = outputWS->dataDx(0);
    }
    for (size_t i = 0; i < outputWS->dataE(0).size(); i++) {
      ws_sumOfCounts->dataE(0)[i] = sqrt(outputWS->dataE(0)[i]);
    }

    MatrixWorkspace_sptr ws_sumOfNormFactors =
        WorkspaceFactory::Instance().create(outputWS);
    ws_sumOfNormFactors->dataX(0) = outputWS->dataX(0);
    if (useQResolution) {
      ws_sumOfNormFactors->dataDx(0) = outputWS->dataDx(0);
    }
    for (size_t i = 0; i < ws_sumOfNormFactors->dataY(0).size(); i++) {
      ws_sumOfNormFactors->dataY(0)[i] = normSum[i];
      ws_sumOfNormFactors->dataE(0)[i] = sqrt(normError2[i]);
    }

    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);
}

/** Creates the output workspace, its size, units, etc.
*  @param binParams the bin boundary specification using the same same syntax as
* param the Rebin algorithm
*  @return A pointer to the newly-created workspace
*/
API::MatrixWorkspace_sptr
Q1D2::setUpOutputWorkspace(const std::vector<double> &binParams) const {
  // Calculate the output binning
  MantidVecPtr XOut;
  size_t sizeOut = static_cast<size_t>(
      VectorHelper::createAxisFromRebinParams(binParams, XOut.access()));

  // Now create the output workspace
  MatrixWorkspace_sptr outputWS =
      WorkspaceFactory::Instance().create(m_dataWS, 1, sizeOut, sizeOut - 1);
  outputWS->getAxis(0)->unit() =
      UnitFactory::Instance().create("MomentumTransfer");
  outputWS->setYUnitLabel("1/cm");

  // Set the X vector for the output workspace
  outputWS->setX(0, XOut);
  outputWS->isDistribution(true);

  outputWS->getSpectrum(0)->clearDetectorIDs();
  outputWS->getSpectrum(0)->setSpectrumNo(1);

  return outputWS;
}

/** Calculate the normalization term for each output bin
*  @param wavStart [in] the index number of the first bin in the input
* wavelengths that is actually being used
*  @param specInd [in] the spectrum to calculate
*  @param pixelAdj [in] if not NULL this is workspace contains single bins with
* the adjustments, e.g. detector efficencies, for the given spectrum index
*  @param wavePixelAdj [in] if not NULL this is workspace that contains the
* adjustments for the pixels and wavelenght dependend values.
*  @param binNorms [in] pointer to a contigious array of doubles that are the
* wavelength correction from waveAdj workspace, can be NULL
*  @param binNormEs [in] pointer to a contigious array of doubles which
* corrospond to the corrections and are their errors, can be NULL
*  @param norm [out] normalization for each bin, including soild angle, pixel
* correction, the proportion that is not masked and the normalization workspace
*  @param normETo2 [out] this pointer must point to the end of the norm array,
* it will be filled with the total of the error on the normalization
*/
void Q1D2::calculateNormalization(const size_t wavStart, const size_t specInd,
                                  API::MatrixWorkspace_const_sptr pixelAdj,
                                  API::MatrixWorkspace_const_sptr wavePixelAdj,
                                  double const *const binNorms,
                                  double const *const binNormEs,
                                  const MantidVec::iterator norm,
                                  const MantidVec::iterator normETo2) const {
  double detectorAdj, detAdjErr;
  pixelWeight(pixelAdj, specInd, detectorAdj, detAdjErr);
  // use that the normalization array ends at the start of the error array
  for (MantidVec::iterator 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->readY(specInd).begin() + wavStart,
                 wavePixelAdj->readE(specInd).begin() + wavStart);
    else
      addWaveAdj(binNorms + wavStart, binNormEs + wavStart, norm, normETo2);
  }
  normToMask(wavStart, specInd, norm, normETo2);
}

/** Calculates the normalisation for the spectrum specified by the index number
* that was passed
*  as the solid anlge multiplied by the pixelAdj that was passed
*  @param[in] pixelAdj if not NULL this is workspace contains single bins with
* the adjustments, e.g. detector efficencies, for the given spectrum index
*  @param[in] specIndex the spectrum index to return the data from
*  @param[out] weight the solid angle or if pixelAdj the solid anlge times the
* pixel adjustment for this spectrum
*  @param[out] error the error on the weight, only non-zero if pixelAdj
*  @throw LogicError if the solid angle is tiny or negative
*/
void Q1D2::pixelWeight(API::MatrixWorkspace_const_sptr pixelAdj,
                       const size_t specIndex, double &weight,
                       double &error) const {
  const V3D samplePos = m_dataWS->getInstrument()->getSample()->getPos();

  if (m_doSolidAngle)
    weight = m_dataWS->getDetector(specIndex)->solidAngle(samplePos);
  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(specIndex)[0];
    error = weight * pixelAdj->readE(specIndex)[0];
  } else {
    error = 0.0;
  }
}

/** Calculates the contribution to the normalization terms from each bin in a
* spectrum
*  @param[in] c pointer to the start of a contigious array of wavelength
* dependent normalization terms
*  @param[in] Dc pointer to the start of a contigious array that corrosponds to
* wavelength dependent term, having its error
*  @param[in,out] bInOut normalization for each bin, this method multiplise this
* by the proportion that is not masked and the normalization workspace
*  @param[in, out] e2InOut this array must follow straight after the
* normalization array and will contain the error on the normalisation term
* before the WavelengthAdj term
*/
void Q1D2::addWaveAdj(const double *c, const double *Dc,
                      MantidVec::iterator bInOut,
                      MantidVec::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 MantidVec::const_iterator 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);
  }
}
/** Calculates the contribution to the normalization terms from each bin in a
* spectrum
*  @param[in] c pointer to the start of a contigious array of wavelength
* dependent normalization terms
*  @param[in] Dc pointer to the start of a contigious array that corrosponds to
* wavelength dependent term, having its error
*  @param[in,out] bInOut normalization for each bin, this method multiplise this
* by the proportion that is not masked and the normalization workspace
*  @param[in, out] e2InOut this array must follow straight after the
* normalization array and will contain the error on the normalisation term
* before the WavelengthAdj term
* @param[in] wavePixelAdjData normalization correction for each bin for each
* detector pixel.
* @param[in] wavePixelAdjError normalization correction incertainty for each bin
* for each detector pixel.
*/
void Q1D2::addWaveAdj(const double *c, const double *Dc,
                      MantidVec::iterator bInOut, MantidVec::iterator e2InOut,
                      MantidVec::const_iterator wavePixelAdjData,
                      MantidVec::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 MantidVec::const_iterator 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);
  }
}

/** Scaled to bin masking, to the normalization
*  @param[in] offSet the inex number of the first bin in the input wavelengths
* that is actually being used
*  @param[in] specIndex the spectrum to calculate
*  @param[in,out] theNorms normalization for each bin, this is multiplied by the
* proportion that is not masked and the normalization workspace
*  @param[in,out] errorSquared the running total of the square of the
* uncertainty in the normalization
*/
void Q1D2::normToMask(const size_t offSet, const size_t specIndex,
                      const MantidVec::iterator theNorms,
                      const MantidVec::iterator errorSquared) const {
  // if any bins are masked it is normally a small proportion
  if (m_dataWS->hasMaskedBins(specIndex)) {
    // Get a reference to the list of masked bins
    const MatrixWorkspace::MaskList &mask = m_dataWS->maskedBins(specIndex);
    // 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;
    }
  }
}

/** Fills a vector with the Q values calculated from the wavelength bin centers
* from the input workspace and
*  the workspace geometry as Q = 4*pi*sin(theta)/lambda
*  @param[in] specInd the spectrum to calculate
*  @param[in] doGravity if to include gravity in the calculation of Q
*  @param[in] offset index number of the first input bin to use
*  @param[in] extraLength for gravitational correction
*  @param[out] Qs points to a preallocated array that is large enough to contain
* all the calculated Q values
*  @throw NotFoundError if the detector associated with the spectrum is not
* found in the instrument definition
*/
void Q1D2::convertWavetoQ(const size_t specInd, const bool doGravity,
                          const size_t offset, MantidVec::iterator Qs,
                          const double extraLength) const {
  static const double FOUR_PI = 4.0 * M_PI;

  IDetector_const_sptr det = m_dataWS->getDetector(specInd);

  // wavelengths (lamda) to be converted to Q
  MantidVec::const_iterator waves = m_dataWS->readX(specInd).begin() + offset;
  // going from bin boundaries to bin centered x-values the size goes down one
  const MantidVec::const_iterator end = m_dataWS->readX(specInd).end() - 1;
  if (doGravity) {
    GravitySANSHelper grav(m_dataWS, det, 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(m_dataWS->detectorTwoTheta(det) / 2.0);
    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));
    }
  }
}
/** This is a slightly "clever" method as it makes some guesses about where is
* best
*  to look for the right Q bin based on the fact that the input Qs (calcualted
* from wavelengths) tend
*  to go down while the output Qs are always in accending order
*  @param[in] OutQs the array of output Q bin boundaries, this finds the bin
* that contains the QIn value
*  @param[in] QToFind the Q value to find the correct bin for
*  @param[in, out] loc points to the bin boundary (in the OutQs array) whos Q is
* higher than QToFind and higher by the smallest amount. Algorithm starts by
* checking the value of loc passed and then all the bins _downwards_ through the
* array
*/
void Q1D2::getQBinPlus1(const MantidVec &OutQs, const double QToFind,
                        MantidVec::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);
}

/** Divides the number of counts in each output Q bin by the wrighting ("number
* that would expected to arrive")
*  The errors are propogated using the uncorrolated error estimate for
* multiplication/division
*  @param[in] normSum the weighting for each bin
*  @param[in] normError2 square of the error on the normalization
*  @param[in, out] counts counts in each bin
*  @param[in, out] errors input the _square_ of the error on each bin, output
* the total error (unsquared)
*/
void Q1D2::normalize(const MantidVec &normSum, const MantidVec &normError2,
                     MantidVec &counts, MantidVec &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 Algorithms
} // namespace Mantid