File indexing completed on 2023-03-17 11:21:50
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "RecoTauTag/RecoTau/interface/TauDiscriminationProducerBase.h"
0011
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Utilities/interface/InputTag.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/ParameterSet/interface/FileInPath.h"
0017
0018 #include "FWCore/Utilities/interface/Exception.h"
0019
0020 #include <FWCore/ParameterSet/interface/ConfigurationDescriptions.h>
0021 #include <FWCore/ParameterSet/interface/ParameterSetDescription.h>
0022
0023 #include "DataFormats/Candidate/interface/Candidate.h"
0024 #include "DataFormats/TauReco/interface/PFTau.h"
0025 #include "DataFormats/TauReco/interface/PFTauFwd.h"
0026 #include "DataFormats/MuonReco/interface/Muon.h"
0027 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0028 #include "DataFormats/Math/interface/deltaR.h"
0029
0030 #include "CondFormats/GBRForest/interface/GBRForest.h"
0031 #include "CondFormats/DataRecord/interface/GBRWrapperRcd.h"
0032
0033 #include <TMath.h>
0034 #include <TFile.h>
0035
0036 #include <iostream>
0037
0038 using namespace reco;
0039
0040 namespace {
0041 const GBRForest* loadMVAfromFile(const edm::FileInPath& inputFileName,
0042 const std::string& mvaName,
0043 std::vector<TFile*>& inputFilesToDelete) {
0044 if (inputFileName.location() == edm::FileInPath::Unknown)
0045 throw cms::Exception("PFRecoTauDiscriminationAgainstMuonMVA::loadMVA")
0046 << " Failed to find File = " << inputFileName << " !!\n";
0047 TFile* inputFile = new TFile(inputFileName.fullPath().data());
0048
0049
0050 const GBRForest* mva = (GBRForest*)inputFile->Get(mvaName.data());
0051 if (!mva)
0052 throw cms::Exception("PFRecoTauDiscriminationAgainstMuonMVA::loadMVA")
0053 << " Failed to load MVA = " << mvaName.data() << " from file = " << inputFileName.fullPath().data()
0054 << " !!\n";
0055
0056 inputFilesToDelete.push_back(inputFile);
0057
0058 return mva;
0059 }
0060
0061 class PFRecoTauDiscriminationAgainstMuonMVA final : public PFTauDiscriminationContainerProducerBase {
0062 public:
0063 explicit PFRecoTauDiscriminationAgainstMuonMVA(const edm::ParameterSet& cfg)
0064 : PFTauDiscriminationContainerProducerBase(cfg),
0065 moduleLabel_(cfg.getParameter<std::string>("@module_label")),
0066 mvaReader_(nullptr),
0067 mvaInput_(nullptr) {
0068 mvaName_ = cfg.getParameter<std::string>("mvaName");
0069 loadMVAfromDB_ = cfg.getParameter<bool>("loadMVAfromDB");
0070 if (!loadMVAfromDB_) {
0071 inputFileName_ = cfg.getParameter<edm::FileInPath>("inputFileName");
0072 } else {
0073 mvaToken_ = esConsumes(edm::ESInputTag{"", mvaName_});
0074 }
0075 mvaInput_ = new float[11];
0076
0077 srcMuons_ = cfg.getParameter<edm::InputTag>("srcMuons");
0078 Muons_token = consumes<reco::MuonCollection>(srcMuons_);
0079 dRmuonMatch_ = cfg.getParameter<double>("dRmuonMatch");
0080
0081 verbosity_ = cfg.getParameter<int>("verbosity");
0082 }
0083
0084 void beginEvent(const edm::Event&, const edm::EventSetup&) override;
0085
0086 reco::SingleTauDiscriminatorContainer discriminate(const PFTauRef&) const override;
0087
0088 ~PFRecoTauDiscriminationAgainstMuonMVA() override {
0089 if (!loadMVAfromDB_)
0090 delete mvaReader_;
0091 delete[] mvaInput_;
0092 for (std::vector<TFile*>::iterator it = inputFilesToDelete_.begin(); it != inputFilesToDelete_.end(); ++it) {
0093 delete (*it);
0094 }
0095 }
0096
0097 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0098
0099 private:
0100 std::string moduleLabel_;
0101
0102 std::string mvaName_;
0103 edm::ESGetToken<GBRForest, GBRWrapperRcd> mvaToken_;
0104 bool loadMVAfromDB_;
0105 edm::FileInPath inputFileName_;
0106 const GBRForest* mvaReader_;
0107 float* mvaInput_;
0108
0109 edm::InputTag srcMuons_;
0110 edm::Handle<reco::MuonCollection> muons_;
0111 edm::EDGetTokenT<reco::MuonCollection> Muons_token;
0112 double dRmuonMatch_;
0113
0114 edm::Handle<TauCollection> taus_;
0115
0116 std::vector<TFile*> inputFilesToDelete_;
0117
0118 int verbosity_;
0119 };
0120
0121 void PFRecoTauDiscriminationAgainstMuonMVA::beginEvent(const edm::Event& evt, const edm::EventSetup& es) {
0122 if (!mvaReader_) {
0123 if (loadMVAfromDB_) {
0124 mvaReader_ = &es.getData(mvaToken_);
0125 } else {
0126 mvaReader_ = loadMVAfromFile(inputFileName_, mvaName_, inputFilesToDelete_);
0127 }
0128 }
0129
0130 evt.getByToken(Muons_token, muons_);
0131
0132 evt.getByToken(Tau_token, taus_);
0133 }
0134
0135 namespace {
0136 void countHits(const reco::Muon& muon,
0137 std::vector<int>& numHitsDT,
0138 std::vector<int>& numHitsCSC,
0139 std::vector<int>& numHitsRPC) {
0140 if (muon.outerTrack().isNonnull()) {
0141 const reco::HitPattern& muonHitPattern = muon.outerTrack()->hitPattern();
0142 for (int iHit = 0; iHit < muonHitPattern.numberOfAllHits(HitPattern::TRACK_HITS); ++iHit) {
0143 uint32_t hit = muonHitPattern.getHitPattern(HitPattern::TRACK_HITS, iHit);
0144 if (hit == 0)
0145 break;
0146 if (muonHitPattern.muonHitFilter(hit) && (muonHitPattern.getHitType(hit) == TrackingRecHit::valid ||
0147 muonHitPattern.getHitType(hit) == TrackingRecHit::bad)) {
0148 int muonStation = muonHitPattern.getMuonStation(hit) - 1;
0149 if (muonStation >= 0 && muonStation < 4) {
0150 if (muonHitPattern.muonDTHitFilter(hit))
0151 ++numHitsDT[muonStation];
0152 else if (muonHitPattern.muonCSCHitFilter(hit))
0153 ++numHitsCSC[muonStation];
0154 else if (muonHitPattern.muonRPCHitFilter(hit))
0155 ++numHitsRPC[muonStation];
0156 }
0157 }
0158 }
0159 }
0160 }
0161 }
0162
0163 reco::SingleTauDiscriminatorContainer PFRecoTauDiscriminationAgainstMuonMVA::discriminate(const PFTauRef& tau) const {
0164 if (verbosity_) {
0165 edm::LogPrint("PFTauAgainstMuonMVA") << "<PFRecoTauDiscriminationAgainstMuonMVA::discriminate>:";
0166 edm::LogPrint("PFTauAgainstMuonMVA") << " moduleLabel = " << moduleLabel_;
0167 }
0168
0169 reco::SingleTauDiscriminatorContainer result;
0170
0171 result.rawValues = {-1., 0.};
0172
0173
0174 if (tau->leadPFChargedHadrCand().isNull())
0175 return 0.;
0176
0177 mvaInput_[0] = TMath::Abs(tau->eta());
0178 double tauCaloEnECAL = 0.;
0179 double tauCaloEnHCAL = 0.;
0180 const std::vector<reco::PFCandidatePtr>& tauSignalPFCands = tau->signalPFCands();
0181 for (std::vector<reco::PFCandidatePtr>::const_iterator tauSignalPFCand = tauSignalPFCands.begin();
0182 tauSignalPFCand != tauSignalPFCands.end();
0183 ++tauSignalPFCand) {
0184 tauCaloEnECAL += (*tauSignalPFCand)->ecalEnergy();
0185 tauCaloEnHCAL += (*tauSignalPFCand)->hcalEnergy();
0186 }
0187 mvaInput_[1] = TMath::Sqrt(TMath::Max(0., tauCaloEnECAL));
0188 mvaInput_[2] = TMath::Sqrt(TMath::Max(0., tauCaloEnHCAL));
0189 mvaInput_[3] = tau->leadPFChargedHadrCand()->pt() / TMath::Max(1., Double_t(tau->pt()));
0190 mvaInput_[4] = TMath::Sqrt(TMath::Max(0., tau->leadPFChargedHadrCand()->ecalEnergy()));
0191 mvaInput_[5] = TMath::Sqrt(TMath::Max(0., tau->leadPFChargedHadrCand()->hcalEnergy()));
0192 int numMatches = 0;
0193 std::vector<int> numHitsDT(4);
0194 std::vector<int> numHitsCSC(4);
0195 std::vector<int> numHitsRPC(4);
0196 for (int iStation = 0; iStation < 4; ++iStation) {
0197 numHitsDT[iStation] = 0;
0198 numHitsCSC[iStation] = 0;
0199 numHitsRPC[iStation] = 0;
0200 }
0201 if (tau->leadPFChargedHadrCand().isNonnull()) {
0202 reco::MuonRef muonRef = tau->leadPFChargedHadrCand()->muonRef();
0203 if (muonRef.isNonnull()) {
0204 numMatches = muonRef->numberOfMatches(reco::Muon::NoArbitration);
0205 countHits(*muonRef, numHitsDT, numHitsCSC, numHitsRPC);
0206 }
0207 }
0208 size_t numMuons = muons_->size();
0209 for (size_t idxMuon = 0; idxMuon < numMuons; ++idxMuon) {
0210 reco::MuonRef muon(muons_, idxMuon);
0211 if (tau->leadPFChargedHadrCand().isNonnull() && tau->leadPFChargedHadrCand()->muonRef().isNonnull() &&
0212 muon == tau->leadPFChargedHadrCand()->muonRef()) {
0213 continue;
0214 }
0215 double dR = deltaR(muon->p4(), tau->p4());
0216 if (dR < dRmuonMatch_) {
0217 numMatches += muon->numberOfMatches(reco::Muon::NoArbitration);
0218 countHits(*muon, numHitsDT, numHitsCSC, numHitsRPC);
0219 }
0220 }
0221 mvaInput_[6] = numMatches;
0222 mvaInput_[7] = numHitsDT[0] + numHitsCSC[0] + numHitsRPC[0];
0223 mvaInput_[8] = numHitsDT[1] + numHitsCSC[1] + numHitsRPC[1];
0224 mvaInput_[9] = numHitsDT[2] + numHitsCSC[2] + numHitsRPC[2];
0225 mvaInput_[10] = numHitsDT[3] + numHitsCSC[3] + numHitsRPC[3];
0226
0227 result.rawValues.at(0) = mvaReader_->GetClassifier(mvaInput_);
0228 if (verbosity_) {
0229 edm::LogPrint("PFTauAgainstMuonMVA") << "mvaValue = " << result.rawValues.at(0);
0230 }
0231 return result.rawValues.at(0);
0232 }
0233 }
0234
0235 void PFRecoTauDiscriminationAgainstMuonMVA::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0236
0237 edm::ParameterSetDescription desc;
0238 desc.add<double>("mvaMin", 0.0);
0239 desc.add<std::string>("mvaName", "againstMuonMVA");
0240 desc.add<edm::InputTag>("PFTauProducer", edm::InputTag("pfTauProducer"));
0241 desc.add<int>("verbosity", 0);
0242 desc.add<bool>("returnMVA", true);
0243 desc.add<edm::FileInPath>("inputFileName", edm::FileInPath("RecoTauTag/RecoTau/data/emptyMVAinputFile"));
0244 desc.add<bool>("loadMVAfromDB", true);
0245 {
0246 edm::ParameterSetDescription psd0;
0247 psd0.add<std::string>("BooleanOperator", "and");
0248 {
0249 edm::ParameterSetDescription psd1;
0250 psd1.add<double>("cut");
0251 psd1.add<edm::InputTag>("Producer");
0252 psd0.addOptional<edm::ParameterSetDescription>("leadTrack", psd1);
0253 }
0254 desc.add<edm::ParameterSetDescription>("Prediscriminants", psd0);
0255 }
0256 desc.add<double>("dRmuonMatch", 0.3);
0257 desc.add<edm::InputTag>("srcMuons", edm::InputTag("muons"));
0258 descriptions.add("pfRecoTauDiscriminationAgainstMuonMVA", desc);
0259 }
0260
0261 DEFINE_FWK_MODULE(PFRecoTauDiscriminationAgainstMuonMVA);
