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  Muons = muonHandle->refVector();
684 
685  edm::RefVector<TrackingParticleCollection> allTPs;
686  for (size_t i = 0; i < nSim; ++i) {
687  allTPs.push_back(TrackingParticleRef(simHandle,i));
688  }
689 
690 
691  muonME_->hNSim_->Fill(nSim);
692  muonME_->hNMuon_->Fill(muonColl.size());
693 
694  MuonAssociatorByHits::MuonToSimCollection muonToSimColl;
695  MuonAssociatorByHits::SimToMuonCollection simToMuonColl;
696 
697 
698  if ( doAssoc_ ) {
699  assoByHits->associateMuons(muonToSimColl, simToMuonColl, Muons, trackType_, allTPs, &event, &eventSetup);
700  } else {
701 
702 /*
703  // SimToMuon associations
704  Handle<reco::RecoToSimCollection> simToTrkMuHandle;
705  event.getByLabel(trkMuAssocLabel_, simToTrkMuHandle);
706  trkSimRecColl = *(simToTrkMuHandle.product());
707 
708  Handle<reco::RecoToSimCollection> simToStaMuHandle;
709  event.getByLabel(staMuAssocLabel_, simToStaMuHandle);
710  staSimRecColl = *(simToStaMuHandle.product());
711 
712  Handle<reco::RecoToSimCollection> simToGlbMuHandle;
713  event.getByLabel(glbMuAssocLabel_, simToGlbMuHandle);
714  glbSimRecColl = *(simToGlbMuHandle.product());
715 
716  // MuonToSim associations
717  Handle<reco::SimToRecoCollection> trkMuToSimHandle;
718  event.getByLabel(trkMuAssocLabel_, trkMuToSimHandle);
719  trkRecSimColl = *(trkMuToSimHandle.product());
720 
721  Handle<reco::SimToRecoCollection> staMuToSimHandle;
722  event.getByLabel(staMuAssocLabel_, staMuToSimHandle);
723  staRecSimColl = *(staMuToSimHandle.product());
724 
725  Handle<reco::SimToRecoCollection> glbMuToSimHandle;
726  event.getByLabel(glbMuAssocLabel_, glbMuToSimHandle);
727  glbRecSimColl = *(glbMuToSimHandle.product());
728 */
729  }
730 
731 
732 
733  int glbNTrackerHits = 0; int trkNTrackerHits = 0;
734  int glbNMuonHits = 0; int staNMuonHits = 0;
735  int NTrackerHits = 0; int NMuonHits = 0;
736 
737 
738  // Analyzer reco::Muon
739  for(View<Muon>::const_iterator iMuon = muonColl.begin();
740  iMuon != muonColl.end(); ++iMuon) {
741 
742  double muonP, muonPt, muonEta, muonPhi;
743  if (usePFMuon_) {
744  muonP = iMuon->pfP4().P();
745  muonPt = iMuon->pfP4().Pt();
746  muonEta = iMuon->pfP4().Eta();
747  muonPhi = iMuon->pfP4().Phi();
748  }
749  else {
750  muonP = iMuon->p();
751  muonPt = iMuon->pt();
752  muonEta = iMuon->eta();
753  muonPhi = iMuon->phi();
754  }
755 
756  //histograms for fractions
757  commonME_->hMuonAllP_->Fill(muonP);
758  commonME_->hMuonAllPt_->Fill(muonPt);
759  commonME_->hMuonAllEta_->Fill(muonEta);
760  commonME_->hMuonAllPhi_->Fill(muonPhi);
761 
762  if (!selector_(*iMuon)) continue;
763  if (wantTightMuon_)
764  {
765  if (!muon::isTightMuon(*iMuon, thePrimaryVertex)) continue;
766  }
767 
768  TrackRef Track = iMuon->track();
769 
770  if (Track.isNonnull()) {
771  commonME_->hMuonTrackP_->Fill(Track->p());
772  commonME_->hMuonTrackPt_->Fill(Track->pt());
773  commonME_->hMuonTrackEta_->Fill(Track->eta());
774  commonME_->hMuonTrackPhi_->Fill(Track->phi());
775 
776  //ip histograms
777  commonME_->hMuonTrackDxy_->Fill(Track->dxy());
778  commonME_->hMuonTrackDz_->Fill(Track->dz());
779  }
780 
781  if (iMuon->isGlobalMuon()) {
782  Track = iMuon->combinedMuon();
783  glbNTrackerHits = countTrackerHits(*Track);
784  glbNMuonHits = countMuonHits(*Track);
785  } else if (iMuon->isTrackerMuon()) {
786  Track = iMuon->track();
787  trkNTrackerHits = countTrackerHits(*Track);
788  } else {
789  Track = iMuon->standAloneMuon();
790  }
791 
792  NTrackerHits = countTrackerHits(*Track);
793  muonME_->hNTrackerHits_->Fill(NTrackerHits);
794  muonME_->hNTrackerHits_vs_Pt_->Fill(Track->pt(), NTrackerHits);
795  muonME_->hNTrackerHits_vs_Eta_->Fill(Track->eta(), NTrackerHits);
796 
797  NMuonHits = countMuonHits(*Track);
798  muonME_->hNMuonHits_->Fill(NMuonHits);
799  muonME_->hNMuonHits_vs_Pt_->Fill(Track->pt(), NMuonHits);
800  muonME_->hNMuonHits_vs_Eta_->Fill(Track->eta(), NMuonHits);
801 
802 //list of histos for each type
803 
804 // muonME_->hNTrks_->Fill();
805  muonME_->hNTrksEta_->Fill(Track->eta());
806  muonME_->hNTrksPt_->Fill(Track->pt());
807 
808  commonME_->hMuonP_->Fill(muonP);
809  commonME_->hMuonPt_->Fill(muonPt);
810  commonME_->hMuonEta_->Fill(muonEta);
811  commonME_->hMuonPhi_->Fill(muonPhi);
812 
813  if (iMuon->isGlobalMuon()) {
814  double gtHitPat = iMuon->globalTrack()->hitPattern().numberOfHits() - iMuon->globalTrack()->hitPattern().numberOfValidHits();
815  double itHitPat = iMuon->innerTrack()->hitPattern().numberOfHits() - iMuon->innerTrack()->hitPattern().numberOfValidHits();
816  double otHitPat = iMuon->outerTrack()->hitPattern().numberOfHits() - iMuon->outerTrack()->hitPattern().numberOfValidHits();
817 
818  commonME_->hNInvalidHitsGTHitPattern_->Fill(gtHitPat);
819  commonME_->hNInvalidHitsITHitPattern_->Fill(itHitPat);
820  commonME_->hNInvalidHitsOTHitPattern_->Fill(otHitPat);
821  commonME_->hNDeltaInvalidHitsHitPattern_->Fill(gtHitPat - itHitPat - otHitPat);
822 
823  //must be global and standalone
824  if (iMuon->isStandAloneMuon()) {
825  commonME_->hStaToGlbDiffNMuonHitsEta_->Fill(Track->eta(),staNMuonHits-glbNMuonHits);
826  commonME_->hStaToGlbDiffNMuonHitsPt_->Fill(Track->pt(),staNMuonHits-glbNMuonHits);
827  commonME_->hStaToGlbDiffNMuonHits_->Fill(staNMuonHits-glbNMuonHits);
828  }
829 
830  //must be global and tracker
831  if (iMuon->isTrackerMuon()) {
832  commonME_->hTrkToGlbDiffNTrackerHitsEta_->Fill(Track->eta(),trkNTrackerHits-glbNTrackerHits);
833  commonME_->hTrkToGlbDiffNTrackerHitsPt_->Fill(Track->pt(),trkNTrackerHits-glbNTrackerHits);
834  commonME_->hTrkToGlbDiffNTrackerHits_->Fill(trkNTrackerHits-glbNTrackerHits);
835  }
836  }
837 
838  }//end of reco muon loop
839 
840 
841  // Associate by hits
842  for(TrackingParticleCollection::size_type i=0; i<nSim; i++) {
843  TrackingParticleRef simRef(simHandle, i);
844  const TrackingParticle* simTP = simRef.get();
845  if ( ! tpSelector_(*simTP) ) continue;
846 
847  //denominators for efficiency plots
848  const double simP = simRef->p();
849  const double simPt = simRef->pt();
850  const double simEta = doAbsEta_ ? fabs(simRef->eta()) : simRef->eta();
851  const double simPhi = simRef->phi();
852 
853  GlobalPoint simVtx(simRef->vertex().x(), simRef->vertex().y(), simRef->vertex().z());
854  GlobalVector simMom(simRef->momentum().x(), simRef->momentum().y(), simRef->momentum().z());
855  const double simDxy = -simVtx.x()*sin(simPhi)+simVtx.y()*cos(simPhi);
856  const double simDz = simVtx.z() - (simVtx.x()*simMom.x()+simVtx.y()*simMom.y())*simMom.z()/simMom.perp2();
857 
858  const unsigned int nSimHits = simRef->numberOfHits();
859 
860  muonME_->hSimP_ ->Fill(simP );
861  muonME_->hSimPt_ ->Fill(simPt );
862  muonME_->hSimEta_->Fill(simEta);
863  muonME_->hSimPhi_->Fill(simPhi);
864  muonME_->hSimDxy_->Fill(simDxy);
865  muonME_->hSimDz_->Fill(simDz);
866  muonME_->hNSimHits_->Fill(nSimHits);
867 
868  // Get sim-reco association for a simRef
869  vector<pair<RefToBase<Muon>, double> > MuRefV;
870  if ( simToMuonColl.find(simRef) != simToMuonColl.end() ) {
871  MuRefV = simToMuonColl[simRef];
872 
873  if ( !MuRefV.empty()) {
874  muonME_->hNSimToReco_->Fill(MuRefV.size());
875  const Muon* Mu = MuRefV.begin()->first.get();
876  if (!selector_(*Mu)) continue;
877  if (wantTightMuon_)
878  {
879  if (!muon::isTightMuon(*Mu, thePrimaryVertex)) continue;
880  }
881 
882  muonME_->fill(&*simTP, Mu);
883  }
884  }
885  }
886 }
887 
888 int
889 RecoMuonValidator::countMuonHits(const reco::Track& track) const {
890  TransientTrackingRecHit::ConstRecHitContainer result;
891 
892  int count = 0;
893 
894  for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
895  if((*hit)->isValid()) {
896  DetId recoid = (*hit)->geographicalId();
897  if ( recoid.det() == DetId::Muon ) count++;
898  }
899  }
900  return count;
901 }
902 
903 int
904 RecoMuonValidator::countTrackerHits(const reco::Track& track) const {
905  TransientTrackingRecHit::ConstRecHitContainer result;
906 
907  int count = 0;
908 
909  for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
910  if((*hit)->isValid()) {
911  DetId recoid = (*hit)->geographicalId();
912  if ( recoid.det() == DetId::Tracker ) count++;
913  }
914  }
915  return count;
916 }
917 /* 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::View::const_iterator
boost::indirect_iterator< typename seq_t::const_iterator > const_iterator
Definition: View.h:81
reco::LeafCandidate::pt
virtual float pt() const
transverse momentum
Definition: LeafCandidate.h:178
TrackingParticleCollection
std::vector< TrackingParticle > TrackingParticleCollection
Definition: HiEvtPlaneFlatProducer.cc:115
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:1391
reco::LeafCandidate::phi
virtual float phi() const
momentum azimuthal angle
Definition: LeafCandidate.h:180
TrackingParticle::vertex
Point vertex() const
Parent vertex position.
Definition: TrackingParticle.cc:221
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:191
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:904
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:166
FTS
FreeTrajectoryState FTS
Definition: RecoMuonValidator.cc:27
RecoMuonValidator::tpSelector_
TrackingParticleSelector tpSelector_
Definition: RecoMuonValidator.h:77
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:250
reco::modules::TrackingParticleSelector
SingleObjectSelector< TrackingParticleCollection,::TrackingParticleSelector > TrackingParticleSelector
Definition: TrackingParticleSelector.cc:18
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
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:889
reco::LeafCandidate::charge
virtual int charge() const
electric charge
Definition: LeafCandidate.h:140
DQMStore::IBooker::book1D
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:113
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:141
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
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:63
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:40
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
DQMStore::IBooker::book2D
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:131
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:62
reco::BeamSpot::BeamWidthY
double BeamWidthY() const
beam width Y
Definition: BeamSpot.h:88
edm::View::size
size_type size() const
TrackingParticle::eta
double eta() const
Momentum pseudorapidity. Note this is taken from the first SimTrack only.
Definition: TrackingParticle.cc:206
edm::Handle::product
T const * product() const
Definition: Handle.h:81
RecoMuonValidator::subsystemname_
std::string subsystemname_
Definition: RecoMuonValidator.h:66
edm::InputTag::label
std::string const & label() const
Definition: InputTag.h:42
DQMStore::save
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", const uint32_t run=0, const uint32_t lumi=0, SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE", const bool resetMEsAfterWriting=false)
Definition: DQMStore.cc:2540
position
static int position[264][3]
Definition: ReadPGInfo.cc:509
TrackingParticle::momentum
Vector momentum() const
spatial momentum vector
Definition: TrackingParticle.cc:131
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
edm::View::begin
const_iterator begin() const
TrackingParticle::charge
int charge() const
Electric charge. Note this is taken from the first SimTrack only.
Definition: TrackingParticle.cc:116
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
edm::Ref::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:242
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:196
edm::Run
Definition: Run.h:41
reco::Track::recHitsEnd
trackingRecHit_iterator recHitsEnd() const
Iterator to last hit on the track.
Definition: Track.h:65
Generated for CMSSW Reference Manual by doxygen 1.8.5