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

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

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

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

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

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

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

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

/** Performs a transformation of the form: Q^A x I^B x Ln(Q^C x I^D x E) v Q^F x
 * I^G x Ln(Q^H x I^I x J).
 *  Uses the 'GeneralFunctionConstants' property where A-J are the 10 (ordered)
 * input constants.
 *  @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::general(API::MatrixWorkspace_sptr ws) {
  MantidVec &X = ws->dataX(0);
  MantidVec &Y = ws->dataY(0);
  MantidVec &E = ws->dataE(0);
  const std::vector<double> C = getProperty("GeneralFunctionConstants");
  // Check for the correct number of elements
  if (C.size() != 10) {
    std::string mess(
        "The General transformation requires 10 values to be provided.");
    g_log.error(mess);
    throw std::invalid_argument(mess);
  }

  for (size_t i = 0; i < Y.size(); ++i) {
    double tmpX = std::pow(X[i], C[7]) * std::pow(Y[i], C[8]) * C[9];
    if (tmpX <= 0.0)
      throw std::range_error(
          "Attempt to take log of a zero or negative number.");
    tmpX = std::pow(X[i], C[5]) * std::pow(Y[i], C[6]) * std::log(tmpX);
    const double tmpY = std::pow(X[i], C[2]) * std::pow(Y[i], C[3]) * C[4];
    if (tmpY <= 0.0)
      throw std::range_error(
          "Attempt to take log of a zero or negative number.");
    const double newY =
        std::pow(X[i], C[0]) * std::pow(Y[i], C[1]) * std::log(tmpY);

    E[i] *= std::pow(X[i], C[0]) *
            (C[1] * std::pow(Y[i], C[1] - 1) * std::log(tmpY) +
             ((std::pow(Y[i], C[1]) * std::pow(X[i], C[2]) * C[4] * C[3] *
               std::pow(Y[i], C[3] - 1)) /
              tmpY));
    X[i] = tmpX;
    Y[i] = newY;
  }

  std::stringstream ylabel;
  ylabel << "Q^" << C[0] << " x I^" << C[1] << " x Ln( Q^" << C[2] << " x I^"
         << C[3] << " x " << C[4] << ")";
  ws->setYUnitLabel(ylabel.str());
  std::stringstream xlabel;
  xlabel << "Q^" << C[5] << " x I^" << C[6] << " x Ln( Q^" << C[7] << " x I^"
         << C[8] << " x " << C[9] << ")";
  m_label->setLabel(xlabel.str());
}

/**
 * Here is the main logic to perform the transformation, to calculate the bin
 *position in degree for each spectrum.
 *
 * The first part of the method is to check if the pixel position is inside the
 *ring defined as minRadio and maxRadio.
 *
 * To do this, it deducts the pixel position. This deduction follows the
 *followin assumption:
 *
 *  - the spectrum_index == row number
 *  - the position in the 'Y' direction is given by getAxis(1)[spectrum_index]
 *  - the position in the 'X' direction is the central point of the bin
 *(dataX[column] + dataX[column+1])/2
 *
 * Having the position of the pixel, as defined above, if the distance is
 *outside the ring defined by minRadio, maxRadio,
 * it defines the bin position as -1.
 *
 * If the pixel is inside the ring, it calculates the angle of the pixel and
 *calls fromAngleToBin to define the bin
 * position.
 * @param ws: pointer to the workspace
 * @param spectrum_index: index of the spectrum
 * @param bins_pos: bin positions (for each column inside the spectrum, the
 *correspondent bin_pos)
 */
void RingProfile::getBinForPixel(const API::MatrixWorkspace_sptr ws,
                                 int spectrum_index,
                                 std::vector<int> &bins_pos) {

  if (bins_pos.size() != ws->dataY(spectrum_index).size())
    throw std::runtime_error("Invalid bin positions vector");

  API::NumericAxis *oldAxis2 = dynamic_cast<API::NumericAxis *>(ws->getAxis(1));
  // assumption y position is the ws->getAxis(1)(spectrum_index)
  if (!oldAxis2) {
    throw std::logic_error("Failed to cast workspace axis to NumericAxis");
  }

  // calculate ypos, the difference of y - centre and  the square of this
  // difference
  double ypos = (*oldAxis2)(spectrum_index);
  double diffy = ypos - centre_y;
  double diffy_quad = pow(diffy, 2.0);

  // the reference to X bins (the limits for each pixel in the horizontal
  // direction)
  auto xvec = ws->dataX(spectrum_index);

  // for each pixel inside this row
  for (size_t i = 0; i < xvec.size() - 1; i++) {

    double xpos = (xvec[i] + xvec[i + 1]) /
                  2.0; // the x position is the centre of the bins boundaries
    double diffx = xpos - centre_x;
    // calculate the distance => norm of pixel position - centre
    double distance = sqrt(pow(diffx, 2.0) + diffy_quad);

    // check if the distance is inside the ring
    if (distance < min_radius || distance > max_radius || distance == 0) {
      bins_pos[i] = -1;
      continue;
    }

    double angle = atan2(diffy, diffx);

    // call fromAngleToBin (radians)
    bins_pos[i] = fromAngleToBin(angle, false);
  }
}

/** Performs the Debye-Bueche transformation: 1/sqrt(I) v Q^2
 *  The output is set to zero for negative input Y values
 *  @param ws The workspace to be transformed
 */
void IQTransform::debyeBueche(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::Squares<double>());
  for (size_t i = 0; i < Y.size(); ++i) {
    if (Y[i] > 0.0) {
      Y[i] = 1.0 / std::sqrt(Y[i]);
      E[i] *= std::pow(Y[i], 3);
    } else {
      Y[i] = 0.0;
      E[i] = 0.0;
    }
  }

  ws->setYUnitLabel("1/sqrt(I)");
  m_label->setLabel("Q^2");
}

  /** 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->dataX(0).size(), inputWS->dataY(0).size());

      std::set<specid_t>::iterator siter;
      // Only copy the X-values from spectra with peaks specified in the table workspace.
      for (siter = m_spectraSet.begin(); siter != m_spectraSet.end(); ++siter)
      {
        specid_t iws = *siter;
        std::copy(inputWS->dataX(iws).begin(), inputWS->dataX(iws).end(), outputWS->dataX(iws).begin());
      }

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

    return outputWS;
  }