Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:13:23 

0001 /// ---------------------------------------------------------------------------
0002 /*! \file   SLambdaJetHunter.ana.h
0003  *  \author Derek Anderson
0004  *  \date   01.25.2024
0005  *
0006  *  A minimal analysis module to find lambda-tagged
0007  *  jets in pythia events.
0008  */
0009 /// ---------------------------------------------------------------------------
0010 
0011 #pragma once
0012 
0013 // make common namespaces implicit
0014 using namespace std;
0015 using namespace fastjet;
0016 using namespace findNode;
0017 
0018 
0019 
0020 namespace SColdQcdCorrelatorAnalysis {
0021 
0022   // analysis methods =========================================================
0023 
0024   // --------------------------------------------------------------------------
0025   //! Grab event-level information
0026   // --------------------------------------------------------------------------
0027   void SLambdaJetHunter::GrabEventInfo(PHCompositeNode* topNode) {
0028 
0029     // print debug statement
0030     if (m_config.isDebugOn) {
0031       cout << "SLambdaJetHunter::GrabEventInfo(PHCompositeNode*) Grabbing event info" << endl;
0032     }
0033 
0034     m_vecSubEvts = Tools::GrabSubevents(topNode);
0035     m_genEvtInfo.SetInfo(topNode, m_config.isEmbed, m_vecSubEvts);
0036     return;
0037 
0038   }  // end 'GrabEventInfo(PHCompositeNode*)'
0039 
0040 
0041 
0042   // --------------------------------------------------------------------------
0043   //! Find all lambds in event
0044   // --------------------------------------------------------------------------
0045   void SLambdaJetHunter::FindLambdas(PHCompositeNode* topNode) {
0046 
0047     // print debug statement
0048     if (m_config.isDebugOn) {
0049       cout << "SLambdaJetHunter::FindLambdas(PHCompositeNode*) Finding all lambdas in event" << endl;
0050     }
0051 
0052     // loop over subevents
0053     for (const int subEvt : m_vecSubEvts) {
0054 
0055       // loop over particles
0056       HepMC::GenEvent* genEvt = Interfaces::GetGenEvent(topNode, subEvt);
0057       for (
0058         HepMC::GenEvent::particle_const_iterator hepPar = genEvt -> particles_begin();
0059         hepPar != genEvt -> particles_end();
0060         ++hepPar
0061       ) {
0062 
0063         // check if lambda and if found yet
0064         const bool isLambda = IsLambda( (*hepPar) -> pdg_id() );
0065         const bool isNew    = IsNewLambda( (*hepPar) -> barcode() );
0066         if (!isLambda || !isNew) continue;
0067 
0068         // grab info
0069         Types::ParInfo lambda( (*hepPar), subEvt );
0070 
0071         // check if good
0072         const bool isGood = IsGoodLambda(lambda);
0073         if (!isGood) continue;
0074 
0075         // if a new good lambda, add to output vector and lambda-jet map
0076         m_lambdaInfo.push_back(lambda);
0077         m_mapLambdaJetAssoc.insert(
0078           make_pair( (*hepPar) -> barcode(), -1 )
0079         );
0080 
0081       }  // end particle loop
0082     }  // end subevent loop
0083     return;
0084 
0085   }  // end 'FindLambdas(PHCompositeNode*)'
0086 
0087 
0088 
0089   // --------------------------------------------------------------------------
0090   //! Reconstruct jets in an event
0091   // --------------------------------------------------------------------------
0092   void SLambdaJetHunter::MakeJets(PHCompositeNode* topNode) {
0093 
0094     // print debug statement
0095     if (m_config.isDebugOn) {
0096       cout << "SLambdaJetHunter::MakeJets() Making jets" << endl;
0097     }
0098 
0099     // for fastjet
0100     vector<PseudoJet> vecPseudoJets;
0101 
0102     // loop over subevents
0103     for (const int subEvt : m_vecSubEvts) {
0104 
0105       // loop over particles
0106       HepMC::GenEvent* genEvt = Interfaces::GetGenEvent(topNode, subEvt);
0107       for (
0108         HepMC::GenEvent::particle_const_iterator hepPar = genEvt -> particles_begin();
0109         hepPar != genEvt -> particles_end();
0110         ++hepPar
0111       ) {
0112 
0113         // grab particle info
0114         Types::ParInfo particle(*hepPar, subEvt);
0115 
0116         // check if particle is good & final state
0117         const bool isParGood = IsGoodParticle(particle);
0118         if (!isParGood || !particle.IsFinalState()) continue;
0119 
0120         // make pseudojet
0121         PseudoJet pseudo(
0122           particle.GetPX(),
0123           particle.GetPY(),
0124           particle.GetPZ(),
0125           particle.GetEne()
0126         );
0127         pseudo.set_user_index(particle.GetBarcode());
0128         vecPseudoJets.push_back(pseudo);
0129 
0130       }  // end particle loop
0131     }  // end subevent loop
0132 
0133     // set jet definition
0134     JetDefinition definition(
0135       Const::MapStringOntoFJAlgo()[m_config.jetAlgo],
0136       m_config.rJet,
0137       Const::MapStringOntoFJRecomb()[m_config.jetRecomb],
0138       fastjet::Best
0139     );
0140 
0141     // run clustering
0142     ClusterSequence* sequence = new ClusterSequence(vecPseudoJets, definition);
0143 
0144     // grab jets and return
0145     m_vecFastJets = sequence -> inclusive_jets();
0146     return;
0147 
0148   }  // end 'MakeJets(PHCompositeNode*)'
0149 
0150 
0151 
0152   // --------------------------------------------------------------------------
0153   //! Set jet output information
0154   // --------------------------------------------------------------------------
0155   /* TODO move to a dedicated interface */
0156   void SLambdaJetHunter::CollectJetOutput(PHCompositeNode* topNode) {
0157 
0158     // print debug statement
0159     if (m_config.isDebugOn) {
0160       cout << "SLambdaJetHunter::CollectOutput() hunting down lambda-taggged jets" << endl;
0161     }
0162 
0163     // reserve space for fastjet output
0164     m_jetInfo.resize( m_vecFastJets.size() );
0165     m_cstInfo.resize( m_vecFastJets.size() );
0166 
0167     // loop over fastjet output
0168     for (size_t iJet = 0; iJet < m_vecFastJets.size(); iJet++) {
0169 
0170       // grab jet info
0171       m_jetInfo[iJet].SetInfo( m_vecFastJets[iJet] );
0172       m_jetInfo[iJet].SetJetID( iJet );
0173 
0174       // loop over constituents
0175       m_cstInfo[iJet].resize( m_vecFastJets[iJet].constituents().size() );
0176       for (size_t iCst = 0; iCst < m_vecFastJets[iJet].constituents().size(); iCst++) {
0177 
0178         // grab cst info
0179         m_cstInfo[iJet][iCst].SetInfo( m_vecFastJets[iJet].constituents()[iCst] );
0180 
0181         // get particle PDG
0182         HepMC::GenParticle* hepCst = Tools::GetHepMCGenParticleFromBarcode(m_cstInfo[iJet][iCst].GetCstID(), topNode);
0183         m_cstInfo[iJet][iCst].SetPID( hepCst -> pdg_id() );
0184 
0185         // run calculations
0186         const float pCst  = hypot( m_cstInfo[iJet][iCst].GetPX(), m_cstInfo[iJet][iCst].GetPY(), m_cstInfo[iJet][iCst].GetPZ() );
0187         const float pJet  = hypot( m_jetInfo[iJet].GetPX(), m_jetInfo[iJet].GetPY(), m_jetInfo[iJet].GetPZ() );
0188         const float dfCst = m_cstInfo[iJet][iCst].GetPhi() - m_jetInfo[iJet].GetPhi();
0189         const float dhCst = m_cstInfo[iJet][iCst].GetEta() - m_jetInfo[iJet].GetEta();
0190 
0191         // grab remaining cst info
0192         m_cstInfo[iJet][iCst].SetJetID( m_jetInfo[iJet].GetJetID() );
0193         m_cstInfo[iJet][iCst].SetZ( pCst / pJet );
0194         m_cstInfo[iJet][iCst].SetDR( hypot(dfCst, dhCst) );
0195 
0196       }  // end cst loop
0197     }  // end jet loop
0198     return;
0199 
0200   }  // end 'CollectJetOutput(PHCompositeNode*)'
0201 
0202 
0203 
0204   // --------------------------------------------------------------------------
0205   //! Associate lambdas to jets using specified method
0206   // --------------------------------------------------------------------------
0207   void SLambdaJetHunter::AssociateLambdasToJets(PHCompositeNode* topNode) {
0208 
0209     // print debug statement
0210     if (m_config.isDebugOn) {
0211       cout << "SLambdaJetHunter::AssociateLambdasToJets() associating found lambdas to jets" << endl;
0212     }
0213 
0214     // loop over lambdas
0215     for (Types::ParInfo& lambda : m_lambdaInfo) {
0216 
0217       optional<int> iAssocJet = nullopt;
0218       switch (m_config.associator) {
0219 
0220         case Associator::Barcode:
0221           iAssocJet = HuntLambdasByBarcode(lambda);
0222           break;
0223 
0224         case Associator::Decay:
0225           iAssocJet = HuntLambdasByDecayChain(lambda, topNode);
0226           break;
0227 
0228         case Associator::Distance:
0229           iAssocJet = HuntLambdasByDistance(lambda);
0230           break;
0231 
0232         default:
0233           iAssocJet = HuntLambdasByDistance(lambda);
0234           break;
0235 
0236       }  // end switch (associator)
0237 
0238       // assign associated jet
0239       if (iAssocJet.has_value()) {
0240         m_mapLambdaJetAssoc[lambda.GetBarcode()] = iAssocJet.value();
0241       } else {
0242         m_mapLambdaJetAssoc[lambda.GetBarcode()] = -1;
0243       }
0244     }  // end lambda loop
0245     return;
0246 
0247   }   // end 'AssociateLambdasToJets(PHCompositeNode*)'
0248 
0249 
0250 
0251   // --------------------------------------------------------------------------
0252   //! Fill output tree
0253   // --------------------------------------------------------------------------
0254   /* TODO move to a dedicated interface */
0255   void SLambdaJetHunter::FillOutputTree() {
0256 
0257     // print debug statement
0258     if (m_config.isDebugOn) {
0259       cout << "SLambdaJetHunter::FillOutputTree() collecting output information and filling tree" << endl;
0260     }
0261 
0262     // collect event info
0263     m_evtNJets       = m_jetInfo.size();
0264     m_evtNLambdas    = m_mapLambdaJetAssoc.size();
0265     m_evtNTaggedJets = GetNTaggedJets();
0266     m_evtNChrgPars   = m_genEvtInfo.GetNChrgPar();
0267     m_evtNNeuPars    = m_genEvtInfo.GetNNeuPar();
0268     m_evtSumEPar     = (m_genEvtInfo.GetESumChrg() + m_genEvtInfo.GetESumNeu());
0269     m_evtVtxX        = m_genEvtInfo.GetPartonA().GetVX();
0270     m_evtVtxY        = m_genEvtInfo.GetPartonA().GetVY();
0271     m_evtVtxZ        = m_genEvtInfo.GetPartonA().GetVZ();
0272 
0273     // collect parton info
0274     m_evtPartID = make_pair( m_genEvtInfo.GetPartonA().GetPID(), m_genEvtInfo.GetPartonB().GetPID() );
0275     m_evtPartPx = make_pair( m_genEvtInfo.GetPartonA().GetPX(),  m_genEvtInfo.GetPartonB().GetPX()  );
0276     m_evtPartPy = make_pair( m_genEvtInfo.GetPartonA().GetPY(),  m_genEvtInfo.GetPartonB().GetPY()  );
0277     m_evtPartPz = make_pair( m_genEvtInfo.GetPartonA().GetPZ(),  m_genEvtInfo.GetPartonB().GetPZ()  );
0278     m_evtPartE  = make_pair( m_genEvtInfo.GetPartonA().GetEne(), m_genEvtInfo.GetPartonB().GetEne() );
0279 
0280     // collect lambda information
0281     for (Types::ParInfo& lambda : m_lambdaInfo) {
0282 
0283       // collect general information
0284       m_lambdaID.push_back( lambda.GetBarcode() );
0285       m_lambdaPID.push_back( lambda.GetPID() );
0286       m_lambdaEmbedID.push_back( lambda.GetEmbedID() );
0287       m_lambdaE.push_back( lambda.GetEne() );
0288       m_lambdaPt.push_back( lambda.GetPT() );
0289       m_lambdaEta.push_back( lambda.GetEta() );
0290       m_lambdaPhi.push_back( lambda.GetPhi() );
0291 
0292       // collect associate jet information
0293       m_lambdaJetID.push_back( m_mapLambdaJetAssoc[lambda.GetBarcode()] );
0294       if (m_mapLambdaJetAssoc[lambda.GetBarcode()] > -1) {
0295         m_lambdaZ.push_back( GetLambdaAssocZ(lambda) );
0296         m_lambdaDr.push_back( GetLambdaAssocDr(lambda) );
0297       } else {
0298         m_lambdaZ.push_back( -1. );
0299         m_lambdaDr.push_back( -1. );
0300       }
0301     }
0302 
0303     // collect jet information
0304     for (Types::JetInfo& jet : m_jetInfo) {
0305       m_jetHasLambda.push_back( HasLambda(jet) );
0306       m_jetNCst.push_back( jet.GetNCsts() );
0307       m_jetID.push_back( jet.GetJetID() );
0308       m_jetE.push_back( jet.GetEne() );
0309       m_jetPt.push_back( jet.GetPT() );
0310       m_jetEta.push_back( jet.GetEta() );
0311       m_jetPhi.push_back( jet.GetPhi() );
0312     }  // end jet loop
0313 
0314     // collect cst information
0315     m_cstID.resize( m_cstInfo.size() );
0316     m_cstPID.resize( m_cstInfo.size() );
0317     m_cstJetID.resize( m_cstInfo.size() );
0318     m_cstEmbedID.resize( m_cstInfo.size() );
0319     m_cstZ.resize( m_cstInfo.size() );
0320     m_cstDr.resize( m_cstInfo.size() );
0321     m_cstE.resize( m_cstInfo.size() );
0322     m_cstPt.resize( m_cstInfo.size() );
0323     m_cstEta.resize( m_cstInfo.size() );
0324     m_cstPhi.resize( m_cstInfo.size() );
0325     for (size_t iJet = 0; iJet < m_cstInfo.size(); iJet++) {
0326       m_cstID[iJet].resize( m_cstInfo[iJet].size() );
0327       m_cstPID[iJet].resize( m_cstInfo[iJet].size() );
0328       m_cstJetID[iJet].resize( m_cstInfo[iJet].size() );
0329       m_cstEmbedID[iJet].resize( m_cstInfo[iJet].size() );
0330       m_cstZ[iJet].resize( m_cstInfo[iJet].size() );
0331       m_cstDr[iJet].resize( m_cstInfo[iJet].size() );
0332       m_cstE[iJet].resize( m_cstInfo[iJet].size() );
0333       m_cstPt[iJet].resize( m_cstInfo[iJet].size() );
0334       m_cstEta[iJet].resize( m_cstInfo[iJet].size() );
0335       m_cstPhi[iJet].resize( m_cstInfo[iJet].size() );
0336       for (size_t iCst = 0; iCst < m_cstInfo[iJet].size(); iCst++) {
0337         m_cstID[iJet][iCst]      = m_cstInfo[iJet][iCst].GetCstID();
0338         m_cstPID[iJet][iCst]     = m_cstInfo[iJet][iCst].GetPID();
0339         m_cstJetID[iJet][iCst]   = m_cstInfo[iJet][iCst].GetJetID();
0340         m_cstEmbedID[iJet][iCst] = m_cstInfo[iJet][iCst].GetEmbedID();
0341         m_cstZ[iJet][iCst]       = m_cstInfo[iJet][iCst].GetZ();
0342         m_cstDr[iJet][iCst]      = m_cstInfo[iJet][iCst].GetDR();
0343         m_cstE[iJet][iCst]       = m_cstInfo[iJet][iCst].GetEne();
0344         m_cstPt[iJet][iCst]      = m_cstInfo[iJet][iCst].GetPT();
0345         m_cstEta[iJet][iCst]     = m_cstInfo[iJet][iCst].GetEta();
0346         m_cstPhi[iJet][iCst]     = m_cstInfo[iJet][iCst].GetPhi();
0347       }
0348     }  // end cst loop
0349 
0350     m_outTree -> Fill();
0351     return;
0352 
0353   }  // end 'FillOutputTree()'
0354 
0355 
0356 
0357   // --------------------------------------------------------------------------
0358   //! Check if a PHG4particle has a parent
0359   // --------------------------------------------------------------------------
0360   bool SLambdaJetHunter::HasParentInfo(const int parent) {
0361 
0362     // print debug statement
0363     if (m_config.isDebugOn && (m_config.verbosity > 5)) {
0364       cout << "SLambdaJetHunter::HasParentInfo(int) checking if PHG4particle has parent info" << endl;
0365     }
0366 
0367     return (parent != 0);
0368 
0369   }  // end 'HasParentInfo(int)'
0370 
0371 
0372 
0373   // --------------------------------------------------------------------------
0374   //! Check if a jet has a lambda associated to it
0375   // --------------------------------------------------------------------------
0376   bool SLambdaJetHunter::HasLambda(Types::JetInfo& jet) {
0377 
0378     // print debug statement
0379     if (m_config.isDebugOn && (m_config.verbosity > 5)) {
0380       cout << "SLambdaJetHunter::HasLambda(JetInfo&) checking if jet has associated lambda" << endl;
0381     }
0382 
0383     bool hasLambda = false;
0384     for (auto& [lamID, jetID]: m_mapLambdaJetAssoc) {
0385       if (jetID == (int) jet.GetJetID()) {
0386         hasLambda = true;
0387         break;
0388       }
0389     }  // end map loop
0390     return hasLambda;
0391 
0392   }  // end 'HasLambda(JetInfo&)'
0393 
0394 
0395 
0396   // --------------------------------------------------------------------------
0397   //! Check if a particle satisfies cuts
0398   // --------------------------------------------------------------------------
0399   bool SLambdaJetHunter::IsGoodParticle(Types::ParInfo& particle) {
0400 
0401     // print debug statement
0402     if (m_config.isDebugOn && (m_config.verbosity > 5)) {
0403       cout << "SLambdaJetHunter::IsGoodParticle(ParInfo&) checking if particle is good" << endl;
0404     }
0405 
0406     // check charge if needed
0407     bool isGoodCharge = true;
0408     if (m_config.isCharged) {
0409       isGoodCharge = (particle.GetCharge() != 0.);
0410     }
0411 
0412     // apply particle cuts & return overall goodness
0413     const bool isInAccept = particle.IsInAcceptance(m_config.parAccept);
0414     return (isGoodCharge && isInAccept);
0415 
0416   }  // end 'IsGoodParticle(ParInfo&)'
0417 
0418 
0419 
0420   // --------------------------------------------------------------------------
0421   //! Check if a lambda satisfies cuts
0422   // --------------------------------------------------------------------------
0423   bool SLambdaJetHunter::IsGoodLambda(Types::ParInfo& lambda) {
0424 
0425     // print debug statement
0426     if (m_config.isDebugOn && (m_config.verbosity > 5)) {
0427       cout << "SLambdaJetHunter::IsGoodLambda(ParInfo&) checking if lambda is good" << endl;
0428     }
0429 
0430     // make sure lambda is in acceptance
0431     const bool isInAccept = lambda.IsInAcceptance(m_config.parAccept);
0432     return isInAccept;
0433 
0434   }  // end 'IsGoodLambda(ParInfo&)'
0435 
0436 
0437 
0438   // --------------------------------------------------------------------------
0439   //! Check if pid corresponds to a lambda
0440   // --------------------------------------------------------------------------
0441   bool SLambdaJetHunter::IsLambda(const int pid) {
0442 
0443     // print debug statement
0444     if (m_config.isDebugOn && (m_config.verbosity > 5)) {
0445       cout << "SLambdaJetHunter::IsLambda(int) checking if particle is a lambda" << endl;
0446     }
0447 
0448     return (abs(pid) == m_const.pidLambda);
0449 
0450   }  // end 'IsLambda(int)'
0451 
0452 
0453 
0454   // --------------------------------------------------------------------------
0455   //! Check if a lambda hasn't been found yet
0456   // --------------------------------------------------------------------------
0457   bool SLambdaJetHunter::IsNewLambda(const int id) {
0458 
0459     // print debug statement
0460     if (m_config.isDebugOn && (m_config.verbosity > 5)) {
0461       cout << "SLambdaJetHunter::IsNewLambda(int) checking if lambda is already found" << endl;
0462     }
0463 
0464     return (m_mapLambdaJetAssoc.count(id) == 0);
0465 
0466   }  // end 'IsNewLambda(int)'
0467 
0468 
0469 
0470   // --------------------------------------------------------------------------
0471   //! Check if a HepMC particle is in a decay chain
0472   // --------------------------------------------------------------------------
0473   bool SLambdaJetHunter::IsInHepMCDecayChain(const int idToFind, HepMC::GenVertex* vtxStart) {
0474 
0475     // print debug statement
0476     if (m_config.isDebugOn && (m_config.verbosity > 5)) {
0477       cout << "SLambdaJetHunter::IsInHepMCDecayChain(int, HepMC::GenVertex*) checking if particle is in HepMC decay chain" << endl;
0478     }
0479 
0480     // make sure vectors are clear and add initial vertex
0481     m_vecVtxChecking.clear();
0482     m_vecVtxToCheck.clear();
0483     m_vecVtxToCheck.push_back(vtxStart);
0484 
0485     // search of all connected vertices for particle w/ barcode 'idToFind'
0486     int  nVtxChecked     = 0;
0487     bool isParticleFound = false;
0488     while (!isParticleFound && (nVtxChecked < m_const.maxVtxToCheck)) {
0489 
0490       // transfer vertices to-check to being-checked list
0491       m_vecVtxChecking.clear();
0492       for (auto vertex : m_vecVtxToCheck) {
0493         m_vecVtxChecking.push_back(vertex);
0494       }
0495       m_vecVtxToCheck.clear();
0496 
0497       // iterate over vertices being checked
0498       for (auto vertex : m_vecVtxChecking) {
0499 
0500         // iterate over particles in vertex
0501         HepMC::GenVertex::particles_out_const_iterator outPar;
0502         for (
0503           outPar = vertex -> particles_out_const_begin();
0504           outPar != vertex -> particles_out_const_end();
0505           ++outPar
0506         ) {
0507           if ((*outPar) -> barcode() == idToFind) {
0508             isParticleFound = true;
0509             break;
0510           } else {
0511             m_vecVtxToCheck.push_back((*outPar) -> end_vertex());
0512           }
0513         }  // end particle loop
0514 
0515         // if found particler, exit
0516         //   - otherwise increment no. of vertices checked
0517         if (isParticleFound) {
0518           break;
0519         } else {
0520           ++nVtxChecked;
0521         }
0522       }  // end vertex loop
0523     }  // end while (!isParticleFound)
0524     return isParticleFound;
0525 
0526   }  // end 'IsInHepMCDecayChain(int, HepMC::GenVertex*)'
0527 
0528 
0529 
0530   // --------------------------------------------------------------------------
0531   //! Check if a PHG4 particle is in a decay chain
0532   // --------------------------------------------------------------------------
0533   bool SLambdaJetHunter::IsInPHG4DecayChain(const int idToFind, const int idLambda, PHCompositeNode* topNode) {
0534 
0535     // print debug statement
0536     if (m_config.isDebugOn && (m_config.verbosity > 5)) {
0537       cout << "SLambdaJetHunter::IsInPHG4DecayChain(int) checking if particle is in Geant decay chain" << endl;
0538     }
0539 
0540     // grab truth information
0541     PHG4TruthInfoContainer*            container = Interfaces::GetTruthContainer(topNode);
0542     PHG4TruthInfoContainer::ConstRange particles = container -> GetParticleRange();
0543 
0544     // loop over particles
0545     m_vecIDToCheck.clear();
0546     for (
0547       PHG4TruthInfoContainer::ConstIterator itPar = particles.first;
0548       itPar != particles.second;
0549       ++itPar
0550     ) {
0551 
0552       // grab current particle
0553       PHG4Particle* phPar = itPar -> second;
0554 
0555       // grab parent if information is available
0556       PHG4Particle* phParent = NULL;
0557       if (phPar -> get_parent_id() != 0) {
0558         phParent = container -> GetParticle(phPar -> get_parent_id());
0559       } else {
0560         continue;
0561       }
0562 
0563       // check if parent matches lambda barcode/pid
0564       const bool hasLambdaBarcode = (phParent -> get_barcode() == idLambda);
0565       const bool hasLambdaPID     = (phParent -> get_pid()     == m_const.pidLambda);
0566       if (hasLambdaBarcode && hasLambdaPID) {
0567             m_vecIDToCheck.push_back(phPar -> get_barcode());
0568       }
0569     }  // end particle loop
0570 
0571     // check if idToFind is among barcodes
0572     bool isParticleFound = false;
0573     for (const int idToCheck : m_vecIDToCheck) {
0574       if (idToCheck == idToFind) {
0575         isParticleFound = true;
0576         break;
0577       }
0578     }
0579     return isParticleFound;
0580 
0581   }  // end 'IsInPHG4DecayChain(int, PHG4Particle*, PHCompositeNode*)'
0582 
0583 
0584 
0585   // --------------------------------------------------------------------------
0586   //! Calculate z between a lambda and its associated jet
0587   // --------------------------------------------------------------------------
0588   double SLambdaJetHunter::GetLambdaAssocZ(Types::ParInfo& lambda) {
0589 
0590     // print debug statement
0591     if (m_config.isDebugOn && (m_config.verbosity > 5)) {
0592       cout << "SLambdaJetHunter::GetLambdaAssocZ(ParInfo&) getting z of lambda relative to associated jet" << endl;
0593     }
0594 
0595     // get index of associated jet
0596     const int iAssoc = m_mapLambdaJetAssoc[lambda.GetBarcode()];
0597 
0598     // get total momenta
0599     const double pJet    = hypot(m_jetInfo.at(iAssoc).GetPX(), m_jetInfo.at(iAssoc).GetPY(), m_jetInfo.at(iAssoc).GetPZ());
0600     const double pLambda = hypot(lambda.GetPX(), lambda.GetPY(), lambda.GetPZ());
0601 
0602     // calculate z and return
0603     const double zLambda = pLambda / pJet;
0604     return zLambda;
0605 
0606   }  // end 'GetLambdaAssocZ(ParInfo& lambda)'
0607 
0608 
0609 
0610   // --------------------------------------------------------------------------
0611   //! Calculate delta-r between a lambda and the axis of its associated jet
0612   // --------------------------------------------------------------------------
0613   double SLambdaJetHunter::GetLambdaAssocDr(Types::ParInfo& lambda) {
0614 
0615     // print debug statement
0616     if (m_config.isDebugOn && (m_config.verbosity > 5)) {
0617       cout << "SLambdaJetHunter::GetLambdaAssocDr(ParInfo&) getting delta-R of lambda relative to associated jet" << endl;
0618     }
0619 
0620     // get index of associated jet
0621     const int iAssoc = m_mapLambdaJetAssoc[lambda.GetBarcode()];
0622 
0623     // calculate delta-phi
0624     double dfLambda = lambda.GetPhi() - m_jetInfo.at(iAssoc).GetPhi();
0625     if (dfLambda < 0.)             dfLambda += TMath::TwoPi();
0626     if (dfLambda > TMath::TwoPi()) dfLambda -= TMath::TwoPi();
0627 
0628     // calculate dr and return
0629     const double dhLambda = lambda.GetEta() - m_jetInfo.at(iAssoc).GetEta();
0630     const double drLambda = hypot(dfLambda, dhLambda);
0631     return drLambda;
0632 
0633   }  // end 'GetLambdaAssocZ(ParInfo& lambda)'
0634 
0635 
0636 
0637   // --------------------------------------------------------------------------
0638   //! Calculate the total no. of tagged jets
0639   // --------------------------------------------------------------------------
0640   uint64_t SLambdaJetHunter::GetNTaggedJets() {
0641 
0642     // print debug statement
0643     if (m_config.isDebugOn && (m_config.verbosity > 5)) {
0644       cout << "SLambdaJetHunter::GetNTaggedJets() getting no. of jets with a lambda in them" << endl;
0645     }
0646 
0647     uint64_t nTaggedJets = 0;
0648     for (auto& [lamID, jetID] : m_mapLambdaJetAssoc) {
0649       if (jetID > -1) ++nTaggedJets;
0650     }
0651     return nTaggedJets;
0652 
0653   }  // end 'GetNTaggedJets()'
0654 
0655 
0656 
0657   // --------------------------------------------------------------------------
0658   //! Associate lambda to a jet by inspecting constituents' barcodes
0659   // --------------------------------------------------------------------------
0660   optional<int> SLambdaJetHunter::HuntLambdasByBarcode(Types::ParInfo& lambda) {
0661 
0662     // print debug statement
0663     if (m_config.isDebugOn) {
0664       cout << "SLambdaJetHunter::HuntLambdasByBarcode(ParInfo&) hunting for lambdas by inspecting barcode" << endl;
0665     }
0666 
0667     // return nullopt by default
0668     optional<int> iAssocJet = nullopt;
0669 
0670     // loop through jet constituents and check barcode
0671     bool foundAssocJet = false;
0672     for (size_t iJet = 0; iJet < m_jetInfo.size(); iJet++) {
0673       for (size_t iCst = 0; iCst < m_cstInfo[iJet].size(); iCst++) {
0674         if (m_cstInfo[iJet][iCst].GetCstID() == lambda.GetBarcode()) {
0675           foundAssocJet = true;
0676           iAssocJet     = iJet;
0677           break;
0678         }
0679       }  // end cst loop
0680 
0681       // if found association, exit loop
0682       if (foundAssocJet) break;
0683 
0684     }  // end jet loop
0685     return iAssocJet;
0686 
0687   }  // end 'HuntLambdasByBarcode(ParInfo&)'
0688 
0689 
0690 
0691   // --------------------------------------------------------------------------
0692   //! Associate lambda to a jet by inspecting constituents' decay chains
0693   // --------------------------------------------------------------------------
0694   optional<int> SLambdaJetHunter::HuntLambdasByDecayChain(Types::ParInfo& lambda, PHCompositeNode* topNode) {
0695 
0696     // print debug statement
0697     if (m_config.isDebugOn) {
0698       cout << "SLambdaJetHunter::HuntLambdasByDecayChain(ParInfo&, PHCompositeNode*) hunting for lambdas by inspecting decay chains" << endl;
0699     }
0700 
0701     // return nullopt by default
0702     optional<int> iAssocJet = nullopt;
0703 
0704     // grab PHG4 and HepMC information
0705     HepMC::GenParticle* hepLambda = Tools::GetHepMCGenParticleFromBarcode(lambda.GetBarcode(), topNode);
0706     HepMC::GenVertex*   hepEndVtx = hepLambda -> end_vertex();
0707 
0708     // loop over constituents and check which cst.s the lambda decayed into
0709     bool foundAssocJet = false;
0710     for (size_t iJet = 0; iJet < m_cstInfo.size(); iJet++) {
0711       for (size_t iCst = 0; iCst < m_cstInfo[iJet].size(); iCst++) {
0712 
0713         // check if constituent is in decay chain,
0714         bool isInDecayChain = false;
0715         if (hepEndVtx) {
0716           isInDecayChain = IsInHepMCDecayChain(m_cstInfo[iJet][iCst].GetCstID(), hepEndVtx);
0717         } else {
0718           isInDecayChain = IsInPHG4DecayChain(m_cstInfo[iJet][iCst].GetCstID(), lambda.GetBarcode(), topNode);
0719         }
0720 
0721         // if so, add lambda to lists
0722         if (isInDecayChain) {
0723           foundAssocJet = true;
0724           iAssocJet     = iJet;
0725           break;
0726         }
0727       } // end cst loop
0728 
0729       // if found association, exit loop
0730       if (foundAssocJet)  break;
0731 
0732     }  // end jet loop
0733     return iAssocJet;
0734 
0735   }  // end 'HuntLambdasByDecayChain(PHCompositeNode*)'
0736 
0737 
0738 
0739   // --------------------------------------------------------------------------
0740   //! Associate lambda to a jet based on distance
0741   // --------------------------------------------------------------------------
0742   optional<int> SLambdaJetHunter::HuntLambdasByDistance(Types::ParInfo& lambda) {
0743 
0744     // print debug statement
0745     if (m_config.isDebugOn) {
0746       cout << "SLambdaJetHunter::HuntLambdasByDistance(ParInfo&) hunting for lambdas by distance" << endl;
0747     }
0748 
0749     // return nullopt by default
0750     optional<int> iAssocJet = nullopt;
0751 
0752     // loop over jets
0753     double drAssocBest = numeric_limits<double>::max();
0754     for (size_t iJet = 0; iJet < m_jetInfo.size(); iJet++) {
0755 
0756       // calculate delta-phi
0757       double dfAssoc = lambda.GetPhi() - m_jetInfo[iJet].GetPhi();
0758       if (dfAssoc < 0.)             dfAssoc += TMath::TwoPi();
0759       if (dfAssoc > TMath::TwoPi()) dfAssoc -= TMath::TwoPi();
0760 
0761       // calculate dr
0762       const double dhAssoc = lambda.GetEta() - m_jetInfo[iJet].GetEta();
0763       const double drAssoc = hypot(dfAssoc, dhAssoc);
0764 
0765       // check if lambda is within rJet
0766       const bool isDistGood = (drAssoc < m_config.rJet);
0767       const bool isDistBest = (drAssoc < drAssocBest);
0768       if (isDistGood && isDistBest) {
0769         drAssocBest = drAssoc;
0770         iAssocJet   = iJet;
0771       }
0772     }  // end jet loop
0773     return iAssocJet;
0774 
0775   }  // end 'HuntLambdasByDistance(ParInfo& lambda)'
0776 
0777 }  // end SColdQcdCorrelatorAnalysis namespace
0778 
0779 // end ------------------------------------------------------------------------
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!