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