packages/trackreco/PHTrackCleaner.cc

0001 #include "PHTrackCleaner.h"
0002
0003 #include "PHTrackCleaner.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 #include <trackbase_historic/SvtxTrack.h>  // for SvtxTrack, SvtxTrack::C...
0011 #include <trackbase_historic/SvtxTrackMap.h>
0012 #include <trackbase_historic/TrackSeedContainer.h>
0013
0014 #include <fun4all/Fun4AllReturnCodes.h>
0015
0016 #include <phool/getClass.h>
0017 #include <phool/phool.h>
0018
0019 #include <cmath>     // for sqrt, fabs, atan2, cos
0020 #include <iostream>  // for operator<<, basic_ostream
0021 #include <map>       // for map
0022 #include <set>       // for _Rb_tree_const_iterator
0023 #include <utility>   // for pair, make_pair
0024
0025 //____________________________________________________________________________..
0026 PHTrackCleaner::PHTrackCleaner(const std::string &name)
0027   : SubsysReco(name)
0028 {
0029 }
0030
0031 //____________________________________________________________________________..
0032 PHTrackCleaner::~PHTrackCleaner() = default;
0033
0034 //____________________________________________________________________________..
0035 int PHTrackCleaner::InitRun(PHCompositeNode *topNode)
0036 {
0037   int ret = GetNodes(topNode);
0038   if (ret != Fun4AllReturnCodes::EVENT_OK)
0039   {
0040     return ret;
0041   }
0042
0043   return ret;
0044 }
0045
0046 //____________________________________________________________________________..
0047 int PHTrackCleaner::process_event(PHCompositeNode * /*topNode*/)
0048 {
0049   if (Verbosity() > 0)
0050   {
0051     std::cout << PHWHERE << " track map size " << _track_map->size() << std::endl;
0052   }
0053
0054   std::set<unsigned int> track_keep_list;
0055   std::set<unsigned int> track_delete_list;
0056
0057   unsigned int good_track = 0;  // for diagnostic output only
0058   unsigned int ok_track = 0;    // tracks to keep
0059
0060   std::multimap<unsigned int, unsigned int> tpcid_track_mmap;
0061   std::set<unsigned int> tpc_id_set;
0062   // loop over the fitted tracks
0063   for (auto &it : *_track_map)
0064   {
0065     auto track_id = it.first;
0066     auto track = it.second;
0067     if (!track)
0068     {
0069       continue;
0070     }
0071
0072     auto tpc_seed = track->get_tpc_seed();
0073     unsigned int tpc_index = _tpc_seed_map->find(tpc_seed);
0074
0075     auto tpc_track_pair = std::make_pair(tpc_index, track_id);
0076
0077     tpc_id_set.insert(tpc_index);
0078     tpcid_track_mmap.insert(tpc_track_pair);
0079   }
0080
0081   if (Verbosity() > 0)
0082   {
0083     std::cout << " tpcid_track_mmap  size " << tpcid_track_mmap.size() << std::endl;
0084   }
0085
0086   // loop over the TPC seed ID's
0087
0088   for (unsigned int tpc_id : tpc_id_set)
0089   {
0090     if (Verbosity() > 1)
0091     {
0092       std::cout << " TPC ID " << tpc_id << std::endl;
0093     }
0094
0095     auto tpc_range = tpcid_track_mmap.equal_range(tpc_id);
0096
0097     unsigned int best_id = 99999;
0098     double min_chisq_df = 99999.0;
0099     unsigned int best_ndf = 1;
0100     for (auto it = tpc_range.first; it != tpc_range.second; ++it)
0101     {
0102       unsigned int track_id = it->second;
0103
0104       // note that the track no longer exists if it failed in the Acts fitter
0105       _track = _track_map->get(track_id);
0106
0107       if (_track)
0108       {
0109         unsigned int si_index = UINT_MAX;
0110         auto si_seed = _track->get_silicon_seed();
0111         if (si_seed)
0112       {
0113         si_index = _silicon_seed_map->find(si_seed);
0114       }
0115     else
0116       {
0117       if(Verbosity() > 1)
0118       {
0119         std::cout << "      no silicon seed found " << std::endl;
0120       }
0121     }
0122
0123         if (_pp_mode)
0124       {
0125         if (!si_seed)
0126           {
0127         // Tracks with no silicon seed in pp mode
0128         if (_track->get_chisq() / _track->get_ndf() < min_chisq_df && _track->get_ndf() > min_ndf && _track->get_ndf() != UINT_MAX)
0129           {
0130             best_id = track_id;
0131             best_ndf = _track->get_ndf();
0132             double qual = _track->get_chisq() / _track->get_ndf();
0133
0134             if (qual < quality_cut * 2)
0135               {
0136             // keep this TPC only track
0137             track_keep_list.insert(best_id);
0138             ok_track++;
0139             if (qual < quality_cut)
0140               {
0141                 good_track++;
0142               }
0143
0144             if (Verbosity() > 1)
0145               {
0146                 std::cout << "        keep unmatched track tpc_id " << tpc_id << " given best_id " << best_id
0147                       << " best_ndf " << best_ndf << " chisq/ndf " << qual << std::endl;
0148               }
0149               }
0150           }
0151
0152         continue;
0153           }
0154
0155       }
0156
0157         // Find the remaining silicon matched track with the best chisq/ndf
0158
0159         if (Verbosity() > 1)
0160       {
0161           std::cout << "        track ID " << track_id << " tpc index " << tpc_id << " si index " << si_index << " crossing " << _track->get_crossing()
0162                     << " chisq " << _track->get_chisq() << " ndf " << _track->get_ndf() << " min_chisq_df " << min_chisq_df << std::endl;
0163         }
0164
0165         // only accept tracks with ndf > min_ndf - very small ndf means something went wrong, as does ndf undefined
0166         if (_track->get_chisq() / _track->get_ndf() < min_chisq_df && _track->get_ndf() > min_ndf && _track->get_ndf() != UINT_MAX)
0167         {
0168           min_chisq_df = _track->get_chisq() / _track->get_ndf();
0169           best_id = track_id;
0170           best_ndf = _track->get_ndf();
0171         }
0172       }
0173     }
0174
0175     if (best_id != 99999)
0176     {
0177       double qual = min_chisq_df;
0178
0179       if (Verbosity() > 1)
0180       {
0181         std::cout << "        best track for tpc_id " << tpc_id << " has track_id " << best_id << " best_ndf " << best_ndf << " chisq/ndf " << qual << std::endl;
0182       }
0183
0184       if (qual < quality_cut * 2)
0185       {
0186         track_keep_list.insert(best_id);
0187         ok_track++;
0188         if (qual < quality_cut)
0189         {
0190           good_track++;
0191         }
0192       }
0193     }
0194     else
0195     {
0196       if (Verbosity() > 1)
0197       {
0198         std::cout << "        no track exists  for tpc_id " << tpc_id << std::endl;
0199       }
0200     }
0201   }
0202
0203   if (Verbosity() > 0)
0204   {
0205     std::cout << " Number of good tracks with qual < " << quality_cut << "  is " << good_track << " OK tracks " << ok_track << std::endl;
0206   }
0207
0208   // make a list of tracks that did not make the keep list
0209   for (auto &track_it : *_track_map)
0210   {
0211     auto id = track_it.first;
0212
0213     auto set_it = track_keep_list.find(id);
0214     if (set_it == track_keep_list.end())
0215     {
0216       if (Verbosity() > 1)
0217       {
0218         std::cout << "    add id " << id << " to track_delete_list " << std::endl;
0219       }
0220       track_delete_list.insert(id);
0221     }
0222   }
0223
0224   if (Verbosity() > 0)
0225   {
0226     std::cout << " track_delete_list size " << track_delete_list.size() << std::endl;
0227   }
0228
0229   // delete failed tracks
0230   for (unsigned int it : track_delete_list)
0231   {
0232     if (Verbosity() > 1)
0233     {
0234       std::cout << " erasing track ID " << it << std::endl;
0235     }
0236     _track_map->erase(it);
0237   }
0238
0239   if (Verbosity() > 0)
0240   {
0241     std::cout << "Track map size after choosing best silicon match: " << _track_map->size() << std::endl;
0242   }
0243
0244   return Fun4AllReturnCodes::EVENT_OK;
0245 }
0246
0247 int PHTrackCleaner::End(PHCompositeNode * /*topNode*/)
0248 {
0249   return Fun4AllReturnCodes::EVENT_OK;
0250 }
0251
0252 int PHTrackCleaner::GetNodes(PHCompositeNode *topNode)
0253 {
0254   _tpc_seed_map = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0255   if (!_tpc_seed_map)
0256   {
0257     std::cout << PHWHERE << " ERROR: Can't find TpcTrackSeedContainer: " << std::endl;
0258     return Fun4AllReturnCodes::ABORTEVENT;
0259   }
0260
0261   _silicon_seed_map = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0262   if (!_silicon_seed_map)
0263   {
0264     std::cout << PHWHERE << " ERROR: Can't find SiliconTrackSeedContainer " << std::endl;
0265     return Fun4AllReturnCodes::ABORTEVENT;
0266   }
0267
0268   _track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0269   if (!_track_map)
0270   {
0271     std::cout << PHWHERE << " ERROR: Can't find SvtxTrackMap: " << std::endl;
0272     return Fun4AllReturnCodes::ABORTEVENT;
0273   }
0274
0275   return Fun4AllReturnCodes::EVENT_OK;
0276 }