IsolationAlgos/plugins/CandIsolatorFromDeposits.cc

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_;  // note: these are a subset of the above. Don't delete twice!
0052     bool skipDefaultVeto_;
0053     edm::Handle<reco::IsoDepositMap> hDeps_;  // transient
0054   };
0055   // datamembers
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 }  // namespace
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 //,vetos_(new AbsVetos())
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         //"Supported modes are 'sum', 'sumRelative', 'max', 'maxRelative', 'count'." << // TODO: on request only
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     //std::cout << "Weight is a simple number, " << toNumber(weight) << std::endl;
0115     weight_ = toNumber(weight);
0116     usesFunction_ = false;
0117   } else {
0118     usesFunction_ = true;
0119     //std::cout << "Weight is a function, this might slow you down... " << std::endl;
0120   }
0121   //std::cout << "CandIsolatorFromDeposits::SingleDeposit::SingleDeposit: Total of " << vetos_.size() << " vetos" << std::endl;
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   // NOTE: we DON'T have to delete the evdepVetos_, they have already been deleted above. We just clear the vectors
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();  // better to center on the deposit direction
0141                                             // that could be, e.g., the impact point at calo
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__;  // avoid gcc warning
0170 }
0171
0172 /// constructor with config
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 /// destructor
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 /// build deposits
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);