Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:08:02

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author S. Bolognesi, Erik - CERN
0005  */
0006 
0007 #include "DQM/Physics/src/BPhysicsOniaDQM.h"
0008 
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 
0013 #include "DQMServices/Core/interface/DQMStore.h"
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 
0016 #include "DataFormats/Common/interface/Handle.h"
0017 #include "DataFormats/MuonReco/interface/Muon.h"
0018 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0019 #include "DataFormats/TrackReco/interface/Track.h"
0020 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0021 
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 
0024 using namespace std;
0025 using namespace edm;
0026 using namespace reco;
0027 
0028 BPhysicsOniaDQM::BPhysicsOniaDQM(const ParameterSet &parameters) {
0029   // Muon Collection Label
0030   vertex_ = consumes<reco::VertexCollection>(parameters.getParameter<InputTag>("vertex"));
0031   theMuonCollectionLabel_ = consumes<reco::MuonCollection>(parameters.getParameter<InputTag>("MuonCollection"));
0032   lumiSummaryToken_ = consumes<LumiSummary, edm::InLumi>(parameters.getParameter<InputTag>("lumiSummary"));
0033 
0034   global_background = nullptr;
0035   diMuonMass_global = nullptr;
0036   tracker_background = nullptr;
0037   diMuonMass_tracker = nullptr;
0038   standalone_background = nullptr;
0039   diMuonMass_standalone = nullptr;
0040 
0041   glbSigCut = nullptr;
0042   glbSigNoCut = nullptr;
0043   glbBkgNoCut = nullptr;
0044   staSigCut = nullptr;
0045   staSigNoCut = nullptr;
0046   staBkgNoCut = nullptr;
0047   trkSigCut = nullptr;
0048   trkSigNoCut = nullptr;
0049   trkBkgNoCut = nullptr;
0050 
0051   metname = "oniaAnalyzer";
0052 }
0053 
0054 BPhysicsOniaDQM::~BPhysicsOniaDQM() {}
0055 
0056 void BPhysicsOniaDQM::bookHistograms(DQMStore::IBooker &iBooker, edm::Run const &, edm::EventSetup const &) {
0057   iBooker.setCurrentFolder("Physics/BPhysics");  // Use folder with name of PAG
0058 
0059   global_background = iBooker.book1D("global_background", "Same-sign global-global dimuon mass", 750, 0, 15);
0060   diMuonMass_global = iBooker.book1D("diMuonMass_global", "Opposite-sign global-global dimuon mass", 750, 0, 15);
0061   tracker_background =
0062       iBooker.book1D("tracker_background", "Same-sign tracker-tracker (arbitrated) dimuon mass", 750, 0, 15);
0063   diMuonMass_tracker =
0064       iBooker.book1D("diMuonMass_tracker", "Opposite-sign tracker-tracker (arbitrated) dimuon mass", 750, 0, 15);
0065   standalone_background =
0066       iBooker.book1D("standalone_background", "Same-sign standalone-standalone dimuon mass", 500, 0, 15);
0067   diMuonMass_standalone =
0068       iBooker.book1D("diMuonMass_standalone", "Opposite-sign standalone-standalone dimuon mass", 500, 0, 15);
0069 
0070   glbSigCut = iBooker.book1D("glbSigCut", "Opposite-sign glb-glb dimuon mass", 650, 0, 130);
0071   glbSigNoCut = iBooker.book1D("glbSigNoCut", "Opposite-sign glb-glb dimuon mass (no cut)", 650, 0, 130);
0072   glbBkgNoCut = iBooker.book1D("glbBkgNoCut", "Same-sign glb-glb dimuon mass (no cut)", 650, 0, 130);
0073   staSigCut = iBooker.book1D("staSigCut", "Opposite-sign sta-sta dimuon mass", 430, 0, 129);
0074   staSigNoCut = iBooker.book1D("staSigNoCut", "Opposite-sign sta-sta dimuon mass (no cut)", 430, 0, 129);
0075   staBkgNoCut = iBooker.book1D("staBkgNoCut", "Same-sign sta-sta dimuon mass (no cut)", 430, 0, 129);
0076   trkSigCut = iBooker.book1D("trkSigCut", "Opposite-sign trk-trk dimuon mass", 650, 0, 130);
0077   trkSigNoCut = iBooker.book1D("trkSigNoCut", "Opposite-sign trk-trk dimuon mass (no cut)", 650, 0, 130);
0078   trkBkgNoCut = iBooker.book1D("trkBkgNoCutt", "Same-sign trk-trk dimuon mass (no cut)", 650, 0, 130);
0079 }
0080 
0081 void BPhysicsOniaDQM::analyze(const Event &iEvent, const EventSetup &iSetup) {
0082   LogTrace(metname) << "[BPhysicsOniaDQM] Analysis of event # ";
0083 
0084   // Take the STA muon container
0085   Handle<MuonCollection> muons;
0086   iEvent.getByToken(theMuonCollectionLabel_, muons);
0087 
0088   Handle<reco::VertexCollection> privtxs;
0089   iEvent.getByToken(vertex_, privtxs);
0090   VertexCollection::const_iterator privtx;
0091 
0092   if (privtxs->begin() != privtxs->end()) {
0093     privtx = privtxs->begin();
0094     RefVtx = privtx->position();
0095   } else {
0096     RefVtx.SetXYZ(0., 0., 0.);
0097   }
0098 
0099   if (muons.isValid()) {
0100     for (MuonCollection::const_iterator recoMu1 = muons->begin(); recoMu1 != muons->end(); ++recoMu1) {
0101       // only loop over the remaining muons if recoMu1 is one of the following
0102       if (recoMu1->isGlobalMuon() || recoMu1->isTrackerMuon() || recoMu1->isStandAloneMuon()) {
0103         for (MuonCollection::const_iterator recoMu2 = recoMu1 + 1; recoMu2 != muons->end(); ++recoMu2) {
0104           // fill the relevant histograms if recoMu2 satisfies one of the
0105           // following
0106           if (recoMu1->isGlobalMuon() && recoMu2->isGlobalMuon()) {
0107             math::XYZVector vec1 = recoMu1->globalTrack()->momentum();
0108             math::XYZVector vec2 = recoMu2->globalTrack()->momentum();
0109             float massJPsi = computeMass(vec1, vec2);
0110 
0111             // if opposite charges, fill glbSig, else fill glbBkg
0112             if (((*recoMu1).charge() * (*recoMu2).charge()) < 0) {
0113               if (diMuonMass_global != nullptr) {  // BPhysicsOniaDQM original one
0114                 diMuonMass_global->Fill(massJPsi);
0115               }
0116 
0117               if (glbSigNoCut != nullptr) {
0118                 glbSigNoCut->Fill(massJPsi);
0119                 if (selGlobalMuon(*recoMu1) && selGlobalMuon(*recoMu2)) {
0120                   if (glbSigCut != nullptr)
0121                     glbSigCut->Fill(massJPsi);
0122                 }
0123               }
0124             } else {
0125               if (global_background != nullptr) {  // BPhysicsOniaDQM original one
0126                 global_background->Fill(massJPsi);
0127               }
0128 
0129               if (glbBkgNoCut != nullptr) {
0130                 glbBkgNoCut->Fill(massJPsi);
0131               }
0132             }
0133           }
0134 
0135           if (recoMu1->isStandAloneMuon() && recoMu2->isStandAloneMuon() && fabs(recoMu1->outerTrack()->d0()) < 5 &&
0136               fabs(recoMu1->outerTrack()->dz()) < 30 && fabs(recoMu2->outerTrack()->d0()) < 5 &&
0137               fabs(recoMu2->outerTrack()->dz()) < 30) {
0138             math::XYZVector vec1 = recoMu1->outerTrack()->momentum();
0139             math::XYZVector vec2 = recoMu2->outerTrack()->momentum();
0140             float massJPsi = computeMass(vec1, vec2);
0141 
0142             // if opposite charges, fill staSig, else fill staBkg
0143             if (((*recoMu1).charge() * (*recoMu2).charge()) < 0) {
0144               if (diMuonMass_standalone != nullptr) {
0145                 diMuonMass_standalone->Fill(massJPsi);
0146               }
0147 
0148               if (staSigNoCut != nullptr) {
0149                 staSigNoCut->Fill(massJPsi);
0150               }
0151             } else {
0152               if (standalone_background != nullptr) {
0153                 standalone_background->Fill(massJPsi);
0154               }
0155 
0156               if (staBkgNoCut != nullptr) {
0157                 staBkgNoCut->Fill(massJPsi);
0158               }
0159             }
0160           }
0161 
0162           if (recoMu1->isTrackerMuon() && recoMu2->isTrackerMuon() &&
0163               muon::isGoodMuon(*recoMu1, muon::TrackerMuonArbitrated) &&
0164               muon::isGoodMuon(*recoMu2, muon::TrackerMuonArbitrated)) {
0165             math::XYZVector vec1 = recoMu1->innerTrack()->momentum();
0166             math::XYZVector vec2 = recoMu2->innerTrack()->momentum();
0167             float massJPsi = computeMass(vec1, vec2);
0168 
0169             // if opposite charges, fill trkSig, else fill trkBkg
0170             if (((*recoMu1).charge() * (*recoMu2).charge()) < 0) {
0171               if (diMuonMass_tracker != nullptr) {
0172                 diMuonMass_tracker->Fill(massJPsi);
0173               }
0174 
0175               if (trkSigNoCut != nullptr) {
0176                 trkSigNoCut->Fill(massJPsi);
0177                 if (selTrackerMuon(*recoMu1) && selTrackerMuon(*recoMu2)) {
0178                   if (trkSigCut != nullptr)
0179                     trkSigCut->Fill(massJPsi);
0180                 }
0181               }
0182             } else {
0183               if (tracker_background != nullptr) {
0184                 tracker_background->Fill(massJPsi);
0185               }
0186 
0187               if (trkBkgNoCut != nullptr) {
0188                 trkBkgNoCut->Fill(massJPsi);
0189               }
0190             }
0191           }
0192         }  // end of 2nd MuonCollection
0193       }    // end of GLB,STA,TRK muon check
0194     }      // end of 1st MuonCollection
0195   }        // Is this MuonCollection vaild?
0196 }
0197 
0198 float BPhysicsOniaDQM::computeMass(const math::XYZVector &vec1, const math::XYZVector &vec2) {
0199   // mass of muon
0200   float massMu = 0.10566;
0201   float eMu1 = -999;
0202   if (massMu * massMu + vec1.Mag2() > 0)
0203     eMu1 = sqrt(massMu * massMu + vec1.Mag2());
0204   float eMu2 = -999;
0205   if (massMu * massMu + vec2.Mag2() > 0)
0206     eMu2 = sqrt(massMu * massMu + vec2.Mag2());
0207 
0208   float pJPsi = -999;
0209   if ((vec1 + vec2).Mag2() > 0)
0210     pJPsi = sqrt((vec1 + vec2).Mag2());
0211   float eJPsi = eMu1 + eMu2;
0212 
0213   float massJPsi = -999;
0214   if ((eJPsi * eJPsi - pJPsi * pJPsi) > 0)
0215     massJPsi = sqrt(eJPsi * eJPsi - pJPsi * pJPsi);
0216 
0217   return massJPsi;
0218 }
0219 
0220 bool BPhysicsOniaDQM::isMuonInAccept(const reco::Muon &recoMu) {
0221   return (fabs(recoMu.eta()) < 2.4 && ((fabs(recoMu.eta()) < 1.3 && recoMu.pt() > 3.3) ||
0222                                        (fabs(recoMu.eta()) > 1.3 && fabs(recoMu.eta()) < 2.2 && recoMu.p() > 2.9) ||
0223                                        (fabs(recoMu.eta()) > 2.2 && recoMu.pt() > 0.8)));
0224 }
0225 
0226 bool BPhysicsOniaDQM::selGlobalMuon(const reco::Muon &recoMu) {
0227   TrackRef iTrack = recoMu.innerTrack();
0228   const reco::HitPattern &p = iTrack->hitPattern();
0229 
0230   TrackRef gTrack = recoMu.globalTrack();
0231   const reco::HitPattern &q = gTrack->hitPattern();
0232 
0233   return (isMuonInAccept(recoMu) && iTrack->found() > 11 && gTrack->chi2() / gTrack->ndof() < 20.0 &&
0234           q.numberOfValidMuonHits() > 0 && iTrack->chi2() / iTrack->ndof() < 4.0 &&
0235           // recoMu.muonID("TrackerMuonArbitrated") &&
0236           // recoMu.muonID("TMLastStationAngTight") &&
0237           p.pixelLayersWithMeasurement() > 1 && fabs(iTrack->dxy(RefVtx)) < 3.0 && fabs(iTrack->dz(RefVtx)) < 15.0);
0238 }
0239 
0240 bool BPhysicsOniaDQM::selTrackerMuon(const reco::Muon &recoMu) {
0241   TrackRef iTrack = recoMu.innerTrack();
0242   const reco::HitPattern &p = iTrack->hitPattern();
0243 
0244   return (isMuonInAccept(recoMu) && iTrack->found() > 11 && iTrack->chi2() / iTrack->ndof() < 4.0 &&
0245           // recoMu.muonID("TrackerMuonArbitrated") &&
0246           // recoMu.muonID("TMLastStationAngTight") &&
0247           p.pixelLayersWithMeasurement() > 1 && fabs(iTrack->dxy(RefVtx)) < 3.0 && fabs(iTrack->dz(RefVtx)) < 15.0);
0248 }
0249 
0250 // Local Variables:
0251 // show-trailing-whitespace: t
0252 // truncate-lines: t
0253 // End: