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
0005 #include "TGraphErrors.h"
0006 #include "TF1.h"
0007 #include "TRandom.h"
0008 #include "TCanvas.h"
0009 #include "TLegend.h"
0010 #include "TMath.h"
0011 #define mpi 0.139  // 1.864 GeV/c^2
0012
0013 void Tracking_Performances(TString filename="tracking_output",TString particle="pi-", double mom=0.1, Double_t pTcut = 0.15, bool truth_seeding=false)
0014 {
0015
0016   // style of the plot
0017    gStyle->SetPalette(1);
0018    gStyle->SetOptTitle(1);
0019    gStyle->SetTitleOffset(.85,"X");gStyle->SetTitleOffset(.85,"Y");
0020    gStyle->SetTitleSize(.05,"X");gStyle->SetTitleSize(.05,"Y");
0021    gStyle->SetLabelSize(.04,"X");gStyle->SetLabelSize(.04,"Y");
0022    gStyle->SetHistLineWidth(2);
0023    gStyle->SetOptFit(1);
0024    gStyle->SetOptStat(1);
0025
0026    bool debug=true;
0027   // Tree with reconstructed tracks
0028    const int nbins_eta = 5;
0029    int theta_val[nbins_eta+1] ={3,50,45,135,130,177};
0030    int nfiles = 100;
0031    double eta[nbins_eta+1]={-3.5,-2.5,-1.0,1.0,2.5,3.5};
0032    double pt[nbins_eta+1]={0.5,1.0,2.0,5.0,10.0,20.1};
0033    TH1D *histp[nbins_eta];
0034
0035     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);
0036     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);
0037
0038    for (int i=0; i<nbins_eta; i++){
0039    histp[i] = new TH1D(Form("hist_etabin%d",i),Form("hist_etabin%d",i),600,-0.3,0.3);
0040    histp[i]->SetTitle(Form("%1.1f < #eta < %1.1f && p = %1.1f ",eta[i],eta[i+1],mom));
0041    histp[i]->SetName(Form("hist_mom_%1.1f_%1.1f_pmax_%1.1f",mom,eta[i],eta[i+1]));
0042    }
0043
0044    TFile* file = TFile::Open(filename.Data());
0045    if (!file) {printf("file not found !!!"); return;}
0046    TTreeReader myReader("events", file); // name of tree and file
0047    if (debug) cout<<"Filename: "<<file->GetName()<<"\t NEvents: "<<myReader.GetEntries()<<endl;
0048
0049    TTree *tree = dynamic_cast<TTree*>(file->Get("events"));
0050    if (!(tree->GetBranch("CentralCKFSeededTrackParameters") || tree->GetBranch("CentralCKFTruthSeededTrackParameters"))) {
0051      cerr << "Found neither CentralCKFSeededTrackParameters nor CentralCKFTruthSeededTrackParameters!" << endl;
0052      return;
0053    }
0054    bool is_old_style = tree->GetBranch("CentralCKFSeededTrackParameters");
0055
0056    TString dir = "";
0057    TString dist_dir_mom = ""; TString dist_dir_dca = "";
0058    TString tag = "";
0059    if (truth_seeding) {
0060       dist_dir_mom = "mom_resol_truth"; dist_dir_dca = "dca_resol_truth"; dir = "truthseed";
0061       if (is_old_style) {
0062         tag = "";
0063       } else {
0064         tag = "TruthSeeded";
0065       }
0066    } else {
0067       dist_dir_mom = "mom_resol_realseed"; dist_dir_dca = "dca_resol_realseed"; dir = "realseed";
0068       if (is_old_style) {
0069         tag = "Seeded";
0070       } else {
0071         tag = "";
0072       }
0073    }
0074
0075    // MC and Reco information
0076    TTreeReaderArray<Float_t> charge(myReader, "MCParticles.charge");
0077    TTreeReaderArray<Double_t> vx_mc(myReader, "MCParticles.vertex.x");
0078    TTreeReaderArray<Double_t> vy_mc(myReader, "MCParticles.vertex.y");
0079    TTreeReaderArray<Double_t> vz_mc(myReader, "MCParticles.vertex.z");
0080    TTreeReaderArray<Float_t> px_mc(myReader, "MCParticles.momentum.x");
0081    TTreeReaderArray<Float_t> py_mc(myReader, "MCParticles.momentum.y");
0082    TTreeReaderArray<Float_t> pz_mc(myReader, "MCParticles.momentum.z");
0083    TTreeReaderArray<Int_t> status(myReader, "MCParticles.generatorStatus");
0084    TTreeReaderArray<Int_t> pdg(myReader, "MCParticles.PDG");
0085    TTreeReaderArray<Int_t> match_flag(myReader, Form("CentralCKF%sTrackParameters.type",tag.Data()));
0086    TTreeReaderArray<Float_t> d0xy(myReader, Form("CentralCKF%sTrackParameters.loc.a",tag.Data()));
0087    TTreeReaderArray<Float_t> d0z(myReader, Form("CentralCKF%sTrackParameters.loc.b",tag.Data()));
0088    TTreeReaderArray<Float_t> theta(myReader, Form("CentralCKF%sTrackParameters.theta",tag.Data()));
0089    TTreeReaderArray<Float_t> phi(myReader, Form("CentralCKF%sTrackParameters.phi",tag.Data()));
0090    TTreeReaderArray<Float_t> qoverp(myReader, Form("CentralCKF%sTrackParameters.qOverP",tag.Data()));
0091
0092   int count =0;
0093   int matchId = 1; // Always matched track assigned the index 0
0094   while (myReader.Next())
0095   {
0096    if (match_flag.GetSize()==0) continue;  // Events with no reco tracks skip them
0097    for (int i = 0; i < matchId; ++i){
0098
0099      for (int j = 0; j < pdg.GetSize(); ++j){
0100
0101       if (status[j] !=1 && pdg.GetSize()!=1) continue;
0102       Double_t ptmc = sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]);
0103
0104       if (fabs(ptmc) < pTcut) continue;
0105
0106       Double_t pmc = (1./charge[j])*sqrt(px_mc[j]*px_mc[j]+py_mc[j]*py_mc[j]+pz_mc[j]*pz_mc[j]); // 1./(q/p); similar to prec
0107       Double_t prec = 1./qoverp[j];
0108
0109       Double_t pzrec = prec*TMath::Cos(theta[j]);  Double_t pt_rec = sqrt(prec*prec-pzrec*pzrec);
0110       Double_t pzmc = pz_mc[j];
0111
0112       Double_t etamc = -1.0*TMath::Log(TMath::Tan((TMath::ACos(pzmc/fabs(pmc)))/2));
0113       Double_t p_resol = (prec-pmc)/pmc;
0114
0115       for (int ibin=0; ibin<nbins_eta; ++ibin){
0116       if(etamc>eta[ibin] && etamc<eta[ibin+1]) histp[ibin]->Fill(p_resol);
0117       }
0118       h_d0xy_3d->Fill(d0xy[j]*0.1, etamc, ptmc); // cm
0119       h_d0z_3d->Fill(d0z[j]*0.1, etamc, ptmc); // cm
0120       } // Generated Tracks
0121      } // Reco Tracks
0122
0123    }// event loop ends
0124
0125    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");
0126    fout_mom->cd();
0127    for (int ibin=0; ibin<nbins_eta; ++ibin) histp[ibin]->Write();
0128    fout_mom->Close();
0129
0130    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");
0131    fout_dca->cd();
0132     h_d0xy_3d->Write();
0133     h_d0z_3d->Write();
0134     fout_dca->Close();
0135 }
0136
0137
0138
0139