File indexing completed on 2023-10-25 10:01:22
0001 #ifndef MuonIsolationProducers_MuIsolatorResultProducer_H
0002 #define MuonIsolationProducers_MuIsolatorResultProducer_H
0003
0004 #include "FWCore/Framework/interface/stream/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
0052 template <typename BT = reco::Candidate>
0053 class MuIsolatorResultProducer : public edm::stream::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(reco::TrackBase::Point const& theBeam,
0089 std::vector<reco::IsoDeposit::Vetos*>& vetos,
0090 CandMap& candMap) const;
0091
0092 template <typename RT>
0093 void writeOutImpl(edm::Event& event, const CandMap& candMapT, const Results& results) const;
0094
0095 void writeOut(edm::Event& event, const CandMap& candMap, const Results& results) const;
0096
0097 edm::ParameterSet theConfig;
0098 std::vector<DepositConf> theDepositConfs;
0099
0100
0101 bool theRemoveOtherVetos;
0102 VetoCuts theVetoCuts;
0103
0104
0105 Isolator* theIsolator;
0106 ResultType theResultType;
0107
0108
0109 std::string theBeamlineOption;
0110 edm::InputTag theBeamSpotLabel;
0111 };
0112
0113
0114 template <typename BT>
0115 template <typename RT>
0116 inline void MuIsolatorResultProducer<BT>::writeOutImpl(edm::Event& event,
0117 const CandMap& candMapT,
0118 const Results& results) const {
0119
0120 std::vector<RT> resV(results.size());
0121 for (unsigned int i = 0; i < resV.size(); ++i)
0122 resV[i] = results[i].val<RT>();
0123 auto outMap = std::make_unique<edm::ValueMap<RT>>();
0124 typename edm::ValueMap<RT>::Filler filler(*outMap);
0125
0126
0127 if (!candMapT.get().empty()) {
0128 filler.insert(candMapT.handle(), resV.begin(), resV.end());
0129 filler.fill();
0130 }
0131
0132 event.put(std::move(outMap));
0133 }
0134
0135
0136 template <typename BT>
0137 inline void MuIsolatorResultProducer<BT>::writeOut(edm::Event& event,
0138 const CandMap& candMapT,
0139 const Results& results) const {
0140 std::string metname = "RecoMuon|MuonIsolationProducers";
0141 LogDebug(metname) << "Before calling writeOutMap with result type " << theIsolator->resultType();
0142
0143 if (theResultType == Isolator::ISOL_INT_TYPE)
0144 writeOutImpl<int>(event, candMapT, results);
0145 if (theResultType == Isolator::ISOL_FLOAT_TYPE)
0146 writeOutImpl<float>(event, candMapT, results);
0147 if (theResultType == Isolator::ISOL_BOOL_TYPE)
0148 writeOutImpl<bool>(event, candMapT, results);
0149 }
0150
0151
0152 template <typename BT>
0153 inline void MuIsolatorResultProducer<BT>::callWhatProduces() {
0154 if (theResultType == Isolator::ISOL_FLOAT_TYPE)
0155 produces<edm::ValueMap<float>>();
0156 if (theResultType == Isolator::ISOL_INT_TYPE)
0157 produces<edm::ValueMap<int>>();
0158 if (theResultType == Isolator::ISOL_BOOL_TYPE)
0159 produces<edm::ValueMap<bool>>();
0160 }
0161
0162
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
0185 template <typename BT>
0186 MuIsolatorResultProducer<BT>::MuIsolatorResultProducer(const edm::ParameterSet& par)
0187 : theConfig(par), theRemoveOtherVetos(par.getParameter<bool>("RemoveOtherVetos")), theIsolator(nullptr) {
0188 LogDebug("RecoMuon|MuonIsolation") << " MuIsolatorResultProducer CTOR";
0189
0190
0191 std::vector<edm::ParameterSet> depositInputs = par.getParameter<std::vector<edm::ParameterSet>>("InputMuIsoDeposits");
0192
0193 std::vector<double> dWeights(depositInputs.size());
0194 std::vector<double> dThresholds(depositInputs.size());
0195
0196 for (unsigned int iDep = 0; iDep < depositInputs.size(); ++iDep) {
0197 DepositConf dConf;
0198 dConf.tag = depositInputs[iDep].getParameter<edm::InputTag>("DepositTag");
0199 dConf.weight = depositInputs[iDep].getParameter<double>("DepositWeight");
0200 dConf.threshold = depositInputs[iDep].getParameter<double>("DepositThreshold");
0201
0202 dWeights[iDep] = dConf.weight;
0203 dThresholds[iDep] = dConf.threshold;
0204
0205 theDepositConfs.push_back(dConf);
0206 }
0207
0208 edm::ParameterSet isoPset = par.getParameter<edm::ParameterSet>("IsolatorPSet");
0209
0210 std::string isolatorType = isoPset.getParameter<std::string>("ComponentName");
0211 if (isolatorType == "IsolatorByDeposit") {
0212 std::string coneSizeType = isoPset.getParameter<std::string>("ConeSizeType");
0213 if (coneSizeType == "FixedConeSize") {
0214 float coneSize(isoPset.getParameter<double>("coneSize"));
0215
0216 theIsolator = new muonisolation::IsolatorByDeposit(coneSize, dWeights, dThresholds);
0217
0218
0219 } else if (coneSizeType == "CutsConeSize") {
0220
0221
0222
0223
0224 }
0225 } else if (isolatorType == "IsolatorByNominalEfficiency") {
0226
0227 theIsolator =
0228 new muonisolation::IsolatorByNominalEfficiency("noname", std::vector<std::string>(1, "8:0.97"), dWeights);
0229 } else if (isolatorType == "IsolatorByDepositCount") {
0230 std::string coneSizeType = isoPset.getParameter<std::string>("ConeSizeType");
0231 if (coneSizeType == "FixedConeSize") {
0232 float coneSize(isoPset.getParameter<double>("coneSize"));
0233
0234 theIsolator = new muonisolation::IsolatorByDepositCount(coneSize, dThresholds);
0235
0236
0237 } else if (coneSizeType == "CutsConeSize") {
0238
0239
0240
0241 }
0242 }
0243
0244 if (theIsolator == nullptr) {
0245 edm::LogError("MuonIsolationProducers") << "Failed to initialize an isolator";
0246 }
0247 theResultType = theIsolator->resultType();
0248
0249 callWhatProduces();
0250
0251 if (theRemoveOtherVetos) {
0252 edm::ParameterSet vetoPSet = par.getParameter<edm::ParameterSet>("VetoPSet");
0253 theVetoCuts.selectAll = vetoPSet.getParameter<bool>("SelectAll");
0254
0255
0256
0257
0258 if (!theVetoCuts.selectAll) {
0259 theVetoCuts.muAbsEtaMax = vetoPSet.getParameter<double>("MuAbsEtaMax");
0260 theVetoCuts.muPtMin = vetoPSet.getParameter<double>("MuPtMin");
0261 theVetoCuts.muAbsZMax = vetoPSet.getParameter<double>("MuAbsZMax");
0262 theVetoCuts.muD0Max = vetoPSet.getParameter<double>("MuD0Max");
0263 theBeamlineOption = par.getParameter<std::string>("BeamlineOption");
0264 theBeamSpotLabel = par.getParameter<edm::InputTag>("BeamSpotLabel");
0265 }
0266 }
0267 }
0268
0269
0270 template <typename BT>
0271 MuIsolatorResultProducer<BT>::~MuIsolatorResultProducer() {
0272 if (theIsolator)
0273 delete theIsolator;
0274 LogDebug("RecoMuon|MuIsolatorResultProducer") << " MuIsolatorResultProducer DTOR";
0275 }
0276
0277 template <typename BT>
0278 void MuIsolatorResultProducer<BT>::produce(edm::Event& event, const edm::EventSetup& eventSetup) {
0279 std::string metname = "RecoMuon|MuonIsolationProducers";
0280 LogDebug(metname) << " Muon Deposit producing..."
0281 << " BEGINING OF EVENT "
0282 << "================================";
0283
0284 reco::TrackBase::Point theBeam = reco::TrackBase::Point(0, 0, 0);
0285
0286
0287 if (theRemoveOtherVetos && !theVetoCuts.selectAll) {
0288 if (theBeamlineOption == "BeamSpotFromEvent") {
0289
0290 reco::BeamSpot beamSpot;
0291 edm::Handle<reco::BeamSpot> beamSpotH;
0292
0293 event.getByLabel(theBeamSpotLabel, beamSpotH);
0294
0295 if (beamSpotH.isValid()) {
0296 theBeam = beamSpotH->position();
0297 LogTrace(metname) << "Extracted beam point at " << theBeam << std::endl;
0298 }
0299 }
0300 }
0301
0302
0303
0304
0305
0306
0307 CandMap candMapT;
0308
0309 unsigned int colSize = initAssociation(event, candMapT);
0310
0311
0312 Results results(colSize);
0313
0314
0315 std::vector<reco::IsoDeposit::Vetos*> vetoDeps(theDepositConfs.size(), nullptr);
0316
0317 if (colSize != 0) {
0318 if (theRemoveOtherVetos) {
0319 initVetos(theBeam, vetoDeps, candMapT);
0320 }
0321
0322
0323 for (unsigned int muI = 0; muI < colSize; ++muI) {
0324 results[muI] = theIsolator->result(candMapT.get()[muI].second, *(candMapT.get()[muI].first));
0325
0326 if (results[muI].typeF() != theIsolator->resultType()) {
0327 edm::LogError("MuonIsolationProducers") << "Failed to get result from the isolator";
0328 }
0329 }
0330 }
0331
0332 LogDebug(metname) << "Ready to write out results of size " << results.size();
0333 writeOut(event, candMapT, results);
0334
0335 for (unsigned int iDep = 0; iDep < vetoDeps.size(); ++iDep) {
0336
0337 if (vetoDeps[iDep]) {
0338 delete vetoDeps[iDep];
0339 vetoDeps[iDep] = nullptr;
0340 }
0341 }
0342 }
0343
0344 template <typename BT>
0345 unsigned int MuIsolatorResultProducer<BT>::initAssociation(edm::Event& event, CandMap& candMapT) const {
0346 std::string metname = "RecoMuon|MuonIsolationProducers";
0347
0348 typedef reco::IsoDepositMap::container CT;
0349
0350 for (unsigned int iMap = 0; iMap < theDepositConfs.size(); ++iMap) {
0351 edm::Handle<reco::IsoDepositMap> depH;
0352 event.getByLabel(theDepositConfs[iMap].tag, depH);
0353 LogDebug(metname) << "Got Deposits of size " << depH->size();
0354 if (depH->empty())
0355 continue;
0356
0357
0358
0359 typename CandMap::handle_type keyH;
0360 event.get(depH->begin().id(), keyH);
0361 candMapT.setHandle(keyH);
0362 typename CT::const_iterator depHCI = depH->begin().begin();
0363 typename CT::const_iterator depEnd = depH->begin().end();
0364 unsigned int keyI = 0;
0365 for (; depHCI != depEnd; ++depHCI, ++keyI) {
0366 typename CandMap::key_type muPtr(keyH->refAt(keyI));
0367
0368 if (iMap == 0)
0369 candMapT.get().push_back(typename CandMap::pair_type(muPtr, DepositContainer(theDepositConfs.size())));
0370 typename CandMap::iterator muI = candMapT.get().begin();
0371 for (; muI != candMapT.get().end(); ++muI) {
0372 if (muI->first == muPtr)
0373 break;
0374 }
0375 if (muI->first != muPtr) {
0376 edm::LogError("MuonIsolationProducers") << "Failed to align muon map";
0377 }
0378 muI->second[iMap].dep = &*depHCI;
0379 }
0380 }
0381
0382 LogDebug(metname) << "Picked and aligned nDeps = " << candMapT.get().size();
0383 return candMapT.get().size();
0384 }
0385
0386 template <typename BT>
0387 void MuIsolatorResultProducer<BT>::initVetos(reco::TrackBase::Point const& theBeam,
0388 std::vector<reco::IsoDeposit::Vetos*>& vetos,
0389 CandMap& candMapT) const {
0390 std::string metname = "RecoMuon|MuonIsolationProducers";
0391
0392 if (theRemoveOtherVetos) {
0393 LogDebug(metname) << "Start checking for vetos based on input/expected vetos.size of " << vetos.size()
0394 << " passed at " << &vetos << " and an input map.size of " << candMapT.get().size();
0395
0396 unsigned int muI = 0;
0397 for (; muI < candMapT.get().size(); ++muI) {
0398 typename CandMap::key_type mu = candMapT.get()[muI].first;
0399 double d0 = ((mu->vx() - theBeam.x()) * mu->py() - (mu->vy() - theBeam.y()) * mu->px()) / mu->pt();
0400 LogDebug(metname) << "Muon at index " << muI;
0401 if (theVetoCuts.selectAll || (fabs(mu->eta()) < theVetoCuts.muAbsEtaMax && mu->pt() > theVetoCuts.muPtMin &&
0402 fabs(mu->vz()) < theVetoCuts.muAbsZMax && fabs(d0) < theVetoCuts.muD0Max)) {
0403 LogDebug(metname) << "muon passes the cuts";
0404 for (unsigned int iDep = 0; iDep < candMapT.get()[muI].second.size(); ++iDep) {
0405 if (vetos[iDep] == nullptr)
0406 vetos[iDep] = new reco::IsoDeposit::Vetos();
0407
0408 vetos[iDep]->push_back(candMapT.get()[muI].second[iDep].dep->veto());
0409 }
0410 }
0411 }
0412
0413 LogDebug(metname) << "Assigning vetos";
0414 muI = 0;
0415 for (; muI < candMapT.get().size(); ++muI) {
0416 for (unsigned int iDep = 0; iDep < candMapT.get()[muI].second.size(); ++iDep) {
0417 candMapT.get()[muI].second[iDep].vetos = vetos[iDep];
0418 }
0419 }
0420 LogDebug(metname) << "Done with vetos";
0421 }
0422 }
0423
0424 #endif