Tracking/LowMomentumTracking/track_eff_fixed_mom.C

0001 //------------------

0002 void track_eff_fixed_mom(std::string gen_part = "muminus"){
0003
0004     //Define Style

0005     gStyle->SetOptStat(0);
0006     gStyle->SetPadBorderMode(0);
0007     gStyle->SetFrameBorderMode(0);
0008     gStyle->SetFrameLineWidth(2);
0009     gStyle->SetLabelSize(0.035,"X");
0010     gStyle->SetLabelSize(0.035,"Y");
0011     //gStyle->SetLabelOffset(0.01,"X");

0012     //gStyle->SetLabelOffset(0.01,"Y");

0013     gStyle->SetTitleXSize(0.04);
0014     gStyle->SetTitleXOffset(0.9);
0015     gStyle->SetTitleYSize(0.04);
0016     gStyle->SetTitleYOffset(0.9);
0017
0018     //Define simulated momentum points

0019     const int N_mom = 6;
0020     double fixed_mom[N_mom] = {0.5,1.0,2.0,5.0,10.0,20.0}; //In GeV/c

0021
0022     //Tracker Efficiency Histograms

0023     TH1 *h_eta_gen[N_mom]; //Generated particle eta distribution

0024     TH1 *h_eta_gen_ts_eff[N_mom]; //Generated particle eta distribution w/ at least 1 truth-seeded track found

0025     TH1 *h_eta_ts_ratio[N_mom]; //Efficiency ratio for truth-seeded tracking

0026     TH1 *h_eta_gen_rs_eff[N_mom]; //Generated particle eta distribution w/ at least 1 real-seeded track found

0027     TH1 *h_eta_rs_ratio[N_mom]; //Efficiency ratio for real-seeded tracking

0028
0029     //Only keep events with -3.5<eta_true<3.5

0030     TH1 *h_phi_gen[N_mom]; //Generated particle phi distribution

0031     TH1 *h_phi_gen_ts_eff[N_mom]; //Generated particle phi distribution w/ at least 1 truth-seeded track found

0032     TH1 *h_phi_ts_ratio[N_mom]; //Efficiency ratio for truth-seeded tracking

0033     TH1 *h_phi_gen_rs_eff[N_mom]; //Generated particle phi distribution w/ at least 1 real-seeded track found

0034     TH1 *h_phi_rs_ratio[N_mom]; //Efficiency ratio for real-seeded tracking

0035
0036     //Loop over momentum points

0037     for(int imom = 0; imom < N_mom; imom++){
0038
0039         cout<<endl<<"Analyzing momentum = "<< fixed_mom[imom] <<"GeV/c point!"<< endl;
0040
0041         //Define histograms

0042         //---Eta---

0043         h_eta_gen[imom] = new TH1D(Form("h_eta_gen[%d]",imom),Form("Tracker Efficiency vs. generated particle #eta: p_{gen} = %.2f GeV/c",fixed_mom[imom]),100,-4,4);
0044         h_eta_gen[imom]->GetXaxis()->SetTitle("#eta");h_eta_gen[imom]->GetXaxis()->CenterTitle();
0045         h_eta_gen[imom]->GetYaxis()->SetTitle("Efficiency");h_eta_gen[imom]->GetYaxis()->CenterTitle();
0046         h_eta_gen[imom]->SetLineColor(kBlue);h_eta_gen[imom]->SetLineWidth(2);
0047
0048         h_eta_gen_ts_eff[imom] = new TH1D(Form("h_eta_gen_ts_eff[%d]",imom),Form("Tracker Efficiency vs. generated particle #eta: p_{gen} = %.2f GeV/c",fixed_mom[imom]),100,-4,4);
0049         h_eta_gen_ts_eff[imom]->GetXaxis()->SetTitle("#eta");h_eta_gen_ts_eff[imom]->GetXaxis()->CenterTitle();
0050         h_eta_gen_ts_eff[imom]->GetYaxis()->SetTitle("Efficiency");h_eta_gen_ts_eff[imom]->GetYaxis()->CenterTitle();
0051         h_eta_gen_ts_eff[imom]->SetLineColor(kBlue);h_eta_gen_ts_eff[imom]->SetLineWidth(2);
0052
0053         h_eta_gen_rs_eff[imom] = new TH1D(Form("h_eta_gen_rs_eff[%d]",imom),Form("Tracker Efficiency vs. generated particle #eta: p_{gen} = %.2f GeV/c",fixed_mom[imom]),100,-4,4);
0054         h_eta_gen_rs_eff[imom]->GetXaxis()->SetTitle("#eta");h_eta_gen_rs_eff[imom]->GetXaxis()->CenterTitle();
0055         h_eta_gen_rs_eff[imom]->GetYaxis()->SetTitle("Efficiency");h_eta_gen_rs_eff[imom]->GetYaxis()->CenterTitle();
0056         h_eta_gen_rs_eff[imom]->SetLineColor(kGreen);h_eta_gen_rs_eff[imom]->SetLineWidth(2);
0057
0058     //---Phi---

0059      h_phi_gen[imom] = new TH1D(Form("h_phi_gen[%d]",imom),Form("Tracker Efficiency vs. generated particle #phi: p_{gen} = %.2f GeV/c",fixed_mom[imom]),100,-3.40,3.40);
0060         h_phi_gen[imom]->GetXaxis()->SetTitle("#phi");h_phi_gen[imom]->GetXaxis()->CenterTitle();
0061         h_phi_gen[imom]->GetYaxis()->SetTitle("Efficiency");h_phi_gen[imom]->GetYaxis()->CenterTitle();
0062         h_phi_gen[imom]->SetLineColor(kBlue);h_phi_gen[imom]->SetLineWidth(2);
0063
0064         h_phi_gen_ts_eff[imom] = new TH1D(Form("h_phi_gen_ts_eff[%d]",imom),Form("Tracker Efficiency vs. generated particle #phi: p_{gen} = %.2f GeV/c",fixed_mom[imom]),100,-3.40,3.40);
0065         h_phi_gen_ts_eff[imom]->GetXaxis()->SetTitle("#phi");h_phi_gen_ts_eff[imom]->GetXaxis()->CenterTitle();
0066         h_phi_gen_ts_eff[imom]->GetYaxis()->SetTitle("Efficiency");h_phi_gen_ts_eff[imom]->GetYaxis()->CenterTitle();
0067         h_phi_gen_ts_eff[imom]->SetLineColor(kBlue);h_phi_gen_ts_eff[imom]->SetLineWidth(2);
0068
0069         h_phi_gen_rs_eff[imom] = new TH1D(Form("h_phi_gen_rs_eff[%d]",imom),Form("Tracker Efficiency vs. generated particle #phi: p_{gen} = %.2f GeV/c",fixed_mom[imom]),100,-3.40,3.40);
0070         h_phi_gen_rs_eff[imom]->GetXaxis()->SetTitle("#phi");h_phi_gen_rs_eff[imom]->GetXaxis()->CenterTitle();
0071         h_phi_gen_rs_eff[imom]->GetYaxis()->SetTitle("Efficiency");h_phi_gen_rs_eff[imom]->GetYaxis()->CenterTitle();
0072         h_phi_gen_rs_eff[imom]->SetLineColor(kGreen);h_phi_gen_rs_eff[imom]->SetLineWidth(2);
0073
0074         //File information

0075         float x_gen;
0076         float y_gen;
0077         float z_gen;
0078         TString run_name;
0079         TString path = "./output/";
0080         int pid_code;
0081
0082         if(gen_part=="muminus"){
0083
0084             pid_code = 13;
0085
0086             if(imom==0){
0087                 run_name = "eicrecon_out_500MeV.root";
0088                 x_gen = 0;
0089                 y_gen = 0;
0090                 z_gen = 0;
0091             }
0092             else if(imom==1){
0093                 run_name = "eicrecon_out_1GeV.root";
0094                 x_gen = 0;
0095                 y_gen = 0;
0096                 z_gen = 0;
0097             }
0098             else if(imom==2){
0099                 run_name = "eicrecon_out_2GeV.root";
0100                 x_gen = 0;
0101                 y_gen = 0;
0102                 z_gen = 0;
0103             }
0104             else if(imom==3){
0105                 run_name = "eicrecon_out_5GeV.root";
0106                 x_gen = 0;
0107                 y_gen = 0;
0108                 z_gen = 0;
0109             }
0110         else if(imom==4){
0111                 run_name = "eicrecon_out_10GeV.root";
0112                 x_gen = 0;
0113                 y_gen = 0;
0114                 z_gen = 0;
0115             }
0116         else if(imom==5){
0117                 run_name = "eicrecon_out_20GeV.root";
0118                 x_gen = 0;
0119                 y_gen = 0;
0120                 z_gen = 0;
0121             }
0122         }
0123
0124         //Open File

0125         TString input = path + run_name;
0126         TFile *f = new TFile(input.Data());
0127         TTree *tree = (TTree*) f->Get("events");
0128
0129         //Create Array Reader

0130         TTreeReader tr(tree);
0131
0132         TTreeReaderArray<int>   gen_status(tr, "MCParticles.generatorStatus");
0133         TTreeReaderArray<int>   gen_pid(tr, "MCParticles.PDG");
0134         TTreeReaderArray<float> gen_px(tr, "MCParticles.momentum.x");
0135         TTreeReaderArray<float> gen_py(tr, "MCParticles.momentum.y");
0136         TTreeReaderArray<float> gen_pz(tr, "MCParticles.momentum.z");
0137         TTreeReaderArray<double> gen_mass(tr, "MCParticles.mass"); //Not important here

0138         TTreeReaderArray<float> gen_charge(tr, "MCParticles.charge");
0139         TTreeReaderArray<double> gen_vx(tr, "MCParticles.vertex.x");
0140         TTreeReaderArray<double> gen_vy(tr, "MCParticles.vertex.y");
0141         TTreeReaderArray<double> gen_vz(tr, "MCParticles.vertex.z");
0142
0143         TTreeReaderArray<float> track_ts_qoverp(tr, "CentralCKFTruthSeededTrackParameters.qOverP");
0144         TTreeReaderArray<float> track_ts_theta(tr, "CentralCKFTruthSeededTrackParameters.theta");
0145         TTreeReaderArray<float> track_ts_phi(tr, "CentralCKFTruthSeededTrackParameters.phi");
0146         TTreeReaderArray<float> track_ts_loca(tr, "CentralCKFTruthSeededTrackParameters.loc.a");
0147         TTreeReaderArray<float> track_ts_locb(tr, "CentralCKFTruthSeededTrackParameters.loc.b");
0148
0149         TTreeReaderArray<float> track_rs_qoverp(tr, "CentralCKFTrackParameters.qOverP");
0150         TTreeReaderArray<float> track_rs_theta(tr, "CentralCKFTrackParameters.theta");
0151         TTreeReaderArray<float> track_rs_phi(tr, "CentralCKFTrackParameters.phi");
0152         TTreeReaderArray<float> track_rs_loca(tr, "CentralCKFTrackParameters.loc.a");
0153         TTreeReaderArray<float> track_rs_locb(tr, "CentralCKFTrackParameters.loc.b");
0154
0155         //Other variables

0156         TLorentzVector gen_vec;
0157         TVector3 gen_vertex;
0158         float charge;
0159         bool found_primary(false);
0160         int counter(0);
0161
0162         //Loop over events

0163         while (tr.Next()) {
0164
0165         if(counter%100==0) cout<<"Analyzing event "<<counter<<endl;
0166         counter++;
0167
0168             //Reset variables

0169             found_primary = false;
0170
0171             //Loop over generated particles, select primary particle (assuming single particle)

0172             for(int igen=0;igen<gen_status.GetSize();igen++){
0173                 //PID requirement so that background particles are not used, but may not be needed because of 'break'

0174                 if(gen_status[igen]==1 && gen_pid[igen]==pid_code){
0175                     gen_vec.SetXYZM(gen_px[igen],gen_py[igen],gen_pz[igen],gen_mass[igen]);
0176                     gen_vertex.SetXYZ(gen_vx[igen],gen_vy[igen],gen_vz[igen]);
0177                     charge = gen_charge[igen];
0178                     found_primary = true;
0179                     break;
0180                 }
0181             }
0182
0183             //Require primary particle for all track results

0184             if(found_primary){
0185
0186                 h_eta_gen[imom]->Fill(gen_vec.Eta());
0187         if( std::abs(gen_vec.Eta())<3.0 ) h_phi_gen[imom]->Fill(gen_vec.Phi());
0188
0189                 //Truth-seeded tracking

0190                 int track_ts_mult = track_ts_qoverp.GetSize();
0191                 if(track_ts_mult>0) {
0192                     h_eta_gen_ts_eff[imom]->Fill(gen_vec.Eta());
0193             if( std::abs(gen_vec.Eta())<3.0 ) h_phi_gen_ts_eff[imom]->Fill(gen_vec.Phi());
0194                 }
0195
0196                 //Real seeded tracking

0197                 int track_rs_mult = track_rs_qoverp.GetSize();
0198                 if(track_rs_mult>0) {
0199                     h_eta_gen_rs_eff[imom]->Fill(gen_vec.Eta());
0200             if( std::abs(gen_vec.Eta())<3.0 ) h_phi_gen_rs_eff[imom]->Fill(gen_vec.Phi());
0201                 }
0202             }
0203         } // End loop over events

0204
0205         //Divide histograms

0206         h_eta_gen_ts_eff[imom] = (TH1*) h_eta_gen_ts_eff[imom]->Clone(Form("h_eta_gen_ts_eff[%d]",imom));
0207         h_eta_gen_ts_eff[imom]->Divide(h_eta_gen[imom]);
0208
0209         h_eta_gen_rs_eff[imom] = (TH1*) h_eta_gen_rs_eff[imom]->Clone(Form("h_eta_gen_rs_eff[%d]",imom));
0210         h_eta_gen_rs_eff[imom]->Divide(h_eta_gen[imom]);
0211
0212     h_phi_gen_ts_eff[imom] = (TH1*) h_phi_gen_ts_eff[imom]->Clone(Form("h_phi_gen_ts_eff[%d]",imom));
0213         h_phi_gen_ts_eff[imom]->Divide(h_phi_gen[imom]);
0214
0215         h_phi_gen_rs_eff[imom] = (TH1*) h_phi_gen_rs_eff[imom]->Clone(Form("h_phi_gen_rs_eff[%d]",imom));
0216         h_phi_gen_rs_eff[imom]->Divide(h_phi_gen[imom]);
0217
0218     } // End loop over momentum points

0219
0220
0221     //Make plots

0222     TCanvas *c1_ts = new TCanvas("c1_ts");
0223     c1_ts->Divide(2,2);
0224     c1_ts->cd(1); h_eta_gen_ts_eff[0]->Draw();
0225     c1_ts->cd(2); h_eta_gen_ts_eff[1]->Draw();
0226     c1_ts->cd(3); h_eta_gen_ts_eff[2]->Draw();
0227     c1_ts->cd(4); h_eta_gen_ts_eff[3]->Draw();
0228
0229     TCanvas *c1_rs = new TCanvas("c1_rs");
0230     c1_rs->Divide(2,2);
0231     c1_rs->cd(1); h_eta_gen_rs_eff[0]->Draw();
0232     c1_rs->cd(2); h_eta_gen_rs_eff[1]->Draw();
0233     c1_rs->cd(3); h_eta_gen_rs_eff[2]->Draw();
0234     c1_rs->cd(4); h_eta_gen_rs_eff[3]->Draw();
0235
0236     TCanvas *c2_ts = new TCanvas("c2_ts");
0237     c2_ts->Divide(2,2);
0238     c2_ts->cd(1); h_eta_gen_ts_eff[4]->Draw();
0239     c2_ts->cd(2); h_eta_gen_ts_eff[5]->Draw();
0240
0241     TCanvas *c2_rs = new TCanvas("c2_rs");
0242     c2_rs->Divide(2,2);
0243     c2_rs->cd(1); h_eta_gen_rs_eff[4]->Draw();
0244     c2_rs->cd(2); h_eta_gen_rs_eff[5]->Draw();
0245
0246     TCanvas *c3_ts = new TCanvas("c3_ts");
0247     c3_ts->Divide(2,2);
0248     c3_ts->cd(1); h_phi_gen_ts_eff[0]->Draw();
0249     c3_ts->cd(2); h_phi_gen_ts_eff[1]->Draw();
0250     c3_ts->cd(3); h_phi_gen_ts_eff[2]->Draw();
0251     c3_ts->cd(4); h_phi_gen_ts_eff[3]->Draw();
0252
0253     TCanvas *c3_rs = new TCanvas("c3_rs");
0254     c3_rs->Divide(2,2);
0255     c3_rs->cd(1); h_phi_gen_rs_eff[0]->Draw();
0256     c3_rs->cd(2); h_phi_gen_rs_eff[1]->Draw();
0257     c3_rs->cd(3); h_phi_gen_rs_eff[2]->Draw();
0258     c3_rs->cd(4); h_phi_gen_rs_eff[3]->Draw();
0259
0260     TCanvas *c4_ts = new TCanvas("c4_ts");
0261     c4_ts->Divide(2,2);
0262     c4_ts->cd(1); h_phi_gen_ts_eff[4]->Draw();
0263     c4_ts->cd(2); h_phi_gen_ts_eff[5]->Draw();
0264
0265     TCanvas *c4_rs = new TCanvas("c4_rs");
0266     c4_rs->Divide(2,2);
0267     c4_rs->cd(1); h_phi_gen_rs_eff[4]->Draw();
0268     c4_rs->cd(2); h_phi_gen_rs_eff[5]->Draw();
0269
0270     //Print plots to file

0271     c1_ts->Print("plots/track_eff_fixed_mom.pdf[");
0272     c1_ts->Print("plots/track_eff_fixed_mom.pdf");
0273     c1_rs->Print("plots/track_eff_fixed_mom.pdf");
0274     c2_ts->Print("plots/track_eff_fixed_mom.pdf");
0275     c2_rs->Print("plots/track_eff_fixed_mom.pdf");
0276     c3_ts->Print("plots/track_eff_fixed_mom.pdf");
0277     c3_rs->Print("plots/track_eff_fixed_mom.pdf");
0278     c4_ts->Print("plots/track_eff_fixed_mom.pdf");
0279     c4_rs->Print("plots/track_eff_fixed_mom.pdf");
0280     c4_rs->Print("plots/track_eff_fixed_mom.pdf]");
0281
0282 }