ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/Jeng/PVStudy/plugins/PVStudy.cc

Comparing UserCode/Jeng/PVStudy/plugins/PVStudy.cc (file contents):
Revision 1.17 by yygao, Tue Jan 26 12:06:49 2010 UTC vs.
Revision 1.23 by yygao, Mon May 31 14:19:45 2010 UTC

#	Line 38 | Line 38 | Implementation:
38		#include "DataFormats/TrackReco/interface/Track.h"
39		#include "DataFormats/TrackReco/interface/TrackFwd.h"
40		#include "FWCore/ServiceRegistry/interface/Service.h"
41	<	#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
41	>	#include "CommonTools/UtilAlgos/interface/TFileService.h"
42		#include "UserCode/PVStudy/interface/PVStudy.h"
43		//
44		#include "DataFormats/VertexReco/interface/Vertex.h"
#	Line 51 | Line 51 | Implementation:
51		#include <SimDataFormats/Track/interface/SimTrackContainer.h>
52		// BeamSpot
53		#include "DataFormats/BeamSpot/interface/BeamSpot.h"
54	–	// Trigger
55	–	#include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
56	–	#include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
57	–	#include "CondFormats/L1TObjects/interface/L1GtPrescaleFactors.h"
58	–	#include "CondFormats/DataRecord/interface/L1GtPrescaleFactorsAlgoTrigRcd.h"
59	–	#include "CondFormats/DataRecord/interface/L1GtPrescaleFactorsTechTrigRcd.h"
60	–	#include "CondFormats/L1TObjects/interface/L1GtTriggerMask.h"
61	–	#include "CondFormats/DataRecord/interface/L1GtTriggerMaskAlgoTrigRcd.h"
62	–	#include "CondFormats/DataRecord/interface/L1GtTriggerMaskTechTrigRcd.h"
63	–	#include "CondFormats/DataRecord/interface/L1GtTriggerMaskVetoAlgoTrigRcd.h"
64	–	#include "CondFormats/DataRecord/interface/L1GtTriggerMaskVetoTechTrigRcd.h"
65	–	#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
66	–	#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
67	–	#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
68	–
54
55		//root
56		#include <TROOT.h>
#	Line 95 | Line 80 | PVStudy::PVStudy(const edm::ParameterSet
80		realData_                  = iConfig.getUntrackedParameter<bool>("realData",false);
81		analyze_                   = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
82		saventuple_                = iConfig.getUntrackedParameter<bool>("saventuple",false);
83	<	outputfilename_            = iConfig.getUntrackedParameter<string>("OutputFileName");
83	>	outputfilename_            = iConfig.getUntrackedParameter<string>("OutputFileName");
84	>	histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
85		ntrkdiffcut_               = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
86		nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
87		nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
88		zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
89	+	useHWTrk_                  = iConfig.getUntrackedParameter<bool>("useHWTrk",false);
90	+	avgTrkPtInPVMin_                  = iConfig.getUntrackedParameter<double>("avgTrkPtInPVMin");
91	+	avgTrkPtInPVMax_                  = iConfig.getUntrackedParameter<double>("avgTrkPtInPVMax");
92		bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
104	–	applyEvtClean_            = iConfig.getUntrackedParameter<bool>("applyEvtClean",false);
105	–	histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
93
94		//now do what ever initialization is needed
95		pvinfo.clear();
96
110	–
97		// Specify the data mode vector
98		if(realData_) datamode.push_back(0);
99		else {
#	Line 127 | Line 113 | PVStudy::PVStudy(const edm::ParameterSet
113
114		// pvtxtree_ analyzing the pvtxs ootb
115		pvtxtree_ = new TTree("pvtxtree","pvtxtree");
130	–	//pvtx using all tracks
131	–	pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
132	–	pvtxtree_->Branch("nTrkPV",&nTrkPV_,"nTrkPV[nrecPV]/I");
133	–	pvtxtree_->Branch("sumptsq",&sumptsq_,"sumptsq[nrecPV]/D");
134	–	pvtxtree_->Branch("isValid",&isValid_,"isValid/I");
135	–	pvtxtree_->Branch("isFake",&isFake_,"isFake/I");
136	–	pvtxtree_->Branch("recx",&recx_,"recx[nrecPV]/D");
137	–	pvtxtree_->Branch("recy",&recy_,"recy[nrecPV]/D");
138	–	pvtxtree_->Branch("recz",&recz_,"recz[nrecPV]/D");
139	–	pvtxtree_->Branch("recx_err",&recx_err_,"recx_err[nrecPV]/D");
140	–	pvtxtree_->Branch("recy_err",&recy_err_,"recy_err[nrecPV]/D");
141	–	pvtxtree_->Branch("recz_err",&recz_err_,"recz_err[nrecPV]/D");
116
117	+	// Event information for the data
118	+	pvtxtree_->Branch("glob_runno",&glob_runno_,"glob_runno/I");
119	+	pvtxtree_->Branch("glob_evtno",&glob_evtno_,"glob_evtno/I");
120	+	pvtxtree_->Branch("glob_ls",&glob_ls_,"glob_ls/I");
121	+	pvtxtree_->Branch("glob_bx",&glob_bx_,"glob_bx/I");
122	+
123	+
124	+	pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
125	+	pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
126	+	pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
127	+	// Event level information
128		pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
129		pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
130		pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
131
132	<	//pvtx using splittracks1
133	<	pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
149	<	pvtxtree_->Branch("nTrkPV1_spl",&nTrkPV1_spl_,"nTrkPV1_spl[nrecPV1_spl]/I");
150	<	pvtxtree_->Branch("sumptsq1_spl",&sumptsq1_spl_,"sumptsq1_spl[nrecPV1_spl]/D");
151	<	pvtxtree_->Branch("isValid1_spl",&isValid1_spl_,"isValid1_spl/I");
152	<	pvtxtree_->Branch("isFake1_spl",&isFake1_spl_,"isFake1_spl/I");
153	<	pvtxtree_->Branch("recx1_spl",&recx1_spl_,"recx1_spl[nrecPV1_spl]/D");
154	<	pvtxtree_->Branch("recy1_spl",&recy1_spl_,"recy1_spl[nrecPV1_spl]/D");
155	<	pvtxtree_->Branch("recz1_spl",&recz1_spl_,"recz1_spl[nrecPV1_spl]/D");
156	<	pvtxtree_->Branch("recx1_err_spl",&recx1_err_spl_,"recx1_err_spl[nrecPV1_spl]/D");
157	<	pvtxtree_->Branch("recy1_err_spl",&recy1_err_spl_,"recy1_err_spl[nrecPV1_spl]/D");
158	<	pvtxtree_->Branch("recz1_err_spl",&recz1_err_spl_,"recz1_err_spl[nrecPV1_spl]/D");
159	<
160	<	//pvtx using splittracks2
161	<	pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
162	<	pvtxtree_->Branch("nTrkPV2_spl",&nTrkPV2_spl_,"nTrkPV2_spl[nrecPV2_spl]/I");
163	<	pvtxtree_->Branch("sumptsq2_spl",&sumptsq2_spl_,"sumptsq2_spl[nrecPV2_spl]/D");
164	<	pvtxtree_->Branch("isValid2_spl",&isValid2_spl_,"isValid2_spl/I");
165	<	pvtxtree_->Branch("isFake2_spl",&isFake2_spl_,"isFake2_spl/I");
166	<	pvtxtree_->Branch("recx2_spl",&recx2_spl_,"recx2_spl[nrecPV2_spl]/D");
167	<	pvtxtree_->Branch("recy2_spl",&recy2_spl_,"recy2_spl[nrecPV2_spl]/D");
168	<	pvtxtree_->Branch("recz2_spl",&recz2_spl_,"recz2_spl[nrecPV2_spl]/D");
169	<	pvtxtree_->Branch("recx2_err_spl",&recx2_err_spl_,"recx2_err_spl[nrecPV2_spl]/D");
170	<	pvtxtree_->Branch("recy2_err_spl",&recy2_err_spl_,"recy2_err_spl[nrecPV2_spl]/D");
171	<	pvtxtree_->Branch("recz2_err_spl",&recz2_err_spl_,"recz2_err_spl[nrecPV2_spl]/D");
172	<
173	<	//pixeVertices
174	<	pvtxtree_->Branch("nrecPV_pxlpvtx",&nrecPV_pxlpvtx_,"nrecPV_pxlpvtx/I");
175	<	pvtxtree_->Branch("nTrkPV_pxlpvtx",&nTrkPV_pxlpvtx_,"nTrkPV_pxlpvtx[nrecPV_pxlpvtx]/I");
176	<	pvtxtree_->Branch("sumptsq_pxlpvtx",&sumptsq_pxlpvtx_,"sumptsq_pxlpvtx[nrecPV_pxlpvtx]/D");
177	<	pvtxtree_->Branch("isValid_pxlpvtx",&isValid_pxlpvtx_,"isValid_pxlpvtx/I");
178	<	pvtxtree_->Branch("isFake_pxlpvtx",&isFake_pxlpvtx_,"isFake_pxlpvtx/I");
179	<	pvtxtree_->Branch("recx_pxlpvtx",&recx_pxlpvtx_,"recx_pxlpvtx[nrecPV_pxlpvtx]/D");
180	<	pvtxtree_->Branch("recy_pxlpvtx",&recy_pxlpvtx_,"recy_pxlpvtx[nrecPV_pxlpvtx]/D");
181	<	pvtxtree_->Branch("recz_pxlpvtx",&recz_pxlpvtx_,"recz_pxlpvtx[nrecPV_pxlpvtx]/D");
182	<	pvtxtree_->Branch("recx_err_pxlpvtx",&recx_err_pxlpvtx_,"recx_err_pxlpvtx[nrecPV_pxlpvtx]/D");
183	<	pvtxtree_->Branch("recy_err_pxlpvtx",&recy_err_pxlpvtx_,"recy_err_pxlpvtx[nrecPV_pxlpvtx]/D");
184	<	pvtxtree_->Branch("recz_err_pxlpvtx",&recz_err_pxlpvtx_,"recz_err_pxlpvtx[nrecPV_pxlpvtx]/D");
132	>	//Fill the variables in the twovtx pair (recvtx1, recvtx2)
133	>	// Information related to the analyzing the two-vertex method
134
186	–	//Fill the variables in the twovtx pair (recvtx1, recvtx2)
135		pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
188	–	pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
189	–	pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
190	–	pvtxtree_->Branch("sumptsq1_spl_twovtx",&sumptsq1_spl_twovtx_,"sumptsq1_spl_twovtx[nrecPV_twovtx]/D");
191	–	pvtxtree_->Branch("sumptsq2_spl_twovtx",&sumptsq2_spl_twovtx_,"sumptsq2_spl_twovtx[nrecPV_twovtx]/D");
136		pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
137		pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
138		pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
#	Line 199 | Line 143 | PVStudy::PVStudy(const edm::ParameterSet
143		pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
144		pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
145		pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
146	+
147	+	// Information for the splitVertexColl1
148	+	pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
149	+	pvtxtree_->Branch("ndofPV1_spl_twovtx",&ndofPV1_spl_twovtx_,"ndofPV1_spl_twovtx[nrecPV_twovtx]/D");
150	+	pvtxtree_->Branch("normchi2PV1_spl_twovtx",&normchi2PV1_spl_twovtx_,"normchi2PV1_spl_twovtx[nrecPV_twovtx]/D");
151	+	pvtxtree_->Branch("avgPtPV1_spl_twovtx",&avgPtPV1_spl_twovtx_,"avgPtPV1_spl_twovtx[nrecPV_twovtx]/D");
152	+	pvtxtree_->Branch("errx1_spl_twovtx",&errx1_spl_twovtx_,"errx1_spl_twovtx[nrecPV_twovtx]/D");
153	+	pvtxtree_->Branch("erry1_spl_twovtx",&erry1_spl_twovtx_,"erry1_spl_twovtx[nrecPV_twovtx]/D");
154	+	pvtxtree_->Branch("errz1_spl_twovtx",&errz1_spl_twovtx_,"errz1_spl_twovtx[nrecPV_twovtx]/D");
155	+
156	+
157	+	// Information for the splitVertexColl2
158	+	pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
159	+	pvtxtree_->Branch("ndofPV2_spl_twovtx",&ndofPV2_spl_twovtx_,"ndofPV2_spl_twovtx[nrecPV_twovtx]/D");
160	+	pvtxtree_->Branch("normchi2PV2_spl_twovtx",&normchi2PV2_spl_twovtx_,"normchi2PV2_spl_twovtx[nrecPV_twovtx]/D");
161	+	pvtxtree_->Branch("avgPtPV2_spl_twovtx",&avgPtPV2_spl_twovtx_,"avgPtPV2_spl_twovtx[nrecPV_twovtx]/D");
162	+	pvtxtree_->Branch("errx2_spl_twovtx",&errx2_spl_twovtx_,"errx2_spl_twovtx[nrecPV_twovtx]/D");
163	+	pvtxtree_->Branch("erry2_spl_twovtx",&erry2_spl_twovtx_,"erry2_spl_twovtx[nrecPV_twovtx]/D");
164	+	pvtxtree_->Branch("errz2_spl_twovtx",&errz2_spl_twovtx_,"errz2_spl_twovtx[nrecPV_twovtx]/D");
165	+
166
167		// MC variables
168		if(!realData_) {
#	Line 216 | Line 180 | PVStudy::PVStudy(const edm::ParameterSet
180		pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
181		pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
182		pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
183	<	pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
183	>	pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
184	>	pvtxtree_->Branch("errx_mct",&errx_mct_,"errx_mct[nrecPV_mct]/D");
185	>	pvtxtree_->Branch("erry_mct",&erry_mct_,"erry_mct[nrecPV_mct]/D");
186	>	pvtxtree_->Branch("errz_mct",&errz_mct_,"errz_mct[nrecPV_mct]/D");
187	>	pvtxtree_->Branch("nTrkPV_mct",&nTrkPV_mct_,"nTrkPV_mct[nrecPV_mct]/I");
188	>	pvtxtree_->Branch("ndofPV_mct",&ndofPV_mct_,"ndofPV_mct[nrecPV_mct]/D");
189	>	pvtxtree_->Branch("normchi2PV_mct",&normchi2PV_mct_,"normchi2PV_mct[nrecPV_mct]/D");
190	>	pvtxtree_->Branch("avgPtPV_mct",&avgPtPV_mct_,"avgPtPV_mct[nrecPV_mct]/D");
191	>
192		// For pvtxs with splittracks1
193		pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
194		pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
#	Line 225 | Line 197 | PVStudy::PVStudy(const edm::ParameterSet
197		pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
198		pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
199		pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
200	+	pvtxtree_->Branch("errx_spl1_mct",&errx_spl1_mct_,"errx_spl1_mct[nrecPV_spl1_mct]/D");
201	+	pvtxtree_->Branch("erry_spl1_mct",&erry_spl1_mct_,"erry_spl1_mct[nrecPV_spl1_mct]/D");
202	+	pvtxtree_->Branch("errz_spl1_mct",&errz_spl1_mct_,"errz_spl1_mct[nrecPV_spl1_mct]/D");
203	+	pvtxtree_->Branch("nTrkPV_spl1_mct",&nTrkPV_spl1_mct_,"nTrkPV_spl1_mct[nrecPV_spl1_mct]/I");
204	+	pvtxtree_->Branch("ndofPV_spl1_mct",&ndofPV_spl1_mct_,"ndofPV_spl1_mct[nrecPV_spl1_mct]/D");
205	+	pvtxtree_->Branch("normchi2PV_spl1_mct",&normchi2PV_spl1_mct_,"normchi2PV_spl1_mct[nrecPV_spl1_mct]/D");
206	+	pvtxtree_->Branch("avgPtPV_spl1_mct",&avgPtPV_spl1_mct_,"avgPtPV_spl1_mct[nrecPV_spl1_mct]/D");
207	+
208		// For pvtxs with splittracks1
209		pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
210		pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
#	Line 232 | Line 212 | PVStudy::PVStudy(const edm::ParameterSet
212		pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
213		pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
214		pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
215	<	pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
215	>	pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
216	>	pvtxtree_->Branch("errx_spl2_mct",&errx_spl2_mct_,"errx_spl2_mct[nrecPV_spl2_mct]/D");
217	>	pvtxtree_->Branch("erry_spl2_mct",&erry_spl2_mct_,"erry_spl2_mct[nrecPV_spl2_mct]/D");
218	>	pvtxtree_->Branch("errz_spl2_mct",&errz_spl2_mct_,"errz_spl2_mct[nrecPV_spl2_mct]/D");
219	>	pvtxtree_->Branch("nTrkPV_spl2_mct",&nTrkPV_spl2_mct_,"nTrkPV_spl2_mct[nrecPV_spl2_mct]/I");
220	>	pvtxtree_->Branch("ndofPV_spl2_mct",&ndofPV_spl2_mct_,"ndofPV_spl2_mct[nrecPV_spl2_mct]/D");
221	>	pvtxtree_->Branch("normchi2PV_spl2_mct",&normchi2PV_spl2_mct_,"normchi2PV_spl2_mct[nrecPV_spl2_mct]/D");
222	>	pvtxtree_->Branch("avgPtPV_spl2_mct",&avgPtPV_spl2_mct_,"avgPtPV_spl2_mct[nrecPV_spl2_mct]/D");
223		}
224		}
225
#	Line 256 | Line 243 | PVStudy::PVStudy(const edm::ParameterSet
243
244		// Book Histograms:
245		h_pvtrk   = new PVHistograms();
259	–	h_pixvtx  = new PVHistograms();
246		h_misc    = new PVHistograms();
247		h_summary = new PVHistograms();
248		h_others  = new PVHistograms();
#	Line 269 | Line 255 | PVStudy::PVStudy(const edm::ParameterSet
255		h_pvtrk->Init("pvTrk", ss.str(),"spl");
256		}
257		}
272	–	h_pixvtx->Init("pixVtx");
258		h_misc->Init("misc");
259
260		// Book MC only plots
#	Line 318 | Line 303 | PVStudy::~PVStudy()
303		theFile->cd();
304		theFile->cd("Summary");
305		h_pvtrk->Save();
321	–	h_pixvtx->Save();
306		h_misc->Save();
307		h_summary->Save();
308		if (!realData_)
#	Line 356 | Line 340 | void PVStudy::setRootStyle() {
340		gStyle->SetTitleSize(0.055, "");    // size for pad title; default is 0.02
341		gStyle->SetLabelSize(0.03, "XYZ");  // size for axis labels; default is 0.04
342		gStyle->SetStatFontSize(0.08);      // size for stat. box
343	<	gStyle->SetTitleFont(32, "XYZ");    // times-bold-italic font (p. 153) for axes
344	<	gStyle->SetTitleFont(32, "");       // same for pad title
345	<	gStyle->SetLabelFont(32, "XYZ");    // same for axis labels
346	<	gStyle->SetStatFont(32);            // same for stat. box
343	>	gStyle->SetTitleFont(42, "XYZ");    // times-bold-italic font (p. 153) for axes
344	>	gStyle->SetTitleFont(42, "");       // same for pad title
345	>	gStyle->SetLabelFont(42, "XYZ");    // same for axis labels
346	>	gStyle->SetStatFont(42);            // same for stat. box
347		gStyle->SetLabelOffset(0.006, "Y"); // default is 0.005
348		//
349		return;
#	Line 546 | Line 530 | PVStudy::analyze(const edm::Event& iEven
530		using namespace edm;
531		using namespace reco;
532
533	+	//========================================================================
534	+	// Step 0: Prepare root variables and get information from the Event
535	+	//========================================================================
536	+
537		edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
550	–	//=======================================================
538		// Initialize Root-tuple variables if needed
539	<	//=======================================================
540	<	//if(saventuple_)
541	<	SetVarToZero();
555	<
556	<	//=======================================================
557	<	// Track accessors
558	<	//=======================================================
559	<	//trackCollection
539	>	SetVarToZero();
540	>
541	>	// ====== TrackCollection
542		static const reco::TrackCollection s_empty_trackColl;
543		const reco::TrackCollection *trackColl = &s_empty_trackColl;
544		edm::Handle<reco::TrackCollection> trackCollectionHandle;
#	Line 566 | Line 548 | PVStudy::analyze(const edm::Event& iEven
548		} else {
549		edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
550		}
551	<	//splitTrackCollection1
551	>	// ====== splitTrackCollection1
552		static const reco::TrackCollection s_empty_splitTrackColl1;
553		const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
554		edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
#	Line 576 | Line 558 | PVStudy::analyze(const edm::Event& iEven
558		} else {
559		edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
560		}
561	<	//splitTrackCollection2
561	>	// ====== splitTrackCollection2
562		static const reco::TrackCollection s_empty_splitTrackColl2;
563		const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
564		edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
#	Line 587 | Line 569 | PVStudy::analyze(const edm::Event& iEven
569		edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
570		}
571
572	<	//=======================================================
591	<	// PVTX accessors
592	<	//=======================================================
593	<	//vertexCollection
572	>	// ======= PrimaryVertexCollection
573		static const reco::VertexCollection s_empty_vertexColl;
574		const reco::VertexCollection *vertexColl = &s_empty_vertexColl;
575		edm::Handle<reco::VertexCollection> vertexCollectionHandle;
#	Line 600 | Line 579 | PVStudy::analyze(const edm::Event& iEven
579		} else {
580		edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
581		}
582	<	//splitVertexCollection1
582	>	// ====== splitVertexCollection1
583		static const reco::VertexCollection s_empty_splitVertexColl1;
584		const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
585		edm::Handle<reco::VertexCollection> splitVertexCollection1Handle;
#	Line 610 | Line 589 | PVStudy::analyze(const edm::Event& iEven
589		} else {
590		edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
591		}
592	<	//splitVertexCollection2
592	>	// ====== splitVertexCollection2
593		static const reco::VertexCollection s_empty_splitVertexColl2;
594		const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
595		edm::Handle<reco::VertexCollection> splitVertexCollection2Handle;
#	Line 622 | Line 601 | PVStudy::analyze(const edm::Event& iEven
601		}
602
603
604	<	//=======================================================
626	<	// GET pixelVertices
627	<	//=======================================================
628	<	static const reco::VertexCollection s_empty_pixelVertexColl;
629	<	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
630	<	edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
631	<	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
632	<	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
633	<	pixelVertexColl = pixelVertexCollectionHandle.product();
634	<	} else {
635	<	edm::LogInfo("Debug") << "[PVStudy] pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
636	<	}
637	<
638	<	//=======================================================
639	<	// BeamSpot accessors
640	<	//=======================================================
641	<
604	>	// ======== BeamSpot accessors
605		edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
606		iEvent.getByLabel(bsSrc,recoBeamSpotHandle);
607		reco::BeamSpot bs = *recoBeamSpotHandle;
608		const Point beamSpot = recoBeamSpotHandle.isValid() ? Point(recoBeamSpotHandle->x0(), recoBeamSpotHandle->y0(), recoBeamSpotHandle->z0()) : Point(0, 0, 0);
609
610	<	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex and pixelvertices collections"<<endl;
610	>	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex collections"<<endl;
611
612	<	//=======================================================
650	<	// MC simvtx accessor
651	<	//=======================================================
612	>	// ========== MC simvtx accessor
613		if (!realData_) {
614		edm::Handle<SimVertexContainer> simVtxs;
615		iEvent.getByLabel( simG4_, simVtxs);
#	Line 657 | Line 618 | PVStudy::analyze(const edm::Event& iEven
618		iEvent.getByLabel( simG4_, simTrks);
619		}
620
621	<	//=======================================================
661	<	// GET PDT
662	<	//=======================================================
621	>	// ========== GET PDT
622		try{
623		iSetup.getData(pdt);
624		}catch(const Exception&){
625		edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
626		}
627	<
628	<	//=======================================================
629	<	// Apply event cleaning for firstcoll data and MC
630	<	//=======================================================
631	<	if(applyEvtClean_) {
632	<	// =====Select on the trigger bits
633	<	edm::ESHandle<L1GtTriggerMenu> menuRcd;
634	<	iSetup.get<L1GtTriggerMenuRcd>().get(menuRcd) ;
635	<	const L1GtTriggerMenu* menu = menuRcd.product();
636	<	edm::Handle< L1GlobalTriggerReadoutRecord > gtRecord;
637	<	iEvent.getByLabel( edm::InputTag("gtDigis"), gtRecord);
638	<
680	<	bool isTechBit0 = false;
681	<	bool isTechBit40 = false;
682	<	bool isBeamHalo = false;
683	<
684	<	TechnicalTriggerWord tw = gtRecord->technicalTriggerWord();
685	<	if ( ! tw.empty() ) {
686	<	// loop over dec. bit to get total rate (no overlap)
687	<	for ( int i = 0; i < 64; ++i ) {
688	<	if ( tw[i] ) {
689	<	//cout<<"technical number "<<i<<"  "<<endl;
690	<	if (i == 0) isTechBit0  = true;
691	<	if (i < 40 && i > 35) isBeamHalo = true;   // The beamHalo bits are 36-39
692	<	if (i == 40 || i == 41) isTechBit40 = true;
693	<	}
694	<	}
695	<	} // =====End select on the trigger bits
696	<	// =====Loop over the trackColl to tet the fraction of HighPurity tracks
697	<	int nTracks = 0;
698	<	int nHighPurityTracks=0;
699	<	double fHighPurity=0;
700	<
701	<	for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
702	<	TrackRef tkref(trackCollectionHandle,i);
703	<	if( (tkref->qualityMask() & 4 ) == 4) ++nHighPurityTracks;
704	<	}
705	<	if(nTracks>0)
706	<	fHighPurity =  double(nHighPurityTracks)/double(nTracks);
707	<	if(verbose_)
708	<	cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;
709	<
710	<	// Apply the event selection for MC events
711	<	if(realData_) {
712	<	if(!isTechBit40 || isBeamHalo || (fHighPurity<0.2&&nTracks>10) || vertexColl->begin()->isFake()) {
713	<	glob_runno_ = iEvent.id().run();
714	<	glob_evtno_ = iEvent.id().event();
715	<	glob_ls_   = iEvent.luminosityBlock();
716	<	glob_bx_  = iEvent.bunchCrossing();
717	<	return;
718	<	}
719	<	else {
720	<	// TechBit 0, beamHalo are not simulted in MC
721	<	if ( !isTechBit0 || !isTechBit40 || isBeamHalo || (fHighPurity<0.2&&nTracks>10) || vertexColl->begin()->isFake()) return;
722	<	}
627	>
628	>	// ======= Get MC information if needed
629	>	bool MC=false;
630	>	Handle<HepMCProduct> evtMC;
631	>	if (!realData_) {
632	>	iEvent.getByLabel("generator",evtMC);
633	>	if (!evtMC.isValid()) {
634	>	MC=false;
635	>	edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
636	>	} else {
637	>	edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
638	>	MC=true;
639		}
640	<	} // End of Apply event cleaning cuts for firstcoll data
725	<
726	<
640	>	}
641
642	<	//=======================================================
643	<	// Fill trackparameters of the input tracks to pvtx fitter
644	<	//=======================================================
642	>	//========================================================================
643	>	// Step 1:  Apply event cleaning for data and MC
644	>	//          WARNING: event selection cut are hard coded!!
645	>	//========================================================================
646	>
647	>	glob_runno_ = iEvent.id().run();
648	>	glob_evtno_ = iEvent.id().event();
649	>	glob_ls_   = iEvent.luminosityBlock();
650	>	glob_bx_  = iEvent.bunchCrossing();
651	>
652	>	//========================================================================
653	>	// Step 2: Fill histograms for the splitting consistency checks
654	>	//========================================================================
655	>
656	>	// === Fill trackparameters of the input tracks to pvtx fitter
657		edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
658		//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs)
659		// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
#	Line 735 | Line 661 | PVStudy::analyze(const edm::Event& iEven
661		fillTrackHisto(splitTrackColl1, 1, beamSpot);
662		fillTrackHisto(splitTrackColl2, 2, beamSpot);
663		edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
664	<
665	<
666	<	//=======================================================
741	<	// Fill pixelVertices related histograms
742	<	//=======================================================
743	<	nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
744	<	if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
745	<	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
746	<	fillTrackHistoInPV(pixelVertexColl, 4, false, true, beamSpot);
747	<	h_pixvtx->Fill1d("dzErr_pxlpvtx", pixelVertexColl->begin()->zError());
748	<	// Get the dZ error of the tracks in the leading pixelVertexColl
749	<	for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
750	<	t!= (pixelVertexColl->begin())->tracks_end(); t++) {
751	<
752	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
753	<	//        h_pvtrk->Fill1d("trkPtPV", 0.);
754	<	}
755	<	else
756	<	h_pixvtx->Fill1d("trkdzErr_pxlpvtx", (**t).dzError());
757	<	}
758	<	// Fill the track->dz() in the PV difference from first pixelpvtx
759	<	for(reco::Vertex::trackRef_iterator t = (vertexColl->begin())->tracks_begin();
760	<	t!= (vertexColl->begin())->tracks_end(); t++) {
761	<	// illegal charge
762	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
763	<	//        h_pvtrk->Fill1d("trkPtPV", 0.);
764	<	}
765	<	else {
766	<	if(pixelVertexColl->size()>0) {
767	<	h_pixvtx->Fill1d("trkdz_pxlpvtxdz", (**t).dz() -  pixelVertexColl->begin()->z());
768	<	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_pxlpvtxdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
769	<	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_trkdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
770	<	h_pixvtx->Fill1d("trkdzErr_pvtx", (**t).dzError());
771	<	}
772	<	}
773	<	}
774	<	}
775	<
776	<	//=======================================================
777	<	// Fill number of reconstructed vertices
778	<	//=======================================================
779	<
664	>
665	>
666	>	// ==== Fill number of reconstructed vertices
667		edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
668		edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
669		edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
#	Line 786 | Line 673 | PVStudy::analyze(const edm::Event& iEven
673		printRecVtxs(splitVertexCollection1Handle);
674		printRecVtxs(splitVertexCollection2Handle);
675		}
676	+
677		nrecPV_ = int (vertexColl->size());
678		nrecPV1_spl_ = int (splitVertexColl1->size());
679		nrecPV2_spl_ = int (splitVertexColl2->size());
#	Line 796 | Line 684 | PVStudy::analyze(const edm::Event& iEven
684		h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
685
686
687	<	//=================================================================
688	<	// Fill track parameter ntuple/hist for tracks used in recoVertices
801	<	//=================================================================
802	<
803	<	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
687	>	// ======= Fill track parameter ntuple/hist for tracks used in recoVertices
688	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
689		if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
690	<	fillTrackHistoInPV(vertexColl, 0, true, true, beamSpot);
806	<
807	<	if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
808	<	fillTrackHistoInPV(splitVertexColl1, 1, false, true, beamSpot);
809	<
810	<	if(splitVertexColl2->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
811	<	fillTrackHistoInPV(splitVertexColl2, 2, false, true, beamSpot);
812	<
813	<	//=======================================================
814	<	// Compare offlinePrimaryVertices with pixelVertices
815	<	//=======================================================
816	<	if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
817	<	&& (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {
818	<	h_pixvtx->Fill1d("nrecPV_minus_nrecPxlPV", double (nrecPV_ - nrecPV_pxlpvtx_));
819	<	// difference in reconstructed position of the leading pvtx
820	<	//edm::LogInfo("Debug")<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;
821	<	//edm::LogInfo("Debug")<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
822	<	h_pixvtx->Fill1d("recxPV_minus_recxPxlPV", recx_[0] - recx_pxlpvtx_[0]);
823	<	h_pixvtx->Fill1d("recyPV_minus_recyPxlPV", recy_[0] - recy_pxlpvtx_[0]);
824	<	h_pixvtx->Fill1d("reczPV_minus_reczPxlPV", recz_[0] - recz_pxlpvtx_[0]);
825	<	}
826	<
827	<
828	<	//==========================================================
829	<	// Fill secondary/primary min separations for multi-vertices
830	<	//==========================================================
690	>	fillTrackHistoInPV(vertexColl, 0, beamSpot);
691
692	+	// ======= Fill secondary/primary min separations for multi-vertices
693		if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {
833	–
694		double min_xsep = 9999.0;
695		double min_ysep = 9999.0;
696		double min_zsep = 9999.0;
#	Line 846 | Line 706 | PVStudy::analyze(const edm::Event& iEven
706		for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
707		v!=vertexColl->end(); ++v) {
708		if(v->isValid() && ! v->isFake()) {
849	–	// xsep
709		if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
710		min_xsep = fabs(v->x()- vertexColl->begin()->x());
852	–	// ysep
711		if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
712		min_ysep = fabs(v->y()- vertexColl->begin()->y());
855	–	// zsep
713		if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
714		min_zsep = fabs(v->z()- vertexColl->begin()->z());
858	–	// xsep_sign
715		double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
716		if(xsep_sign < min_xsep_sign)
717		min_xsep_sign = xsep_sign;
862	–	// ysep_sign
718		double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
719		if(ysep_sign < min_ysep_sign)
720		min_ysep_sign = ysep_sign;
866	–	// zsep_sign
721		double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
722		if(zsep_sign < min_zsep_sign)
723		min_zsep_sign = zsep_sign;
870	–	// ntrksep
724		if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
725		min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
873	–	// sumpt2sep
726		if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
727		min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
728		}
#	Line 887 | Line 739 | PVStudy::analyze(const edm::Event& iEven
739		min_zsep_ = min_zsep;
740		min_ntrksep_ = min_ntrksep;
741		min_sumpt2sep_ = min_sumpt2sep;
890	–
891	–
742		} // End of  if(vertexColl->size()>1) {
743
744		edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
745
746	<
747	<
748	<	//=======================================================
749	<	//           PrimaryVertex Matching
750	<	// 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
751	<	// 2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
902	<	// Assume the first match is the primary vertex,
903	<	// since vertexColl are sorted in decreasing order of sum(pT)
904	<	//=======================================================
746	>	//========================================================================
747	>	// Step 3:  PrimaryVertex Matching
748	>	//   1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
749	>	//   2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
750	>	//   3. The first match is the primary vertex, which has largest sum(pT^2)
751	>	//========================================================================
752
753		edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
754		reco::VertexCollection s_empty_matchedVertexColl1;
#	Line 931 | Line 778 | PVStudy::analyze(const edm::Event& iEven
778		if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
779		edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
780
781	<	int nTrkPV1 = recvtx1->tracksSize();
782	<	int nTrkPV2 = recvtx2->tracksSize();
781	>	int nTrkPV1, nTrkPV2;
782	>	if(useHWTrk_) {
783	>	nTrkPV1 = nHWTrkRecVtx(*recvtx1);
784	>	nTrkPV2 = nHWTrkRecVtx(*recvtx2);
785	>	}
786	>	else {
787	>	nTrkPV1 = recvtx1->tracksSize();
788	>	nTrkPV2 = recvtx2->tracksSize();
789	>	}
790		double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
791		h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
792		h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
#	Line 955 | Line 809 | PVStudy::analyze(const edm::Event& iEven
809		}
810		}
811		edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
812	<
813	<	//=======================================================
814	<	//  Analyze matchedVertexColls
815	<	//=======================================================
812	>
813	>	//========================================================================
814	>	// Step 4: Analyze matchedVertexColls
815	>	//========================================================================
816		edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
817
818	+	// ==== If it is MC, analyze the res/pull of the unsplit vertexColl first
819	+	if(MC){
820	+	// make a list of primary vertices:
821	+	std::vector<simPrimaryVertex> simpv;
822	+	simpv=getSimPVs(evtMC,"");
823	+	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
824	+	h_gen->Fill1d("nsimPV", simpv.size());
825	+	nsimPV_ = simpv.size();
826	+	int isimpv = 0;
827	+	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
828	+	vsim!=simpv.end(); vsim++, isimpv++){
829	+	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
830	+	simx_[isimpv] = vsim->x;
831	+	simy_[isimpv] = vsim->y;
832	+	simz_[isimpv] = vsim->y;
833	+	simptsq_[isimpv] = vsim->ptsq;
834	+	fillMCHisto(vsim, isimpv, vertexColl, 3, avgTrkPtInPVMin_, avgTrkPtInPVMax_);
835	+	}
836	+	}
837	+
838		// Compare the reconstructed vertex position and calculate resolution/pulls
839		if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
840	<	//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
841	<	fillTrackHistoInPV(matchedVertexColl1, 1, true, false, beamSpot);
842	<	fillTrackHistoInPV(matchedVertexColl2, 2, true, false, beamSpot);
840	>	// ==== Analyze the splitted vtx using MC method
841	>	if(MC){
842	>	// make a list of primary vertices:
843	>	std::vector<simPrimaryVertex> simpv;
844	>	simpv=getSimPVs(evtMC,"");
845	>	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
846	>	int isimpv = 0;
847	>	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
848	>	vsim!=simpv.end(); vsim++, isimpv++){
849	>	fillMCHisto(vsim, isimpv, matchedVertexColl1, 1, avgTrkPtInPVMin_, avgTrkPtInPVMax_);
850	>	fillMCHisto(vsim, isimpv, matchedVertexColl2, 2, avgTrkPtInPVMin_, avgTrkPtInPVMax_);
851	>	}
852	>	}
853	>
854	>	// ==== Analyze res/pull two-vertex method
855	>	//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype)
856	>	fillTrackHistoInPV(matchedVertexColl1, 1, beamSpot);
857	>	fillTrackHistoInPV(matchedVertexColl2, 2, beamSpot);
858
859		reco::VertexCollection::const_iterator v1;
860		reco::VertexCollection::const_iterator v2;
861	+
862		nrecPV_twovtx_ = 0;
863		for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
864		v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
865	<	++v1, ++v2) {
865	>	++v1, ++v2) {
866	>
867	>	// ==================================================================
868	>	// With the option to fill the histograms at a given average pT range
869	>	// ==================================================================
870	>	if (avgPtRecVtx(*v1) < avgTrkPtInPVMin_ || avgPtRecVtx(*v1) > avgTrkPtInPVMax_ ) continue;
871	>	if (avgPtRecVtx(*v2) < avgTrkPtInPVMin_ || avgPtRecVtx(*v2) > avgTrkPtInPVMax_ ) continue;
872	>
873		h_misc->Fill1d("ndofPVDiff", v1->ndof() - v2->ndof());
874		h_misc->Fill1d("ndofPVRelDiff", (v1->ndof()-v2->ndof())/(v1->ndof()+v2->ndof()));
875		fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
#	Line 981 | Line 878 | PVStudy::analyze(const edm::Event& iEven
878		ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
879		ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
880		ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
881	<	fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
881	>
882	>	int nTrkPV1, nTrkPV2;
883	>	if(useHWTrk_) {
884	>	nTrkPV1 = nHWTrkRecVtx(*v1);
885	>	nTrkPV2 = nHWTrkRecVtx(*v2);
886	>	}
887	>	else {
888	>	nTrkPV1 = v1->tracksSize();
889	>	nTrkPV2 = v2->tracksSize();
890	>	}
891	>
892	>	fntrk_ = (nTrkPV1 + nTrkPV2)/2;
893
894		h_summary->Fill1d("deltax", fres_[0] );
895		h_summary->Fill1d("deltay", fres_[1] );
#	Line 996 | Line 904 | PVStudy::analyze(const edm::Event& iEven
904		error(ferror_[0], ferror_[1], ferror_[2]),
905		fntrk_));
906		// Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
907	+
908		fillHisto(res(fres_[0], fres_[1], fres_[2]),
909		error(ferror_[0], ferror_[1], ferror_[2]),
910		fntrk_,0);
911	+
912		// Fill the ntuple variables
1003	–	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
1004	–	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
1005	–	sumptsq1_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v1);
1006	–	sumptsq2_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v2);
913		nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
914		deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
915		deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
#	Line 1014 | Line 920 | PVStudy::analyze(const edm::Event& iEven
920		pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
921		pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
922		pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];
1017	–	nrecPV_twovtx_++;
1018	–	} // End of analyzing res/pull
923
1020	–	//=======================================================
1021	–	// Fill simulated vertices
1022	–	//=======================================================
1023	–	if (!realData_) {
1024	–	bool MC=false;
1025	–	Handle<HepMCProduct> evtMC;
1026	–	iEvent.getByLabel("generator",evtMC);
1027	–	if (!evtMC.isValid()) {
1028	–	MC=false;
1029	–	edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
1030	–	} else {
1031	–	edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
1032	–	MC=true;
1033	–	}
1034	–
1035	–	if(MC){
1036	–	// make a list of primary vertices:
1037	–	std::vector<simPrimaryVertex> simpv;
1038	–	simpv=getSimPVs(evtMC,"");
1039	–	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
1040	–	h_gen->Fill1d("nsimPV", simpv.size());
1041	–	nsimPV_ = simpv.size();
1042	–
1043	–	int isimpv = 0;
1044	–	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
1045	–	vsim!=simpv.end(); vsim++, isimpv++){
1046	–	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
1047	–	simx_[isimpv] = vsim->x;
1048	–	simy_[isimpv] = vsim->y;
1049	–	simz_[isimpv] = vsim->y;
1050	–	simptsq_[isimpv] = vsim->ptsq;
1051	–
1052	–	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
1053	–	fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
1054	–	fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
1055	–	fillMCHisto(vsim, isimpv, vertexColl, 3);
1056	–	}
1057	–	} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
1058	–	} // End of   if (!realData_) {
924
925	+	//SplittedVertex
926	+	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = nTrkPV1;
927	+	ndofPV1_spl_twovtx_[nrecPV_twovtx_] = v1->ndof();
928	+	normchi2PV1_spl_twovtx_[nrecPV_twovtx_] = v1->normalizedChi2();
929	+	avgPtPV1_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v1);
930	+	errx1_spl_twovtx_[nrecPV_twovtx_] = v1->xError();
931	+	erry1_spl_twovtx_[nrecPV_twovtx_] = v1->yError();
932	+	errz1_spl_twovtx_[nrecPV_twovtx_] = v1->zError();
933	+
934	+	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = nTrkPV2;
935	+	ndofPV2_spl_twovtx_[nrecPV_twovtx_] = v2->ndof();
936	+	normchi2PV2_spl_twovtx_[nrecPV_twovtx_] = v2->normalizedChi2();
937	+	avgPtPV2_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v2);
938	+	errx2_spl_twovtx_[nrecPV_twovtx_] = v2->xError();
939	+	erry2_spl_twovtx_[nrecPV_twovtx_] = v2->yError();
940	+	errz2_spl_twovtx_[nrecPV_twovtx_] = v2->zError();
941	+
942	+	nrecPV_twovtx_++;
943	+
944	+	// Print some information of the two tracks events
945	+	if(verbose_ && fntrk_ < 4) {
946	+	cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
947	+	printRecVtx(*v1);
948	+	cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
949	+	printRecVtx(*v2);
950	+	}
951	+	} // End of analyzing res/pull
952		} // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
953		else
954		edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
#	Line 1064 | Line 956 | PVStudy::analyze(const edm::Event& iEven
956		edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
957
958
959	<
960	<	// Finally fill the ftree_
959	>	//========================================================================
960	>	// Step 5: Fill ntuple if saventuple_ is on
961	>	//========================================================================
962		if(saventuple_) {
963		ftree_->Fill();
964		pvtxtree_->Fill();
#	Line 1132 | Line 1025 | bool PVStudy::matchVertex(const PVStudy:
1025		{
1026		if(vrec.isFake() || !vrec.isValid()) return false;
1027		else
1028	<	return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1028	>	return true;
1029	>	//return (fabs(vsim.z-vrec.z())<0.0500); // =500um // remove this hard cut
1030		}
1031
1032		// Match two reconstructed vertices
#	Line 1199 | Line 1093 | void PVStudy::printRecVtxs(const Handle<
1093		}
1094		}
1095
1096	+	void PVStudy::printRecVtx(const reco::Vertex & v){
1097	+
1098	+	cout << "#trk " << std::setw(3) << v.tracksSize()
1099	+	<< " chi2 " << std::setw(4) << v.chi2()
1100	+	<< " ndof " << std::setw(3) << v.ndof() << endl
1101	+	<< " x "  << std::setw(8) <<std::fixed << std::setprecision(4) << v.x()
1102	+	<< " dx " << std::setw(8) << v.xError()<< endl
1103	+	<< " y "  << std::setw(8) << v.y()
1104	+	<< " dy " << std::setw(8) << v.yError()<< endl
1105	+	<< " z "  << std::setw(8) << v.z()
1106	+	<< " dz " << std::setw(8) << v.zError()
1107	+	<< endl;
1108	+	}
1109	+
1110		void PVStudy::printSimTrks(const Handle<SimTrackContainer> simTrks){
1111		cout <<  " simTrks   type, (momentum), vertIndex, genpartIndex"  << endl;
1112		int i=1;
#	Line 1217 | Line 1125 | void PVStudy::printSimTrks(const Handle<
1125		}
1126		}
1127
1128	<	void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1128	>	void PVStudy::fillHisto(res r, error er, int ntk, int datamode) {
1129		stringstream ss;
1130		ss << ntk;
1131		string suffix;
#	Line 1370 | Line 1278 | void PVStudy::fillTrackHisto(const reco:
1278		}
1279		}
1280
1281	<	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
1281	>	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
1282		string suffix;
1283		suffix = ""; // for unsplit tracks
1284		if(datatype == 1) suffix = "1_spl"; // for splittracks 1
#	Line 1378 | Line 1286 | void PVStudy::fillTrackHistoInPV(const r
1286		int ivtx = 0;
1287		for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1288		v!=vertexColl->end(); ++v, ++ivtx) {
1289	<	if(v->isFake()) continue;
1290	<	if(fillNtuple) {
1291	<	if(datatype == 0) {
1292	<	nTrkPV_[ivtx] = v->tracksSize();
1293	<	sumptsq_[ivtx] = sumPtSquared(*v);
1294	<	isValid_[ivtx] = v->isValid();
1295	<	isFake_[ivtx] = v->isFake();
1296	<	recx_[ivtx] = v->x();
1297	<	recy_[ivtx] = v->y();
1298	<	recz_[ivtx] = v->z();
1299	<	recx_err_[ivtx] = v->xError();
1300	<	recy_err_[ivtx] = v->yError();
1301	<	recz_err_[ivtx] = v->zError();
1302	<	}
1303	<	if(datatype == 1) {
1304	<	nTrkPV1_spl_[ivtx] = v->tracksSize();
1305	<	sumptsq1_spl_[ivtx] = sumPtSquared(*v);
1306	<	isValid1_spl_[ivtx] = v->isValid();
1307	<	isFake1_spl_[ivtx] = v->isFake();
1400	<	recx1_spl_[ivtx] = v->x();
1401	<	recy1_spl_[ivtx] = v->y();
1402	<	recz1_spl_[ivtx] = v->z();
1403	<	recx1_err_spl_[ivtx] = v->xError();
1404	<	recy1_err_spl_[ivtx] = v->yError();
1405	<	recz1_err_spl_[ivtx] = v->zError();
1406	<
1407	<	}
1408	<	if(datatype == 2) {
1409	<	nTrkPV2_spl_[ivtx] = v->tracksSize();
1410	<	sumptsq2_spl_[ivtx] = sumPtSquared(*v);
1411	<	isValid2_spl_[ivtx] = v->isValid();
1412	<	isFake2_spl_[ivtx] = v->isFake();
1413	<	recx2_spl_[ivtx] = v->x();
1414	<	recy2_spl_[ivtx] = v->y();
1415	<	recz2_spl_[ivtx] = v->z();
1416	<	recx2_err_spl_[ivtx] = v->xError();
1417	<	recy2_err_spl_[ivtx] = v->yError();
1418	<	recz2_err_spl_[ivtx] = v->zError();
1419	<	}
1420	<	if(datatype == 4) { // for pixelVertices
1421	<	nTrkPV_pxlpvtx_[ivtx] = v->tracksSize();
1422	<	sumptsq_pxlpvtx_[ivtx] = sumPtSquared(*v);
1423	<	isValid_pxlpvtx_[ivtx] = v->isValid();
1424	<	isFake_pxlpvtx_[ivtx] = v->isFake();
1425	<	recx_pxlpvtx_[ivtx] = v->x();
1426	<	recy_pxlpvtx_[ivtx] = v->y();
1427	<	recz_pxlpvtx_[ivtx] = v->z();
1428	<	recx_err_pxlpvtx_[ivtx] = v->xError();
1429	<	recy_err_pxlpvtx_[ivtx] = v->yError();
1430	<	recz_err_pxlpvtx_[ivtx] = v->zError();
1431	<	}
1432	<	}
1433	<	}
1434	<	// For histogram only fill the leading pvtx parameters
1435	<	if(fillHisto) {
1436	<	h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1437	<	h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1438	<	int nHWTrkPV_ = 0;
1439	<
1440	<	try {
1441	<	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1442	<	t!=vertexColl->begin()->tracks_end(); t++) {
1443	<	// illegal charge
1444	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1445	<	//      h_pvtrk->Fill1d("trkPtPV", 0.);
1446	<	}
1447	<	else {
1448	<	h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1449	<	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1450	<	h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
1451	<	h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1452	<	h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1453	<	h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1454	<
1455	<	if(vertexColl->begin()->trackWeight(*t) > 0.5)
1456	<	nHWTrkPV_++;
1289	>	if(!v->isFake()) {
1290	>	h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1291	>	h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1292	>	h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkRecVtx(*v));
1293	>	try {
1294	>	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1295	>	t!=vertexColl->begin()->tracks_end(); t++) {
1296	>	// illegal charge
1297	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1298	>	//    h_pvtrk->Fill1d("trkPtPV", 0.);
1299	>	}
1300	>	else {
1301	>	h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1302	>	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1303	>	h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
1304	>	h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1305	>	h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1306	>	h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1307	>	}
1308		}
1309		}
1310	+	catch (...) {
1311	+	// exception thrown when trying to use linked track
1312	+	//        h_pvtrk->Fill1d("trkPtPV", 0.);
1313	+	}
1314		}
1460	–	catch (...) {
1461	–	// exception thrown when trying to use linked track
1462	–	//          h_pvtrk->Fill1d("trkPtPV", 0.);
1463	–	}
1464	–	h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
1315		}
1466	–
1467	–
1316		}
1317
1318	<	void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
1318	>	void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl,
1319	>	int datatype, double avgTrkPtInPVMin_, double avgTrkPtInPVMax_)
1320		{
1321		std::string suffix;
1322		// Set the vertexColl and histogram suffix according to datamode
#	Line 1483 | Line 1332 | void PVStudy::fillMCHisto(std::vector<si
1332		suffix = "_mct";
1333		nrecPV_mct_ = 0;
1334		}
1335	<
1335	>
1336		//========================================================
1337	<	//  look for a matching reconstructed vertex in vertexColl
1338	<	//========================================================
1339	<
1337	>	//  For each simulated vertex, look for a match in the vertexColl
1338	>	//  If more than 1 recVtx is found, use the one with closest in Z
1339	>	//========================================================
1340	>
1341	>	// === look for a matching reconstructed vertex in vertexColl
1342		for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1343		vrec!=vtxColl->end(); ++vrec){
1493	–	int nrectrk = vrec->tracksSize();
1344		vsim->recVtx=NULL;
1495	–
1345		edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1346		edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1347	<
1347	>
1348		if ( matchVertex(*vsim,*vrec)  ) {
1349		vsim->recVtx=&(*vrec);
1350	<
1502	<	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1350	>	//if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1351		//if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
1352		//|| (!vsim->recVtx) )
1353		//vsim->recVtx=&(*vrec);
1506	–	//}
1507	–	edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1508	–
1509	–	double fres_mct[3];
1510	–	double ferror_mct[3];
1511	–
1512	–	fres_mct[0] = vsim->recVtx->x()-vsim->x;
1513	–	fres_mct[1] = vsim->recVtx->y()-vsim->y;
1514	–	fres_mct[2] = vsim->recVtx->z()-vsim->z;
1515	–	ferror_mct[0] = vsim->recVtx->xError();
1516	–	ferror_mct[1] = vsim->recVtx->yError();
1517	–	ferror_mct[2] = vsim->recVtx->zError();
1518	–
1519	–	h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1520	–	h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1521	–	h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1522	–	h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1523	–	h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1524	–	h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1525	–	h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1526	–	h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1527	–	h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1528	–	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1529	–	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1530	–	nrectrk));
1531	–	// Fill histo according to its track multiplicity
1532	–	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1533	–	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1534	–	nrectrk, datatype);
1535	–
1536	–	if(saventuple_) {
1537	–	//Fill the values for variables in ftree_
1538	–	if(datatype == 1) {
1539	–	nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1540	–	deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1541	–	deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1542	–	deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1543	–	pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
1544	–	pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1545	–	pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1546	–	nrecPV_spl1_mct_++;
1547	–	}
1548	–
1549	–	if(datatype == 2) {
1550	–	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1551	–	deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1552	–	deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1553	–	deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1554	–	pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
1555	–	pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1556	–	pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1557	–	nrecPV_spl2_mct_++;
1558	–	}
1559	–	if(datatype == 3) {
1560	–	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1561	–	deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1562	–	deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1563	–	deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1564	–	pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
1565	–	pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1566	–	pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1567	–	nrecPV_mct_++;
1568	–	}
1569	–	} // End of  if(saventuple_) {
1570	–	} //  if ( matchVertex(*vsim,*vrec) ) {
1571	–	else {  // no rec vertex found for this simvertex
1572	–	edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1354		}
1355		}
1356	+
1357	+	// === If match found fill the residual and pulls
1358	+	if(vsim->recVtx) {
1359	+	edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1360	+	double avgPtPV = avgPtRecVtx(*(vsim->recVtx));
1361	+	if(avgPtPV>avgTrkPtInPVMax_ || avgPtPV<avgTrkPtInPVMin_) return;
1362	+
1363	+	// if study the resolution/pull for all tracks in the PVtx
1364	+	int nrectrk;
1365	+	if(useHWTrk_)
1366	+	nrectrk=nHWTrkRecVtx(*(vsim->recVtx));
1367	+	else
1368	+	nrectrk =  int(vsim->recVtx->tracksSize());
1369	+
1370	+	double fres_mct[3];
1371	+	double ferror_mct[3];
1372	+
1373	+	double ndofPV = vsim->recVtx->ndof();
1374	+	double normchi2PV =  vsim->recVtx->normalizedChi2();
1375	+
1376	+	fres_mct[0] = vsim->recVtx->x()-vsim->x;
1377	+	fres_mct[1] = vsim->recVtx->y()-vsim->y;
1378	+	fres_mct[2] = vsim->recVtx->z()-vsim->z;
1379	+	ferror_mct[0] = vsim->recVtx->xError();
1380	+	ferror_mct[1] = vsim->recVtx->yError();
1381	+	ferror_mct[2] = vsim->recVtx->zError();
1382	+
1383	+	h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1384	+	h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1385	+	h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1386	+	h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1387	+	h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1388	+	h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1389	+	h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1390	+	h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1391	+	h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1392	+	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1393	+	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1394	+	nrectrk));
1395	+	// Fill histo according to its track multiplicity
1396	+	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1397	+	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1398	+	nrectrk, datatype);
1399	+
1400	+	if(saventuple_) {
1401	+	//Fill the values for variables in pvtxtree_
1402	+	if(datatype == 1) {
1403	+	nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1404	+	ndofPV_spl1_mct_[nrecPV_spl1_mct_] =  ndofPV;
1405	+	normchi2PV_spl1_mct_[nrecPV_spl1_mct_] =  normchi2PV;
1406	+	avgPtPV_spl1_mct_[nrecPV_spl1_mct_] = avgPtPV;
1407	+	deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1408	+	deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[1];
1409	+	deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2];
1410	+	pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
1411	+	pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1412	+	pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1413	+	errx_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[0];
1414	+	erry_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[1];
1415	+	errz_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[2];
1416	+	nrecPV_spl1_mct_++;
1417	+	}
1418	+	if(datatype == 2) {
1419	+	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1420	+	ndofPV_spl2_mct_[nrecPV_spl2_mct_] =  ndofPV;
1421	+	normchi2PV_spl2_mct_[nrecPV_spl2_mct_] =  normchi2PV;
1422	+	avgPtPV_spl2_mct_[nrecPV_spl2_mct_] = avgPtPV;
1423	+	deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1424	+	deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[1];
1425	+	deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2];
1426	+	pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
1427	+	pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1428	+	pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1429	+	errx_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[0];
1430	+	erry_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[1];
1431	+	errz_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[2];
1432	+	nrecPV_spl2_mct_++;
1433	+	}
1434	+	if(datatype == 3) {
1435	+	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1436	+	ndofPV_mct_[nrecPV_mct_] =  ndofPV;
1437	+	normchi2PV_mct_[nrecPV_mct_] =  normchi2PV;
1438	+	avgPtPV_mct_[nrecPV_mct_] = avgPtPV;
1439	+	deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1440	+	deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1441	+	deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1442	+	pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
1443	+	pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1444	+	pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1445	+	errx_mct_[nrecPV_mct_] =  ferror_mct[0];
1446	+	erry_mct_[nrecPV_mct_] =  ferror_mct[1];
1447	+	errz_mct_[nrecPV_mct_] =  ferror_mct[2];
1448	+	nrecPV_mct_++;
1449	+	}
1450	+	} // End of  if(saventuple_) {
1451	+	} //  if ( matchVertex(*vsim,*vrec) ) {
1452	+	else {  // no rec vertex found for this simvertex
1453	+	edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1454	+	}
1455		}
1456
1457	+
1458		void PVStudy::SetVarToZero() {
1578	–	//Greg's variables
1459		fntrk_ = 0;
1460		//pvtx position (x,y,z) residual and error
1461		for(int i = 0; i<3;i++)    {
#	Line 1583 | Line 1463 | void PVStudy::SetVarToZero() {
1463		ferror_[i] = 0;
1464		}
1465
1466	<	//Yanyan's variables
1467	<	// Number of reconstructed vertices
1466	>	// Event level information
1467	>	glob_runno_ = 0;
1468	>	glob_evtno_ = 0;
1469	>	glob_ls_ = 0;
1470	>	glob_bx_ = 0;
1471		nrecPV_ = 0;
1472		nrecPV1_spl_ = 0;
1473		nrecPV2_spl_ = 0;
#	Line 1600 | Line 1483 | void PVStudy::SetVarToZero() {
1483		min_ntrksep_ = 9999.0;
1484		min_sumpt2sep_ = -9999.0;
1485
1486	<
1486	>	// Variables filled per Vertex
1487		for (int i = 0; i < nMaxPVs_; i++) {
1605	–	// recoVertices with all tracks
1606	–	nTrkPV_[i] = 0; // Number of tracks in the pvtx
1607	–	sumptsq_[i] = 0;
1608	–	isValid_[i] = -1;
1609	–	isFake_[i] = -1;
1610	–	recx_[i] = 0;
1611	–	recy_[i] = 0;
1612	–	recz_[i] = 0;
1613	–	recx_err_[i] = 0;
1614	–	recy_err_[i] = 0;
1615	–	recz_err_[i] = 0;
1616	–
1617	–	// recoVertices with splitTrack1
1618	–	nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx
1619	–	sumptsq1_spl_[i] = 0;
1620	–	isValid1_spl_[i] = -1;
1621	–	isFake1_spl_[i] = -1;
1622	–	recx1_spl_[i] = 0;
1623	–	recy1_spl_[i] = 0;
1624	–	recz1_spl_[i] = 0;
1625	–	recx1_err_spl_[i] = 0;
1626	–	recy1_err_spl_[i] = 0;
1627	–	recz1_err_spl_[i] = 0;
1628	–
1629	–	// recoVertices with splitTrack2
1630	–	nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx
1631	–	sumptsq2_spl_[i] = 0;
1632	–	isValid2_spl_[i] = -1;
1633	–	isFake2_spl_[i] = -1;
1634	–	recx2_spl_[i] = 0;
1635	–	recy2_spl_[i] = 0;
1636	–	recz2_spl_[i] = 0;
1637	–	recx2_err_spl_[i] = 0;
1638	–	recy2_err_spl_[i] = 0;
1639	–	recz2_err_spl_[i] = 0;
1640	–
1641	–	//pixelVertices
1642	–	nTrkPV_pxlpvtx_[i] = 0; // Number of tracks in the pvtx
1643	–	sumptsq_pxlpvtx_[i] = 0;
1644	–	isValid_pxlpvtx_[i] = -1;
1645	–	isFake_pxlpvtx_[i] = -1;
1646	–	recx_pxlpvtx_[i] = 0;
1647	–	recy_pxlpvtx_[i] = 0;
1648	–	recz_pxlpvtx_[i] = 0;
1649	–	recx_err_pxlpvtx_[i] = 0;
1650	–	recy_err_pxlpvtx_[i] = 0;
1651	–	recz_err_pxlpvtx_[i] = 0;
1652	–
1488		// matched two-vertices
1489		nTrkPV1_spl_twovtx_[i] = 0;
1490	+	ndofPV1_spl_twovtx_[i] = 0;
1491	+	normchi2PV1_spl_twovtx_[i] = 0;
1492	+	avgPtPV1_spl_twovtx_[i] = 0;
1493	+	errx1_spl_twovtx_[i] = 0;
1494	+	erry1_spl_twovtx_[i] = 0;
1495	+	errz1_spl_twovtx_[i] = 0;
1496	+
1497		nTrkPV2_spl_twovtx_[i] = 0;
1498	<	nTrkPV_twovtx_[i] = 0;
1499	<	sumptsq1_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1500	<	sumptsq2_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1498	>	ndofPV2_spl_twovtx_[i] = 0;
1499	>	normchi2PV2_spl_twovtx_[i] = 0;
1500	>	avgPtPV2_spl_twovtx_[i] = 0;
1501	>	errx2_spl_twovtx_[i] = 0;
1502	>	erry2_spl_twovtx_[i] = 0;
1503	>	errz2_spl_twovtx_[i] = 0;
1504	>
1505		deltax_twovtx_[i] = 0;
1506		deltay_twovtx_[i] = 0;
1507		deltaz_twovtx_[i] = 0;
#	Line 1680 | Line 1526 | void PVStudy::SetVarToZero() {
1526		pullx_mct_[i] = 0;
1527		pully_mct_[i] = 0;
1528		pullz_mct_[i] = 0;
1529	<
1529	>	nTrkPV_mct_[i] = 0;
1530	>	ndofPV_mct_[i] = 0;
1531	>	normchi2PV_mct_[i] = 0;
1532	>	avgPtPV_mct_[i] = 0;
1533	>	errx_mct_[i] = 0;
1534	>	erry_mct_[i] = 0;
1535	>	errz_mct_[i] = 0;
1536	>
1537		deltax_spl1_mct_[i] = 0;
1538		deltay_spl1_mct_[i] = 0;
1539		deltaz_spl1_mct_[i] = 0;
1540		pullx_spl1_mct_[i] = 0;
1541		pully_spl1_mct_[i] = 0;
1542		pullz_spl1_mct_[i] = 0;
1543	<
1543	>	nTrkPV_spl1_mct_[i] = 0;
1544	>	ndofPV_spl1_mct_[i] = 0;
1545	>	normchi2PV_spl1_mct_[i] = 0;
1546	>	avgPtPV_spl1_mct_[i] = 0;
1547	>	errx_spl1_mct_[i] = 0;
1548	>	erry_spl1_mct_[i] = 0;
1549	>	errz_spl1_mct_[i] = 0;
1550	>
1551		deltax_spl2_mct_[i] = 0;
1552		deltay_spl2_mct_[i] = 0;
1553		deltaz_spl2_mct_[i] = 0;
1554		pullx_spl2_mct_[i] = 0;
1555		pully_spl2_mct_[i] = 0;
1556		pullz_spl2_mct_[i] = 0;
1557	+	nTrkPV_spl2_mct_[i] = 0;
1558	+	ndofPV_spl2_mct_[i] = 0;
1559	+	normchi2PV_spl2_mct_[i] = 0;
1560	+	avgPtPV_spl2_mct_[i] = 0;
1561	+	errx_spl2_mct_[i] = 0;
1562	+	erry_spl2_mct_[i] = 0;
1563	+	errz_spl2_mct_[i] = 0;
1564	+
1565		}
1566		}
1567
1568		double PVStudy::sumPtSquared(const reco::Vertex & v)  {
1569	<	double sum = 0.;
1570	<	double pT;
1571	<	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1572	<	pT = (**it).pt();
1573	<	sum += pT*pT;
1569	>	if(v.isFake() || (!v.isValid()) || v.tracksSize() == 0) return 0.0;
1570	>	else {
1571	>	double sum = 0.;
1572	>	double pT;
1573	>	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1574	>	pT = (**it).pt();
1575	>	sum += pT*pT;
1576	>	}
1577	>	return sum;
1578		}
1707	–	return sum;
1579		}
1580
1581
1582
1583	+	double PVStudy::avgPtRecVtx(const reco::Vertex & v)  {
1584	+
1585	+	if(v.isFake() || !v.isValid() || v.tracksSize()==0 ) return 0;
1586	+	else {
1587	+	int nHWTrk = 0;
1588	+	double sumpT = 0.;
1589	+	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1590	+	if(v.trackWeight(*it) > 0.5 ) {
1591	+	sumpT +=  (**it).pt();
1592	+	nHWTrk++;
1593	+	}
1594	+	}
1595	+	if(nHWTrk > 0)
1596	+	return sumpT/double(nHWTrk);
1597	+	else
1598	+	return 0;
1599	+	}
1600	+	}
1601	+
1602	+	int PVStudy::nHWTrkRecVtx(const reco::Vertex & v)  {
1603	+	int nHWTrkPV = 0;
1604	+	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1605	+	if(v.trackWeight(*it) > 0.5)
1606	+	nHWTrkPV++;
1607	+	}
1608	+	return nHWTrkPV;
1609	+	}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines