Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DrawPrototype4ShowerCalib.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DrawPrototype4ShowerCalib.C
1 // $Id: $
2 
11 #include <TFile.h>
12 #include <TGraphErrors.h>
13 #include <TLatex.h>
14 #include <TLine.h>
15 #include <TString.h>
16 #include <TTree.h>
17 #include <cassert>
18 #include <cmath>
19 #include "SaveCanvas.C"
20 #include "SetOKStyle.C"
21 using namespace std;
22 
23 //#include "Prototype3_DSTReader.h"
24 
25 TCut event_sel;
26 TString cuts;
27 TFile *_file0 = NULL;
28 TTree *T = NULL;
29 
30 class lin_res
31 {
32  public:
33  TString name;
34  TGraphErrors *linearity;
35  TGraphErrors *resolution;
36  TF1 *f_res;
37 };
38 
39 void DrawPrototype4ShowerCalib( //
40  const TString infile =
41  // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2018/ShowerCalib/dst.lst_EMCalCalib.root" //
42  // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2018/ShowerCalib_tilted/dst.lst_EMCalCalib.root"//
43  // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2018/Scan1Block36/dst.lst_EMCalCalib.root" //
44  // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2018/Scan2Block34/dst.lst_EMCalCalib.root" //
45  // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2018/Scan2Block18/dst.lst_EMCalCalib.root" //
46  // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2018/Scan64.28V/dst.lst_EMCalCalib.root" //
47  // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2018/Scan4Block45/dst.lst_EMCalCalib.root" //
48  "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2018/Scan2018b1Tower36/dst.lst_EMCalCalib.root" //
49  )
50 {
51  SetOKStyle();
52  gStyle->SetOptStat(0);
53  gStyle->SetOptFit(1111);
54  gStyle->SetPadGridX(0);
55  gStyle->SetPadGridY(0);
56  TVirtualFitter::SetDefaultFitter("Minuit2");
57 
58  gSystem->Load("libPrototype4.so");
59  gSystem->Load("libProto4ShowCalib.so");
60 
61  // gROOT->LoadMacro("Prototype3_DSTReader.C+");
62 
63  if (!_file0)
64  {
65  TString chian_str = infile;
66  chian_str.ReplaceAll("ALL", "*");
67 
68  TChain *t = new TChain("T");
69  const int n = t->Add(chian_str);
70 
71  cout << "Loaded " << n << " root files with " << chian_str << endl;
72  assert(n > 0);
73 
74  T = t;
75 
76  _file0 = new TFile;
77  _file0->SetName(infile);
78  }
79 
80  assert(_file0);
81 
82  // event_sel = "1";
83  // cuts = "_all_data";
84 // event_sel = "good_e";
85 // cuts = "_good_e";
86 
87  // event_sel = "info.beam_mom == -8 && good_e";
88  // cuts = "_8GeV_good_e";
89 // event_sel = "info.beam_mom == -6 && good_e";
90 // cuts = "_6GeV_good_e";
91  // event_sel = " good_e && info.run == 409";
92  // cuts = "_good_e_run409";
93  // event_sel = "info.beam_mom == -12 && good_e";
94  // cuts = "_12GeV_good_e";
95  // event_sel = "info.beam_mom == -16 && good_e";
96  // cuts = "_16GeV_good_e";
97  // event_sel = "info.beam_mom == -6";
98  // cuts = "_Neg6GeV";
99 
100  // event_sel = "good_e && info.hodo_h==3 && info.hodo_v==5"; // Tower 36
101  // cuts = "_good_e_h3_v5";
102  // event_sel = "good_e && info.hodo_h>=2 && info.hodo_h<=4 && info.hodo_v>=4 && info.hodo_v<=6"; // Tower 36
103  // cuts = "_good_e_h234_v456";
104  // event_sel = "info.beam_mom == -8 && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h==3 && info.hodo_v==5"; // Tower 36
105  // cuts = "_valid_data_h3_v5_8GeV"; // Tower 36
106  // event_sel = "info.beam_mom == -2 && valid_hodo_v && valid_hodo_h&& trigger_veto_pass && info.hodo_h>=2 && info.hodo_h<=4 && info.hodo_v>=4 && info.hodo_v<=6"; // Tower 36
107  // cuts = "_valid_data_h234_v456_2GeV"; // Tower 36
108 
109  // event_sel = "good_e && info.hodo_h==3 && info.hodo_v==4"; // Tower 34/18
110  // cuts = "_good_e_h3_v4";
111  // event_sel = "good_e && info.hodo_h>=2 && info.hodo_h<=4 && info.hodo_v>=3 && info.hodo_v<=5"; // Tower 34/18
112  // cuts = "_good_e_h234_v345";
113  // event_sel = "good_e && info.hodo_h>=1 && info.hodo_h<=5 && info.hodo_v>=2 && info.hodo_v<=6"; // Tower 34/18
114  // cuts = "_good_e_h12345_v23456";
115 
116  event_sel = "good_e && info.hodo_h==3 && info.hodo_v==3"; // Tower 34/18 - 2018b
117  cuts = "_good_e_h3_v3";
118  // event_sel = "good_e && info.hodo_h>=3 && info.hodo_h<=4 && info.hodo_v>=2 && info.hodo_v<=4"; // Tower 34/18
119  // cuts = "_good_e_h34_v234";
120  // event_sel = "good_e && info.hodo_h>=1 && info.hodo_h<=5 && info.hodo_v>=2 && info.hodo_v<=6"; // Tower 34/18
121  // cuts = "_good_e_h12345_v23456";
122 
123  // event_sel = "good_e && info.hodo_h==5 && info.hodo_v==2"; // Tower 34/18 between towers
124  // cuts = "_good_e_h5_v2";
125 
126  T->SetAlias("SimEnergyScale", "1*1");
127  // // based on /phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/UIUC21.lst_EMCalCalib.root_DrawPrototype3ShowerCalib_LineShapeData_Neg8GeV_good_data_h5_v3.svg
128  // T->SetAlias("SimEnergyScale","8.74635e+00/7.60551");
129  // // based on /phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/ShowerCalib/Tilt0.lst_EMCalCalib.root_DrawPrototype3ShowerCalib_LineShapeData_Neg8GeV_quality_h3_v5_col2_row2.root_DrawPrototype3ShowerCalib_SumLineShapeCompare_Electron_8GeV_QGSP_BERT_HP.root
130  // T->SetAlias("SimEnergyScale","8.88178/8.16125e+00");
131  // Tilt0.lst <-> QGSP_BERT_HP Birk 0.151
132  // T->SetAlias("SimEnergyScale","8.88178/8.01845");
133  // Tilt0.lst <-> QGSP_BERT_HP Birk 0.18
134  // T->SetAlias("SimEnergyScale","8.88178/7.94838");
135 
136  // Dummy
137  // T->SetAlias("SimEnergyScale","1*1");
138 
139  cout << "Build event selection of " << (const char *) event_sel << endl;
140 
141  T->Draw(">>EventList", event_sel);
142  TEventList *elist = gDirectory->GetObjectChecked("EventList", "TEventList");
143  cout << elist->GetN() << " / " << T->GetEntriesFast() << " events selected"
144  << endl;
145 
146  T->SetEventList(elist);
147 
148  // // data stuff
149  PositionDependenceData("clus_5x5_prod.sum_E");
150  // PositionDependenceData("clus_5x5_recalib.sum_E");
151  HodoscopeCheck();
152 
153  // LineShapeData("abs(info.C2_sum)<200", "(info.C2_sum)>2000");
154 
155  Get_Res_linear_Summmary("sum_E*.55", 30);
156  // Get_Res_linear_Summmary("sum_E*.13", 10);
157 
158  // simulation stuff
159  // SimPositionCheck(-0); // 0 degree tilted
160  // LineShapeSim();
161 
162  // PositionDependenceSim("clus_5x5_prod.sum_E", -0, 5); // 0 degree tilted
163  // SimPositionCheck(-15); // 10 degree tilted
164  // PositionDependenceSim("clus_5x5_prod.sum_E", -15, 5); // 10 degree tilted
165  // SimPositionCheck(-40+3); // 45 degree tilted
166  // Get_Res_linear_Summmary_Sim();
167 }
168 
169 void PositionDependenceData(TString sTOWER = "clus_5x5_prod.sum_E",
170  const double z_shift = 0, const int n_div = 1)
171 {
172  TH3F *EnergySum_LG3 =
173  new TH3F("EnergySum_LG3",
174  ";Horizontal Hodoscope (5 mm);Vertical Hodoscope (5 mm);5x5 Cluster Energy (GeV)", //
175  8, -.5, 7.5, //
176  8, -.5, 7.5, //
177  200, 1, 20);
178 
179  T->Draw(sTOWER + ":7-hodo_v:hodo_h>>EnergySum_LG3", "", "goff");
180 
181  TProfile2D *EnergySum_LG3_prof_xy = EnergySum_LG3->Project3DProfile("yx");
182  TH2 *EnergySum_LG3_yx = EnergySum_LG3->Project3D("yx");
183  TH2 *EnergySum_LG3_zx = EnergySum_LG3->Project3D("zx");
184  TH2 *EnergySum_LG3_zy = EnergySum_LG3->Project3D("zy");
185 
186  TGraphErrors *ge_EnergySum_LG3_zx = FitProfile(EnergySum_LG3_zx);
187  TGraphErrors *ge_EnergySum_LG3_zy = FitProfile(EnergySum_LG3_zy);
188 
189  TText *t;
190  TCanvas *c1 = new TCanvas(
191  TString(Form("EMCDistributionVSBeam_SUM_NDiv%d_", n_div)) + sTOWER + cuts,
192  TString(Form("EMCDistributionVSBeam_SUM_NDiv%d_", n_div)) + sTOWER + cuts,
193  1000, 960);
194  c1->Divide(2, 2);
195  int idx = 1;
196  TPad *p;
197 
198  p = (TPad *) c1->cd(idx++);
199  c1->Update();
200  // p->SetLogy();
201  p->SetGridx(0);
202  p->SetGridy(0);
203 
204  EnergySum_LG3_prof_xy->SetMinimum(0);
205 
206  EnergySum_LG3_prof_xy->Draw("colz");
207  EnergySum_LG3_prof_xy->SetTitle(
208  "Energy response;Horizontal Hodoscope (5 mm);7 - Vertical Hodoscope (5 mm)");
209 
210  p = (TPad *) c1->cd(idx++);
211  c1->Update();
212  // p->SetLogy();
213  p->SetGridx(0);
214  p->SetGridy(0);
215 
216  EnergySum_LG3_yx->SetMinimum(0);
217  EnergySum_LG3_yx->Draw("colz");
218  EnergySum_LG3_yx->SetTitle(
219  "Event counts;Horizontal Hodoscope (5 mm);7 - Vertical Hodoscope (5 mm)");
220 
221  p = (TPad *) c1->cd(idx++);
222  c1->Update();
223  // p->SetLogy();
224  p->SetGridx(0);
225  p->SetGridy(0);
226 
227  EnergySum_LG3_zx->Draw("colz");
228  EnergySum_LG3_zx->SetTitle(
229  "Position scan;Horizontal Hodoscope (5 mm);5x5 Cluster Energy (GeV)");
230 
231  ge_EnergySum_LG3_zx->SetLineWidth(2);
232  ge_EnergySum_LG3_zx->SetMarkerStyle(kFullCircle);
233  ge_EnergySum_LG3_zx->Draw("pe");
234 
235  p = (TPad *) c1->cd(idx++);
236  c1->Update();
237  // p->SetLogy();
238  p->SetGridx(0);
239  p->SetGridy(0);
240 
241  EnergySum_LG3_zy->Draw("colz");
242  EnergySum_LG3_zy->SetTitle(
243  "Position scan;7 - Vertical Hodoscope (5 mm);5x5 Cluster Energy (GeV)");
244 
245  ge_EnergySum_LG3_zy->SetLineWidth(2);
246  ge_EnergySum_LG3_zy->SetMarkerStyle(kFullCircle);
247  ge_EnergySum_LG3_zy->Draw("pe");
248 
249  // p = (TPad *) c1->cd(idx++);
250  // c1->Update();
252  // p->SetGridx(0);
253  // p->SetGridy(0);
254  //
255  // TH1 * h = (TH1 *) EnergySum_LG3->ProjectionZ();
256  //
257  // TF1 * fgaus = new TF1("fgaus_LG", "gaus", 0, 100);
258  // fgaus->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
259  // h->GetMean() + 2 * h->GetRMS());
260  // h->Fit(fgaus, "M");
261  //
262  // h->Sumw2();
263  // h->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
264  // h->GetMean() + 4 * h->GetRMS());
265  // EnergySum_LG3_zx->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
266  // h->GetMean() + 4 * h->GetRMS());
267  // EnergySum_LG3_zy->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
268  // h->GetMean() + 4 * h->GetRMS());
269  //
270  // h->SetLineWidth(2);
271  // h->SetMarkerStyle(kFullCircle);
272  //
273  // h->SetTitle(
274  // Form("#DeltaE/<E> = %.1f%%",
275  // 100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
276 
277  SaveCanvas(c1,
278  TString(_file0->GetName()) + TString("_DrawPrototype3ShowerCalib_") + TString(c1->GetName()), kTRUE);
279 }
280 
281 void PositionDependenceSim(TString sTOWER = "clus_5x5_prod.sum_E",
282  const double z_shift = 0, const int n_div = 1)
283 {
284  TH3F *EnergySum_LG3 =
285  new TH3F("EnergySum_LG3",
286  ";Beam Horizontal Pos (cm);Beam Vertical Pos (cm);5x5 Cluster Energy (GeV)", //
287  20 * n_div, z_shift - 5, z_shift + 5, //
288  20 * n_div, -5, 5, //
289  200, 1, 20);
290 
291  T->Draw(sTOWER + ":info.truth_y:info.truth_z>>EnergySum_LG3", "", "goff");
292 
293  TProfile2D *EnergySum_LG3_xy = EnergySum_LG3->Project3DProfile("yx");
294  TH2 *EnergySum_LG3_zx = EnergySum_LG3->Project3D("zx");
295  TH2 *EnergySum_LG3_zy = EnergySum_LG3->Project3D("zy");
296 
297  TGraphErrors *ge_EnergySum_LG3_zx = FitProfile(EnergySum_LG3_zx);
298  TGraphErrors *ge_EnergySum_LG3_zy = FitProfile(EnergySum_LG3_zy);
299 
300  TText *t;
301  TCanvas *c1 = new TCanvas(
302  TString(Form("EMCDistributionVSBeam_SUM_NDiv%d_", n_div)) + sTOWER + cuts,
303  TString(Form("EMCDistributionVSBeam_SUM_NDiv%d_", n_div)) + sTOWER + cuts,
304  1000, 960);
305  c1->Divide(2, 2);
306  int idx = 1;
307  TPad *p;
308 
309  p = (TPad *) c1->cd(idx++);
310  c1->Update();
311  // p->SetLogy();
312  p->SetGridx(0);
313  p->SetGridy(0);
314 
315  EnergySum_LG3_xy->Draw("colz");
316  EnergySum_LG3_xy->SetTitle(
317  "Position scan;Beam Horizontal Pos (cm);Beam Vertical Pos (cm)");
318 
319  p = (TPad *) c1->cd(idx++);
320  c1->Update();
321  // p->SetLogy();
322  p->SetGridx(0);
323  p->SetGridy(0);
324 
325  EnergySum_LG3_zx->Draw("colz");
326  EnergySum_LG3_zx->SetTitle(
327  "Position scan;Beam Horizontal Pos (cm);5x5 Cluster Energy (GeV)");
328 
329  ge_EnergySum_LG3_zx->SetLineWidth(2);
330  ge_EnergySum_LG3_zx->SetMarkerStyle(kFullCircle);
331  ge_EnergySum_LG3_zx->Draw("pe");
332 
333  p = (TPad *) c1->cd(idx++);
334  c1->Update();
335  // p->SetLogy();
336  p->SetGridx(0);
337  p->SetGridy(0);
338 
339  EnergySum_LG3_zy->Draw("colz");
340  EnergySum_LG3_zy->SetTitle(
341  "Position scan;Beam Vertical Pos (cm);5x5 Cluster Energy (GeV)");
342 
343  ge_EnergySum_LG3_zy->SetLineWidth(2);
344  ge_EnergySum_LG3_zy->SetMarkerStyle(kFullCircle);
345  ge_EnergySum_LG3_zy->Draw("pe");
346 
347  p = (TPad *) c1->cd(idx++);
348  c1->Update();
349  // p->SetLogy();
350  p->SetGridx(0);
351  p->SetGridy(0);
352 
353  TH1 *h = (TH1 *) EnergySum_LG3->ProjectionZ();
354 
355  TF1 *fgaus = new TF1("fgaus_LG", "gaus", 0, 100);
356  fgaus->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
357  h->GetMean() + 2 * h->GetRMS());
358  h->Fit(fgaus, "M");
359 
360  h->Sumw2();
361  h->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
362  h->GetMean() + 4 * h->GetRMS());
363  EnergySum_LG3_zx->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
364  h->GetMean() + 4 * h->GetRMS());
365  EnergySum_LG3_zy->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
366  h->GetMean() + 4 * h->GetRMS());
367 
368  h->SetLineWidth(2);
369  h->SetMarkerStyle(kFullCircle);
370 
371  h->SetTitle(
372  Form("#DeltaE/<E> = %.1f%%",
373  100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
374 
375  SaveCanvas(c1,
376  TString(_file0->GetName()) + TString("_DrawPrototype3ShowerCalib_") + TString(c1->GetName()), kTRUE);
377 }
378 
379 void LineShapeData(TCut c2_h = "abs(info.C2_sum)<100", TCut c2_e =
380  "(info.C2_sum)>500 && (info.C2_sum)<1300")
381 {
382  vector<double> mom;
383 
384  TText *t;
385  TCanvas *c1 = new TCanvas("LineShapeData" + cuts, "LineShapeData" + cuts,
386  1800, 950);
387 
388  c1->Divide(2, 2);
389  int idx = 1;
390  TPad *p;
391 
392  p = (TPad *) c1->cd(idx++);
393  c1->Update();
394  // p->SetGridx(1);
395  // p->SetGridy(1);
396  p->SetLogy();
397 
398  T->Draw("info.C2_sum>>h_c2_sum(300,-500,22500)");
399  h_c2_sum->SetTitle("Cherenkov Checks;Sum C2 (ADC)");
400  T->Draw("info.C2_sum>>h_c2_h(300,-500,22500)", c2_h, "same");
401  T->Draw("info.C2_sum>>h_c2_e(300,-500,22500)", c2_e, "same");
402 
403  h_c2_h->SetLineColor(kBlue);
404  h_c2_e->SetLineColor(kRed);
405 
406  h_c2_sum->SetLineWidth(2);
407  h_c2_h->SetLineWidth(2);
408  h_c2_e->SetLineWidth(2);
409 
410  p = (TPad *) c1->cd(idx++);
411  c1->Update();
412  p->SetLogy();
413 
414  T->Draw("clus_5x5_prod.sum_E>>h_5x5sum_c2_sum(170,-1,16)");
415  h_5x5sum_c2_sum->SetTitle(
416  "Cluster spectrum decomposition;5x5 cluster energy (GeV)");
417  T->Draw("clus_5x5_prod.sum_E>>h_5x5sum_c2_h(170,-1,16)", c2_h, "same");
418  T->Draw("clus_5x5_prod.sum_E>>h_5x5sum_c2_rej_h(170,-1,16)", !c2_h,
419  "same");
420  T->Draw("clus_5x5_prod.sum_E>>h_5x5sum_c2_e(170,-1,16)", c2_e, "same");
421 
422  h_5x5sum_c2_h->SetLineColor(kBlue);
423  h_5x5sum_c2_rej_h->SetLineColor(kMagenta);
424  h_5x5sum_c2_e->SetLineColor(kRed);
425 
426  h_5x5sum_c2_sum->SetLineWidth(2);
427  h_5x5sum_c2_h->SetLineWidth(2);
428  h_5x5sum_c2_rej_h->SetLineWidth(2);
429  h_5x5sum_c2_e->SetLineWidth(2);
430 
431  p = (TPad *) c1->cd(idx++);
432  c1->Update();
433  p->SetLogy();
434 
435  TH1 *h_5x5sum_c2_h2 = h_5x5sum_c2_h->DrawClone();
436  h_5x5sum_c2_h2->Sumw2();
437  h_5x5sum_c2_h2->SetMarkerColor(kBlue);
438  h_5x5sum_c2_h2->SetMarkerStyle(kFullCircle);
439  h_5x5sum_c2_h2->SetTitle(";5x5 cluster energy (GeV)");
440 
441  TH1 *h_5x5sum_c2_e2 = h_5x5sum_c2_e->DrawClone("same");
442  h_5x5sum_c2_e2->Sumw2();
443  h_5x5sum_c2_e2->SetMarkerColor(kRed);
444  h_5x5sum_c2_e2->SetMarkerStyle(kFullCircle);
445  h_5x5sum_c2_e2->SetTitle(";5x5 cluster energy (GeV)");
446 
447  p = (TPad *) c1->cd(idx++);
448  c1->Update();
449  p->SetLogy();
450 
451  TH1F *h_5x5sum_c2_h3 = h_5x5sum_c2_h->DrawClone();
452  h_5x5sum_c2_h3->SetName("h_5x5sum_c2_h3");
453  h_5x5sum_c2_h3->Sumw2();
454  h_5x5sum_c2_h3->Scale(1. / h_5x5sum_c2_h3->GetSum());
455  h_5x5sum_c2_h3->SetMarkerColor(kBlue);
456  h_5x5sum_c2_h3->SetMarkerStyle(kFullCircle);
457  h_5x5sum_c2_h3->SetTitle(";5x5 cluster energy (GeV);Probability / bin");
458 
459  SaveCanvas(c1,
460  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_" + TString(c1->GetName()), kTRUE);
461 }
462 
463 void LineShapeSim()
464 {
465  vector<double> mom;
466 
467  TText *t;
468  TCanvas *c1 = new TCanvas("LineShapeSim" + cuts, "LineShapeSim" + cuts, 1000,
469  650);
470 
471  // c1->Divide(2,2);
472  int idx = 1;
473  TPad *p;
474 
475  p = (TPad *) c1->cd(idx++);
476  c1->Update();
477  p->SetLogy();
478 
479  T->Draw("clus_5x5_prod.sum_E*SimEnergyScale>>h_5x5sum_c2_sum(170,-1,16)");
480  h_5x5sum_c2_sum->SetTitle(
481  "Cluster spectrum decomposition;5x5 cluster energy (GeV)");
482  h_5x5sum_c2_sum->SetLineWidth(2);
483 
484  SaveCanvas(c1,
485  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_" + TString(c1->GetName()), kTRUE);
486 }
487 
488 void HodoscopeCheck()
489 {
490  vector<double> mom;
491 
492  TText *t;
493  TCanvas *c1 = new TCanvas("HodoscopeCheck" + cuts, "HodoscopeCheck" + cuts,
494  1300, 950);
495 
496  c1->Divide(2, 1);
497  int idx = 1;
498  TPad *p;
499 
500  p = (TPad *) c1->cd(idx++);
501  c1->Update();
502  p->SetGridx(1);
503  p->SetGridy(1);
504  p->SetLogz();
505 
506  T->Draw("clus_5x5_prod.average_col:hodo_h>>h2_h(8,-.5,7.5,160,-.5,7.5)",
507  "good_e", "colz");
508  h2_h->SetTitle(
509  "Horizontal hodoscope check;Horizontal Hodoscope;5x5 cluster mean col");
510 
511  p = (TPad *) c1->cd(idx++);
512  c1->Update();
513  p->SetGridx(1);
514  p->SetGridy(1);
515  p->SetLogz();
516 
517  T->Draw("clus_5x5_prod.average_row:hodo_v>>h2_v(8,-.5,7.5,160,-.5,7.5)",
518  "good_e", "colz");
519  h2_v->SetTitle(
520  "Vertical hodoscope check;Vertical Hodoscope;5x5 cluster mean row");
521 
522  SaveCanvas(c1,
523  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_" + TString(c1->GetName()), kTRUE);
524 
525  return mom;
526 }
527 
528 void SimPositionCheck(const double shift_z = 0)
529 {
530  vector<double> mom;
531 
532  TText *t;
533  TCanvas *c1 = new TCanvas("SimPositionCheck" + cuts,
534  "SimPositionCheck" + cuts, 1300, 950);
535 
536  c1->Divide(2, 1);
537  int idx = 1;
538  TPad *p;
539 
540  p = (TPad *) c1->cd(idx++);
541  c1->Update();
542  p->SetGridx(1);
543  p->SetGridy(1);
544  p->SetLogz();
545 
546  T->Draw(
547  Form("clus_5x5_prod.average_col:truth_z>>h2_h(30,%f,%f,160,-.5,7.5)",
548  shift_z - 1.5, shift_z + 1.5),
549  "1", "colz");
550  h2_h->SetTitle(
551  "Horizontal hodoscope check;Horizontal beam pos;5x5 cluster mean col");
552 
553  p = (TPad *) c1->cd(idx++);
554  c1->Update();
555  p->SetGridx(1);
556  p->SetGridy(1);
557  p->SetLogz();
558 
559  T->Draw("clus_5x5_prod.average_row:truth_y>>h2_v(30,-1.5,1.5,160,-.5,7.5)",
560  "1", "colz");
561  h2_v->SetTitle(
562  "Vertical hodoscope check;Vertical beam pos;5x5 cluster mean row");
563 
564  SaveCanvas(c1,
565  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_" + TString(c1->GetName()), kTRUE);
566 
567  return;
568 }
569 
570 void Get_Res_linear_Summmary(TString e_sum = "sum_E", const double max_E = 32)
571 {
572  vector<double> beam_mom(GetBeamMom());
573 
574  // return;
575 
576  lin_res ges_clus_5x5_prod = GetResolution("clus_5x5_prod", beam_mom, kMagenta + 2, e_sum);
577  lin_res ges_clus_3x3_prod = GetResolution("clus_3x3_prod", beam_mom,
578  kBlue + 3, e_sum);
579  lin_res ges_clus_1x1_prod = GetResolution("clus_1x1_prod", beam_mom,
580  kCyan + 3, e_sum);
581  lin_res ges_clus_5x5_recalib = GetResolution("clus_5x5_recalib", beam_mom,
582  kRed + 3, e_sum);
583 
584  TCanvas *c1 = new TCanvas(Form("Res_linear") + cuts,
585  Form("Res_linear") + cuts, 1300, 600);
586  c1->Divide(2, 1);
587  int idx = 1;
588  TPad *p;
589 
590  p = (TPad *) c1->cd(idx++);
591  c1->Update();
592  // p->SetLogy();
593 
594  TLegend *leg = new TLegend(.15, .7, .6, .85);
595 
596  p->DrawFrame(0, 0, max_E, max_E,
597  Form("Electron Linearity;Input energy (GeV);Measured Energy (GeV)"));
598 
599  TF1 *f_calo_l_sim = new TF1("f_calo_l", "pol2", 0.5, max_E);
600  // f_calo_l_sim->SetParameters(-0.03389, 0.9666, -0.0002822);
601  f_calo_l_sim->SetParameters(-0., 1, -0.);
602  // f_calo_l_sim->SetLineWidth(1);
603  f_calo_l_sim->SetLineColor(kGreen + 2);
604  f_calo_l_sim->SetLineWidth(3);
605 
606  f_calo_l_sim->Draw("same");
607  ges_clus_5x5_prod.linearity->Draw("p");
608  ges_clus_3x3_prod.linearity->Draw("p");
609  ges_clus_1x1_prod.linearity->Draw("p");
610  ges_clus_5x5_recalib.linearity->Draw("p");
611  // ge_linear->Fit(f_calo_l, "RM0");
612  // f_calo_l->Draw("same");
613 
614  leg->AddEntry(ges_clus_5x5_prod.linearity, ges_clus_5x5_prod.name, "ep");
615  leg->AddEntry(ges_clus_3x3_prod.linearity, ges_clus_3x3_prod.name, "ep");
616  leg->AddEntry(ges_clus_1x1_prod.linearity, ges_clus_1x1_prod.name, "ep");
617  leg->AddEntry(ges_clus_5x5_recalib.linearity, "clus_5x5_recalib", "ep");
618  leg->AddEntry(f_calo_l_sim, "Unity", "l");
619  leg->Draw();
620 
621  p = (TPad *) c1->cd(idx++);
622  c1->Update();
623  // p->SetLogy();
624  p->SetGridx(0);
625  p->SetGridy(0);
626 
627  TF1 *f_calo_sim = new TF1("f_calo_sim", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5,
628  30);
629  // TF1 *f_calo_sim = new TF1("f_calo_sim", "sqrt([0]*[0]+[1]*[1]/x+[2]*[2]/x/x)/100", 0.5,
630  // 30);
631  f_calo_sim->SetParameters(3.7, 12.8, 0);
632  f_calo_sim->SetLineWidth(3);
633  f_calo_sim->SetLineColor(kGreen + 2);
634 
635  TH1 *hframe = p->DrawFrame(0, 0, max_E, 0.3,
636  Form("Resolution;Input energy (GeV);#DeltaE/<E>"));
637 
638  TLegend *leg = new TLegend(.2, .6, .85, .9);
639 
640  ges_clus_5x5_prod.f_res->Draw("same");
641  ges_clus_5x5_prod.resolution->Draw("ep");
642  // ges_clus_3x3_prod.f_res->Draw("same");
643  // ges_clus_3x3_prod.resolution->Draw("ep");
644  // ges_clus_1x1_prod.f_res->Draw("same");
645  // ges_clus_1x1_prod.resolution->Draw("ep");
646  ges_clus_5x5_recalib.f_res->Draw("same");
647  ges_clus_5x5_recalib.resolution->Draw("ep");
648  f_calo_sim->Draw("same");
649 
650  leg->AddEntry(ges_clus_5x5_prod.resolution, ges_clus_5x5_prod.name, "ep");
651  leg->AddEntry(ges_clus_5x5_prod.f_res,
652  Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
653  // Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E} ",
654  ges_clus_5x5_prod.f_res->GetParameter(0),
655  ges_clus_5x5_prod.f_res->GetParameter(1)
656  // ges_clus_5x5_prod.f_res->GetParameter(2)
657  ),
658  "l");
659 
660  leg->AddEntry(ges_clus_3x3_prod.resolution, ges_clus_3x3_prod.name, "ep");
661  leg->AddEntry(ges_clus_3x3_prod.f_res,
662  Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
663  // Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E} #oplus %.1f%%/E",
664  ges_clus_3x3_prod.f_res->GetParameter(0),
665  ges_clus_3x3_prod.f_res->GetParameter(1)
666  // ges_clus_3x3_prod.f_res->GetParameter(2)
667  ),
668  "l");
669  //
670  // leg->AddEntry(ges_clus_1x1_prod.resolution, ges_clus_1x1_prod.name, "ep");
671  // leg->AddEntry(ges_clus_1x1_prod.f_res,
672  // Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E} #oplus %.1f%%/E",
673  // ges_clus_1x1_prod.f_res->GetParameter(0),
674  // ges_clus_1x1_prod.f_res->GetParameter(1),
675  // ges_clus_1x1_prod.f_res->GetParameter(2)),
676  // "l");
677  //
678  // leg->AddEntry(ges_clus_5x5_recalib.resolution, "clus_5x5_recalib", "ep");
679  // leg->AddEntry(ges_clus_5x5_recalib.f_res,
680  // Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
681  // // Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E} #oplus %.1f%%/E",
682  // ges_clus_5x5_recalib.f_res->GetParameter(0),
683  // ges_clus_5x5_recalib.f_res->GetParameter(1)
684  // // ges_clus_5x5_recalib.f_res->GetParameter(2)
685  // ),
686  // "l");
687  // leg->AddEntry(new TH1(), "", "l");
688  // leg->AddEntry((TObject *) 0, " ", "");
689 
690  leg->Draw();
691 
692  TLegend *leg = new TLegend(.1, .15, .85, .25);
693 
694  leg->AddEntry(f_calo_sim,
695  Form(
696  "#splitline{Simulation w/ flat light collection}{#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}}",
697  f_calo_sim->GetParameter(0), f_calo_sim->GetParameter(1)),
698  "l");
699  leg->Draw();
700 
701  hframe->SetTitle("Electron Resolution");
702 
703  SaveCanvas(c1,
704  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_" + TString(c1->GetName()), kTRUE);
705 }
706 
707 void Get_Res_linear_Summmary_Sim()
708 {
709  vector<double> beam_mom(GetBeamMom());
710 
711  // return;
712 
713  lin_res ges_clus_5x5_prod = GetResolution("clus_5x5_prod", beam_mom, kBlue, "sum_E");
714  lin_res ges_clus_3x3_prod = GetResolution("clus_3x3_prod", beam_mom,
715  kBlue + 3, "sum_E");
716 
717  TCanvas *c1 = new TCanvas(Form("Res_linear") + cuts,
718  Form("Res_linear") + cuts, 1300, 600);
719  c1->Divide(2, 1);
720  int idx = 1;
721  TPad *p;
722 
723  p = (TPad *) c1->cd(idx++);
724  c1->Update();
725  // p->SetLogy();
726 
727  TLegend *leg = new TLegend(.15, .7, .6, .85);
728 
729  p->DrawFrame(0, 0, 25, 25,
730  Form("Electron Linearity;Input energy (GeV);Measured Energy (GeV)"));
731  TLine *l = new TLine(0, 0, 25, 25);
732  l->SetLineColor(kGray);
733  // l->Draw();
734 
735  TF1 *f_calo_l_sim = new TF1("f_calo_l", "pol2", 0.5, 25);
736  // f_calo_l_sim->SetParameters(-0.03389, 0.9666, -0.0002822);
737  f_calo_l_sim->SetParameters(-0., 1, -0.);
738  // f_calo_l_sim->SetLineWidth(1);
739  f_calo_l_sim->SetLineColor(kGreen + 2);
740  f_calo_l_sim->SetLineWidth(3);
741 
742  f_calo_l_sim->Draw("same");
743  ges_clus_5x5_prod.linearity->Draw("p");
744  ges_clus_3x3_prod.linearity->Draw("p");
745  // ge_linear->Fit(f_calo_l, "RM0");
746  // f_calo_l->Draw("same");
747 
748  leg->AddEntry(ges_clus_5x5_prod.linearity, ges_clus_5x5_prod.name, "ep");
749  leg->AddEntry(ges_clus_3x3_prod.linearity, ges_clus_3x3_prod.name, "ep");
750  leg->AddEntry(f_calo_l_sim, "Unity", "l");
751  leg->Draw();
752 
753  p = (TPad *) c1->cd(idx++);
754  c1->Update();
755  // p->SetLogy();
756  p->SetGridx(0);
757  p->SetGridy(0);
758 
759  TF1 *f_calo_sim = new TF1("f_calo_sim", "sqrt([0]*[0]+[1]*[1]/x)/100", 0.5,
760  25);
761  f_calo_sim->SetParameters(2.4, 11.8);
762  f_calo_sim->SetLineWidth(3);
763  f_calo_sim->SetLineColor(kGreen + 2);
764 
765  TH1 *hframe = p->DrawFrame(0, 0, 25, 0.3,
766  Form("Resolution;Input energy (GeV);#DeltaE/<E>"));
767 
768  TLegend *leg = new TLegend(.2, .6, .85, .9);
769 
770  ges_clus_5x5_prod.f_res->Draw("same");
771  ges_clus_5x5_prod.resolution->Draw("ep");
772  ges_clus_3x3_prod.f_res->Draw("same");
773  ges_clus_3x3_prod.resolution->Draw("ep");
774  f_calo_sim->Draw("same");
775 
776  leg->AddEntry(ges_clus_5x5_prod.resolution, ges_clus_5x5_prod.name, "ep");
777  leg->AddEntry(ges_clus_5x5_prod.f_res,
778  Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
779  ges_clus_5x5_prod.f_res->GetParameter(0),
780  ges_clus_5x5_prod.f_res->GetParameter(1)),
781  "l");
782 
783  leg->AddEntry(ges_clus_3x3_prod.resolution, ges_clus_3x3_prod.name, "ep");
784  leg->AddEntry(ges_clus_3x3_prod.f_res,
785  Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
786  ges_clus_3x3_prod.f_res->GetParameter(0),
787  ges_clus_3x3_prod.f_res->GetParameter(1)),
788  "l");
789 
790  // leg->AddEntry(new TH1(), "", "l");
791  // leg->AddEntry((TObject*) 0, " ", "");
792 
793  leg->Draw();
794 
795  TLegend *leg = new TLegend(.2, .1, .85, .3);
796 
797  leg->AddEntry(f_calo_sim,
798  Form("Prelim. Sim., #DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",
799  f_calo_sim->GetParameter(0), f_calo_sim->GetParameter(1)),
800  "l");
801  leg->Draw();
802 
803  hframe->SetTitle("Electron Resolution");
804 
805  SaveCanvas(c1,
806  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_" + TString(c1->GetName()), kTRUE);
807 }
808 
809 vector<double>
810 GetBeamMom()
811 {
812  vector<double> mom;
813 
814  TH1F *hbeam_mom = new TH1F("hbeam_mom", ";beam momentum (GeV)", 32, .5,
815  32.5);
816 
817  TText *t;
818  TCanvas *c1 = new TCanvas("GetBeamMom" + cuts, "GetBeamMom" + cuts, 1800,
819  900);
820 
821  T->Draw("abs(info.beam_mom)>>hbeam_mom");
822 
823  for (int bin = 1; bin < hbeam_mom->GetNbinsX(); bin++)
824  {
825  if (hbeam_mom->GetBinContent(bin) > 40)
826  {
827  const double momentum = hbeam_mom->GetBinCenter(bin);
828 
829  if (momentum == 1 || momentum == 2 || momentum == 3 || momentum == 4 || momentum == 5 || momentum == 6 || momentum == 8 || momentum == 12 || momentum == 16 || momentum == 24 || momentum == 28 || momentum == 32)
830  {
831  mom.push_back(momentum);
832 
833  cout << "GetBeamMom - " << momentum << " GeV for "
834  << hbeam_mom->GetBinContent(bin) << " event" << endl;
835  }
836  }
837  }
838 
839  SaveCanvas(c1,
840  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_" + TString(c1->GetName()), kTRUE);
841 
842  return mom;
843 }
844 
845 lin_res
846 GetResolution(TString cluster_name, vector<double> beam_mom, Color_t col, TString e_sum = "sum_E")
847 {
848  vector<double> mean;
849  vector<double> mean_err;
850  vector<double> res;
851  vector<double> res_err;
852 
853  TCanvas *c1 = new TCanvas("GetResolution_LineShape_" + cluster_name + cuts,
854  "GetResolution_LineShape_" + cluster_name + cuts, 1800, 900);
855 
856  c1->Divide(4, 2);
857  int idx = 1;
858  TPad *p;
859 
860  for (int i = 0; i < beam_mom.size(); ++i)
861  {
862  p = (TPad *) c1->cd(idx++);
863  c1->Update();
864 
865  const double momemtum = beam_mom[i];
866  const TString histname = Form("hLineShape%.0fGeV_", momemtum) + cluster_name;
867 
868  TH1F *h = new TH1F(histname, histname + ";Observed energy (GeV)",
869  (momemtum < 6 ? 25 : 40), 0.5, momemtum * 1.5);
870  T->Draw(cluster_name + "." + e_sum + ">>" + histname,
871  Form("abs(abs(info.beam_mom)-%f)/%f<.06", momemtum, momemtum));
872 
873  h->Sumw2();
874 
875  TF1 *fgaus_g = new TF1("fgaus_LG_g_" + cluster_name, "gaus",
876  h->GetMean() - 1 * h->GetRMS(), h->GetMean() + 4 * h->GetRMS());
877  fgaus_g->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
878  h->GetMean() + 2 * h->GetRMS());
879  h->Fit(fgaus_g, "MR0N");
880 
881  TF1 *fgaus = new TF1("fgaus_LG_" + cluster_name, "gaus",
882  fgaus_g->GetParameter(1) - 2 * fgaus_g->GetParameter(2),
883  fgaus_g->GetParameter(1) + 3 * fgaus_g->GetParameter(2));
884  fgaus->SetParameters(fgaus_g->GetParameter(0), fgaus_g->GetParameter(1),
885  fgaus_g->GetParameter(2));
886  h->Fit(fgaus, "MR");
887 
888  h->SetLineWidth(2);
889  h->SetMarkerStyle(kFullCircle);
890  fgaus->SetLineWidth(2);
891 
892  h->SetTitle(
893  Form("%.0f GeV/c: #DeltaE/<E> = %.1f%%", momemtum,
894  100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
895 
896  mean.push_back(fgaus->GetParameter(1));
897  mean_err.push_back(fgaus->GetParError(1));
898  res.push_back(fgaus->GetParameter(2) / fgaus->GetParameter(1));
899  res_err.push_back(fgaus->GetParError(2) / fgaus->GetParameter(1));
900  }
901 
902  SaveCanvas(c1,
903  TString(_file0->GetName()) + "_DrawPrototype3ShowerCalib_" + TString(c1->GetName()), kTRUE);
904 
905  TGraphErrors *ge_linear = new TGraphErrors(beam_mom.size(), &beam_mom[0],
906  &mean[0], 0, &mean_err[0]);
907 
908  TGraphErrors *ge_res = new TGraphErrors(beam_mom.size(), &beam_mom[0],
909  &res[0], 0, &res_err[0]);
910  ge_res->GetHistogram()->SetStats(0);
911  ge_res->Print();
912 
913  lin_res ret;
914  ret.name = cluster_name;
915  ret.linearity = ge_linear;
916  ret.resolution = ge_res;
917  ret.f_res = new TF1("f_calo_r_" + cluster_name, "sqrt([0]*[0]+[1]*[1]/x)/100",
918  // ret.f_res = new TF1("f_calo_r_" + cluster_name, "sqrt([0]*[0]+[1]*[1]/x+[2]*[2]/x/x)/100",
919  0.5, 30);
920  ret.f_res->SetParLimits(0, 2, 20);
921  ret.f_res->SetParLimits(1, 10, 40);
922  ret.f_res->SetParLimits(2, 10, 10);
923  ret.f_res->SetParameters(2, 12, 10);
924  ge_res->Fit(ret.f_res, "RM0QN");
925 
926  static int MarkerStyle = kOpenCircle - 1;
927  ++MarkerStyle;
928 
929  ge_linear->SetLineColor(col);
930  ge_linear->SetMarkerColor(col);
931  ge_linear->SetLineWidth(2);
932  ge_linear->SetMarkerStyle(MarkerStyle);
933  ge_linear->SetMarkerSize(2);
934 
935  ge_res->SetLineColor(col);
936  ge_res->SetMarkerColor(col);
937  ge_res->SetLineWidth(2);
938  ge_res->SetMarkerStyle(MarkerStyle);
939  ge_res->SetMarkerSize(2);
940 
941  ge_res->GetHistogram()->SetStats(0);
942 
943  ret.f_res->SetLineColor(col);
944  ret.f_res->SetLineWidth(3);
945 
946  return ret;
947 }
948 
949 TGraphErrors *
950 FitProfile(const TH2 *h2)
951 {
952  TProfile *p2 = h2->ProfileX();
953 
954  int n = 0;
955  double x[1000];
956  double ex[1000];
957  double y[1000];
958  double ey[1000];
959 
960  for (int i = 1; i <= h2->GetNbinsX(); i++)
961  {
962  TH1D *h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
963 
964  if (h1->GetSum() < 30)
965  {
966  cout << "FitProfile - ignore bin " << i << endl;
967  continue;
968  }
969  else
970  {
971  cout << "FitProfile - fit bin " << i << endl;
972  }
973 
974  TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
975  p2->GetBinError(i) * 4);
976 
977  TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
978  -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
979 
980  fgaus.SetParameter(1, p2->GetBinContent(i));
981  fgaus.SetParameter(2, p2->GetBinError(i));
982 
983  h1->Fit(&fgaus, "MQ");
984 
985  f2.SetParameters(fgaus.GetParameter(0) / 2, fgaus.GetParameter(1),
986  fgaus.GetParameter(2), fgaus.GetParameter(0) / 2,
987  fgaus.GetParameter(2) / 4, 0);
988 
989  h1->Fit(&f2, "MQ");
990 
991  // new TCanvas;
992  // h1->Draw();
993  // fgaus.Draw("same");
994  // break;
995 
996  x[n] = p2->GetBinCenter(i);
997  ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
998  y[n] = fgaus.GetParameter(1);
999  ey[n] = fgaus.GetParError(1);
1000 
1001  // p2->SetBinContent(i, fgaus.GetParameter(1));
1002  // p2->SetBinError(i, fgaus.GetParameter(2));
1003 
1004  n++;
1005  delete h1;
1006  }
1007 
1008  return new TGraphErrors(n, x, y, ex, ey);
1009 }