Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DrawEcal_Likelihood.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DrawEcal_Likelihood.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 * _file0 = NULL;
22 TTree * T = NULL;
23 TString cuts = "";
24 
25 void
26 DrawEcal_Likelihood(
27 //
28  TString base_dir =
29  "../..//sPHENIX_work/production_analysis_cemc2x2/emcstudies/pidstudies/spacal2d/fieldmap/",
30 // TString base_dir =
31 // "../../sPHENIX_work/production_analysis_cemc2x2/embedding/emcstudies/pidstudies/spacal2d/fieldmap/",
32 // TString pid = "anti_proton", //
33  TString pid = "e-", //
34  TString kine_config = "eta0_8GeV", int eval_mode = 1)
35 {
36 
37  const TString infile = base_dir + "G4Hits_sPHENIX_" + pid + "_" + kine_config
38  + "-ALL.root_Ana.root.lst_EMCalLikelihood.root";
39 // "../../sPHENIX_work/production_analysis/emcstudies/pidstudies/spacal2d/fieldmap/G4Hits_sPHENIX_pi-_eta0_8GeV-ALL.root_Ana.root.lst_EMCalLikelihood.root"//
40 
41  SetOKStyle();
42  gStyle->SetOptStat(0);
43  gStyle->SetOptFit(1111);
44  TVirtualFitter::SetDefaultFitter("Minuit2");
45 
46  gSystem->Load("libg4eval.so");
47  gSystem->Load("libemcal_ana.so");
48  gSystem->Load("libg4vertex.so");
49 
50  if (!_file0)
51  {
52  TString chian_str = infile;
53  chian_str.ReplaceAll("ALL", "*");
54  chian_str.ReplaceAll("+", "\\+");
55 
56 
57  TChain * t = new TChain("T");
58  const int n = t->Add(chian_str);
59 
60  cout << "Loaded " << n << " root files with " << chian_str << endl;
61  assert(n>0);
62 
63  T = t;
64 
65  _file0 = new TFile;
66  _file0->SetName(infile);
67 
68  fstream flst(infile + ".lst", ios_base::out);
69 
70  TObjArray *fileElements = t->GetListOfFiles();
71  TIter next(fileElements);
72  TChainElement *chEl = 0;
73  while ((chEl = (TChainElement*) next()))
74  {
75  flst << chEl->GetTitle() << endl;
76  }
77  flst.close();
78 
79  cout << "Saved file list to " << infile + ".lst" << endl;
80  }
81 
82  assert(_file0);
83 
84  T->SetAlias("UpsilonPair_trk_gpt",
85  "1*sqrt(DST.UpsilonPair.trk.gpx**2 + DST.UpsilonPair.trk.gpy**2)");
86  T->SetAlias("UpsilonPair_trk_pt",
87  "1*sqrt(DST.UpsilonPair.trk.px**2 + DST.UpsilonPair.trk.py**2)");
88 
89  T->SetAlias("EMCalTrk_pt", "1*sqrt(DST.EMCalTrk.px**2 + DST.EMCalTrk.py**2)");
90  T->SetAlias("EMCalTrk_gpt",
91  "1*sqrt(DST.EMCalTrk.gpx**2 + DST.EMCalTrk.gpy**2)");
92 
93  TCut event_sel = "1*1";
94  cuts = "_all_event";
95  if (eval_mode == 0)
96  {
97  event_sel = "Entry$<50000 && fabs(EMCalTrk_pt/EMCalTrk_gpt - 1)<0.05";
98  cuts = "_ll_sample";
99  }
100  else if (eval_mode == 1)
101  {
102 // event_sel = "Entry$>50000 && fabs(EMCalTrk_pt/EMCalTrk_gpt - 1)<0.05";
103  event_sel = " fabs(EMCalTrk_pt/EMCalTrk_gpt - 1)<0.05";
104  cuts = "_eval_sample";
105  }
106  cout << "Build event selection of " << (const char *) event_sel << endl;
107 
108  T->Draw(">>EventList", event_sel);
109  TEventList * elist = gDirectory->GetObjectChecked("EventList", "TEventList");
110  cout << elist->GetN() << " / " << T->GetEntriesFast() << " events selected"
111  << endl;
112 
113  T->SetEventList(elist);
114 
115  if (eval_mode == 0)
116  {
117  Edep_Distribution(infile);
118  }
119  else if (eval_mode == 1)
120  {
121  Edep_LL_Distribution(infile);
122  EP_LL_Distribution(infile);
123  }
124 
125 // Edep_Checks(infile, "fabs(EMCalTrk_pt/EMCalTrk_gpt - 1)<0.05");
126 // Ep_Checks(infile, "fabs(EMCalTrk_pt/EMCalTrk_gpt - 1)<0.05");
127  Edep_Checks(infile, "1");
128  Ep_Checks(infile, "1");
129  // ShowerShape_Checks(infile, "fabs(EMCalTrk_pt/EMCalTrk_gpt - 1)<0.05 && DST.EMCalTrk.cemc_sum_energy>3");
130 }
131 
132 void
133 EP_LL_Distribution(TString infile)
134 {
135 
136  TCanvas *c1 = new TCanvas("EP_LL_Distribution" + cuts,
137  "EP_LL_Distribution" + cuts, 1900, 900);
138  c1->Divide(2, 2);
139  int idx = 1;
140  TPad * p;
141 
142  p = (TPad *) c1->cd(idx++);
143  c1->Update();
144  p->SetLogy();
145  T->Draw("DST.EMCalTrk.ll_ep_e>>hll_ep_e(300,-30,30)", "", "");
146  hll_ep_e->SetTitle(
147  Form(
148  "EMCal only: Electron likelihood distribution;log(Electron likelihood);A. U."));
149 
150  p = (TPad *) c1->cd(idx++);
151  c1->Update();
152  p->SetLogy();
153  T->Draw("DST.EMCalTrk.ll_ep_h>>hll_ep_h(300,-30,30)", "", "");
154  hll_ep_h->SetTitle(
155  Form(
156  "EMCal only: Hadron likelihood distribution;log(Hadron likelihood);A. U."));
157 
158  p = (TPad *) c1->cd(idx++);
159  c1->Update();
160  p->SetLogy();
161  T->Draw(
162  "DST.EMCalTrk.ll_ep_e - DST.EMCalTrk.ll_ep_h>>hll_ep_diff(300,-30,30)",
163  "", "");
164  TH1F *hll_ep_diff = (TH1F *) gDirectory->GetObjectChecked("hll_ep_diff",
165  "TH1F");
166 
167  hll_ep_diff->SetTitle(
168  Form(
169  "EMCal only: log likelihood cut;log(Electron likelihood) - log(Hadron likelihood);Count / bin"));
170 // hll_ep_diff->Scale(1./hll_ep_diff->GetSum());
171 
172  p = (TPad *) c1->cd(idx++);
173  c1->Update();
174 // p->SetLogy();
175 
176  p->DrawFrame(-30, 1e-4, 30, 1,
177  "EMCal only: Cut Efficiency;Cut on log(Electron likelihood) - log(Hadron likelihood);Cut Efficiency");
178  TGraphErrors * ge = Distribution2Efficiency(hll_ep_diff);
179  ge->SetLineColor(kBlue + 2);
180  ge->SetMarkerColor(kBlue + 21);
181  ge->SetMarkerColor(kFullCircle);
182  ge->SetLineWidth(3);
183  ge->Draw("lp");
184 
185  SaveCanvas(c1,
186  TString(_file0->GetName()) + TString("_DrawEcal_Likelihood_")
187  + TString(c1->GetName()), kFALSE);
188 
189 }
190 
191 void
192 Edep_LL_Distribution(TString infile)
193 {
194 
195  TCanvas *c1 = new TCanvas("Edep_LL_Distribution" + cuts,
196  "Edep_LL_Distribution" + cuts, 1900, 900);
197  c1->Divide(2, 2);
198  int idx = 1;
199  TPad * p;
200 
201  p = (TPad *) c1->cd(idx++);
202  c1->Update();
203  p->SetLogy();
204  T->Draw("DST.EMCalTrk.ll_edep_e>>hll_edep_e(300,-30,30)", "", "");
205  hll_edep_e->SetTitle(
206  Form(
207  "EMCal + HCal_{in}: Electron likelihood distribution;log(Electron likelihood);A. U."));
208 
209  p = (TPad *) c1->cd(idx++);
210  c1->Update();
211  p->SetLogy();
212  T->Draw("DST.EMCalTrk.ll_edep_h>>hll_edep_h(300,-30,30)", "", "");
213  hll_edep_h->SetTitle(
214  Form(
215  "EMCal + HCal_{in}: Hadron likelihood distribution;log(Hadron likelihood);A. U."));
216 
217  p = (TPad *) c1->cd(idx++);
218  c1->Update();
219  p->SetLogy();
220  T->Draw(
221  "DST.EMCalTrk.ll_edep_e - DST.EMCalTrk.ll_edep_h>>hll_edep_diff(300,-30,30)",
222  "", "");
223  TH1F *hll_edep_diff = (TH1F *) gDirectory->GetObjectChecked("hll_edep_diff",
224  "TH1F");
225 
226  hll_edep_diff->SetTitle(
227  Form(
228  "EMCal + HCal_{in}: log likelihood cut;log(Electron likelihood) - log(Hadron likelihood);Count / bin"));
229 // hll_edep_diff->Scale(1./hll_edep_diff->GetSum());
230 
231  p = (TPad *) c1->cd(idx++);
232  c1->Update();
233 // p->SetLogy();
234 
235  p->DrawFrame(-30, 1e-4, 30, 1,
236  "EMCal + HCal_{in}: Cut Efficiency;log(Electron likelihood) - log(Hadron likelihood);Cut Efficiency");
237  TGraphErrors * ge = Distribution2Efficiency(hll_edep_diff);
238  ge->SetLineColor(kBlue + 2);
239  ge->SetMarkerColor(kBlue + 21);
240  ge->SetMarkerColor(kFullCircle);
241  ge->SetLineWidth(3);
242  ge->Draw("lp");
243 
244  SaveCanvas(c1,
245  TString(_file0->GetName()) + TString("_DrawEcal_Likelihood_")
246  + TString(c1->GetName()), kFALSE);
247 
248 }
249 
250 void
251 Edep_Distribution(TString infile)
252 {
253 
254  double N_Event = T->GetEntries();
255 
256  TCanvas *c1 = new TCanvas("Edep_Distribution" + cuts,
257  "Edep_Distribution" + cuts, 1900, 900);
258  c1->Divide(2, 1);
259  int idx = 1;
260  TPad * p;
261 
262  p = (TPad *) c1->cd(idx++);
263  c1->Update();
264  p->SetLogz();
265  T->Draw(
266  "DST.EMCalTrk.hcalin_sum_energy:DST.EMCalTrk.get_ep()>>h2_Edep_Distribution_raw(240,-.0,2, 240,-.0,12)",
267  "", "colz");
268  h2_Edep_Distribution_raw->SetTitle(
269  Form(
270  "Energy distribution;CEMC Cluster Energy in GeV;HCal_{in} Cluster Energy in GeV"));
271  h2_Edep_Distribution_raw->Scale(1. / N_Event);
272 // h2_Edep_Distribution_raw->GetZaxis()->SetRangeUser(1. / N_Event, 1);
273 
274  p = (TPad *) c1->cd(idx++);
275  c1->Update();
276  p->SetLogz();
277 
278  TH2F * h2_Edep_Distribution = (TH2F *) h2_Edep_Distribution_raw->Clone(
279  "h2_Edep_Distribution");
280 
281  h2_Edep_Distribution->Smooth(1, "k5b");
282  h2_Edep_Distribution->Scale(1. / h2_Edep_Distribution->GetSum());
283  h2_Edep_Distribution->Draw("colz");
284 
285  SaveCanvas(c1,
286  TString(_file0->GetName()) + TString("_DrawEcal_Likelihood_")
287  + TString(c1->GetName()), kFALSE);
288 
289 }
290 
291 void
292 ShowerShape_Checks(TString infile, TCut good_track_cut)
293 {
294 
295  double N_Event = T->GetEntries();
296 
297  TCanvas *c1 = new TCanvas("ShowerShape_Checks" + cuts,
298  "ShowerShape_Checks" + cuts, 1900, 900);
299  c1->Divide(2, 1);
300  int idx = 1;
301  TPad * p;
302 
303  p = (TPad *) c1->cd(idx++);
304  c1->Update();
305  p->SetLogz();
306 
307  T->Draw(
308  "DST.EMCalTrk.cemc_energy/DST.EMCalTrk.cemc_sum_energy:DST.EMCalTrk.cemc_radius>>hEMCalTrk_cemc_shape(30,0,2,100,0,1)",
309  good_track_cut, "goff");
310 
311  TH2 * hEMCalTrk_cemc_shape = (TH2 *) gDirectory->GetObjectChecked(
312  "hEMCalTrk_cemc_shape", "TH2");
313  hEMCalTrk_cemc_shape->SetTitle(
314  Form(
315  "CEMC Shower Shape;Distance to track projection (Cluster width);Tower Energy/Cluster Energy"));
316 // hEMCalTrk_cemc_shape->Scale(1. / N_Event);
317 
318  TH1D * hEMCalTrk_cemc_shape_px = hEMCalTrk_cemc_shape->ProjectionX(
319  "hEMCalTrk_cemc_shape_px");
320 
321  for (int r = 1; r <= hEMCalTrk_cemc_shape->GetNbinsX(); r++)
322  for (int e = 1; e <= hEMCalTrk_cemc_shape->GetNbinsY(); e++)
323  {
324  hEMCalTrk_cemc_shape->SetBinContent(r, e,
325  hEMCalTrk_cemc_shape->GetBinContent(r, e)
326  / hEMCalTrk_cemc_shape_px->GetBinContent(r));
327  }
328 
329  hEMCalTrk_cemc_shape->Draw("colz");
330 
331  // return;
332  p = (TPad *) c1->cd(idx++);
333  c1->Update();
334  p->SetLogz();
335 
336  T->Draw(
337  "DST.EMCalTrk.hcalin_energy/DST.EMCalTrk.hcalin_sum_energy:DST.EMCalTrk.hcalin_radius>>hEMCalTrk_hcalin_shape(30,0,2,100,0,1)",
338  good_track_cut, "colz");
339  TH2 * hEMCalTrk_hcalin_shape = (TH2 *) gDirectory->GetObjectChecked(
340  "hEMCalTrk_hcalin_shape", "TH2");
341 
342  hEMCalTrk_hcalin_shape->SetTitle(
343  Form(
344  "HCal_{in} Shower Shape;Distance to track projection (Cluster width);Tower Energy/Cluster Energy"));
345 // hEMCalTrk_hcalin_shape->Scale(1. / N_Event);
346 
347  TH1D * hEMCalTrk_hcalin_shape_px = hEMCalTrk_hcalin_shape->ProjectionX(
348  "hEMCalTrk_hcalin_shape_px");
349 
350  for (int r = 1; r <= hEMCalTrk_hcalin_shape->GetNbinsX(); r++)
351  for (int e = 1; e <= hEMCalTrk_hcalin_shape->GetNbinsY(); e++)
352  {
353  hEMCalTrk_hcalin_shape->SetBinContent(r, e,
354  hEMCalTrk_hcalin_shape->GetBinContent(r, e)
355  / hEMCalTrk_hcalin_shape_px->GetBinContent(r));
356  }
357 
358  hEMCalTrk_hcalin_shape->Draw("colz");
359 
360  SaveCanvas(c1,
361  TString(_file0->GetName()) + TString("_DrawEcal_Likelihood_")
362  + TString(c1->GetName()), kFALSE);
363 
364 }
365 
366 void
367 Edep_Checks(TString infile, TCut good_track_cut)
368 {
369 
370  double N_Event = T->GetEntries();
371 
372  TCanvas *c1 = new TCanvas("Edep_Checks" + cuts, "Edep_Checks" + cuts, 1900,
373  950);
374  c1->Divide(3, 2);
375  int idx = 1;
376  TPad * p;
377 
378  p = (TPad *) c1->cd(idx++);
379  c1->Update();
380  p->SetLogy();
381  T->Draw("DST.EMCalTrk.cemc_sum_n_tower>>hEMCalTrk_cemc_ntower(16,-.5,15.5)",
382  good_track_cut);
383  hEMCalTrk_cemc_ntower->SetTitle(
384  Form("CEMC Cluster Size;Cluster Size (Towers);Probability"));
385  hEMCalTrk_cemc_ntower->Scale(1. / N_Event);
386 
387  p = (TPad *) c1->cd(idx++);
388  c1->Update();
389  p->SetLogy();
390  T->Draw("DST.EMCalTrk.cemc_sum_energy>>hEMCalTrk_cemc_e(240,-.0,12)",
391  good_track_cut);
392  hEMCalTrk_cemc_e->SetTitle(
393  Form("CEMC Cluster Energy;Cluster Energy (GeV);Count/bin"));
394 // hEMCalTrk_cemc_e->Scale(1. / N_Event);
395 
396  p = (TPad *) c1->cd(idx++);
397  c1->Update();
398  p->SetLogy();
399  p->DrawFrame(-.0, 1e-3, 12, 1,
400  "CEMC Cut Eff;Cut on Cluster Energy (GeV);Efficiency");
401 
402  TGraphErrors * ge = Distribution2Efficiency(hEMCalTrk_cemc_e);
403  ge->SetLineColor(kBlue + 2);
404  ge->SetMarkerColor(kBlue + 21);
405  ge->SetMarkerColor(kFullCircle);
406  ge->SetLineWidth(3);
407  ge->Draw("lp");
408 
409  p = (TPad *) c1->cd(idx++);
410  c1->Update();
411  p->SetLogy();
412  T->Draw(
413  "DST.EMCalTrk.hcalin_sum_n_tower>>hEMCalTrk_hcalin_ntower(16,-.5,15.5)",
414  good_track_cut);
415  hEMCalTrk_hcalin_ntower->SetTitle(
416  Form("HCal_{in} Cluster Size;Cluster Size (Towers);Probability"));
417  hEMCalTrk_hcalin_ntower->Scale(1. / N_Event);
418 
419  p = (TPad *) c1->cd(idx++);
420  c1->Update();
421  p->SetLogy();
422  T->Draw("DST.EMCalTrk.hcalin_sum_energy>>hEMCalTrk_hcalin_e(240,-.0,12)",
423  good_track_cut);
424  hEMCalTrk_hcalin_e->SetTitle(
425  Form("HCal_{in} Cluster Energy;Cluster Energy (GeV);Count/bin"));
426 // hEMCalTrk_hcalin_e->Scale(1. / N_Event);
427 
428  p = (TPad *) c1->cd(idx++);
429  c1->Update();
430  p->SetLogy();
431  p->DrawFrame(-.0, 1e-3, 12, 1,
432  "HCal_{in} Cut Eff;Cut on Cluster Energy (GeV);Efficiency");
433 
434  TGraphErrors * ge = Distribution2Efficiency(hEMCalTrk_hcalin_e);
435  ge->SetLineColor(kBlue + 2);
436  ge->SetMarkerColor(kBlue + 21);
437  ge->SetMarkerColor(kFullCircle);
438  ge->SetLineWidth(3);
439  ge->Draw("lp");
440 
441  SaveCanvas(c1,
442  TString(_file0->GetName()) + TString("_DrawEcal_pDST_")
443  + TString(c1->GetName()), kFALSE);
444 
445  TCanvas *c1 = new TCanvas("Edep_Checks_2D" + cuts, "Edep_Checks_2D" + cuts,
446  900, 900);
447 // c1->Divide(2, 2);
448 // int idx = 1;
449 // TPad * p;
450 
451  p = (TPad *) c1->cd(idx++);
452  c1->Update();
453  p->SetLogz();
454  T->Draw(
455  "DST.EMCalTrk.hcalin_sum_energy:DST.EMCalTrk.cemc_sum_energy>>h2_EMCalTrk_hcalin_e_EMCalTrk_cemc_e(240,-.0,8, 240,-.0,8)",
456  good_track_cut, "colz");
457  h2_EMCalTrk_hcalin_e_EMCalTrk_cemc_e->SetTitle(
458  Form(
459  "Energy distribution;CEMC Cluster Energy in GeV;HCal_{in} Cluster Energy in GeV"));
460  h2_EMCalTrk_hcalin_e_EMCalTrk_cemc_e->Scale(1. / N_Event);
461  h2_EMCalTrk_hcalin_e_EMCalTrk_cemc_e->GetZaxis()->SetRangeUser(1. / N_Event,
462  1);
463 
464  SaveCanvas(c1,
465  TString(_file0->GetName()) + TString("_DrawEcal_Likelihood_")
466  + TString(c1->GetName()), kFALSE);
467 
468 }
469 
470 void
471 Ep_Checks(TString infile, TCut good_track_cut)
472 {
473 
474  double N_Event = T->GetEntries();
475 
476  TCanvas *c1 = new TCanvas("Ep_Checks" + cuts, "Ep_Checks" + cuts, 1900, 950);
477  c1->Divide(2, 1);
478  int idx = 1;
479  TPad * p;
480 
481  p = (TPad *) c1->cd(idx++);
482  c1->Update();
483  p->SetLogy();
484  T->Draw("DST.EMCalTrk.get_ep()>>hEMCalTrk_get_ep(240,-.0,2)",
485  good_track_cut);
486  hEMCalTrk_get_ep->SetTitle(
487  Form("CEMC Cluster Energy/Track Momentum;E/p;Count/bin"));
488 // hEMCalTrk_cemc_e->Scale(1. / N_Event);
489 
490  p = (TPad *) c1->cd(idx++);
491  c1->Update();
492 
493  if (infile.Contains("e-") || infile.Contains("e+"))
494  {
495  p->DrawFrame(-.0, 0.8, 1.5, 1,
496  "CEMC E/p Cut Eff;Cut on E/p;Signal Efficiency");
497  }
498  else
499  {
500  p->DrawFrame(-.0, 1e-3, 1.5, 1,
501  "CEMC E/p Cut Eff;Cut on E/p;Background Efficiency or 1/Rejection");
502  p->SetLogy();
503  }
504  TGraphErrors * ge = Distribution2Efficiency(hEMCalTrk_get_ep);
505  ge->SetLineColor(kBlue + 2);
506  ge->SetMarkerColor(kBlue + 21);
507  ge->SetMarkerColor(kFullCircle);
508  ge->SetLineWidth(3);
509  ge->Draw("lp");
510 
511  SaveCanvas(c1,
512  TString(_file0->GetName()) + TString("_DrawEcal_Likelihood_")
513  + TString(c1->GetName()), kFALSE);
514 
515 }
516 
517 TGraphErrors *
518 Distribution2Efficiency(TH1F * hCEMC3_Max)
519 {
520  double threshold[10000] =
521  { 0 };
522  double eff[10000] =
523  { 0 };
524  double eff_err[10000] =
525  { 0 };
526 
527  assert(hCEMC3_Max->GetNbinsX()<10000);
528 
529  const double n = hCEMC3_Max->GetSum();
530  double pass = 0;
531  int cnt = 0;
532  for (int i = hCEMC3_Max->GetNbinsX(); i >= 1; i--)
533  {
534  pass += hCEMC3_Max->GetBinContent(i);
535 
536  const double pp = pass / n;
537 // const double z = 1.96;
538  const double z = 1.;
539 
540  const double A = z * sqrt(1. / n * pp * (1 - pp) + z * z / 4 / n / n);
541  const double B = 1 / (1 + z * z / n);
542 
543  threshold[cnt] = hCEMC3_Max->GetBinCenter(i);
544  eff[cnt] = (pp + z * z / 2 / n) * B;
545  eff_err[cnt] = A * B;
546 
547 // cout << threshold[cnt] << ": " << "CL " << eff[cnt] << "+/-"
548 // << eff_err[cnt] << endl;
549  cnt++;
550  }
551  TGraphErrors * ge = new TGraphErrors(cnt, threshold, eff, NULL, eff_err);
552 
553  return ge;
554 }