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

/** Calculate the derivatives of the given order from the interpolated points
 *
 * @param inputWorkspace :: The input workspace
 * @param outputWorkspace :: The output workspace
 * @param order :: The order of derivatives to calculate
 */
void SplineInterpolation::calculateDerivatives(
    API::MatrixWorkspace_const_sptr inputWorkspace,
    API::MatrixWorkspace_sptr outputWorkspace, int order) const {
  // get x and y parameters from workspaces
  const size_t nData = inputWorkspace->y(0).size();
  const double *xValues = &(inputWorkspace->x(0)[0]);
  double *yValues = &(outputWorkspace->mutableY(order - 1)[0]);

  // calculate the derivatives
  m_cspline->derivative1D(yValues, xValues, nData, order);
}

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

/** 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
  auto &Y = workspace->mutableY(index);
  const size_t nbins = workspace->blocksize();
  int retval;
  for (size_t i = 0; i < nbins; ++i) {
    retval = fscanf(speFile, "%10le", &Y[i]);
    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
  auto &E = workspace->mutableE(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);
}

/** Zeroes all data points outside the X values given
 *  @param outputWorkspace :: The output workspace - data has already been
 * copied
 *  @param inIndex ::         The workspace index of the spectrum in the input
 * workspace
 *  @param outIndex ::        The workspace index of the spectrum in the output
 * workspace
 */
void ExtractSpectra::cropRagged(API::MatrixWorkspace_sptr outputWorkspace,
                                int inIndex, int outIndex) {
  auto &Y = outputWorkspace->mutableY(outIndex);
  auto &E = outputWorkspace->mutableE(outIndex);
  const size_t size = Y.size();
  size_t startX = this->getXMin(inIndex);
  if (startX > size)
    startX = size;
  for (size_t i = 0; i < startX; ++i) {
    Y[i] = 0.0;
    E[i] = 0.0;
  }
  size_t endX = this->getXMax(inIndex);
  if (endX > 0)
    endX -= m_histogram;
  for (size_t i = endX; i < size; ++i) {
    Y[i] = 0.0;
    E[i] = 0.0;
  }
}

/** Generate peaks in the given output workspace
  * @param functionmap :: map to contain the list of functions with key as their
 * spectra
  * @param dataWS :: output matrix workspace
  */
void GeneratePeaks::generatePeaks(
    const std::map<specnum_t,
                   std::vector<std::pair<double, API::IFunction_sptr>>> &
        functionmap,
    API::MatrixWorkspace_sptr dataWS) {
  // Calcualte function
  std::map<specnum_t,
           std::vector<std::pair<double, API::IFunction_sptr>>>::const_iterator
      mapiter;
  for (mapiter = functionmap.begin(); mapiter != functionmap.end(); ++mapiter) {
    // Get spec id and translated to wsindex in the output workspace
    specnum_t specid = mapiter->first;
    specnum_t wsindex;
    if (m_newWSFromParent)
      wsindex = specid;
    else
      wsindex = m_SpectrumMap[specid];

    const std::vector<std::pair<double, API::IFunction_sptr>> &vec_centrefunc =
        mapiter->second;
    size_t numpeaksinspec = mapiter->second.size();

    for (size_t ipeak = 0; ipeak < numpeaksinspec; ++ipeak) {
      const std::pair<double, API::IFunction_sptr> &centrefunc =
          vec_centrefunc[ipeak];

      // Determine boundary
      API::IPeakFunction_sptr thispeak = getPeakFunction(centrefunc.second);
      double centre = centrefunc.first;
      double fwhm = thispeak->fwhm();

      //
      const auto &X = dataWS->x(wsindex);
      double leftbound = centre - m_numPeakWidth * fwhm;
      if (ipeak > 0) {
        // Not left most peak.
        API::IPeakFunction_sptr leftPeak =
            getPeakFunction(vec_centrefunc[ipeak - 1].second);

        double middle = 0.5 * (centre + leftPeak->centre());
        if (leftbound < middle)
          leftbound = middle;
      }
      auto left = std::lower_bound(X.cbegin(), X.cend(), leftbound);
      if (left == X.end())
        left = X.begin();

      double rightbound = centre + m_numPeakWidth * fwhm;
      if (ipeak != numpeaksinspec - 1) {
        // Not the rightmost peak
        IPeakFunction_sptr rightPeak =
            getPeakFunction(vec_centrefunc[ipeak + 1].second);
        double middle = 0.5 * (centre + rightPeak->centre());
        if (rightbound > middle)
          rightbound = middle;
      }
      auto right = std::lower_bound(left + 1, X.cend(), rightbound);

      // Build domain & function
      API::FunctionDomain1DVector domain(left,
                                         right); // dataWS->dataX(wsindex));

      // Evaluate the function
      API::FunctionValues values(domain);
      centrefunc.second->function(domain, values);

      // Put to output
      std::size_t offset = (left - X.begin());
      std::size_t numY = values.size();

      auto &dataY = dataWS->mutableY(wsindex);
      for (std::size_t i = 0; i < numY; i++) {
        dataY[i + offset] += values[i];
      }

    } // ENDFOR(ipeak)
  }
}

API::MatrixWorkspace_sptr
CreateFloodWorkspace::removeBackground(API::MatrixWorkspace_sptr ws) {
  g_log.information() << "Remove background "
                      << getPropertyValue(Prop::BACKGROUND) << '\n';
  auto fitWS = transpose(ws);
  auto const &x = fitWS->x(0);

  // Define the fitting interval
  double startX = getProperty(Prop::START_X);
  double endX = getProperty(Prop::END_X);
  std::vector<double> excludeFromFit;
  if (isDefault(Prop::START_X)) {
    startX = x.front();
  } else {
    excludeFromFit.push_back(x.front());
    excludeFromFit.push_back(startX);
  }
  if (isDefault(Prop::END_X)) {
    endX = x.back();
  } else {
    excludeFromFit.push_back(endX);
    excludeFromFit.push_back(x.back());
  }

  // Exclude any bad detectors.
  for (auto i : m_excludedSpectra) {
    excludeFromFit.push_back(i);
    excludeFromFit.push_back(i);
  }

  std::string const function = getBackgroundFunction();

  // Fit the data to determine unwanted background
  auto alg = createChildAlgorithm("Fit", 0.9, 0.99);
  alg->setProperty("Function", function);
  alg->setProperty("InputWorkspace", fitWS);
  alg->setProperty("WorkspaceIndex", 0);
  if (!excludeFromFit.empty()) {
    alg->setProperty("Exclude", excludeFromFit);
  }
  alg->setProperty("Output", "fit");
  alg->execute();

  IFunction_sptr func = alg->getProperty("Function");
  g_log.information() << "Background function parameters:\n";
  for (size_t i = 0; i < func->nParams(); ++i) {
    g_log.information() << "    " << func->parameterName(i) << ": "
                        << func->getParameter(i) << '\n';
  }

  // Divide the workspace by the fitted curve to remove the background
  // and scale to values around 1
  MatrixWorkspace_sptr bkgWS = alg->getProperty("OutputWorkspace");
  auto const &bkg = bkgWS->y(1);
  auto const nHisto = static_cast<int>(ws->getNumberHistograms());
  PARALLEL_FOR_IF(Kernel::threadSafe(*ws, *bkgWS))
  for (int i = 0; i < nHisto; ++i) {
    PARALLEL_START_INTERUPT_REGION
    auto const xVal = x[i];
    if (isExcludedSpectrum(xVal)) {
      ws->mutableY(i)[0] = VERY_BIG_VALUE;
      ws->mutableE(i)[0] = 0.0;
    } else if (xVal >= startX && xVal <= endX) {
      auto const background = bkg[i];
      if (background <= 0.0) {
        throw std::runtime_error(
            "Background is expected to be positive, found value " +
            std::to_string(background) + " at spectrum with workspace index " +
            std::to_string(i));
      }
      ws->mutableY(i)[0] /= background;
      ws->mutableE(i)[0] /= background;
    } else {
      ws->mutableY(i)[0] = 1.0;
      ws->mutableE(i)[0] = 0.0;
    }
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION

  // Remove the logs
  ws->setSharedRun(make_cow<Run>());

  return ws;
}

/**
 * Takes a single valued histogram workspace and assesses which histograms are
 * within the limits.
 * Those that are not are masked on the input workspace.
 * @param countsWS :: Input/Output Integrated workspace to diagnose.
 * @param medianvec The median value calculated from the current counts.
 * @param indexmap Index map.
 * @param maskWS :: A mask workspace to apply.
 * @return The number of detectors that failed the tests, not including those
 * skipped.
 */
int MedianDetectorTest::doDetectorTests(
    const API::MatrixWorkspace_sptr countsWS,
    const std::vector<double> medianvec,
    std::vector<std::vector<size_t>> indexmap,
    API::MatrixWorkspace_sptr maskWS) {
  g_log.debug("Applying the criteria to find failing detectors");

  // A spectra can't fail if the statistics show its value is consistent with
  // the mean value,
  // check the error and how many errorbars we are away
  const double minSigma = getProperty("SignificanceTest");

  // prepare to report progress
  const int numSpec(m_maxWsIndex - m_minWsIndex);
  const int progStep = static_cast<int>(ceil(numSpec / 30.0));
  int steps(0);

  const double deadValue(1.0);
  int numFailed(0);

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

  PARALLEL_FOR_IF(Kernel::threadSafe(*countsWS, *maskWS))
  for (int j = 0; j < static_cast<int>(indexmap.size()); ++j) {
    std::vector<size_t> hists = indexmap.at(j);
    double median = medianvec.at(j);
    const size_t nhist = hists.size();
    g_log.debug() << "new component with " << nhist << " spectra.\n";
    for (size_t i = 0; i < nhist; ++i) {
      g_log.debug() << "Counts workspace index=" << i
                    << ", Mask workspace index=" << hists.at(i) << '\n';
      PARALLEL_START_INTERUPT_REGION
      ++steps;
      // update the progressbar information
      if (steps % progStep == 0) {
        progress(advanceProgress(progStep * static_cast<double>(RTMarkDetects) /
                                 numSpec));
      }

      if (checkForMask) {
        const auto &detids =
            countsWS->getSpectrum(hists.at(i)).getDetectorIDs();
        if (instrument->isDetectorMasked(detids)) {
          maskWS->mutableY(hists.at(i))[0] = deadValue;
          continue;
        }
        if (instrument->isMonitor(detids)) {
          // Don't include in calculation but don't mask it
          continue;
        }
      }

      const double signal = countsWS->y(hists.at(i))[0];

      // Mask out NaN and infinite
      if (!std::isfinite(signal)) {
        maskWS->mutableY(hists.at(i))[0] = deadValue;
        PARALLEL_ATOMIC
        ++numFailed;
        continue;
      }

      const double error = minSigma * countsWS->e(hists.at(i))[0];

      if ((signal < median * m_loFrac && (signal - median < -error)) ||
          (signal > median * m_hiFrac && (signal - median > error))) {
        maskWS->mutableY(hists.at(i))[0] = deadValue;
        PARALLEL_ATOMIC
        ++numFailed;
      }

      PARALLEL_END_INTERUPT_REGION
    }
    PARALLEL_CHECK_INTERUPT_REGION

    // Log finds
    g_log.information() << numFailed << " spectra failed the median tests.\n";
  }
  return numFailed;
}

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