Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
look_purity.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file look_purity.C
1 #include <TH2D.h>
2 #include <TFile.h>
3 #include <TF1.h>
4 #include <TGraph.h>
5 #include <TCanvas.h>
6 #include <TColor.h>
7 #include <TROOT.h>
8 #include <TStyle.h>
9 #include <TLegend.h>
10 #include <TLatex.h>
11 #include <TLine.h>
12 #include <TEfficiency.h>
13 #include <TGraphAsymmErrors.h>
14 
15 #include "sPhenixStyle.C"
16 
17 void look_purity()
18 {
19  SetsPhenixStyle();
20  /*
21  gROOT->SetStyle("Plain");
22  gStyle->SetOptStat(0);
23  gStyle->SetOptFit(0);
24  gStyle->SetOptTitle(0);
25  */
26 
27  //gStyle->SetOptStat(1);
28 
29  // determines if plots are labeled "HIJING"
30  bool hijing = false;
31 
32  //int n_maps_layer = 3;
33  int n_maps_layer = 0;
34  //int n_intt_layer = 4;
35  int n_intt_layer = 0;
36  int n_tpc_layer = 48;
37 
38  double ptmax = 40.0;
39  double hptmax = 40.0;
40  //double ptstep = 0.25;
41  //double ptstep = 1.0;
42  //double ptstep = 2.0;
43  //double ptstep = 0.5;
44  double ptstep = 0.4;
45 
46  // set to false only to generate pT resolution plots without fits
47  // BEWARE: false means that the 4 sigma cuts are meaningless - they are not done with fitted parameters
48  bool pt_resolution_fit = true;
49  //bool pt_resolution_fit = false;
50 
51  TFile *fin = new TFile("root_files/purity_out.root");
52  //TFile *fin = new TFile("root_files/purity_out_100pions_noINTT_primvert.root");
53 
54 
55  if(!fin)
56  {
57  cout << "Failed to find input file" << endl;
58  exit(1);
59  }
60 
61  TCanvas *c1 = new TCanvas("c1","c1",5,5,900,500);
62  c1->Divide(3,1);
63 
64  //==========================
65  // 2D histo for embedded pions
66  TH2D *hpt_dca2d = 0;
67  fin->GetObject("hpt_dca2d",hpt_dca2d);
68  c1->cd(1);
69  gPad->SetLeftMargin(0.12);
70  hpt_dca2d->GetYaxis()->SetTitleOffset(1.5);
71  hpt_dca2d->SetMarkerStyle(20);
72  hpt_dca2d->SetMarkerSize(0.3);
73  hpt_dca2d->Draw();
74 
75  // 2D histo for embedded pions
76  TH2D *hpt_dcaZ = 0;
77  fin->GetObject("hpt_dcaZ",hpt_dcaZ);
78  c1->cd(2);
79  gPad->SetLeftMargin(0.12);
80  hpt_dcaZ->SetMarkerStyle(20);
81  hpt_dcaZ->SetMarkerSize(0.3);
82  hpt_dcaZ->GetYaxis()->SetTitleOffset(1.5);
83  hpt_dcaZ->Draw();
84 
85  TH2D *hpt_compare = 0;
86  fin->GetObject("hpt_compare",hpt_compare);
87  c1->cd(3);
88  gPad->SetLeftMargin(0.12);
89  hpt_compare->GetYaxis()->SetTitleOffset(1.5);
90  hpt_compare->SetMarkerStyle(20);
91  hpt_compare->SetMarkerSize(0.3);
92  hpt_compare->Draw();
93 
94  //========================================
95  // extract DCA resolution vs pT from 2D histo hpt_dca2d
96  // by fitting NPT pT slices
97 
98  TCanvas *c2 = new TCanvas("c2","c2",20,20,800,600);
99 
100  int NPT = (int) (ptmax/ptstep);
101 
102  double pT[200];
103  double dca2d[200];
104  for(int i = 0;i<NPT;i++)
105  {
106  // divide pT range into bins of width ptstep GeV/c at 0.25, 0.75, 1.25, .....
107  // bins will be (e.g.)from 0-0.5, 0.5-1.0, 1.0-1.5, .....
108  double ptlo = (double) i * ptstep + 0.0;
109  double pthi = (double) i * ptstep + ptstep;
110 
111  int binlo = hpt_dca2d->GetXaxis()->FindBin(ptlo);
112  int binhi = hpt_dca2d->GetXaxis()->FindBin(pthi);
113 
114  std::cout << "ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << std::endl;
115 
116  TH1D *h = new TH1D("h","dca2d resolution",2000, -0.1, 0.1);
117  hpt_dca2d->ProjectionY("h",binlo,binhi);
118  h->GetXaxis()->SetTitle("p_{T} (GeV/c)");
119  h->GetXaxis()->SetTitle("#Delta dca2d (cm)");
120  h->GetXaxis()->SetTitleOffset(1.0);
121  //if(i<8) h->Rebin(4);
122  h->DrawCopy();
123 
124  double mean = h->GetMean();
125  double sigma = h->GetRMS();
126  double yield = h->Integral();
127 
128  pT[i] = (ptlo + pthi) / 2.0;
129 
130  double low = -0.01, high=0.01;
131  if(n_maps_layer == 0)
132  {
133  low = 3.0 * low;
134  high = 3.0 * high;
135  }
136  if(pT[i] < 6.0)
137  {
138  low = 3.0*low;
139  high = 3.0*high;
140  }
141 
142  TF1 *f = new TF1("f","gaus",low, high);
143  f->SetParameter(1, yield/100.0);
144  f->SetParameter(2, 0.0);
145  f->SetParameter(3,0.002);
146  h->Fit(f,"R");
147 
148 
149 
150  dca2d[i] = f->GetParameter(2);
151  cout << " pT " << pT[i] << " dca2d " << dca2d[i] << " counts " << h->Integral() << " hist mean " << h->GetMean() << " hist RMS " << h->GetRMS() << endl;
152  }
153 
154  //============================================
155  // plot the dca2d resolution extracted above for embedded pions
156 
157  TCanvas *c3 = new TCanvas("c3","c3",100,100,800,600);
158  TGraph *grdca2d = new TGraph(NPT,pT,dca2d);
159  grdca2d->SetMarkerStyle(20);
160  grdca2d->SetMarkerSize(1.2);
161  grdca2d->SetMarkerColor(kRed);
162  grdca2d->SetName("dca2d_resolution");
163  grdca2d->SetTitle("dca2d resolution");
164 
165  TH1D *hdummy = new TH1D("hdummy","#Delta dca2d vs p_{T}",100,0.0,hptmax);
166  hdummy->SetMinimum(0);
167  if(n_maps_layer == 0)
168  hdummy->SetMaximum(0.015);
169  else
170  hdummy->SetMaximum(0.0051);
171  hdummy->GetXaxis()->SetTitle("p_{T} (GeV/c)");
172  hdummy->GetYaxis()->SetTitle("DCA(r#phi) resolution (cm)");
173  hdummy->GetYaxis()->SetTitleOffset(1.5);
174  gPad->SetLeftMargin(0.15);
175  hdummy->Draw();
176  grdca2d->Draw("p");
177 
178  char dcalab[500];
179  if(hijing)
180  //sprintf(dcalab, "sPHENIX HIJING b=0-4 fm");
181  sprintf(dcalab, "sPHENIX HIJING b=0-12 fm, 100 kHz pileup");
182  else
183  sprintf(dcalab, "sPHENIX 100 pions");
184  TLegend *ldca = new TLegend(0.3, 0.75, 1.05, 0.90,dcalab,"NDC");
185  ldca->SetBorderSize(0);
186  ldca->SetFillColor(0);
187  ldca->SetFillStyle(0);
188  char lstr1[500];
189  //sprintf(lstr1,"MVTX+INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
190  sprintf(lstr1,"INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
191  ldca->AddEntry(grdca2d, lstr1, "p");
192  ldca->Draw();
193 
194  //===================================
195 
196  //========================================
197  // extract DCA Z resolution vs pT from 2D histo hpt_dcaZ
198  // by fitting NPT pT slices
199 
200  TCanvas *c2a = new TCanvas("c2a","c2a",20,20,800,600);
201 
202  double dcaZ[200];
203  for(int i = 0;i<NPT;i++)
204  {
205  // divide pT range into bins of width 0.5 GeV/c at 0.25, 0.75, 1.25, .....
206  // bins will be from 0-0.5, 0.5-1.0, 1.0-1.5, .....
207  double ptlo = (double) i * ptstep + 0.0;
208  double pthi = (double) i * ptstep + ptstep;
209 
210  int binlo = hpt_dcaZ->GetXaxis()->FindBin(ptlo);
211  int binhi = hpt_dcaZ->GetXaxis()->FindBin(pthi);
212 
213  std::cout << "ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << std::endl;
214  TH1D *h = new TH1D("h","DCA (Z) resolution (cm)",2000, -0.1, 0.1);
215  hpt_dcaZ->ProjectionY("h",binlo,binhi);
216  h->GetXaxis()->SetTitle("p_{T} (GeV/c)");
217  h->GetYaxis()->SetTitle("DCA(Z) resolution (cm)");
218  h->GetYaxis()->SetTitleOffset(1.0);
219  //if(i<8) h->Rebin(4);
220  h->DrawCopy();
221 
222  double mean = h->GetMean();
223  double sigma = h->GetRMS();
224  double yield = h->Integral();
225 
226  pT[i] = (ptlo + pthi) / 2.0;
227 
228  double low = -0.01, high=0.01;
229  if(n_maps_layer == 0)
230  {
231  low = 3.0 * low;
232  high = 3.0 * high;
233  }
234 
235  if(pT[i] < 6.0)
236  {
237  low = 3.0*low;
238  high = 3.0*high;
239  }
240 
241  TF1 *f = new TF1("f","gaus",low, high);
242  f->SetParameter(1, yield/100.0);
243  f->SetParameter(2, 0.0);
244  f->SetParameter(3,0.002);
245  h->Fit(f,"R");
246 
247  dcaZ[i] = f->GetParameter(2);
248  cout << " pT " << pT[i] << " dcaZ " << dcaZ[i] << " counts " << h->Integral() << " hist mean " << h->GetMean() << " hist RMS " << h->GetRMS() << endl;
249  }
250 
251  //============================================
252  // plot the dcaZ resolution extracted above for embedded pions
253 
254  TCanvas *c3a = new TCanvas("c3a","c3a",100,100,800,600);
255  TGraph *grdcaZ = new TGraph(NPT,pT,dcaZ);
256  grdcaZ->SetMarkerStyle(20);
257  grdcaZ->SetMarkerSize(1.2);
258  grdcaZ->SetMarkerColor(kRed);
259  grdcaZ->SetName("dcaZ_resolution");
260  grdcaZ->SetTitle("dcaZ resolution");
261 
262  TH1D *hdummya = new TH1D("hdummya","#Delta dcaZ vs p_{T}",100,0.0,hptmax);
263  hdummya->SetMinimum(0);
264 
265  if(n_maps_layer == 0)
266  hdummya->SetMaximum(0.051);
267  else
268  hdummya->SetMaximum(0.0051);
269  hdummya->GetXaxis()->SetTitle("p_{T} (GeV/c)");
270  hdummya->GetYaxis()->SetTitle("DCA(Z) resolution (cm)");
271  hdummya->GetYaxis()->SetTitleOffset(1.5);
272  gPad->SetLeftMargin(0.15);
273  hdummya->Draw();
274  grdcaZ->Draw("p");
275 
276  TLegend *ldcaZ = new TLegend(0.3, 0.75, 1.05, 0.90, dcalab,"NDC");
277  ldcaZ->SetBorderSize(0);
278  ldcaZ->SetFillColor(0);
279  ldcaZ->SetFillStyle(0);
280  //char lstr1[500];
281  //sprintf(lstr1,"MVTX+INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
282  sprintf(lstr1,"INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
283  ldcaZ->AddEntry(grdcaZ, lstr1, "p");
284  ldcaZ->Draw();
285 
286  //===================================
287 
288 
289  // extract pT resolution vs pT from hpt_compare
290 
291  TCanvas *c4 = new TCanvas("c4","c4",60,60,800,600);
292 
293  double dpT[200];
294 
295  dpT[0] = 1.0; // throw it away, stats are too poor
296  for(int i = 1;i<NPT;i++)
297  //for(int i = 0;i<1;i++)
298  {
299  // divide pT range into bins of width 0.5 GeV/c at 0.25, 0.75, 1.25, .....
300  // bins will be from 0-0.5, 0.5-1.0, 1.0-1.5, .....
301  double ptlo = (double) i * ptstep + 0.0;
302  double pthi = (double) i * ptstep + ptstep;
303 
304  int binlo = hpt_compare->GetXaxis()->FindBin(ptlo);
305  int binhi = hpt_compare->GetXaxis()->FindBin(pthi);
306 
307  TH1D *hpt = new TH1D("hpt","pT resolution ",200, 0, 2.0);
308  hpt_compare->ProjectionY("hpt",binlo,binhi);
309  hpt->GetXaxis()->SetTitle("#Delta p_{T}/p_{T}");
310  hpt->GetXaxis()->SetTitleOffset(1.0);
311  if(i>30) hpt->Rebin(4);
312  hpt->DrawCopy();
313 
314  std::cout << "ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << " integral " << hpt->Integral() << std::endl;
315 
316  TF1 *f = new TF1("f","gaus",0.8,1.2);
317  hpt->Fit(f,"R");
318 
319  pT[i] = (ptlo + pthi) / 2.0;
320 
321  dpT[i] = f->GetParameter(2);
322  cout << " pT " << pT[i] << " dpT " << dpT[i] << " integral " << hpt->Integral() << std::endl;
323  }
324 
325  //==========================================
326  // Plot pT resolution extracted above for embedded pions
327 
328  TCanvas *c5 = new TCanvas("c5","c5",100,100,800,800);
329  c5->SetLeftMargin(0.14);
330  TGraph *grdpt = new TGraph(NPT,pT,dpT);
331  grdpt->SetMarkerStyle(20);
332  grdpt->SetMarkerSize(1.1);
333  grdpt->SetName("pt_resolution");
334  grdpt->SetTitle("pT resolution");
335 
336  TH1D *hdummy2 = new TH1D("hdummy2","#Delta p_{T} vs p_{T}",100,0.0,hptmax);
337  hdummy2->SetMinimum(0);
338  //hdummy2->SetMaximum(0.05);
339  //hdummy2->SetMaximum(0.12);
340  hdummy2->SetMaximum(0.3);
341  hdummy2->GetXaxis()->SetTitle("p_{T}");
342  hdummy2->GetYaxis()->SetTitle("#Delta p_{T}/p_{T}");
343  hdummy2->GetYaxis()->SetTitleOffset(1.5);
344  hdummy2->Draw();
345  grdpt->Draw("p");
346 
347  TLegend *ldpt= new TLegend(0.25, 0.75, 0.95, 0.90, dcalab, "NDC");
348  ldpt->SetBorderSize(0);
349  ldpt->SetFillColor(0);
350  ldpt->SetFillStyle(0);
351  //char lstr1[500];
352  //sprintf(lstr1,"MVTX+INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
353  sprintf(lstr1,"INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
354  ldpt->AddEntry(grdpt, lstr1, "p");
355  ldpt->Draw();
356 
357  // Parameterize pT resolution
358 
359  //TF1 *fpt = new TF1("fpt","sqrt([0]*[0] + [1]*[1]*x*x)", 0, 35.0);
360  TF1 *fpt = new TF1("fpt","[0] + [1]*x", 2.0, hptmax);
361  fpt->SetParameter(0,0.005);
362  //fpt->FixParameter(0,0.005);
363  fpt->SetParameter(1,0.0015);
364  if(pt_resolution_fit)
365  grdpt->Fit(fpt,"R");
366 
367  char lab[1000];
368  //sprintf(lab,"#frac{#Deltap_{T}}{p_{T}} = #sqrt{%.4f^{2} + (%.6f #times p_{T})^{2}}", fpt->GetParameter(0), fpt->GetParameter(1))
369  sprintf(lab,"#frac{#Deltap_{T}}{p_{T}} = %.4f + %.6f #times p_{T}", fpt->GetParameter(0), fpt->GetParameter(1));
370  TLatex *mres = new TLatex(0.45,0.25,lab);
371  //mres->SetTextSize(0.1);
372  mres->SetNDC();
373  if(pt_resolution_fit)
374  mres->Draw();
375 
376  // For making a cut on the momentum difference between rgpT and rpT
377  double pT_sigmas = 6.0;
378  double const_term = fpt->GetParameter(0);
379  double linear_term = fpt->GetParameter(1);
380 
381  TH1D *hpt_truth;
382  fin->GetObject("hpt_truth",hpt_truth);
383  if(!hpt_truth)
384  {
385  cout << "Failed to get hpt_truth, quit!" << endl;
386  exit(1);
387  }
388 
389  //TCanvas *ctruth = new TCanvas("ctruth","ctruth", 5,5,800,600);
390 
391  TH1D *hpt_matched = new TH1D("hpt_matched","hpt_matched", 500, 0.0, hptmax);
392  double eff_pt[200];
393 
394  for(int i = 0;i<NPT;i++)
395  {
396  // divide pT range into bins of width 0.5 GeV/c at 0.25, 0.75, 1.25, .....
397  // bins will be from 0-0.5, 0.5-1.0, 1.0-1.5, .....
398  double ptlo = (double) i * ptstep + 0.0;
399  double pthi = (double) i * ptstep + ptstep;
400  double ptval = (ptlo+pthi)/2.0;
401 
402  int binlo = hpt_dca2d->GetXaxis()->FindBin(ptlo);
403  int binhi = hpt_dca2d->GetXaxis()->FindBin(pthi);
404 
405  TH1D *hpt1 = new TH1D("hpt1","pT resolution ",2000, 0.0, 2.0);
406  hpt_compare->ProjectionY("hpt1",binlo,binhi);
407 
408  double ptres = sqrt(pow(const_term,2) + pow(linear_term * ptval,2));
409  double ptreslo = 1.0 - ptres * pT_sigmas;
410  double ptreshi = 1.0 + ptres * pT_sigmas;
411 
412  int momlo = hpt1->GetXaxis()->FindBin(ptreslo);
413  int momhi = hpt1->GetXaxis()->FindBin(ptreshi);
414 
415  hpt_matched->Fill(ptval, hpt1->Integral(momlo, momhi));
416 
417  int tlo = hpt_truth->FindBin(ptlo);
418  int thi = hpt_truth->FindBin(pthi);
419  double truth_yield = hpt_truth->Integral(tlo, thi);;
420  eff_pt[i] = hpt1->Integral(momlo, momhi) / truth_yield;
421  //eff_pt[i] = hpt1->Integral() / truth_yield;
422  cout << " pT " << ptval << " ptres " << ptres << " ptreslo " << ptreslo << " ptreshi " << ptreshi << " momlo " << momlo << " momhi " << momhi << endl;
423  cout << " truth_yield " << truth_yield << " yield " << hpt1->Integral(momlo,momhi) << " eff_pt " << eff_pt[i] << endl;
424 
425  /*
426  if(i == 20)
427  {
428  ctruth->cd(0);
429  hpt1->DrawCopy();
430  TF1 *fres = new TF1("fres","gaus");
431  fres->SetParameter(0, 20.0);
432  fres->FixParameter(1, 1.0);
433  fres->FixParameter(2, ptres);
434  fres->SetLineColor(kRed);
435  hpt1->Fit(fres);
436  fres->DrawCopy("same");
437  ctruth->Update();
438  cout << hpt1->Integral() << endl;
439  cout << ptres << endl;
440  int k = 0;
441  int y=0;
442  cin >> k >> y;
443  }
444  */
445 
446  cout << " ptval " << ptval
447  << " ptreslo " << ptreslo
448  << " ptreshi " << ptreshi
449  << " total " << hpt1->Integral()
450  << " momlo " << momlo
451  << " momhi " << momhi
452  << " cut " << hpt1->Integral(momlo,momhi)
453  << " tlo " << tlo
454  << " thi " << thi
455  << " truth " << truth_yield
456  << " eff_pt " << eff_pt[i]
457  << endl;
458 
459  delete hpt1;
460 
461  }
462 
463  cout << " create canvas c7" << endl;
464  TCanvas *c7 = new TCanvas("c7","c7",60,60,800,800);
465 
466  TH1F *hd = new TH1F("hd","hd",100, 0.0, hptmax);
467  hd->SetMinimum(0.0);
468  hd->SetMaximum(1.0);
469  hd->GetXaxis()->SetTitle("p_{T} (GeV/c)");
470  hd->GetYaxis()->SetTitle("Single track efficiency");
471  hd->GetYaxis()->SetTitleOffset(1.0);
472  hd->Draw();
473 
474  TGraph *gr_eff = new TGraph(NPT,pT,eff_pt);
475  gr_eff->SetName("single_track_efficiency");
476  gr_eff->SetMarkerStyle(20);
477  gr_eff->SetMarkerSize(1);
478  gr_eff->SetMarkerColor(kRed);
479 
480  gr_eff->Draw("p");
481 
482  TLegend *leff = new TLegend(0.40, 0.35, 0.99, 0.50, dcalab, "NDC");
483  leff->SetBorderSize(0);
484  leff->SetFillColor(0);
485  leff->SetFillStyle(0);
486  char lstr[500];
487  //sprintf(lstr,"MVTX+INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
488  sprintf(lstr,"INTT+TPC",n_maps_layer,n_intt_layer,n_tpc_layer);
489  leff->AddEntry(gr_eff, lstr, "p");
490  //leff->AddEntry(lmax, "max possible efficiency", "l");
491  leff->Draw();
492 
493  /*
494  //=======================
495  // Get track purity for Hijing events by
496  // finding tracks within pt_sigmas of the
497  // truth pT
498 
499  TH2D *hpt_hijing_compare = 0;
500  fin->GetObject("hpt_hijing_compare",hpt_hijing_compare);
501  if(!hpt_hijing_compare)
502  {
503  cout << "Did not get hpt_hijing_compare - quit!" << endl;
504  exit(1);
505  }
506 
507  static const int NVARBINS = 36;
508  double xbins[NVARBINS+1] = {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,
509  1.1, 1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,
510  2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,
511  4.5,5.0,5.5,6.0,7.0,8.0,10.0};
512 
513  TH1D *hpt_hijing_allreco = new TH1D("hpt_hijing_allreco","hpt_hijing_allreco",NVARBINS,xbins);
514  TH1D *hpt_hijing_matched = new TH1D("hpt_hijing_matched","hpt_hijing_matched",NVARBINS,xbins);
515  hpt_hijing_matched->GetYaxis()->SetTitle("Purity");
516 
517  //double purity_pt[NVARBINS];
518 
519  for (int i=0;i<NVARBINS;i++)
520  {
521  double ptlo = xbins[i];
522  double pthi = xbins[i+1];
523  double ptval = (ptlo+pthi)/2.0;
524 
525  int binlo = hpt_hijing_compare->GetXaxis()->FindBin(ptlo);
526  int binhi = hpt_hijing_compare->GetXaxis()->FindBin(pthi);
527 
528  TH1D *hpt1 = new TH1D("hpt1","pT resolution ",2000, 0.0, 2.0);
529  hpt_hijing_compare->ProjectionY("hpt1",binlo,binhi);
530 
531  double ptres = sqrt(pow(const_term,2) + pow(linear_term * ptval,2));
532  double ptreslo = 1.0 - ptres * pT_sigmas;
533  double ptreshi = 1.0 + ptres * pT_sigmas;
534 
535  int momlo = hpt1->GetXaxis()->FindBin(ptreslo);
536  int momhi = hpt1->GetXaxis()->FindBin(ptreshi);
537 
538  // This to avoid TEfficiency's stupid insistence that the histograms cannot be filled with weights
539  for (int j=0;j<(int)hpt1->Integral(momlo,momhi);j++)
540  hpt_hijing_matched->Fill(ptval);
541  for (int j=0;j<(int)hpt1->Integral();j++)
542  hpt_hijing_allreco->Fill(ptval);
543 
544  //purity_pt[i] = hpt1->Integral(momlo, momhi) / hpt1->Integral();
545 
546  delete hpt1;
547  }
548 
549  cout << " create canvas cpur " << endl;
550  TCanvas *cpur = new TCanvas("cpur","cpur",60,60,1000,600);
551 
552  TEfficiency* pEff = 0;
553  //if(TEfficiency::CheckConsistency(*hpt_hijing_matched,*hpt_hijing_allreco))
554  {
555  pEff = new TEfficiency(*hpt_hijing_matched,*hpt_hijing_allreco);
556  char tname[500];
557  sprintf(tname,"Reconstruction efficiency (%.1f#sigma p_{T}) ; p_{T} (GeV/c) ; reco'd tracks within %.0f #sigma p_{T}",pT_sigmas,pT_sigmas);
558  //pEff->SetTitle("Reconstruction efficiency (3.0#sigma p_{T}) ; p_{T} (GeV/c) ; reco'd tracks within 3 #sigma p_{T}");
559  pEff->SetTitle(tname);
560  pEff->SetMarkerStyle(20);
561  pEff->SetMarkerColor(kBlack);
562  pEff->SetMarkerSize(1);
563  pEff->Draw();
564  gPad->Update();
565  pEff->GetPaintedGraph()->GetYaxis()->SetTitleOffset(1.0);
566  pEff->GetPaintedGraph()->GetYaxis()->SetRangeUser(0.0,1.1);
567  pEff->GetPaintedGraph()->GetXaxis()->SetTitleOffset(1.1);
568  }
569  */
570 
571  //=========================
572  // plot quality for Hijing tracks
573 
574  TCanvas *cpq = new TCanvas("cpq","cpq",40,40,1200,600);
575  TH1D * hquality;
576  fin->GetObject("hquality",hquality);
577  hquality->Draw();
578 
579  //======================
580  // Plot number of hits per track
581  TCanvas *c8 = new TCanvas("c8","c8",40,40,1200,600);
582  TH1D * hnhits;
583  fin->GetObject("hnhits",hnhits);
584  hnhits->GetXaxis()->SetTitle("hits per track");
585  hnhits->Draw();
586  cout << "hnhits integral = " << hnhits->Integral() << endl;
587 
588 
589  //================================
590  // extract three DCA 2d distributions inside pT bins
591 
592  TCanvas *cdcadist = new TCanvas("cdcadist","cdcadist",5,20,1200,800);
593  cdcadist->Divide(3,1);
594 
595  double dcaptbinlo[3] = {0.5, 1.0, 2.0};
596  double dcaptbinhi[3] = {1.0, 2.0, 50.0};
597  for(int i = 0;i<3;i++)
598  {
599  double ptlo = dcaptbinlo[i];
600  double pthi = dcaptbinhi[i];
601 
602  int binlo = hpt_dca2d->GetXaxis()->FindBin(ptlo);
603  int binhi = hpt_dca2d->GetXaxis()->FindBin(pthi);
604 
605  std::cout << "DCA distr: ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << std::endl;
606 
607  TH1D *h = new TH1D("h","dca2d resolution",2000, -0.1, 0.1);
608  hpt_dca2d->ProjectionY("h",binlo,binhi);
609  h->GetXaxis()->SetTitle("DCA (r#phi) (cm)");
610  h->GetXaxis()->SetTitleOffset(1.0);
611  h->SetMinimum(0.5);
612 
613  cdcadist->cd(i+1);
614  gPad->SetLogy(1);
615  h->Sumw2();
616  h->Rebin(4);
617  //h->SetMarkerStyle(20);
618  h->SetMarkerSize(0.8);
619  h->GetXaxis()->SetNdivisions(505);
620  h->DrawCopy("E");
621 
622  double mean = h->GetMean();
623  double sigma = h->GetRMS();
624  double yield = h->Integral();
625 
626  double low = -0.01, high=0.01;
627  if(n_maps_layer == 0)
628  {
629  low = 3.0 * low;
630  high = 3.0 * high;
631  }
632 
633  if(pthi < 6.0)
634  {
635  low = 3.0*low;
636  high = 3.0*high;
637  }
638 
639  TF1 *f = new TF1("f","gaus",low, high);
640  f->SetParameter(1, yield/100.0);
641  f->SetParameter(2, 0.0);
642  f->SetParameter(3,0.002);
643  f->SetLineColor(kRed);
644  h->Fit(f,"R");
645 
646  char fitr[500];
647  sprintf(fitr,"p_{T} = %.1f-%.1f GeV/c",ptlo,pthi);
648  TLatex *l1 = new TLatex(0.2,0.971,fitr);
649  l1->SetNDC();
650  l1->SetTextSize(0.07);
651  l1->Draw();
652 
653  sprintf(fitr,"#sigma = %.1f #mum", f->GetParameter(2)*10000);
654  TLatex *l1a = new TLatex(0.2,0.892,fitr);
655  l1a->SetNDC();
656  l1a->SetTextSize(0.07);
657  l1a->Draw();
658 
659  cout << " pT range " << ptlo << " " << pthi << " counts " << h->Integral() << " hist mean " << h->GetMean() << " hist RMS " << h->GetRMS() << endl;
660  }
661 
662  //================================
663  // extract three DCA Z distributions inside pT bins
664 
665  TCanvas *cdcazdist = new TCanvas("cdcazdist","cdcazdist",5,20,1200,800);
666  cdcazdist->Divide(3,1);
667 
668  for(int i = 0;i<3;i++)
669  {
670  double ptlo = dcaptbinlo[i];
671  double pthi = dcaptbinhi[i];
672 
673  int binlo = hpt_dcaZ->GetXaxis()->FindBin(ptlo);
674  int binhi = hpt_dcaZ->GetXaxis()->FindBin(pthi);
675 
676  std::cout << "DCA Z distr: ptlo " << ptlo << " binlo " << binlo << " pthi " << pthi << " binhi " << binhi << std::endl;
677 
678  TH1D *h = new TH1D("h","dcaZ resolution",2000, -0.1, 0.1);
679  hpt_dcaZ->ProjectionY("h",binlo,binhi);
680  h->GetXaxis()->SetTitle("DCA (Z) (cm)");
681  h->GetXaxis()->SetTitleOffset(1.0);
682 
683  cdcazdist->cd(i+1);
684  gPad->SetLogy(1);
685  h->Sumw2();
686  h->Rebin(4);
687  h->SetMarkerSize(0.8);
688  h->GetXaxis()->SetNdivisions(505);
689  h->SetMinimum(0.5);
690  h->DrawCopy("E");
691 
692  double mean = h->GetMean();
693  double sigma = h->GetRMS();
694  double yield = h->Integral();
695 
696  double low = -0.01, high=0.01;
697  if(n_maps_layer == 0)
698  {
699  low = 3.0 * low;
700  high = 3.0 * high;
701  }
702 
703  if(pthi < 6.0)
704  {
705  low = 3.0*low;
706  high = 3.0*high;
707  }
708 
709  TF1 *f = new TF1("f","gaus",low, high);
710  f->SetParameter(1, yield/100.0);
711  f->SetParameter(2, 0.0);
712  f->SetParameter(3,0.002);
713  f->SetLineColor(kRed);
714  h->Fit(f,"R");
715 
716  char fitr[500];
717  sprintf(fitr,"p_{T} = %.1f-%.1f GeV/c",ptlo,pthi);
718  TLatex *l1 = new TLatex(0.2,0.971,fitr);
719  l1->SetNDC();
720  l1->SetTextSize(0.07);
721  l1->Draw();
722 
723  sprintf(fitr,"#sigma = %.1f #mum", f->GetParameter(2)*10000);
724  TLatex *l1a = new TLatex(0.2,0.892,fitr);
725  l1a->SetNDC();
726  l1a->SetTextSize(0.07);
727  l1a->Draw();
728 
729  cout << " pT range " << ptlo << " " << pthi << " counts " << h->Integral() << " hist mean " << h->GetMean() << " hist RMS " << h->GetRMS() << endl;
730  }
731 
732  /*
733 
734  //=========================
735 
736  TCanvas *cdca = new TCanvas("cdca","cdca",5,20,1200,800);
737  //cdca->SetTopMargin(0.2);
738  cdca->Divide(3,1);
739 
740  TH1D *hdca2d[3];
741  fin->GetObject("hdca2d0",hdca2d[0]);
742  fin->GetObject("hdca2d1",hdca2d[1]);
743  fin->GetObject("hdca2d2",hdca2d[2]);
744 
745  cdca->cd(1);
746  gPad->SetLogy(1);
747  //gPad->SetRightMargin(0.1);
748  hdca2d[0]->GetXaxis()->SetNdivisions(505);
749 
750  hdca2d[0]->Draw();
751 
752  cdca->cd(2);
753  gPad->SetLogy(1);
754  hdca2d[1]->GetXaxis()->SetNdivisions(505);
755  hdca2d[1]->Draw();
756 
757  cdca->cd(3);
758  gPad->SetLogy(1);
759  hdca2d[2]->GetXaxis()->SetNdivisions(505);
760  hdca2d[2]->Draw();
761 
762  TF1 *fdca = new TF1("fdca","gaus",-0.1,0.1);
763  fdca->SetLineColor(kRed);
764  cdca->cd(1);
765  hdca2d[0]->Fit(fdca);
766 
767  TLatex *l1a = new TLatex(0.2,0.971,"p_{T} = 0.5-1.0 GeV/c");
768  l1a->SetNDC();
769  l1a->SetTextSize(0.07);
770  l1a->Draw();
771  char fitr[500];
772  sprintf(fitr,"#sigma = %.1f #mum",fdca->GetParameter(2)*10000);
773  TLatex *l1 = new TLatex(0.2,0.892,fitr);
774  l1->SetNDC();
775  l1->SetTextSize(0.07);
776  l1->Draw();
777 
778  cdca->cd(2);
779  hdca2d[1]->Fit(fdca);
780 
781  TLatex *l2a = new TLatex(0.2,0.971,"p_{T} = 1.0-2.0 GeV/c");
782  l2a->SetNDC();
783  l2a->SetTextSize(0.07);
784  l2a->Draw();
785 
786  sprintf(fitr,"#sigma = %.1f #mum",fdca->GetParameter(2)*10000);
787  TLatex *l2 = new TLatex(0.2,0.892,fitr);
788  l2->SetNDC();
789  l2->SetTextSize(0.07);
790  l2->Draw();
791 
792 
793  cdca->cd(3);
794  hdca2d[2]->Fit(fdca);
795  sprintf(fitr,"#splitline{p_{T} > 2.0 GeV/c}{#sigma = %.1f #mum}",fdca->GetParameter(2)*10000);
796 
797  TLatex *l3a = new TLatex(0.2,0.971,"p_{T} = 1.0-2.0 GeV/c");
798  l3a->SetNDC();
799  l3a->SetTextSize(0.07);
800  l3a->Draw();
801 
802  sprintf(fitr,"#sigma = %.1f #mum",fdca->GetParameter(2)*10000);
803  TLatex *l3 = new TLatex(0.2,0.892,fitr);
804  l3->SetNDC();
805  l3->SetTextSize(0.07);
806  l3->Draw();
807  */
808 
809  TCanvas *ceta = new TCanvas("ceta","ceta",10,10,600,600);
810 
811  TH1D *hgeta = 0;
812  fin->GetObject("hgeta",hgeta);
813  if(!hgeta)
814  {
815  cout << "Did not get hgeta" << endl;
816  exit(1);
817  }
818  TH1D *hreta = 0;
819  fin->GetObject("hreta",hreta);
820 
821  hreta->GetXaxis()->SetTitle("#eta");
822  hreta->GetYaxis()->SetNdivisions(505);
823  hreta->SetLineColor(kRed);
824  hreta->SetMinimum(0.0);
825  hreta->Draw();
826 
827  hgeta->Draw("same");
828 
829 
830  TLegend *leta = new TLegend(0.3, 0.35, 1.0, 0.60, dcalab, "NDC");
831  leta->SetBorderSize(0);
832  leta->SetFillColor(0);
833  leta->SetFillStyle(0);
834  leta->AddEntry(hgeta,"Truth","l");
835  leta->AddEntry(hreta,"Reconstructed","l");
836  leta->Draw();
837  ceta->Update();
838 
839  TCanvas *cfake = new TCanvas("cfake","cfake",4,4,800,600);
840  cfake->SetLeftMargin(0.15);
841  TH2D *hpt_nfake = 0;
842  fin->GetObject("hpt_nfake",hpt_nfake);
843  if(!hpt_nfake)
844  {
845  cout << "Did not get hpt_nfake" << endl;
846  exit(1);
847  }
848 
849  int n_noise = 4;
850  double ynoise[10][200];
851  double ynoise_ref[200];
852  //double dptfake = 1.0;
853  double dptfake = 0.2;
854  double fake_ptmax = 20.0;
855  int nptfake = (int) (fake_ptmax/dptfake);
856  double ptfake[200];
857  for(int ipt=0;ipt<nptfake;ipt++)
858  {
859  double ptlo = dptfake * ipt;
860  double pthi = dptfake*ipt + dptfake;
861  ptfake[ipt] = (ptlo+pthi) / 2.0;
862  int binlo = hpt_nfake->GetXaxis()->FindBin(ptlo);
863  int binhi = hpt_nfake->GetXaxis()->FindBin(pthi);
864 
865  TH1D *hnoise = new TH1D("hnoise","",73,-5,68);
866  hpt_nfake->ProjectionY("hnoise",binlo,binhi);
867  double ynoise_ref = hnoise->Integral();
868  for(int j=0;j<n_noise;j++)
869  {
870  int bin = hnoise->FindBin(j);
871  ynoise[j][ipt]= hnoise->GetBinContent(bin);
872  ynoise[j][ipt] /= ynoise_ref;
873  }
874  hnoise->Delete();
875  }
876 
877 
878  TH1D *hfdummy = new TH1D("hfdummy","",100,0.0,fake_ptmax);
879  hfdummy->SetMaximum(1.05);
880  hfdummy->SetMinimum(0.001);
881  hfdummy->GetXaxis()->SetTitle("p_{T} (GeV/c)");
882  hfdummy->GetYaxis()->SetTitle("Fraction of tracks");
883  gPad->SetLogy(1);
884  hfdummy->Draw();
885 
886  TLegend *lfake = new TLegend(0.35,0.55,0.87,0.89,"","NDC");
887  lfake->SetBorderSize(1);
888  TGraph *hnoise[10];
889  int fake_col[5] = {kBlack,kRed,kBlue,kMagenta,kViolet};
890  for(int j=0;j<n_noise;j++)
891  {
892  hnoise[j] = new TGraph(nptfake,ptfake,ynoise[j]);
893  hnoise[j]->SetMarkerStyle(20);
894  hnoise[j]->SetMarkerSize(1.0);
895  hnoise[j]->SetMarkerColor(fake_col[j]);
896  if(j == 0)
897  hnoise[j]->Draw("p");
898  else
899  hnoise[j]->Draw("p same");
900  char lab[500];
901  sprintf(lab,"%i mismatched hits",j);
902  lfake->AddEntry(hnoise[j],lab,"p");
903  }
904  lfake->Draw();
905 
906  TH2D *hcorr_nfake_nmaps = 0;
907  fin->GetObject("hcorr_nfake_nmaps",hcorr_nfake_nmaps);
908  if(hcorr_nfake_nmaps)
909  {
910  TCanvas *cmm = new TCanvas("cmm","cmm",4,4,800,600);
911  cmm->SetLeftMargin(0.15);
912 
913  gPad->SetLogy(1);
914 
915  TLegend *lmm = new TLegend(0.35,0.62,0.9,0.90,"","NDC");
916  lmm->SetBorderSize(1);
917  for(int imaps =0;imaps<n_maps_layer+1;imaps++)
918  {
919  double nmapslo = -0.5 + 1.0 * imaps;
920  double nmapshi = nmapslo + 1.0;
921 
922  int binlo = hcorr_nfake_nmaps->GetYaxis()->FindBin(nmapslo);
923  int binhi = hcorr_nfake_nmaps->GetYaxis()->FindBin(nmapshi);
924 
925  TH1D *h = new TH1D("h","h",40,-1,4);
926  hcorr_nfake_nmaps->ProjectionX("h",binlo,binhi);
927 
928  h->SetMarkerStyle(20);
929  h->SetMarkerSize(1.5);
930  h->SetMarkerColor(fake_col[imaps]);
931  h->SetMinimum(10);
932  h->SetMaximum(5.0e6);
933  char lab[500];
934  sprintf(lab,"%i missed maps layers",imaps);
935  lmm->AddEntry(h,lab,"p");
936 
937  if(imaps == 0)
938  {
939  h->GetXaxis()->SetTitle("Mismatched hits");
940  h->GetYaxis()->SetTitle("Tracks");
941  h->DrawCopy("p");
942 
943  }
944  else
945  h->DrawCopy("same p");
946  }
947 
948  lmm->Draw();
949  }
950 
951  TCanvas *cvtx = new TCanvas("cvtx","cvtx",4,4,800,600);
952  TH1D *hzevt = 0;
953  fin->GetObject("hzevt",hzevt);
954  if(!hzevt)
955  {
956  cout << "Did not get hzevt" << endl;
957  exit(1);
958  }
959  hzevt->Draw();
960 
961  TCanvas *cvtx_res = new TCanvas("cvtx_res","cvtx_res",4,4,1200,800);
962  cvtx_res->Divide(2,2);
963 
964  TH1D *hzvtx_res = 0;
965  fin->GetObject("hzvtx_res",hzvtx_res);
966  if(!hzvtx_res)
967  {
968  cout << "Did not get hzvtx_res" << endl;
969  }
970 
971  if(hzvtx_res)
972  {
973  cvtx_res->cd(1);
974 
975  hzvtx_res->SetNdivisions(505);
976  hzvtx_res->Draw();
977 
978  TF1 *fvtxres = new TF1("fvtxres","gaus");
979  //hzvtx_res->Fit(fvtxres);
980  }
981 
982  TH1D *hxvtx_res = 0;
983  fin->GetObject("hxvtx_res",hxvtx_res);
984  if(!hxvtx_res)
985  {
986  cout << "Did not get hxvtx_res" << endl;
987  }
988 
989  if(hxvtx_res)
990  {
991  cvtx_res->cd(2);
992 
993  hxvtx_res->SetNdivisions(505);
994  hxvtx_res->Draw();
995 
996  TF1 *fvtxres = new TF1("fvtxres","gaus");
997  //hxvtx_res->Fit(fvtxres);
998  }
999 
1000  TH1D *hyvtx_res = 0;
1001  fin->GetObject("hyvtx_res",hyvtx_res);
1002  if(!hyvtx_res)
1003  {
1004  cout << "Did not get hyvtx_res" << endl;
1005  }
1006 
1007  if(hyvtx_res)
1008  {
1009  cvtx_res->cd(3);
1010 
1011  hyvtx_res->SetNdivisions(505);
1012  hyvtx_res->Draw();
1013 
1014  TF1 *fvtxres = new TF1("fvtxres","gaus");
1015  //hyvtx_res->Fit(fvtxres);
1016  }
1017 
1018  TH1D *hdcavtx_res = 0;
1019  fin->GetObject("hdcavtx_res",hdcavtx_res);
1020  if(!hdcavtx_res)
1021  {
1022  cout << "Did not get hdcavtx_res" << endl;
1023  }
1024 
1025  if(hdcavtx_res)
1026  {
1027  cvtx_res->cd(4);
1028 
1029  hdcavtx_res->SetNdivisions(505);
1030  hdcavtx_res->Draw();
1031 
1032  TF1 *fvtxres = new TF1("fvtxres","gaus");
1033  //hdcavtx_res->Fit(fvtxres);
1034  }
1035 
1036 
1037  TCanvas *g4ntr = new TCanvas("g4ntr","g4ntr",5,10,1200,800);
1038  g4ntr->Divide(3,1);
1039  g4ntr->cd(1);
1040  TH2D *hg4ntrack = 0;
1041  fin->GetObject("hg4ntrack",hg4ntrack);
1042  if(hg4ntrack)
1043  {
1044  hg4ntrack->GetXaxis()->SetTitle(" Tracks/event from truth");
1045  hg4ntrack->GetYaxis()->SetTitle(" Tracks/event reconstructed");
1046  hg4ntrack->Draw();
1047  }
1048  else
1049  cout << "Did not get hgn4track" << endl;
1050 
1051  g4ntr->cd(2);
1052  TH1D *hg4ntr = new TH1D("hg4ntr","hg4ntr",200, 0.0, 2000.0);
1053  hg4ntrack->ProjectionX("hg4ntr");
1054  hg4ntr->GetXaxis()->SetTitle(" Tracks/event from truth");
1055  hg4ntr->GetYaxis()->SetTitle(" Yield");
1056  hg4ntr->Draw();
1057 
1058  g4ntr->cd(3);
1059  double cent[20];
1060  double nhits[20];
1061  nhits[0] = 0;
1062  double ytot = hg4ntr->Integral();
1063  for(int i=1;i<20;i++)
1064  {
1065  nhits[i] = i*120.0;
1066  int binhi = hg4ntr->FindBin(nhits[i]);
1067  double y = hg4ntr->Integral(1,binhi);
1068  cent[i] = y / ytot;
1069  cout << "nhits = " << nhits[i] << " binhi = " << binhi << " y = " << y << " ytot = " << ytot << " cent = " << cent[i] << endl;
1070  }
1071  TGraph *gy = new TGraph(20,nhits,cent);
1072  gy->GetXaxis()->SetTitle(" Tracks/event from truth");
1073  gy->GetYaxis()->SetTitle(" integral fraction of events");
1074  gy->Draw();
1075 
1076  TF1 *fgy = new TF1("fgy","[0]*x+[1]*x*x+[2]*x*x*x",0,2300);
1077  gy->Fit(fgy,"R");
1078 
1079  // Output some graphs for comparison plots
1080 
1081  TFile *fout = new TFile("root_files/look_purity_out.root","recreate");
1082  gr_eff->Write();
1083  grdca2d->Write();
1084  grdcaZ->Write();
1085  grdpt->Write();
1086 
1087 }