C++ MatrixWorkspace_const_sptr::getDetectorIDToWorkspaceIndexMap方法代码示例

/**
 * Computed the normalization for the input workspace. Results are stored in
 * m_normWS
 * @param otherValues
 * @param affineTrans
 */
void MDNormDirectSC::calculateNormalization(
    const std::vector<coord_t> &otherValues,
    const Kernel::Matrix<coord_t> &affineTrans) {
  constexpr double energyToK = 8.0 * M_PI * M_PI *
                               PhysicalConstants::NeutronMass *
                               PhysicalConstants::meV * 1e-20 /
                               (PhysicalConstants::h * PhysicalConstants::h);
  const auto &exptInfoZero = *(m_inputWS->getExperimentInfo(0));
  typedef Kernel::PropertyWithValue<std::vector<double>> VectorDoubleProperty;
  auto *rubwLog =
      dynamic_cast<VectorDoubleProperty *>(exptInfoZero.getLog("RUBW_MATRIX"));
  if (!rubwLog) {
    throw std::runtime_error(
        "Wokspace does not contain a log entry for the RUBW matrix."
        "Cannot continue.");
  } else {
    Kernel::DblMatrix rubwValue(
        (*rubwLog)()); // includes the 2*pi factor but not goniometer for now :)
    m_rubw = exptInfoZero.run().getGoniometerMatrix() * rubwValue;
    m_rubw.Invert();
  }
  const double protonCharge = exptInfoZero.run().getProtonCharge();

  auto instrument = exptInfoZero.getInstrument();
  std::vector<detid_t> detIDs = instrument->getDetectorIDs(true);
  // Prune out those that are part of a group and simply leave the head of the
  // group
  detIDs = removeGroupedIDs(exptInfoZero, detIDs);

  // Mapping
  const int64_t ndets = static_cast<int64_t>(detIDs.size());
  bool haveSA = false;
  API::MatrixWorkspace_const_sptr solidAngleWS =
      getProperty("SolidAngleWorkspace");
  detid2index_map solidAngDetToIdx;
  if (solidAngleWS != nullptr) {
    haveSA = true;
    solidAngDetToIdx = solidAngleWS->getDetectorIDToWorkspaceIndexMap();
  }

  auto prog = make_unique<API::Progress>(this, 0.3, 1.0, ndets);
  PARALLEL_FOR_NO_WSP_CHECK()
  for (int64_t i = 0; i < ndets; i++) {
    PARALLEL_START_INTERUPT_REGION

    const auto detID = detIDs[i];
    double theta(0.0), phi(0.0);
    bool skip(false);
    try {
      auto spectrum = getThetaPhi(detID, exptInfoZero, theta, phi);
      if (spectrum->isMonitor() || spectrum->isMasked())
        continue;
    } catch (
        std::exception &) // detector might not exist or has no been included
                          // in grouping
    {
      skip = true; // Intel compiler has a problem with continue inside a catch
                   // inside openmp...
    }
    if (skip)
      continue;

    // Intersections
    auto intersections = calculateIntersections(theta, phi);
    if (intersections.empty())
      continue;

    // Get solid angle for this contribution
    double solid = protonCharge;
    if (haveSA) {
      solid = solidAngleWS->readY(solidAngDetToIdx.find(detID)->second)[0] *
              protonCharge;
    }
    // Compute final position in HKL
    const size_t vmdDims = intersections.front().size();
    // pre-allocate for efficiency and copy non-hkl dim values into place
    std::vector<coord_t> pos(vmdDims + otherValues.size() + 1);
    std::copy(otherValues.begin(), otherValues.end(), pos.begin() + vmdDims);
    pos.push_back(1.);
    auto intersectionsBegin = intersections.begin();
    for (auto it = intersectionsBegin + 1; it != intersections.end(); ++it) {
      const auto &curIntSec = *it;
      const auto &prevIntSec = *(it - 1);
      // the full vector isn't used so compute only what is necessary
      double delta =
          (curIntSec[3] * curIntSec[3] - prevIntSec[3] * prevIntSec[3]) /
          energyToK;
      if (delta < 1e-10)
        continue; // Assume zero contribution if difference is small

      // Average between two intersections for final position
      std::transform(curIntSec.getBareArray(),
                     curIntSec.getBareArray() + vmdDims,
                     prevIntSec.getBareArray(), pos.begin(),
//.........这里部分代码省略.........

/** Get groupings from XML file
*  @param fname :: the full path name of the file to open
*  @param workspace :: a pointer to the input workspace, used to get spectra indexes from numbers
*  @param unUsedSpec :: the list of spectra indexes that have been included in a group (so far)
*  @throw FileError if there's any problem with the file or its format
*/
void GroupDetectors2::processXMLFile(std::string fname,
  API::MatrixWorkspace_const_sptr workspace, std::vector<int64_t> &unUsedSpec)
{
  // 1. Get maps for spectrum ID and detector ID
  spec2index_map specs2index;
  const SpectraAxis* axis = dynamic_cast<const SpectraAxis*>(workspace->getAxis(1));
  if (axis)
  {
    axis->getSpectraIndexMap(specs2index);
  }

  detid2index_map* detIdToWiMap = workspace->getDetectorIDToWorkspaceIndexMap(false);

  // 2. Load XML file
  DataHandling::LoadGroupXMLFile loader;
  loader.setDefaultStartingGroupID(0);
  loader.loadXMLFile(fname);
  std::map<int, std::vector<detid_t> > mGroupDetectorsMap = loader.getGroupDetectorsMap();
  std::map<int, std::vector<int> > mGroupSpectraMap = loader.getGroupSpectraMap();

  // 3. Build m_GroupSpecInds
  std::map<int, std::vector<detid_t> >::iterator dit;
  for (dit = mGroupDetectorsMap.begin(); dit != mGroupDetectorsMap.end(); ++ dit)
  {
    int groupid = dit->first;
    std::vector<size_t> tempv;
    m_GroupSpecInds.insert(std::make_pair(groupid, tempv));
  }

  // 4. Detector IDs
  for (dit = mGroupDetectorsMap.begin(); dit != mGroupDetectorsMap.end(); ++ dit)
  {
    int groupid = dit->first;
    std::vector<detid_t> detids = dit->second;

    storage_map::iterator sit;
    sit = m_GroupSpecInds.find(groupid);
    if (sit == m_GroupSpecInds.end())
      continue;

    std::vector<size_t>& wsindexes = sit->second;

    for (size_t i = 0; i < detids.size(); i++)
    {
      detid_t detid = detids[i];
      detid2index_map::iterator ind =detIdToWiMap->find(detid);
      if ( ind != detIdToWiMap->end() )
      {
        size_t wsid = ind->second;
        wsindexes.push_back(wsid);
        if ( unUsedSpec[wsid] != ( 1000 - INT_MAX ) )
        {
          unUsedSpec[wsid] = ( 1000 - INT_MAX );
        }
      }
      else
      {
        g_log.error() << "Detector with ID " << detid << " is not found in instrument " << std::endl;
      }
    } // for index
  } // for group

  // 5. Spectrum IDs
  std::map<int, std::vector<int> >::iterator pit;
  for (pit = mGroupSpectraMap.begin(); pit != mGroupSpectraMap.end(); ++pit)
  {
    int groupid = pit->first;
    std::vector<int> spectra = pit->second;

    storage_map::iterator sit;
    sit = m_GroupSpecInds.find(groupid);
    if (sit == m_GroupSpecInds.end())
      continue;

    std::vector<size_t>& wsindexes = sit->second;

    for (size_t i = 0; i < spectra.size(); i++)
    {
      int specid = spectra[i];
      spec2index_map::iterator ind = specs2index.find(specid);
      if ( ind != specs2index.end() )
      {
        size_t wsid = ind->second;
        wsindexes.push_back(wsid);
        if ( unUsedSpec[wsid] != ( 1000 - INT_MAX ) )
        {
          unUsedSpec[wsid] = ( 1000 - INT_MAX );
        }
      }
      else
      {
        g_log.error() << "Spectrum with ID " << specid<< " is not found in instrument " << std::endl;
      }
    } // for index
//.........这里部分代码省略.........

/**
 * Computed the normalization for the input workspace. Results are stored in
 * m_normWS
 * @param otherValues
 * @param affineTrans
 */
void MDNormSCD::calculateNormalization(
    const std::vector<coord_t> &otherValues,
    const Kernel::Matrix<coord_t> &affineTrans) {
  API::MatrixWorkspace_const_sptr integrFlux = getProperty("FluxWorkspace");
  integrFlux->getXMinMax(m_kiMin, m_kiMax);
  API::MatrixWorkspace_const_sptr solidAngleWS =
      getProperty("SolidAngleWorkspace");

  const auto &exptInfoZero = *(m_inputWS->getExperimentInfo(0));
  typedef Kernel::PropertyWithValue<std::vector<double>> VectorDoubleProperty;
  auto *rubwLog =
      dynamic_cast<VectorDoubleProperty *>(exptInfoZero.getLog("RUBW_MATRIX"));
  if (!rubwLog) {
    throw std::runtime_error(
        "Wokspace does not contain a log entry for the RUBW matrix."
        "Cannot continue.");
  } else {
    Kernel::DblMatrix rubwValue(
        (*rubwLog)()); // includes the 2*pi factor but not goniometer for now :)
    m_rubw = exptInfoZero.run().getGoniometerMatrix() * rubwValue;
    m_rubw.Invert();
  }
  const double protonCharge = exptInfoZero.run().getProtonCharge();

  auto instrument = exptInfoZero.getInstrument();
  std::vector<detid_t> detIDs = instrument->getDetectorIDs(true);
  // Prune out those that are part of a group and simply leave the head of the
  // group
  detIDs = removeGroupedIDs(exptInfoZero, detIDs);

  // Mappings
  const int64_t ndets = static_cast<int64_t>(detIDs.size());
  const detid2index_map fluxDetToIdx =
      integrFlux->getDetectorIDToWorkspaceIndexMap();
  const detid2index_map solidAngDetToIdx =
      solidAngleWS->getDetectorIDToWorkspaceIndexMap();

  auto prog = make_unique<API::Progress>(this, 0.3, 1.0, ndets);
  PARALLEL_FOR1(integrFlux)
  for (int64_t i = 0; i < ndets; i++) {
    PARALLEL_START_INTERUPT_REGION

    const auto detID = detIDs[i];
    double theta(0.0), phi(0.0);
    bool skip(false);
    try {
      auto spectrum = getThetaPhi(detID, exptInfoZero, theta, phi);
      if (spectrum->isMonitor() || spectrum->isMasked())
        continue;
    } catch (
        std::exception &) // detector might not exist or has no been included
                          // in grouping
    {
      skip = true; // Intel compiler has a problem with continue inside a catch
                   // inside openmp...
    }
    if (skip)
      continue;

    // Intersections
    auto intersections = calculateIntersections(theta, phi);
    if (intersections.empty())
      continue;

    // get the flux spetrum number
    size_t wsIdx = fluxDetToIdx.find(detID)->second;
    // Get solid angle for this contribution
    double solid =
        solidAngleWS->readY(solidAngDetToIdx.find(detID)->second)[0] *
        protonCharge;

    // -- calculate integrals for the intersection --
    // momentum values at intersections
    auto intersectionsBegin = intersections.begin();
    std::vector<double> xValues(intersections.size()),
        yValues(intersections.size());
    {
      // copy momenta to xValues
      auto x = xValues.begin();
      for (auto it = intersectionsBegin; it != intersections.end(); ++it, ++x) {
        *x = (*it)[3];
      }
    }
    // calculate integrals at momenta from xValues by interpolating between
    // points in spectrum sp
    // of workspace integrFlux. The result is stored in yValues
    calcIntegralsForIntersections(xValues, *integrFlux, wsIdx, yValues);

    // Compute final position in HKL
    const size_t vmdDims = intersections.front().size();
    // pre-allocate for efficiency and copy non-hkl dim values into place
    std::vector<coord_t> pos(vmdDims + otherValues.size());
    std::copy(otherValues.begin(), otherValues.end(),
              pos.begin() + vmdDims - 1);
//.........这里部分代码省略.........