16 #include <TPolyLine.h>
20 #include <TLorentzVector.h>
21 #include <TEfficiency.h>
22 #include <TGraphAsymmErrors.h>
29 gROOT->SetStyle(
"Plain");
30 gStyle->SetOptStat(0);
32 gStyle->SetOptTitle(0);
46 int nlayers = n_maps_layer+n_intt_layer+n_tpc_layer;
50 double quality_cut = 3.0;
51 double dca2d_cut = 0.1;
52 double dcaz_cut = 0.1;
53 double gnhits_cut = 20;
55 double truth_hits_cut = 0.8;
69 TH1D *hnhits =
new TH1D(
"hnhits",
"hnhits",100,0,99);
70 TH2D *hpt_nfake =
new TH2D(
"hpt_nfake",
"hpt_nfake",200,0,hptmax, 73, -5, 68.0);
71 TH2D *hpt_nmissed_maps_layers =
new TH2D(
"hpt_nmissed_maps_layers",
"hpt_nmissed_maps_layers",200,0,hptmax, 53, -5, 48.0);
72 TH2D *hpt_nmissed_intt_layers =
new TH2D(
"hpt_nmissed_intt_layers",
"hpt_nmissed_intt_layers",200,0,hptmax, 53, -5, 48.0);
73 TH2D *hpt_nmissed_tpc_layers =
new TH2D(
"hpt_nmissed_tpc_layers",
"hpt_nmissed_tpc_layers",200,0,hptmax, 53, -5, 48.0);
74 TH2D *hcorr_nfake_nmaps =
new TH2D(
"hcorr_nfake_nmaps",
"hcorr_nfake_nmaps",40, -1, 5, 40, -1, 4);
75 TH1D *hzevt =
new TH1D(
"hzevt",
"hzevt",200,-35.0, 35.0);
76 TH1D *hzvtx_res =
new TH1D(
"hzvtx_res",
"hzvtx_res", 1000, -0.1, 0.1);
77 TH1D *hxvtx_res =
new TH1D(
"hxvtx_res",
"hxvtx_res", 1000, -0.04, 0.04);
78 TH1D *hyvtx_res =
new TH1D(
"hyvtx_res",
"hyvtx_res", 1000, -0.04, 0.04);
79 TH1D *hdcavtx_res =
new TH1D(
"hdcavtx_res",
"hdcavtx_res", 1000, -0.04, 0.04);
81 static const int NPTDCA = 3;
84 for(
int ipt=0;ipt<NPTDCA;ipt++)
87 sprintf(hname,
"hdca2d%i",ipt);
89 hdca2d[ipt] =
new TH1D(hname,hname,2000,-0.1,0.1);
90 hdca2d[ipt]->GetXaxis()->SetTitle(
"DCA (cm)");
91 hdca2d[ipt]->GetXaxis()->SetTitleSize(0.055);
92 hdca2d[ipt]->GetXaxis()->SetLabelSize(0.055);
94 hdca2d[ipt]->GetYaxis()->SetTitleSize(0.06);
95 hdca2d[ipt]->GetYaxis()->SetLabelSize(0.055);
98 TH1D *hdca2dsigma =
new TH1D(
"hdca2dsigma",
"hdca2dsigma",2000,-5.0,5.0);
99 hdca2dsigma->GetXaxis()->SetTitle(
"dca2d / dca2dsigma");
101 TH1D *hZdca =
new TH1D(
"hZdca",
"hZdca",2000,-0.5,0.5);
102 hZdca->GetXaxis()->SetTitle(
"Z DCA (cm)");
104 TH1D *hquality =
new TH1D(
"hquality",
"hquality",2000,0.0,20.0);
105 hquality->GetXaxis()->SetTitle(
"quality");
107 TH1D *hgeta =
new TH1D(
"hgeta",
"hgeta",100,-1.2,1.2);
108 TH1D *hreta =
new TH1D(
"hreta",
"hreta",100,-1.2,1.2);
110 static const int NVARBINS = 36;
111 double xbins[NVARBINS+1] = {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,
112 1.1, 1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,
113 2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,
114 4.5,5.0,5.5,6.0,7.0,8.0,10.0};
116 TH1D *hpt_truth =
new TH1D(
"hpt_truth",
"hpt_truth",500, 0.0, hptmax);
117 TH2D *hpt_compare =
new TH2D(
"hpt_compare",
"hpt_compare",500, 0.0, hptmax, 2000, 0.0, 2.0);
118 hpt_compare->GetXaxis()->SetTitle(
"p_{T}");
119 hpt_compare->GetYaxis()->SetTitle(
"#Delta p_{T}/p_{T}");
121 TH2D *hpt_dca2d =
new TH2D(
"hpt_dca2d",
"hpt_dca2d",500, 0.0, hptmax, 2000, -0.1, 0.1);
122 hpt_dca2d->GetXaxis()->SetTitle(
"p_{T}");
123 hpt_dca2d->GetYaxis()->SetTitle(
"dca2d");
125 TH2D *hpt_dcaZ =
new TH2D(
"hpt_dcaZ",
"hpt_dcaZ",500, 0.0, hptmax, 2000, -0.1, 0.1);
126 hpt_dcaZ->GetXaxis()->SetTitle(
"p_{T}");
127 hpt_dcaZ->GetYaxis()->SetTitle(
"dca_Z");
129 TH1D *hpt_hijing_truth =
new TH1D(
"hpt_hijing_truth",
"hpt_hijing_truth",500, 0.0, 10.0);
130 TH2D *hpt_hijing_compare =
new TH2D(
"hpt_hijing_compare",
"hpt_hijing_compare",500, 0.0, 10.0, 2000, 0.0, 2.0);
131 hpt_hijing_compare->GetXaxis()->SetTitle(
"p_{T}");
132 hpt_hijing_compare->GetYaxis()->SetTitle(
"#Delta p_{T}/p_{T}");
133 TH2D *hpt_hijing_dca2d =
new TH2D(
"hpt_hijing_dca2d",
"hpt_hijing_dca2d",500, 0.0, 10.0, 2000, -0.1, 0.1);
134 TH2D *hpt_hijing_dcaZ =
new TH2D(
"hpt_hijing_dcaZ",
"hpt_hijing_dcaZ",500, 0.0, 10.0, 2000, -0.1, 0.1);
136 TH1D *hptreco[NVARBINS];
137 for(
int ipt=0;ipt<NVARBINS;ipt++)
140 sprintf(hn,
"hptreco%i",ipt);
141 hptreco[ipt] =
new TH1D(hn, hn, NVARBINS, xbins);
144 TH2D *hg4ntrack =
new TH2D(
"hg4ntrack",
"hg4ntrack",200,0,3000.0, 200, 0., 2000);
161 for(
int i=0;
i<1200;
i++)
165 TChain* ntp_track =
new TChain(
"ntp_track",
"reco tracks");
166 TChain*
ntp_gtrack =
new TChain(
"ntp_gtrack",
"g4 tracks");
167 TChain* ntp_vertex =
new TChain(
"ntp_vertex",
"events");
168 TChain *ntp_cluster =
new TChain(
"ntp_cluster",
"clusters");
169 TChain *ntp_g4hit =
new TChain(
"ntp_g4hit",
"g4hits");
176 sprintf(name,
"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output/g4svtx_eval_%i.root_g4svtx_eval.root_primary_eval.root",
i);
185 bool checkfiles =
false;
190 TFile*
f=TFile::Open(name);
193 TNtuple*test_ntp_track=(TNtuple*)f->Get(
"ntp_track");
199 if(test_ntp_track->GetEntries() < mintracks)
205 cout <<
"Adding file number " <<
i <<
" with name " << name << endl;
207 ntp_vertex->Add(name);
208 ntp_track->Add(name);
209 ntp_gtrack->Add(name);
216 for(
int ig = 0;ig < ntp_gtrack->GetEntries();ig++)
218 int recoget = ntp_gtrack->GetEntry(ig);
224 for(
int ir = 0;ir < ntp_gtrack->GetEntries();ir++)
226 int recoget = ntp_track->GetEntry(ir);
227 double reta = asinh(
rpz/
rpt);
231 hg4ntrack->Fill(ng4tr, nrtr);
245 <<
" ntp_vertex entries: " << ntp_vertex->GetEntries()
246 <<
" ntp_gtrack entries: " << ntp_gtrack->GetEntries()
247 <<
" ntp_track entries: " << ntp_track->GetEntries()
257 int nev = ntp_vertex->GetEntries();
260 for(
int iev=0;iev<nev;iev++)
263 if(verbose > 0) cout <<
" iev = " << iev <<
" ng " << ng <<
" nr " << nr << endl;
265 int recoget = ntp_vertex->GetEntry(iev);
268 cout <<
"Failed to get ntp_vertex entry " << iev << endl;
293 cout <<
"Get event " << iev
294 <<
" with vertex x " <<
evx
295 <<
" vertex y " <<
evy <<
" vertex z " <<
evz
298 <<
" gtarcks start at " << ng <<
" ntracks start at " << nr
303 cout <<
" ALERT! event vertex is not correct: egvz = " <<
egvz <<
" evz = " <<
evz <<
" evx = " <<
evx <<
" evy = " <<
evy << endl;
309 if(fabs(
egvz) > 10.0)
311 cout <<
" Event vertex is outside 10 cm, reject this event - egvz = " <<
egvz << endl;
331 if(verbose > 0) cout <<
"Process truth tracks:" << endl;
333 int n_embed_gtrack = 0;
335 for(
int ig=ng;ig<ntp_gtrack->GetEntries();ig++)
337 int recoget = ntp_gtrack->GetEntry(ig);
341 cout <<
"Failed to get ntp_gtrack entry " << ig <<
" in ntp_gtrack" << endl;
367 if(verbose>0) cout <<
" -------- failed gnhits cut " << endl;
373 if(verbose>0) cout <<
" ------- failed nfromtruth/nhits cut " << endl;
378 if(verbose>0) cout <<
" ------- failed tgtrackid cut " << endl;
384 if(verbose > 0) cout <<
" -------- failed quality cut - rejected " << endl;
390 if(verbose>0) cout <<
" -------- skip because track failed z vertex cut with dca3dz = " <<
trdca3dz << endl;
396 if(verbose>0) cout <<
" -------- skip because failed dca2d cut, rdca2d = " <<
rdca2d << endl;
410 hpt_truth->Fill(tgpT);
422 hpt_hijing_truth->Fill(tgpT);
426 if(verbose>0) cout <<
"n_embed_gtrack = " << n_embed_gtrack << endl;
434 if(verbose > 0) cout <<
"Process reco tracks:" << endl;
436 int n_embed_rtrack = 0;
439 for(
int ir=nr;ir<ntp_track->GetEntries();ir++)
441 int recoget = ntp_track->GetEntry(ir);
444 if(verbose>0) cout <<
"Failed to get ntp_track entry " << ir << endl;
471 cout <<
" -------- skip because rgnhits too small " << endl;
477 if(verbose>0) cout <<
" ------- failed rnfromtruth/rnhits cut " << endl;
483 if(verbose>0) cout <<
" ------- failed rgtrackid < 0 cut " << endl;
488 if(verbose > 0) cout <<
" -------- failed quality cut - rejected " << endl;
494 if(verbose>0) cout <<
" -------- skip because track " << ir <<
" failed z vertex cut with rvz = " <<
rpcaz <<
" gvz = " <<
rvz << endl;
500 if(verbose>0) cout <<
" -------- skip because failed dca2d cut, rdca2d = " <<
rdca2d << endl;
520 hZdca->Fill(
rpcaz - rvz);
522 if(verbose > 0) cout <<
" accepted: rgembed = " <<
rgembed <<
" rgpt = " <<
rgpt << endl;
526 if(verbose > 0) cout <<
" rdca2d = " <<
rdca2d <<
" rdcaZ = " << rdcaZ << endl;
528 double reta = asinh(
rpz/
rpt);
534 if(
rnmaps == n_maps_layer)
543 hpt_nfake->Fill(
rgpt, nfake);
546 hpt_nmissed_maps_layers->Fill(
rgpt,nmissed_maps_layers);
549 hpt_nmissed_intt_layers->Fill(
rgpt,nmissed_intt_layers);
553 hpt_nmissed_tpc_layers->Fill(
rgpt,nmissed_tpc_layers);
555 hcorr_nfake_nmaps->Fill(nfake,nmissed_maps_layers);
562 if(
rnmaps == n_maps_layer)
565 hpt_hijing_dcaZ->Fill(
rgpt,rdcaZ);
568 for(
int ipt=0;ipt<NVARBINS-1;ipt++)
570 if(
rpt > xbins[ipt] &&
rpt < xbins[ipt+1])
571 hptreco[ipt]->Fill(
rpt);
581 if(verbose > 0) cout <<
" Done with loop: n_embed_gtrack = " << n_embed_gtrack <<
" n_embed_rtrack = " << n_embed_rtrack <<
" naccept = " << naccept << endl;
587 if(verbose > 0) cout <<
" embedded tracks this event = " << n_embed_gtrack <<
" accepted tracks this event = " << naccept << endl;
601 TFile *
fout =
new TFile(
"root_files/purity_out.root",
"recreate");
612 hpt_hijing_truth->Write();
615 hpt_compare->Write();
619 hpt_nmissed_maps_layers->Write();
620 hpt_nmissed_intt_layers->Write();
621 hpt_nmissed_tpc_layers->Write();
622 hcorr_nfake_nmaps->Write();
625 hpt_hijing_compare->Write();
626 hpt_hijing_dca2d->Write();
627 hpt_hijing_dcaZ->Write();
630 for(
int ipt=0;ipt<NVARBINS;ipt++)
632 hptreco[ipt]->Write();
642 hdcavtx_res->Write();