Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DrawHadronShowers.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DrawHadronShowers.C
1 #include <TFile.h>
2 #include <TLine.h>
3 #include <TString.h>
4 #include <TTree.h>
5 #include <cassert>
6 #include <cmath>
7 #include "SaveCanvas.C"
8 #include "sPhenixStyle.C"
9 using namespace std;
10 
11 class lin_res
12 {
13  public:
14  TString name;
15  TGraphErrors *linearity;
16  TGraphErrors *resolution;
17  TF1 *f_res;
18 };
19 
20 
21 
22 //TString base_dataset = "/phenix/u/jinhuang/links/sPHENIX_work/prod_analysis/hadron_shower_res_nightly/";
23 //TString base_dataset = "/phenix/u/jinhuang/links/sPHENIX_work/prod_analysis/hadron_shower_res_hcaldigi/";
24 //TString base_dataset = "/phenix/u/jinhuang/links/sPHENIX_work/prod_analysis/hadron_shower_res_hcaldigi_sampling/";
25 //TString base_dataset ="/phenix/u/jinhuang/links/ePHENIX_work/sPHENIX_work/production_analysis_updates/spacal1d/fieldmap";
26 TString base_dataset = "/sphenix/u/weihuma/analysis/Calorimeter/macros/";
27 
28 void DrawHadronShowers(void)
29 {
30  SetsPhenixStyle();
31  TVirtualFitter::SetDefaultFitter("Minuit2");
32 
33  ClusterSizeScan("eta0", "0<|#eta|<0.1");
34  // ClusterSizeScan("eta0.60","0.6<|#eta|<0.7");
35  // PIDScan();
36  //EtaScan();
37 }
38 
39 void EtaScan(void)
40 {
41  const TString config = "Single #pi^{-}, Track-based 5x5 clusterizer";
42 
43  lin_res *QA_Eta0 =
44  GetResolution("pi-", "eta0", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
45  QA_Eta0->name = "0<|eta|<0.1";
46 
47  lin_res *QA_Eta3 =
48  GetResolution("pi-", "eta0.30", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
49  QA_Eta3->name = "0.3<|eta|<0.4";
50 
51  lin_res *QA_Eta6 =
52  GetResolution("pi-", "eta0.60", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
53  QA_Eta6->name = "0.6<|eta|<0.7";
54 
55  lin_res *QA_Eta9 =
56  GetResolution("pi-", "eta0.90", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
57  QA_Eta9->name = "0.9<|eta|<1.0";
58 
59  // ------------
60  QA_Eta0->resolution->SetMarkerStyle(kFullCircle);
61  QA_Eta0->resolution->SetMarkerColor(kRed + 3);
62  QA_Eta0->resolution->SetLineColor(kRed + 3);
63 
64  QA_Eta0->linearity->SetMarkerStyle(kFullCircle);
65  QA_Eta0->linearity->SetMarkerColor(kRed + 3);
66  QA_Eta0->linearity->SetLineColor(kRed + 3);
67 
68  QA_Eta0->f_res->SetLineColor(kRed + 3);
69 
70  // ------------
71  QA_Eta3->resolution->SetMarkerStyle(kFullSquare);
72  QA_Eta3->resolution->SetMarkerColor(kBlue + 3);
73  QA_Eta3->resolution->SetLineColor(kBlue + 3);
74 
75  QA_Eta3->linearity->SetMarkerStyle(kFullSquare);
76  QA_Eta3->linearity->SetMarkerColor(kBlue + 3);
77  QA_Eta3->linearity->SetLineColor(kBlue + 3);
78 
79  QA_Eta3->f_res->SetLineColor(kBlue + 3);
80 
81  // ------------
82  QA_Eta6->resolution->SetMarkerStyle(kFullCross);
83  QA_Eta6->resolution->SetMarkerColor(kMagenta + 3);
84  QA_Eta6->resolution->SetLineColor(kMagenta + 3);
85 
86  QA_Eta6->linearity->SetMarkerStyle(kFullCross);
87  QA_Eta6->linearity->SetMarkerColor(kMagenta + 3);
88  QA_Eta6->linearity->SetLineColor(kMagenta + 3);
89 
90  QA_Eta6->f_res->SetLineColor(kMagenta + 3);
91  QA_Eta6->f_res->SetLineStyle(kDashed);
92 
93  // ------------
94  QA_Eta9->resolution->SetMarkerStyle(kStar);
95  QA_Eta9->resolution->SetMarkerColor(kCyan + 3);
96  QA_Eta9->resolution->SetLineColor(kCyan + 3);
97 
98  QA_Eta9->linearity->SetMarkerStyle(kStar);
99  QA_Eta9->linearity->SetMarkerColor(kCyan + 3);
100  QA_Eta9->linearity->SetLineColor(kCyan + 3);
101 
102  QA_Eta9->f_res->SetLineColor(kCyan + 3);
103  QA_Eta9->f_res->SetLineStyle(kDashed);
104 
105  TCanvas *c1 = new TCanvas(Form("EtaScan"),
106  Form("EtaScan"), 1300, 600);
107  c1->Divide(2, 1);
108  int idx = 1;
109  TPad *p;
110 
111  p = (TPad *) c1->cd(idx++);
112  c1->Update();
113  // p->SetLogy();
114 
115  p->DrawFrame(0, 0, 60, 1.2,
116  Form("Linearity;Truth energy [GeV];Measured Energy / Truth energy"));
117  TLine *l = new TLine(0, 1, 60, 1);
118  l->SetLineWidth(2);
119  l->SetLineColor(kGray);
120  l->Draw();
121 
122  QA_Eta0->linearity->Draw("ep");
123  QA_Eta3->linearity->Draw("ep");
124  QA_Eta6->linearity->Draw("ep");
125  QA_Eta9->linearity->Draw("ep");
126 
127  TLegend *leg = new TLegend(.2, .2, .85, .45);
128  leg->SetHeader("#splitline{#it{#bf{sPHENIX}} Geant4 Simulation}{" + config + "}");
129  leg->AddEntry("", "", "");
130  leg->AddEntry(QA_Eta0->linearity, QA_Eta0->name, "ep");
131  leg->AddEntry(QA_Eta3->linearity, QA_Eta3->name, "ep");
132  leg->AddEntry(QA_Eta6->linearity, QA_Eta6->name, "ep");
133  leg->AddEntry(QA_Eta9->linearity, QA_Eta9->name, "ep");
134  leg->Draw();
135 
136  p = (TPad *) c1->cd(idx++);
137  c1->Update();
138 
139  TF1 *f_calo_t1044 = new TF1("f_calo_t1044", "sqrt([0]*[0]+[1]*[1]/x)/100", 6,
140  32);
141  f_calo_t1044->SetParameters(13.5, 64.9);
142  f_calo_t1044->SetLineWidth(3);
143  f_calo_t1044->SetLineColor(kGreen + 2);
144 
145  TH1 *hframe = p->DrawFrame(0, 0, 60, 0.7,
146  Form("Resolution;Truth energy [GeV];#DeltaE/<E>"));
147 
148  QA_Eta0->f_res->Draw("same");
149  QA_Eta0->resolution->Draw("ep");
150 
151  QA_Eta3->f_res->Draw("same");
152  QA_Eta3->resolution->Draw("ep");
153 
154  QA_Eta6->f_res->Draw("same");
155  QA_Eta6->resolution->Draw("ep");
156 
157  QA_Eta9->f_res->Draw("same");
158  QA_Eta9->resolution->Draw("ep");
159 
160  f_calo_t1044->Draw("same");
161 
162  TLegend *leg = new TLegend(.2, .6, .85, .9);
163  leg->SetHeader("#splitline{#it{#bf{sPHENIX}} Geant4 Simulation}{" + config + "}");
164  leg->AddEntry("", "", "");
165  leg->AddEntry(QA_Eta0->linearity,
166  QA_Eta0->name + " " +
167  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
168  QA_Eta0->f_res->GetParameter(0),
169  QA_Eta0->f_res->GetParameter(1)),
170  "lpe");
171  leg->AddEntry(f_calo_t1044,
172  TString("T-1044 (#eta=0)") + " " +
173  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
174  f_calo_t1044->GetParameter(0),
175  f_calo_t1044->GetParameter(1)),
176  "l");
177  leg->AddEntry(QA_Eta3->linearity,
178  QA_Eta3->name + " " +
179  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
180  QA_Eta3->f_res->GetParameter(0),
181  QA_Eta3->f_res->GetParameter(1)),
182  "lpe");
183  leg->AddEntry(QA_Eta6->linearity,
184  QA_Eta6->name + " " +
185  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
186  QA_Eta6->f_res->GetParameter(0),
187  QA_Eta6->f_res->GetParameter(1)),
188  "lpe");
189  leg->AddEntry(QA_Eta9->linearity,
190  QA_Eta9->name + " " +
191  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
192  QA_Eta9->f_res->GetParameter(0),
193  QA_Eta9->f_res->GetParameter(1)),
194  "lpe");
195  leg->Draw();
196 
197  SaveCanvas(c1,
198  base_dataset + "DrawHadronShowers_" + TString(c1->GetName()), kTRUE);
199 }
200 
201 void PIDScan(void)
202 {
203  const TString config = "Single Particles, 0<|#eta|<0.1, Track-based 5x5 clusterizer";
204 
205  lin_res *QA_electron =
206  GetResolution("e-", "eta0", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
207  QA_electron->name = "Electron";
208 
209  lin_res *QA_PiNeg =
210  GetResolution("pi-", "eta0", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
211  QA_PiNeg->name = "#pi^{-}";
212 
213  lin_res *QA_PiPos =
214  GetResolution("pi\\+", "eta0", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
215  QA_PiPos->name = "#pi^{+}";
216 
217  lin_res *QA_Proton =
218  GetResolution("proton", "eta0", "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
219  QA_Proton->name = "Proton";
220 
221  // ------------
222  QA_electron->resolution->SetMarkerStyle(kFullCircle);
223  QA_electron->resolution->SetMarkerColor(kRed + 3);
224  QA_electron->resolution->SetLineColor(kRed + 3);
225 
226  QA_electron->linearity->SetMarkerStyle(kFullCircle);
227  QA_electron->linearity->SetMarkerColor(kRed + 3);
228  QA_electron->linearity->SetLineColor(kRed + 3);
229 
230  QA_electron->f_res->SetLineColor(kRed + 3);
231 
232  // ------------
233  QA_PiNeg->resolution->SetMarkerStyle(kFullSquare);
234  QA_PiNeg->resolution->SetMarkerColor(kBlue + 3);
235  QA_PiNeg->resolution->SetLineColor(kBlue + 3);
236 
237  QA_PiNeg->linearity->SetMarkerStyle(kFullSquare);
238  QA_PiNeg->linearity->SetMarkerColor(kBlue + 3);
239  QA_PiNeg->linearity->SetLineColor(kBlue + 3);
240 
241  QA_PiNeg->f_res->SetLineColor(kBlue + 3);
242 
243  // ------------
244  QA_PiPos->resolution->SetMarkerStyle(kFullCross);
245  QA_PiPos->resolution->SetMarkerColor(kMagenta + 3);
246  QA_PiPos->resolution->SetLineColor(kMagenta + 3);
247 
248  QA_PiPos->linearity->SetMarkerStyle(kFullCross);
249  QA_PiPos->linearity->SetMarkerColor(kMagenta + 3);
250  QA_PiPos->linearity->SetLineColor(kMagenta + 3);
251 
252  QA_PiPos->f_res->SetLineColor(kMagenta + 3);
253  QA_PiPos->f_res->SetLineStyle(kDashed);
254 
255  // ------------
256  QA_Proton->resolution->SetMarkerStyle(kStar);
257  QA_Proton->resolution->SetMarkerColor(kCyan + 3);
258  QA_Proton->resolution->SetLineColor(kCyan + 3);
259 
260  QA_Proton->linearity->SetMarkerStyle(kStar);
261  QA_Proton->linearity->SetMarkerColor(kCyan + 3);
262  QA_Proton->linearity->SetLineColor(kCyan + 3);
263 
264  QA_Proton->f_res->SetLineColor(kCyan + 3);
265  QA_Proton->f_res->SetLineStyle(kDashed);
266 
267  TCanvas *c1 = new TCanvas(Form("PIDScan"),
268  Form("PIDScan"), 1300, 600);
269  c1->Divide(2, 1);
270  int idx = 1;
271  TPad *p;
272 
273  p = (TPad *) c1->cd(idx++);
274  c1->Update();
275  // p->SetLogy();
276 
277  p->DrawFrame(0, 0, 60, 1.2,
278  Form("Linearity;Truth energy [GeV];Measured Energy / Truth energy"));
279  TLine *l = new TLine(0, 1, 60, 1);
280  l->SetLineWidth(2);
281  l->SetLineColor(kGray);
282 
283  l->Draw();
284 
285  QA_electron->linearity->Draw("ep");
286  QA_PiNeg->linearity->Draw("ep");
287  QA_PiPos->linearity->Draw("ep");
288  QA_Proton->linearity->Draw("ep");
289 
290  TLegend *leg = new TLegend(.2, .2, 1, .45);
291  leg->SetHeader("#splitline{#it{#bf{sPHENIX}} Geant4 Simulation}{" + config + "}");
292  leg->AddEntry("", "", "");
293  leg->AddEntry(QA_electron->linearity, QA_electron->name, "ep");
294  leg->AddEntry(QA_PiNeg->linearity, QA_PiNeg->name, "ep");
295  leg->AddEntry(QA_PiPos->linearity, QA_PiPos->name, "ep");
296  leg->AddEntry(QA_Proton->linearity, QA_Proton->name, "ep");
297  leg->Draw();
298 
299  p = (TPad *) c1->cd(idx++);
300  c1->Update();
301 
302  TF1 *f_calo_t1044 = new TF1("f_calo_t1044", "sqrt([0]*[0]+[1]*[1]/x)/100", 6,
303  32);
304  f_calo_t1044->SetParameters(13.5, 64.9);
305  f_calo_t1044->SetLineWidth(3);
306  f_calo_t1044->SetLineColor(kGreen + 2);
307 
308  TH1 *hframe = p->DrawFrame(0, 0, 60, 0.7,
309  Form("Resolution;Truth energy [GeV];#DeltaE/<E>"));
310 
311  QA_electron->f_res->Draw("same");
312  QA_electron->resolution->Draw("ep");
313 
314  QA_PiNeg->f_res->Draw("same");
315  QA_PiNeg->resolution->Draw("ep");
316 
317  QA_PiPos->f_res->Draw("same");
318  QA_PiPos->resolution->Draw("ep");
319 
320  QA_Proton->f_res->Draw("same");
321  QA_Proton->resolution->Draw("ep");
322 
323  f_calo_t1044->Draw("same");
324 
325  TLegend *leg = new TLegend(.2, .6, 1, .9);
326  leg->SetHeader("#splitline{#it{#bf{sPHENIX}} Geant4 Simulation}{" + config + "}");
327  leg->AddEntry("", "", "");
328  leg->AddEntry(QA_electron->linearity,
329  QA_electron->name + " " +
330  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
331  QA_electron->f_res->GetParameter(0),
332  QA_electron->f_res->GetParameter(1)),
333  "lpe");
334  leg->AddEntry(f_calo_t1044,
335  TString("T-1044 Data, #pi^{-}") + " " +
336  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
337  f_calo_t1044->GetParameter(0),
338  f_calo_t1044->GetParameter(1)),
339  "l");
340  leg->AddEntry(QA_PiNeg->linearity,
341  QA_PiNeg->name + " " +
342  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
343  QA_PiNeg->f_res->GetParameter(0),
344  QA_PiNeg->f_res->GetParameter(1)),
345  "lpe");
346  leg->AddEntry(QA_PiPos->linearity,
347  QA_PiPos->name + " " +
348  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
349  QA_PiPos->f_res->GetParameter(0),
350  QA_PiPos->f_res->GetParameter(1)),
351  "lpe");
352  leg->AddEntry(QA_Proton->linearity,
353  QA_Proton->name + " " +
354  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
355  QA_Proton->f_res->GetParameter(0),
356  QA_Proton->f_res->GetParameter(1)),
357  "lpe");
358  leg->Draw();
359 
360  SaveCanvas(c1,
361  base_dataset + "DrawHadronShowers_" + TString(c1->GetName()), kTRUE);
362 }
363 
364 void ClusterSizeScan(TString seta = "eta0", TString eta_desc = "0<|#eta|<0.1")
365 {
366  const TString config = "Single #pi^{-}, " + eta_desc;
367 
368  lin_res *QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP =
369  GetResolution("pi-", seta, "h_QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP");
370  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->name = "3x3 cluster";
371 
372  lin_res *QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP =
373  GetResolution("pi-", seta, "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP");
374  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->name = "5x5 cluster";
375 
376  lin_res *QAG4SimJet_AntiKt_Tower_r02_Leading_Et =
377  GetResolution("pi-", seta, "h_QAG4SimJet_AntiKt_Tower_r02_Leading_Et");
378  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->name = "Anti-k_{T} R=0.2";
379 
380  lin_res *QAG4SimJet_AntiKt_Tower_r04_Leading_Et =
381  GetResolution("pi-", seta, "h_QAG4SimJet_AntiKt_Tower_r04_Leading_Et");
382  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->name = "Anti-k_{T} R=0.4";
383 
384  // ------------
385  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->resolution->SetMarkerStyle(kFullCircle);
386  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->resolution->SetMarkerColor(kRed + 3);
387  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->resolution->SetLineColor(kRed + 3);
388 
389  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->linearity->SetMarkerStyle(kFullCircle);
390  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->linearity->SetMarkerColor(kRed + 3);
391  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->linearity->SetLineColor(kRed + 3);
392 
393  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->f_res->SetLineColor(kRed + 3);
394 
395  // ------------
396  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->resolution->SetMarkerStyle(kFullSquare);
397  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->resolution->SetMarkerColor(kBlue + 3);
398  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->resolution->SetLineColor(kBlue + 3);
399 
400  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->linearity->SetMarkerStyle(kFullSquare);
401  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->linearity->SetMarkerColor(kBlue + 3);
402  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->linearity->SetLineColor(kBlue + 3);
403 
404  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->f_res->SetLineColor(kBlue + 3);
405 
406  // ------------
407  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->resolution->SetMarkerStyle(kFullCross);
408  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->resolution->SetMarkerColor(kMagenta + 3);
409  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->resolution->SetLineColor(kMagenta + 3);
410 
411  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->linearity->SetMarkerStyle(kFullCross);
412  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->linearity->SetMarkerColor(kMagenta + 3);
413  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->linearity->SetLineColor(kMagenta + 3);
414 
415  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->f_res->SetLineColor(kMagenta + 3);
416  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->f_res->SetLineStyle(kDashed);
417 
418  // ------------
419  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->resolution->SetMarkerStyle(kStar);
420  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->resolution->SetMarkerColor(kCyan + 3);
421  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->resolution->SetLineColor(kCyan + 3);
422 
423  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->linearity->SetMarkerStyle(kStar);
424  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->linearity->SetMarkerColor(kCyan + 3);
425  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->linearity->SetLineColor(kCyan + 3);
426 
427  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->f_res->SetLineColor(kCyan + 3);
428  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->f_res->SetLineStyle(kDashed);
429 
430  TCanvas *c1 = new TCanvas(Form("ClusterSizeScan_") + seta,
431  Form("ClusterSizeScan_") + seta, 1300, 600);
432  c1->Divide(2, 1);
433  int idx = 1;
434  TPad *p;
435 
436  p = (TPad *) c1->cd(idx++);
437  c1->Update();
438  // p->SetLogy();
439 
440  p->DrawFrame(0, 0, 60, 1.2,
441  Form("Linearity;Truth energy [GeV];Measured Energy / Truth energy"));
442  TLine *l = new TLine(0, 1, 60, 1);
443  l->SetLineWidth(2);
444  l->SetLineColor(kGray);
445  l->Draw();
446 
447  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->linearity->Draw("ep");
448  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->linearity->Draw("ep");
449  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->linearity->Draw("ep");
450  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->linearity->Draw("ep");
451 
452  TLegend *leg = new TLegend(.2, .2, .85, .45);
453  leg->SetHeader("#splitline{#it{#bf{sPHENIX}} Geant4 Simulation}{" + config + "}");
454  leg->AddEntry("", "", "");
455  leg->AddEntry(QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->linearity, QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->name, "ep");
456  leg->AddEntry(QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->linearity, QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->name, "ep");
457  leg->AddEntry(QAG4SimJet_AntiKt_Tower_r02_Leading_Et->linearity, QAG4SimJet_AntiKt_Tower_r02_Leading_Et->name, "ep");
458  leg->AddEntry(QAG4SimJet_AntiKt_Tower_r04_Leading_Et->linearity, QAG4SimJet_AntiKt_Tower_r04_Leading_Et->name, "ep");
459  leg->Draw();
460 
461  p = (TPad *) c1->cd(idx++);
462  c1->Update();
463 
464  TF1 *f_calo_t1044 = new TF1("f_calo_t1044", "sqrt([0]*[0]+[1]*[1]/x)/100", 6,
465  32);
466  f_calo_t1044->SetParameters(13.5, 64.9);
467  f_calo_t1044->SetLineWidth(3);
468  f_calo_t1044->SetLineColor(kGreen + 2);
469 
470  TH1 *hframe = p->DrawFrame(0, 0, 60, 0.7,
471  Form("Resolution;Truth energy [GeV];#DeltaE/<E>"));
472 
473  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->f_res->Draw("same");
474  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->resolution->Draw("ep");
475 
476  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->f_res->Draw("same");
477  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->resolution->Draw("ep");
478 
479  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->f_res->Draw("same");
480  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->resolution->Draw("ep");
481 
482  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->f_res->Draw("same");
483  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->resolution->Draw("ep");
484 
485  f_calo_t1044->Draw("same");
486 
487  TLegend *leg = new TLegend(.2, .6, .85, .9);
488  leg->SetHeader("#splitline{#it{#bf{sPHENIX}} Geant4 Simulation}{" + config + "}");
489  leg->AddEntry("", "", "");
490  leg->AddEntry(QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->linearity,
491  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->name + " " +
492  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
493  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->f_res->GetParameter(0),
494  QAG4Sim_CalorimeterSum_TrackProj_3x3Tower_EP->f_res->GetParameter(1)),
495  "lpe");
496  leg->AddEntry(f_calo_t1044,
497  TString("T-1044 (4x4, #eta=0)") + " " +
498  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
499  f_calo_t1044->GetParameter(0),
500  f_calo_t1044->GetParameter(1)),
501  "l");
502  leg->AddEntry(QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->linearity,
503  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->name + " " +
504  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
505  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->f_res->GetParameter(0),
506  QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP->f_res->GetParameter(1)),
507  "lpe");
508  leg->AddEntry(QAG4SimJet_AntiKt_Tower_r02_Leading_Et->linearity,
509  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->name + " " +
510  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
511  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->f_res->GetParameter(0),
512  QAG4SimJet_AntiKt_Tower_r02_Leading_Et->f_res->GetParameter(1)),
513  "lpe");
514  leg->AddEntry(QAG4SimJet_AntiKt_Tower_r04_Leading_Et->linearity,
515  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->name + " " +
516  Form("%.1f%% #oplus %.1f%%/#sqrt{E}",
517  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->f_res->GetParameter(0),
518  QAG4SimJet_AntiKt_Tower_r04_Leading_Et->f_res->GetParameter(1)),
519  "lpe");
520  leg->Draw();
521 
522  SaveCanvas(c1,
523  base_dataset + "DrawHadronShowers_" + TString(c1->GetName()), kTRUE);
524 }
525 
526 lin_res *
527 GetResolution(TString PID = "pi-", TString eta = "eta0", TString qa_histo = "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP")
528 {
529  TString dataset_name = PID + "_" + eta + "_" + qa_histo;
530 
531  const int N = 6;
532  const double Ps[] = {8, 12, 16, 32, 40, 50, -1, -1};
533 
534  vector<double> mean;
535  vector<double> mean_err;
536  vector<double> res;
537  vector<double> res_err;
538 
539  TCanvas *c1 = new TCanvas("GetResolution_LineShape_" + dataset_name,
540  "GetResolution_LineShape_" + dataset_name, 1800, 900);
541 
542  c1->Divide(3, 2);
543  int idx = 1;
544  TPad *p;
545 
546  // for (int i = 0; i < N; ++i)
547  for (int i = 2; i < 3; ++i) //for 16GeV
548  {
549  // if(i==2){ //for 16GeV
550 
551  p = (TPad *) c1->cd(idx++);
552  c1->Update();
553 
554  const double P = Ps[i];
555 
556  TString sEnergy = Form("%.0fGeV", P);
557 
558  TH1 *h = NULL;
559 
560  //if (qa_histo.Length() == 0)
561  h = GetEvalHisto(PID, eta, sEnergy);
562  // else
563  // h = GetQAHisto(PID, eta, sEnergy, qa_histo);
564 
565  TF1 *fgaus_g = new TF1("fgaus_LG_g_" + dataset_name, "gaus",
566  h->GetMean() - 4 * h->GetRMS(), h->GetMean() + 4 * h->GetRMS());
567  fgaus_g->SetParameters(1, h->GetMean(),
568  h->GetRMS());
569  h->Fit(fgaus_g, "MR0N");
570  fgaus_g->SetLineColor(kGray);
571 
572  TF1 *fgaus = new TF1("fgaus_LG_" + dataset_name, "gaus",
573  fgaus_g->GetParameter(1) - 1.5 * fgaus_g->GetParameter(2),
574  fgaus_g->GetParameter(1) + 1.5 * fgaus_g->GetParameter(2));
575  fgaus->SetParameters(fgaus_g->GetParameter(0), fgaus_g->GetParameter(1),
576  fgaus_g->GetParameter(2));
577  fgaus->SetLineColor(kBlue);
578  h->Fit(fgaus, "MR0N");
579 
580  h->Draw();
581  fgaus_g->Draw("same");
582  fgaus->Draw("same");
583 
584  double scale = 1;
585  if (qa_histo.Contains("Leading_Et")) scale = 1. / P;
586 
587  mean.push_back(fgaus->GetParameter(1) * scale);
588  mean_err.push_back(fgaus->GetParError(1) * scale);
589  res.push_back(fgaus->GetParameter(2) / fgaus->GetParameter(1));
590  res_err.push_back(fgaus->GetParError(2) / fgaus->GetParameter(1));
591 
592 
593  //}// for 16GeV
594  }
595 
596  TGraphErrors *ge_linear = new TGraphErrors(N, Ps,
597  &mean[0], 0, &mean_err[0]);
598 
599  TGraphErrors *ge_res = new TGraphErrors(N, Ps,
600  &res[0], 0, &res_err[0]);
601  ge_res->GetHistogram()->SetStats(0);
602  ge_res->Print();
603 
604  lin_res *ret = new lin_res;
605  ret->name = dataset_name;
606  ret->linearity = ge_linear;
607  ret->resolution = ge_res;
608  ret->f_res = new TF1("f_calo_r_" + dataset_name, "sqrt([0]*[0]+[1]*[1]/x)/100",
609  6, 60);
610  ge_res->Fit(ret->f_res, "RM0N");
611 
612  SaveCanvas(c1,
613  base_dataset + "DrawHadronShowers_" + TString(c1->GetName()), kTRUE);
614 
615  new TCanvas;
616 
617  ret->f_res->Draw();
618  ge_res->Draw("*e");
619 
620  return ret;
621 }
622 
623 TH1 *GetEvalHisto(TString PID = "pi-", TString eta = "eta0", TString energy = "16GeV")
624 {
625  TString dataset_name = PID + "_" + eta + "_" + energy;
626  TH1 *h_E = new TH1F("h_ESumOverp_" + dataset_name, dataset_name + " h_ESumOverp;E/p;Count per bin", 100, 0, 2);
627  //TString infile = base_dataset + "G4Hits_sPHENIX_" + dataset_name + "_*.root_g4svtx_eval.root";
628  TString infile = base_dataset + "G4Hits_sPHENIX_" + dataset_name + ".root_g4svtx_eval.root";
629 
630  TChain *ntp_track = new TChain("ntp_track");
631  const int n = ntp_track->Add(infile);
632 
633  cout << "GetEvalHisto - loaded " << n << " files from " << infile << endl;
634 
635  ntp_track->Draw("(cemce3x3 + hcaloute3x3 + hcaline3x3)/sqrt(gpx**2+gpy**2+gpz**2)>>h_ESumOverp_" + dataset_name,
636  "nfromtruth>50 && gtrackID == 1", "goff");
637 
638  delete ntp_track;
639 
640  return h_E;
641 }
642 
643 TH1 *GetQAHisto(TString PID = "pi-", TString eta = "eta0", TString energy = "16GeV", TString qa_histo = "h_QAG4Sim_CalorimeterSum_TrackProj_5x5Tower_EP")
644 {
645  cout << "Weihu Ma run test " << endl;
646 
647  TString dataset_name = PID + "_" + eta + "_" + energy;
648  TH1 *h = NULL;
649  //TString infile = base_dataset + "G4Hits_sPHENIX_" + dataset_name + "_*.root_qa.root";
650  TString infile = base_dataset + "G4Hits_sPHENIX_" + dataset_name + ".root_qa.root";
651  cout << "GetQAHisto - searching " << infile << endl;
652 
653  TChain *T = new TChain("T");
654  const int n = T->Add(infile);
655 
656  TObjArray *fileElements = T->GetListOfFiles();
657  TIter next(fileElements);
658  TChainElement *chEl = 0;
659  while ((chEl = (TChainElement *) next()))
660  {
661  cout << "GetQAHisto - processing " << chEl->GetTitle() << endl;
662 
663  TFile f(chEl->GetTitle());
664 
665  if (f.IsOpen())
666  {
667  TH1 *hqa = (TH1 *) f.GetObjectChecked(qa_histo, "TH1");
668  cout << "Weihu Ma run test1 " << endl;
669  if (hqa)
670  {
671  cout << "Weihu Ma run test4 " << endl;
672  if (h)
673  {
674  cout << "Weihu Ma run test2 " << endl;
675  h->Add(hqa);
676  }
677  else
678  {
679  cout << "Weihu Ma run test3 " << endl;
680  h = (TH1 *) (hqa->Clone(TString(hqa->GetName()) + "_" + dataset_name));
681  h->SetDirectory(NULL);
682  assert(h);
683  }
684  }
685  }
686  }
687  //assert(h);
688 
689  return h;
690 }