Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-17 01:47:56

0001 
0002 /** \class PFRecoTauDiscriminationAgainstMuonMVA
0003  *
0004  * MVA based discriminator against muon -> tau fakes
0005  * 
0006  * \author Christian Veelken, LLR
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     //const GBRForest* mva = dynamic_cast<GBRForest*>(inputFile->Get(mvaName.data())); // CV: dynamic_cast<GBRForest*> fails for some reason ?!
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;  // CV: map into range 0..3
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   }  // namespace
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     // CV: define dummy category index in order to use RecoTauDiscriminantCutMultiplexer module to apply WP cuts
0171     result.rawValues = {-1., 0.};
0172 
0173     // CV: computation of anti-muon MVA value requires presence of leading charged hadron
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 }  // namespace
0234 
0235 void PFRecoTauDiscriminationAgainstMuonMVA::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0236   // pfRecoTauDiscriminationAgainstMuonMVA
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);