Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-11 04:32:50 

0001 /* This Class Header */
0002 #include "DQMOffline/Muon/interface/DiMuonHistograms.h"
0003 
0004 /* Collaborating Class Header */
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 /* C++ Headers */
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   // initialise parameters:
0033   parameters = pSet;
0034 
0035   // counter
0036   nTightTight = 0;
0037   nMediumMedium = 0;
0038   nLooseLoose = 0;
0039   nGlbGlb = 0;
0040 
0041   // declare consumes:
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& /*iRun*/,
0068                                       edm::EventSetup const& /* iSetup */) {
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     // arround the Z peak
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     //Fraction of bad hits in the tracker track to the total
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     // low-mass resonances
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   // Look for the Primary Vertex (and use the BeamSpot instead, if you can't find it):
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   // Loop on muon collection
0186   TLorentzVector Mu1, Mu2;
0187   float charge = 99.;
0188   float InvMass = -99.;
0189 
0190   //Eta regions
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     // Loop on second muons to fill invariant mass plots
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       // Global-Global Muon
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         // Also Tight-Tight Muon Selection
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         // Also Medium-Medium Muon Selection
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         // Also Loose-Loose Muon Selection
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       // Now check for STA-TRK
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       // TRK-TRK dimuon
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     }  //muon2
0372   }  //Muon1
0373 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!