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

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

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

  int bitDepth = getProperty(PROP_BIT_DEPTH);
  const size_t bytespp = static_cast<size_t>(bitDepth) / 8;

  for (size_t row = 0; row < sizeY; ++row) {
    const auto &dataY = img->readY(row);
    for (size_t col = 0; col < sizeX; ++col) {
      int32_t pixelVal;
      if (8 == bitDepth) {
        pixelVal = static_cast<uint8_t>(dataY[col]);
      } else if (16 == bitDepth) {
        pixelVal = static_cast<uint16_t>(dataY[col]);
      } else if (32 == bitDepth) {
        pixelVal = static_cast<uint32_t>(dataY[col]);
      }

      // change endianness: to sequence of bytes in big-endian
      // this needs revisiting (similarly in LoadFITS)
      // See https://github.com/mantidproject/mantid/pull/15964
      std::array<uint8_t, g_maxBytesPP> bytesPixel;
      uint8_t *iter = reinterpret_cast<uint8_t *>(&pixelVal);
      std::reverse_copy(iter, iter + bytespp, bytesPixel.data());

      file.write(reinterpret_cast<const char *>(bytesPixel.data()), bytespp);
    }
  }
}

/** Reads in the data corresponding to a single spectrum
 *  @param speFile ::   The file handle
 *  @param workspace :: The output workspace
 *  @param index ::     The index of the current spectrum
 */
void LoadSPE::readHistogram(FILE *speFile, API::MatrixWorkspace_sptr workspace,
                            size_t index) {
  // First, there should be a comment line
  char comment[100];
  fgets(comment, 100, speFile);
  if (comment[0] != '#')
    reportFormatError(std::string(comment));

  // Then it's the Y values
  MantidVec &Y = workspace->dataY(index);
  const size_t nbins = workspace->blocksize();
  int retval;
  for (size_t i = 0; i < nbins; ++i) {
    retval = fscanf(speFile, "%10le", &Y[i]);
    // g_log.error() << Y[i] << std::endl;
    if (retval != 1) {
      std::stringstream ss;
      ss << "Reading data value" << i << " of histogram " << index;
      reportFormatError(ss.str());
    }
    // -10^30 is the flag for not a number used in SPE files (from
    // www.mantidproject.org/images/3/3d/Spe_file_format.pdf)
    if (Y[i] == SaveSPE::MASK_FLAG) {
      Y[i] = std::numeric_limits<double>::quiet_NaN();
    }
  }
  // Read to EOL
  fgets(comment, 100, speFile);

  // Another comment line
  fgets(comment, 100, speFile);
  if (comment[0] != '#')
    reportFormatError(std::string(comment));

  // And then the error values
  MantidVec &E = workspace->dataE(index);
  for (size_t i = 0; i < nbins; ++i) {
    retval = fscanf(speFile, "%10le", &E[i]);
    if (retval != 1) {
      std::stringstream ss;
      ss << "Reading error value" << i << " of histogram " << index;
      reportFormatError(ss.str());
    }
  }
  // Read to EOL
  fgets(comment, 100, speFile);

  return;
}

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

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

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

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

  // Only create output workspace if different to input one
  if (outputWorkspace != inputWorkspace )
  {
    if (inputEvent)
    {
      //Make a brand new EventWorkspace
      outputEvent = boost::dynamic_pointer_cast<EventWorkspace>(
          API::WorkspaceFactory::Instance().create("EventWorkspace", inputEvent->getNumberHistograms(), 2, 1));
      //Copy geometry and data
      API::WorkspaceFactory::Instance().initializeFromParent(inputEvent, outputEvent, false);
      outputEvent->copyDataFrom( (*inputEvent) );
      outputWorkspace = boost::dynamic_pointer_cast<MatrixWorkspace>(outputEvent);
    }
    else
      outputWorkspace = WorkspaceFactory::Instance().create(inputWorkspace);
  }

  // Get hold of the monitor spectrum
  const MantidVec& monX = m_monitor->readX(0);
  MantidVec& monY = m_monitor->dataY(0);
  MantidVec& monE = m_monitor->dataE(0);
  // Calculate the overall normalisation just the once if bins are all matching
  if (m_commonBins) this->normalisationFactor(m_monitor->readX(0),&monY,&monE);


  const size_t numHists = inputWorkspace->getNumberHistograms();
  MantidVec::size_type specLength = inputWorkspace->blocksize();
  Progress prog(this,0.0,1.0,numHists);
  // Loop over spectra
  PARALLEL_FOR3(inputWorkspace,outputWorkspace,m_monitor)
  for (int64_t i = 0; i < int64_t(numHists); ++i)
  {
    PARALLEL_START_INTERUPT_REGION
    prog.report();

    const MantidVec& X = inputWorkspace->readX(i);
    // If not rebinning, just point to our monitor spectra, otherwise create new vectors
    MantidVec* Y = ( m_commonBins ? &monY : new MantidVec(specLength) );
    MantidVec* E = ( m_commonBins ? &monE : new MantidVec(specLength) );

    if (!m_commonBins)
    {
      // ConvertUnits can give X vectors of all zeroes - 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,monY,monE,X,*Y,*E,false);
      // Recalculate the overall normalisation factor
      this->normalisationFactor(X,Y,E);
    }

    if (inputEvent)
    {
      // ----------------------------------- EventWorkspace ---------------------------------------
      EventList & outEL = outputEvent->getEventList(i);
      outEL.divide(X, *Y, *E);
    }
    else
    {
      // ----------------------------------- Workspace2D ---------------------------------------
      const MantidVec& inY = inputWorkspace->readY(i);
      const MantidVec& inE = inputWorkspace->readE(i);
      MantidVec& YOut = outputWorkspace->dataY(i);
      MantidVec& EOut = outputWorkspace->dataE(i);
      outputWorkspace->dataX(i) = inputWorkspace->readX(i);
      // The code below comes more or less straight out of Divide.cpp
      for (MantidVec::size_type k = 0; k < specLength; ++k)
      {
        // Get references to the input Y's
        const double& leftY = inY[k];
        const double& rightY = (*Y)[k];

        // 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
        YOut[k] = newY;
      } // end Workspace2D case
    } // end loop over current spectrum

    if (!m_commonBins) { delete Y; delete E; }
    PARALLEL_END_INTERUPT_REGION
  } // end loop over spectra
//.........这里部分代码省略.........

/**
 * Execute smoothing of a single spectrum.
 * @param inputWS :: A workspace to pick a spectrum from.
 * @param wsIndex :: An index of a spectrum to smooth.
 * @return :: A single-spectrum workspace with the smoothed data.
 */
API::MatrixWorkspace_sptr
WienerSmooth::smoothSingleSpectrum(API::MatrixWorkspace_sptr inputWS,
                                   size_t wsIndex) {
  size_t dataSize = inputWS->blocksize();

  // it won't work for very small workspaces
  if (dataSize < 4) {
    g_log.debug() << "No smoothing, spectrum copied." << std::endl;
    return copyInput(inputWS, wsIndex);
  }

  // Due to the way RealFFT works the input should be even-sized
  const bool isOddSize = dataSize % 2 != 0;
  if (isOddSize) {
    // add a fake value to the end to make size even
    inputWS = copyInput(inputWS, wsIndex);
    wsIndex = 0;
    auto &X = inputWS->dataX(wsIndex);
    auto &Y = inputWS->dataY(wsIndex);
    auto &E = inputWS->dataE(wsIndex);
    double dx = X[dataSize - 1] - X[dataSize - 2];
    X.push_back(X.back() + dx);
    Y.push_back(Y.back());
    E.push_back(E.back());
  }

  // the input vectors
  auto &X = inputWS->readX(wsIndex);
  auto &Y = inputWS->readY(wsIndex);
  auto &E = inputWS->readE(wsIndex);

  // Digital fourier transform works best for data oscillating around 0.
  // Fit a spline with a small number of break points to the data.
  // Make sure that the spline passes through the first and the last points
  // of the data.
  // The fitted spline will be subtracted from the data and the difference
  // will be smoothed with the Wiener filter. After that the spline will be
  // added to the smoothed data to produce the output.

  // number of spline break points, must be smaller than the data size but
  // between 2 and 10
  size_t nbreak = 10;
  if (nbreak * 3 > dataSize)
    nbreak = dataSize / 3;

  // NB. The spline mustn't fit too well to the data. If it does smoothing
  // doesn't happen.
  // TODO: it's possible that the spline is unnecessary and a simple linear
  // function will
  //       do a better job.

  g_log.debug() << "Spline break points " << nbreak << std::endl;

  // define the spline
  API::IFunction_sptr spline =
      API::FunctionFactory::Instance().createFunction("BSpline");
  auto xInterval = getStartEnd(X, inputWS->isHistogramData());
  spline->setAttributeValue("StartX", xInterval.first);
  spline->setAttributeValue("EndX", xInterval.second);
  spline->setAttributeValue("NBreak", static_cast<int>(nbreak));
  // fix the first and last parameters to the first and last data values
  spline->setParameter(0, Y.front());
  spline->fix(0);
  size_t lastParamIndex = spline->nParams() - 1;
  spline->setParameter(lastParamIndex, Y.back());
  spline->fix(lastParamIndex);

  // fit the spline to the data
  auto fit = createChildAlgorithm("Fit");
  fit->initialize();
  fit->setProperty("Function", spline);
  fit->setProperty("InputWorkspace", inputWS);
  fit->setProperty("WorkspaceIndex", static_cast<int>(wsIndex));
  fit->setProperty("CreateOutput", true);
  fit->execute();

  // get the fit output workspace; spectrum 2 contains the difference that is to
  // be smoothed
  API::MatrixWorkspace_sptr fitOut = fit->getProperty("OutputWorkspace");

  // Fourier transform the difference spectrum
  auto fourier = createChildAlgorithm("RealFFT");
  fourier->initialize();
  fourier->setProperty("InputWorkspace", fitOut);
  fourier->setProperty("WorkspaceIndex", 2);
  // we don't require bin linearity as we don't need the exact transform
  fourier->setProperty("IgnoreXBins", true);
  fourier->execute();

  API::MatrixWorkspace_sptr fourierOut =
      fourier->getProperty("OutputWorkspace");

  // spectrum 2 of the transformed workspace has the transform modulus which is
  // a square
//.........这里部分代码省略.........

/** Forms the quadrature phase signal (squashogram)
* @param ws :: [input] workspace containing the measured spectra
* @param phase :: [input] table workspace containing the detector phases
* @param n0 :: [input] vector containing the normalization constants
* @return :: workspace containing the quadrature phase signal
*/
API::MatrixWorkspace_sptr
PhaseQuadMuon::squash(const API::MatrixWorkspace_sptr &ws,
                      const API::ITableWorkspace_sptr &phase,
                      const std::vector<double> &n0) {

  // Poisson limit: below this number we consider we don't have enough
  // statistics
  // to apply sqrt(N). This is an arbitrary number used in the original code
  // provided by scientists
  double poissonLimit = 30.;

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

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

  if (n0.size() != nspec) {
    throw std::invalid_argument("Invalid normalization constants");
  }

  // Get the maximum asymmetry
  double maxAsym = 0.;
  for (size_t h = 0; h < nspec; h++) {
    if (phase->Double(h, 1) > maxAsym) {
      maxAsym = phase->Double(h, 1);
    }
  }
  if (maxAsym == 0.0) {
    throw std::invalid_argument("Invalid detector asymmetries");
  }

  std::vector<double> aj, bj;
  {
    // Calculate coefficients aj, bj

    double sxx = 0;
    double syy = 0;
    double sxy = 0;

    for (size_t h = 0; h < nspec; h++) {
      double asym = phase->Double(h, 1) / maxAsym;
      double phi = phase->Double(h, 2);
      double X = n0[h] * asym * cos(phi);
      double Y = n0[h] * asym * sin(phi);
      sxx += X * X;
      syy += Y * Y;
      sxy += X * Y;
    }

    double lam1 = 2 * syy / (sxx * syy - sxy * sxy);
    double mu1 = 2 * sxy / (sxy * sxy - sxx * syy);
    double lam2 = 2 * sxy / (sxy * sxy - sxx * syy);
    double mu2 = 2 * sxx / (sxx * syy - sxy * sxy);
    for (size_t h = 0; h < nspec; h++) {
      double asym = phase->Double(h, 1) / maxAsym;
      double phi = phase->Double(h, 2);
      double X = n0[h] * asym * cos(phi);
      double Y = n0[h] * asym * sin(phi);
      aj.push_back((lam1 * X + mu1 * Y) * 0.5);
      bj.push_back((lam2 * X + mu2 * Y) * 0.5);
    }
  }

  // First X value
  double X0 = ws->x(0).front();

  // Create and populate output workspace
  API::MatrixWorkspace_sptr ows = API::WorkspaceFactory::Instance().create(
      "Workspace2D", 2, npoints + 1, npoints);

  // X
  ows->setSharedX(0, ws->sharedX(0));
  ows->setSharedX(1, ws->sharedX(0));

  // Phase quadrature
  auto &realY = ows->mutableY(0);
  auto &imagY = ows->mutableY(1);
  auto &realE = ows->mutableE(0);
  auto &imagE = ows->mutableE(1);

  for (size_t i = 0; i < npoints; i++) {
    for (size_t h = 0; h < nspec; h++) {

      // (X,Y,E) with exponential decay removed
      const double X = ws->x(h)[i];
      const double Y = ws->y(h)[i] - n0[h] * exp(-(X - X0) / muLife);
      const double E = (ws->y(h)[i] > poissonLimit)
                           ? ws->e(h)[i]
                           : sqrt(n0[h] * exp(-(X - X0) / muLife));

      realY[i] += aj[h] * Y;
      imagY[i] += bj[h] * Y;
      realE[i] += aj[h] * aj[h] * E * E;
//.........这里部分代码省略.........

/** Forms the quadrature phase signal (squashogram)
 * @param ws :: [input] workspace containing the measured spectra
 * @param phase :: [input] table workspace containing the detector phases
 * @param n0 :: [input] vector containing the normalization constants
 * @return :: workspace containing the quadrature phase signal
 */
API::MatrixWorkspace_sptr
PhaseQuadMuon::squash(const API::MatrixWorkspace_sptr &ws,
                      const API::ITableWorkspace_sptr &phase,
                      const std::vector<double> &n0) {

  // Poisson limit: below this number we consider we don't have enough
  // statistics
  // to apply sqrt(N). This is an arbitrary number used in the original code
  // provided by scientists
  const double poissonLimit = 30.;

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

  const size_t nspec = ws->getNumberHistograms();

  if (n0.size() != nspec) {
    throw std::invalid_argument("Invalid normalization constants");
  }

  auto names = phase->getColumnNames();
  for (auto &name : names) {
    std::transform(name.begin(), name.end(), name.begin(), ::tolower);
  }
  auto phaseIndex = findName(phaseNames, names);
  auto asymmetryIndex = findName(asymmNames, names);

  // Get the maximum asymmetry
  double maxAsym = 0.;
  for (size_t h = 0; h < nspec; h++) {
    if (phase->Double(h, asymmetryIndex) > maxAsym &&
        phase->Double(h, asymmetryIndex) != ASYMM_ERROR) {
      maxAsym = phase->Double(h, asymmetryIndex);
    }
  }

  if (maxAsym == 0.0) {
    throw std::invalid_argument("Invalid detector asymmetries");
  }
  std::vector<bool> emptySpectrum;
  emptySpectrum.reserve(nspec);
  std::vector<double> aj, bj;
  {
    // Calculate coefficients aj, bj

    double sxx = 0.;
    double syy = 0.;
    double sxy = 0.;
    for (size_t h = 0; h < nspec; h++) {
      emptySpectrum.push_back(
          std::all_of(ws->y(h).begin(), ws->y(h).end(),
                      [](double value) { return value == 0.; }));
      if (!emptySpectrum[h]) {
        const double asym = phase->Double(h, asymmetryIndex) / maxAsym;
        const double phi = phase->Double(h, phaseIndex);
        const double X = n0[h] * asym * cos(phi);
        const double Y = n0[h] * asym * sin(phi);
        sxx += X * X;
        syy += Y * Y;
        sxy += X * Y;
      }
    }

    const double lam1 = 2 * syy / (sxx * syy - sxy * sxy);
    const double mu1 = 2 * sxy / (sxy * sxy - sxx * syy);
    const double lam2 = 2 * sxy / (sxy * sxy - sxx * syy);
    const double mu2 = 2 * sxx / (sxx * syy - sxy * sxy);
    for (size_t h = 0; h < nspec; h++) {
      if (emptySpectrum[h]) {
        aj.push_back(0.0);
        bj.push_back(0.0);
      } else {
        const double asym = phase->Double(h, asymmetryIndex) / maxAsym;
        const double phi = phase->Double(h, phaseIndex);
        const double X = n0[h] * asym * cos(phi);
        const double Y = n0[h] * asym * sin(phi);
        aj.push_back((lam1 * X + mu1 * Y) * 0.5);
        bj.push_back((lam2 * X + mu2 * Y) * 0.5);
      }
    }
  }

  const size_t npoints = ws->blocksize();
  // Create and populate output workspace
  API::MatrixWorkspace_sptr ows =
      API::WorkspaceFactory::Instance().create(ws, 2, npoints + 1, npoints);

  // X
  ows->setSharedX(0, ws->sharedX(0));
  ows->setSharedX(1, ws->sharedX(0));

  // Phase quadrature
  auto &realY = ows->mutableY(0);
  auto &imagY = ows->mutableY(1);
//.........这里部分代码省略.........

/** 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
//.........这里部分代码省略.........

/** Executes the algorithm
 *
 *  @throw Exception::FileError If the grouping file cannot be opened or read successfully
 *  @throw runtime_error If unable to run one of the sub-algorithms successfully
 */
void DiffractionFocussing::exec()
{
  // retrieve the properties
  std::string groupingFileName=getProperty("GroupingFileName");

  // Get the input workspace
  MatrixWorkspace_sptr inputW = getProperty("InputWorkspace");

  bool dist = inputW->isDistribution();

  //do this first to check that a valid file is available before doing any work
  std::multimap<int64_t,int64_t> detectorGroups;// <group, UDET>
  if (!readGroupingFile(groupingFileName, detectorGroups))
  {
    throw Exception::FileError("Error reading .cal file",groupingFileName);
  }

  //Convert to d-spacing units
  API::MatrixWorkspace_sptr tmpW = convertUnitsToDSpacing(inputW);

  //Rebin to a common set of bins
  RebinWorkspace(tmpW);

  std::set<int64_t> groupNumbers;
  for(std::multimap<int64_t,int64_t>::const_iterator d = detectorGroups.begin();d!=detectorGroups.end();d++)
  {
    if (groupNumbers.find(d->first) == groupNumbers.end())
    {
      groupNumbers.insert(d->first);
    }
  }

  int iprogress = 0;
  int iprogress_count = static_cast<int>(groupNumbers.size());
  int iprogress_step = iprogress_count / 100;
  if (iprogress_step == 0) iprogress_step = 1;
  std::vector<int64_t> resultIndeces;
  for(std::set<int64_t>::const_iterator g = groupNumbers.begin();g!=groupNumbers.end();g++)
  {
    if (iprogress++ % iprogress_step == 0)
    {
      progress(0.68 + double(iprogress)/iprogress_count/3);
    }
    std::multimap<int64_t,int64_t>::const_iterator from = detectorGroups.lower_bound(*g);
    std::multimap<int64_t,int64_t>::const_iterator to =   detectorGroups.upper_bound(*g);
    std::vector<detid_t> detectorList;
    for(std::multimap<int64_t,int64_t>::const_iterator d = from;d!=to;d++)
      detectorList.push_back(static_cast<detid_t>(d->second));
    // Want version 1 of GroupDetectors here
    API::IAlgorithm_sptr childAlg = createSubAlgorithm("GroupDetectors",-1.0,-1.0,true,1);
    childAlg->setProperty("Workspace", tmpW);
    childAlg->setProperty< std::vector<detid_t> >("DetectorList",detectorList);
    childAlg->executeAsSubAlg();
    try
    {
      // get the index of the combined spectrum
      int ri = childAlg->getProperty("ResultIndex");
      if (ri >= 0)
      {
        resultIndeces.push_back(ri);
      }
    }
    catch(...)
    {
      throw std::runtime_error("Unable to get Properties from GroupDetectors sub-algorithm");
    }
  }

  // Discard left-over spectra, but print warning message giving number discarded
  int discarded = 0;
  const int64_t oldHistNumber = tmpW->getNumberHistograms();
  API::Axis *spectraAxis = tmpW->getAxis(1);
  for(int64_t i=0; i < oldHistNumber; i++)
    if ( spectraAxis->spectraNo(i) >= 0 && find(resultIndeces.begin(),resultIndeces.end(),i) == resultIndeces.end())
    {
      ++discarded;
    }
  g_log.warning() << "Discarded " << discarded << " spectra that were not assigned to any group" << std::endl;

  // Running GroupDetectors leads to a load of redundant spectra
  // Create a new workspace that's the right size for the meaningful spectra and copy them in
  int64_t newSize = tmpW->blocksize();
  API::MatrixWorkspace_sptr outputW = API::WorkspaceFactory::Instance().create(tmpW,resultIndeces.size(),newSize+1,newSize);
  // Copy units
  outputW->getAxis(0)->unit() = tmpW->getAxis(0)->unit();
  outputW->getAxis(1)->unit() = tmpW->getAxis(1)->unit();

  API::Axis *spectraAxisNew = outputW->getAxis(1);

  for(int64_t hist=0; hist < static_cast<int64_t>(resultIndeces.size()); hist++)
  {
    int64_t i = resultIndeces[hist];
    double spNo = static_cast<double>(spectraAxis->spectraNo(i));
    MantidVec &tmpE = tmpW->dataE(i);
    MantidVec &outE = outputW->dataE(hist);
//.........这里部分代码省略.........