C++ PeaksWorkspace_sptr::getPeak方法代码示例

/// Fills possibleHKLs with all HKLs from the supplied PeaksWorkspace.
void PredictPeaks::fillPossibleHKLsUsingPeaksWorkspace(
    const PeaksWorkspace_sptr &peaksWorkspace,
    std::vector<V3D> &possibleHKLs) const {
  possibleHKLs.clear();
  possibleHKLs.reserve(peaksWorkspace->getNumberPeaks());

  bool roundHKL = getProperty("RoundHKL");

  /* Q is at the end multiplied with the factor determined in the
   * constructor (-1 for crystallography, 1 otherwise). So to avoid
   * "flippling HKLs" when it's not required, the HKLs of the input
   * workspace are also multiplied by the factor that is appropriate
   * for the convention stored in the workspace.
   */
  double peaks_q_convention_factor =
      get_factor_for_q_convention(peaksWorkspace->getConvention());

  for (int i = 0; i < static_cast<int>(peaksWorkspace->getNumberPeaks()); ++i) {
    IPeak &p = peaksWorkspace->getPeak(i);
    // Get HKL from that peak
    V3D hkl = p.getHKL() * peaks_q_convention_factor;

    if (roundHKL)
      hkl.round();

    possibleHKLs.push_back(hkl);
  } // for each hkl in the workspace
}

/**
 * Updates the map from run number to GoniometerMatrix
 *
 * @param Peaks    The PeaksWorkspace whose peaks contain the run numbers
 *                   along with the corresponding GoniometerMatrix
 *
 * @param OptRuns  A '/' separated "list" of run numbers to include in the
 *                  map. This string must also start and end with a '/'
 *
 * @param Res      The resultant map.
 */
void PeakHKLErrors::getRun2MatMap(
    PeaksWorkspace_sptr &Peaks, const std::string &OptRuns,
    std::map<int, Mantid::Kernel::Matrix<double>> &Res) const {

  for (int i = 0; i < Peaks->getNumberPeaks(); ++i) {
    Geometry::IPeak &peak_old = Peaks->getPeak(i);

    int runNum = peak_old.getRunNumber();
    std::string runNumStr = std::to_string(runNum);
    size_t N = OptRuns.find("/" + runNumStr + "/");
    if (N < OptRuns.size()) {
      double chi =
          getParameter("chi" + boost::lexical_cast<std::string>(runNumStr));
      double phi =
          getParameter("phi" + boost::lexical_cast<std::string>(runNumStr));
      double omega =
          getParameter("omega" + boost::lexical_cast<std::string>(runNumStr));
      Mantid::Geometry::Goniometer uniGonio;
      uniGonio.makeUniversalGoniometer();
      uniGonio.setRotationAngle("phi", phi);
      uniGonio.setRotationAngle("chi", chi);
      uniGonio.setRotationAngle("omega", omega);
      Res[runNum] = uniGonio.getR();
    }
  }
}

/**
 * Checks that a PeaksWorkspace has only one run.
 *
 * @param Peaks   The PeaksWorkspace
 * @param GoniometerMatrix  the goniometer matrix for the run
 */
bool GoniometerAnglesFromPhiRotation::CheckForOneRun(
    const PeaksWorkspace_sptr &Peaks,
    Kernel::Matrix<double> &GoniometerMatrix) const {

  int RunNumber = -1;
  for (int peak = 0; peak < Peaks->getNumberPeaks(); peak++) {
    int thisRunNum = Peaks->getPeak(peak).getRunNumber();
    GoniometerMatrix = Peaks->getPeak(peak).getGoniometerMatrix();

    if (RunNumber < 0)

      RunNumber = thisRunNum;

    else if (thisRunNum != RunNumber)

      return false;
  }

  return true;
}

/**
 * Creates a new parameterized instrument for which the parameter values can be
 *changed
 *
 * @param Peaks - a PeaksWorkspace used to get the original instrument.  The
 *instrument from the 0th peak is
 *                the one that is used.
 *
 * NOTE: All the peaks in the PeaksWorkspace must use the same instrument.
 */
boost::shared_ptr<Geometry::Instrument>
PeakHKLErrors::getNewInstrument(PeaksWorkspace_sptr Peaks) const {
  Geometry::Instrument_const_sptr instSave = Peaks->getPeak(0).getInstrument();
  auto pmap = boost::make_shared<Geometry::ParameterMap>();

  if (!instSave) {
    g_log.error(" Peaks workspace does not have an instrument");
    throw std::invalid_argument(" Not all peaks have an instrument");
  }

  if (!hasParameterMap) {
    pmapSv = instSave->getParameterMap();
    hasParameterMap = true;
    if (!instSave->isParametrized()) {

      boost::shared_ptr<Geometry::Instrument> instClone(instSave->clone());
      auto Pinsta = boost::make_shared<Geometry::Instrument>(instSave, pmap);

      instChange = Pinsta;
      IComponent_const_sptr sample = instChange->getSample();
      sampPos = sample->getRelativePos();
    } else // catch(... )
    {
      auto P1 = boost::make_shared<Geometry::Instrument>(
          instSave->baseInstrument(), instSave->makeLegacyParameterMap());
      instChange = P1;
      IComponent_const_sptr sample = instChange->getSample();
      sampPos = sample->getRelativePos();
    }
  }

  if (!instChange) {
    g_log.error("Cannot 'clone' instrument");
    throw std::logic_error("Cannot clone instrument");
  }
  //------------------"clone" orig instruments pmap -------------------

  cLone(pmap, instSave, pmapSv);
  V3D sampOffsets(getParameter("SampleXOffset"), getParameter("SampleYOffset"),
                  getParameter("SampleZOffset"));

  IComponent_const_sptr sample = instChange->getSample();
  pmap->addPositionCoordinate(sample.get(), std::string("x"),
                              sampPos.X() + sampOffsets.X());
  pmap->addPositionCoordinate(sample.get(), std::string("y"),
                              sampPos.Y() + sampOffsets.Y());
  pmap->addPositionCoordinate(sample.get(), std::string("z"),
                              sampPos.Z() + sampOffsets.Z());

  return instChange;
}

/// Fills possibleHKLs with all HKLs from the supplied PeaksWorkspace.
void PredictPeaks::fillPossibleHKLsUsingPeaksWorkspace(
    const PeaksWorkspace_sptr &peaksWorkspace,
    std::vector<V3D> &possibleHKLs) const {
  possibleHKLs.clear();
  possibleHKLs.reserve(peaksWorkspace->getNumberPeaks());

  bool roundHKL = getProperty("RoundHKL");

  for (int i = 0; i < static_cast<int>(peaksWorkspace->getNumberPeaks()); ++i) {
    IPeak &p = peaksWorkspace->getPeak(i);
    // Get HKL from that peak
    V3D hkl = p.getHKL();

    if (roundHKL)
      hkl.round();

    possibleHKLs.push_back(hkl);
  } // for each hkl in the workspace
}

void PeakHKLErrors::functionDeriv1D(Jacobian *out, const double *xValues,
                                    const size_t nData) {
  PeaksWorkspace_sptr Peaks =
      AnalysisDataService::Instance().retrieveWS<PeaksWorkspace>(
          PeakWorkspaceName);
  boost::shared_ptr<Geometry::Instrument> instNew = getNewInstrument(Peaks);

  const DblMatrix &UB = Peaks->sample().getOrientedLattice().getUB();
  DblMatrix UBinv(UB);
  UBinv.Invert();
  UBinv /= 2 * M_PI;

  double GonRotx = getParameter("GonRotx");
  double GonRoty = getParameter("GonRoty");
  double GonRotz = getParameter("GonRotz");
  Matrix<double> InvGonRotxMat = RotationMatrixAboutRegAxis(GonRotx, 'x');
  Matrix<double> InvGonRotyMat = RotationMatrixAboutRegAxis(GonRoty, 'y');
  Matrix<double> InvGonRotzMat = RotationMatrixAboutRegAxis(GonRotz, 'z');
  Matrix<double> GonRot = InvGonRotxMat * InvGonRotyMat * InvGonRotzMat;

  InvGonRotxMat.Invert();
  InvGonRotyMat.Invert();
  InvGonRotzMat.Invert();

  std::map<int, Kernel::Matrix<double>> RunNums2GonMatrix;
  getRun2MatMap(Peaks, OptRuns, RunNums2GonMatrix);

  g_log.debug()
      << "----------------------------Derivative------------------------\n";

  V3D samplePosition = instNew->getSample()->getPos();
  IPeak &ppeak = Peaks->getPeak(0);
  double L0 = ppeak.getL1();
  double velocity = (L0 + ppeak.getL2()) / ppeak.getTOF();

  double K =
      2 * M_PI / ppeak.getWavelength() / velocity; // 2pi/lambda = K* velocity
  V3D beamDir = instNew->getBeamDirection();

  size_t paramNums[] = {parameterIndex(std::string("SampleXOffset")),
                        parameterIndex(std::string("SampleYOffset")),
                        parameterIndex(std::string("SampleZOffset"))};

  for (size_t i = 0; i < nData; i += 3) {
    int peakNum = boost::math::iround(xValues[i]);
    IPeak &peak_old = Peaks->getPeak(peakNum);
    Peak peak = createNewPeak(peak_old, instNew, 0, peak_old.getL1());

    int runNum = peak_old.getRunNumber();
    std::string runNumStr = std::to_string(runNum);

    for (int kk = 0; kk < static_cast<int>(nParams()); kk++) {
      out->set(i, kk, 0.0);
      out->set(i + 1, kk, 0.0);
      out->set(i + 2, kk, 0.0);
    }

    double chi, phi, omega;
    size_t chiParamNum, phiParamNum, omegaParamNum;

    size_t N = OptRuns.find("/" + runNumStr);
    if (N < OptRuns.size()) {
      chi = getParameter("chi" + (runNumStr));
      phi = getParameter("phi" + (runNumStr));
      omega = getParameter("omega" + (runNumStr));

      peak.setGoniometerMatrix(GonRot * RunNums2GonMatrix[runNum]);

      chiParamNum = parameterIndex("chi" + (runNumStr));
      phiParamNum = parameterIndex("phi" + (runNumStr));
      omegaParamNum = parameterIndex("omega" + (runNumStr));
    } else {

      Geometry::Goniometer Gon(peak.getGoniometerMatrix());
      std::vector<double> phichiOmega = Gon.getEulerAngles("YZY");
      chi = phichiOmega[1];
      phi = phichiOmega[2];
      omega = phichiOmega[0];
      // peak.setGoniometerMatrix( GonRot*Gon.getR());
      chiParamNum = phiParamNum = omegaParamNum = nParams() + 10;
      peak.setGoniometerMatrix(GonRot * peak.getGoniometerMatrix());
    }
    V3D sampOffsets(getParameter("SampleXOffset"),
                    getParameter("SampleYOffset"),
                    getParameter("SampleZOffset"));
    peak.setSamplePos(peak.getSamplePos() + sampOffsets);
    // NOTE:Use getQLabFrame except for below.
    // For parameters the getGoniometerMatrix should remove GonRot, for derivs
    // wrt GonRot*, wrt chi*,phi*,etc.

    // Deriv wrt chi phi and omega
    if (phiParamNum < nParams()) {
      Matrix<double> chiMatrix = RotationMatrixAboutRegAxis(chi, 'z');
      Matrix<double> phiMatrix = RotationMatrixAboutRegAxis(phi, 'y');
      Matrix<double> omegaMatrix = RotationMatrixAboutRegAxis(omega, 'y');

      Matrix<double> dchiMatrix = DerivRotationMatrixAboutRegAxis(chi, 'z');
      Matrix<double> dphiMatrix = DerivRotationMatrixAboutRegAxis(phi, 'y');
      Matrix<double> domegaMatrix = DerivRotationMatrixAboutRegAxis(omega, 'y');

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

    /** Executes the algorithm
    *
    *  @throw runtime_error Thrown if algorithm cannot execute
    */
    void IndexSXPeaks::exec()
    {
      using namespace Mantid::DataObjects;
      std::vector<int> peakindices = getProperty("PeakIndices");

      PeaksWorkspace_sptr ws = boost::dynamic_pointer_cast<PeaksWorkspace>(
        AnalysisDataService::Instance().retrieve(this->getProperty("PeaksWorkspace")) );

      // Need a least two peaks
      std::size_t npeaks=peakindices.size();
      if (npeaks > size_t(ws->getNumberPeaks()))
      {
        throw std::runtime_error("Cannot have more peaks indices than actual peaks");
      }
      if (npeaks == 1 || ws->getNumberPeaks() < 2)
      {
        throw std::runtime_error("At least 2 peaks are required for this algorithm to work");
      }
      if (npeaks == 0)
      {
        //If the user provides no peaks we default to use all the available peaks.
        npeaks = ws->getNumberPeaks();
        peakindices.reserve(npeaks);
        for(int i = 1; i <= int(npeaks); i++) //create indexes corresponding to all peak indexes
        {
          peakindices.push_back(i);
        }
        g_log.information("No peak indexes provided. Algorithm will use all peaks in the workspace for the calculation.");
      }

      //Get the effective unit cell
      double a = getProperty("a");
      double b = getProperty("b");
      double c = getProperty("c");
      double alpha = getProperty("alpha");
      double beta = getProperty("beta");
      double gamma = getProperty("gamma");

      std::vector<int> extents = getProperty("SearchExtents");
      if(extents.size() != 6)
      {
        std::stringstream stream;
        stream << "Expected 6 elements for the extents. Got: " << extents.size();
        throw std::runtime_error(stream.str());      
      }

      // Create the Unit-Cell.
      Mantid::Geometry::UnitCell unitcell(a, b, c, alpha, beta, gamma);

      std::vector<PeakCandidate> peaks;

      for (std::size_t i=0;i<npeaks;i++)
      { 
        int row=peakindices[i]-1;
        if(row < 0)
        {
          throw std::runtime_error("Cannot have a peak index < 0.");
        }
        IPeak& peak = ws->getPeak(row);
        V3D Qs = peak.getQSampleFrame() / (2.0 * M_PI);
        peaks.push_back(PeakCandidate(Qs[0], Qs[1], Qs[2]));
      }

      //Sanity check the generated peaks.
      validateNotColinear(peaks);

      //Generate HKL possibilities for each peak.
      double dtol= getProperty("dTolerance");
      Progress prog(this,0.0,1.0,4);
      for (int h=extents[0]; h<extents[1]; h++)
      {
        for (int k=extents[2]; k<extents[3]; k++)
        {
          for (int l=extents[4]; l<extents[5]; l++)
          {
            double dspacing=unitcell.d(h,k,l); //Create a fictional d spacing
            for (std::size_t p=0;p<npeaks;p++)
            {
              double dSpacingPeaks = peaks[p].getdSpacing();
              if (std::abs(dspacing-dSpacingPeaks)<dtol)
                peaks[p].addHKL(h,k,l); // If the peak position and the fictional d spacing are within tolerance, add it
            }
          }
        }
      }
      prog.report(); //1st Progress report.

      cullHKLs(peaks, unitcell);

      prog.report(); //2nd progress report.
      peaks[0].setFirst(); //On the first peak, now only the first candidate hkl is considered, others are erased,
      //This means the design space of possible peak-hkl alignments has been reduced, will improve future refinements.

      cullHKLs(peaks, unitcell);
      prog.report(); //3rd progress report.

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

int PeakIntegration::fitneighbours(int ipeak, std::string det_name, int x0,
                                   int y0, int idet, double qspan,
                                   PeaksWorkspace_sptr &Peaks,
                                   const detid2index_map &pixel_to_wi) {
  UNUSED_ARG(ipeak);
  UNUSED_ARG(det_name);
  UNUSED_ARG(x0);
  UNUSED_ARG(y0);
  Geometry::IPeak &peak = Peaks->getPeak(ipeak);
  // Number of slices
  int TOFmax = 0;

  IAlgorithm_sptr slice_alg = createChildAlgorithm("IntegratePeakTimeSlices");
  slice_alg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputW);
  std::ostringstream tab_str;
  tab_str << "LogTable" << ipeak;

  slice_alg->setPropertyValue("OutputWorkspace", tab_str.str());
  slice_alg->setProperty<PeaksWorkspace_sptr>("Peaks", Peaks);
  slice_alg->setProperty("PeakIndex", ipeak);
  slice_alg->setProperty("PeakQspan", qspan);

  int nPixels = std::max<int>(0, getProperty("NBadEdgePixels"));

  slice_alg->setProperty("NBadEdgePixels", nPixels);
  slice_alg->executeAsChildAlg();
  Mantid::API::MemoryManager::Instance().releaseFreeMemory();

  MantidVec &Xout = outputW->dataX(idet);
  MantidVec &Yout = outputW->dataY(idet);
  MantidVec &Eout = outputW->dataE(idet);
  TableWorkspace_sptr logtable = slice_alg->getProperty("OutputWorkspace");

  peak.setIntensity(slice_alg->getProperty("Intensity"));
  peak.setSigmaIntensity(slice_alg->getProperty("SigmaIntensity"));

  TOFmax = static_cast<int>(logtable->rowCount());
  for (int iTOF = 0; iTOF < TOFmax; iTOF++) {
    Xout[iTOF] = logtable->getRef<double>(std::string("Time"), iTOF);
    if (m_IC) // Ikeda-Carpenter fit
    {
      Yout[iTOF] = logtable->getRef<double>(std::string("TotIntensity"), iTOF);
      Eout[iTOF] =
          logtable->getRef<double>(std::string("TotIntensityError"), iTOF);
    } else {
      Yout[iTOF] = logtable->getRef<double>(std::string("ISAWIntensity"), iTOF);
      Eout[iTOF] =
          logtable->getRef<double>(std::string("ISAWIntensityError"), iTOF);
    }
  }

  outputW->getSpectrum(idet)->clearDetectorIDs();
  // Find the pixel ID at that XY position on the rectangular detector
  int pixelID = peak.getDetectorID(); // det->getAtXY(x0,y0)->getID();

  // Find the corresponding workspace index, if any
  auto wiEntry = pixel_to_wi.find(pixelID);
  if (wiEntry != pixel_to_wi.end()) {
    size_t wi = wiEntry->second;
    // Set detectorIDs
    outputW->getSpectrum(idet)
        ->addDetectorIDs(inputW->getSpectrum(wi)->getDetectorIDs());
  }

  return TOFmax - 1;
}