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