EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QAG4SimulationCalorimeter.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QAG4SimulationCalorimeter.cc
1 #include "QAG4SimulationCalorimeter.h"
2 
3 #include "QAHistManagerDef.h"
4 
5 #include <g4main/PHG4Hit.h>
6 #include <g4main/PHG4HitContainer.h>
7 #include <g4main/PHG4Particle.h>
8 #include <g4main/PHG4TruthInfoContainer.h>
9 #include <g4main/PHG4VtxPoint.h>
10 
11 #include <calobase/RawCluster.h>
12 #include <calobase/RawClusterContainer.h>
13 #include <calobase/RawTower.h>
14 #include <calobase/RawTowerContainer.h>
15 #include <calobase/RawTowerGeomContainer.h>
16 
17 #include <g4eval/CaloEvalStack.h>
18 #include <g4eval/CaloRawClusterEval.h>
19 
20 #include <fun4all/Fun4AllHistoManager.h>
21 #include <fun4all/Fun4AllReturnCodes.h>
22 #include <fun4all/SubsysReco.h>
23 
24 #include <phool/PHCompositeNode.h>
25 #include <phool/PHNodeIterator.h>
26 #include <phool/getClass.h>
27 #include <phool/phool.h>
28 
29 #include <TAxis.h>
30 #include <TH1.h>
31 #include <TH2.h>
32 #include <TNamed.h>
33 #include <TString.h>
34 #include <TVector3.h>
35 
36 #include <CLHEP/Vector/ThreeVector.h> // for Hep3Vector
37 
38 #include <cassert>
39 #include <cmath>
40 #include <cstdlib>
41 #include <iostream>
42 #include <iterator> // for reverse_iterator
43 #include <map>
44 #include <utility>
45 
46 QAG4SimulationCalorimeter::QAG4SimulationCalorimeter(const std::string &calo_name,
47  QAG4SimulationCalorimeter::enu_flags flags)
48  : SubsysReco("QAG4SimulationCalorimeter_" + calo_name)
49  , _calo_name(calo_name)
50  , _flags(flags)
51  , _calo_hit_container(nullptr)
52  , _calo_abs_hit_container(nullptr)
53  , _truth_container(nullptr)
54 {
55 }
56 
57 int QAG4SimulationCalorimeter::InitRun(PHCompositeNode *topNode)
58 {
59  PHNodeIterator iter(topNode);
60  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
61  assert(dstNode);
62 
63  if (flag(kProcessG4Hit))
64  {
65  _calo_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
66  "G4HIT_" + _calo_name);
67  if (!_calo_hit_container)
68  {
69  std::cout << "QAG4SimulationCalorimeter::InitRun - Fatal Error - "
70  << "unable to find DST node "
71  << "G4HIT_" + _calo_name << std::endl;
72  assert(_calo_hit_container);
73  }
74 
75  _calo_abs_hit_container = findNode::getClass<PHG4HitContainer>(topNode,
76  "G4HIT_ABSORBER_" + _calo_name);
77  if (!_calo_abs_hit_container)
78  {
79  std::cout << "QAG4SimulationCalorimeter::InitRun - Fatal Error - "
80  << "unable to find DST node "
81  << "G4HIT_ABSORBER_" + _calo_name
82  << std::endl;
83  assert(_calo_abs_hit_container);
84  }
85  }
86 
87  _truth_container = findNode::getClass<PHG4TruthInfoContainer>(topNode,
88  "G4TruthInfo");
89  if (!_truth_container)
90  {
91  std::cout << "QAG4SimulationCalorimeter::InitRun - Fatal Error - "
92  << "unable to find DST node "
93  << "G4TruthInfo" << std::endl;
94  assert(_truth_container);
95  }
96 
97  if (flag(kProcessCluster))
98  {
99  if (!_caloevalstack)
100  {
101  _caloevalstack.reset(new CaloEvalStack(topNode, _calo_name));
102  _caloevalstack->set_strict(true);
103  _caloevalstack->set_verbosity(Verbosity() + 1);
104  }
105  else
106  {
107  _caloevalstack->next_event(topNode);
108  }
109  }
110  return Fun4AllReturnCodes::EVENT_OK;
111 }
112 
113 int QAG4SimulationCalorimeter::Init(PHCompositeNode *topNode)
114 {
115  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
116  assert(hm);
117  TH1D *h = new TH1D(TString(get_histo_prefix()) + "_Normalization", //
118  TString(_calo_name) + " Normalization;Items;Count", 10, .5, 10.5);
119  int i = 1;
120  h->GetXaxis()->SetBinLabel(i++, "Event");
121  h->GetXaxis()->SetBinLabel(i++, "G4Hit Active");
122  h->GetXaxis()->SetBinLabel(i++, "G4Hit Absor.");
123  h->GetXaxis()->SetBinLabel(i++, "Tower");
124  h->GetXaxis()->SetBinLabel(i++, "Tower Hit");
125  h->GetXaxis()->SetBinLabel(i++, "Cluster");
126  h->GetXaxis()->LabelsOption("v");
127  hm->registerHisto(h);
128 
129  if (flag(kProcessG4Hit))
130  {
131  if (Verbosity() >= 1)
132  std::cout << "QAG4SimulationCalorimeter::Init - Process sampling fraction"
133  << std::endl;
134  Init_G4Hit(topNode);
135  }
136  if (flag(kProcessTower))
137  {
138  if (Verbosity() >= 1)
139  std::cout << "QAG4SimulationCalorimeter::Init - Process tower occupancies"
140  << std::endl;
141  Init_Tower(topNode);
142  }
143  if (flag(kProcessCluster))
144  {
145  if (Verbosity() >= 1)
146  std::cout << "QAG4SimulationCalorimeter::Init - Process tower occupancies"
147  << std::endl;
148  Init_Cluster(topNode);
149  }
150 
151  return Fun4AllReturnCodes::EVENT_OK;
152 }
153 
154 int QAG4SimulationCalorimeter::process_event(PHCompositeNode *topNode)
155 {
156  if (Verbosity() > 2)
157  std::cout << "QAG4SimulationCalorimeter::process_event() entered" << std::endl;
158 
159  if (_caloevalstack)
160  _caloevalstack->next_event(topNode);
161 
162  if (flag(kProcessG4Hit))
163  {
164  int ret = process_event_G4Hit(topNode);
165 
166  if (ret != Fun4AllReturnCodes::EVENT_OK)
167  return ret;
168  }
169 
170  if (flag(kProcessTower))
171  {
172  int ret = process_event_Tower(topNode);
173 
174  if (ret != Fun4AllReturnCodes::EVENT_OK)
175  return ret;
176  }
177 
178  if (flag(kProcessCluster))
179  {
180  int ret = process_event_Cluster(topNode);
181 
182  if (ret != Fun4AllReturnCodes::EVENT_OK)
183  return ret;
184  }
185 
186  // at the end, count success events
187  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
188  assert(hm);
189  TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(
190  get_histo_prefix() + "_Normalization"));
191  assert(h_norm);
192  h_norm->Fill("Event", 1);
193 
194  return Fun4AllReturnCodes::EVENT_OK;
195 }
196 
197 std::string
198 QAG4SimulationCalorimeter::get_histo_prefix()
199 {
200  return "h_QAG4Sim_" + std::string(_calo_name);
201 }
202 
203 int QAG4SimulationCalorimeter::Init_G4Hit(PHCompositeNode * /*topNode*/)
204 {
205  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
206  assert(hm);
207 
208  hm->registerHisto(
209  new TH2F(TString(get_histo_prefix()) + "_G4Hit_RZ", //
210  TString(_calo_name) + " RZ projection;G4 Hit Z (cm);G4 Hit R (cm)", 1200, -300, 300,
211  600, -000, 300));
212 
213  hm->registerHisto(
214  new TH2F(TString(get_histo_prefix()) + "_G4Hit_XY", //
215  TString(_calo_name) + " XY projection;G4 Hit X (cm);G4 Hit Y (cm)", 1200, -300, 300,
216  1200, -300, 300));
217 
218  hm->registerHisto(
219  new TH2F(TString(get_histo_prefix()) + "_G4Hit_LateralTruthProjection", //
220  TString(_calo_name) + " shower lateral projection (last primary);Polar direction (cm);Azimuthal direction (cm)",
221  200, -30, 30, 200, -30, 30));
222 
223  hm->registerHisto(new TH1F(TString(get_histo_prefix()) + "_G4Hit_SF", //
224  TString(_calo_name) + " sampling fraction;Sampling fraction", 1000, 0, .2));
225 
226  hm->registerHisto(
227  new TH1F(TString(get_histo_prefix()) + "_G4Hit_VSF", //
228  TString(_calo_name) + " visible sampling fraction;Visible sampling fraction", 1000, 0,
229  .2));
230 
231  TH1F *h =
232  new TH1F(TString(get_histo_prefix()) + "_G4Hit_HitTime", //
233  TString(_calo_name) + " hit time (edep weighting);Hit time - T0 (ns);Geant4 energy density",
234  1000, 0.5, 10000);
235  QAHistManagerDef::useLogBins(h->GetXaxis());
236  hm->registerHisto(h);
237 
238  hm->registerHisto(
239  new TH1F(TString(get_histo_prefix()) + "_G4Hit_FractionTruthEnergy", //
240  TString(_calo_name) + " fraction truth energy ;G4 edep / particle energy",
241  1000, 0, 1));
242 
243  hm->registerHisto(
244  new TH1F(TString(get_histo_prefix()) + "_G4Hit_FractionEMVisibleEnergy", //
245  TString(_calo_name) + " fraction visible energy from EM; visible energy from e^{#pm} / total visible energy",
246  100, 0, 1));
247 
248  return Fun4AllReturnCodes::EVENT_OK;
249 }
250 
251 int QAG4SimulationCalorimeter::process_event_G4Hit(PHCompositeNode * /*topNode*/)
252 {
253  if (Verbosity() > 2)
254  std::cout << "QAG4SimulationCalorimeter::process_event_G4Hit() entered" << std::endl;
255 
256  TH1F *h = nullptr;
257 
258  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
259  assert(hm);
260 
261  TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(
262  get_histo_prefix() + "_Normalization"));
263  assert(h_norm);
264 
265  // get primary
266  assert(_truth_container);
267  PHG4TruthInfoContainer::ConstRange primary_range =
268  _truth_container->GetPrimaryParticleRange();
269  double total_primary_energy = 1e-9; //make it zero energy epsilon samll so it can be used for denominator
270  for (PHG4TruthInfoContainer::ConstIterator particle_iter = primary_range.first;
271  particle_iter != primary_range.second; ++particle_iter)
272  {
273  const PHG4Particle *particle = particle_iter->second;
274  assert(particle);
275  total_primary_energy += particle->get_e();
276  }
277 
278  assert(not _truth_container->GetMap().empty());
279  const PHG4Particle *last_primary =
280  _truth_container->GetMap().rbegin()->second;
281  assert(last_primary);
282 
283  if (Verbosity() > 2)
284  {
285  std::cout
286  << "QAG4SimulationCalorimeter::process_event_G4Hit() handle this truth particle"
287  << std::endl;
288  last_primary->identify();
289  }
290  const PHG4VtxPoint *primary_vtx = //
291  _truth_container->GetPrimaryVtx(last_primary->get_vtx_id());
292  assert(primary_vtx);
293  if (Verbosity() > 2)
294  {
295  std::cout
296  << "QAG4SimulationCalorimeter::process_event_G4Hit() handle this vertex"
297  << std::endl;
298  primary_vtx->identify();
299  }
300 
301  const double t0 = primary_vtx->get_t();
302  const TVector3 vertex(primary_vtx->get_x(), primary_vtx->get_y(),
303  primary_vtx->get_z());
304 
305  // projection axis
306  TVector3 axis_proj(last_primary->get_px(), last_primary->get_py(),
307  last_primary->get_pz());
308  if (axis_proj.Mag() == 0)
309  axis_proj.SetXYZ(0, 0, 1);
310  axis_proj = axis_proj.Unit();
311 
312  // azimuthal direction axis
313  TVector3 axis_azimuth = axis_proj.Cross(TVector3(0, 0, 1));
314  if (axis_azimuth.Mag() == 0)
315  axis_azimuth.SetXYZ(1, 0, 0);
316  axis_azimuth = axis_azimuth.Unit();
317 
318  // polar direction axis
319  TVector3 axis_polar = axis_proj.Cross(axis_azimuth);
320  assert(axis_polar.Mag() > 0);
321  axis_polar = axis_polar.Unit();
322 
323  double e_calo = 0.0; // active energy deposition
324  double ev_calo = 0.0; // visible energy
325  double ea_calo = 0.0; // absorber energy
326  double ev_calo_em = 0.0; // EM visible energy
327 
328  if (_calo_hit_container)
329  {
330  TH2F *hrz = dynamic_cast<TH2F *>(hm->getHisto(
331  get_histo_prefix() + "_G4Hit_RZ"));
332  assert(hrz);
333  TH2F *hxy = dynamic_cast<TH2F *>(hm->getHisto(
334  get_histo_prefix() + "_G4Hit_XY"));
335  assert(hxy);
336  TH1F *ht = dynamic_cast<TH1F *>(hm->getHisto(
337  get_histo_prefix() + "_G4Hit_HitTime"));
338  assert(ht);
339  TH2F *hlat = dynamic_cast<TH2F *>(hm->getHisto(
340  get_histo_prefix() + "_G4Hit_LateralTruthProjection"));
341  assert(hlat);
342 
343  h_norm->Fill("G4Hit Active", _calo_hit_container->size());
344  PHG4HitContainer::ConstRange calo_hit_range =
345  _calo_hit_container->getHits();
346  for (PHG4HitContainer::ConstIterator hit_iter = calo_hit_range.first;
347  hit_iter != calo_hit_range.second; hit_iter++)
348  {
349  PHG4Hit *this_hit = hit_iter->second;
350  assert(this_hit);
351 
352  e_calo += this_hit->get_edep();
353  ev_calo += this_hit->get_light_yield();
354 
355  // EM visible energy that is only associated with electron energy deposition
356  PHG4Particle *particle = _truth_container->GetParticle(
357  this_hit->get_trkid());
358  if (!particle)
359  {
360  std::cout << __PRETTY_FUNCTION__ << " - Error - this PHG4hit missing particle: ";
361  this_hit->identify();
362  }
363  assert(particle);
364  if (abs(particle->get_pid()) == 11)
365  ev_calo_em += this_hit->get_light_yield();
366 
367  const TVector3 hit(this_hit->get_avg_x(), this_hit->get_avg_y(),
368  this_hit->get_avg_z());
369 
370  hrz->Fill(hit.Z(), hit.Perp(), this_hit->get_edep());
371  hxy->Fill(hit.X(), hit.Y(), this_hit->get_edep());
372  ht->Fill(this_hit->get_avg_t() - t0, this_hit->get_edep());
373 
374  const double hit_azimuth = axis_azimuth.Dot(hit - vertex);
375  const double hit_polar = axis_polar.Dot(hit - vertex);
376  hlat->Fill(hit_polar, hit_azimuth, this_hit->get_edep());
377  }
378  }
379 
380  if (_calo_abs_hit_container)
381  {
382  h_norm->Fill("G4Hit Absor.", _calo_abs_hit_container->size());
383 
384  PHG4HitContainer::ConstRange calo_abs_hit_range =
385  _calo_abs_hit_container->getHits();
386  for (PHG4HitContainer::ConstIterator hit_iter = calo_abs_hit_range.first;
387  hit_iter != calo_abs_hit_range.second; hit_iter++)
388  {
389  PHG4Hit *this_hit = hit_iter->second;
390  assert(this_hit);
391 
392  ea_calo += this_hit->get_edep();
393  }
394  }
395 
396  if (Verbosity() > 3)
397  std::cout << "QAG4SimulationCalorimeter::process_event_G4Hit::" << _calo_name
398  << " - SF = " << e_calo / (e_calo + ea_calo + 1e-9) << ", VSF = "
399  << ev_calo / (e_calo + ea_calo + 1e-9) << std::endl;
400 
401  if (e_calo + ea_calo > 0)
402  {
403  h = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "_G4Hit_SF"));
404  assert(h);
405  h->Fill(e_calo / (e_calo + ea_calo));
406 
407  h = dynamic_cast<TH1F *>(hm->getHisto(get_histo_prefix() + "_G4Hit_VSF"));
408  assert(h);
409  h->Fill(ev_calo / (e_calo + ea_calo));
410  }
411 
412  h = dynamic_cast<TH1F *>(hm->getHisto(
413  get_histo_prefix() + "_G4Hit_FractionTruthEnergy"));
414  assert(h);
415  h->Fill((e_calo + ea_calo) / total_primary_energy);
416 
417  if (ev_calo > 0)
418  {
419  h = dynamic_cast<TH1F *>(hm->getHisto(
420  get_histo_prefix() + "_G4Hit_FractionEMVisibleEnergy"));
421  assert(h);
422  h->Fill(ev_calo_em / (ev_calo));
423  }
424 
425  if (Verbosity() > 3)
426  std::cout << "QAG4SimulationCalorimeter::process_event_G4Hit::" << _calo_name
427  << " - histogram " << h->GetName() << " Get Sum = " << h->GetSum()
428  << std::endl;
429 
430  return Fun4AllReturnCodes::EVENT_OK;
431 }
432 
433 int QAG4SimulationCalorimeter::Init_Tower(PHCompositeNode * /*topNode*/)
434 {
435  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
436  assert(hm);
437 
438  TH1F *h = new TH1F(TString(get_histo_prefix()) + "_Tower_1x1", //
439  TString(_calo_name) + " 1x1 tower;1x1 TOWER Energy (GeV)", 100, 9e-4, 100);
440  QAHistManagerDef::useLogBins(h->GetXaxis());
441  hm->registerHisto(h);
442 
443  hm->registerHisto(
444  new TH1F(TString(get_histo_prefix()) + "_Tower_1x1_max", //
445  TString(_calo_name) + " 1x1 tower max per event;1x1 tower max per event (GeV)", 5000,
446  0, 50));
447 
448  h = new TH1F(TString(get_histo_prefix()) + "_Tower_2x2", //
449  TString(_calo_name) + " 2x2 tower;2x2 TOWER Energy (GeV)", 100, 9e-4, 100);
450  QAHistManagerDef::useLogBins(h->GetXaxis());
451  hm->registerHisto(h);
452  hm->registerHisto(
453  new TH1F(TString(get_histo_prefix()) + "_Tower_2x2_max", //
454  TString(_calo_name) + " 2x2 tower max per event;2x2 tower max per event (GeV)", 5000,
455  0, 50));
456 
457  h = new TH1F(TString(get_histo_prefix()) + "_Tower_3x3", //
458  TString(_calo_name) + " 3x3 tower;3x3 TOWER Energy (GeV)", 100, 9e-4, 100);
459  QAHistManagerDef::useLogBins(h->GetXaxis());
460  hm->registerHisto(h);
461  hm->registerHisto(
462  new TH1F(TString(get_histo_prefix()) + "_Tower_3x3_max", //
463  TString(_calo_name) + " 3x3 tower max per event;3x3 tower max per event (GeV)", 5000,
464  0, 50));
465 
466  h = new TH1F(TString(get_histo_prefix()) + "_Tower_4x4", //
467  TString(_calo_name) + " 4x4 tower;4x4 TOWER Energy (GeV)", 100, 9e-4, 100);
468  QAHistManagerDef::useLogBins(h->GetXaxis());
469  hm->registerHisto(h);
470  hm->registerHisto(
471  new TH1F(TString(get_histo_prefix()) + "_Tower_4x4_max", //
472  TString(_calo_name) + " 4x4 tower max per event;4x4 tower max per event (GeV)", 5000,
473  0, 50));
474 
475  h = new TH1F(TString(get_histo_prefix()) + "_Tower_5x5", //
476  TString(_calo_name) + " 5x5 tower;5x5 TOWER Energy (GeV)", 100, 9e-4, 100);
477  QAHistManagerDef::useLogBins(h->GetXaxis());
478  hm->registerHisto(h);
479  hm->registerHisto(
480  new TH1F(TString(get_histo_prefix()) + "_Tower_5x5_max", //
481  TString(_calo_name) + " 5x5 tower max per event;5x5 tower max per event (GeV)", 5000,
482  0, 50));
483 
484  return Fun4AllReturnCodes::EVENT_OK;
485 }
486 
487 int QAG4SimulationCalorimeter::process_event_Tower(PHCompositeNode *topNode)
488 {
489  const std::string detector(_calo_name);
490 
491  if (Verbosity() > 2)
492  std::cout << "QAG4SimulationCalorimeter::process_event_Tower() entered" << std::endl;
493 
494  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
495  assert(hm);
496  TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(
497  get_histo_prefix() + "_Normalization"));
498  assert(h_norm);
499 
500  std::string towernodename = "TOWER_CALIB_" + detector;
501  // Grab the towers
502  RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode,
503  towernodename);
504  if (!towers)
505  {
506  std::cout << PHWHERE << ": Could not find node " << towernodename
507  << std::endl;
508  return Fun4AllReturnCodes::ABORTRUN;
509  }
510  std::string towergeomnodename = "TOWERGEOM_" + detector;
511  RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(
512  topNode, towergeomnodename);
513  if (!towergeom)
514  {
515  std::cout << PHWHERE << ": Could not find node " << towergeomnodename
516  << std::endl;
517  return Fun4AllReturnCodes::ABORTRUN;
518  }
519 
520  static const int max_size = 5;
521  std::map<int, std::string> size_label;
522  size_label[1] = "1x1";
523  size_label[2] = "2x2";
524  size_label[3] = "3x3";
525  size_label[4] = "4x4";
526  size_label[5] = "5x5";
527  std::map<int, double> max_energy;
528  std::map<int, TH1F *> energy_hist_list;
529  std::map<int, TH1F *> max_energy_hist_list;
530 
531  for (int size = 1; size <= max_size; ++size)
532  {
533  max_energy[size] = 0;
534 
535  TH1F *h = dynamic_cast<TH1F *>(hm->getHisto(
536  get_histo_prefix() + "_Tower_" + size_label[size]));
537  assert(h);
538  energy_hist_list[size] = h;
539  h = dynamic_cast<TH1F *>(hm->getHisto(
540  get_histo_prefix() + "_Tower_" + size_label[size] + "_max"));
541  assert(h);
542  max_energy_hist_list[size] = h;
543  }
544 
545  h_norm->Fill("Tower", towergeom->size()); // total tower count
546  h_norm->Fill("Tower Hit", towers->size());
547 
548  for (int binphi = 0; binphi < towergeom->get_phibins(); ++binphi)
549  {
550  for (int bineta = 0; bineta < towergeom->get_etabins(); ++bineta)
551  {
552  for (int size = 1; size <= max_size; ++size)
553  {
554  // for 2x2 and 4x4 use slide-2 window as implemented in DAQ
555  if ((size == 2 or size == 4) and ((binphi % 2 != 0) and (bineta % 2 != 0)))
556  continue;
557 
558  double energy = 0;
559 
560  // sliding window made from 2x2 sums
561  for (int iphi = binphi; iphi < binphi + size; ++iphi)
562  {
563  for (int ieta = bineta; ieta < bineta + size; ++ieta)
564  {
565  if (ieta > towergeom->get_etabins())
566  continue;
567 
568  // wrap around
569  int wrapphi = iphi;
570  assert(wrapphi >= 0);
571  if (wrapphi >= towergeom->get_phibins())
572  {
573  wrapphi = wrapphi - towergeom->get_phibins();
574  }
575 
576  RawTower *tower = towers->getTower(ieta, wrapphi);
577 
578  if (tower)
579  {
580  const double e_intput = tower->get_energy();
581 
582  energy += e_intput;
583  }
584  }
585  }
586 
587  energy_hist_list[size]->Fill(energy == 0 ? 9.1e-4 : energy); // trick to fill 0 energy tower to the first bin
588 
589  if (energy > max_energy[size])
590  max_energy[size] = energy;
591 
592  } // for (int size = 1; size <= 4; ++size)
593  }
594  }
595 
596  for (int size = 1; size <= max_size; ++size)
597  {
598  max_energy_hist_list[size]->Fill(max_energy[size]);
599  }
600  return Fun4AllReturnCodes::EVENT_OK;
601 }
602 
603 int QAG4SimulationCalorimeter::Init_Cluster(PHCompositeNode * /*topNode*/)
604 {
605  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
606  assert(hm);
607 
608  hm->registerHisto(
609  new TH1F(TString(get_histo_prefix()) + "_Cluster_BestMatchERatio", //
610  TString(_calo_name) + " best matched cluster E/E_{Truth};E_{Cluster}/E_{Truth}", 150,
611  0, 1.5));
612 
613  hm->registerHisto(
614  new TH2F(TString(get_histo_prefix()) + "_Cluster_LateralTruthProjection", //
615  TString(_calo_name) + " best cluster lateral projection (last primary);Polar direction (cm);Azimuthal direction (cm)",
616  200, -15, 15, 200, -15, 15));
617 
618  return Fun4AllReturnCodes::EVENT_OK;
619 }
620 
621 int QAG4SimulationCalorimeter::process_event_Cluster(PHCompositeNode *topNode)
622 {
623  if (Verbosity() > 2)
624  {
625  std::cout << "QAG4SimulationCalorimeter::process_event_Cluster() entered"
626  << std::endl;
627  }
628 
629  Fun4AllHistoManager *hm = QAHistManagerDef::getHistoManager();
630  assert(hm);
631  std::string towergeomnodename = "TOWERGEOM_" + _calo_name;
632  RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(
633  topNode, towergeomnodename);
634  if (!towergeom)
635  {
636  std::cout << PHWHERE << ": Could not find node " << towergeomnodename << std::endl;
637  return Fun4AllReturnCodes::ABORTRUN;
638  }
639 
640  //get a cluster count
641  TH1D *h_norm = dynamic_cast<TH1D *>(hm->getHisto(get_histo_prefix() + "_Normalization"));
642  assert(h_norm);
643 
644  std::string nodename = "CLUSTER_" + _calo_name;
645  RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, nodename);
646  assert(clusters);
647  h_norm->Fill("Cluster", clusters->size());
648 
649  // get primary
650  assert(_truth_container);
651  assert(not _truth_container->GetMap().empty());
652  PHG4Particle *last_primary = _truth_container->GetMap().rbegin()->second;
653  assert(last_primary);
654 
655  if (Verbosity() > 2)
656  {
657  std::cout
658  << "QAG4SimulationCalorimeter::process_event_Cluster() handle this truth particle"
659  << std::endl;
660  last_primary->identify();
661  }
662 
663  assert(_caloevalstack);
664  CaloRawClusterEval *clustereval = _caloevalstack->get_rawcluster_eval();
665  assert(clustereval);
666 
667  TH1F *h = dynamic_cast<TH1F *>(hm->getHisto(
668  get_histo_prefix() + "_Cluster_BestMatchERatio"));
669  assert(h);
670 
671  RawCluster *cluster = clustereval->best_cluster_from(last_primary);
672  if (cluster)
673  {
674  // has a cluster matched and best cluster selected
675 
676  if (Verbosity() > 3)
677  std::cout << "QAG4SimulationCalorimeter::process_event_Cluster::"
678  << _calo_name << " - get cluster with energy "
679  << cluster->get_energy() << " VS primary energy "
680  << last_primary->get_e() << std::endl;
681 
682  h->Fill(cluster->get_energy() / (last_primary->get_e() + 1e-9)); //avoids divide zero
683 
684  // now work on the projection:
685  const CLHEP::Hep3Vector hit(cluster->get_position());
686 
687  const PHG4VtxPoint *primary_vtx = //
688  _truth_container->GetPrimaryVtx(last_primary->get_vtx_id());
689  assert(primary_vtx);
690  if (Verbosity() > 2)
691  {
692  std::cout
693  << "QAG4SimulationCalorimeter::process_event_Cluster() handle this vertex"
694  << std::endl;
695  primary_vtx->identify();
696  }
697 
698  const CLHEP::Hep3Vector vertex(primary_vtx->get_x(), primary_vtx->get_y(),
699  primary_vtx->get_z());
700 
701  // projection axis
702  CLHEP::Hep3Vector axis_proj(last_primary->get_px(), last_primary->get_py(),
703  last_primary->get_pz());
704  if (axis_proj.mag() == 0)
705  axis_proj.set(0, 0, 1);
706  axis_proj = axis_proj.unit();
707 
708  // azimuthal direction axis
709  CLHEP::Hep3Vector axis_azimuth = axis_proj.cross(CLHEP::Hep3Vector(0, 0, 1));
710  if (axis_azimuth.mag() == 0)
711  axis_azimuth.set(1, 0, 0);
712  axis_azimuth = axis_azimuth.unit();
713 
714  // polar direction axis
715  CLHEP::Hep3Vector axis_polar = axis_proj.cross(axis_azimuth);
716  assert(axis_polar.mag() > 0);
717  axis_polar = axis_polar.unit();
718 
719  TH2F *hlat = dynamic_cast<TH2F *>(hm->getHisto(
720  get_histo_prefix() + "_Cluster_LateralTruthProjection"));
721  assert(hlat);
722 
723  const double hit_azimuth = axis_azimuth.dot(hit - vertex);
724  const double hit_polar = axis_polar.dot(hit - vertex);
725  hlat->Fill(hit_polar, hit_azimuth);
726  }
727  else
728  {
729  if (Verbosity() > 3)
730  std::cout << "QAG4SimulationCalorimeter::process_event_Cluster::"
731  << _calo_name << " - missing cluster !";
732  h->Fill(0); // no cluster matched
733  }
734 
735  return Fun4AllReturnCodes::EVENT_OK;
736 }