Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SCorrelatorJetTree.jet.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SCorrelatorJetTree.jet.h
1 // ----------------------------------------------------------------------------
2 // 'SCorrelatorJetTree.jets.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 fastjet;
16 using namespace findNode;
17 
18 
19 
20 namespace SColdQcdCorrelatorAnalysis {
21 
22  // jet methods --------------------------------------------------------------
23 
24  void SCorrelatorJetTree::FindTrueJets(PHCompositeNode* topNode) {
25 
26  // print debug statement
27  if (m_doDebug) {
28  cout << "SCorrelatorJetTree::FindTrueJets(PHCompositeNode*) Finding truth (inclusive) jets..." << endl;
29  }
30 
31  // define jets
32  m_trueJetDef = new JetDefinition(m_jetAlgo, m_jetR, m_recombScheme, fastjet::Best);
33 
34  // declare constituent vector and jet map
35  vector<fastjet::PseudoJet> particles;
36  map<int, pair<Jet::SRC, int>> fjMapMC;
37 
38  // add constituents
39  AddParticles(topNode, particles, fjMapMC);
40 
41  // cluster jets
42  m_trueClust = new ClusterSequence(particles, *m_trueJetDef);
43  m_trueJets = m_trueClust -> inclusive_jets();
44  return;
45 
46  } // end 'FindTrueJets(PHCompositeNode*)'
47 
48 
49 
50  void SCorrelatorJetTree::FindRecoJets(PHCompositeNode* topNode) {
51 
52  // print debug statement
53  if (m_doDebug) {
54  cout << "SCorrelatorJetTree::FindRecoJets(PHCompositeNode*) Finding jets..." << endl;
55  }
56 
57  // define jets
58  m_recoJetDef = new JetDefinition(m_jetAlgo, m_jetR, m_recombScheme, fastjet::Best);
59 
60  // declare constituent vector and jet map
61  vector<fastjet::PseudoJet> particles;
62  map<int, pair<Jet::SRC, int>> fjMap;
63 
64  // add constitutents
65  if (m_addTracks) AddTracks(topNode, particles, fjMap);
66  if (m_addFlow) AddFlow(topNode, particles, fjMap);
67  if (m_addECal) AddECal(topNode, particles, fjMap);
68  if (m_addHCal) AddHCal(topNode, particles, fjMap);
69 
70  // cluster jets
71  m_recoClust = new ClusterSequence(particles, *m_recoJetDef);
72  m_recoJets = m_recoClust -> inclusive_jets();
73  return;
74 
75  } // end 'FindRecoJets(PHCompositeNode*)'
76 
77 
78 
79  void SCorrelatorJetTree::AddParticles(PHCompositeNode* topNode, vector<PseudoJet>& particles, map<int, pair<Jet::SRC, int>>& fjMap) {
80 
81  // print debug statement
82  if (m_doDebug) {
83  cout << "SCorrelatorJetTree::AddParticles(PHComposite*, vector<PseudoJet>&, map<int, pair<Jet::SRC, int>>&) Adding MC particles..." << endl;
84  }
85 
86  // loop over relevant subevents
87  unsigned int iCst = particles.size();
88  unsigned int nParTot = 0;
89  unsigned int nParAcc = 0;
90  double eParSum = 0.;
91  for (const int evtToGrab : m_vecEvtsToGrab) {
92 
93  // grab subevent
94  HepMC::GenEvent* mcEvt = GetMcEvent(topNode, evtToGrab);
95 
96  // grab embedding ID
97  const int embedID = GetEmbedID(topNode, evtToGrab);
98 
99  // loop over particles in subevent
100  for (HepMC::GenEvent::particle_const_iterator itPar = mcEvt -> particles_begin(); itPar != mcEvt -> particles_end(); ++itPar) {
101 
102  // check if particle is final state
103  const bool isFinalState = ((*itPar) -> status() == 1);
104  if (!isFinalState) {
105  continue;
106  } else {
107  ++nParTot;
108  }
109 
110  // check if particle is good
111  const bool isGoodPar = IsGoodParticle(*itPar);
112  if (!isGoodPar) {
113  continue;
114  } else {
115  ++nParAcc;
116  }
117 
118  // grab particle info
119  const double parPx = (*itPar) -> momentum().px();
120  const double parPy = (*itPar) -> momentum().py();
121  const double parPz = (*itPar) -> momentum().pz();
122  const double parE = (*itPar) -> momentum().e();
123  const int parID = (*itPar) -> barcode();
124 
125  // map barcode onto relevant embeddingID
126  m_mapCstToEmbedID[parID] = embedID;
127 
128  // create pseudojet & add to constituent vector
129  fastjet::PseudoJet fjParticle(parPx, parPy, parPz, parE);
130  fjParticle.set_user_index(parID);
131  particles.push_back(fjParticle);
132 
133 
134  pair<int, pair<Jet::SRC, int>> jetPartPair(iCst, make_pair(Jet::SRC::PARTICLE, parID));
135  fjMap.insert(jetPartPair);
136 
137  // fill QA histograms, increment sums and counters
138  m_hObjectQA[OBJECT::PART][INFO::PT] -> Fill(fjParticle.perp());
139  m_hObjectQA[OBJECT::PART][INFO::ETA] -> Fill(fjParticle.pseudorapidity());
140  m_hObjectQA[OBJECT::PART][INFO::PHI] -> Fill(fjParticle.phi_std());
141  m_hObjectQA[OBJECT::PART][INFO::ENE] -> Fill(fjParticle.E());
142  eParSum += parE;
143  ++iCst;
144  } // end particle loop
145  } // end subevent loop
146 
147  // fill QA histograms
148  m_hNumObject[OBJECT::PART] -> Fill(nParAcc);
149  m_hNumCstAccept[CST_TYPE::PART_CST][0] -> Fill(nParTot);
150  m_hNumCstAccept[CST_TYPE::PART_CST][1] -> Fill(nParAcc);
151  m_hSumCstEne[CST_TYPE::PART_CST] -> Fill(eParSum);
152  return;
153 
154  } // end 'AddParticles(PHCompositeNode*, vector<PseudoJet>&, map<int, pair<Jet::SRC, int>>&)'
155 
156 
157 
158  void SCorrelatorJetTree::AddTracks(PHCompositeNode* topNode, vector<PseudoJet>& particles, map<int, pair<Jet::SRC, int>>& fjMap) {
159 
160  // print debug statement
161  if (m_doDebug) {
162  cout << "SCorrelatorJetTree::AddTracks(PHCompositeNode*, vector<PseudoJet>&, map<int, pair<Jet::SRC, int>>&) Adding tracks..." << endl;
163  }
164 
165  // for weird tracks
166  array<float, 40> arrWeirdTrackLeaves;
167  for (size_t iLeaf = 0; iLeaf < 40; iLeaf++) {
168  arrWeirdTrackLeaves[iLeaf] = 0.;
169  }
170 
171  vector<TrkrDefs::cluskey> clustKeysA;
172  vector<TrkrDefs::cluskey> clustKeysB;
173  clustKeysA.clear();
174  clustKeysB.clear();
175 
176  // loop over tracks
177  unsigned int iCst = particles.size();
178  unsigned int nTrkTot = 0;
179  unsigned int nTrkAcc = 0;
180  double eTrkSum = 0.;
181  SvtxTrack* track = NULL;
182  SvtxTrack* trackB = NULL;
183  SvtxTrackMap* mapTrks = GetTrackMap(topNode);
184  for (SvtxTrackMap::Iter itTrk = mapTrks -> begin(); itTrk != mapTrks -> end(); ++itTrk) {
185 
186  // get track
187  track = itTrk -> second;
188  if (!track) {
189  continue;
190  } else {
191  ++nTrkTot;
192  }
193 
194  // check if good
195  const bool isGoodTrack = IsGoodTrack(track, topNode);
196  if (!isGoodTrack) {
197  continue;
198  } else {
199  ++nTrkAcc;
200  }
201 
202  // create pseudojet and add to constituent vector
203  const int trkID = track -> get_id();
204  const double trkPx = track -> get_px();
205  const double trkPy = track -> get_py();
206  const double trkPz = track -> get_pz();
207  const double trkE = sqrt((trkPx * trkPx) + (trkPy * trkPy) + (trkPz * trkPz) + (0.140 * 0.140)); // FIXME move pion mass to constant in utilities namespace
208 
209  // grab barcode of matching particle
210  int matchID;
211  if (m_isMC) {
212  matchID = GetMatchID(track);
213  } else {
214  matchID = -1;
215  }
216 
217 
218  fastjet::PseudoJet fjTrack(trkPx, trkPy, trkPz, trkE);
219  fjTrack.set_user_index(matchID);
220  particles.push_back(fjTrack);
221 
222  // add track to fastjet map
223  pair<int, pair<Jet::SRC, int>> jetTrkPair(iCst, make_pair(Jet::SRC::TRACK, trkID));
224  fjMap.insert(jetTrkPair);
225 
226  // grab track dca and vertex
227  pair<double, double> trkDcaPair = GetTrackDcaPair(track, topNode);
228  CLHEP::Hep3Vector trkVtx = GetTrackVertex(track, topNode);
229 
230  // grab remaining track info
231  const double trkQuality = track -> get_quality();
232  const double trkDeltaPt = GetTrackDeltaPt(track);
233  const double trkDcaXY = trkDcaPair.first;
234  const double trkDcaZ = trkDcaPair.second;
235  const int trkNumTpc = GetNumLayer(track, SUBSYS::TPC);
236  const int trkNumIntt = GetNumLayer(track, SUBSYS::INTT);
237  const int trkNumMvtx = GetNumLayer(track, SUBSYS::MVTX);
238 
239  const int trkClustTpc = GetNumClust(track, SUBSYS::TPC);
240  const int trkClustIntt = GetNumClust(track, SUBSYS::INTT);
241  const int trkClustMvtx = GetNumClust(track, SUBSYS::MVTX);
242  const double trkPhi = track -> get_phi();
243  const double trkEta = track -> get_eta();
244  const double trkPt = track -> get_pt();
245  if (m_checkWeirdTrks) {
246  for (SvtxTrackMap::Iter itTrkB = mapTrks -> begin(); itTrkB != mapTrks -> end(); ++itTrkB) {
247 
248  // grab track B and skip if bad
249  trackB = itTrkB -> second;
250  if (!trackB) continue;
251 
252  const bool isGoodTrackB = IsGoodTrack(trackB, topNode);
253  if (!isGoodTrackB) continue;
254 
255  // grab track B dca and vertex
256  pair<double, double> trkDcaPairB = GetTrackDcaPair(trackB, topNode);
257  CLHEP::Hep3Vector trkVtxB = GetTrackVertex(trackB, topNode);
258 
259  // get info from track B
260  const int trkIDB = trackB -> get_id();
261  const int trkNumTpcB = GetNumLayer(trackB, SUBSYS::TPC);
262  const int trkNumInttB = GetNumLayer(trackB, SUBSYS::INTT);
263  const int trkNumMvtxB = GetNumLayer(trackB, SUBSYS::MVTX);
264  const int trkClustTpcB = GetNumClust(trackB, SUBSYS::TPC);
265  const int trkClustInttB = GetNumClust(trackB, SUBSYS::INTT);
266  const int trkClustMvtxB = GetNumClust(trackB, SUBSYS::MVTX);
267  const double trkQualityB = trackB -> get_quality();
268  const double trkDeltaPtB = GetTrackDeltaPt(trackB);
269  const double trkDcaXYB = trkDcaPairB.first;
270  const double trkDcaZB = trkDcaPairB.second;
271  const double trkPxB = trackB -> get_px();
272  const double trkPyB = trackB -> get_py();
273  const double trkPzB = trackB -> get_pz();
274  const double trkEB = sqrt((trkPxB * trkPxB) + (trkPyB * trkPyB) + (trkPzB * trkPzB) + (0.140 * 0.140));
275  const double trkPhiB = trackB -> get_phi();
276  const double trkEtaB = trackB -> get_eta();
277  const double trkPtB = trackB -> get_pt();
278 
279  // calculate delta-R
280  const double dfTrkAB = trkPhi - trkPhiB;
281  const double dhTrkAB = trkEta - trkEtaB;
282  const double drTrkAB = sqrt((dfTrkAB * dfTrkAB) + (dhTrkAB * dhTrkAB));
283 
284  // skip if same as track A
285  if (trkID == trkIDB) continue;
286 
287  // clear vectors for checking cluster keys
288  clustKeysA.clear();
289  clustKeysB.clear();
290 
291  // loop over clusters from track A
292  auto seedTpcA = track -> get_tpc_seed();
293  if (seedTpcA) {
294  for (auto local_iterA = seedTpcA -> begin_cluster_keys(); local_iterA != seedTpcA -> end_cluster_keys(); ++local_iterA) {
295  TrkrDefs::cluskey cluster_keyA = *local_iterA;
296  clustKeysA.push_back(cluster_keyA);
297  }
298  }
299 
300  // loop over clusters from track B
301  auto seedTpcB = trackB -> get_tpc_seed();
302  if (seedTpcB) {
303  for (auto local_iterB = seedTpcB -> begin_cluster_keys(); local_iterB != seedTpcB -> end_cluster_keys(); ++local_iterB) {
304  TrkrDefs::cluskey cluster_keyB = *local_iterB;
305  clustKeysB.push_back(cluster_keyB);
306  }
307  }
308 
309  // calculate no. of same cluster keys
310  uint64_t nSameKey = 0;
311  for (auto keyA : clustKeysA) {
312  for (auto keyB : clustKeysB) {
313  if (keyA == keyB) {
314  ++nSameKey;
315  break;
316  }
317  } // end cluster key B loop
318  } // end cluster key A loop
319 
320  // set tuple leaves
321  arrWeirdTrackLeaves[0] = (float) trkID;
322  arrWeirdTrackLeaves[1] = (float) trkIDB;
323  arrWeirdTrackLeaves[2] = (float) trkPt;
324  arrWeirdTrackLeaves[3] = (float) trkPtB;
325  arrWeirdTrackLeaves[4] = (float) trkEta;
326  arrWeirdTrackLeaves[5] = (float) trkEtaB;
327  arrWeirdTrackLeaves[6] = (float) trkPhi;
328  arrWeirdTrackLeaves[7] = (float) trkPhiB;
329  arrWeirdTrackLeaves[8] = (float) trkE;
330  arrWeirdTrackLeaves[9] = (float) trkEB;
331  arrWeirdTrackLeaves[10] = (float) trkDcaXY;
332  arrWeirdTrackLeaves[11] = (float) trkDcaXYB;
333  arrWeirdTrackLeaves[12] = (float) trkDcaZ;
334  arrWeirdTrackLeaves[13] = (float) trkDcaZB;
335  arrWeirdTrackLeaves[14] = (float) trkVtx.x();
336  arrWeirdTrackLeaves[15] = (float) trkVtxB.x();
337  arrWeirdTrackLeaves[16] = (float) trkVtx.y();
338  arrWeirdTrackLeaves[17] = (float) trkVtxB.y();
339  arrWeirdTrackLeaves[18] = (float) trkVtx.z();
340  arrWeirdTrackLeaves[19] = (float) trkVtxB.z();
341  arrWeirdTrackLeaves[20] = (float) trkQuality;
342  arrWeirdTrackLeaves[21] = (float) trkQualityB;
343  arrWeirdTrackLeaves[22] = (float) trkDeltaPt;
344  arrWeirdTrackLeaves[23] = (float) trkDeltaPtB;
345  arrWeirdTrackLeaves[24] = (float) trkNumMvtx;
346  arrWeirdTrackLeaves[25] = (float) trkNumMvtxB;
347  arrWeirdTrackLeaves[26] = (float) trkNumIntt;
348  arrWeirdTrackLeaves[27] = (float) trkNumInttB;
349  arrWeirdTrackLeaves[28] = (float) trkNumTpc;
350  arrWeirdTrackLeaves[29] = (float) trkNumTpcB;
351  arrWeirdTrackLeaves[30] = (float) trkClustMvtx;
352  arrWeirdTrackLeaves[31] = (float) trkClustMvtxB;
353  arrWeirdTrackLeaves[32] = (float) trkClustIntt;
354  arrWeirdTrackLeaves[33] = (float) trkClustInttB;
355  arrWeirdTrackLeaves[34] = (float) trkClustTpc;
356  arrWeirdTrackLeaves[35] = (float) trkClustTpcB;
357  arrWeirdTrackLeaves[36] = (float) clustKeysA.size();
358  arrWeirdTrackLeaves[37] = (float) clustKeysB.size();
359  arrWeirdTrackLeaves[38] = (float) nSameKey;
360  arrWeirdTrackLeaves[39] = (float) drTrkAB;
361 
362  // fill weird track tuple
363  m_ntWeirdTracks -> Fill(arrWeirdTrackLeaves.data());
364 
365  } // end 2nd track loop
366  } // end if (m_checkWeirdTrks)
367 
368  // fill QA histograms
369  m_hObjectQA[OBJECT::TRACK][INFO::PT] -> Fill(fjTrack.perp());
370  m_hObjectQA[OBJECT::TRACK][INFO::ETA] -> Fill(fjTrack.pseudorapidity());
371  m_hObjectQA[OBJECT::TRACK][INFO::PHI] -> Fill(fjTrack.phi_std());
372  m_hObjectQA[OBJECT::TRACK][INFO::ENE] -> Fill(fjTrack.E());
373  m_hObjectQA[OBJECT::TRACK][INFO::QUAL] -> Fill(trkQuality);
374  m_hObjectQA[OBJECT::TRACK][INFO::DCAXY] -> Fill(trkDcaXY);
375  m_hObjectQA[OBJECT::TRACK][INFO::DCAZ] -> Fill(trkDcaZ);
376  m_hObjectQA[OBJECT::TRACK][INFO::DELTAPT] -> Fill(trkDeltaPt);
377  m_hObjectQA[OBJECT::TRACK][INFO::NTPC] -> Fill(trkNumTpc);
378 
379  // fill QA tuple, increment sums and counters
380  m_ntTrkQA -> Fill(
381  (float) fjTrack.perp(),
382  (float) fjTrack.pseudorapidity(),
383  (float) fjTrack.phi_std(),
384  (float) fjTrack.E(),
385  (float) trkDcaXY,
386  (float) trkDcaZ,
387  (float) trkDeltaPt,
388  (float) trkQuality,
389  (float) trkNumMvtx,
390  (float) trkNumIntt,
391  (float) trkNumTpc,
392  (float) trkVtx.x(),
393  (float) trkVtx.y(),
394  (float) trkVtx.z()
395  );
396  eTrkSum += trkE;
397  ++iCst;
398 
399  } // end track loop
400 
401  // fill QA histograms
402  m_hNumObject[OBJECT::TRACK] -> Fill(nTrkAcc);
403  m_hNumCstAccept[CST_TYPE::TRACK_CST][0] -> Fill(nTrkTot);
404  m_hNumCstAccept[CST_TYPE::TRACK_CST][1] -> Fill(nTrkAcc);
405  m_hSumCstEne[CST_TYPE::TRACK_CST] -> Fill(eTrkSum);
406  return;
407 
408  } // end 'AddTracks(PHCompositeNode*, vector<PseudoJet>&, map<int, pair<Jet::SRC, int>>&)'
409 
410 
411 
412  void SCorrelatorJetTree::AddFlow(PHCompositeNode* topNode, vector<PseudoJet>& particles, map<int, pair<Jet::SRC, int>>& fjMap) {
413 
414  // print debug statement
415  if (m_doDebug) {
416  cout << "SCorrelatorJetTree::AddFlow(PHCompositeNode*, vector<PseudoJet>&, map<int, pair<Jet::SRC, int>>&) Adding particle flow elements..." << endl;
417  }
418 
419  // warn if jets should be charged
420  if (m_doDebug && (m_jetType != 1)) {
421  cerr << "SCorrelatorJetTree::AddFlow - Warning - trying to add particle flow elements to charged jets!" << endl;
422  }
423 
424  // loop over pf elements
425  unsigned int iCst = particles.size();
426  unsigned int nFlowTot = 0;
427  unsigned int nFlowAcc = 0;
428  double eFlowSum = 0.;
429  ParticleFlowElementContainer* flowStore = GetFlowStore(topNode);
430  ParticleFlowElementContainer::ConstRange flowRange = flowStore -> getParticleFlowElements();
431  ParticleFlowElementContainer::ConstIterator itFlow;
432  for (itFlow = flowRange.first; itFlow != flowRange.second; ++itFlow) {
433 
434  // get pf element
435  ParticleFlowElement* flow = itFlow -> second;
436  if (!flow) {
437  continue;
438  } else {
439  ++nFlowTot;
440  }
441 
442  // check if good
443  const bool isGoodFlow = IsGoodFlow(flow);
444  if (!isGoodFlow) {
445  continue;
446  } else {
447  ++nFlowAcc;
448  }
449 
450  // create pseudojet and add to constituent vector
451  const int pfID = flow -> get_id();
452  const double pfE = flow -> get_e();
453  const double pfPx = flow -> get_px();
454  const double pfPy = flow -> get_py();
455  const double pfPz = flow -> get_pz();
456 
457  fastjet::PseudoJet fjFlow(pfPx, pfPy, pfPz, pfE);
458  fjFlow.set_user_index(iCst);
459  particles.push_back(fjFlow);
460 
461  // add pf element to fastjet map
462  pair<int, pair<Jet::SRC, int>> jetPartFlowPair(iCst, make_pair(Jet::SRC::PARTICLE, pfID));
463  fjMap.insert(jetPartFlowPair);
464 
465  // fill QA histograms, increment sums and counters
466  m_hObjectQA[OBJECT::FLOW][INFO::PT] -> Fill(fjFlow.perp());
467  m_hObjectQA[OBJECT::FLOW][INFO::ETA] -> Fill(fjFlow.pseudorapidity());
468  m_hObjectQA[OBJECT::FLOW][INFO::PHI] -> Fill(fjFlow.phi_std());
469  m_hObjectQA[OBJECT::FLOW][INFO::ENE] -> Fill(fjFlow.E());
470  eFlowSum += pfE;
471  ++iCst;
472  } // end pf element loop
473 
474  // fill QA histograms
475  m_hNumObject[OBJECT::FLOW] -> Fill(nFlowAcc);
476  m_hNumCstAccept[CST_TYPE::FLOW_CST][0] -> Fill(nFlowTot);
477  m_hNumCstAccept[CST_TYPE::FLOW_CST][1] -> Fill(nFlowAcc);
478  m_hSumCstEne[CST_TYPE::FLOW_CST] -> Fill(eFlowSum);
479  return;
480 
481  } // end 'AddFlow(PHCompositeNode*, vector<PseudoJet>&, map<int, pair<Jet::SRC, int>>&)'
482 
483 
484 
485  void SCorrelatorJetTree::AddECal(PHCompositeNode* topNode, vector<PseudoJet>& particles, map<int, pair<Jet::SRC, int>>& fjMap) {
486 
487  // print debug statement
488  if (m_doDebug) {
489  cout << "SCorrelatorJetTree::AddECal(PHCompositeNode*, vector<PseudoJet>&, map<int, pair<Jet::SRC, int>>&) Adding ECal clusters..." << endl;
490  }
491 
492  // warn if jets should be charged
493  if (m_doDebug && (m_jetType != 1)) {
494  cerr << "SCorrelatorJetTree::AddECal - Warning - trying to add calorimeter clusters to charged jets!" << endl;
495  }
496 
497  // grab vertex and clusters
498  GlobalVertex* vtx = GetGlobalVertex(topNode);
499  RawClusterContainer* emClustStore = GetClusterStore(topNode, "CLUSTER_CEMC");
500 
501  // add emcal clusters if needed
502  unsigned int iCst = particles.size();
503  unsigned int nClustTot = 0;
504  unsigned int nClustAcc = 0;
505  unsigned int nClustEM = 0;
506  double eClustSum = 0.;
507 
508  // loop over em clusters
509  RawClusterContainer::ConstRange emClustRange = emClustStore -> getClusters();
510  RawClusterContainer::ConstIterator itEMClust;
511  for (itEMClust = emClustRange.first; itEMClust != emClustRange.second; ++itEMClust) {
512 
513  // grab cluster
514  const RawCluster* emClust = itEMClust -> second;
515  if (!emClust) {
516  continue;
517  } else {
518  ++nClustTot;
519  }
520 
521  // construct vertex and get 4-momentum
522  const double vX = vtx -> get_x();
523  const double vY = vtx -> get_y();
524  const double vZ = vtx -> get_z();
525 
526  CLHEP::Hep3Vector hepVecVtx = CLHEP::Hep3Vector(vX, vY, vZ);
527  CLHEP::Hep3Vector hepVecEMClust = RawClusterUtility::GetECoreVec(*emClust, hepVecVtx);
528 
529  // check if good
530  const bool isGoodECal = IsGoodECal(hepVecEMClust);
531  if (!isGoodECal) {
532  continue;
533  } else {
534  ++nClustAcc;
535  }
536 
537  // create pseudojet and add to constituent vector
538  const int emClustID = emClust -> get_id();
539  const double emClustE = hepVecEMClust.mag();
540  const double emClustPt = hepVecEMClust.perp();
541  const double emClustPhi = hepVecEMClust.getPhi();
542  const double emClustPx = emClustPt * cos(emClustPhi);
543  const double emClustPy = emClustPt * sin(emClustPhi);
544  const double emClustPz = sqrt((emClustE * emClustE) - (emClustPx * emClustPx) - (emClustPy * emClustPy));
545 
546  fastjet::PseudoJet fjCluster(emClustPx, emClustPy, emClustPz, emClustE);
547  fjCluster.set_user_index(iCst);
548  particles.push_back(fjCluster);
549 
550  // add em cluster to fastjet map
551  pair<int, pair<Jet::SRC, int>> jetEMClustPair(iCst, make_pair(Jet::SRC::CEMC_CLUSTER, emClustID));
552  fjMap.insert(jetEMClustPair);
553 
554  // fill QA histograms, increment sums and counters
555  m_hObjectQA[OBJECT::ECLUST][INFO::PT] -> Fill(fjCluster.perp());
556  m_hObjectQA[OBJECT::ECLUST][INFO::ETA] -> Fill(fjCluster.pseudorapidity());
557  m_hObjectQA[OBJECT::ECLUST][INFO::PHI] -> Fill(fjCluster.phi_std());
558  m_hObjectQA[OBJECT::ECLUST][INFO::ENE] -> Fill(fjCluster.E());
559  eClustSum += emClustE;
560  ++nClustEM;
561  ++iCst;
562  } // end em cluster loop
563 
564  // fill QA histograms
565  m_hNumObject[OBJECT::ECLUST] -> Fill(nClustEM);
566  m_hNumCstAccept[CST_TYPE::ECAL_CST][0] -> Fill(nClustTot);
567  m_hNumCstAccept[CST_TYPE::ECAL_CST][1] -> Fill(nClustAcc);
568  m_hSumCstEne[CST_TYPE::ECAL_CST] -> Fill(eClustSum);
569  return;
570 
571  } // end 'AddECal(PHCompositeNode*, vector<PseudoJet>&, map<int, pair<Jet::SRC, int>>&)'
572 
573 
574 
575  void SCorrelatorJetTree::AddHCal(PHCompositeNode* topNode, vector<PseudoJet>& particles, map<int, pair<Jet::SRC, int>>& fjMap) {
576 
577  // print debug statement
578  if (m_doDebug) {
579  cout << "SCorrelatorJetTree::AddHCal(PHCompositeNode*, vector<PseudoJet>&, map<int, pair<Jet::SRC, int>>&) Adding HCal clusters..." << endl;
580  }
581 
582  // warn if jets should be charged
583  if (m_doDebug && (m_jetType != 1)) {
584  cerr << "SCorrelatorJetTree::AddHCal - Warning - trying to add calorimeter clusters to charged jets!" << endl;
585  }
586 
587  // grab vertex and clusters
588  GlobalVertex* vtx = GetGlobalVertex(topNode);
589  RawClusterContainer* ihClustStore = GetClusterStore(topNode, "CLUSTER_HCALIN");
590  RawClusterContainer* ohClustStore = GetClusterStore(topNode, "CLUSTER_HCALOUT");
591 
592  // add emcal clusters if needed
593  unsigned int iCst = particles.size();
594  unsigned int nClustTot = 0;
595  unsigned int nClustAcc = 0;
596  unsigned int nClustH = 0;
597  double eClustSum = 0.;
598 
599  // Loop over ih clusters
600  RawClusterContainer::ConstRange ihClustRange = ihClustStore -> getClusters();
601  RawClusterContainer::ConstIterator itIHClust;
602  for (itIHClust = ihClustRange.first; itIHClust != ihClustRange.second; ++itIHClust) {
603 
604  // get ih cluster
605  const RawCluster* ihClust = itIHClust -> second;
606  if (!ihClust) {
607  continue;
608  } else {
609  ++nClustTot;
610  }
611 
612  // construct vertex and get 4-momentum
613  const double vX = vtx -> get_x();
614  const double vY = vtx -> get_y();
615  const double vZ = vtx -> get_z();
616 
617  CLHEP::Hep3Vector hepVecVtx = CLHEP::Hep3Vector(vX, vY, vZ);
618  CLHEP::Hep3Vector hepVecIHClust = RawClusterUtility::GetECoreVec(*ihClust, hepVecVtx);
619 
620  // check if good
621  const bool isGoodHCal = IsGoodHCal(hepVecIHClust);
622  if (!isGoodHCal) {
623  continue;
624  } else {
625  ++nClustAcc;
626  }
627 
628  // create pseudojet and add to constituent vector
629  const int ihClustID = ihClust -> get_id();
630  const double ihClustE = hepVecIHClust.mag();
631  const double ihClustPt = hepVecIHClust.perp();
632  const double ihClustPhi = hepVecIHClust.getPhi();
633  const double ihClustPx = ihClustPt * cos(ihClustPhi);
634  const double ihClustPy = ihClustPt * sin(ihClustPhi);
635  const double ihClustPz = sqrt((ihClustE * ihClustE) - (ihClustPx * ihClustPx) - (ihClustPy * ihClustPy));
636 
637  fastjet::PseudoJet fjCluster(ihClustPx, ihClustPy, ihClustPz, ihClustE);
638  fjCluster.set_user_index(iCst);
639  particles.push_back(fjCluster);
640 
641  // add ih cluster to fastjet map
642  pair<int, pair<Jet::SRC, int>> jetIHClustPair(iCst, make_pair(Jet::SRC::HCALIN_CLUSTER, ihClustID));
643  fjMap.insert(jetIHClustPair);
644 
645  // fill QA histograms, increment sums and counters
646  m_hObjectQA[OBJECT::HCLUST][INFO::PT] -> Fill(fjCluster.perp());
647  m_hObjectQA[OBJECT::HCLUST][INFO::ETA] -> Fill(fjCluster.pseudorapidity());
648  m_hObjectQA[OBJECT::HCLUST][INFO::PHI] -> Fill(fjCluster.phi_std());
649  m_hObjectQA[OBJECT::HCLUST][INFO::ENE] -> Fill(fjCluster.E());
650  eClustSum += ihClustE;
651  ++nClustH;
652  ++iCst;
653  } // end ih cluster loop
654 
655  // loop over oh clusters
656  RawClusterContainer::ConstRange ohClustRange = ohClustStore -> getClusters();
657  RawClusterContainer::ConstIterator itOHClust;
658  for (itOHClust = ohClustRange.first; itOHClust != ohClustRange.second; ++itOHClust) {
659 
660  // get oh cluster
661  const RawCluster* ohClust = itOHClust -> second;
662  if (!ohClust) {
663  continue;
664  } else {
665  ++nClustTot;
666  }
667 
668  // construct vertex and get 4-momentum
669  const double vX = vtx -> get_x();
670  const double vY = vtx -> get_y();
671  const double vZ = vtx -> get_z();
672 
673  CLHEP::Hep3Vector hepVecVtx = CLHEP::Hep3Vector(vX, vY, vZ);
674  CLHEP::Hep3Vector hepVecOHClust = RawClusterUtility::GetECoreVec(*ohClust, hepVecVtx);
675 
676  // check if good
677  const bool isGoodHCal = IsGoodHCal(hepVecOHClust);
678  if (!isGoodHCal) {
679  continue;
680  } else {
681  ++nClustAcc;
682  }
683 
684  // create pseudojet and add to constituent vector
685  const int ohClustID = ohClust -> get_id();
686  const double ohClustE = hepVecOHClust.mag();
687  const double ohClustPt = hepVecOHClust.perp();
688  const double ohClustPhi = hepVecOHClust.getPhi();
689  const double ohClustPx = ohClustPt * cos(ohClustPhi);
690  const double ohClustPy = ohClustPt * sin(ohClustPhi);
691  const double ohClustPz = sqrt((ohClustE * ohClustE) - (ohClustPx * ohClustPx) - (ohClustPy * ohClustPy));
692 
693  fastjet::PseudoJet fjCluster(ohClustPx, ohClustPy, ohClustPz, ohClustE);
694  fjCluster.set_user_index(iCst);
695  particles.push_back(fjCluster);
696 
697  // add oh cluster to fastjet map
698  pair<int, pair<Jet::SRC, int>> jetOHClustPair(iCst, make_pair(Jet::SRC::HCALOUT_CLUSTER, ohClustID));
699  fjMap.insert(jetOHClustPair);
700 
701  // fill QA histograms, increment sums and counters
702  m_hObjectQA[OBJECT::HCLUST][INFO::PT] -> Fill(fjCluster.perp());
703  m_hObjectQA[OBJECT::HCLUST][INFO::ETA] -> Fill(fjCluster.pseudorapidity());
704  m_hObjectQA[OBJECT::HCLUST][INFO::PHI] -> Fill(fjCluster.phi_std());
705  m_hObjectQA[OBJECT::HCLUST][INFO::ENE] -> Fill(fjCluster.E());
706  eClustSum += ohClustE;
707  ++nClustH;
708  ++iCst;
709  } // end oh cluster loop
710 
711  // fill QA histograms
712  m_hNumObject[OBJECT::HCLUST] -> Fill(nClustH);
713  m_hSumCstEne[CST_TYPE::HCAL_CST] -> Fill(eClustSum);
714  m_hNumCstAccept[CST_TYPE::HCAL_CST][0] -> Fill(nClustTot);
715  m_hNumCstAccept[CST_TYPE::HCAL_CST][1] -> Fill(nClustAcc);
716  return;
717 
718  } // end 'AddHCal(PHCompositeNode*, vector<PseudoJet>&, map<int, pair<Jet::SRC, int>>&)'
719 
720 } // end SColdQcdCorrelatorAnalysis namespace
721 
722 // end ------------------------------------------------------------------------