All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
CustomPFOAnalysis.cpp
Go to the documentation of this file.
2 
4 {
5  m_pfoVector.clear();
6  m_pfoTargetVector.clear();
7 
8  m_nPfosTotal = 0;
10  m_nPfosPhotons = 0;
11  m_nPfosTracks = 0;
12  m_pfoEnergyTotal = 0.f;
14  m_pfoEnergyPhotons = 0.f;
15 
16  m_pfoEnergyTracks = 0.f;
17  m_pfoECalToEmEnergy = 0.f;
19  m_pfoHCalToEmEnergy = 0.f;
21  m_pfoMuonToEnergy = 0.f;
22  m_pfoOtherEnergy = 0.f;
23  m_pfoMassTotal = 0.f;
24 
25  m_pfoEnergies.clear();
26  m_pfoPx.clear();
27  m_pfoPy.clear();
28  m_pfoPz.clear();
29  m_pfoCosTheta.clear();
30 
31  m_pfoTargetEnergies.clear();
32  m_pfoTargetPx.clear();
33  m_pfoTargetPy.clear();
34  m_pfoTargetPz.clear();
35  m_pfoTargetCosTheta.clear();
36 
37  m_pfoPdgCodes.clear();
38  m_pfoTargetPdgCodes.clear();
39 
44 
49 
50  m_mcEnergyENu = 0.f;
51  m_mcEnergyFwd = 0.f;
52  m_eQQ = -99.f;
53  m_eQ1 = -99.f;
54  m_eQ2 = -99.f;
55  m_costQQ = -99.f;
56  m_costQ1 = -99.f;
57  m_costQ2 = -99.f;
58  m_mQQ = -99.f;
59  m_thrust = -99.f;
60  m_qPdg = -99;
61 
62 }
63 
64 //------------------------------------------------------------------------------------------------------------------------------------------
65 
67 {
68  // Extract reconstructed pfo collection
69  try
70  {
71 
72  ReconstructedParticleVec recos = jet_content->second;
73 
74  for (unsigned int i = 0, nElements = recos.size(); i < nElements; ++i)
75  {
76  const EVENT::ReconstructedParticle *pReconstructedParticle = recos[i];
77 
78  if (NULL == pReconstructedParticle)
79  throw EVENT::Exception("Collection type mismatch");
80 
81  m_pfoVector.push_back(pReconstructedParticle);
82  }
83  }
84  catch (...)
85  {
86  streamlog_out(ERROR) << "Could not extract input particle collection: " << m_inputPfoCollection << std::endl;
87  }
88 
89  // Extract mc pfo collection
90  MCParticleList preliminary_mcPfoCandidates, mcPfoCandidates;
91 
92  try
93  {
94 
95  MCParticleVec mcs = jet_content->first;
96 
97  for (unsigned int i = 0, nElements = mcs.size(); i < nElements; ++i)
98  {
99  MCParticle* pMCParticle = mcs[i];
100 
101  if (pMCParticle == NULL) {
102  throw EVENT::Exception("Collection type mismatch");
103  }
104 
105  // if (!pMCParticle->getParents().empty()) {
106  // continue;
107  // }
108 
109  this->ApplyPfoSelectionRules(pMCParticle, preliminary_mcPfoCandidates);
110  }
111  streamlog_out(DEBUG) << "N leftover mcs: " << preliminary_mcPfoCandidates.size() << " N input: " << mcs.size() << std::endl;
112 
113 
114  std::set<EVENT::MCParticle*>::iterator it;
115  for (it = preliminary_mcPfoCandidates.begin(); it != preliminary_mcPfoCandidates.end(); it++) {
116  if ( hasSomeParentsInMCList( *it, preliminary_mcPfoCandidates ) ) {
117  continue;
118  }
119  mcPfoCandidates.insert(*it);
120  }
121  }
122  catch (...)
123  {
124  streamlog_out(WARNING) << "Could not extract mc particle collection " << m_mcParticleCollection << std::endl;
125  }
126 
127  m_pfoTargetVector.insert(m_pfoTargetVector.begin(), mcPfoCandidates.begin(), mcPfoCandidates.end());
129 }
130 
131 //------------------------------------------------------------------------------------------------------------------------------------------
132 
133 bool TJjetsPFOAnalysisProcessor::hasSomeParentsInMCList(EVENT::MCParticle *pMCParticle, MCParticleList &mcs) const {
134  // Check if the MCparticle has some parents/grandparents/grandgrandparents/... in the given MCParticleList
135  bool found_parent = false;
136  MCParticleVec pMCParticleParents (pMCParticle->getParents());
137  MCParticleVec mcs_vector (mcs.begin(), mcs.end());
138 
139  if ( ! areDisjointVectors( pMCParticleParents, mcs_vector ) ) {
140  found_parent = true;
141  } else { // Check next higher generation
142  for (unsigned int i=0; i<pMCParticleParents.size(); i++) {
143  if ( hasSomeParentsInMCList(pMCParticleParents[i], mcs) ) {
144  found_parent = true;
145  }
146  }
147  }
148 
149  return found_parent;
150 }
151 
152 
153 //------------------------------------------------------------------------------------------------------------------------------------------
154 
155 void TJjetsPFOAnalysisProcessor::ApplyPfoSelectionRules(EVENT::MCParticle *pMCParticle, MCParticleList &mcPfoCandidates) const
156 {
157  const float innerRadius(std::sqrt(pMCParticle->getVertex()[0] * pMCParticle->getVertex()[0] + pMCParticle->getVertex()[1] * pMCParticle->getVertex()[1] + pMCParticle->getVertex()[2] * pMCParticle->getVertex()[2]));
158  const float outerRadius(std::sqrt(pMCParticle->getEndpoint()[0] * pMCParticle->getEndpoint()[0] + pMCParticle->getEndpoint()[1] * pMCParticle->getEndpoint()[1] + pMCParticle->getEndpoint()[2] * pMCParticle->getEndpoint()[2]));
159  const float momentum(std::sqrt(pMCParticle->getMomentum()[0] * pMCParticle->getMomentum()[0] + pMCParticle->getMomentum()[1] * pMCParticle->getMomentum()[1] + pMCParticle->getMomentum()[2] * pMCParticle->getMomentum()[2]));
160 
161  if ((mcPfoCandidates.find(pMCParticle) == mcPfoCandidates.end()) &&
162  (outerRadius > m_mcPfoSelectionRadius) && (innerRadius <= m_mcPfoSelectionRadius) && (momentum > m_mcPfoSelectionMomentum) &&
163  !((pMCParticle->getPDG() == 2212 || pMCParticle->getPDG() == 2112) && (pMCParticle->getEnergy() < m_mcPfoSelectionLowEnergyNPCutOff)))
164  {
165  mcPfoCandidates.insert(pMCParticle);
166  }
167  // else
168  // {
169  // const EVENT::MCParticleVec &daughterVector(pMCParticle->getDaughters());
170  //
171  // for (EVENT::MCParticleVec::const_iterator iter = daughterVector.begin(), iterEnd = daughterVector.end(); iter != iterEnd; ++iter)
172  // {
173  // this->ApplyPfoSelectionRules(*iter, mcPfoCandidates);
174  // }
175  // }
176 }
177 
178 //------------------------------------------------------------------------------------------------------------------------------------------
179 
181 {
182  MCParticleVector mcQuarkVector;
183 
184  try
185  {
186 
187  MCParticleVec mcs = jet_content->first;
188  MCParticleList mcs_list = MCParticleList(mcs.begin(),mcs.end());
189 
190  for (unsigned int i = 0, nElements = mcs.size(); i < nElements; ++i)
191  {
192  EVENT::MCParticle *pMCParticle = mcs[i];
193 
194  if (NULL == pMCParticle)
195  throw EVENT::Exception("Collection type mismatch");
196 
197  const int absPdgCode(std::abs(pMCParticle->getPDG()));
198 
199  // By default, the primary quarks are the ones without parents included in the jet-mc
201  {
202  if ((absPdgCode >= 1) && (absPdgCode <= 6) && !hasSomeParentsInMCList(pMCParticle, mcs_list) )// pMCParticle->getParents().empty()) //!!!!!!!!!!!!!!!!!!!!!!!!!!!
203  mcQuarkVector.push_back(pMCParticle);
204  }
205  else
206  {
207  // For MC files generated in the SLIC environment, the primary quarks have parents; the mother should be the Z-boson
208  if ((absPdgCode >= 1) && (absPdgCode <= 6))
209  {
210  if ((pMCParticle->getParents().size() == 1) && ((pMCParticle->getParents())[0]->getPDG() == 23))
211  mcQuarkVector.push_back(pMCParticle);
212  }
213  }
214  }
215  }
216  catch (...)
217  {
218  streamlog_out(WARNING) << "Could not extract mc quark information" << std::endl;
219  }
220 
221  if (!mcQuarkVector.empty())
222  {
223  m_qPdg = std::abs(mcQuarkVector[0]->getPDG());
224  float energyTot(0.f);
225  float costTot(0.f);
226 
227  for (unsigned int i = 0; i < mcQuarkVector.size(); ++i)
228  {
229  const float px(mcQuarkVector[i]->getMomentum()[0]);
230  const float py(mcQuarkVector[i]->getMomentum()[1]);
231  const float pz(mcQuarkVector[i]->getMomentum()[2]);
232  const float energy(mcQuarkVector[i]->getEnergy());
233  const float p(std::sqrt(px * px + py * py + pz * pz));
234  const float cost(std::fabs(pz) / p);
235  energyTot += energy;
236  costTot += cost * energy;
237  }
238 
239  m_thrust = costTot / energyTot;
240 
241  const float pQ1x = mcQuarkVector[0]->getMomentum()[0];
242  const float pQ1y = mcQuarkVector[0]->getMomentum()[1];
243  const float pQ1z = mcQuarkVector[0]->getMomentum()[2];
244 
245  const float pQ1[3] = {pQ1x, pQ1y, pQ1z};
246  const TLorentzVector q1(pQ1[0], pQ1[1], pQ1[2], mcQuarkVector[0]->getEnergy());
247  const float pQ1Tot(std::sqrt(pQ1[0] * pQ1[0] + pQ1[1] * pQ1[1] + pQ1[2] * pQ1[2]));
248 
249  if (std::fabs(pQ1Tot) > std::numeric_limits<float>::epsilon())
250  m_costQ1 = pQ1[2] / pQ1Tot;
251 
252  m_eQ1 = pQ1Tot;
253 
254  if (mcQuarkVector.size() == 2)
255  {
256 
257  const float pQ2x = mcQuarkVector[1]->getMomentum()[0];
258  const float pQ2y = mcQuarkVector[1]->getMomentum()[1];
259  const float pQ2z = mcQuarkVector[1]->getMomentum()[2];
260 
261  const float pQ2[3] = {pQ2x, pQ2y, pQ2z};
262  const float pQQ[3] = {pQ1[0] + pQ2[0], pQ1[1] + pQ2[1], pQ1[2] + pQ2[2]};
263 
264  const TLorentzVector q2(pQ2[0], pQ2[1], pQ2[2], mcQuarkVector[1]->getEnergy());
265  const TLorentzVector qq = q1 + q2;
266 
267  m_mQQ = qq.M();
268  m_eQQ = mcQuarkVector[0]->getEnergy() + mcQuarkVector[1]->getEnergy();
269 
270  const float pQ2Tot(std::sqrt(pQ2[0] * pQ2[0] + pQ2[1] * pQ2[1] + pQ2[2] * pQ2[2]));
271  const float pQQTot(std::sqrt(pQQ[0] * pQQ[0] + pQQ[1] * pQQ[1] + pQQ[2] * pQQ[2]));
272 
273  if (std::fabs(pQQTot) > std::numeric_limits<float>::epsilon())
274  m_costQQ = pQQ[2] / pQQTot;
275 
276  if (std::fabs(pQ2Tot) > std::numeric_limits<float>::epsilon())
277  m_costQ2 = pQ2[2] / pQ2Tot;
278 
279  m_eQ2 = pQ2Tot;
280 
281  if (m_printing > 0)
282  {
283  std::cout << " eQQ = " << m_eQQ << std::endl
284  << " eQ1 = " << m_eQ1 << std::endl
285  << " eQ2 = " << m_eQ2 << std::endl
286  << " costQQ = " << m_costQQ << std::endl
287  << " costQ1 = " << m_costQ1 << std::endl
288  << " costQ2 = " << m_costQ2 << std::endl
289  << " mQQ = " << m_mQQ << std::endl
290  << " Thrust = " << m_thrust << std::endl
291  << " QPDG = " << m_qPdg << std::endl;
292  }
293  }
294  }
295 }
296 
297 //------------------------------------------------------------------------------------------------------------------------------------------
298 
300 {
301  float momTot[3] = {0.f, 0.f, 0.f};
302 
303  // Extract quantities relating to reconstructed pfos
304  for (ParticleVector::const_iterator iter = m_pfoVector.begin(), iterEnd = m_pfoVector.end(); iter != iterEnd; ++iter)
305  {
306  const EVENT::ReconstructedParticle *pPfo = *iter;
307 
308  ++m_nPfosTotal;
309  m_pfoEnergyTotal += pPfo->getEnergy();
310  m_pfoPdgCodes.push_back(pPfo->getType());
311  m_pfoEnergies.push_back(pPfo->getEnergy());
312 
313  m_pfoPx.push_back(pPfo->getMomentum()[0]);
314  m_pfoPy.push_back(pPfo->getMomentum()[1]);
315  m_pfoPz.push_back(pPfo->getMomentum()[2]);
316 
317  const float momentum(std::sqrt(pPfo->getMomentum()[0] * pPfo->getMomentum()[0] + pPfo->getMomentum()[1] * pPfo->getMomentum()[1] + pPfo->getMomentum()[2] * pPfo->getMomentum()[2]));
318  const float cosTheta((momentum > std::numeric_limits<float>::epsilon()) ? pPfo->getMomentum()[2] / momentum : -999.f);
319  m_pfoCosTheta.push_back(cosTheta);
320 
321  if (!pPfo->getTracks().empty())
322  {
323  // Charged pfos
324  ++m_nPfosTracks;
325  m_pfoEnergyTracks += pPfo->getEnergy();
326  }
327  else
328  {
329  // Neutral pfos
330  float cellEnergySum(0.f);
331  const EVENT::ClusterVec &clusterVec(pPfo->getClusters());
332 
333  for (EVENT::ClusterVec::const_iterator clustIter = clusterVec.begin(), clustIterEnd = clusterVec.end(); clustIter != clustIterEnd; ++clustIter)
334  {
335  const EVENT::CalorimeterHitVec &calorimeterHitVec((*clustIter)->getCalorimeterHits());
336 
337  for (EVENT::CalorimeterHitVec::const_iterator hitIter = calorimeterHitVec.begin(), hitIterEnd = calorimeterHitVec.end(); hitIter != hitIterEnd; ++hitIter)
338  {
339  const EVENT::CalorimeterHit *pCalorimeterHit = *hitIter;
340  cellEnergySum += pCalorimeterHit->getEnergy();
341  }
342  }
343 
344  const float correctionFactor((cellEnergySum < std::numeric_limits<float>::epsilon()) ? 0.f : pPfo->getEnergy() / cellEnergySum);
345 
346  if (22 == pPfo->getType())
347  {
348  // Photons
349  ++m_nPfosPhotons;
350  m_pfoEnergyPhotons += pPfo->getEnergy();
351 
352  for (EVENT::ClusterVec::const_iterator clustIter = clusterVec.begin(), clustIterEnd = clusterVec.end(); clustIter != clustIterEnd; ++clustIter)
353  {
354  const EVENT::CalorimeterHitVec &calorimeterHitVec((*clustIter)->getCalorimeterHits());
355 
356  for (EVENT::CalorimeterHitVec::const_iterator hitIter = calorimeterHitVec.begin(), hitIterEnd = calorimeterHitVec.end(); hitIter != hitIterEnd; ++hitIter)
357  {
358  const EVENT::CalorimeterHit *pCalorimeterHit = *hitIter;
359 
360  const float hitEnergy(correctionFactor * pCalorimeterHit->getEnergy());
361  const CHT cht(pCalorimeterHit->getType());
362 
363  if (cht.is(CHT::ecal))
364  {
365  m_pfoECalToEmEnergy += hitEnergy;
366  }
367  else if (cht.is(CHT::hcal))
368  {
369  m_pfoHCalToEmEnergy += hitEnergy;
370  }
371  else if (cht.is(CHT::muon))
372  {
373  m_pfoMuonToEnergy += hitEnergy;
374  }
375  else
376  {
377  m_pfoOtherEnergy += hitEnergy;
378  }
379  }
380  }
381  }
382  else
383  {
384  // Neutral hadrons
386  m_pfoEnergyNeutralHadrons += pPfo->getEnergy();
387 
388  for (EVENT::ClusterVec::const_iterator clustIter = clusterVec.begin(), clustIterEnd = clusterVec.end(); clustIter != clustIterEnd; ++clustIter)
389  {
390  const EVENT::CalorimeterHitVec &calorimeterHitVec((*clustIter)->getCalorimeterHits());
391 
392  for (EVENT::CalorimeterHitVec::const_iterator hitIter = calorimeterHitVec.begin(), hitIterEnd = calorimeterHitVec.end(); hitIter != hitIterEnd; ++hitIter)
393  {
394  const EVENT::CalorimeterHit *pCalorimeterHit = *hitIter;
395 
396  const float hitEnergy(correctionFactor * pCalorimeterHit->getEnergy());
397  const CHT cht(pCalorimeterHit->getType());
398 
399  if (cht.is(CHT::ecal))
400  {
401  m_pfoECalToHadEnergy += hitEnergy;
402  }
403  else if (cht.is(CHT::hcal))
404  {
405  m_pfoHCalToHadEnergy += hitEnergy;
406  }
407  else if (cht.is(CHT::muon))
408  {
409  m_pfoMuonToEnergy += hitEnergy;
410  }
411  else
412  {
413  m_pfoOtherEnergy += hitEnergy;
414  }
415  }
416  }
417  }
418  }
419 
420  momTot[0] += pPfo->getMomentum()[0];
421  momTot[1] += pPfo->getMomentum()[1];
422  momTot[2] += pPfo->getMomentum()[2];
423 
424  if (m_printing > 0)
425  {
426  std::cout << " RECO PFO, pdg: " << pPfo->getType() << ", E: " << pPfo->getEnergy() << ", nTracks: " << pPfo->getTracks().size()
427  << ", nClusters: " << pPfo->getClusters().size() << ", charge: " << pPfo->getCharge() << std::endl;
428  }
429  }
430 
431  m_pfoMassTotal = std::sqrt(m_pfoEnergyTotal * m_pfoEnergyTotal - momTot[0] * momTot[0] - momTot[1] * momTot[1] - momTot[2] * momTot[2]);
432 
433  // Extract quantities relating to pfo targets
434  for (MCParticleVector::const_iterator iter = m_pfoTargetVector.begin(), iterEnd = m_pfoTargetVector.end(); iter != iterEnd; ++iter)
435  {
436  const EVENT::MCParticle *pMCParticle = *iter;
437 
439  m_pfoTargetsEnergyTotal += pMCParticle->getEnergy();
440  m_pfoTargetPdgCodes.push_back(pMCParticle->getPDG());
441  m_pfoTargetEnergies.push_back(pMCParticle->getEnergy());
442 
443  m_pfoTargetPx.push_back(pMCParticle->getMomentum()[0]);
444  m_pfoTargetPy.push_back(pMCParticle->getMomentum()[1]);
445  m_pfoTargetPz.push_back(pMCParticle->getMomentum()[2]);
446 
447  const float momentum(std::sqrt(pMCParticle->getMomentum()[0] * pMCParticle->getMomentum()[0] + pMCParticle->getMomentum()[1] * pMCParticle->getMomentum()[1] + pMCParticle->getMomentum()[2] * pMCParticle->getMomentum()[2]));
448  const float cosTheta((momentum > std::numeric_limits<float>::epsilon()) ? pMCParticle->getMomentum()[2] / momentum : -999.f);
449  m_pfoTargetCosTheta.push_back(cosTheta);
450 
451  if (std::fabs(cosTheta) > 0.98f)
452  m_mcEnergyFwd += pMCParticle->getEnergy();
453 
454  if ((std::abs(pMCParticle->getPDG()) == 12) || (std::abs(pMCParticle->getPDG()) == 14) || (std::abs(pMCParticle->getPDG()) == 16))
455  m_mcEnergyENu += pMCParticle->getEnergy();
456 
457  if (22 == pMCParticle->getPDG())
458  {
460  m_pfoTargetsEnergyPhotons += pMCParticle->getEnergy();
461  }
462  else if ((11 == std::abs(pMCParticle->getPDG())) || (13 == std::abs(pMCParticle->getPDG())) || (211 == std::abs(pMCParticle->getPDG()))) // TODO, more options here?
463  {
465  m_pfoTargetsEnergyTracks += pMCParticle->getEnergy();
466  }
467  else
468  {
470  m_pfoTargetsEnergyNeutralHadrons += pMCParticle->getEnergy();
471  }
472  }
473 
474  if (m_printing > 0)
475  {
476  std::cout << " EVENT : " << m_nEvt << std::endl
477  << " NPFOs : " << m_nPfosTotal << " (" << m_nPfosTracks << " + " << m_nPfosPhotons + m_nPfosNeutralHadrons << ")" << std::endl
478  << " RECONSTRUCTED ENERGY : " << m_pfoEnergyTotal << std::endl
479  << " RECO ENERGY + eNu : " << m_pfoEnergyTotal + m_mcEnergyENu << std::endl;
480  }
481 
482  // Fill basic histograms
484 
485  if ((m_qPdg >= 1) && (m_qPdg <= 3) && (m_thrust <= 0.7f))
487 }
488 
489 //------------------------------------------------------------------------------------------------------------------------------------------
490 
491 bool TJjetsPFOAnalysisProcessor::SortPfoTargetsByEnergy(const EVENT::MCParticle *const pLhs, const EVENT::MCParticle *const pRhs)
492 {
493  return (pLhs->getEnergy() > pRhs->getEnergy());
494 }
void ApplyPfoSelectionRules(EVENT::MCParticle *pMCParticle, MCParticleList &mcPfoCandidates) const
Apply pfo selection rules, starting with root particles.
bool areDisjointVectors(const std::vector< Object > &v1, const std::vector< Object > &v2) const
Definition: VectorHelper.h:5
std::set< EVENT::MCParticle * > MCParticleList
void MakeQuarkVariables(JetContentPair *jet_content)
Make quark variables.
std::vector< const MCParticle * > MCParticleVector
void PerformPfoAnalysis()
Perform pfo analysis.
void ExtractCollections(JetContentPair *jet_content)
Extract lcio collections.
void Clear()
PFOAnalysis stuff.
static bool SortPfoTargetsByEnergy(const EVENT::MCParticle *const pLhs, const EVENT::MCParticle *const pRhs)
Sort mc pfo targets by decreasing energy.
bool hasSomeParentsInMCList(EVENT::MCParticle *pMCParticle, MCParticleList &mcs) const
std::pair< MCParticleVec, ReconstructedParticleVec > JetContentPair
Definitions using typedef.