Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:33:11 

0001 #include "Validation/RecoMuon/src/RecoMuonValidator.h"
0002 
0003 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
0004 
0005 #include "FWCore/ServiceRegistry/interface/Service.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 //#include "DQMServices/Core/interface/DQMStore.h"
0008 
0009 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0010 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0011 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0012 
0013 // for selection cut
0014 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
0015 
0016 #include "TMath.h"
0017 
0018 using namespace std;
0019 using namespace edm;
0020 using namespace reco;
0021 
0022 typedef TrajectoryStateOnSurface TSOS;
0023 typedef FreeTrajectoryState FTS;
0024 
0025 //
0026 //Struct containing all histograms definitions
0027 //
0028 struct RecoMuonValidator::MuonME {
0029   typedef MonitorElement* MEP;
0030 
0031   //general kinematics
0032   MEP hSimP_, hSimPt_, hSimEta_, hSimPhi_, hSimDxy_, hSimDz_;
0033   //only for efficiencies
0034   MEP hP_, hPt_, hEta_, hPhi_;
0035   MEP hNSim_, hNMuon_;
0036 
0037   //misc vars
0038   MEP hNTrks_, hNTrksEta_, hNTrksPt_;
0039   MEP hMisQPt_, hMisQEta_;
0040 
0041   //resolutions
0042   MEP hErrP_, hErrPt_, hErrEta_, hErrPhi_;
0043   MEP hErrPBarrel_, hErrPOverlap_, hErrPEndcap_;
0044   MEP hErrPtBarrel_, hErrPtOverlap_, hErrPtEndcap_;
0045   MEP hErrDxy_, hErrDz_;
0046 
0047   MEP hErrP_vs_Eta_, hErrPt_vs_Eta_, hErrQPt_vs_Eta_;
0048   MEP hErrP_vs_P_, hErrPt_vs_Pt_, hErrQPt_vs_Pt_, hErrEta_vs_Eta_;
0049 
0050   //PF-RECO event-by-event comparisons
0051   MEP hErrPt_PF_;
0052   MEP hErrQPt_PF_;
0053   MEP hdPt_vs_Eta_;
0054   MEP hdPt_vs_Pt_;
0055   MEP hPFMomAssCorrectness;
0056   MEP hPt_vs_PFMomAssCorrectness;
0057 
0058   //hit pattern
0059   MEP hNSimHits_;
0060   MEP hNSimToReco_, hNRecoToSim_;
0061 
0062   MEP hNHits_, hNLostHits_, hNTrackerHits_, hNMuonHits_;
0063   MEP hNHits_vs_Pt_, hNHits_vs_Eta_;
0064   MEP hNLostHits_vs_Pt_, hNLostHits_vs_Eta_;
0065   MEP hNTrackerHits_vs_Pt_, hNTrackerHits_vs_Eta_;
0066   MEP hNMuonHits_vs_Pt_, hNMuonHits_vs_Eta_;
0067 
0068   //pulls
0069   MEP hPullPt_, hPullEta_, hPullPhi_, hPullQPt_, hPullDxy_, hPullDz_;
0070   MEP hPullPt_vs_Eta_, hPullPt_vs_Pt_, hPullEta_vs_Eta_, hPullPhi_vs_Eta_, hPullEta_vs_Pt_;
0071 
0072   //chi2, ndof
0073   MEP hNDof_, hChi2_, hChi2Norm_, hChi2Prob_;
0074   MEP hNDof_vs_Eta_, hChi2_vs_Eta_, hChi2Norm_vs_Eta_, hChi2Prob_vs_Eta_;
0075 
0076   bool doAbsEta_;
0077   bool usePFMuon_;
0078 
0079   //
0080   //books histograms
0081   //
0082   void bookHistos(DQMStore::IBooker& ibooker, const string& dirName, const HistoDimensions& hDim)
0083 
0084   {
0085     ibooker.cd();
0086     ibooker.setCurrentFolder(dirName);
0087 
0088     doAbsEta_ = hDim.doAbsEta;
0089     usePFMuon_ = hDim.usePFMuon;
0090 
0091     //histograms for efficiency plots
0092     hP_ = ibooker.book1D("P", "p of recoTracks", hDim.nBinP, hDim.minP, hDim.maxP);
0093     hPt_ = ibooker.book1D("Pt", "p_{T} of recoTracks", hDim.nBinPt, hDim.minPt, hDim.maxPt);
0094     hEta_ = ibooker.book1D("Eta", "#eta of recoTracks", hDim.nBinEta, hDim.minEta, hDim.maxEta);
0095     hPhi_ = ibooker.book1D("Phi", "#phi of recoTracks", hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
0096 
0097     hSimP_ = ibooker.book1D("SimP", "p of simTracks", hDim.nBinP, hDim.minP, hDim.maxP);
0098     hSimPt_ = ibooker.book1D("SimPt", "p_{T} of simTracks", hDim.nBinPt, hDim.minPt, hDim.maxPt);
0099     hSimEta_ = ibooker.book1D("SimEta", "#eta of simTracks", hDim.nBinEta, hDim.minEta, hDim.maxEta);
0100     hSimPhi_ = ibooker.book1D("SimPhi", "#phi of simTracks", hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
0101     hSimDxy_ = ibooker.book1D("SimDxy", "Dxy of simTracks", hDim.nBinDxy, hDim.minDxy, hDim.maxDxy);
0102     hSimDz_ = ibooker.book1D("Dz", "Dz of simTracks", hDim.nBinDz, hDim.minDz, hDim.maxDz);
0103 
0104     //track multiplicities
0105     hNSim_ = ibooker.book1D("NSim", "Number of particles per event", hDim.nTrks, -0.5, hDim.nTrks + 0.5);
0106     hNMuon_ = ibooker.book1D("NMuon", "Number of muons per event", hDim.nTrks, -0.5, hDim.nTrks + 0.5);
0107 
0108     // - Misc. variables
0109     hNTrks_ = ibooker.book1D("NTrks", "Number of reco tracks per event", hDim.nTrks, -0.5, hDim.nTrks + 0.5);
0110     hNTrksEta_ = ibooker.book1D("NTrksEta", "Number of reco tracks vs #eta", hDim.nBinEta, hDim.minEta, hDim.maxEta);
0111     hNTrksPt_ = ibooker.book1D("NTrksPt", "Number of reco tracks vs p_{T}", hDim.nBinPt, hDim.minPt, hDim.maxPt);
0112 
0113     hMisQPt_ = ibooker.book1D("MisQPt", "Charge mis-id vs Pt", hDim.nBinPt, hDim.minPt, hDim.maxPt);
0114     hMisQEta_ = ibooker.book1D("MisQEta", "Charge mis-id vs Eta", hDim.nBinEta, hDim.minEta, hDim.maxEta);
0115 
0116     // - Resolutions
0117     hErrP_ = ibooker.book1D("ErrP", "#Delta(p)/p", hDim.nBinErr, hDim.minErrP, hDim.maxErrP);
0118     hErrPBarrel_ = ibooker.book1D("ErrP_barrel", "#Delta(p)/p", hDim.nBinErr, hDim.minErrP, hDim.maxErrP);
0119     hErrPOverlap_ = ibooker.book1D("ErrP_overlap", "#Delta(p)/p", hDim.nBinErr, hDim.minErrP, hDim.maxErrP);
0120     hErrPEndcap_ = ibooker.book1D("ErrP_endcap", "#Delta(p)/p", hDim.nBinErr, hDim.minErrP, hDim.maxErrP);
0121     hErrPt_ = ibooker.book1D("ErrPt", "#Delta(p_{T})/p_{T}", hDim.nBinErr, hDim.minErrPt, hDim.maxErrPt);
0122     hErrPtBarrel_ = ibooker.book1D("ErrPt_barrel", "#Delta(p_{T})/p_{T}", hDim.nBinErr, hDim.minErrPt, hDim.maxErrPt);
0123     hErrPtOverlap_ = ibooker.book1D("ErrPt_overlap", "#Delta(p_{T})/p_{T}", hDim.nBinErr, hDim.minErrPt, hDim.maxErrPt);
0124     hErrPtEndcap_ = ibooker.book1D("ErrPt_endcap", "#Delta(p_{T})/p_{T}", hDim.nBinErr, hDim.minErrPt, hDim.maxErrPt);
0125     hErrEta_ = ibooker.book1D("ErrEta", "#sigma(#eta))", hDim.nBinErr, hDim.minErrEta, hDim.maxErrEta);
0126     hErrPhi_ = ibooker.book1D("ErrPhi", "#sigma(#phi)", hDim.nBinErr, hDim.minErrPhi, hDim.maxErrPhi);
0127     hErrDxy_ = ibooker.book1D("ErrDxy", "#sigma(d_{xy})", hDim.nBinErr, hDim.minErrDxy, hDim.maxErrDxy);
0128     hErrDz_ = ibooker.book1D("ErrDz", "#sigma(d_{z})", hDim.nBinErr, hDim.minErrDz, hDim.maxErrDz);
0129 
0130     //PF-RECO comparisons
0131     if (usePFMuon_) {
0132       hErrPt_PF_ = ibooker.book1D("ErrPt_PF", "#Delta(p_{T})|_{PF}/p_{T}", hDim.nBinErr, hDim.minErrPt, hDim.maxErrPt);
0133       hErrQPt_PF_ =
0134           ibooker.book1D("ErrQPt_PF", "#Delta(q/p_{T})|_{PF}/(q/p_{T})", hDim.nBinErr, hDim.minErrQPt, hDim.maxErrQPt);
0135 
0136       hPFMomAssCorrectness =
0137           ibooker.book1D("hPFMomAssCorrectness", "Corrected momentum assignement PF/RECO", 2, 0.5, 2.5);
0138       hPt_vs_PFMomAssCorrectness = ibooker.book2D("hPt_vs_PFMomAssCorrectness",
0139                                                   "Corrected momentum assignement PF/RECO",
0140                                                   hDim.nBinPt,
0141                                                   hDim.minPt,
0142                                                   hDim.maxP,
0143                                                   2,
0144                                                   0.5,
0145                                                   2.5);
0146 
0147       hdPt_vs_Pt_ = ibooker.book2D("dPt_vs_Pt",
0148                                    "#Delta(p_{T}) vs p_{T}",
0149                                    hDim.nBinPt,
0150                                    hDim.minPt,
0151                                    hDim.maxPt,
0152                                    hDim.nBinErr,
0153                                    hDim.minErrPt,
0154                                    hDim.maxErrPt);
0155       hdPt_vs_Eta_ = ibooker.book2D("dPt_vs_Eta",
0156                                     "#Delta(p_{T}) vs #eta",
0157                                     hDim.nBinEta,
0158                                     hDim.minEta,
0159                                     hDim.maxEta,
0160                                     hDim.nBinErr,
0161                                     hDim.minErrPt,
0162                                     hDim.maxErrPt);
0163     }
0164 
0165     // -- Resolutions vs Eta
0166     hErrP_vs_Eta_ = ibooker.book2D("ErrP_vs_Eta",
0167                                    "#Delta(p)/p vs #eta",
0168                                    hDim.nBinEta,
0169                                    hDim.minEta,
0170                                    hDim.maxEta,
0171                                    hDim.nBinErr,
0172                                    hDim.minErrP,
0173                                    hDim.maxErrP);
0174     hErrPt_vs_Eta_ = ibooker.book2D("ErrPt_vs_Eta",
0175                                     "#Delta(p_{T})/p_{T} vs #eta",
0176                                     hDim.nBinEta,
0177                                     hDim.minEta,
0178                                     hDim.maxEta,
0179                                     hDim.nBinErr,
0180                                     hDim.minErrPt,
0181                                     hDim.maxErrPt);
0182     hErrQPt_vs_Eta_ = ibooker.book2D("ErrQPt_vs_Eta",
0183                                      "#Delta(q/p_{T})/(q/p_{T}) vs #eta",
0184                                      hDim.nBinEta,
0185                                      hDim.minEta,
0186                                      hDim.maxEta,
0187                                      hDim.nBinErr,
0188                                      hDim.minErrQPt,
0189                                      hDim.maxErrQPt);
0190     hErrEta_vs_Eta_ = ibooker.book2D("ErrEta_vs_Eta",
0191                                      "#sigma(#eta) vs #eta",
0192                                      hDim.nBinEta,
0193                                      hDim.minEta,
0194                                      hDim.maxEta,
0195                                      hDim.nBinErr,
0196                                      hDim.minErrEta,
0197                                      hDim.maxErrEta);
0198 
0199     // -- Resolutions vs momentum
0200     hErrP_vs_P_ = ibooker.book2D(
0201         "ErrP_vs_P", "#Delta(p)/p vs p", hDim.nBinP, hDim.minP, hDim.maxP, hDim.nBinErr, hDim.minErrP, hDim.maxErrP);
0202     hErrPt_vs_Pt_ = ibooker.book2D("ErrPt_vs_Pt",
0203                                    "#Delta(p_{T})/p_{T} vs p_{T}",
0204                                    hDim.nBinPt,
0205                                    hDim.minPt,
0206                                    hDim.maxPt,
0207                                    hDim.nBinErr,
0208                                    hDim.minErrPt,
0209                                    hDim.maxErrPt);
0210     hErrQPt_vs_Pt_ = ibooker.book2D("ErrQPt_vs_Pt",
0211                                     "#Delta(q/p_{T})/(q/p_{T}) vs p_{T}",
0212                                     hDim.nBinPt,
0213                                     hDim.minPt,
0214                                     hDim.maxPt,
0215                                     hDim.nBinErr,
0216                                     hDim.minErrQPt,
0217                                     hDim.maxErrQPt);
0218 
0219     // - Pulls
0220     hPullPt_ = ibooker.book1D("PullPt", "Pull(#p_{T})", hDim.nBinPull, -hDim.wPull, hDim.wPull);
0221     hPullEta_ = ibooker.book1D("PullEta", "Pull(#eta)", hDim.nBinPull, -hDim.wPull, hDim.wPull);
0222     hPullPhi_ = ibooker.book1D("PullPhi", "Pull(#phi)", hDim.nBinPull, -hDim.wPull, hDim.wPull);
0223     hPullQPt_ = ibooker.book1D("PullQPt", "Pull(q/p_{T})", hDim.nBinPull, -hDim.wPull, hDim.wPull);
0224     hPullDxy_ = ibooker.book1D("PullDxy", "Pull(D_{xy})", hDim.nBinPull, -hDim.wPull, hDim.wPull);
0225     hPullDz_ = ibooker.book1D("PullDz", "Pull(D_{z})", hDim.nBinPull, -hDim.wPull, hDim.wPull);
0226 
0227     // -- Pulls vs Eta
0228     hPullPt_vs_Eta_ = ibooker.book2D("PullPt_vs_Eta",
0229                                      "Pull(p_{T}) vs #eta",
0230                                      hDim.nBinEta,
0231                                      hDim.minEta,
0232                                      hDim.maxEta,
0233                                      hDim.nBinPull,
0234                                      -hDim.wPull,
0235                                      hDim.wPull);
0236     hPullEta_vs_Eta_ = ibooker.book2D("PullEta_vs_Eta",
0237                                       "Pull(#eta) vs #eta",
0238                                       hDim.nBinEta,
0239                                       hDim.minEta,
0240                                       hDim.maxEta,
0241                                       hDim.nBinPull,
0242                                       -hDim.wPull,
0243                                       hDim.wPull);
0244     hPullPhi_vs_Eta_ = ibooker.book2D("PullPhi_vs_Eta",
0245                                       "Pull(#phi) vs #eta",
0246                                       hDim.nBinEta,
0247                                       hDim.minEta,
0248                                       hDim.maxEta,
0249                                       hDim.nBinPull,
0250                                       -hDim.wPull,
0251                                       hDim.wPull);
0252 
0253     // -- Pulls vs Pt
0254     hPullPt_vs_Pt_ = ibooker.book2D("PullPt_vs_Pt",
0255                                     "Pull(p_{T}) vs p_{T}",
0256                                     hDim.nBinPt,
0257                                     hDim.minPt,
0258                                     hDim.maxPt,
0259                                     hDim.nBinPull,
0260                                     -hDim.wPull,
0261                                     hDim.wPull);
0262     hPullEta_vs_Pt_ = ibooker.book2D("PullEta_vs_Pt",
0263                                      "Pull(#eta) vs p_{T}",
0264                                      hDim.nBinPt,
0265                                      hDim.minPt,
0266                                      hDim.maxPt,
0267                                      hDim.nBinPull,
0268                                      -hDim.wPull,
0269                                      hDim.wPull);
0270 
0271     // -- Number of Hits
0272     const int nHits = 100;
0273     hNHits_ = ibooker.book1D("NHits", "Number of hits", nHits + 1, -0.5, nHits + 0.5);
0274     hNHits_vs_Pt_ = ibooker.book2D("NHits_vs_Pt",
0275                                    "Number of hits vs p_{T}",
0276                                    hDim.nBinPt,
0277                                    hDim.minPt,
0278                                    hDim.maxPt,
0279                                    nHits / 4 + 1,
0280                                    -0.25,
0281                                    nHits + 0.25);
0282     hNHits_vs_Eta_ = ibooker.book2D("NHits_vs_Eta",
0283                                     "Number of hits vs #eta",
0284                                     hDim.nBinEta,
0285                                     hDim.minEta,
0286                                     hDim.maxEta,
0287                                     nHits / 4 + 1,
0288                                     -0.25,
0289                                     nHits + 0.25);
0290     hNSimHits_ = ibooker.book1D("NSimHits", "Number of simHits", nHits + 1, -0.5, nHits + 0.5);
0291 
0292     const int nLostHits = 5;
0293     hNLostHits_ = ibooker.book1D("NLostHits", "Number of Lost hits", nLostHits + 1, -0.5, nLostHits + 0.5);
0294     hNLostHits_vs_Pt_ = ibooker.book2D("NLostHits_vs_Pt",
0295                                        "Number of lost Hits vs p_{T}",
0296                                        hDim.nBinPt,
0297                                        hDim.minPt,
0298                                        hDim.maxPt,
0299                                        nLostHits + 1,
0300                                        -0.5,
0301                                        nLostHits + 0.5);
0302     hNLostHits_vs_Eta_ = ibooker.book2D("NLostHits_vs_Eta",
0303                                         "Number of lost Hits vs #eta",
0304                                         hDim.nBinEta,
0305                                         hDim.minEta,
0306                                         hDim.maxEta,
0307                                         nLostHits + 1,
0308                                         -0.5,
0309                                         nLostHits + 0.5);
0310 
0311     const int nTrackerHits = 40;
0312     hNTrackerHits_ =
0313         ibooker.book1D("NTrackerHits", "Number of valid tracker hits", nTrackerHits + 1, -0.5, nTrackerHits + 0.5);
0314     hNTrackerHits_vs_Pt_ = ibooker.book2D("NTrackerHits_vs_Pt",
0315                                           "Number of valid traker hits vs p_{T}",
0316                                           hDim.nBinPt,
0317                                           hDim.minPt,
0318                                           hDim.maxPt,
0319                                           nTrackerHits / 4 + 1,
0320                                           -0.25,
0321                                           nTrackerHits + 0.25);
0322     hNTrackerHits_vs_Eta_ = ibooker.book2D("NTrackerHits_vs_Eta",
0323                                            "Number of valid tracker hits vs #eta",
0324                                            hDim.nBinEta,
0325                                            hDim.minEta,
0326                                            hDim.maxEta,
0327                                            nTrackerHits / 4 + 1,
0328                                            -0.25,
0329                                            nTrackerHits + 0.25);
0330 
0331     const int nMuonHits = 60;
0332     hNMuonHits_ = ibooker.book1D("NMuonHits", "Number of valid muon hits", nMuonHits + 1, -0.5, nMuonHits + 0.5);
0333     hNMuonHits_vs_Pt_ = ibooker.book2D("NMuonHits_vs_Pt",
0334                                        "Number of valid muon hits vs p_{T}",
0335                                        hDim.nBinPt,
0336                                        hDim.minPt,
0337                                        hDim.maxPt,
0338                                        nMuonHits / 4 + 1,
0339                                        -0.25,
0340                                        nMuonHits + 0.25);
0341     hNMuonHits_vs_Eta_ = ibooker.book2D("NMuonHits_vs_Eta",
0342                                         "Number of valid muon hits vs #eta",
0343                                         hDim.nBinEta,
0344                                         hDim.minEta,
0345                                         hDim.maxEta,
0346                                         nMuonHits / 4 + 1,
0347                                         -0.25,
0348                                         nMuonHits + 0.25);
0349 
0350     hNDof_ = ibooker.book1D("NDof", "Number of DoF", hDim.nDof + 1, -0.5, hDim.nDof + 0.5);
0351     hChi2_ = ibooker.book1D("Chi2", "#Chi^{2}", hDim.nBinErr, 0, 200);
0352     hChi2Norm_ = ibooker.book1D("Chi2Norm", "Normalized #Chi^{2}", hDim.nBinErr, 0, 50);
0353     hChi2Prob_ = ibooker.book1D("Chi2Prob", "Prob(#Chi^{2})", hDim.nBinErr, 0, 1);
0354 
0355     hNDof_vs_Eta_ = ibooker.book2D("NDof_vs_Eta",
0356                                    "Number of DoF vs #eta",
0357                                    hDim.nBinEta,
0358                                    hDim.minEta,
0359                                    hDim.maxEta,
0360                                    hDim.nDof + 1,
0361                                    -0.5,
0362                                    hDim.nDof + 0.5);
0363     hChi2_vs_Eta_ = ibooker.book2D(
0364         "Chi2_vs_Eta", "#Chi^{2} vs #eta", hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinErr, 0., 200.);
0365     hChi2Norm_vs_Eta_ = ibooker.book2D(
0366         "Chi2Norm_vs_Eta", "Normalized #Chi^{2} vs #eta", hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinErr, 0., 50.);
0367     hChi2Prob_vs_Eta_ = ibooker.book2D(
0368         "Chi2Prob_vs_Eta", "Prob(#Chi^{2}) vs #eta", hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinErr, 0., 1.);
0369 
0370     hNSimToReco_ =
0371         ibooker.book1D("NSimToReco", "Number of associated reco tracks", hDim.nAssoc + 1, -0.5, hDim.nAssoc + 0.5);
0372     hNRecoToSim_ =
0373         ibooker.book1D("NRecoToSim", "Number of associated sim TP's", hDim.nAssoc + 1, -0.5, hDim.nAssoc + 0.5);
0374   };
0375 
0376   //
0377   //Fill hists booked, simRef and muonRef are associated by hits
0378   //
0379   void fill(const TrackingParticle* simRef, const Muon* muonRef) {
0380     const double simP = simRef->p();
0381     const double simPt = simRef->pt();
0382     const double simEta = doAbsEta_ ? fabs(simRef->eta()) : simRef->eta();
0383     const double simPhi = simRef->phi();
0384     const double simQ = simRef->charge();
0385     const double simQPt = simQ / simPt;
0386 
0387     GlobalPoint simVtx(simRef->vertex().x(), simRef->vertex().y(), simRef->vertex().z());
0388     GlobalVector simMom(simRef->momentum().x(), simRef->momentum().y(), simRef->momentum().z());
0389     const double simDxy = -simVtx.x() * sin(simPhi) + simVtx.y() * cos(simPhi);
0390     const double simDz = simVtx.z() - (simVtx.x() * simMom.x() + simVtx.y() * simMom.y()) * simMom.z() / simMom.perp2();
0391 
0392     const double recoQ = muonRef->charge();
0393     if (simQ * recoQ < 0) {
0394       hMisQPt_->Fill(simPt);
0395       hMisQEta_->Fill(simEta);
0396     }
0397 
0398     double recoP, recoPt, recoEta, recoPhi, recoQPt;
0399     if (usePFMuon_) {
0400       //      const double origRecoP   = muonRef->p();
0401       const double origRecoPt = muonRef->pt();
0402       //      const double origRecoEta = muonRef->eta();
0403       //      const double origRecoPhi = muonRef->phi();
0404       const double origRecoQPt = recoQ / origRecoPt;
0405       recoP = muonRef->pfP4().P();
0406       recoPt = muonRef->pfP4().Pt();
0407       recoEta = muonRef->pfP4().Eta();
0408       recoPhi = muonRef->pfP4().Phi();
0409       recoQPt = recoQ / recoPt;
0410       hErrPt_PF_->Fill((recoPt - origRecoPt) / origRecoPt);
0411       hErrQPt_PF_->Fill((recoQPt - origRecoQPt) / origRecoQPt);
0412 
0413       hdPt_vs_Eta_->Fill(recoEta, recoPt - origRecoPt);
0414       hdPt_vs_Pt_->Fill(recoPt, recoPt - origRecoPt);
0415 
0416       int theCorrectPFAss = (fabs(recoPt - simPt) < fabs(origRecoPt - simPt)) ? 1 : 2;
0417       hPFMomAssCorrectness->Fill(theCorrectPFAss);
0418       hPt_vs_PFMomAssCorrectness->Fill(simPt, theCorrectPFAss);
0419     }
0420 
0421     else {
0422       recoP = muonRef->p();
0423       recoPt = muonRef->pt();
0424       recoEta = muonRef->eta();
0425       recoPhi = muonRef->phi();
0426       recoQPt = recoQ / recoPt;
0427     }
0428 
0429     const double errP = (recoP - simP) / simP;
0430     const double errPt = (recoPt - simPt) / simPt;
0431     const double errEta = (recoEta - simEta) / simEta;
0432     const double errPhi = (recoPhi - simPhi) / simPhi;
0433     const double errQPt = (recoQPt - simQPt) / simQPt;
0434 
0435     hP_->Fill(simP);
0436     hPt_->Fill(simPt);
0437     hEta_->Fill(simEta);
0438     hPhi_->Fill(simPhi);
0439 
0440     hErrP_->Fill(errP);
0441     hErrPt_->Fill(errPt);
0442     hErrEta_->Fill(errEta);
0443     hErrPhi_->Fill(errPhi);
0444 
0445     if (fabs(simEta) > 0. && fabs(simEta) < 0.8) {
0446       hErrPBarrel_->Fill(errP);
0447       hErrPtBarrel_->Fill(errPt);
0448     } else if (fabs(simEta) > 0.8 && fabs(simEta) < 1.2) {
0449       hErrPOverlap_->Fill(errP);
0450       hErrPtOverlap_->Fill(errPt);
0451     } else if (fabs(simEta) > 1.2) {
0452       hErrPEndcap_->Fill(errP);
0453       hErrPtEndcap_->Fill(errPt);
0454     }
0455 
0456     hErrP_vs_Eta_->Fill(simEta, errP);
0457     hErrPt_vs_Eta_->Fill(simEta, errPt);
0458     hErrQPt_vs_Eta_->Fill(simEta, errQPt);
0459 
0460     hErrP_vs_P_->Fill(simP, errP);
0461     hErrPt_vs_Pt_->Fill(simPt, errPt);
0462     hErrQPt_vs_Pt_->Fill(simQPt, errQPt);
0463 
0464     hErrEta_vs_Eta_->Fill(simEta, errEta);
0465 
0466     //access from track
0467     reco::TrackRef recoRef = muonRef->track();
0468     if (recoRef.isNonnull()) {
0469       // Number of reco-hits
0470       const int nRecoHits = recoRef->numberOfValidHits();
0471       const int nLostHits = recoRef->numberOfLostHits();
0472 
0473       hNHits_->Fill(nRecoHits);
0474       hNHits_vs_Pt_->Fill(simPt, nRecoHits);
0475       hNHits_vs_Eta_->Fill(simEta, nRecoHits);
0476 
0477       hNLostHits_->Fill(nLostHits);
0478       hNLostHits_vs_Pt_->Fill(simPt, nLostHits);
0479       hNLostHits_vs_Eta_->Fill(simEta, nLostHits);
0480 
0481       const double recoNDof = recoRef->ndof();
0482       const double recoChi2 = recoRef->chi2();
0483       const double recoChi2Norm = recoRef->normalizedChi2();
0484       const double recoChi2Prob = TMath::Prob(recoRef->chi2(), static_cast<int>(recoRef->ndof()));
0485 
0486       hNDof_->Fill(recoNDof);
0487       hChi2_->Fill(recoChi2);
0488       hChi2Norm_->Fill(recoChi2Norm);
0489       hChi2Prob_->Fill(recoChi2Prob);
0490 
0491       hNDof_vs_Eta_->Fill(simEta, recoNDof);
0492       hChi2_vs_Eta_->Fill(simEta, recoChi2);
0493       hChi2Norm_vs_Eta_->Fill(simEta, recoChi2Norm);
0494       hChi2Prob_vs_Eta_->Fill(simEta, recoChi2Prob);
0495 
0496       const double recoDxy = recoRef->dxy();
0497       const double recoDz = recoRef->dz();
0498 
0499       const double errDxy = (recoDxy - simDxy) / simDxy;
0500       const double errDz = (recoDz - simDz) / simDz;
0501       hErrDxy_->Fill(errDxy);
0502       hErrDz_->Fill(errDz);
0503 
0504       const double pullPt = (recoPt - simPt) / recoRef->ptError();
0505       const double pullQPt = (recoQPt - simQPt) / recoRef->qoverpError();
0506       const double pullEta = (recoEta - simEta) / recoRef->etaError();
0507       const double pullPhi = (recoPhi - simPhi) / recoRef->phiError();
0508       const double pullDxy = (recoDxy - simDxy) / recoRef->dxyError();
0509       const double pullDz = (recoDz - simDz) / recoRef->dzError();
0510 
0511       hPullPt_->Fill(pullPt);
0512       hPullEta_->Fill(pullEta);
0513       hPullPhi_->Fill(pullPhi);
0514       hPullQPt_->Fill(pullQPt);
0515       hPullDxy_->Fill(pullDxy);
0516       hPullDz_->Fill(pullDz);
0517 
0518       hPullPt_vs_Eta_->Fill(simEta, pullPt);
0519       hPullPt_vs_Pt_->Fill(simPt, pullPt);
0520 
0521       hPullEta_vs_Eta_->Fill(simEta, pullEta);
0522       hPullPhi_vs_Eta_->Fill(simEta, pullPhi);
0523 
0524       hPullEta_vs_Pt_->Fill(simPt, pullEta);
0525     }
0526   };
0527 };
0528 
0529 //
0530 //struct defininiong histograms
0531 //
0532 struct RecoMuonValidator::CommonME {
0533   typedef MonitorElement* MEP;
0534 
0535   //diffs
0536   MEP hTrkToGlbDiffNTrackerHits_, hStaToGlbDiffNMuonHits_;
0537   MEP hTrkToGlbDiffNTrackerHitsEta_, hStaToGlbDiffNMuonHitsEta_;
0538   MEP hTrkToGlbDiffNTrackerHitsPt_, hStaToGlbDiffNMuonHitsPt_;
0539 
0540   //global muon hit pattern
0541   MEP hNInvalidHitsGTHitPattern_, hNInvalidHitsITHitPattern_, hNInvalidHitsOTHitPattern_;
0542   MEP hNDeltaInvalidHitsHitPattern_;
0543 
0544   //muon based momentum assignment
0545   MEP hMuonP_, hMuonPt_, hMuonEta_, hMuonPhi_;
0546   //track based kinematics
0547   MEP hMuonTrackP_, hMuonTrackPt_, hMuonTrackEta_, hMuonTrackPhi_, hMuonTrackDxy_, hMuonTrackDz_;
0548   //histograms for fractions
0549   MEP hMuonAllP_, hMuonAllPt_, hMuonAllEta_, hMuonAllPhi_;
0550 };
0551 
0552 //
0553 //Constructor
0554 //
0555 RecoMuonValidator::RecoMuonValidator(const edm::ParameterSet& pset)
0556     : selector_(pset.getParameter<std::string>("selection")) {
0557   // dump cfg parameters
0558   edm::LogVerbatim("RecoMuonValidator") << "constructing RecoMuonValidator: " << pset.dump();
0559   verbose_ = pset.getUntrackedParameter<unsigned int>("verbose", 0);
0560 
0561   outputFileName_ = pset.getUntrackedParameter<string>("outputFileName", "");
0562 
0563   wantTightMuon_ = pset.getParameter<bool>("wantTightMuon");
0564   beamspotLabel_ = pset.getParameter<edm::InputTag>("beamSpot");
0565   primvertexLabel_ = pset.getParameter<edm::InputTag>("primaryVertex");
0566   beamspotToken_ = consumes<reco::BeamSpot>(beamspotLabel_);
0567   primvertexToken_ = consumes<reco::VertexCollection>(primvertexLabel_);
0568 
0569   // Set histogram dimensions from config
0570 
0571   hDim.nBinP = pset.getUntrackedParameter<unsigned int>("nBinP");
0572   hDim.minP = pset.getUntrackedParameter<double>("minP");
0573   hDim.maxP = pset.getUntrackedParameter<double>("maxP");
0574 
0575   hDim.nBinPt = pset.getUntrackedParameter<unsigned int>("nBinPt");
0576   hDim.minPt = pset.getUntrackedParameter<double>("minPt");
0577   hDim.maxPt = pset.getUntrackedParameter<double>("maxPt");
0578 
0579   doAbsEta_ = pset.getUntrackedParameter<bool>("doAbsEta");
0580   hDim.doAbsEta = doAbsEta_;
0581   hDim.nBinEta = pset.getUntrackedParameter<unsigned int>("nBinEta");
0582   hDim.minEta = pset.getUntrackedParameter<double>("minEta");
0583   hDim.maxEta = pset.getUntrackedParameter<double>("maxEta");
0584 
0585   hDim.nBinDxy = pset.getUntrackedParameter<unsigned int>("nBinDxy");
0586   hDim.minDxy = pset.getUntrackedParameter<double>("minDxy");
0587   hDim.maxDxy = pset.getUntrackedParameter<double>("maxDxy");
0588 
0589   hDim.nBinDz = pset.getUntrackedParameter<unsigned int>("nBinDz");
0590   hDim.minDz = pset.getUntrackedParameter<double>("minDz");
0591   hDim.maxDz = pset.getUntrackedParameter<double>("maxDz");
0592 
0593   hDim.nBinPhi = pset.getUntrackedParameter<unsigned int>("nBinPhi");
0594   hDim.minPhi = pset.getUntrackedParameter<double>("minPhi", -TMath::Pi());
0595   hDim.maxPhi = pset.getUntrackedParameter<double>("maxPhi", TMath::Pi());
0596 
0597   hDim.nBinErr = pset.getUntrackedParameter<unsigned int>("nBinErr");
0598   hDim.nBinPull = pset.getUntrackedParameter<unsigned int>("nBinPull");
0599 
0600   hDim.wPull = pset.getUntrackedParameter<double>("wPull");
0601 
0602   hDim.minErrP = pset.getUntrackedParameter<double>("minErrP");
0603   hDim.maxErrP = pset.getUntrackedParameter<double>("maxErrP");
0604 
0605   hDim.minErrPt = pset.getUntrackedParameter<double>("minErrPt");
0606   hDim.maxErrPt = pset.getUntrackedParameter<double>("maxErrPt");
0607 
0608   hDim.minErrQPt = pset.getUntrackedParameter<double>("minErrQPt");
0609   hDim.maxErrQPt = pset.getUntrackedParameter<double>("maxErrQPt");
0610 
0611   hDim.minErrEta = pset.getUntrackedParameter<double>("minErrEta");
0612   hDim.maxErrEta = pset.getUntrackedParameter<double>("maxErrEta");
0613 
0614   hDim.minErrPhi = pset.getUntrackedParameter<double>("minErrPhi");
0615   hDim.maxErrPhi = pset.getUntrackedParameter<double>("maxErrPhi");
0616 
0617   hDim.minErrDxy = pset.getUntrackedParameter<double>("minErrDxy");
0618   hDim.maxErrDxy = pset.getUntrackedParameter<double>("maxErrDxy");
0619 
0620   hDim.minErrDz = pset.getUntrackedParameter<double>("minErrDz");
0621   hDim.maxErrDz = pset.getUntrackedParameter<double>("maxErrDz");
0622 
0623   hDim.nTrks = pset.getUntrackedParameter<unsigned int>("nTrks");
0624   hDim.nAssoc = pset.getUntrackedParameter<unsigned int>("nAssoc");
0625   hDim.nDof = pset.getUntrackedParameter<unsigned int>("nDof", 55);
0626 
0627   // Labels for simulation and reconstruction tracks
0628   simLabel_ = pset.getParameter<InputTag>("simLabel");
0629   tpRefVector = pset.getParameter<bool>("tpRefVector");
0630   if (tpRefVector)
0631     tpRefVectorToken_ = consumes<TrackingParticleRefVector>(simLabel_);
0632   else
0633     simToken_ = consumes<TrackingParticleCollection>(simLabel_);
0634 
0635   muonLabel_ = pset.getParameter<InputTag>("muonLabel");
0636   muonToken_ = consumes<edm::View<reco::Muon> >(muonLabel_);
0637 
0638   // Labels for sim-reco association
0639   doAssoc_ = pset.getUntrackedParameter<bool>("doAssoc", true);
0640   muAssocLabel_ = pset.getParameter<InputTag>("muAssocLabel");
0641   if (doAssoc_) {
0642     muAssocToken_ = consumes<reco::MuonToTrackingParticleAssociator>(muAssocLabel_);
0643   }
0644 
0645   // Different momentum assignment and additional histos in case of PF muons
0646   usePFMuon_ = pset.getUntrackedParameter<bool>("usePFMuon");
0647   hDim.usePFMuon = usePFMuon_;
0648 
0649   //type of track
0650   std::string trackType = pset.getParameter<std::string>("trackType");
0651   if (trackType == "inner")
0652     trackType_ = reco::InnerTk;
0653   else if (trackType == "outer")
0654     trackType_ = reco::OuterTk;
0655   else if (trackType == "global")
0656     trackType_ = reco::GlobalTk;
0657   else if (trackType == "segments")
0658     trackType_ = reco::Segments;
0659   else
0660     throw cms::Exception("Configuration") << "Track type '" << trackType << "' not supported.\n";
0661 
0662   //  seedPropagatorName_ = pset.getParameter<string>("SeedPropagator");
0663 
0664   ParameterSet tpset = pset.getParameter<ParameterSet>("tpSelector");
0665   tpSelector_ = TrackingParticleSelector(tpset.getParameter<double>("ptMin"),
0666                                          tpset.getParameter<double>("ptMax"),
0667                                          tpset.getParameter<double>("minRapidity"),
0668                                          tpset.getParameter<double>("maxRapidity"),
0669                                          tpset.getParameter<double>("tip"),
0670                                          tpset.getParameter<double>("lip"),
0671                                          tpset.getParameter<int>("minHit"),
0672                                          tpset.getParameter<bool>("signalOnly"),
0673                                          tpset.getParameter<bool>("intimeOnly"),
0674                                          tpset.getParameter<bool>("chargedOnly"),
0675                                          tpset.getParameter<bool>("stableOnly"),
0676                                          tpset.getParameter<std::vector<int> >("pdgId"));
0677 
0678   // the service parameters
0679   ParameterSet serviceParameters = pset.getParameter<ParameterSet>("ServiceParameters");
0680 
0681   // retrieve the instance of DQMService
0682   dbe_ = Service<DQMStore>().operator->();
0683   subsystemname_ = pset.getUntrackedParameter<std::string>("subSystemFolder", "YourSubsystem");
0684 
0685   subDir_ = pset.getUntrackedParameter<string>("subDir");
0686   if (subDir_.empty())
0687     subDir_ = "RecoMuonV";
0688   if (subDir_[subDir_.size() - 1] == '/')
0689     subDir_.erase(subDir_.size() - 1);
0690 }
0691 
0692 void RecoMuonValidator::bookHistograms(DQMStore::IBooker& ibooker,
0693                                        edm::Run const& iRun,
0694                                        edm::EventSetup const& /* iSetup */) {
0695   // book histograms
0696   ibooker.cd();
0697 
0698   ibooker.setCurrentFolder(subDir_);
0699 
0700   commonME_ = new CommonME;
0701   muonME_ = new MuonME;
0702 
0703   //commonME
0704   const int nHits = 100;
0705 
0706   // - diffs
0707   commonME_->hTrkToGlbDiffNTrackerHits_ = ibooker.book1D("TrkGlbDiffNTrackerHits",
0708                                                          "Difference of number of tracker hits (tkMuon - globalMuon)",
0709                                                          2 * nHits + 1,
0710                                                          -nHits - 0.5,
0711                                                          nHits + 0.5);
0712   commonME_->hStaToGlbDiffNMuonHits_ = ibooker.book1D("StaGlbDiffNMuonHits",
0713                                                       "Difference of number of muon hits (staMuon - globalMuon)",
0714                                                       2 * nHits + 1,
0715                                                       -nHits - 0.5,
0716                                                       nHits + 0.5);
0717 
0718   commonME_->hTrkToGlbDiffNTrackerHitsEta_ =
0719       ibooker.book2D("TrkGlbDiffNTrackerHitsEta",
0720                      "Difference of number of tracker hits (tkMuon - globalMuon)",
0721                      hDim.nBinEta,
0722                      hDim.minEta,
0723                      hDim.maxEta,
0724                      2 * nHits + 1,
0725                      -nHits - 0.5,
0726                      nHits + 0.5);
0727   commonME_->hStaToGlbDiffNMuonHitsEta_ = ibooker.book2D("StaGlbDiffNMuonHitsEta",
0728                                                          "Difference of number of muon hits (staMuon - globalMuon)",
0729                                                          hDim.nBinEta,
0730                                                          hDim.minEta,
0731                                                          hDim.maxEta,
0732                                                          2 * nHits + 1,
0733                                                          -nHits - 0.5,
0734                                                          nHits + 0.5);
0735 
0736   commonME_->hTrkToGlbDiffNTrackerHitsPt_ = ibooker.book2D("TrkGlbDiffNTrackerHitsPt",
0737                                                            "Difference of number of tracker hits (tkMuon - globalMuon)",
0738                                                            hDim.nBinPt,
0739                                                            hDim.minPt,
0740                                                            hDim.maxPt,
0741                                                            2 * nHits + 1,
0742                                                            -nHits - 0.5,
0743                                                            nHits + 0.5);
0744   commonME_->hStaToGlbDiffNMuonHitsPt_ = ibooker.book2D("StaGlbDiffNMuonHitsPt",
0745                                                         "Difference of number of muon hits (staMuon - globalMuon)",
0746                                                         hDim.nBinPt,
0747                                                         hDim.minPt,
0748                                                         hDim.maxPt,
0749                                                         2 * nHits + 1,
0750                                                         -nHits - 0.5,
0751                                                         nHits + 0.5);
0752 
0753   // -global muon hit pattern
0754   commonME_->hNInvalidHitsGTHitPattern_ = ibooker.book1D(
0755       "NInvalidHitsGTHitPattern", "Number of invalid hits on a global track", nHits + 1, -0.5, nHits + 0.5);
0756   commonME_->hNInvalidHitsITHitPattern_ = ibooker.book1D(
0757       "NInvalidHitsITHitPattern", "Number of invalid hits on an inner track", nHits + 1, -0.5, nHits + 0.5);
0758   commonME_->hNInvalidHitsOTHitPattern_ = ibooker.book1D(
0759       "NInvalidHitsOTHitPattern", "Number of invalid hits on an outer track", nHits + 1, -0.5, nHits + 0.5);
0760   commonME_->hNDeltaInvalidHitsHitPattern_ =
0761       ibooker.book1D("hNDeltaInvalidHitsHitPattern",
0762                      "The discrepancy for Number of invalid hits on an global track and inner and outer tracks",
0763                      2 * nHits + 1,
0764                      -nHits - 0.5,
0765                      nHits + 0.5);
0766 
0767   //muon based kinematics
0768   commonME_->hMuonP_ = ibooker.book1D("PMuon", "p of muon", hDim.nBinP, hDim.minP, hDim.maxP);
0769   commonME_->hMuonPt_ = ibooker.book1D("PtMuon", "p_{T} of muon", hDim.nBinPt, hDim.minPt, hDim.maxPt);
0770   commonME_->hMuonEta_ = ibooker.book1D("EtaMuon", "#eta of muon", hDim.nBinEta, hDim.minEta, hDim.maxEta);
0771   commonME_->hMuonPhi_ = ibooker.book1D("PhiMuon", "#phi of muon", hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
0772   //track based kinematics
0773   commonME_->hMuonTrackP_ = ibooker.book1D("PMuonTrack", "p of reco muon track", hDim.nBinP, hDim.minP, hDim.maxP);
0774   commonME_->hMuonTrackPt_ =
0775       ibooker.book1D("PtMuonTrack", "p_{T} of reco muon track", hDim.nBinPt, hDim.minPt, hDim.maxPt);
0776   commonME_->hMuonTrackEta_ =
0777       ibooker.book1D("EtaMuonTrack", "#eta of reco muon track", hDim.nBinEta, hDim.minEta, hDim.maxEta);
0778   commonME_->hMuonTrackPhi_ =
0779       ibooker.book1D("PhiMuonTrack", "#phi of reco muon track", hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
0780   commonME_->hMuonTrackDxy_ =
0781       ibooker.book1D("DxyMuonTrack", "Dxy of reco muon track", hDim.nBinDxy, hDim.minDxy, hDim.maxDxy);
0782   commonME_->hMuonTrackDz_ =
0783       ibooker.book1D("DzMuonTrack", "Dz of reco muon track", hDim.nBinDz, hDim.minDz, hDim.maxDz);
0784 
0785   //histograms for fractions
0786   commonME_->hMuonAllP_ = ibooker.book1D("PMuonAll", "p of muons of all types", hDim.nBinP, hDim.minP, hDim.maxP);
0787   commonME_->hMuonAllPt_ =
0788       ibooker.book1D("PtMuonAll", "p_{T} of muon of all types", hDim.nBinPt, hDim.minPt, hDim.maxPt);
0789   commonME_->hMuonAllEta_ =
0790       ibooker.book1D("EtaMuonAll", "#eta of muon of all types", hDim.nBinEta, hDim.minEta, hDim.maxEta);
0791   commonME_->hMuonAllPhi_ =
0792       ibooker.book1D("PhiMuonAll", "#phi of muon of all types", hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
0793 
0794   muonME_->bookHistos(ibooker, subDir_, hDim);
0795 }
0796 
0797 //
0798 //Destructor
0799 //
0800 RecoMuonValidator::~RecoMuonValidator() {}
0801 
0802 //
0803 //Begin run
0804 //
0805 
0806 void RecoMuonValidator::dqmBeginRun(const edm::Run&, const EventSetup& eventSetup) {}
0807 
0808 //
0809 //End run
0810 //
0811 void RecoMuonValidator::dqmEndRun(edm::Run const&, edm::EventSetup const&) {
0812   if (dbe_ && !outputFileName_.empty())
0813     dbe_->save(outputFileName_);
0814 }
0815 
0816 //
0817 //Analyze
0818 //
0819 void RecoMuonValidator::analyze(const Event& event, const EventSetup& eventSetup) {
0820   // Look for the Primary Vertex (and use the BeamSpot instead, if you can't find it):
0821   reco::Vertex::Point posVtx;
0822   reco::Vertex::Error errVtx;
0823   edm::Handle<reco::VertexCollection> recVtxs;
0824   event.getByToken(primvertexToken_, recVtxs);
0825   unsigned int theIndexOfThePrimaryVertex = 999.;
0826   for (unsigned int ind = 0; ind < recVtxs->size(); ++ind) {
0827     if ((*recVtxs)[ind].isValid() && !((*recVtxs)[ind].isFake())) {
0828       theIndexOfThePrimaryVertex = ind;
0829       break;
0830     }
0831   }
0832   if (theIndexOfThePrimaryVertex < 100) {
0833     posVtx = ((*recVtxs)[theIndexOfThePrimaryVertex]).position();
0834     errVtx = ((*recVtxs)[theIndexOfThePrimaryVertex]).error();
0835   } else {
0836     LogInfo("RecoMuonValidator") << "reco::PrimaryVertex not found, use BeamSpot position instead\n";
0837     edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0838     event.getByToken(beamspotToken_, recoBeamSpotHandle);
0839     reco::BeamSpot bs = *recoBeamSpotHandle;
0840     posVtx = bs.position();
0841     errVtx(0, 0) = bs.BeamWidthX();
0842     errVtx(1, 1) = bs.BeamWidthY();
0843     errVtx(2, 2) = bs.sigmaZ();
0844   }
0845   const reco::Vertex thePrimaryVertex(posVtx, errVtx);
0846 
0847   // Get TrackingParticles
0848   TrackingParticleRefVector TPrefV;
0849   const TrackingParticleRefVector* ptr_TPrefV = nullptr;
0850   Handle<TrackingParticleCollection> simHandle;
0851   Handle<TrackingParticleRefVector> TPCollectionRefVector_H;
0852 
0853   if (tpRefVector) {
0854     event.getByToken(tpRefVectorToken_, TPCollectionRefVector_H);
0855     ptr_TPrefV = TPCollectionRefVector_H.product();
0856     TPrefV = *ptr_TPrefV;
0857   } else {
0858     event.getByToken(simToken_, simHandle);
0859     size_t nTP = simHandle->size();
0860     for (size_t i = 0; i < nTP; ++i) {
0861       TPrefV.push_back(TrackingParticleRef(simHandle, i));
0862     }
0863     ptr_TPrefV = &TPrefV;
0864   }
0865 
0866   // Get Muons
0867   Handle<edm::View<Muon> > muonHandle;
0868   event.getByToken(muonToken_, muonHandle);
0869   View<Muon> muonColl = *(muonHandle.product());
0870 
0871   reco::MuonToTrackingParticleAssociator const* assoByHits = nullptr;
0872   if (doAssoc_) {
0873     edm::Handle<reco::MuonToTrackingParticleAssociator> associatorBase;
0874     event.getByToken(muAssocToken_, associatorBase);
0875     assoByHits = associatorBase.product();
0876   }
0877 
0878   const size_t nSim = ptr_TPrefV->size();
0879 
0880   edm::RefToBaseVector<reco::Muon> Muons;
0881   for (size_t i = 0; i < muonHandle->size(); ++i) {
0882     Muons.push_back(muonHandle->refAt(i));
0883   }
0884 
0885   muonME_->hNSim_->Fill(nSim);
0886   muonME_->hNMuon_->Fill(muonColl.size());
0887 
0888   reco::MuonToSimCollection muonToSimColl;
0889   reco::SimToMuonCollection simToMuonColl;
0890 
0891   if (doAssoc_) {
0892     edm::LogVerbatim("RecoMuonValidator")
0893         << "\n >>> MuonToSim association : " << muAssocLabel_ << " <<< \n"
0894         << "     muon collection : " << muonLabel_ << " (size = " << muonHandle->size() << ") \n"
0895         << "     TrackingParticle collection : " << simLabel_ << " (size = " << nSim << ")";
0896 
0897     assoByHits->associateMuons(muonToSimColl, simToMuonColl, Muons, trackType_, TPrefV);
0898   } else {
0899     /*
0900     // SimToMuon associations
0901     Handle<reco::RecoToSimCollection> simToTrkMuHandle;
0902     event.getByLabel(trkMuAssocLabel_, simToTrkMuHandle);
0903     trkSimRecColl = *(simToTrkMuHandle.product());
0904 
0905     Handle<reco::RecoToSimCollection> simToStaMuHandle;
0906     event.getByLabel(staMuAssocLabel_, simToStaMuHandle);
0907     staSimRecColl = *(simToStaMuHandle.product());
0908 
0909     Handle<reco::RecoToSimCollection> simToGlbMuHandle;
0910     event.getByLabel(glbMuAssocLabel_, simToGlbMuHandle);
0911     glbSimRecColl = *(simToGlbMuHandle.product());
0912 
0913     // MuonToSim associations
0914     Handle<reco::SimToRecoCollection> trkMuToSimHandle;
0915     event.getByLabel(trkMuAssocLabel_, trkMuToSimHandle);
0916     trkRecSimColl = *(trkMuToSimHandle.product());
0917 
0918     Handle<reco::SimToRecoCollection> staMuToSimHandle;
0919     event.getByLabel(staMuAssocLabel_, staMuToSimHandle);
0920     staRecSimColl = *(staMuToSimHandle.product());
0921 
0922     Handle<reco::SimToRecoCollection> glbMuToSimHandle;
0923     event.getByLabel(glbMuAssocLabel_, glbMuToSimHandle);
0924     glbRecSimColl = *(glbMuToSimHandle.product());
0925 */
0926   }
0927 
0928   int glbNTrackerHits = 0;
0929   int trkNTrackerHits = 0;
0930   int glbNMuonHits = 0;
0931   int staNMuonHits = 0;
0932   int NTrackerHits = 0;
0933   int NMuonHits = 0;
0934 
0935   // Analyzer reco::Muon
0936   for (View<Muon>::const_iterator iMuon = muonColl.begin(); iMuon != muonColl.end(); ++iMuon) {
0937     double muonP, muonPt, muonEta, muonPhi;
0938     if (usePFMuon_) {
0939       muonP = iMuon->pfP4().P();
0940       muonPt = iMuon->pfP4().Pt();
0941       muonEta = iMuon->pfP4().Eta();
0942       muonPhi = iMuon->pfP4().Phi();
0943     } else {
0944       muonP = iMuon->p();
0945       muonPt = iMuon->pt();
0946       muonEta = iMuon->eta();
0947       muonPhi = iMuon->phi();
0948     }
0949 
0950     //histograms for fractions
0951     commonME_->hMuonAllP_->Fill(muonP);
0952     commonME_->hMuonAllPt_->Fill(muonPt);
0953     commonME_->hMuonAllEta_->Fill(muonEta);
0954     commonME_->hMuonAllPhi_->Fill(muonPhi);
0955 
0956     if (!selector_(*iMuon))
0957       continue;
0958     if (wantTightMuon_) {
0959       if (!muon::isTightMuon(*iMuon, thePrimaryVertex))
0960         continue;
0961     }
0962 
0963     TrackRef Track = iMuon->track();
0964 
0965     if (Track.isNonnull()) {
0966       commonME_->hMuonTrackP_->Fill(Track->p());
0967       commonME_->hMuonTrackPt_->Fill(Track->pt());
0968       commonME_->hMuonTrackEta_->Fill(Track->eta());
0969       commonME_->hMuonTrackPhi_->Fill(Track->phi());
0970 
0971       //ip histograms
0972       commonME_->hMuonTrackDxy_->Fill(Track->dxy());
0973       commonME_->hMuonTrackDz_->Fill(Track->dz());
0974     }
0975 
0976     if (iMuon->isGlobalMuon()) {
0977       Track = iMuon->combinedMuon();
0978       glbNTrackerHits = countTrackerHits(*Track);
0979       glbNMuonHits = countMuonHits(*Track);
0980     } else if (iMuon->isTrackerMuon()) {
0981       Track = iMuon->track();
0982       trkNTrackerHits = countTrackerHits(*Track);
0983     } else {
0984       Track = iMuon->standAloneMuon();
0985     }
0986 
0987     NTrackerHits = countTrackerHits(*Track);
0988     muonME_->hNTrackerHits_->Fill(NTrackerHits);
0989     muonME_->hNTrackerHits_vs_Pt_->Fill(Track->pt(), NTrackerHits);
0990     muonME_->hNTrackerHits_vs_Eta_->Fill(Track->eta(), NTrackerHits);
0991 
0992     NMuonHits = countMuonHits(*Track);
0993     muonME_->hNMuonHits_->Fill(NMuonHits);
0994     muonME_->hNMuonHits_vs_Pt_->Fill(Track->pt(), NMuonHits);
0995     muonME_->hNMuonHits_vs_Eta_->Fill(Track->eta(), NMuonHits);
0996 
0997     //list of histos for each type
0998 
0999     //      muonME_->hNTrks_->Fill();
1000     muonME_->hNTrksEta_->Fill(Track->eta());
1001     muonME_->hNTrksPt_->Fill(Track->pt());
1002 
1003     commonME_->hMuonP_->Fill(muonP);
1004     commonME_->hMuonPt_->Fill(muonPt);
1005     commonME_->hMuonEta_->Fill(muonEta);
1006     commonME_->hMuonPhi_->Fill(muonPhi);
1007 
1008     if (iMuon->isGlobalMuon()) {
1009       double gtHitPat = iMuon->globalTrack()->hitPattern().numberOfAllHits(HitPattern::TRACK_HITS) -
1010                         iMuon->globalTrack()->hitPattern().numberOfValidHits();
1011 
1012       double itHitPat = iMuon->innerTrack()->hitPattern().numberOfAllHits(HitPattern::TRACK_HITS) -
1013                         iMuon->innerTrack()->hitPattern().numberOfValidHits();
1014 
1015       double otHitPat = iMuon->outerTrack()->hitPattern().numberOfAllHits(HitPattern::TRACK_HITS) -
1016                         iMuon->outerTrack()->hitPattern().numberOfValidHits();
1017 
1018       commonME_->hNInvalidHitsGTHitPattern_->Fill(gtHitPat);
1019       commonME_->hNInvalidHitsITHitPattern_->Fill(itHitPat);
1020       commonME_->hNInvalidHitsOTHitPattern_->Fill(otHitPat);
1021       commonME_->hNDeltaInvalidHitsHitPattern_->Fill(gtHitPat - itHitPat - otHitPat);
1022 
1023       //must be global and standalone
1024       if (iMuon->isStandAloneMuon()) {
1025         commonME_->hStaToGlbDiffNMuonHitsEta_->Fill(Track->eta(), staNMuonHits - glbNMuonHits);
1026         commonME_->hStaToGlbDiffNMuonHitsPt_->Fill(Track->pt(), staNMuonHits - glbNMuonHits);
1027         commonME_->hStaToGlbDiffNMuonHits_->Fill(staNMuonHits - glbNMuonHits);
1028       }
1029 
1030       //must be global and tracker
1031       if (iMuon->isTrackerMuon()) {
1032         commonME_->hTrkToGlbDiffNTrackerHitsEta_->Fill(Track->eta(), trkNTrackerHits - glbNTrackerHits);
1033         commonME_->hTrkToGlbDiffNTrackerHitsPt_->Fill(Track->pt(), trkNTrackerHits - glbNTrackerHits);
1034         commonME_->hTrkToGlbDiffNTrackerHits_->Fill(trkNTrackerHits - glbNTrackerHits);
1035       }
1036     }
1037 
1038   }  //end of reco muon loop
1039 
1040   // Associate by hits
1041   for (size_t i = 0; i < nSim; i++) {
1042     TrackingParticleRef simRef = TPrefV[i];
1043     const TrackingParticle* simTP = simRef.get();
1044     if (!tpSelector_(*simTP))
1045       continue;
1046 
1047     //denominators for efficiency plots
1048     const double simP = simRef->p();
1049     const double simPt = simRef->pt();
1050     const double simEta = doAbsEta_ ? fabs(simRef->eta()) : simRef->eta();
1051     const double simPhi = simRef->phi();
1052 
1053     GlobalPoint simVtx(simRef->vertex().x(), simRef->vertex().y(), simRef->vertex().z());
1054     GlobalVector simMom(simRef->momentum().x(), simRef->momentum().y(), simRef->momentum().z());
1055     const double simDxy = -simVtx.x() * sin(simPhi) + simVtx.y() * cos(simPhi);
1056     const double simDz = simVtx.z() - (simVtx.x() * simMom.x() + simVtx.y() * simMom.y()) * simMom.z() / simMom.perp2();
1057 
1058     const unsigned int nSimHits = simRef->numberOfHits();
1059 
1060     muonME_->hSimP_->Fill(simP);
1061     muonME_->hSimPt_->Fill(simPt);
1062     muonME_->hSimEta_->Fill(simEta);
1063     muonME_->hSimPhi_->Fill(simPhi);
1064     muonME_->hSimDxy_->Fill(simDxy);
1065     muonME_->hSimDz_->Fill(simDz);
1066     muonME_->hNSimHits_->Fill(nSimHits);
1067 
1068     // Get sim-reco association for a simRef
1069     vector<pair<RefToBase<Muon>, double> > MuRefV;
1070     if (simToMuonColl.find(simRef) != simToMuonColl.end()) {
1071       MuRefV = simToMuonColl[simRef];
1072 
1073       if (!MuRefV.empty()) {
1074         muonME_->hNSimToReco_->Fill(MuRefV.size());
1075         const Muon* Mu = MuRefV.begin()->first.get();
1076         if (!selector_(*Mu))
1077           continue;
1078         if (wantTightMuon_) {
1079           if (!muon::isTightMuon(*Mu, thePrimaryVertex))
1080             continue;
1081         }
1082 
1083         muonME_->fill(&*simTP, Mu);
1084       }
1085     }
1086   }
1087 }
1088 
1089 int RecoMuonValidator::countMuonHits(const reco::Track& track) const {
1090   TransientTrackingRecHit::ConstRecHitContainer result;
1091 
1092   int count = 0;
1093 
1094   for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
1095     if ((*hit)->isValid()) {
1096       DetId recoid = (*hit)->geographicalId();
1097       if (recoid.det() == DetId::Muon)
1098         count++;
1099     }
1100   }
1101   return count;
1102 }
1103 
1104 int RecoMuonValidator::countTrackerHits(const reco::Track& track) const {
1105   TransientTrackingRecHit::ConstRecHitContainer result;
1106 
1107   int count = 0;
1108 
1109   for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
1110     if ((*hit)->isValid()) {
1111       DetId recoid = (*hit)->geographicalId();
1112       if (recoid.det() == DetId::Tracker)
1113         count++;
1114     }
1115   }
1116   return count;
1117 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!