packages/mvtx/MvtxHitPruner.cc

0001 /**
0002  * @file mvtx/MvtxHitPruner.cc
0003  * @author D. McGlinchey
0004  * @date June 2018
0005  * @brief Implementation of MvtxHitPruner
0006  */
0007 #include "MvtxHitPruner.h"
0008 #include "CylinderGeom_Mvtx.h"
0009
0010 #include <g4detectors/PHG4CylinderGeom.h>
0011 #include <g4detectors/PHG4CylinderGeomContainer.h>
0012
0013 #include <trackbase/MvtxDefs.h>
0014 #include <trackbase/TrkrDefs.h>  // for hitkey, getLayer
0015 #include <trackbase/TrkrHitSet.h>
0016 #include <trackbase/TrkrHitSetContainer.h>
0017 #include <trackbase/TrkrHitv2.h>
0018
0019 #include <fun4all/Fun4AllReturnCodes.h>
0020 #include <fun4all/SubsysReco.h>  // for SubsysReco
0021
0022 #include <phool/PHCompositeNode.h>
0023 #include <phool/PHIODataNode.h>
0024 #include <phool/PHNode.h>  // for PHNode
0025 #include <phool/PHNodeIterator.h>
0026 #include <phool/PHObject.h>  // for PHObject
0027 #include <phool/getClass.h>
0028 #include <phool/phool.h>  // for PHWHERE
0029
0030 #include <TMatrixFfwd.h>    // for TMatrixF
0031 #include <TMatrixT.h>       // for TMatrixT, operator*
0032 #include <TMatrixTUtils.h>  // for TMatrixTRow
0033 #include <TVector3.h>
0034
0035 #pragma GCC diagnostic push
0036 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0037 #include <boost/graph/adjacency_list.hpp>
0038 #pragma GCC diagnostic pop
0039
0040 #include <boost/graph/connected_components.hpp>
0041
0042 #include <array>
0043 #include <cmath>
0044 #include <cstdlib>  // for exit
0045 #include <iostream>
0046 #include <map>  // for multimap<>::iterator
0047 #include <set>  // for set, set<>::iterator
0048 #include <string>
0049 #include <vector>  // for vector
0050
0051 namespace
0052 {
0053   //! range adaptor to be able to use range-based for loop
0054   template<class T> class range_adaptor
0055   {
0056     public:
0057     range_adaptor( const T& range ):m_range(range){}
0058     const typename T::first_type& begin() {return m_range.first;}
0059     const typename T::second_type& end() {return m_range.second;}
0060     private:
0061     T m_range;
0062   };
0063 }
0064
0065 MvtxHitPruner::MvtxHitPruner(const std::string &name)
0066   : SubsysReco(name)
0067 {
0068 }
0069
0070 int MvtxHitPruner::InitRun(PHCompositeNode * /*topNode*/)
0071 { return Fun4AllReturnCodes::EVENT_OK; }
0072
0073 int MvtxHitPruner::process_event(PHCompositeNode *topNode)
0074 {
0075   // get node containing the digitized hits
0076   m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0077   if (!m_hits)
0078   {
0079     std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
0080     return Fun4AllReturnCodes::ABORTRUN;
0081   }
0082
0083   /*
0084    * We want to combine all strobe values for a given hitset
0085    * Start by looping over all MVTX hitsets and making a map of physical sensor
0086    * to hitsetkey-with-strobe
0087    */
0088
0089   // map (bare hitset, hitset with strobe)
0090   std::multimap<TrkrDefs::hitsetkey, TrkrDefs::hitsetkey> hitset_multimap;
0091
0092   // list of all physical sensor hitsetkeys (with strobe set to zero)
0093   std::set<TrkrDefs::hitsetkey> bare_hitset_set;
0094
0095   const auto hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::mvtxId);
0096   for( const auto& [hitsetkey,hitset]:range_adaptor(hitsetrange) )
0097   {
0098
0099     // get strobe, skip if already zero
0100     const int strobe = MvtxDefs::getStrobeId(hitsetkey);
0101     if( strobe == 0 ) continue;
0102
0103     // get the hitsetkey value for strobe 0
0104     const auto bare_hitsetkey = MvtxDefs::resetStrobe(hitsetkey);
0105     hitset_multimap.emplace(bare_hitsetkey, hitsetkey);
0106     bare_hitset_set.insert(bare_hitsetkey);
0107
0108     if (Verbosity())
0109     {
0110       std::cout << " found hitsetkey " << hitsetkey << " for bare_hitsetkey " << bare_hitsetkey << std::endl;
0111     }
0112   }
0113
0114   // Now consolidate all hits into the hitset with strobe 0, and delete the
0115   // other hitsets
0116   //==============================================================
0117   for (const auto& bare_hitsetkey : bare_hitset_set)
0118   {
0119     // find matching hitset of creater
0120     auto bare_hitset = (m_hits->findOrAddHitSet(bare_hitsetkey))->second;
0121     if (Verbosity())
0122     {
0123       std::cout
0124         << "MvtxHitPruner::process_event - bare_hitset " << bare_hitsetkey
0125         << " initially has " << bare_hitset->size() << " hits "
0126         << std::endl;
0127     }
0128
0129     // get all hitsets with non-zero strobe that match the bare hitset key
0130     auto bare_hitsetrange = hitset_multimap.equal_range(bare_hitsetkey);
0131     for( const auto& [unused,hitsetkey]:range_adaptor(bare_hitsetrange) )
0132     {
0133       const int strobe = MvtxDefs::getStrobeId(hitsetkey);
0134       if( strobe == 0 ) continue;
0135
0136       if (Verbosity())
0137       {
0138         std::cout << "MvtxHitPruner::process_event -"
0139           << " process hitsetkey " << hitsetkey
0140           << " from strobe " << strobe
0141           << " for bare_hitsetkey " << bare_hitsetkey
0142           << std::endl;
0143       }
0144
0145       // copy all hits to the hitset with strobe 0
0146       auto hitset = m_hits->findHitSet(hitsetkey);
0147
0148       if (Verbosity())
0149       {
0150         std::cout << "MvtxHitPruner::process_event - hitsetkey " << hitsetkey
0151           << " has strobe " << strobe << " and has " << hitset->size()
0152           << " hits,  so copy it" << std::endl;
0153       }
0154
0155       TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0156       for( const auto& [hitkey,old_hit]:range_adaptor(hitrangei) )
0157       {
0158         if (Verbosity())
0159         {
0160           std::cout << "MvtxHitPruner::process_event - found hitkey " << hitkey << std::endl;
0161         }
0162
0163         // if it is already there, leave it alone, this is a duplicate hit
0164         if (bare_hitset->getHit(hitkey))
0165         {
0166           if (Verbosity())
0167           {
0168             std::cout
0169               << "MvtxHitPruner::process_event - hitkey " << hitkey
0170               << " is already in bare hitsest, do not copy"
0171               << std::endl;
0172           }
0173           continue;
0174         }
0175
0176         // otherwise copy the hit over
0177         if (Verbosity())
0178         {
0179           std::cout
0180             << "MvtxHitPruner::process_event - copying over hitkey "
0181             << hitkey << std::endl;
0182         }
0183
0184         auto new_hit = new TrkrHitv2;
0185         new_hit->CopyFrom(old_hit);
0186         bare_hitset->addHitSpecificKey(hitkey, new_hit);
0187       }
0188
0189       // all hits are copied over to the strobe zero hitset, remove this
0190       // hitset
0191       m_hits->removeHitSet(hitsetkey);
0192     }
0193   }
0194
0195   return Fun4AllReturnCodes::EVENT_OK;
0196 }