Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:26:01

0001 #ifndef MuonIsolationProducers_MuIsolatorResultProducer_H
0002 #define MuonIsolationProducers_MuIsolatorResultProducer_H
0003 
0004 #include "FWCore/Framework/interface/EDProducer.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 
0007 #include "RecoMuon/MuonIsolation/interface/MuIsoBaseIsolator.h"
0008 
0009 #include "FWCore/Framework/interface/Event.h"
0010 
0011 #include "DataFormats/Common/interface/RefToBase.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 
0014 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0015 #include "DataFormats/Common/interface/ValueMap.h"
0016 
0017 #include <string>
0018 
0019 namespace edm {
0020   class EventSetup;
0021 }
0022 
0023 struct muisorhelper {
0024   typedef muonisolation::MuIsoBaseIsolator Isolator;
0025   typedef Isolator::Result Result;
0026   typedef Isolator::ResultType ResultType;
0027   typedef std::vector<Result> Results;
0028   typedef Isolator::DepositContainer DepositContainer;
0029 
0030   template <typename BT>
0031   class CandMap {
0032   public:
0033     typedef typename edm::RefToBase<BT> key_type;
0034     typedef typename edm::Handle<edm::View<BT>> handle_type;
0035     typedef DepositContainer value_type;
0036     typedef std::pair<key_type, value_type> pair_type;
0037     typedef typename std::vector<pair_type> map_type;
0038     typedef typename map_type::iterator iterator;
0039 
0040     map_type& get() { return cMap_; }
0041     const map_type& get() const { return cMap_; }
0042     const handle_type& handle() const { return handle_; }
0043     void setHandle(const handle_type& rhs) { handle_ = rhs; }
0044 
0045   private:
0046     map_type cMap_;
0047     handle_type handle_;
0048   };
0049 };
0050 
0051 //! BT == base type
0052 template <typename BT = reco::Candidate>
0053 class MuIsolatorResultProducer : public edm::EDProducer {
0054 public:
0055   MuIsolatorResultProducer(const edm::ParameterSet&);
0056 
0057   ~MuIsolatorResultProducer() override;
0058 
0059   void produce(edm::Event&, const edm::EventSetup&) override;
0060 
0061 private:
0062   typedef muisorhelper::Isolator Isolator;
0063   typedef muisorhelper::Result Result;
0064   typedef muisorhelper::ResultType ResultType;
0065   typedef muisorhelper::Results Results;
0066   typedef muisorhelper::DepositContainer DepositContainer;
0067 
0068   typedef muisorhelper::CandMap<BT> CandMap;
0069 
0070   struct DepositConf {
0071     edm::InputTag tag;
0072     double weight;
0073     double threshold;
0074   };
0075 
0076   struct VetoCuts {
0077     bool selectAll;
0078     double muAbsEtaMax;
0079     double muPtMin;
0080     double muAbsZMax;
0081     double muD0Max;
0082   };
0083 
0084   void callWhatProduces();
0085 
0086   unsigned int initAssociation(edm::Event& event, CandMap& candMapT) const;
0087 
0088   void initVetos(std::vector<reco::IsoDeposit::Vetos*>& vetos, CandMap& candMap) const;
0089 
0090   template <typename RT>
0091   void writeOutImpl(edm::Event& event, const CandMap& candMapT, const Results& results) const;
0092 
0093   void writeOut(edm::Event& event, const CandMap& candMap, const Results& results) const;
0094 
0095   edm::ParameterSet theConfig;
0096   std::vector<DepositConf> theDepositConfs;
0097 
0098   //!choose which muon vetos should be removed from all deposits
0099   bool theRemoveOtherVetos;
0100   VetoCuts theVetoCuts;
0101 
0102   //!the isolator
0103   Isolator* theIsolator;
0104   ResultType theResultType;
0105 
0106   //! beam spot
0107   std::string theBeamlineOption;
0108   edm::InputTag theBeamSpotLabel;
0109   reco::TrackBase::Point theBeam;
0110 };
0111 
0112 //! actually do the writing here
0113 template <typename BT>
0114 template <typename RT>
0115 inline void MuIsolatorResultProducer<BT>::writeOutImpl(edm::Event& event,
0116                                                        const CandMap& candMapT,
0117                                                        const Results& results) const {
0118   //! make an output vec of what's to be written with a concrete type
0119   std::vector<RT> resV(results.size());
0120   for (unsigned int i = 0; i < resV.size(); ++i)
0121     resV[i] = results[i].val<RT>();
0122   auto outMap = std::make_unique<edm::ValueMap<RT>>();
0123   typename edm::ValueMap<RT>::Filler filler(*outMap);
0124 
0125   //! fill/insert of non-empty values only
0126   if (!candMapT.get().empty()) {
0127     filler.insert(candMapT.handle(), resV.begin(), resV.end());
0128     filler.fill();
0129   }
0130 
0131   event.put(std::move(outMap));
0132 }
0133 
0134 //! choose which result type to write here
0135 template <typename BT>
0136 inline void MuIsolatorResultProducer<BT>::writeOut(edm::Event& event,
0137                                                    const CandMap& candMapT,
0138                                                    const Results& results) const {
0139   std::string metname = "RecoMuon|MuonIsolationProducers";
0140   LogDebug(metname) << "Before calling writeOutMap  with result type " << theIsolator->resultType();
0141 
0142   if (theResultType == Isolator::ISOL_INT_TYPE)
0143     writeOutImpl<int>(event, candMapT, results);
0144   if (theResultType == Isolator::ISOL_FLOAT_TYPE)
0145     writeOutImpl<float>(event, candMapT, results);
0146   if (theResultType == Isolator::ISOL_BOOL_TYPE)
0147     writeOutImpl<bool>(event, candMapT, results);
0148 }
0149 
0150 //! declare what's going to be produced
0151 template <typename BT>
0152 inline void MuIsolatorResultProducer<BT>::callWhatProduces() {
0153   if (theResultType == Isolator::ISOL_FLOAT_TYPE)
0154     produces<edm::ValueMap<float>>();
0155   if (theResultType == Isolator::ISOL_INT_TYPE)
0156     produces<edm::ValueMap<int>>();
0157   if (theResultType == Isolator::ISOL_BOOL_TYPE)
0158     produces<edm::ValueMap<bool>>();
0159 }
0160 
0161 // Framework
0162 #include "FWCore/Framework/interface/EDProducer.h"
0163 #include "FWCore/Framework/interface/EventSetup.h"
0164 #include "DataFormats/Common/interface/Handle.h"
0165 
0166 #include "FWCore/Framework/interface/ESHandle.h"
0167 
0168 #include "DataFormats/MuonReco/interface/Muon.h"
0169 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
0170 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
0171 #include "DataFormats/Common/interface/RefToBase.h"
0172 #include "DataFormats/TrackReco/interface/Track.h"
0173 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0174 
0175 #include "RecoMuon/MuonIsolation/interface/Range.h"
0176 #include "DataFormats/RecoCandidate/interface/IsoDepositDirection.h"
0177 
0178 #include "RecoMuon/MuonIsolation/interface/IsolatorByDeposit.h"
0179 #include "RecoMuon/MuonIsolation/interface/IsolatorByDepositCount.h"
0180 #include "RecoMuon/MuonIsolation/interface/IsolatorByNominalEfficiency.h"
0181 
0182 #include <string>
0183 
0184 //! constructor with config
0185 template <typename BT>
0186 MuIsolatorResultProducer<BT>::MuIsolatorResultProducer(const edm::ParameterSet& par)
0187     : theConfig(par),
0188       theRemoveOtherVetos(par.getParameter<bool>("RemoveOtherVetos")),
0189       theIsolator(nullptr),
0190       theBeam(0, 0, 0) {
0191   LogDebug("RecoMuon|MuonIsolation") << " MuIsolatorResultProducer CTOR";
0192 
0193   //! read input config for deposit types and weights and thresholds to apply to them
0194   std::vector<edm::ParameterSet> depositInputs = par.getParameter<std::vector<edm::ParameterSet>>("InputMuIsoDeposits");
0195 
0196   std::vector<double> dWeights(depositInputs.size());
0197   std::vector<double> dThresholds(depositInputs.size());
0198 
0199   for (unsigned int iDep = 0; iDep < depositInputs.size(); ++iDep) {
0200     DepositConf dConf;
0201     dConf.tag = depositInputs[iDep].getParameter<edm::InputTag>("DepositTag");
0202     dConf.weight = depositInputs[iDep].getParameter<double>("DepositWeight");
0203     dConf.threshold = depositInputs[iDep].getParameter<double>("DepositThreshold");
0204 
0205     dWeights[iDep] = dConf.weight;
0206     dThresholds[iDep] = dConf.threshold;
0207 
0208     theDepositConfs.push_back(dConf);
0209   }
0210 
0211   edm::ParameterSet isoPset = par.getParameter<edm::ParameterSet>("IsolatorPSet");
0212   //! will switch to a factory at some point
0213   std::string isolatorType = isoPset.getParameter<std::string>("ComponentName");
0214   if (isolatorType == "IsolatorByDeposit") {
0215     std::string coneSizeType = isoPset.getParameter<std::string>("ConeSizeType");
0216     if (coneSizeType == "FixedConeSize") {
0217       float coneSize(isoPset.getParameter<double>("coneSize"));
0218 
0219       theIsolator = new muonisolation::IsolatorByDeposit(coneSize, dWeights, dThresholds);
0220 
0221       //      theIsolator = new IsolatorByDeposit(isoPset);
0222     } else if (coneSizeType == "CutsConeSize") {
0223       //! FIXME
0224       //       Cuts cuts(isoPset.getParameter<edm::ParameterSet>("CutsPSet"));
0225 
0226       //       theIsolator = new IsolatorByDeposit(coneSize, dWeights, dThresholds);
0227     }
0228   } else if (isolatorType == "IsolatorByNominalEfficiency") {
0229     //! FIXME: need to get the file name here
0230     theIsolator =
0231         new muonisolation::IsolatorByNominalEfficiency("noname", std::vector<std::string>(1, "8:0.97"), dWeights);
0232   } else if (isolatorType == "IsolatorByDepositCount") {
0233     std::string coneSizeType = isoPset.getParameter<std::string>("ConeSizeType");
0234     if (coneSizeType == "FixedConeSize") {
0235       float coneSize(isoPset.getParameter<double>("coneSize"));
0236 
0237       theIsolator = new muonisolation::IsolatorByDepositCount(coneSize, dThresholds);
0238 
0239       //      theIsolator = new IsolatorByDeposit(isoPset);
0240     } else if (coneSizeType == "CutsConeSize") {
0241       //       Cuts cuts(isoPset.getParameter<edm::ParameterSet>("CutsPSet"));
0242 
0243       //       theIsolator = new IsolatorByDeposit(coneSize, dWeights, dThresholds);
0244     }
0245   }
0246 
0247   if (theIsolator == nullptr) {
0248     edm::LogError("MuonIsolationProducers") << "Failed to initialize an isolator";
0249   }
0250   theResultType = theIsolator->resultType();
0251 
0252   callWhatProduces();
0253 
0254   if (theRemoveOtherVetos) {
0255     edm::ParameterSet vetoPSet = par.getParameter<edm::ParameterSet>("VetoPSet");
0256     theVetoCuts.selectAll = vetoPSet.getParameter<bool>("SelectAll");
0257 
0258     //! "other vetoes" is limited to the same collection now
0259     //! for non-trivial choice an external map with pre-made selection flags
0260     //! can be a better choice
0261     if (!theVetoCuts.selectAll) {
0262       theVetoCuts.muAbsEtaMax = vetoPSet.getParameter<double>("MuAbsEtaMax");
0263       theVetoCuts.muPtMin = vetoPSet.getParameter<double>("MuPtMin");
0264       theVetoCuts.muAbsZMax = vetoPSet.getParameter<double>("MuAbsZMax");
0265       theVetoCuts.muD0Max = vetoPSet.getParameter<double>("MuD0Max");
0266       theBeamlineOption = par.getParameter<std::string>("BeamlineOption");
0267       theBeamSpotLabel = par.getParameter<edm::InputTag>("BeamSpotLabel");
0268     }
0269   }
0270 }
0271 
0272 //! destructor
0273 template <typename BT>
0274 MuIsolatorResultProducer<BT>::~MuIsolatorResultProducer() {
0275   if (theIsolator)
0276     delete theIsolator;
0277   LogDebug("RecoMuon|MuIsolatorResultProducer") << " MuIsolatorResultProducer DTOR";
0278 }
0279 
0280 template <typename BT>
0281 void MuIsolatorResultProducer<BT>::produce(edm::Event& event, const edm::EventSetup& eventSetup) {
0282   std::string metname = "RecoMuon|MuonIsolationProducers";
0283   LogDebug(metname) << " Muon Deposit producing..."
0284                     << " BEGINING OF EVENT "
0285                     << "================================";
0286 
0287   theBeam = reco::TrackBase::Point(0, 0, 0);
0288 
0289   //! do it only if needed
0290   if (theRemoveOtherVetos && !theVetoCuts.selectAll) {
0291     if (theBeamlineOption == "BeamSpotFromEvent") {
0292       //pick beamSpot
0293       reco::BeamSpot beamSpot;
0294       edm::Handle<reco::BeamSpot> beamSpotH;
0295 
0296       event.getByLabel(theBeamSpotLabel, beamSpotH);
0297 
0298       if (beamSpotH.isValid()) {
0299         theBeam = beamSpotH->position();
0300         LogTrace(metname) << "Extracted beam point at " << theBeam << std::endl;
0301       }
0302     }
0303   }
0304 
0305   //! "smart" container
0306   //! used to repackage deposits_type_candIndex into deposits_candIndex_type
0307   //! have to have it for veto removal (could do away without it otherwise)
0308   //! IMPORTANT: ALL THE REFERENCING BUSINESS IS DONE THROUGH POINTERS
0309   //! Access to the mapped values as reference type HAS TO BE AVAILABLE
0310   CandMap candMapT;
0311 
0312   unsigned int colSize = initAssociation(event, candMapT);
0313 
0314   //! isolator results will be here
0315   Results results(colSize);
0316 
0317   //! extra vetos will be filled here
0318   std::vector<reco::IsoDeposit::Vetos*> vetoDeps(theDepositConfs.size(), nullptr);
0319 
0320   if (colSize != 0) {
0321     if (theRemoveOtherVetos) {
0322       initVetos(vetoDeps, candMapT);
0323     }
0324 
0325     //! call the isolator result, passing {[deposit,vetos]_type} set and the candidate
0326     for (unsigned int muI = 0; muI < colSize; ++muI) {
0327       results[muI] = theIsolator->result(candMapT.get()[muI].second, *(candMapT.get()[muI].first));
0328 
0329       if (results[muI].typeF() != theIsolator->resultType()) {
0330         edm::LogError("MuonIsolationProducers") << "Failed to get result from the isolator";
0331       }
0332     }
0333   }
0334 
0335   LogDebug(metname) << "Ready to write out results of size " << results.size();
0336   writeOut(event, candMapT, results);
0337 
0338   for (unsigned int iDep = 0; iDep < vetoDeps.size(); ++iDep) {
0339     //! do cleanup
0340     if (vetoDeps[iDep]) {
0341       delete vetoDeps[iDep];
0342       vetoDeps[iDep] = nullptr;
0343     }
0344   }
0345 }
0346 
0347 template <typename BT>
0348 unsigned int MuIsolatorResultProducer<BT>::initAssociation(edm::Event& event, CandMap& candMapT) const {
0349   std::string metname = "RecoMuon|MuonIsolationProducers";
0350 
0351   typedef reco::IsoDepositMap::container CT;
0352 
0353   for (unsigned int iMap = 0; iMap < theDepositConfs.size(); ++iMap) {
0354     edm::Handle<reco::IsoDepositMap> depH;
0355     event.getByLabel(theDepositConfs[iMap].tag, depH);
0356     LogDebug(metname) << "Got Deposits of size " << depH->size();
0357     if (depH->empty())
0358       continue;
0359 
0360     //! WARNING: the input ValueMaps are better be for a single key product ID
0361     //! no effort is done (FIXME) for more complex cases
0362     typename CandMap::handle_type keyH;
0363     event.get(depH->begin().id(), keyH);
0364     candMapT.setHandle(keyH);
0365     typename CT::const_iterator depHCI = depH->begin().begin();
0366     typename CT::const_iterator depEnd = depH->begin().end();
0367     unsigned int keyI = 0;
0368     for (; depHCI != depEnd; ++depHCI, ++keyI) {
0369       typename CandMap::key_type muPtr(keyH->refAt(keyI));
0370       //! init {muon, {[deposit,veto]_type}} container
0371       if (iMap == 0)
0372         candMapT.get().push_back(typename CandMap::pair_type(muPtr, DepositContainer(theDepositConfs.size())));
0373       typename CandMap::iterator muI = candMapT.get().begin();
0374       for (; muI != candMapT.get().end(); ++muI) {
0375         if (muI->first == muPtr)
0376           break;
0377       }
0378       if (muI->first != muPtr) {
0379         edm::LogError("MuonIsolationProducers") << "Failed to align muon map";
0380       }
0381       muI->second[iMap].dep = &*depHCI;
0382     }
0383   }
0384 
0385   LogDebug(metname) << "Picked and aligned nDeps = " << candMapT.get().size();
0386   return candMapT.get().size();
0387 }
0388 
0389 template <typename BT>
0390 void MuIsolatorResultProducer<BT>::initVetos(std::vector<reco::IsoDeposit::Vetos*>& vetos, CandMap& candMapT) const {
0391   std::string metname = "RecoMuon|MuonIsolationProducers";
0392 
0393   if (theRemoveOtherVetos) {
0394     LogDebug(metname) << "Start checking for vetos based on input/expected vetos.size of " << vetos.size()
0395                       << " passed at " << &vetos << " and an input map.size of " << candMapT.get().size();
0396 
0397     unsigned int muI = 0;
0398     for (; muI < candMapT.get().size(); ++muI) {
0399       typename CandMap::key_type mu = candMapT.get()[muI].first;
0400       double d0 = ((mu->vx() - theBeam.x()) * mu->py() - (mu->vy() - theBeam.y()) * mu->px()) / mu->pt();
0401       LogDebug(metname) << "Muon at index " << muI;
0402       if (theVetoCuts.selectAll || (fabs(mu->eta()) < theVetoCuts.muAbsEtaMax && mu->pt() > theVetoCuts.muPtMin &&
0403                                     fabs(mu->vz()) < theVetoCuts.muAbsZMax && fabs(d0) < theVetoCuts.muD0Max)) {
0404         LogDebug(metname) << "muon passes the cuts";
0405         for (unsigned int iDep = 0; iDep < candMapT.get()[muI].second.size(); ++iDep) {
0406           if (vetos[iDep] == nullptr)
0407             vetos[iDep] = new reco::IsoDeposit::Vetos();
0408 
0409           vetos[iDep]->push_back(candMapT.get()[muI].second[iDep].dep->veto());
0410         }
0411       }
0412     }
0413 
0414     LogDebug(metname) << "Assigning vetos";
0415     muI = 0;
0416     for (; muI < candMapT.get().size(); ++muI) {
0417       for (unsigned int iDep = 0; iDep < candMapT.get()[muI].second.size(); ++iDep) {
0418         candMapT.get()[muI].second[iDep].vetos = vetos[iDep];
0419       }
0420     }
0421     LogDebug(metname) << "Done with vetos";
0422   }
0423 }
0424 
0425 #endif