g4simulation/g4eval/PHG4DstCompressReco.cc

0001 #include "PHG4DstCompressReco.h"
0002
0003 #include <g4main/PHG4Hit.h>
0004 #include <g4main/PHG4HitContainer.h>
0005 #include <g4main/PHG4Particle.h>
0006 #include <g4main/PHG4Shower.h>
0007 #include <g4main/PHG4TruthInfoContainer.h>
0008
0009 #include <g4detectors/PHG4CellContainer.h>
0010
0011 #include <calobase/RawTower.h>
0012 #include <calobase/RawTowerContainer.h>
0013
0014 #include <trackbase_historic/PHG4ParticleSvtxMap.h>
0015 #include <trackbase_historic/SvtxPHG4ParticleMap.h>
0016
0017 #include <fun4all/Fun4AllReturnCodes.h>
0018 #include <fun4all/SubsysReco.h>
0019
0020 #include <phool/PHCompositeNode.h>
0021 #include <phool/PHIODataNode.h>
0022 #include <phool/PHNode.h>
0023 #include <phool/PHNodeIterator.h>
0024 #include <phool/PHPointerListIterator.h>
0025 #include <phool/getClass.h>
0026
0027 #include <iostream>
0028 #include <utility>
0029
0030 PHG4DstCompressReco::PHG4DstCompressReco(const std::string& name)
0031   : SubsysReco(name)
0032 {
0033 }
0034
0035 int PHG4DstCompressReco::InitRun(PHCompositeNode* topNode)
0036 {
0037   _truth_info = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0038   if (!_truth_info)
0039   {
0040     std::cout << "PHG4DstCompressReco::InitRun(): Can't find G4TruthInfo" << std::endl;
0041     return Fun4AllReturnCodes::ABORTRUN;
0042   }
0043
0044   SearchG4HitNodes(topNode);
0045
0046   for (const auto& name : _compress_g4cell_names)
0047   {
0048     PHG4CellContainer* g4cells = findNode::getClass<PHG4CellContainer>(topNode, name);
0049     if (g4cells)
0050     {
0051       _g4cells.insert(g4cells);
0052     }
0053   }
0054
0055   for (const auto& name : _compress_tower_names)
0056   {
0057     RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(topNode, name);
0058     if (towers)
0059     {
0060       _towers.insert(towers);
0061     }
0062   }
0063
0064   if (m_keepRecoTrackMatchedParticles)
0065   {
0066     _recoTruthMap = findNode::getClass<SvtxPHG4ParticleMap>(topNode, "SvtxPHG4ParticleMap");
0067     if (_recoTruthMap == nullptr)
0068     {
0069       std::cout << __PRETTY_FUNCTION__ << " Fatal error: missing SvtxPHG4ParticleMap while m_keepRecoTrackMatchedParticles is set. "
0070                 << "Was PHG4DstCompressReco called before this module?"
0071                 << std::endl;
0072       exit(1);
0073     }
0074
0075     _truthRecoMap = findNode::getClass<PHG4ParticleSvtxMap>(topNode, "PHG4ParticleSvtxMap");
0076     if (_truthRecoMap == nullptr)
0077     {
0078       std::cout << __PRETTY_FUNCTION__ << " Fatal error: missing PHG4ParticleSvtxMap while m_keepRecoTrackMatchedParticles is set. "
0079                 << "Was PHG4DstCompressReco called before this module?"
0080                 << std::endl;
0081       exit(1);
0082     }
0083   }  //  if (m_keepRecoTrackMatchedParticles)
0084
0085   return Fun4AllReturnCodes::EVENT_OK;
0086 }
0087
0088 int PHG4DstCompressReco::process_event(PHCompositeNode* /*topNode*/)
0089 {
0090   if (_g4hits.empty() && _g4cells.empty() && _towers.empty())
0091   {
0092     return Fun4AllReturnCodes::EVENT_OK;
0093   }
0094
0095   //---cells--------------------------------------------------------------------
0096
0097   for (auto* cells : _g4cells)
0098   {
0099     cells->Reset();  // DROP ALL COMPRESSED G4CELLS
0100   }
0101
0102   //---hits---------------------------------------------------------------------
0103
0104   for (auto* hits : _g4hits)
0105   {
0106     hits->Reset();  // DROP ALL COMPRESSED G4HITS
0107   }
0108
0109   //---secondary particles and vertexes-----------------------------------------
0110
0111   std::set<int> keep_particle_ids;
0112   for (auto* hits : _keep_g4hits)
0113   {
0114     for (PHG4HitContainer::ConstIterator jter = hits->getHits().first;
0115          jter != hits->getHits().second;
0116          ++jter)
0117     {
0118       PHG4Hit* hit = jter->second;
0119       keep_particle_ids.insert(hit->get_trkid());
0120       // this will need to include all parents too in a trace back to
0121       // the primary, but let's start here for now
0122     }
0123   }
0124
0125   //---tracker truth map, if set-----------------------------------------
0126   if (_truthRecoMap)
0127   {
0128     for (const auto& [particle_id, map] : *_truthRecoMap)
0129     {
0130       keep_particle_ids.insert(particle_id);
0131     }
0132   }
0133   if (_recoTruthMap)
0134   {
0135     for (const auto& [track_id, weighted_truth_track_map] : *_recoTruthMap)
0136     {
0137       for (const auto& [weight, particle_set] : weighted_truth_track_map)
0138       {
0139         if (weight > 0)
0140         {
0141           for (const auto& particle_id : particle_set)
0142           {
0143             keep_particle_ids.insert(particle_id);
0144           }
0145         }
0146       }
0147     }  //    for (const auto& [track_id, weighted_truth_track_map] : *_recoTruthMap)
0148
0149   }  //  if (_recoTruthMap)
0150
0151   std::set<int> keep_vertex_ids;
0152   PHG4TruthInfoContainer::Range range = _truth_info->GetSecondaryParticleRange();
0153   for (PHG4TruthInfoContainer::Iterator iter = range.first;
0154        iter != range.second;)
0155   {
0156     int id = iter->first;
0157     PHG4Particle* particle = iter->second;
0158
0159     if (keep_particle_ids.contains(id))
0160     {
0161       ++iter;
0162       keep_vertex_ids.insert(particle->get_vtx_id());
0163       continue;
0164     }
0165
0166     _truth_info->delete_particle(iter++);  // DROP PARTICLES NOT ASSOCIATED TO A PRESERVED HIT
0167   }
0168
0169   PHG4TruthInfoContainer::VtxRange vrange = _truth_info->GetSecondaryVtxRange();
0170   for (PHG4TruthInfoContainer::VtxIterator iter = vrange.first;
0171        iter != vrange.second;)
0172   {
0173     int id = iter->first;
0174
0175     if (keep_vertex_ids.contains(id))
0176     {
0177       ++iter;
0178       continue;
0179     }
0180
0181     _truth_info->delete_vtx(iter++);  // DROP VERTEXES NOT ASSOCIATED TO A PRESERVED HIT
0182   }
0183
0184   //---shower entries-----------------------------------------------------------
0185
0186   PHG4TruthInfoContainer::ShowerRange srange = _truth_info->GetShowerRange();
0187   for (PHG4TruthInfoContainer::ShowerIterator iter = srange.first;
0188        iter != srange.second;
0189        ++iter)
0190   {
0191     PHG4Shower* shower = iter->second;
0192
0193     shower->clear_g4particle_id();
0194     shower->clear_g4vertex_id();
0195     shower->clear_g4hit_id();
0196   }
0197
0198   //---tower cell entries-------------------------------------------------------
0199   for (auto* towers : _towers)
0200   {
0201     // loop over all the towers
0202     for (RawTowerContainer::Iterator jter = towers->getTowers().first;
0203          jter != towers->getTowers().second;
0204          ++jter)
0205     {
0206       RawTower* tower = jter->second;
0207       tower->clear_g4cells();
0208     }
0209   }
0210
0211   return Fun4AllReturnCodes::EVENT_OK;
0212 }
0213
0214 void PHG4DstCompressReco::SearchG4HitNodes(PHCompositeNode* top)
0215 {
0216   // fill a lookup map between the g4hit container ids and the containers
0217   // themselves
0218   // without knowing what the container names are in advance, only that they
0219   // begin G4HIT_*
0220
0221   // separate the names into those in the compression list and those not in the
0222   // compression list
0223
0224   PHNodeIterator nodeiter(top);
0225   PHPointerListIterator<PHNode> iter(nodeiter.ls());
0226   PHNode* thisNode;
0227   while ((thisNode = iter()))
0228   {
0229     if (thisNode->getType() == "PHCompositeNode")
0230     {
0231       SearchG4HitNodes(static_cast<PHCompositeNode*>(thisNode));
0232     }
0233     else if (thisNode->getType() == "PHIODataNode")
0234     {
0235       if (thisNode->getName().find("G4HIT_") == 0)
0236       {
0237         PHIODataNode<PHG4HitContainer>* DNode =
0238             static_cast<PHIODataNode<PHG4HitContainer>*>(thisNode);
0239         if (DNode)
0240         {
0241           PHG4HitContainer* object =
0242               dynamic_cast<PHG4HitContainer*>(DNode->getData());
0243           if (object)
0244           {
0245             if (_compress_g4hit_names.contains(thisNode->getName()))
0246             {
0247               _g4hits.insert(object);
0248             }
0249             else
0250             {
0251               _keep_g4hits.insert(object);
0252             }
0253           }
0254         }
0255       }
0256     }
0257   }
0258 }