31 moduleLabel_(pset.getParameter<
std::
string>(
"@module_label")),
35 "ApplyDiscriminationByTrackerIsolation");
37 "ApplyDiscriminationByECALIsolation");
39 calculateWeights_ = pset.
exists(
"ApplyDiscriminationByWeightedECALIsolation") ?
40 pset.
getParameter<
bool>(
"ApplyDiscriminationByWeightedECALIsolation") :
false;
45 weightGammas_ = pset.
exists(
"WeightECALIsolation") ?
49 minPtForNoIso_ = pset.
exists(
"minTauPtForNoIso") ?
52 applyOccupancyCut_ = pset.
getParameter<
bool>(
"applyOccupancyCut");
53 maximumOccupancy_ = pset.
getParameter<uint32_t>(
"maximumOccupancy");
55 applySumPtCut_ = pset.
getParameter<
bool>(
"applySumPtCut");
56 maximumSumPt_ = pset.
getParameter<
double>(
"maximumSumPtCut");
59 "applyRelativeSumPtCut");
62 offsetRelativeSumPt_ = pset.
exists(
"relativeSumPtOffset") ?
65 storeRawOccupancy_ = pset.
exists(
"storeRawOccupancy") ?
67 storeRawSumPt_ = pset.
exists(
"storeRawSumPt") ?
69 storeRawPUsumPt_ = pset.
exists(
"storeRawPUsumPt") ?
71 storeRawFootprintCorrection_ = pset.
exists(
"storeRawFootprintCorrection") ?
72 pset.
getParameter<
bool>(
"storeRawFootprintCorrection") :
false;
73 storeRawPhotonSumPt_outsideSignalCone_ = pset.
exists(
"storeRawPhotonSumPt_outsideSignalCone") ?
74 pset.
getParameter<
bool>(
"storeRawPhotonSumPt_outsideSignalCone") :
false;
78 if ( applySumPtCut_ || applyOccupancyCut_ || applyRelativeSumPtCut_ ) {
79 if ( storeRawSumPt_ || storeRawOccupancy_ || storeRawPUsumPt_ ) {
81 <<
"A 'store raw' and a 'apply cut' option have been set to true " 82 <<
"simultaneously. These options are mutually exclusive.";
87 if ( includeGammas_ && calculateWeights_ ) {
89 <<
"Both 'ApplyDiscriminationByECALIsolation' and 'ApplyDiscriminationByWeightedECALIsolation' " 90 <<
"have been set to true. These options are mutually exclusive.";
94 int numStoreOptions = 0;
95 if ( storeRawSumPt_ ) ++numStoreOptions;
96 if ( storeRawOccupancy_ ) ++numStoreOptions;
97 if ( storeRawPUsumPt_ ) ++numStoreOptions;
98 if ( storeRawFootprintCorrection_ ) ++numStoreOptions;
99 if ( storeRawPhotonSumPt_outsideSignalCone_ ) ++numStoreOptions;
100 if ( numStoreOptions > 1 ) {
102 <<
"Multiple 'store sum pt' and/or 'store occupancy' options are set." 103 <<
" These options are mutually exclusive.";
106 if ( pset.
exists(
"customOuterCone") ) {
107 customIsoCone_ = pset.
getParameter<
double>(
"customOuterCone");
112 applyPhotonPtSumOutsideSignalConeCut_ = ( pset.
exists(
"applyPhotonPtSumOutsideSignalConeCut") ) ?
113 pset.
getParameter<
bool>(
"applyPhotonPtSumOutsideSignalConeCut") :
false;
114 if ( applyPhotonPtSumOutsideSignalConeCut_ ) {
115 maxAbsPhotonSumPt_outsideSignalCone_ = pset.
getParameter<
double>(
"maxAbsPhotonSumPt_outsideSignalCone");
116 maxRelPhotonSumPt_outsideSignalCone_ = pset.
getParameter<
double>(
"maxRelPhotonSumPt_outsideSignalCone");
119 applyFootprintCorrection_ = ( pset.
exists(
"applyFootprintCorrection") ) ?
120 pset.
getParameter<
bool>(
"applyFootprintCorrection") :
false;
121 if ( applyFootprintCorrection_ || storeRawFootprintCorrection_ ) {
123 for ( edm::VParameterSet::const_iterator cfgFootprintCorrection = cfgFootprintCorrections.begin();
124 cfgFootprintCorrection != cfgFootprintCorrections.end(); ++cfgFootprintCorrection ) {
128 footprintCorrections_.push_back(
std::move(footprintCorrection));
134 "isolationQualityCuts");
138 vertexAssociator_.reset(
141 applyDeltaBeta_ = pset.
exists(
"applyDeltaBetaCorrection") ?
142 pset.
getParameter<
bool>(
"applyDeltaBetaCorrection") :
false;
144 if ( applyDeltaBeta_ || calculateWeights_ ) {
147 std::pair<edm::ParameterSet, edm::ParameterSet> puFactorizedIsoQCuts =
152 if ( pset.
exists(
"deltaBetaPUTrackPtCutOverride") ) {
153 puFactorizedIsoQCuts.second.addParameter<
double>(
155 pset.
getParameter<
double>(
"deltaBetaPUTrackPtCutOverride"));
158 puFactorizedIsoQCuts.second.addParameter<
double>(
164 puFactorizedIsoQCuts.first));
167 puFactorizedIsoQCuts.second));
170 pfCand_token = consumes<reco::PFCandidateCollection>(pfCandSrc_);
172 vertex_token = consumes<reco::VertexCollection>(vertexSrc_);
174 "isoConeSizeForDeltaBeta");
177 deltaBetaFormula_.reset(
178 new TFormula(
"DB_corr", deltaBetaFactorFormula.c_str()));
181 applyRhoCorrection_ = pset.
exists(
"applyRhoCorrection") ?
183 if ( applyRhoCorrection_ ) {
185 rho_token=consumes<double>(rhoProducer_);
186 rhoConeSize_ = pset.
getParameter<
double>(
"rhoConeSize");
187 rhoUEOffsetCorrection_ =
190 useAllPFCands_ = pset.
exists(
"UseAllPFCandsForWeights") ?
191 pset.
getParameter<
bool>(
"UseAllPFCandsForWeights") :
false;
193 verbosity_ = ( pset.
exists(
"verbosity") ) ?
202 double discriminate(
const PFTauRef& pfTau)
const override;
204 inline double weightedSum(
const std::vector<PFCandidatePtr>& inColl_,
double eta,
double phi)
const {
206 for (
auto const & inObj_ : inColl_){
207 double sum = (inObj_->pt()*inObj_->pt())/(
deltaR2(eta,phi,inObj_->eta(),inObj_->phi()));
208 if(sum > 1.0) out *= sum;
217 std::auto_ptr<tau::RecoTauQualityCuts>
qcuts_;
248 : selection_(selection),
301 vertexAssociator_->setEvent(event);
305 if ( applyDeltaBeta_ || calculateWeights_ ) {
308 event.getByToken(pfCand_token, pfCandidates);
309 chargedPFCandidatesInEvent_.clear();
310 chargedPFCandidatesInEvent_.reserve(pfCandidates->size());
311 size_t numPFCandidates = pfCandidates->size();
312 for (
size_t i = 0;
i < numPFCandidates; ++
i ) {
314 if ( pfCandidate->
charge() != 0 ) {
315 chargedPFCandidatesInEvent_.push_back(pfCandidate);
321 event.getByToken(vertex_token, vertices);
322 size_t nVtxThisEvent = vertices->size();
323 deltaBetaFactorThisEvent_ = deltaBetaFormula_->Eval(nVtxThisEvent);
326 if ( applyRhoCorrection_ ) {
328 event.getByToken(rho_token, rhoHandle_);
329 rhoThisEvent_ = (*rhoHandle_ - rhoUEOffsetCorrection_)*
330 (3.14159)*rhoConeSize_*rhoConeSize_;
337 LogDebug(
"discriminate") <<
" tau: Pt = " << pfTau->pt() <<
", eta = " << pfTau->eta() <<
", phi = " << pfTau->phi();
338 LogDebug(
"discriminate") << *pfTau ;
341 std::vector<PFCandidatePtr> isoCharged_;
342 std::vector<PFCandidatePtr> isoNeutral_;
343 std::vector<PFCandidatePtr> isoPU_;
345 std::vector<PFCandidatePtr> chPV_;
346 isoCharged_.reserve(pfTau->isolationPFChargedHadrCands().size());
347 isoNeutral_.reserve(pfTau->isolationPFGammaCands().size());
348 isoPU_.reserve(
std::min(100UL, chargedPFCandidatesInEvent_.size()));
349 isoNeutralWeight_.reserve(pfTau->isolationPFGammaCands().size());
351 chPV_.reserve(
std::min(50UL, chargedPFCandidatesInEvent_.size()));
359 LogTrace(
"discriminate") <<
"pv: x = " << pv->position().x() <<
", y = " << pv->position().y() <<
", z = " << pv->position().z() ;
361 LogTrace(
"discriminate") <<
"pv: N/A" ;
363 if ( pfTau->leadPFChargedHadrCand().
isNonnull() ) {
364 LogTrace(
"discriminate") <<
"leadPFChargedHadron:" 365 <<
" Pt = " << pfTau->leadPFChargedHadrCand()->pt() <<
"," 366 <<
" eta = " << pfTau->leadPFChargedHadrCand()->eta() <<
"," 367 <<
" phi = " << pfTau->leadPFChargedHadrCand()->phi() ;
369 LogTrace(
"discriminate") <<
"leadPFChargedHadron: N/A" ;
377 qcuts_->setLeadTrack(pfTau->leadPFChargedHadrCand());
379 if ( applyDeltaBeta_ || calculateWeights_) {
380 pileupQcutsGeneralQCuts_->setPV(pv);
381 pileupQcutsGeneralQCuts_->setLeadTrack(pfTau->leadPFChargedHadrCand());
382 pileupQcutsPUTrackSelection_->setPV(pv);
383 pileupQcutsPUTrackSelection_->setLeadTrack(pfTau->leadPFChargedHadrCand());
387 if ( includeTracks_ ) {
388 for(
auto const &
cand : pfTau->isolationPFChargedHadrCands() ) {
389 if ( qcuts_->filterCandRef(
cand) ) {
390 LogTrace(
"discriminate") <<
"adding charged iso cand with pt " <<
cand->pt() ;
391 isoCharged_.push_back(
cand);
395 if ( includeGammas_ || calculateWeights_ ) {
396 for(
auto const &
cand : pfTau->isolationPFGammaCands() ) {
397 if ( qcuts_->filterCandRef(
cand) ) {
398 LogTrace(
"discriminate") <<
"adding neutral iso cand with pt " <<
cand->pt() ;
399 isoNeutral_.push_back(
cand);
408 if ( applyDeltaBeta_ || calculateWeights_) {
411 std::cout <<
"Initial PFCands: " << chargedPFCandidatesInEvent_.size() << std::endl;
414 std::vector<PFCandidatePtr> allPU =
415 pileupQcutsPUTrackSelection_->filterCandRefs(
416 chargedPFCandidatesInEvent_,
true);
418 std::vector<PFCandidatePtr> allNPU =
419 pileupQcutsPUTrackSelection_->filterCandRefs(
420 chargedPFCandidatesInEvent_);
421 LogTrace(
"discriminate") <<
"After track cuts: " << allPU.size() ;
424 if ( !useAllPFCands_ ) {
425 std::vector<PFCandidatePtr> cleanPU =
426 pileupQcutsGeneralQCuts_->filterCandRefs(allPU);
428 std::vector<PFCandidatePtr> cleanNPU =
429 pileupQcutsGeneralQCuts_->filterCandRefs(allNPU);
431 LogTrace(
"discriminate") <<
"After cleaning cuts: " << cleanPU.size() ;
434 DRFilter deltaBetaFilter(pfTau->p4(), 0, deltaBetaCollectionCone_);
435 for (
auto const &
cand : cleanPU ) {
436 if ( deltaBetaFilter(
cand) ) isoPU_.push_back(
cand);
439 for (
auto const &
cand : cleanNPU ) {
440 if ( deltaBetaFilter(
cand) ) chPV_.push_back(
cand);
442 LogTrace(
"discriminate") <<
"After cone cuts: " << isoPU_.size() <<
" " << chPV_.size() ;
449 if ( calculateWeights_ ) {
450 for (
auto const & isoObject : isoNeutral_ ) {
451 if ( isoObject->charge() != 0 ) {
453 isoNeutralWeight_.push_back(*isoObject);
457 double eta = isoObject->eta();
458 double phi = isoObject->phi();
459 double sumNPU = 0.5*
log(weightedSum(chPV_, eta, phi));
461 double sumPU = 0.5*
log(weightedSum(isoPU_, eta, phi));
463 if ( (sumNPU + sumPU) > 0 ) neutral.
setP4(((sumNPU)/(sumNPU + sumPU))*neutral.
p4());
465 isoNeutralWeight_.push_back(neutral);
470 if ( customIsoCone_ >= 0. ) {
471 DRFilter
filter(pfTau->p4(), 0, customIsoCone_);
472 DRFilter2 filter2(pfTau->p4(), 0, customIsoCone_);
473 std::vector<PFCandidatePtr> isoCharged_filter;
474 std::vector<PFCandidatePtr> isoNeutral_filter;
477 for(
auto const & isoObject : isoCharged_ ) {
478 if (
filter(isoObject) ) isoCharged_filter.push_back(isoObject);
480 if(!calculateWeights_){
481 for(
auto const & isoObject : isoNeutral_ ) {
482 if (
filter(isoObject) ) isoNeutral_filter.push_back(isoObject);
484 isoNeutral_ = isoNeutral_filter;
486 for(
auto const & isoObject : isoNeutralWeight_){
487 if ( filter2(isoObject) ) isoNeutralWeight_filter.push_back(isoObject);
489 isoNeutralWeight_ = isoNeutralWeight_filter;
491 isoCharged_ = isoCharged_filter;
494 bool failsOccupancyCut =
false;
495 bool failsSumPtCut =
false;
496 bool failsRelativeSumPtCut =
false;
499 int neutrals = isoNeutral_.size();
501 if ( applyDeltaBeta_ ) {
502 neutrals -= TMath::Nint(deltaBetaFactorThisEvent_*isoPU_.size());
504 if ( neutrals < 0 ) {
508 size_t nOccupants = isoCharged_.size() + neutrals;
510 failsOccupancyCut = ( nOccupants > maximumOccupancy_ );
512 double footprintCorrection_value = 0.;
513 if ( applyFootprintCorrection_ || storeRawFootprintCorrection_ ) {
514 for ( std::vector<std::unique_ptr<FootprintCorrection> >::const_iterator
footprintCorrection = footprintCorrections_.begin();
516 if ( (*footprintCorrection)->selection_(*pfTau) ) {
517 footprintCorrection_value = (*footprintCorrection)->offset_(*pfTau);
525 if ( applySumPtCut_ || applyRelativeSumPtCut_ || storeRawSumPt_ || storeRawPUsumPt_ ) {
526 double chargedPt = 0.;
527 double neutralPt = 0.;
528 double weightedNeutralPt = 0.;
529 for (
auto const & isoObject : isoCharged_ ) {
530 chargedPt += isoObject->pt();
532 if ( !calculateWeights_ ) {
533 for (
auto const & isoObject : isoNeutral_ ) {
534 neutralPt += isoObject->pt();
537 for (
auto const & isoObject : isoNeutralWeight_ ) {
538 weightedNeutralPt += isoObject.pt();
541 for (
auto const & isoObject : isoPU_ ) {
542 puPt += isoObject->pt();
544 LogTrace(
"discriminate") <<
"chargedPt = " << chargedPt ;
545 LogTrace(
"discriminate") <<
"neutralPt = " << neutralPt ;
546 LogTrace(
"discriminate") <<
"weighted neutral Pt = " << weightedNeutralPt ;
547 LogTrace(
"discriminate") <<
"puPt = " << puPt <<
" (delta-beta corr. = " << (deltaBetaFactorThisEvent_*puPt) <<
")" ;
549 if ( calculateWeights_ ) {
550 neutralPt = weightedNeutralPt;
553 if ( applyDeltaBeta_ ) {
554 neutralPt -= (deltaBetaFactorThisEvent_*puPt);
557 if ( applyFootprintCorrection_ ) {
558 neutralPt -= footprintCorrection_value;
561 if ( applyRhoCorrection_ ) {
562 neutralPt -= rhoThisEvent_;
565 if ( neutralPt < 0. ) {
569 totalPt = chargedPt + weightGammas_ * neutralPt;
570 LogTrace(
"discriminate") <<
"totalPt = " << totalPt <<
" (cut = " << maximumSumPt_ <<
")" ;
572 failsSumPtCut = (totalPt > maximumSumPt_);
575 failsRelativeSumPtCut = (totalPt > ((pfTau->pt() - offsetRelativeSumPt_)*maximumRelativeSumPt_));
578 bool failsPhotonPtSumOutsideSignalConeCut =
false;
579 double photonSumPt_outsideSignalCone = 0.;
580 if ( applyPhotonPtSumOutsideSignalConeCut_ || storeRawPhotonSumPt_outsideSignalCone_ ) {
581 const std::vector<reco::PFCandidatePtr>& signalPFGammas = pfTau->signalPFGammaCands();
582 for ( std::vector<reco::PFCandidatePtr>::const_iterator signalPFGamma = signalPFGammas.begin();
583 signalPFGamma != signalPFGammas.end(); ++signalPFGamma ) {
584 double dR =
deltaR(pfTau->eta(), pfTau->phi(), (*signalPFGamma)->eta(), (*signalPFGamma)->phi());
585 if ( dR > pfTau->signalConeSize() ) photonSumPt_outsideSignalCone += (*signalPFGamma)->pt();
587 if ( photonSumPt_outsideSignalCone > maxAbsPhotonSumPt_outsideSignalCone_ || photonSumPt_outsideSignalCone > (maxRelPhotonSumPt_outsideSignalCone_*pfTau->pt()) ) {
588 failsPhotonPtSumOutsideSignalConeCut =
true;
592 bool fails = (applyOccupancyCut_ && failsOccupancyCut) ||
593 (applySumPtCut_ && failsSumPtCut) ||
594 (applyRelativeSumPtCut_ && failsRelativeSumPtCut) ||
595 (applyPhotonPtSumOutsideSignalConeCut_ && failsPhotonPtSumOutsideSignalConeCut);
598 if (pfTau->pt() > minPtForNoIso_ && minPtForNoIso_ > 0.){
600 LogDebug(
"discriminate") <<
"tau pt = " << pfTau->pt() <<
"\t min cutoff pt = " << minPtForNoIso_ ;
604 if ( storeRawSumPt_ ) {
606 }
else if ( storeRawPUsumPt_ ) {
607 if ( applyDeltaBeta_ )
return puPt;
608 else if ( applyRhoCorrection_ )
return rhoThisEvent_;
610 }
else if ( storeRawOccupancy_ ) {
612 }
else if ( storeRawFootprintCorrection_ ) {
613 return footprintCorrection_value;
614 }
else if ( storeRawPhotonSumPt_outsideSignalCone_ ) {
615 return photonSumPt_outsideSignalCone;
617 return (fails ? 0. : 1.);
PFRecoTauDiscriminationByIsolation(const edm::ParameterSet &pset)
T getParameter(std::string const &) const
edm::EDGetTokenT< double > rho_token
bool isNonnull() const
Checks for non-null.
std::auto_ptr< tau::RecoTauQualityCuts > pileupQcutsGeneralQCuts_
bool applyPhotonPtSumOutsideSignalConeCut_
#define DEFINE_FWK_MODULE(type)
std::vector< ParameterSet > VParameterSet
void beginEvent(const edm::Event &evt, const edm::EventSetup &evtSetup) override
bool exists(std::string const ¶meterName) const
checks if a parameter exists
int charge() const final
electric charge
double weightedSum(const std::vector< PFCandidatePtr > &inColl_, double eta, double phi) const
double rhoUEOffsetCorrection_
double maximumRelativeSumPt_
bool applyRelativeSumPtCut_
~PFRecoTauDiscriminationByIsolation() override
double rhoCorrectionThisEvent_
double deltaBetaFactorThisEvent_
bool applyFootprintCorrection_
double offsetRelativeSumPt_
std::vector< std::unique_ptr< FootprintCorrection > > footprintCorrections_
edm::EDGetTokenT< reco::PFCandidateCollection > pfCand_token
auto deltaR(const T1 &t1, const T2 &t2) -> decltype(t1.eta())
std::vector< reco::PFCandidatePtr > chargedPFCandidatesInEvent_
const LorentzVector & p4() const final
four-momentum Lorentz vector
bool storeRawFootprintCorrection_
std::pair< edm::ParameterSet, edm::ParameterSet > factorizePUQCuts(const edm::ParameterSet &inputSet)
std::auto_ptr< tau::RecoTauQualityCuts > qcuts_
std::auto_ptr< tau::RecoTauQualityCuts > pileupQcutsPUTrackSelection_
double deltaBetaCollectionCone_
std::vector< reco::PFCandidate > PFCandidateCollection
collection of PFCandidates
double discriminate(const PFTauRef &pfTau) const override
ParameterSet const & getParameterSet(std::string const &) const
uint32_t maximumOccupancy_
T1 deltaR2(T1 eta1, T2 phi1, T3 eta2, T4 phi2)
double maxRelPhotonSumPt_outsideSignalCone_
Particle reconstructed by the particle flow algorithm.
std::auto_ptr< tau::RecoTauVertexAssociator > vertexAssociator_
edm::InputTag rhoProducer_
std::auto_ptr< TFormula > deltaBetaFormula_
double maxAbsPhotonSumPt_outsideSignalCone_
bool storeRawPhotonSumPt_outsideSignalCone_
edm::ParameterSet qualityCutsPSet_
void setP4(const LorentzVector &p4) final
set 4-momentum
edm::EDGetTokenT< reco::VertexCollection > vertex_token