File indexing completed on 2024-06-06 04:27:07
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0013 #include "DataFormats/Common/interface/Association.h"
0014 #include "DataFormats/Common/interface/TriggerResults.h"
0015 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0016 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0017 #include "DataFormats/Math/interface/deltaPhi.h"
0018 #include "DataFormats/Math/interface/deltaR.h"
0019 #include "DataFormats/MuonReco/interface/Muon.h"
0020 #include "DataFormats/MuonReco/interface/MuonSimInfo.h"
0021 #include "DataFormats/MuonReco/interface/MuonTimeExtra.h"
0022 #include "DataFormats/ParticleFlowCandidate/interface/IsolatedPFCandidate.h"
0023 #include "DataFormats/PatCandidates/interface/Muon.h"
0024 #include "DataFormats/PatCandidates/interface/PFIsolation.h"
0025 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0026 #include "DataFormats/PatCandidates/interface/TriggerObjectStandAlone.h"
0027 #include "DataFormats/PatCandidates/interface/UserData.h"
0028 #include "DataFormats/VertexReco/interface/Vertex.h"
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/stream/EDProducer.h"
0031 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0032 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0033 #include "FWCore/ParameterSet/interface/EmptyGroupDescription.h"
0034 #include "FWCore/ParameterSet/interface/FileInPath.h"
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0037 #include "FWCore/Utilities/interface/Exception.h"
0038 #include "FWCore/Utilities/interface/transform.h"
0039 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0040 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0041 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0042 #include "PhysicsTools/PatAlgos/interface/EfficiencyLoader.h"
0043 #include "PhysicsTools/PatAlgos/interface/KinResolutionsLoader.h"
0044 #include "PhysicsTools/PatAlgos/interface/MultiIsolator.h"
0045 #include "PhysicsTools/PatAlgos/interface/CalculatePtRatioRel.h"
0046 #include "PhysicsTools/PatAlgos/interface/MuonMvaIDEstimator.h"
0047 #include "PhysicsTools/PatAlgos/interface/PATUserDataHelper.h"
0048 #include "PhysicsTools/PatAlgos/interface/SoftMuonMvaEstimator.h"
0049 #include "PhysicsTools/PatAlgos/interface/SoftMuonMvaRun3Estimator.h"
0050 #include "PhysicsTools/PatUtils/interface/MiniIsolation.h"
0051 #include "PhysicsTools/XGBoost/interface/XGBooster.h"
0052 #include "TrackingTools/IPTools/interface/IPTools.h"
0053 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0054 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0055 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0056
0057 namespace pat {
0058
0059 class PATMuonHeavyObjectCache {
0060 public:
0061 PATMuonHeavyObjectCache(const edm::ParameterSet&);
0062
0063 pat::CalculatePtRatioRel const& calculatePtRatioRel() const { return *calculatePtRatioRel_; }
0064 pat::MuonMvaIDEstimator const& muonMvaIDEstimator() const { return *muonMvaIDEstimator_; }
0065 pat::SoftMuonMvaEstimator const& softMuonMvaEstimator() const { return *softMuonMvaEstimator_; }
0066
0067 private:
0068 std::unique_ptr<const pat::CalculatePtRatioRel> calculatePtRatioRel_;
0069 std::unique_ptr<const pat::MuonMvaIDEstimator> muonMvaIDEstimator_;
0070 std::unique_ptr<const pat::SoftMuonMvaEstimator> softMuonMvaEstimator_;
0071 };
0072
0073
0074 class TrackerIsolationPt;
0075 class CaloIsolationEnergy;
0076
0077
0078 class PATMuonProducer : public edm::stream::EDProducer<edm::GlobalCache<PATMuonHeavyObjectCache>> {
0079 public:
0080
0081 explicit PATMuonProducer(const edm::ParameterSet& iConfig, PATMuonHeavyObjectCache const*);
0082
0083 ~PATMuonProducer() override;
0084
0085 static std::unique_ptr<PATMuonHeavyObjectCache> initializeGlobalCache(const edm::ParameterSet& iConfig) {
0086 return std::make_unique<PATMuonHeavyObjectCache>(iConfig);
0087 }
0088
0089 static void globalEndJob(PATMuonHeavyObjectCache*) {}
0090
0091
0092 void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0093
0094 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0095
0096 private:
0097
0098 typedef edm::RefToBase<reco::Muon> MuonBaseRef;
0099 typedef std::vector<edm::Handle<edm::Association<reco::GenParticleCollection>>> GenAssociations;
0100 typedef std::vector<edm::Handle<edm::ValueMap<IsoDeposit>>> IsoDepositMaps;
0101 typedef std::vector<edm::Handle<edm::ValueMap<double>>> IsolationValueMaps;
0102 typedef std::pair<pat::IsolationKeys, edm::InputTag> IsolationLabel;
0103 typedef std::vector<IsolationLabel> IsolationLabels;
0104
0105
0106 void fillMuon(Muon& aMuon,
0107 const MuonBaseRef& muonRef,
0108 const reco::CandidateBaseRef& baseRef,
0109 const GenAssociations& genMatches,
0110 const IsoDepositMaps& deposits,
0111 const IsolationValueMaps& isolationValues) const;
0112
0113
0114 template <typename T>
0115 void readIsolationLabels(const edm::ParameterSet& iConfig,
0116 const char* psetName,
0117 IsolationLabels& labels,
0118 std::vector<edm::EDGetTokenT<edm::ValueMap<T>>>& tokens);
0119
0120 void setMuonMiniIso(pat::Muon& aMuon, const pat::PackedCandidateCollection* pc);
0121 double getRelMiniIsoPUCorrected(const pat::Muon& muon, double rho, const std::vector<double>& area);
0122
0123 double puppiCombinedIsolation(const pat::Muon& muon, const pat::PackedCandidateCollection* pc);
0124 bool isNeutralHadron(long pdgid);
0125 bool isChargedHadron(long pdgid);
0126 bool isPhoton(long pdgid);
0127
0128
0129
0130 void embedHighLevel(pat::Muon& aMuon,
0131 reco::TrackRef track,
0132 reco::TransientTrack& tt,
0133 reco::Vertex& primaryVertex,
0134 bool primaryVertexIsValid,
0135 reco::BeamSpot& beamspot,
0136 bool beamspotIsValid);
0137 double relMiniIsoPUCorrected(const pat::Muon& aMuon, double rho);
0138 std::optional<GlobalPoint> getMuonDirection(const reco::MuonChamberMatch& chamberMatch,
0139 const edm::ESHandle<GlobalTrackingGeometry>& geometry,
0140 const DetId& chamberId);
0141 void fillL1TriggerInfo(pat::Muon& muon,
0142 edm::Handle<std::vector<pat::TriggerObjectStandAlone>>& triggerObjects,
0143 const edm::TriggerNames& names,
0144 const edm::ESHandle<GlobalTrackingGeometry>& geometry);
0145 void fillHltTriggerInfo(pat::Muon& muon,
0146 edm::Handle<std::vector<pat::TriggerObjectStandAlone>>& triggerObjects,
0147 const edm::TriggerNames& names,
0148 const std::vector<std::string>& collection_names);
0149
0150 private:
0151
0152 edm::EDGetTokenT<edm::View<reco::Muon>> muonToken_;
0153
0154
0155 edm::EDGetTokenT<pat::PackedCandidateCollection> pcToken_;
0156 bool computeMiniIso_;
0157 bool computePuppiCombinedIso_;
0158 std::vector<double> effectiveAreaVec_;
0159 std::vector<double> miniIsoParams_;
0160 double relMiniIsoPUCorrected_;
0161
0162 double mvaIDtightCut_;
0163 double mvaIDmediumCut_;
0164
0165 bool embedBestTrack_;
0166
0167 bool embedTunePBestTrack_;
0168
0169 bool forceEmbedBestTrack_;
0170
0171 bool embedTrack_;
0172
0173 bool embedStandAloneMuon_;
0174
0175 bool embedCombinedMuon_;
0176
0177 bool embedCaloMETMuonCorrs_;
0178
0179 edm::EDGetTokenT<edm::ValueMap<reco::MuonMETCorrectionData>> caloMETMuonCorrsToken_;
0180
0181 bool embedTcMETMuonCorrs_;
0182
0183 edm::EDGetTokenT<edm::ValueMap<reco::MuonMETCorrectionData>> tcMETMuonCorrsToken_;
0184
0185 bool embedPickyMuon_;
0186
0187 bool embedTpfmsMuon_;
0188
0189 bool embedDytMuon_;
0190
0191 bool addInverseBeta_;
0192
0193 edm::EDGetTokenT<edm::ValueMap<reco::MuonTimeExtra>> muonTimeExtraToken_;
0194
0195 bool addGenMatch_;
0196
0197 std::vector<edm::EDGetTokenT<edm::Association<reco::GenParticleCollection>>> genMatchTokens_;
0198
0199 bool embedGenMatch_;
0200
0201 bool addResolutions_;
0202
0203 pat::helper::KinResolutionsLoader resolutionLoader_;
0204
0205 bool useParticleFlow_;
0206
0207 edm::EDGetTokenT<reco::PFCandidateCollection> pfMuonToken_;
0208
0209 bool embedPFCandidate_;
0210
0211 bool embedHighLevelSelection_;
0212
0213 edm::EDGetTokenT<reco::BeamSpot> beamLineToken_;
0214
0215 edm::EDGetTokenT<std::vector<reco::Vertex>> pvToken_;
0216
0217 IsolationLabels isoDepositLabels_;
0218 std::vector<edm::EDGetTokenT<edm::ValueMap<IsoDeposit>>> isoDepositTokens_;
0219
0220 IsolationLabels isolationValueLabels_;
0221 std::vector<edm::EDGetTokenT<edm::ValueMap<double>>> isolationValueTokens_;
0222
0223 bool addEfficiencies_;
0224
0225 bool useUserData_;
0226
0227 bool embedPfEcalEnergy_;
0228
0229 bool addPuppiIsolation_;
0230
0231 edm::EDGetTokenT<edm::ValueMap<float>> PUPPIIsolation_charged_hadrons_;
0232 edm::EDGetTokenT<edm::ValueMap<float>> PUPPIIsolation_neutral_hadrons_;
0233 edm::EDGetTokenT<edm::ValueMap<float>> PUPPIIsolation_photons_;
0234
0235 edm::EDGetTokenT<edm::ValueMap<float>> PUPPINoLeptonsIsolation_charged_hadrons_;
0236 edm::EDGetTokenT<edm::ValueMap<float>> PUPPINoLeptonsIsolation_neutral_hadrons_;
0237 edm::EDGetTokenT<edm::ValueMap<float>> PUPPINoLeptonsIsolation_photons_;
0238
0239 bool computeMuonIDMVA_;
0240 bool computeSoftMuonMVA_;
0241 std::unique_ptr<pat::XGBooster> softMuonMvaRun3Booster_;
0242 bool recomputeBasicSelectors_;
0243 bool useJec_;
0244 edm::EDGetTokenT<reco::JetTagCollection> mvaBTagCollectionTag_;
0245 edm::EDGetTokenT<reco::JetCorrector> mvaL1Corrector_;
0246 edm::EDGetTokenT<reco::JetCorrector> mvaL1L2L3ResCorrector_;
0247 edm::EDGetTokenT<double> rho_;
0248
0249
0250
0251 pat::helper::MultiIsolator isolator_;
0252
0253 pat::helper::MultiIsolator::IsolationValuePairs isolatorTmpStorage_;
0254
0255 pat::helper::EfficiencyLoader efficiencyLoader_;
0256
0257 pat::PATUserDataHelper<pat::Muon> userDataHelper_;
0258
0259
0260 edm::EDGetTokenT<edm::ValueMap<reco::MuonSimInfo>> simInfo_;
0261
0262
0263 bool addTriggerMatching_;
0264 edm::EDGetTokenT<std::vector<pat::TriggerObjectStandAlone>> triggerObjects_;
0265 edm::EDGetTokenT<edm::TriggerResults> triggerResults_;
0266 std::vector<std::string> hltCollectionFilters_;
0267
0268 const edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> geometryToken_;
0269 const edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transientTrackBuilderToken_;
0270
0271 const edm::EDPutTokenT<std::vector<Muon>> patMuonPutToken_;
0272 };
0273
0274 }
0275
0276 template <typename T>
0277 void pat::PATMuonProducer::readIsolationLabels(const edm::ParameterSet& iConfig,
0278 const char* psetName,
0279 pat::PATMuonProducer::IsolationLabels& labels,
0280 std::vector<edm::EDGetTokenT<edm::ValueMap<T>>>& tokens) {
0281 labels.clear();
0282
0283 if (iConfig.exists(psetName)) {
0284 edm::ParameterSet depconf = iConfig.getParameter<edm::ParameterSet>(psetName);
0285
0286 if (depconf.exists("tracker"))
0287 labels.emplace_back(pat::TrackIso, depconf.getParameter<edm::InputTag>("tracker"));
0288 if (depconf.exists("ecal"))
0289 labels.emplace_back(pat::EcalIso, depconf.getParameter<edm::InputTag>("ecal"));
0290 if (depconf.exists("hcal"))
0291 labels.emplace_back(pat::HcalIso, depconf.getParameter<edm::InputTag>("hcal"));
0292 if (depconf.exists("pfAllParticles")) {
0293 labels.emplace_back(pat::PfAllParticleIso, depconf.getParameter<edm::InputTag>("pfAllParticles"));
0294 }
0295 if (depconf.exists("pfChargedHadrons")) {
0296 labels.emplace_back(pat::PfChargedHadronIso, depconf.getParameter<edm::InputTag>("pfChargedHadrons"));
0297 }
0298 if (depconf.exists("pfChargedAll")) {
0299 labels.emplace_back(pat::PfChargedAllIso, depconf.getParameter<edm::InputTag>("pfChargedAll"));
0300 }
0301 if (depconf.exists("pfPUChargedHadrons")) {
0302 labels.emplace_back(pat::PfPUChargedHadronIso, depconf.getParameter<edm::InputTag>("pfPUChargedHadrons"));
0303 }
0304 if (depconf.exists("pfNeutralHadrons")) {
0305 labels.emplace_back(pat::PfNeutralHadronIso, depconf.getParameter<edm::InputTag>("pfNeutralHadrons"));
0306 }
0307 if (depconf.exists("pfPhotons")) {
0308 labels.emplace_back(pat::PfGammaIso, depconf.getParameter<edm::InputTag>("pfPhotons"));
0309 }
0310 if (depconf.exists("user")) {
0311 std::vector<edm::InputTag> userdeps = depconf.getParameter<std::vector<edm::InputTag>>("user");
0312 std::vector<edm::InputTag>::const_iterator it = userdeps.begin(), ed = userdeps.end();
0313 int key = pat::IsolationKeys::UserBaseIso;
0314 for (; it != ed; ++it, ++key) {
0315 labels.push_back(std::make_pair(pat::IsolationKeys(key), *it));
0316 }
0317 }
0318 }
0319 tokens = edm::vector_transform(labels, [this](pat::PATMuonProducer::IsolationLabel const& label) {
0320 return consumes<edm::ValueMap<T>>(label.second);
0321 });
0322 }
0323
0324 using namespace pat;
0325 using namespace std;
0326
0327 PATMuonHeavyObjectCache::PATMuonHeavyObjectCache(const edm::ParameterSet& iConfig) {
0328 if (iConfig.getParameter<bool>("computeMiniIso")) {
0329 float mvaDrMax = iConfig.getParameter<double>("mvaDrMax");
0330 calculatePtRatioRel_ = std::make_unique<CalculatePtRatioRel>(mvaDrMax * mvaDrMax);
0331 }
0332
0333 if (iConfig.getParameter<bool>("computeMuonIDMVA")) {
0334 edm::FileInPath mvaIDTrainingFile = iConfig.getParameter<edm::FileInPath>("mvaIDTrainingFile");
0335 muonMvaIDEstimator_ = std::make_unique<MuonMvaIDEstimator>(mvaIDTrainingFile);
0336 }
0337
0338 if (iConfig.getParameter<bool>("computeSoftMuonMVA")) {
0339 edm::FileInPath softMvaTrainingFile = iConfig.getParameter<edm::FileInPath>("softMvaTrainingFile");
0340 softMuonMvaEstimator_ = std::make_unique<SoftMuonMvaEstimator>(softMvaTrainingFile);
0341 }
0342 }
0343
0344 PATMuonProducer::PATMuonProducer(const edm::ParameterSet& iConfig, PATMuonHeavyObjectCache const*)
0345 : relMiniIsoPUCorrected_(0),
0346 useUserData_(iConfig.exists("userData")),
0347 computeMuonIDMVA_(false),
0348 computeSoftMuonMVA_(false),
0349 recomputeBasicSelectors_(false),
0350 useJec_(false),
0351 isolator_(iConfig.getParameter<edm::ParameterSet>("userIsolation"), consumesCollector(), false),
0352 geometryToken_{esConsumes()},
0353 transientTrackBuilderToken_{esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))},
0354 patMuonPutToken_{produces<std::vector<Muon>>()} {
0355
0356 mvaIDmediumCut_ = iConfig.getParameter<double>("mvaIDwpMedium");
0357 mvaIDtightCut_ = iConfig.getParameter<double>("mvaIDwpTight");
0358
0359
0360 muonToken_ = consumes<edm::View<reco::Muon>>(iConfig.getParameter<edm::InputTag>("muonSource"));
0361
0362 embedBestTrack_ = iConfig.getParameter<bool>("embedMuonBestTrack");
0363 embedTunePBestTrack_ = iConfig.getParameter<bool>("embedTunePMuonBestTrack");
0364 forceEmbedBestTrack_ = iConfig.getParameter<bool>("forceBestTrackEmbedding");
0365 embedTrack_ = iConfig.getParameter<bool>("embedTrack");
0366 embedCombinedMuon_ = iConfig.getParameter<bool>("embedCombinedMuon");
0367 embedStandAloneMuon_ = iConfig.getParameter<bool>("embedStandAloneMuon");
0368
0369 embedCaloMETMuonCorrs_ = iConfig.getParameter<bool>("embedCaloMETMuonCorrs");
0370 embedTcMETMuonCorrs_ = iConfig.getParameter<bool>("embedTcMETMuonCorrs");
0371 caloMETMuonCorrsToken_ =
0372 mayConsume<edm::ValueMap<reco::MuonMETCorrectionData>>(iConfig.getParameter<edm::InputTag>("caloMETMuonCorrs"));
0373 tcMETMuonCorrsToken_ =
0374 mayConsume<edm::ValueMap<reco::MuonMETCorrectionData>>(iConfig.getParameter<edm::InputTag>("tcMETMuonCorrs"));
0375
0376 useParticleFlow_ = iConfig.getParameter<bool>("useParticleFlow");
0377 embedPFCandidate_ = iConfig.getParameter<bool>("embedPFCandidate");
0378 pfMuonToken_ = mayConsume<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("pfMuonSource"));
0379 embedPfEcalEnergy_ = iConfig.getParameter<bool>("embedPfEcalEnergy");
0380
0381 embedPickyMuon_ = iConfig.getParameter<bool>("embedPickyMuon");
0382 embedTpfmsMuon_ = iConfig.getParameter<bool>("embedTpfmsMuon");
0383 embedDytMuon_ = iConfig.getParameter<bool>("embedDytMuon");
0384
0385 addInverseBeta_ = iConfig.getParameter<bool>("addInverseBeta");
0386 if (addInverseBeta_) {
0387 muonTimeExtraToken_ =
0388 consumes<edm::ValueMap<reco::MuonTimeExtra>>(iConfig.getParameter<edm::InputTag>("sourceMuonTimeExtra"));
0389 }
0390
0391 addGenMatch_ = iConfig.getParameter<bool>("addGenMatch");
0392 if (addGenMatch_) {
0393 embedGenMatch_ = iConfig.getParameter<bool>("embedGenMatch");
0394 if (iConfig.existsAs<edm::InputTag>("genParticleMatch")) {
0395 genMatchTokens_.push_back(consumes<edm::Association<reco::GenParticleCollection>>(
0396 iConfig.getParameter<edm::InputTag>("genParticleMatch")));
0397 } else {
0398 genMatchTokens_ = edm::vector_transform(
0399 iConfig.getParameter<std::vector<edm::InputTag>>("genParticleMatch"),
0400 [this](edm::InputTag const& tag) { return consumes<edm::Association<reco::GenParticleCollection>>(tag); });
0401 }
0402 }
0403
0404 addEfficiencies_ = iConfig.getParameter<bool>("addEfficiencies");
0405 if (addEfficiencies_) {
0406 efficiencyLoader_ =
0407 pat::helper::EfficiencyLoader(iConfig.getParameter<edm::ParameterSet>("efficiencies"), consumesCollector());
0408 }
0409
0410 addResolutions_ = iConfig.getParameter<bool>("addResolutions");
0411 if (addResolutions_) {
0412 resolutionLoader_ =
0413 pat::helper::KinResolutionsLoader(iConfig.getParameter<edm::ParameterSet>("resolutions"), consumesCollector());
0414 }
0415
0416 addPuppiIsolation_ = iConfig.getParameter<bool>("addPuppiIsolation");
0417 if (addPuppiIsolation_) {
0418 PUPPIIsolation_charged_hadrons_ =
0419 consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("puppiIsolationChargedHadrons"));
0420 PUPPIIsolation_neutral_hadrons_ =
0421 consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("puppiIsolationNeutralHadrons"));
0422 PUPPIIsolation_photons_ =
0423 consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("puppiIsolationPhotons"));
0424
0425 PUPPINoLeptonsIsolation_charged_hadrons_ =
0426 consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("puppiNoLeptonsIsolationChargedHadrons"));
0427 PUPPINoLeptonsIsolation_neutral_hadrons_ =
0428 consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("puppiNoLeptonsIsolationNeutralHadrons"));
0429 PUPPINoLeptonsIsolation_photons_ =
0430 consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("puppiNoLeptonsIsolationPhotons"));
0431 }
0432
0433 readIsolationLabels(iConfig, "isoDeposits", isoDepositLabels_, isoDepositTokens_);
0434
0435 readIsolationLabels(iConfig, "isolationValues", isolationValueLabels_, isolationValueTokens_);
0436
0437 if (useUserData_) {
0438 userDataHelper_ = PATUserDataHelper<Muon>(iConfig.getParameter<edm::ParameterSet>("userData"), consumesCollector());
0439 }
0440
0441 embedHighLevelSelection_ = iConfig.getParameter<bool>("embedHighLevelSelection");
0442 if (embedHighLevelSelection_) {
0443 beamLineToken_ = consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamLineSrc"));
0444 pvToken_ = consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("pvSrc"));
0445 }
0446
0447
0448 computeMiniIso_ = iConfig.getParameter<bool>("computeMiniIso");
0449
0450 computePuppiCombinedIso_ = iConfig.getParameter<bool>("computePuppiCombinedIso");
0451
0452 effectiveAreaVec_ = iConfig.getParameter<std::vector<double>>("effectiveAreaVec");
0453
0454 miniIsoParams_ = iConfig.getParameter<std::vector<double>>("miniIsoParams");
0455 if (computeMiniIso_ && miniIsoParams_.size() != 9) {
0456 throw cms::Exception("ParameterError") << "miniIsoParams must have exactly 9 elements.\n";
0457 }
0458 if (computeMiniIso_ || computePuppiCombinedIso_)
0459 pcToken_ = consumes<pat::PackedCandidateCollection>(iConfig.getParameter<edm::InputTag>("pfCandsForMiniIso"));
0460
0461
0462 recomputeBasicSelectors_ = iConfig.getParameter<bool>("recomputeBasicSelectors");
0463 computeMuonIDMVA_ = iConfig.getParameter<bool>("computeMuonIDMVA");
0464
0465 if (computeMiniIso_) {
0466
0467 mvaBTagCollectionTag_ = consumes<reco::JetTagCollection>(iConfig.getParameter<edm::InputTag>("mvaJetTag"));
0468 mvaL1Corrector_ = consumes<reco::JetCorrector>(iConfig.getParameter<edm::InputTag>("mvaL1Corrector"));
0469 mvaL1L2L3ResCorrector_ = consumes<reco::JetCorrector>(iConfig.getParameter<edm::InputTag>("mvaL1L2L3ResCorrector"));
0470 rho_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rho"));
0471 useJec_ = iConfig.getParameter<bool>("useJec");
0472 }
0473
0474 computeSoftMuonMVA_ = iConfig.getParameter<bool>("computeSoftMuonMVA");
0475
0476
0477 simInfo_ = consumes<edm::ValueMap<reco::MuonSimInfo>>(iConfig.getParameter<edm::InputTag>("muonSimInfo"));
0478
0479 if (computeSoftMuonMVA_) {
0480 std::string softMvaRun3Model = iConfig.getParameter<string>("softMvaRun3Model");
0481 softMuonMvaRun3Booster_ =
0482 std::make_unique<pat::XGBooster>(edm::FileInPath(softMvaRun3Model + ".model").fullPath(),
0483 edm::FileInPath(softMvaRun3Model + ".features").fullPath());
0484 }
0485
0486 addTriggerMatching_ = iConfig.getParameter<bool>("addTriggerMatching");
0487 if (addTriggerMatching_) {
0488 triggerObjects_ =
0489 consumes<std::vector<pat::TriggerObjectStandAlone>>(iConfig.getParameter<edm::InputTag>("triggerObjects"));
0490 triggerResults_ = consumes<edm::TriggerResults>(iConfig.getParameter<edm::InputTag>("triggerResults"));
0491 }
0492 hltCollectionFilters_ = iConfig.getParameter<std::vector<std::string>>("hltCollectionFilters");
0493 }
0494
0495 PATMuonProducer::~PATMuonProducer() {}
0496
0497 std::optional<GlobalPoint> PATMuonProducer::getMuonDirection(const reco::MuonChamberMatch& chamberMatch,
0498 const edm::ESHandle<GlobalTrackingGeometry>& geometry,
0499 const DetId& chamberId) {
0500 const GeomDet* chamberGeometry = geometry->idToDet(chamberId);
0501 if (chamberGeometry) {
0502 LocalPoint localPosition(chamberMatch.x, chamberMatch.y, 0);
0503 return std::optional<GlobalPoint>(std::in_place, chamberGeometry->toGlobal(localPosition));
0504 }
0505 return std::optional<GlobalPoint>();
0506 }
0507
0508 void PATMuonProducer::fillL1TriggerInfo(pat::Muon& aMuon,
0509 edm::Handle<std::vector<pat::TriggerObjectStandAlone>>& triggerObjects,
0510 const edm::TriggerNames& names,
0511 const edm::ESHandle<GlobalTrackingGeometry>& geometry) {
0512
0513
0514
0515
0516
0517 std::optional<GlobalPoint> muonPosition;
0518
0519
0520
0521
0522 for (const auto& chamberMatch : aMuon.matches()) {
0523 if (chamberMatch.id.subdetId() == MuonSubdetId::DT) {
0524 DTChamberId detId(chamberMatch.id.rawId());
0525 if (abs(detId.station()) > 3)
0526 continue;
0527 muonPosition = getMuonDirection(chamberMatch, geometry, detId);
0528 if (abs(detId.station()) == 2)
0529 break;
0530 }
0531 if (chamberMatch.id.subdetId() == MuonSubdetId::CSC) {
0532 CSCDetId detId(chamberMatch.id.rawId());
0533 if (abs(detId.station()) > 3)
0534 continue;
0535 muonPosition = getMuonDirection(chamberMatch, geometry, detId);
0536 if (abs(detId.station()) == 2)
0537 break;
0538 }
0539 }
0540 if (not muonPosition)
0541 return;
0542 for (const auto& triggerObject : *triggerObjects) {
0543 if (triggerObject.hasTriggerObjectType(trigger::TriggerL1Mu)) {
0544 if (std::abs(triggerObject.eta()) < 0.001) {
0545
0546 if (deltaPhi(triggerObject.phi(), muonPosition->phi()) > 0.1)
0547 continue;
0548 } else {
0549
0550 if (deltaR(triggerObject.p4(), *muonPosition) > 0.15)
0551 continue;
0552 }
0553 pat::TriggerObjectStandAlone obj(triggerObject);
0554 obj.unpackPathNames(names);
0555 aMuon.addTriggerObjectMatch(obj);
0556 }
0557 }
0558 }
0559
0560 void PATMuonProducer::fillHltTriggerInfo(pat::Muon& muon,
0561 edm::Handle<std::vector<pat::TriggerObjectStandAlone>>& triggerObjects,
0562 const edm::TriggerNames& names,
0563 const std::vector<std::string>& collection_filter_names) {
0564
0565
0566 for (const auto& triggerObject : *triggerObjects) {
0567 if (triggerObject.hasTriggerObjectType(trigger::TriggerMuon)) {
0568 bool keepIt = false;
0569 for (const auto& name : collection_filter_names) {
0570 if (triggerObject.hasCollection(name)) {
0571 keepIt = true;
0572 break;
0573 }
0574 }
0575 if (not keepIt)
0576 continue;
0577 if (deltaR(triggerObject.p4(), muon) > 0.1)
0578 continue;
0579 pat::TriggerObjectStandAlone obj(triggerObject);
0580 obj.unpackPathNames(names);
0581 muon.addTriggerObjectMatch(obj);
0582 }
0583 }
0584 }
0585
0586 void PATMuonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0587
0588 auto geometry = iSetup.getHandle(geometryToken_);
0589 if (!geometry.isValid())
0590 throw cms::Exception("FatalError") << "Unable to find GlobalTrackingGeometryRecord in event!\n";
0591
0592
0593 if (iEvent.isRealData()) {
0594 addGenMatch_ = false;
0595 embedGenMatch_ = false;
0596 }
0597
0598 edm::Handle<edm::View<reco::Muon>> muons;
0599 iEvent.getByToken(muonToken_, muons);
0600
0601 edm::Handle<pat::PackedCandidateCollection> pc;
0602 if (computeMiniIso_ || computePuppiCombinedIso_)
0603 iEvent.getByToken(pcToken_, pc);
0604
0605
0606
0607 TransientTrackBuilder const* trackBuilder = nullptr;
0608
0609 if (isolator_.enabled())
0610 isolator_.beginEvent(iEvent, iSetup);
0611 if (efficiencyLoader_.enabled())
0612 efficiencyLoader_.newEvent(iEvent);
0613 if (resolutionLoader_.enabled())
0614 resolutionLoader_.newEvent(iEvent, iSetup);
0615
0616 IsoDepositMaps deposits(isoDepositTokens_.size());
0617 for (size_t j = 0; j < isoDepositTokens_.size(); ++j) {
0618 iEvent.getByToken(isoDepositTokens_[j], deposits[j]);
0619 }
0620
0621 IsolationValueMaps isolationValues(isolationValueTokens_.size());
0622 for (size_t j = 0; j < isolationValueTokens_.size(); ++j) {
0623 iEvent.getByToken(isolationValueTokens_[j], isolationValues[j]);
0624 }
0625
0626
0627 edm::Handle<edm::ValueMap<float>> PUPPIIsolation_charged_hadrons;
0628 edm::Handle<edm::ValueMap<float>> PUPPIIsolation_neutral_hadrons;
0629 edm::Handle<edm::ValueMap<float>> PUPPIIsolation_photons;
0630
0631 edm::Handle<edm::ValueMap<float>> PUPPINoLeptonsIsolation_charged_hadrons;
0632 edm::Handle<edm::ValueMap<float>> PUPPINoLeptonsIsolation_neutral_hadrons;
0633 edm::Handle<edm::ValueMap<float>> PUPPINoLeptonsIsolation_photons;
0634 if (addPuppiIsolation_) {
0635
0636 iEvent.getByToken(PUPPIIsolation_charged_hadrons_, PUPPIIsolation_charged_hadrons);
0637 iEvent.getByToken(PUPPIIsolation_neutral_hadrons_, PUPPIIsolation_neutral_hadrons);
0638 iEvent.getByToken(PUPPIIsolation_photons_, PUPPIIsolation_photons);
0639
0640 iEvent.getByToken(PUPPINoLeptonsIsolation_charged_hadrons_, PUPPINoLeptonsIsolation_charged_hadrons);
0641 iEvent.getByToken(PUPPINoLeptonsIsolation_neutral_hadrons_, PUPPINoLeptonsIsolation_neutral_hadrons);
0642 iEvent.getByToken(PUPPINoLeptonsIsolation_photons_, PUPPINoLeptonsIsolation_photons);
0643 }
0644
0645
0646 edm::Handle<reco::JetTagCollection> mvaBTagCollectionTag;
0647 edm::Handle<reco::JetCorrector> mvaL1Corrector;
0648 edm::Handle<reco::JetCorrector> mvaL1L2L3ResCorrector;
0649 if (computeMiniIso_) {
0650 iEvent.getByToken(mvaBTagCollectionTag_, mvaBTagCollectionTag);
0651 iEvent.getByToken(mvaL1Corrector_, mvaL1Corrector);
0652 iEvent.getByToken(mvaL1L2L3ResCorrector_, mvaL1L2L3ResCorrector);
0653 }
0654
0655
0656 GenAssociations genMatches(genMatchTokens_.size());
0657 if (addGenMatch_) {
0658 for (size_t j = 0, nd = genMatchTokens_.size(); j < nd; ++j) {
0659 iEvent.getByToken(genMatchTokens_[j], genMatches[j]);
0660 }
0661 }
0662
0663
0664
0665 reco::Vertex primaryVertex;
0666 reco::BeamSpot beamSpot;
0667 bool beamSpotIsValid = false;
0668 bool primaryVertexIsValid = false;
0669 if (embedHighLevelSelection_) {
0670
0671 edm::Handle<reco::BeamSpot> beamSpotHandle;
0672 iEvent.getByToken(beamLineToken_, beamSpotHandle);
0673
0674
0675 edm::Handle<std::vector<reco::Vertex>> pvHandle;
0676 iEvent.getByToken(pvToken_, pvHandle);
0677
0678 if (beamSpotHandle.isValid()) {
0679 beamSpot = *beamSpotHandle;
0680 beamSpotIsValid = true;
0681 } else {
0682 edm::LogError("DataNotAvailable") << "No beam spot available from EventSetup, not adding high level selection \n";
0683 }
0684 if (pvHandle.isValid() && !pvHandle->empty()) {
0685 primaryVertex = pvHandle->at(0);
0686 primaryVertexIsValid = true;
0687 } else {
0688 edm::LogError("DataNotAvailable")
0689 << "No primary vertex available from EventSetup, not adding high level selection \n";
0690 }
0691
0692 trackBuilder = &iSetup.getData(transientTrackBuilderToken_);
0693 }
0694
0695
0696 edm::Handle<edm::ValueMap<reco::MuonSimInfo>> simInfo;
0697 bool simInfoIsAvailalbe = iEvent.getByToken(simInfo_, simInfo);
0698
0699
0700 std::vector<Muon> patMuons;
0701
0702 edm::Handle<reco::PFCandidateCollection> pfMuons;
0703 if (useParticleFlow_) {
0704
0705 iEvent.getByToken(pfMuonToken_, pfMuons);
0706
0707 unsigned index = 0;
0708 for (reco::PFCandidateConstIterator i = pfMuons->begin(); i != pfMuons->end(); ++i, ++index) {
0709 const reco::PFCandidate& pfmu = *i;
0710
0711 const reco::MuonRef& muonRef = pfmu.muonRef();
0712 assert(muonRef.isNonnull());
0713
0714 MuonBaseRef muonBaseRef(muonRef);
0715 Muon aMuon(muonBaseRef);
0716
0717 if (useUserData_) {
0718 userDataHelper_.add(aMuon, iEvent, iSetup);
0719 }
0720
0721
0722 if (embedHighLevelSelection_) {
0723
0724 reco::TrackRef innerTrack = muonBaseRef->innerTrack();
0725 reco::TrackRef globalTrack = muonBaseRef->globalTrack();
0726 reco::TrackRef bestTrack = muonBaseRef->muonBestTrack();
0727 reco::TrackRef chosenTrack = innerTrack;
0728
0729 if (bestTrack.isNonnull() && bestTrack.isAvailable())
0730 chosenTrack = bestTrack;
0731
0732 if (chosenTrack.isNonnull() && chosenTrack.isAvailable()) {
0733 unsigned int nhits = chosenTrack->numberOfValidHits();
0734 aMuon.setNumberOfValidHits(nhits);
0735
0736 reco::TransientTrack tt = trackBuilder->build(chosenTrack);
0737 embedHighLevel(aMuon, chosenTrack, tt, primaryVertex, primaryVertexIsValid, beamSpot, beamSpotIsValid);
0738 }
0739
0740 if (globalTrack.isNonnull() && globalTrack.isAvailable() && !embedCombinedMuon_) {
0741 double norm_chi2 = globalTrack->chi2() / globalTrack->ndof();
0742 aMuon.setNormChi2(norm_chi2);
0743 }
0744 }
0745 reco::PFCandidateRef pfRef(pfMuons, index);
0746
0747 reco::CandidateBaseRef pfBaseRef(pfRef);
0748
0749 aMuon.setPFCandidateRef(pfRef);
0750 if (embedPFCandidate_)
0751 aMuon.embedPFCandidate();
0752 fillMuon(aMuon, muonBaseRef, pfBaseRef, genMatches, deposits, isolationValues);
0753
0754 if (computeMiniIso_)
0755 setMuonMiniIso(aMuon, pc.product());
0756
0757 if (addPuppiIsolation_) {
0758 aMuon.setIsolationPUPPI((*PUPPIIsolation_charged_hadrons)[muonBaseRef],
0759 (*PUPPIIsolation_neutral_hadrons)[muonBaseRef],
0760 (*PUPPIIsolation_photons)[muonBaseRef]);
0761
0762 aMuon.setIsolationPUPPINoLeptons((*PUPPINoLeptonsIsolation_charged_hadrons)[muonBaseRef],
0763 (*PUPPINoLeptonsIsolation_neutral_hadrons)[muonBaseRef],
0764 (*PUPPINoLeptonsIsolation_photons)[muonBaseRef]);
0765 } else {
0766 aMuon.setIsolationPUPPI(-999., -999., -999.);
0767 aMuon.setIsolationPUPPINoLeptons(-999., -999., -999.);
0768 }
0769
0770 if (embedPfEcalEnergy_) {
0771 aMuon.setPfEcalEnergy(pfmu.ecalEnergy());
0772 }
0773
0774 patMuons.push_back(aMuon);
0775 }
0776 } else {
0777 edm::Handle<edm::View<reco::Muon>> muons;
0778 iEvent.getByToken(muonToken_, muons);
0779
0780
0781 edm::Handle<edm::ValueMap<reco::MuonMETCorrectionData>> caloMETMuonCorrs;
0782
0783 if (embedCaloMETMuonCorrs_) {
0784 iEvent.getByToken(caloMETMuonCorrsToken_, caloMETMuonCorrs);
0785
0786 }
0787 edm::Handle<edm::ValueMap<reco::MuonMETCorrectionData>> tcMETMuonCorrs;
0788
0789 if (embedTcMETMuonCorrs_) {
0790 iEvent.getByToken(tcMETMuonCorrsToken_, tcMETMuonCorrs);
0791
0792 }
0793
0794 if (embedPfEcalEnergy_ || embedPFCandidate_) {
0795
0796 iEvent.getByToken(pfMuonToken_, pfMuons);
0797 }
0798
0799 edm::Handle<edm::ValueMap<reco::MuonTimeExtra>> muonsTimeExtra;
0800 if (addInverseBeta_) {
0801
0802 iEvent.getByToken(muonTimeExtraToken_, muonsTimeExtra);
0803 }
0804
0805 for (edm::View<reco::Muon>::const_iterator itMuon = muons->begin(); itMuon != muons->end(); ++itMuon) {
0806
0807 unsigned int idx = itMuon - muons->begin();
0808 MuonBaseRef muonRef = muons->refAt(idx);
0809 reco::CandidateBaseRef muonBaseRef(muonRef);
0810
0811 Muon aMuon(muonRef);
0812 fillMuon(aMuon, muonRef, muonBaseRef, genMatches, deposits, isolationValues);
0813 if (computeMiniIso_)
0814 setMuonMiniIso(aMuon, pc.product());
0815 if (addPuppiIsolation_) {
0816 aMuon.setIsolationPUPPI((*PUPPIIsolation_charged_hadrons)[muonRef],
0817 (*PUPPIIsolation_neutral_hadrons)[muonRef],
0818 (*PUPPIIsolation_photons)[muonRef]);
0819 aMuon.setIsolationPUPPINoLeptons((*PUPPINoLeptonsIsolation_charged_hadrons)[muonRef],
0820 (*PUPPINoLeptonsIsolation_neutral_hadrons)[muonRef],
0821 (*PUPPINoLeptonsIsolation_photons)[muonRef]);
0822 } else {
0823 aMuon.setIsolationPUPPI(-999., -999., -999.);
0824 aMuon.setIsolationPUPPINoLeptons(-999., -999., -999.);
0825 }
0826
0827
0828 if (isolator_.enabled()) {
0829
0830 isolator_.fill(*muons, idx, isolatorTmpStorage_);
0831 typedef pat::helper::MultiIsolator::IsolationValuePairs IsolationValuePairs;
0832
0833 for (IsolationValuePairs::const_reverse_iterator it = isolatorTmpStorage_.rbegin(),
0834 ed = isolatorTmpStorage_.rend();
0835 it != ed;
0836 ++it) {
0837 aMuon.setIsolation(it->first, it->second);
0838 }
0839 }
0840
0841
0842
0843
0844
0845
0846
0847 edm::Ptr<reco::Muon> muonsPtr = muons->ptrAt(idx);
0848 if (useUserData_) {
0849 userDataHelper_.add(aMuon, iEvent, iSetup);
0850 }
0851
0852
0853 if (embedHighLevelSelection_) {
0854
0855 reco::TrackRef innerTrack = itMuon->innerTrack();
0856 reco::TrackRef globalTrack = itMuon->globalTrack();
0857 reco::TrackRef bestTrack = itMuon->muonBestTrack();
0858 reco::TrackRef chosenTrack = innerTrack;
0859
0860 if (bestTrack.isNonnull() && bestTrack.isAvailable())
0861 chosenTrack = bestTrack;
0862 if (chosenTrack.isNonnull() && chosenTrack.isAvailable()) {
0863 unsigned int nhits = chosenTrack->numberOfValidHits();
0864 aMuon.setNumberOfValidHits(nhits);
0865
0866 reco::TransientTrack tt = trackBuilder->build(chosenTrack);
0867 embedHighLevel(aMuon, chosenTrack, tt, primaryVertex, primaryVertexIsValid, beamSpot, beamSpotIsValid);
0868 }
0869
0870 if (globalTrack.isNonnull() && globalTrack.isAvailable()) {
0871 double norm_chi2 = globalTrack->chi2() / globalTrack->ndof();
0872 aMuon.setNormChi2(norm_chi2);
0873 }
0874 }
0875
0876
0877 if (embedCaloMETMuonCorrs_)
0878 aMuon.embedCaloMETMuonCorrs((*caloMETMuonCorrs)[muonRef]);
0879 if (embedTcMETMuonCorrs_)
0880 aMuon.embedTcMETMuonCorrs((*tcMETMuonCorrs)[muonRef]);
0881
0882 if (embedPfEcalEnergy_ || embedPFCandidate_) {
0883 if (embedPfEcalEnergy_)
0884 aMuon.setPfEcalEnergy(-99.0);
0885 unsigned index = 0;
0886 for (const reco::PFCandidate& pfmu : *pfMuons) {
0887 if (pfmu.muonRef().isNonnull()) {
0888 if (pfmu.muonRef().id() != muonRef.id())
0889 throw cms::Exception("Configuration")
0890 << "Muon reference within PF candidates does not point to the muon collection." << std::endl;
0891 if (pfmu.muonRef().key() == muonRef.key()) {
0892 reco::PFCandidateRef pfRef(pfMuons, index);
0893 aMuon.setPFCandidateRef(pfRef);
0894 if (embedPfEcalEnergy_)
0895 aMuon.setPfEcalEnergy(pfmu.ecalEnergy());
0896 if (embedPFCandidate_)
0897 aMuon.embedPFCandidate();
0898 break;
0899 }
0900 }
0901 index++;
0902 }
0903 }
0904
0905 if (addInverseBeta_) {
0906 aMuon.readTimeExtra((*muonsTimeExtra)[muonRef]);
0907 }
0908
0909 aMuon.initSimInfo();
0910 if (simInfoIsAvailalbe) {
0911 const auto& msi = (*simInfo)[muonBaseRef];
0912 aMuon.setSimType(msi.primaryClass);
0913 aMuon.setExtSimType(msi.extendedClass);
0914 aMuon.setSimFlavour(msi.flavour);
0915 aMuon.setSimHeaviestMotherFlavour(msi.heaviestMotherFlavour);
0916 aMuon.setSimPdgId(msi.pdgId);
0917 aMuon.setSimMotherPdgId(msi.motherPdgId);
0918 aMuon.setSimBX(msi.tpBX);
0919 aMuon.setSimTpEvent(msi.tpEvent);
0920 aMuon.setSimProdRho(msi.vertex.Rho());
0921 aMuon.setSimProdZ(msi.vertex.Z());
0922 aMuon.setSimPt(msi.p4.pt());
0923 aMuon.setSimEta(msi.p4.eta());
0924 aMuon.setSimPhi(msi.p4.phi());
0925 aMuon.setSimMatchQuality(msi.tpAssoQuality);
0926 }
0927 patMuons.push_back(aMuon);
0928 }
0929 }
0930
0931
0932 std::sort(patMuons.begin(), patMuons.end(), [](auto const& t1, auto const& t2) { return t1.pt() > t2.pt(); });
0933
0934
0935
0936
0937
0938 edm::Handle<double> rho;
0939 if (computeMiniIso_)
0940 iEvent.getByToken(rho_, rho);
0941 const reco::Vertex* pv(nullptr);
0942 if (primaryVertexIsValid)
0943 pv = &primaryVertex;
0944
0945 edm::Handle<std::vector<pat::TriggerObjectStandAlone>> triggerObjects;
0946 edm::Handle<edm::TriggerResults> triggerResults;
0947 bool triggerObjectsAvailable = false;
0948 bool triggerResultsAvailable = false;
0949 if (addTriggerMatching_) {
0950 triggerObjectsAvailable = iEvent.getByToken(triggerObjects_, triggerObjects);
0951 triggerResultsAvailable = iEvent.getByToken(triggerResults_, triggerResults);
0952 }
0953
0954 for (auto& muon : patMuons) {
0955
0956 if (addTriggerMatching_ and triggerObjectsAvailable and triggerResultsAvailable) {
0957 const edm::TriggerNames& triggerNames(iEvent.triggerNames(*triggerResults));
0958 fillL1TriggerInfo(muon, triggerObjects, triggerNames, geometry);
0959 fillHltTriggerInfo(muon, triggerObjects, triggerNames, hltCollectionFilters_);
0960 }
0961
0962 if (recomputeBasicSelectors_) {
0963 muon.setSelectors(0);
0964 bool isRun2016BCDEF = (272728 <= iEvent.run() && iEvent.run() <= 278808);
0965 muon.setSelectors(muon::makeSelectorBitset(muon, pv, isRun2016BCDEF));
0966 }
0967 float miniIsoValue = -1;
0968 if (computeMiniIso_) {
0969
0970
0971 miniIsoValue = getRelMiniIsoPUCorrected(muon, *rho, effectiveAreaVec_);
0972
0973 muon.setSelector(reco::Muon::MiniIsoLoose, miniIsoValue < 0.40);
0974 muon.setSelector(reco::Muon::MiniIsoMedium, miniIsoValue < 0.20);
0975 muon.setSelector(reco::Muon::MiniIsoTight, miniIsoValue < 0.10);
0976 muon.setSelector(reco::Muon::MiniIsoVeryTight, miniIsoValue < 0.05);
0977 }
0978
0979 double puppiCombinedIsolationPAT = -1;
0980 if (computePuppiCombinedIso_) {
0981 puppiCombinedIsolationPAT = puppiCombinedIsolation(muon, pc.product());
0982 muon.setSelector(reco::Muon::PuppiIsoLoose, puppiCombinedIsolationPAT < 0.27);
0983 muon.setSelector(reco::Muon::PuppiIsoMedium, puppiCombinedIsolationPAT < 0.22);
0984 muon.setSelector(reco::Muon::PuppiIsoTight, puppiCombinedIsolationPAT < 0.12);
0985 }
0986
0987 std::array<float, 2> jetPtRatioRel = {{0.0, 0.0}};
0988 if (primaryVertexIsValid && computeMiniIso_) {
0989 if (useJec_) {
0990 jetPtRatioRel = globalCache()->calculatePtRatioRel().computePtRatioRel(
0991 muon, *(mvaBTagCollectionTag.product()), mvaL1Corrector.product(), mvaL1L2L3ResCorrector.product());
0992 } else {
0993 jetPtRatioRel = globalCache()->calculatePtRatioRel().computePtRatioRel(muon, *mvaBTagCollectionTag);
0994 }
0995
0996 muon.setJetPtRatio(jetPtRatioRel[0]);
0997 muon.setJetPtRel(jetPtRatioRel[1]);
0998
0999
1000 if (computeMiniIso_) {
1001 muon.setSelector(reco::Muon::MultiIsoMedium,
1002 miniIsoValue < 0.11 && (muon.jetPtRatio() > 0.74 || muon.jetPtRel() > 6.8));
1003 }
1004 }
1005
1006
1007 float mvaID = 0.0;
1008 constexpr int MVAsentinelValue = -99;
1009 if (computeMuonIDMVA_) {
1010 if (muon.isLooseMuon()) {
1011 mvaID = globalCache()->muonMvaIDEstimator().computeMVAID(muon)[1];
1012 } else {
1013 mvaID = MVAsentinelValue;
1014 }
1015 muon.setMvaIDValue(mvaID);
1016 muon.setSelector(reco::Muon::MvaIDwpMedium, muon.mvaIDValue() > mvaIDmediumCut_);
1017 muon.setSelector(reco::Muon::MvaIDwpTight, muon.mvaIDValue() > mvaIDtightCut_);
1018 }
1019
1020
1021 if (computeSoftMuonMVA_) {
1022 float mva = globalCache()->softMuonMvaEstimator().computeMva(muon);
1023 muon.setSoftMvaValue(mva);
1024
1025 muon.setSelector(reco::Muon::SoftMvaId, muon.softMvaValue() > 0.58);
1026
1027
1028 muon.setSoftMvaRun3Value(computeSoftMvaRun3(*softMuonMvaRun3Booster_, muon));
1029 }
1030 }
1031
1032
1033 iEvent.emplace(patMuonPutToken_, std::move(patMuons));
1034
1035 if (isolator_.enabled())
1036 isolator_.endEvent();
1037 }
1038
1039 void PATMuonProducer::fillMuon(Muon& aMuon,
1040 const MuonBaseRef& muonRef,
1041 const reco::CandidateBaseRef& baseRef,
1042 const GenAssociations& genMatches,
1043 const IsoDepositMaps& deposits,
1044 const IsolationValueMaps& isolationValues) const {
1045
1046
1047
1048
1049
1050 if (useParticleFlow_)
1051 aMuon.setP4(aMuon.pfCandidateRef()->p4());
1052 if (embedTrack_)
1053 aMuon.embedTrack();
1054 if (embedStandAloneMuon_)
1055 aMuon.embedStandAloneMuon();
1056 if (embedCombinedMuon_)
1057 aMuon.embedCombinedMuon();
1058
1059
1060 if (aMuon.isGlobalMuon()) {
1061 if (embedPickyMuon_ && aMuon.isAValidMuonTrack(reco::Muon::Picky))
1062 aMuon.embedPickyMuon();
1063 if (embedTpfmsMuon_ && aMuon.isAValidMuonTrack(reco::Muon::TPFMS))
1064 aMuon.embedTpfmsMuon();
1065 if (embedDytMuon_ && aMuon.isAValidMuonTrack(reco::Muon::DYT))
1066 aMuon.embedDytMuon();
1067 }
1068
1069
1070 if (embedBestTrack_)
1071 aMuon.embedMuonBestTrack(forceEmbedBestTrack_);
1072 if (embedTunePBestTrack_)
1073 aMuon.embedTunePMuonBestTrack(forceEmbedBestTrack_);
1074
1075
1076 if (addGenMatch_) {
1077 for (auto const& genMatch : genMatches) {
1078 reco::GenParticleRef genMuon = (*genMatch)[baseRef];
1079 aMuon.addGenParticleRef(genMuon);
1080 }
1081 if (embedGenMatch_)
1082 aMuon.embedGenParticle();
1083 }
1084 if (efficiencyLoader_.enabled()) {
1085 efficiencyLoader_.setEfficiencies(aMuon, muonRef);
1086 }
1087
1088 for (size_t j = 0, nd = deposits.size(); j < nd; ++j) {
1089 if (useParticleFlow_) {
1090 if (deposits[j]->contains(baseRef.id())) {
1091 aMuon.setIsoDeposit(isoDepositLabels_[j].first, (*deposits[j])[baseRef]);
1092 } else if (deposits[j]->contains(muonRef.id())) {
1093 aMuon.setIsoDeposit(isoDepositLabels_[j].first, (*deposits[j])[muonRef]);
1094 } else {
1095 reco::CandidatePtr source = aMuon.pfCandidateRef()->sourceCandidatePtr(0);
1096 aMuon.setIsoDeposit(isoDepositLabels_[j].first, (*deposits[j])[source]);
1097 }
1098 } else {
1099 aMuon.setIsoDeposit(isoDepositLabels_[j].first, (*deposits[j])[muonRef]);
1100 }
1101 }
1102
1103 for (size_t j = 0; j < isolationValues.size(); ++j) {
1104 if (useParticleFlow_) {
1105 if (isolationValues[j]->contains(baseRef.id())) {
1106 aMuon.setIsolation(isolationValueLabels_[j].first, (*isolationValues[j])[baseRef]);
1107 } else if (isolationValues[j]->contains(muonRef.id())) {
1108 aMuon.setIsolation(isolationValueLabels_[j].first, (*isolationValues[j])[muonRef]);
1109 } else {
1110 reco::CandidatePtr source = aMuon.pfCandidateRef()->sourceCandidatePtr(0);
1111 aMuon.setIsolation(isolationValueLabels_[j].first, (*isolationValues[j])[source]);
1112 }
1113 } else {
1114 aMuon.setIsolation(isolationValueLabels_[j].first, (*isolationValues[j])[muonRef]);
1115 }
1116 }
1117
1118 if (resolutionLoader_.enabled()) {
1119 resolutionLoader_.setResolutions(aMuon);
1120 }
1121 }
1122
1123 void PATMuonProducer::setMuonMiniIso(Muon& aMuon, const PackedCandidateCollection* pc) {
1124 pat::PFIsolation miniiso = pat::getMiniPFIsolation(pc,
1125 aMuon.polarP4(),
1126 miniIsoParams_[0],
1127 miniIsoParams_[1],
1128 miniIsoParams_[2],
1129 miniIsoParams_[3],
1130 miniIsoParams_[4],
1131 miniIsoParams_[5],
1132 miniIsoParams_[6],
1133 miniIsoParams_[7],
1134 miniIsoParams_[8]);
1135 aMuon.setMiniPFIsolation(miniiso);
1136 }
1137
1138 double PATMuonProducer::getRelMiniIsoPUCorrected(const pat::Muon& muon, double rho, const std::vector<double>& area) {
1139 double mindr(miniIsoParams_[0]);
1140 double maxdr(miniIsoParams_[1]);
1141 double kt_scale(miniIsoParams_[2]);
1142 double drcut = pat::miniIsoDr(muon.polarP4(), mindr, maxdr, kt_scale);
1143 return pat::muonRelMiniIsoPUCorrected(muon.miniPFIsolation(), muon.polarP4(), drcut, rho, area);
1144 }
1145
1146 double PATMuonProducer::puppiCombinedIsolation(const pat::Muon& muon, const pat::PackedCandidateCollection* pc) {
1147 constexpr double dR_threshold = 0.4;
1148 constexpr double dR2_threshold = dR_threshold * dR_threshold;
1149 constexpr double mix_fraction = 0.5;
1150 enum particleType { CH = 0, NH = 1, PH = 2, OTHER = 100000 };
1151 double val_PuppiWithLep = 0.0;
1152 double val_PuppiWithoutLep = 0.0;
1153
1154 for (const auto& cand : *pc) {
1155
1156 const particleType pType = isChargedHadron(cand.pdgId()) ? CH
1157 : isNeutralHadron(cand.pdgId()) ? NH
1158 : isPhoton(cand.pdgId()) ? PH
1159 : OTHER;
1160 if (pType == OTHER) {
1161 if (cand.pdgId() != 1 && cand.pdgId() != 2 && abs(cand.pdgId()) != 11 && abs(cand.pdgId()) != 13) {
1162 LogTrace("PATMuonProducer") << "candidate with PDGID = " << cand.pdgId()
1163 << " is not CH/NH/PH/e/mu or 1/2 (and this is removed from isolation calculation)"
1164 << std::endl;
1165 }
1166 continue;
1167 }
1168 double d_eta = std::abs(cand.eta() - muon.eta());
1169 if (d_eta > dR_threshold)
1170 continue;
1171
1172 double d_phi = std::abs(reco::deltaPhi(cand.phi(), muon.phi()));
1173 if (d_phi > dR_threshold)
1174 continue;
1175
1176 double dR2 = reco::deltaR2(cand, muon);
1177 if (dR2 > dR2_threshold)
1178 continue;
1179 if (pType == CH && dR2 < 0.0001 * 0.0001)
1180 continue;
1181 if (pType == NH && dR2 < 0.01 * 0.01)
1182 continue;
1183 if (pType == PH && dR2 < 0.01 * 0.01)
1184 continue;
1185 val_PuppiWithLep += cand.pt() * cand.puppiWeight();
1186 val_PuppiWithoutLep += cand.pt() * cand.puppiWeightNoLep();
1187
1188 }
1189
1190 double reliso_Puppi_withLep = val_PuppiWithLep / muon.pt();
1191 double reliso_Puppi_withoutlep = val_PuppiWithoutLep / muon.pt();
1192 double reliso_Puppi_combined = mix_fraction * reliso_Puppi_withLep + (1.0 - mix_fraction) * reliso_Puppi_withoutlep;
1193 return reliso_Puppi_combined;
1194 }
1195
1196 bool PATMuonProducer::isNeutralHadron(long pdgid) { return std::abs(pdgid) == 130; }
1197
1198 bool PATMuonProducer::isChargedHadron(long pdgid) { return std::abs(pdgid) == 211; }
1199
1200 bool PATMuonProducer::isPhoton(long pdgid) { return pdgid == 22; }
1201
1202
1203 void PATMuonProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1204 edm::ParameterSetDescription iDesc;
1205 iDesc.setComment("PAT muon producer module");
1206
1207
1208 iDesc.add<edm::InputTag>("muonSource", edm::InputTag("no default"))->setComment("input collection");
1209
1210
1211 iDesc.add<bool>("embedMuonBestTrack", true)->setComment("embed muon best track (global pflow)");
1212 iDesc.add<bool>("embedTunePMuonBestTrack", true)->setComment("embed muon best track (muon only)");
1213 iDesc.add<bool>("forceBestTrackEmbedding", true)
1214 ->setComment(
1215 "force embedding separately the best tracks even if they're already embedded e.g. as tracker or global "
1216 "tracks");
1217 iDesc.add<bool>("embedTrack", true)->setComment("embed external track");
1218 iDesc.add<bool>("embedStandAloneMuon", true)->setComment("embed external stand-alone muon");
1219 iDesc.add<bool>("embedCombinedMuon", false)->setComment("embed external combined muon");
1220 iDesc.add<bool>("embedPickyMuon", false)->setComment("embed external picky track");
1221 iDesc.add<bool>("embedTpfmsMuon", false)->setComment("embed external tpfms track");
1222 iDesc.add<bool>("embedDytMuon", false)->setComment("embed external dyt track ");
1223
1224
1225 iDesc.add<bool>("embedCaloMETMuonCorrs", true)->setComment("whether to add MET muon correction for caloMET or not");
1226 iDesc.add<edm::InputTag>("caloMETMuonCorrs", edm::InputTag("muonMETValueMapProducer", "muCorrData"))
1227 ->setComment("source of MET muon corrections for caloMET");
1228 iDesc.add<bool>("embedTcMETMuonCorrs", true)->setComment("whether to add MET muon correction for tcMET or not");
1229 iDesc.add<edm::InputTag>("tcMETMuonCorrs", edm::InputTag("muonTCMETValueMapProducer", "muCorrData"))
1230 ->setComment("source of MET muon corrections for tcMET");
1231
1232
1233 iDesc.add<edm::InputTag>("pfMuonSource", edm::InputTag("pfMuons"))->setComment("particle flow input collection");
1234 iDesc.add<bool>("useParticleFlow", false)->setComment("whether to use particle flow or not");
1235 iDesc.add<bool>("embedPFCandidate", false)->setComment("embed external particle flow object");
1236 iDesc.add<bool>("embedPfEcalEnergy", true)->setComment("add ecal energy as reconstructed by PF");
1237
1238
1239 iDesc.add<bool>("addInverseBeta", true)->setComment("add combined inverse beta");
1240 iDesc.add<edm::InputTag>("sourceInverseBeta", edm::InputTag("muons", "combined"))
1241 ->setComment("source of inverse beta values");
1242
1243
1244 iDesc.add<bool>("addGenMatch", true)->setComment("add MC matching");
1245 iDesc.add<bool>("embedGenMatch", false)->setComment("embed MC matched MC information");
1246 std::vector<edm::InputTag> emptySourceVector;
1247 iDesc
1248 .addNode(edm::ParameterDescription<edm::InputTag>("genParticleMatch", edm::InputTag(), true) xor
1249 edm::ParameterDescription<std::vector<edm::InputTag>>("genParticleMatch", emptySourceVector, true))
1250 ->setComment("input with MC match information");
1251
1252
1253 iDesc.add<bool>("computeMiniIso", false)->setComment("whether or not to compute and store electron mini-isolation");
1254 iDesc.add<bool>("computePuppiCombinedIso", false)
1255 ->setComment("whether or not to compute and store puppi combined isolation");
1256
1257 iDesc.add<edm::InputTag>("pfCandsForMiniIso", edm::InputTag("packedPFCandidates"))
1258 ->setComment("collection to use to compute mini-iso");
1259 iDesc.add<std::vector<double>>("miniIsoParams", std::vector<double>())
1260 ->setComment("mini-iso parameters to use for muons");
1261
1262 iDesc.add<bool>("addTriggerMatching", false)->setComment("add L1 and HLT matching to offline muon");
1263
1264 pat::helper::KinResolutionsLoader::fillDescription(iDesc);
1265
1266
1267 edm::ParameterSetDescription isoDepositsPSet;
1268 isoDepositsPSet.addOptional<edm::InputTag>("tracker");
1269 isoDepositsPSet.addOptional<edm::InputTag>("ecal");
1270 isoDepositsPSet.addOptional<edm::InputTag>("hcal");
1271 isoDepositsPSet.addOptional<edm::InputTag>("particle");
1272 isoDepositsPSet.addOptional<edm::InputTag>("pfChargedHadrons");
1273 isoDepositsPSet.addOptional<edm::InputTag>("pfChargedAll");
1274 isoDepositsPSet.addOptional<edm::InputTag>("pfPUChargedHadrons");
1275 isoDepositsPSet.addOptional<edm::InputTag>("pfNeutralHadrons");
1276 isoDepositsPSet.addOptional<edm::InputTag>("pfPhotons");
1277 isoDepositsPSet.addOptional<std::vector<edm::InputTag>>("user");
1278 iDesc.addOptional("isoDeposits", isoDepositsPSet);
1279
1280
1281 edm::ParameterSetDescription isolationValuesPSet;
1282 isolationValuesPSet.addOptional<edm::InputTag>("tracker");
1283 isolationValuesPSet.addOptional<edm::InputTag>("ecal");
1284 isolationValuesPSet.addOptional<edm::InputTag>("hcal");
1285 isolationValuesPSet.addOptional<edm::InputTag>("particle");
1286 isolationValuesPSet.addOptional<edm::InputTag>("pfChargedHadrons");
1287 isolationValuesPSet.addOptional<edm::InputTag>("pfChargedAll");
1288 isolationValuesPSet.addOptional<edm::InputTag>("pfPUChargedHadrons");
1289 isolationValuesPSet.addOptional<edm::InputTag>("pfNeutralHadrons");
1290 isolationValuesPSet.addOptional<edm::InputTag>("pfPhotons");
1291 iDesc.addOptional("isolationValues", isolationValuesPSet);
1292
1293 iDesc.ifValue(edm::ParameterDescription<bool>("addPuppiIsolation", false, true),
1294 true >> (edm::ParameterDescription<edm::InputTag>(
1295 "puppiIsolationChargedHadrons",
1296 edm::InputTag("muonPUPPIIsolation", "h+-DR030-ThresholdVeto000-ConeVeto000"),
1297 true) and
1298 edm::ParameterDescription<edm::InputTag>(
1299 "puppiIsolationNeutralHadrons",
1300 edm::InputTag("muonPUPPIIsolation", "h0-DR030-ThresholdVeto000-ConeVeto001"),
1301 true) and
1302 edm::ParameterDescription<edm::InputTag>(
1303 "puppiIsolationPhotons",
1304 edm::InputTag("muonPUPPIIsolation", "gamma-DR030-ThresholdVeto000-ConeVeto001"),
1305 true) and
1306 edm::ParameterDescription<edm::InputTag>(
1307 "puppiNoLeptonsIsolationChargedHadrons",
1308 edm::InputTag("muonPUPPINoLeptonsIsolation", "h+-DR030-ThresholdVeto000-ConeVeto000"),
1309 true) and
1310 edm::ParameterDescription<edm::InputTag>(
1311 "puppiNoLeptonsIsolationNeutralHadrons",
1312 edm::InputTag("muonPUPPINoLeptonsIsolation", "h0-DR030-ThresholdVeto000-ConeVeto001"),
1313 true) and
1314 edm::ParameterDescription<edm::InputTag>(
1315 "puppiNoLeptonsIsolationPhotons",
1316 edm::InputTag("muonPUPPINoLeptonsIsolation", "gamma-DR030-ThresholdVeto000-ConeVeto001"),
1317 true)) or
1318 false >> edm::EmptyGroupDescription());
1319
1320
1321 edm::ParameterSetDescription efficienciesPSet;
1322 efficienciesPSet.setAllowAnything();
1323 iDesc.add("efficiencies", efficienciesPSet);
1324 iDesc.add<bool>("addEfficiencies", false);
1325
1326
1327 edm::ParameterSetDescription userDataPSet;
1328 PATUserDataHelper<Muon>::fillDescription(userDataPSet);
1329 iDesc.addOptional("userData", userDataPSet);
1330
1331 edm::ParameterSetDescription isolationPSet;
1332 isolationPSet.setAllowAnything();
1333 iDesc.add("userIsolation", isolationPSet);
1334
1335 iDesc.add<bool>("embedHighLevelSelection", true)->setComment("embed high level selection");
1336 edm::ParameterSetDescription highLevelPSet;
1337 highLevelPSet.setAllowAnything();
1338 iDesc.addNode(edm::ParameterDescription<edm::InputTag>("beamLineSrc", edm::InputTag(), true))
1339 ->setComment("input with high level selection");
1340 iDesc.addNode(edm::ParameterDescription<edm::InputTag>("pvSrc", edm::InputTag(), true))
1341 ->setComment("input with high level selection");
1342
1343
1344 }
1345
1346
1347
1348 void PATMuonProducer::embedHighLevel(pat::Muon& aMuon,
1349 reco::TrackRef track,
1350 reco::TransientTrack& tt,
1351 reco::Vertex& primaryVertex,
1352 bool primaryVertexIsValid,
1353 reco::BeamSpot& beamspot,
1354 bool beamspotIsValid) {
1355
1356
1357
1358 aMuon.setDB(track->dxy(primaryVertex.position()),
1359 track->dxyError(primaryVertex.position(), primaryVertex.covariance()),
1360 pat::Muon::PV2D);
1361
1362
1363 std::pair<bool, Measurement1D> result =
1364 IPTools::signedImpactParameter3D(tt, GlobalVector(track->px(), track->py(), track->pz()), primaryVertex);
1365 double d0_corr = result.second.value();
1366 double d0_err = primaryVertexIsValid ? result.second.error() : -1.0;
1367 aMuon.setDB(d0_corr, d0_err, pat::Muon::PV3D);
1368
1369
1370
1371
1372 aMuon.setDB(track->dxy(beamspot), track->dxyError(beamspot), pat::Muon::BS2D);
1373
1374
1375 reco::Vertex vBeamspot(beamspot.position(), beamspot.rotatedCovariance3D());
1376
1377
1378 result = IPTools::signedImpactParameter3D(tt, GlobalVector(track->px(), track->py(), track->pz()), vBeamspot);
1379 d0_corr = result.second.value();
1380 d0_err = beamspotIsValid ? result.second.error() : -1.0;
1381 aMuon.setDB(d0_corr, d0_err, pat::Muon::BS3D);
1382
1383
1384 aMuon.setDB(
1385 track->dz(primaryVertex.position()), std::hypot(track->dzError(), primaryVertex.zError()), pat::Muon::PVDZ);
1386 }
1387
1388 #include "FWCore/Framework/interface/MakerMacros.h"
1389
1390 DEFINE_FWK_MODULE(PATMuonProducer);