benchmarks/tracking_performances/Tracking_Performances.C

0001 // Code to extract the Tracking Performances
0002 // Shyam Kumar; INFN Bari, Italy
0003 // shyam.kumar@ba.infn.it; shyam.kumar@cern.ch
0004 // Pull distributions: https://indico.bnl.gov/event/28544/contributions/109057/attachments/62799/108633/ePIC_Tracking_Meeting_30June2025_ShyamKumar.pdf
0005
0006 #include "TGraphErrors.h"
0007 #include "TF1.h"
0008 #include "TRandom.h"
0009 #include "TCanvas.h"
0010 #include "TLegend.h"
0011 #include "TMath.h"
0012 #define mpi 0.139  // 1.864 GeV/c^2
0013
0014 void Tracking_Performances(TString filename="tracking_output",TString particle="pi-", double mom=0.1, Double_t pTcut = 0.15, bool truth_seeding=false)
0015 {
0016
0017   // style of the plot
0018    gStyle->SetPalette(1);
0019    gStyle->SetOptTitle(1);
0020    gStyle->SetTitleOffset(.85,"X");gStyle->SetTitleOffset(.85,"Y");
0021    gStyle->SetTitleSize(.05,"X");gStyle->SetTitleSize(.05,"Y");
0022    gStyle->SetLabelSize(.04,"X");gStyle->SetLabelSize(.04,"Y");
0023    gStyle->SetHistLineWidth(2);
0024    gStyle->SetOptFit(1);
0025    gStyle->SetOptStat(1);
0026
0027    bool debug=true;
0028   // Tree with reconstructed tracks
0029    const int nbins_eta = 5;
0030    int theta_val[nbins_eta+1] ={3,50,45,135,130,177};
0031    int nfiles = 100;
0032    double eta[nbins_eta+1]={-3.5,-2.5,-1.0,1.0,2.5,3.5};
0033    double pt[nbins_eta+1]={0.5,1.0,2.0,5.0,10.0,20.1};
0034    TH1D *histp[nbins_eta], *hpull_invp[nbins_eta], *hpull_d0xy[nbins_eta], *hpull_d0z[nbins_eta], *hpull_phi[nbins_eta], *hpull_theta[nbins_eta];
0035
0036     TH3D *h_d0xy_3d= new TH3D("h_d0xy_3d","Transverse Pointing Resolution",500,-0.1,0.1,70,-3.5,3.5,201,0.,20.1);
0037     TH3D *h_d0z_3d= new TH3D("h_d0z_3d","Longitudinal Pointing Resolution",500,-0.1,0.1,70,-3.5,3.5,201,0.,20.1);
0038
0039    for (int i=0; i<nbins_eta; i++){
0040    histp[i] = new TH1D(Form("hist_etabin%d",i),Form("hist_etabin%d",i),600,-0.3,0.3);
0041    histp[i]->SetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f ",eta[i],eta[i+1],mom));
0042    histp[i]->SetName(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1]));
0043
0044    hpull_invp[i] = new TH1D(Form("hpull_invp_etabin%d",i),Form("hist_etabin%d",i),100,-10.0,10.0);
0045    hpull_invp[i]->SetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f (GeV/c);Pull (q/p);Entries (a.u.)",eta[i],eta[i+1],mom));
0046    hpull_invp[i]->SetName(Form("hpull_invp_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1]));
0047
0048    hpull_d0xy[i] = new TH1D(Form("hpull_d0xy_etabin%d",i),Form("hist_etabin%d",i),100,-10.0,10.0);
0049    hpull_d0xy[i]->SetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f (GeV/c);Pull (d0_{xy});Entries (a.u.)",eta[i],eta[i+1],mom));
0050    hpull_d0xy[i]->SetName(Form("hpull_d0xy_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1]));
0051
0052    hpull_d0z[i] = new TH1D(Form("hpull_d0z_etabin%d",i),Form("hist_etabin%d",i),100,-10.0,10.0);
0053    hpull_d0z[i]->SetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f (GeV/c);Pull (d0_{z});Entries (a.u.)",eta[i],eta[i+1],mom));
0054    hpull_d0z[i]->SetName(Form("hpull_d0z_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1]));
0055
0056    hpull_phi[i] = new TH1D(Form("hpull_phi_etabin%d",i),Form("hist_etabin%d",i),100,-10.0,10.0);
0057    hpull_phi[i]->SetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f (GeV/c);Pull (#phi);Entries (a.u.)",eta[i],eta[i+1],mom));
0058    hpull_phi[i]->SetName(Form("hpull_phi_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1]));
0059
0060    hpull_theta[i] = new TH1D(Form("hpull_theta_etabin%d",i),Form("hist_etabin%d",i),100,-10.0,10.0);
0061    hpull_theta[i]->SetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f (GeV/c);Pull (#theta);Entries (a.u.)",eta[i],eta[i+1],mom));
0062    hpull_theta[i]->SetName(Form("hpull_theta_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1]));
0063    }
0064
0065    TFile* file = TFile::Open(filename.Data());
0066    if (!file) {printf("file not found !!!"); return;}
0067    TTreeReader myReader("events", file); // name of tree and file
0068    if (debug) cout<<"Filename: "<<file->GetName()<<"\t NEvents: "<<myReader.GetEntries()<<endl;
0069
0070    TTree *tree = dynamic_cast<TTree*>(file->Get("events"));
0071    if (!(tree->GetBranch("CentralCKFSeededTrackParameters") || tree->GetBranch("CentralCKFTruthSeededTrackParameters"))) {
0072      cerr << "Found neither CentralCKFSeededTrackParameters nor CentralCKFTruthSeededTrackParameters!" << endl;
0073      return;
0074    }
0075    bool is_old_style = tree->GetBranch("CentralCKFSeededTrackParameters");
0076
0077    TString dir = "";
0078    TString dist_dir_mom = ""; TString dist_dir_dca = "";
0079    TString tag = "";
0080    if (truth_seeding) {
0081       dist_dir_mom = "mom_resol_truth"; dist_dir_dca = "dca_resol_truth"; dir = "truthseed";
0082       if (is_old_style) {
0083         tag = "";
0084       } else {
0085         tag = "TruthSeeded";
0086       }
0087    } else {
0088       dist_dir_mom = "mom_resol_realseed"; dist_dir_dca = "dca_resol_realseed"; dir = "realseed";
0089       if (is_old_style) {
0090         tag = "Seeded";
0091       } else {
0092         tag = "";
0093       }
0094    }
0095
0096    // MC and Reco information
0097    TTreeReaderArray<Float_t> charge(myReader, "MCParticles.charge");
0098    TTreeReaderArray<Double_t> vx_mc(myReader, "MCParticles.vertex.x");
0099    TTreeReaderArray<Double_t> vy_mc(myReader, "MCParticles.vertex.y");
0100    TTreeReaderArray<Double_t> vz_mc(myReader, "MCParticles.vertex.z");
0101    TTreeReaderArray<Double_t> px_mc(myReader, "MCParticles.momentum.x");
0102    TTreeReaderArray<Double_t> py_mc(myReader, "MCParticles.momentum.y");
0103    TTreeReaderArray<Double_t> pz_mc(myReader, "MCParticles.momentum.z");
0104    TTreeReaderArray<Int_t> status(myReader, "MCParticles.generatorStatus");
0105    TTreeReaderArray<Int_t> pdg(myReader, "MCParticles.PDG");
0106    TTreeReaderArray<Int_t> match_flag(myReader, Form("CentralCKF%sTrackParameters.type",tag.Data()));
0107    TTreeReaderArray<Float_t> d0xy(myReader, Form("CentralCKF%sTrackParameters.loc.a",tag.Data()));
0108    TTreeReaderArray<Float_t> d0z(myReader, Form("CentralCKF%sTrackParameters.loc.b",tag.Data()));
0109    TTreeReaderArray<Float_t> theta(myReader, Form("CentralCKF%sTrackParameters.theta",tag.Data()));
0110    TTreeReaderArray<Float_t> phi(myReader, Form("CentralCKF%sTrackParameters.phi",tag.Data()));
0111    TTreeReaderArray<Float_t> qoverp(myReader, Form("CentralCKF%sTrackParameters.qOverP",tag.Data()));
0112    TTreeReaderArray<std::array<float, 21>> rcTrkCov(myReader, Form("CentralCKF%sTrackParameters.covariance.covariance[21]",tag.Data()));
0113
0114   int count =0;
0115   int matchId = 1; // Always matched track assigned the index 0
0116   while (myReader.Next())
0117   {
0118    if (match_flag.GetSize()==0) continue;  // Events with no reco tracks skip them
0119    for (int i = 0; i < matchId; ++i){
0120
0121      for (int j = 0; j < pdg.GetSize(); ++j){
0122
0123       if (status[j] !=1 && pdg.GetSize()!=1) continue;
0124       TVector3 mom_MC(px_mc[j],py_mc[j],pz_mc[j]);
0125       Double_t theta_mc = mom_MC.Theta();
0126       Double_t phi_mc = mom_MC.Phi();
0127       Double_t ptmc = mom_MC.Pt();
0128       Double_t etamc = mom_MC.Eta();
0129
0130      if (fabs(ptmc) < pTcut) continue;
0131       Double_t pmc = mom_MC.Mag() / charge[j]; // 1./(q/p); similar to prec
0132       Double_t prec = 1./qoverp[j];
0133
0134       Double_t pzrec = prec*TMath::Cos(theta[j]);  Double_t pt_rec = sqrt(prec*prec-pzrec*pzrec);
0135       Double_t pzmc = pz_mc[j];
0136
0137      // Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pzmc/fabs(pmc)))/2));
0138       Double_t p_resol = (prec-pmc)/pmc;
0139       // Accessing the pull distributions
0140       std::array<float, 21>& cov = rcTrkCov.At(j); // access covariance
0141       Double_t pull_invmom = (fabs(1./prec)-fabs(1./pmc))/sqrt(cov[14]); // cov[14] = sigma_1/p^2
0142       Double_t pull_dcaxy = d0xy[j]/sqrt(cov[0]); // cov[0] = sigma_l0^2
0143       Double_t pull_dcaz = d0z[j]/sqrt(cov[2]);  // cov[2] = sigma_l1^2
0144       Double_t pull_phi = (phi[j]-phi_mc)/sqrt(cov[5]); // cov[5] = sigma_phi^2
0145       Double_t pull_theta = (theta[j]-theta_mc)/sqrt(cov[9]); // cov[9] = sigma_theta^2
0146
0147
0148       for (int ibin=0; ibin<nbins_eta; ++ibin){
0149       if(etamc>eta[ibin] && etamc<eta[ibin+1]){
0150       histp[ibin]->Fill(p_resol);
0151       hpull_invp[ibin]->Fill(pull_invmom);
0152       hpull_d0xy[ibin]->Fill(pull_dcaxy);
0153       hpull_d0z[ibin]->Fill(pull_dcaz);
0154       hpull_phi[ibin]->Fill(pull_phi);
0155       hpull_theta[ibin]->Fill(pull_theta);
0156       }
0157       }
0158       h_d0xy_3d->Fill(d0xy[j]*0.1, etamc, ptmc); // cm
0159       h_d0z_3d->Fill(d0z[j]*0.1, etamc, ptmc); // cm
0160       } // Generated Tracks
0161      } // Reco Tracks
0162
0163    }// event loop ends
0164
0165    TFile *fout_mom = new TFile(Form("%s/%s/mom/Performances_mom_%1.1f_%s_%s.root",dir.Data(),particle.Data(),mom,dist_dir_mom.Data(),particle.Data()),"recreate");
0166    fout_mom->cd();
0167    for (int ibin=0; ibin<nbins_eta; ++ibin){
0168    histp[ibin]->Write();
0169    hpull_invp[ibin]->Write();
0170    hpull_d0xy[ibin]->Write();
0171    hpull_d0z[ibin]->Write();
0172    hpull_phi[ibin]->Write();
0173    hpull_theta[ibin]->Write();
0174    }
0175    fout_mom->Close();
0176
0177    TFile *fout_dca = new TFile(Form("%s/%s/dca/Performances_dca_%1.1f_%s_%s.root",dir.Data(),particle.Data(),mom,dist_dir_dca.Data(),particle.Data()),"recreate");
0178    fout_dca->cd();
0179     h_d0xy_3d->Write();
0180     h_d0z_3d->Write();
0181     fout_dca->Close();
0182 }
0183
0184
0185
0186