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  // dump cfg parameters
558  edm::LogVerbatim("RecoMuonValidator") << "constructing RecoMuonValidator: " << pset.dump();
559  verbose_ = pset.getUntrackedParameter<unsigned int>("verbose", 0);
560 
561  outputFileName_ = pset.getUntrackedParameter<string>("outputFileName", "");
562 
563  wantTightMuon_ = pset.getParameter<bool>("wantTightMuon");
564  beamspotLabel_ = pset.getParameter<edm::InputTag>("beamSpot");
565  primvertexLabel_ = pset.getParameter<edm::InputTag>("primaryVertex");
566  beamspotToken_ = consumes<reco::BeamSpot>(beamspotLabel_);
567  primvertexToken_ = consumes<reco::VertexCollection>(primvertexLabel_);
568 
569  // Set histogram dimensions from config
570 
571  hDim.nBinP = pset.getUntrackedParameter<unsigned int>("nBinP");
572  hDim.minP = pset.getUntrackedParameter<double>("minP");
573  hDim.maxP = pset.getUntrackedParameter<double>("maxP");
574 
575  hDim.nBinPt = pset.getUntrackedParameter<unsigned int>("nBinPt");
576  hDim.minPt = pset.getUntrackedParameter<double>("minPt");
577  hDim.maxPt = pset.getUntrackedParameter<double>("maxPt");
578 
579  doAbsEta_ = pset.getUntrackedParameter<bool>("doAbsEta");
580  hDim.doAbsEta = doAbsEta_;
581  hDim.nBinEta = pset.getUntrackedParameter<unsigned int>("nBinEta");
582  hDim.minEta = pset.getUntrackedParameter<double>("minEta");
583  hDim.maxEta = pset.getUntrackedParameter<double>("maxEta");
584 
585  hDim.nBinDxy = pset.getUntrackedParameter<unsigned int>("nBinDxy");
586  hDim.minDxy = pset.getUntrackedParameter<double>("minDxy");
587  hDim.maxDxy = pset.getUntrackedParameter<double>("maxDxy");
588 
589  hDim.nBinDz = pset.getUntrackedParameter<unsigned int>("nBinDz");
590  hDim.minDz = pset.getUntrackedParameter<double>("minDz");
591  hDim.maxDz = pset.getUntrackedParameter<double>("maxDz");
592 
593  hDim.nBinPhi = pset.getUntrackedParameter<unsigned int>("nBinPhi");
594  hDim.minPhi = pset.getUntrackedParameter<double>("minPhi", -TMath::Pi());
595  hDim.maxPhi = pset.getUntrackedParameter<double>("maxPhi", TMath::Pi());
596 
597  hDim.nBinErr = pset.getUntrackedParameter<unsigned int>("nBinErr");
598  hDim.nBinPull = pset.getUntrackedParameter<unsigned int>("nBinPull");
599 
600  hDim.wPull = pset.getUntrackedParameter<double>("wPull");
601 
602  hDim.minErrP = pset.getUntrackedParameter<double>("minErrP");
603  hDim.maxErrP = pset.getUntrackedParameter<double>("maxErrP");
604 
605  hDim.minErrPt = pset.getUntrackedParameter<double>("minErrPt");
606  hDim.maxErrPt = pset.getUntrackedParameter<double>("maxErrPt");
607 
608  hDim.minErrQPt = pset.getUntrackedParameter<double>("minErrQPt");
609  hDim.maxErrQPt = pset.getUntrackedParameter<double>("maxErrQPt");
610 
611  hDim.minErrEta = pset.getUntrackedParameter<double>("minErrEta");
612  hDim.maxErrEta = pset.getUntrackedParameter<double>("maxErrEta");
613 
614  hDim.minErrPhi = pset.getUntrackedParameter<double>("minErrPhi");
615  hDim.maxErrPhi = pset.getUntrackedParameter<double>("maxErrPhi");
616 
617  hDim.minErrDxy = pset.getUntrackedParameter<double>("minErrDxy");
618  hDim.maxErrDxy = pset.getUntrackedParameter<double>("maxErrDxy");
619 
620  hDim.minErrDz = pset.getUntrackedParameter<double>("minErrDz");
621  hDim.maxErrDz = pset.getUntrackedParameter<double>("maxErrDz");
622 
623  hDim.nTrks = pset.getUntrackedParameter<unsigned int>("nTrks");
624  hDim.nAssoc = pset.getUntrackedParameter<unsigned int>("nAssoc");
625  hDim.nDof = pset.getUntrackedParameter<unsigned int>("nDof", 55);
626 
627  // Labels for simulation and reconstruction tracks
628  simLabel_ = pset.getParameter<InputTag>("simLabel");
629  tpRefVector = pset.getParameter<bool>("tpRefVector");
630  if (tpRefVector)
631  tpRefVectorToken_ = consumes<TrackingParticleRefVector>(simLabel_);
632  else
633  simToken_ = consumes<TrackingParticleCollection>(simLabel_);
634 
635  muonLabel_ = pset.getParameter<InputTag>("muonLabel");
636  muonToken_ = consumes<edm::View<reco::Muon> >(muonLabel_);
637 
638  // Labels for sim-reco association
639  doAssoc_ = pset.getUntrackedParameter<bool>("doAssoc", true);
640  muAssocLabel_ = pset.getParameter<InputTag>("muAssocLabel");
641  if (doAssoc_) {
642  muAssocToken_ = consumes<reco::MuonToTrackingParticleAssociator>(muAssocLabel_);
643  }
644 
645  // Different momentum assignment and additional histos in case of PF muons
646  usePFMuon_ = pset.getUntrackedParameter<bool>("usePFMuon");
647  hDim.usePFMuon = usePFMuon_;
648 
649  //type of track
650  std::string trackType = pset.getParameter<std::string>("trackType");
651  if (trackType == "inner")
652  trackType_ = reco::InnerTk;
653  else if (trackType == "outer")
654  trackType_ = reco::OuterTk;
655  else if (trackType == "global")
656  trackType_ = reco::GlobalTk;
657  else if (trackType == "segments")
658  trackType_ = reco::Segments;
659  else
660  throw cms::Exception("Configuration") << "Track type '" << trackType << "' not supported.\n";
661 
662  // seedPropagatorName_ = pset.getParameter<string>("SeedPropagator");
663 
664  ParameterSet tpset = pset.getParameter<ParameterSet>("tpSelector");
665  tpSelector_ = TrackingParticleSelector(tpset.getParameter<double>("ptMin"),
666  tpset.getParameter<double>("ptMax"),
667  tpset.getParameter<double>("minRapidity"),
668  tpset.getParameter<double>("maxRapidity"),
669  tpset.getParameter<double>("tip"),
670  tpset.getParameter<double>("lip"),
671  tpset.getParameter<int>("minHit"),
672  tpset.getParameter<bool>("signalOnly"),
673  tpset.getParameter<bool>("intimeOnly"),
674  tpset.getParameter<bool>("chargedOnly"),
675  tpset.getParameter<bool>("stableOnly"),
676  tpset.getParameter<std::vector<int> >("pdgId"));
677 
678  // the service parameters
679  ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
680 
681  // retrieve the instance of DQMService
682  dbe_ = Service<DQMStore>().operator->();
683  subsystemname_ = pset.getUntrackedParameter<std::string>("subSystemFolder", "YourSubsystem");
684 
685  subDir_ = pset.getUntrackedParameter<string>("subDir");
686  if (subDir_.empty())
687  subDir_ = "RecoMuonV";
688  if (subDir_[subDir_.size() - 1] == '/')
689  subDir_.erase(subDir_.size() - 1);
690 }
691 
692 void RecoMuonValidator::bookHistograms(DQMStore::IBooker& ibooker,
693  edm::Run const& iRun,
694  edm::EventSetup const& /* iSetup */) {
695  // book histograms
696  ibooker.cd();
697 
698  ibooker.setCurrentFolder(subDir_);
699 
700  commonME_ = new CommonME;
701  muonME_ = new MuonME;
702 
703  //commonME
704  const int nHits = 100;
705 
706  // - diffs
707  commonME_->hTrkToGlbDiffNTrackerHits_ = ibooker.book1D("TrkGlbDiffNTrackerHits",
708  "Difference of number of tracker hits (tkMuon - globalMuon)",
709  2 * nHits + 1,
710  -nHits - 0.5,
711  nHits + 0.5);
712  commonME_->hStaToGlbDiffNMuonHits_ = ibooker.book1D("StaGlbDiffNMuonHits",
713  "Difference of number of muon hits (staMuon - globalMuon)",
714  2 * nHits + 1,
715  -nHits - 0.5,
716  nHits + 0.5);
717 
718  commonME_->hTrkToGlbDiffNTrackerHitsEta_ =
719  ibooker.book2D("TrkGlbDiffNTrackerHitsEta",
720  "Difference of number of tracker hits (tkMuon - globalMuon)",
721  hDim.nBinEta,
722  hDim.minEta,
723  hDim.maxEta,
724  2 * nHits + 1,
725  -nHits - 0.5,
726  nHits + 0.5);
727  commonME_->hStaToGlbDiffNMuonHitsEta_ = ibooker.book2D("StaGlbDiffNMuonHitsEta",
728  "Difference of number of muon hits (staMuon - globalMuon)",
729  hDim.nBinEta,
730  hDim.minEta,
731  hDim.maxEta,
732  2 * nHits + 1,
733  -nHits - 0.5,
734  nHits + 0.5);
735 
736  commonME_->hTrkToGlbDiffNTrackerHitsPt_ = ibooker.book2D("TrkGlbDiffNTrackerHitsPt",
737  "Difference of number of tracker hits (tkMuon - globalMuon)",
738  hDim.nBinPt,
739  hDim.minPt,
740  hDim.maxPt,
741  2 * nHits + 1,
742  -nHits - 0.5,
743  nHits + 0.5);
744  commonME_->hStaToGlbDiffNMuonHitsPt_ = ibooker.book2D("StaGlbDiffNMuonHitsPt",
745  "Difference of number of muon hits (staMuon - globalMuon)",
746  hDim.nBinPt,
747  hDim.minPt,
748  hDim.maxPt,
749  2 * nHits + 1,
750  -nHits - 0.5,
751  nHits + 0.5);
752 
753  // -global muon hit pattern
754  commonME_->hNInvalidHitsGTHitPattern_ = ibooker.book1D(
755  "NInvalidHitsGTHitPattern", "Number of invalid hits on a global track", nHits + 1, -0.5, nHits + 0.5);
756  commonME_->hNInvalidHitsITHitPattern_ = ibooker.book1D(
757  "NInvalidHitsITHitPattern", "Number of invalid hits on an inner track", nHits + 1, -0.5, nHits + 0.5);
758  commonME_->hNInvalidHitsOTHitPattern_ = ibooker.book1D(
759  "NInvalidHitsOTHitPattern", "Number of invalid hits on an outer track", nHits + 1, -0.5, nHits + 0.5);
760  commonME_->hNDeltaInvalidHitsHitPattern_ =
761  ibooker.book1D("hNDeltaInvalidHitsHitPattern",
762  "The discrepancy for Number of invalid hits on an global track and inner and outer tracks",
763  2 * nHits + 1,
764  -nHits - 0.5,
765  nHits + 0.5);
766 
767  //muon based kinematics
768  commonME_->hMuonP_ = ibooker.book1D("PMuon", "p of muon", hDim.nBinP, hDim.minP, hDim.maxP);
769  commonME_->hMuonPt_ = ibooker.book1D("PtMuon", "p_{T} of muon", hDim.nBinPt, hDim.minPt, hDim.maxPt);
770  commonME_->hMuonEta_ = ibooker.book1D("EtaMuon", "#eta of muon", hDim.nBinEta, hDim.minEta, hDim.maxEta);
771  commonME_->hMuonPhi_ = ibooker.book1D("PhiMuon", "#phi of muon", hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
772  //track based kinematics
773  commonME_->hMuonTrackP_ = ibooker.book1D("PMuonTrack", "p of reco muon track", hDim.nBinP, hDim.minP, hDim.maxP);
774  commonME_->hMuonTrackPt_ =
775  ibooker.book1D("PtMuonTrack", "p_{T} of reco muon track", hDim.nBinPt, hDim.minPt, hDim.maxPt);
776  commonME_->hMuonTrackEta_ =
777  ibooker.book1D("EtaMuonTrack", "#eta of reco muon track", hDim.nBinEta, hDim.minEta, hDim.maxEta);
778  commonME_->hMuonTrackPhi_ =
779  ibooker.book1D("PhiMuonTrack", "#phi of reco muon track", hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
780  commonME_->hMuonTrackDxy_ =
781  ibooker.book1D("DxyMuonTrack", "Dxy of reco muon track", hDim.nBinDxy, hDim.minDxy, hDim.maxDxy);
782  commonME_->hMuonTrackDz_ =
783  ibooker.book1D("DzMuonTrack", "Dz of reco muon track", hDim.nBinDz, hDim.minDz, hDim.maxDz);
784 
785  //histograms for fractions
786  commonME_->hMuonAllP_ = ibooker.book1D("PMuonAll", "p of muons of all types", hDim.nBinP, hDim.minP, hDim.maxP);
787  commonME_->hMuonAllPt_ =
788  ibooker.book1D("PtMuonAll", "p_{T} of muon of all types", hDim.nBinPt, hDim.minPt, hDim.maxPt);
789  commonME_->hMuonAllEta_ =
790  ibooker.book1D("EtaMuonAll", "#eta of muon of all types", hDim.nBinEta, hDim.minEta, hDim.maxEta);
791  commonME_->hMuonAllPhi_ =
792  ibooker.book1D("PhiMuonAll", "#phi of muon of all types", hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
793 
794  muonME_->bookHistos(ibooker, subDir_, hDim);
795 }
796 
797 //
798 //Destructor
799 //
800 RecoMuonValidator::~RecoMuonValidator() {}
801 
802 //
803 //Begin run
804 //
805 
806 void RecoMuonValidator::dqmBeginRun(const edm::Run&, const EventSetup& eventSetup) {}
807 
808 //
809 //End run
810 //
811 void RecoMuonValidator::dqmEndRun(edm::Run const&, edm::EventSetup const&) {
812  if (dbe_ && !outputFileName_.empty())
813  dbe_->save(outputFileName_);
814 }
815 
816 //
817 //Analyze
818 //
819 void RecoMuonValidator::analyze(const Event& event, const EventSetup& eventSetup) {
820  // Look for the Primary Vertex (and use the BeamSpot instead, if you can't find it):
821  reco::Vertex::Point posVtx;
822  reco::Vertex::Error errVtx;
823  edm::Handle<reco::VertexCollection> recVtxs;
824  event.getByToken(primvertexToken_, recVtxs);
825  unsigned int theIndexOfThePrimaryVertex = 999.;
826  for (unsigned int ind = 0; ind < recVtxs->size(); ++ind) {
827  if ((*recVtxs)[ind].isValid() && !((*recVtxs)[ind].isFake())) {
828  theIndexOfThePrimaryVertex = ind;
829  break;
830  }
831  }
832  if (theIndexOfThePrimaryVertex < 100) {
833  posVtx = ((*recVtxs)[theIndexOfThePrimaryVertex]).position();
834  errVtx = ((*recVtxs)[theIndexOfThePrimaryVertex]).error();
835  } else {
836  LogInfo("RecoMuonValidator") << "reco::PrimaryVertex not found, use BeamSpot position instead\n";
837  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
838  event.getByToken(beamspotToken_, recoBeamSpotHandle);
839  reco::BeamSpot bs = *recoBeamSpotHandle;
840  posVtx = bs.position();
841  errVtx(0, 0) = bs.BeamWidthX();
842  errVtx(1, 1) = bs.BeamWidthY();
843  errVtx(2, 2) = bs.sigmaZ();
844  }
845  const reco::Vertex thePrimaryVertex(posVtx, errVtx);
846 
847  // Get TrackingParticles
848  TrackingParticleRefVector TPrefV;
849  const TrackingParticleRefVector* ptr_TPrefV = nullptr;
850  Handle<TrackingParticleCollection> simHandle;
851  Handle<TrackingParticleRefVector> TPCollectionRefVector_H;
852 
853  if (tpRefVector) {
854  event.getByToken(tpRefVectorToken_, TPCollectionRefVector_H);
855  ptr_TPrefV = TPCollectionRefVector_H.product();
856  TPrefV = *ptr_TPrefV;
857  } else {
858  event.getByToken(simToken_, simHandle);
859  size_t nTP = simHandle->size();
860  for (size_t i = 0; i < nTP; ++i) {
861  TPrefV.push_back(TrackingParticleRef(simHandle, i));
862  }
863  ptr_TPrefV = &TPrefV;
864  }
865 
866  // Get Muons
867  Handle<edm::View<Muon> > muonHandle;
868  event.getByToken(muonToken_, muonHandle);
869  View<Muon> muonColl = *(muonHandle.product());
870 
871  reco::MuonToTrackingParticleAssociator const* assoByHits = nullptr;
872  if (doAssoc_) {
873  edm::Handle<reco::MuonToTrackingParticleAssociator> associatorBase;
874  event.getByToken(muAssocToken_, associatorBase);
875  assoByHits = associatorBase.product();
876  }
877 
878  const size_t nSim = ptr_TPrefV->size();
879 
880  edm::RefToBaseVector<reco::Muon> Muons;
881  for (size_t i = 0; i < muonHandle->size(); ++i) {
882  Muons.push_back(muonHandle->refAt(i));
883  }
884 
885  muonME_->hNSim_->Fill(nSim);
886  muonME_->hNMuon_->Fill(muonColl.size());
887 
888  reco::MuonToSimCollection muonToSimColl;
889  reco::SimToMuonCollection simToMuonColl;
890 
891  if (doAssoc_) {
892  edm::LogVerbatim("RecoMuonValidator")
893  << "\n >>> MuonToSim association : " << muAssocLabel_ << " <<< \n"
894  << " muon collection : " << muonLabel_ << " (size = " << muonHandle->size() << ") \n"
895  << " TrackingParticle collection : " << simLabel_ << " (size = " << nSim << ")";
896 
897  assoByHits->associateMuons(muonToSimColl, simToMuonColl, Muons, trackType_, TPrefV);
898  } else {
899  /*
900  // SimToMuon associations
901  Handle<reco::RecoToSimCollection> simToTrkMuHandle;
902  event.getByLabel(trkMuAssocLabel_, simToTrkMuHandle);
903  trkSimRecColl = *(simToTrkMuHandle.product());
904 
905  Handle<reco::RecoToSimCollection> simToStaMuHandle;
906  event.getByLabel(staMuAssocLabel_, simToStaMuHandle);
907  staSimRecColl = *(simToStaMuHandle.product());
908 
909  Handle<reco::RecoToSimCollection> simToGlbMuHandle;
910  event.getByLabel(glbMuAssocLabel_, simToGlbMuHandle);
911  glbSimRecColl = *(simToGlbMuHandle.product());
912 
913  // MuonToSim associations
914  Handle<reco::SimToRecoCollection> trkMuToSimHandle;
915  event.getByLabel(trkMuAssocLabel_, trkMuToSimHandle);
916  trkRecSimColl = *(trkMuToSimHandle.product());
917 
918  Handle<reco::SimToRecoCollection> staMuToSimHandle;
919  event.getByLabel(staMuAssocLabel_, staMuToSimHandle);
920  staRecSimColl = *(staMuToSimHandle.product());
921 
922  Handle<reco::SimToRecoCollection> glbMuToSimHandle;
923  event.getByLabel(glbMuAssocLabel_, glbMuToSimHandle);
924  glbRecSimColl = *(glbMuToSimHandle.product());
925 */
926  }
927 
928  int glbNTrackerHits = 0;
929  int trkNTrackerHits = 0;
930  int glbNMuonHits = 0;
931  int staNMuonHits = 0;
932  int NTrackerHits = 0;
933  int NMuonHits = 0;
934 
935  // Analyzer reco::Muon
936  for (View<Muon>::const_iterator iMuon = muonColl.begin(); iMuon != muonColl.end(); ++iMuon) {
937  double muonP, muonPt, muonEta, muonPhi;
938  if (usePFMuon_) {
939  muonP = iMuon->pfP4().P();
940  muonPt = iMuon->pfP4().Pt();
941  muonEta = iMuon->pfP4().Eta();
942  muonPhi = iMuon->pfP4().Phi();
943  } else {
944  muonP = iMuon->p();
945  muonPt = iMuon->pt();
946  muonEta = iMuon->eta();
947  muonPhi = iMuon->phi();
948  }
949 
950  //histograms for fractions
951  commonME_->hMuonAllP_->Fill(muonP);
952  commonME_->hMuonAllPt_->Fill(muonPt);
953  commonME_->hMuonAllEta_->Fill(muonEta);
954  commonME_->hMuonAllPhi_->Fill(muonPhi);
955 
956  if (!selector_(*iMuon))
957  continue;
958  if (wantTightMuon_) {
959  if (!muon::isTightMuon(*iMuon, thePrimaryVertex))
960  continue;
961  }
962 
963  TrackRef Track = iMuon->track();
964 
965  if (Track.isNonnull()) {
966  commonME_->hMuonTrackP_->Fill(Track->p());
967  commonME_->hMuonTrackPt_->Fill(Track->pt());
968  commonME_->hMuonTrackEta_->Fill(Track->eta());
969  commonME_->hMuonTrackPhi_->Fill(Track->phi());
970 
971  //ip histograms
972  commonME_->hMuonTrackDxy_->Fill(Track->dxy());
973  commonME_->hMuonTrackDz_->Fill(Track->dz());
974  }
975 
976  if (iMuon->isGlobalMuon()) {
977  Track = iMuon->combinedMuon();
978  glbNTrackerHits = countTrackerHits(*Track);
979  glbNMuonHits = countMuonHits(*Track);
980  } else if (iMuon->isTrackerMuon()) {
981  Track = iMuon->track();
982  trkNTrackerHits = countTrackerHits(*Track);
983  } else {
984  Track = iMuon->standAloneMuon();
985  }
986 
987  NTrackerHits = countTrackerHits(*Track);
988  muonME_->hNTrackerHits_->Fill(NTrackerHits);
989  muonME_->hNTrackerHits_vs_Pt_->Fill(Track->pt(), NTrackerHits);
990  muonME_->hNTrackerHits_vs_Eta_->Fill(Track->eta(), NTrackerHits);
991 
992  NMuonHits = countMuonHits(*Track);
993  muonME_->hNMuonHits_->Fill(NMuonHits);
994  muonME_->hNMuonHits_vs_Pt_->Fill(Track->pt(), NMuonHits);
995  muonME_->hNMuonHits_vs_Eta_->Fill(Track->eta(), NMuonHits);
996 
997  //list of histos for each type
998 
999  // muonME_->hNTrks_->Fill();
1000  muonME_->hNTrksEta_->Fill(Track->eta());
1001  muonME_->hNTrksPt_->Fill(Track->pt());
1002 
1003  commonME_->hMuonP_->Fill(muonP);
1004  commonME_->hMuonPt_->Fill(muonPt);
1005  commonME_->hMuonEta_->Fill(muonEta);
1006  commonME_->hMuonPhi_->Fill(muonPhi);
1007 
1008  if (iMuon->isGlobalMuon()) {
1009  double gtHitPat = iMuon->globalTrack()->hitPattern().numberOfAllHits(HitPattern::TRACK_HITS) -
1010  iMuon->globalTrack()->hitPattern().numberOfValidHits();
1011 
1012  double itHitPat = iMuon->innerTrack()->hitPattern().numberOfAllHits(HitPattern::TRACK_HITS) -
1013  iMuon->innerTrack()->hitPattern().numberOfValidHits();
1014 
1015  double otHitPat = iMuon->outerTrack()->hitPattern().numberOfAllHits(HitPattern::TRACK_HITS) -
1016  iMuon->outerTrack()->hitPattern().numberOfValidHits();
1017 
1018  commonME_->hNInvalidHitsGTHitPattern_->Fill(gtHitPat);
1019  commonME_->hNInvalidHitsITHitPattern_->Fill(itHitPat);
1020  commonME_->hNInvalidHitsOTHitPattern_->Fill(otHitPat);
1021  commonME_->hNDeltaInvalidHitsHitPattern_->Fill(gtHitPat - itHitPat - otHitPat);
1022 
1023  //must be global and standalone
1024  if (iMuon->isStandAloneMuon()) {
1025  commonME_->hStaToGlbDiffNMuonHitsEta_->Fill(Track->eta(), staNMuonHits - glbNMuonHits);
1026  commonME_->hStaToGlbDiffNMuonHitsPt_->Fill(Track->pt(), staNMuonHits - glbNMuonHits);
1027  commonME_->hStaToGlbDiffNMuonHits_->Fill(staNMuonHits - glbNMuonHits);
1028  }
1029 
1030  //must be global and tracker
1031  if (iMuon->isTrackerMuon()) {
1032  commonME_->hTrkToGlbDiffNTrackerHitsEta_->Fill(Track->eta(), trkNTrackerHits - glbNTrackerHits);
1033  commonME_->hTrkToGlbDiffNTrackerHitsPt_->Fill(Track->pt(), trkNTrackerHits - glbNTrackerHits);
1034  commonME_->hTrkToGlbDiffNTrackerHits_->Fill(trkNTrackerHits - glbNTrackerHits);
1035  }
1036  }
1037 
1038  } //end of reco muon loop
1039 
1040  // Associate by hits
1041  for (size_t i = 0; i < nSim; i++) {
1042  TrackingParticleRef simRef = TPrefV[i];
1043  const TrackingParticle* simTP = simRef.get();
1044  if (!tpSelector_(*simTP))
1045  continue;
1046 
1047  //denominators for efficiency plots
1048  const double simP = simRef->p();
1049  const double simPt = simRef->pt();
1050  const double simEta = doAbsEta_ ? fabs(simRef->eta()) : simRef->eta();
1051  const double simPhi = simRef->phi();
1052 
1053  GlobalPoint simVtx(simRef->vertex().x(), simRef->vertex().y(), simRef->vertex().z());
1054  GlobalVector simMom(simRef->momentum().x(), simRef->momentum().y(), simRef->momentum().z());
1055  const double simDxy = -simVtx.x() * sin(simPhi) + simVtx.y() * cos(simPhi);
1056  const double simDz = simVtx.z() - (simVtx.x() * simMom.x() + simVtx.y() * simMom.y()) * simMom.z() / simMom.perp2();
1057 
1058  const unsigned int nSimHits = simRef->numberOfHits();
1059 
1060  muonME_->hSimP_->Fill(simP);
1061  muonME_->hSimPt_->Fill(simPt);
1062  muonME_->hSimEta_->Fill(simEta);
1063  muonME_->hSimPhi_->Fill(simPhi);
1064  muonME_->hSimDxy_->Fill(simDxy);
1065  muonME_->hSimDz_->Fill(simDz);
1066  muonME_->hNSimHits_->Fill(nSimHits);
1067 
1068  // Get sim-reco association for a simRef
1069  vector<pair<RefToBase<Muon>, double> > MuRefV;
1070  if (simToMuonColl.find(simRef) != simToMuonColl.end()) {
1071  MuRefV = simToMuonColl[simRef];
1072 
1073  if (!MuRefV.empty()) {
1074  muonME_->hNSimToReco_->Fill(MuRefV.size());
1075  const Muon* Mu = MuRefV.begin()->first.get();
1076  if (!selector_(*Mu))
1077  continue;
1078  if (wantTightMuon_) {
1079  if (!muon::isTightMuon(*Mu, thePrimaryVertex))
1080  continue;
1081  }
1082 
1083  muonME_->fill(&*simTP, Mu);
1084  }
1085  }
1086  }
1087 }
1088 
1089 int RecoMuonValidator::countMuonHits(const reco::Track& track) const {
1090  TransientTrackingRecHit::ConstRecHitContainer result;
1091 
1092  int count = 0;
1093 
1094  for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
1095  if ((*hit)->isValid()) {
1096  DetId recoid = (*hit)->geographicalId();
1097  if (recoid.det() == DetId::Muon)
1098  count++;
1099  }
1100  }
1101  return count;
1102 }
1103 
1104 int RecoMuonValidator::countTrackerHits(const reco::Track& track) const {
1105  TransientTrackingRecHit::ConstRecHitContainer result;
1106 
1107  int count = 0;
1108 
1109  for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
1110  if ((*hit)->isValid()) {
1111  DetId recoid = (*hit)->geographicalId();
1112  if (recoid.det() == DetId::Tracker)
1113  count++;
1114  }
1115  }
1116  return count;
1117 }
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:134
RecoMuonValidator::muonToken_
edm::EDGetTokenT< edm::View< reco::Muon > > muonToken_
Definition: RecoMuonValidator.h:51
RecoMuonValidator::simToken_
edm::EDGetTokenT< TrackingParticleCollection > simToken_
Definition: RecoMuonValidator.h:47
reco::MuonToSimCollection
std::map< edm::RefToBase< reco::Muon >, std::vector< std::pair< TrackingParticleRef, double > >, RefToBaseSort > MuonToSimCollection
Definition: MuonTrackType.h:37
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
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:27
mps_fire.i
i
Definition: mps_fire.py:429
TrackingParticle::momentum
Vector momentum() const
spatial momentum vector
Definition: TrackingParticle.h:106
reco::LeafCandidate::pt
double pt() const final
transverse momentum
Definition: LeafCandidate.h:146
mps_fire.result
result
Definition: mps_fire.py:311
TSOS
TrajectoryStateOnSurface TSOS
Definition: RecoMuonValidator.cc:22
dqm::implementation::NavigatorBase::setCurrentFolder
virtual void setCurrentFolder(std::string const &fullpath)
Definition: DQMStore.cc:36
TrackingParticle::p
double p() const
Magnitude of momentum vector. Note this is taken from the first SimTrack only.
Definition: TrackingParticle.h:112
HLTObjectMonitor_cfi.muonEta
muonEta
Definition: HLTObjectMonitor_cfi.py:92
reco::TrackBase::p
double p() const
momentum vector magnitude
Definition: TrackBase.h:631
Mu
Definition: L1GtObject.h:29
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
edm::Handle::product
T const * product() const
Definition: Handle.h:70
RecoMuonValidator::simLabel_
edm::InputTag simLabel_
Definition: RecoMuonValidator.h:44
RecoMuonValidator::selector_
StringCutObjectSelector< reco::Muon > selector_
Definition: RecoMuonValidator.h:129
reco::Vertex::Error
math::Error< dimension >::type Error
covariance error matrix (3x3)
Definition: Vertex.h:45
muonDTDigis_cfi.pset
pset
Definition: muonDTDigis_cfi.py:27
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
HLTObjectMonitor_cfi.muonPhi
muonPhi
Definition: HLTObjectMonitor_cfi.py:102
RecoMuonValidator::countMuonHits
virtual int countMuonHits(const reco::Track &track) const
Definition: RecoMuonValidator.cc:1089
RecoMuonValidator::trackType_
reco::MuonTrackType trackType_
Definition: RecoMuonValidator.h:75
RecoMuonValidator::primvertexToken_
edm::EDGetTokenT< reco::VertexCollection > primvertexToken_
Definition: RecoMuonValidator.h:59
reco::Muon::track
TrackRef track() const override
reference to a Track
Definition: Muon.h:46
HLT_2024v12_cff.track
track
Definition: HLT_2024v12_cff.py:9556
DetId::det
constexpr Detector det() const
get the detector field from this detid
Definition: DetId.h:46
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
RecoMuonValidator::bookHistograms
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
Definition: RecoMuonValidator.cc:692
RecoMuonValidator::countTrackerHits
virtual int countTrackerHits(const reco::Track &track) const
Definition: RecoMuonValidator.cc:1104
RecoMuonValidator::dqmEndRun
void dqmEndRun(edm::Run const &, edm::EventSetup const &) override
Definition: RecoMuonValidator.cc:811
RecoMuonValidator::tpRefVectorToken_
edm::EDGetTokenT< TrackingParticleRefVector > tpRefVectorToken_
Definition: RecoMuonValidator.h:48
edm::View::size
size_type size() const
reco::SimToMuonCollection
std::map< TrackingParticleRef, std::vector< std::pair< edm::RefToBase< reco::Muon >, double > > > SimToMuonCollection
Definition: MuonTrackType.h:38
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
dqm::impl::MonitorElement::Fill
void Fill(long long x)
Definition: MonitorElement.h:292
FTS
FreeTrajectoryState FTS
Definition: RecoMuonValidator.cc:23
RecoMuonValidator::tpSelector_
TrackingParticleSelector tpSelector_
Definition: RecoMuonValidator.h:72
reco::TrackBase::pt
double pt() const
track transverse momentum
Definition: TrackBase.h:637
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
reco::TrackBase::dz
double dz() const
dz parameter (= dsz/cos(lambda)). This is the track z0 w.r.t (0,0,0) only if the refPoint is close to...
Definition: TrackBase.h:622
reco::Muon::pfP4
reco::Candidate::LorentzVector pfP4() const
Definition: Muon.h:97
reco::LeafCandidate::p
double p() const final
magnitude of momentum vector
Definition: LeafCandidate.h:123
Muon
Definition: Muon.py:1
RecoMuonValidator::muAssocLabel_
edm::InputTag muAssocLabel_
Definition: RecoMuonValidator.h:53
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
RecoMuonValidator::hDim
HistoDimensions hDim
Definition: RecoMuonValidator.h:126
RecoMuonValidator::outputFileName_
std::string outputFileName_
Definition: RecoMuonValidator.h:61
RecoMuonValidator::beamspotToken_
edm::EDGetTokenT< reco::BeamSpot > beamspotToken_
Definition: RecoMuonValidator.h:58
reco::Vertex::Point
math::XYZPoint Point
point in the space
Definition: Vertex.h:40
reco::TrackBase::phi
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:649
RecoMuonValidator::verbose_
unsigned int verbose_
Definition: RecoMuonValidator.h:42
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:806
TrackingRecHit::ConstRecHitContainer
std::vector< ConstRecHitPointer > ConstRecHitContainer
Definition: TrackingRecHit.h:32
RecoMuonValidator::muonLabel_
edm::InputTag muonLabel_
Definition: RecoMuonValidator.h:45
TrackingParticle::eta
double eta() const
Momentum pseudorapidity. Note this is taken from the first SimTrack only.
Definition: TrackingParticle.h:154
HLTObjectMonitor_cfi.muonPt
muonPt
Definition: HLTObjectMonitor_cfi.py:82
RecoMuonValidator::MuonME::fill
void fill(const TrackingParticle *simRef, const Muon *muonRef)
Definition: RecoMuonValidator.cc:379
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:131
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:652
DetId
Definition: DetId.h:17
edm::RefVector::size
size_type size() const
Size of the RefVector.
Definition: RefVector.h:102
dqm::implementation::IBooker::book2D
MonitorElement * book2D(TString const &name, TString const &title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, FUNC onbooking=NOOP())
Definition: DQMStore.h:212
RecoMuonValidator::subsystemname_
std::string subsystemname_
Definition: RecoMuonValidator.h:63
TrackingParticle::phi
double phi() const
Momentum azimuthal angle. Note this is taken from the first SimTrack only.
Definition: TrackingParticle.h:148
dqm::implementation::DQMStore::save
DQM_DEPRECATED void save(std::string const &filename, std::string const &path="")
Definition: DQMStore.cc:824
RecoMuonValidator::~RecoMuonValidator
~RecoMuonValidator() override
Definition: RecoMuonValidator.cc:800
reco
fixed size matrix
Definition: AlignmentAlgorithmBase.h:46
edm
HLT enums.
Definition: AlignableModifier.h:19
edm::View::const_iterator
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:88
position
static int position[264][3]
Definition: ReadPGInfo.cc:289
edm::OwnVector::const_iterator
Definition: OwnVector.h:41
RecoMuonValidator::muAssocToken_
edm::EDGetTokenT< reco::MuonToTrackingParticleAssociator > muAssocToken_
Definition: RecoMuonValidator.h:54
edm::ParameterSet::dump
std::string dump(unsigned int indent=0) const
Definition: ParameterSet.cc:888
edm::Ref::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:232
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
PbPb_ZMuSkimMuonDPG_cff.trackType
trackType
Definition: PbPb_ZMuSkimMuonDPG_cff.py:36
edm::View::begin
const_iterator begin() const
RecoMuonValidator::pset
edm::ParameterSet pset
Definition: RecoMuonValidator.h:64
TrackingParticle::pt
double pt() const
Transverse momentum. Note this is taken from the first SimTrack only.
Definition: TrackingParticle.h:145
muon::isTightMuon
bool isTightMuon(const reco::Muon &, const reco::Vertex &)
Definition: MuonSelectors.cc:876
dqm::implementation::IBooker::book1D
MonitorElement * book1D(TString const &name, TString const &title, int const nchX, double const lowX, double const highX, FUNC onbooking=NOOP())
Definition: DQMStore.h:98
nHits
TupleMultiplicity< TrackerTraits > const *__restrict__ uint32_t nHits
Definition: RiemannFitOnGPU.h:27
reco::LeafCandidate::phi
double phi() const final
momentum azimuthal angle
Definition: LeafCandidate.h:148
options_cfi.eventSetup
eventSetup
Definition: options_cfi.py:11
RecoMuonValidator::analyze
void analyze(const edm::Event &event, const edm::EventSetup &eventSetup) override
Definition: RecoMuonValidator.cc:819
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
edm::View::end
const_iterator end() const
TrackingParticleRef
edm::Ref< TrackingParticleCollection > TrackingParticleRef
Definition: TrackingParticleFwd.h:12
RecoMuonValidator::RecoMuonValidator
RecoMuonValidator(const edm::ParameterSet &pset)
Definition: RecoMuonValidator.cc:555
event
Definition: event.py:1
edm::Run
Definition: Run.h:45
reco::LeafCandidate::charge
int charge() const final
electric charge
Definition: LeafCandidate.h:106
TrackingParticle::vertex
Point vertex() const
Parent vertex position.
Definition: TrackingParticle.h:169
HLTTauDQMOffline_cfi.Muons
Muons
Definition: HLTTauDQMOffline_cfi.py:48
TrackingParticle::charge
float charge() const
Electric charge. Note this is taken from the first SimTrack only.
Definition: TrackingParticle.h:98
reco::TrackBase::dxy
double dxy() const
dxy parameter. (This is the transverse impact parameter w.r.t. to (0,0,0) ONLY if refPoint is close t...
Definition: TrackBase.h:608
reco::LeafCandidate::eta
double eta() const final
momentum pseudorapidity
Definition: LeafCandidate.h:152
Generated for CMSSW Reference Manual by doxygen 1.8.14