72 ReconstructedParticleVec recos = jet_content->second;
74 for (
unsigned int i = 0, nElements = recos.size(); i < nElements; ++i)
76 const EVENT::ReconstructedParticle *pReconstructedParticle = recos[i];
78 if (NULL == pReconstructedParticle)
79 throw EVENT::Exception(
"Collection type mismatch");
86 streamlog_out(ERROR) <<
"Could not extract input particle collection: " <<
m_inputPfoCollection << std::endl;
95 MCParticleVec mcs = jet_content->first;
97 for (
unsigned int i = 0, nElements = mcs.size(); i < nElements; ++i)
99 MCParticle* pMCParticle = mcs[i];
101 if (pMCParticle == NULL) {
102 throw EVENT::Exception(
"Collection type mismatch");
111 streamlog_out(DEBUG) <<
"N leftover mcs: " << preliminary_mcPfoCandidates.size() <<
" N input: " << mcs.size() << std::endl;
114 std::set<EVENT::MCParticle*>::iterator it;
115 for (it = preliminary_mcPfoCandidates.begin(); it != preliminary_mcPfoCandidates.end(); it++) {
119 mcPfoCandidates.insert(*it);
124 streamlog_out(WARNING) <<
"Could not extract mc particle collection " <<
m_mcParticleCollection << std::endl;
135 bool found_parent =
false;
136 MCParticleVec pMCParticleParents (pMCParticle->getParents());
137 MCParticleVec mcs_vector (mcs.begin(), mcs.end());
142 for (
unsigned int i=0; i<pMCParticleParents.size(); i++) {
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]));
161 if ((mcPfoCandidates.find(pMCParticle) == mcPfoCandidates.end()) &&
165 mcPfoCandidates.insert(pMCParticle);
187 MCParticleVec mcs = jet_content->first;
190 for (
unsigned int i = 0, nElements = mcs.size(); i < nElements; ++i)
192 EVENT::MCParticle *pMCParticle = mcs[i];
194 if (NULL == pMCParticle)
195 throw EVENT::Exception(
"Collection type mismatch");
197 const int absPdgCode(std::abs(pMCParticle->getPDG()));
203 mcQuarkVector.push_back(pMCParticle);
208 if ((absPdgCode >= 1) && (absPdgCode <= 6))
210 if ((pMCParticle->getParents().size() == 1) && ((pMCParticle->getParents())[0]->getPDG() == 23))
211 mcQuarkVector.push_back(pMCParticle);
218 streamlog_out(WARNING) <<
"Could not extract mc quark information" << std::endl;
221 if (!mcQuarkVector.empty())
223 m_qPdg = std::abs(mcQuarkVector[0]->getPDG());
224 float energyTot(0.f);
227 for (
unsigned int i = 0; i < mcQuarkVector.size(); ++i)
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);
236 costTot += cost * energy;
241 const float pQ1x = mcQuarkVector[0]->getMomentum()[0];
242 const float pQ1y = mcQuarkVector[0]->getMomentum()[1];
243 const float pQ1z = mcQuarkVector[0]->getMomentum()[2];
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]));
249 if (std::fabs(pQ1Tot) > std::numeric_limits<float>::epsilon())
254 if (mcQuarkVector.size() == 2)
257 const float pQ2x = mcQuarkVector[1]->getMomentum()[0];
258 const float pQ2y = mcQuarkVector[1]->getMomentum()[1];
259 const float pQ2z = mcQuarkVector[1]->getMomentum()[2];
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]};
264 const TLorentzVector q2(pQ2[0], pQ2[1], pQ2[2], mcQuarkVector[1]->getEnergy());
265 const TLorentzVector qq = q1 + q2;
268 m_eQQ = mcQuarkVector[0]->getEnergy() + mcQuarkVector[1]->getEnergy();
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]));
273 if (std::fabs(pQQTot) > std::numeric_limits<float>::epsilon())
276 if (std::fabs(pQ2Tot) > std::numeric_limits<float>::epsilon())
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;
301 float momTot[3] = {0.f, 0.f, 0.f};
304 for (ParticleVector::const_iterator iter =
m_pfoVector.begin(), iterEnd =
m_pfoVector.end(); iter != iterEnd; ++iter)
306 const EVENT::ReconstructedParticle *pPfo = *iter;
313 m_pfoPx.push_back(pPfo->getMomentum()[0]);
314 m_pfoPy.push_back(pPfo->getMomentum()[1]);
315 m_pfoPz.push_back(pPfo->getMomentum()[2]);
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);
321 if (!pPfo->getTracks().empty())
330 float cellEnergySum(0.f);
331 const EVENT::ClusterVec &clusterVec(pPfo->getClusters());
333 for (EVENT::ClusterVec::const_iterator clustIter = clusterVec.begin(), clustIterEnd = clusterVec.end(); clustIter != clustIterEnd; ++clustIter)
335 const EVENT::CalorimeterHitVec &calorimeterHitVec((*clustIter)->getCalorimeterHits());
337 for (EVENT::CalorimeterHitVec::const_iterator hitIter = calorimeterHitVec.begin(), hitIterEnd = calorimeterHitVec.end(); hitIter != hitIterEnd; ++hitIter)
339 const EVENT::CalorimeterHit *pCalorimeterHit = *hitIter;
340 cellEnergySum += pCalorimeterHit->getEnergy();
344 const float correctionFactor((cellEnergySum < std::numeric_limits<float>::epsilon()) ? 0.f : pPfo->getEnergy() / cellEnergySum);
346 if (22 == pPfo->getType())
352 for (EVENT::ClusterVec::const_iterator clustIter = clusterVec.begin(), clustIterEnd = clusterVec.end(); clustIter != clustIterEnd; ++clustIter)
354 const EVENT::CalorimeterHitVec &calorimeterHitVec((*clustIter)->getCalorimeterHits());
356 for (EVENT::CalorimeterHitVec::const_iterator hitIter = calorimeterHitVec.begin(), hitIterEnd = calorimeterHitVec.end(); hitIter != hitIterEnd; ++hitIter)
358 const EVENT::CalorimeterHit *pCalorimeterHit = *hitIter;
360 const float hitEnergy(correctionFactor * pCalorimeterHit->getEnergy());
361 const CHT cht(pCalorimeterHit->getType());
363 if (cht.is(CHT::ecal))
367 else if (cht.is(CHT::hcal))
388 for (EVENT::ClusterVec::const_iterator clustIter = clusterVec.begin(), clustIterEnd = clusterVec.end(); clustIter != clustIterEnd; ++clustIter)
390 const EVENT::CalorimeterHitVec &calorimeterHitVec((*clustIter)->getCalorimeterHits());
392 for (EVENT::CalorimeterHitVec::const_iterator hitIter = calorimeterHitVec.begin(), hitIterEnd = calorimeterHitVec.end(); hitIter != hitIterEnd; ++hitIter)
394 const EVENT::CalorimeterHit *pCalorimeterHit = *hitIter;
396 const float hitEnergy(correctionFactor * pCalorimeterHit->getEnergy());
397 const CHT cht(pCalorimeterHit->getType());
399 if (cht.is(CHT::ecal))
403 else if (cht.is(CHT::hcal))
420 momTot[0] += pPfo->getMomentum()[0];
421 momTot[1] += pPfo->getMomentum()[1];
422 momTot[2] += pPfo->getMomentum()[2];
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;
436 const EVENT::MCParticle *pMCParticle = *iter;
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);
451 if (std::fabs(cosTheta) > 0.98f)
454 if ((std::abs(pMCParticle->getPDG()) == 12) || (std::abs(pMCParticle->getPDG()) == 14) || (std::abs(pMCParticle->getPDG()) == 16))
457 if (22 == pMCParticle->getPDG())
462 else if ((11 == std::abs(pMCParticle->getPDG())) || (13 == std::abs(pMCParticle->getPDG())) || (211 == std::abs(pMCParticle->getPDG())))
476 std::cout <<
" EVENT : " <<
m_nEvt << std::endl
493 return (pLhs->getEnergy() > pRhs->getEnergy());
FloatVector m_pfoTargetEnergies
FloatVector m_pfoEnergies
float m_mcPfoSelectionMomentum
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
std::set< EVENT::MCParticle * > MCParticleList
void MakeQuarkVariables(JetContentPair *jet_content)
Make quark variables.
float m_pfoHCalToHadEnergy
float m_mcPfoSelectionLowEnergyNPCutOff
float m_pfoTargetsEnergyTracks
float m_pfoEnergyNeutralHadrons
float m_pfoECalToEmEnergy
int m_nPfoTargetsNeutralHadrons
int m_lookForQuarksWithMotherZ
FloatVector m_pfoCosTheta
float m_pfoHCalToEmEnergy
FloatVector m_pfoTargetPx
MCParticleVector m_pfoTargetVector
float m_pfoTargetsEnergyNeutralHadrons
std::string m_mcParticleCollection
FloatVector m_pfoTargetPz
std::vector< const MCParticle * > MCParticleVector
IntVector m_pfoTargetPdgCodes
float m_pfoTargetsEnergyPhotons
void PerformPfoAnalysis()
Perform pfo analysis.
void ExtractCollections(JetContentPair *jet_content)
Extract lcio collections.
void Clear()
PFOAnalysis stuff.
int m_nPfosNeutralHadrons
TH1F * m_hPfoEnergySumL7A
float m_pfoTargetsEnergyTotal
static bool SortPfoTargetsByEnergy(const EVENT::MCParticle *const pLhs, const EVENT::MCParticle *const pRhs)
Sort mc pfo targets by decreasing energy.
ParticleVector m_pfoVector
FloatVector m_pfoTargetCosTheta
float m_pfoECalToHadEnergy
bool hasSomeParentsInMCList(EVENT::MCParticle *pMCParticle, MCParticleList &mcs) const
std::pair< MCParticleVec, ReconstructedParticleVec > JetContentPair
Definitions using typedef.
FloatVector m_pfoTargetPy
float m_mcPfoSelectionRadius
std::string m_inputPfoCollection