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.14 by jengbou, Fri Dec 4 20:48:04 2009 UTC vs.
Revision 1.19 by yygao, Sat Jan 30 18:27:09 2010 UTC

#	Line 12 | Line 12 | Implementation:
12		*/
13		//
14		// Original Author:  "Geng-yuan Jeng/UC Riverside"
15	<	//                   "Yanyan Gao/Fermilab ygao@fnal.gov"
15	>	//                    "Yanyan Gao/Fermilab ygao@fnal.gov"
16		//         Created:  Thu Aug 20 11:55:40 CDT 2009
17		// $Id$
18		//
#	Line 25 | Line 25 | Implementation:
25		#include <vector>
26		#include <iostream>
27		#include <sstream>
28	–	#include <algorithm>
28
29		// user include files
30		#include "FWCore/Framework/interface/Frameworkfwd.h"
#	Line 50 | Line 49 | Implementation:
49		#include <SimDataFormats/Vertex/interface/SimVertexContainer.h>
50		#include <SimDataFormats/Track/interface/SimTrack.h>
51		#include <SimDataFormats/Track/interface/SimTrackContainer.h>
52	+	// BeamSpot
53	+	#include "DataFormats/BeamSpot/interface/BeamSpot.h"
54
55		//root
56		#include <TROOT.h>
#	Line 60 | Line 61 | Implementation:
61		#include <TPad.h>
62
63		using namespace std;
64	+	typedef math::XYZTLorentzVectorF LorentzVector;
65	+	typedef math::XYZPoint Point;
66
67		PVStudy::PVStudy(const edm::ParameterSet& iConfig)
68		{
#	Line 72 | Line 75 | PVStudy::PVStudy(const edm::ParameterSet
75		splitTrackCollection2Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
76		vertexCollectionTag_       = iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection");
77		splitVertexCollection1Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
78	<	splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
76	<	pixelVertexCollectionTag_  = iConfig.getParameter<edm::InputTag>("pixelVertexCollectionTag");
78	>	splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
79		verbose_                   = iConfig.getUntrackedParameter<bool>("verbose",false);
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");
83	–	zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
86		nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
87		nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
88	<	histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
89	<
88	>	zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
89	>	bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
90	>
91		//now do what ever initialization is needed
92		pvinfo.clear();
93
91	–
94		// Specify the data mode vector
95		if(realData_) datamode.push_back(0);
96		else {
#	Line 97 | Line 99 | PVStudy::PVStudy(const edm::ParameterSet
99		datamode.push_back(2);
100		datamode.push_back(3);
101		}
102	<
101	<	// Create ntuple files if needed
102	>	// Create ntuple files if needed
103		if(saventuple_) {
104		file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
105		ftree_ = new TTree("mytree","mytree");
#	Line 109 | Line 110 | PVStudy::PVStudy(const edm::ParameterSet
110
111		// pvtxtree_ analyzing the pvtxs ootb
112		pvtxtree_ = new TTree("pvtxtree","pvtxtree");
112	–	//pvtx using all tracks
113	–	pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
114	–	pvtxtree_->Branch("nTrkPV",&nTrkPV_,"nTrkPV[nrecPV]/I");
115	–	pvtxtree_->Branch("sumptsq",&sumptsq_,"sumptsq[nrecPV]/D");
116	–	pvtxtree_->Branch("isValid",&isValid_,"isValid/I");
117	–	pvtxtree_->Branch("isFake",&isFake_,"isFake/I");
118	–	pvtxtree_->Branch("recx",&recx_,"recx[nrecPV]/D");
119	–	pvtxtree_->Branch("recy",&recy_,"recy[nrecPV]/D");
120	–	pvtxtree_->Branch("recz",&recz_,"recz[nrecPV]/D");
121	–	pvtxtree_->Branch("recx_err",&recx_err_,"recx_err[nrecPV]/D");
122	–	pvtxtree_->Branch("recy_err",&recy_err_,"recy_err[nrecPV]/D");
123	–	pvtxtree_->Branch("recz_err",&recz_err_,"recz_err[nrecPV]/D");
113
114	+	// Event information for the data
115	+	pvtxtree_->Branch("glob_runno",&glob_runno_,"glob_runno/I");
116	+	pvtxtree_->Branch("glob_evtno",&glob_evtno_,"glob_evtno/I");
117	+	pvtxtree_->Branch("glob_ls",&glob_ls_,"glob_ls/I");
118	+	pvtxtree_->Branch("glob_bx",&glob_bx_,"glob_bx/I");
119	+
120	+
121	+	pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
122	+	pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
123	+	pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
124	+	// Event level information
125		pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
126		pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
127		pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
128
129	<	//pvtx using splittracks1
130	<	pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
131	<	pvtxtree_->Branch("nTrkPV1_spl",&nTrkPV1_spl_,"nTrkPV1_spl[nrecPV1_spl]/I");
132	<	pvtxtree_->Branch("sumptsq1_spl",&sumptsq1_spl_,"sumptsq1_spl[nrecPV1_spl]/D");
133	<	pvtxtree_->Branch("isValid1_spl",&isValid1_spl_,"isValid1_spl/I");
134	<	pvtxtree_->Branch("isFake1_spl",&isFake1_spl_,"isFake1_spl/I");
135	<	pvtxtree_->Branch("recx1_spl",&recx1_spl_,"recx1_spl[nrecPV1_spl]/D");
136	<	pvtxtree_->Branch("recy1_spl",&recy1_spl_,"recy1_spl[nrecPV1_spl]/D");
137	<	pvtxtree_->Branch("recz1_spl",&recz1_spl_,"recz1_spl[nrecPV1_spl]/D");
138	<	pvtxtree_->Branch("recx1_err_spl",&recx1_err_spl_,"recx1_err_spl[nrecPV1_spl]/D");
139	<	pvtxtree_->Branch("recy1_err_spl",&recy1_err_spl_,"recy1_err_spl[nrecPV1_spl]/D");
140	<	pvtxtree_->Branch("recz1_err_spl",&recz1_err_spl_,"recz1_err_spl[nrecPV1_spl]/D");
141	<
142	<	//pvtx using splittracks2
143	<	pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
144	<	pvtxtree_->Branch("nTrkPV2_spl",&nTrkPV2_spl_,"nTrkPV2_spl[nrecPV2_spl]/I");
145	<	pvtxtree_->Branch("sumptsq2_spl",&sumptsq2_spl_,"sumptsq2_spl[nrecPV2_spl]/D");
146	<	pvtxtree_->Branch("isValid2_spl",&isValid2_spl_,"isValid2_spl/I");
147	<	pvtxtree_->Branch("isFake2_spl",&isFake2_spl_,"isFake2_spl/I");
148	<	pvtxtree_->Branch("recx2_spl",&recx2_spl_,"recx2_spl[nrecPV2_spl]/D");
149	<	pvtxtree_->Branch("recy2_spl",&recy2_spl_,"recy2_spl[nrecPV2_spl]/D");
150	<	pvtxtree_->Branch("recz2_spl",&recz2_spl_,"recz2_spl[nrecPV2_spl]/D");
151	<	pvtxtree_->Branch("recx2_err_spl",&recx2_err_spl_,"recx2_err_spl[nrecPV2_spl]/D");
152	<	pvtxtree_->Branch("recy2_err_spl",&recy2_err_spl_,"recy2_err_spl[nrecPV2_spl]/D");
153	<	pvtxtree_->Branch("recz2_err_spl",&recz2_err_spl_,"recz2_err_spl[nrecPV2_spl]/D");
154	<
155	<	//pixeVertices
156	<	pvtxtree_->Branch("nrecPV_pxlpvtx",&nrecPV_pxlpvtx_,"nrecPV_pxlpvtx/I");
157	<	pvtxtree_->Branch("nTrkPV_pxlpvtx",&nTrkPV_pxlpvtx_,"nTrkPV_pxlpvtx[nrecPV_pxlpvtx]/I");
158	<	pvtxtree_->Branch("sumptsq_pxlpvtx",&sumptsq_pxlpvtx_,"sumptsq_pxlpvtx[nrecPV_pxlpvtx]/D");
159	<	pvtxtree_->Branch("isValid_pxlpvtx",&isValid_pxlpvtx_,"isValid_pxlpvtx/I");
160	<	pvtxtree_->Branch("isFake_pxlpvtx",&isFake_pxlpvtx_,"isFake_pxlpvtx/I");
161	<	pvtxtree_->Branch("recx_pxlpvtx",&recx_pxlpvtx_,"recx_pxlpvtx[nrecPV_pxlpvtx]/D");
162	<	pvtxtree_->Branch("recy_pxlpvtx",&recy_pxlpvtx_,"recy_pxlpvtx[nrecPV_pxlpvtx]/D");
163	<	pvtxtree_->Branch("recz_pxlpvtx",&recz_pxlpvtx_,"recz_pxlpvtx[nrecPV_pxlpvtx]/D");
164	<	pvtxtree_->Branch("recx_err_pxlpvtx",&recx_err_pxlpvtx_,"recx_err_pxlpvtx[nrecPV_pxlpvtx]/D");
165	<	pvtxtree_->Branch("recy_err_pxlpvtx",&recy_err_pxlpvtx_,"recy_err_pxlpvtx[nrecPV_pxlpvtx]/D");
166	<	pvtxtree_->Branch("recz_err_pxlpvtx",&recz_err_pxlpvtx_,"recz_err_pxlpvtx[nrecPV_pxlpvtx]/D");
129	>	//Fill the variables in the twovtx pair (recvtx1, recvtx2)
130	>	// Information related to the analyzing the two-vertex method
131	>
132	>
133
168	–	//Fill the variables in the twovtx pair (recvtx1, recvtx2)
134		pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
170	–	pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
171	–	pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
172	–	pvtxtree_->Branch("sumptsq1_spl_twovtx",&sumptsq1_spl_twovtx_,"sumptsq1_spl_twovtx[nrecPV_twovtx]/D");
173	–	pvtxtree_->Branch("sumptsq2_spl_twovtx",&sumptsq2_spl_twovtx_,"sumptsq2_spl_twovtx[nrecPV_twovtx]/D");
135		pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
136		pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
137		pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
#	Line 181 | Line 142 | PVStudy::PVStudy(const edm::ParameterSet
142		pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
143		pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
144		pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
145	+
146	+	// Information for the splitVertexColl1
147	+	pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
148	+	pvtxtree_->Branch("ndofPV1_spl_twovtx",&ndofPV1_spl_twovtx_,"ndofPV1_spl_twovtx[nrecPV_twovtx]/D");
149	+	pvtxtree_->Branch("normchi2PV1_spl_twovtx",&normchi2PV1_spl_twovtx_,"normchi2PV1_spl_twovtx[nrecPV_twovtx]/D");
150	+	pvtxtree_->Branch("avgPtPV1_spl_twovtx",&avgPtPV1_spl_twovtx_,"avgPtPV1_spl_twovtx[nrecPV_twovtx]/D");
151	+	pvtxtree_->Branch("errx1_spl_twovtx",&errx1_spl_twovtx_,"errx1_spl_twovtx[nrecPV_twovtx]/D");
152	+	pvtxtree_->Branch("erry1_spl_twovtx",&erry1_spl_twovtx_,"erry1_spl_twovtx[nrecPV_twovtx]/D");
153	+	pvtxtree_->Branch("errz1_spl_twovtx",&errz1_spl_twovtx_,"errz1_spl_twovtx[nrecPV_twovtx]/D");
154	+
155	+
156	+	// Information for the splitVertexColl2
157	+	pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
158	+	pvtxtree_->Branch("ndofPV2_spl_twovtx",&ndofPV2_spl_twovtx_,"ndofPV2_spl_twovtx[nrecPV_twovtx]/D");
159	+	pvtxtree_->Branch("normchi2PV2_spl_twovtx",&normchi2PV2_spl_twovtx_,"normchi2PV2_spl_twovtx[nrecPV_twovtx]/D");
160	+	pvtxtree_->Branch("avgPtPV2_spl_twovtx",&avgPtPV2_spl_twovtx_,"avgPtPV2_spl_twovtx[nrecPV_twovtx]/D");
161	+	pvtxtree_->Branch("errx2_spl_twovtx",&errx2_spl_twovtx_,"errx2_spl_twovtx[nrecPV_twovtx]/D");
162	+	pvtxtree_->Branch("erry2_spl_twovtx",&erry2_spl_twovtx_,"erry2_spl_twovtx[nrecPV_twovtx]/D");
163	+	pvtxtree_->Branch("errz2_spl_twovtx",&errz2_spl_twovtx_,"errz2_spl_twovtx[nrecPV_twovtx]/D");
164	+
165
166		// MC variables
167		if(!realData_) {
#	Line 198 | Line 179 | PVStudy::PVStudy(const edm::ParameterSet
179		pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
180		pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
181		pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
182	<	pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
182	>	pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
183	>	pvtxtree_->Branch("errx_mct",&errx_mct_,"errx_mct[nrecPV_mct]/D");
184	>	pvtxtree_->Branch("erry_mct",&erry_mct_,"erry_mct[nrecPV_mct]/D");
185	>	pvtxtree_->Branch("errz_mct",&errz_mct_,"errz_mct[nrecPV_mct]/D");
186	>	pvtxtree_->Branch("nTrkPV_mct",&nTrkPV_mct_,"nTrkPV_mct[nrecPV_mct]/I");
187	>	pvtxtree_->Branch("ndofPV_mct",&ndofPV_mct_,"ndofPV_mct[nrecPV_mct]/D");
188	>	pvtxtree_->Branch("normchi2PV_mct",&normchi2PV_mct_,"normchi2PV_mct[nrecPV_mct]/D");
189	>	pvtxtree_->Branch("avgPtPV_mct",&avgPtPV_mct_,"avgPtPV_mct[nrecPV_mct]/D");
190	>
191		// For pvtxs with splittracks1
192		pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
193		pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
#	Line 207 | Line 196 | PVStudy::PVStudy(const edm::ParameterSet
196		pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
197		pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
198		pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
199	+	pvtxtree_->Branch("errx_spl1_mct",&errx_spl1_mct_,"errx_spl1_mct[nrecPV_spl1_mct]/D");
200	+	pvtxtree_->Branch("erry_spl1_mct",&erry_spl1_mct_,"erry_spl1_mct[nrecPV_spl1_mct]/D");
201	+	pvtxtree_->Branch("errz_spl1_mct",&errz_spl1_mct_,"errz_spl1_mct[nrecPV_spl1_mct]/D");
202	+	pvtxtree_->Branch("nTrkPV_spl1_mct",&nTrkPV_spl1_mct_,"nTrkPV_spl1_mct[nrecPV_spl1_mct]/I");
203	+	pvtxtree_->Branch("ndofPV_spl1_mct",&ndofPV_spl1_mct_,"ndofPV_spl1_mct[nrecPV_spl1_mct]/D");
204	+	pvtxtree_->Branch("normchi2PV_spl1_mct",&normchi2PV_spl1_mct_,"normchi2PV_spl1_mct[nrecPV_spl1_mct]/D");
205	+	pvtxtree_->Branch("avgPtPV_spl1_mct",&avgPtPV_spl1_mct_,"avgPtPV_spl1_mct[nrecPV_spl1_mct]/D");
206	+
207		// For pvtxs with splittracks1
208		pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
209		pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
#	Line 214 | Line 211 | PVStudy::PVStudy(const edm::ParameterSet
211		pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
212		pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
213		pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
214	<	pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
214	>	pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
215	>	pvtxtree_->Branch("errx_spl2_mct",&errx_spl2_mct_,"errx_spl2_mct[nrecPV_spl2_mct]/D");
216	>	pvtxtree_->Branch("erry_spl2_mct",&erry_spl2_mct_,"erry_spl2_mct[nrecPV_spl2_mct]/D");
217	>	pvtxtree_->Branch("errz_spl2_mct",&errz_spl2_mct_,"errz_spl2_mct[nrecPV_spl2_mct]/D");
218	>	pvtxtree_->Branch("nTrkPV_spl2_mct",&nTrkPV_spl2_mct_,"nTrkPV_spl2_mct[nrecPV_spl2_mct]/I");
219	>	pvtxtree_->Branch("ndofPV_spl2_mct",&ndofPV_spl2_mct_,"ndofPV_spl2_mct[nrecPV_spl2_mct]/D");
220	>	pvtxtree_->Branch("normchi2PV_spl2_mct",&normchi2PV_spl2_mct_,"normchi2PV_spl2_mct[nrecPV_spl2_mct]/D");
221	>	pvtxtree_->Branch("avgPtPV_spl2_mct",&avgPtPV_spl2_mct_,"avgPtPV_spl2_mct[nrecPV_spl2_mct]/D");
222	>
223		}
224		}
225	<
225	>
226		setRootStyle();
227
228	+	//========================================
229	+	// Booking histograms
230	+	//========================================
231	+
232		// Create a root file for histograms
233		theFile = new TFile(histoFileName_.c_str(), "RECREATE");
234		// make diretories
#	Line 227 | Line 236 | PVStudy::PVStudy(const edm::ParameterSet
236		theFile->cd();
237		theFile->mkdir("Others");
238		theFile->cd();
239	+	if (analyze_) {
240	+	theFile->mkdir("Results");
241	+	theFile->cd();
242	+	}
243
244		// Book Histograms:
245		h_pvtrk   = new PVHistograms();
233	–	h_pixvtx  = new PVHistograms();
246		h_misc    = new PVHistograms();
247		h_summary = new PVHistograms();
248		h_others  = new PVHistograms();
#	Line 240 | Line 252 | PVStudy::PVStudy(const edm::ParameterSet
252		else {
253		stringstream ss;
254		ss << i;
255	<	h_pvtrk->Init("pvTrk",ss.str(),"spl");
255	>	h_pvtrk->Init("pvTrk", ss.str(),"spl");
256		}
257		}
246	–	h_pixvtx->Init("pixVtx");
258		h_misc->Init("misc");
259
260		// Book MC only plots
#	Line 280 | Line 291 | PVStudy::PVStudy(const edm::ParameterSet
291		if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
292		suffix.clear();
293		} // End of Book histograms sensitive to data/mc
294	+
295
296		}
297
298		PVStudy::~PVStudy()
299		{
300	+
301	+	// do anything here that needs to be done at desctruction time
302	+	// (e.g. close files, deallocate resources etc.)
303		theFile->cd();
304		theFile->cd("Summary");
305		h_pvtrk->Save();
291	–	h_pixvtx->Save();
292	–	h_gen->Save();
306		h_misc->Save();
307		h_summary->Save();
308	<	if (analyze_)
308	>	if (!realData_)
309	>	h_gen->Save();
310	>	if (analyze_) {
311	>	theFile->cd();
312	>	theFile->cd("Results");
313		h_ana->Save();
314	+	}
315		theFile->cd();
316		theFile->cd("Others");
317		h_others->Save();
#	Line 322 | 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 335 | Line 353 | void PVStudy::setRootStyle() {
353		//
354		std::vector<PVStudy::simPrimaryVertex> PVStudy::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
355		{
356	<	std::vector<PVStudy::simPrimaryVertex> simpv;
356	>	std::vector<PVStudy::simPrimaryVertex> simpv;
357		const HepMC::GenEvent* evt=evtMC->GetEvent();
358		if (evt) {
359		edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
#	Line 512 | 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;
516	–	//=======================================================
538		// Initialize Root-tuple variables if needed
518	–	//=======================================================
519	–	//if(saventuple_)
539		SetVarToZero();
540
541	<	//=======================================================
523	<	// Track accessors
524	<	//=======================================================
525	<	//trackCollection
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 532 | 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	<
536	<	//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 543 | 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	<
547	<	//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 552 | Line 566 | PVStudy::analyze(const edm::Event& iEven
566		if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
567		splitTrackColl2 = splitTrackCollection2Handle.product();
568		} else {
569	<	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
569	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
570		}
571	<
572	<
559	<	//=======================================================
560	<	// Fill trackparameters of the input tracks to pvtx fitter
561	<	//=======================================================
562	<	edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
563	<	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype)
564	<	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
565	<	fillTrackHisto(trackColl, 0);
566	<	fillTrackHisto(splitTrackColl1, 1);
567	<	fillTrackHisto(splitTrackColl2, 2);
568	<	edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
569	<
570	<	//=======================================================
571	<	// PVTX accessors
572	<	//=======================================================
573	<	//vertexCollection
571	>
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 578 | Line 577 | PVStudy::analyze(const edm::Event& iEven
577		if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
578		vertexColl = vertexCollectionHandle.product();
579		} else {
580	<	edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
580	>	edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
581		}
582	<
584	<	//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 590 | Line 588 | PVStudy::analyze(const edm::Event& iEven
588		splitVertexColl1 = splitVertexCollection1Handle.product();
589		} else {
590		edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
591	<	}
592	<
595	<	//splitVertexCollection2
591	>	}
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 602 | Line 599 | PVStudy::analyze(const edm::Event& iEven
599		} else {
600		edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
601		}
602	+
603
604	<	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track and vertex collections"<<endl;
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	<	//=======================================================
611	<	// MC simvtx accessor
612	<	//=======================================================
610	>	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex collections"<<endl;
611	>
612	>	// ========== MC simvtx accessor
613		if (!realData_) {
614		edm::Handle<SimVertexContainer> simVtxs;
615		iEvent.getByLabel( simG4_, simVtxs);
#	Line 616 | Line 618 | PVStudy::analyze(const edm::Event& iEven
618		iEvent.getByLabel( simG4_, simTrks);
619		}
620
621	<	//=======================================================
620	<	// GET PDT
621	<	//=======================================================
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	<	// GET pixelVertices
630	<	//=======================================================
631	<	static const reco::VertexCollection s_empty_pixelVertexColl;
632	<	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
633	<	edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
634	<	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
635	<	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
636	<	pixelVertexColl = pixelVertexCollectionHandle.product();
637	<	} else {
638	<	edm::LogInfo("Debug") << "[PVStudy] pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
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		}
641
642	<	//=======================================================
643	<	// Fill pixelVertices related histograms
644	<	//=======================================================
645	<	nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
646	<	if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
647	<	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
648	<	fillTrackHistoInPV(pixelVertexColl, 4, false, true);
649	<	h_pixvtx->Fill1d("dzErr_pxlpvtx", pixelVertexColl->begin()->zError());
650	<	// Get the dZ error of the tracks in the leading pixelVertexColl
651	<	for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
652	<	t!= (pixelVertexColl->begin())->tracks_end(); t++) {
653	<
654	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
655	<	//        h_pvtrk->Fill1d("trkPtPV", 0.);
656	<	}
657	<	else
658	<	h_pixvtx->Fill1d("trkdzErr_pxlpvtx", (**t).dzError());
659	<	}
660	<	// Fill the track->dz() in the PV difference from first pixelpvtx
661	<	for(reco::Vertex::trackRef_iterator t = (vertexColl->begin())->tracks_begin();
662	<	t!= (vertexColl->begin())->tracks_end(); t++) {
663	<	// illegal charge
664	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
665	<	//        h_pvtrk->Fill1d("trkPtPV", 0.);
665	<	}
666	<	else {
667	<	if(pixelVertexColl->size()>0) {
668	<	h_pixvtx->Fill1d("trkdz_pxlpvtxdz", (**t).dz() -  pixelVertexColl->begin()->z());
669	<	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_pxlpvtxdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
670	<	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_trkdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
671	<	h_pixvtx->Fill1d("trkdzErr_pvtx", (**t).dzError());
672	<	}
673	<	}
674	<	}
642	>	//========================================================================
643	>	// Step 1:  Apply event cleaning for data and MC
644	>	//          WARNING: event selection cut are hard coded!!
645	>	//========================================================================
646	>	// =====Loop over the trackColl to get the fraction of HighPurity tracks
647	>	int nTracks = 0;
648	>	int nHighPurityTracks=0;
649	>	double fHighPurity=0;
650	>
651	>	for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
652	>	TrackRef tkref(trackCollectionHandle,i);
653	>	if(tkref->quality(reco::TrackBase::highPurity)) ++nHighPurityTracks;
654	>	}
655	>	if(nTracks>0)
656	>	fHighPurity =  double(nHighPurityTracks)/double(nTracks);
657	>	if(verbose_)
658	>	cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;
659	>
660	>	if( (fHighPurity<0.2 && nTracks>10) || vertexColl->begin()->isFake()) {
661	>	glob_runno_ = iEvent.id().run();
662	>	glob_evtno_ = iEvent.id().event();
663	>	glob_ls_   = iEvent.luminosityBlock();
664	>	glob_bx_  = iEvent.bunchCrossing();
665	>	return;
666		}
676	–
677	–	//=======================================================
678	–	// Fill number of reconstructed vertices
679	–	//=======================================================
667
668	+	//========================================================================
669	+	// Step 2: Fill histograms for the splitting consistency checks
670	+	//========================================================================
671	+
672	+	// === Fill trackparameters of the input tracks to pvtx fitter
673	+	edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
674	+	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs)
675	+	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
676	+	fillTrackHisto(trackColl, 0, beamSpot);
677	+	fillTrackHisto(splitTrackColl1, 1, beamSpot);
678	+	fillTrackHisto(splitTrackColl2, 2, beamSpot);
679	+	edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
680	+
681	+
682	+	// ==== Fill number of reconstructed vertices
683		edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
684		edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
685		edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
#	Line 687 | Line 689 | PVStudy::analyze(const edm::Event& iEven
689		printRecVtxs(splitVertexCollection1Handle);
690		printRecVtxs(splitVertexCollection2Handle);
691		}
692	+
693		nrecPV_ = int (vertexColl->size());
694		nrecPV1_spl_ = int (splitVertexColl1->size());
695		nrecPV2_spl_ = int (splitVertexColl2->size());
#	Line 696 | Line 699 | PVStudy::analyze(const edm::Event& iEven
699		h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
700		h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
701
699	–	//=================================================================
700	–	// Fill track parameter ntuple/hist for tracks used in recoVertices
701	–	//=================================================================
702	–
703	–	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
704	–	if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
705	–	fillTrackHistoInPV(vertexColl, 0, true, true);
702
703	<	if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
704	<	fillTrackHistoInPV(splitVertexColl1, 1, false, true);
705	<
706	<	if(splitVertexColl1->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
711	<	fillTrackHistoInPV(splitVertexColl2, 2, false, true);
712	<
713	<
714	<	//=======================================================
715	<	// Compare offlinePrimaryVertices with pixelVertices
716	<	//=======================================================
717	<	if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
718	<	&& (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {
719	<	h_pixvtx->Fill1d("nrecPV_minus_nrecPxlPV", double (nrecPV_ - nrecPV_pxlpvtx_));
720	<	// difference in reconstructed position of the leading pvtx
721	<	//edm::LogInfo("Debug")<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;
722	<	//edm::LogInfo("Debug")<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
723	<	h_pixvtx->Fill1d("recxPV_minus_recxPxlPV", recx_[0] - recx_pxlpvtx_[0]);
724	<	h_pixvtx->Fill1d("recyPV_minus_recyPxlPV", recy_[0] - recy_pxlpvtx_[0]);
725	<	h_pixvtx->Fill1d("reczPV_minus_reczPxlPV", recz_[0] - recz_pxlpvtx_[0]);
726	<	}
727	<
728	<
729	<
730	<	//==========================================================
731	<	// Fill secondary/primary min separations for multi-vertices
732	<	//==========================================================
703	>	// ======= Fill track parameter ntuple/hist for tracks used in recoVertices
704	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
705	>	if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
706	>	fillTrackHistoInPV(vertexColl, 0, beamSpot);
707
708	+	// ======= Fill secondary/primary min separations for multi-vertices
709		if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {
735	–
710		double min_xsep = 9999.0;
711		double min_ysep = 9999.0;
712		double min_zsep = 9999.0;
#	Line 789 | Line 763 | PVStudy::analyze(const edm::Event& iEven
763		min_zsep_ = min_zsep;
764		min_ntrksep_ = min_ntrksep;
765		min_sumpt2sep_ = min_sumpt2sep;
792	–
793	–
766		} // End of  if(vertexColl->size()>1) {
767
768		edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
797	–
798	–	//=======================================================
799	–	//           PrimaryVertex Matching
800	–	// 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
801	–	// 2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
802	–	// Assume the first match is the primary vertex,
803	–	// since vertexColl are sorted in decreasing order of sum(pT)
804	–	//=======================================================
769
770	+	//========================================================================
771	+	// Step 3:  PrimaryVertex Matching
772	+	//   1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
773	+	//   2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
774	+	//   3. The first match is the primary vertex, which has largest sum(pT^2)
775	+	//========================================================================
776	+
777		edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
778		reco::VertexCollection s_empty_matchedVertexColl1;
779		reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
#	Line 825 | Line 796 | PVStudy::analyze(const edm::Event& iEven
796		edm::LogInfo("Debug")<<"[PVStudy] Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
797		edm::LogInfo("Debug")<<"[PVStudy] recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
798		edm::LogInfo("Debug")<<"[PVStudy] recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
799	<
799	>
800	>	double twovtxsig = (recvtx1->z()-recvtx2->z())/max(recvtx1->zError(), recvtx2->zError());
801	>	h_misc->Fill1d("twovtxzsign", twovtxsig);
802		if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
803		edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
804
#	Line 853 | Line 826 | PVStudy::analyze(const edm::Event& iEven
826		}
827		}
828		edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
829	<
830	<	//=======================================================
831	<	//  Analyze matchedVertexColls
832	<	//=======================================================
829	>
830	>	//========================================================================
831	>	// Step 4: Analyze matchedVertexColls
832	>	//========================================================================
833		edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
834
835	+	// ==== If it is MC, analyze the res/pull of the unsplit vertexColl first
836	+	if(MC){
837	+	// make a list of primary vertices:
838	+	std::vector<simPrimaryVertex> simpv;
839	+	simpv=getSimPVs(evtMC,"");
840	+	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
841	+	h_gen->Fill1d("nsimPV", simpv.size());
842	+	nsimPV_ = simpv.size();
843	+	int isimpv = 0;
844	+	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
845	+	vsim!=simpv.end(); vsim++, isimpv++){
846	+	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
847	+	simx_[isimpv] = vsim->x;
848	+	simy_[isimpv] = vsim->y;
849	+	simz_[isimpv] = vsim->y;
850	+	simptsq_[isimpv] = vsim->ptsq;
851	+	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
852	+	fillMCHisto(vsim, isimpv, vertexColl, 3);
853	+	}
854	+	}
855	+
856		// Compare the reconstructed vertex position and calculate resolution/pulls
857		if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
858	<	//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
859	<	fillTrackHistoInPV(matchedVertexColl1, 1, true, false);
860	<	fillTrackHistoInPV(matchedVertexColl2, 2, true, false);
858	>	// ==== Analyze the splitted vtx using MC method
859	>	if(MC){
860	>	// make a list of primary vertices:
861	>	std::vector<simPrimaryVertex> simpv;
862	>	simpv=getSimPVs(evtMC,"");
863	>	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
864	>	int isimpv = 0;
865	>	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
866	>	vsim!=simpv.end(); vsim++, isimpv++){
867	>	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
868	>	fillMCHisto(vsim, isimpv, matchedVertexColl1, 1);
869	>	fillMCHisto(vsim, isimpv, matchedVertexColl2, 2);
870	>	}
871	>	}
872	>
873	>	// ==== Analyze res/pull two-vertex method
874	>	//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype)
875	>	fillTrackHistoInPV(matchedVertexColl1, 1, beamSpot);
876	>	fillTrackHistoInPV(matchedVertexColl2, 2, beamSpot);
877
878		reco::VertexCollection::const_iterator v1;
879		reco::VertexCollection::const_iterator v2;
880	+
881		nrecPV_twovtx_ = 0;
882		for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
883		v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
884	<	++v1, ++v2) {
884	>	++v1, ++v2) {
885	>	h_misc->Fill1d("ndofPVDiff", v1->ndof() - v2->ndof());
886	>	h_misc->Fill1d("ndofPVRelDiff", (v1->ndof()-v2->ndof())/(v1->ndof()+v2->ndof()));
887		fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
888		fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
889		fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
#	Line 895 | Line 908 | PVStudy::analyze(const edm::Event& iEven
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
899	–	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
900	–	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
901	–	sumptsq1_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v1);
902	–	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 910 | 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];
923	+
924	+
925	+	//SplittedVertex
926	+
927	+	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
928	+	ndofPV1_spl_twovtx_[nrecPV_twovtx_] = v1->ndof();
929	+	normchi2PV1_spl_twovtx_[nrecPV_twovtx_] = v1->normalizedChi2();
930	+	avgPtPV1_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v1);
931	+	errx1_spl_twovtx_[nrecPV_twovtx_] = v1->xError();
932	+	erry1_spl_twovtx_[nrecPV_twovtx_] = v1->yError();
933	+	errz1_spl_twovtx_[nrecPV_twovtx_] = v1->zError();
934	+
935	+	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
936	+	ndofPV2_spl_twovtx_[nrecPV_twovtx_] = v2->ndof();
937	+	normchi2PV2_spl_twovtx_[nrecPV_twovtx_] = v2->normalizedChi2();
938	+	avgPtPV2_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v2);
939	+	errx2_spl_twovtx_[nrecPV_twovtx_] = v2->xError();
940	+	erry2_spl_twovtx_[nrecPV_twovtx_] = v2->yError();
941	+	errz2_spl_twovtx_[nrecPV_twovtx_] = v2->zError();
942	+
943		nrecPV_twovtx_++;
944	+
945	+
946	+
947	+	// Print some information of the two tracks events
948	+	if(verbose_ && fntrk_ < 4) {
949	+	cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
950	+	printRecVtx(*v1);
951	+	cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
952	+	printRecVtx(*v2);
953	+	}
954		} // End of analyzing res/pull
955		} // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
956		else
#	Line 918 | Line 958 | PVStudy::analyze(const edm::Event& iEven
958
959		edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
960
961	<	//=======================================================
962	<	// Fill simulated vertices
963	<	//=======================================================
964	<	if (!realData_) {
925	<	bool MC=false;
926	<	Handle<HepMCProduct> evtMC;
927	<	iEvent.getByLabel("generator",evtMC);
928	<	if (!evtMC.isValid()) {
929	<	MC=false;
930	<	edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
931	<	} else {
932	<	edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
933	<	MC=true;
934	<	}
935	<
936	<	if(MC){
937	<	// make a list of primary vertices:
938	<	std::vector<simPrimaryVertex> simpv;
939	<	simpv=getSimPVs(evtMC,"");
940	<	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
941	<	h_gen->Fill1d("nsimPV", simpv.size());
942	<	nsimPV_ = simpv.size();
943	<
944	<	int isimpv = 0;
945	<	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
946	<	vsim!=simpv.end(); vsim++, isimpv++){
947	<	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
948	<	simx_[isimpv] = vsim->x;
949	<	simy_[isimpv] = vsim->y;
950	<	simz_[isimpv] = vsim->y;
951	<	simptsq_[isimpv] = vsim->ptsq;
952	<
953	<	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
954	<	fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
955	<	fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
956	<	fillMCHisto(vsim, isimpv, vertexColl, 3);
957	<	}
958	<	} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
959	<	} // End of   if (!realData_) {
960	<
961	<	// Finally fill the ftree_
961	>
962	>	//========================================================================
963	>	// Step 5: Fill ntuple if saventuple_ is on
964	>	//========================================================================
965		if(saventuple_) {
966		ftree_->Fill();
967		pvtxtree_->Fill();
968		}
969	+
970		}
971
972	+
973	+
974		// ------------ method called once each job just before starting event loop  ------------
975		void
976		PVStudy::beginJob()
#	Line 1010 | Line 1016 | PVStudy::endJob() {
1016		pvtxtree_->Write();
1017		file_->Close();
1018		}
1019	<	}
1014	<
1019	>
1020
1021	+	}
1022
1023		// Match a recovertex with a simvertex
1024		// ? Should the cut be parameter dependent?
#	Line 1020 | Line 1026 | PVStudy::endJob() {
1026		bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
1027		const reco::Vertex & vrec)
1028		{
1029	<	return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1029	>	if(vrec.isFake() || !vrec.isValid()) return false;
1030	>	else
1031	>	return true;
1032	>	//return (fabs(vsim.z-vrec.z())<0.0500); // =500um // remove this hard cut
1033		}
1034
1035		// Match two reconstructed vertices
#	Line 1087 | Line 1096 | void PVStudy::printRecVtxs(const Handle<
1096		}
1097		}
1098
1099	+	void PVStudy::printRecVtx(const reco::Vertex & v){
1100	+
1101	+	cout << "#trk " << std::setw(3) << v.tracksSize()
1102	+	<< " chi2 " << std::setw(4) << v.chi2()
1103	+	<< " ndof " << std::setw(3) << v.ndof() << endl
1104	+	<< " x "  << std::setw(8) <<std::fixed << std::setprecision(4) << v.x()
1105	+	<< " dx " << std::setw(8) << v.xError()<< endl
1106	+	<< " y "  << std::setw(8) << v.y()
1107	+	<< " dy " << std::setw(8) << v.yError()<< endl
1108	+	<< " z "  << std::setw(8) << v.z()
1109	+	<< " dz " << std::setw(8) << v.zError()
1110	+	<< endl;
1111	+	}
1112	+
1113		void PVStudy::printSimTrks(const Handle<SimTrackContainer> simTrks){
1114		cout <<  " simTrks   type, (momentum), vertIndex, genpartIndex"  << endl;
1115		int i=1;
#	Line 1127 | Line 1150 | void PVStudy::fillHisto(res r, error er,
1150		h_others->Fill1d("errPVz"+suffix, er.z());
1151		}
1152
1153	+
1154		PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1155		map<int,double> data;
1156		data.clear();
#	Line 1189 | Line 1213 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1213		ntk);
1214		}
1215
1216	+
1217		void PVStudy::fit(TH1 *hdist, double fitPars[]){
1218		TAxis *axis0 = hdist->GetXaxis();
1219		int nbin = axis0->GetLast();
1220		double nOF = hdist->GetBinContent(nbin+1);
1221		double nUF = hdist->GetBinContent(0);
1222		//   double fitRange = axis0->GetBinUpEdge(nbin);
1223	<	double fitRange = axis0->GetXmax();
1223	>	// double fitRange = axis0->GetXmax();
1224	>	double fitRange = 2*hdist->GetRMS();
1225		double sigMax[2] = {0.,0.};
1226
1227		edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
#	Line 1219 | Line 1245 | void PVStudy::fit(TH1 *hdist, double fit
1245
1246		TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1247		fgaus->SetParameter(1, 0.);
1248	<	fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
1248	>	fgaus->SetParLimits(1, -fitRange/10., fitRange/10.);
1249		fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
1250		fgaus->SetLineColor(4);
1251	<	hdist->Fit(fgaus,"Q");
1251	>	hdist->Fit(fgaus,"QLRM");
1252		gPad->Update();
1253		TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
1254		s->SetOptStat(1111111);
#	Line 1235 | Line 1261 | void PVStudy::fit(TH1 *hdist, double fit
1261		fitPars[0] = -999;
1262		}
1263
1264	<	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype) {
1264	>
1265	>	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs) {
1266		string suffix;
1267		suffix = ""; // for unsplit tracks
1268		if(datatype == 1) suffix = "1_spl"; // for splittracks 1
#	Line 1247 | Line 1274 | void PVStudy::fillTrackHisto(const reco:
1274		++itTrack) {
1275		h_pvtrk->Fill1d("trkPt"+suffix, itTrack->pt());
1276		h_pvtrk->Fill1d("trkDxy"+suffix, itTrack->dxy());
1277	+	h_pvtrk->Fill1d("trkDxyCorr"+suffix, itTrack->dxy(bs));
1278		h_pvtrk->Fill1d("trkDz"+suffix, itTrack->dz());
1279		h_pvtrk->Fill1d("trkEta"+suffix, itTrack->eta());
1280		h_pvtrk->Fill1d("trkPhi"+suffix, itTrack->phi());
1281		}
1282		}
1283
1284	<	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
1284	>	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
1285		string suffix;
1286		suffix = ""; // for unsplit tracks
1287		if(datatype == 1) suffix = "1_spl"; // for splittracks 1
#	Line 1261 | Line 1289 | void PVStudy::fillTrackHistoInPV(const r
1289		int ivtx = 0;
1290		for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1291		v!=vertexColl->end(); ++v, ++ivtx) {
1292	<	if(fillNtuple) {
1293	<	if(datatype == 0) {
1294	<	nTrkPV_[ivtx] = v->tracksSize();
1295	<	sumptsq_[ivtx] = sumPtSquared(*v);
1296	<	isValid_[ivtx] = v->isValid();
1297	<	isFake_[ivtx] = v->isFake();
1298	<	recx_[ivtx] = v->x();
1299	<	recy_[ivtx] = v->y();
1300	<	recz_[ivtx] = v->z();
1301	<	recx_err_[ivtx] = v->xError();
1302	<	recy_err_[ivtx] = v->yError();
1303	<	recz_err_[ivtx] = v->zError();
1304	<	}
1305	<	if(datatype == 1) {
1306	<	nTrkPV1_spl_[ivtx] = v->tracksSize();
1307	<	sumptsq1_spl_[ivtx] = sumPtSquared(*v);
1308	<	isValid1_spl_[ivtx] = v->isValid();
1309	<	isFake1_spl_[ivtx] = v->isFake();
1310	<	recx1_spl_[ivtx] = v->x();
1311	<	recy1_spl_[ivtx] = v->y();
1312	<	recz1_spl_[ivtx] = v->z();
1285	<	recx1_err_spl_[ivtx] = v->xError();
1286	<	recy1_err_spl_[ivtx] = v->yError();
1287	<	recz1_err_spl_[ivtx] = v->zError();
1288	<
1289	<	}
1290	<	if(datatype == 2) {
1291	<	nTrkPV2_spl_[ivtx] = v->tracksSize();
1292	<	sumptsq2_spl_[ivtx] = sumPtSquared(*v);
1293	<	isValid2_spl_[ivtx] = v->isValid();
1294	<	isFake2_spl_[ivtx] = v->isFake();
1295	<	recx2_spl_[ivtx] = v->x();
1296	<	recy2_spl_[ivtx] = v->y();
1297	<	recz2_spl_[ivtx] = v->z();
1298	<	recx2_err_spl_[ivtx] = v->xError();
1299	<	recy2_err_spl_[ivtx] = v->yError();
1300	<	recz2_err_spl_[ivtx] = v->zError();
1301	<	}
1302	<	if(datatype == 4) { // for pixelVertices
1303	<	nTrkPV_pxlpvtx_[ivtx] = v->tracksSize();
1304	<	sumptsq_pxlpvtx_[ivtx] = sumPtSquared(*v);
1305	<	isValid_pxlpvtx_[ivtx] = v->isValid();
1306	<	isFake_pxlpvtx_[ivtx] = v->isFake();
1307	<	recx_pxlpvtx_[ivtx] = v->x();
1308	<	recy_pxlpvtx_[ivtx] = v->y();
1309	<	recz_pxlpvtx_[ivtx] = v->z();
1310	<	recx_err_pxlpvtx_[ivtx] = v->xError();
1311	<	recy_err_pxlpvtx_[ivtx] = v->yError();
1312	<	recz_err_pxlpvtx_[ivtx] = v->zError();
1313	<	}
1314	<	}
1315	<	}
1316	<	// For histogram only fill the leading pvtx parameters
1317	<	if(fillHisto) {
1318	<	h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1319	<	try {
1320	<	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1321	<	t!=vertexColl->begin()->tracks_end(); t++) {
1322	<	// illegal charge
1323	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1324	<	//      h_pvtrk->Fill1d("trkPtPV", 0.);
1325	<	}
1326	<	else {
1327	<	h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1328	<	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1329	<	h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1330	<	h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1331	<	h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1292	>	if(!v->isFake()) {
1293	>	h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1294	>	h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1295	>	int nHWTrkPV_ = 0;
1296	>	try {
1297	>	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1298	>	t!=vertexColl->begin()->tracks_end(); t++) {
1299	>	// illegal charge
1300	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1301	>	//    h_pvtrk->Fill1d("trkPtPV", 0.);
1302	>	}
1303	>	else {
1304	>	h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1305	>	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1306	>	h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
1307	>	h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1308	>	h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1309	>	h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1310	>	if(vertexColl->begin()->trackWeight(*t) > 0.5)
1311	>	nHWTrkPV_++;
1312	>	}
1313		}
1314		}
1315	<	}
1316	<	catch (...) {
1317	<	// exception thrown when trying to use linked track
1318	<	//          h_pvtrk->Fill1d("trkPtPV", 0.);
1315	>	catch (...) {
1316	>	// exception thrown when trying to use linked track
1317	>	//        h_pvtrk->Fill1d("trkPtPV", 0.);
1318	>	}
1319	>	h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
1320		}
1321		}
1340	–
1322		}
1323
1324		void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
#	Line 1363 | Line 1344 | void PVStudy::fillMCHisto(std::vector<si
1344
1345		for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1346		vrec!=vtxColl->end(); ++vrec){
1366	–	int nrectrk = vrec->tracksSize();
1347		vsim->recVtx=NULL;
1368	–
1348		edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1349		edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1350
1351	<	if ( matchVertex(*vsim,*vrec) && vrec->isValid() && !vrec->isFake() ) {
1351	>	if ( matchVertex(*vsim,*vrec)  ) {
1352		vsim->recVtx=&(*vrec);
1353
1354		// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
#	Line 1381 | Line 1360 | void PVStudy::fillMCHisto(std::vector<si
1360
1361		double fres_mct[3];
1362		double ferror_mct[3];
1363	<
1363	>	int nrectrk =  int(vsim->recVtx->tracksSize());
1364	>	double ndofPV = vsim->recVtx->ndof();
1365	>	double normchi2PV =  vsim->recVtx->normalizedChi2();
1366	>	double avgPtPV = avgPtRecVtx(*(vsim->recVtx));
1367	>
1368		fres_mct[0] = vsim->recVtx->x()-vsim->x;
1369		fres_mct[1] = vsim->recVtx->y()-vsim->y;
1370		fres_mct[2] = vsim->recVtx->z()-vsim->z;
#	Line 1407 | Line 1390 | void PVStudy::fillMCHisto(std::vector<si
1390		nrectrk, datatype);
1391
1392		if(saventuple_) {
1393	<	//Fill the values for variables in ftree_
1393	>	//Fill the values for variables in pvtxtree_
1394		if(datatype == 1) {
1395		nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1396	+	ndofPV_spl1_mct_[nrecPV_spl1_mct_] =  ndofPV;
1397	+	normchi2PV_spl1_mct_[nrecPV_spl1_mct_] =  normchi2PV;
1398	+	avgPtPV_spl1_mct_[nrecPV_spl1_mct_] = avgPtPV;
1399		deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1400	<	deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1401	<	deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1400	>	deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[1];
1401	>	deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2];
1402		pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
1403		pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1404		pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1405	+	errx_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[0];
1406	+	erry_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[1];
1407	+	errz_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[2];
1408		nrecPV_spl1_mct_++;
1409		}
1421	–
1410		if(datatype == 2) {
1411	<	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1411	>	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1412	>	ndofPV_spl2_mct_[nrecPV_spl2_mct_] =  ndofPV;
1413	>	normchi2PV_spl2_mct_[nrecPV_spl2_mct_] =  normchi2PV;
1414	>	avgPtPV_spl2_mct_[nrecPV_spl2_mct_] = avgPtPV;
1415		deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1416	<	deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1417	<	deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1416	>	deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[1];
1417	>	deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2];
1418		pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
1419		pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1420	<	pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1420	>	pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1421	>	errx_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[0];
1422	>	erry_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[1];
1423	>	errz_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[2];
1424		nrecPV_spl2_mct_++;
1425		}
1426		if(datatype == 3) {
1427	<	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1427	>	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1428	>	ndofPV_mct_[nrecPV_mct_] =  ndofPV;
1429	>	normchi2PV_mct_[nrecPV_mct_] =  normchi2PV;
1430	>	avgPtPV_mct_[nrecPV_mct_] = avgPtPV;
1431		deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1432		deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1433		deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1434		pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
1435		pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1436	<	pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1436	>	pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1437	>	errx_mct_[nrecPV_mct_] =  ferror_mct[0];
1438	>	erry_mct_[nrecPV_mct_] =  ferror_mct[1];
1439	>	errz_mct_[nrecPV_mct_] =  ferror_mct[2];
1440		nrecPV_mct_++;
1441		}
1442		} // End of  if(saventuple_) {
#	Line 1447 | Line 1447 | void PVStudy::fillMCHisto(std::vector<si
1447		}
1448		}
1449
1450	–
1450		void PVStudy::SetVarToZero() {
1452	–	//Greg's variables
1451		fntrk_ = 0;
1452		//pvtx position (x,y,z) residual and error
1453		for(int i = 0; i<3;i++)    {
#	Line 1457 | Line 1455 | void PVStudy::SetVarToZero() {
1455		ferror_[i] = 0;
1456		}
1457
1458	<	//Yanyan's variables
1459	<	// Number of reconstructed vertices
1458	>	// Event level information
1459	>	glob_runno_ = 0;
1460	>	glob_evtno_ = 0;
1461	>	glob_ls_ = 0;
1462	>	glob_bx_ = 0;
1463		nrecPV_ = 0;
1464		nrecPV1_spl_ = 0;
1465		nrecPV2_spl_ = 0;
#	Line 1474 | Line 1475 | void PVStudy::SetVarToZero() {
1475		min_ntrksep_ = 9999.0;
1476		min_sumpt2sep_ = -9999.0;
1477
1478	<
1478	>	// Variables filled per Vertex
1479		for (int i = 0; i < nMaxPVs_; i++) {
1479	–	// recoVertices with all tracks
1480	–	nTrkPV_[i] = 0; // Number of tracks in the pvtx
1481	–	sumptsq_[i] = 0;
1482	–	isValid_[i] = -1;
1483	–	isFake_[i] = -1;
1484	–	recx_[i] = 0;
1485	–	recy_[i] = 0;
1486	–	recz_[i] = 0;
1487	–	recx_err_[i] = 0;
1488	–	recy_err_[i] = 0;
1489	–	recz_err_[i] = 0;
1490	–
1491	–	// recoVertices with splitTrack1
1492	–	nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx
1493	–	sumptsq1_spl_[i] = 0;
1494	–	isValid1_spl_[i] = -1;
1495	–	isFake1_spl_[i] = -1;
1496	–	recx1_spl_[i] = 0;
1497	–	recy1_spl_[i] = 0;
1498	–	recz1_spl_[i] = 0;
1499	–	recx1_err_spl_[i] = 0;
1500	–	recy1_err_spl_[i] = 0;
1501	–	recz1_err_spl_[i] = 0;
1502	–
1503	–	// recoVertices with splitTrack2
1504	–	nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx
1505	–	sumptsq2_spl_[i] = 0;
1506	–	isValid2_spl_[i] = -1;
1507	–	isFake2_spl_[i] = -1;
1508	–	recx2_spl_[i] = 0;
1509	–	recy2_spl_[i] = 0;
1510	–	recz2_spl_[i] = 0;
1511	–	recx2_err_spl_[i] = 0;
1512	–	recy2_err_spl_[i] = 0;
1513	–	recz2_err_spl_[i] = 0;
1514	–
1515	–	//pixelVertices
1516	–	nTrkPV_pxlpvtx_[i] = 0; // Number of tracks in the pvtx
1517	–	sumptsq_pxlpvtx_[i] = 0;
1518	–	isValid_pxlpvtx_[i] = -1;
1519	–	isFake_pxlpvtx_[i] = -1;
1520	–	recx_pxlpvtx_[i] = 0;
1521	–	recy_pxlpvtx_[i] = 0;
1522	–	recz_pxlpvtx_[i] = 0;
1523	–	recx_err_pxlpvtx_[i] = 0;
1524	–	recy_err_pxlpvtx_[i] = 0;
1525	–	recz_err_pxlpvtx_[i] = 0;
1526	–
1480		// matched two-vertices
1481		nTrkPV1_spl_twovtx_[i] = 0;
1482	+	ndofPV1_spl_twovtx_[i] = 0;
1483	+	normchi2PV1_spl_twovtx_[i] = 0;
1484	+	avgPtPV1_spl_twovtx_[i] = 0;
1485	+	errx1_spl_twovtx_[i] = 0;
1486	+	erry1_spl_twovtx_[i] = 0;
1487	+	errz1_spl_twovtx_[i] = 0;
1488	+
1489		nTrkPV2_spl_twovtx_[i] = 0;
1490	<	nTrkPV_twovtx_[i] = 0;
1491	<	sumptsq1_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1492	<	sumptsq2_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1490	>	ndofPV2_spl_twovtx_[i] = 0;
1491	>	normchi2PV2_spl_twovtx_[i] = 0;
1492	>	avgPtPV2_spl_twovtx_[i] = 0;
1493	>	errx2_spl_twovtx_[i] = 0;
1494	>	erry2_spl_twovtx_[i] = 0;
1495	>	errz2_spl_twovtx_[i] = 0;
1496	>
1497		deltax_twovtx_[i] = 0;
1498		deltay_twovtx_[i] = 0;
1499		deltaz_twovtx_[i] = 0;
#	Line 1554 | Line 1518 | void PVStudy::SetVarToZero() {
1518		pullx_mct_[i] = 0;
1519		pully_mct_[i] = 0;
1520		pullz_mct_[i] = 0;
1521	<
1521	>	nTrkPV_mct_[i] = 0;
1522	>	ndofPV_mct_[i] = 0;
1523	>	normchi2PV_mct_[i] = 0;
1524	>	avgPtPV_mct_[i] = 0;
1525	>	errx_mct_[i] = 0;
1526	>	erry_mct_[i] = 0;
1527	>	errz_mct_[i] = 0;
1528	>
1529		deltax_spl1_mct_[i] = 0;
1530		deltay_spl1_mct_[i] = 0;
1531		deltaz_spl1_mct_[i] = 0;
1532		pullx_spl1_mct_[i] = 0;
1533		pully_spl1_mct_[i] = 0;
1534		pullz_spl1_mct_[i] = 0;
1535	<
1535	>	nTrkPV_spl1_mct_[i] = 0;
1536	>	ndofPV_spl1_mct_[i] = 0;
1537	>	normchi2PV_spl1_mct_[i] = 0;
1538	>	avgPtPV_spl1_mct_[i] = 0;
1539	>	errx_spl1_mct_[i] = 0;
1540	>	erry_spl1_mct_[i] = 0;
1541	>	errz_spl1_mct_[i] = 0;
1542	>
1543		deltax_spl2_mct_[i] = 0;
1544		deltay_spl2_mct_[i] = 0;
1545		deltaz_spl2_mct_[i] = 0;
1546		pullx_spl2_mct_[i] = 0;
1547		pully_spl2_mct_[i] = 0;
1548		pullz_spl2_mct_[i] = 0;
1549	+	nTrkPV_spl2_mct_[i] = 0;
1550	+	ndofPV_spl2_mct_[i] = 0;
1551	+	normchi2PV_spl2_mct_[i] = 0;
1552	+	avgPtPV_spl2_mct_[i] = 0;
1553	+	errx_spl2_mct_[i] = 0;
1554	+	erry_spl2_mct_[i] = 0;
1555	+	errz_spl2_mct_[i] = 0;
1556	+
1557		}
1572	–
1558		}
1559
1560		double PVStudy::sumPtSquared(const reco::Vertex & v)  {
1561	<	double sum = 0.;
1562	<	double pT;
1563	<	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1564	<	pT = (**it).pt();
1565	<	sum += pT*pT;
1561	>	if(v.isFake() || (!v.isValid()) || v.tracksSize() == 0) return 0.0;
1562	>	else {
1563	>	double sum = 0.;
1564	>	double pT;
1565	>	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1566	>	pT = (**it).pt();
1567	>	sum += pT*pT;
1568	>	}
1569	>	return sum;
1570		}
1582	–	return sum;
1571		}
1572
1573	+
1574	+
1575	+	double PVStudy::avgPtRecVtx(const reco::Vertex & v)  {
1576	+
1577	+	if(v.isFake() || !v.isValid() || v.tracksSize()==0 ) return 0;
1578	+	else {
1579	+	double sumpT = 0.;
1580	+	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1581	+	sumpT +=  (**it).pt();
1582	+	}
1583	+	return sumpT/double(v.tracksSize());
1584	+	}
1585	+	}
1586	+
1587	+
1588	+

Diff Legend

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