packages/trackreco/PHSiliconSeedMerger.cc

0001
0002 #include "PHSiliconSeedMerger.h"
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <phool/PHCompositeNode.h>
0006 #include <phool/PHObject.h>
0007 #include <phool/PHTimer.h>
0008 #include <phool/getClass.h>
0009 #include <phool/phool.h>
0010
0011 #include <trackbase/MvtxDefs.h>
0012 #include <trackbase/TrkrDefs.h>
0013
0014 #include <trackbase_historic/TrackSeed.h>
0015 #include <trackbase_historic/TrackSeedContainer.h>
0016
0017 #include <phool/PHCompositeNode.h>
0018
0019 //____________________________________________________________________________..
0020 PHSiliconSeedMerger::PHSiliconSeedMerger(const std::string& name)
0021   : SubsysReco(name)
0022 {
0023 }
0024
0025 //____________________________________________________________________________..
0026 PHSiliconSeedMerger::~PHSiliconSeedMerger()
0027 {
0028 }
0029
0030 //____________________________________________________________________________..
0031 int PHSiliconSeedMerger::Init(PHCompositeNode*)
0032 {
0033   return Fun4AllReturnCodes::EVENT_OK;
0034 }
0035
0036 //____________________________________________________________________________..
0037 int PHSiliconSeedMerger::InitRun(PHCompositeNode* topNode)
0038 {
0039   int ret = getNodes(topNode);
0040   return ret;
0041 }
0042
0043 //____________________________________________________________________________..
0044 int PHSiliconSeedMerger::process_event(PHCompositeNode*)
0045 {
0046   std::multimap<unsigned int, std::set<TrkrDefs::cluskey>> matches;
0047   std::set<unsigned int> seedsToDelete;
0048
0049   if (Verbosity() > 0)
0050   {
0051     std::cout << "Silicon seed track container has " << m_siliconTracks->size() << std::endl;
0052   }
0053
0054   for (unsigned int track1ID = 0;
0055        track1ID != m_siliconTracks->size();
0056        ++track1ID)
0057   {
0058     TrackSeed* track1 = m_siliconTracks->get(track1ID);
0059
0060     if (seedsToDelete.find(track1ID) != seedsToDelete.end())
0061     {
0062       continue;
0063     }
0064
0065     std::set<TrkrDefs::cluskey> mvtx1Keys;
0066     int track1Strobe = std::numeric_limits<int>::quiet_NaN();
0067     for (auto iter = track1->begin_cluster_keys();
0068          iter != track1->end_cluster_keys();
0069          ++iter)
0070     {
0071       TrkrDefs::cluskey ckey = *iter;
0072       auto trkrid = TrkrDefs::getTrkrId(ckey);
0073       if (m_mvtxOnly && trkrid == TrkrDefs::TrkrId::inttId)
0074       {
0075         continue;
0076       }
0077       track1Strobe = MvtxDefs::getStrobeId(ckey);
0078       mvtx1Keys.insert(ckey);
0079     }
0080
0081     /// We can speed up the code by only iterating over the track seeds
0082     /// that are further in the map container from the current track,
0083     /// since the comparison of e.g. track 1 with track 2 doesn't need
0084     /// to be repeated with track 2 to track 1.
0085     for (unsigned int track2ID = track1ID;
0086          track2ID != m_siliconTracks->size();
0087          ++track2ID)
0088     {
0089       if (track1ID == track2ID)
0090       {
0091         continue;
0092       }
0093
0094       TrackSeed* track2 = m_siliconTracks->get(track2ID);
0095       if (track2 == nullptr)
0096       {
0097         continue;
0098       }
0099       int track2Strobe = std::numeric_limits<int>::quiet_NaN();
0100
0101       std::set<TrkrDefs::cluskey> mvtx2Keys;
0102       for (TrackSeed::ConstClusterKeyIter iter = track2->begin_cluster_keys();
0103            iter != track2->end_cluster_keys();
0104            ++iter)
0105       {
0106         TrkrDefs::cluskey ckey = *iter;
0107         auto trkrid = TrkrDefs::getTrkrId(ckey);
0108         if (m_mvtxOnly && trkrid == TrkrDefs::TrkrId::inttId)
0109         {
0110           continue;
0111         }
0112         mvtx2Keys.insert(ckey);
0113         track2Strobe = MvtxDefs::getStrobeId(ckey);
0114       }
0115
0116       std::vector<TrkrDefs::cluskey> intersection;
0117       std::set_intersection(mvtx1Keys.begin(),
0118                             mvtx1Keys.end(),
0119                             mvtx2Keys.begin(),
0120                             mvtx2Keys.end(),
0121                             std::back_inserter(intersection));
0122
0123       /// If we have two clusters in common in the triplet, it is likely
0124       /// from the same track
0125       if (intersection.size() > m_clusterOverlap && track1Strobe == track2Strobe)
0126       {
0127         if (Verbosity() > 2)
0128         {
0129           std::cout << "Track " << track1ID << " keys " << std::endl;
0130           for (auto& key : mvtx1Keys)
0131           {
0132             std::cout << "   ckey: " << key << std::endl;
0133           }
0134           std::cout << "Track " << track2ID << " keys " << std::endl;
0135           for (auto& key : mvtx2Keys)
0136           {
0137             std::cout << "   ckey: " << key << std::endl;
0138           }
0139           std::cout << "Intersection keys " << std::endl;
0140           for (auto& key : intersection)
0141           {
0142             std::cout << "   ckey: " << key << std::endl;
0143           }
0144         }
0145
0146         for (auto& key : mvtx2Keys)
0147         {
0148           mvtx1Keys.insert(key);
0149         }
0150
0151         if (Verbosity() > 2)
0152         {
0153           std::cout << "Match IDed" << std::endl;
0154           for (auto& key : mvtx1Keys)
0155           {
0156             std::cout << "  total track keys " << key << std::endl;
0157           }
0158         }
0159
0160         matches.insert(std::make_pair(track1ID, mvtx1Keys));
0161         seedsToDelete.insert(track2ID);
0162         break;
0163       }
0164     }
0165   }
0166
0167   for (const auto& [trackKey, mvtxKeys] : matches)
0168   {
0169     auto track = m_siliconTracks->get(trackKey);
0170     if (Verbosity() > 2)
0171     {
0172       std::cout << "original track: " << std::endl;
0173       track->identify();
0174     }
0175
0176     for (auto& key : mvtxKeys)
0177     {
0178       if (track->find_cluster_key(key) == track->end_cluster_keys())
0179       {
0180         track->insert_cluster_key(key);
0181         if (Verbosity() > 2)
0182           std::cout << "adding " << key << std::endl;
0183       }
0184     }
0185   }
0186
0187   for (const auto& key : seedsToDelete)
0188   {
0189     if (Verbosity() > 2)
0190     {
0191       std::cout << "Erasing track " << key << std::endl;
0192     }
0193     m_siliconTracks->erase(key);
0194   }
0195
0196   if (Verbosity() > 2)
0197   {
0198     for (const auto& seed : *m_siliconTracks)
0199     {
0200       if (!seed) continue;
0201       seed->identify();
0202     }
0203   }
0204
0205   return Fun4AllReturnCodes::EVENT_OK;
0206 }
0207
0208 //____________________________________________________________________________..
0209 int PHSiliconSeedMerger::ResetEvent(PHCompositeNode*)
0210 {
0211   return Fun4AllReturnCodes::EVENT_OK;
0212 }
0213
0214 //____________________________________________________________________________..
0215 int PHSiliconSeedMerger::End(PHCompositeNode*)
0216 {
0217   return Fun4AllReturnCodes::EVENT_OK;
0218 }
0219
0220 int PHSiliconSeedMerger::getNodes(PHCompositeNode* topNode)
0221 {
0222   m_siliconTracks = findNode::getClass<TrackSeedContainer>(topNode, m_trackMapName.c_str());
0223   if (!m_siliconTracks)
0224   {
0225     std::cout << PHWHERE << "No silicon track container, can't merge seeds"
0226               << std::endl;
0227     return Fun4AllReturnCodes::ABORTEVENT;
0228   }
0229
0230   return Fun4AllReturnCodes::EVENT_OK;
0231 }