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.11 by yygao, Wed Nov 11 19:28:50 2009 UTC vs.
Revision 1.16 by jengbou, Thu Dec 10 17:57:31 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	<	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");
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 217 | Line 218 | PVStudy::PVStudy(const edm::ParameterSet
218		}
219		}
220
220	–	edm::Service<TFileService> fs;
221	–	TFileDirectory subDir = fs->mkdir( "Summary" );
222	–	TFileDirectory subDir1 = fs->mkdir( "Others" );
223	–	//   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
224	–
221		setRootStyle();
226	–
227	–	//Book histograms for track and pvtx parameters for each splitted pvtx
228	–	//i=0 : x (as in unsplit track collection)
229	–	//i=1 : x1_spl_
230	–	//i=2 : x2_spl_
222
223	<
224	<	for(int i=0;i<3;i++)
225	<	{
226	<	string suffix;
227	<	if(i == 0) suffix = "";
228	<	else {
229	<	stringstream ss;
230	<	ss << i;
231	<	suffix =ss.str()+"_spl";
232	<	}
233	<	h["nTrk"+suffix]   = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
234	<	h["trkPt"+suffix]  = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
235	<	h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#eta of rec tracks "+suffix),100,-3,3);
236	<	h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
237	<	h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-0.5,0.5);
238	<	h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
239	<
240	<	h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
241	<	h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
242	<	h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#eta of rec tracks in PV"+suffix),100,-3,3);
243	<	h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
244	<	h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);
245	<	h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
246	<	h["nrecPV"+suffix]   = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);
247	<	suffix.clear();
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	<
253	<	// Histograms on comparing the multi-vertices
254	<	// Difference in reconstructed vtx position
255	<	h["min_xsep"]   = subDir.make<TH1D>("min_xsep", "min x diff of primary and secondary pvtx",300,0,0.1);
256	<	h["min_xsep_sign"]   = subDir.make<TH1D>("min_xsep_sign", "min x diff in signf of primary and secondary pvtx",300,0,5);
257	<	h["min_ysep"]   = subDir.make<TH1D>("min_ysep", "min y diff of primary and secondary pvtx",300,0,0.1);
258	<	h["min_ysep_sign"]   = subDir.make<TH1D>("min_ysep_sign", "min y diff in signf of primary and secondary pvtx",300,0,5);
268	<	h["min_zsep"]   = subDir.make<TH1D>("min_zsep", "min z diff of primary and secondary pvtx",300,0,5);
269	<	h["min_zsep_sign"]   = subDir.make<TH1D>("min_zsep_sign", "min z diff in signf of primary and secondary pvtx",300,0,200);
270	<	// Difference in reconstructed vtx position
271	<	h["min_ntrksep"]   = subDir.make<TH1D>("min_ntrksep", "min nTrk diff of primary and secondary pvtx",201,-50.5,150.5);
272	<	h["min_sumpt2sep"]   = subDir.make<TH1D>("min_sumpt2sep", "min sumpt2 diff of primary and secondary pvtx",300,0,10000);
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 285 | 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);
289	<	h["deltay"+suffix] = subDir.make<TH1D>(TString("deltay"+suffix), TString("y-residual pvtx"+suffix),800,-0.04,0.04);
290	<	h["deltaz"+suffix] = subDir.make<TH1D>(TString("deltaz"+suffix), TString("z-residual pvtx"+suffix),800,-0.04,0.04);
291	<	h["pullx"+suffix]  = subDir.make<TH1D>(TString("pullx"+suffix), TString("x-pull pvtx"+suffix),800,-10,10);
292	<	h["pully"+suffix]  = subDir.make<TH1D>(TString("pully"+suffix), TString("y-pull pvtx"+suffix),800,-10,10);
293	<	h["pullz"+suffix]  = subDir.make<TH1D>(TString("pullz"+suffix), TString("z-pull pvtx"+suffix),800,-10,10);
294	<	h["errPVx"+suffix] = subDir.make<TH1D>(TString("errPVx"+suffix), TString("X"+suffix+" vertex error"),100,0.,0.02);
295	<	h["errPVy"+suffix] = subDir.make<TH1D>(TString("errPVy"+suffix), TString("Y"+suffix+" vertex error"),100,0.,0.02);
296	<	h["errPVz"+suffix] = subDir.make<TH1D>(TString("errPVz"+suffix), TString("Z"+suffix+" vertex error"),100,0.,0.02);
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");
305	<	h["deltay"+suffix2]  = subDir1.make<TH1D>(TString("deltay"+suffix2), TString("y-residual of pvtx"+suffix),100,-0.02,0.02);
306	<	h["deltay"+suffix2]->GetXaxis()->SetTitle("cm");
307	<	h["deltaz"+suffix2]  = subDir1.make<TH1D>(TString("deltaz"+suffix2), TString("z-residual of pvtx"+suffix),100,-0.02,0.02);
308	<	h["deltaz"+suffix2]->GetXaxis()->SetTitle("cm");
309	<	h["pullx"+suffix2] = subDir1.make<TH1D>(TString("pullx"+suffix2), TString("x-pull of pvtx"+suffix),100,-10.,10.);
310	<	h["pully"+suffix2] = subDir1.make<TH1D>(TString("pully"+suffix2), TString("y-pull of pvtx"+suffix),100,-10.,10.);
311	<	h["pullz"+suffix2] = subDir1.make<TH1D>(TString("pullz"+suffix2), TString("z-pull of pvtx"+suffix),100,-10.,10.);
312	<	h["errPVx"+suffix2] = subDir1.make<TH1D>(TString("errPVx"+suffix2), TString("X"+suffix+" vertex error"),100,0.,0.02);
313	<	h["errPVy"+suffix2] = subDir1.make<TH1D>(TString("errPVy"+suffix2), TString("Y"+suffix+" vertex error"),100,0.,0.02);
314	<	h["errPVz"+suffix2] = subDir1.make<TH1D>(TString("errPVz"+suffix2), TString("Z"+suffix+" vertex error"),100,0.,0.02);
315	<	suffix2.clear();
316	<	} // End of   for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
317	<
280	>	h_others->Init("others",suffix,ss.str());
281	>	}
282	>
283		// Book residual and pull histograms only when analyzeOntheFly is enabled
284	<	if(analyze_) {
320	<	h2["resx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resx"+suffix+"_nTrk"), TString("x-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
321	<	h2["resx"+suffix+"_nTrk"]->SetMarkerStyle(21);
322	<	h2["resx"+suffix+"_nTrk"]->SetMarkerColor(4);
323	<	h2["resx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
324	<	h2["resx"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
325	<	h2["resy"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resy"+suffix+"_nTrk"), TString("y-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
326	<	h2["resy"+suffix+"_nTrk"]->SetMarkerStyle(21);
327	<	h2["resy"+suffix+"_nTrk"]->SetMarkerColor(4);
328	<	h2["resy"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
329	<	h2["resy"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
330	<	h2["resz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resz"+suffix+"_nTrk"), TString("z-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
331	<	h2["resz"+suffix+"_nTrk"]->SetMarkerStyle(21);
332	<	h2["resz"+suffix+"_nTrk"]->SetMarkerColor(4);
333	<	h2["resz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
334	<	h2["resz"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
335	<	h2["pullx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullx"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
336	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerStyle(21);
337	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerColor(4);
338	<	h2["pullx"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
339	<	h2["pullx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
340	<	h2["pully"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pully"+suffix+"_nTrk"), TString("y-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
341	<	h2["pully"+suffix+"_nTrk"]->SetMarkerStyle(21);
342	<	h2["pully"+suffix+"_nTrk"]->SetMarkerColor(4);
343	<	h2["pully"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
344	<	h2["pully"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
345	<	h2["pullz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullz"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
346	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerStyle(21);
347	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerColor(4);
348	<	h2["pullz"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
349	<	h2["pullz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
350	<	} // End of   if(analyze_) {
284	>	if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
285		suffix.clear();
286		} // End of Book histograms sensitive to data/mc
353	–
354	–	// Book histograms about pixelVertices
355	–	h["trkdz_pxlpvtxdz"]   = subDir.make<TH1D>("trkdz_pxlpvtxdz", "(Track dz - pixelpvtx dz) in cm",300,-0.5,0.5);
356	–	h["trkdz_pxlpvtxdz_pxlpvtxdzerr"]   = subDir.make<TH1D>("trkdz_pxlpvtxdz_pxlpvtxdzerr", "|Track dz - pixelpvtx dz| / pxlpvtxdzErr",300,0,100);
357	–	h["trkdz_pxlpvtxdz_trkdzerr"]   = subDir.make<TH1D>("trkdz_pxlpvtxdz_trkdzerr", "|Track dz - pixelpvtx dz| / trkdzErr",300,0,50);
358	–	h["trkdzErr_pxlpvtx"]   = subDir.make<TH1D>("trkdzErr_pxlpvtxdz", "Track dzErr of leading pixelpvtx ",300,0,0.5);
359	–	h["trkdzErr_pvtx"]   = subDir.make<TH1D>("trkdzErr_pvtx", "Track dzErr of the leading pvtx ",300,0,0.5);
360	–	h["dzErr_pxlpvtx"]   = subDir.make<TH1D>("dzErr_pxlpvtx", "zError of the leading pvtx ",300,0,0.5);
361	–	// Compare offlinePrimaryVertices with pixelVertices
362	–	h["nrecPV_minus_nrecPxlPV"] =  subDir.make<TH1D>("nrecPV_minus_nrecPxlPV", "nrecPV_minus_nrecPxlPV",21,-10.5,10.5);
363	–	h["recxPV_minus_recxPxlPV"] =  subDir.make<TH1D>("recxPV_minus_recxPxlPV", "recxPV_minus_recxPxlPV",300,-0.02,0.02);
364	–	h["recyPV_minus_recyPxlPV"] =  subDir.make<TH1D>("recyPV_minus_recyPxlPV", "recyPV_minus_recyPxlPV",300,-0.02,0.02);
365	–	h["reczPV_minus_reczPxlPV"] =  subDir.make<TH1D>("reczPV_minus_reczPxlPV", "reczPV_minus_reczPxlPV",300,-0.1,0.1);
287
367	–	// Book MC only plots
368	–	if (!realData_) {
369	–	h["genPart_T"]      = subDir.make<TH1D>("genPart_T","t component of gen particles",300,-0.5,0.5);
370	–	h["genPart_T"]->GetXaxis()->SetTitle("t (nanosecond)");
371	–	h["genPart_cT"]     = subDir.make<TH1D>("genPart_cT","ct component of gen particles",300,-150.,150.);
372	–	h["genPart_cT"]->GetXaxis()->SetTitle("ct (mm)");
373	–	h["nsimPV"]         = subDir.make<TH1D>("nsimPV","Num of sim PV",51,-0.5,50.5);
374	–	}
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 419 | 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;
426	<	cout << "number of vertices " << evt->vertices_size() << endl;
427	<	}
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 435 | 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 449 | 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 480 | 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 508 | 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 542 | 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 591 | Line 520 | PVStudy::analyze(const edm::Event& iEven
520		using namespace edm;
521		using namespace reco;
522
523	<	if (verbose_)
595	<	cout<<"PVStudy::analyze():"<<endl;
523	>	edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
524		//=======================================================
525		// Initialize Root-tuple variables if needed
526		//=======================================================
#	Line 610 | 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
#	Line 621 | 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
#	Line 632 | 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
566
567	<	//=======================================================
567	>	//=======================================================
568		// Fill trackparameters of the input tracks to pvtx fitter
569		//=======================================================
570	<	if(verbose_)
643	<	cout<<"Start filling track parameters of the input tracks to pvtx fitter."<<endl;
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	<	if(verbose_)
650	<	cout<<"End filling track parameters of the input tracks to pvtx fitter."<<endl;
576	>	edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
577
578		//=======================================================
579		// PVTX accessors
#	Line 660 | 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
#	Line 671 | 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;
600	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
601		}
602
603		//splitVertexCollection2
#	Line 682 | 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<<"End accessing the track and vertex collections"<<endl;
614	>	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track and vertex collections"<<endl;
615
616		//=======================================================
617		// MC simvtx accessor
#	Line 704 | 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
636		//=======================================================
#	Line 717 | Line 643 | PVStudy::analyze(const edm::Event& iEven
643		if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
644		pixelVertexColl = pixelVertexCollectionHandle.product();
645		} else {
646	<	cout << "pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
646	>	edm::LogInfo("Debug") << "[PVStudy] pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
647		}
648
649		//=======================================================
#	Line 727 | Line 653 | PVStudy::analyze(const edm::Event& iEven
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["dzErr_pxlpvtx"]->Fill( pixelVertexColl->begin()->zError());
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["trkPtPV"]->Fill(0.);
662	>	//        h_pvtrk->Fill1d("trkPtPV", 0.);
663		}
664		else
665	<	h["trkdzErr_pxlpvtx"]->Fill((**t).dzError());
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["trkPtPV"]->Fill(0.);
672	>	//        h_pvtrk->Fill1d("trkPtPV", 0.);
673		}
674		else {
675		if(pixelVertexColl->size()>0) {
676	<	h["trkdz_pxlpvtxdz"]->Fill((**t).dz() -  pixelVertexColl->begin()->z());
677	<	h["trkdz_pxlpvtxdz_pxlpvtxdzerr"]->Fill(fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
678	<	h["trkdz_pxlpvtxdz_trkdzerr"]->Fill(fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
679	<	h["trkdzErr_pvtx"]->Fill((**t).dzError());
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		}
#	Line 760 | Line 686 | PVStudy::analyze(const edm::Event& iEven
686		// Fill number of reconstructed vertices
687		//=======================================================
688
689	<	if(verbose_)  {
690	<	cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
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);
766	–	cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
695		printRecVtxs(splitVertexCollection1Handle);
768	–	cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
696		printRecVtxs(splitVertexCollection2Handle);
770	–	cout<<"Start filling pvtx parameters."<<endl;
697		}
772	–
698		nrecPV_ = int (vertexColl->size());
699		nrecPV1_spl_ = int (splitVertexColl1->size());
700		nrecPV2_spl_ = int (splitVertexColl2->size());
701
702	<	h["nrecPV"]->Fill(nrecPV_);
703	<	h["nrecPV1_spl"]->Fill(nrecPV1_spl_);
704	<	h["nrecPV2_spl"]->Fill(nrecPV2_spl_);
705	<	h["nrecPVDiff"]->Fill(double(nrecPV1_spl_)-double(nrecPV2_spl_));
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		//=================================================================
708		// Fill track parameter ntuple/hist for tracks used in recoVertices
#	Line 799 | Line 724 | PVStudy::analyze(const edm::Event& iEven
724		//=======================================================
725		if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
726		&& (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {
727	<	h["nrecPV_minus_nrecPxlPV"]->Fill( double (nrecPV_ - nrecPV_pxlpvtx_));
727	>	h_pixvtx->Fill1d("nrecPV_minus_nrecPxlPV", double (nrecPV_ - nrecPV_pxlpvtx_));
728		// difference in reconstructed position of the leading pvtx
729	<	//cout<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;
730	<	//cout<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
731	<	h["recxPV_minus_recxPxlPV"]->Fill (recx_[0] - recx_pxlpvtx_[0]);
732	<	h["recyPV_minus_recyPxlPV"]->Fill (recy_[0] - recy_pxlpvtx_[0]);
733	<	h["reczPV_minus_reczPxlPV"]->Fill (recz_[0] - recz_pxlpvtx_[0]);
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
#	Line 825 | Line 750 | PVStudy::analyze(const edm::Event& iEven
750		double min_ntrksep = 9999.0;
751		double min_sumpt2sep = 9999999.0;
752
753	<	if(verbose_) cout<<"leading pvtx z = "<< vertexColl->begin()->z() <<endl;
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 ;
#	Line 860 | Line 785 | PVStudy::analyze(const edm::Event& iEven
785		min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
786		}
787		}
788	<	h["min_xsep"]->Fill(min_xsep);
789	<	h["min_ysep"]->Fill(min_ysep);
790	<	h["min_zsep"]->Fill(min_zsep);
791	<	h["min_xsep_sign"]->Fill(min_xsep_sign);
792	<	h["min_ysep_sign"]->Fill(min_ysep_sign);
793	<	h["min_zsep_sign"]->Fill(min_zsep_sign);
794	<	h["min_ntrksep"]->Fill(min_ntrksep);
795	<	h["min_sumpt2sep"]->Fill(min_sumpt2sep);
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;
#	Line 876 | Line 801 | PVStudy::analyze(const edm::Event& iEven
801
802		} // End of  if(vertexColl->size()>1) {
803
804	<	if(verbose_)
880	<	cout<<"End filling pvtx parameters."<<endl;
804	>	edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
805
806		//=======================================================
807		//           PrimaryVertex Matching
#	Line 887 | Line 811 | PVStudy::analyze(const edm::Event& iEven
811		// since vertexColl are sorted in decreasing order of sum(pT)
812		//=======================================================
813
814	<	if(verbose_)   cout<<"PVStudy::analyze(): matching pvtxs "<<endl;
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());
#	Line 906 | Line 830 | PVStudy::analyze(const edm::Event& iEven
830		if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
831		reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
832		if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
833	<	if(verbose_) {
834	<	cout<<"Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
835	<	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
836	<	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
913	<	}
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	<	if(verbose_) {
839	<	cout<<"The two splitted vertices match in Z. "<<endl;
917	<	}
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["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
844	<	h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
843	>	h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
844	>	h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
845		if(fabs(ntrkreldiff)<ntrkdiffcut_) {
846	<	if(verbose_) {
847	<	cout<<"(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
926	<	}
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	<	cout<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
854	>	edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
855		}
856		else {
857	<	cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
858	<	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
859	<	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
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		}
861		}
862		}
863	<	if(verbose_)     cout<<"PVStudy::analyze(): End matching pvtxs"<<endl;
863	>	edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
864
865		//=======================================================
866		//  Analyze matchedVertexColls
867		//=======================================================
868	<	if(verbose_) {
948	<	cout<<"Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
949	<	}
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) {
#	Line 968 | Line 887 | PVStudy::analyze(const edm::Event& iEven
887		ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
888		fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
889
890	<	h["deltax"]->Fill( fres_[0] );
891	<	h["deltay"]->Fill( fres_[1] );
892	<	h["deltaz"]->Fill( fres_[2] );
893	<	h["pullx"]->Fill( fres_[0]/ferror_[0]);
894	<	h["pully"]->Fill( fres_[1]/ferror_[1]);
895	<	h["pullz"]->Fill( fres_[2]/ferror_[2]);
896	<	h["errPVx"]->Fill( ferror_[0] );
897	<	h["errPVy"]->Fill( ferror_[1] );
898	<	h["errPVz"]->Fill( ferror_[2] );
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_));
#	Line 1003 | Line 922 | PVStudy::analyze(const edm::Event& iEven
922		} // End of analyzing res/pull
923		} // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
924		else
925	<	cout<<"WARNING: Cannot find matching pvtxs"<<endl;
925	>	edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
926
927	<	if(verbose_)     cout<<"End analyzing the matchedVertexColl"<<endl;
927	>	edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
928
929		//=======================================================
930		// Fill simulated vertices
#	Line 1016 | Line 935 | PVStudy::analyze(const edm::Event& iEven
935		iEvent.getByLabel("generator",evtMC);
936		if (!evtMC.isValid()) {
937		MC=false;
938	<	if(verbose_){
1020	<	cout << "no HepMCProduct found"<< endl;
1021	<	}
938	>	edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
939		} else {
940	<	if(verbose_){
1024	<	cout << "generator HepMCProduct found"<< endl;
1025	<	}
940	>	edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
941		MC=true;
942		}
943
#	Line 1031 | 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());
949	>	h_gen->Fill1d("nsimPV", simpv.size());
950		nsimPV_ = simpv.size();
951
952		int isimpv = 0;
#	Line 1067 | 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 1082 | Line 997 | PVStudy::endJob() {
997		case 3: suffix = "_mct";
998		break;
999		}
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 ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
1006	<	if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
1007	<	if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
1008	<	if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
1009	<	if ( ipv.pull_.y() > 0 ) h2["pully"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
1010	<	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	<	}
1012	>	}
1013		suffix.clear();
1014		}
1015		if(saventuple_) {
#	Line 1120 | Line 1037 | bool PVStudy::matchVertex( const reco::V
1037		int zsigncut)
1038		{
1039		double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
1040	<	if(verbose_)
1124	<	cout<<"The matching criteria in z is "<<cut<<endl;
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 1200 | 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 1236 | 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 1275 | 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 1284 | 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;
1293	<	}
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 1302 | 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.);
1230		fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
1231		fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
1232		fgaus->SetLineColor(4);
1233	<	hdist->Fit(fgaus,"Q");
1233	>	hdist->Fit(fgaus,"QLRM");
1234		gPad->Update();
1235		TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
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
#	Line 1328 | Line 1249 | void PVStudy::fillTrackHisto(const reco:
1249		if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1250		if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1251
1252	<	h["nTrk"+suffix]->Fill(trackColl->size());
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["trkPt"+suffix]->Fill(itTrack->pt());
1257	<	h["trkDxy"+suffix]->Fill(itTrack->dxy());
1258	<	h["trkDz"+suffix]->Fill(itTrack->dz());
1259	<	h["trkEta"+suffix]->Fill(itTrack->eta());
1260	<	h["trkPhi"+suffix]->Fill(itTrack->phi());
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
#	Line 1360 | Line 1281 | void PVStudy::fillTrackHistoInPV(const r
1281		recx_err_[ivtx] = v->xError();
1282		recy_err_[ivtx] = v->yError();
1283		recz_err_[ivtx] = v->zError();
1363	–
1364	–
1284		}
1285		if(datatype == 1) {
1286		nTrkPV1_spl_[ivtx] = v->tracksSize();
#	Line 1404 | Line 1323 | void PVStudy::fillTrackHistoInPV(const r
1323		}
1324		// For histogram only fill the leading pvtx parameters
1325		if(fillHisto) {
1326	<	h["nTrkPV"+suffix]->Fill(vertexColl->begin()->tracksSize());
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["trkPtPV"]->Fill(0.);
1332	>	//      h_pvtrk->Fill1d("trkPtPV", 0.);
1333		}
1334		else {
1335	<	h["trkPtPV"+suffix]->Fill((**t).pt());
1336	<	h["trkDxyPV"+suffix]->Fill((**t).dxy());
1337	<	h["trkDzPV"+suffix]->Fill((**t).dz());
1338	<	h["trkEtaPV"+suffix]->Fill((**t).eta());
1339	<	h["trkPhiPV"+suffix]->Fill((**t).phi());
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["trkPtPV"]->Fill(0.);
1345	>	//          h_pvtrk->Fill1d("trkPtPV", 0.);
1346		}
1347		}
1348
#	Line 1455 | Line 1374 | void PVStudy::fillMCHisto(std::vector<si
1374		int nrectrk = vrec->tracksSize();
1375		vsim->recVtx=NULL;
1376
1377	<	if(verbose_){
1378	<	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1379	<	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1380	<	}
1462	<	if ( matchVertex(*vsim,*vrec) ) {
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
#	Line 1467 | Line 1385 | void PVStudy::fillMCHisto(std::vector<si
1385		//|| (!vsim->recVtx) )
1386		//vsim->recVtx=&(*vrec);
1387		//}
1388	<	if(verbose_)
1471	<	cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
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];
#	Line 1480 | Line 1397 | void PVStudy::fillMCHisto(std::vector<si
1397		ferror_mct[1] = vsim->recVtx->yError();
1398		ferror_mct[2] = vsim->recVtx->zError();
1399
1400	<	h["deltax"+suffix]->Fill( fres_mct[0] );
1401	<	h["deltay"+suffix]->Fill( fres_mct[1] );
1402	<	h["deltaz"+suffix]->Fill( fres_mct[2] );
1403	<	h["pullx"+suffix]->Fill( fres_mct[0]/ferror_mct[0] );
1404	<	h["pully"+suffix]->Fill( fres_mct[1]/ferror_mct[1] );
1405	<	h["pullz"+suffix]->Fill( fres_mct[2]/ferror_mct[2] );
1406	<	h["errPVx"+suffix]->Fill( ferror_mct[0] );
1407	<	h["errPVy"+suffix]->Fill( ferror_mct[1] );
1408	<	h["errPVz"+suffix]->Fill( ferror_mct[2] );
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));
#	Line 1533 | Line 1450 | void PVStudy::fillMCHisto(std::vector<si
1450		} // End of  if(saventuple_) {
1451		} //  if ( matchVertex(*vsim,*vrec) ) {
1452		else {  // no rec vertex found for this simvertex
1453	<	if(verbose_)
1537	<	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1453	>	edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1454		}
1455		}
1456		}
#	Line 1545 | Line 1461 | void PVStudy::SetVarToZero() {
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;
1464	>	fres_[i] = 0;
1465	>	ferror_[i] = 0;
1466		}
1467
1468		//Yanyan's variables

Diff Legend

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