CMS 3D CMS Logo

/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/Validation/RecoMuon/src/RecoMuonValidator.cc

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