diff --git a/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py b/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py index 2c837a10c0e0d..c99a84567762e 100644 --- a/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py +++ b/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py @@ -382,6 +382,34 @@ def miniAOD_customizeCommon(process): era.toReplaceWith(process.deepTauIDTask, deepTauIDTaskNew_) + #-- Adding tauID against dead ECal towers to taus + process.hpsPFTauDiscriminationByDeadECALElectronRejectionForMiniAOD = \ + process.hpsPFTauDiscriminationByDeadECALElectronRejection.clone( + extrapolateToECalEntrance = True + ) + _makePatTausTaskWithDeadECalVeto = process.makePatTausTask.copy() + _makePatTausTaskWithDeadECalVeto.add( + process.hpsPFTauDiscriminationByDeadECALElectronRejectionForMiniAOD + ) + run2_miniAOD_devel.toReplaceWith( + process.makePatTausTask, _makePatTausTaskWithDeadECalVeto + ) + _withDeadEcalTauIDPs = cms.PSet( + process.patTaus.tauIDSources, + againstElectronDeadECAL = cms.InputTag("hpsPFTauDiscriminationByDeadECALElectronRejectionForMiniAOD") + ) + run2_miniAOD_devel.toModify(process.patTaus, + tauIDSources = _withDeadEcalTauIDPs) + #... and to boosted taus + run2_miniAOD_devel.toModify(process.hpsPFTauDiscriminationByDeadECALElectronRejectionBoosted, + extrapolateToECalEntrance = True) + _withDeadEcalTauIDBoostedPs = cms.PSet( + process.patTausBoosted.tauIDSources, + againstElectronDeadECAL = cms.InputTag("hpsPFTauDiscriminationByDeadECALElectronRejectionBoosted") + ) + run2_miniAOD_devel.toModify(process.patTausBoosted, + tauIDSources = _withDeadEcalTauIDBoostedPs) + #-- Adding customization for 80X 2016 legacy reMiniAOD and 2018 heavy ions from Configuration.Eras.Modifier_run2_miniAOD_80XLegacy_cff import run2_miniAOD_80XLegacy from Configuration.Eras.Modifier_pp_on_AA_2018_cff import pp_on_AA_2018 diff --git a/PhysicsTools/PatAlgos/python/tools/tauTools.py b/PhysicsTools/PatAlgos/python/tools/tauTools.py index e9898433a9d62..6d166ab47eeb6 100644 --- a/PhysicsTools/PatAlgos/python/tools/tauTools.py +++ b/PhysicsTools/PatAlgos/python/tools/tauTools.py @@ -175,6 +175,7 @@ def _switchToPFTau(process, ("againstElectronMediumMVA6", "DiscriminationByMVA6MediumElectronRejection"), ("againstElectronTightMVA6", "DiscriminationByMVA6TightElectronRejection"), ("againstElectronVTightMVA6", "DiscriminationByMVA6VTightElectronRejection"), + ("againstElectronDeadECAL", "DiscriminationByDeadECALElectronRejection"), ] # switch to PFTau collection produced for fixed dR = 0.07 signal cone size diff --git a/RecoTauTag/Configuration/python/HPSPFTaus_cff.py b/RecoTauTag/Configuration/python/HPSPFTaus_cff.py index f77ac903a0d44..8fe7b05724c00 100644 --- a/RecoTauTag/Configuration/python/HPSPFTaus_cff.py +++ b/RecoTauTag/Configuration/python/HPSPFTaus_cff.py @@ -13,7 +13,7 @@ from RecoTauTag.RecoTau.PFRecoTauDiscriminationByLeadingTrackFinding_cfi import * from RecoTauTag.RecoTau.PFRecoTauDiscriminationAgainstElectron_cfi import * from RecoTauTag.RecoTau.PFRecoTauDiscriminationAgainstElectronMVA6_cfi import * -from RecoTauTag.RecoTau.PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi import * +from RecoTauTag.RecoTau.pfRecoTauDiscriminationAgainstElectronDeadECAL_cfi import * from RecoTauTag.RecoTau.PFRecoTauDiscriminationAgainstMuon_cfi import * from RecoTauTag.RecoTau.PFRecoTauDiscriminationAgainstMuon2_cfi import * from RecoTauTag.RecoTau.PFRecoTauDiscriminationAgainstMuonMVA_cfi import * @@ -300,7 +300,8 @@ ## ByDeadECALElectronRejection hpsPFTauDiscriminationByDeadECALElectronRejection = pfRecoTauDiscriminationAgainstElectronDeadECAL.clone( PFTauProducer = cms.InputTag('hpsPFTauProducer'), - Prediscriminants = requireDecayMode.clone() + Prediscriminants = requireDecayMode.clone(), + extrapolateToECalEntrance = cms.bool(False) ) ## ByMVA6rawElectronRejection hpsPFTauDiscriminationByMVA6rawElectronRejection = pfRecoTauDiscriminationAgainstElectronMVA6.clone( diff --git a/RecoTauTag/RecoTau/BuildFile.xml b/RecoTauTag/RecoTau/BuildFile.xml index b21ab31b90aba..47921f0d590b0 100644 --- a/RecoTauTag/RecoTau/BuildFile.xml +++ b/RecoTauTag/RecoTau/BuildFile.xml @@ -1,12 +1,16 @@ + + + + diff --git a/RecoTauTag/RecoTau/interface/AntiElectronDeadECAL.h b/RecoTauTag/RecoTau/interface/AntiElectronDeadECAL.h new file mode 100644 index 0000000000000..d37bb1b73b5bd --- /dev/null +++ b/RecoTauTag/RecoTau/interface/AntiElectronDeadECAL.h @@ -0,0 +1,74 @@ +#ifndef RecoTauTag_RecoTau_AntiElectronDeadECAL_h +#define RecoTauTag_RecoTau_AntiElectronDeadECAL_h + +/** \class AntiElectronDeadECAL + * + * Flag tau candidates reconstructed near dead ECAL channels, + * in order to reduce e -> tau fakes not rejected by anti-e MVA discriminator + * + * The motivation for this flag is this presentation: + * https://indico.cern.ch/getFile.py/access?contribId=0&resId=0&materialId=slides&confId=177223 + * + * Code adapted from: + * RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationAgainstElectronDeadECAL.cc + * + * \authors Lauri Andreas Wendland, + * Christian Veelken + * + * + * + */ + +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/ESWatcher.h" +#include "DataFormats/Candidate/interface/Candidate.h" +#include "RecoTauTag/RecoTau/interface/PositionAtECalEntranceComputer.h" + +#include +#include + +class EcalChannelStatusRcd; +class CaloGeometryRecord; +class IdealGeometryRecord; + +class AntiElectronDeadECAL { +public: + explicit AntiElectronDeadECAL(const edm::ParameterSet&); + ~AntiElectronDeadECAL(); + + void beginEvent(const edm::EventSetup&); + + bool operator()(const reco::Candidate* tau) const; + +private: + const unsigned minStatus_; + const double dR2_; + const bool extrapolateToECalEntrance_; + const int verbosity_; + + PositionAtECalEntranceComputer positionAtECalEntrance_; + + void updateBadTowers(const edm::EventSetup&); + + struct TowerInfo { + TowerInfo(uint32_t id, unsigned nBad, unsigned maxStatus, double eta, double phi) + : id_(id), nBad_(nBad), maxStatus_(maxStatus), eta_(eta), phi_(phi) {} + uint32_t id_; + unsigned nBad_; + unsigned maxStatus_; + double eta_; + double phi_; + }; + + std::vector badTowers_; + static const uint16_t statusMask_ = 0x1F; + + edm::ESWatcher channelStatusWatcher_; + edm::ESWatcher caloGeometryWatcher_; + edm::ESWatcher idealGeometryWatcher_; + + bool isFirstEvent_; +}; + +#endif // RecoTauTag_RecoTau_AntiElectronDeadECAL_h diff --git a/RecoTauTag/RecoTau/interface/AntiElectronIDMVA6.h b/RecoTauTag/RecoTau/interface/AntiElectronIDMVA6.h index fe8e9f8232711..0a6473d62d6ad 100644 --- a/RecoTauTag/RecoTau/interface/AntiElectronIDMVA6.h +++ b/RecoTauTag/RecoTau/interface/AntiElectronIDMVA6.h @@ -22,9 +22,8 @@ #include "CondFormats/EgammaObjects/interface/GBRForest.h" #include "DataFormats/PatCandidates/interface/Tau.h" #include "DataFormats/PatCandidates/interface/Electron.h" - -#include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h" #include "DataFormats/PatCandidates/interface/PackedCandidate.h" +#include "RecoTauTag/RecoTau/interface/PositionAtECalEntranceComputer.h" #include "TMVA/Tools.h" #include "TMVA/Reader.h" @@ -33,150 +32,144 @@ #include -class AntiElectronIDMVA6 -{ - public: - - AntiElectronIDMVA6(const edm::ParameterSet&); - ~AntiElectronIDMVA6(); - - void beginEvent(const edm::Event&, const edm::EventSetup&); - - double MVAValue(Float_t TauPt, - Float_t TauEtaAtEcalEntrance, - Float_t TauPhi, - Float_t TauLeadChargedPFCandPt, - Float_t TauLeadChargedPFCandEtaAtEcalEntrance, - Float_t TauEmFraction, - Float_t TauLeadPFChargedHadrHoP, - Float_t TauLeadPFChargedHadrEoP, - Float_t TauVisMassIn, - Float_t TaudCrackEta, - Float_t TaudCrackPhi, - Float_t TauHasGsf, - Int_t TauSignalPFGammaCandsIn, - Int_t TauSignalPFGammaCandsOut, - const std::vector& GammasdEtaInSigCone, - const std::vector& GammasdPhiInSigCone, - const std::vector& GammasPtInSigCone, - const std::vector& GammasdEtaOutSigCone, - const std::vector& GammasdPhiOutSigCone, - const std::vector& GammasPtOutSigCone, - Float_t ElecEta, - Float_t ElecPhi, - Float_t ElecEtotOverPin, - Float_t ElecChi2NormGSF, - Float_t ElecChi2NormKF, - Float_t ElecGSFNumHits, - Float_t ElecKFNumHits, - Float_t ElecGSFTrackResol, - Float_t ElecGSFTracklnPt, - Float_t ElecPin, - Float_t ElecPout, - Float_t ElecEecal, - Float_t ElecDeltaEta, - Float_t ElecDeltaPhi, - Float_t ElecMvaInSigmaEtaEta, - Float_t ElecMvaInHadEnergy, - Float_t ElecMvaInDeltaEta - ); - - double MVAValue(Float_t TauPt, - Float_t TauEtaAtEcalEntrance, - Float_t TauPhi, - Float_t TauLeadChargedPFCandPt, - Float_t TauLeadChargedPFCandEtaAtEcalEntrance, - Float_t TauEmFraction, - Float_t TauLeadPFChargedHadrHoP, - Float_t TauLeadPFChargedHadrEoP, - Float_t TauVisMassIn, - Float_t TaudCrackEta, - Float_t TaudCrackPhi, - Float_t TauHasGsf, - Int_t TauSignalPFGammaCandsIn, - Int_t TauSignalPFGammaCandsOut, - Float_t TauGammaEtaMomIn, - Float_t TauGammaEtaMomOut, - Float_t TauGammaPhiMomIn, - Float_t TauGammaPhiMomOut, - Float_t TauGammaEnFracIn, - Float_t TauGammaEnFracOut, - Float_t ElecEta, - Float_t ElecPhi, - Float_t ElecEtotOverPin, - Float_t ElecChi2NormGSF, - Float_t ElecChi2NormKF, - Float_t ElecGSFNumHits, - Float_t ElecKFNumHits, - Float_t ElecGSFTrackResol, - Float_t ElecGSFTracklnPt, - Float_t ElecPin, - Float_t ElecPout, - Float_t ElecEecal, - Float_t ElecDeltaEta, - Float_t ElecDeltaPhi, - Float_t ElecMvaInSigmaEtaEta, - Float_t ElecMvaInHadEnergy, - Float_t ElecMvaInDeltaEta - ); - - // this function can be called for all categories - double MVAValue(const reco::PFTau& thePFTau, - const reco::GsfElectron& theGsfEle); - // this function can be called for category 1 only !! - double MVAValue(const reco::PFTau& thePFTau); - - // this function can be called for all categories - double MVAValue(const pat::Tau& theTau, - const pat::Electron& theEle); - // this function can be called for category 1 only !! - double MVAValue(const pat::Tau& theTau); - // track extrapolation to ECAL entrance (used to re-calculate varibales that might not be available on miniAOD) - bool atECalEntrance(const reco::Candidate* part, math::XYZPoint &pos); - - private: - - double dCrackEta(double eta); - double minimum(double a,double b); - double dCrackPhi(double phi, double eta); - - bool isInitialized_; - bool loadMVAfromDB_; - edm::FileInPath inputFileName_; - - std::string mvaName_NoEleMatch_woGwoGSF_BL_; - std::string mvaName_NoEleMatch_wGwoGSF_BL_; - std::string mvaName_woGwGSF_BL_; - std::string mvaName_wGwGSF_BL_; - std::string mvaName_NoEleMatch_woGwoGSF_EC_; - std::string mvaName_NoEleMatch_wGwoGSF_EC_; - std::string mvaName_woGwGSF_EC_; - std::string mvaName_wGwGSF_EC_; - - bool usePhiAtEcalEntranceExtrapolation_; - - Float_t* Var_NoEleMatch_woGwoGSF_Barrel_; - Float_t* Var_NoEleMatch_wGwoGSF_Barrel_; - Float_t* Var_woGwGSF_Barrel_; - Float_t* Var_wGwGSF_Barrel_; - Float_t* Var_NoEleMatch_woGwoGSF_Endcap_; - Float_t* Var_NoEleMatch_wGwoGSF_Endcap_; - Float_t* Var_woGwGSF_Endcap_; - Float_t* Var_wGwGSF_Endcap_; - - const GBRForest* mva_NoEleMatch_woGwoGSF_BL_; - const GBRForest* mva_NoEleMatch_wGwoGSF_BL_; - const GBRForest* mva_woGwGSF_BL_; - const GBRForest* mva_wGwGSF_BL_; - const GBRForest* mva_NoEleMatch_woGwoGSF_EC_; - const GBRForest* mva_NoEleMatch_wGwoGSF_EC_; - const GBRForest* mva_woGwGSF_EC_; - const GBRForest* mva_wGwGSF_EC_; - - std::vector inputFilesToDelete_; - - double bField_; - int verbosity_; +class AntiElectronIDMVA6 { +public: + AntiElectronIDMVA6(const edm::ParameterSet&); + ~AntiElectronIDMVA6(); + + void beginEvent(const edm::Event&, const edm::EventSetup&); + + double MVAValue(Float_t TauPt, + Float_t TauEtaAtEcalEntrance, + Float_t TauPhi, + Float_t TauLeadChargedPFCandPt, + Float_t TauLeadChargedPFCandEtaAtEcalEntrance, + Float_t TauEmFraction, + Float_t TauLeadPFChargedHadrHoP, + Float_t TauLeadPFChargedHadrEoP, + Float_t TauVisMassIn, + Float_t TaudCrackEta, + Float_t TaudCrackPhi, + Float_t TauHasGsf, + Int_t TauSignalPFGammaCandsIn, + Int_t TauSignalPFGammaCandsOut, + const std::vector& GammasdEtaInSigCone, + const std::vector& GammasdPhiInSigCone, + const std::vector& GammasPtInSigCone, + const std::vector& GammasdEtaOutSigCone, + const std::vector& GammasdPhiOutSigCone, + const std::vector& GammasPtOutSigCone, + Float_t ElecEta, + Float_t ElecPhi, + Float_t ElecEtotOverPin, + Float_t ElecChi2NormGSF, + Float_t ElecChi2NormKF, + Float_t ElecGSFNumHits, + Float_t ElecKFNumHits, + Float_t ElecGSFTrackResol, + Float_t ElecGSFTracklnPt, + Float_t ElecPin, + Float_t ElecPout, + Float_t ElecEecal, + Float_t ElecDeltaEta, + Float_t ElecDeltaPhi, + Float_t ElecMvaInSigmaEtaEta, + Float_t ElecMvaInHadEnergy, + Float_t ElecMvaInDeltaEta); + + double MVAValue(Float_t TauPt, + Float_t TauEtaAtEcalEntrance, + Float_t TauPhi, + Float_t TauLeadChargedPFCandPt, + Float_t TauLeadChargedPFCandEtaAtEcalEntrance, + Float_t TauEmFraction, + Float_t TauLeadPFChargedHadrHoP, + Float_t TauLeadPFChargedHadrEoP, + Float_t TauVisMassIn, + Float_t TaudCrackEta, + Float_t TaudCrackPhi, + Float_t TauHasGsf, + Int_t TauSignalPFGammaCandsIn, + Int_t TauSignalPFGammaCandsOut, + Float_t TauGammaEtaMomIn, + Float_t TauGammaEtaMomOut, + Float_t TauGammaPhiMomIn, + Float_t TauGammaPhiMomOut, + Float_t TauGammaEnFracIn, + Float_t TauGammaEnFracOut, + Float_t ElecEta, + Float_t ElecPhi, + Float_t ElecEtotOverPin, + Float_t ElecChi2NormGSF, + Float_t ElecChi2NormKF, + Float_t ElecGSFNumHits, + Float_t ElecKFNumHits, + Float_t ElecGSFTrackResol, + Float_t ElecGSFTracklnPt, + Float_t ElecPin, + Float_t ElecPout, + Float_t ElecEecal, + Float_t ElecDeltaEta, + Float_t ElecDeltaPhi, + Float_t ElecMvaInSigmaEtaEta, + Float_t ElecMvaInHadEnergy, + Float_t ElecMvaInDeltaEta); + + // this function can be called for all categories + double MVAValue(const reco::PFTau& thePFTau, const reco::GsfElectron& theGsfEle); + // this function can be called for category 1 only !! + double MVAValue(const reco::PFTau& thePFTau); + + // this function can be called for all categories + double MVAValue(const pat::Tau& theTau, const pat::Electron& theEle); + // this function can be called for category 1 only !! + double MVAValue(const pat::Tau& theTau); + // track extrapolation to ECAL entrance (used to re-calculate variables that might not be available on miniAOD) + bool atECalEntrance(const reco::Candidate* part, math::XYZPoint& pos); + +private: + double dCrackEta(double eta); + double minimum(double a, double b); + double dCrackPhi(double phi, double eta); + + bool isInitialized_; + bool loadMVAfromDB_; + edm::FileInPath inputFileName_; + + std::string mvaName_NoEleMatch_woGwoGSF_BL_; + std::string mvaName_NoEleMatch_wGwoGSF_BL_; + std::string mvaName_woGwGSF_BL_; + std::string mvaName_wGwGSF_BL_; + std::string mvaName_NoEleMatch_woGwoGSF_EC_; + std::string mvaName_NoEleMatch_wGwoGSF_EC_; + std::string mvaName_woGwGSF_EC_; + std::string mvaName_wGwGSF_EC_; + + bool usePhiAtEcalEntranceExtrapolation_; + + Float_t* Var_NoEleMatch_woGwoGSF_Barrel_; + Float_t* Var_NoEleMatch_wGwoGSF_Barrel_; + Float_t* Var_woGwGSF_Barrel_; + Float_t* Var_wGwGSF_Barrel_; + Float_t* Var_NoEleMatch_woGwoGSF_Endcap_; + Float_t* Var_NoEleMatch_wGwoGSF_Endcap_; + Float_t* Var_woGwGSF_Endcap_; + Float_t* Var_wGwGSF_Endcap_; + + const GBRForest* mva_NoEleMatch_woGwoGSF_BL_; + const GBRForest* mva_NoEleMatch_wGwoGSF_BL_; + const GBRForest* mva_woGwGSF_BL_; + const GBRForest* mva_wGwGSF_BL_; + const GBRForest* mva_NoEleMatch_woGwoGSF_EC_; + const GBRForest* mva_NoEleMatch_wGwoGSF_EC_; + const GBRForest* mva_woGwGSF_EC_; + const GBRForest* mva_wGwGSF_EC_; + + std::vector inputFilesToDelete_; + + PositionAtECalEntranceComputer positionAtECalEntrance_; + + int verbosity_; }; #endif diff --git a/RecoTauTag/RecoTau/interface/PositionAtECalEntranceComputer.h b/RecoTauTag/RecoTau/interface/PositionAtECalEntranceComputer.h new file mode 100644 index 0000000000000..0b167c2568549 --- /dev/null +++ b/RecoTauTag/RecoTau/interface/PositionAtECalEntranceComputer.h @@ -0,0 +1,33 @@ +#ifndef RecoTauTag_RecoTau_PositionAtECalEntranceComputer_h +#define RecoTauTag_RecoTau_PositionAtECalEntranceComputer_h + +/** \class PositionAtECalEntranceComputer + * + * Extrapolate particle (charged or neutral) to ECAL entrance, + * in order to compute the distance of the tau to ECAL cracks and/or dead ECAL channels + * + * \authors Fabio Colombo, + * Christian Veelken + * + * + * + */ + +#include "FWCore/Framework/interface/EventSetup.h" +#include "DataFormats/Candidate/interface/Candidate.h" + +class PositionAtECalEntranceComputer { +public: + PositionAtECalEntranceComputer(); + ~PositionAtECalEntranceComputer(); + + void beginEvent(const edm::EventSetup&); + + //To do: it seems to more practical to put this to the ES + reco::Candidate::Point operator()(const reco::Candidate* particle, bool& success) const; + +private: + double bField_z_; +}; + +#endif // RecoTauTag_RecoTau_PositionAtECalEntranceComputer_h diff --git a/RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h b/RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h index 373c29ca4994c..4c3586521319d 100644 --- a/RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h +++ b/RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h @@ -92,6 +92,10 @@ class TauDiscriminationProducerBase : public edm::stream::EDProducer<> { static void fillProducerDescriptions(edm::ParameterSetDescription& desc); + /// helper method to retrieve tau type name, e.g. to build correct cfi getter + //string (i.e. PFTau/PATTauProducer) + static std::string getTauTypeString(); + protected: //value given to taus that fail prediscriminants double prediscriminantFailValue_; @@ -118,14 +122,4 @@ typedef TauDiscriminationProducerBase typedef TauDiscriminationProducerBase CaloTauDiscriminationProducerBase; - -/// helper function retrieve the correct cfi getter string (ie PFTauProducer) -//for this tau type -template std::string getProducerString() -{ - // this generic one shoudl never be called. - // these are specialized in TauDiscriminationProducerBase.cc - throw cms::Exception("TauDiscriminationProducerBase") - << "Unsupported TauType used. You must use either PFTau, PATTau or CaloTaus."; -} #endif diff --git a/RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationAgainstElectronDeadECAL.cc b/RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationAgainstElectronDeadECAL.cc deleted file mode 100644 index f3c1fe6168058..0000000000000 --- a/RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationAgainstElectronDeadECAL.cc +++ /dev/null @@ -1,213 +0,0 @@ - -/** \class PFRecoTauDiscriminationAgainstElectronDeadECAL - * - * Flag tau candidates reconstructed near dead ECAL channels, - * in order to reduce e -> tau fakes not rejected by anti-e MVA discriminator - * - * The motivation for this flag is this presentation: - * https://indico.cern.ch/getFile.py/access?contribId=0&resId=0&materialId=slides&confId=177223 - * - * \authors Lauri Andreas Wendland, - * Christian Veelken - * - * - * - */ - -#include "RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h" -#include "FWCore/Framework/interface/ESHandle.h" -#include -#include - -#include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h" -#include "CondFormats/EcalObjects/interface/EcalChannelStatus.h" -#include "Geometry/Records/interface/CaloGeometryRecord.h" -#include "Geometry/CaloGeometry/interface/CaloGeometry.h" -#include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h" -#include "Geometry/Records/interface/IdealGeometryRecord.h" -#include "DataFormats/DetId/interface/DetId.h" -#include "DataFormats/EcalDetId/interface/EBDetId.h" -#include "DataFormats/EcalDetId/interface/EEDetId.h" -#include "DataFormats/Math/interface/deltaR.h" -#include "DataFormats/TauReco/interface/PFTauDiscriminator.h" - -#include - -using namespace reco; - -class PFRecoTauDiscriminationAgainstElectronDeadECAL : public PFTauDiscriminationProducerBase -{ - public: - explicit PFRecoTauDiscriminationAgainstElectronDeadECAL(const edm::ParameterSet& cfg) - : PFTauDiscriminationProducerBase(cfg), - moduleLabel_(cfg.getParameter("@module_label")), - isFirstEvent_(true) - { - minStatus_ = cfg.getParameter("minStatus"); - dR_ = cfg.getParameter("dR"); - - verbosity_ = cfg.getParameter("verbosity"); - } - ~PFRecoTauDiscriminationAgainstElectronDeadECAL() override {} - - void beginEvent(const edm::Event& evt, const edm::EventSetup& es) override - { - updateBadTowers(es); - } - - double discriminate(const PFTauRef& pfTau) const override - { - if ( verbosity_ ) { - edm::LogPrint("PFTauAgainstEleDeadECAL") << ":" ; - edm::LogPrint("PFTauAgainstEleDeadECAL") << " moduleLabel = " << moduleLabel_ ; - edm::LogPrint("PFTauAgainstEleDeadECAL") << "#badTowers = " << badTowers_.size() ; - edm::LogPrint("PFTauAgainstEleDeadECAL") << "tau: Pt = " << pfTau->pt() << ", eta = " << pfTau->eta() << ", phi = " << pfTau->phi() ; - } - double discriminator = 1.; - for ( std::vector::const_iterator badTower = badTowers_.begin(); - badTower != badTowers_.end(); ++badTower ) { - if ( deltaR(badTower->eta_, badTower->phi_, pfTau->eta(), pfTau->phi()) < dR_ ) { - if ( verbosity_ ) { - edm::LogPrint("PFTauAgainstEleDeadECAL") << " matches badTower: eta = " << badTower->eta_ << ", phi = " << badTower->phi_ ; - } - discriminator = 0.; - } - } - if ( verbosity_ ) { - edm::LogPrint("PFTauAgainstEleDeadECAL") << "--> discriminator = " << discriminator ; - } - return discriminator; - } - - static void fillDescriptions(edm::ConfigurationDescriptions & descriptions); - - private: - void updateBadTowers(const edm::EventSetup& es) - { - // NOTE: modified version of SUSY CAF code - // UserCode/SusyCAF/plugins/SusyCAF_EcalDeadChannels.cc - const uint32_t channelStatusId = es.get().cacheIdentifier(); - const uint32_t caloGeometryId = es.get().cacheIdentifier(); - const uint32_t idealGeometryId = es.get().cacheIdentifier(); - - if ( !isFirstEvent_ && channelStatusId == channelStatusId_cache_ && caloGeometryId == caloGeometryId_cache_ && idealGeometryId == idealGeometryId_cache_ ) return; - - edm::ESHandle channelStatus; - es.get().get(channelStatus); - channelStatusId_cache_ = channelStatusId; - - edm::ESHandle caloGeometry; - es.get().get(caloGeometry); - caloGeometryId_cache_ = caloGeometryId; - - edm::ESHandle ttMap; - es.get().get(ttMap); - idealGeometryId_cache_ = idealGeometryId; - - std::map nBadCrystals, maxStatus; - std::map sumEta, sumPhi; - - loopXtals(nBadCrystals, maxStatus, sumEta, sumPhi, channelStatus.product(), caloGeometry.product(), ttMap.product()); - loopXtals(nBadCrystals, maxStatus, sumEta, sumPhi, channelStatus.product(), caloGeometry.product(), ttMap.product()); - - badTowers_.clear(); - for ( std::map::const_iterator it = nBadCrystals.begin(); - it != nBadCrystals.end(); ++it ) { - uint32_t key = it->first; - badTowers_.push_back(towerInfo(key, it->second, maxStatus[key], sumEta[key]/it->second, sumPhi[key]/it->second)); - } - - isFirstEvent_ = false; - } - - template - void loopXtals(std::map& nBadCrystals, - std::map& maxStatus, - std::map& sumEta, - std::map& sumPhi , - const EcalChannelStatus* channelStatus, - const CaloGeometry* caloGeometry, - const EcalTrigTowerConstituentsMap* ttMap) const - { - // NOTE: modified version of SUSY CAF code - // UserCode/SusyCAF/plugins/SusyCAF_EcalDeadChannels.cc - for ( int i = 0; i < Id::kSizeForDenseIndexing; ++i ) { - Id id = Id::unhashIndex(i); - if ( id == Id(0) ) continue; - EcalChannelStatusMap::const_iterator it = channelStatus->getMap().find(id.rawId()); - unsigned status = ( it == channelStatus->end() ) ? - 0 : (it->getStatusCode() & statusMask_); - if ( status >= minStatus_ ) { - const GlobalPoint& point = caloGeometry->getPosition(id); - uint32_t key = ttMap->towerOf(id); - maxStatus[key] = TMath::Max(status, maxStatus[key]); - ++nBadCrystals[key]; - sumEta[key] += point.eta(); - sumPhi[key] += point.phi(); - } - } - } - - struct towerInfo - { - towerInfo(uint32_t id, unsigned nBad, unsigned maxStatus, double eta, double phi) - : id_(id), - nBad_(nBad), - maxStatus_(maxStatus), - eta_(eta), - phi_(phi) - {} - uint32_t id_; - unsigned nBad_; - unsigned maxStatus_; - double eta_; - double phi_; - }; - typedef ROOT::Math::LorentzVector > PolarLorentzVector; - - std::string moduleLabel_; - unsigned minStatus_; - double dR_; - - std::vector badTowers_; - static const uint16_t statusMask_ = 0x1F; - - uint32_t channelStatusId_cache_; - uint32_t caloGeometryId_cache_; - uint32_t idealGeometryId_cache_; - bool isFirstEvent_; - - int verbosity_; -}; - -void -PFRecoTauDiscriminationAgainstElectronDeadECAL::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { - // pfRecoTauDiscriminationAgainstElectronDeadECAL - edm::ParameterSetDescription desc; - desc.add("verbosity", 0); - { - edm::ParameterSetDescription pset_Prediscriminants; - pset_Prediscriminants.add("BooleanOperator", "and"); - { - edm::ParameterSetDescription psd1; - psd1.add("cut"); - psd1.add("Producer"); - pset_Prediscriminants.addOptional("leadTrack", psd1); - } - { - // encountered this at - // RecoTauTag/Configuration/python/HPSPFTaus_cff.py - edm::ParameterSetDescription psd1; - psd1.add("cut"); - psd1.add("Producer"); - pset_Prediscriminants.addOptional("decayMode", psd1); - } - desc.add("Prediscriminants", pset_Prediscriminants); - } - desc.add("dR", 0.08); - desc.add("PFTauProducer", edm::InputTag("pfTauProducer")); - desc.add("minStatus", 12); - descriptions.add("pfRecoTauDiscriminationAgainstElectronDeadECAL", desc); -} - -DEFINE_FWK_MODULE(PFRecoTauDiscriminationAgainstElectronDeadECAL); diff --git a/RecoTauTag/RecoTau/plugins/TauDiscriminationAgainstElectronDeadECAL.cc b/RecoTauTag/RecoTau/plugins/TauDiscriminationAgainstElectronDeadECAL.cc new file mode 100644 index 0000000000000..578a232bb1d67 --- /dev/null +++ b/RecoTauTag/RecoTau/plugins/TauDiscriminationAgainstElectronDeadECAL.cc @@ -0,0 +1,87 @@ +/** \class TauDiscriminationAgainstElectronDeadECAL + * + * Template class for producing PFTau and PATTau discriminators which + * flag tau candidates reconstructed near dead ECAL channels, + * in order to reduce e -> tau fakes not rejected by anti-e MVA discriminator + * + * The motivation for this flag is this presentation: + * https://indico.cern.ch/getFile.py/access?contribId=0&resId=0&materialId=slides&confId=177223 + * + * \authors Lauri Andreas Wendland, + * Christian Veelken + * + * + * + */ +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h" +#include "RecoTauTag/RecoTau/interface/AntiElectronDeadECAL.h" + +template +class TauDiscriminationAgainstElectronDeadECAL : public TauDiscriminationProducerBase { +public: + typedef std::vector TauCollection; + typedef edm::Ref TauRef; + explicit TauDiscriminationAgainstElectronDeadECAL(const edm::ParameterSet& cfg) + : TauDiscriminationProducerBase::TauDiscriminationProducerBase(cfg), + moduleLabel_(cfg.getParameter("@module_label")), + verbosity_(cfg.getParameter("verbosity")), + antiElectronDeadECAL_(cfg) {} + ~TauDiscriminationAgainstElectronDeadECAL() override {} + + void beginEvent(const edm::Event& evt, const edm::EventSetup& es) override { antiElectronDeadECAL_.beginEvent(es); } + + double discriminate(const TauRef& tau) const override { + if (verbosity_) { + edm::LogPrint(this->getTauTypeString() + "AgainstEleDeadECAL") + << "<" + this->getTauTypeString() + "AgainstElectronDeadECAL::discriminate>:"; + edm::LogPrint(this->getTauTypeString() + "AgainstEleDeadECAL") << " moduleLabel = " << moduleLabel_; + edm::LogPrint(this->getTauTypeString() + "AgainstEleDeadECAL") + << " tau: Pt = " << tau->pt() << ", eta = " << tau->eta() << ", phi = " << tau->phi(); + } + double discriminator = 1.; + if (antiElectronDeadECAL_(tau.get())) { + discriminator = 0.; + } + if (verbosity_) { + edm::LogPrint(this->getTauTypeString() + "AgainstEleDeadECAL") << "--> discriminator = " << discriminator; + } + return discriminator; + } + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + std::string moduleLabel_; + int verbosity_; + + AntiElectronDeadECAL antiElectronDeadECAL_; +}; + +template +void TauDiscriminationAgainstElectronDeadECAL::fillDescriptions( + edm::ConfigurationDescriptions& descriptions) { + // {pfReco,pat}TauDiscriminationAgainstElectronDeadECAL + edm::ParameterSetDescription desc; + + desc.add("dR", 0.08); + desc.add("minStatus", 12); + desc.add("extrapolateToECalEntrance", true); + desc.add("verbosity", 0); + + TauDiscriminationProducerBase::fillProducerDescriptions( + desc); // inherited from the base-class + + descriptions.addWithDefaultLabel(desc); +} + +typedef TauDiscriminationAgainstElectronDeadECAL + PFRecoTauDiscriminationAgainstElectronDeadECAL; +typedef TauDiscriminationAgainstElectronDeadECAL + PATTauDiscriminationAgainstElectronDeadECAL; + +DEFINE_FWK_MODULE(PFRecoTauDiscriminationAgainstElectronDeadECAL); +DEFINE_FWK_MODULE(PATTauDiscriminationAgainstElectronDeadECAL); diff --git a/RecoTauTag/RecoTau/python/PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.py b/RecoTauTag/RecoTau/python/PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.py deleted file mode 100644 index f891f28f305d2..0000000000000 --- a/RecoTauTag/RecoTau/python/PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.py +++ /dev/null @@ -1,21 +0,0 @@ -import FWCore.ParameterSet.Config as cms -from RecoTauTag.RecoTau.TauDiscriminatorTools import requireLeadTrack - -pfRecoTauDiscriminationAgainstElectronDeadECAL = cms.EDProducer( - "PFRecoTauDiscriminationAgainstElectronDeadECAL", - - # tau collection to discriminate - PFTauProducer = cms.InputTag('pfTauProducer'), - - # Require leading pion ensures that: - # 1) these is at least one track above threshold (0.5 GeV) in the signal cone - # 2) a track OR a pi-zero in the signal cone has pT > 5 GeV - Prediscriminants = requireLeadTrack, - - # status flag indicating dead/masked ECAL crystals - minStatus = cms.uint32(12), - - # region around dead/masked ECAL crystals that is to be cut - dR = cms.double(0.08), - verbosity = cms.int32(0) -) diff --git a/RecoTauTag/RecoTau/src/AntiElectronDeadECAL.cc b/RecoTauTag/RecoTau/src/AntiElectronDeadECAL.cc new file mode 100644 index 0000000000000..080c333fcf01b --- /dev/null +++ b/RecoTauTag/RecoTau/src/AntiElectronDeadECAL.cc @@ -0,0 +1,159 @@ +#include "RecoTauTag/RecoTau/interface/AntiElectronDeadECAL.h" + +#include "FWCore/Framework/interface/ESHandle.h" +#include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h" +#include "CondFormats/EcalObjects/interface/EcalChannelStatus.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h" +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "DataFormats/TauReco/interface/PFTau.h" +#include "DataFormats/PatCandidates/interface/Tau.h" +#include "DataFormats/DetId/interface/DetId.h" +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "DataFormats/EcalDetId/interface/EEDetId.h" +#include "DataFormats/Math/interface/deltaR.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +AntiElectronDeadECAL::AntiElectronDeadECAL(const edm::ParameterSet& cfg) + : minStatus_(cfg.getParameter("minStatus")), + dR2_(std::pow(cfg.getParameter("dR"), 2)), + extrapolateToECalEntrance_(cfg.getParameter("extrapolateToECalEntrance")), + verbosity_(cfg.getParameter("verbosity")), + isFirstEvent_(true) {} + +AntiElectronDeadECAL::~AntiElectronDeadECAL() {} + +void AntiElectronDeadECAL::beginEvent(const edm::EventSetup& es) { + updateBadTowers(es); + positionAtECalEntrance_.beginEvent(es); +} + +namespace { + template + void loopXtals(std::map& nBadCrystals, + std::map& maxStatus, + std::map& sumEta, + std::map& sumPhi, + const EcalChannelStatus* channelStatus, + const CaloGeometry* caloGeometry, + const EcalTrigTowerConstituentsMap* ttMap, + unsigned minStatus, + const uint16_t statusMask) { + // NOTE: modified version of SUSY CAF code + // UserCode/SusyCAF/plugins/SusyCAF_EcalDeadChannels.cc + for (int i = 0; i < Id::kSizeForDenseIndexing; ++i) { + Id id = Id::unhashIndex(i); + if (id == Id(0)) + continue; + EcalChannelStatusMap::const_iterator it = channelStatus->getMap().find(id.rawId()); + unsigned status = (it == channelStatus->end()) ? 0 : (it->getStatusCode() & statusMask); + if (status >= minStatus) { + const GlobalPoint& point = caloGeometry->getPosition(id); + uint32_t key = ttMap->towerOf(id); + maxStatus[key] = std::max(status, maxStatus[key]); + ++nBadCrystals[key]; + sumEta[key] += point.eta(); + sumPhi[key] += point.phi(); + } + } + } +} // namespace + +void AntiElectronDeadECAL::updateBadTowers(const edm::EventSetup& es) { + // NOTE: modified version of SUSY CAF code + // UserCode/SusyCAF/plugins/SusyCAF_EcalDeadChannels.cc + + if (!isFirstEvent_ && !channelStatusWatcher_.check(es) && !caloGeometryWatcher_.check(es) && + !idealGeometryWatcher_.check(es)) + return; + + edm::ESHandle channelStatus; + es.get().get(channelStatus); + + edm::ESHandle caloGeometry; + es.get().get(caloGeometry); + + edm::ESHandle ttMap; + es.get().get(ttMap); + + std::map nBadCrystals, maxStatus; + std::map sumEta, sumPhi; + + loopXtals(nBadCrystals, + maxStatus, + sumEta, + sumPhi, + channelStatus.product(), + caloGeometry.product(), + ttMap.product(), + minStatus_, + statusMask_); + loopXtals(nBadCrystals, + maxStatus, + sumEta, + sumPhi, + channelStatus.product(), + caloGeometry.product(), + ttMap.product(), + minStatus_, + statusMask_); + + badTowers_.clear(); + for (auto it : nBadCrystals) { + uint32_t key = it.first; + badTowers_.push_back(TowerInfo(key, it.second, maxStatus[key], sumEta[key] / it.second, sumPhi[key] / it.second)); + } + + isFirstEvent_ = false; +} + +bool AntiElectronDeadECAL::operator()(const reco::Candidate* tau) const { + bool isNearBadTower = false; + double tau_eta = tau->eta(); + double tau_phi = tau->phi(); + const reco::Candidate* leadChargedHadron = nullptr; + if (extrapolateToECalEntrance_) { + const reco::PFTau* pfTau = dynamic_cast(tau); + if (pfTau != nullptr) { + leadChargedHadron = pfTau->leadChargedHadrCand().isNonnull() ? pfTau->leadChargedHadrCand().get() : nullptr; + } else { + const pat::Tau* patTau = dynamic_cast(tau); + if (patTau != nullptr) { + leadChargedHadron = patTau->leadChargedHadrCand().isNonnull() ? patTau->leadChargedHadrCand().get() : nullptr; + } + } + } + if (leadChargedHadron != nullptr) { + bool success = false; + reco::Candidate::Point positionAtECalEntrance = positionAtECalEntrance_(leadChargedHadron, success); + if (success) { + tau_eta = positionAtECalEntrance.eta(); + tau_phi = positionAtECalEntrance.phi(); + } + } + if (verbosity_) { + edm::LogPrint("TauAgainstEleDeadECAL") << ":"; + edm::LogPrint("TauAgainstEleDeadECAL") << " #badTowers = " << badTowers_.size(); + if (leadChargedHadron != nullptr) { + edm::LogPrint("TauAgainstEleDeadECAL") + << " leadChargedHadron (" << leadChargedHadron->pdgId() << "): Pt = " << leadChargedHadron->pt() + << ", eta at ECAL (vtx) = " << tau_eta << " (" << leadChargedHadron->eta() << ")" + << ", phi at ECAL (vtx) = " << tau_phi << " (" << leadChargedHadron->phi() << ")"; + } else { + edm::LogPrint("TauAgainstEleDeadECAL") + << " tau: Pt = " << tau->pt() << ", eta at vtx = " << tau_eta << ", phi at vtx = " << tau_phi; + } + } + for (auto const& badTower : badTowers_) { + if (deltaR2(badTower.eta_, badTower.phi_, tau_eta, tau_phi) < dR2_) { + if (verbosity_) { + edm::LogPrint("TauAgainstEleDeadECAL") + << " matches badTower: eta = " << badTower.eta_ << ", phi = " << badTower.phi_; + } + isNearBadTower = true; + break; + } + } + return isNearBadTower; +} diff --git a/RecoTauTag/RecoTau/src/AntiElectronIDMVA6.cc b/RecoTauTag/RecoTau/src/AntiElectronIDMVA6.cc index 33857cff1d02d..547dc14e07090 100644 --- a/RecoTauTag/RecoTau/src/AntiElectronIDMVA6.cc +++ b/RecoTauTag/RecoTau/src/AntiElectronIDMVA6.cc @@ -11,30 +11,27 @@ #include "CondFormats/DataRecord/interface/GBRWrapperRcd.h" #include "FWCore/Framework/interface/ESHandle.h" -#include "MagneticField/Engine/interface/MagneticField.h" -#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" - #include #include #include AntiElectronIDMVA6::AntiElectronIDMVA6(const edm::ParameterSet& cfg) - : isInitialized_(false), - mva_NoEleMatch_woGwoGSF_BL_(nullptr), - mva_NoEleMatch_wGwoGSF_BL_(nullptr), - mva_woGwGSF_BL_(nullptr), - mva_wGwGSF_BL_(nullptr), - mva_NoEleMatch_woGwoGSF_EC_(nullptr), - mva_NoEleMatch_wGwoGSF_EC_(nullptr), - mva_woGwGSF_EC_(nullptr), - mva_wGwGSF_EC_(nullptr) -{ - loadMVAfromDB_ = cfg.exists("loadMVAfromDB") ? cfg.getParameter("loadMVAfromDB"): false; - if ( !loadMVAfromDB_ ) { - if(cfg.exists("inputFileName")){ + : isInitialized_(false), + mva_NoEleMatch_woGwoGSF_BL_(nullptr), + mva_NoEleMatch_wGwoGSF_BL_(nullptr), + mva_woGwGSF_BL_(nullptr), + mva_wGwGSF_BL_(nullptr), + mva_NoEleMatch_woGwoGSF_EC_(nullptr), + mva_NoEleMatch_wGwoGSF_EC_(nullptr), + mva_woGwGSF_EC_(nullptr), + mva_wGwGSF_EC_(nullptr) { + loadMVAfromDB_ = cfg.exists("loadMVAfromDB") ? cfg.getParameter("loadMVAfromDB") : false; + if (!loadMVAfromDB_) { + if (cfg.exists("inputFileName")) { inputFileName_ = cfg.getParameter("inputFileName"); - }else throw cms::Exception("MVA input not defined") << "Requested to load tau MVA input from ROOT file but no file provided in cfg file"; - + } else + throw cms::Exception("MVA input not defined") + << "Requested to load tau MVA input from ROOT file but no file provided in cfg file"; } mvaName_NoEleMatch_woGwoGSF_BL_ = cfg.getParameter("mvaName_NoEleMatch_woGwoGSF_BL"); @@ -57,22 +54,20 @@ AntiElectronIDMVA6::AntiElectronIDMVA6(const edm::ParameterSet& cfg) Var_woGwGSF_Endcap_ = new Float_t[23]; Var_wGwGSF_Endcap_ = new Float_t[31]; - bField_ = 0; verbosity_ = 0; } -AntiElectronIDMVA6::~AntiElectronIDMVA6() -{ - delete [] Var_NoEleMatch_woGwoGSF_Barrel_; - delete [] Var_NoEleMatch_wGwoGSF_Barrel_; - delete [] Var_woGwGSF_Barrel_; - delete [] Var_wGwGSF_Barrel_; - delete [] Var_NoEleMatch_woGwoGSF_Endcap_; - delete [] Var_NoEleMatch_wGwoGSF_Endcap_; - delete [] Var_woGwGSF_Endcap_; - delete [] Var_wGwGSF_Endcap_; - - if ( !loadMVAfromDB_ ){ +AntiElectronIDMVA6::~AntiElectronIDMVA6() { + delete[] Var_NoEleMatch_woGwoGSF_Barrel_; + delete[] Var_NoEleMatch_wGwoGSF_Barrel_; + delete[] Var_woGwGSF_Barrel_; + delete[] Var_wGwGSF_Barrel_; + delete[] Var_NoEleMatch_woGwoGSF_Endcap_; + delete[] Var_NoEleMatch_wGwoGSF_Endcap_; + delete[] Var_woGwGSF_Endcap_; + delete[] Var_wGwGSF_Endcap_; + + if (!loadMVAfromDB_) { delete mva_NoEleMatch_woGwoGSF_BL_; delete mva_NoEleMatch_wGwoGSF_BL_; delete mva_woGwGSF_BL_; @@ -83,65 +78,59 @@ AntiElectronIDMVA6::~AntiElectronIDMVA6() delete mva_wGwGSF_EC_; } - for ( std::vector::iterator it = inputFilesToDelete_.begin(); - it != inputFilesToDelete_.end(); ++it ) { + for (std::vector::iterator it = inputFilesToDelete_.begin(); it != inputFilesToDelete_.end(); ++it) { delete (*it); } } -namespace -{ - const GBRForest* loadMVAfromFile(TFile* inputFile, const std::string& mvaName) - { +namespace { + const GBRForest* loadMVAfromFile(TFile* inputFile, const std::string& mvaName) { const GBRForest* mva = (GBRForest*)inputFile->Get(mvaName.data()); - if ( !mva ) + if (!mva) throw cms::Exception("PFRecoTauDiscriminationAgainstElectronMVA6::loadMVA") - << " Failed to load MVA = " << mvaName.data() << " from file " << " !!\n"; + << " Failed to load MVA = " << mvaName.data() << " from file " + << " !!\n"; return mva; } - const GBRForest* loadMVAfromDB(const edm::EventSetup& es, const std::string& mvaName) - { + const GBRForest* loadMVAfromDB(const edm::EventSetup& es, const std::string& mvaName) { edm::ESHandle mva; es.get().get(mvaName, mva); return mva.product(); } -} +} // namespace -void AntiElectronIDMVA6::beginEvent(const edm::Event& evt, const edm::EventSetup& es) -{ - if ( !isInitialized_ ) { - if ( loadMVAfromDB_ ) { +void AntiElectronIDMVA6::beginEvent(const edm::Event& evt, const edm::EventSetup& es) { + if (!isInitialized_) { + if (loadMVAfromDB_) { mva_NoEleMatch_woGwoGSF_BL_ = loadMVAfromDB(es, mvaName_NoEleMatch_woGwoGSF_BL_); - mva_NoEleMatch_wGwoGSF_BL_ = loadMVAfromDB(es, mvaName_NoEleMatch_wGwoGSF_BL_); - mva_woGwGSF_BL_ = loadMVAfromDB(es, mvaName_woGwGSF_BL_); - mva_wGwGSF_BL_ = loadMVAfromDB(es, mvaName_wGwGSF_BL_); + mva_NoEleMatch_wGwoGSF_BL_ = loadMVAfromDB(es, mvaName_NoEleMatch_wGwoGSF_BL_); + mva_woGwGSF_BL_ = loadMVAfromDB(es, mvaName_woGwGSF_BL_); + mva_wGwGSF_BL_ = loadMVAfromDB(es, mvaName_wGwGSF_BL_); mva_NoEleMatch_woGwoGSF_EC_ = loadMVAfromDB(es, mvaName_NoEleMatch_woGwoGSF_EC_); - mva_NoEleMatch_wGwoGSF_EC_ = loadMVAfromDB(es, mvaName_NoEleMatch_wGwoGSF_EC_); - mva_woGwGSF_EC_ = loadMVAfromDB(es, mvaName_woGwGSF_EC_); - mva_wGwGSF_EC_ = loadMVAfromDB(es, mvaName_wGwGSF_EC_); + mva_NoEleMatch_wGwoGSF_EC_ = loadMVAfromDB(es, mvaName_NoEleMatch_wGwoGSF_EC_); + mva_woGwGSF_EC_ = loadMVAfromDB(es, mvaName_woGwGSF_EC_); + mva_wGwGSF_EC_ = loadMVAfromDB(es, mvaName_wGwGSF_EC_); } else { - if ( inputFileName_.location() == edm::FileInPath::Unknown ) throw cms::Exception("PFRecoTauDiscriminationAgainstElectronMVA6::loadMVA") - << " Failed to find File = " << inputFileName_ << " !!\n"; - TFile* inputFile = new TFile(inputFileName_.fullPath().data()); + if (inputFileName_.location() == edm::FileInPath::Unknown) + throw cms::Exception("PFRecoTauDiscriminationAgainstElectronMVA6::loadMVA") + << " Failed to find File = " << inputFileName_ << " !!\n"; + TFile* inputFile = new TFile(inputFileName_.fullPath().data()); mva_NoEleMatch_woGwoGSF_BL_ = loadMVAfromFile(inputFile, mvaName_NoEleMatch_woGwoGSF_BL_); - mva_NoEleMatch_wGwoGSF_BL_ = loadMVAfromFile(inputFile, mvaName_NoEleMatch_wGwoGSF_BL_); - mva_woGwGSF_BL_ = loadMVAfromFile(inputFile, mvaName_woGwGSF_BL_); - mva_wGwGSF_BL_ = loadMVAfromFile(inputFile, mvaName_wGwGSF_BL_); + mva_NoEleMatch_wGwoGSF_BL_ = loadMVAfromFile(inputFile, mvaName_NoEleMatch_wGwoGSF_BL_); + mva_woGwGSF_BL_ = loadMVAfromFile(inputFile, mvaName_woGwGSF_BL_); + mva_wGwGSF_BL_ = loadMVAfromFile(inputFile, mvaName_wGwGSF_BL_); mva_NoEleMatch_woGwoGSF_EC_ = loadMVAfromFile(inputFile, mvaName_NoEleMatch_woGwoGSF_EC_); - mva_NoEleMatch_wGwoGSF_EC_ = loadMVAfromFile(inputFile, mvaName_NoEleMatch_wGwoGSF_EC_); - mva_woGwGSF_EC_ = loadMVAfromFile(inputFile, mvaName_woGwGSF_EC_); - mva_wGwGSF_EC_ = loadMVAfromFile(inputFile, mvaName_wGwGSF_EC_); - inputFilesToDelete_.push_back(inputFile); + mva_NoEleMatch_wGwoGSF_EC_ = loadMVAfromFile(inputFile, mvaName_NoEleMatch_wGwoGSF_EC_); + mva_woGwGSF_EC_ = loadMVAfromFile(inputFile, mvaName_woGwGSF_EC_); + mva_wGwGSF_EC_ = loadMVAfromFile(inputFile, mvaName_wGwGSF_EC_); + inputFilesToDelete_.push_back(inputFile); } isInitialized_ = true; } - - edm::ESHandle pSetup; - es.get().get(pSetup); - bField_ = pSetup->inTesla(GlobalPoint(0,0,0)).z(); + positionAtECalEntrance_.beginEvent(es); } double AntiElectronIDMVA6::MVAValue(Float_t TauPt, @@ -180,57 +169,60 @@ double AntiElectronIDMVA6::MVAValue(Float_t TauPt, Float_t ElecDeltaPhi, Float_t ElecMvaInSigmaEtaEta, Float_t ElecMvaInHadEnergy, - Float_t ElecMvaInDeltaEta) -{ - double sumPt = 0.; - double dEta2 = 0.; - double dPhi2 = 0.; + Float_t ElecMvaInDeltaEta) { + double sumPt = 0.; + double dEta2 = 0.; + double dPhi2 = 0.; double sumPt2 = 0.; - for ( unsigned int i = 0 ; i < GammasPtInSigCone.size() ; ++i ) { - double pt_i = GammasPtInSigCone[i]; + for (unsigned int i = 0; i < GammasPtInSigCone.size(); ++i) { + double pt_i = GammasPtInSigCone[i]; double phi_i = GammasdPhiInSigCone[i]; - if ( GammasdPhiInSigCone[i] > M_PI ) phi_i = GammasdPhiInSigCone[i] - 2*M_PI; - else if ( GammasdPhiInSigCone[i] < -M_PI ) phi_i = GammasdPhiInSigCone[i] + 2*M_PI; + if (GammasdPhiInSigCone[i] > M_PI) + phi_i = GammasdPhiInSigCone[i] - 2 * M_PI; + else if (GammasdPhiInSigCone[i] < -M_PI) + phi_i = GammasdPhiInSigCone[i] + 2 * M_PI; double eta_i = GammasdEtaInSigCone[i]; - sumPt += pt_i; - sumPt2 += (pt_i*pt_i); - dEta2 += (pt_i*eta_i*eta_i); - dPhi2 += (pt_i*phi_i*phi_i); + sumPt += pt_i; + sumPt2 += (pt_i * pt_i); + dEta2 += (pt_i * eta_i * eta_i); + dPhi2 += (pt_i * phi_i * phi_i); } Float_t TauGammaEnFracIn = -99.; - if ( TauPt > 0. ) { - TauGammaEnFracIn = sumPt/TauPt; + if (TauPt > 0.) { + TauGammaEnFracIn = sumPt / TauPt; } - if ( sumPt > 0. ) { + if (sumPt > 0.) { dEta2 /= sumPt; dPhi2 /= sumPt; } - Float_t TauGammaEtaMomIn = std::sqrt(dEta2)*std::sqrt(TauGammaEnFracIn)*TauPt; - Float_t TauGammaPhiMomIn = std::sqrt(dPhi2)*std::sqrt(TauGammaEnFracIn)*TauPt; + Float_t TauGammaEtaMomIn = std::sqrt(dEta2) * std::sqrt(TauGammaEnFracIn) * TauPt; + Float_t TauGammaPhiMomIn = std::sqrt(dPhi2) * std::sqrt(TauGammaEnFracIn) * TauPt; - sumPt = 0.; - dEta2 = 0.; - dPhi2 = 0.; + sumPt = 0.; + dEta2 = 0.; + dPhi2 = 0.; sumPt2 = 0.; - for ( unsigned int i = 0 ; i < GammasPtOutSigCone.size() ; ++i ) { - double pt_i = GammasPtOutSigCone[i]; + for (unsigned int i = 0; i < GammasPtOutSigCone.size(); ++i) { + double pt_i = GammasPtOutSigCone[i]; double phi_i = GammasdPhiOutSigCone[i]; - if ( GammasdPhiOutSigCone[i] > M_PI ) phi_i = GammasdPhiOutSigCone[i] - 2*M_PI; - else if ( GammasdPhiOutSigCone[i] < -M_PI ) phi_i = GammasdPhiOutSigCone[i] + 2*M_PI; + if (GammasdPhiOutSigCone[i] > M_PI) + phi_i = GammasdPhiOutSigCone[i] - 2 * M_PI; + else if (GammasdPhiOutSigCone[i] < -M_PI) + phi_i = GammasdPhiOutSigCone[i] + 2 * M_PI; double eta_i = GammasdEtaOutSigCone[i]; - sumPt += pt_i; - sumPt2 += (pt_i*pt_i); - dEta2 += (pt_i*eta_i*eta_i); - dPhi2 += (pt_i*phi_i*phi_i); + sumPt += pt_i; + sumPt2 += (pt_i * pt_i); + dEta2 += (pt_i * eta_i * eta_i); + dPhi2 += (pt_i * phi_i * phi_i); } - Float_t TauGammaEnFracOut = sumPt/TauPt; - if ( sumPt > 0. ) { + Float_t TauGammaEnFracOut = sumPt / TauPt; + if (sumPt > 0.) { dEta2 /= sumPt; dPhi2 /= sumPt; } - Float_t TauGammaEtaMomOut = std::sqrt(dEta2)*std::sqrt(TauGammaEnFracOut)*TauPt; - Float_t TauGammaPhiMomOut = std::sqrt(dPhi2)*std::sqrt(TauGammaEnFracOut)*TauPt; - + Float_t TauGammaEtaMomOut = std::sqrt(dEta2) * std::sqrt(TauGammaEnFracOut) * TauPt; + Float_t TauGammaPhiMomOut = std::sqrt(dPhi2) * std::sqrt(TauGammaEnFracOut) * TauPt; + return MVAValue(TauPt, TauEtaAtEcalEntrance, TauPhi, @@ -306,27 +298,25 @@ double AntiElectronIDMVA6::MVAValue(Float_t TauPt, Float_t ElecDeltaPhi, Float_t ElecMvaInSigmaEtaEta, Float_t ElecMvaInHadEnergy, - Float_t ElecMvaInDeltaEta) -{ - - if ( !isInitialized_ ) { - throw cms::Exception("ClassNotInitialized") - << " AntiElectronMVA not properly initialized !!\n"; + Float_t ElecMvaInDeltaEta) { + if (!isInitialized_) { + throw cms::Exception("ClassNotInitialized") << " AntiElectronMVA not properly initialized !!\n"; } double mvaValue = -99.; - + const float ECALBarrelEndcapEtaBorder = 1.479; - float ElecDeltaPinPoutOverPin = (ElecPin > 0.0) ? (std::abs(ElecPin - ElecPout)/ElecPin) : 1.0; - float ElecEecalOverPout = (ElecPout > 0.0) ? (ElecEecal/ElecPout) : 20.0; - float ElecNumHitsDiffOverSum = ((ElecGSFNumHits + ElecKFNumHits) > 0.0) ? - ((ElecGSFNumHits - ElecKFNumHits)/(ElecGSFNumHits + ElecKFNumHits)) : 1.0; - - if ( deltaR(TauEtaAtEcalEntrance, TauPhi, ElecEta, ElecPhi) > 0.3 && TauSignalPFGammaCandsIn == 0 && TauHasGsf < 0.5) { - if ( std::abs(TauEtaAtEcalEntrance) < ECALBarrelEndcapEtaBorder ){ + float ElecDeltaPinPoutOverPin = (ElecPin > 0.0) ? (std::abs(ElecPin - ElecPout) / ElecPin) : 1.0; + float ElecEecalOverPout = (ElecPout > 0.0) ? (ElecEecal / ElecPout) : 20.0; + float ElecNumHitsDiffOverSum = ((ElecGSFNumHits + ElecKFNumHits) > 0.0) + ? ((ElecGSFNumHits - ElecKFNumHits) / (ElecGSFNumHits + ElecKFNumHits)) + : 1.0; + + if (deltaR(TauEtaAtEcalEntrance, TauPhi, ElecEta, ElecPhi) > 0.3 && TauSignalPFGammaCandsIn == 0 && TauHasGsf < 0.5) { + if (std::abs(TauEtaAtEcalEntrance) < ECALBarrelEndcapEtaBorder) { Var_NoEleMatch_woGwoGSF_Barrel_[0] = TauEtaAtEcalEntrance; Var_NoEleMatch_woGwoGSF_Barrel_[1] = TauLeadChargedPFCandEtaAtEcalEntrance; - Var_NoEleMatch_woGwoGSF_Barrel_[2] = std::min(float(2.), TauLeadChargedPFCandPt/std::max(float(1.), TauPt)); + Var_NoEleMatch_woGwoGSF_Barrel_[2] = std::min(float(2.), TauLeadChargedPFCandPt / std::max(float(1.), TauPt)); Var_NoEleMatch_woGwoGSF_Barrel_[3] = std::log(std::max(float(1.), TauPt)); Var_NoEleMatch_woGwoGSF_Barrel_[4] = TauEmFraction; Var_NoEleMatch_woGwoGSF_Barrel_[5] = TauLeadPFChargedHadrHoP; @@ -338,7 +328,7 @@ double AntiElectronIDMVA6::MVAValue(Float_t TauPt, } else { Var_NoEleMatch_woGwoGSF_Endcap_[0] = TauEtaAtEcalEntrance; Var_NoEleMatch_woGwoGSF_Endcap_[1] = TauLeadChargedPFCandEtaAtEcalEntrance; - Var_NoEleMatch_woGwoGSF_Endcap_[2] = std::min(float(2.), TauLeadChargedPFCandPt/std::max(float(1.), TauPt)); + Var_NoEleMatch_woGwoGSF_Endcap_[2] = std::min(float(2.), TauLeadChargedPFCandPt / std::max(float(1.), TauPt)); Var_NoEleMatch_woGwoGSF_Endcap_[3] = std::log(std::max(float(1.), TauPt)); Var_NoEleMatch_woGwoGSF_Endcap_[4] = TauEmFraction; Var_NoEleMatch_woGwoGSF_Endcap_[5] = TauLeadPFChargedHadrHoP; @@ -347,19 +337,19 @@ double AntiElectronIDMVA6::MVAValue(Float_t TauPt, Var_NoEleMatch_woGwoGSF_Endcap_[8] = TaudCrackEta; mvaValue = mva_NoEleMatch_woGwoGSF_EC_->GetClassifier(Var_NoEleMatch_woGwoGSF_Endcap_); } - } - else if ( deltaR(TauEtaAtEcalEntrance, TauPhi, ElecEta, ElecPhi) > 0.3 && TauSignalPFGammaCandsIn > 0 && TauHasGsf < 0.5 ) { - if ( std::abs(TauEtaAtEcalEntrance) < ECALBarrelEndcapEtaBorder ){ - Var_NoEleMatch_wGwoGSF_Barrel_[0] = TauEtaAtEcalEntrance; - Var_NoEleMatch_wGwoGSF_Barrel_[1] = TauLeadChargedPFCandEtaAtEcalEntrance; - Var_NoEleMatch_wGwoGSF_Barrel_[2] = std::min(float(2.), TauLeadChargedPFCandPt/std::max(float(1.), TauPt)); - Var_NoEleMatch_wGwoGSF_Barrel_[3] = std::log(std::max(float(1.), TauPt)); - Var_NoEleMatch_wGwoGSF_Barrel_[4] = TauEmFraction; - Var_NoEleMatch_wGwoGSF_Barrel_[5] = TauSignalPFGammaCandsIn; - Var_NoEleMatch_wGwoGSF_Barrel_[6] = TauSignalPFGammaCandsOut; - Var_NoEleMatch_wGwoGSF_Barrel_[7] = TauLeadPFChargedHadrHoP; - Var_NoEleMatch_wGwoGSF_Barrel_[8] = TauLeadPFChargedHadrEoP; - Var_NoEleMatch_wGwoGSF_Barrel_[9] = TauVisMassIn; + } else if (deltaR(TauEtaAtEcalEntrance, TauPhi, ElecEta, ElecPhi) > 0.3 && TauSignalPFGammaCandsIn > 0 && + TauHasGsf < 0.5) { + if (std::abs(TauEtaAtEcalEntrance) < ECALBarrelEndcapEtaBorder) { + Var_NoEleMatch_wGwoGSF_Barrel_[0] = TauEtaAtEcalEntrance; + Var_NoEleMatch_wGwoGSF_Barrel_[1] = TauLeadChargedPFCandEtaAtEcalEntrance; + Var_NoEleMatch_wGwoGSF_Barrel_[2] = std::min(float(2.), TauLeadChargedPFCandPt / std::max(float(1.), TauPt)); + Var_NoEleMatch_wGwoGSF_Barrel_[3] = std::log(std::max(float(1.), TauPt)); + Var_NoEleMatch_wGwoGSF_Barrel_[4] = TauEmFraction; + Var_NoEleMatch_wGwoGSF_Barrel_[5] = TauSignalPFGammaCandsIn; + Var_NoEleMatch_wGwoGSF_Barrel_[6] = TauSignalPFGammaCandsOut; + Var_NoEleMatch_wGwoGSF_Barrel_[7] = TauLeadPFChargedHadrHoP; + Var_NoEleMatch_wGwoGSF_Barrel_[8] = TauLeadPFChargedHadrEoP; + Var_NoEleMatch_wGwoGSF_Barrel_[9] = TauVisMassIn; Var_NoEleMatch_wGwoGSF_Barrel_[10] = TauGammaEtaMomIn; Var_NoEleMatch_wGwoGSF_Barrel_[11] = TauGammaEtaMomOut; Var_NoEleMatch_wGwoGSF_Barrel_[12] = TauGammaPhiMomIn; @@ -370,16 +360,16 @@ double AntiElectronIDMVA6::MVAValue(Float_t TauPt, Var_NoEleMatch_wGwoGSF_Barrel_[17] = TaudCrackPhi; mvaValue = mva_NoEleMatch_wGwoGSF_BL_->GetClassifier(Var_NoEleMatch_wGwoGSF_Barrel_); } else { - Var_NoEleMatch_wGwoGSF_Endcap_[0] = TauEtaAtEcalEntrance; - Var_NoEleMatch_wGwoGSF_Endcap_[1] = TauLeadChargedPFCandEtaAtEcalEntrance; - Var_NoEleMatch_wGwoGSF_Endcap_[2] = std::min(float(2.), TauLeadChargedPFCandPt/std::max(float(1.), TauPt)); - Var_NoEleMatch_wGwoGSF_Endcap_[3] = std::log(std::max(float(1.), TauPt)); - Var_NoEleMatch_wGwoGSF_Endcap_[4] = TauEmFraction; - Var_NoEleMatch_wGwoGSF_Endcap_[5] = TauSignalPFGammaCandsIn; - Var_NoEleMatch_wGwoGSF_Endcap_[6] = TauSignalPFGammaCandsOut; - Var_NoEleMatch_wGwoGSF_Endcap_[7] = TauLeadPFChargedHadrHoP; - Var_NoEleMatch_wGwoGSF_Endcap_[8] = TauLeadPFChargedHadrEoP; - Var_NoEleMatch_wGwoGSF_Endcap_[9] = TauVisMassIn; + Var_NoEleMatch_wGwoGSF_Endcap_[0] = TauEtaAtEcalEntrance; + Var_NoEleMatch_wGwoGSF_Endcap_[1] = TauLeadChargedPFCandEtaAtEcalEntrance; + Var_NoEleMatch_wGwoGSF_Endcap_[2] = std::min(float(2.), TauLeadChargedPFCandPt / std::max(float(1.), TauPt)); + Var_NoEleMatch_wGwoGSF_Endcap_[3] = std::log(std::max(float(1.), TauPt)); + Var_NoEleMatch_wGwoGSF_Endcap_[4] = TauEmFraction; + Var_NoEleMatch_wGwoGSF_Endcap_[5] = TauSignalPFGammaCandsIn; + Var_NoEleMatch_wGwoGSF_Endcap_[6] = TauSignalPFGammaCandsOut; + Var_NoEleMatch_wGwoGSF_Endcap_[7] = TauLeadPFChargedHadrHoP; + Var_NoEleMatch_wGwoGSF_Endcap_[8] = TauLeadPFChargedHadrEoP; + Var_NoEleMatch_wGwoGSF_Endcap_[9] = TauVisMassIn; Var_NoEleMatch_wGwoGSF_Endcap_[10] = TauGammaEtaMomIn; Var_NoEleMatch_wGwoGSF_Endcap_[11] = TauGammaEtaMomOut; Var_NoEleMatch_wGwoGSF_Endcap_[12] = TauGammaPhiMomIn; @@ -389,26 +379,25 @@ double AntiElectronIDMVA6::MVAValue(Float_t TauPt, Var_NoEleMatch_wGwoGSF_Endcap_[16] = TaudCrackEta; mvaValue = mva_NoEleMatch_wGwoGSF_EC_->GetClassifier(Var_NoEleMatch_wGwoGSF_Endcap_); } - } - else if ( TauSignalPFGammaCandsIn == 0 && TauHasGsf > 0.5 ) { - if ( std::abs(TauEtaAtEcalEntrance) < ECALBarrelEndcapEtaBorder ) { - Var_woGwGSF_Barrel_[0] = std::max(float(-0.1), ElecEtotOverPin); - Var_woGwGSF_Barrel_[1] = std::log(std::max(float(0.01), ElecChi2NormGSF)); - Var_woGwGSF_Barrel_[2] = ElecGSFNumHits; - Var_woGwGSF_Barrel_[3] = std::log(std::max(float(0.01), ElecGSFTrackResol)); - Var_woGwGSF_Barrel_[4] = ElecGSFTracklnPt; - Var_woGwGSF_Barrel_[5] = ElecNumHitsDiffOverSum; - Var_woGwGSF_Barrel_[6] = std::log(std::max(float(0.01), ElecChi2NormKF)); - Var_woGwGSF_Barrel_[7] = std::min(ElecDeltaPinPoutOverPin, float(1.)); - Var_woGwGSF_Barrel_[8] = std::min(ElecEecalOverPout, float(20.)); - Var_woGwGSF_Barrel_[9] = ElecDeltaEta; + } else if (TauSignalPFGammaCandsIn == 0 && TauHasGsf > 0.5) { + if (std::abs(TauEtaAtEcalEntrance) < ECALBarrelEndcapEtaBorder) { + Var_woGwGSF_Barrel_[0] = std::max(float(-0.1), ElecEtotOverPin); + Var_woGwGSF_Barrel_[1] = std::log(std::max(float(0.01), ElecChi2NormGSF)); + Var_woGwGSF_Barrel_[2] = ElecGSFNumHits; + Var_woGwGSF_Barrel_[3] = std::log(std::max(float(0.01), ElecGSFTrackResol)); + Var_woGwGSF_Barrel_[4] = ElecGSFTracklnPt; + Var_woGwGSF_Barrel_[5] = ElecNumHitsDiffOverSum; + Var_woGwGSF_Barrel_[6] = std::log(std::max(float(0.01), ElecChi2NormKF)); + Var_woGwGSF_Barrel_[7] = std::min(ElecDeltaPinPoutOverPin, float(1.)); + Var_woGwGSF_Barrel_[8] = std::min(ElecEecalOverPout, float(20.)); + Var_woGwGSF_Barrel_[9] = ElecDeltaEta; Var_woGwGSF_Barrel_[10] = ElecDeltaPhi; Var_woGwGSF_Barrel_[11] = std::min(ElecMvaInSigmaEtaEta, float(0.01)); Var_woGwGSF_Barrel_[12] = std::min(ElecMvaInHadEnergy, float(20.)); Var_woGwGSF_Barrel_[13] = std::min(ElecMvaInDeltaEta, float(0.1)); Var_woGwGSF_Barrel_[14] = TauEtaAtEcalEntrance; Var_woGwGSF_Barrel_[15] = TauLeadChargedPFCandEtaAtEcalEntrance; - Var_woGwGSF_Barrel_[16] = std::min(float(2.), TauLeadChargedPFCandPt/std::max(float(1.), TauPt)); + Var_woGwGSF_Barrel_[16] = std::min(float(2.), TauLeadChargedPFCandPt / std::max(float(1.), TauPt)); Var_woGwGSF_Barrel_[17] = std::log(std::max(float(1.), TauPt)); Var_woGwGSF_Barrel_[18] = TauEmFraction; Var_woGwGSF_Barrel_[19] = TauLeadPFChargedHadrHoP; @@ -418,23 +407,23 @@ double AntiElectronIDMVA6::MVAValue(Float_t TauPt, Var_woGwGSF_Barrel_[23] = TaudCrackPhi; mvaValue = mva_woGwGSF_BL_->GetClassifier(Var_woGwGSF_Barrel_); } else { - Var_woGwGSF_Endcap_[0] = std::max(float(-0.1), ElecEtotOverPin); - Var_woGwGSF_Endcap_[1] = std::log(std::max(float(0.01), ElecChi2NormGSF)); - Var_woGwGSF_Endcap_[2] = ElecGSFNumHits; - Var_woGwGSF_Endcap_[3] = std::log(std::max(float(0.01), ElecGSFTrackResol)); - Var_woGwGSF_Endcap_[4] = ElecGSFTracklnPt; - Var_woGwGSF_Endcap_[5] = ElecNumHitsDiffOverSum; - Var_woGwGSF_Endcap_[6] = std::log(std::max(float(0.01), ElecChi2NormKF)); - Var_woGwGSF_Endcap_[7] = std::min(ElecDeltaPinPoutOverPin, float(1.)); - Var_woGwGSF_Endcap_[8] = std::min(ElecEecalOverPout, float(20.)); - Var_woGwGSF_Endcap_[9] = ElecDeltaEta; + Var_woGwGSF_Endcap_[0] = std::max(float(-0.1), ElecEtotOverPin); + Var_woGwGSF_Endcap_[1] = std::log(std::max(float(0.01), ElecChi2NormGSF)); + Var_woGwGSF_Endcap_[2] = ElecGSFNumHits; + Var_woGwGSF_Endcap_[3] = std::log(std::max(float(0.01), ElecGSFTrackResol)); + Var_woGwGSF_Endcap_[4] = ElecGSFTracklnPt; + Var_woGwGSF_Endcap_[5] = ElecNumHitsDiffOverSum; + Var_woGwGSF_Endcap_[6] = std::log(std::max(float(0.01), ElecChi2NormKF)); + Var_woGwGSF_Endcap_[7] = std::min(ElecDeltaPinPoutOverPin, float(1.)); + Var_woGwGSF_Endcap_[8] = std::min(ElecEecalOverPout, float(20.)); + Var_woGwGSF_Endcap_[9] = ElecDeltaEta; Var_woGwGSF_Endcap_[10] = ElecDeltaPhi; Var_woGwGSF_Endcap_[11] = std::min(ElecMvaInSigmaEtaEta, float(0.01)); Var_woGwGSF_Endcap_[12] = std::min(ElecMvaInHadEnergy, float(20.)); Var_woGwGSF_Endcap_[13] = std::min(ElecMvaInDeltaEta, float(0.1)); Var_woGwGSF_Endcap_[14] = TauEtaAtEcalEntrance; Var_woGwGSF_Endcap_[15] = TauLeadChargedPFCandEtaAtEcalEntrance; - Var_woGwGSF_Endcap_[16] = std::min(float(2.), TauLeadChargedPFCandPt/std::max(float(1.), TauPt)); + Var_woGwGSF_Endcap_[16] = std::min(float(2.), TauLeadChargedPFCandPt / std::max(float(1.), TauPt)); Var_woGwGSF_Endcap_[17] = std::log(std::max(float(1.), TauPt)); Var_woGwGSF_Endcap_[18] = TauEmFraction; Var_woGwGSF_Endcap_[19] = TauLeadPFChargedHadrHoP; @@ -442,27 +431,26 @@ double AntiElectronIDMVA6::MVAValue(Float_t TauPt, Var_woGwGSF_Endcap_[21] = TauVisMassIn; Var_woGwGSF_Endcap_[22] = TaudCrackEta; mvaValue = mva_woGwGSF_EC_->GetClassifier(Var_woGwGSF_Endcap_); - } - } - else if ( TauSignalPFGammaCandsIn > 0 && TauHasGsf > 0.5 ) { - if ( std::abs(TauEtaAtEcalEntrance) < ECALBarrelEndcapEtaBorder ) { - Var_wGwGSF_Barrel_[0] = std::max(float(-0.1), ElecEtotOverPin); - Var_wGwGSF_Barrel_[1] = std::log(std::max(float(0.01), ElecChi2NormGSF)); - Var_wGwGSF_Barrel_[2] = ElecGSFNumHits; - Var_wGwGSF_Barrel_[3] = std::log(std::max(float(0.01), ElecGSFTrackResol)); - Var_wGwGSF_Barrel_[4] = ElecGSFTracklnPt; - Var_wGwGSF_Barrel_[5] = ElecNumHitsDiffOverSum; - Var_wGwGSF_Barrel_[6] = std::log(std::max(float(0.01), ElecChi2NormKF)); - Var_wGwGSF_Barrel_[7] = std::min(ElecDeltaPinPoutOverPin, float(1.)); - Var_wGwGSF_Barrel_[8] = std::min(ElecEecalOverPout, float(20.)); - Var_wGwGSF_Barrel_[9] = ElecDeltaEta; + } + } else if (TauSignalPFGammaCandsIn > 0 && TauHasGsf > 0.5) { + if (std::abs(TauEtaAtEcalEntrance) < ECALBarrelEndcapEtaBorder) { + Var_wGwGSF_Barrel_[0] = std::max(float(-0.1), ElecEtotOverPin); + Var_wGwGSF_Barrel_[1] = std::log(std::max(float(0.01), ElecChi2NormGSF)); + Var_wGwGSF_Barrel_[2] = ElecGSFNumHits; + Var_wGwGSF_Barrel_[3] = std::log(std::max(float(0.01), ElecGSFTrackResol)); + Var_wGwGSF_Barrel_[4] = ElecGSFTracklnPt; + Var_wGwGSF_Barrel_[5] = ElecNumHitsDiffOverSum; + Var_wGwGSF_Barrel_[6] = std::log(std::max(float(0.01), ElecChi2NormKF)); + Var_wGwGSF_Barrel_[7] = std::min(ElecDeltaPinPoutOverPin, float(1.)); + Var_wGwGSF_Barrel_[8] = std::min(ElecEecalOverPout, float(20.)); + Var_wGwGSF_Barrel_[9] = ElecDeltaEta; Var_wGwGSF_Barrel_[10] = ElecDeltaPhi; Var_wGwGSF_Barrel_[11] = std::min(ElecMvaInSigmaEtaEta, float(0.01)); Var_wGwGSF_Barrel_[12] = std::min(ElecMvaInHadEnergy, float(20.)); Var_wGwGSF_Barrel_[13] = std::min(ElecMvaInDeltaEta, float(0.1)); Var_wGwGSF_Barrel_[14] = TauEtaAtEcalEntrance; Var_wGwGSF_Barrel_[15] = TauLeadChargedPFCandEtaAtEcalEntrance; - Var_wGwGSF_Barrel_[16] = std::min(float(2.), TauLeadChargedPFCandPt/std::max(float(1.), TauPt)); + Var_wGwGSF_Barrel_[16] = std::min(float(2.), TauLeadChargedPFCandPt / std::max(float(1.), TauPt)); Var_wGwGSF_Barrel_[17] = std::log(std::max(float(1.), TauPt)); Var_wGwGSF_Barrel_[18] = TauEmFraction; Var_wGwGSF_Barrel_[19] = TauSignalPFGammaCandsIn; @@ -480,23 +468,23 @@ double AntiElectronIDMVA6::MVAValue(Float_t TauPt, Var_wGwGSF_Barrel_[31] = TaudCrackPhi; mvaValue = mva_wGwGSF_BL_->GetClassifier(Var_wGwGSF_Barrel_); } else { - Var_wGwGSF_Endcap_[0] = std::max(float(-0.1), ElecEtotOverPin); - Var_wGwGSF_Endcap_[1] = std::log(std::max(float(0.01), ElecChi2NormGSF)); - Var_wGwGSF_Endcap_[2] = ElecGSFNumHits; - Var_wGwGSF_Endcap_[3] = std::log(std::max(float(0.01), ElecGSFTrackResol)); - Var_wGwGSF_Endcap_[4] = ElecGSFTracklnPt; - Var_wGwGSF_Endcap_[5] = ElecNumHitsDiffOverSum; - Var_wGwGSF_Endcap_[6] = std::log(std::max(float(0.01), ElecChi2NormKF)); - Var_wGwGSF_Endcap_[7] = std::min(ElecDeltaPinPoutOverPin, float(1.)); - Var_wGwGSF_Endcap_[8] = std::min(ElecEecalOverPout, float(20.)); - Var_wGwGSF_Endcap_[9] = ElecDeltaEta; + Var_wGwGSF_Endcap_[0] = std::max(float(-0.1), ElecEtotOverPin); + Var_wGwGSF_Endcap_[1] = std::log(std::max(float(0.01), ElecChi2NormGSF)); + Var_wGwGSF_Endcap_[2] = ElecGSFNumHits; + Var_wGwGSF_Endcap_[3] = std::log(std::max(float(0.01), ElecGSFTrackResol)); + Var_wGwGSF_Endcap_[4] = ElecGSFTracklnPt; + Var_wGwGSF_Endcap_[5] = ElecNumHitsDiffOverSum; + Var_wGwGSF_Endcap_[6] = std::log(std::max(float(0.01), ElecChi2NormKF)); + Var_wGwGSF_Endcap_[7] = std::min(ElecDeltaPinPoutOverPin, float(1.)); + Var_wGwGSF_Endcap_[8] = std::min(ElecEecalOverPout, float(20.)); + Var_wGwGSF_Endcap_[9] = ElecDeltaEta; Var_wGwGSF_Endcap_[10] = ElecDeltaPhi; Var_wGwGSF_Endcap_[11] = std::min(ElecMvaInSigmaEtaEta, float(0.01)); Var_wGwGSF_Endcap_[12] = std::min(ElecMvaInHadEnergy, float(20.)); Var_wGwGSF_Endcap_[13] = std::min(ElecMvaInDeltaEta, float(0.1)); Var_wGwGSF_Endcap_[14] = TauEtaAtEcalEntrance; Var_wGwGSF_Endcap_[15] = TauLeadChargedPFCandEtaAtEcalEntrance; - Var_wGwGSF_Endcap_[16] = std::min(float(2.), TauLeadChargedPFCandPt/std::max(float(1.), TauPt)); + Var_wGwGSF_Endcap_[16] = std::min(float(2.), TauLeadChargedPFCandPt / std::max(float(1.), TauPt)); Var_wGwGSF_Endcap_[17] = std::log(std::max(float(1.), TauPt)); Var_wGwGSF_Endcap_[18] = TauEmFraction; Var_wGwGSF_Endcap_[19] = TauSignalPFGammaCandsIn; @@ -512,13 +500,12 @@ double AntiElectronIDMVA6::MVAValue(Float_t TauPt, Var_wGwGSF_Endcap_[29] = TauGammaEnFracOut; Var_wGwGSF_Endcap_[30] = TaudCrackEta; mvaValue = mva_wGwGSF_EC_->GetClassifier(Var_wGwGSF_Endcap_); - } + } } return mvaValue; } -double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau, - const reco::GsfElectron& theGsfEle) +double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau, const reco::GsfElectron& theGsfEle) { // === tau variables === @@ -526,27 +513,32 @@ double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau, float sumEtaTimesEnergy = 0.; float sumEnergy = 0.; const std::vector& signalPFCands = thePFTau.signalPFCands(); - for ( const auto & pfCandidate : signalPFCands ) { - sumEtaTimesEnergy += pfCandidate->positionAtECALEntrance().eta()*pfCandidate->energy(); + for (const auto& pfCandidate : signalPFCands) { + sumEtaTimesEnergy += pfCandidate->positionAtECALEntrance().eta() * pfCandidate->energy(); sumEnergy += pfCandidate->energy(); } - if ( sumEnergy > 0. ) { - TauEtaAtEcalEntrance = sumEtaTimesEnergy/sumEnergy; + if (sumEnergy > 0.) { + TauEtaAtEcalEntrance = sumEtaTimesEnergy / sumEnergy; } - + float TauLeadChargedPFCandEtaAtEcalEntrance = -99.; float TauLeadChargedPFCandPt = -99.; - for ( const auto & pfCandidate : signalPFCands ) { + for (const auto& pfCandidate : signalPFCands) { const reco::Track* track = nullptr; - if ( pfCandidate->trackRef().isNonnull() ) track = pfCandidate->trackRef().get(); - else if ( pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->innerTrack().isNonnull() ) track = pfCandidate->muonRef()->innerTrack().get(); - else if ( pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->globalTrack().isNonnull() ) track = pfCandidate->muonRef()->globalTrack().get(); - else if ( pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->outerTrack().isNonnull() ) track = pfCandidate->muonRef()->outerTrack().get(); - else if ( pfCandidate->gsfTrackRef().isNonnull() ) track = pfCandidate->gsfTrackRef().get(); - if ( track ) { - if ( track->pt() > TauLeadChargedPFCandPt ) { - TauLeadChargedPFCandEtaAtEcalEntrance = pfCandidate->positionAtECALEntrance().eta(); - TauLeadChargedPFCandPt = track->pt(); + if (pfCandidate->trackRef().isNonnull()) + track = pfCandidate->trackRef().get(); + else if (pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->innerTrack().isNonnull()) + track = pfCandidate->muonRef()->innerTrack().get(); + else if (pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->globalTrack().isNonnull()) + track = pfCandidate->muonRef()->globalTrack().get(); + else if (pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->outerTrack().isNonnull()) + track = pfCandidate->muonRef()->outerTrack().get(); + else if (pfCandidate->gsfTrackRef().isNonnull()) + track = pfCandidate->gsfTrackRef().get(); + if (track) { + if (track->pt() > TauLeadChargedPFCandPt) { + TauLeadChargedPFCandEtaAtEcalEntrance = pfCandidate->positionAtECALEntrance().eta(); + TauLeadChargedPFCandPt = track->pt(); } } } @@ -555,9 +547,9 @@ double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau, Float_t TauEmFraction = std::max(thePFTau.emFraction(), (Float_t)0.); Float_t TauLeadPFChargedHadrHoP = 0.; Float_t TauLeadPFChargedHadrEoP = 0.; - if ( thePFTau.leadPFChargedHadrCand()->p() > 0. ) { - TauLeadPFChargedHadrHoP = thePFTau.leadPFChargedHadrCand()->hcalEnergy()/thePFTau.leadPFChargedHadrCand()->p(); - TauLeadPFChargedHadrEoP = thePFTau.leadPFChargedHadrCand()->ecalEnergy()/thePFTau.leadPFChargedHadrCand()->p(); + if (thePFTau.leadPFChargedHadrCand()->p() > 0.) { + TauLeadPFChargedHadrHoP = thePFTau.leadPFChargedHadrCand()->hcalEnergy() / thePFTau.leadPFChargedHadrCand()->p(); + TauLeadPFChargedHadrEoP = thePFTau.leadPFChargedHadrCand()->ecalEnergy() / thePFTau.leadPFChargedHadrCand()->p(); } std::vector GammasdEtaInSigCone; @@ -566,20 +558,19 @@ double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau, std::vector GammasdEtaOutSigCone; std::vector GammasdPhiOutSigCone; std::vector GammasPtOutSigCone; - reco::Candidate::LorentzVector pfGammaSum(0,0,0,0); - reco::Candidate::LorentzVector pfChargedSum(0,0,0,0); - - for ( const auto & gamma : thePFTau.signalGammaCands() ) { + reco::Candidate::LorentzVector pfGammaSum(0, 0, 0, 0); + reco::Candidate::LorentzVector pfChargedSum(0, 0, 0, 0); + + for (const auto& gamma : thePFTau.signalGammaCands()) { float dR = deltaR(gamma->p4(), thePFTau.leadChargedHadrCand()->p4()); - float signalrad = std::max(0.05, std::min(0.10, 3.0/std::max(1.0, thePFTau.pt()))); + float signalrad = std::max(0.05, std::min(0.10, 3.0 / std::max(1.0, thePFTau.pt()))); // pfGammas inside the tau signal cone if (dR < signalrad) { - if ( thePFTau.leadChargedHadrCand().isNonnull() ) { + if (thePFTau.leadChargedHadrCand().isNonnull()) { GammasdEtaInSigCone.push_back(gamma->eta() - thePFTau.leadChargedHadrCand()->eta()); GammasdPhiInSigCone.push_back(gamma->phi() - thePFTau.leadChargedHadrCand()->phi()); - } - else { + } else { GammasdEtaInSigCone.push_back(gamma->eta() - thePFTau.eta()); GammasdPhiInSigCone.push_back(gamma->phi() - thePFTau.phi()); } @@ -588,28 +579,27 @@ double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau, } // pfGammas outside the tau signal cone else { - if ( thePFTau.leadChargedHadrCand().isNonnull() ) { + if (thePFTau.leadChargedHadrCand().isNonnull()) { GammasdEtaOutSigCone.push_back(gamma->eta() - thePFTau.leadChargedHadrCand()->eta()); GammasdPhiOutSigCone.push_back(gamma->phi() - thePFTau.leadChargedHadrCand()->phi()); - } - else { + } else { GammasdEtaOutSigCone.push_back(gamma->eta() - thePFTau.eta()); GammasdPhiOutSigCone.push_back(gamma->phi() - thePFTau.phi()); } GammasPtOutSigCone.push_back(gamma->pt()); } } - - for ( const auto & charged : thePFTau.signalChargedHadrCands() ) { + + for (const auto& charged : thePFTau.signalChargedHadrCands()) { float dR = deltaR(charged->p4(), thePFTau.leadChargedHadrCand()->p4()); - float signalrad = std::max(0.05, std::min(0.10, 3.0/std::max(1.0, thePFTau.pt()))); - + float signalrad = std::max(0.05, std::min(0.10, 3.0 / std::max(1.0, thePFTau.pt()))); + // charged particles inside the tau signal cone if (dR < signalrad) { - pfChargedSum += charged->p4(); + pfChargedSum += charged->p4(); } } - + Int_t TauSignalPFGammaCandsIn = GammasPtInSigCone.size(); Int_t TauSignalPFGammaCandsOut = GammasPtOutSigCone.size(); Float_t TauVisMassIn = (pfGammaSum + pfChargedSum).mass(); @@ -617,75 +607,80 @@ double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau, Float_t TauPhi = thePFTau.phi(); float sumPhiTimesEnergy = 0.; float sumEnergyPhi = 0.; - if ( !usePhiAtEcalEntranceExtrapolation_ ) { - for ( const auto & pfc : signalPFCands ) { - sumPhiTimesEnergy += pfc->positionAtECALEntrance().phi()*pfc->energy(); + if (!usePhiAtEcalEntranceExtrapolation_) { + for (const auto& pfc : signalPFCands) { + sumPhiTimesEnergy += pfc->positionAtECALEntrance().phi() * pfc->energy(); sumEnergyPhi += pfc->energy(); } - } - else{ - TauPhi= -99.; - for ( const auto & signalPFCand : signalPFCands ) { - reco::Candidate const* signalCand = signalPFCand.get(); + } else { + TauPhi = -99.; + for (const auto& signalPFCand : signalPFCands) { + reco::Candidate const* signalCand = signalPFCand.get(); float phi = thePFTau.phi(); - math::XYZPoint aPos; - if ( atECalEntrance(signalCand, aPos) ) phi = aPos.Phi(); - sumPhiTimesEnergy += phi*signalCand->energy(); + bool success = false; + reco::Candidate::Point aPos = positionAtECalEntrance_(signalCand, success); + if (success) { + phi = aPos.Phi(); + } + sumPhiTimesEnergy += phi * signalCand->energy(); sumEnergy += signalCand->energy(); } } - if ( sumEnergyPhi > 0. ) { - TauPhi = sumPhiTimesEnergy/sumEnergyPhi; + if (sumEnergyPhi > 0.) { + TauPhi = sumPhiTimesEnergy / sumEnergyPhi; } Float_t TaudCrackPhi = dCrackPhi(TauPhi, TauEtaAtEcalEntrance); Float_t TaudCrackEta = dCrackEta(TauEtaAtEcalEntrance); Float_t TauHasGsf = thePFTau.leadPFChargedHadrCand()->gsfTrackRef().isNonnull(); - + // === electron variables === Float_t ElecEta = theGsfEle.eta(); Float_t ElecPhi = theGsfEle.phi(); - + //Variables related to the electron Cluster Float_t ElecEe = 0.; Float_t ElecEgamma = 0.; reco::SuperClusterRef pfSuperCluster = theGsfEle.superCluster(); - if ( pfSuperCluster.isNonnull() && pfSuperCluster.isAvailable() ) { - for ( reco::CaloCluster_iterator pfCluster = pfSuperCluster->clustersBegin(); - pfCluster != pfSuperCluster->clustersEnd(); ++pfCluster ) { + if (pfSuperCluster.isNonnull() && pfSuperCluster.isAvailable()) { + for (reco::CaloCluster_iterator pfCluster = pfSuperCluster->clustersBegin(); + pfCluster != pfSuperCluster->clustersEnd(); + ++pfCluster) { double pfClusterEn = (*pfCluster)->energy(); - if ( pfCluster == pfSuperCluster->clustersBegin() ) ElecEe += pfClusterEn; - else ElecEgamma += pfClusterEn; + if (pfCluster == pfSuperCluster->clustersBegin()) + ElecEe += pfClusterEn; + else + ElecEgamma += pfClusterEn; } } - + Float_t ElecPin = std::sqrt(theGsfEle.trackMomentumAtVtx().Mag2()); - Float_t ElecPout = std::sqrt(theGsfEle.trackMomentumOut().Mag2()); - Float_t ElecEtotOverPin = (ElecPin > 0.0) ? ((ElecEe + ElecEgamma)/ElecPin) : -0.1; + Float_t ElecPout = std::sqrt(theGsfEle.trackMomentumOut().Mag2()); + Float_t ElecEtotOverPin = (ElecPin > 0.0) ? ((ElecEe + ElecEgamma) / ElecPin) : -0.1; Float_t ElecEecal = theGsfEle.ecalEnergy(); Float_t ElecDeltaEta = theGsfEle.deltaEtaSeedClusterTrackAtCalo(); Float_t ElecDeltaPhi = theGsfEle.deltaPhiSeedClusterTrackAtCalo(); Float_t ElecMvaInSigmaEtaEta = (theGsfEle).mvaInput().sigmaEtaEta; Float_t ElecMvaInHadEnergy = (theGsfEle).mvaInput().hadEnergy; Float_t ElecMvaInDeltaEta = (theGsfEle).mvaInput().deltaEta; - + //Variables related to the GsfTrack Float_t ElecChi2NormGSF = -99.; Float_t ElecGSFNumHits = -99.; Float_t ElecGSFTrackResol = -99.; Float_t ElecGSFTracklnPt = -99.; - if ( theGsfEle.gsfTrack().isNonnull() ) { + if (theGsfEle.gsfTrack().isNonnull()) { ElecChi2NormGSF = (theGsfEle).gsfTrack()->normalizedChi2(); ElecGSFNumHits = (theGsfEle).gsfTrack()->numberOfValidHits(); - if ( theGsfEle.gsfTrack()->pt() > 0. ) { - ElecGSFTrackResol = theGsfEle.gsfTrack()->ptError()/theGsfEle.gsfTrack()->pt(); - ElecGSFTracklnPt = log(theGsfEle.gsfTrack()->pt())*M_LN10; + if (theGsfEle.gsfTrack()->pt() > 0.) { + ElecGSFTrackResol = theGsfEle.gsfTrack()->ptError() / theGsfEle.gsfTrack()->pt(); + ElecGSFTracklnPt = log(theGsfEle.gsfTrack()->pt()) * M_LN10; } } //Variables related to the CtfTrack Float_t ElecChi2NormKF = -99.; Float_t ElecKFNumHits = -99.; - if ( theGsfEle.closestCtfTrackRef().isNonnull() ) { + if (theGsfEle.closestCtfTrackRef().isNonnull()) { ElecChi2NormKF = (theGsfEle).closestCtfTrackRef()->normalizedChi2(); ElecKFNumHits = (theGsfEle).closestCtfTrackRef()->numberOfValidHits(); } @@ -729,34 +724,38 @@ double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau, ElecMvaInDeltaEta); } -double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau) -{ +double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau) { // === tau variables === float TauEtaAtEcalEntrance = -99.; float sumEtaTimesEnergy = 0.; float sumEnergy = 0.; const std::vector& signalPFCands = thePFTau.signalPFCands(); - for ( const auto & pfCandidate : signalPFCands ) { - sumEtaTimesEnergy += pfCandidate->positionAtECALEntrance().eta()*pfCandidate->energy(); + for (const auto& pfCandidate : signalPFCands) { + sumEtaTimesEnergy += pfCandidate->positionAtECALEntrance().eta() * pfCandidate->energy(); sumEnergy += pfCandidate->energy(); } - if ( sumEnergy > 0. ) { - TauEtaAtEcalEntrance = sumEtaTimesEnergy/sumEnergy; + if (sumEnergy > 0.) { + TauEtaAtEcalEntrance = sumEtaTimesEnergy / sumEnergy; } - + float TauLeadChargedPFCandEtaAtEcalEntrance = -99.; float TauLeadChargedPFCandPt = -99.; - for ( const auto & pfCandidate : signalPFCands ) { + for (const auto& pfCandidate : signalPFCands) { const reco::Track* track = nullptr; - if ( pfCandidate->trackRef().isNonnull() ) track = pfCandidate->trackRef().get(); - else if ( pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->innerTrack().isNonnull() ) track = pfCandidate->muonRef()->innerTrack().get(); - else if ( pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->globalTrack().isNonnull() ) track = pfCandidate->muonRef()->globalTrack().get(); - else if ( pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->outerTrack().isNonnull() ) track = pfCandidate->muonRef()->outerTrack().get(); - else if ( pfCandidate->gsfTrackRef().isNonnull() ) track = pfCandidate->gsfTrackRef().get(); - if ( track ) { - if ( track->pt() > TauLeadChargedPFCandPt ) { - TauLeadChargedPFCandEtaAtEcalEntrance = pfCandidate->positionAtECALEntrance().eta(); - TauLeadChargedPFCandPt = track->pt(); + if (pfCandidate->trackRef().isNonnull()) + track = pfCandidate->trackRef().get(); + else if (pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->innerTrack().isNonnull()) + track = pfCandidate->muonRef()->innerTrack().get(); + else if (pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->globalTrack().isNonnull()) + track = pfCandidate->muonRef()->globalTrack().get(); + else if (pfCandidate->muonRef().isNonnull() && pfCandidate->muonRef()->outerTrack().isNonnull()) + track = pfCandidate->muonRef()->outerTrack().get(); + else if (pfCandidate->gsfTrackRef().isNonnull()) + track = pfCandidate->gsfTrackRef().get(); + if (track) { + if (track->pt() > TauLeadChargedPFCandPt) { + TauLeadChargedPFCandEtaAtEcalEntrance = pfCandidate->positionAtECALEntrance().eta(); + TauLeadChargedPFCandPt = track->pt(); } } } @@ -765,9 +764,9 @@ double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau) Float_t TauEmFraction = std::max(thePFTau.emFraction(), (Float_t)0.); Float_t TauLeadPFChargedHadrHoP = 0.; Float_t TauLeadPFChargedHadrEoP = 0.; - if ( thePFTau.leadPFChargedHadrCand()->p() > 0. ) { - TauLeadPFChargedHadrHoP = thePFTau.leadPFChargedHadrCand()->hcalEnergy()/thePFTau.leadPFChargedHadrCand()->p(); - TauLeadPFChargedHadrEoP = thePFTau.leadPFChargedHadrCand()->ecalEnergy()/thePFTau.leadPFChargedHadrCand()->p(); + if (thePFTau.leadPFChargedHadrCand()->p() > 0.) { + TauLeadPFChargedHadrHoP = thePFTau.leadPFChargedHadrCand()->hcalEnergy() / thePFTau.leadPFChargedHadrCand()->p(); + TauLeadPFChargedHadrEoP = thePFTau.leadPFChargedHadrCand()->ecalEnergy() / thePFTau.leadPFChargedHadrCand()->p(); } std::vector GammasdEtaInSigCone; @@ -776,20 +775,19 @@ double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau) std::vector GammasdEtaOutSigCone; std::vector GammasdPhiOutSigCone; std::vector GammasPtOutSigCone; - reco::Candidate::LorentzVector pfGammaSum(0,0,0,0); - reco::Candidate::LorentzVector pfChargedSum(0,0,0,0); - - for ( const auto & gamma : thePFTau.signalGammaCands() ) { + reco::Candidate::LorentzVector pfGammaSum(0, 0, 0, 0); + reco::Candidate::LorentzVector pfChargedSum(0, 0, 0, 0); + + for (const auto& gamma : thePFTau.signalGammaCands()) { float dR = deltaR(gamma->p4(), thePFTau.leadChargedHadrCand()->p4()); - float signalrad = std::max(0.05, std::min(0.10, 3.0/std::max(1.0, thePFTau.pt()))); + float signalrad = std::max(0.05, std::min(0.10, 3.0 / std::max(1.0, thePFTau.pt()))); // pfGammas inside the tau signal cone if (dR < signalrad) { - if ( thePFTau.leadChargedHadrCand().isNonnull() ) { + if (thePFTau.leadChargedHadrCand().isNonnull()) { GammasdEtaInSigCone.push_back(gamma->eta() - thePFTau.leadChargedHadrCand()->eta()); GammasdPhiInSigCone.push_back(gamma->phi() - thePFTau.leadChargedHadrCand()->phi()); - } - else { + } else { GammasdEtaInSigCone.push_back(gamma->eta() - thePFTau.eta()); GammasdPhiInSigCone.push_back(gamma->phi() - thePFTau.phi()); } @@ -798,28 +796,27 @@ double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau) } // pfGammas outside the tau signal cone else { - if ( thePFTau.leadChargedHadrCand().isNonnull() ) { + if (thePFTau.leadChargedHadrCand().isNonnull()) { GammasdEtaOutSigCone.push_back(gamma->eta() - thePFTau.leadChargedHadrCand()->eta()); GammasdPhiOutSigCone.push_back(gamma->phi() - thePFTau.leadChargedHadrCand()->phi()); - } - else { + } else { GammasdEtaOutSigCone.push_back(gamma->eta() - thePFTau.eta()); GammasdPhiOutSigCone.push_back(gamma->phi() - thePFTau.phi()); } GammasPtOutSigCone.push_back(gamma->pt()); } } - - for ( const auto & charged : thePFTau.signalChargedHadrCands() ) { + + for (const auto& charged : thePFTau.signalChargedHadrCands()) { float dR = deltaR(charged->p4(), thePFTau.leadChargedHadrCand()->p4()); - float signalrad = std::max(0.05, std::min(0.10, 3.0/std::max(1.0, thePFTau.pt()))); - + float signalrad = std::max(0.05, std::min(0.10, 3.0 / std::max(1.0, thePFTau.pt()))); + // charged particles inside the tau signal cone if (dR < signalrad) { - pfChargedSum += charged->p4(); + pfChargedSum += charged->p4(); } } - + Int_t TauSignalPFGammaCandsIn = GammasPtInSigCone.size(); Int_t TauSignalPFGammaCandsOut = GammasPtOutSigCone.size(); Float_t TauVisMassIn = (pfGammaSum + pfChargedSum).mass(); @@ -827,31 +824,32 @@ double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau) Float_t TauPhi = thePFTau.phi(); float sumPhiTimesEnergy = 0.; float sumEnergyPhi = 0.; - if ( !usePhiAtEcalEntranceExtrapolation_ ){ - for ( const auto & pfCandidate : signalPFCands ) { - sumPhiTimesEnergy += pfCandidate->positionAtECALEntrance().phi()*pfCandidate->energy(); + if (!usePhiAtEcalEntranceExtrapolation_) { + for (const auto& pfCandidate : signalPFCands) { + sumPhiTimesEnergy += pfCandidate->positionAtECALEntrance().phi() * pfCandidate->energy(); sumEnergyPhi += pfCandidate->energy(); } - } - else{ - TauPhi= -99.; - for ( const auto & signalPFCand : signalPFCands ) { - reco::Candidate const* signalCand = signalPFCand.get(); + } else { + TauPhi = -99.; + for (const auto& signalPFCand : signalPFCands) { + reco::Candidate const* signalCand = signalPFCand.get(); float phi = thePFTau.phi(); - math::XYZPoint aPos; - if ( atECalEntrance(signalCand, aPos) == true ) phi = aPos.Phi(); - sumPhiTimesEnergy += phi*signalCand->energy(); + bool success = false; + reco::Candidate::Point aPos = positionAtECalEntrance_(signalCand, success); + if (success) { + phi = aPos.Phi(); + } + sumPhiTimesEnergy += phi * signalCand->energy(); sumEnergy += signalCand->energy(); } } - if ( sumEnergyPhi > 0. ) { - TauPhi = sumPhiTimesEnergy/sumEnergyPhi; + if (sumEnergyPhi > 0.) { + TauPhi = sumPhiTimesEnergy / sumEnergyPhi; } Float_t TaudCrackPhi = dCrackPhi(TauPhi, TauEtaAtEcalEntrance); Float_t TaudCrackEta = dCrackEta(TauEtaAtEcalEntrance); Float_t TauHasGsf = thePFTau.leadPFChargedHadrCand()->gsfTrackRef().isNonnull(); - // === electron variables === Float_t dummyElecEta = 9.9; @@ -894,11 +892,10 @@ double AntiElectronIDMVA6::MVAValue(const reco::PFTau& thePFTau) 0.); } -double AntiElectronIDMVA6::MVAValue(const pat::Tau& theTau, const pat::Electron& theEle) -{ +double AntiElectronIDMVA6::MVAValue(const pat::Tau& theTau, const pat::Electron& theEle) { // === tau variables === float TauEtaAtEcalEntrance = theTau.etaAtEcalEntrance(); - + float TauLeadChargedPFCandEtaAtEcalEntrance = theTau.etaAtEcalEntranceLeadChargedCand(); float TauLeadChargedPFCandPt = theTau.ptLeadChargedCand(); @@ -907,9 +904,9 @@ double AntiElectronIDMVA6::MVAValue(const pat::Tau& theTau, const pat::Electron& Float_t TauEmFraction = std::max(theTau.emFraction_MVA(), (Float_t)0.); Float_t TauLeadPFChargedHadrHoP = 0.; Float_t TauLeadPFChargedHadrEoP = 0.; - if ( theTau.leadChargedHadrCand()->p() > 0. ) { - TauLeadPFChargedHadrHoP = theTau.hcalEnergyLeadChargedHadrCand()/theTau.leadChargedHadrCand()->p(); - TauLeadPFChargedHadrEoP = theTau.ecalEnergyLeadChargedHadrCand()/theTau.leadChargedHadrCand()->p(); + if (theTau.leadChargedHadrCand()->p() > 0.) { + TauLeadPFChargedHadrHoP = theTau.hcalEnergyLeadChargedHadrCand() / theTau.leadChargedHadrCand()->p(); + TauLeadPFChargedHadrEoP = theTau.ecalEnergyLeadChargedHadrCand() / theTau.leadChargedHadrCand()->p(); } std::vector GammasdEtaInSigCone; @@ -918,26 +915,25 @@ double AntiElectronIDMVA6::MVAValue(const pat::Tau& theTau, const pat::Electron& std::vector GammasdEtaOutSigCone; std::vector GammasdPhiOutSigCone; std::vector GammasPtOutSigCone; - reco::Candidate::LorentzVector pfGammaSum(0,0,0,0); - reco::Candidate::LorentzVector pfChargedSum(0,0,0,0); - + reco::Candidate::LorentzVector pfGammaSum(0, 0, 0, 0); + reco::Candidate::LorentzVector pfChargedSum(0, 0, 0, 0); + const reco::CandidatePtrVector signalGammaCands = theTau.signalGammaCands(); - for ( const auto & gamma : signalGammaCands ) { + for (const auto& gamma : signalGammaCands) { float dR = deltaR(gamma->p4(), theTau.leadChargedHadrCand()->p4()); - float signalrad = std::max(0.05, std::min(0.10, 3.0/std::max(1.0, theTau.pt()))); + float signalrad = std::max(0.05, std::min(0.10, 3.0 / std::max(1.0, theTau.pt()))); // pfGammas inside the tau signal cone if (dR < signalrad) { - if ( theTau.leadChargedHadrCand().isNonnull() ) { + if (theTau.leadChargedHadrCand().isNonnull()) { GammasdEtaInSigCone.push_back(gamma->eta() - theTau.leadChargedHadrCand()->eta()); GammasdPhiInSigCone.push_back(gamma->phi() - theTau.leadChargedHadrCand()->phi()); - //A.-C. please check whether this change is safe against future trainings + //A.-C. please check whether this change is safe against future trainings //GammasdPhiInSigCone.push_back(deltaPhi((*gamma)->phi(), theTau.leadChargedHadrCand()->phi())); - } - else { + } else { GammasdEtaInSigCone.push_back(gamma->eta() - theTau.eta()); GammasdPhiInSigCone.push_back(gamma->phi() - theTau.phi()); - //A.-C. please check whether this change is safe against future trainings + //A.-C. please check whether this change is safe against future trainings //GammasdPhiInSigCone.push_back(deltaPhi(gamma->phi(), theTau.phi())); } GammasPtInSigCone.push_back(gamma->pt()); @@ -945,108 +941,115 @@ double AntiElectronIDMVA6::MVAValue(const pat::Tau& theTau, const pat::Electron& } // pfGammas outside the tau signal cone else { - if ( theTau.leadChargedHadrCand().isNonnull() ) { + if (theTau.leadChargedHadrCand().isNonnull()) { GammasdEtaOutSigCone.push_back(gamma->eta() - theTau.leadChargedHadrCand()->eta()); GammasdPhiOutSigCone.push_back(gamma->phi() - theTau.leadChargedHadrCand()->phi()); - //A.-C. please check whether this change is safe against future trainings + //A.-C. please check whether this change is safe against future trainings //GammasdPhiOutSigCone.push_back(deltaPhi(gamma->phi(), theTau.leadChargedHadrCand()->phi())); - } - else { + } else { GammasdEtaOutSigCone.push_back(gamma->eta() - theTau.eta()); GammasdPhiOutSigCone.push_back(gamma->phi() - theTau.phi()); - //A.-C. please chaekc whether this change is safe against future trainings + //A.-C. please chaekc whether this change is safe against future trainings //GammasdPhiOutSigCone.push_back(deltaPhi(gamma->phi(), theTau.phi())); } GammasPtOutSigCone.push_back(gamma->pt()); } } - + const reco::CandidatePtrVector signalChargedCands = theTau.signalChargedHadrCands(); - for ( const auto & charged : signalChargedCands ) { + for (const auto& charged : signalChargedCands) { float dR = deltaR(charged->p4(), theTau.leadChargedHadrCand()->p4()); - float signalrad = std::max(0.05, std::min(0.10, 3.0/std::max(1.0, theTau.pt()))); - + float signalrad = std::max(0.05, std::min(0.10, 3.0 / std::max(1.0, theTau.pt()))); + // charged particles inside the tau signal cone if (dR < signalrad) { pfChargedSum += charged->p4(); } } - + Int_t TauSignalPFGammaCandsIn = GammasPtInSigCone.size(); Int_t TauSignalPFGammaCandsOut = GammasPtOutSigCone.size(); Float_t TauVisMassIn = (pfGammaSum + pfChargedSum).mass(); Float_t TauPhi = -99.; - if ( usePhiAtEcalEntranceExtrapolation_ ) { + if (usePhiAtEcalEntranceExtrapolation_) { float sumPhiTimesEnergy = 0.; float sumEnergy = 0.; const reco::CandidatePtrVector signalCands = theTau.signalCands(); - for ( const auto & signalCandPtr : signalCands ) { + for (const auto& signalCandPtr : signalCands) { reco::Candidate const* signalCand = signalCandPtr.get(); float phi = theTau.phi(); - math::XYZPoint aPos; - if ( atECalEntrance(signalCand, aPos) == true ) phi = aPos.Phi(); - sumPhiTimesEnergy += phi*signalCand->energy(); + bool success = false; + reco::Candidate::Point aPos = positionAtECalEntrance_(signalCand, success); + if (success) { + phi = aPos.Phi(); + } + sumPhiTimesEnergy += phi * signalCand->energy(); sumEnergy += signalCand->energy(); } - if ( sumEnergy > 0. ) { - TauPhi = sumPhiTimesEnergy/sumEnergy; + if (sumEnergy > 0.) { + TauPhi = sumPhiTimesEnergy / sumEnergy; } - } - else { + } else { TauPhi = theTau.phiAtEcalEntrance(); - } + } Float_t TaudCrackPhi = dCrackPhi(TauPhi, TauEtaAtEcalEntrance); - Float_t TaudCrackEta = dCrackEta(TauEtaAtEcalEntrance); - + Float_t TaudCrackEta = dCrackEta(TauEtaAtEcalEntrance); + Float_t TauHasGsf = 0; - pat::PackedCandidate const* packedLeadTauCand = dynamic_cast(theTau.leadChargedHadrCand().get()); - if( abs(packedLeadTauCand->pdgId()) == 11 ) TauHasGsf = 1; - + pat::PackedCandidate const* packedLeadTauCand = + dynamic_cast(theTau.leadChargedHadrCand().get()); + if (abs(packedLeadTauCand->pdgId()) == 11) + TauHasGsf = 1; + // === electron variables === Float_t ElecEta = theEle.eta(); Float_t ElecPhi = theEle.phi(); - + //Variables related to the electron Cluster Float_t ElecEe = 0.; Float_t ElecEgamma = 0.; reco::SuperClusterRef pfSuperCluster = theEle.superCluster(); - if ( pfSuperCluster.isNonnull() && pfSuperCluster.isAvailable() ) { - for ( reco::CaloCluster_iterator pfCluster = pfSuperCluster->clustersBegin(); pfCluster != pfSuperCluster->clustersEnd(); ++pfCluster ) { + if (pfSuperCluster.isNonnull() && pfSuperCluster.isAvailable()) { + for (reco::CaloCluster_iterator pfCluster = pfSuperCluster->clustersBegin(); + pfCluster != pfSuperCluster->clustersEnd(); + ++pfCluster) { double pfClusterEn = (*pfCluster)->energy(); - if ( pfCluster == pfSuperCluster->clustersBegin() ) ElecEe += pfClusterEn; - else ElecEgamma += pfClusterEn; + if (pfCluster == pfSuperCluster->clustersBegin()) + ElecEe += pfClusterEn; + else + ElecEgamma += pfClusterEn; } } - + Float_t ElecPin = std::sqrt(theEle.trackMomentumAtVtx().Mag2()); - Float_t ElecPout = std::sqrt(theEle.trackMomentumOut().Mag2()); - Float_t ElecEtotOverPin = (ElecPin > 0.0) ? ((ElecEe + ElecEgamma)/ElecPin) : -0.1; + Float_t ElecPout = std::sqrt(theEle.trackMomentumOut().Mag2()); + Float_t ElecEtotOverPin = (ElecPin > 0.0) ? ((ElecEe + ElecEgamma) / ElecPin) : -0.1; Float_t ElecEecal = theEle.ecalEnergy(); Float_t ElecDeltaEta = theEle.deltaEtaSeedClusterTrackAtCalo(); Float_t ElecDeltaPhi = theEle.deltaPhiSeedClusterTrackAtCalo(); Float_t ElecMvaInSigmaEtaEta = (theEle).mvaInput().sigmaEtaEta; Float_t ElecMvaInHadEnergy = (theEle).mvaInput().hadEnergy; Float_t ElecMvaInDeltaEta = (theEle).mvaInput().deltaEta; - + //Variables related to the GsfTrack Float_t ElecChi2NormGSF = -99.; Float_t ElecGSFNumHits = -99.; Float_t ElecGSFTrackResol = -99.; Float_t ElecGSFTracklnPt = -99.; - if ( theEle.gsfTrack().isNonnull() ) { + if (theEle.gsfTrack().isNonnull()) { ElecChi2NormGSF = (theEle).gsfTrack()->normalizedChi2(); ElecGSFNumHits = (theEle).gsfTrack()->numberOfValidHits(); - if ( theEle.gsfTrack()->pt() > 0. ) { - ElecGSFTrackResol = theEle.gsfTrack()->ptError()/theEle.gsfTrack()->pt(); - ElecGSFTracklnPt = log(theEle.gsfTrack()->pt())*M_LN10; + if (theEle.gsfTrack()->pt() > 0.) { + ElecGSFTrackResol = theEle.gsfTrack()->ptError() / theEle.gsfTrack()->pt(); + ElecGSFTracklnPt = log(theEle.gsfTrack()->pt()) * M_LN10; } } //Variables related to the CtfTrack Float_t ElecChi2NormKF = -99.; Float_t ElecKFNumHits = -99.; - if ( theEle.closestCtfTrackRef().isNonnull() ) { + if (theEle.closestCtfTrackRef().isNonnull()) { ElecChi2NormKF = (theEle).closestCtfTrackRef()->normalizedChi2(); ElecKFNumHits = (theEle).closestCtfTrackRef()->numberOfValidHits(); } @@ -1090,11 +1093,10 @@ double AntiElectronIDMVA6::MVAValue(const pat::Tau& theTau, const pat::Electron& ElecMvaInDeltaEta); } -double AntiElectronIDMVA6::MVAValue(const pat::Tau& theTau) -{ +double AntiElectronIDMVA6::MVAValue(const pat::Tau& theTau) { // === tau variables === float TauEtaAtEcalEntrance = theTau.etaAtEcalEntrance(); - + float TauLeadChargedPFCandEtaAtEcalEntrance = theTau.etaAtEcalEntranceLeadChargedCand(); float TauLeadChargedPFCandPt = theTau.ptLeadChargedCand(); @@ -1103,9 +1105,9 @@ double AntiElectronIDMVA6::MVAValue(const pat::Tau& theTau) Float_t TauEmFraction = std::max(theTau.emFraction_MVA(), (Float_t)0.); Float_t TauLeadPFChargedHadrHoP = 0.; Float_t TauLeadPFChargedHadrEoP = 0.; - if ( theTau.leadChargedHadrCand()->p() > 0. ) { - TauLeadPFChargedHadrHoP = theTau.hcalEnergyLeadChargedHadrCand()/theTau.leadChargedHadrCand()->p(); - TauLeadPFChargedHadrEoP = theTau.ecalEnergyLeadChargedHadrCand()/theTau.leadChargedHadrCand()->p(); + if (theTau.leadChargedHadrCand()->p() > 0.) { + TauLeadPFChargedHadrHoP = theTau.hcalEnergyLeadChargedHadrCand() / theTau.leadChargedHadrCand()->p(); + TauLeadPFChargedHadrEoP = theTau.ecalEnergyLeadChargedHadrCand() / theTau.leadChargedHadrCand()->p(); } std::vector GammasdEtaInSigCone; @@ -1114,21 +1116,20 @@ double AntiElectronIDMVA6::MVAValue(const pat::Tau& theTau) std::vector GammasdEtaOutSigCone; std::vector GammasdPhiOutSigCone; std::vector GammasPtOutSigCone; - reco::Candidate::LorentzVector pfGammaSum(0,0,0,0); - reco::Candidate::LorentzVector pfChargedSum(0,0,0,0); - + reco::Candidate::LorentzVector pfGammaSum(0, 0, 0, 0); + reco::Candidate::LorentzVector pfChargedSum(0, 0, 0, 0); + const reco::CandidatePtrVector signalGammaCands = theTau.signalGammaCands(); - for ( const auto & gamma : signalGammaCands ) { + for (const auto& gamma : signalGammaCands) { float dR = deltaR(gamma->p4(), theTau.leadChargedHadrCand()->p4()); - float signalrad = std::max(0.05, std::min(0.10, 3.0/std::max(1.0, theTau.pt()))); + float signalrad = std::max(0.05, std::min(0.10, 3.0 / std::max(1.0, theTau.pt()))); // pfGammas inside the tau signal cone if (dR < signalrad) { - if ( theTau.leadChargedHadrCand().isNonnull() ) { + if (theTau.leadChargedHadrCand().isNonnull()) { GammasdEtaInSigCone.push_back(gamma->eta() - theTau.leadChargedHadrCand()->eta()); GammasdPhiInSigCone.push_back(gamma->phi() - theTau.leadChargedHadrCand()->phi()); - } - else { + } else { GammasdEtaInSigCone.push_back(gamma->eta() - theTau.eta()); GammasdPhiInSigCone.push_back(gamma->phi() - theTau.phi()); } @@ -1137,61 +1138,64 @@ double AntiElectronIDMVA6::MVAValue(const pat::Tau& theTau) } // pfGammas outside the tau signal cone else { - if ( theTau.leadChargedHadrCand().isNonnull() ) { + if (theTau.leadChargedHadrCand().isNonnull()) { GammasdEtaOutSigCone.push_back(gamma->eta() - theTau.leadChargedHadrCand()->eta()); GammasdPhiOutSigCone.push_back(gamma->phi() - theTau.leadChargedHadrCand()->phi()); - } - else { + } else { GammasdEtaOutSigCone.push_back(gamma->eta() - theTau.eta()); GammasdPhiOutSigCone.push_back(gamma->phi() - theTau.phi()); } GammasPtOutSigCone.push_back(gamma->pt()); } } - + const reco::CandidatePtrVector signalChargedCands = theTau.signalChargedHadrCands(); - for ( const auto & charged : signalChargedCands ) { + for (const auto& charged : signalChargedCands) { float dR = deltaR(charged->p4(), theTau.leadChargedHadrCand()->p4()); - float signalrad = std::max(0.05, std::min(0.10, 3.0/std::max(1.0, theTau.pt()))); - + float signalrad = std::max(0.05, std::min(0.10, 3.0 / std::max(1.0, theTau.pt()))); + // charged particles inside the tau signal cone if (dR < signalrad) { - pfChargedSum += charged->p4(); + pfChargedSum += charged->p4(); } } - + Int_t TauSignalPFGammaCandsIn = GammasPtInSigCone.size(); Int_t TauSignalPFGammaCandsOut = GammasPtOutSigCone.size(); Float_t TauVisMassIn = (pfGammaSum + pfChargedSum).mass(); Float_t TauPhi = -99.; - if ( usePhiAtEcalEntranceExtrapolation_ ) { + if (usePhiAtEcalEntranceExtrapolation_) { float sumPhiTimesEnergy = 0.; float sumEnergy = 0.; const reco::CandidatePtrVector signalCands = theTau.signalCands(); - for ( const auto & signalCandPtr : signalCands ) { + for (const auto& signalCandPtr : signalCands) { reco::Candidate const* signalCand = signalCandPtr.get(); float phi = theTau.phi(); - math::XYZPoint aPos; - if ( atECalEntrance(signalCand, aPos) == true ) phi = aPos.Phi(); - sumPhiTimesEnergy += phi*signalCand->energy(); + bool success = false; + reco::Candidate::Point aPos = positionAtECalEntrance_(signalCand, success); + if (success) { + phi = aPos.Phi(); + } + sumPhiTimesEnergy += phi * signalCand->energy(); sumEnergy += signalCand->energy(); } - if ( sumEnergy > 0. ) { - TauPhi = sumPhiTimesEnergy/sumEnergy; + if (sumEnergy > 0.) { + TauPhi = sumPhiTimesEnergy / sumEnergy; } - } - else { + } else { TauPhi = theTau.phiAtEcalEntrance(); - } + } Float_t TaudCrackPhi = dCrackPhi(TauPhi, TauEtaAtEcalEntrance); - Float_t TaudCrackEta = dCrackEta(TauEtaAtEcalEntrance); - + Float_t TaudCrackEta = dCrackEta(TauEtaAtEcalEntrance); + Float_t TauHasGsf = 0; - pat::PackedCandidate const* packedLeadTauCand = dynamic_cast(theTau.leadChargedHadrCand().get()); + pat::PackedCandidate const* packedLeadTauCand = + dynamic_cast(theTau.leadChargedHadrCand().get()); //const reco::Track & pseudoTrack = packedLeadTauCand->pseudoTrack(); - if( abs(packedLeadTauCand->pdgId()) == 11 ) TauHasGsf = 1; - + if (abs(packedLeadTauCand->pdgId()) == 11) + TauHasGsf = 1; + // === electron variables === Float_t dummyElecEta = 9.9; @@ -1234,119 +1238,95 @@ double AntiElectronIDMVA6::MVAValue(const pat::Tau& theTau) 0.); } -double AntiElectronIDMVA6::minimum(double a, double b) -{ - if ( std::abs(b) < std::abs(a) ) return b; - else return a; +double AntiElectronIDMVA6::minimum(double a, double b) { + if (std::abs(b) < std::abs(a)) + return b; + else + return a; } namespace { // IN: define locations of the 18 phi-cracks - std::array fill_cPhi() { - constexpr double pi = M_PI; // 3.14159265358979323846; - std::array cPhi; - // IN: define locations of the 18 phi-cracks - cPhi[0] = 2.97025; - for ( unsigned iCrack = 1; iCrack <= 17; ++iCrack ) { - cPhi[iCrack] = cPhi[0] - 2.*iCrack*pi/18; - } - return cPhi; - } - - const std::array cPhi = fill_cPhi(); + std::array fill_cPhi() { + constexpr double pi = M_PI; // 3.14159265358979323846; + std::array cPhi; + // IN: define locations of the 18 phi-cracks + cPhi[0] = 2.97025; + for (unsigned iCrack = 1; iCrack <= 17; ++iCrack) { + cPhi[iCrack] = cPhi[0] - 2. * iCrack * pi / 18; + } + return cPhi; + } -} + const std::array cPhi = fill_cPhi(); -double AntiElectronIDMVA6::dCrackPhi(double phi, double eta) -{ -//--- compute the (unsigned) distance to the closest phi-crack in the ECAL barrel +} // namespace + +double AntiElectronIDMVA6::dCrackPhi(double phi, double eta) { + //--- compute the (unsigned) distance to the closest phi-crack in the ECAL barrel - constexpr double pi = M_PI; // 3.14159265358979323846; + constexpr double pi = M_PI; // 3.14159265358979323846; // IN: shift of this location if eta < 0 constexpr double delta_cPhi = 0.00638; - double retVal = 99.; - - if ( eta >= -1.47464 && eta <= 1.47464 ) { + double retVal = 99.; + if (eta >= -1.47464 && eta <= 1.47464) { // the location is shifted - if ( eta < 0. ) phi += delta_cPhi; + if (eta < 0.) + phi += delta_cPhi; // CV: need to bring-back phi into interval [-pi,+pi] - if ( phi > pi ) phi -= 2.*pi; - if ( phi < -pi ) phi += 2.*pi; - - if ( phi >= -pi && phi <= pi ) { + if (phi > pi) + phi -= 2. * pi; + if (phi < -pi) + phi += 2. * pi; + if (phi >= -pi && phi <= pi) { // the problem of the extrema: - if ( phi < cPhi[17] || phi >= cPhi[0] ) { - if ( phi < 0. ) phi += 2.*pi; - retVal = minimum(phi - cPhi[0], phi - cPhi[17] - 2.*pi); + if (phi < cPhi[17] || phi >= cPhi[0]) { + if (phi < 0.) + phi += 2. * pi; + retVal = minimum(phi - cPhi[0], phi - cPhi[17] - 2. * pi); } else { - // between these extrema... - bool OK = false; - unsigned iCrack = 16; - while( !OK ) { - if ( phi < cPhi[iCrack] ) { - retVal = minimum(phi - cPhi[iCrack + 1], phi - cPhi[iCrack]); - OK = true; - } else { - iCrack -= 1; - } - } + // between these extrema... + bool OK = false; + unsigned iCrack = 16; + while (!OK) { + if (phi < cPhi[iCrack]) { + retVal = minimum(phi - cPhi[iCrack + 1], phi - cPhi[iCrack]); + OK = true; + } else { + iCrack -= 1; + } + } } } else { - retVal = 0.; // IN: if there is a problem, we assume that we are in a crack + retVal = 0.; // IN: if there is a problem, we assume that we are in a crack } } else { - return -99.; + return -99.; } - + return std::abs(retVal); } -double AntiElectronIDMVA6::dCrackEta(double eta) -{ -//--- compute the (unsigned) distance to the closest eta-crack in the ECAL barrel - +double AntiElectronIDMVA6::dCrackEta(double eta) { + //--- compute the (unsigned) distance to the closest eta-crack in the ECAL barrel + // IN: define locations of the eta-cracks - double cracks[5] = { 0., 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00 }; - + double cracks[5] = {0., 4.44747e-01, 7.92824e-01, 1.14090e+00, 1.47464e+00}; + double retVal = 99.; - - for ( int iCrack = 0; iCrack < 5 ; ++iCrack ) { + + for (int iCrack = 0; iCrack < 5; ++iCrack) { double d = minimum(eta - cracks[iCrack], eta + cracks[iCrack]); - if ( std::abs(d) < std::abs(retVal) ) { + if (std::abs(d) < std::abs(retVal)) { retVal = d; } } return std::abs(retVal); } - -bool AntiElectronIDMVA6::atECalEntrance(const reco::Candidate* part, math::XYZPoint &pos) -{ - bool result = false; - BaseParticlePropagator theParticle = - BaseParticlePropagator(RawParticle(math::XYZTLorentzVector(part->px(), - part->py(), - part->pz(), - part->energy()), - math::XYZTLorentzVector(part->vertex().x(), - part->vertex().y(), - part->vertex().z(), - 0.), - part->charge()), - 0.,0.,bField_); - theParticle.propagateToEcalEntrance(false); - if(theParticle.getSuccess()!=0){ - pos = math::XYZPoint(theParticle.particle().vertex()); - result = true; - } - else { - result = false; - } - return result; -} diff --git a/RecoTauTag/RecoTau/src/PositionAtECalEntranceComputer.cc b/RecoTauTag/RecoTau/src/PositionAtECalEntranceComputer.cc new file mode 100644 index 0000000000000..aedab0dab18fb --- /dev/null +++ b/RecoTauTag/RecoTau/src/PositionAtECalEntranceComputer.cc @@ -0,0 +1,37 @@ +#include "RecoTauTag/RecoTau/interface/PositionAtECalEntranceComputer.h" + +#include "MagneticField/Engine/interface/MagneticField.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h" + +#include + +PositionAtECalEntranceComputer::PositionAtECalEntranceComputer() : bField_z_(-1.) {} + +PositionAtECalEntranceComputer::~PositionAtECalEntranceComputer() {} + +void PositionAtECalEntranceComputer::beginEvent(const edm::EventSetup& es) { + edm::ESHandle bField; + es.get().get(bField); + bField_z_ = bField->inTesla(GlobalPoint(0., 0., 0.)).z(); +} + +reco::Candidate::Point PositionAtECalEntranceComputer::operator()(const reco::Candidate* particle, bool& success) const { + assert(bField_z_ != -1.); + BaseParticlePropagator propagator = BaseParticlePropagator( + RawParticle(particle->p4(), + math::XYZTLorentzVector(particle->vertex().x(), particle->vertex().y(), particle->vertex().z(), 0.), + particle->charge()), + 0., + 0., + bField_z_); + propagator.propagateToEcalEntrance(false); + reco::Candidate::Point position; + if (propagator.getSuccess() != 0) { + position = propagator.particle().vertex().Vect(); + success = true; + } else { + success = false; + } + return position; +} diff --git a/RecoTauTag/RecoTau/src/TauDiscriminationProducerBase.cc b/RecoTauTag/RecoTau/src/TauDiscriminationProducerBase.cc index 81de9ffad163b..1e36ac828ea70 100644 --- a/RecoTauTag/RecoTau/src/TauDiscriminationProducerBase.cc +++ b/RecoTauTag/RecoTau/src/TauDiscriminationProducerBase.cc @@ -18,8 +18,8 @@ TauDiscriminationProducerBase::TauDiscriminationProdu : moduleLabel_(iConfig.getParameter("@module_label")) { // tau collection to discriminate - TauProducer_ = iConfig.getParameter(getProducerString()); - Tau_token= consumes(TauProducer_); + TauProducer_ = iConfig.getParameter(getTauTypeString() + "Producer"); + Tau_token = consumes(TauProducer_); // prediscriminant operator // require the tau to pass the following prediscriminants @@ -159,7 +159,7 @@ void TauDiscriminationProducerBase::produce(edm::Even template void TauDiscriminationProducerBase::fillProducerDescriptions(edm::ParameterSetDescription& desc) { // helper function, it fills the description of the Producers parameter - desc.add(getProducerString(),edm::InputTag("fixme")); + desc.add(getTauTypeString() + "Producer", edm::InputTag("fixme")); { edm::ParameterSetDescription pset_prediscriminants; pset_prediscriminants.add("BooleanOperator","AND"); @@ -178,10 +178,17 @@ void TauDiscriminationProducerBase::fillProducerDescr } } -// template specialiazation to get the correct (Calo/PF)TauProducer names -template<> std::string getProducerString() { return "PFTauProducer"; } -template<> std::string getProducerString() { return "CaloTauProducer"; } -template<> std::string getProducerString() { return "PATTauProducer"; } +template +std::string TauDiscriminationProducerBase::getTauTypeString() { + if(std::is_same::value) + return "PFTau"; + if(std::is_same::value) + return "CaloTau"; + if(std::is_same::value) + return "PATTau"; + throw cms::Exception("TauDiscriminationProducerBase") + << "Unsupported TauType used. You must use either PFTau, PATTau or CaloTau."; +} // compile our desired types and make available to linker template class TauDiscriminationProducerBase;