Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GenTools.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GenTools.h
1 // ----------------------------------------------------------------------------
2 // 'GenTools.h'
3 // Derek Anderson
4 // 10.30.2023
5 //
6 // Collection of frequent generator/MC-related methods utilized
7 // in the sPHENIX Cold QCD Energy-Energy Correlator analysis.
8 // ----------------------------------------------------------------------------
9 
10 #pragma once
11 
12 #pragma GCC diagnostic push
13 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
14 
15 // c++ utilities
16 #include <cmath>
17 #include <limits>
18 #include <string>
19 #include <vector>
20 #include <utility>
21 #include <optional>
22 // phool libraries
23 #include <phool/phool.h>
24 #include <phool/getClass.h>
25 #include <phool/PHIODataNode.h>
26 #include <phool/PHNodeIterator.h>
27 #include <phool/PHCompositeNode.h>
28 // hepmc includes
29 #include <HepMC/GenEvent.h>
30 #include <HepMC/GenVertex.h>
31 #include <HepMC/GenParticle.h>
32 #include <phhepmc/PHHepMCGenEvent.h>
33 #include <phhepmc/PHHepMCGenEventMap.h>
34 // analysis utilities
35 #include "Constants.h"
36 
37 #pragma GCC diagnostic pop
38 
39 // make common namespaces implicit
40 using namespace std;
41 using namespace findNode;
42 
43 
44 
45 namespace SColdQcdCorrelatorAnalysis {
46  namespace SCorrelatorUtilities {
47 
48  // ParInfo definition -----------------------------------------------------
49 
50  struct ParInfo {
51 
52  // data members
53  int pid = numeric_limits<int>::max();
54  int status = numeric_limits<int>::max();
55  int barcode = numeric_limits<int>::max();
56  int embedID = numeric_limits<int>::max();
57  float charge = numeric_limits<float>::max();
58  double mass = numeric_limits<double>::max();
59  double eta = numeric_limits<double>::max();
60  double phi = numeric_limits<double>::max();
61  double ene = numeric_limits<double>::max();
62  double px = numeric_limits<double>::max();
63  double py = numeric_limits<double>::max();
64  double pz = numeric_limits<double>::max();
65  double pt = numeric_limits<double>::max();
66  double vx = numeric_limits<double>::max();
67  double vy = numeric_limits<double>::max();
68  double vz = numeric_limits<double>::max();
69  double vr = numeric_limits<double>::max();
70 
71  void SetInfo(const HepMC::GenParticle* particle, const int event) {
72  pid = particle -> pdg_id();
73  status = particle -> status();
74  barcode = particle -> barcode();
75  embedID = event;
76  charge = MapPidOntoCharge[pid];
77  mass = particle -> momentum().m();
78  eta = particle -> momentum().eta();
79  phi = particle -> momentum().phi();
80  ene = particle -> momentum().e();
81  px = particle -> momentum().px();
82  py = particle -> momentum().py();
83  pz = particle -> momentum().pz();
84  pt = particle -> momentum().perp();
85  vx = particle -> production_vertex() -> position().x();
86  vy = particle -> production_vertex() -> position().y();
87  vz = particle -> production_vertex() -> position().z();
88  vr = hypot(vx, vy);
89  return;
90  };
91 
92  void Reset() {
93  pid = numeric_limits<int>::max();
94  status = numeric_limits<int>::max();
95  barcode = numeric_limits<int>::max();
96  embedID = numeric_limits<int>::max();
97  charge = numeric_limits<float>::max();
98  mass = numeric_limits<double>::max();
99  eta = numeric_limits<double>::max();
100  phi = numeric_limits<double>::max();
101  ene = numeric_limits<double>::max();
102  px = numeric_limits<double>::max();
103  py = numeric_limits<double>::max();
104  pz = numeric_limits<double>::max();
105  pt = numeric_limits<double>::max();
106  vx = numeric_limits<double>::max();
107  vy = numeric_limits<double>::max();
108  vz = numeric_limits<double>::max();
109  vr = numeric_limits<double>::max();
110  return;
111  };
112 
113  static vector<string> GetListOfMembers() {
114  vector<string> members = {
115  "pid",
116  "status",
117  "barcode",
118  "embedID",
119  "charge",
120  "mass",
121  "eta",
122  "phi",
123  "ene",
124  "px",
125  "py",
126  "pz",
127  "pt",
128  "vx",
129  "vy",
130  "vz",
131  "vr"
132  };
133  return members;
134  } // end 'GetListOfMembers()'
135 
136  // overloaded < operator
137  friend bool operator<(const ParInfo& lhs, const ParInfo& rhs) {
138 
139  // note that some quantities aren't relevant for this comparison
140  const bool isLessThan = (
141  (lhs.eta < rhs.eta) &&
142  (lhs.phi < rhs.phi) &&
143  (lhs.ene < rhs.ene) &&
144  (lhs.px < rhs.px) &&
145  (lhs.py < rhs.py) &&
146  (lhs.pz < rhs.pz) &&
147  (lhs.pt < rhs.pt)
148  );
149  return isLessThan;
150 
151  } // end 'operator>(ParInfo&, ParInfo&)'
152 
153  // overloaded > operator
154  friend bool operator>(const ParInfo& lhs, const ParInfo& rhs) {
155 
156  // note that some quantities aren't relevant for this comparison
157  const bool isGreaterThan = (
158  (lhs.eta > rhs.eta) &&
159  (lhs.phi > rhs.phi) &&
160  (lhs.ene > rhs.ene) &&
161  (lhs.px > rhs.px) &&
162  (lhs.py > rhs.py) &&
163  (lhs.pz > rhs.pz) &&
164  (lhs.pt > rhs.pt)
165  );
166  return isGreaterThan;
167 
168  } // end 'operator>(ParInfo&, ParInfo&)'
169 
170  // overloaded, <=, >= operators
171  inline friend bool operator<=(const ParInfo& lhs, const ParInfo& rhs) {return !(lhs > rhs);}
172  inline friend bool operator>=(const ParInfo& lhs, const ParInfo& rhs) {return !(lhs < rhs);}
173 
174  // default ctor/dtor
175  ParInfo() {};
176  ~ParInfo() {};
177 
178  ParInfo(HepMC::GenParticle* particle, const int event) {
179  SetInfo(particle, event);
180  };
181 
182  }; // end ParInfo definition
183 
184 
185 
186  // generator/mc methods ---------------------------------------------------
187 
188  PHHepMCGenEventMap* GetMcEventMap(PHCompositeNode* topNode) {
189 
190  PHHepMCGenEventMap* mapMcEvts = getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
191  if (!mapMcEvts) {
192  cerr << PHWHERE
193  << "PANIC: HEPMC event map node is missing!"
194  << endl;
195  assert(mapMcEvts);
196  }
197  return mapMcEvts;
198 
199  } // end 'GetMcEventMap(PHCompositeNode*)'
200 
201 
202 
203  PHHepMCGenEvent* GetMcEvent(PHCompositeNode* topNode, const int iEvtToGrab) {
204 
205  PHHepMCGenEventMap* mcEvtMap = GetMcEventMap(topNode);
206  PHHepMCGenEvent* mcEvt = mcEvtMap -> get(iEvtToGrab);
207  if (!mcEvt) {
208  cerr << PHWHERE
209  << "PANIC: Couldn't grab mc event!"
210  << endl;
211  assert(mcEvt);
212  }
213  return mcEvt;
214 
215  } // end 'GetMcEvent(PHCompositeNode*, int)'
216 
217 
218 
219  HepMC::GenEvent* GetGenEvent(PHCompositeNode* topNode, const int iEvtToGrab) {
220 
221  PHHepMCGenEvent* mcEvt = GetMcEvent(topNode, iEvtToGrab);
222  HepMC::GenEvent* genEvt = mcEvt -> getEvent();
223  if (!genEvt) {
224  cerr << PHWHERE
225  << "PANIC: Couldn't grab HepMC event!"
226  << endl;
227  assert(mcEvt);
228  }
229  return genEvt;
230 
231  } // end 'GetGenEvent(PHCompositeNode*, int)'
232 
233 
234 
235  int GetEmbedID(PHCompositeNode* topNode, const int iEvtToGrab) {
236 
237  // grab mc event & return embedding id
238  PHHepMCGenEvent* mcEvt = GetMcEvent(topNode, iEvtToGrab);
239  return mcEvt -> get_embedding_id();
240 
241  } // end 'GetEmbedID(PHCompositeNode*, int)'
242 
243 
244 
245  bool IsInAcceptance(const ParInfo& particle, const ParInfo& minimum, const ParInfo& maximum) {
246 
247  return ((particle >= minimum) && (particle <= maximum));
248 
249  } // end 'IsInAcceptance(ParInfo&, ParInfo&, ParInfo&)'
250 
251 
252 
253  bool IsFinalState(const int status) {
254 
255  return (status == 1);
256 
257  } // end 'IsFinalState(int)'
258 
259 
260 
261  bool IsHardScatterProduct(const int status) {
262 
263  return ((status == HardScatterStatus::First) || (status == HardScatterStatus::Second));
264 
265  } // end 'IsHardScatterProduct(int)'
266 
267 
268 
269  bool IsParton(const int pid) {
270 
271  const bool isLightQuark = ((pid == Parton::Down) || (pid == Parton::Up));
272  const bool isStrangeQuark = ((pid == Parton::Strange) || (pid == Parton::Charm));
273  const bool isHeavyQuark = ((pid == Parton::Bottom) || (pid == Parton::Top));
274  const bool isGluon = (pid == Parton::Gluon);
275  return (isLightQuark || isStrangeQuark || isHeavyQuark || isGluon);
276 
277  } // end 'IsParton(int)'
278 
279 
280 
281  bool IsOutgoingParton(const HepMC::GenParticle* par) {
282 
283  // grab particle info
284  const int pid = par -> pdg_id();
285  const int status = par -> status();
286 
287  // check if is outgoing parton
288  const bool isOutgoingParton = (IsHardScatterProduct(status) && IsParton(pid));
289  return isOutgoingParton;
290 
291  } // end 'IsOutgoingParton(HepMC::GenParticle*)'
292 
293 
294 
295  float GetParticleCharge(const int pid) {
296 
297  // particle charge
298  float charge = MapPidOntoCharge[abs(pid)];
299 
300  // if antiparticle, flip charge and return
301  if (pid < 0) {
302  charge *= -1.;
303  }
304  return charge;
305 
306  } // end 'GetParticleCharge(int)'
307 
308 
309 
310  vector<int> GrabSubevents(PHCompositeNode* topNode, optional<vector<int>> evtsToGrab = nullopt) {
311 
312  // instantiate vector to hold subevents
313  vector<int> subevents;
314 
315  PHHepMCGenEventMap* mcEvtMap = GetMcEventMap(topNode);
316  for (
317  PHHepMCGenEventMap::ConstIter itEvt = mcEvtMap -> begin();
318  itEvt != mcEvtMap -> end();
319  ++itEvt
320  ) {
321 
322  // grab event id
323  const int embedID = itEvt -> second -> get_embedding_id();
324 
325  // if selecting certain subevents, check if matched
326  bool addToList = false;
327  if (evtsToGrab.has_value()) {
328  for (const int idToCheck : evtsToGrab.value()) {
329  if (embedID == idToCheck) {
330  addToList = true;
331  break;
332  }
333  } // end evtsToGrab loop
334  } else {
335  addToList = true;
336  }
337 
338  // add id to list if needed
339  if (addToList) subevents.push_back(embedID);
340  }
341  return subevents;
342 
343  } // end 'GrabSubevents(PHCompositeNode*, optional<vector<int>>)'
344 
345 
346 
347  bool IsSubEvtGood(const int embedID, const int option, const bool isEmbed) {
348 
349  // set ID of signal
350  int signalID = SubEvt::NotEmbedSignal;
351  if (isEmbed) {
352  signalID = SubEvt::EmbedSignal;
353  }
354 
355  bool isSubEvtGood = true;
356  switch (option) {
357 
358  // consider everything
359  case SubEvtOpt::Everything:
360  isSubEvtGood = true;
361  break;
362 
363  // only consider signal event
364  case SubEvtOpt::OnlySignal:
365  isSubEvtGood = (embedID == signalID);
366  break;
367 
368  // only consider background events
369  case SubEvtOpt::AllBkgd:
370  isSubEvtGood = (embedID <= SubEvt::Background);
371  break;
372 
373  // only consider primary background event
374  case SubEvtOpt::PrimaryBkgd:
375  isSubEvtGood = (embedID == SubEvt::Background);
376  break;
377 
378  // only consider pileup events
379  case SubEvtOpt::Pileup:
380  isSubEvtGood = (embedID < SubEvt::Background);
381  break;
382 
383  // by default do nothing
384  default:
385  isSubEvtGood = true;
386  break;
387  }
388  return isSubEvtGood;
389 
390  } // end 'IsSubEvtGood(int, int, bool)'
391 
392 
393 
394  bool IsSubEvtGood(const int embedID, vector<int> subEvtsToUse) {
395 
396  bool isSubEvtGood = false;
397  for (const int evtToUse : subEvtsToUse) {
398  if (embedID == evtToUse) {
399  isSubEvtGood = true;
400  break;
401  }
402  }
403  return isSubEvtGood;
404 
405  } // end 'IsSubEvtGood(int, vector<int>)'
406 
407  } // end SCorrelatorUtilities namespace
408 } // end SColdQcdCorrealtorAnalysis namespace
409 
410 // end ------------------------------------------------------------------------