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