QA/SimulationModules/QAG4SimulationMicromegas.cc

0001 #include "QAG4SimulationMicromegas.h"
0002
0003 #include <qautils/QAHistManagerDef.h>
0004
0005 #include <g4detectors/PHG4CylinderGeomContainer.h>
0006
0007 #include <g4main/PHG4Hit.h>
0008 #include <g4main/PHG4HitContainer.h>
0009 #include <g4main/PHG4Particle.h>
0010 #include <g4main/PHG4TruthInfoContainer.h>
0011
0012 #include <micromegas/CylinderGeomMicromegas.h>
0013 #include <micromegas/MicromegasDefs.h>  // for SegmentationType
0014
0015 #include <trackbase_historic/ActsTransformations.h>
0016
0017 #include <trackbase/ActsGeometry.h>
0018 #include <trackbase/ClusterErrorPara.h>
0019 #include <trackbase/TrackFitUtils.h>
0020 #include <trackbase/TrkrCluster.h>
0021 #include <trackbase/TrkrClusterContainer.h>
0022 #include <trackbase/TrkrClusterHitAssoc.h>
0023 #include <trackbase/TrkrDefs.h>  // for getTrkrId, getHit...
0024 #include <trackbase/TrkrHit.h>
0025 #include <trackbase/TrkrHitSet.h>
0026 #include <trackbase/TrkrHitSetContainer.h>
0027 #include <trackbase/TrkrHitTruthAssoc.h>
0028
0029 #include <g4eval/SvtxClusterEval.h>  // for SvtxClusterEval
0030 #include <g4eval/SvtxEvalStack.h>
0031 #include <g4eval/SvtxTruthEval.h>  // for SvtxTruthEval
0032
0033 #include <fun4all/Fun4AllHistoManager.h>
0034 #include <fun4all/Fun4AllReturnCodes.h>
0035 #include <fun4all/SubsysReco.h>  // for SubsysReco
0036
0037 #include <phool/getClass.h>
0038 #include <phool/phool.h>  // for PHWHERE
0039
0040 #include <TAxis.h>  // for TAxis
0041 #include <TH1.h>
0042 #include <TVector3.h>
0043
0044 #include <cassert>
0045 #include <cmath>  // for cos, sin
0046 #include <format>
0047 #include <iostream>  // for operator<<, basic...
0048 #include <iterator>  // for distance
0049 #include <limits>
0050 #include <map>  // for map
0051 #include <string>
0052 #include <utility>  // for pair, make_pair
0053 #include <vector>   // for vector
0054
0055 //_____________________________________________________________________
0056 namespace
0057 {
0058   //! square
0059   template <class T>
0060   constexpr T square(T x)
0061   {
0062     return x * x;
0063   }
0064
0065   //! needed for weighted linear interpolation
0066   struct interpolation_data_t
0067   {
0068     using list = std::vector<interpolation_data_t>;
0069     double x() const { return position.x(); }
0070     double y() const { return position.y(); }
0071     double z() const { return position.z(); }
0072
0073     double px() const { return momentum.x(); }
0074     double py() const { return momentum.y(); }
0075     double pz() const { return momentum.z(); }
0076
0077     // NOLINTNEXTLINE(misc-non-private-member-variables-in-classes)
0078     TVector3 position;
0079     // NOLINTNEXTLINE(misc-non-private-member-variables-in-classes)
0080     TVector3 momentum;
0081     // NOLINTNEXTLINE(misc-non-private-member-variables-in-classes)
0082     double weight = 1;
0083   };
0084
0085   //! calculate the interpolation of member function called on all members in collection to the provided y_extrap
0086   template <double (interpolation_data_t::*accessor)() const>
0087   double interpolate(const interpolation_data_t::list& hits, double z_extrap)
0088   {
0089     // calculate all terms needed for the interpolation
0090     // need to use double everywhere here due to numerical divergences
0091     double sw = 0;
0092     double swz = 0;
0093     double swz2 = 0;
0094     double swx = 0;
0095     double swzx = 0;
0096
0097     bool valid(false);
0098     for (const auto& hit : hits)
0099     {
0100       const double x = (hit.*accessor)();
0101       const double w = hit.weight;
0102       if (w <= 0)
0103       {
0104         continue;
0105       }
0106
0107       valid = true;
0108       const double z = hit.z();
0109
0110       sw += w;
0111       swz += w * z;
0112       swz2 += w * square(z);
0113       swx += w * x;
0114       swzx += w * x * z;
0115     }
0116
0117     if (!valid)
0118     {
0119       return std::numeric_limits<double>::quiet_NaN();
0120     }
0121
0122     const auto alpha = (sw * swzx - swz * swx);
0123     const auto beta = (swz2 * swx - swz * swzx);
0124     const auto denom = (sw * swz2 - square(swz));
0125
0126     return (alpha * z_extrap + beta) / denom;
0127   }
0128
0129 }  // namespace
0130
0131 //________________________________________________________________________
0132 QAG4SimulationMicromegas::QAG4SimulationMicromegas(const std::string& name)
0133   : SubsysReco(name)
0134   , m_truthContainer(nullptr)
0135 {
0136 }
0137
0138 //________________________________________________________________________
0139 int QAG4SimulationMicromegas::InitRun(PHCompositeNode* topNode)
0140 {
0141   // prevent multiple creations of histograms
0142   if (m_initialized)
0143   {
0144     return Fun4AllReturnCodes::EVENT_OK;
0145   }
0146
0147   m_initialized = true;
0148
0149   // find mvtx geometry
0150   auto* geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
0151   if (!geom_container)
0152   {
0153     std::cout << PHWHERE << " unable to find DST node CYLINDERGEOM_MICROMEGAS_FULL" << std::endl;
0154     return Fun4AllReturnCodes::ABORTRUN;
0155   }
0156
0157   if (!m_svtxEvalStack)
0158   {
0159     m_svtxEvalStack.reset(new SvtxEvalStack(topNode));
0160     //    m_svtxEvalStack->set_strict(true);
0161     m_svtxEvalStack->set_strict(false);
0162     m_svtxEvalStack->set_verbosity(Verbosity());
0163   }
0164
0165   m_truthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode,
0166                                                                 "G4TruthInfo");
0167   if (!m_truthContainer)
0168   {
0169     std::cout << "QAG4SimulationTpc::InitRun - Fatal Error - "
0170               << "unable to find node G4TruthInfo" << std::endl;
0171     assert(m_truthContainer);
0172   }
0173
0174   // get layers from mvtx geometry
0175   const auto range = geom_container->get_begin_end();
0176   for (auto iter = range.first; iter != range.second; ++iter)
0177   {
0178     m_layers.insert(iter->first);
0179   }
0180
0181   // histogram manager
0182   auto* hm = QAHistManagerDef::getHistoManager();
0183   assert(hm);
0184
0185   // create histograms
0186   for (const auto& layer : m_layers)
0187   {
0188     if (Verbosity())
0189     {
0190       std::cout << PHWHERE << " adding layer " << layer << std::endl;
0191     }
0192
0193     // get layer geometry
0194     auto* const layergeom = dynamic_cast<CylinderGeomMicromegas*>(geom_container->GetLayerGeom(layer));
0195     assert(layergeom);
0196
0197     // get segmentation type
0198     const auto segmentation_type = layergeom->get_segmentation_type();
0199     const bool is_segmentation_phi = (segmentation_type == MicromegasDefs::SegmentationType::SEGMENTATION_PHI);
0200
0201     {
0202       // ADC distributions
0203       auto* h = new TH1F(std::format("{}adc_{}", get_histo_prefix(), layer).c_str(),
0204                          std::format("micromegas ADC distribution layer_{}", layer).c_str(),
0205                          1024, 0, 1024);
0206       h->GetXaxis()->SetTitle("ADC");
0207       hm->registerHisto(h);
0208     }
0209
0210     {
0211       // residuals (cluster - truth)
0212       const double max_residual = is_segmentation_phi ? 0.04 : 0.08;
0213       auto* h = new TH1F(std::format("{}residual_{}", get_histo_prefix(), layer).c_str(),
0214                          std::format("micromegas {}_{{cluster-truth}} layer_{}", (is_segmentation_phi ? "r#Delta#phi" : "#Deltaz"), layer).c_str(),
0215                          100, -max_residual, max_residual);
0216       h->GetXaxis()->SetTitle(std::format("{}_{{cluster-truth}} (cm)", (is_segmentation_phi ? "r#Delta#phi" : "#Deltaz")).c_str());
0217       hm->registerHisto(h);
0218     }
0219
0220     {
0221       // cluster errors
0222       const double max_error = is_segmentation_phi ? 0.04 : 0.08;
0223       auto* h = new TH1F(std::format("{}residual_error_{}", get_histo_prefix(), layer).c_str(),
0224                          std::format("micromegas {} error layer_{}", (is_segmentation_phi ? "r#Delta#phi" : "#Deltaz"), layer).c_str(),
0225                          100, 0, max_error);
0226       h->GetXaxis()->SetTitle(std::format("{} error (cm)", (is_segmentation_phi ? "r#Delta#phi" : "#Deltaz")).c_str());
0227       hm->registerHisto(h);
0228     }
0229
0230     {
0231       // pulls (cluster - truth)
0232       auto* h = new TH1F(std::format("{}cluster_pulls_{}", get_histo_prefix(), layer).c_str(),
0233                          std::format("micromegas {} layer_{}", (is_segmentation_phi ? "#Delta#phi/#sigma#phi" : "#Deltaz/#sigmaz"), layer).c_str(),
0234                          100, -5, 5);
0235       h->GetXaxis()->SetTitle(is_segmentation_phi ? "#Delta#phi/#sigma#phi" : "#Deltaz/#sigmaz");
0236       hm->registerHisto(h);
0237     }
0238
0239     {
0240       // cluster size
0241       auto* h = new TH1F(std::format("{}clus_size_{}", get_histo_prefix(), layer).c_str(), std::format("micromegas cluster size layer_{}", layer).c_str(), 20, 0, 20);
0242       h->GetXaxis()->SetTitle("csize");
0243       hm->registerHisto(h);
0244     }
0245   }
0246
0247   return Fun4AllReturnCodes::EVENT_OK;
0248 }
0249
0250 //_____________________________________________________________________
0251 int QAG4SimulationMicromegas::process_event(PHCompositeNode* topNode)
0252 {
0253   // load nodes
0254   auto res = load_nodes(topNode);
0255   if (res != Fun4AllReturnCodes::EVENT_OK)
0256   {
0257     return res;
0258   }
0259   m_svtxEvalStack->next_event(topNode);
0260   // run evaluation
0261   evaluate_hits();
0262   evaluate_clusters();
0263   return Fun4AllReturnCodes::EVENT_OK;
0264 }
0265
0266 //________________________________________________________________________
0267 std::string QAG4SimulationMicromegas::get_histo_prefix() const
0268 {
0269   return std::string("h_") + Name() + std::string("_");
0270 }
0271
0272 //________________________________________________________________________
0273 int QAG4SimulationMicromegas::load_nodes(PHCompositeNode* topNode)
0274 {
0275   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0276   if (!m_tGeometry)
0277   {
0278     std::cout << PHWHERE << "No acts tracking geometry, exiting."
0279               << std::endl;
0280     return Fun4AllReturnCodes::ABORTEVENT;
0281   }
0282
0283   m_micromegas_geonode = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
0284   if (!m_micromegas_geonode)
0285   {
0286     std::cout << PHWHERE << " ERROR: Can't find CYLINDERGEOM_MICROMEGAS_FULL." << std::endl;
0287     return Fun4AllReturnCodes::ABORTEVENT;
0288   }
0289
0290   m_hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0291   if (!m_hitsets)
0292   {
0293     std::cout << PHWHERE << " ERROR: Can't find TrkrHitSetContainer." << std::endl;
0294   }
0295
0296   m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0297   if (!m_cluster_map)
0298   {
0299     std::cout << PHWHERE << " unable to find DST node TRKR_CLUSTER" << std::endl;
0300     return Fun4AllReturnCodes::ABORTEVENT;
0301   }
0302
0303   m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0304   if (!m_cluster_hit_map)
0305   {
0306     std::cout << PHWHERE << " unable to find DST node TRKR_CLUSTERHITASSOC" << std::endl;
0307     return Fun4AllReturnCodes::ABORTEVENT;
0308   }
0309
0310   m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
0311   if (!m_hit_truth_map)
0312   {
0313     std::cout << PHWHERE << " unable to find DST node TRKR_HITTRUTHASSOC" << std::endl;
0314     return Fun4AllReturnCodes::ABORTEVENT;
0315   }
0316
0317   m_g4hits_micromegas = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MICROMEGAS");
0318   if (!m_g4hits_micromegas)
0319   {
0320     std::cout << PHWHERE << " unable to find DST node G4HIT_MVTX" << std::endl;
0321     return Fun4AllReturnCodes::ABORTEVENT;
0322   }
0323
0324   return Fun4AllReturnCodes::EVENT_OK;
0325 }
0326
0327 //________________________________________________________________________
0328 void QAG4SimulationMicromegas::evaluate_hits()
0329 {
0330   // do nothing if hitsets are not there
0331   if (!m_hitsets)
0332   {
0333     return;
0334   }
0335
0336   // histogram manager
0337   auto* hm = QAHistManagerDef::getHistoManager();
0338   assert(hm);
0339
0340   // load relevant histograms
0341   struct HistogramList
0342   {
0343     TH1* adc = nullptr;
0344   };
0345
0346   using HistogramMap = std::map<int, HistogramList>;
0347   HistogramMap histograms;
0348
0349   for (const auto& layer : m_layers)
0350   {
0351     HistogramList h;
0352     h.adc = dynamic_cast<TH1*>(hm->getHisto(std::format("{}adc_{}", get_histo_prefix(), layer)));
0353     histograms.insert(std::make_pair(layer, h));
0354   }
0355
0356   // loop over hitsets
0357   const auto hitsetrange = m_hitsets->getHitSets(TrkrDefs::TrkrId::micromegasId);
0358   for (auto hitsetiter = hitsetrange.first; hitsetiter != hitsetrange.second; ++hitsetiter)
0359   {
0360     // get hitsetkey and layer
0361     const auto hitsetkey = hitsetiter->first;
0362     const auto layer = TrkrDefs::getLayer(hitsetkey);
0363
0364     // get relevant histograms
0365     const auto hiter = histograms.find(layer);
0366     if (hiter == histograms.end())
0367     {
0368       continue;
0369     }
0370
0371     // get all of the hits from this hitset
0372     TrkrHitSet* hitset = hitsetiter->second;
0373     TrkrHitSet::ConstRange const hitrange = hitset->getHits();
0374
0375     // loop over hits
0376     for (auto hit_it = hitrange.first; hit_it != hitrange.second; ++hit_it)
0377     {
0378       // store ADC
0379       auto fill = [](TH1* h, float value)
0380       { if( h ) { h->Fill( value );
0381 } };
0382       fill(hiter->second.adc, hit_it->second->getAdc());
0383     }
0384   }
0385 }
0386
0387 //________________________________________________________________________
0388 void QAG4SimulationMicromegas::evaluate_clusters()
0389 {
0390   SvtxTruthEval* trutheval = m_svtxEvalStack->get_truth_eval();
0391   assert(trutheval);
0392   SvtxClusterEval* clustereval = m_svtxEvalStack->get_cluster_eval();
0393   assert(clustereval);
0394   // histogram manager
0395   auto* hm = QAHistManagerDef::getHistoManager();
0396   assert(hm);
0397   // load relevant histograms
0398   struct HistogramList
0399   {
0400     TH1* residual = nullptr;
0401     TH1* residual_error = nullptr;
0402     TH1* pulls = nullptr;
0403
0404     TH1* csize = nullptr;
0405   };
0406
0407   using HistogramMap = std::map<int, HistogramList>;
0408   HistogramMap histograms;
0409
0410   for (const auto& layer : m_layers)
0411   {
0412     HistogramList h;
0413     h.residual = dynamic_cast<TH1*>(hm->getHisto(std::format("{}residual_{}", get_histo_prefix(), layer)));
0414     h.residual_error = dynamic_cast<TH1*>(hm->getHisto(std::format("{}residual_error_{}", get_histo_prefix(), layer)));
0415     h.pulls = dynamic_cast<TH1*>(hm->getHisto(std::format("{}cluster_pulls_{}", get_histo_prefix(), layer)));
0416     h.csize = dynamic_cast<TH1*>(hm->getHisto(std::format("{}clus_size_{}", get_histo_prefix(), layer)));
0417
0418     histograms.insert(std::make_pair(layer, h));
0419   }
0420   // get primary truth particles
0421
0422   PHG4TruthInfoContainer::ConstRange const range = m_truthContainer->GetPrimaryParticleRange();  // only from primary particles
0423
0424   for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0425        iter != range.second;
0426        ++iter)
0427   {
0428     PHG4Particle* g4particle = iter->second;
0429     float const gtrackID = g4particle->get_track_id();
0430     float const gflavor = g4particle->get_pid();
0431     float const gembed = trutheval->get_embed(g4particle);
0432     float const gprimary = trutheval->is_primary(g4particle);
0433
0434     if (Verbosity() > 0)
0435     {
0436       std::cout << PHWHERE << " PHG4Particle ID " << gtrackID << " gembed " << gembed << " gflavor " << gflavor << " gprimary " << gprimary << std::endl;
0437     }
0438     const auto ckeyset = clustereval->all_clusters_from(g4particle);
0439     // Get the truth clusters from this particle
0440     const auto truth_clusters = trutheval->all_truth_clusters(g4particle);
0441
0442     // get circle fit parameters first
0443     TrackFitUtils::position_vector_t xy_pts;
0444     TrackFitUtils::position_vector_t rz_pts;
0445
0446     for (const auto& [gkey, gclus] : truth_clusters)
0447     {
0448       const auto layer = TrkrDefs::getLayer(gkey);
0449       if (layer < 7)
0450       {
0451         continue;
0452       }
0453
0454       float const gx = gclus->getX();
0455       float const gy = gclus->getY();
0456       float const gz = gclus->getZ();
0457
0458       xy_pts.emplace_back(gx, gy);
0459       rz_pts.emplace_back(std::sqrt(gx * gx + gy * gy), gz);
0460     }
0461
0462     // fit a circle through x,y coordinates
0463     const auto [R, X0, Y0] = TrackFitUtils::circle_fit_by_taubin(xy_pts);
0464
0465     // skip chain entirely if fit fails
0466     if (std::isnan(R))
0467     {
0468       continue;
0469     }
0470
0471     // process residuals and pulls
0472     // get reco clusters
0473
0474     for (const auto ckey : ckeyset)
0475     {
0476       const auto layer = TrkrDefs::getLayer(ckey);
0477       const auto detID = TrkrDefs::getTrkrId(ckey);
0478       if (detID != TrkrDefs::TrkrId::micromegasId)
0479       {
0480         continue;
0481       }
0482
0483       // get tileid
0484       const auto tileid = MicromegasDefs::getTileId(ckey);
0485
0486       // load geometry
0487       auto* const layergeom = dynamic_cast<CylinderGeomMicromegas*>(m_micromegas_geonode->GetLayerGeom(layer));
0488       if (!layergeom)
0489       {
0490         continue;
0491       }
0492
0493       // get relevant histograms
0494       const auto hiter = histograms.find(layer);
0495       if (hiter == histograms.end())
0496       {
0497         continue;
0498       }
0499
0500       // get cluster
0501       auto* const cluster = m_cluster_map->findCluster(ckey);
0502       const auto global = m_tGeometry->getGlobalPosition(ckey, cluster);
0503
0504       // get segmentation type
0505       const auto segmentation_type = MicromegasDefs::getSegmentationType(ckey);
0506
0507       // get relevant cluster information
0508       double const rphi_error = cluster->getRPhiError();
0509       double const z_error = cluster->getZError();
0510
0511       // convert cluster position to local tile coordinates
0512       const TVector3 cluster_world(global(0), global(1), global(2));
0513       const auto cluster_local = layergeom->get_local_from_world_coords(tileid, m_tGeometry, cluster_world);
0514
0515       // find associated g4hits
0516       const auto g4hits = find_g4hits(ckey);
0517       // convert hits to list of interpolation_data_t
0518       /* TODO: should actually use the same code as in TrackEvaluation */
0519       interpolation_data_t::list hits;
0520       for (const auto& g4hit : g4hits)
0521       {
0522         const auto weight = g4hit->get_edep();
0523         for (int i = 0; i < 2; ++i)
0524         {
0525           // convert position to local
0526           TVector3 const g4hit_world(g4hit->get_x(i), g4hit->get_y(i), g4hit->get_z(i));
0527           TVector3 const g4hit_local = layergeom->get_local_from_world_coords(tileid, m_tGeometry, g4hit_world);
0528
0529           // convert momentum to local
0530           TVector3 const momentum_world(g4hit->get_px(i), g4hit->get_py(i), g4hit->get_pz(i));
0531           TVector3 const momentum_local = layergeom->get_local_from_world_vect(tileid, m_tGeometry, momentum_world);
0532
0533           hits.push_back({.position = g4hit_local, .momentum = momentum_local, .weight = weight});
0534         }
0535       }
0536
0537       // do position interpolation along z
0538       const auto z_extrap = cluster_local.z();
0539       const TVector3 interpolation_local(
0540           interpolate<&interpolation_data_t::x>(hits, z_extrap),
0541           interpolate<&interpolation_data_t::y>(hits, z_extrap),
0542           interpolate<&interpolation_data_t::z>(hits, z_extrap));
0543
0544       // fill phi residuals, errors and pulls
0545       auto fill = [](TH1* h, float value)
0546       { if( h ) { h->Fill( value );
0547 } };
0548       switch (segmentation_type)
0549       {
0550       case MicromegasDefs::SegmentationType::SEGMENTATION_PHI:
0551       {
0552         const auto drphi = cluster_local.x() - interpolation_local.x();
0553         fill(hiter->second.residual, drphi);
0554         fill(hiter->second.residual_error, rphi_error);
0555         fill(hiter->second.pulls, drphi / rphi_error);
0556         break;
0557       }
0558
0559       case MicromegasDefs::SegmentationType::SEGMENTATION_Z:
0560       {
0561         const auto dz = cluster_local.y() - interpolation_local.y();
0562         fill(hiter->second.residual, dz);
0563         fill(hiter->second.residual_error, z_error);
0564         fill(hiter->second.pulls, dz / z_error);
0565         break;
0566       }
0567       }
0568
0569       // cluster size
0570       // get associated hits
0571       const auto hit_range = m_cluster_hit_map->getHits(ckey);
0572       fill(hiter->second.csize, std::distance(hit_range.first, hit_range.second));
0573     }
0574   }
0575 }
0576
0577 //_____________________________________________________________________
0578 QAG4SimulationMicromegas::G4HitSet QAG4SimulationMicromegas::find_g4hits(TrkrDefs::cluskey cluster_key) const
0579 {
0580   // find hitset associated to cluster
0581   G4HitSet out;
0582   const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
0583
0584   // loop over hits associated to clusters
0585   const auto range = m_cluster_hit_map->getHits(cluster_key);
0586   for (auto iter = range.first; iter != range.second; ++iter)
0587   {
0588     // hit key
0589     const auto& hit_key = iter->second;
0590
0591     // store hits to g4hit associations
0592     TrkrHitTruthAssoc::MMap g4hit_map;
0593     m_hit_truth_map->getG4Hits(hitset_key, hit_key, g4hit_map);
0594
0595     // find corresponding g4 hist
0596     for (auto& truth_iter : g4hit_map)
0597     {
0598       // g4hit key
0599       const auto g4hit_key = truth_iter.second.second;
0600
0601       // g4 hit
0602       PHG4Hit* g4hit = (TrkrDefs::getTrkrId(hitset_key) == TrkrDefs::micromegasId) ? m_g4hits_micromegas->findHit(g4hit_key) : nullptr;
0603
0604       // insert in set
0605       if (g4hit)
0606       {
0607         out.insert(g4hit);
0608       }
0609     }
0610   }
0611
0612   return out;
0613 }