Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-06 08:18:44

0001 /// ===========================================================================
0002 /*! \file   TrksInJetQAInJetFiller.cc
0003  *  \author Derek Anderson
0004  *  \date   04.03.2024
0005  *
0006  *  A submodule for the TrksInJetsQAM F4A module to produce
0007  *  QA histograms for tracks and more in jets
0008  */
0009 /// ===========================================================================
0010 
0011 #define TRKSINJETQAINJETFILLER_CC
0012 
0013 // submodule definition
0014 #include "TrksInJetQAInJetFiller.h"
0015 
0016 // inherited public methods ===================================================
0017 
0018 // ----------------------------------------------------------------------------
0019 //! Run fill for relevant histograms
0020 // ----------------------------------------------------------------------------
0021 void TrksInJetQAInJetFiller::Fill(PHCompositeNode* topNode)
0022 {
0023   GetNodes(topNode);
0024 
0025   FillJetAndTrackQAHists(topNode);
0026 }  // end 'Fill(PHCompositeNode* topNode)'
0027 
0028 // private methods ============================================================
0029 
0030 // ----------------------------------------------------------------------------
0031 //! Grab additional input nodes
0032 // ----------------------------------------------------------------------------
0033 void TrksInJetQAInJetFiller::GetNode(const int node, PHCompositeNode* topNode)
0034 {
0035   // jump to relevant node
0036   // NOLINTNEXTLINE(hicpp-multiway-paths-covered)
0037   switch (node)
0038   {
0039     case Node::Flow:
0040       m_flowStore = findNode::getClass<ParticleFlowElementContainer>(topNode, "ParticleFlowElements");
0041       if (!m_flowStore)
0042       {
0043         std::cerr << PHWHERE << ": PANIC: Couldn't grab particle flow container from node tree!" << std::endl;
0044         assert(m_flowStore);
0045       }
0046       break;
0047     default:
0048       std::cerr << PHWHERE << ": WARNING: trying to grab unkown additional node..." << std::endl;
0049       break;
0050   }
0051 }  // end 'GetNode(int, PHCompositeNode*)'
0052 
0053 // ----------------------------------------------------------------------------
0054 //! Fill histograms for jet and in-jet track histograms
0055 // ----------------------------------------------------------------------------
0056 /*! This defines how to (1) loop over jets, and (2) extract
0057  *  tracks inside them. Some care is required here since
0058  *  jets can be made
0059  *    a. completely without tracks (e.g. from calo towers),
0060  *    b. completely with tracks (i.e. track jets),
0061  *    c. or with a mix of tracks and other objects (e.g.
0062  *       particle flow jets).
0063  *  Once tracks have been identified, the relevant hit and
0064  *  cluster populations are extracted from the tracks, and
0065  *  their histograms are filled.
0066  */
0067 void TrksInJetQAInJetFiller::FillJetAndTrackQAHists(PHCompositeNode* topNode)
0068 {
0069   // loop over jets
0070   for (
0071       uint64_t iJet = 0;
0072       iJet < m_jetMap->size();
0073       ++iJet)
0074   {
0075     // grab jet and make sure track vector is clear
0076     Jet* jet = m_jetMap->get_jet(iJet);
0077     m_trksInJet.clear();
0078 
0079     // get all tracks "in" a jet
0080     GetCstTracks(jet, topNode);
0081     GetNonCstTracks(jet);
0082 
0083     // grab jet info and fill histograms
0084     if (m_config.doJetQA)
0085     {
0086       m_jetManager->GetInfo(jet, m_trksInJet);
0087     }
0088 
0089     // loop over tracks in the jet
0090     for (SvtxTrack* track : m_trksInJet)
0091     {
0092       // grab track info and fill histograms
0093       if (m_config.doTrackQA)
0094       {
0095         m_trackManager->GetInfo(track, jet);
0096       }
0097 
0098       // fill cluster and hit histograms as needed
0099       if (m_config.doClustQA || m_config.doHitQA)
0100       {
0101         FillClustAndHitQAHists(track);
0102       }
0103     }  // end track loop
0104   }  // end jet loop
0105 }  // end 'FillJetAndTrackQAHists(PHCompositeNode*)'
0106 
0107 // ----------------------------------------------------------------------------
0108 //! Fill histograms for in-jet track hits and clusters
0109 // ----------------------------------------------------------------------------
0110 void TrksInJetQAInJetFiller::FillClustAndHitQAHists(SvtxTrack* track)
0111 {
0112   // get cluster keys
0113   for (auto clustKey : ClusKeyIter(track))
0114   {
0115     // grab cluster and its info
0116     if (m_config.doClustQA)
0117     {
0118       m_clustManager->GetInfo(
0119           m_clustMap->findCluster(clustKey),
0120           clustKey,
0121           m_actsGeom);
0122     }
0123 
0124     // get hits if needed
0125     if (m_config.doHitQA)
0126     {
0127       // grab hit set and key associated with cluster key
0128       TrkrDefs::hitsetkey setKey = TrkrDefs::getHitSetKeyFromClusKey(clustKey);
0129       TrkrHitSet* set = m_hitMap->findHitSet(setKey);
0130       if (!set || !(set->size() > 0))
0131       {
0132         return;
0133       }
0134 
0135       // loop over all hits in hit set
0136       TrkrHitSet::ConstRange hits = set->getHits();
0137       for (
0138           TrkrHitSet::ConstIterator itHit = hits.first;
0139           itHit != hits.second;
0140           ++itHit)
0141       {
0142         // grab hit
0143         TrkrDefs::hitkey hitKey = itHit->first;
0144         TrkrHit* hit = itHit->second;
0145 
0146         // grab info and fill histograms
0147         m_hitManager->GetInfo(hit, setKey, hitKey);
0148 
0149       }  // end hit loop
0150     }
0151   }  // end cluster key loop
0152 }  // end 'FillClustQAHists(SvtxTrack*)'
0153 
0154 // ----------------------------------------------------------------------------
0155 //! Identify constituent tracks in a jet
0156 // ----------------------------------------------------------------------------
0157 void TrksInJetQAInJetFiller::GetCstTracks(Jet* jet, PHCompositeNode* topNode)
0158 {
0159   // loop over consituents
0160   auto csts = jet->get_comp_vec();
0161   for (auto& cst : csts)
0162   {
0163     // ignore cst if non-relevent type
0164     const uint32_t src = cst.first;
0165     if (IsCstNotRelevant(src))
0166     {
0167       continue;
0168     }
0169 
0170     // if cst is track, add to list
0171     if (src == Jet::SRC::TRACK)
0172     {
0173       m_trksInJet.push_back(m_trkMap->get(cst.second));
0174     }
0175 
0176     // if pfo, grab track if needed
0177     if (src == Jet::SRC::PARTICLE)
0178     {
0179       TrksInJetQADefs::PFObject* pfo = GetPFObject(cst.second, topNode);
0180       SvtxTrack* track = GetTrkFromPFO(pfo);
0181       if (track)
0182       {
0183         m_trksInJet.push_back(track);
0184       }
0185     }
0186   }  // end cst loop
0187 }  // end 'GetCstTracks(Jet* jet, PHCompositeNode* topNode)'
0188 
0189 // ----------------------------------------------------------------------------
0190 //! Identify tracks in jet cone that are NOT constituents
0191 // ----------------------------------------------------------------------------
0192 /*! This method comes into play if, e.g., the jets being
0193  *  analyzed were made with calorimeter information only.
0194  *  In this case you can still have tracks lying inside
0195  *  the jet cone despite them not being consituents of
0196  *  the constructed jet.
0197  */
0198 void TrksInJetQAInJetFiller::GetNonCstTracks(Jet* jet)
0199 {
0200   // loop over tracks
0201   for (auto& itTrk : *m_trkMap)
0202   {
0203     // grab track
0204     SvtxTrack* track = itTrk.second;
0205 
0206     // ignore tracks we've already added to the list
0207     if (IsTrkInList(track->get_id()))
0208     {
0209       continue;
0210     }
0211 
0212     // FIXME this can be improved!
0213     //   - jets don't necessarily have areas of
0214     //     pi*(Rjet)^2
0215     //   - it may be better to instead check
0216     //     if a track projection falls in
0217     //     a constituent tower/cluster
0218     //   - Also track projections to a calo
0219     //     would be better to use than just
0220     //     the track
0221     const double dr = GetTrackJetDist(track, jet);
0222     if (dr < m_config.rJet)
0223     {
0224       m_trksInJet.push_back(track);
0225     }
0226   }  // end track loop
0227 }  // end 'GetNonCstTracks(Jet* jet)'
0228 
0229 // ----------------------------------------------------------------------------
0230 //! Check if a consituent is not a relevant type
0231 // ----------------------------------------------------------------------------
0232 bool TrksInJetQAInJetFiller::IsCstNotRelevant(const uint32_t type)
0233 {
0234   const bool isVoid = (type == Jet::SRC::VOID);
0235   const bool isImport = (type == Jet::SRC::HEPMC_IMPORT);
0236   const bool isProbe = (type == Jet::SRC::JET_PROBE);
0237   return (isVoid || isImport || isProbe);
0238 }  // end 'IsCstNotRelevant(uint32_t)'
0239 
0240 // ----------------------------------------------------------------------------
0241 //! Check if a track has already been added to the list of tracks in the jet
0242 // ----------------------------------------------------------------------------
0243 bool TrksInJetQAInJetFiller::IsTrkInList(const uint32_t id)
0244 {
0245   bool isAdded = false;
0246   for (SvtxTrack* trkInJet : m_trksInJet)
0247   {
0248     if (id == trkInJet->get_id())
0249     {
0250       isAdded = true;
0251       break;
0252     }
0253   }
0254   return isAdded;
0255 }  // end 'IsTrkInList(uint32_t)'
0256 
0257 // ----------------------------------------------------------------------------
0258 //! Get eta-phi distance between track & jet axis
0259 // ----------------------------------------------------------------------------
0260 double TrksInJetQAInJetFiller::GetTrackJetDist(SvtxTrack* track, Jet* jet)
0261 {
0262   // get delta eta
0263   const double dEta = (track->get_eta()) - (jet->get_eta());
0264 
0265   // get delta phi
0266   double dPhi = (track->get_phi()) - (jet->get_phi());
0267   if (dPhi < (-1. * M_PI))
0268   {
0269     dPhi += 2. * M_PI;
0270   }
0271   if (dPhi > (1. * M_PI))
0272   {
0273     dPhi -= 2. * M_PI;
0274   }
0275 
0276   // return distance
0277   return std::hypot(dEta, dPhi);
0278 
0279 }  // end 'GetTrackJetDist(SvtxTrack*, Jet*)'
0280 
0281 // ----------------------------------------------------------------------------
0282 //! Find a particle flow element in node
0283 // ----------------------------------------------------------------------------
0284 TrksInJetQADefs::PFObject* TrksInJetQAInJetFiller::GetPFObject(const uint32_t id,
0285                                                                PHCompositeNode* topNode)
0286 {
0287   // pointer to pfo
0288   TrksInJetQADefs::PFObject* pfoToFind = nullptr;
0289 
0290   // grab pf node if needed
0291   if (!m_flowStore)
0292   {
0293     GetNode(Node::Flow, topNode);
0294   }
0295 
0296   // loop over pfos
0297   for (
0298       ParticleFlowElementContainer::ConstIterator itFlow = m_flowStore->getParticleFlowElements().first;
0299       itFlow != m_flowStore->getParticleFlowElements().second;
0300       ++itFlow)
0301   {
0302     // get pfo
0303     TrksInJetQADefs::PFObject* pfo = itFlow->second;
0304 
0305     // if has provided id, set pointer and exit
0306     if (id == pfo->get_id())
0307     {
0308       pfoToFind = pfo;
0309       break;
0310     }
0311   }  // end pfo loop
0312   return pfoToFind;
0313 }  // end 'GetPFObject(uint32_t, PHCompositeNode*)'
0314 
0315 // ----------------------------------------------------------------------------
0316 //! Get track from a particle flow element
0317 // ----------------------------------------------------------------------------
0318 SvtxTrack* TrksInJetQAInJetFiller::GetTrkFromPFO(TrksInJetQADefs::PFObject* pfo)
0319 {
0320   // pointer to track
0321   SvtxTrack* track = nullptr;
0322 
0323   // if pfo has track, try to grab
0324   const auto type = pfo->get_type();
0325   if (
0326       (type == ParticleFlowElement::PFLOWTYPE::MATCHED_CHARGED_HADRON) ||
0327       (type == ParticleFlowElement::PFLOWTYPE::UNMATCHED_CHARGED_HADRON))
0328   {
0329     track = pfo->get_track();
0330   }
0331   return track;
0332 }  // end 'GetTrkFromPFO(TrksInJetQADefs::PFObject*)'
0333 
0334 // end ========================================================================