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

/** 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;
    }
  }
}

/** Loads, checks and passes back the values passed to the algorithm
 * @param whiteBeam1 :: A white beam vanadium spectrum that will be used to
 * check detector efficiency variations
 * @param whiteBeam2 :: The other white beam vanadium spectrum from the same
 * instrument to use for comparison
 * @param variation :: The maximum fractional variation above the median that is
 * allowed for god detectors
 * @param startWsIndex :: Index number of the first spectrum to use
 * @param endWsIndex :: Index number of the last spectrum to use
 * @throw invalid_argument if there is an incapatible property value and so the
 * algorithm can't continue
 */
void DetectorEfficiencyVariation::retrieveProperties(
    API::MatrixWorkspace_sptr &whiteBeam1,
    API::MatrixWorkspace_sptr &whiteBeam2, double &variation, int &startWsIndex,
    int &endWsIndex) {
  whiteBeam1 = getProperty("WhiteBeamBase");
  whiteBeam2 = getProperty("WhiteBeamCompare");
  if (whiteBeam1->getInstrument()->getName() !=
      whiteBeam2->getInstrument()->getName()) {
    throw std::invalid_argument("The two input white beam vanadium workspaces "
                                "must be from the same instrument");
  }
  int maxWsIndex = static_cast<int>(whiteBeam1->getNumberHistograms()) - 1;
  if (maxWsIndex !=
      static_cast<int>(whiteBeam2->getNumberHistograms()) -
          1) { // we would get a crash later on if this were not true
    throw std::invalid_argument("The input white beam vanadium workspaces must "
                                "be have the same number of histograms");
  }

  variation = getProperty("Variation");

  startWsIndex = getProperty("StartWorkspaceIndex");
  if ((startWsIndex < 0) || (startWsIndex > maxWsIndex)) {
    g_log.warning("StartWorkspaceIndex out of range, changed to 0");
    startWsIndex = 0;
  }
  endWsIndex = getProperty("EndWorkspaceIndex");
  if (endWsIndex == Mantid::EMPTY_INT())
    endWsIndex = maxWsIndex;
  if ((endWsIndex < 0) || (endWsIndex > maxWsIndex)) {
    g_log.warning(
        "EndWorkspaceIndex out of range, changed to max Workspace number");
    endWsIndex = maxWsIndex;
  }
  if ((endWsIndex < startWsIndex)) {
    g_log.warning(
        "EndWorkspaceIndex can not be less than the StartWorkspaceIndex, "
        "changed to max Workspace number");
    endWsIndex = maxWsIndex;
  }
}

    std::vector<std::vector<size_t> > DetectorDiagnostic::makeInstrumentMap(API::MatrixWorkspace_sptr countsWS)
    {
      std::vector<std::vector<size_t> > mymap;
      std::vector<size_t> single;

      for(size_t i=0;i < countsWS->getNumberHistograms();i++)
      {
        single.push_back(i);
      }
      mymap.push_back(single);
      return mymap;
    }

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

/** Returns a given spectrum as a complex number
* @param inWS :: [input] The input workspace containing all the spectra
* @param spec :: [input] The spectrum of interest
* @param errors :: [input] If true, returns the errors, otherwise returns the
* counts
* @return : Spectrum 'spec' as a complex vector
*/
std::vector<double> MaxEnt::toComplex(const API::MatrixWorkspace_sptr &inWS,
                                      size_t spec, bool errors) {

  std::vector<double> result(inWS->blocksize() * 2);

  if (inWS->getNumberHistograms() % 2)
    throw std::invalid_argument(
        "Cannot convert input workspace to complex data");

  size_t nspec = inWS->getNumberHistograms() / 2;

  if (!errors) {
    for (size_t i = 0; i < inWS->blocksize(); i++) {
      result[2 * i] = inWS->readY(spec)[i];
      result[2 * i + 1] = inWS->readY(spec + nspec)[i];
    }
  } else {
    for (size_t i = 0; i < inWS->blocksize(); i++) {
      result[2 * i] = inWS->readE(spec)[i];
      result[2 * i + 1] = inWS->readE(spec + nspec)[i];
    }
  }
  return result;
}

/** Divide each bin by the width of its q bin.
 *  @param outputWS :: The output workspace
 *  @param qAxis ::    A vector of the q bin boundaries
 */
void SofQWCentre::makeDistribution(API::MatrixWorkspace_sptr outputWS,
                                   const std::vector<double> qAxis) {
  std::vector<double> widths(qAxis.size());
  std::adjacent_difference(qAxis.begin(), qAxis.end(), widths.begin());

  const size_t numQBins = outputWS->getNumberHistograms();
  for (size_t i = 0; i < numQBins; ++i) {
    auto &Y = outputWS->mutableY(i);
    auto &E = outputWS->mutableE(i);
    std::transform(Y.begin(), Y.end(), Y.begin(),
                   std::bind2nd(std::divides<double>(), widths[i + 1]));
    std::transform(E.begin(), E.end(), E.begin(),
                   std::bind2nd(std::divides<double>(), widths[i + 1]));
  }
}

void SaveFITS::writeFITSHeaderAxesSizes(const API::MatrixWorkspace_sptr img,
                                        std::ofstream &file) {
  const std::string sizeX = std::to_string(img->blocksize());
  const std::string sizeY = std::to_string(img->getNumberHistograms());

  const size_t fieldWidth = 20;
  std::stringstream axis1;
  axis1 << "NAXIS1  = " << std::setw(fieldWidth) << sizeX
        << " / length of data axis 1";
  writeFITSHeaderEntry(axis1.str(), file);

  std::stringstream axis2;
  axis2 << "NAXIS2  = " << std::setw(fieldWidth) << sizeY
        << " / length of data axis 2";
  writeFITSHeaderEntry(axis2.str(), file);
}

/** 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());
  }
}

/**
 * Validation of the inputs of the RingProfile algorithm. 
 *
 * Inside this method, the Workspace is considered a 2D Matrix, where each spectrum is
 * the rows of the matrix and have the variation in axis0. The columns of the matrix 
 * is the position of dataX(0)
 * 
 * The main validation are: 
 *  - the centre of the ring is inside the image it self.
 *  - The minimum ring is smaller than the limits of the image to allow 
 * @param inputWS: pointer to the input workspace
*/
void RingProfile::checkInputsForNumericWorkspace(const API::MatrixWorkspace_sptr inputWS){
  g_log.notice() << "CheckingInputs For Numeric Workspace" << std::endl; 

  // The Axis0 is defined by the values of readX inside the spectra of the workspace. 
  // The limits of this axis will be get by inspection of the readX vector taking the first 
  // and the last value. 
  
  // check that centre is inside the range available for the instrument
  const MantidVec & refX  = inputWS->readX(inputWS->getNumberHistograms()/2);  
  // get the limits of the axis 0 (X)
  double min_v_x, max_v_x;
  min_v_x = std::min(refX[0], refX[refX.size() -1]); 
  max_v_x = std::max(refX[0], refX[refX.size() -1]); 
  g_log.notice() << "Limits X = " << min_v_x << " " << max_v_x << std::endl;   
  // check centre is inside the X domain
  if ( centre_x < min_v_x || centre_x > max_v_x){
    std::stringstream s;
    s << "The input value for centre (X="<< centre_x << ") is outside the limits of the instrument ["
      << min_v_x << ", "<< max_v_x << "]" ; 
    throw std::invalid_argument(s.str());
  }        
  

  // The Axis1 is defined by the spectra inside the workspace. Its limits and values are given by the 
  // ws->getAxis(1)

  // get the limits of the axis1 (Y)
  API::NumericAxis *oldAxis2 = dynamic_cast<API::NumericAxis*>(inputWS->getAxis(1) );
  // we cannot have the positions in Y direction without a NumericAxis  
  if( !oldAxis2 ) 
    throw std::invalid_argument("Vertical axis is not a numeric axis. If it is a spectra axis try running ConvertSpectrumAxis first.");
  double min_v_y = std::min(oldAxis2->getMin(), oldAxis2->getMax()); 
  double max_v_y = std::max(oldAxis2->getMin(), oldAxis2->getMax()); 
  g_log.notice() << "Limits Y = " << min_v_y << " " << max_v_y << std::endl;   
  // check centre is inside the Y domain
  if (centre_y < min_v_y || centre_y > max_v_y){
    std::stringstream s; 
    s << "The input value for centre (Y=" << centre_y << ") is outside the limits of the instrument ["
      << min_v_y << ", " << max_v_y << "]" ;
    throw std::invalid_argument(s.str()); 
  }
  g_log.notice() << "Centre: " << centre_x << "  " << centre_y << std::endl; 
  // check minradius is inside the limits of the region of the instrument
  if (centre_x - min_radius > max_v_x || centre_x + min_radius < min_v_x || centre_y - min_radius > max_v_y || centre_y + min_radius < min_v_y)
    throw std::invalid_argument("The minimun radius is outside the region of the instrument");

}

/** 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()));
    auto anc = instrument->getDetector(d)->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.emplace(anc[m_parents - 1]->getComponentID(), i);
  }

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

  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;
    auto keyRange = mymap.equal_range(theKey);
    // Iterate over all map elements with key == theKey
    std::vector<size_t> speclistsingle;
    for (s_it = keyRange.first; s_it != keyRange.second; ++s_it) {
      speclistsingle.push_back(s_it->second);
    }
    speclist.push_back(std::move(speclistsingle));
  }

  return speclist;
}

/** Assuming that the input workspace is a Numeric Image where the pixel
 *positions depend on their
 *  relative position inside the workspace, this function extracts the position
 *of the first and last
 *  pixel of the image.
 *
 *  It is important that the input workspace must be a numeric image, and not an
 *instrument related workspace.
 *  The function will raise exception (std::invalid_argument) if an invalid
 *input is give.
 *
 *  @see RadiusSum::inputWorkspaceHasInstrumentAssociated for reference.
 *
 *  @param inWS reference to the workspace
 *  @return a list of values that defines the limits of the image in this order:
 *Xmin, Xmax, Ymin, Ymax
 */
std::vector<double>
RadiusSum::getBoundariesOfNumericImage(API::MatrixWorkspace_sptr inWS) {

  // horizontal axis

  // get the pixel position in the horizontal axis from the middle of the image.
  const MantidVec &refX = inWS->readX(inWS->getNumberHistograms() / 2);

  double min_x, max_x;

  const double &first_x(refX[0]);
  const double &last_x(refX[refX.size() - 1]);
  if (first_x < last_x) {
    min_x = first_x;
    max_x = last_x;
  } else {
    min_x = last_x;
    max_x = first_x;
  }

  // vertical axis
  API::NumericAxis *verticalAxis =
      dynamic_cast<API::NumericAxis *>(inWS->getAxis(1));
  if (!verticalAxis)
    throw std::invalid_argument("Vertical axis is not a numeric axis. Can not "
                                "find the limits of the image.");

  double min_y, max_y;
  min_y = verticalAxis->getMin();
  max_y = verticalAxis->getMax();

  // check the assumption made that verticalAxis will provide the correct
  // answer.
  if (min_y > max_y) {
    throw std::logic_error("Failure to get the boundaries of this image. "
                           "Internal logic error. Please, inform MantidHelp");
  }

  std::vector<double> output(4); // output = {min_x, max_x, min_y, max_y}; not
                                 // supported in all compilers
  output[0] = min_x;
  output[1] = max_x;
  output[2] = min_y;
  output[3] = max_y;
  return output;
}

  /** 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;
    }
  }
}

/** Executes the algorithm
 *
 */
void SplineBackground::exec()
{

  API::MatrixWorkspace_sptr inWS = getProperty("InputWorkspace");
  int spec = getProperty("WorkspaceIndex");

  if (spec > static_cast<int>(inWS->getNumberHistograms()))
    throw std::out_of_range("WorkspaceIndex is out of range.");

  const MantidVec& X = inWS->readX(spec);
  const MantidVec& Y = inWS->readY(spec);
  const MantidVec& E = inWS->readE(spec);
  const bool isHistogram = inWS->isHistogramData();

  const int ncoeffs = getProperty("NCoeff");
  const int k = 4; // order of the spline + 1 (cubic)
  const int nbreak = ncoeffs - (k - 2);

  if (nbreak <= 0)
    throw std::out_of_range("Too low NCoeff");

  gsl_bspline_workspace *bw;
  gsl_vector *B;

  gsl_vector *c, *w, *x, *y;
  gsl_matrix *Z, *cov;
  gsl_multifit_linear_workspace *mw;
  double chisq;

  int n = static_cast<int>(Y.size());
  bool isMasked = inWS->hasMaskedBins(spec);
  std::vector<int> masked(Y.size());
  if (isMasked)
  {
    for(API::MatrixWorkspace::MaskList::const_iterator it=inWS->maskedBins(spec).begin();it!=inWS->maskedBins(spec).end();++it)
      masked[it->first] = 1;
    n -= static_cast<int>(inWS->maskedBins(spec).size());
  }

  if (n < ncoeffs)
  {
    g_log.error("Too many basis functions (NCoeff)");
    throw std::out_of_range("Too many basis functions (NCoeff)");
  }

  /* allocate a cubic bspline workspace (k = 4) */
  bw = gsl_bspline_alloc(k, nbreak);
  B = gsl_vector_alloc(ncoeffs);

  x = gsl_vector_alloc(n);
  y = gsl_vector_alloc(n);
  Z = gsl_matrix_alloc(n, ncoeffs);
  c = gsl_vector_alloc(ncoeffs);
  w = gsl_vector_alloc(n);
  cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
  mw = gsl_multifit_linear_alloc(n, ncoeffs);

  /* this is the data to be fitted */
  int j = 0;
  for (MantidVec::size_type i = 0; i < Y.size(); ++i)
  {
    if (isMasked && masked[i]) continue;
    gsl_vector_set(x, j, (isHistogram ? (0.5*(X[i]+X[i+1])) : X[i])); // Middle of the bins, if a histogram
    gsl_vector_set(y, j, Y[i]);
    gsl_vector_set(w, j, E[i]>0.?1./(E[i]*E[i]):0.);

    ++j;
  }

  if (n != j)
  {
    gsl_bspline_free(bw);
    gsl_vector_free(B);
    gsl_vector_free(x);
    gsl_vector_free(y);
    gsl_matrix_free(Z);
    gsl_vector_free(c);
    gsl_vector_free(w);
    gsl_matrix_free(cov);
    gsl_multifit_linear_free(mw);

    throw std::runtime_error("Assertion failed: n != j");
  }

  double xStart = X.front();
  double xEnd =   X.back();

  /* use uniform breakpoints */
  gsl_bspline_knots_uniform(xStart, xEnd, bw);

  /* construct the fit matrix X */
  for (int i = 0; i < n; ++i)
  {
    double xi=gsl_vector_get(x, i);

    /* compute B_j(xi) for all j */
    gsl_bspline_eval(xi, B, bw);
//.........这里部分代码省略.........

/** Carries out the bin-by-bin normalization
 *  @param inputWorkspace The input workspace
 *  @param outputWorkspace The result workspace
 */
void NormaliseToMonitor::normaliseBinByBin(
    const API::MatrixWorkspace_sptr &inputWorkspace,
    API::MatrixWorkspace_sptr &outputWorkspace) {
  EventWorkspace_sptr inputEvent =
      boost::dynamic_pointer_cast<EventWorkspace>(inputWorkspace);

  // Only create output workspace if different to input one
  if (outputWorkspace != inputWorkspace) {
    if (inputEvent) {
      outputWorkspace = inputWorkspace->clone();
    } else
      outputWorkspace = WorkspaceFactory::Instance().create(inputWorkspace);
  }
  auto outputEvent =
      boost::dynamic_pointer_cast<EventWorkspace>(outputWorkspace);

  // Get hold of the monitor spectrum
  const auto &monX = m_monitor->binEdges(0);
  auto monY = m_monitor->counts(0);
  auto monE = m_monitor->countStandardDeviations(0);
  // Calculate the overall normalization just the once if bins are all matching
  if (m_commonBins)
    this->normalisationFactor(monX, monY, monE);

  const size_t numHists = inputWorkspace->getNumberHistograms();
  auto specLength = inputWorkspace->blocksize();
  // Flag set when a division by 0 is found
  bool hasZeroDivision = false;
  Progress prog(this, 0.0, 1.0, numHists);
  // Loop over spectra
  PARALLEL_FOR_IF(
      Kernel::threadSafe(*inputWorkspace, *outputWorkspace, *m_monitor))
  for (int64_t i = 0; i < int64_t(numHists); ++i) {
    PARALLEL_START_INTERUPT_REGION
    prog.report();

    const auto &X = inputWorkspace->binEdges(i);
    // If not rebinning, just point to our monitor spectra, otherwise create new
    // vectors

    auto Y = (m_commonBins ? monY : Counts(specLength));
    auto E = (m_commonBins ? monE : CountStandardDeviations(specLength));

    if (!m_commonBins) {
      // ConvertUnits can give X vectors of all zeros - skip these, they cause
      // problems
      if (X.back() == 0.0 && X.front() == 0.0)
        continue;
      // Rebin the monitor spectrum to match the binning of the current data
      // spectrum
      VectorHelper::rebinHistogram(
          monX.rawData(), monY.mutableRawData(), monE.mutableRawData(),
          X.rawData(), Y.mutableRawData(), E.mutableRawData(), false);
      // Recalculate the overall normalization factor
      this->normalisationFactor(X, Y, E);
    }

    if (inputEvent) {
      // ----------------------------------- EventWorkspace
      // ---------------------------------------
      EventList &outEL = outputEvent->getSpectrum(i);
      outEL.divide(X.rawData(), Y.mutableRawData(), E.mutableRawData());
    } else {
      // ----------------------------------- Workspace2D
      // ---------------------------------------
      auto &YOut = outputWorkspace->mutableY(i);
      auto &EOut = outputWorkspace->mutableE(i);
      const auto &inY = inputWorkspace->y(i);
      const auto &inE = inputWorkspace->e(i);
      outputWorkspace->mutableX(i) = inputWorkspace->x(i);

      // The code below comes more or less straight out of Divide.cpp
      for (size_t k = 0; k < specLength; ++k) {
        // Get the input Y's
        const double leftY = inY[k];
        const double rightY = Y[k];

        if (rightY == 0.0) {
          hasZeroDivision = true;
        }

        // Calculate result and store in local variable to avoid overwriting
        // original data if
        // output workspace is same as one of the input ones
        const double newY = leftY / rightY;

        if (fabs(rightY) > 1.0e-12 && fabs(newY) > 1.0e-12) {
          const double lhsFactor = (inE[k] < 1.0e-12 || fabs(leftY) < 1.0e-12)
                                       ? 0.0
                                       : pow((inE[k] / leftY), 2);
          const double rhsFactor =
              E[k] < 1.0e-12 ? 0.0 : pow((E[k] / rightY), 2);
          EOut[k] = std::abs(newY) * sqrt(lhsFactor + rhsFactor);
        }

        // Now store the result
//.........这里部分代码省略.........