plugins/cuts/GsfEleMVACut.cc

0001 #include "PhysicsTools/SelectorUtils/interface/CutApplicatorWithEventContentBase.h"
0002 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0003 #include "CommonTools/Utils/interface/StringObjectFunction.h"
0004
0005 class GsfEleMVACut : public CutApplicatorWithEventContentBase {
0006 public:
0007   GsfEleMVACut(const edm::ParameterSet& c);
0008
0009   result_type operator()(const reco::GsfElectronPtr&) const final;
0010
0011   void setConsumes(edm::ConsumesCollector&) final;
0012   void getEventContent(const edm::EventBase&) final;
0013
0014   CandidateType candidateType() const final { return ELECTRON; }
0015
0016 private:
0017   double value(const reco::CandidatePtr& cand) const final;
0018
0019   // Cut formulas
0020   const std::vector<std::string> mvaCutStrings_;
0021   std::vector<StringObjectFunction<reco::GsfElectron>> cutFormula_;
0022
0023   const int nCuts_;
0024
0025   // Pre-computed MVA value map
0026   edm::Handle<edm::ValueMap<float>> mvaValueMap_;
0027   edm::Handle<edm::ValueMap<int>> mvaCategoriesMap_;
0028 };
0029
0030 DEFINE_EDM_PLUGIN(CutApplicatorFactory, GsfEleMVACut, "GsfEleMVACut");
0031
0032 GsfEleMVACut::GsfEleMVACut(const edm::ParameterSet& c)
0033     : CutApplicatorWithEventContentBase(c),
0034       mvaCutStrings_(c.getParameter<std::vector<std::string>>("mvaCuts")),
0035       nCuts_(mvaCutStrings_.size()) {
0036   edm::InputTag mvaValTag = c.getParameter<edm::InputTag>("mvaValueMapName");
0037   contentTags_.emplace("mvaVal", mvaValTag);
0038
0039   edm::InputTag mvaCatTag = c.getParameter<edm::InputTag>("mvaCategoriesMapName");
0040   contentTags_.emplace("mvaCat", mvaCatTag);
0041
0042   for (auto& cutString : mvaCutStrings_) {
0043     cutFormula_.push_back(StringObjectFunction<reco::GsfElectron>(cutString));
0044   }
0045 }
0046
0047 void GsfEleMVACut::setConsumes(edm::ConsumesCollector& cc) {
0048   auto mvaVal = cc.consumes<edm::ValueMap<float>>(contentTags_["mvaVal"]);
0049   contentTokens_.emplace("mvaVal", mvaVal);
0050
0051   auto mvaCat = cc.consumes<edm::ValueMap<int>>(contentTags_["mvaCat"]);
0052   contentTokens_.emplace("mvaCat", mvaCat);
0053 }
0054
0055 void GsfEleMVACut::getEventContent(const edm::EventBase& ev) {
0056   ev.getByLabel(contentTags_["mvaVal"], mvaValueMap_);
0057   ev.getByLabel(contentTags_["mvaCat"], mvaCategoriesMap_);
0058 }
0059
0060 CutApplicatorBase::result_type GsfEleMVACut::operator()(const reco::GsfElectronPtr& cand) const {
0061   // in case we are by-value
0062   const std::string& val_name = contentTags_.find("mvaVal")->second.instance();
0063   const std::string& cat_name = contentTags_.find("mvaCat")->second.instance();
0064   edm::Ptr<pat::Electron> pat(cand);
0065   float val = -1.0;
0066   int cat = -1;
0067   if (mvaCategoriesMap_.isValid() && mvaCategoriesMap_->contains(cand.id()) && mvaValueMap_.isValid() &&
0068       mvaValueMap_->contains(cand.id())) {
0069     cat = (*mvaCategoriesMap_)[cand];
0070     val = (*mvaValueMap_)[cand];
0071   } else if (mvaCategoriesMap_.isValid() && mvaValueMap_.isValid() && mvaCategoriesMap_->idSize() == 1 &&
0072              mvaValueMap_->idSize() == 1 && cand.id() == edm::ProductID()) {
0073     // in case we have spoofed a ptr
0074     //note this must be a 1:1 valuemap (only one product input)
0075     cat = mvaCategoriesMap_->begin()[cand.key()];
0076     val = mvaValueMap_->begin()[cand.key()];
0077   } else if (mvaCategoriesMap_.isValid() && mvaValueMap_.isValid()) {  // throw an exception
0078     cat = (*mvaCategoriesMap_)[cand];
0079     val = (*mvaValueMap_)[cand];
0080   }
0081
0082   // Find the cut formula
0083   const int iCategory = mvaCategoriesMap_.isValid() ? cat : pat->userInt(cat_name);
0084   if (iCategory >= nCuts_)
0085     throw cms::Exception(" Error in MVA categories: ")
0086         << " found a particle with a category larger than max configured " << std::endl;
0087
0088   // Look up the MVA value for this particle
0089   const float mvaValue = mvaValueMap_.isValid() ? val : pat->userFloat(val_name);
0090
0091   // Apply the cut and return the result
0092   return mvaValue > cutFormula_[iCategory](*cand);
0093 }
0094
0095 double GsfEleMVACut::value(const reco::CandidatePtr& cand) const {
0096   // in case we are by-value
0097   const std::string& val_name = contentTags_.find("mvaVal")->second.instance();
0098   edm::Ptr<pat::Electron> pat(cand);
0099   float val = 0.0;
0100   if (mvaCategoriesMap_.isValid() && mvaCategoriesMap_->contains(cand.id()) && mvaValueMap_.isValid() &&
0101       mvaValueMap_->contains(cand.id())) {
0102     val = (*mvaValueMap_)[cand];
0103   } else if (mvaCategoriesMap_.isValid() && mvaValueMap_.isValid() && mvaCategoriesMap_->idSize() == 1 &&
0104              mvaValueMap_->idSize() == 1 && cand.id() == edm::ProductID()) {
0105     // in case we have spoofed a ptr
0106     //note this must be a 1:1 valuemap (only one product input)
0107     val = mvaValueMap_->begin()[cand.key()];
0108   } else if (mvaCategoriesMap_.isValid() && mvaValueMap_.isValid()) {  // throw an exception
0109     val = (*mvaValueMap_)[cand];
0110   }
0111
0112   const float mvaValue = mvaValueMap_.isValid() ? val : pat->userFloat(val_name);
0113   return mvaValue;
0114 }