CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonTrackValidator.cc
Go to the documentation of this file.
1 #include "Validation/RecoMuon/plugins/MuonTrackValidator.h"
2 #include "DQMServices/ClientConfig/interface/FitSlicesYTool.h"
3 
4 #include "FWCore/Framework/interface/MakerMacros.h"
5 #include "FWCore/MessageLogger/interface/MessageLogger.h"
6 
7 #include "DataFormats/TrackReco/interface/Track.h"
8 #include "DataFormats/TrackReco/interface/TrackFwd.h"
9 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
10 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
11 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
12 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
13 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
14 #include "SimTracker/Records/interface/TrackAssociatorRecord.h"
15 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
16 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
17 #include "SimTracker/TrackAssociation/plugins/ParametersDefinerForTPESProducer.h"
18 #include "SimTracker/TrackAssociation/plugins/CosmicParametersDefinerForTPESProducer.h"
19 
20 #include "TMath.h"
21 #include <TF1.h>
22 
23 using namespace std;
24 using namespace edm;
25 
26 void MuonTrackValidator::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const&, edm::EventSetup const& setup) {
27 
28  int j=0;
29  for (unsigned int ww=0;ww<associators.size();ww++){
30  for (unsigned int www=0;www<label.size();www++){
31 
32  ibooker.cd();
33  InputTag algo = label[www];
34  string dirName=dirName_;
35  if (algo.process()!="")
36  dirName+=algo.process()+"_";
37  if(algo.label()!="")
38  dirName+=algo.label();
39  if(algo.instance()!="")
40  dirName+=("_"+algo.instance());
41  if (dirName.find("Tracks")<dirName.length()){
42  dirName.replace(dirName.find("Tracks"),6,"Trks");
43  }
44  if (dirName.find("UpdatedAtVtx")<dirName.length()){
45  dirName.replace(dirName.find("UpdatedAtVtx"),12,"UpdAtVtx");
46  }
47  string assoc= associators[ww];
48  if (assoc.find("tpToTkmuTrackAssociation")<assoc.length()){
49  dirName+="_TkAsso";
50  }
51  std::replace(dirName.begin(), dirName.end(), ':', '_');
52  ibooker.setCurrentFolder(dirName.c_str());
53 
54  setUpVectors();
55 
56  ibooker.goUp();
57  string subDirName = dirName + "/simulation";
58  ibooker.setCurrentFolder(subDirName.c_str());
59  h_ptSIM.push_back( ibooker.book1D("ptSIM", "generated p_{t}", 5500, 0, 110 ) );
60  h_etaSIM.push_back( ibooker.book1D("etaSIM", "generated pseudorapidity", 500, -2.5, 2.5 ) );
61  h_tracksSIM.push_back( ibooker.book1D("tracksSIM","number of simulated tracks",200,-0.5,99.5) );
62  h_vertposSIM.push_back( ibooker.book1D("vertposSIM","Transverse position of sim vertices",100,0.,120.) );
63 
64  ibooker.cd();
65  ibooker.setCurrentFolder(dirName.c_str());
66  h_tracks.push_back( ibooker.book1D("tracks","number of reconstructed tracks",200,-0.5,19.5) );
67  h_fakes.push_back( ibooker.book1D("fakes","number of fake reco tracks",20,-0.5,19.5) );
68  h_charge.push_back( ibooker.book1D("charge","charge",3,-1.5,1.5) );
69  h_hits.push_back( ibooker.book1D("hits", "number of hits per track", nintHit,minHit,maxHit ) );
70  h_losthits.push_back( ibooker.book1D("losthits", "number of lost hits per track", nintHit,minHit,maxHit) );
71  h_nchi2.push_back( ibooker.book1D("chi2", "normalized #chi^{2}", 200, 0, 20 ) );
72  h_nchi2_prob.push_back( ibooker.book1D("chi2_prob", "normalized #chi^{2} probability",100,0,1));
73 
75  h_recoeta.push_back( ibooker.book1D("num_reco_eta","N of reco track vs eta",nint,min,max) );
76  h_assoceta.push_back( ibooker.book1D("num_assoc(simToReco)_eta","N of associated tracks (simToReco) vs eta",nint,min,max) );
77  h_assoc2eta.push_back( ibooker.book1D("num_assoc(recoToSim)_eta","N of associated (recoToSim) tracks vs eta",nint,min,max) );
78  h_simuleta.push_back( ibooker.book1D("num_simul_eta","N of simulated tracks vs eta",nint,min,max) );
79  h_recopT.push_back( ibooker.book1D("num_reco_pT","N of reco track vs pT",nintpT,minpT,maxpT) );
80  h_assocpT.push_back( ibooker.book1D("num_assoc(simToReco)_pT","N of associated tracks (simToReco) vs pT",nintpT,minpT,maxpT) );
81  h_assoc2pT.push_back( ibooker.book1D("num_assoc(recoToSim)_pT","N of associated (recoToSim) tracks vs pT",nintpT,minpT,maxpT) );
82  h_simulpT.push_back( ibooker.book1D("num_simul_pT","N of simulated tracks vs pT",nintpT,minpT,maxpT) );
83  //
84  h_recohit.push_back( ibooker.book1D("num_reco_hit","N of reco track vs hit",nintHit,minHit,maxHit) );
85  h_assochit.push_back( ibooker.book1D("num_assoc(simToReco)_hit","N of associated tracks (simToReco) vs hit",nintHit,minHit,maxHit) );
86  h_assoc2hit.push_back( ibooker.book1D("num_assoc(recoToSim)_hit","N of associated (recoToSim) tracks vs hit",nintHit,minHit,maxHit) );
87  h_simulhit.push_back( ibooker.book1D("num_simul_hit","N of simulated tracks vs hit",nintHit,minHit,maxHit) );
88  //
89  h_recophi.push_back( ibooker.book1D("num_reco_phi","N of reco track vs phi",nintPhi,minPhi,maxPhi) );
90  h_assocphi.push_back( ibooker.book1D("num_assoc(simToReco)_phi","N of associated tracks (simToReco) vs phi",nintPhi,minPhi,maxPhi) );
91  h_assoc2phi.push_back( ibooker.book1D("num_assoc(recoToSim)_phi","N of associated (recoToSim) tracks vs phi",nintPhi,minPhi,maxPhi) );
92  h_simulphi.push_back( ibooker.book1D("num_simul_phi","N of simulated tracks vs phi",nintPhi,minPhi,maxPhi) );
93 
94  h_recodxy.push_back( ibooker.book1D("num_reco_dxy","N of reco track vs dxy",nintDxy,minDxy,maxDxy) );
95  h_assocdxy.push_back( ibooker.book1D("num_assoc(simToReco)_dxy","N of associated tracks (simToReco) vs dxy",nintDxy,minDxy,maxDxy) );
96  h_assoc2dxy.push_back( ibooker.book1D("num_assoc(recoToSim)_dxy","N of associated (recoToSim) tracks vs dxy",nintDxy,minDxy,maxDxy) );
97  h_simuldxy.push_back( ibooker.book1D("num_simul_dxy","N of simulated tracks vs dxy",nintDxy,minDxy,maxDxy) );
98 
99  h_recodz.push_back( ibooker.book1D("num_reco_dz","N of reco track vs dz",nintDz,minDz,maxDz) );
100  h_assocdz.push_back( ibooker.book1D("num_assoc(simToReco)_dz","N of associated tracks (simToReco) vs dz",nintDz,minDz,maxDz) );
101  h_assoc2dz.push_back( ibooker.book1D("num_assoc(recoToSim)_dz","N of associated (recoToSim) tracks vs dz",nintDz,minDz,maxDz) );
102  h_simuldz.push_back( ibooker.book1D("num_simul_dz","N of simulated tracks vs dz",nintDz,minDz,maxDz) );
103 
104  h_assocvertpos.push_back( ibooker.book1D("num_assoc(simToReco)_vertpos","N of associated tracks (simToReco) vs transverse vert position",nintVertpos,minVertpos,maxVertpos) );
105  h_simulvertpos.push_back( ibooker.book1D("num_simul_vertpos","N of simulated tracks vs transverse vert position",nintVertpos,minVertpos,maxVertpos) );
106 
107  h_assoczpos.push_back( ibooker.book1D("num_assoc(simToReco)_zpos","N of associated tracks (simToReco) vs z vert position",nintZpos,minZpos,maxZpos) );
108  h_simulzpos.push_back( ibooker.book1D("num_simul_zpos","N of simulated tracks vs z vert position",nintZpos,minZpos,maxZpos) );
109 
111 
112  h_eta.push_back( ibooker.book1D("eta", "pseudorapidity residue", 1000, -0.1, 0.1 ) );
113  h_pt.push_back( ibooker.book1D("pullPt", "pull of p_{t}", 100, -10, 10 ) );
114  h_pullTheta.push_back( ibooker.book1D("pullTheta","pull of #theta parameter",250,-25,25) );
115  h_pullPhi.push_back( ibooker.book1D("pullPhi","pull of #phi parameter",250,-25,25) );
116  h_pullDxy.push_back( ibooker.book1D("pullDxy","pull of dxy parameter",250,-25,25) );
117  h_pullDz.push_back( ibooker.book1D("pullDz","pull of dz parameter",250,-25,25) );
118  h_pullQoverp.push_back( ibooker.book1D("pullQoverp","pull of qoverp parameter",250,-25,25) );
119 
120  if (associators[ww]=="trackAssociatorByChi2"){
121  h_assochi2.push_back( ibooker.book1D("assocChi2","track association #chi^{2}",1000000,0,100000) );
122  h_assochi2_prob.push_back(ibooker.book1D("assocChi2_prob","probability of association #chi^{2}",100,0,1));
123  } else if (associators[ww]=="trackAssociatorByHits"){
124  h_assocFraction.push_back( ibooker.book1D("assocFraction","fraction of shared hits",200,0,2) );
125  h_assocSharedHit.push_back(ibooker.book1D("assocSharedHit","number of shared hits",20,0,20));
126  }
127 
128  chi2_vs_nhits.push_back( ibooker.book2D("chi2_vs_nhits","#chi^{2} vs nhits",25,0,25,100,0,10) );
129  h_chi2meanhitsh.push_back( ibooker.bookProfile("chi2mean_vs_nhits","mean #chi^{2} vs nhits",25,0,25,100,0,10) );
130 
131  etares_vs_eta.push_back( ibooker.book2D("etares_vs_eta","etaresidue vs eta",nint,min,max,200,-0.1,0.1) );
132  nrec_vs_nsim.push_back( ibooker.book2D("nrec_vs_nsim","nrec vs nsim",20,-0.5,19.5,20,-0.5,19.5) );
133 
134  chi2_vs_eta.push_back( ibooker.book2D("chi2_vs_eta","chi2_vs_eta",nint,min,max, 200, 0, 20 ));
135  h_chi2meanh.push_back( ibooker.bookProfile("chi2mean","mean #chi^{2} vs #eta",nint,min,max, 200, 0, 20) );
136  chi2_vs_phi.push_back( ibooker.book2D("chi2_vs_phi","#chi^{2} vs #phi",nintPhi,minPhi,maxPhi, 200, 0, 20 ) );
137  h_chi2mean_vs_phi.push_back( ibooker.bookProfile("chi2mean_vs_phi","mean of #chi^{2} vs #phi",nintPhi,minPhi,maxPhi, 200, 0, 20) );
138 
139  nhits_vs_eta.push_back( ibooker.book2D("nhits_vs_eta","nhits vs eta",nint,min,max,nintHit,minHit,maxHit) );
140  nDThits_vs_eta.push_back( ibooker.book2D("nDThits_vs_eta","# DT hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
141  nCSChits_vs_eta.push_back( ibooker.book2D("nCSChits_vs_eta","# CSC hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
142  nRPChits_vs_eta.push_back( ibooker.book2D("nRPChits_vs_eta","# RPC hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
143  if(useGEMs_) nGEMhits_vs_eta.push_back( ibooker.book2D("nGEMhits_vs_eta","# GEM hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
144  if(useME0_) nME0hits_vs_eta.push_back( ibooker.book2D("nME0hits_vs_eta","# ME0 hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
145 
146  h_DThits_eta.push_back( ibooker.bookProfile("DThits_eta","mean # DT hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
147  h_CSChits_eta.push_back( ibooker.bookProfile("CSChits_eta","mean # CSC hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
148  h_RPChits_eta.push_back( ibooker.bookProfile("RPChits_eta","mean # RPC hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
149  if(useGEMs_) h_GEMhits_eta.push_back( ibooker.bookProfile("GEMhits_eta","mean # GEM hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
150  if(useME0_) h_ME0hits_eta.push_back( ibooker.bookProfile("ME0hits_eta","mean # ME0 hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
151 
152  h_hits_eta.push_back( ibooker.bookProfile("hits_eta","mean #hits vs eta",nint,min,max,nintHit,minHit,maxHit) );
153  nhits_vs_phi.push_back( ibooker.book2D("nhits_vs_phi","#hits vs #phi",nintPhi,minPhi,maxPhi,nintHit,minHit,maxHit) );
154  h_hits_phi.push_back( ibooker.bookProfile("hits_phi","mean #hits vs #phi",nintPhi,minPhi,maxPhi, nintHit,minHit,maxHit) );
155 
156  nlosthits_vs_eta.push_back( ibooker.book2D("nlosthits_vs_eta","nlosthits vs eta",nint,min,max,nintHit,minHit,maxHit) );
157  h_losthits_eta.push_back( ibooker.bookProfile("losthits_eta","losthits_eta",nint,min,max,nintHit,minHit,maxHit) );
158 
159  ptres_vs_eta.push_back(ibooker.book2D("ptres_vs_eta","ptres_vs_eta",nint,min,max, ptRes_nbin, ptRes_rangeMin, ptRes_rangeMax));
160  ptres_vs_phi.push_back( ibooker.book2D("ptres_vs_phi","p_{t} res vs #phi",nintPhi,minPhi,maxPhi, ptRes_nbin, ptRes_rangeMin, ptRes_rangeMax));
161  ptres_vs_pt.push_back(ibooker.book2D("ptres_vs_pt","ptres_vs_pt",nintpT,minpT,maxpT, ptRes_nbin, ptRes_rangeMin, ptRes_rangeMax));
162 
163  cotThetares_vs_eta.push_back(ibooker.book2D("cotThetares_vs_eta","cotThetares_vs_eta",nint,min,max,cotThetaRes_nbin, cotThetaRes_rangeMin, cotThetaRes_rangeMax));
164  cotThetares_vs_pt.push_back(ibooker.book2D("cotThetares_vs_pt","cotThetares_vs_pt",nintpT,minpT,maxpT, cotThetaRes_nbin, cotThetaRes_rangeMin, cotThetaRes_rangeMax));
165 
166  phires_vs_eta.push_back(ibooker.book2D("phires_vs_eta","phires_vs_eta",nint,min,max, phiRes_nbin, phiRes_rangeMin, phiRes_rangeMax));
167  phires_vs_pt.push_back(ibooker.book2D("phires_vs_pt","phires_vs_pt",nintpT,minpT,maxpT, phiRes_nbin, phiRes_rangeMin, phiRes_rangeMax));
168  phires_vs_phi.push_back(ibooker.book2D("phires_vs_phi","#phi res vs #phi",nintPhi,minPhi,maxPhi,phiRes_nbin, phiRes_rangeMin, phiRes_rangeMax));
169 
170  dxyres_vs_eta.push_back(ibooker.book2D("dxyres_vs_eta","dxyres_vs_eta",nint,min,max,dxyRes_nbin, dxyRes_rangeMin, dxyRes_rangeMax));
171  dxyres_vs_pt.push_back( ibooker.book2D("dxyres_vs_pt","dxyres_vs_pt",nintpT,minpT,maxpT,dxyRes_nbin, dxyRes_rangeMin, dxyRes_rangeMax));
172 
173  dzres_vs_eta.push_back(ibooker.book2D("dzres_vs_eta","dzres_vs_eta",nint,min,max,dzRes_nbin, dzRes_rangeMin, dzRes_rangeMax));
174  dzres_vs_pt.push_back(ibooker.book2D("dzres_vs_pt","dzres_vs_pt",nintpT,minpT,maxpT,dzRes_nbin, dzRes_rangeMin, dzRes_rangeMax));
175 
176  ptmean_vs_eta_phi.push_back(ibooker.bookProfile2D("ptmean_vs_eta_phi","mean p_{t} vs #eta and #phi",nintPhi,minPhi,maxPhi,nint,min,max,1000,0,1000));
177  phimean_vs_eta_phi.push_back(ibooker.bookProfile2D("phimean_vs_eta_phi","mean #phi vs #eta and #phi",nintPhi,minPhi,maxPhi,nint,min,max,nintPhi,minPhi,maxPhi));
178 
179  //pulls of track params vs eta: to be used with fitslicesytool
180  dxypull_vs_eta.push_back(ibooker.book2D("dxypull_vs_eta","dxypull_vs_eta",nint,min,max,100,-10,10));
181  ptpull_vs_eta.push_back(ibooker.book2D("ptpull_vs_eta","ptpull_vs_eta",nint,min,max,100,-10,10));
182  dzpull_vs_eta.push_back(ibooker.book2D("dzpull_vs_eta","dzpull_vs_eta",nint,min,max,100,-10,10));
183  phipull_vs_eta.push_back(ibooker.book2D("phipull_vs_eta","phipull_vs_eta",nint,min,max,100,-10,10));
184  thetapull_vs_eta.push_back(ibooker.book2D("thetapull_vs_eta","thetapull_vs_eta",nint,min,max,100,-10,10));
185 
186  //pulls of track params vs phi
187  ptpull_vs_phi.push_back(ibooker.book2D("ptpull_vs_phi","p_{t} pull vs #phi",nintPhi,minPhi,maxPhi,100,-10,10));
188  phipull_vs_phi.push_back(ibooker.book2D("phipull_vs_phi","#phi pull vs #phi",nintPhi,minPhi,maxPhi,100,-10,10));
189  thetapull_vs_phi.push_back(ibooker.book2D("thetapull_vs_phi","#theta pull vs #phi",nintPhi,minPhi,maxPhi,100,-10,10));
190 
191  nrecHit_vs_nsimHit_sim2rec.push_back( ibooker.book2D("nrecHit_vs_nsimHit_sim2rec","nrecHit vs nsimHit (Sim2RecAssoc)",nintHit,minHit,maxHit, nintHit,minHit,maxHit ));
192  nrecHit_vs_nsimHit_rec2sim.push_back( ibooker.book2D("nrecHit_vs_nsimHit_rec2sim","nrecHit vs nsimHit (Rec2simAssoc)",nintHit,minHit,maxHit, nintHit,minHit,maxHit ));
193 
194  if (MABH) {
195  h_PurityVsQuality.push_back
196  (ibooker.book2D("PurityVsQuality","Purity vs Quality (MABH)",20,0.01,1.01,20,0.01,1.01) );
197  h_assoceta_Quality05.push_back
198  (ibooker.book1D("num_assoc(simToReco)_eta_Q05","N of associated tracks (simToReco) vs eta (Quality>0.5)",nint,min,max) );
199  h_assoceta_Quality075.push_back
200  (ibooker.book1D("num_assoc(simToReco)_eta_Q075","N of associated tracks (simToReco) vs eta (Quality>0.75)",nint,min,max) );
201  h_assocpT_Quality05.push_back
202  (ibooker.book1D("num_assoc(simToReco)_pT_Q05","N of associated tracks (simToReco) vs pT (Quality>0.5)",nintpT,minpT,maxpT) );
203  h_assocpT_Quality075.push_back
204  (ibooker.book1D("num_assoc(simToReco)_pT_Q075","N of associated tracks (simToReco) vs pT (Quality>0.75)",nintpT,minpT,maxpT) );
205  h_assocphi_Quality05.push_back
206  (ibooker.book1D("num_assoc(simToReco)_phi_Q05","N of associated tracks (simToReco) vs phi (Quality>0.5)",nintPhi,minPhi,maxPhi) );
207  h_assocphi_Quality075.push_back
208  (ibooker.book1D("num_assoc(simToReco)_phi_Q075","N of associated tracks (simToReco) vs phi (Quality>0.75)",nintPhi,minPhi,maxPhi) );
209  }
210 
211  if(useLogPt){
212  BinLogX(dzres_vs_pt[j]->getTH2F());
213  BinLogX(dxyres_vs_pt[j]->getTH2F());
214  BinLogX(phires_vs_pt[j]->getTH2F());
215  BinLogX(cotThetares_vs_pt[j]->getTH2F());
216  BinLogX(ptres_vs_pt[j]->getTH2F());
217  BinLogX(h_recopT[j]->getTH1F());
218  BinLogX(h_assocpT[j]->getTH1F());
219  BinLogX(h_assoc2pT[j]->getTH1F());
220  BinLogX(h_simulpT[j]->getTH1F());
221  if (MABH) {
222  BinLogX(h_assocpT_Quality05[j]->getTH1F());
223  BinLogX(h_assocpT_Quality075[j]->getTH1F());
224  }
225  j++;
226  }
227 
228  }
229  }
230 }
231 
232 void MuonTrackValidator::analyze(const edm::Event& event, const edm::EventSetup& setup){
233  using namespace reco;
234 
235  edm::LogInfo("MuonTrackValidator") << "\n====================================================" << "\n"
236  << "Analyzing new event" << "\n"
237  << "====================================================\n" << "\n";
238 
239  edm::ESHandle<ParametersDefinerForTP> Lhc_parametersDefinerTP;
240  std::unique_ptr<ParametersDefinerForTP> Cosmic_parametersDefinerTP;
241 
242  if(parametersDefiner=="LhcParametersDefinerForTP") {
243  setup.get<TrackAssociatorRecord>().get(parametersDefiner, Lhc_parametersDefinerTP);
244  }
245  else if(parametersDefiner=="CosmicParametersDefinerForTP") {
246  edm::ESHandle<CosmicParametersDefinerForTP> _Cosmic_parametersDefinerTP;
247  setup.get<TrackAssociatorRecord>().get(parametersDefiner, _Cosmic_parametersDefinerTP);
248 
249  //Since we modify the object, we must clone it
250  Cosmic_parametersDefinerTP = _Cosmic_parametersDefinerTP->clone();
251 
252  edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
253  //warning: make sure the TP collection used in the map is the same used here
254  event.getByToken(_simHitTpMapTag,simHitsTPAssoc);
255  Cosmic_parametersDefinerTP->initEvent(simHitsTPAssoc);
256  cosmictpSelector.initEvent(simHitsTPAssoc);
257  }
258  else {
259  edm::LogError("MuonTrackValidator")
260  << "Unexpected label: parametersDefiner = "<< parametersDefiner.c_str() << "\n";
261  }
262 
263  edm::Handle<TrackingParticleCollection> TPCollectionHeff ;
264  event.getByToken(tp_effic_Token,TPCollectionHeff);
265  TrackingParticleCollection const & tPCeff = *(TPCollectionHeff.product());
266 
267  edm::Handle<TrackingParticleCollection> TPCollectionHfake ;
268  event.getByToken(tp_fake_Token,TPCollectionHfake);
269 
270  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
271  event.getByToken(bsSrc_Token,recoBeamSpotHandle);
272  reco::BeamSpot bs = *recoBeamSpotHandle;
273 
274  std::vector<const reco::TrackToTrackingParticleAssociator*> associator;
275  if (UseAssociators) {
276  edm::Handle<reco::TrackToTrackingParticleAssociator> theAssociator;
277  for (unsigned int w=0;w<associators.size();w++) {
278  event.getByLabel(associators[w],theAssociator);
279  associator.push_back( theAssociator.product() );
280  }
281  }
282 
283  int w=0;
284  for (unsigned int ww=0;ww<associators.size();ww++){
285  for (unsigned int www=0;www<label.size();www++){
286  //
287  //get collections from the event
288  //
289  edm::Handle<edm::View<Track> > trackCollection;
290 
291  reco::RecoToSimCollection recSimColl;
292  reco::SimToRecoCollection simRecColl;
293  unsigned int trackCollectionSize = 0;
294 
295  if(!event.getByToken(track_Collection_Token[www], trackCollection)&&ignoremissingtkcollection_) {
296 
297  recSimColl.post_insert();
298  simRecColl.post_insert();
299 
300  }
301 
302  else {
303 
304  trackCollectionSize = trackCollection->size();
305 
306  //associate tracks
307  if(UseAssociators){
308  edm::LogVerbatim("MuonTrackValidator") << "Analyzing "
309  << label[www].process()<<":"
310  << label[www].label()<<":"
311  << label[www].instance()<<" with "
312  << associators[ww].c_str() <<"\n";
313 
314  LogTrace("MuonTrackValidator") << "Calling associateRecoToSim method" << "\n";
315  recSimColl=associator[ww]->associateRecoToSim(trackCollection,
316  TPCollectionHfake);
317  LogTrace("MuonTrackValidator") << "Calling associateSimToReco method" << "\n";
318  simRecColl=associator[ww]->associateSimToReco(trackCollection,
319  TPCollectionHeff);
320  }
321  else{
322  edm::LogVerbatim("MuonTrackValidator") << "Analyzing "
323  << label[www].process()<<":"
324  << label[www].label()<<":"
325  << label[www].instance()<<" with "
326  << associatormap.process()<<":"
327  << associatormap.label()<<":"
328  << associatormap.instance()<<"\n";
329 
330  Handle<reco::SimToRecoCollection > simtorecoCollectionH;
331  event.getByToken(simToRecoCollection_Token,simtorecoCollectionH);
332  simRecColl= *(simtorecoCollectionH.product());
333 
334  Handle<reco::RecoToSimCollection > recotosimCollectionH;
335  event.getByToken(recoToSimCollection_Token,recotosimCollectionH);
336  recSimColl= *(recotosimCollectionH.product());
337  }
338 
339  }
340 
341  //
342  //fill simulation histograms
343  //compute number of tracks per eta interval
344  //
345  edm::LogVerbatim("MuonTrackValidator") << "\n# of TrackingParticles: " << tPCeff.size() << "\n";
346  int ats = 0;
347  int st = 0;
348  for (TrackingParticleCollection::size_type i=0; i<tPCeff.size(); i++){
349  bool TP_is_matched = false;
350  double quality = 0.;
351  bool Quality05 = false;
352  bool Quality075 = false;
353 
354  TrackingParticleRef tpr(TPCollectionHeff, i);
355  TrackingParticle* tp = const_cast<TrackingParticle*>(tpr.get());
356 
357  TrackingParticle::Vector momentumTP;
358  TrackingParticle::Point vertexTP;
359  double dxySim = 0;
360  double dzSim = 0;
361 
362  //If the TrackingParticle is collision-like, get the momentum and vertex at production state
363  //and the impact parameters w.r.t. PCA
364  if(parametersDefiner=="LhcParametersDefinerForTP")
365  {
366  LogTrace("MuonTrackValidator") <<"TrackingParticle "<< i;
367  if(! tpSelector(*tp)) continue;
368  momentumTP = tp->momentum();
369  vertexTP = tp->vertex();
370  TrackingParticle::Vector momentum = Lhc_parametersDefinerTP->momentum(event,setup,tpr);
371  TrackingParticle::Point vertex = Lhc_parametersDefinerTP->vertex(event,setup,tpr);
372  dxySim = (-vertex.x()*sin(momentum.phi())+vertex.y()*cos(momentum.phi()));
373  dzSim = vertex.z() - (vertex.x()*momentum.x()+vertex.y()*momentum.y()) /
374  sqrt(momentum.perp2()) * momentum.z()/sqrt(momentum.perp2());
375  }
376  //for cosmics get the momentum and vertex at PCA
377  else if(parametersDefiner=="CosmicParametersDefinerForTP")
378  {
379  edm::LogVerbatim("MuonTrackValidator") <<"TrackingParticle "<< i;
380  if(! cosmictpSelector(tpr,&bs,event,setup)) continue;
381  momentumTP = Cosmic_parametersDefinerTP->momentum(event,setup,tpr);
382  vertexTP = Cosmic_parametersDefinerTP->vertex(event,setup,tpr);
383  dxySim = (-vertexTP.x()*sin(momentumTP.phi())+vertexTP.y()*cos(momentumTP.phi()));
384  dzSim = vertexTP.z() - (vertexTP.x()*momentumTP.x()+vertexTP.y()*momentumTP.y()) /
385  sqrt(momentumTP.perp2()) * momentumTP.z()/sqrt(momentumTP.perp2());
386  }
387  edm::LogVerbatim("MuonTrackValidator") <<"--------------------Selected TrackingParticle #"<<tpr.key();
388  edm::LogVerbatim("MuonTrackValidator") <<"momentumTP: pt = "<<sqrt(momentumTP.perp2())<<", pz = "<<momentumTP.z()
389  <<", \t vertexTP: radius = "<<sqrt(vertexTP.perp2())<< ", z = "<<vertexTP.z() <<"\n";
390  st++;
391 
392  h_ptSIM[w]->Fill(sqrt(momentumTP.perp2()));
393  h_etaSIM[w]->Fill(momentumTP.eta());
394  h_vertposSIM[w]->Fill(sqrt(vertexTP.perp2()));
395 
396  std::vector<std::pair<RefToBase<Track>, double> > rt;
397  if(simRecColl.find(tpr) != simRecColl.end()){
398  rt = (std::vector<std::pair<RefToBase<Track>, double> >) simRecColl[tpr];
399  if (rt.size()!=0) {
400  RefToBase<Track> assoc_recoTrack = rt.begin()->first;
401  edm::LogVerbatim("MuonTrackValidator")<<"-----------------------------associated Track #"<<assoc_recoTrack.key();
402  TP_is_matched = true;
403  ats++;
404  quality = rt.begin()->second;
405  edm::LogVerbatim("MuonTrackValidator") << "TrackingParticle #" <<tpr.key()
406  << " with pt=" << sqrt(momentumTP.perp2())
407  << " associated with quality:" << quality <<"\n";
408  if (MABH) {
409  if (quality > 0.75) {
410  Quality075 = true;
411  Quality05 = true;
412  }
413  else if (quality > 0.5) {
414  Quality05 = true;
415  }
416  }
417  }
418  }else{
419  edm::LogVerbatim("MuonTrackValidator")
420  << "TrackingParticle #" << tpr.key()
421  << " with pt,eta,phi: "
422  << sqrt(momentumTP.perp2()) << " , "
423  << momentumTP.eta() << " , "
424  << momentumTP.phi() << " , "
425  << " NOT associated to any reco::Track" << "\n";
426  }
427 
428  for (unsigned int f=0; f<etaintervals[w].size()-1; f++){
429  if (getEta(momentumTP.eta())>etaintervals[w][f]&&
430  getEta(momentumTP.eta())<etaintervals[w][f+1]) {
431  totSIMeta[w][f]++;
432  if (TP_is_matched) {
433  totASSeta[w][f]++;
434 
435  if (MABH) {
436  if (Quality075) {
437  totASSeta_Quality075[w][f]++;
438  totASSeta_Quality05[w][f]++;
439  }
440  else if (Quality05) {
441  totASSeta_Quality05[w][f]++;
442  }
443  }
444  }
445  }
446  } // END for (unsigned int f=0; f<etaintervals[w].size()-1; f++){
447 
448  for (unsigned int f=0; f<phiintervals[w].size()-1; f++){
449  if (momentumTP.phi() > phiintervals[w][f]&&
450  momentumTP.phi() <phiintervals[w][f+1]) {
451  totSIM_phi[w][f]++;
452  if (TP_is_matched) {
453  totASS_phi[w][f]++;
454 
455  if (MABH) {
456  if (Quality075) {
457  totASS_phi_Quality075[w][f]++;
458  totASS_phi_Quality05[w][f]++;
459  }
460  else if (Quality05) {
461  totASS_phi_Quality05[w][f]++;
462  }
463  }
464  }
465  }
466  } // END for (unsigned int f=0; f<phiintervals[w].size()-1; f++){
467 
468  for (unsigned int f=0; f<pTintervals[w].size()-1; f++){
469  if (getPt(sqrt(momentumTP.perp2()))>pTintervals[w][f]&&
470  getPt(sqrt(momentumTP.perp2()))<pTintervals[w][f+1]) {
471  totSIMpT[w][f]++;
472  if (TP_is_matched) {
473  totASSpT[w][f]++;
474 
475  if (MABH) {
476  if (Quality075) {
477  totASSpT_Quality075[w][f]++;
478  totASSpT_Quality05[w][f]++;
479  }
480  else if (Quality05) {
481  totASSpT_Quality05[w][f]++;
482  }
483  }
484  }
485  }
486  } // END for (unsigned int f=0; f<pTintervals[w].size()-1; f++){
487 
488  for (unsigned int f=0; f<dxyintervals[w].size()-1; f++){
489  if (dxySim>dxyintervals[w][f]&&
490  dxySim<dxyintervals[w][f+1]) {
491  totSIM_dxy[w][f]++;
492  if (TP_is_matched) {
493  totASS_dxy[w][f]++;
494  }
495  }
496  } // END for (unsigned int f=0; f<dxyintervals[w].size()-1; f++){
497 
498  for (unsigned int f=0; f<dzintervals[w].size()-1; f++){
499  if (dzSim>dzintervals[w][f]&&
500  dzSim<dzintervals[w][f+1]) {
501  totSIM_dz[w][f]++;
502  if (TP_is_matched) {
503  totASS_dz[w][f]++;
504  }
505  }
506  } // END for (unsigned int f=0; f<dzintervals[w].size()-1; f++){
507 
508  for (unsigned int f=0; f<vertposintervals[w].size()-1; f++){
509  if (sqrt(vertexTP.perp2())>vertposintervals[w][f]&&
510  sqrt(vertexTP.perp2())<vertposintervals[w][f+1]) {
511  totSIM_vertpos[w][f]++;
512  if (TP_is_matched) {
513  totASS_vertpos[w][f]++;
514  }
515  }
516  } // END for (unsigned int f=0; f<vertposintervals[w].size()-1; f++){
517 
518  for (unsigned int f=0; f<zposintervals[w].size()-1; f++){
519  if (vertexTP.z()>zposintervals[w][f]&&
520  vertexTP.z()<zposintervals[w][f+1]) {
521  totSIM_zpos[w][f]++;
522  if (TP_is_matched) {
523  totASS_zpos[w][f]++;
524  }
525  }
526  } // END for (unsigned int f=0; f<zposintervals[w].size()-1; f++){
527 
528  int nSimHits = 0;
529  if (usetracker && usemuon) {
530  nSimHits= tpr.get()->numberOfHits();
531  }
532  else if (!usetracker && usemuon) {
533  nSimHits= tpr.get()->numberOfHits() - tpr.get()->numberOfTrackerHits();
534  }
535  else if (usetracker && !usemuon) {
536  nSimHits=tpr.get()->numberOfTrackerHits();
537  }
538 
539  int tmp = std::min(nSimHits,int(maxHit-1));
540  edm::LogVerbatim("MuonTrackValidator") << "\t N simhits = "<< nSimHits<<"\n";
541 
542  totSIM_hit[w][tmp]++;
543  if (TP_is_matched) totASS_hit[w][tmp]++;
544 
545  if (TP_is_matched)
546  {
547  RefToBase<Track> assoctrack = rt.begin()->first;
548  nrecHit_vs_nsimHit_sim2rec[w]->Fill( assoctrack->numberOfValidHits(),nSimHits);
549  }
550  } // End for (TrackingParticleCollection::size_type i=0; i<tPCeff.size(); i++){
551  if (st!=0) h_tracksSIM[w]->Fill(st);
552 
553  //
554  //fill reconstructed track histograms
555  //
556  edm::LogVerbatim("MuonTrackValidator") << "\n# of reco::Tracks with "
557  << label[www].process()<<":"
558  << label[www].label()<<":"
559  << label[www].instance()
560  << ": " << trackCollectionSize << "\n";
561  int at = 0;
562  int rT = 0;
563  for(edm::View<Track>::size_type i=0; i<trackCollectionSize; ++i){
564  bool Track_is_matched = false;
565  RefToBase<Track> track(trackCollection, i);
566  rT++;
567 
568  std::vector<std::pair<TrackingParticleRef, double> > tp;
569  TrackingParticleRef tpr;
570 
571  // new logic (bidirectional)
572  if (BiDirectional_RecoToSim_association) {
573  edm::LogVerbatim("MuonTrackValidator")<<"----------------------------------------Track #"<< track.key();
574 
575  if(recSimColl.find(track) != recSimColl.end()) {
576  tp = recSimColl[track];
577  if (tp.size() != 0) {
578  tpr = tp.begin()->first;
579  // RtS and StR must associate the same pair !
580  if(simRecColl.find(tpr) != simRecColl.end()) {
581  std::vector<std::pair<RefToBase<Track>, double> > track_checkback = simRecColl[tpr];
582  RefToBase<Track> assoc_track_checkback;
583  assoc_track_checkback = track_checkback.begin()->first;
584 
585  if ( assoc_track_checkback.key() == track.key() ) {
586  edm::LogVerbatim("MuonTrackValidator")<<"------------------associated TrackingParticle #"<<tpr.key();
587  Track_is_matched = true;
588  at++;
589  double Purity = tp.begin()->second;
590  double Quality = track_checkback.begin()->second;
591  edm::LogVerbatim("MuonTrackValidator") << "reco::Track #" << track.key() << " with pt=" << track->pt()
592  << " associated with quality:" << Purity <<"\n";
593  if (MABH) h_PurityVsQuality[w]->Fill(Quality,Purity);
594  }
595  }
596  }
597  }
598 
599  if (!Track_is_matched)
600  edm::LogVerbatim("MuonTrackValidator")
601  << "reco::Track #" << track.key() << " with pt=" << track->pt() << " NOT associated to any TrackingParticle" << "\n";
602  }
603  // old logic (bugged for collision scenario, still valid for cosmics 2 legs reco)
604  else {
605  if(recSimColl.find(track) != recSimColl.end()){
606  tp = recSimColl[track];
607  if (tp.size()!=0) {
608  Track_is_matched = true;
609  tpr = tp.begin()->first;
610  at++;
611  edm::LogVerbatim("MuonTrackValidator") << "reco::Track #" << track.key() << " with pt=" << track->pt()
612  << " associated with quality:" << tp.begin()->second <<"\n";
613  }
614  } else {
615  edm::LogVerbatim("MuonTrackValidator") << "reco::Track #" << track.key() << " with pt=" << track->pt()
616  << " NOT associated to any TrackingParticle" << "\n";
617  }
618  }
619 
620  //Compute fake rate vs eta
621  for (unsigned int f=0; f<etaintervals[w].size()-1; f++){
622  if (getEta(track->momentum().eta())>etaintervals[w][f]&&
623  getEta(track->momentum().eta())<etaintervals[w][f+1]) {
624  totRECeta[w][f]++;
625  if (Track_is_matched) {
626  totASS2eta[w][f]++;
627  }
628  }
629  } // End for (unsigned int f=0; f<etaintervals[w].size()-1; f++){
630 
631  for (unsigned int f=0; f<phiintervals[w].size()-1; f++){
632  if (track->momentum().phi()>phiintervals[w][f]&&
633  track->momentum().phi()<phiintervals[w][f+1]) {
634  totREC_phi[w][f]++;
635  if (Track_is_matched) {
636  totASS2_phi[w][f]++;
637  }
638  }
639  } // End for (unsigned int f=0; f<phiintervals[w].size()-1; f++){
640 
641  for (unsigned int f=0; f<pTintervals[w].size()-1; f++){
642  if (getPt(sqrt(track->momentum().perp2()))>pTintervals[w][f]&&
643  getPt(sqrt(track->momentum().perp2()))<pTintervals[w][f+1]) {
644  totRECpT[w][f]++;
645  if (Track_is_matched) {
646  totASS2pT[w][f]++;
647  }
648  }
649  } // End for (unsigned int f=0; f<pTintervals[w].size()-1; f++){
650 
651  for (unsigned int f=0; f<dxyintervals[w].size()-1; f++){
652  if (track->dxy(bs.position())>dxyintervals[w][f]&&
653  track->dxy(bs.position())<dxyintervals[w][f+1]) {
654  totREC_dxy[w][f]++;
655  if (Track_is_matched) {
656  totASS2_dxy[w][f]++;
657  }
658  }
659  } // End for (unsigned int f=0; f<dxyintervals[w].size()-1; f++){
660 
661  for (unsigned int f=0; f<dzintervals[w].size()-1; f++){
662  if (track->dz(bs.position())>dzintervals[w][f]&&
663  track->dz(bs.position())<dzintervals[w][f+1]) {
664  totREC_dz[w][f]++;
665  if (Track_is_matched) {
666  totASS2_dz[w][f]++;
667  }
668  }
669  } // End for (unsigned int f=0; f<dzintervals[w].size()-1; f++){
670 
671  int tmp = std::min((int)track->found(),int(maxHit-1));
672  totREC_hit[w][tmp]++;
673  if (Track_is_matched) totASS2_hit[w][tmp]++;
674 
675  edm::LogVerbatim("MuonTrackValidator") << "\t N valid rechits = "<< (int)track->found() <<"\n";
676 
677  // Fill other histos
678  TrackingParticle* tpp = const_cast<TrackingParticle*>(tpr.get());
679  // TrackingParticle parameters at point of closest approach to the beamline
680  TrackingParticle::Vector momentumTP;
681  TrackingParticle::Point vertexTP;
682 
683  if (parametersDefiner=="LhcParametersDefinerForTP") {
684  // following reco plots are made only from tracks associated to selected signal TPs
685  if (! (Track_is_matched && tpSelector(*tpp)) ) continue;
686  else {
687  momentumTP = Lhc_parametersDefinerTP->momentum(event,setup,tpr) ;
688  vertexTP = Lhc_parametersDefinerTP->vertex(event,setup,tpr);
689  }
690  }
691  else if (parametersDefiner=="CosmicParametersDefinerForTP") {
692  // following reco plots are made only from tracks associated to selected signal TPs
693  if (! (Track_is_matched && cosmictpSelector(tpr,&bs,event,setup)) ) continue;
694  else {
695  momentumTP = Cosmic_parametersDefinerTP->momentum(event,setup,tpr) ;
696  vertexTP = Cosmic_parametersDefinerTP->vertex(event,setup,tpr);
697  }
698  }
699 
700  if (associators[ww]=="trackAssociatorByChi2"){
701  //association chi2
702  double assocChi2 = -tp.begin()->second;//in association map is stored -chi2
703  h_assochi2[www]->Fill(assocChi2);
704  h_assochi2_prob[www]->Fill(TMath::Prob((assocChi2)*5,5));
705  }
706  else if (associators[ww]=="trackAssociatorByHits"){
707  double fraction = tp.begin()->second;
708  h_assocFraction[www]->Fill(fraction);
709  h_assocSharedHit[www]->Fill(fraction*track->numberOfValidHits());
710  }
711 
712  //nchi2 and hits global distributions
713  h_nchi2[w]->Fill(track->normalizedChi2());
714  h_nchi2_prob[w]->Fill(TMath::Prob(track->chi2(),(int)track->ndof()));
715  h_hits[w]->Fill(track->numberOfValidHits());
716  h_losthits[w]->Fill(track->numberOfLostHits());
717  chi2_vs_nhits[w]->Fill(track->numberOfValidHits(),track->normalizedChi2());
718  h_charge[w]->Fill( track->charge() );
719 
720  double ptSim = sqrt(momentumTP.perp2());
721  double qoverpSim = tpr->charge()/sqrt(momentumTP.x()*momentumTP.x()+momentumTP.y()*momentumTP.y()+momentumTP.z()*momentumTP.z());
722  double thetaSim = momentumTP.theta();
723  double lambdaSim = M_PI/2-momentumTP.theta();
724  double phiSim = momentumTP.phi();
725  double dxySim = (-vertexTP.x()*sin(momentumTP.phi())+vertexTP.y()*cos(momentumTP.phi()));
726  double dzSim = vertexTP.z() - (vertexTP.x()*momentumTP.x()+vertexTP.y()*momentumTP.y()) /
727  sqrt(momentumTP.perp2()) * momentumTP.z()/sqrt(momentumTP.perp2());
728 
729  // removed unused variable, left this in case it has side effects
730  track->parameters();
731 
732  double qoverpRec(0);
733  double qoverpErrorRec(0);
734  double ptRec(0);
735  double ptErrorRec(0);
736  double lambdaRec(0);
737  double lambdaErrorRec(0);
738  double phiRec(0);
739  double phiErrorRec(0);
740 
741  //loop to decide whether to take gsfTrack (utilisation of mode-function) or common track
742  const GsfTrack* gsfTrack(0);
743  if(useGsf){
744  gsfTrack = dynamic_cast<const GsfTrack*>(&(*track));
745  if (gsfTrack==0) edm::LogInfo("MuonTrackValidator") << "Trying to access mode for a non-GsfTrack";
746  }
747 
748  if (gsfTrack) {
749  // get values from mode
750  getRecoMomentum(*gsfTrack, ptRec, ptErrorRec, qoverpRec, qoverpErrorRec,
751  lambdaRec,lambdaErrorRec, phiRec, phiErrorRec);
752  }
753 
754  else {
755  // get values from track (without mode)
756  getRecoMomentum(*track, ptRec, ptErrorRec, qoverpRec, qoverpErrorRec,
757  lambdaRec,lambdaErrorRec, phiRec, phiErrorRec);
758  }
759 
760  double thetaRec = track->theta();
761  double ptError = ptErrorRec;
762  double ptres = ptRec - ptSim;
763  double etares = track->eta()-momentumTP.Eta();
764  double dxyRec = track->dxy(bs.position());
765  double dzRec = track->dz(bs.position());
766  // eta residue; pt, k, theta, phi, dxy, dz pulls
767  double qoverpPull=(qoverpRec-qoverpSim)/qoverpErrorRec;
768  double thetaPull=(lambdaRec-lambdaSim)/lambdaErrorRec;
769  double phiDiff = phiRec - phiSim;
770  if (abs(phiDiff) > M_PI) {
771  if (phiDiff >0.) phiDiff = phiDiff - 2.*M_PI;
772  else phiDiff = phiDiff + 2.*M_PI;
773  }
774  double phiPull=phiDiff/phiErrorRec;
775  double dxyPull=(dxyRec-dxySim)/track->dxyError();
776  double dzPull=(dzRec-dzSim)/track->dzError();
777 
778  double contrib_Qoverp = ((qoverpRec-qoverpSim)/qoverpErrorRec)*((qoverpRec-qoverpSim)/qoverpErrorRec)/5;
779  double contrib_dxy = ((dxyRec-dxySim)/track->dxyError())*((dxyRec-dxySim)/track->dxyError())/5;
780  double contrib_dz = ((dzRec-dzSim)/track->dzError())*((dzRec-dzSim)/track->dzError())/5;
781  double contrib_theta = ((lambdaRec-lambdaSim)/lambdaErrorRec)*((lambdaRec-lambdaSim)/lambdaErrorRec)/5;
782  double contrib_phi = (phiDiff/phiErrorRec)*(phiDiff/phiErrorRec)/5;
783 
784  edm::LogVerbatim("MuonTrackValidator") << "assocChi2=" << tp.begin()->second << "\n"
785  << "" << "\n"
786  << "ptREC=" << ptRec << "\n"
787  << "etaREC=" << track->eta() << "\n"
788  << "qoverpREC=" << qoverpRec << "\n"
789  << "dxyREC=" << dxyRec << "\n"
790  << "dzREC=" << dzRec << "\n"
791  << "thetaREC=" << track->theta() << "\n"
792  << "phiREC=" << phiRec << "\n"
793  << "" << "\n"
794  << "qoverpError()=" << qoverpErrorRec << "\n"
795  << "dxyError()=" << track->dxyError() << "\n"
796  << "dzError()=" << track->dzError() << "\n"
797  << "thetaError()=" << lambdaErrorRec << "\n"
798  << "phiError()=" << phiErrorRec << "\n"
799  << "" << "\n"
800  << "ptSIM=" << ptSim << "\n"
801  << "etaSIM=" << momentumTP.Eta() << "\n"
802  << "qoverpSIM=" << qoverpSim << "\n"
803  << "dxySIM=" << dxySim << "\n"
804  << "dzSIM=" << dzSim << "\n"
805  << "thetaSIM=" << M_PI/2-lambdaSim << "\n"
806  << "phiSIM=" << phiSim << "\n"
807  << "" << "\n"
808  << "contrib_Qoverp=" << contrib_Qoverp << "\n"
809  << "contrib_dxy=" << contrib_dxy << "\n"
810  << "contrib_dz=" << contrib_dz << "\n"
811  << "contrib_theta=" << contrib_theta << "\n"
812  << "contrib_phi=" << contrib_phi << "\n"
813  << "" << "\n"
814  <<"chi2PULL="<<contrib_Qoverp+contrib_dxy+contrib_dz+contrib_theta+contrib_phi<<"\n";
815 
816  h_pullQoverp[w]->Fill(qoverpPull);
817  h_pullTheta[w]->Fill(thetaPull);
818  h_pullPhi[w]->Fill(phiPull);
819  h_pullDxy[w]->Fill(dxyPull);
820  h_pullDz[w]->Fill(dzPull);
821 
822  h_pt[w]->Fill(ptres/ptError);
823  h_eta[w]->Fill(etares);
824  etares_vs_eta[w]->Fill(getEta(track->eta()),etares);
825 
826  //chi2 and #hit vs eta: fill 2D histos
827  chi2_vs_eta[w]->Fill(getEta(track->eta()),track->normalizedChi2());
828  nhits_vs_eta[w]->Fill(getEta(track->eta()),track->numberOfValidHits());
829  nDThits_vs_eta[w]->Fill(getEta(track->eta()),track->hitPattern().numberOfValidMuonDTHits());
830  nCSChits_vs_eta[w]->Fill(getEta(track->eta()),track->hitPattern().numberOfValidMuonCSCHits());
831  nRPChits_vs_eta[w]->Fill(getEta(track->eta()),track->hitPattern().numberOfValidMuonRPCHits());
832  if(useGEMs_) nGEMhits_vs_eta[w]->Fill(getEta(track->eta()),track->hitPattern().numberOfValidMuonGEMHits());
833  if(useME0_) nME0hits_vs_eta[w]->Fill(getEta(track->eta()),track->hitPattern().numberOfValidMuonME0Hits());
834  nlosthits_vs_eta[w]->Fill(getEta(track->eta()),track->numberOfLostHits());
835 
836  //resolution of track params: fill 2D histos
837  dxyres_vs_eta[w]->Fill(getEta(track->eta()),dxyRec-dxySim);
838  ptres_vs_eta[w]->Fill(getEta(track->eta()),(ptRec-ptSim)/ptRec);
839  dzres_vs_eta[w]->Fill(getEta(track->eta()),dzRec-dzSim);
840  phires_vs_eta[w]->Fill(getEta(track->eta()),phiDiff);
841  cotThetares_vs_eta[w]->Fill(getEta(track->eta()), cos(thetaRec)/sin(thetaRec) - cos(thetaSim)/sin(thetaSim));
842 
843  //same as before but vs pT
844  dxyres_vs_pt[w]->Fill(getPt(ptRec),dxyRec-dxySim);
845  ptres_vs_pt[w]->Fill(getPt(ptRec),(ptRec-ptSim)/ptRec);
846  dzres_vs_pt[w]->Fill(getPt(ptRec),dzRec-dzSim);
847  phires_vs_pt[w]->Fill(getPt(ptRec),phiDiff);
848  cotThetares_vs_pt[w]->Fill(getPt(ptRec), cos(thetaRec)/sin(thetaRec) - cos(thetaSim)/sin(thetaSim));
849 
850  //pulls of track params vs eta: fill 2D histos
851  dxypull_vs_eta[w]->Fill(getEta(track->eta()),dxyPull);
852  ptpull_vs_eta[w]->Fill(getEta(track->eta()),ptres/ptError);
853  dzpull_vs_eta[w]->Fill(getEta(track->eta()),dzPull);
854  phipull_vs_eta[w]->Fill(getEta(track->eta()),phiPull);
855  thetapull_vs_eta[w]->Fill(getEta(track->eta()),thetaPull);
856 
857  //plots vs phi
858  nhits_vs_phi[w]->Fill(phiRec,track->numberOfValidHits());
859  chi2_vs_phi[w]->Fill(phiRec,track->normalizedChi2());
860  ptmean_vs_eta_phi[w]->Fill(phiRec,getEta(track->eta()),ptRec);
861  phimean_vs_eta_phi[w]->Fill(phiRec,getEta(track->eta()),phiRec);
862  ptres_vs_phi[w]->Fill(phiRec,(ptRec-ptSim)/ptRec);
863  phires_vs_phi[w]->Fill(phiRec,phiDiff);
864  ptpull_vs_phi[w]->Fill(phiRec,ptres/ptError);
865  phipull_vs_phi[w]->Fill(phiRec,phiPull);
866  thetapull_vs_phi[w]->Fill(phiRec,thetaPull);
867 
868  int nSimHits = 0;
869  if (usetracker && usemuon) {
870  nSimHits= tpr.get()->numberOfHits();
871  }
872  else if (!usetracker && usemuon) {
873  nSimHits= tpr.get()->numberOfHits() - tpr.get()->numberOfTrackerHits();
874  }
875  else if (usetracker && !usemuon) {
876  nSimHits=tpr.get()->numberOfTrackerHits();
877  }
878 
879  nrecHit_vs_nsimHit_rec2sim[w]->Fill(track->numberOfValidHits(), nSimHits);
880 
881  } // End of for(View<Track>::size_type i=0; i<trackCollectionSize; ++i){
882 
883  if (at!=0) h_tracks[w]->Fill(at);
884  h_fakes[w]->Fill(rT-at);
885  edm::LogVerbatim("MuonTrackValidator") << "Total Simulated: " << st << "\n"
886  << "Total Associated (simToReco): " << ats << "\n"
887  << "Total Reconstructed: " << rT << "\n"
888  << "Total Associated (recoToSim): " << at << "\n"
889  << "Total Fakes: " << rT-at << "\n";
890  nrec_vs_nsim[w]->Fill(rT,st);
891  w++;
892  } // End of for (unsigned int www=0;www<label.size();www++){
893  } //END of for (unsigned int ww=0;ww<associators.size();ww++){
894 }
895 
896 void MuonTrackValidator::endRun(Run const&, EventSetup const&) {
897 
898  int w=0;
899  for (unsigned int ww=0;ww<associators.size();ww++){
900  for (unsigned int www=0;www<label.size();www++){
901 
902  //chi2 and #hit vs eta: get mean from 2D histos
903  doProfileX(chi2_vs_eta[w],h_chi2meanh[w]);
904  doProfileX(nhits_vs_eta[w],h_hits_eta[w]);
905  doProfileX(nDThits_vs_eta[w],h_DThits_eta[w]);
906  doProfileX(nCSChits_vs_eta[w],h_CSChits_eta[w]);
907  doProfileX(nRPChits_vs_eta[w],h_RPChits_eta[w]);
908  if (useGEMs_) doProfileX(nGEMhits_vs_eta[w],h_GEMhits_eta[w]);
909  if (useME0_) doProfileX(nME0hits_vs_eta[w],h_ME0hits_eta[w]);
910 
911  doProfileX(nlosthits_vs_eta[w],h_losthits_eta[w]);
912  //vs phi
913  doProfileX(chi2_vs_nhits[w],h_chi2meanhitsh[w]);
914  doProfileX(chi2_vs_phi[w],h_chi2mean_vs_phi[w]);
915  doProfileX(nhits_vs_phi[w],h_hits_phi[w]);
916 
917  fillPlotFromVector(h_recoeta[w],totRECeta[w]);
918  fillPlotFromVector(h_simuleta[w],totSIMeta[w]);
919  fillPlotFromVector(h_assoceta[w],totASSeta[w]);
920  fillPlotFromVector(h_assoc2eta[w],totASS2eta[w]);
921 
922  fillPlotFromVector(h_recopT[w],totRECpT[w]);
923  fillPlotFromVector(h_simulpT[w],totSIMpT[w]);
924  fillPlotFromVector(h_assocpT[w],totASSpT[w]);
925  fillPlotFromVector(h_assoc2pT[w],totASS2pT[w]);
926 
927  fillPlotFromVector(h_recohit[w],totREC_hit[w]);
928  fillPlotFromVector(h_simulhit[w],totSIM_hit[w]);
929  fillPlotFromVector(h_assochit[w],totASS_hit[w]);
930  fillPlotFromVector(h_assoc2hit[w],totASS2_hit[w]);
931 
932  fillPlotFromVector(h_recophi[w],totREC_phi[w]);
933  fillPlotFromVector(h_simulphi[w],totSIM_phi[w]);
934  fillPlotFromVector(h_assocphi[w],totASS_phi[w]);
935  fillPlotFromVector(h_assoc2phi[w],totASS2_phi[w]);
936 
937  fillPlotFromVector(h_recodxy[w],totREC_dxy[w]);
938  fillPlotFromVector(h_simuldxy[w],totSIM_dxy[w]);
939  fillPlotFromVector(h_assocdxy[w],totASS_dxy[w]);
940  fillPlotFromVector(h_assoc2dxy[w],totASS2_dxy[w]);
941 
942  fillPlotFromVector(h_recodz[w],totREC_dz[w]);
943  fillPlotFromVector(h_simuldz[w],totSIM_dz[w]);
944  fillPlotFromVector(h_assocdz[w],totASS_dz[w]);
945  fillPlotFromVector(h_assoc2dz[w],totASS2_dz[w]);
946 
947  fillPlotFromVector(h_simulvertpos[w],totSIM_vertpos[w]);
948  fillPlotFromVector(h_assocvertpos[w],totASS_vertpos[w]);
949 
950  fillPlotFromVector(h_simulzpos[w],totSIM_zpos[w]);
951  fillPlotFromVector(h_assoczpos[w],totASS_zpos[w]);
952 
953  if (MABH) {
954  fillPlotFromVector(h_assoceta_Quality05[w] ,totASSeta_Quality05[w]);
955  fillPlotFromVector(h_assoceta_Quality075[w],totASSeta_Quality075[w]);
956  fillPlotFromVector(h_assocpT_Quality05[w] ,totASSpT_Quality05[w]);
957  fillPlotFromVector(h_assocpT_Quality075[w],totASSpT_Quality075[w]);
958  fillPlotFromVector(h_assocphi_Quality05[w] ,totASS_phi_Quality05[w]);
959  fillPlotFromVector(h_assocphi_Quality075[w],totASS_phi_Quality075[w]);
960  }
961 
962  w++;
963  }
964  }
965 
966  if ( out.size() != 0 && dbe_ ) dbe_->save(out);
967 }
968 
969 void MuonTrackValidator::getRecoMomentum (const reco::Track& track,
970  double& pt, double& ptError, double& qoverp, double& qoverpError,
971  double& lambda,double& lambdaError, double& phi, double& phiError ) const {
972  pt = track.pt();
973  ptError = track.ptError();
974  qoverp = track.qoverp();
975  qoverpError = track.qoverpError();
976  lambda = track.lambda();
977  lambdaError = track.lambdaError();
978  phi = track.phi();
979  phiError = track.phiError();
980 }
981 
982 void MuonTrackValidator::getRecoMomentum (const reco::GsfTrack& gsfTrack,
983  double& pt, double& ptError, double& qoverp, double& qoverpError,
984  double& lambda,double& lambdaError, double& phi, double& phiError ) const {
985  pt = gsfTrack.ptMode();
986  ptError = gsfTrack.ptModeError();
987  qoverp = gsfTrack.qoverpMode();
988  qoverpError = gsfTrack.qoverpModeError();
989  lambda = gsfTrack.lambdaMode();
990  lambdaError = gsfTrack.lambdaModeError();
991  phi = gsfTrack.phiMode();
992  phiError = gsfTrack.phiModeError();
993 }
reco::TrackBase::qoverp
double qoverp() const
q / p
Definition: TrackBase.h:568
reco::GsfTrack::phiModeError
double phiModeError() const
error on phi from mode
Definition: GsfTrack.h:94
edm::View::size_type
unsigned int size_type
Definition: View.h:85
i
int i
Definition: DBlmapReader.cc:9
w
const double w
Definition: UKUtility.cc:23
TrackingParticleCollection
std::vector< TrackingParticle > TrackingParticleCollection
Definition: HiEvtPlaneFlatProducer.cc:111
reco::GsfTrack::lambdaMode
double lambdaMode() const
Lambda angle from mode.
Definition: GsfTrack.h:45
TrackingParticle::momentum
Vector momentum() const
spatial momentum vector
Definition: TrackingParticle.h:109
edm::AssociationMap::end
const_iterator end() const
last iterator over the map (read only)
Definition: AssociationMap.h:160
reco::TrackBase::normalizedChi2
double normalizedChi2() const
chi-squared divided by n.d.o.f. (or chi-squared * 1e6 if n.d.o.f. is zero)
Definition: TrackBase.h:556
DQMStore::IBooker::bookProfile
MonitorElement * bookProfile(Args &&...args)
Definition: DQMStore.h:157
gather_cfg.maxDxy
tuple maxDxy
Definition: gather_cfg.py:52
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:457
DQMStore::IBooker::cd
void cd(void)
Definition: DQMStore.cc:269
reco::TrackBase::theta
double theta() const
polar angle
Definition: TrackBase.h:574
reco::TrackBase::dxyError
double dxyError() const
error on dxy
Definition: TrackBase.h:791
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
edm::AssociationMap::find
const_iterator find(const key_type &k) const
find element with specified reference key
Definition: AssociationMap.h:162
getEta
static float getEta(float eta, int bits=5)
Definition: PtAssignment.cc:386
reco::TrackBase::phi
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:640
reco::TrackBase::numberOfLostHits
unsigned short numberOfLostHits() const
number of cases where track crossed a layer without getting a hit.
Definition: TrackBase.h:821
HLT_FULL_cff.fraction
tuple fraction
Definition: HLT_FULL_cff.py:43957
ecalcalib_dqm_sourceclient-live_cfg.algo
tuple algo
Definition: ecalcalib_dqm_sourceclient-live_cfg.py:65
edm::Ref::key
key_type key() const
Accessor for product key.
Definition: Ref.h:264
reco::TrackBase::momentum
const Vector & momentum() const
track momentum vector
Definition: TrackBase.h:670
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:107
MuonTrackValidator::getRecoMomentum
void getRecoMomentum(const reco::Track &track, double &pt, double &ptError, double &qoverp, double &qoverpError, double &lambda, double &lambdaError, double &phi, double &phiError) const
retrieval of reconstructed momentum components from reco::Track (== mean values for GSF) ...
Definition: MuonTrackValidator.cc:969
trigger::size_type
uint16_t size_type
Definition: TriggerTypeDefs.h:19
reco::HitPattern::numberOfValidMuonCSCHits
int numberOfValidMuonCSCHits() const
Definition: HitPattern.h:883
reco::GsfTrack::qoverpMode
double qoverpMode() const
q/p from mode
Definition: GsfTrack.h:41
reco::HitPattern::numberOfValidMuonRPCHits
int numberOfValidMuonRPCHits() const
Definition: HitPattern.h:888
reco::GsfTrack::ptModeError
double ptModeError() const
error on Pt (set to 1000 TeV if charge==0 for safety) from mode
Definition: GsfTrack.h:81
TrackingParticle::Point
math::XYZPointD Point
point in the space
Definition: TrackingParticle.h:36
HLT_FULL_cff.minHit
tuple minHit
Definition: HLT_FULL_cff.py:48660
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:646
TrackCollections2monitor_cff.quality
tuple quality
[pTError/pT]*max(1,normChi2) &lt;= ptErrorCut
Definition: TrackCollections2monitor_cff.py:205
reco::TrackBase::chi2
double chi2() const
chi-squared of the fit
Definition: TrackBase.h:544
DQMStore::IBooker::bookProfile2D
MonitorElement * bookProfile2D(Args &&...args)
Definition: DQMStore.h:163
TrackerOfflineValidation_Dqm_cff.dirName
tuple dirName
Definition: TrackerOfflineValidation_Dqm_cff.py:56
edm::RefToBase::key
size_t key() const
Definition: RefToBase.h:250
edm::AssociationMap::post_insert
void post_insert()
post insert action
Definition: AssociationMap.h:219
reco::TrackBase::ndof
double ndof() const
number of degrees of freedom of the fit
Definition: TrackBase.h:550
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
reco::TrackBase::pt
double pt() const
track transverse momentum
Definition: TrackBase.h:616
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
reco::TrackBase::ptError
double ptError() const
error on Pt (set to 1000 TeV if charge==0 for safety)
Definition: TrackBase.h:758
reco::TrackBase::phiError
double phiError() const
error on phi
Definition: TrackBase.h:785
DQMStore::IBooker::book1D
MonitorElement * book1D(Args &&...args)
Definition: DQMStore.h:115
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
j
int j
Definition: DBlmapReader.cc:9
reco::TrackBase::lambda
double lambda() const
Lambda angle.
Definition: TrackBase.h:580
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
fileinputsource_cfi.out
tuple out
Definition: fileinputsource_cfi.py:62
MuonTrackValidator::analyze
void analyze(const edm::Event &, const edm::EventSetup &)
Method called before the event loop.
Definition: MuonTrackValidator.cc:232
reco::TrackBase::numberOfValidHits
unsigned short numberOfValidHits() const
number of valid hits found
Definition: TrackBase.h:815
edm::Ref::get
T const * get() const
Returns C++ pointer to the item.
Definition: Ref.h:244
min
T min(T a, T b)
Definition: MathUtil.h:58
event
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
LogTrace
#define LogTrace(id)
Definition: PFTrackAlgoTools.cc:502
dbe_
DQMStore * dbe_
Definition: PFJetBenchmarkAnalyzer.cc:79
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:51
nint
int nint(float a)
Return the nearest integer - analogous to the FORTRAN intrinsic NINT.
Definition: nint.h:8
reco::TrackBase::qoverpError
double qoverpError() const
error on signed transverse curvature
Definition: TrackBase.h:752
reco::TrackBase::dz
double dz() const
dz parameter (= dsz/cos(lambda)). This is the track z0 w.r.t (0,0,0) only if the refPoint is close to...
Definition: TrackBase.h:604
reco::TrackBase::dzError
double dzError() const
error on dz
Definition: TrackBase.h:809
HLT_FULL_cff.trackCollection
tuple trackCollection
Definition: HLT_FULL_cff.py:62148
DQMStore::IBooker::setCurrentFolder
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:277
DQMStore::IBooker::goUp
void goUp(void)
Definition: DQMStore.cc:281
edm::Handle::product
T const * product() const
Definition: Handle.h:81
DQMStore::IBooker::book2D
MonitorElement * book2D(Args &&...args)
Definition: DQMStore.h:133
reco::TrackBase::hitPattern
const HitPattern & hitPattern() const
Access the hit pattern, indicating in which Tracker layers the track has hits.
Definition: TrackBase.h:445
reco::GsfTrack::qoverpModeError
double qoverpModeError() const
error on signed transverse curvature from mode
Definition: GsfTrack.h:79
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:56
edm::InputTag::label
std::string const & label() const
Definition: InputTag.h:36
TrackingParticle::vertex
Point vertex() const
Parent vertex position.
Definition: TrackingParticle.h:199
reco::TrackBase::parameters
ParameterVector parameters() const
Track parameters with one-to-one correspondence to the covariance matrix.
Definition: TrackBase.h:720
edm::InputTag::process
std::string const & process() const
Definition: InputTag.h:40
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
reco::TrackBase::lambdaError
double lambdaError() const
error on lambda
Definition: TrackBase.h:773
reco::HitPattern::numberOfValidMuonGEMHits
int numberOfValidMuonGEMHits() const
Definition: HitPattern.h:893
MuonTrackValidator::bookHistograms
void bookHistograms(DQMStore::IBooker &, edm::Run const &, edm::EventSetup const &)
Definition: MuonTrackValidator.cc:26
reco::HitPattern::numberOfValidMuonME0Hits
int numberOfValidMuonME0Hits() const
Definition: HitPattern.h:898
reco::Track::found
unsigned short found() const
Number of valid hits on track.
Definition: Track.h:194
reco::HitPattern::numberOfValidMuonDTHits
int numberOfValidMuonDTHits() const
Definition: HitPattern.h:878
reco::GsfTrack::phiMode
double phiMode() const
azimuthal angle of momentum vector from mode
Definition: GsfTrack.h:57
TrackingParticle
Monte Carlo truth information used for tracking validation.
Definition: TrackingParticle.h:29
reco::TrackBase::charge
int charge() const
track electric charge
Definition: TrackBase.h:562
reco::GsfTrack::lambdaModeError
double lambdaModeError() const
error on lambda from mode
Definition: GsfTrack.h:90
TrackingParticle::Vector
math::XYZVectorD Vector
point in the space
Definition: TrackingParticle.h:37
reco::TrackBase::dxy
double dxy() const
dxy parameter. (This is the transverse impact parameter w.r.t. to (0,0,0) ONLY if refPoint is close t...
Definition: TrackBase.h:586
reco::GsfTrack::ptMode
double ptMode() const
track transverse momentum from mode
Definition: GsfTrack.h:49
edm::InputTag::instance
std::string const & instance() const
Definition: InputTag.h:37
edm::Run
Definition: Run.h:42
MuonTrackValidator::endRun
void endRun(edm::Run const &, edm::EventSetup const &)
Method called at the end of the event loop.
Definition: MuonTrackValidator.cc:896
Generated for CMSSW Reference Manual by doxygen 1.8.5