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