C++ MatrixWorkspace_sptr::getSpectrum方法代码示例

/***
 * This will add the maximum spectrum number from ws1 to the data coming from ws2.
 *
 * @param ws1 The first workspace supplied to the algorithm.
 * @param ws2 The second workspace supplied to the algorithm.
 * @param output The workspace that is going to be returned by the algorithm.
 */
void ConjoinWorkspaces::fixSpectrumNumbers(API::MatrixWorkspace_const_sptr ws1, API::MatrixWorkspace_const_sptr /*ws2*/,
                                      API::MatrixWorkspace_sptr output)
{
  // make sure we should bother. If you don't check for overlap, you don't need to
  if (this->getProperty("CheckOverlapping"))
    return;

  // is everything possibly ok?
  specid_t min;
  specid_t max;
  getMinMax(output, min, max);
  if (max - min >= static_cast<specid_t>(output->getNumberHistograms())) // nothing to do then
    return;

  // information for remapping the spectra numbers
  specid_t ws1min;
  specid_t ws1max;
  getMinMax(ws1, ws1min, ws1max);

  // change the axis by adding the maximum existing spectrum number to the current value
  for (size_t i = ws1->getNumberHistograms(); i < output->getNumberHistograms(); i++)
  {
    specid_t origid;
    origid = output->getSpectrum(i)->getSpectrumNo();
    output->getSpectrum(i)->setSpectrumNo(origid + ws1max);
  }
  // To be deprecated:
  output->generateSpectraMap();
}

/** Calculates the normalization constant for the exponential decay
* @param ws :: [input] Workspace containing the spectra to remove exponential
* from
* @return :: Vector containing the normalization constants, N0, for each
* spectrum
*/
std::vector<double>
PhaseQuadMuon::getExponentialDecay(const API::MatrixWorkspace_sptr &ws) {

  size_t nspec = ws->getNumberHistograms();
  size_t npoints = ws->blocksize();

  // Muon life time in microseconds
  double muLife = PhysicalConstants::MuonLifetime * 1e6;

  std::vector<double> n0(nspec, 0.);

  for (size_t h = 0; h < ws->getNumberHistograms(); h++) {

    const auto &X = ws->getSpectrum(h).x();
    const auto &Y = ws->getSpectrum(h).y();
    const auto &E = ws->getSpectrum(h).e();

    double s, sx, sy;
    s = sx = sy = 0;
    for (size_t i = 0; i < npoints; i++) {

      if (Y[i] > 0) {
        double sig = E[i] * E[i] / Y[i] / Y[i];
        s += 1. / sig;
        sx += (X[i] - X[0]) / sig;
        sy += log(Y[i]) / sig;
      }
    }
    n0[h] = exp((sy + sx / muLife) / s);
  }

  return n0;
}

/** Add workspace2 to workspace1 by adding spectrum.
  */
MatrixWorkspace_sptr
AlignAndFocusPowder::conjoinWorkspaces(API::MatrixWorkspace_sptr ws1,
                                       API::MatrixWorkspace_sptr ws2,
                                       size_t offset) {
  // Get information from ws1: maximum spectrum number, and store original
  // spectrum Nos
  size_t nspec1 = ws1->getNumberHistograms();
  specnum_t maxspecNo1 = 0;
  std::vector<specnum_t> origspecNos;
  for (size_t i = 0; i < nspec1; ++i) {
    specnum_t tmpspecNo = ws1->getSpectrum(i).getSpectrumNo();
    origspecNos.push_back(tmpspecNo);
    if (tmpspecNo > maxspecNo1)
      maxspecNo1 = tmpspecNo;
  }

  g_log.information() << "[DBx536] Max spectrum number of ws1 = " << maxspecNo1
                      << ", Offset = " << offset << ".\n";

  size_t nspec2 = ws2->getNumberHistograms();

  // Conjoin 2 workspaces
  Algorithm_sptr alg = this->createChildAlgorithm("AppendSpectra");
  alg->initialize();
  ;

  alg->setProperty("InputWorkspace1", ws1);
  alg->setProperty("InputWorkspace2", ws2);
  alg->setProperty("OutputWorkspace", ws1);
  alg->setProperty("ValidateInputs", false);

  alg->executeAsChildAlg();

  API::MatrixWorkspace_sptr outws = alg->getProperty("OutputWorkspace");

  // FIXED : Restore the original spectrum Nos to spectra from ws1
  for (size_t i = 0; i < nspec1; ++i) {
    specnum_t tmpspecNo = outws->getSpectrum(i).getSpectrumNo();
    outws->getSpectrum(i).setSpectrumNo(origspecNos[i]);

    g_log.information() << "[DBx540] Conjoined spectrum " << i
                        << ": restore spectrum number to "
                        << outws->getSpectrum(i).getSpectrumNo()
                        << " from spectrum number = " << tmpspecNo << ".\n";
  }

  // Rename spectrum number
  if (offset >= 1) {
    for (size_t i = 0; i < nspec2; ++i) {
      specnum_t newspecid = maxspecNo1 + static_cast<specnum_t>((i) + offset);
      outws->getSpectrum(nspec1 + i).setSpectrumNo(newspecid);
      // ISpectrum* spec = outws->getSpectrum(nspec1+i);
      // if (spec)
      // spec->setSpectrumNo(3);
    }
  }

  return outws;
}

/**
 * Read spectra from the DAE
 * @param period :: Current period index
 * @param index :: First spectrum index
 * @param count :: Number of spectra to read
 * @param workspace :: Workspace to store the data
 * @param workspaceIndex :: index in workspace to store data
 */
void ISISHistoDataListener::getData(int period, int index, int count,
                                    API::MatrixWorkspace_sptr workspace,
                                    size_t workspaceIndex) {
  const int numberOfBins = m_numberOfBins[m_timeRegime];
  const size_t bufferSize = count * (numberOfBins + 1) * sizeof(int);
  std::vector<int> dataBuffer(bufferSize);
  // Read in spectra from DAE
  int ndims = 2, dims[2];
  dims[0] = count;
  dims[1] = numberOfBins + 1;

  int spectrumIndex = index + period * (m_totalNumberOfSpectra + 1);
  if (IDCgetdat(m_daeHandle, spectrumIndex, count, dataBuffer.data(), dims,
                &ndims) != 0) {
    g_log.error("Unable to read DATA from DAE " + m_daeName);
    throw Kernel::Exception::FileError("Unable to read DATA from DAE ",
                                       m_daeName);
  }

  for (size_t i = 0; i < static_cast<size_t>(count); ++i) {
    size_t wi = workspaceIndex + i;
    workspace->setX(wi, m_bins[m_timeRegime]);
    MantidVec &y = workspace->dataY(wi);
    MantidVec &e = workspace->dataE(wi);
    workspace->getSpectrum(wi)->setSpectrumNo(index + static_cast<specid_t>(i));
    size_t shift = i * (numberOfBins + 1) + 1;
    y.assign(dataBuffer.begin() + shift, dataBuffer.begin() + shift + y.size());
    std::transform(y.begin(), y.end(), e.begin(), dblSqrt);
  }
}

    /** This function will check how to group spectra when calculating median
     *
     *
     */
    std::vector<std::vector<size_t> > DetectorDiagnostic::makeMap(API::MatrixWorkspace_sptr countsWS)
    {
      std::multimap<Mantid::Geometry::ComponentID,size_t> mymap;

      Geometry::Instrument_const_sptr instrument = countsWS->getInstrument();
      if (m_parents==0)
      {
        return makeInstrumentMap(countsWS);
      }
      if (!instrument)
      {
        g_log.warning("Workspace has no instrument. LevelsUP is ignored");
        return makeInstrumentMap(countsWS);
      }

      //check if not grouped. If grouped, it will throw
      if ( countsWS->hasGroupedDetectors() )
      {
        throw std::runtime_error("Median detector test: not able to create detector to spectra map. Try with LevelUp=0.");
      }

      for(size_t i=0;i < countsWS->getNumberHistograms();i++)
      {
        detid_t d=(*((countsWS->getSpectrum(i))->getDetectorIDs().begin()));
        std::vector<boost::shared_ptr<const Mantid::Geometry::IComponent> > anc=instrument->getDetector(d)->getAncestors();
        //std::vector<boost::shared_ptr<const IComponent> > anc=(*(countsWS->getSpectrum(i)->getDetectorIDs().begin()))->getAncestors();
        if (anc.size()<static_cast<size_t>(m_parents))
        {
          g_log.warning("Too many levels up. Will ignore LevelsUp");
          m_parents=0;
          return makeInstrumentMap(countsWS);
        }
        mymap.insert(std::pair<Mantid::Geometry::ComponentID,size_t>(anc[m_parents-1]->getComponentID(),i));
      }

      std::vector<std::vector<size_t> > speclist;
      std::vector<size_t>  speclistsingle;

      std::multimap<Mantid::Geometry::ComponentID,size_t>::iterator m_it, s_it;

      for (m_it = mymap.begin();  m_it != mymap.end();  m_it = s_it)
      {
        Mantid::Geometry::ComponentID theKey = (*m_it).first;

        std::pair<std::multimap<Mantid::Geometry::ComponentID,size_t>::iterator,std::multimap<Mantid::Geometry::ComponentID,size_t>::iterator> keyRange = mymap.equal_range(theKey);

        // Iterate over all map elements with key == theKey
        speclistsingle.clear();
        for (s_it = keyRange.first;  s_it != keyRange.second;  ++s_it)
        {
          speclistsingle.push_back( (*s_it).second );
        }
        speclist.push_back(speclistsingle);
      }

      return speclist;
    }

/** Convert X units from nano-secs to micro-secs and shift to start at m_tMin
* @param inputWs :: [input/output] workspace to convert
*/
void PhaseQuadMuon::convertToMicroSecs (API::MatrixWorkspace_sptr inputWs)
{
  for (size_t h=0; h<inputWs->getNumberHistograms(); h++)
  {
    auto spec = inputWs->getSpectrum(h);
    for (int t=0; t<m_nData+1; t++)
    {
     spec->dataX()[t] = spec->dataX()[t]/1000+m_tMin;
    }
  }
}

/***
 * This will ensure the spectrum numbers do not overlap by starting the second
 *on at the first + 1
 *
 * @param ws1 The first workspace supplied to the algorithm.
 * @param ws2 The second workspace supplied to the algorithm.
 * @param output The workspace that is going to be returned by the algorithm.
 */
void ConjoinWorkspaces::fixSpectrumNumbers(API::MatrixWorkspace_const_sptr ws1,
                                           API::MatrixWorkspace_const_sptr ws2,
                                           API::MatrixWorkspace_sptr output) {
  bool needsFix(false);

  if (this->getProperty("CheckOverlapping")) {
    // If CheckOverlapping is required, then either skip fixing spectrum number
    // or get stopped by an exception
    if (!m_overlapChecked)
      checkForOverlap(ws1, ws2, true);
    needsFix = false;
  } else {
    // It will be determined later whether spectrum number needs to be fixed.
    needsFix = true;
  }
  if (!needsFix)
    return;

  // is everything possibly ok?
  specid_t min;
  specid_t max;
  getMinMax(output, min, max);
  if (max - min >= static_cast<specid_t>(
                       output->getNumberHistograms())) // nothing to do then
    return;

  // information for remapping the spectra numbers
  specid_t ws1min;
  specid_t ws1max;
  getMinMax(ws1, ws1min, ws1max);

  // change the axis by adding the maximum existing spectrum number to the
  // current value
  for (size_t i = ws1->getNumberHistograms(); i < output->getNumberHistograms();
       i++) {
    specid_t origid;
    origid = output->getSpectrum(i)->getSpectrumNo();
    output->getSpectrum(i)->setSpectrumNo(origid + ws1max);
  }
}

    /**
     * Mask the outlier values to get a better median value.
     * @param medianvec The median value calculated from the current counts.
     * @param countsWS The counts workspace. Any outliers will be masked here.
     * @param indexmap Index map.
     * @returns The number failed.
     */
    int MedianDetectorTest::maskOutliers(const std::vector<double> medianvec, API::MatrixWorkspace_sptr countsWS,std::vector<std::vector<size_t> > indexmap)
    {

      // Fractions of the median
      const double out_lo = getProperty("LowOutlier");
      const double out_hi = getProperty("HighOutlier");

      int numFailed(0);

      bool checkForMask = false;
      Geometry::Instrument_const_sptr instrument = countsWS->getInstrument();
      if (instrument != NULL)
      {
        checkForMask = ((instrument->getSource() != NULL) && (instrument->getSample() != NULL));
      }

      for (size_t i=0; i<indexmap.size();  ++i)
      {
        std::vector<size_t> hists=indexmap.at(i);
        double median=medianvec.at(i);

        PARALLEL_FOR1(countsWS)
        for(int j = 0; j < static_cast<int>(hists.size()); ++j)
        {
          const double value = countsWS->readY(hists.at(j))[0];
          if ((value == 0.) && checkForMask)
          {
            const std::set<detid_t>& detids = countsWS->getSpectrum(hists.at(j))->getDetectorIDs();
            if (instrument->isDetectorMasked(detids))
            {
              numFailed -= 1; // it was already masked
            }
          }
          if( (value < out_lo*median) && (value > 0.0) )
          {
            countsWS->maskWorkspaceIndex(hists.at(j));
            PARALLEL_ATOMIC
            ++numFailed;
          }
          else if( value > out_hi*median )
          {
            countsWS->maskWorkspaceIndex(hists.at(j));
            PARALLEL_ATOMIC
            ++numFailed;
          }
        }
        PARALLEL_CHECK_INTERUPT_REGION
      }

      return numFailed;
    }

/** Apply dead time correction to spectra in inputWs and create temporary workspace with corrected spectra
* @param inputWs :: [input] input workspace containing spectra to correct
* @param deadTimeTable :: [input] table containing dead times
* @param tempWs :: [output] workspace containing corrected spectra
*/
void PhaseQuadMuon::deadTimeCorrection(API::MatrixWorkspace_sptr inputWs, API::ITableWorkspace_sptr deadTimeTable, API::MatrixWorkspace_sptr& tempWs)
{

  // Apply correction only from t = m_tPulseOver
  // To do so, we first apply corrections to the whole spectrum (ApplyDeadTimeCorr
  // does not allow to select a range in the spectrum)
  // Then we recover counts from 0 to m_tPulseOver

  auto alg = this->createChildAlgorithm("ApplyDeadTimeCorr",-1,-1);
  alg->initialize();
  alg->setProperty("DeadTimeTable", deadTimeTable);
  alg->setPropertyValue("InputWorkspace", inputWs->getName());
  alg->setPropertyValue("OutputWorkspace", inputWs->getName()+"_deadtime");
  bool sucDeadTime = alg->execute();
  if (!sucDeadTime)
  {
    g_log.error() << "PhaseQuad: Unable to apply dead time corrections" << std::endl;
    throw std::runtime_error("PhaseQuad: Unable to apply dead time corrections");
  }
  tempWs = alg->getProperty("OutputWorkspace");

  // Now recover counts from t=0 to m_tPulseOver
  // Errors are set to m_bigNumber
  for (int h=0; h<m_nHist; h++)
  {
    auto specOld = inputWs->getSpectrum(h);
    auto specNew = tempWs->getSpectrum(h);

    for (int t=0; t<m_tPulseOver; t++)
    {
      specNew->dataY()[t] = specOld->dataY()[t];
      specNew->dataE()[t] = m_bigNumber;
    }
  }

}

/**
*  Only to be used if the KeepUnGrouped property is true, moves the spectra that were not selected
*  to be in a group to the end of the output spectrum
*  @param unGroupedSet :: list of WORKSPACE indexes that were included in a group
*  @param inputWS :: user selected input workspace for the algorithm
*  @param outputWS :: user selected output workspace for the algorithm
*  @param outIndex :: the next spectra index available after the grouped spectra
*/
void GroupDetectors2::moveOthers(const std::set<int64_t> &unGroupedSet, API::MatrixWorkspace_const_sptr inputWS, API::MatrixWorkspace_sptr outputWS,
         size_t outIndex)
{
  g_log.debug() << "Starting to copy the ungrouped spectra" << std::endl;
  double prog4Copy = (1. - 1.*static_cast<double>(m_FracCompl))/static_cast<double>(unGroupedSet.size());

  std::set<int64_t>::const_iterator copyFrIt = unGroupedSet.begin();
  // go thorugh all the spectra in the input workspace
  for ( ; copyFrIt != unGroupedSet.end(); ++copyFrIt )
  {
    if( *copyFrIt == USED ) continue; //Marked as not to be used
    size_t sourceIndex = static_cast<size_t>(*copyFrIt);

    // The input spectrum we'll copy
    const ISpectrum * inputSpec = inputWS->getSpectrum(sourceIndex);

    // Destination of the copying
    ISpectrum * outputSpec = outputWS->getSpectrum(outIndex);

    // Copy the data
    outputSpec->dataX() = inputSpec->dataX();
    outputSpec->dataY() = inputSpec->dataY();
    outputSpec->dataE() = inputSpec->dataE();

    // Spectrum numbers etc.
    outputSpec->setSpectrumNo(inputSpec->getSpectrumNo());
    outputSpec->clearDetectorIDs();
    outputSpec->addDetectorIDs( inputSpec->getDetectorIDs() );

    // go to the next free index in the output workspace
    outIndex ++;
    // make regular progress reports and check for cancelling the algorithm
    if ( outIndex % INTERVAL == 0 )
    {
      m_FracCompl += INTERVAL*prog4Copy;
      if ( m_FracCompl > 1.0 )
      {
        m_FracCompl = 1.0;
      }
      progress(m_FracCompl);
      interruption_point();
    }
  }
  // Refresh the spectraDetectorMap
  outputWS->generateSpectraMap();

  g_log.debug() << name() << " copied " << unGroupedSet.size()-1 << " ungrouped spectra\n";
}

/** Fits each spectrum in the workspace to f(x) = A * sin( w * x + p)
* @param ws :: [input] The workspace to fit
* @param freq :: [input] Hint for the frequency (w)
* @param groupName :: [input] The name of the output workspace group
* @param resTab :: [output] Table workspace storing the asymmetries and phases
* @param resGroup :: [output] Workspace group storing the fitting results
*/
void CalMuonDetectorPhases::fitWorkspace(const API::MatrixWorkspace_sptr &ws,
                                         double freq, std::string groupName,
                                         API::ITableWorkspace_sptr &resTab,
                                         API::WorkspaceGroup_sptr &resGroup) {

  int nhist = static_cast<int>(ws->getNumberHistograms());

  // Create the fitting function f(x) = A * sin ( w * x + p )
  // The same function and initial parameters are used for each fit
  std::string funcStr = createFittingFunction(freq, true);

  // Set up results table
  resTab->addColumn("int", "Spectrum number");
  resTab->addColumn("double", "Asymmetry");
  resTab->addColumn("double", "Phase");

  // Loop through fitting all spectra individually
  const static std::string success = "success";
  for (int wsIndex = 0; wsIndex < nhist; wsIndex++) {
    reportProgress(wsIndex, nhist);
    auto fit = createChildAlgorithm("Fit");
    fit->initialize();
    fit->setPropertyValue("Function", funcStr);
    fit->setProperty("InputWorkspace", ws);
    fit->setProperty("WorkspaceIndex", wsIndex);
    fit->setProperty("CreateOutput", true);
    fit->setPropertyValue("Output", groupName);
    fit->execute();

    std::string status = fit->getProperty("OutputStatus");
    if (!fit->isExecuted() || status != success) {
      std::ostringstream error;
      error << "Fit failed for spectrum at workspace index " << wsIndex;
      error << ": " << status;
      throw std::runtime_error(error.str());
    }

    API::MatrixWorkspace_sptr fitOut = fit->getProperty("OutputWorkspace");
    resGroup->addWorkspace(fitOut);
    API::ITableWorkspace_sptr tab = fit->getProperty("OutputParameters");
    // Now we have our fitting results stored in tab
    // but we need to extract the relevant information, i.e.
    // the detector phases (parameter 'p') and asymmetries ('A')
    const auto &spectrum = ws->getSpectrum(static_cast<size_t>(wsIndex));
    extractDetectorInfo(tab, resTab, spectrum.getSpectrumNo());
  }
}

  /** If there is an overlap in spectrum numbers between ws1 and ws2,
   * then the spectrum numbers are reset as a simple 1-1 correspondence
   * with the workspace index.
   *
   * @param ws1 The first workspace supplied to the algorithm.
   * @param ws2 The second workspace supplied to the algorithm.
   * @param output The workspace that is going to be returned by the algorithm.
   */
  void AppendSpectra::fixSpectrumNumbers(API::MatrixWorkspace_const_sptr ws1, API::MatrixWorkspace_const_sptr ws2,
                                             API::MatrixWorkspace_sptr output)
  {
    specid_t ws1min;
    specid_t ws1max;
    getMinMax(ws1, ws1min, ws1max);

    specid_t ws2min;
    specid_t ws2max;
    getMinMax(ws2, ws2min, ws2max);

    // is everything possibly ok?
    if (ws2min > ws1max)
      return;

    // change the axis by adding the maximum existing spectrum number to the current value
    for (size_t i = 0; i < output->getNumberHistograms(); i++)
      output->getSpectrum(i)->setSpectrumNo( specid_t(i) );
  }

/**
 * The main method to calculate the ring profile for workspaces based on
 *instruments.
 *
 * It will iterate over all the spectrum inside the workspace.
 * For each spectrum, it will use the RingProfile::getBinForPixel method to
 *identify
 * where, in the output_bins, the sum of all the spectrum values should be
 *placed in.
 *
 * @param inputWS: pointer to the input workspace
 * @param output_bins: the reference to the vector to be filled with the
 *integration values
 */
void RingProfile::processInstrumentRingProfile(
    const API::MatrixWorkspace_sptr inputWS, std::vector<double> &output_bins) {

  for (int i = 0; i < static_cast<int>(inputWS->getNumberHistograms()); i++) {
    m_progress->report("Computing ring bins positions for detectors");
    // for the detector based, the positions will be taken from the detector
    // itself.
    try {
      Mantid::Geometry::IDetector_const_sptr det = inputWS->getDetector(i);

      // skip monitors
      if (det->isMonitor()) {
        continue;
      }

      // this part will be executed if the instrument is attached to the
      // workspace

      // get the bin position
      int bin_n = getBinForPixel(det);

      if (bin_n < 0) // -1 is the agreement for an invalid bin, or outside the
                     // ring being integrated
        continue;

      g_log.debug() << "Bin for the index " << i << " = " << bin_n
                    << " Pos = " << det->getPos() << std::endl;

      // get the reference to the spectrum
      auto spectrum_pt = inputWS->getSpectrum(i);
      const MantidVec &refY = spectrum_pt->dataY();
      // accumulate the values of this spectrum inside this bin
      for (size_t sp_ind = 0; sp_ind < inputWS->blocksize(); sp_ind++)
        output_bins[bin_n] += refY[sp_ind];

    } catch (Kernel::Exception::NotFoundError &ex) {
      g_log.information() << "It found that detector for " << i
                          << " is not valid. " << ex.what() << std::endl;
      continue;
    }
  }
}

/** Create the final output workspace after converting the X axis
* @returns the final output workspace
*
* @param progress :: Progress indicator
* @param targetUnit :: Target conversion unit
* @param inputWS :: Input workspace
* @param nHist :: Stores the number of histograms
*/
MatrixWorkspace_sptr ConvertSpectrumAxis2::createOutputWorkspace(
    API::Progress &progress, const std::string &targetUnit,
    API::MatrixWorkspace_sptr &inputWS, size_t nHist) {
  // Create the output workspace. Can not re-use the input one because the
  // spectra are re-ordered.
  MatrixWorkspace_sptr outputWorkspace = WorkspaceFactory::Instance().create(
      inputWS, m_indexMap.size(), inputWS->x(0).size(), inputWS->y(0).size());

  // Now set up a new numeric axis holding the theta values corresponding to
  // each spectrum.
  auto const newAxis = new NumericAxis(m_indexMap.size());
  outputWorkspace->replaceAxis(1, newAxis);

  progress.setNumSteps(nHist + m_indexMap.size());

  // Set the units of the axis.
  if (targetUnit == "theta" || targetUnit == "Theta" ||
      targetUnit == "signed_theta" || targetUnit == "SignedTheta") {
    newAxis->unit() = boost::make_shared<Units::Degrees>();
  } else if (targetUnit == "ElasticQ") {
    newAxis->unit() = UnitFactory::Instance().create("MomentumTransfer");
  } else if (targetUnit == "ElasticQSquared") {
    newAxis->unit() = UnitFactory::Instance().create("QSquared");
  }

  std::multimap<double, size_t>::const_iterator it;
  size_t currentIndex = 0;
  for (it = m_indexMap.begin(); it != m_indexMap.end(); ++it) {
    // Set the axis value.
    newAxis->setValue(currentIndex, it->first);
    // Copy over the data.
    outputWorkspace->setHistogram(currentIndex, inputWS->histogram(it->second));
    // We can keep the spectrum numbers etc.
    outputWorkspace->getSpectrum(currentIndex)
        .copyInfoFrom(inputWS->getSpectrum(it->second));
    ++currentIndex;

    progress.report("Creating output workspace...");
  }
  return outputWorkspace;
}

/**
 * Construct an "empty" output workspace in virtual-lambda for summation in Q.
 *
 * @param detectorWS [in] :: the input workspace
 * @param indices [in] :: the workspace indices of the foreground histograms
 * @param refAngles [in] :: the reference angles
 * @return :: a 1D workspace where y values are all zero
 */
API::MatrixWorkspace_sptr ReflectometrySumInQ::constructIvsLamWS(
    const API::MatrixWorkspace &detectorWS,
    const Indexing::SpectrumIndexSet &indices, const Angles &refAngles) {

  // Calculate the number of bins based on the min/max wavelength, using
  // the same bin width as the input workspace
  const auto &edges = detectorWS.binEdges(refAngles.referenceWSIndex);
  const double binWidth =
      (edges.back() - edges.front()) / static_cast<double>(edges.size());
  const auto wavelengthRange =
      findWavelengthMinMax(detectorWS, indices, refAngles);
  if (std::abs(wavelengthRange.max - wavelengthRange.min) < binWidth) {
    throw std::runtime_error("Given wavelength range too small.");
  }
  const int numBins = static_cast<int>(
      std::ceil((wavelengthRange.max - wavelengthRange.min) / binWidth));
  // Construct the histogram with these X values. Y and E values are zero.
  const HistogramData::BinEdges bins(
      numBins + 1,
      HistogramData::LinearGenerator(wavelengthRange.min, binWidth));
  const HistogramData::Counts counts(numBins, 0.);
  const HistogramData::Histogram modelHistogram(std::move(bins),
                                                std::move(counts));
  // Create the output workspace
  API::MatrixWorkspace_sptr outputWS =
      DataObjects::create<DataObjects::Workspace2D>(detectorWS, 1,
                                                    std::move(modelHistogram));

  // Set the detector IDs and specturm number from the twoThetaR detector.
  const auto &thetaSpec = detectorWS.getSpectrum(refAngles.referenceWSIndex);
  auto &outSpec = outputWS->getSpectrum(0);
  outSpec.clearDetectorIDs();
  outSpec.addDetectorIDs(thetaSpec.getDetectorIDs());
  outSpec.setSpectrumNo(thetaSpec.getSpectrumNo());

  return outputWS;
}