EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LiteCaloEval.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LiteCaloEval.cc
1 #include "LiteCaloEval.h"
2 
3 #include <calobase/RawTower.h>
4 #include <calobase/RawTowerContainer.h>
5 #include <calobase/RawTowerGeom.h>
6 #include <calobase/RawTowerGeomContainer.h>
7 
9 #include <fun4all/SubsysReco.h>
10 
11 #include <phool/PHCompositeNode.h>
12 #include <phool/getClass.h>
13 #include <phool/phool.h>
14 
15 #include <TF1.h>
16 #include <TFile.h>
17 #include <TGraph.h>
18 #include <TH1.h>
19 #include <TH2.h>
20 #include <TH3.h>
21 #include <TStyle.h>
22 #include <TSystem.h>
23 
24 #include <cstdlib>
25 #include <iostream>
26 #include <cmath> // for abs
27 #include <map> // for _Rb_tree_const_iterator
28 #include <memory>
29 #include <utility> // for pair
30 
31 //____________________________________________________________________________..
32 LiteCaloEval::LiteCaloEval(const std::string& name, const std::string& caloname, const std::string& filename)
33  : SubsysReco(name)
34  , _caloname(caloname)
35  , _filename(filename)
36 {
37 }
38 
39 //____________________________________________________________________________..
41 {
42  // just quit if we forgot to set the calorimeter type
44  {
45  std::cout << Name() << ": No Calo Type set" << std::endl;
46  gSystem->Exit(1);
47  }
48 
49  _ievent = 0;
50 
51  cal_output = new TFile(_filename.c_str(), "RECREATE");
52 
54  {
55  hcalin_energy_eta = new TH2F("hcalin_energy_eta", "hcalin energy eta", 1000, 0, 10, 240, -1.1, 1.1);
56  hcalin_e_eta_phi = new TH3F("hcalin_e_eta_phi", "hcalin e eta phi", 50, 0, 10, 24, -1.1, 1.1, 64, -3.14159, 3.14159);
57  for (int i = 0; i < 24; i++)
58  {
59  for (int j = 0; j < 64; j++)
60  {
61  std::string hist_name = "hcal_in_eta_" + std::to_string(i) + "_phi_" + std::to_string(j);
62 
63  hcal_in_eta_phi[i][j] = new TH1F(hist_name.c_str(), "Hcal_in_energy", 1000, 0, 10);
64  }
65  }
66 
67  for (int i = 0; i < 24; i++)
68  {
69  std::string hist_name = "hcalin_eta_" + std::to_string(i);
70 
71  hcalin_eta[i] = new TH1F(hist_name.c_str(), "hcalin eta's", 1000, 0, 10);
72  }
73  }
74  else if (calotype == LiteCaloEval::HCALOUT)
75  {
76  hcalout_energy_eta = new TH2F("hcalout_energy_eta", "hcalout energy eta", 10, 0, 10, 24000, -1.1, 1.1);
77  hcalout_e_eta_phi = new TH3F("hcalout_e_eta_phi", "hcalout e eta phi", 50, 0, 10, 24, -1.1, 1.1, 64, -3.14159, 3.14159);
78  for (int i = 0; i < 24; i++)
79  {
80  for (int j = 0; j < 64; j++)
81  {
82  std::string hist_name = "hcal_out_eta_" + std::to_string(i) + "_phi_" + std::to_string(j);
83 
84  hcal_out_eta_phi[i][j] = new TH1F(hist_name.c_str(), "Hcal_out energy", 1000, 0, 10);
85  }
86  }
87 
88  for (int i = 0; i < 24; i++)
89  {
90  std::string hist_name = "hcalout_eta_" + std::to_string(i);
91 
92  hcalout_eta[i] = new TH1F(hist_name.c_str(), "hcalout eta's", 1000, 0, 10);
93  }
94  }
95  else if (calotype == LiteCaloEval::CEMC)
96  {
97  for (int i = 0; i < 96; i++)
98  {
99  for (int j = 0; j < 258; j++)
100  {
101  std::string hist_name = "emc_ieta" + std::to_string(i) + "_phi" + std::to_string(j);
102 
103  cemc_hist_eta_phi[i][j] = new TH1F(hist_name.c_str(), "Hist_ieta_phi_leaf(e)", 1000, 0, 10);
104  }
105  }
106 
107  //create the eta 1d histos
108  for (int i = 0; i < 96; i++)
109  {
110  gStyle->SetOptFit(1);
111  std::string b = "eta_" + std::to_string(i);
112 
113  eta_hist[i] = new TH1F(b.c_str(), "eta and all phi's", 1000, 0, 10);
114  }
115 
116  //make 2d histo
117  energy_eta_hist = new TH2F("energy_eta_hist", "energy eta and all phi", 10, 0, 10, 9600, -1, 1);
118 
119  //make 3d histo
120  e_eta_phi = new TH3F("e_eta_phi", "e v eta v phi", 50, 0, 10, 100, -1, 1, 256, -3.14159, 3.14159);
121  }
122 
124 }
125 
126 //____________________________________________________________________________..
128 {
129  if (_ievent % 100 == 0) std::cout << "LiteCaloEval::process_event(PHCompositeNode *topNode) Processing Event " << _ievent << std::endl;
130 
131  std::string towernode = "TOWER_CALIB_" + _caloname;
132  RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(topNode, towernode);
133  if (!towers)
134  {
135  std::cout << PHWHERE << " ERROR: Can't find " << towernode << std::endl;
136  exit(-1);
137  }
138 
139  std::string towergeomnode = "TOWERGEOM_" + _caloname;
140  RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnode);
141  if (!towergeom)
142  {
143  std::cout << PHWHERE << " ERROR: Can't find " << towergeomnode << std::endl;
144  exit(-1);
145  }
146 
147  RawTowerContainer::ConstRange begin_end = towers->getTowers();
149  for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
150  {
151  RawTower* tower = rtiter->second;
152 
153  if (tower->get_energy() < 0.0005) continue;
154 
155  RawTowerGeom* tower_geom = towergeom->get_tower_geometry(tower->get_id());
156  if (!tower_geom)
157  {
158  std::cout << PHWHERE << " ERROR: Can't find tower geometry for this tower hit: ";
159  tower->identify();
160  exit(-1);
161  }
162 
163  const int towerid = tower->get_id();
164  const int ieta = tower->get_bineta();
165  const int iphi = tower->get_binphi();
166  const float eta = tower_geom->get_eta();
167  const float phi = tower_geom->get_phi();
168  const float e = tower->get_energy();
170  {
171  if (towerid < 0)
172  {
173  std::cout << "a towerid was less than 0 " << std::endl;
174  break;
175  }
176 
177  //fill the hist with energy data
178  cemc_hist_eta_phi[ieta][iphi]->Fill(e);
179 
180  //fill and fit the 1d hist eta and all phi
181  eta_hist[ieta]->Fill(e);
182 
183  //fill the 2d histo eta, energy and all phi
184  energy_eta_hist->Fill(e, eta);
185 
186  //fill 3d histo e_eta_phid
187  e_eta_phi->Fill(e, eta, phi);
188  }
189  else if (calotype == LiteCaloEval::HCALOUT)
190  {
191  //fill the hist with energy data
192  //std::cout << ieta << " " << iphi << std::endl;
193 
194  hcal_out_eta_phi[ieta][iphi]->Fill(e);
195 
196  hcalout_eta[ieta]->Fill(e);
197 
198  hcalout_energy_eta->Fill(e, eta);
199 
200  hcalout_e_eta_phi->Fill(e, eta, phi);
201  }
202  else if (calotype == LiteCaloEval::HCALIN)
203  {
204  //fill the hist with energy data
205 
206  hcal_in_eta_phi[ieta][iphi]->Fill(e);
207 
208  hcalin_eta[ieta]->Fill(e);
209 
210  hcalin_energy_eta->Fill(e, eta);
211 
212  hcalin_e_eta_phi->Fill(e, eta, phi);
213  }
214  }
215 
216  _ievent++;
217 
219 }
220 
221 //____________________________________________________________________________..
223 {
224  cal_output->cd();
226  {
227  double par_value[24];
228  double eta_value[24];
229  double par_err[24];
230  double rel_err[24];
231 
232  //same as above but for the second fit performed below in code
233  double par_value2[24];
234  double par_err2[24];
235  double rel_err2[24];
236 
237  for (int i = 0; i < 24; i++)
238  {
239  //a,b,c,d hold the par value, error, par2 value, par2 error respectivley. The 2 refers to the second fit below
240  double a;
241  double b;
242  double c;
243  double d;
244 
245  //make functions to fit the eta histos
246  std::unique_ptr<TF1> f1(new TF1("f1", "expo", 0.02, 0.1));
247  std::unique_ptr<TF1> f2(new TF1("f2", "expo", 0.1, 1));
248  f2->SetLineColor(7);
249 
250  hcalin_eta[i]->Fit("f1", "R");
251  hcalin_eta[i]->Fit("f2", "R+");
252 
253  a = f1->GetParameter(1);
254  b = f1->GetParError(1);
255 
256  c = f2->GetParameter(1);
257  d = f2->GetParError(1);
258 
259  //assign retreived parameter values
260  par_value[i] = abs(a);
261  par_err[i] = b;
262  eta_value[i] = i;
263  rel_err[i] = par_err[i] / par_value[i];
264 
265  par_value2[i] = abs(c);
266  par_err2[i] = d;
267  rel_err2[i] = par_err2[i] / par_value2[i];
268  }
269 
270  //create graphs
271  TGraph g1(24, eta_value, par_value);
272  g1.SetTitle("HCal In (0.02-0.1 GeV); eta; p1");
273  g1.SetMarkerStyle(20);
274  g1.Draw("ap");
275  g1.SetName("Fit1_hcalin");
276  g1.Write();
277 
278  TGraph g2(24, eta_value, rel_err);
279  g2.SetTitle("HCal In Error (0.02-0.1 GeV); eta; p1 rel error");
280  g2.SetMarkerStyle(20);
281  g2.Draw("ap");
282  g2.SetName("Fit1_err_hcalin");
283  g2.Write();
284 
285  TGraph g3(24, eta_value, par_value2);
286  g3.SetTitle("HCal In (0.1-1 GeV); eta; p1");
287  g3.SetMarkerStyle(20);
288  g3.Draw("ap");
289  g3.SetName("Fit2_hcalin");
290  g3.Write();
291 
292  TGraph g4(24, eta_value, rel_err2);
293  g4.SetTitle("HCal In Error (0.1-1 GeV); eta; p1 rel error");
294  g4.SetMarkerStyle(20);
295  g4.Draw("ap");
296  g4.SetName("Fit2_err_hcalin");
297  g4.Write();
298  }
299  else if (calotype == LiteCaloEval::HCALOUT)
300  {
301  double par_value[24];
302  double eta_value[24];
303  double par_err[24];
304  double rel_err[24];
305 
306  //same as above but for the second/third fit performed below in code
307  double par_value2[24];
308  double par_err2[24];
309  double rel_err2[24];
310 
311  double par_value3[24];
312  double par_err3[24];
313  double rel_err3[24];
314 
315  for (int i = 0; i < 24; i++)
316  {
317  //a,b hold the p1 value, p1 error, and then repeats for the second fit (c,d) and third fit (e,f)
318  double a;
319  double b;
320  double c;
321  double d;
322  double e;
323  double f;
324 
325  //make functions to fit the eta histos
326  std::unique_ptr<TF1> f1(new TF1("f1", "expo", 0.05, 0.2));
327  std::unique_ptr<TF1> f2(new TF1("f2", "expo", 0.2, 1));
328  std::unique_ptr<TF1> f3(new TF1("f3", "expo", 1, 2));
329  f2->SetLineColor(7);
330  f3->SetLineColor(1);
331 
332  hcalout_eta[i]->Fit("f1", "R");
333  hcalout_eta[i]->Fit("f2", "R+");
334  hcalout_eta[i]->Fit("f3", "R+");
335 
336  a = f1->GetParameter(1);
337  b = f1->GetParError(1);
338 
339  c = f2->GetParameter(1);
340  d = f2->GetParError(1);
341 
342  e = f3->GetParameter(1);
343  f = f3->GetParError(1);
344 
345  //assign retreived parameter values
346  par_value[i] = abs(a);
347  par_err[i] = b;
348  eta_value[i] = i;
349  rel_err[i] = par_err[i] / par_value[i];
350 
351  par_value2[i] = abs(c);
352  par_err2[i] = d;
353  rel_err2[i] = par_err2[i] / par_value2[i];
354 
355  par_value3[i] = abs(e);
356  par_err3[i] = f;
357  rel_err3[i] = par_err3[i] / par_value3[i];
358  }
359 
360  //create graphs
361  TGraph g1(24, eta_value, par_value);
362  g1.SetTitle("HCal Out (0.05-0.2 GeV); eta; p1");
363  g1.SetMarkerStyle(20);
364  g1.Draw("ap");
365  g1.SetName("Fit1_hcalout");
366  g1.Write();
367 
368  TGraph g2(24, eta_value, rel_err);
369  g2.SetTitle("HCal Out Error (0.05-0.2 GeV); eta; p1 rel err");
370  g2.SetMarkerStyle(20);
371  g2.Draw("ap");
372  g2.SetName("Fit1_err_hcalout");
373  g2.Write();
374 
375  TGraph g3(24, eta_value, par_value2);
376  g3.SetTitle("HCal Out (0.2-1 GeV); eta; p1");
377  g3.SetMarkerStyle(20);
378  g3.Draw("ap");
379  g3.SetName("Fit2_hcalout");
380  g3.Write();
381 
382  TGraph g4(24, eta_value, rel_err2);
383  g4.SetTitle("HCal Out Error (0.2-1 GeV); eta; p1 rel err");
384  g4.SetMarkerStyle(20);
385  g4.Draw("ap");
386  g4.SetName("Fit2_err_hcalout");
387  g4.Write();
388 
389  TGraph g5(24, eta_value, par_value3);
390  g5.SetTitle("HCal Out (1-2 GeV); eta; p1");
391  g5.SetMarkerStyle(20);
392  g5.Draw("ap");
393  g5.SetName("Fit3_hcalout");
394  g5.Write();
395 
396  TGraph g6(24, eta_value, rel_err3);
397  g6.SetTitle("HCal Out Error (1-2 GeV); eta; p1 rel err");
398  g6.SetMarkerStyle(20);
399  g6.Draw("ap");
400  g6.SetName("Fit3_err_hcalout");
401  g6.Write();
402  }
403  else if (calotype == LiteCaloEval::CEMC)
404  {
405  //create arrays that holds parameter (p1) value and error
406  double par_value[96];
407  double eta_value[96];
408  double par_err[96];
409  double rel_err[96];
410 
411  //same as above but for the second fit performed below in code
412  double par_value2[96];
413  double par_err2[96];
414  double rel_err2[96];
415 
416  double par_value3[96];
417  double par_err3[96];
418  double rel_err3[96];
419 
420  //create graphs for parameter (p1) vs eta
421  for (int i = 0; i < 96; i++)
422  {
423  //a,b,c,d hold the par value, error, par2 value, par2 error respectivley. The 2 refers to the second fit below
424  double a;
425  double b;
426  double c;
427  double d;
428  double e;
429  double f;
430 
431  //make functions to fit the eta histos
432  std::unique_ptr<TF1> f1(new TF1("f1", "expo", 0.04, 0.1));
433  std::unique_ptr<TF1> f2(new TF1("f2", "expo", 0.1, 0.4));
434  std::unique_ptr<TF1> f3(new TF1("f3", "expo", 0.4, 2));
435  f2->SetLineColor(7);
436  f3->SetLineColor(1);
437 
438  eta_hist[i]->Fit("f1", "R");
439  eta_hist[i]->Fit("f2", "R+");
440  eta_hist[i]->Fit("f3", "R+");
441 
442  a = f1->GetParameter(1);
443  b = f1->GetParError(1);
444 
445  c = f2->GetParameter(1);
446  d = f2->GetParError(1);
447 
448  e = f3->GetParameter(1);
449  f = f3->GetParError(1);
450 
451  //assign retreived parameter values
452  par_value[i] = abs(a);
453  par_err[i] = b;
454  eta_value[i] = i;
455  rel_err[i] = par_err[i] / par_value[i];
456 
457  par_value2[i] = abs(c);
458  par_err2[i] = d;
459  rel_err2[i] = par_err2[i] / par_value2[i];
460 
461  par_value3[i] = abs(e);
462  par_err3[i] = f;
463  rel_err3[i] = par_err3[i] / par_value3[i];
464  }
465 
466  //create graphs
467  TGraph g1(96, eta_value, par_value);
468  g1.SetTitle("EMCal (0.04-0.1 GeV); eta; p1");
469  g1.SetMarkerStyle(20);
470  g1.Draw("ap");
471  g1.SetName("Fit1_emc");
472  g1.Write();
473 
474  TGraph g2(96, eta_value, rel_err);
475  g2.SetTitle("EMCal Error (0.04-0.1 GeV); eta; p1 rel error");
476  g2.SetMarkerStyle(20);
477  g2.Draw("ap");
478  g2.SetName("Fit1_err_emc");
479  g2.Write();
480 
481  TGraph g3(96, eta_value, par_value2);
482  g3.SetTitle("EMCal (0.1-0.4 GeV); eta; p1");
483  g3.SetMarkerStyle(20);
484  g3.Draw("ap");
485  g3.SetName("Fit2_emc");
486  g3.Write();
487 
488  TGraph g4(96, eta_value, rel_err2);
489  g4.SetTitle("EMCal Error (0.1-0.4 GeV); eta; p1 rel error");
490  g4.SetMarkerStyle(20);
491  g4.Draw("ap");
492  g4.SetName("Fit2_err_emc");
493  g4.Write();
494 
495  TGraph g5(96, eta_value, par_value3);
496  g5.SetTitle("EMCal (0.4-2 GeV); eta; p1");
497  g5.SetMarkerStyle(20);
498  g5.Draw("ap");
499  g5.SetName("Fit3_emc");
500  g5.Write();
501 
502  TGraph g6(96, eta_value, rel_err3);
503  g6.SetTitle("EMCal Error (0.4-2 GeV); eta; p1 rel err");
504  g6.SetMarkerStyle(20);
505  g6.Draw("ap");
506  g6.SetName("Fit3_err_emc");
507  g6.Write();
508  }
509 
510  cal_output->Write();
511 
513 }