Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DrawEMCalTower.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DrawEMCalTower.C
1 // $Id: $
2 
11 #include <cmath>
12 #include <TFile.h>
13 #include <TString.h>
14 #include <TLine.h>
15 #include <TTree.h>
16 #include <TLatex.h>
17 #include <TGraphErrors.h>
18 #include <cassert>
19 #include "SaveCanvas.C"
20 #include "SetOKStyle.C"
21 using namespace std;
22 
23 TFile * _file0 = NULL;
24 TTree * T = NULL;
25 TString cuts = "";
26 
27 void
28 DrawEMCalTower( //
29  const TString infile = // "data/Prototype2_CEMC.root_DSTReader.root" //
30 // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/15Degree_1Col_LightCollection//Prototype_e-_16_SegALL_DSTReader.root"//
31  "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./0Degree_1Col_LightCollectionSeanStoll_CrackScan/Prototype_e-_8_SegALL_DSTReader.root"//
32 // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/EMCal_sim/./10DegreeRot_1Col_LightCollectionSeanStoll_CrackScan/Prototype_e-_8_SegALL_DSTReader.root"//
33  )
34 {
35  SetOKStyle();
36  gStyle->SetOptStat(0);
37  gStyle->SetOptFit(1111);
38  TVirtualFitter::SetDefaultFitter("Minuit2");
39  gSystem->Load("libg4eval.so");
40  gSystem->Load("libqa_modules.so");
41 
42  if (!_file0)
43  {
44  TString chian_str = infile;
45  chian_str.ReplaceAll("ALL", "*");
46 
47  TChain * t = new TChain("T");
48  const int n = t->Add(chian_str);
49 
50  cout << "Loaded " << n << " root files with " << chian_str << endl;
51  assert(n > 0);
52 
53  T = t;
54 
55  _file0 = new TFile;
56  _file0->SetName(infile);
57  }
58 
59  assert(_file0);
60 
61  T->SetAlias("ActiveTower_LG",
62  "TOWER_CALIB_LG_CEMC[].get_binphi()<8 && TOWER_CALIB_LG_CEMC[].get_bineta()<8");
63  T->SetAlias("EnergySum_LG",
64  "1*Sum$(TOWER_CALIB_LG_CEMC[].get_energy() * ActiveTower_LG)");
65 
66  T->SetAlias("ActiveTower_HG",
67  "TOWER_CALIB_HG_CEMC[].get_binphi()<8 && TOWER_CALIB_HG_CEMC[].get_bineta()<8");
68  T->SetAlias("EnergySum_HG",
69  "1*Sum$(TOWER_CALIB_HG_CEMC[].get_energy() * ActiveTower_HG)");
70 
71 //
72  const TCut event_sel = "1*1";
73  cuts = "_all_event";
74 // const TCut event_sel = "Entry$==3";
75 // cuts = "ev3";
76 // const TCut event_sel = "Entry$<2";
77 // cuts = "2ev";
78 // const TCut event_sel = "Entry$<1000";
79 // cuts = "1000ev";
80 
81  cout << "Build event selection of " << (const char *) event_sel << endl;
82 
83  T->Draw(">>EventList", event_sel);
84  TEventList * elist = gDirectory->GetObjectChecked("EventList", "TEventList");
85  cout << elist->GetN() << " / " << T->GetEntriesFast() << " events selected"
86  << endl;
87 
88  T->SetEventList(elist);
89 
90 // EMCDistribution_Fast("HG");
91 // EMCDistribution_Fast("LG", true);
92 // EMCDistributionVSBeam_SUM();
93  EMCDistributionVSBeam_SUM("TOWER_CEMC", -0, 5); // 0 degree tilted
94 // EMCDistributionVSBeam_SUM("TOWER_CEMC", -15); // 10 degree tilted
95 // EMCDistributionVSBeam_SUM("TOWER_CEMC", -15, 5); // 10 degree tilted
96 // EMCDistribution_SUM();
97 }
98 
99 void
100 EMCDistribution_Fast(TString gain = "LG", bool log_scale = false)
101 {
102  TString hname = "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "")
103  + cuts;
104 
105  TH2 * h2 = NULL;
106  if (log_scale)
107  {
108  h2 = new TH2F(hname,
109  Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 300, 10e-3,
110  100, 64, -.5, 63.5);
111  QAHistManagerDef::useLogBins(h2->GetXaxis());
112  }
113  else
114  {
115  h2 = new TH2F(hname,
116  Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 100, -.050,
117  .5, 64, -.5, 63.5);
118  }
119 
120  T->Draw(
121  "TOWER_CALIB_" + gain + "_CEMC[].get_bineta() + 8* TOWER_CALIB_" + gain
122  + "_CEMC[].get_binphi():TOWER_CALIB_" + gain
123  + "_CEMC[].get_energy()>>" + hname, "", "goff");
124  TText * t;
125  TCanvas *c1 = new TCanvas(
126  "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "") + cuts,
127  "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "") + cuts, 1800,
128  950);
129  c1->Divide(8, 8, 0., 0.01);
130  int idx = 1;
131  TPad * p;
132 
133  for (int iphi = 8 - 1; iphi >= 0; iphi--)
134  {
135  for (int ieta = 0; ieta < 8; ieta++)
136  {
137  p = (TPad *) c1->cd(idx++);
138  c1->Update();
139 
140  if (log_scale)
141  {
142  p->SetLogy();
143  p->SetLogx();
144  }
145  p->SetGridx(0);
146  p->SetGridy(0);
147 
148  TString hname = Form("hEnergy_ieta%d_iphi%d", ieta, iphi)
149  + TString(log_scale ? "_Log" : "");
150 
151  TH1 * h = h2->ProjectionX(hname, ieta + 8 * iphi + 1,
152  ieta + 8 * iphi + 1); // axis bin number is encoded as ieta+8*iphi+1
153 
154  h->SetLineWidth(0);
155  h->SetLineColor(kBlue + 3);
156  h->SetFillColor(kBlue + 3);
157  h->GetXaxis()->SetTitleSize(.09);
158  h->GetXaxis()->SetLabelSize(.08);
159  h->GetYaxis()->SetLabelSize(.08);
160 
161  h->Draw();
162 
163  TText *t = new TText(.9, .9, Form("Col%d Row%d", ieta, iphi));
164  t->SetTextAlign(33);
165  t->SetTextSize(.15);
166  t->SetNDC();
167  t->Draw();
168  }
169  }
170 
171  SaveCanvas(c1,
172  TString(_file0->GetName()) + TString("_DrawEMCalTower_")
173  + TString(c1->GetName()), kTRUE);
174 
175 }
176 
177 void
178 EMCDistribution(TString gain = "LG", bool log_scale = false)
179 {
180 
181  TText * t;
182  TCanvas *c1 = new TCanvas(
183  "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "") + cuts,
184  "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "") + cuts, 1800,
185  1000);
186  c1->Divide(8, 8, 0., 0.01);
187  int idx = 1;
188  TPad * p;
189 
190  for (int iphi = 8 - 1; iphi >= 0; iphi--)
191  {
192  for (int ieta = 0; ieta < 8; ieta++)
193  {
194  p = (TPad *) c1->cd(idx++);
195  c1->Update();
196 
197  if (log_scale)
198  {
199  p->SetLogy();
200  p->SetLogx();
201  }
202  p->SetGridx(0);
203  p->SetGridy(0);
204 
205  TString hname = Form("hEnergy_ieta%d_iphi%d", ieta, iphi)
206  + TString(log_scale ? "_Log" : "");
207 
208  TH1 * h = NULL;
209 
210  if (log_scale)
211  h = new TH1F(hname,
212  Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 300,
213  5e-3, 100);
214  else
215  h = new TH1F(hname,
216  Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 100,
217  -.050, .5);
218 
219  h->SetLineWidth(0);
220  h->SetLineColor(kBlue + 3);
221  h->SetFillColor(kBlue + 3);
222  h->GetXaxis()->SetTitleSize(.09);
223  h->GetXaxis()->SetLabelSize(.08);
224  h->GetYaxis()->SetLabelSize(.08);
225 
226  if (log_scale)
227  QAHistManagerDef::useLogBins(h->GetXaxis());
228 
229  T->Draw("TOWER_CALIB_" + gain + "_CEMC[].get_energy()>>" + hname,
230  Form(
231  "TOWER_CALIB_%s_CEMC[].get_bineta()==%d && TOWER_CALIB_%s_CEMC[].get_binphi()==%d",
232  gain.Data(), ieta, gain.Data(), iphi), "");
233 
234  TText *t = new TText(.9, .9, Form("Col%d Row%d", ieta, iphi));
235  t->SetTextAlign(33);
236  t->SetTextSize(.15);
237  t->SetNDC();
238  t->Draw();
239  }
240  }
241 
242  SaveCanvas(c1,
243  TString(_file0->GetName()) + TString("_DrawEMCalTower_")
244  + TString(c1->GetName()), kTRUE);
245 
246 }
247 void
248 EMCDistribution_SUM(TString sTOWER = "TOWER_CEMC")
249 {
250  TH1 * EnergySum_LG = new TH1F("EnergySum_LG",
251  ";Low-gain Tower Energy Sum (GeV);Count / bin", 1500, 0, 100);
252  TH1 * EnergySum_HG = new TH1F("EnergySum_HG",
253  ";High-gain Tower Energy Sum (GeV);Count / bin", 300, 0, 3);
254 
255  T->Draw("EnergySum_LG>>EnergySum_LG", "", "goff");
256  T->Draw("EnergySum_HG>>EnergySum_HG", "", "goff");
257 
258  TText * t;
259  TCanvas *c1 = new TCanvas("EMCDistribution_SUM" + cuts,
260  "EMCDistribution_SUM" + cuts, 1000, 960);
261  c1->Divide(2, 2);
262  int idx = 1;
263  TPad * p;
264 
265  p = (TPad *) c1->cd(idx++);
266  c1->Update();
267  p->SetLogy();
268  p->SetGridx(0);
269  p->SetGridy(0);
270 
271  TH1 * h = (TH1 *) EnergySum_LG->DrawClone();
272  h->SetLineWidth(2);
273  h->SetLineColor(kBlue + 3);
274 // h->Sumw2();
275  h->GetXaxis()->SetRangeUser(0, h->GetMean() + 5 * h->GetRMS());
276 
277  p = (TPad *) c1->cd(idx++);
278  c1->Update();
279 // p->SetLogy();
280  p->SetGridx(0);
281  p->SetGridy(0);
282 
283  TH1 * h = (TH1 *) EnergySum_LG->DrawClone();
284 
285  TF1 * fgaus = new TF1("fgaus_LG", "gaus", 0, 100);
286  fgaus->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
287  h->GetMean() + 2 * h->GetRMS());
288  h->Fit(fgaus, "M");
289 
290  h->Sumw2();
291  h->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
292  h->GetMean() + 4 * h->GetRMS());
293 
294  h->SetLineWidth(2);
295  h->SetMarkerStyle(kFullCircle);
296 
297  h->SetTitle(
298  Form("#DeltaE/<E> = %.1f%%",
299  100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
300 
301  p = (TPad *) c1->cd(idx++);
302  c1->Update();
303  p->SetLogy();
304  p->SetGridx(0);
305  p->SetGridy(0);
306 
307  TH1 * h = (TH1 *) EnergySum_HG->DrawClone();
308  h->SetLineWidth(2);
309  h->SetLineColor(kBlue + 3);
310 // h->Sumw2();
311  h->GetXaxis()->SetRangeUser(0, h->GetMean() + 5 * h->GetRMS());
312 
313  p = (TPad *) c1->cd(idx++);
314  c1->Update();
315 // p->SetLogy();
316  p->SetGridx(0);
317  p->SetGridy(0);
318 
319  TH1 * h = (TH1 *) EnergySum_HG->DrawClone();
320 
321  TF1 * fgaus = new TF1("fgaus_HG", "gaus", 0, 100);
322  fgaus->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
323  h->GetMean() + 2 * h->GetRMS());
324  h->Fit(fgaus, "M");
325 
326  h->Sumw2();
327  h->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
328  h->GetMean() + 4 * h->GetRMS());
329 
330  h->SetLineWidth(2);
331  h->SetMarkerStyle(kFullCircle);
332 
333  h->SetTitle(
334  Form("#DeltaE/<E> = %.1f%%",
335  100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
336 
337  SaveCanvas(c1,
338  TString(_file0->GetName()) + TString("_DrawEMCalTower_")
339  + TString(c1->GetName()), kTRUE);
340 
341 }
342 
343 void
344 EMCDistributionVSBeam_SUM(TString sTOWER = "TOWER_CEMC", const double z_shift =
345  0, const int n_div = 1)
346 {
347  TH3F * EnergySum_LG3 =
348  new TH3F("EnergySum_LG3",
349  ";Beam Horizontal Pos (cm);Beam Horizontal Vertical (cm);Low-gain Tower Energy Sum (GeV)", //
350  20 * n_div, z_shift - 5, z_shift + 5, //
351  20 * n_div, -5, 5, //
352  200, 0, 20);
353 
354  T->Draw("EnergySum_LG:PHG4VtxPoint.vy:PHG4VtxPoint.vz>>EnergySum_LG3", "",
355  "goff");
356 
357  TProfile2D * EnergySum_LG3_xy = EnergySum_LG3->Project3DProfile("yx");
358  TH2 * EnergySum_LG3_zx = EnergySum_LG3->Project3D("zx");
359  TH2 * EnergySum_LG3_zy = EnergySum_LG3->Project3D("zy");
360 
361  TGraphErrors * ge_EnergySum_LG3_zx = FitProfile(EnergySum_LG3_zx);
362  TGraphErrors * ge_EnergySum_LG3_zy = FitProfile(EnergySum_LG3_zy);
363 
364  TText * t;
365  TCanvas *c1 = new TCanvas(
366  TString(Form("EMCDistributionVSBeam_SUM_NDiv%d", n_div)) + cuts,
367  TString(Form("EMCDistributionVSBeam_SUM_NDiv%d", n_div)) + cuts, 1000,
368  960);
369  c1->Divide(2, 2);
370  int idx = 1;
371  TPad * p;
372 
373  p = (TPad *) c1->cd(idx++);
374  c1->Update();
375 // p->SetLogy();
376  p->SetGridx(0);
377  p->SetGridy(0);
378 
379  EnergySum_LG3_xy->Draw("colz");
380  EnergySum_LG3_xy->SetTitle(
381  "Position scan;Beam Horizontal Pos (cm);Beam Horizontal Vertical (cm)");
382 
383  p = (TPad *) c1->cd(idx++);
384  c1->Update();
385 // p->SetLogy();
386  p->SetGridx(0);
387  p->SetGridy(0);
388 
389  EnergySum_LG3_zx->Draw("colz");
390  EnergySum_LG3_zx->SetTitle(
391  "Position scan;Beam Horizontal Pos (cm);Energy Sum (GeV)");
392 
393  ge_EnergySum_LG3_zx->SetLineWidth(2);
394  ge_EnergySum_LG3_zx->SetMarkerStyle(kFullCircle);
395  ge_EnergySum_LG3_zx->Draw("pe");
396 
397  p = (TPad *) c1->cd(idx++);
398  c1->Update();
399 // p->SetLogy();
400  p->SetGridx(0);
401  p->SetGridy(0);
402 
403  EnergySum_LG3_zy->Draw("colz");
404  EnergySum_LG3_zy->SetTitle(
405  "Position scan;Beam Vertical Pos (cm);Energy Sum (GeV)");
406 
407  ge_EnergySum_LG3_zy->SetLineWidth(2);
408  ge_EnergySum_LG3_zy->SetMarkerStyle(kFullCircle);
409  ge_EnergySum_LG3_zy->Draw("pe");
410 
411  p = (TPad *) c1->cd(idx++);
412  c1->Update();
413 // p->SetLogy();
414  p->SetGridx(0);
415  p->SetGridy(0);
416 
417  TH1 * h = (TH1 *) EnergySum_LG3->ProjectionZ();
418 
419  TF1 * fgaus = new TF1("fgaus_LG", "gaus", 0, 100);
420  fgaus->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
421  h->GetMean() + 2 * h->GetRMS());
422  h->Fit(fgaus, "M");
423 
424  h->Sumw2();
425  h->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
426  h->GetMean() + 4 * h->GetRMS());
427  EnergySum_LG3_zx->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
428  h->GetMean() + 4 * h->GetRMS());
429  EnergySum_LG3_zy->GetYaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
430  h->GetMean() + 4 * h->GetRMS());
431 
432  h->SetLineWidth(2);
433  h->SetMarkerStyle(kFullCircle);
434 
435  h->SetTitle(
436  Form("#DeltaE/<E> = %.1f%%",
437  100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
438 
439  SaveCanvas(c1,
440  TString(_file0->GetName()) + TString("_DrawEMCalTower_")
441  + TString(c1->GetName()), kTRUE);
442 
443 }
444 
445 TGraphErrors *
446 FitProfile(const TH2 * h2)
447 {
448 
449  TProfile * p2 = h2->ProfileX();
450 
451  int n = 0;
452  double x[1000];
453  double ex[1000];
454  double y[1000];
455  double ey[1000];
456 
457  for (int i = 1; i <= h2->GetNbinsX(); i++)
458  {
459  TH1D * h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
460 
461  if (h1->GetSum() < 30)
462  {
463  cout << "FitProfile - ignore bin " << i << endl;
464  continue;
465  }
466  else
467  {
468  cout << "FitProfile - fit bin " << i << endl;
469  }
470 
471  TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
472  p2->GetBinError(i) * 4);
473 
474  TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
475  -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
476 
477  fgaus.SetParameter(1, p2->GetBinContent(i));
478  fgaus.SetParameter(2, p2->GetBinError(i));
479 
480  h1->Fit(&fgaus, "MQ");
481 
482  f2.SetParameters(fgaus.GetParameter(0) / 2, fgaus.GetParameter(1),
483  fgaus.GetParameter(2), fgaus.GetParameter(0) / 2,
484  fgaus.GetParameter(2) / 4, 0);
485 
486  h1->Fit(&f2, "MQ");
487 
488 // new TCanvas;
489 // h1->Draw();
490 // fgaus.Draw("same");
491 // break;
492 
493  x[n] = p2->GetBinCenter(i);
494  ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
495  y[n] = fgaus.GetParameter(1);
496  ey[n] = fgaus.GetParError(1);
497 
498 // p2->SetBinContent(i, fgaus.GetParameter(1));
499 // p2->SetBinError(i, fgaus.GetParameter(2));
500 
501  n++;
502  delete h1;
503  }
504 
505  return new TGraphErrors(n, x, y, ex, ey);
506 }
507