ParticleFlow/plugins/PFCandIsolatorFromDeposit.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/Common/interface/OwnVector.h"
0006 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0007 #include "DataFormats/MuonReco/interface/Muon.h"
0008 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0009 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0010 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
0011 #include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
0012 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
0013 #include "DataFormats/RecoCandidate/interface/IsoDepositVetos.h"
0014 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0015 #include "DataFormats/TrackReco/interface/Track.h"
0016 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0017 #include "FWCore/Framework/interface/ConsumesCollector.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Framework/interface/MakerMacros.h"
0021 #include "FWCore/Framework/interface/stream/EDProducer.h"
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "PhysicsTools/IsolationAlgos/interface/EventDependentAbsVeto.h"
0025 #include "PhysicsTools/IsolationAlgos/interface/IsoDepositVetoFactory.h"
0026
0027 #include <memory>
0028 #include <regex>
0029 #include <string>
0030
0031 class PFCandIsolatorFromDeposits : public edm::stream::EDProducer<> {
0032 public:
0033   typedef edm::ValueMap<double> CandDoubleMap;
0034
0035   enum Mode { Sum, SumRelative, Sum2, Sum2Relative, Max, MaxRelative, Count, NearestDR };
0036   PFCandIsolatorFromDeposits(const edm::ParameterSet &);
0037
0038   ~PFCandIsolatorFromDeposits() override;
0039
0040   void produce(edm::Event &, const edm::EventSetup &) override;
0041
0042 private:
0043   class SingleDeposit {
0044   public:
0045     SingleDeposit(const edm::ParameterSet &, edm::ConsumesCollector &&iC);
0046     void cleanup();
0047     void open(const edm::Event &iEvent, const edm::EventSetup &iSetup);
0048     double compute(const reco::CandidateBaseRef &cand);
0049     const reco::IsoDepositMap &map() { return *hDeps_; }
0050
0051   private:
0052     Mode mode_;
0053     edm::EDGetTokenT<reco::IsoDepositMap> srcToken_;
0054     double deltaR_;
0055     bool usesFunction_;
0056     double weight_;
0057
0058     StringObjectFunction<reco::Candidate> weightExpr_;
0059     reco::isodeposit::AbsVetos barrelVetos_;
0060     reco::isodeposit::AbsVetos endcapVetos_;
0061     reco::isodeposit::EventDependentAbsVetos evdepVetos_;  // note: these are a subset of the above. Don't delete twice!
0062     bool skipDefaultVeto_;
0063     bool usePivotForBarrelEndcaps_;
0064     edm::Handle<reco::IsoDepositMap> hDeps_;  // transient
0065
0066     bool isNumber(const std::string &str) const;
0067     double toNumber(const std::string &str) const;
0068   };
0069   // datamembers
0070   std::vector<SingleDeposit> sources_;
0071 };
0072
0073 using namespace edm;
0074 using namespace reco;
0075 using namespace reco::isodeposit;
0076
0077 bool PFCandIsolatorFromDeposits::SingleDeposit::isNumber(const std::string &str) const {
0078   static const std::regex re("^[+-]?(\\d+\\.?|\\d*\\.\\d*)$");
0079   return regex_match(str.c_str(), re);
0080 }
0081 double PFCandIsolatorFromDeposits::SingleDeposit::toNumber(const std::string &str) const { return atof(str.c_str()); }
0082
0083 PFCandIsolatorFromDeposits::SingleDeposit::SingleDeposit(const edm::ParameterSet &iConfig, edm::ConsumesCollector &&iC)
0084     : srcToken_(iC.consumes<reco::IsoDepositMap>(iConfig.getParameter<edm::InputTag>("src"))),
0085       deltaR_(iConfig.getParameter<double>("deltaR")),
0086       weightExpr_(iConfig.getParameter<std::string>("weight")),
0087       skipDefaultVeto_(iConfig.getParameter<bool>("skipDefaultVeto")),
0088       usePivotForBarrelEndcaps_(iConfig.getParameter<bool>("PivotCoordinatesForEBEE"))
0089 //,vetos_(new AbsVetos())
0090 {
0091   std::string mode = iConfig.getParameter<std::string>("mode");
0092   if (mode == "sum")
0093     mode_ = Sum;
0094   else if (mode == "sumRelative")
0095     mode_ = SumRelative;
0096   else if (mode == "sum2")
0097     mode_ = Sum2;
0098   else if (mode == "sum2Relative")
0099     mode_ = Sum2Relative;
0100   else if (mode == "max")
0101     mode_ = Max;
0102   else if (mode == "maxRelative")
0103     mode_ = MaxRelative;
0104   else if (mode == "nearestDR")
0105     mode_ = NearestDR;
0106   else if (mode == "count")
0107     mode_ = Count;
0108   else
0109     throw cms::Exception("Not Implemented") << "Mode '" << mode << "' not implemented. "
0110                                             << "Supported modes are 'sum', 'sumRelative', 'count'." <<
0111         //"Supported modes are 'sum', 'sumRelative', 'max', 'maxRelative', 'count'." << // TODO: on request only
0112         "New methods can be easily implemented if requested.";
0113   typedef std::vector<std::string> vstring;
0114   vstring vetos = iConfig.getParameter<vstring>("vetos");
0115   reco::isodeposit::EventDependentAbsVeto *evdep = nullptr;
0116   static const std::regex ecalSwitch("^Ecal(Barrel|Endcaps):(.*)");
0117
0118   for (vstring::const_iterator it = vetos.begin(), ed = vetos.end(); it != ed; ++it) {
0119     std::cmatch match;
0120     // in that case, make two series of vetoes
0121     if (usePivotForBarrelEndcaps_) {
0122       if (regex_match(it->c_str(), match, ecalSwitch)) {
0123         if (match[1] == "Barrel") {
0124           //        std::cout << " Adding Barrel veto " << std::string(match[2]) << std::endl;
0125           barrelVetos_.push_back(
0126               IsoDepositVetoFactory::make(std::string(match[2]).c_str(), evdep, iC));  // I don't know a better syntax
0127         }
0128         if (match[1] == "Endcaps") {
0129           //        std::cout << " Adding Endcap veto " << std::string(match[2]) << std::endl;
0130           endcapVetos_.push_back(IsoDepositVetoFactory::make(std::string(match[2]).c_str(), evdep, iC));
0131         }
0132       } else {
0133         barrelVetos_.push_back(IsoDepositVetoFactory::make(it->c_str(), evdep, iC));
0134         endcapVetos_.push_back(IsoDepositVetoFactory::make(it->c_str(), evdep, iC));
0135       }
0136     } else {
0137       //only one serie of vetoes, just barrel
0138       barrelVetos_.push_back(IsoDepositVetoFactory::make(it->c_str(), evdep, iC));
0139     }
0140     if (evdep)
0141       evdepVetos_.push_back(evdep);
0142   }
0143
0144   std::string weight = iConfig.getParameter<std::string>("weight");
0145   if (isNumber(weight)) {
0146     //std::cout << "Weight is a simple number, " << toNumber(weight) << std::endl;
0147     weight_ = toNumber(weight);
0148     usesFunction_ = false;
0149   } else {
0150     usesFunction_ = true;
0151     //std::cout << "Weight is a function, this might slow you down... " << std::endl;
0152   }
0153   //std::cout << "PFCandIsolatorFromDeposits::SingleDeposit::SingleDeposit: Total of " << vetos_.size() << " vetos" << std::endl;
0154 }
0155 void PFCandIsolatorFromDeposits::SingleDeposit::cleanup() {
0156   for (AbsVetos::iterator it = barrelVetos_.begin(), ed = barrelVetos_.end(); it != ed; ++it) {
0157     delete *it;
0158   }
0159   for (AbsVetos::iterator it = endcapVetos_.begin(), ed = endcapVetos_.end(); it != ed; ++it) {
0160     delete *it;
0161   }
0162   barrelVetos_.clear();
0163   endcapVetos_.clear();
0164   // NOTE: we DON'T have to delete the evdepVetos_, they have already been deleted above. We just clear the vectors
0165   evdepVetos_.clear();
0166 }
0167 void PFCandIsolatorFromDeposits::SingleDeposit::open(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0168   iEvent.getByToken(srcToken_, hDeps_);
0169   for (EventDependentAbsVetos::iterator it = evdepVetos_.begin(), ed = evdepVetos_.end(); it != ed; ++it) {
0170     (*it)->setEvent(iEvent, iSetup);
0171   }
0172 }
0173
0174 double PFCandIsolatorFromDeposits::SingleDeposit::compute(const reco::CandidateBaseRef &cand) {
0175   const IsoDeposit &dep = (*hDeps_)[cand];
0176   double eta = dep.eta(), phi = dep.phi();  // better to center on the deposit direction
0177                                             // that could be, e.g., the impact point at calo
0178   bool barrel = true;
0179   if (usePivotForBarrelEndcaps_) {
0180     const reco::PFCandidate *myPFCand = dynamic_cast<const reco::PFCandidate *>(&(*cand));
0181     if (myPFCand) {
0182       // exact barrel boundary
0183       barrel = fabs(myPFCand->positionAtECALEntrance().eta()) < 1.479;
0184     } else {
0185       const reco::RecoCandidate *myRecoCand = dynamic_cast<const reco::RecoCandidate *>(&(*cand));
0186       if (myRecoCand) {
0187         // not optimal. isEB should be used.
0188         barrel = (fabs(myRecoCand->superCluster()->eta()) < 1.479);
0189       }
0190     }
0191   }
0192   // if ! usePivotForBarrelEndcaps_ only the barrel series is used, which does not prevent the vetoes do be different in barrel & endcaps
0193   reco::isodeposit::AbsVetos *vetos = (barrel) ? &barrelVetos_ : &endcapVetos_;
0194
0195   for (AbsVetos::iterator it = vetos->begin(), ed = vetos->end(); it != ed; ++it) {
0196     (*it)->centerOn(eta, phi);
0197   }
0198   double weight = (usesFunction_ ? weightExpr_(*cand) : weight_);
0199   switch (mode_) {
0200     case Count:
0201       return weight * dep.countWithin(deltaR_, *vetos, skipDefaultVeto_);
0202     case Sum:
0203       return weight * dep.sumWithin(deltaR_, *vetos, skipDefaultVeto_);
0204     case SumRelative:
0205       return weight * dep.sumWithin(deltaR_, *vetos, skipDefaultVeto_) / dep.candEnergy();
0206     case Sum2:
0207       return weight * dep.sum2Within(deltaR_, *vetos, skipDefaultVeto_);
0208     case Sum2Relative:
0209       return weight * dep.sum2Within(deltaR_, *vetos, skipDefaultVeto_) / (dep.candEnergy() * dep.candEnergy());
0210     case Max:
0211       return weight * dep.maxWithin(deltaR_, *vetos, skipDefaultVeto_);
0212     case NearestDR:
0213       return weight * dep.nearestDR(deltaR_, *vetos, skipDefaultVeto_);
0214     case MaxRelative:
0215       return weight * dep.maxWithin(deltaR_, *vetos, skipDefaultVeto_) / dep.candEnergy();
0216   }
0217   throw cms::Exception("Logic error") << "Should not happen at " << __FILE__ << ", line "
0218                                       << __LINE__;  // avoid gcc warning
0219 }
0220
0221 /// constructor with config
0222 PFCandIsolatorFromDeposits::PFCandIsolatorFromDeposits(const ParameterSet &par) {
0223   typedef std::vector<edm::ParameterSet> VPSet;
0224   VPSet depPSets = par.getParameter<VPSet>("deposits");
0225   for (VPSet::const_iterator it = depPSets.begin(), ed = depPSets.end(); it != ed; ++it) {
0226     sources_.push_back(SingleDeposit(*it, consumesCollector()));
0227   }
0228   if (sources_.empty())
0229     throw cms::Exception("Configuration Error") << "Please specify at least one deposit!";
0230   produces<CandDoubleMap>();
0231 }
0232
0233 /// destructor
0234 PFCandIsolatorFromDeposits::~PFCandIsolatorFromDeposits() {
0235   std::vector<SingleDeposit>::iterator it, begin = sources_.begin(), end = sources_.end();
0236   for (it = begin; it != end; ++it)
0237     it->cleanup();
0238 }
0239
0240 /// build deposits
0241 void PFCandIsolatorFromDeposits::produce(Event &event, const EventSetup &eventSetup) {
0242   std::vector<SingleDeposit>::iterator it, begin = sources_.begin(), end = sources_.end();
0243   for (it = begin; it != end; ++it)
0244     it->open(event, eventSetup);
0245
0246   const IsoDepositMap &map = begin->map();
0247
0248   if (map.empty()) {  // !!???
0249     event.put(std::make_unique<CandDoubleMap>());
0250     return;
0251   }
0252   std::unique_ptr<CandDoubleMap> ret(new CandDoubleMap());
0253   CandDoubleMap::Filler filler(*ret);
0254
0255   typedef reco::IsoDepositMap::const_iterator iterator_i;
0256   typedef reco::IsoDepositMap::container::const_iterator iterator_ii;
0257   iterator_i depI = map.begin();
0258   iterator_i depIEnd = map.end();
0259   for (; depI != depIEnd; ++depI) {
0260     std::vector<double> retV(depI.size(), 0);
0261     edm::Handle<edm::View<reco::Candidate> > candH;
0262     event.get(depI.id(), candH);
0263     const edm::View<reco::Candidate> &candV = *candH;
0264
0265     iterator_ii depII = depI.begin();
0266     iterator_ii depIIEnd = depI.end();
0267     size_t iRet = 0;
0268     for (; depII != depIIEnd; ++depII, ++iRet) {
0269       double sum = 0;
0270       for (it = begin; it != end; ++it)
0271         sum += it->compute(candV.refAt(iRet));
0272       retV[iRet] = sum;
0273     }
0274     filler.insert(candH, retV.begin(), retV.end());
0275   }
0276   filler.fill();
0277   event.put(std::move(ret));
0278 }
0279
0280 DEFINE_FWK_MODULE(PFCandIsolatorFromDeposits);