packages/trackreco/PHSiliconSeedPruner.cc

0001 #include "PHSiliconSeedPruner.h"
0002 #include "PHSiliconSeedPrunerHelper.h"
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005
0006 #include <phool/PHCompositeNode.h>
0007 #include <phool/PHIODataNode.h>
0008 #include <phool/PHNodeIterator.h>
0009 #include <phool/PHObject.h>
0010 #include <phool/PHRandomSeed.h>
0011 #include <phool/getClass.h>
0012 #include <phool/phool.h>
0013
0014 #include <trackbase_historic/TrackSeed.h>
0015 #include <trackbase_historic/TrackSeedContainer.h>
0016 #include <trackbase_historic/TrackSeedContainer_v1.h>
0017
0018 #include <climits>
0019 #include <iostream>
0020 #include <map>
0021 #include <set>
0022 #include <vector>
0023
0024 PHSiliconSeedPruner::PHSiliconSeedPruner(const std::string &name)
0025   : SubsysReco(name)
0026   , m_rng(gsl_rng_alloc(gsl_rng_mt19937))
0027 {
0028   set_random_seed(PHRandomSeed());  // default random seed
0029 }
0030
0031 PHSiliconSeedPruner::~PHSiliconSeedPruner()
0032 {
0033   gsl_rng_free(m_rng);
0034 }
0035
0036 void PHSiliconSeedPruner::set_random_seed(unsigned int seed)
0037 {
0038   m_randomSeed = seed;
0039   gsl_rng_set(m_rng, m_randomSeed);
0040 }
0041
0042 int PHSiliconSeedPruner::InitRun(PHCompositeNode *topNode)
0043 {
0044   if (m_trackMapName == m_outputTrackMapName)
0045   {
0046     std::cout << PHWHERE << " ERROR: Input and output TrackSeedContainer names are both " << m_trackMapName << std::endl;
0047     return Fun4AllReturnCodes::ABORTRUN;
0048   }
0049
0050   m_siliconSeeds = findNode::getClass<TrackSeedContainer>(topNode, m_trackMapName);
0051   if (!m_siliconSeeds)
0052   {
0053     std::cout << PHWHERE << " ERROR: Can't find " << m_trackMapName << std::endl;
0054     return Fun4AllReturnCodes::ABORTRUN;
0055   }
0056
0057   m_prunedSiliconSeeds = findNode::getClass<TrackSeedContainer>(topNode, m_outputTrackMapName);
0058   if (!m_prunedSiliconSeeds)
0059   {
0060     PHNodeIterator iter(topNode);
0061     PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0062     if (!dstNode)
0063     {
0064       std::cout << PHWHERE << " ERROR: Can't find DST node" << std::endl;
0065       return Fun4AllReturnCodes::ABORTRUN;
0066     }
0067
0068     PHNodeIterator dstIter(dstNode);
0069     PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0070     if (!svtxNode)
0071     {
0072       svtxNode = new PHCompositeNode("SVTX");
0073       dstNode->addNode(svtxNode);
0074     }
0075
0076     m_prunedSiliconSeeds = new TrackSeedContainer_v1;
0077     auto *outputNode = new PHIODataNode<PHObject>(m_prunedSiliconSeeds, m_outputTrackMapName, "PHObject");
0078     svtxNode->addNode(outputNode);
0079   }
0080
0081   return Fun4AllReturnCodes::EVENT_OK;
0082 }
0083
0084 int PHSiliconSeedPruner::process_event(PHCompositeNode * /*topNode*/)
0085 {
0086   m_prunedSiliconSeeds->Reset();
0087
0088   std::map<short int, std::vector<size_t>> seedIndicesByCrossing;
0089   std::set<size_t> selectedSeedIndices;
0090
0091   size_t inputSeedCount = 0;
0092   size_t invalidCrossingCount = 0;
0093   for (size_t seedIndex = 0; seedIndex < m_siliconSeeds->size(); ++seedIndex)
0094   {
0095     TrackSeed *seed = m_siliconSeeds->get(seedIndex);
0096     if (!seed)
0097     {
0098       continue;
0099     }
0100
0101     ++inputSeedCount;
0102     const short int crossing = seed->get_crossing();
0103     if (crossing == SHRT_MAX)
0104     {
0105       selectedSeedIndices.insert(seedIndex);
0106       ++invalidCrossingCount;
0107       continue;
0108     }
0109
0110     seedIndicesByCrossing[crossing].push_back(seedIndex);
0111   }
0112
0113   size_t uncertifiedSeedCount = 0;
0114   for (const auto &[crossing, seedIndices] : seedIndicesByCrossing)
0115   {
0116     const PHSiliconSeedPrunerHelper::Result result = PHSiliconSeedPrunerHelper::SelectSeeds(
0117         *m_siliconSeeds,
0118         seedIndices,
0119         m_rng,
0120         kMvtxLayerCount,
0121         kMaxRepresentatives,
0122         kSearchBudget,
0123         kHeuristicRestarts,
0124         kDebugMinimumGroupSize,
0125         Verbosity() > 0,
0126         Verbosity() > 1);
0127
0128     selectedSeedIndices.insert(result.selectedSeedIndices.begin(), result.selectedSeedIndices.end());
0129     uncertifiedSeedCount += result.uncertifiedSeedIndices.size();
0130
0131     if (Verbosity() > 0)
0132     {
0133       std::cout << Name() << ": crossing " << crossing
0134                 << " input seeds " << seedIndices.size()
0135                 << " selected seeds " << result.selectedSeedIndices.size()
0136                 << " uncertified seeds " << result.uncertifiedSeedIndices.size()
0137                 << std::endl;
0138     }
0139   }
0140
0141   for (const size_t seedIndex : selectedSeedIndices)
0142   {
0143     m_prunedSiliconSeeds->insert(m_siliconSeeds->get(seedIndex));
0144   }
0145
0146   if (Verbosity() > 0)
0147   {
0148     std::cout << Name()
0149               << ": input seeds " << inputSeedCount
0150               << ", retained seeds " << selectedSeedIndices.size()
0151               << ", retained invalid-crossing seeds " << invalidCrossingCount
0152               << ", uncertified selected seeds " << uncertifiedSeedCount
0153               << ", random seed " << m_randomSeed
0154               << std::endl;
0155   }
0156
0157   return Fun4AllReturnCodes::EVENT_OK;
0158 }