Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EvtTools.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EvtTools.h
1 // ----------------------------------------------------------------------------
2 // 'EvtTools.h'
3 // Derek Anderson
4 // 10.30.2023
5 //
6 // Collection of frequent event-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 <vector>
19 #include <cassert>
20 #include <utility>
21 #include <optional>
22 // root libraries
23 #include <Math/Vector3D.h>
24 // phool libraries
25 #include <phool/phool.h>
26 #include <phool/getClass.h>
27 #include <phool/PHIODataNode.h>
28 #include <phool/PHNodeIterator.h>
29 #include <phool/PHCompositeNode.h>
30 // tracking includes
31 #include <trackbase_historic/SvtxTrack.h>
32 #include <trackbase_historic/SvtxTrackMap.h>
33 // calo includes
34 #include <calobase/RawCluster.h>
35 #include <calobase/RawClusterUtility.h>
36 #include <calobase/RawClusterContainer.h>
37 // hepmc includes
38 #include <HepMC/GenEvent.h>
39 #include <HepMC/GenVertex.h>
40 #include <HepMC/GenParticle.h>
41 #include <phhepmc/PHHepMCGenEvent.h>
42 #include <phhepmc/PHHepMCGenEventMap.h>
43 // vertex libraries
44 #include <globalvertex/GlobalVertex.h>
45 #include <globalvertex/GlobalVertexMap.h>
46 // analysis utilities
47 #include "Constants.h"
48 #include "Interfaces.h"
49 #include "CalTools.h"
50 #include "TrkTools.h"
51 #include "VtxTools.h"
52 #include "GenTools.h"
53 
54 #pragma GCC diagnostic pop
55 
56 // make common namespaces implicit
57 using namespace std;
58 using namespace findNode;
59 
60 
61 
62 namespace SColdQcdCorrelatorAnalysis {
63  namespace SCorrelatorUtilities {
64 
65  // forward declarations used in Reco/GenInfo
66  long GetNumTrks(PHCompositeNode* topNode);
67  long GetNumFinalStatePars(PHCompositeNode* topNode, const vector<int> evtsToGrab, const int subset = 0, optional<float> chargeToGrab = nullopt);
68  double GetSumTrkMomentum(PHCompositeNode* topNode);
69  double GetSumCaloEne(PHCompositeNode* topNode, const string store);
70  double GetSumFinalStateParEne(PHCompositeNode* topNode, const vector<int> evtsToGrab, const int subset, optional<float> chargeToGrab = nullopt);
71  ParInfo GetPartonInfo(PHCompositeNode* topNode, const int event, const int status);
72  ROOT::Math::XYZVector GetRecoVtx(PHCompositeNode* topNode);
73 
74 
75  // RecoInfo definition ----------------------------------------------------
76 
77  struct RecoInfo {
78 
79  // data members
80  int nTrks = numeric_limits<int>::max();
81  double pSumTrks = numeric_limits<double>::max();
82  double eSumEMCal = numeric_limits<double>::max();
83  double eSumIHCal = numeric_limits<double>::max();
84  double eSumOHCal = numeric_limits<double>::max();
85  double vx = numeric_limits<double>::max();
86  double vy = numeric_limits<double>::max();
87  double vz = numeric_limits<double>::max();
88  double vr = numeric_limits<double>::max();
89 
90  void SetInfo(PHCompositeNode* topNode) {
91 
92  // get sums
93  nTrks = GetNumTrks(topNode);
94  pSumTrks = GetSumTrkMomentum(topNode);
95  eSumEMCal = GetSumCaloEne(topNode, "CLUSTER_CEMC");
96  eSumIHCal = GetSumCaloEne(topNode, "CLUSTER_HCALIN");
97  eSumOHCal = GetSumCaloEne(topNode, "CLUSTER_HCALOUT");
98 
99  // get vertex
100  ROOT::Math::XYZVector vtx = GetRecoVtx(topNode);
101  vx = vtx.x();
102  vy = vtx.y();
103  vz = vtx.z();
104  vr = hypot(vx, vy);
105  return;
106  } // end 'SetInfo(PHCompositeNode*)'
107 
108  void Reset() {
109  nTrks = numeric_limits<int>::max();
110  pSumTrks = numeric_limits<double>::max();
111  eSumEMCal = numeric_limits<double>::max();
112  eSumIHCal = numeric_limits<double>::max();
113  eSumOHCal = numeric_limits<double>::max();
114  vx = numeric_limits<double>::max();
115  vy = numeric_limits<double>::max();
116  vz = numeric_limits<double>::max();
117  vr = numeric_limits<double>::max();
118  return;
119  } // end 'Reset()'
120 
121  static vector<string> GetListOfMembers() {
122  vector<string> members = {
123  "nTrks",
124  "pSumTrks",
125  "eSumEMCal",
126  "eSumIHCal",
127  "eSumOHCal",
128  "vx",
129  "vy",
130  "vz",
131  "vr"
132  };
133  return members;
134  } // end 'GetListOfMembers()'
135 
136  // default ctor/dtor
137  RecoInfo() {};
138  ~RecoInfo() {};
139 
140  // ctor accepting PHCompositeNode*
141  RecoInfo(PHCompositeNode* topNode) {
142  SetInfo(topNode);
143  };
144 
145  }; // end RecoInfo definition
146 
147 
148 
149  // GenInfo definition ------------------------------------------------------
150 
151  struct GenInfo {
152 
153  // atomic data members
154  int nChrgPar = numeric_limits<int>::max();
155  int nNeuPar = numeric_limits<int>::max();
156  bool isEmbed = false;
157  double eSumChrg = numeric_limits<double>::max();
158  double eSumNeu = numeric_limits<double>::max();
159 
160  // hard scatter products
161  pair<ParInfo, ParInfo> partons;
162 
163  void SetInfo(PHCompositeNode* topNode, const bool embed, const vector<int> evtsToGrab) {
164 
165  // set embed flag
166  isEmbed = embed;
167 
168  // set parton info
169  isEmbed = embed;
170  if (isEmbed) {
171  partons.first = GetPartonInfo(topNode, SubEvt::EmbedSignal, HardScatterStatus::First);
172  partons.second = GetPartonInfo(topNode, SubEvt::EmbedSignal, HardScatterStatus::Second);
173 
174  } else {
175  partons.first = GetPartonInfo(topNode, SubEvt::NotEmbedSignal, HardScatterStatus::First);
176  partons.second = GetPartonInfo(topNode, SubEvt::NotEmbedSignal, HardScatterStatus::Second);
177  }
178 
179  // get sums
180  nChrgPar = GetNumFinalStatePars(topNode, evtsToGrab, Subset::Charged);
181  nNeuPar = GetNumFinalStatePars(topNode, evtsToGrab, Subset::Neutral);
182  eSumChrg = GetSumFinalStateParEne(topNode, evtsToGrab, Subset::Charged);
183  eSumNeu = GetSumFinalStateParEne(topNode, evtsToGrab, Subset::Neutral);
184  return;
185  } // end 'SetInfo(PHCompositeNode*, vector<int>)'
186 
187  void Reset() {
188  nChrgPar = numeric_limits<int>::max();
189  nNeuPar = numeric_limits<int>::max();
190  isEmbed = false;
191  eSumChrg = numeric_limits<double>::max();
192  eSumNeu = numeric_limits<double>::max();
193  return;
194  } // end 'Reset()'
195 
196  static vector<string> GetListOfMembers() {
197  vector<string> members = {
198  "nChrgPar",
199  "nNeuPar",
200  "isEmbed",
201  "eSumChrg",
202  "eSumNeu"
203  };
204  AddLeavesToVector<ParInfo>(members, "PartonA");
205  AddLeavesToVector<ParInfo>(members, "PartonB");
206  return members;
207  } // end 'GetListOfMembers()'
208 
209  // default ctor/dtor
210  GenInfo() {};
211  ~GenInfo() {};
212 
213  // ctor accepting PHCompositeNode* and list of subevents
214  GenInfo(PHCompositeNode* topNode, const bool embed, vector<int> evtsToGrab) {
215  SetInfo(topNode, embed, evtsToGrab);
216  };
217  }; // end GenInfo definition
218 
219 
220 
221  // event methods ----------------------------------------------------------
222 
223  long GetNumTrks(PHCompositeNode* topNode) {
224 
225  // grab size of track map
226  SvtxTrackMap* mapTrks = GetTrackMap(topNode);
227  return mapTrks -> size();
228 
229  } // end 'GetNumTrks(PHCompositeNode*)'
230 
231 
232 
233  long GetNumFinalStatePars(PHCompositeNode* topNode, const vector<int> evtsToGrab, const int subset, optional<float> chargeToGrab) {
234 
235  // loop over subevents
236  long nPar = 0;
237  for (const int evtToGrab : evtsToGrab) {
238 
239  // loop over particles
240  HepMC::GenEvent* genEvt = GetGenEvent(topNode, evtToGrab);
241  for (
242  HepMC::GenEvent::particle_const_iterator particle = genEvt -> particles_begin();
243  particle != genEvt -> particles_end();
244  ++particle
245  ) {
246 
247  // check if particle is final state
248  if (!IsFinalState((*particle) -> status())) continue;
249 
250  // if chargeToGrab is set, select only particle with charge
251  const float charge = GetParticleCharge((*particle) -> pdg_id());
252  if (chargeToGrab.has_value()) {
253  if (charge == chargeToGrab.value()) {
254  ++nPar;
255  }
256  } else {
257  switch (subset) {
258 
259  // everything
260  case Subset::All:
261  ++nPar;
262  break;
263 
264  // all charged
265  case Subset::Charged:
266  if (charge != 0.) {
267  ++nPar;
268  }
269  break;
270 
271  // only neutral
272  case Subset::Neutral:
273  if (charge == 0.) {
274  ++nPar;
275  }
276  break;
277 
278  default:
279  ++nPar;
280  break;
281  }
282  } // end if-else chargeToGrab.has_value()
283  } // end particle loop
284  } // end subevent loop
285  return nPar;
286 
287  } // end 'GetNumFinalStatePars(PHCompositeNode*, vector<int>, optional<float>, optional<float> bool)'
288 
289 
290 
291  double GetSumTrkMomentum(PHCompositeNode* topNode) {
292 
293  // grab track map
294  SvtxTrackMap* mapTrks = GetTrackMap(topNode);
295 
296  // loop over tracks
297  double pSum = 0.;
298  for (
299  SvtxTrackMap::Iter itTrk = mapTrks -> begin();
300  itTrk != mapTrks -> end();
301  ++itTrk
302  ) {
303 
304  // grab track
305  SvtxTrack* track = itTrk -> second;
306  if (!track) continue;
307 
308  // momentum to sum
309  const double pTrk = std::hypot(track -> get_px(), track -> get_py(), track -> get_pz());
310  pSum += pTrk;
311  }
312  return pSum;
313 
314  } // end 'GetSumTrkMomentum(PHCompositeNode*)'
315 
316 
317 
318  double GetSumCaloEne(PHCompositeNode* topNode, const string store) {
319 
320  // grab clusters
321  RawClusterContainer::ConstRange clusters = GetClusters(topNode, store);
322 
323  // loop over clusters
324  double eSum = 0.;
325  for (
326  RawClusterContainer::ConstIterator itClust = clusters.first;
327  itClust != clusters.second;
328  itClust++
329  ) {
330 
331  // grab cluster and increment sum
332  const RawCluster* cluster = itClust -> second;
333  if (!cluster) {
334  continue;
335  } else {
336  eSum += cluster -> get_energy();
337  }
338  } // end cluster loop
339  return eSum;
340 
341  } // end 'GetSumCaloEne(PHCompositeNode*, string)'
342 
343 
344 
345  double GetSumFinalStateParEne(PHCompositeNode* topNode, const vector<int> evtsToGrab, const int subset, optional<float> chargeToGrab) {
346 
347  // loop over subevents
348  double eSum = 0.;
349  for (const int evtToGrab : evtsToGrab) {
350 
351  HepMC::GenEvent* genEvt = GetGenEvent(topNode, evtToGrab);
352  for (
353  HepMC::GenEvent::particle_const_iterator particle = genEvt -> particles_begin();
354  particle != genEvt -> particles_end();
355  ++particle
356  ) {
357 
358  // check if particle is final state
359  if (!IsFinalState((*particle) -> status())) continue;
360 
361  // if chargeToGrab is set, select only particle with charge
362  const float charge = GetParticleCharge((*particle) -> pdg_id());
363  const float energy = (*particle) -> momentum().e();
364  if (chargeToGrab.has_value()) {
365  if (charge == chargeToGrab.value()) {
366  eSum += energy;
367  }
368  } else {
369  switch (subset) {
370 
371  // everything
372  case Subset::All:
373  eSum += energy;
374  break;
375 
376  // all charged
377  case Subset::Charged:
378  if (charge != 0.) {
379  eSum += energy;
380  }
381  break;
382 
383  // only neutral
384  case Subset::Neutral:
385  if (charge == 0.) {
386  eSum += energy;
387  }
388  break;
389 
390  default:
391  eSum += energy;
392  break;
393  }
394  } // end if-else chargeToGrab.has_value()
395  } // end particle loop
396  } // end subevent loop
397  return eSum;
398 
399  } // end 'GetSumFinalStateParEne(PHCompositeNode*, vector<int>, int, optional<float>)'
400 
401 
402 
403  ParInfo GetPartonInfo(PHCompositeNode* topNode, const int event, const int status) {
404 
405  // pick out relevant sub-sevent to grab
406  HepMC::GenEvent* genEvt = GetGenEvent(topNode, event);
407 
408  // loop over particles
409  ParInfo parton;
410  for (
411  HepMC::GenEvent::particle_const_iterator particle = genEvt -> particles_begin();
412  particle != genEvt -> particles_end();
413  ++particle
414  ) {
415 
416  // ignore all non-partons
417  if (!IsParton((*particle) -> pdg_id())) continue;
418 
419  // set info if parton is desired status
420  if (((*particle) -> status()) == status) {
421  parton.SetInfo(*particle, event);
422  }
423  } // end particle loop
424  return parton;
425 
426  } // end 'GetPartonInfo(PHCompositeNode*, int, int)'
427 
428  } // end SCorrelatorUtilities namespace
429 } // end SColdQcdCorrealtorAnalysis namespace
430 
431 // end ------------------------------------------------------------------------