32 moduleLabel_(
pset.getParameter<
std::
string>(
"@module_label")),
34 includeTracks_ =
pset.getParameter<
bool>(
"ApplyDiscriminationByTrackerIsolation");
35 includeGammas_ =
pset.getParameter<
bool>(
"ApplyDiscriminationByECALIsolation");
37 calculateWeights_ =
pset.getParameter<
bool>(
"ApplyDiscriminationByWeightedECALIsolation");
39 enableHGCalWorkaround_ =
pset.getParameter<
bool>(
"enableHGCalWorkaround");
44 weightGammas_ =
pset.getParameter<
double>(
"WeightECALIsolation");
47 minPtForNoIso_ =
pset.getParameter<
double>(
"minTauPtForNoIso");
49 applyOccupancyCut_ =
pset.getParameter<
bool>(
"applyOccupancyCut");
50 maximumOccupancy_ =
pset.getParameter<uint32_t>(
"maximumOccupancy");
52 applySumPtCut_ =
pset.getParameter<
bool>(
"applySumPtCut");
53 maximumSumPt_ =
pset.getParameter<
double>(
"maximumSumPtCut");
55 applyRelativeSumPtCut_ =
pset.getParameter<
bool>(
"applyRelativeSumPtCut");
56 maximumRelativeSumPt_ =
pset.getParameter<
double>(
"relativeSumPtCut");
57 offsetRelativeSumPt_ =
pset.getParameter<
double>(
"relativeSumPtOffset");
59 storeRawOccupancy_ =
pset.getParameter<
bool>(
"storeRawOccupancy");
60 storeRawSumPt_ =
pset.getParameter<
bool>(
"storeRawSumPt");
61 storeRawPUsumPt_ =
pset.getParameter<
bool>(
"storeRawPUsumPt");
62 storeRawFootprintCorrection_ =
pset.getParameter<
bool>(
"storeRawFootprintCorrection");
63 storeRawPhotonSumPt_outsideSignalCone_ =
pset.getParameter<
bool>(
"storeRawPhotonSumPt_outsideSignalCone");
67 if (applySumPtCut_ || applyOccupancyCut_ || applyRelativeSumPtCut_) {
68 if (storeRawSumPt_ || storeRawOccupancy_ || storeRawPUsumPt_) {
69 throw cms::Exception(
"BadIsoConfig") <<
"A 'store raw' and a 'apply cut' option have been set to true "
70 <<
"simultaneously. These options are mutually exclusive.";
75 if (includeGammas_ && calculateWeights_) {
77 <<
"Both 'ApplyDiscriminationByECALIsolation' and 'ApplyDiscriminationByWeightedECALIsolation' "
78 <<
"have been set to true. These options are mutually exclusive.";
82 int numStoreOptions = 0;
85 if (storeRawOccupancy_)
89 if (storeRawFootprintCorrection_)
91 if (storeRawPhotonSumPt_outsideSignalCone_)
93 if (numStoreOptions > 1) {
94 throw cms::Exception(
"BadIsoConfig") <<
"Multiple 'store sum pt' and/or 'store occupancy' options are set."
95 <<
" These options are mutually exclusive.";
98 customIsoCone_ =
pset.getParameter<
double>(
"customOuterCone");
100 applyPhotonPtSumOutsideSignalConeCut_ =
pset.getParameter<
bool>(
"applyPhotonPtSumOutsideSignalConeCut");
101 if (applyPhotonPtSumOutsideSignalConeCut_) {
102 maxAbsPhotonSumPt_outsideSignalCone_ =
pset.getParameter<
double>(
"maxAbsPhotonSumPt_outsideSignalCone");
103 maxRelPhotonSumPt_outsideSignalCone_ =
pset.getParameter<
double>(
"maxRelPhotonSumPt_outsideSignalCone");
106 applyFootprintCorrection_ =
pset.getParameter<
bool>(
"applyFootprintCorrection");
107 if (applyFootprintCorrection_ || storeRawFootprintCorrection_) {
109 for (edm::VParameterSet::const_iterator cfgFootprintCorrection = cfgFootprintCorrections.begin();
110 cfgFootprintCorrection != cfgFootprintCorrections.end();
111 ++cfgFootprintCorrection) {
115 footprintCorrections_.push_back(
std::move(footprintCorrection));
126 applyDeltaBeta_ =
pset.getParameter<
bool>(
"applyDeltaBetaCorrection");
128 if (applyDeltaBeta_ || calculateWeights_) {
131 std::pair<edm::ParameterSet, edm::ParameterSet> puFactorizedIsoQCuts =
143 puFactorizedIsoQCuts.second.addParameter<
double>(
"minTrackPt",
152 pfCand_token = consumes<edm::View<reco::Candidate> >(pfCandSrc_);
154 vertex_token = consumes<reco::VertexCollection>(vertexSrc_);
155 deltaBetaCollectionCone_ =
pset.getParameter<
double>(
"isoConeSizeForDeltaBeta");
156 std::string deltaBetaFactorFormula =
pset.getParameter<
string>(
"deltaBetaFactor");
157 deltaBetaFormula_.reset(
new TFormula(
"DB_corr", deltaBetaFactorFormula.c_str()));
160 applyRhoCorrection_ =
pset.getParameter<
bool>(
"applyRhoCorrection");
161 if (applyRhoCorrection_) {
163 rho_token = consumes<double>(rhoProducer_);
164 rhoConeSize_ =
pset.getParameter<
double>(
"rhoConeSize");
165 rhoUEOffsetCorrection_ =
pset.getParameter<
double>(
"rhoUEOffsetCorrection");
167 useAllPFCands_ =
pset.getParameter<
bool>(
"UseAllPFCandsForWeights");
169 verbosity_ =
pset.getParameter<
int>(
"verbosity");
175 double discriminate(
const PFTauRef& pfTau)
const override;
177 inline double weightedSum(
const std::vector<CandidatePtr>& inColl_,
double eta,
double phi)
const {
179 for (
auto const& inObj_ : inColl_) {
180 double sum = (inObj_->pt() * inObj_->pt()) / (
deltaR2(
eta, phi, inObj_->eta(), inObj_->phi()));
193 std::unique_ptr<tau::RecoTauQualityCuts>
qcuts_;
274 vertexAssociator_->setEvent(
event);
278 if (applyDeltaBeta_ || calculateWeights_) {
282 chargedPFCandidatesInEvent_.clear();
283 chargedPFCandidatesInEvent_.reserve(
pfCandidates->size());
285 for (
size_t i = 0;
i < numPFCandidates; ++
i) {
287 if (pfCandidate->
charge() != 0) {
288 chargedPFCandidatesInEvent_.push_back(pfCandidate);
294 event.getByToken(vertex_token,
vertices);
295 size_t nVtxThisEvent =
vertices->size();
296 deltaBetaFactorThisEvent_ = deltaBetaFormula_->Eval(nVtxThisEvent);
299 if (applyRhoCorrection_) {
301 event.getByToken(rho_token, rhoHandle_);
302 rhoThisEvent_ = (*rhoHandle_ - rhoUEOffsetCorrection_) * (3.14159) * rhoConeSize_ * rhoConeSize_;
307 LogDebug(
"discriminate") <<
" tau: Pt = " << pfTau->pt() <<
", eta = " << pfTau->eta() <<
", phi = " << pfTau->phi();
311 std::vector<CandidatePtr> isoCharged_;
312 std::vector<CandidatePtr> isoNeutral_;
313 std::vector<CandidatePtr> isoPU_;
315 std::vector<CandidatePtr> chPV_;
316 isoCharged_.
reserve(pfTau->isolationChargedHadrCands().size());
317 isoNeutral_.reserve(pfTau->isolationGammaCands().size());
318 isoPU_.reserve(
std::min(100UL, chargedPFCandidatesInEvent_.size()));
319 isoNeutralWeight_.
reserve(pfTau->isolationGammaCands().size());
321 chPV_.reserve(
std::min(50UL, chargedPFCandidatesInEvent_.size()));
328 if (
pv.isNonnull()) {
329 LogTrace(
"discriminate") <<
"pv: x = " <<
pv->position().x() <<
", y = " <<
pv->position().y()
330 <<
", z = " <<
pv->position().z();
332 LogTrace(
"discriminate") <<
"pv: N/A";
334 if (pfTau->leadChargedHadrCand().
isNonnull()) {
335 LogTrace(
"discriminate") <<
"leadPFChargedHadron:"
336 <<
" Pt = " << pfTau->leadChargedHadrCand()->pt() <<
","
337 <<
" eta = " << pfTau->leadChargedHadrCand()->eta() <<
","
338 <<
" phi = " << pfTau->leadChargedHadrCand()->phi();
340 LogTrace(
"discriminate") <<
"leadPFChargedHadron: N/A";
345 if (!(
pv.isNonnull() && pfTau->leadChargedHadrCand().
isNonnull()))
349 qcuts_->setLeadTrack(*pfTau->leadChargedHadrCand());
351 if (applyDeltaBeta_ || calculateWeights_) {
352 pileupQcutsGeneralQCuts_->setPV(
pv);
353 pileupQcutsGeneralQCuts_->setLeadTrack(*pfTau->leadChargedHadrCand());
354 pileupQcutsPUTrackSelection_->setPV(
pv);
355 pileupQcutsPUTrackSelection_->setLeadTrack(*pfTau->leadChargedHadrCand());
359 if (includeTracks_) {
360 for (
auto const&
cand : pfTau->isolationChargedHadrCands()) {
361 if (qcuts_->filterCandRef(
cand)) {
362 LogTrace(
"discriminate") <<
"adding charged iso cand with pt " <<
cand->pt();
363 isoCharged_.push_back(
cand);
367 if (includeGammas_ || calculateWeights_) {
368 for (
auto const&
cand : pfTau->isolationGammaCands()) {
369 if (qcuts_->filterCandRef(
cand)) {
370 LogTrace(
"discriminate") <<
"adding neutral iso cand with pt " <<
cand->pt();
371 isoNeutral_.push_back(
cand);
380 if (applyDeltaBeta_ || calculateWeights_) {
383 std::cout <<
"Initial PFCands: " << chargedPFCandidatesInEvent_.size() << std::endl;
386 std::vector<CandidatePtr> allPU = pileupQcutsPUTrackSelection_->filterCandRefs(chargedPFCandidatesInEvent_,
true);
388 std::vector<CandidatePtr> allNPU = pileupQcutsPUTrackSelection_->filterCandRefs(chargedPFCandidatesInEvent_);
389 LogTrace(
"discriminate") <<
"After track cuts: " << allPU.size();
392 if (!useAllPFCands_) {
393 std::vector<CandidatePtr> cleanPU = pileupQcutsGeneralQCuts_->filterCandRefs(allPU);
395 std::vector<CandidatePtr> cleanNPU = pileupQcutsGeneralQCuts_->filterCandRefs(allNPU);
397 LogTrace(
"discriminate") <<
"After cleaning cuts: " << cleanPU.size();
400 DRFilter deltaBetaFilter(pfTau->p4(), 0, deltaBetaCollectionCone_);
401 for (
auto const&
cand : cleanPU) {
402 if (deltaBetaFilter(
cand))
403 isoPU_.push_back(
cand);
406 for (
auto const&
cand : cleanNPU) {
407 if (deltaBetaFilter(
cand))
408 chPV_.push_back(
cand);
410 LogTrace(
"discriminate") <<
"After cone cuts: " << isoPU_.size() <<
" " << chPV_.size();
417 if (calculateWeights_) {
418 for (
auto const& isoObject : isoNeutral_) {
419 if (isoObject->charge() != 0) {
425 double eta = isoObject->eta();
426 double phi = isoObject->phi();
427 double sumNPU = 0.5 *
log(weightedSum(chPV_,
eta, phi));
429 double sumPU = 0.5 *
log(weightedSum(isoPU_,
eta, phi));
431 if ((sumNPU + sumPU) > 0)
432 neutral.
setP4(((sumNPU) / (sumNPU + sumPU)) * neutral.
p4());
439 if (customIsoCone_ >= 0.) {
440 DRFilter
filter(pfTau->p4(), 0, customIsoCone_);
441 DRFilter2 filter2(pfTau->p4(), 0, customIsoCone_);
442 std::vector<CandidatePtr> isoCharged_filter;
443 std::vector<CandidatePtr> isoNeutral_filter;
446 for (
auto const& isoObject : isoCharged_) {
450 if (!calculateWeights_) {
451 for (
auto const& isoObject : isoNeutral_) {
453 isoNeutral_filter.push_back(isoObject);
455 isoNeutral_ = isoNeutral_filter;
457 for (
auto const& isoObject : isoNeutralWeight_) {
458 if (filter2(isoObject))
459 isoNeutralWeight_filter.
push_back(isoObject);
461 isoNeutralWeight_ = isoNeutralWeight_filter;
463 isoCharged_ = isoCharged_filter;
466 bool failsOccupancyCut =
false;
467 bool failsSumPtCut =
false;
468 bool failsRelativeSumPtCut =
false;
471 int neutrals = isoNeutral_.
size();
473 if (applyDeltaBeta_) {
474 neutrals -= TMath::Nint(deltaBetaFactorThisEvent_ * isoPU_.size());
480 size_t nOccupants = isoCharged_.size() + neutrals;
482 failsOccupancyCut = (nOccupants > maximumOccupancy_);
484 double footprintCorrection_value = 0.;
485 if (applyFootprintCorrection_ || storeRawFootprintCorrection_) {
486 for (std::vector<std::unique_ptr<FootprintCorrection> >::const_iterator footprintCorrection =
487 footprintCorrections_.begin();
488 footprintCorrection != footprintCorrections_.end();
489 ++footprintCorrection) {
490 if ((*footprintCorrection)->selection_(*pfTau)) {
491 footprintCorrection_value = (*footprintCorrection)->offset_(*pfTau);
499 if (applySumPtCut_ || applyRelativeSumPtCut_ || storeRawSumPt_ || storeRawPUsumPt_) {
500 double chargedPt = 0.;
501 double neutralPt = 0.;
502 double weightedNeutralPt = 0.;
503 for (
auto const& isoObject : isoCharged_) {
507 if (enableHGCalWorkaround_) {
508 double trackPt = (isoObject->bestTrack()) ? isoObject->bestTrack()->pt() : 0.;
509 double pfCandPt = isoObject->pt();
512 neutralPt += pfCandPt -
trackPt;
514 chargedPt += isoObject->pt();
517 chargedPt += isoObject->pt();
521 if (!calculateWeights_) {
522 for (
auto const& isoObject : isoNeutral_) {
523 neutralPt += isoObject->pt();
526 for (
auto const& isoObject : isoNeutralWeight_) {
527 weightedNeutralPt += isoObject.pt();
530 for (
auto const& isoObject : isoPU_) {
531 puPt += isoObject->pt();
533 LogTrace(
"discriminate") <<
"chargedPt = " << chargedPt;
534 LogTrace(
"discriminate") <<
"neutralPt = " << neutralPt;
535 LogTrace(
"discriminate") <<
"weighted neutral Pt = " << weightedNeutralPt;
536 LogTrace(
"discriminate") <<
"puPt = " << puPt <<
" (delta-beta corr. = " << (deltaBetaFactorThisEvent_ * puPt)
539 if (calculateWeights_) {
540 neutralPt = weightedNeutralPt;
543 if (applyDeltaBeta_) {
544 neutralPt -= (deltaBetaFactorThisEvent_ * puPt);
547 if (applyFootprintCorrection_) {
548 neutralPt -= footprintCorrection_value;
551 if (applyRhoCorrection_) {
552 neutralPt -= rhoThisEvent_;
555 if (neutralPt < 0.) {
559 totalPt = chargedPt + weightGammas_ * neutralPt;
560 LogTrace(
"discriminate") <<
"totalPt = " << totalPt <<
" (cut = " << maximumSumPt_ <<
")";
562 failsSumPtCut = (totalPt > maximumSumPt_);
565 failsRelativeSumPtCut = (totalPt > ((pfTau->pt() - offsetRelativeSumPt_) * maximumRelativeSumPt_));
568 bool failsPhotonPtSumOutsideSignalConeCut =
false;
569 double photonSumPt_outsideSignalCone = 0.;
570 if (applyPhotonPtSumOutsideSignalConeCut_ || storeRawPhotonSumPt_outsideSignalCone_) {
571 const std::vector<reco::CandidatePtr>& signalGammas = pfTau->signalGammaCands();
572 for (std::vector<reco::CandidatePtr>::const_iterator signalGamma = signalGammas.begin();
573 signalGamma != signalGammas.end();
575 double dR =
deltaR(pfTau->eta(), pfTau->phi(), (*signalGamma)->eta(), (*signalGamma)->phi());
576 if (
dR > pfTau->signalConeSize())
577 photonSumPt_outsideSignalCone += (*signalGamma)->pt();
579 if (photonSumPt_outsideSignalCone > maxAbsPhotonSumPt_outsideSignalCone_ ||
580 photonSumPt_outsideSignalCone > (maxRelPhotonSumPt_outsideSignalCone_ * pfTau->pt())) {
581 failsPhotonPtSumOutsideSignalConeCut =
true;
585 bool fails = (applyOccupancyCut_ && failsOccupancyCut) || (applySumPtCut_ && failsSumPtCut) ||
586 (applyRelativeSumPtCut_ && failsRelativeSumPtCut) ||
587 (applyPhotonPtSumOutsideSignalConeCut_ && failsPhotonPtSumOutsideSignalConeCut);
589 if (pfTau->pt() > minPtForNoIso_ && minPtForNoIso_ > 0.) {
591 LogDebug(
"discriminate") <<
"tau pt = " << pfTau->pt() <<
"\t min cutoff pt = " << minPtForNoIso_;
595 if (storeRawSumPt_) {
597 }
else if (storeRawPUsumPt_) {
600 else if (applyRhoCorrection_)
601 return rhoThisEvent_;
604 }
else if (storeRawOccupancy_) {
606 }
else if (storeRawFootprintCorrection_) {
607 return footprintCorrection_value;
608 }
else if (storeRawPhotonSumPt_outsideSignalCone_) {
609 return photonSumPt_outsideSignalCone;
611 return (fails ? 0. : 1.);
618 desc.
add<
bool>(
"storeRawFootprintCorrection",
false);
620 desc.
add<
bool>(
"storeRawOccupancy",
false);
621 desc.
add<
double>(
"maximumSumPtCut", 6.0);
625 pset_signalQualityCuts.
add<
double>(
"maxDeltaZ", 0.4);
626 pset_signalQualityCuts.
add<
double>(
"minTrackPt", 0.5);
627 pset_signalQualityCuts.
add<
double>(
"minTrackVertexWeight", -1.0);
628 pset_signalQualityCuts.
add<
double>(
"maxTrackChi2", 100.0);
629 pset_signalQualityCuts.
add<
unsigned int>(
"minTrackPixelHits", 0);
630 pset_signalQualityCuts.
add<
double>(
"minGammaEt", 1.0);
631 pset_signalQualityCuts.
add<
unsigned int>(
"minTrackHits", 3);
632 pset_signalQualityCuts.
add<
double>(
"minNeutralHadronEt", 30.0);
633 pset_signalQualityCuts.
add<
double>(
"maxTransverseImpactParameter", 0.1);
634 pset_signalQualityCuts.
addOptional<
bool>(
"useTracksInsteadOfPFHadrons");
637 pset_vxAssocQualityCuts.
add<
double>(
"minTrackPt", 0.5);
638 pset_vxAssocQualityCuts.add<
double>(
"minTrackVertexWeight", -1.0);
639 pset_vxAssocQualityCuts.add<
double>(
"maxTrackChi2", 100.0);
640 pset_vxAssocQualityCuts.add<
unsigned int>(
"minTrackPixelHits", 0);
641 pset_vxAssocQualityCuts.add<
double>(
"minGammaEt", 1.0);
642 pset_vxAssocQualityCuts.add<
unsigned int>(
"minTrackHits", 3);
643 pset_vxAssocQualityCuts.add<
double>(
"maxTransverseImpactParameter", 0.1);
644 pset_vxAssocQualityCuts.addOptional<
bool>(
"useTracksInsteadOfPFHadrons");
647 pset_isolationQualityCuts.
add<
double>(
"maxDeltaZ", 0.2);
648 pset_isolationQualityCuts.add<
double>(
"minTrackPt", 1.0);
649 pset_isolationQualityCuts.add<
double>(
"minTrackVertexWeight", -1.0);
650 pset_isolationQualityCuts.add<
double>(
"maxTrackChi2", 100.0);
651 pset_isolationQualityCuts.add<
unsigned int>(
"minTrackPixelHits", 0);
652 pset_isolationQualityCuts.add<
double>(
"minGammaEt", 1.5);
653 pset_isolationQualityCuts.add<
unsigned int>(
"minTrackHits", 8);
654 pset_isolationQualityCuts.add<
double>(
"maxTransverseImpactParameter", 0.03);
655 pset_isolationQualityCuts.addOptional<
bool>(
"useTracksInsteadOfPFHadrons");
661 pset_qualityCuts.
add<
std::string>(
"leadingTrkOrPFCandOption",
"leadPFCand");
662 pset_qualityCuts.add<
std::string>(
"pvFindingAlgo",
"closestInDeltaZ");
664 pset_qualityCuts.add<
bool>(
"vertexTrackFiltering",
false);
665 pset_qualityCuts.add<
bool>(
"recoverLeadingTrk",
false);
670 desc.
add<
double>(
"minTauPtForNoIso", -99.0);
671 desc.
add<
double>(
"maxAbsPhotonSumPt_outsideSignalCone", 1000000000.0);
673 desc.
add<
bool>(
"applySumPtCut",
false);
674 desc.
add<
double>(
"rhoConeSize", 0.5);
675 desc.
add<
bool>(
"ApplyDiscriminationByTrackerIsolation",
true);
676 desc.
add<
bool>(
"storeRawPhotonSumPt_outsideSignalCone",
false);
678 desc.
add<
bool>(
"enableHGCalWorkaround",
false);
684 desc.
addVPSet(
"footprintCorrections", vpsd1, {});
688 desc.
add<
bool>(
"applyFootprintCorrection",
false);
689 desc.
add<
bool>(
"UseAllPFCandsForWeights",
false);
690 desc.
add<
double>(
"relativeSumPtCut", 0.0);
696 psd1.
add<
double>(
"cut");
712 psd1.
add<
double>(
"cut");
725 psd1.
add<
double>(
"cut");
732 desc.
add<
unsigned int>(
"maximumOccupancy", 0);
733 desc.
add<
int>(
"verbosity", 0);
735 desc.
add<
bool>(
"applyOccupancyCut",
true);
736 desc.
add<
bool>(
"applyDeltaBetaCorrection",
false);
737 desc.
add<
bool>(
"applyRelativeSumPtCut",
false);
738 desc.
add<
bool>(
"storeRawPUsumPt",
false);
739 desc.
add<
bool>(
"applyPhotonPtSumOutsideSignalConeCut",
false);
740 desc.
add<
bool>(
"deltaBetaPUTrackPtCutOverride",
false);
741 desc.
add<
bool>(
"ApplyDiscriminationByWeightedECALIsolation",
false);
742 desc.
add<
bool>(
"storeRawSumPt",
false);
743 desc.
add<
bool>(
"ApplyDiscriminationByECALIsolation",
true);
744 desc.
add<
bool>(
"applyRhoCorrection",
false);
746 desc.
add<
double>(
"WeightECALIsolation", 1.0);
747 desc.
add<
double>(
"rhoUEOffsetCorrection", 1.0);
748 desc.
add<
double>(
"maxRelPhotonSumPt_outsideSignalCone", 0.1);
749 desc.
add<
double>(
"deltaBetaPUTrackPtCutOverride_val", -1.5);
750 desc.
add<
double>(
"isoConeSizeForDeltaBeta", 0.5);
751 desc.
add<
double>(
"relativeSumPtOffset", 0.0);
752 desc.
add<
double>(
"customOuterCone", -1.0);
755 descriptions.
add(
"pfRecoTauDiscriminationByIsolation", desc);