File indexing completed on 2023-03-17 11:15:51
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 namespace {
0064 bool isNumber(const std::string &str) {
0065 static const std::regex re("^[+-]?(\\d+\\.?|\\d*\\.\\d*)$");
0066 return regex_match(str.c_str(), re);
0067 }
0068 double toNumber(const std::string &str) { return atof(str.c_str()); }
0069 }
0070
0071 CandIsolatorFromDeposits::SingleDeposit::SingleDeposit(const edm::ParameterSet &iConfig, edm::ConsumesCollector &&iC)
0072 : srcToken_(iC.consumes<reco::IsoDepositMap>(iConfig.getParameter<edm::InputTag>("src"))),
0073 deltaR_(iConfig.getParameter<double>("deltaR")),
0074 weightExpr_(iConfig.getParameter<std::string>("weight")),
0075 skipDefaultVeto_(iConfig.getParameter<bool>("skipDefaultVeto"))
0076
0077 {
0078 std::string mode = iConfig.getParameter<std::string>("mode");
0079 if (mode == "sum")
0080 mode_ = Sum;
0081 else if (mode == "sumRelative")
0082 mode_ = SumRelative;
0083 else if (mode == "sum2")
0084 mode_ = Sum2;
0085 else if (mode == "sum2Relative")
0086 mode_ = Sum2Relative;
0087 else if (mode == "max")
0088 mode_ = Max;
0089 else if (mode == "maxRelative")
0090 mode_ = MaxRelative;
0091 else if (mode == "nearestDR")
0092 mode_ = NearestDR;
0093 else if (mode == "count")
0094 mode_ = Count;
0095 else if (mode == "meanDR")
0096 mode_ = MeanDR;
0097 else if (mode == "sumDR")
0098 mode_ = SumDR;
0099 else
0100 throw cms::Exception("Not Implemented") << "Mode '" << mode << "' not implemented. "
0101 << "Supported modes are 'sum', 'sumRelative', 'count'." <<
0102
0103 "New methods can be easily implemented if requested.";
0104 typedef std::vector<std::string> vstring;
0105 vstring vetos = iConfig.getParameter<vstring>("vetos");
0106 reco::isodeposit::EventDependentAbsVeto *evdep;
0107 for (vstring::const_iterator it = vetos.begin(), ed = vetos.end(); it != ed; ++it) {
0108 vetos_.push_back(IsoDepositVetoFactory::make(it->c_str(), evdep, iC));
0109 if (evdep)
0110 evdepVetos_.push_back(evdep);
0111 }
0112 std::string weight = iConfig.getParameter<std::string>("weight");
0113 if (isNumber(weight)) {
0114
0115 weight_ = toNumber(weight);
0116 usesFunction_ = false;
0117 } else {
0118 usesFunction_ = true;
0119
0120 }
0121
0122 }
0123 void CandIsolatorFromDeposits::SingleDeposit::cleanup() {
0124 for (AbsVetos::iterator it = vetos_.begin(), ed = vetos_.end(); it != ed; ++it) {
0125 delete *it;
0126 }
0127 vetos_.clear();
0128
0129 evdepVetos_.clear();
0130 }
0131 void CandIsolatorFromDeposits::SingleDeposit::open(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0132 iEvent.getByToken(srcToken_, hDeps_);
0133 for (EventDependentAbsVetos::iterator it = evdepVetos_.begin(), ed = evdepVetos_.end(); it != ed; ++it) {
0134 (*it)->setEvent(iEvent, iSetup);
0135 }
0136 }
0137
0138 double CandIsolatorFromDeposits::SingleDeposit::compute(const reco::CandidateBaseRef &cand) {
0139 const IsoDeposit &dep = (*hDeps_)[cand];
0140 double eta = dep.eta(), phi = dep.phi();
0141
0142 for (AbsVetos::iterator it = vetos_.begin(), ed = vetos_.end(); it != ed; ++it) {
0143 (*it)->centerOn(eta, phi);
0144 }
0145 double weight = (usesFunction_ ? weightExpr_(*cand) : weight_);
0146 switch (mode_) {
0147 case Count:
0148 return weight * dep.countWithin(deltaR_, vetos_, skipDefaultVeto_);
0149 case Sum:
0150 return weight * dep.sumWithin(deltaR_, vetos_, skipDefaultVeto_);
0151 case SumRelative:
0152 return weight * dep.sumWithin(deltaR_, vetos_, skipDefaultVeto_) / dep.candEnergy();
0153 case Sum2:
0154 return weight * dep.sum2Within(deltaR_, vetos_, skipDefaultVeto_);
0155 case Sum2Relative:
0156 return weight * dep.sum2Within(deltaR_, vetos_, skipDefaultVeto_) / (dep.candEnergy() * dep.candEnergy());
0157 case Max:
0158 return weight * dep.maxWithin(deltaR_, vetos_, skipDefaultVeto_);
0159 case NearestDR:
0160 return weight * dep.nearestDR(deltaR_, vetos_, skipDefaultVeto_);
0161 case MaxRelative:
0162 return weight * dep.maxWithin(deltaR_, vetos_, skipDefaultVeto_) / dep.candEnergy();
0163 case MeanDR:
0164 return weight * dep.algoWithin<reco::IsoDeposit::MeanDRAlgo>(deltaR_, vetos_, skipDefaultVeto_);
0165 case SumDR:
0166 return weight * dep.algoWithin<reco::IsoDeposit::SumDRAlgo>(deltaR_, vetos_, skipDefaultVeto_);
0167 }
0168 throw cms::Exception("Logic error") << "Should not happen at " << __FILE__ << ", line "
0169 << __LINE__;
0170 }
0171
0172
0173 CandIsolatorFromDeposits::CandIsolatorFromDeposits(const ParameterSet &par) {
0174 typedef std::vector<edm::ParameterSet> VPSet;
0175 VPSet depPSets = par.getParameter<VPSet>("deposits");
0176 for (VPSet::const_iterator it = depPSets.begin(), ed = depPSets.end(); it != ed; ++it) {
0177 sources_.push_back(SingleDeposit(*it, consumesCollector()));
0178 }
0179 if (sources_.empty())
0180 throw cms::Exception("Configuration Error") << "Please specify at least one deposit!";
0181 produces<CandDoubleMap>();
0182 }
0183
0184
0185 CandIsolatorFromDeposits::~CandIsolatorFromDeposits() {
0186 std::vector<SingleDeposit>::iterator it, begin = sources_.begin(), end = sources_.end();
0187 for (it = begin; it != end; ++it)
0188 it->cleanup();
0189 }
0190
0191
0192 void CandIsolatorFromDeposits::produce(Event &event, const EventSetup &eventSetup) {
0193 std::vector<SingleDeposit>::iterator it, begin = sources_.begin(), end = sources_.end();
0194 for (it = begin; it != end; ++it)
0195 it->open(event, eventSetup);
0196
0197 const IsoDepositMap &map = begin->map();
0198
0199 if (map.empty()) {
0200 event.put(std::make_unique<CandDoubleMap>());
0201 return;
0202 }
0203 auto ret = std::make_unique<CandDoubleMap>();
0204 CandDoubleMap::Filler filler(*ret);
0205
0206 typedef reco::IsoDepositMap::const_iterator iterator_i;
0207 typedef reco::IsoDepositMap::container::const_iterator iterator_ii;
0208 iterator_i depI = map.begin();
0209 iterator_i depIEnd = map.end();
0210 for (; depI != depIEnd; ++depI) {
0211 std::vector<double> retV(depI.size(), 0);
0212 edm::Handle<edm::View<reco::Candidate> > candH;
0213 event.get(depI.id(), candH);
0214 const edm::View<reco::Candidate> &candV = *candH;
0215
0216 iterator_ii depII = depI.begin();
0217 iterator_ii depIIEnd = depI.end();
0218 size_t iRet = 0;
0219 for (; depII != depIIEnd; ++depII, ++iRet) {
0220 double sum = 0;
0221 for (it = begin; it != end; ++it)
0222 sum += it->compute(candV.refAt(iRet));
0223 retV[iRet] = sum;
0224 }
0225 filler.insert(candH, retV.begin(), retV.end());
0226 }
0227 filler.fill();
0228 event.put(std::move(ret));
0229 }
0230
0231 #include "FWCore/Framework/interface/MakerMacros.h"
0232 DEFINE_FWK_MODULE(CandIsolatorFromDeposits);