Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 /* class TauDiscriminationAgainstElectronMVA6
0002  * created : Nov 2 2015,
0003  * revised : May 29 2020,
0004  * Authors : Fabio Colombo (KIT)
0005  *           Anne-Catherine Le Bihan (IPHC),
0006  *           Michal Bluj (NCBJ)
0007  */
0008 
0009 #include "RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h"
0010 #include "RecoTauTag/RecoTau/interface/AntiElectronIDMVA6.h"
0011 #include "RecoTauTag/RecoTau/interface/PositionAtECalEntranceComputer.h"
0012 #include "DataFormats/Math/interface/deltaR.h"
0013 
0014 template <class TauType, class TauDiscriminator, class ElectronType>
0015 class TauDiscriminationAgainstElectronMVA6 : public TauDiscriminationProducerBase<TauType,
0016                                                                                   reco::TauDiscriminatorContainer,
0017                                                                                   reco::SingleTauDiscriminatorContainer,
0018                                                                                   TauDiscriminator> {
0019 public:
0020   typedef std::vector<TauType> TauCollection;
0021   typedef edm::Ref<TauCollection> TauRef;
0022   typedef std::vector<ElectronType> ElectronCollection;
0023 
0024   explicit TauDiscriminationAgainstElectronMVA6(const edm::ParameterSet& cfg)
0025       : TauDiscriminationProducerBase<TauType,
0026                                       reco::TauDiscriminatorContainer,
0027                                       reco::SingleTauDiscriminatorContainer,
0028                                       TauDiscriminator>::TauDiscriminationProducerBase(cfg),
0029         moduleLabel_(cfg.getParameter<std::string>("@module_label")),
0030         mva_(
0031             std::make_unique<AntiElectronIDMVA6<TauType, ElectronType>>(cfg, edm::EDConsumerBase::consumesCollector())),
0032         Electron_token(edm::EDConsumerBase::consumes<ElectronCollection>(
0033             cfg.getParameter<edm::InputTag>("srcElectrons"))),  // MB: full specification with prefix mandatory
0034         positionAtECalEntrance_(PositionAtECalEntranceComputer(edm::EDConsumerBase::consumesCollector(),
0035                                                                cfg.getParameter<bool>("isPhase2"))),
0036         vetoEcalCracks_(cfg.getParameter<bool>("vetoEcalCracks")),
0037         isPhase2_(cfg.getParameter<bool>("isPhase2")),
0038         verbosity_(cfg.getParameter<int>("verbosity")) {
0039     deltaREleTauMax_ = (isPhase2_ ? 0.2 : 0.3);
0040   }
0041 
0042   void beginEvent(const edm::Event& evt, const edm::EventSetup& es) override {
0043     mva_->beginEvent(evt, es);
0044     positionAtECalEntrance_.beginEvent(es);
0045     evt.getByToken(this->Tau_token, taus_);
0046     evt.getByToken(Electron_token, electrons_);
0047   }
0048 
0049   reco::SingleTauDiscriminatorContainer discriminate(const TauRef&) const override;
0050 
0051   ~TauDiscriminationAgainstElectronMVA6() override {}
0052 
0053   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0054 
0055 private:
0056   bool isInEcalCrack(double) const;
0057 
0058   // Overloaded method with explicit type specification to avoid partial
0059   //implementation of full class
0060   std::pair<float, float> getTauEtaAtECalEntrance(const reco::PFTauRef& theTauRef) const;
0061   std::pair<float, float> getTauEtaAtECalEntrance(const pat::TauRef& theTauRef) const;
0062 
0063   std::string moduleLabel_;
0064   std::unique_ptr<AntiElectronIDMVA6<TauType, ElectronType>> mva_;
0065 
0066   edm::EDGetTokenT<ElectronCollection> Electron_token;
0067   edm::Handle<ElectronCollection> electrons_;
0068   edm::Handle<TauCollection> taus_;
0069 
0070   PositionAtECalEntranceComputer positionAtECalEntrance_;
0071 
0072   static constexpr float ecalBarrelEndcapEtaBorder_ = 1.479;
0073   static constexpr float ecalEndcapVFEndcapEtaBorder_ = 2.4;
0074 
0075   bool vetoEcalCracks_;
0076 
0077   bool isPhase2_;
0078   float deltaREleTauMax_;
0079 
0080   int verbosity_;
0081 };
0082 
0083 template <class TauType, class TauDiscriminator, class ElectronType>
0084 reco::SingleTauDiscriminatorContainer
0085 TauDiscriminationAgainstElectronMVA6<TauType, TauDiscriminator, ElectronType>::discriminate(
0086     const TauRef& theTauRef) const {
0087   reco::SingleTauDiscriminatorContainer result;
0088   result.rawValues = {1., -1.};
0089   double category = -1.;
0090   bool isGsfElectronMatched = false;
0091 
0092   double deltaRDummy = 9.9;
0093 
0094   std::pair<float, float> tauEtaAtECalEntrance;
0095   if (std::is_same<TauType, reco::PFTau>::value || std::is_same<TauType, pat::Tau>::value)
0096     tauEtaAtECalEntrance = getTauEtaAtECalEntrance(theTauRef);
0097   else
0098     throw cms::Exception("TauDiscriminationAgainstElectronMVA6")
0099         << "Unsupported TauType used. You must use either reco::PFTau or pat::Tau.";
0100 
0101   if ((*theTauRef).leadChargedHadrCand().isNonnull()) {
0102     int numSignalGammaCandsInSigCone = 0;
0103     double signalrad = std::clamp(3.0 / std::max(1.0, theTauRef->pt()), 0.05, 0.10);
0104     for (const auto& gamma : theTauRef->signalGammaCands()) {
0105       double dR = deltaR(gamma->p4(), theTauRef->leadChargedHadrCand()->p4());
0106       // pfGammas inside the tau signal cone
0107       if (dR < signalrad) {
0108         numSignalGammaCandsInSigCone += 1;
0109       }
0110     }
0111 
0112     bool hasGsfTrack = false;
0113     const reco::CandidatePtr& leadChCand = theTauRef->leadChargedHadrCand();
0114     if (leadChCand.isNonnull()) {
0115       if (isPhase2_) {
0116         //MB: for phase-2 has gsf-track reads lead charged cand is pf-electron
0117         hasGsfTrack = (std::abs(leadChCand->pdgId()) == 11);
0118       } else {
0119         const pat::PackedCandidate* packedLeadChCand = dynamic_cast<const pat::PackedCandidate*>(leadChCand.get());
0120         if (packedLeadChCand != nullptr) {
0121           hasGsfTrack = (std::abs(packedLeadChCand->pdgId()) == 11);
0122         } else {
0123           const reco::PFCandidate* pfLeadChCand = dynamic_cast<const reco::PFCandidate*>(leadChCand.get());
0124           //pfLeadChCand can not be a nullptr here as it would be imply taus not built either with PFCandidates or PackedCandidates
0125           hasGsfTrack = pfLeadChCand->gsfTrackRef().isNonnull();
0126         }
0127       }
0128     }
0129 
0130     // loop over the electrons
0131     size_t iElec = 0;
0132     for (const auto& theElectron : *electrons_) {
0133       edm::Ref<ElectronCollection> theElecRef(electrons_, iElec);
0134       iElec++;
0135       if (theElectron.pt() > 10.) {  // CV: only take electrons above some minimal energy/Pt into account...
0136         double deltaREleTau = deltaR(theElectron.p4(), theTauRef->p4());
0137         deltaRDummy = std::min(deltaREleTau, deltaRDummy);
0138         if (deltaREleTau < deltaREleTauMax_) {
0139           double mva_match = mva_->mvaValue(*theTauRef, theElecRef);
0140           if (!hasGsfTrack)
0141             hasGsfTrack = theElectron.gsfTrack().isNonnull();
0142 
0143           // veto taus that go to ECal crack
0144           if (vetoEcalCracks_ &&
0145               (isInEcalCrack(tauEtaAtECalEntrance.first) || isInEcalCrack(tauEtaAtECalEntrance.second))) {
0146             // add category index
0147             result.rawValues.at(1) = category;
0148             // return MVA output value
0149             result.rawValues.at(0) = -99.;
0150             return result;
0151           }
0152           // veto taus that go to ECal crack
0153 
0154           if (std::abs(tauEtaAtECalEntrance.first) < ecalBarrelEndcapEtaBorder_) {  // Barrel
0155             if (numSignalGammaCandsInSigCone == 0 && hasGsfTrack) {
0156               category = 5.;
0157             } else if (numSignalGammaCandsInSigCone >= 1 && hasGsfTrack) {
0158               category = 7.;
0159             }
0160           } else if (!isPhase2_ || std::abs(tauEtaAtECalEntrance.first) < ecalEndcapVFEndcapEtaBorder_) {  // Endcap
0161             if (numSignalGammaCandsInSigCone == 0 && hasGsfTrack) {
0162               category = 13.;
0163             } else if (numSignalGammaCandsInSigCone >= 1 && hasGsfTrack) {
0164               category = 15.;
0165             }
0166           } else {  // VeryForwardEndcap
0167             if (numSignalGammaCandsInSigCone == 0 && hasGsfTrack) {
0168               category = 14.;
0169             } else if (numSignalGammaCandsInSigCone >= 1 && hasGsfTrack) {
0170               category = 16.;
0171             }
0172           }
0173 
0174           result.rawValues.at(0) = std::min(result.rawValues.at(0), float(mva_match));
0175           isGsfElectronMatched = true;
0176         }  // deltaR < deltaREleTauMax_
0177       }    // electron pt > 10
0178     }      // end of loop over electrons
0179 
0180     if (!isGsfElectronMatched) {
0181       double mva_nomatch = mva_->mvaValue(*theTauRef);
0182 
0183       // veto taus that go to ECal crack
0184       if (vetoEcalCracks_ &&
0185           (isInEcalCrack(tauEtaAtECalEntrance.first) || isInEcalCrack(tauEtaAtECalEntrance.second))) {
0186         // add category index
0187         result.rawValues.at(1) = category;
0188         // return MVA output value
0189         result.rawValues.at(0) = -99.;
0190         return result;
0191       }
0192       // veto taus that go to ECal crack
0193 
0194       if (std::abs(tauEtaAtECalEntrance.first) < ecalBarrelEndcapEtaBorder_) {  // Barrel
0195         if (numSignalGammaCandsInSigCone == 0 && !hasGsfTrack) {
0196           category = 0.;
0197         } else if (numSignalGammaCandsInSigCone >= 1 && !hasGsfTrack) {
0198           category = 2.;
0199         }
0200       } else if (!isPhase2_ || std::abs(tauEtaAtECalEntrance.first) < ecalEndcapVFEndcapEtaBorder_) {  // Endcap
0201         if (numSignalGammaCandsInSigCone == 0 && !hasGsfTrack) {
0202           category = 8.;
0203         } else if (numSignalGammaCandsInSigCone >= 1 && !hasGsfTrack) {
0204           category = 10.;
0205         }
0206       } else {  // VeryForwardEndcap
0207         if (numSignalGammaCandsInSigCone == 0 && !hasGsfTrack) {
0208           category = 9.;
0209         } else if (numSignalGammaCandsInSigCone >= 1 && !hasGsfTrack) {
0210           category = 11.;
0211         }
0212       }
0213 
0214       result.rawValues.at(0) = std::min(result.rawValues.at(0), float(mva_nomatch));
0215     }
0216   }
0217 
0218   if (verbosity_) {
0219     edm::LogPrint(this->getTauTypeString() + "AgainstEleMVA6")
0220         << "<" + this->getTauTypeString() + "AgainstElectronMVA6::discriminate>:";
0221     edm::LogPrint(this->getTauTypeString() + "AgainstEleMVA6")
0222         << " tau: Pt = " << theTauRef->pt() << ", eta = " << theTauRef->eta() << ", phi = " << theTauRef->phi();
0223     edm::LogPrint(this->getTauTypeString() + "AgainstEleMVA6")
0224         << " deltaREleTau = " << deltaRDummy << ", isGsfElectronMatched = " << isGsfElectronMatched;
0225     edm::LogPrint(this->getTauTypeString() + "AgainstEleMVA6")
0226         << " #Prongs = " << theTauRef->signalChargedHadrCands().size();
0227     edm::LogPrint(this->getTauTypeString() + "AgainstEleMVA6")
0228         << " MVA = " << result.rawValues.at(0) << ", category = " << category;
0229   }
0230 
0231   // add category index
0232   result.rawValues.at(1) = category;
0233   // return MVA output value
0234   return result;
0235 }
0236 
0237 template <class TauType, class TauDiscriminator, class ElectronType>
0238 bool TauDiscriminationAgainstElectronMVA6<TauType, TauDiscriminator, ElectronType>::isInEcalCrack(double eta) const {
0239   double absEta = std::abs(eta);
0240   return (absEta > 1.460 && absEta < 1.558);
0241 }
0242 
0243 template <class TauType, class TauDiscriminator, class ElectronType>
0244 std::pair<float, float>
0245 TauDiscriminationAgainstElectronMVA6<TauType, TauDiscriminator, ElectronType>::getTauEtaAtECalEntrance(
0246     const reco::PFTauRef& theTauRef) const {
0247   float tauEtaAtECalEntrance = -99;
0248   float leadChargedCandEtaAtECalEntrance = -99;
0249   float sumEtaTimesEnergy = 0;
0250   float sumEnergy = 0;
0251   float leadChargedCandPt = -99;
0252 
0253   for (const auto& candidate : theTauRef->signalCands()) {
0254     float etaAtECalEntrance = candidate->eta();
0255     const reco::Track* track = nullptr;
0256     const reco::PFCandidate* pfCandidate = dynamic_cast<const reco::PFCandidate*>(candidate.get());
0257     if (pfCandidate != nullptr) {
0258       if (!isPhase2_ || std::abs(theTauRef->eta()) < ecalBarrelEndcapEtaBorder_) {  // ECal
0259         etaAtECalEntrance = pfCandidate->positionAtECALEntrance().eta();
0260       } else {  // HGCal
0261         bool success = false;
0262         reco::Candidate::Point posAtECal = positionAtECalEntrance_(candidate.get(), success);
0263         if (success) {
0264           etaAtECalEntrance = posAtECal.eta();
0265         }
0266       }
0267       if (pfCandidate->trackRef().isNonnull())
0268         track = pfCandidate->trackRef().get();
0269       else if (pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->innerTrack().isNonnull())
0270         track = pfCandidate->muonRef()->innerTrack().get();
0271       else if (pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->globalTrack().isNonnull())
0272         track = pfCandidate->muonRef()->globalTrack().get();
0273       else if (pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->outerTrack().isNonnull())
0274         track = pfCandidate->muonRef()->outerTrack().get();
0275       else if (pfCandidate->gsfTrackRef().isNonnull())
0276         track = pfCandidate->gsfTrackRef().get();
0277     } else {
0278       bool success = false;
0279       reco::Candidate::Point posAtECal = positionAtECalEntrance_(candidate.get(), success);
0280       if (success) {
0281         etaAtECalEntrance = posAtECal.eta();
0282       }
0283       track = candidate->bestTrack();
0284     }
0285     if (track != nullptr) {
0286       if (track->pt() > leadChargedCandPt) {
0287         leadChargedCandEtaAtECalEntrance = etaAtECalEntrance;
0288         leadChargedCandPt = track->pt();
0289       }
0290     }
0291     sumEtaTimesEnergy += etaAtECalEntrance * candidate->energy();
0292     sumEnergy += candidate->energy();
0293   }
0294   if (sumEnergy > 0.) {
0295     tauEtaAtECalEntrance = sumEtaTimesEnergy / sumEnergy;
0296   }
0297   return std::pair<float, float>(tauEtaAtECalEntrance, leadChargedCandEtaAtECalEntrance);
0298 }
0299 
0300 template <class TauType, class TauDiscriminator, class ElectronType>
0301 std::pair<float, float>
0302 TauDiscriminationAgainstElectronMVA6<TauType, TauDiscriminator, ElectronType>::getTauEtaAtECalEntrance(
0303     const pat::TauRef& theTauRef) const {
0304   if (!isPhase2_ || std::abs(theTauRef->eta()) < ecalBarrelEndcapEtaBorder_) {  // ECal
0305     return std::pair<float, float>(theTauRef->etaAtEcalEntrance(), theTauRef->etaAtEcalEntranceLeadChargedCand());
0306   } else {  // HGCal
0307     float tauEtaAtECalEntrance = -99;
0308     float leadChargedCandEtaAtECalEntrance = -99;
0309     float sumEtaTimesEnergy = 0.;
0310     float sumEnergy = 0.;
0311     float leadChargedCandPt = -99;
0312 
0313     for (const auto& candidate : theTauRef->signalCands()) {
0314       float etaAtECalEntrance = candidate->eta();
0315       bool success = false;
0316       reco::Candidate::Point posAtECal = positionAtECalEntrance_(candidate.get(), success);
0317       if (success) {
0318         etaAtECalEntrance = posAtECal.eta();
0319       }
0320       const reco::Track* track = candidate->bestTrack();
0321       if (track != nullptr) {
0322         if (track->pt() > leadChargedCandPt) {
0323           leadChargedCandEtaAtECalEntrance = etaAtECalEntrance;
0324           leadChargedCandPt = track->pt();
0325         }
0326       }
0327       sumEtaTimesEnergy += etaAtECalEntrance * candidate->energy();
0328       sumEnergy += candidate->energy();
0329     }
0330     if (sumEnergy > 0.) {
0331       tauEtaAtECalEntrance = sumEtaTimesEnergy / sumEnergy;
0332     }
0333     return std::pair<float, float>(tauEtaAtECalEntrance, leadChargedCandEtaAtECalEntrance);
0334   }
0335 }
0336 
0337 template <class TauType, class TauDiscriminator, class ElectronType>
0338 void TauDiscriminationAgainstElectronMVA6<TauType, TauDiscriminator, ElectronType>::fillDescriptions(
0339     edm::ConfigurationDescriptions& descriptions) {
0340   // {pfReco,pat}TauDiscriminationAgainstElectronMVA6
0341   edm::ParameterSetDescription desc;
0342 
0343   desc.add<std::string>("method", "BDTG");
0344   desc.add<bool>("loadMVAfromDB", true);
0345   desc.add<bool>("returnMVA", true);
0346 
0347   desc.add<std::string>("mvaName_NoEleMatch_woGwoGSF_BL", "gbr_NoEleMatch_woGwoGSF_BL");
0348   desc.add<std::string>("mvaName_NoEleMatch_wGwoGSF_BL", "gbr_NoEleMatch_wGwoGSF_BL");
0349   desc.add<std::string>("mvaName_woGwGSF_BL", "gbr_woGwGSF_BL");
0350   desc.add<std::string>("mvaName_wGwGSF_BL", "gbr_wGwGSF_BL");
0351   desc.add<std::string>("mvaName_NoEleMatch_woGwoGSF_EC", "gbr_NoEleMatch_woGwoGSF_EC");
0352   desc.add<std::string>("mvaName_NoEleMatch_wGwoGSF_EC", "gbr_NoEleMatch_wGwoGSF_EC");
0353   desc.add<std::string>("mvaName_woGwGSF_EC", "gbr_woGwGSF_EC");
0354   desc.add<std::string>("mvaName_wGwGSF_EC", "gbr_wGwGSF_EC");
0355 
0356   desc.add<double>("minMVANoEleMatchWOgWOgsfBL", 0.0);
0357   desc.add<double>("minMVANoEleMatchWgWOgsfBL", 0.0);
0358   desc.add<double>("minMVAWOgWgsfBL", 0.0);
0359   desc.add<double>("minMVAWgWgsfBL", 0.0);
0360   desc.add<double>("minMVANoEleMatchWOgWOgsfEC", 0.0);
0361   desc.add<double>("minMVANoEleMatchWgWOgsfEC", 0.0);
0362   desc.add<double>("minMVAWOgWgsfEC", 0.0);
0363   desc.add<double>("minMVAWgWgsfEC", 0.0);
0364 
0365   desc.ifValue(
0366       edm::ParameterDescription<bool>("isPhase2", false, true),
0367       // MB: "srcElectrons" present for both phase-2 and non-phase2 to have a non-empy case for default, i.e. isPhase2=false
0368       false >> (edm::ParameterDescription<edm::InputTag>("srcElectrons", edm::InputTag("fixme"), true)) or
0369           // The following used only for Phase2
0370           true >> (edm::ParameterDescription<edm::InputTag>("srcElectrons", edm::InputTag("fixme"), true) and
0371                    edm::ParameterDescription<std::string>("mvaName_wGwGSF_VFEC", "gbr_wGwGSF_VFEC", true) and
0372                    edm::ParameterDescription<std::string>("mvaName_woGwGSF_VFEC", "gbr_woGwGSF_VFEC", true) and
0373                    edm::ParameterDescription<std::string>(
0374                        "mvaName_NoEleMatch_wGwoGSF_VFEC", "gbr_NoEleMatch_wGwoGSF_VFEC", true) and
0375                    edm::ParameterDescription<std::string>(
0376                        "mvaName_NoEleMatch_woGwoGSF_VFEC", "gbr_NoEleMatch_woGwoGSF_VFEC", true) and
0377                    edm::ParameterDescription<double>("minMVAWOgWgsfVFEC", 0.0, true) and
0378                    edm::ParameterDescription<double>("minMVAWgWgsfVFEC", 0.0, true) and
0379                    edm::ParameterDescription<double>("minMVANoEleMatchWgWOgsfVFEC", 0.0, true) and
0380                    edm::ParameterDescription<double>("minMVANoEleMatchWOgWOgsfVFEC", 0.0, true)));
0381 
0382   // Relevant only for gsfElectrons for Phase2
0383   if (std::is_same<ElectronType, reco::GsfElectron>::value) {
0384     desc.add<std::vector<edm::InputTag>>("hgcalElectronIDs", std::vector<edm::InputTag>())
0385         ->setComment("Relevant only for Phase-2");
0386   }
0387   desc.add<bool>("vetoEcalCracks", true);
0388   desc.add<bool>("usePhiAtEcalEntranceExtrapolation", false);
0389   desc.add<int>("verbosity", 0);
0390 
0391   TauDiscriminationProducerBase<TauType,
0392                                 reco::TauDiscriminatorContainer,
0393                                 reco::SingleTauDiscriminatorContainer,
0394                                 TauDiscriminator>::fillProducerDescriptions(desc);  // inherited from the base-class
0395 
0396   descriptions.addWithDefaultLabel(desc);
0397 }
0398 
0399 typedef TauDiscriminationAgainstElectronMVA6<reco::PFTau, reco::PFTauDiscriminator, reco::GsfElectron>
0400     PFRecoTauDiscriminationAgainstElectronMVA6;
0401 typedef TauDiscriminationAgainstElectronMVA6<pat::Tau, pat::PATTauDiscriminator, pat::Electron>
0402     PATTauDiscriminationAgainstElectronMVA6;
0403 
0404 DEFINE_FWK_MODULE(PFRecoTauDiscriminationAgainstElectronMVA6);
0405 DEFINE_FWK_MODULE(PATTauDiscriminationAgainstElectronMVA6);
[ 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!