Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TrkTools.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TrkTools.h
1 // ----------------------------------------------------------------------------
2 // 'TrkTools.h'
3 // Derek Anderson
4 // 10.30.2023
5 //
6 // Definition of TrkInfo class and collection of frequent track-
7 // related methods utilized in the sPHENIX Cold QCD Energy-Energy
8 // Correlator analysis.
9 // ----------------------------------------------------------------------------
10 
11 #pragma once
12 
13 // c++ utilities
14 #include <array>
15 #include <limits>
16 #include <string>
17 #include <vector>
18 #include <cassert>
19 #include <utility>
20 // root utilities
21 #include <TF1.h>
22 #include <Math/Vector3D.h>
23 // phool libraries
24 #include <phool/phool.h>
25 #include <phool/getClass.h>
26 #include <phool/PHIODataNode.h>
27 #include <phool/PHNodeIterator.h>
28 #include <phool/PHCompositeNode.h>
29 // tracking libraries
30 #include <trackbase_historic/SvtxTrack.h>
31 #include <trackbase_historic/SvtxTrackMap.h>
32 #include <trackbase_historic/TrackAnalysisUtils.h>
33 // track evaluator utilities
34 #include <g4eval/SvtxTrackEval.h>
35 // vertex libraries
36 #include <globalvertex/GlobalVertex.h>
37 #include <globalvertex/GlobalVertexMap.h>
38 // phenix Geant4 utilities
39 #include <g4main/PHG4Particle.h>
40 // analysis utilities
41 #include "Constants.h"
42 #include "VtxTools.h"
43 
44 // make common namespaces implicit
45 using namespace std;
46 using namespace findNode;
47 
48 
49 
50 namespace SColdQcdCorrelatorAnalysis {
51  namespace SCorrelatorUtilities {
52 
53  // forward declarations used in TrkInfo
54  pair<double, double> GetTrackDcaPair(SvtxTrack* track, PHCompositeNode* topNode);
55  ROOT::Math::XYZVector GetTrackVertex(SvtxTrack* track, PHCompositeNode* topNode);
56  double GetTrackDeltaPt(SvtxTrack* track);
57  int GetNumLayer(SvtxTrack* track, const uint8_t sys);
58  int GetNumClust(SvtxTrack* track, const uint8_t sys);
59 
60  // TrkInfo definition -----------------------------------------------------
61 
62  struct TrkInfo {
63 
64  // data members
65  int id = numeric_limits<int>::max();
66  int nMvtxLayer = numeric_limits<int>::max();
67  int nInttLayer = numeric_limits<int>::max();
68  int nTpcLayer = numeric_limits<int>::max();
69  int nMvtxClust = numeric_limits<int>::max();
70  int nInttClust = numeric_limits<int>::max();
71  int nTpcClust = numeric_limits<int>::max();
72  double phi = numeric_limits<double>::max();
73  double ene = numeric_limits<double>::max();
74  double px = numeric_limits<double>::max();
75  double py = numeric_limits<double>::max();
76  double pz = numeric_limits<double>::max();
77  double pt = numeric_limits<double>::max();
78  double eta = numeric_limits<double>::max();
79  double dcaXY = numeric_limits<double>::max();
80  double dcaZ = numeric_limits<double>::max();
81  double ptErr = numeric_limits<double>::max();
82  double quality = numeric_limits<double>::max();
83  double vtxX = numeric_limits<double>::max();
84  double vtxY = numeric_limits<double>::max();
85  double vtxZ = numeric_limits<double>::max();
86 
87  void SetInfo(SvtxTrack* track, PHCompositeNode* topNode) {
88 
89  // do relevant calculations
90  const ROOT::Math::XYZVector trkVtx = GetTrackVertex(track, topNode);
91  const pair<double, double> trkDcaPair = GetTrackDcaPair(track, topNode);
92 
93  // set track info
94  id = track -> get_id();
95  quality = track -> get_quality();
96  eta = track -> get_eta();
97  phi = track -> get_phi();
98  px = track -> get_px();
99  py = track -> get_py();
100  pz = track -> get_pz();
101  pt = track -> get_pt();
102  ene = sqrt((px * px) + (py * py) + (pz * pz) + (MassPion * MassPion));
103  vtxX = trkVtx.x();
104  vtxY = trkVtx.y();
105  vtxZ = trkVtx.z();
106  dcaXY = trkDcaPair.first;
107  dcaZ = trkDcaPair.second;
108  nMvtxLayer = GetNumLayer(track, Subsys::Mvtx);
109  nInttLayer = GetNumLayer(track, Subsys::Intt);
110  nTpcLayer = GetNumLayer(track, Subsys::Tpc);
111  nMvtxClust = GetNumClust(track, Subsys::Mvtx);
112  nInttClust = GetNumClust(track, Subsys::Intt);
113  nTpcClust = GetNumClust(track, Subsys::Tpc);
114  ptErr = GetTrackDeltaPt(track);
115  return;
116  } // end 'SetInfo(SvtxTrack*, PHCompositeNode*)'
117 
118  void Reset() {
119  id = numeric_limits<int>::max();
120  nMvtxLayer = numeric_limits<int>::max();
121  nInttLayer = numeric_limits<int>::max();
122  nTpcLayer = numeric_limits<int>::max();
123  nMvtxClust = numeric_limits<int>::max();
124  nInttClust = numeric_limits<int>::max();
125  nTpcClust = numeric_limits<int>::max();
126  eta = numeric_limits<double>::max();
127  phi = numeric_limits<double>::max();
128  px = numeric_limits<double>::max();
129  py = numeric_limits<double>::max();
130  pz = numeric_limits<double>::max();
131  pt = numeric_limits<double>::max();
132  ene = numeric_limits<double>::max();
133  dcaXY = numeric_limits<double>::max();
134  dcaZ = numeric_limits<double>::max();
135  ptErr = numeric_limits<double>::max();
136  quality = numeric_limits<double>::max();
137  vtxX = numeric_limits<double>::max();
138  vtxY = numeric_limits<double>::max();
139  vtxZ = numeric_limits<double>::max();
140  return;
141  } // end 'Reset()'
142 
143  static vector<string> GetListOfMembers() {
144  vector<string> members = {
145  "id",
146  "nMvtxLayer",
147  "nInttLayer",
148  "nTpcLayer",
149  "nMvtxClust",
150  "nInttClust",
151  "nTpcClust",
152  "eta",
153  "phi",
154  "px",
155  "py",
156  "pz",
157  "pt",
158  "ene",
159  "dcaXY",
160  "dcaZ",
161  "ptErr",
162  "quality",
163  "vtxX",
164  "vtxY",
165  "vtxZ"
166  };
167  return members;
168  } // end 'GetListOfMembers()'
169 
170  // overloaded < operator
171  friend bool operator<(const TrkInfo& lhs, const TrkInfo& rhs) {
172 
173  // note that some quantities aren't relevant for this comparison
174  const bool isLessThan = (
175  (lhs.nMvtxLayer < rhs.nMvtxLayer) &&
176  (lhs.nInttLayer < rhs.nInttLayer) &&
177  (lhs.nTpcLayer < rhs.nTpcLayer) &&
178  (lhs.nMvtxClust < rhs.nMvtxClust) &&
179  (lhs.nInttClust < rhs.nInttClust) &&
180  (lhs.nTpcClust < rhs.nTpcClust) &&
181  (lhs.eta < rhs.eta) &&
182  (lhs.phi < rhs.phi) &&
183  (lhs.px < rhs.px) &&
184  (lhs.py < rhs.py) &&
185  (lhs.pz < rhs.pz) &&
186  (lhs.pt < rhs.pt) &&
187  (lhs.ene < rhs.ene) &&
188  (lhs.dcaXY < rhs.dcaXY) &&
189  (lhs.dcaZ < rhs.dcaZ) &&
190  (lhs.ptErr < rhs.ptErr) &&
191  (lhs.quality < rhs.quality) &&
192  (lhs.vtxX < rhs.vtxX) &&
193  (lhs.vtxY < rhs.vtxY) &&
194  (lhs.vtxZ < rhs.vtxZ)
195  );
196  return isLessThan;
197 
198  } // end 'operator<(TrkInfo&, TrkInfo&)'
199 
200  // overloaded > operator
201  friend bool operator>(const TrkInfo& lhs, const TrkInfo& rhs) {
202 
203  // note that some quantities aren't relevant for this comparison
204  const bool isGreaterThan = (
205  (lhs.nMvtxLayer > rhs.nMvtxLayer) &&
206  (lhs.nInttLayer > rhs.nInttLayer) &&
207  (lhs.nTpcLayer > rhs.nTpcLayer) &&
208  (lhs.nMvtxClust > rhs.nMvtxClust) &&
209  (lhs.nInttClust > rhs.nInttClust) &&
210  (lhs.nTpcClust > rhs.nTpcClust) &&
211  (lhs.eta > rhs.eta) &&
212  (lhs.phi > rhs.phi) &&
213  (lhs.px > rhs.px) &&
214  (lhs.py > rhs.py) &&
215  (lhs.pz > rhs.pz) &&
216  (lhs.pt > rhs.pt) &&
217  (lhs.ene > rhs.ene) &&
218  (lhs.dcaXY > rhs.dcaXY) &&
219  (lhs.dcaZ > rhs.dcaZ) &&
220  (lhs.ptErr > rhs.ptErr) &&
221  (lhs.quality > rhs.quality) &&
222  (lhs.vtxX > rhs.vtxX) &&
223  (lhs.vtxY > rhs.vtxY) &&
224  (lhs.vtxZ > rhs.vtxZ)
225  );
226  return isGreaterThan;
227 
228  } // end 'operator>(TrkInfo&, TrkInfo&)'
229 
230  // overloaded, <=, >= operators
231  inline friend bool operator<=(const TrkInfo& lhs, const TrkInfo& rhs) {return !(lhs > rhs);}
232  inline friend bool operator>=(const TrkInfo& lhs, const TrkInfo& rhs) {return !(lhs < rhs);}
233 
234  // default ctor/dtor
235  TrkInfo() {};
236  ~TrkInfo() {};
237 
238  // ctor accepting SvtxTracks
239  TrkInfo(SvtxTrack* track, PHCompositeNode* topNode) {
240  SetInfo(track, topNode);
241  }
242 
243  }; // end TrkInfo def
244 
245 
246 
247  // track methods ----------------------------------------------------------
248 
249  SvtxTrackMap* GetTrackMap(PHCompositeNode* topNode) {
250 
251  // grab track map
252  SvtxTrackMap* mapTrks = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
253  if (!mapTrks) {
254  cerr << PHWHERE
255  << "PANIC: SvtxTrackMap node is missing!"
256  << endl;
257  assert(mapTrks);
258  }
259  return mapTrks;
260 
261  } // end 'GetTrackMap(PHCompositeNode*)'
262 
263 
264 
265  bool IsInAcceptance(const TrkInfo& trk, const TrkInfo& minimum, const TrkInfo& maximum) {
266 
267  return ((trk >= minimum) && (trk <= maximum));
268 
269  } // end 'IsInAcceptance(TrkInfo&, TrkInfo&, TrkInfo&)'
270 
271 
272 
273  bool IsInSigmaDcaCut(const TrkInfo& trk, const pair<float, float> nSigCut, const pair<float, float> ptFitMax, const pair<TF1*, TF1*> fSigmaDca) {
274 
275  // if above max pt used to fit dca width, use value of fit at max pt
276  const double ptEvalXY = (trk.pt > ptFitMax.first) ? ptFitMax.first : trk.pt;
277  const double ptEvalZ = (trk.pt > ptFitMax.second) ? ptFitMax.second : trk.pt;
278 
279  // check if dca is in cut
280  const bool isInDcaRangeXY = (abs(trk.dcaXY) < (nSigCut.first * (fSigmaDca.first -> Eval(ptEvalXY))));
281  const bool isInDcaRangeZ = (abs(trk.dcaZ) < (nSigCut.second * (fSigmaDca.second -> Eval(ptEvalZ))));
282  const bool isInSigmaDcaCut = (isInDcaRangeXY && isInDcaRangeZ);
283  return isInSigmaDcaCut;
284 
285  } // end 'IsInSigmaDcaCut(TrkInfo&, pair<float, float>, pair<float, float>, pair<TF1*, TF1*>)'
286 
287 
288 
289  bool IsGoodTrackSeed(SvtxTrack* track, const bool requireSiSeeds = true) {
290 
291  // get track seeds
292  TrackSeed* trkSiSeed = track -> get_silicon_seed();
293  TrackSeed* trkTpcSeed = track -> get_tpc_seed();
294 
295  // check if one or both seeds are present as needed
296  bool isSeedGood = (trkSiSeed && trkTpcSeed);
297  if (!requireSiSeeds) {
298  isSeedGood = (trkSiSeed || trkTpcSeed);
299  }
300  return isSeedGood;
301 
302  } // end 'IsGoodSeed(SvtxTrack*)'
303 
304 
305 
306  bool IsFromPrimaryVtx(SvtxTrack* track, PHCompositeNode* topNode) {
307 
308  // get id of vertex associated with track
309  const int vtxID = (int) track -> get_vertex_id();
310 
311  // get id of primary vertex
312  GlobalVertex* primVtx = GetGlobalVertex(topNode);
313  const int primVtxID = primVtx -> get_id();
314 
315  // check if from vertex and return
316  const bool isFromPrimVtx = (vtxID == primVtxID);
317  return isFromPrimVtx;
318 
319  } // end 'IsFromPrimaryVtx(SvtTrack*, PHCompsiteNode*)'
320 
321 
322 
323  pair<double, double> GetTrackDcaPair(SvtxTrack* track, PHCompositeNode* topNode) {
324 
325  // get global vertex and convert to acts vector
326  GlobalVertex* sphxVtx = GetGlobalVertex(topNode);
327  Acts::Vector3 actsVtx = Acts::Vector3(sphxVtx -> get_x(), sphxVtx -> get_y(), sphxVtx -> get_z());
328 
329  // return dca
330  const auto dcaAndErr = TrackAnalysisUtils::get_dca(track, actsVtx);
331  return make_pair(dcaAndErr.first.first, dcaAndErr.second.first);
332 
333  } // end 'GetTrackDcaPair(SvtxTrack*, PHCompositeNode*)'
334 
335 
336 
337  ROOT::Math::XYZVector GetTrackVertex(SvtxTrack* track, PHCompositeNode* topNode) {
338 
339  // get vertex associated with track
340  const int vtxID = (int) track -> get_vertex_id();
341  GlobalVertex* vtx = GetGlobalVertex(topNode, vtxID);
342 
343  // return vertex 3-vector
344  ROOT::Math::XYZVector xyzVtx = ROOT::Math::XYZVector(vtx -> get_x(), vtx -> get_y(), vtx -> get_z());
345  return xyzVtx;
346 
347  } // end 'GetTrackVertex(SvtxTrack*, PHCompositeNode*)'
348 
349 
350 
351  double GetTrackDeltaPt(SvtxTrack* track) {
352 
353  // grab covariances
354  const float trkCovXX = track -> get_error(3, 3);
355  const float trkCovXY = track -> get_error(3, 4);
356  const float trkCovYY = track -> get_error(4, 4);
357 
358  // grab momentum
359  const float trkPx = track -> get_px();
360  const float trkPy = track -> get_py();
361  const float trkPt = track -> get_pt();
362 
363  // calculate delta-pt
364  const float numer = (trkCovXX * trkPx * trkPx) + 2 * (trkCovXY * trkPx * trkPy) + (trkCovYY * trkPy * trkPy);
365  const float denom = (trkPx * trkPx) + (trkPy * trkPy);
366  const float ptDelta2 = numer / denom;
367  const float ptDelta = sqrt(ptDelta2) / trkPt;
368  return ptDelta;
369 
370  } // end 'GetTrackDeltaPt(SvtxTrack*)'
371 
372 
373 
374  TrackSeed* GetTrackSeed(SvtxTrack* track, const uint8_t sys) {
375 
376  // get both track seeds
377  TrackSeed* trkSiSeed = track -> get_silicon_seed();
378  TrackSeed* trkTpcSeed = track -> get_tpc_seed();
379 
380  // select relevant seed
381  const bool hasBothSeeds = (trkSiSeed && trkTpcSeed);
382  const bool hasOnlyTpcSeed = (!trkSiSeed && trkTpcSeed);
383 
384  TrackSeed* seed = NULL;
385  switch (sys) {
386  case Subsys::Mvtx:
387  if (hasBothSeeds) seed = trkSiSeed;
388  if (hasOnlyTpcSeed) seed = trkTpcSeed;
389  break;
390  case Subsys::Intt:
391  if (hasBothSeeds) seed = trkSiSeed;
392  if (hasOnlyTpcSeed) seed = trkTpcSeed;
393  break;
394  case Subsys::Tpc:
395  seed = trkTpcSeed;
396  break;
397  default:
398  seed = NULL;
399  break;
400  }
401  return seed;
402 
403  } // end 'GetTrackSeed(SvtxTrack*, uint8_t)'
404 
405 
406 
407  int GetNumLayer(SvtxTrack* track, const uint8_t sys = 0) {
408 
409  // grab track seed
410  TrackSeed* seed = GetTrackSeed(track, sys);
411 
412  // set min no. of layers
413  const int minInttLayer = NMvtxLayer;
414  const int minTpcLayer = NMvtxLayer + NInttLayer;
415 
416  // initialize hit flags
417  array<bool, NMvtxLayer> isMvtxLayerHit = {false};
418  array<bool, NInttLayer> isInttLayerHit = {false};
419  array<bool, NTpcLayer> isTpcLayerHit = {false};
420 
421  // determine which layers were hit
422  unsigned int layer = 0;
423  unsigned int mvtxLayer = 0;
424  unsigned int inttLayer = 0;
425  unsigned int tpcLayer = 0;
426  for (auto itClustKey = (seed -> begin_cluster_keys()); itClustKey != (seed -> end_cluster_keys()); ++itClustKey) {
427 
428  // grab layer number
429  layer = TrkrDefs::getLayer(*itClustKey);
430 
431  // increment accordingly
432  switch (sys) {
433  case Subsys::Mvtx:
434  if (layer < NMvtxLayer) {
435  mvtxLayer = layer;
436  isMvtxLayerHit[mvtxLayer] = true;
437  }
438  break;
439  case Subsys::Intt:
440  if ((layer >= minInttLayer) && (layer < minTpcLayer)) {
441  inttLayer = layer - minInttLayer;
442  isInttLayerHit[inttLayer] = true;
443  }
444  break;
445  case Subsys::Tpc:
446  if (layer >= minTpcLayer) {
447  tpcLayer = layer - minTpcLayer;
448  isTpcLayerHit[tpcLayer] = true;
449  }
450  break;
451  default:
452  break;
453  }
454  } // end cluster loop
455 
456  // get the relevant no. of layers
457  int nLayer = 0;
458  switch (sys) {
459  case Subsys::Mvtx:
460  for (int iMvtxLayer = 0; iMvtxLayer < NMvtxLayer; iMvtxLayer++) {
461  if (isMvtxLayerHit[iMvtxLayer]) ++nLayer;
462  }
463  break;
464  case Subsys::Intt:
465  for (int iInttLayer = 0; iInttLayer < NInttLayer; iInttLayer++) {
466  if (isInttLayerHit[iInttLayer]) ++nLayer;
467  }
468  break;
469  case Subsys::Tpc:
470  for (int iTpcLayer = 0; iTpcLayer < NTpcLayer; iTpcLayer++) {
471  if (isTpcLayerHit[iTpcLayer]) ++nLayer;
472  }
473  break;
474  default:
475  break;
476  }
477  return nLayer;
478 
479  } // end 'GetNumLayer(SvtxTrack*, uint8_t)'
480 
481 
482 
483  int GetNumClust(SvtxTrack* track, const uint8_t sys = 0) {
484 
485  // grab track seed
486  TrackSeed* seed = GetTrackSeed(track, sys);
487 
488  // set min no. of layers
489  const int minInttLayer = NMvtxLayer;
490  const int minTpcLayer = NMvtxLayer + NInttLayer;
491 
492  // determine no. of clusters for a given layer
493  unsigned int layer = 0;
494  unsigned int nCluster = 0;
495  for (auto itClustKey = (seed -> begin_cluster_keys()); itClustKey != (seed -> end_cluster_keys()); ++itClustKey) {
496 
497  // grab layer number
498  layer = TrkrDefs::getLayer(*itClustKey);
499 
500  // increment accordingly
501  switch (sys) {
502  case Subsys::Mvtx:
503  if (layer < NMvtxLayer) {
504  ++nCluster;
505  }
506  break;
507  case Subsys::Intt:
508  if ((layer >= minInttLayer) && (layer < minTpcLayer)) {
509  ++nCluster;
510  }
511  break;
512  case Subsys::Tpc:
513  if (layer >= minTpcLayer) {
514  ++nCluster;
515  }
516  break;
517  default:
518  break;
519  }
520  } // end cluster loop
521  return nCluster;
522 
523  } // end 'GetNumClust(SvtxTrack*, uint8_t)'
524 
525 
526 
527  int GetMatchID(SvtxTrack* track, SvtxTrackEval* trackEval) {
528 
529  // get best match from truth particles
530  PHG4Particle* bestMatch = trackEval -> max_truth_particle_by_nclusters(track);
531 
532  // grab barcode of best match
533  int matchID;
534  if (bestMatch) {
535  matchID = bestMatch -> get_barcode();
536  } else {
537  matchID = numeric_limits<int>::max();
538  }
539  return matchID;
540 
541  } // end 'GetMatchID(SvtxTrack*)'
542 
543  } // end SCorrelatorUtilities namespace
544 } // end SColdQcdCorrealtorAnalysis namespace
545 
546 // end ------------------------------------------------------------------------