Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TracksInJets.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TracksInJets.cc
1 
2 #include "TracksInJets.h"
3 
4 #include <fun4all/Fun4AllReturnCodes.h>
5 #include <Acts/Definitions/Algebra.hpp>
6 
7 #include <phool/PHCompositeNode.h>
8 #include <phool/getClass.h>
9 
10 #include <g4main/PHG4Hit.h>
11 #include <g4main/PHG4Particle.h>
12 #include <g4main/PHG4TruthInfoContainer.h>
13 #include <g4main/PHG4VtxPoint.h>
14 
15 #include <g4jets/Jet.h>
16 #include <g4jets/JetMap.h>
17 
18 #include <trackbase_historic/SvtxTrack.h>
19 #include <trackbase_historic/SvtxTrackMap.h>
20 #include <trackbase_historic/SvtxVertex.h>
21 #include <trackbase_historic/SvtxVertexMap.h>
22 
23 #include <particleflowreco/ParticleFlowElementContainer.h>
24 #include <particleflowreco/ParticleFlowElement.h>
25 
26 #include <g4eval/JetEvalStack.h>
27 #include <g4eval/SvtxEvalStack.h>
28 
29 #include <HepMC/GenEvent.h>
30 #include <HepMC/GenVertex.h>
31 #include <phhepmc/PHHepMCGenEvent.h>
32 #include <phhepmc/PHHepMCGenEventMap.h>
33 
34 
35 //____________________________________________________________________________..
36 TracksInJets::TracksInJets(const std::string& name,
37  const std::string& out):
38  SubsysReco(name)
39  , m_outfilename(out)
40 {
41 }
42 
43 //____________________________________________________________________________..
44 TracksInJets::~TracksInJets()
45 {
46 }
47 
48 //____________________________________________________________________________..
49 int TracksInJets::Init(PHCompositeNode *topNode)
50 {
51  return Fun4AllReturnCodes::EVENT_OK;
52 }
53 
54 //____________________________________________________________________________..
55 int TracksInJets::InitRun(PHCompositeNode *topNode)
56 {
57  m_outfile = new TFile(m_outfilename.c_str(), "RECREATE");
58  m_tree = new TTree("tree","A tree with tracks in jets");
59  setBranches();
60 
61  return Fun4AllReturnCodes::EVENT_OK;
62 }
63 
64 //____________________________________________________________________________..
65 int TracksInJets::process_event(PHCompositeNode *topNode)
66 {
67  PHHepMCGenEventMap *geneventmap = findNode::getClass<PHHepMCGenEventMap>(
68  topNode, "PHHepMCGenEventMap");
69  if (!geneventmap)
70  {
71  std::cout << PHWHERE
72  << " - Fatal error - missing node PHHepMCGenEventMap"
73  << std::endl;
74  return Fun4AllReturnCodes::ABORTEVENT;
75  }
76 
77  PHHepMCGenEvent *genevt = geneventmap->get(m_embeddingid);
78  if (!genevt)
79  {
80  std::cout << PHWHERE
81  << " - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID "
82  << m_embeddingid;
83  std::cout << ". Print PHHepMCGenEventMap:";
84  geneventmap->identify();
85  return Fun4AllReturnCodes::ABORTEVENT;
86  }
87 
88  JetEvalStack* jetevalstack = new JetEvalStack(topNode, m_recojetmapname, m_truthjetmapname);
89  SvtxEvalStack* svtxevalstack = new SvtxEvalStack(topNode);
90  if(!jetevalstack or !svtxevalstack)
91  {
92  std::cout << "No eval stacks, can't continue." << std::endl;
93  return Fun4AllReturnCodes::ABORTEVENT;
94  }
95 
96  svtxevalstack->next_event(topNode);
97  //SvtxTrackEval* trackeval = svtxevalstack->get_track_eval();
98 
99  JetMap *truth_jets = findNode::getClass<JetMap>(topNode, m_truthjetmapname);
100  JetMap *reco_jets = findNode::getClass<JetMap>(topNode, m_recojetmapname);
101 
102  JetTruthEval* jettrutheval = jetevalstack->get_truth_eval();
103 
104  ParticleFlowElementContainer* pflowcontainer = findNode::getClass<ParticleFlowElementContainer>(topNode, "ParticleFlowElements");
105  if(!pflowcontainer)
106  {
107  std::cout << "No particle flow elements on node tree, can't continue."
108  << std::endl;
109  return Fun4AllReturnCodes::ABORTEVENT;
110  }
111 
112  //HepMC::GenEvent *theEvent = genevt->getEvent();
113 
114 
115  PHG4TruthInfoContainer *truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
116 
117  PHG4VtxPoint *first_point = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
118 
119  float truth_vx = first_point->get_x();
120  float truth_vy = first_point->get_y();
121  float truth_vz = first_point->get_z();
122 
123 
124  for (JetMap::Iter iter = truth_jets->begin(); iter != truth_jets->end();
125  ++iter)
126  {
127  Jet *truthjet = iter->second;
128  if(truthjet->get_pt() < m_minjettruthpt)
129  { continue; }
130  if(fabs(truthjet->get_eta()) > 2)
131  { continue; }
132 
133  m_truthjetpx = truthjet->get_px();
134  m_truthjetpy = truthjet->get_py();
135  m_truthjetpz = truthjet->get_pz();
136  m_truthjete = truthjet->get_e();
137 
138  std::set<PHG4Particle*> truthjetconst =
139  jettrutheval->all_truth_particles(truthjet);
140  m_truthjetconst = 0;
141  for(const PHG4Particle* truthpart : truthjetconst)
142  {
143  int truthpid = truthpart->get_pid();
144  int fabstruthpid = fabs(truthpid);
145  if(! (fabstruthpid == 211 or fabstruthpid == 321 or fabstruthpid == 2212 or fabstruthpid == 11 or fabstruthpid == 13))
146  { continue; }
147  m_truthjettrackpx.push_back(truthpart->get_px());
148  m_truthjettrackpy.push_back(truthpart->get_py());
149  m_truthjettrackpz.push_back(truthpart->get_pz());
150  m_truthjettrackvx.push_back(truth_vx);
151  m_truthjettrackvy.push_back(truth_vy);
152  m_truthjettrackvz.push_back(truth_vz);
153 
154  m_truthjetconst++;
155  }
156 
157  Jet* matchedrecojet = nullptr;
158  float mindr = std::numeric_limits<float>::max();
159  for (JetMap::Iter riter = reco_jets->begin(); riter != reco_jets->end();
160  ++riter)
161  {
162  Jet* mjet = riter->second;
163  if(mjet->get_pt() < 5)
164  { continue; }
165  float dr = dR(mjet->get_eta(), truthjet->get_eta(),
166  mjet->get_phi(), truthjet->get_phi());
167  if(dr < mindr)
168  {
169  mindr = dr;
170  matchedrecojet = mjet;
171  }
172  }
173 
174  if(!matchedrecojet)
175  {
177  m_recojetpx = std::numeric_limits<float>::max();
178  m_recojetpy = std::numeric_limits<float>::max();
179  m_recojetpz = std::numeric_limits<float>::max();
180  m_recojete = std::numeric_limits<float>::max();
181  m_recojetconst = 0;
182  m_tree->Fill();
183  continue;
184  }
185 
187  m_recojetpx = matchedrecojet->get_px();
188  m_recojetpy = matchedrecojet->get_py();
189  m_recojetpz = matchedrecojet->get_pz();
190  m_recojete = matchedrecojet->get_e();
191 
192  m_recojetconst = matchedrecojet->size_comp();
193  for(auto citer = matchedrecojet->begin_comp(); citer != matchedrecojet->end_comp(); ++citer)
194  {
195  ParticleFlowElement *pflow = pflowcontainer->getParticleFlowElement(citer->second);
196  ParticleFlowElement::PFLOWTYPE type = pflow->get_type();
197 
198  if(!(type == ParticleFlowElement::PFLOWTYPE::MATCHED_CHARGED_HADRON or
199  type == ParticleFlowElement::PFLOWTYPE::UNMATCHED_CHARGED_HADRON))
200  { continue; }
201 
202  m_recojettrackpx.push_back(pflow->get_px());
203  m_recojettrackpy.push_back(pflow->get_py());
204  m_recojettrackpz.push_back(pflow->get_pz());
205  m_recojettracke.push_back(pflow->get_e());
206  m_recojettracktype.push_back(pflow->get_type());
207  }
208 
209  m_tree->Fill();
210 
211  }
212 
213  return Fun4AllReturnCodes::EVENT_OK;
214 }
215 
216 //____________________________________________________________________________..
217 int TracksInJets::ResetEvent(PHCompositeNode *topNode)
218 {
219  m_truthjettrackpx.clear();
220  m_truthjettrackpy.clear();
221  m_truthjettrackpz.clear();
222  m_truthjettrackpid.clear();
223  m_truthjettrackvx.clear();
224  m_truthjettrackvy.clear();
225  m_truthjettrackvz.clear();
226  m_recojettrackpx.clear();
227  m_recojettrackpy.clear();
228  m_recojettrackpz.clear();
229  m_recojettracke.clear();
230  m_recojettrackvx.clear();
231  m_recojettrackvy.clear();
232  m_recojettrackvz.clear();
233  m_recojettracktype.clear();
234  m_recojettrackpcax.clear();
235  m_recojettrackpcay.clear();
236  m_recojettrackpcaz.clear();
237  m_recojettracktype.clear();
238 
239  return Fun4AllReturnCodes::EVENT_OK;
240 }
241 
242 //____________________________________________________________________________..
243 int TracksInJets::End(PHCompositeNode *topNode)
244 {
245  m_outfile->Write();
246  m_outfile->Close();
247  return Fun4AllReturnCodes::EVENT_OK;
248 }
249 
250 void TracksInJets::setBranches()
251 {
252  m_tree->Branch("truthjetpx", &m_truthjetpx, "truthjetpx/F");
253  m_tree->Branch("truthjetpy", &m_truthjetpy, "truthjetpy/F");
254  m_tree->Branch("truthjetpz", &m_truthjetpz, "truthjetpz/F");
255  m_tree->Branch("truthjete", &m_truthjete, "truthjete/F");
256  m_tree->Branch("recojetpx", &m_recojetpx, "recojetpx/F");
257  m_tree->Branch("recojetpy", &m_recojetpy, "recojetpy/F");
258  m_tree->Branch("recojetpz", &m_recojetpz, "recojetpz/F");
259  m_tree->Branch("recojete", &m_recojete, "recojete/F");
260  m_tree->Branch("truthjetconst", &m_truthjetconst, "truthjetconst/I");
261  m_tree->Branch("recojetconst", &m_recojetconst, "recojetconst/I");
262 
263  m_tree->Branch("truthjettrackpx",&m_truthjettrackpx);
264  m_tree->Branch("truthjettrackpy", &m_truthjettrackpy);
265  m_tree->Branch("truthjettrackpz", &m_truthjettrackpz);
266  m_tree->Branch("truthjettrackpid", &m_truthjettrackpid);
267  m_tree->Branch("truthjettrackvx", &m_truthjettrackvx);
268  m_tree->Branch("truthjettrackvy", &m_truthjettrackvy);
269  m_tree->Branch("truthjettrackvz", &m_truthjettrackvz);
270 
271  m_tree->Branch("recojettrackpx",&m_recojettrackpx);
272  m_tree->Branch("recojettrackpy", &m_recojettrackpy);
273  m_tree->Branch("recojettrackpz", &m_recojettrackpz);
274  m_tree->Branch("recojettracke", &m_recojettracke);
275  m_tree->Branch("recojettrackvx", &m_recojettrackvx);
276  m_tree->Branch("recojettrackvy", &m_recojettrackvy);
277  m_tree->Branch("recojettrackvz", &m_recojettrackvz);
278  m_tree->Branch("recojettrackpcax", &m_recojettrackpcax);
279  m_tree->Branch("recojettrackpcay", &m_recojettrackpcay);
280  m_tree->Branch("recojettrackpcaz", &m_recojettrackpcaz);
281  m_tree->Branch("recojettracktype", &m_recojettracktype);
282 
283 
284 }
285 
286 
287 float TracksInJets::dR(const float& eta1, const float& eta2,
288  const float& phi1, const float& phi2)
289 {
290  float deta = eta1-eta2;
291  float dphi = phi1-phi2;
292  if(dphi > M_PI)
293  { dphi -= 2. * M_PI; }
294  if(dphi < -M_PI)
295  { dphi += 2. * M_PI; }
296 
297  return sqrt(pow(deta,2) + pow(dphi,2));
298 
299 }