File indexing completed on 2024-09-11 04:32:50
0001
0002 #include "DQMOffline/Muon/interface/DiMuonHistograms.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/MuonReco/interface/MuonSelectors.h"
0012 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0013 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0014 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0015 #include "DataFormats/TrackReco/interface/Track.h"
0016
0017 #include "DQMServices/Core/interface/DQMStore.h"
0018
0019 #include "TLorentzVector.h"
0020 #include "TFile.h"
0021 #include <vector>
0022 #include <cmath>
0023
0024
0025 #include <iostream>
0026 #include <fstream>
0027 #include <cmath>
0028 using namespace std;
0029 using namespace edm;
0030
0031 DiMuonHistograms::DiMuonHistograms(const edm::ParameterSet& pSet) {
0032
0033 parameters = pSet;
0034
0035
0036 nTightTight = 0;
0037 nMediumMedium = 0;
0038 nLooseLoose = 0;
0039 nGlbGlb = 0;
0040
0041
0042 theMuonCollectionLabel_ = consumes<edm::View<reco::Muon> >(parameters.getParameter<edm::InputTag>("MuonCollection"));
0043 theVertexLabel_ = consumes<reco::VertexCollection>(parameters.getParameter<edm::InputTag>("VertexLabel"));
0044
0045 theBeamSpotLabel_ = mayConsume<reco::BeamSpot>(parameters.getParameter<edm::InputTag>("BeamSpotLabel"));
0046
0047 etaBin = parameters.getParameter<int>("etaBin");
0048 etaBBin = parameters.getParameter<int>("etaBBin");
0049 etaEBin = parameters.getParameter<int>("etaEBin");
0050
0051 etaBMin = parameters.getParameter<double>("etaBMin");
0052 etaBMax = parameters.getParameter<double>("etaBMax");
0053 etaECMin = parameters.getParameter<double>("etaECMin");
0054 etaECMax = parameters.getParameter<double>("etaECMax");
0055
0056 lowMassMin = parameters.getParameter<double>("lowMassMin");
0057 lowMassMax = parameters.getParameter<double>("lowMassMax");
0058 highMassMin = parameters.getParameter<double>("highMassMin");
0059 highMassMax = parameters.getParameter<double>("highMassMax");
0060
0061 theFolder = parameters.getParameter<string>("folder");
0062 }
0063
0064 DiMuonHistograms::~DiMuonHistograms() {}
0065
0066 void DiMuonHistograms::bookHistograms(DQMStore::IBooker& ibooker,
0067 edm::Run const& ,
0068 edm::EventSetup const& ) {
0069 ibooker.cd();
0070 ibooker.setCurrentFolder(theFolder);
0071
0072 int nBin[3] = {etaBin, etaBBin, etaEBin};
0073 EtaName[0] = "";
0074 EtaName[1] = "_Barrel";
0075 EtaName[2] = "_EndCap";
0076 test = ibooker.book1D("test", "InvMass_{Tight,Tight}", 100, 0., 200.);
0077 for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
0078 GlbGlbMuon_LM.push_back(ibooker.book1D("GlbGlbMuon_LM" + EtaName[iEtaRegion],
0079 "InvMass_{GLB,GLB}" + EtaName[iEtaRegion],
0080 nBin[iEtaRegion],
0081 lowMassMin,
0082 lowMassMax));
0083 TrkTrkMuon_LM.push_back(ibooker.book1D("TrkTrkMuon_LM" + EtaName[iEtaRegion],
0084 "InvMass_{TRK,TRK}" + EtaName[iEtaRegion],
0085 nBin[iEtaRegion],
0086 lowMassMin,
0087 lowMassMax));
0088 StaTrkMuon_LM.push_back(ibooker.book1D("StaTrkMuon_LM" + EtaName[iEtaRegion],
0089 "InvMass_{STA,TRK}" + EtaName[iEtaRegion],
0090 nBin[iEtaRegion],
0091 lowMassMin,
0092 lowMassMax));
0093
0094 GlbGlbMuon_HM.push_back(ibooker.book1D("GlbGlbMuon_HM" + EtaName[iEtaRegion],
0095 "InvMass_{GLB,GLB}" + EtaName[iEtaRegion],
0096 nBin[iEtaRegion],
0097 highMassMin,
0098 highMassMax));
0099 TrkTrkMuon_HM.push_back(ibooker.book1D("TrkTrkMuon_HM" + EtaName[iEtaRegion],
0100 "InvMass_{TRK,TRK}" + EtaName[iEtaRegion],
0101 nBin[iEtaRegion],
0102 highMassMin,
0103 highMassMax));
0104 StaTrkMuon_HM.push_back(ibooker.book1D("StaTrkMuon_HM" + EtaName[iEtaRegion],
0105 "InvMass_{STA,TRK}" + EtaName[iEtaRegion],
0106 nBin[iEtaRegion],
0107 highMassMin,
0108 highMassMax));
0109
0110
0111 TightTightMuon.push_back(ibooker.book1D("TightTightMuon" + EtaName[iEtaRegion],
0112 "InvMass_{Tight,Tight}" + EtaName[iEtaRegion],
0113 nBin[iEtaRegion],
0114 highMassMin,
0115 highMassMax));
0116 MediumMediumMuon.push_back(ibooker.book1D("MediumMediumMuon" + EtaName[iEtaRegion],
0117 "InvMass_{Medium,Medium}" + EtaName[iEtaRegion],
0118 nBin[iEtaRegion],
0119 highMassMin,
0120 highMassMax));
0121 LooseLooseMuon.push_back(ibooker.book1D("LooseLooseMuon" + EtaName[iEtaRegion],
0122 "InvMass_{Loose,Loose}" + EtaName[iEtaRegion],
0123 nBin[iEtaRegion],
0124 highMassMin,
0125 highMassMax));
0126
0127 TightTightMuonBadFrac.push_back(ibooker.book1D(
0128 "TightTightMuonBadFrac" + EtaName[iEtaRegion], "BadFrac_{Tight,Tight}" + EtaName[iEtaRegion], 10, 0, 0.4));
0129 MediumMediumMuonBadFrac.push_back(ibooker.book1D(
0130 "MediumMediumMuonBadFrac" + EtaName[iEtaRegion], "BadFrac_{Medium,Medium}" + EtaName[iEtaRegion], 10, 0, 0.4));
0131 LooseLooseMuonBadFrac.push_back(ibooker.book1D(
0132 "LooseLooseMuonBadFrac" + EtaName[iEtaRegion], "BadFrac_{Loose,Loose}" + EtaName[iEtaRegion], 10, 0, 0.4));
0133
0134
0135 SoftSoftMuon.push_back(ibooker.book1D(
0136 "SoftSoftMuon" + EtaName[iEtaRegion], "InvMass_{Soft,Soft}" + EtaName[iEtaRegion], nBin[iEtaRegion], 0.0, 55.0));
0137 SoftSoftMuonBadFrac.push_back(ibooker.book1D(
0138 "SoftSoftMuonBadFrac" + EtaName[iEtaRegion], "BadFrac_{Soft,Soft}" + EtaName[iEtaRegion], 10, 0, 0.4));
0139 }
0140 }
0141
0142 void DiMuonHistograms::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0143 LogTrace(metname) << "[DiMuonHistograms] Analyze the mu in different eta regions";
0144 edm::Handle<edm::View<reco::Muon> > muons;
0145 iEvent.getByToken(theMuonCollectionLabel_, muons);
0146
0147
0148
0149 reco::Vertex::Point posVtx;
0150 reco::Vertex::Error errVtx;
0151 unsigned int theIndexOfThePrimaryVertex = 999.;
0152
0153 edm::Handle<reco::VertexCollection> vertex;
0154 iEvent.getByToken(theVertexLabel_, vertex);
0155 if (vertex.isValid()) {
0156 for (unsigned int ind = 0; ind < vertex->size(); ++ind) {
0157 if ((*vertex)[ind].isValid() && !((*vertex)[ind].isFake())) {
0158 theIndexOfThePrimaryVertex = ind;
0159 break;
0160 }
0161 }
0162 }
0163
0164 if (theIndexOfThePrimaryVertex < 100) {
0165 posVtx = ((*vertex)[theIndexOfThePrimaryVertex]).position();
0166 errVtx = ((*vertex)[theIndexOfThePrimaryVertex]).error();
0167 } else {
0168 LogInfo("RecoMuonValidator") << "reco::PrimaryVertex not found, use BeamSpot position instead\n";
0169
0170 edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0171 iEvent.getByToken(theBeamSpotLabel_, recoBeamSpotHandle);
0172 reco::BeamSpot bs = *recoBeamSpotHandle;
0173
0174 posVtx = bs.position();
0175 errVtx(0, 0) = bs.BeamWidthX();
0176 errVtx(1, 1) = bs.BeamWidthY();
0177 errVtx(2, 2) = bs.sigmaZ();
0178 }
0179
0180 const reco::Vertex vtx(posVtx, errVtx);
0181
0182 if (!muons.isValid())
0183 return;
0184
0185
0186 TLorentzVector Mu1, Mu2;
0187 float charge = 99.;
0188 float InvMass = -99.;
0189
0190
0191 double EtaCutMin[] = {0, etaBMin, etaECMin};
0192 double EtaCutMax[] = {2.4, etaBMax, etaECMax};
0193
0194 for (edm::View<reco::Muon>::const_iterator muon1 = muons->begin(); muon1 != muons->end(); ++muon1) {
0195 LogTrace(metname) << "[DiMuonHistograms] loop over 1st muon" << endl;
0196
0197
0198 for (edm::View<reco::Muon>::const_iterator muon2 = muon1; muon2 != muons->end(); ++muon2) {
0199 LogTrace(metname) << "[DiMuonHistograms] loop over 2nd muon" << endl;
0200 if (muon1 == muon2)
0201 continue;
0202
0203
0204 if (muon1->isGlobalMuon() && muon2->isGlobalMuon()) {
0205 LogTrace(metname) << "[DiMuonHistograms] Glb-Glb pair" << endl;
0206 reco::TrackRef recoCombinedGlbTrack1 = muon1->combinedMuon();
0207 reco::TrackRef recoCombinedGlbTrack2 = muon2->combinedMuon();
0208 Mu1.SetPxPyPzE(recoCombinedGlbTrack1->px(),
0209 recoCombinedGlbTrack1->py(),
0210 recoCombinedGlbTrack1->pz(),
0211 recoCombinedGlbTrack1->p());
0212 Mu2.SetPxPyPzE(recoCombinedGlbTrack2->px(),
0213 recoCombinedGlbTrack2->py(),
0214 recoCombinedGlbTrack2->pz(),
0215 recoCombinedGlbTrack2->p());
0216
0217 charge = recoCombinedGlbTrack1->charge() * recoCombinedGlbTrack2->charge();
0218 if (charge < 0) {
0219 InvMass = (Mu1 + Mu2).M();
0220 for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
0221 if (fabs(recoCombinedGlbTrack1->eta()) > EtaCutMin[iEtaRegion] &&
0222 fabs(recoCombinedGlbTrack1->eta()) < EtaCutMax[iEtaRegion] &&
0223 fabs(recoCombinedGlbTrack2->eta()) > EtaCutMin[iEtaRegion] &&
0224 fabs(recoCombinedGlbTrack2->eta()) < EtaCutMax[iEtaRegion]) {
0225 if (InvMass < lowMassMax)
0226 GlbGlbMuon_LM[iEtaRegion]->Fill(InvMass);
0227 if (InvMass > highMassMin)
0228 GlbGlbMuon_HM[iEtaRegion]->Fill(InvMass);
0229 }
0230 }
0231 }
0232
0233 if (muon::isTightMuon(*muon1, vtx) && muon::isTightMuon(*muon2, vtx)) {
0234 test->Fill(InvMass);
0235 LogTrace(metname) << "[DiMuonHistograms] Tight-Tight pair" << endl;
0236 if (charge < 0) {
0237 for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
0238 if (fabs(recoCombinedGlbTrack1->eta()) > EtaCutMin[iEtaRegion] &&
0239 fabs(recoCombinedGlbTrack1->eta()) < EtaCutMax[iEtaRegion] &&
0240 fabs(recoCombinedGlbTrack2->eta()) > EtaCutMin[iEtaRegion] &&
0241 fabs(recoCombinedGlbTrack2->eta()) < EtaCutMax[iEtaRegion]) {
0242 if (InvMass > 55. && InvMass < 125.) {
0243 TightTightMuon[iEtaRegion]->Fill(InvMass);
0244 TightTightMuonBadFrac[iEtaRegion]->Fill(muon1->innerTrack()->lost() / muon1->innerTrack()->found());
0245 }
0246 }
0247 }
0248 }
0249 }
0250
0251 if (muon::isMediumMuon(*muon1) && muon::isMediumMuon(*muon2)) {
0252 test->Fill(InvMass);
0253 LogTrace(metname) << "[DiMuonHistograms] Medium-Medium pair" << endl;
0254 if (charge < 0) {
0255 for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
0256 if (fabs(recoCombinedGlbTrack1->eta()) > EtaCutMin[iEtaRegion] &&
0257 fabs(recoCombinedGlbTrack1->eta()) < EtaCutMax[iEtaRegion] &&
0258 fabs(recoCombinedGlbTrack2->eta()) > EtaCutMin[iEtaRegion] &&
0259 fabs(recoCombinedGlbTrack2->eta()) < EtaCutMax[iEtaRegion]) {
0260 if (InvMass > 55. && InvMass < 125.) {
0261 MediumMediumMuon[iEtaRegion]->Fill(InvMass);
0262 MediumMediumMuonBadFrac[iEtaRegion]->Fill(muon1->innerTrack()->lost() / muon1->innerTrack()->found());
0263 }
0264 }
0265 }
0266 }
0267 }
0268
0269 if (muon::isLooseMuon(*muon1) && muon::isLooseMuon(*muon2)) {
0270 test->Fill(InvMass);
0271 LogTrace(metname) << "[DiMuonHistograms] Loose-Loose pair" << endl;
0272 if (charge < 0) {
0273 for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
0274 if (fabs(recoCombinedGlbTrack1->eta()) > EtaCutMin[iEtaRegion] &&
0275 fabs(recoCombinedGlbTrack1->eta()) < EtaCutMax[iEtaRegion] &&
0276 fabs(recoCombinedGlbTrack2->eta()) > EtaCutMin[iEtaRegion] &&
0277 fabs(recoCombinedGlbTrack2->eta()) < EtaCutMax[iEtaRegion]) {
0278 if (InvMass > 55. && InvMass < 125.) {
0279 LooseLooseMuon[iEtaRegion]->Fill(InvMass);
0280 LooseLooseMuonBadFrac[iEtaRegion]->Fill(muon1->innerTrack()->lost() / muon1->innerTrack()->found());
0281 }
0282 }
0283 }
0284 }
0285 }
0286 }
0287
0288
0289 if (muon2->isStandAloneMuon() && muon1->isTrackerMuon()) {
0290 LogTrace(metname) << "[DiMuonHistograms] STA-Trk pair" << endl;
0291 reco::TrackRef recoStaTrack = muon2->standAloneMuon();
0292 reco::TrackRef recoTrack = muon1->track();
0293 Mu2.SetPxPyPzE(recoStaTrack->px(), recoStaTrack->py(), recoStaTrack->pz(), recoStaTrack->p());
0294 Mu1.SetPxPyPzE(recoTrack->px(), recoTrack->py(), recoTrack->pz(), recoTrack->p());
0295
0296 charge = recoStaTrack->charge() * recoTrack->charge();
0297 if (charge < 0) {
0298 InvMass = (Mu1 + Mu2).M();
0299 for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
0300 if (fabs(recoStaTrack->eta()) > EtaCutMin[iEtaRegion] &&
0301 fabs(recoStaTrack->eta()) < EtaCutMax[iEtaRegion] && fabs(recoTrack->eta()) > EtaCutMin[iEtaRegion] &&
0302 fabs(recoTrack->eta()) < EtaCutMax[iEtaRegion]) {
0303 if (InvMass < lowMassMax)
0304 StaTrkMuon_LM[iEtaRegion]->Fill(InvMass);
0305 if (InvMass > highMassMin)
0306 StaTrkMuon_HM[iEtaRegion]->Fill(InvMass);
0307 }
0308 }
0309 }
0310 }
0311 if (muon1->isStandAloneMuon() && muon2->isTrackerMuon()) {
0312 LogTrace(metname) << "[DiMuonHistograms] STA-Trk pair" << endl;
0313 reco::TrackRef recoStaTrack = muon1->standAloneMuon();
0314 reco::TrackRef recoTrack = muon2->track();
0315 Mu1.SetPxPyPzE(recoStaTrack->px(), recoStaTrack->py(), recoStaTrack->pz(), recoStaTrack->p());
0316 Mu2.SetPxPyPzE(recoTrack->px(), recoTrack->py(), recoTrack->pz(), recoTrack->p());
0317
0318 charge = recoStaTrack->charge() * recoTrack->charge();
0319 if (charge < 0) {
0320 InvMass = (Mu1 + Mu2).M();
0321 for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
0322 if (fabs(recoStaTrack->eta()) > EtaCutMin[iEtaRegion] &&
0323 fabs(recoStaTrack->eta()) < EtaCutMax[iEtaRegion] && fabs(recoTrack->eta()) > EtaCutMin[iEtaRegion] &&
0324 fabs(recoTrack->eta()) < EtaCutMax[iEtaRegion]) {
0325 if (InvMass < lowMassMax)
0326 StaTrkMuon_LM[iEtaRegion]->Fill(InvMass);
0327 if (InvMass > highMassMin)
0328 StaTrkMuon_HM[iEtaRegion]->Fill(InvMass);
0329 }
0330 }
0331 }
0332 }
0333
0334
0335 if (muon1->isTrackerMuon() && muon2->isTrackerMuon()) {
0336 LogTrace(metname) << "[DiMuonHistograms] Trk-Trk dimuon pair" << endl;
0337 reco::TrackRef recoTrack2 = muon2->track();
0338 reco::TrackRef recoTrack1 = muon1->track();
0339 Mu2.SetPxPyPzE(recoTrack2->px(), recoTrack2->py(), recoTrack2->pz(), recoTrack2->p());
0340 Mu1.SetPxPyPzE(recoTrack1->px(), recoTrack1->py(), recoTrack1->pz(), recoTrack1->p());
0341
0342 charge = recoTrack1->charge() * recoTrack2->charge();
0343 if (charge < 0) {
0344 InvMass = (Mu1 + Mu2).M();
0345 for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
0346 if (fabs(recoTrack1->eta()) > EtaCutMin[iEtaRegion] && fabs(recoTrack1->eta()) < EtaCutMax[iEtaRegion] &&
0347 fabs(recoTrack2->eta()) > EtaCutMin[iEtaRegion] && fabs(recoTrack2->eta()) < EtaCutMax[iEtaRegion]) {
0348 if (InvMass < lowMassMax)
0349 TrkTrkMuon_LM[iEtaRegion]->Fill(InvMass);
0350 if (InvMass > highMassMin)
0351 TrkTrkMuon_HM[iEtaRegion]->Fill(InvMass);
0352 }
0353 }
0354 }
0355
0356 LogTrace(metname) << "[DiMuonHistograms] Soft-Soft pair" << endl;
0357
0358 if (muon::isSoftMuon(*muon1, vtx) && muon::isSoftMuon(*muon2, vtx)) {
0359 if (charge < 0) {
0360 InvMass = (Mu1 + Mu2).M();
0361 for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
0362 if (fabs(recoTrack1->eta()) > EtaCutMin[iEtaRegion] && fabs(recoTrack1->eta()) < EtaCutMax[iEtaRegion] &&
0363 fabs(recoTrack2->eta()) > EtaCutMin[iEtaRegion] && fabs(recoTrack2->eta()) < EtaCutMax[iEtaRegion]) {
0364 SoftSoftMuon[iEtaRegion]->Fill(InvMass);
0365 SoftSoftMuonBadFrac[iEtaRegion]->Fill(muon1->innerTrack()->lost() / muon1->innerTrack()->found());
0366 }
0367 }
0368 }
0369 }
0370 }
0371 }
0372 }
0373 }