void CountPrimaries(TH1F *hMultCount) {

  if (hMultCount==0) hMultCount = new TH1F("mult","averaged multiplicity (charg. prim)",80,-4.,4.);
  
  AliRunLoader *rl = AliRunLoader::Open("galice.root");
  rl->SetKineFileName("Kinematics.root");
  rl->LoadHeader();
  rl->LoadKinematics(); 
  Int_t nEvents = rl->GetNumberOfEvents();
  cout<< "N events "<<nEvents<<endl;
  for(Int_t iEv=0; iEv<nEvents; iEv++){
    rl->GetEvent(iEv);
    AliStack *s = rl->Stack();
    for(Int_t iP=0; iP<s->GetNtrack(); iP++ ){
      TParticle *p = s->Particle(iP);
      if (!(s->IsPhysicalPrimary(iP))) continue;
      Float_t eta = p->Eta();
      if (p->Pt()>0.06) {
	hMultCount->Fill(eta);
      }
    }
  }

  hMultCount->DrawCopy();
  rl->UnloadHeader();
  rl->UnloadKinematics();
  delete rl;



}

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

//_____________________________________________________________________________
Int_t  amsvmc_MCStack::GetCurrentParentTrackNumber() const 
{
/// \return  The current track parent ID.

  TParticle* current = GetCurrentTrack();

  if (current) 
    return current->GetFirstMother();
  else 
    return -1;
}  

void kine_daughters(IlcEveTrack* parent,  IlcStack* stack,
		    Double_t     min_pt,  Double_t  min_p,
		    Bool_t       pdg_col, Bool_t    recurse)
{
  TParticle *p = stack->Particle(parent->GetLabel());
  if (p->GetNDaughters() > 0)
  {
    TEveTrackPropagator* rs = parent->GetPropagator();
    for (int d=p->GetFirstDaughter(); d>0 && d<=p->GetLastDaughter(); ++d)
    {
      TParticle* dp = stack->Particle(d);
      if (dp->Pt() < min_pt && dp->P() < min_p) continue;

      IlcEveTrack* dtrack = new IlcEveTrack(dp, d, rs);
      char form[1000];
      sprintf(form,"%s [%d]", dp->GetName(), d);
      dtrack->SetName(form);
      dtrack->SetStdTitle();
      set_track_color(dtrack, pdg_col);

      gEve->AddElement(dtrack, parent);

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

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

}

//_____________________________________________________________________________
TParticle* amsvmc_MCStack::PopNextTrack(Int_t& itrack)
{
/// Get next particle for tracking from the stack.
/// \return       The popped particle object
/// \param track  The index of the popped track
//  cout<<"DEBUG: in PopNextTrack:"<<endl;
  itrack = -1;
  if  (fStack.empty()) return 0;
  TParticle* particle = fStack.top();
  fStack.pop();
  if (!particle) return 0;
  fCurrentTrack = particle->GetSecondMother();
  itrack = fCurrentTrack;
  //  cout << "fCurrentTrack:   " <<  fCurrentTrack << endl;
  //  fCurrentTrack = particle->GetID();
  //  itrack = fCurrentTrack;
  return particle;
}    

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 TVolumeSource::FindPotentialMinimum(){
	cout << "Sampling phase space ";
	const int N = 100000;
	int percent = 0;
	for (int i = 0; i < N; i++){
		PrintPercent((double)i/N, percent);
		double t = fmc->UniformDist(0, fActiveTime); // dice start time
		int polarisation = fmc->DicePolarisation(fParticleName); // dice polarisation
		double x, y, z;
		RandomPointInSourceVolume(x, y, z); // dice point in source volume
		TParticle *p = TParticleSource::CreateParticle(t, x, y, z, 0, 0, 0, polarisation); // create dummy particle with Ekin = 0
		double V = p->Hstart(); // potential at particle position equals its total energy
		if (V < MinPot)
			MinPot = V; // remember minimal potential
		delete p;
		ParticleCounter--;
	}
	cout << " minimal potential = " << MinPot << "eV\n";
}

void GetFinalDecayProducts(Int_t ind, IlcStack & stack , TArrayI & ar){

  // Recursive algorithm to get the final decay products of a particle
  //
  // ind is the index of the particle in the IlcStack
  // stack is the particle stack from the generator
  // ar contains the indexes of the final decay products
  // ar[0] is the number of final decay products

  if (ind<0 || ind>stack.GetNtrack()) {
    cerr << "Invalid index of the particle " << ind << endl;
    return;
  } 
  if (ar.GetSize()==0) {
    ar.Set(10);
    ar[0] = 0;
  }

  TParticle * part = stack.Particle(ind);

  Int_t iFirstDaughter = part->GetFirstDaughter();
  if( iFirstDaughter<0) {
    // This particle is a final decay product, add its index to the array
    ar[0]++;
    if (ar.GetSize() <= ar[0]) ar.Set(ar.GetSize()+10); // resize if needed
    ar[ar[0]] = ind;
    return;
  } 

  Int_t iLastDaughter = part->GetLastDaughter();

  for (Int_t id=iFirstDaughter; id<=iLastDaughter;id++) {
    // Now search for final decay products of the daughters
    GetFinalDecayProducts(id,stack,ar);
  }
}

//.........这里部分代码省略.........
  TH1F* hPtMUON = CreateHisto("hPtMUON", "MUON", 100, 0, 20, 
			      "p_{t} [GeV/c]", "N");

  // V0s and cascades
  TH1F* hMassK0 = CreateHisto("hMassK0", "K^{0}", 100, 0.4, 0.6, 
			      "M(#pi^{+}#pi^{-}) [GeV/c^{2}]", "N");
  TH1F* hMassLambda = CreateHisto("hMassLambda", "#Lambda", 100, 1.0, 1.2, 
				  "M(p#pi^{-}) [GeV/c^{2}]", "N");
  TH1F* hMassLambdaBar = CreateHisto("hMassLambdaBar", "#bar{#Lambda}", 
				     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) {

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

//.........这里部分代码省略.........
  std::ofstream myfile;
  myfile.open("_txt/generation_time.txt");
  int execution_time[Nevts/1000];

  Int_t n_lepton_total=0; //counter for the total number of charged leptons (sum on the events)

  for (int ievt = 0; ievt < Nevts; ievt++) {

    if (ievt ==0) cout << "Generating event #: " << flush;
    if (ievt % 1000 == 0) {
      cout << ievt << " " << flush;
      execution_time[ievt/1000]=time(NULL);
      if (ievt !=0)
	myfile << execution_time[ievt/1000] << "\t" << ievt << "\t" << 1000/(float)(execution_time[ievt/1000]-execution_time[ievt/1000 -1])<< endl;
    }

    // ******************************************************************
    // Generate event

    pythia->GenerateEvent();

    if ( passEvtSelection(pythia)==false ) { ///Cut on W!!!!!!!!!!!!!!!
      exclEvts++;
      continue;
    }

    // Loop over particles and fill histos

    particlesArray = (TClonesArray*)pythia->GetListOfParticles();
    int N = particlesArray->GetEntriesFast();
    //cout << "<I> PYTHIA particles in event: " << N << " " << flush;

    TMCParticle *part;
    TParticle *partic;

    Int_t n_leptons = 0;

    for (int ip = 0; ip < N; ip++ ) {

      part = (TMCParticle*)particlesArray->At(ip);
      int status = part->GetKS();
      if (status!=1) continue; // Only consider final-state particles

      int parent_id = part->GetParent();
      //if (parent_id != 22) continue; //Only consider particles coming from photons

      int pdg = part->GetKF();
      double charge = PDG->GetParticle(pdg)->Charge();
      if(charge==0) continue; // only charged particles

      //if ( abs(pdg)!=11 && abs(pdg)!=13 ) continue;// only leptons
      if ( abs(pdg)!=13 ) continue;// only muons

      partic = (TParticle*)TMCParticle2TParticle(part);
      double ptPartic = partic->Pt();
      double etaPartic = partic->Eta();
      //phiPartic = partic->Phi();

      //if ( ptPartic<0.15 ) continue; //cut on Pt!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      //if (abs(etaPartic)>1.5) continue; //Fill only if |eta|<1.5!!!!!!!!!!!!!!!!!!!!!   


      // Histo filling
      h_id_part->Fill(abs(pdg));
      hdsigmadeta->Fill(etaPartic);
      

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));
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  // Get pointers to Alice detectors and Hits containers
  tofl->LoadHits();
  TTree *TH = tofl->TreeH();
  tof->SetTreeAddress();
  if (!TH) {
    cout << "<TOFquickanal> No hit tree found" << endl;
    rc = 4;
    return rc;
  }
  
  // Import the Kine Tree for the event eventNumber in the file  
  rl->LoadHeader();
  rl->LoadKinematics();
  //AliStack * stack = rl->Stack();
  
  Int_t ntracks = TH->GetEntries();
  cout<<" ntracks = "<<ntracks<<endl;
  
  AliTOFhitT0 *tofHit;
  
  // Start loop on tracks in the hits containers
  for (Int_t track=0; track<ntracks;track++) {
    
    tof->ResetHits();
    TH->GetEvent(track);
    
    for(tofHit=(AliTOFhitT0*)tof->FirstHit(track); tofHit; tofHit=(AliTOFhitT0*)tof->NextHit()) {
      
      Float_t toflight = tofHit->GetTof();
      toflight        *= 1.E+09;  // conversion from s to ns
      Double_t tofmom  = tofHit->GetMom();
      
      Int_t ipart = tofHit->Track();
      TParticle *particle = gAlice->Particle(ipart);
      if (particle->GetFirstMother() < 0) {
	htofprim->Fill(toflight);
	htofmom->Fill(tofmom); 
      } else {
	htofsec->Fill(toflight); 
      }
      
      Double_t zcoor = tofHit->Z();
      hzcoor->Fill(zcoor);
      
      Double_t incangle = tofHit->GetIncA();
      hincangle->Fill(incangle);
      
      Double_t xcoor  = particle->Vx();
      Double_t ycoor  = particle->Vy();
      Double_t radius = TMath::Sqrt(xcoor*xcoor+ycoor*ycoor);
      if (particle->GetFirstMother() >= 0) hradius->Fill(radius);
      
      Double_t prodmom = particle->P();        
      if (prodmom!=0.) {
	Double_t dummy = (particle->Pz())/prodmom;
	Double_t prodthe = TMath::ACos(dummy);
	prodthe *= 57.29578; // conversion from rad to deg
	if (particle->GetFirstMother() < 0) hprodthe->Fill(prodthe);
      } // theta at production vertex
      
      if (particle->GetFirstMother() < 0) {         
	hprodmom->Fill(prodmom);
	Double_t dummypx = particle->Px();
	Double_t dummypy = particle->Py();
	Double_t prodphi = TMath::ATan2(dummypy,dummypx);
	prodphi *= 57.29578; // conversion from rad to deg

//.........这里部分代码省略.........
  esd->ReadFromTree(tree);
  
  Int_t nTrackTotalMC = 0;
  Int_t nTrackFound = 0;
  Int_t nTrackImperfect = 0;
  Int_t nTrackPerfect = 0;
  Int_t nNoMCTrack = 0;

  
  for(Int_t iEv =0; iEv<tree->GetEntries(); iEv++){
    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;

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

//!PG main function
int 
  selector (TChain * tree, histos & plots, int if_signal)
{
 
 plots.v_hardTAGPt = -99;
 plots.v_softTAGPt = -99;
 plots.v_TAGDProdEta = -99;
 plots.v_TAGDeta = -99;
 plots.v_TAGMinv = -99;
 plots.v_LepLep = -99;
 plots.v_hardLEPPt = -99;
 plots.v_softLEPPt = -99;
 plots.v_LEPDPhi = -99;
 plots.v_LEPDEta = -99;
 plots.v_LEPDR = -99;
 plots.v_LEPMinv = -99;
 plots.v_LEPProdCharge = -99;
 plots.v_hardLEPCharge = -99;
 plots.v_softLEPCharge = -99;
 plots.v_MET = -99;
 
 plots.v_ojets = -99 ;
 plots.v_ojetsCJV = -99 ;
 plots.v_ojetsRegionalCJV = -99 ;
 
 plots.v_ojetsZepp_01 = -99 ;
 plots.v_ojetsZepp_02 = -99 ;
 plots.v_ojetsZepp_03 = -99 ;
 plots.v_ojetsZepp_04 = -99 ;
 plots.v_ojetsZepp_05 = -99 ;
 plots.v_ojetsZepp_06 = -99 ;
 plots.v_ojetsZepp_07 = -99 ;
 plots.v_ojetsZepp_08 = -99 ;
 plots.v_ojetsZepp_09 = -99 ;
 plots.v_ojetsZepp_10 = -99 ;
 plots.v_ojetsZepp_11 = -99 ;
 plots.v_ojetsZepp_12 = -99 ;
 plots.v_ojetsZepp_13 = -99 ;
 plots.v_ojetsZepp_14 = -99 ;
 
 plots.v_decay_Channel_e = -99 ;
 plots.v_decay_Channel_mu = -99 ;
 plots.v_decay_Channel_tau = -99 ;
 
 
 TClonesArray * genParticles = new TClonesArray ("TParticle") ;
 tree->SetBranchAddress ("genParticles", &genParticles) ;
 
 
//  TClonesArray * tagJets = new TClonesArray ("TLorentzVector") ; 
//  tree->SetBranchAddress ("tagJets", &tagJets) ;
 TClonesArray * otherJets_temp = new TClonesArray ("TLorentzVector") ;
 tree->SetBranchAddress (g_KindOfJet.c_str(), &otherJets_temp) ;
//  tree->SetBranchAddress ("otherJets", &otherJets_temp) ;
 
 
 TClonesArray * electrons = new TClonesArray ("TLorentzVector") ;
 tree->SetBranchAddress ("electrons", &electrons) ;
 TClonesArray * muons = new TClonesArray ("TLorentzVector") ;
 tree->SetBranchAddress ("muons", &muons) ;
 TClonesArray * MET = new TClonesArray ("TLorentzVector") ;
 tree->SetBranchAddress ("MET", &MET) ;
 TClonesArray * tracks = new TClonesArray ("TLorentzVector") ;
 tree->SetBranchAddress ("tracks", &tracks) ;
 
 TClonesArray * tagJets = new TClonesArray ("TLorentzVector") ; 
 TClonesArray * otherJets = new TClonesArray ("TLorentzVector") ;
 
  
 int EleId[100];
 float IsolEleSumPt_VBF[100];
 int nEle;
 int EleCharge[30];
 tree->SetBranchAddress ("nEle", &nEle) ;
 tree->SetBranchAddress ("EleId",EleId ) ;
 tree->SetBranchAddress ("IsolEleSumPt_VBF",IsolEleSumPt_VBF ) ;
 tree->SetBranchAddress ("EleCharge",EleCharge ) ;
 
 float IsolMuTr[100];
 int nMu ;
 int MuCharge[30];
 tree->SetBranchAddress ("nMu", &nMu) ;
 tree->SetBranchAddress ("IsolMuTr",IsolMuTr ) ;
 tree->SetBranchAddress ("MuCharge", MuCharge) ;
 

 int IdEvent;
 tree->SetBranchAddress ("IdEvent", &IdEvent) ;
 
 
 int nentries = (int) tree->GetEntries () ;

 
 plots.passedJetAndLepNumberSelections = 0;
 plots.analyzed = 0;
 
 plots.analyzed_ee = 0;
 plots.analyzed_mumu = 0;
 plots.analyzed_tautau = 0;
//.........这里部分代码省略.........

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

//.........这里部分代码省略.........
          cSum.nRowPatt = cl-> GetPatternRowSpan();
          cSum.nColPatt = cl-> GetPatternColSpan();
	  DB.AccountTopology(*cl, cSum.dX, cSum.dZ, cSum.alpha, cSum.beta);
          
          GetHistoClSize(clsize,kDR,&histoArr)->Fill((rht-rcl)*1e4);
          if (cl-