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