ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

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

Comparing UserCode/Jeng/PVStudy/plugins/PVStudy.cc (file contents):
Revision 1.5 by jengbou, Wed Sep 9 05:56:50 2009 UTC vs.
Revision 1.11 by yygao, Wed Nov 11 19:28:50 2009 UTC

#	Line 12 | Line 12 | Implementation:
12		*/
13		//
14		// Original Author:  "Geng-yuan Jeng/UC Riverside"
15	+	//                    "Yanyan Gao/Fermilab ygao@fnal.gov"
16		//         Created:  Thu Aug 20 11:55:40 CDT 2009
17		// $Id$
18		//
#	Line 24 | Line 25 | Implementation:
25		#include <vector>
26		#include <iostream>
27		#include <sstream>
28	+	#include <algorithm>
29
30		// user include files
31		#include "FWCore/Framework/interface/Frameworkfwd.h"
#	Line 59 | Line 61 | Implementation:
61
62		using namespace std;
63
64	<	PVStudy::PVStudy(const edm::ParameterSet& iConfig):
63	<	nTrkMin_(10),nTrkMax_(100)
64	>	PVStudy::PVStudy(const edm::ParameterSet& iConfig)
65		{
66	+	//=======================================================================
67	+	// Get configuration for TrackTupleMaker
68	+	//=======================================================================
69		simG4_           = iConfig.getParameter<edm::InputTag>( "simG4" );
70	<	tracksLabel_     = iConfig.getUntrackedParameter<edm::InputTag>("tracks");
71	<	vtxSample_       = iConfig.getUntrackedParameter<edm::InputTag>("vtxSample");
72	<	verbose_         = iConfig.getUntrackedParameter<bool>("verbose", false);
70	>	trackCollectionTag_     = iConfig.getUntrackedParameter<edm::InputTag>("trackCollection");
71	>	splitTrackCollection1Tag_     = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection1");
72	>	splitTrackCollection2Tag_     = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
73	>	vertexCollectionTag_       = iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection");
74	>	splitVertexCollection1Tag_       = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
75	>	splitVertexCollection2Tag_       = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
76	>	pixelVertexCollectionTag_      = iConfig.getParameter<edm::InputTag>("pixelVertexCollectionTag");
77	>	verbose_         = iConfig.getUntrackedParameter<bool>("verbose",false);
78		realData_        = iConfig.getUntrackedParameter<bool>("realData",false);
79		analyze_         = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
80	<	outputfilename_  = iConfig.getUntrackedParameter<string>("OutputFileName");
80	>	saventuple_         = iConfig.getUntrackedParameter<bool>("saventuple",false);
81	>	outputfilename_  = iConfig.getUntrackedParameter<string>("OutputFileName");
82	>	ntrkdiffcut_         = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
83	>	zsigncut_         = iConfig.getUntrackedParameter<int>("zsigncut");
84	>	nTrkMin_ =  iConfig.getUntrackedParameter<int>("nTrkMin");
85	>	nTrkMax_ =  iConfig.getUntrackedParameter<int>("nTrkMax");
86
87		//now do what ever initialization is needed
88		pvinfo.clear();
89
76	–	file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
77	–	ftree_ = new TTree("mytree","mytree");
78	–	ftree_->AutoSave();
79	–	ftree_->Branch("residual",&fres,"fres[3]/D");
80	–	ftree_->Branch("error",&ferror,"ferror[3]/D");
81	–	ftree_->Branch("nTrk",&fntrk,"fntrk/I");
90
91	+	// Specify the data mode vector
92	+	if(realData_) datamode.push_back(0);
93	+	else {
94	+	datamode.push_back(0);
95	+	datamode.push_back(1);
96	+	datamode.push_back(2);
97	+	datamode.push_back(3);
98	+	}
99	+
100	+	// Create ntuple files if needed
101	+	if(saventuple_) {
102	+	file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
103	+	ftree_ = new TTree("mytree","mytree");
104	+	ftree_->AutoSave();
105	+	ftree_->Branch("residual",&fres_,"fres_[3]/D");
106	+	ftree_->Branch("error",&ferror_,"ferror_[3]/D");
107	+	ftree_->Branch("nTrkPV",&fntrk_,"fntrk_/I");
108	+
109	+	// pvtxtree_ analyzing the pvtxs ootb
110	+	pvtxtree_ = new TTree("pvtxtree","pvtxtree");
111	+	//pvtx using all tracks
112	+	pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
113	+	pvtxtree_->Branch("nTrkPV",&nTrkPV_,"nTrkPV[nrecPV]/I");
114	+	pvtxtree_->Branch("sumptsq",&sumptsq_,"sumptsq[nrecPV]/D");
115	+	pvtxtree_->Branch("isValid",&isValid_,"isValid/I");
116	+	pvtxtree_->Branch("isFake",&isFake_,"isFake/I");
117	+	pvtxtree_->Branch("recx",&recx_,"recx[nrecPV]/D");
118	+	pvtxtree_->Branch("recy",&recy_,"recy[nrecPV]/D");
119	+	pvtxtree_->Branch("recz",&recz_,"recz[nrecPV]/D");
120	+	pvtxtree_->Branch("recx_err",&recx_err_,"recx_err[nrecPV]/D");
121	+	pvtxtree_->Branch("recy_err",&recy_err_,"recy_err[nrecPV]/D");
122	+	pvtxtree_->Branch("recz_err",&recz_err_,"recz_err[nrecPV]/D");
123	+
124	+	pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
125	+	pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
126	+	pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
127	+
128	+	//pvtx using splittracks1
129	+	pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
130	+	pvtxtree_->Branch("nTrkPV1_spl",&nTrkPV1_spl_,"nTrkPV1_spl[nrecPV1_spl]/I");
131	+	pvtxtree_->Branch("sumptsq1_spl",&sumptsq1_spl_,"sumptsq1_spl[nrecPV1_spl]/D");
132	+	pvtxtree_->Branch("isValid1_spl",&isValid1_spl_,"isValid1_spl/I");
133	+	pvtxtree_->Branch("isFake1_spl",&isFake1_spl_,"isFake1_spl/I");
134	+	pvtxtree_->Branch("recx1_spl",&recx1_spl_,"recx1_spl[nrecPV1_spl]/D");
135	+	pvtxtree_->Branch("recy1_spl",&recy1_spl_,"recy1_spl[nrecPV1_spl]/D");
136	+	pvtxtree_->Branch("recz1_spl",&recz1_spl_,"recz1_spl[nrecPV1_spl]/D");
137	+	pvtxtree_->Branch("recx1_err_spl",&recx1_err_spl_,"recx1_err_spl[nrecPV1_spl]/D");
138	+	pvtxtree_->Branch("recy1_err_spl",&recy1_err_spl_,"recy1_err_spl[nrecPV1_spl]/D");
139	+	pvtxtree_->Branch("recz1_err_spl",&recz1_err_spl_,"recz1_err_spl[nrecPV1_spl]/D");
140	+
141	+	//pvtx using splittracks2
142	+	pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
143	+	pvtxtree_->Branch("nTrkPV2_spl",&nTrkPV2_spl_,"nTrkPV2_spl[nrecPV2_spl]/I");
144	+	pvtxtree_->Branch("sumptsq2_spl",&sumptsq2_spl_,"sumptsq2_spl[nrecPV2_spl]/D");
145	+	pvtxtree_->Branch("isValid2_spl",&isValid2_spl_,"isValid2_spl/I");
146	+	pvtxtree_->Branch("isFake2_spl",&isFake2_spl_,"isFake2_spl/I");
147	+	pvtxtree_->Branch("recx2_spl",&recx2_spl_,"recx2_spl[nrecPV2_spl]/D");
148	+	pvtxtree_->Branch("recy2_spl",&recy2_spl_,"recy2_spl[nrecPV2_spl]/D");
149	+	pvtxtree_->Branch("recz2_spl",&recz2_spl_,"recz2_spl[nrecPV2_spl]/D");
150	+	pvtxtree_->Branch("recx2_err_spl",&recx2_err_spl_,"recx2_err_spl[nrecPV2_spl]/D");
151	+	pvtxtree_->Branch("recy2_err_spl",&recy2_err_spl_,"recy2_err_spl[nrecPV2_spl]/D");
152	+	pvtxtree_->Branch("recz2_err_spl",&recz2_err_spl_,"recz2_err_spl[nrecPV2_spl]/D");
153	+
154	+	//pixeVertices
155	+	pvtxtree_->Branch("nrecPV_pxlpvtx",&nrecPV_pxlpvtx_,"nrecPV_pxlpvtx/I");
156	+	pvtxtree_->Branch("nTrkPV_pxlpvtx",&nTrkPV_pxlpvtx_,"nTrkPV_pxlpvtx[nrecPV_pxlpvtx]/I");
157	+	pvtxtree_->Branch("sumptsq_pxlpvtx",&sumptsq_pxlpvtx_,"sumptsq_pxlpvtx[nrecPV_pxlpvtx]/D");
158	+	pvtxtree_->Branch("isValid_pxlpvtx",&isValid_pxlpvtx_,"isValid_pxlpvtx/I");
159	+	pvtxtree_->Branch("isFake_pxlpvtx",&isFake_pxlpvtx_,"isFake_pxlpvtx/I");
160	+	pvtxtree_->Branch("recx_pxlpvtx",&recx_pxlpvtx_,"recx_pxlpvtx[nrecPV_pxlpvtx]/D");
161	+	pvtxtree_->Branch("recy_pxlpvtx",&recy_pxlpvtx_,"recy_pxlpvtx[nrecPV_pxlpvtx]/D");
162	+	pvtxtree_->Branch("recz_pxlpvtx",&recz_pxlpvtx_,"recz_pxlpvtx[nrecPV_pxlpvtx]/D");
163	+	pvtxtree_->Branch("recx_err_pxlpvtx",&recx_err_pxlpvtx_,"recx_err_pxlpvtx[nrecPV_pxlpvtx]/D");
164	+	pvtxtree_->Branch("recy_err_pxlpvtx",&recy_err_pxlpvtx_,"recy_err_pxlpvtx[nrecPV_pxlpvtx]/D");
165	+	pvtxtree_->Branch("recz_err_pxlpvtx",&recz_err_pxlpvtx_,"recz_err_pxlpvtx[nrecPV_pxlpvtx]/D");
166	+
167	+	//Fill the variables in the twovtx pair (recvtx1, recvtx2)
168	+	pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
169	+	pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
170	+	pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
171	+	pvtxtree_->Branch("sumptsq1_spl_twovtx",&sumptsq1_spl_twovtx_,"sumptsq1_spl_twovtx[nrecPV_twovtx]/D");
172	+	pvtxtree_->Branch("sumptsq2_spl_twovtx",&sumptsq2_spl_twovtx_,"sumptsq2_spl_twovtx[nrecPV_twovtx]/D");
173	+	pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
174	+	pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
175	+	pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
176	+	pvtxtree_->Branch("deltaz_twovtx",&deltaz_twovtx_,"deltaz_twovtx[nrecPV_twovtx]/D");
177	+	pvtxtree_->Branch("errx_twovtx",&errx_twovtx_,"errx_twovtx[nrecPV_twovtx]/D");
178	+	pvtxtree_->Branch("erry_twovtx",&erry_twovtx_,"erry_twovtx[nrecPV_twovtx]/D");
179	+	pvtxtree_->Branch("errz_twovtx",&errz_twovtx_,"errz_twovtx[nrecPV_twovtx]/D");
180	+	pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
181	+	pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
182	+	pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
183	+
184	+	// MC variables
185	+	if(!realData_) {
186	+	pvtxtree_->Branch("nsimPV",&nsimPV_,"nsimPV/I");
187	+	pvtxtree_->Branch("nsimTrkPV",&nsimTrkPV_,"nsimTrkPV[nsimPV]/I");
188	+	pvtxtree_->Branch("simx",&simx_,"simx[nsimPV]/D");
189	+	pvtxtree_->Branch("simy",&simy_,"simy[nsimPV]/D");
190	+	pvtxtree_->Branch("simz",&simz_,"simz[nsimPV]/D");
191	+	pvtxtree_->Branch("simptsq",&simptsq_,"simptsq[nsimPV]/D");
192	+
193	+	// For pvtxs with all the tracks
194	+	pvtxtree_->Branch("nrecPV_mct",&nrecPV_mct_,"nrecPV_mct/I");
195	+	pvtxtree_->Branch("deltax_mct",&deltax_mct_,"deltax_mct[nrecPV_mct]/D");
196	+	pvtxtree_->Branch("deltay_mct",&deltay_mct_,"deltay_mct[nrecPV_mct]/D");
197	+	pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
198	+	pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
199	+	pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
200	+	pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
201	+	// For pvtxs with splittracks1
202	+	pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
203	+	pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
204	+	pvtxtree_->Branch("deltay_spl1_mct",&deltay_spl1_mct_,"deltay_spl1_mct[nrecPV_spl1_mct]/D");
205	+	pvtxtree_->Branch("deltaz_spl1_mct",&deltaz_spl1_mct_,"deltaz_spl1_mct[nrecPV_spl1_mct]/D");
206	+	pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
207	+	pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
208	+	pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
209	+	// For pvtxs with splittracks1
210	+	pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
211	+	pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
212	+	pvtxtree_->Branch("deltay_spl2_mct",&deltay_spl2_mct_,"deltay_spl2_mct[nrecPV_spl2_mct]/D");
213	+	pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
214	+	pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
215	+	pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
216	+	pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
217	+	}
218	+	}
219	+
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();
226	<
227	<	h["nTrk"]           = subDir.make<TH1F>("nTrk","Num of total tracks",300,0,300);
228	<	h["nTrkPV"]         = subDir.make<TH1F>("nTrkPV","Num of Tracks from PV",300,0,300);
229	<	h["trkPt"]          = subDir.make<TH1D>("trkPt","Pt of all tracks",100,0,100);
230	<	h["trkPtPV"]        = subDir.make<TH1D>("trkPtPV","Pt dist of tracks from PV",100,0,100);
231	<
232	<	if (!realData_) {
233	<	h["deltax"]       = subDir.make<TH1D>("deltax","x-distance (Rec - Sim)",800,-0.04,0.04);
234	<	h["deltay"]       = subDir.make<TH1D>("deltay","y-distance (Rec - Sim)",800,-0.04,0.04);
235	<	h["deltaz"]       = subDir.make<TH1D>("deltaz","z-distance (Rec - Sim)",800,-0.04,0.04);
236	<	h["pullx"]        = subDir.make<TH1D>("pullx","x-pull (Rec - Sim)",1000,-25.,25.);
237	<	h["pully"]        = subDir.make<TH1D>("pully","y-pull (Rec - Sim)",1000,-25.,25.);
238	<	h["pullz"]        = subDir.make<TH1D>("pullz","z-pull (Rec - Sim)",1000,-25.,25.);
239	<	h["errPVx"]       = subDir.make<TH1D>("errPVx","X vertex error",100,0.,0.02);
240	<	h["errPVy"]       = subDir.make<TH1D>("errPVy","Y vertex error",100,0.,0.02);
241	<	h["errPVz"]       = subDir.make<TH1D>("errPVz","Z vertex error",100,0.,0.02);
242	<
243	<	// Summary of analysis:
244	<	if (analyze_) {
245	<	h2["resx_nTrk"]= subDir.make<TH2D>("resx_ntrk","x-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
246	<	h2["resx_nTrk"]->SetMarkerStyle(21);
247	<	h2["resx_nTrk"]->SetMarkerColor(4);
248	<	h2["resx_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
249	<	h2["resx_nTrk"]->GetYaxis()->SetTitle("#mum");
250	<	h2["resy_nTrk"]= subDir.make<TH2D>("resy_ntrk","y-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
251	<	h2["resy_nTrk"]->SetMarkerStyle(21);
252	<	h2["resy_nTrk"]->SetMarkerColor(4);
253	<	h2["resy_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
254	<	h2["resy_nTrk"]->GetYaxis()->SetTitle("#mum");
255	<	h2["resz_nTrk"]= subDir.make<TH2D>("resz_ntrk","z-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
256	<	h2["resz_nTrk"]->SetMarkerStyle(21);
120	<	h2["resz_nTrk"]->SetMarkerColor(4);
121	<	h2["resz_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
122	<	h2["resz_nTrk"]->GetYaxis()->SetTitle("#mum");
123	<	h2["pullx_nTrk"]= subDir.make<TH2D>("pullx_ntrk","x-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
124	<	h2["pullx_nTrk"]->SetMarkerStyle(21);
125	<	h2["pullx_nTrk"]->SetMarkerColor(4);
126	<	h2["pullx_nTrk"]->SetBit(TH1::kCanRebin);
127	<	h2["pullx_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
128	<	h2["pully_nTrk"]= subDir.make<TH2D>("pully_ntrk","y-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
129	<	h2["pully_nTrk"]->SetMarkerStyle(21);
130	<	h2["pully_nTrk"]->SetMarkerColor(4);
131	<	h2["pully_nTrk"]->SetBit(TH1::kCanRebin);
132	<	h2["pully_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
133	<	h2["pullz_nTrk"]= subDir.make<TH2D>("pullz_ntrk","z-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
134	<	h2["pullz_nTrk"]->SetMarkerStyle(21);
135	<	h2["pullz_nTrk"]->SetMarkerColor(4);
136	<	h2["pullz_nTrk"]->SetBit(TH1::kCanRebin);
137	<	h2["pullz_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
226	>
227	>	//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_
231	>
232	>
233	>	for(int i=0;i<3;i++)
234	>	{
235	>	string suffix;
236	>	if(i == 0) suffix = "";
237	>	else {
238	>	stringstream ss;
239	>	ss << i;
240	>	suffix =ss.str()+"_spl";
241	>	}
242	>	h["nTrk"+suffix]   = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
243	>	h["trkPt"+suffix]  = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
244	>	h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#eta of rec tracks "+suffix),100,-3,3);
245	>	h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
246	>	h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-0.5,0.5);
247	>	h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
248	>
249	>	h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
250	>	h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
251	>	h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#eta of rec tracks in PV"+suffix),100,-3,3);
252	>	h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
253	>	h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);
254	>	h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
255	>	h["nrecPV"+suffix]   = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);
256	>	suffix.clear();
257		}
258	+	h["nrecPVDiff"]     = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
259	+	h["nTrkPVDiff"]     = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
260	+	h["nTrkPVRelDiff"]  = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
261	+
262	+	// Histograms on comparing the multi-vertices
263	+	// Difference in reconstructed vtx position
264	+	h["min_xsep"]   = subDir.make<TH1D>("min_xsep", "min x diff of primary and secondary pvtx",300,0,0.1);
265	+	h["min_xsep_sign"]   = subDir.make<TH1D>("min_xsep_sign", "min x diff in signf of primary and secondary pvtx",300,0,5);
266	+	h["min_ysep"]   = subDir.make<TH1D>("min_ysep", "min y diff of primary and secondary pvtx",300,0,0.1);
267	+	h["min_ysep_sign"]   = subDir.make<TH1D>("min_ysep_sign", "min y diff in signf of primary and secondary pvtx",300,0,5);
268	+	h["min_zsep"]   = subDir.make<TH1D>("min_zsep", "min z diff of primary and secondary pvtx",300,0,5);
269	+	h["min_zsep_sign"]   = subDir.make<TH1D>("min_zsep_sign", "min z diff in signf of primary and secondary pvtx",300,0,200);
270	+	// Difference in reconstructed vtx position
271	+	h["min_ntrksep"]   = subDir.make<TH1D>("min_ntrksep", "min nTrk diff of primary and secondary pvtx",201,-50.5,150.5);
272	+	h["min_sumpt2sep"]   = subDir.make<TH1D>("min_sumpt2sep", "min sumpt2 diff of primary and secondary pvtx",300,0,10000);
273	+
274	+	//Book histograms sensitive to data/mc
275	+	for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
276	+	string suffix;
277	+	if(verbose_)  cout<<"datamode = "<< *it<<endl;
278	+	switch(*it) {
279	+	case 0: suffix = "";
280	+	break;
281	+	case 1: suffix = "_spl1_mct";
282	+	break;
283	+	case 2: suffix = "_spl2_mct";
284	+	break;
285	+	case 3: suffix = "_mct";
286	+	break;
287	+	}
288	+	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);
297	+
298	+	// Book residual, error and pull histograms for each nTrk bin
299		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
300		stringstream ss;
301		ss << ntrk;
302	<	string suffix = ss.str();
303	<	h["resx_"+suffix]  = fs->make<TH1D>(TString("deltax_"+suffix),"x-distance (Rec - Sim)",100,-0.02,0.02);
304	<	h["resx_"+suffix]->GetXaxis()->SetTitle("cm");
305	<	h["resy_"+suffix]  = fs->make<TH1D>(TString("deltay_"+suffix),"y-distance (Rec - Sim)",100,-0.02,0.02);
306	<	h["resy_"+suffix]->GetXaxis()->SetTitle("cm");
307	<	h["resz_"+suffix]  = fs->make<TH1D>(TString("deltaz_"+suffix),"z-distance (Rec - Sim)",100,-0.02,0.02);
308	<	h["resz_"+suffix]->GetXaxis()->SetTitle("cm");
309	<	h["pullx_"+suffix] = fs->make<TH1D>(TString("pullx_"+suffix),"x-pull (Rec - Sim)",100,-10.,10.);
310	<	h["pully_"+suffix] = fs->make<TH1D>(TString("pully_"+suffix),"y-pull (Rec - Sim)",100,-10.,10.);
311	<	h["pullz_"+suffix] = fs->make<TH1D>(TString("pullz_"+suffix),"z-pull (Rec - Sim)",100,-10.,10.);
312	<	h["errPVx_"+suffix]  = subDir1.make<TH1D>(TString("errPVx_"+suffix),"X vertex error",100,0.,0.02);
313	<	h["errPVy_"+suffix]  = subDir1.make<TH1D>(TString("errPVy_"+suffix),"Y vertex error",100,0.,0.02);
314	<	h["errPVz_"+suffix]  = subDir1.make<TH1D>(TString("errPVz_"+suffix),"Z vertex error",100,0.,0.02);
315	<	suffix.clear();
316	<	}
302	>	string suffix2 = suffix+"_"+ss.str();
303	>	h["deltax"+suffix2]  = subDir1.make<TH1D>(TString("deltax"+suffix2), TString("x-residual of pvtx"+suffix),100,-0.02,0.02);
304	>	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	>
318	>	// Book residual and pull histograms only when analyzeOntheFly is enabled
319	>	if(analyze_) {
320	>	h2["resx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resx"+suffix+"_nTrk"), TString("x-resolution vs number of tracks in pvtx"+suffix),150,0,150,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),150,0,150,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),150,0,150,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),150,0,150,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),150,0,150,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("x-pull vs number of tracks"+suffix),150,0,150,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_) {
351	>	suffix.clear();
352	>	} // End of Book histograms sensitive to data/mc
353	>
354	>	// Book histograms about pixelVertices
355	>	h["trkdz_pxlpvtxdz"]   = subDir.make<TH1D>("trkdz_pxlpvtxdz", "(Track dz - pixelpvtx dz) in cm",300,-0.5,0.5);
356	>	h["trkdz_pxlpvtxdz_pxlpvtxdzerr"]   = subDir.make<TH1D>("trkdz_pxlpvtxdz_pxlpvtxdzerr", "|Track dz - pixelpvtx dz| / pxlpvtxdzErr",300,0,100);
357	>	h["trkdz_pxlpvtxdz_trkdzerr"]   = subDir.make<TH1D>("trkdz_pxlpvtxdz_trkdzerr", "|Track dz - pixelpvtx dz| / trkdzErr",300,0,50);
358	>	h["trkdzErr_pxlpvtx"]   = subDir.make<TH1D>("trkdzErr_pxlpvtxdz", "Track dzErr of leading pixelpvtx ",300,0,0.5);
359	>	h["trkdzErr_pvtx"]   = subDir.make<TH1D>("trkdzErr_pvtx", "Track dzErr of the leading pvtx ",300,0,0.5);
360	>	h["dzErr_pxlpvtx"]   = subDir.make<TH1D>("dzErr_pxlpvtx", "zError of the leading pvtx ",300,0,0.5);
361	>	// Compare offlinePrimaryVertices with pixelVertices
362	>	h["nrecPV_minus_nrecPxlPV"] =  subDir.make<TH1D>("nrecPV_minus_nrecPxlPV", "nrecPV_minus_nrecPxlPV",21,-10.5,10.5);
363	>	h["recxPV_minus_recxPxlPV"] =  subDir.make<TH1D>("recxPV_minus_recxPxlPV", "recxPV_minus_recxPxlPV",300,-0.02,0.02);
364	>	h["recyPV_minus_recyPxlPV"] =  subDir.make<TH1D>("recyPV_minus_recyPxlPV", "recyPV_minus_recyPxlPV",300,-0.02,0.02);
365	>	h["reczPV_minus_reczPxlPV"] =  subDir.make<TH1D>("reczPV_minus_reczPxlPV", "reczPV_minus_reczPxlPV",300,-0.1,0.1);
366	>
367	>	// Book MC only plots
368	>	if (!realData_) {
369	>	h["genPart_T"]      = subDir.make<TH1D>("genPart_T","t component of gen particles",300,-0.5,0.5);
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		}
375		}
376
160	–
377		PVStudy::~PVStudy()
378		{
379
#	Line 210 | Line 426 | std::vector<PVStudy::simPrimaryVertex> P
426		cout << "number of vertices " << evt->vertices_size() << endl;
427		}
428
429	<	int idx=0;
429	>	int idx=0; int npv=0;
430		for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
431	<	vitr != evt->vertices_end(); ++vitr )
432	<	{ // loop for vertex ...
433	<	HepMC::FourVector pos = (*vitr)->position();
434	<	//HepLorentzVector pos = (*vitr)->position();
435	<	if (pos.t()>0) { continue;} // ?
436	<
437	<	bool hasMotherVertex=false;
438	<	//      if(verbose_) cout << "mothers" << endl;
439	<	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
440	<	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
441	<	HepMC::GenVertex * mv=(*mother)->production_vertex();
442	<	if (mv) {hasMotherVertex=true;}
443	<	//      if(verbose_) {
444	<	//        cout << "\t";
229	<	//        (*mother)->print();
230	<	//      }
231	<	}
232	<
233	<	if(hasMotherVertex) {continue;}
234	<
235	<	// could be a new vertex, check  all primaries found so far to avoid multiple entries
236	<	const double mm2cm=0.1;
237	<	simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
238	<	simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
239	<	for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
240	<	v0!=simpv.end(); v0++){
241	<	if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
242	<	vp=&(*v0);
243	<	break;
431	>	vitr != evt->vertices_end(); ++vitr ) { // loop for vertex ...
432	>	HepMC::FourVector pos = (*vitr)->position();
433	>	//HepLorentzVector pos = (*vitr)->position();
434	>
435	>	// t component of PV:
436	>	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
437	>	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
438	>	//        std::cout << "Status = " << (*mother)->status() << std::endl;
439	>	HepMC::GenVertex * mv=(*mother)->production_vertex();
440	>	if( ((*mother)->status() == 3) && (!mv)) {
441	>	//        std::cout << "npv= " << npv << std::endl;
442	>	if (npv == 0) {
443	>	h["genPart_cT"]->Fill(pos.t()); // mm
444	>	h["genPart_T"]->Fill(pos.t()/299.792458); // ns
445		}
446	+	npv++;
447	+	}
448	+	}
449	+	//       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
450	+
451	+	bool hasMotherVertex=false;
452	+	if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl;
453	+	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
454	+	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
455	+	HepMC::GenVertex * mv=(*mother)->production_vertex();
456	+	//        std::cout << "Status = " << (*mother)->status() << std::endl;
457	+	if (mv) {
458	+	hasMotherVertex=true;
459	+	if(!verbose_) break; //if verbose_, print all particles of gen vertices
460	+	}
461	+	if(verbose_) {
462	+	cout << "\t";
463	+	(*mother)->print();
464		}
465	+	}
466
467	<	if(!vp){
468	<	// this is a new vertex
469	<	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
470	<	simpv.push_back(sv);
471	<	vp=&simpv.back();
472	<	}else{
473	<	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
474	<	}
475	<	vp->genVertex.push_back((*vitr)->barcode());
476	<	// collect final state descendants
477	<	for ( HepMC::GenVertex::particle_iterator
478	<	daughter  = (*vitr)->particles_begin(HepMC::descendants);
479	<	daughter != (*vitr)->particles_end(HepMC::descendants);
480	<	++daughter ) {
481	<	if (isFinalstateParticle(*daughter)){
482	<	if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
483	<	== vp->finalstateParticles.end()){
484	<	vp->finalstateParticles.push_back((*daughter)->barcode());
485	<	HepMC::FourVector m=(*daughter)->momentum();
486	<	// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
487	<	// but adding FourVectors seems not to be foreseen
488	<	vp->ptot.setPx(vp->ptot.px()+m.px());
489	<	vp->ptot.setPy(vp->ptot.py()+m.py());
490	<	vp->ptot.setPz(vp->ptot.pz()+m.pz());
491	<	vp->ptot.setE(vp->ptot.e()+m.e());
492	<	vp->ptsq+=(m.perp())*(m.perp());
493	<	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
494	<	vp->nGenTrk++;
495	<	}
467	>	if(hasMotherVertex) continue;
468	>
469	>	// could be a new vertex, check  all primaries found so far to avoid multiple entries
470	>	const double mm2cm=0.1;
471	>	simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
472	>	simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
473	>	for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
474	>	v0!=simpv.end(); v0++){
475	>	if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
476	>	vp=&(*v0);
477	>	break;
478	>	}
479	>	}
480	>
481	>	if(!vp){
482	>	// this is a new vertex
483	>	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
484	>	simpv.push_back(sv);
485	>	vp=&simpv.back();
486	>	}else{
487	>	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
488	>	}
489	>	vp->genVertex.push_back((*vitr)->barcode());
490	>	// collect final state descendants
491	>	for ( HepMC::GenVertex::particle_iterator daughter  = (*vitr)->particles_begin(HepMC::descendants);
492	>	daughter != (*vitr)->particles_end(HepMC::descendants);
493	>	++daughter ) {
494	>	if (isFinalstateParticle(*daughter)){
495	>	if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
496	>	== vp->finalstateParticles.end()){
497	>	vp->finalstateParticles.push_back((*daughter)->barcode());
498	>	HepMC::FourVector m=(*daughter)->momentum();
499	>	// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
500	>	// but adding FourVectors seems not to be foreseen
501	>	vp->ptot.setPx(vp->ptot.px()+m.px());
502	>	vp->ptot.setPy(vp->ptot.py()+m.py());
503	>	vp->ptot.setPz(vp->ptot.pz()+m.pz());
504	>	vp->ptot.setE(vp->ptot.e()+m.e());
505	>	vp->ptsq+=(m.perp())*(m.perp());
506	>	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
507	>	vp->nGenTrk++;
508		}
509		}
510		}
279	–	idx++;
511		}
512	+	}
513		}
514		return simpv;
515		}
#	Line 357 | Line 589 | void
589		PVStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
590		{
591		using namespace edm;
592	<	using reco::TrackCollection;
592	>	using namespace reco;
593	>
594	>	if (verbose_)
595	>	cout<<"PVStudy::analyze():"<<endl;
596	>	//=======================================================
597	>	// Initialize Root-tuple variables if needed
598	>	//=======================================================
599	>	//if(saventuple_)
600	>	SetVarToZero();
601	>
602	>	//=======================================================
603	>	// Track accessors
604	>	//=======================================================
605	>	//trackCollection
606	>	static const reco::TrackCollection s_empty_trackColl;
607	>	const reco::TrackCollection *trackColl = &s_empty_trackColl;
608	>	edm::Handle<reco::TrackCollection> trackCollectionHandle;
609	>	iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle);
610	>	if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
611	>	trackColl = trackCollectionHandle.product();
612	>	} else {
613	>	cout << "trackCollection cannot be found -> using empty collection of same type." <<endl;
614	>	}
615
616	<	ftree_->SetBranchAddress("residual",&fres);
617	<	ftree_->SetBranchAddress("error",&ferror);
618	<	ftree_->SetBranchAddress("nTrk",&fntrk);
616	>	//splitTrackCollection1
617	>	static const reco::TrackCollection s_empty_splitTrackColl1;
618	>	const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
619	>	edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
620	>	iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle);
621	>	if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
622	>	splitTrackColl1 = splitTrackCollection1Handle.product();
623	>	} else {
624	>	cout << "splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
625	>	}
626
627	<	Handle<TrackCollection> tracks;
628	<	iEvent.getByLabel(tracksLabel_,tracks);
627	>	//splitTrackCollection2
628	>	static const reco::TrackCollection s_empty_splitTrackColl2;
629	>	const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
630	>	edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
631	>	iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle);
632	>	if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
633	>	splitTrackColl2 = splitTrackCollection2Handle.product();
634	>	} else {
635	>	cout << "splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
636	>	}
637
369	–	Handle<reco::VertexCollection> recVtxs;
370	–	iEvent.getByLabel(vtxSample_, recVtxs);
638
639	+	//=======================================================
640	+	// Fill trackparameters of the input tracks to pvtx fitter
641	+	//=======================================================
642	+	if(verbose_)
643	+	cout<<"Start filling track parameters of the input tracks to pvtx fitter."<<endl;
644	+	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype)
645	+	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
646	+	fillTrackHisto(trackColl, 0);
647	+	fillTrackHisto(splitTrackColl1, 1);
648	+	fillTrackHisto(splitTrackColl2, 2);
649	+	if(verbose_)
650	+	cout<<"End filling track parameters of the input tracks to pvtx fitter."<<endl;
651	+
652	+	//=======================================================
653	+	// PVTX accessors
654	+	//=======================================================
655	+	//vertexCollection
656	+	static const reco::VertexCollection s_empty_vertexColl;
657	+	const reco::VertexCollection *vertexColl = &s_empty_vertexColl;
658	+	edm::Handle<reco::VertexCollection> vertexCollectionHandle;
659	+	iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle);
660	+	if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
661	+	vertexColl = vertexCollectionHandle.product();
662	+	} else {
663	+	cout << "vertexCollection cannot be found -> using empty collection of same type." <<endl;
664	+	}
665	+
666	+	//splitVertexCollection1
667	+	static const reco::VertexCollection s_empty_splitVertexColl1;
668	+	const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
669	+	edm::Handle<reco::VertexCollection> splitVertexCollection1Handle;
670	+	iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle);
671	+	if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
672	+	splitVertexColl1 = splitVertexCollection1Handle.product();
673	+	} else {
674	+	cout << "splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
675	+	}
676	+
677	+	//splitVertexCollection2
678	+	static const reco::VertexCollection s_empty_splitVertexColl2;
679	+	const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
680	+	edm::Handle<reco::VertexCollection> splitVertexCollection2Handle;
681	+	iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle);
682	+	if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
683	+	splitVertexColl2 = splitVertexCollection2Handle.product();
684	+	} else {
685	+	cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
686	+	}
687	+
688	+	if(verbose_) cout<<"End accessing the track and vertex collections"<<endl;
689	+
690	+	//=======================================================
691	+	// MC simvtx accessor
692	+	//=======================================================
693		if (!realData_) {
694	<	Handle<SimVertexContainer> simVtxs;
694	>	edm::Handle<SimVertexContainer> simVtxs;
695		iEvent.getByLabel( simG4_, simVtxs);
696
697	<	Handle<SimTrackContainer> simTrks;
697	>	edm::Handle<SimTrackContainer> simTrks;
698		iEvent.getByLabel( simG4_, simTrks);
699		}
700
701	<	if(verbose_) printRecVtxs(recVtxs);
702	<
701	>	//=======================================================
702	>	// GET PDT
703	>	//=======================================================
704		try{
705		iSetup.getData(pdt);
706		}catch(const Exception&){
707		cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
708		}
709	+
710	+	//=======================================================
711	+	// GET pixelVertices
712	+	//=======================================================
713	+	static const reco::VertexCollection s_empty_pixelVertexColl;
714	+	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
715	+	edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
716	+	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
717	+	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
718	+	pixelVertexColl = pixelVertexCollectionHandle.product();
719	+	} else {
720	+	cout << "pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
721	+	}
722	+
723	+	//=======================================================
724	+	// Fill pixelVertices related histograms
725	+	//=======================================================
726	+	nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
727	+	if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
728	+	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
729	+	fillTrackHistoInPV(pixelVertexColl, 4, false, true);
730	+	h["dzErr_pxlpvtx"]->Fill( pixelVertexColl->begin()->zError());
731	+	// Get the dZ error of the tracks in the leading pixelVertexColl
732	+	for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
733	+	t!= (pixelVertexColl->begin())->tracks_end(); t++) {
734	+
735	+	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
736	+	//        h["trkPtPV"]->Fill(0.);
737	+	}
738	+	else
739	+	h["trkdzErr_pxlpvtx"]->Fill((**t).dzError());
740	+	}
741	+	// Fill the track->dz() in the PV difference from first pixelpvtx
742	+	for(reco::Vertex::trackRef_iterator t = (vertexColl->begin())->tracks_begin();
743	+	t!= (vertexColl->begin())->tracks_end(); t++) {
744	+	// illegal charge
745	+	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
746	+	//        h["trkPtPV"]->Fill(0.);
747	+	}
748	+	else {
749	+	if(pixelVertexColl->size()>0) {
750	+	h["trkdz_pxlpvtxdz"]->Fill((**t).dz() -  pixelVertexColl->begin()->z());
751	+	h["trkdz_pxlpvtxdz_pxlpvtxdzerr"]->Fill(fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
752	+	h["trkdz_pxlpvtxdz_trkdzerr"]->Fill(fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
753	+	h["trkdzErr_pvtx"]->Fill((**t).dzError());
754	+	}
755	+	}
756	+	}
757	+	}
758	+
759	+	//=======================================================
760	+	// Fill number of reconstructed vertices
761	+	//=======================================================
762	+
763	+	if(verbose_)  {
764	+	cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
765	+	printRecVtxs(vertexCollectionHandle);
766	+	cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
767	+	printRecVtxs(splitVertexCollection1Handle);
768	+	cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
769	+	printRecVtxs(splitVertexCollection2Handle);
770	+	cout<<"Start filling pvtx parameters."<<endl;
771	+	}
772	+
773	+	nrecPV_ = int (vertexColl->size());
774	+	nrecPV1_spl_ = int (splitVertexColl1->size());
775	+	nrecPV2_spl_ = int (splitVertexColl2->size());
776	+
777	+	h["nrecPV"]->Fill(nrecPV_);
778	+	h["nrecPV1_spl"]->Fill(nrecPV1_spl_);
779	+	h["nrecPV2_spl"]->Fill(nrecPV2_spl_);
780	+	h["nrecPVDiff"]->Fill(double(nrecPV1_spl_)-double(nrecPV2_spl_));
781	+
782	+	//=================================================================
783	+	// Fill track parameter ntuple/hist for tracks used in recoVertices
784	+	//=================================================================
785	+
786	+	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
787	+	if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
788	+	fillTrackHistoInPV(vertexColl, 0, true, true);
789	+
790	+	if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
791	+	fillTrackHistoInPV(splitVertexColl1, 1, false, true);
792	+
793	+	if(splitVertexColl1->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
794	+	fillTrackHistoInPV(splitVertexColl2, 2, false, true);
795	+
796	+
797	+	//=======================================================
798	+	// Compare offlinePrimaryVertices with pixelVertices
799	+	//=======================================================
800	+	if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
801	+	&& (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {
802	+	h["nrecPV_minus_nrecPxlPV"]->Fill( double (nrecPV_ - nrecPV_pxlpvtx_));
803	+	// difference in reconstructed position of the leading pvtx
804	+	//cout<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;
805	+	//cout<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
806	+	h["recxPV_minus_recxPxlPV"]->Fill (recx_[0] - recx_pxlpvtx_[0]);
807	+	h["recyPV_minus_recyPxlPV"]->Fill (recy_[0] - recy_pxlpvtx_[0]);
808	+	h["reczPV_minus_reczPxlPV"]->Fill (recz_[0] - recz_pxlpvtx_[0]);
809	+	}
810	+
811
812	+
813	+	//==========================================================
814	+	// Fill secondary/primary min separations for multi-vertices
815	+	//==========================================================
816	+
817	+	if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {
818	+
819	+	double min_xsep = 9999.0;
820	+	double min_ysep = 9999.0;
821	+	double min_zsep = 9999.0;
822	+	double min_xsep_sign = 9999.0;
823	+	double min_ysep_sign = 9999.0;
824	+	double min_zsep_sign = 9999.0;
825	+	double min_ntrksep = 9999.0;
826	+	double min_sumpt2sep = 9999999.0;
827	+
828	+	if(verbose_) cout<<"leading pvtx z = "<< vertexColl->begin()->z() <<endl;
829	+
830	+	// Looping through the secondary vertices to find the mininum separation to the primary
831	+	for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
832	+	v!=vertexColl->end(); ++v) {
833	+	if(v->isValid() && ! v->isFake()) {
834	+	// xsep
835	+	if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
836	+	min_xsep = fabs(v->x()- vertexColl->begin()->x());
837	+	// ysep
838	+	if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
839	+	min_ysep = fabs(v->y()- vertexColl->begin()->y());
840	+	// zsep
841	+	if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
842	+	min_zsep = fabs(v->z()- vertexColl->begin()->z());
843	+	// xsep_sign
844	+	double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
845	+	if(xsep_sign < min_xsep_sign)
846	+	min_xsep_sign = xsep_sign;
847	+	// ysep_sign
848	+	double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
849	+	if(ysep_sign < min_ysep_sign)
850	+	min_ysep_sign = ysep_sign;
851	+	// zsep_sign
852	+	double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
853	+	if(zsep_sign < min_zsep_sign)
854	+	min_zsep_sign = zsep_sign;
855	+	// ntrksep
856	+	if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
857	+	min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
858	+	// sumpt2sep
859	+	if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
860	+	min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
861	+	}
862	+	}
863	+	h["min_xsep"]->Fill(min_xsep);
864	+	h["min_ysep"]->Fill(min_ysep);
865	+	h["min_zsep"]->Fill(min_zsep);
866	+	h["min_xsep_sign"]->Fill(min_xsep_sign);
867	+	h["min_ysep_sign"]->Fill(min_ysep_sign);
868	+	h["min_zsep_sign"]->Fill(min_zsep_sign);
869	+	h["min_ntrksep"]->Fill(min_ntrksep);
870	+	h["min_sumpt2sep"]->Fill(min_sumpt2sep);
871	+
872	+	min_zsep_ = min_zsep;
873	+	min_ntrksep_ = min_ntrksep;
874	+	min_sumpt2sep_ = min_sumpt2sep;
875	+
876	+
877	+	} // End of  if(vertexColl->size()>1) {
878	+
879	+	if(verbose_)
880	+	cout<<"End filling pvtx parameters."<<endl;
881	+
882	+	//=======================================================
883	+	//           PrimaryVertex Matching
884	+	// 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
885	+	// 2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
886	+	// Assume the first match is the primary vertex,
887	+	// since vertexColl are sorted in decreasing order of sum(pT)
888	+	//=======================================================
889	+
890	+	if(verbose_)   cout<<"PVStudy::analyze(): matching pvtxs "<<endl;
891	+	reco::VertexCollection s_empty_matchedVertexColl1;
892	+	reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
893	+	matchedVertexColl1->reserve(splitVertexColl1->size());
894	+	reco::VertexCollection s_empty_matchedVertexColl2;
895	+	reco::VertexCollection *matchedVertexColl2 = &s_empty_matchedVertexColl2;
896	+	matchedVertexColl2->reserve(splitVertexColl2->size());
897	+
898	+	bool stopmatching = false;
899	+	for (size_t ivtx = 0; ivtx<splitVertexCollection1Handle->size() && !stopmatching; ++ivtx) {
900	+	reco::VertexRef recvtx1(splitVertexCollection1Handle, ivtx);
901	+	if( !(recvtx1->isValid()) || recvtx1->isFake()) break;
902	+	for (size_t jvtx = ivtx; jvtx < splitVertexCollection2Handle->size(); ++jvtx) {
903	+	//------------------------------------------------------------------------
904	+	//== If only considering matching the first vertex of the splitVertexColl
905	+	//------------------------------------------------------------------------
906	+	if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
907	+	reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
908	+	if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
909	+	if(verbose_) {
910	+	cout<<"Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
911	+	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
912	+	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
913	+	}
914	+	if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
915	+	if(verbose_) {
916	+	cout<<"The two splitted vertices match in Z. "<<endl;
917	+	}
918	+	int nTrkPV1 = recvtx1->tracksSize();
919	+	int nTrkPV2 = recvtx2->tracksSize();
920	+	double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
921	+	h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
922	+	h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
923	+	if(fabs(ntrkreldiff)<ntrkdiffcut_) {
924	+	if(verbose_) {
925	+	cout<<"(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
926	+	}
927	+	matchedVertexColl1->push_back(*recvtx1);
928	+	matchedVertexColl2->push_back(*recvtx2);
929	+	stopmatching = true; // stop the matching when the first match is found!
930	+	break;
931	+	}
932	+	else
933	+	cout<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
934	+	}
935	+	else {
936	+	cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
937	+	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
938	+	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
939	+	}
940	+	}
941	+	}
942	+	if(verbose_)     cout<<"PVStudy::analyze(): End matching pvtxs"<<endl;
943	+
944	+	//=======================================================
945	+	//  Analyze matchedVertexColls
946	+	//=======================================================
947	+	if(verbose_) {
948	+	cout<<"Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
949	+	}
950	+
951	+	// Compare the reconstructed vertex position and calculate resolution/pulls
952	+	if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
953	+	//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
954	+	fillTrackHistoInPV(matchedVertexColl1, 1, true, false);
955	+	fillTrackHistoInPV(matchedVertexColl2, 2, true, false);
956	+
957	+	reco::VertexCollection::const_iterator v1;
958	+	reco::VertexCollection::const_iterator v2;
959	+	nrecPV_twovtx_ = 0;
960	+	for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
961	+	v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
962	+	++v1, ++v2) {
963	+	fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
964	+	fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
965	+	fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
966	+	ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
967	+	ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
968	+	ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
969	+	fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
970	+
971	+	h["deltax"]->Fill( fres_[0] );
972	+	h["deltay"]->Fill( fres_[1] );
973	+	h["deltaz"]->Fill( fres_[2] );
974	+	h["pullx"]->Fill( fres_[0]/ferror_[0]);
975	+	h["pully"]->Fill( fres_[1]/ferror_[1]);
976	+	h["pullz"]->Fill( fres_[2]/ferror_[2]);
977	+	h["errPVx"]->Fill( ferror_[0] );
978	+	h["errPVy"]->Fill( ferror_[1] );
979	+	h["errPVz"]->Fill( ferror_[2] );
980	+	pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
981	+	error(ferror_[0], ferror_[1], ferror_[2]),
982	+	fntrk_));
983	+	// Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
984	+	fillHisto(res(fres_[0], fres_[1], fres_[2]),
985	+	error(ferror_[0], ferror_[1], ferror_[2]),
986	+	fntrk_,0);
987	+	// Fill the ntuple variables
988	+	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
989	+	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
990	+	sumptsq1_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v1);
991	+	sumptsq2_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v2);
992	+	nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
993	+	deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
994	+	deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
995	+	deltaz_twovtx_[nrecPV_twovtx_] = fres_[2];
996	+	errx_twovtx_[nrecPV_twovtx_] = ferror_[0];
997	+	erry_twovtx_[nrecPV_twovtx_] = ferror_[1];
998	+	errz_twovtx_[nrecPV_twovtx_] = ferror_[2];
999	+	pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
1000	+	pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
1001	+	pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];
1002	+	nrecPV_twovtx_++;
1003	+	} // End of analyzing res/pull
1004	+	} // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
1005	+	else
1006	+	cout<<"WARNING: Cannot find matching pvtxs"<<endl;
1007	+
1008	+	if(verbose_)     cout<<"End analyzing the matchedVertexColl"<<endl;
1009	+
1010	+	//=======================================================
1011	+	// Fill simulated vertices
1012	+	//=======================================================
1013		if (!realData_) {
1014		bool MC=false;
1015		Handle<HepMCProduct> evtMC;
#	Line 400 | Line 1025 | PVStudy::analyze(const edm::Event& iEven
1025		}
1026		MC=true;
1027		}
1028	<
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	<
1035	<	int nsimtrk=0;
411	<	int nrectrk=0;
412	<	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
413	<	vsim!=simpv.end(); vsim++){
414	<	nsimtrk+=vsim->nGenTrk;
415	<	// look for a matching reconstructed vertex
416	<	vsim->recVtx=NULL;
417	<	for(reco::VertexCollection::const_iterator vrec=recVtxs->begin();
418	<	vrec!=recVtxs->end(); ++vrec){
419	<	nrectrk=vrec->tracksSize();
420	<	if(verbose_){
421	<	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
422	<	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
423	<	}
424	<	if ( matchVertex(*vsim,*vrec) ){
425	<	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
426	<	if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
427	<	|| (!vsim->recVtx) )
428	<	{
429	<	vsim->recVtx=&(*vrec);
430	<	}
431	<	}
432	<	}
433	<	if (vsim->recVtx){
434	<	if(verbose_) {
435	<	cout <<"primary matched " << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
436	<	}
437	<	fres[0] = vsim->recVtx->x()-vsim->x;
438	<	fres[1] = vsim->recVtx->y()-vsim->y;
439	<	fres[2] = vsim->recVtx->z()-vsim->z;
440	<	ferror[0] = vsim->recVtx->xError();
441	<	ferror[1] = vsim->recVtx->yError();
442	<	ferror[2] = vsim->recVtx->zError();
443	<	fntrk = nrectrk;
444	<
445	<	h["deltax"]->Fill( fres[0] );
446	<	h["deltay"]->Fill( fres[1] );
447	<	h["deltaz"]->Fill( fres[2] );
448	<	h["pullx"]->Fill( fres[0]/ferror[0] );
449	<	h["pully"]->Fill( fres[1]/ferror[1] );
450	<	h["pullz"]->Fill( fres[2]/ferror[2] );
451	<	h["errPVx"]->Fill( ferror[0] );
452	<	h["errPVy"]->Fill( ferror[1] );
453	<	h["errPVz"]->Fill( ferror[2] );
454	<	pvinfo.push_back(PVStudy::PVInfo(res(fres[0], fres[1], fres[2]),
455	<	error(ferror[0], ferror[1], ferror[2]),
456	<	nrectrk));
457	<	// Fill histo according to its track multiplicity
458	<	fillHisto(res(fres[0], fres[1], fres[2]),
459	<	error(ferror[0], ferror[1], ferror[2]),
460	<	nrectrk);
1034	>	h["nsimPV"]->Fill(simpv.size());
1035	>	nsimPV_ = simpv.size();
1036
1037	<	ftree_->Fill();
1038	<	}
1039	<	else{  // no rec vertex found for this simvertex
1040	<	if(verbose_) {
1041	<	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1042	<	}
1043	<	}
1044	<	}
1045	<	}
1046	<	}// With MC
1047	<
1048	<	// Loop RecVtxs and find matched SimVtxs
1049	<	for(reco::VertexCollection::const_iterator v=recVtxs->begin();
1050	<	v!=recVtxs->end(); ++v){
1051	<	try {
1052	<	for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
1053	<	t!=v->tracks_end(); t++) {
1054	<	// illegal charge
1055	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1056	<	//      h["trkPtPV"]->Fill(0.);
1057	<	}
483	<	else {
484	<	h["trkPtPV"]->Fill((**t).pt());
485	<	}
486	<	}
487	<	}
488	<	catch (...) {
489	<	// exception thrown when trying to use linked track
490	<	//          h["trkPtPV"]->Fill(0.);
491	<	}
492	<
493	<	h["nTrkPV"]->Fill(v->tracksSize());
494	<
495	<	}
496	<
497	<	h["nTrk"]->Fill(tracks->size());
498	<
499	<	for(TrackCollection::const_iterator itTrack = tracks->begin();
500	<	itTrack != tracks->end();
501	<	++itTrack) {
502	<	int charge = 0;
503	<	charge = itTrack->charge();
504	<	h["trkPt"]->Fill(itTrack->pt());
1037	>	int isimpv = 0;
1038	>	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
1039	>	vsim!=simpv.end(); vsim++, isimpv++){
1040	>	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
1041	>	simx_[isimpv] = vsim->x;
1042	>	simy_[isimpv] = vsim->y;
1043	>	simz_[isimpv] = vsim->y;
1044	>	simptsq_[isimpv] = vsim->ptsq;
1045	>
1046	>	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
1047	>	fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
1048	>	fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
1049	>	fillMCHisto(vsim, isimpv, vertexColl, 3);
1050	>	}
1051	>	} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
1052	>	} // End of   if (!realData_) {
1053	>
1054	>	// Finally fill the ftree_
1055	>	if(saventuple_) {
1056	>	ftree_->Fill();
1057	>	pvtxtree_->Fill();
1058		}
506	–
1059		}
1060
509	–
1061		// ------------ method called once each job just before starting event loop  ------------
1062		void
1063		PVStudy::beginJob()
#	Line 516 | Line 1067 | PVStudy::beginJob()
1067		// ------------ method called once each job just after ending the event loop  ------------
1068		void
1069		PVStudy::endJob() {
1070	<	double sqt2 = sqrt(2.);
1071	<	if (!realData_) {
1072	<	if (verbose_) cout << "Analyzing PV info" << endl;
1073	<	for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1074	<	if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
1075	<	PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk);
1076	<	if (analyze_) {
1077	<	//      if ( ipv.res_.x() > 0 ) h2["resx_nTrk"]->Fill(ntrk,ipv.res_.x()/sqt2, 1.);
1078	<	//      if ( ipv.res_.y() > 0 ) h2["resy_nTrk"]->Fill(ntrk,ipv.res_.y()/sqt2, 1.);
1079	<	//      if ( ipv.res_.z() > 0 ) h2["resz_nTrk"]->Fill(ntrk,ipv.res_.z()/sqt2, 1.);
1080	<	if ( ipv.res_.x() > 0 ) h2["resx_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
1081	<	if ( ipv.res_.y() > 0 ) h2["resy_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
1082	<	if ( ipv.res_.z() > 0 ) h2["resz_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
1083	<	if ( ipv.pull_.x() > 0 ) h2["pullx_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
533	<	if ( ipv.pull_.y() > 0 ) h2["pully_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
534	<	if ( ipv.pull_.z() > 0 ) h2["pullz_nTrk"]->Fill(ntrk,ipv.pull_.z(), 1.);
535	<	}
1070	>	if (verbose_) cout << "Analyzing PV info" << endl;
1071	>	// Looping through the datamodes available
1072	>	for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
1073	>	string suffix;
1074	>	if(verbose_)  cout<<"datamode = "<< *it<<endl;
1075	>	switch(*it) {
1076	>	case 0: suffix = "";
1077	>	break;
1078	>	case 1: suffix = "_spl1_mct";
1079	>	break;
1080	>	case 2: suffix = "_spl2_mct";
1081	>	break;
1082	>	case 3: suffix = "_mct";
1083	>	break;
1084		}
1085	+	if(analyze_) {
1086	+	for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1087	+	PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
1088	+	if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
1089	+	if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
1090	+	if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
1091	+	if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
1092	+	if ( ipv.pull_.y() > 0 ) h2["pully"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
1093	+	if ( ipv.pull_.z() > 0 ) h2["pullz"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.z(), 1.);
1094	+	}
1095	+	}
1096	+	suffix.clear();
1097	+	}
1098	+	if(saventuple_) {
1099	+	file_->cd();
1100	+	ftree_->Write();
1101	+	pvtxtree_->Write();
1102	+	file_->Close();
1103		}
538	–	file_->cd();
539	–	ftree_->Write();
540	–	file_->Close();
1104		}
1105
1106	+
1107	+
1108	+	// Match a recovertex with a simvertex
1109	+	// ? Should the cut be parameter dependent?
1110	+	// cut can be within n sigma of the vertex.zerror()
1111		bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
1112		const reco::Vertex & vrec)
1113		{
1114		return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1115		}
1116
1117	+	// Match two reconstructed vertices
1118	+	bool PVStudy::matchVertex( const reco::VertexRef & vrec1,
1119	+	const reco::VertexRef & vrec2,
1120	+	int zsigncut)
1121	+	{
1122	+	double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
1123	+	if(verbose_)
1124	+	cout<<"The matching criteria in z is "<<cut<<endl;
1125	+	return (fabs(vrec1->z()-vrec2->z())<cut);
1126	+	}
1127	+
1128	+
1129		bool PVStudy::isResonance(const HepMC::GenParticle * p){
1130		double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
1131		if(verbose_) cout << "isResonance   " << p->pdg_id() << " " << ctau << endl;
#	Line 617 | Line 1197 | void PVStudy::printSimTrks(const Handle<
1197		}
1198		}
1199
1200	<	void PVStudy::fillHisto(res r, error er, int ntk){
1200	>	void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1201		stringstream ss;
1202		ss << ntk;
1203		string suffix = ss.str();
1204	+	if(datamode == 0 )  suffix = "_" + suffix;
1205	+	if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1206	+	if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1207	+	if(datamode == 3 )  suffix = "_mct_" + suffix;
1208	+
1209		if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
1210		// Fill histograms of res and pull of ntk
1211	<	h["resx_"+suffix]->Fill(r.x());
1212	<	h["resy_"+suffix]->Fill(r.y());
1213	<	h["resz_"+suffix]->Fill(r.z());
1214	<	h["pullx_"+suffix]->Fill(r.x()/er.x());
1215	<	h["pully_"+suffix]->Fill(r.y()/er.y());
1216	<	h["pullz_"+suffix]->Fill(r.z()/er.z());
1217	<	h["errPVx_"+suffix]->Fill(er.x());
1218	<	h["errPVy_"+suffix]->Fill(er.y());
1219	<	h["errPVz_"+suffix]->Fill(er.z());
1211	>	h["deltax"+suffix]->Fill(r.x());
1212	>	h["deltay"+suffix]->Fill(r.y());
1213	>	h["deltaz"+suffix]->Fill(r.z());
1214	>	h["pullx"+suffix]->Fill(r.x()/er.x());
1215	>	h["pully"+suffix]->Fill(r.y()/er.y());
1216	>	h["pullz"+suffix]->Fill(r.z()/er.z());
1217	>	h["errPVx"+suffix]->Fill(er.x());
1218	>	h["errPVy"+suffix]->Fill(er.y());
1219	>	h["errPVz"+suffix]->Fill(er.z());
1220		}
1221
1222	<	PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk) {
1222	>	PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1223		map<int,double> data;
1224		data.clear();
1225		double fpar[2] = {-999,-999};
1226		double cm2um = 10000;
1227		stringstream ss;
1228		ss << ntk;
1229	<	string suffix = ss.str();
1230	<
1229	>	string suffix = ss.str();
1230	>	if(datamode == 0 )  suffix = "_" + suffix;
1231	>	if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1232	>	if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1233	>	if(datamode == 3 )  suffix = "_mct_" + suffix;
1234	>
1235		if(analyze_) {
1236		for ( int i = 0; i < 6; ++i) {
1237		switch (i) {
1238		case 0:
1239	<	fit(h["resx_"+suffix], fpar);
1239	>	fit(h["deltax"+suffix], fpar);
1240		data[i] = fpar[0]*cm2um;
1241		break;
1242		case 1:
1243	<	fit(h["resy_"+suffix], fpar);
1243	>	fit(h["deltay"+suffix], fpar);
1244		data[i] = fpar[0]*cm2um;
1245		break;
1246		case 2:
1247	<	fit(h["resz_"+suffix], fpar);
1247	>	fit(h["deltaz"+suffix], fpar);
1248		data[i] = fpar[0]*cm2um;
1249		break;
1250		case 3:
1251	<	fit(h["pullx_"+suffix], fpar);
1251	>	fit(h["pullx"+suffix], fpar);
1252		data[i] = fpar[0];
1253		break;
1254		case 4:
1255	<	fit(h["pully_"+suffix], fpar);
1255	>	fit(h["pully"+suffix], fpar);
1256		data[i] = fpar[0];
1257		break;
1258		case 5:
1259	<	fit(h["pullz_"+suffix], fpar);
1259	>	fit(h["pullz"+suffix], fpar);
1260		data[i] = fpar[0];
1261		break;
1262		}
#	Line 733 | Line 1322 | void PVStudy::fit(TH1 *&hdist, double fi
1322		fitPars[0] = -999;
1323		}
1324
1325	<	//define this as a plug-in
1326	<	DEFINE_FWK_MODULE(PVStudy);
1325	>	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype) {
1326	>	string suffix;
1327	>	suffix = ""; // for unsplit tracks
1328	>	if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1329	>	if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1330	>
1331	>	h["nTrk"+suffix]->Fill(trackColl->size());
1332	>	for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1333	>	itTrack != trackColl->end();
1334	>	++itTrack) {
1335	>	h["trkPt"+suffix]->Fill(itTrack->pt());
1336	>	h["trkDxy"+suffix]->Fill(itTrack->dxy());
1337	>	h["trkDz"+suffix]->Fill(itTrack->dz());
1338	>	h["trkEta"+suffix]->Fill(itTrack->eta());
1339	>	h["trkPhi"+suffix]->Fill(itTrack->phi());
1340	>	}
1341	>	}
1342	>
1343	>	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
1344	>	string suffix;
1345	>	suffix = ""; // for unsplit tracks
1346	>	if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1347	>	if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1348	>	int ivtx = 0;
1349	>	for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1350	>	v!=vertexColl->end(); ++v, ++ivtx) {
1351	>	if(fillNtuple) {
1352	>	if(datatype == 0) {
1353	>	nTrkPV_[ivtx] = v->tracksSize();
1354	>	sumptsq_[ivtx] = sumPtSquared(*v);
1355	>	isValid_[ivtx] = v->isValid();
1356	>	isFake_[ivtx] = v->isFake();
1357	>	recx_[ivtx] = v->x();
1358	>	recy_[ivtx] = v->y();
1359	>	recz_[ivtx] = v->z();
1360	>	recx_err_[ivtx] = v->xError();
1361	>	recy_err_[ivtx] = v->yError();
1362	>	recz_err_[ivtx] = v->zError();
1363	>
1364	>
1365	>	}
1366	>	if(datatype == 1) {
1367	>	nTrkPV1_spl_[ivtx] = v->tracksSize();
1368	>	sumptsq1_spl_[ivtx] = sumPtSquared(*v);
1369	>	isValid1_spl_[ivtx] = v->isValid();
1370	>	isFake1_spl_[ivtx] = v->isFake();
1371	>	recx1_spl_[ivtx] = v->x();
1372	>	recy1_spl_[ivtx] = v->y();
1373	>	recz1_spl_[ivtx] = v->z();
1374	>	recx1_err_spl_[ivtx] = v->xError();
1375	>	recy1_err_spl_[ivtx] = v->yError();
1376	>	recz1_err_spl_[ivtx] = v->zError();
1377	>
1378	>	}
1379	>	if(datatype == 2) {
1380	>	nTrkPV2_spl_[ivtx] = v->tracksSize();
1381	>	sumptsq2_spl_[ivtx] = sumPtSquared(*v);
1382	>	isValid2_spl_[ivtx] = v->isValid();
1383	>	isFake2_spl_[ivtx] = v->isFake();
1384	>	recx2_spl_[ivtx] = v->x();
1385	>	recy2_spl_[ivtx] = v->y();
1386	>	recz2_spl_[ivtx] = v->z();
1387	>	recx2_err_spl_[ivtx] = v->xError();
1388	>	recy2_err_spl_[ivtx] = v->yError();
1389	>	recz2_err_spl_[ivtx] = v->zError();
1390	>	}
1391	>	if(datatype == 4) { // for pixelVertices
1392	>	nTrkPV_pxlpvtx_[ivtx] = v->tracksSize();
1393	>	sumptsq_pxlpvtx_[ivtx] = sumPtSquared(*v);
1394	>	isValid_pxlpvtx_[ivtx] = v->isValid();
1395	>	isFake_pxlpvtx_[ivtx] = v->isFake();
1396	>	recx_pxlpvtx_[ivtx] = v->x();
1397	>	recy_pxlpvtx_[ivtx] = v->y();
1398	>	recz_pxlpvtx_[ivtx] = v->z();
1399	>	recx_err_pxlpvtx_[ivtx] = v->xError();
1400	>	recy_err_pxlpvtx_[ivtx] = v->yError();
1401	>	recz_err_pxlpvtx_[ivtx] = v->zError();
1402	>	}
1403	>	}
1404	>	}
1405	>	// For histogram only fill the leading pvtx parameters
1406	>	if(fillHisto) {
1407	>	h["nTrkPV"+suffix]->Fill(vertexColl->begin()->tracksSize());
1408	>	try {
1409	>	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1410	>	t!=vertexColl->begin()->tracks_end(); t++) {
1411	>	// illegal charge
1412	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1413	>	//      h["trkPtPV"]->Fill(0.);
1414	>	}
1415	>	else {
1416	>	h["trkPtPV"+suffix]->Fill((**t).pt());
1417	>	h["trkDxyPV"+suffix]->Fill((**t).dxy());
1418	>	h["trkDzPV"+suffix]->Fill((**t).dz());
1419	>	h["trkEtaPV"+suffix]->Fill((**t).eta());
1420	>	h["trkPhiPV"+suffix]->Fill((**t).phi());
1421	>	}
1422	>	}
1423	>	}
1424	>	catch (...) {
1425	>	// exception thrown when trying to use linked track
1426	>	//          h["trkPtPV"]->Fill(0.);
1427	>	}
1428	>	}
1429	>
1430	>	}
1431	>
1432	>	void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
1433	>	{
1434	>	std::string suffix;
1435	>	// Set the vertexColl and histogram suffix according to datamode
1436	>	if (datatype == 1) {
1437	>	suffix = "_spl1_mct";
1438	>	nrecPV_spl1_mct_ = 0;
1439	>	}
1440	>	if (datatype == 2) {
1441	>	suffix = "_spl2_mct";
1442	>	nrecPV_spl2_mct_ = 0;
1443	>	}
1444	>	if (datatype == 3) {
1445	>	suffix = "_mct";
1446	>	nrecPV_mct_ = 0;
1447	>	}
1448	>
1449	>	//========================================================
1450	>	//  look for a matching reconstructed vertex in vertexColl
1451	>	//========================================================
1452	>
1453	>	for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1454	>	vrec!=vtxColl->end(); ++vrec){
1455	>	int nrectrk = vrec->tracksSize();
1456	>	vsim->recVtx=NULL;
1457	>
1458	>	if(verbose_){
1459	>	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1460	>	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1461	>	}
1462	>	if ( matchVertex(*vsim,*vrec) ) {
1463	>	vsim->recVtx=&(*vrec);
1464	>
1465	>	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1466	>	//if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
1467	>	//|| (!vsim->recVtx) )
1468	>	//vsim->recVtx=&(*vrec);
1469	>	//}
1470	>	if(verbose_)
1471	>	cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1472	>
1473	>	double fres_mct[3];
1474	>	double ferror_mct[3];
1475	>
1476	>	fres_mct[0] = vsim->recVtx->x()-vsim->x;
1477	>	fres_mct[1] = vsim->recVtx->y()-vsim->y;
1478	>	fres_mct[2] = vsim->recVtx->z()-vsim->z;
1479	>	ferror_mct[0] = vsim->recVtx->xError();
1480	>	ferror_mct[1] = vsim->recVtx->yError();
1481	>	ferror_mct[2] = vsim->recVtx->zError();
1482	>
1483	>	h["deltax"+suffix]->Fill( fres_mct[0] );
1484	>	h["deltay"+suffix]->Fill( fres_mct[1] );
1485	>	h["deltaz"+suffix]->Fill( fres_mct[2] );
1486	>	h["pullx"+suffix]->Fill( fres_mct[0]/ferror_mct[0] );
1487	>	h["pully"+suffix]->Fill( fres_mct[1]/ferror_mct[1] );
1488	>	h["pullz"+suffix]->Fill( fres_mct[2]/ferror_mct[2] );
1489	>	h["errPVx"+suffix]->Fill( ferror_mct[0] );
1490	>	h["errPVy"+suffix]->Fill( ferror_mct[1] );
1491	>	h["errPVz"+suffix]->Fill( ferror_mct[2] );
1492	>	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1493	>	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1494	>	nrectrk));
1495	>	// Fill histo according to its track multiplicity
1496	>	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1497	>	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1498	>	nrectrk, datatype);
1499	>
1500	>	if(saventuple_) {
1501	>	//Fill the values for variables in ftree_
1502	>	if(datatype == 1) {
1503	>	nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1504	>	deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1505	>	deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1506	>	deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1507	>	pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
1508	>	pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1509	>	pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1510	>	nrecPV_spl1_mct_++;
1511	>	}
1512	>
1513	>	if(datatype == 2) {
1514	>	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1515	>	deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1516	>	deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1517	>	deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1518	>	pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
1519	>	pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1520	>	pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1521	>	nrecPV_spl2_mct_++;
1522	>	}
1523	>	if(datatype == 3) {
1524	>	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1525	>	deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1526	>	deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1527	>	deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1528	>	pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
1529	>	pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1530	>	pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1531	>	nrecPV_mct_++;
1532	>	}
1533	>	} // End of  if(saventuple_) {
1534	>	} //  if ( matchVertex(*vsim,*vrec) ) {
1535	>	else {  // no rec vertex found for this simvertex
1536	>	if(verbose_)
1537	>	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1538	>	}
1539	>	}
1540	>	}
1541	>
1542	>
1543	>	void PVStudy::SetVarToZero() {
1544	>	//Greg's variables
1545	>	fntrk_ = 0;
1546	>	//pvtx position (x,y,z) residual and error
1547	>	for(int i = 0; i<3;i++)    {
1548	>	fres_[i] = 0;
1549	>	ferror_[i] = 0;
1550	>	}
1551	>
1552	>	//Yanyan's variables
1553	>	// Number of reconstructed vertices
1554	>	nrecPV_ = 0;
1555	>	nrecPV1_spl_ = 0;
1556	>	nrecPV2_spl_ = 0;
1557	>	nrecPV_twovtx_ = 0;
1558	>	nsimPV_ = 0;
1559	>
1560	>	nrecPV_mct_ = 0;
1561	>	nrecPV_spl1_mct_ = 0;
1562	>	nrecPV_spl2_mct_ = 0;
1563	>
1564	>	// Mininum separation between the secondary pvtxes and leading pvtx
1565	>	min_zsep_ = 9999.0;
1566	>	min_ntrksep_ = 9999.0;
1567	>	min_sumpt2sep_ = -9999.0;
1568	>
1569	>
1570	>	for (int i = 0; i < nMaxPVs_; i++) {
1571	>	// recoVertices with all tracks
1572	>	nTrkPV_[i] = 0; // Number of tracks in the pvtx
1573	>	sumptsq_[i] = 0;
1574	>	isValid_[i] = -1;
1575	>	isFake_[i] = -1;
1576	>	recx_[i] = 0;
1577	>	recy_[i] = 0;
1578	>	recz_[i] = 0;
1579	>	recx_err_[i] = 0;
1580	>	recy_err_[i] = 0;
1581	>	recz_err_[i] = 0;
1582	>
1583	>	// recoVertices with splitTrack1
1584	>	nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx
1585	>	sumptsq1_spl_[i] = 0;
1586	>	isValid1_spl_[i] = -1;
1587	>	isFake1_spl_[i] = -1;
1588	>	recx1_spl_[i] = 0;
1589	>	recy1_spl_[i] = 0;
1590	>	recz1_spl_[i] = 0;
1591	>	recx1_err_spl_[i] = 0;
1592	>	recy1_err_spl_[i] = 0;
1593	>	recz1_err_spl_[i] = 0;
1594	>
1595	>	// recoVertices with splitTrack2
1596	>	nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx
1597	>	sumptsq2_spl_[i] = 0;
1598	>	isValid2_spl_[i] = -1;
1599	>	isFake2_spl_[i] = -1;
1600	>	recx2_spl_[i] = 0;
1601	>	recy2_spl_[i] = 0;
1602	>	recz2_spl_[i] = 0;
1603	>	recx2_err_spl_[i] = 0;
1604	>	recy2_err_spl_[i] = 0;
1605	>	recz2_err_spl_[i] = 0;
1606	>
1607	>	//pixelVertices
1608	>	nTrkPV_pxlpvtx_[i] = 0; // Number of tracks in the pvtx
1609	>	sumptsq_pxlpvtx_[i] = 0;
1610	>	isValid_pxlpvtx_[i] = -1;
1611	>	isFake_pxlpvtx_[i] = -1;
1612	>	recx_pxlpvtx_[i] = 0;
1613	>	recy_pxlpvtx_[i] = 0;
1614	>	recz_pxlpvtx_[i] = 0;
1615	>	recx_err_pxlpvtx_[i] = 0;
1616	>	recy_err_pxlpvtx_[i] = 0;
1617	>	recz_err_pxlpvtx_[i] = 0;
1618	>
1619	>	// matched two-vertices
1620	>	nTrkPV1_spl_twovtx_[i] = 0;
1621	>	nTrkPV2_spl_twovtx_[i] = 0;
1622	>	nTrkPV_twovtx_[i] = 0;
1623	>	sumptsq1_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1624	>	sumptsq2_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1625	>	deltax_twovtx_[i] = 0;
1626	>	deltay_twovtx_[i] = 0;
1627	>	deltaz_twovtx_[i] = 0;
1628	>	errx_twovtx_[i] = 0;
1629	>	erry_twovtx_[i] = 0;
1630	>	errz_twovtx_[i] = 0;
1631	>	pullx_twovtx_[i] = 0;
1632	>	pully_twovtx_[i] = 0;
1633	>	pullz_twovtx_[i] = 0;
1634	>
1635	>	//simpv
1636	>	nsimTrkPV_[i] = 0;
1637	>	simx_[i] = 0;
1638	>	simy_[i] = 0;
1639	>	simz_[i] = 0;
1640	>	simptsq_[i] = 0;
1641	>
1642	>	// residual and pulls with mc truth required
1643	>	deltax_mct_[i] = 0;
1644	>	deltay_mct_[i] = 0;
1645	>	deltaz_mct_[i] = 0;
1646	>	pullx_mct_[i] = 0;
1647	>	pully_mct_[i] = 0;
1648	>	pullz_mct_[i] = 0;
1649	>
1650	>	deltax_spl1_mct_[i] = 0;
1651	>	deltay_spl1_mct_[i] = 0;
1652	>	deltaz_spl1_mct_[i] = 0;
1653	>	pullx_spl1_mct_[i] = 0;
1654	>	pully_spl1_mct_[i] = 0;
1655	>	pullz_spl1_mct_[i] = 0;
1656	>
1657	>	deltax_spl2_mct_[i] = 0;
1658	>	deltay_spl2_mct_[i] = 0;
1659	>	deltaz_spl2_mct_[i] = 0;
1660	>	pullx_spl2_mct_[i] = 0;
1661	>	pully_spl2_mct_[i] = 0;
1662	>	pullz_spl2_mct_[i] = 0;
1663	>	}
1664	>
1665	>	}
1666	>
1667	>	double PVStudy::sumPtSquared(const reco::Vertex & v)  {
1668	>	double sum = 0.;
1669	>	double pT;
1670	>	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1671	>	pT = (**it).pt();
1672	>	sum += pT*pT;
1673	>	}
1674	>	return sum;
1675	>	}
1676	>

Diff Legend

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