32 MEP hSimP_,
hSimPt_, hSimEta_, hSimPhi_, hSimDxy_, hSimDz_;
69 MEP hPullPt_, hPullEta_, hPullPhi_,
hPullQPt_, hPullDxy_, hPullDz_;
105 hNSim_ = ibooker.
book1D(
"NSim",
"Number of particles per event", hDim.
nTrks, -0.5, hDim.
nTrks + 0.5);
106 hNMuon_ = ibooker.
book1D(
"NMuon",
"Number of muons per event", hDim.
nTrks, -0.5, hDim.
nTrks + 0.5);
109 hNTrks_ = ibooker.
book1D(
"NTrks",
"Number of reco tracks per event", hDim.
nTrks, -0.5, hDim.
nTrks + 0.5);
111 hNTrksPt_ = ibooker.
book1D(
"NTrksPt",
"Number of reco tracks vs p_{T}", hDim.
nBinPt, hDim.
minPt, hDim.
maxPt);
136 hPFMomAssCorrectness =
137 ibooker.
book1D(
"hPFMomAssCorrectness",
"Corrected momentum assignement PF/RECO", 2, 0.5, 2.5);
138 hPt_vs_PFMomAssCorrectness = ibooker.
book2D(
"hPt_vs_PFMomAssCorrectness",
139 "Corrected momentum assignement PF/RECO",
147 hdPt_vs_Pt_ = ibooker.
book2D(
"dPt_vs_Pt",
148 "#Delta(p_{T}) vs p_{T}",
155 hdPt_vs_Eta_ = ibooker.
book2D(
"dPt_vs_Eta",
156 "#Delta(p_{T}) vs #eta",
166 hErrP_vs_Eta_ = ibooker.
book2D(
"ErrP_vs_Eta",
167 "#Delta(p)/p vs #eta",
174 hErrPt_vs_Eta_ = ibooker.
book2D(
"ErrPt_vs_Eta",
175 "#Delta(p_{T})/p_{T} vs #eta",
182 hErrQPt_vs_Eta_ = ibooker.
book2D(
"ErrQPt_vs_Eta",
183 "#Delta(q/p_{T})/(q/p_{T}) vs #eta",
190 hErrEta_vs_Eta_ = ibooker.
book2D(
"ErrEta_vs_Eta",
191 "#sigma(#eta) vs #eta",
200 hErrP_vs_P_ = ibooker.
book2D(
202 hErrPt_vs_Pt_ = ibooker.
book2D(
"ErrPt_vs_Pt",
203 "#Delta(p_{T})/p_{T} vs p_{T}",
210 hErrQPt_vs_Pt_ = ibooker.
book2D(
"ErrQPt_vs_Pt",
211 "#Delta(q/p_{T})/(q/p_{T}) vs p_{T}",
228 hPullPt_vs_Eta_ = ibooker.
book2D(
"PullPt_vs_Eta",
229 "Pull(p_{T}) vs #eta",
236 hPullEta_vs_Eta_ = ibooker.
book2D(
"PullEta_vs_Eta",
237 "Pull(#eta) vs #eta",
244 hPullPhi_vs_Eta_ = ibooker.
book2D(
"PullPhi_vs_Eta",
245 "Pull(#phi) vs #eta",
254 hPullPt_vs_Pt_ = ibooker.
book2D(
"PullPt_vs_Pt",
255 "Pull(p_{T}) vs p_{T}",
262 hPullEta_vs_Pt_ = ibooker.
book2D(
"PullEta_vs_Pt",
263 "Pull(#eta) vs p_{T}",
272 const int nHits = 100;
273 hNHits_ = ibooker.
book1D(
"NHits",
"Number of hits", nHits + 1, -0.5, nHits + 0.5);
274 hNHits_vs_Pt_ = ibooker.
book2D(
"NHits_vs_Pt",
275 "Number of hits vs p_{T}",
282 hNHits_vs_Eta_ = ibooker.
book2D(
"NHits_vs_Eta",
283 "Number of hits vs #eta",
290 hNSimHits_ = ibooker.
book1D(
"NSimHits",
"Number of simHits", nHits + 1, -0.5, nHits + 0.5);
292 const int nLostHits = 5;
293 hNLostHits_ = ibooker.
book1D(
"NLostHits",
"Number of Lost hits", nLostHits + 1, -0.5, nLostHits + 0.5);
294 hNLostHits_vs_Pt_ = ibooker.
book2D(
"NLostHits_vs_Pt",
295 "Number of lost Hits vs p_{T}",
302 hNLostHits_vs_Eta_ = ibooker.
book2D(
"NLostHits_vs_Eta",
303 "Number of lost Hits vs #eta",
311 const int nTrackerHits = 40;
313 ibooker.
book1D(
"NTrackerHits",
"Number of valid tracker hits", nTrackerHits + 1, -0.5, nTrackerHits + 0.5);
314 hNTrackerHits_vs_Pt_ = ibooker.
book2D(
"NTrackerHits_vs_Pt",
315 "Number of valid traker hits vs p_{T}",
319 nTrackerHits / 4 + 1,
321 nTrackerHits + 0.25);
322 hNTrackerHits_vs_Eta_ = ibooker.
book2D(
"NTrackerHits_vs_Eta",
323 "Number of valid tracker hits vs #eta",
327 nTrackerHits / 4 + 1,
329 nTrackerHits + 0.25);
331 const int nMuonHits = 60;
332 hNMuonHits_ = ibooker.
book1D(
"NMuonHits",
"Number of valid muon hits", nMuonHits + 1, -0.5, nMuonHits + 0.5);
333 hNMuonHits_vs_Pt_ = ibooker.
book2D(
"NMuonHits_vs_Pt",
334 "Number of valid muon hits vs p_{T}",
341 hNMuonHits_vs_Eta_ = ibooker.
book2D(
"NMuonHits_vs_Eta",
342 "Number of valid muon hits vs #eta",
350 hNDof_ = ibooker.
book1D(
"NDof",
"Number of DoF", hDim.
nDof + 1, -0.5, hDim.
nDof + 0.5);
351 hChi2_ = ibooker.
book1D(
"Chi2",
"#Chi^{2}", hDim.
nBinErr, 0, 200);
352 hChi2Norm_ = ibooker.
book1D(
"Chi2Norm",
"Normalized #Chi^{2}", hDim.
nBinErr, 0, 50);
353 hChi2Prob_ = ibooker.
book1D(
"Chi2Prob",
"Prob(#Chi^{2})", hDim.
nBinErr, 0, 1);
355 hNDof_vs_Eta_ = ibooker.
book2D(
"NDof_vs_Eta",
356 "Number of DoF vs #eta",
363 hChi2_vs_Eta_ = ibooker.
book2D(
365 hChi2Norm_vs_Eta_ = ibooker.
book2D(
367 hChi2Prob_vs_Eta_ = ibooker.
book2D(
371 ibooker.
book1D(
"NSimToReco",
"Number of associated reco tracks", hDim.
nAssoc + 1, -0.5, hDim.
nAssoc + 0.5);
373 ibooker.
book1D(
"NRecoToSim",
"Number of associated sim TP's", hDim.
nAssoc + 1, -0.5, hDim.
nAssoc + 0.5);
380 const double simP = simRef->
p();
381 const double simPt = simRef->
pt();
382 const double simEta = doAbsEta_ ? fabs(simRef->
eta()) : simRef->
eta();
383 const double simPhi = simRef->
phi();
384 const double simQ = simRef->
charge();
385 const double simQPt = simQ / simPt;
389 const double simDxy = -simVtx.
x() *
sin(simPhi) + simVtx.y() *
cos(simPhi);
390 const double simDz = simVtx.z() - (simVtx.x() * simMom.x() + simVtx.y() * simMom.y()) * simMom.z() / simMom.perp2();
392 const double recoQ = muonRef->
charge();
393 if (simQ * recoQ < 0) {
394 hMisQPt_->
Fill(simPt);
395 hMisQEta_->
Fill(simEta);
398 double recoP, recoPt, recoEta, recoPhi, recoQPt;
401 const double origRecoPt = muonRef->
pt();
404 const double origRecoQPt = recoQ / origRecoPt;
405 recoP = muonRef->
pfP4().P();
406 recoPt = muonRef->
pfP4().Pt();
407 recoEta = muonRef->
pfP4().Eta();
408 recoPhi = muonRef->
pfP4().Phi();
409 recoQPt = recoQ / recoPt;
410 hErrPt_PF_->
Fill((recoPt - origRecoPt) / origRecoPt);
411 hErrQPt_PF_->
Fill((recoQPt - origRecoQPt) / origRecoQPt);
413 hdPt_vs_Eta_->
Fill(recoEta, recoPt - origRecoPt);
414 hdPt_vs_Pt_->
Fill(recoPt, recoPt - origRecoPt);
416 int theCorrectPFAss = (fabs(recoPt - simPt) < fabs(origRecoPt - simPt)) ? 1 : 2;
417 hPFMomAssCorrectness->
Fill(theCorrectPFAss);
418 hPt_vs_PFMomAssCorrectness->
Fill(simPt, theCorrectPFAss);
422 recoP = muonRef->
p();
423 recoPt = muonRef->
pt();
424 recoEta = muonRef->
eta();
425 recoPhi = muonRef->
phi();
426 recoQPt = recoQ / recoPt;
429 const double errP = (recoP - simP) / simP;
430 const double errPt = (recoPt - simPt) / simPt;
431 const double errEta = (recoEta - simEta) / simEta;
432 const double errPhi = (recoPhi - simPhi) / simPhi;
433 const double errQPt = (recoQPt - simQPt) / simQPt;
441 hErrPt_->
Fill(errPt);
442 hErrEta_->
Fill(errEta);
443 hErrPhi_->
Fill(errPhi);
445 if (fabs(simEta) > 0. && fabs(simEta) < 0.8) {
446 hErrPBarrel_->
Fill(errP);
447 hErrPtBarrel_->
Fill(errPt);
448 }
else if (fabs(simEta) > 0.8 && fabs(simEta) < 1.2) {
449 hErrPOverlap_->
Fill(errP);
450 hErrPtOverlap_->
Fill(errPt);
451 }
else if (fabs(simEta) > 1.2) {
452 hErrPEndcap_->
Fill(errP);
453 hErrPtEndcap_->
Fill(errPt);
456 hErrP_vs_Eta_->
Fill(simEta, errP);
457 hErrPt_vs_Eta_->
Fill(simEta, errPt);
458 hErrQPt_vs_Eta_->
Fill(simEta, errQPt);
460 hErrP_vs_P_->
Fill(simP, errP);
461 hErrPt_vs_Pt_->
Fill(simPt, errPt);
462 hErrQPt_vs_Pt_->
Fill(simQPt, errQPt);
464 hErrEta_vs_Eta_->
Fill(simEta, errEta);
470 const int nRecoHits = recoRef->numberOfValidHits();
471 const int nLostHits = recoRef->numberOfLostHits();
473 hNHits_->
Fill(nRecoHits);
474 hNHits_vs_Pt_->
Fill(simPt, nRecoHits);
475 hNHits_vs_Eta_->
Fill(simEta, nRecoHits);
477 hNLostHits_->
Fill(nLostHits);
478 hNLostHits_vs_Pt_->
Fill(simPt, nLostHits);
479 hNLostHits_vs_Eta_->
Fill(simEta, nLostHits);
481 const double recoNDof = recoRef->ndof();
482 const double recoChi2 = recoRef->chi2();
483 const double recoChi2Norm = recoRef->normalizedChi2();
484 const double recoChi2Prob = TMath::Prob(recoRef->chi2(), static_cast<int>(recoRef->ndof()));
486 hNDof_->
Fill(recoNDof);
487 hChi2_->
Fill(recoChi2);
488 hChi2Norm_->
Fill(recoChi2Norm);
489 hChi2Prob_->
Fill(recoChi2Prob);
491 hNDof_vs_Eta_->
Fill(simEta, recoNDof);
492 hChi2_vs_Eta_->
Fill(simEta, recoChi2);
493 hChi2Norm_vs_Eta_->
Fill(simEta, recoChi2Norm);
494 hChi2Prob_vs_Eta_->
Fill(simEta, recoChi2Prob);
496 const double recoDxy = recoRef->dxy();
497 const double recoDz = recoRef->dz();
499 const double errDxy = (recoDxy - simDxy) / simDxy;
500 const double errDz = (recoDz - simDz) / simDz;
501 hErrDxy_->
Fill(errDxy);
502 hErrDz_->
Fill(errDz);
504 const double pullPt = (recoPt - simPt) / recoRef->ptError();
505 const double pullQPt = (recoQPt - simQPt) / recoRef->qoverpError();
506 const double pullEta = (recoEta - simEta) / recoRef->etaError();
507 const double pullPhi = (recoPhi - simPhi) / recoRef->phiError();
508 const double pullDxy = (recoDxy - simDxy) / recoRef->dxyError();
509 const double pullDz = (recoDz - simDz) / recoRef->dzError();
511 hPullPt_->
Fill(pullPt);
512 hPullEta_->
Fill(pullEta);
513 hPullPhi_->
Fill(pullPhi);
514 hPullQPt_->
Fill(pullQPt);
515 hPullDxy_->
Fill(pullDxy);
516 hPullDz_->
Fill(pullDz);
518 hPullPt_vs_Eta_->
Fill(simEta, pullPt);
519 hPullPt_vs_Pt_->
Fill(simPt, pullPt);
521 hPullEta_vs_Eta_->
Fill(simEta, pullEta);
522 hPullPhi_vs_Eta_->
Fill(simEta, pullPhi);
524 hPullEta_vs_Pt_->
Fill(simPt, pullEta);
547 MEP hMuonTrackP_,
hMuonTrackPt_, hMuonTrackEta_, hMuonTrackPhi_, hMuonTrackDxy_, hMuonTrackDz_;
698 const int nHits = 100;
702 "Difference of number of tracker hits (tkMuon - globalMuon)",
707 "Difference of number of muon hits (staMuon - globalMuon)",
713 ibooker.
book2D(
"TrkGlbDiffNTrackerHitsEta",
714 "Difference of number of tracker hits (tkMuon - globalMuon)",
722 "Difference of number of muon hits (staMuon - globalMuon)",
731 "Difference of number of tracker hits (tkMuon - globalMuon)",
739 "Difference of number of muon hits (staMuon - globalMuon)",
749 "NInvalidHitsGTHitPattern",
"Number of invalid hits on a global track", nHits + 1, -0.5, nHits + 0.5);
751 "NInvalidHitsITHitPattern",
"Number of invalid hits on an inner track", nHits + 1, -0.5, nHits + 0.5);
753 "NInvalidHitsOTHitPattern",
"Number of invalid hits on an outer track", nHits + 1, -0.5, nHits + 0.5);
755 ibooker.
book1D(
"hNDeltaInvalidHitsHitPattern",
756 "The discrepancy for Number of invalid hits on an global track and inner and outer tracks",
819 unsigned int theIndexOfThePrimaryVertex = 999.;
820 for (
unsigned int ind = 0; ind < recVtxs->size(); ++ind) {
821 if ((*recVtxs)[ind].isValid() && !((*recVtxs)[ind].isFake())) {
822 theIndexOfThePrimaryVertex = ind;
826 if (theIndexOfThePrimaryVertex < 100) {
827 posVtx = ((*recVtxs)[theIndexOfThePrimaryVertex]).
position();
828 errVtx = ((*recVtxs)[theIndexOfThePrimaryVertex]).
error();
830 LogInfo(
"RecoMuonValidator") <<
"reco::PrimaryVertex not found, use BeamSpot position instead\n";
834 posVtx =
bs.position();
835 errVtx(0, 0) =
bs.BeamWidthX();
836 errVtx(1, 1) =
bs.BeamWidthY();
837 errVtx(2, 2) =
bs.sigmaZ();
855 assoByHits = associatorBase.
product();
861 for (
size_t i = 0;
i < muonHandle->size(); ++
i) {
862 Muons.push_back(muonHandle->refAt(
i));
866 for (
size_t i = 0;
i < nSim; ++
i) {
908 int glbNTrackerHits = 0;
909 int trkNTrackerHits = 0;
910 int glbNMuonHits = 0;
911 int staNMuonHits = 0;
912 int NTrackerHits = 0;
919 muonP = iMuon->pfP4().P();
920 muonPt = iMuon->pfP4().Pt();
945 if (
Track.isNonnull()) {
956 if (iMuon->isGlobalMuon()) {
957 Track = iMuon->combinedMuon();
960 }
else if (iMuon->isTrackerMuon()) {
961 Track = iMuon->track();
988 if (iMuon->isGlobalMuon()) {
989 double gtHitPat = iMuon->globalTrack()->hitPattern().numberOfAllHits(HitPattern::TRACK_HITS) -
990 iMuon->globalTrack()->hitPattern().numberOfValidHits();
992 double itHitPat = iMuon->innerTrack()->hitPattern().numberOfAllHits(HitPattern::TRACK_HITS) -
993 iMuon->innerTrack()->hitPattern().numberOfValidHits();
995 double otHitPat = iMuon->outerTrack()->hitPattern().numberOfAllHits(HitPattern::TRACK_HITS) -
996 iMuon->outerTrack()->hitPattern().numberOfValidHits();
1004 if (iMuon->isStandAloneMuon()) {
1011 if (iMuon->isTrackerMuon()) {
1028 const double simP = simRef->p();
1029 const double simPt = simRef->pt();
1030 const double simEta =
doAbsEta_ ? fabs(simRef->eta()) : simRef->eta();
1031 const double simPhi = simRef->phi();
1033 GlobalPoint simVtx(simRef->vertex().x(), simRef->vertex().y(), simRef->vertex().z());
1034 GlobalVector simMom(simRef->momentum().x(), simRef->momentum().y(), simRef->momentum().z());
1035 const double simDxy = -simVtx.
x() *
sin(simPhi) + simVtx.y() *
cos(simPhi);
1036 const double simDz = simVtx.z() - (simVtx.x() * simMom.x() + simVtx.y() * simMom.y()) * simMom.z() / simMom.perp2();
1038 const unsigned int nSimHits = simRef->numberOfHits();
1049 vector<pair<RefToBase<Muon>,
double> > MuRefV;
1050 if (simToMuonColl.find(simRef) != simToMuonColl.end()) {
1051 MuRefV = simToMuonColl[simRef];
1053 if (!MuRefV.empty()) {
1055 const Muon*
Mu = MuRefV.begin()->first.get();
1075 if ((*hit)->isValid()) {
1076 DetId recoid = (*hit)->geographicalId();
1090 if ((*hit)->isValid()) {
1091 DetId recoid = (*hit)->geographicalId();