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

/**
 * Outputs message to log if the detector at the given index is not a monitor in
 *both input workspaces.
 *
 * @param sampleWS :: the input sample workspace
 * @param directWS :: the input direct workspace
 * @param index    :: the index of the detector to checked
 */
void CalculateTransmission::logIfNotMonitor(API::MatrixWorkspace_sptr sampleWS,
                                            API::MatrixWorkspace_sptr directWS,
                                            size_t index) {
  const std::string message = "The detector at index " +
                              boost::lexical_cast<std::string>(index) +
                              " is not a monitor in the ";
  if (!sampleWS->getDetector(index)->isMonitor())
    g_log.information(message + "sample workspace.");
  if (!directWS->getDetector(index)->isMonitor())
    g_log.information(message + "direct workspace.");
}

void FindDetectorsPar::populate_values_from_file(
    const API::MatrixWorkspace_sptr &inputWS) {
  size_t nHist = inputWS->getNumberHistograms();

  if (this->current_ASCII_file.Type == PAR_type) {
    // in this case data in azimuthal width and polar width are in fact real
    // sizes in meters; have to transform it in into angular values
    for (size_t i = 0; i < nHist; i++) {
      azimuthalWidth[i] =
          atan2(azimuthalWidth[i], secondaryFlightpath[i]) * rad2deg;
      polarWidth[i] = atan2(polarWidth[i], secondaryFlightpath[i]) * rad2deg;
    }
    m_SizesAreLinear = false;
  } else {

    Geometry::IComponent_const_sptr sample =
        inputWS->getInstrument()->getSample();
    secondaryFlightpath.resize(nHist);
    // Loop over the spectra
    for (size_t i = 0; i < nHist; i++) {
      Geometry::IDetector_const_sptr spDet;
      try {
        spDet = inputWS->getDetector(i);
      } catch (Kernel::Exception::NotFoundError &) {
        continue;
      }
      // Check that we aren't writing a monitor...
      if (spDet->isMonitor())
        continue;
      /// this is the only value, which is not defined in phx file, so we
      /// calculate it
      secondaryFlightpath[i] = spDet->getDistance(*sample);
    }
  }
}

/** Converts X axis to theta representation
* @param progress :: Progress indicator
* @param targetUnit :: Target conversion unit
* @param inputWS :: Input Workspace
* @param nHist :: Stores the number of histograms
*/
void ConvertSpectrumAxis2::createThetaMap(API::Progress &progress,
                                          const std::string &targetUnit,
                                          API::MatrixWorkspace_sptr &inputWS,
                                          size_t nHist) {
  // Set up binding to member funtion. Avoids condition as part of loop over
  // nHistograms.
  boost::function<double(const IDetector &)> thetaFunction;
  if (targetUnit.compare("signed_theta") == 0 ||
      targetUnit.compare("SignedTheta") == 0) {
    thetaFunction =
        boost::bind(&MatrixWorkspace::detectorSignedTwoTheta, inputWS, _1);
  } else if (targetUnit == "theta" || targetUnit == "Theta") {
    thetaFunction =
        boost::bind(&MatrixWorkspace::detectorTwoTheta, inputWS, _1);
  }

  bool warningGiven = false;

  for (size_t i = 0; i < nHist; ++i) {
    try {
      IDetector_const_sptr det = inputWS->getDetector(i);
      // Invoke relevant member function.
      m_indexMap.emplace(thetaFunction(*det) * rad2deg, i);
    } catch (Exception::NotFoundError &) {
      if (!warningGiven)
        g_log.warning("The instrument definition is incomplete - spectra "
                      "dropped from output");
      warningGiven = true;
    }

    progress.report("Converting to theta...");
  }
}

double RadiusSum::getMinBinSizeForInstrument(API::MatrixWorkspace_sptr inWS) {
  // Assumption made: the detectors are placed one after the other, so the
  // minimum
  // reasonalbe size for the bin is the width of one detector.

  double width;
  size_t i = 0;
  while (true) {
    i++;

    // this should never happen because it was done in
    // getBoundariesOfInstrument,
    // but it is here to avoid risk of infiniti loop
    if (i >= inWS->getNumberHistograms())
      throw std::invalid_argument(
          "Did not find any non monitor detector position");

    auto det = inWS->getDetector(i);
    if (det->isMonitor())
      continue;

    Geometry::BoundingBox bbox;
    det->getBoundingBox(bbox);

    width = bbox.width().norm();
    break;
  }

  return width;
}

 double LoadHelper::getL2(const API::MatrixWorkspace_sptr& workspace, int detId)
 {
   // Get a pointer to the instrument contained in the workspace
   Geometry::Instrument_const_sptr instrument = workspace->getInstrument();
   // Get the distance between the source and the sample (assume in metres)
   Geometry::IComponent_const_sptr sample = instrument->getSample();
   // Get the sample-detector distance for this detector (in metres)
   double l2 = workspace->getDetector(detId)->getPos().distance(sample->getPos());
   return l2;
 }

/** Extracts a single spectrum from a Workspace2D into a new workspaces. Uses CropWorkspace to do this.
 *  @param WS ::    The workspace containing the spectrum to extract
 *  @param index :: The workspace index of the spectrum to extract
 *  @return A Workspace2D containing the extracted spectrum
 */
API::MatrixWorkspace_sptr CalculateTransmissionBeamSpreader::extractSpectrum(API::MatrixWorkspace_sptr WS, const size_t index)
{
  // Check that given spectra are monitors
  if ( !WS->getDetector(index)->isMonitor() )
  {
    g_log.information("The Incident Beam Monitor UDET provided is not marked as a monitor");
  }

  Algorithm_sptr childAlg = createSubAlgorithm("ExtractSingleSpectrum",0.0,0.4);
  childAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", WS);
  childAlg->setProperty<int>("WorkspaceIndex", static_cast<int>(index));
  childAlg->executeAsSubAlg();
  return childAlg->getProperty("OutputWorkspace");
}

/**
 * The main method to calculate the ring profile for workspaces based on
 *instruments.
 *
 * 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 sum of all the spectrum values 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::processInstrumentRingProfile(
    const API::MatrixWorkspace_sptr inputWS, std::vector<double> &output_bins) {

  for (int i = 0; i < static_cast<int>(inputWS->getNumberHistograms()); i++) {
    m_progress->report("Computing ring bins positions for detectors");
    // for the detector based, the positions will be taken from the detector
    // itself.
    try {
      Mantid::Geometry::IDetector_const_sptr det = inputWS->getDetector(i);

      // skip monitors
      if (det->isMonitor()) {
        continue;
      }

      // this part will be executed if the instrument is attached to the
      // workspace

      // get the bin position
      int bin_n = getBinForPixel(det);

      if (bin_n < 0) // -1 is the agreement for an invalid bin, or outside the
                     // ring being integrated
        continue;

      g_log.debug() << "Bin for the index " << i << " = " << bin_n
                    << " Pos = " << det->getPos() << std::endl;

      // get the reference to the spectrum
      auto spectrum_pt = inputWS->getSpectrum(i);
      const MantidVec &refY = spectrum_pt->dataY();
      // accumulate the values of this spectrum inside this bin
      for (size_t sp_ind = 0; sp_ind < inputWS->blocksize(); sp_ind++)
        output_bins[bin_n] += refY[sp_ind];

    } catch (Kernel::Exception::NotFoundError &ex) {
      g_log.information() << "It found that detector for " << i
                          << " is not valid. " << ex.what() << std::endl;
      continue;
    }
  }
}

/**
 * Validation of the inputs of the RingProfile algorithm.
 *
 * Inside this method, the Workspace is considered an instrument based
 *instrument. Each spectrum
 * has a detector associated which has a position in the 3D space.
 *
 * The main validation are:
 *  - the centre of the ring is inside the image it self.
 *  - The minimum ring is smaller than the limits of the image to allow
 *
 * @param inputWS: the input workspace
*/
void RingProfile::checkInputsForSpectraWorkspace(
    const API::MatrixWorkspace_sptr inputWS) {
  try {
    // finding the limits of the instrument
    double first_x, first_y, first_z;
    size_t i = 0;
    while (true) {
      i++;
      if (i >= inputWS->getNumberHistograms())
        throw std::invalid_argument(
            "Did not find any non monitor detector position");

      auto det = inputWS->getDetector(i);
      if (det->isMonitor())
        continue;
      first_x = det->getPos().X();
      first_y = det->getPos().Y();
      first_z = det->getPos().Z();
      break;
    }

    double last_x, last_y, last_z;
    i = inputWS->getNumberHistograms() - 1;
    while (true) {
      i--;
      if (i == 0)
        throw std::invalid_argument(
            "There is no region defined for the instrument of this workspace");

      auto det = inputWS->getDetector(i);
      if (det->isMonitor())
        continue;
      last_x = det->getPos().X();
      last_y = det->getPos().Y();
      last_z = det->getPos().Z();
      break;
    }

    double xMax, yMax, zMax;
    double xMin, yMin, zMin;
    xMax = std::max(first_x, last_x);
    yMax = std::max(first_y, last_y);
    zMax = std::max(first_z, last_z);
    xMin = std::min(first_x, last_x);
    yMin = std::min(first_y, last_y);
    zMin = std::min(first_z, last_z);

    std::stringstream limits_s;
    limits_s << "([" << xMin << ", " << xMax << "], [" << yMin << ", " << yMax
             << "], [" << zMin << ", " << zMax << "])";
    g_log.debug() << "The limits for the instrument is : " << limits_s.str()
                  << std::endl;
    int xOutside = 0, yOutside = 0, zOutside = 0;
    if (centre_x < xMin || centre_x > xMax)
      xOutside = 1;
    if (centre_y < yMin || centre_y > yMax)
      yOutside = 1;
    if (centre_z < zMin || centre_z > zMax)
      zOutside = 1;
    int summed = xOutside + yOutside + zOutside;
    // if at least 2 are outside, the centre is considered outside the box.
    if (summed >= 2) {
      std::stringstream s;
      s << "The defined centre (" << centre_x << ", " << centre_y << ", "
        << centre_z
        << ") is outside the limits of the detectors inside this instrument: "
        << limits_s.str();
      throw std::invalid_argument(s.str());
    }

    xOutside = yOutside = zOutside = 0;
    if (centre_x - min_radius > xMax || centre_x + min_radius < xMin)
      xOutside = 1;
    if (centre_y - min_radius > yMax || centre_y + min_radius < yMin)
      yOutside = 1;
    if (centre_z - min_radius > zMax || centre_z + min_radius < zMin)
      zOutside = 1;

    summed = xOutside + yOutside + zOutside;

    if (summed >= 2) {
      std::stringstream s;
      s << "The defined minRadius make the inner ring outside the limits of "
           "the detectors inside this instrument: " << limits_s.str();
      throw std::invalid_argument(s.str());
    }

//.........这里部分代码省略.........

/** Convert the workspace units using TOF as an intermediate step in the
 * conversion
 * @param fromUnit :: The unit of the input workspace
 * @param outputWS :: The output workspace
 */
void ConvertUnitsUsingDetectorTable::convertViaTOF(
    Kernel::Unit_const_sptr fromUnit, API::MatrixWorkspace_sptr outputWS) {
  using namespace Geometry;

  // Let's see if we are using a TableWorkspace to override parameters
  TableWorkspace_sptr paramWS = getProperty("DetectorParameters");

  // See if we have supplied a DetectorParameters Workspace
  // TODO: Check if paramWS is NULL and if so throw an exception

  //      const std::string l1ColumnLabel("l1");

  // Let's check all the columns exist and are readable
  try {
    auto spectraColumnTmp = paramWS->getColumn("spectra");
    auto l1ColumnTmp = paramWS->getColumn("l1");
    auto l2ColumnTmp = paramWS->getColumn("l2");
    auto twoThetaColumnTmp = paramWS->getColumn("twotheta");
    auto efixedColumnTmp = paramWS->getColumn("efixed");
    auto emodeColumnTmp = paramWS->getColumn("emode");
  } catch (...) {
    throw Exception::InstrumentDefinitionError(
        "DetectorParameter TableWorkspace is not defined correctly.");
  }

  // Now let's read them into some vectors.
  auto l1Column = paramWS->getColVector<double>("l1");
  auto l2Column = paramWS->getColVector<double>("l2");
  auto twoThetaColumn = paramWS->getColVector<double>("twotheta");
  auto efixedColumn = paramWS->getColVector<double>("efixed");
  auto emodeColumn = paramWS->getColVector<int>("emode");
  auto spectraColumn = paramWS->getColVector<int>("spectra");

  EventWorkspace_sptr eventWS =
      boost::dynamic_pointer_cast<EventWorkspace>(outputWS);
  assert(static_cast<bool>(eventWS) == m_inputEvents); // Sanity check

  Progress prog(this, 0.2, 1.0, m_numberOfSpectra);
  int64_t numberOfSpectra_i =
      static_cast<int64_t>(m_numberOfSpectra); // cast to make openmp happy

  // Get the unit object for each workspace
  Kernel::Unit_const_sptr outputUnit = outputWS->getAxis(0)->unit();

  std::vector<double> emptyVec;
  int failedDetectorCount = 0;

  // ConstColumnVector<int> spectraNumber = paramWS->getVector("spectra");

  // TODO: Check why this parallel stuff breaks
  // Loop over the histograms (detector spectra)
  // PARALLEL_FOR1(outputWS)
  for (int64_t i = 0; i < numberOfSpectra_i; ++i) {

    // Lets find what row this spectrum ID appears in our detector table.

    // PARALLEL_START_INTERUPT_REGION

    std::size_t wsid = i;

    try {

      double deg2rad = M_PI / 180.;

      auto det = outputWS->getDetector(i);
      int specid = det->getID();

      // int spectraNumber = static_cast<int>(spectraColumn->toDouble(i));
      // wsid = outputWS->getIndexFromSpectrumNumber(spectraNumber);
      g_log.debug() << "###### Spectra #" << specid
                    << " ==> Workspace ID:" << wsid << std::endl;

      // Now we need to find the row that contains this spectrum
      std::vector<int>::iterator specIter;

      specIter = std::find(spectraColumn.begin(), spectraColumn.end(), specid);
      if (specIter != spectraColumn.end()) {
        size_t detectorRow = std::distance(spectraColumn.begin(), specIter);
        double l1 = l1Column[detectorRow];
        double l2 = l2Column[detectorRow];
        double twoTheta = twoThetaColumn[detectorRow] * deg2rad;
        double efixed = efixedColumn[detectorRow];
        int emode = emodeColumn[detectorRow];

        g_log.debug() << "specId from detector table = "
                      << spectraColumn[detectorRow] << std::endl;

        // l1 = l1Column->toDouble(detectorRow);
        // l2 = l2Column->toDouble(detectorRow);
        // twoTheta = deg2rad * twoThetaColumn->toDouble(detectorRow);
        // efixed = efixedColumn->toDouble(detectorRow);
        // emode = static_cast<int>(emodeColumn->toDouble(detectorRow));

        g_log.debug() << "###### Spectra #" << specid
                      << " ==> Det Table Row:" << detectorRow << std::endl;
//.........这里部分代码省略.........

/** Assuming that the workspace has an instrument associated with it from which
 *the pixel positions has to be taken,
 *  this function extracts the position of the first and last valid pixel
 *(detector) and return a list of values
 *  giving the boundaries of the instrument.
 *
 * @param inWS Input Workspace
 * @return a list of values that defines the limits of the image in this order:
 *Xmin, Xmax, Ymin, Ymax, Zmin, Zmax
 */
std::vector<double>
RadiusSum::getBoundariesOfInstrument(API::MatrixWorkspace_sptr inWS) {

  // This function is implemented based in the following assumption:
  //   - The workspace is composed by spectrum with associated spectrum No which
  //   is associated to one detector or monitor
  //   - The first spectrum No (non monitor) is associated with one detector
  //   while the last spectrum No (non monitor)
  //     is associated with one detector.
  //   - They are in complete oposite direction.
  //
  //   Consider the following 'image' (where the ID is the number and the
  //   position is where it is displayed)
  //
  //    1  2  3
  //    4  5  6
  //    7  8  9
  //   10 11 12
  //
  //    In this image, the assumption is true, because, we can derive the
  //    boundaries of the image looking just to the
  //    ids 1 and 12.
  //
  //    But the following image:
  //
  //   1  2  3       6  5  4
  //   6  5  4       1  2  3
  //   7  8  9      12 11 10
  //  12 11 12       7  8  9
  //
  //   Although valid 'IDF' instrument, fail the assumption, and will return
  //   wrong values.
  //   Bear in mind these words if you face problems with the return of the
  //   boundaries of one instrument
  //

  double first_x, first_y, first_z;
  size_t i = 0;
  while (true) {
    i++;
    if (i >= inWS->getNumberHistograms())
      throw std::invalid_argument("Did not find any non monitor detector. "
                                  "Failed to identify the boundaries of this "
                                  "instrument.");

    auto det = inWS->getDetector(i);
    if (det->isMonitor())
      continue;
    // get the position of the first valid (non-monitor) detector.
    first_x = det->getPos().X();
    first_y = det->getPos().Y();
    first_z = det->getPos().Z();
    break;
  }

  double last_x, last_y, last_z;
  i = inWS->getNumberHistograms() - 1;
  while (true) {
    i--;
    if (i == 0)
      throw std::invalid_argument("There is no region defined for the "
                                  "instrument of this workspace. Failed to "
                                  "identify the boundaries of this instrument");

    auto det = inWS->getDetector(i);
    if (det->isMonitor())
      continue;
    // get the last valid detector position
    last_x = det->getPos().X();
    last_y = det->getPos().Y();
    last_z = det->getPos().Z();
    break;
  }

  // order the values
  double xMax, yMax, zMax;
  double xMin, yMin, zMin;
  xMax = std::max(first_x, last_x);
  yMax = std::max(first_y, last_y);
  zMax = std::max(first_z, last_z);
  xMin = std::min(first_x, last_x);
  yMin = std::min(first_y, last_y);
  zMin = std::min(first_z, last_z);

  std::vector<double> output(6); // output  = {xMin, xMax, yMin, yMax, zMin,
                                 // zMax }; not supported in all compilers
  output[0] = xMin;
  output[1] = xMax;
  output[2] = yMin;
  output[3] = yMax;
//.........这里部分代码省略.........