Reference Manual

Go to the documentation of this file.
00001 #include "Validation/RecoMuon/plugins/MuonTrackValidator.h"
00002 #include "DQMServices/ClientConfig/interface/FitSlicesYTool.h"
00003 
00004 #include "FWCore/Framework/interface/MakerMacros.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 
00007 #include "DataFormats/TrackReco/interface/Track.h"
00008 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00009 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00010 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
00011 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00012 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
00013 #include "SimTracker/TrackAssociation/interface/TrackAssociatorByChi2.h"
00014 #include "SimTracker/TrackAssociation/interface/TrackAssociatorByHits.h"
00015 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
00016 #include "SimTracker/Records/interface/TrackAssociatorRecord.h"
00017 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
00018 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00019 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
00020 #include "SimTracker/TrackAssociation/plugins/ParametersDefinerForTPESProducer.h"
00021 #include "SimTracker/TrackAssociation/plugins/CosmicParametersDefinerForTPESProducer.h"
00022 
00023 #include "TMath.h"
00024 #include <TF1.h>
00025 
00026 using namespace std;
00027 using namespace edm;
00028 
00029 void MuonTrackValidator::beginRun(Run const&, EventSetup const& setup) {
00030 
00031   //  dbe_->showDirStructure();
00032 
00033   int j=0;
00034   for (unsigned int ww=0;ww<associators.size();ww++){
00035     for (unsigned int www=0;www<label.size();www++){
00036 
00037       dbe_->cd();
00038       InputTag algo = label[www];
00039       string dirName=dirName_;
00040       if (algo.process()!="")
00041         dirName+=algo.process()+"_";
00042       if(algo.label()!="")
00043         dirName+=algo.label()+"_";
00044       if(algo.instance()!="")
00045         dirName+=algo.instance()+"_";      
00046       if (dirName.find("Tracks")<dirName.length()){
00047         dirName.replace(dirName.find("Tracks"),6,"");
00048       }
00049       string assoc= associators[ww];
00050       if (assoc.find("Track")<assoc.length()){
00051         assoc.replace(assoc.find("Track"),5,"");
00052       }
00053       dirName+=assoc;
00054       std::replace(dirName.begin(), dirName.end(), ':', '_');
00055       dbe_->setCurrentFolder(dirName.c_str());
00056 
00057       setUpVectors();
00058 
00059       dbe_->goUp();
00060       string subDirName = dirName + "/simulation";
00061       dbe_->setCurrentFolder(subDirName.c_str());
00062       h_ptSIM.push_back( dbe_->book1D("ptSIM", "generated p_{t}", 5500, 0, 110 ) );
00063       h_etaSIM.push_back( dbe_->book1D("etaSIM", "generated pseudorapidity", 500, -2.5, 2.5 ) );
00064       h_tracksSIM.push_back( dbe_->book1D("tracksSIM","number of simulated tracks",200,-0.5,99.5) );
00065       h_vertposSIM.push_back( dbe_->book1D("vertposSIM","Transverse position of sim vertices",100,0.,120.) );
00066       
00067       dbe_->cd();
00068       dbe_->setCurrentFolder(dirName.c_str());
00069       h_tracks.push_back( dbe_->book1D("tracks","number of reconstructed tracks",200,-0.5,19.5) );
00070       h_fakes.push_back( dbe_->book1D("fakes","number of fake reco tracks",20,-0.5,19.5) );
00071       h_charge.push_back( dbe_->book1D("charge","charge",3,-1.5,1.5) );
00072       h_hits.push_back( dbe_->book1D("hits", "number of hits per track", nintHit,minHit,maxHit ) );
00073       h_losthits.push_back( dbe_->book1D("losthits", "number of lost hits per track", nintHit,minHit,maxHit) );
00074       h_nchi2.push_back( dbe_->book1D("chi2", "normalized #chi^{2}", 200, 0, 20 ) );
00075       h_nchi2_prob.push_back( dbe_->book1D("chi2_prob", "normalized #chi^{2} probability",100,0,1));
00076 
00078       h_recoeta.push_back( dbe_->book1D("num_reco_eta","N of reco track vs eta",nint,min,max) );
00079       h_assoceta.push_back( dbe_->book1D("num_assoc(simToReco)_eta","N of associated tracks (simToReco) vs eta",nint,min,max) );
00080       h_assoc2eta.push_back( dbe_->book1D("num_assoc(recoToSim)_eta","N of associated (recoToSim) tracks vs eta",nint,min,max) );
00081       h_simuleta.push_back( dbe_->book1D("num_simul_eta","N of simulated tracks vs eta",nint,min,max) );
00082       h_recopT.push_back( dbe_->book1D("num_reco_pT","N of reco track vs pT",nintpT,minpT,maxpT) );
00083       h_assocpT.push_back( dbe_->book1D("num_assoc(simToReco)_pT","N of associated tracks (simToReco) vs pT",nintpT,minpT,maxpT) );
00084       h_assoc2pT.push_back( dbe_->book1D("num_assoc(recoToSim)_pT","N of associated (recoToSim) tracks vs pT",nintpT,minpT,maxpT) );
00085       h_simulpT.push_back( dbe_->book1D("num_simul_pT","N of simulated tracks vs pT",nintpT,minpT,maxpT) );
00086       //
00087       h_recohit.push_back( dbe_->book1D("num_reco_hit","N of reco track vs hit",nintHit,minHit,maxHit) );
00088       h_assochit.push_back( dbe_->book1D("num_assoc(simToReco)_hit","N of associated tracks (simToReco) vs hit",nintHit,minHit,maxHit) );
00089       h_assoc2hit.push_back( dbe_->book1D("num_assoc(recoToSim)_hit","N of associated (recoToSim) tracks vs hit",nintHit,minHit,maxHit) );
00090       h_simulhit.push_back( dbe_->book1D("num_simul_hit","N of simulated tracks vs hit",nintHit,minHit,maxHit) );
00091       //
00092       h_recophi.push_back( dbe_->book1D("num_reco_phi","N of reco track vs phi",nintPhi,minPhi,maxPhi) );
00093       h_assocphi.push_back( dbe_->book1D("num_assoc(simToReco)_phi","N of associated tracks (simToReco) vs phi",nintPhi,minPhi,maxPhi) );
00094       h_assoc2phi.push_back( dbe_->book1D("num_assoc(recoToSim)_phi","N of associated (recoToSim) tracks vs phi",nintPhi,minPhi,maxPhi) );
00095       h_simulphi.push_back( dbe_->book1D("num_simul_phi","N of simulated tracks vs phi",nintPhi,minPhi,maxPhi) );
00096 
00097       h_recodxy.push_back( dbe_->book1D("num_reco_dxy","N of reco track vs dxy",nintDxy,minDxy,maxDxy) );
00098       h_assocdxy.push_back( dbe_->book1D("num_assoc(simToReco)_dxy","N of associated tracks (simToReco) vs dxy",nintDxy,minDxy,maxDxy) );
00099       h_assoc2dxy.push_back( dbe_->book1D("num_assoc(recoToSim)_dxy","N of associated (recoToSim) tracks vs dxy",nintDxy,minDxy,maxDxy) );
00100       h_simuldxy.push_back( dbe_->book1D("num_simul_dxy","N of simulated tracks vs dxy",nintDxy,minDxy,maxDxy) );
00101       
00102       h_recodz.push_back( dbe_->book1D("num_reco_dz","N of reco track vs dz",nintDz,minDz,maxDz) );
00103       h_assocdz.push_back( dbe_->book1D("num_assoc(simToReco)_dz","N of associated tracks (simToReco) vs dz",nintDz,minDz,maxDz) );
00104       h_assoc2dz.push_back( dbe_->book1D("num_assoc(recoToSim)_dz","N of associated (recoToSim) tracks vs dz",nintDz,minDz,maxDz) );
00105       h_simuldz.push_back( dbe_->book1D("num_simul_dz","N of simulated tracks vs dz",nintDz,minDz,maxDz) );
00106 
00107       h_assocvertpos.push_back( dbe_->book1D("num_assoc(simToReco)_vertpos","N of associated tracks (simToReco) vs transverse vert position",nintVertpos,minVertpos,maxVertpos) );
00108       h_simulvertpos.push_back( dbe_->book1D("num_simul_vertpos","N of simulated tracks vs transverse vert position",nintVertpos,minVertpos,maxVertpos) );
00109 
00110       h_assoczpos.push_back( dbe_->book1D("num_assoc(simToReco)_zpos","N of associated tracks (simToReco) vs z vert position",nintZpos,minZpos,maxZpos) );
00111       h_simulzpos.push_back( dbe_->book1D("num_simul_zpos","N of simulated tracks vs z vert position",nintZpos,minZpos,maxZpos) );
00112 
00113 
00115 
00116       h_eta.push_back( dbe_->book1D("eta", "pseudorapidity residue", 1000, -0.1, 0.1 ) );
00117       h_pt.push_back( dbe_->book1D("pullPt", "pull of p_{t}", 100, -10, 10 ) );
00118       h_pullTheta.push_back( dbe_->book1D("pullTheta","pull of #theta parameter",250,-25,25) );
00119       h_pullPhi.push_back( dbe_->book1D("pullPhi","pull of #phi parameter",250,-25,25) );
00120       h_pullDxy.push_back( dbe_->book1D("pullDxy","pull of dxy parameter",250,-25,25) );
00121       h_pullDz.push_back( dbe_->book1D("pullDz","pull of dz parameter",250,-25,25) );
00122       h_pullQoverp.push_back( dbe_->book1D("pullQoverp","pull of qoverp parameter",250,-25,25) );
00123       
00124       if (associators[ww]=="TrackAssociatorByChi2"){
00125         h_assochi2.push_back( dbe_->book1D("assocChi2","track association #chi^{2}",1000000,0,100000) );
00126         h_assochi2_prob.push_back(dbe_->book1D("assocChi2_prob","probability of association #chi^{2}",100,0,1));
00127       } else if (associators[ww]=="TrackAssociatorByHits"){
00128         h_assocFraction.push_back( dbe_->book1D("assocFraction","fraction of shared hits",200,0,2) );
00129         h_assocSharedHit.push_back(dbe_->book1D("assocSharedHit","number of shared hits",20,0,20));
00130       }
00131 
00132       chi2_vs_nhits.push_back( dbe_->book2D("chi2_vs_nhits","#chi^{2} vs nhits",25,0,25,100,0,10) );
00133       h_chi2meanhitsh.push_back( dbe_->bookProfile("chi2mean_vs_nhits","mean #chi^{2} vs nhits",25,0,25,100,0,10) );
00134 
00135       etares_vs_eta.push_back( dbe_->book2D("etares_vs_eta","etaresidue vs eta",nint,min,max,200,-0.1,0.1) );
00136       nrec_vs_nsim.push_back( dbe_->book2D("nrec_vs_nsim","nrec vs nsim",20,-0.5,19.5,20,-0.5,19.5) );
00137 
00138       chi2_vs_eta.push_back( dbe_->book2D("chi2_vs_eta","chi2_vs_eta",nint,min,max, 200, 0, 20 ));
00139       h_chi2meanh.push_back( dbe_->bookProfile("chi2mean","mean #chi^{2} vs #eta",nint,min,max, 200, 0, 20) );
00140       chi2_vs_phi.push_back( dbe_->book2D("chi2_vs_phi","#chi^{2} vs #phi",nintPhi,minPhi,maxPhi, 200, 0, 20 ) );
00141       h_chi2mean_vs_phi.push_back( dbe_->bookProfile("chi2mean_vs_phi","mean of #chi^{2} vs #phi",nintPhi,minPhi,maxPhi, 200, 0, 20) );
00142 
00143       nhits_vs_eta.push_back( dbe_->book2D("nhits_vs_eta","nhits vs eta",nint,min,max,nintHit,minHit,maxHit) );
00144       nDThits_vs_eta.push_back( dbe_->book2D("nDThits_vs_eta","# DT hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
00145       nCSChits_vs_eta.push_back( dbe_->book2D("nCSChits_vs_eta","# CSC hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
00146       nRPChits_vs_eta.push_back( dbe_->book2D("nRPChits_vs_eta","# RPC hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
00147 
00148       h_DThits_eta.push_back( dbe_->bookProfile("DThits_eta","mean # DT hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
00149       h_CSChits_eta.push_back( dbe_->bookProfile("CSChits_eta","mean # CSC hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
00150       h_RPChits_eta.push_back( dbe_->bookProfile("RPChits_eta","mean # RPC hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
00151       h_hits_eta.push_back( dbe_->bookProfile("hits_eta","mean #hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
00152       nhits_vs_phi.push_back( dbe_->book2D("nhits_vs_phi","#hits vs #phi",nintPhi,minPhi,maxPhi,nintHit,minHit,maxHit) );
00153       h_hits_phi.push_back( dbe_->bookProfile("hits_phi","mean #hits vs #phi",nintPhi,minPhi,maxPhi, nintHit,minHit,maxHit) );
00154 
00155       nlosthits_vs_eta.push_back( dbe_->book2D("nlosthits_vs_eta","nlosthits vs eta",nint,min,max,nintHit,minHit,maxHit) );
00156       h_losthits_eta.push_back( dbe_->bookProfile("losthits_eta","losthits_eta",nint,min,max,nintHit,minHit,maxHit) );
00157 
00158       ptres_vs_eta.push_back(dbe_->book2D("ptres_vs_eta","ptres_vs_eta",nint,min,max, ptRes_nbin, ptRes_rangeMin, ptRes_rangeMax));
00159       ptres_vs_phi.push_back( dbe_->book2D("ptres_vs_phi","p_{t} res vs #phi",nintPhi,minPhi,maxPhi, ptRes_nbin, ptRes_rangeMin, ptRes_rangeMax));
00160       ptres_vs_pt.push_back(dbe_->book2D("ptres_vs_pt","ptres_vs_pt",nintpT,minpT,maxpT, ptRes_nbin, ptRes_rangeMin, ptRes_rangeMax));
00161 
00162       cotThetares_vs_eta.push_back(dbe_->book2D("cotThetares_vs_eta","cotThetares_vs_eta",nint,min,max,cotThetaRes_nbin, cotThetaRes_rangeMin, cotThetaRes_rangeMax));
00163       cotThetares_vs_pt.push_back(dbe_->book2D("cotThetares_vs_pt","cotThetares_vs_pt",nintpT,minpT,maxpT, cotThetaRes_nbin, cotThetaRes_rangeMin, cotThetaRes_rangeMax));
00164 
00165       phires_vs_eta.push_back(dbe_->book2D("phires_vs_eta","phires_vs_eta",nint,min,max, phiRes_nbin, phiRes_rangeMin, phiRes_rangeMax));
00166       phires_vs_pt.push_back(dbe_->book2D("phires_vs_pt","phires_vs_pt",nintpT,minpT,maxpT, phiRes_nbin, phiRes_rangeMin, phiRes_rangeMax));
00167       phires_vs_phi.push_back(dbe_->book2D("phires_vs_phi","#phi res vs #phi",nintPhi,minPhi,maxPhi,phiRes_nbin, phiRes_rangeMin, phiRes_rangeMax));
00168 
00169       dxyres_vs_eta.push_back(dbe_->book2D("dxyres_vs_eta","dxyres_vs_eta",nint,min,max,dxyRes_nbin, dxyRes_rangeMin, dxyRes_rangeMax));
00170       dxyres_vs_pt.push_back( dbe_->book2D("dxyres_vs_pt","dxyres_vs_pt",nintpT,minpT,maxpT,dxyRes_nbin, dxyRes_rangeMin, dxyRes_rangeMax));
00171 
00172       dzres_vs_eta.push_back(dbe_->book2D("dzres_vs_eta","dzres_vs_eta",nint,min,max,dzRes_nbin, dzRes_rangeMin, dzRes_rangeMax));
00173       dzres_vs_pt.push_back(dbe_->book2D("dzres_vs_pt","dzres_vs_pt",nintpT,minpT,maxpT,dzRes_nbin, dzRes_rangeMin, dzRes_rangeMax));
00174 
00175       ptmean_vs_eta_phi.push_back(dbe_->bookProfile2D("ptmean_vs_eta_phi","mean p_{t} vs #eta and #phi",nintPhi,minPhi,maxPhi,nint,min,max,1000,0,1000));
00176       phimean_vs_eta_phi.push_back(dbe_->bookProfile2D("phimean_vs_eta_phi","mean #phi vs #eta and #phi",nintPhi,minPhi,maxPhi,nint,min,max,nintPhi,minPhi,maxPhi));
00177 
00178       //pulls of track params vs eta: to be used with fitslicesytool
00179       dxypull_vs_eta.push_back(dbe_->book2D("dxypull_vs_eta","dxypull_vs_eta",nint,min,max,100,-10,10));
00180       ptpull_vs_eta.push_back(dbe_->book2D("ptpull_vs_eta","ptpull_vs_eta",nint,min,max,100,-10,10)); 
00181       dzpull_vs_eta.push_back(dbe_->book2D("dzpull_vs_eta","dzpull_vs_eta",nint,min,max,100,-10,10)); 
00182       phipull_vs_eta.push_back(dbe_->book2D("phipull_vs_eta","phipull_vs_eta",nint,min,max,100,-10,10)); 
00183       thetapull_vs_eta.push_back(dbe_->book2D("thetapull_vs_eta","thetapull_vs_eta",nint,min,max,100,-10,10));
00184 
00185       //pulls of track params vs phi
00186       ptpull_vs_phi.push_back(dbe_->book2D("ptpull_vs_phi","p_{t} pull vs #phi",nintPhi,minPhi,maxPhi,100,-10,10)); 
00187       phipull_vs_phi.push_back(dbe_->book2D("phipull_vs_phi","#phi pull vs #phi",nintPhi,minPhi,maxPhi,100,-10,10)); 
00188       thetapull_vs_phi.push_back(dbe_->book2D("thetapull_vs_phi","#theta pull vs #phi",nintPhi,minPhi,maxPhi,100,-10,10));
00189 
00190       nrecHit_vs_nsimHit_sim2rec.push_back( dbe_->book2D("nrecHit_vs_nsimHit_sim2rec","nrecHit vs nsimHit (Sim2RecAssoc)",nintHit,minHit,maxHit, nintHit,minHit,maxHit ));
00191       nrecHit_vs_nsimHit_rec2sim.push_back( dbe_->book2D("nrecHit_vs_nsimHit_rec2sim","nrecHit vs nsimHit (Rec2simAssoc)",nintHit,minHit,maxHit, nintHit,minHit,maxHit ));
00192       
00193       if (MABH) {
00194         h_PurityVsQuality.push_back( dbe_->book2D("PurityVsQuality","Purity vs Quality (MABH)",20,0.01,1.01,20,0.01,1.01) );
00195         h_assoceta_Quality05.push_back( dbe_->book1D("num_assoc(simToReco)_eta_Q05","N of associated tracks (simToReco) vs eta (Quality>0.5)",nint,min,max) );
00196         h_assoceta_Quality075.push_back( dbe_->book1D("num_assoc(simToReco)_eta_Q075","N of associated tracks (simToReco) vs eta (Quality>0.75)",nint,min,max) );
00197         h_assocpT_Quality05.push_back( dbe_->book1D("num_assoc(simToReco)_pT_Q05","N of associated tracks (simToReco) vs pT (Quality>0.5)",nintpT,minpT,maxpT) );
00198         h_assocpT_Quality075.push_back( dbe_->book1D("num_assoc(simToReco)_pT_Q075","N of associated tracks (simToReco) vs pT (Quality>0.75)",nintpT,minpT,maxpT) );
00199         h_assocphi_Quality05.push_back( dbe_->book1D("num_assoc(simToReco)_phi_Q05","N of associated tracks (simToReco) vs phi (Quality>0.5)",nintPhi,minPhi,maxPhi) );
00200         h_assocphi_Quality075.push_back( dbe_->book1D("num_assoc(simToReco)_phi_Q075","N of associated tracks (simToReco) vs phi (Quality>0.75)",nintPhi,minPhi,maxPhi) );
00201       }
00202 
00203       if(useLogPt){
00204         BinLogX(dzres_vs_pt[j]->getTH2F());
00205         BinLogX(dxyres_vs_pt[j]->getTH2F());
00206         BinLogX(phires_vs_pt[j]->getTH2F());
00207         BinLogX(cotThetares_vs_pt[j]->getTH2F());
00208         BinLogX(ptres_vs_pt[j]->getTH2F());
00209         BinLogX(h_recopT[j]->getTH1F());
00210         BinLogX(h_assocpT[j]->getTH1F());
00211         BinLogX(h_assoc2pT[j]->getTH1F());
00212         BinLogX(h_simulpT[j]->getTH1F());
00213         if (MABH)       {
00214           BinLogX(h_assocpT_Quality05[j]->getTH1F());
00215           BinLogX(h_assocpT_Quality075[j]->getTH1F());
00216         }
00217         j++;
00218       }
00219 
00220     }
00221   }
00222   if (UseAssociators) {
00223     edm::ESHandle<TrackAssociatorBase> theAssociator;
00224     for (unsigned int w=0;w<associators.size();w++) {
00225       setup.get<TrackAssociatorRecord>().get(associators[w],theAssociator);
00226       associator.push_back( theAssociator.product() );
00227     }
00228   }
00229 }
00230 
00231 void MuonTrackValidator::analyze(const edm::Event& event, const edm::EventSetup& setup){
00232   using namespace reco;
00233   
00234   edm::LogInfo("MuonTrackValidator") << "\n====================================================" << "\n"
00235                                      << "Analyzing new event" << "\n"
00236                                      << "====================================================\n" << "\n";
00237   edm::ESHandle<ParametersDefinerForTP> parametersDefinerTP; 
00238   setup.get<TrackAssociatorRecord>().get(parametersDefiner,parametersDefinerTP);    
00239   
00240   edm::Handle<TrackingParticleCollection>  TPCollectionHeff ;
00241   event.getByLabel(label_tp_effic,TPCollectionHeff);
00242   const TrackingParticleCollection tPCeff = *(TPCollectionHeff.product());
00243   
00244   edm::Handle<TrackingParticleCollection>  TPCollectionHfake ;
00245   event.getByLabel(label_tp_fake,TPCollectionHfake);
00246   const TrackingParticleCollection tPCfake = *(TPCollectionHfake.product());
00247   
00248   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
00249   event.getByLabel(bsSrc,recoBeamSpotHandle);
00250   reco::BeamSpot bs = *recoBeamSpotHandle;      
00251   
00252   int w=0;
00253   for (unsigned int ww=0;ww<associators.size();ww++){
00254     for (unsigned int www=0;www<label.size();www++){
00255       //
00256       //get collections from the event
00257       //
00258       edm::Handle<View<Track> >  trackCollection;
00259 
00260       reco::RecoToSimCollection recSimColl;
00261       reco::SimToRecoCollection simRecColl;
00262       unsigned int trackCollectionSize = 0;
00263 
00264       //      if(!event.getByLabel(label[www], trackCollection)&&ignoremissingtkcollection_) continue;
00265       if(!event.getByLabel(label[www], trackCollection)&&ignoremissingtkcollection_) {
00266 
00267         recSimColl.post_insert();
00268         simRecColl.post_insert();
00269 
00270       }
00271 
00272       else {
00273 
00274         trackCollectionSize = trackCollection->size();
00275         //associate tracks
00276         if(UseAssociators){
00277           edm::LogVerbatim("MuonTrackValidator") << "Analyzing " 
00278                                                  << label[www].process()<<":"
00279                                                  << label[www].label()<<":"
00280                                                  << label[www].instance()<<" with "
00281                                                  << associators[ww].c_str() <<"\n";
00282         
00283           LogTrace("MuonTrackValidator") << "Calling associateRecoToSim method" << "\n";
00284           recSimColl=associator[ww]->associateRecoToSim(trackCollection,
00285                                                         TPCollectionHfake,
00286                                                         &event);
00287           LogTrace("MuonTrackValidator") << "Calling associateSimToReco method" << "\n";
00288           simRecColl=associator[ww]->associateSimToReco(trackCollection,
00289                                                         TPCollectionHeff, 
00290                                                         &event);
00291         }
00292         else{
00293           edm::LogVerbatim("MuonTrackValidator") << "Analyzing " 
00294                                                  << label[www].process()<<":"
00295                                                  << label[www].label()<<":"
00296                                                  << label[www].instance()<<" with "
00297                                                  << associatormap.process()<<":"
00298                                                  << associatormap.label()<<":"
00299                                                  << associatormap.instance()<<"\n";
00300         
00301           Handle<reco::SimToRecoCollection > simtorecoCollectionH;
00302           event.getByLabel(associatormap,simtorecoCollectionH);
00303           simRecColl= *(simtorecoCollectionH.product()); 
00304         
00305           Handle<reco::RecoToSimCollection > recotosimCollectionH;
00306           event.getByLabel(associatormap,recotosimCollectionH);
00307           recSimColl= *(recotosimCollectionH.product()); 
00308         }
00309 
00310       }
00311 
00312       
00313       //
00314       //fill simulation histograms
00315       //compute number of tracks per eta interval
00316       //
00317       edm::LogVerbatim("MuonTrackValidator") << "\n# of TrackingParticles: " << tPCeff.size() << "\n";
00318       int ats = 0;
00319       int st=0;
00320       for (TrackingParticleCollection::size_type i=0; i<tPCeff.size(); i++){
00321         bool TP_is_matched = false;
00322         double quality = 0.;
00323         bool Quality05  = false;
00324         bool Quality075 = false;
00325 
00326         TrackingParticleRef tpr(TPCollectionHeff, i);
00327         TrackingParticle* tp=const_cast<TrackingParticle*>(tpr.get());
00328         ParticleBase::Vector momentumTP; 
00329         ParticleBase::Point vertexTP;
00330         double dxySim = 0;
00331         double dzSim = 0; 
00332 
00333         //If the TrackingParticle is collison like, get the momentum and vertex at production state
00334         if(parametersDefiner=="LhcParametersDefinerForTP")
00335           {
00336             if(! tpSelector(*tp)) continue;
00337             momentumTP = tp->momentum();
00338             vertexTP = tp->vertex();
00339             //Calcualte the impact parameters w.r.t. PCA
00340             ParticleBase::Vector momentum = parametersDefinerTP->momentum(event,setup,*tp);
00341             ParticleBase::Point vertex = parametersDefinerTP->vertex(event,setup,*tp);
00342             dxySim = (-vertex.x()*sin(momentum.phi())+vertex.y()*cos(momentum.phi()));
00343             dzSim = vertex.z() - (vertex.x()*momentum.x()+vertex.y()*momentum.y())/sqrt(momentum.perp2()) * momentum.z()/sqrt(momentum.perp2());
00344           }
00345         //If the TrackingParticle is comics, get the momentum and vertex at PCA
00346         if(parametersDefiner=="CosmicParametersDefinerForTP")
00347           {
00348             if(! cosmictpSelector(*tp,&bs,event,setup)) continue;       
00349             momentumTP = parametersDefinerTP->momentum(event,setup,*tp);
00350             vertexTP = parametersDefinerTP->vertex(event,setup,*tp);
00351             dxySim = (-vertexTP.x()*sin(momentumTP.phi())+vertexTP.y()*cos(momentumTP.phi()));
00352             dzSim = vertexTP.z() - (vertexTP.x()*momentumTP.x()+vertexTP.y()*momentumTP.y())/sqrt(momentumTP.perp2()) * momentumTP.z()/sqrt(momentumTP.perp2());
00353           }
00354         edm::LogVerbatim("MuonTrackValidator") <<"--------------------Selected TrackingParticle #"<<tpr.key();
00355         st++;
00356 
00357         h_ptSIM[w]->Fill(sqrt(momentumTP.perp2()));
00358         h_etaSIM[w]->Fill(momentumTP.eta());
00359         h_vertposSIM[w]->Fill(sqrt(vertexTP.perp2()));
00360         
00361         std::vector<std::pair<RefToBase<Track>, double> > rt;
00362         if(simRecColl.find(tpr) != simRecColl.end()){
00363           rt = (std::vector<std::pair<RefToBase<Track>, double> >) simRecColl[tpr];
00364           if (rt.size()!=0) {
00365             RefToBase<Track> assoc_recoTrack = rt.begin()->first;
00366             edm::LogVerbatim("MuonTrackValidator")<<"-----------------------------associated Track #"<<assoc_recoTrack.key();
00367             TP_is_matched = true;
00368             ats++;
00369             quality = rt.begin()->second;
00370             edm::LogVerbatim("MuonTrackValidator") << "TrackingParticle #" <<tpr.key()  
00371                                                    << " with pt=" << sqrt(momentumTP.perp2()) 
00372                                                    << " associated with quality:" << quality <<"\n";
00373             if (MABH) {
00374               if (quality > 0.75) {
00375                 Quality075 = true;
00376                 Quality05  = true;
00377               } 
00378               else if (quality > 0.5) {
00379                 Quality05  = true;
00380               }
00381             }       
00382           }
00383         }else{
00384           edm::LogVerbatim("MuonTrackValidator") 
00385             << "TrackingParticle #" << tpr.key()
00386             << " with pt,eta,phi: " 
00387             << sqrt(momentumTP.perp2()) << " , "
00388             << momentumTP.eta() << " , "
00389             << momentumTP.phi() << " , "
00390             << " NOT associated to any reco::Track" << "\n";
00391         }
00392         
00393         for (unsigned int f=0; f<etaintervals[w].size()-1; f++){
00394           if (getEta(momentumTP.eta())>etaintervals[w][f]&&
00395               getEta(momentumTP.eta())<etaintervals[w][f+1]) {
00396             totSIMeta[w][f]++;
00397             if (TP_is_matched) {
00398               totASSeta[w][f]++;
00399 
00400               if (MABH) {
00401                 if (Quality075) {
00402                   totASSeta_Quality075[w][f]++;
00403                   totASSeta_Quality05[w][f]++;
00404                 }
00405                 else if (Quality05) {
00406                   totASSeta_Quality05[w][f]++;
00407                 }
00408               }
00409             }
00410           }
00411         } // END for (unsigned int f=0; f<etaintervals[w].size()-1; f++){
00412         
00413         for (unsigned int f=0; f<phiintervals[w].size()-1; f++){
00414           if (momentumTP.phi() > phiintervals[w][f]&&
00415               momentumTP.phi() <phiintervals[w][f+1]) {
00416             totSIM_phi[w][f]++;
00417             if (TP_is_matched) {
00418               totASS_phi[w][f]++;
00419               
00420               if (MABH) {
00421                 if (Quality075) {
00422                   totASS_phi_Quality075[w][f]++;
00423                   totASS_phi_Quality05[w][f]++;
00424                 }
00425                 else if (Quality05) {
00426                   totASS_phi_Quality05[w][f]++;
00427                 }
00428               }
00429             }
00430           }
00431         } // END for (unsigned int f=0; f<phiintervals[w].size()-1; f++){
00432         
00433         
00434         for (unsigned int f=0; f<pTintervals[w].size()-1; f++){
00435           if (getPt(sqrt(momentumTP.perp2()))>pTintervals[w][f]&&
00436               getPt(sqrt(momentumTP.perp2()))<pTintervals[w][f+1]) {
00437             totSIMpT[w][f]++;
00438             if (TP_is_matched) {
00439               totASSpT[w][f]++;
00440               
00441               if (MABH) {
00442                 if (Quality075) {
00443                   totASSpT_Quality075[w][f]++;
00444                   totASSpT_Quality05[w][f]++;
00445                 }
00446                 else if (Quality05) { 
00447                   totASSpT_Quality05[w][f]++;
00448                 }
00449               }
00450             }
00451           }
00452         } // END for (unsigned int f=0; f<pTintervals[w].size()-1; f++){
00453         
00454         for (unsigned int f=0; f<dxyintervals[w].size()-1; f++){
00455           if (dxySim>dxyintervals[w][f]&&
00456               dxySim<dxyintervals[w][f+1]) {
00457             totSIM_dxy[w][f]++;
00458             if (TP_is_matched) {
00459               totASS_dxy[w][f]++;
00460             }
00461           }
00462         } // END for (unsigned int f=0; f<dxyintervals[w].size()-1; f++){
00463 
00464         for (unsigned int f=0; f<dzintervals[w].size()-1; f++){
00465           if (dzSim>dzintervals[w][f]&&
00466               dzSim<dzintervals[w][f+1]) {
00467             totSIM_dz[w][f]++;
00468             if (TP_is_matched) {
00469               totASS_dz[w][f]++;
00470             }
00471           }
00472         } // END for (unsigned int f=0; f<dzintervals[w].size()-1; f++){
00473 
00474         for (unsigned int f=0; f<vertposintervals[w].size()-1; f++){
00475           if (sqrt(vertexTP.perp2())>vertposintervals[w][f]&&
00476               sqrt(vertexTP.perp2())<vertposintervals[w][f+1]) {
00477             totSIM_vertpos[w][f]++;
00478             if (TP_is_matched) {
00479               totASS_vertpos[w][f]++;
00480             }
00481           }
00482         } // END for (unsigned int f=0; f<vertposintervals[w].size()-1; f++){
00483 
00484         for (unsigned int f=0; f<zposintervals[w].size()-1; f++){
00485           if (vertexTP.z()>zposintervals[w][f]&&
00486               vertexTP.z()<zposintervals[w][f+1]) {
00487             totSIM_zpos[w][f]++;
00488             if (TP_is_matched) {
00489               totASS_zpos[w][f]++;
00490             }
00491           }
00492         } // END for (unsigned int f=0; f<zposintervals[w].size()-1; f++){
00493         
00494         std::vector<PSimHit> simhits;
00495         
00496         if (usetracker && usemuon) {
00497           simhits=tp->trackPSimHit();
00498         } 
00499         else if (!usetracker && usemuon) {
00500           simhits=tp->trackPSimHit(DetId::Muon);
00501         }
00502         else if (usetracker && !usemuon) {
00503           simhits=tp->trackPSimHit(DetId::Tracker);
00504         }
00505         
00506         int tmp = std::min((int)(simhits.end()-simhits.begin()),int(maxHit-1));
00507         edm::LogVerbatim("MuonTrackValidator") << "\t N simhits = "<< (int)(simhits.end()-simhits.begin())<<"\n";
00508 
00509         totSIM_hit[w][tmp]++;
00510         if (TP_is_matched) totASS_hit[w][tmp]++;
00511 
00512         if (TP_is_matched)       
00513           {
00514             RefToBase<Track> assoctrack = rt.begin()->first; 
00515             nrecHit_vs_nsimHit_sim2rec[w]->Fill( assoctrack->numberOfValidHits(),(int)(simhits.end()-simhits.begin() ));
00516           }
00517       } // End  for (TrackingParticleCollection::size_type i=0; i<tPCeff.size(); i++){
00518       if (st!=0) h_tracksSIM[w]->Fill(st);
00519       
00520 
00521       //
00522       //fill reconstructed track histograms
00523       // 
00524       edm::LogVerbatim("MuonTrackValidator") << "\n# of reco::Tracks with "
00525                                          << label[www].process()<<":"
00526                                          << label[www].label()<<":"
00527                                          << label[www].instance()
00528                                          << ": " << trackCollectionSize << "\n";
00529       int at=0;
00530       int rT=0;
00531       for(View<Track>::size_type i=0; i<trackCollectionSize; ++i){
00532         bool Track_is_matched = false; 
00533         RefToBase<Track> track(trackCollection, i);
00534         rT++;
00535 
00536         std::vector<std::pair<TrackingParticleRef, double> > tp;
00537         TrackingParticleRef tpr;
00538 
00539         // new logic (bidirectional)
00540         if (BiDirectional_RecoToSim_association) {        
00541           edm::LogVerbatim("MuonTrackValidator")<<"----------------------------------------Track #"<< track.key();
00542 
00543           if(recSimColl.find(track) != recSimColl.end()) {
00544             tp = recSimColl[track];         
00545             if (tp.size() != 0) {
00546               tpr = tp.begin()->first;        
00547               // RtS and StR must associate the same pair !
00548               if(simRecColl.find(tpr) != simRecColl.end()) {
00549                 std::vector<std::pair<RefToBase<Track>, double> > track_checkback  = simRecColl[tpr];
00550                 RefToBase<Track> assoc_track_checkback;
00551                 assoc_track_checkback = track_checkback.begin()->first;
00552 
00553                 if ( assoc_track_checkback.key() == track.key() ) {
00554                   edm::LogVerbatim("MuonTrackValidator")<<"------------------associated TrackingParticle #"<<tpr.key();
00555                   Track_is_matched = true;
00556                   at++;
00557                   double Purity = tp.begin()->second;
00558                   double Quality = track_checkback.begin()->second;
00559                   edm::LogVerbatim("MuonTrackValidator") << "reco::Track #" << track.key() << " with pt=" << track->pt() 
00560                                                          << " associated with quality:" << Purity <<"\n";
00561                   if (MABH) h_PurityVsQuality[w]->Fill(Quality,Purity);
00562                 }
00563               }
00564             }
00565           }
00566 
00567           if (!Track_is_matched) edm::LogVerbatim("MuonTrackValidator") 
00568             << "reco::Track #" << track.key() << " with pt=" << track->pt() << " NOT associated to any TrackingParticle" << "\n";
00569         }
00570         // old logic (bugged)
00571         else {
00572           if(recSimColl.find(track) != recSimColl.end()){
00573             tp = recSimColl[track];
00574             if (tp.size()!=0) {
00575               Track_is_matched = true;
00576               tpr = tp.begin()->first;
00577               at++;
00578               edm::LogVerbatim("MuonTrackValidator") << "reco::Track #" << track.key() << " with pt=" << track->pt() 
00579                                                      << " associated with quality:" << tp.begin()->second <<"\n";
00580             }
00581           } else {
00582             edm::LogVerbatim("MuonTrackValidator") << "reco::Track #" << track.key() << " with pt=" << track->pt()
00583                                                    << " NOT associated to any TrackingParticle" << "\n";                  
00584           }
00585         }
00586         
00587         //Compute fake rate vs eta
00588         for (unsigned int f=0; f<etaintervals[w].size()-1; f++){
00589           if (getEta(track->momentum().eta())>etaintervals[w][f]&&
00590               getEta(track->momentum().eta())<etaintervals[w][f+1]) {
00591             totRECeta[w][f]++; 
00592             if (Track_is_matched) {
00593               totASS2eta[w][f]++;
00594             }           
00595           }
00596         } // End for (unsigned int f=0; f<etaintervals[w].size()-1; f++){
00597 
00598         for (unsigned int f=0; f<phiintervals[w].size()-1; f++){
00599           if (track->momentum().phi()>phiintervals[w][f]&&
00600               track->momentum().phi()<phiintervals[w][f+1]) {
00601             totREC_phi[w][f]++; 
00602             if (Track_is_matched) {
00603               totASS2_phi[w][f]++;
00604             }           
00605           }
00606         } // End for (unsigned int f=0; f<phiintervals[w].size()-1; f++){
00607 
00608         
00609         for (unsigned int f=0; f<pTintervals[w].size()-1; f++){
00610           if (getPt(sqrt(track->momentum().perp2()))>pTintervals[w][f]&&
00611               getPt(sqrt(track->momentum().perp2()))<pTintervals[w][f+1]) {
00612             totRECpT[w][f]++; 
00613             if (Track_is_matched) {
00614               totASS2pT[w][f]++;
00615             }         
00616           }
00617         } // End for (unsigned int f=0; f<pTintervals[w].size()-1; f++){
00618 
00619         for (unsigned int f=0; f<dxyintervals[w].size()-1; f++){
00620           if (track->dxy(bs.position())>dxyintervals[w][f]&&
00621               track->dxy(bs.position())<dxyintervals[w][f+1]) {
00622             totREC_dxy[w][f]++; 
00623             if (Track_is_matched) {
00624               totASS2_dxy[w][f]++;
00625             }         
00626           }
00627         } // End for (unsigned int f=0; f<dxyintervals[w].size()-1; f++){
00628 
00629         for (unsigned int f=0; f<dzintervals[w].size()-1; f++){
00630           if (track->dz(bs.position())>dzintervals[w][f]&&
00631               track->dz(bs.position())<dzintervals[w][f+1]) {
00632             totREC_dz[w][f]++; 
00633             if (Track_is_matched) {
00634               totASS2_dz[w][f]++;
00635             }         
00636           }
00637         } // End for (unsigned int f=0; f<dzintervals[w].size()-1; f++){
00638 
00639         int tmp = std::min((int)track->found(),int(maxHit-1));
00640         totREC_hit[w][tmp]++;
00641         if (Track_is_matched) totASS2_hit[w][tmp]++;
00642 
00643         edm::LogVerbatim("MuonTrackValidator") << "\t N valid rechits = "<< (int)track->found() <<"\n";
00644 
00645         //Fill other histos
00646         try{
00647           if (!Track_is_matched) continue;
00648 
00649           if (associators[ww]=="TrackAssociatorByChi2"){
00650             //association chi2
00651             double assocChi2 = -tp.begin()->second;//in association map is stored -chi2
00652             h_assochi2[www]->Fill(assocChi2);
00653             h_assochi2_prob[www]->Fill(TMath::Prob((assocChi2)*5,5));
00654           }
00655           else if (associators[ww]=="TrackAssociatorByHits"){
00656             double fraction = tp.begin()->second;
00657             h_assocFraction[www]->Fill(fraction);
00658             h_assocSharedHit[www]->Fill(fraction*track->numberOfValidHits());
00659           }
00660     
00661           //nchi2 and hits global distributions
00662           h_nchi2[w]->Fill(track->normalizedChi2());
00663           h_nchi2_prob[w]->Fill(TMath::Prob(track->chi2(),(int)track->ndof()));
00664           h_hits[w]->Fill(track->numberOfValidHits());
00665           h_losthits[w]->Fill(track->numberOfLostHits());
00666           chi2_vs_nhits[w]->Fill(track->numberOfValidHits(),track->normalizedChi2());
00667           h_charge[w]->Fill( track->charge() );
00668           
00669           //Get tracking particle parameters at point of closest approach to the beamline
00670           ParticleBase::Vector momentumTP = parametersDefinerTP->momentum(event,setup,*(tpr.get()));
00671           ParticleBase::Point vertexTP = parametersDefinerTP->vertex(event,setup,*(tpr.get()));
00672           double ptSim = sqrt(momentumTP.perp2());
00673           double qoverpSim = tpr->charge()/sqrt(momentumTP.x()*momentumTP.x()+momentumTP.y()*momentumTP.y()+momentumTP.z()*momentumTP.z());
00674           double thetaSim = momentumTP.theta();
00675           double lambdaSim = M_PI/2-momentumTP.theta();
00676           double phiSim    = momentumTP.phi();
00677           double dxySim    = (-vertexTP.x()*sin(momentumTP.phi())+vertexTP.y()*cos(momentumTP.phi()));
00678           double dzSim     = vertexTP.z() - (vertexTP.x()*momentumTP.x()+vertexTP.y()*momentumTP.y())/sqrt(momentumTP.perp2()) * momentumTP.z()/sqrt(momentumTP.perp2());
00679           
00680           // removed unused variable, left this in case it has side effects
00681           track->parameters();
00682 
00683           double qoverpRec(0);
00684           double qoverpErrorRec(0); 
00685           double ptRec(0);
00686           double ptErrorRec(0);
00687           double lambdaRec(0); 
00688           double lambdaErrorRec(0);
00689           double phiRec(0);
00690           double phiErrorRec(0);
00691 
00692 
00693           //loop to decide whether to take gsfTrack (utilisation of mode-function) or common track
00694           const GsfTrack* gsfTrack(0);
00695           if(useGsf){
00696             gsfTrack = dynamic_cast<const GsfTrack*>(&(*track));
00697             if (gsfTrack==0) edm::LogInfo("MuonTrackValidator") << "Trying to access mode for a non-GsfTrack";
00698           }
00699           
00700           if (gsfTrack) {
00701             // get values from mode
00702             getRecoMomentum(*gsfTrack, ptRec, ptErrorRec, qoverpRec, qoverpErrorRec, 
00703                             lambdaRec,lambdaErrorRec, phiRec, phiErrorRec); 
00704           }
00705          
00706           else {
00707             // get values from track (without mode) 
00708             getRecoMomentum(*track, ptRec, ptErrorRec, qoverpRec, qoverpErrorRec, 
00709                             lambdaRec,lambdaErrorRec, phiRec, phiErrorRec); 
00710           }
00711          
00712           double thetaRec = track->theta();
00713           double ptError = ptErrorRec;
00714           double ptres = ptRec - ptSim; 
00715           double etares = track->eta()-momentumTP.Eta();
00716           double dxyRec    = track->dxy(bs.position());
00717           double dzRec     = track->dz(bs.position());
00718           // eta residue; pt, k, theta, phi, dxy, dz pulls
00719           double qoverpPull=(qoverpRec-qoverpSim)/qoverpErrorRec;
00720           double thetaPull=(lambdaRec-lambdaSim)/lambdaErrorRec;
00721           double phiDiff = phiRec - phiSim;
00722           if (abs(phiDiff) > M_PI) {
00723             if (phiDiff >0.) phiDiff = phiDiff - 2.*M_PI;
00724             else phiDiff = phiDiff + 2.*M_PI;
00725           }
00726           double phiPull=phiDiff/phiErrorRec;
00727           double dxyPull=(dxyRec-dxySim)/track->dxyError();
00728           double dzPull=(dzRec-dzSim)/track->dzError();
00729 
00730           double contrib_Qoverp = ((qoverpRec-qoverpSim)/qoverpErrorRec)*
00731             ((qoverpRec-qoverpSim)/qoverpErrorRec)/5;
00732           double contrib_dxy = ((dxyRec-dxySim)/track->dxyError())*((dxyRec-dxySim)/track->dxyError())/5;
00733           double contrib_dz = ((dzRec-dzSim)/track->dzError())*((dzRec-dzSim)/track->dzError())/5;
00734           double contrib_theta = ((lambdaRec-lambdaSim)/lambdaErrorRec)*
00735             ((lambdaRec-lambdaSim)/lambdaErrorRec)/5;
00736           double contrib_phi = (phiDiff/phiErrorRec)*(phiDiff/phiErrorRec)/5;
00737           
00738           LogTrace("MuonTrackValidator") << "assocChi2=" << tp.begin()->second << "\n"
00739                                          << "" <<  "\n"
00740                                          << "ptREC=" << ptRec << "\n"
00741                                          << "etaREC=" << track->eta() << "\n"
00742                                          << "qoverpREC=" << qoverpRec << "\n"
00743                                          << "dxyREC=" << dxyRec << "\n"
00744                                          << "dzREC=" << dzRec << "\n"
00745                                          << "thetaREC=" << track->theta() << "\n"
00746                                          << "phiREC=" << phiRec << "\n"
00747                                          << "" <<  "\n"
00748                                          << "qoverpError()=" << qoverpErrorRec << "\n"
00749                                          << "dxyError()=" << track->dxyError() << "\n"
00750                                          << "dzError()=" << track->dzError() << "\n"
00751                                          << "thetaError()=" << lambdaErrorRec << "\n"
00752                                          << "phiError()=" << phiErrorRec << "\n"
00753                                          << "" <<  "\n"
00754                                          << "ptSIM=" << ptSim << "\n"
00755                                          << "etaSIM=" << momentumTP.Eta() << "\n"    
00756                                          << "qoverpSIM=" << qoverpSim << "\n"
00757                                          << "dxySIM=" << dxySim << "\n"
00758                                          << "dzSIM=" << dzSim << "\n"
00759                                          << "thetaSIM=" << M_PI/2-lambdaSim << "\n"
00760                                          << "phiSIM=" << phiSim << "\n"
00761                                          << "" << "\n"
00762                                          << "contrib_Qoverp=" << contrib_Qoverp << "\n"
00763                                          << "contrib_dxy=" << contrib_dxy << "\n"
00764                                          << "contrib_dz=" << contrib_dz << "\n"
00765                                          << "contrib_theta=" << contrib_theta << "\n"
00766                                          << "contrib_phi=" << contrib_phi << "\n"
00767                                          << "" << "\n"
00768                                          <<"chi2PULL="<<contrib_Qoverp+contrib_dxy+contrib_dz+contrib_theta+contrib_phi<<"\n";
00769           
00770           h_pullQoverp[w]->Fill(qoverpPull);
00771           h_pullTheta[w]->Fill(thetaPull);
00772           h_pullPhi[w]->Fill(phiPull);
00773           h_pullDxy[w]->Fill(dxyPull);
00774           h_pullDz[w]->Fill(dzPull);
00775 
00776 
00777           h_pt[w]->Fill(ptres/ptError);
00778           h_eta[w]->Fill(etares);
00779           etares_vs_eta[w]->Fill(getEta(track->eta()),etares);
00780  
00781 
00782           //chi2 and #hit vs eta: fill 2D histos
00783           chi2_vs_eta[w]->Fill(getEta(track->eta()),track->normalizedChi2());
00784           nhits_vs_eta[w]->Fill(getEta(track->eta()),track->numberOfValidHits());
00785           nDThits_vs_eta[w]->Fill(getEta(track->eta()),track->hitPattern().numberOfValidMuonDTHits());
00786           nCSChits_vs_eta[w]->Fill(getEta(track->eta()),track->hitPattern().numberOfValidMuonCSCHits());
00787           nRPChits_vs_eta[w]->Fill(getEta(track->eta()),track->hitPattern().numberOfValidMuonRPCHits());
00788 
00789           nlosthits_vs_eta[w]->Fill(getEta(track->eta()),track->numberOfLostHits());
00790 
00791           //resolution of track params: fill 2D histos
00792           dxyres_vs_eta[w]->Fill(getEta(track->eta()),dxyRec-dxySim);
00793           ptres_vs_eta[w]->Fill(getEta(track->eta()),(ptRec-ptSim)/ptRec);
00794           dzres_vs_eta[w]->Fill(getEta(track->eta()),dzRec-dzSim);
00795           phires_vs_eta[w]->Fill(getEta(track->eta()),phiDiff);
00796           cotThetares_vs_eta[w]->Fill(getEta(track->eta()), cos(thetaRec)/sin(thetaRec) - cos(thetaSim)/sin(thetaSim));
00797           
00798           //same as before but vs pT
00799           dxyres_vs_pt[w]->Fill(getPt(ptRec),dxyRec-dxySim);
00800           ptres_vs_pt[w]->Fill(getPt(ptRec),(ptRec-ptSim)/ptRec);
00801           dzres_vs_pt[w]->Fill(getPt(ptRec),dzRec-dzSim);
00802           phires_vs_pt[w]->Fill(getPt(ptRec),phiDiff);
00803           cotThetares_vs_pt[w]->Fill(getPt(ptRec), cos(thetaRec)/sin(thetaRec) - cos(thetaSim)/sin(thetaSim));
00804                  
00805           //pulls of track params vs eta: fill 2D histos
00806           dxypull_vs_eta[w]->Fill(getEta(track->eta()),dxyPull);
00807           ptpull_vs_eta[w]->Fill(getEta(track->eta()),ptres/ptError);
00808           dzpull_vs_eta[w]->Fill(getEta(track->eta()),dzPull);
00809           phipull_vs_eta[w]->Fill(getEta(track->eta()),phiPull);
00810           thetapull_vs_eta[w]->Fill(getEta(track->eta()),thetaPull);
00811 
00812           //plots vs phi
00813           nhits_vs_phi[w]->Fill(phiRec,track->numberOfValidHits());
00814           chi2_vs_phi[w]->Fill(phiRec,track->normalizedChi2());
00815           ptmean_vs_eta_phi[w]->Fill(phiRec,getEta(track->eta()),ptRec);
00816           phimean_vs_eta_phi[w]->Fill(phiRec,getEta(track->eta()),phiRec);
00817           ptres_vs_phi[w]->Fill(phiRec,(ptRec-ptSim)/ptRec);
00818           phires_vs_phi[w]->Fill(phiRec,phiDiff);
00819           ptpull_vs_phi[w]->Fill(phiRec,ptres/ptError);
00820           phipull_vs_phi[w]->Fill(phiRec,phiPull); 
00821           thetapull_vs_phi[w]->Fill(phiRec,thetaPull); 
00822           
00823           std::vector<PSimHit> simhits;
00824           
00825           if (usetracker && usemuon) {
00826             simhits=tpr.get()->trackPSimHit();
00827           } 
00828           else if (!usetracker && usemuon) {
00829             simhits=tpr.get()->trackPSimHit(DetId::Muon);
00830           }
00831           else if (usetracker && !usemuon) {
00832             simhits=tpr.get()->trackPSimHit(DetId::Tracker);
00833           }
00834           
00835           nrecHit_vs_nsimHit_rec2sim[w]->Fill(track->numberOfValidHits(), (int)(simhits.end()-simhits.begin() ));
00836           
00837         } // End of try{
00838         catch (cms::Exception e){
00839           LogTrace("MuonTrackValidator") << "exception found: " << e.what() << "\n";
00840         }
00841       } // End of for(View<Track>::size_type i=0; i<trackCollectionSize; ++i){
00842       if (at!=0) h_tracks[w]->Fill(at);
00843       h_fakes[w]->Fill(rT-at);
00844       edm::LogVerbatim("MuonTrackValidator") << "Total Simulated: " << st << "\n"
00845                                              << "Total Associated (simToReco): " << ats << "\n"
00846                                              << "Total Reconstructed: " << rT << "\n"
00847                                              << "Total Associated (recoToSim): " << at << "\n"
00848                                              << "Total Fakes: " << rT-at << "\n";
00849       nrec_vs_nsim[w]->Fill(rT,st);
00850       w++;
00851     } // End of  for (unsigned int www=0;www<label.size();www++){
00852   } //END of for (unsigned int ww=0;ww<associators.size();ww++){
00853 }
00854 
00855 void MuonTrackValidator::endRun(Run const&, EventSetup const&) {
00856 
00857   int w=0;
00858   for (unsigned int ww=0;ww<associators.size();ww++){
00859     for (unsigned int www=0;www<label.size();www++){
00860 
00861       //chi2 and #hit vs eta: get mean from 2D histos
00862       doProfileX(chi2_vs_eta[w],h_chi2meanh[w]);
00863       doProfileX(nhits_vs_eta[w],h_hits_eta[w]);    
00864       doProfileX(nDThits_vs_eta[w],h_DThits_eta[w]);    
00865       doProfileX(nCSChits_vs_eta[w],h_CSChits_eta[w]);    
00866       doProfileX(nRPChits_vs_eta[w],h_RPChits_eta[w]);    
00867 
00868       doProfileX(nlosthits_vs_eta[w],h_losthits_eta[w]);    
00869       //vs phi
00870       doProfileX(chi2_vs_nhits[w],h_chi2meanhitsh[w]); 
00871       doProfileX(chi2_vs_phi[w],h_chi2mean_vs_phi[w]);
00872       doProfileX(nhits_vs_phi[w],h_hits_phi[w]);
00873 
00874       fillPlotFromVector(h_recoeta[w],totRECeta[w]);
00875       fillPlotFromVector(h_simuleta[w],totSIMeta[w]);
00876       fillPlotFromVector(h_assoceta[w],totASSeta[w]);
00877       fillPlotFromVector(h_assoc2eta[w],totASS2eta[w]);
00878 
00879       fillPlotFromVector(h_recopT[w],totRECpT[w]);
00880       fillPlotFromVector(h_simulpT[w],totSIMpT[w]);
00881       fillPlotFromVector(h_assocpT[w],totASSpT[w]);
00882       fillPlotFromVector(h_assoc2pT[w],totASS2pT[w]);
00883 
00884       fillPlotFromVector(h_recohit[w],totREC_hit[w]);
00885       fillPlotFromVector(h_simulhit[w],totSIM_hit[w]);
00886       fillPlotFromVector(h_assochit[w],totASS_hit[w]);
00887       fillPlotFromVector(h_assoc2hit[w],totASS2_hit[w]);
00888 
00889       fillPlotFromVector(h_recophi[w],totREC_phi[w]);
00890       fillPlotFromVector(h_simulphi[w],totSIM_phi[w]);
00891       fillPlotFromVector(h_assocphi[w],totASS_phi[w]);
00892       fillPlotFromVector(h_assoc2phi[w],totASS2_phi[w]);
00893 
00894       fillPlotFromVector(h_recodxy[w],totREC_dxy[w]);
00895       fillPlotFromVector(h_simuldxy[w],totSIM_dxy[w]);
00896       fillPlotFromVector(h_assocdxy[w],totASS_dxy[w]);
00897       fillPlotFromVector(h_assoc2dxy[w],totASS2_dxy[w]);
00898 
00899       fillPlotFromVector(h_recodz[w],totREC_dz[w]);
00900       fillPlotFromVector(h_simuldz[w],totSIM_dz[w]);
00901       fillPlotFromVector(h_assocdz[w],totASS_dz[w]);
00902       fillPlotFromVector(h_assoc2dz[w],totASS2_dz[w]);
00903 
00904       fillPlotFromVector(h_simulvertpos[w],totSIM_vertpos[w]);
00905       fillPlotFromVector(h_assocvertpos[w],totASS_vertpos[w]);
00906 
00907       fillPlotFromVector(h_simulzpos[w],totSIM_zpos[w]);
00908       fillPlotFromVector(h_assoczpos[w],totASS_zpos[w]);
00909       
00910       if (MABH) {
00911         fillPlotFromVector(h_assoceta_Quality05[w] ,totASSeta_Quality05[w]);
00912         fillPlotFromVector(h_assoceta_Quality075[w],totASSeta_Quality075[w]);
00913         fillPlotFromVector(h_assocpT_Quality05[w] ,totASSpT_Quality05[w]);
00914         fillPlotFromVector(h_assocpT_Quality075[w],totASSpT_Quality075[w]);
00915         fillPlotFromVector(h_assocphi_Quality05[w] ,totASS_phi_Quality05[w]);
00916         fillPlotFromVector(h_assocphi_Quality075[w],totASS_phi_Quality075[w]);
00917       }
00918       
00919       w++;
00920     }
00921   }
00922   
00923   if ( out.size() != 0 && dbe_ ) dbe_->save(out);
00924 }
00925 
00926 
00927 void 
00928 MuonTrackValidator::getRecoMomentum (const reco::Track& track, double& pt, double& ptError,
00929                                       double& qoverp, double& qoverpError, double& lambda,double& lambdaError,  double& phi, double& phiError ) const {
00930   pt = track.pt();
00931   ptError = track.ptError();
00932   qoverp = track.qoverp();
00933   qoverpError = track.qoverpError();
00934   lambda = track.lambda();
00935   lambdaError = track.lambdaError(); 
00936   phi = track.phi(); 
00937   phiError = track.phiError();
00938 
00939 }
00940 
00941 void 
00942 MuonTrackValidator::getRecoMomentum (const reco::GsfTrack& gsfTrack, double& pt, double& ptError,
00943                                       double& qoverp, double& qoverpError, double& lambda,double& lambdaError,  double& phi, double& phiError  ) const {
00944 
00945   pt = gsfTrack.ptMode();
00946   ptError = gsfTrack.ptModeError();
00947   qoverp = gsfTrack.qoverpMode();
00948   qoverpError = gsfTrack.qoverpModeError();
00949   lambda = gsfTrack.lambdaMode();
00950   lambdaError = gsfTrack.lambdaModeError(); 
00951   phi = gsfTrack.phiMode(); 
00952   phiError = gsfTrack.phiModeError();
00953 
00954 }
00955
/data/refman/pasoursint/CMSSW_5_3_10/src/Validation/RecoMuon/plugins/MuonTrackValidator.cc
