CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TkConvValidator.cc
Go to the documentation of this file.
1 #include <iostream>
2 //
5 //
7 
8 //
12 //
18 //
22 //
29 //
33 //
36 #include "CLHEP/Units/GlobalPhysicalConstants.h"
38 //
48 
49 //
73 
74 //
83 //
84 //
85 #include "TFile.h"
86 #include "TH1.h"
87 #include "TH2.h"
88 #include "TTree.h"
89 #include "TVector3.h"
90 #include "TProfile.h"
91 //
102 using namespace std;
103 
104 
106  {
107 
108  fName_ = pset.getUntrackedParameter<std::string>("Name");
109  verbosity_ = pset.getUntrackedParameter<int>("Verbosity");
110  parameters_ = pset;
111 
112  conversionCollectionProducer_ = pset.getParameter<std::string>("convProducer");
113  conversionCollection_ = pset.getParameter<std::string>("conversionCollection");
114  // conversionTrackProducer_ = pset.getParameter<std::string>("trackProducer");
115  minPhoEtCut_ = pset.getParameter<double>("minPhoEtCut");
116  mergedTracks_ = pset.getParameter<bool>("mergedTracks");
117  isRunCentrally_= pset.getParameter<bool>("isRunCentrally");
118 
119  }
120 
121 
122 
123 
124 
126 
127 
128 
129 
131 
132  nEvt_=0;
133  nEntry_=0;
134  nRecConv_=0;
135  nRecConvAss_=0;
136  nRecConvAssWithEcal_=0;
137 
138  nInvalidPCA_=0;
139 
140  dbe_ = 0;
142 
143 
144  double etMin = parameters_.getParameter<double>("etMin");
145  double etMax = parameters_.getParameter<double>("etMax");
146  int etBin = parameters_.getParameter<int>("etBin");
147 
148 
149  double resMin = parameters_.getParameter<double>("resMin");
150  double resMax = parameters_.getParameter<double>("resMax");
151  int resBin = parameters_.getParameter<int>("resBin");
152 
153  double etaMin = parameters_.getParameter<double>("etaMin");
154  double etaMax = parameters_.getParameter<double>("etaMax");
155  int etaBin = parameters_.getParameter<int>("etaBin");
156  int etaBin2 = parameters_.getParameter<int>("etaBin2");
157 
158 
159  double phiMin = parameters_.getParameter<double>("phiMin");
160  double phiMax = parameters_.getParameter<double>("phiMax");
161  int phiBin = parameters_.getParameter<int>("phiBin");
162 
163 
164  double rMin = parameters_.getParameter<double>("rMin");
165  double rMax = parameters_.getParameter<double>("rMax");
166  int rBin = parameters_.getParameter<int>("rBin");
167 
168  double zMin = parameters_.getParameter<double>("zMin");
169  double zMax = parameters_.getParameter<double>("zMax");
170  int zBin = parameters_.getParameter<int>("zBin");
171 
172  double dPhiTracksMin = parameters_.getParameter<double>("dPhiTracksMin");
173  double dPhiTracksMax = parameters_.getParameter<double>("dPhiTracksMax");
174  int dPhiTracksBin = parameters_.getParameter<int>("dPhiTracksBin");
175 
176  // double dEtaTracksMin = parameters_.getParameter<double>("dEtaTracksMin"); // unused
177  // double dEtaTracksMax = parameters_.getParameter<double>("dEtaTracksMax"); // unused
178  // int dEtaTracksBin = parameters_.getParameter<int>("dEtaTracksBin"); // unused
179 
180  double dCotTracksMin = parameters_.getParameter<double>("dCotTracksMin");
181  double dCotTracksMax = parameters_.getParameter<double>("dCotTracksMax");
182  int dCotTracksBin = parameters_.getParameter<int>("dCotTracksBin");
183 
184 
185  double chi2Min = parameters_.getParameter<double>("chi2Min");
186  double chi2Max = parameters_.getParameter<double>("chi2Max");
187 
188 
189  double rMinForXray = parameters_.getParameter<double>("rMinForXray");
190  double rMaxForXray = parameters_.getParameter<double>("rMaxForXray");
191  int rBinForXray = parameters_.getParameter<int>("rBinForXray");
192  double zMinForXray = parameters_.getParameter<double>("zMinForXray");
193  double zMaxForXray = parameters_.getParameter<double>("zMaxForXray");
194  int zBinForXray = parameters_.getParameter<int>("zBinForXray");
195  int zBin2ForXray = parameters_.getParameter<int>("zBin2ForXray");
196 
197  minPhoPtForEffic = parameters_.getParameter<double>("minPhoPtForEffic");
198  maxPhoEtaForEffic = parameters_.getParameter<double>("maxPhoEtaForEffic");
199  maxPhoZForEffic = parameters_.getParameter<double>("maxPhoZForEffic");
200  maxPhoRForEffic = parameters_.getParameter<double>("maxPhoRForEffic");
201  minPhoPtForPurity = parameters_.getParameter<double>("minPhoPtForPurity");
202  maxPhoEtaForPurity = parameters_.getParameter<double>("maxPhoEtaForPurity");
203  maxPhoZForPurity = parameters_.getParameter<double>("maxPhoZForPurity");
204  maxPhoRForPurity = parameters_.getParameter<double>("maxPhoRForPurity");
205 
206  if (dbe_) {
207 
209  // SC from reco photons
210 
211  dbe_->setCurrentFolder("EgammaV/ConversionValidator/SimulationInfo");
212  //
213  // simulation information about conversions
215  std::string histname = "nOfSimConversions";
216  h_nSimConv_[0] = dbe_->book1D(histname,"# of Sim conversions per event ",20,-0.5,19.5);
218  histname = "h_AllSimConvEta";
219  h_AllSimConv_[0] = dbe_->book1D(histname," All conversions: simulated #eta",etaBin2,etaMin,etaMax);
220  histname = "h_AllSimConvPhi";
221  h_AllSimConv_[1] = dbe_->book1D(histname," All conversions: simulated #phi",phiBin,phiMin,phiMax);
222  histname = "h_AllSimConvR";
223  h_AllSimConv_[2] = dbe_->book1D(histname," All conversions: simulated R",rBin,rMin,rMax);
224  histname = "h_AllSimConvZ";
225  h_AllSimConv_[3] = dbe_->book1D(histname," All conversions: simulated Z",zBin,zMin,zMax);
226  histname = "h_AllSimConvEt";
227  h_AllSimConv_[4] = dbe_->book1D(histname," All conversions: simulated Et",etBin,etMin,etMax);
228  //
229  histname = "nOfVisSimConversions";
230  h_nSimConv_[1] = dbe_->book1D(histname,"# of Sim conversions per event ",20,-0.5,19.5);
231  histname = "h_VisSimConvEta";
232  h_VisSimConv_[0] = dbe_->book1D(histname," All vis conversions: simulated #eta",etaBin2,etaMin, etaMax);
233  histname = "h_VisSimConvPhi";
234  h_VisSimConv_[1] = dbe_->book1D(histname," All vis conversions: simulated #phi",phiBin,phiMin, phiMax);
235  histname = "h_VisSimConvR";
236  h_VisSimConv_[2] = dbe_->book1D(histname," All vis conversions: simulated R",rBin,rMin,rMax);
237  histname = "h_VisSimConvZ";
238  h_VisSimConv_[3] = dbe_->book1D(histname," All vis conversions: simulated Z",zBin,zMin, zMax);
239  histname = "h_VisSimConvEt";
240  h_VisSimConv_[4] = dbe_->book1D(histname," All vis conversions: simulated Et",etBin,etMin, etMax);
241 
242  //
243  histname = "h_SimConvTwoMTracksEta";
244  h_SimConvTwoMTracks_[0] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks: simulated #eta",etaBin2,etaMin, etaMax);
245  histname = "h_SimConvTwoMTracksPhi";
246  h_SimConvTwoMTracks_[1] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks: simulated #phi",phiBin,phiMin, phiMax);
247  histname = "h_SimConvTwoMTracksR";
248  h_SimConvTwoMTracks_[2] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks: simulated R",rBin,rMin, rMax);
249  histname = "h_SimConvTwoMTracksZ";
250  h_SimConvTwoMTracks_[3] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks: simulated Z",zBin,zMin, zMax);
251  histname = "h_SimConvTwoMTracksEt";
252  h_SimConvTwoMTracks_[4] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks: simulated Et",etBin,etMin, etMax);
253  //
254  histname = "h_SimConvTwoTracksEta";
255  h_SimConvTwoTracks_[0] = dbe_->book1D(histname," All vis conversions with 2 reco tracks: simulated #eta",etaBin2,etaMin, etaMax);
256  histname = "h_SimConvTwoTracksPhi";
257  h_SimConvTwoTracks_[1] = dbe_->book1D(histname," All vis conversions with 2 reco tracks: simulated #phi",phiBin,phiMin, phiMax);
258  histname = "h_SimConvTwoTracksR";
259  h_SimConvTwoTracks_[2] = dbe_->book1D(histname," All vis conversions with 2 reco tracks: simulated R",rBin,rMin, rMax);
260  histname = "h_SimConvTwoTracksZ";
261  h_SimConvTwoTracks_[3] = dbe_->book1D(histname," All vis conversions with 2 reco tracks: simulated Z",zBin,zMin, zMax);
262  histname = "h_SimConvTwoTracksEt";
263  h_SimConvTwoTracks_[4] = dbe_->book1D(histname," All vis conversions with 2 reco tracks: simulated Et",etBin,etMin, etMax);
264  //
265  histname = "h_SimConvTwoMTracksEtaAndVtxPGT0";
266  h_SimConvTwoMTracksAndVtxPGT0_[0] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks + vertex: simulated #eta",etaBin2,etaMin, etaMax);
267  histname = "h_SimConvTwoMTracksPhiAndVtxPGT0";
268  h_SimConvTwoMTracksAndVtxPGT0_[1] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks + vertex: simulated #phi",phiBin,phiMin, phiMax);
269  histname = "h_SimConvTwoMTracksRAndVtxPGT0";
270  h_SimConvTwoMTracksAndVtxPGT0_[2] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks + vertex: simulated R",rBin,rMin, rMax);
271  histname = "h_SimConvTwoMTracksZAndVtxPGT0";
272  h_SimConvTwoMTracksAndVtxPGT0_[3] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks + vertex: simulated Z",zBin,zMin, zMax);
273  histname = "h_SimConvTwoMTracksEtAndVtxPGT0";
274  h_SimConvTwoMTracksAndVtxPGT0_[4] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks + vertex: simulated Et",etBin,etMin, etMax);
275 
276  //
277  histname = "h_SimConvTwoMTracksEtaAndVtxPGT0005";
278  h_SimConvTwoMTracksAndVtxPGT0005_[0] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks + vertex: simulated #eta",etaBin2,etaMin, etaMax);
279  histname = "h_SimConvTwoMTracksPhiAndVtxPGT0005";
280  h_SimConvTwoMTracksAndVtxPGT0005_[1] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks + vertex: simulated #phi",phiBin,phiMin, phiMax);
281  histname = "h_SimConvTwoMTracksRAndVtxPGT0005";
282  h_SimConvTwoMTracksAndVtxPGT0005_[2] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks + vertex: simulated R",rBin,rMin, rMax);
283  histname = "h_SimConvTwoMTracksZAndVtxPGT0005";
284  h_SimConvTwoMTracksAndVtxPGT0005_[3] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks + vertex: simulated Z",zBin,zMin, zMax);
285  histname = "h_SimConvTwoMTracksEtAndVtxPGT0005";
286  h_SimConvTwoMTracksAndVtxPGT0005_[4] = dbe_->book1D(histname," All vis conversions with 2 reco-matching tracks + vertex: simulated Et",etBin,etMin, etMax);
287 
288 
289  h_SimConvEtaPix_[0] = dbe_->book1D("simConvEtaPix"," sim converted Photon Eta: Pix ",etaBin,etaMin, etaMax) ;
290  h_simTkPt_ = dbe_->book1D("simTkPt","Sim conversion tracks pt ",etBin*3,0.,etMax);
291  h_simTkEta_ = dbe_->book1D("simTkEta","Sim conversion tracks eta ",etaBin,etaMin,etaMax);
292 
293  h_simConvVtxRvsZ_[0] = dbe_->book2D("simConvVtxRvsZAll"," Photon Sim conversion vtx position",zBinForXray, zMinForXray, zMaxForXray, rBinForXray, rMinForXray, rMaxForXray);
294  h_simConvVtxRvsZ_[1] = dbe_->book2D("simConvVtxRvsZBarrel"," Photon Sim conversion vtx position",zBinForXray, zMinForXray, zMaxForXray, rBinForXray, rMinForXray, rMaxForXray);
295  h_simConvVtxRvsZ_[2] = dbe_->book2D("simConvVtxRvsZEndcap"," Photon Sim conversion vtx position",zBin2ForXray, zMinForXray, zMaxForXray, rBinForXray, rMinForXray, rMaxForXray);
296  h_simConvVtxRvsZ_[3] = dbe_->book2D("simConvVtxRvsZBarrel2"," Photon Sim conversion vtx position when reco R<4cm",zBinForXray, zMinForXray, zMaxForXray, rBinForXray, rMinForXray, rMaxForXray);
297  h_simConvVtxYvsX_ = dbe_->book2D("simConvVtxYvsXTrkBarrel"," Photon Sim conversion vtx position, (x,y) eta<1 ",100, -80., 80., 100, -80., 80.);
298 
299  dbe_->setCurrentFolder("EgammaV/ConversionValidator/ConversionInfo");
300 
301  histname="nConv";
302  h_nConv_[0][0] = dbe_->book1D(histname+"All","Number Of Conversions per isolated candidates per events: All Ecal ",10,-0.5, 9.5);
303  h_nConv_[0][1] = dbe_->book1D(histname+"Barrel","Number Of Conversions per isolated candidates per events: Ecal Barrel ",10,-0.5, 9.5);
304  h_nConv_[0][2] = dbe_->book1D(histname+"Endcap","Number Of Conversions per isolated candidates per events: Ecal Endcap ",10,-0.5, 9.5);
305  h_nConv_[1][0] = dbe_->book1D(histname+"All_Ass","Number Of associated Conversions per isolated candidates per events: All Ecal ",10,-0.5, 9.5);
306 
307  h_convEta_[0][0] = dbe_->book1D("convEta"," converted Photon Eta ",etaBin,etaMin, etaMax) ;
308  h_convPhi_[0][0] = dbe_->book1D("convPhi"," converted Photon Phi ",phiBin,phiMin,phiMax) ;
309  h_convR_[0][0] = dbe_->book1D("convR"," converted photon R",rBin,rMin, rMax);
310  h_convZ_[0][0] = dbe_->book1D("convZ"," converted photon Z",zBin,zMin, zMax);
311  h_convPt_[0][0] = dbe_->book1D("convPt"," conversions Transverse Energy: all eta ", etBin,etMin, etMax);
312 
313  h_convEta_[1][0] = dbe_->book1D("convEtaAss"," Matched converted Photon Eta ",etaBin,etaMin, etaMax) ;
314  h_convPhi_[1][0] = dbe_->book1D("convPhiAss"," Matched converted Photon Phi ",phiBin,phiMin,phiMax) ;
315  h_convR_[1][0] = dbe_->book1D("convRAss"," Matched converted photon R",rBin,rMin, rMax);
316  h_convZ_[1][0] = dbe_->book1D("convZAss"," Matched converted photon Z",zBin,zMin, zMax);
317  h_convPt_[1][0] = dbe_->book1D("convPtAss","Matched conversions Transverse Energy: all eta ", etBin,etMin, etMax);
318 
319  h_convEta_[2][0] = dbe_->book1D("convEtaFake"," Fake converted Photon Eta ",etaBin,etaMin, etaMax) ;
320  h_convPhi_[2][0] = dbe_->book1D("convPhiFake"," Fake converted Photon Phi ",phiBin,phiMin,phiMax) ;
321  h_convR_[2][0] = dbe_->book1D("convRFake"," Fake converted photon R",rBin,rMin, rMax);
322  h_convZ_[2][0] = dbe_->book1D("convZFake"," Fake converted photon Z",zBin,zMin, zMax);
323  h_convPt_[2][0] = dbe_->book1D("convPtFake","Fake conversions Transverse Energy: all eta ", etBin,etMin, etMax);
324 
325  h_convRplot_ = dbe_->book1D("convRplot"," converted photon R",400, 0.,80.);
326  h_convZplot_ = dbe_->book1D("convZplot"," converted photon Z",320,-160.,160.);
327 
328 
329  histname = "convPtRes";
330  h_convPtRes_[0] = dbe_->book1D(histname+"All"," Conversion Pt rec/true : All ecal ", resBin,resMin, resMax);
331  h_convPtRes_[1] = dbe_->book1D(histname+"Barrel"," Conversion Pt rec/true : Barrel ",resBin,resMin, resMax);
332  h_convPtRes_[2] = dbe_->book1D(histname+"Endcap"," Conversion Pt rec/true : Endcap ",resBin,resMin, resMax);
333 
334 
335  histname="hInvMass";
336  h_invMass_[0][0]= dbe_->book1D(histname+"All_AllTracks"," Photons:Tracks from conversion: Pair invariant mass: all Ecal ",100, 0., 1.5);
337  h_invMass_[0][1]= dbe_->book1D(histname+"Barrel_AllTracks"," Photons:Tracks from conversion: Pair invariant mass: Barrel Ecal ",100, 0., 1.5);
338  h_invMass_[0][2]= dbe_->book1D(histname+"Endcap_AllTracks"," Photons:Tracks from conversion: Pair invariant mass: Endcap Ecal ",100, 0., 1.5);
339  //
340  h_invMass_[1][0]= dbe_->book1D(histname+"All_AssTracks"," Photons:Tracks from conversion: Pair invariant mass: all Ecal ",100, 0., 1.5);
341  h_invMass_[1][1]= dbe_->book1D(histname+"Barrel_AssTracks"," Photons:Tracks from conversion: Pair invariant mass: Barrel Ecal ",100, 0., 1.5);
342  h_invMass_[1][2]= dbe_->book1D(histname+"Endcap_AssTracks"," Photons:Tracks from conversion: Pair invariant mass: Endcap Ecal ",100, 0., 1.5);
343  //
344  h_invMass_[2][0]= dbe_->book1D(histname+"All_FakeTracks"," Photons:Tracks from conversion: Pair invariant mass: all Ecal ",100, 0., 1.5);
345  h_invMass_[2][1]= dbe_->book1D(histname+"Barrel_FakeTracks"," Photons:Tracks from conversion: Pair invariant mass: Barrel Ecal ",100, 0., 1.5);
346  h_invMass_[2][2]= dbe_->book1D(histname+"Endcap_FaleTracks"," Photons:Tracks from conversion: Pair invariant mass: Endcap Ecal ",100, 0., 1.5);
347 
348 
349 
350  histname="hDPhiTracksAtVtx";
351  h_DPhiTracksAtVtx_[0][0] =dbe_->book1D(histname+"All", " Photons:Tracks from conversions: #delta#phi Tracks at vertex: all Ecal",dPhiTracksBin,dPhiTracksMin,dPhiTracksMax);
352  h_DPhiTracksAtVtx_[0][1] =dbe_->book1D(histname+"Barrel", " Photons:Tracks from conversions: #delta#phi Tracks at vertex: Barrel Ecal",dPhiTracksBin,dPhiTracksMin,dPhiTracksMax);
353  h_DPhiTracksAtVtx_[0][2] =dbe_->book1D(histname+"Endcap", " Photons:Tracks from conversions: #delta#phi Tracks at vertex: Endcap Ecal",dPhiTracksBin,dPhiTracksMin,dPhiTracksMax);
354  h_DPhiTracksAtVtx_[1][0] =dbe_->book1D(histname+"All_Ass", " Photons:Tracks from conversions: #delta#phi Tracks at vertex: all Ecal",dPhiTracksBin,dPhiTracksMin,dPhiTracksMax);
355  h_DPhiTracksAtVtx_[1][1] =dbe_->book1D(histname+"Barrel_Ass", " Photons:Tracks from conversions: #delta#phi Tracks at vertex: Barrel Ecal",dPhiTracksBin,dPhiTracksMin,dPhiTracksMax);
356  h_DPhiTracksAtVtx_[1][2] =dbe_->book1D(histname+"Endcap_Ass", " Photons:Tracks from conversions: #delta#phi Tracks at vertex: Endcap Ecal",dPhiTracksBin,dPhiTracksMin,dPhiTracksMax);
357  h_DPhiTracksAtVtx_[2][0] =dbe_->book1D(histname+"All_Fakes", " Photons:Tracks from conversions: #delta#phi Tracks at vertex: all Ecal",dPhiTracksBin,dPhiTracksMin,dPhiTracksMax);
358  h_DPhiTracksAtVtx_[2][1] =dbe_->book1D(histname+"Barrel_Fakes", " Photons:Tracks from conversions: #delta#phi Tracks at vertex: Barrel Ecal",dPhiTracksBin,dPhiTracksMin,dPhiTracksMax);
359  h_DPhiTracksAtVtx_[2][2] =dbe_->book1D(histname+"Endcap_Fakes", " Photons:Tracks from conversions: #delta#phi Tracks at vertex: Endcap Ecal",dPhiTracksBin,dPhiTracksMin,dPhiTracksMax);
360 
361 
362 
363  histname="hDPhiTracksAtVtxVsEta";
364  h2_DPhiTracksAtVtxVsEta_ = dbe_->book2D(histname+"All"," Photons:Tracks from conversions: #delta#phi Tracks at vertex vs #eta",etaBin2,etaMin, etaMax,100, -0.5, 0.5);
365  histname="pDPhiTracksAtVtxVsEta";
366  p_DPhiTracksAtVtxVsEta_ = dbe_->bookProfile(histname+"All"," Photons:Tracks from conversions: #delta#phi Tracks at vertex vs #eta ",etaBin2,etaMin, etaMax, 100, -0.5, 0.5,"");
367 
368  histname="hDPhiTracksAtVtxVsR";
369  h2_DPhiTracksAtVtxVsR_ = dbe_->book2D(histname+"All"," Photons:Tracks from conversions: #delta#phi Tracks at vertex vs R",rBin,rMin, rMax,100, -0.5, 0.5);
370  histname="pDPhiTracksAtVtxVsR";
371  p_DPhiTracksAtVtxVsR_ = dbe_->bookProfile(histname+"All"," Photons:Tracks from conversions: #delta#phi Tracks at vertex vs R ",rBin,rMin, rMax,100, -0.5, 0.5,"");
372 
373 
374  histname="hDCotTracks";
375  h_DCotTracks_[0][0]= dbe_->book1D(histname+"All"," Photons:Tracks from conversions #delta cotg(#Theta) Tracks: all Ecal ",dCotTracksBin,dCotTracksMin,dCotTracksMax);
376  h_DCotTracks_[0][1]= dbe_->book1D(histname+"Barrel"," Photons:Tracks from conversions #delta cotg(#Theta) Tracks: Barrel Ecal ",dCotTracksBin,dCotTracksMin,dCotTracksMax);
377  h_DCotTracks_[0][2]= dbe_->book1D(histname+"Endcap"," Photons:Tracks from conversions #delta cotg(#Theta) Tracks: Endcap Ecal ",dCotTracksBin,dCotTracksMin,dCotTracksMax);
378  h_DCotTracks_[1][0]= dbe_->book1D(histname+"All_Ass"," Photons:Tracks from conversions #delta cotg(#Theta) Tracks: all Ecal ",dCotTracksBin,dCotTracksMin,dCotTracksMax);
379  h_DCotTracks_[1][1]= dbe_->book1D(histname+"Barrel_Ass"," Photons:Tracks from conversions #delta cotg(#Theta) Tracks: Barrel Ecal ",dCotTracksBin,dCotTracksMin,dCotTracksMax);
380  h_DCotTracks_[1][2]= dbe_->book1D(histname+"Endcap_Ass"," Photons:Tracks from conversions #delta cotg(#Theta) Tracks: Endcap Ecal ",dCotTracksBin,dCotTracksMin,dCotTracksMax);
381  h_DCotTracks_[2][0]= dbe_->book1D(histname+"All_Fakes"," Photons:Tracks from conversions #delta cotg(#Theta) Tracks: all Ecal ",dCotTracksBin,dCotTracksMin,dCotTracksMax);
382  h_DCotTracks_[2][1]= dbe_->book1D(histname+"Barrel_Fakes"," Photons:Tracks from conversions #delta cotg(#Theta) Tracks: Barrel Ecal ",dCotTracksBin,dCotTracksMin,dCotTracksMax);
383  h_DCotTracks_[2][2]= dbe_->book1D(histname+"Endcap_Fakes"," Photons:Tracks from conversions #delta cotg(#Theta) Tracks: Endcap Ecal ",dCotTracksBin,dCotTracksMin,dCotTracksMax);
384 
385 
386  histname="hDCotTracksVsEta";
387  h2_DCotTracksVsEta_ = dbe_->book2D(histname+"All"," Photons:Tracks from conversions: #delta cotg(#Theta) Tracks vs #eta",etaBin2,etaMin, etaMax,100, -0.2, 0.2);
388  histname="pDCotTracksVsEta";
389  p_DCotTracksVsEta_ = dbe_->bookProfile(histname+"All"," Photons:Tracks from conversions: #delta cotg(#Theta) Tracks vs #eta ",etaBin2,etaMin, etaMax, 100, -0.2, 0.2,"");
390 
391  histname="hDCotTracksVsR";
392  h2_DCotTracksVsR_ = dbe_->book2D(histname+"All"," Photons:Tracks from conversions: #delta cotg(#Theta) Tracks at vertex vs R",rBin,rMin, rMax,100, -0.2, 0.2);
393  histname="pDCotTracksVsR";
394  p_DCotTracksVsR_ = dbe_->bookProfile(histname+"All"," Photons:Tracks from conversions: #delta cotg(#Theta) Tracks at vertex vs R ",rBin,rMin, rMax,100, -0.2, 0.2,"");
395 
396 
397  histname="hDistMinAppTracks";
398  h_distMinAppTracks_[0][0]= dbe_->book1D(histname+"All"," Photons:Tracks from conversions Min Approach Dist Tracks: all Ecal ",120, -0.5, 1.0);
399  h_distMinAppTracks_[0][1]= dbe_->book1D(histname+"Barrel"," Photons:Tracks from conversions Min Approach Dist Tracks: Barrel Ecal ",120, -0.5, 1.0);
400  h_distMinAppTracks_[0][2]= dbe_->book1D(histname+"Endcap"," Photons:Tracks from conversions Min Approach Dist Tracks: Endcap Ecal ",120, -0.5, 1.0);
401  h_distMinAppTracks_[1][0]= dbe_->book1D(histname+"All_Ass"," Photons:Tracks from conversions Min Approach Dist Tracks: all Ecal ",120, -0.5, 1.0);
402  h_distMinAppTracks_[1][1]= dbe_->book1D(histname+"Barrel_Ass"," Photons:Tracks from conversions Min Approach Dist Tracks: Barrel Ecal ",120, -0.5, 1.0);
403  h_distMinAppTracks_[1][2]= dbe_->book1D(histname+"Endcap_Ass"," Photons:Tracks from conversions Min Approach Dist Tracks: Endcap Ecal ",120, -0.5, 1.0);
404  h_distMinAppTracks_[2][0]= dbe_->book1D(histname+"All_Fakes"," Photons:Tracks from conversions Min Approach Dist Tracks: all Ecal ",120, -0.5, 1.0);
405  h_distMinAppTracks_[2][1]= dbe_->book1D(histname+"Barrel_Fakes"," Photons:Tracks from conversions Min Approach Dist Tracks: Barrel Ecal ",120, -0.5, 1.0);
406  h_distMinAppTracks_[2][2]= dbe_->book1D(histname+"Endcap_Fakes"," Photons:Tracks from conversions Min Approach Dist Tracks: Endcap Ecal ",120, -0.5, 1.0);
407 
408 
409  h_convVtxRvsZ_[0] = dbe_->book2D("convVtxRvsZAll"," Photon Reco conversion vtx position",zBinForXray, zMinForXray, zMaxForXray, rBinForXray, rMinForXray, rMaxForXray);
410  h_convVtxRvsZ_[1] = dbe_->book2D("convVtxRvsZBarrel"," Photon Reco conversion vtx position",zBinForXray, zMinForXray, zMaxForXray, rBinForXray, rMinForXray, rMaxForXray);
411  h_convVtxRvsZ_[2] = dbe_->book2D("convVtxRvsZEndcap"," Photon Reco conversion vtx position",zBin2ForXray, zMinForXray, zMaxForXray, rBinForXray, rMinForXray, rMaxForXray);
412  h_convVtxYvsX_ = dbe_->book2D("convVtxYvsXTrkBarrel"," Photon Reco conversion vtx position, (x,y) eta<1 ", 1000, -60., 60., 1000, -60., 60.);
414  h_convVtxRvsZ_zoom_[0] = dbe_->book2D("convVtxRvsZBarrelZoom1"," Photon Reco conversion vtx position",zBinForXray, zMinForXray, zMaxForXray, rBinForXray, -10., 40.);
415  h_convVtxRvsZ_zoom_[1] = dbe_->book2D("convVtxRvsZBarrelZoom2"," Photon Reco conversion vtx position",zBinForXray, zMinForXray, zMaxForXray, rBinForXray, -10., 20.);
416  h_convVtxYvsX_zoom_[0] = dbe_->book2D("convVtxYvsXTrkBarrelZoom1"," Photon Reco conversion vtx position, (x,y) eta<1 ",100, -40., 40., 100, -40., 40.);
417  h_convVtxYvsX_zoom_[1] = dbe_->book2D("convVtxYvsXTrkBarrelZoom2"," Photon Reco conversion vtx position, (x,y) eta<1 ",100, -20., 20., 100, -20., 20.);
418 
419  h_convVtxdR_ = dbe_->book1D("convVtxdR"," Photon Reco conversion vtx dR",100, -10.,10.);
420  h_convVtxdX_ = dbe_->book1D("convVtxdX"," Photon Reco conversion vtx dX",100, -10.,10.);
421  h_convVtxdY_ = dbe_->book1D("convVtxdY"," Photon Reco conversion vtx dY",100, -10.,10.);
422  h_convVtxdZ_ = dbe_->book1D("convVtxdZ"," Photon Reco conversion vtx dZ",100, -20.,20.);
423 
424  h_convVtxdPhi_ = dbe_->book1D("convVtxdPhi"," Photon Reco conversion vtx dPhi",100, -0.01,0.01);
425  h_convVtxdEta_ = dbe_->book1D("convVtxdEta"," Photon Reco conversion vtx dEta",100, -0.5,0.5);
426 
427  h_convVtxdR_barrel_ = dbe_->book1D("convVtxdR_barrel"," Photon Reco conversion vtx dR, |eta|<=1.2",100, -10.,10.);
428  h_convVtxdX_barrel_ = dbe_->book1D("convVtxdX_barrel"," Photon Reco conversion vtx dX, |eta|<=1.2",100, -10.,10.);
429  h_convVtxdY_barrel_ = dbe_->book1D("convVtxdY_barrel"," Photon Reco conversion vtx dY, |eta|<=1.2 ",100, -10.,10.);
430  h_convVtxdZ_barrel_ = dbe_->book1D("convVtxdZ_barrel"," Photon Reco conversion vtx dZ, |eta|<=1.2,",100, -20.,20.);
431 
432  h_convVtxdR_endcap_ = dbe_->book1D("convVtxdR_endcap"," Photon Reco conversion vtx dR, |eta|>1.2 ",100, -10.,10.);
433  h_convVtxdX_endcap_ = dbe_->book1D("convVtxdX_endcap"," Photon Reco conversion vtx dX, |eta|>1.2",100, -10.,10.);
434  h_convVtxdY_endcap_ = dbe_->book1D("convVtxdY_endcap"," Photon Reco conversion vtx dY, |eta|>1.2",100, -10.,10.);
435  h_convVtxdZ_endcap_ = dbe_->book1D("convVtxdZ_endcap"," Photon Reco conversion vtx dZ, |eta|>1.2",100, -20.,20.);
436 
437 
438 
439  h2_convVtxdRVsR_ = dbe_->book2D("h2ConvVtxdRVsR"," Conversion vtx dR vsR" ,rBin,rMin, rMax,100, -20.,20.);
440  h2_convVtxdRVsEta_ = dbe_->book2D("h2ConvVtxdRVsEta","Conversion vtx dR vs Eta" ,etaBin2,etaMin, etaMax,100, -20.,20.);
441 
442  p_convVtxdRVsR_ = dbe_->bookProfile("pConvVtxdRVsR"," Conversion vtx dR vsR" ,rBin,rMin, rMax ,100, -20.,20., "");
443  p_convVtxdRVsEta_ = dbe_->bookProfile("pConvVtxdRVsEta","Conversion vtx dR vs Eta" ,etaBin2,etaMin, etaMax, 100, -20.,20., "");
444  p_convVtxdXVsX_ = dbe_->bookProfile("pConvVtxdXVsX","Conversion vtx dX vs X" ,120,-60, 60 ,100, -20.,20., "");
445  p_convVtxdYVsY_ = dbe_->bookProfile("pConvVtxdYVsY","Conversion vtx dY vs Y" ,120,-60, 60 ,100, -20.,20., "");
446  p_convVtxdZVsZ_ = dbe_->bookProfile("pConvVtxdZVsZ","Conversion vtx dZ vs Z" ,zBin,zMin,zMax ,100, -20.,20., "");
447 
448  p_convVtxdZVsR_ = dbe_->bookProfile("pConvVtxdZVsR","Conversion vtx dZ vs R" ,rBin,rMin,rMax ,100, -20.,20., "");
449  p2_convVtxdRVsRZ_ = dbe_->bookProfile2D("p2ConvVtxdRVsRZ","Conversion vtx dR vs RZ" ,zBin,zMin, zMax,rBin,rMin,rMax,100, 0.,20.,"s");
450  p2_convVtxdZVsRZ_ = dbe_->bookProfile2D("p2ConvVtxdZVsRZ","Conversion vtx dZ vs RZ" ,zBin,zMin, zMax,rBin,rMin,rMax,100, 0.,20.,"s");
451 
452 
453 
454  h2_convVtxRrecVsTrue_ = dbe_->book2D("h2ConvVtxRrecVsTrue","Photon Reco conversion vtx R rec vs true" ,rBin,rMin, rMax,rBin,rMin, rMax);
455 
456  histname="vtxChi2Prob";
457  h_vtxChi2Prob_[0][0] = dbe_->book1D(histname+"All","vertex #chi^{2} all", 100, 0., 1.);
458  h_vtxChi2Prob_[0][1] = dbe_->book1D(histname+"Barrel","vertex #chi^{2} barrel", 100, 0., 1.);
459  h_vtxChi2Prob_[0][2] = dbe_->book1D(histname+"Endcap","vertex #chi^{2} endcap", 100, 0., 1.);
460  h_vtxChi2Prob_[1][0] = dbe_->book1D(histname+"All_Ass","vertex #chi^{2} all", 100, 0., 1.);
461  h_vtxChi2Prob_[1][1] = dbe_->book1D(histname+"Barrel_Ass","vertex #chi^{2} barrel", 100, 0., 1.);
462  h_vtxChi2Prob_[1][2] = dbe_->book1D(histname+"Endcap_Ass","vertex #chi^{2} endcap", 100, 0., 1.);
463  h_vtxChi2Prob_[2][0] = dbe_->book1D(histname+"All_Fakes","vertex #chi^{2} all", 100, 0., 1.);
464  h_vtxChi2Prob_[2][1] = dbe_->book1D(histname+"Barrel_Fakes","vertex #chi^{2} barrel", 100, 0., 1.);
465  h_vtxChi2Prob_[2][2] = dbe_->book1D(histname+"Endcap_Fakes","vertex #chi^{2} endcap", 100, 0., 1.);
466 
467 
468  h_zPVFromTracks_[1] = dbe_->book1D("zPVFromTracks"," Photons: PV z from conversion tracks",100, -25., 25.);
469  h_dzPVFromTracks_[1] = dbe_->book1D("dzPVFromTracks"," Photons: PV Z_rec - Z_true from conversion tracks",100, -5., 5.);
470  h2_dzPVVsR_ = dbe_->book2D("h2dzPVVsR","Photon Reco conversions: dz(PV) vs R" ,rBin,rMin, rMax,100, -3.,3.);
471  p_dzPVVsR_ = dbe_->bookProfile("pdzPVVsR","Photon Reco conversions: dz(PV) vs R" ,rBin,rMin, rMax, 100, -3.,3.,"");
472 
473 
475  histname="nHits";
476  nHits_[0] = dbe_->book2D(histname+"AllTracks","Photons:Tracks from conversions: # of hits all tracks",etaBin,etaMin, etaMax,30,0., 30.);
477  nHits_[1] = dbe_->book2D(histname+"AllTracks_Ass","Photons:Tracks from conversions: # of hits all tracks ass",etaBin,etaMin, etaMax,30,0., 30.);
478  nHits_[2] = dbe_->book2D(histname+"AllTracks_Fakes","Photons:Tracks from conversions: # of hits all tracks fakes",etaBin,etaMin, etaMax,30,0., 30.);
479 
480 
481  histname="nHitsVsEta";
482  nHitsVsEta_[0] = dbe_->book2D(histname+"AllTracks","Photons:Tracks from conversions: # of hits vs #eta all tracks",etaBin,etaMin, etaMax,30,0., 30.);
483  nHitsVsEta_[1] = dbe_->book2D(histname+"AllTracks_Ass","Photons:Tracks from conversions: # of hits vs #eta all tracks",etaBin,etaMin, etaMax,30,0., 30.);
484  nHitsVsEta_[2] = dbe_->book2D(histname+"AllTracks_Fakes","Photons:Tracks from conversions: # of hits vs #eta all tracks",etaBin,etaMin, etaMax,30,0., 30.);
485  histname="h_nHitsVsEta";
486  p_nHitsVsEta_[0] = dbe_->bookProfile(histname+"AllTracks","Photons:Tracks from conversions: # of hits vs #eta all tracks",etaBin,etaMin, etaMax, 30,-0.5, 29.5,"");
487  p_nHitsVsEta_[1] = dbe_->bookProfile(histname+"AllTracks_Ass","Photons:Tracks from conversions: # of hits vs #eta all tracks",etaBin,etaMin, etaMax, 30,-0.5, 29.5,"");
488  p_nHitsVsEta_[2] = dbe_->bookProfile(histname+"AllTracks_Fakes","Photons:Tracks from conversions: # of hits vs #eta all tracks",etaBin,etaMin, etaMax, 30,-0.5, 29.5,"");
489 
490 
491  histname="nHitsVsR";
492  nHitsVsR_[0] = dbe_->book2D(histname+"AllTracks","Photons:Tracks from conversions: # of hits vs radius all tracks" ,rBin,rMin, rMax,30,0., 30.);
493  nHitsVsR_[1] = dbe_->book2D(histname+"AllTracks_Ass","Photons:Tracks from conversions: # of hits vs radius all tracks" ,rBin,rMin, rMax,30,0., 30.);
494  nHitsVsR_[2] = dbe_->book2D(histname+"AllTracks_Fakes","Photons:Tracks from conversions: # of hits vs radius all tracks" ,rBin,rMin, rMax,30,0., 30.);
495 
496  histname="h_nHitsVsR";
497  p_nHitsVsR_[0] = dbe_->bookProfile(histname+"AllTracks","Photons:Tracks from conversions: # of hits vs radius all tracks",rBin,rMin, rMax, 30,-0.5, 29.5,"");
498  p_nHitsVsR_[1] = dbe_->bookProfile(histname+"AllTracks_Ass","Photons:Tracks from conversions: # of hits vs radius all tracks",rBin,rMin, rMax, 30,-0.5, 29.5,"");
499  p_nHitsVsR_[2] = dbe_->bookProfile(histname+"AllTracks_Fakes","Photons:Tracks from conversions: # of hits vs radius all tracks",rBin,rMin, rMax, 30,-0.5, 29.5,"");
500 
501  histname="tkChi2";
502  h_tkChi2_[0] = dbe_->book1D(histname+"AllTracks","Photons:Tracks from conversions: #chi^{2} of all tracks", 100, chi2Min, chi2Max);
503  h_tkChi2_[1] = dbe_->book1D(histname+"AllTracks_Ass","Photons:Tracks from conversions: #chi^{2} of all tracks", 100, chi2Min, chi2Max);
504  h_tkChi2_[2] = dbe_->book1D(histname+"AllTracks_Fakes","Photons:Tracks from conversions: #chi^{2} of all tracks", 100, chi2Min, chi2Max);
505 
506  histname="tkChi2Large";
507  h_tkChi2Large_[0] = dbe_->book1D(histname+"AllTracks","Photons:Tracks from conversions: #chi^{2} of all tracks", 1000, 0., 5000.0);
508  h_tkChi2Large_[1] = dbe_->book1D(histname+"AllTracks_Ass","Photons:Tracks from conversions: #chi^{2} of all tracks", 1000, 0., 5000.0);
509  h_tkChi2Large_[2] = dbe_->book1D(histname+"AllTracks_Fakes","Photons:Tracks from conversions: #chi^{2} of all tracks", 1000, 0., 5000.0);
510 
511 
512  histname="h2Chi2VsEta";
513  h2_Chi2VsEta_[0]=dbe_->book2D(histname+"All"," Reco Track #chi^{2} vs #eta: All ",etaBin2,etaMin, etaMax,100, chi2Min, chi2Max);
514  h2_Chi2VsEta_[1]=dbe_->book2D(histname+"All_Ass"," Reco Track #chi^{2} vs #eta: All ",etaBin2,etaMin, etaMax,100, chi2Min, chi2Max);
515  h2_Chi2VsEta_[2]=dbe_->book2D(histname+"All_Fakes"," Reco Track #chi^{2} vs #eta: All ",etaBin2,etaMin, etaMax,100, chi2Min, chi2Max);
516  histname="pChi2VsEta";
517  p_Chi2VsEta_[0]=dbe_->bookProfile(histname+"All"," Reco Track #chi^{2} vs #eta : All ",etaBin2,etaMin, etaMax, 100, chi2Min, chi2Max,"");
518  p_Chi2VsEta_[1]=dbe_->bookProfile(histname+"All_Ass"," Reco Track #chi^{2} vs #eta : All ",etaBin2,etaMin, etaMax, 100, chi2Min, chi2Max,"");
519  p_Chi2VsEta_[2]=dbe_->bookProfile(histname+"All_Fakes"," Reco Track #chi^{2} vs #eta : All ",etaBin2,etaMin, etaMax, 100, chi2Min, chi2Max,"");
520 
521  histname="h2Chi2VsR";
522  h2_Chi2VsR_[0]=dbe_->book2D(histname+"All"," Reco Track #chi^{2} vs R: All ",rBin,rMin, rMax,100,chi2Min, chi2Max);
523  h2_Chi2VsR_[1]=dbe_->book2D(histname+"All_Ass"," Reco Track #chi^{2} vs R: All ",rBin,rMin, rMax,100,chi2Min, chi2Max);
524  h2_Chi2VsR_[2]=dbe_->book2D(histname+"All_Fakes"," Reco Track #chi^{2} vs R: All ",rBin,rMin, rMax,100,chi2Min, chi2Max);
525  histname="pChi2VsR";
526  p_Chi2VsR_[0]=dbe_->bookProfile(histname+"All"," Reco Track #chi^{2} vas R : All ",rBin,rMin,rMax, 100,chi2Min, chi2Max,"");
527  p_Chi2VsR_[1]=dbe_->bookProfile(histname+"All_Ass"," Reco Track #chi^{2} vas R : All ",rBin,rMin,rMax, 100,chi2Min, chi2Max,"");
528  p_Chi2VsR_[2]=dbe_->bookProfile(histname+"All_Fakes"," Reco Track #chi^{2} vas R : All ",rBin,rMin,rMax, 100,chi2Min, chi2Max,"");
529 
530  histname="hTkD0";
531  h_TkD0_[0]=dbe_->book1D(histname+"All"," Reco Track D0*q: All ",200,-0.1,60);
532  h_TkD0_[1]=dbe_->book1D(histname+"All_Ass"," Reco Track D0*q: Barrel ",200,-0.1,60);
533  h_TkD0_[2]=dbe_->book1D(histname+"All_Fakes"," Reco Track D0*q: Endcap ",200,-0.1,60);
534 
535 
536 
537  histname="hTkPtPull";
538  h_TkPtPull_[0]=dbe_->book1D(histname+"All"," Reco Track Pt pull: All ",100, -20., 10.);
539  histname="hTkPtPull";
540  h_TkPtPull_[1]=dbe_->book1D(histname+"Barrel"," Reco Track Pt pull: Barrel ",100, -20., 10.);
541  histname="hTkPtPull";
542  h_TkPtPull_[2]=dbe_->book1D(histname+"Endcap"," Reco Track Pt pull: Endcap ",100, -20., 10.);
543 
544  histname="h2TkPtPullEta";
545  h2_TkPtPull_[0]=dbe_->book2D(histname+"All"," Reco Track Pt pull: All ",etaBin2,etaMin, etaMax,100, -20., 10.);
546  histname="pTkPtPullEta";
547  p_TkPtPull_[0]=dbe_->bookProfile(histname+"All"," Reco Track Pt pull: All ",etaBin2,etaMin, etaMax, 100, -20., 10., " ");
548 
549 
550  histname="PtRecVsPtSim";
551  h2_PtRecVsPtSim_[0]=dbe_->book2D(histname+"All", "Pt Rec vs Pt sim: All ", etBin,etMin,etMax,etBin,etMin, etMax);
552  h2_PtRecVsPtSim_[1]=dbe_->book2D(histname+"Barrel", "Pt Rec vs Pt sim: Barrel ", etBin,etMin,etMax,etBin,etMin, etMax);
553  h2_PtRecVsPtSim_[2]=dbe_->book2D(histname+"Endcap", "Pt Rec vs Pt sim: Endcap ", etBin,etMin,etMax,etBin,etMin, etMax);
554 
555  histname="photonPtRecVsPtSim";
556  h2_photonPtRecVsPtSim_=dbe_->book2D(histname+"All", "Pt Rec vs Pt sim: All ", etBin,etMin,etMax,etBin,etMin, etMax);
557 
558 
559 
560  h_match_= dbe_->book1D("h_match"," ", 3, -0.5,2.5);
561 
562 
563  } // if DQM
564 
565 
566 
567 }
568 
569 
570 
571  void TkConvValidator::beginRun (edm::Run const & r, edm::EventSetup const & theEventSetup) {
572 
573  //get magnetic field
574  edm::LogInfo("ConvertedPhotonProducer") << " get magnetic field" << "\n";
575  theEventSetup.get<IdealMagneticFieldRecord>().get(theMF_);
576 
577 
578  edm::ESHandle<TrackAssociatorBase> theHitsAssociator;
579  theEventSetup.get<TrackAssociatorRecord>().get("TrackAssociatorByHits",theHitsAssociator);
580  theTrackAssociator_ = (TrackAssociatorBase *) theHitsAssociator.product();
581 
582 
583 
584 
585  thePhotonMCTruthFinder_ = new PhotonMCTruthFinder();
586 
587 }
588 
589 void TkConvValidator::endRun (edm::Run& r, edm::EventSetup const & theEventSetup) {
590 
591  delete thePhotonMCTruthFinder_;
592 
593 }
594 
595 
596 
598 
599  using namespace edm;
600  // const float etaPhiDistance=0.01;
601  // Fiducial region
602  // const float TRK_BARL =0.9;
603  const float BARL = 1.4442; // DAQ TDR p.290
604  // const float END_LO = 1.566; // unused
605  const float END_HI = 2.5;
606  // Electron mass
607  // const Float_t mElec= 0.000511; // unused
608 
609 
610  nEvt_++;
611  LogInfo("TkConvValidator") << "TkConvValidator Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
612  // std::cout << "TkConvValidator Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
613 
614 
615  // get the geometry from the event setup:
616  esup.get<CaloGeometryRecord>().get(theCaloGeom_);
617 
618 
619  // Transform Track into TransientTrack (needed by the Vertex fitter)
621  esup.get<TransientTrackRecord>().get("TransientTrackBuilder",theTTB);
622 
623 
626  e.getByLabel(conversionCollectionProducer_, conversionCollection_ , convHandle);
627  const reco::ConversionCollection convCollection = *(convHandle.product());
628  if (!convHandle.isValid()) {
629  edm::LogError("ConversionsProducer") << "Error! Can't get the collection "<< std::endl;
630  return;
631  }
632 
633  // offline Primary vertex
636  e.getByLabel("offlinePrimaryVertices", vertexHandle);
637  if (!vertexHandle.isValid()) {
638  edm::LogError("TrackerOnlyConversionProducer") << "Error! Can't get the product primary Vertex Collection "<< "\n";
639  } else {
640  vertexCollection = *(vertexHandle.product());
641  }
642  reco::Vertex the_pvtx;
643  bool valid_pvtx = false;
644  if (!vertexCollection.empty()){
645  the_pvtx = *(vertexCollection.begin());
646  //asking for one good vertex
647  if (the_pvtx.isValid() && fabs(the_pvtx.position().z())<=15 && the_pvtx.position().Rho()<=2){
648  valid_pvtx = true;
649  }
650  }
651 
652 
653  //get tracker geometry for hits positions
656  const TrackerGeometry* trackerGeom = tracker.product();
657 
658 
659 
661  //get simtrack info
662  std::vector<SimTrack> theSimTracks;
663  std::vector<SimVertex> theSimVertices;
664 
667  e.getByLabel("g4SimHits",SimTk);
668  e.getByLabel("g4SimHits",SimVtx);
669 
670  bool useTP= parameters_.getParameter<bool>("useTP");
671  TrackingParticleCollection tpForEfficiency;
672  TrackingParticleCollection tpForFakeRate;
674  edm::Handle<TrackingParticleCollection> TPHandleForFakeRate;
675  if ( useTP) {
676  e.getByLabel("tpSelecForEfficiency",TPHandleForEff);
677  tpForEfficiency = *(TPHandleForEff.product());
678  e.getByLabel("tpSelecForFakeRate",TPHandleForFakeRate);
679  tpForFakeRate = *(TPHandleForFakeRate.product());
680  }
681 
682 
683 
684  theSimTracks.insert(theSimTracks.end(),SimTk->begin(),SimTk->end());
685  theSimVertices.insert(theSimVertices.end(),SimVtx->begin(),SimVtx->end());
686  std::vector<PhotonMCTruth> mcPhotons=thePhotonMCTruthFinder_->find (theSimTracks, theSimVertices);
687 
689  e.getByLabel("generator",hepMC);
690  // const HepMC::GenEvent *myGenEvent = hepMC->GetEvent(); // unused
691 
692 
693  // get generated jets
695  e.getByLabel("iterativeCone5GenJets","",GenJetsHandle);
696  reco::GenJetCollection genJetCollection = *(GenJetsHandle.product());
697 
698 
699 
700  // ################ SIM to RECO ######################### //
701  std::map<const reco::Track*,TrackingParticleRef> myAss;
702  std::map<const reco::Track*,TrackingParticleRef>::const_iterator itAss;
703 
704  for ( std::vector<PhotonMCTruth>::const_iterator mcPho=mcPhotons.begin(); mcPho !=mcPhotons.end(); mcPho++) {
705 
706  mcConvPt_= (*mcPho).fourMomentum().et();
707  float mcPhi= (*mcPho).fourMomentum().phi();
708  mcPhi_= phiNormalization(mcPhi);
709  mcEta_= (*mcPho).fourMomentum().pseudoRapidity();
710  mcEta_ = etaTransformation(mcEta_, (*mcPho).primaryVertex().z() );
711  mcConvR_= (*mcPho).vertex().perp();
712  mcConvX_= (*mcPho).vertex().x();
713  mcConvY_= (*mcPho).vertex().y();
714  mcConvZ_= (*mcPho).vertex().z();
715  mcConvEta_= (*mcPho).vertex().eta();
716  mcConvPhi_= (*mcPho).vertex().phi();
717 
718  if ( fabs(mcEta_) > END_HI ) continue;
719 
720  if (mcConvPt_<minPhoPtForEffic) continue;
721  if (fabs(mcEta_)>maxPhoEtaForEffic) continue;
722  if (fabs(mcConvZ_)>maxPhoZForEffic) continue;
723  if (mcConvR_>maxPhoRForEffic) continue;
725 
726  bool goodSimConversion=false;
727  bool visibleConversion=false;
728  bool visibleConversionsWithTwoSimTracks=false;
729  if ( (*mcPho).isAConversion() == 1 ) {
730  nSimConv_[0]++;
731  h_AllSimConv_[0]->Fill( mcEta_ ) ;
732  h_AllSimConv_[1]->Fill( mcPhi_ );
733  h_AllSimConv_[2]->Fill( mcConvR_ );
734  h_AllSimConv_[3]->Fill( mcConvZ_ );
735  h_AllSimConv_[4]->Fill( (*mcPho).fourMomentum().et());
736 
737  if ( mcConvR_ <15) h_SimConvEtaPix_[0]->Fill( mcEta_ ) ;
738 
739  if ( ( fabs(mcEta_) <= BARL && mcConvR_ <85 ) ||
740  ( fabs(mcEta_) > BARL && fabs(mcEta_) <=END_HI && fabs( (*mcPho).vertex().z() ) < 210 ) ) visibleConversion=true;
741 
742  theConvTP_.clear();
743  // std::cout << " TkConvValidator TrackingParticles TrackingParticleCollection size "<< trackingParticles.size() << "\n";
744  //duplicated TP collections for two associations
745  for(size_t i = 0; i < tpForEfficiency.size(); ++i){
746  TrackingParticleRef tp (TPHandleForEff,i);
747  if ( fabs( tp->vx() - (*mcPho).vertex().x() ) < 0.0001 &&
748  fabs( tp->vy() - (*mcPho).vertex().y() ) < 0.0001 &&
749  fabs( tp->vz() - (*mcPho).vertex().z() ) < 0.0001) {
750  theConvTP_.push_back( tp );
751  }
752  }
753  //std::cout << " TkConvValidator theConvTP_ size " << theConvTP_.size() << std::endl;
754 
755  if ( theConvTP_.size() == 2 ) visibleConversionsWithTwoSimTracks=true;
756  goodSimConversion=false;
757 
758  if ( visibleConversion && visibleConversionsWithTwoSimTracks ) goodSimConversion=true;
759  if ( goodSimConversion ) {
760  nSimConv_[1]++;
761  h_VisSimConv_[0]->Fill( mcEta_ ) ;
762  h_VisSimConv_[1]->Fill( mcPhi_ );
763  h_VisSimConv_[2]->Fill( mcConvR_ );
764  h_VisSimConv_[3]->Fill( mcConvZ_ );
765  h_VisSimConv_[4]->Fill( (*mcPho).fourMomentum().et());
766 
767  }
768 
769  for ( edm::RefVector<TrackingParticleCollection>::iterator iTrk=theConvTP_.begin(); iTrk!=theConvTP_.end(); ++iTrk) {
770  h_simTkPt_ -> Fill ( (*iTrk)->pt() );
771  h_simTkEta_ -> Fill ( (*iTrk)->eta() );
772  }
773 
774 
775  }
776 
777  if ( ! (visibleConversion && visibleConversionsWithTwoSimTracks ) ) continue;
778 
779  h_simConvVtxRvsZ_[0] ->Fill ( fabs (mcConvZ_), mcConvR_ ) ;
780  if ( fabs(mcEta_) <=1.) {
781  h_simConvVtxRvsZ_[1] ->Fill ( fabs (mcConvZ_), mcConvR_ ) ;
782  h_simConvVtxYvsX_ ->Fill ( mcConvX_, mcConvY_ ) ;
783  }
784  else
785  h_simConvVtxRvsZ_[2] ->Fill ( fabs (mcConvZ_), mcConvR_ ) ;
786 
787  //std::cout << " TkConvValidator theConvTP_ size " << theConvTP_.size() << std::endl;
788  for ( edm::RefVector<TrackingParticleCollection>::iterator iTP= theConvTP_.begin(); iTP!=theConvTP_.end(); iTP++)
789  {
790  // std::cout << " SIM to RECO TP vertex " << (*iTP)->vx() << " " << (*iTP)->vy() << " " << (*iTP)->vz() << " pt " << (*iTP)->pt() << std::endl;
791  }
792 
793 
795  for (reco::ConversionCollection::const_iterator conv = convHandle->begin();conv!=convHandle->end();++conv) {
796 
797  const reco::Conversion aConv = (*conv);
798  if ( mergedTracks_ ) {
800  } else {
802  }
803 
804 
805  //problematic?
806  std::vector<edm::RefToBase<reco::Track> > tracks = aConv.tracks();
807 
808  const reco::Vertex& vtx = aConv.conversionVertex();
809  //requires two tracks and a valid vertex
810  if (tracks.size() !=2 || !(vtx.isValid())) continue;
811 
812  // bool phoIsInBarrel=false; // unused
813  // bool phoIsInEndcap=false; // unused
814  RefToBase<reco::Track> tfrb1 = aConv.tracks().front();
815  RefToBase<reco::Track> tfrb2 = aConv.tracks().back();
816  //reco::TrackRef tk1 = aConv.tracks().front();
817  //reco::TrackRef tk2 = aConv.tracks().back();
818  //std::cout << "SIM to RECO conversion track pt " << tk1->pt() << " " << tk2->pt() << endl;
819  //
820  //Use two RefToBaseVector and do two association actions to avoid that if two tracks from different collection
822  tc1.push_back(tfrb1);
823  tc2.push_back(tfrb2);
824  bool isAssociated = false;
825  reco::SimToRecoCollection q1 = theTrackAssociator_->associateSimToReco(tc1,theConvTP_,&e);
826  reco::SimToRecoCollection q2 = theTrackAssociator_->associateSimToReco(tc2,theConvTP_,&e);
827  //try {
828  std::vector<std::pair<RefToBase<reco::Track>, double> > trackV1, trackV2;
829 
830  int tp_1 = 0, tp_2 = 1;//the index of associated tp in theConvTP_ for two tracks
831  if (q1.find(theConvTP_[0])!=q1.end()){
832  trackV1 = (std::vector<std::pair<RefToBase<reco::Track>, double> >) q1[theConvTP_[0]];
833  } else if (q1.find(theConvTP_[1])!=q1.end()){
834  trackV1 = (std::vector<std::pair<RefToBase<reco::Track>, double> >) q1[theConvTP_[1]];
835  tp_1 = 1;
836  }
837  if (q2.find(theConvTP_[1])!=q2.end()){
838  trackV2 = (std::vector<std::pair<RefToBase<reco::Track>, double> >) q2[theConvTP_[1]];
839  } else if (q2.find(theConvTP_[0])!=q2.end()){
840  trackV2 = (std::vector<std::pair<RefToBase<reco::Track>, double> >) q2[theConvTP_[0]];
841  tp_2 = 0;
842  }
843  if (!(trackV1.size()&&trackV2.size()))
844  continue;
845  if (tp_1 == tp_2) continue;
846 
847  edm::RefToBase<reco::Track> tr1 = trackV1.front().first;
848  edm::RefToBase<reco::Track> tr2 = trackV2.front().first;
849  //std::cout << "associated tp1 " <<theConvTP_[0]->pt() << " to track with pT=" << tr1->pt() << " " << (tr1.get())->pt() << endl;
850  //std::cout << "associated tp2 " <<theConvTP_[1]->pt() << " to track with pT=" << tr2->pt() << " " << (tr2.get())->pt() << endl;
851  myAss.insert( std::make_pair (tr1.get(),theConvTP_[tp_1] ) );
852  myAss.insert( std::make_pair (tr2.get(),theConvTP_[tp_2]) );
853 
854  //} catch (Exception event) {
855  //cout << "continue: " << event.what() << endl;
856  // continue;
857  //}
858 
859 
860  isAssociated = true;
861 
863  h_SimConvTwoMTracks_[0]->Fill( mcEta_ ) ;
864  h_SimConvTwoMTracks_[1]->Fill( mcPhi_ );
865  h_SimConvTwoMTracks_[2]->Fill( mcConvR_ );
866  h_SimConvTwoMTracks_[3]->Fill( mcConvZ_ );
867  h_SimConvTwoMTracks_[4]->Fill( (*mcPho).fourMomentum().et());
868 
869  float chi2Prob = ChiSquaredProbability( aConv.conversionVertex().chi2(), aConv.conversionVertex().ndof() );
870  if ( chi2Prob > 0) {
871  h_SimConvTwoMTracksAndVtxPGT0_[0]->Fill( mcEta_ ) ;
872  h_SimConvTwoMTracksAndVtxPGT0_[1]->Fill( mcPhi_ );
873  h_SimConvTwoMTracksAndVtxPGT0_[2]->Fill( mcConvR_ );
874  h_SimConvTwoMTracksAndVtxPGT0_[3]->Fill( mcConvZ_ );
875  h_SimConvTwoMTracksAndVtxPGT0_[4]->Fill( (*mcPho).fourMomentum().et());
876  }
877  if ( chi2Prob > 0.0005) {
878  h_SimConvTwoMTracksAndVtxPGT0005_[0]->Fill( mcEta_ ) ;
879  h_SimConvTwoMTracksAndVtxPGT0005_[1]->Fill( mcPhi_ );
880  h_SimConvTwoMTracksAndVtxPGT0005_[2]->Fill( mcConvR_ );
881  h_SimConvTwoMTracksAndVtxPGT0005_[3]->Fill( mcConvZ_ );
882  h_SimConvTwoMTracksAndVtxPGT0005_[4]->Fill( (*mcPho).fourMomentum().et());
883 
884  }
885 
886 
887  } // loop over reco conversions
888  } //End loop over simulated conversions
889 
890 
891  // ########################### RECO to SIM ############################## //
892 
893  for (reco::ConversionCollection::const_iterator conv = convHandle->begin();conv!=convHandle->end();++conv) {
894  const reco::Conversion aConv = (*conv);
895  if ( mergedTracks_ ) {
897 
898  } else {
900  }
901 
902 
903 
904  std::vector<edm::RefToBase<reco::Track> > tracks = aConv.tracks();
905  const reco::Vertex& vtx = aConv.conversionVertex();
906  //requires two tracks and a valid vertex
907  if (tracks.size() !=2 || !(vtx.isValid())) continue;
908 
909  bool phoIsInBarrel=false;
910  bool phoIsInEndcap=false;
911  RefToBase<reco::Track> tk1 = aConv.tracks().front();
912  RefToBase<reco::Track> tk2 = aConv.tracks().back();
914  tc1.push_back(tk1);
915  tc2.push_back(tk2);
916 
917  //std::cout << " RECO to SIM conversion track pt " << tk1->pt() << " " << tk2->pt() << endl;
918  const reco::Track refTk1 = aConv.conversionVertex().refittedTracks().front();
919  const reco::Track refTk2 = aConv.conversionVertex().refittedTracks().back();
920 
921  //TODO replace it with phi at vertex
922  float dPhiTracksAtVtx = aConv.dPhiTracksAtVtx();
923  // override with the phi calculated at the vertex
924  math::XYZVector p1AtVtx= recalculateMomentumAtFittedVertex ( (*theMF_), *trackerGeom, tk1, aConv.conversionVertex() );
925  math::XYZVector p2AtVtx= recalculateMomentumAtFittedVertex ( (*theMF_), *trackerGeom, tk2, aConv.conversionVertex() );
926  if ( sqrt(p1AtVtx.perp2()) > sqrt(p2AtVtx.perp2()) )
927  dPhiTracksAtVtx = p1AtVtx.phi() - p2AtVtx.phi();
928  else
929  dPhiTracksAtVtx = p2AtVtx.phi() - p1AtVtx.phi();
930 
931 
932  math::XYZVector convMom = tk1->momentum() + tk2->momentum();
933  math::XYZVector refittedMom = aConv.refittedPairMomentum();
934 
935 
936  if (fabs(refittedMom.eta())< 1.479 ) {
937  phoIsInBarrel=true;
938  } else {
939  phoIsInEndcap=true;
940  }
941 
942  nRecConv_++;
943 
945  int match =0;
946  float invM=aConv.pairInvariantMass();
947  float chi2Prob = ChiSquaredProbability( aConv.conversionVertex().chi2(), aConv.conversionVertex().ndof() );
948 
949  h_convEta_[match][0]->Fill( refittedMom.eta() );
950  h_convPhi_[match][0]->Fill( refittedMom.phi() );
951  h_convR_[match][0]->Fill( sqrt(aConv.conversionVertex().position().perp2()) );
952  h_convRplot_->Fill( sqrt(aConv.conversionVertex().position().perp2()) );
953  h_convZ_[match][0]->Fill( aConv.conversionVertex().position().z() );
954  h_convZplot_->Fill( aConv.conversionVertex().position().z() );
955  h_convPt_[match][0]->Fill( sqrt(refittedMom.perp2()) );
956  h_invMass_[match][0] ->Fill( invM);
957  h_vtxChi2Prob_[match][0] ->Fill (chi2Prob);
958 
959 
960 
961  h_distMinAppTracks_[match][0] ->Fill (aConv.distOfMinimumApproach());
962  h_DPhiTracksAtVtx_[match][0]->Fill( dPhiTracksAtVtx);
963  h2_DPhiTracksAtVtxVsEta_->Fill( mcEta_, dPhiTracksAtVtx);
964  h2_DPhiTracksAtVtxVsR_->Fill( mcConvR_, dPhiTracksAtVtx);
965  p_DPhiTracksAtVtxVsEta_->Fill( mcEta_, dPhiTracksAtVtx);
966  p_DPhiTracksAtVtxVsR_->Fill( mcConvR_, dPhiTracksAtVtx);
967 
968  h_DCotTracks_[match][0] ->Fill ( aConv.pairCotThetaSeparation() );
969  h2_DCotTracksVsEta_->Fill( mcEta_, aConv.pairCotThetaSeparation() );
970  h2_DCotTracksVsR_->Fill( mcConvR_, aConv.pairCotThetaSeparation() );
971  p_DCotTracksVsEta_->Fill( mcEta_, aConv.pairCotThetaSeparation() );
972  p_DCotTracksVsR_->Fill( mcConvR_, aConv.pairCotThetaSeparation() );
973 
974  if ( phoIsInBarrel ) {
975  h_invMass_[match][1] ->Fill(invM);
976  h_vtxChi2Prob_[match][1] ->Fill (chi2Prob);
977  h_distMinAppTracks_[match][1] ->Fill (aConv.distOfMinimumApproach());
978  h_DPhiTracksAtVtx_[match][1]->Fill( dPhiTracksAtVtx);
979  h_DCotTracks_[match][1] ->Fill ( aConv.pairCotThetaSeparation() );
980  }
981 
982 
983  if ( phoIsInEndcap ) {
984  h_invMass_[match][2] ->Fill(invM);
985  h_vtxChi2Prob_[match][2] ->Fill (chi2Prob);
986  h_distMinAppTracks_[match][2] ->Fill (aConv.distOfMinimumApproach());
987  h_DPhiTracksAtVtx_[match][2]->Fill( dPhiTracksAtVtx);
988  h_DCotTracks_[match][2] ->Fill ( aConv.pairCotThetaSeparation() );
989  }
990 
991  h_convVtxRvsZ_[0] ->Fill ( fabs (aConv.conversionVertex().position().z() ), sqrt(aConv.conversionVertex().position().perp2()) ) ;
992  h_convVtxYvsX_ ->Fill ( aConv.conversionVertex().position().x(), aConv.conversionVertex().position().y() );
993  h_convVtxYvsX_zoom_[0] ->Fill ( aConv.conversionVertex().position().x(), aConv.conversionVertex().position().y() );
994  h_convVtxYvsX_zoom_[1] ->Fill ( aConv.conversionVertex().position().x(), aConv.conversionVertex().position().y() );
995 
996 
997  // quantities per track: all conversions
998  for (unsigned int i=0; i<tracks.size(); i++) {
999  double d0;
1000  if (valid_pvtx){
1001  d0 = - tracks[i]->dxy(the_pvtx.position());
1002  } else {
1003  d0 = tracks[i]->d0();
1004  }
1005  h_TkD0_[match]->Fill ( d0* tracks[i]->charge() );
1006  nHitsVsEta_[match] ->Fill (mcEta_, float(tracks[i]->numberOfValidHits()) );
1007  nHitsVsR_[match] ->Fill (mcConvR_, float(tracks[i]->numberOfValidHits()) );
1008  p_nHitsVsEta_[match] ->Fill (mcEta_, float(tracks[i]->numberOfValidHits()) -0.0001);
1009  p_nHitsVsR_[match] ->Fill (mcConvR_, float(tracks[i]->numberOfValidHits()) -0.0001);
1010  h_tkChi2_[match] ->Fill (tracks[i]->normalizedChi2() );
1011  h_tkChi2Large_[match] ->Fill (tracks[i]->normalizedChi2() );
1012  h2_Chi2VsEta_[match] ->Fill( mcEta_, tracks[i]->normalizedChi2() );
1013  h2_Chi2VsR_[match] ->Fill( mcConvR_, tracks[i]->normalizedChi2() );
1014  p_Chi2VsEta_[match] ->Fill( mcEta_, tracks[i]->normalizedChi2() );
1015  p_Chi2VsR_[match] ->Fill( mcConvR_, tracks[i]->normalizedChi2() );
1016 
1017  }
1018 
1019  bool associated = false;
1020  float mcConvPt_= -99999999;
1021  // float mcPhi= 0; // unused
1022  float simPV_Z=0;
1023  for ( std::vector<PhotonMCTruth>::const_iterator mcPho=mcPhotons.begin(); mcPho !=mcPhotons.end(); mcPho++) {
1024  mcConvPt_= (*mcPho).fourMomentum().et();
1025  float mcPhi= (*mcPho).fourMomentum().phi();
1026  simPV_Z = (*mcPho).primaryVertex().z();
1027  mcPhi_= phiNormalization(mcPhi);
1028  mcEta_= (*mcPho).fourMomentum().pseudoRapidity();
1029  mcEta_ = etaTransformation(mcEta_, (*mcPho).primaryVertex().z() );
1030  mcConvR_= (*mcPho).vertex().perp();
1031  mcConvX_= (*mcPho).vertex().x();
1032  mcConvY_= (*mcPho).vertex().y();
1033  mcConvZ_= (*mcPho).vertex().z();
1034  mcConvEta_= (*mcPho).vertex().eta();
1035  mcConvPhi_= (*mcPho).vertex().phi();
1036  if ( fabs(mcEta_) > END_HI ) continue;
1037  if (mcConvPt_<minPhoPtForPurity) continue;
1038  if (fabs(mcEta_)>maxPhoEtaForPurity) continue;
1039  if (fabs(mcConvZ_)>maxPhoZForPurity) continue;
1040  if (mcConvR_>maxPhoRForEffic) continue;
1041 
1042  if ( (*mcPho).isAConversion() != 1 ) continue;
1043  if (!( ( fabs(mcEta_) <= BARL && mcConvR_ <85 ) ||
1044  ( fabs(mcEta_) > BARL && fabs(mcEta_) <=END_HI && fabs( (*mcPho).vertex().z() ) < 210 ) ) )
1045  continue;
1046 
1047 
1048  theConvTP_.clear();
1049  for(size_t i = 0; i < tpForFakeRate.size(); ++i){
1050  TrackingParticleRef tp (TPHandleForFakeRate,i);
1051  if ( fabs( tp->vx() - (*mcPho).vertex().x() ) < 0.0001 &&
1052  fabs( tp->vy() - (*mcPho).vertex().y() ) < 0.0001 &&
1053  fabs( tp->vz() - (*mcPho).vertex().z() ) < 0.0001) {
1054  theConvTP_.push_back( tp );
1055 
1056 
1057  }
1058  }
1059 
1060  if ( theConvTP_.size() < 2 ) continue;
1061 
1062  associated = false;
1063  reco::RecoToSimCollection p1 = theTrackAssociator_->associateRecoToSim(tc1,theConvTP_,&e);
1064  reco::RecoToSimCollection p2 = theTrackAssociator_->associateRecoToSim(tc2,theConvTP_,&e);
1065  try{
1066  std::vector<std::pair<TrackingParticleRef, double> > tp1 = p1[tk1];
1067  std::vector<std::pair<TrackingParticleRef, double> > tp2 = p2[tk2];
1068  if (!(tp1.size()&&tp2.size())){
1069  tp1 = p1[tk2];
1070  tp2 = p2[tk1];
1071  }
1072  if (tp1.size()&&tp2.size()) {
1073  TrackingParticleRef tpr1 = tp1.front().first;
1074  TrackingParticleRef tpr2 = tp2.front().first;
1075  if (abs(tpr1->pdgId())==11&&abs(tpr2->pdgId())==11) {
1076  if ( (tpr1->parentVertex()->sourceTracks_end()-tpr1->parentVertex()->sourceTracks_begin()==1) &&
1077  (tpr2->parentVertex()->sourceTracks_end()-tpr2->parentVertex()->sourceTracks_begin()==1)) {
1078  if (tpr1->parentVertex().key()==tpr2->parentVertex().key() && ((*tpr1->parentVertex()->sourceTracks_begin())->pdgId()==22)) {
1079  mcConvR_ = sqrt(tpr1->parentVertex()->position().Perp2());
1080  mcConvZ_ = tpr1->parentVertex()->position().z();
1081  mcConvX_ = tpr1->parentVertex()->position().x();
1082  mcConvY_ = tpr1->parentVertex()->position().y();
1083  mcConvEta_ = tpr1->parentVertex()->position().eta();
1084  mcConvPhi_ = tpr1->parentVertex()->position().phi();
1085  mcConvPt_ = sqrt((*tpr1->parentVertex()->sourceTracks_begin())->momentum().Perp2());
1086  //std::cout << " Reco to Sim mcconvpt " << mcConvPt_ << std::endl;
1087  //cout << "associated track1 to " << tpr1->pdgId() << " with p=" << tpr1->p4() << " with pT=" << tpr1->pt() << endl;
1088  //cout << "associated track2 to " << tpr2->pdgId() << " with p=" << tpr2->p4() << " with pT=" << tpr2->pt() << endl;
1089  associated = true;
1090  break;
1091  }
1092  }
1093  }
1094  }
1095  } catch (Exception event) {
1096  //cout << "do not continue: " << event.what() << endl;
1097  //continue;
1098  }
1099 
1100  }// end loop on sim photons
1101 
1102  if ( associated ) match=1;
1103  else
1104  match=2;
1105 
1106  h_match_->Fill(float(match));
1108  if ( match == 1) nRecConvAss_++;
1109  h_convEta_[match][0]->Fill( refittedMom.eta() );
1110  h_convPhi_[match][0]->Fill( refittedMom.phi() );
1111  h_convR_[match][0]->Fill( sqrt(aConv.conversionVertex().position().perp2()) );
1112  h_convZ_[match][0]->Fill( aConv.conversionVertex().position().z() );
1113  h_convPt_[match][0]->Fill( sqrt(refittedMom.perp2()) );
1114  h_invMass_[match][0] ->Fill( invM);
1115  h_vtxChi2Prob_[match][0] ->Fill (chi2Prob);
1116  h_DPhiTracksAtVtx_[match][0]->Fill( dPhiTracksAtVtx);
1117  h_DCotTracks_[match][0] ->Fill ( aConv.pairCotThetaSeparation() );
1118  h_distMinAppTracks_[match][0] ->Fill (aConv.distOfMinimumApproach());
1119 
1120  if ( match==1) {
1121  h2_photonPtRecVsPtSim_->Fill ( mcConvPt_, sqrt(refittedMom.perp2()) );
1122  h_convPtRes_[0]->Fill ( sqrt(refittedMom.perp2())/mcConvPt_);
1123  }
1124 
1125  if ( phoIsInBarrel ) {
1126  h_invMass_[match][1] ->Fill(invM);
1127  h_vtxChi2Prob_[match][1] ->Fill (chi2Prob);
1128  h_DPhiTracksAtVtx_[match][1]->Fill( dPhiTracksAtVtx);
1129  h_DCotTracks_[match][1] ->Fill ( aConv.pairCotThetaSeparation() );
1130  h_distMinAppTracks_[match][1] ->Fill (aConv.distOfMinimumApproach());
1131  if ( match==1) h_convPtRes_[1]->Fill ( sqrt(refittedMom.perp2())/mcConvPt_);
1132  }
1133 
1134 
1135  if ( phoIsInEndcap ) {
1136  h_invMass_[match][2] ->Fill(invM);
1137  h_vtxChi2Prob_[match][2] ->Fill (chi2Prob);
1138  h_DPhiTracksAtVtx_[match][2]->Fill( dPhiTracksAtVtx);
1139  h_DCotTracks_[match][2] ->Fill ( aConv.pairCotThetaSeparation() );
1140  h_distMinAppTracks_[match][2] ->Fill (aConv.distOfMinimumApproach());
1141  if ( match==1) h_convPtRes_[2]->Fill ( sqrt(refittedMom.perp2())/mcConvPt_);
1142  }
1143 
1144 
1145  if ( match == 1 ) {
1146  h_convVtxdX_ ->Fill ( aConv.conversionVertex().position().x() - mcConvX_);
1147  h_convVtxdY_ ->Fill ( aConv.conversionVertex().position().y() - mcConvY_);
1148  h_convVtxdZ_ ->Fill ( aConv.conversionVertex().position().z() - mcConvZ_);
1149  h_convVtxdR_ ->Fill ( sqrt(aConv.conversionVertex().position().perp2()) - mcConvR_);
1150  h_convVtxdPhi_ ->Fill ( aConv.conversionVertex().position().phi() - mcConvPhi_);
1151  h_convVtxdEta_ ->Fill ( aConv.conversionVertex().position().eta() - mcConvEta_);
1152  h2_convVtxdRVsR_ ->Fill (mcConvR_, sqrt(aConv.conversionVertex().position().perp2()) - mcConvR_ );
1153  h2_convVtxdRVsEta_ ->Fill (mcEta_, sqrt(aConv.conversionVertex().position().perp2()) - mcConvR_ );
1154  p_convVtxdRVsR_ ->Fill (mcConvR_, sqrt(aConv.conversionVertex().position().perp2()) - mcConvR_ );
1155  p_convVtxdRVsEta_ ->Fill (mcEta_, sqrt(aConv.conversionVertex().position().perp2()) - mcConvR_ );
1156  p_convVtxdXVsX_ ->Fill (mcConvX_, aConv.conversionVertex().position().x() - mcConvX_ );
1157  p_convVtxdYVsY_ ->Fill (mcConvY_, aConv.conversionVertex().position().y() - mcConvY_ );
1158  p_convVtxdZVsZ_ ->Fill (mcConvZ_, aConv.conversionVertex().position().z() - mcConvZ_ );
1159  p_convVtxdZVsR_ ->Fill (mcConvR_, aConv.conversionVertex().position().z() - mcConvZ_ );
1160 
1161  float dR=sqrt(aConv.conversionVertex().position().perp2()) - mcConvR_;
1162  float dZ=aConv.conversionVertex().position().z() - mcConvZ_;
1163  p2_convVtxdRVsRZ_ ->Fill (mcConvZ_,mcConvR_, dR );
1164  p2_convVtxdZVsRZ_ ->Fill (mcConvZ_,mcConvR_, dZ );
1165 
1166 
1167 
1168 
1169  h2_convVtxRrecVsTrue_ -> Fill (mcConvR_, sqrt(aConv.conversionVertex().position().perp2()) );
1170 
1171 
1172  h_zPVFromTracks_[match]->Fill ( aConv.zOfPrimaryVertexFromTracks() );
1173  h_dzPVFromTracks_[match]->Fill ( aConv.zOfPrimaryVertexFromTracks() - simPV_Z );
1174  h2_dzPVVsR_ ->Fill(mcConvR_, aConv.zOfPrimaryVertexFromTracks() - simPV_Z );
1175  p_dzPVVsR_ ->Fill(mcConvR_, aConv.zOfPrimaryVertexFromTracks() - simPV_Z );
1176 
1177  if ( phoIsInBarrel ) {
1178  h_convVtxdX_barrel_ ->Fill ( aConv.conversionVertex().position().x() - mcConvX_);
1179  h_convVtxdY_barrel_ ->Fill ( aConv.conversionVertex().position().y() - mcConvY_);
1180  h_convVtxdZ_barrel_ ->Fill ( aConv.conversionVertex().position().z() - mcConvZ_);
1181  h_convVtxdR_barrel_ ->Fill ( sqrt(aConv.conversionVertex().position().perp2()) - mcConvR_);
1182 
1183  }
1184  if ( phoIsInEndcap ) {
1185  h_convVtxdX_endcap_ ->Fill ( aConv.conversionVertex().position().x() - mcConvX_);
1186  h_convVtxdY_endcap_ ->Fill ( aConv.conversionVertex().position().y() - mcConvY_);
1187  h_convVtxdZ_endcap_ ->Fill ( aConv.conversionVertex().position().z() - mcConvZ_);
1188  h_convVtxdR_endcap_ ->Fill ( sqrt(aConv.conversionVertex().position().perp2()) - mcConvR_);
1189 
1190  }
1191 
1192 
1193  }
1194 
1196  for (unsigned int i=0; i<tracks.size(); i++) {
1197  //std::cout << " Loop over tracks pt " << tracks[i]->pt() << std::endl;
1198  RefToBase<reco::Track> tfrb(aConv.tracks()[i] );
1199  itAss= myAss.find( tfrb.get() );
1200 
1201  nHitsVsEta_[match] ->Fill (mcEta_, float(tracks[i]->numberOfValidHits()) );
1202  nHitsVsR_[match] ->Fill (mcConvR_, float(tracks[i]->numberOfValidHits()) );
1203  p_nHitsVsEta_[match] ->Fill (mcEta_, float(tracks[i]->numberOfValidHits()) -0.0001);
1204  p_nHitsVsR_[match] ->Fill (mcConvR_, float(tracks[i]->numberOfValidHits()) -0.0001);
1205  h_tkChi2_[match] ->Fill (tracks[i]->normalizedChi2() );
1206  h_tkChi2Large_[match] ->Fill (tracks[i]->normalizedChi2() );
1207  h2_Chi2VsEta_[match] ->Fill( mcEta_, tracks[i]->normalizedChi2() );
1208  h2_Chi2VsR_[match] ->Fill( mcConvR_, tracks[i]->normalizedChi2() );
1209  p_Chi2VsEta_[match] ->Fill( mcEta_, tracks[i]->normalizedChi2() );
1210  p_Chi2VsR_[match] ->Fill( mcConvR_, tracks[i]->normalizedChi2() );
1211  double d0;
1212  if (valid_pvtx){
1213  d0 = - tracks[i]->dxy(the_pvtx.position());
1214  } else {
1215  d0 = tracks[i]->d0();
1216  }
1217  h_TkD0_[match]->Fill (d0* tracks[i]->charge() );
1218 
1219 
1220  if ( itAss == myAss.end() ) continue;
1221  reco::Track refTrack= aConv.conversionVertex().refittedTracks()[i];
1222 
1223  float simPt = sqrt( ((*itAss).second)->momentum().perp2() );
1224  float recPt = refTrack.pt();
1225  float ptres= recPt - simPt ;
1226  //float pterror = aConv.tracks()[i]->ptError();
1227  float pterror = aConv.conversionVertex().refittedTracks()[i].ptError();
1228  h2_PtRecVsPtSim_[0]->Fill ( simPt, recPt);
1229  h_TkPtPull_[0] ->Fill(ptres/pterror);
1230  h2_TkPtPull_[0] ->Fill(mcEta_, ptres/pterror);
1231 
1232  if ( phoIsInBarrel ) {
1233  h_TkPtPull_[1] ->Fill(ptres/pterror);
1234  h2_PtRecVsPtSim_[1]->Fill ( simPt, recPt);
1235  }
1236  if ( phoIsInEndcap ) {
1237  h_TkPtPull_[2] ->Fill(ptres/pterror);
1238  h2_PtRecVsPtSim_[2]->Fill ( simPt, recPt);
1239  }
1240  } // end loop over track
1241 
1242 
1243 
1244  } // loop over reco conversions
1245 
1246 
1247  h_nConv_[0][0]->Fill (float(nRecConv_));
1248  h_nConv_[1][0]->Fill (float(nRecConvAss_));
1249 
1250 
1251 
1252 }
1253 
1254 
1255 
1256 
1257 
1259 
1260 
1261  std::string outputFileName = parameters_.getParameter<std::string>("OutputFileName");
1262  if ( ! isRunCentrally_ ) {
1263  dbe_->save(outputFileName);
1264  }
1265 
1266  edm::LogInfo("TkConvValidator") << "Analyzed " << nEvt_ << "\n";
1267  // std::cout << "::endJob Analyzed " << nEvt_ << " events " << " with total " << nPho_ << " Photons " << "\n";
1268  // std::cout << "TkConvValidator::endJob Analyzed " << nEvt_ << " events " << "\n";
1269 
1270  return ;
1271 }
1272 
1273 
1275 
1279  SimpleCylinderBounds( sqrt(vtx.position().perp2())-0.001,
1280  sqrt(vtx.position().perp2())+0.001,
1281  -fabs(vtx.position().z()),
1282  fabs(vtx.position().z()))));
1283 
1284  ReferenceCountingPointer<BoundDisk> theDisk_(new BoundDisk( Surface::PositionType( 0, 0, vtx.position().z()), rot,
1285  SimpleDiskBounds( 0, sqrt(vtx.position().perp2()), -0.001, 0.001)));
1286 
1287  TrajectoryStateTransform transformer;
1288  const TrajectoryStateOnSurface myTSOS = transformer.innerStateOnSurface(*tk, trackerGeom, &mf);
1289  PropagatorWithMaterial propag( anyDirection, 0.000511, &mf );
1290  TrajectoryStateOnSurface stateAtVtx;
1291  stateAtVtx = propag.propagate(myTSOS, *theBarrel_);
1292  if (!stateAtVtx.isValid() ) {
1293  stateAtVtx = propag.propagate(myTSOS, *theDisk_);
1294  }
1295  if (stateAtVtx.isValid()){
1296  return math::XYZVector ( double(stateAtVtx.globalMomentum().x()), double(stateAtVtx.globalMomentum().y()), double(stateAtVtx.globalMomentum().z()));
1297  } else {
1298  return math::XYZVector(0.,0.,0.);
1299  }
1300 
1301 
1302 
1303 }
1304 
1305 
1307 {
1308  //---Definitions
1309  const float PI = 3.1415927;
1310  const float TWOPI = 2.0*PI;
1311 
1312 
1313  if(phi > PI) {phi = phi - TWOPI;}
1314  if(phi < -PI) {phi = phi + TWOPI;}
1315 
1316  // cout << " Float_t PHInormalization out " << PHI << endl;
1317  return phi;
1318 
1319 }
1320 
1321 
1322 float TkConvValidator::etaTransformation( float EtaParticle , float Zvertex) {
1323 
1324  //---Definitions
1325  const float PI = 3.1415927;
1326 
1327  //---Definitions for ECAL
1328  const float R_ECAL = 136.5;
1329  const float Z_Endcap = 328.0;
1330  const float etaBarrelEndcap = 1.479;
1331 
1332  //---ETA correction
1333 
1334  float Theta = 0.0 ;
1335  float ZEcal = R_ECAL*sinh(EtaParticle)+Zvertex;
1336 
1337  if(ZEcal != 0.0) Theta = atan(R_ECAL/ZEcal);
1338  if(Theta<0.0) Theta = Theta+PI ;
1339  float ETA = - log(tan(0.5*Theta));
1340 
1341  if( fabs(ETA) > etaBarrelEndcap )
1342  {
1343  float Zend = Z_Endcap ;
1344  if(EtaParticle<0.0 ) Zend = -Zend ;
1345  float Zlen = Zend - Zvertex ;
1346  float RR = Zlen/sinh(EtaParticle);
1347  Theta = atan(RR/Zend);
1348  if(Theta<0.0) Theta = Theta+PI ;
1349  ETA = - log(tan(0.5*Theta));
1350  }
1351  //---Return the result
1352  return ETA;
1353  //---end
1354 }
1355 
1356 
const reco::Vertex & conversionVertex() const
returns the reco conversion vertex
Definition: Conversion.h:100
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
virtual void endRun(edm::Run &r, edm::EventSetup const &es)
const Vector & momentum() const
track momentum vector
Definition: TrackBase.h:149
math::XYZVector refittedPairMomentum() const
Conversion tracks momentum from the tracks refitted with vertex constraint.
Definition: Conversion.cc:277
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:519
#define PI
const_iterator end() const
last iterator over the map (read only)
tuple d0
Definition: debug_cff.py:3
static HepMC::IO_HEPEVT conv
bool quality(ConversionQuality q) const
Definition: Conversion.h:168
math::XYZVector recalculateMomentumAtFittedVertex(const MagneticField &mf, const TrackerGeometry &trackerGeom, const edm::RefToBase< reco::Track > &tk, const reco::Vertex &vtx)
bool isValid() const
Tells whether the vertex is valid.
Definition: Vertex.h:61
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE")
Definition: DQMStore.cc:1883
std::vector< GenJet > GenJetCollection
collection of GenJet objects
const_iterator find(const key_type &k) const
find element with specified reference key
T y() const
Definition: PV3DBase.h:57
virtual void analyze(const edm::Event &, const edm::EventSetup &)
double distOfMinimumApproach() const
Definition: Conversion.h:130
#define abs(x)
Definition: mlp_lapack.h:159
double pairCotThetaSeparation() const
Delta cot(Theta) where Theta is the angle in the (y,z) plane between the two tracks. Original tracks are used.
Definition: Conversion.cc:238
#define TWOPI
EgammaCoreTools.
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
const std::vector< Track > & refittedTracks() const
Returns the container of refitted tracks.
Definition: Vertex.h:136
double pairInvariantMass() const
if nTracks=2 returns the pair invariant mass. Original tracks are used here
Definition: Conversion.cc:219
const Point & position() const
position
Definition: Vertex.h:93
double charge(const std::vector< uint8_t > &Ampls)
double q2[4]
Definition: TauolaWrapper.h:88
std::vector< Conversion > ConversionCollection
collectin of Conversion objects
Definition: ConversionFwd.h:9
return((rh^lh)&mask)
virtual ~TkConvValidator()
std::vector< edm::RefToBase< reco::Track > > tracks() const
vector of track to base references
Definition: Conversion.cc:180
double zOfPrimaryVertexFromTracks() const
z coordinate of the photon origin calculated from the track-pair direction and the position of the ve...
Definition: Conversion.cc:205
#define ETA
float etaTransformation(float a, float b)
float phiNormalization(float &a)
T sqrt(T t)
Definition: SSEVec.h:28
double pt() const
track transverse momentum
Definition: TrackBase.h:130
T z() const
Definition: PV3DBase.h:58
tuple result
Definition: query.py:137
TkConvValidator(const edm::ParameterSet &)
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
double chi2() const
chi-squares
Definition: Vertex.h:82
float ChiSquaredProbability(double chiSquared, double nrDOF)
tuple pset
Definition: CrabTask.py:85
MonitorElement * bookProfile(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, const char *option="s")
Definition: DQMStore.cc:833
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
double p2[4]
Definition: TauolaWrapper.h:90
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:359
DQMStore * dbe_
double ndof() const
Definition: Vertex.h:89
virtual TrajectoryStateOnSurface propagate(const TrajectoryStateOnSurface &tsos, const Plane &plane) const
Log< T >::type log(const T &t)
Definition: Log.h:22
double q1[4]
Definition: TauolaWrapper.h:87
tuple tracks
Definition: testEve_cfg.py:39
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
const T & get() const
Definition: EventSetup.h:55
key_type key() const
Accessor for product key.
Definition: Ref.h:264
T const * product() const
Definition: ESHandle.h:62
virtual void endJob()
static const float etaBarrelEndcap
TrajectoryStateOnSurface innerStateOnSurface(const reco::Track &tk, const TrackingGeometry &geom, const MagneticField *field) const
static const float Z_Endcap
edm::EventID id() const
Definition: EventBase.h:56
double p1[4]
Definition: TauolaWrapper.h:89
GlobalVector globalMomentum() const
void push_back(const RefToBase< T > &)
static const float R_ECAL
double recPt
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:6
virtual void beginJob()
std::vector< TrackingParticle > TrackingParticleCollection
virtual void beginRun(edm::Run const &r, edm::EventSetup const &theEventSetup)
MonitorElement * book2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY)
Book 2D histogram.
Definition: DQMStore.cc:647
T x() const
Definition: PV3DBase.h:56
value_type const * get() const
Definition: RefToBase.h:207
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:237
double dPhiTracksAtVtx() const
Definition: Conversion.cc:343
Definition: Run.h:31
Definition: DDAxes.h:10
MonitorElement * bookProfile2D(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, int nchZ, double lowZ, double highZ, const char *option="s")
Definition: DQMStore.cc:977