CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
DiMuonHistograms.cc
Go to the documentation of this file.
1 /* This Class Header */
3 
4 /* Collaborating Class Header */
16 
18 
19 #include "TLorentzVector.h"
20 #include "TFile.h"
21 #include <vector>
22 #include <cmath>
23 
24 /* C++ Headers */
25 #include <iostream>
26 #include <fstream>
27 #include <cmath>
28 using namespace std;
29 using namespace edm;
30 
32  // initialise parameters:
33  parameters = pSet;
34 
35  // counter
36  nTightTight = 0;
37  nMediumMedium = 0;
38  nLooseLoose = 0;
39  nGlbGlb = 0;
40 
41  // declare consumes:
42  theMuonCollectionLabel_ = consumes<edm::View<reco::Muon> >(parameters.getParameter<edm::InputTag>("MuonCollection"));
43  theVertexLabel_ = consumes<reco::VertexCollection>(parameters.getParameter<edm::InputTag>("VertexLabel"));
44 
45  theBeamSpotLabel_ = mayConsume<reco::BeamSpot>(parameters.getParameter<edm::InputTag>("BeamSpotLabel"));
46 
47  etaBin = parameters.getParameter<int>("etaBin");
48  etaBBin = parameters.getParameter<int>("etaBBin");
49  etaEBin = parameters.getParameter<int>("etaEBin");
50 
51  etaBMin = parameters.getParameter<double>("etaBMin");
52  etaBMax = parameters.getParameter<double>("etaBMax");
53  etaECMin = parameters.getParameter<double>("etaECMin");
54  etaECMax = parameters.getParameter<double>("etaECMax");
55 
56  LowMassMin = parameters.getParameter<double>("LowMassMin");
57  LowMassMax = parameters.getParameter<double>("LowMassMax");
58  HighMassMin = parameters.getParameter<double>("HighMassMin");
59  HighMassMax = parameters.getParameter<double>("HighMassMax");
60 
61  theFolder = parameters.getParameter<string>("folder");
62 }
63 
65 
67  edm::Run const& /*iRun*/,
68  edm::EventSetup const& /* iSetup */) {
69  ibooker.cd();
70  ibooker.setCurrentFolder(theFolder);
71 
72  int nBin[3] = {etaBin, etaBBin, etaEBin};
73  EtaName[0] = "";
74  EtaName[1] = "_Barrel";
75  EtaName[2] = "_EndCap";
76  test = ibooker.book1D("test", "InvMass_{Tight,Tight}", 100, 0., 200.);
77  for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
78  GlbGlbMuon_LM.push_back(ibooker.book1D("GlbGlbMuon_LM" + EtaName[iEtaRegion],
79  "InvMass_{GLB,GLB}" + EtaName[iEtaRegion],
80  nBin[iEtaRegion],
81  LowMassMin,
82  LowMassMax));
83  TrkTrkMuon_LM.push_back(ibooker.book1D("TrkTrkMuon_LM" + EtaName[iEtaRegion],
84  "InvMass_{TRK,TRK}" + EtaName[iEtaRegion],
85  nBin[iEtaRegion],
86  LowMassMin,
87  LowMassMax));
88  StaTrkMuon_LM.push_back(ibooker.book1D("StaTrkMuon_LM" + EtaName[iEtaRegion],
89  "InvMass_{STA,TRK}" + EtaName[iEtaRegion],
90  nBin[iEtaRegion],
91  LowMassMin,
92  LowMassMax));
93 
94  GlbGlbMuon_HM.push_back(ibooker.book1D("GlbGlbMuon_HM" + EtaName[iEtaRegion],
95  "InvMass_{GLB,GLB}" + EtaName[iEtaRegion],
96  nBin[iEtaRegion],
98  HighMassMax));
99  TrkTrkMuon_HM.push_back(ibooker.book1D("TrkTrkMuon_HM" + EtaName[iEtaRegion],
100  "InvMass_{TRK,TRK}" + EtaName[iEtaRegion],
101  nBin[iEtaRegion],
102  HighMassMin,
103  HighMassMax));
104  StaTrkMuon_HM.push_back(ibooker.book1D("StaTrkMuon_HM" + EtaName[iEtaRegion],
105  "InvMass_{STA,TRK}" + EtaName[iEtaRegion],
106  nBin[iEtaRegion],
107  HighMassMin,
108  HighMassMax));
109 
110  // arround the Z peak
111  TightTightMuon.push_back(ibooker.book1D("TightTightMuon" + EtaName[iEtaRegion],
112  "InvMass_{Tight,Tight}" + EtaName[iEtaRegion],
113  nBin[iEtaRegion],
114  LowMassMin,
115  LowMassMax));
116  MediumMediumMuon.push_back(ibooker.book1D("MediumMediumMuon" + EtaName[iEtaRegion],
117  "InvMass_{Medium,Medium}" + EtaName[iEtaRegion],
118  nBin[iEtaRegion],
119  LowMassMin,
120  LowMassMax));
121  LooseLooseMuon.push_back(ibooker.book1D("LooseLooseMuon" + EtaName[iEtaRegion],
122  "InvMass_{Loose,Loose}" + EtaName[iEtaRegion],
123  nBin[iEtaRegion],
124  LowMassMin,
125  LowMassMax));
126  //Fraction of bad hits in the tracker track to the total
127  TightTightMuonBadFrac.push_back(ibooker.book1D(
128  "TightTightMuonBadFrac" + EtaName[iEtaRegion], "BadFrac_{Tight,Tight}" + EtaName[iEtaRegion], 10, 0, 0.4));
129  MediumMediumMuonBadFrac.push_back(ibooker.book1D(
130  "MediumMediumMuonBadFrac" + EtaName[iEtaRegion], "BadFrac_{Medium,Medium}" + EtaName[iEtaRegion], 10, 0, 0.4));
131  LooseLooseMuonBadFrac.push_back(ibooker.book1D(
132  "LooseLooseMuonBadFrac" + EtaName[iEtaRegion], "BadFrac_{Loose,Loose}" + EtaName[iEtaRegion], 10, 0, 0.4));
133 
134  // low-mass resonances
135  SoftSoftMuon.push_back(ibooker.book1D(
136  "SoftSoftMuon" + EtaName[iEtaRegion], "InvMass_{Soft,Soft}" + EtaName[iEtaRegion], nBin[iEtaRegion], 0.0, 55.0));
137  SoftSoftMuonBadFrac.push_back(ibooker.book1D(
138  "SoftSoftMuonBadFrac" + EtaName[iEtaRegion], "BadFrac_{Soft,Soft}" + EtaName[iEtaRegion], 10, 0, 0.4));
139  }
140 }
141 
143  LogTrace(metname) << "[DiMuonHistograms] Analyze the mu in different eta regions";
145  iEvent.getByToken(theMuonCollectionLabel_, muons);
146 
147  // =================================================================================
148  // Look for the Primary Vertex (and use the BeamSpot instead, if you can't find it):
149  reco::Vertex::Point posVtx;
150  reco::Vertex::Error errVtx;
151  unsigned int theIndexOfThePrimaryVertex = 999.;
152 
154  iEvent.getByToken(theVertexLabel_, vertex);
155  if (vertex.isValid()) {
156  for (unsigned int ind = 0; ind < vertex->size(); ++ind) {
157  if ((*vertex)[ind].isValid() && !((*vertex)[ind].isFake())) {
158  theIndexOfThePrimaryVertex = ind;
159  break;
160  }
161  }
162  }
163 
164  if (theIndexOfThePrimaryVertex < 100) {
165  posVtx = ((*vertex)[theIndexOfThePrimaryVertex]).position();
166  errVtx = ((*vertex)[theIndexOfThePrimaryVertex]).error();
167  } else {
168  LogInfo("RecoMuonValidator") << "reco::PrimaryVertex not found, use BeamSpot position instead\n";
169 
170  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
171  iEvent.getByToken(theBeamSpotLabel_, recoBeamSpotHandle);
172  reco::BeamSpot bs = *recoBeamSpotHandle;
173 
174  posVtx = bs.position();
175  errVtx(0, 0) = bs.BeamWidthX();
176  errVtx(1, 1) = bs.BeamWidthY();
177  errVtx(2, 2) = bs.sigmaZ();
178  }
179 
180  const reco::Vertex vtx(posVtx, errVtx);
181 
182  if (!muons.isValid())
183  return;
184 
185  // Loop on muon collection
186  TLorentzVector Mu1, Mu2;
187  float charge = 99.;
188  float InvMass = -99.;
189 
190  //Eta regions
191  double EtaCutMin[] = {0, etaBMin, etaECMin};
192  double EtaCutMax[] = {2.4, etaBMax, etaECMax};
193 
194  for (edm::View<reco::Muon>::const_iterator muon1 = muons->begin(); muon1 != muons->end(); ++muon1) {
195  LogTrace(metname) << "[DiMuonHistograms] loop over 1st muon" << endl;
196 
197  // Loop on second muons to fill invariant mass plots
198  for (edm::View<reco::Muon>::const_iterator muon2 = muon1; muon2 != muons->end(); ++muon2) {
199  LogTrace(metname) << "[DiMuonHistograms] loop over 2nd muon" << endl;
200  if (muon1 == muon2)
201  continue;
202 
203  // Global-Global Muon
204  if (muon1->isGlobalMuon() && muon2->isGlobalMuon()) {
205  LogTrace(metname) << "[DiMuonHistograms] Glb-Glb pair" << endl;
206  reco::TrackRef recoCombinedGlbTrack1 = muon1->combinedMuon();
207  reco::TrackRef recoCombinedGlbTrack2 = muon2->combinedMuon();
208  Mu1.SetPxPyPzE(recoCombinedGlbTrack1->px(),
209  recoCombinedGlbTrack1->py(),
210  recoCombinedGlbTrack1->pz(),
211  recoCombinedGlbTrack1->p());
212  Mu2.SetPxPyPzE(recoCombinedGlbTrack2->px(),
213  recoCombinedGlbTrack2->py(),
214  recoCombinedGlbTrack2->pz(),
215  recoCombinedGlbTrack2->p());
216 
217  charge = recoCombinedGlbTrack1->charge() * recoCombinedGlbTrack2->charge();
218  if (charge < 0) {
219  InvMass = (Mu1 + Mu2).M();
220  for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
221  if (fabs(recoCombinedGlbTrack1->eta()) > EtaCutMin[iEtaRegion] &&
222  fabs(recoCombinedGlbTrack1->eta()) < EtaCutMax[iEtaRegion] &&
223  fabs(recoCombinedGlbTrack2->eta()) > EtaCutMin[iEtaRegion] &&
224  fabs(recoCombinedGlbTrack2->eta()) < EtaCutMax[iEtaRegion]) {
225  if (InvMass < LowMassMax)
226  GlbGlbMuon_LM[iEtaRegion]->Fill(InvMass);
227  if (InvMass > HighMassMin)
228  GlbGlbMuon_HM[iEtaRegion]->Fill(InvMass);
229  }
230  }
231  }
232  // Also Tight-Tight Muon Selection
233  if (muon::isTightMuon(*muon1, vtx) && muon::isTightMuon(*muon2, vtx)) {
234  test->Fill(InvMass);
235  LogTrace(metname) << "[DiMuonHistograms] Tight-Tight pair" << endl;
236  if (charge < 0) {
237  for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
238  if (fabs(recoCombinedGlbTrack1->eta()) > EtaCutMin[iEtaRegion] &&
239  fabs(recoCombinedGlbTrack1->eta()) < EtaCutMax[iEtaRegion] &&
240  fabs(recoCombinedGlbTrack2->eta()) > EtaCutMin[iEtaRegion] &&
241  fabs(recoCombinedGlbTrack2->eta()) < EtaCutMax[iEtaRegion]) {
242  if (InvMass > 55. && InvMass < 125.) {
243  TightTightMuon[iEtaRegion]->Fill(InvMass);
244  TightTightMuonBadFrac[iEtaRegion]->Fill(muon1->innerTrack()->lost() / muon1->innerTrack()->found());
245  }
246  }
247  }
248  }
249  }
250  // Also Medium-Medium Muon Selection
251  if (muon::isMediumMuon(*muon1) && muon::isMediumMuon(*muon2)) {
252  test->Fill(InvMass);
253  LogTrace(metname) << "[DiMuonHistograms] Medium-Medium pair" << endl;
254  if (charge < 0) {
255  for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
256  if (fabs(recoCombinedGlbTrack1->eta()) > EtaCutMin[iEtaRegion] &&
257  fabs(recoCombinedGlbTrack1->eta()) < EtaCutMax[iEtaRegion] &&
258  fabs(recoCombinedGlbTrack2->eta()) > EtaCutMin[iEtaRegion] &&
259  fabs(recoCombinedGlbTrack2->eta()) < EtaCutMax[iEtaRegion]) {
260  if (InvMass > 55. && InvMass < 125.) {
261  MediumMediumMuon[iEtaRegion]->Fill(InvMass);
262  MediumMediumMuonBadFrac[iEtaRegion]->Fill(muon1->innerTrack()->lost() / muon1->innerTrack()->found());
263  }
264  }
265  }
266  }
267  }
268  // Also Loose-Loose Muon Selection
269  if (muon::isLooseMuon(*muon1) && muon::isLooseMuon(*muon2)) {
270  test->Fill(InvMass);
271  LogTrace(metname) << "[DiMuonHistograms] Loose-Loose pair" << endl;
272  if (charge < 0) {
273  for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
274  if (fabs(recoCombinedGlbTrack1->eta()) > EtaCutMin[iEtaRegion] &&
275  fabs(recoCombinedGlbTrack1->eta()) < EtaCutMax[iEtaRegion] &&
276  fabs(recoCombinedGlbTrack2->eta()) > EtaCutMin[iEtaRegion] &&
277  fabs(recoCombinedGlbTrack2->eta()) < EtaCutMax[iEtaRegion]) {
278  if (InvMass > 55. && InvMass < 125.) {
279  LooseLooseMuon[iEtaRegion]->Fill(InvMass);
280  LooseLooseMuonBadFrac[iEtaRegion]->Fill(muon1->innerTrack()->lost() / muon1->innerTrack()->found());
281  }
282  }
283  }
284  }
285  }
286  }
287 
288  // Now check for STA-TRK
289  if (muon2->isStandAloneMuon() && muon1->isTrackerMuon()) {
290  LogTrace(metname) << "[DiMuonHistograms] STA-Trk pair" << endl;
291  reco::TrackRef recoStaTrack = muon2->standAloneMuon();
292  reco::TrackRef recoTrack = muon1->track();
293  Mu2.SetPxPyPzE(recoStaTrack->px(), recoStaTrack->py(), recoStaTrack->pz(), recoStaTrack->p());
294  Mu1.SetPxPyPzE(recoTrack->px(), recoTrack->py(), recoTrack->pz(), recoTrack->p());
295 
296  charge = recoStaTrack->charge() * recoTrack->charge();
297  if (charge < 0) {
298  InvMass = (Mu1 + Mu2).M();
299  for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
300  if (fabs(recoStaTrack->eta()) > EtaCutMin[iEtaRegion] &&
301  fabs(recoStaTrack->eta()) < EtaCutMax[iEtaRegion] && fabs(recoTrack->eta()) > EtaCutMin[iEtaRegion] &&
302  fabs(recoTrack->eta()) < EtaCutMax[iEtaRegion]) {
303  if (InvMass < LowMassMax)
304  StaTrkMuon_LM[iEtaRegion]->Fill(InvMass);
305  if (InvMass > HighMassMin)
306  StaTrkMuon_HM[iEtaRegion]->Fill(InvMass);
307  }
308  }
309  }
310  }
311  if (muon1->isStandAloneMuon() && muon2->isTrackerMuon()) {
312  LogTrace(metname) << "[DiMuonHistograms] STA-Trk pair" << endl;
313  reco::TrackRef recoStaTrack = muon1->standAloneMuon();
314  reco::TrackRef recoTrack = muon2->track();
315  Mu1.SetPxPyPzE(recoStaTrack->px(), recoStaTrack->py(), recoStaTrack->pz(), recoStaTrack->p());
316  Mu2.SetPxPyPzE(recoTrack->px(), recoTrack->py(), recoTrack->pz(), recoTrack->p());
317 
318  charge = recoStaTrack->charge() * recoTrack->charge();
319  if (charge < 0) {
320  InvMass = (Mu1 + Mu2).M();
321  for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
322  if (fabs(recoStaTrack->eta()) > EtaCutMin[iEtaRegion] &&
323  fabs(recoStaTrack->eta()) < EtaCutMax[iEtaRegion] && fabs(recoTrack->eta()) > EtaCutMin[iEtaRegion] &&
324  fabs(recoTrack->eta()) < EtaCutMax[iEtaRegion]) {
325  if (InvMass < LowMassMax)
326  StaTrkMuon_LM[iEtaRegion]->Fill(InvMass);
327  if (InvMass > HighMassMin)
328  StaTrkMuon_HM[iEtaRegion]->Fill(InvMass);
329  }
330  }
331  }
332  }
333 
334  // TRK-TRK dimuon
335  if (muon1->isTrackerMuon() && muon2->isTrackerMuon()) {
336  LogTrace(metname) << "[DiMuonHistograms] Trk-Trk dimuon pair" << endl;
337  reco::TrackRef recoTrack2 = muon2->track();
338  reco::TrackRef recoTrack1 = muon1->track();
339  Mu2.SetPxPyPzE(recoTrack2->px(), recoTrack2->py(), recoTrack2->pz(), recoTrack2->p());
340  Mu1.SetPxPyPzE(recoTrack1->px(), recoTrack1->py(), recoTrack1->pz(), recoTrack1->p());
341 
342  charge = recoTrack1->charge() * recoTrack2->charge();
343  if (charge < 0) {
344  InvMass = (Mu1 + Mu2).M();
345  for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
346  if (fabs(recoTrack1->eta()) > EtaCutMin[iEtaRegion] && fabs(recoTrack1->eta()) < EtaCutMax[iEtaRegion] &&
347  fabs(recoTrack2->eta()) > EtaCutMin[iEtaRegion] && fabs(recoTrack2->eta()) < EtaCutMax[iEtaRegion]) {
348  if (InvMass < LowMassMax)
349  TrkTrkMuon_LM[iEtaRegion]->Fill(InvMass);
350  if (InvMass > HighMassMin)
351  TrkTrkMuon_HM[iEtaRegion]->Fill(InvMass);
352  }
353  }
354  }
355 
356  LogTrace(metname) << "[DiMuonHistograms] Soft-Soft pair" << endl;
357 
358  if (muon::isSoftMuon(*muon1, vtx) && muon::isSoftMuon(*muon2, vtx)) {
359  if (charge < 0) {
360  InvMass = (Mu1 + Mu2).M();
361  for (unsigned int iEtaRegion = 0; iEtaRegion < 3; iEtaRegion++) {
362  if (fabs(recoTrack1->eta()) > EtaCutMin[iEtaRegion] && fabs(recoTrack1->eta()) < EtaCutMax[iEtaRegion] &&
363  fabs(recoTrack2->eta()) > EtaCutMin[iEtaRegion] && fabs(recoTrack2->eta()) < EtaCutMax[iEtaRegion]) {
364  SoftSoftMuon[iEtaRegion]->Fill(InvMass);
365  SoftSoftMuonBadFrac[iEtaRegion]->Fill(muon1->innerTrack()->lost() / muon1->innerTrack()->found());
366  }
367  }
368  }
369  }
370  }
371  } //muon2
372  } //Muon1
373 }
~DiMuonHistograms() override
void analyze(const edm::Event &, const edm::EventSetup &) override
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:32
bool isMediumMuon(const reco::Muon &, bool run2016_hip_mitigation=false)
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:539
const std::string metname
math::Error< dimension >::type Error
covariance error matrix (3x3)
Definition: Vertex.h:44
#define LogTrace(id)
bool isLooseMuon(const reco::Muon &)
int iEvent
Definition: GenABIO.cc:224
math::XYZPoint Point
point in the space
Definition: Vertex.h:40
double BeamWidthX() const
beam width X
Definition: BeamSpot.h:82
DiMuonHistograms(const edm::ParameterSet &pset)
bool isSoftMuon(const reco::Muon &, const reco::Vertex &, bool run2016_hip_mitigation=false)
bool isValid() const
Definition: HandleBase.h:70
Log< level::Info, false > LogInfo
double sigmaZ() const
sigma z
Definition: BeamSpot.h:76
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:84
tuple muons
Definition: patZpeak.py:39
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:86
static int position[264][3]
Definition: ReadPGInfo.cc:289
const Point & position() const
position
Definition: BeamSpot.h:59
bool isTightMuon(const reco::Muon &, const reco::Vertex &)
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
int etaBin(const l1t::HGCalMulticluster *cl)
Definition: Run.h:45