37 , _clustermap(nullptr)
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()
67 cout <<
"SvtxClusterEval::~SvtxClusterEval() - Error Count: " <<
_errors << endl;
98 std::map<TrkrDefs::cluskey, std::set<std::shared_ptr<TrkrCluster> > >::iterator iter =
106 std::set<std::shared_ptr<TrkrCluster> > truth_clusters;
111 for (std::set<PHG4Particle*>::iterator iter = particles.begin();
112 iter != particles.end();
118 for (std::map<
unsigned int, std::shared_ptr<TrkrCluster> >::iterator citer = gclusters.begin();
119 citer != gclusters.end();
122 if(citer->first == cluster_layer)
124 truth_clusters.insert(citer->second);
129 return truth_clusters;
136 std::map<TrkrDefs::cluskey, std::shared_ptr<TrkrCluster> >::iterator iter =
144 std::shared_ptr<TrkrCluster> truth_cluster = 0;
150 return truth_cluster;
153 cout <<
" max truth particle by cluster energy has trackID " << max_particle->
get_track_id() << endl;
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);
163 double reco_rphi = r*atan2(reco_y, reco_x);
166 for (std::map<
unsigned int, std::shared_ptr<TrkrCluster> >::iterator citer = gclusters.begin();
167 citer != gclusters.end();
170 if(citer->first == cluster_layer)
172 std::shared_ptr<TrkrCluster> candidate_truth_cluster = citer->second;
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);
182 double dz = reco_z - gz;
183 double drphi = reco_rphi - grphi;
186 if(cluster_layer > 6 && cluster_layer < 23)
191 return candidate_truth_cluster;;
194 if(cluster_layer > 22 && cluster_layer < 39)
199 return candidate_truth_cluster;;
202 if(cluster_layer > 38 && cluster_layer < 55)
207 return candidate_truth_cluster;;
210 else if(cluster_layer < 3)
215 return candidate_truth_cluster;;
218 else if(cluster_layer == 55)
222 return candidate_truth_cluster;;
225 else if(cluster_layer == 56)
229 return candidate_truth_cluster;;
237 return candidate_truth_cluster;;
243 return truth_cluster;
250 std::map<std::shared_ptr<TrkrCluster>,
TrkrCluster* >::iterator iter =
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);
269 std::set<TrkrDefs::cluskey> reco_cluskeys;
271 for(
auto it = contributing_hits.begin();
it != contributing_hits.end(); ++
it)
277 cout <<
" contributing g4hitID " << cont_g4hit->
get_hit_id() <<
" g4trackID " << cont_g4hit->
get_trkid() << endl;
279 for (std::set<TrkrDefs::cluskey>::iterator iter = cluskeys.begin();
280 iter != cluskeys.end();
285 if(clus_layer != truth_layer)
continue;
287 reco_cluskeys.insert(*iter);
291 unsigned int nreco = reco_cluskeys.size();
296 for(std::set<TrkrDefs::cluskey>::iterator
it = reco_cluskeys.begin();
it != reco_cluskeys.end(); ++
it)
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);
308 double dz = this_z - gz;
309 double drphi = this_rphi - grphi;
312 if(truth_layer > 6 && truth_layer < 23)
317 return this_cluster;;
320 if(truth_layer > 22 && truth_layer < 39)
325 return this_cluster;;
328 if(truth_layer > 38 && truth_layer < 55)
333 return this_cluster;;
336 else if(truth_layer < 3)
341 return this_cluster;;
344 else if(truth_layer == 55)
348 return this_cluster;;
351 else if(truth_layer == 56)
355 return this_cluster;;
363 return this_cluster;;
377 return std::set<PHG4Hit*>();
382 std::map<TrkrDefs::cluskey, std::set<PHG4Hit*> >::iterator iter =
390 std::set<PHG4Hit*> truth_hits;
394 std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
397 for(std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
398 clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
405 std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> > temp_map;
409 for(std::multimap<
TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> >::iterator htiter = temp_map.begin(); htiter != temp_map.end(); ++htiter)
424 if( g4hit ) truth_hits.insert(g4hit);
464 TVector3 cvec(glob(0), glob(1), glob(2));
466 std::set<PHG4Hit*> truth_hits;
468 std::multimap<PHG4HitDefs::keytype,TrkrDefs::hitkey> g4keyperhit;
469 std::vector<PHG4HitDefs::keytype> g4hitkeys;
474 std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
476 for(std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
477 clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
483 std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> > temp_map;
485 for(std::multimap<
TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> >::iterator htiter = temp_map.begin(); htiter != temp_map.end(); ++htiter)
490 cout <<
" g4key: " << g4hitkey <<
" layer: " << layer << endl;
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);
505 unsigned int n_max = 0;
507 if(g4hitkeys.size()==1 ){
508 std::vector<PHG4HitDefs::keytype>::iterator
it = g4hitkeys.begin();
511 for(std::vector<PHG4HitDefs::keytype>::iterator
it = g4hitkeys.begin();
it != g4hitkeys.end(); ++
it){
512 unsigned int ng4hit = g4keyperhit.count(*
it);
516 if(this_g4hit!=NULL){
517 unsigned int glayer = this_g4hit->
get_layer();
518 if(layer != glayer)
continue;
546 if( g4hit ) truth_hits.insert(g4hit);
556 return make_pair(0,0);
580 std::pair<int, int> out_pair;
582 out_pair.second = -1;
587 TVector3 cvec(global(0), global(1), global(2));
590 std::multimap<PHG4HitDefs::keytype,TrkrDefs::hitkey> g4keyperhit;
591 std::vector<PHG4HitDefs::keytype> g4hitkeys;
596 std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
598 for(std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
599 clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
605 std::multimap< TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> > temp_map;
608 for(std::multimap<
TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype> >::iterator htiter = temp_map.begin(); htiter != temp_map.end(); ++htiter)
613 cout <<
" g4key: " << g4hitkey <<
" layer: " << layer << endl;
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);
627 unsigned int n_max = 0;
629 cout <<
" n matches found: " << g4hitkeys.size() <<
" phi: " << cvec.Phi() <<
" z: " << cvec.Z() <<
" ckey: " << cluster_key << endl;
631 if(g4hitkeys.size()==1 ){
632 std::vector<PHG4HitDefs::keytype>::iterator
it = g4hitkeys.begin();
635 for(std::vector<PHG4HitDefs::keytype>::iterator
it = g4hitkeys.begin();
it != g4hitkeys.end(); ++
it){
636 unsigned int ng4hit = g4keyperhit.count(*
it);
641 if(this_g4hit!=NULL){
642 unsigned int glayer = this_g4hit->
get_layer();
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;
658 cout <<
"found in layer: " << layer <<
" n_max: " << n_max <<
" max_key: " << max_key <<
" ckey: " << cluster_key << endl;
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();
681 TVector3 gv(gpx, gpy, gpz);
686 double deta = TMath::Abs(gpeta - this_vec.Eta());
692 if(deta>0.1) is_loop = 1;
697 out_pair.second =
layer;
712 std::map<TrkrDefs::cluskey, PHG4Hit*>::iterator iter =
723 for (std::set<PHG4Hit*>::iterator iter = hits.begin();
745 return std::set<PHG4Particle*>();
750 std::map<TrkrDefs::cluskey, std::set<PHG4Particle*> >::iterator iter =
758 std::set<PHG4Particle*> truth_particles;
762 for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
763 iter != g4hits.end();
780 truth_particles.insert(particle);
785 return truth_particles;
798 std::map<TrkrDefs::cluskey, PHG4Particle*>::iterator iter =
813 for (std::set<PHG4Particle*>::iterator iter = particles.begin();
814 iter != particles.end();
819 for(
auto it = truth_clus.begin();
it != truth_clus.end(); ++
it)
821 if(
it->first == layer)
823 float e =
it->second->getError(0,0);
851 std::map<TrkrDefs::cluskey, PHG4Particle*>::iterator iter =
864 for (std::set<PHG4Particle*>::iterator iter = particles.begin();
865 iter != particles.end();
887 return std::set<TrkrDefs::cluskey>();
892 assert(truthparticle);
894 else if (!truthparticle)
897 return std::set<TrkrDefs::cluskey>();
907 std::map<PHG4Particle*, std::set<TrkrDefs::cluskey> >::iterator iter =
914 std::set<TrkrDefs::cluskey> clusters;
919 auto Mytimer = std::make_unique<PHTimer>(
"ReCl_timer");
923 std::multimap<PHG4Particle*, TrkrDefs::cluskey> temp_clusters_from_particles;
926 for (
auto hitsetitr = hitsetrange.first;
927 hitsetitr != hitsetrange.second;
930 for(
auto iter = range.first; iter != range.second; ++iter ){
935 for (std::set<PHG4Particle*>::iterator jter = particles.begin();
936 jter != particles.end();
939 temp_clusters_from_particles.insert(make_pair(candidate, cluster_key));
946 iter != range.second; ++iter){
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){
954 clusters.insert(cluster_key);
968 return std::set<TrkrDefs::cluskey>();
978 return std::set<TrkrDefs::cluskey>();
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;
990 if(
_verbosity > 1) cout <<
"all_clusters_from_g4hit: list all reco clusters " << endl;
993 for (
auto hitsetitr = hitsetrange.first;
994 hitsetitr != hitsetrange.second;
997 for(
auto iter = range.first; iter != range.second; ++iter ){
1003 cout <<
" layer " << layer <<
" cluster_key " << cluster_key <<
" adc " << clus->
getAdc()
1007 cout <<
" associated hits:";
1008 std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
1010 for(std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
1011 clushititer = hitrange.first; clushititer != hitrange.second; ++clushititer)
1021 for (std::set<PHG4Hit*>::iterator jter = hits.begin();
1031 cout <<
" adding cluster with cluster_key " << cluster_key <<
" g4hit with g4hit_key " << g4hit_key
1032 <<
" gtrackID " << gtrackID
1036 all_g4hits_set.insert(g4hit_key);
1037 all_g4hits_map.insert(std::make_pair(g4hit_key, candidate));
1039 truth_cluster_map.insert(std::make_pair(g4hit_key, cluster_key));
1046 for(std::set<PHG4HitDefs::keytype>::iterator iter = all_g4hits_set.begin(); iter != all_g4hits_set.end(); ++iter)
1049 if(
_verbosity > 5) cout <<
" associations for g4hit_key " << g4hit_key << endl;
1051 std::map<PHG4HitDefs::keytype, PHG4Hit*>::iterator
it = all_g4hits_map.find(g4hit_key);
1054 std::set<TrkrDefs::cluskey> assoc_clusters;
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)
1060 assoc_clusters.insert(jter->second);
1062 if(
_verbosity > 5) cout <<
" g4hit_key " << g4hit_key <<
" associated with cluster_key " << jter->second << endl;
1071 std::set<TrkrDefs::cluskey> clusters;
1072 std::map<PHG4Hit*, std::set<TrkrDefs::cluskey> >::iterator iter =
1076 return iter->second;
1115 for (
auto hitsetitr = hitsetrange.first;
1116 hitsetitr != hitsetrange.second;
1119 for(
auto iter = range.first; iter != range.second; ++iter ){
1123 if(layer_in<0)
continue;
1131 if(gid_lay.second >=0)
1135 if(
_verbosity > 2){ cout <<
"found doublematch" << endl;
1136 cout <<
"ckey: " << cluster_key <<
" gtrackID: " << gid_lay.first <<
" layer: " << gid_lay.second << endl;
1150 return iter->second;
1153 return best_cluster;
1176 std::map<PHG4Hit*, TrkrDefs::cluskey>::iterator iter =
1180 return iter->second;
1185 float best_energy = 0.0;
1187 for (std::set<TrkrDefs::cluskey>::iterator iter = clusters.begin();
1188 iter != clusters.end();
1193 if (energy > best_energy)
1195 best_cluster = cluster_key;
1202 return best_cluster;
1222 else if ( !particle)
1230 std::map<std::pair<TrkrDefs::cluskey, PHG4Particle*>,
float>::iterator iter =
1234 return iter->second;
1240 for (std::set<PHG4Hit*>::iterator iter = hits.begin();
1287 for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
1288 iter != g4hits.end();
1305 _clustermap = findNode::getClass<TrkrClusterContainer>(topNode,
"CORRECTED_TRKR_CLUSTER");
1307 _clustermap = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
1314 _clustermap = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
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");
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");
1338 for (
auto hitsetitr = hitsetrange.first;
1339 hitsetitr != hitsetrange.second;
1342 for(
auto iter = range.first; iter != range.second; ++iter ){
1347 float clus_phi = atan2(glob(1), glob(0));
1354 it->second.insert(make_pair(clus_phi, cluster_key));
1383 if (fabsf(x) > fabsf(y))
1385 const float z = y /
x;
1404 const float z = x /
y;
1435 const float n1 = 0.97239411f;
1436 const float n2 = -0.19194795f;
1437 return (n1 + n2 * z * z) *
z;