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.7 by yygao, Thu Oct 8 14:57:46 2009 UTC vs.
Revision 1.11 by yygao, Wed Nov 11 19:28:50 2009 UTC

#	Line 25 | Line 25 | Implementation:
25		#include <vector>
26		#include <iostream>
27		#include <sstream>
28	+	#include <algorithm>
29
30		// user include files
31		#include "FWCore/Framework/interface/Frameworkfwd.h"
#	Line 71 | Line 72 | PVStudy::PVStudy(const edm::ParameterSet
72		splitTrackCollection2Tag_     = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
73		vertexCollectionTag_       = iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection");
74		splitVertexCollection1Tag_       = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
75	<	splitVertexCollection2Tag_       = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
75	>	splitVertexCollection2Tag_       = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
76	>	pixelVertexCollectionTag_      = iConfig.getParameter<edm::InputTag>("pixelVertexCollectionTag");
77		verbose_         = iConfig.getUntrackedParameter<bool>("verbose",false);
78		realData_        = iConfig.getUntrackedParameter<bool>("realData",false);
79		analyze_         = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
80		saventuple_         = iConfig.getUntrackedParameter<bool>("saventuple",false);
81		outputfilename_  = iConfig.getUntrackedParameter<string>("OutputFileName");
82		ntrkdiffcut_         = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
83	+	zsigncut_         = iConfig.getUntrackedParameter<int>("zsigncut");
84		nTrkMin_ =  iConfig.getUntrackedParameter<int>("nTrkMin");
85		nTrkMax_ =  iConfig.getUntrackedParameter<int>("nTrkMax");
86
87		//now do what ever initialization is needed
88		pvinfo.clear();
89
90	<
90	>
91		// Specify the data mode vector
92		if(realData_) datamode.push_back(0);
93		else {
#	Line 102 | Line 105 | PVStudy::PVStudy(const edm::ParameterSet
105		ftree_->Branch("residual",&fres_,"fres_[3]/D");
106		ftree_->Branch("error",&ferror_,"ferror_[3]/D");
107		ftree_->Branch("nTrkPV",&fntrk_,"fntrk_/I");
108	<	ftree_->Branch("nrecPV",&nrecPV_,"nrecPV_/I");
109	<	ftree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl_/I");
110	<	ftree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl_/I");
111	<
112	<	// If it is MC, get the resA, resB
113	<	if(!realData_) {
114	<	ftree_->Branch("residual_spl1_mct",&fres_spl1_mct_,"fres_spl1_mct_[3]/D");
115	<	ftree_->Branch("error_spl1_mct",&ferror_spl1_mct_,"ferror_spl1_mct_[3]/D");
116	<	ftree_->Branch("nTrkPV_spl1_mct",&fntrk_spl1_mct_,"fntrk_spl1_mct_/I");
117	<	ftree_->Branch("residual_spl2_mct",&fres_spl2_mct_,"fres_spl2_mct_[3]/D");
118	<	ftree_->Branch("error_spl2_mct",&ferror_spl2_mct_,"ferror_spl2_mct_[3]/D");
119	<	ftree_->Branch("nTrkPV_spl2_mct",&fntrk_spl2_mct_,"fntrk_spl2_mct_/I");
120	<	ftree_->Branch("residual_mct",&fres_mct_,"fres_mct_[3]/D");
121	<	ftree_->Branch("error_mct",&ferror_mct_,"ferror_mct_[3]/D");
122	<	ftree_->Branch("nTrkPV_mct",&fntrk_mct_,"fntrk_mct_/I");
108	>
109	>	// pvtxtree_ analyzing the pvtxs ootb
110	>	pvtxtree_ = new TTree("pvtxtree","pvtxtree");
111	>	//pvtx using all tracks
112	>	pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
113	>	pvtxtree_->Branch("nTrkPV",&nTrkPV_,"nTrkPV[nrecPV]/I");
114	>	pvtxtree_->Branch("sumptsq",&sumptsq_,"sumptsq[nrecPV]/D");
115	>	pvtxtree_->Branch("isValid",&isValid_,"isValid/I");
116	>	pvtxtree_->Branch("isFake",&isFake_,"isFake/I");
117	>	pvtxtree_->Branch("recx",&recx_,"recx[nrecPV]/D");
118	>	pvtxtree_->Branch("recy",&recy_,"recy[nrecPV]/D");
119	>	pvtxtree_->Branch("recz",&recz_,"recz[nrecPV]/D");
120	>	pvtxtree_->Branch("recx_err",&recx_err_,"recx_err[nrecPV]/D");
121	>	pvtxtree_->Branch("recy_err",&recy_err_,"recy_err[nrecPV]/D");
122	>	pvtxtree_->Branch("recz_err",&recz_err_,"recz_err[nrecPV]/D");
123	>
124	>	pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
125	>	pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
126	>	pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
127	>
128	>	//pvtx using splittracks1
129	>	pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
130	>	pvtxtree_->Branch("nTrkPV1_spl",&nTrkPV1_spl_,"nTrkPV1_spl[nrecPV1_spl]/I");
131	>	pvtxtree_->Branch("sumptsq1_spl",&sumptsq1_spl_,"sumptsq1_spl[nrecPV1_spl]/D");
132	>	pvtxtree_->Branch("isValid1_spl",&isValid1_spl_,"isValid1_spl/I");
133	>	pvtxtree_->Branch("isFake1_spl",&isFake1_spl_,"isFake1_spl/I");
134	>	pvtxtree_->Branch("recx1_spl",&recx1_spl_,"recx1_spl[nrecPV1_spl]/D");
135	>	pvtxtree_->Branch("recy1_spl",&recy1_spl_,"recy1_spl[nrecPV1_spl]/D");
136	>	pvtxtree_->Branch("recz1_spl",&recz1_spl_,"recz1_spl[nrecPV1_spl]/D");
137	>	pvtxtree_->Branch("recx1_err_spl",&recx1_err_spl_,"recx1_err_spl[nrecPV1_spl]/D");
138	>	pvtxtree_->Branch("recy1_err_spl",&recy1_err_spl_,"recy1_err_spl[nrecPV1_spl]/D");
139	>	pvtxtree_->Branch("recz1_err_spl",&recz1_err_spl_,"recz1_err_spl[nrecPV1_spl]/D");
140	>
141	>	//pvtx using splittracks2
142	>	pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
143	>	pvtxtree_->Branch("nTrkPV2_spl",&nTrkPV2_spl_,"nTrkPV2_spl[nrecPV2_spl]/I");
144	>	pvtxtree_->Branch("sumptsq2_spl",&sumptsq2_spl_,"sumptsq2_spl[nrecPV2_spl]/D");
145	>	pvtxtree_->Branch("isValid2_spl",&isValid2_spl_,"isValid2_spl/I");
146	>	pvtxtree_->Branch("isFake2_spl",&isFake2_spl_,"isFake2_spl/I");
147	>	pvtxtree_->Branch("recx2_spl",&recx2_spl_,"recx2_spl[nrecPV2_spl]/D");
148	>	pvtxtree_->Branch("recy2_spl",&recy2_spl_,"recy2_spl[nrecPV2_spl]/D");
149	>	pvtxtree_->Branch("recz2_spl",&recz2_spl_,"recz2_spl[nrecPV2_spl]/D");
150	>	pvtxtree_->Branch("recx2_err_spl",&recx2_err_spl_,"recx2_err_spl[nrecPV2_spl]/D");
151	>	pvtxtree_->Branch("recy2_err_spl",&recy2_err_spl_,"recy2_err_spl[nrecPV2_spl]/D");
152	>	pvtxtree_->Branch("recz2_err_spl",&recz2_err_spl_,"recz2_err_spl[nrecPV2_spl]/D");
153	>
154	>	//pixeVertices
155	>	pvtxtree_->Branch("nrecPV_pxlpvtx",&nrecPV_pxlpvtx_,"nrecPV_pxlpvtx/I");
156	>	pvtxtree_->Branch("nTrkPV_pxlpvtx",&nTrkPV_pxlpvtx_,"nTrkPV_pxlpvtx[nrecPV_pxlpvtx]/I");
157	>	pvtxtree_->Branch("sumptsq_pxlpvtx",&sumptsq_pxlpvtx_,"sumptsq_pxlpvtx[nrecPV_pxlpvtx]/D");
158	>	pvtxtree_->Branch("isValid_pxlpvtx",&isValid_pxlpvtx_,"isValid_pxlpvtx/I");
159	>	pvtxtree_->Branch("isFake_pxlpvtx",&isFake_pxlpvtx_,"isFake_pxlpvtx/I");
160	>	pvtxtree_->Branch("recx_pxlpvtx",&recx_pxlpvtx_,"recx_pxlpvtx[nrecPV_pxlpvtx]/D");
161	>	pvtxtree_->Branch("recy_pxlpvtx",&recy_pxlpvtx_,"recy_pxlpvtx[nrecPV_pxlpvtx]/D");
162	>	pvtxtree_->Branch("recz_pxlpvtx",&recz_pxlpvtx_,"recz_pxlpvtx[nrecPV_pxlpvtx]/D");
163	>	pvtxtree_->Branch("recx_err_pxlpvtx",&recx_err_pxlpvtx_,"recx_err_pxlpvtx[nrecPV_pxlpvtx]/D");
164	>	pvtxtree_->Branch("recy_err_pxlpvtx",&recy_err_pxlpvtx_,"recy_err_pxlpvtx[nrecPV_pxlpvtx]/D");
165	>	pvtxtree_->Branch("recz_err_pxlpvtx",&recz_err_pxlpvtx_,"recz_err_pxlpvtx[nrecPV_pxlpvtx]/D");
166	>
167	>	//Fill the variables in the twovtx pair (recvtx1, recvtx2)
168	>	pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
169	>	pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
170	>	pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
171	>	pvtxtree_->Branch("sumptsq1_spl_twovtx",&sumptsq1_spl_twovtx_,"sumptsq1_spl_twovtx[nrecPV_twovtx]/D");
172	>	pvtxtree_->Branch("sumptsq2_spl_twovtx",&sumptsq2_spl_twovtx_,"sumptsq2_spl_twovtx[nrecPV_twovtx]/D");
173	>	pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
174	>	pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
175	>	pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
176	>	pvtxtree_->Branch("deltaz_twovtx",&deltaz_twovtx_,"deltaz_twovtx[nrecPV_twovtx]/D");
177	>	pvtxtree_->Branch("errx_twovtx",&errx_twovtx_,"errx_twovtx[nrecPV_twovtx]/D");
178	>	pvtxtree_->Branch("erry_twovtx",&erry_twovtx_,"erry_twovtx[nrecPV_twovtx]/D");
179	>	pvtxtree_->Branch("errz_twovtx",&errz_twovtx_,"errz_twovtx[nrecPV_twovtx]/D");
180	>	pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
181	>	pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
182	>	pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
183	>
184	>	// MC variables
185	>	if(!realData_) {
186	>	pvtxtree_->Branch("nsimPV",&nsimPV_,"nsimPV/I");
187	>	pvtxtree_->Branch("nsimTrkPV",&nsimTrkPV_,"nsimTrkPV[nsimPV]/I");
188	>	pvtxtree_->Branch("simx",&simx_,"simx[nsimPV]/D");
189	>	pvtxtree_->Branch("simy",&simy_,"simy[nsimPV]/D");
190	>	pvtxtree_->Branch("simz",&simz_,"simz[nsimPV]/D");
191	>	pvtxtree_->Branch("simptsq",&simptsq_,"simptsq[nsimPV]/D");
192	>
193	>	// For pvtxs with all the tracks
194	>	pvtxtree_->Branch("nrecPV_mct",&nrecPV_mct_,"nrecPV_mct/I");
195	>	pvtxtree_->Branch("deltax_mct",&deltax_mct_,"deltax_mct[nrecPV_mct]/D");
196	>	pvtxtree_->Branch("deltay_mct",&deltay_mct_,"deltay_mct[nrecPV_mct]/D");
197	>	pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
198	>	pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
199	>	pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
200	>	pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
201	>	// For pvtxs with splittracks1
202	>	pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
203	>	pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
204	>	pvtxtree_->Branch("deltay_spl1_mct",&deltay_spl1_mct_,"deltay_spl1_mct[nrecPV_spl1_mct]/D");
205	>	pvtxtree_->Branch("deltaz_spl1_mct",&deltaz_spl1_mct_,"deltaz_spl1_mct[nrecPV_spl1_mct]/D");
206	>	pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
207	>	pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
208	>	pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
209	>	// For pvtxs with splittracks1
210	>	pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
211	>	pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
212	>	pvtxtree_->Branch("deltay_spl2_mct",&deltay_spl2_mct_,"deltay_spl2_mct[nrecPV_spl2_mct]/D");
213	>	pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
214	>	pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
215	>	pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
216	>	pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
217		}
218		}
219
#	Line 131 | Line 228 | PVStudy::PVStudy(const edm::ParameterSet
228		//i=0 : x (as in unsplit track collection)
229		//i=1 : x1_spl_
230		//i=2 : x2_spl_
231	<
231	>
232	>
233		for(int i=0;i<3;i++)
234		{
235		string suffix;
#	Line 143 | Line 241 | PVStudy::PVStudy(const edm::ParameterSet
241		}
242		h["nTrk"+suffix]   = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
243		h["trkPt"+suffix]  = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
244	<	h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#Eta of rec tracks "+suffix),100,-3,3);
244	>	h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#eta of rec tracks "+suffix),100,-3,3);
245		h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
246	<	h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-5,5);
246	>	h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-0.5,0.5);
247		h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
248
249		h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
250		h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
251	<	h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#Eta of rec tracks in PV"+suffix),100,-3,3);
251	>	h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#eta of rec tracks in PV"+suffix),100,-3,3);
252		h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
253		h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);
254		h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
#	Line 160 | Line 258 | PVStudy::PVStudy(const edm::ParameterSet
258		h["nrecPVDiff"]     = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
259		h["nTrkPVDiff"]     = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
260		h["nTrkPVRelDiff"]  = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
261	<
261	>
262	>	// Histograms on comparing the multi-vertices
263	>	// Difference in reconstructed vtx position
264	>	h["min_xsep"]   = subDir.make<TH1D>("min_xsep", "min x diff of primary and secondary pvtx",300,0,0.1);
265	>	h["min_xsep_sign"]   = subDir.make<TH1D>("min_xsep_sign", "min x diff in signf of primary and secondary pvtx",300,0,5);
266	>	h["min_ysep"]   = subDir.make<TH1D>("min_ysep", "min y diff of primary and secondary pvtx",300,0,0.1);
267	>	h["min_ysep_sign"]   = subDir.make<TH1D>("min_ysep_sign", "min y diff in signf of primary and secondary pvtx",300,0,5);
268	>	h["min_zsep"]   = subDir.make<TH1D>("min_zsep", "min z diff of primary and secondary pvtx",300,0,5);
269	>	h["min_zsep_sign"]   = subDir.make<TH1D>("min_zsep_sign", "min z diff in signf of primary and secondary pvtx",300,0,200);
270	>	// Difference in reconstructed vtx position
271	>	h["min_ntrksep"]   = subDir.make<TH1D>("min_ntrksep", "min nTrk diff of primary and secondary pvtx",201,-50.5,150.5);
272	>	h["min_sumpt2sep"]   = subDir.make<TH1D>("min_sumpt2sep", "min sumpt2 diff of primary and secondary pvtx",300,0,10000);
273	>
274		//Book histograms sensitive to data/mc
275		for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
276		string suffix;
#	Line 239 | Line 349 | PVStudy::PVStudy(const edm::ParameterSet
349		h2["pullz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
350		} // End of   if(analyze_) {
351		suffix.clear();
352	<	} // End of Book histograms sensitive to data/mc
352	>	} // End of Book histograms sensitive to data/mc
353	>
354	>	// Book histograms about pixelVertices
355	>	h["trkdz_pxlpvtxdz"]   = subDir.make<TH1D>("trkdz_pxlpvtxdz", "(Track dz - pixelpvtx dz) in cm",300,-0.5,0.5);
356	>	h["trkdz_pxlpvtxdz_pxlpvtxdzerr"]   = subDir.make<TH1D>("trkdz_pxlpvtxdz_pxlpvtxdzerr", "|Track dz - pixelpvtx dz| / pxlpvtxdzErr",300,0,100);
357	>	h["trkdz_pxlpvtxdz_trkdzerr"]   = subDir.make<TH1D>("trkdz_pxlpvtxdz_trkdzerr", "|Track dz - pixelpvtx dz| / trkdzErr",300,0,50);
358	>	h["trkdzErr_pxlpvtx"]   = subDir.make<TH1D>("trkdzErr_pxlpvtxdz", "Track dzErr of leading pixelpvtx ",300,0,0.5);
359	>	h["trkdzErr_pvtx"]   = subDir.make<TH1D>("trkdzErr_pvtx", "Track dzErr of the leading pvtx ",300,0,0.5);
360	>	h["dzErr_pxlpvtx"]   = subDir.make<TH1D>("dzErr_pxlpvtx", "zError of the leading pvtx ",300,0,0.5);
361	>	// Compare offlinePrimaryVertices with pixelVertices
362	>	h["nrecPV_minus_nrecPxlPV"] =  subDir.make<TH1D>("nrecPV_minus_nrecPxlPV", "nrecPV_minus_nrecPxlPV",21,-10.5,10.5);
363	>	h["recxPV_minus_recxPxlPV"] =  subDir.make<TH1D>("recxPV_minus_recxPxlPV", "recxPV_minus_recxPxlPV",300,-0.02,0.02);
364	>	h["recyPV_minus_recyPxlPV"] =  subDir.make<TH1D>("recyPV_minus_recyPxlPV", "recyPV_minus_recyPxlPV",300,-0.02,0.02);
365	>	h["reczPV_minus_reczPxlPV"] =  subDir.make<TH1D>("reczPV_minus_reczPxlPV", "reczPV_minus_reczPxlPV",300,-0.1,0.1);
366	>
367		// Book MC only plots
368		if (!realData_) {
369		h["genPart_T"]      = subDir.make<TH1D>("genPart_T","t component of gen particles",300,-0.5,0.5);
#	Line 472 | Line 596 | PVStudy::analyze(const edm::Event& iEven
596		//=======================================================
597		// Initialize Root-tuple variables if needed
598		//=======================================================
599	<	if(saventuple_)
600	<	SetVarToZero();
599	>	//if(saventuple_)
600	>	SetVarToZero();
601
602		//=======================================================
603		// Track accessors
#	Line 488 | Line 612 | PVStudy::analyze(const edm::Event& iEven
612		} else {
613		cout << "trackCollection cannot be found -> using empty collection of same type." <<endl;
614		}
615	+
616		//splitTrackCollection1
617		static const reco::TrackCollection s_empty_splitTrackColl1;
618		const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
#	Line 498 | Line 623 | PVStudy::analyze(const edm::Event& iEven
623		} else {
624		cout << "splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
625		}
626	+
627		//splitTrackCollection2
628		static const reco::TrackCollection s_empty_splitTrackColl2;
629		const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
#	Line 508 | Line 634 | PVStudy::analyze(const edm::Event& iEven
634		} else {
635		cout << "splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
636		}
637	<
637	>
638	>
639	>	//=======================================================
640	>	// Fill trackparameters of the input tracks to pvtx fitter
641	>	//=======================================================
642	>	if(verbose_)
643	>	cout<<"Start filling track parameters of the input tracks to pvtx fitter."<<endl;
644	>	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype)
645	>	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
646	>	fillTrackHisto(trackColl, 0);
647	>	fillTrackHisto(splitTrackColl1, 1);
648	>	fillTrackHisto(splitTrackColl2, 2);
649	>	if(verbose_)
650	>	cout<<"End filling track parameters of the input tracks to pvtx fitter."<<endl;
651	>
652		//=======================================================
653		// PVTX accessors
654		//=======================================================
#	Line 522 | Line 662 | PVStudy::analyze(const edm::Event& iEven
662		} else {
663		cout << "vertexCollection cannot be found -> using empty collection of same type." <<endl;
664		}
665	+
666		//splitVertexCollection1
667		static const reco::VertexCollection s_empty_splitVertexColl1;
668		const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
#	Line 531 | Line 672 | PVStudy::analyze(const edm::Event& iEven
672		splitVertexColl1 = splitVertexCollection1Handle.product();
673		} else {
674		cout << "splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
675	<	}
675	>	}
676	>
677		//splitVertexCollection2
678		static const reco::VertexCollection s_empty_splitVertexColl2;
679		const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
#	Line 542 | Line 684 | PVStudy::analyze(const edm::Event& iEven
684		} else {
685		cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
686		}
687	<
688	<	if(verbose_) cout<<"Done accessing the track and vertex collections"<<endl;
687	>
688	>	if(verbose_) cout<<"End accessing the track and vertex collections"<<endl;
689
690		//=======================================================
691		// MC simvtx accessor
#	Line 564 | Line 706 | PVStudy::analyze(const edm::Event& iEven
706		}catch(const Exception&){
707		cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
708		}
709	<
709	>
710		//=======================================================
711	<	// Fill trackparameters of the input tracks to pvtx fitter
711	>	// GET pixelVertices
712		//=======================================================
713	<
714	<	// GeneralTracks
715	<	h["nTrk"]->Fill(trackColl->size());
716	<	for(TrackCollection::const_iterator itTrack = trackColl->begin();
717	<	itTrack != trackColl->end();
718	<	++itTrack) {
719	<	h["trkPt"]->Fill(itTrack->pt());
720	<	h["trkDxy"]->Fill(itTrack->dxy());
579	<	h["trkDz"]->Fill(itTrack->dz());
580	<	h["trkEta"]->Fill(itTrack->eta());
581	<	h["trkPhi"]->Fill(itTrack->phi());
713	>	static const reco::VertexCollection s_empty_pixelVertexColl;
714	>	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
715	>	edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
716	>	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
717	>	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
718	>	pixelVertexColl = pixelVertexCollectionHandle.product();
719	>	} else {
720	>	cout << "pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
721		}
722	<
723	<	//SplittedTracks1
724	<	h["nTrk1_spl"]->Fill(splitTrackColl1->size());
725	<	for(TrackCollection::const_iterator itTrack = splitTrackColl1->begin();
726	<	itTrack != splitTrackColl1->end();
727	<	++itTrack) {
728	<	h["trkPt1_spl"]->Fill(itTrack->pt());
729	<	h["trkDxy1_spl"]->Fill(itTrack->dxy());
730	<	h["trkDz1_spl"]->Fill(itTrack->dz());
731	<	h["trkEta1_spl"]->Fill(itTrack->eta());
732	<	h["trkPhi1_spl"]->Fill(itTrack->phi());
733	<	}
734	<	//SplittedTracks2
735	<	h["nTrk2_spl"]->Fill(splitTrackColl2->size());
736	<	for(TrackCollection::const_iterator itTrack = splitTrackColl2->begin();
737	<	itTrack != splitTrackColl2->end();
738	<	++itTrack) {
739	<	h["trkPt2_spl"]->Fill(itTrack->pt());
740	<	h["trkDxy2_spl"]->Fill(itTrack->dxy());
741	<	h["trkDz2_spl"]->Fill(itTrack->dz());
742	<	h["trkEta2_spl"]->Fill(itTrack->eta());
743	<	h["trkPhi2_spl"]->Fill(itTrack->phi());
722	>
723	>	//=======================================================
724	>	// Fill pixelVertices related histograms
725	>	//=======================================================
726	>	nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
727	>	if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
728	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
729	>	fillTrackHistoInPV(pixelVertexColl, 4, false, true);
730	>	h["dzErr_pxlpvtx"]->Fill( pixelVertexColl->begin()->zError());
731	>	// Get the dZ error of the tracks in the leading pixelVertexColl
732	>	for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
733	>	t!= (pixelVertexColl->begin())->tracks_end(); t++) {
734	>
735	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
736	>	//        h["trkPtPV"]->Fill(0.);
737	>	}
738	>	else
739	>	h["trkdzErr_pxlpvtx"]->Fill((**t).dzError());
740	>	}
741	>	// Fill the track->dz() in the PV difference from first pixelpvtx
742	>	for(reco::Vertex::trackRef_iterator t = (vertexColl->begin())->tracks_begin();
743	>	t!= (vertexColl->begin())->tracks_end(); t++) {
744	>	// illegal charge
745	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
746	>	//        h["trkPtPV"]->Fill(0.);
747	>	}
748	>	else {
749	>	if(pixelVertexColl->size()>0) {
750	>	h["trkdz_pxlpvtxdz"]->Fill((**t).dz() -  pixelVertexColl->begin()->z());
751	>	h["trkdz_pxlpvtxdz_pxlpvtxdzerr"]->Fill(fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
752	>	h["trkdz_pxlpvtxdz_trkdzerr"]->Fill(fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
753	>	h["trkdzErr_pvtx"]->Fill((**t).dzError());
754	>	}
755	>	}
756	>	}
757		}
606	–
607	–	if(verbose_)
608	–	cout<<"Done filling track parameters of the input tracks to pvtx fitter."<<endl;
758
759		//=======================================================
760	<	// Fill reconstructed vertices
760	>	// Fill number of reconstructed vertices
761		//=======================================================
762	+
763		if(verbose_)  {
764		cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
765		printRecVtxs(vertexCollectionHandle);
#	Line 617 | Line 767 | PVStudy::analyze(const edm::Event& iEven
767		printRecVtxs(splitVertexCollection1Handle);
768		cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
769		printRecVtxs(splitVertexCollection2Handle);
770	+	cout<<"Start filling pvtx parameters."<<endl;
771		}
772
773	<	h["nrecPV"]->Fill(vertexColl->size());
774	<	h["nrecPV1_spl"]->Fill(splitVertexColl1->size());
775	<	h["nrecPV2_spl"]->Fill(splitVertexColl2->size());
776	<	h["nrecPVDiff"]->Fill(double(splitVertexColl1->size())-double(splitVertexColl2->size()));
777	<	//Fill the tree variables
778	<	nrecPV_ = vertexColl->size();
779	<	nrecPV1_spl_ = splitVertexColl1->size();
780	<	nrecPV2_spl_ = splitVertexColl2->size();
773	>	nrecPV_ = int (vertexColl->size());
774	>	nrecPV1_spl_ = int (splitVertexColl1->size());
775	>	nrecPV2_spl_ = int (splitVertexColl2->size());
776	>
777	>	h["nrecPV"]->Fill(nrecPV_);
778	>	h["nrecPV1_spl"]->Fill(nrecPV1_spl_);
779	>	h["nrecPV2_spl"]->Fill(nrecPV2_spl_);
780	>	h["nrecPVDiff"]->Fill(double(nrecPV1_spl_)-double(nrecPV2_spl_));
781
782	<	// Fill the track paramters for the vertexColl
783	<	for(reco::VertexCollection::const_iterator v=vertexColl->begin();
784	<	v!=vertexColl->end(); ++v){
785	<	try {
786	<	for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
787	<	t!=v->tracks_end(); t++) {
788	<	// illegal charge
789	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
790	<	//      h["trkPtPV"]->Fill(0.);
791	<	}
792	<	else {
793	<	h["trkPtPV"]->Fill((**t).pt());
794	<	h["trkDxyPV"]->Fill((**t).dxy());
795	<	h["trkDzPV"]->Fill((**t).dz());
796	<	h["trkEtaPV"]->Fill((**t).eta());
797	<	h["trkPhiPV"]->Fill((**t).phi());
798	<	}
799	<	}
800	<	}
801	<	catch (...) {
802	<	// exception thrown when trying to use linked track
803	<	//          h["trkPtPV"]->Fill(0.);
804	<	}
782	>	//=================================================================
783	>	// Fill track parameter ntuple/hist for tracks used in recoVertices
784	>	//=================================================================
785	>
786	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
787	>	if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
788	>	fillTrackHistoInPV(vertexColl, 0, true, true);
789	>
790	>	if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
791	>	fillTrackHistoInPV(splitVertexColl1, 1, false, true);
792	>
793	>	if(splitVertexColl1->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
794	>	fillTrackHistoInPV(splitVertexColl2, 2, false, true);
795	>
796	>
797	>	//=======================================================
798	>	// Compare offlinePrimaryVertices with pixelVertices
799	>	//=======================================================
800	>	if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
801	>	&& (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {
802	>	h["nrecPV_minus_nrecPxlPV"]->Fill( double (nrecPV_ - nrecPV_pxlpvtx_));
803	>	// difference in reconstructed position of the leading pvtx
804	>	//cout<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;
805	>	//cout<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
806	>	h["recxPV_minus_recxPxlPV"]->Fill (recx_[0] - recx_pxlpvtx_[0]);
807	>	h["recyPV_minus_recyPxlPV"]->Fill (recy_[0] - recy_pxlpvtx_[0]);
808	>	h["reczPV_minus_reczPxlPV"]->Fill (recz_[0] - recz_pxlpvtx_[0]);
809		}
810
811	+
812	+
813	+	//==========================================================
814	+	// Fill secondary/primary min separations for multi-vertices
815	+	//==========================================================
816
817	<	// == ignore multi-vertices (for now)
818	<	//take single pvt collection and compare
819	<	if(splitVertexColl1->size() == 1 && splitVertexColl2->size() == 1 ) {
820	<	const reco::Vertex recvtx1 = *(splitVertexColl1->begin());
821	<	const reco::Vertex recvtx2 = *(splitVertexColl2->begin());
822	<	if (matchVertex(recvtx1, recvtx2))  {
823	<	int nTrkPV1 = recvtx1.tracksSize();
824	<	int nTrkPV2 = recvtx2.tracksSize();
825	<	h["nTrkPV1_spl"]->Fill(nTrkPV1);
826	<	h["nTrkPV2_spl"]->Fill(nTrkPV2);
827	<	double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
828	<	h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
829	<	h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
830	<
831	<	// Fill the track paramters for the splitVertexColl1
832	<	try {
833	<	for(reco::Vertex::trackRef_iterator t = recvtx1.tracks_begin();
834	<	t!=recvtx1.tracks_end(); t++) {
835	<	// illegal charge
836	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
837	<	//    h["trkPtPV"]->Fill(0.);
838	<	}
839	<	else {
840	<	h["trkPtPV1_spl"]->Fill((**t).pt());
841	<	h["trkDxyPV1_spl"]->Fill((**t).dxy());
842	<	h["trkDzPV1_spl"]->Fill((**t).dz());
843	<	h["trkEtaPV1_spl"]->Fill((**t).eta());
844	<	h["trkPhiPV1_spl"]->Fill((**t).phi());
845	<	}
846	<	}
817	>	if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {
818	>
819	>	double min_xsep = 9999.0;
820	>	double min_ysep = 9999.0;
821	>	double min_zsep = 9999.0;
822	>	double min_xsep_sign = 9999.0;
823	>	double min_ysep_sign = 9999.0;
824	>	double min_zsep_sign = 9999.0;
825	>	double min_ntrksep = 9999.0;
826	>	double min_sumpt2sep = 9999999.0;
827	>
828	>	if(verbose_) cout<<"leading pvtx z = "<< vertexColl->begin()->z() <<endl;
829	>
830	>	// Looping through the secondary vertices to find the mininum separation to the primary
831	>	for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
832	>	v!=vertexColl->end(); ++v) {
833	>	if(v->isValid() && ! v->isFake()) {
834	>	// xsep
835	>	if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
836	>	min_xsep = fabs(v->x()- vertexColl->begin()->x());
837	>	// ysep
838	>	if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
839	>	min_ysep = fabs(v->y()- vertexColl->begin()->y());
840	>	// zsep
841	>	if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
842	>	min_zsep = fabs(v->z()- vertexColl->begin()->z());
843	>	// xsep_sign
844	>	double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
845	>	if(xsep_sign < min_xsep_sign)
846	>	min_xsep_sign = xsep_sign;
847	>	// ysep_sign
848	>	double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
849	>	if(ysep_sign < min_ysep_sign)
850	>	min_ysep_sign = ysep_sign;
851	>	// zsep_sign
852	>	double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
853	>	if(zsep_sign < min_zsep_sign)
854	>	min_zsep_sign = zsep_sign;
855	>	// ntrksep
856	>	if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
857	>	min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
858	>	// sumpt2sep
859	>	if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
860	>	min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
861		}
862	<	catch (...) {
863	<	// exception thrown when trying to use linked track
864	<	//        h["trkPtPV"]->Fill(0.);
862	>	}
863	>	h["min_xsep"]->Fill(min_xsep);
864	>	h["min_ysep"]->Fill(min_ysep);
865	>	h["min_zsep"]->Fill(min_zsep);
866	>	h["min_xsep_sign"]->Fill(min_xsep_sign);
867	>	h["min_ysep_sign"]->Fill(min_ysep_sign);
868	>	h["min_zsep_sign"]->Fill(min_zsep_sign);
869	>	h["min_ntrksep"]->Fill(min_ntrksep);
870	>	h["min_sumpt2sep"]->Fill(min_sumpt2sep);
871	>
872	>	min_zsep_ = min_zsep;
873	>	min_ntrksep_ = min_ntrksep;
874	>	min_sumpt2sep_ = min_sumpt2sep;
875	>
876	>
877	>	} // End of  if(vertexColl->size()>1) {
878	>
879	>	if(verbose_)
880	>	cout<<"End filling pvtx parameters."<<endl;
881	>
882	>	//=======================================================
883	>	//           PrimaryVertex Matching
884	>	// 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
885	>	// 2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
886	>	// Assume the first match is the primary vertex,
887	>	// since vertexColl are sorted in decreasing order of sum(pT)
888	>	//=======================================================
889	>
890	>	if(verbose_)   cout<<"PVStudy::analyze(): matching pvtxs "<<endl;
891	>	reco::VertexCollection s_empty_matchedVertexColl1;
892	>	reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
893	>	matchedVertexColl1->reserve(splitVertexColl1->size());
894	>	reco::VertexCollection s_empty_matchedVertexColl2;
895	>	reco::VertexCollection *matchedVertexColl2 = &s_empty_matchedVertexColl2;
896	>	matchedVertexColl2->reserve(splitVertexColl2->size());
897	>
898	>	bool stopmatching = false;
899	>	for (size_t ivtx = 0; ivtx<splitVertexCollection1Handle->size() && !stopmatching; ++ivtx) {
900	>	reco::VertexRef recvtx1(splitVertexCollection1Handle, ivtx);
901	>	if( !(recvtx1->isValid()) || recvtx1->isFake()) break;
902	>	for (size_t jvtx = ivtx; jvtx < splitVertexCollection2Handle->size(); ++jvtx) {
903	>	//------------------------------------------------------------------------
904	>	//== If only considering matching the first vertex of the splitVertexColl
905	>	//------------------------------------------------------------------------
906	>	if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
907	>	reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
908	>	if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
909	>	if(verbose_) {
910	>	cout<<"Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
911	>	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
912	>	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
913		}
914	<
915	<	// Fill the track paramters for the splitVertexColl2
916	<	try {
917	<	for(reco::Vertex::trackRef_iterator t = recvtx2.tracks_begin();
918	<	t!=recvtx2.tracks_end(); t++) {
919	<	// illegal charge
920	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
921	<	//    h["trkPtPV"]->Fill(0.);
922	<	}
923	<	else {
924	<	h["trkPtPV2_spl"]->Fill((**t).pt());
925	<	h["trkDxyPV2_spl"]->Fill((**t).dxy());
704	<	h["trkDzPV2_spl"]->Fill((**t).dz());
705	<	h["trkEtaPV2_spl"]->Fill((**t).eta());
706	<	h["trkPhiPV2_spl"]->Fill((**t).phi());
914	>	if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
915	>	if(verbose_) {
916	>	cout<<"The two splitted vertices match in Z. "<<endl;
917	>	}
918	>	int nTrkPV1 = recvtx1->tracksSize();
919	>	int nTrkPV2 = recvtx2->tracksSize();
920	>	double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
921	>	h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
922	>	h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
923	>	if(fabs(ntrkreldiff)<ntrkdiffcut_) {
924	>	if(verbose_) {
925	>	cout<<"(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
926		}
927	+	matchedVertexColl1->push_back(*recvtx1);
928	+	matchedVertexColl2->push_back(*recvtx2);
929	+	stopmatching = true; // stop the matching when the first match is found!
930	+	break;
931		}
932	+	else
933	+	cout<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
934		}
935	<	catch (...) {
936	<	// exception thrown when trying to use linked track
937	<	//        h["trkPtPV"]->Fill(0.);
935	>	else {
936	>	cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
937	>	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
938	>	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
939		}
714	–
715	–	if(verbose_)
716	–	cout<<"Done filling track parameters in reconstruced vertices."<<endl;
717	–
718	–	if(fabs(ntrkreldiff)<ntrkdiffcut_) {
719	–	fres_[0] = (recvtx1.x()-recvtx2.x())/sqrt(2.0);
720	–	fres_[1] = (recvtx1.y()-recvtx2.y())/sqrt(2.0);
721	–	fres_[2] = (recvtx1.z()-recvtx2.z())/sqrt(2.0);
722	–	ferror_[0] = sqrt(pow(recvtx1.xError(),2)+pow(recvtx2.xError(),2))/sqrt(2);
723	–	ferror_[1] = sqrt(pow(recvtx1.yError(),2)+pow(recvtx2.yError(),2))/sqrt(2);
724	–	ferror_[2] = sqrt(pow(recvtx1.zError(),2)+pow(recvtx2.zError(),2))/sqrt(2);
725	–	fntrk_ =  (nTrkPV1+nTrkPV2)/2;
726	–
727	–	h["deltax"]->Fill( fres_[0] );
728	–	h["deltay"]->Fill( fres_[1] );
729	–	h["deltaz"]->Fill( fres_[2] );
730	–	h["pullx"]->Fill( fres_[0]/ferror_[0]);
731	–	h["pully"]->Fill( fres_[1]/ferror_[1]);
732	–	h["pullz"]->Fill( fres_[2]/ferror_[2]);
733	–	h["errPVx"]->Fill( ferror_[0] );
734	–	h["errPVy"]->Fill( ferror_[1] );
735	–	h["errPVz"]->Fill( ferror_[2] );
736	–	pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
737	–	error(ferror_[0], ferror_[1], ferror_[2]),
738	–	fntrk_));
739	–	// Fill histo according to its track multiplicity, set datamode = 0
740	–	fillHisto(res(fres_[0], fres_[1], fres_[2]),
741	–	error(ferror_[0], ferror_[1], ferror_[2]),
742	–	fntrk_,0);
743	–
744	–	ftree_->Fill();
745	–	} // End of  if(fabs(ntrkreldiff)<ntrkdiffcut_) {
746	–	} // End of  if (matchVertex(recvtx1, recvtx2))  {
747	–	else   {
748	–	cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
749	–	cout<<"recvtx1.z() = " << recvtx1.z() <<"; recvtx2.z() = " <<  recvtx2.z() <<endl;
940		}
941	<	} // End of fill histograms for two-vertices Method
942	<
943	<	if(verbose_)
944	<	cout<<"Done filling res/pull with two-vertices method"<<endl;
941	>	}
942	>	if(verbose_)     cout<<"PVStudy::analyze(): End matching pvtxs"<<endl;
943	>
944	>	//=======================================================
945	>	//  Analyze matchedVertexColls
946	>	//=======================================================
947	>	if(verbose_) {
948	>	cout<<"Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
949	>	}
950	>
951	>	// Compare the reconstructed vertex position and calculate resolution/pulls
952	>	if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
953	>	//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
954	>	fillTrackHistoInPV(matchedVertexColl1, 1, true, false);
955	>	fillTrackHistoInPV(matchedVertexColl2, 2, true, false);
956	>
957	>	reco::VertexCollection::const_iterator v1;
958	>	reco::VertexCollection::const_iterator v2;
959	>	nrecPV_twovtx_ = 0;
960	>	for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
961	>	v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
962	>	++v1, ++v2) {
963	>	fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
964	>	fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
965	>	fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
966	>	ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
967	>	ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
968	>	ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
969	>	fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
970	>
971	>	h["deltax"]->Fill( fres_[0] );
972	>	h["deltay"]->Fill( fres_[1] );
973	>	h["deltaz"]->Fill( fres_[2] );
974	>	h["pullx"]->Fill( fres_[0]/ferror_[0]);
975	>	h["pully"]->Fill( fres_[1]/ferror_[1]);
976	>	h["pullz"]->Fill( fres_[2]/ferror_[2]);
977	>	h["errPVx"]->Fill( ferror_[0] );
978	>	h["errPVy"]->Fill( ferror_[1] );
979	>	h["errPVz"]->Fill( ferror_[2] );
980	>	pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
981	>	error(ferror_[0], ferror_[1], ferror_[2]),
982	>	fntrk_));
983	>	// Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
984	>	fillHisto(res(fres_[0], fres_[1], fres_[2]),
985	>	error(ferror_[0], ferror_[1], ferror_[2]),
986	>	fntrk_,0);
987	>	// Fill the ntuple variables
988	>	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
989	>	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
990	>	sumptsq1_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v1);
991	>	sumptsq2_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v2);
992	>	nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
993	>	deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
994	>	deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
995	>	deltaz_twovtx_[nrecPV_twovtx_] = fres_[2];
996	>	errx_twovtx_[nrecPV_twovtx_] = ferror_[0];
997	>	erry_twovtx_[nrecPV_twovtx_] = ferror_[1];
998	>	errz_twovtx_[nrecPV_twovtx_] = ferror_[2];
999	>	pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
1000	>	pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
1001	>	pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];
1002	>	nrecPV_twovtx_++;
1003	>	} // End of analyzing res/pull
1004	>	} // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
1005	>	else
1006	>	cout<<"WARNING: Cannot find matching pvtxs"<<endl;
1007	>
1008	>	if(verbose_)     cout<<"End analyzing the matchedVertexColl"<<endl;
1009
1010		//=======================================================
1011		// Fill simulated vertices
#	Line 778 | Line 1032 | PVStudy::analyze(const edm::Event& iEven
1032		simpv=getSimPVs(evtMC,"");
1033		//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
1034		h["nsimPV"]->Fill(simpv.size());
1035	<	int nsimtrk=0;
1036	<	int nrectrk=0;
1037	<
1035	>	nsimPV_ = simpv.size();
1036	>
1037	>	int isimpv = 0;
1038		for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
1039	<	vsim!=simpv.end(); vsim++){
1040	<	nsimtrk+=vsim->nGenTrk;
1041	<	// Loop over the datamodes, for MC, fill only 1, 2, 3
1042	<	for (int i=1;i<=3;i++) {
1043	<	static const reco::VertexCollection s_empty_vtxColl;
1044	<	const reco::VertexCollection *vtxColl = &s_empty_vtxColl;
1045	<	std::string suffix;
1046	<	// Set the vertexColl and histogram suffix according to datamode
1047	<	if (i == 1) {
1048	<	vtxColl = splitVertexColl1;
1049	<	suffix = "_spl1_mct";
796	<	}
797	<	if (i == 2) {
798	<	vtxColl =  splitVertexColl2;
799	<	suffix = "_spl2_mct";
800	<	}
801	<	if (i == 3) {
802	<	vtxColl =  vertexColl;
803	<	suffix = "_mct";
804	<	}
805	<	//========================================================
806	<	//  look for a matching reconstructed vertex in vertexColl
807	<	//========================================================
808	<	vsim->recVtx=NULL;
809	<	for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
810	<	vrec!=vtxColl->end(); ++vrec){
811	<	nrectrk=vrec->tracksSize();
812	<	if(verbose_){
813	<	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
814	<	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
815	<	}
816	<	if ( matchVertex(*vsim,*vrec) ){
817	<	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
818	<	if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
819	<	|| (!vsim->recVtx) )
820	<	vsim->recVtx=&(*vrec);
821	<	}// End of finding the matched reco_vertex
822	<
823	<	if (vsim->recVtx){
824	<	if(verbose_) {
825	<	cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
826	<	}
827	<	double fres_mct[3];
828	<	double ferror_mct[3];
829	<	int fntrk_mct;
830	<
831	<	fres_mct[0] = vsim->recVtx->x()-vsim->x;
832	<	fres_mct[1] = vsim->recVtx->y()-vsim->y;
833	<	fres_mct[2] = vsim->recVtx->z()-vsim->z;
834	<	ferror_mct[0] = vsim->recVtx->xError();
835	<	ferror_mct[1] = vsim->recVtx->yError();
836	<	ferror_mct[2] = vsim->recVtx->zError();
837	<	fntrk_mct = nrectrk;
838	<
839	<	//Fill the values for variables in ftree_
840	<	if(i == 1) {
841	<	for(int j = 0;j<3;j++)
842	<	fres_spl1_mct_[j] = fres_mct[j];
843	<	fntrk_spl1_mct_ =  fntrk_mct;
844	<	}
845	<	if(i == 2) {
846	<	for(int j = 0;j<3;j++)
847	<	fres_spl2_mct_[j] = fres_mct[j];
848	<	fntrk_spl2_mct_ =  fntrk_mct;
849	<	}
850	<	if(i == 3) {
851	<	for(int j=0;j<3;j++)
852	<	fres_mct_[j] = fres_mct[j];
853	<	fntrk_mct_ =  fntrk_mct;
854	<	}
855	<
856	<	h["deltax"+suffix]->Fill( fres_mct[0] );
857	<	h["deltay"+suffix]->Fill( fres_mct[1] );
858	<	h["deltaz"+suffix]->Fill( fres_mct[2] );
859	<	h["pullx"+suffix]->Fill( fres_mct[0]/ferror_mct[0] );
860	<	h["pully"+suffix]->Fill( fres_mct[1]/ferror_mct[1] );
861	<	h["pullz"+suffix]->Fill( fres_mct[2]/ferror_mct[2] );
862	<	h["errPVx"+suffix]->Fill( ferror_mct[0] );
863	<	h["errPVy"+suffix]->Fill( ferror_mct[1] );
864	<	h["errPVz"+suffix]->Fill( ferror_mct[2] );
865	<	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
866	<	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
867	<	fntrk_mct));
868	<	// Fill histo according to its track multiplicity
869	<	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
870	<	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
871	<	fntrk_mct, i);
872	<	ftree_->Fill();
873	<	suffix.clear();
874	<	}
875	<	else {  // no rec vertex found for this simvertex
876	<	if(verbose_)
877	<	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
878	<	}
879	<	}
880	<	} // === End of Looping through vertexColl
1039	>	vsim!=simpv.end(); vsim++, isimpv++){
1040	>	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
1041	>	simx_[isimpv] = vsim->x;
1042	>	simy_[isimpv] = vsim->y;
1043	>	simz_[isimpv] = vsim->y;
1044	>	simptsq_[isimpv] = vsim->ptsq;
1045	>
1046	>	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
1047	>	fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
1048	>	fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
1049	>	fillMCHisto(vsim, isimpv, vertexColl, 3);
1050		}
1051		} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
1052	<	} // End of Filling MC variables
1052	>	} // End of   if (!realData_) {
1053	>
1054	>	// Finally fill the ftree_
1055	>	if(saventuple_) {
1056	>	ftree_->Fill();
1057	>	pvtxtree_->Fill();
1058	>	}
1059		}
1060
1061		// ------------ method called once each job just before starting event loop  ------------
#	Line 907 | Line 1082 | PVStudy::endJob() {
1082		case 3: suffix = "_mct";
1083		break;
1084		}
1085	<	for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1086	<	if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
1087	<	PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
913	<	if (analyze_) {
1085	>	if(analyze_) {
1086	>	for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1087	>	PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
1088		if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
1089		if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
1090		if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
#	Line 918 | Line 1092 | PVStudy::endJob() {
1092		if ( ipv.pull_.y() > 0 ) h2["pully"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
1093		if ( ipv.pull_.z() > 0 ) h2["pullz"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.z(), 1.);
1094		}
1095	<	} // End of  for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1095	>	}
1096		suffix.clear();
1097		}
1098	<	file_->cd();
1099	<	ftree_->Write();
1100	<	file_->Close();
1098	>	if(saventuple_) {
1099	>	file_->cd();
1100	>	ftree_->Write();
1101	>	pvtxtree_->Write();
1102	>	file_->Close();
1103	>	}
1104		}
1105
1106
#	Line 938 | Line 1115 | bool PVStudy::matchVertex(const PVStudy:
1115		}
1116
1117		// Match two reconstructed vertices
1118	<	bool PVStudy::matchVertex( const reco::Vertex & vrec1,
1119	<	const reco::Vertex & vrec2)
1118	>	bool PVStudy::matchVertex( const reco::VertexRef & vrec1,
1119	>	const reco::VertexRef & vrec2,
1120	>	int zsigncut)
1121		{
1122	<	return (fabs(vrec1.z()-vrec2.z())<0.0500); // =500um
1122	>	double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
1123	>	if(verbose_)
1124	>	cout<<"The matching criteria in z is "<<cut<<endl;
1125	>	return (fabs(vrec1->z()-vrec2->z())<cut);
1126		}
1127
1128
#	Line 1141 | Line 1322 | void PVStudy::fit(TH1 *&hdist, double fi
1322		fitPars[0] = -999;
1323		}
1324
1325	<	void PVStudy::SetVarToZero() {
1325	>	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype) {
1326	>	string suffix;
1327	>	suffix = ""; // for unsplit tracks
1328	>	if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1329	>	if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1330
1331	<	// number of reconstructed vertices
1332	<	nrecPV_ = 0;
1333	<	nrecPV1_spl_ = 0;
1334	<	nrecPV2_spl_ = 0;
1335	<	// number of tracks in the vertex
1331	>	h["nTrk"+suffix]->Fill(trackColl->size());
1332	>	for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1333	>	itTrack != trackColl->end();
1334	>	++itTrack) {
1335	>	h["trkPt"+suffix]->Fill(itTrack->pt());
1336	>	h["trkDxy"+suffix]->Fill(itTrack->dxy());
1337	>	h["trkDz"+suffix]->Fill(itTrack->dz());
1338	>	h["trkEta"+suffix]->Fill(itTrack->eta());
1339	>	h["trkPhi"+suffix]->Fill(itTrack->phi());
1340	>	}
1341	>	}
1342	>
1343	>	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
1344	>	string suffix;
1345	>	suffix = ""; // for unsplit tracks
1346	>	if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1347	>	if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1348	>	int ivtx = 0;
1349	>	for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1350	>	v!=vertexColl->end(); ++v, ++ivtx) {
1351	>	if(fillNtuple) {
1352	>	if(datatype == 0) {
1353	>	nTrkPV_[ivtx] = v->tracksSize();
1354	>	sumptsq_[ivtx] = sumPtSquared(*v);
1355	>	isValid_[ivtx] = v->isValid();
1356	>	isFake_[ivtx] = v->isFake();
1357	>	recx_[ivtx] = v->x();
1358	>	recy_[ivtx] = v->y();
1359	>	recz_[ivtx] = v->z();
1360	>	recx_err_[ivtx] = v->xError();
1361	>	recy_err_[ivtx] = v->yError();
1362	>	recz_err_[ivtx] = v->zError();
1363	>
1364	>
1365	>	}
1366	>	if(datatype == 1) {
1367	>	nTrkPV1_spl_[ivtx] = v->tracksSize();
1368	>	sumptsq1_spl_[ivtx] = sumPtSquared(*v);
1369	>	isValid1_spl_[ivtx] = v->isValid();
1370	>	isFake1_spl_[ivtx] = v->isFake();
1371	>	recx1_spl_[ivtx] = v->x();
1372	>	recy1_spl_[ivtx] = v->y();
1373	>	recz1_spl_[ivtx] = v->z();
1374	>	recx1_err_spl_[ivtx] = v->xError();
1375	>	recy1_err_spl_[ivtx] = v->yError();
1376	>	recz1_err_spl_[ivtx] = v->zError();
1377	>
1378	>	}
1379	>	if(datatype == 2) {
1380	>	nTrkPV2_spl_[ivtx] = v->tracksSize();
1381	>	sumptsq2_spl_[ivtx] = sumPtSquared(*v);
1382	>	isValid2_spl_[ivtx] = v->isValid();
1383	>	isFake2_spl_[ivtx] = v->isFake();
1384	>	recx2_spl_[ivtx] = v->x();
1385	>	recy2_spl_[ivtx] = v->y();
1386	>	recz2_spl_[ivtx] = v->z();
1387	>	recx2_err_spl_[ivtx] = v->xError();
1388	>	recy2_err_spl_[ivtx] = v->yError();
1389	>	recz2_err_spl_[ivtx] = v->zError();
1390	>	}
1391	>	if(datatype == 4) { // for pixelVertices
1392	>	nTrkPV_pxlpvtx_[ivtx] = v->tracksSize();
1393	>	sumptsq_pxlpvtx_[ivtx] = sumPtSquared(*v);
1394	>	isValid_pxlpvtx_[ivtx] = v->isValid();
1395	>	isFake_pxlpvtx_[ivtx] = v->isFake();
1396	>	recx_pxlpvtx_[ivtx] = v->x();
1397	>	recy_pxlpvtx_[ivtx] = v->y();
1398	>	recz_pxlpvtx_[ivtx] = v->z();
1399	>	recx_err_pxlpvtx_[ivtx] = v->xError();
1400	>	recy_err_pxlpvtx_[ivtx] = v->yError();
1401	>	recz_err_pxlpvtx_[ivtx] = v->zError();
1402	>	}
1403	>	}
1404	>	}
1405	>	// For histogram only fill the leading pvtx parameters
1406	>	if(fillHisto) {
1407	>	h["nTrkPV"+suffix]->Fill(vertexColl->begin()->tracksSize());
1408	>	try {
1409	>	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1410	>	t!=vertexColl->begin()->tracks_end(); t++) {
1411	>	// illegal charge
1412	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1413	>	//      h["trkPtPV"]->Fill(0.);
1414	>	}
1415	>	else {
1416	>	h["trkPtPV"+suffix]->Fill((**t).pt());
1417	>	h["trkDxyPV"+suffix]->Fill((**t).dxy());
1418	>	h["trkDzPV"+suffix]->Fill((**t).dz());
1419	>	h["trkEtaPV"+suffix]->Fill((**t).eta());
1420	>	h["trkPhiPV"+suffix]->Fill((**t).phi());
1421	>	}
1422	>	}
1423	>	}
1424	>	catch (...) {
1425	>	// exception thrown when trying to use linked track
1426	>	//          h["trkPtPV"]->Fill(0.);
1427	>	}
1428	>	}
1429	>
1430	>	}
1431	>
1432	>	void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
1433	>	{
1434	>	std::string suffix;
1435	>	// Set the vertexColl and histogram suffix according to datamode
1436	>	if (datatype == 1) {
1437	>	suffix = "_spl1_mct";
1438	>	nrecPV_spl1_mct_ = 0;
1439	>	}
1440	>	if (datatype == 2) {
1441	>	suffix = "_spl2_mct";
1442	>	nrecPV_spl2_mct_ = 0;
1443	>	}
1444	>	if (datatype == 3) {
1445	>	suffix = "_mct";
1446	>	nrecPV_mct_ = 0;
1447	>	}
1448	>
1449	>	//========================================================
1450	>	//  look for a matching reconstructed vertex in vertexColl
1451	>	//========================================================
1452	>
1453	>	for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1454	>	vrec!=vtxColl->end(); ++vrec){
1455	>	int nrectrk = vrec->tracksSize();
1456	>	vsim->recVtx=NULL;
1457	>
1458	>	if(verbose_){
1459	>	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1460	>	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1461	>	}
1462	>	if ( matchVertex(*vsim,*vrec) ) {
1463	>	vsim->recVtx=&(*vrec);
1464	>
1465	>	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1466	>	//if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
1467	>	//|| (!vsim->recVtx) )
1468	>	//vsim->recVtx=&(*vrec);
1469	>	//}
1470	>	if(verbose_)
1471	>	cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1472	>
1473	>	double fres_mct[3];
1474	>	double ferror_mct[3];
1475	>
1476	>	fres_mct[0] = vsim->recVtx->x()-vsim->x;
1477	>	fres_mct[1] = vsim->recVtx->y()-vsim->y;
1478	>	fres_mct[2] = vsim->recVtx->z()-vsim->z;
1479	>	ferror_mct[0] = vsim->recVtx->xError();
1480	>	ferror_mct[1] = vsim->recVtx->yError();
1481	>	ferror_mct[2] = vsim->recVtx->zError();
1482	>
1483	>	h["deltax"+suffix]->Fill( fres_mct[0] );
1484	>	h["deltay"+suffix]->Fill( fres_mct[1] );
1485	>	h["deltaz"+suffix]->Fill( fres_mct[2] );
1486	>	h["pullx"+suffix]->Fill( fres_mct[0]/ferror_mct[0] );
1487	>	h["pully"+suffix]->Fill( fres_mct[1]/ferror_mct[1] );
1488	>	h["pullz"+suffix]->Fill( fres_mct[2]/ferror_mct[2] );
1489	>	h["errPVx"+suffix]->Fill( ferror_mct[0] );
1490	>	h["errPVy"+suffix]->Fill( ferror_mct[1] );
1491	>	h["errPVz"+suffix]->Fill( ferror_mct[2] );
1492	>	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1493	>	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1494	>	nrectrk));
1495	>	// Fill histo according to its track multiplicity
1496	>	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1497	>	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1498	>	nrectrk, datatype);
1499	>
1500	>	if(saventuple_) {
1501	>	//Fill the values for variables in ftree_
1502	>	if(datatype == 1) {
1503	>	nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1504	>	deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1505	>	deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1506	>	deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1507	>	pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
1508	>	pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1509	>	pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1510	>	nrecPV_spl1_mct_++;
1511	>	}
1512	>
1513	>	if(datatype == 2) {
1514	>	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1515	>	deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1516	>	deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1517	>	deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1518	>	pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
1519	>	pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1520	>	pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1521	>	nrecPV_spl2_mct_++;
1522	>	}
1523	>	if(datatype == 3) {
1524	>	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1525	>	deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1526	>	deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1527	>	deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1528	>	pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
1529	>	pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1530	>	pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1531	>	nrecPV_mct_++;
1532	>	}
1533	>	} // End of  if(saventuple_) {
1534	>	} //  if ( matchVertex(*vsim,*vrec) ) {
1535	>	else {  // no rec vertex found for this simvertex
1536	>	if(verbose_)
1537	>	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1538	>	}
1539	>	}
1540	>	}
1541	>
1542	>
1543	>	void PVStudy::SetVarToZero() {
1544	>	//Greg's variables
1545		fntrk_ = 0;
1152	–	fntrk_spl1_mct_ = 0;
1153	–	fntrk_spl2_mct_ = 0;
1154	–	fntrk_mct_ = 0;
1546		//pvtx position (x,y,z) residual and error
1547	<	for(int i = 0; i<3;i++)
1157	<	{
1547	>	for(int i = 0; i<3;i++)    {
1548		fres_[i] = 0;
1549		ferror_[i] = 0;
1550	<	fres_spl1_mct_[i] = 0;
1551	<	ferror_spl1_mct_[i] = 0;
1552	<	fres_spl2_mct_[i] = 0;
1553	<	ferror_spl2_mct_[i] = 0;
1554	<	fres_mct_[i] = 0;
1555	<	ferror_mct_[i] = 0;
1556	<	}
1550	>	}
1551	>
1552	>	//Yanyan's variables
1553	>	// Number of reconstructed vertices
1554	>	nrecPV_ = 0;
1555	>	nrecPV1_spl_ = 0;
1556	>	nrecPV2_spl_ = 0;
1557	>	nrecPV_twovtx_ = 0;
1558	>	nsimPV_ = 0;
1559	>
1560	>	nrecPV_mct_ = 0;
1561	>	nrecPV_spl1_mct_ = 0;
1562	>	nrecPV_spl2_mct_ = 0;
1563	>
1564	>	// Mininum separation between the secondary pvtxes and leading pvtx
1565	>	min_zsep_ = 9999.0;
1566	>	min_ntrksep_ = 9999.0;
1567	>	min_sumpt2sep_ = -9999.0;
1568	>
1569	>
1570	>	for (int i = 0; i < nMaxPVs_; i++) {
1571	>	// recoVertices with all tracks
1572	>	nTrkPV_[i] = 0; // Number of tracks in the pvtx
1573	>	sumptsq_[i] = 0;
1574	>	isValid_[i] = -1;
1575	>	isFake_[i] = -1;
1576	>	recx_[i] = 0;
1577	>	recy_[i] = 0;
1578	>	recz_[i] = 0;
1579	>	recx_err_[i] = 0;
1580	>	recy_err_[i] = 0;
1581	>	recz_err_[i] = 0;
1582	>
1583	>	// recoVertices with splitTrack1
1584	>	nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx
1585	>	sumptsq1_spl_[i] = 0;
1586	>	isValid1_spl_[i] = -1;
1587	>	isFake1_spl_[i] = -1;
1588	>	recx1_spl_[i] = 0;
1589	>	recy1_spl_[i] = 0;
1590	>	recz1_spl_[i] = 0;
1591	>	recx1_err_spl_[i] = 0;
1592	>	recy1_err_spl_[i] = 0;
1593	>	recz1_err_spl_[i] = 0;
1594	>
1595	>	// recoVertices with splitTrack2
1596	>	nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx
1597	>	sumptsq2_spl_[i] = 0;
1598	>	isValid2_spl_[i] = -1;
1599	>	isFake2_spl_[i] = -1;
1600	>	recx2_spl_[i] = 0;
1601	>	recy2_spl_[i] = 0;
1602	>	recz2_spl_[i] = 0;
1603	>	recx2_err_spl_[i] = 0;
1604	>	recy2_err_spl_[i] = 0;
1605	>	recz2_err_spl_[i] = 0;
1606	>
1607	>	//pixelVertices
1608	>	nTrkPV_pxlpvtx_[i] = 0; // Number of tracks in the pvtx
1609	>	sumptsq_pxlpvtx_[i] = 0;
1610	>	isValid_pxlpvtx_[i] = -1;
1611	>	isFake_pxlpvtx_[i] = -1;
1612	>	recx_pxlpvtx_[i] = 0;
1613	>	recy_pxlpvtx_[i] = 0;
1614	>	recz_pxlpvtx_[i] = 0;
1615	>	recx_err_pxlpvtx_[i] = 0;
1616	>	recy_err_pxlpvtx_[i] = 0;
1617	>	recz_err_pxlpvtx_[i] = 0;
1618	>
1619	>	// matched two-vertices
1620	>	nTrkPV1_spl_twovtx_[i] = 0;
1621	>	nTrkPV2_spl_twovtx_[i] = 0;
1622	>	nTrkPV_twovtx_[i] = 0;
1623	>	sumptsq1_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1624	>	sumptsq2_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1625	>	deltax_twovtx_[i] = 0;
1626	>	deltay_twovtx_[i] = 0;
1627	>	deltaz_twovtx_[i] = 0;
1628	>	errx_twovtx_[i] = 0;
1629	>	erry_twovtx_[i] = 0;
1630	>	errz_twovtx_[i] = 0;
1631	>	pullx_twovtx_[i] = 0;
1632	>	pully_twovtx_[i] = 0;
1633	>	pullz_twovtx_[i] = 0;
1634	>
1635	>	//simpv
1636	>	nsimTrkPV_[i] = 0;
1637	>	simx_[i] = 0;
1638	>	simy_[i] = 0;
1639	>	simz_[i] = 0;
1640	>	simptsq_[i] = 0;
1641	>
1642	>	// residual and pulls with mc truth required
1643	>	deltax_mct_[i] = 0;
1644	>	deltay_mct_[i] = 0;
1645	>	deltaz_mct_[i] = 0;
1646	>	pullx_mct_[i] = 0;
1647	>	pully_mct_[i] = 0;
1648	>	pullz_mct_[i] = 0;
1649	>
1650	>	deltax_spl1_mct_[i] = 0;
1651	>	deltay_spl1_mct_[i] = 0;
1652	>	deltaz_spl1_mct_[i] = 0;
1653	>	pullx_spl1_mct_[i] = 0;
1654	>	pully_spl1_mct_[i] = 0;
1655	>	pullz_spl1_mct_[i] = 0;
1656	>
1657	>	deltax_spl2_mct_[i] = 0;
1658	>	deltay_spl2_mct_[i] = 0;
1659	>	deltaz_spl2_mct_[i] = 0;
1660	>	pullx_spl2_mct_[i] = 0;
1661	>	pully_spl2_mct_[i] = 0;
1662	>	pullz_spl2_mct_[i] = 0;
1663	>	}
1664	>
1665		}
1666
1667	+	double PVStudy::sumPtSquared(const reco::Vertex & v)  {
1668	+	double sum = 0.;
1669	+	double pT;
1670	+	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1671	+	pT = (**it).pt();
1672	+	sum += pT*pT;
1673	+	}
1674	+	return sum;
1675	+	}
1676

Diff Legend

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