Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Draw_PHG4DSTReader.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Draw_PHG4DSTReader.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 <cassert>
17 #include "SaveCanvas.C"
18 #include "SetOKStyle.C"
19 using namespace std;
20 
21 TFile * _infile = NULL;
22 TTree * T = NULL;
23 TString cuts = "";
24 
25 void
26 Draw_PHG4DSTReader( //
27 // const TString infile = "G4Hits_sPHENIX_e-_eta0_8GeV.root_DSTReader.root"
28 // const TString infile = "G4Hits_sPHENIX_e-_eta0_8GeV.root_Ana.root_DSTReader.root"//
29 // const TString infile =
30 // "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/single_particle/spacal2d/zerofield/G4Hits_sPHENIX_gamma_etaALL_50GeV.root_DSTReader.root" //
31 // const TString infile = "G4sPHENIXCells.root_DSTReader.root" //
32 // const TString infile = "G4sPHENIXCells_2DSpacal_100e10GeV.root_Ana.root_DSTReader.root"//
33 // const TString infile =
34 // "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_8GeV_CALICEBirk/Sim1096_eneg_DSTReader.root" //
35 // const TString infile =
36 // "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_8GeV/Sum_muonneg.lst_Process.root_DSTReader.root" //
37 // const TString infile = "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/Spacal2D_Leakage/Job_ALL_e-_DSTReader.root"//
38 // const TString infile = "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/Spacal2D_InsPHENIX/Sim1393ALL_eneg_DSTReader.root"//
39  )
40 {
41  SetOKStyle();
42  gStyle->SetOptStat(0);
43  gStyle->SetOptFit(1111);
44  TVirtualFitter::SetDefaultFitter("Minuit2");
45  gSystem->Load("libg4eval.so");
46 
47  if (!_infile)
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  _infile = new TFile;
61  _infile->SetName(infile);
62  }
63 
64  assert(_infile);
65 
66  T->SetAlias("CEMC_Sample",
67  "Sum$(G4HIT_CEMC.light_yield)/(Sum$(G4HIT_CEMC.edep) + Sum$(G4HIT_ABSORBER_CEMC.edep))");
68  T->SetAlias("HCALIN_Sample",
69  "Sum$(G4HIT_HCALIN.light_yield)/(Sum$(G4HIT_HCALIN.edep) + Sum$(G4HIT_ABSORBER_HCALIN.edep))");
70  T->SetAlias("HCALOUT_Sample",
71  "Sum$(G4HIT_HCALOUT.light_yield)/(Sum$(G4HIT_HCALOUT.edep) + Sum$(G4HIT_ABSORBER_HCALOUT.edep))");
72  T->SetAlias("FHCAL_Sample",
73  "Sum$(G4HIT_FHCAL.light_yield)/(Sum$(G4HIT_FHCAL.edep) + Sum$(G4HIT_ABSORBER_FHCAL.edep))");
74 
75  T->SetAlias("PHG4Particle0_e", "0+(PHG4Particle[0].fe)");
76  T->SetAlias("PHG4Particle0_p",
77  "1*sqrt(PHG4Particle[0].fpx**2+PHG4Particle[0].fpy**2+PHG4Particle[0].fpz**2)");
78  T->SetAlias("PHG4Particle0_eta",
79  "0.5*log((PHG4Particle0_p+PHG4Particle[0].fpz)/(PHG4Particle0_p-PHG4Particle[0].fpz))");
80 
81  T->SetAlias("PHG4Particle_p",
82  "1*sqrt(PHG4Particle.fpx**2+PHG4Particle.fpy**2+PHG4Particle.fpz**2)");
83  T->SetAlias("PHG4Particle_eta",
84  "0.5*log((PHG4Particle_p+PHG4Particle.fpz)/(PHG4Particle_p-PHG4Particle.fpz))");
85 
86  T->SetAlias("Entrace_CEMC_x", "G4HIT_ABSORBER_CEMC[0].get_avg_x()+0");
87  T->SetAlias("Entrace_CEMC_y", "G4HIT_ABSORBER_CEMC[0].get_avg_y()+0");
88  T->SetAlias("Entrace_CEMC_z", "G4HIT_ABSORBER_CEMC[0].get_avg_z()+0");
89  T->SetAlias("Entrace_CEMC_azimuthal",
90  "atan2(Entrace_CEMC_y, Entrace_CEMC_x)*95.");
91  T->SetAlias("Average_CEMC_z",
92  "Sum$(G4HIT_ABSORBER_CEMC.get_avg_z()*G4HIT_ABSORBER_CEMC.edep)/Sum$(G4HIT_ABSORBER_CEMC.edep)");
93  T->SetAlias("Average_BH_1_z",
94  "Sum$(G4HIT_BH_1.get_avg_z())/Sum$(G4HIT_BH_1.hit_id>-1)");
95 
96  T->SetAlias("CEMC_E", "Sum$(G4HIT_CEMC.light_yield)+0");
97  T->SetAlias("TOWER_CEMC_E",
98  "Sum$(TOWER_CEMC.light_yield * (TOWER_CEMC.light_yield>32))");
99  T->SetAlias("ABSORBER_CEMC_E", "Sum$(G4HIT_ABSORBER_CEMC.edep)+0");
100  T->SetAlias("Total_CEMC_E",
101  "Sum$(G4HIT_ABSORBER_CEMC.edep)+Sum$(G4HIT_CEMC.edep)");
102  T->SetAlias("HCALIN_E", "Sum$(G4HIT_HCALIN.light_yield)+0");
103  T->SetAlias("HCALOUT_E", "Sum$(G4HIT_HCALOUT.light_yield)+0");
104  T->SetAlias("Total_HCALIN_E",
105  "Sum$(G4HIT_ABSORBER_HCALIN.edep)+Sum$(G4HIT_HCALIN.edep)");
106  T->SetAlias("Total_HCALOUT_E",
107  "Sum$(G4HIT_ABSORBER_HCALOUT.edep)+Sum$(G4HIT_HCALOUT.edep)");
108  T->SetAlias("BH_1_E", "Sum$(G4HIT_BH_1.edep)+0");
109  T->SetAlias("EMCScale", "1*1");
110 
111  T->SetAlias("PHG4Particle0_pT",
112  "1*sqrt(PHG4Particle[0].fpx**2+PHG4Particle[0].fpy**2)");
113  T->SetAlias("PHG4Particle0_theta",
114  "1*atan2(PHG4Particle0_pT,PHG4Particle[0].fpz)");
115  T->SetAlias("PHG4Particle0_phi",
116  "1*atan2(PHG4Particle[0].fpy, PHG4Particle[0].fpx)");
117 
118  T->SetAlias("PHG4Particle0_xhx",
119  "cos(PHG4Particle0_theta)*cos(PHG4Particle0_phi)");
120  T->SetAlias("PHG4Particle0_xhy",
121  "cos(PHG4Particle0_theta)*sin(PHG4Particle0_phi)");
122  T->SetAlias("PHG4Particle0_xhz", "-sin(PHG4Particle0_theta)*1");
123  T->SetAlias("PHG4Particle0_yhx", "1*cos(PHG4Particle0_phi + pi/2)");
124  T->SetAlias("PHG4Particle0_yhy", "1*sin(PHG4Particle0_phi + pi/2)");
125  T->SetAlias("PHG4Particle0_yhz", "0*1");
126  T->SetAlias("PHG4Particle0_zhz", "cos(PHG4Particle0_theta)*1");
127  T->SetAlias("PHG4Particle0_zhy",
128  "sin(PHG4Particle0_theta)*cos(PHG4Particle0_phi)");
129  T->SetAlias("PHG4Particle0_zhz",
130  "sin(PHG4Particle0_theta)*sin(PHG4Particle0_phi)");
131 
132  T->SetAlias("Field", "0*0");
133  TString field_cor = "";
134 // T->SetAlias("Field", "1*1");
135 // TString field_cor = "_field_cor";
136  T->SetAlias("CEMC_yp_cor",
137  "(-19.984 + 5.64862 * PHG4Particle0_pT -0.531359* PHG4Particle0_pT * PHG4Particle0_pT) * (-Field) ");
138  T->SetAlias("ABSORBER_CEMC_yp_cor", "CEMC_yp_cor + 0");
139  T->SetAlias("HCALIN_yp_cor",
140  "(-25.6969 + 6.24025 * PHG4Particle0_pT -0.472037 * PHG4Particle0_pT * PHG4Particle0_pT) * (-Field) ");
141  T->SetAlias("HCALOUT_yp_cor",
142  "(-40.6969 + 6.24025 * PHG4Particle0_pT -0.472037 * PHG4Particle0_pT * PHG4Particle0_pT) * (-Field) ");
143 
144  T->SetAlias("CEMC_xp",
145  "G4HIT_CEMC.get_avg_x() * PHG4Particle0_xhx + G4HIT_CEMC.get_avg_y() * PHG4Particle0_xhy + G4HIT_CEMC.get_avg_z() * PHG4Particle0_xhz");
146  T->SetAlias("CEMC_yp",
147  "CEMC_yp_cor + G4HIT_CEMC.get_avg_x() * PHG4Particle0_yhx + G4HIT_CEMC.get_avg_y() * PHG4Particle0_yhy + G4HIT_CEMC.get_avg_z() * PHG4Particle0_yhz");
148  T->SetAlias("CEMC_zp",
149  "G4HIT_CEMC.get_avg_x() * PHG4Particle0_zhx + G4HIT_CEMC.get_avg_y() * PHG4Particle0_zhy + G4HIT_CEMC.get_avg_z() * PHG4Particle0_zhz");
150  T->SetAlias("CEMC_R", "sqrt(CEMC_xp**2 + CEMC_yp**2)");
151 
152  T->SetAlias("ABSORBER_CEMC_xp",
153  "G4HIT_ABSORBER_CEMC.get_avg_x() * PHG4Particle0_xhx + G4HIT_ABSORBER_CEMC.get_avg_y() * PHG4Particle0_xhy + G4HIT_ABSORBER_CEMC.get_avg_z() * PHG4Particle0_xhz");
154  T->SetAlias("ABSORBER_CEMC_yp",
155  "ABSORBER_CEMC_yp_cor + G4HIT_ABSORBER_CEMC.get_avg_x() * PHG4Particle0_yhx + G4HIT_ABSORBER_CEMC.get_avg_y() * PHG4Particle0_yhy + G4HIT_ABSORBER_CEMC.get_avg_z() * PHG4Particle0_yhz");
156  T->SetAlias("ABSORBER_CEMC_zp",
157  "G4HIT_ABSORBER_CEMC.get_avg_x() * PHG4Particle0_zhx + G4HIT_ABSORBER_CEMC.get_avg_y() * PHG4Particle0_zhy + G4HIT_ABSORBER_CEMC.get_avg_z() * PHG4Particle0_zhz");
158  T->SetAlias("ABSORBER_CEMC_R",
159  "sqrt(ABSORBER_CEMC_xp**2 + ABSORBER_CEMC_yp**2)");
160 
161  T->SetAlias("HCALIN_xp",
162  "G4HIT_HCALIN.get_avg_x() * PHG4Particle0_xhx + G4HIT_HCALIN.get_avg_y() * PHG4Particle0_xhy + G4HIT_HCALIN.get_avg_z() * PHG4Particle0_xhz");
163  T->SetAlias("HCALIN_yp",
164  "HCALIN_yp_cor + G4HIT_HCALIN.get_avg_x() * PHG4Particle0_yhx + G4HIT_HCALIN.get_avg_y() * PHG4Particle0_yhy + G4HIT_HCALIN.get_avg_z() * PHG4Particle0_yhz");
165  T->SetAlias("HCALIN_zp",
166  "G4HIT_HCALIN.get_avg_x() * PHG4Particle0_zhx + G4HIT_HCALIN.get_avg_y() * PHG4Particle0_zhy + G4HIT_HCALIN.get_avg_z() * PHG4Particle0_zhz");
167  T->SetAlias("HCALIN_R", "sqrt(HCALIN_xp**2 + HCALIN_yp**2)");
168 
169  T->SetAlias("HCALOUT_xp",
170  "G4HIT_HCALOUT.get_avg_x() * PHG4Particle0_xhx + G4HIT_HCALOUT.get_avg_y() * PHG4Particle0_xhy + G4HIT_HCALOUT.get_avg_z() * PHG4Particle0_xhz");
171  T->SetAlias("HCALOUT_yp",
172  "HCALOUT_yp_cor + G4HIT_HCALOUT.get_avg_x() * PHG4Particle0_yhx + G4HIT_HCALOUT.get_avg_y() * PHG4Particle0_yhy + G4HIT_HCALOUT.get_avg_z() * PHG4Particle0_yhz");
173  T->SetAlias("HCALOUT_zp",
174  "G4HIT_HCALOUT.get_avg_x() * PHG4Particle0_zhx + G4HIT_HCALOUT.get_avg_y() * PHG4Particle0_zhy + G4HIT_HCALOUT.get_avg_z() * PHG4Particle0_zhz");
175  T->SetAlias("HCALOUT_R", "sqrt(HCALOUT_xp**2 + HCALOUT_yp**2)");
176 
177  T->SetAlias("TOWER_CEMC_eta_mean",
178  "Sum$(TOWER_CEMC.bineta*TOWER_CEMC.get_energy())/Sum$(TOWER_CEMC.get_energy())");
179  T->SetAlias("G4HIT_CEMC_eta_mean",
180  "Sum$(-log(tan(0.5*atan2(sqrt(G4HIT_CEMC.get_avg_x()*G4HIT_CEMC.get_avg_x() + G4HIT_CEMC.get_avg_y()*G4HIT_CEMC.get_avg_y()) , G4HIT_CEMC.get_avg_z()))) * G4HIT_CEMC.edep)/Sum$(G4HIT_CEMC.edep)");
181 
182  T->SetAlias("TOWER_CEMC_phi_mean",
183  "Sum$(TOWER_CEMC.binphi*TOWER_CEMC.get_energy())/Sum$(TOWER_CEMC.get_energy())");
184  T->SetAlias("G4HIT_CEMC_phi_mean",
185  "Sum$(atan2(G4HIT_CEMC.get_avg_y(),G4HIT_CEMC.get_avg_x()) * G4HIT_CEMC.edep)/Sum$(G4HIT_CEMC.edep)");
186 
187  const TCut event_sel = "1*1";
188  cuts = "_all_event";
189 
190  cout << "Build event selection of " << (const char *) event_sel << endl;
191 
192  T->Draw(">>EventList", event_sel);
193  TEventList * elist = gDirectory->GetObjectChecked("EventList", "TEventList");
194  cout << elist->GetN() << " / " << T->GetEntriesFast() << " events selected"
195  << endl;
196 
197  T->SetEventList(elist);
199 // DrawDist(infile);
200 // DrawLeakage(infile);
201 // DrawLeakage_Wide(infile);
202 // DrawLeakage_Phi(infile);
203  Sampling(infile);
204 // DrawTower_Raw_E(infile);
205 
206 // DrawCalibratedE(infile, 0.02206);
207 // DrawCalibratedE(infile, 0.02206*9.71762904052642762e-01);
208 //
209 // DrawCalibratedE(infile, 0.02206*7.71244e+00/7.86120e+00);
210 
211 // DrawCalibratedE_Tower(infile, 500.);
212 // DrawCalibratedE_Tower(infile, 500.*9.71762904052642762e-01);
213 // DrawCalibratedE_Tower(infile, 500.*7.71244e+00/7.86120e+00);
214 
215 // DrawLeakage_LY(infile);
216 //
217 // DrawTowerIDCheck(infile);
218 }
219 
220 void
221 DrawCalibratedE(TString infile, const double SF = 0.021)
222 {
223 
224  TCanvas *c1 = new TCanvas("DrawCalibratedE" + cuts, "DrawCalibratedE" + cuts,
225  1800, 900);
226  c1->Divide(2, 1);
227  int idx = 1;
228  TPad * p;
229 
230  p = (TPad *) c1->cd(idx++);
231  c1->Update();
232  T->Draw("CEMC_E>>hCEMC_E(1000,0,.5)");
233 
234  p = (TPad *) c1->cd(idx++);
235  c1->Update();
236  T->Draw(Form("CEMC_E/%f>>hCEMC_E_SF(1000,0,15)", SF));
237 
238  SaveCanvas(c1,
239  TString(_infile->GetName()) + TString("_DrawJet_")
240  + TString(c1->GetName()), kFALSE);
241 }
242 
243 void
244 DrawCalibratedE_Tower(TString infile, const double SF = 0.021)
245 {
246 
247  TCanvas *c1 = new TCanvas("DrawCalibratedE" + cuts, "DrawCalibratedE" + cuts,
248  1800, 900);
249  c1->Divide(2, 1);
250  int idx = 1;
251  TPad * p;
252 
253  p = (TPad *) c1->cd(idx++);
254  c1->Update();
255  T->Draw("TOWER_CEMC_E>>hCEMC_E(1000,0,8000)");
256 
257  p = (TPad *) c1->cd(idx++);
258  c1->Update();
259  T->Draw(Form("TOWER_CEMC_E/%f>>hCEMC_E_SF(1000,0,15)", SF));
260 
261  SaveCanvas(c1,
262  TString(_infile->GetName()) + TString("_DrawJet_")
263  + TString(c1->GetName()), kFALSE);
264 }
265 
266 void
267 DrawTower_Raw_E(TString infile)
268 {
269  SetOKStyle();
270  gStyle->SetOptStat(0);
271  gStyle->SetOptFit(0);
272  TVirtualFitter::SetDefaultFitter("Minuit2");
273  gSystem->Load("libg4eval.so");
274 
275  TCanvas *c1 = new TCanvas("DrawTower_Raw_E" + cuts, "DrawTower_Raw_E" + cuts,
276  1100, 900);
277  c1->Divide(1, 1);
278  int idx = 1;
279  TPad * p;
280 
281  p = (TPad *) c1->cd(idx++);
282  c1->Update();
283  p->SetLogy();
284  p->SetGridx(0);
285  p->SetGridy(0);
286 
287  T->Draw("TOWER_RAW_CEMC.energy>>he(100,0,25000)");
288 
289  he->SetTitle(";Photo-electron count per tower; A. U.");
290  he->SetLineWidth(4);
291  he->SetLineColor(kBlue + 1);
292 
293  SaveCanvas(c1,
294  TString(_infile->GetName()) + "_Draw_PHG4DSTReader_"
295  + TString(c1->GetName()), true);
296 }
297 
298 void
299 DrawCalibratedE_PlotTestBeam(
300  const TString base_name =
301  "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_Data")
302 {
303  SetOKStyle();
304  gStyle->SetOptStat(0);
305  gStyle->SetOptFit(0);
306  TVirtualFitter::SetDefaultFitter("Minuit2");
307  gSystem->Load("libg4eval.so");
308 
309  const TString config = "nocut";
310 
311  const int N = 8;
312  const double beam_res = 2.7e-2;
313 
314  const double beam_E[8] =
315  { 0.995, 2.026, 3.0, 4.12, 6.0, 8.0, 10.0, 12.0 };
316  double res[100] =
317  { 0 };
318  double res_sub[100] =
319  { 0 };
320  double res_err[100] =
321  { 0 };
322 
323  TCanvas *c1 = new TCanvas(config, config);
324  c1->Print(base_name + "/input_" + config + ".pdf[");
325 
326  for (int i = 0; i < N; i++)
327  {
328  TFile *f = new TFile(
329  base_name + Form("/%.0fGeV_", beam_E[i]) + config + ".root");
330  assert(f);
331 
332  TH1D * ECalSumElectron = (TH1D *) f->GetObjectChecked("ECalSumElectron",
333  "TH1D");
334  TH1D * ECalSumNonElectron = (TH1D *) f->GetObjectChecked(
335  "ECalSumNonElectron", "TH1D");
336  assert(ECalSumElectron);
337  assert(ECalSumNonElectron);
338 // ECalSumElectron->SetStats(false);
339 // ECalSumNonElectron->SetStats(false);
340  TF1 * fgaus =
341  ECalSumElectron->GetListOfFunctions()->At(0);
342  assert(fgaus);
343 
344  ECalSumElectron->Draw();
345  c1->Print(base_name + "/input_" + config + ".pdf");
346 
347  res[i] = fgaus->GetParameter(2)/fgaus->GetParameter(1);
348  res_sub[i] = sqrt(res[i]*res[i] - beam_res*beam_res);
349  res_err[i] = fgaus->GetParError(1)/fgaus->GetParameter(1);
350  }
351  c1->Print(base_name + "/input_" + config + ".pdf]");
352 
353  TGraphErrors * gamma_eta0 = new TGraphErrors(N, beam_E, res_sub, 0, res_err);
354  gamma_eta0->Print();
355  TGraphErrors * gamma_eta9 = new TGraphErrors(N, beam_E, res, 0, res_err);
356  gamma_eta0->Print();
357 
358  TCanvas *c1 = new TCanvas("DrawCalibratedE_PlotTestBeam", "DrawCalibratedE_PlotTestBeam",
359  1100, 900);
360  c1->Divide(1, 1);
361  int idx = 1;
362  TPad * p;
363 
364  p = (TPad *) c1->cd(idx++);
365  c1->Update();
366 
367  p->SetGridx(0);
368  p->SetGridy(0);
369 
370 
371  p->DrawFrame(0, 0, 14, 20e-2,
372  ";Beam Energy (GeV);Relative energy resolution, #DeltaE/E");
373 
374  TLegend * lg = new TLegend(2.*14./35., 9.6e-2, 17.*14./35., 19.6e-2, NULL, "br");
375  TLegend * lg2 = new TLegend(4, 9e-2, 33.*14./35., 19e-2, NULL, "br");
376 
377  TF1 * f_calo = new TF1("f_calo_gamma_eta0", "sqrt([0]*[0]+[1]*[1]/x)/100",
378  0.5, 60);
379  TF1 * f_calo_l = new TF1("f_calo_l_gamma_eta0", "([0]+[1]/sqrt(x))/100", 0.5,
380  60);
381  gamma_eta0->Fit(f_calo, "RM0");
382  f_calo->Print();
383  gamma_eta0->Fit(f_calo_l, "RM0");
384  f_calo_l->Print();
385 
386  TF1 * f_caloR = new TF1("f_calo_gamma_eta0R", "sqrt([0]*[0]+[1]*[1]/x)/100",
387  1.5, 60);
388  TF1 * f_calo_lR = new TF1("f_calo_l_gamma_eta0R", "([0]+[1]/sqrt(x))/100", 1.5,
389  60);
390  gamma_eta0->Fit(f_caloR, "RM0");
391  f_caloR->Print();
392  gamma_eta0->Fit(f_calo_lR, "RM0");
393  f_calo_lR->Print();
394 
395 
396  gamma_eta0->SetLineColor(kBlack);
397  gamma_eta0->SetMarkerColor(kBlack);
398  gamma_eta0->SetLineWidth(2);
399  gamma_eta0->SetMarkerStyle(kFullSquare);
400  gamma_eta0->SetMarkerSize(2);
401 
402  gamma_eta9->SetLineColor(kBlack);
403  gamma_eta9->SetMarkerColor(kBlack);
404  gamma_eta9->SetLineWidth(2);
405  gamma_eta9->SetMarkerStyle(kOpenSquare);
406  gamma_eta9->SetMarkerSize(2);
407 
408  f_calo->SetLineColor(kRed + 1);
409  f_calo->SetLineWidth(2);
410  f_calo_l->SetLineColor(kRed + 1);
411  f_calo_l->SetLineWidth(2);
412  f_calo_l->SetLineStyle(kDashed);
413 
414  f_caloR->SetLineColor(kBlue + 3);
415  f_caloR->SetLineWidth(2);
416  f_calo_lR->SetLineColor(kBlue + 3);
417  f_calo_lR->SetLineWidth(2);
418  f_calo_lR->SetLineStyle(kDashed);
419 
420  f_calo->Draw("same");
421  f_calo_l->Draw("same");
422  f_caloR->Draw("same");
423  f_calo_lR->Draw("same");
424  gamma_eta0->Draw("p");
425  gamma_eta9->Draw("p");
426 
427  lg2->AddEntry(gamma_eta9,
428  Form("Electrons Data, #eta = 0.3-0.4", abs(f_calo->GetParameter(0)),
429  f_calo->GetParameter(1)), "p");
430  lg2->AddEntry(gamma_eta0,
431  Form("Electrons Data - 2.7%% Beam #DeltaE", f_calo->GetParameter(0),
432  f_calo->GetParameter(1)), "p");
433  lg2->AddEntry(f_calo,
434  Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}",abs( f_calo->GetParameter(0)),
435  f_calo->GetParameter(1)), "l");
436  TLegendEntry * entry = lg2->AddEntry(f_calo_l,
437  Form("#DeltaE/E = %.1f%% + %.1f%%/#sqrt{E}", abs(f_calo_l->GetParameter(0)),
438  f_calo_l->GetParameter(1)), "l");
439  entry->SetTextColor(kGray + 1);
440  lg2->AddEntry(f_caloR,
441  Form("#DeltaE/E = %.1f%% #oplus %.1f%%/#sqrt{E}, E #geq 2 GeV", abs(f_caloR->GetParameter(0)),
442  f_caloR->GetParameter(1)), "l");
443  TLegendEntry * entry = lg2->AddEntry(f_calo_lR,
444  Form("#DeltaE/E = %.1f%% + %.1f%%/#sqrt{E}, E #geq 2 GeV",abs( f_calo_lR->GetParameter(0)),
445  f_calo_lR->GetParameter(1)), "l");
446  entry->SetTextColor(kGray + 1);
447 
448 // lg->Draw();
449  lg2->Draw();
450 
451  SaveCanvas(c1, base_name + "/" + TString(c1->GetName()), true);
452 }
453 
454 void
455 //DrawCalibratedE_Plot(TString ID = "4.12GeV", TString config = "")
456 //DrawCalibratedE_Plot(TString ID = "8GeV", TString config = "")
457 DrawCalibratedE_Plot(TString ID = "12GeV",TString config = "")
458 //DrawCalibratedE_Plot(TString ID = "8GeV",TString config = "_Range1um")
459 //DrawCalibratedE_Plot(TString ID = "8GeV",TString config = "_CALICEBirk")
460 //DrawCalibratedE_Plot(TString ID = "4GeV",TString config = "")
461 //DrawCalibratedE_Plot(TString ID = "4GeV",TString config = "_Range1um")
462 //DrawCalibratedE_Plot(TString ID = "4GeV",TString config = "_CALICEBirk")
463 {
464  SetOKStyle();
465  gStyle->SetOptStat(0);
466  gStyle->SetOptFit(0);
467  TVirtualFitter::SetDefaultFitter("Minuit2");
468  gSystem->Load("libg4eval.so");
469 
470 // /direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_8GeV//SimALL_gamma_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root
471 
472 // TFile *f = new TFile("/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"+ID+config+"//SimALL_eneg_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
473 // assert(f);
474 // TH1F * h_eneg = (TH1F *)f->GetObjectChecked("hCEMC_E_SF","TH1F");
475 // assert(h_eneg);
476 // TFile *f = new TFile("/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"+ID+config+"//SimALL_pineg_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
477 // assert(f);
478 // TH1F * h_pineg = (TH1F *)f->GetObjectChecked("hCEMC_E_SF","TH1F");
479 // assert(h_pineg);
480 // TFile *f = new TFile("/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"+ID+config+"//SimALL_kaonneg_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
481 // assert(f);
482 // TH1F * h_kaonneg = (TH1F *)f->GetObjectChecked("hCEMC_E_SF","TH1F");
483 // assert(h_kaonneg);
484 
485  TFile *f =
486  new TFile(
487  "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"
488  + ID + config
489  + "//Sum_eneg.lst_Process.root_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
490  assert(f);
491  TH1F * h_eneg = (TH1F *) f->GetObjectChecked("hCEMC_E_SF", "TH1F");
492  assert(h_eneg);
493  TFile *f =
494  new TFile(
495  "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"
496  + ID + config
497  + "//Sum_pineg.lst_Process.root_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
498  assert(f);
499  TH1F * h_pineg = (TH1F *) f->GetObjectChecked("hCEMC_E_SF", "TH1F");
500  assert(h_pineg);
501  TFile *f =
502  new TFile(
503  "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"
504  + ID + config
505  + "//Sum_kaonneg.lst_Process.root_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
506  assert(f);
507  TH1F * h_kaonneg = (TH1F *) f->GetObjectChecked("hCEMC_E_SF", "TH1F");
508  assert(h_kaonneg);
509  TFile *f =
510  new TFile(
511  "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"
512  + ID + config
513  + "//Sum_muonneg.lst_Process.root_DSTReader.root_DrawJet_DrawCalibratedE_all_event.root");
514  assert(f);
515  TH1F * h_muonneg = (TH1F *) f->GetObjectChecked("hCEMC_E_SF", "TH1F");
516  assert(h_muonneg);
517 
518  TH1F * h_eneg_scale_hadron_data_tail = (TH1F *) h_eneg->Clone(
519  "h_eneg_scale_hadron_data_tail");
520 
521  TFile *f =
522  new TFile(
523  "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_8GeV/beam_data/"
524  + ID + ".root");
525  assert(f);
526 
527  TH1D * ECalSumElectron = (TH1D *) f->GetObjectChecked("ECalSumElectron",
528  "TH1D");
529  TH1D * ECalSumNonElectron = (TH1D *) f->GetObjectChecked("ECalSumNonElectron",
530  "TH1D");
531  assert(ECalSumElectron);
532  assert(ECalSumNonElectron);
533  ECalSumElectron->SetStats(false);
534  ECalSumNonElectron->SetStats(false);
535  ECalSumElectron->GetListOfFunctions()->RemoveAt(0);
536 
537  const double sim_E = 7.86131e+00;
538  const double data_E = 1.09505e+03;
539 
540  ECalSumElectron->GetXaxis()->Set(ECalSumElectron->GetXaxis()->GetNbins(),
541  ECalSumElectron->GetXaxis()->GetXmin() * sim_E / data_E,
542  ECalSumElectron->GetXaxis()->GetXmax() * sim_E / data_E //
543  );
544  ECalSumNonElectron->GetXaxis()->Set(
545  ECalSumNonElectron->GetXaxis()->GetNbins(),
546  ECalSumNonElectron->GetXaxis()->GetXmin() * sim_E / data_E,
547  ECalSumNonElectron->GetXaxis()->GetXmax() * sim_E / data_E //
548  );
549 
550  const double energy = h_eneg->GetMean();
551 
552  double eneg_scale = ECalSumElectron->Integral(
553  ECalSumElectron->GetXaxis()->FindBin(energy * .7),
554  ECalSumElectron->GetXaxis()->FindBin(energy * 1.3))
555  / h_eneg->Integral(h_eneg->GetXaxis()->FindBin(energy * .7),
556  h_eneg->GetXaxis()->FindBin(energy * 1.3));
557 
558  const double hadron_norm_max_range = 0.7;
559 
560  double muonneg_scale = 0.1 * //10% muon in hadron sample
561  ECalSumNonElectron->Integral(ECalSumNonElectron->GetXaxis()->FindBin(.01),
562  ECalSumNonElectron->GetXaxis()->FindBin(energy * 1.3))
563  / h_muonneg->Integral(h_muonneg->GetXaxis()->FindBin(.01),
564  h_muonneg->GetXaxis()->FindBin(energy * 1.3));
565 
566  const double muon_ratio_in_MIP = h_muonneg->Integral(
567  h_muonneg->GetXaxis()->FindBin(.01),
568  h_muonneg->GetXaxis()->FindBin(hadron_norm_max_range))
569  / h_muonneg->Integral(h_muonneg->GetXaxis()->FindBin(.01),
570  h_muonneg->GetXaxis()->FindBin(energy * 1.3));
571  const double muon_count_in_MIP = muon_ratio_in_MIP //
572  * 0.1 * //10% muon in hadron sample
573  ECalSumNonElectron->Integral(ECalSumNonElectron->GetXaxis()->FindBin(.01),
574  ECalSumNonElectron->GetXaxis()->FindBin(energy * 1.3));
575 
576  cout << "DrawCalibratedE_Plot: muon_ratio_in_MIP = " << muon_ratio_in_MIP
577  << ", muon_count_in_MIP = " << muon_count_in_MIP << endl;
578 
579  double pineg_scale = (ECalSumNonElectron->Integral(
580  ECalSumNonElectron->GetXaxis()->FindBin(.01),
581  ECalSumNonElectron->GetXaxis()->FindBin(hadron_norm_max_range))
582  - muon_count_in_MIP)
583  / h_pineg->Integral(h_pineg->GetXaxis()->FindBin(.01),
584  h_pineg->GetXaxis()->FindBin(hadron_norm_max_range));
585 
586  double kaonneg_scale = (ECalSumNonElectron->Integral(
587  ECalSumNonElectron->GetXaxis()->FindBin(.01),
588  ECalSumNonElectron->GetXaxis()->FindBin(hadron_norm_max_range))
589  - muon_count_in_MIP)
590  / h_kaonneg->Integral(h_kaonneg->GetXaxis()->FindBin(.01),
591  h_kaonneg->GetXaxis()->FindBin(hadron_norm_max_range));
592 
593  double eneg_scale_hadron_data_tail = ECalSumNonElectron->Integral(
594  ECalSumNonElectron->GetXaxis()->FindBin(energy * 0.95),
595  ECalSumNonElectron->GetXaxis()->FindBin(energy * 1.15))
596  / h_eneg_scale_hadron_data_tail->Integral(
597  h_eneg_scale_hadron_data_tail->GetXaxis()->FindBin(energy * 0.95),
598  h_eneg_scale_hadron_data_tail->GetXaxis()->FindBin(energy * 1.15));
599 
600  h_eneg->Rebin(5);
601  h_pineg->Rebin(5);
602  h_kaonneg->Rebin(5);
603  h_muonneg->Rebin(5);
604  h_eneg_scale_hadron_data_tail->Rebin(5);
605 
606 // ECalSumElectron->Rebin();
607  ECalSumNonElectron->Rebin(40);
608 
609  ECalSumElectron->Sumw2();
610  ECalSumNonElectron->Sumw2();
611 
612  ECalSumElectron->SetMarkerStyle(kFullCircle);
613  ECalSumElectron->SetMarkerSize(1);
614  ECalSumNonElectron->SetMarkerStyle(kFullCircle);
615  ECalSumNonElectron->SetMarkerSize(1);
616 
617  eneg_scale *= ECalSumElectron->GetBinWidth(1) / h_eneg->GetBinWidth(1);
618  pineg_scale *= ECalSumNonElectron->GetBinWidth(1) / h_pineg->GetBinWidth(1);
619  kaonneg_scale *= ECalSumNonElectron->GetBinWidth(1)
620  / h_kaonneg->GetBinWidth(1);
621  muonneg_scale *= ECalSumNonElectron->GetBinWidth(1)
622  / h_muonneg->GetBinWidth(1);
623  eneg_scale_hadron_data_tail *= ECalSumNonElectron->GetBinWidth(1)
624  / h_eneg_scale_hadron_data_tail->GetBinWidth(1);
625 
626  //Resolution for electron sample
627 
628  TF1 * f_gaus_sim = new TF1("f_gaus_sim", "gaus",
629  h_eneg->GetMean() - 2 * h_eneg->GetRMS(),
630  h_eneg->GetMean() + 2 * h_eneg->GetRMS());
631  f_gaus_sim->SetParameters(1, h_eneg->GetMean(), h_eneg->GetRMS());
632  new TCanvas;
633  h_eneg->Fit(f_gaus_sim, "RM0");
634  h_eneg->DrawClone();
635  f_gaus_sim->Draw("same");
636 
637  TF1 * f_gaus_data = new TF1("f_gaus_data", "gaus",
638  h_eneg->GetMean() - 2 * h_eneg->GetRMS(),
639  h_eneg->GetMean() + 2 * h_eneg->GetRMS());
640  f_gaus_data->SetParameters(1, h_eneg->GetMean(), h_eneg->GetRMS());
641  new TCanvas;
642  ECalSumElectron->Fit(f_gaus_data, "RM0");
643  ECalSumElectron->DrawClone();
644  f_gaus_data->Draw("same");
645 
646  h_eneg->Scale(eneg_scale);
647  h_pineg->Scale(pineg_scale);
648  h_kaonneg->Scale(kaonneg_scale);
649  h_muonneg->Scale(muonneg_scale);
650  h_eneg_scale_hadron_data_tail->Scale(eneg_scale_hadron_data_tail);
651 
652 // h_pineg -> Add(h_eneg_scale_hadron_data_tail);
653 // h_kaonneg -> Add(h_eneg_scale_hadron_data_tail);
654 
655  TCanvas *c1 = new TCanvas("DrawCalibratedE_Plot" + cuts,
656  "DrawCalibratedE_Plot" + cuts, 1800, 700);
657  c1->Divide(2, 1);
658  int idx = 1;
659  TPad * p;
660 
661  p = (TPad *) c1->cd(idx++);
662  c1->Update();
663  p->SetLogy();
664  p->SetGridx(0);
665  p->SetGridy(0);
666 
667  h_eneg->SetLineWidth(2);
668 
669  ECalSumElectron->GetXaxis()->SetRangeUser(0, energy * 1.5);
670  ECalSumElectron->Draw();
671  h_eneg->SetLineWidth(3);
672  h_eneg->SetLineColor(kGreen + 3);
673  h_eneg->SetFillColor(kGreen - 7);
674  h_eneg->SetFillStyle(1001);
675  h_eneg->Draw("same");
676  ECalSumElectron->Draw("same");
677 
678  ECalSumElectron->SetTitle(";Calibrated Energy (GeV);Count/bin");
679 
680  TLegend * lg = new TLegend(0.13, 0.65, 0.55, 0.85);
681  lg->AddEntry(ECalSumElectron,
682  Form("#splitline{2014 eRD1 beam test (e-cut)}{#DeltaE/E = %.1f%%}",
683  f_gaus_data->GetParameter(2) / f_gaus_data->GetParameter(1) * 100),
684  "pe");
685  lg->AddEntry(h_eneg,
686  Form("#splitline{e^{-} sim using sPHENIX Geant4,}{#DeltaE/E = %.1f%%}",
687  f_gaus_sim->GetParameter(2) / f_gaus_sim->GetParameter(1) * 100),
688  "lf");
689  lg->Draw();
690 
691  p = (TPad *) c1->cd(idx++);
692  c1->Update();
693  p->SetLogy();
694  p->SetGridx(0);
695  p->SetGridy(0);
696 
697  h_pineg->SetLineColor(kRed);
698  h_pineg->SetLineWidth(3);
699  h_kaonneg->SetLineColor(kBlue);
700  h_kaonneg->SetLineWidth(3);
701  h_muonneg->SetLineColor(kBlack);
702  h_muonneg->SetFillColor(kGray);
703  h_muonneg->SetLineWidth(3);
704  h_muonneg->SetFillStyle(1001);
705  h_eneg_scale_hadron_data_tail->SetLineColor(kGreen + 3);
706  h_eneg_scale_hadron_data_tail->SetLineWidth(2);
707  h_eneg_scale_hadron_data_tail->SetFillColor(kGreen - 7);
708  h_eneg_scale_hadron_data_tail->SetFillStyle(1001);
709 
710  ECalSumNonElectron->GetXaxis()->SetRangeUser(0, energy * 1.5);
711 
712  ECalSumNonElectron->Draw();
713  h_muonneg->Draw("same");
714 // h_eneg_scale_hadron_data_tail->Draw("same"); ////////////// Guess of electron contamination
715  h_pineg->Draw("same");
716  h_kaonneg->Draw("same");
717  ECalSumNonElectron->Draw("same");
718 
719  ECalSumNonElectron->SetTitle(";Calibrated Energy (GeV);Count/bin");
720 
721  TLegend * lg = new TLegend(0.35, 0.65, 0.85, 0.9);
722  lg->AddEntry(ECalSumNonElectron, Form("2014 eRD1 beam test (#bar{e-cut})"),
723  "pe");
724  lg->AddEntry(h_pineg, Form("#pi^{-} sim using sPHENIX Geant4"), "l");
725  lg->AddEntry(h_kaonneg, Form("K^{-} sim using sPHENIX Geant4"), "l");
726 // lg->AddEntry(h_eneg_scale_hadron_data_tail, Form("e^{-} sim using sPHENIX Geant4"),////////////// Guess of electron contamination
727 // "lf");////////////// Guess of electron contamination
728  lg->AddEntry(h_muonneg, Form("#mu^{-} (~10%% data) sim using sPHENIX Geant4"), "lf");
729  lg->Draw();
730 
731  SaveCanvas(c1,
732  TString(
733  "/direct/phenix+sim02/phnxreco/ePHENIX/jinhuang/sPHENIX_work/SPACAL_TestBeam_"
734  + ID + config + "/DrawCalibratedE_Plot_DrawJet_")
735  + TString(c1->GetName()), true);
736 }
737 
738 void
739 DrawTowerIDCheck(TString infile)
740 {
741 
742  TCanvas *c1 = new TCanvas("DrawTowerIDCheck" + cuts,
743  "DrawTowerIDCheck" + cuts, 1800, 900);
744  c1->Divide(2, 1);
745  int idx = 1;
746  TPad * p;
747 
748  p = (TPad *) c1->cd(idx++);
749  c1->Update();
750 
751  T->Draw("TOWER_CEMC_eta_mean:G4HIT_CEMC_eta_mean", "", "*");
752 
753  p = (TPad *) c1->cd(idx++);
754  c1->Update();
755 
756  T->Draw("TOWER_CEMC_phi_mean:G4HIT_CEMC_phi_mean", "", "*");
757 
758  SaveCanvas(c1,
759  TString(_infile->GetName()) + TString("_DrawJet_")
760  + TString(c1->GetName()), kFALSE);
761 }
762 
763 void
764 DrawDist(TString infile)
765 {
766 
767  TCanvas *c1 = new TCanvas("DrawDist" + cuts, "DrawDist" + cuts, 1800, 900);
768  c1->Divide(2, 1);
769  int idx = 1;
770  TPad * p;
771 
772  p = (TPad *) c1->cd(idx++);
773  c1->Update();
774 
775  T->Draw(
776  "G4HIT_CEMC.get_avg_y():G4HIT_CEMC.get_avg_x()>>h_EMC_Azim(600,-300,300,600,-300,300)",
777  "G4HIT_CEMC.edep", "goff");
778 // T->Draw(
779 // "G4HIT_HCALIN.get_avg_y():G4HIT_HCALIN.get_avg_x()>>h_HCALIN_Azim(600,-300,300,600,-300,300)",
780 // "G4HIT_HCALIN.edep", "goff");
781 // T->Draw(
782 // "G4HIT_HCALOUT.get_avg_y():G4HIT_HCALOUT.get_avg_x()>>h_HCALOUT_Azim(600,-300,300,600,-300,300)",
783 // "G4HIT_HCALOUT.edep", "goff");
784 // TGraph * g_EMC =
785 // (TGraph *) gPad->GetListOfPrimitives()->FindObject("Graph")->Clone(
786 // "g_EMC");
787 // assert(g_EMC);
788 // g_EMC->SetName("g_EMC");
789 // g_EMC->SetMarkerColor(kRed);
790 // g_EMC->SetMarkerStyle(kDot);
791 
792  p->DrawFrame(-300, -300, 300, 300,
793  "Azimuthal distribution of hits;X (cm);Y (cm)");
794  h_EMC_Azim->Draw("colz same");
795 // h_HCALIN_Azim->Draw("colz same");
796 // h_HCALOUT_Azim->Draw("colz same");
797 
798  TEllipse * cBaBar = new TEllipse(0, 0, 140);
799 
800  cBaBar->SetFillStyle(0);
801  cBaBar->SetLineWidth(3);
802  cBaBar->SetLineColor(kMagenta);
803  cBaBar->Draw();
804  TEllipse * cBaBar = new TEllipse(0, 0, 170);
805 
806  cBaBar->SetFillStyle(0);
807  cBaBar->SetLineWidth(3);
808  cBaBar->SetLineColor(kMagenta);
809  cBaBar->Draw();
810 
811  p = (TPad *) c1->cd(idx++);
812  c1->Update();
813 
814  T->Draw(
815  "sqrt(G4HIT_CEMC.get_avg_x()**2 + G4HIT_CEMC.get_avg_y()**2):G4HIT_CEMC.get_avg_z()>>h_CEMC_Rz(600,-300,300,300,0,300)",
816  "G4HIT_CEMC.edep", "goff");
817 // T->Draw(
818 // "sqrt(G4HIT_HCALIN.get_avg_x()**2 + G4HIT_HCALIN.get_avg_y()**2):G4HIT_HCALIN.get_avg_z()>>h_HCALIN_Rz(600,-300,300,300,0,300)",
819 // "G4HIT_HCALIN.edep", "goff");
820 // T->Draw(
821 // "sqrt(G4HIT_HCALOUT.get_avg_x()**2 + G4HIT_HCALOUT.get_avg_y()**2):G4HIT_HCALOUT.get_avg_z()>>h_HCALOUT_Rz(600,-300,300,300,0,300)",
822 // "G4HIT_HCALOUT.edep", "goff");
823 
824 // TGraph * g_EMC =
825 // (TGraph *) gPad->GetListOfPrimitives()->FindObject("Graph")->Clone(
826 // "g_EMC");
827 // assert(g_EMC);
828 // g_EMC->SetName("g_EMC");
829 // g_EMC->SetMarkerColor(kRed);
830 // g_EMC->SetMarkerStyle(kDot);
831 
832  p->DrawFrame(-300, 00, 300, 300,
833  "Azimuthal distribution of hits;Z (cm);R (cm)");
834  h_CEMC_Rz->Draw("colz same");
835 // h_HCALIN_Rz->Draw("colz same");
836 // h_HCALOUT_Rz->Draw("colz same");
837 
838  TLine * lBaBar = new TLine(-174.5, 140, 174.5, 140);
839 
840  lBaBar->SetLineWidth(3);
841  lBaBar->SetLineColor(kMagenta);
842  lBaBar->Draw();
843  TLine * lBaBar = new TLine(-174.5, 170, 174.5, 170);
844 
845  lBaBar->SetLineWidth(3);
846  lBaBar->SetLineColor(kMagenta);
847  lBaBar->Draw();
848 
849  SaveCanvas(c1,
850  TString(_infile->GetName()) + TString("_DrawJet_")
851  + TString(c1->GetName()), kFALSE);
852 }
853 
854 void
855 DrawLeakage(TString infile)
856 {
857 
858  TCanvas *c1 = new TCanvas("DrawLeakage" + cuts, "DrawLeakage" + cuts, 1800,
859  900);
860  c1->Divide(2, 3);
861  int idx = 1;
862  TPad * p;
863 
864  p = (TPad *) c1->cd(idx++);
865  c1->Update();
866  T->Draw(
867  "Average_CEMC_z:PHG4Particle0_eta>>hAverage_CEMC_z_PHG4Particle0_eta(120,0.8,0.9,120,75,125)",
868  "", "colz");
869 
870  p = (TPad *) c1->cd(idx++);
871  c1->Update();
872  T->Draw(
873  "Entrace_CEMC_z:PHG4Particle0_eta>>hEntrace_CEMC_z_PHG4Particle0_eta(120,0.8,0.9,120,75,125)",
874  "", "colz");
875 
876  p = (TPad *) c1->cd(idx++);
877  c1->Update();
878 
879  T->Draw(
880  "Total_HCALIN_E/PHG4Particle0_e:Average_CEMC_z>>hLeakage_Average_CEMC_z(120,75,125,100,0,0.25)",
881  "", "colz");
882 
883  p = (TPad *) c1->cd(idx++);
884  c1->Update();
885 
886  T->Draw(
887  "Total_HCALIN_E/PHG4Particle0_e:PHG4Particle0_eta>>hLeakage_PHG4Particle0_eta(120,0.8,0.9,100,0,0.25)",
888  "", "colz");
889 
890  p = (TPad *) c1->cd(idx++);
891  c1->Update();
892  T->Draw(
893  "Total_CEMC_E/PHG4Particle0_e:Average_CEMC_z>>hTotal_CEMC_E_Entrance_Z_Entrace_CEMC_z(120,75,125,100,0,1.2)",
894  "", "colz");
895 
896  p = (TPad *) c1->cd(idx++);
897  c1->Update();
898  T->Draw(
899  "Total_CEMC_E/PHG4Particle0_e:PHG4Particle0_eta>>hTotal_CEMC_E_PHG4Particle0_eta(120,0.8,0.9,100,0,1.2)",
900  "", "colz");
901 
902  SaveCanvas(c1,
903  TString(_infile->GetName()) + TString("_DrawJet_")
904  + TString(c1->GetName()), kFALSE);
905 }
906 void
907 DrawLeakage_LY(TString infile)
908 {
909 
910  TCanvas *c1 = new TCanvas("DrawLeakage_LY" + cuts, "DrawLeakage_LY" + cuts,
911  1800, 900);
912  c1->Divide(2, 3);
913  int idx = 1;
914  TPad * p;
915 
916  p = (TPad *) c1->cd(idx++);
917  c1->Update();
918  T->Draw(
919  "Sum$(G4HIT_CEMC.light_yield)/PHG4Particle0_e:Average_CEMC_z>>hTotal_CEMC_E_Entrance_Z_Entrace_CEMC_z(120,75,125,100,0,.05)",
920  "", "colz");
921 
922  p = (TPad *) c1->cd(idx++);
923  c1->Update();
924  T->Draw(
925  "Sum$(G4HIT_CEMC.light_yield)/PHG4Particle0_e:PHG4Particle0_eta>>hTotal_CEMC_E_PHG4Particle0_eta(120,0.8,0.9,100,0,.05)",
926  "", "colz");
927 
928  SaveCanvas(c1,
929  TString(_infile->GetName()) + TString("_DrawJet_")
930  + TString(c1->GetName()), kFALSE);
931 }
932 
933 void
934 DrawLeakage_Wide(TString infile)
935 {
936 
937  TCanvas *c1 = new TCanvas("DrawLeakage_Wide" + cuts,
938  "DrawLeakage_Wide" + cuts, 1800, 900);
939  c1->Divide(2, 3);
940  int idx = 1;
941  TPad * p;
942 
943  p = (TPad *) c1->cd(idx++);
944  c1->Update();
945  T->Draw(
946  "Average_CEMC_z:PHG4Particle0_eta>>hAverage_CEMC_z_PHG4Particle0_eta(300,0,1.1,300,0,150)",
947  "", "colz");
948 
949  p = (TPad *) c1->cd(idx++);
950  c1->Update();
951  T->Draw(
952  "Entrace_CEMC_z:PHG4Particle0_eta>>hEntrace_CEMC_z_PHG4Particle0_eta(300,0,1.1,300,0,150)",
953  "", "colz");
954 
955  p = (TPad *) c1->cd(idx++);
956  c1->Update();
957 
958  T->Draw(
959  "Total_HCALIN_E/PHG4Particle0_e:Average_CEMC_z>>hLeakage_Average_CEMC_z(300,0,150,100,0,0.25)",
960  "", "colz");
961 
962  p = (TPad *) c1->cd(idx++);
963  c1->Update();
964 
965  T->Draw(
966  "Total_HCALIN_E/PHG4Particle0_e:PHG4Particle0_eta>>hLeakage_PHG4Particle0_eta(300,0,1.1,100,0,0.25)",
967  "", "colz");
968 
969  p = (TPad *) c1->cd(idx++);
970  c1->Update();
971  T->Draw(
972  "Total_CEMC_E/PHG4Particle0_e:Average_CEMC_z>>hTotal_CEMC_E_Entrance_Z_Entrace_CEMC_z(300,0,150,100,0,1.2)",
973  "", "colz");
974 
975  p = (TPad *) c1->cd(idx++);
976  c1->Update();
977  T->Draw(
978  "Total_CEMC_E/PHG4Particle0_e:PHG4Particle0_eta>>hTotal_CEMC_E_PHG4Particle0_eta(300,0,1.1,100,0,1.2)",
979  "", "colz");
980 
981  SaveCanvas(c1,
982  TString(_infile->GetName()) + TString("_DrawJet_")
983  + TString(c1->GetName()), kFALSE);
984 }
985 
986 void
987 DrawLeakage_Phi(TString infile)
988 {
989 
990  TCanvas *c1 = new TCanvas("DrawLeakage_Phi" + cuts, "DrawLeakage_Phi" + cuts,
991  1800, 900);
992  c1->Divide(1, 2);
993  int idx = 1;
994  TPad * p;
995 
996  p = (TPad *) c1->cd(idx++);
997  c1->Update();
998 
999  T->Draw(
1000  "Total_HCALIN_E/PHG4Particle0_e:PHG4Particle0_phi>>hLeakage_PHG4Particle0_e(320,-3.14159,3.14159,100,0,0.25)",
1001  "", "colz");
1002 
1003  p = (TPad *) c1->cd(idx++);
1004  c1->Update();
1005  T->Draw(
1006  "Total_CEMC_E/PHG4Particle0_e:PHG4Particle0_phi>>hTotal_CEMC_E_PHG4Particle0_e(320,-3.14159,3.14159,100,0,1.2)",
1007  "", "colz");
1008 
1009  SaveCanvas(c1,
1010  TString(_infile->GetName()) + TString("_DrawJet_")
1011  + TString(c1->GetName()), kFALSE);
1012 }
1013 
1014 void
1015 DrawLeakage_Phi_Folding(const int N_fold = 16)
1016 {
1017  SetOKStyle();
1018  gStyle->SetOptStat(0);
1019  gStyle->SetOptFit(1111);
1020  TVirtualFitter::SetDefaultFitter("Minuit2");
1021  gSystem->Load("libg4eval.so");
1022 
1023  TH2F * hLeakage_PHG4Particle0_e = (TH2F *) gDirectory->GetObjectChecked(
1024  "hLeakage_PHG4Particle0_e", "TH2F");
1025  TH2F * hTotal_CEMC_E_PHG4Particle0_e = (TH2F *) gDirectory->GetObjectChecked(
1026  "hTotal_CEMC_E_PHG4Particle0_e", "TH2F");
1027 
1028  TH2F * hLeakage_PHG4Particle0_e_Fold = new TH2F(
1029  "hLeakage_PHG4Particle0_e_Fold", ";Phi (rad);Leakage Ratio", 320 / N_fold,
1030  -3.14159 / N_fold, 3.14159 / N_fold, 100, 0, 0.25);
1031  TH2F * hTotal_CEMC_E_PHG4Particle0_e_Fold = new TH2F(
1032  "hTotal_CEMC_E_PHG4Particle0_e_Fold", ";Phi (rad);EMCal Energy Ratio",
1033  320 / N_fold, -3.14159 / N_fold, 3.14159 / N_fold, 100, 0, 1.2);
1034 
1035  for (int biny = 1; biny <= hLeakage_PHG4Particle0_e->GetNbinsY(); biny++)
1036  for (int bin = 1; bin <= hLeakage_PHG4Particle0_e->GetNbinsX(); bin++)
1037  {
1038  const int folded_bin = (bin - 1) % (320 / N_fold) + 1;
1039 
1040  hLeakage_PHG4Particle0_e_Fold->SetBinContent(folded_bin, biny,
1041  hLeakage_PHG4Particle0_e_Fold->GetBinContent(folded_bin, biny)
1042  + hLeakage_PHG4Particle0_e->GetBinContent(bin, biny));
1043  hTotal_CEMC_E_PHG4Particle0_e_Fold->SetBinContent(folded_bin, biny,
1044  hTotal_CEMC_E_PHG4Particle0_e_Fold->GetBinContent(folded_bin, biny)
1045  + hTotal_CEMC_E_PHG4Particle0_e->GetBinContent(bin, biny));
1046 
1047  }
1048 
1049  TCanvas *c1 = new TCanvas("DrawLeakage_Phi_Folding" + cuts,
1050  "DrawLeakage_Phi_Folding" + cuts, 1000, 900);
1051  c1->Divide(1, 2);
1052  int idx = 1;
1053  TPad * p;
1054 
1055  p = (TPad *) c1->cd(idx++);
1056  c1->Update();
1057  hLeakage_PHG4Particle0_e_Fold->Draw("colz");
1058 
1059  p = (TPad *) c1->cd(idx++);
1060  c1->Update();
1061  hTotal_CEMC_E_PHG4Particle0_e_Fold->Draw("colz");
1062 
1063  SaveCanvas(c1, TString("_DrawJet_") + TString(c1->GetName()), kFALSE);
1064 }
1065 
1066 void
1067 Sampling(TString infile)
1068 {
1069 
1070  TCanvas *c1 = new TCanvas("Sampling" + cuts, "Sampling" + cuts, 900, 900);
1071  c1->Divide(1, 1);
1072  int idx = 1;
1073  TPad * p = gPad;
1074 
1075  T->Draw(
1076  "Sum$(G4HIT_CEMC.light_yield):(Sum$(G4HIT_CEMC.edep) + Sum$(G4HIT_ABSORBER_CEMC.edep))",
1077  "", "*");
1078  TGraph * gh2D =
1079  (TGraph *) gPad->GetListOfPrimitives()->FindObject("Graph")->Clone(
1080  "gh2D");
1081  assert(gh2D);
1082  gh2D->SetName("gh2D");
1083  gh2D->SetMarkerStyle(kDot);
1084 
1085  p->DrawFrame(0, 0, 10, 1,
1086  "EMC Energy Deposition;Absorber+Scintillator (GeV);Scintillator (GeV)");
1087  gh2D->Draw("p");
1088 
1089  TF1 * fsampling = new TF1("fsampling", "[0]*x", 0.4, 100);
1090  gh2D->Fit(fsampling, "MR");
1091 
1092  SaveCanvas(c1,
1093  TString(_infile->GetName()) + TString("_DrawJet_")
1094  + TString(c1->GetName()), kFALSE);
1095 }