EgammaTools/interface/MVAValueMapProducer.h

0001 #ifndef __RecoEgamma_EgammaTools_MVAValueMapProducer_H__
0002 #define __RecoEgamma_EgammaTools_MVAValueMapProducer_H__
0003
0004 #include "DataFormats/Common/interface/ValueMap.h"
0005 #include "DataFormats/Common/interface/View.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/global/EDProducer.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "RecoEgamma/EgammaTools/interface/AnyMVAEstimatorRun2Base.h"
0012 #include "RecoEgamma/EgammaTools/interface/AnyMVAEstimatorRun2Factory.h"
0013 #include "FWCore/Utilities/interface/EDGetToken.h"
0014 #include "RecoEgamma/EgammaTools/interface/MVAVariableHelper.h"
0015 #include "DataFormats/Common/interface/Handle.h"
0016 #include "DataFormats/Common/interface/ValueMap.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "RecoEgamma/EgammaTools/interface/validateEgammaCandidate.h"
0019 #include "FWCore/Framework/interface/ConsumesCollector.h"
0020 #include "FWCore/Utilities/interface/thread_safety_macros.h"
0021
0022 #include <atomic>
0023 #include <cmath>
0024 #include <memory>
0025 #include <string>
0026 #include <vector>
0027
0028 template <class ParticleType>
0029 class MVAValueMapProducer : public edm::global::EDProducer<> {
0030 public:
0031   MVAValueMapProducer(const edm::ParameterSet&);
0032
0033   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0034
0035 private:
0036   static auto getMVAEstimators(const edm::VParameterSet& vConfig) {
0037     std::vector<std::unique_ptr<AnyMVAEstimatorRun2Base>> mvaEstimators;
0038
0039     // Loop over the list of MVA configurations passed here from python and
0040     // construct all requested MVA estimators.
0041     for (auto& imva : vConfig) {
0042       // The factory below constructs the MVA of the appropriate type based
0043       // on the "mvaName" which is the name of the derived MVA class (plugin)
0044       if (!imva.empty()) {
0045         mvaEstimators.emplace_back(
0046             AnyMVAEstimatorRun2Factory::get()->create(imva.getParameter<std::string>("mvaName"), imva));
0047
0048       } else
0049         throw cms::Exception(" MVA configuration not found: ")
0050             << " failed to find proper configuration for one of the MVAs in the main python script " << std::endl;
0051     }
0052
0053     return mvaEstimators;
0054   }
0055
0056   static std::vector<std::string> getValueMapNames(const edm::VParameterSet& vConfig, std::string&& suffix) {
0057     std::vector<std::string> names;
0058     for (auto& imva : vConfig) {
0059       names.push_back(imva.getParameter<std::string>("mvaName") + imva.getParameter<std::string>("mvaTag") + suffix);
0060     }
0061
0062     return names;
0063   }
0064
0065   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0066
0067   const edm::EDGetTokenT<edm::View<ParticleType>> srcToken_;
0068   const edm::EDGetTokenT<edm::View<ParticleType>> keysForValueMapsToken_;
0069
0070   // MVA estimators
0071   const std::vector<std::unique_ptr<AnyMVAEstimatorRun2Base>> mvaEstimators_;
0072
0073   // Value map names
0074   const std::vector<std::string> mvaValueMapNames_;
0075   const std::vector<std::string> mvaRawValueMapNames_;
0076   const std::vector<std::string> mvaCategoriesMapNames_;
0077
0078   // To get the auxiliary MVA variables
0079   const MVAVariableHelper variableHelper_;
0080
0081   CMS_THREAD_SAFE mutable std::atomic<bool> validated_ = false;
0082 };
0083
0084 namespace {
0085
0086   template <typename ValueType, class HandleType>
0087   void writeValueMap(edm::Event& iEvent,
0088                      const edm::Handle<HandleType>& handle,
0089                      const std::vector<ValueType>& values,
0090                      const std::string& label) {
0091     auto valMap = std::make_unique<edm::ValueMap<ValueType>>();
0092     typename edm::ValueMap<ValueType>::Filler filler(*valMap);
0093     filler.insert(handle, values.begin(), values.end());
0094     filler.fill();
0095     iEvent.put(std::move(valMap), label);
0096   }
0097
0098   template <class ParticleType>
0099   auto getKeysForValueMapsToken(edm::InputTag const& keysForValueMapsTag, edm::ConsumesCollector&& cc) {
0100     const bool tagGiven = !keysForValueMapsTag.label().empty();
0101     return tagGiven ? cc.consumes<edm::View<ParticleType>>(keysForValueMapsTag)
0102                     : edm::EDGetTokenT<edm::View<ParticleType>>{};
0103   }
0104
0105 }  // namespace
0106
0107 template <class ParticleType>
0108 MVAValueMapProducer<ParticleType>::MVAValueMapProducer(const edm::ParameterSet& iConfig)
0109     : srcToken_(consumes<edm::View<ParticleType>>(iConfig.getParameter<edm::InputTag>("src"))),
0110       keysForValueMapsToken_(getKeysForValueMapsToken<ParticleType>(
0111           iConfig.getParameter<edm::InputTag>("keysForValueMaps"), consumesCollector())),
0112       mvaEstimators_(getMVAEstimators(iConfig.getParameterSetVector("mvaConfigurations"))),
0113       mvaValueMapNames_(getValueMapNames(iConfig.getParameterSetVector("mvaConfigurations"), "Values")),
0114       mvaRawValueMapNames_(getValueMapNames(iConfig.getParameterSetVector("mvaConfigurations"), "RawValues")),
0115       mvaCategoriesMapNames_(getValueMapNames(iConfig.getParameterSetVector("mvaConfigurations"), "Categories")),
0116       variableHelper_(consumesCollector()) {
0117   for (auto const& name : mvaValueMapNames_)
0118     produces<edm::ValueMap<float>>(name);
0119   for (auto const& name : mvaRawValueMapNames_)
0120     produces<edm::ValueMap<float>>(name);
0121   for (auto const& name : mvaCategoriesMapNames_)
0122     produces<edm::ValueMap<int>>(name);
0123 }
0124
0125 template <class ParticleType>
0126 void MVAValueMapProducer<ParticleType>::produce(edm::StreamID,
0127                                                 edm::Event& iEvent,
0128                                                 const edm::EventSetup& iSetup) const {
0129   std::vector<float> auxVariables = variableHelper_.getAuxVariables(iEvent);
0130
0131   auto srcHandle = iEvent.getHandle(srcToken_);
0132   auto keysForValueMapsHandle =
0133       keysForValueMapsToken_.isUninitialized() ? srcHandle : iEvent.getHandle(keysForValueMapsToken_);
0134
0135   // check if nothing is wrong with the data format of the candidates
0136   if (!validated_ && !srcHandle->empty()) {
0137     egammaTools::validateEgammaCandidate((*srcHandle)[0]);
0138     validated_ = true;
0139   }
0140
0141   // Loop over MVA estimators
0142   for (unsigned iEstimator = 0; iEstimator < mvaEstimators_.size(); iEstimator++) {
0143     std::vector<float> mvaValues;
0144     std::vector<float> mvaRawValues;
0145     std::vector<int> mvaCategories;
0146
0147     // Loop over particles
0148     for (auto const& cand : *srcHandle) {
0149       int cat = -1;  // Passed by reference to the mvaValue function to store the category
0150       const float response = mvaEstimators_[iEstimator]->mvaValue(&cand, auxVariables, cat);
0151       mvaRawValues.push_back(response);                             // The MVA score
0152       mvaValues.push_back(2.0 / (1.0 + exp(-2.0 * response)) - 1);  // MVA output between -1 and 1
0153       mvaCategories.push_back(cat);
0154     }  // end loop over particles
0155
0156     writeValueMap(iEvent, keysForValueMapsHandle, mvaValues, mvaValueMapNames_[iEstimator]);
0157     writeValueMap(iEvent, keysForValueMapsHandle, mvaRawValues, mvaRawValueMapNames_[iEstimator]);
0158     writeValueMap(iEvent, keysForValueMapsHandle, mvaCategories, mvaCategoriesMapNames_[iEstimator]);
0159
0160   }  // end loop over estimators
0161 }
0162
0163 template <class ParticleType>
0164 void MVAValueMapProducer<ParticleType>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0165   edm::ParameterSetDescription desc;
0166   desc.add<edm::InputTag>("src", {});
0167   desc.add<edm::InputTag>("keysForValueMaps", {});
0168   {
0169     //The following says we do not know what parameters are allowed so do no validation
0170     // Please change this to state exactly what you do use, even if it is no parameters
0171     edm::ParameterSetDescription mvaConfigurations;
0172     mvaConfigurations.setUnknown();
0173     desc.addVPSet("mvaConfigurations", mvaConfigurations);
0174   }
0175   descriptions.addDefault(desc);
0176 }
0177
0178 #endif