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

    /**
     * A map detector ID and Q ranges
     * This method looks unnecessary as it could be calculated on the fly but
     * the parallelization means that lazy instantation slows it down due to the
     * necessary CRITICAL sections required to update the cache. The Q range
     * values are required very frequently so the total time is more than
     * offset by this precaching step
     */
    void SofQW2::initThetaCache(API::MatrixWorkspace_const_sptr workspace)
    {
      const size_t nhist = workspace->getNumberHistograms();
      m_thetaPts = std::vector<double>(nhist);
      size_t ndets(0);
      double minTheta(DBL_MAX), maxTheta(-DBL_MAX);

      for(int64_t i = 0 ; i < (int64_t)nhist; ++i) //signed for OpenMP
      {

        m_progress->report("Calculating detector angles");
        IDetector_const_sptr det;
        try
        {
          det = workspace->getDetector(i);
          // Check to see if there is an EFixed, if not skip it
          try
          {
            m_EmodeProperties.getEFixed(det);
          }
          catch(std::runtime_error&)
          {
            det.reset();
          }
        }
        catch(Kernel::Exception::NotFoundError&)
        {
          // Catch if no detector. Next line tests whether this happened - test placed
          // outside here because Mac Intel compiler doesn't like 'continue' in a catch
          // in an openmp block.
        }
        // If no detector found, skip onto the next spectrum
        if( !det || det->isMonitor() )
        {
          m_thetaPts[i] = -1.0; // Indicates a detector to skip
        }
        else
        {
          ++ndets;
          const double theta = workspace->detectorTwoTheta(det);
          m_thetaPts[i] = theta;
          if( theta < minTheta )
          {
            minTheta = theta;
          }
          else if( theta > maxTheta )
          {
            maxTheta = theta;
          }
        }
      }

      m_thetaWidth = (maxTheta - minTheta)/static_cast<double>(ndets);
      g_log.information() << "Calculated detector width in theta=" << (m_thetaWidth*180.0/M_PI) << " degrees.\n";
    }

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

//.........这里部分代码省略.........
  auto &L2 = targWS->getColVector<double>("L2");
  auto &TwoTheta = targWS->getColVector<double>("TwoTheta");
  auto &Azimuthal = targWS->getColVector<double>("Azimuthal");
  auto &detDir = targWS->getColVector<Kernel::V3D>("DetDirections");

  // Efixed; do we need one and does one exist?
  double Efi = targWS->getLogs()->getPropertyValueAsType<double>("Ei");
  float *pEfixedArray(nullptr);
  const Geometry::ParameterMap &pmap = inputWS->constInstrumentParameters();
  if (m_getEFixed)
    pEfixedArray = targWS->getColDataArray<float>("eFixed");

  // check if one needs to generate masked detectors column.
  int *pMasksArray(nullptr);
  if (m_getIsMasked)
    pMasksArray = targWS->getColDataArray<int>("detMask");

  //// progress message appearance
  size_t div = 100;
  size_t nHist = targWS->rowCount();
  Mantid::API::Progress theProgress(this, 0, 1, nHist);
  //// Loop over the spectra
  uint32_t liveDetectorsCount(0);
  for (size_t i = 0; i < nHist; i++) {
    sp2detMap[i] = std::numeric_limits<uint64_t>::quiet_NaN();
    detId[i] = std::numeric_limits<int32_t>::quiet_NaN();
    detIDMap[i] = std::numeric_limits<uint64_t>::quiet_NaN();
    L2[i] = std::numeric_limits<double>::quiet_NaN();
    TwoTheta[i] = std::numeric_limits<double>::quiet_NaN();
    Azimuthal[i] = std::numeric_limits<double>::quiet_NaN();
    //     detMask[i]  = true;

    // get detector or detector group which corresponds to the spectra i
    Geometry::IDetector_const_sptr spDet;
    try {
      spDet = inputWS->getDetector(i);
    } catch (Kernel::Exception::NotFoundError &) {
      continue;
    }

    // Check that we aren't dealing with monitor...
    if (spDet->isMonitor())
      continue;

    // if masked detectors state is not used, masked detectors just ignored;
    bool maskDetector = spDet->isMasked();
    if (m_getIsMasked)
      *(pMasksArray + liveDetectorsCount) = maskDetector ? 1 : 0;
    else if (maskDetector)
      continue;

    // calculate the requested values;
    sp2detMap[i] = liveDetectorsCount;
    detId[liveDetectorsCount] = int32_t(spDet->getID());
    detIDMap[liveDetectorsCount] = i;
    L2[liveDetectorsCount] = spDet->getDistance(*sample);

    double polar = inputWS->detectorTwoTheta(spDet);
    double azim = spDet->getPhi();
    TwoTheta[liveDetectorsCount] = polar;
    Azimuthal[liveDetectorsCount] = azim;

    double sPhi = sin(polar);
    double ez = cos(polar);
    double ex = sPhi * cos(azim);
    double ey = sPhi * sin(azim);

    detDir[liveDetectorsCount].setX(ex);
    detDir[liveDetectorsCount].setY(ey);
    detDir[liveDetectorsCount].setZ(ez);

    // double sinTheta=sin(0.5*polar);
    // this->SinThetaSq[liveDetectorsCount]  = sinTheta*sinTheta;

    // specific code which should work and makes sense
    // for indirect instrument but may be deployed on any code with Ei property
    // defined;
    if (pEfixedArray) {
      try {
        Geometry::Parameter_sptr par = pmap.getRecursive(spDet.get(), "eFixed");
        if (par)
          Efi = par->value<double>();
      } catch (std::runtime_error &) {
      }
      // set efixed for each existing detector
      *(pEfixedArray + liveDetectorsCount) = static_cast<float>(Efi);
    }

    liveDetectorsCount++;
    if (i % div == 0)
      theProgress.report(i, "Preprocessing detectors");
  }
  targWS->logs()->addProperty<uint32_t>("ActualDetectorsNum",
                                        liveDetectorsCount, true);

  theProgress.report();
  g_log.information() << "Finished preprocessing detector locations. Found: "
                      << liveDetectorsCount << " detectors out of: " << nHist
                      << " histograms\n";
}