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.8 by jengbou, Fri Oct 16 00:11:30 2009 UTC vs.
Revision 1.15 by jengbou, Fri Dec 4 21:20:26 2009 UTC

#	Line 12 | Line 12 | Implementation:
12		*/
13		//
14		// Original Author:  "Geng-yuan Jeng/UC Riverside"
15	<	//                    "Yanyan Gao/Fermilab ygao@fnal.gov"
15	>	//                   "Yanyan Gao/Fermilab ygao@fnal.gov"
16		//         Created:  Thu Aug 20 11:55:40 CDT 2009
17		// $Id$
18		//
#	Line 25 | Line 25 | Implementation:
25		#include <vector>
26		#include <iostream>
27		#include <sstream>
28	+	#include <algorithm>
29
30		// user include files
31		#include "FWCore/Framework/interface/Frameworkfwd.h"
#	Line 71 | Line 72 | PVStudy::PVStudy(const edm::ParameterSet
72		splitTrackCollection2Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
73		vertexCollectionTag_       = iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection");
74		splitVertexCollection1Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
75	<	splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
75	>	splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
76	>	pixelVertexCollectionTag_  = iConfig.getParameter<edm::InputTag>("pixelVertexCollectionTag");
77		verbose_                   = iConfig.getUntrackedParameter<bool>("verbose",false);
78		realData_                  = iConfig.getUntrackedParameter<bool>("realData",false);
79		analyze_                   = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
80		saventuple_                = iConfig.getUntrackedParameter<bool>("saventuple",false);
81		outputfilename_            = iConfig.getUntrackedParameter<string>("OutputFileName");
82		ntrkdiffcut_               = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
83	+	zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
84		nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
85		nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
86	+	histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
87
88		//now do what ever initialization is needed
89		pvinfo.clear();
90
91	<
91	>
92		// Specify the data mode vector
93		if(realData_) datamode.push_back(0);
94		else {
#	Line 102 | Line 106 | PVStudy::PVStudy(const edm::ParameterSet
106		ftree_->Branch("residual",&fres_,"fres_[3]/D");
107		ftree_->Branch("error",&ferror_,"ferror_[3]/D");
108		ftree_->Branch("nTrkPV",&fntrk_,"fntrk_/I");
109	<	ftree_->Branch("nrecPV",&nrecPV_,"nrecPV_/I");
110	<	ftree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl_/I");
111	<	ftree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl_/I");
112	<
113	<	// If it is MC, get the resA, resB
114	<	if(!realData_) {
115	<	ftree_->Branch("residual_spl1_mct",&fres_spl1_mct_,"fres_spl1_mct_[3]/D");
116	<	ftree_->Branch("error_spl1_mct",&ferror_spl1_mct_,"ferror_spl1_mct_[3]/D");
117	<	ftree_->Branch("nTrkPV_spl1_mct",&fntrk_spl1_mct_,"fntrk_spl1_mct_/I");
118	<	ftree_->Branch("residual_spl2_mct",&fres_spl2_mct_,"fres_spl2_mct_[3]/D");
119	<	ftree_->Branch("error_spl2_mct",&ferror_spl2_mct_,"ferror_spl2_mct_[3]/D");
120	<	ftree_->Branch("nTrkPV_spl2_mct",&fntrk_spl2_mct_,"fntrk_spl2_mct_/I");
121	<	ftree_->Branch("residual_mct",&fres_mct_,"fres_mct_[3]/D");
122	<	ftree_->Branch("error_mct",&ferror_mct_,"ferror_mct_[3]/D");
123	<	ftree_->Branch("nTrkPV_mct",&fntrk_mct_,"fntrk_mct_/I");
109	>
110	>	// pvtxtree_ analyzing the pvtxs ootb
111	>	pvtxtree_ = new TTree("pvtxtree","pvtxtree");
112	>	//pvtx using all tracks
113	>	pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
114	>	pvtxtree_->Branch("nTrkPV",&nTrkPV_,"nTrkPV[nrecPV]/I");
115	>	pvtxtree_->Branch("sumptsq",&sumptsq_,"sumptsq[nrecPV]/D");
116	>	pvtxtree_->Branch("isValid",&isValid_,"isValid/I");
117	>	pvtxtree_->Branch("isFake",&isFake_,"isFake/I");
118	>	pvtxtree_->Branch("recx",&recx_,"recx[nrecPV]/D");
119	>	pvtxtree_->Branch("recy",&recy_,"recy[nrecPV]/D");
120	>	pvtxtree_->Branch("recz",&recz_,"recz[nrecPV]/D");
121	>	pvtxtree_->Branch("recx_err",&recx_err_,"recx_err[nrecPV]/D");
122	>	pvtxtree_->Branch("recy_err",&recy_err_,"recy_err[nrecPV]/D");
123	>	pvtxtree_->Branch("recz_err",&recz_err_,"recz_err[nrecPV]/D");
124	>
125	>	pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
126	>	pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
127	>	pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
128	>
129	>	//pvtx using splittracks1
130	>	pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
131	>	pvtxtree_->Branch("nTrkPV1_spl",&nTrkPV1_spl_,"nTrkPV1_spl[nrecPV1_spl]/I");
132	>	pvtxtree_->Branch("sumptsq1_spl",&sumptsq1_spl_,"sumptsq1_spl[nrecPV1_spl]/D");
133	>	pvtxtree_->Branch("isValid1_spl",&isValid1_spl_,"isValid1_spl/I");
134	>	pvtxtree_->Branch("isFake1_spl",&isFake1_spl_,"isFake1_spl/I");
135	>	pvtxtree_->Branch("recx1_spl",&recx1_spl_,"recx1_spl[nrecPV1_spl]/D");
136	>	pvtxtree_->Branch("recy1_spl",&recy1_spl_,"recy1_spl[nrecPV1_spl]/D");
137	>	pvtxtree_->Branch("recz1_spl",&recz1_spl_,"recz1_spl[nrecPV1_spl]/D");
138	>	pvtxtree_->Branch("recx1_err_spl",&recx1_err_spl_,"recx1_err_spl[nrecPV1_spl]/D");
139	>	pvtxtree_->Branch("recy1_err_spl",&recy1_err_spl_,"recy1_err_spl[nrecPV1_spl]/D");
140	>	pvtxtree_->Branch("recz1_err_spl",&recz1_err_spl_,"recz1_err_spl[nrecPV1_spl]/D");
141	>
142	>	//pvtx using splittracks2
143	>	pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
144	>	pvtxtree_->Branch("nTrkPV2_spl",&nTrkPV2_spl_,"nTrkPV2_spl[nrecPV2_spl]/I");
145	>	pvtxtree_->Branch("sumptsq2_spl",&sumptsq2_spl_,"sumptsq2_spl[nrecPV2_spl]/D");
146	>	pvtxtree_->Branch("isValid2_spl",&isValid2_spl_,"isValid2_spl/I");
147	>	pvtxtree_->Branch("isFake2_spl",&isFake2_spl_,"isFake2_spl/I");
148	>	pvtxtree_->Branch("recx2_spl",&recx2_spl_,"recx2_spl[nrecPV2_spl]/D");
149	>	pvtxtree_->Branch("recy2_spl",&recy2_spl_,"recy2_spl[nrecPV2_spl]/D");
150	>	pvtxtree_->Branch("recz2_spl",&recz2_spl_,"recz2_spl[nrecPV2_spl]/D");
151	>	pvtxtree_->Branch("recx2_err_spl",&recx2_err_spl_,"recx2_err_spl[nrecPV2_spl]/D");
152	>	pvtxtree_->Branch("recy2_err_spl",&recy2_err_spl_,"recy2_err_spl[nrecPV2_spl]/D");
153	>	pvtxtree_->Branch("recz2_err_spl",&recz2_err_spl_,"recz2_err_spl[nrecPV2_spl]/D");
154	>
155	>	//pixeVertices
156	>	pvtxtree_->Branch("nrecPV_pxlpvtx",&nrecPV_pxlpvtx_,"nrecPV_pxlpvtx/I");
157	>	pvtxtree_->Branch("nTrkPV_pxlpvtx",&nTrkPV_pxlpvtx_,"nTrkPV_pxlpvtx[nrecPV_pxlpvtx]/I");
158	>	pvtxtree_->Branch("sumptsq_pxlpvtx",&sumptsq_pxlpvtx_,"sumptsq_pxlpvtx[nrecPV_pxlpvtx]/D");
159	>	pvtxtree_->Branch("isValid_pxlpvtx",&isValid_pxlpvtx_,"isValid_pxlpvtx/I");
160	>	pvtxtree_->Branch("isFake_pxlpvtx",&isFake_pxlpvtx_,"isFake_pxlpvtx/I");
161	>	pvtxtree_->Branch("recx_pxlpvtx",&recx_pxlpvtx_,"recx_pxlpvtx[nrecPV_pxlpvtx]/D");
162	>	pvtxtree_->Branch("recy_pxlpvtx",&recy_pxlpvtx_,"recy_pxlpvtx[nrecPV_pxlpvtx]/D");
163	>	pvtxtree_->Branch("recz_pxlpvtx",&recz_pxlpvtx_,"recz_pxlpvtx[nrecPV_pxlpvtx]/D");
164	>	pvtxtree_->Branch("recx_err_pxlpvtx",&recx_err_pxlpvtx_,"recx_err_pxlpvtx[nrecPV_pxlpvtx]/D");
165	>	pvtxtree_->Branch("recy_err_pxlpvtx",&recy_err_pxlpvtx_,"recy_err_pxlpvtx[nrecPV_pxlpvtx]/D");
166	>	pvtxtree_->Branch("recz_err_pxlpvtx",&recz_err_pxlpvtx_,"recz_err_pxlpvtx[nrecPV_pxlpvtx]/D");
167	>
168	>	//Fill the variables in the twovtx pair (recvtx1, recvtx2)
169	>	pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
170	>	pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
171	>	pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
172	>	pvtxtree_->Branch("sumptsq1_spl_twovtx",&sumptsq1_spl_twovtx_,"sumptsq1_spl_twovtx[nrecPV_twovtx]/D");
173	>	pvtxtree_->Branch("sumptsq2_spl_twovtx",&sumptsq2_spl_twovtx_,"sumptsq2_spl_twovtx[nrecPV_twovtx]/D");
174	>	pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
175	>	pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
176	>	pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
177	>	pvtxtree_->Branch("deltaz_twovtx",&deltaz_twovtx_,"deltaz_twovtx[nrecPV_twovtx]/D");
178	>	pvtxtree_->Branch("errx_twovtx",&errx_twovtx_,"errx_twovtx[nrecPV_twovtx]/D");
179	>	pvtxtree_->Branch("erry_twovtx",&erry_twovtx_,"erry_twovtx[nrecPV_twovtx]/D");
180	>	pvtxtree_->Branch("errz_twovtx",&errz_twovtx_,"errz_twovtx[nrecPV_twovtx]/D");
181	>	pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
182	>	pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
183	>	pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
184	>
185	>	// MC variables
186	>	if(!realData_) {
187	>	pvtxtree_->Branch("nsimPV",&nsimPV_,"nsimPV/I");
188	>	pvtxtree_->Branch("nsimTrkPV",&nsimTrkPV_,"nsimTrkPV[nsimPV]/I");
189	>	pvtxtree_->Branch("simx",&simx_,"simx[nsimPV]/D");
190	>	pvtxtree_->Branch("simy",&simy_,"simy[nsimPV]/D");
191	>	pvtxtree_->Branch("simz",&simz_,"simz[nsimPV]/D");
192	>	pvtxtree_->Branch("simptsq",&simptsq_,"simptsq[nsimPV]/D");
193	>
194	>	// For pvtxs with all the tracks
195	>	pvtxtree_->Branch("nrecPV_mct",&nrecPV_mct_,"nrecPV_mct/I");
196	>	pvtxtree_->Branch("deltax_mct",&deltax_mct_,"deltax_mct[nrecPV_mct]/D");
197	>	pvtxtree_->Branch("deltay_mct",&deltay_mct_,"deltay_mct[nrecPV_mct]/D");
198	>	pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
199	>	pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
200	>	pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
201	>	pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
202	>	// For pvtxs with splittracks1
203	>	pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
204	>	pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
205	>	pvtxtree_->Branch("deltay_spl1_mct",&deltay_spl1_mct_,"deltay_spl1_mct[nrecPV_spl1_mct]/D");
206	>	pvtxtree_->Branch("deltaz_spl1_mct",&deltaz_spl1_mct_,"deltaz_spl1_mct[nrecPV_spl1_mct]/D");
207	>	pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
208	>	pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
209	>	pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
210	>	// For pvtxs with splittracks1
211	>	pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
212	>	pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
213	>	pvtxtree_->Branch("deltay_spl2_mct",&deltay_spl2_mct_,"deltay_spl2_mct[nrecPV_spl2_mct]/D");
214	>	pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
215	>	pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
216	>	pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
217	>	pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
218		}
219		}
220
123	–	edm::Service<TFileService> fs;
124	–	TFileDirectory subDir = fs->mkdir( "Summary" );
125	–	TFileDirectory subDir1 = fs->mkdir( "Others" );
126	–	//   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
127	–
221		setRootStyle();
222	<
223	<	//Book histograms for track and pvtx parameters for each splitted pvtx
224	<	//i=0 : x (as in unsplit track collection)
225	<	//i=1 : x1_spl_
226	<	//i=2 : x2_spl_
227	<
228	<	for(int i=0;i<3;i++)
229	<	{
230	<	string suffix;
231	<	if(i == 0) suffix = "";
232	<	else {
233	<	stringstream ss;
234	<	ss << i;
235	<	suffix =ss.str()+"_spl";
236	<	}
237	<	h["nTrk"+suffix]   = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
238	<	h["trkPt"+suffix]  = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
239	<	h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#Eta of rec tracks "+suffix),100,-3,3);
240	<	h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
241	<	h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-5,5);
242	<	h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
243	<
244	<	h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
245	<	h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
246	<	h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#Eta of rec tracks in PV"+suffix),100,-3,3);
247	<	h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
155	<	h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);
156	<	h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
157	<	h["nrecPV"+suffix]   = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);
158	<	suffix.clear();
222	>
223	>	// Create a root file for histograms
224	>	theFile = new TFile(histoFileName_.c_str(), "RECREATE");
225	>	// make diretories
226	>	theFile->mkdir("Summary");
227	>	theFile->cd();
228	>	theFile->mkdir("Others");
229	>	theFile->cd();
230	>	if (analyze_) {
231	>	theFile->mkdir("Results");
232	>	theFile->cd();
233	>	}
234	>
235	>	// Book Histograms:
236	>	h_pvtrk   = new PVHistograms();
237	>	h_pixvtx  = new PVHistograms();
238	>	h_misc    = new PVHistograms();
239	>	h_summary = new PVHistograms();
240	>	h_others  = new PVHistograms();
241	>
242	>	for(int i=0;i<3;i++) {
243	>	if(i == 0) h_pvtrk->Init("pvTrk");
244	>	else {
245	>	stringstream ss;
246	>	ss << i;
247	>	h_pvtrk->Init("pvTrk",ss.str(),"spl");
248		}
249	<	h["nrecPVDiff"]     = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
250	<	h["nTrkPVDiff"]     = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
251	<	h["nTrkPVRelDiff"]  = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
252	<
249	>	}
250	>	h_pixvtx->Init("pixVtx");
251	>	h_misc->Init("misc");
252	>
253	>	// Book MC only plots
254	>	if (!realData_) {
255	>	h_gen = new PVHistograms();
256	>	h_gen->Init("generator");
257	>	}
258	>	if (analyze_) h_ana = new PVHistograms();
259	>
260		//Book histograms sensitive to data/mc
261		for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
262		string suffix;
263	<	if(verbose_)  cout<<"datamode = "<< *it<<endl;
263	>	edm::LogInfo("Debug")<<"datamode = "<< *it<<endl;
264		switch(*it) {
265		case 0: suffix = "";
266		break;
#	Line 175 | Line 271 | PVStudy::PVStudy(const edm::ParameterSet
271		case 3: suffix = "_mct";
272		break;
273		}
274	<	h["deltax"+suffix] = subDir.make<TH1D>(TString("deltax"+suffix), TString("x-residual pvtx"+suffix),800,-0.04,0.04);
179	<	h["deltay"+suffix] = subDir.make<TH1D>(TString("deltay"+suffix), TString("y-residual pvtx"+suffix),800,-0.04,0.04);
180	<	h["deltaz"+suffix] = subDir.make<TH1D>(TString("deltaz"+suffix), TString("z-residual pvtx"+suffix),800,-0.04,0.04);
181	<	h["pullx"+suffix]  = subDir.make<TH1D>(TString("pullx"+suffix), TString("x-pull pvtx"+suffix),800,-10,10);
182	<	h["pully"+suffix]  = subDir.make<TH1D>(TString("pully"+suffix), TString("y-pull pvtx"+suffix),800,-10,10);
183	<	h["pullz"+suffix]  = subDir.make<TH1D>(TString("pullz"+suffix), TString("z-pull pvtx"+suffix),800,-10,10);
184	<	h["errPVx"+suffix] = subDir.make<TH1D>(TString("errPVx"+suffix), TString("X"+suffix+" vertex error"),100,0.,0.02);
185	<	h["errPVy"+suffix] = subDir.make<TH1D>(TString("errPVy"+suffix), TString("Y"+suffix+" vertex error"),100,0.,0.02);
186	<	h["errPVz"+suffix] = subDir.make<TH1D>(TString("errPVz"+suffix), TString("Z"+suffix+" vertex error"),100,0.,0.02);
274	>	h_summary->Init("summary",suffix);
275
276		// Book residual, error and pull histograms for each nTrk bin
277		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
278		stringstream ss;
279		ss << ntrk;
280	<	string suffix2 = suffix+"_"+ss.str();
281	<	h["deltax"+suffix2]  = subDir1.make<TH1D>(TString("deltax"+suffix2), TString("x-residual of pvtx"+suffix),100,-0.02,0.02);
282	<	h["deltax"+suffix2]->GetXaxis()->SetTitle("cm");
195	<	h["deltay"+suffix2]  = subDir1.make<TH1D>(TString("deltay"+suffix2), TString("y-residual of pvtx"+suffix),100,-0.02,0.02);
196	<	h["deltay"+suffix2]->GetXaxis()->SetTitle("cm");
197	<	h["deltaz"+suffix2]  = subDir1.make<TH1D>(TString("deltaz"+suffix2), TString("z-residual of pvtx"+suffix),100,-0.02,0.02);
198	<	h["deltaz"+suffix2]->GetXaxis()->SetTitle("cm");
199	<	h["pullx"+suffix2] = subDir1.make<TH1D>(TString("pullx"+suffix2), TString("x-pull of pvtx"+suffix),100,-10.,10.);
200	<	h["pully"+suffix2] = subDir1.make<TH1D>(TString("pully"+suffix2), TString("y-pull of pvtx"+suffix),100,-10.,10.);
201	<	h["pullz"+suffix2] = subDir1.make<TH1D>(TString("pullz"+suffix2), TString("z-pull of pvtx"+suffix),100,-10.,10.);
202	<	h["errPVx"+suffix2] = subDir1.make<TH1D>(TString("errPVx"+suffix2), TString("X"+suffix+" vertex error"),100,0.,0.02);
203	<	h["errPVy"+suffix2] = subDir1.make<TH1D>(TString("errPVy"+suffix2), TString("Y"+suffix+" vertex error"),100,0.,0.02);
204	<	h["errPVz"+suffix2] = subDir1.make<TH1D>(TString("errPVz"+suffix2), TString("Z"+suffix+" vertex error"),100,0.,0.02);
205	<	suffix2.clear();
206	<	} // End of   for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
207	<
280	>	h_others->Init("others",suffix,ss.str());
281	>	}
282	>
283		// Book residual and pull histograms only when analyzeOntheFly is enabled
284	<	if(analyze_) {
210	<	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);
211	<	h2["resx"+suffix+"_nTrk"]->SetMarkerStyle(21);
212	<	h2["resx"+suffix+"_nTrk"]->SetMarkerColor(4);
213	<	h2["resx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
214	<	h2["resx"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
215	<	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);
216	<	h2["resy"+suffix+"_nTrk"]->SetMarkerStyle(21);
217	<	h2["resy"+suffix+"_nTrk"]->SetMarkerColor(4);
218	<	h2["resy"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
219	<	h2["resy"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
220	<	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);
221	<	h2["resz"+suffix+"_nTrk"]->SetMarkerStyle(21);
222	<	h2["resz"+suffix+"_nTrk"]->SetMarkerColor(4);
223	<	h2["resz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
224	<	h2["resz"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
225	<	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.);
226	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerStyle(21);
227	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerColor(4);
228	<	h2["pullx"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
229	<	h2["pullx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
230	<	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.);
231	<	h2["pully"+suffix+"_nTrk"]->SetMarkerStyle(21);
232	<	h2["pully"+suffix+"_nTrk"]->SetMarkerColor(4);
233	<	h2["pully"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
234	<	h2["pully"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
235	<	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.);
236	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerStyle(21);
237	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerColor(4);
238	<	h2["pullz"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
239	<	h2["pullz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
240	<	} // End of   if(analyze_) {
284	>	if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
285		suffix.clear();
286	<	} // End of Book histograms sensitive to data/mc
287	<	// Book MC only plots
244	<	if (!realData_) {
245	<	h["genPart_T"]      = subDir.make<TH1D>("genPart_T","t component of gen particles",300,-0.5,0.5);
246	<	h["genPart_T"]->GetXaxis()->SetTitle("t (nanosecond)");
247	<	h["genPart_cT"]     = subDir.make<TH1D>("genPart_cT","ct component of gen particles",300,-150.,150.);
248	<	h["genPart_cT"]->GetXaxis()->SetTitle("ct (mm)");
249	<	h["nsimPV"]         = subDir.make<TH1D>("nsimPV","Num of sim PV",51,-0.5,50.5);
250	<	}
286	>	} // End of Book histograms sensitive to data/mc
287	>
288		}
289
290		PVStudy::~PVStudy()
291		{
292	<
293	<	// do anything here that needs to be done at desctruction time
294	<	// (e.g. close files, deallocate resources etc.)
292	>	theFile->cd();
293	>	theFile->cd("Summary");
294	>	h_pvtrk->Save();
295	>	h_pixvtx->Save();
296	>	h_misc->Save();
297	>	h_summary->Save();
298	>	if (!realData_)
299	>	h_gen->Save();
300	>	if (analyze_) {
301	>	theFile->cd();
302	>	theFile->cd("Results");
303	>	h_ana->Save();
304	>	}
305	>	theFile->cd();
306	>	theFile->cd("Others");
307	>	h_others->Save();
308
309	+	theFile->Close();
310		}
311
312		void PVStudy::setRootStyle() {
#	Line 295 | Line 346 | std::vector<PVStudy::simPrimaryVertex> P
346		std::vector<PVStudy::simPrimaryVertex> simpv;
347		const HepMC::GenEvent* evt=evtMC->GetEvent();
348		if (evt) {
349	<	if(verbose_) {
350	<	cout << "process id " << evt->signal_process_id()<<endl;
351	<	cout << "signal process vertex " << ( evt->signal_process_vertex() ?
352	<	evt->signal_process_vertex()->barcode() : 0 )   <<endl;
302	<	cout << "number of vertices " << evt->vertices_size() << endl;
303	<	}
349	>	edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
350	>	edm::LogInfo("SimPVs") << "[getSimPVs] signal process vertex " << ( evt->signal_process_vertex() ?
351	>	evt->signal_process_vertex()->barcode() : 0 ) <<endl;
352	>	edm::LogInfo("SimPVs") << "[getSimPVs] number of vertices " << evt->vertices_size() << endl;
353
354		int idx=0; int npv=0;
355		for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
#	Line 311 | Line 360 | std::vector<PVStudy::simPrimaryVertex> P
360		// t component of PV:
361		for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
362		mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
363	<	//        std::cout << "Status = " << (*mother)->status() << std::endl;
363	>	//        edm::LogInfo("SimPVs") << "Status = " << (*mother)->status() << endl;
364		HepMC::GenVertex * mv=(*mother)->production_vertex();
365		if( ((*mother)->status() == 3) && (!mv)) {
366	<	//        std::cout << "npv= " << npv << std::endl;
366	>	//        edm::LogInfo("SimPVs") << "npv= " << npv << endl;
367		if (npv == 0) {
368	<	h["genPart_cT"]->Fill(pos.t()); // mm
369	<	h["genPart_T"]->Fill(pos.t()/299.792458); // ns
368	>	h_gen->Fill1d("genPart_cT", pos.t()); // mm
369	>	h_gen->Fill1d("genPart_T", pos.t()/299.792458); // ns
370		}
371		npv++;
372		}
#	Line 325 | Line 374 | std::vector<PVStudy::simPrimaryVertex> P
374		//       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
375
376		bool hasMotherVertex=false;
377	<	if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl;
377	>	if (verbose_) cout << "[getSimPVs] mothers of vertex[" << ++idx << "]: " << endl;
378		for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
379		mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
380		HepMC::GenVertex * mv=(*mother)->production_vertex();
381	<	//        std::cout << "Status = " << (*mother)->status() << std::endl;
381	>	// if (verbose_) cout << "Status = " << (*mother)->status() << endl;
382		if (mv) {
383		hasMotherVertex=true;
384		if(!verbose_) break; //if verbose_, print all particles of gen vertices
#	Line 356 | Line 405 | std::vector<PVStudy::simPrimaryVertex> P
405
406		if(!vp){
407		// this is a new vertex
408	<	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
408	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
409		simpv.push_back(sv);
410		vp=&simpv.back();
411		}else{
412	<	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
412	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
413		}
414		vp->genVertex.push_back((*vitr)->barcode());
415		// collect final state descendants
#	Line 384 | Line 433 | std::vector<PVStudy::simPrimaryVertex> P
433		}
434		}
435		}
436	<	}
437	<	}
436	>	}//loop MC vertices daughters
437	>	}//loop MC vertices
438		}
439		return simpv;
440		}
#	Line 418 | Line 467 | std::vector<PVStudy::simPrimaryVertex> P
467		ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
468		}
469		}
470	<	primary =  ( gv->position().t()==0);
470	>	//primary =  ( gv->position().t()==0 );
471	>	primary = true;
472	>	h_gen->Fill1d("genPart_cT", gv->position().t()); // mm
473	>	h_gen->Fill1d("genPart_T", gv->position().t()/299.792458); // ns
474	>
475		//resonance= ( gp->mother() && isResonance(gp->mother()));  // in CLHEP/HepMC days
476		// no more mother pointer in the improved HepMC GenParticle
477		resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
#	Line 467 | Line 520 | PVStudy::analyze(const edm::Event& iEven
520		using namespace edm;
521		using namespace reco;
522
523	<	if (verbose_)
471	<	cout<<"PVStudy::analyze():"<<endl;
523	>	edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
524		//=======================================================
525		// Initialize Root-tuple variables if needed
526		//=======================================================
527		//if(saventuple_)
528	<	SetVarToZero();
528	>	SetVarToZero();
529
530		//=======================================================
531		// Track accessors
#	Line 486 | Line 538 | PVStudy::analyze(const edm::Event& iEven
538		if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
539		trackColl = trackCollectionHandle.product();
540		} else {
541	<	cout << "trackCollection cannot be found -> using empty collection of same type." <<endl;
541	>	edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
542		}
543	+
544		//splitTrackCollection1
545		static const reco::TrackCollection s_empty_splitTrackColl1;
546		const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
#	Line 496 | Line 549 | PVStudy::analyze(const edm::Event& iEven
549		if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
550		splitTrackColl1 = splitTrackCollection1Handle.product();
551		} else {
552	<	cout << "splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
552	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
553		}
554	+
555		//splitTrackCollection2
556		static const reco::TrackCollection s_empty_splitTrackColl2;
557		const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
#	Line 506 | Line 560 | PVStudy::analyze(const edm::Event& iEven
560		if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
561		splitTrackColl2 = splitTrackCollection2Handle.product();
562		} else {
563	<	cout << "splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
563	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
564		}
565	<
565	>
566	>
567	>	//=======================================================
568	>	// Fill trackparameters of the input tracks to pvtx fitter
569	>	//=======================================================
570	>	edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
571	>	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype)
572	>	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
573	>	fillTrackHisto(trackColl, 0);
574	>	fillTrackHisto(splitTrackColl1, 1);
575	>	fillTrackHisto(splitTrackColl2, 2);
576	>	edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
577	>
578		//=======================================================
579		// PVTX accessors
580		//=======================================================
#	Line 520 | Line 586 | PVStudy::analyze(const edm::Event& iEven
586		if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
587		vertexColl = vertexCollectionHandle.product();
588		} else {
589	<	cout << "vertexCollection cannot be found -> using empty collection of same type." <<endl;
589	>	edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
590		}
591	+
592		//splitVertexCollection1
593		static const reco::VertexCollection s_empty_splitVertexColl1;
594		const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
#	Line 530 | Line 597 | PVStudy::analyze(const edm::Event& iEven
597		if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
598		splitVertexColl1 = splitVertexCollection1Handle.product();
599		} else {
600	<	cout << "splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
601	<	}
600	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
601	>	}
602	>
603		//splitVertexCollection2
604		static const reco::VertexCollection s_empty_splitVertexColl2;
605		const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
#	Line 540 | Line 608 | PVStudy::analyze(const edm::Event& iEven
608		if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
609		splitVertexColl2 = splitVertexCollection2Handle.product();
610		} else {
611	<	cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
611	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
612		}
613	<
614	<	if(verbose_) cout<<"Done accessing the track and vertex collections"<<endl;
613	>
614	>	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track and vertex collections"<<endl;
615
616		//=======================================================
617		// MC simvtx accessor
#	Line 562 | Line 630 | PVStudy::analyze(const edm::Event& iEven
630		try{
631		iSetup.getData(pdt);
632		}catch(const Exception&){
633	<	cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
633	>	edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
634		}
635	<
635	>
636		//=======================================================
637	<	// Fill trackparameters of the input tracks to pvtx fitter
637	>	// GET pixelVertices
638		//=======================================================
639	<
640	<	// GeneralTracks
641	<	h["nTrk"]->Fill(trackColl->size());
642	<	for(TrackCollection::const_iterator itTrack = trackColl->begin();
643	<	itTrack != trackColl->end();
644	<	++itTrack) {
645	<	h["trkPt"]->Fill(itTrack->pt());
646	<	h["trkDxy"]->Fill(itTrack->dxy());
579	<	h["trkDz"]->Fill(itTrack->dz());
580	<	h["trkEta"]->Fill(itTrack->eta());
581	<	h["trkPhi"]->Fill(itTrack->phi());
639	>	static const reco::VertexCollection s_empty_pixelVertexColl;
640	>	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
641	>	edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
642	>	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
643	>	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
644	>	pixelVertexColl = pixelVertexCollectionHandle.product();
645	>	} else {
646	>	edm::LogInfo("Debug") << "[PVStudy] pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
647		}
648	<
649	<	//SplittedTracks1
650	<	h["nTrk1_spl"]->Fill(splitTrackColl1->size());
651	<	for(TrackCollection::const_iterator itTrack = splitTrackColl1->begin();
652	<	itTrack != splitTrackColl1->end();
653	<	++itTrack) {
654	<	h["trkPt1_spl"]->Fill(itTrack->pt());
655	<	h["trkDxy1_spl"]->Fill(itTrack->dxy());
656	<	h["trkDz1_spl"]->Fill(itTrack->dz());
657	<	h["trkEta1_spl"]->Fill(itTrack->eta());
658	<	h["trkPhi1_spl"]->Fill(itTrack->phi());
659	<	}
660	<	//SplittedTracks2
661	<	h["nTrk2_spl"]->Fill(splitTrackColl2->size());
662	<	for(TrackCollection::const_iterator itTrack = splitTrackColl2->begin();
663	<	itTrack != splitTrackColl2->end();
664	<	++itTrack) {
665	<	h["trkPt2_spl"]->Fill(itTrack->pt());
666	<	h["trkDxy2_spl"]->Fill(itTrack->dxy());
667	<	h["trkDz2_spl"]->Fill(itTrack->dz());
668	<	h["trkEta2_spl"]->Fill(itTrack->eta());
669	<	h["trkPhi2_spl"]->Fill(itTrack->phi());
648	>
649	>	//=======================================================
650	>	// Fill pixelVertices related histograms
651	>	//=======================================================
652	>	nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
653	>	if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
654	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
655	>	fillTrackHistoInPV(pixelVertexColl, 4, false, true);
656	>	h_pixvtx->Fill1d("dzErr_pxlpvtx", pixelVertexColl->begin()->zError());
657	>	// Get the dZ error of the tracks in the leading pixelVertexColl
658	>	for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
659	>	t!= (pixelVertexColl->begin())->tracks_end(); t++) {
660	>
661	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
662	>	//        h_pvtrk->Fill1d("trkPtPV", 0.);
663	>	}
664	>	else
665	>	h_pixvtx->Fill1d("trkdzErr_pxlpvtx", (**t).dzError());
666	>	}
667	>	// Fill the track->dz() in the PV difference from first pixelpvtx
668	>	for(reco::Vertex::trackRef_iterator t = (vertexColl->begin())->tracks_begin();
669	>	t!= (vertexColl->begin())->tracks_end(); t++) {
670	>	// illegal charge
671	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
672	>	//        h_pvtrk->Fill1d("trkPtPV", 0.);
673	>	}
674	>	else {
675	>	if(pixelVertexColl->size()>0) {
676	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz", (**t).dz() -  pixelVertexColl->begin()->z());
677	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_pxlpvtxdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
678	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_trkdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
679	>	h_pixvtx->Fill1d("trkdzErr_pvtx", (**t).dzError());
680	>	}
681	>	}
682	>	}
683		}
606	–
607	–	if(verbose_)
608	–	cout<<"Done filling track parameters of the input tracks to pvtx fitter."<<endl;
684
685		//=======================================================
686	<	// Fill reconstructed vertices
686	>	// Fill number of reconstructed vertices
687		//=======================================================
688	<	if(verbose_)  {
689	<	cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
688	>
689	>	edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
690	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
691	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
692	>	edm::LogInfo("Debug")<<"[PVStudy] Start filling pvtx parameters."<<endl;
693	>	if (verbose_) {
694		printRecVtxs(vertexCollectionHandle);
616	–	cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
695		printRecVtxs(splitVertexCollection1Handle);
618	–	cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
696		printRecVtxs(splitVertexCollection2Handle);
697		}
698	+	nrecPV_ = int (vertexColl->size());
699	+	nrecPV1_spl_ = int (splitVertexColl1->size());
700	+	nrecPV2_spl_ = int (splitVertexColl2->size());
701	+
702	+	h_pvtrk->Fill1d("nrecPV", nrecPV_);
703	+	h_pvtrk->Fill1d("nrecPV1_spl", nrecPV1_spl_);
704	+	h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
705	+	h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
706
707	<	//Fill the tree variables
708	<	nrecPV_ = vertexColl->size();
709	<	nrecPV1_spl_ = splitVertexColl1->size();
710	<	nrecPV2_spl_ = splitVertexColl2->size();
711	<
712	<	h["nrecPV"]->Fill(nrecPV_);
713	<	h["nrecPV1_spl"]->Fill(nrecPV1_spl_);
629	<	h["nrecPV2_spl"]->Fill(nrecPV2_spl_);
630	<	h["nrecPVDiff"]->Fill(double(nrecPV1_spl_)-double(nrecPV2_spl_));
707	>	//=================================================================
708	>	// Fill track parameter ntuple/hist for tracks used in recoVertices
709	>	//=================================================================
710	>
711	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
712	>	if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
713	>	fillTrackHistoInPV(vertexColl, 0, true, true);
714
715	<	// Fill the track paramters for the vertexColl
716	<	for(reco::VertexCollection::const_iterator v=vertexColl->begin();
717	<	v!=vertexColl->end(); ++v){
718	<	h["nTrkPV"]->Fill(v->tracksSize());
719	<	try {
720	<	for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
721	<	t!=v->tracks_end(); t++) {
722	<	// illegal charge
723	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
724	<	//      h["trkPtPV"]->Fill(0.);
725	<	}
726	<	else {
727	<	h["trkPtPV"]->Fill((**t).pt());
728	<	h["trkDxyPV"]->Fill((**t).dxy());
729	<	h["trkDzPV"]->Fill((**t).dz());
730	<	h["trkEtaPV"]->Fill((**t).eta());
731	<	h["trkPhiPV"]->Fill((**t).phi());
732	<	}
733	<	}
651	<	}
652	<	catch (...) {
653	<	// exception thrown when trying to use linked track
654	<	//          h["trkPtPV"]->Fill(0.);
655	<	}
715	>	if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
716	>	fillTrackHistoInPV(splitVertexColl1, 1, false, true);
717	>
718	>	if(splitVertexColl1->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
719	>	fillTrackHistoInPV(splitVertexColl2, 2, false, true);
720	>
721	>
722	>	//=======================================================
723	>	// Compare offlinePrimaryVertices with pixelVertices
724	>	//=======================================================
725	>	if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
726	>	&& (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {
727	>	h_pixvtx->Fill1d("nrecPV_minus_nrecPxlPV", double (nrecPV_ - nrecPV_pxlpvtx_));
728	>	// difference in reconstructed position of the leading pvtx
729	>	//edm::LogInfo("Debug")<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;
730	>	//edm::LogInfo("Debug")<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
731	>	h_pixvtx->Fill1d("recxPV_minus_recxPxlPV", recx_[0] - recx_pxlpvtx_[0]);
732	>	h_pixvtx->Fill1d("recyPV_minus_recyPxlPV", recy_[0] - recy_pxlpvtx_[0]);
733	>	h_pixvtx->Fill1d("reczPV_minus_reczPxlPV", recz_[0] - recz_pxlpvtx_[0]);
734		}
735
736	+
737	+
738	+	//==========================================================
739	+	// Fill secondary/primary min separations for multi-vertices
740	+	//==========================================================
741
742	<	// == ignore multi-vertices (for now)
743	<	//take single pvt collection and compare
744	<	if(splitVertexColl1->size() == 1 && splitVertexColl2->size() == 1 ) {
745	<	const reco::Vertex recvtx1 = *(splitVertexColl1->begin());
746	<	const reco::Vertex recvtx2 = *(splitVertexColl2->begin());
747	<	if (matchVertex(recvtx1, recvtx2))  {
748	<	int nTrkPV1 = recvtx1.tracksSize();
749	<	int nTrkPV2 = recvtx2.tracksSize();
750	<	h["nTrkPV1_spl"]->Fill(nTrkPV1);
751	<	h["nTrkPV2_spl"]->Fill(nTrkPV2);
752	<	double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
753	<	h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
754	<	h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
755	<
756	<	// Fill the track paramters for the splitVertexColl1
757	<	try {
758	<	for(reco::Vertex::trackRef_iterator t = recvtx1.tracks_begin();
759	<	t!=recvtx1.tracks_end(); t++) {
760	<	// illegal charge
761	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
762	<	//    h["trkPtPV"]->Fill(0.);
763	<	}
764	<	else {
765	<	h["trkPtPV1_spl"]->Fill((**t).pt());
766	<	h["trkDxyPV1_spl"]->Fill((**t).dxy());
767	<	h["trkDzPV1_spl"]->Fill((**t).dz());
768	<	h["trkEtaPV1_spl"]->Fill((**t).eta());
769	<	h["trkPhiPV1_spl"]->Fill((**t).phi());
770	<	}
771	<	}
772	<	}
773	<	catch (...) {
774	<	// exception thrown when trying to use linked track
775	<	//        h["trkPtPV"]->Fill(0.);
742	>	if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {
743	>
744	>	double min_xsep = 9999.0;
745	>	double min_ysep = 9999.0;
746	>	double min_zsep = 9999.0;
747	>	double min_xsep_sign = 9999.0;
748	>	double min_ysep_sign = 9999.0;
749	>	double min_zsep_sign = 9999.0;
750	>	double min_ntrksep = 9999.0;
751	>	double min_sumpt2sep = 9999999.0;
752	>
753	>	edm::LogInfo("Debug")<<"[PVStudy] leading pvtx z = "<< vertexColl->begin()->z() <<endl;
754	>
755	>	// Looping through the secondary vertices to find the mininum separation to the primary
756	>	for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
757	>	v!=vertexColl->end(); ++v) {
758	>	if(v->isValid() && ! v->isFake()) {
759	>	// xsep
760	>	if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
761	>	min_xsep = fabs(v->x()- vertexColl->begin()->x());
762	>	// ysep
763	>	if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
764	>	min_ysep = fabs(v->y()- vertexColl->begin()->y());
765	>	// zsep
766	>	if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
767	>	min_zsep = fabs(v->z()- vertexColl->begin()->z());
768	>	// xsep_sign
769	>	double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
770	>	if(xsep_sign < min_xsep_sign)
771	>	min_xsep_sign = xsep_sign;
772	>	// ysep_sign
773	>	double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
774	>	if(ysep_sign < min_ysep_sign)
775	>	min_ysep_sign = ysep_sign;
776	>	// zsep_sign
777	>	double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
778	>	if(zsep_sign < min_zsep_sign)
779	>	min_zsep_sign = zsep_sign;
780	>	// ntrksep
781	>	if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
782	>	min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
783	>	// sumpt2sep
784	>	if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
785	>	min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
786		}
787	<
788	<	// Fill the track paramters for the splitVertexColl2
789	<	try {
790	<	for(reco::Vertex::trackRef_iterator t = recvtx2.tracks_begin();
791	<	t!=recvtx2.tracks_end(); t++) {
792	<	// illegal charge
793	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
794	<	//    h["trkPtPV"]->Fill(0.);
795	<	}
796	<	else {
797	<	h["trkPtPV2_spl"]->Fill((**t).pt());
798	<	h["trkDxyPV2_spl"]->Fill((**t).dxy());
799	<	h["trkDzPV2_spl"]->Fill((**t).dz());
800	<	h["trkEtaPV2_spl"]->Fill((**t).eta());
801	<	h["trkPhiPV2_spl"]->Fill((**t).phi());
802	<	}
787	>	}
788	>	h_misc->Fill1d("min_xsep", min_xsep);
789	>	h_misc->Fill1d("min_ysep", min_ysep);
790	>	h_misc->Fill1d("min_zsep", min_zsep);
791	>	h_misc->Fill1d("min_xsep_sign", min_xsep_sign);
792	>	h_misc->Fill1d("min_ysep_sign", min_ysep_sign);
793	>	h_misc->Fill1d("min_zsep_sign", min_zsep_sign);
794	>	h_misc->Fill1d("min_ntrksep", min_ntrksep);
795	>	h_misc->Fill1d("min_sumpt2sep", min_sumpt2sep);
796	>
797	>	min_zsep_ = min_zsep;
798	>	min_ntrksep_ = min_ntrksep;
799	>	min_sumpt2sep_ = min_sumpt2sep;
800	>
801	>
802	>	} // End of  if(vertexColl->size()>1) {
803	>
804	>	edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
805	>
806	>	//=======================================================
807	>	//           PrimaryVertex Matching
808	>	// 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
809	>	// 2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
810	>	// Assume the first match is the primary vertex,
811	>	// since vertexColl are sorted in decreasing order of sum(pT)
812	>	//=======================================================
813	>
814	>	edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
815	>	reco::VertexCollection s_empty_matchedVertexColl1;
816	>	reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
817	>	matchedVertexColl1->reserve(splitVertexColl1->size());
818	>	reco::VertexCollection s_empty_matchedVertexColl2;
819	>	reco::VertexCollection *matchedVertexColl2 = &s_empty_matchedVertexColl2;
820	>	matchedVertexColl2->reserve(splitVertexColl2->size());
821	>
822	>	bool stopmatching = false;
823	>	for (size_t ivtx = 0; ivtx<splitVertexCollection1Handle->size() && !stopmatching; ++ivtx) {
824	>	reco::VertexRef recvtx1(splitVertexCollection1Handle, ivtx);
825	>	if( !(recvtx1->isValid()) || recvtx1->isFake()) break;
826	>	for (size_t jvtx = ivtx; jvtx < splitVertexCollection2Handle->size(); ++jvtx) {
827	>	//------------------------------------------------------------------------
828	>	//== If only considering matching the first vertex of the splitVertexColl
829	>	//------------------------------------------------------------------------
830	>	if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
831	>	reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
832	>	if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
833	>	edm::LogInfo("Debug")<<"[PVStudy] Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
834	>	edm::LogInfo("Debug")<<"[PVStudy] recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
835	>	edm::LogInfo("Debug")<<"[PVStudy] recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
836	>
837	>	if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
838	>	edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
839	>
840	>	int nTrkPV1 = recvtx1->tracksSize();
841	>	int nTrkPV2 = recvtx2->tracksSize();
842	>	double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
843	>	h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
844	>	h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
845	>	if(fabs(ntrkreldiff)<ntrkdiffcut_) {
846	>	edm::LogInfo("Debug")<<"[PVStudy] (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
847	>
848	>	matchedVertexColl1->push_back(*recvtx1);
849	>	matchedVertexColl2->push_back(*recvtx2);
850	>	stopmatching = true; // stop the matching when the first match is found!
851	>	break;
852		}
853	+	else
854	+	edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
855		}
856	<	catch (...) {
857	<	// exception thrown when trying to use linked track
858	<	//        h["trkPtPV"]->Fill(0.);
856	>	else {
857	>	edm::LogInfo("Debug")<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
858	>	edm::LogInfo("Debug")<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
859	>	edm::LogInfo("Debug")<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
860		}
716	–
717	–	if(verbose_)
718	–	cout<<"Done filling track parameters in reconstruced vertices."<<endl;
719	–
720	–	if(fabs(ntrkreldiff)<ntrkdiffcut_) {
721	–	fres_[0] = (recvtx1.x()-recvtx2.x())/sqrt(2.0);
722	–	fres_[1] = (recvtx1.y()-recvtx2.y())/sqrt(2.0);
723	–	fres_[2] = (recvtx1.z()-recvtx2.z())/sqrt(2.0);
724	–	ferror_[0] = sqrt(pow(recvtx1.xError(),2)+pow(recvtx2.xError(),2))/sqrt(2);
725	–	ferror_[1] = sqrt(pow(recvtx1.yError(),2)+pow(recvtx2.yError(),2))/sqrt(2);
726	–	ferror_[2] = sqrt(pow(recvtx1.zError(),2)+pow(recvtx2.zError(),2))/sqrt(2);
727	–	fntrk_ =  (nTrkPV1+nTrkPV2)/2;
728	–
729	–	h["deltax"]->Fill( fres_[0] );
730	–	h["deltay"]->Fill( fres_[1] );
731	–	h["deltaz"]->Fill( fres_[2] );
732	–	h["pullx"]->Fill( fres_[0]/ferror_[0]);
733	–	h["pully"]->Fill( fres_[1]/ferror_[1]);
734	–	h["pullz"]->Fill( fres_[2]/ferror_[2]);
735	–	h["errPVx"]->Fill( ferror_[0] );
736	–	h["errPVy"]->Fill( ferror_[1] );
737	–	h["errPVz"]->Fill( ferror_[2] );
738	–	pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
739	–	error(ferror_[0], ferror_[1], ferror_[2]),
740	–	fntrk_));
741	–	// Fill histo according to its track multiplicity, set datamode = 0
742	–	fillHisto(res(fres_[0], fres_[1], fres_[2]),
743	–	error(ferror_[0], ferror_[1], ferror_[2]),
744	–	fntrk_,0);
745	–
746	–	if (saventuple_) ftree_->Fill();
747	–	} // End of  if(fabs(ntrkreldiff)<ntrkdiffcut_) {
748	–	} // End of  if (matchVertex(recvtx1, recvtx2))  {
749	–	else   {
750	–	cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
751	–	cout<<"recvtx1.z() = " << recvtx1.z() <<"; recvtx2.z() = " <<  recvtx2.z() <<endl;
861		}
862	<	} // End of fill histograms for two-vertices Method
863	<
864	<	if(verbose_)
865	<	cout<<"Done filling res/pull with two-vertices method"<<endl;
862	>	}
863	>	edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
864	>
865	>	//=======================================================
866	>	//  Analyze matchedVertexColls
867	>	//=======================================================
868	>	edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
869	>
870	>	// Compare the reconstructed vertex position and calculate resolution/pulls
871	>	if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
872	>	//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
873	>	fillTrackHistoInPV(matchedVertexColl1, 1, true, false);
874	>	fillTrackHistoInPV(matchedVertexColl2, 2, true, false);
875	>
876	>	reco::VertexCollection::const_iterator v1;
877	>	reco::VertexCollection::const_iterator v2;
878	>	nrecPV_twovtx_ = 0;
879	>	for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
880	>	v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
881	>	++v1, ++v2) {
882	>	fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
883	>	fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
884	>	fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
885	>	ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
886	>	ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
887	>	ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
888	>	fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
889	>
890	>	h_summary->Fill1d("deltax", fres_[0] );
891	>	h_summary->Fill1d("deltay", fres_[1] );
892	>	h_summary->Fill1d("deltaz", fres_[2] );
893	>	h_summary->Fill1d("pullx", fres_[0]/ferror_[0]);
894	>	h_summary->Fill1d("pully", fres_[1]/ferror_[1]);
895	>	h_summary->Fill1d("pullz", fres_[2]/ferror_[2]);
896	>	h_summary->Fill1d("errPVx", ferror_[0] );
897	>	h_summary->Fill1d("errPVy", ferror_[1] );
898	>	h_summary->Fill1d("errPVz", ferror_[2] );
899	>	pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
900	>	error(ferror_[0], ferror_[1], ferror_[2]),
901	>	fntrk_));
902	>	// Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
903	>	fillHisto(res(fres_[0], fres_[1], fres_[2]),
904	>	error(ferror_[0], ferror_[1], ferror_[2]),
905	>	fntrk_,0);
906	>	// Fill the ntuple variables
907	>	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
908	>	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
909	>	sumptsq1_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v1);
910	>	sumptsq2_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v2);
911	>	nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
912	>	deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
913	>	deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
914	>	deltaz_twovtx_[nrecPV_twovtx_] = fres_[2];
915	>	errx_twovtx_[nrecPV_twovtx_] = ferror_[0];
916	>	erry_twovtx_[nrecPV_twovtx_] = ferror_[1];
917	>	errz_twovtx_[nrecPV_twovtx_] = ferror_[2];
918	>	pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
919	>	pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
920	>	pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];
921	>	nrecPV_twovtx_++;
922	>	} // End of analyzing res/pull
923	>	} // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
924	>	else
925	>	edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
926	>
927	>	edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
928
929		//=======================================================
930		// Fill simulated vertices
#	Line 764 | Line 935 | PVStudy::analyze(const edm::Event& iEven
935		iEvent.getByLabel("generator",evtMC);
936		if (!evtMC.isValid()) {
937		MC=false;
938	<	if(verbose_){
768	<	cout << "no HepMCProduct found"<< endl;
769	<	}
938	>	edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
939		} else {
940	<	if(verbose_){
772	<	cout << "generator HepMCProduct found"<< endl;
773	<	}
940	>	edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
941		MC=true;
942		}
943
#	Line 779 | Line 946 | PVStudy::analyze(const edm::Event& iEven
946		std::vector<simPrimaryVertex> simpv;
947		simpv=getSimPVs(evtMC,"");
948		//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
949	<	h["nsimPV"]->Fill(simpv.size());
950	<	int nsimtrk=0;
951	<	int nrectrk=0;
952	<
949	>	h_gen->Fill1d("nsimPV", simpv.size());
950	>	nsimPV_ = simpv.size();
951	>
952	>	int isimpv = 0;
953		for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
954	<	vsim!=simpv.end(); vsim++){
955	<	nsimtrk+=vsim->nGenTrk;
956	<	// Loop over the datamodes, for MC, fill only 1, 2, 3
957	<	for (int i=1;i<=3;i++) {
958	<	static const reco::VertexCollection s_empty_vtxColl;
959	<	const reco::VertexCollection *vtxColl = &s_empty_vtxColl;
960	<	std::string suffix;
961	<	// Set the vertexColl and histogram suffix according to datamode
962	<	if (i == 1) {
963	<	vtxColl = splitVertexColl1;
964	<	suffix = "_spl1_mct";
798	<	}
799	<	if (i == 2) {
800	<	vtxColl =  splitVertexColl2;
801	<	suffix = "_spl2_mct";
802	<	}
803	<	if (i == 3) {
804	<	vtxColl =  vertexColl;
805	<	suffix = "_mct";
806	<	}
807	<	//========================================================
808	<	//  look for a matching reconstructed vertex in vertexColl
809	<	//========================================================
810	<	vsim->recVtx=NULL;
811	<	for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
812	<	vrec!=vtxColl->end(); ++vrec){
813	<	nrectrk=vrec->tracksSize();
814	<	if(verbose_){
815	<	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
816	<	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
817	<	}
818	<	if ( matchVertex(*vsim,*vrec) ){
819	<	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
820	<	if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
821	<	|| (!vsim->recVtx) )
822	<	vsim->recVtx=&(*vrec);
823	<	}// End of finding the matched reco_vertex
824	<
825	<	if (vsim->recVtx){
826	<	if(verbose_) {
827	<	cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
828	<	}
829	<	double fres_mct[3];
830	<	double ferror_mct[3];
831	<	int fntrk_mct;
832	<
833	<	fres_mct[0] = vsim->recVtx->x()-vsim->x;
834	<	fres_mct[1] = vsim->recVtx->y()-vsim->y;
835	<	fres_mct[2] = vsim->recVtx->z()-vsim->z;
836	<	ferror_mct[0] = vsim->recVtx->xError();
837	<	ferror_mct[1] = vsim->recVtx->yError();
838	<	ferror_mct[2] = vsim->recVtx->zError();
839	<	fntrk_mct = nrectrk;
840	<
841	<	//Fill the values for variables in ftree_
842	<	if(i == 1) {
843	<	for(int j = 0;j<3;j++)
844	<	fres_spl1_mct_[j] = fres_mct[j];
845	<	fntrk_spl1_mct_ =  fntrk_mct;
846	<	}
847	<	if(i == 2) {
848	<	for(int j = 0;j<3;j++)
849	<	fres_spl2_mct_[j] = fres_mct[j];
850	<	fntrk_spl2_mct_ =  fntrk_mct;
851	<	}
852	<	if(i == 3) {
853	<	for(int j=0;j<3;j++)
854	<	fres_mct_[j] = fres_mct[j];
855	<	fntrk_mct_ =  fntrk_mct;
856	<	}
857	<
858	<	h["deltax"+suffix]->Fill( fres_mct[0] );
859	<	h["deltay"+suffix]->Fill( fres_mct[1] );
860	<	h["deltaz"+suffix]->Fill( fres_mct[2] );
861	<	h["pullx"+suffix]->Fill( fres_mct[0]/ferror_mct[0] );
862	<	h["pully"+suffix]->Fill( fres_mct[1]/ferror_mct[1] );
863	<	h["pullz"+suffix]->Fill( fres_mct[2]/ferror_mct[2] );
864	<	h["errPVx"+suffix]->Fill( ferror_mct[0] );
865	<	h["errPVy"+suffix]->Fill( ferror_mct[1] );
866	<	h["errPVz"+suffix]->Fill( ferror_mct[2] );
867	<	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
868	<	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
869	<	fntrk_mct));
870	<	// Fill histo according to its track multiplicity
871	<	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
872	<	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
873	<	fntrk_mct, i);
874	<	if (saventuple_) ftree_->Fill();
875	<	suffix.clear();
876	<	}
877	<	else {  // no rec vertex found for this simvertex
878	<	if(verbose_)
879	<	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
880	<	}
881	<	}
882	<	} // === End of Looping through vertexColl
954	>	vsim!=simpv.end(); vsim++, isimpv++){
955	>	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
956	>	simx_[isimpv] = vsim->x;
957	>	simy_[isimpv] = vsim->y;
958	>	simz_[isimpv] = vsim->y;
959	>	simptsq_[isimpv] = vsim->ptsq;
960	>
961	>	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
962	>	fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
963	>	fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
964	>	fillMCHisto(vsim, isimpv, vertexColl, 3);
965		}
966		} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
967	<	} // End of Filling MC variables
967	>	} // End of   if (!realData_) {
968	>
969	>	// Finally fill the ftree_
970	>	if(saventuple_) {
971	>	ftree_->Fill();
972	>	pvtxtree_->Fill();
973	>	}
974		}
975
976		// ------------ method called once each job just before starting event loop  ------------
#	Line 894 | Line 982 | PVStudy::beginJob()
982		// ------------ method called once each job just after ending the event loop  ------------
983		void
984		PVStudy::endJob() {
985	<	if (verbose_) cout << "Analyzing PV info" << endl;
985	>	edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl;
986		// Looping through the datamodes available
987		for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
988		string suffix;
989	<	if(verbose_)  cout<<"datamode = "<< *it<<endl;
989	>	edm::LogInfo("Analysis")<<"[endJob] datamode = "<< *it<<endl;
990		switch(*it) {
991		case 0: suffix = "";
992		break;
#	Line 909 | Line 997 | PVStudy::endJob() {
997		case 3: suffix = "_mct";
998		break;
999		}
1000	<	if (analyze_) {
1000	>	suffix += "_nTrk";
1001	>	if(analyze_) {
1002		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1003	+	edm::LogInfo("Analysis") << "[endJob] ntrk = " << ntrk << endl;
1004		PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
1005	<	if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
1006	<	if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
1007	<	if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
1008	<	if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
1009	<	if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
1010	<	if ( ipv.pull_.y() > 0 ) h2["pully"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
921	<	if ( ipv.pull_.z() > 0 ) h2["pullz"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.z(), 1.);
1005	>	if ( ipv.res_.x() > 0 ) h_ana->Fill2d("resx"+suffix, ntrk,ipv.res_.x());
1006	>	if ( ipv.res_.y() > 0 ) h_ana->Fill2d("resy"+suffix, ntrk,ipv.res_.y());
1007	>	if ( ipv.res_.z() > 0 ) h_ana->Fill2d("resz"+suffix, ntrk,ipv.res_.z());
1008	>	if ( ipv.pull_.x() > 0 ) h_ana->Fill2d("pullx"+suffix, ntrk,ipv.pull_.x());
1009	>	if ( ipv.pull_.y() > 0 ) h_ana->Fill2d("pully"+suffix, ntrk,ipv.pull_.y());
1010	>	if ( ipv.pull_.z() > 0 ) h_ana->Fill2d("pullz"+suffix, ntrk,ipv.pull_.z());
1011		}
1012		}
1013		suffix.clear();
#	Line 926 | Line 1015 | PVStudy::endJob() {
1015		if(saventuple_) {
1016		file_->cd();
1017		ftree_->Write();
1018	+	pvtxtree_->Write();
1019		file_->Close();
1020		}
1021		}
#	Line 942 | Line 1032 | bool PVStudy::matchVertex(const PVStudy:
1032		}
1033
1034		// Match two reconstructed vertices
1035	<	bool PVStudy::matchVertex( const reco::Vertex & vrec1,
1036	<	const reco::Vertex & vrec2)
1035	>	bool PVStudy::matchVertex( const reco::VertexRef & vrec1,
1036	>	const reco::VertexRef & vrec2,
1037	>	int zsigncut)
1038		{
1039	<	return (fabs(vrec1.z()-vrec2.z())<0.0500); // =500um
1039	>	double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
1040	>	edm::LogInfo("Debug")<<"[matchVertex] The matching criteria in z is "<<cut<<endl;
1041	>	return (fabs(vrec1->z()-vrec2->z())<cut);
1042		}
1043
1044
1045		bool PVStudy::isResonance(const HepMC::GenParticle * p){
1046		double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
1047	<	if(verbose_) cout << "isResonance   " << p->pdg_id() << " " << ctau << endl;
1047	>	edm::LogInfo("Debug") << "[isResonance] isResonance   " << p->pdg_id() << " " << ctau << endl;
1048		return  ctau >0 && ctau <1e-6;
1049		}
1050
#	Line 1023 | Line 1116 | void PVStudy::printSimTrks(const Handle<
1116		void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1117		stringstream ss;
1118		ss << ntk;
1119	<	string suffix = ss.str();
1120	<	if(datamode == 0 )  suffix = "_" + suffix;
1121	<	if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1122	<	if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1123	<	if(datamode == 3 )  suffix = "_mct_" + suffix;
1119	>	string suffix;
1120	>	if(datamode == 0 )  suffix = "_" + ss.str();
1121	>	if(datamode == 1 )  suffix = "_spl1_mct_" + ss.str();
1122	>	if(datamode == 2 )  suffix = "_spl2_mct_" + ss.str();
1123	>	if(datamode == 3 )  suffix = "_mct_" + ss.str();
1124
1125		if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
1126		// Fill histograms of res and pull of ntk
1127	<	h["deltax"+suffix]->Fill(r.x());
1128	<	h["deltay"+suffix]->Fill(r.y());
1129	<	h["deltaz"+suffix]->Fill(r.z());
1130	<	h["pullx"+suffix]->Fill(r.x()/er.x());
1131	<	h["pully"+suffix]->Fill(r.y()/er.y());
1132	<	h["pullz"+suffix]->Fill(r.z()/er.z());
1133	<	h["errPVx"+suffix]->Fill(er.x());
1134	<	h["errPVy"+suffix]->Fill(er.y());
1135	<	h["errPVz"+suffix]->Fill(er.z());
1127	>	h_others->Fill1d("deltax"+suffix, r.x());
1128	>	h_others->Fill1d("deltay"+suffix, r.y());
1129	>	h_others->Fill1d("deltaz"+suffix, r.z());
1130	>	h_others->Fill1d("pullx"+suffix, r.x()/er.x());
1131	>	h_others->Fill1d("pully"+suffix, r.y()/er.y());
1132	>	h_others->Fill1d("pullz"+suffix, r.z()/er.z());
1133	>	h_others->Fill1d("errPVx"+suffix, er.x());
1134	>	h_others->Fill1d("errPVy"+suffix, er.y());
1135	>	h_others->Fill1d("errPVz"+suffix, er.z());
1136		}
1137
1138		PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
#	Line 1059 | Line 1152 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1152		for ( int i = 0; i < 6; ++i) {
1153		switch (i) {
1154		case 0:
1155	<	fit(h["deltax"+suffix], fpar);
1155	>	if (!h_others->ReadHisto1D("deltax"+suffix)) break;
1156	>	fit(h_others->ReadHisto1D("deltax"+suffix), fpar);
1157		data[i] = fpar[0]*cm2um;
1158		break;
1159		case 1:
1160	<	fit(h["deltay"+suffix], fpar);
1160	>	if (!h_others->ReadHisto1D("deltay"+suffix)) break;
1161	>	fit(h_others->ReadHisto1D("deltay"+suffix), fpar);
1162		data[i] = fpar[0]*cm2um;
1163		break;
1164		case 2:
1165	<	fit(h["deltaz"+suffix], fpar);
1165	>	if (!h_others->ReadHisto1D("deltaz"+suffix)) break;
1166	>	fit(h_others->ReadHisto1D("deltaz"+suffix), fpar);
1167		data[i] = fpar[0]*cm2um;
1168		break;
1169		case 3:
1170	<	fit(h["pullx"+suffix], fpar);
1170	>	if (!h_others->ReadHisto1D("pullx"+suffix)) break;
1171	>	fit(h_others->ReadHisto1D("pullx"+suffix), fpar);
1172		data[i] = fpar[0];
1173		break;
1174		case 4:
1175	<	fit(h["pully"+suffix], fpar);
1175	>	if (!h_others->ReadHisto1D("pully"+suffix)) break;
1176	>	fit(h_others->ReadHisto1D("pully"+suffix), fpar);
1177		data[i] = fpar[0];
1178		break;
1179		case 5:
1180	<	fit(h["pullz"+suffix], fpar);
1180	>	if (!h_others->ReadHisto1D("pullz"+suffix)) break;
1181	>	fit(h_others->ReadHisto1D("pullz"+suffix), fpar);
1182		data[i] = fpar[0];
1183		break;
1184		}
1185	<	if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1185	>	if (data[i] > 0) edm::LogInfo("Analysis") << "[Analysis] data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1186		}
1187		}
1188		else
#	Line 1098 | Line 1197 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1197		ntk);
1198		}
1199
1200	<	void PVStudy::fit(TH1 *&hdist, double fitPars[]){
1200	>	void PVStudy::fit(TH1 *hdist, double fitPars[]){
1201		TAxis *axis0 = hdist->GetXaxis();
1202		int nbin = axis0->GetLast();
1203		double nOF = hdist->GetBinContent(nbin+1);
#	Line 1107 | Line 1206 | void PVStudy::fit(TH1 *&hdist, double fi
1206		double fitRange = axis0->GetXmax();
1207		double sigMax[2] = {0.,0.};
1208
1209	<	if ( verbose_ ){
1210	<	cout << "Last bin = " << nbin;
1211	<	cout << "; hdist: Overflow Entries = " << nOF;
1212	<	cout << "; Underflow Entries = " << nUF;
1213	<	cout << "; hdist->GetEntries() = " << hdist->GetEntries();
1214	<	cout << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1116	<	}
1209	>	edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
1210	>	<< "; hdist: Overflow Entries = " << nOF
1211	>	<< "; Underflow Entries = " << nUF
1212	>	<< "; hdist->GetEntries() = " << hdist->GetEntries()
1213	>	<< "; fitting range = " << -fitRange << " to " << fitRange << endl;
1214	>
1215		if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit
1216		for (int bi = 0; bi < nbin; bi++) {
1217		if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) {
#	Line 1125 | Line 1223 | void PVStudy::fit(TH1 *&hdist, double fi
1223		if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1224		}
1225		}
1226	<	if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1226	>	//edm::LogInfo("Analysis") << "[Fit] Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1227
1228		TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1229		fgaus->SetParameter(1, 0.);
#	Line 1138 | Line 1236 | void PVStudy::fit(TH1 *&hdist, double fi
1236		s->SetOptStat(1111111);
1237		s->SetOptFit(0111);
1238		gPad->Update();
1239	<	if (verbose_) cout << "got function" << endl;
1239	>	//edm::LogInfo("Analysis") << "[Fit] got function" << endl;
1240		fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1241		}
1242		else
1243		fitPars[0] = -999;
1244		}
1245
1246	+	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype) {
1247	+	string suffix;
1248	+	suffix = ""; // for unsplit tracks
1249	+	if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1250	+	if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1251	+
1252	+	h_pvtrk->Fill1d("nTrk"+suffix, trackColl->size());
1253	+	for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1254	+	itTrack != trackColl->end();
1255	+	++itTrack) {
1256	+	h_pvtrk->Fill1d("trkPt"+suffix, itTrack->pt());
1257	+	h_pvtrk->Fill1d("trkDxy"+suffix, itTrack->dxy());
1258	+	h_pvtrk->Fill1d("trkDz"+suffix, itTrack->dz());
1259	+	h_pvtrk->Fill1d("trkEta"+suffix, itTrack->eta());
1260	+	h_pvtrk->Fill1d("trkPhi"+suffix, itTrack->phi());
1261	+	}
1262	+	}
1263	+
1264	+	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
1265	+	string suffix;
1266	+	suffix = ""; // for unsplit tracks
1267	+	if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1268	+	if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1269	+	int ivtx = 0;
1270	+	for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1271	+	v!=vertexColl->end(); ++v, ++ivtx) {
1272	+	if(fillNtuple) {
1273	+	if(datatype == 0) {
1274	+	nTrkPV_[ivtx] = v->tracksSize();
1275	+	sumptsq_[ivtx] = sumPtSquared(*v);
1276	+	isValid_[ivtx] = v->isValid();
1277	+	isFake_[ivtx] = v->isFake();
1278	+	recx_[ivtx] = v->x();
1279	+	recy_[ivtx] = v->y();
1280	+	recz_[ivtx] = v->z();
1281	+	recx_err_[ivtx] = v->xError();
1282	+	recy_err_[ivtx] = v->yError();
1283	+	recz_err_[ivtx] = v->zError();
1284	+	}
1285	+	if(datatype == 1) {
1286	+	nTrkPV1_spl_[ivtx] = v->tracksSize();
1287	+	sumptsq1_spl_[ivtx] = sumPtSquared(*v);
1288	+	isValid1_spl_[ivtx] = v->isValid();
1289	+	isFake1_spl_[ivtx] = v->isFake();
1290	+	recx1_spl_[ivtx] = v->x();
1291	+	recy1_spl_[ivtx] = v->y();
1292	+	recz1_spl_[ivtx] = v->z();
1293	+	recx1_err_spl_[ivtx] = v->xError();
1294	+	recy1_err_spl_[ivtx] = v->yError();
1295	+	recz1_err_spl_[ivtx] = v->zError();
1296	+
1297	+	}
1298	+	if(datatype == 2) {
1299	+	nTrkPV2_spl_[ivtx] = v->tracksSize();
1300	+	sumptsq2_spl_[ivtx] = sumPtSquared(*v);
1301	+	isValid2_spl_[ivtx] = v->isValid();
1302	+	isFake2_spl_[ivtx] = v->isFake();
1303	+	recx2_spl_[ivtx] = v->x();
1304	+	recy2_spl_[ivtx] = v->y();
1305	+	recz2_spl_[ivtx] = v->z();
1306	+	recx2_err_spl_[ivtx] = v->xError();
1307	+	recy2_err_spl_[ivtx] = v->yError();
1308	+	recz2_err_spl_[ivtx] = v->zError();
1309	+	}
1310	+	if(datatype == 4) { // for pixelVertices
1311	+	nTrkPV_pxlpvtx_[ivtx] = v->tracksSize();
1312	+	sumptsq_pxlpvtx_[ivtx] = sumPtSquared(*v);
1313	+	isValid_pxlpvtx_[ivtx] = v->isValid();
1314	+	isFake_pxlpvtx_[ivtx] = v->isFake();
1315	+	recx_pxlpvtx_[ivtx] = v->x();
1316	+	recy_pxlpvtx_[ivtx] = v->y();
1317	+	recz_pxlpvtx_[ivtx] = v->z();
1318	+	recx_err_pxlpvtx_[ivtx] = v->xError();
1319	+	recy_err_pxlpvtx_[ivtx] = v->yError();
1320	+	recz_err_pxlpvtx_[ivtx] = v->zError();
1321	+	}
1322	+	}
1323	+	}
1324	+	// For histogram only fill the leading pvtx parameters
1325	+	if(fillHisto) {
1326	+	h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1327	+	try {
1328	+	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1329	+	t!=vertexColl->begin()->tracks_end(); t++) {
1330	+	// illegal charge
1331	+	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1332	+	//      h_pvtrk->Fill1d("trkPtPV", 0.);
1333	+	}
1334	+	else {
1335	+	h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1336	+	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1337	+	h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1338	+	h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1339	+	h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1340	+	}
1341	+	}
1342	+	}
1343	+	catch (...) {
1344	+	// exception thrown when trying to use linked track
1345	+	//          h_pvtrk->Fill1d("trkPtPV", 0.);
1346	+	}
1347	+	}
1348	+
1349	+	}
1350	+
1351	+	void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
1352	+	{
1353	+	std::string suffix;
1354	+	// Set the vertexColl and histogram suffix according to datamode
1355	+	if (datatype == 1) {
1356	+	suffix = "_spl1_mct";
1357	+	nrecPV_spl1_mct_ = 0;
1358	+	}
1359	+	if (datatype == 2) {
1360	+	suffix = "_spl2_mct";
1361	+	nrecPV_spl2_mct_ = 0;
1362	+	}
1363	+	if (datatype == 3) {
1364	+	suffix = "_mct";
1365	+	nrecPV_mct_ = 0;
1366	+	}
1367	+
1368	+	//========================================================
1369	+	//  look for a matching reconstructed vertex in vertexColl
1370	+	//========================================================
1371	+
1372	+	for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1373	+	vrec!=vtxColl->end(); ++vrec){
1374	+	int nrectrk = vrec->tracksSize();
1375	+	vsim->recVtx=NULL;
1376	+
1377	+	edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1378	+	edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1379	+
1380	+	if ( matchVertex(*vsim,*vrec) && vrec->isValid() && !vrec->isFake() ) {
1381	+	vsim->recVtx=&(*vrec);
1382	+
1383	+	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1384	+	//if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
1385	+	//|| (!vsim->recVtx) )
1386	+	//vsim->recVtx=&(*vrec);
1387	+	//}
1388	+	edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1389	+
1390	+	double fres_mct[3];
1391	+	double ferror_mct[3];
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 ftree_
1419	+	if(datatype == 1) {
1420	+	nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1421	+	deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1422	+	deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1423	+	deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1424	+	pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
1425	+	pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1426	+	pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1427	+	nrecPV_spl1_mct_++;
1428	+	}
1429	+
1430	+	if(datatype == 2) {
1431	+	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1432	+	deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1433	+	deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1434	+	deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1435	+	pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
1436	+	pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1437	+	pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1438	+	nrecPV_spl2_mct_++;
1439	+	}
1440	+	if(datatype == 3) {
1441	+	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1442	+	deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1443	+	deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1444	+	deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1445	+	pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
1446	+	pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1447	+	pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1448	+	nrecPV_mct_++;
1449	+	}
1450	+	} // End of  if(saventuple_) {
1451	+	} //  if ( matchVertex(*vsim,*vrec) ) {
1452	+	else {  // no rec vertex found for this simvertex
1453	+	edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1454	+	}
1455	+	}
1456	+	}
1457	+
1458	+
1459		void PVStudy::SetVarToZero() {
1460	+	//Greg's variables
1461	+	fntrk_ = 0;
1462	+	//pvtx position (x,y,z) residual and error
1463	+	for(int i = 0; i<3;i++)    {
1464	+	fres_[i] = 0;
1465	+	ferror_[i] = 0;
1466	+	}
1467
1468	<	// number of reconstructed vertices
1468	>	//Yanyan's variables
1469	>	// Number of reconstructed vertices
1470		nrecPV_ = 0;
1471	<	nrecPV1_spl_ = 0;
1471	>	nrecPV1_spl_ = 0;
1472		nrecPV2_spl_ = 0;
1473	<	// number of tracks in the vertex
1474	<	fntrk_ = 0;
1475	<	fntrk_spl1_mct_ = 0;
1476	<	fntrk_spl2_mct_ = 0;
1477	<	fntrk_mct_ = 0;
1478	<	//pvtx position (x,y,z) residual and error
1479	<	for(int i = 0; i<3;i++)
1480	<	{
1481	<	fres_[i] = 0;
1482	<	ferror_[i] = 0;
1483	<	fres_spl1_mct_[i] = 0;
1484	<	ferror_spl1_mct_[i] = 0;
1485	<	fres_spl2_mct_[i] = 0;
1486	<	ferror_spl2_mct_[i] = 0;
1487	<	fres_mct_[i] = 0;
1488	<	ferror_mct_[i] = 0;
1489	<	}
1473	>	nrecPV_twovtx_ = 0;
1474	>	nsimPV_ = 0;
1475	>
1476	>	nrecPV_mct_ = 0;
1477	>	nrecPV_spl1_mct_ = 0;
1478	>	nrecPV_spl2_mct_ = 0;
1479	>
1480	>	// Mininum separation between the secondary pvtxes and leading pvtx
1481	>	min_zsep_ = 9999.0;
1482	>	min_ntrksep_ = 9999.0;
1483	>	min_sumpt2sep_ = -9999.0;
1484	>
1485	>
1486	>	for (int i = 0; i < nMaxPVs_; i++) {
1487	>	// recoVertices with all tracks
1488	>	nTrkPV_[i] = 0; // Number of tracks in the pvtx
1489	>	sumptsq_[i] = 0;
1490	>	isValid_[i] = -1;
1491	>	isFake_[i] = -1;
1492	>	recx_[i] = 0;
1493	>	recy_[i] = 0;
1494	>	recz_[i] = 0;
1495	>	recx_err_[i] = 0;
1496	>	recy_err_[i] = 0;
1497	>	recz_err_[i] = 0;
1498	>
1499	>	// recoVertices with splitTrack1
1500	>	nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx
1501	>	sumptsq1_spl_[i] = 0;
1502	>	isValid1_spl_[i] = -1;
1503	>	isFake1_spl_[i] = -1;
1504	>	recx1_spl_[i] = 0;
1505	>	recy1_spl_[i] = 0;
1506	>	recz1_spl_[i] = 0;
1507	>	recx1_err_spl_[i] = 0;
1508	>	recy1_err_spl_[i] = 0;
1509	>	recz1_err_spl_[i] = 0;
1510	>
1511	>	// recoVertices with splitTrack2
1512	>	nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx
1513	>	sumptsq2_spl_[i] = 0;
1514	>	isValid2_spl_[i] = -1;
1515	>	isFake2_spl_[i] = -1;
1516	>	recx2_spl_[i] = 0;
1517	>	recy2_spl_[i] = 0;
1518	>	recz2_spl_[i] = 0;
1519	>	recx2_err_spl_[i] = 0;
1520	>	recy2_err_spl_[i] = 0;
1521	>	recz2_err_spl_[i] = 0;
1522	>
1523	>	//pixelVertices
1524	>	nTrkPV_pxlpvtx_[i] = 0; // Number of tracks in the pvtx
1525	>	sumptsq_pxlpvtx_[i] = 0;
1526	>	isValid_pxlpvtx_[i] = -1;
1527	>	isFake_pxlpvtx_[i] = -1;
1528	>	recx_pxlpvtx_[i] = 0;
1529	>	recy_pxlpvtx_[i] = 0;
1530	>	recz_pxlpvtx_[i] = 0;
1531	>	recx_err_pxlpvtx_[i] = 0;
1532	>	recy_err_pxlpvtx_[i] = 0;
1533	>	recz_err_pxlpvtx_[i] = 0;
1534	>
1535	>	// matched two-vertices
1536	>	nTrkPV1_spl_twovtx_[i] = 0;
1537	>	nTrkPV2_spl_twovtx_[i] = 0;
1538	>	nTrkPV_twovtx_[i] = 0;
1539	>	sumptsq1_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1540	>	sumptsq2_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1541	>	deltax_twovtx_[i] = 0;
1542	>	deltay_twovtx_[i] = 0;
1543	>	deltaz_twovtx_[i] = 0;
1544	>	errx_twovtx_[i] = 0;
1545	>	erry_twovtx_[i] = 0;
1546	>	errz_twovtx_[i] = 0;
1547	>	pullx_twovtx_[i] = 0;
1548	>	pully_twovtx_[i] = 0;
1549	>	pullz_twovtx_[i] = 0;
1550	>
1551	>	//simpv
1552	>	nsimTrkPV_[i] = 0;
1553	>	simx_[i] = 0;
1554	>	simy_[i] = 0;
1555	>	simz_[i] = 0;
1556	>	simptsq_[i] = 0;
1557	>
1558	>	// residual and pulls with mc truth required
1559	>	deltax_mct_[i] = 0;
1560	>	deltay_mct_[i] = 0;
1561	>	deltaz_mct_[i] = 0;
1562	>	pullx_mct_[i] = 0;
1563	>	pully_mct_[i] = 0;
1564	>	pullz_mct_[i] = 0;
1565	>
1566	>	deltax_spl1_mct_[i] = 0;
1567	>	deltay_spl1_mct_[i] = 0;
1568	>	deltaz_spl1_mct_[i] = 0;
1569	>	pullx_spl1_mct_[i] = 0;
1570	>	pully_spl1_mct_[i] = 0;
1571	>	pullz_spl1_mct_[i] = 0;
1572	>
1573	>	deltax_spl2_mct_[i] = 0;
1574	>	deltay_spl2_mct_[i] = 0;
1575	>	deltaz_spl2_mct_[i] = 0;
1576	>	pullx_spl2_mct_[i] = 0;
1577	>	pully_spl2_mct_[i] = 0;
1578	>	pullz_spl2_mct_[i] = 0;
1579	>	}
1580	>
1581		}
1582
1583	+	double PVStudy::sumPtSquared(const reco::Vertex & v)  {
1584	+	double sum = 0.;
1585	+	double pT;
1586	+	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1587	+	pT = (**it).pt();
1588	+	sum += pT*pT;
1589	+	}
1590	+	return sum;
1591	+	}
1592

Diff Legend

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