Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DrawTrackingEval.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DrawTrackingEval.C
1 #include "SaveCanvas.C"
2 #include "sPhenixStyle.C"
3 
4 #include <TChain.h>
5 #include <TCut.h>
6 #include <TEfficiency.h>
7 #include <TF1.h>
8 #include <TGraphAsymmErrors.h>
9 #include <TGraphErrors.h>
10 #include <TH2.h>
11 #include <TH3.h>
12 
13 #include <TFile.h>
14 
15 #include <TColor.h>
16 #include <TLatex.h>
17 #include <TLegend.h>
18 #include <TLine.h>
19 
20 #include <TMath.h>
21 #include <TString.h>
22 #include <TTree.h>
23 #include <TVirtualFitter.h>
24 
25 #include <cmath>
26 #include <iostream>
27 
28 using namespace std;
29 
30 // ROOT6 disabled assert. Well....
31 #ifdef assert
32 #undef assert
33 #endif
34 #define assert(exp) \
35  { \
36  if (!(exp)) \
37  { \
38  cout << "Assert (" << #exp << ") failed at " << __FILE__ << " line " << __LINE__ << endl; \
39  exit(1); \
40  } \
41  }
42 
43 TTree *ntp_gtrack(nullptr);
44 TFile *_file0 = NULL;
45 TString description;
46 TCut primPartCut("gpt>.1 && abs(geta)<1 && abs(gvt)<1e-10"); // primary particle in acceptance
47 TCut recoCut("nlmaps>=2 && pt>0");
48 const double pTMin(.3);
49 
50 //#define ASSERT(exp) {{exit(1);} }
51 
53 void useLogBins(TAxis *axis)
54 {
55  assert(axis);
56  assert(axis->GetXmin() > 0);
57  assert(axis->GetXmax() > 0);
58 
59  const int bins = axis->GetNbins();
60 
61  Axis_t from = log10(axis->GetXmin());
62  Axis_t to = log10(axis->GetXmax());
63  Axis_t width = (to - from) / bins;
64  vector<Axis_t> new_bins(bins + 1);
65 
66  for (int i = 0; i <= bins; i++)
67  {
68  new_bins[i] = TMath::Power(10, from + i * width);
69  }
70  axis->Set(bins, new_bins.data());
71 }
72 
73 TGraphErrors *
74 FitProfile(const TH2 *h2)
75 {
76  TProfile *p2 = h2->ProfileX();
77 
78  int n = 0;
79  double x[1000];
80  double ex[1000];
81  double y[1000];
82  double ey[1000];
83 
84  for (int i = 1; i <= h2->GetNbinsX(); i++)
85  {
86  TH1D *h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
87 
88  if (h1->GetSum() < 30)
89  {
90  cout << "FitProfile - ignore bin " << i << endl;
91  continue;
92  }
93  else
94  {
95  // cout << "FitProfile - fit bin " << i << endl;
96  }
97 
98  TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
99  p2->GetBinError(i) * 4);
100 
101  // TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
102  // -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
103 
104  fgaus.SetParameter(1, p2->GetBinContent(i));
105  fgaus.SetParameter(2, p2->GetBinError(i));
106 
107  h1->Fit(&fgaus, "MQ");
108 
109  TF1 f2(Form("dgaus"), "gaus + [3]",
110  -fgaus.GetParameter(2) * 1.5, fgaus.GetParameter(2) * 1.5);
111 
112  f2.SetParameters(fgaus.GetParameter(0), fgaus.GetParameter(1),
113  fgaus.GetParameter(2), 0);
114 
115  h1->Fit(&f2, "MQR");
116 
117  // new TCanvas;
118  // h1->Draw();
119  // fgaus.Draw("same");
120  // break;
121 
122  x[n] = p2->GetBinCenter(i);
123  ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
124  y[n] = fgaus.GetParameter(2);
125  ey[n] = fgaus.GetParError(2);
126 
127  // p2->SetBinContent(i, fgaus.GetParameter(1));
128  // p2->SetBinError(i, fgaus.GetParameter(2));
129 
130  n++;
131  delete h1;
132  }
133 
134  return new TGraphErrors(n, x, y, ex, ey);
135 }
136 
137 void RecoCheck(const TString &name = "pi", const TCut &cut = "abs(gflavor)==211", const int nBins = 100)
138 {
139  assert(_file0);
140 
141  TCanvas *c1 = new TCanvas("RecoCheck_" + name, "RecoCheck_" + name, 1800, 960);
142  c1->Divide(3, 2);
143  int idx = 1;
144  TPad *p;
145 
146  p = (TPad *) c1->cd(idx++);
147  c1->Update();
148  p->SetLogx();
149 
150  p->DrawFrame(pTMin, 0, 20, 1.1, "Efficiency;Truth p_{T} [GeV/c];Tracking efficiency");
151  TLine *l = new TLine(pTMin, 1, 20, 1);
152  l->Draw();
153 
154  TH1 *hTruthpT = new TH1F("hTruthpT", "hTruthpT;Truth p_{T} [GeV/c]", nBins, pTMin, 20);
155  useLogBins(hTruthpT->GetXaxis());
156  TH1 *hRecopT = new TH1F("hRecopT", "hRecopT;Truth p_{T} [GeV/c]", nBins, pTMin, 20);
157  useLogBins(hRecopT->GetXaxis());
158  TH1 *hRecopTMVTX = new TH1F("hRecopTMVTX", "hRecopT;Truth p_{T} [GeV/c]", nBins, pTMin, 20);
159  useLogBins(hRecopTMVTX->GetXaxis());
160 
161  ntp_gtrack->Draw("gpt>>hTruthpT", cut && primPartCut, "goff");
162  ntp_gtrack->Draw("gpt>>hRecopT", cut && primPartCut && "pt>0", "goff");
163  ntp_gtrack->Draw("gpt>>hRecopTMVTX", cut && primPartCut && recoCut, "goff");
164 
165  TEfficiency *hEffpTMVTX = new TEfficiency(*hRecopTMVTX, *hTruthpT);
166  hEffpTMVTX->SetName("hEffpTMVTX");
167  hEffpTMVTX->SetTitle("Efficiency;Truth p_{T} [GeV/c];Tracking efficiency");
168  hEffpTMVTX->Draw("p same");
169  hEffpTMVTX->SetLineColor(kBlue + 2);
170  hEffpTMVTX->SetMarkerColor(kBlue + 2);
171 
172  // TH1 * hDCAxy = new TH1F("hDCAxy","")
173  TEfficiency *hEffpT = new TEfficiency(*hRecopT, *hTruthpT);
174  hEffpT->SetName("hEffpT");
175  hEffpT->SetTitle("Efficiency;Truth p_{T} [GeV/c];Tracking efficiency");
176  hEffpT->Draw("P same");
177  hEffpT->SetLineColor(kRed + 2);
178  hEffpT->SetMarkerColor(kRed + 2);
179 
180  // TH2F * hDCA;
181  // ntp_gtrack -> Draw();
182 
183  p = (TPad *) c1->cd(idx++);
184  c1->Update();
185  p->SetLogx();
186  p->SetLogz();
187 
188  TH2 *hpTResolution = new TH2F("hpTResolution", "hpTResolution;Truth p_{T} [GeV/c];#Deltap_{T}/p_{T}", nBins, pTMin, 20, 1000, -.2, .2);
189  useLogBins(hpTResolution->GetXaxis());
190  hpTResolution->GetYaxis()->SetRangeUser(-.1, .1);
191 
192  ntp_gtrack->Draw("pt/gpt - 1 : gpt>>hpTResolution", cut && primPartCut && recoCut, "goff");
193 
194  TGraphErrors *gepTResolution = FitProfile(hpTResolution);
195  gepTResolution->SetName("gepTResolution");
196 
197  hpTResolution->Draw("colz");
198  gepTResolution->Draw("p");
199 
200  p = (TPad *) c1->cd(idx++);
201  c1->Update();
202  p->SetLogx();
203  p->SetLogz();
204 
205  TH2 *hcotThetaResolution = new TH2F("hcotThetaResolution", "hcotThetaResolution;Truth p_{T} [GeV/c];#Delta[p_{z}/p_{T}]", nBins, pTMin, 20, 1000, -.1, .1);
206  useLogBins(hcotThetaResolution->GetXaxis());
207  hcotThetaResolution->GetYaxis()->SetRangeUser(-.01, .01);
208 
209  ntp_gtrack->Draw("pz/pt - gpz/gpt: gpt>>hcotThetaResolution", cut && primPartCut && recoCut, "goff");
210 
211  TGraphErrors *gecotThetaResolution = FitProfile(hcotThetaResolution);
212  gecotThetaResolution->SetName("gecotThetaResolution");
213 
214  hcotThetaResolution->Draw("colz");
215  gecotThetaResolution->Draw("p");
216 
217  p = (TPad *) c1->cd(idx++);
218  c1->Update();
219  p->SetLogx();
220  p->SetLogz();
221 
222  TH2 *hsinPhiResolution = new TH2F("hsinPhiResolution", "hsinPhiResolution;Truth p_{T} [GeV/c];sin(#Delta#phi)", nBins, pTMin, 20, 1000, -.1, .1);
223  useLogBins(hsinPhiResolution->GetXaxis());
224  hsinPhiResolution->GetYaxis()->SetRangeUser(-.01, .01);
225 
226  ntp_gtrack->Draw("sin(atan2(py,px) - atan2(gpy,gpx)) : gpt>>hsinPhiResolution", cut && primPartCut && recoCut, "goff");
227 
228  TGraphErrors *gesinPhiResolution = FitProfile(hsinPhiResolution);
229  gesinPhiResolution->SetName("gesinPhiResolution");
230 
231  hsinPhiResolution->Draw("colz");
232  gesinPhiResolution->Draw("p");
233 
234  p = (TPad *) c1->cd(idx++);
235  c1->Update();
236  p->SetLogx();
237  p->SetLogz();
238 
239  TH2 *hDCA3DXYResolution = new TH2F("hDCA3DXYResolution", "hDCA3DXYResolution;Truth p_{T} [GeV/c];DCA_{xy} [cm]", nBins, pTMin, 20, 1000, -.2, .2);
240  useLogBins(hDCA3DXYResolution->GetXaxis());
241  hDCA3DXYResolution->GetYaxis()->SetRangeUser(-.02, .02);
242 
243  ntp_gtrack->Draw("dca3dxy : gpt>>hDCA3DXYResolution", cut && primPartCut && recoCut, "goff");
244 
245  TGraphErrors *geDCA3DXYResolution = FitProfile(hDCA3DXYResolution);
246  geDCA3DXYResolution->SetName("geDCA3DXYResolution");
247 
248  hDCA3DXYResolution->Draw("colz");
249  geDCA3DXYResolution->Draw("p");
250 
251  p = (TPad *) c1->cd(idx++);
252  c1->Update();
253  p->SetLogx();
254  p->SetLogz();
255 
256  TH2 *hDCA3DZResolution = new TH2F("hDCA3DZResolution", "hDCA3DZResolution;Truth p_{T} [GeV/c];DCA_{z} [cm]", nBins, pTMin, 20, 1000, -.2, .2);
257  useLogBins(hDCA3DZResolution->GetXaxis());
258  hDCA3DZResolution->GetYaxis()->SetRangeUser(-.02, .02);
259 
260  ntp_gtrack->Draw("dca3dz : gpt>>hDCA3DZResolution", cut && primPartCut && recoCut, "goff");
261 
262  TGraphErrors *geDCA3DZResolution = FitProfile(hDCA3DZResolution);
263  geDCA3DZResolution->SetName("geDCA3DZResolution");
264 
265  hDCA3DZResolution->Draw("colz");
266  geDCA3DZResolution->Draw("p");
267 
268  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), false);
269 }
270 
271 void DCACheck(const TString &name = "pi", const TCut &cut = "abs(gflavor)==211", const int nBins = 14)
272 {
273  assert(_file0);
274 
275  const int nDCAbins = 2 * (100 / 1 * 2 + (1000 - 100) / 10 + (2000 - 1000) / 100);
276 
277  TH3F *hDCA3Dpz = new TH3F("hDCA3Dpz", "hDCA3Dpz;DCA_{xy} [cm];DCA_{z} [cm];Truth p_{T} [GeV/c]", //
278  nDCAbins, -.2, .2, //
279  nDCAbins, -.2, .2, //
280  nBins, pTMin, 20);
281 
282  vector<Axis_t> new_bins(nDCAbins + 1);
283  int binedge = -20000;
284  for (int i = 0; i <= nDCAbins; i++)
285  {
286  new_bins[i] = binedge * 1e-5; // convert .1um to cm
287 
288  double dedge = 1000;
289  if (abs(binedge + .5) <= 10000)
290  dedge = 100;
291  if (abs(binedge + .5) <= 1000)
292  dedge = 5;
293 
294  // cout << "new_bins[" << i << "] = " << binedge << "* 1e-4 = " << new_bins[i] << ", dedge = " << dedge << endl;
295 
296  binedge += dedge;
297  }
298  assert(binedge == 21000);
299  hDCA3Dpz->GetXaxis()->Set(nDCAbins, new_bins.data());
300  hDCA3Dpz->GetYaxis()->Set(nDCAbins, new_bins.data());
301  // cout << "hDCA3Dpz->GetXaxis()->GetBinCenter(nDCAbins) = " << hDCA3Dpz->GetXaxis()->GetBinCenter(nDCAbins) << endl;
302 
303  // ntp_gtrack->Draw("dca3dxy : dca3dz : gpt >> hDCA3Dpz", cut && primPartCut && recoCut, "goff");
304  ntp_gtrack->Draw("gpt : dca3dz : dca3dxy >> hDCA3Dpz", cut && primPartCut && recoCut, "goff");
305 
306  // plot
307  TCanvas *c1 = new TCanvas("DCACheck_" + name, "DCACheck_" + name, 1800, 960);
308 
309  c1->Divide(5, 3);
310  int idx = 1;
311  TPad *p;
312 
313  p = (TPad *) c1->cd(idx++);
314  c1->Update();
315  p->SetLogx();
316  p->SetLogz();
317 
318  TH2 *hDCA3DXYCheck = new TH2F("hDCA3DXYCheck", "hDCA3DXYCheck;Truth p_{T} [GeV/c];DCA_{xy}", nBins, pTMin, 20, 1000, -.2, .2);
319  useLogBins(hDCA3DXYCheck->GetXaxis());
320  hDCA3DXYCheck->GetYaxis()->SetRangeUser(-.02, .02);
321 
322  ntp_gtrack->Draw("dca3dxy : gpt>>hDCA3DXYCheck", cut && primPartCut && recoCut, "goff");
323 
324  TGraphErrors *geDCA3DXYCheck = FitProfile(hDCA3DXYCheck);
325  geDCA3DXYCheck->SetName("geDCA3DXYCheck");
326 
327  hDCA3DXYCheck->Draw("colz");
328  geDCA3DXYCheck->Draw("p");
329 
330  for (int bin = 1; bin <= nBins; bin++)
331  {
332  p = (TPad *) c1->cd(idx++);
333  c1->Update();
334  // p->SetLogx();
335  p->SetLogz();
336  p->SetRightMargin(.15);
337 
338  hDCA3Dpz->GetZaxis()->SetRange(bin, bin);
339  TH1 *h2D = hDCA3Dpz->Project3D("yx");
340  h2D->SetName(Form("hDCA3Dpz%d", bin));
341  h2D->Draw("colz");
342  }
343 
344  // ping the 3D histogram for saving
345  c1->GetListOfPrimitives()->Add(hDCA3Dpz);
346 
347  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), false);
348 }
349 
350 void DrawTrackingEval(
351  const TString infile = "/gpfs/mnt/gpfs02/sphenix/user/jinhuang/HF-jet/HF-production-piKp-pp200-dataset1_333ALL.cfg_ALL_g4svtx_eval.root",
352  const TString &name = "pi", const TCut &cut = "abs(gflavor)==211" //
353 )
354 {
355  SetsPhenixStyle();
356  TVirtualFitter::SetDefaultFitter("Minuit2");
357 
358 // description = disc;
359  gSystem->Load("libg4eval.so");
360 
361  if (!_file0)
362  {
363  TString chian_str = infile;
364  chian_str.ReplaceAll("ALL", "*");
365 
366  TChain *t = new TChain("ntp_gtrack");
367  const int n = t->Add(chian_str);
368 
369  cout << "Loaded " << n << " root files with " << chian_str << endl;
370  assert(n > 0);
371  ntp_gtrack = t;
372 
373  _file0 = new TFile;
374  _file0->SetName(infile);
375  }
376 
377 // gDirectory->mkdir("/pi");
378 // gDirectory->Cd("/pi");
379  RecoCheck(name, cut);
380  DCACheck(name, cut);
381 
382  // gDirectory->mkdir("/kaon");
383  // gDirectory->Cd("/kaon");
384  // RecoCheck("kaon", "abs(gflavor)==321");
385  // DCACheck("kaon", "abs(gflavor)==321");
386  //
387  // gDirectory->mkdir("/proton");
388  // gDirectory->Cd("/proton");
389  // RecoCheck("proton", "abs(gflavor)==2212");
390  // DCACheck("proton", "abs(gflavor)==2212");
391 }