File indexing completed on 2024-04-06 12:09:36
0001
0002
0003
0004
0005
0006
0007
0008 #include "DQMOffline/L1Trigger/interface/L1TPhase2MuonOffline.h"
0009
0010
0011 #include "DataFormats/L1TMuonPhase2/interface/Constants.h"
0012
0013 using namespace reco;
0014 using namespace edm;
0015 using namespace std;
0016 using namespace l1t;
0017
0018
0019 GenMuonGMTPair::GenMuonGMTPair(const reco::GenParticle* muon, const l1t::L1Candidate* gmtmu)
0020 : mu_(muon), gmtmu_(gmtmu) {
0021 if (gmtmu) {
0022 gmtEta_ = gmtmu_->eta();
0023 gmtPhi_ = gmtmu_->phi();
0024 } else {
0025 gmtEta_ = -5.;
0026 gmtPhi_ = -5.;
0027 }
0028 if (mu_) {
0029 muEta_ = mu_->eta();
0030 muPhi_ = mu_->phi();
0031 } else {
0032 muEta_ = 999.;
0033 muPhi_ = 999.;
0034 }
0035 };
0036
0037 GenMuonGMTPair::GenMuonGMTPair(const GenMuonGMTPair& muonGmtPair) {
0038 mu_ = muonGmtPair.mu_;
0039 gmtmu_ = muonGmtPair.gmtmu_;
0040
0041 gmtEta_ = muonGmtPair.gmtEta_;
0042 gmtPhi_ = muonGmtPair.gmtPhi_;
0043
0044 muEta_ = muonGmtPair.muEta_;
0045 muPhi_ = muonGmtPair.muPhi_;
0046 }
0047
0048 float GenMuonGMTPair::dR2() {
0049 if (!gmtmu_)
0050 return 999.;
0051 float dEta = gmtEta_ - muEta_;
0052 float dPhi = reco::deltaPhi(gmtPhi_, muPhi_);
0053 return dEta * dEta + dPhi * dPhi;
0054 }
0055
0056 L1TPhase2MuonOffline::EtaRegion GenMuonGMTPair::etaRegion() const {
0057 if (std::abs(muEta_) < 0.83)
0058 return L1TPhase2MuonOffline::kEtaRegionBmtf;
0059 if (std::abs(muEta_) < 1.24)
0060 return L1TPhase2MuonOffline::kEtaRegionOmtf;
0061 if (std::abs(muEta_) < 2.4)
0062 return L1TPhase2MuonOffline::kEtaRegionEmtf;
0063 return L1TPhase2MuonOffline::kEtaRegionAll;
0064 }
0065
0066 double GenMuonGMTPair::getDeltaVar(const L1TPhase2MuonOffline::ResType type) const {
0067 if (type == L1TPhase2MuonOffline::kResPt)
0068 return (gmtPt() - pt()) / pt();
0069 if (type == L1TPhase2MuonOffline::kRes1OverPt)
0070 return (pt() - gmtPt()) / gmtPt();
0071 if (type == L1TPhase2MuonOffline::kResQOverPt)
0072 return (gmtCharge() * charge() * pt() - gmtPt()) /
0073 gmtPt();
0074 if (type == L1TPhase2MuonOffline::kResPhi)
0075 return reco::deltaPhi(gmtPhi(), muPhi_);
0076 if (type == L1TPhase2MuonOffline::kResEta)
0077 return gmtEta() - muEta_;
0078 if (type == L1TPhase2MuonOffline::kResCh)
0079 return gmtCharge() - charge();
0080 return -999.;
0081 }
0082
0083 double GenMuonGMTPair::getVar(const L1TPhase2MuonOffline::EffType type) const {
0084 if (type == L1TPhase2MuonOffline::kEffPt)
0085 return pt();
0086 if (type == L1TPhase2MuonOffline::kEffPhi)
0087 return muPhi_;
0088 if (type == L1TPhase2MuonOffline::kEffEta)
0089 return muEta_;
0090 return -999.;
0091 }
0092
0093
0094 L1TPhase2MuonOffline::L1TPhase2MuonOffline(const ParameterSet& ps)
0095 : gmtMuonToken_(consumes<l1t::SAMuonCollection>(ps.getParameter<edm::InputTag>("gmtMuonToken"))),
0096 gmtTkMuonToken_(consumes<l1t::TrackerMuonCollection>(ps.getParameter<edm::InputTag>("gmtTkMuonToken"))),
0097 genParticleToken_(
0098 consumes<std::vector<reco::GenParticle>>(ps.getUntrackedParameter<edm::InputTag>("genParticlesInputTag"))),
0099 muonTypes_({kSAMuon, kTkMuon}),
0100 effTypes_({kEffPt, kEffPhi, kEffEta}),
0101 resTypes_({kResPt, kResQOverPt, kResPhi, kResEta}),
0102 etaRegions_({kEtaRegionAll, kEtaRegionBmtf, kEtaRegionOmtf, kEtaRegionEmtf}),
0103 qualLevels_({kQualOpen, kQualDouble, kQualSingle}),
0104 resNames_({{kResPt, "pt"},
0105 {kRes1OverPt, "1overpt"},
0106 {kResQOverPt, "qoverpt"},
0107 {kResPhi, "phi"},
0108 {kResEta, "eta"},
0109 {kResCh, "charge"}}),
0110 resLabels_({{kResPt, "(p_{T}^{L1} - p_{T}^{reco})/p_{T}^{reco}"},
0111 {kRes1OverPt, "(p_{T}^{reco} - p_{T}^{L1})/p_{T}^{L1}"},
0112 {kResQOverPt, "(q^{L1}*q^{reco}*p_{T}^{reco} - p_{T}^{L1})/p_{T}^{L1}"},
0113 {kResPhi, "#phi_{L1} - #phi_{reco}"},
0114 {kResEta, "#eta_{L1} - #eta_{reco}"},
0115 {kResCh, "charge^{L1} - charge^{reco}"}}),
0116 etaNames_({{kEtaRegionAll, "etaMin0_etaMax2p4"},
0117 {kEtaRegionBmtf, "etaMin0_etaMax0p83"},
0118 {kEtaRegionOmtf, "etaMin0p83_etaMax1p24"},
0119 {kEtaRegionEmtf, "etaMin1p24_etaMax2p4"}}),
0120 qualNames_({{kQualOpen, "qualOpen"}, {kQualDouble, "qualDouble"}, {kQualSingle, "qualSingle"}}),
0121 muonNames_({{kSAMuon, "SAMuon"}, {kTkMuon, "TkMuon"}}),
0122 histFolder_(ps.getUntrackedParameter<string>("histFolder")),
0123 cutsVPSet_(ps.getUntrackedParameter<std::vector<edm::ParameterSet>>("cuts")),
0124 effVsPtBins_(ps.getUntrackedParameter<std::vector<double>>("efficiencyVsPtBins")),
0125 effVsPhiBins_(ps.getUntrackedParameter<std::vector<double>>("efficiencyVsPhiBins")),
0126 effVsEtaBins_(ps.getUntrackedParameter<std::vector<double>>("efficiencyVsEtaBins")),
0127 maxGmtMuonDR_(ps.getUntrackedParameter<double>("maxDR")) {
0128 edm::LogInfo("L1TPhase2MuonOffline") << "L1TPhase2MuonOffline::L1TPhase2MuonOffline()" << endl;
0129
0130 for (const auto& c : cutsVPSet_) {
0131 const auto qCut = c.getUntrackedParameter<int>("qualCut");
0132 QualLevel qLevel = kQualOpen;
0133 if (qCut > 11) {
0134 qLevel = kQualSingle;
0135 } else if (qCut > 7) {
0136 qLevel = kQualDouble;
0137 } else if (qCut > 3) {
0138 qLevel = kQualOpen;
0139 }
0140 cuts_.emplace_back(std::make_pair(c.getUntrackedParameter<int>("ptCut"), qLevel));
0141 }
0142 }
0143
0144
0145 void L1TPhase2MuonOffline::dqmBeginRun(const edm::Run& run, const edm::EventSetup& iSetup) {
0146 edm::LogInfo("L1TPhase2MuonOFfline") << "L1TPhase2MuonOffline::dqmBeginRun" << endl;
0147 }
0148
0149
0150 void L1TPhase2MuonOffline::bookHistograms(DQMStore::IBooker& ibooker,
0151 const edm::Run& run,
0152 const edm::EventSetup& iSetup) {
0153 edm::LogInfo("L1TPhase2MuonOFfline") << "L1TPhase2MuonOffline::bookHistograms" << endl;
0154
0155
0156 for (const auto mutype : muonTypes_) {
0157 bookControlHistos(ibooker, mutype);
0158 bookEfficiencyHistos(ibooker, mutype);
0159 bookResolutionHistos(ibooker, mutype);
0160 }
0161 }
0162
0163
0164 void L1TPhase2MuonOffline::analyze(const Event& iEvent, const EventSetup& eventSetup) {
0165 edm::LogInfo("L1TPhase2MuonOffline") << "L1TPhase2MuonOffline::analyze() " << endl;
0166
0167
0168 iEvent.getByToken(genParticleToken_, genparticles_);
0169
0170 std::vector<const reco::GenParticle*> genmus;
0171 for (const reco::GenParticle& gen : *genparticles_) {
0172 if (std::abs(gen.pdgId()) != 13)
0173 continue;
0174 genmus.push_back(&gen);
0175 }
0176 edm::LogInfo("L1TPhase2MuonOffline") << "L1TPhase2MuonOffline::analyze() N of genmus: " << genmus.size() << endl;
0177
0178
0179 iEvent.getByToken(gmtMuonToken_, gmtSAMuon_);
0180 iEvent.getByToken(gmtTkMuonToken_, gmtTkMuon_);
0181
0182
0183 edm::LogInfo("L1TPhase2MuonOffline") << "Fill Control histograms for GMT Muons" << endl;
0184 fillControlHistos();
0185
0186
0187 if (genmus.empty())
0188 return;
0189 edm::LogInfo("L1TPhase2MuonOffline") << "L1TPhase2MuonOffline::analyze() calling matchMuonsToGen() " << endl;
0190 matchMuonsToGen(genmus);
0191
0192
0193 fillEfficiencyHistos();
0194 fillResolutionHistos();
0195 edm::LogInfo("L1TPhase2MuonOffline") << "L1TPhase2MuonOffline::analyze() Computation finished" << endl;
0196 }
0197
0198
0199 void L1TPhase2MuonOffline::bookControlHistos(DQMStore::IBooker& ibooker, MuType mutype) {
0200 edm::LogInfo("L1TPhase2MuonOffline") << "L1TPhase2MuonOffline::bookControlHistos()" << endl;
0201
0202 ibooker.setCurrentFolder(histFolder_ + "/" + muonNames_[mutype] + "/control_variables");
0203
0204 controlHistos_[mutype][kPt] = ibooker.book1D(muonNames_[mutype] + "Pt", "MuonPt; p_{T}", 50, 0., 100.);
0205 controlHistos_[mutype][kPhi] = ibooker.book1D(muonNames_[mutype] + "Phi", "MuonPhi; #phi", 66, -3.3, 3.3);
0206 controlHistos_[mutype][kEta] = ibooker.book1D(muonNames_[mutype] + "Eta", "MuonEta; #eta", 50, -2.5, 2.5);
0207 controlHistos_[mutype][kIso] = ibooker.book1D(muonNames_[mutype] + "Iso", "MuonIso; RelIso", 50, 0, 1.0);
0208 controlHistos_[mutype][kQual] = ibooker.book1D(muonNames_[mutype] + "Qual", "MuonQual; Quality", 15, 0.5, 15.5);
0209 controlHistos_[mutype][kZ0] = ibooker.book1D(muonNames_[mutype] + "Z0", "MuonZ0; Z_{0}", 50, 0, 50.0);
0210 controlHistos_[mutype][kD0] = ibooker.book1D(muonNames_[mutype] + "D0", "MuonD0; D_{0}", 50, 0, 200.);
0211 }
0212
0213 void L1TPhase2MuonOffline::bookEfficiencyHistos(DQMStore::IBooker& ibooker, MuType mutype) {
0214 edm::LogInfo("L1TPhase2MuonOffline") << "L1TPhase2MuonOffline::bookEfficiencyHistos()" << endl;
0215
0216 ibooker.setCurrentFolder(histFolder_ + "/" + muonNames_[mutype] + "/nums_and_dens");
0217
0218 std::string histoname = "";
0219 for (const auto eta : etaRegions_) {
0220 for (const auto q : qualLevels_) {
0221 histoname = "Eff_" + muonNames_[mutype] + "_" + etaNames_[eta] + "_" + qualNames_[q];
0222
0223 auto histBins = getHistBinsEff(kEffPt);
0224 efficiencyNum_[mutype][eta][q][kEffPt] =
0225 ibooker.book1D(histoname + "_Pt_Num", "MuonPt; p_{T} ;", histBins.size() - 1, &histBins[0]);
0226 efficiencyDen_[mutype][eta][q][kEffPt] =
0227 ibooker.book1D(histoname + "_Pt_Den", "MuonPt; p_{T} ;", histBins.size() - 1, &histBins[0]);
0228
0229 histBins = getHistBinsEff(kEffEta);
0230 efficiencyNum_[mutype][eta][q][kEffEta] =
0231 ibooker.book1D(histoname + "_Eta_Num", "MuonEta; #eta ;", histBins.size() - 1, &histBins[0]);
0232 efficiencyDen_[mutype][eta][q][kEffEta] =
0233 ibooker.book1D(histoname + "_Eta_Den", "MuonEta; #eta ;", histBins.size() - 1, &histBins[0]);
0234
0235 histBins = getHistBinsEff(kEffPhi);
0236 efficiencyNum_[mutype][eta][q][kEffPhi] =
0237 ibooker.book1D(histoname + "_Phi_Num", "MuonPhi; #phi ;", histBins.size() - 1, &histBins[0]);
0238 efficiencyDen_[mutype][eta][q][kEffPhi] =
0239 ibooker.book1D(histoname + "_Phi_Den", "MuonPhi; #phi ;", histBins.size() - 1, &histBins[0]);
0240 }
0241 }
0242 }
0243
0244 void L1TPhase2MuonOffline::bookResolutionHistos(DQMStore::IBooker& ibooker, MuType mutype) {
0245 edm::LogInfo("L1TPhase2MuonOffline") << "L1TPhase2MuonOffline::bookResolutionHistos()" << endl;
0246
0247 ibooker.setCurrentFolder(histFolder_ + "/" + muonNames_[mutype] + "/resolution");
0248 std::string histoname = "";
0249 for (const auto eta : etaRegions_) {
0250 for (const auto q : qualLevels_) {
0251 for (const auto var : resTypes_) {
0252 histoname = "Res_" + muonNames_[mutype] + "_" + etaNames_[eta] + "_" + qualNames_[q] + "_" + resNames_[var];
0253 auto nbins = std::get<0>(getHistBinsRes(var));
0254 auto xmin = std::get<1>(getHistBinsRes(var));
0255 auto xmax = std::get<2>(getHistBinsRes(var));
0256 resolutionHistos_[mutype][eta][q][var] =
0257 ibooker.book1D(histoname, resNames_[var] + ";" + resLabels_[var], nbins, xmin, xmax);
0258 }
0259 }
0260 }
0261 }
0262
0263
0264 void L1TPhase2MuonOffline::fillControlHistos() {
0265 for (auto& muIt : *gmtSAMuon_) {
0266 controlHistos_[kSAMuon][kPt]->Fill(lsb_pt * muIt.hwPt());
0267 controlHistos_[kSAMuon][kPhi]->Fill(lsb_phi * muIt.hwPhi());
0268 controlHistos_[kSAMuon][kEta]->Fill(lsb_eta * muIt.hwEta());
0269 controlHistos_[kSAMuon][kIso]->Fill(muIt.hwIso());
0270 controlHistos_[kSAMuon][kQual]->Fill(muIt.hwQual());
0271 controlHistos_[kSAMuon][kZ0]->Fill(lsb_z0 * muIt.hwZ0());
0272 controlHistos_[kSAMuon][kD0]->Fill(lsb_d0 * muIt.hwD0());
0273 }
0274
0275 for (auto& muIt : *gmtTkMuon_) {
0276 controlHistos_[kTkMuon][kPt]->Fill(lsb_pt * muIt.hwPt());
0277 controlHistos_[kTkMuon][kPhi]->Fill(lsb_phi * muIt.hwPhi());
0278 controlHistos_[kTkMuon][kEta]->Fill(lsb_eta * muIt.hwEta());
0279 controlHistos_[kTkMuon][kIso]->Fill(muIt.hwIso());
0280 controlHistos_[kTkMuon][kQual]->Fill(muIt.hwQual());
0281 controlHistos_[kTkMuon][kZ0]->Fill(lsb_z0 * muIt.hwZ0());
0282 controlHistos_[kTkMuon][kD0]->Fill(lsb_d0 * muIt.hwD0());
0283 }
0284 }
0285
0286 void L1TPhase2MuonOffline::fillEfficiencyHistos() {
0287 for (const auto& muIt : gmtSAMuonPairs_) {
0288 auto eta = muIt.etaRegion();
0289 for (const auto var : effTypes_) {
0290 auto varToFill = muIt.getVar(var);
0291 for (const auto& cut : cuts_) {
0292 const auto q = cut.second;
0293 efficiencyDen_[kSAMuon][eta][q][var]->Fill(varToFill);
0294 if (muIt.gmtPt() < 0)
0295 continue;
0296 if (muIt.gmtQual() < q * 4)
0297 continue;
0298 const auto gmtPtCut = cut.first;
0299 if (var != kEffPt && muIt.gmtPt() < gmtPtCut)
0300 continue;
0301 efficiencyNum_[kSAMuon][eta][q][var]->Fill(varToFill);
0302 }
0303 }
0304 }
0305
0306
0307 for (const auto& muIt : gmtTkMuonPairs_) {
0308 auto eta = muIt.etaRegion();
0309 for (const auto var : effTypes_) {
0310 auto varToFill = muIt.getVar(var);
0311 for (const auto& cut : cuts_) {
0312 const auto q = cut.second;
0313 efficiencyDen_[kTkMuon][eta][q][var]->Fill(varToFill);
0314 if (muIt.gmtPt() < 0)
0315 continue;
0316 if (muIt.gmtQual() < q * 4)
0317 continue;
0318 const auto gmtPtCut = cut.first;
0319 if (var != kEffPt && muIt.gmtPt() < gmtPtCut)
0320 continue;
0321 efficiencyNum_[kTkMuon][eta][q][var]->Fill(varToFill);
0322 }
0323 }
0324 }
0325 }
0326
0327 void L1TPhase2MuonOffline::fillResolutionHistos() {
0328 for (const auto& muIt : gmtSAMuonPairs_) {
0329 if (muIt.gmtPt() < 0)
0330 continue;
0331
0332 auto eta = muIt.etaRegion();
0333 for (const auto q : qualLevels_) {
0334 if (muIt.gmtQual() < q * 4)
0335 continue;
0336 for (const auto var : resTypes_) {
0337 auto varToFill = muIt.getDeltaVar(var);
0338
0339 resolutionHistos_[kSAMuon][eta][q][var]->Fill(varToFill);
0340 }
0341 }
0342 }
0343
0344 for (const auto& muIt : gmtTkMuonPairs_) {
0345 if (muIt.gmtPt() < 0)
0346 continue;
0347
0348 auto eta = muIt.etaRegion();
0349 for (const auto q : qualLevels_) {
0350 if (muIt.gmtQual() < q * 4)
0351 continue;
0352 for (const auto var : resTypes_) {
0353 auto varToFill = muIt.getDeltaVar(var);
0354
0355 resolutionHistos_[kTkMuon][eta][q][var]->Fill(varToFill);
0356 }
0357 }
0358 }
0359 }
0360
0361
0362 void L1TPhase2MuonOffline::matchMuonsToGen(std::vector<const reco::GenParticle*> genmus) {
0363 gmtSAMuonPairs_.clear();
0364 gmtTkMuonPairs_.clear();
0365
0366 edm::LogInfo("L1TPhase2MuonOffline") << "L1TPhase2MuonOffline::matchMuonsToGen() " << endl;
0367
0368 for (const reco::GenParticle* gen : genmus) {
0369 edm::LogInfo("L1TPhase2MuonOffline") << "Looping on genmus: " << gen << endl;
0370 GenMuonGMTPair pairBestCand(&(*gen), nullptr);
0371 float dr2Best = maxGmtMuonDR_ * maxGmtMuonDR_;
0372 for (auto& muIt : *gmtSAMuon_) {
0373 GenMuonGMTPair pairTmpCand(&(*gen), &(muIt));
0374 float dr2Tmp = pairTmpCand.dR2();
0375 if (dr2Tmp < dr2Best) {
0376 dr2Best = dr2Tmp;
0377 pairBestCand = pairTmpCand;
0378 }
0379 }
0380 gmtSAMuonPairs_.emplace_back(pairBestCand);
0381
0382 GenMuonGMTPair pairBestCand2(&(*gen), nullptr);
0383 dr2Best = maxGmtMuonDR_ * maxGmtMuonDR_;
0384 for (auto& tkmuIt : *gmtTkMuon_) {
0385 GenMuonGMTPair pairTmpCand(&(*gen), &(tkmuIt));
0386 float dr2Tmp = pairTmpCand.dR2();
0387 if (dr2Tmp < dr2Best) {
0388 dr2Best = dr2Tmp;
0389 pairBestCand2 = pairTmpCand;
0390 }
0391 }
0392 gmtTkMuonPairs_.emplace_back(pairBestCand2);
0393 }
0394 edm::LogInfo("L1TPhase2MuonOffline") << "L1TPhase2MuonOffline::matchMuonsToGen() gmtSAMuons: "
0395 << gmtSAMuonPairs_.size() << endl;
0396 edm::LogInfo("L1TPhase2MuonOffline") << "L1TPhase2MuonOffline::matchMuonsToGen() gmtTkMuons: "
0397 << gmtTkMuonPairs_.size() << endl;
0398 edm::LogInfo("L1TPhase2MuonOffline") << "L1TPhase2MuonOffline::matchMuonsToGen() END " << endl;
0399 }
0400
0401 std::vector<float> L1TPhase2MuonOffline::getHistBinsEff(EffType eff) {
0402 if (eff == kEffPt) {
0403 std::vector<float> effVsPtBins(effVsPtBins_.begin(), effVsPtBins_.end());
0404 return effVsPtBins;
0405 }
0406 if (eff == kEffPhi) {
0407 std::vector<float> effVsPhiBins(effVsPhiBins_.begin(), effVsPhiBins_.end());
0408 return effVsPhiBins;
0409 }
0410 if (eff == kEffEta) {
0411 std::vector<float> effVsEtaBins(effVsEtaBins_.begin(), effVsEtaBins_.end());
0412 return effVsEtaBins;
0413 }
0414 return {0., 1.};
0415 }
0416
0417 std::tuple<int, double, double> L1TPhase2MuonOffline::getHistBinsRes(ResType res) {
0418 if (res == kResPt)
0419 return {50, -2., 2.};
0420 if (res == kRes1OverPt)
0421 return {50, -2., 2.};
0422 if (res == kResQOverPt)
0423 return {50, -2., 2.};
0424 if (res == kResPhi)
0425 return {96, -0.2, 0.2};
0426 if (res == kResEta)
0427 return {100, -0.1, 0.1};
0428 if (res == kResCh)
0429 return {5, -2, 3};
0430 return {1, 0, 1};
0431 }
0432
0433
0434 DEFINE_FWK_MODULE(L1TPhase2MuonOffline);