g4simulation/g4eval/SvtxClusterEval.cc

0001 #include "SvtxClusterEval.h"
0002
0003 #include "SvtxHitEval.h"
0004 #include "SvtxTruthEval.h"
0005
0006 #include <trackbase/TrkrCluster.h>
0007 #include <trackbase/TrkrClusterContainer.h>
0008 #include <trackbase/TrkrClusterHitAssoc.h>
0009 #include <trackbase/TrkrDefs.h>
0010 #include <trackbase/TrkrHitSet.h>
0011 #include <trackbase/TrkrHitTruthAssoc.h>
0012
0013 #include <trackbase_historic/ActsTransformations.h>
0014
0015 #include <g4main/PHG4Hit.h>
0016 #include <g4main/PHG4HitContainer.h>
0017 #include <g4main/PHG4HitDefs.h>
0018 #include <g4main/PHG4Particle.h>
0019 #include <g4main/PHG4TruthInfoContainer.h>
0020 #include <g4main/PHG4VtxPoint.h>
0021
0022 #include <phool/PHTimer.h>
0023 #include <phool/getClass.h>
0024
0025 #include <TVector3.h>
0026
0027 #include <cassert>
0028 #include <cfloat>
0029 #include <cmath>
0030 #include <iostream>  // for operator<<, basic_ostream
0031 #include <map>
0032 #include <set>
0033
0034 SvtxClusterEval::SvtxClusterEval(PHCompositeNode* topNode)
0035   : _hiteval(topNode)
0036 {
0037   get_node_pointers(topNode);
0038 }
0039
0040 SvtxClusterEval::~SvtxClusterEval()
0041 {
0042   if (_verbosity > 0)
0043   {
0044     if ((_errors > 0) || (_verbosity > 1))
0045     {
0046       std::cout << "SvtxClusterEval::~SvtxClusterEval() - Error Count: " << _errors << std::endl;
0047     }
0048   }
0049 }
0050
0051 void SvtxClusterEval::next_event(PHCompositeNode* topNode)
0052 {
0053   _cache_all_truth_hits.clear();
0054   _cache_all_truth_clusters.clear();
0055   _cache_max_truth_hit_by_energy.clear();
0056   _cache_max_truth_cluster_by_energy.clear();
0057   _cache_all_truth_particles.clear();
0058   _cache_max_truth_particle_by_energy.clear();
0059   _cache_max_truth_particle_by_cluster_energy.clear();
0060   _cache_all_clusters_from_particle.clear();
0061   _cache_all_clusters_from_g4hit.clear();
0062   _cache_best_cluster_from_g4hit.clear();
0063   _cache_get_energy_contribution_g4particle.clear();
0064   _cache_get_energy_contribution_g4hit.clear();
0065   _cache_best_cluster_from_gtrackid_layer.clear();
0066   _clusters_per_layer.clear();
0067   //  _g4hits_per_layer.clear();
0068   _hiteval.next_event(topNode);
0069
0070   get_node_pointers(topNode);
0071 }
0072
0073 std::map<TrkrDefs::cluskey, std::shared_ptr<TrkrCluster>> SvtxClusterEval::all_truth_clusters(TrkrDefs::cluskey cluster_key)
0074 {
0075   if (_do_cache)
0076   {
0077     const auto iter = _cache_all_truth_clusters.find(cluster_key);
0078     if (iter != _cache_all_truth_clusters.end())
0079     {
0080       return iter->second;
0081     }
0082   }
0083
0084   std::map<TrkrDefs::cluskey, std::shared_ptr<TrkrCluster>> truth_clusters;
0085
0086   unsigned int cluster_layer = TrkrDefs::getLayer(cluster_key);
0087   std::set<PHG4Particle*> particles = all_truth_particles(cluster_key);
0088   for (auto* particle : particles)
0089   {
0090     for (const auto& [ckey, cluster] : get_truth_eval()->all_truth_clusters(particle))
0091     {
0092       if (TrkrDefs::getLayer(ckey) == cluster_layer)
0093       {
0094         truth_clusters.insert(std::make_pair(ckey, cluster));
0095       }
0096     }
0097   }
0098
0099   return _cache_all_truth_clusters.insert(std::make_pair(cluster_key, truth_clusters)).first->second;
0100 }
0101
0102 std::pair<TrkrDefs::cluskey, std::shared_ptr<TrkrCluster>> SvtxClusterEval::max_truth_cluster_by_energy(TrkrDefs::cluskey cluster_key)
0103 {
0104   if (_do_cache)
0105   {
0106     const auto iter = _cache_max_truth_cluster_by_energy.find(cluster_key);
0107     if (iter != _cache_max_truth_cluster_by_energy.end())
0108     {
0109       return iter->second;
0110     }
0111   }
0112
0113   unsigned int cluster_layer = TrkrDefs::getLayer(cluster_key);
0114
0115   PHG4Particle* max_particle = max_truth_particle_by_cluster_energy(cluster_key);
0116   if (!max_particle)
0117   {
0118     return std::make_pair(0, nullptr);
0119   }
0120
0121   if (_verbosity > 0)
0122   {
0123     std::cout << "         max truth particle by cluster energy has  trackID  " << max_particle->get_track_id() << std::endl;
0124   }
0125
0126   TrkrCluster* reco_cluster = _clustermap->findCluster(cluster_key);
0127
0128   auto global = getGlobalPosition(cluster_key, reco_cluster);
0129   double reco_x = global(0);
0130   double reco_y = global(1);
0131   double reco_z = global(2);
0132   double r = sqrt(reco_x * reco_x + reco_y * reco_y);
0133   // double reco_rphi = r*fast_approx_atan2(reco_y, reco_x);
0134   double reco_rphi = r * atan2(reco_y, reco_x);
0135
0136   const std::map<TrkrDefs::cluskey, std::shared_ptr<TrkrCluster>> gclusters = get_truth_eval()->all_truth_clusters(max_particle);
0137   for (const auto& [ckey, candidate_truth_cluster] : gclusters)
0138   {
0139     if (TrkrDefs::getLayer(ckey) != cluster_layer)
0140     {
0141       continue;
0142     }
0143
0144     double gx = candidate_truth_cluster->getX();
0145     double gy = candidate_truth_cluster->getY();
0146     double gz = candidate_truth_cluster->getZ();
0147     double gr = sqrt(gx * gx + gy * gy);
0148     double grphi = gr * atan2(gy, gx);
0149     // double grphi = gr*fast_approx_atan2(gy, gx);
0150
0151     // Find the difference in position from the reco cluster
0152     double dz = reco_z - gz;
0153     double drphi = reco_rphi - grphi;
0154
0155     // approximate 4 sigmas cut
0156     if (cluster_layer > 6 && cluster_layer < 23)
0157     {
0158       if (fabs(drphi) < 4.0 * sig_tpc_rphi_inner &&
0159           fabs(dz) < 4.0 * sig_tpc_z)
0160       {
0161         return std::make_pair(ckey, candidate_truth_cluster);
0162       }
0163     }
0164     if (cluster_layer > 22 && cluster_layer < 39)
0165     {
0166       if (fabs(drphi) < 4.0 * sig_tpc_rphi_mid &&
0167           fabs(dz) < 4.0 * sig_tpc_z)
0168       {
0169         return std::make_pair(ckey, candidate_truth_cluster);
0170       }
0171     }
0172     if (cluster_layer > 38 && cluster_layer < 55)
0173     {
0174       if (fabs(drphi) < 4.0 * sig_tpc_rphi_outer &&
0175           fabs(dz) < 4.0 * sig_tpc_z)
0176       {
0177         return std::make_pair(ckey, candidate_truth_cluster);
0178       }
0179     }
0180     else if (cluster_layer < 3)
0181     {
0182       if (fabs(drphi) < 4.0 * sig_mvtx_rphi &&
0183           fabs(dz) < 4.0 * sig_mvtx_z)
0184       {
0185         return std::make_pair(ckey, candidate_truth_cluster);
0186       }
0187     }
0188     else if (cluster_layer == 55)
0189     {
0190       if (fabs(drphi) < 4.0 * sig_mms_rphi_55)
0191       {
0192         return std::make_pair(ckey, candidate_truth_cluster);
0193       }
0194     }
0195     else if (cluster_layer == 56)
0196     {
0197       if (fabs(dz) < 4.0 * sig_mms_z_56)
0198       {
0199         return std::make_pair(ckey, candidate_truth_cluster);
0200       }
0201     }
0202     else
0203     {
0204       if (fabs(drphi) < 4.0 * sig_intt_rphi &&
0205           fabs(dz) < range_intt_z)
0206       {
0207         return std::make_pair(ckey, candidate_truth_cluster);
0208       }
0209     }
0210   }
0211
0212   return std::make_pair(0, nullptr);
0213 }
0214
0215 std::pair<TrkrDefs::cluskey, TrkrCluster*> SvtxClusterEval::reco_cluster_from_truth_cluster(TrkrDefs::cluskey ckey, const std::shared_ptr<TrkrCluster>& gclus)
0216 {
0217   if (_do_cache)
0218   {
0219     /* this does not work. Cache is not filled in the code below, so always remains empty */
0220     const auto iter = _cache_reco_cluster_from_truth_cluster.find(gclus);
0221     if (iter != _cache_reco_cluster_from_truth_cluster.end())
0222     {
0223       return iter->second;
0224     }
0225   }
0226
0227   double gx = gclus->getX();
0228   double gy = gclus->getY();
0229   double gz = gclus->getZ();
0230   double gr = sqrt(gx * gx + gy * gy);
0231   double grphi = gr * atan2(gy, gx);
0232   // double grphi = gr*fast_approx_atan2(gy, gx);
0233
0234   unsigned int truth_layer = TrkrDefs::getLayer(ckey);
0235
0236   std::set<TrkrDefs::cluskey> reco_cluskeys;
0237   std::set<PHG4Hit*> contributing_hits = get_truth_eval()->get_truth_hits_from_truth_cluster(ckey);
0238   for (auto* cont_g4hit : contributing_hits)
0239   {
0240     std::set<TrkrDefs::cluskey> cluskeys = all_clusters_from(cont_g4hit);  // this returns clusters from this hit in any layer using TrkrAssoc maps
0241
0242     if (_verbosity > 0)
0243     {
0244       std::cout << "       contributing g4hitID " << cont_g4hit->get_hit_id() << " g4trackID " << cont_g4hit->get_trkid() << std::endl;
0245     }
0246
0247     for (unsigned long iter : cluskeys)
0248     {
0249       unsigned int clus_layer = TrkrDefs::getLayer(iter);
0250       if (clus_layer != truth_layer)
0251       {
0252         continue;
0253       }
0254
0255       reco_cluskeys.insert(iter);
0256     }
0257   }
0258
0259   unsigned int nreco = reco_cluskeys.size();
0260   if (nreco > 0)
0261   {
0262     // Find a matching reco cluster with position inside 4 sigmas, and replace reco_cluskey
0263     for (const auto& this_ckey : reco_cluskeys)
0264     {
0265       // get the cluster
0266       TrkrCluster* this_cluster = _clustermap->findCluster(this_ckey);
0267       auto global = getGlobalPosition(this_ckey, this_cluster);
0268       double this_x = global(0);
0269       double this_y = global(1);
0270       double this_z = global(2);
0271       double this_rphi = gr * atan2(this_y, this_x);
0272       // double this_rphi = gr*fast_approx_atan2(this_y, this_x);
0273
0274       // Find the difference in position from the g4cluster
0275       double dz = this_z - gz;
0276       double drphi = this_rphi - grphi;
0277
0278       // approximate 4 sigmas cut
0279       if (truth_layer > 6 && truth_layer < 23)
0280       {
0281         if (fabs(drphi) < 4.0 * sig_tpc_rphi_inner &&
0282             fabs(dz) < 4.0 * sig_tpc_z)
0283         {
0284           return std::make_pair(this_ckey, this_cluster);
0285         }
0286       }
0287       if (truth_layer > 22 && truth_layer < 39)
0288       {
0289         if (fabs(drphi) < 4.0 * sig_tpc_rphi_mid &&
0290             fabs(dz) < 4.0 * sig_tpc_z)
0291         {
0292           return std::make_pair(this_ckey, this_cluster);
0293         }
0294       }
0295       if (truth_layer > 38 && truth_layer < 55)
0296       {
0297         if (fabs(drphi) < 4.0 * sig_tpc_rphi_outer &&
0298             fabs(dz) < 4.0 * sig_tpc_z)
0299         {
0300           return std::make_pair(this_ckey, this_cluster);
0301         }
0302       }
0303       else if (truth_layer < 3)
0304       {
0305         if (fabs(drphi) < 4.0 * sig_mvtx_rphi &&
0306             fabs(dz) < 4.0 * sig_mvtx_z)
0307         {
0308           return std::make_pair(this_ckey, this_cluster);
0309         }
0310       }
0311       else if (truth_layer == 55)
0312       {
0313         if (fabs(drphi) < 4.0 * sig_mms_rphi_55)
0314         {
0315           return std::make_pair(this_ckey, this_cluster);
0316         }
0317       }
0318       else if (truth_layer == 56)
0319       {
0320         if (fabs(dz) < 4.0 * sig_mms_z_56)
0321         {
0322           return std::make_pair(this_ckey, this_cluster);
0323         }
0324       }
0325       else
0326       {
0327         if (fabs(drphi) < 4.0 * sig_intt_rphi &&
0328             fabs(dz) < range_intt_z)
0329         {
0330           return std::make_pair(this_ckey, this_cluster);
0331         }
0332       }
0333     }
0334   }
0335
0336   return std::make_pair(0, nullptr);
0337 }
0338
0339 std::set<PHG4Hit*> SvtxClusterEval::all_truth_hits(TrkrDefs::cluskey cluster_key)
0340 {
0341   if (!has_node_pointers())
0342   {
0343     ++_errors;
0344     return std::set<PHG4Hit*>();
0345   }
0346
0347   if (_do_cache)
0348   {
0349     std::map<TrkrDefs::cluskey, std::set<PHG4Hit*>>::iterator iter =
0350         _cache_all_truth_hits.find(cluster_key);
0351     if (iter != _cache_all_truth_hits.end())
0352     {
0353       return iter->second;
0354     }
0355   }
0356
0357   std::set<PHG4Hit*> truth_hits;
0358
0359   // get all truth hits for this cluster
0360   //_cluster_hit_map->identify();
0361   std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
0362       hitrange = _cluster_hit_map->getHits(cluster_key);  // returns range of pairs {cluster key, hit key} for this cluskey
0363
0364   for (std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
0365            clushititer = hitrange.first;
0366        clushititer != hitrange.second; ++clushititer)
0367   {
0368     TrkrDefs::hitkey hitkey = clushititer->second;
0369     // TrkrHitTruthAssoc uses a map with (hitsetkey, std::pair(hitkey, g4hitkey)) - get the hitsetkey from the cluskey
0370     TrkrDefs::hitsetkey hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
0371
0372     // get all of the g4hits for this hitkey
0373     std::multimap<TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype>> temp_map;
0374     _hit_truth_map->getG4Hits(hitsetkey, hitkey, temp_map);
0375     // returns pairs (hitsetkey, std::pair(hitkey, g4hitkey)) for this hitkey only
0376
0377     for (auto& htiter : temp_map)
0378     {
0379       // extract the g4 hit key here and add the hits to the set
0380       PHG4HitDefs::keytype g4hitkey = htiter.second.second;
0381       PHG4Hit* g4hit = nullptr;
0382       unsigned int trkrid = TrkrDefs::getTrkrId(hitsetkey);
0383       switch (trkrid)
0384       {
0385       case TrkrDefs::tpcId:
0386         g4hit = _g4hits_tpc->findHit(g4hitkey);
0387         break;
0388       case TrkrDefs::inttId:
0389         g4hit = _g4hits_intt->findHit(g4hitkey);
0390         break;
0391       case TrkrDefs::mvtxId:
0392         g4hit = _g4hits_mvtx->findHit(g4hitkey);
0393         break;
0394       case TrkrDefs::micromegasId:
0395         g4hit = _g4hits_mms->findHit(g4hitkey);
0396         break;
0397       default:
0398         break;
0399       }
0400       if (g4hit)
0401       {
0402         truth_hits.insert(g4hit);
0403       }
0404     }  // end loop over g4hits associated with hitsetkey and hitkey
0405   }  // end loop over hits associated with cluskey
0406
0407   if (_do_cache)
0408   {
0409     _cache_all_truth_hits.insert(std::make_pair(cluster_key, truth_hits));
0410   }
0411
0412   return truth_hits;
0413 }
0414
0415 PHG4Hit* SvtxClusterEval::all_truth_hits_by_nhit(TrkrDefs::cluskey cluster_key)
0416 {
0417   if (!has_node_pointers())
0418   {
0419     ++_errors;
0420     return nullptr;
0421   }
0422
0423   //  if (_strict)
0424   //    {
0425   //      assert(cluster_key);
0426   //    }
0427   //  else if (!cluster_key)
0428   //    {
0429   //      ++_errors;
0430   //      return std::set<PHG4Hit*>();
0431   //    }
0432   /*
0433   if (_do_cache)
0434     {
0435       std::map<TrkrDefs::cluskey, std::set<PHG4Hit*> >::iterator iter =
0436         _cache_all_truth_hits.find(cluster_key);
0437       if (iter != _cache_all_truth_hits.end())
0438         {
0439           return iter->second;
0440         }
0441     }
0442   */
0443   TrkrCluster* cluster = _clustermap->findCluster(cluster_key);
0444   auto glob = getGlobalPosition(cluster_key, cluster);
0445   TVector3 cvec(glob(0), glob(1), glob(2));
0446   unsigned int layer = TrkrDefs::getLayer(cluster_key);
0447   std::set<PHG4Hit*> truth_hits;
0448
0449   std::multimap<PHG4HitDefs::keytype, TrkrDefs::hitkey> g4keyperhit;
0450   std::vector<PHG4HitDefs::keytype> g4hitkeys;
0451   // get all truth hits for this cluster
0452   //_cluster_hit_map->identify();
0453   TrkrDefs::hitsetkey hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
0454
0455   std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
0456       hitrange = _cluster_hit_map->getHits(cluster_key);  // returns range of pairs {cluster key, hit key} for this cluskey
0457   for (std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
0458            clushititer = hitrange.first;
0459        clushititer != hitrange.second; ++clushititer)
0460   {
0461     TrkrDefs::hitkey hitkey = clushititer->second;
0462     // TrkrHitTruthAssoc uses a map with (hitsetkey, std::pair(hitkey, g4hitkey)) - get the hitsetkey from the cluskey
0463
0464     // get all of the g4hits for this hitkey
0465     std::multimap<TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype>> temp_map;
0466     _hit_truth_map->getG4Hits(hitsetkey, hitkey, temp_map);  // returns pairs (hitsetkey, std::pair(hitkey, g4hitkey)) for this hitkey only
0467     for (auto& htiter : temp_map)
0468     {
0469       // extract the g4 hit key here and add the hits to the set
0470       PHG4HitDefs::keytype g4hitkey = htiter.second.second;
0471       if (_verbosity > 2)
0472       {
0473         std::cout << " g4key:  " << g4hitkey << " layer: " << layer << std::endl;
0474       }
0475       TrkrDefs::hitkey local_hitkey = htiter.second.first;
0476       /*      if(layer>=7){
0477         PHG4Hit *match_g4hit = _g4hits_tpc->findHit(g4hitkey);
0478         if(layer != match_g4hit->get_layer() ) continue;
0479         }
0480       */
0481       g4keyperhit.insert(std::pair<PHG4HitDefs::keytype, TrkrDefs::hitkey>(g4hitkey, local_hitkey));
0482       std::vector<PHG4HitDefs::keytype>::iterator itg4keys = find(g4hitkeys.begin(), g4hitkeys.end(), g4hitkey);
0483       if (itg4keys == g4hitkeys.end())
0484       {
0485         g4hitkeys.push_back(g4hitkey);
0486       }
0487     }  // end loop over g4hits associated with hitsetkey and hitkey
0488   }  // end loop over hits associated with cluskey
0489
0490   //  if (_do_cache) _cache_all_truth_hits.insert(std::make_pair(cluster_key, truth_hits));
0491   PHG4HitDefs::keytype max_key = 0;
0492   unsigned int n_max = 0;
0493
0494   if (g4hitkeys.size() == 1)
0495   {
0496     std::vector<PHG4HitDefs::keytype>::iterator it = g4hitkeys.begin();
0497     max_key = *it;
0498   }
0499   else
0500   {
0501     for (unsigned long long& g4hitkey : g4hitkeys)
0502     {
0503       unsigned int ng4hit = g4keyperhit.count(g4hitkey);
0504       PHG4Hit* this_g4hit = _g4hits_tpc->findHit(g4hitkey);
0505
0506       if (layer >= 7)
0507       {  // in tpc
0508         if (this_g4hit != nullptr)
0509         {
0510           unsigned int glayer = this_g4hit->get_layer();
0511           if (layer != glayer)
0512           {
0513             continue;
0514           }
0515
0516           TVector3 vec(this_g4hit->get_avg_x(), this_g4hit->get_avg_y(), this_g4hit->get_avg_z());
0517           // std::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 << std::endl; //<< " keyrec: "<< *it.second << std::endl;
0518         }
0519         /*else{
0520           std::cout << "g4hit == NULL " << std::endl;
0521         }
0522         */
0523       }
0524       if (ng4hit > n_max)
0525       {
0526         max_key = g4hitkey;
0527         n_max = ng4hit;
0528       }
0529     }
0530   }
0531
0532   PHG4Hit* g4hit = nullptr;
0533   unsigned int trkrid = TrkrDefs::getTrkrId(hitsetkey);
0534   switch (trkrid)
0535   {
0536   case TrkrDefs::tpcId:
0537     g4hit = _g4hits_tpc->findHit(max_key);
0538     break;
0539   case TrkrDefs::inttId:
0540     g4hit = _g4hits_intt->findHit(max_key);
0541     break;
0542   case TrkrDefs::mvtxId:
0543     g4hit = _g4hits_mvtx->findHit(max_key);
0544     break;
0545   case TrkrDefs::micromegasId:
0546     g4hit = _g4hits_mms->findHit(max_key);
0547     break;
0548   default:
0549     break;
0550   }
0551   if (g4hit)
0552   {
0553     truth_hits.insert(g4hit);
0554   }
0555
0556   return g4hit;
0557 }
0558
0559 std::pair<int, int> SvtxClusterEval::gtrackid_and_layer_by_nhit(TrkrDefs::cluskey cluster_key)
0560 {
0561   if (!has_node_pointers())
0562   {
0563     ++_errors;
0564     return std::make_pair(0, 0);
0565   }
0566
0567   //  if (_strict)
0568   //    {
0569   //      assert(cluster_key);
0570   //    }
0571   //  else if (!cluster_key)
0572   //    {
0573   //      ++_errors;
0574   //      return std::set<PHG4Hit*>();
0575   //    }
0576   /*
0577   if (_do_cache)
0578     {
0579       std::map<TrkrDefs::cluskey, std::set<PHG4Hit*> >::iterator iter =
0580         _cache_all_truth_hits.find(cluster_key);
0581       if (iter != _cache_all_truth_hits.end())
0582         {
0583           return iter->second;
0584         }
0585     }
0586   */
0587
0588   std::pair<int, int> out_pair;
0589   out_pair.first = 0;
0590   out_pair.second = -1;
0591
0592   TrkrCluster* cluster = _clustermap->findCluster(cluster_key);
0593   auto global = getGlobalPosition(cluster_key, cluster);
0594   TVector3 cvec(global(0), global(1), global(2));
0595   unsigned int layer = TrkrDefs::getLayer(cluster_key);
0596
0597   std::multimap<PHG4HitDefs::keytype, TrkrDefs::hitkey> g4keyperhit;
0598   std::vector<PHG4HitDefs::keytype> g4hitkeys;
0599   // get all truth hits for this cluster
0600   //_cluster_hit_map->identify();
0601   TrkrDefs::hitsetkey hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
0602
0603   std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
0604       hitrange = _cluster_hit_map->getHits(cluster_key);  // returns range of pairs {cluster key, hit key} for this cluskey
0605   for (std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
0606            clushititer = hitrange.first;
0607        clushititer != hitrange.second; ++clushititer)
0608   {
0609     TrkrDefs::hitkey hitkey = clushititer->second;
0610     // TrkrHitTruthAssoc uses a map with (hitsetkey, std::pair(hitkey, g4hitkey)) - get the hitsetkey from the cluskey
0611
0612     // get all of the g4hits for this hitkey
0613     std::multimap<TrkrDefs::hitsetkey, std::pair<TrkrDefs::hitkey, PHG4HitDefs::keytype>> temp_map;
0614     _hit_truth_map->getG4Hits(hitsetkey, hitkey, temp_map);  // returns pairs (hitsetkey, std::pair(hitkey, g4hitkey)) for this hitkey only
0615
0616     for (auto& htiter : temp_map)
0617     {
0618       // extract the g4 hit key here and add the hits to the set
0619       PHG4HitDefs::keytype g4hitkey = htiter.second.second;
0620       if (_verbosity > 2)
0621       {
0622         std::cout << " g4key:  " << g4hitkey << " layer: " << layer << std::endl;
0623       }
0624       TrkrDefs::hitkey local_hitkey = htiter.second.first;
0625       /*      if(layer>=7){
0626         PHG4Hit *match_g4hit = _g4hits_tpc->findHit(g4hitkey);
0627         if(layer != match_g4hit->get_layer() ) continue;
0628         }
0629       */
0630       g4keyperhit.insert(std::pair<PHG4HitDefs::keytype, TrkrDefs::hitkey>(g4hitkey, local_hitkey));
0631       std::vector<PHG4HitDefs::keytype>::iterator itg4keys = find(g4hitkeys.begin(), g4hitkeys.end(), g4hitkey);
0632       if (itg4keys == g4hitkeys.end())
0633       {
0634         g4hitkeys.push_back(g4hitkey);
0635       }
0636     }  // end loop over g4hits associated with hitsetkey and hitkey
0637   }  // end loop over hits associated with cluskey
0638
0639   PHG4HitDefs::keytype max_key = 0;
0640   unsigned int n_max = 0;
0641   if (_verbosity > 2)
0642   {
0643     std::cout << " n matches found: " << g4hitkeys.size() << " phi: " << cvec.Phi() << " z: " << cvec.Z() << " ckey: " << cluster_key << std::endl;
0644   }
0645
0646   if (g4hitkeys.size() == 1)
0647   {
0648     std::vector<PHG4HitDefs::keytype>::iterator it = g4hitkeys.begin();
0649     max_key = *it;
0650   }
0651   else
0652   {
0653     for (unsigned long long& g4hitkey : g4hitkeys)
0654     {
0655       unsigned int ng4hit = g4keyperhit.count(g4hitkey);
0656       PHG4Hit* this_g4hit = _g4hits_tpc->findHit(g4hitkey);
0657
0658       if (layer >= 7)
0659       {  // in tpc
0660         if (this_g4hit != nullptr)
0661         {
0662           unsigned int glayer = this_g4hit->get_layer();
0663           //  if(layer != glayer) continue;
0664
0665           TVector3 vec(this_g4hit->get_avg_x(), this_g4hit->get_avg_y(), this_g4hit->get_avg_z());
0666           if (_verbosity > 2)
0667           {
0668             std::cout << "layer: " << layer << " (" << glayer << ") "
0669                       << " gtrackID: " << this_g4hit->get_trkid() << " novlp: " << ng4hit << " phi: " << vec.Phi() << " z: " << this_g4hit->get_avg_z() << " r: " << vec.Perp() << " keyg4: " << g4hitkey << " cz: " << cluster->getZ() << std::endl;  //<< " keyrec: "<< *it.second << std::endl;
0670           }
0671         }
0672       }
0673       if (ng4hit > n_max)
0674       {
0675         max_key = g4hitkey;
0676         n_max = ng4hit;
0677       }
0678     }
0679   }
0680   if (_verbosity > 2)
0681   {
0682     std::cout << "found in layer: " << layer << " n_max: " << n_max << " max_key: " << max_key << " ckey: " << cluster_key << std::endl;
0683   }
0684   if (max_key != 0)
0685   {
0686     PHG4Hit* g4hit = nullptr;
0687     unsigned int trkrid = TrkrDefs::getTrkrId(hitsetkey);
0688     switch (trkrid)
0689     {
0690     case TrkrDefs::tpcId:
0691       g4hit = _g4hits_tpc->findHit(max_key);
0692       break;
0693     case TrkrDefs::inttId:
0694       g4hit = _g4hits_intt->findHit(max_key);
0695       break;
0696     case TrkrDefs::mvtxId:
0697       g4hit = _g4hits_mvtx->findHit(max_key);
0698       break;
0699     case TrkrDefs::micromegasId:
0700       g4hit = _g4hits_mms->findHit(max_key);
0701       break;
0702     default:
0703       break;
0704     }
0705
0706     // check if we on a looper
0707     PHG4Particle* g4particle = _truthinfo->GetParticle(g4hit->get_trkid());
0708
0709     PHG4VtxPoint* vtx = _truthinfo->GetVtx(g4particle->get_vtx_id());
0710     float vtx_z = vtx->get_z();
0711     float gpx = g4particle->get_px();
0712     float gpy = g4particle->get_py();
0713     float gpz = g4particle->get_pz();
0714     float gpeta = std::numeric_limits<float>::quiet_NaN();
0715
0716     TVector3 gv(gpx, gpy, gpz);
0717     gpeta = gv.Eta();
0718     TVector3 this_vec(g4hit->get_avg_x(),
0719                       g4hit->get_avg_y(),
0720                       g4hit->get_avg_z() - vtx_z);
0721     double deta = std::abs(gpeta - this_vec.Eta());
0722
0723     int is_loop = 0;
0724
0725     if (layer >= 7)
0726     {
0727       //        std::cout << " in tpc " << std::endl;
0728       if (deta > 0.1)
0729       {
0730         is_loop = 1;
0731       }
0732     }
0733
0734     out_pair.first = g4hit->get_trkid();
0735     if (!is_loop)
0736     {
0737       out_pair.second = layer;
0738     }
0739   }
0740   return out_pair;
0741 }
0742
0743 PHG4Hit* SvtxClusterEval::max_truth_hit_by_energy(TrkrDefs::cluskey cluster_key)
0744 {
0745   if (!has_node_pointers())
0746   {
0747     ++_errors;
0748     return nullptr;
0749   }
0750
0751   if (_do_cache)
0752   {
0753     std::map<TrkrDefs::cluskey, PHG4Hit*>::iterator iter =
0754         _cache_max_truth_hit_by_energy.find(cluster_key);
0755     if (iter != _cache_max_truth_hit_by_energy.end())
0756     {
0757       return iter->second;
0758     }
0759   }
0760
0761   std::set<PHG4Hit*> hits = all_truth_hits(cluster_key);
0762   PHG4Hit* max_hit = nullptr;
0763   float max_e = std::numeric_limits<float>::min();
0764   for (auto* hit : hits)
0765   {
0766     if (hit->get_edep() > max_e)
0767     {
0768       max_e = hit->get_edep();
0769       max_hit = hit;
0770     }
0771   }
0772
0773   if (_do_cache)
0774   {
0775     _cache_max_truth_hit_by_energy.insert(std::make_pair(cluster_key, max_hit));
0776   }
0777
0778   return max_hit;
0779 }
0780
0781 std::set<PHG4Particle*> SvtxClusterEval::all_truth_particles(TrkrDefs::cluskey cluster_key)
0782 {
0783   if (!has_node_pointers())
0784   {
0785     ++_errors;
0786     return std::set<PHG4Particle*>();
0787   }
0788
0789   if (_do_cache)
0790   {
0791     std::map<TrkrDefs::cluskey, std::set<PHG4Particle*>>::iterator iter =
0792         _cache_all_truth_particles.find(cluster_key);
0793     if (iter != _cache_all_truth_particles.end())
0794     {
0795       return iter->second;
0796     }
0797   }
0798
0799   std::set<PHG4Particle*> truth_particles;
0800
0801   std::set<PHG4Hit*> g4hits = all_truth_hits(cluster_key);
0802
0803   for (auto* hit : g4hits)
0804   {
0805     PHG4Particle* particle = get_truth_eval()->get_particle(hit);
0806     // std::cout << "cluster key " << cluster_key << " has hit " << hit->get_hit_id() << " and has particle " << particle->get_track_id() << std::endl;
0807
0808     if (_strict)
0809     {
0810       assert(particle);
0811     }
0812     else if (!particle)
0813     {
0814       ++_errors;
0815       continue;
0816     }
0817
0818     truth_particles.insert(particle);
0819   }
0820
0821   if (_do_cache)
0822   {
0823     _cache_all_truth_particles.insert(std::make_pair(cluster_key, truth_particles));
0824   }
0825
0826   return truth_particles;
0827 }
0828
0829 PHG4Particle* SvtxClusterEval::max_truth_particle_by_cluster_energy(TrkrDefs::cluskey cluster_key)
0830 {
0831   if (!has_node_pointers())
0832   {
0833     ++_errors;
0834     return nullptr;
0835   }
0836
0837   if (_do_cache)
0838   {
0839     std::map<TrkrDefs::cluskey, PHG4Particle*>::iterator iter =
0840         _cache_max_truth_particle_by_cluster_energy.find(cluster_key);
0841     if (iter != _cache_max_truth_particle_by_cluster_energy.end())
0842     {
0843       return iter->second;
0844     }
0845   }
0846
0847   unsigned int layer = TrkrDefs::getLayer(cluster_key);
0848
0849   // loop over all particles associated with this cluster and
0850   // get the energy contribution for each one, record the max
0851   PHG4Particle* max_particle = nullptr;
0852   float max_e = std::numeric_limits<float>::min();
0853   std::set<PHG4Particle*> particles = all_truth_particles(cluster_key);
0854   for (auto* particle : particles)
0855   {
0856     std::map<TrkrDefs::cluskey, std::shared_ptr<TrkrCluster>> truth_clus = get_truth_eval()->all_truth_clusters(particle);
0857     for (const auto& [ckey, cluster] : truth_clus)
0858     {
0859       if (TrkrDefs::getLayer(ckey) == layer)
0860       {
0861         float e = cluster->getError(0, 0);
0862         if (e > max_e)
0863         {
0864           max_e = e;
0865           max_particle = particle;
0866         }
0867       }
0868     }
0869   }
0870
0871   if (_do_cache)
0872   {
0873     _cache_max_truth_particle_by_cluster_energy.insert(std::make_pair(cluster_key, max_particle));
0874   }
0875
0876   return max_particle;
0877 }
0878
0879 PHG4Particle* SvtxClusterEval::max_truth_particle_by_energy(TrkrDefs::cluskey cluster_key)
0880 {
0881   // Note: this does not quite work correctly for the TPC - it assumes one g4hit per layer
0882   // use max_truth_particle_by_cluster_energy instead
0883
0884   if (!has_node_pointers())
0885   {
0886     ++_errors;
0887     return nullptr;
0888   }
0889
0890   if (_do_cache)
0891   {
0892     std::map<TrkrDefs::cluskey, PHG4Particle*>::iterator iter =
0893         _cache_max_truth_particle_by_energy.find(cluster_key);
0894     if (iter != _cache_max_truth_particle_by_energy.end())
0895     {
0896       return iter->second;
0897     }
0898   }
0899
0900   // loop over all particles associated with this cluster and
0901   // get the energy contribution for each one, record the max
0902   PHG4Particle* max_particle = nullptr;
0903   float max_e = std::numeric_limits<float>::min();
0904   std::set<PHG4Particle*> particles = all_truth_particles(cluster_key);
0905   for (auto* particle : particles)
0906   {
0907     float e = get_energy_contribution(cluster_key, particle);
0908     if (e > max_e)
0909     {
0910       max_e = e;
0911       max_particle = particle;
0912     }
0913   }
0914
0915   if (_do_cache)
0916   {
0917     _cache_max_truth_particle_by_energy.insert(std::make_pair(cluster_key, max_particle));
0918   }
0919
0920   return max_particle;
0921 }
0922
0923 std::set<TrkrDefs::cluskey> SvtxClusterEval::all_clusters_from(PHG4Particle* truthparticle)
0924 {
0925   if (!has_node_pointers())
0926   {
0927     ++_errors;
0928     return std::set<TrkrDefs::cluskey>();
0929   }
0930
0931   if (_strict)
0932   {
0933     assert(truthparticle);
0934   }
0935   else if (!truthparticle)
0936   {
0937     ++_errors;
0938     return std::set<TrkrDefs::cluskey>();
0939   }
0940   // check if cache is filled, if not fill it.
0941   //   if(_cache_all_clusters_from_particle.count(truthparticle)==0){
0942   if (_cache_all_clusters_from_particle.empty())
0943   {
0944     FillRecoClusterFromG4HitCache();
0945   }
0946
0947   if (_do_cache)
0948   {
0949     std::map<PHG4Particle*, std::set<TrkrDefs::cluskey>>::iterator iter =
0950         _cache_all_clusters_from_particle.find(truthparticle);
0951     if (iter != _cache_all_clusters_from_particle.end())
0952     {
0953       return iter->second;
0954     }
0955   }
0956   std::set<TrkrDefs::cluskey> clusters;
0957   return clusters;
0958 }
0959
0960 void SvtxClusterEval::FillRecoClusterFromG4HitCache()
0961 {
0962   auto Mytimer = std::make_unique<PHTimer>("ReCl_timer");
0963   Mytimer->stop();
0964   Mytimer->restart();
0965
0966   std::multimap<PHG4Particle*, TrkrDefs::cluskey> temp_clusters_from_particles;
0967   // loop over all the clusters
0968   for (const auto& hitsetkey : _clustermap->getHitSetKeys())
0969   {
0970     auto range = _clustermap->getClusters(hitsetkey);
0971     for (auto iter = range.first; iter != range.second; ++iter)
0972     {
0973       TrkrDefs::cluskey cluster_key = iter->first;
0974
0975       // loop over all truth particles connected to this cluster
0976       std::set<PHG4Particle*> particles = all_truth_particles(cluster_key);
0977       for (auto* candidate : particles)
0978       {
0979         temp_clusters_from_particles.insert(std::make_pair(candidate, cluster_key));
0980       }
0981     }
0982   }
0983   // Loop over particles and fill cache
0984   PHG4TruthInfoContainer::ConstRange range = _truthinfo->GetParticleRange();
0985   for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0986        iter != range.second; ++iter)
0987   {
0988     PHG4Particle* g4particle = iter->second;
0989     std::set<TrkrDefs::cluskey> clusters;
0990     std::multimap<PHG4Particle*, TrkrDefs::cluskey>::const_iterator lower_bound = temp_clusters_from_particles.lower_bound(g4particle);
0991     std::multimap<PHG4Particle*, TrkrDefs::cluskey>::const_iterator upper_bound = temp_clusters_from_particles.upper_bound(g4particle);
0992     std::multimap<PHG4Particle*, TrkrDefs::cluskey>::const_iterator cfp_iter;
0993     for (cfp_iter = lower_bound; cfp_iter != upper_bound; ++cfp_iter)
0994     {
0995       TrkrDefs::cluskey cluster_key = cfp_iter->second;
0996       clusters.insert(cluster_key);
0997     }
0998     _cache_all_clusters_from_particle.insert(std::make_pair(g4particle, clusters));
0999   }
1000
1001   Mytimer->stop();
1002 }
1003
1004 std::set<TrkrDefs::cluskey> SvtxClusterEval::all_clusters_from(PHG4Hit* truthhit)
1005 {
1006   if (!has_node_pointers())
1007   {
1008     ++_errors;
1009     return std::set<TrkrDefs::cluskey>();
1010   }
1011
1012   if (_strict)
1013   {
1014     assert(truthhit);
1015   }
1016   else if (!truthhit)
1017   {
1018     ++_errors;
1019     return std::set<TrkrDefs::cluskey>();
1020   }
1021
1022   // one time, fill cache of g4hit/cluster pairs
1023   if (_cache_all_clusters_from_g4hit.empty())
1024   {
1025     // make a map of truthhit, cluster_key inside this loop
1026     std::multimap<PHG4HitDefs::keytype, TrkrDefs::cluskey> truth_cluster_map;
1027     std::set<PHG4HitDefs::keytype> all_g4hits_set;
1028     std::map<PHG4HitDefs::keytype, PHG4Hit*> all_g4hits_map;
1029
1030     // get all reco clusters
1031     if (_verbosity > 1)
1032     {
1033       std::cout << "all_clusters_from_g4hit: list all reco clusters " << std::endl;
1034     }
1035
1036     for (const auto& hitsetkey : _clustermap->getHitSetKeys())
1037     {
1038       auto range = _clustermap->getClusters(hitsetkey);
1039       for (auto iter = range.first; iter != range.second; ++iter)
1040       {
1041         TrkrDefs::cluskey cluster_key = iter->first;
1042         int layer = TrkrDefs::getLayer(cluster_key);
1043         TrkrCluster* clus = iter->second;
1044         if (_verbosity > 1)
1045         {
1046           std::cout << " layer " << layer << " cluster_key " << cluster_key << " adc " << clus->getAdc()
1047                     << " localx " << clus->getLocalX()
1048                     << " localy " << clus->getLocalY()
1049                     << std::endl;
1050           std::cout << "  associated hits:";
1051           std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
1052               hitrange = _cluster_hit_map->getHits(cluster_key);  // returns range of pairs {cluster key, hit key} for this cluskey
1053           for (std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
1054                    clushititer = hitrange.first;
1055                clushititer != hitrange.second; ++clushititer)
1056           {
1057             TrkrDefs::hitkey hitkey = clushititer->second;
1058             std::cout << " " << hitkey;
1059           }
1060           std::cout << std::endl;
1061         }
1062
1063         // the returned truth hits were obtained from TrkrAssoc maps
1064         std::set<PHG4Hit*> hits = all_truth_hits(cluster_key);
1065         for (auto* candidate : hits)
1066         {
1067           PHG4HitDefs::keytype g4hit_key = candidate->get_hit_id();
1068
1069           if (_verbosity > 5)
1070           {
1071             int gtrackID = candidate->get_trkid();
1072             std::cout << "   adding cluster with cluster_key " << cluster_key << " g4hit with g4hit_key " << g4hit_key
1073                       << " gtrackID " << gtrackID
1074                       << std::endl;
1075           }
1076
1077           all_g4hits_set.insert(g4hit_key);
1078           all_g4hits_map.insert(std::make_pair(g4hit_key, candidate));
1079
1080           truth_cluster_map.insert(std::make_pair(g4hit_key, cluster_key));
1081         }
1082       }
1083     }
1084
1085     // now fill the cache
1086     // loop over all entries in all_g4hits
1087     for (unsigned long long g4hit_key : all_g4hits_set)
1088     {
1089       if (_verbosity > 5)
1090       {
1091         std::cout << " associations for g4hit_key " << g4hit_key << std::endl;
1092       }
1093
1094       std::map<PHG4HitDefs::keytype, PHG4Hit*>::iterator it = all_g4hits_map.find(g4hit_key);
1095       PHG4Hit* g4hit = it->second;
1096
1097       std::set<TrkrDefs::cluskey> assoc_clusters;
1098
1099       std::pair<std::multimap<PHG4HitDefs::keytype, TrkrDefs::cluskey>::iterator,
1100                 std::multimap<PHG4HitDefs::keytype, TrkrDefs::cluskey>::iterator>
1101           ret = truth_cluster_map.equal_range(g4hit_key);
1102       for (std::multimap<PHG4HitDefs::keytype, TrkrDefs::cluskey>::iterator jter = ret.first; jter != ret.second; ++jter)
1103       {
1104         assoc_clusters.insert(jter->second);
1105
1106         if (_verbosity > 5)
1107         {
1108           std::cout << "             g4hit_key " << g4hit_key << " associated with cluster_key " << jter->second << std::endl;
1109         }
1110       }
1111       // done with this g4hit
1112       _cache_all_clusters_from_g4hit.insert(std::make_pair(g4hit, assoc_clusters));
1113     }
1114   }
1115
1116   // get the clusters
1117   std::set<TrkrDefs::cluskey> clusters;
1118   std::map<PHG4Hit*, std::set<TrkrDefs::cluskey>>::iterator iter =
1119       _cache_all_clusters_from_g4hit.find(truthhit);
1120   if (iter != _cache_all_clusters_from_g4hit.end())
1121   {
1122     return iter->second;
1123   }
1124
1125   if (_clusters_per_layer.empty())
1126   {
1127     fill_cluster_layer_map();
1128   }
1129
1130   return clusters;
1131 }
1132
1133 TrkrDefs::cluskey SvtxClusterEval::best_cluster_by_nhit(int gid, int layer)
1134 {
1135   TrkrDefs::cluskey val = 0;
1136   if (!has_node_pointers())
1137   {
1138     ++_errors;
1139     return val;
1140   }
1141
1142   /*  if (_strict)
1143   {
1144     assert(truthhit);
1145   }
1146   else if (!truthhit)
1147   {
1148     ++_errors;
1149     return val;
1150   }
1151   */
1152   // one time, fill cache of g4hit/cluster pairs
1153   if (_cache_best_cluster_from_gtrackid_layer.empty())
1154   {
1155     // get all reco clusters
1156     // std::cout << "cache size ==0" << std::endl;
1157     if (_verbosity > 1)
1158     {
1159       std::cout << "all_clusters: found # " << _clustermap->size() << std::endl;
1160     }
1161
1162     for (const auto& hitsetkey : _clustermap->getHitSetKeys())
1163     {
1164       auto range = _clustermap->getClusters(hitsetkey);
1165       for (auto iter = range.first; iter != range.second; ++iter)
1166       {
1167         TrkrDefs::cluskey cluster_key = iter->first;
1168         int layer_in = TrkrDefs::getLayer(cluster_key);
1169
1170         if (layer_in < 0)
1171         {
1172           continue;
1173         }
1174         // TrkrCluster *clus = iter->second;
1175
1176         std::pair<int, int> gid_lay = gtrackid_and_layer_by_nhit(cluster_key);
1177
1178         //      std::map<std::pair<int, unsigned int>, TrkrDefs::cluskey>::iterator it_exists;
1179         //      it_exists =
1180         if (!_cache_best_cluster_from_gtrackid_layer.contains(gid_lay))
1181         {
1182           if (gid_lay.second >= 0)
1183           {
1184             _cache_best_cluster_from_gtrackid_layer.insert(std::make_pair(gid_lay, cluster_key));
1185           }
1186         }
1187         else if (_verbosity > 2)
1188         {
1189           std::cout << "found doublematch" << std::endl;
1190           std::cout << "ckey: " << cluster_key << " gtrackID: " << gid_lay.first << " layer: " << gid_lay.second << std::endl;
1191         }
1192       }
1193     }
1194   }
1195
1196   // get the clusters
1197   TrkrDefs::cluskey best_cluster = 0;
1198   //  PHG4Hit*, std::set<TrkrDefs::cluskey> >::iterator iter =
1199
1200   std::map<std::pair<int, int>, TrkrDefs::cluskey>::iterator iter = _cache_best_cluster_from_gtrackid_layer.find(std::make_pair(gid, layer));
1201   if (iter != _cache_best_cluster_from_gtrackid_layer.end())
1202   {
1203     return iter->second;
1204   }
1205
1206   return best_cluster;
1207 }
1208
1209 TrkrDefs::cluskey SvtxClusterEval::best_cluster_from(PHG4Hit* truthhit)
1210 {
1211   if (!has_node_pointers())
1212   {
1213     ++_errors;
1214     return 0;
1215   }
1216
1217   if (_strict)
1218   {
1219     assert(truthhit);
1220   }
1221   else if (!truthhit)
1222   {
1223     ++_errors;
1224     return 0;
1225   }
1226
1227   if (_do_cache)
1228   {
1229     std::map<PHG4Hit*, TrkrDefs::cluskey>::iterator iter =
1230         _cache_best_cluster_from_g4hit.find(truthhit);
1231     if (iter != _cache_best_cluster_from_g4hit.end())
1232     {
1233       return iter->second;
1234     }
1235   }
1236
1237   TrkrDefs::cluskey best_cluster = 0;
1238   float best_energy = 0.0;
1239   std::set<TrkrDefs::cluskey> clusters = all_clusters_from(truthhit);
1240   for (unsigned long cluster_key : clusters)
1241   {
1242     float energy = get_energy_contribution(cluster_key, truthhit);
1243     if (energy > best_energy)
1244     {
1245       best_cluster = cluster_key;
1246       best_energy = energy;
1247     }
1248   }
1249
1250   if (_do_cache)
1251   {
1252     _cache_best_cluster_from_g4hit.insert(std::make_pair(truthhit, best_cluster));
1253   }
1254
1255   return best_cluster;
1256 }
1257
1258 // overlap calculations
1259 float SvtxClusterEval::get_energy_contribution(TrkrDefs::cluskey cluster_key, PHG4Particle* particle)
1260 {
1261   // Note: this does not work correctly for the TPC
1262   // It assumes one g4hit per layer. Use the truth cluster energy instead.
1263
1264   if (!has_node_pointers())
1265   {
1266     ++_errors;
1267     return std::numeric_limits<float>::quiet_NaN();
1268   }
1269
1270   if (_strict)
1271   {
1272     //    assert(cluster_key);
1273     assert(particle);
1274   }
1275   else if (!particle)
1276   {
1277     ++_errors;
1278     return std::numeric_limits<float>::quiet_NaN();
1279   }
1280
1281   if (_do_cache)
1282   {
1283     std::map<std::pair<TrkrDefs::cluskey, PHG4Particle*>, float>::iterator iter =
1284         _cache_get_energy_contribution_g4particle.find(std::make_pair(cluster_key, particle));
1285     if (iter != _cache_get_energy_contribution_g4particle.end())
1286     {
1287       return iter->second;
1288     }
1289   }
1290
1291   float energy = 0.0;
1292   std::set<PHG4Hit*> hits = all_truth_hits(cluster_key);
1293   for (auto* hit : hits)
1294   {
1295     if (get_truth_eval()->is_g4hit_from_particle(hit, particle))
1296     {
1297       energy += hit->get_edep();
1298     }
1299   }
1300
1301   if (_do_cache)
1302   {
1303     _cache_get_energy_contribution_g4particle.insert(std::make_pair(std::make_pair(cluster_key, particle), energy));
1304   }
1305
1306   return energy;
1307 }
1308
1309 float SvtxClusterEval::get_energy_contribution(TrkrDefs::cluskey cluster_key, PHG4Hit* g4hit)
1310 {
1311   if (!has_node_pointers())
1312   {
1313     ++_errors;
1314     return std::numeric_limits<float>::quiet_NaN();
1315   }
1316
1317   if (_strict)
1318   {
1319     //    assert(cluster_key);
1320     assert(g4hit);
1321   }
1322   else if (!g4hit)
1323   {
1324     ++_errors;
1325     return std::numeric_limits<float>::quiet_NaN();
1326   }
1327
1328   if ((_do_cache) &&
1329       (_cache_get_energy_contribution_g4hit.contains(std::make_pair(cluster_key, g4hit))))
1330   {
1331     return _cache_get_energy_contribution_g4hit[std::make_pair(cluster_key, g4hit)];
1332   }
1333
1334   // this is a fairly simple existance check right now, but might be more
1335   // complex in the future, so this is here mostly as future-proofing.
1336
1337   float energy = 0.0;
1338   std::set<PHG4Hit*> g4hits = all_truth_hits(cluster_key);
1339   for (auto* candidate : g4hits)
1340   {
1341     if (candidate->get_hit_id() != g4hit->get_hit_id())
1342     {
1343       continue;
1344     }
1345     energy += candidate->get_edep();
1346   }
1347
1348   if (_do_cache)
1349   {
1350     _cache_get_energy_contribution_g4hit.insert(std::make_pair(std::make_pair(cluster_key, g4hit), energy));
1351   }
1352
1353   return energy;
1354 }
1355
1356 void SvtxClusterEval::get_node_pointers(PHCompositeNode* topNode)
1357 {
1358   // need things off of the DST...
1359
1360   _clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
1361   if (!_clustermap)
1362   {
1363     _clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
1364   }
1365   else
1366   {
1367     /// Need a catch for if the node corrected cluster node exists but hasn't been filled
1368     /// yet, in e.g. the case of truth track seeding
1369     if (_clustermap->size() == 0)
1370     {
1371       _clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
1372     }
1373   }
1374
1375   _cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
1376   _hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
1377   _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1378
1379   _g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
1380   _g4hits_intt = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_INTT");
1381   _g4hits_mvtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
1382   _g4hits_mms = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MICROMEGAS");
1383   _tgeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
1384
1385   return;
1386 }
1387
1388 void SvtxClusterEval::fill_cluster_layer_map()
1389 {
1390   // loop over all the clusters
1391   for (const auto& hitsetkey : _clustermap->getHitSetKeys())
1392   {
1393     auto range = _clustermap->getClusters(hitsetkey);
1394     for (auto iter = range.first; iter != range.second; ++iter)
1395     {
1396       TrkrDefs::cluskey cluster_key = iter->first;
1397       unsigned int ilayer = TrkrDefs::getLayer(cluster_key);
1398       TrkrCluster* cluster = iter->second;
1399       auto glob = getGlobalPosition(cluster_key, cluster);
1400       float clus_phi = atan2(glob(1), glob(0));
1401
1402       std::multimap<unsigned int, innerMap>::iterator it = _clusters_per_layer.find(ilayer);
1403       if (it == _clusters_per_layer.end())
1404       {
1405         it = _clusters_per_layer.insert(std::make_pair(ilayer, innerMap()));
1406       }
1407       it->second.insert(std::make_pair(clus_phi, cluster_key));
1408
1409       // address wrapping along +/-PI by filling larger area of the map
1410       if (clus_phi - (-M_PI) < _clusters_searching_window)
1411       {
1412         it->second.insert(std::make_pair(clus_phi + 2 * M_PI, cluster_key));
1413       }
1414       if (M_PI - clus_phi < _clusters_searching_window)
1415       {
1416         it->second.insert(std::make_pair(clus_phi - 2 * M_PI, cluster_key));
1417       }
1418     }
1419   }
1420   return;
1421 }
1422
1423 bool SvtxClusterEval::has_node_pointers()
1424 {
1425   if (_strict)
1426   {
1427     assert(_clustermap);
1428   }
1429   else if (!_clustermap)
1430   {
1431     return false;
1432   }
1433
1434   if (_strict)
1435   {
1436     assert(_truthinfo);
1437   }
1438   else if (!_truthinfo)
1439   {
1440     return false;
1441   }
1442
1443   return true;
1444 }
1445
1446 float SvtxClusterEval::fast_approx_atan2(float y, float x)
1447 {
1448   if (x != 0.0F)
1449   {
1450     if (std::abs(x) > std::abs(y))
1451     {
1452       const float z = y / x;
1453       if (x > 0.0)
1454       {
1455         // atan2(y,x) = atan(y/x) if x > 0
1456         return fast_approx_atan2(z);
1457       }
1458       if (y >= 0.0)
1459       {
1460         // atan2(y,x) = atan(y/x) + PI if x < 0, y >= 0
1461         return fast_approx_atan2(z) + M_PI;
1462       }
1463
1464       // atan2(y,x) = atan(y/x) - PI if x < 0, y < 0
1465       return fast_approx_atan2(z) - M_PI;
1466     }
1467     const float z = x / y;
1468     if (y > 0.0)
1469     {
1470       // atan2(y,x) = PI/2 - atan(x/y) if |y/x| > 1, y > 0
1471       return -fast_approx_atan2(z) + M_PI_2;
1472     }
1473
1474     // atan2(y,x) = -PI/2 - atan(x/y) if |y/x| > 1, y < 0
1475     return -fast_approx_atan2(z) - M_PI_2;
1476   }
1477
1478   if (y > 0.0F)  // x = 0, y > 0
1479   {
1480     return M_PI_2;
1481   }
1482   if (y < 0.0F)  // x = 0, y < 0
1483   {
1484     return -M_PI_2;
1485   }
1486
1487   return 0.0F;  // x,y = 0. Could return NaN instead.
1488 }
1489
1490 float SvtxClusterEval::fast_approx_atan2(float z)
1491 {
1492   // Polynomial approximating arctangenet on the range -1,1.
1493   // Max error < 0.005 (or 0.29 degrees)
1494   const float n1 = 0.97239411F;
1495   const float n2 = -0.19194795F;
1496   return (n1 + n2 * z * z) * z;
1497 }
1498
1499 Acts::Vector3 SvtxClusterEval::getGlobalPosition(TrkrDefs::cluskey cluster_key, TrkrCluster* cluster)
1500 {
1501   Acts::Vector3 glob;
1502   glob = _tgeometry->getGlobalPosition(cluster_key, cluster);
1503
1504   return glob;
1505 }