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.9 by yygao, Thu Oct 29 02:05:54 2009 UTC vs.
Revision 1.14 by jengbou, Fri Dec 4 20:48:04 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 66 | Line 66 | PVStudy::PVStudy(const edm::ParameterSet
66		//=======================================================================
67		// Get configuration for TrackTupleMaker
68		//=======================================================================
69	<	simG4_           = iConfig.getParameter<edm::InputTag>( "simG4" );
70	<	trackCollectionTag_     = iConfig.getUntrackedParameter<edm::InputTag>("trackCollection");
71	<	splitTrackCollection1Tag_     = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection1");
72	<	splitTrackCollection2Tag_     = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
69	>	simG4_                     = iConfig.getParameter<edm::InputTag>( "simG4" );
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	<	verbose_         = iConfig.getUntrackedParameter<bool>("verbose",false);
77	<	realData_        = iConfig.getUntrackedParameter<bool>("realData",false);
78	<	analyze_         = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
79	<	saventuple_         = iConfig.getUntrackedParameter<bool>("saventuple",false);
80	<	outputfilename_  = iConfig.getUntrackedParameter<string>("OutputFileName");
81	<	ntrkdiffcut_         = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
82	<	zsigncut_         = iConfig.getUntrackedParameter<int>("zsigncut");
83	<	nTrkMin_ =  iConfig.getUntrackedParameter<int>("nTrkMin");
84	<	nTrkMax_ =  iConfig.getUntrackedParameter<int>("nTrkMax");
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	>	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();
#	Line 110 | Line 112 | PVStudy::PVStudy(const edm::ParameterSet
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");
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");
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");
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");
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");
#	Line 187 | Line 218 | PVStudy::PVStudy(const edm::ParameterSet
218		}
219		}
220
190	–	edm::Service<TFileService> fs;
191	–	TFileDirectory subDir = fs->mkdir( "Summary" );
192	–	TFileDirectory subDir1 = fs->mkdir( "Others" );
193	–	//   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
194	–
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	<
218	<	h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
219	<	h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
220	<	h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#eta of rec tracks in PV"+suffix),100,-3,3);
221	<	h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
222	<	h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);
223	<	h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
224	<	h["nrecPV"+suffix]   = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);
225	<	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	>
231	>	// Book Histograms:
232	>	h_pvtrk   = new PVHistograms();
233	>	h_pixvtx  = new PVHistograms();
234	>	h_misc    = new PVHistograms();
235	>	h_summary = new PVHistograms();
236	>	h_others  = new PVHistograms();
237	>
238	>	for(int i=0;i<3;i++) {
239	>	if(i == 0) h_pvtrk->Init("pvTrk");
240	>	else {
241	>	stringstream ss;
242	>	ss << i;
243	>	h_pvtrk->Init("pvTrk",ss.str(),"spl");
244		}
245	<	h["nrecPVDiff"]     = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
246	<	h["nTrkPVDiff"]     = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
247	<	h["nTrkPVRelDiff"]  = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
248	<
245	>	}
246	>	h_pixvtx->Init("pixVtx");
247	>	h_misc->Init("misc");
248	>
249	>	// Book MC only plots
250	>	if (!realData_) {
251	>	h_gen = new PVHistograms();
252	>	h_gen->Init("generator");
253	>	}
254	>	if (analyze_) h_ana = new PVHistograms();
255	>
256		//Book histograms sensitive to data/mc
257		for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
258		string suffix;
259	<	if(verbose_)  cout<<"datamode = "<< *it<<endl;
259	>	edm::LogInfo("Debug")<<"datamode = "<< *it<<endl;
260		switch(*it) {
261		case 0: suffix = "";
262		break;
#	Line 242 | Line 267 | PVStudy::PVStudy(const edm::ParameterSet
267		case 3: suffix = "_mct";
268		break;
269		}
270	<	h["deltax"+suffix] = subDir.make<TH1D>(TString("deltax"+suffix), TString("x-residual pvtx"+suffix),800,-0.04,0.04);
246	<	h["deltay"+suffix] = subDir.make<TH1D>(TString("deltay"+suffix), TString("y-residual pvtx"+suffix),800,-0.04,0.04);
247	<	h["deltaz"+suffix] = subDir.make<TH1D>(TString("deltaz"+suffix), TString("z-residual pvtx"+suffix),800,-0.04,0.04);
248	<	h["pullx"+suffix]  = subDir.make<TH1D>(TString("pullx"+suffix), TString("x-pull pvtx"+suffix),800,-10,10);
249	<	h["pully"+suffix]  = subDir.make<TH1D>(TString("pully"+suffix), TString("y-pull pvtx"+suffix),800,-10,10);
250	<	h["pullz"+suffix]  = subDir.make<TH1D>(TString("pullz"+suffix), TString("z-pull pvtx"+suffix),800,-10,10);
251	<	h["errPVx"+suffix] = subDir.make<TH1D>(TString("errPVx"+suffix), TString("X"+suffix+" vertex error"),100,0.,0.02);
252	<	h["errPVy"+suffix] = subDir.make<TH1D>(TString("errPVy"+suffix), TString("Y"+suffix+" vertex error"),100,0.,0.02);
253	<	h["errPVz"+suffix] = subDir.make<TH1D>(TString("errPVz"+suffix), TString("Z"+suffix+" vertex error"),100,0.,0.02);
270	>	h_summary->Init("summary",suffix);
271
272		// Book residual, error and pull histograms for each nTrk bin
273		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
274		stringstream ss;
275		ss << ntrk;
276	<	string suffix2 = suffix+"_"+ss.str();
277	<	h["deltax"+suffix2]  = subDir1.make<TH1D>(TString("deltax"+suffix2), TString("x-residual of pvtx"+suffix),100,-0.02,0.02);
278	<	h["deltax"+suffix2]->GetXaxis()->SetTitle("cm");
262	<	h["deltay"+suffix2]  = subDir1.make<TH1D>(TString("deltay"+suffix2), TString("y-residual of pvtx"+suffix),100,-0.02,0.02);
263	<	h["deltay"+suffix2]->GetXaxis()->SetTitle("cm");
264	<	h["deltaz"+suffix2]  = subDir1.make<TH1D>(TString("deltaz"+suffix2), TString("z-residual of pvtx"+suffix),100,-0.02,0.02);
265	<	h["deltaz"+suffix2]->GetXaxis()->SetTitle("cm");
266	<	h["pullx"+suffix2] = subDir1.make<TH1D>(TString("pullx"+suffix2), TString("x-pull of pvtx"+suffix),100,-10.,10.);
267	<	h["pully"+suffix2] = subDir1.make<TH1D>(TString("pully"+suffix2), TString("y-pull of pvtx"+suffix),100,-10.,10.);
268	<	h["pullz"+suffix2] = subDir1.make<TH1D>(TString("pullz"+suffix2), TString("z-pull of pvtx"+suffix),100,-10.,10.);
269	<	h["errPVx"+suffix2] = subDir1.make<TH1D>(TString("errPVx"+suffix2), TString("X"+suffix+" vertex error"),100,0.,0.02);
270	<	h["errPVy"+suffix2] = subDir1.make<TH1D>(TString("errPVy"+suffix2), TString("Y"+suffix+" vertex error"),100,0.,0.02);
271	<	h["errPVz"+suffix2] = subDir1.make<TH1D>(TString("errPVz"+suffix2), TString("Z"+suffix+" vertex error"),100,0.,0.02);
272	<	suffix2.clear();
273	<	} // End of   for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
274	<
276	>	h_others->Init("others",suffix,ss.str());
277	>	}
278	>
279		// Book residual and pull histograms only when analyzeOntheFly is enabled
280	<	if(analyze_) {
277	<	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);
278	<	h2["resx"+suffix+"_nTrk"]->SetMarkerStyle(21);
279	<	h2["resx"+suffix+"_nTrk"]->SetMarkerColor(4);
280	<	h2["resx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
281	<	h2["resx"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
282	<	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);
283	<	h2["resy"+suffix+"_nTrk"]->SetMarkerStyle(21);
284	<	h2["resy"+suffix+"_nTrk"]->SetMarkerColor(4);
285	<	h2["resy"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
286	<	h2["resy"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
287	<	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);
288	<	h2["resz"+suffix+"_nTrk"]->SetMarkerStyle(21);
289	<	h2["resz"+suffix+"_nTrk"]->SetMarkerColor(4);
290	<	h2["resz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
291	<	h2["resz"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
292	<	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.);
293	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerStyle(21);
294	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerColor(4);
295	<	h2["pullx"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
296	<	h2["pullx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
297	<	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.);
298	<	h2["pully"+suffix+"_nTrk"]->SetMarkerStyle(21);
299	<	h2["pully"+suffix+"_nTrk"]->SetMarkerColor(4);
300	<	h2["pully"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
301	<	h2["pully"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
302	<	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.);
303	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerStyle(21);
304	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerColor(4);
305	<	h2["pullz"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
306	<	h2["pullz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
307	<	} // End of   if(analyze_) {
280	>	if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
281		suffix.clear();
282	<	} // End of Book histograms sensitive to data/mc
283	<	// Book MC only plots
311	<	if (!realData_) {
312	<	h["genPart_T"]      = subDir.make<TH1D>("genPart_T","t component of gen particles",300,-0.5,0.5);
313	<	h["genPart_T"]->GetXaxis()->SetTitle("t (nanosecond)");
314	<	h["genPart_cT"]     = subDir.make<TH1D>("genPart_cT","ct component of gen particles",300,-150.,150.);
315	<	h["genPart_cT"]->GetXaxis()->SetTitle("ct (mm)");
316	<	h["nsimPV"]         = subDir.make<TH1D>("nsimPV","Num of sim PV",51,-0.5,50.5);
317	<	}
282	>	} // End of Book histograms sensitive to data/mc
283	>
284		}
285
286		PVStudy::~PVStudy()
287		{
288	<
289	<	// do anything here that needs to be done at desctruction time
290	<	// (e.g. close files, deallocate resources etc.)
288	>	theFile->cd();
289	>	theFile->cd("Summary");
290	>	h_pvtrk->Save();
291	>	h_pixvtx->Save();
292	>	h_gen->Save();
293	>	h_misc->Save();
294	>	h_summary->Save();
295	>	if (analyze_)
296	>	h_ana->Save();
297	>	theFile->cd();
298	>	theFile->cd("Others");
299	>	h_others->Save();
300
301	+	theFile->Close();
302		}
303
304		void PVStudy::setRootStyle() {
#	Line 362 | Line 338 | std::vector<PVStudy::simPrimaryVertex> P
338		std::vector<PVStudy::simPrimaryVertex> simpv;
339		const HepMC::GenEvent* evt=evtMC->GetEvent();
340		if (evt) {
341	<	if(verbose_) {
342	<	cout << "process id " << evt->signal_process_id()<<endl;
343	<	cout << "signal process vertex " << ( evt->signal_process_vertex() ?
344	<	evt->signal_process_vertex()->barcode() : 0 )   <<endl;
369	<	cout << "number of vertices " << evt->vertices_size() << endl;
370	<	}
341	>	edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
342	>	edm::LogInfo("SimPVs") << "[getSimPVs] signal process vertex " << ( evt->signal_process_vertex() ?
343	>	evt->signal_process_vertex()->barcode() : 0 ) <<endl;
344	>	edm::LogInfo("SimPVs") << "[getSimPVs] number of vertices " << evt->vertices_size() << endl;
345
346		int idx=0; int npv=0;
347		for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
#	Line 378 | Line 352 | std::vector<PVStudy::simPrimaryVertex> P
352		// t component of PV:
353		for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
354		mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
355	<	//        std::cout << "Status = " << (*mother)->status() << std::endl;
355	>	//        edm::LogInfo("SimPVs") << "Status = " << (*mother)->status() << endl;
356		HepMC::GenVertex * mv=(*mother)->production_vertex();
357		if( ((*mother)->status() == 3) && (!mv)) {
358	<	//        std::cout << "npv= " << npv << std::endl;
358	>	//        edm::LogInfo("SimPVs") << "npv= " << npv << endl;
359		if (npv == 0) {
360	<	h["genPart_cT"]->Fill(pos.t()); // mm
361	<	h["genPart_T"]->Fill(pos.t()/299.792458); // ns
360	>	h_gen->Fill1d("genPart_cT", pos.t()); // mm
361	>	h_gen->Fill1d("genPart_T", pos.t()/299.792458); // ns
362		}
363		npv++;
364		}
#	Line 392 | Line 366 | std::vector<PVStudy::simPrimaryVertex> P
366		//       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
367
368		bool hasMotherVertex=false;
369	<	if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl;
369	>	if (verbose_) cout << "[getSimPVs] mothers of vertex[" << ++idx << "]: " << endl;
370		for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
371		mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
372		HepMC::GenVertex * mv=(*mother)->production_vertex();
373	<	//        std::cout << "Status = " << (*mother)->status() << std::endl;
373	>	// if (verbose_) cout << "Status = " << (*mother)->status() << endl;
374		if (mv) {
375		hasMotherVertex=true;
376		if(!verbose_) break; //if verbose_, print all particles of gen vertices
#	Line 423 | Line 397 | std::vector<PVStudy::simPrimaryVertex> P
397
398		if(!vp){
399		// this is a new vertex
400	<	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
400	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
401		simpv.push_back(sv);
402		vp=&simpv.back();
403		}else{
404	<	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
404	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
405		}
406		vp->genVertex.push_back((*vitr)->barcode());
407		// collect final state descendants
#	Line 451 | Line 425 | std::vector<PVStudy::simPrimaryVertex> P
425		}
426		}
427		}
428	<	}
429	<	}
428	>	}//loop MC vertices daughters
429	>	}//loop MC vertices
430		}
431		return simpv;
432		}
#	Line 485 | Line 459 | std::vector<PVStudy::simPrimaryVertex> P
459		ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
460		}
461		}
462	<	primary =  ( gv->position().t()==0);
462	>	//primary =  ( gv->position().t()==0 );
463	>	primary = true;
464	>	h_gen->Fill1d("genPart_cT", gv->position().t()); // mm
465	>	h_gen->Fill1d("genPart_T", gv->position().t()/299.792458); // ns
466	>
467		//resonance= ( gp->mother() && isResonance(gp->mother()));  // in CLHEP/HepMC days
468		// no more mother pointer in the improved HepMC GenParticle
469		resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
#	Line 534 | Line 512 | PVStudy::analyze(const edm::Event& iEven
512		using namespace edm;
513		using namespace reco;
514
515	<	if (verbose_)
538	<	cout<<"PVStudy::analyze():"<<endl;
515	>	edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
516		//=======================================================
517		// Initialize Root-tuple variables if needed
518		//=======================================================
#	Line 553 | Line 530 | PVStudy::analyze(const edm::Event& iEven
530		if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
531		trackColl = trackCollectionHandle.product();
532		} else {
533	<	cout << "trackCollection cannot be found -> using empty collection of same type." <<endl;
533	>	edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
534		}
535	+
536		//splitTrackCollection1
537		static const reco::TrackCollection s_empty_splitTrackColl1;
538		const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
#	Line 563 | Line 541 | PVStudy::analyze(const edm::Event& iEven
541		if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
542		splitTrackColl1 = splitTrackCollection1Handle.product();
543		} else {
544	<	cout << "splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
544	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
545		}
546	+
547		//splitTrackCollection2
548		static const reco::TrackCollection s_empty_splitTrackColl2;
549		const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
#	Line 573 | Line 552 | PVStudy::analyze(const edm::Event& iEven
552		if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
553		splitTrackColl2 = splitTrackCollection2Handle.product();
554		} else {
555	<	cout << "splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
555	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
556		}
557	<
557	>
558	>
559	>	//=======================================================
560	>	// Fill trackparameters of the input tracks to pvtx fitter
561	>	//=======================================================
562	>	edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
563	>	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype)
564	>	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
565	>	fillTrackHisto(trackColl, 0);
566	>	fillTrackHisto(splitTrackColl1, 1);
567	>	fillTrackHisto(splitTrackColl2, 2);
568	>	edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
569	>
570		//=======================================================
571		// PVTX accessors
572		//=======================================================
#	Line 587 | Line 578 | PVStudy::analyze(const edm::Event& iEven
578		if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
579		vertexColl = vertexCollectionHandle.product();
580		} else {
581	<	cout << "vertexCollection cannot be found -> using empty collection of same type." <<endl;
581	>	edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
582		}
583	+
584		//splitVertexCollection1
585		static const reco::VertexCollection s_empty_splitVertexColl1;
586		const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
#	Line 597 | Line 589 | PVStudy::analyze(const edm::Event& iEven
589		if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
590		splitVertexColl1 = splitVertexCollection1Handle.product();
591		} else {
592	<	cout << "splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
593	<	}
592	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
593	>	}
594	>
595		//splitVertexCollection2
596		static const reco::VertexCollection s_empty_splitVertexColl2;
597		const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
#	Line 607 | Line 600 | PVStudy::analyze(const edm::Event& iEven
600		if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
601		splitVertexColl2 = splitVertexCollection2Handle.product();
602		} else {
603	<	cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
603	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
604		}
605	<
606	<	if(verbose_) cout<<"End accessing the track and vertex collections"<<endl;
605	>
606	>	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track and vertex collections"<<endl;
607
608		//=======================================================
609		// MC simvtx accessor
#	Line 629 | Line 622 | PVStudy::analyze(const edm::Event& iEven
622		try{
623		iSetup.getData(pdt);
624		}catch(const Exception&){
625	<	cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
625	>	edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
626	>	}
627	>
628	>	//=======================================================
629	>	// GET pixelVertices
630	>	//=======================================================
631	>	static const reco::VertexCollection s_empty_pixelVertexColl;
632	>	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
633	>	edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
634	>	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
635	>	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
636	>	pixelVertexColl = pixelVertexCollectionHandle.product();
637	>	} else {
638	>	edm::LogInfo("Debug") << "[PVStudy] pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
639		}
640
641		//=======================================================
642	<	// Fill trackparameters of the input tracks to pvtx fitter
642	>	// Fill pixelVertices related histograms
643		//=======================================================
644	<	if(verbose_)
645	<	cout<<"Start filling track parameters of the input tracks to pvtx fitter."<<endl;
646	<	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype)
647	<	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
648	<	fillTrackHisto(trackColl, 0);
649	<	fillTrackHisto(splitTrackColl1, 1);
650	<	fillTrackHisto(splitTrackColl2, 2);
651	<	if(verbose_)
652	<	cout<<"End filling track parameters of the input tracks to pvtx fitter."<<endl;
644	>	nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
645	>	if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
646	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
647	>	fillTrackHistoInPV(pixelVertexColl, 4, false, true);
648	>	h_pixvtx->Fill1d("dzErr_pxlpvtx", pixelVertexColl->begin()->zError());
649	>	// Get the dZ error of the tracks in the leading pixelVertexColl
650	>	for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
651	>	t!= (pixelVertexColl->begin())->tracks_end(); t++) {
652	>
653	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
654	>	//        h_pvtrk->Fill1d("trkPtPV", 0.);
655	>	}
656	>	else
657	>	h_pixvtx->Fill1d("trkdzErr_pxlpvtx", (**t).dzError());
658	>	}
659	>	// Fill the track->dz() in the PV difference from first pixelpvtx
660	>	for(reco::Vertex::trackRef_iterator t = (vertexColl->begin())->tracks_begin();
661	>	t!= (vertexColl->begin())->tracks_end(); t++) {
662	>	// illegal charge
663	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
664	>	//        h_pvtrk->Fill1d("trkPtPV", 0.);
665	>	}
666	>	else {
667	>	if(pixelVertexColl->size()>0) {
668	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz", (**t).dz() -  pixelVertexColl->begin()->z());
669	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_pxlpvtxdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
670	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_trkdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
671	>	h_pixvtx->Fill1d("trkdzErr_pvtx", (**t).dzError());
672	>	}
673	>	}
674	>	}
675	>	}
676
677		//=======================================================
678	<	// Fill reconstructed vertices ootb
678	>	// Fill number of reconstructed vertices
679		//=======================================================
680	<	if(verbose_)  {
681	<	cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
680	>
681	>	edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
682	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
683	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
684	>	edm::LogInfo("Debug")<<"[PVStudy] Start filling pvtx parameters."<<endl;
685	>	if (verbose_) {
686		printRecVtxs(vertexCollectionHandle);
654	–	cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
687		printRecVtxs(splitVertexCollection1Handle);
656	–	cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
688		printRecVtxs(splitVertexCollection2Handle);
658	–	cout<<"Start filling pvtx parameters."<<endl;
689		}
660	–
690		nrecPV_ = int (vertexColl->size());
691		nrecPV1_spl_ = int (splitVertexColl1->size());
692		nrecPV2_spl_ = int (splitVertexColl2->size());
693
694	<	h["nrecPV"]->Fill(nrecPV_);
695	<	h["nrecPV1_spl"]->Fill(nrecPV1_spl_);
696	<	h["nrecPV2_spl"]->Fill(nrecPV2_spl_);
697	<	h["nrecPVDiff"]->Fill(double(nrecPV1_spl_)-double(nrecPV2_spl_));
698	<
694	>	h_pvtrk->Fill1d("nrecPV", nrecPV_);
695	>	h_pvtrk->Fill1d("nrecPV1_spl", nrecPV1_spl_);
696	>	h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
697	>	h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
698	>
699	>	//=================================================================
700	>	// Fill track parameter ntuple/hist for tracks used in recoVertices
701	>	//=================================================================
702	>
703		//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
704	<	fillTrackHistoInPV(vertexColl, 0, true, true);
705	<	fillTrackHistoInPV(splitVertexColl1, 1, false, true);
706	<	fillTrackHistoInPV(splitVertexColl2, 2, false, true);
707	<	if(verbose_)
708	<	cout<<"End filling pvtx parameters."<<endl;
704	>	if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
705	>	fillTrackHistoInPV(vertexColl, 0, true, true);
706	>
707	>	if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
708	>	fillTrackHistoInPV(splitVertexColl1, 1, false, true);
709	>
710	>	if(splitVertexColl1->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
711	>	fillTrackHistoInPV(splitVertexColl2, 2, false, true);
712	>
713	>
714	>	//=======================================================
715	>	// Compare offlinePrimaryVertices with pixelVertices
716	>	//=======================================================
717	>	if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
718	>	&& (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {
719	>	h_pixvtx->Fill1d("nrecPV_minus_nrecPxlPV", double (nrecPV_ - nrecPV_pxlpvtx_));
720	>	// difference in reconstructed position of the leading pvtx
721	>	//edm::LogInfo("Debug")<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;
722	>	//edm::LogInfo("Debug")<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
723	>	h_pixvtx->Fill1d("recxPV_minus_recxPxlPV", recx_[0] - recx_pxlpvtx_[0]);
724	>	h_pixvtx->Fill1d("recyPV_minus_recyPxlPV", recy_[0] - recy_pxlpvtx_[0]);
725	>	h_pixvtx->Fill1d("reczPV_minus_reczPxlPV", recz_[0] - recz_pxlpvtx_[0]);
726	>	}
727	>
728	>
729	>
730	>	//==========================================================
731	>	// Fill secondary/primary min separations for multi-vertices
732	>	//==========================================================
733	>
734	>	if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {
735	>
736	>	double min_xsep = 9999.0;
737	>	double min_ysep = 9999.0;
738	>	double min_zsep = 9999.0;
739	>	double min_xsep_sign = 9999.0;
740	>	double min_ysep_sign = 9999.0;
741	>	double min_zsep_sign = 9999.0;
742	>	double min_ntrksep = 9999.0;
743	>	double min_sumpt2sep = 9999999.0;
744	>
745	>	edm::LogInfo("Debug")<<"[PVStudy] leading pvtx z = "<< vertexColl->begin()->z() <<endl;
746	>
747	>	// Looping through the secondary vertices to find the mininum separation to the primary
748	>	for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
749	>	v!=vertexColl->end(); ++v) {
750	>	if(v->isValid() && ! v->isFake()) {
751	>	// xsep
752	>	if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
753	>	min_xsep = fabs(v->x()- vertexColl->begin()->x());
754	>	// ysep
755	>	if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
756	>	min_ysep = fabs(v->y()- vertexColl->begin()->y());
757	>	// zsep
758	>	if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
759	>	min_zsep = fabs(v->z()- vertexColl->begin()->z());
760	>	// xsep_sign
761	>	double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
762	>	if(xsep_sign < min_xsep_sign)
763	>	min_xsep_sign = xsep_sign;
764	>	// ysep_sign
765	>	double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
766	>	if(ysep_sign < min_ysep_sign)
767	>	min_ysep_sign = ysep_sign;
768	>	// zsep_sign
769	>	double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
770	>	if(zsep_sign < min_zsep_sign)
771	>	min_zsep_sign = zsep_sign;
772	>	// ntrksep
773	>	if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
774	>	min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
775	>	// sumpt2sep
776	>	if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
777	>	min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
778	>	}
779	>	}
780	>	h_misc->Fill1d("min_xsep", min_xsep);
781	>	h_misc->Fill1d("min_ysep", min_ysep);
782	>	h_misc->Fill1d("min_zsep", min_zsep);
783	>	h_misc->Fill1d("min_xsep_sign", min_xsep_sign);
784	>	h_misc->Fill1d("min_ysep_sign", min_ysep_sign);
785	>	h_misc->Fill1d("min_zsep_sign", min_zsep_sign);
786	>	h_misc->Fill1d("min_ntrksep", min_ntrksep);
787	>	h_misc->Fill1d("min_sumpt2sep", min_sumpt2sep);
788	>
789	>	min_zsep_ = min_zsep;
790	>	min_ntrksep_ = min_ntrksep;
791	>	min_sumpt2sep_ = min_sumpt2sep;
792	>
793	>
794	>	} // End of  if(vertexColl->size()>1) {
795	>
796	>	edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
797
798		//=======================================================
799		//           PrimaryVertex Matching
#	Line 682 | Line 803 | PVStudy::analyze(const edm::Event& iEven
803		// since vertexColl are sorted in decreasing order of sum(pT)
804		//=======================================================
805
806	<	if(verbose_)   cout<<"PVStudy::analyze(): matching pvtxs "<<endl;
806	>	edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
807		reco::VertexCollection s_empty_matchedVertexColl1;
808		reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
809		matchedVertexColl1->reserve(splitVertexColl1->size());
#	Line 690 | Line 811 | PVStudy::analyze(const edm::Event& iEven
811		reco::VertexCollection *matchedVertexColl2 = &s_empty_matchedVertexColl2;
812		matchedVertexColl2->reserve(splitVertexColl2->size());
813
814	<	for (size_t ivtx = 0; ivtx<splitVertexCollection1Handle->size(); ++ivtx) {
814	>	bool stopmatching = false;
815	>	for (size_t ivtx = 0; ivtx<splitVertexCollection1Handle->size() && !stopmatching; ++ivtx) {
816		reco::VertexRef recvtx1(splitVertexCollection1Handle, ivtx);
817	<	for (size_t jvtx = ivtx; jvtx < splitVertexCollection2Handle->size(); ++jvtx) {
817	>	if( !(recvtx1->isValid()) || recvtx1->isFake()) break;
818	>	for (size_t jvtx = ivtx; jvtx < splitVertexCollection2Handle->size(); ++jvtx) {
819	>	//------------------------------------------------------------------------
820	>	//== If only considering matching the first vertex of the splitVertexColl
821	>	//------------------------------------------------------------------------
822	>	if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
823		reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
824	+	if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
825	+	edm::LogInfo("Debug")<<"[PVStudy] Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
826	+	edm::LogInfo("Debug")<<"[PVStudy] recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
827	+	edm::LogInfo("Debug")<<"[PVStudy] recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
828	+
829		if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
830	<	if(verbose_) {
831	<	cout<<"The two splitted vertices match in Z: "<<endl;
700	<	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
701	<	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
702	<	}
830	>	edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
831	>
832		int nTrkPV1 = recvtx1->tracksSize();
833		int nTrkPV2 = recvtx2->tracksSize();
834		double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
835	<	h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
836	<	h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
835	>	h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
836	>	h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
837		if(fabs(ntrkreldiff)<ntrkdiffcut_) {
838	<	if(verbose_) {
839	<	cout<<"(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
711	<	}
838	>	edm::LogInfo("Debug")<<"[PVStudy] (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
839	>
840		matchedVertexColl1->push_back(*recvtx1);
841		matchedVertexColl2->push_back(*recvtx2);
842	<	continue;
842	>	stopmatching = true; // stop the matching when the first match is found!
843	>	break;
844		}
845		else
846	<	cout<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
846	>	edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
847		}
848		else {
849	<	cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
850	<	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
851	<	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
849	>	edm::LogInfo("Debug")<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
850	>	edm::LogInfo("Debug")<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
851	>	edm::LogInfo("Debug")<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
852		}
853		}
854		}
855	<	if(verbose_)     cout<<"PVStudy::analyze(): End matching pvtxs"<<endl;
855	>	edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
856
857		//=======================================================
858		//  Analyze matchedVertexColls
859		//=======================================================
860	<	if(verbose_) {
732	<	cout<<"Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
733	<	}
734	<
860	>	edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
861
862		// Compare the reconstructed vertex position and calculate resolution/pulls
863		if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
#	Line 753 | Line 879 | PVStudy::analyze(const edm::Event& iEven
879		ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
880		fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
881
882	<	h["deltax"]->Fill( fres_[0] );
883	<	h["deltay"]->Fill( fres_[1] );
884	<	h["deltaz"]->Fill( fres_[2] );
885	<	h["pullx"]->Fill( fres_[0]/ferror_[0]);
886	<	h["pully"]->Fill( fres_[1]/ferror_[1]);
887	<	h["pullz"]->Fill( fres_[2]/ferror_[2]);
888	<	h["errPVx"]->Fill( ferror_[0] );
889	<	h["errPVy"]->Fill( ferror_[1] );
890	<	h["errPVz"]->Fill( ferror_[2] );
882	>	h_summary->Fill1d("deltax", fres_[0] );
883	>	h_summary->Fill1d("deltay", fres_[1] );
884	>	h_summary->Fill1d("deltaz", fres_[2] );
885	>	h_summary->Fill1d("pullx", fres_[0]/ferror_[0]);
886	>	h_summary->Fill1d("pully", fres_[1]/ferror_[1]);
887	>	h_summary->Fill1d("pullz", fres_[2]/ferror_[2]);
888	>	h_summary->Fill1d("errPVx", ferror_[0] );
889	>	h_summary->Fill1d("errPVy", ferror_[1] );
890	>	h_summary->Fill1d("errPVz", ferror_[2] );
891		pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
892		error(ferror_[0], ferror_[1], ferror_[2]),
893		fntrk_));
#	Line 788 | Line 914 | PVStudy::analyze(const edm::Event& iEven
914		} // End of analyzing res/pull
915		} // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
916		else
917	<	cout<<"WARNING: Cannot find matching pvtxs"<<endl;
917	>	edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
918
919	<	if(verbose_)     cout<<"End analyzing the matchedVertexColl"<<endl;
919	>	edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
920
921		//=======================================================
922		// Fill simulated vertices
#	Line 801 | Line 927 | PVStudy::analyze(const edm::Event& iEven
927		iEvent.getByLabel("generator",evtMC);
928		if (!evtMC.isValid()) {
929		MC=false;
930	<	if(verbose_){
805	<	cout << "no HepMCProduct found"<< endl;
806	<	}
930	>	edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
931		} else {
932	<	if(verbose_){
809	<	cout << "generator HepMCProduct found"<< endl;
810	<	}
932	>	edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
933		MC=true;
934		}
935
#	Line 816 | Line 938 | PVStudy::analyze(const edm::Event& iEven
938		std::vector<simPrimaryVertex> simpv;
939		simpv=getSimPVs(evtMC,"");
940		//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
941	<	h["nsimPV"]->Fill(simpv.size());
941	>	h_gen->Fill1d("nsimPV", simpv.size());
942		nsimPV_ = simpv.size();
943
944		int isimpv = 0;
#	Line 852 | Line 974 | PVStudy::beginJob()
974		// ------------ method called once each job just after ending the event loop  ------------
975		void
976		PVStudy::endJob() {
977	<	if (verbose_) cout << "Analyzing PV info" << endl;
977	>	edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl;
978		// Looping through the datamodes available
979		for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
980		string suffix;
981	<	if(verbose_)  cout<<"datamode = "<< *it<<endl;
981	>	edm::LogInfo("Analysis")<<"[endJob] datamode = "<< *it<<endl;
982		switch(*it) {
983		case 0: suffix = "";
984		break;
#	Line 867 | Line 989 | PVStudy::endJob() {
989		case 3: suffix = "_mct";
990		break;
991		}
992	+	suffix += "_nTrk";
993		if(analyze_) {
994		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
995	+	edm::LogInfo("Analysis") << "[endJob] ntrk = " << ntrk << endl;
996		PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
997	<	if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
998	<	if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
999	<	if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
1000	<	if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
1001	<	if ( ipv.pull_.y() > 0 ) h2["pully"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
1002	<	if ( ipv.pull_.z() > 0 ) h2["pullz"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.z(), 1.);
997	>	if ( ipv.res_.x() > 0 ) h_ana->Fill2d("resx"+suffix, ntrk,ipv.res_.x());
998	>	if ( ipv.res_.y() > 0 ) h_ana->Fill2d("resy"+suffix, ntrk,ipv.res_.y());
999	>	if ( ipv.res_.z() > 0 ) h_ana->Fill2d("resz"+suffix, ntrk,ipv.res_.z());
1000	>	if ( ipv.pull_.x() > 0 ) h_ana->Fill2d("pullx"+suffix, ntrk,ipv.pull_.x());
1001	>	if ( ipv.pull_.y() > 0 ) h_ana->Fill2d("pully"+suffix, ntrk,ipv.pull_.y());
1002	>	if ( ipv.pull_.z() > 0 ) h_ana->Fill2d("pullz"+suffix, ntrk,ipv.pull_.z());
1003		}
1004	<	}
1004	>	}
1005		suffix.clear();
1006		}
1007		if(saventuple_) {
#	Line 905 | Line 1029 | bool PVStudy::matchVertex( const reco::V
1029		int zsigncut)
1030		{
1031		double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
1032	<	if(verbose_)
909	<	cout<<"The matching criteria in z is "<<cut<<endl;
1032	>	edm::LogInfo("Debug")<<"[matchVertex] The matching criteria in z is "<<cut<<endl;
1033		return (fabs(vrec1->z()-vrec2->z())<cut);
1034		}
1035
1036
1037		bool PVStudy::isResonance(const HepMC::GenParticle * p){
1038		double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
1039	<	if(verbose_) cout << "isResonance   " << p->pdg_id() << " " << ctau << endl;
1039	>	edm::LogInfo("Debug") << "[isResonance] isResonance   " << p->pdg_id() << " " << ctau << endl;
1040		return  ctau >0 && ctau <1e-6;
1041		}
1042
#	Line 985 | Line 1108 | void PVStudy::printSimTrks(const Handle<
1108		void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1109		stringstream ss;
1110		ss << ntk;
1111	<	string suffix = ss.str();
1112	<	if(datamode == 0 )  suffix = "_" + suffix;
1113	<	if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1114	<	if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1115	<	if(datamode == 3 )  suffix = "_mct_" + suffix;
1111	>	string suffix;
1112	>	if(datamode == 0 )  suffix = "_" + ss.str();
1113	>	if(datamode == 1 )  suffix = "_spl1_mct_" + ss.str();
1114	>	if(datamode == 2 )  suffix = "_spl2_mct_" + ss.str();
1115	>	if(datamode == 3 )  suffix = "_mct_" + ss.str();
1116
1117		if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
1118		// Fill histograms of res and pull of ntk
1119	<	h["deltax"+suffix]->Fill(r.x());
1120	<	h["deltay"+suffix]->Fill(r.y());
1121	<	h["deltaz"+suffix]->Fill(r.z());
1122	<	h["pullx"+suffix]->Fill(r.x()/er.x());
1123	<	h["pully"+suffix]->Fill(r.y()/er.y());
1124	<	h["pullz"+suffix]->Fill(r.z()/er.z());
1125	<	h["errPVx"+suffix]->Fill(er.x());
1126	<	h["errPVy"+suffix]->Fill(er.y());
1127	<	h["errPVz"+suffix]->Fill(er.z());
1119	>	h_others->Fill1d("deltax"+suffix, r.x());
1120	>	h_others->Fill1d("deltay"+suffix, r.y());
1121	>	h_others->Fill1d("deltaz"+suffix, r.z());
1122	>	h_others->Fill1d("pullx"+suffix, r.x()/er.x());
1123	>	h_others->Fill1d("pully"+suffix, r.y()/er.y());
1124	>	h_others->Fill1d("pullz"+suffix, r.z()/er.z());
1125	>	h_others->Fill1d("errPVx"+suffix, er.x());
1126	>	h_others->Fill1d("errPVy"+suffix, er.y());
1127	>	h_others->Fill1d("errPVz"+suffix, er.z());
1128		}
1129
1130		PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
#	Line 1021 | Line 1144 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1144		for ( int i = 0; i < 6; ++i) {
1145		switch (i) {
1146		case 0:
1147	<	fit(h["deltax"+suffix], fpar);
1147	>	if (!h_others->ReadHisto1D("deltax"+suffix)) break;
1148	>	fit(h_others->ReadHisto1D("deltax"+suffix), fpar);
1149		data[i] = fpar[0]*cm2um;
1150		break;
1151		case 1:
1152	<	fit(h["deltay"+suffix], fpar);
1152	>	if (!h_others->ReadHisto1D("deltay"+suffix)) break;
1153	>	fit(h_others->ReadHisto1D("deltay"+suffix), fpar);
1154		data[i] = fpar[0]*cm2um;
1155		break;
1156		case 2:
1157	<	fit(h["deltaz"+suffix], fpar);
1157	>	if (!h_others->ReadHisto1D("deltaz"+suffix)) break;
1158	>	fit(h_others->ReadHisto1D("deltaz"+suffix), fpar);
1159		data[i] = fpar[0]*cm2um;
1160		break;
1161		case 3:
1162	<	fit(h["pullx"+suffix], fpar);
1162	>	if (!h_others->ReadHisto1D("pullx"+suffix)) break;
1163	>	fit(h_others->ReadHisto1D("pullx"+suffix), fpar);
1164		data[i] = fpar[0];
1165		break;
1166		case 4:
1167	<	fit(h["pully"+suffix], fpar);
1167	>	if (!h_others->ReadHisto1D("pully"+suffix)) break;
1168	>	fit(h_others->ReadHisto1D("pully"+suffix), fpar);
1169		data[i] = fpar[0];
1170		break;
1171		case 5:
1172	<	fit(h["pullz"+suffix], fpar);
1172	>	if (!h_others->ReadHisto1D("pullz"+suffix)) break;
1173	>	fit(h_others->ReadHisto1D("pullz"+suffix), fpar);
1174		data[i] = fpar[0];
1175		break;
1176		}
1177	<	if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1177	>	if (data[i] > 0) edm::LogInfo("Analysis") << "[Analysis] data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1178		}
1179		}
1180		else
#	Line 1060 | Line 1189 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1189		ntk);
1190		}
1191
1192	<	void PVStudy::fit(TH1 *&hdist, double fitPars[]){
1192	>	void PVStudy::fit(TH1 *hdist, double fitPars[]){
1193		TAxis *axis0 = hdist->GetXaxis();
1194		int nbin = axis0->GetLast();
1195		double nOF = hdist->GetBinContent(nbin+1);
#	Line 1069 | Line 1198 | void PVStudy::fit(TH1 *&hdist, double fi
1198		double fitRange = axis0->GetXmax();
1199		double sigMax[2] = {0.,0.};
1200
1201	<	if ( verbose_ ){
1202	<	cout << "Last bin = " << nbin;
1203	<	cout << "; hdist: Overflow Entries = " << nOF;
1204	<	cout << "; Underflow Entries = " << nUF;
1205	<	cout << "; hdist->GetEntries() = " << hdist->GetEntries();
1206	<	cout << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1078	<	}
1201	>	edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
1202	>	<< "; hdist: Overflow Entries = " << nOF
1203	>	<< "; Underflow Entries = " << nUF
1204	>	<< "; hdist->GetEntries() = " << hdist->GetEntries()
1205	>	<< "; fitting range = " << -fitRange << " to " << fitRange << endl;
1206	>
1207		if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit
1208		for (int bi = 0; bi < nbin; bi++) {
1209		if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) {
#	Line 1087 | Line 1215 | void PVStudy::fit(TH1 *&hdist, double fi
1215		if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1216		}
1217		}
1218	<	if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1218	>	//edm::LogInfo("Analysis") << "[Fit] Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1219
1220		TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1221		fgaus->SetParameter(1, 0.);
#	Line 1100 | Line 1228 | void PVStudy::fit(TH1 *&hdist, double fi
1228		s->SetOptStat(1111111);
1229		s->SetOptFit(0111);
1230		gPad->Update();
1231	<	if (verbose_) cout << "got function" << endl;
1231	>	//edm::LogInfo("Analysis") << "[Fit] got function" << endl;
1232		fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1233		}
1234		else
#	Line 1113 | Line 1241 | void PVStudy::fillTrackHisto(const reco:
1241		if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1242		if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1243
1244	<	h["nTrk"+suffix]->Fill(trackColl->size());
1244	>	h_pvtrk->Fill1d("nTrk"+suffix, trackColl->size());
1245		for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1246		itTrack != trackColl->end();
1247		++itTrack) {
1248	<	h["trkPt"+suffix]->Fill(itTrack->pt());
1249	<	h["trkDxy"+suffix]->Fill(itTrack->dxy());
1250	<	h["trkDz"+suffix]->Fill(itTrack->dz());
1251	<	h["trkEta"+suffix]->Fill(itTrack->eta());
1252	<	h["trkPhi"+suffix]->Fill(itTrack->phi());
1248	>	h_pvtrk->Fill1d("trkPt"+suffix, itTrack->pt());
1249	>	h_pvtrk->Fill1d("trkDxy"+suffix, itTrack->dxy());
1250	>	h_pvtrk->Fill1d("trkDz"+suffix, itTrack->dz());
1251	>	h_pvtrk->Fill1d("trkEta"+suffix, itTrack->eta());
1252	>	h_pvtrk->Fill1d("trkPhi"+suffix, itTrack->phi());
1253		}
1254		}
1255
#	Line 1135 | Line 1263 | void PVStudy::fillTrackHistoInPV(const r
1263		v!=vertexColl->end(); ++v, ++ivtx) {
1264		if(fillNtuple) {
1265		if(datatype == 0) {
1266	<	nTrkPV_[ivtx] = v->tracksSize();
1266	>	nTrkPV_[ivtx] = v->tracksSize();
1267	>	sumptsq_[ivtx] = sumPtSquared(*v);
1268	>	isValid_[ivtx] = v->isValid();
1269	>	isFake_[ivtx] = v->isFake();
1270		recx_[ivtx] = v->x();
1271		recy_[ivtx] = v->y();
1272		recz_[ivtx] = v->z();
#	Line 1144 | Line 1275 | void PVStudy::fillTrackHistoInPV(const r
1275		recz_err_[ivtx] = v->zError();
1276		}
1277		if(datatype == 1) {
1278	<	nTrkPV1_spl_[ivtx] = v->tracksSize();
1278	>	nTrkPV1_spl_[ivtx] = v->tracksSize();
1279	>	sumptsq1_spl_[ivtx] = sumPtSquared(*v);
1280	>	isValid1_spl_[ivtx] = v->isValid();
1281	>	isFake1_spl_[ivtx] = v->isFake();
1282		recx1_spl_[ivtx] = v->x();
1283		recy1_spl_[ivtx] = v->y();
1284		recz1_spl_[ivtx] = v->z();
1285		recx1_err_spl_[ivtx] = v->xError();
1286		recy1_err_spl_[ivtx] = v->yError();
1287		recz1_err_spl_[ivtx] = v->zError();
1288	+
1289		}
1290		if(datatype == 2) {
1291	<	nTrkPV2_spl_[ivtx] = v->tracksSize();
1291	>	nTrkPV2_spl_[ivtx] = v->tracksSize();
1292	>	sumptsq2_spl_[ivtx] = sumPtSquared(*v);
1293	>	isValid2_spl_[ivtx] = v->isValid();
1294	>	isFake2_spl_[ivtx] = v->isFake();
1295		recx2_spl_[ivtx] = v->x();
1296		recy2_spl_[ivtx] = v->y();
1297		recz2_spl_[ivtx] = v->z();
#	Line 1161 | Line 1299 | void PVStudy::fillTrackHistoInPV(const r
1299		recy2_err_spl_[ivtx] = v->yError();
1300		recz2_err_spl_[ivtx] = v->zError();
1301		}
1302	+	if(datatype == 4) { // for pixelVertices
1303	+	nTrkPV_pxlpvtx_[ivtx] = v->tracksSize();
1304	+	sumptsq_pxlpvtx_[ivtx] = sumPtSquared(*v);
1305	+	isValid_pxlpvtx_[ivtx] = v->isValid();
1306	+	isFake_pxlpvtx_[ivtx] = v->isFake();
1307	+	recx_pxlpvtx_[ivtx] = v->x();
1308	+	recy_pxlpvtx_[ivtx] = v->y();
1309	+	recz_pxlpvtx_[ivtx] = v->z();
1310	+	recx_err_pxlpvtx_[ivtx] = v->xError();
1311	+	recy_err_pxlpvtx_[ivtx] = v->yError();
1312	+	recz_err_pxlpvtx_[ivtx] = v->zError();
1313	+	}
1314		}
1315	<
1316	<	if(fillHisto) {
1317	<	h["nTrkPV"+suffix]->Fill(v->tracksSize());
1318	<	try {
1319	<	for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
1320	<	t!=v->tracks_end(); t++) {
1321	<	// illegal charge
1322	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1323	<	//    h["trkPtPV"]->Fill(0.);
1324	<	}
1325	<	else {
1326	<	h["trkPtPV"+suffix]->Fill((**t).pt());
1327	<	h["trkDxyPV"+suffix]->Fill((**t).dxy());
1328	<	h["trkDzPV"+suffix]->Fill((**t).dz());
1329	<	h["trkEtaPV"+suffix]->Fill((**t).eta());
1330	<	h["trkPhiPV"+suffix]->Fill((**t).phi());
1331	<	}
1315	>	}
1316	>	// For histogram only fill the leading pvtx parameters
1317	>	if(fillHisto) {
1318	>	h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1319	>	try {
1320	>	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1321	>	t!=vertexColl->begin()->tracks_end(); t++) {
1322	>	// illegal charge
1323	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1324	>	//      h_pvtrk->Fill1d("trkPtPV", 0.);
1325	>	}
1326	>	else {
1327	>	h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1328	>	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1329	>	h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1330	>	h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1331	>	h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1332		}
1183	–	}
1184	–	catch (...) {
1185	–	// exception thrown when trying to use linked track
1186	–	//        h["trkPtPV"]->Fill(0.);
1333		}
1334		}
1335	+	catch (...) {
1336	+	// exception thrown when trying to use linked track
1337	+	//          h_pvtrk->Fill1d("trkPtPV", 0.);
1338	+	}
1339		}
1340	+
1341		}
1342
1343		void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
#	Line 1215 | Line 1366 | void PVStudy::fillMCHisto(std::vector<si
1366		int nrectrk = vrec->tracksSize();
1367		vsim->recVtx=NULL;
1368
1369	<	if(verbose_){
1370	<	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1371	<	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1372	<	}
1222	<	if ( matchVertex(*vsim,*vrec) ) {
1369	>	edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1370	>	edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1371	>
1372	>	if ( matchVertex(*vsim,*vrec) && vrec->isValid() && !vrec->isFake() ) {
1373		vsim->recVtx=&(*vrec);
1374
1375		// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
#	Line 1227 | Line 1377 | void PVStudy::fillMCHisto(std::vector<si
1377		//|| (!vsim->recVtx) )
1378		//vsim->recVtx=&(*vrec);
1379		//}
1380	<	if(verbose_)
1231	<	cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1380	>	edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1381
1382		double fres_mct[3];
1383		double ferror_mct[3];
#	Line 1240 | Line 1389 | void PVStudy::fillMCHisto(std::vector<si
1389		ferror_mct[1] = vsim->recVtx->yError();
1390		ferror_mct[2] = vsim->recVtx->zError();
1391
1392	<	h["deltax"+suffix]->Fill( fres_mct[0] );
1393	<	h["deltay"+suffix]->Fill( fres_mct[1] );
1394	<	h["deltaz"+suffix]->Fill( fres_mct[2] );
1395	<	h["pullx"+suffix]->Fill( fres_mct[0]/ferror_mct[0] );
1396	<	h["pully"+suffix]->Fill( fres_mct[1]/ferror_mct[1] );
1397	<	h["pullz"+suffix]->Fill( fres_mct[2]/ferror_mct[2] );
1398	<	h["errPVx"+suffix]->Fill( ferror_mct[0] );
1399	<	h["errPVy"+suffix]->Fill( ferror_mct[1] );
1400	<	h["errPVz"+suffix]->Fill( ferror_mct[2] );
1392	>	h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1393	>	h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1394	>	h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1395	>	h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1396	>	h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1397	>	h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1398	>	h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1399	>	h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1400	>	h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1401		pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1402		error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1403		nrectrk));
#	Line 1293 | Line 1442 | void PVStudy::fillMCHisto(std::vector<si
1442		} // End of  if(saventuple_) {
1443		} //  if ( matchVertex(*vsim,*vrec) ) {
1444		else {  // no rec vertex found for this simvertex
1445	<	if(verbose_)
1297	<	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1445	>	edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1446		}
1447		}
1448		}
#	Line 1305 | Line 1453 | void PVStudy::SetVarToZero() {
1453		fntrk_ = 0;
1454		//pvtx position (x,y,z) residual and error
1455		for(int i = 0; i<3;i++)    {
1456	<	fres_[i] = 0;
1457	<	ferror_[i] = 0;
1456	>	fres_[i] = 0;
1457	>	ferror_[i] = 0;
1458		}
1459
1460		//Yanyan's variables
1461	<	// Number of reconstructed vertices  // for GeneralTracks
1461	>	// Number of reconstructed vertices
1462		nrecPV_ = 0;
1463		nrecPV1_spl_ = 0;
1464		nrecPV2_spl_ = 0;
1465		nrecPV_twovtx_ = 0;
1466		nsimPV_ = 0;
1467	+
1468		nrecPV_mct_ = 0;
1469		nrecPV_spl1_mct_ = 0;
1470		nrecPV_spl2_mct_ = 0;
1471
1472	+	// Mininum separation between the secondary pvtxes and leading pvtx
1473	+	min_zsep_ = 9999.0;
1474	+	min_ntrksep_ = 9999.0;
1475	+	min_sumpt2sep_ = -9999.0;
1476	+
1477	+
1478		for (int i = 0; i < nMaxPVs_; i++) {
1479	<	nTrkPV_[i] = 0; // Number of tracks in the pvtx
1479	>	// recoVertices with all tracks
1480	>	nTrkPV_[i] = 0; // Number of tracks in the pvtx
1481	>	sumptsq_[i] = 0;
1482	>	isValid_[i] = -1;
1483	>	isFake_[i] = -1;
1484		recx_[i] = 0;
1485		recy_[i] = 0;
1486		recz_[i] = 0;
1487		recx_err_[i] = 0;
1488		recy_err_[i] = 0;
1489		recz_err_[i] = 0;
1490	<
1491	<	nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx
1490	>
1491	>	// recoVertices with splitTrack1
1492	>	nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx
1493	>	sumptsq1_spl_[i] = 0;
1494	>	isValid1_spl_[i] = -1;
1495	>	isFake1_spl_[i] = -1;
1496		recx1_spl_[i] = 0;
1497		recy1_spl_[i] = 0;
1498		recz1_spl_[i] = 0;
1499		recx1_err_spl_[i] = 0;
1500		recy1_err_spl_[i] = 0;
1501		recz1_err_spl_[i] = 0;
1339	–
1340	–	nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx
1502
1503	+	// recoVertices with splitTrack2
1504	+	nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx
1505	+	sumptsq2_spl_[i] = 0;
1506	+	isValid2_spl_[i] = -1;
1507	+	isFake2_spl_[i] = -1;
1508		recx2_spl_[i] = 0;
1509		recy2_spl_[i] = 0;
1510		recz2_spl_[i] = 0;
#	Line 1346 | Line 1512 | void PVStudy::SetVarToZero() {
1512		recy2_err_spl_[i] = 0;
1513		recz2_err_spl_[i] = 0;
1514
1515	+	//pixelVertices
1516	+	nTrkPV_pxlpvtx_[i] = 0; // Number of tracks in the pvtx
1517	+	sumptsq_pxlpvtx_[i] = 0;
1518	+	isValid_pxlpvtx_[i] = -1;
1519	+	isFake_pxlpvtx_[i] = -1;
1520	+	recx_pxlpvtx_[i] = 0;
1521	+	recy_pxlpvtx_[i] = 0;
1522	+	recz_pxlpvtx_[i] = 0;
1523	+	recx_err_pxlpvtx_[i] = 0;
1524	+	recy_err_pxlpvtx_[i] = 0;
1525	+	recz_err_pxlpvtx_[i] = 0;
1526	+
1527	+	// matched two-vertices
1528		nTrkPV1_spl_twovtx_[i] = 0;
1529		nTrkPV2_spl_twovtx_[i] = 0;
1530		nTrkPV_twovtx_[i] = 0;
#	Line 1361 | Line 1540 | void PVStudy::SetVarToZero() {
1540		pully_twovtx_[i] = 0;
1541		pullz_twovtx_[i] = 0;
1542
1543	+	//simpv
1544		nsimTrkPV_[i] = 0;
1545		simx_[i] = 0;
1546		simy_[i] = 0;
1547		simz_[i] = 0;
1548		simptsq_[i] = 0;
1549
1550	+	// residual and pulls with mc truth required
1551		deltax_mct_[i] = 0;
1552		deltay_mct_[i] = 0;
1553		deltaz_mct_[i] = 0;
#	Line 1400 | Line 1581 | double PVStudy::sumPtSquared(const reco:
1581		}
1582		return sum;
1583		}
1584	+

Diff Legend

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