C++ TParticle::Momentum方法代码示例

void testTDime(Int_t nev = 100) {

  gSystem->Load("libEVGEN");
  gSystem->Load("libTDime");
  gSystem->Load("libdime");

  TDime* dime = new TDime();
  dime->SetEnergyCMS(7000.0);
  dime->SetYRange(-2.0, 2.0);   // Set rapidity range of mesons
  dime->SetMinPt(0.1);          // Minimum pT of mesons
  dime->Initialize();

  // (pi+pi-) histograms
  TH1* hM = new TH1D("hM", "DIME #pi^{+}#pi^{-};M_{#pi^{+}#pi^{-}} #[]{GeV/#it{c}^{2}}", 100,  0.0, 5.0);

  TClonesArray* particles = new TClonesArray("TParticle", 6);
  TParticle* part = NULL;
  TLorentzVector v[2];
  TLorentzVector vSum;

  // Event loop
  for (Int_t i = 0; i < nev; ++i) {

    dime->GenerateEvent();
    Int_t np = dime->ImportParticles(particles, "All");
    printf("\n DIME Event %d: Imported %3d particles \n", i, np);

    Int_t nPrimary = 0;

    // Loop over pion (j = 4,5) tracks
    for (Int_t j = 4; j < 6; ++j) {
      part = (TParticle*) particles->At(j); // Choose the particle
      part->Print();
      part->Momentum(v[nPrimary]);          // Copy content to v
      nPrimary++;
    }
    //particles.Clear();

    // 4-vector sum
    vSum = v[0] + v[1];

    // Fill pi+pi- histograms
    hM->Fill(vSum.M());
  }

  // Save plots as pdf
  hM->Draw();    c1->SaveAs("massTDime.pdf");

}

void Kin2Txt() {
  AliRunLoader* rl = AliRunLoader::Open("galice.root");
  rl->LoadKinematics();
  rl->LoadHeader();

  TH1* hM  = new TH1D("hM",  "TreeK;M#(){#pi^{+}#pi^{-}} #(){GeV/#it{c}^{2}}", 100, 0., 2.);
  TH1* hPt = new TH1D("hPt", "TreeK;P_{T}#(){#pi^{+}#pi^{-}} #(){GeV/#it{c}}", 100, 0., 1.);
  TH1* hY  = new TH1D("hY",  "TreeK;Y#(){#pi^{+}#pi^{-}}",                     160,-8., 8.);
  std::ofstream ofs("rho0.txt");
  
  AliStack *stack = NULL;
  TParticle *part = NULL;
  TLorentzVector v[2], vSum;
  for (Int_t i=0, n(rl->GetNumberOfEvents()); i<n; ++i) {
    rl->GetEvent(i);
    stack = rl->Stack();

    Int_t nPrimary(0);
    for (Int_t j(0), m(stack->GetNtrack()); j<m; ++j) {
      part = stack->Particle(j);
      if (part->IsPrimary())
	part->Momentum(v[nPrimary++]);
    }
    if (nPrimary != 2) {
      Printf("Error: nPrimary=%d != 2", nPrimary);
      continue;
    }
    vSum = v[0] + v[1];
    hM->Fill(vSum.M());
    hPt->Fill(vSum.Perp());
    hY->Fill(vSum.Rapidity());
    ofs << std::fixed << std::setprecision(4)
	<< vSum.M() << " " << vSum.Perp() << " " << vSum.Rapidity() << " "
	<< v[0].Eta() << " " << v[0].Px() << " " << v[0].Py() << " " << v[0].Pz() << " "
	<< v[1].Eta() << " " << v[1].Px() << " " << v[1].Py() << " " << v[1].Pz()
	<< std::endl;
  }
  hM->Draw();
  c1->SaveAs("TreeK.pdf(");
  hPt->Draw();
  c1->SaveAs("TreeK.pdf");
  hY->Draw();
  c1->SaveAs("TreeK.pdf)");
}

void plots() {

  gSystem->Load("libEVGEN"); // Needs to be!

  AliRunLoader* rl = AliRunLoader::Open("galice.root");
  rl->LoadKinematics();
  rl->LoadHeader();

  // 4pi histograms
  TH1* hM    = new TH1D("hM",  "DIME #rho#rho;M_{4#pi} #(){GeV/#it{c}^{2}}", 100,  1.0, 3.0);
  TH1* hPt   = new TH1D("hPt", "DIME #rho#rho;p_{T}#(){4#pi} #(){GeV/#it{c}}", 100,  0.0, 3.0);

  // pi+- histograms
  TH1* hPt1   = new TH1D("hPt1", "DIME #rho#rho;p_{T}#(){#pi^{#pm}} #(){Gev/#it{c}}",  100, 0.0, 3.0);

  AliStack* stack = NULL;
  TParticle* part = NULL;
  TLorentzVector v[4];
  TLorentzVector vSum;

  // Loop over events
  for (Int_t i = 0; i < rl->GetNumberOfEvents(); ++i) {

    rl->GetEvent(i);
    stack = rl->Stack();
    Int_t nPrimary = 0;

    // Loop over all particles
    for (Int_t j = 0; j < stack->GetNtrack(); ++j) {
      part = stack->Particle(j);           // Get particle
      part->Print();                       // Print contents

      if (abs(part->GetPdgCode()) == 211   // Is pi+ or pi-
          & part->GetStatusCode() == 1     // Is stable final state
          & stack->IsPhysicalPrimary(j)) { // Is from generator level

        part->Momentum(v[nPrimary]);       // Set content of v
        ++nPrimary;
      }
    }
    if (nPrimary != 4) {
      printf("Error: nPrimary=%d != 4 \n", nPrimary);
      continue;
    }

    // 4-vector sum
    vSum = v[0] + v[1] + v[2] + v[3];

    // Fill 4pi histograms
    hM->Fill(vSum.M());
    hPt->Fill(vSum.Perp());

    // Fill pi+- histograms
    for (Int_t k = 0; k < 4; ++k) {
      hPt1->Fill(v[k].Perp());
    }
    printf("\n");
  }

  // Save plots as pdf
  hM->Draw();    c1->SaveAs("plotM.pdf");
  hPt->Draw();   c1->SaveAs("plotPt.pdf");
  hPt1->Draw();  c1->SaveAs("plotPt1.pdf");

}

//.........这里部分代码省略.........
  plots.v_Jet_8_DR = -99;
 
  plots.v_numEle = -99;
  plots.v_numMu = -99;
  plots.v_numJets = -99;
  plots.v_totNumJets = -99;
  
  tree->GetEntry (evt) ;
  
  tagJets -> Clear () ;  
  otherJets -> Clear () ;    
  
  
  
    //---- check if signal ----
  if (if_signal && (IdEvent!=123 && IdEvent!=124)) continue;
  plots.analyzed++;
  
  
  //---- MC data ----
  std::vector<TLorentzVector*> MCJets ;
  
  TLorentzVector* MCJets_temp[6] ;
  int counter = 0;
  
  if (if_signal && (IdEvent==123 || IdEvent==124)){
   for(int ii=0; ii<9; ii++){
//     if (ii==0 || ii==1){
//     if (ii!=0 && ii!=1 && ii!=2 && ii!=3 && ii!=6){
    if (ii!=2 && ii!=3 && ii!=6){
     TParticle* myparticle = (TParticle*) HiggsParticle->At(ii);
//      std::cerr << "pdg = " << ii << " = " << myparticle->GetPdgCode() << std::endl;
     MCJets_temp[counter] = new TLorentzVector;
     myparticle->Momentum(*(MCJets_temp[counter]));
     MCJets.push_back((MCJets_temp[counter]));
     counter++;
    }
   }
  }
  

  
   //---- find Tagging Jets ----
 
  double m_jetPtMin = 15.;
  double m_jetEtaMax = 5.;
  double m_jetDEtaMin = -1;
  double m_jetMjjMin = -1;
 
  

  std::vector<myJet> goodJets ;

//   std::cerr << std::endl << std::endl << std::endl << std::endl << std::endl;
  
  for (int l=0; l<otherJets_temp->GetEntries (); l++ ){
   TLorentzVector* jet_temp = (TLorentzVector*) otherJets_temp->At(l);
   if (jet_temp->Pt()<m_jetPtMin) continue;
   //---- Eta max threshold ----
   if (jet_temp->Eta()>m_jetEtaMax) continue;
   //---- pt min threshold ----
   myJet dummy (jet_temp, 0, 0) ;
   goodJets.push_back (dummy) ;
  }
  
//   for (int gg=0; gg<goodJets.size(); gg++ ) std::cerr << " goodJets[" << gg << "] = " << &(goodJets.at(gg)) << std::endl;

void testsl() {
  gSystem->Load("libStarLight");
  gSystem->Load("libAliStarLight");

  TStarLight* sl = new TStarLight("starlight generator", "title", "");

  sl->SetParameter("baseFileName = slight	#suite of output files will be saved with this base name");
  sl->SetParameter("BEAM_1_Z = 82    #Z of projectile");
  sl->SetParameter("BEAM_1_A = 208   #A of projectile");
  sl->SetParameter("BEAM_2_Z = 82   #Z of target");
  sl->SetParameter("BEAM_2_A = 208   #A of target");
  sl->SetParameter("BEAM_1_GAMMA = 1470.0 #Gamma of the colliding ion 1");
  sl->SetParameter("BEAM_2_GAMMA = 1470.0 #Gamma of the colliding ion 2");
  sl->SetParameter("W_MAX = -1   #Max value of w");
  sl->SetParameter("W_MIN = -1    #Min value of w");
  sl->SetParameter("W_N_BINS = 50    #Bins i w");
  sl->SetParameter("RAP_MAX = 9.    #max y");
  sl->SetParameter("RAP_N_BINS = 200    #Bins i y");
  sl->SetParameter("CUT_PT = 0 #Cut in pT? 0 = (no, 1 = yes)");
  sl->SetParameter("PT_MIN = 1.0 #Minimum pT in GeV");
  sl->SetParameter("PT_MAX = 3.0 #Maximum pT in GeV");
  sl->SetParameter("CUT_ETA = 0 #Cut in pseudorapidity? (0 = no, 1 = yes)");
  sl->SetParameter("ETA_MIN = -10 #Minimum pseudorapidity");
  sl->SetParameter("ETA_MAX = 10 #Maximum pseudorapidity");
  sl->SetParameter("PROD_MODE = 2     #gg or gP switch (1 = 2-photon, 2 = coherent vector meson (narrow), 3 = coherent vector meson (wide), 4 = incoherent vector meson)");
  sl->SetParameter("N_EVENTS = 1000   #Number of events");
  sl->SetParameter("PROD_PID = 443013   #Channel of interest; this is j/psi --> mu+ mu-");
  sl->SetParameter("RND_SEED = 5574533 #Random number seed");
  sl->SetParameter("BREAKUP_MODE = 5     #Controls the nuclear breakup; a 5 here makes no requirement on the breakup of the ions");
  sl->SetParameter("INTERFERENCE = 0     #Interference (0 = off, 1 = on)");
  sl->SetParameter("IF_STRENGTH = 1.    #percent of intefernce (0.0 - 0.1)");
  sl->SetParameter("INT_PT_MAX = 0.24  #Maximum pt considered, when interference is turned on");
  sl->SetParameter("INT_PT_N_BINS =120   #Number of pt bins when interference is turned on");
  sl->SetParameter("XSEC_METHOD = 1 # Set to 0 to use old method for calculating gamma-gamma luminosity");
  sl->SetParameter("PYTHIA_FULL_EVENTRECORD = 0 # Write full pythia information to output (vertex, parents, daughter etc).");

  sl->InitStarLight();
  sl->PrintInputs(std::cout);
  TClonesArray tca("TParticle", 100);

  TLorentzVector v[2], vSum;
  TH1* hM  = new TH1D("hM",  "STARLIGHT;M#(){#pi^{+}#pi^{-}}",    200, 3.0, 3.2);
  TH1* hPt = new TH1D("hPt", "STARLIGHT;P_{T}#(){#pi^{+}#pi^{-}}", 80, 0., 2.);
  TH1* hY  = new TH1D("hY",  "STARLIGHT;Y#(){#pi^{+}#pi^{-}}",    100,-10., 10.);

  std::ofstream ofs("sl.txt");
  TParticle *p;
  for (Int_t counter(0); counter<20000; ) {
    sl->GenerateEvent();
    sl->BoostEvent();
    sl->ImportParticles(&tca, "ALL");
    Bool_t genOK = kTRUE;
    TLorentzVector vSum;
    for (Int_t i=0; i<tca.GetEntries() && genOK; ++i) {
      p = (TParticle*)tca.At(i);
      p->Momentum(v[i]);
      vSum += v[i];
//       genOK = TMath::Abs(v[i].Rapidity()) <= 1.5;
    }
    tca.Clear();
    if (!genOK) continue;
    Printf("%5d %d", counter, genOK);
    ++counter;
    vSum = v[0] + v[1];
    ofs << std::fixed << std::setprecision(4)
	<< vSum.M() << " " << vSum.Perp() << " " << vSum.Rapidity() << " "
	<< v[0].Eta() << " " << v[0].Px() << " " << v[0].Py() << " " << v[0].Pz() << " "
	<< v[1].Eta() << " " << v[1].Px() << " " << v[1].Py() << " " << v[1].Pz()
	<< std::endl;
    hM->Fill(vSum.M());
    hPt->Fill(vSum.Perp());
    hY->Fill(vSum.Rapidity());
  }
  TFile::Open("sl.root", "RECREATE");
  sl->Write();
  gFile->Write();

  hM->Draw();
  c1->SaveAs("SL.pdf(");
  hPt->Draw();
  c1->SaveAs("SL.pdf");
  hY->Draw();
  c1->SaveAs("SL.pdf)");
}