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

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

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

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

//.........这里部分代码省略.........
    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);
  auto &realE = ows->mutableE(0);
  auto &imagE = ows->mutableE(1);

  const auto xPointData = ws->histogram(0).points();
  // First X value
  const double X0 = xPointData.front();

  // calculate exponential decay outside of the loop
  std::vector<double> expDecay = xPointData.rawData();
  std::transform(expDecay.begin(), expDecay.end(), expDecay.begin(),
                 [X0, muLife](double x) { return exp(-(x - X0) / muLife); });

  for (size_t i = 0; i < npoints; i++) {
    for (size_t h = 0; h < nspec; h++) {
      if (!emptySpectrum[h]) {
        // (X,Y,E) with exponential decay removed
        const double X = ws->x(h)[i];
        const double exponential = n0[h] * exp(-(X - X0) / muLife);
        const double Y = ws->y(h)[i] - exponential;
        const double E =
            (ws->y(h)[i] > poissonLimit) ? ws->e(h)[i] : sqrt(exponential);

        realY[i] += aj[h] * Y;
        imagY[i] += bj[h] * Y;
        realE[i] += aj[h] * aj[h] * E * E;
        imagE[i] += bj[h] * bj[h] * E * E;
      }
    }
    realE[i] = sqrt(realE[i]);
    imagE[i] = sqrt(imagE[i]);

    // Regain exponential decay
    realY[i] /= expDecay[i];
    imagY[i] /= expDecay[i];
    realE[i] /= expDecay[i];
    imagE[i] /= expDecay[i];
  }

  // New Y axis label
  ows->setYUnit("Asymmetry");

  return ows;
}

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