packages/mvtx/MvtxClusterPruner.cc

0001 /**
0002  * @file mvtx/MvtxClusterPruner.cc
0003  * @author Hugo Pereira Da Costa <hugo.pereira-da-costa@lanl.gov>
0004  * @date May 2025
0005  * @brief Implementation of MvtxClusterPruner
0006  */
0007
0008 #include "MvtxClusterPruner.h"
0009
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011 #include <phool/getClass.h>
0012
0013 #include <trackbase/MvtxDefs.h>
0014 #include <trackbase/TrkrCluster.h>
0015 #include <trackbase/TrkrClusterHitAssoc.h>
0016 #include <trackbase/TrkrClusterContainer.h>
0017
0018 #include <set>
0019 #include <algorithm>
0020
0021 namespace
0022 {
0023   //! range adaptor to be able to use range-based for loop
0024   template<class T> class range_adaptor
0025   {
0026     public:
0027     range_adaptor( const T& range ):m_range(range){}
0028     const typename T::first_type& begin() {return m_range.first;}
0029     const typename T::second_type& end() {return m_range.second;}
0030     private:
0031     T m_range;
0032   };
0033
0034   // print cluster information
0035   void print_cluster_information( TrkrDefs::cluskey ckey, TrkrCluster* cluster )
0036   {
0037     if( cluster )
0038     {
0039       std::cout << " MVTX cluster: " << ckey
0040         << " position: (" << cluster->getLocalX() << ", " << cluster->getLocalY() << ")"
0041         << " size: " << (int)cluster->getSize()
0042         << " layer: " << (int)TrkrDefs::getLayer(ckey)
0043         << " stave: " << (int) MvtxDefs::getStaveId(ckey)
0044         << " chip: " << (int)MvtxDefs::getChipId(ckey)
0045         << " strobe: " << (int)MvtxDefs::getStrobeId(ckey)
0046         << " index: " << (int)TrkrDefs::getClusIndex(ckey)
0047         << std::endl;
0048     } else {
0049       std::cout << " MVTX cluster: " << ckey
0050         << " layer: " << (int)TrkrDefs::getLayer(ckey)
0051         << " stave: " << (int) MvtxDefs::getStaveId(ckey)
0052         << " chip: " << (int)MvtxDefs::getChipId(ckey)
0053         << " strobe: " << (int)MvtxDefs::getStrobeId(ckey)
0054         << " index: " << (int)TrkrDefs::getClusIndex(ckey)
0055         << std::endl;
0056     }
0057   }
0058
0059   using hitkeyset_t = std::set<TrkrDefs::hitkey>;
0060   using clustermap_t = std::map<TrkrDefs::cluskey,hitkeyset_t>;
0061
0062 }
0063
0064 //_____________________________________________________________________________
0065 MvtxClusterPruner::MvtxClusterPruner(const std::string &name)
0066   : SubsysReco(name)
0067 {
0068 }
0069
0070 //_____________________________________________________________________________
0071 int MvtxClusterPruner::InitRun(PHCompositeNode * /*topNode*/)
0072 {
0073   std::cout << "MvtxClusterPruner::InitRun - m_use_strict_matching: " << m_use_strict_matching << std::endl;
0074   return Fun4AllReturnCodes::EVENT_OK;
0075 }
0076
0077 //_____________________________________________________________________________
0078 int MvtxClusterPruner::process_event(PHCompositeNode *topNode)
0079 {
0080   // load relevant nodes
0081   auto trkrclusters = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0082   if( !trkrclusters )
0083   {
0084     std::cout << "MvtxClusterPruner::process_event - TRKR_CLUSTER not found. Doing nothing" << std::endl;
0085     return Fun4AllReturnCodes::EVENT_OK;
0086   }
0087
0088   auto clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0089   if( !clusterhitassoc )
0090   {
0091     std::cout << "MvtxClusterPruner::process_event - TRKR_CLUSTERHITASSOC not found. Doing nothing" << std::endl;
0092     return Fun4AllReturnCodes::EVENT_OK;
0093   }
0094
0095   // lambda method to create map of cluster keys and associated hits
0096   auto get_cluster_map = [trkrclusters,clusterhitassoc]( TrkrDefs::hitsetkey key )
0097   {
0098     clustermap_t out;
0099
0100     // get all clusters for this hitsetkey
0101     const auto cluster_range= trkrclusters->getClusters(key);
0102     for( const auto& [ckey,cluster]:range_adaptor(cluster_range) )
0103     {
0104       // get associated hits
0105       const auto& hit_range = clusterhitassoc->getHits(ckey);
0106       hitkeyset_t hitkeys;
0107       std::transform(hit_range.first, hit_range.second, std::inserter(hitkeys,hitkeys.end()),
0108         [](const TrkrClusterHitAssoc::Map::value_type& pair ){ return pair.second; });
0109       out.emplace(ckey,std::move(hitkeys));
0110     }
0111     return out;
0112   };
0113
0114   // loop over MVTX hitset keys
0115   const auto hitsetkeys = trkrclusters->getHitSetKeys(TrkrDefs::mvtxId);
0116   for( const auto& hitsetkey:hitsetkeys )
0117   {
0118     // get layer, stave, chip and current strobe
0119     const auto layer = TrkrDefs::getLayer(hitsetkey);
0120     const auto stave = MvtxDefs::getStaveId(hitsetkey);
0121     const auto chip = MvtxDefs::getChipId(hitsetkey);
0122     const auto current_strobe = MvtxDefs::getStrobeId(hitsetkey);
0123
0124     // get clusters for this hitsetkey
0125     const auto cluster_map1 = get_cluster_map(hitsetkey);
0126
0127     // get clusters for the next strobe
0128     int next_strobe = current_strobe+1;
0129     const auto hitsetkey_next_strobe = MvtxDefs::genHitSetKey(layer, stave, chip, next_strobe);
0130     const auto clusterk_map2 = get_cluster_map(hitsetkey_next_strobe);
0131
0132     // loop over clusters from first range
0133     for( auto [ckey1,hitkeys1]:cluster_map1)
0134     {
0135       // increment counter
0136       ++m_cluster_counter_total;
0137
0138       // get correcponding cluser
0139       auto cluster1 = Verbosity() ? trkrclusters->findCluster(ckey1):nullptr;
0140
0141       // loop over clusters from second range
0142       for( auto [ckey2,hitkeys2]:clusterk_map2)
0143       {
0144         auto cluster2 = Verbosity() ? trkrclusters->findCluster(ckey2):nullptr;
0145
0146         if( m_use_strict_matching )
0147         {
0148           // see if hitsets are identical
0149           if(hitkeys1 == hitkeys2)
0150           {
0151             // increment counter
0152             ++m_cluster_counter_deleted;
0153
0154             if( Verbosity() )
0155             {
0156               std::cout << "Removing cluster ";
0157               print_cluster_information( ckey2, cluster2);
0158
0159               std::cout << "Keeping  cluster ";
0160               print_cluster_information( ckey1, cluster1);
0161             }
0162
0163             // always remove second cluster
0164             trkrclusters->removeCluster(ckey2);
0165             break;
0166           }
0167
0168         } else {
0169
0170           // make sure first set is larger than second
0171           const bool swapped = hitkeys2.size() > hitkeys1.size();
0172           if( swapped ) { std::swap(hitkeys2,hitkeys1); }
0173
0174           // see if hitkeys2 is a subset of hitkeys1
0175           if( std::includes(hitkeys1.begin(), hitkeys1.end(), hitkeys2.begin(), hitkeys2.end()) )
0176           {
0177             // increment counter
0178             ++m_cluster_counter_deleted;
0179
0180             if( swapped )
0181             {
0182
0183               if( Verbosity() )
0184               {
0185                 std::cout << "Removing cluster ";
0186                 print_cluster_information( ckey1, cluster1);
0187
0188                 std::cout << "Keeping  cluster ";
0189                 print_cluster_information( ckey2, cluster2);
0190               }
0191
0192               // remove first cluster
0193               trkrclusters->removeCluster(ckey1);
0194               break;
0195             } else {
0196
0197               if( Verbosity() )
0198               {
0199                 std::cout << "Removing cluster ";
0200                 print_cluster_information( ckey2, cluster2);
0201
0202                 std::cout << "Keeping  cluster ";
0203                 print_cluster_information( ckey1, cluster1);
0204               }
0205
0206               // remove second cluster
0207               trkrclusters->removeCluster(ckey2);
0208             }
0209           }
0210
0211         } // strict matching
0212
0213       } // second cluster loop
0214     } // first cluster loop
0215   } // hitsetkey loop
0216
0217   return Fun4AllReturnCodes::EVENT_OK;
0218
0219 }
0220
0221 //_____________________________________________________________________________
0222 int MvtxClusterPruner::End(PHCompositeNode * /*topNode*/)
0223 {
0224     std::cout << "MvtxClusterPruner::End -"
0225     << " m_cluster_counter_total: " << m_cluster_counter_total
0226     << std::endl;
0227
0228   std::cout << "MvtxClusterPruner::End -"
0229     << " m_cluster_counter_deleted: " << m_cluster_counter_deleted
0230     << " fraction: " << double( m_cluster_counter_deleted )/m_cluster_counter_total
0231     << std::endl;
0232
0233
0234   return Fun4AllReturnCodes::EVENT_OK;
0235 }