113 std::cout <<
"<PFRecoTauEnergyAlgorithmPlugin::operator()>:" << std::endl;
114 std::cout <<
"tau: Pt = " << tau.
pt() <<
", eta = " << tau.
eta() <<
", phi = " << tau.
phi() <<
" (En = " << tau.
energy() <<
", decayMode = " << tau.
decayMode() <<
")" << std::endl;
118 std::vector<reco::CandidatePtr> addNeutrals;
120 const auto& jetConstituents = tau.
jetRef()->daughterPtrVector();
121 for (
const auto& jetConstituent : jetConstituents) {
123 int jetConstituentPdgId =
std::abs(jetConstituent->pdgId());
125 (jetConstituentPdgId == 22 && jetConstituent->et() >
minGammaEt_ )) )
continue;
127 bool isSignalCand =
false;
129 for (
const auto& signalCand : signalCands) {
130 if (
isPtrEqual(jetConstituent, signalCand) ) isSignalCand =
true;
132 if ( isSignalCand )
continue;
134 double dR =
deltaR(jetConstituent->p4(), tau.
p4());
137 else if ( jetConstituentPdgId == 22 ) dRadd =
dRaddPhoton_;
139 addNeutrals.push_back(jetConstituent);
140 addNeutralsSumP4 += jetConstituent->p4();
144 std::cout <<
"addNeutralsSumP4: En = " << addNeutralsSumP4.energy() << std::endl;
147 unsigned numNonPFCandTracks = 0;
148 double nonPFCandTracksSumP = 0.;
149 double nonPFCandTracksSumPerr2 = 0.;
151 for ( std::vector<PFRecoTauChargedHadron>::const_iterator
chargedHadron = chargedHadrons.begin();
154 ++numNonPFCandTracks;
156 if ( chargedHadronTrack !=
nullptr ) {
157 nonPFCandTracksSumP += chargedHadronTrack->
p();
158 nonPFCandTracksSumPerr2 += getTrackPerr2(*chargedHadronTrack);
163 std::cout <<
"nonPFCandTracksSumP = " << nonPFCandTracksSumP <<
" +/- " <<
sqrt(nonPFCandTracksSumPerr2)
164 <<
" (numNonPFCandTracks = " << numNonPFCandTracks <<
")" << std::endl;
167 if ( numNonPFCandTracks == 0 ) {
172 std::cout <<
"easy case: all tracks are associated to PFCandidates --> leaving tau momentum untouched." << std::endl;
174 updateTauP4(tau, 1., addNeutralsSumP4);
183 if ( nonPFCandTracksSumP < addNeutralsSumP4.energy() ) {
184 double scaleFactor = 1. - nonPFCandTracksSumP/addNeutralsSumP4.energy();
185 if ( !(scaleFactor >= 0. && scaleFactor <= 1.) ) {
187 <<
"Failed to compute tau energy --> killing tau candidate !!" << std::endl;
192 std::cout <<
"case (2): addNeutralsSumEn > nonPFCandTracksSumP --> adjusting tau momentum." << std::endl;
194 updateTauP4(tau, scaleFactor, addNeutralsSumP4);
200 std::vector<reco::CandidatePtr> mergedNeutrals;
202 for ( std::vector<PFRecoTauChargedHadron>::const_iterator
chargedHadron = chargedHadrons.begin();
205 const std::vector<reco::CandidatePtr>& neutralPFCands =
chargedHadron->getNeutralPFCandidates();
206 for ( std::vector<reco::CandidatePtr>::const_iterator neutralPFCand = neutralPFCands.begin();
207 neutralPFCand != neutralPFCands.end(); ++neutralPFCand ) {
208 mergedNeutrals.push_back(*neutralPFCand);
209 mergedNeutralsSumP4 += (*neutralPFCand)->p4();
214 std::cout <<
"mergedNeutralsSumP4: En = " << mergedNeutralsSumP4.energy() << std::endl;
219 if ( nonPFCandTracksSumP < (addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy()) ) {
220 double scaleFactor = ((addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy()) - nonPFCandTracksSumP)/mergedNeutralsSumP4.energy();
221 if ( !(scaleFactor >= 0. && scaleFactor <= 1.) ) {
223 <<
"Failed to compute tau energy --> killing tau candidate !!" << std::endl;
228 size_t numChargedHadrons = chargedHadrons.size();
229 for (
size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron ) {
230 const PFRecoTauChargedHadron&
chargedHadron = chargedHadrons[iChargedHadron];
231 if ( !chargedHadron.getNeutralPFCandidates().empty() ) {
232 PFRecoTauChargedHadron chargedHadron_modified =
chargedHadron;
235 diffP4 += (chargedHadron.p4() - chargedHadron_modified.p4());
239 std::cout <<
"case (3): (addNeutralsSumEn + mergedNeutralsSumEn) > nonPFCandTracksSumP --> adjusting tau momentum." << std::endl;
241 updateTauP4(tau, -1., diffP4);
246 unsigned numChargedHadronNeutrals = 0;
247 std::vector<reco::CandidatePtr> chargedHadronNeutrals;
249 for ( std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();
252 ++numChargedHadronNeutrals;
253 chargedHadronNeutrals.push_back(chargedHadron->getChargedPFCandidate());
254 chargedHadronNeutralsSumP4 += chargedHadron->getChargedPFCandidate()->p4();
258 std::cout <<
"chargedHadronNeutralsSumP4: En = " << chargedHadronNeutralsSumP4.energy()
259 <<
" (numChargedHadronNeutrals = " << numChargedHadronNeutrals <<
")" << std::endl;
264 if ( nonPFCandTracksSumP < (addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy() + chargedHadronNeutralsSumP4.energy()) ) {
265 double scaleFactor = ((addNeutralsSumP4.energy() + mergedNeutralsSumP4.energy() + chargedHadronNeutralsSumP4.energy()) - nonPFCandTracksSumP)/chargedHadronNeutralsSumP4.energy();
266 if ( !(scaleFactor >= 0. && scaleFactor <= 1.) ) {
268 <<
"Failed to compute tau energy --> killing tau candidate !!" << std::endl;
273 size_t numChargedHadrons = chargedHadrons.size();
274 for (
size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron ) {
275 const PFRecoTauChargedHadron& chargedHadron = chargedHadrons[iChargedHadron];
277 PFRecoTauChargedHadron chargedHadron_modified =
chargedHadron;
278 chargedHadron_modified.neutralPFCandidates_.clear();
279 const CandidatePtr& chargedPFCand = chargedHadron.getChargedPFCandidate();
280 double chargedHadronPx_modified = scaleFactor*chargedPFCand->px();
281 double chargedHadronPy_modified = scaleFactor*chargedPFCand->py();
282 double chargedHadronPz_modified = scaleFactor*chargedPFCand->pz();
284 chargedHadron_modified.setP4(chargedHadronP4_modified);
286 diffP4 += (chargedHadron.p4() - chargedHadron_modified.p4());
290 std::cout <<
"case (4): (addNeutralsSumEn + mergedNeutralsSumEn + chargedHadronNeutralsSumEn) > nonPFCandTracksSumP --> adjusting momenta of tau and chargedHadrons." << std::endl;
292 updateTauP4(tau, -1., diffP4);
295 double allTracksSumP = 0.;
296 double allTracksSumPerr2 = 0.;
298 for ( std::vector<PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();
302 if ( chargedHadronTrack !=
nullptr ) {
303 allTracksSumP += chargedHadronTrack->
p();
304 allTracksSumPerr2 += getTrackPerr2(*chargedHadronTrack);
306 edm::LogInfo(
"PFRecoTauEnergyAlgorithmPlugin::operator()")
307 <<
"PFRecoTauChargedHadron has no associated reco::Track !!";
309 chargedHadron->print();
315 std::cout <<
"allTracksSumP = " << allTracksSumP <<
" +/- " <<
sqrt(allTracksSumPerr2) << std::endl;
317 double allNeutralsSumEn = 0.;
319 for (
const auto& signalCand : signalCands) {
321 std::cout <<
"Candidate #" << signalCand.id() <<
":" << signalCand.key() <<
":" 322 <<
" Pt = " << (signalCand)->
pt() <<
", eta = " << (signalCand)->
eta() <<
", phi = " << (signalCand)->
phi() << std::endl;
328 <<
" ECAL = " << (pfCand)->ecalEnergy() <<
"," 329 <<
" HCAL = " << (pfCand)->hcalEnergy() <<
"," 330 <<
" HO = " << (pfCand)->hoEnergy() << std::endl;
333 if (
edm::isFinite(pfCand->ecalEnergy()) ) allNeutralsSumEn += pfCand->ecalEnergy();
334 if (
edm::isFinite(pfCand->hcalEnergy()) ) allNeutralsSumEn += pfCand->hcalEnergy();
335 if (
edm::isFinite(pfCand->hoEnergy()) ) allNeutralsSumEn += pfCand->hoEnergy();
338 allNeutralsSumEn += addNeutralsSumP4.energy();
339 if ( allNeutralsSumEn < 0. ) allNeutralsSumEn = 0.;
341 std::cout <<
"allNeutralsSumEn = " << allNeutralsSumEn << std::endl;
343 if ( allNeutralsSumEn > allTracksSumP ) {
345 size_t numChargedHadrons = chargedHadrons.size();
346 for (
size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron ) {
347 const PFRecoTauChargedHadron& chargedHadron = chargedHadrons[iChargedHadron];
349 PFRecoTauChargedHadron chargedHadron_modified =
chargedHadron;
350 chargedHadron_modified.neutralPFCandidates_.clear();
351 chargedHadron_modified.setP4(chargedHadron.getChargedPFCandidate()->p4());
353 std::cout <<
"chargedHadron #" << iChargedHadron <<
": changing En = " << chargedHadron.energy() <<
" to " << chargedHadron_modified.energy() << std::endl;
357 PFRecoTauChargedHadron chargedHadron_modified =
chargedHadron;
358 chargedHadron_modified.neutralPFCandidates_.clear();
361 if ( chTrack !=
nullptr ) {
362 double chargedHadronPx_modified = chTrack->
px();
363 double chargedHadronPy_modified = chTrack->
py();
364 double chargedHadronPz_modified = chTrack->
pz();
368 <<
"PFRecoTauChargedHadron has no associated reco::Track !!" << std::endl;
370 chargedHadron.print();
373 chargedHadron_modified.setP4(chargedHadronP4_modified);
375 std::cout <<
"chargedHadron #" << iChargedHadron <<
": changing En = " << chargedHadron.energy() <<
" to " << chargedHadron_modified.energy() << std::endl;
380 double scaleFactor = allNeutralsSumEn/tau.
energy();
382 std::cout <<
"case (5): allNeutralsSumEn > allTracksSumP --> adjusting momenta of tau and chargedHadrons." << std::endl;
384 updateTauP4(tau, scaleFactor - 1., tau.
p4());
389 size_t numChargedHadrons = chargedHadrons.size();
390 for (
size_t iChargedHadron = 0; iChargedHadron < numChargedHadrons; ++iChargedHadron ) {
391 const PFRecoTauChargedHadron& chargedHadron = chargedHadrons[iChargedHadron];
393 PFRecoTauChargedHadron chargedHadron_modified =
chargedHadron;
394 chargedHadron_modified.neutralPFCandidates_.clear();
397 if ( chargedHadronTrack !=
nullptr ) {
398 double trackP = chargedHadronTrack->
p();
399 double trackPerr2 = getTrackPerr2(*chargedHadronTrack);
401 std::cout <<
"trackP = " << trackP <<
" +/- " <<
sqrt(trackPerr2) << std::endl;
405 double trackP_modified =
406 (trackP*(allTracksSumPerr2 - trackPerr2)
407 + trackPerr2*(allNeutralsSumEn - (allTracksSumP - trackP)))/
412 if ( trackP_modified < 1.
e-1 ) trackP_modified = 1.e-1;
414 std::cout <<
"trackP (modified) = " << trackP_modified << std::endl;
416 double scaleFactor = trackP_modified/trackP;
417 if ( !(scaleFactor >= 0. && scaleFactor <= 1.) ) {
419 <<
"Failed to compute tau energy --> killing tau candidate !!" << std::endl;
423 double chargedHadronPx_modified = scaleFactor*chargedHadronTrack->
px();
424 double chargedHadronPy_modified = scaleFactor*chargedHadronTrack->
py();
425 double chargedHadronPz_modified = scaleFactor*chargedHadronTrack->
pz();
428 edm::LogInfo(
"PFRecoTauEnergyAlgorithmPlugin::operator()")
429 <<
"PFRecoTauChargedHadron has no associated reco::Track !!";
431 chargedHadron.print();
434 chargedHadron_modified.setP4(chargedHadronP4_modified);
436 std::cout <<
"chargedHadron #" << iChargedHadron <<
": changing En = " << chargedHadron.energy() <<
" to " << chargedHadron_modified.energy() << std::endl;
441 double scaleFactor = allNeutralsSumEn/tau.
energy();
443 std::cout <<
"case (6): allNeutralsSumEn < allTracksSumP --> adjusting momenta of tau and chargedHadrons." << std::endl;
445 updateTauP4(tau, scaleFactor - 1., tau.
p4());
452 std::cout <<
"case (7): allNeutralsSumEn < allTracksSumP not compatible with tau decay mode hypothesis --> killing tau candidate." << std::endl;
463 std::cout <<
"undefined case: you should never come here !!" << std::endl;
double p() const
momentum vector magnitude
double eta() const final
momentum pseudorapidity
double pt() const final
transverse momentum
double px() const
x coordinate of momentum vector
constexpr bool isFinite(T x)
double minNeutralHadronEt_
hadronicDecayMode decayMode() const
const std::vector< RecoTauPiZero > & signalPiZeroCandidates() const
Retrieve the association of signal region gamma candidates into candidate PiZeros.
const JetBaseRef & jetRef() const
double energy() const final
energy
Abs< T >::type abs(const T &t)
const LorentzVector & p4() const final
four-momentum Lorentz vector
double pz() const
z coordinate of momentum vector
constexpr auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
bool isPtrEqual(const edm::Ptr< Base > &b, const edm::Ptr< Der > &d)
edm::Ptr< Candidate > CandidatePtr
persistent reference to an object in a collection of Candidate objects
double dRaddNeutralHadron_
math::XYZTLorentzVector LorentzVector
Lorentz vector.
const reco::Track * getTrackFromChargedHadron(const reco::PFRecoTauChargedHadron &chargedHadron)
std::vector< PFRecoTauChargedHadron > & signalTauChargedHadronCandidatesRestricted()
const std::vector< reco::CandidatePtr > & signalCands() const
Candidates in signal region.
void setChargedHadronP4(reco::PFRecoTauChargedHadron &chargedHadron, double scaleFactor_neutralPFCands=1.0)
double phi() const final
momentum azimuthal angle
double py() const
y coordinate of momentum vector
reco::Candidate::LorentzVector compChargedHadronP4fromPxPyPz(double, double, double)