packages/trackreco/PHGhostRejection.cc

0001 #include "PHGhostRejection.h"
0002
0003 #include "PHGhostRejection.h"
0004
0005 /// Tracking includes
0006
0007 #include <trackbase/TrkrCluster.h>  // for TrkrCluster
0008 #include <trackbase/TrkrClusterContainer.h>
0009 #include <trackbase/TrkrDefs.h>  // for cluskey, getLayer, TrkrId
0010
0011 #include <trackbase_historic/TrackSeed.h>
0012 #include <trackbase_historic/TrackSeed_v2.h>
0013 #include <trackbase_historic/TrackSeedContainer.h>
0014 #include <trackbase_historic/TrackSeedHelper.h>
0015
0016 #include <cmath>     // for sqrt, fabs, atan2, cos
0017 #include <iostream>  // for operator<<, basic_ostream
0018 #include <map>       // for map
0019 #include <set>       // for _Rb_tree_const_iterator
0020 #include <utility>   // for pair, make_pair
0021
0022 //____________________________________________________________________________..
0023 bool PHGhostRejection::cut_from_clusters(int itrack) {
0024   const auto& track = seeds[itrack];
0025   unsigned int nclusters = track.size_cluster_keys();
0026   if (nclusters < _min_clusters)
0027   {
0028     m_rejected[itrack] = true;
0029     if (m_verbosity > 3) {
0030       std::cout << " rejecting track ID (" << ((int)itrack)
0031         <<") because n-clusters(" << ((int)nclusters) <<")" << std::endl;
0032     }
0033     return true;
0034   }
0035   if (_must_span_sectors)
0036   {
0037     // check that there are clusters in at least 2 of the three layers of sectors
0038     bool in_two_sectors = false;
0039     bool in_0 = false;
0040     bool in_1 = false;
0041     bool in_2 = false;
0042     if (m_verbosity > 3)
0043     {
0044       std::cout << " layers in track: ";
0045     }
0046     for (auto key = track.begin_cluster_keys();
0047         key != track.end_cluster_keys();
0048         ++key)
0049     {
0050       unsigned int layer = TrkrDefs::getLayer(*key);
0051       if (m_verbosity > 4)
0052       {
0053         std::cout << ((int) layer) << " ";
0054       }
0055
0056       if (layer < 23)
0057       { in_0 = true; }
0058       else if (layer < 40)
0059       { in_1 = true; }
0060       else
0061       { in_2 = true; }
0062
0063       if ((in_0+in_1+in_2)>1)
0064       {
0065         in_two_sectors = true;
0066         break;
0067       }
0068     }
0069     if (m_verbosity > 3)
0070     { std::cout << std::endl; }
0071     if (!in_two_sectors)
0072     {
0073       m_rejected[itrack] = true;
0074       if (m_verbosity > 1) {
0075         std::cout << " Cutting track ID " << ((int)itrack)
0076           << "  because only has clusters in "
0077           << (in_0 ? "inner sector  (layers <23)"
0078             : in_1 ? "middle sector  (layers 23-39)"
0079             : "outer   (layers >40) ")
0080           << std::endl;
0081       }
0082       return true;
0083     }
0084   }
0085   return false;
0086 }
0087
0088 void PHGhostRejection::find_ghosts(const std::vector<float>& trackChi2)
0089 {
0090   if (m_verbosity > 0)
0091   {
0092     std::cout << "PHGhostRejection beginning track map size " << seeds.size() << std::endl;
0093   }
0094
0095   // cut now pt tracks
0096   if (_min_pt>0.) {
0097     for (unsigned int i=0;i<seeds.size();++i) {
0098       if (seeds[i].get_pt() < _min_pt) {
0099         m_rejected[i] = true;
0100       }
0101     }
0102   }
0103
0104   // Elimate low-interest track, and try to eliminate repeated tracks
0105   std::set<unsigned int> matches_set;
0106   std::multimap<unsigned int, unsigned int> matches;
0107   for (size_t trid1 = 0; trid1 < seeds.size(); ++trid1)
0108   {
0109     if (m_rejected[trid1]) { continue; }
0110     const auto& track1 = seeds[trid1];
0111     const float track1phi = track1.get_phi();
0112
0113     const auto track1_pos = TrackSeedHelper::get_xyz(&track1);
0114     const float track1eta = track1.get_eta();
0115     for (size_t trid2 = trid1+1; trid2 < seeds.size(); ++trid2)
0116     {
0117       if (m_rejected[trid2])
0118       {
0119         continue;
0120       }
0121
0122       const auto& track2 = seeds[trid2];
0123       const auto track2_pos = TrackSeedHelper::get_xyz(&track2);
0124       const float track2eta = track2.get_eta();
0125       auto delta_phi = std::abs(track1phi - track2.get_phi());
0126
0127       if (delta_phi > 2 * M_PI) {
0128         delta_phi = delta_phi - 2*M_PI;
0129       }
0130       if (delta_phi < _phi_cut &&
0131           std::fabs(track1eta - track2eta) < _eta_cut &&
0132           std::fabs(track1_pos.x() - track2_pos.x()) < _x_cut &&
0133           std::fabs(track1_pos.y() - track2_pos.y()) < _y_cut &&
0134           std::fabs(track1_pos.z() - track2_pos.z()) < _z_cut)
0135       {
0136         matches_set.insert(trid1);
0137         matches.insert(std::pair(trid1, trid2));
0138
0139         if (m_verbosity > 1)
0140         {
0141           std::cout << "Found match for tracks " << trid1 << " and " << trid2 << std::endl;
0142         }
0143       }
0144     }
0145   }
0146
0147   for (auto set_it : matches_set)
0148   {
0149     if (m_rejected[set_it]) { continue; } // already rejected
0150     auto match_list = matches.equal_range(set_it);
0151
0152     auto tr1 = seeds[set_it];
0153     double best_qual = trackChi2.at(set_it);
0154     unsigned int best_track = set_it;
0155
0156     if (m_verbosity > 1)
0157     {
0158       std::cout << " ****** start checking track " << set_it << " with best quality " << best_qual << " best_track " << best_track << std::endl;
0159     }
0160
0161     for (auto it = match_list.first; it != match_list.second; ++it)
0162     {
0163       if (m_verbosity > 1)
0164       {
0165         std::cout << "    match of track " << it->first << " to track " << it->second << std::endl;
0166       }
0167
0168       auto tr2 = seeds[it->second];
0169
0170       // Check that these two tracks actually share the same clusters, if not skip this pair
0171       bool is_same_track = checkClusterSharing(tr1, tr2);
0172       if (!is_same_track)
0173       {
0174         continue;
0175       }
0176
0177       // which one has the best quality?
0178       double tr2_qual = trackChi2.at(it->second);
0179       if (m_verbosity > 1)
0180       {
0181         std::cout << "       Compare: best quality " << best_qual << " track 2 quality " << tr2_qual << std::endl;
0182
0183         const auto track1_pos = TrackSeedHelper::get_xyz(&tr1);
0184         std::cout << "       tr1: phi " << tr1.get_phi() << " eta " << tr1.get_eta()
0185                   << " x " << track1_pos.x() << " y " << track1_pos.y() << " z " << track1_pos.z() << std::endl;
0186
0187         const auto track2_pos = TrackSeedHelper::get_xyz(&tr2);
0188         std::cout << "       tr2: phi " << tr2.get_phi() << " eta " << tr2.get_eta()
0189                   << " x " << track2_pos.x() << " y " << track2_pos.y() << " z " << track2_pos.z() << std::endl;
0190       }
0191
0192       if (tr2_qual < best_qual)
0193       {
0194         if (m_verbosity > 1)
0195         {
0196           std::cout << "       --------- Track " << it->second << " has better quality, erase track " << best_track << std::endl;
0197           std::cout << " rejecting track ID " << ((int)best_track) << "  because it is a ghost " << std::endl;
0198         }
0199         m_rejected[best_track] = true;
0200         best_qual = tr2_qual;
0201         best_track = it->second;
0202       }
0203       else
0204       {
0205         if (m_verbosity > 1)
0206         {
0207           std::cout << "       --------- Track " << best_track << " has better quality, erase track " << it->second << std::endl;
0208           std::cout << " rejecting track ID " << ((int)best_track) << "  because it is a ghost " << std::endl;
0209         }
0210         m_rejected[it->second] = true;
0211       }
0212     }
0213     if (m_verbosity > 1)
0214     {
0215       std::cout << " best track " << best_track << " best_qual " << best_qual << std::endl;
0216     }
0217   }
0218
0219   // delete ghost tracks
0220   if (m_verbosity > 1)
0221   {
0222     for (unsigned int it=0 ; it<m_rejected.size(); ++it)
0223     {
0224       if (m_rejected[it]) {
0225         std::cout << " rejecting track ID " << it << std::endl;
0226       }
0227     }
0228   }
0229
0230   if (m_verbosity > 0)
0231   {
0232     int n_ghost = std::count(m_rejected.begin(), m_rejected.end(), true);
0233     std::cout << " Track list sizes: n_init(" << ((int)m_rejected.size()) <<") - n_ghost("<<n_ghost << ") = n_good(" << (m_rejected.size()-n_ghost) << ")" << std::endl;
0234   }
0235 }
0236
0237 // there is no check, at this point, about which is the best chi2 track
0238 bool PHGhostRejection::checkClusterSharing(const TrackSeed& tr1, const TrackSeed& tr2) const
0239 {
0240   // count shared clusters that tr1 and tr2 share many clusters
0241   size_t nclus_tr1 = tr1.size_cluster_keys();
0242   size_t nclus_tr2 = tr2.size_cluster_keys();
0243   size_t n_shared_clus = 0;
0244
0245   for (auto key_tr1 = tr1.begin_cluster_keys();
0246        key_tr1 != tr1.end_cluster_keys();
0247        ++key_tr1)
0248   {
0249     if (tr2.find_cluster_key(*key_tr1) != tr2.end_cluster_keys())
0250     {
0251       ++n_shared_clus;
0252     }
0253   }
0254
0255   if (m_verbosity > 2)
0256   {
0257     std::cout << " N-clusters tr1: " << nclus_tr1 << " N-clusters tr2: " << nclus_tr2 << " N-clusters shared: " << n_shared_clus << std::endl;
0258   }
0259   size_t nreq = 2 * n_shared_clus + 1;
0260   return (nreq > nclus_tr1) || (nreq > nclus_tr2);
0261 }