Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:09:36

0001 /**
0002  * \file L1TPhase2MuonOffline.cc
0003  *
0004  * \author S. Folgueras 
0005  *
0006  */
0007 
0008 #include "DQMOffline/L1Trigger/interface/L1TPhase2MuonOffline.h"
0009 
0010 // To convert from HW to Physical coordinates
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 //__________RECO-GMT Muon Pair Helper Class____________________________
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();  // (1/gmtPt - 1/pt) / (1/pt)
0071   if (type == L1TPhase2MuonOffline::kResQOverPt)
0072     return (gmtCharge() * charge() * pt() - gmtPt()) /
0073            gmtPt();  // (gmtCharge/gmtPt - charge/pt) / (charge/pt) with gmtCharge/charge = gmtCharge*charge
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 //__________DQM_base_class_______________________________________________
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   //book histos
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   // COLLECT GEN MUONS
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   // Collect both muon collection:
0179   iEvent.getByToken(gmtMuonToken_, gmtSAMuon_);
0180   iEvent.getByToken(gmtTkMuonToken_, gmtTkMuon_);
0181 
0182   // Fill Control histograms
0183   edm::LogInfo("L1TPhase2MuonOffline") << "Fill Control histograms for GMT Muons" << endl;
0184   fillControlHistos();
0185 
0186   // Match each muon to a gen muon, if possible.
0187   if (genmus.empty())
0188     return;
0189   edm::LogInfo("L1TPhase2MuonOffline") << "L1TPhase2MuonOffline::analyze() calling matchMuonsToGen() " << endl;
0190   matchMuonsToGen(genmus);
0191 
0192   // Fill efficiency and resolution once, matching has been done...
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;  // there is not an assciated gmt muon
0296         if (muIt.gmtQual() < q * 4)
0297           continue;  //quality requirements
0298         const auto gmtPtCut = cut.first;
0299         if (var != kEffPt && muIt.gmtPt() < gmtPtCut)
0300           continue;  // pt requirement
0301         efficiencyNum_[kSAMuon][eta][q][var]->Fill(varToFill);
0302       }
0303     }
0304   }
0305 
0306   /// FOR TK MUONS
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;  // there is not an assciated gmt muon
0316         if (muIt.gmtQual() < q * 4)
0317           continue;  //quality requirements
0318         const auto gmtPtCut = cut.first;
0319         if (var != kEffPt && muIt.gmtPt() < gmtPtCut)
0320           continue;  // pt requirement
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 //define this as a plug-in
0434 DEFINE_FWK_MODULE(L1TPhase2MuonOffline);