File indexing completed on 2022-04-29 23:10:54
0001 #include "CommonTools/Utils/interface/StringObjectFunction.h"
0002 #include "DataFormats/Candidate/interface/CandAssociation.h"
0003 #include "DataFormats/Candidate/interface/Candidate.h"
0004 #include "DataFormats/Common/interface/Handle.h"
0005 #include "DataFormats/MuonReco/interface/Muon.h"
0006 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
0007 #include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
0008 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
0009 #include "DataFormats/RecoCandidate/interface/IsoDepositVetos.h"
0010 #include "DataFormats/TrackReco/interface/Track.h"
0011 #include "FWCore/Framework/interface/ConsumesCollector.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Framework/interface/stream/EDProducer.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "PhysicsTools/IsolationAlgos/interface/EventDependentAbsVeto.h"
0018 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositVetoFactory.h"
0019
0020 #include <string>
0021 #include <regex>
0022
0023 class CandIsolatorFromDeposits : public edm::stream::EDProducer<> {
0024 public:
0025 typedef edm::ValueMap<double> CandDoubleMap;
0026
0027 enum Mode { Sum, SumRelative, Sum2, Sum2Relative, Max, MaxRelative, Count, NearestDR, MeanDR, SumDR };
0028 CandIsolatorFromDeposits(const edm::ParameterSet &);
0029
0030 ~CandIsolatorFromDeposits() override;
0031
0032 void produce(edm::Event &, const edm::EventSetup &) override;
0033
0034 private:
0035 class SingleDeposit {
0036 public:
0037 SingleDeposit(const edm::ParameterSet &, edm::ConsumesCollector &&iC);
0038 void cleanup();
0039 void open(const edm::Event &iEvent, const edm::EventSetup &iSetup);
0040 double compute(const reco::CandidateBaseRef &cand);
0041 const reco::IsoDepositMap &map() { return *hDeps_; }
0042
0043 private:
0044 Mode mode_;
0045 edm::EDGetTokenT<reco::IsoDepositMap> srcToken_;
0046 double deltaR_;
0047 bool usesFunction_;
0048 double weight_;
0049 StringObjectFunction<reco::Candidate> weightExpr_;
0050 reco::isodeposit::AbsVetos vetos_;
0051 reco::isodeposit::EventDependentAbsVetos evdepVetos_;
0052 bool skipDefaultVeto_;
0053 edm::Handle<reco::IsoDepositMap> hDeps_;
0054 };
0055
0056 std::vector<SingleDeposit> sources_;
0057 };
0058
0059 using namespace edm;
0060 using namespace reco;
0061 using namespace reco::isodeposit;
0062
0063 bool isNumber(const std::string &str) {
0064 static const std::regex re("^[+-]?(\\d+\\.?|\\d*\\.\\d*)$");
0065 return regex_match(str.c_str(), re);
0066 }
0067 double toNumber(const std::string &str) { return atof(str.c_str()); }
0068
0069 CandIsolatorFromDeposits::SingleDeposit::SingleDeposit(const edm::ParameterSet &iConfig, edm::ConsumesCollector &&iC)
0070 : srcToken_(iC.consumes<reco::IsoDepositMap>(iConfig.getParameter<edm::InputTag>("src"))),
0071 deltaR_(iConfig.getParameter<double>("deltaR")),
0072 weightExpr_(iConfig.getParameter<std::string>("weight")),
0073 skipDefaultVeto_(iConfig.getParameter<bool>("skipDefaultVeto"))
0074
0075 {
0076 std::string mode = iConfig.getParameter<std::string>("mode");
0077 if (mode == "sum")
0078 mode_ = Sum;
0079 else if (mode == "sumRelative")
0080 mode_ = SumRelative;
0081 else if (mode == "sum2")
0082 mode_ = Sum2;
0083 else if (mode == "sum2Relative")
0084 mode_ = Sum2Relative;
0085 else if (mode == "max")
0086 mode_ = Max;
0087 else if (mode == "maxRelative")
0088 mode_ = MaxRelative;
0089 else if (mode == "nearestDR")
0090 mode_ = NearestDR;
0091 else if (mode == "count")
0092 mode_ = Count;
0093 else if (mode == "meanDR")
0094 mode_ = MeanDR;
0095 else if (mode == "sumDR")
0096 mode_ = SumDR;
0097 else
0098 throw cms::Exception("Not Implemented") << "Mode '" << mode << "' not implemented. "
0099 << "Supported modes are 'sum', 'sumRelative', 'count'." <<
0100
0101 "New methods can be easily implemented if requested.";
0102 typedef std::vector<std::string> vstring;
0103 vstring vetos = iConfig.getParameter<vstring>("vetos");
0104 reco::isodeposit::EventDependentAbsVeto *evdep;
0105 for (vstring::const_iterator it = vetos.begin(), ed = vetos.end(); it != ed; ++it) {
0106 vetos_.push_back(IsoDepositVetoFactory::make(it->c_str(), evdep, iC));
0107 if (evdep)
0108 evdepVetos_.push_back(evdep);
0109 }
0110 std::string weight = iConfig.getParameter<std::string>("weight");
0111 if (isNumber(weight)) {
0112
0113 weight_ = toNumber(weight);
0114 usesFunction_ = false;
0115 } else {
0116 usesFunction_ = true;
0117
0118 }
0119
0120 }
0121 void CandIsolatorFromDeposits::SingleDeposit::cleanup() {
0122 for (AbsVetos::iterator it = vetos_.begin(), ed = vetos_.end(); it != ed; ++it) {
0123 delete *it;
0124 }
0125 vetos_.clear();
0126
0127 evdepVetos_.clear();
0128 }
0129 void CandIsolatorFromDeposits::SingleDeposit::open(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0130 iEvent.getByToken(srcToken_, hDeps_);
0131 for (EventDependentAbsVetos::iterator it = evdepVetos_.begin(), ed = evdepVetos_.end(); it != ed; ++it) {
0132 (*it)->setEvent(iEvent, iSetup);
0133 }
0134 }
0135
0136 double CandIsolatorFromDeposits::SingleDeposit::compute(const reco::CandidateBaseRef &cand) {
0137 const IsoDeposit &dep = (*hDeps_)[cand];
0138 double eta = dep.eta(), phi = dep.phi();
0139
0140 for (AbsVetos::iterator it = vetos_.begin(), ed = vetos_.end(); it != ed; ++it) {
0141 (*it)->centerOn(eta, phi);
0142 }
0143 double weight = (usesFunction_ ? weightExpr_(*cand) : weight_);
0144 switch (mode_) {
0145 case Count:
0146 return weight * dep.countWithin(deltaR_, vetos_, skipDefaultVeto_);
0147 case Sum:
0148 return weight * dep.sumWithin(deltaR_, vetos_, skipDefaultVeto_);
0149 case SumRelative:
0150 return weight * dep.sumWithin(deltaR_, vetos_, skipDefaultVeto_) / dep.candEnergy();
0151 case Sum2:
0152 return weight * dep.sum2Within(deltaR_, vetos_, skipDefaultVeto_);
0153 case Sum2Relative:
0154 return weight * dep.sum2Within(deltaR_, vetos_, skipDefaultVeto_) / (dep.candEnergy() * dep.candEnergy());
0155 case Max:
0156 return weight * dep.maxWithin(deltaR_, vetos_, skipDefaultVeto_);
0157 case NearestDR:
0158 return weight * dep.nearestDR(deltaR_, vetos_, skipDefaultVeto_);
0159 case MaxRelative:
0160 return weight * dep.maxWithin(deltaR_, vetos_, skipDefaultVeto_) / dep.candEnergy();
0161 case MeanDR:
0162 return weight * dep.algoWithin<reco::IsoDeposit::MeanDRAlgo>(deltaR_, vetos_, skipDefaultVeto_);
0163 case SumDR:
0164 return weight * dep.algoWithin<reco::IsoDeposit::SumDRAlgo>(deltaR_, vetos_, skipDefaultVeto_);
0165 }
0166 throw cms::Exception("Logic error") << "Should not happen at " << __FILE__ << ", line "
0167 << __LINE__;
0168 }
0169
0170
0171 CandIsolatorFromDeposits::CandIsolatorFromDeposits(const ParameterSet &par) {
0172 typedef std::vector<edm::ParameterSet> VPSet;
0173 VPSet depPSets = par.getParameter<VPSet>("deposits");
0174 for (VPSet::const_iterator it = depPSets.begin(), ed = depPSets.end(); it != ed; ++it) {
0175 sources_.push_back(SingleDeposit(*it, consumesCollector()));
0176 }
0177 if (sources_.empty())
0178 throw cms::Exception("Configuration Error") << "Please specify at least one deposit!";
0179 produces<CandDoubleMap>();
0180 }
0181
0182
0183 CandIsolatorFromDeposits::~CandIsolatorFromDeposits() {
0184 std::vector<SingleDeposit>::iterator it, begin = sources_.begin(), end = sources_.end();
0185 for (it = begin; it != end; ++it)
0186 it->cleanup();
0187 }
0188
0189
0190 void CandIsolatorFromDeposits::produce(Event &event, const EventSetup &eventSetup) {
0191 std::vector<SingleDeposit>::iterator it, begin = sources_.begin(), end = sources_.end();
0192 for (it = begin; it != end; ++it)
0193 it->open(event, eventSetup);
0194
0195 const IsoDepositMap &map = begin->map();
0196
0197 if (map.empty()) {
0198 event.put(std::make_unique<CandDoubleMap>());
0199 return;
0200 }
0201 auto ret = std::make_unique<CandDoubleMap>();
0202 CandDoubleMap::Filler filler(*ret);
0203
0204 typedef reco::IsoDepositMap::const_iterator iterator_i;
0205 typedef reco::IsoDepositMap::container::const_iterator iterator_ii;
0206 iterator_i depI = map.begin();
0207 iterator_i depIEnd = map.end();
0208 for (; depI != depIEnd; ++depI) {
0209 std::vector<double> retV(depI.size(), 0);
0210 edm::Handle<edm::View<reco::Candidate> > candH;
0211 event.get(depI.id(), candH);
0212 const edm::View<reco::Candidate> &candV = *candH;
0213
0214 iterator_ii depII = depI.begin();
0215 iterator_ii depIIEnd = depI.end();
0216 size_t iRet = 0;
0217 for (; depII != depIIEnd; ++depII, ++iRet) {
0218 double sum = 0;
0219 for (it = begin; it != end; ++it)
0220 sum += it->compute(candV.refAt(iRet));
0221 retV[iRet] = sum;
0222 }
0223 filler.insert(candH, retV.begin(), retV.end());
0224 }
0225 filler.fill();
0226 event.put(std::move(ret));
0227 }
0228
0229 #include "FWCore/Framework/interface/MakerMacros.h"
0230 DEFINE_FWK_MODULE(CandIsolatorFromDeposits);