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

//_____________________________________________________________________________
Bool_t ProofPythia::Process(Long64_t entry)
{
   // Main event loop  

   fPythia->GenerateEvent();
   if (entry < 2)
      fPythia->EventListing();
   fPythia->ImportParticles(fP, "All");
   Int_t nTot = fPythia->GetN();
   fPythia->ImportParticles(fP, "All");
   Int_t np = fP->GetEntriesFast();
   // Particle loop
   Int_t nCharged = 0;
   for (Int_t ip = 0; ip < np; ip++) {
      TParticle* part = (TParticle*) fP->At(ip);
      Int_t ist = part->GetStatusCode();
      Int_t pdg = part->GetPdgCode();
      if (ist != 1) continue;
      Float_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
      if (charge == 0.) continue;
      nCharged++;
      Float_t eta = part->Eta();
      Float_t pt  = part->Pt();
      if (pt > 0.) fPt->Fill(pt);
      if ((eta > -10) && (eta < 10)) fEta->Fill(eta);
   }
   fHist->Fill(nCharged);
   fTot->Fill(nTot);

   return kTRUE;
}

void pythia8(Int_t nev  = 100, Int_t ndeb = 1)
{
   const char *p8dataenv = gSystem->Getenv("PYTHIA8DATA");
   if (!p8dataenv) {
      const char *p8env = gSystem->Getenv("PYTHIA8");
      if (!p8env) {
         Error("pythia8.C",
               "Environment variable PYTHIA8 must contain path to pythia directory!");
         return;
      }
      TString p8d = p8env;
      p8d += "/xmldoc";
      gSystem->Setenv("PYTHIA8DATA", p8d);
   }

   const char* path = gSystem->ExpandPathName("$PYTHIA8DATA");
   if (gSystem->AccessPathName(path)) {
         Error("pythia8.C",
               "Environment variable PYTHIA8DATA must contain path to $PYTHIA8/xmldoc directory !");
      return;
   }

// Load libraries
#ifndef G__WIN32 // Pythia8 is a static library on Windows
   if (gSystem->Getenv("PYTHIA8")) {
      gSystem->Load("$PYTHIA8/lib/libpythia8");
   } else {
      gSystem->Load("libpythia8");
   }
#endif
   gSystem->Load("libEG");
   gSystem->Load("libEGPythia8");
// Histograms
   TH1F* etaH = new TH1F("etaH", "Pseudorapidity", 120, -12., 12.);
   TH1F* ptH  = new TH1F("ptH",  "pt",              50,   0., 10.);


// Array of particles
   TClonesArray* particles = new TClonesArray("TParticle", 1000);
// Create pythia8 object
   TPythia8* pythia8 = new TPythia8();

// Configure
   pythia8->ReadString("HardQCD:all = on");


// Initialize

   pythia8->Initialize(2212 /* p */, 2212 /* p */, 14000. /* TeV */);

// Event loop
   for (Int_t iev = 0; iev < nev; iev++) {
      pythia8->GenerateEvent();
      if (iev < ndeb) pythia8->EventListing();
      pythia8->ImportParticles(particles,"All");
      Int_t np = particles->GetEntriesFast();
// Particle loop
      for (Int_t ip = 0; ip < np; ip++) {
         TParticle* part = (TParticle*) particles->At(ip);
         Int_t ist = part->GetStatusCode();
         // Positive codes are final particles.
         if (ist <= 0) continue;
         Int_t pdg = part->GetPdgCode();
         Float_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
         if (charge == 0.) continue;
         Float_t eta = part->Eta();
         Float_t pt  = part->Pt();

         etaH->Fill(eta);
         if (pt > 0.) ptH->Fill(pt, 1./(2. * pt));
      }
   }

   pythia8->PrintStatistics();

   TCanvas* c1 = new TCanvas("c1","Pythia8 test example",800,800);
   c1->Divide(1, 2);
   c1->cd(1);
   etaH->Scale(5./Float_t(nev));
   etaH->Draw();
   etaH->SetXTitle("#eta");
   etaH->SetYTitle("dN/d#eta");

   c1->cd(2);
   gPad->SetLogy();
   ptH->Scale(5./Float_t(nev));
   ptH->Draw();
   ptH->SetXTitle("p_{t} [GeV/c]");
   ptH->SetYTitle("dN/dp_{t}^{2} [GeV/c]^{-2}");
 }

void pythia8_susy() {

  Int_t maxEvts = 100; // Maximo numero de eventos

  char* path = gSystem->ExpandPathName("$PYTHIA8DATA");
  if (gSystem->AccessPathName(path)) {
    Warning("pythia8.C", 
            "Environment variable PYTHIA8DATA must contain path to pythi8100/xmldoc directory !");
    return;
  }
  
  // Load libraries
  gSystem->Load("$PYTHIA8/lib/libpythia8");
  gSystem->Load("$PYTHIA8/lib/liblhapdfdummy");
  
  gSystem->Load("libEG");
  gSystem->Load("libEGPythia8");
  
  //Definir archivo de salida
  TFile * outfile = new TFile("eventos_pythia8_SUSY.root","RECREATE");
  
  // Array of particles
  TClonesArray* particles = new TClonesArray("TParticle", 5000);
  
  //Definir el TTree
  TTree*tree= new TTree("tree","Arbol con particulas segun Pythia8");
  tree->Branch("particles",&particles);
  
  // Create pythia8 object
  TPythia8* pythia8 = new TPythia8();
  
  //*Configurar: Aqui seleccione el proceso que quiere simular    
  pythia8->ReadString("SUSY:all = on"); //Todos los procesos susy posibles
  //pythia8->ReadString("SUSY:qqbar2chi+-chi0 = on"); //Un proceso en especial
  
  //Importante: pasar a Pythia8 el nombre del archivo SLHA
  pythia8->ReadString("SLHA:file = SUSY_LM2_sftsht.slha"); //insertar aqui el nombre del archivo SLHA

  // Initialize 
  
  pythia8->Initialize(2212 /* p */, 2212 /* p */, 7000. /* TeV */);
  
  int iev = 0;
  
  // Event loop
  
  while( iev < maxEvts ) {
    
    pythia8->GenerateEvent();
    if (iev < 1) pythia8->EventListing();
    pythia8->ImportParticles(particles,"All");
    
    Int_t np = particles->GetEntriesFast();
    
    // Particle loop
    
    for (Int_t ip = 0; ip < np; ip++) {
      
      TParticle* part = (TParticle*) particles->At(ip);
      Int_t ist = part->GetStatusCode();
      Int_t pdg = part->GetPdgCode();
      
    }

   tree->Fill();
   ++iev;
     	 
  }
  
  pythia8->PrintStatistics();
  
  outfile->Write();
  outfile->Close();
  
   
}

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");

}

void Pythia8(const Int_t nEvents = 10)
{
  gROOT->LoadMacro("TUtils.h");

  if (LoadRootLibs()) return;
  if (LoadPythia8())  return;
  if (LoadThermalClass()) return;
//=============================================================================

  TPythia8 *pythia8 = new TPythia8();

  pythia8->ReadString("SoftQCD:all = on");
  pythia8->ReadString("SoftQCD:singleDiffractive = on");
  pythia8->ReadString("SoftQCD:doubleDiffractive = on");

  pythia8->Initialize(2212, 2212, 14000.);
//=============================================================================

  TGenThermalParticles *thermal = new TGenThermalParticles("Boltzmann");

  thermal->SetMultiplicity(2000);

  thermal->SetMeanPt(0.7);
  thermal->SetPtRange(0.15, 200.);

  thermal->SetEtaRange(-0.8, 0.8);
  thermal->SetPhiRange(0., TMath::TwoPi());
//=============================================================================

  TClonesArray *particles = new TClonesArray("TParticle", 1000);

  for (Int_t iEvent=0; iEvent<nEvents; iEvent++) {
    pythia8->GenerateEvent();
    if (iEvent==0) pythia8->EventListing();

    pythia8->ImportParticles(particles, "Final");
    Int_t nb = particles->GetEntriesFast();
    cout << "iEvent = "<< iEvent << ", np before = " << nb;

    thermal->ImportParticles(particles, "Boltzmann");
    Int_t na = particles->GetEntriesFast();
    cout << ", np after = " << na << endl;

    TParticle *part = 0;
    for (Int_t i=0; i<na; i++) {
      part = (TParticle*)particles->At(i); if (!part) continue;
      Bool_t bThermalBkg = (part->GetStatusCode()==-1);

      if (!bThermalBkg) {
        Int_t kPDG = part->GetPdgCode();
        Float_t dCharge = TDatabasePDG::Instance()->GetParticle(kPDG)->Charge();
      }

      part = 0;
    }
  }
//=============================================================================

   pythia8->PrintStatistics();
//=============================================================================

  return;
}

TEveTrackList*
kine_tracks(Double_t min_pt,  Double_t min_p,
	    Bool_t   pdg_col, Bool_t   recurse,
	    Bool_t   use_track_refs)
{
  IlcRunLoader* rl =  IlcEveEventManager::AssertRunLoader();
  rl->LoadKinematics();
  IlcStack* stack = rl->Stack();
  if (!stack)
  {
    Error("kine_tracks", "can not get kinematics.");
    return 0;
  }

  gEve->DisableRedraw();

  TEveTrackList* cont = new TEveTrackList("Kine Tracks");
  cont->SetMainColor(3);
  TEveTrackPropagator* trkProp = cont->GetPropagator();

  kine_track_propagator_setup(trkProp);

  gEve->AddElement(cont);
  Int_t count = 0;
  Int_t Np = stack->GetNprimary();
  for (Int_t i = 0; i < Np; ++i)
  {
    TParticle* p = stack->Particle(i);
    if (p->GetStatusCode() <= 1)
    {
      if (p->Pt() < min_pt && p->P() < min_p) continue;

      ++count;
      IlcEveTrack* track = new IlcEveTrack(p, i, trkProp);

      //PH The line below is replaced waiting for a fix in Root
      //PH which permits to use variable siza arguments in CINT
      //PH on some platforms (alphalinuxgcc, solariscc5, etc.)
      //PH    track->SetName(Form("%s [%d]", p->GetName(), i));
      char form[1000];
      sprintf(form,"%s [%d]", p->GetName(), i);
      track->SetName(form);
      track->SetStdTitle();
      Int_t ml = p->GetMother(0);
      if (ml != -1)
      {
        track->SetTitle(Form("%s\nMother label=%d\nMother Pdg=%d",
                             track->GetElementTitle(),
                             ml, stack->Particle(ml)->GetPdgCode()));
      }
      set_track_color(track, pdg_col);

      gEve->AddElement(track, cont);

      if (recurse)
	kine_daughters(track, stack, min_pt, min_p, pdg_col, recurse);
    }
  }

  // set path marks
  IlcEveKineTools kt;
  kt.SetDaughterPathMarks(cont, stack, recurse);
  if (use_track_refs && rl->LoadTrackRefs() == 0)
  {
    kt.SetTrackReferences(cont, rl->TreeTR(), recurse);
    trkProp->SetEditPathMarks(kTRUE);
  }
  kt.SortPathMarks(cont, recurse);

  //PH  const Text_t* tooltip = Form("min pT=%.2lf, min P=%.2lf), N=%d", min_pt, min_p, count);
  char tooltip[1000];
  sprintf(tooltip,"min pT=%.2lf, min P=%.2lf), N=%d", min_pt, min_p, count);
  cont->SetTitle(tooltip); // Not broadcasted automatically ...

  cont->MakeTracks(recurse);
  gEve->EnableRedraw();
  gEve->Redraw3D();

  return cont;
}