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

//________________________________________________________________________________
void StarMCHits::FinishEvent() {
  static const Double_t pEMax = 1 - 1.e-10;
  TDataSet *m_DataSet = StarMCHits::instance()->GetHitHolder();
  if (! m_DataSet) return;
  St_g2t_event *g2t_event = new St_g2t_event("g2t_event",1);  
  m_DataSet->Add(g2t_event);
  g2t_event_st event;
  memset (&event, 0, sizeof(g2t_event_st));
  fEventNumber++;
  event.n_event            = fEventNumber;//IHEAD(2)
  event.ge_rndm[0]         =        fSeed;//IHEAD(3)
  event.ge_rndm[1]         =            0;//IHEAD(4)
  event.n_run              =            1;
  event.n_track_eg_fs      = StarVMCApplication::Instance()->GetStack()->GetNtrack();
  event.n_track_prim       = StarVMCApplication::Instance()->GetStack()->GetNprimary();
  event.prim_vertex_p      =            1;
  event.b_impact           =           99;
  event.phi_impact         =          0.5;
  g2t_event->AddAt(&event);
  Int_t NoVertex = 1;
  St_g2t_vertex  *g2t_vertex  = new St_g2t_vertex("g2t_vertex",NoVertex);
  m_DataSet->Add(g2t_vertex); 
  g2t_vertex_st vertex;
  Int_t NTracks = StarVMCApplication::Instance()->GetStack()->GetNtrack();
  St_g2t_track   *g2t_track   = new St_g2t_track ("g2t_track",NTracks);
  m_DataSet->Add(g2t_track);
  g2t_track_st track;
  StarMCParticle  *particle = 0;   
  Int_t iv = 0;
  TLorentzVector oldV(0,0,0,0);
  TLorentzVector newV(0,0,0,0);
  TLorentzVector devV(0,0,0,0);
  for (Int_t it = 0; it <NTracks; it++) {
    memset(&track, 0, sizeof(g2t_track_st));
    particle = (StarMCParticle*) StarVMCApplication::Instance()->GetStack()->GetParticle(it);
    TParticle  *part = (TParticle *) particle->GetParticle();
    part->ProductionVertex(newV);
    devV = newV - oldV;
    if (iv == 0 || devV.Mag() > 1.e-7) {
      if (iv > 0) g2t_vertex->AddAt(&vertex);
      memset (&vertex, 0, sizeof(g2t_vertex_st));
      iv++;
      vertex.id           = iv             ;// primary key 
      vertex.event_p      = 0              ;// pointer to event
      vertex.eg_label     = 0              ;// generator label (0 if GEANT)
      vertex.eg_tof       = 0              ;// vertex production time
      vertex.eg_proc      = 0              ;// event generator mechanism
      memcpy(vertex.ge_volume,"   ",4);    ;// GEANT volume name
      vertex.ge_medium    = 0              ;// GEANT Medium
      vertex.ge_tof       = 0              ;// GEANT vertex production time
      vertex.ge_proc      = 0              ;// GEANT mechanism (0 if eg)
      vertex.ge_x[0]      = newV.X()       ;// GEANT vertex coordinate
      vertex.ge_x[1]      = newV.Y()       ;
      vertex.ge_x[2]      = newV.Z()       ;
      vertex.ge_tof       = newV.T()       ;
      vertex.n_parent     = 0              ;// number of parent tracks
      vertex.parent_p     = 0              ;// first parent track
      vertex.is_itrmd     = 0              ;// flags intermediate vertex
      vertex.next_itrmd_p = 0              ;// next intermedate vertex 
      vertex.next_prim_v_p= 0              ;// next primary vertex
      oldV                = newV;
    }
    vertex.n_daughter++;
    track.id             = it+1;
    track.eg_label       = particle->GetIdGen();
    track.eg_pid         = part->GetPdgCode();
    track.ge_pid         = gMC->IdFromPDG(track.eg_pid);
    track.start_vertex_p = iv;
    track.p[0]           = part->Px();
    track.p[1]           = part->Py();
    track.p[2]           = part->Pz();
    track.ptot           = part->P();
    track.e              = part->Energy();
    track.charge         = part->GetPDG()->Charge()/3;
    Double_t   ratio     = part->Pz()/part->Energy();
    ratio                = TMath::Min(1.-1e-10,TMath::Max(-1.+1e-10, ratio));
    track.rapidity       = TMath::ATanH(ratio);
    track.pt             = part->Pt();
    ratio                = part->Pz()/part->P();
    ratio                = TMath::Min(pEMax,TMath::Max(-pEMax, ratio));
    track.eta            = TMath::ATanH(ratio);
    g2t_track->AddAt(&track);
  }
  g2t_vertex->AddAt(&vertex);   
}

void createGlauberTree(Int_t nEvents,
                       const char *outFileName) 
{
  AliPDG::AddParticlesToPdgDataBase();
  TDatabasePDG::Instance();

  // Run loader
  TFolder *folder = new TFolder("myfolder","myfolder");
  AliRunLoader* rl = new AliRunLoader(folder);
  rl->MakeHeader();
  rl->MakeStack();
  AliStack* stack = rl->Stack();
  //AliHeader* rheader = rl->GetHeader();

  AliGenHijing *genHi = new AliGenHijing(-1);
  genHi->SetStack(stack);
  genHi->SetEnergyCMS(2760);
  genHi->SetReferenceFrame("CMS");
  genHi->SetProjectile("A", 208, 82);
  genHi->SetTarget    ("A", 208, 82);
  genHi->SetPtHardMin (2.3);
  genHi->SetImpactParameterRange(0.,30);
  genHi->SetJetQuenching(0); // enable jet quenching
  genHi->SetShadowing(1);    // enable shadowing
  genHi->SetDecaysOff(1);    // neutral pion and heavy particle decays switched off
  genHi->Init();

  MyHeader  *myheader = new MyHeader;
  MyResponse *myresp  = new MyResponse;

  TFile *outFile = TFile::Open(outFileName, "RECREATE");
  outFile->SetCompressionLevel(5);
  TDirectory::TContext context(outFile);

  TTree *tree = new TTree("glaubertree", "Glauber tree");
  tree->Branch("header",&myheader, 32*1024, 99);
  tree->Branch("response",&myresp, 32*1024, 99);

  TNtuple *ntuple = new TNtuple("gnt", "Glauber ntuple", "npart:ncoll:b");

  Double_t etas[] = {-10,-5,-4,-3,-2,-1,0,1,2,3,4,5,10};
  TH1D *hNEta = new TH1D("hNeta","",12,etas);
  TH1D *hEtEta = new TH1D("hEteta","",12,etas);

  // create events and fill them
  for (Int_t iEvent = 0; iEvent < nEvents; ++iEvent) {

    cout << "Event " << iEvent+1 << "/" << nEvents << endl;;
    stack->Reset();
    hNEta->Reset();
    hEtEta->Reset();
    genHi->Generate();
  
    AliStack *s = genHi->GetStack();
    const TObjArray *parts = s->Particles();
    Int_t nents = parts->GetEntries();
    for (Int_t i = 0; i<nents; ++i) {
      TParticle *p = (TParticle*)parts->At(i);
      //p->Print();
      TParticlePDG *pdg = p->GetPDG(1);
      Int_t c = (Int_t)(TMath::Abs(pdg->Charge()));
      if (c!=0) {
        hNEta->Fill(p->Eta());
        hEtEta->Fill(p->Eta(),p->Pt());
      }
    }

    AliGenHijingEventHeader *h = (AliGenHijingEventHeader*)genHi->CollisionGeometry();
    myheader->fNATT = nents;
    myheader->fEATT = h->TotalEnergy();
    myheader->fJATT = h->HardScatters();
    myheader->fNT   = h->TargetParticipants();
    myheader->fNP   = h->ProjectileParticipants();
    myheader->fN00  = h->NwNw();
    myheader->fN01  = h->NwN();
    myheader->fN10  = h->NNw();
    myheader->fN11  = h->NN();
    myheader->fBB   = h->ImpactParameter();
    myheader->fRP   = h->ReactionPlaneAngle();
    myheader->fPSn  = h->ProjSpectatorsn();
    myheader->fPSp  = h->ProjSpectatorsp();
    myheader->fTSn  = h->TargSpectatorsn();
    myheader->fTSp  = h->TargSpectatorsn();

    myresp->fEtch0p = hEtEta->GetBinContent(hEtEta->FindBin(0.5));
    myresp->fEtch1p = hEtEta->GetBinContent(hEtEta->FindBin(1.5));
    myresp->fEtch2p = hEtEta->GetBinContent(hEtEta->FindBin(2.5));
    myresp->fEtch3p = hEtEta->GetBinContent(hEtEta->FindBin(3.5));
    myresp->fEtch4p = hEtEta->GetBinContent(hEtEta->FindBin(4.5));
    myresp->fEtch5p = hEtEta->GetBinContent(hEtEta->FindBin(5.5));
    myresp->fEtchrp = hEtEta->GetBinContent(hEtEta->FindBin(10.5));
    myresp->fEtch0n = hEtEta->GetBinContent(hEtEta->FindBin(-0.5));
    myresp->fEtch1n = hEtEta->GetBinContent(hEtEta->FindBin(-1.5));
    myresp->fEtch2n = hEtEta->GetBinContent(hEtEta->FindBin(-2.5));
    myresp->fEtch3n = hEtEta->GetBinContent(hEtEta->FindBin(-3.5));
    myresp->fEtch4n = hEtEta->GetBinContent(hEtEta->FindBin(-4.5));
    myresp->fEtch5n = hEtEta->GetBinContent(hEtEta->FindBin(-5.5));
    myresp->fEtchrn = hEtEta->GetBinContent(hEtEta->FindBin(-10.5));
    myresp->fNch0p  = hNEta->GetBinContent(hNEta->FindBin(0.5));
    myresp->fNch1p  = hNEta->GetBinContent(hNEta->FindBin(1.5));
//.........这里部分代码省略.........