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.12 by yygao, Fri Nov 27 12:31:05 2009 UTC vs.
Revision 1.18 by jengbou, Fri Jan 29 21:36:22 2010 UTC

#	Line 25 | Line 25 | Implementation:
25		#include <vector>
26		#include <iostream>
27		#include <sstream>
28	–	#include <algorithm>
28
29		// user include files
30		#include "FWCore/Framework/interface/Frameworkfwd.h"
#	Line 50 | Line 49 | Implementation:
49		#include <SimDataFormats/Vertex/interface/SimVertexContainer.h>
50		#include <SimDataFormats/Track/interface/SimTrack.h>
51		#include <SimDataFormats/Track/interface/SimTrackContainer.h>
52	+	// BeamSpot
53	+	#include "DataFormats/BeamSpot/interface/BeamSpot.h"
54	+
55
56		//root
57		#include <TROOT.h>
#	Line 60 | Line 62 | Implementation:
62		#include <TPad.h>
63
64		using namespace std;
65	+	typedef math::XYZTLorentzVectorF LorentzVector;
66	+	typedef math::XYZPoint Point;
67
68		PVStudy::PVStudy(const edm::ParameterSet& iConfig)
69		{
70		//=======================================================================
71		// Get configuration for TrackTupleMaker
72		//=======================================================================
73	<	simG4_           = iConfig.getParameter<edm::InputTag>( "simG4" );
74	<	trackCollectionTag_     = iConfig.getUntrackedParameter<edm::InputTag>("trackCollection");
75	<	splitTrackCollection1Tag_     = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection1");
76	<	splitTrackCollection2Tag_     = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
73	>	simG4_                     = iConfig.getParameter<edm::InputTag>( "simG4" );
74	>	trackCollectionTag_        = iConfig.getUntrackedParameter<edm::InputTag>("trackCollection");
75	>	splitTrackCollection1Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection1");
76	>	splitTrackCollection2Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
77		vertexCollectionTag_       = iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection");
78	<	splitVertexCollection1Tag_       = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
79	<	splitVertexCollection2Tag_       = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
80	<	pixelVertexCollectionTag_      = iConfig.getParameter<edm::InputTag>("pixelVertexCollectionTag");
81	<	verbose_         = iConfig.getUntrackedParameter<bool>("verbose",false);
82	<	realData_        = iConfig.getUntrackedParameter<bool>("realData",false);
83	<	analyze_         = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
84	<	saventuple_         = iConfig.getUntrackedParameter<bool>("saventuple",false);
85	<	outputfilename_  = iConfig.getUntrackedParameter<string>("OutputFileName");
86	<	ntrkdiffcut_         = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
87	<	zsigncut_         = iConfig.getUntrackedParameter<int>("zsigncut");
88	<	nTrkMin_ =  iConfig.getUntrackedParameter<int>("nTrkMin");
89	<	nTrkMax_ =  iConfig.getUntrackedParameter<int>("nTrkMax");
90	<
78	>	splitVertexCollection1Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
79	>	splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
80	>	verbose_                   = iConfig.getUntrackedParameter<bool>("verbose",false);
81	>	realData_                  = iConfig.getUntrackedParameter<bool>("realData",false);
82	>	analyze_                   = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
83	>	saventuple_                = iConfig.getUntrackedParameter<bool>("saventuple",false);
84	>	outputfilename_            = iConfig.getUntrackedParameter<string>("OutputFileName");
85	>	ntrkdiffcut_               = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
86	>	nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
87	>	nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
88	>	zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
89	>	bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
90	>	//   applyEvtClean_             = iConfig.getUntrackedParameter<bool>("applyEvtClean",false);
91	>	histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
92	>
93		//now do what ever initialization is needed
94		pvinfo.clear();
95
90	–
96		// Specify the data mode vector
97		if(realData_) datamode.push_back(0);
98		else {
#	Line 96 | Line 101 | PVStudy::PVStudy(const edm::ParameterSet
101		datamode.push_back(2);
102		datamode.push_back(3);
103		}
104	<
100	<	// Create ntuple files if needed
104	>	// Create ntuple files if needed
105		if(saventuple_) {
106		file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
107		ftree_ = new TTree("mytree","mytree");
#	Line 216 | Line 220 | PVStudy::PVStudy(const edm::ParameterSet
220		pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
221		}
222		}
219	–
220	–	edm::Service<TFileService> fs;
221	–	TFileDirectory subDir = fs->mkdir( "Summary" );
222	–	TFileDirectory subDir1 = fs->mkdir( "Others" );
223	–	//   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
224	–
225	–	setRootStyle();
223
224	<	//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_
224	>	setRootStyle();
225
226	<
227	<	for(int i=0;i<3;i++)
228	<	{
229	<	string suffix;
230	<	if(i == 0) suffix = "";
231	<	else {
232	<	stringstream ss;
233	<	ss << i;
234	<	suffix =ss.str()+"_spl";
235	<	}
236	<	h["nTrk"+suffix]   = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
237	<	h["trkPt"+suffix]  = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
238	<	h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#eta of rec tracks "+suffix),100,-3,3);
239	<	h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
240	<	h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-0.5,0.5);
241	<	h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
242	<
243	<	h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
244	<	h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
245	<	h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#eta of rec tracks in PV"+suffix),100,-3,3);
246	<	h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
247	<	h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);
248	<	h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
249	<	h["nrecPV"+suffix]   = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);
250	<	suffix.clear();
226	>	//========================================
227	>	// Booking histograms
228	>	//========================================
229	>
230	>	// Create a root file for histograms
231	>	theFile = new TFile(histoFileName_.c_str(), "RECREATE");
232	>	// make diretories
233	>	theFile->mkdir("Summary");
234	>	theFile->cd();
235	>	theFile->mkdir("Others");
236	>	theFile->cd();
237	>	if (analyze_) {
238	>	theFile->mkdir("Results");
239	>	theFile->cd();
240	>	}
241	>
242	>	// Book Histograms:
243	>	h_pvtrk   = new PVHistograms();
244	>	h_pixvtx  = new PVHistograms();
245	>	h_misc    = new PVHistograms();
246	>	h_summary = new PVHistograms();
247	>	h_others  = new PVHistograms();
248	>
249	>	for(int i=0;i<3;i++) {
250	>	if(i == 0) h_pvtrk->Init("pvTrk");
251	>	else {
252	>	stringstream ss;
253	>	ss << i;
254	>	h_pvtrk->Init("pvTrk", ss.str(),"spl");
255		}
256	<	h["nrecPVDiff"]     = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
257	<	h["nTrkPVDiff"]     = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
258	<	h["nTrkPVRelDiff"]  = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
259	<
260	<	// Histograms on comparing the multi-vertices
261	<	// Difference in reconstructed vtx position
262	<	h["min_xsep"]   = subDir.make<TH1D>("min_xsep", "min x diff of primary and secondary pvtx",300,0,0.1);
263	<	h["min_xsep_sign"]   = subDir.make<TH1D>("min_xsep_sign", "min x diff in signf of primary and secondary pvtx",300,0,5);
264	<	h["min_ysep"]   = subDir.make<TH1D>("min_ysep", "min y diff of primary and secondary pvtx",300,0,0.1);
265	<	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);
256	>	}
257	>	h_pixvtx->Init("pixVtx");
258	>	h_misc->Init("misc");
259	>
260	>	// Book MC only plots
261	>	if (!realData_) {
262	>	h_gen = new PVHistograms();
263	>	h_gen->Init("generator");
264	>	}
265	>	if (analyze_) h_ana = new PVHistograms();
266
267		//Book histograms sensitive to data/mc
268		for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
269		string suffix;
270	<	if(verbose_)  cout<<"datamode = "<< *it<<endl;
270	>	edm::LogInfo("Debug")<<"datamode = "<< *it<<endl;
271		switch(*it) {
272		case 0: suffix = "";
273		break;
#	Line 285 | Line 278 | PVStudy::PVStudy(const edm::ParameterSet
278		case 3: suffix = "_mct";
279		break;
280		}
281	<	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);
281	>	h_summary->Init("summary",suffix);
282
283		// Book residual, error and pull histograms for each nTrk bin
284		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
285		stringstream ss;
286		ss << ntrk;
287	<	string suffix2 = suffix+"_"+ss.str();
288	<	h["deltax"+suffix2]  = subDir1.make<TH1D>(TString("deltax"+suffix2), TString("x-residual of pvtx"+suffix),100,-0.02,0.02);
289	<	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	<
287	>	h_others->Init("others",suffix,ss.str());
288	>	}
289	>
290		// Book residual and pull histograms only when analyzeOntheFly is enabled
291	<	if(analyze_) {
320	<	h2["resx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resx"+suffix+"_nTrk"), TString("x-resolution vs number of tracks in pvtx"+suffix),nTrkMax_-nTrkMin_,nTrkMin_-0.5,nTrkMax_-0.5,400,0.,200);
321	<	h2["resx"+suffix+"_nTrk"]->SetMarkerStyle(21);
322	<	h2["resx"+suffix+"_nTrk"]->SetMarkerColor(4);
323	<	h2["resx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
324	<	h2["resx"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
325	<	h2["resy"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resy"+suffix+"_nTrk"), TString("y-resolution vs number of tracks in pvtx"+suffix),nTrkMax_-nTrkMin_,nTrkMin_-0.5,nTrkMax_-0.5,400,0.,200);
326	<	h2["resy"+suffix+"_nTrk"]->SetMarkerStyle(21);
327	<	h2["resy"+suffix+"_nTrk"]->SetMarkerColor(4);
328	<	h2["resy"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
329	<	h2["resy"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
330	<	h2["resz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resz"+suffix+"_nTrk"), TString("z-resolution vs number of tracks in pvtx"+suffix),nTrkMax_-nTrkMin_,nTrkMin_-0.5,nTrkMax_-0.5,400,0.,200);
331	<	h2["resz"+suffix+"_nTrk"]->SetMarkerStyle(21);
332	<	h2["resz"+suffix+"_nTrk"]->SetMarkerColor(4);
333	<	h2["resz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
334	<	h2["resz"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
335	<	h2["pullx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullx"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),nTrkMax_-nTrkMin_,nTrkMin_-0.5,nTrkMax_-0.5,100,0.,2.);
336	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerStyle(21);
337	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerColor(4);
338	<	h2["pullx"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
339	<	h2["pullx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
340	<	h2["pully"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pully"+suffix+"_nTrk"), TString("y-pull vs number of tracks"+suffix),nTrkMax_-nTrkMin_,nTrkMin_-0.5,nTrkMax_-0.5,100,0.,2.);
341	<	h2["pully"+suffix+"_nTrk"]->SetMarkerStyle(21);
342	<	h2["pully"+suffix+"_nTrk"]->SetMarkerColor(4);
343	<	h2["pully"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
344	<	h2["pully"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
345	<	h2["pullz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullz"+suffix+"_nTrk"), TString("z-pull vs number of tracks"+suffix),nTrkMax_-nTrkMin_,nTrkMin_-0.5,nTrkMax_-0.5,100,0.,2.);
346	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerStyle(21);
347	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerColor(4);
348	<	h2["pullz"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
349	<	h2["pullz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
350	<	} // End of   if(analyze_) {
291	>	if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
292		suffix.clear();
293		} // End of Book histograms sensitive to data/mc
294
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);
295
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	–	}
296		}
297
298		PVStudy::~PVStudy()
#	Line 379 | Line 300 | PVStudy::~PVStudy()
300
301		// do anything here that needs to be done at desctruction time
302		// (e.g. close files, deallocate resources etc.)
303	+	theFile->cd();
304	+	theFile->cd("Summary");
305	+	h_pvtrk->Save();
306	+	h_pixvtx->Save();
307	+	h_misc->Save();
308	+	h_summary->Save();
309	+	if (!realData_)
310	+	h_gen->Save();
311	+	if (analyze_) {
312	+	theFile->cd();
313	+	theFile->cd("Results");
314	+	h_ana->Save();
315	+	}
316	+	theFile->cd();
317	+	theFile->cd("Others");
318	+	h_others->Save();
319
320	+	theFile->Close();
321		}
322
323		void PVStudy::setRootStyle() {
#	Line 416 | Line 354 | void PVStudy::setRootStyle() {
354		//
355		std::vector<PVStudy::simPrimaryVertex> PVStudy::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
356		{
357	<	std::vector<PVStudy::simPrimaryVertex> simpv;
357	>	std::vector<PVStudy::simPrimaryVertex> simpv;
358		const HepMC::GenEvent* evt=evtMC->GetEvent();
359		if (evt) {
360	<	if(verbose_) {
361	<	cout << "process id " << evt->signal_process_id()<<endl;
362	<	cout << "signal process vertex " << ( evt->signal_process_vertex() ?
363	<	evt->signal_process_vertex()->barcode() : 0 )   <<endl;
426	<	cout << "number of vertices " << evt->vertices_size() << endl;
427	<	}
360	>	edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
361	>	edm::LogInfo("SimPVs") << "[getSimPVs] signal process vertex " << ( evt->signal_process_vertex() ?
362	>	evt->signal_process_vertex()->barcode() : 0 ) <<endl;
363	>	edm::LogInfo("SimPVs") << "[getSimPVs] number of vertices " << evt->vertices_size() << endl;
364
365		int idx=0; int npv=0;
366		for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
#	Line 435 | Line 371 | std::vector<PVStudy::simPrimaryVertex> P
371		// t component of PV:
372		for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
373		mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
374	<	//        std::cout << "Status = " << (*mother)->status() << std::endl;
374	>	//        edm::LogInfo("SimPVs") << "Status = " << (*mother)->status() << endl;
375		HepMC::GenVertex * mv=(*mother)->production_vertex();
376		if( ((*mother)->status() == 3) && (!mv)) {
377	<	//        std::cout << "npv= " << npv << std::endl;
377	>	//        edm::LogInfo("SimPVs") << "npv= " << npv << endl;
378		if (npv == 0) {
379	<	h["genPart_cT"]->Fill(pos.t()); // mm
380	<	h["genPart_T"]->Fill(pos.t()/299.792458); // ns
379	>	h_gen->Fill1d("genPart_cT", pos.t()); // mm
380	>	h_gen->Fill1d("genPart_T", pos.t()/299.792458); // ns
381		}
382		npv++;
383		}
#	Line 449 | Line 385 | std::vector<PVStudy::simPrimaryVertex> P
385		//       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
386
387		bool hasMotherVertex=false;
388	<	if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl;
388	>	if (verbose_) cout << "[getSimPVs] mothers of vertex[" << ++idx << "]: " << endl;
389		for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
390		mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
391		HepMC::GenVertex * mv=(*mother)->production_vertex();
392	<	//        std::cout << "Status = " << (*mother)->status() << std::endl;
392	>	// if (verbose_) cout << "Status = " << (*mother)->status() << endl;
393		if (mv) {
394		hasMotherVertex=true;
395		if(!verbose_) break; //if verbose_, print all particles of gen vertices
#	Line 480 | Line 416 | std::vector<PVStudy::simPrimaryVertex> P
416
417		if(!vp){
418		// this is a new vertex
419	<	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
419	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
420		simpv.push_back(sv);
421		vp=&simpv.back();
422		}else{
423	<	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
423	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
424		}
425		vp->genVertex.push_back((*vitr)->barcode());
426		// collect final state descendants
#	Line 508 | Line 444 | std::vector<PVStudy::simPrimaryVertex> P
444		}
445		}
446		}
447	<	}
448	<	}
447	>	}//loop MC vertices daughters
448	>	}//loop MC vertices
449		}
450		return simpv;
451		}
#	Line 542 | Line 478 | std::vector<PVStudy::simPrimaryVertex> P
478		ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
479		}
480		}
481	<	primary =  ( gv->position().t()==0);
481	>	//primary =  ( gv->position().t()==0 );
482	>	primary = true;
483	>	h_gen->Fill1d("genPart_cT", gv->position().t()); // mm
484	>	h_gen->Fill1d("genPart_T", gv->position().t()/299.792458); // ns
485	>
486		//resonance= ( gp->mother() && isResonance(gp->mother()));  // in CLHEP/HepMC days
487		// no more mother pointer in the improved HepMC GenParticle
488		resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
#	Line 591 | Line 531 | PVStudy::analyze(const edm::Event& iEven
531		using namespace edm;
532		using namespace reco;
533
534	<	if (verbose_)
595	<	cout<<"PVStudy::analyze():"<<endl;
534	>	edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
535		//=======================================================
536		// Initialize Root-tuple variables if needed
537		//=======================================================
538		//if(saventuple_)
539	<	SetVarToZero();
539	>	SetVarToZero();
540
541		//=======================================================
542		// Track accessors
#	Line 610 | Line 549 | PVStudy::analyze(const edm::Event& iEven
549		if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
550		trackColl = trackCollectionHandle.product();
551		} else {
552	<	cout << "trackCollection cannot be found -> using empty collection of same type." <<endl;
552	>	edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
553		}
615	–
554		//splitTrackCollection1
555		static const reco::TrackCollection s_empty_splitTrackColl1;
556		const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
#	Line 621 | Line 559 | PVStudy::analyze(const edm::Event& iEven
559		if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
560		splitTrackColl1 = splitTrackCollection1Handle.product();
561		} else {
562	<	cout << "splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
562	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
563		}
626	–
564		//splitTrackCollection2
565		static const reco::TrackCollection s_empty_splitTrackColl2;
566		const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
#	Line 632 | Line 569 | PVStudy::analyze(const edm::Event& iEven
569		if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
570		splitTrackColl2 = splitTrackCollection2Handle.product();
571		} else {
572	<	cout << "splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
572	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
573		}
574	<
638	<
639	<	//=======================================================
640	<	// Fill trackparameters of the input tracks to pvtx fitter
641	<	//=======================================================
642	<	if(verbose_)
643	<	cout<<"Start filling track parameters of the input tracks to pvtx fitter."<<endl;
644	<	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype)
645	<	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
646	<	fillTrackHisto(trackColl, 0);
647	<	fillTrackHisto(splitTrackColl1, 1);
648	<	fillTrackHisto(splitTrackColl2, 2);
649	<	if(verbose_)
650	<	cout<<"End filling track parameters of the input tracks to pvtx fitter."<<endl;
651	<
574	>
575		//=======================================================
576		// PVTX accessors
577		//=======================================================
#	Line 660 | Line 583 | PVStudy::analyze(const edm::Event& iEven
583		if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
584		vertexColl = vertexCollectionHandle.product();
585		} else {
586	<	cout << "vertexCollection cannot be found -> using empty collection of same type." <<endl;
586	>	edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
587		}
665	–
588		//splitVertexCollection1
589		static const reco::VertexCollection s_empty_splitVertexColl1;
590		const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
#	Line 671 | Line 593 | PVStudy::analyze(const edm::Event& iEven
593		if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
594		splitVertexColl1 = splitVertexCollection1Handle.product();
595		} else {
596	<	cout << "splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
597	<	}
676	<
596	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
597	>	}
598		//splitVertexCollection2
599		static const reco::VertexCollection s_empty_splitVertexColl2;
600		const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
#	Line 682 | Line 603 | PVStudy::analyze(const edm::Event& iEven
603		if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
604		splitVertexColl2 = splitVertexCollection2Handle.product();
605		} else {
606	<	cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
606	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
607	>	}
608	>
609	>
610	>	//=======================================================
611	>	// GET pixelVertices
612	>	//=======================================================
613	>	static const reco::VertexCollection s_empty_pixelVertexColl;
614	>	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
615	>	edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
616	>	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
617	>	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
618	>	pixelVertexColl = pixelVertexCollectionHandle.product();
619	>	} else {
620	>	edm::LogInfo("Debug") << "[PVStudy] pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
621		}
622	+
623	+	//=======================================================
624	+	// BeamSpot accessors
625	+	//=======================================================
626	+
627	+	edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
628	+	iEvent.getByLabel(bsSrc,recoBeamSpotHandle);
629	+	reco::BeamSpot bs = *recoBeamSpotHandle;
630	+	const Point beamSpot = recoBeamSpotHandle.isValid() ? Point(recoBeamSpotHandle->x0(), recoBeamSpotHandle->y0(), recoBeamSpotHandle->z0()) : Point(0, 0, 0);
631
632	<	if(verbose_) cout<<"End accessing the track and vertex collections"<<endl;
632	>	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex and pixelvertices collections"<<endl;
633
634		//=======================================================
635		// MC simvtx accessor
#	Line 704 | Line 648 | PVStudy::analyze(const edm::Event& iEven
648		try{
649		iSetup.getData(pdt);
650		}catch(const Exception&){
651	<	cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
651	>	edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
652		}
653	<
653	>
654		//=======================================================
655	<	// GET pixelVertices
655	>	// Apply event cleaning for firstcoll data and MC
656		//=======================================================
657	<	static const reco::VertexCollection s_empty_pixelVertexColl;
658	<	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
659	<	edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
660	<	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
661	<	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
662	<	pixelVertexColl = pixelVertexCollectionHandle.product();
663	<	} else {
664	<	cout << "pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
657	>
658	>	// =====Loop over the trackColl to tet the fraction of HighPurity tracks
659	>	int nTracks = 0;
660	>	int nHighPurityTracks=0;
661	>	double fHighPurity=0;
662	>
663	>	for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
664	>	TrackRef tkref(trackCollectionHandle,i);
665	>	if(tkref->quality(reco::TrackBase::highPurity)) ++nHighPurityTracks;
666		}
667	+	if(nTracks>0)
668	+	fHighPurity =  double(nHighPurityTracks)/double(nTracks);
669	+	if(verbose_)
670	+	cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;
671	+
672	+	if( (fHighPurity<0.2 && nTracks>10) || vertexColl->begin()->isFake()) {
673	+	//  glob_runno_ = iEvent.id().run();
674	+	//  glob_evtno_ = iEvent.id().event();
675	+	//  glob_ls_   = iEvent.luminosityBlock();
676	+	//  glob_bx_  = iEvent.bunchCrossing();
677	+	return;
678	+	}
679	+
680	+	//=======================================================
681	+	// Fill trackparameters of the input tracks to pvtx fitter
682	+	//=======================================================
683	+	edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
684	+	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs)
685	+	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
686	+	fillTrackHisto(trackColl, 0, beamSpot);
687	+	fillTrackHisto(splitTrackColl1, 1, beamSpot);
688	+	fillTrackHisto(splitTrackColl2, 2, beamSpot);
689	+	edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
690	+
691
692		//=======================================================
693		// Fill pixelVertices related histograms
694		//=======================================================
695		nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
696		if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
697	<	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
698	<	fillTrackHistoInPV(pixelVertexColl, 4, false, true);
699	<	h["dzErr_pxlpvtx"]->Fill( pixelVertexColl->begin()->zError());
697	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
698	>	fillTrackHistoInPV(pixelVertexColl, 4, false, true, beamSpot);
699	>	h_pixvtx->Fill1d("dzErr_pxlpvtx", pixelVertexColl->begin()->zError());
700		// Get the dZ error of the tracks in the leading pixelVertexColl
701		for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
702		t!= (pixelVertexColl->begin())->tracks_end(); t++) {
703
704		if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
705	<	//        h["trkPtPV"]->Fill(0.);
705	>	//        h_pvtrk->Fill1d("trkPtPV", 0.);
706		}
707		else
708	<	h["trkdzErr_pxlpvtx"]->Fill((**t).dzError());
708	>	h_pixvtx->Fill1d("trkdzErr_pxlpvtx", (**t).dzError());
709		}
710		// Fill the track->dz() in the PV difference from first pixelpvtx
711		for(reco::Vertex::trackRef_iterator t = (vertexColl->begin())->tracks_begin();
712		t!= (vertexColl->begin())->tracks_end(); t++) {
713		// illegal charge
714		if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
715	<	//        h["trkPtPV"]->Fill(0.);
715	>	//        h_pvtrk->Fill1d("trkPtPV", 0.);
716		}
717		else {
718		if(pixelVertexColl->size()>0) {
719	<	h["trkdz_pxlpvtxdz"]->Fill((**t).dz() -  pixelVertexColl->begin()->z());
720	<	h["trkdz_pxlpvtxdz_pxlpvtxdzerr"]->Fill(fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
721	<	h["trkdz_pxlpvtxdz_trkdzerr"]->Fill(fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
722	<	h["trkdzErr_pvtx"]->Fill((**t).dzError());
719	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz", (**t).dz() -  pixelVertexColl->begin()->z());
720	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_pxlpvtxdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
721	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_trkdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
722	>	h_pixvtx->Fill1d("trkdzErr_pvtx", (**t).dzError());
723		}
724		}
725		}
#	Line 760 | Line 729 | PVStudy::analyze(const edm::Event& iEven
729		// Fill number of reconstructed vertices
730		//=======================================================
731
732	<	if(verbose_)  {
733	<	cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
732	>	edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
733	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
734	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
735	>	edm::LogInfo("Debug")<<"[PVStudy] Start filling pvtx parameters."<<endl;
736	>	if (verbose_) {
737		printRecVtxs(vertexCollectionHandle);
766	–	cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
738		printRecVtxs(splitVertexCollection1Handle);
768	–	cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
739		printRecVtxs(splitVertexCollection2Handle);
770	–	cout<<"Start filling pvtx parameters."<<endl;
740		}
772	–
741		nrecPV_ = int (vertexColl->size());
742		nrecPV1_spl_ = int (splitVertexColl1->size());
743		nrecPV2_spl_ = int (splitVertexColl2->size());
744
745	<	h["nrecPV"]->Fill(nrecPV_);
746	<	h["nrecPV1_spl"]->Fill(nrecPV1_spl_);
747	<	h["nrecPV2_spl"]->Fill(nrecPV2_spl_);
748	<	h["nrecPVDiff"]->Fill(double(nrecPV1_spl_)-double(nrecPV2_spl_));
745	>	h_pvtrk->Fill1d("nrecPV", nrecPV_);
746	>	h_pvtrk->Fill1d("nrecPV1_spl", nrecPV1_spl_);
747	>	h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
748	>	h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
749
750	<	//=================================================================
750	>
751	>	//=================================================================
752		// Fill track parameter ntuple/hist for tracks used in recoVertices
753		//=================================================================
754
755	<	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
755	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
756		if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
757	<	fillTrackHistoInPV(vertexColl, 0, true, true);
757	>	fillTrackHistoInPV(vertexColl, 0, true, true, beamSpot);
758
759		if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
760	<	fillTrackHistoInPV(splitVertexColl1, 1, false, true);
792	<
793	<	if(splitVertexColl1->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
794	<	fillTrackHistoInPV(splitVertexColl2, 2, false, true);
760	>	fillTrackHistoInPV(splitVertexColl1, 1, false, true, beamSpot);
761
762	+	if(splitVertexColl2->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
763	+	fillTrackHistoInPV(splitVertexColl2, 2, false, true, beamSpot);
764
765		//=======================================================
766		// Compare offlinePrimaryVertices with pixelVertices
767		//=======================================================
768		if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
769		&& (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {
770	<	h["nrecPV_minus_nrecPxlPV"]->Fill( double (nrecPV_ - nrecPV_pxlpvtx_));
770	>	h_pixvtx->Fill1d("nrecPV_minus_nrecPxlPV", double (nrecPV_ - nrecPV_pxlpvtx_));
771		// difference in reconstructed position of the leading pvtx
772	<	//cout<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;
773	<	//cout<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
774	<	h["recxPV_minus_recxPxlPV"]->Fill (recx_[0] - recx_pxlpvtx_[0]);
775	<	h["recyPV_minus_recyPxlPV"]->Fill (recy_[0] - recy_pxlpvtx_[0]);
776	<	h["reczPV_minus_reczPxlPV"]->Fill (recz_[0] - recz_pxlpvtx_[0]);
772	>	//edm::LogInfo("Debug")<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;
773	>	//edm::LogInfo("Debug")<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
774	>	h_pixvtx->Fill1d("recxPV_minus_recxPxlPV", recx_[0] - recx_pxlpvtx_[0]);
775	>	h_pixvtx->Fill1d("recyPV_minus_recyPxlPV", recy_[0] - recy_pxlpvtx_[0]);
776	>	h_pixvtx->Fill1d("reczPV_minus_reczPxlPV", recz_[0] - recz_pxlpvtx_[0]);
777		}
778
779	<
780	<
813	<	//==========================================================
779	>
780	>	//==========================================================
781		// Fill secondary/primary min separations for multi-vertices
782		//==========================================================
783
#	Line 825 | Line 792 | PVStudy::analyze(const edm::Event& iEven
792		double min_ntrksep = 9999.0;
793		double min_sumpt2sep = 9999999.0;
794
795	<	if(verbose_) cout<<"leading pvtx z = "<< vertexColl->begin()->z() <<endl;
795	>	edm::LogInfo("Debug")<<"[PVStudy] leading pvtx z = "<< vertexColl->begin()->z() <<endl;
796
797		// Looping through the secondary vertices to find the mininum separation to the primary
798		for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
#	Line 860 | Line 827 | PVStudy::analyze(const edm::Event& iEven
827		min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
828		}
829		}
830	<	h["min_xsep"]->Fill(min_xsep);
831	<	h["min_ysep"]->Fill(min_ysep);
832	<	h["min_zsep"]->Fill(min_zsep);
833	<	h["min_xsep_sign"]->Fill(min_xsep_sign);
834	<	h["min_ysep_sign"]->Fill(min_ysep_sign);
835	<	h["min_zsep_sign"]->Fill(min_zsep_sign);
836	<	h["min_ntrksep"]->Fill(min_ntrksep);
837	<	h["min_sumpt2sep"]->Fill(min_sumpt2sep);
830	>	h_misc->Fill1d("min_xsep", min_xsep);
831	>	h_misc->Fill1d("min_ysep", min_ysep);
832	>	h_misc->Fill1d("min_zsep", min_zsep);
833	>	h_misc->Fill1d("min_xsep_sign", min_xsep_sign);
834	>	h_misc->Fill1d("min_ysep_sign", min_ysep_sign);
835	>	h_misc->Fill1d("min_zsep_sign", min_zsep_sign);
836	>	h_misc->Fill1d("min_ntrksep", min_ntrksep);
837	>	h_misc->Fill1d("min_sumpt2sep", min_sumpt2sep);
838
839		min_zsep_ = min_zsep;
840		min_ntrksep_ = min_ntrksep;
#	Line 876 | Line 843 | PVStudy::analyze(const edm::Event& iEven
843
844		} // End of  if(vertexColl->size()>1) {
845
846	<	if(verbose_)
847	<	cout<<"End filling pvtx parameters."<<endl;
846	>	edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
847	>
848
849	+
850		//=======================================================
851		//           PrimaryVertex Matching
852		// 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
#	Line 886 | Line 854 | PVStudy::analyze(const edm::Event& iEven
854		// Assume the first match is the primary vertex,
855		// since vertexColl are sorted in decreasing order of sum(pT)
856		//=======================================================
857	<
858	<	if(verbose_)   cout<<"PVStudy::analyze(): matching pvtxs "<<endl;
857	>
858	>	edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
859		reco::VertexCollection s_empty_matchedVertexColl1;
860		reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
861		matchedVertexColl1->reserve(splitVertexColl1->size());
#	Line 906 | Line 874 | PVStudy::analyze(const edm::Event& iEven
874		if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
875		reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
876		if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
877	<	if(verbose_) {
878	<	cout<<"Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
879	<	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
880	<	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
881	<	}
877	>	edm::LogInfo("Debug")<<"[PVStudy] Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
878	>	edm::LogInfo("Debug")<<"[PVStudy] recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
879	>	edm::LogInfo("Debug")<<"[PVStudy] recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
880	>
881	>	double twovtxsig = (recvtx1->z()-recvtx2->z())/max(recvtx1->zError(), recvtx2->zError());
882	>	h_misc->Fill1d("twovtxzsign", twovtxsig);
883		if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
884	<	if(verbose_) {
885	<	cout<<"The two splitted vertices match in Z. "<<endl;
917	<	}
884	>	edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
885	>
886		int nTrkPV1 = recvtx1->tracksSize();
887		int nTrkPV2 = recvtx2->tracksSize();
888		double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
889	<	h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
890	<	h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
889	>	h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
890	>	h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
891		if(fabs(ntrkreldiff)<ntrkdiffcut_) {
892	<	if(verbose_) {
893	<	cout<<"(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
926	<	}
892	>	edm::LogInfo("Debug")<<"[PVStudy] (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
893	>
894		matchedVertexColl1->push_back(*recvtx1);
895		matchedVertexColl2->push_back(*recvtx2);
896		stopmatching = true; // stop the matching when the first match is found!
897		break;
898		}
899		else
900	<	cout<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
900	>	edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
901		}
902		else {
903	<	cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
904	<	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
905	<	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
903	>	edm::LogInfo("Debug")<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
904	>	edm::LogInfo("Debug")<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
905	>	edm::LogInfo("Debug")<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
906		}
907		}
908		}
909	<	if(verbose_)     cout<<"PVStudy::analyze(): End matching pvtxs"<<endl;
909	>	edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
910
911		//=======================================================
912		//  Analyze matchedVertexColls
913		//=======================================================
914	<	if(verbose_) {
948	<	cout<<"Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
949	<	}
914	>	edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
915
916		// Compare the reconstructed vertex position and calculate resolution/pulls
917		if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
918		//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
919	<	fillTrackHistoInPV(matchedVertexColl1, 1, true, false);
920	<	fillTrackHistoInPV(matchedVertexColl2, 2, true, false);
919	>	fillTrackHistoInPV(matchedVertexColl1, 1, true, false, beamSpot);
920	>	fillTrackHistoInPV(matchedVertexColl2, 2, true, false, beamSpot);
921
922		reco::VertexCollection::const_iterator v1;
923		reco::VertexCollection::const_iterator v2;
924		nrecPV_twovtx_ = 0;
925		for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
926		v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
927	<	++v1, ++v2) {
927	>	++v1, ++v2) {
928	>	h_misc->Fill1d("ndofPVDiff", v1->ndof() - v2->ndof());
929	>	h_misc->Fill1d("ndofPVRelDiff", (v1->ndof()-v2->ndof())/(v1->ndof()+v2->ndof()));
930		fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
931		fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
932		fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
#	Line 968 | Line 935 | PVStudy::analyze(const edm::Event& iEven
935		ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
936		fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
937
938	<	h["deltax"]->Fill( fres_[0] );
939	<	h["deltay"]->Fill( fres_[1] );
940	<	h["deltaz"]->Fill( fres_[2] );
941	<	h["pullx"]->Fill( fres_[0]/ferror_[0]);
942	<	h["pully"]->Fill( fres_[1]/ferror_[1]);
943	<	h["pullz"]->Fill( fres_[2]/ferror_[2]);
944	<	h["errPVx"]->Fill( ferror_[0] );
945	<	h["errPVy"]->Fill( ferror_[1] );
946	<	h["errPVz"]->Fill( ferror_[2] );
938	>	h_summary->Fill1d("deltax", fres_[0] );
939	>	h_summary->Fill1d("deltay", fres_[1] );
940	>	h_summary->Fill1d("deltaz", fres_[2] );
941	>	h_summary->Fill1d("pullx", fres_[0]/ferror_[0]);
942	>	h_summary->Fill1d("pully", fres_[1]/ferror_[1]);
943	>	h_summary->Fill1d("pullz", fres_[2]/ferror_[2]);
944	>	h_summary->Fill1d("errPVx", ferror_[0] );
945	>	h_summary->Fill1d("errPVy", ferror_[1] );
946	>	h_summary->Fill1d("errPVz", ferror_[2] );
947		pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
948		error(ferror_[0], ferror_[1], ferror_[2]),
949		fntrk_));
#	Line 1001 | Line 968 | PVStudy::analyze(const edm::Event& iEven
968		pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];
969		nrecPV_twovtx_++;
970		} // End of analyzing res/pull
971	+
972	+	//=======================================================
973	+	// Fill simulated vertices
974	+	//=======================================================
975	+	if (!realData_) {
976	+	bool MC=false;
977	+	Handle<HepMCProduct> evtMC;
978	+	iEvent.getByLabel("generator",evtMC);
979	+	if (!evtMC.isValid()) {
980	+	MC=false;
981	+	edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
982	+	} else {
983	+	edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
984	+	MC=true;
985	+	}
986	+
987	+	if(MC){
988	+	// make a list of primary vertices:
989	+	std::vector<simPrimaryVertex> simpv;
990	+	simpv=getSimPVs(evtMC,"");
991	+	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
992	+	h_gen->Fill1d("nsimPV", simpv.size());
993	+	nsimPV_ = simpv.size();
994	+
995	+	int isimpv = 0;
996	+	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
997	+	vsim!=simpv.end(); vsim++, isimpv++){
998	+	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
999	+	simx_[isimpv] = vsim->x;
1000	+	simy_[isimpv] = vsim->y;
1001	+	simz_[isimpv] = vsim->y;
1002	+	simptsq_[isimpv] = vsim->ptsq;
1003	+
1004	+	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
1005	+	fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
1006	+	fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
1007	+	fillMCHisto(vsim, isimpv, vertexColl, 3);
1008	+	}
1009	+	} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
1010	+	} // End of   if (!realData_) {
1011	+
1012		} // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
1013		else
1014	<	cout<<"WARNING: Cannot find matching pvtxs"<<endl;
1014	>	edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
1015
1016	<	if(verbose_)     cout<<"End analyzing the matchedVertexColl"<<endl;
1016	>	edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
1017
1018	<	//=======================================================
1011	<	// Fill simulated vertices
1012	<	//=======================================================
1013	<	if (!realData_) {
1014	<	bool MC=false;
1015	<	Handle<HepMCProduct> evtMC;
1016	<	iEvent.getByLabel("generator",evtMC);
1017	<	if (!evtMC.isValid()) {
1018	<	MC=false;
1019	<	if(verbose_){
1020	<	cout << "no HepMCProduct found"<< endl;
1021	<	}
1022	<	} else {
1023	<	if(verbose_){
1024	<	cout << "generator HepMCProduct found"<< endl;
1025	<	}
1026	<	MC=true;
1027	<	}
1028	<
1029	<	if(MC){
1030	<	// make a list of primary vertices:
1031	<	std::vector<simPrimaryVertex> simpv;
1032	<	simpv=getSimPVs(evtMC,"");
1033	<	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
1034	<	h["nsimPV"]->Fill(simpv.size());
1035	<	nsimPV_ = simpv.size();
1036	<
1037	<	int isimpv = 0;
1038	<	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
1039	<	vsim!=simpv.end(); vsim++, isimpv++){
1040	<	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
1041	<	simx_[isimpv] = vsim->x;
1042	<	simy_[isimpv] = vsim->y;
1043	<	simz_[isimpv] = vsim->y;
1044	<	simptsq_[isimpv] = vsim->ptsq;
1045	<
1046	<	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
1047	<	fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
1048	<	fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
1049	<	fillMCHisto(vsim, isimpv, vertexColl, 3);
1050	<	}
1051	<	} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
1052	<	} // End of   if (!realData_) {
1018	>
1019
1020		// Finally fill the ftree_
1021		if(saventuple_) {
1022		ftree_->Fill();
1023		pvtxtree_->Fill();
1024		}
1025	+
1026		}
1027
1028	+
1029	+
1030		// ------------ method called once each job just before starting event loop  ------------
1031		void
1032		PVStudy::beginJob()
#	Line 1067 | Line 1036 | PVStudy::beginJob()
1036		// ------------ method called once each job just after ending the event loop  ------------
1037		void
1038		PVStudy::endJob() {
1039	<	if (verbose_) cout << "Analyzing PV info" << endl;
1039	>	edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl;
1040		// Looping through the datamodes available
1041		for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
1042		string suffix;
1043	<	if(verbose_)  cout<<"datamode = "<< *it<<endl;
1043	>	edm::LogInfo("Analysis")<<"[endJob] datamode = "<< *it<<endl;
1044		switch(*it) {
1045		case 0: suffix = "";
1046		break;
#	Line 1082 | Line 1051 | PVStudy::endJob() {
1051		case 3: suffix = "_mct";
1052		break;
1053		}
1054	+	suffix += "_nTrk";
1055		if(analyze_) {
1056		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1057	+	edm::LogInfo("Analysis") << "[endJob] ntrk = " << ntrk << endl;
1058		PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
1059	<	if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
1060	<	if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
1061	<	if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
1062	<	if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
1063	<	if ( ipv.pull_.y() > 0 ) h2["pully"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
1064	<	if ( ipv.pull_.z() > 0 ) h2["pullz"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.z());//, 1.);
1059	>	if ( ipv.res_.x() > 0 ) h_ana->Fill2d("resx"+suffix, ntrk,ipv.res_.x());
1060	>	if ( ipv.res_.y() > 0 ) h_ana->Fill2d("resy"+suffix, ntrk,ipv.res_.y());
1061	>	if ( ipv.res_.z() > 0 ) h_ana->Fill2d("resz"+suffix, ntrk,ipv.res_.z());
1062	>	if ( ipv.pull_.x() > 0 ) h_ana->Fill2d("pullx"+suffix, ntrk,ipv.pull_.x());
1063	>	if ( ipv.pull_.y() > 0 ) h_ana->Fill2d("pully"+suffix, ntrk,ipv.pull_.y());
1064	>	if ( ipv.pull_.z() > 0 ) h_ana->Fill2d("pullz"+suffix, ntrk,ipv.pull_.z());
1065		}
1066	<	}
1066	>	}
1067		suffix.clear();
1068		}
1069		if(saventuple_) {
#	Line 1101 | Line 1072 | PVStudy::endJob() {
1072		pvtxtree_->Write();
1073		file_->Close();
1074		}
1075	<	}
1105	<
1075	>
1076
1077	+	}
1078
1079		// Match a recovertex with a simvertex
1080		// ? Should the cut be parameter dependent?
#	Line 1111 | Line 1082 | PVStudy::endJob() {
1082		bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
1083		const reco::Vertex & vrec)
1084		{
1085	<	return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1085	>	if(vrec.isFake() || !vrec.isValid()) return false;
1086	>	else
1087	>	return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1088		}
1089
1090		// Match two reconstructed vertices
#	Line 1120 | Line 1093 | bool PVStudy::matchVertex( const reco::V
1093		int zsigncut)
1094		{
1095		double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
1096	<	if(verbose_)
1124	<	cout<<"The matching criteria in z is "<<cut<<endl;
1096	>	edm::LogInfo("Debug")<<"[matchVertex] The matching criteria in z is "<<cut<<endl;
1097		return (fabs(vrec1->z()-vrec2->z())<cut);
1098		}
1099
1100
1101		bool PVStudy::isResonance(const HepMC::GenParticle * p){
1102		double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
1103	<	if(verbose_) cout << "isResonance   " << p->pdg_id() << " " << ctau << endl;
1103	>	edm::LogInfo("Debug") << "[isResonance] isResonance   " << p->pdg_id() << " " << ctau << endl;
1104		return  ctau >0 && ctau <1e-6;
1105		}
1106
#	Line 1200 | Line 1172 | void PVStudy::printSimTrks(const Handle<
1172		void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1173		stringstream ss;
1174		ss << ntk;
1175	<	string suffix = ss.str();
1176	<	if(datamode == 0 )  suffix = "_" + suffix;
1177	<	if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1178	<	if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1179	<	if(datamode == 3 )  suffix = "_mct_" + suffix;
1175	>	string suffix;
1176	>	if(datamode == 0 )  suffix = "_" + ss.str();
1177	>	if(datamode == 1 )  suffix = "_spl1_mct_" + ss.str();
1178	>	if(datamode == 2 )  suffix = "_spl2_mct_" + ss.str();
1179	>	if(datamode == 3 )  suffix = "_mct_" + ss.str();
1180
1181		if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
1182		// Fill histograms of res and pull of ntk
1183	<	h["deltax"+suffix]->Fill(r.x());
1184	<	h["deltay"+suffix]->Fill(r.y());
1185	<	h["deltaz"+suffix]->Fill(r.z());
1186	<	h["pullx"+suffix]->Fill(r.x()/er.x());
1187	<	h["pully"+suffix]->Fill(r.y()/er.y());
1188	<	h["pullz"+suffix]->Fill(r.z()/er.z());
1189	<	h["errPVx"+suffix]->Fill(er.x());
1190	<	h["errPVy"+suffix]->Fill(er.y());
1191	<	h["errPVz"+suffix]->Fill(er.z());
1183	>	h_others->Fill1d("deltax"+suffix, r.x());
1184	>	h_others->Fill1d("deltay"+suffix, r.y());
1185	>	h_others->Fill1d("deltaz"+suffix, r.z());
1186	>	h_others->Fill1d("pullx"+suffix, r.x()/er.x());
1187	>	h_others->Fill1d("pully"+suffix, r.y()/er.y());
1188	>	h_others->Fill1d("pullz"+suffix, r.z()/er.z());
1189	>	h_others->Fill1d("errPVx"+suffix, er.x());
1190	>	h_others->Fill1d("errPVy"+suffix, er.y());
1191	>	h_others->Fill1d("errPVz"+suffix, er.z());
1192		}
1193
1194	+
1195		PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1196		map<int,double> data;
1197		data.clear();
#	Line 1236 | Line 1209 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1209		for ( int i = 0; i < 6; ++i) {
1210		switch (i) {
1211		case 0:
1212	<	fit(h["deltax"+suffix], fpar);
1212	>	if (!h_others->ReadHisto1D("deltax"+suffix)) break;
1213	>	fit(h_others->ReadHisto1D("deltax"+suffix), fpar);
1214		data[i] = fpar[0]*cm2um;
1215		break;
1216		case 1:
1217	<	fit(h["deltay"+suffix], fpar);
1217	>	if (!h_others->ReadHisto1D("deltay"+suffix)) break;
1218	>	fit(h_others->ReadHisto1D("deltay"+suffix), fpar);
1219		data[i] = fpar[0]*cm2um;
1220		break;
1221		case 2:
1222	<	fit(h["deltaz"+suffix], fpar);
1222	>	if (!h_others->ReadHisto1D("deltaz"+suffix)) break;
1223	>	fit(h_others->ReadHisto1D("deltaz"+suffix), fpar);
1224		data[i] = fpar[0]*cm2um;
1225		break;
1226		case 3:
1227	<	fit(h["pullx"+suffix], fpar);
1227	>	if (!h_others->ReadHisto1D("pullx"+suffix)) break;
1228	>	fit(h_others->ReadHisto1D("pullx"+suffix), fpar);
1229		data[i] = fpar[0];
1230		break;
1231		case 4:
1232	<	fit(h["pully"+suffix], fpar);
1232	>	if (!h_others->ReadHisto1D("pully"+suffix)) break;
1233	>	fit(h_others->ReadHisto1D("pully"+suffix), fpar);
1234		data[i] = fpar[0];
1235		break;
1236		case 5:
1237	<	fit(h["pullz"+suffix], fpar);
1237	>	if (!h_others->ReadHisto1D("pullz"+suffix)) break;
1238	>	fit(h_others->ReadHisto1D("pullz"+suffix), fpar);
1239		data[i] = fpar[0];
1240		break;
1241		}
1242	<	if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1242	>	if (data[i] > 0) edm::LogInfo("Analysis") << "[Analysis] data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1243		}
1244		}
1245		else
#	Line 1275 | Line 1254 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1254		ntk);
1255		}
1256
1257	<	void PVStudy::fit(TH1 *&hdist, double fitPars[]){
1257	>
1258	>	void PVStudy::fit(TH1 *hdist, double fitPars[]){
1259		TAxis *axis0 = hdist->GetXaxis();
1260		int nbin = axis0->GetLast();
1261		double nOF = hdist->GetBinContent(nbin+1);
1262		double nUF = hdist->GetBinContent(0);
1263		//   double fitRange = axis0->GetBinUpEdge(nbin);
1264	<	double fitRange = axis0->GetXmax();
1264	>	// double fitRange = axis0->GetXmax();
1265	>	double fitRange = 2*hdist->GetRMS();
1266		double sigMax[2] = {0.,0.};
1267
1268	<	if ( verbose_ ){
1269	<	cout << "Last bin = " << nbin;
1270	<	cout << "; hdist: Overflow Entries = " << nOF;
1271	<	cout << "; Underflow Entries = " << nUF;
1272	<	cout << "; hdist->GetEntries() = " << hdist->GetEntries();
1273	<	cout << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1293	<	}
1268	>	edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
1269	>	<< "; hdist: Overflow Entries = " << nOF
1270	>	<< "; Underflow Entries = " << nUF
1271	>	<< "; hdist->GetEntries() = " << hdist->GetEntries()
1272	>	<< "; fitting range = " << -fitRange << " to " << fitRange << endl;
1273	>
1274		if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit
1275		for (int bi = 0; bi < nbin; bi++) {
1276		if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) {
#	Line 1302 | Line 1282 | void PVStudy::fit(TH1 *&hdist, double fi
1282		if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1283		}
1284		}
1285	<	if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1285	>	//edm::LogInfo("Analysis") << "[Fit] Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1286
1287		TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1288		fgaus->SetParameter(1, 0.);
1289	<	fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
1289	>	fgaus->SetParLimits(1, -fitRange/10., fitRange/10.);
1290		fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
1291		fgaus->SetLineColor(4);
1292	<	hdist->Fit(fgaus,"Q");
1292	>	hdist->Fit(fgaus,"QLRM");
1293		gPad->Update();
1294		TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
1295		s->SetOptStat(1111111);
1296		s->SetOptFit(0111);
1297		gPad->Update();
1298	<	if (verbose_) cout << "got function" << endl;
1298	>	//edm::LogInfo("Analysis") << "[Fit] got function" << endl;
1299		fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1300		}
1301		else
1302		fitPars[0] = -999;
1303		}
1304
1305	<	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype) {
1305	>
1306	>	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs) {
1307		string suffix;
1308		suffix = ""; // for unsplit tracks
1309		if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1310		if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1311
1312	<	h["nTrk"+suffix]->Fill(trackColl->size());
1312	>	h_pvtrk->Fill1d("nTrk"+suffix, trackColl->size());
1313		for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1314		itTrack != trackColl->end();
1315		++itTrack) {
1316	<	h["trkPt"+suffix]->Fill(itTrack->pt());
1317	<	h["trkDxy"+suffix]->Fill(itTrack->dxy());
1318	<	h["trkDz"+suffix]->Fill(itTrack->dz());
1319	<	h["trkEta"+suffix]->Fill(itTrack->eta());
1320	<	h["trkPhi"+suffix]->Fill(itTrack->phi());
1316	>	h_pvtrk->Fill1d("trkPt"+suffix, itTrack->pt());
1317	>	h_pvtrk->Fill1d("trkDxy"+suffix, itTrack->dxy());
1318	>	h_pvtrk->Fill1d("trkDxyCorr"+suffix, itTrack->dxy(bs));
1319	>	h_pvtrk->Fill1d("trkDz"+suffix, itTrack->dz());
1320	>	h_pvtrk->Fill1d("trkEta"+suffix, itTrack->eta());
1321	>	h_pvtrk->Fill1d("trkPhi"+suffix, itTrack->phi());
1322		}
1323		}
1324
1325	<	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
1325	>	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
1326		string suffix;
1327		suffix = ""; // for unsplit tracks
1328		if(datatype == 1) suffix = "1_spl"; // for splittracks 1
#	Line 1348 | Line 1330 | void PVStudy::fillTrackHistoInPV(const r
1330		int ivtx = 0;
1331		for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1332		v!=vertexColl->end(); ++v, ++ivtx) {
1333	+	if(v->isFake()) continue;
1334		if(fillNtuple) {
1335		if(datatype == 0) {
1336		nTrkPV_[ivtx] = v->tracksSize();
#	Line 1360 | Line 1343 | void PVStudy::fillTrackHistoInPV(const r
1343		recx_err_[ivtx] = v->xError();
1344		recy_err_[ivtx] = v->yError();
1345		recz_err_[ivtx] = v->zError();
1363	–
1364	–
1346		}
1347		if(datatype == 1) {
1348		nTrkPV1_spl_[ivtx] = v->tracksSize();
#	Line 1404 | Line 1385 | void PVStudy::fillTrackHistoInPV(const r
1385		}
1386		// For histogram only fill the leading pvtx parameters
1387		if(fillHisto) {
1388	<	h["nTrkPV"+suffix]->Fill(vertexColl->begin()->tracksSize());
1388	>	h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1389	>	h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1390	>	int nHWTrkPV_ = 0;
1391	>
1392		try {
1393		for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1394		t!=vertexColl->begin()->tracks_end(); t++) {
1395		// illegal charge
1396		if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1397	<	//      h["trkPtPV"]->Fill(0.);
1397	>	//      h_pvtrk->Fill1d("trkPtPV", 0.);
1398		}
1399		else {
1400	<	h["trkPtPV"+suffix]->Fill((**t).pt());
1401	<	h["trkDxyPV"+suffix]->Fill((**t).dxy());
1402	<	h["trkDzPV"+suffix]->Fill((**t).dz());
1403	<	h["trkEtaPV"+suffix]->Fill((**t).eta());
1404	<	h["trkPhiPV"+suffix]->Fill((**t).phi());
1400	>	h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1401	>	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1402	>	h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
1403	>	h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1404	>	h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1405	>	h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1406	>
1407	>	if(vertexColl->begin()->trackWeight(*t) > 0.5)
1408	>	nHWTrkPV_++;
1409		}
1410		}
1411		}
1412		catch (...) {
1413		// exception thrown when trying to use linked track
1414	<	//          h["trkPtPV"]->Fill(0.);
1414	>	//          h_pvtrk->Fill1d("trkPtPV", 0.);
1415		}
1416	+	h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
1417		}
1418	+
1419
1420		}
1421
#	Line 1455 | Line 1445 | void PVStudy::fillMCHisto(std::vector<si
1445		int nrectrk = vrec->tracksSize();
1446		vsim->recVtx=NULL;
1447
1448	<	if(verbose_){
1449	<	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1450	<	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1451	<	}
1462	<	if ( matchVertex(*vsim,*vrec) && vrec->isValid() && !vrec->isFake() ) {
1448	>	edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1449	>	edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1450	>
1451	>	if ( matchVertex(*vsim,*vrec)  ) {
1452		vsim->recVtx=&(*vrec);
1453
1454		// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
#	Line 1467 | Line 1456 | void PVStudy::fillMCHisto(std::vector<si
1456		//|| (!vsim->recVtx) )
1457		//vsim->recVtx=&(*vrec);
1458		//}
1459	<	if(verbose_)
1471	<	cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1459	>	edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1460
1461		double fres_mct[3];
1462		double ferror_mct[3];
#	Line 1480 | Line 1468 | void PVStudy::fillMCHisto(std::vector<si
1468		ferror_mct[1] = vsim->recVtx->yError();
1469		ferror_mct[2] = vsim->recVtx->zError();
1470
1471	<	h["deltax"+suffix]->Fill( fres_mct[0] );
1472	<	h["deltay"+suffix]->Fill( fres_mct[1] );
1473	<	h["deltaz"+suffix]->Fill( fres_mct[2] );
1474	<	h["pullx"+suffix]->Fill( fres_mct[0]/ferror_mct[0] );
1475	<	h["pully"+suffix]->Fill( fres_mct[1]/ferror_mct[1] );
1476	<	h["pullz"+suffix]->Fill( fres_mct[2]/ferror_mct[2] );
1477	<	h["errPVx"+suffix]->Fill( ferror_mct[0] );
1478	<	h["errPVy"+suffix]->Fill( ferror_mct[1] );
1479	<	h["errPVz"+suffix]->Fill( ferror_mct[2] );
1471	>	h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1472	>	h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1473	>	h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1474	>	h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1475	>	h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1476	>	h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1477	>	h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1478	>	h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1479	>	h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1480		pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1481		error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1482		nrectrk));
#	Line 1533 | Line 1521 | void PVStudy::fillMCHisto(std::vector<si
1521		} // End of  if(saventuple_) {
1522		} //  if ( matchVertex(*vsim,*vrec) ) {
1523		else {  // no rec vertex found for this simvertex
1524	<	if(verbose_)
1537	<	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1524	>	edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1525		}
1526		}
1527		}
1528
1542	–
1529		void PVStudy::SetVarToZero() {
1530		//Greg's variables
1531		fntrk_ = 0;
1532		//pvtx position (x,y,z) residual and error
1533		for(int i = 0; i<3;i++)    {
1534	<	fres_[i] = 0;
1535	<	ferror_[i] = 0;
1534	>	fres_[i] = 0;
1535	>	ferror_[i] = 0;
1536		}
1537
1538		//Yanyan's variables
#	Line 1661 | Line 1647 | void PVStudy::SetVarToZero() {
1647		pully_spl2_mct_[i] = 0;
1648		pullz_spl2_mct_[i] = 0;
1649		}
1664	–
1650		}
1651
1652		double PVStudy::sumPtSquared(const reco::Vertex & v)  {
#	Line 1674 | Line 1659 | double PVStudy::sumPtSquared(const reco:
1659		return sum;
1660		}
1661
1662	+
1663	+

Diff Legend

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