File indexing completed on 2023-03-17 11:21:51
0001 #include "RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h"
0002 #include "FWCore/Utilities/interface/InputTag.h"
0003
0004 #include <FWCore/ParameterSet/interface/ConfigurationDescriptions.h>
0005 #include <FWCore/ParameterSet/interface/ParameterSetDescription.h>
0006
0007 #include "CommonTools/Utils/interface/StringObjectFunction.h"
0008 #include "DataFormats/Math/interface/deltaR.h"
0009
0010 #include "DataFormats/TauReco/interface/PFRecoTauChargedHadron.h"
0011 #include "DataFormats/TauReco/interface/PFRecoTauChargedHadronFwd.h"
0012 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0013
0014 namespace {
0015
0016 math::XYZTLorentzVector applyMassConstraint(const math::XYZTLorentzVector& vec, double mass) {
0017 double factor = sqrt(vec.energy() * vec.energy() - mass * mass) / vec.P();
0018 return math::XYZTLorentzVector(vec.px() * factor, vec.py() * factor, vec.pz() * factor, vec.energy());
0019 }
0020 }
0021
0022 class PFRecoTauDiscriminationByHPSSelection : public PFTauDiscriminationProducerBase {
0023 public:
0024 explicit PFRecoTauDiscriminationByHPSSelection(const edm::ParameterSet&);
0025 ~PFRecoTauDiscriminationByHPSSelection() override;
0026 double discriminate(const reco::PFTauRef&) const override;
0027
0028 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0029
0030 private:
0031 typedef StringObjectFunction<reco::PFTau> TauFunc;
0032
0033 struct DecayModeCuts {
0034 DecayModeCuts() : maxMass_(nullptr) {}
0035 ~DecayModeCuts() {}
0036 unsigned nTracksMin_;
0037 unsigned nChargedPFCandsMin_;
0038 double minMass_;
0039 TauFunc* maxMass_;
0040 bool applyBendCorrection_mass_;
0041 bool applyBendCorrection_eta_;
0042 bool applyBendCorrection_phi_;
0043 double minPi0Mass_;
0044 double maxPi0Mass_;
0045 double assumeStripMass_;
0046 };
0047
0048 typedef std::pair<unsigned int, unsigned int> IntPair;
0049 typedef std::pair<double, double> DoublePair;
0050 typedef std::map<IntPair, DecayModeCuts> DecayModeCutMap;
0051
0052 DecayModeCutMap decayModeCuts_;
0053 double matchingCone_;
0054 double minPt_;
0055
0056 bool requireTauChargedHadronsToBeChargedPFCands_;
0057
0058 int minPixelHits_;
0059
0060 int verbosity_;
0061 };
0062
0063 PFRecoTauDiscriminationByHPSSelection::PFRecoTauDiscriminationByHPSSelection(const edm::ParameterSet& pset)
0064 : PFTauDiscriminationProducerBase(pset) {
0065
0066 matchingCone_ = pset.getParameter<double>("matchingCone");
0067 minPt_ = pset.getParameter<double>("minTauPt");
0068
0069 typedef std::vector<edm::ParameterSet> VPSet;
0070 const VPSet& decayModes = pset.getParameter<VPSet>("decayModes");
0071 for (auto const& decayMode : decayModes) {
0072
0073 DecayModeCuts cuts;
0074 cuts.nTracksMin_ = decayMode.getParameter<unsigned>("nTracksMin");
0075 cuts.nChargedPFCandsMin_ = decayMode.getParameter<unsigned>("nChargedPFCandsMin");
0076 cuts.minMass_ = decayMode.getParameter<double>("minMass");
0077 cuts.maxMass_ = new TauFunc(decayMode.getParameter<std::string>("maxMass"));
0078 edm::ParameterSet applyBendCorrection = decayMode.getParameter<edm::ParameterSet>("applyBendCorrection");
0079 cuts.applyBendCorrection_eta_ = applyBendCorrection.getParameter<bool>("eta");
0080 cuts.applyBendCorrection_phi_ = applyBendCorrection.getParameter<bool>("phi");
0081 cuts.applyBendCorrection_mass_ = applyBendCorrection.getParameter<bool>("mass");
0082 cuts.minPi0Mass_ = decayMode.getParameter<double>("minPi0Mass");
0083 cuts.maxPi0Mass_ = decayMode.getParameter<double>("maxPi0Mass");
0084 cuts.assumeStripMass_ = decayMode.getParameter<double>("assumeStripMass");
0085 decayModeCuts_.insert(std::make_pair(
0086
0087 std::make_pair(decayMode.getParameter<uint32_t>("nCharged"), decayMode.getParameter<uint32_t>("nPiZeros")),
0088 cuts));
0089 }
0090 requireTauChargedHadronsToBeChargedPFCands_ = pset.getParameter<bool>("requireTauChargedHadronsToBeChargedPFCands");
0091 minPixelHits_ = pset.getParameter<int>("minPixelHits");
0092 verbosity_ = pset.getParameter<int>("verbosity");
0093 }
0094
0095 PFRecoTauDiscriminationByHPSSelection::~PFRecoTauDiscriminationByHPSSelection() {
0096 for (DecayModeCutMap::iterator it = decayModeCuts_.begin(); it != decayModeCuts_.end(); ++it) {
0097 delete it->second.maxMass_;
0098 }
0099 }
0100
0101 namespace {
0102 inline const reco::Track* getTrack(const reco::Candidate& cand) {
0103 const reco::PFCandidate* pfCandPtr = dynamic_cast<const reco::PFCandidate*>(&cand);
0104 if (pfCandPtr) {
0105 if (pfCandPtr->trackRef().isNonnull())
0106 return pfCandPtr->trackRef().get();
0107 else if (pfCandPtr->gsfTrackRef().isNonnull())
0108 return pfCandPtr->gsfTrackRef().get();
0109 else
0110 return nullptr;
0111 }
0112 const pat::PackedCandidate* packedCand = dynamic_cast<const pat::PackedCandidate*>(&cand);
0113 if (packedCand && packedCand->hasTrackDetails())
0114 return &packedCand->pseudoTrack();
0115
0116 return nullptr;
0117 }
0118 }
0119
0120 double PFRecoTauDiscriminationByHPSSelection::discriminate(const reco::PFTauRef& tau) const {
0121 if (verbosity_) {
0122 edm::LogPrint("PFTauByHPSSelect") << "<PFRecoTauDiscriminationByHPSSelection::discriminate>:";
0123 edm::LogPrint("PFTauByHPSSelect") << " nCharged = " << tau->signalTauChargedHadronCandidates().size();
0124 edm::LogPrint("PFTauByHPSSelect") << " nPiZeros = " << tau->signalPiZeroCandidates().size();
0125 }
0126
0127
0128 if (tau->pt() < minPt_) {
0129 if (verbosity_) {
0130 edm::LogPrint("PFTauByHPSSelect") << " fails minPt cut.";
0131 }
0132 return 0.0;
0133 }
0134
0135
0136 DecayModeCutMap::const_iterator massWindowIter = decayModeCuts_.find(
0137 std::make_pair(tau->signalTauChargedHadronCandidates().size(), tau->signalPiZeroCandidates().size()));
0138
0139 if (massWindowIter == decayModeCuts_.end()) {
0140 if (verbosity_) {
0141 edm::LogPrint("PFTauByHPSSelect") << " fails mass-window definition requirement.";
0142 }
0143 return 0.0;
0144 }
0145 const DecayModeCuts& massWindow = massWindowIter->second;
0146
0147 if (massWindow.nTracksMin_ > 0) {
0148 unsigned int nTracks = 0;
0149 const std::vector<reco::PFRecoTauChargedHadron>& chargedHadrons = tau->signalTauChargedHadronCandidates();
0150 for (std::vector<reco::PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();
0151 chargedHadron != chargedHadrons.end();
0152 ++chargedHadron) {
0153 if (chargedHadron->algoIs(reco::PFRecoTauChargedHadron::kChargedPFCandidate) ||
0154 chargedHadron->algoIs(reco::PFRecoTauChargedHadron::kTrack)) {
0155 ++nTracks;
0156 }
0157 }
0158 if (verbosity_) {
0159 edm::LogPrint("PFTauByHPSSelect") << "nTracks = " << nTracks << " (min = " << massWindow.nTracksMin_ << ")";
0160 }
0161 if (nTracks < massWindow.nTracksMin_) {
0162 if (verbosity_) {
0163 edm::LogPrint("PFTauByHPSSelect") << " fails nTracks requirement for mass-window.";
0164 }
0165 return 0.0;
0166 }
0167 }
0168 if (massWindow.nChargedPFCandsMin_ > 0) {
0169 unsigned int nChargedPFCands = 0;
0170 const std::vector<reco::PFRecoTauChargedHadron>& chargedHadrons = tau->signalTauChargedHadronCandidates();
0171 for (std::vector<reco::PFRecoTauChargedHadron>::const_iterator chargedHadron = chargedHadrons.begin();
0172 chargedHadron != chargedHadrons.end();
0173 ++chargedHadron) {
0174 if (chargedHadron->algoIs(reco::PFRecoTauChargedHadron::kChargedPFCandidate)) {
0175 ++nChargedPFCands;
0176 }
0177 }
0178 if (verbosity_) {
0179 edm::LogPrint("PFTauByHPSSelect") << "nChargedPFCands = " << nChargedPFCands
0180 << " (min = " << massWindow.nChargedPFCandsMin_ << ")";
0181 }
0182 if (nChargedPFCands < massWindow.nChargedPFCandsMin_) {
0183 if (verbosity_) {
0184 edm::LogPrint("PFTauByHPSSelect") << " fails nChargedPFCands requirement for mass-window.";
0185 }
0186 return 0.0;
0187 }
0188 }
0189
0190 math::XYZTLorentzVector tauP4 = tau->p4();
0191 if (verbosity_) {
0192 edm::LogPrint("PFTauByHPSSelect") << "tau: Pt = " << tauP4.pt() << ", eta = " << tauP4.eta()
0193 << ", phi = " << tauP4.phi() << ", mass = " << tauP4.mass();
0194 }
0195
0196 reco::Candidate::LorentzVector stripsP4;
0197 for (auto const& cand : tau->signalPiZeroCandidates()) {
0198 const math::XYZTLorentzVector& candP4 = cand.p4();
0199 stripsP4 += candP4;
0200 }
0201
0202
0203 if (massWindow.assumeStripMass_ >= 0) {
0204 for (auto const& cand : tau->signalPiZeroCandidates()) {
0205 const math::XYZTLorentzVector& uncorrected = cand.p4();
0206 math::XYZTLorentzVector corrected = applyMassConstraint(uncorrected, massWindow.assumeStripMass_);
0207 math::XYZTLorentzVector correction = corrected - uncorrected;
0208 tauP4 += correction;
0209 stripsP4 += correction;
0210 }
0211 }
0212 if (verbosity_) {
0213 edm::LogPrint("PFTauByHPSSelect") << "strips: Pt = " << stripsP4.pt() << ", eta = " << stripsP4.eta()
0214 << ", phi = " << stripsP4.phi() << ", mass = " << stripsP4.mass();
0215 }
0216
0217
0218 double maxMass_value = (*massWindow.maxMass_)(*tau);
0219 double bendCorrection_mass = (massWindow.applyBendCorrection_mass_) ? tau->bendCorrMass() : 0.;
0220 if (!((tauP4.M() - bendCorrection_mass) < maxMass_value && (tauP4.M() + bendCorrection_mass) > massWindow.minMass_)) {
0221 if (verbosity_) {
0222 edm::LogPrint("PFTauByHPSSelect") << " fails tau mass-window cut.";
0223 }
0224 return 0.0;
0225 }
0226
0227
0228 reco::Candidate::LorentzVector tauP4_charged;
0229 const std::vector<reco::PFRecoTauChargedHadron>& signalChargedHadrons = tau->signalTauChargedHadronCandidates();
0230 for (std::vector<reco::PFRecoTauChargedHadron>::const_iterator signalChargedHadron = signalChargedHadrons.begin();
0231 signalChargedHadron != signalChargedHadrons.end();
0232 ++signalChargedHadron) {
0233 tauP4_charged += signalChargedHadron->p4();
0234 }
0235 if (!(tauP4_charged.mass() < maxMass_value)) {
0236 if (verbosity_) {
0237 edm::LogPrint("PFTauByHPSSelect") << " fails tau mass-window cut.";
0238 }
0239 return 0.0;
0240 }
0241
0242
0243 if (stripsP4.M() > massWindow.maxPi0Mass_ || stripsP4.M() < massWindow.minPi0Mass_) {
0244 if (verbosity_) {
0245 edm::LogPrint("PFTauByHPSSelect") << " fails strip mass-window cut.";
0246 }
0247 return 0.0;
0248 }
0249
0250
0251
0252 if (deltaR(tauP4, tau->jetRef()->p4()) > matchingCone_) {
0253 if (verbosity_) {
0254 edm::LogPrint("PFTauByHPSSelect") << " fails matching-cone cut.";
0255 }
0256 return 0.0;
0257 }
0258
0259
0260 double cone_size = tau->signalConeSize();
0261
0262 for (auto const& cand : tau->signalTauChargedHadronCandidates()) {
0263 if (verbosity_) {
0264 edm::LogPrint("PFTauByHPSSelect") << "dR(tau, signalPFChargedHadr) = " << deltaR(cand.p4(), tauP4);
0265 }
0266 if (deltaR(cand.p4(), tauP4) > cone_size) {
0267 if (verbosity_) {
0268 edm::LogPrint("PFTauByHPSSelect") << " fails signal-cone cut for charged hadron(s).";
0269 }
0270 return 0.0;
0271 }
0272 }
0273
0274 for (auto const& cand : tau->signalPiZeroCandidates()) {
0275 double bendCorrection_eta = (massWindow.applyBendCorrection_eta_) ? cand.bendCorrEta() : 0.;
0276 double dEta = std::max(0., fabs(cand.eta() - tauP4.eta()) - bendCorrection_eta);
0277 double bendCorrection_phi = (massWindow.applyBendCorrection_phi_) ? cand.bendCorrPhi() : 0.;
0278 double dPhi = std::max(0., std::abs(reco::deltaPhi(cand.phi(), tauP4.phi())) - bendCorrection_phi);
0279 double dR2 = dEta * dEta + dPhi * dPhi;
0280 if (verbosity_) {
0281 edm::LogPrint("PFTauByHPSSelect") << "dR2(tau, signalPiZero) = " << dR2;
0282 }
0283 if (dR2 > cone_size * cone_size) {
0284 if (verbosity_) {
0285 edm::LogPrint("PFTauByHPSSelect") << " fails signal-cone cut for strip(s).";
0286 }
0287 return 0.0;
0288 }
0289 }
0290
0291 if (requireTauChargedHadronsToBeChargedPFCands_) {
0292 for (auto const& cand : tau->signalTauChargedHadronCandidates()) {
0293 if (verbosity_) {
0294 std::string algo_string;
0295 if (cand.algo() == reco::PFRecoTauChargedHadron::kChargedPFCandidate)
0296 algo_string = "ChargedPFCandidate";
0297 else if (cand.algo() == reco::PFRecoTauChargedHadron::kTrack)
0298 algo_string = "Track";
0299 else if (cand.algo() == reco::PFRecoTauChargedHadron::kPFNeutralHadron)
0300 algo_string = "PFNeutralHadron";
0301 else
0302 algo_string = "Undefined";
0303 edm::LogPrint("PFTauByHPSSelect") << "algo(signalPFChargedHadr) = " << algo_string;
0304 }
0305 if (!(cand.algo() == reco::PFRecoTauChargedHadron::kChargedPFCandidate)) {
0306 if (verbosity_) {
0307 edm::LogPrint("PFTauByHPSSelect") << " fails cut on PFRecoTauChargedHadron algo.";
0308 }
0309 return 0.0;
0310 }
0311 }
0312 }
0313
0314 if (minPixelHits_ > 0) {
0315 int numPixelHits = 0;
0316 const int nProngs = tau->decayMode() / 5 + 1;
0317 int nTracks = 0;
0318 for (const auto& chargedHadrCand : tau->signalChargedHadrCands()) {
0319 const reco::Track* track = getTrack(*chargedHadrCand);
0320 if (track != nullptr) {
0321 numPixelHits += track->hitPattern().numberOfValidPixelHits();
0322 nTracks++;
0323 }
0324 }
0325
0326 if (nTracks < nProngs) {
0327 for (const auto& track : tau->signalTracks()) {
0328 if (track.isNonnull()) {
0329 numPixelHits += track->hitPattern().numberOfValidPixelHits();
0330 nTracks++;
0331 }
0332 }
0333 }
0334 if (!(numPixelHits >= minPixelHits_)) {
0335 if (verbosity_) {
0336 edm::LogPrint("PFTauByHPSSelect") << " fails cut on sum of pixel hits.";
0337 }
0338 return 0.0;
0339 }
0340 }
0341
0342
0343 if (verbosity_) {
0344 edm::LogPrint("PFTauByHPSSelect") << " passes all cuts.";
0345 }
0346 return 1.0;
0347 }
0348
0349 void PFRecoTauDiscriminationByHPSSelection::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0350
0351 edm::ParameterSetDescription desc;
0352 desc.add<edm::InputTag>("PFTauProducer", edm::InputTag("combinatoricRecoTaus"));
0353 desc.add<int>("verbosity", 0);
0354 desc.add<double>("minTauPt", 0.0);
0355 {
0356 edm::ParameterSetDescription psd0;
0357 psd0.add<std::string>("BooleanOperator", "and");
0358 desc.add<edm::ParameterSetDescription>("Prediscriminants", psd0);
0359 }
0360
0361 {
0362 edm::ParameterSetDescription vpset_decayModes;
0363 vpset_decayModes.add<double>("minPi0Mass", -1.e3);
0364 vpset_decayModes.add<std::string>("maxMass");
0365 vpset_decayModes.add<double>("maxPi0Mass", 1.e9);
0366 vpset_decayModes.add<unsigned int>("nPiZeros");
0367 vpset_decayModes.add<double>("minMass");
0368 vpset_decayModes.add<unsigned int>("nChargedPFCandsMin", 0);
0369 vpset_decayModes.add<unsigned int>("nTracksMin", 0);
0370 vpset_decayModes.add<unsigned int>("nCharged");
0371 {
0372 edm::ParameterSetDescription psd0;
0373 psd0.add<bool>("phi");
0374 psd0.add<bool>("eta");
0375 psd0.add<bool>("mass");
0376 vpset_decayModes.add<edm::ParameterSetDescription>("applyBendCorrection", psd0);
0377 }
0378 vpset_decayModes.add<double>("assumeStripMass", -1.0);
0379 std::vector<edm::ParameterSet> vpset_default;
0380 {
0381 edm::ParameterSet pset;
0382 pset.addParameter<double>("minPi0Mass", -1.e3);
0383 pset.addParameter<double>("maxPi0Mass", 1.e9);
0384 pset.addParameter<unsigned int>("nChargedPFCandsMin", 0);
0385 pset.addParameter<unsigned int>("nTracksMin", 0);
0386 pset.addParameter<double>("assumeStripMass", -1.0);
0387 vpset_default.push_back(pset);
0388 }
0389 desc.addVPSet("decayModes", vpset_decayModes, vpset_default);
0390 }
0391
0392 desc.add<double>("matchingCone", 0.5);
0393 desc.add<int>("minPixelHits", 1);
0394 desc.add<bool>("requireTauChargedHadronsToBeChargedPFCands", false);
0395 descriptions.add("hpsSelectionDiscriminator", desc);
0396 }
0397
0398 DEFINE_FWK_MODULE(PFRecoTauDiscriminationByHPSSelection);