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>("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()
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: */
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
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:252
TrackingParticleCollection
std::vector< TrackingParticle > TrackingParticleCollection
Definition: HiEvtPlaneFlatProducer.cc:111
TrackingParticle::momentum
Vector momentum() const
spatial momentum vector
Definition: TrackingParticle.h:109
DQMStore::IBooker::cd
void cd(void)
Definition: DQMStore.cc:268
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
TrackValidation_cff.pset
tuple pset
Definition: TrackValidation_cff.py:374
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::LeafCandidate::phi
virtual double phi() const final
momentum azimuthal angle
Definition: LeafCandidate.h:133
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
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:910
reco::Muon::pfP4
reco::Candidate::LorentzVector pfP4() const
Definition: Muon.h:102
trigger::size_type
uint16_t size_type
Definition: TriggerTypeDefs.h:19
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:549
mps_fire.result
tuple result
Definition: mps_fire.py:83
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
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:99
RecoMuonValidator::analyze
virtual void analyze(const edm::Event &event, const edm::EventSetup &eventSetup)
Definition: RecoMuonValidator.cc:631
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
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
reco::LeafCandidate::charge
virtual int charge() const final
electric charge
Definition: LeafCandidate.h:91
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
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
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:615
DetId
Definition: DetId.h:18
TSOS
TrajectoryStateOnSurface TSOS
Definition: TestHits.cc:19
DQMStore::IBooker::setCurrentFolder
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:276
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.h:184
reco::LeafCandidate::p
virtual double p() const final
magnitude of momentum vector
Definition: LeafCandidate.h:108
RecoMuonValidator::subsystemname_
std::string subsystemname_
Definition: RecoMuonValidator.h:65
TrackingParticle::vertex
Point vertex() const
Parent vertex position.
Definition: TrackingParticle.h:199
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
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&lt;C, T&gt; to the RefVector.
Definition: RefVector.h:62
RecoMuonValidator::endRun
virtual void endRun()
Definition: RecoMuonValidator.cc:623
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
reco::LeafCandidate::eta
virtual double eta() const final
momentum pseudorapidity
Definition: LeafCandidate.h:137
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
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
edm::Run
Definition: Run.h:43
reco::LeafCandidate::pt
virtual double pt() const final
transverse momentum
Definition: LeafCandidate.h:131
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