diff --git a/RecoTauTag/HLTProducers/python/deepTauAtHLT.py b/RecoTauTag/HLTProducers/python/deepTauAtHLT.py new file mode 100644 index 0000000000000..5bc006fc64035 --- /dev/null +++ b/RecoTauTag/HLTProducers/python/deepTauAtHLT.py @@ -0,0 +1,140 @@ +import FWCore.ParameterSet.Config as cms + +from RecoTauTag.RecoTau.PFRecoTauDiscriminationByIsolation_cfi import * +from RecoTauTag.RecoTau.PFRecoTauQualityCuts_cfi import PFTauQualityCuts +from RecoTauTag.RecoTau.PFRecoTauDiscriminationByHPSSelection_cfi import hpsSelectionDiscriminator + +from RecoTauTag.RecoTau.PFTauPrimaryVertexProducer_cfi import * +from RecoTauTag.RecoTau.PFTauSecondaryVertexProducer_cfi import * +from RecoTauTag.RecoTau.PFTauTransverseImpactParameters_cfi import * + +from RecoTauTag.RecoTau.DeepTau_cfi import * + +from RecoTauTag.RecoTau.PFRecoTauPFJetInputs_cfi import PFRecoTauPFJetInputs +## DeltaBeta correction factor +ak4dBetaCorrection = 0.20 + +def update(process): + process.options.wantSummary = cms.untracked.bool(True) + + process.hltFixedGridRhoFastjetAll = cms.EDProducer( "FixedGridRhoProducerFastjet", + gridSpacing = cms.double( 0.55 ), + maxRapidity = cms.double( 5.0 ), + pfCandidatesTag = cms.InputTag( "hltParticleFlowReg" ) + ) + + PFTauQualityCuts.primaryVertexSrc = cms.InputTag("hltPixelVertices") + + ## Decay mode prediscriminant + requireDecayMode = cms.PSet( + BooleanOperator = cms.string("and"), + decayMode = cms.PSet( + Producer = cms.InputTag('hltHpsPFTauDiscriminationByDecayModeFindingNewDMsReg'), + cut = cms.double(0.5) + ) + ) + + ## Cut based isolations dR=0.5 + process.hpsPFTauBasicDiscriminators = pfRecoTauDiscriminationByIsolation.clone( + PFTauProducer = 'hltHpsPFTauProducerReg', + Prediscriminants = requireDecayMode.clone(), + deltaBetaPUTrackPtCutOverride = True, # Set the boolean = True to override. + deltaBetaPUTrackPtCutOverride_val = 0.5, # Set the value for new value. + particleFlowSrc = 'hltParticleFlowReg', + vertexSrc = PFTauQualityCuts.primaryVertexSrc, + customOuterCone = PFRecoTauPFJetInputs.isolationConeSize, + isoConeSizeForDeltaBeta = 0.8, + deltaBetaFactor = "%0.4f"%(ak4dBetaCorrection), + qualityCuts = dict(isolationQualityCuts = dict(minTrackHits = 3, minGammaEt = 1.0, minTrackPt = 0.5)), + IDdefinitions = [ + cms.PSet( + IDname = cms.string("ChargedIsoPtSum"), + ApplyDiscriminationByTrackerIsolation = cms.bool(True), + storeRawSumPt = cms.bool(True) + ), + cms.PSet( + IDname = cms.string("NeutralIsoPtSum"), + ApplyDiscriminationByECALIsolation = cms.bool(True), + storeRawSumPt = cms.bool(True) + ), + cms.PSet( + IDname = cms.string("NeutralIsoPtSumWeight"), + ApplyDiscriminationByWeightedECALIsolation = cms.bool(True), + storeRawSumPt = cms.bool(True), + UseAllPFCandsForWeights = cms.bool(True) + ), + cms.PSet( + IDname = cms.string("TauFootprintCorrection"), + storeRawFootprintCorrection = cms.bool(True) + ), + cms.PSet( + IDname = cms.string("PhotonPtSumOutsideSignalCone"), + storeRawPhotonSumPt_outsideSignalCone = cms.bool(True) + ), + cms.PSet( + IDname = cms.string("PUcorrPtSum"), + applyDeltaBetaCorrection = cms.bool(True), + storeRawPUsumPt = cms.bool(True) + ), + ], + ) + + ## Cut based isolations dR=0.3 + process.hpsPFTauBasicDiscriminatorsdR03 = process.hpsPFTauBasicDiscriminators.clone( + customOuterCone = 0.3 + ) + + process.hpsPFTauPrimaryVertexProducer = PFTauPrimaryVertexProducer.clone( + PFTauTag = "hltHpsPFTauProducerReg", + ElectronTag = "hltEgammaCandidates", + MuonTag = "hltMuonsReg", + PVTag = "hltPixelVertices", + beamSpot = "hltOnlineBeamSpot", + discriminators = [ + cms.PSet( + discriminator = cms.InputTag('hltHpsPFTauDiscriminationByDecayModeFindingNewDMsReg'), + selectionCut = cms.double(0.5) + ) + ], + cut = "pt > 18.0 & abs(eta) < 2.4", + qualityCuts = PFTauQualityCuts + ) + + process.hpsPFTauSecondaryVertexProducer = PFTauSecondaryVertexProducer.clone( + PFTauTag = "hltHpsPFTauProducerReg" + ) + process.hpsPFTauTransverseImpactParameters = PFTauTransverseImpactParameters.clone( + PFTauTag = "hltHpsPFTauProducerReg", + PFTauPVATag = "hpsPFTauPrimaryVertexProducer", + PFTauSVATag = "hpsPFTauSecondaryVertexProducer", + useFullCalculation = True + ) + + file_names = [ + 'core:RecoTauTag/TrainingFiles/data/DeepTauId/deepTau_2017v2p6_e6_core.pb', + 'inner:RecoTauTag/TrainingFiles/data/DeepTauId/deepTau_2017v2p6_e6_inner.pb', + 'outer:RecoTauTag/TrainingFiles/data/DeepTauId/deepTau_2017v2p6_e6_outer.pb', + ] + + working_points = ["0.", "0.92"] + + process.deepTauProducer = DeepTau.clone( + taus = 'hltHpsPFTauProducerReg', + pfcands = 'hltParticleFlowReg', + vertices = 'hltPixelVertices', + rho = 'hltFixedGridRhoFastjetAll', + graph_file = file_names, + disable_dxy_pca = cms.bool(True), + is_online = cms.bool(True), + basicTauDiscriminators = 'hpsPFTauBasicDiscriminators', + basicTauDiscriminatorsdR03 = 'hpsPFTauBasicDiscriminatorsdR03', + Prediscriminants = requireDecayMode.clone(), + VSeWP = working_points, + VSmuWP = working_points, + VSjetWP = working_points + ) + + # Add DeepTauProducer + process.HLTHPSMediumChargedIsoPFTauSequenceReg += (process.hpsPFTauPrimaryVertexProducer + process.hpsPFTauSecondaryVertexProducer + process.hpsPFTauTransverseImpactParameters + process.hltFixedGridRhoFastjetAll + process.hpsPFTauBasicDiscriminators + process.hpsPFTauBasicDiscriminatorsdR03 + process.deepTauProducer) + + return process diff --git a/RecoTauTag/RecoTau/interface/DeepTauBase.h b/RecoTauTag/RecoTau/interface/DeepTauBase.h index 6dba2956316cb..6f487d102d880 100644 --- a/RecoTauTag/RecoTau/interface/DeepTauBase.h +++ b/RecoTauTag/RecoTau/interface/DeepTauBase.h @@ -19,18 +19,26 @@ #include "DataFormats/PatCandidates/interface/Muon.h" #include "DataFormats/PatCandidates/interface/Tau.h" #include "DataFormats/TauReco/interface/TauDiscriminatorContainer.h" +#include "DataFormats/TauReco/interface/PFTauDiscriminator.h" +#include "DataFormats/PatCandidates/interface/PATTauDiscriminator.h" #include "CommonTools/Utils/interface/StringObjectFunction.h" #include "RecoTauTag/RecoTau/interface/PFRecoTauClusterVariables.h" #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" #include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "DataFormats/Common/interface/View.h" +#include "DataFormats/Common/interface/RefToBase.h" +#include "DataFormats/Provenance/interface/ProductProvenance.h" +#include "DataFormats/Provenance/interface/ProcessHistoryID.h" +#include "FWCore/Common/interface/Provenance.h" #include +#include namespace deep_tau { class TauWPThreshold { public: explicit TauWPThreshold(const std::string& cut_str); - double operator()(const pat::Tau& tau) const; + double operator()(const reco::BaseTau& tau, bool isPFTau) const; private: std::unique_ptr fn_; @@ -57,9 +65,9 @@ namespace deep_tau { class DeepTauBase : public edm::stream::EDProducer> { public: - using TauType = pat::Tau; using TauDiscriminator = reco::TauDiscriminatorContainer; - using TauCollection = std::vector; + using TauCollection = edm::View; + using CandidateCollection = edm::View; using TauRef = edm::Ref; using TauRefProd = edm::RefProd; using ElectronCollection = pat::ElectronCollection; @@ -76,7 +84,8 @@ namespace deep_tau { std::unique_ptr get_value(const edm::Handle& taus, const tensorflow::Tensor& pred, - const WPList& working_points) const; + const WPList* working_points, + bool is_online) const; }; using OutputCollection = std::map; @@ -89,19 +98,45 @@ namespace deep_tau { static std::unique_ptr initializeGlobalCache(const edm::ParameterSet& cfg); static void globalEndJob(const DeepTauCache* cache) {} + template + struct TauDiscInfo { + edm::InputTag label; + edm::Handle handle; + edm::EDGetTokenT disc_token; + double cut; + void fill(const edm::Event& evt) { evt.getByToken(disc_token, handle); } + }; + + // select boolean operation on prediscriminants (and = 0x01, or = 0x00) + uint8_t andPrediscriminants_; + std::vector> patPrediscriminants_; + std::vector> recoPrediscriminants_; + + enum BasicDiscriminator { + ChargedIsoPtSum, + NeutralIsoPtSum, + NeutralIsoPtSumWeight, + FootprintCorrection, + PhotonPtSumOutsideSignalCone, + PUcorrPtSum + }; + private: - virtual tensorflow::Tensor getPredictions(edm::Event& event, - const edm::EventSetup& es, - edm::Handle taus) = 0; + virtual tensorflow::Tensor getPredictions(edm::Event& event, edm::Handle taus) = 0; virtual void createOutputs(edm::Event& event, const tensorflow::Tensor& pred, edm::Handle taus); protected: edm::EDGetTokenT tausToken_; - edm::EDGetTokenT pfcandToken_; + edm::EDGetTokenT pfcandToken_; edm::EDGetTokenT vtxToken_; std::map workingPoints_; + const bool is_online_; OutputCollection outputs_; const DeepTauCache* cache_; + + static const std::map stringFromDiscriminator_; + static const std::vector requiredBasicDiscriminators_; + static const std::vector requiredBasicDiscriminatorsdR03_; }; } // namespace deep_tau diff --git a/RecoTauTag/RecoTau/plugins/DPFIsolation.cc b/RecoTauTag/RecoTau/plugins/DPFIsolation.cc index 0b87390d32a2e..13b281cc948bf 100644 --- a/RecoTauTag/RecoTau/plugins/DPFIsolation.cc +++ b/RecoTauTag/RecoTau/plugins/DPFIsolation.cc @@ -74,9 +74,7 @@ class DPFIsolation : public deep_tau::DeepTauBase { } private: - tensorflow::Tensor getPredictions(edm::Event& event, - const edm::EventSetup& es, - edm::Handle taus) override { + tensorflow::Tensor getPredictions(edm::Event& event, edm::Handle taus) override { edm::Handle pfcands; event.getByToken(pfcandToken_, pfcands); diff --git a/RecoTauTag/RecoTau/plugins/DeepTauId.cc b/RecoTauTag/RecoTau/plugins/DeepTauId.cc index 9e707b50955e1..b33fabd88d1b3 100644 --- a/RecoTauTag/RecoTau/plugins/DeepTauId.cc +++ b/RecoTauTag/RecoTau/plugins/DeepTauId.cc @@ -4,10 +4,15 @@ * Tau identification using Deep NN. * * \author Konstantin Androsov, INFN Pisa + * Christian Veelken, Tallinn */ #include "RecoTauTag/RecoTau/interface/DeepTauBase.h" #include "FWCore/Utilities/interface/isFinite.h" +#include "DataFormats/TauReco/interface/PFTauTransverseImpactParameterAssociation.h" + +#include +#include "tbb/concurrent_unordered_set.h" namespace deep_tau { constexpr int NumberOfOutputs = 4; @@ -414,6 +419,260 @@ namespace { } } // namespace dnn_inputs_2017_v2 + float getTauID(const pat::Tau& tau, const std::string& tauID, float default_value = -999., bool assert_input = true) { + static tbb::concurrent_unordered_set isFirstWarning; + if (tau.isTauIDAvailable(tauID)) { + return tau.tauID(tauID); + } else { + if (assert_input) { + throw cms::Exception("DeepTauId") + << "Exception in : No tauID '" << tauID << "' available in pat::Tau given as function argument."; + } + if (isFirstWarning.insert(tauID).second) { + edm::LogWarning("DeepTauID") << "Warning in : No tauID '" << tauID + << "' available in pat::Tau given as function argument." + << " Using default_value = " << default_value << " instead." << std::endl; + } + return default_value; + } + } + + struct TauFunc { + const reco::TauDiscriminatorContainer* basicTauDiscriminatorCollection; + const reco::TauDiscriminatorContainer* basicTauDiscriminatordR03Collection; + const edm::AssociationVector>* + pfTauTransverseImpactParameters; + + using BasicDiscr = deep_tau::DeepTauBase::BasicDiscriminator; + std::map indexMap; + std::map indexMapdR03; + + const float getChargedIsoPtSum(const reco::PFTau& tau, const edm::RefToBase tau_ref) const { + return (*basicTauDiscriminatorCollection)[tau_ref].rawValues.at(indexMap.at(BasicDiscr::ChargedIsoPtSum)); + } + const float getChargedIsoPtSum(const pat::Tau& tau, const edm::RefToBase tau_ref) const { + return getTauID(tau, "chargedIsoPtSum"); + } + const float getChargedIsoPtSumdR03(const reco::PFTau& tau, const edm::RefToBase tau_ref) const { + return (*basicTauDiscriminatordR03Collection)[tau_ref].rawValues.at(indexMapdR03.at(BasicDiscr::ChargedIsoPtSum)); + } + const float getChargedIsoPtSumdR03(const pat::Tau& tau, const edm::RefToBase tau_ref) const { + return getTauID(tau, "chargedIsoPtSumdR03"); + } + const float getFootprintCorrectiondR03(const reco::PFTau& tau, const edm::RefToBase tau_ref) const { + return (*basicTauDiscriminatordR03Collection)[tau_ref].rawValues.at( + indexMapdR03.at(BasicDiscr::FootprintCorrection)); + } + const float getFootprintCorrectiondR03(const pat::Tau& tau, const edm::RefToBase tau_ref) const { + return getTauID(tau, "footprintCorrectiondR03"); + } + const float getNeutralIsoPtSum(const reco::PFTau& tau, const edm::RefToBase tau_ref) const { + return (*basicTauDiscriminatorCollection)[tau_ref].rawValues.at(indexMap.at(BasicDiscr::NeutralIsoPtSum)); + } + const float getNeutralIsoPtSum(const pat::Tau& tau, const edm::RefToBase tau_ref) const { + return getTauID(tau, "neutralIsoPtSum"); + } + const float getNeutralIsoPtSumdR03(const reco::PFTau& tau, const edm::RefToBase tau_ref) const { + return (*basicTauDiscriminatordR03Collection)[tau_ref].rawValues.at(indexMapdR03.at(BasicDiscr::NeutralIsoPtSum)); + } + const float getNeutralIsoPtSumdR03(const pat::Tau& tau, const edm::RefToBase tau_ref) const { + return getTauID(tau, "neutralIsoPtSumdR03"); + } + const float getNeutralIsoPtSumWeight(const reco::PFTau& tau, const edm::RefToBase tau_ref) const { + return (*basicTauDiscriminatorCollection)[tau_ref].rawValues.at(indexMap.at(BasicDiscr::NeutralIsoPtSumWeight)); + } + const float getNeutralIsoPtSumWeight(const pat::Tau& tau, const edm::RefToBase tau_ref) const { + return getTauID(tau, "neutralIsoPtSumWeight"); + } + const float getNeutralIsoPtSumdR03Weight(const reco::PFTau& tau, + const edm::RefToBase tau_ref) const { + return (*basicTauDiscriminatordR03Collection)[tau_ref].rawValues.at( + indexMapdR03.at(BasicDiscr::NeutralIsoPtSumWeight)); + } + const float getNeutralIsoPtSumdR03Weight(const pat::Tau& tau, const edm::RefToBase tau_ref) const { + return getTauID(tau, "neutralIsoPtSumWeightdR03"); + } + const float getPhotonPtSumOutsideSignalCone(const reco::PFTau& tau, + const edm::RefToBase tau_ref) const { + return (*basicTauDiscriminatorCollection)[tau_ref].rawValues.at( + indexMap.at(BasicDiscr::PhotonPtSumOutsideSignalCone)); + } + const float getPhotonPtSumOutsideSignalCone(const pat::Tau& tau, + const edm::RefToBase tau_ref) const { + return getTauID(tau, "photonPtSumOutsideSignalCone"); + } + const float getPhotonPtSumOutsideSignalConedR03(const reco::PFTau& tau, + const edm::RefToBase tau_ref) const { + return (*basicTauDiscriminatordR03Collection)[tau_ref].rawValues.at( + indexMapdR03.at(BasicDiscr::PhotonPtSumOutsideSignalCone)); + } + const float getPhotonPtSumOutsideSignalConedR03(const pat::Tau& tau, + const edm::RefToBase tau_ref) const { + return getTauID(tau, "photonPtSumOutsideSignalConedR03"); + } + const float getPuCorrPtSum(const reco::PFTau& tau, const edm::RefToBase tau_ref) const { + return (*basicTauDiscriminatorCollection)[tau_ref].rawValues.at(indexMap.at(BasicDiscr::PUcorrPtSum)); + } + const float getPuCorrPtSum(const pat::Tau& tau, const edm::RefToBase tau_ref) const { + return getTauID(tau, "puCorrPtSum"); + } + + auto getdxyPCA(const reco::PFTau& tau, const size_t tau_index) const { + return pfTauTransverseImpactParameters->value(tau_index)->dxy_PCA(); + } + auto getdxyPCA(const pat::Tau& tau, const size_t tau_index) const { return tau.dxy_PCA(); } + auto getdxy(const reco::PFTau& tau, const size_t tau_index) const { + return pfTauTransverseImpactParameters->value(tau_index)->dxy(); + } + auto getdxy(const pat::Tau& tau, const size_t tau_index) const { return tau.dxy(); } + auto getdxyError(const reco::PFTau& tau, const size_t tau_index) const { + return pfTauTransverseImpactParameters->value(tau_index)->dxy_error(); + } + auto getdxyError(const pat::Tau& tau, const size_t tau_index) const { return tau.dxy_error(); } + auto getdxySig(const reco::PFTau& tau, const size_t tau_index) const { + return pfTauTransverseImpactParameters->value(tau_index)->dxy_Sig(); + } + auto getdxySig(const pat::Tau& tau, const size_t tau_index) const { return tau.dxy_Sig(); } + auto getip3d(const reco::PFTau& tau, const size_t tau_index) const { + return pfTauTransverseImpactParameters->value(tau_index)->ip3d(); + } + auto getip3d(const pat::Tau& tau, const size_t tau_index) const { return tau.ip3d(); } + auto getip3dError(const reco::PFTau& tau, const size_t tau_index) const { + return pfTauTransverseImpactParameters->value(tau_index)->ip3d_error(); + } + auto getip3dError(const pat::Tau& tau, const size_t tau_index) const { return tau.ip3d_error(); } + auto getip3dSig(const reco::PFTau& tau, const size_t tau_index) const { + return pfTauTransverseImpactParameters->value(tau_index)->ip3d_Sig(); + } + auto getip3dSig(const pat::Tau& tau, const size_t tau_index) const { return tau.ip3d_Sig(); } + auto getHasSecondaryVertex(const reco::PFTau& tau, const size_t tau_index) const { + return pfTauTransverseImpactParameters->value(tau_index)->hasSecondaryVertex(); + } + auto getHasSecondaryVertex(const pat::Tau& tau, const size_t tau_index) const { return tau.hasSecondaryVertex(); } + auto getFlightLength(const reco::PFTau& tau, const size_t tau_index) const { + return pfTauTransverseImpactParameters->value(tau_index)->flightLength(); + } + auto getFlightLength(const pat::Tau& tau, const size_t tau_index) const { return tau.flightLength(); } + auto getFlightLengthSig(const reco::PFTau& tau, const size_t tau_index) const { + return pfTauTransverseImpactParameters->value(tau_index)->flightLengthSig(); + } + auto getFlightLengthSig(const pat::Tau& tau, const size_t tau_index) const { return tau.flightLengthSig(); } + + auto getLeadingTrackNormChi2(const reco::PFTau& tau) { return reco::tau::lead_track_chi2(tau); } + auto getLeadingTrackNormChi2(const pat::Tau& tau) { return tau.leadingTrackNormChi2(); } + auto getEmFraction(const pat::Tau& tau) { return tau.emFraction_MVA(); } + auto getEmFraction(const reco::PFTau& tau) { return tau.emFraction(); } + auto getEtaAtEcalEntrance(const pat::Tau& tau) { return tau.etaAtEcalEntranceLeadChargedCand(); } + auto getEtaAtEcalEntrance(const reco::PFTau& tau) { + return tau.leadPFChargedHadrCand()->positionAtECALEntrance().eta(); + } + auto getEcalEnergyLeadingChargedHadr(const reco::PFTau& tau) { return tau.leadPFChargedHadrCand()->ecalEnergy(); } + auto getEcalEnergyLeadingChargedHadr(const pat::Tau& tau) { return tau.ecalEnergyLeadChargedHadrCand(); } + auto getHcalEnergyLeadingChargedHadr(const reco::PFTau& tau) { return tau.leadPFChargedHadrCand()->hcalEnergy(); } + auto getHcalEnergyLeadingChargedHadr(const pat::Tau& tau) { return tau.hcalEnergyLeadChargedHadrCand(); } + + template + bool passPrediscriminants(const PreDiscrType prediscriminants, + const size_t andPrediscriminants, + const edm::RefToBase tau_ref) { + bool passesPrediscriminants = (andPrediscriminants ? 1 : 0); + // check tau passes prediscriminants + size_t nPrediscriminants = prediscriminants.size(); + for (size_t iDisc = 0; iDisc < nPrediscriminants; ++iDisc) { + // current discriminant result for this tau + double discResult = (*prediscriminants[iDisc].handle)[tau_ref]; + uint8_t thisPasses = (discResult > prediscriminants[iDisc].cut) ? 1 : 0; + + // if we are using the AND option, as soon as one fails, + // the result is FAIL and we can quit looping. + // if we are using the OR option as soon as one passes, + // the result is pass and we can quit looping + + // truth table + // | result (thisPasses) + // | F | T + //----------------------------------- + // AND(T) | res=fails | continue + // | break | + //----------------------------------- + // OR (F) | continue | res=passes + // | | break + + if (thisPasses ^ andPrediscriminants) //XOR + { + passesPrediscriminants = (andPrediscriminants ? 0 : 1); //NOR + break; + } + } + return passesPrediscriminants; + } + }; + + namespace candFunc { + auto getTauDz(const reco::PFCandidate& cand) { return cand.bestTrack()->dz(); } + auto getTauDz(const pat::PackedCandidate& cand) { return cand.dz(); } + auto getTauDZSigValid(const reco::PFCandidate& cand) { + return cand.bestTrack() != nullptr && std::isnormal(cand.bestTrack()->dz()) && std::isnormal(cand.dzError()) && + cand.dzError() > 0; + } + auto getTauDZSigValid(const pat::PackedCandidate& cand) { + return cand.hasTrackDetails() && std::isnormal(cand.dz()) && std::isnormal(cand.dzError()) && cand.dzError() > 0; + } + auto getTauDxy(const reco::PFCandidate& cand) { return cand.bestTrack()->dxy(); } + auto getTauDxy(const pat::PackedCandidate& cand) { return cand.dxy(); } + auto getPvAssocationQuality(const reco::PFCandidate& cand) { return 0.7013f; } + auto getPvAssocationQuality(const pat::PackedCandidate& cand) { return cand.pvAssociationQuality(); } + auto getPuppiWeight(const reco::PFCandidate& cand, const float aod_value) { return aod_value; } + auto getPuppiWeight(const pat::PackedCandidate& cand, const float aod_value) { return cand.puppiWeight(); } + auto getPuppiWeightNoLep(const reco::PFCandidate& cand, const float aod_value) { return aod_value; } + auto getPuppiWeightNoLep(const pat::PackedCandidate& cand, const float aod_value) { + return cand.puppiWeightNoLep(); + } + auto getLostInnerHits(const reco::PFCandidate& cand, float default_value) { + return cand.bestTrack() != nullptr + ? cand.bestTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) + : default_value; + } + auto getLostInnerHits(const pat::PackedCandidate& cand, float default_value) { return cand.lostInnerHits(); } + auto getNumberOfPixelHits(const reco::PFCandidate& cand, float default_value) { + return cand.bestTrack() != nullptr + ? cand.bestTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) + : default_value; + } + auto getNumberOfPixelHits(const pat::PackedCandidate& cand, float default_value) { + return cand.numberOfPixelHits(); + } + auto getHasTrackDetails(const reco::PFCandidate& cand) { return cand.bestTrack() != nullptr; } + auto getHasTrackDetails(const pat::PackedCandidate& cand) { return cand.hasTrackDetails(); } + auto getPseudoTrack(const reco::PFCandidate& cand) { return *cand.bestTrack(); } + auto getPseudoTrack(const pat::PackedCandidate& cand) { return cand.pseudoTrack(); } + auto getFromPV(const reco::PFCandidate& cand) { return 0.9994f; } + auto getFromPV(const pat::PackedCandidate& cand) { return cand.fromPV(); } + auto getHCalFraction(const reco::PFCandidate& cand, bool disable_hcalFraction_workaround) { + return cand.rawHcalEnergy() / (cand.rawHcalEnergy() + cand.rawEcalEnergy()); + } + auto getHCalFraction(const pat::PackedCandidate& cand, bool disable_hcalFraction_workaround) { + float hcal_fraction = 0.; + if (disable_hcalFraction_workaround) { + // CV: use consistent definition for pfCand_chHad_hcalFraction + // in DeepTauId.cc code and in TauMLTools/Production/plugins/TauTupleProducer.cc + hcal_fraction = cand.hcalFraction(); + } else { + // CV: backwards compatibility with DeepTau training v2p1 used during Run 2 + if (cand.pdgId() == 1 || cand.pdgId() == 130) { + hcal_fraction = cand.hcalFraction(); + } else if (cand.isIsolatedChargedHadron()) { + hcal_fraction = cand.rawHcalFraction(); + } + } + return hcal_fraction; + } + auto getRawCaloFraction(const reco::PFCandidate& cand) { + return (cand.rawEcalEnergy() + cand.rawHcalEnergy()) / cand.energy(); + } + auto getRawCaloFraction(const pat::PackedCandidate& cand) { return cand.rawCaloFraction(); } + }; // namespace candFunc + template float dEta(const LVector1& p4, const LVector2& tau_p4) { return static_cast(p4.eta() - tau_p4.eta()); @@ -441,7 +700,7 @@ namespace { n_hits[MuonSubdetId::RPC].assign(n_muon_stations, 0); } - void addMatchedMuon(const pat::Muon& muon, const pat::Tau& tau) { + void addMatchedMuon(const pat::Muon& muon, reco::BaseTau const& tau) { static constexpr int n_stations = 4; ++n_muons; @@ -484,8 +743,9 @@ namespace { } } - static std::vector findMatchedMuons(const pat::Tau& tau, - const pat::MuonCollection& muons, + template + static std::vector findMatchedMuons(const TauCastType& tau, + const std::vector* muons, double deltaR, double minPt) { const reco::Muon* hadr_cand_muon = nullptr; @@ -493,7 +753,7 @@ namespace { hadr_cand_muon = tau.leadPFChargedHadrCand()->muonRef().get(); std::vector matched_muons; const double dR2 = deltaR * deltaR; - for (const pat::Muon& muon : muons) { + for (const pat::Muon& muon : *muons) { const reco::Muon* reco_muon = &muon; if (muon.pt() <= minPt) continue; @@ -506,8 +766,8 @@ namespace { return matched_muons; } - template - void fillTensor(const TensorElemGet& get, const pat::Tau& tau, float default_value) const { + template + void fillTensor(const TensorElemGet& get, const TauCastType& tau, float default_value) const { get(dnn::n_matched_muons) = n_muons; get(dnn::muon_pt) = best_matched_muon != nullptr ? best_matched_muon->p4().pt() : default_value; get(dnn::muon_dEta) = best_matched_muon != nullptr ? dEta(best_matched_muon->p4(), tau.p4()) : default_value; @@ -723,7 +983,7 @@ namespace { } template <> - CellObjectType GetCellObjectType(const pat::PackedCandidate& cand) { + CellObjectType GetCellObjectType(reco::Candidate const& cand) { static const std::map obj_types = {{11, CellObjectType::PfCand_electron}, {13, CellObjectType::PfCand_muon}, {22, CellObjectType::PfCand_gamma}, @@ -752,12 +1012,17 @@ namespace { using Map = std::map; using const_iterator = Map::const_iterator; - CellGrid(unsigned n_cells_eta, unsigned n_cells_phi, double cell_size_eta, double cell_size_phi) + CellGrid(unsigned n_cells_eta, + unsigned n_cells_phi, + double cell_size_eta, + double cell_size_phi, + bool disable_CellIndex_workaround) : nCellsEta(n_cells_eta), nCellsPhi(n_cells_phi), nTotal(nCellsEta * nCellsPhi), cellSizeEta(cell_size_eta), - cellSizePhi(cell_size_phi) { + cellSizePhi(cell_size_phi), + disable_CellIndex_workaround_(disable_CellIndex_workaround) { if (nCellsEta % 2 != 1 || nCellsEta < 1) throw cms::Exception("DeepTauId") << "Invalid number of eta cells."; if (nCellsPhi % 2 != 1 || nCellsPhi < 1) @@ -774,11 +1039,19 @@ namespace { int getPhiTensorIndex(const CellIndex& cellIndex) const { return cellIndex.phi + maxPhiIndex(); } bool tryGetCellIndex(double deltaEta, double deltaPhi, CellIndex& cellIndex) const { - static auto getCellIndex = [](double x, double maxX, double size, int& index) { + const auto getCellIndex = [this](double x, double maxX, double size, int& index) { const double absX = std::abs(x); if (absX > maxX) return false; - const double absIndex = std::floor(std::abs(absX / size - 0.5)); + double absIndex; + if (disable_CellIndex_workaround_) { + // CV: use consistent definition for CellIndex + // in DeepTauId.cc code and new DeepTau trainings + absIndex = std::floor(absX / size + 0.5); + } else { + // CV: backwards compatibility with DeepTau training v2p1 used during Run 2 + absIndex = std::floor(std::abs(absX / size - 0.5)); + } index = static_cast(std::copysign(absIndex, x)); return true; }; @@ -801,10 +1074,29 @@ namespace { private: std::map cells; + const bool disable_CellIndex_workaround_; }; - } // anonymous namespace +using bd = deep_tau::DeepTauBase::BasicDiscriminator; +const std::map deep_tau::DeepTauBase::stringFromDiscriminator_{ + {bd::ChargedIsoPtSum, "ChargedIsoPtSum"}, + {bd::NeutralIsoPtSum, "NeutralIsoPtSum"}, + {bd::NeutralIsoPtSumWeight, "NeutralIsoPtSumWeight"}, + {bd::FootprintCorrection, "TauFootprintCorrection"}, + {bd::PhotonPtSumOutsideSignalCone, "PhotonPtSumOutsideSignalCone"}, + {bd::PUcorrPtSum, "PUcorrPtSum"}}; +const std::vector deep_tau::DeepTauBase::requiredBasicDiscriminators_ = {bd::ChargedIsoPtSum, + bd::NeutralIsoPtSum, + bd::NeutralIsoPtSumWeight, + bd::PhotonPtSumOutsideSignalCone, + bd::PUcorrPtSum}; +const std::vector deep_tau::DeepTauBase::requiredBasicDiscriminatorsdR03_ = {bd::ChargedIsoPtSum, + bd::NeutralIsoPtSum, + bd::NeutralIsoPtSumWeight, + bd::PhotonPtSumOutsideSignalCone, + bd::FootprintCorrection}; + class DeepTauId : public deep_tau::DeepTauBase { public: static constexpr float default_value = -999.; @@ -819,6 +1111,38 @@ class DeepTauId : public deep_tau::DeepTauBase { return outputs_; } + const std::map matchDiscriminatorIndices( + edm::Event& event, + edm::EDGetTokenT discriminatorContainerToken, + std::vector requiredDiscr) { + std::map discrIndexMapStr; + auto const aHandle = event.getHandle(discriminatorContainerToken); + auto const aProv = aHandle.provenance(); + if (aProv == nullptr) + aHandle.whyFailed()->raise(); + const auto& psetsFromProvenance = edm::parameterSet(*aProv, event.processHistory()); + auto const idlist = psetsFromProvenance.getParameter>("IDdefinitions"); + for (size_t j = 0; j < idlist.size(); ++j) { + std::string idname = idlist[j].getParameter("IDname"); + if (discrIndexMapStr.count(idname)) { + throw cms::Exception("DeepTauId") + << "basic discriminator " << idname << " appears more than once in the input."; + } + discrIndexMapStr[idname] = j; + } + + //translate to a map of and check if all discriminators are present + std::map discrIndexMap; + for (size_t i = 0; i < requiredDiscr.size(); i++) { + if (discrIndexMapStr.find(stringFromDiscriminator_.at(requiredDiscr[i])) == discrIndexMapStr.end()) + throw cms::Exception("DeepTauId") << "Basic Discriminator " << stringFromDiscriminator_.at(requiredDiscr[i]) + << " was not provided in the config file."; + else + discrIndexMap[requiredDiscr[i]] = discrIndexMapStr[stringFromDiscriminator_.at(requiredDiscr[i])]; + } + return discrIndexMap; + } + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; desc.add("electrons", edm::InputTag("slimmedElectrons")); @@ -833,30 +1157,63 @@ class DeepTauId : public deep_tau::DeepTauBase { desc.add("version", 2); desc.add("debug_level", 0); desc.add("disable_dxy_pca", false); + desc.add("disable_hcalFraction_workaround", false); + desc.add("disable_CellIndex_workaround", false); + desc.add("save_inputs", false); + desc.add("is_online", false); desc.add>("VSeWP"); desc.add>("VSmuWP"); desc.add>("VSjetWP"); + + desc.addUntracked("basicTauDiscriminators", edm::InputTag("basicTauDiscriminators")); + desc.addUntracked("basicTauDiscriminatorsdR03", edm::InputTag("basicTauDiscriminatorsdR03")); + desc.add("pfTauTransverseImpactParameters", edm::InputTag("hpsPFTauTransverseImpactParameters")); + + { + edm::ParameterSetDescription pset_Prediscriminants; + pset_Prediscriminants.add("BooleanOperator", "and"); + { + edm::ParameterSetDescription psd1; + psd1.add("cut"); + psd1.add("Producer"); + pset_Prediscriminants.addOptional("decayMode", psd1); + } + desc.add("Prediscriminants", pset_Prediscriminants); + } + descriptions.add("DeepTau", desc); } public: explicit DeepTauId(const edm::ParameterSet& cfg, const deep_tau::DeepTauCache* cache) : DeepTauBase(cfg, GetOutputs(), cache), - electrons_token_(consumes(cfg.getParameter("electrons"))), - muons_token_(consumes(cfg.getParameter("muons"))), + electrons_token_(consumes>(cfg.getParameter("electrons"))), + muons_token_(consumes>(cfg.getParameter("muons"))), rho_token_(consumes(cfg.getParameter("rho"))), - version(cfg.getParameter("version")), + basicTauDiscriminators_inputToken_(consumes( + cfg.getUntrackedParameter("basicTauDiscriminators"))), + basicTauDiscriminatorsdR03_inputToken_(consumes( + cfg.getUntrackedParameter("basicTauDiscriminatorsdR03"))), + pfTauTransverseImpactParameters_token_( + consumes>>( + cfg.getParameter("pfTauTransverseImpactParameters"))), + version_(cfg.getParameter("version")), debug_level(cfg.getParameter("debug_level")), - disable_dxy_pca_(cfg.getParameter("disable_dxy_pca")) { - if (version == 1) { + disable_dxy_pca_(cfg.getParameter("disable_dxy_pca")), + disable_hcalFraction_workaround_(cfg.getParameter("disable_hcalFraction_workaround")), + disable_CellIndex_workaround_(cfg.getParameter("disable_CellIndex_workaround")), + save_inputs_(cfg.getParameter("save_inputs")), + json_file_(nullptr), + file_counter_(0) { + if (version_ == 1) { input_layer_ = cache_->getGraph().node(0).name(); output_layer_ = cache_->getGraph().node(cache_->getGraph().node_size() - 1).name(); const auto& shape = cache_->getGraph().node(0).attr().at("shape").shape(); if (shape.dim(1).size() != dnn_inputs_2017v1::NumberOfInputs) throw cms::Exception("DeepTauId") << "number of inputs does not match the expected inputs for the given version"; - } else if (version == 2) { + } else if (version_ == 2) { tauBlockTensor_ = std::make_unique( tensorflow::DT_FLOAT, tensorflow::TensorShape{1, dnn_inputs_2017_v2::TauBlockInputs::NumberOfInputs}); for (size_t n = 0; n < 2; ++n) { @@ -885,7 +1242,7 @@ class DeepTauId : public deep_tau::DeepTauBase { setCellConvFeatures(*zeroOutputTensor_[is_inner], getPartialPredictions(is_inner), 0, 0, 0); } } else { - throw cms::Exception("DeepTauId") << "version " << version << " is not supported."; + throw cms::Exception("DeepTauId") << "version " << version_ << " is not supported."; } } @@ -920,6 +1277,8 @@ class DeepTauId : public deep_tau::DeepTauBase { return std::clamp(norm_value, -n_sigmas_max, n_sigmas_max); } + static bool isAbove(double value, double min) { return std::isnormal(value) && value > min; } + static bool calculateElectronClusterVarsV2(const pat::Electron& ele, float& cc_ele_energy, float& cc_gamma_energy, @@ -943,38 +1302,169 @@ class DeepTauId : public deep_tau::DeepTauBase { return false; } - inline void checkInputs( - const tensorflow::Tensor& inputs, const char* block_name, int n_inputs, int n_eta = 1, int n_phi = 1) const { + inline void checkInputs(const tensorflow::Tensor& inputs, + const std::string& block_name, + int n_inputs, + const CellGrid* grid = nullptr) const { if (debug_level >= 1) { - for (int eta = 0; eta < n_eta; ++eta) { - for (int phi = 0; phi < n_phi; phi++) { - for (int k = 0; k < n_inputs; ++k) { - const float input = - n_eta == 1 && n_phi == 1 ? inputs.matrix()(0, k) : inputs.tensor()(0, eta, phi, k); - if (edm::isNotFinite(input)) - throw cms::Exception("DeepTauId") - << "in the " << block_name << ", input is not finite, i.e. infinite or NaN, for eta_index = " << n_eta - << ", phi_index = " << n_phi << ", input_index = " << k; - if (debug_level >= 2) - std::cout << block_name << "," << eta << "," << phi << "," << k << "," << std::setprecision(5) - << std::fixed << input << '\n'; + std::cout << ": block_name = " << block_name << std::endl; + if (block_name == "input_tau") { + for (int input_index = 0; input_index < n_inputs; ++input_index) { + float input = inputs.matrix()(0, input_index); + if (edm::isNotFinite(input)) { + throw cms::Exception("DeepTauId") + << "in the " << block_name + << ", input is not finite, i.e. infinite or NaN, for input_index = " << input_index; + } + if (debug_level >= 2) { + std::cout << block_name << "[var = " << input_index << "] = " << std::setprecision(5) << std::fixed << input + << std::endl; + } + } + } else { + assert(grid); + int n_eta, n_phi; + if (block_name.find("input_inner") != std::string::npos) { + n_eta = 5; + n_phi = 5; + } else if (block_name.find("input_outer") != std::string::npos) { + n_eta = 10; + n_phi = 10; + } else + assert(0); + int eta_phi_index = 0; + for (int eta = -n_eta; eta <= n_eta; ++eta) { + for (int phi = -n_phi; phi <= n_phi; ++phi) { + const CellIndex cell_index{eta, phi}; + const auto cell_iter = grid->find(cell_index); + if (cell_iter != grid->end()) { + for (int input_index = 0; input_index < n_inputs; ++input_index) { + float input = inputs.tensor()(eta_phi_index, 0, 0, input_index); + if (edm::isNotFinite(input)) { + throw cms::Exception("DeepTauId") + << "in the " << block_name << ", input is not finite, i.e. infinite or NaN, for eta = " << eta + << ", phi = " << phi << ", input_index = " << input_index; + } + if (debug_level >= 2) { + std::cout << block_name << "[eta = " << eta << "][phi = " << phi << "][var = " << input_index + << "] = " << std::setprecision(5) << std::fixed << input << std::endl; + } + } + eta_phi_index += 1; + } } } } } } + inline void saveInputs(const tensorflow::Tensor& inputs, + const std::string& block_name, + int n_inputs, + const CellGrid* grid = nullptr) { + if (debug_level >= 1) { + std::cout << ": block_name = " << block_name << std::endl; + } + if (!is_first_block_) + (*json_file_) << ", "; + (*json_file_) << "\"" << block_name << "\": ["; + if (block_name == "input_tau") { + for (int input_index = 0; input_index < n_inputs; ++input_index) { + float input = inputs.matrix()(0, input_index); + if (input_index != 0) + (*json_file_) << ", "; + (*json_file_) << input; + } + } else { + assert(grid); + int n_eta, n_phi; + if (block_name.find("input_inner") != std::string::npos) { + n_eta = 5; + n_phi = 5; + } else if (block_name.find("input_outer") != std::string::npos) { + n_eta = 10; + n_phi = 10; + } else + assert(0); + int eta_phi_index = 0; + for (int eta = -n_eta; eta <= n_eta; ++eta) { + if (eta != -n_eta) + (*json_file_) << ", "; + (*json_file_) << "["; + for (int phi = -n_phi; phi <= n_phi; ++phi) { + if (phi != -n_phi) + (*json_file_) << ", "; + (*json_file_) << "["; + const CellIndex cell_index{eta, phi}; + const auto cell_iter = grid->find(cell_index); + for (int input_index = 0; input_index < n_inputs; ++input_index) { + float input = 0.; + if (cell_iter != grid->end()) { + input = inputs.tensor()(eta_phi_index, 0, 0, input_index); + } + if (input_index != 0) + (*json_file_) << ", "; + (*json_file_) << input; + } + if (cell_iter != grid->end()) { + eta_phi_index += 1; + } + (*json_file_) << "]"; + } + (*json_file_) << "]"; + } + } + (*json_file_) << "]"; + is_first_block_ = false; + } + private: - tensorflow::Tensor getPredictions(edm::Event& event, - const edm::EventSetup& es, - edm::Handle taus) override { - edm::Handle electrons; - event.getByToken(electrons_token_, electrons); + tensorflow::Tensor getPredictions(edm::Event& event, edm::Handle taus) override { + // Empty dummy vectors + const std::vector electron_collection_default; + const std::vector muon_collection_default; + const reco::TauDiscriminatorContainer basicTauDiscriminators_default; + const reco::TauDiscriminatorContainer basicTauDiscriminatorsdR03_default; + const edm::AssociationVector> + pfTauTransverseImpactParameters_default; + + const std::vector* electron_collection; + const std::vector* muon_collection; + const reco::TauDiscriminatorContainer* basicTauDiscriminators; + const reco::TauDiscriminatorContainer* basicTauDiscriminatorsdR03; + const edm::AssociationVector>* + pfTauTransverseImpactParameters; + + if (!is_online_) { + electron_collection = &event.get(electrons_token_); + muon_collection = &event.get(muons_token_); + pfTauTransverseImpactParameters = &pfTauTransverseImpactParameters_default; + basicTauDiscriminators = &basicTauDiscriminators_default; + basicTauDiscriminatorsdR03 = &basicTauDiscriminatorsdR03_default; + } else { + electron_collection = &electron_collection_default; + muon_collection = &muon_collection_default; + pfTauTransverseImpactParameters = &event.get(pfTauTransverseImpactParameters_token_); + basicTauDiscriminators = &event.get(basicTauDiscriminators_inputToken_); + basicTauDiscriminatorsdR03 = &event.get(basicTauDiscriminatorsdR03_inputToken_); + + // Get indices for discriminators + if (!discrIndicesMapped_) { + basicDiscrIndexMap_ = + matchDiscriminatorIndices(event, basicTauDiscriminators_inputToken_, requiredBasicDiscriminators_); + basicDiscrdR03IndexMap_ = + matchDiscriminatorIndices(event, basicTauDiscriminatorsdR03_inputToken_, requiredBasicDiscriminatorsdR03_); + discrIndicesMapped_ = true; + } + } - edm::Handle muons; - event.getByToken(muons_token_, muons); + TauFunc tauIDs = {basicTauDiscriminators, + basicTauDiscriminatorsdR03, + pfTauTransverseImpactParameters, + basicDiscrIndexMap_, + basicDiscrdR03IndexMap_}; - edm::Handle pfCands; + edm::Handle> pfCands; event.getByToken(pfcandToken_, pfCands); edm::Handle vertices; @@ -984,49 +1474,174 @@ class DeepTauId : public deep_tau::DeepTauBase { event.getByToken(rho_token_, rho); tensorflow::Tensor predictions(tensorflow::DT_FLOAT, {static_cast(taus->size()), deep_tau::NumberOfOutputs}); + for (size_t tau_index = 0; tau_index < taus->size(); ++tau_index) { + const edm::RefToBase tauRef = taus->refAt(tau_index); + std::vector pred_vector; - if (version == 1) - getPredictionsV1(taus->at(tau_index), *electrons, *muons, pred_vector); - else if (version == 2) - getPredictionsV2(taus->at(tau_index), *electrons, *muons, *pfCands, vertices->at(0), *rho, pred_vector); - else - throw cms::Exception("DeepTauId") << "version " << version << " is not supported."; - for (int k = 0; k < deep_tau::NumberOfOutputs; ++k) { - const float pred = pred_vector[0].flat()(k); - if (!(pred >= 0 && pred <= 1)) - throw cms::Exception("DeepTauId") - << "invalid prediction = " << pred << " for tau_index = " << tau_index << ", pred_index = " << k; - predictions.matrix()(tau_index, k) = pred; + + bool passesPrediscriminants; + if (is_online_) { + passesPrediscriminants = tauIDs.passPrediscriminants>>( + recoPrediscriminants_, andPrediscriminants_, tauRef); + } else { + passesPrediscriminants = tauIDs.passPrediscriminants>>( + patPrediscriminants_, andPrediscriminants_, tauRef); + } + + if (passesPrediscriminants) { + if (version_ == 1) { + if (is_online_) + getPredictionsV1( + taus->at(tau_index), tau_index, tauRef, electron_collection, muon_collection, pred_vector, tauIDs); + else + getPredictionsV1( + taus->at(tau_index), tau_index, tauRef, electron_collection, muon_collection, pred_vector, tauIDs); + } else if (version_ == 2) { + if (is_online_) { + getPredictionsV2(taus->at(tau_index), + tau_index, + tauRef, + electron_collection, + muon_collection, + *pfCands, + vertices->at(0), + *rho, + pred_vector, + tauIDs); + } else + getPredictionsV2(taus->at(tau_index), + tau_index, + tauRef, + electron_collection, + muon_collection, + *pfCands, + vertices->at(0), + *rho, + pred_vector, + tauIDs); + } else { + throw cms::Exception("DeepTauId") << "version " << version_ << " is not supported."; + } + + for (int k = 0; k < deep_tau::NumberOfOutputs; ++k) { + const float pred = pred_vector[0].flat()(k); + if (!(pred >= 0 && pred <= 1)) + throw cms::Exception("DeepTauId") + << "invalid prediction = " << pred << " for tau_index = " << tau_index << ", pred_index = " << k; + predictions.matrix()(tau_index, k) = pred; + } } } return predictions; } - void getPredictionsV1(const TauType& tau, - const pat::ElectronCollection& electrons, - const pat::MuonCollection& muons, - std::vector& pred_vector) { - const tensorflow::Tensor& inputs = createInputsV1(tau, electrons, muons); + template + void getPredictionsV1(TauCollection::const_reference& tau, + const size_t tau_index, + const edm::RefToBase tau_ref, + const std::vector* electrons, + const std::vector* muons, + std::vector& pred_vector, + TauFunc tau_funcs) { + const tensorflow::Tensor& inputs = createInputsV1( + dynamic_cast(tau), tau_index, tau_ref, electrons, muons, tau_funcs); tensorflow::run(&(cache_->getSession()), {{input_layer_, inputs}}, {output_layer_}, &pred_vector); } - void getPredictionsV2(const TauType& tau, - const pat::ElectronCollection& electrons, - const pat::MuonCollection& muons, - const pat::PackedCandidateCollection& pfCands, + template + void getPredictionsV2(TauCollection::const_reference& tau, + const size_t tau_index, + const edm::RefToBase tau_ref, + const std::vector* electrons, + const std::vector* muons, + const edm::View& pfCands, const reco::Vertex& pv, double rho, - std::vector& pred_vector) { - CellGrid inner_grid(dnn_inputs_2017_v2::number_of_inner_cell, dnn_inputs_2017_v2::number_of_inner_cell, 0.02, 0.02); - CellGrid outer_grid(dnn_inputs_2017_v2::number_of_outer_cell, dnn_inputs_2017_v2::number_of_outer_cell, 0.05, 0.05); - fillGrids(tau, electrons, inner_grid, outer_grid); - fillGrids(tau, muons, inner_grid, outer_grid); - fillGrids(tau, pfCands, inner_grid, outer_grid); - - createTauBlockInputs(tau, pv, rho); - createConvFeatures(tau, pv, rho, electrons, muons, pfCands, inner_grid, true); - createConvFeatures(tau, pv, rho, electrons, muons, pfCands, outer_grid, false); + std::vector& pred_vector, + TauFunc tau_funcs) { + if (debug_level >= 2) { + std::cout << ":" << std::endl; + std::cout << " tau: pT = " << tau.pt() << ", eta = " << tau.eta() << ", phi = " << tau.phi() << std::endl; + } + CellGrid inner_grid(dnn_inputs_2017_v2::number_of_inner_cell, + dnn_inputs_2017_v2::number_of_inner_cell, + 0.02, + 0.02, + disable_CellIndex_workaround_); + CellGrid outer_grid(dnn_inputs_2017_v2::number_of_outer_cell, + dnn_inputs_2017_v2::number_of_outer_cell, + 0.05, + 0.05, + disable_CellIndex_workaround_); + fillGrids(dynamic_cast(tau), *electrons, inner_grid, outer_grid); + fillGrids(dynamic_cast(tau), *muons, inner_grid, outer_grid); + fillGrids(dynamic_cast(tau), pfCands, inner_grid, outer_grid); + + createTauBlockInputs( + dynamic_cast(tau), tau_index, tau_ref, pv, rho, tau_funcs); + using namespace dnn_inputs_2017_v2; + checkInputs(*tauBlockTensor_, "input_tau", TauBlockInputs::NumberOfInputs); + createConvFeatures(dynamic_cast(tau), + tau_index, + tau_ref, + pv, + rho, + electrons, + muons, + pfCands, + inner_grid, + tau_funcs, + true); + checkInputs(*eGammaTensor_[true], "input_inner_egamma", EgammaBlockInputs::NumberOfInputs, &inner_grid); + checkInputs(*muonTensor_[true], "input_inner_muon", MuonBlockInputs::NumberOfInputs, &inner_grid); + checkInputs(*hadronsTensor_[true], "input_inner_hadrons", HadronBlockInputs::NumberOfInputs, &inner_grid); + createConvFeatures(dynamic_cast(tau), + tau_index, + tau_ref, + pv, + rho, + electrons, + muons, + pfCands, + outer_grid, + tau_funcs, + false); + checkInputs(*eGammaTensor_[false], "input_outer_egamma", EgammaBlockInputs::NumberOfInputs, &outer_grid); + checkInputs(*muonTensor_[false], "input_outer_muon", MuonBlockInputs::NumberOfInputs, &outer_grid); + checkInputs(*hadronsTensor_[false], "input_outer_hadrons", HadronBlockInputs::NumberOfInputs, &outer_grid); + + if (save_inputs_) { + std::string json_file_name = "DeepTauId_" + std::to_string(file_counter_) + ".json"; + json_file_ = new std::ofstream(json_file_name.data()); + is_first_block_ = true; + (*json_file_) << "{"; + saveInputs(*tauBlockTensor_, "input_tau", dnn_inputs_2017_v2::TauBlockInputs::NumberOfInputs); + saveInputs(*eGammaTensor_[true], + "input_inner_egamma", + dnn_inputs_2017_v2::EgammaBlockInputs::NumberOfInputs, + &inner_grid); + saveInputs( + *muonTensor_[true], "input_inner_muon", dnn_inputs_2017_v2::MuonBlockInputs::NumberOfInputs, &inner_grid); + saveInputs(*hadronsTensor_[true], + "input_inner_hadrons", + dnn_inputs_2017_v2::HadronBlockInputs::NumberOfInputs, + &inner_grid); + saveInputs(*eGammaTensor_[false], + "input_outer_egamma", + dnn_inputs_2017_v2::EgammaBlockInputs::NumberOfInputs, + &outer_grid); + saveInputs( + *muonTensor_[false], "input_outer_muon", dnn_inputs_2017_v2::MuonBlockInputs::NumberOfInputs, &outer_grid); + saveInputs(*hadronsTensor_[false], + "input_outer_hadrons", + dnn_inputs_2017_v2::HadronBlockInputs::NumberOfInputs, + &outer_grid); + (*json_file_) << "}"; + delete json_file_; + ++file_counter_; + } tensorflow::run(&(cache_->getSession("core")), {{"input_tau", *tauBlockTensor_}, @@ -1034,10 +1649,30 @@ class DeepTauId : public deep_tau::DeepTauBase { {"input_outer", *convTensor_.at(false)}}, {"main_output/Softmax"}, &pred_vector); + if (debug_level >= 1) { + std::cout << "output = { "; + for (int idx = 0; idx < deep_tau::NumberOfOutputs; ++idx) { + if (idx > 0) + std::cout << ", "; + std::string label; + if (idx == 0) + label = "e"; + else if (idx == 1) + label = "mu"; + else if (idx == 2) + label = "tau"; + else if (idx == 3) + label = "jet"; + else + assert(0); + std::cout << label << " = " << pred_vector[0].flat()(idx); + } + std::cout << " }" << std::endl; + } } - template - void fillGrids(const TauType& tau, const Collection& objects, CellGrid& inner_grid, CellGrid& outer_grid) { + template + void fillGrids(const TauCastType& tau, const Collection& objects, CellGrid& inner_grid, CellGrid& outer_grid) { static constexpr double outer_dR2 = 0.25; //0.5^2 const double inner_radius = getInnerSignalConeRadius(tau.polarP4().pt()); const double inner_dR2 = std::pow(inner_radius, 2); @@ -1097,16 +1732,22 @@ class DeepTauId : public deep_tau::DeepTauBase { return pred_vector.at(0); } - void createConvFeatures(const TauType& tau, + template + void createConvFeatures(const TauCastType& tau, + const size_t tau_index, + const edm::RefToBase tau_ref, const reco::Vertex& pv, double rho, - const pat::ElectronCollection& electrons, - const pat::MuonCollection& muons, - const pat::PackedCandidateCollection& pfCands, + const std::vector* electrons, + const std::vector* muons, + const edm::View& pfCands, const CellGrid& grid, + TauFunc tau_funcs, bool is_inner) { + if (debug_level >= 2) { + std::cout << ":" << std::endl; + } tensorflow::Tensor& convTensor = *convTensor_.at(is_inner); - eGammaTensor_[is_inner] = std::make_unique( tensorflow::DT_FLOAT, tensorflow::TensorShape{ @@ -1127,20 +1768,34 @@ class DeepTauId : public deep_tau::DeepTauBase { unsigned idx = 0; for (int eta = -grid.maxEtaIndex(); eta <= grid.maxEtaIndex(); ++eta) { for (int phi = -grid.maxPhiIndex(); phi <= grid.maxPhiIndex(); ++phi) { + if (debug_level >= 2) { + std::cout << "processing ( eta = " << eta << ", phi = " << phi << " )" << std::endl; + } const CellIndex cell_index{eta, phi}; const auto cell_iter = grid.find(cell_index); if (cell_iter != grid.end()) { + if (debug_level >= 2) { + std::cout << " creating inputs for ( eta = " << eta << ", phi = " << phi << " ): idx = " << idx + << std::endl; + } const Cell& cell = cell_iter->second; - createEgammaBlockInputs(idx, tau, pv, rho, electrons, pfCands, cell, is_inner); - createMuonBlockInputs(idx, tau, pv, rho, muons, pfCands, cell, is_inner); - createHadronsBlockInputs(idx, tau, pv, rho, pfCands, cell, is_inner); + createEgammaBlockInputs( + idx, tau, tau_index, tau_ref, pv, rho, electrons, pfCands, cell, tau_funcs, is_inner); + createMuonBlockInputs( + idx, tau, tau_index, tau_ref, pv, rho, muons, pfCands, cell, tau_funcs, is_inner); + createHadronsBlockInputs( + idx, tau, tau_index, tau_ref, pv, rho, pfCands, cell, tau_funcs, is_inner); idx += 1; + } else { + if (debug_level >= 2) { + std::cout << " skipping creation of inputs, because ( eta = " << eta << ", phi = " << phi + << " ) is not in the grid !!" << std::endl; + } } } } const auto predTensor = getPartialPredictions(is_inner); - idx = 0; for (int eta = -grid.maxEtaIndex(); eta <= grid.maxEtaIndex(); ++eta) { for (int phi = -grid.maxPhiIndex(); phi <= grid.maxPhiIndex(); ++phi) { @@ -1164,11 +1819,18 @@ class DeepTauId : public deep_tau::DeepTauBase { unsigned batch_idx, int eta_index, int phi_index) { - for (int n = 0; n < dnn_inputs_2017_v2::number_of_conv_features; ++n) + for (int n = 0; n < dnn_inputs_2017_v2::number_of_conv_features; ++n) { convTensor.tensor()(0, eta_index, phi_index, n) = features.tensor()(batch_idx, 0, 0, n); + } } - void createTauBlockInputs(const TauType& tau, const reco::Vertex& pv, double rho) { + template + void createTauBlockInputs(const TauCastType& tau, + const size_t& tau_index, + const edm::RefToBase tau_ref, + const reco::Vertex& pv, + double rho, + TauFunc tau_funcs) { namespace dnn = dnn_inputs_2017_v2::TauBlockInputs; tensorflow::Tensor& inputs = *tauBlockTensor_; @@ -1176,7 +1838,7 @@ class DeepTauId : public deep_tau::DeepTauBase { const auto& get = [&](int var_index) -> float& { return inputs.matrix()(0, var_index); }; - auto leadChargedHadrCand = dynamic_cast(tau.leadChargedHadrCand().get()); + auto leadChargedHadrCand = dynamic_cast(tau.leadChargedHadrCand().get()); get(dnn::rho) = getValueNorm(rho, 21.49f, 9.713f); get(dnn::tau_pt) = getValueLinear(tau.polarP4().pt(), 20.f, 1000.f, true); @@ -1187,59 +1849,61 @@ class DeepTauId : public deep_tau::DeepTauBase { get(dnn::tau_charge) = getValue(tau.charge()); get(dnn::tau_n_charged_prongs) = getValueLinear(tau.decayMode() / 5 + 1, 1, 3, true); get(dnn::tau_n_neutral_prongs) = getValueLinear(tau.decayMode() % 5, 0, 2, true); - get(dnn::chargedIsoPtSum) = getValueNorm(tau.tauID("chargedIsoPtSum"), 47.78f, 123.5f); - get(dnn::chargedIsoPtSumdR03_over_dR05) = getValue(tau.tauID("chargedIsoPtSumdR03") / tau.tauID("chargedIsoPtSum")); - get(dnn::footprintCorrection) = getValueNorm(tau.tauID("footprintCorrectiondR03"), 9.029f, 26.42f); - get(dnn::neutralIsoPtSum) = getValueNorm(tau.tauID("neutralIsoPtSum"), 57.59f, 155.3f); + get(dnn::chargedIsoPtSum) = getValueNorm(tau_funcs.getChargedIsoPtSum(tau, tau_ref), 47.78f, 123.5f); + get(dnn::chargedIsoPtSumdR03_over_dR05) = + getValue(tau_funcs.getChargedIsoPtSumdR03(tau, tau_ref) / tau_funcs.getChargedIsoPtSum(tau, tau_ref)); + get(dnn::footprintCorrection) = getValueNorm(tau_funcs.getFootprintCorrectiondR03(tau, tau_ref), 9.029f, 26.42f); + get(dnn::neutralIsoPtSum) = getValueNorm(tau_funcs.getNeutralIsoPtSum(tau, tau_ref), 57.59f, 155.3f); get(dnn::neutralIsoPtSumWeight_over_neutralIsoPtSum) = - getValue(tau.tauID("neutralIsoPtSumWeight") / tau.tauID("neutralIsoPtSum")); + getValue(tau_funcs.getNeutralIsoPtSumWeight(tau, tau_ref) / tau_funcs.getNeutralIsoPtSum(tau, tau_ref)); get(dnn::neutralIsoPtSumWeightdR03_over_neutralIsoPtSum) = - getValue(tau.tauID("neutralIsoPtSumWeightdR03") / tau.tauID("neutralIsoPtSum")); - get(dnn::neutralIsoPtSumdR03_over_dR05) = getValue(tau.tauID("neutralIsoPtSumdR03") / tau.tauID("neutralIsoPtSum")); + getValue(tau_funcs.getNeutralIsoPtSumdR03Weight(tau, tau_ref) / tau_funcs.getNeutralIsoPtSum(tau, tau_ref)); + get(dnn::neutralIsoPtSumdR03_over_dR05) = + getValue(tau_funcs.getNeutralIsoPtSumdR03(tau, tau_ref) / tau_funcs.getNeutralIsoPtSum(tau, tau_ref)); get(dnn::photonPtSumOutsideSignalCone) = - getValueNorm(tau.tauID("photonPtSumOutsideSignalConedR03"), 1.731f, 6.846f); - get(dnn::puCorrPtSum) = getValueNorm(tau.tauID("puCorrPtSum"), 22.38f, 16.34f); + getValueNorm(tau_funcs.getPhotonPtSumOutsideSignalCone(tau, tau_ref), 1.731f, 6.846f); + get(dnn::puCorrPtSum) = getValueNorm(tau_funcs.getPuCorrPtSum(tau, tau_ref), 22.38f, 16.34f); // The global PCA coordinates were used as inputs during the NN training, but it was decided to disable // them for the inference, because modeling of dxy_PCA in MC poorly describes the data, and x and y coordinates // in data results outside of the expected 5 std. dev. input validity range. On the other hand, // these coordinates are strongly era-dependent. Kept as comment to document what NN expects. if (!disable_dxy_pca_) { - get(dnn::tau_dxy_pca_x) = getValueNorm(tau.dxy_PCA().x(), -0.0241f, 0.0074f); - get(dnn::tau_dxy_pca_y) = getValueNorm(tau.dxy_PCA().y(), 0.0675f, 0.0128f); - get(dnn::tau_dxy_pca_z) = getValueNorm(tau.dxy_PCA().z(), 0.7973f, 3.456f); + auto const pca = tau_funcs.getdxyPCA(tau, tau_index); + get(dnn::tau_dxy_pca_x) = getValueNorm(pca.x(), -0.0241f, 0.0074f); + get(dnn::tau_dxy_pca_y) = getValueNorm(pca.y(), 0.0675f, 0.0128f); + get(dnn::tau_dxy_pca_z) = getValueNorm(pca.z(), 0.7973f, 3.456f); } else { get(dnn::tau_dxy_pca_x) = 0; get(dnn::tau_dxy_pca_y) = 0; get(dnn::tau_dxy_pca_z) = 0; } - const bool tau_dxy_valid = - std::isnormal(tau.dxy()) && tau.dxy() > -10 && std::isnormal(tau.dxy_error()) && tau.dxy_error() > 0; + isAbove(tau_funcs.getdxy(tau, tau_index), -10) && isAbove(tau_funcs.getdxyError(tau, tau_index), 0); if (tau_dxy_valid) { get(dnn::tau_dxy_valid) = tau_dxy_valid; - get(dnn::tau_dxy) = getValueNorm(tau.dxy(), 0.0018f, 0.0085f); - get(dnn::tau_dxy_sig) = getValueNorm(std::abs(tau.dxy()) / tau.dxy_error(), 2.26f, 4.191f); + get(dnn::tau_dxy) = getValueNorm(tau_funcs.getdxy(tau, tau_index), 0.0018f, 0.0085f); + get(dnn::tau_dxy_sig) = getValueNorm( + std::abs(tau_funcs.getdxy(tau, tau_index)) / tau_funcs.getdxyError(tau, tau_index), 2.26f, 4.191f); } const bool tau_ip3d_valid = - std::isnormal(tau.ip3d()) && tau.ip3d() > -10 && std::isnormal(tau.ip3d_error()) && tau.ip3d_error() > 0; + isAbove(tau_funcs.getip3d(tau, tau_index), -10) && isAbove(tau_funcs.getip3dError(tau, tau_index), 0); if (tau_ip3d_valid) { get(dnn::tau_ip3d_valid) = tau_ip3d_valid; - get(dnn::tau_ip3d) = getValueNorm(tau.ip3d(), 0.0026f, 0.0114f); - get(dnn::tau_ip3d_sig) = getValueNorm(std::abs(tau.ip3d()) / tau.ip3d_error(), 2.928f, 4.466f); + get(dnn::tau_ip3d) = getValueNorm(tau_funcs.getip3d(tau, tau_index), 0.0026f, 0.0114f); + get(dnn::tau_ip3d_sig) = getValueNorm( + std::abs(tau_funcs.getip3d(tau, tau_index)) / tau_funcs.getip3dError(tau, tau_index), 2.928f, 4.466f); } if (leadChargedHadrCand) { - get(dnn::tau_dz) = getValueNorm(leadChargedHadrCand->dz(), 0.f, 0.0190f); - const bool tau_dz_sig_valid = leadChargedHadrCand->hasTrackDetails() && - std::isnormal(leadChargedHadrCand->dz()) && - std::isnormal(leadChargedHadrCand->dzError()) && leadChargedHadrCand->dzError() > 0; - get(dnn::tau_dz_sig_valid) = tau_dz_sig_valid; - const double dzError = leadChargedHadrCand->hasTrackDetails() ? leadChargedHadrCand->dzError() : default_value; - get(dnn::tau_dz_sig) = getValueNorm(std::abs(leadChargedHadrCand->dz()) / dzError, 4.717f, 11.78f); - } - get(dnn::tau_flightLength_x) = getValueNorm(tau.flightLength().x(), -0.0003f, 0.7362f); - get(dnn::tau_flightLength_y) = getValueNorm(tau.flightLength().y(), -0.0009f, 0.7354f); - get(dnn::tau_flightLength_z) = getValueNorm(tau.flightLength().z(), -0.0022f, 1.993f); - // get(dnn::tau_flightLength_sig) = getValueNorm(tau.flightLengthSig(), -4.78f, 9.573f); + const bool hasTrackDetails = candFunc::getHasTrackDetails(*leadChargedHadrCand); + const float tau_dz = (is_online_ && !hasTrackDetails) ? 0 : candFunc::getTauDz(*leadChargedHadrCand); + get(dnn::tau_dz) = getValueNorm(tau_dz, 0.f, 0.0190f); + get(dnn::tau_dz_sig_valid) = candFunc::getTauDZSigValid(*leadChargedHadrCand); + const double dzError = hasTrackDetails ? leadChargedHadrCand->dzError() : -999.; + get(dnn::tau_dz_sig) = getValueNorm(std::abs(tau_dz) / dzError, 4.717f, 11.78f); + } + get(dnn::tau_flightLength_x) = getValueNorm(tau_funcs.getFlightLength(tau, tau_index).x(), -0.0003f, 0.7362f); + get(dnn::tau_flightLength_y) = getValueNorm(tau_funcs.getFlightLength(tau, tau_index).y(), -0.0009f, 0.7354f); + get(dnn::tau_flightLength_z) = getValueNorm(tau_funcs.getFlightLength(tau, tau_index).z(), -0.0022f, 1.993f); get(dnn::tau_flightLength_sig) = 0.55756444; //This value is set due to a bug in the training get(dnn::tau_pt_weighted_deta_strip) = getValueLinear(reco::tau::pt_weighted_deta_strip(tau, tau.decayMode()), 0, 1, true); @@ -1249,31 +1913,34 @@ class DeepTauId : public deep_tau::DeepTauBase { get(dnn::tau_pt_weighted_dr_signal) = getValueNorm(reco::tau::pt_weighted_dr_signal(tau, tau.decayMode()), 0.0052f, 0.01433f); get(dnn::tau_pt_weighted_dr_iso) = getValueLinear(reco::tau::pt_weighted_dr_iso(tau, tau.decayMode()), 0, 1, true); - get(dnn::tau_leadingTrackNormChi2) = getValueNorm(tau.leadingTrackNormChi2(), 1.538f, 4.401f); + get(dnn::tau_leadingTrackNormChi2) = getValueNorm(tau_funcs.getLeadingTrackNormChi2(tau), 1.538f, 4.401f); const auto eratio = reco::tau::eratio(tau); const bool tau_e_ratio_valid = std::isnormal(eratio) && eratio > 0.f; get(dnn::tau_e_ratio_valid) = tau_e_ratio_valid; get(dnn::tau_e_ratio) = tau_e_ratio_valid ? getValueLinear(eratio, 0, 1, true) : 0.f; - const double gj_angle_diff = calculateGottfriedJacksonAngleDifference(tau); + const double gj_angle_diff = calculateGottfriedJacksonAngleDifference(tau, tau_index, tau_funcs); const bool tau_gj_angle_diff_valid = (std::isnormal(gj_angle_diff) || gj_angle_diff == 0) && gj_angle_diff >= 0; get(dnn::tau_gj_angle_diff_valid) = tau_gj_angle_diff_valid; get(dnn::tau_gj_angle_diff) = tau_gj_angle_diff_valid ? getValueLinear(gj_angle_diff, 0, pi, true) : 0; get(dnn::tau_n_photons) = getValueNorm(reco::tau::n_photons_total(tau), 2.95f, 3.927f); - get(dnn::tau_emFraction) = getValueLinear(tau.emFraction_MVA(), -1, 1, false); + get(dnn::tau_emFraction) = getValueLinear(tau_funcs.getEmFraction(tau), -1, 1, false); + get(dnn::tau_inside_ecal_crack) = getValue(isInEcalCrack(tau.p4().eta())); get(dnn::leadChargedCand_etaAtEcalEntrance_minus_tau_eta) = - getValueNorm(tau.etaAtEcalEntranceLeadChargedCand() - tau.p4().eta(), 0.0042f, 0.0323f); - - checkInputs(inputs, "tau_block", dnn::NumberOfInputs); + getValueNorm(tau_funcs.getEtaAtEcalEntrance(tau) - tau.p4().eta(), 0.0042f, 0.0323f); } + template void createEgammaBlockInputs(unsigned idx, - const TauType& tau, + const TauCastType& tau, + const size_t tau_index, + const edm::RefToBase tau_ref, const reco::Vertex& pv, double rho, - const pat::ElectronCollection& electrons, - const pat::PackedCandidateCollection& pfCands, + const std::vector* electrons, + const edm::View& pfCands, const Cell& cell_map, + TauFunc tau_funcs, bool is_inner) { namespace dnn = dnn_inputs_2017_v2::EgammaBlockInputs; @@ -1293,6 +1960,7 @@ class DeepTauId : public deep_tau::DeepTauBase { } if (valid_index_pf_ele) { size_t index_pf_ele = cell_map.at(CellObjectType::PfCand_electron); + const auto& ele_cand = dynamic_cast(pfCands.at(index_pf_ele)); get(dnn::pfCand_ele_valid) = valid_index_pf_ele; get(dnn::pfCand_ele_rel_pt) = getValueNorm(pfCands.at(index_pf_ele).polarP4().pt() / tau.polarP4().pt(), @@ -1307,12 +1975,12 @@ class DeepTauId : public deep_tau::DeepTauBase { is_inner ? 0.1f : 0.5f, false); get(dnn::pfCand_ele_pvAssociationQuality) = - getValueLinear(pfCands.at(index_pf_ele).pvAssociationQuality(), 0, 7, true); - get(dnn::pfCand_ele_puppiWeight) = getValue(pfCands.at(index_pf_ele).puppiWeight()); - get(dnn::pfCand_ele_charge) = getValue(pfCands.at(index_pf_ele).charge()); - get(dnn::pfCand_ele_lostInnerHits) = getValue(pfCands.at(index_pf_ele).lostInnerHits()); - get(dnn::pfCand_ele_numberOfPixelHits) = - getValueLinear(pfCands.at(index_pf_ele).numberOfPixelHits(), 0, 10, true); + getValueLinear(candFunc::getPvAssocationQuality(ele_cand), 0, 7, true); + get(dnn::pfCand_ele_puppiWeight) = is_inner ? getValue(candFunc::getPuppiWeight(ele_cand, 0.9906834f)) + : getValue(candFunc::getPuppiWeight(ele_cand, 0.9669586f)); + get(dnn::pfCand_ele_charge) = getValue(ele_cand.charge()); + get(dnn::pfCand_ele_lostInnerHits) = getValue(candFunc::getLostInnerHits(ele_cand, 0)); + get(dnn::pfCand_ele_numberOfPixelHits) = getValueLinear(candFunc::getNumberOfPixelHits(ele_cand, 0), 0, 10, true); get(dnn::pfCand_ele_vertex_dx) = getValueNorm(pfCands.at(index_pf_ele).vertex().x() - pv.position().x(), 0.f, 0.1221f); get(dnn::pfCand_ele_vertex_dy) = @@ -1320,30 +1988,36 @@ class DeepTauId : public deep_tau::DeepTauBase { get(dnn::pfCand_ele_vertex_dz) = getValueNorm(pfCands.at(index_pf_ele).vertex().z() - pv.position().z(), 0.001f, 1.024f); get(dnn::pfCand_ele_vertex_dx_tauFL) = getValueNorm( - pfCands.at(index_pf_ele).vertex().x() - pv.position().x() - tau.flightLength().x(), 0.f, 0.3411f); + pfCands.at(index_pf_ele).vertex().x() - pv.position().x() - tau_funcs.getFlightLength(tau, tau_index).x(), + 0.f, + 0.3411f); get(dnn::pfCand_ele_vertex_dy_tauFL) = getValueNorm( - pfCands.at(index_pf_ele).vertex().y() - pv.position().y() - tau.flightLength().y(), 0.0003f, 0.3385f); - get(dnn::pfCand_ele_vertex_dz_tauFL) = - getValueNorm(pfCands.at(index_pf_ele).vertex().z() - pv.position().z() - tau.flightLength().z(), 0.f, 1.307f); - - const bool hasTrackDetails = pfCands.at(index_pf_ele).hasTrackDetails(); + pfCands.at(index_pf_ele).vertex().y() - pv.position().y() - tau_funcs.getFlightLength(tau, tau_index).y(), + 0.0003f, + 0.3385f); + get(dnn::pfCand_ele_vertex_dz_tauFL) = getValueNorm( + pfCands.at(index_pf_ele).vertex().z() - pv.position().z() - tau_funcs.getFlightLength(tau, tau_index).z(), + 0.f, + 1.307f); + + const bool hasTrackDetails = candFunc::getHasTrackDetails(ele_cand); if (hasTrackDetails) { get(dnn::pfCand_ele_hasTrackDetails) = hasTrackDetails; - get(dnn::pfCand_ele_dxy) = getValueNorm(pfCands.at(index_pf_ele).dxy(), 0.f, 0.171f); + get(dnn::pfCand_ele_dxy) = getValueNorm(candFunc::getTauDxy(ele_cand), 0.f, 0.171f); get(dnn::pfCand_ele_dxy_sig) = - getValueNorm(std::abs(pfCands.at(index_pf_ele).dxy()) / pfCands.at(index_pf_ele).dxyError(), 1.634f, 6.45f); - get(dnn::pfCand_ele_dz) = getValueNorm(pfCands.at(index_pf_ele).dz(), 0.001f, 1.02f); + getValueNorm(std::abs(candFunc::getTauDxy(ele_cand)) / pfCands.at(index_pf_ele).dxyError(), 1.634f, 6.45f); + get(dnn::pfCand_ele_dz) = getValueNorm(candFunc::getTauDz(ele_cand), 0.001f, 1.02f); get(dnn::pfCand_ele_dz_sig) = - getValueNorm(std::abs(pfCands.at(index_pf_ele).dz()) / pfCands.at(index_pf_ele).dzError(), 24.56f, 210.4f); - get(dnn::pfCand_ele_track_chi2_ndof) = - getValueNorm(pfCands.at(index_pf_ele).pseudoTrack().chi2() / pfCands.at(index_pf_ele).pseudoTrack().ndof(), - 2.272f, - 8.439f); - get(dnn::pfCand_ele_track_ndof) = getValueNorm(pfCands.at(index_pf_ele).pseudoTrack().ndof(), 15.18f, 3.203f); + getValueNorm(std::abs(candFunc::getTauDz(ele_cand)) / ele_cand.dzError(), 24.56f, 210.4f); + get(dnn::pfCand_ele_track_chi2_ndof) = getValueNorm( + candFunc::getPseudoTrack(ele_cand).chi2() / candFunc::getPseudoTrack(ele_cand).ndof(), 2.272f, 8.439f); + get(dnn::pfCand_ele_track_ndof) = getValueNorm(candFunc::getPseudoTrack(ele_cand).ndof(), 15.18f, 3.203f); } } if (valid_index_pf_gamma) { size_t index_pf_gamma = cell_map.at(CellObjectType::PfCand_gamma); + const auto& gamma_cand = dynamic_cast(pfCands.at(index_pf_gamma)); + get(dnn::pfCand_gamma_valid) = valid_index_pf_gamma; get(dnn::pfCand_gamma_rel_pt) = getValueNorm(pfCands.at(index_pf_gamma).polarP4().pt() / tau.polarP4().pt(), is_inner ? 0.6048f : 0.02576f, @@ -1357,13 +2031,16 @@ class DeepTauId : public deep_tau::DeepTauBase { is_inner ? 0.1f : 0.5f, false); get(dnn::pfCand_gamma_pvAssociationQuality) = - getValueLinear(pfCands.at(index_pf_gamma).pvAssociationQuality(), 0, 7, true); - get(dnn::pfCand_gamma_fromPV) = getValueLinear(pfCands.at(index_pf_gamma).fromPV(), 0, 3, true); - get(dnn::pfCand_gamma_puppiWeight) = getValue(pfCands.at(index_pf_gamma).puppiWeight()); - get(dnn::pfCand_gamma_puppiWeightNoLep) = getValue(pfCands.at(index_pf_gamma).puppiWeightNoLep()); - get(dnn::pfCand_gamma_lostInnerHits) = getValue(pfCands.at(index_pf_gamma).lostInnerHits()); + getValueLinear(candFunc::getPvAssocationQuality(gamma_cand), 0, 7, true); + get(dnn::pfCand_gamma_fromPV) = getValueLinear(candFunc::getFromPV(gamma_cand), 0, 3, true); + get(dnn::pfCand_gamma_puppiWeight) = is_inner ? getValue(candFunc::getPuppiWeight(gamma_cand, 0.9084110f)) + : getValue(candFunc::getPuppiWeight(gamma_cand, 0.4211567f)); + get(dnn::pfCand_gamma_puppiWeightNoLep) = is_inner + ? getValue(candFunc::getPuppiWeightNoLep(gamma_cand, 0.8857716f)) + : getValue(candFunc::getPuppiWeightNoLep(gamma_cand, 0.3822604f)); + get(dnn::pfCand_gamma_lostInnerHits) = getValue(candFunc::getLostInnerHits(gamma_cand, 0)); get(dnn::pfCand_gamma_numberOfPixelHits) = - getValueLinear(pfCands.at(index_pf_gamma).numberOfPixelHits(), 0, 7, true); + getValueLinear(candFunc::getNumberOfPixelHits(gamma_cand, 0), 0, 7, true); get(dnn::pfCand_gamma_vertex_dx) = getValueNorm(pfCands.at(index_pf_gamma).vertex().x() - pv.position().x(), 0.f, 0.0067f); get(dnn::pfCand_gamma_vertex_dy) = @@ -1371,30 +2048,35 @@ class DeepTauId : public deep_tau::DeepTauBase { get(dnn::pfCand_gamma_vertex_dz) = getValueNorm(pfCands.at(index_pf_gamma).vertex().z() - pv.position().z(), 0.f, 0.0578f); get(dnn::pfCand_gamma_vertex_dx_tauFL) = getValueNorm( - pfCands.at(index_pf_gamma).vertex().x() - pv.position().x() - tau.flightLength().x(), 0.001f, 0.9565f); + pfCands.at(index_pf_gamma).vertex().x() - pv.position().x() - tau_funcs.getFlightLength(tau, tau_index).x(), + 0.001f, + 0.9565f); get(dnn::pfCand_gamma_vertex_dy_tauFL) = getValueNorm( - pfCands.at(index_pf_gamma).vertex().y() - pv.position().y() - tau.flightLength().y(), 0.0008f, 0.9592f); + pfCands.at(index_pf_gamma).vertex().y() - pv.position().y() - tau_funcs.getFlightLength(tau, tau_index).y(), + 0.0008f, + 0.9592f); get(dnn::pfCand_gamma_vertex_dz_tauFL) = getValueNorm( - pfCands.at(index_pf_gamma).vertex().z() - pv.position().z() - tau.flightLength().z(), 0.0038f, 2.154f); - - const bool hasTrackDetails = pfCands.at(index_pf_gamma).hasTrackDetails(); + pfCands.at(index_pf_gamma).vertex().z() - pv.position().z() - tau_funcs.getFlightLength(tau, tau_index).z(), + 0.0038f, + 2.154f); + const bool hasTrackDetails = candFunc::getHasTrackDetails(gamma_cand); if (hasTrackDetails) { get(dnn::pfCand_gamma_hasTrackDetails) = hasTrackDetails; - get(dnn::pfCand_gamma_dxy) = getValueNorm(pfCands.at(index_pf_gamma).dxy(), 0.0004f, 0.882f); - get(dnn::pfCand_gamma_dxy_sig) = getValueNorm( - std::abs(pfCands.at(index_pf_gamma).dxy()) / pfCands.at(index_pf_gamma).dxyError(), 4.271f, 63.78f); - get(dnn::pfCand_gamma_dz) = getValueNorm(pfCands.at(index_pf_gamma).dz(), 0.0071f, 5.285f); - get(dnn::pfCand_gamma_dz_sig) = getValueNorm( - std::abs(pfCands.at(index_pf_gamma).dz()) / pfCands.at(index_pf_gamma).dzError(), 162.1f, 622.4f); - get(dnn::pfCand_gamma_track_chi2_ndof) = pfCands.at(index_pf_gamma).pseudoTrack().ndof() > 0 - ? getValueNorm(pfCands.at(index_pf_gamma).pseudoTrack().chi2() / - pfCands.at(index_pf_gamma).pseudoTrack().ndof(), + get(dnn::pfCand_gamma_dxy) = getValueNorm(candFunc::getTauDxy(gamma_cand), 0.0004f, 0.882f); + get(dnn::pfCand_gamma_dxy_sig) = + getValueNorm(std::abs(candFunc::getTauDxy(gamma_cand)) / gamma_cand.dxyError(), 4.271f, 63.78f); + get(dnn::pfCand_gamma_dz) = getValueNorm(candFunc::getTauDz(gamma_cand), 0.0071f, 5.285f); + get(dnn::pfCand_gamma_dz_sig) = + getValueNorm(std::abs(candFunc::getTauDz(gamma_cand)) / gamma_cand.dzError(), 162.1f, 622.4f); + get(dnn::pfCand_gamma_track_chi2_ndof) = candFunc::getPseudoTrack(gamma_cand).ndof() > 0 + ? getValueNorm(candFunc::getPseudoTrack(gamma_cand).chi2() / + candFunc::getPseudoTrack(gamma_cand).ndof(), 4.268f, 15.47f) : 0; get(dnn::pfCand_gamma_track_ndof) = - pfCands.at(index_pf_gamma).pseudoTrack().ndof() > 0 - ? getValueNorm(pfCands.at(index_pf_gamma).pseudoTrack().ndof(), 12.25f, 4.774f) + candFunc::getPseudoTrack(gamma_cand).ndof() > 0 + ? getValueNorm(candFunc::getPseudoTrack(gamma_cand).ndof(), 12.25f, 4.774f) : 0; } } @@ -1402,14 +2084,14 @@ class DeepTauId : public deep_tau::DeepTauBase { size_t index_ele = cell_map.at(CellObjectType::Electron); get(dnn::ele_valid) = valid_index_ele; - get(dnn::ele_rel_pt) = getValueNorm(electrons.at(index_ele).polarP4().pt() / tau.polarP4().pt(), + get(dnn::ele_rel_pt) = getValueNorm(electrons->at(index_ele).polarP4().pt() / tau.polarP4().pt(), is_inner ? 1.067f : 0.5111f, is_inner ? 1.521f : 2.765f); - get(dnn::ele_deta) = getValueLinear(electrons.at(index_ele).polarP4().eta() - tau.polarP4().eta(), + get(dnn::ele_deta) = getValueLinear(electrons->at(index_ele).polarP4().eta() - tau.polarP4().eta(), is_inner ? -0.1f : -0.5f, is_inner ? 0.1f : 0.5f, false); - get(dnn::ele_dphi) = getValueLinear(dPhi(tau.polarP4(), electrons.at(index_ele).polarP4()), + get(dnn::ele_dphi) = getValueLinear(dPhi(tau.polarP4(), electrons->at(index_ele).polarP4()), is_inner ? -0.1f : -0.5f, is_inner ? 0.1f : 0.5f, false); @@ -1417,64 +2099,63 @@ class DeepTauId : public deep_tau::DeepTauBase { float cc_ele_energy, cc_gamma_energy; int cc_n_gamma; const bool cc_valid = - calculateElectronClusterVarsV2(electrons.at(index_ele), cc_ele_energy, cc_gamma_energy, cc_n_gamma); + calculateElectronClusterVarsV2(electrons->at(index_ele), cc_ele_energy, cc_gamma_energy, cc_n_gamma); if (cc_valid) { get(dnn::ele_cc_valid) = cc_valid; get(dnn::ele_cc_ele_rel_energy) = - getValueNorm(cc_ele_energy / electrons.at(index_ele).polarP4().pt(), 1.729f, 1.644f); + getValueNorm(cc_ele_energy / electrons->at(index_ele).polarP4().pt(), 1.729f, 1.644f); get(dnn::ele_cc_gamma_rel_energy) = getValueNorm(cc_gamma_energy / cc_ele_energy, 0.1439f, 0.3284f); get(dnn::ele_cc_n_gamma) = getValueNorm(cc_n_gamma, 1.794f, 2.079f); } get(dnn::ele_rel_trackMomentumAtVtx) = getValueNorm( - electrons.at(index_ele).trackMomentumAtVtx().R() / electrons.at(index_ele).polarP4().pt(), 1.531f, 1.424f); + electrons->at(index_ele).trackMomentumAtVtx().R() / electrons->at(index_ele).polarP4().pt(), 1.531f, 1.424f); get(dnn::ele_rel_trackMomentumAtCalo) = getValueNorm( - electrons.at(index_ele).trackMomentumAtCalo().R() / electrons.at(index_ele).polarP4().pt(), 1.531f, 1.424f); + electrons->at(index_ele).trackMomentumAtCalo().R() / electrons->at(index_ele).polarP4().pt(), 1.531f, 1.424f); get(dnn::ele_rel_trackMomentumOut) = getValueNorm( - electrons.at(index_ele).trackMomentumOut().R() / electrons.at(index_ele).polarP4().pt(), 0.7735f, 0.935f); + electrons->at(index_ele).trackMomentumOut().R() / electrons->at(index_ele).polarP4().pt(), 0.7735f, 0.935f); get(dnn::ele_rel_trackMomentumAtEleClus) = - getValueNorm(electrons.at(index_ele).trackMomentumAtEleClus().R() / electrons.at(index_ele).polarP4().pt(), + getValueNorm(electrons->at(index_ele).trackMomentumAtEleClus().R() / electrons->at(index_ele).polarP4().pt(), 0.7735f, 0.935f); get(dnn::ele_rel_trackMomentumAtVtxWithConstraint) = getValueNorm( - electrons.at(index_ele).trackMomentumAtVtxWithConstraint().R() / electrons.at(index_ele).polarP4().pt(), + electrons->at(index_ele).trackMomentumAtVtxWithConstraint().R() / electrons->at(index_ele).polarP4().pt(), 1.625f, 1.581f); get(dnn::ele_rel_ecalEnergy) = - getValueNorm(electrons.at(index_ele).ecalEnergy() / electrons.at(index_ele).polarP4().pt(), 1.993f, 1.308f); + getValueNorm(electrons->at(index_ele).ecalEnergy() / electrons->at(index_ele).polarP4().pt(), 1.993f, 1.308f); get(dnn::ele_ecalEnergy_sig) = getValueNorm( - electrons.at(index_ele).ecalEnergy() / electrons.at(index_ele).ecalEnergyError(), 70.25f, 58.16f); - get(dnn::ele_eSuperClusterOverP) = getValueNorm(electrons.at(index_ele).eSuperClusterOverP(), 2.432f, 15.13f); - get(dnn::ele_eSeedClusterOverP) = getValueNorm(electrons.at(index_ele).eSeedClusterOverP(), 2.034f, 13.96f); - get(dnn::ele_eSeedClusterOverPout) = getValueNorm(electrons.at(index_ele).eSeedClusterOverPout(), 6.64f, 36.8f); - get(dnn::ele_eEleClusterOverPout) = getValueNorm(electrons.at(index_ele).eEleClusterOverPout(), 4.183f, 20.63f); + electrons->at(index_ele).ecalEnergy() / electrons->at(index_ele).ecalEnergyError(), 70.25f, 58.16f); + get(dnn::ele_eSuperClusterOverP) = getValueNorm(electrons->at(index_ele).eSuperClusterOverP(), 2.432f, 15.13f); + get(dnn::ele_eSeedClusterOverP) = getValueNorm(electrons->at(index_ele).eSeedClusterOverP(), 2.034f, 13.96f); + get(dnn::ele_eSeedClusterOverPout) = getValueNorm(electrons->at(index_ele).eSeedClusterOverPout(), 6.64f, 36.8f); + get(dnn::ele_eEleClusterOverPout) = getValueNorm(electrons->at(index_ele).eEleClusterOverPout(), 4.183f, 20.63f); get(dnn::ele_deltaEtaSuperClusterTrackAtVtx) = - getValueNorm(electrons.at(index_ele).deltaEtaSuperClusterTrackAtVtx(), 0.f, 0.0363f); + getValueNorm(electrons->at(index_ele).deltaEtaSuperClusterTrackAtVtx(), 0.f, 0.0363f); get(dnn::ele_deltaEtaSeedClusterTrackAtCalo) = - getValueNorm(electrons.at(index_ele).deltaEtaSeedClusterTrackAtCalo(), -0.0001f, 0.0512f); + getValueNorm(electrons->at(index_ele).deltaEtaSeedClusterTrackAtCalo(), -0.0001f, 0.0512f); get(dnn::ele_deltaEtaEleClusterTrackAtCalo) = - getValueNorm(electrons.at(index_ele).deltaEtaEleClusterTrackAtCalo(), -0.0001f, 0.0541f); + getValueNorm(electrons->at(index_ele).deltaEtaEleClusterTrackAtCalo(), -0.0001f, 0.0541f); get(dnn::ele_deltaPhiEleClusterTrackAtCalo) = - getValueNorm(electrons.at(index_ele).deltaPhiEleClusterTrackAtCalo(), 0.0002f, 0.0553f); + getValueNorm(electrons->at(index_ele).deltaPhiEleClusterTrackAtCalo(), 0.0002f, 0.0553f); get(dnn::ele_deltaPhiSuperClusterTrackAtVtx) = - getValueNorm(electrons.at(index_ele).deltaPhiSuperClusterTrackAtVtx(), 0.0001f, 0.0523f); + getValueNorm(electrons->at(index_ele).deltaPhiSuperClusterTrackAtVtx(), 0.0001f, 0.0523f); get(dnn::ele_deltaPhiSeedClusterTrackAtCalo) = - getValueNorm(electrons.at(index_ele).deltaPhiSeedClusterTrackAtCalo(), 0.0004f, 0.0777f); - get(dnn::ele_mvaInput_earlyBrem) = getValue(electrons.at(index_ele).mvaInput().earlyBrem); - get(dnn::ele_mvaInput_lateBrem) = getValue(electrons.at(index_ele).mvaInput().lateBrem); + getValueNorm(electrons->at(index_ele).deltaPhiSeedClusterTrackAtCalo(), 0.0004f, 0.0777f); + get(dnn::ele_mvaInput_earlyBrem) = getValue(electrons->at(index_ele).mvaInput().earlyBrem); + get(dnn::ele_mvaInput_lateBrem) = getValue(electrons->at(index_ele).mvaInput().lateBrem); get(dnn::ele_mvaInput_sigmaEtaEta) = - getValueNorm(electrons.at(index_ele).mvaInput().sigmaEtaEta, 0.0008f, 0.0052f); - get(dnn::ele_mvaInput_hadEnergy) = getValueNorm(electrons.at(index_ele).mvaInput().hadEnergy, 14.04f, 69.48f); - get(dnn::ele_mvaInput_deltaEta) = getValueNorm(electrons.at(index_ele).mvaInput().deltaEta, 0.0099f, 0.0851f); - - const auto& gsfTrack = electrons.at(index_ele).gsfTrack(); + getValueNorm(electrons->at(index_ele).mvaInput().sigmaEtaEta, 0.0008f, 0.0052f); + get(dnn::ele_mvaInput_hadEnergy) = getValueNorm(electrons->at(index_ele).mvaInput().hadEnergy, 14.04f, 69.48f); + get(dnn::ele_mvaInput_deltaEta) = getValueNorm(electrons->at(index_ele).mvaInput().deltaEta, 0.0099f, 0.0851f); + const auto& gsfTrack = electrons->at(index_ele).gsfTrack(); if (gsfTrack.isNonnull()) { get(dnn::ele_gsfTrack_normalizedChi2) = getValueNorm(gsfTrack->normalizedChi2(), 3.049f, 10.39f); get(dnn::ele_gsfTrack_numberOfValidHits) = getValueNorm(gsfTrack->numberOfValidHits(), 16.52f, 2.806f); get(dnn::ele_rel_gsfTrack_pt) = - getValueNorm(gsfTrack->pt() / electrons.at(index_ele).polarP4().pt(), 1.355f, 16.81f); + getValueNorm(gsfTrack->pt() / electrons->at(index_ele).polarP4().pt(), 1.355f, 16.81f); get(dnn::ele_gsfTrack_pt_sig) = getValueNorm(gsfTrack->pt() / gsfTrack->ptError(), 5.046f, 3.119f); } - const auto& closestCtfTrack = electrons.at(index_ele).closestCtfTrackRef(); + const auto& closestCtfTrack = electrons->at(index_ele).closestCtfTrackRef(); const bool has_closestCtfTrack = closestCtfTrack.isNonnull(); if (has_closestCtfTrack) { get(dnn::ele_has_closestCtfTrack) = has_closestCtfTrack; @@ -1483,16 +2164,19 @@ class DeepTauId : public deep_tau::DeepTauBase { getValueNorm(closestCtfTrack->numberOfValidHits(), 15.16f, 5.26f); } } - checkInputs(inputs, is_inner ? "egamma_inner_block" : "egamma_outer_block", dnn::NumberOfInputs); } + template void createMuonBlockInputs(unsigned idx, - const TauType& tau, + const TauCastType& tau, + const size_t tau_index, + const edm::RefToBase tau_ref, const reco::Vertex& pv, double rho, - const pat::MuonCollection& muons, - const pat::PackedCandidateCollection& pfCands, + const std::vector* muons, + const edm::View& pfCands, const Cell& cell_map, + TauFunc tau_funcs, bool is_inner) { namespace dnn = dnn_inputs_2017_v2::MuonBlockInputs; @@ -1511,6 +2195,7 @@ class DeepTauId : public deep_tau::DeepTauBase { } if (valid_index_pf_muon) { size_t index_pf_muon = cell_map.at(CellObjectType::PfCand_muon); + const auto& muon_cand = dynamic_cast(pfCands.at(index_pf_muon)); get(dnn::pfCand_muon_valid) = valid_index_pf_muon; get(dnn::pfCand_muon_rel_pt) = getValueNorm(pfCands.at(index_pf_muon).polarP4().pt() / tau.polarP4().pt(), @@ -1525,13 +2210,14 @@ class DeepTauId : public deep_tau::DeepTauBase { is_inner ? 0.1f : 0.5f, false); get(dnn::pfCand_muon_pvAssociationQuality) = - getValueLinear(pfCands.at(index_pf_muon).pvAssociationQuality(), 0, 7, true); - get(dnn::pfCand_muon_fromPV) = getValueLinear(pfCands.at(index_pf_muon).fromPV(), 0, 3, true); - get(dnn::pfCand_muon_puppiWeight) = getValue(pfCands.at(index_pf_muon).puppiWeight()); - get(dnn::pfCand_muon_charge) = getValue(pfCands.at(index_pf_muon).charge()); - get(dnn::pfCand_muon_lostInnerHits) = getValue(pfCands.at(index_pf_muon).lostInnerHits()); + getValueLinear(candFunc::getPvAssocationQuality(muon_cand), 0, 7, true); + get(dnn::pfCand_muon_fromPV) = getValueLinear(candFunc::getFromPV(muon_cand), 0, 3, true); + get(dnn::pfCand_muon_puppiWeight) = is_inner ? getValue(candFunc::getPuppiWeight(muon_cand, 0.9786588f)) + : getValue(candFunc::getPuppiWeight(muon_cand, 0.8132477f)); + get(dnn::pfCand_muon_charge) = getValue(muon_cand.charge()); + get(dnn::pfCand_muon_lostInnerHits) = getValue(candFunc::getLostInnerHits(muon_cand, 0)); get(dnn::pfCand_muon_numberOfPixelHits) = - getValueLinear(pfCands.at(index_pf_muon).numberOfPixelHits(), 0, 11, true); + getValueLinear(candFunc::getNumberOfPixelHits(muon_cand, 0), 0, 11, true); get(dnn::pfCand_muon_vertex_dx) = getValueNorm(pfCands.at(index_pf_muon).vertex().x() - pv.position().x(), -0.0007f, 0.6869f); get(dnn::pfCand_muon_vertex_dy) = @@ -1539,66 +2225,72 @@ class DeepTauId : public deep_tau::DeepTauBase { get(dnn::pfCand_muon_vertex_dz) = getValueNorm(pfCands.at(index_pf_muon).vertex().z() - pv.position().z(), -0.0117f, 4.097f); get(dnn::pfCand_muon_vertex_dx_tauFL) = getValueNorm( - pfCands.at(index_pf_muon).vertex().x() - pv.position().x() - tau.flightLength().x(), -0.0001f, 0.8642f); + pfCands.at(index_pf_muon).vertex().x() - pv.position().x() - tau_funcs.getFlightLength(tau, tau_index).x(), + -0.0001f, + 0.8642f); get(dnn::pfCand_muon_vertex_dy_tauFL) = getValueNorm( - pfCands.at(index_pf_muon).vertex().y() - pv.position().y() - tau.flightLength().y(), 0.0004f, 0.8561f); + pfCands.at(index_pf_muon).vertex().y() - pv.position().y() - tau_funcs.getFlightLength(tau, tau_index).y(), + 0.0004f, + 0.8561f); get(dnn::pfCand_muon_vertex_dz_tauFL) = getValueNorm( - pfCands.at(index_pf_muon).vertex().z() - pv.position().z() - tau.flightLength().z(), -0.0118f, 4.405f); + pfCands.at(index_pf_muon).vertex().z() - pv.position().z() - tau_funcs.getFlightLength(tau, tau_index).z(), + -0.0118f, + 4.405f); - const bool hasTrackDetails = pfCands.at(index_pf_muon).hasTrackDetails(); + const bool hasTrackDetails = candFunc::getHasTrackDetails(muon_cand); if (hasTrackDetails) { get(dnn::pfCand_muon_hasTrackDetails) = hasTrackDetails; - get(dnn::pfCand_muon_dxy) = getValueNorm(pfCands.at(index_pf_muon).dxy(), -0.0045f, 0.9655f); - get(dnn::pfCand_muon_dxy_sig) = getValueNorm( - std::abs(pfCands.at(index_pf_muon).dxy()) / pfCands.at(index_pf_muon).dxyError(), 4.575f, 42.36f); - get(dnn::pfCand_muon_dz) = getValueNorm(pfCands.at(index_pf_muon).dz(), -0.0117f, 4.097f); - get(dnn::pfCand_muon_dz_sig) = getValueNorm( - std::abs(pfCands.at(index_pf_muon).dz()) / pfCands.at(index_pf_muon).dzError(), 80.37f, 343.3f); + get(dnn::pfCand_muon_dxy) = getValueNorm(candFunc::getTauDxy(muon_cand), -0.0045f, 0.9655f); + get(dnn::pfCand_muon_dxy_sig) = + getValueNorm(std::abs(candFunc::getTauDxy(muon_cand)) / muon_cand.dxyError(), 4.575f, 42.36f); + get(dnn::pfCand_muon_dz) = getValueNorm(candFunc::getTauDz(muon_cand), -0.0117f, 4.097f); + get(dnn::pfCand_muon_dz_sig) = + getValueNorm(std::abs(candFunc::getTauDz(muon_cand)) / muon_cand.dzError(), 80.37f, 343.3f); get(dnn::pfCand_muon_track_chi2_ndof) = getValueNorm( - pfCands.at(index_pf_muon).pseudoTrack().chi2() / pfCands.at(index_pf_muon).pseudoTrack().ndof(), - 0.69f, - 1.711f); - get(dnn::pfCand_muon_track_ndof) = getValueNorm(pfCands.at(index_pf_muon).pseudoTrack().ndof(), 17.5f, 5.11f); + candFunc::getPseudoTrack(muon_cand).chi2() / candFunc::getPseudoTrack(muon_cand).ndof(), 0.69f, 1.711f); + get(dnn::pfCand_muon_track_ndof) = getValueNorm(candFunc::getPseudoTrack(muon_cand).ndof(), 17.5f, 5.11f); } } if (valid_index_muon) { size_t index_muon = cell_map.at(CellObjectType::Muon); get(dnn::muon_valid) = valid_index_muon; - get(dnn::muon_rel_pt) = getValueNorm(muons.at(index_muon).polarP4().pt() / tau.polarP4().pt(), + get(dnn::muon_rel_pt) = getValueNorm(muons->at(index_muon).polarP4().pt() / tau.polarP4().pt(), is_inner ? 0.7966f : 0.2678f, is_inner ? 3.402f : 3.592f); - get(dnn::muon_deta) = getValueLinear(muons.at(index_muon).polarP4().eta() - tau.polarP4().eta(), + get(dnn::muon_deta) = getValueLinear(muons->at(index_muon).polarP4().eta() - tau.polarP4().eta(), + is_inner ? -0.1f : -0.5f, + is_inner ? 0.1f : 0.5f, + false); + get(dnn::muon_dphi) = getValueLinear(dPhi(tau.polarP4(), muons->at(index_muon).polarP4()), is_inner ? -0.1f : -0.5f, is_inner ? 0.1f : 0.5f, false); - get(dnn::muon_dphi) = getValueLinear( - dPhi(tau.polarP4(), muons.at(index_muon).polarP4()), is_inner ? -0.1f : -0.5f, is_inner ? 0.1f : 0.5f, false); - get(dnn::muon_dxy) = getValueNorm(muons.at(index_muon).dB(pat::Muon::PV2D), 0.0019f, 1.039f); + get(dnn::muon_dxy) = getValueNorm(muons->at(index_muon).dB(pat::Muon::PV2D), 0.0019f, 1.039f); get(dnn::muon_dxy_sig) = - getValueNorm(std::abs(muons.at(index_muon).dB(pat::Muon::PV2D)) / muons.at(index_muon).edB(pat::Muon::PV2D), + getValueNorm(std::abs(muons->at(index_muon).dB(pat::Muon::PV2D)) / muons->at(index_muon).edB(pat::Muon::PV2D), 8.98f, 71.17f); const bool normalizedChi2_valid = - muons.at(index_muon).globalTrack().isNonnull() && muons.at(index_muon).normChi2() >= 0; + muons->at(index_muon).globalTrack().isNonnull() && muons->at(index_muon).normChi2() >= 0; if (normalizedChi2_valid) { get(dnn::muon_normalizedChi2_valid) = normalizedChi2_valid; - get(dnn::muon_normalizedChi2) = getValueNorm(muons.at(index_muon).normChi2(), 21.52f, 265.8f); - if (muons.at(index_muon).innerTrack().isNonnull()) - get(dnn::muon_numberOfValidHits) = getValueNorm(muons.at(index_muon).numberOfValidHits(), 21.84f, 10.59f); + get(dnn::muon_normalizedChi2) = getValueNorm(muons->at(index_muon).normChi2(), 21.52f, 265.8f); + if (muons->at(index_muon).innerTrack().isNonnull()) + get(dnn::muon_numberOfValidHits) = getValueNorm(muons->at(index_muon).numberOfValidHits(), 21.84f, 10.59f); } - get(dnn::muon_segmentCompatibility) = getValue(muons.at(index_muon).segmentCompatibility()); - get(dnn::muon_caloCompatibility) = getValue(muons.at(index_muon).caloCompatibility()); + get(dnn::muon_segmentCompatibility) = getValue(muons->at(index_muon).segmentCompatibility()); + get(dnn::muon_caloCompatibility) = getValue(muons->at(index_muon).caloCompatibility()); - const bool pfEcalEnergy_valid = muons.at(index_muon).pfEcalEnergy() >= 0; + const bool pfEcalEnergy_valid = muons->at(index_muon).pfEcalEnergy() >= 0; if (pfEcalEnergy_valid) { get(dnn::muon_pfEcalEnergy_valid) = pfEcalEnergy_valid; get(dnn::muon_rel_pfEcalEnergy) = - getValueNorm(muons.at(index_muon).pfEcalEnergy() / muons.at(index_muon).polarP4().pt(), 0.2273f, 0.4865f); + getValueNorm(muons->at(index_muon).pfEcalEnergy() / muons->at(index_muon).polarP4().pt(), 0.2273f, 0.4865f); } - MuonHitMatchV2 hit_match(muons.at(index_muon)); + MuonHitMatchV2 hit_match(muons->at(index_muon)); static const std::map> muonMatchHitVars = { {MuonSubdetId::DT, {dnn::muon_n_matches_DT_1, dnn::muon_n_hits_DT_1}}, {MuonSubdetId::CSC, {dnn::muon_n_matches_CSC_1, dnn::muon_n_hits_CSC_1}}, @@ -1623,15 +2315,18 @@ class DeepTauId : public deep_tau::DeepTauBase { } } } - checkInputs(inputs, is_inner ? "muon_inner_block" : "muon_outer_block", dnn::NumberOfInputs); } + template void createHadronsBlockInputs(unsigned idx, - const TauType& tau, + const TauCastType& tau, + const size_t tau_index, + const edm::RefToBase tau_ref, const reco::Vertex& pv, double rho, - const pat::PackedCandidateCollection& pfCands, + const edm::View& pfCands, const Cell& cell_map, + TauFunc tau_funcs, bool is_inner) { namespace dnn = dnn_inputs_2017_v2::HadronBlockInputs; @@ -1650,6 +2345,7 @@ class DeepTauId : public deep_tau::DeepTauBase { } if (valid_chH) { size_t index_chH = cell_map.at(CellObjectType::PfCand_chargedHadron); + const auto& chH_cand = dynamic_cast(pfCands.at(index_chH)); get(dnn::pfCand_chHad_valid) = valid_chH; get(dnn::pfCand_chHad_rel_pt) = getValueNorm(pfCands.at(index_chH).polarP4().pt() / tau.polarP4().pt(), @@ -1663,16 +2359,23 @@ class DeepTauId : public deep_tau::DeepTauBase { is_inner ? -0.1f : -0.5f, is_inner ? 0.1f : 0.5f, false); - get(dnn::pfCand_chHad_leadChargedHadrCand) = getValue( - &pfCands.at(index_chH) == dynamic_cast(tau.leadChargedHadrCand().get())); + get(dnn::pfCand_chHad_leadChargedHadrCand) = + getValue(&chH_cand == dynamic_cast(tau.leadChargedHadrCand().get())); get(dnn::pfCand_chHad_pvAssociationQuality) = - getValueLinear(pfCands.at(index_chH).pvAssociationQuality(), 0, 7, true); - get(dnn::pfCand_chHad_fromPV) = getValueLinear(pfCands.at(index_chH).fromPV(), 0, 3, true); - get(dnn::pfCand_chHad_puppiWeight) = getValue(pfCands.at(index_chH).puppiWeight()); - get(dnn::pfCand_chHad_puppiWeightNoLep) = getValue(pfCands.at(index_chH).puppiWeightNoLep()); - get(dnn::pfCand_chHad_charge) = getValue(pfCands.at(index_chH).charge()); - get(dnn::pfCand_chHad_lostInnerHits) = getValue(pfCands.at(index_chH).lostInnerHits()); - get(dnn::pfCand_chHad_numberOfPixelHits) = getValueLinear(pfCands.at(index_chH).numberOfPixelHits(), 0, 12, true); + getValueLinear(candFunc::getPvAssocationQuality(chH_cand), 0, 7, true); + get(dnn::pfCand_chHad_fromPV) = getValueLinear(candFunc::getFromPV(chH_cand), 0, 3, true); + const float default_chH_pw_inner = 0.7614090f; + const float default_chH_pw_outer = 0.1974930f; + get(dnn::pfCand_chHad_puppiWeight) = is_inner + ? getValue(candFunc::getPuppiWeight(chH_cand, default_chH_pw_inner)) + : getValue(candFunc::getPuppiWeight(chH_cand, default_chH_pw_outer)); + get(dnn::pfCand_chHad_puppiWeightNoLep) = + is_inner ? getValue(candFunc::getPuppiWeightNoLep(chH_cand, default_chH_pw_inner)) + : getValue(candFunc::getPuppiWeightNoLep(chH_cand, default_chH_pw_outer)); + get(dnn::pfCand_chHad_charge) = getValue(chH_cand.charge()); + get(dnn::pfCand_chHad_lostInnerHits) = getValue(candFunc::getLostInnerHits(chH_cand, 0)); + get(dnn::pfCand_chHad_numberOfPixelHits) = + getValueLinear(candFunc::getNumberOfPixelHits(chH_cand, 0), 0, 12, true); get(dnn::pfCand_chHad_vertex_dx) = getValueNorm(pfCands.at(index_chH).vertex().x() - pv.position().x(), 0.0005f, 1.735f); get(dnn::pfCand_chHad_vertex_dy) = @@ -1680,43 +2383,45 @@ class DeepTauId : public deep_tau::DeepTauBase { get(dnn::pfCand_chHad_vertex_dz) = getValueNorm(pfCands.at(index_chH).vertex().z() - pv.position().z(), -0.0201f, 8.333f); get(dnn::pfCand_chHad_vertex_dx_tauFL) = getValueNorm( - pfCands.at(index_chH).vertex().x() - pv.position().x() - tau.flightLength().x(), -0.0014f, 1.93f); + pfCands.at(index_chH).vertex().x() - pv.position().x() - tau_funcs.getFlightLength(tau, tau_index).x(), + -0.0014f, + 1.93f); get(dnn::pfCand_chHad_vertex_dy_tauFL) = getValueNorm( - pfCands.at(index_chH).vertex().y() - pv.position().y() - tau.flightLength().y(), 0.0022f, 1.948f); + pfCands.at(index_chH).vertex().y() - pv.position().y() - tau_funcs.getFlightLength(tau, tau_index).y(), + 0.0022f, + 1.948f); get(dnn::pfCand_chHad_vertex_dz_tauFL) = getValueNorm( - pfCands.at(index_chH).vertex().z() - pv.position().z() - tau.flightLength().z(), -0.0138f, 8.622f); + pfCands.at(index_chH).vertex().z() - pv.position().z() - tau_funcs.getFlightLength(tau, tau_index).z(), + -0.0138f, + 8.622f); - const bool hasTrackDetails = pfCands.at(index_chH).hasTrackDetails(); + const bool hasTrackDetails = candFunc::getHasTrackDetails(chH_cand); if (hasTrackDetails) { get(dnn::pfCand_chHad_hasTrackDetails) = hasTrackDetails; - get(dnn::pfCand_chHad_dxy) = getValueNorm(pfCands.at(index_chH).dxy(), -0.012f, 2.386f); + get(dnn::pfCand_chHad_dxy) = getValueNorm(candFunc::getTauDxy(chH_cand), -0.012f, 2.386f); get(dnn::pfCand_chHad_dxy_sig) = - getValueNorm(std::abs(pfCands.at(index_chH).dxy()) / pfCands.at(index_chH).dxyError(), 6.417f, 36.28f); - get(dnn::pfCand_chHad_dz) = getValueNorm(pfCands.at(index_chH).dz(), -0.0246f, 7.618f); + getValueNorm(std::abs(candFunc::getTauDxy(chH_cand)) / chH_cand.dxyError(), 6.417f, 36.28f); + get(dnn::pfCand_chHad_dz) = getValueNorm(candFunc::getTauDz(chH_cand), -0.0246f, 7.618f); get(dnn::pfCand_chHad_dz_sig) = - getValueNorm(std::abs(pfCands.at(index_chH).dz()) / pfCands.at(index_chH).dzError(), 301.3f, 491.1f); + getValueNorm(std::abs(candFunc::getTauDz(chH_cand)) / chH_cand.dzError(), 301.3f, 491.1f); get(dnn::pfCand_chHad_track_chi2_ndof) = - pfCands.at(index_chH).pseudoTrack().ndof() > 0 - ? getValueNorm(pfCands.at(index_chH).pseudoTrack().chi2() / pfCands.at(index_chH).pseudoTrack().ndof(), + candFunc::getPseudoTrack(chH_cand).ndof() > 0 + ? getValueNorm(candFunc::getPseudoTrack(chH_cand).chi2() / candFunc::getPseudoTrack(chH_cand).ndof(), 0.7876f, 3.694f) : 0; get(dnn::pfCand_chHad_track_ndof) = - pfCands.at(index_chH).pseudoTrack().ndof() > 0 - ? getValueNorm(pfCands.at(index_chH).pseudoTrack().ndof(), 13.92f, 6.581f) + candFunc::getPseudoTrack(chH_cand).ndof() > 0 + ? getValueNorm(candFunc::getPseudoTrack(chH_cand).ndof(), 13.92f, 6.581f) : 0; } - float hcal_fraction = 0.; - if (pfCands.at(index_chH).pdgId() == 1 || pfCands.at(index_chH).pdgId() == 130) { - hcal_fraction = pfCands.at(index_chH).hcalFraction(); - } else if (pfCands.at(index_chH).isIsolatedChargedHadron()) { - hcal_fraction = pfCands.at(index_chH).rawHcalFraction(); - } + float hcal_fraction = candFunc::getHCalFraction(chH_cand, disable_hcalFraction_workaround_); get(dnn::pfCand_chHad_hcalFraction) = getValue(hcal_fraction); - get(dnn::pfCand_chHad_rawCaloFraction) = getValueLinear(pfCands.at(index_chH).rawCaloFraction(), 0.f, 2.6f, true); + get(dnn::pfCand_chHad_rawCaloFraction) = getValueLinear(candFunc::getRawCaloFraction(chH_cand), 0.f, 2.6f, true); } if (valid_nH) { size_t index_nH = cell_map.at(CellObjectType::PfCand_neutralHadron); + const auto& nH_cand = dynamic_cast(pfCands.at(index_nH)); get(dnn::pfCand_nHad_valid) = valid_nH; get(dnn::pfCand_nHad_rel_pt) = getValueNorm(pfCands.at(index_nH).polarP4().pt() / tau.polarP4().pt(), @@ -1728,29 +2433,28 @@ class DeepTauId : public deep_tau::DeepTauBase { false); get(dnn::pfCand_nHad_dphi) = getValueLinear( dPhi(tau.polarP4(), pfCands.at(index_nH).polarP4()), is_inner ? -0.1f : -0.5f, is_inner ? 0.1f : 0.5f, false); - get(dnn::pfCand_nHad_puppiWeight) = getValue(pfCands.at(index_nH).puppiWeight()); - get(dnn::pfCand_nHad_puppiWeightNoLep) = getValue(pfCands.at(index_nH).puppiWeightNoLep()); - float hcal_fraction = 0.; - if (pfCands.at(index_nH).pdgId() == 1 || pfCands.at(index_nH).pdgId() == 130) { - hcal_fraction = pfCands.at(index_nH).hcalFraction(); - } else if (pfCands.at(index_nH).isIsolatedChargedHadron()) { - hcal_fraction = pfCands.at(index_nH).rawHcalFraction(); - } + get(dnn::pfCand_nHad_puppiWeight) = is_inner ? getValue(candFunc::getPuppiWeight(nH_cand, 0.9798355f)) + : getValue(candFunc::getPuppiWeight(nH_cand, 0.7813260f)); + get(dnn::pfCand_nHad_puppiWeightNoLep) = is_inner ? getValue(candFunc::getPuppiWeightNoLep(nH_cand, 0.9046796f)) + : getValue(candFunc::getPuppiWeightNoLep(nH_cand, 0.6554860f)); + float hcal_fraction = candFunc::getHCalFraction(nH_cand, disable_hcalFraction_workaround_); get(dnn::pfCand_nHad_hcalFraction) = getValue(hcal_fraction); } - checkInputs(inputs, is_inner ? "hadron_inner_block" : "hadron_outer_block", dnn::NumberOfInputs); } - template - tensorflow::Tensor createInputsV1(const TauType& tau, - const ElectronCollection& electrons, - const MuonCollection& muons) const { + template + tensorflow::Tensor createInputsV1(const TauCastType& tau, + const size_t tau_index, + const edm::RefToBase tau_ref, + const std::vector* electrons, + const std::vector* muons, + TauFunc tau_funcs) const { static constexpr bool check_all_set = false; static constexpr float default_value_for_set_check = -42; tensorflow::Tensor inputs(tensorflow::DT_FLOAT, {1, dnn_inputs_2017v1::NumberOfInputs}); const auto& get = [&](int var_index) -> float& { return inputs.matrix()(0, var_index); }; - auto leadChargedHadrCand = dynamic_cast(tau.leadChargedHadrCand().get()); + auto leadChargedHadrCand = dynamic_cast(tau.leadChargedHadrCand().get()); if (check_all_set) { for (int var_index = 0; var_index < dnn::NumberOfInputs; ++var_index) { @@ -1762,21 +2466,21 @@ class DeepTauId : public deep_tau::DeepTauBase { get(dnn::eta) = tau.p4().eta(); get(dnn::mass) = tau.p4().mass(); get(dnn::decayMode) = tau.decayMode(); - get(dnn::chargedIsoPtSum) = tau.tauID("chargedIsoPtSum"); - get(dnn::neutralIsoPtSum) = tau.tauID("neutralIsoPtSum"); - get(dnn::neutralIsoPtSumWeight) = tau.tauID("neutralIsoPtSumWeight"); - get(dnn::photonPtSumOutsideSignalCone) = tau.tauID("photonPtSumOutsideSignalCone"); - get(dnn::puCorrPtSum) = tau.tauID("puCorrPtSum"); - get(dnn::dxy) = tau.dxy(); - get(dnn::dxy_sig) = tau.dxy_Sig(); - get(dnn::dz) = leadChargedHadrCand ? leadChargedHadrCand->dz() : default_value; - get(dnn::ip3d) = tau.ip3d(); - get(dnn::ip3d_sig) = tau.ip3d_Sig(); - get(dnn::hasSecondaryVertex) = tau.hasSecondaryVertex(); - get(dnn::flightLength_r) = tau.flightLength().R(); - get(dnn::flightLength_dEta) = dEta(tau.flightLength(), tau.p4()); - get(dnn::flightLength_dPhi) = dPhi(tau.flightLength(), tau.p4()); - get(dnn::flightLength_sig) = tau.flightLengthSig(); + get(dnn::chargedIsoPtSum) = tau_funcs.getChargedIsoPtSum(tau, tau_ref); + get(dnn::neutralIsoPtSum) = tau_funcs.getNeutralIsoPtSum(tau, tau_ref); + get(dnn::neutralIsoPtSumWeight) = tau_funcs.getNeutralIsoPtSumWeight(tau, tau_ref); + get(dnn::photonPtSumOutsideSignalCone) = tau_funcs.getPhotonPtSumOutsideSignalCone(tau, tau_ref); + get(dnn::puCorrPtSum) = tau_funcs.getPuCorrPtSum(tau, tau_ref); + get(dnn::dxy) = tau_funcs.getdxy(tau, tau_index); + get(dnn::dxy_sig) = tau_funcs.getdxySig(tau, tau_index); + get(dnn::dz) = leadChargedHadrCand ? candFunc::getTauDz(*leadChargedHadrCand) : default_value; + get(dnn::ip3d) = tau_funcs.getip3d(tau, tau_index); + get(dnn::ip3d_sig) = tau_funcs.getip3dSig(tau, tau_index); + get(dnn::hasSecondaryVertex) = tau_funcs.getHasSecondaryVertex(tau, tau_index); + get(dnn::flightLength_r) = tau_funcs.getFlightLength(tau, tau_index).R(); + get(dnn::flightLength_dEta) = dEta(tau_funcs.getFlightLength(tau, tau_index), tau.p4()); + get(dnn::flightLength_dPhi) = dPhi(tau_funcs.getFlightLength(tau, tau_index), tau.p4()); + get(dnn::flightLength_sig) = tau_funcs.getFlightLengthSig(tau, tau_index); get(dnn::leadChargedHadrCand_pt) = leadChargedHadrCand ? leadChargedHadrCand->p4().Pt() : default_value; get(dnn::leadChargedHadrCand_dEta) = leadChargedHadrCand ? dEta(leadChargedHadrCand->p4(), tau.p4()) : default_value; @@ -1787,11 +2491,11 @@ class DeepTauId : public deep_tau::DeepTauBase { get(dnn::pt_weighted_dphi_strip) = reco::tau::pt_weighted_dphi_strip(tau, tau.decayMode()); get(dnn::pt_weighted_dr_signal) = reco::tau::pt_weighted_dr_signal(tau, tau.decayMode()); get(dnn::pt_weighted_dr_iso) = reco::tau::pt_weighted_dr_iso(tau, tau.decayMode()); - get(dnn::leadingTrackNormChi2) = tau.leadingTrackNormChi2(); + get(dnn::leadingTrackNormChi2) = tau_funcs.getLeadingTrackNormChi2(tau); get(dnn::e_ratio) = reco::tau::eratio(tau); - get(dnn::gj_angle_diff) = calculateGottfriedJacksonAngleDifference(tau); + get(dnn::gj_angle_diff) = calculateGottfriedJacksonAngleDifference(tau, tau_index, tau_funcs); get(dnn::n_photons) = reco::tau::n_photons_total(tau); - get(dnn::emFraction) = tau.emFraction_MVA(); + get(dnn::emFraction) = tau_funcs.getEmFraction(tau); get(dnn::has_gsf_track) = leadChargedHadrCand && std::abs(leadChargedHadrCand->pdgId()) == 11; get(dnn::inside_ecal_crack) = isInEcalCrack(tau.p4().Eta()); auto gsf_ele = findMatchedElectron(tau, electrons, 0.3); @@ -1834,14 +2538,14 @@ class DeepTauId : public deep_tau::DeepTauBase { get(dnn::gsf_ele_Chi2NormKF) = gsf_ele->closestCtfTrackRef()->normalizedChi2(); get(dnn::gsf_ele_KFNumHits) = gsf_ele->closestCtfTrackRef()->numberOfValidHits(); } - get(dnn::leadChargedCand_etaAtEcalEntrance) = tau.etaAtEcalEntranceLeadChargedCand(); - get(dnn::leadChargedCand_pt) = tau.ptLeadChargedCand(); + get(dnn::leadChargedCand_etaAtEcalEntrance) = tau_funcs.getEtaAtEcalEntrance(tau); + get(dnn::leadChargedCand_pt) = leadChargedHadrCand->pt(); get(dnn::leadChargedHadrCand_HoP) = default_value; get(dnn::leadChargedHadrCand_EoP) = default_value; - if (tau.leadChargedHadrCand()->pt() > 0) { - get(dnn::leadChargedHadrCand_HoP) = tau.hcalEnergyLeadChargedHadrCand() / tau.leadChargedHadrCand()->pt(); - get(dnn::leadChargedHadrCand_EoP) = tau.ecalEnergyLeadChargedHadrCand() / tau.leadChargedHadrCand()->pt(); + if (leadChargedHadrCand->pt() > 0) { + get(dnn::leadChargedHadrCand_HoP) = tau_funcs.getEcalEnergyLeadingChargedHadr(tau) / leadChargedHadrCand->pt(); + get(dnn::leadChargedHadrCand_EoP) = tau_funcs.getHcalEnergyLeadingChargedHadr(tau) / leadChargedHadrCand->pt(); } MuonHitMatchV1 muon_hit_match; @@ -1962,8 +2666,8 @@ class DeepTauId : public deep_tau::DeepTauBase { } } - template - static void processSignalPFComponents(const pat::Tau& tau, + template + static void processSignalPFComponents(const TauCastType& tau, const CandidateCollection& candidates, LorentzVectorXYZ& p4_inner, LorentzVectorXYZ& p4_outer, @@ -2006,8 +2710,8 @@ class DeepTauId : public deep_tau::DeepTauBase { m_outer = n_outer != 0 ? p4_outer.mass() : default_value; } - template - static void processIsolationPFComponents(const pat::Tau& tau, + template + static void processIsolationPFComponents(const TauCastType& tau, const CandidateCollection& candidates, LorentzVectorXYZ& p4, float& pt, @@ -2036,14 +2740,18 @@ class DeepTauId : public deep_tau::DeepTauBase { } // Copied from https://github.com/cms-sw/cmssw/blob/CMSSW_9_4_X/RecoTauTag/RecoTau/plugins/PATTauDiscriminationByMVAIsolationRun2.cc#L218 - static bool calculateGottfriedJacksonAngleDifference(const pat::Tau& tau, double& gj_diff) { - if (tau.hasSecondaryVertex()) { + template + static bool calculateGottfriedJacksonAngleDifference(const TauCastType& tau, + const size_t tau_index, + double& gj_diff, + TauFunc tau_funcs) { + if (tau_funcs.getHasSecondaryVertex(tau, tau_index)) { static constexpr double mTau = 1.77682; const double mAOne = tau.p4().M(); const double pAOneMag = tau.p(); const double argumentThetaGJmax = (std::pow(mTau, 2) - std::pow(mAOne, 2)) / (2 * mTau * pAOneMag); - const double argumentThetaGJmeasured = - tau.p4().Vect().Dot(tau.flightLength()) / (pAOneMag * tau.flightLength().R()); + const double argumentThetaGJmeasured = tau.p4().Vect().Dot(tau_funcs.getFlightLength(tau, tau_index)) / + (pAOneMag * tau_funcs.getFlightLength(tau, tau_index).R()); if (std::abs(argumentThetaGJmax) <= 1. && std::abs(argumentThetaGJmeasured) <= 1.) { double thetaGJmax = std::asin(argumentThetaGJmax); double thetaGJmeasured = std::acos(argumentThetaGJmeasured); @@ -2054,9 +2762,12 @@ class DeepTauId : public deep_tau::DeepTauBase { return false; } - static float calculateGottfriedJacksonAngleDifference(const pat::Tau& tau) { + template + static float calculateGottfriedJacksonAngleDifference(const TauCastType& tau, + const size_t tau_index, + TauFunc tau_funcs) { double gj_diff; - if (calculateGottfriedJacksonAngleDifference(tau, gj_diff)) + if (calculateGottfriedJacksonAngleDifference(tau, tau_index, gj_diff, tau_funcs)) return static_cast(gj_diff); return default_value; } @@ -2066,12 +2777,13 @@ class DeepTauId : public deep_tau::DeepTauBase { return abs_eta > 1.46 && abs_eta < 1.558; } - static const pat::Electron* findMatchedElectron(const pat::Tau& tau, - const pat::ElectronCollection& electrons, + template + static const pat::Electron* findMatchedElectron(const TauCastType& tau, + const std::vector* electrons, double deltaR) { const double dR2 = deltaR * deltaR; const pat::Electron* matched_ele = nullptr; - for (const auto& ele : electrons) { + for (const auto& ele : *electrons) { if (reco::deltaR2(tau.p4(), ele.p4()) < dR2 && (!matched_ele || matched_ele->pt() < ele.pt())) { matched_ele = &ele; } @@ -2080,16 +2792,31 @@ class DeepTauId : public deep_tau::DeepTauBase { } private: - edm::EDGetTokenT electrons_token_; - edm::EDGetTokenT muons_token_; + edm::EDGetTokenT> electrons_token_; + edm::EDGetTokenT> muons_token_; edm::EDGetTokenT rho_token_; + edm::EDGetTokenT basicTauDiscriminators_inputToken_; + edm::EDGetTokenT basicTauDiscriminatorsdR03_inputToken_; + edm::EDGetTokenT>> + pfTauTransverseImpactParameters_token_; std::string input_layer_, output_layer_; - const unsigned version; + const unsigned version_; const int debug_level; const bool disable_dxy_pca_; + const bool disable_hcalFraction_workaround_; + const bool disable_CellIndex_workaround_; std::unique_ptr tauBlockTensor_; std::array, 2> eGammaTensor_, muonTensor_, hadronsTensor_, convTensor_, zeroOutputTensor_; + const bool save_inputs_; + std::ofstream* json_file_; + bool is_first_block_; + int file_counter_; + + //boolean to check if discriminator indices are already mapped + bool discrIndicesMapped_ = false; + std::map basicDiscrIndexMap_; + std::map basicDiscrdR03IndexMap_; }; #include "FWCore/Framework/interface/MakerMacros.h" diff --git a/RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationByIsolation.cc b/RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationByIsolation.cc index 2fc1a15a147b0..5c9e8d78fe204 100644 --- a/RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationByIsolation.cc +++ b/RecoTauTag/RecoTau/plugins/PFRecoTauDiscriminationByIsolation.cc @@ -36,6 +36,8 @@ class PFRecoTauDiscriminationByIsolation : public PFTauDiscriminationProducerBas calculateWeights_ = pset.getParameter("ApplyDiscriminationByWeightedECALIsolation"); + enableHGCalWorkaround_ = pset.getParameter("enableHGCalWorkaround"); + // RIC: multiply neutral isolation by a flat factor. // Useful, for instance, to combine charged and neutral isolations // with different relative weights @@ -199,6 +201,7 @@ class PFRecoTauDiscriminationByIsolation : public PFTauDiscriminationProducerBas bool includeTracks_; bool includeGammas_; bool calculateWeights_; + bool enableHGCalWorkaround_; double weightGammas_; bool applyOccupancyCut_; uint32_t maximumOccupancy_; @@ -498,7 +501,22 @@ double PFRecoTauDiscriminationByIsolation::discriminate(const PFTauRef& pfTau) c double neutralPt = 0.; double weightedNeutralPt = 0.; for (auto const& isoObject : isoCharged_) { - chargedPt += isoObject->pt(); + //------------------------------------------------------------------------- + // CV: fix for Phase-2 HLT tau trigger studies + // (pT of PFCandidates within HGCal acceptance is significantly higher than track pT !!) + if (enableHGCalWorkaround_) { + double trackPt = (isoObject->bestTrack()) ? isoObject->bestTrack()->pt() : 0.; + double pfCandPt = isoObject->pt(); + if (pfCandPt > trackPt) { + chargedPt += trackPt; + neutralPt += pfCandPt - trackPt; + } else { + chargedPt += isoObject->pt(); + } + } else { + chargedPt += isoObject->pt(); + } + //------------------------------------------------------------------------- } if (!calculateWeights_) { for (auto const& isoObject : isoNeutral_) { @@ -657,6 +675,7 @@ void PFRecoTauDiscriminationByIsolation::fillDescriptions(edm::ConfigurationDesc desc.add("ApplyDiscriminationByTrackerIsolation", true); desc.add("storeRawPhotonSumPt_outsideSignalCone", false); desc.add("rhoProducer", edm::InputTag("fixedGridRhoFastjetAll")); + desc.add("enableHGCalWorkaround", false); { edm::ParameterSetDescription vpsd1; diff --git a/RecoTauTag/RecoTau/python/tools/runTauIdMVA.py b/RecoTauTag/RecoTau/python/tools/runTauIdMVA.py index 311ad4065215c..60dcb42c3468f 100644 --- a/RecoTauTag/RecoTau/python/tools/runTauIdMVA.py +++ b/RecoTauTag/RecoTau/python/tools/runTauIdMVA.py @@ -701,7 +701,8 @@ def runTauID(self): mem_mapped = cms.bool(False), version = cms.uint32(self.getDeepTauVersion(file_names[0])[1]), debug_level = cms.int32(0), - disable_dxy_pca = cms.bool(True) + disable_dxy_pca = cms.bool(True), + is_online = cms.bool(False) ) self.processDeepProducer('deepTau2017v2p1', tauIDSources, workingPoints_) diff --git a/RecoTauTag/RecoTau/src/DeepTauBase.cc b/RecoTauTag/RecoTau/src/DeepTauBase.cc index e3ec63b3da2a4..4db4610cf2706 100644 --- a/RecoTauTag/RecoTau/src/DeepTauBase.cc +++ b/RecoTauTag/RecoTau/src/DeepTauBase.cc @@ -7,6 +7,9 @@ * \author Maria Rosaria Di Domenico, University of Siena & INFN Pisa */ +//TODO: port to offline RECO/AOD inputs to allow usage with offline AOD +//TODO: Take into account that PFTaus can also be build with pat::PackedCandidates + #include "RecoTauTag/RecoTau/interface/DeepTauBase.h" namespace deep_tau { @@ -36,10 +39,14 @@ namespace deep_tau { } } - double TauWPThreshold::operator()(const pat::Tau& tau) const { + double TauWPThreshold::operator()(const reco::BaseTau& tau, bool isPFTau) const { if (!fn_) return value_; - fn_->SetParameter(0, tau.decayMode()); + + if (isPFTau) + fn_->SetParameter(0, dynamic_cast(tau).decayMode()); + else + fn_->SetParameter(0, dynamic_cast(tau).decayMode()); fn_->SetParameter(1, tau.pt()); fn_->SetParameter(2, tau.eta()); return fn_->Eval(0); @@ -47,7 +54,8 @@ namespace deep_tau { std::unique_ptr DeepTauBase::Output::get_value(const edm::Handle& taus, const tensorflow::Tensor& pred, - const WPList& working_points) const { + const WPList* working_points, + bool is_online) const { std::vector outputbuffer(taus->size()); for (size_t tau_index = 0; tau_index < taus->size(); ++tau_index) { @@ -61,9 +69,11 @@ namespace deep_tau { x = den_val != 0 ? x / den_val : std::numeric_limits::max(); } outputbuffer[tau_index].rawValues.push_back(x); - for (const auto& wp : working_points) { - const bool pass = x > (*wp)(taus->at(tau_index)); - outputbuffer[tau_index].workingPoints.push_back(pass); + if (working_points) { + for (const auto& wp : *working_points) { + const bool pass = x > (*wp)(taus->at(tau_index), is_online); + outputbuffer[tau_index].workingPoints.push_back(pass); + } } } std::unique_ptr output = std::make_unique(); @@ -77,8 +87,9 @@ namespace deep_tau { const OutputCollection& outputCollection, const DeepTauCache* cache) : tausToken_(consumes(cfg.getParameter("taus"))), - pfcandToken_(consumes(cfg.getParameter("pfcands"))), + pfcandToken_(consumes(cfg.getParameter("pfcands"))), vtxToken_(consumes(cfg.getParameter("vertices"))), + is_online_(cfg.getParameter("is_online")), outputs_(outputCollection), cache_(cache) { for (const auto& output_desc : outputs_) { @@ -88,19 +99,88 @@ namespace deep_tau { workingPoints_[output_desc.first].push_back(std::make_unique(cut_str)); } } + + // prediscriminant operator + // require the tau to pass the following prediscriminants + const edm::ParameterSet& prediscriminantConfig = cfg.getParameter("Prediscriminants"); + + // determine boolean operator used on the prediscriminants + std::string pdBoolOperator = prediscriminantConfig.getParameter("BooleanOperator"); + // convert string to lowercase + transform(pdBoolOperator.begin(), pdBoolOperator.end(), pdBoolOperator.begin(), ::tolower); + + if (pdBoolOperator == "and") { + andPrediscriminants_ = 0x1; //use chars instead of bools so we can do a bitwise trick later + } else if (pdBoolOperator == "or") { + andPrediscriminants_ = 0x0; + } else { + throw cms::Exception("TauDiscriminationProducerBase") + << "PrediscriminantBooleanOperator defined incorrectly, options are: AND,OR"; + } + + // get the list of prediscriminants + std::vector prediscriminantsNames = + prediscriminantConfig.getParameterNamesForType(); + + for (auto const& iDisc : prediscriminantsNames) { + const edm::ParameterSet& iPredisc = prediscriminantConfig.getParameter(iDisc); + const edm::InputTag& label = iPredisc.getParameter("Producer"); + double cut = iPredisc.getParameter("cut"); + + if (is_online_) { + TauDiscInfo thisDiscriminator; + thisDiscriminator.label = label; + thisDiscriminator.cut = cut; + thisDiscriminator.disc_token = consumes(label); + recoPrediscriminants_.push_back(thisDiscriminator); + } else { + TauDiscInfo thisDiscriminator; + thisDiscriminator.label = label; + thisDiscriminator.cut = cut; + thisDiscriminator.disc_token = consumes(label); + patPrediscriminants_.push_back(thisDiscriminator); + } + } } void DeepTauBase::produce(edm::Event& event, const edm::EventSetup& es) { edm::Handle taus; event.getByToken(tausToken_, taus); + edm::ProductID tauProductID = taus.id(); + + // load prediscriminators + size_t nPrediscriminants = + patPrediscriminants_.empty() ? recoPrediscriminants_.size() : patPrediscriminants_.size(); + for (size_t iDisc = 0; iDisc < nPrediscriminants; ++iDisc) { + edm::ProductID discKeyId; + if (is_online_) { + recoPrediscriminants_[iDisc].fill(event); + discKeyId = recoPrediscriminants_[iDisc].handle->keyProduct().id(); + } else { + patPrediscriminants_[iDisc].fill(event); + discKeyId = patPrediscriminants_[iDisc].handle->keyProduct().id(); + } - const tensorflow::Tensor& pred = getPredictions(event, es, taus); + // Check to make sure the product is correct for the discriminator. + // If not, throw a more informative exception. + if (tauProductID != discKeyId) { + throw cms::Exception("MisconfiguredPrediscriminant") + << "The tau collection has product ID: " << tauProductID + << " but the pre-discriminator is keyed with product ID: " << discKeyId << std::endl; + } + } + + const tensorflow::Tensor& pred = getPredictions(event, taus); createOutputs(event, pred, taus); } void DeepTauBase::createOutputs(edm::Event& event, const tensorflow::Tensor& pred, edm::Handle taus) { for (const auto& output_desc : outputs_) { - auto result = output_desc.second.get_value(taus, pred, workingPoints_.at(output_desc.first)); + const WPList* working_points = nullptr; + if (workingPoints_.find(output_desc.first) != workingPoints_.end()) { + working_points = &workingPoints_.at(output_desc.first); + } + auto result = output_desc.second.get_value(taus, pred, working_points, is_online_); event.put(std::move(result), output_desc.first); } } diff --git a/RecoTauTag/RecoTau/src/PositionAtECalEntranceComputer.cc b/RecoTauTag/RecoTau/src/PositionAtECalEntranceComputer.cc index aedab0dab18fb..391ac598bb414 100644 --- a/RecoTauTag/RecoTau/src/PositionAtECalEntranceComputer.cc +++ b/RecoTauTag/RecoTau/src/PositionAtECalEntranceComputer.cc @@ -16,7 +16,8 @@ void PositionAtECalEntranceComputer::beginEvent(const edm::EventSetup& es) { bField_z_ = bField->inTesla(GlobalPoint(0., 0., 0.)).z(); } -reco::Candidate::Point PositionAtECalEntranceComputer::operator()(const reco::Candidate* particle, bool& success) const { +reco::Candidate::Point PositionAtECalEntranceComputer::operator()(const reco::Candidate* particle, + bool& success) const { assert(bField_z_ != -1.); BaseParticlePropagator propagator = BaseParticlePropagator( RawParticle(particle->p4(),