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

/** Do the initial copy of the data from the input to the output workspace for
 * histogram workspaces.
 *  Takes out the bin width if necessary.
 *  @param inputWS  The input workspace
 *  @param outputWS The output workspace
 */
void ConvertUnitsUsingDetectorTable::fillOutputHist(
    const API::MatrixWorkspace_const_sptr inputWS,
    const API::MatrixWorkspace_sptr outputWS) {
  const int size = static_cast<int>(inputWS->blocksize());

  // Loop over the histograms (detector spectra)
  Progress prog(this, 0.0, 0.2, m_numberOfSpectra);
  int64_t numberOfSpectra_i =
      static_cast<int64_t>(m_numberOfSpectra); // cast to make openmp happy
  PARALLEL_FOR2(inputWS, outputWS)
  for (int64_t i = 0; i < numberOfSpectra_i; ++i) {
    PARALLEL_START_INTERUPT_REGION
    // Take the bin width dependency out of the Y & E data
    if (m_distribution) {
      for (int j = 0; j < size; ++j) {
        const double width =
            std::abs(inputWS->dataX(i)[j + 1] - inputWS->dataX(i)[j]);
        outputWS->dataY(i)[j] = inputWS->dataY(i)[j] * width;
        outputWS->dataE(i)[j] = inputWS->dataE(i)[j] * width;
      }
    } else {
      // Just copy over
      outputWS->dataY(i) = inputWS->readY(i);
      outputWS->dataE(i) = inputWS->readE(i);
    }
    // Copy over the X data
    outputWS->setX(i, inputWS->refX(i));

    prog.report("Convert to " + m_outputUnit->unitID());
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION
}

void LoadNexusProcessed::loadBlock(NXDataSetTyped<double> & data, NXDataSetTyped<double> & errors,
    int64_t blocksize, int64_t nchannels, int64_t &hist,int64_t& wsIndex,
    API::MatrixWorkspace_sptr local_workspace)
{
  data.load(static_cast<int>(blocksize),static_cast<int>(hist));
  errors.load(static_cast<int>(blocksize),static_cast<int>(hist));
  double *data_start = data();
  double *data_end = data_start + nchannels;
  double *err_start = errors();
  double *err_end = err_start + nchannels;
  int64_t final(hist + blocksize);
  while( hist < final )
  {
    MantidVec& Y = local_workspace->dataY(wsIndex);
    Y.assign(data_start, data_end);
    data_start += nchannels; data_end += nchannels;
    MantidVec& E = local_workspace->dataE(wsIndex);
    E.assign(err_start, err_end);
    err_start += nchannels; err_end += nchannels;
    local_workspace->setX(wsIndex, m_xbins);
    ++hist;
    ++wsIndex;

  }
}

void ConvertEmptyToTof::setTofInWS(const std::vector<double> &tofAxis,
                                   API::MatrixWorkspace_sptr outputWS) {

  const size_t numberOfSpectra = m_inputWS->getNumberHistograms();
  int64_t numberOfSpectraInt64 =
      static_cast<int64_t>(numberOfSpectra); // cast to make openmp happy

  g_log.debug() << "Setting the TOF X Axis for numberOfSpectra="
                << numberOfSpectra << '\n';

  Progress prog(this, 0.0, 0.2, numberOfSpectra);
  PARALLEL_FOR2(m_inputWS, outputWS)
  for (int64_t i = 0; i < numberOfSpectraInt64; ++i) {
    PARALLEL_START_INTERUPT_REGION
    // Just copy over
    outputWS->dataY(i) = m_inputWS->readY(i);
    outputWS->dataE(i) = m_inputWS->readE(i);
    // copy
    outputWS->setX(i, tofAxis);

    prog.report();
    PARALLEL_END_INTERUPT_REGION
  } // end for i
  PARALLEL_CHECK_INTERUPT_REGION
  outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
}

/**
 * Read spectra from the DAE
 * @param period :: Current period index
 * @param index :: First spectrum index
 * @param count :: Number of spectra to read
 * @param workspace :: Workspace to store the data
 * @param workspaceIndex :: index in workspace to store data
 */
void ISISHistoDataListener::getData(int period, int index, int count,
                                    API::MatrixWorkspace_sptr workspace,
                                    size_t workspaceIndex) {
  const int numberOfBins = m_numberOfBins[m_timeRegime];
  const size_t bufferSize = count * (numberOfBins + 1) * sizeof(int);
  std::vector<int> dataBuffer(bufferSize);
  // Read in spectra from DAE
  int ndims = 2, dims[2];
  dims[0] = count;
  dims[1] = numberOfBins + 1;

  int spectrumIndex = index + period * (m_totalNumberOfSpectra + 1);
  if (IDCgetdat(m_daeHandle, spectrumIndex, count, dataBuffer.data(), dims,
                &ndims) != 0) {
    g_log.error("Unable to read DATA from DAE " + m_daeName);
    throw Kernel::Exception::FileError("Unable to read DATA from DAE ",
                                       m_daeName);
  }

  for (size_t i = 0; i < static_cast<size_t>(count); ++i) {
    size_t wi = workspaceIndex + i;
    workspace->setX(wi, m_bins[m_timeRegime]);
    MantidVec &y = workspace->dataY(wi);
    MantidVec &e = workspace->dataE(wi);
    workspace->getSpectrum(wi)->setSpectrumNo(index + static_cast<specid_t>(i));
    size_t shift = i * (numberOfBins + 1) + 1;
    y.assign(dataBuffer.begin() + shift, dataBuffer.begin() + shift + y.size());
    std::transform(y.begin(), y.end(), e.begin(), dblSqrt);
  }
}

/** Convert the workspace units according to a simple output = a * (input^b) relationship
 *  @param outputWS :: the output workspace
 *  @param factor :: the conversion factor a to apply
 *  @param power :: the Power b to apply to the conversion
 */
void ConvertUnits::convertQuickly(API::MatrixWorkspace_sptr outputWS, const double& factor, const double& power)
{
  Progress prog(this,0.2,1.0,m_numberOfSpectra);
  int64_t numberOfSpectra_i = static_cast<int64_t>(m_numberOfSpectra); // cast to make openmp happy

  // See if the workspace has common bins - if so the X vector can be common
  // First a quick check using the validator
  CommonBinsValidator sameBins;
  bool commonBoundaries = false;
  if ( sameBins.isValid(outputWS) == "" )
  {
    commonBoundaries =  WorkspaceHelpers::commonBoundaries(outputWS);
    // Only do the full check if the quick one passes
    if (commonBoundaries)
    {
      // Calculate the new (common) X values
      MantidVec::iterator iter;
      for (iter = outputWS->dataX(0).begin(); iter != outputWS->dataX(0).end(); ++iter)
      {
        *iter = factor * std::pow(*iter,power);
      }

      MantidVecPtr xVals;
      xVals.access() = outputWS->dataX(0);

      PARALLEL_FOR1(outputWS)
      for (int64_t j = 1; j < numberOfSpectra_i; ++j)
      {
        PARALLEL_START_INTERUPT_REGION
        outputWS->setX(j,xVals);
        prog.report("Convert to " + m_outputUnit->unitID());
        PARALLEL_END_INTERUPT_REGION
      }
      PARALLEL_CHECK_INTERUPT_REGION
      if (!m_inputEvents) // if in event mode the work is done
        return;
    }
  }

//.........这里部分代码省略.........
  size_t i0 = 0;
  // intermediate variables
  double xx = 0.0;
  double xy = 0.0;
  double ym = 0.0;
  // low frequency filter values: the higher the power spectrum the closer the
  // filter to 1.0
  for (size_t i = 0; i < n2; ++i) {
    double cd1 = powerSpec[i] / noise;
    if (cd1 < 1.0 && i > imax) {
      i0 = i;
      break;
    }
    double cd2 = log(cd1);
    wf[i] = cd1 / (1.0 + cd1);
    double j = static_cast<double>(i + 1);
    xx += j * j;
    xy += j * cd2;
    ym += cd2;
  }

  // i0 should always be > 0 but in case something goes wrong make a check
  if (i0 > 0) {
    g_log.debug() << "Noise start index " << i0 << std::endl;

    // high frequency filter values: smooth decreasing function
    double ri0f = static_cast<double>(i0 + 1);
    double xm = (1.0 + ri0f) / 2;
    ym /= ri0f;
    double a1 = (xy - ri0f * xm * ym) / (xx - ri0f * xm * xm);
    double b1 = ym - a1 * xm;

    g_log.debug() << "(a1,b1) = (" << a1 << ',' << b1 << ')' << std::endl;

    const double dblev = -20.0;
    // cut-off index
    double ri1 = floor((dblev / 4 - b1) / a1);
    if (ri1 < static_cast<double>(i0)) {
      g_log.warning() << "Failed to build Wiener filter: no smoothing."
                      << std::endl;
      ri1 = static_cast<double>(i0);
    }
    size_t i1 = static_cast<size_t>(ri1);
    if (i1 > n2)
      i1 = n2;
    for (size_t i = i0; i < i1; ++i) {
      double s = exp(a1 * static_cast<double>(i + 1) + b1);
      wf[i] = s / (1.0 + s);
    }
    // wf[i] for i1 <= i < n2 are 0.0

    g_log.debug() << "Cut-off index " << i1 << std::endl;
  } else {
    g_log.warning() << "Power spectrum has an unexpected shape: no smoothing"
                    << std::endl;
    return copyInput(inputWS, wsIndex);
  }

  // multiply the fourier transform by the filter
  auto &re = fourierOut->dataY(0);
  auto &im = fourierOut->dataY(1);

  std::transform(re.begin(), re.end(), wf.begin(), re.begin(),
                 std::multiplies<double>());
  std::transform(im.begin(), im.end(), wf.begin(), im.begin(),
                 std::multiplies<double>());

  // inverse fourier transform
  fourier = createChildAlgorithm("RealFFT");
  fourier->initialize();
  fourier->setProperty("InputWorkspace", fourierOut);
  fourier->setProperty("IgnoreXBins", true);
  fourier->setPropertyValue("Transform", "Backward");
  fourier->execute();

  API::MatrixWorkspace_sptr out = fourier->getProperty("OutputWorkspace");
  auto &background = fitOut->readY(1);
  auto &y = out->dataY(0);

  if (y.size() != background.size()) {
    throw std::logic_error("Logic error: inconsistent arrays");
  }

  // add the spline "background" to the smoothed data
  std::transform(y.begin(), y.end(), background.begin(), y.begin(),
                 std::plus<double>());

  // copy the x-values and errors from the original spectrum
  // remove the last values for odd-sized inputs
  if (isOddSize) {
    out->dataX(0).assign(X.begin(), X.end() - 1);
    out->dataE(0).assign(E.begin(), E.end() - 1);
    out->dataY(0).resize(Y.size() - 1);
  } else {
    out->setX(0, X);
    out->dataE(0).assign(E.begin(), E.end());
  }

  return out;
}

/** Executes the regroup algorithm
 *
 *  @throw runtime_error Thrown if
 */
void Regroup::exec()
{
  // retrieve the properties
  std::vector<double> rb_params=getProperty("Params");

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

  // can work only if all histograms have the same boundaries
  if (!API::WorkspaceHelpers::commonBoundaries(inputW))
  {
    g_log.error("Histograms with different boundaries");
    throw std::runtime_error("Histograms with different boundaries");
  }

  bool dist = inputW->isDistribution();

  int histnumber = static_cast<int>(inputW->getNumberHistograms());
  MantidVecPtr XValues_new;
  const MantidVec & XValues_old = inputW->readX(0);
  std::vector<int> xoldIndex;// indeces of new x in XValues_old
  // create new output X axis
  int ntcnew = newAxis(rb_params,XValues_old,XValues_new.access(),xoldIndex);

  // make output Workspace the same type is the input, but with new length of signal array
  API::MatrixWorkspace_sptr outputW = API::WorkspaceFactory::Instance().create(inputW,histnumber,ntcnew,ntcnew-1);

  int progress_step = histnumber / 100;
  if (progress_step == 0) progress_step = 1;
  for (int hist=0; hist <  histnumber;hist++)
  {
    // get const references to input Workspace arrays (no copying)
    const MantidVec& XValues = inputW->readX(hist);
    const MantidVec& YValues = inputW->readY(hist);
    const MantidVec& YErrors = inputW->readE(hist);

    //get references to output workspace data (no copying)
    MantidVec& YValues_new=outputW->dataY(hist);
    MantidVec& YErrors_new=outputW->dataE(hist);

    // output data arrays are implicitly filled by function
    rebin(XValues,YValues,YErrors,xoldIndex,YValues_new,YErrors_new, dist);

    outputW->setX(hist,XValues_new);

    if (hist % progress_step == 0)
    {
        progress(double(hist)/histnumber);
        interruption_point();
    }
  }

  outputW->isDistribution(dist);

  // Copy units
  if (outputW->getAxis(0)->unit().get())
    outputW->getAxis(0)->unit() = inputW->getAxis(0)->unit();
  try
  {
    if (inputW->getAxis(1)->unit().get())
      outputW->getAxis(1)->unit() = inputW->getAxis(1)->unit();
   }
  catch(Exception::IndexError) {
    // OK, so this isn't a Workspace2D
  }

  // Assign it to the output workspace property
  setProperty("OutputWorkspace",outputW);

  return;
}