Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
purity.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file purity.C
1 #include <TChain.h>
2 #include <TNtuple.h>
3 #include <TCanvas.h>
4 #include <TH1.h>
5 #include <TH1F.h>
6 #include <TH2D.h>
7 #include <TF1.h>
8 #include <TFile.h>
9 #include <TGraph.h>
10 #include <TStyle.h>
11 #include <TROOT.h>
12 #include <TLegend.h>
13 #include <TLine.h>
14 #include <TLatex.h>
15 #include <TRandom1.h>
16 #include <TPolyLine.h>
17 #include <iostream>
18 #include <fstream>
19 #include <TMath.h>
20 #include <TLorentzVector.h>
21 #include <TEfficiency.h>
22 #include <TGraphAsymmErrors.h>
23 
24 using namespace std;
25 #define use_events
26 
27 void purity()
28 {
29  gROOT->SetStyle("Plain");
30  gStyle->SetOptStat(0);
31  gStyle->SetOptFit(0);
32  gStyle->SetOptTitle(0);
33 
34  int verbose = -1;
35 
36  // Setup parameters
37  //=================
38  int embed_flag = 2; // set to match value used when events were generated
39 
40  // these should match the values used when the files were created
41  //int n_maps_layer = 3;
42  int n_maps_layer = 0;
43  //int n_intt_layer = 4;
44  int n_intt_layer = 0;
45  int n_tpc_layer = 48; // 16 inner
46  int nlayers = n_maps_layer+n_intt_layer+n_tpc_layer; // maximum number of tracking layers for tpc+intt+maps
47 
48  // track cuts
49  //static const int nmissed = 20; // maximum number of missed layers to allow for a track
50  double quality_cut = 3.0;
51  double dca2d_cut = 0.1; // cm
52  double dcaz_cut = 0.1; // cm
53  double gnhits_cut = 20;
54  //double gnhits_cut = 10;
55  double truth_hits_cut = 0.8;
56 
57  // other parameters
58  double ptmax = 12.2; // used for Hijing tracks only
59  //double hptmax = 10.0;
60  double hptmax = 50.0;
61 
62  //===============================
63 
64 
65  //=======================================
66  // define histograms
67 
68 
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);
80 
81  static const int NPTDCA = 3;
82 
83  TH1D *hdca2d[NPTDCA];
84  for(int ipt=0;ipt<NPTDCA;ipt++)
85  {
86  char hname[500];
87  sprintf(hname,"hdca2d%i",ipt);
88 
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);
93 
94  hdca2d[ipt]->GetYaxis()->SetTitleSize(0.06);
95  hdca2d[ipt]->GetYaxis()->SetLabelSize(0.055);
96  }
97 
98  TH1D *hdca2dsigma = new TH1D("hdca2dsigma","hdca2dsigma",2000,-5.0,5.0);
99  hdca2dsigma->GetXaxis()->SetTitle("dca2d / dca2dsigma");
100 
101  TH1D *hZdca = new TH1D("hZdca","hZdca",2000,-0.5,0.5);
102  hZdca->GetXaxis()->SetTitle("Z DCA (cm)");
103 
104  TH1D *hquality = new TH1D("hquality","hquality",2000,0.0,20.0);
105  hquality->GetXaxis()->SetTitle("quality");
106 
107  TH1D *hgeta = new TH1D("hgeta","hgeta",100,-1.2,1.2);
108  TH1D *hreta = new TH1D("hreta","hreta",100,-1.2,1.2);
109 
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};
115 
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}");
120 
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");
124 
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");
128 
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);
135 
136  TH1D *hptreco[NVARBINS];
137  for(int ipt=0;ipt<NVARBINS;ipt++)
138  {
139  char hn[1000];
140  sprintf(hn,"hptreco%i",ipt);
141  hptreco[ipt] = new TH1D(hn, hn, NVARBINS, xbins);
142  }
143 
144  TH2D *hg4ntrack = new TH2D("hg4ntrack","hg4ntrack",200,0,3000.0, 200, 0., 2000);
145 
146  // end of histogram definitions
147  //===================================================
148 
149  //============================================================
150  // Loop over events
151  // Loop over reco'd tracks
152  // Make quality cut
153  // Fill dca2d histos
154  // Make loose dca2d cut (usually < 1 mm)
155  // Fill Z dca histo
156  // Drop tracks outside 1 mm in Z dca from evt vertex
157  // Add this reco'd track 2D histo of (rpT-pT) vs pT
158  //============================================================
159 
160  // The condor job output files
161  for(int i=0;i<1200;i++)
162  {
163  // Open the evaluator output file
164 
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");
170 
171  // This include file contains the definitions of the ntuple variables, and the chain definitions
172 #include "ntuple_variables.C"
173 
174  char name[1000];
175  // latest files
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);
177 
178  // for CD1 plots with no MVTX
179  //sprintf(name,"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output_noINTT_80ns_100pions_clusterizer_fixed/g4svtx_eval_%i.root_g4svtx_eval.root_primary_eval.root",i);
180  //sprintf(name,"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output_noINTT_80ns_100pions_clusterizer_fixed/g4svtx_eval_%i.root_g4svtx_eval.root",i);
181 
182  //sprintf(name,"/sphenix/user/frawley/fresh_mar8_testing/macros/macros/g4simulations/eval_output_ups1s_100pions_80ns/g4svtx_eval_%i.root_g4svtx_eval.root_primary_eval.root",i);
183 
184 
185  bool checkfiles = false;
186  int mintracks = 10;
187  if(checkfiles)
188  {
189  // some QA on the input files
190  TFile*f=TFile::Open(name);
191  if(!f)
192  continue;
193  TNtuple*test_ntp_track=(TNtuple*)f->Get("ntp_track");
194  if(!test_ntp_track)
195  {
196  f->Close();
197  continue;
198  }
199  if(test_ntp_track->GetEntries() < mintracks)
200  continue;
201 
202  f->Close();
203  }
204 
205  cout << "Adding file number " << i << " with name " << name << endl;
206 
207  ntp_vertex->Add(name);
208  ntp_track->Add(name);
209  ntp_gtrack->Add(name);
210 
211 
212 
213  // Use the number of g4 tracks (with some cuts) as a measure of centrality
214  // Examination of the output file shows that the 10% most central events have 1700 or more g4 tracks with these cuts
215  double ng4tr = 0;
216  for(int ig = 0;ig < ntp_gtrack->GetEntries();ig++)
217  {
218  int recoget = ntp_gtrack->GetEntry(ig);
219  if(tembed == 0 && tprimary == 1 && fabs(teta) < 1.0 && sqrt(tpx*tpx+tpy*tpy) > 0.2)
220  ng4tr++;
221  }
222 
223  double nrtr = 0;
224  for(int ir = 0;ir < ntp_gtrack->GetEntries();ir++)
225  {
226  int recoget = ntp_track->GetEntry(ir);
227  double reta = asinh(rpz/rpt);
228  if(rgembed == 0 && rprimary == 1 && fabs(reta) < 1.0 && rpt > 0.2)
229  nrtr++;
230  }
231  hg4ntrack->Fill(ng4tr, nrtr);
232 
233  //if(ng4tr < 1250) // keep 0-20%
234  //if(ng4tr > 700) // keep 60-100%
235  if(ng4tr > 10000) // keep 0-100%
236  {
237  delete ntp_gtrack;
238  delete ntp_track;
239  delete ntp_cluster;
240  delete ntp_vertex;
241  continue;
242  }
243  if(verbose> -1)
244  cout
245  << " ntp_vertex entries: " << ntp_vertex->GetEntries()
246  << " ntp_gtrack entries: " << ntp_gtrack->GetEntries()
247  << " ntp_track entries: " << ntp_track->GetEntries()
248  << endl;
249 
250  //==================================
251  // Begin event loop
252  //==================================
253 
254  // These keep track of the starting position in the track ntuples for each event
255  int nr = 0;
256  int ng = 0;
257  int nev = ntp_vertex->GetEntries();
258  //int nev = 1;
259 
260  for(int iev=0;iev<nev;iev++)
261  {
262 
263  if(verbose > 0) cout << " iev = " << iev << " ng " << ng << " nr " << nr << endl;
264 
265  int recoget = ntp_vertex->GetEntry(iev);
266  if(!recoget)
267  {
268  cout << "Failed to get ntp_vertex entry " << iev << endl;
269  continue;
270  }
271 
272  /*
273  if(egvt != 0)
274  {
275  // this is a pileup event, we want only the triggered event, so skip it
276  // we need to advance the pointers to the track ntuples
277  if(isnan(ntracks))
278  ntracks = 0;
279  nr += ntracks;
280  ng += ngtracks;
281 
282  if(verbose>0)
283  cout << " skip pileup event with ngtracks = " << ngtracks << " and ntracks = " << ntracks
284  << " advance ng to " << ng << " and nr to " << nr
285  << endl;
286 
287  continue;
288  }
289  */
290 
291  if(iev%1 == 0)
292  if(verbose>0)
293  cout << "Get event " << iev
294  << " with vertex x " << evx
295  << " vertex y " << evy << " vertex z " << evz
296  << " egvz " << egvz
297  << " ngtracks " << ngtracks << " ntracks " << ntracks
298  << " gtarcks start at " << ng << " ntracks start at " << nr
299  << endl;
300 
301  if(fabs(evz - egvz) > 0.1)
302  {
303  cout << " ALERT! event vertex is not correct: egvz = " << egvz << " evz = " << evz << " evx = " << evx << " evy = " << evy << endl;
304  continue;
305  }
306 
307 
308  // do not consider events that are not within the full acceptance
309  if(fabs(egvz) > 10.0)
310  {
311  cout << " Event vertex is outside 10 cm, reject this event - egvz = " << egvz << endl;
312  if(isnan(ntracks))
313  ntracks = 0;
314  nr += ntracks;
315  ng += ngtracks;
316  continue;
317  }
318 
319  hzevt->Fill(evz);
320  hzvtx_res->Fill(evz-egvz);
321  hxvtx_res->Fill(evx-egvx);
322  hyvtx_res->Fill(evy-egvy);
323  hdcavtx_res->Fill(sqrt(evx*evx+evy*evy) - sqrt(egvx*egvx+egvy*egvy));
324 
325  //====================================================
326  // ntp_gtracks
327  //====================================================
328  // ngtracks is defined in ntuple_variables.C and is the number of g4 tracks
329  // ntracks is defined in ntuple_variables.C and is the number of reco'd tracks
330 
331  if(verbose > 0) cout << "Process truth tracks:" << endl;
332 
333  int n_embed_gtrack = 0;
334  //for(int ig=ng;ig<ng+ngtracks;ig++)
335  for(int ig=ng;ig<ntp_gtrack->GetEntries();ig++)
336  {
337  int recoget = ntp_gtrack->GetEntry(ig);
338  if(recoget == 0)
339  {
340  //if(verbose) cout << "Failed to get ntp_gtrack entry " << ig << " in ntp_gtrack" << endl;
341  cout << "Failed to get ntp_gtrack entry " << ig << " in ntp_gtrack" << endl;
342  break;
343  }
344 
345  //if(tembed != embed_flag)
346  // continue;
347 
348  /*
349  // if the scan_for_embedded flag is set to true, the number of reco tracks recorded in ntp_vertex will be all reco'd tracks, but only embedded tracks will be present
350  // check for change of event number to detect this
351  if(tevent != iev)
352  {
353  if(verbose > 0) cout << "Change of event number from " << iev << " to " << tevent << " go to next event" << endl;
354  ngtracks = ig - ng;
355  break;
356  }
357 
358  if(verbose > 0) cout << " ig = " << ig << " tevent " << tevent << " tgtrackid " << tgtrackid
359  << " ttrackid " << ttrackid << " tgnhits " << tgnhits << " tnhits " << tnhits << " tnfromtruth " << tnfromtruth << " tembed " << tembed << endl;
360  */
361 
362  // ntp_gtrack track cuts
363  //===================
364  // skip tracks that do not pass through enough layers (judged using truth track value of nhits)
365  if(tgnhits <= gnhits_cut)
366  {
367  if(verbose>0) cout << " -------- failed gnhits cut " << endl;
368  continue;
369  }
370 
371  if(tnfromtruth / tnhits <= truth_hits_cut)
372  {
373  if(verbose>0) cout << " ------- failed nfromtruth/nhits cut " << endl;
374  continue;
375  }
376  if(tgtrackid < 0)
377  {
378  if(verbose>0) cout << " ------- failed tgtrackid cut " << endl;
379  continue;
380  }
381 
382  if(tquality > quality_cut)
383  {
384  if(verbose > 0) cout << " -------- failed quality cut - rejected " << endl;
385  continue;
386  }
387 
388  if(fabs(trdca3dz) > dcaz_cut)
389  {
390  if(verbose>0) cout << " -------- skip because track failed z vertex cut with dca3dz = " << trdca3dz << endl;
391  continue;
392  }
393 
394  if(trdca2d > dca2d_cut)
395  {
396  if(verbose>0) cout << " -------- skip because failed dca2d cut, rdca2d = " << rdca2d << endl;
397  continue;
398  }
399  // get the truth pT
400  double tgpT = sqrt(tpx*tpx+tpy*tpy);
401 
402  if(tembed == embed_flag) // take only embedded pions
403  {
404  double geta = asinh(tpz/sqrt(tpx*tpx+tpy*tpy));
405  hgeta->Fill(geta); // optional cut
406  if(fabs(geta) < 1.0)
407  {
408  n_embed_gtrack++;
409  //cout << "fill htp_truth with tgpT = " << tgpT << endl;
410  hpt_truth->Fill(tgpT);
411 
412  // Capture hpt_compare from ntp_gtrack
413  hpt_compare->Fill(tpt,trpt/tpt);
414  hpt_dca2d->Fill(tpt, trdca2d);
415  hpt_dcaZ->Fill(tpt, trdca3dz);
416  }
417  }
418 
419  // record embedded pions and Hijing tracks separately
420  if(tembed == 0)
421  {
422  hpt_hijing_truth->Fill(tgpT);
423  }
424 
425  } // end loop over truth tracks
426  if(verbose>0) cout << "n_embed_gtrack = " << n_embed_gtrack << endl;
427  //===================================================
428 
429 
430  //====================================================
431  // ntp_tracks
432  //====================================================
433 
434  if(verbose > 0) cout << "Process reco tracks:" << endl;
435 
436  int n_embed_rtrack = 0;
437  int naccept = 0;
438  //for(int ir=nr;ir<nr+ntracks;ir++)
439  for(int ir=nr;ir<ntp_track->GetEntries();ir++)
440  {
441  int recoget = ntp_track->GetEntry(ir);
442  if(!recoget)
443  {
444  if(verbose>0) cout << "Failed to get ntp_track entry " << ir << endl;
445  break;
446  }
447 
448  if(rgembed != embed_flag)
449  continue;
450 
451  /*
452  // if the scan_for_embedded flag is set to true, the number of reco tracks recorded in ntp_vertex will be all reco'd tracks, but only embedded tracks will be present
453  // check for change of event number to detect this
454  if(revent != iev)
455  {
456  if(verbose > 0) cout << "Change of event number from " << iev << " to " << revent << " go to next event" << endl;
457  ntracks = ir - nr;
458  break;
459  }
460 
461  if(verbose > 0) cout << " ir = " << ir << " rtrackid " << rtrackid << " revent " << revent
462  << " rquality " << rquality << " rgtrackid = " << rgtrackid << " rgnhits " << rgnhits
463  << " rnhits " << rnhits << " rnfromtruth " << rnfromtruth << " rgembed " << rgembed << endl;
464  */
465 
466  // ntp_track track cuts
467  //================================
468  if(rgnhits <= gnhits_cut)
469  {
470  //skip tracks that do not pass through enough layers in truth, same cut made on truth tracks earlier
471  cout << " -------- skip because rgnhits too small " << endl;
472  continue;
473  }
474 
475  if(rnfromtruth / rnhits <= truth_hits_cut)
476  {
477  if(verbose>0) cout << " ------- failed rnfromtruth/rnhits cut " << endl;
478  continue;
479  }
480 
481  if(rgtrackid < 0)
482  {
483  if(verbose>0) cout << " ------- failed rgtrackid < 0 cut " << endl;
484  continue;
485  }
486  if(rquality > quality_cut)
487  {
488  if(verbose > 0) cout << " -------- failed quality cut - rejected " << endl;
489  continue;
490  }
491 
492  if(fabs(rpcaz - rvz) > dcaz_cut)
493  {
494  if(verbose>0) cout << " -------- skip because track " << ir << " failed z vertex cut with rvz = " << rpcaz << " gvz = " << rvz << endl;
495  continue;
496  }
497 
498  if(rdca2d > dca2d_cut)
499  {
500  if(verbose>0) cout << " -------- skip because failed dca2d cut, rdca2d = " << rdca2d << endl;
501  continue;
502  }
503 
504  /*
505  if(rnmaps != n_maps_layer)
506  {
507  if(verbose>0) cout << "skip track " << ir << " because nmaps = " << rnmaps << endl;
508  continue;
509  }
510  */
511 
512  //=============================
513 
514  double rdcaZ = rpcaz - rvz;
515 
516  if(rgembed == embed_flag) // take only embedded pions
517  {
518  // get the quality and Zdca histos before vertex cuts
519  hquality->Fill(rquality);
520  hZdca->Fill(rpcaz - rvz);
521 
522  if(verbose > 0) cout << " accepted: rgembed = " << rgembed << " rgpt = " << rgpt << endl;
523 
524  naccept++;
525 
526  if(verbose > 0) cout << " rdca2d = " << rdca2d << " rdcaZ = " << rdcaZ << endl;
527 
528  double reta = asinh(rpz/rpt);
529  hreta->Fill(reta);
530  if(fabs(reta) < 1.0) // optional cut
531  {
532  n_embed_rtrack++;
533  //hpt_compare->Fill(rgpt,rpt/rgpt);
534  if(rnmaps == n_maps_layer) // require all maps layers
535  {
536  //hpt_dca2d->Fill(rgpt, rdca2d);
537  //hpt_dcaZ->Fill(rgpt, rdcaZ);
538  }
539  }
540  hnhits->Fill(rnhits);
541 
542  double nfake = rnhits - rnfromtruth;
543  hpt_nfake->Fill(rgpt, nfake);
544 
545  double nmissed_maps_layers = rgnlmaps - rnlmaps;
546  hpt_nmissed_maps_layers->Fill(rgpt,nmissed_maps_layers);
547 
548  double nmissed_intt_layers = rgnlintt - rnlintt;
549  hpt_nmissed_intt_layers->Fill(rgpt,nmissed_intt_layers);
550 
551 
552  double nmissed_tpc_layers = rgnltpc - rnltpc;
553  hpt_nmissed_tpc_layers->Fill(rgpt,nmissed_tpc_layers);
554 
555  hcorr_nfake_nmaps->Fill(nfake,nmissed_maps_layers);
556  }
557 
558  if(rgembed == 1)
559  {
560  // for the non-embedded tracks from Hijing, we want to be able to do pt_sigma cuts later
561  hpt_hijing_compare->Fill(rgpt,rpt/rgpt);
562  if(rnmaps == n_maps_layer) // require all maps layers
563  {
564  hpt_hijing_dca2d->Fill(rgpt,rdca2d);
565  hpt_hijing_dcaZ->Fill(rgpt,rdcaZ);
566  }
567 
568  for(int ipt=0;ipt<NVARBINS-1;ipt++)
569  {
570  if(rpt > xbins[ipt] && rpt < xbins[ipt+1])
571  hptreco[ipt]->Fill(rpt);
572  }
573 
574  // Add to the 3 panel dca2d histos
575  if(rgpt > 0.5 && rgpt <= 1.0) hdca2d[0]->Fill(rdca2d);
576  if(rgpt > 1.0 && rgpt <= 2.0) hdca2d[1]->Fill(rdca2d);
577  if(rgpt > 2.0) hdca2d[2]->Fill(rdca2d);
578  }
579 
580  } // end loop over reco'd tracks
581  if(verbose > 0) cout << " Done with loop: n_embed_gtrack = " << n_embed_gtrack << " n_embed_rtrack = " << n_embed_rtrack << " naccept = " << naccept << endl;
582  //====================================================
583 
584  // set the starting positions in the track ntuples for the next event
585  nr += ntracks;
586  ng += ngtracks;
587  if(verbose > 0) cout << " embedded tracks this event = " << n_embed_gtrack << " accepted tracks this event = " << naccept << endl;
588  } // end loop over events
589 
590 
591  delete ntp_gtrack;
592  delete ntp_track;
593  delete ntp_cluster;
594  delete ntp_vertex;
595 
596  } // end loop over files
597 
598  //==============
599  // Output the results
600 
601  TFile *fout = new TFile("root_files/purity_out.root","recreate");
602 
603  hnhits->Write();
604 
605  // Three dca2d histos made from Hijing tracks
606  hdca2d[0]->Write();
607  hdca2d[1]->Write();
608  hdca2d[2]->Write();
609 
610  // truth pT distributions for embedded and non-embedded particles
611  hpt_truth->Write();
612  hpt_hijing_truth->Write();
613 
614  // 2D histos made with embedded pions
615  hpt_compare->Write();
616  hpt_dca2d->Write();
617  hpt_dcaZ->Write();
618  hpt_nfake->Write();
619  hpt_nmissed_maps_layers->Write();
620  hpt_nmissed_intt_layers->Write();
621  hpt_nmissed_tpc_layers->Write();
622  hcorr_nfake_nmaps->Write();
623 
624  // 2d histo made with hijing tracks only
625  hpt_hijing_compare->Write();
626  hpt_hijing_dca2d->Write();
627  hpt_hijing_dcaZ->Write();
628 
629  // efficiency plot made with Hijing tracks
630  for(int ipt=0;ipt<NVARBINS;ipt++)
631  {
632  hptreco[ipt]->Write();
633  }
634 
635  hquality->Write();
636  hZdca->Write();
637 
638  hzevt->Write();
639  hzvtx_res->Write();
640  hxvtx_res->Write();
641  hyvtx_res->Write();
642  hdcavtx_res->Write();
643 
644  hreta->Write();
645  hgeta->Write();
646  hg4ntrack->Write();
647 
648  fout->Close();
649 
650 
651 }
652