CMS 3D CMS Logo

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