Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TruthJetTagging.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TruthJetTagging.cc
1 //____________________________________________________________________________..
2 
3 #include "TruthJetTagging.h"
4 
5 #include <fun4all/Fun4AllReturnCodes.h>
6 #include <fun4all/Fun4AllServer.h>
7 #include <fun4all/PHTFileServer.h>
8 #include <phool/PHCompositeNode.h>
9 #include <phool/PHIODataNode.h> // for PHIODataNode
10 #include <phool/PHNode.h> // for PHNode
11 #include <phool/PHNodeIterator.h> // for PHNodeIterator
12 #include <phool/PHObject.h> // for PHObject
13 #include <phool/getClass.h>
14 #include <phool/phool.h> // for PHWHERE
15 #include <g4main/PHG4Utils.h>
16 
17 #pragma GCC diagnostic push
18 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
19 #include <HepMC/GenEvent.h>
20 #include <HepMC/GenVertex.h>
21 #pragma GCC diagnostic pop
22 
23 #include <phhepmc/PHHepMCGenEvent.h>
24 #include <phhepmc/PHHepMCGenEventMap.h>
25 #include <TDatabasePDG.h>
26 
27 #include <g4main/PHG4Particle.h>
28 #include <g4main/PHG4TruthInfoContainer.h>
29 #include <TMath.h>
30 
31 #include <g4jets/JetMap.h>
32 #include <g4jets/Jet.h>
33 #include <g4jets/Jetv1.h>
34 #include <sstream>
35 
36 //____________________________________________________________________________..
37 TruthJetTagging::TruthJetTagging(const std::string &name):
38  SubsysReco(name)
39  , _algorithms()
40  , _radii()
41  , _do_photon_tagging()
42  , _do_hf_tagging()
43  , _embedding_id(1)
44 {
45 
46 }
47 
48 //____________________________________________________________________________..
49 TruthJetTagging::~TruthJetTagging()
50 {
51 
52 }
53 
54 //____________________________________________________________________________..
55 int TruthJetTagging::Init(PHCompositeNode *topNode)
56 {
57  return Fun4AllReturnCodes::EVENT_OK;
58 }
59 
60 //____________________________________________________________________________..
61 int TruthJetTagging::InitRun(PHCompositeNode *topNode)
62 {
63  return Fun4AllReturnCodes::EVENT_OK;
64 }
65 
66 //____________________________________________________________________________..
67 int TruthJetTagging::process_event(PHCompositeNode *topNode)
68 {
69 PHG4TruthInfoContainer* truthcontainer = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
70 if (!truthcontainer)
71  {
72  std::cout
73  << "MyJetAnalysis::process_event - Error can not find DST truthcontainer node "
74  << std::endl;
75 exit(-1);
76 }
77  PHHepMCGenEventMap* geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
78  if (!geneventmap)
79  {
80  std::cout << PHWHERE << " - Fatal error - missing node PHHepMCGenEventMap" << std::endl;
81  return Fun4AllReturnCodes::ABORTRUN;
82  }
83 
84  PHHepMCGenEvent* genevt = geneventmap->get(_embedding_id);
85  if (!genevt)
86  {
87  std::cout << PHWHERE << " - Fatal error - node PHHepMCGenEventMap missing subevent with embedding ID " << _embedding_id;
88  std::cout << ". Print PHHepMCGenEventMap:";
89  geneventmap->identify();
90  return Fun4AllReturnCodes::ABORTRUN;
91  }
92 
93 
94 
95 HepMC::GenEvent* theEvent = genevt->getEvent();
96 
97  int n_algo = _algorithms.size();
98  int n_radii = _radii.size();
99  if (_do_hf_tagging)
100  {
101  if (n_radii != n_algo)
102  {
103  std::cout << "TruthJetTagging::process_event - errorr unequal number of jet radii and algoirthms specified" << std::endl;
104  exit(-1);
105  }
106  }
107  for (int algoiter = 0;algoiter < n_algo;algoiter++)
108  {
109  JetMap* tjets = findNode::getClass<JetMap>(topNode, _algorithms.at(algoiter).c_str());
110  if (!tjets)
111  {
112  std::cout
113  << "MyJetAnalysis::process_event - Error can not find DST JetMap node "
114  << std::endl;
115  exit(-1);
116  }
117  for (JetMap::Iter iter = tjets->begin(); iter != tjets->end(); ++iter)
118  {
119  Jet* tjet = iter->second;
120  assert(tjet);
121 
122  if (_do_photon_tagging)
123  {
124  float photon_fraction = TruthJetTagging::TruthPhotonTagging ( truthcontainer , tjet );
125  tjet->set_property(Jet::prop_gamma,photon_fraction);
126  }
127  if (_do_hf_tagging)
128  {
129 
130  float jet_radius = 0.4;
131  int jet_flavor = TruthJetTagging::hadron_tagging(tjet, theEvent, jet_radius);
132  tjet->set_property(Jet::prop_JetHadronFlavor,jet_flavor);
133  }
134 
135  }
136  }
137  return Fun4AllReturnCodes::EVENT_OK;
138 }
139 
140  float TruthJetTagging::TruthPhotonTagging ( PHG4TruthInfoContainer* truthnode, Jet* tjet)
141  {
142  assert(truthnode);
143  assert(tjet);
144  float jetpt = tjet->get_pt();
145  float max_gamma_pt = 0;
146  for (Jet::ConstIter comp = tjet->begin_comp(); comp != tjet->end_comp(); ++comp)
147  {
148  PHG4Particle* particle = truthnode->GetParticle((*comp).second);
149  if (particle->get_pid() == 22)
150  {
151  float particle_pt = TMath::Sqrt(TMath::Power(particle->get_px(),2) + TMath::Power(particle->get_py(),2));
152  if (particle_pt > max_gamma_pt)
153  {
154  max_gamma_pt = particle_pt;
155  }
156  }
157  }
158  float ratio = max_gamma_pt/jetpt;
159  return ratio;
160  }
161 
162 
163 
164 
165 
166 
167 int TruthJetTagging::hadron_tagging(Jet* this_jet, HepMC::GenEvent* theEvent,
168  const double match_radius)
169 {
170  float this_pt = this_jet->get_pt();
171  float this_phi = this_jet->get_phi();
172  float this_eta = this_jet->get_eta();
173 
174  int jet_flavor = 0;
175  double jet_parton_zt = 0;
176 
177  for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin();
178  p != theEvent->particles_end(); ++p)
179  {
180  float dR = deltaR((*p)->momentum().pseudoRapidity(), this_eta,
181  (*p)->momentum().phi(), this_phi);
182  if (dR > match_radius)
183  continue;
184 
185  int pidabs = 0;
186  TParticlePDG* pdg_p = TDatabasePDG::Instance()->GetParticle((*p)->pdg_id());
187  if (pdg_p)
188  {
189  if (TString(pdg_p->ParticleClass()).BeginsWith("B-"))
190  {
191  pidabs = TDatabasePDG::Instance()->GetParticle("b")->PdgCode();
192  }
193  else if (TString(pdg_p->ParticleClass()).BeginsWith("Charmed"))
194  {
195  pidabs = TDatabasePDG::Instance()->GetParticle("c")->PdgCode();
196  }
197  }
198 
199  const double zt = (*p)->momentum().perp() / this_pt;
200 
201  if (pidabs == TDatabasePDG::Instance()->GetParticle("c")->PdgCode() //
202  or pidabs == TDatabasePDG::Instance()->GetParticle("b")->PdgCode()) // handle heavy quarks only. All other favor tagged as default 0
203  {
204  if (pidabs > abs(jet_flavor)) // heavy quark found
205  {
206  jet_parton_zt = zt;
207  jet_flavor = pidabs;
208  }
209  else if (pidabs == abs(jet_flavor)) // same quark mass. next compare zt
210  {
211  if (zt > jet_parton_zt)
212  {
213  jet_parton_zt = zt;
214  jet_flavor = pidabs;
215  }
216  }
217 
218  // if (pidabs == TDatabasePDG::Instance()->GetParticle("b")->PdgCode())
219  // {
220  // // if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
221  // // {
222  // // std::cout << __PRETTY_FUNCTION__
223  // // << " --BOTTOM--> pt / eta / phi = "
224  // // << (*p)->momentum().perp() << " / "
225  // // << (*p)->momentum().pseudoRapidity() << " / "
226  // // << (*p)->momentum().phi() << " with hadron ";
227  // // pdg_p->Print();
228  // // }
229  // }
230  // else if (pidabs == TDatabasePDG::Instance()->GetParticle("c")->PdgCode())
231  // {
232  // // if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
233  // // {
234  // // std::cout << __PRETTY_FUNCTION__
235  // // << " --CHARM --> pt / eta / phi = "
236  // // << (*p)->momentum().perp() << " / "
237  // // << (*p)->momentum().pseudoRapidity() << " / "
238  // // << (*p)->momentum().phi() << " with hadron ";
239  // // pdg_p->Print();
240  // // }
241  // }
242  }
243  } // for (HepMC::GenEvent::particle_const_iterator p =
244 
245 
246  this_jet->set_property(Jet::prop_JetHadronZT,
247  jet_parton_zt);
248  // this_jet->identify();
249 
250  // if (Verbosity() >= HFJetTruthTrigger::VERBOSITY_MORE)
251  // std::cout << __PRETTY_FUNCTION__ << " jet_flavor = " << jet_flavor
252  // << std::endl;
253 
254  return jet_flavor;
255 }
256 
257 
258 
259 
260 
261 
262 
263 
264 
265 
266 
267 
268 
269 
270 
271 
272 
273 
274 
275 
276 
277 //____________________________________________________________________________..
278 int TruthJetTagging::ResetEvent(PHCompositeNode *topNode)
279 {
280  return Fun4AllReturnCodes::EVENT_OK;
281 }
282 
283 //____________________________________________________________________________..
284 int TruthJetTagging::EndRun(const int runnumber)
285 {
286  return Fun4AllReturnCodes::EVENT_OK;
287 }
288 
289 //____________________________________________________________________________..
290 int TruthJetTagging::End(PHCompositeNode *topNode)
291 {
292  return Fun4AllReturnCodes::EVENT_OK;
293 }
294 
295 //____________________________________________________________________________..
296 int TruthJetTagging::Reset(PHCompositeNode *topNode)
297 {
298  return Fun4AllReturnCodes::EVENT_OK;
299 }
300 
301 //____________________________________________________________________________..
302 void TruthJetTagging::Print(const std::string &what) const
303 {
304  std::cout << "TruthJetTagging::Print(const std::string &what) const Printing info for " << what << std::endl;
305 }