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

/** Calculate the derivatives for each spectrum in the input workspace
 *
 * @param index :: index of the spectrum
 * @param order :: order of derivatives to calculate
 */
void SplineSmoothing::calculateSpectrumDerivatives(const int index,
        const int order) {
    if (order > 0) {
        API::MatrixWorkspace_sptr derivs =
            setupOutputWorkspace(m_inputWorkspace, order);

        for (int j = 0; j < order; ++j) {
            derivs->setSharedX(j, m_inputWorkspace->sharedX(index));
            calculateDerivatives(*m_inputWorkspacePointData, *derivs, j + 1, index);
        }

        m_derivativeWorkspaceGroup->addWorkspace(derivs);
    }
}

/** Create output workspace
 */
API::MatrixWorkspace_sptr GeneratePeaks::createOutputWorkspace() {
  // Reference workspace and output workspace
  API::MatrixWorkspace_sptr outputWS;

  m_newWSFromParent = true;
  if (!inputWS && binParameters.empty()) {
    // Error! Neither bin parameters or reference workspace is given.
    std::string errmsg(
        "Must define either InputWorkspace or BinningParameters.");
    g_log.error(errmsg);
    throw std::invalid_argument(errmsg);
  } else if (inputWS) {
    // Generate Workspace2D from input workspace
    if (!binParameters.empty())
      g_log.notice()
          << "Both binning parameters and input workspace are given. "
          << "Using input worksapce to generate output workspace!\n";

    outputWS = API::WorkspaceFactory::Instance().create(
        inputWS, inputWS->getNumberHistograms(), inputWS->x(0).size(),
        inputWS->y(0).size());

    // Only copy the X-values from spectra with peaks specified in the table
    // workspace.
    for (const auto &iws : m_spectraSet) {
      outputWS->setSharedX(iws, inputWS->sharedX(iws));
    }

    m_newWSFromParent = true;
  } else {
    // Generate a one-spectrum Workspace2D from binning
    outputWS = createDataWorkspace(binParameters);
    m_newWSFromParent = false;
  }

  return outputWS;
}

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