EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QAG4SimulationCalorimeterSum.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QAG4SimulationCalorimeterSum.cc
1 #include "QAG4SimulationCalorimeterSum.h"
2 #include "QAHistManagerDef.h"
3 
4 #include <g4eval/CaloEvalStack.h>
5 #include <g4eval/CaloRawClusterEval.h>
6 #include <g4eval/SvtxEvalStack.h>
7 
8 #include <g4main/PHG4Particle.h>
9 #include <g4main/PHG4TruthInfoContainer.h>
10 
11 #include <calobase/RawCluster.h>
12 #include <calobase/RawTower.h>
13 #include <calobase/RawTowerContainer.h>
14 #include <calobase/RawTowerGeomContainer.h>
15 
16 #include <trackbase_historic/SvtxTrack.h>
17 
18 #include <g4eval/SvtxTrackEval.h> // for SvtxTrackEval
19 
20 #include <fun4all/Fun4AllHistoManager.h>
21 #include <fun4all/Fun4AllReturnCodes.h>
22 #include <fun4all/SubsysReco.h>
23 
24 #include <phool/getClass.h>
25 
26 #include <TAxis.h>
27 #include <TH1.h>
28 #include <TH2.h>
29 #include <TNamed.h>
30 #include <TString.h>
31 
32 #include <cassert>
33 #include <cmath>
34 #include <iostream>
35 #include <iterator> // for reverse_iterator
36 #include <utility> // for pair
37 #include <vector>
38 
39 QAG4SimulationCalorimeterSum::QAG4SimulationCalorimeterSum(
40  QAG4SimulationCalorimeterSum::enu_flags flags)
41  : SubsysReco("QAG4SimulationCalorimeterSum")
42  , _flags(flags)
43  , _calo_name_cemc("CEMC")
44  , _calo_name_hcalin("HCALIN")
45  , _calo_name_hcalout("HCALOUT")
46  , _truth_container(nullptr)
47  , _magField(+1.4)
48 {
49 }
50 
51 int QAG4SimulationCalorimeterSum::InitRun(PHCompositeNode *topNode)
52 {
53  _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode,
54  "G4TruthInfo");
55  if (!_truth_container)
56  {
57  std::cout << "QAG4SimulationCalorimeterSum::InitRun - Fatal Error - "
58  << "unable to find DST node "
59  << "G4TruthInfo" << std::endl;
60  assert(_truth_container);
61  }
62 
63  if (flag(kProcessCluster))
64  {
65  if (!_caloevalstack_cemc)
66  {
67  _caloevalstack_cemc.reset(
68  new CaloEvalStack(topNode, _calo_name_cemc));
69  _caloevalstack_cemc->set_strict(true);
70  _caloevalstack_cemc->set_verbosity(Verbosity() + 1);
71  }
72  if (!_caloevalstack_hcalin)
73  {
74  _caloevalstack_hcalin.reset(
75  new CaloEvalStack(topNode, _calo_name_hcalin));
76  _caloevalstack_hcalin->set_strict(true);
77  _caloevalstack_hcalin->set_verbosity(Verbosity() + 1);
78  }
79  if (!_caloevalstack_hcalout)
80  {
81  _caloevalstack_hcalout.reset(
82  new CaloEvalStack(topNode, _calo_name_hcalout));
83  _caloevalstack_hcalout->set_strict(true);
84  _caloevalstack_hcalout->set_verbosity(Verbosity() + 1);
85  }
86  }
87 
88  if (flag(kProcessTrackProj))
89  {
90  if (!_svtxevalstack)
91  {
92  _svtxevalstack.reset(new SvtxEvalStack(topNode));
93  _svtxevalstack->set_strict(true);
94  _svtxevalstack->set_verbosity(Verbosity() + 1);
95  }
96  }
97  return Fun4AllReturnCodes::EVENT_OK;
98 }
99 
100 int QAG4SimulationCalorimeterSum::Init(PHCompositeNode *topNode)
101 {
102  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
103  assert(hm);
104  TH1D *h = new TH1D(TString(get_histo_prefix()) + "Normalization", //
105  TString(get_histo_prefix()) + " Normalization;Items;Count", 10, .5, 10.5);
106  int i = 1;
107  h->GetXaxis()->SetBinLabel(i++, "Event");
108  h->GetXaxis()->SetBinLabel(i++, (_calo_name_cemc + " Tower").c_str());
109  h->GetXaxis()->SetBinLabel(i++, (_calo_name_hcalin + " Tower").c_str());
110  h->GetXaxis()->SetBinLabel(i++, (_calo_name_hcalout + " Tower").c_str());
111  h->GetXaxis()->SetBinLabel(i++, (_calo_name_cemc + " Cluster").c_str());
112  h->GetXaxis()->SetBinLabel(i++, (_calo_name_hcalin + " Cluster").c_str());
113  h->GetXaxis()->SetBinLabel(i++, (_calo_name_hcalout + " Cluster").c_str());
114  h->GetXaxis()->SetBinLabel(i++, "Track");
115  h->GetXaxis()->LabelsOption("v");
116  hm->registerHisto(h);
117 
118  // if (flag(kProcessTower))
119  // {
120  // if (Verbosity() >= 1)
121  // std::cout << "QAG4SimulationCalorimeterSum::Init - Process tower occupancies"
122  // << std::endl;
123  // Init_Tower(topNode);
124  // }
125  if (flag(kProcessCluster))
126  {
127  if (Verbosity() >= 1)
128  std::cout << "QAG4SimulationCalorimeterSum::Init - Process tower occupancies"
129  << std::endl;
130  Init_Cluster(topNode);
131  }
132 
133  if (flag(kProcessTrackProj))
134  {
135  if (Verbosity() >= 1)
136  std::cout << "QAG4SimulationCalorimeterSum::Init - Process sampling fraction"
137  << std::endl;
138  Init_TrackProj(topNode);
139  }
140  return Fun4AllReturnCodes::EVENT_OK;
141 }
142 
143 int QAG4SimulationCalorimeterSum::process_event(PHCompositeNode *topNode)
144 {
145  if (Verbosity() > 2)
146  std::cout << "QAG4SimulationCalorimeterSum::process_event() entered" << std::endl;
147 
148  if (_caloevalstack_cemc)
149  _caloevalstack_cemc->next_event(topNode);
150  if (_caloevalstack_hcalin)
151  _caloevalstack_hcalin->next_event(topNode);
152  if (_caloevalstack_hcalout)
153  _caloevalstack_hcalout->next_event(topNode);
154  if (_svtxevalstack)
155  _svtxevalstack->next_event(topNode);
156 
157  // if (flag(kProcessTower))
158  // {
159  // int ret = process_event_Tower(topNode);
160  //
161  // if (ret != Fun4AllReturnCodes::EVENT_OK)
162  // return ret;
163  // }
164 
165  if (flag(kProcessCluster))
166  {
167  int ret = process_event_Cluster(topNode);
168 
169  if (ret != Fun4AllReturnCodes::EVENT_OK)
170  return ret;
171  }
172 
173  if (flag(kProcessTrackProj))
174  {
175  int ret = process_event_TrackProj(topNode);
176 
177  if (ret != Fun4AllReturnCodes::EVENT_OK)
178  return ret;
179  }
180 
181  // at the end, count success events
182  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
183  assert(hm);
184  TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(
185  get_histo_prefix() + "Normalization"));
186  assert(h_norm);
187  h_norm->Fill("Event", 1);
188 
189  return Fun4AllReturnCodes::EVENT_OK;
190 }
191 
192 std::string
193 QAG4SimulationCalorimeterSum::get_histo_prefix()
194 {
195  return "h_QAG4Sim_CalorimeterSum_";
196 }
197 
198 PHG4Particle *
199 QAG4SimulationCalorimeterSum::get_truth_particle()
200 {
201  // get last primary
202  assert(_truth_container);
203  assert(not _truth_container->GetMap().empty());
204  PHG4Particle *last_primary = _truth_container->GetMap().rbegin()->second;
205  assert(last_primary);
206 
207  return last_primary;
208 }
209 
210 int QAG4SimulationCalorimeterSum::Init_TrackProj(PHCompositeNode * /*topNode*/)
211 {
212  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
213  assert(hm);
214 
215  hm->registerHisto(
216  new TH2F(
217  TString(get_histo_prefix()) + TString(_calo_name_cemc.c_str()) + "_TrackProj", //
218  TString(_calo_name_cemc.c_str()) + " Tower Energy Distr. around Track Proj.;Polar distance / Tower width;Azimuthal distance / Tower width",
219  (Max_N_Tower - 1) * 10, -Max_N_Tower / 2, Max_N_Tower / 2,
220  (Max_N_Tower - 1) * 10, -Max_N_Tower / 2, Max_N_Tower / 2));
221 
222  hm->registerHisto(
223  new TH2F(
224  TString(get_histo_prefix()) + TString(_calo_name_hcalin.c_str()) + "_TrackProj", //
225  TString(_calo_name_hcalin.c_str()) + " Tower Energy Distr. around Track Proj.;Polar distance / Tower width;Azimuthal distance / Tower width",
226  (Max_N_Tower - 1) * 10, -Max_N_Tower / 2, Max_N_Tower / 2,
227  (Max_N_Tower - 1) * 10, -Max_N_Tower / 2, Max_N_Tower / 2));
228 
229  hm->registerHisto(
230  new TH2F(
231  TString(get_histo_prefix()) + TString(_calo_name_hcalout.c_str()) + "_TrackProj", //
232  TString(_calo_name_hcalout.c_str()) + " Tower Energy Distr. around Track Proj.;Polar distance / Tower width;Azimuthal distance / Tower width",
233  (Max_N_Tower - 1) * 10, -Max_N_Tower / 2, Max_N_Tower / 2,
234  (Max_N_Tower - 1) * 10, -Max_N_Tower / 2, Max_N_Tower / 2));
235 
236  hm->registerHisto(
237  new TH1F(TString(get_histo_prefix()) + "TrackProj_3x3Tower_EP", //
238  "Tower 3x3 sum /E_{Truth};#Sigma_{3x3}[E_{Tower}] / total truth energy",
239  150, 0, 1.5));
240 
241  hm->registerHisto(
242  new TH1F(TString(get_histo_prefix()) + "TrackProj_5x5Tower_EP", //
243  "Tower 5x5 sum /E_{Truth};#Sigma_{5x5}[E_{Tower}] / total truth energy",
244  150, 0, 1.5));
245 
246  return Fun4AllReturnCodes::EVENT_OK;
247 }
248 
249 int QAG4SimulationCalorimeterSum::process_event_TrackProj(PHCompositeNode *topNode)
250 {
251  if (Verbosity() > 2)
252  std::cout << "QAG4SimulationCalorimeterSum::process_event_TrackProj() entered"
253  << std::endl;
254 
255  PHG4Particle *primary = get_truth_particle();
256  if (!primary)
257  return Fun4AllReturnCodes::DISCARDEVENT;
258 
259  SvtxTrackEval *trackeval = _svtxevalstack->get_track_eval();
260  assert(trackeval);
261  SvtxTrack *track = trackeval->best_track_from(primary);
262  if (!track)
263  return Fun4AllReturnCodes::EVENT_OK; // not through the whole event for missing track.
264 
265  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
266  assert(hm);
267  TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(
268  get_histo_prefix() + "Normalization"));
269  assert(h_norm);
270  h_norm->Fill("Track", 1);
271 
272  {
273  TH1F *hsum = dynamic_cast<TH1F *>(hm->getHisto(
274  (get_histo_prefix()) + "TrackProj_3x3Tower_EP"));
275  assert(hsum);
276 
277  hsum->Fill(
278  (track->get_cal_energy_3x3(SvtxTrack::CEMC) + track->get_cal_energy_3x3(SvtxTrack::HCALIN) + track->get_cal_energy_3x3(SvtxTrack::HCALOUT)) / (primary->get_e() + 1e-9));
279  }
280  {
281  TH1F *hsum = dynamic_cast<TH1F *>(hm->getHisto(
282  (get_histo_prefix()) + "TrackProj_5x5Tower_EP"));
283  assert(hsum);
284 
285  hsum->Fill(
286  (track->get_cal_energy_5x5(SvtxTrack::CEMC) + track->get_cal_energy_5x5(SvtxTrack::HCALIN) + track->get_cal_energy_5x5(SvtxTrack::HCALOUT)) / (primary->get_e() + 1e-9));
287  }
288 
289  eval_trk_proj(_calo_name_cemc, track, topNode);
290  eval_trk_proj(_calo_name_hcalin, track, topNode);
291  eval_trk_proj(_calo_name_hcalout, track, topNode);
292 
293  return Fun4AllReturnCodes::EVENT_OK;
294 }
295 
296 bool QAG4SimulationCalorimeterSum::eval_trk_proj(const std::string &detector, SvtxTrack *track,
297  PHCompositeNode *topNode)
298 // Track projections
299 {
300  assert(track);
301  assert(topNode);
302 
303  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
304  assert(hm);
305 
306  TH2F *h2_proj = dynamic_cast<TH2F *>(hm->getHisto(
307  (get_histo_prefix()) + detector + "_TrackProj"));
308  assert(h2_proj);
309 
310  // pull the tower geometry
311  std::string towergeonodename = "TOWERGEOM_" + detector;
312  RawTowerGeomContainer *towergeo = findNode::getClass<RawTowerGeomContainer>(
313  topNode, towergeonodename);
314  assert(towergeo);
315 
316  if (Verbosity() > 2)
317  {
318  towergeo->identify();
319  }
320  // pull the towers
321  std::string towernodename = "TOWER_CALIB_" + detector;
322  RawTowerContainer *towerList = findNode::getClass<RawTowerContainer>(topNode,
323  towernodename);
324  assert(towerList);
325 
326  if (Verbosity() > 3)
327  {
328  std::cout << __PRETTY_FUNCTION__ << " - info - handling track track p = ("
329  << track->get_px() << ", " << track->get_py() << ", "
330  << track->get_pz() << "), v = (" << track->get_x() << ", "
331  << track->get_z() << ", " << track->get_z() << ")"
332  << " size_states "
333  << track->size_states() << std::endl;
334  }
335 
336  // curved tracks inside mag field
337  // straight projections thereafter
338 
339  std::vector<double> point;
340  point.assign(3, NAN);
341 
342  // const double radius = towergeo->get_radius() + towergeo->get_thickness() * 0.5;
343 
344  // PHG4HoughTransform::projectToRadius(track, _magField, radius, point);
345 
346  if (std::isnan(point[0]) or std::isnan(point[1]) or std::isnan(point[2]))
347  {
348  // std::cout << __PRETTY_FUNCTION__ << "::" << Name()
349  // << " - Error - track extrapolation failure:";
350  // track->identify();
351  return false;
352  }
353 
354  assert(not std::isnan(point[0]));
355  assert(not std::isnan(point[1]));
356  assert(not std::isnan(point[2]));
357 
358  double x = point[0];
359  double y = point[1];
360  double z = point[2];
361 
362  double phi = atan2(y, x);
363  double eta = asinh(z / sqrt(x * x + y * y));
364 
365  // calculate 3x3 tower energy
366  int binphi = towergeo->get_phibin(phi);
367  int bineta = towergeo->get_etabin(eta);
368 
369  double etabin_width = towergeo->get_etabounds(bineta).second - towergeo->get_etabounds(bineta).first;
370  if (bineta > 1 and bineta < towergeo->get_etabins() - 1)
371  etabin_width = (towergeo->get_etacenter(bineta + 1) - towergeo->get_etacenter(bineta - 1)) / 2.;
372 
373  double phibin_width = towergeo->get_phibounds(binphi).second - towergeo->get_phibounds(binphi).first;
374 
375  assert(etabin_width > 0);
376  assert(phibin_width > 0);
377 
378  const double etabin_shift = (towergeo->get_etacenter(bineta) - eta) / etabin_width;
379  const double phibin_shift = (fmod(
380  towergeo->get_phicenter(binphi) - phi + 5 * M_PI, 2 * M_PI) -
381  M_PI) /
382  phibin_width;
383 
384  if (Verbosity() > 3)
385  std::cout << __PRETTY_FUNCTION__ << " - info - handling track proj. (" << x
386  << ", " << y << ", " << z << ")"
387  << ", eta " << eta << ", phi" << phi
388  << ", bineta " << bineta << " - ["
389  << towergeo->get_etabounds(bineta).first << ", "
390  << towergeo->get_etabounds(bineta).second << "], etabin_shift = "
391  << etabin_shift << ", binphi " << binphi << " - ["
392  << towergeo->get_phibounds(binphi).first << ", "
393  << towergeo->get_phibounds(binphi).second << "], phibin_shift = "
394  << phibin_shift << std::endl;
395 
396  const int bin_search_range = (Max_N_Tower - 1) / 2;
397  for (int iphi = binphi - bin_search_range; iphi <= binphi + bin_search_range;
398  ++iphi)
399 
400  for (int ieta = bineta - bin_search_range;
401  ieta <= bineta + bin_search_range; ++ieta)
402  {
403  // wrap around
404  int wrapphi = iphi;
405  if (wrapphi < 0)
406  {
407  wrapphi = towergeo->get_phibins() + wrapphi;
408  }
409  if (wrapphi >= towergeo->get_phibins())
410  {
411  wrapphi = wrapphi - towergeo->get_phibins();
412  }
413 
414  // edges
415  if (ieta < 0)
416  continue;
417  if (ieta >= towergeo->get_etabins())
418  continue;
419 
420  if (Verbosity() > 3)
421  std::cout << __PRETTY_FUNCTION__ << " - info - processing tower"
422  << ", bineta " << ieta << " - ["
423  << towergeo->get_etabounds(ieta).first << ", "
424  << towergeo->get_etabounds(ieta).second << "]"
425  << ", binphi "
426  << wrapphi << " - [" << towergeo->get_phibounds(wrapphi).first
427  << ", " << towergeo->get_phibounds(wrapphi).second << "]" << std::endl;
428 
429  RawTower *tower = towerList->getTower(ieta, wrapphi);
430  double energy = 0;
431  if (tower)
432  {
433  if (Verbosity() > 1)
434  std::cout << __PRETTY_FUNCTION__ << " - info - tower " << ieta << " "
435  << wrapphi << " energy = " << tower->get_energy() << std::endl;
436 
437  energy = tower->get_energy();
438  }
439 
440  h2_proj->Fill(ieta - bineta + etabin_shift,
441  iphi - binphi + phibin_shift, energy);
442 
443  } // for (int ieta = bineta-1; ieta < bineta+2; ++ieta) {
444 
445  return true;
446 } // // Track projections
447 
448 //int
449 //QAG4SimulationCalorimeterSum::Init_Tower(PHCompositeNode *topNode)
450 //{
451 //
452 // Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
453 // assert(hm);
454 //
456 //
457 // return Fun4AllReturnCodes::EVENT_OK;
458 //}
459 //
460 //int
461 //QAG4SimulationCalorimeterSum::process_event_Tower(PHCompositeNode *topNode)
462 //{
463 // if (Verbosity() > 2)
464 // std::cout << "QAG4SimulationCalorimeterSum::process_event_Tower() entered"
465 // << std::endl;
466 //
467 // return Fun4AllReturnCodes::EVENT_OK;
468 //}
469 
470 int QAG4SimulationCalorimeterSum::Init_Cluster(PHCompositeNode * /*topNode*/)
471 {
472  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
473  assert(hm);
474 
475  hm->registerHisto(
476  new TH2F(
477  TString(get_histo_prefix()) + "Cluster_" + _calo_name_cemc.c_str() + "_" + _calo_name_hcalin.c_str(), //
478  TString(_calo_name_hcalin.c_str()) + " VS " + TString(_calo_name_cemc.c_str()) + ": best cluster energy;" + TString(_calo_name_cemc.c_str()) + " cluster energy (GeV);" + TString(_calo_name_hcalin.c_str()) + " cluster energy (GeV)",
479  70, 0, 70, 70, 0, 70));
480 
481  hm->registerHisto(
482  new TH2F(
483  TString(get_histo_prefix()) + "Cluster_" + _calo_name_cemc.c_str() + "_" + _calo_name_hcalin.c_str() + "_" + _calo_name_hcalout.c_str(), //
484  TString(_calo_name_cemc.c_str()) + " + " + TString(_calo_name_hcalin.c_str()) + " VS " + TString(_calo_name_hcalout.c_str()) + ": best cluster energy;" + TString(_calo_name_cemc.c_str()) + " + " + TString(_calo_name_hcalin.c_str()) + " cluster energy (GeV);" + TString(_calo_name_hcalout.c_str()) + " cluster energy (GeV)",
485  70, 0, 70, 70, 0, 70));
486 
487  hm->registerHisto(
488  new TH1F(TString(get_histo_prefix()) + "Cluster_EP", //
489  "Total Cluster E_{Reco}/E_{Truth};Reco cluster energy sum / total truth energy",
490  150, 0, 1.5));
491 
492  hm->registerHisto(
493  new TH1F(
494  TString(get_histo_prefix()) + "Cluster_Ratio_" + _calo_name_cemc.c_str() + "_" + _calo_name_hcalin.c_str(), //
495  "Energy ratio " + TString(_calo_name_cemc.c_str()) + " VS " + TString(_calo_name_hcalin.c_str()) + ";Best cluster " + TString(_calo_name_cemc.c_str()) + " / (" + TString(_calo_name_cemc.c_str()) + " + " + TString(_calo_name_hcalin.c_str()) + ")", 110, 0, 1.1));
496 
497  hm->registerHisto(
498  new TH1F(
499  TString(get_histo_prefix()) + "Cluster_Ratio_" + _calo_name_cemc.c_str() + "_" + _calo_name_hcalin.c_str() + "_" + TString(_calo_name_hcalout.c_str()), //
500  "Energy ratio " + TString(_calo_name_cemc.c_str()) + " + " + TString(_calo_name_hcalin.c_str()) + " VS " + TString(_calo_name_hcalout.c_str()) + ";Best cluster (" + TString(_calo_name_cemc.c_str()) + " + " + TString(_calo_name_hcalin.c_str()) + ") / (" + TString(_calo_name_cemc.c_str()) + " + " + TString(_calo_name_hcalin.c_str()) + " + " + TString(_calo_name_hcalout.c_str()) + ")", 110, 0, 1.1));
501 
502  return Fun4AllReturnCodes::EVENT_OK;
503 }
504 
505 int QAG4SimulationCalorimeterSum::process_event_Cluster(PHCompositeNode * /*topNode*/)
506 {
507  if (Verbosity() > 2)
508  std::cout << "QAG4SimulationCalorimeterSum::process_event_Cluster() entered"
509  << std::endl;
510 
511  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
512  assert(hm);
513 
514  PHG4Particle *primary = get_truth_particle();
515  if (!primary)
516  return Fun4AllReturnCodes::DISCARDEVENT;
517 
518  RawCluster *cluster_cemc =
519  _caloevalstack_cemc->get_rawcluster_eval()->best_cluster_from(primary);
520  RawCluster *cluster_hcalin =
521  _caloevalstack_hcalin->get_rawcluster_eval()->best_cluster_from(primary);
522  RawCluster *cluster_hcalout =
523  _caloevalstack_hcalout->get_rawcluster_eval()->best_cluster_from(primary);
524 
525  const double cluster_cemc_e = cluster_cemc ? cluster_cemc->get_energy() : 0;
526  const double cluster_hcalin_e =
527  cluster_hcalin ? cluster_hcalin->get_energy() : 0;
528  const double cluster_hcalout_e =
529  cluster_hcalout ? cluster_hcalout->get_energy() : 0;
530 
531  if (cluster_cemc_e + cluster_hcalin_e > 0)
532  {
533  TH2F *h2 = dynamic_cast<TH2F *>(hm->getHisto(
534  (get_histo_prefix()) + "Cluster_" + _calo_name_cemc + "_" + _calo_name_hcalin));
535  assert(h2);
536 
537  h2->Fill(cluster_cemc_e, cluster_hcalin_e);
538 
539  TH1F *hr = dynamic_cast<TH1F *>(hm->getHisto(
540  (get_histo_prefix()) + "Cluster_Ratio_" + _calo_name_cemc + "_" + _calo_name_hcalin));
541  assert(hr);
542 
543  hr->Fill(cluster_cemc_e / (cluster_cemc_e + cluster_hcalin_e));
544  }
545 
546  if (cluster_cemc_e + cluster_hcalin_e + cluster_hcalout_e > 0)
547  {
548  if (Verbosity())
549  {
550  std::cout << "QAG4SimulationCalorimeterSum::process_event_Cluster - "
551  << " cluster_cemc_e = " << cluster_cemc_e
552  << " cluster_hcalin_e = " << cluster_hcalin_e
553  << " cluster_hcalout_e = " << cluster_hcalout_e << " hr = "
554  << (cluster_cemc_e + cluster_hcalin_e) / (cluster_cemc_e + cluster_hcalin_e + cluster_hcalout_e)
555  << std::endl;
556  }
557 
558  TH2F *h2 = dynamic_cast<TH2F *>(hm->getHisto(
559  (get_histo_prefix()) + "Cluster_" + _calo_name_cemc + "_" + _calo_name_hcalin + "_" + _calo_name_hcalout));
560  assert(h2);
561 
562  h2->Fill((cluster_cemc_e + cluster_hcalin_e), cluster_hcalout_e);
563 
564  TH1F *hr = dynamic_cast<TH1F *>(hm->getHisto(
565  (get_histo_prefix()) + "Cluster_Ratio_" + _calo_name_cemc + "_" + _calo_name_hcalin + "_" + _calo_name_hcalout));
566  assert(hr);
567 
568  hr->Fill(
569  (cluster_cemc_e + cluster_hcalin_e) / (cluster_cemc_e + cluster_hcalin_e + cluster_hcalout_e));
570 
571  TH1F *hsum = dynamic_cast<TH1F *>(hm->getHisto(
572  (get_histo_prefix()) + "Cluster_EP"));
573  assert(hsum);
574 
575  hsum->Fill(
576  (cluster_cemc_e + cluster_hcalin_e + cluster_hcalout_e) / (primary->get_e() + 1e-9));
577  }
578 
579  return Fun4AllReturnCodes::EVENT_OK;
580 }