ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SvtxClusterEval.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SvtxClusterEval.cc
1 #include "SvtxClusterEval.h"
2 
3 #include "SvtxHitEval.h"
4 #include "SvtxTruthEval.h"
5 
7 #include <trackbase/TrkrDefs.h>
11 #include <trackbase/TrkrHitSet.h>
14 
15 #include <g4main/PHG4Hit.h>
16 #include <g4main/PHG4Particle.h>
17 #include <g4main/PHG4VtxPoint.h>
19 #include <g4main/PHG4HitDefs.h>
21 
22 #include <phool/getClass.h>
23 #include <phool/PHTimer.h>
24 
25 #include <cassert>
26 #include <cfloat>
27 #include <cmath>
28 #include <iostream> // for operator<<, basic_ostream
29 #include <map>
30 #include <set>
31 #include <TVector3.h>
32 
33 using namespace std;
34 
36  : _hiteval(topNode)
37  , _clustermap(nullptr)
38  , _hitsets(nullptr)
39  , _truthinfo(nullptr)
40  , _strict(false)
41  , _verbosity(0)
42  , _errors(0)
43  , _do_cache(true)
44  , _cache_all_truth_hits()
45  , _cache_all_truth_clusters()
46  , _cache_max_truth_hit_by_energy()
47  , _cache_max_truth_cluster_by_energy()
48  , _cache_all_truth_particles()
49  , _cache_max_truth_particle_by_energy()
50  , _cache_max_truth_particle_by_cluster_energy()
51  , _cache_all_clusters_from_particle()
52  , _cache_all_clusters_from_g4hit()
53  , _cache_best_cluster_from_g4hit()
54  , _cache_get_energy_contribution_g4particle()
55  , _cache_get_energy_contribution_g4hit()
56  , _cache_reco_cluster_from_truth_cluster()
57 {
58  get_node_pointers(topNode);
59 }
60 
62 {
63  if (_verbosity > 0)
64  {
65  if ((_errors > 0) || (_verbosity > 1))
66  {
67  cout << "SvtxClusterEval::~SvtxClusterEval() - Error Count: " << _errors << endl;
68  }
69  }
70 }
71 
73 {
74  _cache_all_truth_hits.clear();
87  _clusters_per_layer.clear();
88  // _g4hits_per_layer.clear();
89  _hiteval.next_event(topNode);
90 
91  get_node_pointers(topNode);
92 }
93 
94 std::set<std::shared_ptr<TrkrCluster> > SvtxClusterEval::all_truth_clusters(TrkrDefs::cluskey cluster_key)
95 {
96  if (_do_cache)
97  {
98  std::map<TrkrDefs::cluskey, std::set<std::shared_ptr<TrkrCluster> > >::iterator iter =
99  _cache_all_truth_clusters.find(cluster_key);
100  if (iter != _cache_all_truth_clusters.end())
101  {
102  return iter->second;
103  }
104  }
105 
106  std::set<std::shared_ptr<TrkrCluster> > truth_clusters;
107 
108  unsigned int cluster_layer = TrkrDefs::getLayer(cluster_key);
109 
110  std::set<PHG4Particle*> particles = all_truth_particles(cluster_key);
111  for (std::set<PHG4Particle*>::iterator iter = particles.begin();
112  iter != particles.end();
113  ++iter)
114  {
115  PHG4Particle* particle = *iter;
116 
117  std::map<unsigned int, std::shared_ptr<TrkrCluster> > gclusters = get_truth_eval()->all_truth_clusters(particle);
118  for (std::map<unsigned int, std::shared_ptr<TrkrCluster> >::iterator citer = gclusters.begin();
119  citer != gclusters.end();
120  ++citer)
121  {
122  if(citer->first == cluster_layer)
123  {
124  truth_clusters.insert(citer->second);
125  }
126  }
127  }
128 
129  return truth_clusters;
130 }
131 
132 std::shared_ptr<TrkrCluster> SvtxClusterEval::max_truth_cluster_by_energy(TrkrDefs::cluskey cluster_key)
133 {
134  if (_do_cache)
135  {
136  std::map<TrkrDefs::cluskey, std::shared_ptr<TrkrCluster> >::iterator iter =
137  _cache_max_truth_cluster_by_energy.find(cluster_key);
138  if (iter != _cache_max_truth_cluster_by_energy.end())
139  {
140  return iter->second;
141  }
142  }
143 
144  std::shared_ptr<TrkrCluster> truth_cluster = 0;
145 
146  unsigned int cluster_layer = TrkrDefs::getLayer(cluster_key);
147 
148  PHG4Particle* max_particle = max_truth_particle_by_cluster_energy(cluster_key);
149  if(!max_particle)
150  return truth_cluster;
151 
152  if(_verbosity > 0)
153  cout << " max truth particle by cluster energy has trackID " << max_particle->get_track_id() << endl;
154 
155  TrkrCluster* reco_cluster = _clustermap->findCluster(cluster_key);
156  ActsTransformations transformer;
157  auto global = transformer.getGlobalPosition(reco_cluster,_surfmaps,_tgeometry);
158  double reco_x = global(0);
159  double reco_y = global(1);
160  double reco_z = global(2);
161  double r = sqrt(reco_x*reco_x + reco_y*reco_y);
162  //double reco_rphi = r*fast_approx_atan2(reco_y, reco_x);
163  double reco_rphi = r*atan2(reco_y, reco_x);
164 
165  std::map<unsigned int, std::shared_ptr<TrkrCluster> > gclusters = get_truth_eval()->all_truth_clusters(max_particle);
166  for (std::map<unsigned int, std::shared_ptr<TrkrCluster> >::iterator citer = gclusters.begin();
167  citer != gclusters.end();
168  ++citer)
169  {
170  if(citer->first == cluster_layer)
171  {
172  std::shared_ptr<TrkrCluster> candidate_truth_cluster = citer->second;
173 
174  double gx = candidate_truth_cluster->getX();
175  double gy = candidate_truth_cluster->getY();
176  double gz = candidate_truth_cluster->getZ();
177  double gr = sqrt(gx*gx+gy*gy);
178  double grphi = gr*atan2(gy, gx);
179  //double grphi = gr*fast_approx_atan2(gy, gx);
180 
181  // Find the difference in position from the reco cluster
182  double dz = reco_z - gz;
183  double drphi = reco_rphi - grphi;
184 
185  // approximate 4 sigmas cut
186  if(cluster_layer > 6 && cluster_layer < 23)
187  {
188  if(fabs(drphi) < 4.0 * sig_tpc_rphi_inner &&
189  fabs(dz) < 4.0 * sig_tpc_z)
190  {
191  return candidate_truth_cluster;;
192  }
193  }
194  if(cluster_layer > 22 && cluster_layer < 39)
195  {
196  if(fabs(drphi) < 4.0 * sig_tpc_rphi_mid &&
197  fabs(dz) < 4.0 * sig_tpc_z)
198  {
199  return candidate_truth_cluster;;
200  }
201  }
202  if(cluster_layer > 38 && cluster_layer < 55)
203  {
204  if(fabs(drphi) < 4.0 * sig_tpc_rphi_outer &&
205  fabs(dz) < 4.0 * sig_tpc_z)
206  {
207  return candidate_truth_cluster;;
208  }
209  }
210  else if(cluster_layer < 3)
211  {
212  if(fabs(drphi) < 4.0 * sig_mvtx_rphi &&
213  fabs(dz) < 4.0 * sig_mvtx_z)
214  {
215  return candidate_truth_cluster;;
216  }
217  }
218  else if(cluster_layer == 55)
219  {
220  if(fabs(drphi) < 4.0 * sig_mms_rphi_55)
221  {
222  return candidate_truth_cluster;;
223  }
224  }
225  else if(cluster_layer == 56)
226  {
227  if(fabs(dz) < 4.0 * sig_mms_z_56)
228  {
229  return candidate_truth_cluster;;
230  }
231  }
232  else
233  {
234  if(fabs(drphi) < 4.0 * sig_intt_rphi &&
235  fabs(dz) < range_intt_z)
236  {
237  return candidate_truth_cluster;;
238  }
239  }
240  }
241  }
242 
243  return truth_cluster;
244 }
245 
247 {
248  if (_do_cache)
249  {
250  std::map<std::shared_ptr<TrkrCluster>, TrkrCluster* >::iterator iter =
252  if (iter != _cache_reco_cluster_from_truth_cluster.end())
253  {
254  return iter->second;
255  }
256  }
257 
258  TrkrDefs::cluskey ckey = gclus->getClusKey();
259  double gx = gclus->getX();
260  double gy = gclus->getY();
261  double gz = gclus->getZ();
262  double gr = sqrt(gx*gx+gy*gy);
263  double grphi = gr*atan2(gy, gx);
264  //double grphi = gr*fast_approx_atan2(gy, gx);
265 
266  unsigned int truth_layer = TrkrDefs::getLayer(ckey);
267  TrkrCluster *reco_cluster = 0;
268 
269  std::set<TrkrDefs::cluskey> reco_cluskeys;
270  std::set<PHG4Hit*> contributing_hits = get_truth_eval()->get_truth_hits_from_truth_cluster(ckey);
271  for(auto it = contributing_hits.begin(); it != contributing_hits.end(); ++it)
272  {
273  PHG4Hit* cont_g4hit = *it;
274  std::set<TrkrDefs::cluskey> cluskeys = all_clusters_from(cont_g4hit); // this returns clusters from this hit in any layer using TrkrAssoc maps
275 
276  if(_verbosity > 0)
277  cout << " contributing g4hitID " << cont_g4hit->get_hit_id() << " g4trackID " << cont_g4hit->get_trkid() << endl;
278 
279  for (std::set<TrkrDefs::cluskey>::iterator iter = cluskeys.begin();
280  iter != cluskeys.end();
281  ++iter)
282  {
283  unsigned int clus_layer = TrkrDefs::getLayer(*iter);
284  // discard if reco cluster is in the wrong layer
285  if(clus_layer != truth_layer) continue;
286 
287  reco_cluskeys.insert(*iter);
288  }
289  }
290 
291  unsigned int nreco = reco_cluskeys.size();
292  if(nreco > 0)
293  {
294  // Find a matching reco cluster with position inside 4 sigmas, and replace reco_cluskey
296  for(std::set<TrkrDefs::cluskey>::iterator it = reco_cluskeys.begin(); it != reco_cluskeys.end(); ++it)
297  {
298  // get the cluster
299  TrkrCluster* this_cluster = _clustermap->findCluster(*it);
300  auto global = transform.getGlobalPosition(this_cluster,_surfmaps,_tgeometry);
301  double this_x = global(0);
302  double this_y = global(1);
303  double this_z = global(2);
304  double this_rphi = gr*atan2(this_y, this_x);
305  //double this_rphi = gr*fast_approx_atan2(this_y, this_x);
306 
307  // Find the difference in position from the g4cluster
308  double dz = this_z - gz;
309  double drphi = this_rphi - grphi;
310 
311  // approximate 4 sigmas cut
312  if(truth_layer > 6 && truth_layer < 23)
313  {
314  if(fabs(drphi) < 4.0 * sig_tpc_rphi_inner &&
315  fabs(dz) < 4.0 * sig_tpc_z)
316  {
317  return this_cluster;;
318  }
319  }
320  if(truth_layer > 22 && truth_layer < 39)
321  {
322  if(fabs(drphi) < 4.0 * sig_tpc_rphi_mid &&
323  fabs(dz) < 4.0 * sig_tpc_z)
324  {
325  return this_cluster;;
326  }
327  }
328  if(truth_layer > 38 && truth_layer < 55)
329  {
330  if(fabs(drphi) < 4.0 * sig_tpc_rphi_outer &&
331  fabs(dz) < 4.0 * sig_tpc_z)
332  {
333  return this_cluster;;
334  }
335  }
336  else if(truth_layer < 3)
337  {
338  if(fabs(drphi) < 4.0 * sig_mvtx_rphi &&
339  fabs(dz) < 4.0 * sig_mvtx_z)
340  {
341  return this_cluster;;
342  }
343  }
344  else if(truth_layer == 55)
345  {
346  if(fabs(drphi) < 4.0 * sig_mms_rphi_55)
347  {
348  return this_cluster;;
349  }
350  }
351  else if(truth_layer == 56)
352  {
353  if(fabs(dz) < 4.0 * sig_mms_z_56)
354  {
355  return this_cluster;;
356  }
357  }
358  else
359  {
360  if(fabs(drphi) < 4.0 * sig_intt_rphi &&
361  fabs(dz) < range_intt_z)
362  {
363  return this_cluster;;
364  }
365  }
366  }
367  }
368 
369  return reco_cluster;
370 }
371 
372 std::set<PHG4Hit*> SvtxClusterEval::all_truth_hits(TrkrDefs::cluskey cluster_key)
373 {
374  if (!has_node_pointers())
375  {
376  ++_errors;
377  return std::set<PHG4Hit*>();
378  }
379 
380  if (_do_cache)
381  {
382  std::map<TrkrDefs::cluskey, std::set<PHG4Hit*> >::iterator iter =
383  _cache_all_truth_hits.find(cluster_key);
384  if (iter != _cache_all_truth_hits.end())
385  {
386  return iter->second;
387  }
388  }
389 
390  std::set<PHG4Hit*> truth_hits;
391 
392  // get all truth hits for this cluster
393  //_cluster_hit_map->identify();
394  std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
395  hitrange = _cluster_hit_map->getHits(cluster_key); // returns range of pairs {cluster key, hit key} for this cluskey
396 
397  for(std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
398  clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
399  {
400  TrkrDefs::hitkey hitkey = clushititer->second;
401  // TrkrHitTruthAssoc uses a map with (hitsetkey, std::pair(hitkey, g4hitkey)) - get the hitsetkey from the cluskey
403 
404  // get all of the g4hits for this hitkey
405  std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> > temp_map;
406  _hit_truth_map->getG4Hits(hitsetkey, hitkey, temp_map);
407  // returns pairs (hitsetkey, std::pair(hitkey, g4hitkey)) for this hitkey only
408 
409  for(std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> >::iterator htiter = temp_map.begin(); htiter != temp_map.end(); ++htiter)
410  {
411 
412  // extract the g4 hit key here and add the hits to the set
413  PHG4HitDefs::keytype g4hitkey = htiter->second.second;
414  PHG4Hit * g4hit = nullptr;
415  unsigned int trkrid = TrkrDefs::getTrkrId(hitsetkey);
416  switch( trkrid )
417  {
418  case TrkrDefs::tpcId: g4hit = _g4hits_tpc->findHit(g4hitkey); break;
419  case TrkrDefs::inttId: g4hit = _g4hits_intt->findHit(g4hitkey); break;
420  case TrkrDefs::mvtxId: g4hit = _g4hits_mvtx->findHit(g4hitkey); break;
421  case TrkrDefs::micromegasId: g4hit = _g4hits_mms->findHit(g4hitkey); break;
422  default: break;
423  }
424  if( g4hit ) truth_hits.insert(g4hit);
425  } // end loop over g4hits associated with hitsetkey and hitkey
426  } // end loop over hits associated with cluskey
427 
428  if (_do_cache) _cache_all_truth_hits.insert(make_pair(cluster_key, truth_hits));
429 
430  return truth_hits;
431 }
432 
434 {
435  if (!has_node_pointers())
436  {
437  ++_errors;
438  return nullptr;
439  }
440 
441 // if (_strict)
442 // {
443 // assert(cluster_key);
444 // }
445 // else if (!cluster_key)
446 // {
447 // ++_errors;
448 // return std::set<PHG4Hit*>();
449 // }
450  /*
451  if (_do_cache)
452  {
453  std::map<TrkrDefs::cluskey, std::set<PHG4Hit*> >::iterator iter =
454  _cache_all_truth_hits.find(cluster_key);
455  if (iter != _cache_all_truth_hits.end())
456  {
457  return iter->second;
458  }
459  }
460  */
461  ActsTransformations transformer;
462  TrkrCluster* cluster = _clustermap->findCluster(cluster_key);
463  auto glob = transformer.getGlobalPosition(cluster,_surfmaps,_tgeometry);
464  TVector3 cvec(glob(0), glob(1), glob(2));
465  unsigned int layer = TrkrDefs::getLayer(cluster_key);
466  std::set<PHG4Hit*> truth_hits;
467 
468  std::multimap<PHG4HitDefs::keytype,TrkrDefs::hitkey> g4keyperhit;
469  std::vector<PHG4HitDefs::keytype> g4hitkeys;
470  // get all truth hits for this cluster
471  //_cluster_hit_map->identify();
473 
474  std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
475  hitrange = _cluster_hit_map->getHits(cluster_key); // returns range of pairs {cluster key, hit key} for this cluskey
476  for(std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
477  clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
478  {
479  TrkrDefs::hitkey hitkey = clushititer->second;
480  // TrkrHitTruthAssoc uses a map with (hitsetkey, std::pair(hitkey, g4hitkey)) - get the hitsetkey from the cluskey
481 
482  // get all of the g4hits for this hitkey
483  std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> > temp_map;
484  _hit_truth_map->getG4Hits(hitsetkey, hitkey, temp_map); // returns pairs (hitsetkey, std::pair(hitkey, g4hitkey)) for this hitkey only
485  for(std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> >::iterator htiter = temp_map.begin(); htiter != temp_map.end(); ++htiter)
486  {
487  // extract the g4 hit key here and add the hits to the set
488  PHG4HitDefs::keytype g4hitkey = htiter->second.second;
489  if(_verbosity > 2)
490  cout << " g4key: " << g4hitkey << " layer: " << layer << endl;
491  TrkrDefs::hitkey hitkey = htiter->second.first;
492  /* if(layer>=7){
493  PHG4Hit *match_g4hit = _g4hits_tpc->findHit(g4hitkey);
494  if(layer != match_g4hit->get_layer() ) continue;
495  }
496  */
497  g4keyperhit.insert(std::pair<PHG4HitDefs::keytype,TrkrDefs::hitkey>(g4hitkey,hitkey));
498  std::vector<PHG4HitDefs::keytype>::iterator itg4keys = find(g4hitkeys.begin(),g4hitkeys.end(),g4hitkey);
499  if(itg4keys==g4hitkeys.end()) g4hitkeys.push_back(g4hitkey);
500  } // end loop over g4hits associated with hitsetkey and hitkey
501  } // end loop over hits associated with cluskey
502 
503  // if (_do_cache) _cache_all_truth_hits.insert(make_pair(cluster_key, truth_hits));
504  PHG4HitDefs::keytype max_key = 0;
505  unsigned int n_max = 0;
506 
507  if(g4hitkeys.size()==1 ){
508  std::vector<PHG4HitDefs::keytype>::iterator it = g4hitkeys.begin();
509  max_key = *it;
510  }else{
511  for(std::vector<PHG4HitDefs::keytype>::iterator it = g4hitkeys.begin(); it != g4hitkeys.end(); ++it){
512  unsigned int ng4hit = g4keyperhit.count(*it);
513  PHG4Hit * this_g4hit = _g4hits_tpc->findHit(*it);
514 
515  if(layer >= 7 ){ //in tpc
516  if(this_g4hit!=NULL){
517  unsigned int glayer = this_g4hit->get_layer();
518  if(layer != glayer) continue;
519 
520  TVector3 vec(this_g4hit->get_avg_x(), this_g4hit->get_avg_y(), this_g4hit->get_avg_z());
521  //cout << "layer: " << layer << " (" << glayer << ") " << " gtrackID: " << this_g4hit->get_trkid() << " novlp: " << ng4hit << " phi: " << vec.Phi() << " z: " << this_g4hit->get_avg_z() << " r: " << vec.Perp() << " keyg4: " << *it << endl; //<< " keyrec: "<< *it.second << endl;
522  }
523  /*else{
524  cout << "g4hit == NULL " << endl;
525  }
526  */
527  }
528  if(ng4hit>n_max){
529  max_key = *it;
530  n_max = ng4hit;
531  }
532 
533  }
534  }
535 
536  PHG4Hit * g4hit = nullptr;
537  unsigned int trkrid = TrkrDefs::getTrkrId(hitsetkey);
538  switch( trkrid )
539  {
540  case TrkrDefs::tpcId: g4hit = _g4hits_tpc->findHit(max_key); break;
541  case TrkrDefs::inttId: g4hit = _g4hits_intt->findHit(max_key); break;
542  case TrkrDefs::mvtxId: g4hit = _g4hits_mvtx->findHit(max_key); break;
543  case TrkrDefs::micromegasId: g4hit = _g4hits_mms->findHit(max_key); break;
544  default: break;
545  }
546  if( g4hit ) truth_hits.insert(g4hit);
547 
548  return g4hit;
549 }
550 
552 {
553  if (!has_node_pointers())
554  {
555  ++_errors;
556  return make_pair(0,0);
557  }
558 
559 // if (_strict)
560 // {
561 // assert(cluster_key);
562 // }
563 // else if (!cluster_key)
564 // {
565 // ++_errors;
566 // return std::set<PHG4Hit*>();
567 // }
568  /*
569  if (_do_cache)
570  {
571  std::map<TrkrDefs::cluskey, std::set<PHG4Hit*> >::iterator iter =
572  _cache_all_truth_hits.find(cluster_key);
573  if (iter != _cache_all_truth_hits.end())
574  {
575  return iter->second;
576  }
577  }
578  */
579 
580  std::pair<int, int> out_pair;
581  out_pair.first = 0;
582  out_pair.second = -1;
584 
585  TrkrCluster* cluster = _clustermap->findCluster(cluster_key);
586  auto global = transform.getGlobalPosition(cluster,_surfmaps,_tgeometry);
587  TVector3 cvec(global(0), global(1), global(2));
588  unsigned int layer = TrkrDefs::getLayer(cluster_key);
589 
590  std::multimap<PHG4HitDefs::keytype,TrkrDefs::hitkey> g4keyperhit;
591  std::vector<PHG4HitDefs::keytype> g4hitkeys;
592  // get all truth hits for this cluster
593  //_cluster_hit_map->identify();
595 
596  std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
597  hitrange = _cluster_hit_map->getHits(cluster_key); // returns range of pairs {cluster key, hit key} for this cluskey
598  for(std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
599  clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
600  {
601  TrkrDefs::hitkey hitkey = clushititer->second;
602  // TrkrHitTruthAssoc uses a map with (hitsetkey, std::pair(hitkey, g4hitkey)) - get the hitsetkey from the cluskey
603 
604  // get all of the g4hits for this hitkey
605  std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> > temp_map;
606  _hit_truth_map->getG4Hits(hitsetkey, hitkey, temp_map); // returns pairs (hitsetkey, std::pair(hitkey, g4hitkey)) for this hitkey only
607 
608  for(std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> >::iterator htiter = temp_map.begin(); htiter != temp_map.end(); ++htiter)
609  {
610  // extract the g4 hit key here and add the hits to the set
611  PHG4HitDefs::keytype g4hitkey = htiter->second.second;
612  if(_verbosity > 2)
613  cout << " g4key: " << g4hitkey << " layer: " << layer << endl;
614  TrkrDefs::hitkey hitkey = htiter->second.first;
615  /* if(layer>=7){
616  PHG4Hit *match_g4hit = _g4hits_tpc->findHit(g4hitkey);
617  if(layer != match_g4hit->get_layer() ) continue;
618  }
619  */
620  g4keyperhit.insert(std::pair<PHG4HitDefs::keytype,TrkrDefs::hitkey>(g4hitkey,hitkey));
621  std::vector<PHG4HitDefs::keytype>::iterator itg4keys = find(g4hitkeys.begin(),g4hitkeys.end(),g4hitkey);
622  if(itg4keys==g4hitkeys.end()) g4hitkeys.push_back(g4hitkey);
623  } // end loop over g4hits associated with hitsetkey and hitkey
624  } // end loop over hits associated with cluskey
625 
626  PHG4HitDefs::keytype max_key = 0;
627  unsigned int n_max = 0;
628  if(_verbosity > 2)
629  cout << " n matches found: " << g4hitkeys.size() << " phi: " << cvec.Phi() << " z: " << cvec.Z() << " ckey: " << cluster_key << endl;
630 
631  if(g4hitkeys.size()==1 ){
632  std::vector<PHG4HitDefs::keytype>::iterator it = g4hitkeys.begin();
633  max_key = *it;
634  }else{
635  for(std::vector<PHG4HitDefs::keytype>::iterator it = g4hitkeys.begin(); it != g4hitkeys.end(); ++it){
636  unsigned int ng4hit = g4keyperhit.count(*it);
637  PHG4Hit * this_g4hit = _g4hits_tpc->findHit(*it);
638 
639 
640  if(layer >= 7 ){ //in tpc
641  if(this_g4hit!=NULL){
642  unsigned int glayer = this_g4hit->get_layer();
643  // if(layer != glayer) continue;
644 
645  TVector3 vec(this_g4hit->get_avg_x(), this_g4hit->get_avg_y(), this_g4hit->get_avg_z());
646  if(_verbosity > 2)
647  cout << "layer: " << layer << " (" << glayer << ") " << " gtrackID: " << this_g4hit->get_trkid() << " novlp: " << ng4hit << " phi: " << vec.Phi() << " z: " << this_g4hit->get_avg_z() << " r: " << vec.Perp() << " keyg4: " << *it << " cz: " << cluster->getZ() << endl; //<< " keyrec: "<< *it.second << endl;
648  }
649  }
650  if(ng4hit>n_max){
651  max_key = *it;
652  n_max = ng4hit;
653  }
654 
655  }
656  }
657  if(_verbosity > 2)
658  cout << "found in layer: " << layer << " n_max: " << n_max << " max_key: " << max_key << " ckey: " << cluster_key << endl;
659  if(max_key !=0){
660  PHG4Hit * g4hit = nullptr;
661  unsigned int trkrid = TrkrDefs::getTrkrId(hitsetkey);
662  switch( trkrid )
663  {
664  case TrkrDefs::tpcId: g4hit = _g4hits_tpc->findHit(max_key); break;
665  case TrkrDefs::inttId: g4hit = _g4hits_intt->findHit(max_key); break;
666  case TrkrDefs::mvtxId: g4hit = _g4hits_mvtx->findHit(max_key); break;
667  case TrkrDefs::micromegasId: g4hit = _g4hits_mms->findHit(max_key); break;
668  default: break;
669  }
670 
671  //check if we on a looper
672  PHG4Particle* g4particle = _truthinfo->GetParticle(g4hit->get_trkid());
673 
674  PHG4VtxPoint* vtx = _truthinfo->GetVtx(g4particle->get_vtx_id());
675  float vtx_z = vtx->get_z();
676  float gpx = g4particle->get_px();
677  float gpy = g4particle->get_py();
678  float gpz = g4particle->get_pz();
679  float gpeta = NAN;
680 
681  TVector3 gv(gpx, gpy, gpz);
682  gpeta = gv.Eta();
683  TVector3 this_vec( g4hit->get_avg_x() ,
684  g4hit->get_avg_y() ,
685  g4hit->get_avg_z() - vtx_z);
686  double deta = TMath::Abs(gpeta - this_vec.Eta());
687 
688  int is_loop = 0;
689 
690  if(layer >= 7){
691  // cout << " in tpc " << endl;
692  if(deta>0.1) is_loop = 1;
693  }
694 
695  out_pair.first = g4hit->get_trkid();
696  if(!is_loop)
697  out_pair.second = layer;
698  }
699  return out_pair;
700 }
701 
703 {
704  if (!has_node_pointers())
705  {
706  ++_errors;
707  return nullptr;
708  }
709 
710  if (_do_cache)
711  {
712  std::map<TrkrDefs::cluskey, PHG4Hit*>::iterator iter =
713  _cache_max_truth_hit_by_energy.find(cluster_key);
714  if (iter != _cache_max_truth_hit_by_energy.end())
715  {
716  return iter->second;
717  }
718  }
719 
720  std::set<PHG4Hit*> hits = all_truth_hits(cluster_key);
721  PHG4Hit* max_hit = nullptr;
722  float max_e = FLT_MAX * -1.0;
723  for (std::set<PHG4Hit*>::iterator iter = hits.begin();
724  iter != hits.end();
725  ++iter)
726  {
727  PHG4Hit* hit = *iter;
728  if (hit->get_edep() > max_e)
729  {
730  max_e = hit->get_edep();
731  max_hit = hit;
732  }
733  }
734 
735  if (_do_cache) _cache_max_truth_hit_by_energy.insert(make_pair(cluster_key, max_hit));
736 
737  return max_hit;
738 }
739 
740 std::set<PHG4Particle*> SvtxClusterEval::all_truth_particles(TrkrDefs::cluskey cluster_key)
741 {
742  if (!has_node_pointers())
743  {
744  ++_errors;
745  return std::set<PHG4Particle*>();
746  }
747 
748  if (_do_cache)
749  {
750  std::map<TrkrDefs::cluskey, std::set<PHG4Particle*> >::iterator iter =
751  _cache_all_truth_particles.find(cluster_key);
752  if (iter != _cache_all_truth_particles.end())
753  {
754  return iter->second;
755  }
756  }
757 
758  std::set<PHG4Particle*> truth_particles;
759 
760  std::set<PHG4Hit*> g4hits = all_truth_hits(cluster_key);
761 
762  for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
763  iter != g4hits.end();
764  ++iter)
765  {
766  PHG4Hit* hit = *iter;
768  //cout << "cluster key " << cluster_key << " has hit " << hit->get_hit_id() << " and has particle " << particle->get_track_id() << endl;
769 
770  if (_strict)
771  {
772  assert(particle);
773  }
774  else if (!particle)
775  {
776  ++_errors;
777  continue;
778  }
779 
780  truth_particles.insert(particle);
781  }
782 
783  if (_do_cache) _cache_all_truth_particles.insert(make_pair(cluster_key, truth_particles));
784 
785  return truth_particles;
786 }
787 
789 {
790  if (!has_node_pointers())
791  {
792  ++_errors;
793  return nullptr;
794  }
795 
796  if (_do_cache)
797  {
798  std::map<TrkrDefs::cluskey, PHG4Particle*>::iterator iter =
801  {
802  return iter->second;
803  }
804  }
805 
806  unsigned int layer = TrkrDefs::getLayer(cluster_key);
807 
808  // loop over all particles associated with this cluster and
809  // get the energy contribution for each one, record the max
810  PHG4Particle* max_particle = nullptr;
811  float max_e = FLT_MAX * -1.0;
812  std::set<PHG4Particle*> particles = all_truth_particles(cluster_key);
813  for (std::set<PHG4Particle*>::iterator iter = particles.begin();
814  iter != particles.end();
815  ++iter)
816  {
817  PHG4Particle* particle = *iter;
818  std::map<unsigned int, std::shared_ptr<TrkrCluster> > truth_clus = get_truth_eval()->all_truth_clusters(particle);
819  for(auto it = truth_clus.begin(); it != truth_clus.end(); ++it)
820  {
821  if(it->first == layer)
822  {
823  float e = it->second->getError(0,0);
824  if (e > max_e)
825  {
826  max_e = e;
827  max_particle = particle;
828  }
829  }
830  }
831  }
832 
833  if (_do_cache) _cache_max_truth_particle_by_cluster_energy.insert(make_pair(cluster_key, max_particle));
834 
835  return max_particle;
836 }
837 
839 {
840  // Note: this does not quite work correctly for the TPC - it assumes one g4hit per layer
841  // use max_truth_particle_by_cluster_energy instead
842 
843  if (!has_node_pointers())
844  {
845  ++_errors;
846  return nullptr;
847  }
848 
849  if (_do_cache)
850  {
851  std::map<TrkrDefs::cluskey, PHG4Particle*>::iterator iter =
852  _cache_max_truth_particle_by_energy.find(cluster_key);
853  if (iter != _cache_max_truth_particle_by_energy.end())
854  {
855  return iter->second;
856  }
857  }
858 
859  // loop over all particles associated with this cluster and
860  // get the energy contribution for each one, record the max
861  PHG4Particle* max_particle = nullptr;
862  float max_e = FLT_MAX * -1.0;
863  std::set<PHG4Particle*> particles = all_truth_particles(cluster_key);
864  for (std::set<PHG4Particle*>::iterator iter = particles.begin();
865  iter != particles.end();
866  ++iter)
867  {
868  PHG4Particle* particle = *iter;
869  float e = get_energy_contribution(cluster_key, particle);
870  if (e > max_e)
871  {
872  max_e = e;
873  max_particle = particle;
874  }
875  }
876 
877  if (_do_cache) _cache_max_truth_particle_by_energy.insert(make_pair(cluster_key, max_particle));
878 
879  return max_particle;
880 }
881 
882 std::set<TrkrDefs::cluskey> SvtxClusterEval::all_clusters_from(PHG4Particle* truthparticle)
883 {
884  if (!has_node_pointers())
885  {
886  ++_errors;
887  return std::set<TrkrDefs::cluskey>();
888  }
889 
890  if (_strict)
891  {
892  assert(truthparticle);
893  }
894  else if (!truthparticle)
895  {
896  ++_errors;
897  return std::set<TrkrDefs::cluskey>();
898  }
899  //check if cache is filled, if not fill it.
900  // if(_cache_all_clusters_from_particle.count(truthparticle)==0){
903  }
904 
905  if (_do_cache)
906  {
907  std::map<PHG4Particle*, std::set<TrkrDefs::cluskey> >::iterator iter =
908  _cache_all_clusters_from_particle.find(truthparticle);
909  if (iter != _cache_all_clusters_from_particle.end())
910  {
911  return iter->second;
912  }
913  }
914  std::set<TrkrDefs::cluskey> clusters;
915  return clusters;
916 }
917 
919  auto Mytimer = std::make_unique<PHTimer>("ReCl_timer");
920  Mytimer->stop();
921  Mytimer->restart();
922 
923  std::multimap<PHG4Particle*, TrkrDefs::cluskey> temp_clusters_from_particles;
924  // loop over all the clusters
925  auto hitsetrange = _hitsets->getHitSets();
926  for (auto hitsetitr = hitsetrange.first;
927  hitsetitr != hitsetrange.second;
928  ++hitsetitr){
929  auto range = _clustermap->getClusters(hitsetitr->first);
930  for( auto iter = range.first; iter != range.second; ++iter ){
931  TrkrDefs::cluskey cluster_key = iter->first;
932 
933  // loop over all truth particles connected to this cluster
934  std::set<PHG4Particle*> particles = all_truth_particles(cluster_key);
935  for (std::set<PHG4Particle*>::iterator jter = particles.begin();
936  jter != particles.end();
937  ++jter){
938  PHG4Particle* candidate = *jter;
939  temp_clusters_from_particles.insert(make_pair(candidate, cluster_key));
940  }
941  }
942  }
943  //Loop over particles and fill cache
945  for(PHG4TruthInfoContainer::ConstIterator iter = range.first;
946  iter != range.second; ++iter){
947  PHG4Particle* g4particle = iter->second;
948  std::set<TrkrDefs::cluskey> clusters;
949  std::multimap<PHG4Particle*, TrkrDefs::cluskey>::const_iterator lower_bound = temp_clusters_from_particles.lower_bound(g4particle);
950  std::multimap<PHG4Particle*, TrkrDefs::cluskey>::const_iterator upper_bound = temp_clusters_from_particles.upper_bound(g4particle);
951  std::multimap<PHG4Particle*, TrkrDefs::cluskey>::const_iterator cfp_iter;
952  for(cfp_iter = lower_bound;cfp_iter != upper_bound;++cfp_iter){
953  TrkrDefs::cluskey cluster_key = cfp_iter->second;
954  clusters.insert(cluster_key);
955  }
956  _cache_all_clusters_from_particle.insert(make_pair(g4particle, clusters));
957  }
958 
959  Mytimer->stop();
960 
961 }
962 
963 std::set<TrkrDefs::cluskey> SvtxClusterEval::all_clusters_from(PHG4Hit* truthhit)
964 {
965  if (!has_node_pointers())
966  {
967  ++_errors;
968  return std::set<TrkrDefs::cluskey>();
969  }
970 
971  if (_strict)
972  {
973  assert(truthhit);
974  }
975  else if (!truthhit)
976  {
977  ++_errors;
978  return std::set<TrkrDefs::cluskey>();
979  }
980 
981  // one time, fill cache of g4hit/cluster pairs
982  if(_cache_all_clusters_from_g4hit.size() == 0)
983  {
984  // make a map of truthhit, cluster_key inside this loop
985  std::multimap<PHG4HitDefs::keytype, TrkrDefs::cluskey> truth_cluster_map;
986  std::set<PHG4HitDefs::keytype> all_g4hits_set;
987  std::map<PHG4HitDefs::keytype, PHG4Hit*> all_g4hits_map;
988 
989  // get all reco clusters
990  if(_verbosity > 1) cout << "all_clusters_from_g4hit: list all reco clusters " << endl;
991 
992  auto hitsetrange = _hitsets->getHitSets();
993  for (auto hitsetitr = hitsetrange.first;
994  hitsetitr != hitsetrange.second;
995  ++hitsetitr){
996  auto range = _clustermap->getClusters(hitsetitr->first);
997  for( auto iter = range.first; iter != range.second; ++iter ){
998  TrkrDefs::cluskey cluster_key = iter->first;
999  int layer = TrkrDefs::getLayer(cluster_key);
1000  TrkrCluster *clus = iter->second;
1001  if(_verbosity > 1)
1002  {
1003  cout << " layer " << layer << " cluster_key " << cluster_key << " adc " << clus->getAdc()
1004  << " localx " << clus->getLocalX()
1005  << " localy " << clus->getLocalY()
1006  << endl;
1007  cout << " associated hits:";
1008  std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
1009  hitrange = _cluster_hit_map->getHits(cluster_key); // returns range of pairs {cluster key, hit key} for this cluskey
1010  for(std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
1011  clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
1012  {
1013  TrkrDefs::hitkey hitkey = clushititer->second;
1014  cout << " " << hitkey;
1015  }
1016  cout << endl;
1017  }
1018 
1019  // the returned truth hits were obtained from TrkrAssoc maps
1020  std::set<PHG4Hit*> hits = all_truth_hits(cluster_key);
1021  for (std::set<PHG4Hit*>::iterator jter = hits.begin();
1022  jter != hits.end();
1023  ++jter)
1024  {
1025  PHG4Hit* candidate = *jter;
1026  PHG4HitDefs::keytype g4hit_key = candidate->get_hit_id();
1027 
1028  if(_verbosity > 5)
1029  {
1030  int gtrackID = candidate->get_trkid();
1031  cout << " adding cluster with cluster_key " << cluster_key << " g4hit with g4hit_key " << g4hit_key
1032  << " gtrackID " << gtrackID
1033  << endl;
1034  }
1035 
1036  all_g4hits_set.insert(g4hit_key);
1037  all_g4hits_map.insert(std::make_pair(g4hit_key, candidate));
1038 
1039  truth_cluster_map.insert(std::make_pair(g4hit_key, cluster_key));
1040  }
1041  }
1042  }
1043 
1044  // now fill the cache
1045  // loop over all entries in all_g4hits
1046  for(std::set<PHG4HitDefs::keytype>::iterator iter = all_g4hits_set.begin(); iter != all_g4hits_set.end(); ++iter)
1047  {
1048  PHG4HitDefs::keytype g4hit_key = *iter;
1049  if(_verbosity > 5) cout << " associations for g4hit_key " << g4hit_key << endl;
1050 
1051  std::map<PHG4HitDefs::keytype, PHG4Hit*>::iterator it = all_g4hits_map.find(g4hit_key);
1052  PHG4Hit *g4hit = it->second;
1053 
1054  std::set<TrkrDefs::cluskey> assoc_clusters;
1055 
1056  std::pair<std::multimap<PHG4HitDefs::keytype, TrkrDefs::cluskey>::iterator,
1057  std::multimap<PHG4HitDefs::keytype, TrkrDefs::cluskey>::iterator> ret = truth_cluster_map.equal_range(g4hit_key);
1058  for(std::multimap<PHG4HitDefs::keytype, TrkrDefs::cluskey>::iterator jter = ret.first; jter != ret.second; ++jter)
1059  {
1060  assoc_clusters.insert(jter->second);
1061 
1062  if(_verbosity > 5) cout << " g4hit_key " << g4hit_key << " associated with cluster_key " << jter->second << endl;
1063  }
1064  // done with this g4hit
1065  _cache_all_clusters_from_g4hit.insert(make_pair(g4hit, assoc_clusters));
1066  }
1067 
1068  }
1069 
1070  // get the clusters
1071  std::set<TrkrDefs::cluskey> clusters;
1072  std::map<PHG4Hit*, std::set<TrkrDefs::cluskey> >::iterator iter =
1073  _cache_all_clusters_from_g4hit.find(truthhit);
1074  if (iter != _cache_all_clusters_from_g4hit.end())
1075  {
1076  return iter->second;
1077  }
1078 
1079 
1080  if (_clusters_per_layer.size() == 0)
1081  {
1083  }
1084 
1085  return clusters;
1086 }
1087 
1089 {
1090  TrkrDefs::cluskey val =0;
1091  if (!has_node_pointers())
1092  {
1093  ++_errors;
1094  return val;
1095  }
1096 
1097  /* if (_strict)
1098  {
1099  assert(truthhit);
1100  }
1101  else if (!truthhit)
1102  {
1103  ++_errors;
1104  return val;
1105  }
1106  */
1107  // one time, fill cache of g4hit/cluster pairs
1109  // get all reco clusters
1110  // cout << "cache size ==0" << endl;
1111  if(_verbosity > 1) cout << "all_clusters: found # " << _clustermap->size() << endl;
1112 
1113  // loop over clusters and get all contributing hits
1114  auto hitsetrange = _hitsets->getHitSets();
1115  for (auto hitsetitr = hitsetrange.first;
1116  hitsetitr != hitsetrange.second;
1117  ++hitsetitr){
1118  auto range = _clustermap->getClusters(hitsetitr->first);
1119  for( auto iter = range.first; iter != range.second; ++iter ){
1120  TrkrDefs::cluskey cluster_key = iter->first;
1121  int layer_in = TrkrDefs::getLayer(cluster_key);
1122 
1123  if(layer_in<0) continue;
1124  // TrkrCluster *clus = iter->second;
1125 
1126  std::pair<int, int> gid_lay = gtrackid_and_layer_by_nhit(cluster_key);
1127 
1128  // std::map<std::pair<int, unsigned int>, TrkrDefs::cluskey>::iterator it_exists;
1129  // it_exists =
1130  if(_cache_best_cluster_from_gtrackid_layer.count(gid_lay)==0){
1131  if(gid_lay.second >=0)
1132  _cache_best_cluster_from_gtrackid_layer.insert(make_pair(gid_lay, cluster_key));
1133  }
1134  else
1135  if(_verbosity > 2){ cout << "found doublematch" << endl;
1136  cout << "ckey: " << cluster_key << " gtrackID: " << gid_lay.first << " layer: " << gid_lay.second << endl;
1137  }
1138  }
1139  }
1140  }
1141 
1142 
1143  // get the clusters
1144  TrkrDefs::cluskey best_cluster = 0;
1145  // PHG4Hit*, std::set<TrkrDefs::cluskey> >::iterator iter =
1146 
1147  std::map<std::pair<int, int>, TrkrDefs::cluskey>::iterator iter = _cache_best_cluster_from_gtrackid_layer.find(make_pair(gid,layer));
1148  if (iter != _cache_best_cluster_from_gtrackid_layer.end())
1149  {
1150  return iter->second;
1151  }
1152 
1153  return best_cluster;
1154 }
1155 
1157 {
1158  if (!has_node_pointers())
1159  {
1160  ++_errors;
1161  return 0;
1162  }
1163 
1164  if (_strict)
1165  {
1166  assert(truthhit);
1167  }
1168  else if (!truthhit)
1169  {
1170  ++_errors;
1171  return 0;
1172  }
1173 
1174  if (_do_cache)
1175  {
1176  std::map<PHG4Hit*, TrkrDefs::cluskey>::iterator iter =
1177  _cache_best_cluster_from_g4hit.find(truthhit);
1178  if (iter != _cache_best_cluster_from_g4hit.end())
1179  {
1180  return iter->second;
1181  }
1182  }
1183 
1184  TrkrDefs::cluskey best_cluster = 0;
1185  float best_energy = 0.0;
1186  std::set<TrkrDefs::cluskey> clusters = all_clusters_from(truthhit);
1187  for (std::set<TrkrDefs::cluskey>::iterator iter = clusters.begin();
1188  iter != clusters.end();
1189  ++iter)
1190  {
1191  TrkrDefs::cluskey cluster_key = *iter;
1192  float energy = get_energy_contribution(cluster_key, truthhit);
1193  if (energy > best_energy)
1194  {
1195  best_cluster = cluster_key;
1196  best_energy = energy;
1197  }
1198  }
1199 
1200  if (_do_cache) _cache_best_cluster_from_g4hit.insert(make_pair(truthhit, best_cluster));
1201 
1202  return best_cluster;
1203 }
1204 
1205 // overlap calculations
1207 {
1208  // Note: this does not work correctly for the TPC
1209  // It assumes one g4hit per layer. Use the truth cluster energy instead.
1210 
1211  if (!has_node_pointers())
1212  {
1213  ++_errors;
1214  return NAN;
1215  }
1216 
1217  if (_strict)
1218  {
1219 // assert(cluster_key);
1220  assert(particle);
1221  }
1222  else if ( !particle)
1223  {
1224  ++_errors;
1225  return NAN;
1226  }
1227 
1228  if (_do_cache)
1229  {
1230  std::map<std::pair<TrkrDefs::cluskey, PHG4Particle*>, float>::iterator iter =
1231  _cache_get_energy_contribution_g4particle.find(make_pair(cluster_key, particle));
1233  {
1234  return iter->second;
1235  }
1236  }
1237 
1238  float energy = 0.0;
1239  std::set<PHG4Hit*> hits = all_truth_hits(cluster_key);
1240  for (std::set<PHG4Hit*>::iterator iter = hits.begin();
1241  iter != hits.end();
1242  ++iter)
1243  {
1244  PHG4Hit* hit = *iter;
1245  if (get_truth_eval()->is_g4hit_from_particle(hit, particle))
1246  {
1247  energy += hit->get_edep();
1248  }
1249  }
1250 
1251  if (_do_cache) _cache_get_energy_contribution_g4particle.insert(make_pair(make_pair(cluster_key, particle), energy));
1252 
1253  return energy;
1254 }
1255 
1257 {
1258  if (!has_node_pointers())
1259  {
1260  ++_errors;
1261  return NAN;
1262  }
1263 
1264  if (_strict)
1265  {
1266 // assert(cluster_key);
1267  assert(g4hit);
1268  }
1269  else if ( !g4hit)
1270  {
1271  ++_errors;
1272  return NAN;
1273  }
1274 
1275  if ((_do_cache) &&
1276  (_cache_get_energy_contribution_g4hit.find(make_pair(cluster_key, g4hit)) !=
1278  {
1279  return _cache_get_energy_contribution_g4hit[make_pair(cluster_key, g4hit)];
1280  }
1281 
1282  // this is a fairly simple existance check right now, but might be more
1283  // complex in the future, so this is here mostly as future-proofing.
1284 
1285  float energy = 0.0;
1286  std::set<PHG4Hit*> g4hits = all_truth_hits(cluster_key);
1287  for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
1288  iter != g4hits.end();
1289  ++iter)
1290  {
1291  PHG4Hit* candidate = *iter;
1292  if (candidate->get_hit_id() != g4hit->get_hit_id()) continue;
1293  energy += candidate->get_edep();
1294  }
1295 
1296  if (_do_cache) _cache_get_energy_contribution_g4hit.insert(make_pair(make_pair(cluster_key, g4hit), energy));
1297 
1298  return energy;
1299 }
1300 
1302 {
1303  // need things off of the DST...
1304 
1305  _clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
1306  if(!_clustermap) {
1307  _clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
1308  }
1309  else
1310  {
1313  if(_clustermap->size() == 0)
1314  _clustermap = findNode::getClass<TrkrClusterContainer>(topNode,"TRKR_CLUSTER");
1315  }
1316 
1317 
1318  _hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
1319  _cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
1320  _hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode,"TRKR_HITTRUTHASSOC");
1321  _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1322 
1323  _g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
1324  _g4hits_intt = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_INTT");
1325  _g4hits_mvtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
1326  _g4hits_mms = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MICROMEGAS");
1327  _surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode, "ActsSurfaceMaps");
1328  _tgeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
1329 
1330  return;
1331 }
1332 
1334 {
1335  ActsTransformations transformer;
1336  // loop over all the clusters
1337  auto hitsetrange = _hitsets->getHitSets(TrkrDefs::TrkrId::inttId);
1338  for (auto hitsetitr = hitsetrange.first;
1339  hitsetitr != hitsetrange.second;
1340  ++hitsetitr){
1341  auto range = _clustermap->getClusters(hitsetitr->first);
1342  for( auto iter = range.first; iter != range.second; ++iter ){
1343  TrkrDefs::cluskey cluster_key = iter->first;
1344  unsigned int ilayer = TrkrDefs::getLayer(cluster_key);
1345  TrkrCluster *cluster = iter->second;
1346  auto glob = transformer.getGlobalPosition(cluster, _surfmaps,_tgeometry);
1347  float clus_phi = atan2(glob(1), glob(0));
1348 
1349  multimap<unsigned int, innerMap>::iterator it = _clusters_per_layer.find(ilayer);
1350  if (it == _clusters_per_layer.end())
1351  {
1352  it = _clusters_per_layer.insert(make_pair(ilayer, innerMap()));
1353  }
1354  it->second.insert(make_pair(clus_phi, cluster_key));
1355 
1356  //address wrapping along +/-PI by filling larger area of the map
1357  if (clus_phi - (-M_PI) < _clusters_searching_window) it->second.insert(make_pair(clus_phi + 2 * M_PI, cluster_key));
1358  if (M_PI - clus_phi < _clusters_searching_window) it->second.insert(make_pair(clus_phi - 2 * M_PI, cluster_key));
1359  }
1360  }
1361  return;
1362 }
1363 
1365 {
1366  if (_strict)
1367  assert(_clustermap);
1368  else if (!_clustermap)
1369  return false;
1370 
1371  if (_strict)
1372  assert(_truthinfo);
1373  else if (!_truthinfo)
1374  return false;
1375 
1376  return true;
1377 }
1378 
1380 {
1381  if (x != 0.0f)
1382  {
1383  if (fabsf(x) > fabsf(y))
1384  {
1385  const float z = y / x;
1386  if (x > 0.0)
1387  {
1388  // atan2(y,x) = atan(y/x) if x > 0
1389  return fast_approx_atan2(z);
1390  }
1391  else if (y >= 0.0)
1392  {
1393  // atan2(y,x) = atan(y/x) + PI if x < 0, y >= 0
1394  return fast_approx_atan2(z) + M_PI;
1395  }
1396  else
1397  {
1398  // atan2(y,x) = atan(y/x) - PI if x < 0, y < 0
1399  return fast_approx_atan2(z) - M_PI;
1400  }
1401  }
1402  else // Use property atan(y/x) = PI/2 - atan(x/y) if |y/x| > 1.
1403  {
1404  const float z = x / y;
1405  if (y > 0.0)
1406  {
1407  // atan2(y,x) = PI/2 - atan(x/y) if |y/x| > 1, y > 0
1408  return -fast_approx_atan2(z) + M_PI_2;
1409  }
1410  else
1411  {
1412  // atan2(y,x) = -PI/2 - atan(x/y) if |y/x| > 1, y < 0
1413  return -fast_approx_atan2(z) - M_PI_2;
1414  }
1415  }
1416  }
1417  else
1418  {
1419  if (y > 0.0f) // x = 0, y > 0
1420  {
1421  return M_PI_2;
1422  }
1423  else if (y < 0.0f) // x = 0, y < 0
1424  {
1425  return -M_PI_2;
1426  }
1427  }
1428  return 0.0f; // x,y = 0. Could return NaN instead.
1429 }
1430 
1432 {
1433  // Polynomial approximating arctangenet on the range -1,1.
1434  // Max error < 0.005 (or 0.29 degrees)
1435  const float n1 = 0.97239411f;
1436  const float n2 = -0.19194795f;
1437  return (n1 + n2 * z * z) * z;
1438 }