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.7 by yygao, Thu Oct 8 14:57:46 2009 UTC vs.
Revision 1.20 by yygao, Thu Feb 18 10:59:48 2010 UTC

#	Line 49 | 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		//root
56		#include <TROOT.h>
#	Line 59 | Line 61 | Implementation:
61		#include <TPad.h>
62
63		using namespace std;
64	+	typedef math::XYZTLorentzVectorF LorentzVector;
65	+	typedef math::XYZPoint Point;
66
67		PVStudy::PVStudy(const edm::ParameterSet& iConfig)
68		{
69		//=======================================================================
70		// Get configuration for TrackTupleMaker
71		//=======================================================================
72	<	simG4_           = iConfig.getParameter<edm::InputTag>( "simG4" );
73	<	trackCollectionTag_     = iConfig.getUntrackedParameter<edm::InputTag>("trackCollection");
74	<	splitTrackCollection1Tag_     = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection1");
75	<	splitTrackCollection2Tag_     = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
72	>	simG4_                     = iConfig.getParameter<edm::InputTag>( "simG4" );
73	>	trackCollectionTag_        = iConfig.getUntrackedParameter<edm::InputTag>("trackCollection");
74	>	splitTrackCollection1Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection1");
75	>	splitTrackCollection2Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
76		vertexCollectionTag_       = iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection");
77	<	splitVertexCollection1Tag_       = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
78	<	splitVertexCollection2Tag_       = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
79	<	verbose_         = iConfig.getUntrackedParameter<bool>("verbose",false);
80	<	realData_        = iConfig.getUntrackedParameter<bool>("realData",false);
81	<	analyze_         = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
82	<	saventuple_         = iConfig.getUntrackedParameter<bool>("saventuple",false);
83	<	outputfilename_  = iConfig.getUntrackedParameter<string>("OutputFileName");
84	<	ntrkdiffcut_         = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
85	<	nTrkMin_ =  iConfig.getUntrackedParameter<int>("nTrkMin");
86	<	nTrkMax_ =  iConfig.getUntrackedParameter<int>("nTrkMax");
87	<
77	>	splitVertexCollection1Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
78	>	splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
79	>	verbose_                   = iConfig.getUntrackedParameter<bool>("verbose",false);
80	>	realData_                  = iConfig.getUntrackedParameter<bool>("realData",false);
81	>	analyze_                   = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
82	>	saventuple_                = iConfig.getUntrackedParameter<bool>("saventuple",false);
83	>	outputfilename_            = iConfig.getUntrackedParameter<string>("OutputFileName");
84	>	histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
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	>
91		//now do what ever initialization is needed
92		pvinfo.clear();
93
87	–
94		// Specify the data mode vector
95		if(realData_) datamode.push_back(0);
96		else {
#	Line 93 | Line 99 | PVStudy::PVStudy(const edm::ParameterSet
99		datamode.push_back(2);
100		datamode.push_back(3);
101		}
102	<
97	<	// Create ntuple files if needed
102	>	// Create ntuple files if needed
103		if(saventuple_) {
104		file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
105		ftree_ = new TTree("mytree","mytree");
#	Line 102 | Line 107 | PVStudy::PVStudy(const edm::ParameterSet
107		ftree_->Branch("residual",&fres_,"fres_[3]/D");
108		ftree_->Branch("error",&ferror_,"ferror_[3]/D");
109		ftree_->Branch("nTrkPV",&fntrk_,"fntrk_/I");
105	–	ftree_->Branch("nrecPV",&nrecPV_,"nrecPV_/I");
106	–	ftree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl_/I");
107	–	ftree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl_/I");
108	–
109	–	// If it is MC, get the resA, resB
110	–	if(!realData_) {
111	–	ftree_->Branch("residual_spl1_mct",&fres_spl1_mct_,"fres_spl1_mct_[3]/D");
112	–	ftree_->Branch("error_spl1_mct",&ferror_spl1_mct_,"ferror_spl1_mct_[3]/D");
113	–	ftree_->Branch("nTrkPV_spl1_mct",&fntrk_spl1_mct_,"fntrk_spl1_mct_/I");
114	–	ftree_->Branch("residual_spl2_mct",&fres_spl2_mct_,"fres_spl2_mct_[3]/D");
115	–	ftree_->Branch("error_spl2_mct",&ferror_spl2_mct_,"ferror_spl2_mct_[3]/D");
116	–	ftree_->Branch("nTrkPV_spl2_mct",&fntrk_spl2_mct_,"fntrk_spl2_mct_/I");
117	–	ftree_->Branch("residual_mct",&fres_mct_,"fres_mct_[3]/D");
118	–	ftree_->Branch("error_mct",&ferror_mct_,"ferror_mct_[3]/D");
119	–	ftree_->Branch("nTrkPV_mct",&fntrk_mct_,"fntrk_mct_/I");
120	–	}
121	–	}
110
111	<	edm::Service<TFileService> fs;
112	<	TFileDirectory subDir = fs->mkdir( "Summary" );
125	<	TFileDirectory subDir1 = fs->mkdir( "Others" );
126	<	//   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
111	>	// pvtxtree_ analyzing the pvtxs ootb
112	>	pvtxtree_ = new TTree("pvtxtree","pvtxtree");
113
114	<	setRootStyle();
114	>	// Event information for the data
115	>	pvtxtree_->Branch("glob_runno",&glob_runno_,"glob_runno/I");
116	>	pvtxtree_->Branch("glob_evtno",&glob_evtno_,"glob_evtno/I");
117	>	pvtxtree_->Branch("glob_ls",&glob_ls_,"glob_ls/I");
118	>	pvtxtree_->Branch("glob_bx",&glob_bx_,"glob_bx/I");
119	>
120	>
121	>	pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
122	>	pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
123	>	pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
124	>	// Event level information
125	>	pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
126	>	pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
127	>	pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
128	>
129	>	//Fill the variables in the twovtx pair (recvtx1, recvtx2)
130	>	// Information related to the analyzing the two-vertex method
131	>
132	>
133	>
134	>	pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
135	>	pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
136	>	pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
137	>	pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
138	>	pvtxtree_->Branch("deltaz_twovtx",&deltaz_twovtx_,"deltaz_twovtx[nrecPV_twovtx]/D");
139	>	pvtxtree_->Branch("errx_twovtx",&errx_twovtx_,"errx_twovtx[nrecPV_twovtx]/D");
140	>	pvtxtree_->Branch("erry_twovtx",&erry_twovtx_,"erry_twovtx[nrecPV_twovtx]/D");
141	>	pvtxtree_->Branch("errz_twovtx",&errz_twovtx_,"errz_twovtx[nrecPV_twovtx]/D");
142	>	pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
143	>	pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
144	>	pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
145	>
146	>	// Information for the splitVertexColl1
147	>	pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
148	>	pvtxtree_->Branch("ndofPV1_spl_twovtx",&ndofPV1_spl_twovtx_,"ndofPV1_spl_twovtx[nrecPV_twovtx]/D");
149	>	pvtxtree_->Branch("normchi2PV1_spl_twovtx",&normchi2PV1_spl_twovtx_,"normchi2PV1_spl_twovtx[nrecPV_twovtx]/D");
150	>	pvtxtree_->Branch("avgPtPV1_spl_twovtx",&avgPtPV1_spl_twovtx_,"avgPtPV1_spl_twovtx[nrecPV_twovtx]/D");
151	>	pvtxtree_->Branch("errx1_spl_twovtx",&errx1_spl_twovtx_,"errx1_spl_twovtx[nrecPV_twovtx]/D");
152	>	pvtxtree_->Branch("erry1_spl_twovtx",&erry1_spl_twovtx_,"erry1_spl_twovtx[nrecPV_twovtx]/D");
153	>	pvtxtree_->Branch("errz1_spl_twovtx",&errz1_spl_twovtx_,"errz1_spl_twovtx[nrecPV_twovtx]/D");
154
155	<	//Book histograms for track and pvtx parameters for each splitted pvtx
156	<	//i=0 : x (as in unsplit track collection)
157	<	//i=1 : x1_spl_
158	<	//i=2 : x2_spl_
159	<
160	<	for(int i=0;i<3;i++)
161	<	{
162	<	string suffix;
163	<	if(i == 0) suffix = "";
164	<	else {
165	<	stringstream ss;
166	<	ss << i;
167	<	suffix =ss.str()+"_spl";
168	<	}
169	<	h["nTrk"+suffix]   = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
170	<	h["trkPt"+suffix]  = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
171	<	h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#Eta of rec tracks "+suffix),100,-3,3);
172	<	h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
173	<	h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-5,5);
174	<	h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
175	<
176	<	h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
177	<	h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
178	<	h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#Eta of rec tracks in PV"+suffix),100,-3,3);
179	<	h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
180	<	h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);
181	<	h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
182	<	h["nrecPV"+suffix]   = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);
183	<	suffix.clear();
155	>
156	>	// Information for the splitVertexColl2
157	>	pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
158	>	pvtxtree_->Branch("ndofPV2_spl_twovtx",&ndofPV2_spl_twovtx_,"ndofPV2_spl_twovtx[nrecPV_twovtx]/D");
159	>	pvtxtree_->Branch("normchi2PV2_spl_twovtx",&normchi2PV2_spl_twovtx_,"normchi2PV2_spl_twovtx[nrecPV_twovtx]/D");
160	>	pvtxtree_->Branch("avgPtPV2_spl_twovtx",&avgPtPV2_spl_twovtx_,"avgPtPV2_spl_twovtx[nrecPV_twovtx]/D");
161	>	pvtxtree_->Branch("errx2_spl_twovtx",&errx2_spl_twovtx_,"errx2_spl_twovtx[nrecPV_twovtx]/D");
162	>	pvtxtree_->Branch("erry2_spl_twovtx",&erry2_spl_twovtx_,"erry2_spl_twovtx[nrecPV_twovtx]/D");
163	>	pvtxtree_->Branch("errz2_spl_twovtx",&errz2_spl_twovtx_,"errz2_spl_twovtx[nrecPV_twovtx]/D");
164	>
165	>
166	>	// MC variables
167	>	if(!realData_) {
168	>	pvtxtree_->Branch("nsimPV",&nsimPV_,"nsimPV/I");
169	>	pvtxtree_->Branch("nsimTrkPV",&nsimTrkPV_,"nsimTrkPV[nsimPV]/I");
170	>	pvtxtree_->Branch("simx",&simx_,"simx[nsimPV]/D");
171	>	pvtxtree_->Branch("simy",&simy_,"simy[nsimPV]/D");
172	>	pvtxtree_->Branch("simz",&simz_,"simz[nsimPV]/D");
173	>	pvtxtree_->Branch("simptsq",&simptsq_,"simptsq[nsimPV]/D");
174	>
175	>	// For pvtxs with all the tracks
176	>	pvtxtree_->Branch("nrecPV_mct",&nrecPV_mct_,"nrecPV_mct/I");
177	>	pvtxtree_->Branch("deltax_mct",&deltax_mct_,"deltax_mct[nrecPV_mct]/D");
178	>	pvtxtree_->Branch("deltay_mct",&deltay_mct_,"deltay_mct[nrecPV_mct]/D");
179	>	pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
180	>	pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
181	>	pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
182	>	pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
183	>	pvtxtree_->Branch("errx_mct",&errx_mct_,"errx_mct[nrecPV_mct]/D");
184	>	pvtxtree_->Branch("erry_mct",&erry_mct_,"erry_mct[nrecPV_mct]/D");
185	>	pvtxtree_->Branch("errz_mct",&errz_mct_,"errz_mct[nrecPV_mct]/D");
186	>	pvtxtree_->Branch("nTrkPV_mct",&nTrkPV_mct_,"nTrkPV_mct[nrecPV_mct]/I");
187	>	pvtxtree_->Branch("ndofPV_mct",&ndofPV_mct_,"ndofPV_mct[nrecPV_mct]/D");
188	>	pvtxtree_->Branch("normchi2PV_mct",&normchi2PV_mct_,"normchi2PV_mct[nrecPV_mct]/D");
189	>	pvtxtree_->Branch("avgPtPV_mct",&avgPtPV_mct_,"avgPtPV_mct[nrecPV_mct]/D");
190	>
191	>	// For pvtxs with splittracks1
192	>	pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
193	>	pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
194	>	pvtxtree_->Branch("deltay_spl1_mct",&deltay_spl1_mct_,"deltay_spl1_mct[nrecPV_spl1_mct]/D");
195	>	pvtxtree_->Branch("deltaz_spl1_mct",&deltaz_spl1_mct_,"deltaz_spl1_mct[nrecPV_spl1_mct]/D");
196	>	pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
197	>	pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
198	>	pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
199	>	pvtxtree_->Branch("errx_spl1_mct",&errx_spl1_mct_,"errx_spl1_mct[nrecPV_spl1_mct]/D");
200	>	pvtxtree_->Branch("erry_spl1_mct",&erry_spl1_mct_,"erry_spl1_mct[nrecPV_spl1_mct]/D");
201	>	pvtxtree_->Branch("errz_spl1_mct",&errz_spl1_mct_,"errz_spl1_mct[nrecPV_spl1_mct]/D");
202	>	pvtxtree_->Branch("nTrkPV_spl1_mct",&nTrkPV_spl1_mct_,"nTrkPV_spl1_mct[nrecPV_spl1_mct]/I");
203	>	pvtxtree_->Branch("ndofPV_spl1_mct",&ndofPV_spl1_mct_,"ndofPV_spl1_mct[nrecPV_spl1_mct]/D");
204	>	pvtxtree_->Branch("normchi2PV_spl1_mct",&normchi2PV_spl1_mct_,"normchi2PV_spl1_mct[nrecPV_spl1_mct]/D");
205	>	pvtxtree_->Branch("avgPtPV_spl1_mct",&avgPtPV_spl1_mct_,"avgPtPV_spl1_mct[nrecPV_spl1_mct]/D");
206	>
207	>	// For pvtxs with splittracks1
208	>	pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
209	>	pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
210	>	pvtxtree_->Branch("deltay_spl2_mct",&deltay_spl2_mct_,"deltay_spl2_mct[nrecPV_spl2_mct]/D");
211	>	pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
212	>	pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
213	>	pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
214	>	pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
215	>	pvtxtree_->Branch("errx_spl2_mct",&errx_spl2_mct_,"errx_spl2_mct[nrecPV_spl2_mct]/D");
216	>	pvtxtree_->Branch("erry_spl2_mct",&erry_spl2_mct_,"erry_spl2_mct[nrecPV_spl2_mct]/D");
217	>	pvtxtree_->Branch("errz_spl2_mct",&errz_spl2_mct_,"errz_spl2_mct[nrecPV_spl2_mct]/D");
218	>	pvtxtree_->Branch("nTrkPV_spl2_mct",&nTrkPV_spl2_mct_,"nTrkPV_spl2_mct[nrecPV_spl2_mct]/I");
219	>	pvtxtree_->Branch("ndofPV_spl2_mct",&ndofPV_spl2_mct_,"ndofPV_spl2_mct[nrecPV_spl2_mct]/D");
220	>	pvtxtree_->Branch("normchi2PV_spl2_mct",&normchi2PV_spl2_mct_,"normchi2PV_spl2_mct[nrecPV_spl2_mct]/D");
221	>	pvtxtree_->Branch("avgPtPV_spl2_mct",&avgPtPV_spl2_mct_,"avgPtPV_spl2_mct[nrecPV_spl2_mct]/D");
222	>
223		}
224	<	h["nrecPVDiff"]     = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
225	<	h["nTrkPVDiff"]     = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
226	<	h["nTrkPVRelDiff"]  = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
224	>	}
225	>
226	>	setRootStyle();
227	>
228	>	//========================================
229	>	// Booking histograms
230	>	//========================================
231
232	+	// Create a root file for histograms
233	+	theFile = new TFile(histoFileName_.c_str(), "RECREATE");
234	+	// make diretories
235	+	theFile->mkdir("Summary");
236	+	theFile->cd();
237	+	theFile->mkdir("Others");
238	+	theFile->cd();
239	+	if (analyze_) {
240	+	theFile->mkdir("Results");
241	+	theFile->cd();
242	+	}
243	+
244	+	// Book Histograms:
245	+	h_pvtrk   = new PVHistograms();
246	+	h_misc    = new PVHistograms();
247	+	h_summary = new PVHistograms();
248	+	h_others  = new PVHistograms();
249	+
250	+	for(int i=0;i<3;i++) {
251	+	if(i == 0) h_pvtrk->Init("pvTrk");
252	+	else {
253	+	stringstream ss;
254	+	ss << i;
255	+	h_pvtrk->Init("pvTrk", ss.str(),"spl");
256	+	}
257	+	}
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 175 | 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);
179	<	h["deltay"+suffix] = subDir.make<TH1D>(TString("deltay"+suffix), TString("y-residual pvtx"+suffix),800,-0.04,0.04);
180	<	h["deltaz"+suffix] = subDir.make<TH1D>(TString("deltaz"+suffix), TString("z-residual pvtx"+suffix),800,-0.04,0.04);
181	<	h["pullx"+suffix]  = subDir.make<TH1D>(TString("pullx"+suffix), TString("x-pull pvtx"+suffix),800,-10,10);
182	<	h["pully"+suffix]  = subDir.make<TH1D>(TString("pully"+suffix), TString("y-pull pvtx"+suffix),800,-10,10);
183	<	h["pullz"+suffix]  = subDir.make<TH1D>(TString("pullz"+suffix), TString("z-pull pvtx"+suffix),800,-10,10);
184	<	h["errPVx"+suffix] = subDir.make<TH1D>(TString("errPVx"+suffix), TString("X"+suffix+" vertex error"),100,0.,0.02);
185	<	h["errPVy"+suffix] = subDir.make<TH1D>(TString("errPVy"+suffix), TString("Y"+suffix+" vertex error"),100,0.,0.02);
186	<	h["errPVz"+suffix] = subDir.make<TH1D>(TString("errPVz"+suffix), TString("Z"+suffix+" vertex error"),100,0.,0.02);
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");
195	<	h["deltay"+suffix2]  = subDir1.make<TH1D>(TString("deltay"+suffix2), TString("y-residual of pvtx"+suffix),100,-0.02,0.02);
196	<	h["deltay"+suffix2]->GetXaxis()->SetTitle("cm");
197	<	h["deltaz"+suffix2]  = subDir1.make<TH1D>(TString("deltaz"+suffix2), TString("z-residual of pvtx"+suffix),100,-0.02,0.02);
198	<	h["deltaz"+suffix2]->GetXaxis()->SetTitle("cm");
199	<	h["pullx"+suffix2] = subDir1.make<TH1D>(TString("pullx"+suffix2), TString("x-pull of pvtx"+suffix),100,-10.,10.);
200	<	h["pully"+suffix2] = subDir1.make<TH1D>(TString("pully"+suffix2), TString("y-pull of pvtx"+suffix),100,-10.,10.);
201	<	h["pullz"+suffix2] = subDir1.make<TH1D>(TString("pullz"+suffix2), TString("z-pull of pvtx"+suffix),100,-10.,10.);
202	<	h["errPVx"+suffix2] = subDir1.make<TH1D>(TString("errPVx"+suffix2), TString("X"+suffix+" vertex error"),100,0.,0.02);
203	<	h["errPVy"+suffix2] = subDir1.make<TH1D>(TString("errPVy"+suffix2), TString("Y"+suffix+" vertex error"),100,0.,0.02);
204	<	h["errPVz"+suffix2] = subDir1.make<TH1D>(TString("errPVz"+suffix2), TString("Z"+suffix+" vertex error"),100,0.,0.02);
205	<	suffix2.clear();
206	<	} // End of   for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
207	<
287	>	h_others->Init("others",suffix,ss.str());
288	>	}
289	>
290		// Book residual and pull histograms only when analyzeOntheFly is enabled
291	<	if(analyze_) {
210	<	h2["resx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resx"+suffix+"_nTrk"), TString("x-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
211	<	h2["resx"+suffix+"_nTrk"]->SetMarkerStyle(21);
212	<	h2["resx"+suffix+"_nTrk"]->SetMarkerColor(4);
213	<	h2["resx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
214	<	h2["resx"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
215	<	h2["resy"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resy"+suffix+"_nTrk"), TString("y-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
216	<	h2["resy"+suffix+"_nTrk"]->SetMarkerStyle(21);
217	<	h2["resy"+suffix+"_nTrk"]->SetMarkerColor(4);
218	<	h2["resy"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
219	<	h2["resy"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
220	<	h2["resz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resz"+suffix+"_nTrk"), TString("z-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
221	<	h2["resz"+suffix+"_nTrk"]->SetMarkerStyle(21);
222	<	h2["resz"+suffix+"_nTrk"]->SetMarkerColor(4);
223	<	h2["resz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
224	<	h2["resz"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
225	<	h2["pullx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullx"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
226	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerStyle(21);
227	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerColor(4);
228	<	h2["pullx"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
229	<	h2["pullx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
230	<	h2["pully"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pully"+suffix+"_nTrk"), TString("y-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
231	<	h2["pully"+suffix+"_nTrk"]->SetMarkerStyle(21);
232	<	h2["pully"+suffix+"_nTrk"]->SetMarkerColor(4);
233	<	h2["pully"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
234	<	h2["pully"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
235	<	h2["pullz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullz"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
236	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerStyle(21);
237	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerColor(4);
238	<	h2["pullz"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
239	<	h2["pullz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
240	<	} // End of   if(analyze_) {
291	>	if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
292		suffix.clear();
293	<	} // End of Book histograms sensitive to data/mc
294	<	// Book MC only plots
295	<	if (!realData_) {
245	<	h["genPart_T"]      = subDir.make<TH1D>("genPart_T","t component of gen particles",300,-0.5,0.5);
246	<	h["genPart_T"]->GetXaxis()->SetTitle("t (nanosecond)");
247	<	h["genPart_cT"]     = subDir.make<TH1D>("genPart_cT","ct component of gen particles",300,-150.,150.);
248	<	h["genPart_cT"]->GetXaxis()->SetTitle("ct (mm)");
249	<	h["nsimPV"]         = subDir.make<TH1D>("nsimPV","Num of sim PV",51,-0.5,50.5);
250	<	}
293	>	} // End of Book histograms sensitive to data/mc
294	>
295	>
296		}
297
298		PVStudy::~PVStudy()
#	Line 255 | 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_misc->Save();
307	+	h_summary->Save();
308	+	if (!realData_)
309	+	h_gen->Save();
310	+	if (analyze_) {
311	+	theFile->cd();
312	+	theFile->cd("Results");
313	+	h_ana->Save();
314	+	}
315	+	theFile->cd();
316	+	theFile->cd("Others");
317	+	h_others->Save();
318
319	+	theFile->Close();
320		}
321
322		void PVStudy::setRootStyle() {
#	Line 279 | Line 340 | void PVStudy::setRootStyle() {
340		gStyle->SetTitleSize(0.055, "");    // size for pad title; default is 0.02
341		gStyle->SetLabelSize(0.03, "XYZ");  // size for axis labels; default is 0.04
342		gStyle->SetStatFontSize(0.08);      // size for stat. box
343	<	gStyle->SetTitleFont(32, "XYZ");    // times-bold-italic font (p. 153) for axes
344	<	gStyle->SetTitleFont(32, "");       // same for pad title
345	<	gStyle->SetLabelFont(32, "XYZ");    // same for axis labels
346	<	gStyle->SetStatFont(32);            // same for stat. box
343	>	gStyle->SetTitleFont(42, "XYZ");    // times-bold-italic font (p. 153) for axes
344	>	gStyle->SetTitleFont(42, "");       // same for pad title
345	>	gStyle->SetLabelFont(42, "XYZ");    // same for axis labels
346	>	gStyle->SetStatFont(42);            // same for stat. box
347		gStyle->SetLabelOffset(0.006, "Y"); // default is 0.005
348		//
349		return;
#	Line 292 | Line 353 | void PVStudy::setRootStyle() {
353		//
354		std::vector<PVStudy::simPrimaryVertex> PVStudy::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
355		{
356	<	std::vector<PVStudy::simPrimaryVertex> simpv;
356	>	std::vector<PVStudy::simPrimaryVertex> simpv;
357		const HepMC::GenEvent* evt=evtMC->GetEvent();
358		if (evt) {
359	<	if(verbose_) {
360	<	cout << "process id " << evt->signal_process_id()<<endl;
361	<	cout << "signal process vertex " << ( evt->signal_process_vertex() ?
362	<	evt->signal_process_vertex()->barcode() : 0 )   <<endl;
302	<	cout << "number of vertices " << evt->vertices_size() << endl;
303	<	}
359	>	edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
360	>	edm::LogInfo("SimPVs") << "[getSimPVs] signal process vertex " << ( evt->signal_process_vertex() ?
361	>	evt->signal_process_vertex()->barcode() : 0 ) <<endl;
362	>	edm::LogInfo("SimPVs") << "[getSimPVs] number of vertices " << evt->vertices_size() << endl;
363
364		int idx=0; int npv=0;
365		for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
#	Line 311 | Line 370 | std::vector<PVStudy::simPrimaryVertex> P
370		// t component of PV:
371		for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
372		mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
373	<	//        std::cout << "Status = " << (*mother)->status() << std::endl;
373	>	//        edm::LogInfo("SimPVs") << "Status = " << (*mother)->status() << endl;
374		HepMC::GenVertex * mv=(*mother)->production_vertex();
375		if( ((*mother)->status() == 3) && (!mv)) {
376	<	//        std::cout << "npv= " << npv << std::endl;
376	>	//        edm::LogInfo("SimPVs") << "npv= " << npv << endl;
377		if (npv == 0) {
378	<	h["genPart_cT"]->Fill(pos.t()); // mm
379	<	h["genPart_T"]->Fill(pos.t()/299.792458); // ns
378	>	h_gen->Fill1d("genPart_cT", pos.t()); // mm
379	>	h_gen->Fill1d("genPart_T", pos.t()/299.792458); // ns
380		}
381		npv++;
382		}
#	Line 325 | Line 384 | std::vector<PVStudy::simPrimaryVertex> P
384		//       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
385
386		bool hasMotherVertex=false;
387	<	if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl;
387	>	if (verbose_) cout << "[getSimPVs] mothers of vertex[" << ++idx << "]: " << endl;
388		for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
389		mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
390		HepMC::GenVertex * mv=(*mother)->production_vertex();
391	<	//        std::cout << "Status = " << (*mother)->status() << std::endl;
391	>	// if (verbose_) cout << "Status = " << (*mother)->status() << endl;
392		if (mv) {
393		hasMotherVertex=true;
394		if(!verbose_) break; //if verbose_, print all particles of gen vertices
#	Line 356 | Line 415 | std::vector<PVStudy::simPrimaryVertex> P
415
416		if(!vp){
417		// this is a new vertex
418	<	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
418	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
419		simpv.push_back(sv);
420		vp=&simpv.back();
421		}else{
422	<	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
422	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
423		}
424		vp->genVertex.push_back((*vitr)->barcode());
425		// collect final state descendants
#	Line 384 | Line 443 | std::vector<PVStudy::simPrimaryVertex> P
443		}
444		}
445		}
446	<	}
447	<	}
446	>	}//loop MC vertices daughters
447	>	}//loop MC vertices
448		}
449		return simpv;
450		}
#	Line 418 | Line 477 | std::vector<PVStudy::simPrimaryVertex> P
477		ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
478		}
479		}
480	<	primary =  ( gv->position().t()==0);
480	>	//primary =  ( gv->position().t()==0 );
481	>	primary = true;
482	>	h_gen->Fill1d("genPart_cT", gv->position().t()); // mm
483	>	h_gen->Fill1d("genPart_T", gv->position().t()/299.792458); // ns
484	>
485		//resonance= ( gp->mother() && isResonance(gp->mother()));  // in CLHEP/HepMC days
486		// no more mother pointer in the improved HepMC GenParticle
487		resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
#	Line 467 | Line 530 | PVStudy::analyze(const edm::Event& iEven
530		using namespace edm;
531		using namespace reco;
532
533	<	if (verbose_)
534	<	cout<<"PVStudy::analyze():"<<endl;
535	<	//=======================================================
533	>	//========================================================================
534	>	// Step 0: Prepare root variables and get information from the Event
535	>	//========================================================================
536	>
537	>	edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
538		// Initialize Root-tuple variables if needed
539	<	//=======================================================
540	<	if(saventuple_)
541	<	SetVarToZero();
477	<
478	<	//=======================================================
479	<	// Track accessors
480	<	//=======================================================
481	<	//trackCollection
539	>	SetVarToZero();
540	>
541	>	// ====== TrackCollection
542		static const reco::TrackCollection s_empty_trackColl;
543		const reco::TrackCollection *trackColl = &s_empty_trackColl;
544		edm::Handle<reco::TrackCollection> trackCollectionHandle;
#	Line 486 | Line 546 | PVStudy::analyze(const edm::Event& iEven
546		if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
547		trackColl = trackCollectionHandle.product();
548		} else {
549	<	cout << "trackCollection cannot be found -> using empty collection of same type." <<endl;
549	>	edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
550		}
551	<	//splitTrackCollection1
551	>	// ====== splitTrackCollection1
552		static const reco::TrackCollection s_empty_splitTrackColl1;
553		const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
554		edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
#	Line 496 | Line 556 | PVStudy::analyze(const edm::Event& iEven
556		if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
557		splitTrackColl1 = splitTrackCollection1Handle.product();
558		} else {
559	<	cout << "splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
559	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
560		}
561	<	//splitTrackCollection2
561	>	// ====== splitTrackCollection2
562		static const reco::TrackCollection s_empty_splitTrackColl2;
563		const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
564		edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
#	Line 506 | Line 566 | PVStudy::analyze(const edm::Event& iEven
566		if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
567		splitTrackColl2 = splitTrackCollection2Handle.product();
568		} else {
569	<	cout << "splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
569	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
570		}
571
572	<	//=======================================================
513	<	// PVTX accessors
514	<	//=======================================================
515	<	//vertexCollection
572	>	// ======= PrimaryVertexCollection
573		static const reco::VertexCollection s_empty_vertexColl;
574		const reco::VertexCollection *vertexColl = &s_empty_vertexColl;
575		edm::Handle<reco::VertexCollection> vertexCollectionHandle;
#	Line 520 | Line 577 | PVStudy::analyze(const edm::Event& iEven
577		if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
578		vertexColl = vertexCollectionHandle.product();
579		} else {
580	<	cout << "vertexCollection cannot be found -> using empty collection of same type." <<endl;
580	>	edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
581		}
582	<	//splitVertexCollection1
582	>	// ====== splitVertexCollection1
583		static const reco::VertexCollection s_empty_splitVertexColl1;
584		const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
585		edm::Handle<reco::VertexCollection> splitVertexCollection1Handle;
#	Line 530 | Line 587 | PVStudy::analyze(const edm::Event& iEven
587		if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
588		splitVertexColl1 = splitVertexCollection1Handle.product();
589		} else {
590	<	cout << "splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
590	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
591		}
592	<	//splitVertexCollection2
592	>	// ====== splitVertexCollection2
593		static const reco::VertexCollection s_empty_splitVertexColl2;
594		const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
595		edm::Handle<reco::VertexCollection> splitVertexCollection2Handle;
#	Line 540 | Line 597 | PVStudy::analyze(const edm::Event& iEven
597		if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
598		splitVertexColl2 = splitVertexCollection2Handle.product();
599		} else {
600	<	cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
600	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
601		}
602
546	–	if(verbose_) cout<<"Done accessing the track and vertex collections"<<endl;
603
604	<	//=======================================================
605	<	// MC simvtx accessor
606	<	//=======================================================
604	>	// ======== BeamSpot accessors
605	>	edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
606	>	iEvent.getByLabel(bsSrc,recoBeamSpotHandle);
607	>	reco::BeamSpot bs = *recoBeamSpotHandle;
608	>	const Point beamSpot = recoBeamSpotHandle.isValid() ? Point(recoBeamSpotHandle->x0(), recoBeamSpotHandle->y0(), recoBeamSpotHandle->z0()) : Point(0, 0, 0);
609	>
610	>	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex collections"<<endl;
611	>
612	>	// ========== MC simvtx accessor
613		if (!realData_) {
614		edm::Handle<SimVertexContainer> simVtxs;
615		iEvent.getByLabel( simG4_, simVtxs);
#	Line 556 | Line 618 | PVStudy::analyze(const edm::Event& iEven
618		iEvent.getByLabel( simG4_, simTrks);
619		}
620
621	<	//=======================================================
560	<	// GET PDT
561	<	//=======================================================
621	>	// ========== GET PDT
622		try{
623		iSetup.getData(pdt);
624		}catch(const Exception&){
625	<	cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
625	>	edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
626	>	}
627	>
628	>	// ======= Get MC information if needed
629	>	bool MC=false;
630	>	Handle<HepMCProduct> evtMC;
631	>	if (!realData_) {
632	>	iEvent.getByLabel("generator",evtMC);
633	>	if (!evtMC.isValid()) {
634	>	MC=false;
635	>	edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
636	>	} else {
637	>	edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
638	>	MC=true;
639	>	}
640		}
641
642	<	//=======================================================
643	<	// Fill trackparameters of the input tracks to pvtx fitter
644	<	//=======================================================
645	<
646	<	// GeneralTracks
647	<	h["nTrk"]->Fill(trackColl->size());
648	<	for(TrackCollection::const_iterator itTrack = trackColl->begin();
649	<	itTrack != trackColl->end();
650	<	++itTrack) {
651	<	h["trkPt"]->Fill(itTrack->pt());
652	<	h["trkDxy"]->Fill(itTrack->dxy());
653	<	h["trkDz"]->Fill(itTrack->dz());
580	<	h["trkEta"]->Fill(itTrack->eta());
581	<	h["trkPhi"]->Fill(itTrack->phi());
642	>	//========================================================================
643	>	// Step 1:  Apply event cleaning for data and MC
644	>	//          WARNING: event selection cut are hard coded!!
645	>	//========================================================================
646	>	// =====Loop over the trackColl to get the fraction of HighPurity tracks
647	>	int nTracks = 0;
648	>	int nHighPurityTracks=0;
649	>	double fHighPurity=0;
650	>
651	>	for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
652	>	TrackRef tkref(trackCollectionHandle,i);
653	>	if(tkref->quality(reco::TrackBase::highPurity)) ++nHighPurityTracks;
654		}
655	+	if(nTracks>0)
656	+	fHighPurity =  double(nHighPurityTracks)/double(nTracks);
657	+	if(verbose_)
658	+	cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;
659	+
660	+	if( (fHighPurity<0.2 && nTracks>10) || vertexColl->begin()->isFake()) return;
661	+
662	+	glob_runno_ = iEvent.id().run();
663	+	glob_evtno_ = iEvent.id().event();
664	+	glob_ls_   = iEvent.luminosityBlock();
665	+	glob_bx_  = iEvent.bunchCrossing();
666	+	return;
667
584	–	//SplittedTracks1
585	–	h["nTrk1_spl"]->Fill(splitTrackColl1->size());
586	–	for(TrackCollection::const_iterator itTrack = splitTrackColl1->begin();
587	–	itTrack != splitTrackColl1->end();
588	–	++itTrack) {
589	–	h["trkPt1_spl"]->Fill(itTrack->pt());
590	–	h["trkDxy1_spl"]->Fill(itTrack->dxy());
591	–	h["trkDz1_spl"]->Fill(itTrack->dz());
592	–	h["trkEta1_spl"]->Fill(itTrack->eta());
593	–	h["trkPhi1_spl"]->Fill(itTrack->phi());
594	–	}
595	–	//SplittedTracks2
596	–	h["nTrk2_spl"]->Fill(splitTrackColl2->size());
597	–	for(TrackCollection::const_iterator itTrack = splitTrackColl2->begin();
598	–	itTrack != splitTrackColl2->end();
599	–	++itTrack) {
600	–	h["trkPt2_spl"]->Fill(itTrack->pt());
601	–	h["trkDxy2_spl"]->Fill(itTrack->dxy());
602	–	h["trkDz2_spl"]->Fill(itTrack->dz());
603	–	h["trkEta2_spl"]->Fill(itTrack->eta());
604	–	h["trkPhi2_spl"]->Fill(itTrack->phi());
605	–	}
668
669	<	if(verbose_)
670	<	cout<<"Done filling track parameters of the input tracks to pvtx fitter."<<endl;
671	<
672	<	//=======================================================
673	<	// Fill reconstructed vertices
674	<	//=======================================================
675	<	if(verbose_)  {
676	<	cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
669	>	//========================================================================
670	>	// Step 2: Fill histograms for the splitting consistency checks
671	>	//========================================================================
672	>
673	>	// === Fill trackparameters of the input tracks to pvtx fitter
674	>	edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
675	>	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs)
676	>	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
677	>	fillTrackHisto(trackColl, 0, beamSpot);
678	>	fillTrackHisto(splitTrackColl1, 1, beamSpot);
679	>	fillTrackHisto(splitTrackColl2, 2, beamSpot);
680	>	edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
681	>
682	>
683	>	// ==== Fill number of reconstructed vertices
684	>	edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
685	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
686	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
687	>	edm::LogInfo("Debug")<<"[PVStudy] Start filling pvtx parameters."<<endl;
688	>	if (verbose_) {
689		printRecVtxs(vertexCollectionHandle);
616	–	cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
690		printRecVtxs(splitVertexCollection1Handle);
618	–	cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
691		printRecVtxs(splitVertexCollection2Handle);
692		}
693	+
694	+	nrecPV_ = int (vertexColl->size());
695	+	nrecPV1_spl_ = int (splitVertexColl1->size());
696	+	nrecPV2_spl_ = int (splitVertexColl2->size());
697	+
698	+	h_pvtrk->Fill1d("nrecPV", nrecPV_);
699	+	h_pvtrk->Fill1d("nrecPV1_spl", nrecPV1_spl_);
700	+	h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
701	+	h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
702
622	–	h["nrecPV"]->Fill(vertexColl->size());
623	–	h["nrecPV1_spl"]->Fill(splitVertexColl1->size());
624	–	h["nrecPV2_spl"]->Fill(splitVertexColl2->size());
625	–	h["nrecPVDiff"]->Fill(double(splitVertexColl1->size())-double(splitVertexColl2->size()));
626	–	//Fill the tree variables
627	–	nrecPV_ = vertexColl->size();
628	–	nrecPV1_spl_ = splitVertexColl1->size();
629	–	nrecPV2_spl_ = splitVertexColl2->size();
703
704	<	// Fill the track paramters for the vertexColl
705	<	for(reco::VertexCollection::const_iterator v=vertexColl->begin();
706	<	v!=vertexColl->end(); ++v){
707	<	try {
708	<	for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
709	<	t!=v->tracks_end(); t++) {
710	<	// illegal charge
711	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
712	<	//      h["trkPtPV"]->Fill(0.);
713	<	}
714	<	else {
715	<	h["trkPtPV"]->Fill((**t).pt());
716	<	h["trkDxyPV"]->Fill((**t).dxy());
717	<	h["trkDzPV"]->Fill((**t).dz());
718	<	h["trkEtaPV"]->Fill((**t).eta());
719	<	h["trkPhiPV"]->Fill((**t).phi());
720	<	}
704	>	// ======= Fill track parameter ntuple/hist for tracks used in recoVertices
705	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
706	>	if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
707	>	fillTrackHistoInPV(vertexColl, 0, beamSpot);
708	>
709	>	// ======= Fill secondary/primary min separations for multi-vertices
710	>	if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {
711	>	double min_xsep = 9999.0;
712	>	double min_ysep = 9999.0;
713	>	double min_zsep = 9999.0;
714	>	double min_xsep_sign = 9999.0;
715	>	double min_ysep_sign = 9999.0;
716	>	double min_zsep_sign = 9999.0;
717	>	double min_ntrksep = 9999.0;
718	>	double min_sumpt2sep = 9999999.0;
719	>
720	>	edm::LogInfo("Debug")<<"[PVStudy] leading pvtx z = "<< vertexColl->begin()->z() <<endl;
721	>
722	>	// Looping through the secondary vertices to find the mininum separation to the primary
723	>	for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
724	>	v!=vertexColl->end(); ++v) {
725	>	if(v->isValid() && ! v->isFake()) {
726	>	// xsep
727	>	if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
728	>	min_xsep = fabs(v->x()- vertexColl->begin()->x());
729	>	// ysep
730	>	if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
731	>	min_ysep = fabs(v->y()- vertexColl->begin()->y());
732	>	// zsep
733	>	if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
734	>	min_zsep = fabs(v->z()- vertexColl->begin()->z());
735	>	// xsep_sign
736	>	double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
737	>	if(xsep_sign < min_xsep_sign)
738	>	min_xsep_sign = xsep_sign;
739	>	// ysep_sign
740	>	double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
741	>	if(ysep_sign < min_ysep_sign)
742	>	min_ysep_sign = ysep_sign;
743	>	// zsep_sign
744	>	double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
745	>	if(zsep_sign < min_zsep_sign)
746	>	min_zsep_sign = zsep_sign;
747	>	// ntrksep
748	>	if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
749	>	min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
750	>	// sumpt2sep
751	>	if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
752	>	min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
753		}
754		}
755	<	catch (...) {
756	<	// exception thrown when trying to use linked track
757	<	//          h["trkPtPV"]->Fill(0.);
758	<	}
759	<	}
755	>	h_misc->Fill1d("min_xsep", min_xsep);
756	>	h_misc->Fill1d("min_ysep", min_ysep);
757	>	h_misc->Fill1d("min_zsep", min_zsep);
758	>	h_misc->Fill1d("min_xsep_sign", min_xsep_sign);
759	>	h_misc->Fill1d("min_ysep_sign", min_ysep_sign);
760	>	h_misc->Fill1d("min_zsep_sign", min_zsep_sign);
761	>	h_misc->Fill1d("min_ntrksep", min_ntrksep);
762	>	h_misc->Fill1d("min_sumpt2sep", min_sumpt2sep);
763	>
764	>	min_zsep_ = min_zsep;
765	>	min_ntrksep_ = min_ntrksep;
766	>	min_sumpt2sep_ = min_sumpt2sep;
767	>	} // End of  if(vertexColl->size()>1) {
768
769	+	edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
770	+
771	+	//========================================================================
772	+	// Step 3:  PrimaryVertex Matching
773	+	//   1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
774	+	//   2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
775	+	//   3. The first match is the primary vertex, which has largest sum(pT^2)
776	+	//========================================================================
777
778	<	// == ignore multi-vertices (for now)
779	<	//take single pvt collection and compare
780	<	if(splitVertexColl1->size() == 1 && splitVertexColl2->size() == 1 ) {
781	<	const reco::Vertex recvtx1 = *(splitVertexColl1->begin());
782	<	const reco::Vertex recvtx2 = *(splitVertexColl2->begin());
783	<	if (matchVertex(recvtx1, recvtx2))  {
784	<	int nTrkPV1 = recvtx1.tracksSize();
785	<	int nTrkPV2 = recvtx2.tracksSize();
786	<	h["nTrkPV1_spl"]->Fill(nTrkPV1);
787	<	h["nTrkPV2_spl"]->Fill(nTrkPV2);
788	<	double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
789	<	h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
790	<	h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
791	<
792	<	// Fill the track paramters for the splitVertexColl1
793	<	try {
794	<	for(reco::Vertex::trackRef_iterator t = recvtx1.tracks_begin();
795	<	t!=recvtx1.tracks_end(); t++) {
796	<	// illegal charge
797	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
798	<	//    h["trkPtPV"]->Fill(0.);
799	<	}
679	<	else {
680	<	h["trkPtPV1_spl"]->Fill((**t).pt());
681	<	h["trkDxyPV1_spl"]->Fill((**t).dxy());
682	<	h["trkDzPV1_spl"]->Fill((**t).dz());
683	<	h["trkEtaPV1_spl"]->Fill((**t).eta());
684	<	h["trkPhiPV1_spl"]->Fill((**t).phi());
685	<	}
686	<	}
687	<	}
688	<	catch (...) {
689	<	// exception thrown when trying to use linked track
690	<	//        h["trkPtPV"]->Fill(0.);
691	<	}
778	>	edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
779	>	reco::VertexCollection s_empty_matchedVertexColl1;
780	>	reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
781	>	matchedVertexColl1->reserve(splitVertexColl1->size());
782	>	reco::VertexCollection s_empty_matchedVertexColl2;
783	>	reco::VertexCollection *matchedVertexColl2 = &s_empty_matchedVertexColl2;
784	>	matchedVertexColl2->reserve(splitVertexColl2->size());
785	>
786	>	bool stopmatching = false;
787	>	for (size_t ivtx = 0; ivtx<splitVertexCollection1Handle->size() && !stopmatching; ++ivtx) {
788	>	reco::VertexRef recvtx1(splitVertexCollection1Handle, ivtx);
789	>	if( !(recvtx1->isValid()) || recvtx1->isFake()) break;
790	>	for (size_t jvtx = ivtx; jvtx < splitVertexCollection2Handle->size(); ++jvtx) {
791	>	//------------------------------------------------------------------------
792	>	//== If only considering matching the first vertex of the splitVertexColl
793	>	//------------------------------------------------------------------------
794	>	if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
795	>	reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
796	>	if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
797	>	edm::LogInfo("Debug")<<"[PVStudy] Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
798	>	edm::LogInfo("Debug")<<"[PVStudy] recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
799	>	edm::LogInfo("Debug")<<"[PVStudy] recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
800
801	<	// Fill the track paramters for the splitVertexColl2
802	<	try {
803	<	for(reco::Vertex::trackRef_iterator t = recvtx2.tracks_begin();
804	<	t!=recvtx2.tracks_end(); t++) {
805	<	// illegal charge
806	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
807	<	//    h["trkPtPV"]->Fill(0.);
808	<	}
809	<	else {
810	<	h["trkPtPV2_spl"]->Fill((**t).pt());
811	<	h["trkDxyPV2_spl"]->Fill((**t).dxy());
812	<	h["trkDzPV2_spl"]->Fill((**t).dz());
813	<	h["trkEtaPV2_spl"]->Fill((**t).eta());
814	<	h["trkPhiPV2_spl"]->Fill((**t).phi());
815	<	}
801	>	double twovtxsig = (recvtx1->z()-recvtx2->z())/max(recvtx1->zError(), recvtx2->zError());
802	>	h_misc->Fill1d("twovtxzsign", twovtxsig);
803	>	if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
804	>	edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
805	>
806	>	int nTrkPV1 = recvtx1->tracksSize();
807	>	int nTrkPV2 = recvtx2->tracksSize();
808	>	double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
809	>	h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
810	>	h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
811	>	if(fabs(ntrkreldiff)<ntrkdiffcut_) {
812	>	edm::LogInfo("Debug")<<"[PVStudy] (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
813	>
814	>	matchedVertexColl1->push_back(*recvtx1);
815	>	matchedVertexColl2->push_back(*recvtx2);
816	>	stopmatching = true; // stop the matching when the first match is found!
817	>	break;
818		}
819	+	else
820	+	edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
821		}
822	<	catch (...) {
823	<	// exception thrown when trying to use linked track
824	<	//        h["trkPtPV"]->Fill(0.);
822	>	else {
823	>	edm::LogInfo("Debug")<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
824	>	edm::LogInfo("Debug")<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
825	>	edm::LogInfo("Debug")<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
826		}
714	–
715	–	if(verbose_)
716	–	cout<<"Done filling track parameters in reconstruced vertices."<<endl;
717	–
718	–	if(fabs(ntrkreldiff)<ntrkdiffcut_) {
719	–	fres_[0] = (recvtx1.x()-recvtx2.x())/sqrt(2.0);
720	–	fres_[1] = (recvtx1.y()-recvtx2.y())/sqrt(2.0);
721	–	fres_[2] = (recvtx1.z()-recvtx2.z())/sqrt(2.0);
722	–	ferror_[0] = sqrt(pow(recvtx1.xError(),2)+pow(recvtx2.xError(),2))/sqrt(2);
723	–	ferror_[1] = sqrt(pow(recvtx1.yError(),2)+pow(recvtx2.yError(),2))/sqrt(2);
724	–	ferror_[2] = sqrt(pow(recvtx1.zError(),2)+pow(recvtx2.zError(),2))/sqrt(2);
725	–	fntrk_ =  (nTrkPV1+nTrkPV2)/2;
726	–
727	–	h["deltax"]->Fill( fres_[0] );
728	–	h["deltay"]->Fill( fres_[1] );
729	–	h["deltaz"]->Fill( fres_[2] );
730	–	h["pullx"]->Fill( fres_[0]/ferror_[0]);
731	–	h["pully"]->Fill( fres_[1]/ferror_[1]);
732	–	h["pullz"]->Fill( fres_[2]/ferror_[2]);
733	–	h["errPVx"]->Fill( ferror_[0] );
734	–	h["errPVy"]->Fill( ferror_[1] );
735	–	h["errPVz"]->Fill( ferror_[2] );
736	–	pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
737	–	error(ferror_[0], ferror_[1], ferror_[2]),
738	–	fntrk_));
739	–	// Fill histo according to its track multiplicity, set datamode = 0
740	–	fillHisto(res(fres_[0], fres_[1], fres_[2]),
741	–	error(ferror_[0], ferror_[1], ferror_[2]),
742	–	fntrk_,0);
743	–
744	–	ftree_->Fill();
745	–	} // End of  if(fabs(ntrkreldiff)<ntrkdiffcut_) {
746	–	} // End of  if (matchVertex(recvtx1, recvtx2))  {
747	–	else   {
748	–	cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
749	–	cout<<"recvtx1.z() = " << recvtx1.z() <<"; recvtx2.z() = " <<  recvtx2.z() <<endl;
827		}
828	<	} // End of fill histograms for two-vertices Method
828	>	}
829	>	edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
830
831	<	if(verbose_)
832	<	cout<<"Done filling res/pull with two-vertices method"<<endl;
833	<
834	<	//=======================================================
835	<	// Fill simulated vertices
836	<	//=======================================================
837	<	if (!realData_) {
838	<	bool MC=false;
839	<	Handle<HepMCProduct> evtMC;
840	<	iEvent.getByLabel("generator",evtMC);
841	<	if (!evtMC.isValid()) {
842	<	MC=false;
843	<	if(verbose_){
844	<	cout << "no HepMCProduct found"<< endl;
845	<	}
846	<	} else {
847	<	if(verbose_){
848	<	cout << "generator HepMCProduct found"<< endl;
849	<	}
850	<	MC=true;
831	>	//========================================================================
832	>	// Step 4: Analyze matchedVertexColls
833	>	//========================================================================
834	>	edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
835	>
836	>	// ==== If it is MC, analyze the res/pull of the unsplit vertexColl first
837	>	if(MC){
838	>	// make a list of primary vertices:
839	>	std::vector<simPrimaryVertex> simpv;
840	>	simpv=getSimPVs(evtMC,"");
841	>	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
842	>	h_gen->Fill1d("nsimPV", simpv.size());
843	>	nsimPV_ = simpv.size();
844	>	int isimpv = 0;
845	>	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
846	>	vsim!=simpv.end(); vsim++, isimpv++){
847	>	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
848	>	simx_[isimpv] = vsim->x;
849	>	simy_[isimpv] = vsim->y;
850	>	simz_[isimpv] = vsim->y;
851	>	simptsq_[isimpv] = vsim->ptsq;
852	>	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
853	>	fillMCHisto(vsim, isimpv, vertexColl, 3);
854		}
855	<
855	>	}
856	>
857	>	// Compare the reconstructed vertex position and calculate resolution/pulls
858	>	if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
859	>	// ==== Analyze the splitted vtx using MC method
860		if(MC){
861		// make a list of primary vertices:
862		std::vector<simPrimaryVertex> simpv;
863		simpv=getSimPVs(evtMC,"");
864		//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
865	<	h["nsimPV"]->Fill(simpv.size());
781	<	int nsimtrk=0;
782	<	int nrectrk=0;
783	<
865	>	int isimpv = 0;
866		for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
867	<	vsim!=simpv.end(); vsim++){
868	<	nsimtrk+=vsim->nGenTrk;
869	<	// Loop over the datamodes, for MC, fill only 1, 2, 3
870	<	for (int i=1;i<=3;i++) {
789	<	static const reco::VertexCollection s_empty_vtxColl;
790	<	const reco::VertexCollection *vtxColl = &s_empty_vtxColl;
791	<	std::string suffix;
792	<	// Set the vertexColl and histogram suffix according to datamode
793	<	if (i == 1) {
794	<	vtxColl = splitVertexColl1;
795	<	suffix = "_spl1_mct";
796	<	}
797	<	if (i == 2) {
798	<	vtxColl =  splitVertexColl2;
799	<	suffix = "_spl2_mct";
800	<	}
801	<	if (i == 3) {
802	<	vtxColl =  vertexColl;
803	<	suffix = "_mct";
804	<	}
805	<	//========================================================
806	<	//  look for a matching reconstructed vertex in vertexColl
807	<	//========================================================
808	<	vsim->recVtx=NULL;
809	<	for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
810	<	vrec!=vtxColl->end(); ++vrec){
811	<	nrectrk=vrec->tracksSize();
812	<	if(verbose_){
813	<	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
814	<	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
815	<	}
816	<	if ( matchVertex(*vsim,*vrec) ){
817	<	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
818	<	if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
819	<	|| (!vsim->recVtx) )
820	<	vsim->recVtx=&(*vrec);
821	<	}// End of finding the matched reco_vertex
822	<
823	<	if (vsim->recVtx){
824	<	if(verbose_) {
825	<	cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
826	<	}
827	<	double fres_mct[3];
828	<	double ferror_mct[3];
829	<	int fntrk_mct;
830	<
831	<	fres_mct[0] = vsim->recVtx->x()-vsim->x;
832	<	fres_mct[1] = vsim->recVtx->y()-vsim->y;
833	<	fres_mct[2] = vsim->recVtx->z()-vsim->z;
834	<	ferror_mct[0] = vsim->recVtx->xError();
835	<	ferror_mct[1] = vsim->recVtx->yError();
836	<	ferror_mct[2] = vsim->recVtx->zError();
837	<	fntrk_mct = nrectrk;
838	<
839	<	//Fill the values for variables in ftree_
840	<	if(i == 1) {
841	<	for(int j = 0;j<3;j++)
842	<	fres_spl1_mct_[j] = fres_mct[j];
843	<	fntrk_spl1_mct_ =  fntrk_mct;
844	<	}
845	<	if(i == 2) {
846	<	for(int j = 0;j<3;j++)
847	<	fres_spl2_mct_[j] = fres_mct[j];
848	<	fntrk_spl2_mct_ =  fntrk_mct;
849	<	}
850	<	if(i == 3) {
851	<	for(int j=0;j<3;j++)
852	<	fres_mct_[j] = fres_mct[j];
853	<	fntrk_mct_ =  fntrk_mct;
854	<	}
855	<
856	<	h["deltax"+suffix]->Fill( fres_mct[0] );
857	<	h["deltay"+suffix]->Fill( fres_mct[1] );
858	<	h["deltaz"+suffix]->Fill( fres_mct[2] );
859	<	h["pullx"+suffix]->Fill( fres_mct[0]/ferror_mct[0] );
860	<	h["pully"+suffix]->Fill( fres_mct[1]/ferror_mct[1] );
861	<	h["pullz"+suffix]->Fill( fres_mct[2]/ferror_mct[2] );
862	<	h["errPVx"+suffix]->Fill( ferror_mct[0] );
863	<	h["errPVy"+suffix]->Fill( ferror_mct[1] );
864	<	h["errPVz"+suffix]->Fill( ferror_mct[2] );
865	<	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
866	<	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
867	<	fntrk_mct));
868	<	// Fill histo according to its track multiplicity
869	<	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
870	<	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
871	<	fntrk_mct, i);
872	<	ftree_->Fill();
873	<	suffix.clear();
874	<	}
875	<	else {  // no rec vertex found for this simvertex
876	<	if(verbose_)
877	<	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
878	<	}
879	<	}
880	<	} // === End of Looping through vertexColl
867	>	vsim!=simpv.end(); vsim++, isimpv++){
868	>	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
869	>	fillMCHisto(vsim, isimpv, matchedVertexColl1, 1);
870	>	fillMCHisto(vsim, isimpv, matchedVertexColl2, 2);
871		}
872	<	} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
873	<	} // End of Filling MC variables
872	>	}
873	>
874	>	// ==== Analyze res/pull two-vertex method
875	>	//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype)
876	>	fillTrackHistoInPV(matchedVertexColl1, 1, beamSpot);
877	>	fillTrackHistoInPV(matchedVertexColl2, 2, beamSpot);
878	>
879	>	reco::VertexCollection::const_iterator v1;
880	>	reco::VertexCollection::const_iterator v2;
881	>
882	>	nrecPV_twovtx_ = 0;
883	>	for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
884	>	v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
885	>	++v1, ++v2) {
886	>	h_misc->Fill1d("ndofPVDiff", v1->ndof() - v2->ndof());
887	>	h_misc->Fill1d("ndofPVRelDiff", (v1->ndof()-v2->ndof())/(v1->ndof()+v2->ndof()));
888	>	fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
889	>	fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
890	>	fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
891	>	ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
892	>	ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
893	>	ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
894	>	fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
895	>
896	>	h_summary->Fill1d("deltax", fres_[0] );
897	>	h_summary->Fill1d("deltay", fres_[1] );
898	>	h_summary->Fill1d("deltaz", fres_[2] );
899	>	h_summary->Fill1d("pullx", fres_[0]/ferror_[0]);
900	>	h_summary->Fill1d("pully", fres_[1]/ferror_[1]);
901	>	h_summary->Fill1d("pullz", fres_[2]/ferror_[2]);
902	>	h_summary->Fill1d("errPVx", ferror_[0] );
903	>	h_summary->Fill1d("errPVy", ferror_[1] );
904	>	h_summary->Fill1d("errPVz", ferror_[2] );
905	>	pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
906	>	error(ferror_[0], ferror_[1], ferror_[2]),
907	>	fntrk_));
908	>	// Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
909	>	fillHisto(res(fres_[0], fres_[1], fres_[2]),
910	>	error(ferror_[0], ferror_[1], ferror_[2]),
911	>	fntrk_,0);
912	>
913	>	// Fill the ntuple variables
914	>	nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
915	>	deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
916	>	deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
917	>	deltaz_twovtx_[nrecPV_twovtx_] = fres_[2];
918	>	errx_twovtx_[nrecPV_twovtx_] = ferror_[0];
919	>	erry_twovtx_[nrecPV_twovtx_] = ferror_[1];
920	>	errz_twovtx_[nrecPV_twovtx_] = ferror_[2];
921	>	pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
922	>	pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
923	>	pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];
924	>
925	>
926	>	//SplittedVertex
927	>
928	>	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
929	>	ndofPV1_spl_twovtx_[nrecPV_twovtx_] = v1->ndof();
930	>	normchi2PV1_spl_twovtx_[nrecPV_twovtx_] = v1->normalizedChi2();
931	>	avgPtPV1_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v1);
932	>	errx1_spl_twovtx_[nrecPV_twovtx_] = v1->xError();
933	>	erry1_spl_twovtx_[nrecPV_twovtx_] = v1->yError();
934	>	errz1_spl_twovtx_[nrecPV_twovtx_] = v1->zError();
935	>
936	>	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
937	>	ndofPV2_spl_twovtx_[nrecPV_twovtx_] = v2->ndof();
938	>	normchi2PV2_spl_twovtx_[nrecPV_twovtx_] = v2->normalizedChi2();
939	>	avgPtPV2_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v2);
940	>	errx2_spl_twovtx_[nrecPV_twovtx_] = v2->xError();
941	>	erry2_spl_twovtx_[nrecPV_twovtx_] = v2->yError();
942	>	errz2_spl_twovtx_[nrecPV_twovtx_] = v2->zError();
943	>
944	>	nrecPV_twovtx_++;
945	>
946	>
947	>
948	>	// Print some information of the two tracks events
949	>	if(verbose_ && fntrk_ < 4) {
950	>	cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
951	>	printRecVtx(*v1);
952	>	cout<<"Printing matchedVertexColl1 for low ntrack bins"<<endl;
953	>	printRecVtx(*v2);
954	>	}
955	>	} // End of analyzing res/pull
956	>	} // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
957	>	else
958	>	edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
959	>
960	>	edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
961	>
962	>
963	>	//========================================================================
964	>	// Step 5: Fill ntuple if saventuple_ is on
965	>	//========================================================================
966	>	if(saventuple_) {
967	>	ftree_->Fill();
968	>	pvtxtree_->Fill();
969	>	}
970	>
971		}
972
973	+
974	+
975		// ------------ method called once each job just before starting event loop  ------------
976		void
977		PVStudy::beginJob()
#	Line 892 | Line 981 | PVStudy::beginJob()
981		// ------------ method called once each job just after ending the event loop  ------------
982		void
983		PVStudy::endJob() {
984	<	if (verbose_) cout << "Analyzing PV info" << endl;
984	>	edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl;
985		// Looping through the datamodes available
986		for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
987		string suffix;
988	<	if(verbose_)  cout<<"datamode = "<< *it<<endl;
988	>	edm::LogInfo("Analysis")<<"[endJob] datamode = "<< *it<<endl;
989		switch(*it) {
990		case 0: suffix = "";
991		break;
#	Line 907 | Line 996 | PVStudy::endJob() {
996		case 3: suffix = "_mct";
997		break;
998		}
999	<	for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1000	<	if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
1001	<	PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
1002	<	if (analyze_) {
1003	<	if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
1004	<	if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
1005	<	if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
1006	<	if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
1007	<	if ( ipv.pull_.y() > 0 ) h2["pully"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
1008	<	if ( ipv.pull_.z() > 0 ) h2["pullz"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.z(), 1.);
999	>	suffix += "_nTrk";
1000	>	if(analyze_) {
1001	>	for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1002	>	edm::LogInfo("Analysis") << "[endJob] ntrk = " << ntrk << endl;
1003	>	PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
1004	>	if ( ipv.res_.x() > 0 ) h_ana->Fill2d("resx"+suffix, ntrk,ipv.res_.x());
1005	>	if ( ipv.res_.y() > 0 ) h_ana->Fill2d("resy"+suffix, ntrk,ipv.res_.y());
1006	>	if ( ipv.res_.z() > 0 ) h_ana->Fill2d("resz"+suffix, ntrk,ipv.res_.z());
1007	>	if ( ipv.pull_.x() > 0 ) h_ana->Fill2d("pullx"+suffix, ntrk,ipv.pull_.x());
1008	>	if ( ipv.pull_.y() > 0 ) h_ana->Fill2d("pully"+suffix, ntrk,ipv.pull_.y());
1009	>	if ( ipv.pull_.z() > 0 ) h_ana->Fill2d("pullz"+suffix, ntrk,ipv.pull_.z());
1010		}
1011	<	} // End of  for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1011	>	}
1012		suffix.clear();
1013		}
1014	<	file_->cd();
1015	<	ftree_->Write();
1016	<	file_->Close();
1017	<	}
1018	<
1014	>	if(saventuple_) {
1015	>	file_->cd();
1016	>	ftree_->Write();
1017	>	pvtxtree_->Write();
1018	>	file_->Close();
1019	>	}
1020	>
1021
1022	+	}
1023
1024		// Match a recovertex with a simvertex
1025		// ? Should the cut be parameter dependent?
#	Line 934 | Line 1027 | PVStudy::endJob() {
1027		bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
1028		const reco::Vertex & vrec)
1029		{
1030	<	return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1030	>	if(vrec.isFake() || !vrec.isValid()) return false;
1031	>	else
1032	>	return true;
1033	>	//return (fabs(vsim.z-vrec.z())<0.0500); // =500um // remove this hard cut
1034		}
1035
1036		// Match two reconstructed vertices
1037	<	bool PVStudy::matchVertex( const reco::Vertex & vrec1,
1038	<	const reco::Vertex & vrec2)
1037	>	bool PVStudy::matchVertex( const reco::VertexRef & vrec1,
1038	>	const reco::VertexRef & vrec2,
1039	>	int zsigncut)
1040		{
1041	<	return (fabs(vrec1.z()-vrec2.z())<0.0500); // =500um
1041	>	double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
1042	>	edm::LogInfo("Debug")<<"[matchVertex] The matching criteria in z is "<<cut<<endl;
1043	>	return (fabs(vrec1->z()-vrec2->z())<cut);
1044		}
1045
1046
1047		bool PVStudy::isResonance(const HepMC::GenParticle * p){
1048		double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
1049	<	if(verbose_) cout << "isResonance   " << p->pdg_id() << " " << ctau << endl;
1049	>	edm::LogInfo("Debug") << "[isResonance] isResonance   " << p->pdg_id() << " " << ctau << endl;
1050		return  ctau >0 && ctau <1e-6;
1051		}
1052
#	Line 998 | Line 1097 | void PVStudy::printRecVtxs(const Handle<
1097		}
1098		}
1099
1100	+	void PVStudy::printRecVtx(const reco::Vertex & v){
1101	+
1102	+	cout << "#trk " << std::setw(3) << v.tracksSize()
1103	+	<< " chi2 " << std::setw(4) << v.chi2()
1104	+	<< " ndof " << std::setw(3) << v.ndof() << endl
1105	+	<< " x "  << std::setw(8) <<std::fixed << std::setprecision(4) << v.x()
1106	+	<< " dx " << std::setw(8) << v.xError()<< endl
1107	+	<< " y "  << std::setw(8) << v.y()
1108	+	<< " dy " << std::setw(8) << v.yError()<< endl
1109	+	<< " z "  << std::setw(8) << v.z()
1110	+	<< " dz " << std::setw(8) << v.zError()
1111	+	<< endl;
1112	+	}
1113	+
1114		void PVStudy::printSimTrks(const Handle<SimTrackContainer> simTrks){
1115		cout <<  " simTrks   type, (momentum), vertIndex, genpartIndex"  << endl;
1116		int i=1;
#	Line 1019 | Line 1132 | void PVStudy::printSimTrks(const Handle<
1132		void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1133		stringstream ss;
1134		ss << ntk;
1135	<	string suffix = ss.str();
1136	<	if(datamode == 0 )  suffix = "_" + suffix;
1137	<	if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1138	<	if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1139	<	if(datamode == 3 )  suffix = "_mct_" + suffix;
1135	>	string suffix;
1136	>	if(datamode == 0 )  suffix = "_" + ss.str();
1137	>	if(datamode == 1 )  suffix = "_spl1_mct_" + ss.str();
1138	>	if(datamode == 2 )  suffix = "_spl2_mct_" + ss.str();
1139	>	if(datamode == 3 )  suffix = "_mct_" + ss.str();
1140
1141		if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
1142		// Fill histograms of res and pull of ntk
1143	<	h["deltax"+suffix]->Fill(r.x());
1144	<	h["deltay"+suffix]->Fill(r.y());
1145	<	h["deltaz"+suffix]->Fill(r.z());
1146	<	h["pullx"+suffix]->Fill(r.x()/er.x());
1147	<	h["pully"+suffix]->Fill(r.y()/er.y());
1148	<	h["pullz"+suffix]->Fill(r.z()/er.z());
1149	<	h["errPVx"+suffix]->Fill(er.x());
1150	<	h["errPVy"+suffix]->Fill(er.y());
1151	<	h["errPVz"+suffix]->Fill(er.z());
1143	>	h_others->Fill1d("deltax"+suffix, r.x());
1144	>	h_others->Fill1d("deltay"+suffix, r.y());
1145	>	h_others->Fill1d("deltaz"+suffix, r.z());
1146	>	h_others->Fill1d("pullx"+suffix, r.x()/er.x());
1147	>	h_others->Fill1d("pully"+suffix, r.y()/er.y());
1148	>	h_others->Fill1d("pullz"+suffix, r.z()/er.z());
1149	>	h_others->Fill1d("errPVx"+suffix, er.x());
1150	>	h_others->Fill1d("errPVy"+suffix, er.y());
1151	>	h_others->Fill1d("errPVz"+suffix, er.z());
1152		}
1153
1154	+
1155		PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1156		map<int,double> data;
1157		data.clear();
#	Line 1055 | Line 1169 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1169		for ( int i = 0; i < 6; ++i) {
1170		switch (i) {
1171		case 0:
1172	<	fit(h["deltax"+suffix], fpar);
1172	>	if (!h_others->ReadHisto1D("deltax"+suffix)) break;
1173	>	fit(h_others->ReadHisto1D("deltax"+suffix), fpar);
1174		data[i] = fpar[0]*cm2um;
1175		break;
1176		case 1:
1177	<	fit(h["deltay"+suffix], fpar);
1177	>	if (!h_others->ReadHisto1D("deltay"+suffix)) break;
1178	>	fit(h_others->ReadHisto1D("deltay"+suffix), fpar);
1179		data[i] = fpar[0]*cm2um;
1180		break;
1181		case 2:
1182	<	fit(h["deltaz"+suffix], fpar);
1182	>	if (!h_others->ReadHisto1D("deltaz"+suffix)) break;
1183	>	fit(h_others->ReadHisto1D("deltaz"+suffix), fpar);
1184		data[i] = fpar[0]*cm2um;
1185		break;
1186		case 3:
1187	<	fit(h["pullx"+suffix], fpar);
1187	>	if (!h_others->ReadHisto1D("pullx"+suffix)) break;
1188	>	fit(h_others->ReadHisto1D("pullx"+suffix), fpar);
1189		data[i] = fpar[0];
1190		break;
1191		case 4:
1192	<	fit(h["pully"+suffix], fpar);
1192	>	if (!h_others->ReadHisto1D("pully"+suffix)) break;
1193	>	fit(h_others->ReadHisto1D("pully"+suffix), fpar);
1194		data[i] = fpar[0];
1195		break;
1196		case 5:
1197	<	fit(h["pullz"+suffix], fpar);
1197	>	if (!h_others->ReadHisto1D("pullz"+suffix)) break;
1198	>	fit(h_others->ReadHisto1D("pullz"+suffix), fpar);
1199		data[i] = fpar[0];
1200		break;
1201		}
1202	<	if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1202	>	if (data[i] > 0) edm::LogInfo("Analysis") << "[Analysis] data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1203		}
1204		}
1205		else
#	Line 1094 | Line 1214 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1214		ntk);
1215		}
1216
1217	<	void PVStudy::fit(TH1 *&hdist, double fitPars[]){
1217	>
1218	>	void PVStudy::fit(TH1 *hdist, double fitPars[]){
1219		TAxis *axis0 = hdist->GetXaxis();
1220		int nbin = axis0->GetLast();
1221		double nOF = hdist->GetBinContent(nbin+1);
1222		double nUF = hdist->GetBinContent(0);
1223		//   double fitRange = axis0->GetBinUpEdge(nbin);
1224	<	double fitRange = axis0->GetXmax();
1224	>	// double fitRange = axis0->GetXmax();
1225	>	double fitRange = 2*hdist->GetRMS();
1226		double sigMax[2] = {0.,0.};
1227
1228	<	if ( verbose_ ){
1229	<	cout << "Last bin = " << nbin;
1230	<	cout << "; hdist: Overflow Entries = " << nOF;
1231	<	cout << "; Underflow Entries = " << nUF;
1232	<	cout << "; hdist->GetEntries() = " << hdist->GetEntries();
1233	<	cout << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1112	<	}
1228	>	edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
1229	>	<< "; hdist: Overflow Entries = " << nOF
1230	>	<< "; Underflow Entries = " << nUF
1231	>	<< "; hdist->GetEntries() = " << hdist->GetEntries()
1232	>	<< "; fitting range = " << -fitRange << " to " << fitRange << endl;
1233	>
1234		if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit
1235		for (int bi = 0; bi < nbin; bi++) {
1236		if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) {
#	Line 1121 | Line 1242 | void PVStudy::fit(TH1 *&hdist, double fi
1242		if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1243		}
1244		}
1245	<	if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1245	>	//edm::LogInfo("Analysis") << "[Fit] Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1246
1247		TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1248		fgaus->SetParameter(1, 0.);
1249	<	fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
1249	>	fgaus->SetParLimits(1, -fitRange/10., fitRange/10.);
1250		fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
1251		fgaus->SetLineColor(4);
1252	<	hdist->Fit(fgaus,"Q");
1252	>	hdist->Fit(fgaus,"QLRM");
1253		gPad->Update();
1254		TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
1255		s->SetOptStat(1111111);
1256		s->SetOptFit(0111);
1257		gPad->Update();
1258	<	if (verbose_) cout << "got function" << endl;
1258	>	//edm::LogInfo("Analysis") << "[Fit] got function" << endl;
1259		fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1260		}
1261		else
1262		fitPars[0] = -999;
1263		}
1264
1265	+
1266	+	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs) {
1267	+	string suffix;
1268	+	suffix = ""; // for unsplit tracks
1269	+	if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1270	+	if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1271	+
1272	+	h_pvtrk->Fill1d("nTrk"+suffix, trackColl->size());
1273	+	for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1274	+	itTrack != trackColl->end();
1275	+	++itTrack) {
1276	+	h_pvtrk->Fill1d("trkPt"+suffix, itTrack->pt());
1277	+	h_pvtrk->Fill1d("trkDxy"+suffix, itTrack->dxy());
1278	+	h_pvtrk->Fill1d("trkDxyCorr"+suffix, itTrack->dxy(bs));
1279	+	h_pvtrk->Fill1d("trkDz"+suffix, itTrack->dz());
1280	+	h_pvtrk->Fill1d("trkEta"+suffix, itTrack->eta());
1281	+	h_pvtrk->Fill1d("trkPhi"+suffix, itTrack->phi());
1282	+	}
1283	+	}
1284	+
1285	+	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, const Point & bs) {
1286	+	string suffix;
1287	+	suffix = ""; // for unsplit tracks
1288	+	if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1289	+	if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1290	+	int ivtx = 0;
1291	+	for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1292	+	v!=vertexColl->end(); ++v, ++ivtx) {
1293	+	if(!v->isFake()) {
1294	+	h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1295	+	h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1296	+	int nHWTrkPV_ = 0;
1297	+	try {
1298	+	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1299	+	t!=vertexColl->begin()->tracks_end(); t++) {
1300	+	// illegal charge
1301	+	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1302	+	//    h_pvtrk->Fill1d("trkPtPV", 0.);
1303	+	}
1304	+	else {
1305	+	h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1306	+	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1307	+	h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
1308	+	h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1309	+	h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1310	+	h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1311	+	if(vertexColl->begin()->trackWeight(*t) > 0.5)
1312	+	nHWTrkPV_++;
1313	+	}
1314	+	}
1315	+	}
1316	+	catch (...) {
1317	+	// exception thrown when trying to use linked track
1318	+	//        h_pvtrk->Fill1d("trkPtPV", 0.);
1319	+	}
1320	+	h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
1321	+	}
1322	+	}
1323	+	}
1324	+
1325	+	void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
1326	+	{
1327	+	std::string suffix;
1328	+	// Set the vertexColl and histogram suffix according to datamode
1329	+	if (datatype == 1) {
1330	+	suffix = "_spl1_mct";
1331	+	nrecPV_spl1_mct_ = 0;
1332	+	}
1333	+	if (datatype == 2) {
1334	+	suffix = "_spl2_mct";
1335	+	nrecPV_spl2_mct_ = 0;
1336	+	}
1337	+	if (datatype == 3) {
1338	+	suffix = "_mct";
1339	+	nrecPV_mct_ = 0;
1340	+	}
1341	+
1342	+	//========================================================
1343	+	//  look for a matching reconstructed vertex in vertexColl
1344	+	//========================================================
1345	+
1346	+	for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1347	+	vrec!=vtxColl->end(); ++vrec){
1348	+	vsim->recVtx=NULL;
1349	+	edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1350	+	edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1351	+
1352	+	if ( matchVertex(*vsim,*vrec)  ) {
1353	+	vsim->recVtx=&(*vrec);
1354	+
1355	+	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1356	+	//if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
1357	+	//|| (!vsim->recVtx) )
1358	+	//vsim->recVtx=&(*vrec);
1359	+	//}
1360	+	edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1361	+
1362	+	double fres_mct[3];
1363	+	double ferror_mct[3];
1364	+	int nrectrk =  int(vsim->recVtx->tracksSize());
1365	+	double ndofPV = vsim->recVtx->ndof();
1366	+	double normchi2PV =  vsim->recVtx->normalizedChi2();
1367	+	double avgPtPV = avgPtRecVtx(*(vsim->recVtx));
1368	+
1369	+	fres_mct[0] = vsim->recVtx->x()-vsim->x;
1370	+	fres_mct[1] = vsim->recVtx->y()-vsim->y;
1371	+	fres_mct[2] = vsim->recVtx->z()-vsim->z;
1372	+	ferror_mct[0] = vsim->recVtx->xError();
1373	+	ferror_mct[1] = vsim->recVtx->yError();
1374	+	ferror_mct[2] = vsim->recVtx->zError();
1375	+
1376	+	h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1377	+	h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1378	+	h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1379	+	h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1380	+	h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1381	+	h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1382	+	h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1383	+	h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1384	+	h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1385	+	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1386	+	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1387	+	nrectrk));
1388	+	// Fill histo according to its track multiplicity
1389	+	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1390	+	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1391	+	nrectrk, datatype);
1392	+
1393	+	if(saventuple_) {
1394	+	//Fill the values for variables in pvtxtree_
1395	+	if(datatype == 1) {
1396	+	nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1397	+	ndofPV_spl1_mct_[nrecPV_spl1_mct_] =  ndofPV;
1398	+	normchi2PV_spl1_mct_[nrecPV_spl1_mct_] =  normchi2PV;
1399	+	avgPtPV_spl1_mct_[nrecPV_spl1_mct_] = avgPtPV;
1400	+	deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1401	+	deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[1];
1402	+	deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2];
1403	+	pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
1404	+	pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1405	+	pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1406	+	errx_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[0];
1407	+	erry_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[1];
1408	+	errz_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[2];
1409	+	nrecPV_spl1_mct_++;
1410	+	}
1411	+	if(datatype == 2) {
1412	+	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1413	+	ndofPV_spl2_mct_[nrecPV_spl2_mct_] =  ndofPV;
1414	+	normchi2PV_spl2_mct_[nrecPV_spl2_mct_] =  normchi2PV;
1415	+	avgPtPV_spl2_mct_[nrecPV_spl2_mct_] = avgPtPV;
1416	+	deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1417	+	deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[1];
1418	+	deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2];
1419	+	pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
1420	+	pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1421	+	pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1422	+	errx_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[0];
1423	+	erry_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[1];
1424	+	errz_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[2];
1425	+	nrecPV_spl2_mct_++;
1426	+	}
1427	+	if(datatype == 3) {
1428	+	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1429	+	ndofPV_mct_[nrecPV_mct_] =  ndofPV;
1430	+	normchi2PV_mct_[nrecPV_mct_] =  normchi2PV;
1431	+	avgPtPV_mct_[nrecPV_mct_] = avgPtPV;
1432	+	deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1433	+	deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1434	+	deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1435	+	pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
1436	+	pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1437	+	pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1438	+	errx_mct_[nrecPV_mct_] =  ferror_mct[0];
1439	+	erry_mct_[nrecPV_mct_] =  ferror_mct[1];
1440	+	errz_mct_[nrecPV_mct_] =  ferror_mct[2];
1441	+	nrecPV_mct_++;
1442	+	}
1443	+	} // End of  if(saventuple_) {
1444	+	} //  if ( matchVertex(*vsim,*vrec) ) {
1445	+	else {  // no rec vertex found for this simvertex
1446	+	edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1447	+	}
1448	+	}
1449	+	}
1450	+
1451		void PVStudy::SetVarToZero() {
1452	+	fntrk_ = 0;
1453	+	//pvtx position (x,y,z) residual and error
1454	+	for(int i = 0; i<3;i++)    {
1455	+	fres_[i] = 0;
1456	+	ferror_[i] = 0;
1457	+	}
1458
1459	<	// number of reconstructed vertices
1459	>	// Event level information
1460	>	glob_runno_ = 0;
1461	>	glob_evtno_ = 0;
1462	>	glob_ls_ = 0;
1463	>	glob_bx_ = 0;
1464		nrecPV_ = 0;
1465	<	nrecPV1_spl_ = 0;
1465	>	nrecPV1_spl_ = 0;
1466		nrecPV2_spl_ = 0;
1467	<	// number of tracks in the vertex
1468	<	fntrk_ = 0;
1469	<	fntrk_spl1_mct_ = 0;
1470	<	fntrk_spl2_mct_ = 0;
1471	<	fntrk_mct_ = 0;
1472	<	//pvtx position (x,y,z) residual and error
1473	<	for(int i = 0; i<3;i++)
1474	<	{
1475	<	fres_[i] = 0;
1476	<	ferror_[i] = 0;
1477	<	fres_spl1_mct_[i] = 0;
1478	<	ferror_spl1_mct_[i] = 0;
1479	<	fres_spl2_mct_[i] = 0;
1480	<	ferror_spl2_mct_[i] = 0;
1481	<	fres_mct_[i] = 0;
1482	<	ferror_mct_[i] = 0;
1467	>	nrecPV_twovtx_ = 0;
1468	>	nsimPV_ = 0;
1469	>
1470	>	nrecPV_mct_ = 0;
1471	>	nrecPV_spl1_mct_ = 0;
1472	>	nrecPV_spl2_mct_ = 0;
1473	>
1474	>	// Mininum separation between the secondary pvtxes and leading pvtx
1475	>	min_zsep_ = 9999.0;
1476	>	min_ntrksep_ = 9999.0;
1477	>	min_sumpt2sep_ = -9999.0;
1478	>
1479	>	// Variables filled per Vertex
1480	>	for (int i = 0; i < nMaxPVs_; i++) {
1481	>	// matched two-vertices
1482	>	nTrkPV1_spl_twovtx_[i] = 0;
1483	>	ndofPV1_spl_twovtx_[i] = 0;
1484	>	normchi2PV1_spl_twovtx_[i] = 0;
1485	>	avgPtPV1_spl_twovtx_[i] = 0;
1486	>	errx1_spl_twovtx_[i] = 0;
1487	>	erry1_spl_twovtx_[i] = 0;
1488	>	errz1_spl_twovtx_[i] = 0;
1489	>
1490	>	nTrkPV2_spl_twovtx_[i] = 0;
1491	>	ndofPV2_spl_twovtx_[i] = 0;
1492	>	normchi2PV2_spl_twovtx_[i] = 0;
1493	>	avgPtPV2_spl_twovtx_[i] = 0;
1494	>	errx2_spl_twovtx_[i] = 0;
1495	>	erry2_spl_twovtx_[i] = 0;
1496	>	errz2_spl_twovtx_[i] = 0;
1497	>
1498	>	deltax_twovtx_[i] = 0;
1499	>	deltay_twovtx_[i] = 0;
1500	>	deltaz_twovtx_[i] = 0;
1501	>	errx_twovtx_[i] = 0;
1502	>	erry_twovtx_[i] = 0;
1503	>	errz_twovtx_[i] = 0;
1504	>	pullx_twovtx_[i] = 0;
1505	>	pully_twovtx_[i] = 0;
1506	>	pullz_twovtx_[i] = 0;
1507	>
1508	>	//simpv
1509	>	nsimTrkPV_[i] = 0;
1510	>	simx_[i] = 0;
1511	>	simy_[i] = 0;
1512	>	simz_[i] = 0;
1513	>	simptsq_[i] = 0;
1514	>
1515	>	// residual and pulls with mc truth required
1516	>	deltax_mct_[i] = 0;
1517	>	deltay_mct_[i] = 0;
1518	>	deltaz_mct_[i] = 0;
1519	>	pullx_mct_[i] = 0;
1520	>	pully_mct_[i] = 0;
1521	>	pullz_mct_[i] = 0;
1522	>	nTrkPV_mct_[i] = 0;
1523	>	ndofPV_mct_[i] = 0;
1524	>	normchi2PV_mct_[i] = 0;
1525	>	avgPtPV_mct_[i] = 0;
1526	>	errx_mct_[i] = 0;
1527	>	erry_mct_[i] = 0;
1528	>	errz_mct_[i] = 0;
1529	>
1530	>	deltax_spl1_mct_[i] = 0;
1531	>	deltay_spl1_mct_[i] = 0;
1532	>	deltaz_spl1_mct_[i] = 0;
1533	>	pullx_spl1_mct_[i] = 0;
1534	>	pully_spl1_mct_[i] = 0;
1535	>	pullz_spl1_mct_[i] = 0;
1536	>	nTrkPV_spl1_mct_[i] = 0;
1537	>	ndofPV_spl1_mct_[i] = 0;
1538	>	normchi2PV_spl1_mct_[i] = 0;
1539	>	avgPtPV_spl1_mct_[i] = 0;
1540	>	errx_spl1_mct_[i] = 0;
1541	>	erry_spl1_mct_[i] = 0;
1542	>	errz_spl1_mct_[i] = 0;
1543	>
1544	>	deltax_spl2_mct_[i] = 0;
1545	>	deltay_spl2_mct_[i] = 0;
1546	>	deltaz_spl2_mct_[i] = 0;
1547	>	pullx_spl2_mct_[i] = 0;
1548	>	pully_spl2_mct_[i] = 0;
1549	>	pullz_spl2_mct_[i] = 0;
1550	>	nTrkPV_spl2_mct_[i] = 0;
1551	>	ndofPV_spl2_mct_[i] = 0;
1552	>	normchi2PV_spl2_mct_[i] = 0;
1553	>	avgPtPV_spl2_mct_[i] = 0;
1554	>	errx_spl2_mct_[i] = 0;
1555	>	erry_spl2_mct_[i] = 0;
1556	>	errz_spl2_mct_[i] = 0;
1557	>
1558	>	}
1559	>	}
1560	>
1561	>	double PVStudy::sumPtSquared(const reco::Vertex & v)  {
1562	>	if(v.isFake() || (!v.isValid()) || v.tracksSize() == 0) return 0.0;
1563	>	else {
1564	>	double sum = 0.;
1565	>	double pT;
1566	>	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1567	>	pT = (**it).pt();
1568	>	sum += pT*pT;
1569		}
1570	+	return sum;
1571	+	}
1572	+	}
1573	+
1574	+
1575	+
1576	+	double PVStudy::avgPtRecVtx(const reco::Vertex & v)  {
1577	+
1578	+	if(v.isFake() || !v.isValid() || v.tracksSize()==0 ) return 0;
1579	+	else {
1580	+	double sumpT = 0.;
1581	+	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1582	+	sumpT +=  (**it).pt();
1583	+	}
1584	+	return sumpT/double(v.tracksSize());
1585	+	}
1586		}
1587
1588
1589	+

Diff Legend

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