ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CaloRawTowerEval.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CaloRawTowerEval.cc
1 
2 #include "CaloRawTowerEval.h"
3 #include "CaloTruthEval.h"
4 
5 #include <calobase/RawTower.h>
6 #include <calobase/RawTowerContainer.h>
7 #include <g4detectors/PHG4Cell.h>
10 #include <g4main/PHG4Shower.h>
12 
13 #include <phool/getClass.h>
14 
15 #include <cfloat>
16 #include <cassert>
17 #include <cmath>
18 #include <iostream>
19 #include <map>
20 #include <set>
21 #include <string>
22 
23 class PHG4Hit;
24 
25 using namespace std;
26 
27 CaloRawTowerEval::CaloRawTowerEval(PHCompositeNode* topNode, const std::string& caloname)
28  : _caloname(caloname)
29  , _trutheval(topNode, caloname)
30  , _towers(nullptr)
31  , _g4cells(nullptr)
32  , _g4hits(nullptr)
33  , _truthinfo(nullptr)
34  , _strict(false)
35  , _verbosity(1)
36  , _errors(0)
37  , _do_cache(true)
38  , _cache_all_truth_primary_showers()
39  , _cache_max_truth_primary_shower_by_energy()
40  , _cache_all_towers_from_primary_shower()
41  , _cache_best_tower_from_primary_shower()
42  , _cache_get_energy_contribution_primary_shower()
43  , _cache_all_truth_primary_particles()
44  , _cache_max_truth_primary_particle_by_energy()
45  , _cache_all_towers_from_primary_particle()
46  , _cache_best_tower_from_primary_particle()
47  , _cache_get_energy_contribution_primary_particle()
48  , _cache_all_truth_hits()
49 {
50  get_node_pointers(topNode);
51 }
52 
54 {
55  if (_verbosity > 0)
56  {
57  if ((_errors > 0) || (_verbosity > 1))
58  {
59  cout << "CaloRawTowerEval::~CaloRawTowerEval() - Error Count: " << _errors << endl;
60  }
61  }
62 }
63 
65 {
71 
77 
78  _cache_all_truth_hits.clear();
79 
80  _trutheval.next_event(topNode);
81 
82  get_node_pointers(topNode);
83 }
84 
86 {
87  if (!get_truth_eval()->has_reduced_node_pointers()) return false;
88 
89  if (_strict)
90  assert(_towers);
91  else if (!_towers)
92  return false;
93 
94  if (_strict)
95  assert(_truthinfo);
96  else if (!_truthinfo)
97  return false;
98 
99  return true;
100 }
101 
103 {
105  {
106  ++_errors;
107  return std::set<PHG4Shower*>();
108  }
109 
110  if (_strict)
111  {
112  assert(tower);
113  }
114  else if (!tower)
115  {
116  ++_errors;
117  return std::set<PHG4Shower*>();
118  }
119 
120  if (_do_cache)
121  {
122  std::map<RawTower*, std::set<PHG4Shower*> >::iterator iter =
124  if (iter != _cache_all_truth_primary_showers.end())
125  {
126  return iter->second;
127  }
128  }
129 
130  std::set<PHG4Shower*> showers;
131 
132  RawTower::ShowerConstRange shower_range = tower->get_g4showers();
133  for (RawTower::ShowerConstIterator iter = shower_range.first;
134  iter != shower_range.second;
135  ++iter)
136  {
137  PHG4Shower* shower = _truthinfo->GetShower(iter->first);
138 
139  if (_strict)
140  {
141  assert(shower);
142  }
143  else if (!shower)
144  {
145  ++_errors;
146  continue;
147  }
148 
149  showers.insert(shower);
150  }
151 
152  if (_do_cache) _cache_all_truth_primary_showers.insert(make_pair(tower, showers));
153 
154  return showers;
155 }
156 
158 {
160  {
161  ++_errors;
162  return nullptr;
163  }
164 
165  if (_strict)
166  {
167  assert(tower);
168  }
169  else if (!tower)
170  {
171  ++_errors;
172  return nullptr;
173  }
174 
175  if (_do_cache)
176  {
177  std::map<RawTower*, PHG4Shower*>::iterator iter =
180  {
181  return iter->second;
182  }
183  }
184 
185  PHG4Shower* max_shower = nullptr;
186  float max_e = FLT_MAX * -1.0;
187  std::set<PHG4Shower*> showers = all_truth_primary_showers(tower);
188 
189  for (std::set<PHG4Shower*>::iterator iter = showers.begin();
190  iter != showers.end();
191  ++iter)
192  {
193  PHG4Shower* shower = *iter;
194 
195  if (_strict)
196  {
197  assert(shower);
198  }
199  else if (!shower)
200  {
201  ++_errors;
202  continue;
203  }
204 
205  float e = get_energy_contribution(tower, shower);
206  if (isnan(e)) continue;
207  if (e > max_e)
208  {
209  max_e = e;
210  max_shower = shower;
211  }
212  }
213 
214  if (_do_cache) _cache_max_truth_primary_shower_by_energy.insert(make_pair(tower, max_shower));
215 
216  return max_shower;
217 }
218 
220 {
222  {
223  ++_errors;
224  return nullptr;
225  }
226 
227  if (_strict)
228  {
229  assert(shower);
230  }
231  else if (!shower)
232  {
233  ++_errors;
234  return nullptr;
235  }
236 
237  if (!_trutheval.is_primary(shower)) return nullptr;
238 
239  if (_do_cache)
240  {
241  std::map<PHG4Shower*, RawTower*>::iterator iter =
243  if (iter != _cache_best_tower_from_primary_shower.end())
244  {
245  return iter->second;
246  }
247  }
248 
249  RawTower* best_tower = nullptr;
250  float best_energy = FLT_MAX * -1.0;
251  std::set<RawTower*> towers = all_towers_from(shower);
252  for (std::set<RawTower*>::iterator iter = towers.begin();
253  iter != towers.end();
254  ++iter)
255  {
256  RawTower* tower = *iter;
257 
258  if (_strict)
259  {
260  assert(tower);
261  }
262  else if (!tower)
263  {
264  ++_errors;
265  continue;
266  }
267 
268  float energy = get_energy_contribution(tower, shower);
269  if (isnan(energy)) continue;
270  if (energy > best_energy)
271  {
272  best_tower = tower;
273  best_energy = energy;
274  }
275  }
276 
277  if (_do_cache) _cache_best_tower_from_primary_shower.insert(make_pair(shower, best_tower));
278 
279  return best_tower;
280 }
281 
282 std::set<RawTower*> CaloRawTowerEval::all_towers_from(PHG4Shower* shower)
283 {
285  {
286  ++_errors;
287  return std::set<RawTower*>();
288  }
289 
290  if (_strict)
291  {
292  assert(shower);
293  }
294  else if (!shower)
295  {
296  ++_errors;
297  return std::set<RawTower*>();
298  }
299 
300  if (!_trutheval.is_primary(shower)) return std::set<RawTower*>();
301 
302  if (_do_cache)
303  {
304  std::map<PHG4Shower*, std::set<RawTower*> >::iterator iter =
306  if (iter != _cache_all_towers_from_primary_shower.end())
307  {
308  return iter->second;
309  }
310  }
311 
312  std::set<RawTower*> towers;
313 
314  // loop over all the towers
315  for (RawTowerContainer::Iterator iter = _towers->getTowers().first;
316  iter != _towers->getTowers().second;
317  ++iter)
318  {
319  RawTower* tower = iter->second;
320 
321  std::set<PHG4Shower*> showers = all_truth_primary_showers(tower);
322  for (std::set<PHG4Shower*>::iterator jter = showers.begin();
323  jter != showers.end();
324  ++jter)
325  {
326  PHG4Shower* candidate = *jter;
327 
328  if (_strict)
329  {
330  assert(candidate);
331  }
332  else if (!candidate)
333  {
334  ++_errors;
335  continue;
336  }
337 
338  if (candidate->get_id() == shower->get_id())
339  {
340  towers.insert(tower);
341  }
342  }
343  }
344 
345  if (_do_cache) _cache_all_towers_from_primary_shower.insert(make_pair(shower, towers));
346 
347  return towers;
348 }
349 
351 {
353  {
354  ++_errors;
355  return NAN;
356  }
357 
358  if (_strict)
359  {
360  assert(tower);
361  assert(shower);
362  }
363  else if (!tower || !shower)
364  {
365  ++_errors;
366  return NAN;
367  }
368 
369  if (!_trutheval.is_primary(shower)) return NAN;
370 
371  if (_do_cache)
372  {
373  std::map<std::pair<RawTower*, PHG4Shower*>, float>::iterator iter =
374  _cache_get_energy_contribution_primary_shower.find(make_pair(tower, shower));
376  {
377  return iter->second;
378  }
379  }
380 
381  // loop over the tower shower entries
382  float energy = 0.0;
383  RawTower::ShowerConstRange range = tower->get_g4showers();
384  RawTower::ShowerConstIterator iter = tower->find_g4shower(shower->get_id());
385  if (iter != range.second)
386  {
387  energy = iter->second;
388  }
389 
390  if (_do_cache) _cache_get_energy_contribution_primary_shower.insert(make_pair(make_pair(tower, shower), energy));
391 
392  return energy;
393 }
394 
396 {
398  {
399  ++_errors;
400  return std::set<PHG4Particle*>();
401  }
402 
403  if (_strict)
404  {
405  assert(tower);
406  }
407  else if (!tower)
408  {
409  ++_errors;
410  return std::set<PHG4Particle*>();
411  }
412 
413  if (_do_cache)
414  {
415  std::map<RawTower*, std::set<PHG4Particle*> >::iterator iter =
417  if (iter != _cache_all_truth_primary_particles.end())
418  {
419  return iter->second;
420  }
421  }
422 
423  std::set<PHG4Particle*> truth_primaries;
424 
425  std::set<PHG4Shower*> showers = all_truth_primary_showers(tower);
426 
427  for (std::set<PHG4Shower*>::iterator iter = showers.begin();
428  iter != showers.end();
429  ++iter)
430  {
431  PHG4Shower* shower = *iter;
432  PHG4Particle* primary = get_truth_eval()->get_primary_particle(shower);
433 
434  if (_strict)
435  {
436  assert(primary);
437  }
438  else if (!primary)
439  {
440  ++_errors;
441  continue;
442  }
443 
444  truth_primaries.insert(primary);
445  }
446 
447  if (_do_cache) _cache_all_truth_primary_particles.insert(make_pair(tower, truth_primaries));
448 
449  return truth_primaries;
450 }
451 
453 {
455  {
456  ++_errors;
457  return nullptr;
458  }
459 
460  if (_strict)
461  {
462  assert(tower);
463  }
464  else if (!tower)
465  {
466  ++_errors;
467  return nullptr;
468  }
469 
470  if (_do_cache)
471  {
472  std::map<RawTower*, PHG4Particle*>::iterator iter =
475  {
476  return iter->second;
477  }
478  }
479 
480  PHG4Particle* max_primary = nullptr;
481  PHG4Shower* max_shower = max_truth_primary_shower_by_energy(tower);
482 
483  if (max_shower)
484  {
485  max_primary = get_truth_eval()->get_primary_particle(max_shower);
486  }
487 
488  if (_do_cache) _cache_max_truth_primary_particle_by_energy.insert(make_pair(tower, max_primary));
489 
490  return max_primary;
491 }
492 
493 std::set<RawTower*> CaloRawTowerEval::all_towers_from(PHG4Particle* primary)
494 {
496  {
497  ++_errors;
498  return std::set<RawTower*>();
499  }
500 
501  if (_strict)
502  {
503  assert(primary);
504  }
505  else if (!primary)
506  {
507  ++_errors;
508  return std::set<RawTower*>();
509  }
510 
511  if (!_trutheval.is_primary(primary)) return std::set<RawTower*>();
512 
513  // use primary map pointer
514  primary = get_truth_eval()->get_primary_particle(primary);
515 
516  if (_strict)
517  {
518  assert(primary);
519  }
520  else if (!primary)
521  {
522  ++_errors;
523  return std::set<RawTower*>();
524  }
525 
526  if (_do_cache)
527  {
528  std::map<PHG4Particle*, std::set<RawTower*> >::iterator iter =
530  if (iter != _cache_all_towers_from_primary_particle.end())
531  {
532  return iter->second;
533  }
534  }
535 
536  std::set<RawTower*> towers;
537 
538  PHG4Shower* shower = get_truth_eval()->get_primary_shower(primary);
539 
540  if (shower)
541  {
542  towers = all_towers_from(shower);
543  }
544 
545  if (_do_cache) _cache_all_towers_from_primary_particle.insert(make_pair(primary, towers));
546 
547  return towers;
548 }
549 
551 {
553  {
554  ++_errors;
555  return nullptr;
556  }
557 
558  if (_strict)
559  {
560  assert(primary);
561  }
562  else if (!primary)
563  {
564  ++_errors;
565  return nullptr;
566  }
567 
568  if (!_trutheval.is_primary(primary)) return nullptr;
569 
570  primary = get_truth_eval()->get_primary_particle(primary);
571 
572  if (_strict)
573  {
574  assert(primary);
575  }
576  else if (!primary)
577  {
578  ++_errors;
579  return nullptr;
580  }
581 
582  if (_do_cache)
583  {
584  std::map<PHG4Particle*, RawTower*>::iterator iter =
586  if (iter != _cache_best_tower_from_primary_particle.end())
587  {
588  return iter->second;
589  }
590  }
591 
592  RawTower* best_tower = nullptr;
593  PHG4Shower* shower = get_truth_eval()->get_primary_shower(primary);
594  if (shower)
595  {
596  best_tower = best_tower_from(shower);
597  }
598 
599  if (_do_cache) _cache_best_tower_from_primary_particle.insert(make_pair(primary, best_tower));
600 
601  return best_tower;
602 }
603 
604 // overlap calculations
606 {
608  {
609  ++_errors;
610  return NAN;
611  }
612 
613  if (_strict)
614  {
615  assert(tower);
616  assert(primary);
617  }
618  else if (!tower || !primary)
619  {
620  ++_errors;
621  return NAN;
622  }
623 
624  if (!_trutheval.is_primary(primary)) return NAN;
625 
626  // reduce cache misses by using only pointer from PrimaryMap
627  primary = get_truth_eval()->get_primary_particle(primary);
628 
629  if (_strict)
630  {
631  assert(primary);
632  }
633  else if (!primary)
634  {
635  ++_errors;
636  return NAN;
637  }
638 
639  if (_do_cache)
640  {
641  std::map<std::pair<RawTower*, PHG4Particle*>, float>::iterator iter =
642  _cache_get_energy_contribution_primary_particle.find(make_pair(tower, primary));
644  {
645  return iter->second;
646  }
647  }
648 
649  float energy = 0.0;
650 
651  PHG4Shower* shower = get_truth_eval()->get_primary_shower(primary);
652 
653  if (shower)
654  {
655  energy = get_energy_contribution(tower, shower);
656  }
657 
658  if (_do_cache) _cache_get_energy_contribution_primary_particle.insert(make_pair(make_pair(tower, primary), energy));
659 
660  return energy;
661 }
662 
664 {
665  if (!get_truth_eval()->has_full_node_pointers()) return false;
666 
667  if (_strict)
668  assert(_towers);
669  else if (!_towers)
670  return false;
671 
672  if (_strict)
673  assert(_g4cells);
674  else if (!_g4cells)
675  return false;
676 
677  if (_strict)
678  assert(_g4hits);
679  else if (!_g4hits)
680  return false;
681 
682  if (_strict)
683  assert(_truthinfo);
684  else if (!_truthinfo)
685  return false;
686 
687  return true;
688 }
689 
690 std::set<PHG4Hit*> CaloRawTowerEval::all_truth_hits(RawTower* tower)
691 {
692  if (!has_full_node_pointers())
693  {
694  ++_errors;
695  return std::set<PHG4Hit*>();
696  }
697 
698  if (_strict)
699  {
700  assert(tower);
701  }
702  else if (!tower)
703  {
704  ++_errors;
705  return std::set<PHG4Hit*>();
706  }
707 
708  if (_do_cache)
709  {
710  std::map<RawTower*, std::set<PHG4Hit*> >::iterator iter =
711  _cache_all_truth_hits.find(tower);
712  if (iter != _cache_all_truth_hits.end())
713  {
714  return iter->second;
715  }
716  }
717 
718  std::set<PHG4Hit*> truth_hits;
719  // loop over all the towered cells
720  RawTower::CellConstRange cell_range = tower->get_g4cells();
721  for (RawTower::CellConstIterator cell_iter = cell_range.first;
722  cell_iter != cell_range.second; ++cell_iter)
723  {
724  unsigned int cell_id = cell_iter->first;
725  PHG4Cell* cell = _g4cells->findCell(cell_id);
726 
727  if (_strict)
728  {
729  assert(cell);
730  }
731  else if (!cell)
732  {
733  ++_errors;
734  continue;
735  }
736 
737  // loop over all the g4hits in this cell
738  for (PHG4Cell::EdepConstIterator hit_iter = cell->get_g4hits().first;
739  hit_iter != cell->get_g4hits().second;
740  ++hit_iter)
741  {
742  PHG4Hit* g4hit = _g4hits->findHit(hit_iter->first);
743 
744  if (_strict)
745  {
746  assert(g4hit);
747  }
748  else if (!g4hit)
749  {
750  ++_errors;
751  continue;
752  }
753 
754  // fill output set
755  truth_hits.insert(g4hit);
756  }
757  }
758 
759  if (_do_cache) _cache_all_truth_hits.insert(make_pair(tower, truth_hits));
760 
761  return truth_hits;
762 }
763 
765 {
766  // need things off of the DST...
767  std::string towername = "TOWER_CALIB_" + _caloname;
768  _towers = findNode::getClass<RawTowerContainer>(topNode, towername.c_str());
769 
770  std::string cellname = "G4CELL_" + _caloname;
771  _g4cells = findNode::getClass<PHG4CellContainer>(topNode, cellname.c_str());
772 
773  std::string hitname = "G4HIT_" + _caloname;
774  _g4hits = findNode::getClass<PHG4HitContainer>(topNode, hitname.c_str());
775 
776  _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
777 
778  return;
779 }