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

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

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

//.........这里部分代码省略.........
  // 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) {
	  nGenCascades++;
	  selCascades.Add(particle);
	}

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

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

      n_leptons++; //Count the total number of charged hadrons in this event
      delete partic;
      // if (TMath::Abs(etaPartic)<1.) Nch++;
    } // End loop over all particles in event

    h_lep_per_ev->Fill(n_leptons);

    n_lepton_total += n_leptons;
    
    
  } // END EVENT GENERATION LOOP

  myfile.close();

  
  
  // **********************************************************************************  
  // FINAL GENERATION MESSAGES:
  
  pythia->Pystat(1);

void ExtractOutputHistos(Bool_t onlyPrims=0,Bool_t onlyPion=0,Int_t plotFlag=0) {

  //  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1);

  const Int_t nbins=20;
  Double_t ptmin=0.06;//04;
  Double_t ptmax=2.0;//GeV
  Double_t logxmin = TMath::Log10(ptmin);
  Double_t logxmax = TMath::Log10(ptmax);
  Double_t binwidth = (logxmax-logxmin)/(nbins+1);
  enum {nb=nbins+1};
  Double_t xbins[nb];
  xbins[0] = ptmin;
  for (Int_t i=1;i<=nbins;i++) {
    xbins[i] = ptmin + TMath::Power(10,logxmin+(i)*binwidth);
    //    cout<<xbins[i]<<endl;
  }
  //  TH1F *h = new TH1F("h","hist with log x axis",nbins,xbins);

  TH1F *hMultCount = new TH1F("mult","averaged multiplicity (charg. prim)",80,-4.,4.);
  hMultCount->GetXaxis()->SetTitle("eta");
  hMultCount->GetYaxis()->SetTitle("N/d#eta");

  TH1F *hAllMC = new TH1F("allMC","All Tracks MC primaries",nbins,xbins);
  TH1F *hAllFound = new TH1F("allFound","All Tracks found",nbins,xbins);
  TH1F *hImperfect = new TH1F("imperfect","Imperfect tracks",nbins,xbins);
  TH1F *hPerfect = new TH1F("perfect","Perfect tracks",nbins,xbins);
  TH1F *hEff = new TH1F("efficiency","Efficiency (Perfect tracks in \"ALL MC\")",nbins,xbins);
  TH1F *hFake = new TH1F("fake","Fake tracks (Inperfect tracks in \"ALL MC\")",nbins,xbins);
  TH1F *hPurity = new TH1F("purity","Purity (Perfect tracks in \"All Found\")",nbins,xbins);
  TH1F *hAnna = new TH1F("annaEff","AnnalisaEff ",nbins,xbins);
  TH1F *hNoMCTrack = new TH1F("noMCtrack","noMCtrack ",nbins,xbins);

  TH1F *hEta = new TH1F("","",50,-2,2);
  //  TH1F *hEtaMC = new TH1F("","",50,-2,2);

  TH2D *h2Ddca = new TH2D("dca2D","DCAvsPt2D",nbins,xbins,50,-0.05,0.05);
  TH2D *h2Dpt = new TH2D("dPt2D","dPtdvsPt2D",nbins,xbins,50,-25,25);

  // 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->LoadgAlice();
  gAlice = runLoader->GetAliRun();
  if (!gAlice) {
    Error("Check kine", "no galice object found");
    return;
  }
  runLoader->LoadHeader();
  runLoader->LoadKinematics();

  TFile* esdFile = TFile::Open("AliESDs.root");
  if (!esdFile || !esdFile->IsOpen()) {
    Error("CheckESD", "opening ESD file %s failed", "AliESDs.root");
    return;
  }
  AliESDEvent *esd = new AliESDEvent();
  TTree* tree = (TTree*) esdFile->Get("esdTree");
  if (!tree) {
    Error("CheckESD", "no ESD tree found");
    return;
  }
  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();
    

//.........这里部分代码省略.........

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

//.........这里部分代码省略.........
          
          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]);
          }
          */
          //
        }
      }
    }
    
    //    layerClus.Clear();

int main(int argc, char* argv[])
{
  TApplication theApp(srcName.Data(), &argc, argv);
//=============================================================================

  for (int i=0; i<argc; i++) cout << i << ", " << argv[i] << endl;
//=============================================================================

  if (argc<5) return -1;
  TString sPath = argv[1]; if (sPath.IsNull()) return -1;
  TString sFile = argv[2]; if (sFile.IsNull()) return -1;
  TString sJetR = argv[3]; if (sJetR.IsNull()) return -1;
  TString sSjeR = argv[4]; if (sSjeR.IsNull()) return -1;
//=============================================================================

  sPath.ReplaceAll("#", "/");
//=============================================================================

  double dJetR = -1.;
  if (sJetR=="JetR02") dJetR = 0.2;
  if (sJetR=="JetR03") dJetR = 0.3;
  if (sJetR=="JetR04") dJetR = 0.4;
  if (sJetR=="JetR05") dJetR = 0.5;

  if (dJetR<0.) return -1;
  cout << "Jet R = " << dJetR << endl;
//=============================================================================

  double dSjeR = -1.;
  if (sSjeR=="SjeR01") dSjeR = 0.1;
  if (sSjeR=="SjeR02") dSjeR = 0.2;
  if (sSjeR=="SjeR03") dSjeR = 0.3;
  if (sSjeR=="SjeR04") dSjeR = 0.4;

  if (dSjeR<0.) return -1;
  cout << "Sub-jet R = " << dSjeR << endl;
//=============================================================================

  const double dJetsPtMin  = 0.001;
  const double dCutEtaMax  = 1.6;
  const double dJetEtaMax  = 1.;
  const double dJetAreaRef = TMath::Pi() * dJetR * dJetR;

  fastjet::GhostedAreaSpec areaSpc(dCutEtaMax);
  fastjet::JetDefinition   jetsDef(fastjet::antikt_algorithm, dJetR, fastjet::BIpt_scheme, fastjet::Best);

//fastjet::AreaDefinition  areaDef(fastjet::active_area,areaSpc);
  fastjet::AreaDefinition  areaDef(fastjet::active_area_explicit_ghosts,areaSpc);

//fastjet::JetDefinition   bkgsDef(fastjet::kt_algorithm, 0.2, fastjet::BIpt_scheme, fastjet::Best);
//fastjet::AreaDefinition  aBkgDef(fastjet::active_area_explicit_ghosts, areaSpc);

  fastjet::Selector selectJet = fastjet::SelectorAbsEtaMax(dJetEtaMax);
//fastjet::Selector selectRho = fastjet::SelectorAbsEtaMax(dCutEtaMax-0.2);
//fastjet::Selector selecHard = fastjet::SelectorNHardest(2);
//fastjet::Selector selectBkg = selectRho * (!(selecHard));
//fastjet::JetMedianBackgroundEstimator bkgsEstimator(selectBkg, bkgsDef, aBkgDef);
//fastjet::Subtractor                   bkgSubtractor(&bkgsEstimator);

  fastjet::JetDefinition subjDef(fastjet::antikt_algorithm, dSjeR, fastjet::BIpt_scheme, fastjet::Best);
//=============================================================================

  std::vector<fastjet::PseudoJet> fjInput;
//=============================================================================

  TList *list = new TList();
  TH1D *hPtHat = new TH1D("hPtHat", "", 1000, 0., 1000.);

  TH1D *hJet = new TH1D("hJet", "", 1000, 0., 1000.); hJet->Sumw2(); list->Add(hJet);
  TH2D *hJetNsj = new TH2D("hJetNsj", "", 1000, 0., 1000., 101, -0.5, 100.5); hJetNsj->Sumw2(); list->Add(hJetNsj);

  TH2D *hJetIsj = new TH2D("hJetIsj", "", 1000, 0., 1000., 1000, 0., 1000.); hJetIsj->Sumw2(); list->Add(hJetIsj);
  TH2D *hJet1sj = new TH2D("hJet1sj", "", 1000, 0., 1000., 1000, 0., 1000.); hJet1sj->Sumw2(); list->Add(hJet1sj);
  TH2D *hJet2sj = new TH2D("hJet2sj", "", 1000, 0., 1000., 1000, 0., 1000.); hJet2sj->Sumw2(); list->Add(hJet2sj);
  TH2D *hJetDsj = new TH2D("hJetDsj", "", 1000, 0., 1000., 1000, 0., 1000.); hJetDsj->Sumw2(); list->Add(hJetDsj);

  TH2D *hJetIsz = new TH2D("hJetIsz", "", 1000, 0., 1000., 120, 0., 1.2); hJetIsz->Sumw2(); list->Add(hJetIsz);
  TH2D *hJet1sz = new TH2D("hJet1sz", "", 1000, 0., 1000., 120, 0., 1.2); hJet1sz->Sumw2(); list->Add(hJet1sz);
  TH2D *hJet2sz = new TH2D("hJet2sz", "", 1000, 0., 1000., 120, 0., 1.2); hJet2sz->Sumw2(); list->Add(hJet2sz);
  TH2D *hJetDsz = new TH2D("hJetDsz", "", 1000, 0., 1000., 120, 0., 1.2); hJetDsz->Sumw2(); list->Add(hJetDsz);
//=============================================================================

  AliRunLoader *rl = AliRunLoader::Open(Form("%s/galice.root",sPath.Data())); if (!rl) return -1;

  if (rl->LoadHeader()) return -1;
  if (rl->LoadKinematics("READ")) return -1;
//=============================================================================

  for (Int_t iEvent=0; iEvent<rl->GetNumberOfEvents(); iEvent++) {
    fjInput.resize(0);
    if (rl->GetEvent(iEvent)) continue;
//=============================================================================

    AliStack  *pStack  = rl->Stack();     if (!pStack)  continue;
    AliHeader *pHeader = rl->GetHeader(); if (!pHeader) continue;
//=============================================================================

    AliGenPythiaEventHeader *pHeadPy = (AliGenPythiaEventHeader*)pHeader->GenEventHeader();

    if (!pHeadPy) continue;
//.........这里部分代码省略.........