Back to home page

Project CMSSW displayed by LXR

 
 

    


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_;  // 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 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 //,vetos_(new AbsVetos())
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         //"Supported modes are 'sum', 'sumRelative', 'max', 'maxRelative', 'count'." << // TODO: on request only
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     //std::cout << "Weight is a simple number, " << toNumber(weight) << std::endl;
0113     weight_ = toNumber(weight);
0114     usesFunction_ = false;
0115   } else {
0116     usesFunction_ = true;
0117     //std::cout << "Weight is a function, this might slow you down... " << std::endl;
0118   }
0119   //std::cout << "CandIsolatorFromDeposits::SingleDeposit::SingleDeposit: Total of " << vetos_.size() << " vetos" << std::endl;
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   // NOTE: we DON'T have to delete the evdepVetos_, they have already been deleted above. We just clear the vectors
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();  // better to center on the deposit direction
0139                                             // that could be, e.g., the impact point at calo
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__;  // avoid gcc warning
0168 }
0169 
0170 /// constructor with config
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 /// destructor
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 /// build deposits
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);