RecoTau/plugins/PFRecoTauDiscriminationByMVAIsolation2.cc

0001
0002 /** \class PFRecoTauDiscriminationByMVAIsolation2
0003  *
0004  * MVA based discriminator against jet -> 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/TauReco/interface/PFTauDiscriminator.h"
0027 #include "DataFormats/TauReco/interface/PFTauTransverseImpactParameterAssociation.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("PFRecoTauDiscriminationByIsolationMVA2::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("PFRecoTauDiscriminationByIsolationMVA2::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 }  // namespace
0061
0062 class PFRecoTauDiscriminationByIsolationMVA2 : public PFTauDiscriminationContainerProducerBase {
0063 public:
0064   explicit PFRecoTauDiscriminationByIsolationMVA2(const edm::ParameterSet& cfg)
0065       : PFTauDiscriminationContainerProducerBase(cfg),
0066         moduleLabel_(cfg.getParameter<std::string>("@module_label")),
0067         mvaReader_(nullptr),
0068         mvaInput_(nullptr) {
0069     mvaName_ = cfg.getParameter<std::string>("mvaName");
0070     loadMVAfromDB_ = cfg.getParameter<bool>("loadMVAfromDB");
0071     if (!loadMVAfromDB_) {
0072       inputFileName_ = cfg.getParameter<edm::FileInPath>("inputFileName");
0073     } else {
0074       mvaToken_ = esConsumes(edm::ESInputTag{"", mvaName_});
0075     }
0076     std::string mvaOpt_string = cfg.getParameter<std::string>("mvaOpt");
0077     if (mvaOpt_string == "oldDMwoLT")
0078       mvaOpt_ = kOldDMwoLT;
0079     else if (mvaOpt_string == "oldDMwLT")
0080       mvaOpt_ = kOldDMwLT;
0081     else if (mvaOpt_string == "newDMwoLT")
0082       mvaOpt_ = kNewDMwoLT;
0083     else if (mvaOpt_string == "newDMwLT")
0084       mvaOpt_ = kNewDMwLT;
0085     else
0086       throw cms::Exception("PFRecoTauDiscriminationByIsolationMVA2")
0087           << " Invalid Configuration Parameter 'mvaOpt' = " << mvaOpt_string << " !!\n";
0088
0089     if (mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT)
0090       mvaInput_ = new float[6];
0091     else if (mvaOpt_ == kOldDMwLT || mvaOpt_ == kNewDMwLT)
0092       mvaInput_ = new float[12];
0093     else
0094       assert(0);
0095
0096     tauTransverseImpactParameters_token_ =
0097         consumes<PFTauTIPAssociationByRef>(cfg.getParameter<edm::InputTag>("srcTauTransverseImpactParameters"));
0098
0099     basicTauDiscriminators_token_ =
0100         consumes<reco::TauDiscriminatorContainer>(cfg.getParameter<edm::InputTag>("srcBasicTauDiscriminators"));
0101     chargedIsoPtSum_index_ = cfg.getParameter<int>("srcChargedIsoPtSumIndex");
0102     neutralIsoPtSum_index_ = cfg.getParameter<int>("srcNeutralIsoPtSumIndex");
0103     pucorrPtSum_index_ = cfg.getParameter<int>("srcPUcorrPtSumIndex");
0104
0105     verbosity_ = cfg.getParameter<int>("verbosity");
0106   }
0107
0108   void beginEvent(const edm::Event&, const edm::EventSetup&) override;
0109
0110   reco::SingleTauDiscriminatorContainer discriminate(const PFTauRef&) const override;
0111
0112   ~PFRecoTauDiscriminationByIsolationMVA2() override {
0113     if (!loadMVAfromDB_)
0114       delete mvaReader_;
0115     delete[] mvaInput_;
0116     for (std::vector<TFile*>::iterator it = inputFilesToDelete_.begin(); it != inputFilesToDelete_.end(); ++it) {
0117       delete (*it);
0118     }
0119   }
0120
0121   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0122
0123 private:
0124   std::string moduleLabel_;
0125
0126   std::string mvaName_;
0127   edm::ESGetToken<GBRForest, GBRWrapperRcd> mvaToken_;
0128   bool loadMVAfromDB_;
0129   edm::FileInPath inputFileName_;
0130   const GBRForest* mvaReader_;
0131   enum { kOldDMwoLT, kOldDMwLT, kNewDMwoLT, kNewDMwLT };
0132   int mvaOpt_;
0133   float* mvaInput_;
0134
0135   typedef edm::AssociationVector<reco::PFTauRefProd, std::vector<reco::PFTauTransverseImpactParameterRef> >
0136       PFTauTIPAssociationByRef;
0137   edm::EDGetTokenT<PFTauTIPAssociationByRef> tauTransverseImpactParameters_token_;
0138   edm::Handle<PFTauTIPAssociationByRef> tauLifetimeInfos_;
0139
0140   edm::EDGetTokenT<reco::TauDiscriminatorContainer> basicTauDiscriminators_token_;
0141   edm::Handle<reco::TauDiscriminatorContainer> basicTauDiscriminators_;
0142   int chargedIsoPtSum_index_;
0143   int neutralIsoPtSum_index_;
0144   int pucorrPtSum_index_;
0145
0146   edm::Handle<TauCollection> taus_;
0147
0148   std::vector<TFile*> inputFilesToDelete_;
0149
0150   int verbosity_;
0151 };
0152
0153 void PFRecoTauDiscriminationByIsolationMVA2::beginEvent(const edm::Event& evt, const edm::EventSetup& es) {
0154   if (!mvaReader_) {
0155     if (loadMVAfromDB_) {
0156       mvaReader_ = &es.getData(mvaToken_);
0157     } else {
0158       mvaReader_ = loadMVAfromFile(inputFileName_, mvaName_, inputFilesToDelete_);
0159     }
0160   }
0161
0162   evt.getByToken(tauTransverseImpactParameters_token_, tauLifetimeInfos_);
0163
0164   evt.getByToken(basicTauDiscriminators_token_, basicTauDiscriminators_);
0165
0166   evt.getByToken(Tau_token, taus_);
0167 }
0168
0169 reco::SingleTauDiscriminatorContainer PFRecoTauDiscriminationByIsolationMVA2::discriminate(const PFTauRef& tau) const {
0170   reco::SingleTauDiscriminatorContainer result;
0171   // CV: define dummy category index in order to use RecoTauDiscriminantCutMultiplexer module to apply WP cuts
0172   result.rawValues = {-1., 0.};
0173
0174   // CV: computation of MVA value requires presence of leading charged hadron
0175   if (tau->leadChargedHadrCand().isNull())
0176     return 0.;
0177
0178   int tauDecayMode = tau->decayMode();
0179
0180   if (((mvaOpt_ == kOldDMwoLT || mvaOpt_ == kOldDMwLT) &&
0181        (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 10)) ||
0182       ((mvaOpt_ == kNewDMwoLT || mvaOpt_ == kNewDMwLT) &&
0183        (tauDecayMode == 0 || tauDecayMode == 1 || tauDecayMode == 2 || tauDecayMode == 5 || tauDecayMode == 6 ||
0184         tauDecayMode == 10))) {
0185     double chargedIsoPtSum = (*basicTauDiscriminators_)[tau].rawValues.at(chargedIsoPtSum_index_);
0186     double neutralIsoPtSum = (*basicTauDiscriminators_)[tau].rawValues.at(neutralIsoPtSum_index_);
0187     double puCorrPtSum = (*basicTauDiscriminators_)[tau].rawValues.at(pucorrPtSum_index_);
0188
0189     const reco::PFTauTransverseImpactParameter& tauLifetimeInfo = *(*tauLifetimeInfos_)[tau];
0190
0191     double decayDistX = tauLifetimeInfo.flightLength().x();
0192     double decayDistY = tauLifetimeInfo.flightLength().y();
0193     double decayDistZ = tauLifetimeInfo.flightLength().z();
0194     double decayDistMag = TMath::Sqrt(decayDistX * decayDistX + decayDistY * decayDistY + decayDistZ * decayDistZ);
0195
0196     if (mvaOpt_ == kOldDMwoLT || mvaOpt_ == kNewDMwoLT) {
0197       mvaInput_[0] = TMath::Log(TMath::Max(1., Double_t(tau->pt())));
0198       mvaInput_[1] = TMath::Abs(tau->eta());
0199       mvaInput_[2] = TMath::Log(TMath::Max(1.e-2, chargedIsoPtSum));
0200       mvaInput_[3] = TMath::Log(TMath::Max(1.e-2, neutralIsoPtSum - 0.125 * puCorrPtSum));
0201       mvaInput_[4] = TMath::Log(TMath::Max(1.e-2, puCorrPtSum));
0202       mvaInput_[5] = tauDecayMode;
0203     } else if (mvaOpt_ == kOldDMwLT || mvaOpt_ == kNewDMwLT) {
0204       mvaInput_[0] = TMath::Log(TMath::Max(1., Double_t(tau->pt())));
0205       mvaInput_[1] = TMath::Abs(tau->eta());
0206       mvaInput_[2] = TMath::Log(TMath::Max(1.e-2, chargedIsoPtSum));
0207       mvaInput_[3] = TMath::Log(TMath::Max(1.e-2, neutralIsoPtSum - 0.125 * puCorrPtSum));
0208       mvaInput_[4] = TMath::Log(TMath::Max(1.e-2, puCorrPtSum));
0209       mvaInput_[5] = tauDecayMode;
0210       mvaInput_[6] = TMath::Sign(+1., tauLifetimeInfo.dxy());
0211       mvaInput_[7] = TMath::Sqrt(TMath::Abs(TMath::Min(1., tauLifetimeInfo.dxy())));
0212       mvaInput_[8] = TMath::Min(10., TMath::Abs(tauLifetimeInfo.dxy_Sig()));
0213       mvaInput_[9] = (tauLifetimeInfo.hasSecondaryVertex()) ? 1. : 0.;
0214       mvaInput_[10] = TMath::Sqrt(decayDistMag);
0215       mvaInput_[11] = TMath::Min(10., tauLifetimeInfo.flightLengthSig());
0216     }
0217
0218     double mvaValue = mvaReader_->GetClassifier(mvaInput_);
0219     if (verbosity_) {
0220       edm::LogPrint("PFTauDiscByMVAIsol2") << "<PFRecoTauDiscriminationByIsolationMVA2::discriminate>:";
0221       edm::LogPrint("PFTauDiscByMVAIsol2") << " tau: Pt = " << tau->pt() << ", eta = " << tau->eta();
0222       edm::LogPrint("PFTauDiscByMVAIsol2") << " isolation: charged = " << chargedIsoPtSum
0223                                            << ", neutral = " << neutralIsoPtSum << ", PUcorr = " << puCorrPtSum;
0224       edm::LogPrint("PFTauDiscByMVAIsol2") << " decay mode = " << tauDecayMode;
0225       edm::LogPrint("PFTauDiscByMVAIsol2") << " impact parameter: distance = " << tauLifetimeInfo.dxy()
0226                                            << ", significance = " << tauLifetimeInfo.dxy_Sig();
0227       edm::LogPrint("PFTauDiscByMVAIsol2")
0228           << " has decay vertex = " << tauLifetimeInfo.hasSecondaryVertex() << ":"
0229           << " distance = " << decayDistMag << ", significance = " << tauLifetimeInfo.flightLengthSig();
0230       edm::LogPrint("PFTauDiscByMVAIsol2") << "--> mvaValue = " << mvaValue;
0231     }
0232     result.rawValues.at(0) = mvaValue;
0233   }
0234   return result;
0235 }
0236
0237 void PFRecoTauDiscriminationByIsolationMVA2::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0238   // pfRecoTauDiscriminationByIsolationMVA2
0239   edm::ParameterSetDescription desc;
0240
0241   desc.add<std::string>("mvaName");
0242   desc.add<bool>("loadMVAfromDB");
0243   desc.addOptional<edm::FileInPath>("inputFileName");
0244   desc.add<std::string>("mvaOpt");
0245
0246   desc.add<edm::InputTag>("srcTauTransverseImpactParameters");
0247   desc.add<edm::InputTag>("srcBasicTauDiscriminators");
0248   desc.add<int>("srcChargedIsoPtSumIndex");
0249   desc.add<int>("srcNeutralIsoPtSumIndex");
0250   desc.add<int>("srcPUcorrPtSumIndex");
0251   desc.add<int>("verbosity", 0);
0252
0253   fillProducerDescriptions(desc);  // inherited from the base
0254
0255   descriptions.add("pfRecoTauDiscriminationByIsolationMVA2", desc);
0256 }
0257
0258 DEFINE_FWK_MODULE(PFRecoTauDiscriminationByIsolationMVA2);