[ViewVC] Diff of: cvsroot/UserCode/Jeng/PVStudy/plugins/PVStudy.cc

Comparing UserCode/Jeng/PVStudy/plugins/PVStudy.cc (file contents):
Revision 1.12 by yygao, Fri Nov 27 12:31:05 2009 UTC vs.
Revision 1.21 by yygao, Sun Feb 28 00:07:09 2010 UTC

#	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		{
69		//=======================================================================
70		// Get configuration for TrackTupleMaker
71		//=======================================================================
72	<	simG4_ = iConfig.getParameter<edm::InputTag>( "simG4" );
73	<	trackCollectionTag_ = iConfig.getUntrackedParameter<edm::InputTag>("trackCollection");
74	<	splitTrackCollection1Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection1");
75	<	splitTrackCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
72	>	simG4_ = iConfig.getParameter<edm::InputTag>( "simG4" );
73	>	trackCollectionTag_ = iConfig.getUntrackedParameter<edm::InputTag>("trackCollection");
74	>	splitTrackCollection1Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection1");
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");
79	<	pixelVertexCollectionTag_ = iConfig.getParameter<edm::InputTag>("pixelVertexCollectionTag");
80	<	verbose_ = iConfig.getUntrackedParameter<bool>("verbose",false);
81	<	realData_ = iConfig.getUntrackedParameter<bool>("realData",false);
82	<	analyze_ = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
83	<	saventuple_ = iConfig.getUntrackedParameter<bool>("saventuple",false);
84	<	outputfilename_ = iConfig.getUntrackedParameter<string>("OutputFileName");
85	<	ntrkdiffcut_ = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
86	<	zsigncut_ = iConfig.getUntrackedParameter<int>("zsigncut");
87	<	nTrkMin_ = iConfig.getUntrackedParameter<int>("nTrkMin");
88	<	nTrkMax_ = iConfig.getUntrackedParameter<int>("nTrkMax");
89	<
77	>	splitVertexCollection1Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
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");
84	>	histoFileName_ = iConfig.getUntrackedParameter<std::string> ("histoFileName");
85	>	ntrkdiffcut_ = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
86	>	nTrkMin_ = iConfig.getUntrackedParameter<int>("nTrkMin");
87	>	nTrkMax_ = iConfig.getUntrackedParameter<int>("nTrkMax");
88	>	zsigncut_ = iConfig.getUntrackedParameter<int>("zsigncut");
89	>	fHighPuritycut_ = iConfig.getUntrackedParameter<double>("fHighPuritycut");
90	>	useHWTrk_ = iConfig.getUntrackedParameter<bool>("useHWTrk",false);
91	>	bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
92	>
93		//now do what ever initialization is needed
94		pvinfo.clear();
95
90	–
96		// Specify the data mode vector
97		if(realData_) datamode.push_back(0);
98		else {
#	Line 96 \| Line 101 \| PVStudy::PVStudy(const edm::ParameterSet
101		datamode.push_back(2);
102		datamode.push_back(3);
103		}
104	<
100	<	// Create ntuple files if needed
104	>	// Create ntuple files if needed
105		if(saventuple_) {
106		file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
107		ftree_ = new TTree("mytree","mytree");
#	Line 108 \| Line 112 \| PVStudy::PVStudy(const edm::ParameterSet
112
113		// pvtxtree_ analyzing the pvtxs ootb
114		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");
115
116	+	// Event information for the data
117	+	pvtxtree_->Branch("glob_runno",&glob_runno_,"glob_runno/I");
118	+	pvtxtree_->Branch("glob_evtno",&glob_evtno_,"glob_evtno/I");
119	+	pvtxtree_->Branch("glob_ls",&glob_ls_,"glob_ls/I");
120	+	pvtxtree_->Branch("glob_bx",&glob_bx_,"glob_bx/I");
121	+
122	+
123	+	pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
124	+	pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
125	+	pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
126	+	// Event level information
127		pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
128		pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
129		pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
130
131	<	//pvtx using splittracks1
132	<	pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
133	<	pvtxtree_->Branch("nTrkPV1_spl",&nTrkPV1_spl_,"nTrkPV1_spl[nrecPV1_spl]/I");
134	<	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");
131	>	//Fill the variables in the twovtx pair (recvtx1, recvtx2)
132	>	// Information related to the analyzing the two-vertex method
133	>
134	>
135
167	–	//Fill the variables in the twovtx pair (recvtx1, recvtx2)
136		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");
137		pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
138		pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
139		pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
#	Line 180 \| Line 144 \| PVStudy::PVStudy(const edm::ParameterSet
144		pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
145		pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
146		pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
147	+
148	+	// Information for the splitVertexColl1
149	+	pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
150	+	pvtxtree_->Branch("ndofPV1_spl_twovtx",&ndofPV1_spl_twovtx_,"ndofPV1_spl_twovtx[nrecPV_twovtx]/D");
151	+	pvtxtree_->Branch("normchi2PV1_spl_twovtx",&normchi2PV1_spl_twovtx_,"normchi2PV1_spl_twovtx[nrecPV_twovtx]/D");
152	+	pvtxtree_->Branch("avgPtPV1_spl_twovtx",&avgPtPV1_spl_twovtx_,"avgPtPV1_spl_twovtx[nrecPV_twovtx]/D");
153	+	pvtxtree_->Branch("errx1_spl_twovtx",&errx1_spl_twovtx_,"errx1_spl_twovtx[nrecPV_twovtx]/D");
154	+	pvtxtree_->Branch("erry1_spl_twovtx",&erry1_spl_twovtx_,"erry1_spl_twovtx[nrecPV_twovtx]/D");
155	+	pvtxtree_->Branch("errz1_spl_twovtx",&errz1_spl_twovtx_,"errz1_spl_twovtx[nrecPV_twovtx]/D");
156	+
157	+
158	+	// Information for the splitVertexColl2
159	+	pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
160	+	pvtxtree_->Branch("ndofPV2_spl_twovtx",&ndofPV2_spl_twovtx_,"ndofPV2_spl_twovtx[nrecPV_twovtx]/D");
161	+	pvtxtree_->Branch("normchi2PV2_spl_twovtx",&normchi2PV2_spl_twovtx_,"normchi2PV2_spl_twovtx[nrecPV_twovtx]/D");
162	+	pvtxtree_->Branch("avgPtPV2_spl_twovtx",&avgPtPV2_spl_twovtx_,"avgPtPV2_spl_twovtx[nrecPV_twovtx]/D");
163	+	pvtxtree_->Branch("errx2_spl_twovtx",&errx2_spl_twovtx_,"errx2_spl_twovtx[nrecPV_twovtx]/D");
164	+	pvtxtree_->Branch("erry2_spl_twovtx",&erry2_spl_twovtx_,"erry2_spl_twovtx[nrecPV_twovtx]/D");
165	+	pvtxtree_->Branch("errz2_spl_twovtx",&errz2_spl_twovtx_,"errz2_spl_twovtx[nrecPV_twovtx]/D");
166	+
167
168		// MC variables
169		if(!realData_) {
#	Line 197 \| Line 181 \| PVStudy::PVStudy(const edm::ParameterSet
181		pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
182		pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
183		pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
184	<	pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
184	>	pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
185	>	pvtxtree_->Branch("errx_mct",&errx_mct_,"errx_mct[nrecPV_mct]/D");
186	>	pvtxtree_->Branch("erry_mct",&erry_mct_,"erry_mct[nrecPV_mct]/D");
187	>	pvtxtree_->Branch("errz_mct",&errz_mct_,"errz_mct[nrecPV_mct]/D");
188	>	pvtxtree_->Branch("nTrkPV_mct",&nTrkPV_mct_,"nTrkPV_mct[nrecPV_mct]/I");
189	>	pvtxtree_->Branch("ndofPV_mct",&ndofPV_mct_,"ndofPV_mct[nrecPV_mct]/D");
190	>	pvtxtree_->Branch("normchi2PV_mct",&normchi2PV_mct_,"normchi2PV_mct[nrecPV_mct]/D");
191	>	pvtxtree_->Branch("avgPtPV_mct",&avgPtPV_mct_,"avgPtPV_mct[nrecPV_mct]/D");
192	>
193		// For pvtxs with splittracks1
194		pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
195		pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
#	Line 206 \| Line 198 \| PVStudy::PVStudy(const edm::ParameterSet
198		pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
199		pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
200		pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
201	+	pvtxtree_->Branch("errx_spl1_mct",&errx_spl1_mct_,"errx_spl1_mct[nrecPV_spl1_mct]/D");
202	+	pvtxtree_->Branch("erry_spl1_mct",&erry_spl1_mct_,"erry_spl1_mct[nrecPV_spl1_mct]/D");
203	+	pvtxtree_->Branch("errz_spl1_mct",&errz_spl1_mct_,"errz_spl1_mct[nrecPV_spl1_mct]/D");
204	+	pvtxtree_->Branch("nTrkPV_spl1_mct",&nTrkPV_spl1_mct_,"nTrkPV_spl1_mct[nrecPV_spl1_mct]/I");
205	+	pvtxtree_->Branch("ndofPV_spl1_mct",&ndofPV_spl1_mct_,"ndofPV_spl1_mct[nrecPV_spl1_mct]/D");
206	+	pvtxtree_->Branch("normchi2PV_spl1_mct",&normchi2PV_spl1_mct_,"normchi2PV_spl1_mct[nrecPV_spl1_mct]/D");
207	+	pvtxtree_->Branch("avgPtPV_spl1_mct",&avgPtPV_spl1_mct_,"avgPtPV_spl1_mct[nrecPV_spl1_mct]/D");
208	+
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");
#	Line 213 \| Line 213 \| PVStudy::PVStudy(const edm::ParameterSet
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");
216	>	pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
217	>	pvtxtree_->Branch("errx_spl2_mct",&errx_spl2_mct_,"errx_spl2_mct[nrecPV_spl2_mct]/D");
218	>	pvtxtree_->Branch("erry_spl2_mct",&erry_spl2_mct_,"erry_spl2_mct[nrecPV_spl2_mct]/D");
219	>	pvtxtree_->Branch("errz_spl2_mct",&errz_spl2_mct_,"errz_spl2_mct[nrecPV_spl2_mct]/D");
220	>	pvtxtree_->Branch("nTrkPV_spl2_mct",&nTrkPV_spl2_mct_,"nTrkPV_spl2_mct[nrecPV_spl2_mct]/I");
221	>	pvtxtree_->Branch("ndofPV_spl2_mct",&ndofPV_spl2_mct_,"ndofPV_spl2_mct[nrecPV_spl2_mct]/D");
222	>	pvtxtree_->Branch("normchi2PV_spl2_mct",&normchi2PV_spl2_mct_,"normchi2PV_spl2_mct[nrecPV_spl2_mct]/D");
223	>	pvtxtree_->Branch("avgPtPV_spl2_mct",&avgPtPV_spl2_mct_,"avgPtPV_spl2_mct[nrecPV_spl2_mct]/D");
224	>
225		}
226		}
219	–
220	–	edm::Service<TFileService> fs;
221	–	TFileDirectory subDir = fs->mkdir( "Summary" );
222	–	TFileDirectory subDir1 = fs->mkdir( "Others" );
223	–	// TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
224	–
225	–	setRootStyle();
227
228	<	//Book histograms for track and pvtx parameters for each splitted pvtx
228	<	//i=0 : x (as in unsplit track collection)
229	<	//i=1 : x1_spl_
230	<	//i=2 : x2_spl_
228	>	setRootStyle();
229
230	<
231	<	for(int i=0;i<3;i++)
232	<	{
233	<	string suffix;
234	<	if(i == 0) suffix = "";
235	<	else {
236	<	stringstream ss;
237	<	ss << i;
238	<	suffix =ss.str()+"_spl";
239	<	}
240	<	h["nTrk"+suffix] = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
241	<	h["trkPt"+suffix] = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
242	<	h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#eta of rec tracks "+suffix),100,-3,3);
243	<	h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
244	<	h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-0.5,0.5);
245	<	h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
246	<
247	<	h["nTrkPV"+suffix] = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
248	<	h["trkPtPV"+suffix] = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
249	<	h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#eta of rec tracks in PV"+suffix),100,-3,3);
250	<	h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
251	<	h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);
252	<	h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
253	<	h["nrecPV"+suffix] = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);
254	<	suffix.clear();
230	>	//========================================
231	>	// Booking histograms
232	>	//========================================
233	>
234	>	// Create a root file for histograms
235	>	theFile = new TFile(histoFileName_.c_str(), "RECREATE");
236	>	// make diretories
237	>	theFile->mkdir("Summary");
238	>	theFile->cd();
239	>	theFile->mkdir("Others");
240	>	theFile->cd();
241	>	if (analyze_) {
242	>	theFile->mkdir("Results");
243	>	theFile->cd();
244	>	}
245	>
246	>	// Book Histograms:
247	>	h_pvtrk = new PVHistograms();
248	>	h_misc = new PVHistograms();
249	>	h_summary = new PVHistograms();
250	>	h_others = new PVHistograms();
251	>
252	>	for(int i=0;i<3;i++) {
253	>	if(i == 0) h_pvtrk->Init("pvTrk");
254	>	else {
255	>	stringstream ss;
256	>	ss << i;
257	>	h_pvtrk->Init("pvTrk", ss.str(),"spl");
258		}
259	<	h["nrecPVDiff"] = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
260	<	h["nTrkPVDiff"] = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
261	<	h["nTrkPVRelDiff"] = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
262	<
263	<	// Histograms on comparing the multi-vertices
264	<	// Difference in reconstructed vtx position
265	<	h["min_xsep"] = subDir.make<TH1D>("min_xsep", "min x diff of primary and secondary pvtx",300,0,0.1);
266	<	h["min_xsep_sign"] = subDir.make<TH1D>("min_xsep_sign", "min x diff in signf of primary and secondary pvtx",300,0,5);
267	<	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);
259	>	}
260	>	h_misc->Init("misc");
261	>
262	>	// Book MC only plots
263	>	if (!realData_) {
264	>	h_gen = new PVHistograms();
265	>	h_gen->Init("generator");
266	>	}
267	>	if (analyze_) h_ana = new PVHistograms();
268
269		//Book histograms sensitive to data/mc
270		for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
271		string suffix;
272	<	if(verbose_) cout<<"datamode = "<< *it<<endl;
272	>	edm::LogInfo("Debug")<<"datamode = "<< *it<<endl;
273		switch(*it) {
274		case 0: suffix = "";
275		break;
#	Line 285 \| Line 280 \| PVStudy::PVStudy(const edm::ParameterSet
280		case 3: suffix = "_mct";
281		break;
282		}
283	<	h["deltax"+suffix] = subDir.make<TH1D>(TString("deltax"+suffix), TString("x-residual pvtx"+suffix),800,-0.04,0.04);
289	<	h["deltay"+suffix] = subDir.make<TH1D>(TString("deltay"+suffix), TString("y-residual pvtx"+suffix),800,-0.04,0.04);
290	<	h["deltaz"+suffix] = subDir.make<TH1D>(TString("deltaz"+suffix), TString("z-residual pvtx"+suffix),800,-0.04,0.04);
291	<	h["pullx"+suffix] = subDir.make<TH1D>(TString("pullx"+suffix), TString("x-pull pvtx"+suffix),800,-10,10);
292	<	h["pully"+suffix] = subDir.make<TH1D>(TString("pully"+suffix), TString("y-pull pvtx"+suffix),800,-10,10);
293	<	h["pullz"+suffix] = subDir.make<TH1D>(TString("pullz"+suffix), TString("z-pull pvtx"+suffix),800,-10,10);
294	<	h["errPVx"+suffix] = subDir.make<TH1D>(TString("errPVx"+suffix), TString("X"+suffix+" vertex error"),100,0.,0.02);
295	<	h["errPVy"+suffix] = subDir.make<TH1D>(TString("errPVy"+suffix), TString("Y"+suffix+" vertex error"),100,0.,0.02);
296	<	h["errPVz"+suffix] = subDir.make<TH1D>(TString("errPVz"+suffix), TString("Z"+suffix+" vertex error"),100,0.,0.02);
283	>	h_summary->Init("summary",suffix);
284
285		// Book residual, error and pull histograms for each nTrk bin
286		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
287		stringstream ss;
288		ss << ntrk;
289	<	string suffix2 = suffix+"_"+ss.str();
290	<	h["deltax"+suffix2] = subDir1.make<TH1D>(TString("deltax"+suffix2), TString("x-residual of pvtx"+suffix),100,-0.02,0.02);
291	<	h["deltax"+suffix2]->GetXaxis()->SetTitle("cm");
305	<	h["deltay"+suffix2] = subDir1.make<TH1D>(TString("deltay"+suffix2), TString("y-residual of pvtx"+suffix),100,-0.02,0.02);
306	<	h["deltay"+suffix2]->GetXaxis()->SetTitle("cm");
307	<	h["deltaz"+suffix2] = subDir1.make<TH1D>(TString("deltaz"+suffix2), TString("z-residual of pvtx"+suffix),100,-0.02,0.02);
308	<	h["deltaz"+suffix2]->GetXaxis()->SetTitle("cm");
309	<	h["pullx"+suffix2] = subDir1.make<TH1D>(TString("pullx"+suffix2), TString("x-pull of pvtx"+suffix),100,-10.,10.);
310	<	h["pully"+suffix2] = subDir1.make<TH1D>(TString("pully"+suffix2), TString("y-pull of pvtx"+suffix),100,-10.,10.);
311	<	h["pullz"+suffix2] = subDir1.make<TH1D>(TString("pullz"+suffix2), TString("z-pull of pvtx"+suffix),100,-10.,10.);
312	<	h["errPVx"+suffix2] = subDir1.make<TH1D>(TString("errPVx"+suffix2), TString("X"+suffix+" vertex error"),100,0.,0.02);
313	<	h["errPVy"+suffix2] = subDir1.make<TH1D>(TString("errPVy"+suffix2), TString("Y"+suffix+" vertex error"),100,0.,0.02);
314	<	h["errPVz"+suffix2] = subDir1.make<TH1D>(TString("errPVz"+suffix2), TString("Z"+suffix+" vertex error"),100,0.,0.02);
315	<	suffix2.clear();
316	<	} // End of for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
317	<
289	>	h_others->Init("others",suffix,ss.str());
290	>	}
291	>
292		// Book residual and pull histograms only when analyzeOntheFly is enabled
293	<	if(analyze_) {
320	<	h2["resx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resx"+suffix+"_nTrk"), TString("x-resolution vs number of tracks in pvtx"+suffix),nTrkMax_-nTrkMin_,nTrkMin_-0.5,nTrkMax_-0.5,400,0.,200);
321	<	h2["resx"+suffix+"_nTrk"]->SetMarkerStyle(21);
322	<	h2["resx"+suffix+"_nTrk"]->SetMarkerColor(4);
323	<	h2["resx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
324	<	h2["resx"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
325	<	h2["resy"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resy"+suffix+"_nTrk"), TString("y-resolution vs number of tracks in pvtx"+suffix),nTrkMax_-nTrkMin_,nTrkMin_-0.5,nTrkMax_-0.5,400,0.,200);
326	<	h2["resy"+suffix+"_nTrk"]->SetMarkerStyle(21);
327	<	h2["resy"+suffix+"_nTrk"]->SetMarkerColor(4);
328	<	h2["resy"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
329	<	h2["resy"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
330	<	h2["resz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resz"+suffix+"_nTrk"), TString("z-resolution vs number of tracks in pvtx"+suffix),nTrkMax_-nTrkMin_,nTrkMin_-0.5,nTrkMax_-0.5,400,0.,200);
331	<	h2["resz"+suffix+"_nTrk"]->SetMarkerStyle(21);
332	<	h2["resz"+suffix+"_nTrk"]->SetMarkerColor(4);
333	<	h2["resz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
334	<	h2["resz"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
335	<	h2["pullx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullx"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),nTrkMax_-nTrkMin_,nTrkMin_-0.5,nTrkMax_-0.5,100,0.,2.);
336	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerStyle(21);
337	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerColor(4);
338	<	h2["pullx"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
339	<	h2["pullx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
340	<	h2["pully"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pully"+suffix+"_nTrk"), TString("y-pull vs number of tracks"+suffix),nTrkMax_-nTrkMin_,nTrkMin_-0.5,nTrkMax_-0.5,100,0.,2.);
341	<	h2["pully"+suffix+"_nTrk"]->SetMarkerStyle(21);
342	<	h2["pully"+suffix+"_nTrk"]->SetMarkerColor(4);
343	<	h2["pully"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
344	<	h2["pully"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
345	<	h2["pullz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullz"+suffix+"_nTrk"), TString("z-pull vs number of tracks"+suffix),nTrkMax_-nTrkMin_,nTrkMin_-0.5,nTrkMax_-0.5,100,0.,2.);
346	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerStyle(21);
347	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerColor(4);
348	<	h2["pullz"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
349	<	h2["pullz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
350	<	} // End of if(analyze_) {
293	>	if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
294		suffix.clear();
295		} // End of Book histograms sensitive to data/mc
296
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);
297
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);
370	–	h["genPart_T"]->GetXaxis()->SetTitle("t (nanosecond)");
371	–	h["genPart_cT"] = subDir.make<TH1D>("genPart_cT","ct component of gen particles",300,-150.,150.);
372	–	h["genPart_cT"]->GetXaxis()->SetTitle("ct (mm)");
373	–	h["nsimPV"] = subDir.make<TH1D>("nsimPV","Num of sim PV",51,-0.5,50.5);
374	–	}
298		}
299
300		PVStudy::~PVStudy()
#	Line 379 \| Line 302 \| PVStudy::~PVStudy()
302
303		// do anything here that needs to be done at desctruction time
304		// (e.g. close files, deallocate resources etc.)
305	+	theFile->cd();
306	+	theFile->cd("Summary");
307	+	h_pvtrk->Save();
308	+	h_misc->Save();
309	+	h_summary->Save();
310	+	if (!realData_)
311	+	h_gen->Save();
312	+	if (analyze_) {
313	+	theFile->cd();
314	+	theFile->cd("Results");
315	+	h_ana->Save();
316	+	}
317	+	theFile->cd();
318	+	theFile->cd("Others");
319	+	h_others->Save();
320
321	+	theFile->Close();
322		}
323
324		void PVStudy::setRootStyle() {
#	Line 403 \| Line 342 \| void PVStudy::setRootStyle() {
342		gStyle->SetTitleSize(0.055, ""); // size for pad title; default is 0.02
343		gStyle->SetLabelSize(0.03, "XYZ"); // size for axis labels; default is 0.04
344		gStyle->SetStatFontSize(0.08); // size for stat. box
345	<	gStyle->SetTitleFont(32, "XYZ"); // times-bold-italic font (p. 153) for axes
346	<	gStyle->SetTitleFont(32, ""); // same for pad title
347	<	gStyle->SetLabelFont(32, "XYZ"); // same for axis labels
348	<	gStyle->SetStatFont(32); // same for stat. box
345	>	gStyle->SetTitleFont(42, "XYZ"); // times-bold-italic font (p. 153) for axes
346	>	gStyle->SetTitleFont(42, ""); // same for pad title
347	>	gStyle->SetLabelFont(42, "XYZ"); // same for axis labels
348	>	gStyle->SetStatFont(42); // same for stat. box
349		gStyle->SetLabelOffset(0.006, "Y"); // default is 0.005
350		//
351		return;
#	Line 416 \| Line 355 \| void PVStudy::setRootStyle() {
355		//
356		std::vector<PVStudy::simPrimaryVertex> PVStudy::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
357		{
358	<	std::vector<PVStudy::simPrimaryVertex> simpv;
358	>	std::vector<PVStudy::simPrimaryVertex> simpv;
359		const HepMC::GenEvent* evt=evtMC->GetEvent();
360		if (evt) {
361	<	if(verbose_) {
362	<	cout << "process id " << evt->signal_process_id()<<endl;
363	<	cout << "signal process vertex " << ( evt->signal_process_vertex() ?
364	<	evt->signal_process_vertex()->barcode() : 0 ) <<endl;
426	<	cout << "number of vertices " << evt->vertices_size() << endl;
427	<	}
361	>	edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
362	>	edm::LogInfo("SimPVs") << "[getSimPVs] signal process vertex " << ( evt->signal_process_vertex() ?
363	>	evt->signal_process_vertex()->barcode() : 0 ) <<endl;
364	>	edm::LogInfo("SimPVs") << "[getSimPVs] number of vertices " << evt->vertices_size() << endl;
365
366		int idx=0; int npv=0;
367		for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
#	Line 435 \| Line 372 \| std::vector<PVStudy::simPrimaryVertex> P
372		// t component of PV:
373		for ( HepMC::GenVertex::particle_iterator mother = (*vitr)->particles_begin(HepMC::parents);
374		mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
375	<	// std::cout << "Status = " << (*mother)->status() << std::endl;
375	>	// edm::LogInfo("SimPVs") << "Status = " << (*mother)->status() << endl;
376		HepMC::GenVertex * mv=(*mother)->production_vertex();
377		if( ((*mother)->status() == 3) && (!mv)) {
378	<	// std::cout << "npv= " << npv << std::endl;
378	>	// edm::LogInfo("SimPVs") << "npv= " << npv << endl;
379		if (npv == 0) {
380	<	h["genPart_cT"]->Fill(pos.t()); // mm
381	<	h["genPart_T"]->Fill(pos.t()/299.792458); // ns
380	>	h_gen->Fill1d("genPart_cT", pos.t()); // mm
381	>	h_gen->Fill1d("genPart_T", pos.t()/299.792458); // ns
382		}
383		npv++;
384		}
#	Line 449 \| Line 386 \| std::vector<PVStudy::simPrimaryVertex> P
386		// if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
387
388		bool hasMotherVertex=false;
389	<	if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl;
389	>	if (verbose_) cout << "[getSimPVs] mothers of vertex[" << ++idx << "]: " << endl;
390		for ( HepMC::GenVertex::particle_iterator mother = (*vitr)->particles_begin(HepMC::parents);
391		mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
392		HepMC::GenVertex * mv=(*mother)->production_vertex();
393	<	// std::cout << "Status = " << (*mother)->status() << std::endl;
393	>	// if (verbose_) cout << "Status = " << (*mother)->status() << endl;
394		if (mv) {
395		hasMotherVertex=true;
396		if(!verbose_) break; //if verbose_, print all particles of gen vertices
#	Line 480 \| Line 417 \| std::vector<PVStudy::simPrimaryVertex> P
417
418		if(!vp){
419		// this is a new vertex
420	<	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
420	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
421		simpv.push_back(sv);
422		vp=&simpv.back();
423		}else{
424	<	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
424	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
425		}
426		vp->genVertex.push_back((*vitr)->barcode());
427		// collect final state descendants
#	Line 508 \| Line 445 \| std::vector<PVStudy::simPrimaryVertex> P
445		}
446		}
447		}
448	<	}
449	<	}
448	>	}//loop MC vertices daughters
449	>	}//loop MC vertices
450		}
451		return simpv;
452		}
#	Line 542 \| Line 479 \| std::vector<PVStudy::simPrimaryVertex> P
479		ptsq+=(daughter)->momentum().perp()(*daughter)->momentum().perp();
480		}
481		}
482	<	primary = ( gv->position().t()==0);
482	>	//primary = ( gv->position().t()==0 );
483	>	primary = true;
484	>	h_gen->Fill1d("genPart_cT", gv->position().t()); // mm
485	>	h_gen->Fill1d("genPart_T", gv->position().t()/299.792458); // ns
486	>
487		//resonance= ( gp->mother() && isResonance(gp->mother())); // in CLHEP/HepMC days
488		// no more mother pointer in the improved HepMC GenParticle
489		resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
#	Line 591 \| Line 532 \| PVStudy::analyze(const edm::Event& iEven
532		using namespace edm;
533		using namespace reco;
534
535	<	if (verbose_)
536	<	cout<<"PVStudy::analyze():"<<endl;
537	<	//=======================================================
535	>	//========================================================================
536	>	// Step 0: Prepare root variables and get information from the Event
537	>	//========================================================================
538	>
539	>	edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
540		// Initialize Root-tuple variables if needed
598	–	//=======================================================
599	–	//if(saventuple_)
541		SetVarToZero();
542
543	<	//=======================================================
603	<	// Track accessors
604	<	//=======================================================
605	<	//trackCollection
543	>	// ====== TrackCollection
544		static const reco::TrackCollection s_empty_trackColl;
545		const reco::TrackCollection *trackColl = &s_empty_trackColl;
546		edm::Handle<reco::TrackCollection> trackCollectionHandle;
#	Line 610 \| Line 548 \| PVStudy::analyze(const edm::Event& iEven
548		if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
549		trackColl = trackCollectionHandle.product();
550		} else {
551	<	cout << "trackCollection cannot be found -> using empty collection of same type." <<endl;
551	>	edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
552		}
553	<
616	<	//splitTrackCollection1
553	>	// ====== splitTrackCollection1
554		static const reco::TrackCollection s_empty_splitTrackColl1;
555		const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
556		edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
#	Line 621 \| Line 558 \| PVStudy::analyze(const edm::Event& iEven
558		if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
559		splitTrackColl1 = splitTrackCollection1Handle.product();
560		} else {
561	<	cout << "splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
561	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
562		}
563	<
627	<	//splitTrackCollection2
563	>	// ====== splitTrackCollection2
564		static const reco::TrackCollection s_empty_splitTrackColl2;
565		const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
566		edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
#	Line 632 \| Line 568 \| PVStudy::analyze(const edm::Event& iEven
568		if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
569		splitTrackColl2 = splitTrackCollection2Handle.product();
570		} else {
571	<	cout << "splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
571	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
572		}
573	<
574	<
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	<	//=======================================================
655	<	//vertexCollection
573	>
574	>	// ======= PrimaryVertexCollection
575		static const reco::VertexCollection s_empty_vertexColl;
576		const reco::VertexCollection *vertexColl = &s_empty_vertexColl;
577		edm::Handle<reco::VertexCollection> vertexCollectionHandle;
#	Line 660 \| Line 579 \| PVStudy::analyze(const edm::Event& iEven
579		if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
580		vertexColl = vertexCollectionHandle.product();
581		} else {
582	<	cout << "vertexCollection cannot be found -> using empty collection of same type." <<endl;
582	>	edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
583		}
584	<
666	<	//splitVertexCollection1
584	>	// ====== splitVertexCollection1
585		static const reco::VertexCollection s_empty_splitVertexColl1;
586		const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
587		edm::Handle<reco::VertexCollection> splitVertexCollection1Handle;
#	Line 671 \| Line 589 \| PVStudy::analyze(const edm::Event& iEven
589		if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
590		splitVertexColl1 = splitVertexCollection1Handle.product();
591		} else {
592	<	cout << "splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
593	<	}
594	<
677	<	//splitVertexCollection2
592	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
593	>	}
594	>	// ====== splitVertexCollection2
595		static const reco::VertexCollection s_empty_splitVertexColl2;
596		const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
597		edm::Handle<reco::VertexCollection> splitVertexCollection2Handle;
#	Line 682 \| Line 599 \| PVStudy::analyze(const edm::Event& iEven
599		if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
600		splitVertexColl2 = splitVertexCollection2Handle.product();
601		} else {
602	<	cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
602	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
603		}
604	+
605	+
606	+	// ======== BeamSpot accessors
607	+	edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
608	+	iEvent.getByLabel(bsSrc,recoBeamSpotHandle);
609	+	reco::BeamSpot bs = *recoBeamSpotHandle;
610	+	const Point beamSpot = recoBeamSpotHandle.isValid() ? Point(recoBeamSpotHandle->x0(), recoBeamSpotHandle->y0(), recoBeamSpotHandle->z0()) : Point(0, 0, 0);
611
612	<	if(verbose_) cout<<"End accessing the track and vertex collections"<<endl;
612	>	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex collections"<<endl;
613
614	<	//=======================================================
691	<	// MC simvtx accessor
692	<	//=======================================================
614	>	// ========== MC simvtx accessor
615		if (!realData_) {
616		edm::Handle<SimVertexContainer> simVtxs;
617		iEvent.getByLabel( simG4_, simVtxs);
#	Line 698 \| Line 620 \| PVStudy::analyze(const edm::Event& iEven
620		iEvent.getByLabel( simG4_, simTrks);
621		}
622
623	<	//=======================================================
702	<	// GET PDT
703	<	//=======================================================
623	>	// ========== GET PDT
624		try{
625		iSetup.getData(pdt);
626		}catch(const Exception&){
627	<	cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
627	>	edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
628		}
629
630	<	//=======================================================
631	<	// GET pixelVertices
632	<	//=======================================================
633	<	static const reco::VertexCollection s_empty_pixelVertexColl;
634	<	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
635	<	edm::Handle<reco::VertexCollection> pixelVertexCollectionHandle;
636	<	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
637	<	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
638	<	pixelVertexColl = pixelVertexCollectionHandle.product();
639	<	} else {
640	<	cout << "pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
721	<	}
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	<	}
630	>	// ======= Get MC information if needed
631	>	bool MC=false;
632	>	Handle<HepMCProduct> evtMC;
633	>	if (!realData_) {
634	>	iEvent.getByLabel("generator",evtMC);
635	>	if (!evtMC.isValid()) {
636	>	MC=false;
637	>	edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
638	>	} else {
639	>	edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
640	>	MC=true;
641		}
642		}
758	–
759	–	//=======================================================
760	–	// Fill number of reconstructed vertices
761	–	//=======================================================
643
644	<	if(verbose_) {
645	<	cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
644	>	//========================================================================
645	>	// Step 1: Apply event cleaning for data and MC
646	>	// WARNING: event selection cut are hard coded!!
647	>	//========================================================================
648	>	// =====Loop over the trackColl to get the fraction of HighPurity tracks
649	>	int nTracks = 0;
650	>	int nHighPurityTracks=0;
651	>	double fHighPurity=0;
652	>
653	>	for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
654	>	TrackRef tkref(trackCollectionHandle,i);
655	>	if(tkref->quality(reco::TrackBase::highPurity)) ++nHighPurityTracks;
656	>	}
657	>	if(nTracks>0)
658	>	fHighPurity = double(nHighPurityTracks)/double(nTracks);
659	>	if(verbose_)
660	>	cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;
661	>
662	>	if( (fHighPurity<fHighPuritycut_ && nTracks>10) \|\| vertexColl->begin()->isFake()) return;
663	>
664	>	glob_runno_ = iEvent.id().run();
665	>	glob_evtno_ = iEvent.id().event();
666	>	glob_ls_ = iEvent.luminosityBlock();
667	>	glob_bx_ = iEvent.bunchCrossing();
668	>
669	>	//========================================================================
670	>	// Step 2: Fill histograms for the splitting consistency checks
671	>	//========================================================================
672	>
673	>	// === Fill trackparameters of the input tracks to pvtx fitter
674	>	edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
675	>	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs)
676	>	// datatype: unsplittracks (0); splittracks1 (1); splittracks2 (2);
677	>	fillTrackHisto(trackColl, 0, beamSpot);
678	>	fillTrackHisto(splitTrackColl1, 1, beamSpot);
679	>	fillTrackHisto(splitTrackColl2, 2, beamSpot);
680	>	edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
681	>
682	>
683	>	// ==== Fill number of reconstructed vertices
684	>	edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
685	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
686	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
687	>	edm::LogInfo("Debug")<<"[PVStudy] Start filling pvtx parameters."<<endl;
688	>	if (verbose_) {
689		printRecVtxs(vertexCollectionHandle);
766	–	cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
690		printRecVtxs(splitVertexCollection1Handle);
768	–	cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
691		printRecVtxs(splitVertexCollection2Handle);
770	–	cout<<"Start filling pvtx parameters."<<endl;
692		}
693	<
693	>
694		nrecPV_ = int (vertexColl->size());
695		nrecPV1_spl_ = int (splitVertexColl1->size());
696		nrecPV2_spl_ = int (splitVertexColl2->size());
697
698	<	h["nrecPV"]->Fill(nrecPV_);
699	<	h["nrecPV1_spl"]->Fill(nrecPV1_spl_);
700	<	h["nrecPV2_spl"]->Fill(nrecPV2_spl_);
701	<	h["nrecPVDiff"]->Fill(double(nrecPV1_spl_)-double(nrecPV2_spl_));
781	<
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);
698	>	h_pvtrk->Fill1d("nrecPV", nrecPV_);
699	>	h_pvtrk->Fill1d("nrecPV1_spl", nrecPV1_spl_);
700	>	h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
701	>	h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
702
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	–	}
703
704	<
705	<
706	<	//==========================================================
707	<	// Fill secondary/primary min separations for multi-vertices
815	<	//==========================================================
704	>	// ======= Fill track parameter ntuple/hist for tracks used in recoVertices
705	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
706	>	if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
707	>	fillTrackHistoInPV(vertexColl, 0, beamSpot);
708
709	+	// ======= Fill secondary/primary min separations for multi-vertices
710		if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {
818	–
711		double min_xsep = 9999.0;
712		double min_ysep = 9999.0;
713		double min_zsep = 9999.0;
#	Line 825 \| Line 717 \| PVStudy::analyze(const edm::Event& iEven
717		double min_ntrksep = 9999.0;
718		double min_sumpt2sep = 9999999.0;
719
720	<	if(verbose_) cout<<"leading pvtx z = "<< vertexColl->begin()->z() <<endl;
720	>	edm::LogInfo("Debug")<<"[PVStudy] leading pvtx z = "<< vertexColl->begin()->z() <<endl;
721
722		// Looping through the secondary vertices to find the mininum separation to the primary
723		for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
#	Line 860 \| Line 752 \| PVStudy::analyze(const edm::Event& iEven
752		min_sumpt2sep = fabs(sumPtSquared(v) - sumPtSquared((vertexColl->begin())));
753		}
754		}
755	<	h["min_xsep"]->Fill(min_xsep);
756	<	h["min_ysep"]->Fill(min_ysep);
757	<	h["min_zsep"]->Fill(min_zsep);
758	<	h["min_xsep_sign"]->Fill(min_xsep_sign);
759	<	h["min_ysep_sign"]->Fill(min_ysep_sign);
760	<	h["min_zsep_sign"]->Fill(min_zsep_sign);
761	<	h["min_ntrksep"]->Fill(min_ntrksep);
762	<	h["min_sumpt2sep"]->Fill(min_sumpt2sep);
755	>	h_misc->Fill1d("min_xsep", min_xsep);
756	>	h_misc->Fill1d("min_ysep", min_ysep);
757	>	h_misc->Fill1d("min_zsep", min_zsep);
758	>	h_misc->Fill1d("min_xsep_sign", min_xsep_sign);
759	>	h_misc->Fill1d("min_ysep_sign", min_ysep_sign);
760	>	h_misc->Fill1d("min_zsep_sign", min_zsep_sign);
761	>	h_misc->Fill1d("min_ntrksep", min_ntrksep);
762	>	h_misc->Fill1d("min_sumpt2sep", min_sumpt2sep);
763
764		min_zsep_ = min_zsep;
765		min_ntrksep_ = min_ntrksep;
766		min_sumpt2sep_ = min_sumpt2sep;
875	–
876	–
767		} // End of if(vertexColl->size()>1) {
768
769	<	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	<	//=======================================================
769	>	edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
770
771	<	if(verbose_) cout<<"PVStudy::analyze(): matching pvtxs "<<endl;
771	>	//========================================================================
772	>	// Step 3: PrimaryVertex Matching
773	>	// 1. In z \|z1-z2\|< zsigncut * max(sigmaz1, sigmaz2)
774	>	// 2. \|(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)\|<ntrkdiffcut_
775	>	// 3. The first match is the primary vertex, which has largest sum(pT^2)
776	>	//========================================================================
777	>
778	>	edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
779		reco::VertexCollection s_empty_matchedVertexColl1;
780		reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
781		matchedVertexColl1->reserve(splitVertexColl1->size());
#	Line 906 \| Line 794 \| PVStudy::analyze(const edm::Event& iEven
794		if (ivtx!=0 \|\| jvtx!=0) { stopmatching = true; break; }
795		reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
796		if( !(recvtx2->isValid()) \|\| recvtx2->isFake()) break;
797	<	if(verbose_) {
798	<	cout<<"Matching splitVertexColl1: # "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
799	<	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
800	<	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
801	<	}
797	>	edm::LogInfo("Debug")<<"[PVStudy] Matching splitVertexColl1: # "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
798	>	edm::LogInfo("Debug")<<"[PVStudy] recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
799	>	edm::LogInfo("Debug")<<"[PVStudy] recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
800	>
801	>	double twovtxsig = (recvtx1->z()-recvtx2->z())/max(recvtx1->zError(), recvtx2->zError());
802	>	h_misc->Fill1d("twovtxzsign", twovtxsig);
803		if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
804	<	if(verbose_) {
805	<	cout<<"The two splitted vertices match in Z. "<<endl;
804	>	edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
805	>
806	>	int nTrkPV1, nTrkPV2;
807	>	if(useHWTrk_) {
808	>	nTrkPV1 = nHWTrkRecVtx(*recvtx1);
809	>	nTrkPV2 = nHWTrkRecVtx(*recvtx2);
810	>	}
811	>	else {
812	>	nTrkPV1 = recvtx1->tracksSize();
813	>	nTrkPV2 = recvtx2->tracksSize();
814		}
918	–	int nTrkPV1 = recvtx1->tracksSize();
919	–	int nTrkPV2 = recvtx2->tracksSize();
815		double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
816	<	h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
817	<	h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
816	>	h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
817	>	h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
818		if(fabs(ntrkreldiff)<ntrkdiffcut_) {
819	<	if(verbose_) {
820	<	cout<<"(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
926	<	}
819	>	edm::LogInfo("Debug")<<"[PVStudy] (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
820	>
821		matchedVertexColl1->push_back(*recvtx1);
822		matchedVertexColl2->push_back(*recvtx2);
823		stopmatching = true; // stop the matching when the first match is found!
824		break;
825		}
826		else
827	<	cout<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
827	>	edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
828		}
829		else {
830	<	cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
831	<	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
832	<	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
830	>	edm::LogInfo("Debug")<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
831	>	edm::LogInfo("Debug")<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
832	>	edm::LogInfo("Debug")<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
833		}
834		}
835		}
836	<	if(verbose_) cout<<"PVStudy::analyze(): End matching pvtxs"<<endl;
837	<
838	<	//=======================================================
839	<	// Analyze matchedVertexColls
840	<	//=======================================================
841	<	if(verbose_) {
842	<	cout<<"Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
843	<	}
836	>	edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
837	>
838	>	//========================================================================
839	>	// Step 4: Analyze matchedVertexColls
840	>	//========================================================================
841	>	edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
842	>
843	>	// ==== If it is MC, analyze the res/pull of the unsplit vertexColl first
844	>	if(MC){
845	>	// make a list of primary vertices:
846	>	std::vector<simPrimaryVertex> simpv;
847	>	simpv=getSimPVs(evtMC,"");
848	>	// simpv=getSimPVs(evtMC, simVtxs, simTrks);
849	>	h_gen->Fill1d("nsimPV", simpv.size());
850	>	nsimPV_ = simpv.size();
851	>	int isimpv = 0;
852	>	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
853	>	vsim!=simpv.end(); vsim++, isimpv++){
854	>	nsimTrkPV_[isimpv] =vsim->nGenTrk;
855	>	simx_[isimpv] = vsim->x;
856	>	simy_[isimpv] = vsim->y;
857	>	simz_[isimpv] = vsim->y;
858	>	simptsq_[isimpv] = vsim->ptsq;
859	>	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
860	>	fillMCHisto(vsim, isimpv, vertexColl, 3);
861	>	}
862	>	}
863
864		// Compare the reconstructed vertex position and calculate resolution/pulls
865		if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
866	<	//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
867	<	fillTrackHistoInPV(matchedVertexColl1, 1, true, false);
868	<	fillTrackHistoInPV(matchedVertexColl2, 2, true, false);
866	>	// ==== Analyze the splitted vtx using MC method
867	>	if(MC){
868	>	// make a list of primary vertices:
869	>	std::vector<simPrimaryVertex> simpv;
870	>	simpv=getSimPVs(evtMC,"");
871	>	// simpv=getSimPVs(evtMC, simVtxs, simTrks);
872	>	int isimpv = 0;
873	>	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
874	>	vsim!=simpv.end(); vsim++, isimpv++){
875	>	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
876	>	fillMCHisto(vsim, isimpv, matchedVertexColl1, 1);
877	>	fillMCHisto(vsim, isimpv, matchedVertexColl2, 2);
878	>	}
879	>	}
880	>
881	>	// ==== Analyze res/pull two-vertex method
882	>	//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype)
883	>	fillTrackHistoInPV(matchedVertexColl1, 1, beamSpot);
884	>	fillTrackHistoInPV(matchedVertexColl2, 2, beamSpot);
885
886		reco::VertexCollection::const_iterator v1;
887		reco::VertexCollection::const_iterator v2;
888	+
889		nrecPV_twovtx_ = 0;
890		for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
891		v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
892	<	++v1, ++v2) {
892	>	++v1, ++v2) {
893	>	h_misc->Fill1d("ndofPVDiff", v1->ndof() - v2->ndof());
894	>	h_misc->Fill1d("ndofPVRelDiff", (v1->ndof()-v2->ndof())/(v1->ndof()+v2->ndof()));
895		fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
896		fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
897		fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
898		ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
899		ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
900		ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
969	–	fntrk_ = (v1->tracksSize()+v2->tracksSize())/2;
901
902	<	h["deltax"]->Fill( fres_[0] );
903	<	h["deltay"]->Fill( fres_[1] );
904	<	h["deltaz"]->Fill( fres_[2] );
905	<	h["pullx"]->Fill( fres_[0]/ferror_[0]);
906	<	h["pully"]->Fill( fres_[1]/ferror_[1]);
907	<	h["pullz"]->Fill( fres_[2]/ferror_[2]);
908	<	h["errPVx"]->Fill( ferror_[0] );
909	<	h["errPVy"]->Fill( ferror_[1] );
910	<	h["errPVz"]->Fill( ferror_[2] );
902	>	int nTrkPV1, nTrkPV2;
903	>	if(useHWTrk_) {
904	>	nTrkPV1 = nHWTrkRecVtx(*v1);
905	>	nTrkPV2 = nHWTrkRecVtx(*v2);
906	>	}
907	>	else {
908	>	nTrkPV1 = v1->tracksSize();
909	>	nTrkPV2 = v2->tracksSize();
910	>	}
911	>
912	>	fntrk_ = (nTrkPV1 + nTrkPV2)/2;
913	>
914	>	h_summary->Fill1d("deltax", fres_[0] );
915	>	h_summary->Fill1d("deltay", fres_[1] );
916	>	h_summary->Fill1d("deltaz", fres_[2] );
917	>	h_summary->Fill1d("pullx", fres_[0]/ferror_[0]);
918	>	h_summary->Fill1d("pully", fres_[1]/ferror_[1]);
919	>	h_summary->Fill1d("pullz", fres_[2]/ferror_[2]);
920	>	h_summary->Fill1d("errPVx", ferror_[0] );
921	>	h_summary->Fill1d("errPVy", ferror_[1] );
922	>	h_summary->Fill1d("errPVz", ferror_[2] );
923		pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
924		error(ferror_[0], ferror_[1], ferror_[2]),
925		fntrk_));
#	Line 984 \| Line 927 \| PVStudy::analyze(const edm::Event& iEven
927		fillHisto(res(fres_[0], fres_[1], fres_[2]),
928		error(ferror_[0], ferror_[1], ferror_[2]),
929		fntrk_,0);
930	+
931		// 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);
932		nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
933		deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
934		deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
#	Line 999 \| Line 939 \| PVStudy::analyze(const edm::Event& iEven
939		pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
940		pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
941		pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];
942	+
943	+
944	+	//SplittedVertex
945	+	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = nTrkPV1;
946	+	ndofPV1_spl_twovtx_[nrecPV_twovtx_] = v1->ndof();
947	+	normchi2PV1_spl_twovtx_[nrecPV_twovtx_] = v1->normalizedChi2();
948	+	avgPtPV1_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v1);
949	+	errx1_spl_twovtx_[nrecPV_twovtx_] = v1->xError();
950	+	erry1_spl_twovtx_[nrecPV_twovtx_] = v1->yError();
951	+	errz1_spl_twovtx_[nrecPV_twovtx_] = v1->zError();
952	+
953	+	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = nTrkPV2;
954	+	ndofPV2_spl_twovtx_[nrecPV_twovtx_] = v2->ndof();
955	+	normchi2PV2_spl_twovtx_[nrecPV_twovtx_] = v2->normalizedChi2();
956	+	avgPtPV2_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v2);
957	+	errx2_spl_twovtx_[nrecPV_twovtx_] = v2->xError();
958	+	erry2_spl_twovtx_[nrecPV_twovtx_] = v2->yError();
959	+	errz2_spl_twovtx_[nrecPV_twovtx_] = v2->zError();
960	+
961		nrecPV_twovtx_++;
962	+
963	+	// Print some information of the two tracks events
964	+	if(verbose_ && fntrk_ < 4) {
965	+	cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
966	+	printRecVtx(*v1);
967	+	cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
968	+	printRecVtx(*v2);
969	+	}
970		} // End of analyzing res/pull
971		} // End of if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
972		else
973	<	cout<<"WARNING: Cannot find matching pvtxs"<<endl;
973	>	edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
974
975	<	if(verbose_) cout<<"End analyzing the matchedVertexColl"<<endl;
1009	<
1010	<	//=======================================================
1011	<	// Fill simulated vertices
1012	<	//=======================================================
1013	<	if (!realData_) {
1014	<	bool MC=false;
1015	<	Handle<HepMCProduct> evtMC;
1016	<	iEvent.getByLabel("generator",evtMC);
1017	<	if (!evtMC.isValid()) {
1018	<	MC=false;
1019	<	if(verbose_){
1020	<	cout << "no HepMCProduct found"<< endl;
1021	<	}
1022	<	} else {
1023	<	if(verbose_){
1024	<	cout << "generator HepMCProduct found"<< endl;
1025	<	}
1026	<	MC=true;
1027	<	}
1028	<
1029	<	if(MC){
1030	<	// make a list of primary vertices:
1031	<	std::vector<simPrimaryVertex> simpv;
1032	<	simpv=getSimPVs(evtMC,"");
1033	<	// simpv=getSimPVs(evtMC, simVtxs, simTrks);
1034	<	h["nsimPV"]->Fill(simpv.size());
1035	<	nsimPV_ = simpv.size();
975	>	edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
976
977	<	int isimpv = 0;
978	<	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
979	<	vsim!=simpv.end(); vsim++, isimpv++){
980	<	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 if (!realData_) {
1053	<
1054	<	// Finally fill the ftree_
977	>
978	>	//========================================================================
979	>	// Step 5: Fill ntuple if saventuple_ is on
980	>	//========================================================================
981		if(saventuple_) {
982		ftree_->Fill();
983		pvtxtree_->Fill();
984		}
985	+
986		}
987
988	+
989	+
990		// ------------ method called once each job just before starting event loop ------------
991		void
992		PVStudy::beginJob()
#	Line 1067 \| Line 996 \| PVStudy::beginJob()
996		// ------------ method called once each job just after ending the event loop ------------
997		void
998		PVStudy::endJob() {
999	<	if (verbose_) cout << "Analyzing PV info" << endl;
999	>	edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl;
1000		// Looping through the datamodes available
1001		for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
1002		string suffix;
1003	<	if(verbose_) cout<<"datamode = "<< *it<<endl;
1003	>	edm::LogInfo("Analysis")<<"[endJob] datamode = "<< *it<<endl;
1004		switch(*it) {
1005		case 0: suffix = "";
1006		break;
#	Line 1082 \| Line 1011 \| PVStudy::endJob() {
1011		case 3: suffix = "_mct";
1012		break;
1013		}
1014	+	suffix += "_nTrk";
1015		if(analyze_) {
1016		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1017	+	edm::LogInfo("Analysis") << "[endJob] ntrk = " << ntrk << endl;
1018		PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
1019	<	if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
1020	<	if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
1021	<	if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
1022	<	if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
1023	<	if ( ipv.pull_.y() > 0 ) h2["pully"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
1024	<	if ( ipv.pull_.z() > 0 ) h2["pullz"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.z());//, 1.);
1019	>	if ( ipv.res_.x() > 0 ) h_ana->Fill2d("resx"+suffix, ntrk,ipv.res_.x());
1020	>	if ( ipv.res_.y() > 0 ) h_ana->Fill2d("resy"+suffix, ntrk,ipv.res_.y());
1021	>	if ( ipv.res_.z() > 0 ) h_ana->Fill2d("resz"+suffix, ntrk,ipv.res_.z());
1022	>	if ( ipv.pull_.x() > 0 ) h_ana->Fill2d("pullx"+suffix, ntrk,ipv.pull_.x());
1023	>	if ( ipv.pull_.y() > 0 ) h_ana->Fill2d("pully"+suffix, ntrk,ipv.pull_.y());
1024	>	if ( ipv.pull_.z() > 0 ) h_ana->Fill2d("pullz"+suffix, ntrk,ipv.pull_.z());
1025		}
1026	<	}
1026	>	}
1027		suffix.clear();
1028		}
1029		if(saventuple_) {
#	Line 1101 \| Line 1032 \| PVStudy::endJob() {
1032		pvtxtree_->Write();
1033		file_->Close();
1034		}
1035	<	}
1105	<
1035	>
1036
1037	+	}
1038
1039		// Match a recovertex with a simvertex
1040		// ? Should the cut be parameter dependent?
#	Line 1111 \| Line 1042 \| PVStudy::endJob() {
1042		bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
1043		const reco::Vertex & vrec)
1044		{
1045	<	return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1045	>	if(vrec.isFake() \|\| !vrec.isValid()) return false;
1046	>	else
1047	>	return true;
1048	>	//return (fabs(vsim.z-vrec.z())<0.0500); // =500um // remove this hard cut
1049		}
1050
1051		// Match two reconstructed vertices
#	Line 1120 \| Line 1054 \| bool PVStudy::matchVertex( const reco::V
1054		int zsigncut)
1055		{
1056		double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
1057	<	if(verbose_)
1124	<	cout<<"The matching criteria in z is "<<cut<<endl;
1057	>	edm::LogInfo("Debug")<<"[matchVertex] The matching criteria in z is "<<cut<<endl;
1058		return (fabs(vrec1->z()-vrec2->z())<cut);
1059		}
1060
1061
1062		bool PVStudy::isResonance(const HepMC::GenParticle * p){
1063		double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
1064	<	if(verbose_) cout << "isResonance " << p->pdg_id() << " " << ctau << endl;
1064	>	edm::LogInfo("Debug") << "[isResonance] isResonance " << p->pdg_id() << " " << ctau << endl;
1065		return ctau >0 && ctau <1e-6;
1066		}
1067
#	Line 1179 \| Line 1112 \| void PVStudy::printRecVtxs(const Handle<
1112		}
1113		}
1114
1115	+	void PVStudy::printRecVtx(const reco::Vertex & v){
1116	+
1117	+	cout << "#trk " << std::setw(3) << v.tracksSize()
1118	+	<< " chi2 " << std::setw(4) << v.chi2()
1119	+	<< " ndof " << std::setw(3) << v.ndof() << endl
1120	+	<< " x " << std::setw(8) <<std::fixed << std::setprecision(4) << v.x()
1121	+	<< " dx " << std::setw(8) << v.xError()<< endl
1122	+	<< " y " << std::setw(8) << v.y()
1123	+	<< " dy " << std::setw(8) << v.yError()<< endl
1124	+	<< " z " << std::setw(8) << v.z()
1125	+	<< " dz " << std::setw(8) << v.zError()
1126	+	<< endl;
1127	+	}
1128	+
1129		void PVStudy::printSimTrks(const Handle<SimTrackContainer> simTrks){
1130		cout << " simTrks type, (momentum), vertIndex, genpartIndex" << endl;
1131		int i=1;
#	Line 1200 \| Line 1147 \| void PVStudy::printSimTrks(const Handle<
1147		void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1148		stringstream ss;
1149		ss << ntk;
1150	<	string suffix = ss.str();
1151	<	if(datamode == 0 ) suffix = "_" + suffix;
1152	<	if(datamode == 1 ) suffix = "_spl1_mct_" + suffix;
1153	<	if(datamode == 2 ) suffix = "_spl2_mct_" + suffix;
1154	<	if(datamode == 3 ) suffix = "_mct_" + suffix;
1150	>	string suffix;
1151	>	if(datamode == 0 ) suffix = "_" + ss.str();
1152	>	if(datamode == 1 ) suffix = "_spl1_mct_" + ss.str();
1153	>	if(datamode == 2 ) suffix = "_spl2_mct_" + ss.str();
1154	>	if(datamode == 3 ) suffix = "_mct_" + ss.str();
1155
1156		if (ntk < nTrkMin_ \|\| ntk > nTrkMax_ ) return;
1157		// Fill histograms of res and pull of ntk
1158	<	h["deltax"+suffix]->Fill(r.x());
1159	<	h["deltay"+suffix]->Fill(r.y());
1160	<	h["deltaz"+suffix]->Fill(r.z());
1161	<	h["pullx"+suffix]->Fill(r.x()/er.x());
1162	<	h["pully"+suffix]->Fill(r.y()/er.y());
1163	<	h["pullz"+suffix]->Fill(r.z()/er.z());
1164	<	h["errPVx"+suffix]->Fill(er.x());
1165	<	h["errPVy"+suffix]->Fill(er.y());
1166	<	h["errPVz"+suffix]->Fill(er.z());
1158	>	h_others->Fill1d("deltax"+suffix, r.x());
1159	>	h_others->Fill1d("deltay"+suffix, r.y());
1160	>	h_others->Fill1d("deltaz"+suffix, r.z());
1161	>	h_others->Fill1d("pullx"+suffix, r.x()/er.x());
1162	>	h_others->Fill1d("pully"+suffix, r.y()/er.y());
1163	>	h_others->Fill1d("pullz"+suffix, r.z()/er.z());
1164	>	h_others->Fill1d("errPVx"+suffix, er.x());
1165	>	h_others->Fill1d("errPVy"+suffix, er.y());
1166	>	h_others->Fill1d("errPVz"+suffix, er.z());
1167		}
1168
1169	+
1170		PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1171		map<int,double> data;
1172		data.clear();
#	Line 1236 \| Line 1184 \| PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1184		for ( int i = 0; i < 6; ++i) {
1185		switch (i) {
1186		case 0:
1187	<	fit(h["deltax"+suffix], fpar);
1187	>	if (!h_others->ReadHisto1D("deltax"+suffix)) break;
1188	>	fit(h_others->ReadHisto1D("deltax"+suffix), fpar);
1189		data[i] = fpar[0]*cm2um;
1190		break;
1191		case 1:
1192	<	fit(h["deltay"+suffix], fpar);
1192	>	if (!h_others->ReadHisto1D("deltay"+suffix)) break;
1193	>	fit(h_others->ReadHisto1D("deltay"+suffix), fpar);
1194		data[i] = fpar[0]*cm2um;
1195		break;
1196		case 2:
1197	<	fit(h["deltaz"+suffix], fpar);
1197	>	if (!h_others->ReadHisto1D("deltaz"+suffix)) break;
1198	>	fit(h_others->ReadHisto1D("deltaz"+suffix), fpar);
1199		data[i] = fpar[0]*cm2um;
1200		break;
1201		case 3:
1202	<	fit(h["pullx"+suffix], fpar);
1202	>	if (!h_others->ReadHisto1D("pullx"+suffix)) break;
1203	>	fit(h_others->ReadHisto1D("pullx"+suffix), fpar);
1204		data[i] = fpar[0];
1205		break;
1206		case 4:
1207	<	fit(h["pully"+suffix], fpar);
1207	>	if (!h_others->ReadHisto1D("pully"+suffix)) break;
1208	>	fit(h_others->ReadHisto1D("pully"+suffix), fpar);
1209		data[i] = fpar[0];
1210		break;
1211		case 5:
1212	<	fit(h["pullz"+suffix], fpar);
1212	>	if (!h_others->ReadHisto1D("pullz"+suffix)) break;
1213	>	fit(h_others->ReadHisto1D("pullz"+suffix), fpar);
1214		data[i] = fpar[0];
1215		break;
1216		}
1217	<	if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1217	>	if (data[i] > 0) edm::LogInfo("Analysis") << "[Analysis] data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1218		}
1219		}
1220		else
#	Line 1275 \| Line 1229 \| PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1229		ntk);
1230		}
1231
1232	<	void PVStudy::fit(TH1 *&hdist, double fitPars[]){
1232	>
1233	>	void PVStudy::fit(TH1 *hdist, double fitPars[]){
1234		TAxis *axis0 = hdist->GetXaxis();
1235		int nbin = axis0->GetLast();
1236		double nOF = hdist->GetBinContent(nbin+1);
1237		double nUF = hdist->GetBinContent(0);
1238		// double fitRange = axis0->GetBinUpEdge(nbin);
1239	<	double fitRange = axis0->GetXmax();
1239	>	// double fitRange = axis0->GetXmax();
1240	>	double fitRange = 2*hdist->GetRMS();
1241		double sigMax[2] = {0.,0.};
1242
1243	<	if ( verbose_ ){
1244	<	cout << "Last bin = " << nbin;
1245	<	cout << "; hdist: Overflow Entries = " << nOF;
1246	<	cout << "; Underflow Entries = " << nUF;
1247	<	cout << "; hdist->GetEntries() = " << hdist->GetEntries();
1248	<	cout << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1293	<	}
1243	>	edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
1244	>	<< "; hdist: Overflow Entries = " << nOF
1245	>	<< "; Underflow Entries = " << nUF
1246	>	<< "; hdist->GetEntries() = " << hdist->GetEntries()
1247	>	<< "; fitting range = " << -fitRange << " to " << fitRange << endl;
1248	>
1249		if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit
1250		for (int bi = 0; bi < nbin; bi++) {
1251		if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) {
#	Line 1302 \| Line 1257 \| *void PVStudy::fit(TH1 &hdist, double fi**
1257		if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1258		}
1259		}
1260	<	if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1260	>	//edm::LogInfo("Analysis") << "[Fit] Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1261
1262		TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1263		fgaus->SetParameter(1, 0.);
1264	<	fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
1264	>	fgaus->SetParLimits(1, -fitRange/10., fitRange/10.);
1265		fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
1266		fgaus->SetLineColor(4);
1267	<	hdist->Fit(fgaus,"Q");
1267	>	hdist->Fit(fgaus,"QLRM");
1268		gPad->Update();
1269		TPaveStats s = (TPaveStats) hdist->GetListOfFunctions()->FindObject("stats");
1270		s->SetOptStat(1111111);
1271		s->SetOptFit(0111);
1272		gPad->Update();
1273	<	if (verbose_) cout << "got function" << endl;
1273	>	//edm::LogInfo("Analysis") << "[Fit] got function" << endl;
1274		fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1275		}
1276		else
1277		fitPars[0] = -999;
1278		}
1279
1280	<	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype) {
1280	>
1281	>	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs) {
1282		string suffix;
1283		suffix = ""; // for unsplit tracks
1284		if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1285		if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1286
1287	<	h["nTrk"+suffix]->Fill(trackColl->size());
1287	>	h_pvtrk->Fill1d("nTrk"+suffix, trackColl->size());
1288		for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1289		itTrack != trackColl->end();
1290		++itTrack) {
1291	<	h["trkPt"+suffix]->Fill(itTrack->pt());
1292	<	h["trkDxy"+suffix]->Fill(itTrack->dxy());
1293	<	h["trkDz"+suffix]->Fill(itTrack->dz());
1294	<	h["trkEta"+suffix]->Fill(itTrack->eta());
1295	<	h["trkPhi"+suffix]->Fill(itTrack->phi());
1291	>	h_pvtrk->Fill1d("trkPt"+suffix, itTrack->pt());
1292	>	h_pvtrk->Fill1d("trkDxy"+suffix, itTrack->dxy());
1293	>	h_pvtrk->Fill1d("trkDxyCorr"+suffix, itTrack->dxy(bs));
1294	>	h_pvtrk->Fill1d("trkDz"+suffix, itTrack->dz());
1295	>	h_pvtrk->Fill1d("trkEta"+suffix, itTrack->eta());
1296	>	h_pvtrk->Fill1d("trkPhi"+suffix, itTrack->phi());
1297		}
1298		}
1299
1300	<	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
1300	>	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
1301		string suffix;
1302		suffix = ""; // for unsplit tracks
1303		if(datatype == 1) suffix = "1_spl"; // for splittracks 1
#	Line 1348 \| Line 1305 \| void PVStudy::fillTrackHistoInPV(const r
1305		int ivtx = 0;
1306		for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1307		v!=vertexColl->end(); ++v, ++ivtx) {
1308	<	if(fillNtuple) {
1309	<	if(datatype == 0) {
1310	<	nTrkPV_[ivtx] = v->tracksSize();
1311	<	sumptsq_[ivtx] = sumPtSquared(*v);
1312	<	isValid_[ivtx] = v->isValid();
1313	<	isFake_[ivtx] = v->isFake();
1314	<	recx_[ivtx] = v->x();
1315	<	recy_[ivtx] = v->y();
1316	<	recz_[ivtx] = v->z();
1317	<	recx_err_[ivtx] = v->xError();
1318	<	recy_err_[ivtx] = v->yError();
1319	<	recz_err_[ivtx] = v->zError();
1320	<
1321	<
1322	<	}
1323	<	if(datatype == 1) {
1324	<	nTrkPV1_spl_[ivtx] = v->tracksSize();
1325	<	sumptsq1_spl_[ivtx] = sumPtSquared(*v);
1326	<	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());
1308	>	if(!v->isFake()) {
1309	>	h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1310	>	h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1311	>	h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkRecVtx(*v));
1312	>	try {
1313	>	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1314	>	t!=vertexColl->begin()->tracks_end(); t++) {
1315	>	// illegal charge
1316	>	if ( (t).charge() < -1 \|\| (t).charge() > 1 ) {
1317	>	// h_pvtrk->Fill1d("trkPtPV", 0.);
1318	>	}
1319	>	else {
1320	>	h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1321	>	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1322	>	h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
1323	>	h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1324	>	h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1325	>	h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1326	>	}
1327		}
1328		}
1329	<	}
1330	<	catch (...) {
1331	<	// exception thrown when trying to use linked track
1332	<	// h["trkPtPV"]->Fill(0.);
1329	>	catch (...) {
1330	>	// exception thrown when trying to use linked track
1331	>	// h_pvtrk->Fill1d("trkPtPV", 0.);
1332	>	}
1333		}
1334		}
1429	–
1335		}
1336
1337		void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
#	Line 1445 \| Line 1350 \| void PVStudy::fillMCHisto(std::vector<si
1350		suffix = "_mct";
1351		nrecPV_mct_ = 0;
1352		}
1353	<
1353	>
1354		//========================================================
1355	<	// look for a matching reconstructed vertex in vertexColl
1356	<	//========================================================
1357	<
1355	>	// For each simulated vertex, look for a match in the vertexColl
1356	>	// If more than 1 recVtx is found, use the one with closest in Z
1357	>	//========================================================
1358	>
1359	>	// === look for a matching reconstructed vertex in vertexColl
1360		for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1361		vrec!=vtxColl->end(); ++vrec){
1455	–	int nrectrk = vrec->tracksSize();
1362		vsim->recVtx=NULL;
1363	<
1364	<	if(verbose_){
1365	<	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z << endl;
1366	<	cout << "Is matched? " << (matchVertex(vsim,vrec)?"Yes":"No") << endl;
1461	<	}
1462	<	if ( matchVertex(vsim,vrec) && vrec->isValid() && !vrec->isFake() ) {
1363	>	edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z << endl;
1364	>	edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(vsim,vrec)?"Yes":"No") << endl;
1365	>
1366	>	if ( matchVertex(vsim,vrec) ) {
1367		vsim->recVtx=&(*vrec);
1368	<
1465	<	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1368	>	//if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1369		//if( ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
1370		//\|\| (!vsim->recVtx) )
1371		//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;
1372		}
1373		}
1374	<	}
1374	>
1375	>	// === If match found fill the residual and pulls
1376	>	if(vsim->recVtx) {
1377	>	edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z << endl;
1378	>	// if study the resolution/pull for all tracks in the PVtx
1379	>	int nrectrk;
1380	>
1381	>	if(useHWTrk_)
1382	>	nrectrk=nHWTrkRecVtx(*(vsim->recVtx));
1383	>	else
1384	>	nrectrk = int(vsim->recVtx->tracksSize());
1385
1386	+	double fres_mct[3];
1387	+	double ferror_mct[3];
1388	+
1389	+	double ndofPV = vsim->recVtx->ndof();
1390	+	double normchi2PV = vsim->recVtx->normalizedChi2();
1391	+	double avgPtPV = avgPtRecVtx(*(vsim->recVtx));
1392	+
1393	+	fres_mct[0] = vsim->recVtx->x()-vsim->x;
1394	+	fres_mct[1] = vsim->recVtx->y()-vsim->y;
1395	+	fres_mct[2] = vsim->recVtx->z()-vsim->z;
1396	+	ferror_mct[0] = vsim->recVtx->xError();
1397	+	ferror_mct[1] = vsim->recVtx->yError();
1398	+	ferror_mct[2] = vsim->recVtx->zError();
1399	+
1400	+	h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1401	+	h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1402	+	h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1403	+	h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1404	+	h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1405	+	h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1406	+	h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1407	+	h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1408	+	h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1409	+	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1410	+	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1411	+	nrectrk));
1412	+	// Fill histo according to its track multiplicity
1413	+	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1414	+	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1415	+	nrectrk, datatype);
1416	+
1417	+	if(saventuple_) {
1418	+	//Fill the values for variables in pvtxtree_
1419	+	if(datatype == 1) {
1420	+	nTrkPV_spl1_mct_[nrecPV_spl1_mct_] = nrectrk;
1421	+	ndofPV_spl1_mct_[nrecPV_spl1_mct_] = ndofPV;
1422	+	normchi2PV_spl1_mct_[nrecPV_spl1_mct_] = normchi2PV;
1423	+	avgPtPV_spl1_mct_[nrecPV_spl1_mct_] = avgPtPV;
1424	+	deltax_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0];
1425	+	deltay_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1];
1426	+	deltaz_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[2];
1427	+	pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
1428	+	pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1429	+	pullz_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[2]/ferror_mct[2];
1430	+	errx_spl1_mct_[nrecPV_spl1_mct_] = ferror_mct[0];
1431	+	erry_spl1_mct_[nrecPV_spl1_mct_] = ferror_mct[1];
1432	+	errz_spl1_mct_[nrecPV_spl1_mct_] = ferror_mct[2];
1433	+	nrecPV_spl1_mct_++;
1434	+	}
1435	+	if(datatype == 2) {
1436	+	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1437	+	ndofPV_spl2_mct_[nrecPV_spl2_mct_] = ndofPV;
1438	+	normchi2PV_spl2_mct_[nrecPV_spl2_mct_] = normchi2PV;
1439	+	avgPtPV_spl2_mct_[nrecPV_spl2_mct_] = avgPtPV;
1440	+	deltax_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0];
1441	+	deltay_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1];
1442	+	deltaz_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[2];
1443	+	pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
1444	+	pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1445	+	pullz_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[2]/ferror_mct[2];
1446	+	errx_spl2_mct_[nrecPV_spl2_mct_] = ferror_mct[0];
1447	+	erry_spl2_mct_[nrecPV_spl2_mct_] = ferror_mct[1];
1448	+	errz_spl2_mct_[nrecPV_spl2_mct_] = ferror_mct[2];
1449	+	nrecPV_spl2_mct_++;
1450	+	}
1451	+	if(datatype == 3) {
1452	+	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1453	+	ndofPV_mct_[nrecPV_mct_] = ndofPV;
1454	+	normchi2PV_mct_[nrecPV_mct_] = normchi2PV;
1455	+	avgPtPV_mct_[nrecPV_mct_] = avgPtPV;
1456	+	deltax_mct_[nrecPV_mct_] = fres_mct[0];
1457	+	deltay_mct_[nrecPV_mct_] = fres_mct[0];
1458	+	deltaz_mct_[nrecPV_mct_] = fres_mct[0];
1459	+	pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
1460	+	pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1461	+	pullz_mct_[nrecPV_mct_] = fres_mct[2]/ferror_mct[2];
1462	+	errx_mct_[nrecPV_mct_] = ferror_mct[0];
1463	+	erry_mct_[nrecPV_mct_] = ferror_mct[1];
1464	+	errz_mct_[nrecPV_mct_] = ferror_mct[2];
1465	+	nrecPV_mct_++;
1466	+	}
1467	+	} // End of if(saventuple_) {
1468	+	} // if ( matchVertex(vsim,vrec) ) {
1469	+	else { // no rec vertex found for this simvertex
1470	+	edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1471	+	}
1472	+	}
1473
1474	+
1475		void PVStudy::SetVarToZero() {
1544	–	//Greg's variables
1476		fntrk_ = 0;
1477		//pvtx position (x,y,z) residual and error
1478		for(int i = 0; i<3;i++) {
1479	<	fres_[i] = 0;
1480	<	ferror_[i] = 0;
1479	>	fres_[i] = 0;
1480	>	ferror_[i] = 0;
1481		}
1482
1483	<	//Yanyan's variables
1484	<	// Number of reconstructed vertices
1483	>	// Event level information
1484	>	glob_runno_ = 0;
1485	>	glob_evtno_ = 0;
1486	>	glob_ls_ = 0;
1487	>	glob_bx_ = 0;
1488		nrecPV_ = 0;
1489		nrecPV1_spl_ = 0;
1490		nrecPV2_spl_ = 0;
#	Line 1566 \| Line 1500 \| void PVStudy::SetVarToZero() {
1500		min_ntrksep_ = 9999.0;
1501		min_sumpt2sep_ = -9999.0;
1502
1503	<
1503	>	// Variables filled per Vertex
1504		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	–
1505		// matched two-vertices
1506		nTrkPV1_spl_twovtx_[i] = 0;
1507	+	ndofPV1_spl_twovtx_[i] = 0;
1508	+	normchi2PV1_spl_twovtx_[i] = 0;
1509	+	avgPtPV1_spl_twovtx_[i] = 0;
1510	+	errx1_spl_twovtx_[i] = 0;
1511	+	erry1_spl_twovtx_[i] = 0;
1512	+	errz1_spl_twovtx_[i] = 0;
1513	+
1514		nTrkPV2_spl_twovtx_[i] = 0;
1515	<	nTrkPV_twovtx_[i] = 0;
1516	<	sumptsq1_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1517	<	sumptsq2_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1515	>	ndofPV2_spl_twovtx_[i] = 0;
1516	>	normchi2PV2_spl_twovtx_[i] = 0;
1517	>	avgPtPV2_spl_twovtx_[i] = 0;
1518	>	errx2_spl_twovtx_[i] = 0;
1519	>	erry2_spl_twovtx_[i] = 0;
1520	>	errz2_spl_twovtx_[i] = 0;
1521	>
1522		deltax_twovtx_[i] = 0;
1523		deltay_twovtx_[i] = 0;
1524		deltaz_twovtx_[i] = 0;
#	Line 1646 \| Line 1543 \| void PVStudy::SetVarToZero() {
1543		pullx_mct_[i] = 0;
1544		pully_mct_[i] = 0;
1545		pullz_mct_[i] = 0;
1546	<
1546	>	nTrkPV_mct_[i] = 0;
1547	>	ndofPV_mct_[i] = 0;
1548	>	normchi2PV_mct_[i] = 0;
1549	>	avgPtPV_mct_[i] = 0;
1550	>	errx_mct_[i] = 0;
1551	>	erry_mct_[i] = 0;
1552	>	errz_mct_[i] = 0;
1553	>
1554		deltax_spl1_mct_[i] = 0;
1555		deltay_spl1_mct_[i] = 0;
1556		deltaz_spl1_mct_[i] = 0;
1557		pullx_spl1_mct_[i] = 0;
1558		pully_spl1_mct_[i] = 0;
1559		pullz_spl1_mct_[i] = 0;
1560	<
1560	>	nTrkPV_spl1_mct_[i] = 0;
1561	>	ndofPV_spl1_mct_[i] = 0;
1562	>	normchi2PV_spl1_mct_[i] = 0;
1563	>	avgPtPV_spl1_mct_[i] = 0;
1564	>	errx_spl1_mct_[i] = 0;
1565	>	erry_spl1_mct_[i] = 0;
1566	>	errz_spl1_mct_[i] = 0;
1567	>
1568		deltax_spl2_mct_[i] = 0;
1569		deltay_spl2_mct_[i] = 0;
1570		deltaz_spl2_mct_[i] = 0;
1571		pullx_spl2_mct_[i] = 0;
1572		pully_spl2_mct_[i] = 0;
1573		pullz_spl2_mct_[i] = 0;
1574	+	nTrkPV_spl2_mct_[i] = 0;
1575	+	ndofPV_spl2_mct_[i] = 0;
1576	+	normchi2PV_spl2_mct_[i] = 0;
1577	+	avgPtPV_spl2_mct_[i] = 0;
1578	+	errx_spl2_mct_[i] = 0;
1579	+	erry_spl2_mct_[i] = 0;
1580	+	errz_spl2_mct_[i] = 0;
1581	+
1582		}
1664	–
1583		}
1584
1585		double PVStudy::sumPtSquared(const reco::Vertex & v) {
1586	<	double sum = 0.;
1587	<	double pT;
1588	<	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1589	<	pT = (**it).pt();
1590	<	sum += pT*pT;
1586	>	if(v.isFake() \|\| (!v.isValid()) \|\| v.tracksSize() == 0) return 0.0;
1587	>	else {
1588	>	double sum = 0.;
1589	>	double pT;
1590	>	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1591	>	pT = (**it).pt();
1592	>	sum += pT*pT;
1593	>	}
1594	>	return sum;
1595	>	}
1596	>	}
1597	>
1598	>
1599	>
1600	>	double PVStudy::avgPtRecVtx(const reco::Vertex & v) {
1601	>
1602	>	if(v.isFake() \|\| !v.isValid() \|\| v.tracksSize()==0 ) return 0;
1603	>	else {
1604	>	int nHWTrk = 0;
1605	>	double sumpT = 0.;
1606	>	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1607	>	if(v.trackWeight(*it) > 0.5 ) {
1608	>	sumpT += (**it).pt();
1609	>	nHWTrk++;
1610	>	}
1611	>	}
1612	>	if(nHWTrk > 0)
1613	>	return sumpT/double(nHWTrk);
1614	>	else
1615	>	return 0;
1616		}
1674	–	return sum;
1617		}
1618
1619	+	int PVStudy::nHWTrkRecVtx(const reco::Vertex & v) {
1620	+	int nHWTrkPV = 0;
1621	+	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1622	+	if(v.trackWeight(*it) > 0.5)
1623	+	nHWTrkPV++;
1624	+	}
1625	+	return nHWTrkPV;
1626	+	}

Diff Legend

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

Comparing UserCode/Jeng/PVStudy/plugins/PVStudy.cc (file contents): Revision 1.12 by yygao, Fri Nov 27 12:31:05 2009 UTC vs. Revision 1.21 by yygao, Sun Feb 28 00:07:09 2010 UTC

Diff Legend

Comparing UserCode/Jeng/PVStudy/plugins/PVStudy.cc (file contents):
Revision 1.12 by yygao, Fri Nov 27 12:31:05 2009 UTC vs.
Revision 1.21 by yygao, Sun Feb 28 00:07:09 2010 UTC