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