File indexing completed on 2024-04-06 12:09:40
0001
0002 #include "DQMOffline/Muon/interface/EfficiencyAnalyzer.h"
0003
0004
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Framework/interface/ESHandle.h"
0010 #include "FWCore/Utilities/interface/Exception.h"
0011 #include "DataFormats/Math/interface/deltaR.h"
0012 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0013 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0014 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0015
0016 #include "TLorentzVector.h"
0017 #include "TFile.h"
0018 #include <vector>
0019 #include <cmath>
0020 #include <algorithm>
0021
0022
0023 #include <iostream>
0024 #include <fstream>
0025 #include <cmath>
0026 using namespace std;
0027 using namespace edm;
0028
0029 EfficiencyAnalyzer::EfficiencyAnalyzer(const edm::ParameterSet& pSet) {
0030 parameters = pSet;
0031
0032
0033 theMuonCollectionLabel_ = consumes<edm::View<reco::Muon> >(parameters.getParameter<edm::InputTag>("MuonCollection"));
0034 theTrackCollectionLabel_ = consumes<reco::TrackCollection>(parameters.getParameter<edm::InputTag>("TrackCollection"));
0035 theVertexLabel_ = consumes<reco::VertexCollection>(parameters.getParameter<edm::InputTag>("VertexLabel"));
0036 theBeamSpotLabel_ = mayConsume<reco::BeamSpot>(parameters.getParameter<edm::InputTag>("BeamSpotLabel"));
0037
0038
0039 doPVCheck_ = parameters.getParameter<bool>("doPrimaryVertexCheck");
0040
0041 ptBin_ = parameters.getParameter<int>("ptBin");
0042 ptMin_ = parameters.getParameter<double>("ptMin");
0043 ptMax_ = parameters.getParameter<double>("ptMax");
0044
0045 etaBin_ = parameters.getParameter<int>("etaBin");
0046 etaMin_ = parameters.getParameter<double>("etaMin");
0047 etaMax_ = parameters.getParameter<double>("etaMax");
0048
0049 phiBin_ = parameters.getParameter<int>("phiBin");
0050 phiMin_ = parameters.getParameter<double>("phiMin");
0051 phiMax_ = parameters.getParameter<double>("phiMax");
0052
0053 vtxBin_ = parameters.getParameter<int>("vtxBin");
0054 vtxMin_ = parameters.getParameter<double>("vtxMin");
0055 vtxMax_ = parameters.getParameter<double>("vtxMax");
0056
0057 ID_ = parameters.getParameter<string>("ID");
0058 theFolder = parameters.getParameter<string>("folder");
0059 }
0060
0061 EfficiencyAnalyzer::~EfficiencyAnalyzer() {}
0062
0063 void EfficiencyAnalyzer::bookHistograms(DQMStore::IBooker& ibooker,
0064 edm::Run const& ,
0065 edm::EventSetup const& ) {
0066 ibooker.cd();
0067 ibooker.setCurrentFolder(theFolder + ID_);
0068
0069 h_allProbes_pt = ibooker.book1D("allProbes_pt", "All Probes Pt", ptBin_, ptMin_, ptMax_);
0070 h_allProbes_inner_pt = ibooker.book1D("allProbes_inner_pt", "All Probes inner Pt", ptBin_, ptMin_, ptMax_);
0071 h_allProbes_inner_eta = ibooker.book1D("allProbes_inner_eta", "All Probes inner eta", etaBin_, etaMin_, etaMax_);
0072 h_allProbes_inner_phi = ibooker.book1D("allProbes_inner_phi", "All Probes inner phi", phiBin_, phiMin_, phiMax_);
0073 h_allProbes_EB_pt = ibooker.book1D("allProbes_EB_pt", "Barrel: all Probes Pt", ptBin_, ptMin_, ptMax_);
0074 h_allProbes_EE_pt = ibooker.book1D("allProbes_EE_pt", "Endcap: all Probes Pt", ptBin_, ptMin_, ptMax_);
0075 h_allProbes_eta = ibooker.book1D("allProbes_eta", "All Probes Eta", etaBin_, etaMin_, etaMax_);
0076 h_allProbes_hp_eta = ibooker.book1D("allProbes_hp_eta", "High Pt all Probes Eta", etaBin_, etaMin_, etaMax_);
0077 h_allProbes_phi = ibooker.book1D("allProbes_phi", "All Probes Phi", phiBin_, phiMin_, phiMax_);
0078
0079 h_allProbes_ID_pt = ibooker.book1D("allProbes_ID_pt", "All ID Probes Pt", ptBin_, ptMin_, ptMax_);
0080 h_allProbes_EB_ID_pt = ibooker.book1D("allProbes_EB_ID_pt", "Barrel: all ID Probes Pt", ptBin_, ptMin_, ptMax_);
0081 h_allProbes_EE_ID_pt = ibooker.book1D("allProbes_EE_ID_pt", "Endcap: all ID Probes Pt", ptBin_, ptMin_, ptMax_);
0082 h_allProbes_ID_nVtx = ibooker.book1D("allProbes_ID_nVtx", "All Probes (ID) nVtx", vtxBin_, vtxMin_, vtxMax_);
0083 h_allProbes_EB_ID_nVtx =
0084 ibooker.book1D("allProbes_EB_ID_nVtx", "Barrel: All Probes (ID) nVtx", vtxBin_, vtxMin_, vtxMax_);
0085 h_allProbes_EE_ID_nVtx =
0086 ibooker.book1D("allProbes_EE_ID_nVtx", "Endcap: All Probes (ID) nVtx", vtxBin_, vtxMin_, vtxMax_);
0087
0088 h_passProbes_ID_pt = ibooker.book1D("passProbes_ID_pt", "ID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
0089 h_passProbes_ID_inner_pt =
0090 ibooker.book1D("passProbes_ID_inner_pt", "ID Passing Probes inner Pt", ptBin_, ptMin_, ptMax_);
0091 h_passProbes_ID_inner_eta =
0092 ibooker.book1D("passProbes_ID_inner_eta", "ID Passing Probes inner eta", etaBin_, etaMin_, etaMax_);
0093 h_passProbes_ID_inner_phi =
0094 ibooker.book1D("passProbes_ID_inner_phi", "ID Passing Probes inner phi", phiBin_, phiMin_, phiMax_);
0095 h_passProbes_ID_EB_pt = ibooker.book1D("passProbes_ID_EB_pt", "Barrel: ID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
0096 h_passProbes_ID_EE_pt = ibooker.book1D("passProbes_ID_EE_pt", "Endcap: ID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
0097 h_passProbes_ID_eta = ibooker.book1D("passProbes_ID_eta", "ID Passing Probes #eta", etaBin_, etaMin_, etaMax_);
0098 h_passProbes_ID_hp_eta =
0099 ibooker.book1D("passProbes_ID_hp_eta", "High Pt ID Passing Probes #eta", etaBin_, etaMin_, etaMax_);
0100 h_passProbes_ID_phi = ibooker.book1D("passProbes_ID_phi", "ID Passing Probes #phi", phiBin_, phiMin_, phiMax_);
0101
0102 h_passProbes_detIsoID_pt =
0103 ibooker.book1D("passProbes_detIsoID_pt", "detIsoID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
0104 h_passProbes_EB_detIsoID_pt =
0105 ibooker.book1D("passProbes_EB_detIsoID_pt", "Barrel: detIsoID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
0106 h_passProbes_EE_detIsoID_pt =
0107 ibooker.book1D("passProbes_EE_detIsoID_pt", "Endcap: detIsoID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
0108
0109 h_passProbes_pfIsoID_pt =
0110 ibooker.book1D("passProbes_pfIsoID_pt", "pfIsoID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
0111 h_passProbes_EB_pfIsoID_pt =
0112 ibooker.book1D("passProbes_EB_pfIsoID_pt", "Barrel: pfIsoID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
0113 h_passProbes_EE_pfIsoID_pt =
0114 ibooker.book1D("passProbes_EE_pfIsoID_pt", "Endcap: pfIsoID Passing Probes Pt", ptBin_, ptMin_, ptMax_);
0115
0116 h_passProbes_detIsoID_nVtx =
0117 ibooker.book1D("passProbes_detIsoID_nVtx", "detIsoID Passing Probes nVtx (R03)", vtxBin_, vtxMin_, vtxMax_);
0118 h_passProbes_pfIsoID_nVtx =
0119 ibooker.book1D("passProbes_pfIsoID_nVtx", "pfIsoID Passing Probes nVtx (R04)", vtxBin_, vtxMin_, vtxMax_);
0120 h_passProbes_EB_detIsoID_nVtx = ibooker.book1D(
0121 "passProbes_EB_detIsoID_nVtx", "Barrel: detIsoID Passing Probes nVtx (R03)", vtxBin_, vtxMin_, vtxMax_);
0122 h_passProbes_EE_detIsoID_nVtx = ibooker.book1D(
0123 "passProbes_EE_detIsoID_nVtx", "Endcap: detIsoID Passing Probes nVtx (R03)", vtxBin_, vtxMin_, vtxMax_);
0124 h_passProbes_EB_pfIsoID_nVtx = ibooker.book1D(
0125 "passProbes_EB_pfIsoID_nVtx", "Barrel: pfIsoID Passing Probes nVtx (R04)", vtxBin_, vtxMin_, vtxMax_);
0126 h_passProbes_EE_pfIsoID_nVtx = ibooker.book1D(
0127 "passProbes_EE_pfIsoID_nVtx", "Endcap: pfIsoID Passing Probes nVtx (R04)", vtxBin_, vtxMin_, vtxMax_);
0128
0129
0130
0131 h_passProbes_pfIsodBID_pt = ibooker.book1D(
0132 "passProbes_pfIsodBID_pt", "pfIsoID Passing Probes Pt (deltaB PU correction)", ptBin_, ptMin_, ptMax_);
0133 h_passProbes_EB_pfIsodBID_pt = ibooker.book1D(
0134 "passProbes_EB_pfIsodBID_pt", "Barrel: pfIsoID Passing Probes Pt (deltaB PU correction)", ptBin_, ptMin_, ptMax_);
0135 h_passProbes_EE_pfIsodBID_pt = ibooker.book1D(
0136 "passProbes_EE_pfIsodBID_pt", "Endcap: pfIsoID Passing Probes Pt (deltaB PU correction)", ptBin_, ptMin_, ptMax_);
0137 h_passProbes_pfIsodBID_nVtx = ibooker.book1D("passProbes_pfIsodBID_nVtx",
0138 "pfIsoID Passing Probes nVtx (R04) (deltaB PU correction)",
0139 vtxBin_,
0140 vtxMin_,
0141 vtxMax_);
0142 h_passProbes_EB_pfIsodBID_nVtx = ibooker.book1D("passProbes_EB_pfIsodBID_nVtx",
0143 "Barrel: pfIsoID Passing Probes nVtx (R04) (deltaB PU correction)",
0144 vtxBin_,
0145 vtxMin_,
0146 vtxMax_);
0147 h_passProbes_EE_pfIsodBID_nVtx = ibooker.book1D("passProbes_EE_pfIsodBID_nVtx",
0148 "Endcap: pfIsoID Passing Probes nVtx (R04) (deltaB PU correction)",
0149 vtxBin_,
0150 vtxMin_,
0151 vtxMax_);
0152
0153 #ifdef DEBUG
0154 cout << "[EfficiencyAnalyzer] Parameters initialization DONE" << endl;
0155 #endif
0156 }
0157
0158 void EfficiencyAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0159 LogTrace(metname) << "[EfficiencyAnalyzer] Analyze the mu in different eta regions";
0160
0161
0162
0163 edm::Handle<edm::View<reco::Muon> > muons;
0164 iEvent.getByToken(theMuonCollectionLabel_, muons);
0165
0166
0167 edm::Handle<reco::TrackCollection> tracks;
0168 iEvent.getByToken(theTrackCollectionLabel_, tracks);
0169
0170
0171 edm::Handle<reco::VertexCollection> vertex;
0172 iEvent.getByToken(theVertexLabel_, vertex);
0173
0174
0175
0176 _numPV = 0;
0177 bool bPrimaryVertex = true;
0178 if (doPVCheck_) {
0179 bPrimaryVertex = false;
0180
0181 if (!vertex.isValid()) {
0182 LogTrace(metname) << "[EfficiencyAnalyzer] Could not find vertex collection" << std::endl;
0183 bPrimaryVertex = false;
0184 }
0185
0186 if (vertex.isValid()) {
0187 const reco::VertexCollection& vertexCollection = *(vertex.product());
0188 int vertex_number = vertexCollection.size();
0189
0190 reco::VertexCollection::const_iterator v = vertexCollection.begin();
0191 for (; v != vertexCollection.end(); ++v) {
0192 double vertex_chi2 = v->normalizedChi2();
0193 double vertex_ndof = v->ndof();
0194 bool fakeVtx = v->isFake();
0195 double vertex_Z = v->z();
0196
0197 if (!fakeVtx && vertex_number >= 1 && vertex_ndof > 4 && vertex_chi2 < 999 && fabs(vertex_Z) < 24.) {
0198 bPrimaryVertex = true;
0199 ++_numPV;
0200 }
0201 }
0202 }
0203 }
0204
0205
0206
0207 reco::Vertex::Point posVtx;
0208 reco::Vertex::Error errVtx;
0209 unsigned int theIndexOfThePrimaryVertex = 999.;
0210 if (vertex.isValid()) {
0211 for (unsigned int ind = 0; ind < vertex->size(); ++ind) {
0212 if ((*vertex)[ind].isValid() && !((*vertex)[ind].isFake())) {
0213 theIndexOfThePrimaryVertex = ind;
0214 break;
0215 }
0216 }
0217 }
0218
0219 if (theIndexOfThePrimaryVertex < 100) {
0220 posVtx = ((*vertex)[theIndexOfThePrimaryVertex]).position();
0221 errVtx = ((*vertex)[theIndexOfThePrimaryVertex]).error();
0222 } else {
0223 LogInfo("RecoMuonValidator") << "reco::PrimaryVertex not found, use BeamSpot position instead\n";
0224
0225 edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0226 iEvent.getByToken(theBeamSpotLabel_, recoBeamSpotHandle);
0227 reco::BeamSpot bs = *recoBeamSpotHandle;
0228
0229 posVtx = bs.position();
0230 errVtx(0, 0) = bs.BeamWidthX();
0231 errVtx(1, 1) = bs.BeamWidthY();
0232 errVtx(2, 2) = bs.sigmaZ();
0233 }
0234
0235 const reco::Vertex thePrimaryVertex(posVtx, errVtx);
0236
0237
0238 if (!muons.isValid())
0239 return;
0240
0241
0242 TLorentzVector Mu1, Mu2;
0243
0244 bool isMB = false;
0245 bool isME = false;
0246
0247 for (edm::View<reco::Muon>::const_iterator muon1 = muons->begin(); muon1 != muons->end(); ++muon1) {
0248 LogTrace(metname) << "[EfficiencyAnalyzer] loop over first muons" << endl;
0249
0250
0251 reco::MuonIsolation Iso_muon = muon1->isolationR03();
0252 float combIso = (Iso_muon.emEt + Iso_muon.hadEt + Iso_muon.sumPt);
0253
0254
0255 if (!muon1->isGlobalMuon())
0256 continue;
0257
0258
0259 reco::TrackRef recoCombinedGlbTrack1 = muon1->combinedMuon();
0260 float muPt1 = recoCombinedGlbTrack1->pt();
0261 Mu1.SetPxPyPzE(recoCombinedGlbTrack1->px(),
0262 recoCombinedGlbTrack1->py(),
0263 recoCombinedGlbTrack1->pz(),
0264 recoCombinedGlbTrack1->p());
0265
0266
0267
0268 if (ID_ == "Loose" && !muon::isLooseMuon(*muon1))
0269 continue;
0270 if (ID_ == "Medium" && !muon::isMediumMuon(*muon1))
0271 continue;
0272 if (ID_ == "Tight" && !muon::isTightMuon(*muon1, thePrimaryVertex))
0273 continue;
0274
0275
0276 if (muPt1 <= 15)
0277 continue;
0278 if (combIso / muPt1 > 0.1)
0279 continue;
0280
0281 for (edm::View<reco::Muon>::const_iterator muon2 = muons->begin(); muon2 != muons->end(); ++muon2) {
0282 LogTrace(metname) << "[EfficiencyAnalyzer] loop over second muon" << endl;
0283 if (muon2 == muon1)
0284 continue;
0285
0286 if (muon2->eta() < 1.479)
0287 isMB = true;
0288 if (muon2->eta() >= 1.479)
0289 isME = true;
0290
0291
0292 Mu2.SetPxPyPzE(muon2->px(), muon2->py(), muon2->pz(), muon2->p());
0293
0294 float Minv = (Mu1 + Mu2).M();
0295 if (!muon2->isTrackerMuon())
0296 continue;
0297 if (muon2->pt() < 5)
0298 continue;
0299 if ((muon1->charge()) * (muon2->charge()) > 0)
0300 continue;
0301 if (Minv < 70 || Minv > 110)
0302 continue;
0303
0304 h_allProbes_pt->Fill(muon2->pt());
0305 h_allProbes_eta->Fill(muon2->eta());
0306 h_allProbes_phi->Fill(muon2->phi());
0307 if (muon2->innerTrack()->extra().isAvailable()) {
0308 h_allProbes_inner_pt->Fill(muon2->innerTrack()->innerMomentum().Rho());
0309 h_allProbes_inner_eta->Fill(muon2->innerTrack()->innerPosition().Eta());
0310 h_allProbes_inner_phi->Fill(muon2->innerTrack()->innerPosition().Phi());
0311 }
0312 if (isMB)
0313 h_allProbes_EB_pt->Fill(muon2->pt());
0314 if (isME)
0315 h_allProbes_EE_pt->Fill(muon2->pt());
0316 if (muon2->pt() > 20)
0317 h_allProbes_hp_eta->Fill(muon2->eta());
0318
0319
0320 if (ID_ == "Loose" && !muon::isLooseMuon(*muon2))
0321 continue;
0322 if (ID_ == "Medium" && !muon::isMediumMuon(*muon2))
0323 continue;
0324 if (ID_ == "Tight" && !muon::isTightMuon(*muon2, thePrimaryVertex))
0325 continue;
0326
0327 h_passProbes_ID_pt->Fill(muon2->pt());
0328 h_passProbes_ID_eta->Fill(muon2->eta());
0329 h_passProbes_ID_phi->Fill(muon2->phi());
0330 if (muon2->innerTrack()->extra().isAvailable()) {
0331 h_passProbes_ID_inner_pt->Fill(muon2->innerTrack()->innerMomentum().Rho());
0332 h_passProbes_ID_inner_eta->Fill(muon2->innerTrack()->innerPosition().Eta());
0333 h_passProbes_ID_inner_phi->Fill(muon2->innerTrack()->innerPosition().Phi());
0334 }
0335
0336 if (isMB)
0337 h_passProbes_ID_EB_pt->Fill(muon2->pt());
0338 if (isME)
0339 h_passProbes_ID_EE_pt->Fill(muon2->pt());
0340 if (muon2->pt() > 20)
0341 h_passProbes_ID_hp_eta->Fill(muon2->eta());
0342
0343 h_allProbes_ID_pt->Fill(muon2->pt());
0344 if (isMB)
0345 h_allProbes_EB_ID_pt->Fill(muon2->pt());
0346 if (isME)
0347 h_allProbes_EE_ID_pt->Fill(muon2->pt());
0348
0349
0350 if (bPrimaryVertex)
0351 h_allProbes_ID_nVtx->Fill(_numPV);
0352 if (bPrimaryVertex && isMB)
0353 h_allProbes_EB_ID_nVtx->Fill(_numPV);
0354 if (bPrimaryVertex && isME)
0355 h_allProbes_EE_ID_nVtx->Fill(_numPV);
0356
0357
0358 float tkIso = muon2->isolationR03().sumPt;
0359 float emIso = muon2->isolationR03().emEt;
0360 float hadIso = muon2->isolationR03().hadEt + muon2->isolationR03().hoEt;
0361 float relDetIso = (tkIso + emIso + hadIso) / (muon2->pt());
0362
0363 if (relDetIso < 0.05) {
0364 h_passProbes_detIsoID_pt->Fill(muon2->pt());
0365 if (isMB)
0366 h_passProbes_EB_detIsoID_pt->Fill(muon2->pt());
0367 if (isME)
0368 h_passProbes_EE_detIsoID_pt->Fill(muon2->pt());
0369
0370 if (bPrimaryVertex)
0371 h_passProbes_detIsoID_nVtx->Fill(_numPV);
0372 if (bPrimaryVertex && isMB)
0373 h_passProbes_EB_detIsoID_nVtx->Fill(_numPV);
0374 if (bPrimaryVertex && isME)
0375 h_passProbes_EE_detIsoID_nVtx->Fill(_numPV);
0376 }
0377
0378
0379 float chargedIso = muon2->pfIsolationR04().sumChargedHadronPt;
0380 float neutralIso = muon2->pfIsolationR04().sumNeutralHadronEt;
0381 float photonIso = muon2->pfIsolationR04().sumPhotonEt;
0382 float relPFIso = (chargedIso + neutralIso + photonIso) / (muon2->pt());
0383
0384 float pu = muon2->pfIsolationR04().sumPUPt;
0385 float neutralphotonPUCorrected = std::max(0.0, (neutralIso + photonIso - 0.5 * pu));
0386 float relPFIsoPUCorrected = (chargedIso + neutralphotonPUCorrected) / (muon2->pt());
0387
0388 if (relPFIso < 0.12) {
0389 h_passProbes_pfIsoID_pt->Fill(muon2->pt());
0390 if (isMB)
0391 h_passProbes_EB_pfIsoID_pt->Fill(muon2->pt());
0392 if (isME)
0393 h_passProbes_EE_pfIsoID_pt->Fill(muon2->pt());
0394
0395 if (bPrimaryVertex)
0396 h_passProbes_pfIsoID_nVtx->Fill(_numPV);
0397 if (bPrimaryVertex && isMB)
0398 h_passProbes_EB_pfIsoID_nVtx->Fill(_numPV);
0399 if (bPrimaryVertex && isME)
0400 h_passProbes_EE_pfIsoID_nVtx->Fill(_numPV);
0401 }
0402
0403
0404 if (relPFIsoPUCorrected < 0.12) {
0405 h_passProbes_pfIsodBID_pt->Fill(muon2->pt());
0406 if (isMB)
0407 h_passProbes_EB_pfIsodBID_pt->Fill(muon2->pt());
0408 if (isME)
0409 h_passProbes_EE_pfIsodBID_pt->Fill(muon2->pt());
0410
0411 if (bPrimaryVertex)
0412 h_passProbes_pfIsodBID_nVtx->Fill(_numPV);
0413 if (bPrimaryVertex && isMB)
0414 h_passProbes_EB_pfIsodBID_nVtx->Fill(_numPV);
0415 if (bPrimaryVertex && isME)
0416 h_passProbes_EE_pfIsodBID_nVtx->Fill(_numPV);
0417 }
0418 }
0419 }
0420 }