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