File indexing completed on 2025-08-05 08:12:46
0001
0002 #include "TracksInJets.h"
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <Acts/Definitions/Algebra.hpp>
0006
0007 #include <phool/PHCompositeNode.h>
0008 #include <phool/getClass.h>
0009
0010 #include <g4main/PHG4Hit.h>
0011 #include <g4main/PHG4Particle.h>
0012 #include <g4main/PHG4TruthInfoContainer.h>
0013 #include <g4main/PHG4VtxPoint.h>
0014
0015 #include <g4jets/Jet.h>
0016 #include <g4jets/JetMap.h>
0017
0018 #include <trackbase_historic/SvtxTrack.h>
0019 #include <trackbase_historic/SvtxTrackMap.h>
0020 #include <trackbase_historic/SvtxVertex.h>
0021 #include <trackbase_historic/SvtxVertexMap.h>
0022
0023 #include <particleflowreco/ParticleFlowElementContainer.h>
0024 #include <particleflowreco/ParticleFlowElement.h>
0025
0026 #include <g4eval/JetEvalStack.h>
0027 #include <g4eval/SvtxEvalStack.h>
0028
0029 #include <HepMC/GenEvent.h>
0030 #include <HepMC/GenVertex.h>
0031 #include <phhepmc/PHHepMCGenEvent.h>
0032 #include <phhepmc/PHHepMCGenEventMap.h>
0033
0034
0035
0036 TracksInJets::TracksInJets(const std::string& name,
0037 const std::string& out):
0038 SubsysReco(name)
0039 , m_outfilename(out)
0040 {
0041 }
0042
0043
0044 TracksInJets::~TracksInJets()
0045 {
0046 }
0047
0048
0049 int TracksInJets::Init(PHCompositeNode *topNode)
0050 {
0051 return Fun4AllReturnCodes::EVENT_OK;
0052 }
0053
0054
0055 int TracksInJets::InitRun(PHCompositeNode *topNode)
0056 {
0057 m_outfile = new TFile(m_outfilename.c_str(), "RECREATE");
0058 m_tree = new TTree("tree","A tree with tracks in jets");
0059 setBranches();
0060
0061 return Fun4AllReturnCodes::EVENT_OK;
0062 }
0063
0064
0065 int TracksInJets::process_event(PHCompositeNode *topNode)
0066 {
0067 PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(
0068 topNode, "PHHepMCGenEventMap");
0069 if (!geneventmap)
0070 {
0071 std::cout << PHWHERE
0072 << " - Fatal error - missing node PHHepMCGenEventMap"
0073 << std::endl;
0074 return Fun4AllReturnCodes::ABORTEVENT;
0075 }
0076
0077 PHHepMCGenEvent *genevt = geneventmap->get(m_embeddingid);
0078 if (!genevt)
0079 {
0080 std::cout << PHWHERE
0081 << " - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID "
0082 << m_embeddingid;
0083 std::cout << ". Print PHHepMCGenEventMap:";
0084 geneventmap->identify();
0085 return Fun4AllReturnCodes::ABORTEVENT;
0086 }
0087
0088 JetEvalStack* jetevalstack = new JetEvalStack(topNode, m_recojetmapname, m_truthjetmapname);
0089 SvtxEvalStack* svtxevalstack = new SvtxEvalStack(topNode);
0090 if(!jetevalstack or !svtxevalstack)
0091 {
0092 std::cout << "No eval stacks, can't continue." << std::endl;
0093 return Fun4AllReturnCodes::ABORTEVENT;
0094 }
0095
0096 svtxevalstack->next_event(topNode);
0097
0098
0099 JetMap *truth_jets = findNode::getClass<JetMap>(topNode, m_truthjetmapname);
0100 JetMap *reco_jets = findNode::getClass<JetMap>(topNode, m_recojetmapname);
0101
0102 JetTruthEval* jettrutheval = jetevalstack->get_truth_eval();
0103
0104 ParticleFlowElementContainer* pflowcontainer = findNode::getClass<ParticleFlowElementContainer>(topNode, "ParticleFlowElements");
0105 if(!pflowcontainer)
0106 {
0107 std::cout << "No particle flow elements on node tree, can't continue."
0108 << std::endl;
0109 return Fun4AllReturnCodes::ABORTEVENT;
0110 }
0111
0112
0113
0114
0115 PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0116
0117 PHG4VtxPoint *first_point = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
0118
0119 float truth_vx = first_point->get_x();
0120 float truth_vy = first_point->get_y();
0121 float truth_vz = first_point->get_z();
0122
0123
0124 for (JetMap::Iter iter = truth_jets->begin(); iter != truth_jets->end();
0125 ++iter)
0126 {
0127 Jet *truthjet = iter->second;
0128 if(truthjet->get_pt() < m_minjettruthpt)
0129 { continue; }
0130 if(fabs(truthjet->get_eta()) > 2)
0131 { continue; }
0132
0133 m_truthjetpx = truthjet->get_px();
0134 m_truthjetpy = truthjet->get_py();
0135 m_truthjetpz = truthjet->get_pz();
0136 m_truthjete = truthjet->get_e();
0137
0138 std::set<PHG4Particle*> truthjetconst =
0139 jettrutheval->all_truth_particles(truthjet);
0140 m_truthjetconst = 0;
0141 for(const PHG4Particle* truthpart : truthjetconst)
0142 {
0143 int truthpid = truthpart->get_pid();
0144 int fabstruthpid = fabs(truthpid);
0145 if(! (fabstruthpid == 211 or fabstruthpid == 321 or fabstruthpid == 2212 or fabstruthpid == 11 or fabstruthpid == 13))
0146 { continue; }
0147 m_truthjettrackpx.push_back(truthpart->get_px());
0148 m_truthjettrackpy.push_back(truthpart->get_py());
0149 m_truthjettrackpz.push_back(truthpart->get_pz());
0150 m_truthjettrackvx.push_back(truth_vx);
0151 m_truthjettrackvy.push_back(truth_vy);
0152 m_truthjettrackvz.push_back(truth_vz);
0153
0154 m_truthjetconst++;
0155 }
0156
0157 Jet* matchedrecojet = nullptr;
0158 float mindr = std::numeric_limits<float>::max();
0159 for (JetMap::Iter riter = reco_jets->begin(); riter != reco_jets->end();
0160 ++riter)
0161 {
0162 Jet* mjet = riter->second;
0163 if(mjet->get_pt() < 5)
0164 { continue; }
0165 float dr = dR(mjet->get_eta(), truthjet->get_eta(),
0166 mjet->get_phi(), truthjet->get_phi());
0167 if(dr < mindr)
0168 {
0169 mindr = dr;
0170 matchedrecojet = mjet;
0171 }
0172 }
0173
0174 if(!matchedrecojet)
0175 {
0176
0177 m_recojetpx = std::numeric_limits<float>::max();
0178 m_recojetpy = std::numeric_limits<float>::max();
0179 m_recojetpz = std::numeric_limits<float>::max();
0180 m_recojete = std::numeric_limits<float>::max();
0181 m_recojetconst = 0;
0182 m_tree->Fill();
0183 continue;
0184 }
0185
0186
0187 m_recojetpx = matchedrecojet->get_px();
0188 m_recojetpy = matchedrecojet->get_py();
0189 m_recojetpz = matchedrecojet->get_pz();
0190 m_recojete = matchedrecojet->get_e();
0191
0192 m_recojetconst = matchedrecojet->size_comp();
0193 for(auto citer = matchedrecojet->begin_comp(); citer != matchedrecojet->end_comp(); ++citer)
0194 {
0195 ParticleFlowElement *pflow = pflowcontainer->getParticleFlowElement(citer->second);
0196 ParticleFlowElement::PFLOWTYPE type = pflow->get_type();
0197
0198 if(!(type == ParticleFlowElement::PFLOWTYPE::MATCHED_CHARGED_HADRON or
0199 type == ParticleFlowElement::PFLOWTYPE::UNMATCHED_CHARGED_HADRON))
0200 { continue; }
0201
0202 m_recojettrackpx.push_back(pflow->get_px());
0203 m_recojettrackpy.push_back(pflow->get_py());
0204 m_recojettrackpz.push_back(pflow->get_pz());
0205 m_recojettracke.push_back(pflow->get_e());
0206 m_recojettracktype.push_back(pflow->get_type());
0207 }
0208
0209 m_tree->Fill();
0210
0211 }
0212
0213 return Fun4AllReturnCodes::EVENT_OK;
0214 }
0215
0216
0217 int TracksInJets::ResetEvent(PHCompositeNode *topNode)
0218 {
0219 m_truthjettrackpx.clear();
0220 m_truthjettrackpy.clear();
0221 m_truthjettrackpz.clear();
0222 m_truthjettrackpid.clear();
0223 m_truthjettrackvx.clear();
0224 m_truthjettrackvy.clear();
0225 m_truthjettrackvz.clear();
0226 m_recojettrackpx.clear();
0227 m_recojettrackpy.clear();
0228 m_recojettrackpz.clear();
0229 m_recojettracke.clear();
0230 m_recojettrackvx.clear();
0231 m_recojettrackvy.clear();
0232 m_recojettrackvz.clear();
0233 m_recojettracktype.clear();
0234 m_recojettrackpcax.clear();
0235 m_recojettrackpcay.clear();
0236 m_recojettrackpcaz.clear();
0237 m_recojettracktype.clear();
0238
0239 return Fun4AllReturnCodes::EVENT_OK;
0240 }
0241
0242
0243 int TracksInJets::End(PHCompositeNode *topNode)
0244 {
0245 m_outfile->Write();
0246 m_outfile->Close();
0247 return Fun4AllReturnCodes::EVENT_OK;
0248 }
0249
0250 void TracksInJets::setBranches()
0251 {
0252 m_tree->Branch("truthjetpx", &m_truthjetpx, "truthjetpx/F");
0253 m_tree->Branch("truthjetpy", &m_truthjetpy, "truthjetpy/F");
0254 m_tree->Branch("truthjetpz", &m_truthjetpz, "truthjetpz/F");
0255 m_tree->Branch("truthjete", &m_truthjete, "truthjete/F");
0256 m_tree->Branch("recojetpx", &m_recojetpx, "recojetpx/F");
0257 m_tree->Branch("recojetpy", &m_recojetpy, "recojetpy/F");
0258 m_tree->Branch("recojetpz", &m_recojetpz, "recojetpz/F");
0259 m_tree->Branch("recojete", &m_recojete, "recojete/F");
0260 m_tree->Branch("truthjetconst", &m_truthjetconst, "truthjetconst/I");
0261 m_tree->Branch("recojetconst", &m_recojetconst, "recojetconst/I");
0262
0263 m_tree->Branch("truthjettrackpx",&m_truthjettrackpx);
0264 m_tree->Branch("truthjettrackpy", &m_truthjettrackpy);
0265 m_tree->Branch("truthjettrackpz", &m_truthjettrackpz);
0266 m_tree->Branch("truthjettrackpid", &m_truthjettrackpid);
0267 m_tree->Branch("truthjettrackvx", &m_truthjettrackvx);
0268 m_tree->Branch("truthjettrackvy", &m_truthjettrackvy);
0269 m_tree->Branch("truthjettrackvz", &m_truthjettrackvz);
0270
0271 m_tree->Branch("recojettrackpx",&m_recojettrackpx);
0272 m_tree->Branch("recojettrackpy", &m_recojettrackpy);
0273 m_tree->Branch("recojettrackpz", &m_recojettrackpz);
0274 m_tree->Branch("recojettracke", &m_recojettracke);
0275 m_tree->Branch("recojettrackvx", &m_recojettrackvx);
0276 m_tree->Branch("recojettrackvy", &m_recojettrackvy);
0277 m_tree->Branch("recojettrackvz", &m_recojettrackvz);
0278 m_tree->Branch("recojettrackpcax", &m_recojettrackpcax);
0279 m_tree->Branch("recojettrackpcay", &m_recojettrackpcay);
0280 m_tree->Branch("recojettrackpcaz", &m_recojettrackpcaz);
0281 m_tree->Branch("recojettracktype", &m_recojettracktype);
0282
0283
0284 }
0285
0286
0287 float TracksInJets::dR(const float& eta1, const float& eta2,
0288 const float& phi1, const float& phi2)
0289 {
0290 float deta = eta1-eta2;
0291 float dphi = phi1-phi2;
0292 if(dphi > M_PI)
0293 { dphi -= 2. * M_PI; }
0294 if(dphi < -M_PI)
0295 { dphi += 2. * M_PI; }
0296
0297 return sqrt(pow(deta,2) + pow(dphi,2));
0298
0299 }
