CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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>("chargedOnly"),
528  tpset.getParameter<bool>("stableOnly"),
529  tpset.getParameter<std::vector<int> >("pdgId"));
530 
531 
532  // the service parameters
533  ParameterSet serviceParameters
534  = pset.getParameter<ParameterSet>("ServiceParameters");
535  theMuonService = new MuonServiceProxy(serviceParameters);
536 
537  // retrieve the instance of DQMService
538  dbe_ = Service<DQMStore>().operator->();
539  subsystemname_ = pset.getUntrackedParameter<std::string>("subSystemFolder", "YourSubsystem") ;
540 
541  subDir_ = pset.getUntrackedParameter<string>("subDir");
542  if ( subDir_.empty() ) subDir_ = "RecoMuonV";
543  if ( subDir_[subDir_.size()-1] == '/' ) subDir_.erase(subDir_.size()-1);
544 
545 
546 }
547 
548 void RecoMuonValidator::bookHistograms(DQMStore::IBooker & ibooker,
549  edm::Run const & iRun,
550  edm::EventSetup const & /* iSetup */){
551 
552 
553  // book histograms
554  ibooker.cd();
555 
556  ibooker.setCurrentFolder(subDir_);
557 
558  commonME_ = new CommonME;
559  muonME_ = new MuonME;
560 
561  //commonME
562  const int nHits = 100;
563 
564  // - diffs
565  commonME_->hTrkToGlbDiffNTrackerHits_ = ibooker.book1D("TrkGlbDiffNTrackerHits", "Difference of number of tracker hits (tkMuon - globalMuon)", 2*nHits+1, -nHits-0.5, nHits+0.5);
566  commonME_->hStaToGlbDiffNMuonHits_ = ibooker.book1D("StaGlbDiffNMuonHits", "Difference of number of muon hits (staMuon - globalMuon)", 2*nHits+1, -nHits-0.5, nHits+0.5);
567 
568  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);
569  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);
570 
571  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);
572  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);
573 
574  // -global muon hit pattern
575  commonME_->hNInvalidHitsGTHitPattern_ = ibooker.book1D("NInvalidHitsGTHitPattern", "Number of invalid hits on a global track", nHits+1, -0.5, nHits+0.5);
576  commonME_->hNInvalidHitsITHitPattern_ = ibooker.book1D("NInvalidHitsITHitPattern", "Number of invalid hits on an inner track", nHits+1, -0.5, nHits+0.5);
577  commonME_->hNInvalidHitsOTHitPattern_ = ibooker.book1D("NInvalidHitsOTHitPattern", "Number of invalid hits on an outer track", nHits+1, -0.5, nHits+0.5);
578  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);
579 
580  //muon based kinematics
581  commonME_->hMuonP_ = ibooker.book1D("PMuon" , "p of muon" , hDim.nBinP , hDim.minP , hDim.maxP );
582  commonME_->hMuonPt_ = ibooker.book1D("PtMuon" , "p_{T} of muon", hDim.nBinPt , hDim.minPt , hDim.maxPt );
583  commonME_->hMuonEta_ = ibooker.book1D("EtaMuon", "#eta of muon" , hDim.nBinEta, hDim.minEta, hDim.maxEta);
584  commonME_->hMuonPhi_ = ibooker.book1D("PhiMuon", "#phi of muon" , hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
585  //track based kinematics
586  commonME_->hMuonTrackP_ = ibooker.book1D("PMuonTrack" , "p of reco muon track" , hDim.nBinP , hDim.minP , hDim.maxP );
587  commonME_->hMuonTrackPt_ = ibooker.book1D("PtMuonTrack" , "p_{T} of reco muon track", hDim.nBinPt , hDim.minPt , hDim.maxPt );
588  commonME_->hMuonTrackEta_ = ibooker.book1D("EtaMuonTrack", "#eta of reco muon track" , hDim.nBinEta, hDim.minEta, hDim.maxEta);
589  commonME_->hMuonTrackPhi_ = ibooker.book1D("PhiMuonTrack", "#phi of reco muon track" , hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
590  commonME_->hMuonTrackDxy_ = ibooker.book1D("DxyMuonTrack", "Dxy of reco muon track" , hDim.nBinDxy, hDim.minDxy, hDim.maxDxy);
591  commonME_->hMuonTrackDz_ = ibooker.book1D("DzMuonTrack", "Dz of reco muon track" , hDim.nBinDz, hDim.minDz, hDim.maxDz);
592 
593  //histograms for fractions
594  commonME_->hMuonAllP_ = ibooker.book1D("PMuonAll" , "p of muons of all types" , hDim.nBinP , hDim.minP , hDim.maxP );
595  commonME_->hMuonAllPt_ = ibooker.book1D("PtMuonAll" , "p_{T} of muon of all types", hDim.nBinPt , hDim.minPt , hDim.maxPt );
596  commonME_->hMuonAllEta_ = ibooker.book1D("EtaMuonAll", "#eta of muon of all types" , hDim.nBinEta, hDim.minEta, hDim.maxEta);
597  commonME_->hMuonAllPhi_ = ibooker.book1D("PhiMuonAll", "#phi of muon of all types" , hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
598 
599  muonME_->bookHistos(ibooker, subDir_, hDim);
600 }
601 
602 //
603 //Destructor
604 //
605 RecoMuonValidator::~RecoMuonValidator()
606 {
607  if ( theMuonService ) delete theMuonService;
608 }
609 
610 //
611 //Begin run
612 //
613 
614 void RecoMuonValidator::dqmBeginRun(const edm::Run& , const EventSetup& eventSetup)
615 {
616  if ( theMuonService ) theMuonService->update(eventSetup);
617 }
618 
619 //
620 //End run
621 //
622 void RecoMuonValidator::endRun()
623 {
624  if ( dbe_ && ! outputFileName_.empty() ) dbe_->save(outputFileName_);
625 }
626 
627 //
628 //Analyze
629 //
630 void RecoMuonValidator::analyze(const Event& event, const EventSetup& eventSetup)
631 {
632  // Look for the Primary Vertex (and use the BeamSpot instead, if you can't find it):
633  reco::Vertex::Point posVtx;
634  reco::Vertex::Error errVtx;
635  edm::Handle<reco::VertexCollection> recVtxs;
636  event.getByToken(primvertexToken_,recVtxs);
637  unsigned int theIndexOfThePrimaryVertex = 999.;
638  for (unsigned int ind=0; ind<recVtxs->size(); ++ind) {
639  if ( (*recVtxs)[ind].isValid() && !((*recVtxs)[ind].isFake()) ) {
640  theIndexOfThePrimaryVertex = ind;
641  break;
642  }
643  }
644  if (theIndexOfThePrimaryVertex<100) {
645  posVtx = ((*recVtxs)[theIndexOfThePrimaryVertex]).position();
646  errVtx = ((*recVtxs)[theIndexOfThePrimaryVertex]).error();
647  }
648  else {
649  LogInfo("RecoMuonValidator") << "reco::PrimaryVertex not found, use BeamSpot position instead\n";
650  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
651  event.getByToken(beamspotToken_,recoBeamSpotHandle);
652  reco::BeamSpot bs = *recoBeamSpotHandle;
653  posVtx = bs.position();
654  errVtx(0,0) = bs.BeamWidthX();
655  errVtx(1,1) = bs.BeamWidthY();
656  errVtx(2,2) = bs.sigmaZ();
657  }
658  const reco::Vertex thePrimaryVertex(posVtx,errVtx);
659 
660 
661  // Get TrackingParticles
662  Handle<TrackingParticleCollection> simHandle;
663  event.getByToken(simToken_, simHandle);
664  const TrackingParticleCollection simColl = *(simHandle.product());
665 
666  // Get Muons
667  Handle<edm::View<Muon> > muonHandle;
668  event.getByToken(muonToken_, muonHandle);
669  View<Muon> muonColl = *(muonHandle.product());
670 
671  reco::MuonToTrackingParticleAssociator const* assoByHits = nullptr;
672  if ( doAssoc_ ) {
673  edm::Handle<reco::MuonToTrackingParticleAssociator> associatorBase;
674  event.getByToken(muAssocToken_, associatorBase);
675  assoByHits = associatorBase.product();
676  }
677 
678  const TrackingParticleCollection::size_type nSim = simColl.size();
679 
680  edm::RefToBaseVector<reco::Muon> Muons;
681  for (size_t i = 0; i < muonHandle->size(); ++i) {
682  Muons.push_back(muonHandle->refAt(i));
683  }
684 
685  edm::RefVector<TrackingParticleCollection> allTPs;
686  for (size_t i = 0; i < nSim; ++i) {
687  allTPs.push_back(TrackingParticleRef(simHandle,i));
688  }
689 
690 
691  muonME_->hNSim_->Fill(nSim);
692  muonME_->hNMuon_->Fill(muonColl.size());
693 
694  reco::MuonToSimCollection muonToSimColl;
695  reco::SimToMuonCollection simToMuonColl;
696 
697 
698  if ( doAssoc_ ) {
699  assoByHits->associateMuons(muonToSimColl, simToMuonColl, Muons, trackType_, allTPs);
700  } else {
701 
702 /*
703  // SimToMuon associations
704  Handle<reco::RecoToSimCollection> simToTrkMuHandle;
705  event.getByLabel(trkMuAssocLabel_, simToTrkMuHandle);
706  trkSimRecColl = *(simToTrkMuHandle.product());
707 
708  Handle<reco::RecoToSimCollection> simToStaMuHandle;
709  event.getByLabel(staMuAssocLabel_, simToStaMuHandle);
710  staSimRecColl = *(simToStaMuHandle.product());
711 
712  Handle<reco::RecoToSimCollection> simToGlbMuHandle;
713  event.getByLabel(glbMuAssocLabel_, simToGlbMuHandle);
714  glbSimRecColl = *(simToGlbMuHandle.product());
715 
716  // MuonToSim associations
717  Handle<reco::SimToRecoCollection> trkMuToSimHandle;
718  event.getByLabel(trkMuAssocLabel_, trkMuToSimHandle);
719  trkRecSimColl = *(trkMuToSimHandle.product());
720 
721  Handle<reco::SimToRecoCollection> staMuToSimHandle;
722  event.getByLabel(staMuAssocLabel_, staMuToSimHandle);
723  staRecSimColl = *(staMuToSimHandle.product());
724 
725  Handle<reco::SimToRecoCollection> glbMuToSimHandle;
726  event.getByLabel(glbMuAssocLabel_, glbMuToSimHandle);
727  glbRecSimColl = *(glbMuToSimHandle.product());
728 */
729  }
730 
731 
732 
733  int glbNTrackerHits = 0; int trkNTrackerHits = 0;
734  int glbNMuonHits = 0; int staNMuonHits = 0;
735  int NTrackerHits = 0; int NMuonHits = 0;
736 
737 
738  // Analyzer reco::Muon
739  for(View<Muon>::const_iterator iMuon = muonColl.begin();
740  iMuon != muonColl.end(); ++iMuon) {
741 
742  double muonP, muonPt, muonEta, muonPhi;
743  if (usePFMuon_) {
744  muonP = iMuon->pfP4().P();
745  muonPt = iMuon->pfP4().Pt();
746  muonEta = iMuon->pfP4().Eta();
747  muonPhi = iMuon->pfP4().Phi();
748  }
749  else {
750  muonP = iMuon->p();
751  muonPt = iMuon->pt();
752  muonEta = iMuon->eta();
753  muonPhi = iMuon->phi();
754  }
755 
756  //histograms for fractions
757  commonME_->hMuonAllP_->Fill(muonP);
758  commonME_->hMuonAllPt_->Fill(muonPt);
759  commonME_->hMuonAllEta_->Fill(muonEta);
760  commonME_->hMuonAllPhi_->Fill(muonPhi);
761 
762  if (!selector_(*iMuon)) continue;
763  if (wantTightMuon_)
764  {
765  if (!muon::isTightMuon(*iMuon, thePrimaryVertex)) continue;
766  }
767 
768  TrackRef Track = iMuon->track();
769 
770  if (Track.isNonnull()) {
771  commonME_->hMuonTrackP_->Fill(Track->p());
772  commonME_->hMuonTrackPt_->Fill(Track->pt());
773  commonME_->hMuonTrackEta_->Fill(Track->eta());
774  commonME_->hMuonTrackPhi_->Fill(Track->phi());
775 
776  //ip histograms
777  commonME_->hMuonTrackDxy_->Fill(Track->dxy());
778  commonME_->hMuonTrackDz_->Fill(Track->dz());
779  }
780 
781  if (iMuon->isGlobalMuon()) {
782  Track = iMuon->combinedMuon();
783  glbNTrackerHits = countTrackerHits(*Track);
784  glbNMuonHits = countMuonHits(*Track);
785  } else if (iMuon->isTrackerMuon()) {
786  Track = iMuon->track();
787  trkNTrackerHits = countTrackerHits(*Track);
788  } else {
789  Track = iMuon->standAloneMuon();
790  }
791 
792  NTrackerHits = countTrackerHits(*Track);
793  muonME_->hNTrackerHits_->Fill(NTrackerHits);
794  muonME_->hNTrackerHits_vs_Pt_->Fill(Track->pt(), NTrackerHits);
795  muonME_->hNTrackerHits_vs_Eta_->Fill(Track->eta(), NTrackerHits);
796 
797  NMuonHits = countMuonHits(*Track);
798  muonME_->hNMuonHits_->Fill(NMuonHits);
799  muonME_->hNMuonHits_vs_Pt_->Fill(Track->pt(), NMuonHits);
800  muonME_->hNMuonHits_vs_Eta_->Fill(Track->eta(), NMuonHits);
801 
802 //list of histos for each type
803 
804 // muonME_->hNTrks_->Fill();
805  muonME_->hNTrksEta_->Fill(Track->eta());
806  muonME_->hNTrksPt_->Fill(Track->pt());
807 
808  commonME_->hMuonP_->Fill(muonP);
809  commonME_->hMuonPt_->Fill(muonPt);
810  commonME_->hMuonEta_->Fill(muonEta);
811  commonME_->hMuonPhi_->Fill(muonPhi);
812 
813  if (iMuon->isGlobalMuon()) {
814  double gtHitPat = iMuon->globalTrack()->hitPattern().numberOfHits(HitPattern::TRACK_HITS)
815  - iMuon->globalTrack()->hitPattern().numberOfValidHits();
816 
817  double itHitPat = iMuon->innerTrack()->hitPattern().numberOfHits(HitPattern::TRACK_HITS)
818  - iMuon->innerTrack()->hitPattern().numberOfValidHits();
819 
820  double otHitPat = iMuon->outerTrack()->hitPattern().numberOfHits(HitPattern::TRACK_HITS)
821  - iMuon->outerTrack()->hitPattern().numberOfValidHits();
822 
823  commonME_->hNInvalidHitsGTHitPattern_->Fill(gtHitPat);
824  commonME_->hNInvalidHitsITHitPattern_->Fill(itHitPat);
825  commonME_->hNInvalidHitsOTHitPattern_->Fill(otHitPat);
826  commonME_->hNDeltaInvalidHitsHitPattern_->Fill(gtHitPat - itHitPat - otHitPat);
827 
828  //must be global and standalone
829  if (iMuon->isStandAloneMuon()) {
830  commonME_->hStaToGlbDiffNMuonHitsEta_->Fill(Track->eta(),staNMuonHits-glbNMuonHits);
831  commonME_->hStaToGlbDiffNMuonHitsPt_->Fill(Track->pt(),staNMuonHits-glbNMuonHits);
832  commonME_->hStaToGlbDiffNMuonHits_->Fill(staNMuonHits-glbNMuonHits);
833  }
834 
835  //must be global and tracker
836  if (iMuon->isTrackerMuon()) {
837  commonME_->hTrkToGlbDiffNTrackerHitsEta_->Fill(Track->eta(),trkNTrackerHits-glbNTrackerHits);
838  commonME_->hTrkToGlbDiffNTrackerHitsPt_->Fill(Track->pt(),trkNTrackerHits-glbNTrackerHits);
839  commonME_->hTrkToGlbDiffNTrackerHits_->Fill(trkNTrackerHits-glbNTrackerHits);
840  }
841  }
842 
843  }//end of reco muon loop
844 
845 
846  // Associate by hits
847  for(TrackingParticleCollection::size_type i=0; i<nSim; i++) {
848  TrackingParticleRef simRef(simHandle, i);
849  const TrackingParticle* simTP = simRef.get();
850  if ( ! tpSelector_(*simTP) ) continue;
851 
852  //denominators for efficiency plots
853  const double simP = simRef->p();
854  const double simPt = simRef->pt();
855  const double simEta = doAbsEta_ ? fabs(simRef->eta()) : simRef->eta();
856  const double simPhi = simRef->phi();
857 
858  GlobalPoint simVtx(simRef->vertex().x(), simRef->vertex().y(), simRef->vertex().z());
859  GlobalVector simMom(simRef->momentum().x(), simRef->momentum().y(), simRef->momentum().z());
860  const double simDxy = -simVtx.x()*sin(simPhi)+simVtx.y()*cos(simPhi);
861  const double simDz = simVtx.z() - (simVtx.x()*simMom.x()+simVtx.y()*simMom.y())*simMom.z()/simMom.perp2();
862 
863  const unsigned int nSimHits = simRef->numberOfHits();
864 
865  muonME_->hSimP_ ->Fill(simP );
866  muonME_->hSimPt_ ->Fill(simPt );
867  muonME_->hSimEta_->Fill(simEta);
868  muonME_->hSimPhi_->Fill(simPhi);
869  muonME_->hSimDxy_->Fill(simDxy);
870  muonME_->hSimDz_->Fill(simDz);
871  muonME_->hNSimHits_->Fill(nSimHits);
872 
873  // Get sim-reco association for a simRef
874  vector<pair<RefToBase<Muon>, double> > MuRefV;
875  if ( simToMuonColl.find(simRef) != simToMuonColl.end() ) {
876  MuRefV = simToMuonColl[simRef];
877 
878  if ( !MuRefV.empty()) {
879  muonME_->hNSimToReco_->Fill(MuRefV.size());
880  const Muon* Mu = MuRefV.begin()->first.get();
881  if (!selector_(*Mu)) continue;
882  if (wantTightMuon_)
883  {
884  if (!muon::isTightMuon(*Mu, thePrimaryVertex)) continue;
885  }
886 
887  muonME_->fill(&*simTP, Mu);
888  }
889  }
890  }
891 }
892 
893 int
894 RecoMuonValidator::countMuonHits(const reco::Track& track) const {
895  TransientTrackingRecHit::ConstRecHitContainer result;
896 
897  int count = 0;
898 
899  for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
900  if((*hit)->isValid()) {
901  DetId recoid = (*hit)->geographicalId();
902  if ( recoid.det() == DetId::Muon ) count++;
903  }
904  }
905  return count;
906 }
907 
908 int
909 RecoMuonValidator::countTrackerHits(const reco::Track& track) const {
910  TransientTrackingRecHit::ConstRecHitContainer result;
911 
912  int count = 0;
913 
914  for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
915  if((*hit)->isValid()) {
916  DetId recoid = (*hit)->geographicalId();
917  if ( recoid.det() == DetId::Tracker ) count++;
918  }
919  }
920  return count;
921 }
922 /* vim:set ts=2 sts=2 sw=2 expandtab: */
MuonServiceProxy_cff.MuonServiceProxy
tuple MuonServiceProxy
Definition: MuonServiceProxy_cff.py:15
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
RecoMuonValidator::muonToken_
edm::EDGetTokenT< edm::View< reco::Muon > > muonToken_
Definition: RecoMuonValidator.h:53
RecoMuonValidator::simToken_
edm::EDGetTokenT< TrackingParticleCollection > simToken_
Definition: RecoMuonValidator.h:52
i
int i
Definition: DBlmapReader.cc:9
reco::MuonToSimCollection
std::map< edm::RefToBase< reco::Muon >, std::vector< std::pair< TrackingParticleRef, double > >, RefToBaseSort > MuonToSimCollection
Definition: MuonTrackType.h:34
reco::LeafCandidate::p
virtual double p() const
magnitude of momentum vector
Definition: LeafCandidate.h:108
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
TrackingParticleCollection
std::vector< TrackingParticle > TrackingParticleCollection
Definition: HiEvtPlaneFlatProducer.cc:102
DQMStore::IBooker::cd
void cd(void)
Definition: DQMStore.cc:266
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
RecoMuonValidator::simLabel_
edm::InputTag simLabel_
Definition: RecoMuonValidator.h:49
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.cc:182
edm::View::size
size_type size() const
RecoMuonValidator::trackType_
reco::MuonTrackType trackType_
Definition: RecoMuonValidator.h:79
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:107
RecoMuonValidator::primvertexToken_
edm::EDGetTokenT< reco::VertexCollection > primvertexToken_
Definition: RecoMuonValidator.h:61
RecoMuonValidator::countTrackerHits
virtual int countTrackerHits(const reco::Track &track) const
Definition: RecoMuonValidator.cc:909
reco::Muon::pfP4
reco::Candidate::LorentzVector pfP4() const
Definition: Muon.h:102
trigger::size_type
uint16_t size_type
Definition: TriggerTypeDefs.h:19
reco::LeafCandidate::eta
virtual double eta() const
momentum pseudorapidity
Definition: LeafCandidate.h:137
reco::LeafCandidate::pt
virtual double pt() const
transverse momentum
Definition: LeafCandidate.h:131
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
RecoMuonValidator::bookHistograms
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &) override
Definition: RecoMuonValidator.cc:548
MonitorElement::Fill
void Fill(long long x)
Definition: MonitorElement.h:164
reco::SimToMuonCollection
std::map< TrackingParticleRef, std::vector< std::pair< edm::RefToBase< reco::Muon >, double > > > SimToMuonCollection
Definition: MuonTrackType.h:35
FTS
FreeTrajectoryState FTS
Definition: RecoMuonValidator.cc:24
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:92
RecoMuonValidator::analyze
virtual void analyze(const edm::Event &event, const edm::EventSetup &eventSetup)
Definition: RecoMuonValidator.cc:630
edm::View::begin
const_iterator begin() const
RecoMuonValidator::theMuonService
MuonServiceProxy * theMuonService
Definition: RecoMuonValidator.h:68
TrackerOfflineValidation_Dqm_cff.dirName
tuple dirName
Definition: TrackerOfflineValidation_Dqm_cff.py:56
RecoMuonValidator::muAssocLabel_
edm::InputTag muAssocLabel_
Definition: RecoMuonValidator.h:55
query.result
tuple result
Definition: query.py:137
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:894
reco::LeafCandidate::charge
virtual int charge() const
electric charge
Definition: LeafCandidate.h:91
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
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.cc:132
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:242
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
event
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
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
RecoMuonValidator::dqmBeginRun
virtual void dqmBeginRun(const edm::Run &, const edm::EventSetup &eventSetup)
Definition: RecoMuonValidator.cc:614
DetId
Definition: DetId.h:18
TSOS
TrajectoryStateOnSurface TSOS
Definition: TestHits.cc:19
DQMStore::IBooker::setCurrentFolder
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:274
Mu
Definition: L1GtObject.h:30
edm::Handle::product
T const * product() const
Definition: Handle.h:81
DQMStore::IBooker::book2D
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:133
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.cc:197
RecoMuonValidator::subsystemname_
std::string subsystemname_
Definition: RecoMuonValidator.h:65
TrackingParticle::vertex
Point vertex() const
Parent vertex position.
Definition: TrackingParticle.h:120
edm::View::const_iterator
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:81
position
static int position[264][3]
Definition: ReadPGInfo.cc:509
TrackingParticle::momentum
Vector momentum() const
spatial momentum vector
Definition: TrackingParticle.cc:122
edm::RefToBaseVector::push_back
void push_back(const RefToBase< T > &)
Definition: RefToBaseVector.h:258
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&lt;C, T&gt; to the RefVector.
Definition: RefVector.h:64
RecoMuonValidator::endRun
virtual void endRun()
Definition: RecoMuonValidator.cc:622
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:745
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
reco::LeafCandidate::phi
virtual double phi() const
momentum azimuthal angle
Definition: LeafCandidate.h:133
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.cc:187
trackingRecHit_iterator
TrackingRecHitCollection::base::const_iterator trackingRecHit_iterator
iterator over a vector of reference to TrackingRecHit in the same collection
Definition: TrackingRecHitFwd.h:19
edm::Run
Definition: Run.h:41
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.5