File indexing completed on 2023-03-17 11:21:59
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
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
0047 template <typename T>
0048 reco::TrackBaseRef convertRef(const T& ref) {
0049 return reco::TrackBaseRef(ref);
0050 }
0051 }
0052
0053
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
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;
0067 else if (quality == pat::PackedCandidate::UsedInFitLoose)
0068 weight = 0.1;
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 }
0078
0079 RecoTauQualityCuts::RecoTauQualityCuts(const edm::ParameterSet& qcuts) {
0080
0081 CandQCutFuncCollection chargedHadronCuts;
0082 CandQCutFuncCollection gammaCuts;
0083 CandQCutFuncCollection neutralHadronCuts;
0084
0085
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
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
0120 minTrackVertexWeight_ = getDouble("minTrackVertexWeight");
0121
0122
0123 checkHitPattern_ = (minTrackHits_ > 0) || (minTrackPixelHits_ > 0);
0124 checkPV_ = (maxTransverseImpactParameter_ >= 0) || (maxDeltaZ_ >= 0) || (maxDeltaZToLeadTrack_ >= 0) ||
0125 (minTrackVertexWeight_ >= 0);
0126
0127
0128 minGammaEt_ = getDouble("minGammaEt");
0129
0130
0131 minNeutralHadronEt_ = getDouble("minNeutralHadronEt");
0132
0133
0134 if (!passedOptionSet.empty()) {
0135 std::string unParsedOptions;
0136 bool thereIsABadParameter = false;
0137 for (auto const& option : passedOptionSet) {
0138
0139 if (option == "useTracksInsteadOfPFHadrons") {
0140
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
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
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
0245
0246 auto track = getTrack(cand);
0247 if (track != nullptr) {
0248 if (!filterTrack(*track))
0249 return false;
0250 } else {
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
0300 case 211:
0301 case 11:
0302 case 13:
0303
0304 return true;
0305
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
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);
0353 desc_signalQualityCuts.add<double>("minTrackVertexWeight", -1.0);
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);
0366 desc_isolationQualityCuts.add<double>("minTrackVertexWeight", -1.0);
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 }