Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:48 

0001 /** \class PFRecoTauDiscriminationAgainstMuonSimple
0002  *
0003  * Compute tau Id. discriminator against muons for MiniAOD.
0004  * 
0005  * \author Michal Bluj, NCBJ Warsaw
0006  * based on PFRecoTauDiscriminationAgainstMuon2 by Christian Veelken
0007  *
0008  * Note: it is not granted that information on muon track is/will be always 
0009  * accesible with MiniAOD, if not one can consider to veto muons which are 
0010  * not only Trk by also STA or RPC muons, i.e. have many (>=2) good muon segments 
0011  */
0012 
0013 #include "RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h"
0014 
0015 #include "FWCore/Utilities/interface/Exception.h"
0016 
0017 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0018 #include "DataFormats/PatCandidates/interface/Muon.h"
0019 #include "DataFormats/Common/interface/Handle.h"
0020 #include "DataFormats/TrackReco/interface/HitPattern.h"
0021 #include "DataFormats/TrackReco/interface/Track.h"
0022 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0023 #include "DataFormats/Math/interface/deltaR.h"
0024 
0025 #include "RecoTauTag/RecoTau/interface/RecoTauMuonTools.h"
0026 
0027 #include <vector>
0028 #include <string>
0029 #include <iostream>
0030 #include <atomic>
0031 
0032 using reco::tau::format_vint;
0033 
0034 namespace {
0035 
0036   struct PFRecoTauDiscriminationAgainstMuonSimpleConfigSet {
0037     PFRecoTauDiscriminationAgainstMuonSimpleConfigSet(double hop, int mNOM, bool doCMV, int mNHL2S, int mNSTA, int mNRPC)
0038         : hop(hop),
0039           maxNumberOfMatches(mNOM),
0040           doCaloMuonVeto(doCMV),
0041           maxNumberOfHitsLast2Stations(mNHL2S),
0042           maxNumberOfSTAMuons(mNSTA),
0043           maxNumberOfRPCMuons(mNRPC) {}
0044 
0045     double hop;
0046     int maxNumberOfMatches;
0047     bool doCaloMuonVeto;
0048     int maxNumberOfHitsLast2Stations;
0049     int maxNumberOfSTAMuons, maxNumberOfRPCMuons;
0050   };
0051 
0052   class PFRecoTauDiscriminationAgainstMuonSimple final : public PFTauDiscriminationContainerProducerBase {
0053   public:
0054     explicit PFRecoTauDiscriminationAgainstMuonSimple(const edm::ParameterSet& cfg)
0055         : PFTauDiscriminationContainerProducerBase(cfg), moduleLabel_(cfg.getParameter<std::string>("@module_label")) {
0056       auto const wpDefs = cfg.getParameter<std::vector<edm::ParameterSet>>("IDWPdefinitions");
0057       for (auto& wpDef : wpDefs) {
0058         wpDefs_.push_back(
0059             PFRecoTauDiscriminationAgainstMuonSimpleConfigSet(wpDef.getParameter<double>("HoPMin"),
0060                                                               wpDef.getParameter<int>("maxNumberOfMatches"),
0061                                                               wpDef.getParameter<bool>("doCaloMuonVeto"),
0062                                                               wpDef.getParameter<int>("maxNumberOfHitsLast2Stations"),
0063                                                               wpDef.getParameter<int>("maxNumberOfSTAMuons"),
0064                                                               wpDef.getParameter<int>("maxNumberOfRPCMuons")));
0065       }
0066       srcPatMuons_ = cfg.getParameter<edm::InputTag>("srcPatMuons");
0067       muons_token = consumes<pat::MuonCollection>(srcPatMuons_);
0068       dRmuonMatch_ = cfg.getParameter<double>("dRmuonMatch");
0069       dRmuonMatchLimitedToJetArea_ = cfg.getParameter<bool>("dRmuonMatchLimitedToJetArea");
0070       minPtMatchedMuon_ = cfg.getParameter<double>("minPtMatchedMuon");
0071       typedef std::vector<int> vint;
0072       maskMatchesDT_ = cfg.getParameter<vint>("maskMatchesDT");
0073       maskMatchesCSC_ = cfg.getParameter<vint>("maskMatchesCSC");
0074       maskMatchesRPC_ = cfg.getParameter<vint>("maskMatchesRPC");
0075       maskHitsDT_ = cfg.getParameter<vint>("maskHitsDT");
0076       maskHitsCSC_ = cfg.getParameter<vint>("maskHitsCSC");
0077       maskHitsRPC_ = cfg.getParameter<vint>("maskHitsRPC");
0078       verbosity_ = cfg.getParameter<int>("verbosity");
0079     }
0080     ~PFRecoTauDiscriminationAgainstMuonSimple() override {}
0081 
0082     void beginEvent(const edm::Event&, const edm::EventSetup&) override;
0083 
0084     reco::SingleTauDiscriminatorContainer discriminate(const reco::PFTauRef&) const override;
0085 
0086     static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0087 
0088   private:
0089     std::string moduleLabel_;
0090     std::vector<PFRecoTauDiscriminationAgainstMuonSimpleConfigSet> wpDefs_;
0091     edm::InputTag srcPatMuons_;
0092     edm::Handle<pat::MuonCollection> muons_;
0093     edm::EDGetTokenT<pat::MuonCollection> muons_token;
0094     double dRmuonMatch_;
0095     bool dRmuonMatchLimitedToJetArea_;
0096     double minPtMatchedMuon_;
0097     std::vector<int> maskMatchesDT_;
0098     std::vector<int> maskMatchesCSC_;
0099     std::vector<int> maskMatchesRPC_;
0100     std::vector<int> maskHitsDT_;
0101     std::vector<int> maskHitsCSC_;
0102     std::vector<int> maskHitsRPC_;
0103 
0104     static std::atomic<unsigned int> numWarnings_;
0105     static constexpr unsigned int maxWarnings_ = 3;
0106     int verbosity_;
0107   };
0108 
0109   std::atomic<unsigned int> PFRecoTauDiscriminationAgainstMuonSimple::numWarnings_{0};
0110 
0111   void PFRecoTauDiscriminationAgainstMuonSimple::beginEvent(const edm::Event& evt, const edm::EventSetup& es) {
0112     evt.getByToken(muons_token, muons_);
0113   }
0114 
0115   reco::SingleTauDiscriminatorContainer PFRecoTauDiscriminationAgainstMuonSimple::discriminate(
0116       const reco::PFTauRef& pfTau) const {
0117     if (verbosity_) {
0118       edm::LogPrint("PFTauAgainstMuonSimple") << "<PFRecoTauDiscriminationAgainstMuonSimple::discriminate>:";
0119       edm::LogPrint("PFTauAgainstMuonSimple") << " moduleLabel = " << moduleLabel_;
0120       edm::LogPrint("PFTauAgainstMuonSimple")
0121           << "tau #" << pfTau.key() << ": Pt = " << pfTau->pt() << ", eta = " << pfTau->eta()
0122           << ", phi = " << pfTau->phi() << ", decay mode = " << pfTau->decayMode();
0123     }
0124 
0125     //energy fraction carried by leading charged hadron
0126     const reco::CandidatePtr& pfLeadChargedHadron = pfTau->leadChargedHadrCand();
0127     double caloEnergyFraction = 99;
0128     if (pfLeadChargedHadron.isNonnull()) {
0129       const pat::PackedCandidate* pCand = dynamic_cast<const pat::PackedCandidate*>(pfLeadChargedHadron.get());
0130       if (pCand != nullptr) {
0131         caloEnergyFraction = pCand->caloFraction();
0132         if (pCand->hasTrackDetails())  //MB: relate energy fraction to track momentum rather than to energy of candidate
0133           caloEnergyFraction *= pCand->energy() / pCand->bestTrack()->p();
0134         if (verbosity_) {
0135           edm::LogPrint("PFTauAgainstMuonSimple")
0136               << "decayMode = " << pfTau->decayMode() << ", caloEnergy(ECAL+HCAL)Fraction = " << caloEnergyFraction
0137               << ", leadPFChargedHadronP = " << pCand->p() << ", leadPFChargedHadron pdgId = " << pCand->pdgId();
0138         }
0139       }
0140     }
0141     //iterate over muons to find ones to use for discrimination
0142     std::vector<const pat::Muon*> muonsToCheck;
0143     size_t numMuons = muons_->size();
0144     for (size_t idxMuon = 0; idxMuon < numMuons; ++idxMuon) {
0145       bool matched = false;
0146       const pat::MuonRef muon(muons_, idxMuon);
0147       if (verbosity_)
0148         edm::LogPrint("PFTauAgainstMuonSimple") << "muon #" << muon.key() << ": Pt = " << muon->pt()
0149                                                 << ", eta = " << muon->eta() << ", phi = " << muon->phi();
0150       //firsty check if the muon corrsponds with the leading tau particle
0151       for (size_t iCand = 0; iCand < muon->numberOfSourceCandidatePtrs(); ++iCand) {
0152         const reco::CandidatePtr& srcCand = muon->sourceCandidatePtr(iCand);
0153         if (srcCand.isNonnull() && pfLeadChargedHadron.isNonnull() && srcCand.id() == pfLeadChargedHadron.id() &&
0154             srcCand.key() == pfLeadChargedHadron.key()) {
0155           muonsToCheck.push_back(muon.get());
0156           matched = true;
0157           if (verbosity_)
0158             edm::LogPrint("PFTauAgainstMuonSimple") << " tau has muonRef.";
0159           break;
0160         }
0161       }
0162       if (matched)
0163         continue;
0164       //check dR matching
0165       if (!(muon->pt() > minPtMatchedMuon_)) {
0166         if (verbosity_) {
0167           edm::LogPrint("PFTauAgainstMuonSimple") << " fails Pt cut --> skipping it.";
0168         }
0169         continue;
0170       }
0171       double dR = deltaR(muon->p4(), pfTau->p4());
0172       double dRmatch = dRmuonMatch_;
0173       if (dRmuonMatchLimitedToJetArea_) {
0174         double jetArea = 0.;
0175         if (pfTau->jetRef().isNonnull())
0176           jetArea = pfTau->jetRef()->jetArea();
0177         if (jetArea > 0.) {
0178           dRmatch = std::min(dRmatch, std::sqrt(jetArea / M_PI));
0179         } else {
0180           if (numWarnings_ < maxWarnings_) {
0181             edm::LogInfo("PFRecoTauDiscriminationAgainstMuonSimple::discriminate")
0182                 << "Jet associated to Tau: Pt = " << pfTau->pt() << ", eta = " << pfTau->eta()
0183                 << ", phi = " << pfTau->phi() << " has area = " << jetArea << " !!";
0184             ++numWarnings_;
0185           }
0186           dRmatch = pfTau->signalConeSize();
0187         }
0188         dRmatch = std::max(dRmatch, pfTau->signalConeSize());
0189         if (dR < dRmatch) {
0190           if (verbosity_)
0191             edm::LogPrint("PFTauAgainstMuonSimple") << " overlaps with tau, dR = " << dR;
0192           muonsToCheck.push_back(muon.get());
0193         }
0194       }
0195     }
0196     //now examine selected muons
0197     std::vector<int> numMatchesDT(4);
0198     std::vector<int> numMatchesCSC(4);
0199     std::vector<int> numMatchesRPC(4);
0200     std::vector<int> numHitsDT(4);
0201     std::vector<int> numHitsCSC(4);
0202     std::vector<int> numHitsRPC(4);
0203     //MB: clear counters of matched segments and hits globally as in the AgainstMuon2 discriminant, but note that they will sum for all matched muons. However it is not likely that there is more than one matched muon.
0204     for (int iStation = 0; iStation < 4; ++iStation) {
0205       numMatchesDT[iStation] = 0;
0206       numMatchesCSC[iStation] = 0;
0207       numMatchesRPC[iStation] = 0;
0208       numHitsDT[iStation] = 0;
0209       numHitsCSC[iStation] = 0;
0210       numHitsRPC[iStation] = 0;
0211     }
0212     int numSTAMuons = 0, numRPCMuons = 0;
0213     int numStationsWithMatches = 0;
0214     int numLast2StationsWithHits = 0;
0215     if (verbosity_ && !muonsToCheck.empty())
0216       edm::LogPrint("PFTauAgainstMuonSimple") << "Muons to check (" << muonsToCheck.size() << "):";
0217     size_t iMu = 0;
0218     for (const auto& mu : muonsToCheck) {
0219       if (mu->isStandAloneMuon())
0220         numSTAMuons++;
0221       if (mu->muonID("RPCMuLoose"))
0222         numRPCMuons++;
0223       reco::tau::countMatches(*mu, numMatchesDT, numMatchesCSC, numMatchesRPC);
0224       for (int iStation = 0; iStation < 4; ++iStation) {
0225         if (numMatchesDT[iStation] > 0 && !maskMatchesDT_[iStation])
0226           ++numStationsWithMatches;
0227         if (numMatchesCSC[iStation] > 0 && !maskMatchesCSC_[iStation])
0228           ++numStationsWithMatches;
0229         if (numMatchesRPC[iStation] > 0 && !maskMatchesRPC_[iStation])
0230           ++numStationsWithMatches;
0231       }
0232       reco::tau::countHits(*mu, numHitsDT, numHitsCSC, numHitsRPC);
0233       for (int iStation = 2; iStation < 4; ++iStation) {
0234         if (numHitsDT[iStation] > 0 && !maskHitsDT_[iStation])
0235           ++numLast2StationsWithHits;
0236         if (numHitsCSC[iStation] > 0 && !maskHitsCSC_[iStation])
0237           ++numLast2StationsWithHits;
0238         if (numHitsRPC[iStation] > 0 && !maskHitsRPC_[iStation])
0239           ++numLast2StationsWithHits;
0240       }
0241       if (verbosity_)
0242         edm::LogPrint("PFTauAgainstMuonSimple")
0243             << "\t" << iMu << ": Pt = " << mu->pt() << ", eta = " << mu->eta() << ", phi = " << mu->phi() << "\n\t"
0244             << "   isSTA: " << mu->isStandAloneMuon() << ", isRPCLoose: " << mu->muonID("RPCMuLoose")
0245             << "\n\t   numMatchesDT  = " << format_vint(numMatchesDT)
0246             << "\n\t   numMatchesCSC = " << format_vint(numMatchesCSC)
0247             << "\n\t   numMatchesRPC = " << format_vint(numMatchesRPC)
0248             << "\n\t   --> numStationsWithMatches = " << numStationsWithMatches
0249             << "\n\t   numHitsDT  = " << format_vint(numHitsDT) << "\n\t   numHitsCSC = " << format_vint(numHitsCSC)
0250             << "\n\t   numHitsRPC = " << format_vint(numHitsRPC)
0251             << "\n\t   --> numLast2StationsWithHits = " << numLast2StationsWithHits;
0252       iMu++;
0253     }
0254 
0255     reco::SingleTauDiscriminatorContainer result;
0256     for (auto const& wpDef : wpDefs_) {
0257       bool discriminatorValue = true;
0258       if (wpDef.maxNumberOfMatches >= 0 && numStationsWithMatches > wpDef.maxNumberOfMatches)
0259         discriminatorValue = false;
0260       if (wpDef.maxNumberOfHitsLast2Stations >= 0 && numLast2StationsWithHits > wpDef.maxNumberOfHitsLast2Stations)
0261         discriminatorValue = false;
0262       if (wpDef.maxNumberOfSTAMuons >= 0 && numSTAMuons > wpDef.maxNumberOfSTAMuons) {
0263         discriminatorValue = false;
0264       }
0265       if (wpDef.maxNumberOfRPCMuons >= 0 && numRPCMuons > wpDef.maxNumberOfRPCMuons) {
0266         discriminatorValue = false;
0267       }
0268       bool passesCaloMuonVeto = true;
0269       if (pfTau->decayMode() == 0 && caloEnergyFraction < wpDef.hop) {
0270         passesCaloMuonVeto = false;
0271       }
0272       if (wpDef.doCaloMuonVeto && !passesCaloMuonVeto) {
0273         discriminatorValue = false;
0274       }
0275       result.workingPoints.push_back(discriminatorValue);
0276       if (verbosity_)
0277         edm::LogPrint("PFTauAgainstMuonSimple") << "--> returning discriminatorValue = " << result.workingPoints.back();
0278     }
0279     return result;
0280   }
0281 
0282   void PFRecoTauDiscriminationAgainstMuonSimple::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0283     // pfRecoTauDiscriminationAgainstMuonSimple
0284     edm::ParameterSetDescription desc;
0285     desc.add<edm::InputTag>("PFTauProducer", edm::InputTag("pfRecoTauProducer"));
0286     {
0287       edm::ParameterSetDescription psd0;
0288       psd0.add<std::string>("BooleanOperator", "and");
0289       {
0290         edm::ParameterSetDescription psd1;
0291         psd1.add<double>("cut");              //, 0.5);
0292         psd1.add<edm::InputTag>("Producer");  //, edm::InputTag("pfRecoTauDiscriminationByLeadingTrackFinding"));
0293         psd0.addOptional<edm::ParameterSetDescription>("leadTrack", psd1);  // optional with default?
0294       }
0295       // Prediscriminants can be
0296       // Prediscriminants = noPrediscriminants,
0297       // as in RecoTauTag/Configuration/python/HPSPFTaus_cff.py
0298       //
0299       // and the definition is:
0300       // RecoTauTag/RecoTau/python/TauDiscriminatorTools.py
0301       // # Require no prediscriminants
0302       // noPrediscriminants = cms.PSet(
0303       //       BooleanOperator = cms.string("and"),
0304       //       )
0305       // -- so this is the minimum required definition
0306       // otherwise it inserts the leadTrack with Producer = InpuTag(...) and does not find the corresponding output during the run
0307       desc.add<edm::ParameterSetDescription>("Prediscriminants", psd0);
0308     }
0309     {
0310       edm::ParameterSetDescription desc_wp;
0311       desc_wp.add<std::string>("IDname");
0312       desc_wp.add<double>("HoPMin");
0313       desc_wp.add<int>("maxNumberOfMatches")->setComment("negative value would turn off this cut");
0314       desc_wp.add<bool>("doCaloMuonVeto");
0315       desc_wp.add<int>("maxNumberOfHitsLast2Stations")->setComment("negative value would turn off this cut");
0316       desc_wp.add<int>("maxNumberOfSTAMuons")->setComment("negative value would turn off this cut");
0317       desc_wp.add<int>("maxNumberOfRPCMuons")->setComment("negative value would turn off this cut");
0318       edm::ParameterSet pset_wp;
0319       pset_wp.addParameter<std::string>("IDname", "ByLooseMuonRejectionSimple");
0320       pset_wp.addParameter<double>("HoPMin", 0.2);
0321       pset_wp.addParameter<int>("maxNumberOfMatches", 1);
0322       pset_wp.addParameter<bool>("doCaloMuonVeto", true);
0323       pset_wp.addParameter<int>("maxNumberOfHitsLast2Stations", -1);
0324       pset_wp.addParameter<int>("maxNumberOfSTAMuons", -1);
0325       pset_wp.addParameter<int>("maxNumberOfRPCMuons", -1);
0326       std::vector<edm::ParameterSet> vpsd_wp;
0327       vpsd_wp.push_back(pset_wp);
0328       desc.addVPSet("IDWPdefinitions", desc_wp, vpsd_wp);
0329     }
0330     //add empty raw value config to simplify subsequent provenance searches
0331     desc.addVPSet("IDdefinitions", edm::ParameterSetDescription(), {});
0332 
0333     desc.add<edm::InputTag>("srcPatMuons", edm::InputTag("slimmedMuons"));
0334     desc.add<double>("dRmuonMatch", 0.3);
0335     desc.add<bool>("dRmuonMatchLimitedToJetArea", false);
0336     desc.add<double>("minPtMatchedMuon", 5.0);
0337     desc.add<std::vector<int>>("maskMatchesDT", {0, 0, 0, 0});
0338     desc.add<std::vector<int>>("maskMatchesCSC", {1, 0, 0, 0})
0339         ->setComment(
0340             "flags to mask/unmask DT, CSC and RPC chambers in individual muon stations. Segments and hits that are "
0341             "present in that muon station are ignored in case the 'mask' is set to 1. Per default only the innermost "
0342             "CSC "
0343             "chamber is ignored, as it is affected by spurious hits in high pile-up events.");
0344     desc.add<std::vector<int>>("maskMatchesRPC", {0, 0, 0, 0});
0345     desc.add<std::vector<int>>("maskHitsDT", {0, 0, 0, 0});
0346     desc.add<std::vector<int>>("maskHitsCSC", {0, 0, 0, 0});
0347     desc.add<std::vector<int>>("maskHitsRPC", {0, 0, 0, 0});
0348     desc.add<int>("verbosity", 0);
0349     descriptions.add("pfRecoTauDiscriminationAgainstMuonSimple", desc);
0350   }
0351 
0352 }  // namespace
0353 
0354 DEFINE_FWK_MODULE(PFRecoTauDiscriminationAgainstMuonSimple);
[ 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!