ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CaloRawClusterEval.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CaloRawClusterEval.cc
1 
2 #include "CaloRawClusterEval.h"
3 #include "CaloTruthEval.h"
4 
5 #include <calobase/RawCluster.h>
6 #include <calobase/RawClusterContainer.h>
7 #include <calobase/RawTowerContainer.h>
8 
9 #include <phool/getClass.h>
10 
11 #include <cassert>
12 #include <cfloat>
13 #include <cmath>
14 #include <iostream>
15 #include <map>
16 #include <set>
17 #include <string>
18 
19 class RawTower;
20 
21 using namespace std;
22 
23 CaloRawClusterEval::CaloRawClusterEval(PHCompositeNode* topNode, const std::string& caloname)
24  : _caloname(caloname)
25  , _towereval(topNode, caloname)
26  , _clusters(nullptr)
27  , _towers(nullptr)
28  , _strict(false)
29  , _verbosity(1)
30  , _errors(0)
31  , _do_cache(true)
32  , _cache_all_truth_primary_showers()
33  , _cache_max_truth_primary_shower_by_energy()
34  , _cache_all_clusters_from_primary_shower()
35  , _cache_best_cluster_from_primary_shower()
36  , _cache_get_energy_contribution_primary_shower()
37  , _cache_all_truth_primary_particles()
38  , _cache_max_truth_primary_particle_by_energy()
39  , _cache_all_clusters_from_primary_particle()
40  , _cache_best_cluster_from_primary_particle()
41  , _cache_get_energy_contribution_primary_particle()
42  , _cache_all_truth_hits()
43 {
44  get_node_pointers(topNode);
45 }
46 
48 {
49  if (_verbosity > 0)
50  {
51  if ((_errors > 0) || (_verbosity > 1))
52  {
53  cout << "CaloRawClusterEval::~CaloRawClusterEval() - Error Count: " << _errors << endl;
54  }
55  }
56 }
57 
59 {
65 
71 
72  _cache_all_truth_hits.clear();
73 
74  _towereval.next_event(topNode);
75 
76  get_node_pointers(topNode);
77 }
78 
80 {
81  if (!get_rawtower_eval()->has_reduced_node_pointers()) return false;
82 
83  if (_strict)
84  assert(_clusters);
85  else if (!_clusters)
86  return false;
87 
88  if (_strict)
89  assert(_towers);
90  else if (!_towers)
91  return false;
92 
93  return true;
94 }
95 
97 {
99  {
100  ++_errors;
101  return std::set<PHG4Shower*>();
102  }
103 
104  if (_strict)
105  {
106  assert(cluster);
107  }
108  else if (!cluster)
109  {
110  ++_errors;
111  return std::set<PHG4Shower*>();
112  }
113 
114  if (_do_cache)
115  {
116  std::map<RawCluster*, std::set<PHG4Shower*> >::iterator iter =
117  _cache_all_truth_primary_showers.find(cluster);
118  if (iter != _cache_all_truth_primary_showers.end())
119  {
120  return iter->second;
121  }
122  }
123 
124  std::set<PHG4Shower*> truth_primary_showers;
125 
126  // loop over all the clustered towers
127  RawCluster::TowerConstRange begin_end = cluster->get_towers();
128  for (RawCluster::TowerConstIterator iter = begin_end.first;
129  iter != begin_end.second;
130  ++iter)
131  {
132  RawTower* tower = _towers->getTower(iter->first);
133 
134  if (_strict)
135  {
136  assert(tower);
137  }
138  else if (!tower)
139  {
140  ++_errors;
141  continue;
142  }
143 
144  std::set<PHG4Shower*> new_primary_showers = _towereval.all_truth_primary_showers(tower);
145 
146  for (std::set<PHG4Shower*>::iterator iter = new_primary_showers.begin();
147  iter != new_primary_showers.end();
148  ++iter)
149  {
150  PHG4Shower* shower = *iter;
151 
152  if (_strict)
153  {
154  assert(shower);
155  }
156  else if (!shower)
157  {
158  ++_errors;
159  continue;
160  }
161 
162  truth_primary_showers.insert(shower);
163  }
164  }
165 
166  if (_do_cache) _cache_all_truth_primary_showers.insert(make_pair(cluster, truth_primary_showers));
167 
168  return truth_primary_showers;
169 }
170 
172 {
174  {
175  ++_errors;
176  return nullptr;
177  }
178 
179  if (_strict)
180  {
181  assert(cluster);
182  }
183  else if (!cluster)
184  {
185  ++_errors;
186  return nullptr;
187  }
188 
189  if (_do_cache)
190  {
191  std::map<RawCluster*, PHG4Shower*>::iterator iter =
194  {
195  return iter->second;
196  }
197  }
198 
199  // loop over all primaries associated with this cluster and
200  // get the energy contribution for each one, record the max
201  PHG4Shower* max_primary = nullptr;
202  float max_e = FLT_MAX * -1.0;
203  std::set<PHG4Shower*> primary_showers = all_truth_primary_showers(cluster);
204  for (std::set<PHG4Shower*>::iterator iter = primary_showers.begin();
205  iter != primary_showers.end();
206  ++iter)
207  {
208  PHG4Shower* primary = *iter;
209 
210  if (_strict)
211  {
212  assert(primary);
213  }
214  else if (!primary)
215  {
216  ++_errors;
217  continue;
218  }
219 
220  float e = get_energy_contribution(cluster, primary);
221  if (isnan(e)) continue;
222  if (e > max_e)
223  {
224  max_e = e;
225  max_primary = primary;
226  }
227  }
228 
229  if (_do_cache) _cache_max_truth_primary_shower_by_energy.insert(make_pair(cluster, max_primary));
230 
231  return max_primary;
232 }
233 
234 std::set<RawCluster*> CaloRawClusterEval::all_clusters_from(PHG4Shower* primary)
235 {
237  {
238  ++_errors;
239  return std::set<RawCluster*>();
240  }
241 
242  if (_strict)
243  {
244  assert(primary);
245  }
246  else if (!primary)
247  {
248  ++_errors;
249  return std::set<RawCluster*>();
250  }
251 
252  if (!get_truth_eval()->is_primary(primary)) return std::set<RawCluster*>();
253 
254  primary = get_truth_eval()->get_primary_shower(primary);
255 
256  if (_strict)
257  {
258  assert(primary);
259  }
260  else if (!primary)
261  {
262  ++_errors;
263  return std::set<RawCluster*>();
264  }
265 
266  if (_do_cache)
267  {
268  std::map<PHG4Shower*, std::set<RawCluster*> >::iterator iter =
270  if (iter != _cache_all_clusters_from_primary_shower.end())
271  {
272  return iter->second;
273  }
274  }
275 
276  std::set<RawCluster*> clusters;
277 
278  // loop over all the clusters
280  iter != _clusters->getClusters().second;
281  ++iter)
282  {
283  RawCluster* cluster = iter->second;
284 
285  std::set<PHG4Shower*> primary_showers = all_truth_primary_showers(cluster);
286  for (std::set<PHG4Shower*>::iterator jter = primary_showers.begin();
287  jter != primary_showers.end();
288  ++jter)
289  {
290  PHG4Shower* candidate = *jter;
291 
292  if (_strict)
293  {
294  assert(candidate);
295  }
296  else if (!candidate)
297  {
298  ++_errors;
299  continue;
300  }
301 
302  if (get_truth_eval()->are_same_shower(candidate, primary))
303  {
304  clusters.insert(cluster);
305  }
306  }
307  }
308 
309  if (_do_cache) _cache_all_clusters_from_primary_shower.insert(make_pair(primary, clusters));
310 
311  return clusters;
312 }
313 
315 {
317  {
318  ++_errors;
319  return nullptr;
320  }
321 
322  if (_strict)
323  {
324  assert(primary);
325  }
326  else if (!primary)
327  {
328  ++_errors;
329  return nullptr;
330  }
331 
332  if (!get_truth_eval()->is_primary(primary)) return nullptr;
333 
334  primary = get_truth_eval()->get_primary_shower(primary);
335 
336  if (_strict)
337  {
338  assert(primary);
339  }
340  else if (!primary)
341  {
342  ++_errors;
343  return nullptr;
344  }
345 
346  if (_do_cache)
347  {
348  std::map<PHG4Shower*, RawCluster*>::iterator iter =
350  if (iter != _cache_best_cluster_from_primary_shower.end())
351  {
352  return iter->second;
353  }
354  }
355 
356  RawCluster* best_cluster = nullptr;
357  float best_energy = FLT_MAX * -1.0;
358  std::set<RawCluster*> clusters = all_clusters_from(primary);
359  for (std::set<RawCluster*>::iterator iter = clusters.begin();
360  iter != clusters.end();
361  ++iter)
362  {
363  RawCluster* cluster = *iter;
364 
365  if (_strict)
366  {
367  assert(cluster);
368  }
369  else if (!cluster)
370  {
371  ++_errors;
372  continue;
373  }
374 
375  float energy = get_energy_contribution(cluster, primary);
376  if (isnan(energy)) continue;
377  if (energy > best_energy)
378  {
379  best_cluster = cluster;
380  best_energy = energy;
381  }
382  }
383 
384  if (_do_cache) _cache_best_cluster_from_primary_shower.insert(make_pair(primary, best_cluster));
385 
386  return best_cluster;
387 }
388 
390 {
392  {
393  ++_errors;
394  return NAN;
395  }
396 
397  if (_strict)
398  {
399  assert(cluster);
400  assert(primary);
401  }
402  else if (!cluster || !primary)
403  {
404  ++_errors;
405  return NAN;
406  }
407 
408  if (!get_truth_eval()->is_primary(primary)) return NAN;
409 
410  // reduce cache misses by using only pointer from PrimaryMap
411  primary = get_truth_eval()->get_primary_shower(primary);
412 
413  if (_strict)
414  {
415  assert(primary);
416  }
417  else if (!primary)
418  {
419  ++_errors;
420  return NAN;
421  }
422 
423  if (_do_cache)
424  {
425  std::map<std::pair<RawCluster*, PHG4Shower*>, float>::iterator iter =
426  _cache_get_energy_contribution_primary_shower.find(make_pair(cluster, primary));
428  {
429  return iter->second;
430  }
431  }
432 
433  float energy = 0.0;
434 
435  // loop over all the clustered towers
436  RawCluster::TowerConstRange begin_end = cluster->get_towers();
437  for (RawCluster::TowerConstIterator iter = begin_end.first;
438  iter != begin_end.second;
439  ++iter)
440  {
441  RawTower* tower = _towers->getTower(iter->first);
442 
443  if (_strict)
444  {
445  assert(tower);
446  }
447  else if (!tower)
448  {
449  ++_errors;
450  continue;
451  }
452 
453  float edep = get_rawtower_eval()->get_energy_contribution(tower, primary);
454  if (!isnan(edep)) energy += edep;
455  }
456 
457  if (_do_cache) _cache_get_energy_contribution_primary_shower.insert(make_pair(make_pair(cluster, primary), energy));
458 
459  return energy;
460 }
461 
463 {
465  {
466  ++_errors;
467  return std::set<PHG4Particle*>();
468  }
469 
470  if (_strict)
471  {
472  assert(cluster);
473  }
474  else if (!cluster)
475  {
476  ++_errors;
477  return std::set<PHG4Particle*>();
478  }
479 
480  if (_do_cache)
481  {
482  std::map<RawCluster*, std::set<PHG4Particle*> >::iterator iter =
484  if (iter != _cache_all_truth_primary_particles.end())
485  {
486  return iter->second;
487  }
488  }
489 
490  std::set<PHG4Particle*> truth_primary_particles;
491 
492  std::set<PHG4Shower*> primary_showers = all_truth_primary_showers(cluster);
493 
494  for (std::set<PHG4Shower*>::iterator iter = primary_showers.begin();
495  iter != primary_showers.end();
496  ++iter)
497  {
498  PHG4Shower* shower = *iter;
499 
500  if (_strict)
501  assert(shower);
502  else if (!shower)
503  {
504  ++_errors;
505  continue;
506  }
507 
509 
510  if (_strict)
511  assert(particle);
512  else if (!particle)
513  {
514  ++_errors;
515  continue;
516  }
517 
518  truth_primary_particles.insert(particle);
519  }
520 
521  if (_do_cache) _cache_all_truth_primary_particles.insert(make_pair(cluster, truth_primary_particles));
522 
523  return truth_primary_particles;
524 }
525 
527 {
529  {
530  ++_errors;
531  return nullptr;
532  }
533 
534  if (_strict)
535  {
536  assert(cluster);
537  }
538  else if (!cluster)
539  {
540  ++_errors;
541  return nullptr;
542  }
543 
544  if (_do_cache)
545  {
546  std::map<RawCluster*, PHG4Particle*>::iterator iter =
549  {
550  return iter->second;
551  }
552  }
553 
554  PHG4Particle* max_primary = nullptr;
555  PHG4Shower* max_shower = max_truth_primary_shower_by_energy(cluster);
556 
557  if (max_shower)
558  {
559  max_primary = get_truth_eval()->get_primary_particle(max_shower);
560  }
561 
562  if (_do_cache) _cache_max_truth_primary_particle_by_energy.insert(make_pair(cluster, max_primary));
563 
564  return max_primary;
565 }
566 
568 {
570  {
571  ++_errors;
572  return std::set<RawCluster*>();
573  }
574 
575  if (_strict)
576  {
577  assert(primary);
578  }
579  else if (!primary)
580  {
581  ++_errors;
582  return std::set<RawCluster*>();
583  }
584 
585  if (!get_truth_eval()->is_primary(primary)) return std::set<RawCluster*>();
586 
587  primary = get_truth_eval()->get_primary_particle(primary);
588 
589  if (_strict)
590  {
591  assert(primary);
592  }
593  else if (!primary)
594  {
595  ++_errors;
596  return std::set<RawCluster*>();
597  }
598 
599  if (_do_cache)
600  {
601  std::map<PHG4Particle*, std::set<RawCluster*> >::iterator iter =
604  {
605  return iter->second;
606  }
607  }
608 
609  std::set<RawCluster*> clusters;
610 
611  PHG4Shower* shower = get_truth_eval()->get_primary_shower(primary);
612 
613  if (shower) clusters = all_clusters_from(shower);
614 
615  if (_do_cache) _cache_all_clusters_from_primary_particle.insert(make_pair(primary, clusters));
616 
617  return clusters;
618 }
619 
621 {
623  {
624  ++_errors;
625  return nullptr;
626  }
627 
628  if (_strict)
629  {
630  assert(primary);
631  }
632  else if (!primary)
633  {
634  ++_errors;
635  return nullptr;
636  }
637 
638  if (!get_truth_eval()->is_primary(primary)) return nullptr;
639 
640  primary = get_truth_eval()->get_primary_particle(primary);
641 
642  if (_strict)
643  {
644  assert(primary);
645  }
646  else if (!primary)
647  {
648  ++_errors;
649  return nullptr;
650  }
651 
652  if (_do_cache)
653  {
654  std::map<PHG4Particle*, RawCluster*>::iterator iter =
657  {
658  return iter->second;
659  }
660  }
661 
662  RawCluster* best_cluster = nullptr;
663 
664  PHG4Shower* shower = get_truth_eval()->get_primary_shower(primary);
665  if (shower) best_cluster = best_cluster_from(shower);
666 
667  if (_do_cache) _cache_best_cluster_from_primary_particle.insert(make_pair(primary, best_cluster));
668 
669  return best_cluster;
670 }
671 
673 {
675  {
676  ++_errors;
677  return NAN;
678  }
679 
680  if (_strict)
681  {
682  assert(cluster);
683  assert(primary);
684  }
685  else if (!cluster || !primary)
686  {
687  ++_errors;
688  return NAN;
689  }
690 
691  if (!get_truth_eval()->is_primary(primary)) return NAN;
692 
693  // reduce cache misses by using only pointer from PrimaryMap
694  primary = get_truth_eval()->get_primary_particle(primary);
695 
696  if (_strict)
697  {
698  assert(primary);
699  }
700  else if (!primary)
701  {
702  ++_errors;
703  return 0.;
704  }
705 
706  if (_do_cache)
707  {
708  std::map<std::pair<RawCluster*, PHG4Particle*>, float>::iterator iter =
709  _cache_get_energy_contribution_primary_particle.find(make_pair(cluster, primary));
711  {
712  return iter->second;
713  }
714  }
715 
716  float energy = 0.0;
717  PHG4Shower* shower = get_truth_eval()->get_primary_shower(primary);
718  if (shower)
719  {
720  float edep = get_energy_contribution(cluster, shower);
721  if (!isnan(edep)) energy += edep;
722  }
723 
724  if (_do_cache) _cache_get_energy_contribution_primary_particle.insert(make_pair(make_pair(cluster, primary), energy));
725 
726  return energy;
727 }
728 
730 {
731  if (!get_rawtower_eval()->has_full_node_pointers()) return false;
732 
733  if (_strict)
734  assert(_clusters);
735  else if (!_clusters)
736  return false;
737 
738  if (_strict)
739  assert(_towers);
740  else if (!_towers)
741  return false;
742 
743  return true;
744 }
745 
746 std::set<PHG4Hit*> CaloRawClusterEval::all_truth_hits(RawCluster* cluster)
747 {
748  if (!has_full_node_pointers())
749  {
750  ++_errors;
751  return std::set<PHG4Hit*>();
752  }
753 
754  if (_strict)
755  {
756  assert(cluster);
757  }
758  else if (!cluster)
759  {
760  ++_errors;
761  return std::set<PHG4Hit*>();
762  }
763 
764  if (_do_cache)
765  {
766  std::map<RawCluster*, std::set<PHG4Hit*> >::iterator iter =
767  _cache_all_truth_hits.find(cluster);
768  if (iter != _cache_all_truth_hits.end())
769  {
770  return iter->second;
771  }
772  }
773 
774  std::set<PHG4Hit*> truth_hits;
775 
776  // loop over all the clustered towers
777  RawCluster::TowerConstRange begin_end = cluster->get_towers();
778  for (RawCluster::TowerConstIterator iter = begin_end.first;
779  iter != begin_end.second;
780  ++iter)
781  {
782  RawTower* tower = _towers->getTower(iter->first);
783 
784  if (_strict)
785  {
786  assert(tower);
787  }
788  else if (!tower)
789  {
790  ++_errors;
791  continue;
792  }
793 
794  std::set<PHG4Hit*> new_hits = get_rawtower_eval()->all_truth_hits(tower);
795 
796  for (std::set<PHG4Hit*>::iterator iter = new_hits.begin();
797  iter != new_hits.end();
798  ++iter)
799  {
800  PHG4Hit* g4hit = *iter;
801 
802  if (_strict)
803  {
804  assert(g4hit);
805  }
806  else if (!g4hit)
807  {
808  ++_errors;
809  continue;
810  }
811 
812  truth_hits.insert(g4hit);
813  }
814  }
815 
816  if (_do_cache) _cache_all_truth_hits.insert(make_pair(cluster, truth_hits));
817 
818  return truth_hits;
819 }
820 
822 {
823  // need things off of the DST...
824  std::string nodename = "CLUSTER_" + _caloname;
825  _clusters = findNode::getClass<RawClusterContainer>(topNode, nodename.c_str());
826 
827  std::string towername = "TOWER_CALIB_" + _caloname;
828  _towers = findNode::getClass<RawTowerContainer>(topNode, towername.c_str());
829 
830  return;
831 }