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