Number of ADC Bins histogram QA

offline/QA/Tpc/TpcRawHitQA.cc

@@ -99,6 +99,8 @@ int TpcRawHitQA::InitRun(PHCompositeNode *topNode)
 	}
     }
 
+  h_nadc_bins_event = dynamic_cast<TH1 *>(hm->getHisto(std::string(getHistoPrefix() + "nadc_bins_event")));
+
   h_xy_N = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "xyPos_North")));
   h_xy_S = dynamic_cast<TH2 *>(hm->getHisto(std::string(getHistoPrefix() + "xyPos_South")));
 
@@ -120,6 +122,8 @@ int TpcRawHitQA::process_event(PHCompositeNode * /*unused*/)
   float nhit_sectors_fees_sampas[24][26][8] = {{{0}}};
 
   unsigned int raw_hit_num = 0;
+  unsigned int nadc_bins_event = 0;
+
   for (TpcRawHitContainer *&rawhitcont : rawhitcont_vec)
     {
       raw_hit_num = rawhitcont->get_nhits();
@@ -215,6 +219,9 @@ int TpcRawHitQA::process_event(PHCompositeNode * /*unused*/)
 	    {
 	      const uint16_t sampleN = adc_iterator->CurrentTimeBin();
 	      const uint16_t adc = adc_iterator->CurrentAdc();
+
+        ++nadc_bins_event;
+
 	      if (adc - median <= (std::max(5 * stdDev, (float) 20.)))
 		{
 		  continue;
@@ -249,6 +256,8 @@ int TpcRawHitQA::process_event(PHCompositeNode * /*unused*/)
       return Fun4AllReturnCodes::EVENT_OK;
     }
 
+  h_nadc_bins_event->Fill(nadc_bins_event);
+
   for (int s = 0; s < 24; s++)
     {
       h_nhits_sectors[s]->Fill(nhit_sectors[s]);
@@ -309,6 +318,8 @@ void TpcRawHitQA::createHistos()
 	}
     }
 
+  hm->registerHisto(new TH1F(std::string(getHistoPrefix() + "nadc_bins_event").c_str(), "Number of ADC Bins per Event;Number of ADC Bins/Event;Entries", 2500, 0, 1000000));
+
   hm->registerHisto(new TH2F(std::string(getHistoPrefix() + "xyPos_North").c_str(), "Hit XY distribution (North);X [mm];Y [mm]", 400, -800, 800, 400, -800, 800));
   hm->registerHisto(new TH2F(std::string(getHistoPrefix() + "xyPos_South").c_str(), "Hit XY distribution (South);X [mm];Y [mm]", 400, -800, 800, 400, -800, 800));
 }

offline/QA/Tpc/TpcRawHitQA.h

@@ -45,6 +45,8 @@ class TpcRawHitQA : public SubsysReco
   TH2* h_nhits_sectors_fees_sampas[24][26]{{nullptr}};
   TH1* h_nhits_sam[24][3]{{nullptr}};
   TH1* h_adc[24][3]{{nullptr}};
+  TH1* h_nadc_bins_event{nullptr};
+
   TH2* h_xy_N{nullptr};
   TH2* h_xy_S{nullptr};
 

offline/packages/trackreco/Makefile.am

@@ -64,6 +64,7 @@ pkginclude_HEADERS = \
   PHRaveVertexing.h \
   PHSiliconHelicalPropagator.h \
   PHSiliconSeedMerger.h \
+  PHSiliconSeedPruner.h \
   PHSimpleKFProp.h \
   PHStreakFinder.h \
   PHTpcClusterMover.h \
@@ -87,6 +88,9 @@ pkginclude_HEADERS = \
   TpcSeedFilter.h \
   WeightedFitter.h
 
+noinst_HEADERS = \
+  PHSiliconSeedPrunerHelper.h
+
 ROOTDICTS = \
   AssocInfoContainer_Dict.cc \
   AssocInfoContainerv1_Dict.cc
@@ -145,6 +149,8 @@ libtrack_reco_la_SOURCES = \
   PHRaveVertexing.cc \
   PHSiliconHelicalPropagator.cc \
   PHSiliconSeedMerger.cc \
+  PHSiliconSeedPruner.cc \
+  PHSiliconSeedPrunerHelper.cc \
   PHSiliconTpcTrackMatching.cc \
   PHSiliconCosmicSeeding.cc \
   PHSimpleKFProp.cc \

offline/packages/trackreco/PHSiliconSeedPruner.cc

@@ -0,0 +1,158 @@
+#include "PHSiliconSeedPruner.h"
+#include "PHSiliconSeedPrunerHelper.h"
+
+#include <fun4all/Fun4AllReturnCodes.h>
+
+#include <phool/PHCompositeNode.h>
+#include <phool/PHIODataNode.h>
+#include <phool/PHNodeIterator.h>
+#include <phool/PHObject.h>
+#include <phool/PHRandomSeed.h>
+#include <phool/getClass.h>
+#include <phool/phool.h>
+
+#include <trackbase_historic/TrackSeed.h>
+#include <trackbase_historic/TrackSeedContainer.h>
+#include <trackbase_historic/TrackSeedContainer_v1.h>
+
+#include <climits>
+#include <iostream>
+#include <map>
+#include <set>
+#include <vector>
+
+PHSiliconSeedPruner::PHSiliconSeedPruner(const std::string &name)
+  : SubsysReco(name)
+  , m_rng(gsl_rng_alloc(gsl_rng_mt19937))
+{
+  set_random_seed(PHRandomSeed());  // default random seed
+}
+
+PHSiliconSeedPruner::~PHSiliconSeedPruner()
+{
+  gsl_rng_free(m_rng);
+}
+
+void PHSiliconSeedPruner::set_random_seed(unsigned int seed)
+{
+  m_randomSeed = seed;
+  gsl_rng_set(m_rng, m_randomSeed);
+}
+
+int PHSiliconSeedPruner::InitRun(PHCompositeNode *topNode)
+{
+  if (m_trackMapName == m_outputTrackMapName)
+  {
+    std::cout << PHWHERE << " ERROR: Input and output TrackSeedContainer names are both " << m_trackMapName << std::endl;
+    return Fun4AllReturnCodes::ABORTRUN;
+  }
+
+  m_siliconSeeds = findNode::getClass<TrackSeedContainer>(topNode, m_trackMapName);
+  if (!m_siliconSeeds)
+  {
+    std::cout << PHWHERE << " ERROR: Can't find " << m_trackMapName << std::endl;
+    return Fun4AllReturnCodes::ABORTRUN;
+  }
+
+  m_prunedSiliconSeeds = findNode::getClass<TrackSeedContainer>(topNode, m_outputTrackMapName);
+  if (!m_prunedSiliconSeeds)
+  {
+    PHNodeIterator iter(topNode);
+    PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
+    if (!dstNode)
+    {
+      std::cout << PHWHERE << " ERROR: Can't find DST node" << std::endl;
+      return Fun4AllReturnCodes::ABORTRUN;
+    }
+
+    PHNodeIterator dstIter(dstNode);
+    PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(dstIter.findFirst("PHCompositeNode", "SVTX"));
+    if (!svtxNode)
+    {
+      svtxNode = new PHCompositeNode("SVTX");
+      dstNode->addNode(svtxNode);
+    }
+
+    m_prunedSiliconSeeds = new TrackSeedContainer_v1;
+    auto *outputNode = new PHIODataNode<PHObject>(m_prunedSiliconSeeds, m_outputTrackMapName, "PHObject");
+    svtxNode->addNode(outputNode);
+  }
+
+  return Fun4AllReturnCodes::EVENT_OK;
+}
+
+int PHSiliconSeedPruner::process_event(PHCompositeNode * /*topNode*/)
+{
+  m_prunedSiliconSeeds->Reset();
+
+  std::map<short int, std::vector<size_t>> seedIndicesByCrossing;
+  std::set<size_t> selectedSeedIndices;
+
+  size_t inputSeedCount = 0;
+  size_t invalidCrossingCount = 0;
+  for (size_t seedIndex = 0; seedIndex < m_siliconSeeds->size(); ++seedIndex)
+  {
+    TrackSeed *seed = m_siliconSeeds->get(seedIndex);
+    if (!seed)
+    {
+      continue;
+    }
+
+    ++inputSeedCount;
+    const short int crossing = seed->get_crossing();
+    if (crossing == SHRT_MAX)
+    {
+      selectedSeedIndices.insert(seedIndex);
+      ++invalidCrossingCount;
+      continue;
+    }
+
+    seedIndicesByCrossing[crossing].push_back(seedIndex);
+  }
+
+  size_t uncertifiedSeedCount = 0;
+  for (const auto &[crossing, seedIndices] : seedIndicesByCrossing)
+  {
+    const PHSiliconSeedPrunerHelper::Result result = PHSiliconSeedPrunerHelper::SelectSeeds(
+        *m_siliconSeeds,
+        seedIndices,
+        m_rng,
+        kMvtxLayerCount,
+        kMaxRepresentatives,
+        kSearchBudget,
+        kHeuristicRestarts,
+        kDebugMinimumGroupSize,
+        Verbosity() > 0,
+        Verbosity() > 1);
+
+    selectedSeedIndices.insert(result.selectedSeedIndices.begin(), result.selectedSeedIndices.end());
+    uncertifiedSeedCount += result.uncertifiedSeedIndices.size();
+
+    if (Verbosity() > 0)
+    {
+      std::cout << Name() << ": crossing " << crossing
+                << " input seeds " << seedIndices.size()
+                << " selected seeds " << result.selectedSeedIndices.size()
+                << " uncertified seeds " << result.uncertifiedSeedIndices.size()
+                << std::endl;
+    }
+  }
+
+  for (const size_t seedIndex : selectedSeedIndices)
+  {
+    m_prunedSiliconSeeds->insert(m_siliconSeeds->get(seedIndex));
+  }
+
+  if (Verbosity() > 0)
+  {
+    std::cout << Name()
+              << ": input seeds " << inputSeedCount
+              << ", retained seeds " << selectedSeedIndices.size()
+              << ", retained invalid-crossing seeds " << invalidCrossingCount
+              << ", uncertified selected seeds " << uncertifiedSeedCount
+              << ", random seed " << m_randomSeed
+              << std::endl;
+  }
+
+  return Fun4AllReturnCodes::EVENT_OK;
+}

offline/packages/trackreco/PHSiliconSeedPruner.h

@@ -0,0 +1,51 @@
+// Tell emacs that this is a C++ source
+//  -*- C++ -*-.
+
+#ifndef PHSILICONSEEDPRUNER_H
+#define PHSILICONSEEDPRUNER_H
+
+#include <fun4all/SubsysReco.h>
+
+#include <gsl/gsl_rng.h>
+
+#include <cstddef>
+#include <string>
+
+class PHCompositeNode;
+class TrackSeedContainer;
+
+class PHSiliconSeedPruner : public SubsysReco
+{
+ public:
+  PHSiliconSeedPruner(const std::string& name = "PHSiliconSeedPruner");
+  ~PHSiliconSeedPruner() override;
+
+  PHSiliconSeedPruner(const PHSiliconSeedPruner&) = delete;
+  PHSiliconSeedPruner& operator=(const PHSiliconSeedPruner&) = delete;
+  PHSiliconSeedPruner(PHSiliconSeedPruner&&) = delete;
+  PHSiliconSeedPruner& operator=(PHSiliconSeedPruner&&) = delete;
+
+  int InitRun(PHCompositeNode* topNode) override;
+  int process_event(PHCompositeNode* topNode) override;
+
+  void set_track_map_name(const std::string& name) { m_trackMapName = name; }
+  void set_output_track_map_name(const std::string& name) { m_outputTrackMapName = name; }
+  void set_random_seed(unsigned int seed);
+
+ private:
+  static constexpr std::size_t kMvtxLayerCount = 3;
+  static constexpr std::size_t kMaxRepresentatives = 10;
+  static constexpr std::size_t kSearchBudget = 2000000;
+  static constexpr std::size_t kHeuristicRestarts = 256;
+  static constexpr std::size_t kDebugMinimumGroupSize = 10;
+
+  TrackSeedContainer* m_siliconSeeds{nullptr};
+  TrackSeedContainer* m_prunedSiliconSeeds{nullptr};
+  std::string m_trackMapName{"SiliconTrackSeedContainer"};              // the default input node name for the TrackSeedContainer
+  std::string m_outputTrackMapName{"PrunedSiliconTrackSeedContainer"};  // the default output node name for the TrackSeedContainer
+
+  unsigned int m_randomSeed{0};  // the default random seed
+  gsl_rng* m_rng{nullptr};
+};
+
+#endif  // PHSILICONSEEDPRUNER_H