ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SvtxHitEval.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SvtxHitEval.cc
1 #include "SvtxHitEval.h"
2 
3 #include <trackbase/TrkrDefs.h>
4 #include <trackbase/TrkrHitSet.h>
8 
9 #include <g4main/PHG4Hit.h>
11 #include <g4main/PHG4HitDefs.h>
12 #include <g4main/PHG4Particle.h>
14 
15 #include <phool/getClass.h>
16 #include <phool/phool.h>
17 
18 #include <cassert>
19 #include <cfloat>
20 #include <cmath>
21 #include <iostream> // for operator<<, endl, basic_...
22 #include <map>
23 #include <set>
24 
25 class TrkrHit;
26 
27 using namespace std;
28 
30  : _trutheval(topNode)
31  , _hitmap(nullptr)
32  // , _g4cells_svtx(nullptr)
33  // , _g4cells_tracker(nullptr)
34  //, _g4cells_maps(nullptr)
35  , _g4hits_tpc(nullptr)
36  , _g4hits_intt(nullptr)
37  , _g4hits_mvtx(nullptr)
38  , _g4hits_mms(nullptr)
39  , _truthinfo(nullptr)
40  , _strict(false)
41  , _verbosity(0)
42  , _errors(0)
43  , _do_cache(true)
44  , _cache_all_truth_hits()
45  , _cache_max_truth_hit_by_energy()
46  , _cache_all_truth_particles()
47  , _cache_max_truth_particle_by_energy()
48  , _cache_all_hits_from_particle()
49  , _cache_all_hits_from_g4hit()
50  , _cache_best_hit_from_g4hit()
51  , _cache_get_energy_contribution_g4particle()
52  , _cache_get_energy_contribution_g4hit()
53 {
54  get_node_pointers(topNode);
55 }
56 
58 {
59  if (_verbosity > 0)
60  {
61  if ((_errors > 0) || (_verbosity > 1))
62  {
63  cout << "SvtxHitEval::~SvtxHitEval() - Error Count: " << _errors << endl;
64  }
65  }
66 }
67 
69 {
70  _cache_all_truth_hits.clear();
79 
80  _trutheval.next_event(topNode);
81 
82  get_node_pointers(topNode);
83 }
84 
85 /*
86 PHG4Cell* SvtxHitEval::get_cell(SvtxHit* hit)
87 {
88  if (!has_node_pointers())
89  {
90  ++_errors;
91  cout << PHWHERE << " nerr: " << _errors << endl;
92  return nullptr;
93  }
94 
95  if (_strict)
96  {
97  assert(hit);
98  }
99  else if (!hit)
100  {
101  ++_errors;
102  cout << PHWHERE << " nerr: " << _errors << endl;
103  return nullptr;
104  }
105 
106  // hop from reco hit to g4cell
107  PHG4Cell* cell = nullptr;
108  if (_g4cells_svtx) cell = _g4cells_svtx->findCell(hit->get_cellid());
109  if (!cell && _g4cells_tracker) cell = _g4cells_tracker->findCell(hit->get_cellid());
110  if (!cell && _g4cells_maps) cell = _g4cells_maps->findCell(hit->get_cellid());
111 
112  // only noise hits (cellid left at default value) should not trace
113  if ((_strict) && (hit->get_cellid() != 0xFFFFFFFF))
114  {
115  assert(cell);
116  }
117  else if (!cell)
118  {
119  ++_errors;
120  cout << PHWHERE << " nerr: " << _errors << endl;
121  }
122 
123  return cell;
124 }
125 */
126 
128 {
129  if (!has_node_pointers())
130  {
131  ++_errors;
132  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
133  return std::set<PHG4Hit*>();
134  }
135 
136  if (_strict)
137  {
138  assert(hit_key);
139  }
140  else if (!hit_key)
141  {
142  ++_errors;
143  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
144  return std::set<PHG4Hit*>();
145  }
146 
147  if (_do_cache)
148  {
149  std::map<TrkrDefs::hitkey, std::set<PHG4Hit*> >::iterator iter =
150  _cache_all_truth_hits.find(hit_key);
151  if (iter != _cache_all_truth_hits.end())
152  {
153  return iter->second;
154  }
155  }
156 
157  std::set<PHG4Hit*> truth_hits;
158 
159  /*
160  // hop from reco hit to g4cell
161  PHG4Cell* cell = nullptr;
162  if (_g4cells_svtx) cell = _g4cells_svtx->findCell(hit->get_cellid());
163  if (!cell && _g4cells_tracker) cell = _g4cells_tracker->findCell(hit->get_cellid());
164  if (!cell && _g4cells_maps) cell = _g4cells_maps->findCell(hit->get_cellid());
165 
166  // only noise hits (cellid left at default value) should not trace
167  if ((_strict) && (hit->get_cellid() != 0xFFFFFFFF))
168  assert(cell);
169  else if (!cell)
170  {
171  ++_errors;
172  cout << PHWHERE << " nerr: " << _errors << endl;
173  return truth_hits;
174  }
175 
176  //cout << "Eval: hitid " << hit->get_id() << " cellid " << cell->get_cellid() << endl;
177  // loop over all the g4hits in this cell
178  for (PHG4Cell::EdepConstIterator g4iter = cell->get_g4hits().first;
179  g4iter != cell->get_g4hits().second;
180  ++g4iter)
181  {
182  //cout << " Looking for hit " << g4iter->first << " in layer " << cell->get_layer() << " with edep " << g4iter->second << endl;
183  PHG4Hit* g4hit = nullptr;
184  if (_g4hits_svtx) g4hit = _g4hits_svtx->findHit(g4iter->first);
185  if (!g4hit && _g4hits_tracker) g4hit = _g4hits_tracker->findHit(g4iter->first);
186  if (!g4hit && _g4hits_maps) g4hit = _g4hits_maps->findHit(g4iter->first);
187  if (!g4hit) cout << " Failed to find g4hit " << g4iter->first << " with edep " << g4iter->second << endl;
188  if (_strict)
189  assert(g4hit);
190  else if (!g4hit)
191  {
192  ++_errors;
193  cout << PHWHERE << " nerr: " << _errors << endl;
194  continue;
195  }
196  */
197 
198  // get all of the g4hits for this hit_key
199  // have to start with all hitsets, unfortunately
201  for(TrkrHitSetContainer::ConstIterator iter = all_hitsets.first; iter != all_hitsets.second; ++iter)
202  {
203  TrkrDefs::hitsetkey hitset_key = iter->first;
204  unsigned int trkrid = TrkrDefs::getTrkrId(hitset_key);
205  TrkrHitSet *hitset = iter->second;
206 
207  // does this hitset contain our hitkey?
208  TrkrHit *hit = nullptr;
209  hit = hitset->getHit(hit_key);
210  if(hit)
211  {
212  // get g4hits for this hit
213 
214  std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> > temp_map;
215  _hit_truth_map->getG4Hits(hitset_key, hit_key, temp_map); // returns pairs (hitsetkey, std::pair(hitkey, g4hitkey)) for this hitkey only
216  for(std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> >::iterator htiter = temp_map.begin();
217  htiter != temp_map.end(); ++htiter)
218  {
219  // extract the g4 hit key here and add the g4hit to the set
220  PHG4HitDefs::keytype g4hitkey = htiter->second.second;
221  //cout << " hitkey " << hitkey << " g4hitkey " << g4hitkey << endl;
222  PHG4Hit * g4hit = nullptr;
223  switch( trkrid )
224  {
225  case TrkrDefs::tpcId: g4hit = _g4hits_tpc->findHit(g4hitkey); break;
226  case TrkrDefs::inttId: g4hit = _g4hits_intt->findHit(g4hitkey); break;
227  case TrkrDefs::mvtxId: g4hit = _g4hits_mvtx->findHit(g4hitkey); break;
228  case TrkrDefs::micromegasId: g4hit = _g4hits_mms->findHit(g4hitkey); break;
229  default: break;
230  }
231  // fill output set
232  if( g4hit ) truth_hits.insert(g4hit);
233  }
234  }
235  }
236 
237  if (_do_cache) _cache_all_truth_hits.insert(make_pair(hit_key, truth_hits));
238 
239  return truth_hits;
240 }
241 
243 {
244  if (!has_node_pointers())
245  {
246  ++_errors;
247  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
248  return nullptr;
249  }
250 
251  if (_strict)
252  {
253  assert(hit_key);
254  }
255  else if (!hit_key)
256  {
257  ++_errors;
258  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
259  return nullptr;
260  }
261 
262  if (_do_cache)
263  {
264  std::map<TrkrDefs::hitkey, PHG4Hit*>::iterator iter =
265  _cache_max_truth_hit_by_energy.find(hit_key);
266  if (iter != _cache_max_truth_hit_by_energy.end())
267  {
268  return iter->second;
269  }
270  }
271 
272  std::set<PHG4Hit*> hits = all_truth_hits(hit_key);
273  PHG4Hit* max_hit = nullptr;
274  float max_e = FLT_MAX * -1.0;
275  for (std::set<PHG4Hit*>::iterator iter = hits.begin();
276  iter != hits.end();
277  ++iter)
278  {
279  PHG4Hit* hit = *iter;
280  if (hit->get_edep() > max_e)
281  {
282  max_e = hit->get_edep();
283  max_hit = hit;
284  }
285  }
286 
287  if (_do_cache) _cache_max_truth_hit_by_energy.insert(make_pair(hit_key, max_hit));
288 
289  return max_hit;
290 }
291 
292 std::set<PHG4Particle*> SvtxHitEval::all_truth_particles(TrkrDefs::hitkey hit_key)
293 {
294  if (!has_node_pointers())
295  {
296  ++_errors;
297  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
298  return std::set<PHG4Particle*>();
299  }
300 
301  if (_strict)
302  {
303  assert(hit_key);
304  }
305  else if (!hit_key)
306  {
307  ++_errors;
308  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
309  return std::set<PHG4Particle*>();
310  }
311 
312  if (_do_cache)
313  {
314  std::map<TrkrDefs::hitkey, std::set<PHG4Particle*> >::iterator iter =
315  _cache_all_truth_particles.find(hit_key);
316  if (iter != _cache_all_truth_particles.end())
317  {
318  return iter->second;
319  }
320  }
321 
322  std::set<PHG4Particle*> truth_particles;
323 
324  std::set<PHG4Hit*> g4hits = all_truth_hits(hit_key);
325 
326  for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
327  iter != g4hits.end();
328  ++iter)
329  {
330  PHG4Hit* g4hit = *iter;
332 
333  if (_strict)
334  assert(particle);
335  else if (!particle)
336  {
337  ++_errors;
338  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
339  continue;
340  }
341 
342  truth_particles.insert(particle);
343  }
344 
345  if (_do_cache) _cache_all_truth_particles.insert(make_pair(hit_key, truth_particles));
346 
347  return truth_particles;
348 }
349 
351 {
352  if (!has_node_pointers())
353  {
354  ++_errors;
355  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
356  return nullptr;
357  }
358 
359  if (_strict)
360  {
361  assert(hit_key);
362  }
363  else if (!hit_key)
364  {
365  ++_errors;
366  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
367  return nullptr;
368  }
369 
370  if (_do_cache)
371  {
372  std::map<TrkrDefs::hitkey, PHG4Particle*>::iterator iter =
374  if (iter != _cache_max_truth_particle_by_energy.end())
375  {
376  return iter->second;
377  }
378  }
379 
380  // loop over all particles associated with this hit and
381  // get the energy contribution for each one, record the max
382  PHG4Particle* max_particle = nullptr;
383  float max_e = FLT_MAX * -1.0;
384  std::set<PHG4Particle*> particles = all_truth_particles(hit_key);
385  for (std::set<PHG4Particle*>::iterator iter = particles.begin();
386  iter != particles.end();
387  ++iter)
388  {
389  PHG4Particle* particle = *iter;
390  float e = get_energy_contribution(hit_key, particle);
391  if (e > max_e)
392  {
393  max_e = e;
394  max_particle = particle;
395  }
396  }
397 
398  if (_do_cache) _cache_max_truth_particle_by_energy.insert(make_pair(hit_key, max_particle));
399 
400  return max_particle;
401 }
402 
403 std::set<TrkrDefs::hitkey> SvtxHitEval::all_hits_from(PHG4Particle* g4particle)
404 {
405  if (!has_node_pointers())
406  {
407  ++_errors;
408  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
409  return std::set<TrkrDefs::hitkey>();
410  }
411 
412  if (_strict)
413  {
414  assert(g4particle);
415  }
416  else if (!g4particle)
417  {
418  ++_errors;
419  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
420  return std::set<TrkrDefs::hitkey>();
421  }
422 
423  if (_do_cache)
424  {
425  std::map<PHG4Particle*, std::set<TrkrDefs::hitkey> >::iterator iter =
426  _cache_all_hits_from_particle.find(g4particle);
427  if (iter != _cache_all_hits_from_particle.end())
428  {
429  return iter->second;
430  }
431  }
432 
433  std::set<TrkrDefs::hitkey> hits;
434 
435  // loop over all the hits
437  for (TrkrHitSetContainer::ConstIterator iter = all_hitsets.first;
438  iter != all_hitsets.second;
439  ++iter)
440  {
441  TrkrHitSet::ConstRange range = iter->second->getHits();
442  for(TrkrHitSet::ConstIterator hitr = range.first; hitr != range.second; ++hitr)
443  {
444  TrkrDefs::hitkey hit_key = hitr->first;
445 
446  // loop over all truth hits connected to this hit
447  std::set<PHG4Hit*> g4hits = all_truth_hits(hit_key);
448  for (std::set<PHG4Hit*>::iterator jter = g4hits.begin();
449  jter != g4hits.end();
450  ++jter)
451  {
452  PHG4Hit* candidate = *jter;
454  if (g4particle->get_track_id() == particle->get_track_id())
455  {
456  hits.insert(hit_key);
457  }
458  }
459  }
460  }
461 
462  if (_do_cache) _cache_all_hits_from_particle.insert(make_pair(g4particle, hits));
463 
464  return hits;
465 }
466 
467 std::set<TrkrDefs::hitkey> SvtxHitEval::all_hits_from(PHG4Hit* g4hit)
468 {
469  if (!has_node_pointers())
470  {
471  ++_errors;
472  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
473  return std::set<TrkrDefs::hitkey>();
474  }
475 
476  if (_strict)
477  {
478  assert(g4hit);
479  }
480  else if (!g4hit)
481  {
482  ++_errors;
483  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
484  return std::set<TrkrDefs::hitkey>();
485  }
486 
487  if (_do_cache)
488  {
489  std::map<PHG4Hit*, std::set<TrkrDefs::hitkey> >::iterator iter =
490  _cache_all_hits_from_g4hit.find(g4hit);
491  if (iter != _cache_all_hits_from_g4hit.end())
492  {
493  return iter->second;
494  }
495  }
496 
497  std::set<TrkrDefs::hitkey> hits;
498 
499  unsigned int hit_layer = g4hit->get_layer();
500 
501  // loop over all the hits
503  for (TrkrHitSetContainer::ConstIterator iter = all_hitsets.first;
504  iter != all_hitsets.second;
505  ++iter)
506  {
507  TrkrHitSet::ConstRange range = iter->second->getHits();
508  for(TrkrHitSet::ConstIterator hitr = range.first; hitr != range.second; ++hitr)
509  {
510  TrkrDefs::hitkey hit_key = hitr->first;
511 
512  if (TrkrDefs::getLayer(hit_key) != hit_layer) continue;
513 
514  // loop over all truth hits connected to this hit
515  std::set<PHG4Hit*> g4hits = all_truth_hits(hit_key);
516  for (std::set<PHG4Hit*>::iterator jter = g4hits.begin();
517  jter != g4hits.end();
518  ++jter)
519  {
520  PHG4Hit* candidate = *jter;
521  if (candidate->get_hit_id() == g4hit->get_hit_id())
522  {
523  hits.insert(hit_key);
524  }
525  }
526  }
527  }
528 
529  if (_do_cache) _cache_all_hits_from_g4hit.insert(make_pair(g4hit, hits));
530 
531  return hits;
532 }
533 
535 {
536  if (!has_node_pointers())
537  {
538  ++_errors;
539  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
540  return 0;
541  }
542 
543  if (_strict)
544  {
545  assert(g4hit);
546  }
547  else if (!g4hit)
548  {
549  ++_errors;
550  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
551  return 0;
552  }
553 
554  if (_do_cache)
555  {
556  std::map<PHG4Hit*, TrkrDefs::hitkey>::iterator iter =
557  _cache_best_hit_from_g4hit.find(g4hit);
558  if (iter != _cache_best_hit_from_g4hit.end())
559  {
560  return iter->second;
561  }
562  }
563 
564  TrkrDefs::hitkey best_hit = 0;
565  float best_energy = 0.0;
566  std::set<TrkrDefs::hitkey> hits = all_hits_from(g4hit);
567  for (std::set<TrkrDefs::hitkey>::iterator iter = hits.begin();
568  iter != hits.end();
569  ++iter)
570  {
571  TrkrDefs::hitkey hit_key = *iter;
572  float energy = get_energy_contribution(hit_key, g4hit);
573  if (energy > best_energy)
574  {
575  best_hit = hit_key;
576  best_energy = energy;
577  }
578  }
579 
580  if (_do_cache) _cache_best_hit_from_g4hit.insert(make_pair(g4hit, best_hit));
581 
582  return best_hit;
583 }
584 
585 // overlap calculations
587 {
588  if (!has_node_pointers())
589  {
590  ++_errors;
591  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
592  return NAN;
593  }
594 
595  if (_strict)
596  {
597  assert(hit_key);
598  assert(particle);
599  }
600  else if (!hit_key || !particle)
601  {
602  ++_errors;
603  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
604  return NAN;
605  }
606 
607  if (_do_cache)
608  {
609  std::map<std::pair<TrkrDefs::hitkey, PHG4Particle*>, float>::iterator iter =
610  _cache_get_energy_contribution_g4particle.find(make_pair(hit_key, particle));
612  {
613  return iter->second;
614  }
615  }
616 
617  float energy = 0.0;
618  std::set<PHG4Hit*> g4hits = all_truth_hits(hit_key);
619  for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
620  iter != g4hits.end();
621  ++iter)
622  {
623  PHG4Hit* g4hit = *iter;
624  if (get_truth_eval()->is_g4hit_from_particle(g4hit, particle))
625  {
626  energy += g4hit->get_edep();
627  }
628  }
629 
630  if (_do_cache) _cache_get_energy_contribution_g4particle.insert(make_pair(make_pair(hit_key, particle), energy));
631 
632  return energy;
633 }
634 
636 {
637  if (!has_node_pointers())
638  {
639  ++_errors;
640  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
641  return NAN;
642  }
643 
644  if (_strict)
645  {
646  assert(hit_key);
647  assert(g4hit);
648  }
649  else if (!hit_key || !g4hit)
650  {
651  ++_errors;
652  if(_verbosity > 0) cout << PHWHERE << " nerr: " << _errors << endl;
653  return NAN;
654  }
655 
656  if (_do_cache)
657  {
658  std::map<std::pair<TrkrDefs::hitkey, PHG4Hit*>, float>::iterator iter =
659  _cache_get_energy_contribution_g4hit.find(make_pair(hit_key, g4hit));
660  if (iter != _cache_get_energy_contribution_g4hit.end())
661  {
662  return iter->second;
663  }
664  }
665 
666  // this is a fairly simple existance check right now, but might be more
667  // complex in the future, so this is here mostly as future-proofing.
668 
669  float energy = 0.0;
670  std::set<PHG4Hit*> g4hits = all_truth_hits(hit_key);
671  for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
672  iter != g4hits.end();
673  ++iter)
674  {
675  PHG4Hit* candidate = *iter;
676  if (candidate->get_hit_id() != g4hit->get_hit_id()) continue;
677  energy += candidate->get_edep();
678  }
679 
680  if (_do_cache) _cache_get_energy_contribution_g4hit.insert(make_pair(make_pair(hit_key, g4hit), energy));
681 
682  return energy;
683 }
684 
686 {
687  // need things off of the DST...
688  _hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
689 
690  _clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
691  if(!_clustermap)
692  _clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
693 
694  _hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode,"TRKR_HITTRUTHASSOC");
695 
696  // need things off of the DST...
697  _g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
698  _g4hits_intt = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_INTT");
699  _g4hits_mvtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
700  _g4hits_mms = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MICROMEGAS");
701 
702  _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
703 
704  return;
705 }
706 
708 {
709  if (_strict)
710  assert(_hitmap);
711  else if (!_hitmap)
712  {
713  return false;
714  }
715 
716  if (_strict)
718  else if (!_g4hits_mms && !_g4hits_tpc && !_g4hits_intt && !_g4hits_mvtx)
719  {
720  cout << "no hits" << endl;
721  return false;
722  }
723  if (_strict)
724  assert(_truthinfo);
725  else if (!_truthinfo)
726  {
727  cout << " no truth" << endl;
728  return false;
729  }
730 
731  return true;
732 }