Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SCorrelatorJetTree.cst.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SCorrelatorJetTree.cst.h
1 // ----------------------------------------------------------------------------
2 // 'SCorrelatorJetTree.csts.h'
3 // Derek Anderson
4 // 01.18.2023
5 //
6 // A module to produce a tree of jets for the sPHENIX
7 // Cold QCD Energy-Energy Correlator analysis.
8 //
9 // Derived from code by Antonio Silva (thanks!!)
10 // ----------------------------------------------------------------------------
11 
12 #pragma once
13 
14 using namespace std;
15 using namespace findNode;
16 
17 
18 
19 namespace SColdQcdCorrelatorAnalysis {
20 
21  // constituent methods ------------------------------------------------------
22 
23  bool SCorrelatorJetTree::IsGoodParticle(HepMC::GenParticle* par, const bool ignoreCharge) {
24 
25  // print debug statement
26  if (m_doDebug && (Verbosity() > 1)) {
27  cout << "SCorrelatorJetTree::IsGoodParticle(HepMC::GenParticle*) Checking if MC particle is good..." << endl;
28  }
29 
30  // check charge if needed
31  const bool isJetCharged = (m_jetType != 1);
32  const bool doChargeCheck = (isJetCharged && !ignoreCharge);
33 
34  int parID;
35  bool isGoodCharge;
36  float parChrg;
37  if (doChargeCheck) {
38  parID = par -> pdg_id();
39  parChrg = GetParticleCharge(parID);
40  isGoodCharge = (parChrg != 0.);
41  } else {
42  isGoodCharge = true;
43  }
44 
45  const double parEta = par -> momentum().eta();
46  const double parPx = par -> momentum().px();
47  const double parPy = par -> momentum().py();
48  const double parPt = sqrt((parPx * parPx) + (parPy * parPy));
49  const bool isInPtRange = ((parPt > m_parPtRange[0]) && (parPt < m_parPtRange[1]));
50  const bool isInEtaRange = ((parEta > m_parEtaRange[0]) && (parEta < m_parEtaRange[1]));
51  const bool isGoodPar = (isGoodCharge && isInPtRange && isInEtaRange);
52  return isGoodPar;
53 
54  } // end 'IsGoodParticle(HepMC::GenParticle*, bool)'
55 
56 
57 
58  bool SCorrelatorJetTree::IsGoodTrack(SvtxTrack* track, PHCompositeNode* topNode) {
59 
60  // print debug statement
61  if (m_doDebug && (Verbosity() > 1)) {
62  cout << "SCorrelatorJetTree::IsGoodTrack(SvtxTrack*) Checking if track is good..." << endl;
63  }
64 
65  // grab track info
66  const double trkPt = track -> get_pt();
67  const double trkEta = track -> get_eta();
68  const double trkQual = track -> get_quality();
69  const double trkDeltaPt = GetTrackDeltaPt(track);
70  const int trkNMvtx = GetNumLayer(track, SUBSYS::MVTX);
71  const int trkNIntt = GetNumLayer(track, SUBSYS::INTT);
72  const int trkNTpc = GetNumLayer(track, SUBSYS::TPC);
73 
74  // grab track dca
75  const auto trkDca = GetTrackDcaPair(track, topNode);
76  const double trkDcaXY = trkDca.first;
77  const double trkDcaZ = trkDca.second;
78 
79 
80  // if above max pt used to fit dca width,
81  // use value of fit at max pt
82  double ptEvalXY = (trkPt > m_dcaPtFitMaxXY) ? m_dcaPtFitMaxXY : trkPt;
83  double ptEvalZ = (trkPt > m_dcaPtFitMaxZ) ? m_dcaPtFitMaxZ : trkPt;
84 
85  // check if dca is good
86  bool isInDcaRangeXY = false;
87  bool isInDcaRangeZ = false;
88  if (m_doDcaSigmaCut) {
89  isInDcaRangeXY = (abs(trkDcaXY) < (m_nSigCutXY * (m_fSigDcaXY -> Eval(ptEvalXY))));
90  isInDcaRangeZ = (abs(trkDcaZ) < (m_nSigCutZ * (m_fSigDcaZ -> Eval(ptEvalZ))));
91  } else {
92  isInDcaRangeXY = ((trkDcaXY > m_trkDcaRangeXY[0]) && (trkDcaXY < m_trkDcaRangeXY[1]));
93  isInDcaRangeZ = ((trkDcaZ > m_trkDcaRangeZ[0]) && (trkDcaZ < m_trkDcaRangeZ[1]));
94  }
95 
96  // if applying vertex cuts, grab track
97  // vertex and check if good
98  bool isInVtxRange = true;
99  if (m_doVtxCut) {
100  CLHEP::Hep3Vector trkVtx = GetTrackVertex(track, topNode);
101  isInVtxRange = IsGoodVertex(trkVtx);
102  }
103 
104  // if using only primary vertex,
105  // ignore tracks from other vertices
106  if (m_useOnlyPrimVtx) {
107  const bool isFromPrimVtx = IsFromPrimaryVtx(track, topNode);
108  if (!isFromPrimVtx) {
109  isInVtxRange = false;
110  }
111  }
112 
113  // if masking tpc sector boundaries,
114  // ignore tracks near boundaries
115  bool isGoodPhi = true;
116  if (m_maskTpcSectors) {
117  isGoodPhi = IsGoodTrackPhi(track);
118  }
119 
120  // apply cuts
121  const bool isSeedGood = IsGoodTrackSeed(track);
122  const bool isInPtRange = ((trkPt > m_trkPtRange[0]) && (trkPt < m_trkPtRange[1]));
123  const bool isInEtaRange = ((trkEta > m_trkEtaRange[0]) && (trkEta < m_trkEtaRange[1]));
124  const bool isInQualRange = ((trkQual > m_trkQualRange[0]) && (trkQual < m_trkQualRange[1]));
125  const bool isInNMvtxRange = ((trkNMvtx > m_trkNMvtxRange[0]) && (trkNMvtx <= m_trkNMvtxRange[1]));
126  const bool isInNInttRange = ((trkNIntt > m_trkNInttRange[0]) && (trkNIntt <= m_trkNInttRange[1]));
127  const bool isInNTpcRange = ((trkNTpc > m_trkNTpcRange[0]) && (trkNTpc <= m_trkNTpcRange[1]));
128  const bool isInDeltaPtRange = ((trkDeltaPt > m_trkDeltaPtRange[0]) && (trkDeltaPt < m_trkDeltaPtRange[1]));
129  const bool isInNumRange = (isInNMvtxRange && isInNInttRange && isInNTpcRange);
130  const bool isInDcaRange = (isInDcaRangeXY && isInDcaRangeZ);
131  const bool isGoodTrack = (isSeedGood && isGoodPhi && isInPtRange && isInEtaRange && isInQualRange && isInNumRange && isInDcaRange && isInDeltaPtRange && isInVtxRange);
132  return isGoodTrack;
133 
134  } // end 'IsGoodTrack(SvtxTrack*)'
135 
136 
137 
138  bool SCorrelatorJetTree::IsGoodFlow(ParticleFlowElement* flow) {
139 
140  // print debug statement
141  if (m_doDebug && (Verbosity() > 1)) {
142  cout << "SCorrelatorJetTree::IsGoodFlow(ParticleFlowElement*) Checking if particle flow element is good..." << endl;
143  }
144 
145  // TODO: explore particle flow cuts
146  const double pfEta = flow -> get_eta();
147  const bool isInEtaRange = ((pfEta > m_flowEtaRange[0]) && (pfEta < m_flowEtaRange[1]));
148  const bool isGoodFlow = isInEtaRange;
149  return isGoodFlow;
150 
151  } // end 'IsGoodFlow(ParticleFlowElement*)'
152 
153 
154 
155  bool SCorrelatorJetTree::IsGoodECal(CLHEP::Hep3Vector& hepVecECal) {
156 
157  // print debug statement
158  if (m_doDebug && (Verbosity() > 1)) {
159  cout << "SCorrelatorJetTree::IsGoodECal(CLHEP::Hep3Vector&) Checking if ECal cluster is good..." << endl;
160  }
161 
162  const double clustPt = hepVecECal.perp();
163  const double clustEta = hepVecECal.pseudoRapidity();
164  const bool isInPtRange = ((clustPt > m_ecalPtRange[0]) && (clustPt < m_ecalPtRange[1]));
165  const bool isInEtaRange = ((clustEta > m_ecalEtaRange[0]) && (clustEta < m_ecalEtaRange[1]));
166  const bool isGoodClust = (isInPtRange && isInEtaRange);
167  return isGoodClust;
168 
169  } // end 'IsGoodECal(CLHEP::Hep3Vector&)'
170 
171 
172 
173  bool SCorrelatorJetTree::IsGoodHCal(CLHEP::Hep3Vector& hepVecHCal) {
174 
175  // print debug statement
176  if (m_doDebug && (Verbosity() > 1)) {
177  cout << "SCorrelatorJetTree::IsGoodHCal(CLHEP::Hep3Vector&) Checking if HCal cluster is good..." << endl;
178  }
179 
180  // TODO: explore particle cuts. These should vary with particle charge/species.
181  const double clustPt = hepVecHCal.perp();
182  const double clustEta = hepVecHCal.pseudoRapidity();
183  const bool isInPtRange = ((clustPt > m_hcalPtRange[0]) && (clustPt < m_hcalPtRange[1]));
184  const bool isInEtaRange = ((clustEta > m_hcalEtaRange[0]) && (clustEta < m_hcalEtaRange[1]));
185  const bool isGoodClust = (isInPtRange && isInEtaRange);
186  return isGoodClust;
187 
188  } // end 'IsGoodHCal(CLHEP::Hep3Vector&)'
189 
190 
191 
192  bool SCorrelatorJetTree::IsGoodTrackSeed(SvtxTrack* track) {
193 
194  // print debug statement
195  if (m_doDebug && (Verbosity() > 2)) {
196  cout << "SCorrelatorJetTree::IsGoodSeedTrack(SvtxTrack*) Checking if track seed is good..." << endl;
197  }
198 
199  // get track seeds
200  TrackSeed* trkSiSeed = track -> get_silicon_seed();
201  TrackSeed* trkTpcSeed = track -> get_tpc_seed();
202 
203  // check if one or both seeds are present as needed
204  bool isSeedGood = (trkSiSeed && trkTpcSeed);
205  if (!m_requireSiSeeds) {
206  isSeedGood = (trkSiSeed || trkTpcSeed);
207  }
208  return isSeedGood;
209 
210  } // end 'IsGoodSeed(SvtxTrack*)'
211 
212 
213 
214  bool SCorrelatorJetTree::IsGoodTrackPhi(SvtxTrack* track, const float phiMaskSize) {
215 
216  // print debug statement
217  if (m_doDebug && (Verbosity() > 2)) {
218  cout << "SCorrelatorJetTree::IsGoodTrackPhi(SvtxTrack*) Checking if track phi is good..." << endl;
219  }
220 
221  // TPC sector boundaries:
222  // 12 sectors --> ~0.523 rad/sector,
223  // assumed to be symmetric about phi = 0
224  // FIXME move to constant in utilities namespace
225  const array<float, NTpcSector> phiSectorBoundaries = {
226  -2.877,
227  -2.354,
228  -1.831,
229  -1.308,
230  -0.785,
231  -0.262,
232  0.262,
233  0.785,
234  1.308,
235  1.831,
236  2.354,
237  2.877
238  };
239 
240  // flag phi as bad if within boundary +- (phiMaskSize / 2)
241  const double halfMaskSize = phiMaskSize / 2.;
242  const double trkPhi = track -> get_phi();
243 
244  // loop over sector boundaries and check phi
245  bool isGoodPhi = true;
246  for (const float boundary : phiSectorBoundaries) {
247  if ((trkPhi > (boundary - halfMaskSize)) && (trkPhi < (boundary + halfMaskSize))) {
248  isGoodPhi = false;
249  break;
250  }
251  }
252  return isGoodPhi;
253 
254  } // end 'IsGoodTrackPhi(SvtxTrack*, float)'
255 
256 
257 
258  bool SCorrelatorJetTree::IsOutgoingParton(HepMC::GenParticle* par) {
259 
260  // print debug statement
261  if (m_doDebug && (Verbosity() > 2)) {
262  cout << "SCorrelatorJetTree::IsParton(HepMC::GenParticle*) Checking if particle is a parton..." << endl;
263  }
264 
265  // grab particle info
266  const int pid = par -> pdg_id();
267  const int status = par -> status();
268 
269  // check if is outgoing parton
270  const bool isStatusGood = ((status == 23) || (status == 24));
271  const bool isLightQuark = ((pid == 1) || (pid == 2));
272  const bool isStrangeQuark = ((pid == 3) || (pid == 4));
273  const bool isHeavyQuark = ((pid == 5) || (pid == 6));
274  const bool isGluon = (pid == 21);
275  const bool isParton = (isLightQuark || isStrangeQuark || isHeavyQuark || isGluon);
276  const bool isOutgoingParton = (isStatusGood && isParton);
277  return isOutgoingParton;
278 
279  } // end 'IsOutgoingParton(HepMC::GenParticle*)'
280 
281 
282 
283  bool SCorrelatorJetTree::IsFromPrimaryVtx(SvtxTrack* track, PHCompositeNode* topNode) {
284 
285  // print debug statement
286  if (m_doDebug) {
287  cout << "SCorrelatorJetTree::IsFromPrimaryVtx(SvtxTrack*, PHCompositeNode*) Checking if track is from primary vertex..." << endl;
288  }
289 
290  // get id of vertex associated with track
291  const int vtxID = (int) track -> get_vertex_id();
292 
293  // get id of primary vertex
294  GlobalVertex* primVtx = GetGlobalVertex(topNode);
295  const int primVtxID = primVtx -> get_id();
296 
297  // check if from vertex and return
298  const bool isFromPrimVtx = (vtxID == primVtxID);
299  return isFromPrimVtx;
300 
301  } // end 'IsFromPrimaryVtx(SvtTrack*, PHCompsiteNode*)'
302 
303 
304 
305  pair<double, double> SCorrelatorJetTree::GetTrackDcaPair(SvtxTrack* track, PHCompositeNode* topNode) {
306 
307  // print debug statement
308  if (m_doDebug) {
309  cout << "SCorrelatorJetTree::GetTrackDcaPair(SvtxTrack*, PHCompositeNode*) Getting track dca values..." << endl;
310  }
311 
312  // get global vertex and convert to acts vector
313  GlobalVertex* sphxVtx = GetGlobalVertex(topNode);
314  Acts::Vector3 actsVtx = Acts::Vector3(sphxVtx -> get_x(), sphxVtx -> get_y(), sphxVtx -> get_z());
315 
316  // return dca
317  const auto dcaAndErr = TrackAnalysisUtils::get_dca(track, actsVtx);
318  return make_pair(dcaAndErr.first.first, dcaAndErr.second.first);
319 
320  } // end 'GetTrackDcaPair(SvtxTrack*, PHCompositeNode*)'
321 
322 
323 
324  CLHEP::Hep3Vector SCorrelatorJetTree::GetTrackVertex(SvtxTrack* track, PHCompositeNode* topNode) {
325 
326  // print debug statement
327  if (m_doDebug) {
328  cout << "SCorrelatorJetTree::GetTrackVertex(SvtxTrack*, PHCompositeNode*) Getting track vertex..." << endl;
329  }
330 
331  // get vertex associated with track
332  const int vtxID = (int) track -> get_vertex_id();
333  GlobalVertex* vtx = GetGlobalVertex(topNode, vtxID);
334 
335  // return vertex 3-vector
336  CLHEP::Hep3Vector hepVecVtx = CLHEP::Hep3Vector(vtx -> get_x(), vtx -> get_y(), vtx -> get_z());
337  return hepVecVtx;
338 
339  } // end 'GetTrackVertex(SvtxTrack*, PHCompositeNode*)'
340 
341 
342 
343  double SCorrelatorJetTree::GetTrackDeltaPt(SvtxTrack* track) {
344 
345  // print debug statement
346  if (m_doDebug) {
347  cout << "SCorrelatorJetTree::GetTrackDeltaPt(SvtxTrack*) Getting track delta pt..." << endl;
348  }
349 
350  // grab covariances
351  const float trkCovXX = track -> get_error(3, 3);
352  const float trkCovXY = track -> get_error(3, 4);
353  const float trkCovYY = track -> get_error(4, 4);
354 
355  // grab momentum
356  const float trkPx = track -> get_px();
357  const float trkPy = track -> get_py();
358  const float trkPt = track -> get_pt();
359 
360  // calculate delta-pt
361  const float numer = (trkCovXX * trkPx * trkPx) + 2 * (trkCovXY * trkPx * trkPy) + (trkCovYY * trkPy * trkPy);
362  const float denom = (trkPx * trkPx) + (trkPy * trkPy);
363  const float ptDelta2 = numer / denom;
364  const float ptDelta = sqrt(ptDelta2) / trkPt;
365  return ptDelta;
366 
367  } // end 'GetTrackDeltaPt(SvtxTrack*)'
368 
369 
370 
371  float SCorrelatorJetTree::GetParticleCharge(const int pid) {
372 
373  // print debug statement
374  if (m_doDebug && (Verbosity() > 2)) {
375  cout << "SCorrelatorJetTree::GetParticleCharge(int) Grabbing MC particle charge..." << endl;
376  }
377 
378  // particle charge
379  float charge(-999.);
380 
381  switch (abs(pid)) {
382  // e+/e-
383  case 11:
384  charge = 1.;
385  break;
386  // e neutrinos
387  case 12:
388  charge = 0.;
389  break;
390  // mu-/mu+
391  case 13:
392  charge = -1.;
393  break;
394  // mu neutrinos
395  case 14:
396  charge = 0.;
397  break;
398  // tau-/tau+
399  case 15:
400  charge = -1.;
401  break;
402  // tau neutrinos
403  case 16:
404  charge = 0.;
405  break;
406  // photon
407  case 22:
408  charge = 0.;
409  break;
410  // Z0
411  case 23:
412  charge = 0.;
413  break;
414  // W+/W-
415  case 24:
416  charge = 1.;
417  break;
418  // pi0
419  case 111:
420  charge = 0.;
421  break;
422  // pi+/pi-
423  case 211:
424  charge = 1.;
425  break;
426  // K0 (long)
427  case 130:
428  charge = 0.;
429  break;
430  // K0 (short)
431  case 310:
432  charge = 0.;
433  break;
434  // K+/K-
435  case 321:
436  charge = 1.;
437  break;
438  // D+/D-
439  case 441:
440  charge = 1.;
441  break;
442  // D0
443  case 421:
444  charge = 0.;
445  break;
446  // DS+/DS-
447  case 431:
448  charge = 1.;
449  break;
450  // eta
451  case 221:
452  charge = 0.;
453  break;
454  // proton/antiproton
455  case 2212:
456  charge = 1.;
457  break;
458  // neutron
459  case 2112:
460  charge = 0.;
461  break;
462  // lambda
463  case 3122:
464  charge = 0.;
465  break;
466  // sigma+/antisigma+
467  case 3222:
468  charge = 1.;
469  break;
470  // sigma0
471  case 3212:
472  charge = 0.;
473  break;
474  // sigma-/antisigma-
475  case 3112:
476  charge = -1.;
477  break;
478  // xi0
479  case 3322:
480  charge = 0.;
481  break;
482  // deuteron
483  case 700201:
484  charge = 0.;
485  break;
486  // alpha
487  case 700202:
488  charge = 2.;
489  break;
490  // triton
491  case 700301:
492  charge = 0.;
493  break;
494  // he3
495  case 700302:
496  charge = 3.;
497  break;
498  default:
499  charge = 0.;
500  if (m_doDebug && (Verbosity() > 1)) {
501  cerr << "SCorrelatorJetTree::GetParticleCharge(int) WARNING: trying to determine charge of unknown particle (PID = " << pid << ")..." << endl;
502  }
503  break;
504  } // end switch (abs(pid))
505 
506  // if antiparticle, flip charge and return
507  if (pid < 0) {
508  charge *= -1.;
509  }
510  return charge;
511 
512  } // end 'GetParticleCharge(int)'
513 
514 
515 
516  int SCorrelatorJetTree::GetNumLayer(SvtxTrack* track, const uint8_t subsys) {
517 
518  // print debug statement
519  if (m_doDebug && (Verbosity() > 3)) {
520  cout << "SCorrelatorJetTree::GetNumLayer(SvtxTrack*, uint8_t) Grabbing number of track clusters..." << endl;
521  }
522 
523  // issue warning if subsys is not set correctly
524  const bool isSubsysWrong = (subsys > 2);
525  if (isSubsysWrong && m_doDebug && (Verbosity() > 3)) {
526  cerr << "SCorrelatorJetTree::GetNumLayer(SvtxTrack*, uint8_t) PANIC: trying to determine no. of clusters for an unknown subsystem (subsys = " << subsys << ")! Aborting!" << endl;
527  assert(subsys <= 2);
528  }
529 
530  // check if seed is good
531  const bool isSeedGood = IsGoodTrackSeed(track);
532  if (!isSeedGood && m_doDebug && (Verbosity() > 3)) {
533  cerr << "SCorrelatorJetTree::GetNumLayer(SvtxTrack*, uint8_t) PANIC: track seed(s) is (are) no good! Aborting!" << endl;
534  assert(isSeedGood);
535  }
536 
537  // get both track seeds
538  TrackSeed* trkSiSeed = track -> get_silicon_seed();
539  TrackSeed* trkTpcSeed = track -> get_tpc_seed();
540 
541  // select relevant seed
542  const bool hasBothSeeds = (trkSiSeed && trkTpcSeed);
543  const bool hasOnlyTpcSeed = (!trkSiSeed && trkTpcSeed);
544 
545  TrackSeed* seed = NULL;
546  switch (subsys) {
547  case SUBSYS::MVTX:
548  if (hasBothSeeds) seed = trkSiSeed;
549  if (hasOnlyTpcSeed) seed = trkTpcSeed;
550  break;
551  case SUBSYS::INTT:
552  if (hasBothSeeds) seed = trkSiSeed;
553  if (hasOnlyTpcSeed) seed = trkTpcSeed;
554  break;
555  case SUBSYS::TPC:
556  seed = trkTpcSeed;
557  break;
558  }
559  if (!seed && m_doDebug && (Verbosity() > 3)) {
560  cerr << "SCorrelatorJetTree::GetNumLayer(SvtxTrack*, uint8_t) PANIC: couldn't set seed! Aborting!" << endl;
561  assert(seed);
562  }
563 
564  // set min no. of layers
565  const int minInttLayer = CONST::NMvtxLayer;
566  const int minTpcLayer = CONST::NMvtxLayer + CONST::NInttLayer;
567 
568  // reset hit flags
569  switch (subsys) {
570  case SUBSYS::MVTX:
571  for (int iMvtxLayer = 0; iMvtxLayer < CONST::NMvtxLayer; iMvtxLayer++) {
572  isMvtxLayerHit[iMvtxLayer] = false;
573  }
574  break;
575  case SUBSYS::INTT:
576  for (int iInttLayer = 0; iInttLayer < CONST::NInttLayer; iInttLayer++) {
577  isInttLayerHit[iInttLayer] = false;
578  }
579  break;
580  case SUBSYS::TPC:
581  for (int iTpcLayer = 0; iTpcLayer < CONST::NTpcLayer; iTpcLayer++) {
582  isTpcLayerHit[iTpcLayer] = false;
583  }
584  break;
585  }
586 
587  // determine which layers were hit
588  unsigned int layer = 0;
589  unsigned int mvtxLayer = 0;
590  unsigned int inttLayer = 0;
591  unsigned int tpcLayer = 0;
592  for (auto itClustKey = (seed -> begin_cluster_keys()); itClustKey != (seed -> end_cluster_keys()); ++itClustKey) {
593 
594  // grab layer number
595  layer = TrkrDefs::getLayer(*itClustKey);
596 
597  // increment accordingly
598  switch (subsys) {
599  case SUBSYS::MVTX:
600  if (layer < CONST::NMvtxLayer) {
601  mvtxLayer = layer;
602  isMvtxLayerHit[mvtxLayer] = true;
603  }
604  break;
605  case SUBSYS::INTT:
606  if ((layer >= minInttLayer) && (layer < minTpcLayer)) {
607  inttLayer = layer - minInttLayer;
608  isInttLayerHit[inttLayer] = true;
609  }
610  break;
611  case SUBSYS::TPC:
612  if (layer >= minTpcLayer) {
613  tpcLayer = layer - minTpcLayer;
614  isTpcLayerHit[tpcLayer] = true;
615  }
616  break;
617  }
618  } // end cluster loop
619 
620  // get the relevant no. of layers
621  int nLayer = 0;
622  switch (subsys) {
623  case SUBSYS::MVTX:
624  for (int iMvtxLayer = 0; iMvtxLayer < CONST::NMvtxLayer; iMvtxLayer++) {
625  if (isMvtxLayerHit[iMvtxLayer]) ++nLayer;
626  }
627  break;
628  case SUBSYS::INTT:
629  for (int iInttLayer = 0; iInttLayer < CONST::NInttLayer; iInttLayer++) {
630  if (isInttLayerHit[iInttLayer]) ++nLayer;
631  }
632  break;
633  case SUBSYS::TPC:
634  for (int iTpcLayer = 0; iTpcLayer < CONST::NTpcLayer; iTpcLayer++) {
635  if (isTpcLayerHit[iTpcLayer]) ++nLayer;
636  }
637  break;
638  }
639  return nLayer;
640 
641  } // end 'GetNumLayer(SvtxTrack*, uint8_t)'
642 
643 
644 
645  int SCorrelatorJetTree::GetNumClust(SvtxTrack* track, const uint8_t subsys) {
646 
647  // issue warning if subsys is not set correctly
648  const bool isSubsysWrong = (subsys > 2);
649  if (isSubsysWrong && m_doDebug && (Verbosity() > 3)) {
650  cerr << "SCorrelatorJetTree::GetNumLClust(SvtxTrack*, uint8_t) PANIC: trying to determine no. of clusters for an unknown subsystem (subsys = " << subsys << ")! Aborting!" << endl;
651  assert(subsys <= 2);
652  }
653 
654  // check if seed is good
655  const bool isSeedGood = IsGoodTrackSeed(track);
656  if (!isSeedGood && m_doDebug && (Verbosity() > 3)) {
657  cerr << "SCorrelatorJetTree::GetNumLayer(SvtxTrack*, uint8_t) PANIC: track seed(s) is (are) no good! Aborting!" << endl;
658  assert(isSeedGood);
659  }
660 
661  // get both track seeds
662  TrackSeed* trkSiSeed = track -> get_silicon_seed();
663  TrackSeed* trkTpcSeed = track -> get_tpc_seed();
664 
665  // select relevant seed
666  const bool hasBothSeeds = (trkSiSeed && trkTpcSeed);
667  const bool hasOnlyTpcSeed = (!trkSiSeed && trkTpcSeed);
668 
669  TrackSeed* seed = NULL;
670  switch (subsys) {
671  case SUBSYS::MVTX:
672  if (hasBothSeeds) seed = trkSiSeed;
673  if (hasOnlyTpcSeed) seed = trkTpcSeed;
674  break;
675  case SUBSYS::INTT:
676  if (hasBothSeeds) seed = trkSiSeed;
677  if (hasOnlyTpcSeed) seed = trkTpcSeed;
678  break;
679  case SUBSYS::TPC:
680  seed = trkTpcSeed;
681  break;
682  }
683  if (!seed && m_doDebug && (Verbosity() > 3)) {
684  cerr << "SCorrelatorJetTree::GetNumClust(SvtxTrack*, uint8_t) PANIC: couldn't set seed! Aborting!" << endl;
685  assert(seed);
686  }
687 
688  // set min no. of layers
689  const int minInttLayer = CONST::NMvtxLayer;
690  const int minTpcLayer = CONST::NMvtxLayer + CONST::NInttLayer;
691 
692  // determine no. of clusters for a given layer
693  unsigned int layer = 0;
694  unsigned int nCluster = 0;
695  for (auto itClustKey = (seed -> begin_cluster_keys()); itClustKey != (seed -> end_cluster_keys()); ++itClustKey) {
696 
697  // grab layer number
698  layer = TrkrDefs::getLayer(*itClustKey);
699 
700  // increment accordingly
701  switch (subsys) {
702  case SUBSYS::MVTX:
703  if (layer < CONST::NMvtxLayer) {
704  ++nCluster;
705  }
706  break;
707  case SUBSYS::INTT:
708  if ((layer >= minInttLayer) && (layer < minTpcLayer)) {
709  ++nCluster;
710  }
711  break;
712  case SUBSYS::TPC:
713  if (layer >= minTpcLayer) {
714  ++nCluster;
715  }
716  break;
717  }
718  } // end cluster loop
719  return nCluster;
720 
721  } // end 'GetNumClust(SvtxTrack*, uint8_t)'
722 
723 
724 
725  int SCorrelatorJetTree::GetMatchID(SvtxTrack* track) {
726 
727  // print debug statement
728  if (m_doDebug && (Verbosity() > 1)) {
729  cout << "SCorrelatorJetTree::GetMatchID(SvtxTrack*) Grabbing barcode of matching particle..." << endl;
730  }
731 
732  // get best match from truth particles
733  PHG4Particle *bestMatch = m_trackEval -> max_truth_particle_by_nclusters(track);
734 
735  // grab barcode of best match
736  int matchID;
737  if (bestMatch) {
738  matchID = bestMatch -> get_barcode();
739  } else {
740  matchID = -1;
741  }
742  return matchID;
743 
744  } // end 'GetMatchID(SvtxTrack*)'
745 
746 } // end SColdQcdCorrelatorAnalysis namespace
747 
748 // end ------------------------------------------------------------------------