CMS 3D CMS Logo
RecoMuonValidator.cc
Go to the documentation of this file.
1 #include "Validation/RecoMuon/src/RecoMuonValidator.h"
2 
3 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
4 
5 #include "FWCore/ServiceRegistry/interface/Service.h"
6 #include "FWCore/MessageLogger/interface/MessageLogger.h"
7 //#include "DQMServices/Core/interface/DQMStore.h"
8 
9 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
10 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
11 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
12 
13 // for selection cut
14 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
15 
16 #include "TMath.h"
17 
18 using namespace std;
19 using namespace edm;
20 using namespace reco;
21 
22 typedef TrajectoryStateOnSurface TSOS;
23 typedef FreeTrajectoryState FTS;
24 
25 //
26 //Struct containing all histograms definitions
27 //
28 struct RecoMuonValidator::MuonME {
29  typedef MonitorElement* MEP;
30 
31  //general kinematics
32  MEP hSimP_, hSimPt_, hSimEta_, hSimPhi_, hSimDxy_, hSimDz_;
33  //only for efficiencies
34  MEP hP_, hPt_, hEta_, hPhi_;
35  MEP hNSim_, hNMuon_;
36 
37  //misc vars
38  MEP hNTrks_, hNTrksEta_, hNTrksPt_;
39  MEP hMisQPt_, hMisQEta_;
40 
41  //resolutions
42  MEP hErrP_, hErrPt_, hErrEta_, hErrPhi_;
43  MEP hErrPBarrel_, hErrPOverlap_, hErrPEndcap_;
44  MEP hErrPtBarrel_, hErrPtOverlap_, hErrPtEndcap_;
45  MEP hErrDxy_, hErrDz_;
46 
47  MEP hErrP_vs_Eta_, hErrPt_vs_Eta_, hErrQPt_vs_Eta_;
48  MEP hErrP_vs_P_, hErrPt_vs_Pt_, hErrQPt_vs_Pt_, hErrEta_vs_Eta_;
49 
50  //PF-RECO event-by-event comparisons
51  MEP hErrPt_PF_;
52  MEP hErrQPt_PF_;
53  MEP hdPt_vs_Eta_;
54  MEP hdPt_vs_Pt_;
55  MEP hPFMomAssCorrectness;
56  MEP hPt_vs_PFMomAssCorrectness;
57 
58  //hit pattern
59  MEP hNSimHits_;
60  MEP hNSimToReco_, hNRecoToSim_;
61 
62  MEP hNHits_, hNLostHits_, hNTrackerHits_, hNMuonHits_;
63  MEP hNHits_vs_Pt_, hNHits_vs_Eta_;
64  MEP hNLostHits_vs_Pt_, hNLostHits_vs_Eta_;
65  MEP hNTrackerHits_vs_Pt_, hNTrackerHits_vs_Eta_;
66  MEP hNMuonHits_vs_Pt_, hNMuonHits_vs_Eta_;
67 
68  //pulls
69  MEP hPullPt_, hPullEta_, hPullPhi_, hPullQPt_, hPullDxy_, hPullDz_;
70  MEP hPullPt_vs_Eta_, hPullPt_vs_Pt_, hPullEta_vs_Eta_, hPullPhi_vs_Eta_, hPullEta_vs_Pt_;
71 
72  //chi2, ndof
73  MEP hNDof_, hChi2_, hChi2Norm_, hChi2Prob_;
74  MEP hNDof_vs_Eta_, hChi2_vs_Eta_, hChi2Norm_vs_Eta_, hChi2Prob_vs_Eta_;
75 
76  bool doAbsEta_;
77  bool usePFMuon_;
78 
79  //
80  //books histograms
81  //
82  void bookHistos(DQMStore::IBooker& ibooker, const string& dirName, const HistoDimensions& hDim)
83 
84  {
85  ibooker.cd();
86  ibooker.setCurrentFolder(dirName);
87 
88  doAbsEta_ = hDim.doAbsEta;
89  usePFMuon_ = hDim.usePFMuon;
90 
91  //histograms for efficiency plots
92  hP_ = ibooker.book1D("P", "p of recoTracks", hDim.nBinP, hDim.minP, hDim.maxP);
93  hPt_ = ibooker.book1D("Pt", "p_{T} of recoTracks", hDim.nBinPt, hDim.minPt, hDim.maxPt);
94  hEta_ = ibooker.book1D("Eta", "#eta of recoTracks", hDim.nBinEta, hDim.minEta, hDim.maxEta);
95  hPhi_ = ibooker.book1D("Phi", "#phi of recoTracks", hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
96 
97  hSimP_ = ibooker.book1D("SimP", "p of simTracks", hDim.nBinP, hDim.minP, hDim.maxP);
98  hSimPt_ = ibooker.book1D("SimPt", "p_{T} of simTracks", hDim.nBinPt, hDim.minPt, hDim.maxPt);
99  hSimEta_ = ibooker.book1D("SimEta", "#eta of simTracks", hDim.nBinEta, hDim.minEta, hDim.maxEta);
100  hSimPhi_ = ibooker.book1D("SimPhi", "#phi of simTracks", hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
101  hSimDxy_ = ibooker.book1D("SimDxy", "Dxy of simTracks", hDim.nBinDxy, hDim.minDxy, hDim.maxDxy);
102  hSimDz_ = ibooker.book1D("Dz", "Dz of simTracks", hDim.nBinDz, hDim.minDz, hDim.maxDz);
103 
104  //track multiplicities
105  hNSim_ = ibooker.book1D("NSim", "Number of particles per event", hDim.nTrks, -0.5, hDim.nTrks + 0.5);
106  hNMuon_ = ibooker.book1D("NMuon", "Number of muons per event", hDim.nTrks, -0.5, hDim.nTrks + 0.5);
107 
108  // - Misc. variables
109  hNTrks_ = ibooker.book1D("NTrks", "Number of reco tracks per event", hDim.nTrks, -0.5, hDim.nTrks + 0.5);
110  hNTrksEta_ = ibooker.book1D("NTrksEta", "Number of reco tracks vs #eta", hDim.nBinEta, hDim.minEta, hDim.maxEta);
111  hNTrksPt_ = ibooker.book1D("NTrksPt", "Number of reco tracks vs p_{T}", hDim.nBinPt, hDim.minPt, hDim.maxPt);
112 
113  hMisQPt_ = ibooker.book1D("MisQPt", "Charge mis-id vs Pt", hDim.nBinPt, hDim.minPt, hDim.maxPt);
114  hMisQEta_ = ibooker.book1D("MisQEta", "Charge mis-id vs Eta", hDim.nBinEta, hDim.minEta, hDim.maxEta);
115 
116  // - Resolutions
117  hErrP_ = ibooker.book1D("ErrP", "#Delta(p)/p", hDim.nBinErr, hDim.minErrP, hDim.maxErrP);
118  hErrPBarrel_ = ibooker.book1D("ErrP_barrel", "#Delta(p)/p", hDim.nBinErr, hDim.minErrP, hDim.maxErrP);
119  hErrPOverlap_ = ibooker.book1D("ErrP_overlap", "#Delta(p)/p", hDim.nBinErr, hDim.minErrP, hDim.maxErrP);
120  hErrPEndcap_ = ibooker.book1D("ErrP_endcap", "#Delta(p)/p", hDim.nBinErr, hDim.minErrP, hDim.maxErrP);
121  hErrPt_ = ibooker.book1D("ErrPt", "#Delta(p_{T})/p_{T}", hDim.nBinErr, hDim.minErrPt, hDim.maxErrPt);
122  hErrPtBarrel_ = ibooker.book1D("ErrPt_barrel", "#Delta(p_{T})/p_{T}", hDim.nBinErr, hDim.minErrPt, hDim.maxErrPt);
123  hErrPtOverlap_ = ibooker.book1D("ErrPt_overlap", "#Delta(p_{T})/p_{T}", hDim.nBinErr, hDim.minErrPt, hDim.maxErrPt);
124  hErrPtEndcap_ = ibooker.book1D("ErrPt_endcap", "#Delta(p_{T})/p_{T}", hDim.nBinErr, hDim.minErrPt, hDim.maxErrPt);
125  hErrEta_ = ibooker.book1D("ErrEta", "#sigma(#eta))", hDim.nBinErr, hDim.minErrEta, hDim.maxErrEta);
126  hErrPhi_ = ibooker.book1D("ErrPhi", "#sigma(#phi)", hDim.nBinErr, hDim.minErrPhi, hDim.maxErrPhi);
127  hErrDxy_ = ibooker.book1D("ErrDxy", "#sigma(d_{xy})", hDim.nBinErr, hDim.minErrDxy, hDim.maxErrDxy);
128  hErrDz_ = ibooker.book1D("ErrDz", "#sigma(d_{z})", hDim.nBinErr, hDim.minErrDz, hDim.maxErrDz);
129 
130  //PF-RECO comparisons
131  if (usePFMuon_) {
132  hErrPt_PF_ = ibooker.book1D("ErrPt_PF", "#Delta(p_{T})|_{PF}/p_{T}", hDim.nBinErr, hDim.minErrPt, hDim.maxErrPt);
133  hErrQPt_PF_ =
134  ibooker.book1D("ErrQPt_PF", "#Delta(q/p_{T})|_{PF}/(q/p_{T})", hDim.nBinErr, hDim.minErrQPt, hDim.maxErrQPt);
135 
136  hPFMomAssCorrectness =
137  ibooker.book1D("hPFMomAssCorrectness", "Corrected momentum assignement PF/RECO", 2, 0.5, 2.5);
138  hPt_vs_PFMomAssCorrectness = ibooker.book2D("hPt_vs_PFMomAssCorrectness",
139  "Corrected momentum assignement PF/RECO",
140  hDim.nBinPt,
141  hDim.minPt,
142  hDim.maxP,
143  2,
144  0.5,
145  2.5);
146 
147  hdPt_vs_Pt_ = ibooker.book2D("dPt_vs_Pt",
148  "#Delta(p_{T}) vs p_{T}",
149  hDim.nBinPt,
150  hDim.minPt,
151  hDim.maxPt,
152  hDim.nBinErr,
153  hDim.minErrPt,
154  hDim.maxErrPt);
155  hdPt_vs_Eta_ = ibooker.book2D("dPt_vs_Eta",
156  "#Delta(p_{T}) vs #eta",
157  hDim.nBinEta,
158  hDim.minEta,
159  hDim.maxEta,
160  hDim.nBinErr,
161  hDim.minErrPt,
162  hDim.maxErrPt);
163  }
164 
165  // -- Resolutions vs Eta
166  hErrP_vs_Eta_ = ibooker.book2D("ErrP_vs_Eta",
167  "#Delta(p)/p vs #eta",
168  hDim.nBinEta,
169  hDim.minEta,
170  hDim.maxEta,
171  hDim.nBinErr,
172  hDim.minErrP,
173  hDim.maxErrP);
174  hErrPt_vs_Eta_ = ibooker.book2D("ErrPt_vs_Eta",
175  "#Delta(p_{T})/p_{T} vs #eta",
176  hDim.nBinEta,
177  hDim.minEta,
178  hDim.maxEta,
179  hDim.nBinErr,
180  hDim.minErrPt,
181  hDim.maxErrPt);
182  hErrQPt_vs_Eta_ = ibooker.book2D("ErrQPt_vs_Eta",
183  "#Delta(q/p_{T})/(q/p_{T}) vs #eta",
184  hDim.nBinEta,
185  hDim.minEta,
186  hDim.maxEta,
187  hDim.nBinErr,
188  hDim.minErrQPt,
189  hDim.maxErrQPt);
190  hErrEta_vs_Eta_ = ibooker.book2D("ErrEta_vs_Eta",
191  "#sigma(#eta) vs #eta",
192  hDim.nBinEta,
193  hDim.minEta,
194  hDim.maxEta,
195  hDim.nBinErr,
196  hDim.minErrEta,
197  hDim.maxErrEta);
198 
199  // -- Resolutions vs momentum
200  hErrP_vs_P_ = ibooker.book2D(
201  "ErrP_vs_P", "#Delta(p)/p vs p", hDim.nBinP, hDim.minP, hDim.maxP, hDim.nBinErr, hDim.minErrP, hDim.maxErrP);
202  hErrPt_vs_Pt_ = ibooker.book2D("ErrPt_vs_Pt",
203  "#Delta(p_{T})/p_{T} vs p_{T}",
204  hDim.nBinPt,
205  hDim.minPt,
206  hDim.maxPt,
207  hDim.nBinErr,
208  hDim.minErrPt,
209  hDim.maxErrPt);
210  hErrQPt_vs_Pt_ = ibooker.book2D("ErrQPt_vs_Pt",
211  "#Delta(q/p_{T})/(q/p_{T}) vs p_{T}",
212  hDim.nBinPt,
213  hDim.minPt,
214  hDim.maxPt,
215  hDim.nBinErr,
216  hDim.minErrQPt,
217  hDim.maxErrQPt);
218 
219  // - Pulls
220  hPullPt_ = ibooker.book1D("PullPt", "Pull(#p_{T})", hDim.nBinPull, -hDim.wPull, hDim.wPull);
221  hPullEta_ = ibooker.book1D("PullEta", "Pull(#eta)", hDim.nBinPull, -hDim.wPull, hDim.wPull);
222  hPullPhi_ = ibooker.book1D("PullPhi", "Pull(#phi)", hDim.nBinPull, -hDim.wPull, hDim.wPull);
223  hPullQPt_ = ibooker.book1D("PullQPt", "Pull(q/p_{T})", hDim.nBinPull, -hDim.wPull, hDim.wPull);
224  hPullDxy_ = ibooker.book1D("PullDxy", "Pull(D_{xy})", hDim.nBinPull, -hDim.wPull, hDim.wPull);
225  hPullDz_ = ibooker.book1D("PullDz", "Pull(D_{z})", hDim.nBinPull, -hDim.wPull, hDim.wPull);
226 
227  // -- Pulls vs Eta
228  hPullPt_vs_Eta_ = ibooker.book2D("PullPt_vs_Eta",
229  "Pull(p_{T}) vs #eta",
230  hDim.nBinEta,
231  hDim.minEta,
232  hDim.maxEta,
233  hDim.nBinPull,
234  -hDim.wPull,
235  hDim.wPull);
236  hPullEta_vs_Eta_ = ibooker.book2D("PullEta_vs_Eta",
237  "Pull(#eta) vs #eta",
238  hDim.nBinEta,
239  hDim.minEta,
240  hDim.maxEta,
241  hDim.nBinPull,
242  -hDim.wPull,
243  hDim.wPull);
244  hPullPhi_vs_Eta_ = ibooker.book2D("PullPhi_vs_Eta",
245  "Pull(#phi) vs #eta",
246  hDim.nBinEta,
247  hDim.minEta,
248  hDim.maxEta,
249  hDim.nBinPull,
250  -hDim.wPull,
251  hDim.wPull);
252 
253  // -- Pulls vs Pt
254  hPullPt_vs_Pt_ = ibooker.book2D("PullPt_vs_Pt",
255  "Pull(p_{T}) vs p_{T}",
256  hDim.nBinPt,
257  hDim.minPt,
258  hDim.maxPt,
259  hDim.nBinPull,
260  -hDim.wPull,
261  hDim.wPull);
262  hPullEta_vs_Pt_ = ibooker.book2D("PullEta_vs_Pt",
263  "Pull(#eta) vs p_{T}",
264  hDim.nBinPt,
265  hDim.minPt,
266  hDim.maxPt,
267  hDim.nBinPull,
268  -hDim.wPull,
269  hDim.wPull);
270 
271  // -- Number of Hits
272  const int nHits = 100;
273  hNHits_ = ibooker.book1D("NHits", "Number of hits", nHits + 1, -0.5, nHits + 0.5);
274  hNHits_vs_Pt_ = ibooker.book2D("NHits_vs_Pt",
275  "Number of hits vs p_{T}",
276  hDim.nBinPt,
277  hDim.minPt,
278  hDim.maxPt,
279  nHits / 4 + 1,
280  -0.25,
281  nHits + 0.25);
282  hNHits_vs_Eta_ = ibooker.book2D("NHits_vs_Eta",
283  "Number of hits vs #eta",
284  hDim.nBinEta,
285  hDim.minEta,
286  hDim.maxEta,
287  nHits / 4 + 1,
288  -0.25,
289  nHits + 0.25);
290  hNSimHits_ = ibooker.book1D("NSimHits", "Number of simHits", nHits + 1, -0.5, nHits + 0.5);
291 
292  const int nLostHits = 5;
293  hNLostHits_ = ibooker.book1D("NLostHits", "Number of Lost hits", nLostHits + 1, -0.5, nLostHits + 0.5);
294  hNLostHits_vs_Pt_ = ibooker.book2D("NLostHits_vs_Pt",
295  "Number of lost Hits vs p_{T}",
296  hDim.nBinPt,
297  hDim.minPt,
298  hDim.maxPt,
299  nLostHits + 1,
300  -0.5,
301  nLostHits + 0.5);
302  hNLostHits_vs_Eta_ = ibooker.book2D("NLostHits_vs_Eta",
303  "Number of lost Hits vs #eta",
304  hDim.nBinEta,
305  hDim.minEta,
306  hDim.maxEta,
307  nLostHits + 1,
308  -0.5,
309  nLostHits + 0.5);
310 
311  const int nTrackerHits = 40;
312  hNTrackerHits_ =
313  ibooker.book1D("NTrackerHits", "Number of valid tracker hits", nTrackerHits + 1, -0.5, nTrackerHits + 0.5);
314  hNTrackerHits_vs_Pt_ = ibooker.book2D("NTrackerHits_vs_Pt",
315  "Number of valid traker hits vs p_{T}",
316  hDim.nBinPt,
317  hDim.minPt,
318  hDim.maxPt,
319  nTrackerHits / 4 + 1,
320  -0.25,
321  nTrackerHits + 0.25);
322  hNTrackerHits_vs_Eta_ = ibooker.book2D("NTrackerHits_vs_Eta",
323  "Number of valid tracker hits vs #eta",
324  hDim.nBinEta,
325  hDim.minEta,
326  hDim.maxEta,
327  nTrackerHits / 4 + 1,
328  -0.25,
329  nTrackerHits + 0.25);
330 
331  const int nMuonHits = 60;
332  hNMuonHits_ = ibooker.book1D("NMuonHits", "Number of valid muon hits", nMuonHits + 1, -0.5, nMuonHits + 0.5);
333  hNMuonHits_vs_Pt_ = ibooker.book2D("NMuonHits_vs_Pt",
334  "Number of valid muon hits vs p_{T}",
335  hDim.nBinPt,
336  hDim.minPt,
337  hDim.maxPt,
338  nMuonHits / 4 + 1,
339  -0.25,
340  nMuonHits + 0.25);
341  hNMuonHits_vs_Eta_ = ibooker.book2D("NMuonHits_vs_Eta",
342  "Number of valid muon hits vs #eta",
343  hDim.nBinEta,
344  hDim.minEta,
345  hDim.maxEta,
346  nMuonHits / 4 + 1,
347  -0.25,
348  nMuonHits + 0.25);
349 
350  hNDof_ = ibooker.book1D("NDof", "Number of DoF", hDim.nDof + 1, -0.5, hDim.nDof + 0.5);
351  hChi2_ = ibooker.book1D("Chi2", "#Chi^{2}", hDim.nBinErr, 0, 200);
352  hChi2Norm_ = ibooker.book1D("Chi2Norm", "Normalized #Chi^{2}", hDim.nBinErr, 0, 50);
353  hChi2Prob_ = ibooker.book1D("Chi2Prob", "Prob(#Chi^{2})", hDim.nBinErr, 0, 1);
354 
355  hNDof_vs_Eta_ = ibooker.book2D("NDof_vs_Eta",
356  "Number of DoF vs #eta",
357  hDim.nBinEta,
358  hDim.minEta,
359  hDim.maxEta,
360  hDim.nDof + 1,
361  -0.5,
362  hDim.nDof + 0.5);
363  hChi2_vs_Eta_ = ibooker.book2D(
364  "Chi2_vs_Eta", "#Chi^{2} vs #eta", hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinErr, 0., 200.);
365  hChi2Norm_vs_Eta_ = ibooker.book2D(
366  "Chi2Norm_vs_Eta", "Normalized #Chi^{2} vs #eta", hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinErr, 0., 50.);
367  hChi2Prob_vs_Eta_ = ibooker.book2D(
368  "Chi2Prob_vs_Eta", "Prob(#Chi^{2}) vs #eta", hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinErr, 0., 1.);
369 
370  hNSimToReco_ =
371  ibooker.book1D("NSimToReco", "Number of associated reco tracks", hDim.nAssoc + 1, -0.5, hDim.nAssoc + 0.5);
372  hNRecoToSim_ =
373  ibooker.book1D("NRecoToSim", "Number of associated sim TP's", hDim.nAssoc + 1, -0.5, hDim.nAssoc + 0.5);
374  };
375 
376  //
377  //Fill hists booked, simRef and muonRef are associated by hits
378  //
379  void fill(const TrackingParticle* simRef, const Muon* muonRef) {
380  const double simP = simRef->p();
381  const double simPt = simRef->pt();
382  const double simEta = doAbsEta_ ? fabs(simRef->eta()) : simRef->eta();
383  const double simPhi = simRef->phi();
384  const double simQ = simRef->charge();
385  const double simQPt = simQ / simPt;
386 
387  GlobalPoint simVtx(simRef->vertex().x(), simRef->vertex().y(), simRef->vertex().z());
388  GlobalVector simMom(simRef->momentum().x(), simRef->momentum().y(), simRef->momentum().z());
389  const double simDxy = -simVtx.x() * sin(simPhi) + simVtx.y() * cos(simPhi);
390  const double simDz = simVtx.z() - (simVtx.x() * simMom.x() + simVtx.y() * simMom.y()) * simMom.z() / simMom.perp2();
391 
392  const double recoQ = muonRef->charge();
393  if (simQ * recoQ < 0) {
394  hMisQPt_->Fill(simPt);
395  hMisQEta_->Fill(simEta);
396  }
397 
398  double recoP, recoPt, recoEta, recoPhi, recoQPt;
399  if (usePFMuon_) {
400  // const double origRecoP = muonRef->p();
401  const double origRecoPt = muonRef->pt();
402  // const double origRecoEta = muonRef->eta();
403  // const double origRecoPhi = muonRef->phi();
404  const double origRecoQPt = recoQ / origRecoPt;
405  recoP = muonRef->pfP4().P();
406  recoPt = muonRef->pfP4().Pt();
407  recoEta = muonRef->pfP4().Eta();
408  recoPhi = muonRef->pfP4().Phi();
409  recoQPt = recoQ / recoPt;
410  hErrPt_PF_->Fill((recoPt - origRecoPt) / origRecoPt);
411  hErrQPt_PF_->Fill((recoQPt - origRecoQPt) / origRecoQPt);
412 
413  hdPt_vs_Eta_->Fill(recoEta, recoPt - origRecoPt);
414  hdPt_vs_Pt_->Fill(recoPt, recoPt - origRecoPt);
415 
416  int theCorrectPFAss = (fabs(recoPt - simPt) < fabs(origRecoPt - simPt)) ? 1 : 2;
417  hPFMomAssCorrectness->Fill(theCorrectPFAss);
418  hPt_vs_PFMomAssCorrectness->Fill(simPt, theCorrectPFAss);
419  }
420 
421  else {
422  recoP = muonRef->p();
423  recoPt = muonRef->pt();
424  recoEta = muonRef->eta();
425  recoPhi = muonRef->phi();
426  recoQPt = recoQ / recoPt;
427  }
428 
429  const double errP = (recoP - simP) / simP;
430  const double errPt = (recoPt - simPt) / simPt;
431  const double errEta = (recoEta - simEta) / simEta;
432  const double errPhi = (recoPhi - simPhi) / simPhi;
433  const double errQPt = (recoQPt - simQPt) / simQPt;
434 
435  hP_->Fill(simP);
436  hPt_->Fill(simPt);
437  hEta_->Fill(simEta);
438  hPhi_->Fill(simPhi);
439 
440  hErrP_->Fill(errP);
441  hErrPt_->Fill(errPt);
442  hErrEta_->Fill(errEta);
443  hErrPhi_->Fill(errPhi);
444 
445  if (fabs(simEta) > 0. && fabs(simEta) < 0.8) {
446  hErrPBarrel_->Fill(errP);
447  hErrPtBarrel_->Fill(errPt);
448  } else if (fabs(simEta) > 0.8 && fabs(simEta) < 1.2) {
449  hErrPOverlap_->Fill(errP);
450  hErrPtOverlap_->Fill(errPt);
451  } else if (fabs(simEta) > 1.2) {
452  hErrPEndcap_->Fill(errP);
453  hErrPtEndcap_->Fill(errPt);
454  }
455 
456  hErrP_vs_Eta_->Fill(simEta, errP);
457  hErrPt_vs_Eta_->Fill(simEta, errPt);
458  hErrQPt_vs_Eta_->Fill(simEta, errQPt);
459 
460  hErrP_vs_P_->Fill(simP, errP);
461  hErrPt_vs_Pt_->Fill(simPt, errPt);
462  hErrQPt_vs_Pt_->Fill(simQPt, errQPt);
463 
464  hErrEta_vs_Eta_->Fill(simEta, errEta);
465 
466  //access from track
467  reco::TrackRef recoRef = muonRef->track();
468  if (recoRef.isNonnull()) {
469  // Number of reco-hits
470  const int nRecoHits = recoRef->numberOfValidHits();
471  const int nLostHits = recoRef->numberOfLostHits();
472 
473  hNHits_->Fill(nRecoHits);
474  hNHits_vs_Pt_->Fill(simPt, nRecoHits);
475  hNHits_vs_Eta_->Fill(simEta, nRecoHits);
476 
477  hNLostHits_->Fill(nLostHits);
478  hNLostHits_vs_Pt_->Fill(simPt, nLostHits);
479  hNLostHits_vs_Eta_->Fill(simEta, nLostHits);
480 
481  const double recoNDof = recoRef->ndof();
482  const double recoChi2 = recoRef->chi2();
483  const double recoChi2Norm = recoRef->normalizedChi2();
484  const double recoChi2Prob = TMath::Prob(recoRef->chi2(), static_cast<int>(recoRef->ndof()));
485 
486  hNDof_->Fill(recoNDof);
487  hChi2_->Fill(recoChi2);
488  hChi2Norm_->Fill(recoChi2Norm);
489  hChi2Prob_->Fill(recoChi2Prob);
490 
491  hNDof_vs_Eta_->Fill(simEta, recoNDof);
492  hChi2_vs_Eta_->Fill(simEta, recoChi2);
493  hChi2Norm_vs_Eta_->Fill(simEta, recoChi2Norm);
494  hChi2Prob_vs_Eta_->Fill(simEta, recoChi2Prob);
495 
496  const double recoDxy = recoRef->dxy();
497  const double recoDz = recoRef->dz();
498 
499  const double errDxy = (recoDxy - simDxy) / simDxy;
500  const double errDz = (recoDz - simDz) / simDz;
501  hErrDxy_->Fill(errDxy);
502  hErrDz_->Fill(errDz);
503 
504  const double pullPt = (recoPt - simPt) / recoRef->ptError();
505  const double pullQPt = (recoQPt - simQPt) / recoRef->qoverpError();
506  const double pullEta = (recoEta - simEta) / recoRef->etaError();
507  const double pullPhi = (recoPhi - simPhi) / recoRef->phiError();
508  const double pullDxy = (recoDxy - simDxy) / recoRef->dxyError();
509  const double pullDz = (recoDz - simDz) / recoRef->dzError();
510 
511  hPullPt_->Fill(pullPt);
512  hPullEta_->Fill(pullEta);
513  hPullPhi_->Fill(pullPhi);
514  hPullQPt_->Fill(pullQPt);
515  hPullDxy_->Fill(pullDxy);
516  hPullDz_->Fill(pullDz);
517 
518  hPullPt_vs_Eta_->Fill(simEta, pullPt);
519  hPullPt_vs_Pt_->Fill(simPt, pullPt);
520 
521  hPullEta_vs_Eta_->Fill(simEta, pullEta);
522  hPullPhi_vs_Eta_->Fill(simEta, pullPhi);
523 
524  hPullEta_vs_Pt_->Fill(simPt, pullEta);
525  }
526  };
527 };
528 
529 //
530 //struct defininiong histograms
531 //
532 struct RecoMuonValidator::CommonME {
533  typedef MonitorElement* MEP;
534 
535  //diffs
536  MEP hTrkToGlbDiffNTrackerHits_, hStaToGlbDiffNMuonHits_;
537  MEP hTrkToGlbDiffNTrackerHitsEta_, hStaToGlbDiffNMuonHitsEta_;
538  MEP hTrkToGlbDiffNTrackerHitsPt_, hStaToGlbDiffNMuonHitsPt_;
539 
540  //global muon hit pattern
541  MEP hNInvalidHitsGTHitPattern_, hNInvalidHitsITHitPattern_, hNInvalidHitsOTHitPattern_;
542  MEP hNDeltaInvalidHitsHitPattern_;
543 
544  //muon based momentum assignment
545  MEP hMuonP_, hMuonPt_, hMuonEta_, hMuonPhi_;
546  //track based kinematics
547  MEP hMuonTrackP_, hMuonTrackPt_, hMuonTrackEta_, hMuonTrackPhi_, hMuonTrackDxy_, hMuonTrackDz_;
548  //histograms for fractions
549  MEP hMuonAllP_, hMuonAllPt_, hMuonAllEta_, hMuonAllPhi_;
550 };
551 
552 //
553 //Constructor
554 //
555 RecoMuonValidator::RecoMuonValidator(const edm::ParameterSet& pset)
556  : selector_(pset.getParameter<std::string>("selection")) {
557  // pset=ps;
558  verbose_ = pset.getUntrackedParameter<unsigned int>("verbose", 0);
559 
560  outputFileName_ = pset.getUntrackedParameter<string>("outputFileName", "");
561 
562  wantTightMuon_ = pset.getParameter<bool>("wantTightMuon");
563  beamspotLabel_ = pset.getParameter<edm::InputTag>("beamSpot");
564  primvertexLabel_ = pset.getParameter<edm::InputTag>("primaryVertex");
565  beamspotToken_ = consumes<reco::BeamSpot>(beamspotLabel_);
566  primvertexToken_ = consumes<reco::VertexCollection>(primvertexLabel_);
567 
568  // Set histogram dimensions from config
569 
570  hDim.nBinP = pset.getUntrackedParameter<unsigned int>("nBinP");
571  hDim.minP = pset.getUntrackedParameter<double>("minP");
572  hDim.maxP = pset.getUntrackedParameter<double>("maxP");
573 
574  hDim.nBinPt = pset.getUntrackedParameter<unsigned int>("nBinPt");
575  hDim.minPt = pset.getUntrackedParameter<double>("minPt");
576  hDim.maxPt = pset.getUntrackedParameter<double>("maxPt");
577 
578  doAbsEta_ = pset.getUntrackedParameter<bool>("doAbsEta");
579  hDim.doAbsEta = doAbsEta_;
580  hDim.nBinEta = pset.getUntrackedParameter<unsigned int>("nBinEta");
581  hDim.minEta = pset.getUntrackedParameter<double>("minEta");
582  hDim.maxEta = pset.getUntrackedParameter<double>("maxEta");
583 
584  hDim.nBinDxy = pset.getUntrackedParameter<unsigned int>("nBinDxy");
585  hDim.minDxy = pset.getUntrackedParameter<double>("minDxy");
586  hDim.maxDxy = pset.getUntrackedParameter<double>("maxDxy");
587 
588  hDim.nBinDz = pset.getUntrackedParameter<unsigned int>("nBinDz");
589  hDim.minDz = pset.getUntrackedParameter<double>("minDz");
590  hDim.maxDz = pset.getUntrackedParameter<double>("maxDz");
591 
592  hDim.nBinPhi = pset.getUntrackedParameter<unsigned int>("nBinPhi");
593  hDim.minPhi = pset.getUntrackedParameter<double>("minPhi", -TMath::Pi());
594  hDim.maxPhi = pset.getUntrackedParameter<double>("maxPhi", TMath::Pi());
595 
596  hDim.nBinErr = pset.getUntrackedParameter<unsigned int>("nBinErr");
597  hDim.nBinPull = pset.getUntrackedParameter<unsigned int>("nBinPull");
598 
599  hDim.wPull = pset.getUntrackedParameter<double>("wPull");
600 
601  hDim.minErrP = pset.getUntrackedParameter<double>("minErrP");
602  hDim.maxErrP = pset.getUntrackedParameter<double>("maxErrP");
603 
604  hDim.minErrPt = pset.getUntrackedParameter<double>("minErrPt");
605  hDim.maxErrPt = pset.getUntrackedParameter<double>("maxErrPt");
606 
607  hDim.minErrQPt = pset.getUntrackedParameter<double>("minErrQPt");
608  hDim.maxErrQPt = pset.getUntrackedParameter<double>("maxErrQPt");
609 
610  hDim.minErrEta = pset.getUntrackedParameter<double>("minErrEta");
611  hDim.maxErrEta = pset.getUntrackedParameter<double>("maxErrEta");
612 
613  hDim.minErrPhi = pset.getUntrackedParameter<double>("minErrPhi");
614  hDim.maxErrPhi = pset.getUntrackedParameter<double>("maxErrPhi");
615 
616  hDim.minErrDxy = pset.getUntrackedParameter<double>("minErrDxy");
617  hDim.maxErrDxy = pset.getUntrackedParameter<double>("maxErrDxy");
618 
619  hDim.minErrDz = pset.getUntrackedParameter<double>("minErrDz");
620  hDim.maxErrDz = pset.getUntrackedParameter<double>("maxErrDz");
621 
622  hDim.nTrks = pset.getUntrackedParameter<unsigned int>("nTrks");
623  hDim.nAssoc = pset.getUntrackedParameter<unsigned int>("nAssoc");
624  hDim.nDof = pset.getUntrackedParameter<unsigned int>("nDof", 55);
625 
626  // Labels for simulation and reconstruction tracks
627  simLabel_ = pset.getParameter<InputTag>("simLabel");
628  muonLabel_ = pset.getParameter<InputTag>("muonLabel");
629  simToken_ = consumes<TrackingParticleCollection>(simLabel_);
630  muonToken_ = consumes<edm::View<reco::Muon> >(muonLabel_);
631 
632  // Labels for sim-reco association
633  doAssoc_ = pset.getUntrackedParameter<bool>("doAssoc", true);
634  muAssocLabel_ = pset.getParameter<InputTag>("muAssocLabel");
635  if (doAssoc_) {
636  muAssocToken_ = consumes<reco::MuonToTrackingParticleAssociator>(muAssocLabel_);
637  }
638 
639  // Different momentum assignment and additional histos in case of PF muons
640  usePFMuon_ = pset.getUntrackedParameter<bool>("usePFMuon");
641  hDim.usePFMuon = usePFMuon_;
642 
643  //type of track
644  std::string trackType = pset.getParameter<std::string>("trackType");
645  if (trackType == "inner")
646  trackType_ = reco::InnerTk;
647  else if (trackType == "outer")
648  trackType_ = reco::OuterTk;
649  else if (trackType == "global")
650  trackType_ = reco::GlobalTk;
651  else if (trackType == "segments")
652  trackType_ = reco::Segments;
653  else
654  throw cms::Exception("Configuration") << "Track type '" << trackType << "' not supported.\n";
655 
656  // seedPropagatorName_ = pset.getParameter<string>("SeedPropagator");
657 
658  ParameterSet tpset = pset.getParameter<ParameterSet>("tpSelector");
659  tpSelector_ = TrackingParticleSelector(tpset.getParameter<double>("ptMin"),
660  tpset.getParameter<double>("ptMax"),
661  tpset.getParameter<double>("minRapidity"),
662  tpset.getParameter<double>("maxRapidity"),
663  tpset.getParameter<double>("tip"),
664  tpset.getParameter<double>("lip"),
665  tpset.getParameter<int>("minHit"),
666  tpset.getParameter<bool>("signalOnly"),
667  tpset.getParameter<bool>("intimeOnly"),
668  tpset.getParameter<bool>("chargedOnly"),
669  tpset.getParameter<bool>("stableOnly"),
670  tpset.getParameter<std::vector<int> >("pdgId"));
671 
672  // the service parameters
673  ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
674  theMuonService = new MuonServiceProxy(serviceParameters);
675 
676  // retrieve the instance of DQMService
677  dbe_ = Service<DQMStore>().operator->();
678  subsystemname_ = pset.getUntrackedParameter<std::string>("subSystemFolder", "YourSubsystem");
679 
680  subDir_ = pset.getUntrackedParameter<string>("subDir");
681  if (subDir_.empty())
682  subDir_ = "RecoMuonV";
683  if (subDir_[subDir_.size() - 1] == '/')
684  subDir_.erase(subDir_.size() - 1);
685 }
686 
687 void RecoMuonValidator::bookHistograms(DQMStore::IBooker& ibooker,
688  edm::Run const& iRun,
689  edm::EventSetup const& /* iSetup */) {
690  // book histograms
691  ibooker.cd();
692 
693  ibooker.setCurrentFolder(subDir_);
694 
695  commonME_ = new CommonME;
696  muonME_ = new MuonME;
697 
698  //commonME
699  const int nHits = 100;
700 
701  // - diffs
702  commonME_->hTrkToGlbDiffNTrackerHits_ = ibooker.book1D("TrkGlbDiffNTrackerHits",
703  "Difference of number of tracker hits (tkMuon - globalMuon)",
704  2 * nHits + 1,
705  -nHits - 0.5,
706  nHits + 0.5);
707  commonME_->hStaToGlbDiffNMuonHits_ = ibooker.book1D("StaGlbDiffNMuonHits",
708  "Difference of number of muon hits (staMuon - globalMuon)",
709  2 * nHits + 1,
710  -nHits - 0.5,
711  nHits + 0.5);
712 
713  commonME_->hTrkToGlbDiffNTrackerHitsEta_ =
714  ibooker.book2D("TrkGlbDiffNTrackerHitsEta",
715  "Difference of number of tracker hits (tkMuon - globalMuon)",
716  hDim.nBinEta,
717  hDim.minEta,
718  hDim.maxEta,
719  2 * nHits + 1,
720  -nHits - 0.5,
721  nHits + 0.5);
722  commonME_->hStaToGlbDiffNMuonHitsEta_ = ibooker.book2D("StaGlbDiffNMuonHitsEta",
723  "Difference of number of muon hits (staMuon - globalMuon)",
724  hDim.nBinEta,
725  hDim.minEta,
726  hDim.maxEta,
727  2 * nHits + 1,
728  -nHits - 0.5,
729  nHits + 0.5);
730 
731  commonME_->hTrkToGlbDiffNTrackerHitsPt_ = ibooker.book2D("TrkGlbDiffNTrackerHitsPt",
732  "Difference of number of tracker hits (tkMuon - globalMuon)",
733  hDim.nBinPt,
734  hDim.minPt,
735  hDim.maxPt,
736  2 * nHits + 1,
737  -nHits - 0.5,
738  nHits + 0.5);
739  commonME_->hStaToGlbDiffNMuonHitsPt_ = ibooker.book2D("StaGlbDiffNMuonHitsPt",
740  "Difference of number of muon hits (staMuon - globalMuon)",
741  hDim.nBinPt,
742  hDim.minPt,
743  hDim.maxPt,
744  2 * nHits + 1,
745  -nHits - 0.5,
746  nHits + 0.5);
747 
748  // -global muon hit pattern
749  commonME_->hNInvalidHitsGTHitPattern_ = ibooker.book1D(
750  "NInvalidHitsGTHitPattern", "Number of invalid hits on a global track", nHits + 1, -0.5, nHits + 0.5);
751  commonME_->hNInvalidHitsITHitPattern_ = ibooker.book1D(
752  "NInvalidHitsITHitPattern", "Number of invalid hits on an inner track", nHits + 1, -0.5, nHits + 0.5);
753  commonME_->hNInvalidHitsOTHitPattern_ = ibooker.book1D(
754  "NInvalidHitsOTHitPattern", "Number of invalid hits on an outer track", nHits + 1, -0.5, nHits + 0.5);
755  commonME_->hNDeltaInvalidHitsHitPattern_ =
756  ibooker.book1D("hNDeltaInvalidHitsHitPattern",
757  "The discrepancy for Number of invalid hits on an global track and inner and outer tracks",
758  2 * nHits + 1,
759  -nHits - 0.5,
760  nHits + 0.5);
761 
762  //muon based kinematics
763  commonME_->hMuonP_ = ibooker.book1D("PMuon", "p of muon", hDim.nBinP, hDim.minP, hDim.maxP);
764  commonME_->hMuonPt_ = ibooker.book1D("PtMuon", "p_{T} of muon", hDim.nBinPt, hDim.minPt, hDim.maxPt);
765  commonME_->hMuonEta_ = ibooker.book1D("EtaMuon", "#eta of muon", hDim.nBinEta, hDim.minEta, hDim.maxEta);
766  commonME_->hMuonPhi_ = ibooker.book1D("PhiMuon", "#phi of muon", hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
767  //track based kinematics
768  commonME_->hMuonTrackP_ = ibooker.book1D("PMuonTrack", "p of reco muon track", hDim.nBinP, hDim.minP, hDim.maxP);
769  commonME_->hMuonTrackPt_ =
770  ibooker.book1D("PtMuonTrack", "p_{T} of reco muon track", hDim.nBinPt, hDim.minPt, hDim.maxPt);
771  commonME_->hMuonTrackEta_ =
772  ibooker.book1D("EtaMuonTrack", "#eta of reco muon track", hDim.nBinEta, hDim.minEta, hDim.maxEta);
773  commonME_->hMuonTrackPhi_ =
774  ibooker.book1D("PhiMuonTrack", "#phi of reco muon track", hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
775  commonME_->hMuonTrackDxy_ =
776  ibooker.book1D("DxyMuonTrack", "Dxy of reco muon track", hDim.nBinDxy, hDim.minDxy, hDim.maxDxy);
777  commonME_->hMuonTrackDz_ =
778  ibooker.book1D("DzMuonTrack", "Dz of reco muon track", hDim.nBinDz, hDim.minDz, hDim.maxDz);
779 
780  //histograms for fractions
781  commonME_->hMuonAllP_ = ibooker.book1D("PMuonAll", "p of muons of all types", hDim.nBinP, hDim.minP, hDim.maxP);
782  commonME_->hMuonAllPt_ =
783  ibooker.book1D("PtMuonAll", "p_{T} of muon of all types", hDim.nBinPt, hDim.minPt, hDim.maxPt);
784  commonME_->hMuonAllEta_ =
785  ibooker.book1D("EtaMuonAll", "#eta of muon of all types", hDim.nBinEta, hDim.minEta, hDim.maxEta);
786  commonME_->hMuonAllPhi_ =
787  ibooker.book1D("PhiMuonAll", "#phi of muon of all types", hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
788 
789  muonME_->bookHistos(ibooker, subDir_, hDim);
790 }
791 
792 //
793 //Destructor
794 //
795 RecoMuonValidator::~RecoMuonValidator() {
796  if (theMuonService)
797  delete theMuonService;
798 }
799 
800 //
801 //Begin run
802 //
803 
804 void RecoMuonValidator::dqmBeginRun(const edm::Run&, const EventSetup& eventSetup) {
805  if (theMuonService)
806  theMuonService->update(eventSetup);
807 }
808 
809 //
810 //End run
811 //
812 void RecoMuonValidator::dqmEndRun(edm::Run const&, edm::EventSetup const&) {
813  if (dbe_ && !outputFileName_.empty())
814  dbe_->save(outputFileName_);
815 }
816 
817 //
818 //Analyze
819 //
820 void RecoMuonValidator::analyze(const Event& event, const EventSetup& eventSetup) {
821  // Look for the Primary Vertex (and use the BeamSpot instead, if you can't find it):
822  reco::Vertex::Point posVtx;
823  reco::Vertex::Error errVtx;
824  edm::Handle<reco::VertexCollection> recVtxs;
825  event.getByToken(primvertexToken_, recVtxs);
826  unsigned int theIndexOfThePrimaryVertex = 999.;
827  for (unsigned int ind = 0; ind < recVtxs->size(); ++ind) {
828  if ((*recVtxs)[ind].isValid() && !((*recVtxs)[ind].isFake())) {
829  theIndexOfThePrimaryVertex = ind;
830  break;
831  }
832  }
833  if (theIndexOfThePrimaryVertex < 100) {
834  posVtx = ((*recVtxs)[theIndexOfThePrimaryVertex]).position();
835  errVtx = ((*recVtxs)[theIndexOfThePrimaryVertex]).error();
836  } else {
837  LogInfo("RecoMuonValidator") << "reco::PrimaryVertex not found, use BeamSpot position instead\n";
838  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
839  event.getByToken(beamspotToken_, recoBeamSpotHandle);
840  reco::BeamSpot bs = *recoBeamSpotHandle;
841  posVtx = bs.position();
842  errVtx(0, 0) = bs.BeamWidthX();
843  errVtx(1, 1) = bs.BeamWidthY();
844  errVtx(2, 2) = bs.sigmaZ();
845  }
846  const reco::Vertex thePrimaryVertex(posVtx, errVtx);
847 
848  // Get TrackingParticles
849  Handle<TrackingParticleCollection> simHandle;
850  event.getByToken(simToken_, simHandle);
851  const TrackingParticleCollection simColl = *(simHandle.product());
852 
853  // Get Muons
854  Handle<edm::View<Muon> > muonHandle;
855  event.getByToken(muonToken_, muonHandle);
856  View<Muon> muonColl = *(muonHandle.product());
857 
858  reco::MuonToTrackingParticleAssociator const* assoByHits = nullptr;
859  if (doAssoc_) {
860  edm::Handle<reco::MuonToTrackingParticleAssociator> associatorBase;
861  event.getByToken(muAssocToken_, associatorBase);
862  assoByHits = associatorBase.product();
863  }
864 
865  const TrackingParticleCollection::size_type nSim = simColl.size();
866 
867  edm::RefToBaseVector<reco::Muon> Muons;
868  for (size_t i = 0; i < muonHandle->size(); ++i) {
869  Muons.push_back(muonHandle->refAt(i));
870  }
871 
872  edm::RefVector<TrackingParticleCollection> allTPs;
873  for (size_t i = 0; i < nSim; ++i) {
874  allTPs.push_back(TrackingParticleRef(simHandle, i));
875  }
876 
877  muonME_->hNSim_->Fill(nSim);
878  muonME_->hNMuon_->Fill(muonColl.size());
879 
880  reco::MuonToSimCollection muonToSimColl;
881  reco::SimToMuonCollection simToMuonColl;
882 
883  if (doAssoc_) {
884  assoByHits->associateMuons(muonToSimColl, simToMuonColl, Muons, trackType_, allTPs);
885  } else {
886  /*
887  // SimToMuon associations
888  Handle<reco::RecoToSimCollection> simToTrkMuHandle;
889  event.getByLabel(trkMuAssocLabel_, simToTrkMuHandle);
890  trkSimRecColl = *(simToTrkMuHandle.product());
891 
892  Handle<reco::RecoToSimCollection> simToStaMuHandle;
893  event.getByLabel(staMuAssocLabel_, simToStaMuHandle);
894  staSimRecColl = *(simToStaMuHandle.product());
895 
896  Handle<reco::RecoToSimCollection> simToGlbMuHandle;
897  event.getByLabel(glbMuAssocLabel_, simToGlbMuHandle);
898  glbSimRecColl = *(simToGlbMuHandle.product());
899 
900  // MuonToSim associations
901  Handle<reco::SimToRecoCollection> trkMuToSimHandle;
902  event.getByLabel(trkMuAssocLabel_, trkMuToSimHandle);
903  trkRecSimColl = *(trkMuToSimHandle.product());
904 
905  Handle<reco::SimToRecoCollection> staMuToSimHandle;
906  event.getByLabel(staMuAssocLabel_, staMuToSimHandle);
907  staRecSimColl = *(staMuToSimHandle.product());
908 
909  Handle<reco::SimToRecoCollection> glbMuToSimHandle;
910  event.getByLabel(glbMuAssocLabel_, glbMuToSimHandle);
911  glbRecSimColl = *(glbMuToSimHandle.product());
912 */
913  }
914 
915  int glbNTrackerHits = 0;
916  int trkNTrackerHits = 0;
917  int glbNMuonHits = 0;
918  int staNMuonHits = 0;
919  int NTrackerHits = 0;
920  int NMuonHits = 0;
921 
922  // Analyzer reco::Muon
923  for (View<Muon>::const_iterator iMuon = muonColl.begin(); iMuon != muonColl.end(); ++iMuon) {
924  double muonP, muonPt, muonEta, muonPhi;
925  if (usePFMuon_) {
926  muonP = iMuon->pfP4().P();
927  muonPt = iMuon->pfP4().Pt();
928  muonEta = iMuon->pfP4().Eta();
929  muonPhi = iMuon->pfP4().Phi();
930  } else {
931  muonP = iMuon->p();
932  muonPt = iMuon->pt();
933  muonEta = iMuon->eta();
934  muonPhi = iMuon->phi();
935  }
936 
937  //histograms for fractions
938  commonME_->hMuonAllP_->Fill(muonP);
939  commonME_->hMuonAllPt_->Fill(muonPt);
940  commonME_->hMuonAllEta_->Fill(muonEta);
941  commonME_->hMuonAllPhi_->Fill(muonPhi);
942 
943  if (!selector_(*iMuon))
944  continue;
945  if (wantTightMuon_) {
946  if (!muon::isTightMuon(*iMuon, thePrimaryVertex))
947  continue;
948  }
949 
950  TrackRef Track = iMuon->track();
951 
952  if (Track.isNonnull()) {
953  commonME_->hMuonTrackP_->Fill(Track->p());
954  commonME_->hMuonTrackPt_->Fill(Track->pt());
955  commonME_->hMuonTrackEta_->Fill(Track->eta());
956  commonME_->hMuonTrackPhi_->Fill(Track->phi());
957 
958  //ip histograms
959  commonME_->hMuonTrackDxy_->Fill(Track->dxy());
960  commonME_->hMuonTrackDz_->Fill(Track->dz());
961  }
962 
963  if (iMuon->isGlobalMuon()) {
964  Track = iMuon->combinedMuon();
965  glbNTrackerHits = countTrackerHits(*Track);
966  glbNMuonHits = countMuonHits(*Track);
967  } else if (iMuon->isTrackerMuon()) {
968  Track = iMuon->track();
969  trkNTrackerHits = countTrackerHits(*Track);
970  } else {
971  Track = iMuon->standAloneMuon();
972  }
973 
974  NTrackerHits = countTrackerHits(*Track);
975  muonME_->hNTrackerHits_->Fill(NTrackerHits);
976  muonME_->hNTrackerHits_vs_Pt_->Fill(Track->pt(), NTrackerHits);
977  muonME_->hNTrackerHits_vs_Eta_->Fill(Track->eta(), NTrackerHits);
978 
979  NMuonHits = countMuonHits(*Track);
980  muonME_->hNMuonHits_->Fill(NMuonHits);
981  muonME_->hNMuonHits_vs_Pt_->Fill(Track->pt(), NMuonHits);
982  muonME_->hNMuonHits_vs_Eta_->Fill(Track->eta(), NMuonHits);
983 
984  //list of histos for each type
985 
986  // muonME_->hNTrks_->Fill();
987  muonME_->hNTrksEta_->Fill(Track->eta());
988  muonME_->hNTrksPt_->Fill(Track->pt());
989 
990  commonME_->hMuonP_->Fill(muonP);
991  commonME_->hMuonPt_->Fill(muonPt);
992  commonME_->hMuonEta_->Fill(muonEta);
993  commonME_->hMuonPhi_->Fill(muonPhi);
994 
995  if (iMuon->isGlobalMuon()) {
996  double gtHitPat = iMuon->globalTrack()->hitPattern().numberOfAllHits(HitPattern::TRACK_HITS) -
997  iMuon->globalTrack()->hitPattern().numberOfValidHits();
998 
999  double itHitPat = iMuon->innerTrack()->hitPattern().numberOfAllHits(HitPattern::TRACK_HITS) -
1000  iMuon->innerTrack()->hitPattern().numberOfValidHits();
1001 
1002  double otHitPat = iMuon->outerTrack()->hitPattern().numberOfAllHits(HitPattern::TRACK_HITS) -
1003  iMuon->outerTrack()->hitPattern().numberOfValidHits();
1004 
1005  commonME_->hNInvalidHitsGTHitPattern_->Fill(gtHitPat);
1006  commonME_->hNInvalidHitsITHitPattern_->Fill(itHitPat);
1007  commonME_->hNInvalidHitsOTHitPattern_->Fill(otHitPat);
1008  commonME_->hNDeltaInvalidHitsHitPattern_->Fill(gtHitPat - itHitPat - otHitPat);
1009 
1010  //must be global and standalone
1011  if (iMuon->isStandAloneMuon()) {
1012  commonME_->hStaToGlbDiffNMuonHitsEta_->Fill(Track->eta(), staNMuonHits - glbNMuonHits);
1013  commonME_->hStaToGlbDiffNMuonHitsPt_->Fill(Track->pt(), staNMuonHits - glbNMuonHits);
1014  commonME_->hStaToGlbDiffNMuonHits_->Fill(staNMuonHits - glbNMuonHits);
1015  }
1016 
1017  //must be global and tracker
1018  if (iMuon->isTrackerMuon()) {
1019  commonME_->hTrkToGlbDiffNTrackerHitsEta_->Fill(Track->eta(), trkNTrackerHits - glbNTrackerHits);
1020  commonME_->hTrkToGlbDiffNTrackerHitsPt_->Fill(Track->pt(), trkNTrackerHits - glbNTrackerHits);
1021  commonME_->hTrkToGlbDiffNTrackerHits_->Fill(trkNTrackerHits - glbNTrackerHits);
1022  }
1023  }
1024 
1025  } //end of reco muon loop
1026 
1027  // Associate by hits
1028  for (TrackingParticleCollection::size_type i = 0; i < nSim; i++) {
1029  TrackingParticleRef simRef(simHandle, i);
1030  const TrackingParticle* simTP = simRef.get();
1031  if (!tpSelector_(*simTP))
1032  continue;
1033 
1034  //denominators for efficiency plots
1035  const double simP = simRef->p();
1036  const double simPt = simRef->pt();
1037  const double simEta = doAbsEta_ ? fabs(simRef->eta()) : simRef->eta();
1038  const double simPhi = simRef->phi();
1039 
1040  GlobalPoint simVtx(simRef->vertex().x(), simRef->vertex().y(), simRef->vertex().z());
1041  GlobalVector simMom(simRef->momentum().x(), simRef->momentum().y(), simRef->momentum().z());
1042  const double simDxy = -simVtx.x() * sin(simPhi) + simVtx.y() * cos(simPhi);
1043  const double simDz = simVtx.z() - (simVtx.x() * simMom.x() + simVtx.y() * simMom.y()) * simMom.z() / simMom.perp2();
1044 
1045  const unsigned int nSimHits = simRef->numberOfHits();
1046 
1047  muonME_->hSimP_->Fill(simP);
1048  muonME_->hSimPt_->Fill(simPt);
1049  muonME_->hSimEta_->Fill(simEta);
1050  muonME_->hSimPhi_->Fill(simPhi);
1051  muonME_->hSimDxy_->Fill(simDxy);
1052  muonME_->hSimDz_->Fill(simDz);
1053  muonME_->hNSimHits_->Fill(nSimHits);
1054 
1055  // Get sim-reco association for a simRef
1056  vector<pair<RefToBase<Muon>, double> > MuRefV;
1057  if (simToMuonColl.find(simRef) != simToMuonColl.end()) {
1058  MuRefV = simToMuonColl[simRef];
1059 
1060  if (!MuRefV.empty()) {
1061  muonME_->hNSimToReco_->Fill(MuRefV.size());
1062  const Muon* Mu = MuRefV.begin()->first.get();
1063  if (!selector_(*Mu))
1064  continue;
1065  if (wantTightMuon_) {
1066  if (!muon::isTightMuon(*Mu, thePrimaryVertex))
1067  continue;
1068  }
1069 
1070  muonME_->fill(&*simTP, Mu);
1071  }
1072  }
1073  }
1074 }
1075 
1076 int RecoMuonValidator::countMuonHits(const reco::Track& track) const {
1077  TransientTrackingRecHit::ConstRecHitContainer result;
1078 
1079  int count = 0;
1080 
1081  for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
1082  if ((*hit)->isValid()) {
1083  DetId recoid = (*hit)->geographicalId();
1084  if (recoid.det() == DetId::Muon)
1085  count++;
1086  }
1087  }
1088  return count;
1089 }
1090 
1091 int RecoMuonValidator::countTrackerHits(const reco::Track& track) const {
1092  TransientTrackingRecHit::ConstRecHitContainer result;
1093 
1094  int count = 0;
1095 
1096  for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
1097  if ((*hit)->isValid()) {
1098  DetId recoid = (*hit)->geographicalId();
1099  if (recoid.det() == DetId::Tracker)
1100  count++;
1101  }
1102  }
1103  return count;
1104 }
1105 /* vim:set ts=2 sts=2 sw=2 expandtab: */
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
dqm::dqmstoreimpl::DQMStore::IBooker::book1D
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX)
Definition: DQMStore.cc:239
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
RecoMuonValidator::muonToken_
edm::EDGetTokenT< edm::View< reco::Muon > > muonToken_
Definition: RecoMuonValidator.h:51
RecoMuonValidator::simToken_
edm::EDGetTokenT< TrackingParticleCollection > simToken_
Definition: RecoMuonValidator.h:50
reco::MuonToSimCollection
std::map< edm::RefToBase< reco::Muon >, std::vector< std::pair< TrackingParticleRef, double > >, RefToBaseSort > MuonToSimCollection
Definition: MuonTrackType.h:37
reco::Muon::track
TrackRef track() const override
reference to a Track
Definition: Muon.h:46
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
mps_fire.i
i
Definition: mps_fire.py:341
reco::LeafCandidate::eta
double eta() const final
momentum pseudorapidity
Definition: LeafCandidate.h:152
mps_fire.result
result
Definition: mps_fire.py:294
TrackingParticleCollection
std::vector< TrackingParticle > TrackingParticleCollection
Definition: HiEvtPlaneFlatProducer.cc:110
TSOS
TrajectoryStateOnSurface TSOS
Definition: RecoMuonValidator.cc:22
TrackingParticle::momentum
Vector momentum() const
spatial momentum vector
Definition: TrackingParticle.h:106
HLTObjectMonitor_cfi.muonEta
muonEta
Definition: HLTObjectMonitor_cfi.py:92
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
dqm::dqmstoreimpl::DQMStore::IBooker::setCurrentFolder
void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:418
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
RecoMuonValidator::simLabel_
edm::InputTag simLabel_
Definition: RecoMuonValidator.h:47
RecoMuonValidator::selector_
StringCutObjectSelector< reco::Muon > selector_
Definition: RecoMuonValidator.h:130
reco::Vertex::Error
math::Error< dimension >::type Error
covariance error matrix (3x3)
Definition: Vertex.h:43
TrackingParticle::pt
double pt() const
Transverse momentum. Note this is taken from the first SimTrack only.
Definition: TrackingParticle.h:142
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
edm::View::size
size_type size() const
reco::LeafCandidate::pt
double pt() const final
transverse momentum
Definition: LeafCandidate.h:146
HLTObjectMonitor_cfi.muonPhi
muonPhi
Definition: HLTObjectMonitor_cfi.py:102
reco::LeafCandidate::charge
int charge() const final
electric charge
Definition: LeafCandidate.h:106
RecoMuonValidator::trackType_
reco::MuonTrackType trackType_
Definition: RecoMuonValidator.h:76
RecoMuonValidator::primvertexToken_
edm::EDGetTokenT< reco::VertexCollection > primvertexToken_
Definition: RecoMuonValidator.h:59
RecoMuonValidator::countTrackerHits
virtual int countTrackerHits(const reco::Track &track) const
Definition: RecoMuonValidator.cc:1091
reco::Muon::pfP4
reco::Candidate::LorentzVector pfP4() const
Definition: Muon.h:97
trigger::size_type
uint16_t size_type
Definition: TriggerTypeDefs.h:18
RecoMuonValidator::bookHistograms
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
Definition: RecoMuonValidator.cc:687
RecoMuonValidator::dqmEndRun
void dqmEndRun(edm::Run const &, edm::EventSetup const &) override
Definition: RecoMuonValidator.cc:812
reco::SimToMuonCollection
std::map< TrackingParticleRef, std::vector< std::pair< edm::RefToBase< reco::Muon >, double > > > SimToMuonCollection
Definition: MuonTrackType.h:38
dqm::impl::MonitorElement::Fill
void Fill(long long x)
Definition: MonitorElement.h:242
FTS
FreeTrajectoryState FTS
Definition: RecoMuonValidator.cc:23
RecoMuonValidator::tpSelector_
TrackingParticleSelector tpSelector_
Definition: RecoMuonValidator.h:73
TrackingParticle::charge
float charge() const
Electric charge. Note this is taken from the first SimTrack only.
Definition: TrackingParticle.h:98
edm::View::begin
const_iterator begin() const
Muon
Definition: Muon.py:1
RecoMuonValidator::theMuonService
MuonServiceProxy * theMuonService
Definition: RecoMuonValidator.h:66
RecoMuonValidator::muAssocLabel_
edm::InputTag muAssocLabel_
Definition: RecoMuonValidator.h:53
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
reco::MuonToTrackingParticleAssociator::associateMuons
void associateMuons(MuonToSimCollection &recoToSim, SimToMuonCollection &simToReco, const edm::RefToBaseVector< reco::Muon > &muons, MuonTrackType type, const edm::RefVector< TrackingParticleCollection > &tpColl) const
Definition: MuonToTrackingParticleAssociator.h:25
RecoMuonValidator::countMuonHits
virtual int countMuonHits(const reco::Track &track) const
Definition: RecoMuonValidator.cc:1076
RecoMuonValidator::hDim
HistoDimensions hDim
Definition: RecoMuonValidator.h:127
RecoMuonValidator::outputFileName_
std::string outputFileName_
Definition: RecoMuonValidator.h:61
RecoMuonValidator::beamspotToken_
edm::EDGetTokenT< reco::BeamSpot > beamspotToken_
Definition: RecoMuonValidator.h:58
TrackingParticle::p
double p() const
Magnitude of momentum vector. Note this is taken from the first SimTrack only.
Definition: TrackingParticle.h:112
reco::Vertex::Point
math::XYZPoint Point
point in the space
Definition: Vertex.h:39
reco::BeamSpot::BeamWidthX
double BeamWidthX() const
beam width X
Definition: BeamSpot.h:82
edm::Ref::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:232
RecoMuonValidator::verbose_
unsigned int verbose_
Definition: RecoMuonValidator.h:45
reco::Track::recHitsBegin
trackingRecHit_iterator recHitsBegin() const
Iterator to first hit on the track.
Definition: Track.h:88
TrackerOfflineValidation_Dqm_cff.dirName
dirName
Definition: TrackerOfflineValidation_Dqm_cff.py:55
RecoMuonValidator::dqmBeginRun
void dqmBeginRun(const edm::Run &, const edm::EventSetup &eventSetup) override
Definition: RecoMuonValidator.cc:804
TrackingRecHit::ConstRecHitContainer
std::vector< ConstRecHitPointer > ConstRecHitContainer
Definition: TrackingRecHit.h:33
RecoMuonValidator::muonLabel_
edm::InputTag muonLabel_
Definition: RecoMuonValidator.h:48
HLTObjectMonitor_cfi.muonPt
muonPt
Definition: HLTObjectMonitor_cfi.py:82
RecoMuonValidator::MuonME::fill
void fill(const TrackingParticle *simRef, const Muon *muonRef)
Definition: RecoMuonValidator.cc:379
DetId
Definition: DetId.h:17
reco::LeafCandidate::p
double p() const final
magnitude of momentum vector
Definition: LeafCandidate.h:123
Mu
Definition: L1GtObject.h:29
MuonServiceProxy_cff.MuonServiceProxy
MuonServiceProxy
Definition: MuonServiceProxy_cff.py:15
edm::Handle::product
T const * product() const
Definition: Handle.h:69
reco::BeamSpot::sigmaZ
double sigmaZ() const
sigma z
Definition: BeamSpot.h:76
reco::BeamSpot::BeamWidthY
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:84
TrackingParticle::eta
double eta() const
Momentum pseudorapidity. Note this is taken from the first SimTrack only.
Definition: TrackingParticle.h:151
RecoMuonValidator::subsystemname_
std::string subsystemname_
Definition: RecoMuonValidator.h:63
TrackingParticle::vertex
Point vertex() const
Parent vertex position.
Definition: TrackingParticle.h:160
RecoMuonValidator::~RecoMuonValidator
~RecoMuonValidator() override
Definition: RecoMuonValidator.cc:795
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:45
edm
HLT enums.
Definition: AlignableModifier.h:19
edm::View::const_iterator
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:86
dqm::dqmstoreimpl::DQMStore::IBooker::book2D
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Definition: DQMStore.cc:266
position
static int position[264][3]
Definition: ReadPGInfo.cc:289
edm::OwnVector::const_iterator
Definition: OwnVector.h:41
edm::RefToBaseVector::push_back
void push_back(const RefToBase< T > &)
Definition: RefToBaseVector.h:217
RecoMuonValidator::muAssocToken_
edm::EDGetTokenT< reco::MuonToTrackingParticleAssociator > muAssocToken_
Definition: RecoMuonValidator.h:54
edm::RefVector::push_back
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
Definition: RefVector.h:67
RecoMuonValidator::beamspotLabel_
edm::InputTag beamspotLabel_
Definition: RecoMuonValidator.h:56
TrackingParticle
Monte Carlo truth information used for tracking validation.
Definition: TrackingParticle.h:29
RecoMuonValidator::primvertexLabel_
edm::InputTag primvertexLabel_
Definition: RecoMuonValidator.h:57
reco::BeamSpot::position
const Point & position() const
position
Definition: BeamSpot.h:59
edm::View::end
const_iterator end() const
PbPb_ZMuSkimMuonDPG_cff.trackType
trackType
Definition: PbPb_ZMuSkimMuonDPG_cff.py:36
muon::isTightMuon
bool isTightMuon(const reco::Muon &, const reco::Vertex &)
Definition: MuonSelectors.cc:913
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
RecoMuonValidator::analyze
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup) override
Definition: RecoMuonValidator.cc:820
reco::LeafCandidate::phi
double phi() const final
momentum azimuthal angle
Definition: LeafCandidate.h:148
reco::modules::TrackingParticleSelector
SingleObjectSelector< TrackingParticleCollection,::TrackingParticleSelector > TrackingParticleSelector
Definition: TrackingParticleSelector.cc:17
RecoMuonValidator::MuonME::bookHistos
void bookHistos(DQMStore::IBooker &ibooker, const string &dirName, const HistoDimensions &hDim)
Definition: RecoMuonValidator.cc:82
TrackingParticleRef
edm::Ref< TrackingParticleCollection > TrackingParticleRef
Definition: TrackingParticleFwd.h:10
RecoMuonValidator::RecoMuonValidator
RecoMuonValidator(const edm::ParameterSet &pset)
Definition: RecoMuonValidator.cc:555
dqm::dqmstoreimpl::DQMStore::save
void save(std::string const &filename, std::string const &path="", std::string const &pattern="", std::string const &rewrite="", uint32_t run=0, uint32_t lumi=0, SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, std::string const &fileupdate="RECREATE")
Definition: DQMStore.cc:2244
TrackingParticle::phi
double phi() const
Momentum azimuthal angle. Note this is taken from the first SimTrack only.
Definition: TrackingParticle.h:145
event
Definition: event.py:1
edm::Run
Definition: Run.h:45
HLT_2018_cff.track
track
Definition: HLT_2018_cff.py:10339
HLTTauDQMOffline_cfi.Muons
Muons
Definition: HLTTauDQMOffline_cfi.py:49
DetId::det
constexpr Detector det() const
get the detector field from this detid
Definition: DetId.h:46
reco::Track::recHitsEnd
trackingRecHit_iterator recHitsEnd() const
Iterator to last hit on the track.
Definition: Track.h:91
Generated for CMSSW Reference Manual by doxygen 1.8.11