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

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

/**
 * The main method to calculate the ring profile for 2d image based workspace.
 *
 * 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 elements of the spectrum 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::processNumericImageRingProfile(
    const API::MatrixWorkspace_sptr inputWS, std::vector<double> &output_bins) {
  // allocate the bin positions vector
  std::vector<int> bin_n(inputWS->dataY(0).size(), -1);

  // consider that each spectrum is a row in the image
  for (int i = 0; i < static_cast<int>(inputWS->getNumberHistograms()); i++) {
    m_progress->report("Computing ring bins positions for pixels");
    // get bin for the pixels inside this spectrum
    // for each column of the image
    getBinForPixel(inputWS, i, bin_n);

    // accumulate the values from the spectrum to the target bin
    // each column has it correspondend bin_position inside bin_n
    const MantidVec &refY = inputWS->dataY(i);
    for (size_t j = 0; j < bin_n.size(); j++) {

      // is valid bin? No -> skip
      if (bin_n[j] < 0)
        continue;

      // accumulate the values of this spectrum inside this bin
      output_bins[bin_n[j]] += refY[j];
    }
  }
}

/**
 * Perform a call to nxgetslab, via the NexusClasses wrapped methods for a given blocksize. The xbins are read along with
 * each call to the data/error loading
 * @param data :: The NXDataSet object of y values
 * @param errors :: The NXDataSet object of error values
 * @param xbins :: The xbin NXDataSet
 * @param blocksize :: The blocksize to use
 * @param nchannels :: The number of channels for the block
 * @param hist :: The workspace index to start reading into
 * @param wsIndex :: The workspace index to save data into
 * @param local_workspace :: A pointer to the workspace
 */
void LoadNexusProcessed::loadBlock(NXDataSetTyped<double> & data, NXDataSetTyped<double> & errors, NXDouble & xbins,
    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));
  double *data_start = data();
  double *data_end = data_start + nchannels;
  errors.load(static_cast<int>(blocksize),static_cast<int>(hist));
  double *err_start = errors();
  double *err_end = err_start + nchannels;
  xbins.load(static_cast<int>(blocksize),static_cast<int>(hist));
  const int64_t nxbins(nchannels + 1);
  double *xbin_start = xbins();
  double *xbin_end = xbin_start + nxbins;
  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;
    MantidVec& X = local_workspace->dataX(wsIndex);
    X.assign(xbin_start, xbin_end);
    xbin_start += nxbins; xbin_end += nxbins;
    ++hist;
    ++wsIndex;
  }
}

/** loadData
*  Load the counts data from an NXInt into a workspace
*/
void LoadMuonNexus2::loadData(const Mantid::NeXus::NXInt &counts,
                              const std::vector<double> &timeBins, int wsIndex,
                              int period, int spec,
                              API::MatrixWorkspace_sptr localWorkspace) {
  MantidVec &X = localWorkspace->dataX(wsIndex);
  MantidVec &Y = localWorkspace->dataY(wsIndex);
  MantidVec &E = localWorkspace->dataE(wsIndex);
  X.assign(timeBins.begin(), timeBins.end());

  int nBins = 0;
  int *data = nullptr;

  if (counts.rank() == 3) {
    nBins = counts.dim2();
    data = &counts(period, spec, 0);
  } else if (counts.rank() == 2) {
    nBins = counts.dim1();
    data = &counts(spec, 0);
  } else {
    throw std::runtime_error("Data have unsupported dimansionality");
  }
  assert(nBins + 1 == static_cast<int>(timeBins.size()));

  Y.assign(data, data + nBins);
  typedef double (*uf)(double);
  uf dblSqrt = std::sqrt;
  std::transform(Y.begin(), Y.end(), E.begin(), dblSqrt);
}

/** Smoothing by zeroing.
 *  @param n :: The order of truncation
 *  @param unfilteredWS :: workspace for storing the unfiltered Fourier
 * transform of the input spectrum
 *  @param filteredWS :: workspace for storing the filtered spectrum
 */
void FFTSmooth2::zero(int n, API::MatrixWorkspace_sptr &unfilteredWS,
                      API::MatrixWorkspace_sptr &filteredWS) {
  int mx = static_cast<int>(unfilteredWS->readX(0).size());
  int my = static_cast<int>(unfilteredWS->readY(0).size());
  int ny = my / n;

  if (ny == 0)
    ny = 1;

  filteredWS =
      API::WorkspaceFactory::Instance().create(unfilteredWS, 2, mx, my);

  const Mantid::MantidVec &Yr = unfilteredWS->readY(0);
  const Mantid::MantidVec &Yi = unfilteredWS->readY(1);
  const Mantid::MantidVec &X = unfilteredWS->readX(0);

  Mantid::MantidVec &yr = filteredWS->dataY(0);
  Mantid::MantidVec &yi = filteredWS->dataY(1);
  Mantid::MantidVec &xr = filteredWS->dataX(0);
  Mantid::MantidVec &xi = filteredWS->dataX(1);

  xr.assign(X.begin(), X.end());
  xi.assign(X.begin(), X.end());
  yr.assign(Yr.size(), 0);
  yi.assign(Yr.size(), 0);

  for (int i = 0; i < ny; i++) {
    yr[i] = Yr[i];
    yi[i] = Yi[i];
  }
}

void RadiusSum::setUpOutputWorkspace(std::vector<double> &values) {

  g_log.debug() << "Output calculated, setting up the output workspace\n";

  API::MatrixWorkspace_sptr outputWS = API::WorkspaceFactory::Instance().create(
      inputWS, 1, values.size() + 1, values.size());

  g_log.debug() << "Set the data\n";
  MantidVec &refY = outputWS->dataY(0);
  std::copy(values.begin(), values.end(), refY.begin());

  g_log.debug() << "Set the bins limits\n";
  MantidVec &refX = outputWS->dataX(0);
  double bin_size = (max_radius - min_radius) / num_bins;

  for (int i = 0; i < (static_cast<int>(refX.size())) - 1; i++)
    refX[i] = min_radius + i * bin_size;
  refX.back() = max_radius;

  // configure the axis:
  // for numeric images, the axis are the same as the input workspace, and are
  // copied in the creation.

  // for instrument related, the axis Y (1) continues to be the same.
  // it is necessary to change only the axis X. We have to change it to radius.
  if (inputWorkspaceHasInstrumentAssociated(inputWS)) {
    API::Axis *const horizontal = new API::NumericAxis(refX.size());
    auto labelX = UnitFactory::Instance().create("Label");
    boost::dynamic_pointer_cast<Units::Label>(labelX)->setLabel("Radius");
    horizontal->unit() = labelX;
    outputWS->replaceAxis(0, horizontal);
  }

  setProperty("OutputWorkspace", outputWS);
}

/** Smoothing using Butterworth filter.
 *  @param n ::     The cutoff frequency control parameter.
 *               Cutoff frequency = my/n where my is the
 *               number of sample points in the data.
 *               As with the "Zeroing" case, the cutoff
 *               frequency is truncated to an integer value
 *               and set to 1 if the truncated value was zero.
 *  @param order :: The order of the Butterworth filter, 1, 2, etc.
 *               This must be a positive integer.
 *  @param unfilteredWS :: workspace for storing the unfiltered Fourier
 * transform of the input spectrum
 *  @param filteredWS :: workspace for storing the filtered spectrum
 */
void FFTSmooth2::Butterworth(int n, int order,
                             API::MatrixWorkspace_sptr &unfilteredWS,
                             API::MatrixWorkspace_sptr &filteredWS) {
  int mx = static_cast<int>(unfilteredWS->readX(0).size());
  int my = static_cast<int>(unfilteredWS->readY(0).size());
  int ny = my / n;

  if (ny == 0)
    ny = 1;

  filteredWS =
      API::WorkspaceFactory::Instance().create(unfilteredWS, 2, mx, my);

  const Mantid::MantidVec &Yr = unfilteredWS->readY(0);
  const Mantid::MantidVec &Yi = unfilteredWS->readY(1);
  const Mantid::MantidVec &X = unfilteredWS->readX(0);

  Mantid::MantidVec &yr = filteredWS->dataY(0);
  Mantid::MantidVec &yi = filteredWS->dataY(1);
  Mantid::MantidVec &xr = filteredWS->dataX(0);
  Mantid::MantidVec &xi = filteredWS->dataX(1);

  xr.assign(X.begin(), X.end());
  xi.assign(X.begin(), X.end());
  yr.assign(Yr.size(), 0);
  yi.assign(Yr.size(), 0);

  double cutoff = ny;

  for (int i = 0; i < my; i++) {
    double scale = 1.0 / (1.0 + pow(i / cutoff, 2 * order));
    yr[i] = scale * Yr[i];
    yi[i] = scale * Yi[i];
  }
}

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

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

/** Execute the algorithm.
 */
void DampSq::exec()
{
    // TODO Auto-generated execute stub

    // 1. Generate new workspace
    API::MatrixWorkspace_const_sptr isqspace = getProperty("InputWorkspace");
    API::MatrixWorkspace_sptr osqspace = WorkspaceFactory::Instance().create(isqspace, 1, isqspace->size(), isqspace->size());

    int mode = getProperty("Mode");
    double qmax = getProperty("QMax");

    if (mode < 1 || mode > 4) {
        g_log.error("Damp mode can only be 1, 2, 3, or 4");
        return;
    }

    // 2. Get access to all
    const MantidVec& iQVec = isqspace->dataX(0);
    const MantidVec& iSVec = isqspace->dataY(0);
    const MantidVec& iEVec = isqspace->dataE(0);

    MantidVec& oQVec = osqspace->dataX(0);
    MantidVec& oSVec = osqspace->dataY(0);
    MantidVec& oEVec = osqspace->dataE(0);

    // 3. Calculation
    double dqmax = qmax - iQVec[0];

    double damp;
    for (unsigned int i = 0; i < iQVec.size(); i ++) {
        // a) calculate damp coefficient
        switch (mode) {
        case 1:
            damp = dampcoeff1(iQVec[i], qmax, dqmax);
            break;
        case 2:
            damp = dampcoeff2(iQVec[i], qmax, dqmax);;
            break;
        case 3:
            damp = dampcoeff3(iQVec[i], qmax, dqmax);;
            break;
        case 4:
            damp = dampcoeff4(iQVec[i], qmax, dqmax);;
            break;
        default:
            damp = 0;
            break;
        }
        // b) calculate new S(q)
        oQVec[i] = iQVec[i];
        oSVec[i] = 1 + damp*(iSVec[i]-1);
        oEVec[i] = damp*iEVec[i];
    }  // i

    // 4. Over
    setProperty("OutputWorkspace", osqspace);

    return;
}

/** Performs the Holtzer transformation: IQ v Q
 *  @param ws The workspace to be transformed
 */
void IQTransform::holtzer(API::MatrixWorkspace_sptr ws)
{
  MantidVec& X = ws->dataX(0);
  MantidVec& Y = ws->dataY(0);
  MantidVec& E = ws->dataE(0);
  std::transform(Y.begin(),Y.end(),X.begin(),Y.begin(),std::multiplies<double>());
  std::transform(E.begin(),E.end(),X.begin(),E.begin(),std::multiplies<double>());

  ws->setYUnitLabel("I x Q");
}

/** Performs the Porod transformation: IQ^4 v Q
 *  @param ws The workspace to be transformed
 */
void IQTransform::porod(API::MatrixWorkspace_sptr ws)
{
  MantidVec& X = ws->dataX(0);
  MantidVec& Y = ws->dataY(0);
  MantidVec& E = ws->dataE(0);
  MantidVec Q4(X.size());
  std::transform(X.begin(),X.end(),X.begin(),Q4.begin(),VectorHelper::TimesSquares<double>());
  std::transform(Y.begin(),Y.end(),Q4.begin(),Y.begin(),std::multiplies<double>());
  std::transform(E.begin(),E.end(),Q4.begin(),E.begin(),std::multiplies<double>());

  ws->setYUnitLabel("I x Q^4");
}

/** Performs a log-log transformation: Ln(I) v Ln(Q)
 *  @param ws The workspace to be transformed
 *  @throw std::range_error if an attempt is made to take log of a negative number
 */
void IQTransform::logLog(API::MatrixWorkspace_sptr ws)
{
  MantidVec& X = ws->dataX(0);
  MantidVec& Y = ws->dataY(0);
  MantidVec& E = ws->dataE(0);
  std::transform(X.begin(),X.end(),X.begin(),VectorHelper::Log<double>());
  std::transform(E.begin(),E.end(),Y.begin(),E.begin(),std::divides<double>());
  std::transform(Y.begin(),Y.end(),Y.begin(),VectorHelper::LogNoThrow<double>());

  ws->setYUnitLabel("Ln(I)");
  m_label->setLabel("Ln(Q)");
}

/** Performs the Guinier (sheets) transformation: Ln(IQ^2) v Q^2
 *  @param ws The workspace to be transformed
 *  @throw std::range_error if an attempt is made to take log of a negative number
 */
void IQTransform::guinierSheets(API::MatrixWorkspace_sptr ws)
{
  MantidVec& X = ws->dataX(0);
  MantidVec& Y = ws->dataY(0);
  MantidVec& E = ws->dataE(0);
  std::transform(E.begin(),E.end(),Y.begin(),E.begin(),std::divides<double>());
  std::transform(X.begin(),X.end(),X.begin(),VectorHelper::Squares<double>());
  std::transform(Y.begin(),Y.end(),X.begin(),Y.begin(),std::multiplies<double>());
  std::transform(Y.begin(),Y.end(),Y.begin(),VectorHelper::LogNoThrow<double>());

  ws->setYUnitLabel("Ln(I x Q^2)");
  m_label->setLabel("Q^2");
}

/** Remove background per pixel
 * @brief ConvertCWSDExpToMomentum::removeBackground
 * @param dataws
 */
void ConvertCWSDExpToMomentum::removeBackground(
    API::MatrixWorkspace_sptr dataws) {
  if (dataws->getNumberHistograms() != m_backgroundWS->getNumberHistograms())
    throw std::runtime_error("Impossible to have this situation");

  size_t numhist = dataws->getNumberHistograms();
  for (size_t i = 0; i < numhist; ++i) {
    double bkgd_y = m_backgroundWS->readY(i)[0];
    if (fabs(bkgd_y) > 1.E-2) {
      dataws->dataY(i)[0] -= bkgd_y;
      dataws->dataE(i)[0] = std::sqrt(dataws->readY(i)[0]);
    }
  }
}