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

/** Gets the distances between the source and detectors whose IDs you pass to it
*  @param WS :: the input workspace
*  @param mon0Spec :: Spectrum number of the output from the first monitor
*  @param mon1Spec :: Spectrum number of the output from the second monitor
*  @param monitor0Dist :: the calculated distance to the detector whose ID was
* passed to this function first
*  @param monitor1Dist :: calculated distance to the detector whose ID was
* passed to this function second
*  @throw NotFoundError if no detector is found for the detector ID given
*/
void GetEi::getGeometry(API::MatrixWorkspace_const_sptr WS, specid_t mon0Spec,
                        specid_t mon1Spec, double &monitor0Dist,
                        double &monitor1Dist) const {
    const IComponent_const_sptr source = WS->getInstrument()->getSource();

    // retrieve a pointer to the first detector and get its distance
    size_t monWI = 0;
    try {
        monWI = WS->getIndexFromSpectrumNumber(mon0Spec);
    } catch (std::runtime_error &) {
        g_log.error()
                << "Could not find the workspace index for the monitor at spectrum "
                << mon0Spec << "\n";
        g_log.error() << "Error retrieving data for the first monitor" << std::endl;
        throw std::bad_cast();
    }
    const std::set<detid_t> &dets = WS->getSpectrum(monWI)->getDetectorIDs();

    if (dets.size() != 1) {
        g_log.error() << "The detector for spectrum number " << mon0Spec
                      << " was either not found or is a group, grouped monitors "
                      "are not supported by this algorithm\n";
        g_log.error() << "Error retrieving data for the first monitor" << std::endl;
        throw std::bad_cast();
    }
    IDetector_const_sptr det = WS->getInstrument()->getDetector(*dets.begin());
    monitor0Dist = det->getDistance(*(source.get()));

    // repeat for the second detector
    try {
        monWI = WS->getIndexFromSpectrumNumber(mon0Spec);
    } catch (std::runtime_error &) {
        g_log.error()
                << "Could not find the workspace index for the monitor at spectrum "
                << mon0Spec << "\n";
        g_log.error() << "Error retrieving data for the second monitor\n";
        throw std::bad_cast();
    }
    const std::set<detid_t> &dets2 = WS->getSpectrum(monWI)->getDetectorIDs();
    if (dets2.size() != 1) {
        g_log.error() << "The detector for spectrum number " << mon1Spec
                      << " was either not found or is a group, grouped monitors "
                      "are not supported by this algorithm\n";
        g_log.error() << "Error retrieving data for the second monitor\n";
        throw std::bad_cast();
    }
    det = WS->getInstrument()->getDetector(*dets2.begin());
    monitor1Dist = det->getDistance(*(source.get()));
}

/** Checks that the two input workspaces have non-overlapping spectra numbers and contributing detectors
 *  @param ws1 :: The first input workspace
 *  @param ws2 :: The second input workspace
 *  @param checkSpectra :: set to true to check for overlapping spectra numbers (non-sensical for event workspaces)
 *  @throw std::invalid_argument If there is some overlap
 */
void ConjoinWorkspaces::checkForOverlap(API::MatrixWorkspace_const_sptr ws1, API::MatrixWorkspace_const_sptr ws2, bool checkSpectra) const
{
  // make sure we should bother checking
  if (!this->getProperty("CheckOverlapping"))
    return;
  // Loop through the first workspace adding all the spectrum numbers & UDETS to a set
  std::set<specid_t> spectra;
  std::set<detid_t> detectors;
  const size_t& nhist1 = ws1->getNumberHistograms();
  for (size_t i = 0; i < nhist1; ++i)
  {
    const ISpectrum * spec = ws1->getSpectrum(i);
    const specid_t spectrum = spec->getSpectrumNo();
    spectra.insert(spectrum);
    const std::set<detid_t> & dets = spec->getDetectorIDs();
    std::set<detid_t>::const_iterator it;
    for (it = dets.begin(); it != dets.end(); ++it)
    {
      detectors.insert(*it);
    }
  }

  // Now go throught the spectrum numbers & UDETS in the 2nd workspace, making sure that there's no overlap
  const size_t& nhist2 = ws2->getNumberHistograms();
  for (size_t j = 0; j < nhist2; ++j)
  {
    const ISpectrum * spec = ws2->getSpectrum(j);
    const specid_t spectrum = spec->getSpectrumNo();
    if (checkSpectra)
    {
      if ( spectrum > 0 && spectra.find(spectrum) != spectra.end() )
      {
        g_log.error("The input workspaces have overlapping spectrum numbers");
        throw std::invalid_argument("The input workspaces have overlapping spectrum numbers");
      }
    }
    const std::set<detid_t> & dets = spec->getDetectorIDs();
    std::set<detid_t>::const_iterator it;
    for (it = dets.begin(); it != dets.end(); ++it)
    {
      if ( detectors.find(*it) != detectors.end() )
      {
        g_log.error("The input workspaces have common detectors");
        throw std::invalid_argument("The input workspaces have common detectors");
      }
    }
  }
}

/** Checks that the two input workspaces have non-overlapping spectra numbers
 * and contributing detectors
 *  @param ws1 :: The first input workspace
 *  @param ws2 :: The second input workspace
 *  @param checkSpectra :: set to true to check for overlapping spectra numbers
 * (non-sensical for event workspaces)
 *  @throw std::invalid_argument If there is some overlap
 */
void ConjoinWorkspaces::checkForOverlap(API::MatrixWorkspace_const_sptr ws1,
                                        API::MatrixWorkspace_const_sptr ws2,
                                        bool checkSpectra) const {
  // Loop through the first workspace adding all the spectrum numbers & UDETS to
  // a set
  std::set<specnum_t> spectra;
  std::set<detid_t> detectors;
  const size_t &nhist1 = ws1->getNumberHistograms();
  for (size_t i = 0; i < nhist1; ++i) {
    const ISpectrum *spec = ws1->getSpectrum(i);
    const specnum_t spectrum = spec->getSpectrumNo();
    spectra.insert(spectrum);
    const auto &dets = spec->getDetectorIDs();
    for (auto const &det : dets) {
      detectors.insert(det);
    }
  }

  // Now go throught the spectrum numbers & UDETS in the 2nd workspace, making
  // sure that there's no overlap
  const size_t &nhist2 = ws2->getNumberHistograms();
  for (size_t j = 0; j < nhist2; ++j) {
    const ISpectrum *spec = ws2->getSpectrum(j);
    const specnum_t spectrum = spec->getSpectrumNo();
    if (checkSpectra) {
      if (spectrum > 0 && spectra.find(spectrum) != spectra.end()) {
        g_log.error()
            << "The input workspaces have overlapping spectrum numbers "
            << spectrum << "\n";
        throw std::invalid_argument(
            "The input workspaces have overlapping spectrum numbers");
      }
    }
    const auto &dets = spec->getDetectorIDs();
    for (const auto &det : dets) {
      if (detectors.find(det) != detectors.end()) {
        g_log.error() << "The input workspaces have common detectors: " << (det)
                      << "\n";
        throw std::invalid_argument(
            "The input workspaces have common detectors");
      }
    }
  }
}

/**
*  Only to be used if the KeepUnGrouped property is true, moves the spectra that were not selected
*  to be in a group to the end of the output spectrum
*  @param unGroupedSet :: list of WORKSPACE indexes that were included in a group
*  @param inputWS :: user selected input workspace for the algorithm
*  @param outputWS :: user selected output workspace for the algorithm
*  @param outIndex :: the next spectra index available after the grouped spectra
*/
void GroupDetectors2::moveOthers(const std::set<int64_t> &unGroupedSet, API::MatrixWorkspace_const_sptr inputWS, API::MatrixWorkspace_sptr outputWS,
         size_t outIndex)
{
  g_log.debug() << "Starting to copy the ungrouped spectra" << std::endl;
  double prog4Copy = (1. - 1.*static_cast<double>(m_FracCompl))/static_cast<double>(unGroupedSet.size());

  std::set<int64_t>::const_iterator copyFrIt = unGroupedSet.begin();
  // go thorugh all the spectra in the input workspace
  for ( ; copyFrIt != unGroupedSet.end(); ++copyFrIt )
  {
    if( *copyFrIt == USED ) continue; //Marked as not to be used
    size_t sourceIndex = static_cast<size_t>(*copyFrIt);

    // The input spectrum we'll copy
    const ISpectrum * inputSpec = inputWS->getSpectrum(sourceIndex);

    // Destination of the copying
    ISpectrum * outputSpec = outputWS->getSpectrum(outIndex);

    // Copy the data
    outputSpec->dataX() = inputSpec->dataX();
    outputSpec->dataY() = inputSpec->dataY();
    outputSpec->dataE() = inputSpec->dataE();

    // Spectrum numbers etc.
    outputSpec->setSpectrumNo(inputSpec->getSpectrumNo());
    outputSpec->clearDetectorIDs();
    outputSpec->addDetectorIDs( inputSpec->getDetectorIDs() );

    // go to the next free index in the output workspace
    outIndex ++;
    // make regular progress reports and check for cancelling the algorithm
    if ( outIndex % INTERVAL == 0 )
    {
      m_FracCompl += INTERVAL*prog4Copy;
      if ( m_FracCompl > 1.0 )
      {
        m_FracCompl = 1.0;
      }
      progress(m_FracCompl);
      interruption_point();
    }
  }
  // Refresh the spectraDetectorMap
  outputWS->generateSpectraMap();

  g_log.debug() << name() << " copied " << unGroupedSet.size()-1 << " ungrouped spectra\n";
}

/**
 * If the InstrumentParameter property is set then it attempts to retrieve the parameter
 * from the component, else it returns the value of the Factor property
 * @param inputWS A pointer to the input workspace
 * @param index The current index to inspect
 * @return Value for the scale factor
 */
double ScaleX::getScaleFactor(const API::MatrixWorkspace_const_sptr & inputWS, const size_t index)
{
  if(m_parname.empty()) return m_algFactor;

  // Try and get factor from component. If we see a DetectorGroup use this will use the first component
  Geometry::IDetector_const_sptr det;
  auto inst = inputWS->getInstrument();

  auto *spec = inputWS->getSpectrum(index);
  const auto & ids = spec->getDetectorIDs();
  const size_t ndets(ids.size());
  if(ndets > 0)
  {
    try
    {
      det = inst->getDetector(*ids.begin());
    }
    catch(Exception::NotFoundError&)
    {
      return 0.0;
    }
  }
  else return 0.0;

  const auto & pmap = inputWS->constInstrumentParameters();
  auto par = pmap.getRecursive(det->getComponentID(), m_parname);
  if(par)
  {
    if(!m_combine) return par->value<double>();
    else return m_binOp(m_algFactor,par->value<double>());
  }
  else
  {
    std::ostringstream os;
    os << "Spectrum at index '" << index << "' has no parameter named '" << m_parname << "'\n";
    throw std::runtime_error(os.str());
  }
}

/**
*  Move the user selected spectra in the input workspace into groups in the output workspace
*  @param inputWS :: user selected input workspace for the algorithm
*  @param outputWS :: user selected output workspace for the algorithm
*  @param prog4Copy :: the amount of algorithm progress to attribute to moving a single spectra
*  @return number of new grouped spectra
*/
size_t GroupDetectors2::formGroups( API::MatrixWorkspace_const_sptr inputWS, API::MatrixWorkspace_sptr outputWS, 
            const double prog4Copy)
{
  // get "Behaviour" string
  const std::string behaviour = getProperty("Behaviour");
  int bhv = 0;
  if ( behaviour == "Average" ) bhv = 1;

  API::MatrixWorkspace_sptr beh = API::WorkspaceFactory::Instance().create(
    "Workspace2D", static_cast<int>(m_GroupSpecInds.size()), 1, 1);

  g_log.debug() << name() << ": Preparing to group spectra into " << m_GroupSpecInds.size() << " groups\n";

  // where we are copying spectra to, we start copying to the start of the output workspace
  size_t outIndex = 0;
  // Only used for averaging behaviour. We may have a 1:1 map where a Divide would be waste as it would be just dividing by 1
  bool requireDivide(false);
  for ( storage_map::const_iterator it = m_GroupSpecInds.begin(); it != m_GroupSpecInds.end() ; ++it )
  {
    // This is the grouped spectrum
    ISpectrum * outSpec = outputWS->getSpectrum(outIndex);

    // The spectrum number of the group is the key
    outSpec->setSpectrumNo(it->first);
    // Start fresh with no detector IDs
    outSpec->clearDetectorIDs();

    // Copy over X data from first spectrum, the bin boundaries for all spectra are assumed to be the same here
    outSpec->dataX() = inputWS->readX(0);

    // the Y values and errors from spectra being grouped are combined in the output spectrum
    // Keep track of number of detectors required for masking
    size_t nonMaskedSpectra(0);
    beh->dataX(outIndex)[0] = 0.0;
    beh->dataE(outIndex)[0] = 0.0;
    for( std::vector<size_t>::const_iterator wsIter = it->second.begin(); wsIter != it->second.end(); ++wsIter)
    {
      const size_t originalWI = *wsIter;

      // detectors to add to firstSpecNum
      const ISpectrum * fromSpectrum = inputWS->getSpectrum(originalWI);

      // Add up all the Y spectra and store the result in the first one
      // Need to keep the next 3 lines inside loop for now until ManagedWorkspace mru-list works properly
      MantidVec &firstY = outSpec->dataY();
      MantidVec::iterator fYit;
      MantidVec::iterator fEit = outSpec->dataE().begin();
      MantidVec::const_iterator Yit = fromSpectrum->dataY().begin();
      MantidVec::const_iterator Eit = fromSpectrum->dataE().begin();
      for (fYit = firstY.begin(); fYit != firstY.end(); ++fYit, ++fEit, ++Yit, ++Eit)
      {
        *fYit += *Yit;
        // Assume 'normal' (i.e. Gaussian) combination of errors
        *fEit = std::sqrt( (*fEit)*(*fEit) + (*Eit)*(*Eit) );
      }

      // detectors to add to the output spectrum
      outSpec->addDetectorIDs(fromSpectrum->getDetectorIDs() );
      try
      {
        Geometry::IDetector_const_sptr det = inputWS->getDetector(originalWI);
        if( !det->isMasked() ) ++nonMaskedSpectra;
      }
      catch(Exception::NotFoundError&)
      {
        // If a detector cannot be found, it cannot be masked
        ++nonMaskedSpectra;
      }
    }
    if( nonMaskedSpectra == 0 ) ++nonMaskedSpectra; // Avoid possible divide by zero
    if(!requireDivide) requireDivide = (nonMaskedSpectra > 1);
    beh->dataY(outIndex)[0] = static_cast<double>(nonMaskedSpectra);

    // make regular progress reports and check for cancelling the algorithm
    if ( outIndex % INTERVAL == 0 )
    {
      m_FracCompl += INTERVAL*prog4Copy;
      if ( m_FracCompl > 1.0 )
        m_FracCompl = 1.0;
      progress(m_FracCompl);
      interruption_point();
    }
    outIndex ++;
  }
  
  // Refresh the spectraDetectorMap
  outputWS->generateSpectraMap();

  if ( bhv == 1 && requireDivide )
  {
    g_log.debug() << "Running Divide algorithm to perform averaging.\n";
    Mantid::API::IAlgorithm_sptr divide = createChildAlgorithm("Divide");
    divide->initialize();
//.........这里部分代码省略.........

  /**
   * Function that retrieves the two-theta and azimuthal angles from a given
   * detector. It then looks up the nearest neighbours. Using those detectors,
   * it calculates the two-theta and azimuthal angle widths.
   * @param workspace : the workspace containing the needed detector information
   */
  void SofQW3::getValuesAndWidths(API::MatrixWorkspace_const_sptr workspace)
  {
    // Trigger a build of the nearst neighbors outside the OpenMP loop
    const int numNeighbours = 4;
    const size_t nHistos = workspace->getNumberHistograms();
    g_log.debug() << "Number of Histograms: " << nHistos << std::endl;

    this->m_theta = std::vector<double>(nHistos);
    this->m_thetaWidths = std::vector<double>(nHistos);
    this->m_phi = std::vector<double>(nHistos);
    this->m_phiWidths = std::vector<double>(nHistos);

    for (size_t i = 0; i < nHistos; ++i)
    {
      m_progress->report("Calculating detector angular widths");
      DetConstPtr detector = workspace->getDetector(i);
      g_log.debug() << "Current histogram: " << i << std::endl;
      specid_t inSpec = workspace->getSpectrum(i)->getSpectrumNo();
      SpectraDistanceMap neighbours = workspace->getNeighboursExact(inSpec,
                                                                    numNeighbours,
                                                                    true);

      g_log.debug() << "Current ID: " << inSpec << std::endl;
      // Convert from spectrum numbers to workspace indices
      double thetaWidth = -DBL_MAX;
      double phiWidth = -DBL_MAX;

      // Find theta and phi widths
      double theta = workspace->detectorTwoTheta(detector);
      double phi = detector->getPhi();

      specid_t deltaPlus1 = inSpec + 1;
      specid_t deltaMinus1 = inSpec - 1;
      specid_t deltaPlusT = inSpec + this->m_detNeighbourOffset;
      specid_t deltaMinusT = inSpec - this->m_detNeighbourOffset;

      for (SpectraDistanceMap::iterator it = neighbours.begin();
           it != neighbours.end(); ++it)
      {
        specid_t spec = it->first;
        g_log.debug() << "Neighbor ID: " << spec << std::endl;
        if (spec == deltaPlus1 || spec == deltaMinus1 ||
            spec == deltaPlusT || spec == deltaMinusT)
        {
          DetConstPtr detector_n = workspace->getDetector(spec - 1);
          double theta_n = workspace->detectorTwoTheta(detector_n);
          double phi_n = detector_n->getPhi();

          double dTheta = std::fabs(theta - theta_n);
          double dPhi = std::fabs(phi - phi_n);
          if (dTheta > thetaWidth)
          {
            thetaWidth = dTheta;
            //g_log.debug() << "Current ThetaWidth: " << thetaWidth << std::endl;
          }
          if (dPhi > phiWidth)
          {
            phiWidth = dPhi;
            //g_log.debug() << "Current PhiWidth: " << phiWidth << std::endl;
          }
        }
      }
      this->m_theta[i] = theta;
      this->m_phi[i] = phi;
      this->m_thetaWidths[i] = thetaWidth;
      this->m_phiWidths[i] = phiWidth;
    }
  }

/**
 * Apply the detector test criterion
 * @param counts1 :: A workspace containing the integrated counts of the first
 * white beam run
 * @param counts2 :: A workspace containing the integrated counts of the first
 * white beam run
 * @param average :: The computed median
 * @param variation :: The allowed variation in terms of number of medians, i.e
 * those spectra where
 * the ratio of the counts outside this range will fail the tests and will be
 * masked on counts1
 * @return number of detectors for which tests failed
 */
int DetectorEfficiencyVariation::doDetectorTests(
    API::MatrixWorkspace_const_sptr counts1,
    API::MatrixWorkspace_const_sptr counts2, const double average,
    double variation) {
  // DIAG in libISIS did this.  A variation of less than 1 doesn't make sense in
  // this algorithm
  if (variation < 1) {
    variation = 1.0 / variation;
  }
  // criterion for if the the first spectrum is larger than expected
  double largest = average * variation;
  // criterion for if the the first spectrum is lower than expected
  double lowest = average / variation;

  const int numSpec = static_cast<int>(counts1->getNumberHistograms());
  const int progStep = static_cast<int>(std::ceil(numSpec / 30.0));

  // Create a workspace for the output
  MaskWorkspace_sptr maskWS = this->generateEmptyMask(counts1);

  bool checkForMask = false;
  Geometry::Instrument_const_sptr instrument = counts1->getInstrument();
  if (instrument != nullptr) {
    checkForMask = ((instrument->getSource() != nullptr) &&
                    (instrument->getSample() != nullptr));
  }

  const double deadValue(1.0);
  int numFailed(0);
  PARALLEL_FOR3(counts1, counts2, maskWS)
  for (int i = 0; i < numSpec; ++i) {
    PARALLEL_START_INTERUPT_REGION
    // move progress bar
    if (i % progStep == 0) {
      advanceProgress(progStep * static_cast<double>(RTMarkDetects) / numSpec);
      progress(m_fracDone);
      interruption_point();
    }

    if (checkForMask) {
      const std::set<detid_t> &detids =
          counts1->getSpectrum(i)->getDetectorIDs();
      if (instrument->isMonitor(detids))
        continue;
      if (instrument->isDetectorMasked(detids)) {
        // Ensure it is masked on the output
        maskWS->dataY(i)[0] = deadValue;
        continue;
      }
    }

    const double signal1 = counts1->readY(i)[0];
    const double signal2 = counts2->readY(i)[0];

    // Mask out NaN and infinite
    if (boost::math::isinf(signal1) || boost::math::isnan(signal1) ||
        boost::math::isinf(signal2) || boost::math::isnan(signal2)) {
      maskWS->dataY(i)[0] = deadValue;
      PARALLEL_ATOMIC
      ++numFailed;
      continue;
    }

    // Check the ratio is within the given range
    const double ratio = signal1 / signal2;
    if (ratio < lowest || ratio > largest) {
      maskWS->dataY(i)[0] = deadValue;
      PARALLEL_ATOMIC
      ++numFailed;
    }

    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION

  // Register the results with the ADS
  setProperty("OutputWorkspace", maskWS);

  return numFailed;
}