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.5 by jengbou, Wed Sep 9 05:56:50 2009 UTC vs.
Revision 1.18 by jengbou, Fri Jan 29 21:36:22 2010 UTC

#	Line 12 | Line 12 | Implementation:
12		*/
13		//
14		// Original Author:  "Geng-yuan Jeng/UC Riverside"
15	+	//                    "Yanyan Gao/Fermilab ygao@fnal.gov"
16		//         Created:  Thu Aug 20 11:55:40 CDT 2009
17		// $Id$
18		//
#	Line 48 | 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 58 | 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):
63	<	nTrkMin_(10),nTrkMax_(100)
68	>	PVStudy::PVStudy(const edm::ParameterSet& iConfig)
69		{
70	<	simG4_           = iConfig.getParameter<edm::InputTag>( "simG4" );
71	<	tracksLabel_     = iConfig.getUntrackedParameter<edm::InputTag>("tracks");
72	<	vtxSample_       = iConfig.getUntrackedParameter<edm::InputTag>("vtxSample");
73	<	verbose_         = iConfig.getUntrackedParameter<bool>("verbose", false);
74	<	realData_        = iConfig.getUntrackedParameter<bool>("realData",false);
75	<	analyze_         = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
76	<	outputfilename_  = iConfig.getUntrackedParameter<string>("OutputFileName");
77	<
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");
77	>	vertexCollectionTag_       = iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection");
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
96	<	file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
97	<	ftree_ = new TTree("mytree","mytree");
98	<	ftree_->AutoSave();
99	<	ftree_->Branch("residual",&fres,"fres[3]/D");
100	<	ftree_->Branch("error",&ferror,"ferror[3]/D");
101	<	ftree_->Branch("nTrk",&fntrk,"fntrk/I");
102	<
103	<	edm::Service<TFileService> fs;
104	<	TFileDirectory subDir = fs->mkdir( "Summary" );
105	<	TFileDirectory subDir1 = fs->mkdir( "Others" );
106	<	//   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
107	<
96	>	// Specify the data mode vector
97	>	if(realData_) datamode.push_back(0);
98	>	else {
99	>	datamode.push_back(0);
100	>	datamode.push_back(1);
101	>	datamode.push_back(2);
102	>	datamode.push_back(3);
103	>	}
104	>	// Create ntuple files if needed
105	>	if(saventuple_) {
106	>	file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
107	>	ftree_ = new TTree("mytree","mytree");
108	>	ftree_->AutoSave();
109	>	ftree_->Branch("residual",&fres_,"fres_[3]/D");
110	>	ftree_->Branch("error",&ferror_,"ferror_[3]/D");
111	>	ftree_->Branch("nTrkPV",&fntrk_,"fntrk_/I");
112	>
113	>	// pvtxtree_ analyzing the pvtxs ootb
114	>	pvtxtree_ = new TTree("pvtxtree","pvtxtree");
115	>	//pvtx using all tracks
116	>	pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
117	>	pvtxtree_->Branch("nTrkPV",&nTrkPV_,"nTrkPV[nrecPV]/I");
118	>	pvtxtree_->Branch("sumptsq",&sumptsq_,"sumptsq[nrecPV]/D");
119	>	pvtxtree_->Branch("isValid",&isValid_,"isValid/I");
120	>	pvtxtree_->Branch("isFake",&isFake_,"isFake/I");
121	>	pvtxtree_->Branch("recx",&recx_,"recx[nrecPV]/D");
122	>	pvtxtree_->Branch("recy",&recy_,"recy[nrecPV]/D");
123	>	pvtxtree_->Branch("recz",&recz_,"recz[nrecPV]/D");
124	>	pvtxtree_->Branch("recx_err",&recx_err_,"recx_err[nrecPV]/D");
125	>	pvtxtree_->Branch("recy_err",&recy_err_,"recy_err[nrecPV]/D");
126	>	pvtxtree_->Branch("recz_err",&recz_err_,"recz_err[nrecPV]/D");
127	>
128	>	pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
129	>	pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
130	>	pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
131	>
132	>	//pvtx using splittracks1
133	>	pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
134	>	pvtxtree_->Branch("nTrkPV1_spl",&nTrkPV1_spl_,"nTrkPV1_spl[nrecPV1_spl]/I");
135	>	pvtxtree_->Branch("sumptsq1_spl",&sumptsq1_spl_,"sumptsq1_spl[nrecPV1_spl]/D");
136	>	pvtxtree_->Branch("isValid1_spl",&isValid1_spl_,"isValid1_spl/I");
137	>	pvtxtree_->Branch("isFake1_spl",&isFake1_spl_,"isFake1_spl/I");
138	>	pvtxtree_->Branch("recx1_spl",&recx1_spl_,"recx1_spl[nrecPV1_spl]/D");
139	>	pvtxtree_->Branch("recy1_spl",&recy1_spl_,"recy1_spl[nrecPV1_spl]/D");
140	>	pvtxtree_->Branch("recz1_spl",&recz1_spl_,"recz1_spl[nrecPV1_spl]/D");
141	>	pvtxtree_->Branch("recx1_err_spl",&recx1_err_spl_,"recx1_err_spl[nrecPV1_spl]/D");
142	>	pvtxtree_->Branch("recy1_err_spl",&recy1_err_spl_,"recy1_err_spl[nrecPV1_spl]/D");
143	>	pvtxtree_->Branch("recz1_err_spl",&recz1_err_spl_,"recz1_err_spl[nrecPV1_spl]/D");
144	>
145	>	//pvtx using splittracks2
146	>	pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
147	>	pvtxtree_->Branch("nTrkPV2_spl",&nTrkPV2_spl_,"nTrkPV2_spl[nrecPV2_spl]/I");
148	>	pvtxtree_->Branch("sumptsq2_spl",&sumptsq2_spl_,"sumptsq2_spl[nrecPV2_spl]/D");
149	>	pvtxtree_->Branch("isValid2_spl",&isValid2_spl_,"isValid2_spl/I");
150	>	pvtxtree_->Branch("isFake2_spl",&isFake2_spl_,"isFake2_spl/I");
151	>	pvtxtree_->Branch("recx2_spl",&recx2_spl_,"recx2_spl[nrecPV2_spl]/D");
152	>	pvtxtree_->Branch("recy2_spl",&recy2_spl_,"recy2_spl[nrecPV2_spl]/D");
153	>	pvtxtree_->Branch("recz2_spl",&recz2_spl_,"recz2_spl[nrecPV2_spl]/D");
154	>	pvtxtree_->Branch("recx2_err_spl",&recx2_err_spl_,"recx2_err_spl[nrecPV2_spl]/D");
155	>	pvtxtree_->Branch("recy2_err_spl",&recy2_err_spl_,"recy2_err_spl[nrecPV2_spl]/D");
156	>	pvtxtree_->Branch("recz2_err_spl",&recz2_err_spl_,"recz2_err_spl[nrecPV2_spl]/D");
157	>
158	>	//pixeVertices
159	>	pvtxtree_->Branch("nrecPV_pxlpvtx",&nrecPV_pxlpvtx_,"nrecPV_pxlpvtx/I");
160	>	pvtxtree_->Branch("nTrkPV_pxlpvtx",&nTrkPV_pxlpvtx_,"nTrkPV_pxlpvtx[nrecPV_pxlpvtx]/I");
161	>	pvtxtree_->Branch("sumptsq_pxlpvtx",&sumptsq_pxlpvtx_,"sumptsq_pxlpvtx[nrecPV_pxlpvtx]/D");
162	>	pvtxtree_->Branch("isValid_pxlpvtx",&isValid_pxlpvtx_,"isValid_pxlpvtx/I");
163	>	pvtxtree_->Branch("isFake_pxlpvtx",&isFake_pxlpvtx_,"isFake_pxlpvtx/I");
164	>	pvtxtree_->Branch("recx_pxlpvtx",&recx_pxlpvtx_,"recx_pxlpvtx[nrecPV_pxlpvtx]/D");
165	>	pvtxtree_->Branch("recy_pxlpvtx",&recy_pxlpvtx_,"recy_pxlpvtx[nrecPV_pxlpvtx]/D");
166	>	pvtxtree_->Branch("recz_pxlpvtx",&recz_pxlpvtx_,"recz_pxlpvtx[nrecPV_pxlpvtx]/D");
167	>	pvtxtree_->Branch("recx_err_pxlpvtx",&recx_err_pxlpvtx_,"recx_err_pxlpvtx[nrecPV_pxlpvtx]/D");
168	>	pvtxtree_->Branch("recy_err_pxlpvtx",&recy_err_pxlpvtx_,"recy_err_pxlpvtx[nrecPV_pxlpvtx]/D");
169	>	pvtxtree_->Branch("recz_err_pxlpvtx",&recz_err_pxlpvtx_,"recz_err_pxlpvtx[nrecPV_pxlpvtx]/D");
170	>
171	>	//Fill the variables in the twovtx pair (recvtx1, recvtx2)
172	>	pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
173	>	pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
174	>	pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
175	>	pvtxtree_->Branch("sumptsq1_spl_twovtx",&sumptsq1_spl_twovtx_,"sumptsq1_spl_twovtx[nrecPV_twovtx]/D");
176	>	pvtxtree_->Branch("sumptsq2_spl_twovtx",&sumptsq2_spl_twovtx_,"sumptsq2_spl_twovtx[nrecPV_twovtx]/D");
177	>	pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
178	>	pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
179	>	pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
180	>	pvtxtree_->Branch("deltaz_twovtx",&deltaz_twovtx_,"deltaz_twovtx[nrecPV_twovtx]/D");
181	>	pvtxtree_->Branch("errx_twovtx",&errx_twovtx_,"errx_twovtx[nrecPV_twovtx]/D");
182	>	pvtxtree_->Branch("erry_twovtx",&erry_twovtx_,"erry_twovtx[nrecPV_twovtx]/D");
183	>	pvtxtree_->Branch("errz_twovtx",&errz_twovtx_,"errz_twovtx[nrecPV_twovtx]/D");
184	>	pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
185	>	pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
186	>	pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
187	>
188	>	// MC variables
189	>	if(!realData_) {
190	>	pvtxtree_->Branch("nsimPV",&nsimPV_,"nsimPV/I");
191	>	pvtxtree_->Branch("nsimTrkPV",&nsimTrkPV_,"nsimTrkPV[nsimPV]/I");
192	>	pvtxtree_->Branch("simx",&simx_,"simx[nsimPV]/D");
193	>	pvtxtree_->Branch("simy",&simy_,"simy[nsimPV]/D");
194	>	pvtxtree_->Branch("simz",&simz_,"simz[nsimPV]/D");
195	>	pvtxtree_->Branch("simptsq",&simptsq_,"simptsq[nsimPV]/D");
196	>
197	>	// For pvtxs with all the tracks
198	>	pvtxtree_->Branch("nrecPV_mct",&nrecPV_mct_,"nrecPV_mct/I");
199	>	pvtxtree_->Branch("deltax_mct",&deltax_mct_,"deltax_mct[nrecPV_mct]/D");
200	>	pvtxtree_->Branch("deltay_mct",&deltay_mct_,"deltay_mct[nrecPV_mct]/D");
201	>	pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
202	>	pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
203	>	pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
204	>	pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
205	>	// For pvtxs with splittracks1
206	>	pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
207	>	pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
208	>	pvtxtree_->Branch("deltay_spl1_mct",&deltay_spl1_mct_,"deltay_spl1_mct[nrecPV_spl1_mct]/D");
209	>	pvtxtree_->Branch("deltaz_spl1_mct",&deltaz_spl1_mct_,"deltaz_spl1_mct[nrecPV_spl1_mct]/D");
210	>	pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
211	>	pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
212	>	pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
213	>	// For pvtxs with splittracks1
214	>	pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
215	>	pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
216	>	pvtxtree_->Branch("deltay_spl2_mct",&deltay_spl2_mct_,"deltay_spl2_mct[nrecPV_spl2_mct]/D");
217	>	pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
218	>	pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
219	>	pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
220	>	pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
221	>	}
222	>	}
223	>
224		setRootStyle();
225
226	<	h["nTrk"]           = subDir.make<TH1F>("nTrk","Num of total tracks",300,0,300);
227	<	h["nTrkPV"]         = subDir.make<TH1F>("nTrkPV","Num of Tracks from PV",300,0,300);
228	<	h["trkPt"]          = subDir.make<TH1D>("trkPt","Pt of all tracks",100,0,100);
229	<	h["trkPtPV"]        = subDir.make<TH1D>("trkPtPV","Pt dist of tracks from PV",100,0,100);
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	+	}
257	+	h_pixvtx->Init("pixVtx");
258	+	h_misc->Init("misc");
259	+
260	+	// Book MC only plots
261		if (!realData_) {
262	<	h["deltax"]       = subDir.make<TH1D>("deltax","x-distance (Rec - Sim)",800,-0.04,0.04);
263	<	h["deltay"]       = subDir.make<TH1D>("deltay","y-distance (Rec - Sim)",800,-0.04,0.04);
264	<	h["deltaz"]       = subDir.make<TH1D>("deltaz","z-distance (Rec - Sim)",800,-0.04,0.04);
265	<	h["pullx"]        = subDir.make<TH1D>("pullx","x-pull (Rec - Sim)",1000,-25.,25.);
266	<	h["pully"]        = subDir.make<TH1D>("pully","y-pull (Rec - Sim)",1000,-25.,25.);
267	<	h["pullz"]        = subDir.make<TH1D>("pullz","z-pull (Rec - Sim)",1000,-25.,25.);
268	<	h["errPVx"]       = subDir.make<TH1D>("errPVx","X vertex error",100,0.,0.02);
269	<	h["errPVy"]       = subDir.make<TH1D>("errPVy","Y vertex error",100,0.,0.02);
270	<	h["errPVz"]       = subDir.make<TH1D>("errPVz","Z vertex error",100,0.,0.02);
271	<
272	<	// Summary of analysis:
273	<	if (analyze_) {
274	<	h2["resx_nTrk"]= subDir.make<TH2D>("resx_ntrk","x-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
275	<	h2["resx_nTrk"]->SetMarkerStyle(21);
276	<	h2["resx_nTrk"]->SetMarkerColor(4);
277	<	h2["resx_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
278	<	h2["resx_nTrk"]->GetYaxis()->SetTitle("#mum");
279	<	h2["resy_nTrk"]= subDir.make<TH2D>("resy_ntrk","y-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
114	<	h2["resy_nTrk"]->SetMarkerStyle(21);
115	<	h2["resy_nTrk"]->SetMarkerColor(4);
116	<	h2["resy_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
117	<	h2["resy_nTrk"]->GetYaxis()->SetTitle("#mum");
118	<	h2["resz_nTrk"]= subDir.make<TH2D>("resz_ntrk","z-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
119	<	h2["resz_nTrk"]->SetMarkerStyle(21);
120	<	h2["resz_nTrk"]->SetMarkerColor(4);
121	<	h2["resz_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
122	<	h2["resz_nTrk"]->GetYaxis()->SetTitle("#mum");
123	<	h2["pullx_nTrk"]= subDir.make<TH2D>("pullx_ntrk","x-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
124	<	h2["pullx_nTrk"]->SetMarkerStyle(21);
125	<	h2["pullx_nTrk"]->SetMarkerColor(4);
126	<	h2["pullx_nTrk"]->SetBit(TH1::kCanRebin);
127	<	h2["pullx_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
128	<	h2["pully_nTrk"]= subDir.make<TH2D>("pully_ntrk","y-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
129	<	h2["pully_nTrk"]->SetMarkerStyle(21);
130	<	h2["pully_nTrk"]->SetMarkerColor(4);
131	<	h2["pully_nTrk"]->SetBit(TH1::kCanRebin);
132	<	h2["pully_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
133	<	h2["pullz_nTrk"]= subDir.make<TH2D>("pullz_ntrk","z-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
134	<	h2["pullz_nTrk"]->SetMarkerStyle(21);
135	<	h2["pullz_nTrk"]->SetMarkerColor(4);
136	<	h2["pullz_nTrk"]->SetBit(TH1::kCanRebin);
137	<	h2["pullz_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
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	>	edm::LogInfo("Debug")<<"datamode = "<< *it<<endl;
271	>	switch(*it) {
272	>	case 0: suffix = "";
273	>	break;
274	>	case 1: suffix = "_spl1_mct";
275	>	break;
276	>	case 2: suffix = "_spl2_mct";
277	>	break;
278	>	case 3: suffix = "_mct";
279	>	break;
280		}
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 suffix = ss.str();
143	<	h["resx_"+suffix]  = fs->make<TH1D>(TString("deltax_"+suffix),"x-distance (Rec - Sim)",100,-0.02,0.02);
144	<	h["resx_"+suffix]->GetXaxis()->SetTitle("cm");
145	<	h["resy_"+suffix]  = fs->make<TH1D>(TString("deltay_"+suffix),"y-distance (Rec - Sim)",100,-0.02,0.02);
146	<	h["resy_"+suffix]->GetXaxis()->SetTitle("cm");
147	<	h["resz_"+suffix]  = fs->make<TH1D>(TString("deltaz_"+suffix),"z-distance (Rec - Sim)",100,-0.02,0.02);
148	<	h["resz_"+suffix]->GetXaxis()->SetTitle("cm");
149	<	h["pullx_"+suffix] = fs->make<TH1D>(TString("pullx_"+suffix),"x-pull (Rec - Sim)",100,-10.,10.);
150	<	h["pully_"+suffix] = fs->make<TH1D>(TString("pully_"+suffix),"y-pull (Rec - Sim)",100,-10.,10.);
151	<	h["pullz_"+suffix] = fs->make<TH1D>(TString("pullz_"+suffix),"z-pull (Rec - Sim)",100,-10.,10.);
152	<	h["errPVx_"+suffix]  = subDir1.make<TH1D>(TString("errPVx_"+suffix),"X vertex error",100,0.,0.02);
153	<	h["errPVy_"+suffix]  = subDir1.make<TH1D>(TString("errPVy_"+suffix),"Y vertex error",100,0.,0.02);
154	<	h["errPVz_"+suffix]  = subDir1.make<TH1D>(TString("errPVz_"+suffix),"Z vertex error",100,0.,0.02);
155	<	suffix.clear();
287	>	h_others->Init("others",suffix,ss.str());
288		}
157	–	}
158	–	}
289
290	+	// Book residual and pull histograms only when analyzeOntheFly is enabled
291	+	if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
292	+	suffix.clear();
293	+	} // End of Book histograms sensitive to data/mc
294	+
295	+
296	+	}
297
298		PVStudy::~PVStudy()
299		{
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 200 | 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;
210	<	cout << "number of vertices " << evt->vertices_size() << endl;
211	<	}
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;
365	>	int idx=0; int npv=0;
366		for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
367	<	vitr != evt->vertices_end(); ++vitr )
368	<	{ // loop for vertex ...
369	<	HepMC::FourVector pos = (*vitr)->position();
370	<	//HepLorentzVector pos = (*vitr)->position();
371	<	if (pos.t()>0) { continue;} // ?
372	<
373	<	bool hasMotherVertex=false;
374	<	//      if(verbose_) cout << "mothers" << endl;
375	<	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
376	<	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
377	<	HepMC::GenVertex * mv=(*mother)->production_vertex();
378	<	if (mv) {hasMotherVertex=true;}
379	<	//      if(verbose_) {
380	<	//        cout << "\t";
229	<	//        (*mother)->print();
230	<	//      }
231	<	}
232	<
233	<	if(hasMotherVertex) {continue;}
234	<
235	<	// could be a new vertex, check  all primaries found so far to avoid multiple entries
236	<	const double mm2cm=0.1;
237	<	simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
238	<	simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
239	<	for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
240	<	v0!=simpv.end(); v0++){
241	<	if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
242	<	vp=&(*v0);
243	<	break;
367	>	vitr != evt->vertices_end(); ++vitr ) { // loop for vertex ...
368	>	HepMC::FourVector pos = (*vitr)->position();
369	>	//HepLorentzVector pos = (*vitr)->position();
370	>
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	>	//        edm::LogInfo("SimPVs") << "Status = " << (*mother)->status() << endl;
375	>	HepMC::GenVertex * mv=(*mother)->production_vertex();
376	>	if( ((*mother)->status() == 3) && (!mv)) {
377	>	//        edm::LogInfo("SimPVs") << "npv= " << npv << endl;
378	>	if (npv == 0) {
379	>	h_gen->Fill1d("genPart_cT", pos.t()); // mm
380	>	h_gen->Fill1d("genPart_T", pos.t()/299.792458); // ns
381		}
382	+	npv++;
383	+	}
384	+	}
385	+	//       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
386	+
387	+	bool hasMotherVertex=false;
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	+	// 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
396		}
397	+	if(verbose_) {
398	+	cout << "\t";
399	+	(*mother)->print();
400	+	}
401	+	}
402
403	<	if(!vp){
404	<	// this is a new vertex
405	<	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
406	<	simpv.push_back(sv);
407	<	vp=&simpv.back();
408	<	}else{
409	<	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
410	<	}
411	<	vp->genVertex.push_back((*vitr)->barcode());
412	<	// collect final state descendants
413	<	for ( HepMC::GenVertex::particle_iterator
414	<	daughter  = (*vitr)->particles_begin(HepMC::descendants);
415	<	daughter != (*vitr)->particles_end(HepMC::descendants);
416	<	++daughter ) {
417	<	if (isFinalstateParticle(*daughter)){
418	<	if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
419	<	== vp->finalstateParticles.end()){
420	<	vp->finalstateParticles.push_back((*daughter)->barcode());
421	<	HepMC::FourVector m=(*daughter)->momentum();
422	<	// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
423	<	// but adding FourVectors seems not to be foreseen
424	<	vp->ptot.setPx(vp->ptot.px()+m.px());
425	<	vp->ptot.setPy(vp->ptot.py()+m.py());
426	<	vp->ptot.setPz(vp->ptot.pz()+m.pz());
427	<	vp->ptot.setE(vp->ptot.e()+m.e());
428	<	vp->ptsq+=(m.perp())*(m.perp());
429	<	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
430	<	vp->nGenTrk++;
431	<	}
403	>	if(hasMotherVertex) continue;
404	>
405	>	// could be a new vertex, check  all primaries found so far to avoid multiple entries
406	>	const double mm2cm=0.1;
407	>	simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
408	>	simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
409	>	for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
410	>	v0!=simpv.end(); v0++){
411	>	if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
412	>	vp=&(*v0);
413	>	break;
414	>	}
415	>	}
416	>
417	>	if(!vp){
418	>	// this is a new vertex
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	>	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
427	>	for ( HepMC::GenVertex::particle_iterator daughter  = (*vitr)->particles_begin(HepMC::descendants);
428	>	daughter != (*vitr)->particles_end(HepMC::descendants);
429	>	++daughter ) {
430	>	if (isFinalstateParticle(*daughter)){
431	>	if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
432	>	== vp->finalstateParticles.end()){
433	>	vp->finalstateParticles.push_back((*daughter)->barcode());
434	>	HepMC::FourVector m=(*daughter)->momentum();
435	>	// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
436	>	// but adding FourVectors seems not to be foreseen
437	>	vp->ptot.setPx(vp->ptot.px()+m.px());
438	>	vp->ptot.setPy(vp->ptot.py()+m.py());
439	>	vp->ptot.setPz(vp->ptot.pz()+m.pz());
440	>	vp->ptot.setE(vp->ptot.e()+m.e());
441	>	vp->ptsq+=(m.perp())*(m.perp());
442	>	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
443	>	vp->nGenTrk++;
444		}
445		}
446		}
447	<	idx++;
448	<	}
447	>	}//loop MC vertices daughters
448	>	}//loop MC vertices
449		}
450		return simpv;
451		}
#	Line 310 | 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 357 | Line 529 | void
529		PVStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
530		{
531		using namespace edm;
532	<	using reco::TrackCollection;
533	<
534	<	ftree_->SetBranchAddress("residual",&fres);
535	<	ftree_->SetBranchAddress("error",&ferror);
536	<	ftree_->SetBranchAddress("nTrk",&fntrk);
537	<
538	<	Handle<TrackCollection> tracks;
539	<	iEvent.getByLabel(tracksLabel_,tracks);
532	>	using namespace reco;
533	>
534	>	edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
535	>	//=======================================================
536	>	// Initialize Root-tuple variables if needed
537	>	//=======================================================
538	>	//if(saventuple_)
539	>	SetVarToZero();
540	>
541	>	//=======================================================
542	>	// Track accessors
543	>	//=======================================================
544	>	//trackCollection
545	>	static const reco::TrackCollection s_empty_trackColl;
546	>	const reco::TrackCollection *trackColl = &s_empty_trackColl;
547	>	edm::Handle<reco::TrackCollection> trackCollectionHandle;
548	>	iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle);
549	>	if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
550	>	trackColl = trackCollectionHandle.product();
551	>	} else {
552	>	edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
553	>	}
554	>	//splitTrackCollection1
555	>	static const reco::TrackCollection s_empty_splitTrackColl1;
556	>	const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
557	>	edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
558	>	iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle);
559	>	if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
560	>	splitTrackColl1 = splitTrackCollection1Handle.product();
561	>	} else {
562	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
563	>	}
564	>	//splitTrackCollection2
565	>	static const reco::TrackCollection s_empty_splitTrackColl2;
566	>	const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
567	>	edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
568	>	iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle);
569	>	if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
570	>	splitTrackColl2 = splitTrackCollection2Handle.product();
571	>	} else {
572	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
573	>	}
574	>
575	>	//=======================================================
576	>	// PVTX accessors
577	>	//=======================================================
578	>	//vertexCollection
579	>	static const reco::VertexCollection s_empty_vertexColl;
580	>	const reco::VertexCollection *vertexColl = &s_empty_vertexColl;
581	>	edm::Handle<reco::VertexCollection> vertexCollectionHandle;
582	>	iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle);
583	>	if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
584	>	vertexColl = vertexCollectionHandle.product();
585	>	} else {
586	>	edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
587	>	}
588	>	//splitVertexCollection1
589	>	static const reco::VertexCollection s_empty_splitVertexColl1;
590	>	const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
591	>	edm::Handle<reco::VertexCollection> splitVertexCollection1Handle;
592	>	iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle);
593	>	if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
594	>	splitVertexColl1 = splitVertexCollection1Handle.product();
595	>	} else {
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;
601	>	edm::Handle<reco::VertexCollection> splitVertexCollection2Handle;
602	>	iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle);
603	>	if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
604	>	splitVertexColl2 = splitVertexCollection2Handle.product();
605	>	} else {
606	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
607	>	}
608
609	<	Handle<reco::VertexCollection> recVtxs;
610	<	iEvent.getByLabel(vtxSample_, recVtxs);
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	+	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex and pixelvertices collections"<<endl;
633	+
634	+	//=======================================================
635	+	// MC simvtx accessor
636	+	//=======================================================
637		if (!realData_) {
638	<	Handle<SimVertexContainer> simVtxs;
638	>	edm::Handle<SimVertexContainer> simVtxs;
639		iEvent.getByLabel( simG4_, simVtxs);
640
641	<	Handle<SimTrackContainer> simTrks;
641	>	edm::Handle<SimTrackContainer> simTrks;
642		iEvent.getByLabel( simG4_, simTrks);
643		}
644
645	<	if(verbose_) printRecVtxs(recVtxs);
646	<
645	>	//=======================================================
646	>	// GET PDT
647	>	//=======================================================
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
654	<	if (!realData_) {
655	<	bool MC=false;
656	<	Handle<HepMCProduct> evtMC;
657	<	iEvent.getByLabel("generator",evtMC);
658	<	if (!evtMC.isValid()) {
659	<	MC=false;
660	<	if(verbose_){
661	<	cout << "no HepMCProduct found"<< endl;
662	<	}
663	<	} else {
664	<	if(verbose_){
665	<	cout << "generator HepMCProduct found"<< endl;
666	<	}
667	<	MC=true;
668	<	}
669	<
670	<	if(MC){
671	<	// make a list of primary vertices:
672	<	std::vector<simPrimaryVertex> simpv;
673	<	simpv=getSimPVs(evtMC,"");
674	<	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
675	<
676	<	int nsimtrk=0;
677	<	int nrectrk=0;
678	<	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
679	<	vsim!=simpv.end(); vsim++){
680	<	nsimtrk+=vsim->nGenTrk;
681	<	// look for a matching reconstructed vertex
682	<	vsim->recVtx=NULL;
683	<	for(reco::VertexCollection::const_iterator vrec=recVtxs->begin();
684	<	vrec!=recVtxs->end(); ++vrec){
685	<	nrectrk=vrec->tracksSize();
686	<	if(verbose_){
687	<	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
688	<	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
689	<	}
690	<	if ( matchVertex(*vsim,*vrec) ){
691	<	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
692	<	if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
693	<	|| (!vsim->recVtx) )
694	<	{
695	<	vsim->recVtx=&(*vrec);
696	<	}
697	<	}
698	<	}
699	<	if (vsim->recVtx){
700	<	if(verbose_) {
701	<	cout <<"primary matched " << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
702	<	}
703	<	fres[0] = vsim->recVtx->x()-vsim->x;
704	<	fres[1] = vsim->recVtx->y()-vsim->y;
705	<	fres[2] = vsim->recVtx->z()-vsim->z;
440	<	ferror[0] = vsim->recVtx->xError();
441	<	ferror[1] = vsim->recVtx->yError();
442	<	ferror[2] = vsim->recVtx->zError();
443	<	fntrk = nrectrk;
444	<
445	<	h["deltax"]->Fill( fres[0] );
446	<	h["deltay"]->Fill( fres[1] );
447	<	h["deltaz"]->Fill( fres[2] );
448	<	h["pullx"]->Fill( fres[0]/ferror[0] );
449	<	h["pully"]->Fill( fres[1]/ferror[1] );
450	<	h["pullz"]->Fill( fres[2]/ferror[2] );
451	<	h["errPVx"]->Fill( ferror[0] );
452	<	h["errPVy"]->Fill( ferror[1] );
453	<	h["errPVz"]->Fill( ferror[2] );
454	<	pvinfo.push_back(PVStudy::PVInfo(res(fres[0], fres[1], fres[2]),
455	<	error(ferror[0], ferror[1], ferror[2]),
456	<	nrectrk));
457	<	// Fill histo according to its track multiplicity
458	<	fillHisto(res(fres[0], fres[1], fres[2]),
459	<	error(ferror[0], ferror[1], ferror[2]),
460	<	nrectrk);
461	<
462	<	ftree_->Fill();
463	<	}
464	<	else{  // no rec vertex found for this simvertex
465	<	if(verbose_) {
466	<	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
467	<	}
468	<	}
654	>	//=======================================================
655	>	// Apply event cleaning for firstcoll data and MC
656	>	//=======================================================
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, 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_pvtrk->Fill1d("trkPtPV", 0.);
706		}
707	+	else
708	+	h_pixvtx->Fill1d("trkdzErr_pxlpvtx", (**t).dzError());
709		}
710	<	}// With MC
711	<
712	<	// Loop RecVtxs and find matched SimVtxs
713	<	for(reco::VertexCollection::const_iterator v=recVtxs->begin();
714	<	v!=recVtxs->end(); ++v){
715	<	try {
716	<	for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
717	<	t!=v->tracks_end(); t++) {
718	<	// illegal charge
719	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
720	<	//      h["trkPtPV"]->Fill(0.);
721	<	}
722	<	else {
484	<	h["trkPtPV"]->Fill((**t).pt());
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_pvtrk->Fill1d("trkPtPV", 0.);
716	>	}
717	>	else {
718	>	if(pixelVertexColl->size()>0) {
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		}
488	–	catch (...) {
489	–	// exception thrown when trying to use linked track
490	–	//          h["trkPtPV"]->Fill(0.);
491	–	}
492	–
493	–	h["nTrkPV"]->Fill(v->tracksSize());
494	–
726		}
727	+
728	+	//=======================================================
729	+	// Fill number of reconstructed vertices
730	+	//=======================================================
731	+
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);
738	+	printRecVtxs(splitVertexCollection1Handle);
739	+	printRecVtxs(splitVertexCollection2Handle);
740	+	}
741	+	nrecPV_ = int (vertexColl->size());
742	+	nrecPV1_spl_ = int (splitVertexColl1->size());
743	+	nrecPV2_spl_ = int (splitVertexColl2->size());
744	+
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
497	–	h["nTrk"]->Fill(tracks->size());
750
751	<	for(TrackCollection::const_iterator itTrack = tracks->begin();
752	<	itTrack != tracks->end();
753	<	++itTrack) {
754	<	int charge = 0;
755	<	charge = itTrack->charge();
756	<	h["trkPt"]->Fill(itTrack->pt());
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, const Point & bs) {
756	>	if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
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, 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_pixvtx->Fill1d("nrecPV_minus_nrecPxlPV", double (nrecPV_ - nrecPV_pxlpvtx_));
771	>	// difference in reconstructed position of the leading pvtx
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	+	//==========================================================
781	+	// Fill secondary/primary min separations for multi-vertices
782	+	//==========================================================
783	+
784	+	if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {
785	+
786	+	double min_xsep = 9999.0;
787	+	double min_ysep = 9999.0;
788	+	double min_zsep = 9999.0;
789	+	double min_xsep_sign = 9999.0;
790	+	double min_ysep_sign = 9999.0;
791	+	double min_zsep_sign = 9999.0;
792	+	double min_ntrksep = 9999.0;
793	+	double min_sumpt2sep = 9999999.0;
794	+
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 ;
799	+	v!=vertexColl->end(); ++v) {
800	+	if(v->isValid() && ! v->isFake()) {
801	+	// xsep
802	+	if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
803	+	min_xsep = fabs(v->x()- vertexColl->begin()->x());
804	+	// ysep
805	+	if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
806	+	min_ysep = fabs(v->y()- vertexColl->begin()->y());
807	+	// zsep
808	+	if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
809	+	min_zsep = fabs(v->z()- vertexColl->begin()->z());
810	+	// xsep_sign
811	+	double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
812	+	if(xsep_sign < min_xsep_sign)
813	+	min_xsep_sign = xsep_sign;
814	+	// ysep_sign
815	+	double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
816	+	if(ysep_sign < min_ysep_sign)
817	+	min_ysep_sign = ysep_sign;
818	+	// zsep_sign
819	+	double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
820	+	if(zsep_sign < min_zsep_sign)
821	+	min_zsep_sign = zsep_sign;
822	+	// ntrksep
823	+	if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
824	+	min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
825	+	// sumpt2sep
826	+	if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
827	+	min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
828	+	}
829	+	}
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;
841	+	min_sumpt2sep_ = min_sumpt2sep;
842	+
843	+
844	+	} // End of  if(vertexColl->size()>1) {
845	+
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)
853	+	// 2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
854	+	// Assume the first match is the primary vertex,
855	+	// since vertexColl are sorted in decreasing order of sum(pT)
856	+	//=======================================================
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());
862	+	reco::VertexCollection s_empty_matchedVertexColl2;
863	+	reco::VertexCollection *matchedVertexColl2 = &s_empty_matchedVertexColl2;
864	+	matchedVertexColl2->reserve(splitVertexColl2->size());
865	+
866	+	bool stopmatching = false;
867	+	for (size_t ivtx = 0; ivtx<splitVertexCollection1Handle->size() && !stopmatching; ++ivtx) {
868	+	reco::VertexRef recvtx1(splitVertexCollection1Handle, ivtx);
869	+	if( !(recvtx1->isValid()) || recvtx1->isFake()) break;
870	+	for (size_t jvtx = ivtx; jvtx < splitVertexCollection2Handle->size(); ++jvtx) {
871	+	//------------------------------------------------------------------------
872	+	//== If only considering matching the first vertex of the splitVertexColl
873	+	//------------------------------------------------------------------------
874	+	if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
875	+	reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
876	+	if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
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	+	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_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
890	+	h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
891	+	if(fabs(ntrkreldiff)<ntrkdiffcut_) {
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	+	edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
901	+	}
902	+	else {
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	+	edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
910	+
911	+	//=======================================================
912	+	//  Analyze matchedVertexColls
913	+	//=======================================================
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, 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) {
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);
933	+	ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
934	+	ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
935	+	ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
936	+	fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
937	+
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_));
950	+	// Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
951	+	fillHisto(res(fres_[0], fres_[1], fres_[2]),
952	+	error(ferror_[0], ferror_[1], ferror_[2]),
953	+	fntrk_,0);
954	+	// Fill the ntuple variables
955	+	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
956	+	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
957	+	sumptsq1_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v1);
958	+	sumptsq2_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v2);
959	+	nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
960	+	deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
961	+	deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
962	+	deltaz_twovtx_[nrecPV_twovtx_] = fres_[2];
963	+	errx_twovtx_[nrecPV_twovtx_] = ferror_[0];
964	+	erry_twovtx_[nrecPV_twovtx_] = ferror_[1];
965	+	errz_twovtx_[nrecPV_twovtx_] = ferror_[2];
966	+	pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
967	+	pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
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	+	edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
1015	+
1016	+	edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
1017	+
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 516 | Line 1036 | PVStudy::beginJob()
1036		// ------------ method called once each job just after ending the event loop  ------------
1037		void
1038		PVStudy::endJob() {
1039	<	double sqt2 = sqrt(2.);
1040	<	if (!realData_) {
1041	<	if (verbose_) cout << "Analyzing PV info" << endl;
1042	<	for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1043	<	if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
1044	<	PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk);
1045	<	if (analyze_) {
1046	<	//      if ( ipv.res_.x() > 0 ) h2["resx_nTrk"]->Fill(ntrk,ipv.res_.x()/sqt2, 1.);
1047	<	//      if ( ipv.res_.y() > 0 ) h2["resy_nTrk"]->Fill(ntrk,ipv.res_.y()/sqt2, 1.);
1048	<	//      if ( ipv.res_.z() > 0 ) h2["resz_nTrk"]->Fill(ntrk,ipv.res_.z()/sqt2, 1.);
1049	<	if ( ipv.res_.x() > 0 ) h2["resx_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
1050	<	if ( ipv.res_.y() > 0 ) h2["resy_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
1051	<	if ( ipv.res_.z() > 0 ) h2["resz_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
1052	<	if ( ipv.pull_.x() > 0 ) h2["pullx_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
1053	<	if ( ipv.pull_.y() > 0 ) h2["pully_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
1054	<	if ( ipv.pull_.z() > 0 ) h2["pullz_nTrk"]->Fill(ntrk,ipv.pull_.z(), 1.);
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	>	edm::LogInfo("Analysis")<<"[endJob] datamode = "<< *it<<endl;
1044	>	switch(*it) {
1045	>	case 0: suffix = "";
1046	>	break;
1047	>	case 1: suffix = "_spl1_mct";
1048	>	break;
1049	>	case 2: suffix = "_spl2_mct";
1050	>	break;
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 ) 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		}
1067	+	suffix.clear();
1068	+	}
1069	+	if(saventuple_) {
1070	+	file_->cd();
1071	+	ftree_->Write();
1072	+	pvtxtree_->Write();
1073	+	file_->Close();
1074		}
1075	<	file_->cd();
1076	<	ftree_->Write();
540	<	file_->Close();
1075	>
1076	>
1077		}
1078
1079	+	// Match a recovertex with a simvertex
1080	+	// ? Should the cut be parameter dependent?
1081	+	// cut can be within n sigma of the vertex.zerror()
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
1091	+	bool PVStudy::matchVertex( const reco::VertexRef & vrec1,
1092	+	const reco::VertexRef & vrec2,
1093	+	int zsigncut)
1094	+	{
1095	+	double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
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 617 | Line 1169 | void PVStudy::printSimTrks(const Handle<
1169		}
1170		}
1171
1172	<	void PVStudy::fillHisto(res r, error er, int ntk){
1172	>	void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1173		stringstream ss;
1174		ss << ntk;
1175	<	string suffix = ss.str();
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["resx_"+suffix]->Fill(r.x());
1184	<	h["resy_"+suffix]->Fill(r.y());
1185	<	h["resz_"+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	<	PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk) {
1194	>
1195	>	PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1196		map<int,double> data;
1197		data.clear();
1198		double fpar[2] = {-999,-999};
1199		double cm2um = 10000;
1200		stringstream ss;
1201		ss << ntk;
1202	<	string suffix = ss.str();
1203	<
1202	>	string suffix = ss.str();
1203	>	if(datamode == 0 )  suffix = "_" + suffix;
1204	>	if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1205	>	if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1206	>	if(datamode == 3 )  suffix = "_mct_" + suffix;
1207	>
1208		if(analyze_) {
1209		for ( int i = 0; i < 6; ++i) {
1210		switch (i) {
1211		case 0:
1212	<	fit(h["resx_"+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["resy_"+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["resz_"+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 686 | 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;
704	<	}
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 713 | 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	<	//define this as a plug-in
1306	<	DEFINE_FWK_MODULE(PVStudy);
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_pvtrk->Fill1d("nTrk"+suffix, trackColl->size());
1313	>	for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1314	>	itTrack != trackColl->end();
1315	>	++itTrack) {
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, const Point & bs) {
1326	>	string suffix;
1327	>	suffix = ""; // for unsplit tracks
1328	>	if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1329	>	if(datatype == 2) suffix = "2_spl"; // for splittracks 2
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();
1337	>	sumptsq_[ivtx] = sumPtSquared(*v);
1338	>	isValid_[ivtx] = v->isValid();
1339	>	isFake_[ivtx] = v->isFake();
1340	>	recx_[ivtx] = v->x();
1341	>	recy_[ivtx] = v->y();
1342	>	recz_[ivtx] = v->z();
1343	>	recx_err_[ivtx] = v->xError();
1344	>	recy_err_[ivtx] = v->yError();
1345	>	recz_err_[ivtx] = v->zError();
1346	>	}
1347	>	if(datatype == 1) {
1348	>	nTrkPV1_spl_[ivtx] = v->tracksSize();
1349	>	sumptsq1_spl_[ivtx] = sumPtSquared(*v);
1350	>	isValid1_spl_[ivtx] = v->isValid();
1351	>	isFake1_spl_[ivtx] = v->isFake();
1352	>	recx1_spl_[ivtx] = v->x();
1353	>	recy1_spl_[ivtx] = v->y();
1354	>	recz1_spl_[ivtx] = v->z();
1355	>	recx1_err_spl_[ivtx] = v->xError();
1356	>	recy1_err_spl_[ivtx] = v->yError();
1357	>	recz1_err_spl_[ivtx] = v->zError();
1358	>
1359	>	}
1360	>	if(datatype == 2) {
1361	>	nTrkPV2_spl_[ivtx] = v->tracksSize();
1362	>	sumptsq2_spl_[ivtx] = sumPtSquared(*v);
1363	>	isValid2_spl_[ivtx] = v->isValid();
1364	>	isFake2_spl_[ivtx] = v->isFake();
1365	>	recx2_spl_[ivtx] = v->x();
1366	>	recy2_spl_[ivtx] = v->y();
1367	>	recz2_spl_[ivtx] = v->z();
1368	>	recx2_err_spl_[ivtx] = v->xError();
1369	>	recy2_err_spl_[ivtx] = v->yError();
1370	>	recz2_err_spl_[ivtx] = v->zError();
1371	>	}
1372	>	if(datatype == 4) { // for pixelVertices
1373	>	nTrkPV_pxlpvtx_[ivtx] = v->tracksSize();
1374	>	sumptsq_pxlpvtx_[ivtx] = sumPtSquared(*v);
1375	>	isValid_pxlpvtx_[ivtx] = v->isValid();
1376	>	isFake_pxlpvtx_[ivtx] = v->isFake();
1377	>	recx_pxlpvtx_[ivtx] = v->x();
1378	>	recy_pxlpvtx_[ivtx] = v->y();
1379	>	recz_pxlpvtx_[ivtx] = v->z();
1380	>	recx_err_pxlpvtx_[ivtx] = v->xError();
1381	>	recy_err_pxlpvtx_[ivtx] = v->yError();
1382	>	recz_err_pxlpvtx_[ivtx] = v->zError();
1383	>	}
1384	>	}
1385	>	}
1386	>	// For histogram only fill the leading pvtx parameters
1387	>	if(fillHisto) {
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_pvtrk->Fill1d("trkPtPV", 0.);
1398	>	}
1399	>	else {
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_pvtrk->Fill1d("trkPtPV", 0.);
1415	>	}
1416	>	h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
1417	>	}
1418	>
1419	>
1420	>	}
1421	>
1422	>	void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
1423	>	{
1424	>	std::string suffix;
1425	>	// Set the vertexColl and histogram suffix according to datamode
1426	>	if (datatype == 1) {
1427	>	suffix = "_spl1_mct";
1428	>	nrecPV_spl1_mct_ = 0;
1429	>	}
1430	>	if (datatype == 2) {
1431	>	suffix = "_spl2_mct";
1432	>	nrecPV_spl2_mct_ = 0;
1433	>	}
1434	>	if (datatype == 3) {
1435	>	suffix = "_mct";
1436	>	nrecPV_mct_ = 0;
1437	>	}
1438	>
1439	>	//========================================================
1440	>	//  look for a matching reconstructed vertex in vertexColl
1441	>	//========================================================
1442	>
1443	>	for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1444	>	vrec!=vtxColl->end(); ++vrec){
1445	>	int nrectrk = vrec->tracksSize();
1446	>	vsim->recVtx=NULL;
1447	>
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
1455	>	//if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
1456	>	//|| (!vsim->recVtx) )
1457	>	//vsim->recVtx=&(*vrec);
1458	>	//}
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];
1463	>
1464	>	fres_mct[0] = vsim->recVtx->x()-vsim->x;
1465	>	fres_mct[1] = vsim->recVtx->y()-vsim->y;
1466	>	fres_mct[2] = vsim->recVtx->z()-vsim->z;
1467	>	ferror_mct[0] = vsim->recVtx->xError();
1468	>	ferror_mct[1] = vsim->recVtx->yError();
1469	>	ferror_mct[2] = vsim->recVtx->zError();
1470	>
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));
1483	>	// Fill histo according to its track multiplicity
1484	>	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1485	>	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1486	>	nrectrk, datatype);
1487	>
1488	>	if(saventuple_) {
1489	>	//Fill the values for variables in ftree_
1490	>	if(datatype == 1) {
1491	>	nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1492	>	deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1493	>	deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1494	>	deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1495	>	pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
1496	>	pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1497	>	pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1498	>	nrecPV_spl1_mct_++;
1499	>	}
1500	>
1501	>	if(datatype == 2) {
1502	>	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1503	>	deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1504	>	deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1505	>	deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1506	>	pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
1507	>	pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1508	>	pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1509	>	nrecPV_spl2_mct_++;
1510	>	}
1511	>	if(datatype == 3) {
1512	>	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1513	>	deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1514	>	deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1515	>	deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1516	>	pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
1517	>	pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1518	>	pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1519	>	nrecPV_mct_++;
1520	>	}
1521	>	} // End of  if(saventuple_) {
1522	>	} //  if ( matchVertex(*vsim,*vrec) ) {
1523	>	else {  // no rec vertex found for this simvertex
1524	>	edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1525	>	}
1526	>	}
1527	>	}
1528	>
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;
1536	>	}
1537	>
1538	>	//Yanyan's variables
1539	>	// Number of reconstructed vertices
1540	>	nrecPV_ = 0;
1541	>	nrecPV1_spl_ = 0;
1542	>	nrecPV2_spl_ = 0;
1543	>	nrecPV_twovtx_ = 0;
1544	>	nsimPV_ = 0;
1545	>
1546	>	nrecPV_mct_ = 0;
1547	>	nrecPV_spl1_mct_ = 0;
1548	>	nrecPV_spl2_mct_ = 0;
1549	>
1550	>	// Mininum separation between the secondary pvtxes and leading pvtx
1551	>	min_zsep_ = 9999.0;
1552	>	min_ntrksep_ = 9999.0;
1553	>	min_sumpt2sep_ = -9999.0;
1554	>
1555	>
1556	>	for (int i = 0; i < nMaxPVs_; i++) {
1557	>	// recoVertices with all tracks
1558	>	nTrkPV_[i] = 0; // Number of tracks in the pvtx
1559	>	sumptsq_[i] = 0;
1560	>	isValid_[i] = -1;
1561	>	isFake_[i] = -1;
1562	>	recx_[i] = 0;
1563	>	recy_[i] = 0;
1564	>	recz_[i] = 0;
1565	>	recx_err_[i] = 0;
1566	>	recy_err_[i] = 0;
1567	>	recz_err_[i] = 0;
1568	>
1569	>	// recoVertices with splitTrack1
1570	>	nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx
1571	>	sumptsq1_spl_[i] = 0;
1572	>	isValid1_spl_[i] = -1;
1573	>	isFake1_spl_[i] = -1;
1574	>	recx1_spl_[i] = 0;
1575	>	recy1_spl_[i] = 0;
1576	>	recz1_spl_[i] = 0;
1577	>	recx1_err_spl_[i] = 0;
1578	>	recy1_err_spl_[i] = 0;
1579	>	recz1_err_spl_[i] = 0;
1580	>
1581	>	// recoVertices with splitTrack2
1582	>	nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx
1583	>	sumptsq2_spl_[i] = 0;
1584	>	isValid2_spl_[i] = -1;
1585	>	isFake2_spl_[i] = -1;
1586	>	recx2_spl_[i] = 0;
1587	>	recy2_spl_[i] = 0;
1588	>	recz2_spl_[i] = 0;
1589	>	recx2_err_spl_[i] = 0;
1590	>	recy2_err_spl_[i] = 0;
1591	>	recz2_err_spl_[i] = 0;
1592	>
1593	>	//pixelVertices
1594	>	nTrkPV_pxlpvtx_[i] = 0; // Number of tracks in the pvtx
1595	>	sumptsq_pxlpvtx_[i] = 0;
1596	>	isValid_pxlpvtx_[i] = -1;
1597	>	isFake_pxlpvtx_[i] = -1;
1598	>	recx_pxlpvtx_[i] = 0;
1599	>	recy_pxlpvtx_[i] = 0;
1600	>	recz_pxlpvtx_[i] = 0;
1601	>	recx_err_pxlpvtx_[i] = 0;
1602	>	recy_err_pxlpvtx_[i] = 0;
1603	>	recz_err_pxlpvtx_[i] = 0;
1604	>
1605	>	// matched two-vertices
1606	>	nTrkPV1_spl_twovtx_[i] = 0;
1607	>	nTrkPV2_spl_twovtx_[i] = 0;
1608	>	nTrkPV_twovtx_[i] = 0;
1609	>	sumptsq1_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1610	>	sumptsq2_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1611	>	deltax_twovtx_[i] = 0;
1612	>	deltay_twovtx_[i] = 0;
1613	>	deltaz_twovtx_[i] = 0;
1614	>	errx_twovtx_[i] = 0;
1615	>	erry_twovtx_[i] = 0;
1616	>	errz_twovtx_[i] = 0;
1617	>	pullx_twovtx_[i] = 0;
1618	>	pully_twovtx_[i] = 0;
1619	>	pullz_twovtx_[i] = 0;
1620	>
1621	>	//simpv
1622	>	nsimTrkPV_[i] = 0;
1623	>	simx_[i] = 0;
1624	>	simy_[i] = 0;
1625	>	simz_[i] = 0;
1626	>	simptsq_[i] = 0;
1627	>
1628	>	// residual and pulls with mc truth required
1629	>	deltax_mct_[i] = 0;
1630	>	deltay_mct_[i] = 0;
1631	>	deltaz_mct_[i] = 0;
1632	>	pullx_mct_[i] = 0;
1633	>	pully_mct_[i] = 0;
1634	>	pullz_mct_[i] = 0;
1635	>
1636	>	deltax_spl1_mct_[i] = 0;
1637	>	deltay_spl1_mct_[i] = 0;
1638	>	deltaz_spl1_mct_[i] = 0;
1639	>	pullx_spl1_mct_[i] = 0;
1640	>	pully_spl1_mct_[i] = 0;
1641	>	pullz_spl1_mct_[i] = 0;
1642	>
1643	>	deltax_spl2_mct_[i] = 0;
1644	>	deltay_spl2_mct_[i] = 0;
1645	>	deltaz_spl2_mct_[i] = 0;
1646	>	pullx_spl2_mct_[i] = 0;
1647	>	pully_spl2_mct_[i] = 0;
1648	>	pullz_spl2_mct_[i] = 0;
1649	>	}
1650	>	}
1651	>
1652	>	double PVStudy::sumPtSquared(const reco::Vertex & v)  {
1653	>	double sum = 0.;
1654	>	double pT;
1655	>	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1656	>	pT = (**it).pt();
1657	>	sum += pT*pT;
1658	>	}
1659	>	return sum;
1660	>	}
1661	>
1662	>
1663	>

Diff Legend

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