Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DrawPrototype2EMCalTower.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DrawPrototype2EMCalTower.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 //#include "Prototype2_DSTReader.h"
24 
25 TFile * _file0 = NULL;
26 TTree * T = NULL;
27 TString cuts = "";
28 double beam_momentum_selection = -16;
29 
30 void
31 DrawPrototype2EMCalTower( //
32  const TString infile = "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/Production_0715_EMCalSet2_HCalPR12/beam_00002237-0000_DSTReader.root", //
33  bool plot_all = false, const double momentum = -16)
34 {
35  beam_momentum_selection = momentum;
36 
37  SetOKStyle();
38  gStyle->SetOptStat(0);
39  gStyle->SetOptFit(1111);
40  TVirtualFitter::SetDefaultFitter("Minuit2");
41  gSystem->Load("libg4eval.so");
42  gSystem->Load("libqa_modules.so");
43  gSystem->Load("libPrototype2.so");
44 
45 // gROOT->LoadMacro("Prototype2_DSTReader.C+");
46 
47  if (!_file0)
48  {
49  TString chian_str = infile;
50  chian_str.ReplaceAll("ALL", "*");
51 
52  TChain * t = new TChain("T");
53  const int n = t->Add(chian_str);
54 
55  cout << "Loaded " << n << " root files with " << chian_str << endl;
56  assert(n > 0);
57 
58  T = t;
59 
60  _file0 = new TFile;
61  _file0->SetName(infile);
62  }
63 
64  assert(_file0);
65 
66  T->SetAlias("ActiveTower_LG",
67  "TOWER_LG_CEMC[].get_binphi()<8 && TOWER_LG_CEMC[].get_bineta()<8");
68  T->SetAlias("EnergySum_LG",
69  "1*Sum$(TOWER_LG_CEMC[].get_energy() * ActiveTower_LG)");
70 
71  T->SetAlias("ActiveTower_HG",
72  "TOWER_HG_CEMC[].get_binphi()<8 && TOWER_HG_CEMC[].get_bineta()<8");
73  T->SetAlias("EnergySum_HG",
74  "1*Sum$(TOWER_HG_CEMC[].get_energy() * ActiveTower_HG)");
75 
76 // T->SetAlias("C2_Inner_e", "1*TOWER_RAW_C2[0].energy");
77  T->SetAlias("C2_Inner_e", "1*abs(TOWER_RAW_C2[2].energy)");
78  T->SetAlias("C2_Outer_e", "1*abs(TOWER_RAW_C2[3].energy)");
79  T->SetAlias("C2_Sum_e", "C2_Inner_e + C2_Outer_e");
80 
81 // "TOWER_CALIB_CEMC.energy * ( Sum$( TOWER_CALIB_CEMC.get_column()==2 && TOWER_CALIB_CEMC.get_row()==1
82 
83  T->SetAlias("Average_column",
84  "Sum$(TOWER_CALIB_CEMC.get_column() * TOWER_CALIB_CEMC.get_energy())/Sum$(TOWER_CALIB_CEMC.get_energy())");
85  T->SetAlias("Average_row",
86  "Sum$(TOWER_CALIB_CEMC.get_row() * TOWER_CALIB_CEMC.get_energy())/Sum$(TOWER_CALIB_CEMC.get_energy())");
87 
88  T->SetAlias("Average_HODO_VERTICAL",
89  "Sum$(TOWER_CALIB_HODO_VERTICAL.towerid * (abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) * abs(TOWER_CALIB_HODO_VERTICAL.energy))/Sum$((abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) * abs(TOWER_CALIB_HODO_VERTICAL.energy))");
90  T->SetAlias("Valid_HODO_VERTICAL",
91  "Sum$(abs(TOWER_CALIB_HODO_VERTICAL.energy)>30) > 0");
92 
93  T->SetAlias("Average_HODO_HORIZONTAL",
94  "Sum$(TOWER_CALIB_HODO_HORIZONTAL.towerid * (abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) * abs(TOWER_CALIB_HODO_HORIZONTAL.energy))/Sum$((abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) * abs(TOWER_CALIB_HODO_HORIZONTAL.energy))");
95  T->SetAlias("Valid_HODO_HORIZONTAL",
96  "Sum$(abs(TOWER_CALIB_HODO_HORIZONTAL.energy)>30) > 0");
97 
98  T->SetAlias("No_Triger_VETO",
99  "Sum$(abs(TOWER_RAW_TRIGGER_VETO.energy)>15)==0");
100 
101  T->SetAlias("Energy_Sum_col1_row2_3x3",
102  "Sum$( (abs(TOWER_CALIB_CEMC.get_column()-1)<=1 && abs(TOWER_CALIB_CEMC.get_row()-2)<=1 ) * TOWER_CALIB_CEMC.get_energy())");
103  T->SetAlias("Energy_Sum_col1_row2_5x5",
104  "Sum$( (abs(TOWER_CALIB_CEMC.get_column()-1)<=2 && abs(TOWER_CALIB_CEMC.get_row()-2)<=2 ) * TOWER_CALIB_CEMC.get_energy())");
105  T->SetAlias("Energy_Sum_col2_row2_5x5",
106  "Sum$( (abs(TOWER_CALIB_CEMC.get_column()-2)<=2 && abs(TOWER_CALIB_CEMC.get_row()-2)<=2 ) * TOWER_CALIB_CEMC.get_energy())");
107  T->SetAlias("Energy_Sum_CEMC", "1*Sum$(TOWER_CALIB_CEMC.get_energy())");
108 
109  // 12 GeV calibration
110 // EDM=9.83335e-18 STRATEGY= 1 ERROR MATRIX ACCURATE
111 //EXT PARAMETER STEP FIRST
112 //NO. NAME VALUE ERROR SIZE DERIVATIVE
113 //1 p0 1.19768e+01 7.30605e-02 3.76799e-05 4.24290e-09
114 //2 p1 8.71776e+00 6.82987e-02 3.52240e-05 6.80808e-08
115 
116 // T->SetAlias("Energy_Sum_Hadron_CEMC",
117 // "1*Sum$(TOWER_CALIB_CEMC.get_energy())"); // full bias
118  T->SetAlias("Energy_Sum_Hadron_CEMC",
119  "1.14*12./8.71776e+00*Sum$(TOWER_CALIB_CEMC.get_energy())"); // full bias
120 // T->SetAlias("Energy_Sum_Hadron_CEMC",
121 // "1.14*12./8.71776e+00*(16./6.93250e+00)*(28/33.3405)*Sum$(TOWER_CALIB_CEMC.get_energy())"); // half-gain bias
122  T->SetAlias("CEMC_MIP", "Energy_Sum_Hadron_CEMC<0.7");
123 
124  // 12 GeV calibration
125 // FCN=9.63681 FROM HESSE STATUS=OK 14 CALLS 56 TOTAL
126 // EDM=1.49963e-17 STRATEGY= 1 ERROR MATRIX ACCURATE
127 // EXT PARAMETER STEP FIRST
128 // NO. NAME VALUE ERROR SIZE DERIVATIVE
129 // 1 p0 9.50430e+00 8.42691e-02 3.83666e-06 1.41105e-08
130 // 2 p1 6.99727e+00 1.06583e-01 4.85258e-06 5.85631e-08
131 
132 // T->SetAlias("Energy_Sum_Hadron_HCALIN",
133 // "1*Sum$(TOWER_CALIB_LG_HCALIN.get_energy())");
134 // T->SetAlias("HCALIN_MIP", "Energy_Sum_Hadron_HCALIN<0.5");
135 // T->SetAlias("Energy_Sum_Hadron_HCALOUT",
136 // "1*Sum$(TOWER_CALIB_LG_HCALOUT.get_energy())");
137  T->SetAlias("Energy_Sum_Hadron_HCALIN",
138  "12./6.99727e+00*Sum$(TOWER_CALIB_LG_HCALIN.get_energy())");
139  T->SetAlias("HCALIN_MIP", "Energy_Sum_Hadron_HCALIN<0.5");
140  T->SetAlias("Energy_Sum_Hadron_HCALOUT",
141  "12./9.50430e+00*Sum$(TOWER_CALIB_LG_HCALOUT.get_energy())");
142 
143  T->SetAlias("MIP_Count_Col2",
144  "Sum$( abs( TOWER_RAW_CEMC.get_energy() )>20 && abs( TOWER_RAW_CEMC.get_energy() )<200 && TOWER_CALIB_CEMC.get_column() == 2 )");
145  T->SetAlias("Pedestal_Count_AllCEMC",
146  "Sum$( abs( TOWER_RAW_CEMC.get_energy() )<20)");
147 
148 //
149  TCut event_sel = "1*1";
150 
151  if (plot_all)
152  {
153 // event_sel = "1*1";
154 // cuts = "_all_event";
155 // event_sel = "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL";
156 // cuts = "_Valid_HODO";
157  event_sel =
158  "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO";
159  cuts = "_Valid_HODO_Trigger_VETO";
160  }
161  else
162  {
163  if (0)
164  {// energy selection
165  event_sel = Form("(beam_MTNRG_GeV == %f)", beam_momentum_selection);
166  cuts = Form("_%.0fGeV", beam_momentum_selection);
167 
168  cout << "Build event selection of " << (const char *) event_sel << endl;
169 
170  T->Draw(">>EventListRunCut", event_sel);
171  TEventList * elist = gDirectory->GetObjectChecked("EventListRunCut",
172  "TEventList");
173  cout << elist->GetN() << " / " << T->GetEntriesFast()
174  << " events selected" << endl;
175  T->SetEventList(elist);
176  }
177 
178 // event_sel = "1*1";
179 // cuts = "_all_event";
180 // event_sel = "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL";
181 // cuts = "_Valid_HODO";
182  event_sel =
183  "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO";
184  cuts = "_Valid_HODO_Trigger_VETO";
185 
186 // event_sel =
187 // event_sel
188 // + "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO && (C2_Sum_e<50)";
189 // cuts = cuts + "_Valid_HODO_Trigger_VETO_C2_Sum_Hadron";
190 
191 // event_sel =
192 // event_sel
193 // + "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && No_Triger_VETO && CEMC_MIP && (C2_Sum_e<100)";
194 // cuts = cuts + "_Valid_HODO_Trigger_VETO_CEMC_MIP_C2_Sum_Hadron";
195 
196 // event_sel = "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && (MIP_Count_Col2 == 8) && (Pedestal_Count_AllCEMC >= 64 - 8)";
197 // cuts = "_Valid_HODO_MIP_Col2_PedestalOther";
198 
199 // event_sel = "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && Average_HODO_VERTICAL>1.5 && Average_HODO_VERTICAL<4.5 && Average_HODO_HORIZONTAL>3.5 && Average_HODO_HORIZONTAL<6.5";
200 // cuts = "_Valid_HODO_center_col1_row2";
201 
202 // event_sel =
203 // "C2_Inner_e>100 && Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && Average_HODO_HORIZONTAL>1.5 && Average_HODO_HORIZONTAL<4.5 && Average_HODO_VERTICAL>3.5 && Average_HODO_VERTICAL<6.5";
204 // cuts = "_Valid_HODO_center_col1_row2_C2";
205 
206 // event_sel = "abs(Average_column - 1)<.5 && abs(Average_row - 2) < .5";
207 // cuts = "_tower_center_col1_row2";
208 // event_sel = "Valid_HODO_HORIZONTAL && Valid_HODO_VERTICAL && abs(C2_Inner_e)<10";
209 // cuts = "_Valid_HODO_VetoC2";
210  }
211 
212  cout << "Build event selection of " << (const char *) event_sel << endl;
213 
214  T->Draw(">>EventList", event_sel);
215  TEventList * elist = gDirectory->GetObjectChecked("EventList", "TEventList");
216  cout << elist->GetN() << " / " << T->GetEntriesFast() << " events selected"
217  << endl;
218 
219  T->SetEventList(elist);
220 
221 // int rnd = rand();
222 // gDirectory->mkdir(Form("dir_%d", rnd));
223 // gDirectory->cd(Form("dir_%d", rnd));
224 // if (plot_all)
225 // EMCDistribution_SUM("Energy_Sum_col1_row2_5x5");
226 
227  int rnd = rand();
228  gDirectory->mkdir(Form("dir_%d", rnd));
229  gDirectory->cd(Form("dir_%d", rnd));
230  if (plot_all)
231  EMCDistribution_ShowShape("C2_Sum_e");
232 
233  int rnd = rand();
234  gDirectory->mkdir(Form("dir_%d", rnd));
235  gDirectory->cd(Form("dir_%d", rnd));
236  if (plot_all)
237  EMCDistribution_SUM("Energy_Sum_CEMC", "C2_Sum_e");
238 
239  int rnd = rand();
240  gDirectory->mkdir(Form("dir_%d", rnd));
241  gDirectory->cd(Form("dir_%d", rnd));
242  if (plot_all)
243  EMCDistribution_SUM("Energy_Sum_col2_row2_5x5", "C2_Sum_e");
244 
245  int rnd = rand();
246  gDirectory->mkdir(Form("dir_%d", rnd));
247  gDirectory->cd(Form("dir_%d", rnd));
248  if (plot_all)
249  EMCDistribution_HCalCalibration();
250 
251  int rnd = rand();
252  gDirectory->mkdir(Form("dir_%d", rnd));
253  gDirectory->cd(Form("dir_%d", rnd));
254  if (plot_all)
255  EMCDistribution_Fast();
256 
257  int rnd = rand();
258  gDirectory->mkdir(Form("dir_%d", rnd));
259  gDirectory->cd(Form("dir_%d", rnd));
260  if (plot_all)
261  EMCDistribution_Fast("RAW");
262 
263  int rnd = rand();
264  gDirectory->mkdir(Form("dir_%d", rnd));
265  gDirectory->cd(Form("dir_%d", rnd));
266  if (plot_all)
267  EMCDistribution_PeakSample_Fast();
268 
269  int rnd = rand();
270  gDirectory->mkdir(Form("dir_%d", rnd));
271  gDirectory->cd(Form("dir_%d", rnd));
272  if (plot_all)
273  EMCDistribution_Fast("CALIB", true);
274 
275  int rnd = rand();
276  gDirectory->mkdir(Form("dir_%d", rnd));
277  gDirectory->cd(Form("dir_%d", rnd));
278  if (plot_all)
279  EMCDistribution_ADC();
280 
281 // if (!plot_all)
282 // T->Process("Prototype2_DSTReader.C+",
283 // TString(_file0->GetName())
284 // + TString("_DrawPrototype2EMCalTower_Prototype2_DSTReader") + cuts
285 // + TString(".dat"));
286 }
287 
288 void
289 EMCDistribution_HCalCalibration()
290 {
291  gStyle->SetOptStat(0);
292  gStyle->SetOptFit(1111);
293 
294  TH1F * h1_HCalOut = new TH1F("h1_HCalOut",
295  "(CEMC MIP && HCal_{IN} MIP);EMCal + HCal_{IN} + HCal_{OUT} (GeV)", 100,
296  0, abs(beam_momentum_selection) * 2.);
297 
298  TH2 * h2_HCalOut_HCalIn =
299  new TH2F("h2_HCalOut_HCalIn",
300  "Energy Balance HCal_{OUT} and HCal_{IN} (CEMC MIP && NOT HCal_{IN} MIP);(HCal_{OUT} - HCal_{IN})/HCal Sum;EMCal + HCal_{IN} + HCal_{OUT} (GeV)",
301  10, -1, 1, 100, 0, abs(beam_momentum_selection) * 2.);
302 
303  TH2 * h2_HCal_EMCal =
304  new TH2F("h2_HCal_EMCal",
305  "Energy Balance HCal_{SUM} and EMCal (NOT CEMC MIP);(HCal_{SUM} - EMCal)/Sum;EMCal + HCal_{IN} + HCal_{OUT} (GeV)",
306  10, -1, 1, 100, 0, abs(beam_momentum_selection) * 2.);
307 
308  TText * t;
309  TCanvas *c1 = new TCanvas("EMCDistribution_HCalCalibration" + cuts,
310  "EMCDistribution_HCalCalibration" + cuts, 1800, 600);
311  c1->Divide(3, 1);
312  int idx = 1;
313  TPad * p;
314 
315  p = (TPad *) c1->cd(idx++);
316  c1->Update();
317 
318  T->Draw(
319  "Energy_Sum_Hadron_CEMC + Energy_Sum_Hadron_HCALIN + Energy_Sum_Hadron_HCALOUT:(Energy_Sum_Hadron_HCALOUT - Energy_Sum_Hadron_HCALIN)/(Energy_Sum_Hadron_HCALIN+Energy_Sum_Hadron_HCALOUT)>>h2_HCalOut_HCalIn",
320  "(!HCALIN_MIP) && CEMC_MIP", "goff");
321  h2_HCalOut_HCalIn->FitSlicesY();
322 
323  TH1 * h2_HCalOut_HCalIn_1 = (TH1 *) gDirectory->GetObjectChecked(
324  "h2_HCalOut_HCalIn_1", "TH1D");
325  assert(h2_HCalOut_HCalIn_1);
326  h2_HCalOut_HCalIn_1->SetMarkerStyle(kFullCircle);
327 
328  TF1 * f_HCalOut_HCalIn = new TF1("f_HCalOut_HCalIn",
329  "[0] * (x+1)/2 + [1]*(1-x)/2", -.8, 0.8);
330  f_HCalOut_HCalIn->SetLineColor(kMagenta);
331  f_HCalOut_HCalIn->SetLineWidth(3);
332  h2_HCalOut_HCalIn_1->Fit(f_HCalOut_HCalIn, "MR0");
333 
334  h2_HCalOut_HCalIn->DrawClone("colz");
335  h2_HCalOut_HCalIn_1->DrawClone("SAME");
336  f_HCalOut_HCalIn->DrawClone("same");
337 
338  p = (TPad *) c1->cd(idx++);
339  c1->Update();
340 
341  T->Draw(
342  "Energy_Sum_Hadron_CEMC + Energy_Sum_Hadron_HCALIN + Energy_Sum_Hadron_HCALOUT:(Energy_Sum_Hadron_HCALIN + Energy_Sum_Hadron_HCALOUT - Energy_Sum_Hadron_CEMC)/(Energy_Sum_Hadron_CEMC + Energy_Sum_Hadron_HCALIN + Energy_Sum_Hadron_HCALOUT)>>h2_HCal_EMCal",
343  "! CEMC_MIP", "goff");
344  h2_HCal_EMCal->FitSlicesY();
345 
346  TH1 * h2_HCal_EMCal_1 = (TH1 *) gDirectory->GetObjectChecked(
347  "h2_HCal_EMCal_1", "TH1D");
348  assert(h2_HCal_EMCal_1);
349  h2_HCal_EMCal_1->SetMarkerStyle(kFullCircle);
350 
351  TF1 * f_HCal_EMCal = new TF1("f_HCal_EMCal", "[0] * (x+1)/2 + [1]*(1-x)/2",
352  -.8, 0.8);
353  f_HCal_EMCal->SetLineColor(kMagenta);
354  f_HCal_EMCal->SetLineWidth(3);
355  h2_HCal_EMCal_1->Fit(f_HCal_EMCal, "MR0");
356 
357  h2_HCal_EMCal->DrawClone("colz");
358  h2_HCal_EMCal_1->DrawClone("SAME");
359  f_HCal_EMCal->DrawClone("same");
360 
361  p = (TPad *) c1->cd(idx++);
362  c1->Update();
363 
364  T->Draw(
365  "Energy_Sum_Hadron_CEMC + Energy_Sum_Hadron_HCALIN + Energy_Sum_Hadron_HCALOUT>>h1_HCalOut",
366  "HCALIN_MIP && CEMC_MIP", "goff");
367 
368  TH1 * h2_HCal_EMCal_ProjY = (TH1 *) h2_HCal_EMCal->ProjectionY(
369  "h2_HCal_EMCal_ProjY",2,9)->DrawClone();
370  TH1 * h2_HCalOut_HCalIn_ProjY = (TH1 *) h2_HCalOut_HCalIn->ProjectionY(
371  "h2_HCalOut_HCalIn_ProjY")->DrawClone("same");
372  h1_HCalOut->Draw("same");
373 
374  h2_HCal_EMCal_ProjY->SetTitle("Red: 3 Calo shower, Blue: MIP EMCal, Green: MIP EMCal & HCal_{IN}");
375 
376  h2_HCal_EMCal_ProjY->SetLineColor(kRed + 2);
377  h2_HCal_EMCal_ProjY->SetLineWidth(3);
378  h2_HCalOut_HCalIn_ProjY->SetLineColor(kBlue + 2);
379  h2_HCalOut_HCalIn_ProjY->SetLineWidth(2);
380  h1_HCalOut->SetLineColor(kGreen+3);
381  h1_HCalOut->SetLineWidth(3);
382 
383  SaveCanvas(c1,
384  TString(_file0->GetName()) + TString("_DrawPrototype2EMCalTower_")
385  + TString(c1->GetName()), false);
386 
387 }
388 
389 void
390 EMCDistribution_ShowShape(TString CherenkovSignal = "C2_Inner",
391  const double che_cut = 100)
392 {
393  TString cut_pass = CherenkovSignal + Form(">%.1f", che_cut);
394  TString cut_rej = CherenkovSignal + Form("<%.1f", che_cut);
395 
396  const double event_pass = T->GetEntries(cut_pass);
397  const double event_rej = T->GetEntries(cut_rej);
398 
399  TH2 * EnergyDist_pass = new TH2F("EnergyDist_pass",
400  cut_pass + ";Column;Row;<Energy> / Event / Tower", 8, -.5, 7.5, 8, -.5,
401  7.5);
402  TH2 * EnergyDist_rej = new TH2F("EnergyDist_rej",
403  cut_rej + ";Column;Row;<Energy> / Event / Tower", 8, -.5, 7.5, 8, -.5,
404  7.5);
405 
406  TH1 * Che_full = new TH1F("Che_full",
407  ";" + CherenkovSignal + " Signal (ADC);Count / bin", 200, 0, 2000);
408  TH1 * Che_pass = new TH1F("Che_pass",
409  ";" + CherenkovSignal + " Signal (ADC);Count / bin", 200, 0, 2000);
410  TH1 * Che_rej = new TH1F("Che_rej",
411  ";" + CherenkovSignal + " Signal (ADC);Count / bin", 200, 0, 2000);
412 
413  Che_full->SetLineColor(kBlue + 3);
414  Che_full->SetLineWidth(2);
415  Che_pass->SetLineColor(kGreen + 3);
416  Che_pass->SetLineWidth(2);
417  Che_pass->SetFillColor(kGreen + 3);
418  Che_pass->SetFillStyle(1);
419  Che_rej->SetLineColor(kBlue + 3);
420  Che_rej->SetLineWidth(2);
421  Che_rej->SetFillColor(kBlue + 3);
422  Che_rej->SetFillStyle(1);
423 
424  T->Draw(
425  "TOWER_CALIB_CEMC[].get_binphi():TOWER_CALIB_CEMC[].get_bineta()>>EnergyDist_pass",
426  Form("(%s) * (TOWER_CALIB_CEMC[].get_energy())", cut_pass.Data()),
427  "goff");
428  T->Draw(
429  "TOWER_CALIB_CEMC[].get_binphi():TOWER_CALIB_CEMC[].get_bineta()>>EnergyDist_rej",
430  Form("(%s) * (TOWER_CALIB_CEMC[].get_energy())", cut_rej.Data()), "goff");
431 
432  T->Draw(CherenkovSignal + ">>Che_full", NULL, "goff");
433  T->Draw(CherenkovSignal + ">>Che_pass", cut_pass, "goff");
434  T->Draw(CherenkovSignal + ">>Che_rej", cut_rej, "goff");
435 
436  EnergyDist_pass->Scale(1. / event_pass);
437  EnergyDist_rej->Scale(1. / event_rej);
438 
439  TText * t;
440  TCanvas *c1 = new TCanvas("EMCDistribution_ShowShape" + cuts,
441  "EMCDistribution_ShowShape" + cuts, 1800, 600);
442  c1->Divide(3, 1);
443  int idx = 1;
444  TPad * p;
445 
446  p = (TPad *) c1->cd(idx++);
447  c1->Update();
448  p->SetLogy();
449  p->SetGridx(0);
450  p->SetGridy(0);
451 
452  Che_full->DrawClone();
453  Che_pass->DrawClone("same");
454  Che_rej->DrawClone("same");
455 
456  p = (TPad *) c1->cd(idx++);
457  c1->Update();
458 // p->SetLogy();
459 
460  EnergyDist_pass->DrawClone("LEGO2Z");
461 
462  p = (TPad *) c1->cd(idx++);
463  c1->Update();
464 // p->SetLogy();
465 // p->SetGridx(0);
466 // p->SetGridy(0);
467 
468  EnergyDist_rej->DrawClone("LEGO2Z");
469 
470  SaveCanvas(c1,
471  TString(_file0->GetName()) + TString("_DrawPrototype2EMCalTower_")
472  + TString(c1->GetName()), false);
473 
474 }
475 
476 void
477 EMCDistribution_SUM(TString sTOWER = "Energy_Sum_col1_row2_5x5",
478  TString CherenkovSignal = "C2_Inner")
479 {
480  TH1 * EnergySum_LG_full = new TH1F("EnergySum_LG_full",
481  ";Full range Tower Energy Sum (GeV);Count / bin", 300, 0, 40);
482  TH1 * EnergySum_LG = new TH1F("EnergySum_LG",
483  ";Full range Tower Energy Sum (GeV);Count / bin", 300, 0, 40);
484 // TH1 * EnergySum_HG = new TH1F("EnergySum_HG",
485 // ";Low range Tower Energy Sum (ADC);Count / bin", 50, 0, 500);
486 
487  TH1 * C2_Inner_full = new TH1F("C2_Inner_full",
488  CherenkovSignal + ";Cherenkov Signal (ADC);Count / bin", 200, 0, 2000);
489  TH1 * C2_Inner = new TH1F("C2_Inner",
490  CherenkovSignal + ";Cherenkov Inner Signal (ADC);Count / bin", 200, 0,
491  2000);
492 
493  EnergySum_LG_full->SetLineColor(kBlue + 3);
494  EnergySum_LG_full->SetLineWidth(2);
495 
496  EnergySum_LG->SetLineColor(kGreen + 3);
497  EnergySum_LG->SetLineWidth(3);
498  EnergySum_LG->SetMarkerColor(kGreen + 3);
499 
500  C2_Inner_full->SetLineColor(kBlue + 3);
501  C2_Inner_full->SetLineWidth(2);
502 
503  C2_Inner->SetLineColor(kGreen + 3);
504  C2_Inner->SetLineWidth(3);
505  C2_Inner->SetMarkerColor(kGreen + 3);
506 
507  TCut c2 = CherenkovSignal + ">100";
508 
509  T->Draw(sTOWER + ">>EnergySum_LG_full", "", "goff");
510  T->Draw(sTOWER + ">>EnergySum_LG", c2, "goff");
511  T->Draw(CherenkovSignal + ">>C2_Inner_full", "", "goff");
512  T->Draw(CherenkovSignal + ">>C2_Inner", c2, "goff");
513 
514  TText * t;
515  TCanvas *c1 = new TCanvas(
516  "EMCDistribution_SUM_" + sTOWER + "_" + CherenkovSignal + cuts,
517  "EMCDistribution_SUM_" + sTOWER + "_" + CherenkovSignal + cuts, 1800,
518  600);
519  c1->Divide(3, 1);
520  int idx = 1;
521  TPad * p;
522 
523  p = (TPad *) c1->cd(idx++);
524  c1->Update();
525  p->SetLogy();
526  p->SetGridx(0);
527  p->SetGridy(0);
528 
529  C2_Inner_full->DrawClone();
530  C2_Inner->DrawClone("same");
531 
532  p = (TPad *) c1->cd(idx++);
533  c1->Update();
534  p->SetLogy();
535  p->SetGridx(0);
536  p->SetGridy(0);
537 
538  TH1 * h = (TH1 *) EnergySum_LG_full->DrawClone();
539  h->GetXaxis()->SetRangeUser(0, h->GetMean() + 5 * h->GetRMS());
540  (TH1 *) EnergySum_LG->DrawClone("same");
541 
542  p = (TPad *) c1->cd(idx++);
543  c1->Update();
544 // p->SetLogy();
545  p->SetGridx(0);
546  p->SetGridy(0);
547 
548  TH1 * h_full = (TH1 *) EnergySum_LG_full->DrawClone();
549  TH1 * h = (TH1 *) EnergySum_LG->DrawClone("same");
550 
551  TF1 * fgaus_g = new TF1("fgaus_LG_g", "gaus", h->GetMean() - 1 * h->GetRMS(),
552  h->GetMean() + 4 * h->GetRMS());
553  fgaus_g->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
554  h->GetMean() + 2 * h->GetRMS());
555  h->Fit(fgaus_g, "MR0N");
556 
557  TF1 * fgaus = new TF1("fgaus_LG", "gaus",
558  fgaus_g->GetParameter(1) - 1 * fgaus_g->GetParameter(2),
559  fgaus_g->GetParameter(1) + 4 * fgaus_g->GetParameter(2));
560  fgaus->SetParameters(fgaus_g->GetParameter(0), fgaus_g->GetParameter(1),
561  fgaus_g->GetParameter(2));
562  h->Fit(fgaus, "MR");
563 
564  h->Sumw2();
565  h_full->Sumw2();
566  h_full->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
567  h->GetMean() + 4 * h->GetRMS());
568 
569  h->SetLineWidth(2);
570  h->SetMarkerStyle(kFullCircle);
571 
572  h_full->SetTitle(
573  Form("#DeltaE/<E> = %.1f%%",
574  100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
575 
576 // p = (TPad *) c1->cd(idx++);
577 // c1->Update();
578 // p->SetLogy();
579 // p->SetGridx(0);
580 // p->SetGridy(0);
581 //
582 // TH1 * h = (TH1 *) EnergySum_LG->DrawClone();
583 // h->GetXaxis()->SetRangeUser(0,500);
584 // h->SetLineWidth(2);
585 // h->SetLineColor(kBlue + 3);
587 // h->GetXaxis()->SetRangeUser(0, h->GetMean() + 5 * h->GetRMS());
588 //
589 // p = (TPad *) c1->cd(idx++);
590 // c1->Update();
592 // p->SetGridx(0);
593 // p->SetGridy(0);
594 //
595 // TH1 * h = (TH1 *) EnergySum_LG->DrawClone();
596 // h->GetXaxis()->SetRangeUser(0,500);
597 //
598 // TF1 * fgaus = new TF1("fgaus_HG", "gaus", 0, 100);
599 // fgaus->SetParameters(1, h->GetMean() - 2 * h->GetRMS(),
600 // h->GetMean() + 2 * h->GetRMS());
601 // h->Fit(fgaus, "M");
602 //
603 // h->Sumw2();
604 // h->GetXaxis()->SetRangeUser(h->GetMean() - 4 * h->GetRMS(),
605 // h->GetMean() + 4 * h->GetRMS());
606 //
607 // h->SetLineWidth(2);
608 // h->SetMarkerStyle(kFullCircle);
609 //
610 // h->SetTitle(
611 // Form("#DeltaE/<E> = %.1f%%",
612 // 100 * fgaus->GetParameter(2) / fgaus->GetParameter(1)));
613 
614  SaveCanvas(c1,
615  TString(_file0->GetName()) + TString("_DrawPrototype2EMCalTower_")
616  + TString(c1->GetName()), false);
617 
618 }
619 
620 void
621 EMCDistribution_Fast(TString gain = "CALIB", bool full_gain = false)
622 {
623  TString hname = "EMCDistribution_" + gain
624  + TString(full_gain ? "_FullGain" : "") + cuts;
625 
626  TH2 * h2 = NULL;
627  if (gain.Contains("CALIB"))
628  {
629  if (full_gain)
630  {
631  h2 = new TH2F(hname,
632  Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 100, .05,
633  25, 64, -.5, 63.5);
634  QAHistManagerDef::useLogBins(h2->GetXaxis());
635  }
636  else
637  {
638  h2 = new TH2F(hname,
639  Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 260, -.2,
640  5, 64, -.5, 63.5);
641  }
642  T->Draw(
643  "TOWER_" + gain + "_CEMC[].get_bineta() + 8* TOWER_" + gain
644  + "_CEMC[].get_binphi():TOWER_" + gain + "_CEMC[].get_energy()>>"
645  + hname, "", "goff");
646  }
647  else if (gain.Contains("RAW"))
648  {
649  if (full_gain)
650  {
651  h2 = new TH2F(hname,
652  Form(";Calibrated Tower Energy Sum (ADC);Count / bin"), 100,
653  .05 * 100, 25 * 100, 64, -.5, 63.5);
654  QAHistManagerDef::useLogBins(h2->GetXaxis());
655  }
656  else
657  {
658  h2 = new TH2F(hname,
659  Form(";Calibrated Tower Energy Sum (ADC);Count / bin"), 260,
660  -.2 * 100, 5 * 100, 64, -.5, 63.5);
661  }
662  T->Draw(
663  "TOWER_" + gain + "_CEMC[].get_bineta() + 8* TOWER_" + gain
664  + "_CEMC[].get_binphi():TOWER_" + gain
665  + "_CEMC[].get_energy()*(-1)>>" + hname, "", "goff");
666  }
667 
668  TText * t;
669  TCanvas *c1 = new TCanvas(
670  "EMCDistribution_" + gain + TString(full_gain ? "_FullGain" : "") + cuts,
671  "EMCDistribution_" + gain + TString(full_gain ? "_FullGain" : "") + cuts,
672  1800, 950);
673  c1->Divide(8, 8, 0., 0.01);
674  int idx = 1;
675  TPad * p;
676 
677  for (int iphi = 8 - 1; iphi >= 0; iphi--)
678  {
679  for (int ieta = 0; ieta < 8; ieta++)
680  {
681  p = (TPad *) c1->cd(idx++);
682  c1->Update();
683 
684  p->SetLogy();
685  if (full_gain)
686  {
687  p->SetLogx();
688  }
689  p->SetGridx(0);
690  p->SetGridy(0);
691 
692  TString hname = Form("hEnergy_ieta%d_iphi%d", ieta, iphi)
693  + TString(full_gain ? "_FullGain" : "");
694 
695  TH1 * h = h2->ProjectionX(hname, ieta + 8 * iphi + 1,
696  ieta + 8 * iphi + 1); // axis bin number is encoded as ieta+8*iphi+1
697 
698  h->SetLineWidth(0);
699  h->SetLineColor(kBlue + 3);
700  h->SetFillColor(kBlue + 3);
701 
702  h->GetXaxis()->SetTitleSize(.09);
703  h->GetXaxis()->SetLabelSize(.08);
704  h->GetYaxis()->SetLabelSize(.08);
705 
706  h->Draw();
707 
708  if (full_gain)
709  h->Fit("x*gaus", "M");
710  else
711  h->Fit("landau", "M");
712 
713  double peak = -1;
714 
715  TF1 * fit = ((TF1 *) (h->GetListOfFunctions()->At(0)));
716  if (fit)
717  {
718 
719  fit->SetLineColor(kRed);
720  peak = fit->GetParameter(1);
721 
722  }
723 
724  cout << Form("Finished <Col%d Row%d> = %.1f", ieta, iphi, peak)
725  << endl;
726 
727  TText *t = new TText(.9, .9,
728  Form("<Col%d Row%d> = %.1f", ieta, iphi, peak));
729  t->SetTextAlign(33);
730  t->SetTextSize(.15);
731  t->SetNDC();
732  t->Draw();
733  }
734  }
735 
736  SaveCanvas(c1,
737  TString(_file0->GetName()) + TString("_DrawPrototype2EMCalTower_")
738  + TString(c1->GetName()), false);
739 
740 }
741 
742 void
743 EMCDistribution_PeakSample_Fast(bool full_gain = false)
744 {
745 
746  const TString gain = "RAW";
747 
748  TString hname = "EMCDistribution_" + gain
749  + TString(full_gain ? "_FullGain" : "") + cuts;
750 
751  TH2 * h2 = NULL;
752  {
753  if (full_gain)
754  {
755  h2 = new TH2F(hname,
756  Form(";Calibrated Tower Energy Sum (ADC);Count / bin"), 100,
757  .05 * 100, 25 * 100, 64, -.5, 63.5);
758  QAHistManagerDef::useLogBins(h2->GetXaxis());
759  }
760  else
761  {
762  h2 = new TH2F(hname,
763  Form(";Calibrated Tower Energy Sum (ADC);Count / bin"), 260,
764  -.2 * 100, 5 * 100, 64, -.5, 63.5);
765  }
766  T->Draw(
767  "TOWER_" + gain + "_CEMC[].get_bineta() + 8* TOWER_" + gain
768  + "_CEMC[].get_binphi():(TOWER_RAW_CEMC[].signal_samples[10] - TOWER_RAW_CEMC[].signal_samples[0])*(-1)>>" + hname, "", "goff");
769  }
770 
771  TText * t;
772  TCanvas *c1 = new TCanvas(
773  "EMCDistribution_PeakSample_Fast_" + TString(full_gain ? "_FullGain" : "") + cuts,
774  "EMCDistribution_PeakSample_Fast_" + TString(full_gain ? "_FullGain" : "") + cuts,
775  1800, 950);
776  c1->Divide(8, 8, 0., 0.01);
777  int idx = 1;
778  TPad * p;
779 
780  for (int iphi = 8 - 1; iphi >= 0; iphi--)
781  {
782  for (int ieta = 0; ieta < 8; ieta++)
783  {
784  p = (TPad *) c1->cd(idx++);
785  c1->Update();
786 
787  p->SetLogy();
788  if (full_gain)
789  {
790  p->SetLogx();
791  }
792  p->SetGridx(0);
793  p->SetGridy(0);
794 
795  TString hname = Form("hEnergy_ieta%d_iphi%d", ieta, iphi)
796  + TString(full_gain ? "_FullGain" : "");
797 
798  TH1 * h = h2->ProjectionX(hname, ieta + 8 * iphi + 1,
799  ieta + 8 * iphi + 1); // axis bin number is encoded as ieta+8*iphi+1
800 
801  h->SetLineWidth(0);
802  h->SetLineColor(kBlue + 3);
803  h->SetFillColor(kBlue + 3);
804 
805  h->GetXaxis()->SetTitleSize(.09);
806  h->GetXaxis()->SetLabelSize(.08);
807  h->GetYaxis()->SetLabelSize(.08);
808 
809  h->Draw();
810 
811  if (full_gain)
812  h->Fit("x*gaus", "M");
813  else
814  h->Fit("landau", "M");
815 
816  double peak = -1;
817 
818  TF1 * fit = ((TF1 *) (h->GetListOfFunctions()->At(0)));
819  if (fit)
820  {
821 
822  fit->SetLineColor(kRed);
823  peak = fit->GetParameter(1);
824 
825  }
826 
827  cout << Form("Finished <Col%d Row%d> = %.1f", ieta, iphi, peak)
828  << endl;
829 
830  TText *t = new TText(.9, .9,
831  Form("<Col%d Row%d> = %.1f", ieta, iphi, peak));
832  t->SetTextAlign(33);
833  t->SetTextSize(.15);
834  t->SetNDC();
835  t->Draw();
836  }
837  }
838 
839  SaveCanvas(c1,
840  TString(_file0->GetName()) + TString("_DrawPrototype2EMCalTower_")
841  + TString(c1->GetName()), false);
842 
843 }
844 
845 
846 
847 void
848 EMCDistribution(TString gain = "CALIB", bool log_scale = false)
849 {
850 
851  TText * t;
852  TCanvas *c1 = new TCanvas(
853  "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "") + cuts,
854  "EMCDistribution_" + gain + TString(log_scale ? "_Log" : "") + cuts, 1800,
855  1000);
856  c1->Divide(8, 8, 0., 0.01);
857  int idx = 1;
858  TPad * p;
859 
860  for (int iphi = 8 - 1; iphi >= 0; iphi--)
861  {
862  for (int ieta = 0; ieta < 8; ieta++)
863  {
864  p = (TPad *) c1->cd(idx++);
865  c1->Update();
866 
867  p->SetLogy();
868  p->SetGridx(0);
869  p->SetGridy(0);
870 
871  TString hname = Form("hEnergy_ieta%d_iphi%d", ieta, iphi)
872  + TString(log_scale ? "_Log" : "");
873 
874  TH1 * h = NULL;
875 
876  if (log_scale)
877  h = new TH1F(hname,
878  Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 300,
879  5e-3, 3096);
880  else
881 // h = new TH1F(hname,
882 // Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 196,
883 // 1900, 2096);
884  h = new TH1F(hname,
885  Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 596,
886  -96, 500);
887 
888  h->SetLineWidth(0);
889  h->SetLineColor(kBlue + 3);
890  h->SetFillColor(kBlue + 3);
891  h->GetXaxis()->SetTitleSize(.09);
892  h->GetXaxis()->SetLabelSize(.08);
893  h->GetYaxis()->SetLabelSize(.08);
894 
895  if (log_scale)
896  QAHistManagerDef::useLogBins(h->GetXaxis());
897 
898  T->Draw(
899  "TOWER_" + gain + "_CEMC[].get_energy_power_law_exp()>>" + hname,
900  Form(
901  "TOWER_%s_CEMC[].get_bineta()==%d && TOWER_%s_CEMC[].get_binphi()==%d",
902  gain.Data(), ieta, gain.Data(), iphi), "");
903 
904  TText *t = new TText(.9, .9, Form("Col%d Row%d", ieta, iphi));
905  t->SetTextAlign(33);
906  t->SetTextSize(.15);
907  t->SetNDC();
908  t->Draw();
909 
910 // return;
911  }
912  }
913 
914  SaveCanvas(c1,
915  TString(_file0->GetName()) + TString("_DrawPrototype2EMCalTower_")
916  + TString(c1->GetName()), false);
917 
918 }
919 
920 void
921 EMCDistribution_ADC(bool log_scale = true)
922 {
923  TString gain = "RAW";
924 
925  TText * t;
926  TCanvas *c1 = new TCanvas(
927  "EMCDistribution_ADC_" + gain + TString(log_scale ? "_Log" : "") + cuts,
928  "EMCDistribution_ADC_" + gain + TString(log_scale ? "_Log" : "") + cuts,
929  1800, 1000);
930  c1->Divide(8, 8, 0., 0.01);
931  int idx = 1;
932  TPad * p;
933 
934  for (int iphi = 8 - 1; iphi >= 0; iphi--)
935  {
936  for (int ieta = 0; ieta < 8; ieta++)
937  {
938  p = (TPad *) c1->cd(idx++);
939  c1->Update();
940 
941  if (log_scale)
942  {
943  p->SetLogz();
944  }
945  p->SetGridx(0);
946  p->SetGridy(0);
947 
948  TString hname = Form("hEnergy_ieta%d_iphi%d", ieta, iphi)
949  + TString(log_scale ? "_Log" : "");
950 
951  TH1 * h = NULL;
952 
953  if (log_scale)
954  h = new TH2F(hname,
955  Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 24, -.5,
956  23.5,
957 // 128+64, 0, 3096);
958  550, 1500, 2050);
959 // else
960 // h = new TH2F(hname,
961 // Form(";Calibrated Tower Energy Sum (GeV);Count / bin"), 100,
962 // -.050, .5,128,0,2048);
963 
964  h->SetLineWidth(0);
965  h->SetLineColor(kBlue + 3);
966  h->SetFillColor(kBlue + 3);
967  h->GetXaxis()->SetTitleSize(.09);
968  h->GetXaxis()->SetLabelSize(.08);
969  h->GetYaxis()->SetLabelSize(.08);
970 
971 // if (log_scale)
972 // QAHistManagerDef::useLogBins(h->GetYaxis());
973 
974  TString sdraw = "TOWER_" + gain
975  + "_CEMC[].signal_samples[]:fmod(Iteration$,24)>>" + hname;
976  TString scut =
977  Form(
978  "TOWER_%s_CEMC[].get_bineta()==%d && TOWER_%s_CEMC[].get_binphi()==%d",
979  gain.Data(), ieta, gain.Data(), iphi);
980 
981  cout << "T->Draw(\"" << sdraw << "\",\"" << scut << "\");" << endl;
982 
983  T->Draw(sdraw, scut, "colz");
984 
985  TText *t = new TText(.9, .9, Form("Col%d Row%d", ieta, iphi));
986  t->SetTextAlign(33);
987  t->SetTextSize(.15);
988  t->SetNDC();
989  t->Draw();
990 
991 // return;
992  }
993  }
994 
995  SaveCanvas(c1,
996  TString(_file0->GetName()) + TString("_DrawPrototype2EMCalTower_")
997  + TString(c1->GetName()), false);
998 
999 }
1000