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

/// 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();
    }
  }
}

  /** Count the peaks from a .peaks file and compare with the workspace
   * @param outWS :: the workspace in which to place the information
   * @param filename :: path to the .peaks file
   */
  void LoadIsawPeaks::checkNumberPeaks( PeaksWorkspace_sptr outWS, std::string filename )
  {

    // Open the file
    std::ifstream in( filename.c_str() );
    std::string first;
    int NumberPeaks = 0;
    while (getline(in,first))
    {
    	if (first[0] == '3')NumberPeaks++;
    }
    if(NumberPeaks != outWS->getNumberPeaks())
    {
      g_log.error()<<"Number of peaks in file is " << NumberPeaks << " but only read "
    		  <<outWS->getNumberPeaks() << std::endl;
      throw std::length_error("Wrong number of peaks read");
    }
  }

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

/**
  @param  inname       Name of workspace containing peaks
  @param  params       optimized cell parameters
  @param  out          residuals from optimization
*/
void OptimizeLatticeForCellType::optLattice(std::string inname,
                                            std::vector<double> &params,
                                            double *out) {
  PeaksWorkspace_sptr ws = boost::dynamic_pointer_cast<PeaksWorkspace>(
      AnalysisDataService::Instance().retrieve(inname));
  const std::vector<Peak> &peaks = ws->getPeaks();
  size_t n_peaks = ws->getNumberPeaks();
  std::vector<V3D> q_vector;
  std::vector<V3D> hkl_vector;

  for (size_t i = 0; i < params.size(); i++)
    params[i] = std::abs(params[i]);
  for (size_t i = 0; i < n_peaks; i++) {
    q_vector.push_back(peaks[i].getQSampleFrame());
    hkl_vector.push_back(peaks[i].getHKL());
  }

  Mantid::API::IAlgorithm_sptr alg = createChildAlgorithm("CalculateUMatrix");
  alg->setPropertyValue("PeaksWorkspace", inname);
  alg->setProperty("a", params[0]);
  alg->setProperty("b", params[1]);
  alg->setProperty("c", params[2]);
  alg->setProperty("alpha", params[3]);
  alg->setProperty("beta", params[4]);
  alg->setProperty("gamma", params[5]);
  alg->executeAsChildAlg();

  ws = alg->getProperty("PeaksWorkspace");
  OrientedLattice latt = ws->mutableSample().getOrientedLattice();
  DblMatrix UB = latt.getUB();
  DblMatrix A = aMatrix(params);
  DblMatrix Bc = A;
  Bc.Invert();
  DblMatrix U1_B1 = UB * A;
  OrientedLattice o_lattice;
  o_lattice.setUB(U1_B1);
  DblMatrix U1 = o_lattice.getU();
  DblMatrix U1_Bc = U1 * Bc;

  for (size_t i = 0; i < hkl_vector.size(); i++) {
    V3D error = U1_Bc * hkl_vector[i] - q_vector[i] / (2.0 * M_PI);
    out[i] = error.norm2();
  }

  return;
}

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

  /** Execute the algorithm.
   */
  void SelectCellWithForm::exec()
  {
    PeaksWorkspace_sptr ws = this->getProperty("PeaksWorkspace");
    if (!ws) 
    { 
      throw std::runtime_error("Could not read the peaks workspace");
    }

    OrientedLattice o_lattice = ws->mutableSample().getOrientedLattice();
    Matrix<double> UB = o_lattice.getUB();

    bool   allowPermutations        = this->getProperty("AllowPermutations");

    if ( ! IndexingUtils::CheckUB( UB ) )
    {
       throw std::runtime_error(
             "ERROR: The stored UB is not a valid orientation matrix");
    }

    int    form_num  = this->getProperty("FormNumber");
    bool   apply     = this->getProperty("Apply");
    double tolerance = this->getProperty("Tolerance");

    ConventionalCell info = ScalarUtils::GetCellForForm( UB, form_num, allowPermutations );

    DblMatrix newUB = info.GetNewUB();

    std::string message = info.GetDescription() + " Lat Par:" +
                          IndexingUtils::GetLatticeParameterString( newUB );

    g_log.notice( std::string(message) );

    Kernel::Matrix<double>T(UB);
    T.Invert();
    T = newUB * T;
    g_log.notice() << "Transformation Matrix =  " << T.str() << std::endl;

    if ( apply )
    {
    //----------------------------------- Try to optimize(LSQ) to find lattice errors ------------------------
    //                       UB matrix may NOT have been found by unconstrained least squares optimization

     //----------------------------------------------
    o_lattice.setUB( newUB );
    std::vector<double> sigabc(6);
    DetermineErrors(sigabc,newUB,ws, tolerance);

    o_lattice.setError( sigabc[0],sigabc[1],sigabc[2],sigabc[3],sigabc[4],sigabc[5]);

    ws->mutableSample().setOrientedLattice( &o_lattice );

    std::vector<Peak> &peaks = ws->getPeaks();
    size_t n_peaks = ws->getNumberPeaks();

    int    num_indexed   = 0;
    double average_error = 0.0;
    std::vector<V3D> miller_indices;
    std::vector<V3D> q_vectors;
    for ( size_t i = 0; i < n_peaks; i++ )
    {
      q_vectors.push_back( peaks[i].getQSampleFrame() );
    }

    num_indexed = IndexingUtils::CalculateMillerIndices( newUB, q_vectors,
                                                         tolerance,
                                                         miller_indices,
                                                         average_error );

    for ( size_t i = 0; i < n_peaks; i++ )
    {
      peaks[i].setHKL( miller_indices[i] );
    }

    // Tell the user what happened.
    g_log.notice() << "Re-indexed the peaks with the new UB. " << std::endl;
    g_log.notice() << "Now, " << num_indexed << " are indexed with average error " << average_error << std::endl;

    // Save output properties

    this->setProperty("NumIndexed", num_indexed);
    this->setProperty("AverageError", average_error);
    }
  }

  /** Execute the algorithm.
   */
  void CalculateUMatrix::exec()
  {
    double a=this->getProperty("a");
    double b=this->getProperty("b");
    double c=this->getProperty("c");
    double alpha=this->getProperty("alpha");
    double beta=this->getProperty("beta");
    double gamma=this->getProperty("gamma");
    OrientedLattice o(a,b,c,alpha,beta,gamma);
    Matrix<double> B=o.getB();

    double H,K,L;

    PeaksWorkspace_sptr ws;
    ws = AnalysisDataService::Instance().retrieveWS<PeaksWorkspace>(this->getProperty("PeaksWorkspace") );
    if (!ws) throw std::runtime_error("Problems reading the peaks workspace");

    size_t nIndexedpeaks=0;
    bool found2nc=false;
    V3D old(0,0,0);
    Matrix<double> Hi(4,4),Si(4,4),HS(4,4),zero(4,4);
    for (int i=0;i<ws->getNumberPeaks();i++)
    {
      Peak p=ws->getPeaks()[i];
      H=p.getH();
      K=p.getK();
      L=p.getL();
      if(H*H+K*K+L*L>0)
      {
        nIndexedpeaks++;
        if (!found2nc)
        {
          if (nIndexedpeaks==1)
          {
            old=V3D(H,K,L);
          }
          else
          {
            if (!old.coLinear(V3D(0,0,0),V3D(H,K,L))) found2nc=true;
          }
        }
        V3D Qhkl=B*V3D(H,K,L);
        Hi[0][0]=0.;        Hi[0][1]=-Qhkl.X(); Hi[0][2]=-Qhkl.Y(); Hi[0][3]=-Qhkl.Z();
        Hi[1][0]=Qhkl.X();  Hi[1][1]=0.;        Hi[1][2]=Qhkl.Z();  Hi[1][3]=-Qhkl.Y();
        Hi[2][0]=Qhkl.Y();  Hi[2][1]=-Qhkl.Z(); Hi[2][2]=0.;        Hi[2][3]=Qhkl.X();
        Hi[3][0]=Qhkl.Z();  Hi[3][1]=Qhkl.Y();  Hi[3][2]=-Qhkl.X(); Hi[3][3]=0.;

        V3D Qgon=p.getQSampleFrame();
        Si[0][0]=0.;        Si[0][1]=-Qgon.X(); Si[0][2]=-Qgon.Y(); Si[0][3]=-Qgon.Z();
        Si[1][0]=Qgon.X();  Si[1][1]=0.;        Si[1][2]=-Qgon.Z(); Si[1][3]=Qgon.Y();
        Si[2][0]=Qgon.Y();  Si[2][1]=Qgon.Z();  Si[2][2]=0.;        Si[2][3]=-Qgon.X();
        Si[3][0]=Qgon.Z();  Si[3][1]=-Qgon.Y(); Si[3][2]=Qgon.X(); Si[3][3]=0.;

        HS+=(Hi*Si);
      }
    }
    //check if enough peaks are indexed or if HS is 0
    if ((nIndexedpeaks<2) || (found2nc==false)) throw std::invalid_argument("Less then two non-colinear peaks indexed");
    if (HS==zero) throw std::invalid_argument("Something really bad happened");

    Matrix<double> Eval;
    Matrix<double> Diag;
    HS.Diagonalise(Eval,Diag);
    Eval.sortEigen(Diag);
    Mantid::Kernel::Quat qR(Eval[0][0],Eval[1][0],Eval[2][0],Eval[3][0]);//the first column corresponds to the highest eigenvalue
    DblMatrix U(qR.getRotation());
    o.setU(U);

    ws->mutableSample().setOrientedLattice(new OrientedLattice(o));

  }

    /** 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.

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

/**
  @param  inname       Name of Filename containing peaks
  @param  cell_type    cell type to optimize
  @param  params       optimized cell parameters
  @return  chisq of optimization
*/
double OptimizeLatticeForCellType::optLatticeSum(std::string inname,
                                                 std::string cell_type,
                                                 std::vector<double> &params) {
  std::vector<double> lattice_parameters;
  lattice_parameters.assign(6, 0);
  if (cell_type == ReducedCell::CUBIC()) {
    lattice_parameters[0] = params[0];
    lattice_parameters[1] = params[0];
    lattice_parameters[2] = params[0];

    lattice_parameters[3] = 90;
    lattice_parameters[4] = 90;
    lattice_parameters[5] = 90;
  } else if (cell_type == ReducedCell::TETRAGONAL()) {
    lattice_parameters[0] = params[0];
    lattice_parameters[1] = params[0];
    lattice_parameters[2] = params[1];

    lattice_parameters[3] = 90;
    lattice_parameters[4] = 90;
    lattice_parameters[5] = 90;
  } else if (cell_type == ReducedCell::ORTHORHOMBIC()) {
    lattice_parameters[0] = params[0];
    lattice_parameters[1] = params[1];
    lattice_parameters[2] = params[2];

    lattice_parameters[3] = 90;
    lattice_parameters[4] = 90;
    lattice_parameters[5] = 90;
  } else if (cell_type == ReducedCell::RHOMBOHEDRAL()) {
    lattice_parameters[0] = params[0];
    lattice_parameters[1] = params[0];
    lattice_parameters[2] = params[0];

    lattice_parameters[3] = params[1];
    lattice_parameters[4] = params[1];
    lattice_parameters[5] = params[1];
  } else if (cell_type == ReducedCell::HEXAGONAL()) {
    lattice_parameters[0] = params[0];
    lattice_parameters[1] = params[0];
    lattice_parameters[2] = params[1];

    lattice_parameters[3] = 90;
    lattice_parameters[4] = 90;
    lattice_parameters[5] = 120;
  } else if (cell_type == "Monoclinic ( a unique )") {
    lattice_parameters[0] = params[0];
    lattice_parameters[1] = params[1];
    lattice_parameters[2] = params[2];

    lattice_parameters[3] = params[3];
    lattice_parameters[4] = 90;
    lattice_parameters[5] = 90;
  } else if (cell_type == ReducedCell::MONOCLINIC() ||
             cell_type == "Monoclinic ( b unique )") {
    lattice_parameters[0] = params[0];
    lattice_parameters[1] = params[1];
    lattice_parameters[2] = params[2];

    lattice_parameters[3] = 90;
    lattice_parameters[4] = params[3];
    lattice_parameters[5] = 90;
  } else if (cell_type == "Monoclinic ( c unique )") {
    lattice_parameters[0] = params[0];
    lattice_parameters[1] = params[1];
    lattice_parameters[2] = params[2];

    lattice_parameters[3] = 90;
    lattice_parameters[4] = 90;
    lattice_parameters[5] = params[3];
  } else if (cell_type == ReducedCell::TRICLINIC()) {
    lattice_parameters[0] = params[0];
    lattice_parameters[1] = params[1];
    lattice_parameters[2] = params[2];

    lattice_parameters[3] = params[3];
    lattice_parameters[4] = params[4];
    lattice_parameters[5] = params[5];
  }

  PeaksWorkspace_sptr ws = boost::dynamic_pointer_cast<PeaksWorkspace>(
      AnalysisDataService::Instance().retrieve(inname));
  size_t n_peaks = ws->getNumberPeaks();
  double *out = new double[n_peaks];
  optLattice(inname, lattice_parameters, out);
  double ChiSqTot = 0;
  for (size_t i = 0; i < n_peaks; i++)
    ChiSqTot += out[i];
  delete[] out;
  return ChiSqTot;
}

  /** Execute the algorithm.
   */
  void FindUBUsingIndexedPeaks::exec()
  {
    PeaksWorkspace_sptr ws;
    ws = boost::dynamic_pointer_cast<PeaksWorkspace>(
         AnalysisDataService::Instance().retrieve(this->getProperty("PeaksWorkspace")) );

    if (!ws)
    {
      throw std::runtime_error("Could not read the peaks workspace");
    }

    std::vector<Peak> peaks = ws->getPeaks();
    size_t n_peaks = ws->getNumberPeaks();

    std::vector<V3D>  q_vectors;
    std::vector<V3D>  hkl_vectors;

    q_vectors.reserve( n_peaks );
    hkl_vectors.reserve( n_peaks );

    size_t indexed_count = 0;
    for ( size_t i = 0; i < n_peaks; i++ )
    {
      V3D hkl( peaks[i].getH(), peaks[i].getK(), peaks[i].getL() );    // ##### KEEP
      if ( IndexingUtils::ValidIndex( hkl, 1.0 ) )    // use tolerance == 1 to 
                                                      // just check for (0,0,0) 
      {
        q_vectors.push_back( peaks[i].getQSampleFrame() );
        V3D miller_ind( round(hkl[0]), round(hkl[1]), round(hkl[2]) );
        hkl_vectors.push_back( V3D(miller_ind) );
        indexed_count++;
      }
    }

    if ( indexed_count < 3 ) 
    { 
      throw std::runtime_error(
            "At least three linearly independent indexed peaks are needed.");
    }

    Matrix<double> UB(3,3,false);
    double error = IndexingUtils::Optimize_UB( UB, hkl_vectors, q_vectors );

    std::cout << "Error = " << error << std::endl;
    std::cout << "UB = " << UB << std::endl;

    if ( ! IndexingUtils::CheckUB( UB ) ) // UB not found correctly
    {
      g_log.notice( std::string(
         "Found Invalid UB...peaks used might not be linearly independent") );
      g_log.notice( std::string(
         "UB NOT SAVED.") );
    }
    else                                 // tell user how many would be indexed
    {                                    // from the full list of peaks, and  
      q_vectors.clear();                 // save the UB in the sample
      q_vectors.reserve( n_peaks );
      for ( size_t i = 0; i < n_peaks; i++ )
      {
        q_vectors.push_back( peaks[i].getQSampleFrame() );
      }
      double tolerance = 0.1;
      int num_indexed = IndexingUtils::NumberIndexed(UB, q_vectors, tolerance);

      char logInfo[200];
      sprintf( logInfo,
               std::string("New UB will index %1d Peaks out of %1d with tolerance %5.3f").c_str(),
               num_indexed, n_peaks, tolerance);
      g_log.notice( std::string(logInfo) );

      OrientedLattice o_lattice;
      o_lattice.setUB( UB );
      double calc_a = o_lattice.a();
      double calc_b = o_lattice.b();
      double calc_c = o_lattice.c();
      double calc_alpha = o_lattice.alpha();
      double calc_beta  = o_lattice.beta();
      double calc_gamma = o_lattice.gamma();
                                       // Show the modified lattice parameters
      sprintf( logInfo, 
               std::string("Lattice Parameters: %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f").c_str(),
               calc_a, calc_b, calc_c, calc_alpha, calc_beta, calc_gamma);
      g_log.notice( std::string(logInfo) );

      ws->mutableSample().setOrientedLattice( new OrientedLattice(o_lattice) );
    }
  }