Reference Manual

Go to the documentation of this file.
00001 #include "Validation/RecoMuon/src/RecoMuonValidator.h"
00002 
00003 #include "DataFormats/TrackReco/interface/Track.h"
00004 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00005 #include "DataFormats/MuonReco/interface/Muon.h"
00006 #include "DataFormats/MuonReco/interface/MuonFwd.h"
00007 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
00008 
00009 #include "SimTracker/Records/interface/TrackAssociatorRecord.h"
00010 #include "SimTracker/TrackAssociation/interface/TrackAssociatorByChi2.h"
00011 #include "SimTracker/TrackAssociation/interface/TrackAssociatorByHits.h"
00012 
00013 #include "SimMuon/MCTruth/interface/MuonAssociatorByHits.h"
00014 
00015 #include "FWCore/ServiceRegistry/interface/Service.h"
00016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00017 #include "DQMServices/Core/interface/DQMStore.h"
00018 #include "DQMServices/Core/interface/MonitorElement.h"
00019 
00020 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00021 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
00022 
00023 // for selection cut
00024 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
00025 
00026 #include "TMath.h"
00027 
00028 using namespace std;
00029 using namespace edm;
00030 using namespace reco;
00031 
00032 typedef TrajectoryStateOnSurface TSOS;
00033 typedef FreeTrajectoryState FTS;
00034 
00035 //
00036 //struct for histogram dimensions
00037 //
00038 struct HistoDimensions {
00039   //p
00040   unsigned int nBinP;
00041   double minP, maxP;
00042   //pt
00043   unsigned int nBinPt;
00044   double minPt, maxPt;
00045   //if abs eta
00046   bool doAbsEta;
00047   //eta
00048   unsigned int nBinEta;
00049   double minEta, maxEta;
00050   //phi
00051   unsigned int nBinPhi;
00052   double minPhi, maxPhi;
00053   //dxy
00054   unsigned int nBinDxy;
00055   double minDxy, maxDxy;
00056   //dz
00057   unsigned int nBinDz; 
00058   double minDz, maxDz;
00059   //pulls
00060   unsigned int nBinPull;
00061   double wPull;
00062   //resolustions
00063   unsigned int nBinErr;
00064   double minErrP, maxErrP;
00065   double minErrPt, maxErrPt;
00066   double minErrQPt, maxErrQPt;
00067   double minErrEta, maxErrEta;
00068   double minErrPhi, maxErrPhi;
00069   double minErrDxy, maxErrDxy;
00070   double minErrDz, maxErrDz;
00071   //track multiplicities
00072   unsigned int nTrks, nAssoc;
00073   unsigned int nDof;
00074   // for PF muons
00075   bool usePFMuon;
00076 };
00077 
00078 //
00079 //Struct containing all histograms definitions
00080 //
00081 struct RecoMuonValidator::MuonME {
00082   typedef MonitorElement* MEP;
00083 
00084 //general kinematics
00085   MEP hSimP_, hSimPt_, hSimEta_, hSimPhi_, hSimDxy_, hSimDz_;
00086 //only for efficiencies
00087   MEP hP_, hPt_, hEta_, hPhi_;
00088   MEP hNSim_, hNMuon_;
00089 
00090 //misc vars
00091   MEP hNTrks_, hNTrksEta_,  hNTrksPt_;
00092   MEP hMisQPt_, hMisQEta_;
00093 
00094 //resolutions
00095   MEP hErrP_, hErrPt_, hErrEta_, hErrPhi_;
00096   MEP hErrPBarrel_, hErrPOverlap_, hErrPEndcap_;
00097   MEP hErrPtBarrel_, hErrPtOverlap_, hErrPtEndcap_;
00098   MEP hErrDxy_, hErrDz_;
00099 
00100   MEP hErrP_vs_Eta_, hErrPt_vs_Eta_, hErrQPt_vs_Eta_;
00101   MEP hErrP_vs_P_, hErrPt_vs_Pt_, hErrQPt_vs_Pt_, hErrEta_vs_Eta_;
00102 
00103 //PF-RECO event-by-event comparisons
00104   MEP hErrPt_PF_;
00105   MEP hErrQPt_PF_;
00106   MEP hdPt_vs_Eta_;
00107   MEP hdPt_vs_Pt_;
00108   MEP hPFMomAssCorrectness;
00109   MEP hPt_vs_PFMomAssCorrectness;
00110 
00111 //hit pattern
00112   MEP hNSimHits_;
00113   MEP hNSimToReco_, hNRecoToSim_;
00114 
00115   MEP hNHits_, hNLostHits_, hNTrackerHits_, hNMuonHits_;
00116   MEP hNHits_vs_Pt_, hNHits_vs_Eta_;
00117   MEP hNLostHits_vs_Pt_, hNLostHits_vs_Eta_;
00118   MEP hNTrackerHits_vs_Pt_, hNTrackerHits_vs_Eta_;
00119   MEP hNMuonHits_vs_Pt_, hNMuonHits_vs_Eta_;
00120 
00121 //pulls
00122   MEP hPullPt_, hPullEta_, hPullPhi_, hPullQPt_, hPullDxy_, hPullDz_;
00123   MEP hPullPt_vs_Eta_, hPullPt_vs_Pt_, hPullEta_vs_Eta_, hPullPhi_vs_Eta_, hPullEta_vs_Pt_;
00124 
00125 //chi2, ndof
00126   MEP hNDof_, hChi2_, hChi2Norm_, hChi2Prob_;
00127   MEP hNDof_vs_Eta_, hChi2_vs_Eta_, hChi2Norm_vs_Eta_, hChi2Prob_vs_Eta_;
00128   
00129   bool doAbsEta_;
00130   bool usePFMuon_;
00131 
00132 
00133 //
00134 //books histograms
00135 //
00136   void bookHistograms(DQMStore* dqm, const string& dirName, const HistoDimensions& hDim)
00137   {
00138     dqm->cd();
00139     dqm->setCurrentFolder(dirName.c_str());
00140 
00141     doAbsEta_ = hDim.doAbsEta;
00142     usePFMuon_ = hDim.usePFMuon;
00143 
00144     //histograms for efficiency plots
00145     hP_   = dqm->book1D("P"  , "p of recoTracks"    , hDim.nBinP  , hDim.minP  , hDim.maxP  );
00146     hPt_  = dqm->book1D("Pt" , "p_{T} of recoTracks", hDim.nBinPt , hDim.minPt , hDim.maxPt );
00147     hEta_ = dqm->book1D("Eta", "#eta of recoTracks" , hDim.nBinEta, hDim.minEta, hDim.maxEta);
00148     hPhi_ = dqm->book1D("Phi", "#phi of recoTracks" , hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
00149 
00150     hSimP_   = dqm->book1D("SimP"  , "p of simTracks"    , hDim.nBinP  , hDim.minP  , hDim.maxP  );
00151     hSimPt_  = dqm->book1D("SimPt" , "p_{T} of simTracks", hDim.nBinPt , hDim.minPt , hDim.maxPt );
00152     hSimEta_ = dqm->book1D("SimEta", "#eta of simTracks" , hDim.nBinEta, hDim.minEta, hDim.maxEta);
00153     hSimPhi_ = dqm->book1D("SimPhi", "#phi of simTracks" , hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
00154     hSimDxy_ = dqm->book1D("SimDxy", "Dxy of simTracks" , hDim.nBinDxy, hDim.minDxy, hDim.maxDxy);
00155     hSimDz_ = dqm->book1D("Dz", "Dz of simTracks" , hDim.nBinDz, hDim.minDz, hDim.maxDz);
00156 
00157     //track multiplicities
00158     hNSim_  = dqm->book1D("NSim" , "Number of particles per event", hDim.nTrks, -0.5, hDim.nTrks+0.5);
00159     hNMuon_ = dqm->book1D("NMuon", "Number of muons per event"    , hDim.nTrks, -0.5, hDim.nTrks+0.5);
00160 
00161     // - Misc. variables
00162     hNTrks_ = dqm->book1D("NTrks", "Number of reco tracks per event", hDim.nTrks, -0.5, hDim.nTrks+0.5);
00163     hNTrksEta_ = dqm->book1D("NTrksEta", "Number of reco tracks vs #eta", hDim.nBinEta, hDim.minEta, hDim.maxEta);
00164     hNTrksPt_ = dqm->book1D("NTrksPt", "Number of reco tracks vs p_{T}", hDim.nBinPt, hDim.minPt, hDim.maxPt);
00165 
00166     hMisQPt_  = dqm->book1D("MisQPt" , "Charge mis-id vs Pt" , hDim.nBinPt , hDim.minPt , hDim.maxPt );
00167     hMisQEta_ = dqm->book1D("MisQEta", "Charge mis-id vs Eta", hDim.nBinEta, hDim.minEta, hDim.maxEta);
00168 
00169     // - Resolutions
00170     hErrP_   = dqm->book1D("ErrP"  , "#Delta(p)/p"        , hDim.nBinErr, hDim.minErrP  , hDim.maxErrP  );
00171     hErrPBarrel_   = dqm->book1D("ErrP_barrel"  , "#Delta(p)/p"        , hDim.nBinErr, hDim.minErrP  , hDim.maxErrP  );
00172     hErrPOverlap_   = dqm->book1D("ErrP_overlap"  , "#Delta(p)/p"        , hDim.nBinErr, hDim.minErrP  , hDim.maxErrP  );
00173     hErrPEndcap_   = dqm->book1D("ErrP_endcap"  , "#Delta(p)/p"        , hDim.nBinErr, hDim.minErrP  , hDim.maxErrP  );
00174     hErrPt_  = dqm->book1D("ErrPt" , "#Delta(p_{T})/p_{T}", hDim.nBinErr, hDim.minErrPt , hDim.maxErrPt );
00175     hErrPtBarrel_  = dqm->book1D("ErrPt_barrel" , "#Delta(p_{T})/p_{T}", hDim.nBinErr, hDim.minErrPt , hDim.maxErrPt );
00176     hErrPtOverlap_  = dqm->book1D("ErrPt_overlap" , "#Delta(p_{T})/p_{T}", hDim.nBinErr, hDim.minErrPt , hDim.maxErrPt );
00177     hErrPtEndcap_  = dqm->book1D("ErrPt_endcap" , "#Delta(p_{T})/p_{T}", hDim.nBinErr, hDim.minErrPt , hDim.maxErrPt );
00178     hErrEta_ = dqm->book1D("ErrEta", "#sigma(#eta))"      , hDim.nBinErr, hDim.minErrEta, hDim.maxErrEta);
00179     hErrPhi_ = dqm->book1D("ErrPhi", "#sigma(#phi)"       , hDim.nBinErr, hDim.minErrPhi, hDim.maxErrPhi);
00180     hErrDxy_ = dqm->book1D("ErrDxy", "#sigma(d_{xy})"     , hDim.nBinErr, hDim.minErrDxy, hDim.maxErrDxy);
00181     hErrDz_  = dqm->book1D("ErrDz" , "#sigma(d_{z})"      , hDim.nBinErr, hDim.minErrDz , hDim.maxErrDz );
00182 
00183     //PF-RECO comparisons
00184     if (usePFMuon_) {
00185       hErrPt_PF_  = dqm->book1D("ErrPt_PF" , "#Delta(p_{T})|_{PF}/p_{T}", hDim.nBinErr, hDim.minErrPt, hDim.maxErrPt );
00186       hErrQPt_PF_  = dqm->book1D("ErrQPt_PF" , "#Delta(q/p_{T})|_{PF}/(q/p_{T})", hDim.nBinErr, hDim.minErrQPt, hDim.maxErrQPt);
00187       
00188       hPFMomAssCorrectness  = dqm->book1D("hPFMomAssCorrectness", "Corrected momentum assignement PF/RECO",2,0.5,2.5);
00189       hPt_vs_PFMomAssCorrectness  = dqm->book2D("hPt_vs_PFMomAssCorrectness", "Corrected momentum assignement PF/RECO", hDim.nBinPt, hDim.minPt, hDim.maxP, 2, 0.5, 2.5);
00190       
00191       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);
00192       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);
00193     }
00194 
00195     // -- Resolutions vs Eta
00196     hErrP_vs_Eta_   = dqm->book2D("ErrP_vs_Eta", "#Delta(p)/p vs #eta",
00197                                   hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinErr, hDim.minErrP, hDim.maxErrP);
00198     hErrPt_vs_Eta_  = dqm->book2D("ErrPt_vs_Eta", "#Delta(p_{T})/p_{T} vs #eta",
00199                                   hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinErr, hDim.minErrPt, hDim.maxErrPt);
00200     hErrQPt_vs_Eta_ = dqm->book2D("ErrQPt_vs_Eta", "#Delta(q/p_{T})/(q/p_{T}) vs #eta",
00201                                   hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinErr, hDim.minErrQPt, hDim.maxErrQPt);
00202     hErrEta_vs_Eta_ = dqm->book2D("ErrEta_vs_Eta", "#sigma(#eta) vs #eta",
00203                                   hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinErr, hDim.minErrEta, hDim.maxErrEta);
00204 
00205     // -- Resolutions vs momentum
00206     hErrP_vs_P_    = dqm->book2D("ErrP_vs_P", "#Delta(p)/p vs p",
00207                                  hDim.nBinP, hDim.minP, hDim.maxP, hDim.nBinErr, hDim.minErrP, hDim.maxErrP);
00208     hErrPt_vs_Pt_  = dqm->book2D("ErrPt_vs_Pt", "#Delta(p_{T})/p_{T} vs p_{T}",
00209                                  hDim.nBinPt, hDim.minPt, hDim.maxPt, hDim.nBinErr, hDim.minErrPt, hDim.maxErrPt);
00210     hErrQPt_vs_Pt_ = dqm->book2D("ErrQPt_vs_Pt", "#Delta(q/p_{T})/(q/p_{T}) vs p_{T}",
00211                                  hDim.nBinPt, hDim.minPt, hDim.maxPt, hDim.nBinErr, hDim.minErrQPt, hDim.maxErrQPt);
00212 
00213     // - Pulls
00214     hPullPt_  = dqm->book1D("PullPt" , "Pull(#p_{T})" , hDim.nBinPull, -hDim.wPull, hDim.wPull);
00215     hPullEta_ = dqm->book1D("PullEta", "Pull(#eta)"   , hDim.nBinPull, -hDim.wPull, hDim.wPull);
00216     hPullPhi_ = dqm->book1D("PullPhi", "Pull(#phi)"   , hDim.nBinPull, -hDim.wPull, hDim.wPull);
00217     hPullQPt_ = dqm->book1D("PullQPt", "Pull(q/p_{T})", hDim.nBinPull, -hDim.wPull, hDim.wPull);
00218     hPullDxy_ = dqm->book1D("PullDxy", "Pull(D_{xy})" , hDim.nBinPull, -hDim.wPull, hDim.wPull);
00219     hPullDz_  = dqm->book1D("PullDz" , "Pull(D_{z})"  , hDim.nBinPull, -hDim.wPull, hDim.wPull);
00220 
00221     // -- Pulls vs Eta
00222     hPullPt_vs_Eta_  = dqm->book2D("PullPt_vs_Eta", "Pull(p_{T}) vs #eta",
00223                                    hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinPull, -hDim.wPull, hDim.wPull);
00224     hPullEta_vs_Eta_ = dqm->book2D("PullEta_vs_Eta", "Pull(#eta) vs #eta",
00225                                    hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinPull, -hDim.wPull, hDim.wPull);
00226     hPullPhi_vs_Eta_ = dqm->book2D("PullPhi_vs_Eta", "Pull(#phi) vs #eta",
00227                                    hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinPull, -hDim.wPull, hDim.wPull);
00228 
00229     // -- Pulls vs Pt
00230     hPullPt_vs_Pt_ = dqm->book2D("PullPt_vs_Pt", "Pull(p_{T}) vs p_{T}",
00231                                  hDim.nBinPt, hDim.minPt, hDim.maxPt, hDim.nBinPull, -hDim.wPull, hDim.wPull);
00232     hPullEta_vs_Pt_ = dqm->book2D("PullEta_vs_Pt", "Pull(#eta) vs p_{T}",
00233                                   hDim.nBinPt, hDim.minPt, hDim.maxPt, hDim.nBinPull, -hDim.wPull, hDim.wPull);
00234 
00235     // -- Number of Hits
00236     const int nHits = 100;
00237     hNHits_ = dqm->book1D("NHits", "Number of hits", nHits+1, -0.5, nHits+0.5);
00238     hNHits_vs_Pt_  = dqm->book2D("NHits_vs_Pt", "Number of hits vs p_{T}",
00239                                  hDim.nBinPt, hDim.minPt, hDim.maxPt, nHits/4+1, -0.25, nHits+0.25);
00240     hNHits_vs_Eta_ = dqm->book2D("NHits_vs_Eta", "Number of hits vs #eta",
00241                                  hDim.nBinEta, hDim.minEta, hDim.maxEta, nHits/4+1, -0.25, nHits+0.25);
00242     hNSimHits_ = dqm->book1D("NSimHits", "Number of simHits", nHits+1, -0.5, nHits+0.5);
00243 
00244     const int nLostHits = 5;
00245     hNLostHits_ = dqm->book1D("NLostHits", "Number of Lost hits", nLostHits+1, -0.5, nLostHits+0.5);
00246     hNLostHits_vs_Pt_  = dqm->book2D("NLostHits_vs_Pt", "Number of lost Hits vs p_{T}",
00247                                      hDim.nBinPt, hDim.minPt, hDim.maxPt, nLostHits+1, -0.5, nLostHits+0.5);
00248     hNLostHits_vs_Eta_ = dqm->book2D("NLostHits_vs_Eta", "Number of lost Hits vs #eta",
00249                                      hDim.nBinEta, hDim.minEta, hDim.maxEta, nLostHits+1, -0.5, nLostHits+0.5);
00250 
00251     const int nTrackerHits = 40;
00252     hNTrackerHits_ = dqm->book1D("NTrackerHits", "Number of valid tracker hits", nTrackerHits+1, -0.5, nTrackerHits+0.5);
00253     hNTrackerHits_vs_Pt_ = dqm->book2D("NTrackerHits_vs_Pt", "Number of valid traker hits vs p_{T}",
00254                                        hDim.nBinPt, hDim.minPt, hDim.maxPt, nTrackerHits/4+1, -0.25, nTrackerHits+0.25);
00255     hNTrackerHits_vs_Eta_ = dqm->book2D("NTrackerHits_vs_Eta", "Number of valid tracker hits vs #eta",
00256                                         hDim.nBinEta, hDim.minEta, hDim.maxEta, nTrackerHits/4+1, -0.25, nTrackerHits+0.25);
00257 
00258     const int nMuonHits = 60;
00259     hNMuonHits_ = dqm->book1D("NMuonHits", "Number of valid muon hits", nMuonHits+1, -0.5, nMuonHits+0.5);
00260     hNMuonHits_vs_Pt_  = dqm->book2D("NMuonHits_vs_Pt", "Number of valid muon hits vs p_{T}",
00261                                      hDim.nBinPt, hDim.minPt, hDim.maxPt, nMuonHits/4+1, -0.25, nMuonHits+0.25);
00262     hNMuonHits_vs_Eta_ = dqm->book2D("NMuonHits_vs_Eta", "Number of valid muon hits vs #eta",
00263                                      hDim.nBinEta, hDim.minEta, hDim.maxEta, nMuonHits/4+1, -0.25, nMuonHits+0.25);
00264 
00265     hNDof_ = dqm->book1D("NDof", "Number of DoF", hDim.nDof+1, -0.5, hDim.nDof+0.5);
00266     hChi2_ = dqm->book1D("Chi2", "#Chi^{2}", hDim.nBinErr, 0, 200);
00267     hChi2Norm_ = dqm->book1D("Chi2Norm", "Normalized #Chi^{2}", hDim.nBinErr, 0, 50);
00268     hChi2Prob_ = dqm->book1D("Chi2Prob", "Prob(#Chi^{2})", hDim.nBinErr, 0, 1);
00269 
00270     hNDof_vs_Eta_ = dqm->book2D("NDof_vs_Eta", "Number of DoF vs #eta",
00271                                 hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nDof+1, -0.5, hDim.nDof+0.5);
00272     hChi2_vs_Eta_ = dqm->book2D("Chi2_vs_Eta", "#Chi^{2} vs #eta",
00273                                 hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinErr, 0., 200.);
00274     hChi2Norm_vs_Eta_ = dqm->book2D("Chi2Norm_vs_Eta", "Normalized #Chi^{2} vs #eta",
00275                                     hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinErr, 0., 50.);
00276     hChi2Prob_vs_Eta_ = dqm->book2D("Chi2Prob_vs_Eta", "Prob(#Chi^{2}) vs #eta",
00277                                     hDim.nBinEta, hDim.minEta, hDim.maxEta, hDim.nBinErr, 0., 1.);
00278 
00279     hNSimToReco_ = dqm->book1D("NSimToReco", "Number of associated reco tracks", hDim.nAssoc+1, -0.5, hDim.nAssoc+0.5);
00280     hNRecoToSim_ = dqm->book1D("NRecoToSim", "Number of associated sim TP's", hDim.nAssoc+1, -0.5, hDim.nAssoc+0.5);
00281 
00282   };
00283 
00284 //
00285 //Fill hists booked, simRef and muonRef are associated by hits
00286 //
00287   void fill(const TrackingParticle* simRef, const Muon* muonRef)
00288   {
00289 
00290     const double simP   = simRef->p();
00291     const double simPt  = simRef->pt();
00292     const double simEta = doAbsEta_ ? fabs(simRef->eta()) : simRef->eta();
00293     const double simPhi = simRef->phi();
00294     const double simQ   = simRef->charge();
00295     const double simQPt = simQ/simPt;
00296 
00297     GlobalPoint  simVtx(simRef->vertex().x(), simRef->vertex().y(), simRef->vertex().z());
00298     GlobalVector simMom(simRef->momentum().x(), simRef->momentum().y(), simRef->momentum().z());
00299     const double simDxy = -simVtx.x()*sin(simPhi)+simVtx.y()*cos(simPhi);
00300     const double simDz  = simVtx.z() - (simVtx.x()*simMom.x()+simVtx.y()*simMom.y())*simMom.z()/simMom.perp2();
00301 
00302     const double recoQ   = muonRef->charge();
00303     if ( simQ*recoQ < 0 ) {
00304       hMisQPt_ ->Fill(simPt );
00305       hMisQEta_->Fill(simEta);
00306     }
00307 
00308     double recoP, recoPt, recoEta, recoPhi, recoQPt;
00309     if (usePFMuon_) {
00310       //      const double origRecoP   = muonRef->p(); 
00311       const double origRecoPt  = muonRef->pt();
00312       //      const double origRecoEta = muonRef->eta();
00313       //      const double origRecoPhi = muonRef->phi();
00314       const double origRecoQPt = recoQ/origRecoPt;
00315       recoP = muonRef->pfP4().P();
00316       recoPt = muonRef->pfP4().Pt();
00317       recoEta = muonRef->pfP4().Eta(); 
00318       recoPhi = muonRef->pfP4().Phi();
00319       recoQPt = recoQ/recoPt;
00320       hErrPt_PF_->Fill((recoPt-origRecoPt)/origRecoPt);
00321       hErrQPt_PF_->Fill((recoQPt-origRecoQPt)/origRecoQPt);
00322 
00323       hdPt_vs_Eta_->Fill(recoEta,recoPt-origRecoPt);
00324       hdPt_vs_Pt_->Fill(recoPt,recoPt-origRecoPt);
00325 
00326       int theCorrectPFAss = (fabs(recoPt-simPt) < fabs(origRecoPt - simPt))? 1 : 2;
00327       hPFMomAssCorrectness->Fill(theCorrectPFAss);
00328       hPt_vs_PFMomAssCorrectness->Fill(simPt,theCorrectPFAss);
00329     }
00330   
00331     else {
00332       recoP   = muonRef->p(); 
00333       recoPt  = muonRef->pt();
00334       recoEta = muonRef->eta();
00335       recoPhi = muonRef->phi();
00336       recoQPt = recoQ/recoPt;
00337     }
00338 
00339     const double errP   = (recoP-simP)/simP;
00340     const double errPt  = (recoPt-simPt)/simPt;
00341     const double errEta = (recoEta-simEta)/simEta;
00342     const double errPhi = (recoPhi-simPhi)/simPhi;
00343     const double errQPt = (recoQPt-simQPt)/simQPt;
00344 
00345     hP_  ->Fill(simP);
00346     hPt_ ->Fill(simPt );
00347     hEta_->Fill(simEta);   
00348     hPhi_->Fill(simPhi);
00349 
00350     hErrP_  ->Fill(errP  );
00351     hErrPt_ ->Fill(errPt );
00352     hErrEta_->Fill(errEta);
00353     hErrPhi_->Fill(errPhi);
00354 
00355     if(fabs(simEta) > 0. && fabs(simEta) < 0.8) {
00356       hErrPBarrel_->Fill(errP);
00357       hErrPtBarrel_->Fill(errPt);
00358     } else if (fabs(simEta) > 0.8 && fabs(simEta) < 1.2) {
00359       hErrPOverlap_->Fill(errP);
00360       hErrPtOverlap_->Fill(errPt);
00361     } else if (fabs(simEta) > 1.2 ){
00362       hErrPEndcap_->Fill(errP);
00363       hErrPtEndcap_->Fill(errPt);
00364     }
00365 
00366     hErrP_vs_Eta_  ->Fill(simEta, errP  );
00367     hErrPt_vs_Eta_ ->Fill(simEta, errPt );
00368     hErrQPt_vs_Eta_->Fill(simEta, errQPt);
00369 
00370     hErrP_vs_P_   ->Fill(simP  , errP  );
00371     hErrPt_vs_Pt_ ->Fill(simPt , errPt );
00372     hErrQPt_vs_Pt_->Fill(simQPt, errQPt);
00373 
00374     hErrEta_vs_Eta_->Fill(simEta, errEta);
00375 
00376     //access from track
00377     reco::TrackRef recoRef = muonRef->track();
00378     if (recoRef.isNonnull()) {
00379 
00380     // Number of reco-hits
00381     const int nRecoHits = recoRef->numberOfValidHits();
00382     const int nLostHits = recoRef->numberOfLostHits();
00383 
00384     hNHits_->Fill(nRecoHits);
00385     hNHits_vs_Pt_ ->Fill(simPt , nRecoHits);
00386     hNHits_vs_Eta_->Fill(simEta, nRecoHits);
00387 
00388     hNLostHits_->Fill(nLostHits);
00389     hNLostHits_vs_Pt_ ->Fill(simPt , nLostHits);
00390     hNLostHits_vs_Eta_->Fill(simEta, nLostHits);
00391 
00392     const double recoNDof = recoRef->ndof();
00393     const double recoChi2 = recoRef->chi2();
00394     const double recoChi2Norm = recoRef->normalizedChi2();
00395     const double recoChi2Prob = TMath::Prob(recoRef->chi2(), static_cast<int>(recoRef->ndof()));
00396 
00397     hNDof_->Fill(recoNDof);
00398     hChi2_->Fill(recoChi2);
00399     hChi2Norm_->Fill(recoChi2Norm);
00400     hChi2Prob_->Fill(recoChi2Prob);
00401 
00402     hNDof_vs_Eta_->Fill(simEta, recoNDof);
00403     hChi2_vs_Eta_->Fill(simEta, recoChi2);
00404     hChi2Norm_vs_Eta_->Fill(simEta, recoChi2Norm);
00405     hChi2Prob_vs_Eta_->Fill(simEta, recoChi2Prob);
00406 
00407     const double recoDxy = recoRef->dxy();
00408     const double recoDz  = recoRef->dz();
00409 
00410     const double errDxy = (recoDxy-simDxy)/simDxy;
00411     const double errDz  = (recoDz-simDz)/simDz;
00412     hErrDxy_->Fill(errDxy);
00413     hErrDz_ ->Fill(errDz );
00414 
00415     const double pullPt  = (recoPt-simPt)/recoRef->ptError();
00416     const double pullQPt = (recoQPt-simQPt)/recoRef->qoverpError();
00417     const double pullEta = (recoEta-simEta)/recoRef->etaError();
00418     const double pullPhi = (recoPhi-simPhi)/recoRef->phiError();
00419     const double pullDxy = (recoDxy-simDxy)/recoRef->dxyError();
00420     const double pullDz  = (recoDz-simDz)/recoRef->dzError();
00421 
00422     hPullPt_ ->Fill(pullPt );
00423     hPullEta_->Fill(pullEta);
00424     hPullPhi_->Fill(pullPhi);
00425     hPullQPt_->Fill(pullQPt);
00426     hPullDxy_->Fill(pullDxy);
00427     hPullDz_ ->Fill(pullDz );
00428 
00429     hPullPt_vs_Eta_->Fill(simEta, pullPt);
00430     hPullPt_vs_Pt_ ->Fill(simPt, pullPt);
00431 
00432     hPullEta_vs_Eta_->Fill(simEta, pullEta);
00433     hPullPhi_vs_Eta_->Fill(simEta, pullPhi);
00434 
00435     hPullEta_vs_Pt_->Fill(simPt, pullEta);
00436     }
00437   };
00438 };
00439 
00440 //
00441 //struct defininiong histograms
00442 //
00443 struct RecoMuonValidator::CommonME {
00444   typedef MonitorElement* MEP;
00445 
00446 //diffs
00447   MEP hTrkToGlbDiffNTrackerHits_, hStaToGlbDiffNMuonHits_;
00448   MEP hTrkToGlbDiffNTrackerHitsEta_, hStaToGlbDiffNMuonHitsEta_;
00449   MEP hTrkToGlbDiffNTrackerHitsPt_, hStaToGlbDiffNMuonHitsPt_;
00450 
00451 //global muon hit pattern
00452   MEP hNInvalidHitsGTHitPattern_, hNInvalidHitsITHitPattern_, hNInvalidHitsOTHitPattern_;
00453   MEP hNDeltaInvalidHitsHitPattern_;
00454 
00455 //muon based momentum assignment
00456   MEP hMuonP_, hMuonPt_, hMuonEta_, hMuonPhi_;
00457 //track based kinematics
00458   MEP hMuonTrackP_, hMuonTrackPt_, hMuonTrackEta_, hMuonTrackPhi_, hMuonTrackDxy_, hMuonTrackDz_;
00459 //histograms for fractions
00460   MEP hMuonAllP_, hMuonAllPt_, hMuonAllEta_, hMuonAllPhi_;
00461 };
00462 
00463 //
00464 //Constructor
00465 //
00466 RecoMuonValidator::RecoMuonValidator(const ParameterSet& pset):
00467   selector_(pset.getParameter<std::string>("selection"))
00468 {
00469   verbose_ = pset.getUntrackedParameter<unsigned int>("verbose", 0);
00470 
00471   outputFileName_ = pset.getUntrackedParameter<string>("outputFileName", "");
00472 
00473 
00474   // Set histogram dimensions from config
00475   HistoDimensions hDim;
00476   
00477   hDim.nBinP = pset.getUntrackedParameter<unsigned int>("nBinP");
00478   hDim.minP = pset.getUntrackedParameter<double>("minP");
00479   hDim.maxP = pset.getUntrackedParameter<double>("maxP");
00480 
00481   hDim.nBinPt = pset.getUntrackedParameter<unsigned int>("nBinPt");
00482   hDim.minPt = pset.getUntrackedParameter<double>("minPt");
00483   hDim.maxPt = pset.getUntrackedParameter<double>("maxPt");
00484 
00485   doAbsEta_ = pset.getUntrackedParameter<bool>("doAbsEta");
00486   hDim.doAbsEta = doAbsEta_;
00487   hDim.nBinEta  = pset.getUntrackedParameter<unsigned int>("nBinEta");
00488   hDim.minEta = pset.getUntrackedParameter<double>("minEta");
00489   hDim.maxEta = pset.getUntrackedParameter<double>("maxEta");
00490 
00491   hDim.nBinDxy  = pset.getUntrackedParameter<unsigned int>("nBinDxy");
00492   hDim.minDxy = pset.getUntrackedParameter<double>("minDxy");
00493   hDim.maxDxy = pset.getUntrackedParameter<double>("maxDxy");
00494 
00495   hDim.nBinDz  = pset.getUntrackedParameter<unsigned int>("nBinDz");
00496   hDim.minDz = pset.getUntrackedParameter<double>("minDz");
00497   hDim.maxDz = pset.getUntrackedParameter<double>("maxDz");
00498 
00499   hDim.nBinPhi  = pset.getUntrackedParameter<unsigned int>("nBinPhi");
00500   hDim.minPhi = pset.getUntrackedParameter<double>("minPhi", -TMath::Pi());
00501   hDim.maxPhi = pset.getUntrackedParameter<double>("maxPhi",  TMath::Pi());
00502 
00503   hDim.nBinErr  = pset.getUntrackedParameter<unsigned int>("nBinErr");
00504   hDim.nBinPull = pset.getUntrackedParameter<unsigned int>("nBinPull");
00505 
00506   hDim.wPull = pset.getUntrackedParameter<double>("wPull");
00507 
00508   hDim.minErrP = pset.getUntrackedParameter<double>("minErrP");
00509   hDim.maxErrP = pset.getUntrackedParameter<double>("maxErrP");
00510 
00511   hDim.minErrPt = pset.getUntrackedParameter<double>("minErrPt");
00512   hDim.maxErrPt = pset.getUntrackedParameter<double>("maxErrPt");
00513 
00514   hDim.minErrQPt = pset.getUntrackedParameter<double>("minErrQPt");
00515   hDim.maxErrQPt = pset.getUntrackedParameter<double>("maxErrQPt");
00516 
00517   hDim.minErrEta = pset.getUntrackedParameter<double>("minErrEta");
00518   hDim.maxErrEta = pset.getUntrackedParameter<double>("maxErrEta");
00519 
00520   hDim.minErrPhi = pset.getUntrackedParameter<double>("minErrPhi");
00521   hDim.maxErrPhi = pset.getUntrackedParameter<double>("maxErrPhi");
00522 
00523   hDim.minErrDxy = pset.getUntrackedParameter<double>("minErrDxy");
00524   hDim.maxErrDxy = pset.getUntrackedParameter<double>("maxErrDxy");
00525 
00526   hDim.minErrDz  = pset.getUntrackedParameter<double>("minErrDz" );
00527   hDim.maxErrDz  = pset.getUntrackedParameter<double>("maxErrDz" );
00528 
00529   hDim.nTrks = pset.getUntrackedParameter<unsigned int>("nTrks");
00530   hDim.nAssoc = pset.getUntrackedParameter<unsigned int>("nAssoc");
00531   hDim.nDof = pset.getUntrackedParameter<unsigned int>("nDof", 55);
00532 
00533   // Labels for simulation and reconstruction tracks
00534   simLabel_  = pset.getParameter<InputTag>("simLabel" );
00535   muonLabel_ = pset.getParameter<InputTag>("muonLabel");
00536 
00537   // Labels for sim-reco association
00538   doAssoc_ = pset.getUntrackedParameter<bool>("doAssoc", true);
00539   muAssocLabel_ = pset.getParameter<InputTag>("muAssocLabel");
00540 
00541   // Different momentum assignment and additional histos in case of PF muons
00542   usePFMuon_ = pset.getUntrackedParameter<bool>("usePFMuon");
00543   hDim.usePFMuon = usePFMuon_;
00544 
00545   //type of track
00546   std::string trackType = pset.getParameter< std::string >("trackType");
00547   if (trackType == "inner") trackType_ = MuonAssociatorByHits::InnerTk;
00548   else if (trackType == "outer") trackType_ = MuonAssociatorByHits::OuterTk;
00549   else if (trackType == "global") trackType_ = MuonAssociatorByHits::GlobalTk;
00550   else if (trackType == "segments") trackType_ = MuonAssociatorByHits::Segments;
00551   else throw cms::Exception("Configuration") << "Track type '" << trackType << "' not supported.\n";
00552 
00553 //  seedPropagatorName_ = pset.getParameter<string>("SeedPropagator");
00554 
00555   ParameterSet tpset = pset.getParameter<ParameterSet>("tpSelector");
00556   tpSelector_ = TrackingParticleSelector(tpset.getParameter<double>("ptMin"),
00557                                          tpset.getParameter<double>("minRapidity"),
00558                                          tpset.getParameter<double>("maxRapidity"),
00559                                          tpset.getParameter<double>("tip"),
00560                                          tpset.getParameter<double>("lip"),
00561                                          tpset.getParameter<int>("minHit"),
00562                                          tpset.getParameter<bool>("signalOnly"),
00563                                          tpset.getParameter<bool>("chargedOnly"),
00564                                          tpset.getParameter<bool>("stableOnly"),
00565                                          tpset.getParameter<std::vector<int> >("pdgId"));
00566 
00567   // the service parameters
00568   ParameterSet serviceParameters 
00569     = pset.getParameter<ParameterSet>("ServiceParameters");
00570   theMuonService = new MuonServiceProxy(serviceParameters);
00571 
00572   // retrieve the instance of DQMService
00573   theDQM = 0;
00574   theDQM = Service<DQMStore>().operator->();
00575 
00576   if ( ! theDQM ) {
00577     LogError("RecoMuonValidator") << "DQMService not initialized\n";
00578     return;
00579   }
00580 
00581   subDir_ = pset.getUntrackedParameter<string>("subDir");
00582   if ( subDir_.empty() ) subDir_ = "RecoMuonV";
00583   if ( subDir_[subDir_.size()-1] == '/' ) subDir_.erase(subDir_.size()-1);
00584 
00585   // book histograms
00586   theDQM->cd();
00587 
00588   theDQM->setCurrentFolder(subDir_);
00589 
00590   commonME_ = new CommonME;
00591   muonME_  = new MuonME;
00592 
00593   //commonME
00594   const int nHits = 100;
00595 
00596   // - diffs
00597   commonME_->hTrkToGlbDiffNTrackerHits_ = theDQM->book1D("TrkGlbDiffNTrackerHits", "Difference of number of tracker hits (tkMuon - globalMuon)", 2*nHits+1, -nHits-0.5, nHits+0.5);
00598   commonME_->hStaToGlbDiffNMuonHits_ = theDQM->book1D("StaGlbDiffNMuonHits", "Difference of number of muon hits (staMuon - globalMuon)", 2*nHits+1, -nHits-0.5, nHits+0.5);
00599 
00600   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);
00601   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);
00602 
00603   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);
00604   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);
00605 
00606   // -global muon hit pattern
00607   commonME_->hNInvalidHitsGTHitPattern_ = theDQM->book1D("NInvalidHitsGTHitPattern", "Number of invalid hits on a global track", nHits+1, -0.5, nHits+0.5);
00608   commonME_->hNInvalidHitsITHitPattern_ = theDQM->book1D("NInvalidHitsITHitPattern", "Number of invalid hits on an inner track", nHits+1, -0.5, nHits+0.5);
00609   commonME_->hNInvalidHitsOTHitPattern_ = theDQM->book1D("NInvalidHitsOTHitPattern", "Number of invalid hits on an outer track", nHits+1, -0.5, nHits+0.5);
00610   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);
00611 
00612   //muon based kinematics
00613   commonME_->hMuonP_   = theDQM->book1D("PMuon"  , "p of muon"    , hDim.nBinP  , hDim.minP  , hDim.maxP  );
00614   commonME_->hMuonPt_  = theDQM->book1D("PtMuon" , "p_{T} of muon", hDim.nBinPt , hDim.minPt , hDim.maxPt );
00615   commonME_->hMuonEta_ = theDQM->book1D("EtaMuon", "#eta of muon" , hDim.nBinEta, hDim.minEta, hDim.maxEta);
00616   commonME_->hMuonPhi_ = theDQM->book1D("PhiMuon", "#phi of muon" , hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
00617   //track based kinematics
00618   commonME_->hMuonTrackP_   = theDQM->book1D("PMuonTrack"  , "p of reco muon track"    , hDim.nBinP  , hDim.minP  , hDim.maxP  );
00619   commonME_->hMuonTrackPt_  = theDQM->book1D("PtMuonTrack" , "p_{T} of reco muon track", hDim.nBinPt , hDim.minPt , hDim.maxPt );
00620   commonME_->hMuonTrackEta_ = theDQM->book1D("EtaMuonTrack", "#eta of reco muon track" , hDim.nBinEta, hDim.minEta, hDim.maxEta);
00621   commonME_->hMuonTrackPhi_ = theDQM->book1D("PhiMuonTrack", "#phi of reco muon track" , hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
00622   commonME_->hMuonTrackDxy_ = theDQM->book1D("DxyMuonTrack", "Dxy of reco muon track" , hDim.nBinDxy, hDim.minDxy, hDim.maxDxy);
00623   commonME_->hMuonTrackDz_ = theDQM->book1D("DzMuonTrack", "Dz of reco muon track" , hDim.nBinDz, hDim.minDz, hDim.maxDz);
00624 
00625   //histograms for fractions
00626   commonME_->hMuonAllP_   = theDQM->book1D("PMuonAll"  , "p of muons of all types"    , hDim.nBinP  , hDim.minP  , hDim.maxP  );
00627   commonME_->hMuonAllPt_  = theDQM->book1D("PtMuonAll" , "p_{T} of muon of all types", hDim.nBinPt , hDim.minPt , hDim.maxPt );
00628   commonME_->hMuonAllEta_ = theDQM->book1D("EtaMuonAll", "#eta of muon of all types" , hDim.nBinEta, hDim.minEta, hDim.maxEta);
00629   commonME_->hMuonAllPhi_ = theDQM->book1D("PhiMuonAll", "#phi of muon of all types" , hDim.nBinPhi, hDim.minPhi, hDim.maxPhi);
00630 
00631   muonME_->bookHistograms(theDQM, subDir_, hDim);
00632 
00633   if ( verbose_ > 0 ) theDQM->showDirStructure();
00634 }
00635 
00636 //
00637 //Destructor
00638 //
00639 RecoMuonValidator::~RecoMuonValidator()
00640 {
00641   if ( theMuonService ) delete theMuonService;
00642 }
00643 
00644 //
00645 //Begin run
00646 //
00647 void RecoMuonValidator::beginRun(const edm::Run& , const EventSetup& eventSetup)
00648 {
00649   if ( theMuonService ) theMuonService->update(eventSetup);
00650 
00651   if ( doAssoc_ ) {
00652     edm::ESHandle<TrackAssociatorBase> associatorBase;
00653     eventSetup.get<TrackAssociatorRecord>().get(muAssocLabel_.label(), associatorBase);
00654     assoByHits = dynamic_cast<const MuonAssociatorByHits *>(associatorBase.product());
00655     if (assoByHits == 0) throw cms::Exception("Configuration") << "The Track Associator with label '" << muAssocLabel_.label() << "' is not a MuonAssociatorByHits.\n";
00656     }
00657 
00658 }
00659 
00660 //
00661 //End run
00662 //
00663 void RecoMuonValidator::endRun()
00664 {
00665   if ( theDQM && ! outputFileName_.empty() ) theDQM->save(outputFileName_);
00666 }
00667 
00668 //
00669 //Analyze
00670 //
00671 void RecoMuonValidator::analyze(const Event& event, const EventSetup& eventSetup)
00672 {
00673   if ( ! theDQM ) {
00674     LogError("RecoMuonValidator") << "DQMService not initialized\n";
00675     return;
00676   }
00677 
00678   // Get TrackingParticles
00679   Handle<TrackingParticleCollection> simHandle;
00680   event.getByLabel(simLabel_, simHandle);
00681   const TrackingParticleCollection simColl = *(simHandle.product());
00682 
00683   // Get Muons
00684   Handle<View<Muon> > muonHandle;
00685   event.getByLabel(muonLabel_, muonHandle);
00686   View<Muon> muonColl = *(muonHandle.product());
00687 
00688   const TrackingParticleCollection::size_type nSim = simColl.size();
00689 
00690   edm::RefToBaseVector<reco::Muon> Muons;
00691   Muons = muonHandle->refVector();
00692 
00693   edm::RefVector<TrackingParticleCollection> allTPs;
00694   for (size_t i = 0; i < nSim; ++i) {
00695       allTPs.push_back(TrackingParticleRef(simHandle,i));
00696   }
00697 
00698   muonME_->hNSim_->Fill(nSim);
00699   muonME_->hNMuon_->Fill(muonColl.size());
00700 
00701   MuonAssociatorByHits::MuonToSimCollection muonToSimColl;
00702   MuonAssociatorByHits::SimToMuonCollection simToMuonColl;
00703 
00704   if ( doAssoc_ ) {
00705   assoByHits->associateMuons(muonToSimColl, simToMuonColl, Muons, trackType_, allTPs, &event, &eventSetup);
00706   } else {
00707 /*
00708     // SimToMuon associations
00709     Handle<reco::RecoToSimCollection> simToTrkMuHandle;
00710     event.getByLabel(trkMuAssocLabel_, simToTrkMuHandle);
00711     trkSimRecColl = *(simToTrkMuHandle.product());
00712 
00713     Handle<reco::RecoToSimCollection> simToStaMuHandle;
00714     event.getByLabel(staMuAssocLabel_, simToStaMuHandle);
00715     staSimRecColl = *(simToStaMuHandle.product());
00716 
00717     Handle<reco::RecoToSimCollection> simToGlbMuHandle;
00718     event.getByLabel(glbMuAssocLabel_, simToGlbMuHandle);
00719     glbSimRecColl = *(simToGlbMuHandle.product());
00720 
00721     // MuonToSim associations
00722     Handle<reco::SimToRecoCollection> trkMuToSimHandle;
00723     event.getByLabel(trkMuAssocLabel_, trkMuToSimHandle);
00724     trkRecSimColl = *(trkMuToSimHandle.product());
00725 
00726     Handle<reco::SimToRecoCollection> staMuToSimHandle;
00727     event.getByLabel(staMuAssocLabel_, staMuToSimHandle);
00728     staRecSimColl = *(staMuToSimHandle.product());
00729 
00730     Handle<reco::SimToRecoCollection> glbMuToSimHandle;
00731     event.getByLabel(glbMuAssocLabel_, glbMuToSimHandle);
00732     glbRecSimColl = *(glbMuToSimHandle.product());
00733 */
00734   }
00735 
00736   int glbNTrackerHits = 0; int trkNTrackerHits = 0;
00737   int glbNMuonHits = 0; int staNMuonHits = 0;
00738   int NTrackerHits = 0; int NMuonHits = 0;
00739   
00740   // Analyzer reco::Muon  
00741   for(View<Muon>::const_iterator iMuon = muonColl.begin();
00742       iMuon != muonColl.end(); ++iMuon) {
00743       
00744     double muonP, muonPt, muonEta, muonPhi;
00745     if (usePFMuon_) {
00746       muonP   = iMuon->pfP4().P();
00747       muonPt  = iMuon->pfP4().Pt();
00748       muonEta = iMuon->pfP4().Eta(); 
00749       muonPhi = iMuon->pfP4().Phi();
00750     }
00751     else {
00752       muonP   = iMuon->p(); 
00753       muonPt  = iMuon->pt();
00754       muonEta = iMuon->eta();
00755       muonPhi = iMuon->phi();
00756     }
00757 
00758     //histograms for fractions
00759     commonME_->hMuonAllP_->Fill(muonP);
00760     commonME_->hMuonAllPt_->Fill(muonPt);
00761     commonME_->hMuonAllEta_->Fill(muonEta);
00762     commonME_->hMuonAllPhi_->Fill(muonPhi);
00763 
00764 
00765     if (!selector_(*iMuon)) continue;
00766 
00767 
00768     TrackRef Track = iMuon->track();
00769 
00770     if (Track.isNonnull()) {
00771       commonME_->hMuonTrackP_->Fill(Track->p());
00772       commonME_->hMuonTrackPt_->Fill(Track->pt());
00773       commonME_->hMuonTrackEta_->Fill(Track->eta());
00774       commonME_->hMuonTrackPhi_->Fill(Track->phi());
00775       
00776       //ip histograms
00777       commonME_->hMuonTrackDxy_->Fill(Track->dxy());
00778       commonME_->hMuonTrackDz_->Fill(Track->dz());
00779     }
00780     
00781     if (iMuon->isGlobalMuon()) {
00782       Track = iMuon->combinedMuon();
00783       glbNTrackerHits = countTrackerHits(*Track);
00784       glbNMuonHits = countMuonHits(*Track);
00785     } else if (iMuon->isTrackerMuon()) {
00786       Track = iMuon->track();
00787       trkNTrackerHits = countTrackerHits(*Track);
00788     } else {
00789       Track = iMuon->standAloneMuon();
00790     } 
00791     
00792     NTrackerHits = countTrackerHits(*Track);
00793     muonME_->hNTrackerHits_->Fill(NTrackerHits);
00794     muonME_->hNTrackerHits_vs_Pt_->Fill(Track->pt(), NTrackerHits);
00795     muonME_->hNTrackerHits_vs_Eta_->Fill(Track->eta(), NTrackerHits);
00796 
00797     NMuonHits = countMuonHits(*Track);
00798     muonME_->hNMuonHits_->Fill(NMuonHits);
00799     muonME_->hNMuonHits_vs_Pt_->Fill(Track->pt(), NMuonHits);
00800     muonME_->hNMuonHits_vs_Eta_->Fill(Track->eta(), NMuonHits);
00801 
00802 //list of histos for each type
00803 
00804 //      muonME_->hNTrks_->Fill();
00805     muonME_->hNTrksEta_->Fill(Track->eta());
00806     muonME_->hNTrksPt_->Fill(Track->pt());
00807 
00808     commonME_->hMuonP_->Fill(muonP);
00809     commonME_->hMuonPt_->Fill(muonPt);
00810     commonME_->hMuonEta_->Fill(muonEta);
00811     commonME_->hMuonPhi_->Fill(muonPhi);
00812 
00813     if (iMuon->isGlobalMuon()) {
00814       double gtHitPat = iMuon->globalTrack()->hitPattern().numberOfHits() - iMuon->globalTrack()->hitPattern().numberOfValidHits();
00815       double itHitPat = iMuon->innerTrack()->hitPattern().numberOfHits() - iMuon->innerTrack()->hitPattern().numberOfValidHits();
00816       double otHitPat = iMuon->outerTrack()->hitPattern().numberOfHits() - iMuon->outerTrack()->hitPattern().numberOfValidHits();
00817       
00818       commonME_->hNInvalidHitsGTHitPattern_->Fill(gtHitPat);
00819       commonME_->hNInvalidHitsITHitPattern_->Fill(itHitPat);
00820       commonME_->hNInvalidHitsOTHitPattern_->Fill(otHitPat);
00821       commonME_->hNDeltaInvalidHitsHitPattern_->Fill(gtHitPat - itHitPat - otHitPat);
00822 
00823       //must be global and standalone
00824       if (iMuon->isStandAloneMuon()) { 
00825         commonME_->hStaToGlbDiffNMuonHitsEta_->Fill(Track->eta(),staNMuonHits-glbNMuonHits);
00826         commonME_->hStaToGlbDiffNMuonHitsPt_->Fill(Track->pt(),staNMuonHits-glbNMuonHits);
00827         commonME_->hStaToGlbDiffNMuonHits_->Fill(staNMuonHits-glbNMuonHits);
00828       }
00829       
00830       //must be global and tracker
00831       if (iMuon->isTrackerMuon()) {
00832         commonME_->hTrkToGlbDiffNTrackerHitsEta_->Fill(Track->eta(),trkNTrackerHits-glbNTrackerHits);
00833         commonME_->hTrkToGlbDiffNTrackerHitsPt_->Fill(Track->pt(),trkNTrackerHits-glbNTrackerHits);
00834         commonME_->hTrkToGlbDiffNTrackerHits_->Fill(trkNTrackerHits-glbNTrackerHits);
00835       }
00836     }
00837 
00838   }//end of reco muon loop
00839 
00840   // Associate by hits
00841   for(TrackingParticleCollection::size_type i=0; i<nSim; i++) {
00842     TrackingParticleRef simRef(simHandle, i);
00843     const TrackingParticle* simTP = simRef.get();
00844     if ( ! tpSelector_(*simTP) ) continue;
00845 
00846       //denominators for efficiency plots
00847       const double simP   = simRef->p();
00848       const double simPt  = simRef->pt();
00849       const double simEta = doAbsEta_ ? fabs(simRef->eta()) : simRef->eta();
00850       const double simPhi = simRef->phi();
00851 
00852       GlobalPoint  simVtx(simRef->vertex().x(), simRef->vertex().y(), simRef->vertex().z());
00853       GlobalVector simMom(simRef->momentum().x(), simRef->momentum().y(), simRef->momentum().z());
00854       const double simDxy = -simVtx.x()*sin(simPhi)+simVtx.y()*cos(simPhi);
00855       const double simDz  = simVtx.z() - (simVtx.x()*simMom.x()+simVtx.y()*simMom.y())*simMom.z()/simMom.perp2();
00856       
00857       const unsigned int nSimHits = simRef->pSimHit_end() - simRef->pSimHit_begin();
00858 
00859       muonME_->hSimP_  ->Fill(simP  );
00860       muonME_->hSimPt_ ->Fill(simPt );
00861       muonME_->hSimEta_->Fill(simEta);
00862       muonME_->hSimPhi_->Fill(simPhi);
00863       muonME_->hSimDxy_->Fill(simDxy);
00864       muonME_->hSimDz_->Fill(simDz);
00865       muonME_->hNSimHits_->Fill(nSimHits);
00866 
00867     // Get sim-reco association for a simRef
00868     vector<pair<RefToBase<Muon>, double> > MuRefV;
00869     if ( simToMuonColl.find(simRef) != simToMuonColl.end() ) {
00870       MuRefV = simToMuonColl[simRef];
00871 
00872       if ( !MuRefV.empty()) {
00873         muonME_->hNSimToReco_->Fill(MuRefV.size());
00874         const Muon* Mu = MuRefV.begin()->first.get();
00875         if (!selector_(*Mu)) continue;
00876 
00877         muonME_->fill(&*simTP, Mu); 
00878       }
00879     }
00880   }
00881 }
00882 
00883 int
00884 RecoMuonValidator::countMuonHits(const reco::Track& track) const {
00885   TransientTrackingRecHit::ConstRecHitContainer result;
00886   
00887   int count = 0;
00888 
00889   for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
00890     if((*hit)->isValid()) {
00891       DetId recoid = (*hit)->geographicalId();
00892       if ( recoid.det() == DetId::Muon ) count++;
00893     }
00894   }
00895   return count;
00896 }
00897 
00898 int
00899 RecoMuonValidator::countTrackerHits(const reco::Track& track) const {
00900   TransientTrackingRecHit::ConstRecHitContainer result;
00901   
00902   int count = 0;
00903 
00904   for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
00905     if((*hit)->isValid()) {
00906       DetId recoid = (*hit)->geographicalId();
00907       if ( recoid.det() == DetId::Tracker ) count++;
00908     }
00909   }
00910   return count;
00911 }
00912 /* vim:set ts=2 sts=2 sw=2 expandtab: */
/data/refman/pasoursint/CMSSW_5_2_9/src/Validation/RecoMuon/src/RecoMuonValidator.cc
