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