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

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

//.........这里部分代码省略.........
 plots.analyzed_mutau = 0;
  
 plots.passedJetAndLepNumberSelections_ee = 0;
 plots.passedJetAndLepNumberSelections_mumu = 0;
 plots.passedJetAndLepNumberSelections_tautau = 0;
 plots.passedJetAndLepNumberSelections_emu = 0;
 plots.passedJetAndLepNumberSelections_etau = 0;
 plots.passedJetAndLepNumberSelections_mutau = 0;

 
 
 //PG loop over the events
 for (int evt = 0 ; evt < nentries ; ++evt)
 {

  tree->GetEntry (evt) ;
  
  tagJets -> Clear () ;  
  otherJets -> Clear () ;    
  
  
    //---- check if signal ----
  if (if_signal && (IdEvent!=123 && IdEvent!=124)) continue;
  plots.analyzed++;
  
   //!---- MC ----

  if (IdEvent==123 || IdEvent==124) { //---- VBF event ----
   plots.v_decay_Channel_e = 0;
   plots.v_decay_Channel_mu = 0;
   plots.v_decay_Channel_tau = 0;
   for (int iGen = 0; iGen < genParticles->GetEntries() ; ++iGen){
    TParticle* myparticle = (TParticle*) genParticles->At(iGen);
    if (abs(myparticle->GetPdgCode()) == 24) { //---- W
     Int_t mother1 = 0;
     mother1 = myparticle->GetMother(0);
     if (mother1 == 25) { //---- mother is higgs ----
      for (int iDaughter = 0; iDaughter<2; iDaughter++){
       if (abs(myparticle->GetDaughter(iDaughter)) == 11) {//---- W -> e
        plots.v_decay_Channel_e++;
       }
       if (abs(myparticle->GetDaughter(iDaughter)) == 13) {//---- W -> mu    
        plots.v_decay_Channel_mu++;
       }
       if (abs(myparticle->GetDaughter(iDaughter)) == 15) {//---- W -> tau
        plots.v_decay_Channel_tau++;
       }       
      }
     }
    }
   }
  }

  if (plots.v_decay_Channel_e == 2) plots.analyzed_ee++;
  if (plots.v_decay_Channel_mu == 2) plots.analyzed_mumu++;
  if (plots.v_decay_Channel_tau == 2) plots.analyzed_tautau++;
  if (plots.v_decay_Channel_e == 1 && plots.v_decay_Channel_mu == 1) plots.analyzed_emu++;
  if (plots.v_decay_Channel_e == 1 && plots.v_decay_Channel_tau == 1) plots.analyzed_etau++;
  if (plots.v_decay_Channel_mu == 1 && plots.v_decay_Channel_tau == 1) plots.analyzed_mutau++;
  
  
  
  
  
  
  

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

//.........这里部分代码省略.........
				     100, 1.0, 1.2, 
				     "M(#bar{p}#pi^{+}) [GeV/c^{2}]", "N");
  Int_t nGenV0s = 0;
  Int_t nRecV0s = 0;
  TH1F* hMassXi = CreateHisto("hMassXi", "#Xi", 100, 1.2, 1.5, 
			      "M(#Lambda#pi) [GeV/c^{2}]", "N");
  TH1F* hMassOmega = CreateHisto("hMassOmega", "#Omega", 100, 1.5, 1.8, 
				 "M(#LambdaK) [GeV/c^{2}]", "N");
  Int_t nGenCascades = 0;
  Int_t nRecCascades = 0;

  // loop over events
  for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
    runLoader->GetEvent(iEvent);

    // select simulated primary particles, V0s and cascades
    AliStack* stack = runLoader->Stack();
    Int_t nParticles = stack->GetNtrack();
    TArrayF vertex(3);
    runLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
    TObjArray selParticles;
    TObjArray selV0s;
    TObjArray selCascades;
    for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
      TParticle* particle = stack->Particle(iParticle);
      if (!particle) continue;
      if (particle->Pt() < 0.001) continue;
      if (TMath::Abs(particle->Eta()) > 0.9) continue;
      TVector3 dVertex(particle->Vx() - vertex[0], 
		       particle->Vy() - vertex[1],
		       particle->Vz() - vertex[2]);
      if (dVertex.Mag() > 0.0001) continue;

      switch (TMath::Abs(particle->GetPdgCode())) {
      case kElectron:
      case kMuonMinus:
      case kPiPlus:
      case kKPlus:
      case kProton: {
	if (particle->Pt() > minPt) {
	  selParticles.Add(particle);
	  nGen++;
	  hGen->Fill(particle->Pt());
	}
	break;
      }
      case kK0Short:
      case kLambda0: {
	if (particle->Pt() > cutPtV0) {
	  nGenV0s++;
	  selV0s.Add(particle);
	}
	break;
      }
      case kXiMinus:
      case kOmegaMinus: {
	if (particle->Pt() > cutPtCascade) {
	  nGenCascades++;
	  selCascades.Add(particle);
	}
	break;
      }
      default: break;
      }
    }

//.........这里部分代码省略.........
    tree->GetEvent(iEv);
    runLoader->GetEvent(iEv);
    
    printf("+++ event %i (of %lld) +++++++++++++++++++++++  # ESDtracks: %d \n",iEv,tree->GetEntries()-1,esd->GetNumberOfTracks());
    Int_t nESDtracks = esd->GetNumberOfTracks();
    for (Int_t iTrack = 0; iTrack < nESDtracks; iTrack++) {
      AliESDtrack* track = esd->GetTrack(iTrack);
      if (!(iTrack%1000)) printf("event %i: ESD track count %d (of %d)\n",iEv,iTrack,nESDtracks);

      Int_t label = track->GetLabel();
  
      Int_t idx[12];
      //      Int_t ncl = track->GetITSclusters(idx);
   
      if(label<0) {
	//	cout<< " ESD track label " << label;
	//	cout<<"  ---> imperfect track (label "<<label<<"<0) !! -> track Pt: "<< track->Pt() << endl;
      }

      AliStack* stack = runLoader->Stack();
      //     nTrackTotalMC += stack->GetNprimary();
    

      TParticle* particle = stack->Particle(TMath::Abs(label)); 
      Double_t pt = track->Pt();
      
      if(particle) {

	if (TMath::Abs(particle->Eta())>etaCut) continue;

	Double_t ptMC = particle->Pt();

	// Efficiencies
	if (onlyPion && TMath::Abs(particle->GetPdgCode())!=211) continue;

	if ( (!onlyPrims) || stack->IsPhysicalPrimary(TMath::Abs(label))) {
	  //  cout<<" # clusters "<<ncl<<endl;

	  nTrackFound++;
	  hAllFound->Fill(ptMC);
	  hEta->Fill(track->Eta());
	  
	  if (label<0) {
	    nTrackImperfect++;
	    hImperfect->Fill(ptMC);
	  } else {
	    nTrackPerfect++;
	    hPerfect->Fill(ptMC);
	  }

	}


	// following only for "true tracks, pions

	if(particle->Pt() < 0.001)continue;
	if (TMath::Abs(particle->GetPdgCode())!=211) continue;
	if (label>0) {
	  
	  // Impact parameters for Pions only
	  Double_t dca = track->GetD(0,0,0.5);
	  h2Ddca->Fill(ptMC,dca);
	  
	  // Pt resolution for Pions only
	  Double_t dPt = (pt-ptMC)/ptMC*100;
	  h2Dpt->Fill(ptMC,dPt);

void CountTrackableMCs(TH1F *hAllMC, Bool_t onlyPrims,Bool_t onlyPion) {
  
  gSystem->Load("libITSUpgradeBase");
  gSystem->Load("libITSUpgradeSim");

 // open run loader and load gAlice, kinematics and header
  AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
  if (!runLoader) {
    Error("Check kine", "getting run loader from file %s failed",
          "galice.root");
    return;
  }
  runLoader->LoadHeader();
  runLoader->LoadKinematics();
  runLoader->LoadTrackRefs();

  AliLoader *dl = runLoader->GetDetectorLoader("ITS");

  //Trackf
  TTree *trackRefTree = 0x0; 
  TClonesArray *trackRef = new TClonesArray("AliTrackReference",1000);

  //  TH1F *hRef = new TH1F("","",100,0,100);
  TH1F *hR = new TH1F("","",100,0,100);
  if (hAllMC==0) hAllMC = new TH1F("","",100,0.1,2);
  Float_t ptmin = hAllMC->GetBinCenter(1)-hAllMC->GetBinWidth(1)/2;
  Float_t ptmax = hAllMC->GetBinCenter(hAllMC->GetNbinsX())+hAllMC->GetBinWidth(hAllMC->GetNbinsX())/2;
  //  Int_t nAllMC = 0;

  // Detector geometry
  TArrayD rmin(0);   TArrayD rmax(0); 
  GetDetectorRadii(&rmin,&rmax);
  TArrayI nLaySigs(rmin.GetSize());

  printf("Counting trackable MC tracks ...\n");
  
  for(Int_t iEv =0; iEv<runLoader->GetNumberOfEvents(); iEv++){
    Int_t nTrackableTracks = 0;
    runLoader->GetEvent(iEv);
    AliStack* stack = runLoader->Stack();  
    printf("+++ event %i (of %d) +++++++++++++++++++++++  # total MCtracks: %d \n",iEv,runLoader->GetNumberOfEvents()-1,stack->GetNtrack());

    trackRefTree=runLoader->TreeTR();
    TBranch *br = trackRefTree->GetBranch("TrackReferences");
    if(!br) {
      printf("no TR branch available , exiting \n");
      return;
    }
    br->SetAddress(&trackRef);

    // init the trackRef tree 
    trackRefTree=runLoader->TreeTR();
    trackRefTree->SetBranchAddress("TrackReferences",&trackRef);
 
    // Count trackable MC tracks
    for (Int_t iMC=0; iMC<stack->GetNtrack(); iMC++) {

      TParticle* particle = stack->Particle(iMC); 
      if (TMath::Abs(particle->Eta())>etaCut) continue;
      if (onlyPrims && !stack->IsPhysicalPrimary(iMC)) continue;
      if (onlyPion && TMath::Abs(particle->GetPdgCode())!=211) continue;


      Bool_t isTrackable = 0;
      nLaySigs.Reset(0);
 
      trackRefTree->GetEntry(stack->TreeKEntry(iMC));
      Int_t nref=trackRef->GetEntriesFast();
      for(Int_t iref =0; iref<nref; iref++){
	AliTrackReference *trR = (AliTrackReference*)trackRef->At(iref);
	if(!trR) continue;
	if(trR->DetectorId()!=AliTrackReference::kITS) continue;
	Float_t radPos = trR->R();
	hR->Fill(radPos);
	for (Int_t il=0; il<rmin.GetSize();il++) {
	  if (radPos>=rmin.At(il)-0.1 && radPos<=rmax.At(il)+0.1) {
	    //	    cout<<"  in Layer "<<il<<" "<<radPos;
	    nLaySigs.AddAt(1.,il);
	    //	    cout<<" "<<nLaySigs.At(il)<<endl;
	  }
	}
      }

      if (nLaySigs.GetSum()>=3) {
	isTrackable =1;
	//	cout<<nLaySigs.GetSum()<<endl;
      }
      
      if (isTrackable) {
	Double_t ptMC = particle->Pt();
	//	Double_t etaMC = particle->Eta();
	//	if (ptMC>ptmin&&ptMC<ptmax) {nTrackableTracks++;hAllMC->Fill(ptMC);}
	if (ptMC>ptmin) {nTrackableTracks++;hAllMC->Fill(ptMC);}

      }

      
    } // entries tracks MC
    printf(" -> trackable MC tracks: %d (%d)\n",nTrackableTracks,hAllMC->GetEntries());
  }//entries Events
//.........这里部分代码省略.........

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

void fastGen(Int_t nev = 1, char* filename = "gilc.root")
{
  IlcPDG::AddParticlesToPdgDataBase();
  TDatabasePDG::Instance();
 


  // Run loader
  IlcRunLoader* rl = IlcRunLoader::Open("gilc.root","FASTRUN","recreate");
  
  rl->SetCompressionLevel(2);
  rl->SetNumberOfEventsPerFile(nev);
  rl->LoadKinematics("RECREATE");
  rl->MakeTree("E");
  gIlc->SetRunLoader(rl);
  
  //  Create stack
  rl->MakeStack();
  IlcStack* stack = rl->Stack();
  
  //  Header
  IlcHeader* header = rl->GetHeader();
  
  //  Create and Initialize Generator
 
  // Example of charm generation taken from Config_PythiaHeavyFlavours.C
  IlcGenPythia *gener = new IlcGenPythia(-1);
  gener->SetEnergyCMS(14000.);
  gener->SetMomentumRange(0,999999);
  gener->SetPhiRange(0., 360.);
  gener->SetThetaRange(0.,180.);
  //  gener->SetProcess(kPyCharmppMNR); // Correct Pt distribution, wrong mult
  gener->SetProcess(kPyMb); // Correct multiplicity, wrong Pt
  gener->SetStrucFunc(kCTEQ4L);
  gener->SetPtHard(2.1,-1.0);
  gener->SetFeedDownHigherFamily(kFALSE);
  gener->SetStack(stack);
  gener->Init();

  // Go to gilc.root
  rl->CdGAFile();

  // Forbid some decays. Do it after gener->Init(0, because
  // the initialization of the generator includes reading of the decay table.

  IlcPythia * py= IlcPythia::Instance();
  py->SetMDME(737,1,0); //forbid D*+->D+ + pi0
  py->SetMDME(738,1,0);//forbid D*+->D+ + gamma

  // Forbid all D0 decays except D0->K- pi+
  for(Int_t d=747; d<=762; d++){ 
    py->SetMDME(d,1,0);
  }
  // decay 763 is D0->K- pi+
  for(Int_t d=764; d<=807; d++){ 
    py->SetMDME(d,1,0);
  }

  //
  //                        Event Loop
  //
  
  TStopwatch timer;
  timer.Start();
  for (Int_t iev = 0; iev < nev; iev++) {
    
    cout <<"Event number "<< iev << endl;
    
    //  Initialize event
    header->Reset(0,iev);
    rl->SetEventNumber(iev);
    stack->Reset();
    rl->MakeTree("K");
    
    //  Generate event
    Int_t nprim = 0;
    Int_t ntrial = 0;
    Int_t ndstar = 0;
   
   

    //-------------------------------------------------------------------------------------

    while(!ndstar) {
      // Selection of events with D*
      stack->Reset();
      stack->ConnectTree(rl->TreeK());
      gener->Generate();
      ntrial++;
      nprim = stack->GetNprimary();
      
      for(Int_t ipart =0; ipart < nprim; ipart++){
        TParticle * part = stack->Particle(ipart);
        if(part)    {
          
          if (TMath::Abs(part->GetPdgCode())== 413) {

	    TArrayI daughtersId;

	    GetFinalDecayProducts(ipart,*stack,daughtersId);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	  DB.AccountTopology(*cl, cSum.dX, cSum.dZ, cSum.alpha, cSum.beta);
          
          GetHistoClSize(clsize,kDR,&histoArr)->Fill((rht-rcl)*1e4);
          if (cl->TestBit(kSplit)) {
            if (col%2) GetHistoClSize(clsize,kDTXoddSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
            else       GetHistoClSize(clsize,kDTXevenSPL,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
            GetHistoClSize(clsize,kDTZSPL,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4);
            GetHistoClSize(0,kNPixSPL,&histoArr)->Fill(clsize);
          }
          if (col%2) GetHistoClSize(clsize,kDTXodd,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
          else       GetHistoClSize(clsize,kDTXeven,&histoArr)->Fill((txyzH[0]-xyzClTr[0])*1e4);
          GetHistoClSize(clsize,kDTZ,&histoArr)->Fill((txyzH[2]-xyzClTr[2])*1e4);
          GetHistoClSize(0,kNPixAll,&histoArr)->Fill(clsize);
          //
          cSum.evID = iEvent;
          cSum.volID = cl->GetVolumeId();
          cSum.lrID = ilr;
          cSum.clID = icl;
          cSum.nPix = cl->GetNPix();
          cSum.nX   = cl->GetNx();
          cSum.nZ   = cl->GetNz();
          cSum.q    = cl->GetQ();
          cSum.split = cl->TestBit(kSplit);
          cSum.dX = (txyzH[0]-xyzClTr[0])*1e4;
          cSum.dY = (txyzH[1]-xyzClTr[1])*1e4;
          cSum.dZ = (txyzH[2]-xyzClTr[2])*1e4;
          cSum.nRowPatt = cl-> GetPatternRowSpan();
          cSum.nColPatt = cl-> GetPatternColSpan();
 
                    
          int label = cl->GetLabel(0);
          TParticle* part = 0;
          if (label>=0 && (part=stack->Particle(label)) ) {
            cSum.pdg = part->GetPdgCode();
            cSum.eta = part->Eta();
            cSum.pt  = part->Pt();
            cSum.phi = part->Phi();
            cSum.prim = stack->IsPhysicalPrimary(label);
          } 
          cSum.ntr = 0;
          for (int ilb=0;ilb<3;ilb++) if (cl->GetLabel(ilb)>=0) cSum.ntr++;
          for (int i=0;i<3;i++) cSum.xyz[i] = xyzClGloF[i];
          //
          trOut->Fill();
          /*
          if (clsize==5) {
            printf("\nL%d(%c) Mod%d, Cl:%d | %+5.1f %+5.1f (%d/%d)|H:%e %e %e | C:%e %e %e\n",ilr,cl->TestBit(kSplit) ? 'S':'N',
             modID,icl,(txyzH[0]-xyzClTr[0])*1e4,(txyzH[2]-xyzClTr[2])*1e4, row,col,
             gxyzH[0],gxyzH[1],gxyzH[2],xyzClGlo[0],xyzClGlo[1],xyzClGlo[2]);
            cl->Print();
            pHit->Print();
            //
            double a0,b0,c0,a1,b1,c1,e0;
            pHit->GetPositionL0(a0,b0,c0,e0);
            pHit->GetPositionL(a1,b1,c1);
            float cloc[3];
            cl->GetLocalXYZ(cloc);
            printf("LocH: %e %e %e | %e %e %e\n",a0,b0,c0,a1,b1,c1);
            printf("LocC: %e %e %e | %e %e %e\n",cloc[0],cloc[1],cloc[2],xyzClTr[0],xyzClTr[1],xyzClTr[2]);
          }
          */
          //
        }
      }
    }
    

void fastGenPA(Int_t nev = 1, char* filename = "gilc.root")
{
//  Runloader
    IlcRunLoader* rl = IlcRunLoader::Open("gilc.root", "FASTRUN", "recreate");

    rl->SetCompressionLevel(2);
    rl->SetNumberOfEventsPerFile(10000);
    rl->LoadKinematics("RECREATE");
    rl->MakeTree("E");
    gIlc->SetRunLoader(rl);


//  Create stack
    rl->MakeStack();
    IlcStack* stack      = rl->Stack();
 
//  Header
    IlcHeader* header = rl->GetHeader();


//  Create and Initialize Generator
    IlcGenerator *gener = CreateGenerator();
    gener->Init();
    gener->SetStack(stack);
    
//
//                        Event Loop
//
    Int_t iev;
     
    for (iev = 0; iev < nev; iev++) {

	printf("\n \n Event number %d \n \n", iev);

	//  Initialize event

	header->Reset(0,iev);
	rl->SetEventNumber(iev);
	stack->Reset();
	rl->MakeTree("K");

//  Generate event

	gener->Generate();

//  Analysis
	Int_t npart = stack->GetNprimary();
	printf("Analyse %d Particles\n", npart);
	for (Int_t part=0; part<npart; part++) {
	    TParticle *MPart = stack->Particle(part);
	    Int_t mpart  = MPart->GetPdgCode();
	}
	
//  Finish event
	header->SetNprimary(stack->GetNprimary());
	header->SetNtrack(stack->GetNtrack());  
//      I/O
//	
	stack->FinishEvent();
	header->SetStack(stack);
	rl->TreeE()->Fill();
	rl->WriteKinematics("OVERWRITE");
    } // event loop
//
//                         Termination
//  Generator
    gener->FinishRun();

    //  Write file
    rl->WriteHeader("OVERWRITE");
    gener->Write();
    rl->Write();
}