RecoTau/src/RecoTauQualityCuts.cc

0001 #include "RecoTauTag/RecoTau/interface/RecoTauQualityCuts.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0005 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0006 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0007
0008 namespace reco::tau {
0009
0010   namespace {
0011     const reco::Track* getTrack(const Candidate& cand) {
0012       const PFCandidate* pfCandPtr = dynamic_cast<const PFCandidate*>(&cand);
0013       if (pfCandPtr) {
0014         // Get the KF track if it exists.  Otherwise, see if PFCandidate has a GSF track.
0015         if (pfCandPtr->trackRef().isNonnull())
0016           return pfCandPtr->trackRef().get();
0017         else if (pfCandPtr->gsfTrackRef().isNonnull())
0018           return pfCandPtr->gsfTrackRef().get();
0019         else
0020           return nullptr;
0021       }
0022
0023       const pat::PackedCandidate* packedCand = dynamic_cast<const pat::PackedCandidate*>(&cand);
0024       if (packedCand && packedCand->hasTrackDetails())
0025         return &packedCand->pseudoTrack();
0026
0027       return nullptr;
0028     }
0029
0030     const reco::TrackRef getTrackRef(const Candidate& cand) {
0031       const PFCandidate* pfCandPtr = dynamic_cast<const PFCandidate*>(&cand);
0032       if (pfCandPtr)
0033         return pfCandPtr->trackRef();
0034
0035       return reco::TrackRef();
0036     }
0037
0038     const reco::TrackBaseRef getGsfTrackRef(const Candidate& cand) {
0039       const PFCandidate* pfCandPtr = dynamic_cast<const PFCandidate*>(&cand);
0040       if (pfCandPtr) {
0041         return reco::TrackBaseRef(pfCandPtr->gsfTrackRef());
0042       }
0043       return reco::TrackBaseRef();
0044     }
0045
0046     // Translate GsfTrackRef to TrackBaseRef
0047     template <typename T>
0048     reco::TrackBaseRef convertRef(const T& ref) {
0049       return reco::TrackBaseRef(ref);
0050     }
0051   }  // namespace
0052
0053   // Quality cut implementations
0054   namespace qcuts {
0055     bool minPackedCandVertexWeight(const pat::PackedCandidate& pCand, const reco::VertexRef* pv, double cut) {
0056       if (pv->isNull()) {
0057         edm::LogError("QCutsNoPrimaryVertex") << "Primary vertex Ref in "
0058                                               << "RecoTauQualityCuts is invalid. - minPackedCandVertexWeight";
0059         return false;
0060       }
0061       //there is some low granular information on track weight in the vertex available with packed cands
0062       double weight = -9.9;
0063       if (pCand.vertexRef().isNonnull() && pCand.vertexRef().key() == pv->key()) {
0064         int quality = pCand.pvAssociationQuality();
0065         if (quality == pat::PackedCandidate::UsedInFitTight)
0066           weight = 0.6;  //0.6 as proxy for weight above 0.5
0067         else if (quality == pat::PackedCandidate::UsedInFitLoose)
0068           weight = 0.1;  //0.6 as proxy for weight below 0.5
0069       }
0070       LogDebug("TauQCuts") << " packedCand: Pt = " << pCand.pt() << ", eta = " << pCand.eta()
0071                            << ", phi = " << pCand.phi();
0072       LogDebug("TauQCuts") << " vertex: x = " << (*pv)->position().x() << ", y = " << (*pv)->position().y()
0073                            << ", z = " << (*pv)->position().z();
0074       LogDebug("TauQCuts") << "--> trackWeight from packedCand = " << weight << " (cut = " << cut << ")";
0075       return (weight >= cut);
0076     }
0077   }  // namespace qcuts
0078
0079   RecoTauQualityCuts::RecoTauQualityCuts(const edm::ParameterSet& qcuts) {
0080     // Setup all of our predicates
0081     CandQCutFuncCollection chargedHadronCuts;
0082     CandQCutFuncCollection gammaCuts;
0083     CandQCutFuncCollection neutralHadronCuts;
0084
0085     // Make sure there are no extra passed options
0086     std::set<std::string> passedOptionSet;
0087     std::vector<std::string> passedOptions = qcuts.getParameterNames();
0088
0089     for (auto const& option : passedOptions) {
0090       passedOptionSet.insert(option);
0091     }
0092
0093     unsigned int nCuts = 0;
0094     auto getDouble = [&qcuts, &passedOptionSet, &nCuts](const std::string& name) {
0095       if (qcuts.exists(name)) {
0096         ++nCuts;
0097         passedOptionSet.erase(name);
0098         return qcuts.getParameter<double>(name);
0099       }
0100       return -1.0;
0101     };
0102     auto getUint = [&qcuts, &passedOptionSet, &nCuts](const std::string& name) -> unsigned int {
0103       if (qcuts.exists(name)) {
0104         ++nCuts;
0105         passedOptionSet.erase(name);
0106         return qcuts.getParameter<unsigned int>(name);
0107       }
0108       return 0;
0109     };
0110
0111     // Build all the QCuts for tracks
0112     minTrackPt_ = getDouble("minTrackPt");
0113     maxTrackChi2_ = getDouble("maxTrackChi2");
0114     minTrackPixelHits_ = getUint("minTrackPixelHits");
0115     minTrackHits_ = getUint("minTrackHits");
0116     maxTransverseImpactParameter_ = getDouble("maxTransverseImpactParameter");
0117     maxDeltaZ_ = getDouble("maxDeltaZ");
0118     maxDeltaZToLeadTrack_ = getDouble("maxDeltaZToLeadTrack");
0119     // Require tracks to contribute a minimum weight to the associated vertex.
0120     minTrackVertexWeight_ = getDouble("minTrackVertexWeight");
0121
0122     // Use bit-wise & to avoid conditional code
0123     checkHitPattern_ = (minTrackHits_ > 0) || (minTrackPixelHits_ > 0);
0124     checkPV_ = (maxTransverseImpactParameter_ >= 0) || (maxDeltaZ_ >= 0) || (maxDeltaZToLeadTrack_ >= 0) ||
0125                (minTrackVertexWeight_ >= 0);
0126
0127     // Build the QCuts for gammas
0128     minGammaEt_ = getDouble("minGammaEt");
0129
0130     // Build QCuts for netural hadrons
0131     minNeutralHadronEt_ = getDouble("minNeutralHadronEt");
0132
0133     // Check if there are any remaining unparsed QCuts
0134     if (!passedOptionSet.empty()) {
0135       std::string unParsedOptions;
0136       bool thereIsABadParameter = false;
0137       for (auto const& option : passedOptionSet) {
0138         // Workaround for HLT - TODO FIXME
0139         if (option == "useTracksInsteadOfPFHadrons") {
0140           // Crash if true - no one should have this option enabled.
0141           if (qcuts.getParameter<bool>("useTracksInsteadOfPFHadrons")) {
0142             throw cms::Exception("DontUseTracksInQcuts") << "The obsolete exception useTracksInsteadOfPFHadrons "
0143                                                          << "is set to true in the quality cut config." << std::endl;
0144           }
0145           continue;
0146         }
0147
0148         // If we get to this point, there is a real unknown parameter
0149         thereIsABadParameter = true;
0150
0151         unParsedOptions += option;
0152         unParsedOptions += "\n";
0153       }
0154       if (thereIsABadParameter) {
0155         throw cms::Exception("BadQualityCutConfig") << " The PSet passed to the RecoTauQualityCuts class had"
0156                                                     << " the following unrecognized options: " << std::endl
0157                                                     << unParsedOptions;
0158       }
0159     }
0160
0161     // Make sure there are at least some quality cuts
0162     if (!nCuts) {
0163       throw cms::Exception("BadQualityCutConfig") << " No options were passed to the quality cut class!" << std::endl;
0164     }
0165   }
0166
0167   std::pair<edm::ParameterSet, edm::ParameterSet> factorizePUQCuts(const edm::ParameterSet& input) {
0168     edm::ParameterSet puCuts;
0169     edm::ParameterSet nonPUCuts;
0170
0171     std::vector<std::string> inputNames = input.getParameterNames();
0172     for (auto const& cut : inputNames) {
0173       if (cut == "minTrackVertexWeight" || cut == "maxDeltaZ" || cut == "maxDeltaZToLeadTrack") {
0174         puCuts.copyFrom(input, cut);
0175       } else {
0176         nonPUCuts.copyFrom(input, cut);
0177       }
0178     }
0179     return std::make_pair(puCuts, nonPUCuts);
0180   }
0181
0182   bool RecoTauQualityCuts::filterTrack(const reco::TrackBaseRef& track) const {
0183     if (!filterTrack_(track.get()))
0184       return false;
0185     if (minTrackVertexWeight_ >= 0. && !(pv_->trackWeight(convertRef(track)) >= minTrackVertexWeight_))
0186       return false;
0187     return true;
0188   }
0189
0190   bool RecoTauQualityCuts::filterTrack(const reco::TrackRef& track) const {
0191     if (!filterTrack_(track.get()))
0192       return false;
0193     if (minTrackVertexWeight_ >= 0. && !(pv_->trackWeight(convertRef(track)) >= minTrackVertexWeight_))
0194       return false;
0195     return true;
0196   }
0197
0198   bool RecoTauQualityCuts::filterTrack(const reco::Track& track) const { return filterTrack_(&track); }
0199
0200   bool RecoTauQualityCuts::filterTrack_(const reco::Track* track) const {
0201     if (minTrackPt_ >= 0 && !(track->pt() > minTrackPt_))
0202       return false;
0203     if (maxTrackChi2_ >= 0 && !(track->normalizedChi2() <= maxTrackChi2_))
0204       return false;
0205     if (checkHitPattern_) {
0206       const reco::HitPattern& hitPattern = track->hitPattern();
0207       if (minTrackPixelHits_ > 0 && !(hitPattern.numberOfValidPixelHits() >= minTrackPixelHits_))
0208         return false;
0209       if (minTrackHits_ > 0 && !(hitPattern.numberOfValidHits() >= minTrackHits_))
0210         return false;
0211     }
0212     if (checkPV_ && pv_.isNull()) {
0213       edm::LogError("QCutsNoPrimaryVertex") << "Primary vertex Ref in "
0214                                             << "RecoTauQualityCuts is invalid. - filterTrack";
0215       return false;
0216     }
0217
0218     if (maxTransverseImpactParameter_ >= 0 &&
0219         !(std::fabs(track->dxy(pv_->position())) <= maxTransverseImpactParameter_))
0220       return false;
0221     if (maxDeltaZ_ >= 0 && !(std::fabs(track->dz(pv_->position())) <= maxDeltaZ_))
0222       return false;
0223     if (maxDeltaZToLeadTrack_ >= 0) {
0224       if (!leadTrack_) {
0225         edm::LogError("QCutsNoValidLeadTrack") << "Lead track Ref in "
0226                                                << "RecoTauQualityCuts is invalid. - filterTrack";
0227         return false;
0228       }
0229
0230       if (!(std::fabs(track->dz(pv_->position()) - leadTrack_->dz(pv_->position())) <= maxDeltaZToLeadTrack_))
0231         return false;
0232     }
0233
0234     return true;
0235   }
0236
0237   bool RecoTauQualityCuts::filterChargedCand(const reco::Candidate& cand) const {
0238     if (cand.charge() == 0)
0239       return true;
0240     const pat::PackedCandidate* pCand = dynamic_cast<const pat::PackedCandidate*>(&cand);
0241     if (pCand == nullptr)
0242       return true;
0243
0244     //Get track, it should be present for cands with pT(charged)>0.5GeV
0245     //and check track quality critera other than vertex weight
0246     auto track = getTrack(cand);
0247     if (track != nullptr) {
0248       if (!filterTrack(*track))
0249         return false;
0250     } else {  //Candidates without track (pT(charged)<0.5GeV): Can still check pT and calculate dxy and dz
0251       if (minTrackPt_ >= 0 && !(pCand->pt() > minTrackPt_))
0252         return false;
0253       if (checkPV_ && pv_.isNull()) {
0254         edm::LogError("QCutsNoPrimaryVertex") << "Primary vertex Ref in "
0255                                               << "RecoTauQualityCuts is invalid. - filterChargedCand";
0256         return false;
0257       }
0258
0259       if (maxTransverseImpactParameter_ >= 0 &&
0260           !(std::fabs(pCand->dxy(pv_->position())) <= maxTransverseImpactParameter_))
0261         return false;
0262       if (maxDeltaZ_ >= 0 && !(std::fabs(pCand->dz(pv_->position())) <= maxDeltaZ_))
0263         return false;
0264       if (maxDeltaZToLeadTrack_ >= 0) {
0265         if (leadTrack_ == nullptr) {
0266           edm::LogError("QCutsNoValidLeadTrack") << "Lead track Ref in "
0267                                                  << "RecoTauQualityCuts is invalid. - filterChargedCand";
0268           return false;
0269         }
0270
0271         if (!(std::fabs(pCand->dz(pv_->position()) - leadTrack_->dz(pv_->position())) <= maxDeltaZToLeadTrack_))
0272           return false;
0273       }
0274     }
0275     if (minTrackVertexWeight_ >= 0. && !(qcuts::minPackedCandVertexWeight(*pCand, &pv_, minTrackVertexWeight_)))
0276       return false;
0277
0278     return true;
0279   }
0280
0281   bool RecoTauQualityCuts::filterGammaCand(const reco::Candidate& cand) const {
0282     if (minGammaEt_ >= 0 && !(cand.et() > minGammaEt_))
0283       return false;
0284     return true;
0285   }
0286
0287   bool RecoTauQualityCuts::filterNeutralHadronCand(const reco::Candidate& cand) const {
0288     if (minNeutralHadronEt_ >= 0 && !(cand.et() > minNeutralHadronEt_))
0289       return false;
0290     return true;
0291   }
0292
0293   bool RecoTauQualityCuts::filterCandByType(const reco::Candidate& cand) const {
0294     switch (std::abs(cand.pdgId())) {
0295       case 22:
0296         return filterGammaCand(cand);
0297       case 130:
0298         return filterNeutralHadronCand(cand);
0299       // We use the same qcuts for muons/electrons and charged hadrons.
0300       case 211:
0301       case 11:
0302       case 13:
0303         // no cuts ATM (track cuts applied in filterCand)
0304         return true;
0305       // Return false if we dont' know how to deal with this particle type
0306       default:
0307         return false;
0308     };
0309     return false;
0310   }
0311
0312   bool RecoTauQualityCuts::filterCand(const reco::Candidate& cand) const {
0313     auto trackRef = getTrackRef(cand);
0314     bool result = true;
0315
0316     if (trackRef.isNonnull()) {
0317       result = filterTrack(trackRef);
0318     } else {
0319       auto gsfTrackRef = getGsfTrackRef(cand);
0320       if (gsfTrackRef.isNonnull())
0321         result = filterTrack(gsfTrackRef);
0322       else if (cand.charge() != 0) {
0323         result = filterChargedCand(cand);
0324       }
0325     }
0326
0327     if (result)
0328       result = filterCandByType(cand);
0329
0330     return result;
0331   }
0332
0333   void RecoTauQualityCuts::setLeadTrack(const reco::Track& leadTrack) { leadTrack_ = &leadTrack; }
0334
0335   void RecoTauQualityCuts::setLeadTrack(const reco::Candidate& leadCand) { leadTrack_ = getTrack(leadCand); }
0336
0337   void RecoTauQualityCuts::setLeadTrack(const reco::CandidateRef& leadCand) {
0338     if (leadCand.isNonnull()) {
0339       leadTrack_ = getTrack(*leadCand);
0340     } else {
0341       // Set null
0342       leadTrack_ = nullptr;
0343     }
0344   }
0345
0346   void RecoTauQualityCuts::fillDescriptions(edm::ParameterSetDescription& desc_qualityCuts) {
0347     edm::ParameterSetDescription desc_signalQualityCuts;
0348     desc_signalQualityCuts.add<double>("minTrackPt", 0.5);
0349     desc_signalQualityCuts.add<double>("maxTrackChi2", 100.0);
0350     desc_signalQualityCuts.add<double>("maxTransverseImpactParameter", 0.1);
0351     desc_signalQualityCuts.add<double>("maxDeltaZ", 0.4);
0352     desc_signalQualityCuts.add<double>("maxDeltaZToLeadTrack", -1.0);  // by default disabled
0353     desc_signalQualityCuts.add<double>("minTrackVertexWeight", -1.0);  // by default disabled
0354     desc_signalQualityCuts.add<unsigned int>("minTrackPixelHits", 0);
0355     desc_signalQualityCuts.add<unsigned int>("minTrackHits", 3);
0356     desc_signalQualityCuts.add<double>("minGammaEt", 1.0);
0357     desc_signalQualityCuts.addOptional<bool>("useTracksInsteadOfPFHadrons");
0358     desc_signalQualityCuts.add<double>("minNeutralHadronEt", 30.0);
0359
0360     edm::ParameterSetDescription desc_isolationQualityCuts;
0361     desc_isolationQualityCuts.add<double>("minTrackPt", 1.0);
0362     desc_isolationQualityCuts.add<double>("maxTrackChi2", 100.0);
0363     desc_isolationQualityCuts.add<double>("maxTransverseImpactParameter", 0.03);
0364     desc_isolationQualityCuts.add<double>("maxDeltaZ", 0.2);
0365     desc_isolationQualityCuts.add<double>("maxDeltaZToLeadTrack", -1.0);  // by default disabled
0366     desc_isolationQualityCuts.add<double>("minTrackVertexWeight", -1.0);  // by default disabled
0367     desc_isolationQualityCuts.add<unsigned int>("minTrackPixelHits", 0);
0368     desc_isolationQualityCuts.add<unsigned int>("minTrackHits", 8);
0369     desc_isolationQualityCuts.add<double>("minGammaEt", 1.5);
0370     desc_isolationQualityCuts.addOptional<bool>("useTracksInsteadOfPFHadrons");
0371
0372     edm::ParameterSetDescription desc_vxAssocQualityCuts;
0373     desc_vxAssocQualityCuts.add<double>("minTrackPt", 0.5);
0374     desc_vxAssocQualityCuts.add<double>("maxTrackChi2", 100.0);
0375     desc_vxAssocQualityCuts.add<double>("maxTransverseImpactParameter", 0.1);
0376     desc_vxAssocQualityCuts.add<double>("minTrackVertexWeight", -1.0);
0377     desc_vxAssocQualityCuts.add<unsigned int>("minTrackPixelHits", 0);
0378     desc_vxAssocQualityCuts.add<unsigned int>("minTrackHits", 3);
0379     desc_vxAssocQualityCuts.add<double>("minGammaEt", 1.0);
0380     desc_vxAssocQualityCuts.addOptional<bool>("useTracksInsteadOfPFHadrons");
0381
0382     desc_qualityCuts.add<edm::ParameterSetDescription>("signalQualityCuts", desc_signalQualityCuts);
0383     desc_qualityCuts.add<edm::ParameterSetDescription>("isolationQualityCuts", desc_isolationQualityCuts);
0384     desc_qualityCuts.add<edm::ParameterSetDescription>("vxAssocQualityCuts", desc_vxAssocQualityCuts);
0385     desc_qualityCuts.add<edm::InputTag>("primaryVertexSrc", edm::InputTag("offlinePrimaryVertices"));
0386     desc_qualityCuts.add<std::string>("pvFindingAlgo", "closestInDeltaZ");
0387     desc_qualityCuts.add<bool>("vertexTrackFiltering", false);
0388     desc_qualityCuts.add<bool>("recoverLeadingTrk", false);
0389     desc_qualityCuts.add<std::string>("leadingTrkOrPFCandOption", "leadPFCand");
0390   }
0391
0392 }  // end namespace reco::tau