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

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

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

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

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

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

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

/** Unwraps the Y & E vectors of a spectrum according to the ranges found in
 * unwrapX.
 *  @param tempWS ::      A pointer to the temporary workspace in which the
 * results are being stored
 *  @param spectrum ::    The workspace index
 *  @param rangeBounds :: The upper and lower ranges for the unwrapping
 */
void UnwrapMonitor::unwrapYandE(const API::MatrixWorkspace_sptr &tempWS,
                                const int &spectrum,
                                const std::vector<int> &rangeBounds) {
  // Copy over the relevant ranges of Y & E data
  MantidVec &Y = tempWS->dataY(spectrum);
  MantidVec &E = tempWS->dataE(spectrum);
  // Get references to the input data
  const MantidVec &YIn = m_inputWS->dataY(spectrum);
  const MantidVec &EIn = m_inputWS->dataE(spectrum);
  if (rangeBounds[2] != -1) {
    // Copy in the upper range
    Y.assign(YIn.begin() + rangeBounds[2], YIn.end());
    E.assign(EIn.begin() + rangeBounds[2], EIn.end());
    // Propagate masking, if necessary
    if (m_inputWS->hasMaskedBins(spectrum)) {
      const MatrixWorkspace::MaskList &inputMasks =
          m_inputWS->maskedBins(spectrum);
      MatrixWorkspace::MaskList::const_iterator it;
      for (it = inputMasks.begin(); it != inputMasks.end(); ++it) {
        if (static_cast<int>((*it).first) >= rangeBounds[2])
          tempWS->flagMasked(spectrum, (*it).first - rangeBounds[2],
                             (*it).second);
      }
    }
  } else {
    // Y & E are references to existing vector. Assign above clears them, so
    // need to explicitly here
    Y.clear();
    E.clear();
  }
  if (rangeBounds[0] != -1 && rangeBounds[1] > 0) {
    // Now append the lower range
    MantidVec::const_iterator YStart = YIn.begin();
    MantidVec::const_iterator EStart = EIn.begin();
    Y.insert(Y.end(), YStart + rangeBounds[0], YStart + rangeBounds[1]);
    E.insert(E.end(), EStart + rangeBounds[0], EStart + rangeBounds[1]);
    // Propagate masking, if necessary
    if (m_inputWS->hasMaskedBins(spectrum)) {
      const MatrixWorkspace::MaskList &inputMasks =
          m_inputWS->maskedBins(spectrum);
      MatrixWorkspace::MaskList::const_iterator it;
      for (it = inputMasks.begin(); it != inputMasks.end(); ++it) {
        const int maskIndex = static_cast<int>((*it).first);
        if (maskIndex >= rangeBounds[0] && maskIndex < rangeBounds[1])
          tempWS->flagMasked(spectrum, maskIndex - rangeBounds[0],
                             (*it).second);
      }
    }
  }
}

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

/** 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
  MantidVec &Y = workspace->dataY(index);
  const size_t nbins = workspace->blocksize();
  int retval;
  for (size_t i = 0; i < nbins; ++i) {
    retval = fscanf(speFile, "%10le", &Y[i]);
    // g_log.error() << Y[i] << std::endl;
    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
  MantidVec &E = workspace->dataE(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);

  return;
}

/** 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) {
    MantidVec &Y = outputWS->dataY(i);
    MantidVec &E = outputWS->dataE(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]));
  }
}

/**
 * Sum counts from the input workspace in lambda along lines of constant Q by
 * projecting to "virtual lambda" at a reference angle.
 *
 * @param detectorWS [in] :: the input workspace in wavelength
 * @param indices [in] :: an index set defining the foreground histograms
 * @return :: the single histogram output workspace in wavelength
 */
API::MatrixWorkspace_sptr
ReflectometrySumInQ::sumInQ(const API::MatrixWorkspace &detectorWS,
                            const Indexing::SpectrumIndexSet &indices) {

  const auto spectrumInfo = detectorWS.spectrumInfo();
  const auto refAngles = referenceAngles(spectrumInfo);
  // Construct the output workspace in virtual lambda
  API::MatrixWorkspace_sptr IvsLam =
      constructIvsLamWS(detectorWS, indices, refAngles);
  auto &outputE = IvsLam->dataE(0);
  // Loop through each spectrum in the detector group
  for (auto spIdx : indices) {
    if (spectrumInfo.isMasked(spIdx) || spectrumInfo.isMonitor(spIdx)) {
      continue;
    }
    // Get the size of this detector in twoTheta
    const auto twoThetaRange = twoThetaWidth(spIdx, spectrumInfo);
    // Check X length is Y length + 1
    const auto inputBinEdges = detectorWS.binEdges(spIdx);
    const auto inputCounts = detectorWS.counts(spIdx);
    const auto inputStdDevs = detectorWS.countStandardDeviations(spIdx);
    // Create a vector for the projected errors for this spectrum.
    // (Output Y values can simply be accumulated directly into the output
    // workspace, but for error values we need to create a separate error
    // vector for the projected errors from each input spectrum and then
    // do an overall sum in quadrature.)
    std::vector<double> projectedE(outputE.size(), 0.0);
    // Process each value in the spectrum
    const int ySize = static_cast<int>(inputCounts.size());
    for (int inputIdx = 0; inputIdx < ySize; ++inputIdx) {
      // Do the summation in Q
      processValue(inputIdx, twoThetaRange, refAngles, inputBinEdges,
                   inputCounts, inputStdDevs, *IvsLam, projectedE);
    }
    // Sum errors in quadrature
    const int eSize = static_cast<int>(outputE.size());
    for (int outIdx = 0; outIdx < eSize; ++outIdx) {
      outputE[outIdx] += projectedE[outIdx] * projectedE[outIdx];
    }
  }

  // Take the square root of all the accumulated squared errors for this
  // detector group. Assumes Gaussian errors
  double (*rs)(double) = std::sqrt;
  std::transform(outputE.begin(), outputE.end(), outputE.begin(), rs);

  return IvsLam;
}