Back to home page

Project CMSSW displayed by LXR
 
 

    

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   // Apply a hypothesis on the mass of the strips.
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 }  // namespace
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() {}  // CV: maxMass object gets deleted by PFRecoTauDiscriminationByHPSSelection destructor
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   // Get the matchign cut
0066   matchingCone_ = pset.getParameter<double>("matchingCone");
0067   minPt_ = pset.getParameter<double>("minTauPt");
0068   // Get the mass cuts for each decay mode
0069   typedef std::vector<edm::ParameterSet> VPSet;
0070   const VPSet& decayModes = pset.getParameter<VPSet>("decayModes");
0071   for (auto const& decayMode : decayModes) {
0072     // The mass window(s)
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         // The decay mode as a key
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 }  // namespace
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   // Check if we pass the min pt
0128   if (tau->pt() < minPt_) {
0129     if (verbosity_) {
0130       edm::LogPrint("PFTauByHPSSelect") << " fails minPt cut.";
0131     }
0132     return 0.0;
0133   }
0134 
0135   // See if we select this decay mode
0136   DecayModeCutMap::const_iterator massWindowIter = decayModeCuts_.find(
0137       std::make_pair(tau->signalTauChargedHadronCandidates().size(), tau->signalPiZeroCandidates().size()));
0138   // Check if decay mode is supported
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   // Find the total pizero p4
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   // Apply strip mass assumption corrections
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   // Check if tau fails mass cut
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   // CV: require that mass of charged tau decay products is always within specified mass window,
0227   //     irrespective of bendCorrection_mass
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   // Check if it fails the pi0 IM cut
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   // Check if tau passes matching cone cut
0251   //edm::LogPrint("PFTauByHPSSelect") << "dR(tau, jet) = " << deltaR(tauP4, tau->jetRef()->p4()) ;
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   // Check if tau passes cone cut
0260   double cone_size = tau->signalConeSize();
0261   // Check if any charged objects fail the signal cone cut
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   // Now check the pizeros
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;  //no. of charged prongs expected for this tau
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     //MB: how to deal with tau@miniAOD case?
0326     if (nTracks < nProngs) {  //"lost track" probably used to build this tau
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   // Otherwise, we pass!
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   // hpsSelectionDiscriminator
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!