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.3 by jengbou, Tue Sep 8 07:15:38 2009 UTC vs.
Revision 1.14 by jengbou, Fri Dec 4 20:48:04 2009 UTC

#	Line 12 | Line 12 | Implementation:
12		*/
13		//
14		// Original Author:  "Geng-yuan Jeng/UC Riverside"
15	+	//                   "Yanyan Gao/Fermilab ygao@fnal.gov"
16		//         Created:  Thu Aug 20 11:55:40 CDT 2009
17		// $Id$
18		//
#	Line 24 | Line 25 | Implementation:
25		#include <vector>
26		#include <iostream>
27		#include <sstream>
28	+	#include <algorithm>
29
30		// user include files
31		#include "FWCore/Framework/interface/Frameworkfwd.h"
#	Line 50 | Line 52 | Implementation:
52		#include <SimDataFormats/Track/interface/SimTrackContainer.h>
53
54		//root
55	<	#include "TROOT.h"
56	<	#include "TF1.h"
57	<	#include "TString.h"
58	<	#include "TStyle.h"
59	<
55	>	#include <TROOT.h>
56	>	#include <TF1.h>
57	>	#include <TString.h>
58	>	#include <TStyle.h>
59	>	#include <TPaveStats.h>
60	>	#include <TPad.h>
61
62		using namespace std;
63
64	<	PVStudy::PVStudy(const edm::ParameterSet& iConfig):
62	<	nTrkMin_(10),nTrkMax_(100)
64	>	PVStudy::PVStudy(const edm::ParameterSet& iConfig)
65		{
66	<	simG4_           = iConfig.getParameter<edm::InputTag>( "simG4" );
67	<	tracksLabel_     = iConfig.getUntrackedParameter<edm::InputTag>("tracks");
68	<	vtxSample_       = iConfig.getUntrackedParameter<edm::InputTag>("vtxSample");
69	<	verbose_         = iConfig.getUntrackedParameter<bool>("verbose", false);
70	<	realData_        = iConfig.getUntrackedParameter<bool>("realData",false);
71	<	analyze_         = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
72	<	outputfilename_  = iConfig.getUntrackedParameter<string>("OutputFileName");
66	>	//=======================================================================
67	>	// Get configuration for TrackTupleMaker
68	>	//=======================================================================
69	>	simG4_                     = iConfig.getParameter<edm::InputTag>( "simG4" );
70	>	trackCollectionTag_        = iConfig.getUntrackedParameter<edm::InputTag>("trackCollection");
71	>	splitTrackCollection1Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection1");
72	>	splitTrackCollection2Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
73	>	vertexCollectionTag_       = iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection");
74	>	splitVertexCollection1Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
75	>	splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
76	>	pixelVertexCollectionTag_  = iConfig.getParameter<edm::InputTag>("pixelVertexCollectionTag");
77	>	verbose_                   = iConfig.getUntrackedParameter<bool>("verbose",false);
78	>	realData_                  = iConfig.getUntrackedParameter<bool>("realData",false);
79	>	analyze_                   = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
80	>	saventuple_                = iConfig.getUntrackedParameter<bool>("saventuple",false);
81	>	outputfilename_            = iConfig.getUntrackedParameter<string>("OutputFileName");
82	>	ntrkdiffcut_               = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
83	>	zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
84	>	nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
85	>	nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
86	>	histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
87
88		//now do what ever initialization is needed
89		pvinfo.clear();
90
91	<	file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
92	<	ftree_ = new TTree("mytree","mytree");
93	<	ftree_->AutoSave();
94	<	ftree_->Branch("residual",&fres,"fres[3]/D");
95	<	ftree_->Branch("error",&ferror,"ferror[3]/D");
96	<	ftree_->Branch("nTrk",&fntrk,"fntrk/I");
97	<
98	<	edm::Service<TFileService> fs;
99	<	//   TFileDirectory subDir = fs->mkdir( "Summary" );
100	<	//   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
101	<
102	<	h_nTrk           = fs->make<TH1F>("nTrk","Num of total tracks",300,0,300);
103	<	h_nTrkPV         = fs->make<TH1F>("nTrkPV","Num of Tracks from PV",300,0,300);
104	<	h_trkPt          = fs->make<TH1D>("trkPt","Pt of all tracks",100,0,100);
105	<	h_trkPtPV        = fs->make<TH1D>("trkPtPV","Pt dist of tracks from PV",100,0,100);
91	>
92	>	// Specify the data mode vector
93	>	if(realData_) datamode.push_back(0);
94	>	else {
95	>	datamode.push_back(0);
96	>	datamode.push_back(1);
97	>	datamode.push_back(2);
98	>	datamode.push_back(3);
99	>	}
100	>
101	>	// Create ntuple files if needed
102	>	if(saventuple_) {
103	>	file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
104	>	ftree_ = new TTree("mytree","mytree");
105	>	ftree_->AutoSave();
106	>	ftree_->Branch("residual",&fres_,"fres_[3]/D");
107	>	ftree_->Branch("error",&ferror_,"ferror_[3]/D");
108	>	ftree_->Branch("nTrkPV",&fntrk_,"fntrk_/I");
109	>
110	>	// pvtxtree_ analyzing the pvtxs ootb
111	>	pvtxtree_ = new TTree("pvtxtree","pvtxtree");
112	>	//pvtx using all tracks
113	>	pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
114	>	pvtxtree_->Branch("nTrkPV",&nTrkPV_,"nTrkPV[nrecPV]/I");
115	>	pvtxtree_->Branch("sumptsq",&sumptsq_,"sumptsq[nrecPV]/D");
116	>	pvtxtree_->Branch("isValid",&isValid_,"isValid/I");
117	>	pvtxtree_->Branch("isFake",&isFake_,"isFake/I");
118	>	pvtxtree_->Branch("recx",&recx_,"recx[nrecPV]/D");
119	>	pvtxtree_->Branch("recy",&recy_,"recy[nrecPV]/D");
120	>	pvtxtree_->Branch("recz",&recz_,"recz[nrecPV]/D");
121	>	pvtxtree_->Branch("recx_err",&recx_err_,"recx_err[nrecPV]/D");
122	>	pvtxtree_->Branch("recy_err",&recy_err_,"recy_err[nrecPV]/D");
123	>	pvtxtree_->Branch("recz_err",&recz_err_,"recz_err[nrecPV]/D");
124	>
125	>	pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
126	>	pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
127	>	pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
128	>
129	>	//pvtx using splittracks1
130	>	pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
131	>	pvtxtree_->Branch("nTrkPV1_spl",&nTrkPV1_spl_,"nTrkPV1_spl[nrecPV1_spl]/I");
132	>	pvtxtree_->Branch("sumptsq1_spl",&sumptsq1_spl_,"sumptsq1_spl[nrecPV1_spl]/D");
133	>	pvtxtree_->Branch("isValid1_spl",&isValid1_spl_,"isValid1_spl/I");
134	>	pvtxtree_->Branch("isFake1_spl",&isFake1_spl_,"isFake1_spl/I");
135	>	pvtxtree_->Branch("recx1_spl",&recx1_spl_,"recx1_spl[nrecPV1_spl]/D");
136	>	pvtxtree_->Branch("recy1_spl",&recy1_spl_,"recy1_spl[nrecPV1_spl]/D");
137	>	pvtxtree_->Branch("recz1_spl",&recz1_spl_,"recz1_spl[nrecPV1_spl]/D");
138	>	pvtxtree_->Branch("recx1_err_spl",&recx1_err_spl_,"recx1_err_spl[nrecPV1_spl]/D");
139	>	pvtxtree_->Branch("recy1_err_spl",&recy1_err_spl_,"recy1_err_spl[nrecPV1_spl]/D");
140	>	pvtxtree_->Branch("recz1_err_spl",&recz1_err_spl_,"recz1_err_spl[nrecPV1_spl]/D");
141	>
142	>	//pvtx using splittracks2
143	>	pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
144	>	pvtxtree_->Branch("nTrkPV2_spl",&nTrkPV2_spl_,"nTrkPV2_spl[nrecPV2_spl]/I");
145	>	pvtxtree_->Branch("sumptsq2_spl",&sumptsq2_spl_,"sumptsq2_spl[nrecPV2_spl]/D");
146	>	pvtxtree_->Branch("isValid2_spl",&isValid2_spl_,"isValid2_spl/I");
147	>	pvtxtree_->Branch("isFake2_spl",&isFake2_spl_,"isFake2_spl/I");
148	>	pvtxtree_->Branch("recx2_spl",&recx2_spl_,"recx2_spl[nrecPV2_spl]/D");
149	>	pvtxtree_->Branch("recy2_spl",&recy2_spl_,"recy2_spl[nrecPV2_spl]/D");
150	>	pvtxtree_->Branch("recz2_spl",&recz2_spl_,"recz2_spl[nrecPV2_spl]/D");
151	>	pvtxtree_->Branch("recx2_err_spl",&recx2_err_spl_,"recx2_err_spl[nrecPV2_spl]/D");
152	>	pvtxtree_->Branch("recy2_err_spl",&recy2_err_spl_,"recy2_err_spl[nrecPV2_spl]/D");
153	>	pvtxtree_->Branch("recz2_err_spl",&recz2_err_spl_,"recz2_err_spl[nrecPV2_spl]/D");
154	>
155	>	//pixeVertices
156	>	pvtxtree_->Branch("nrecPV_pxlpvtx",&nrecPV_pxlpvtx_,"nrecPV_pxlpvtx/I");
157	>	pvtxtree_->Branch("nTrkPV_pxlpvtx",&nTrkPV_pxlpvtx_,"nTrkPV_pxlpvtx[nrecPV_pxlpvtx]/I");
158	>	pvtxtree_->Branch("sumptsq_pxlpvtx",&sumptsq_pxlpvtx_,"sumptsq_pxlpvtx[nrecPV_pxlpvtx]/D");
159	>	pvtxtree_->Branch("isValid_pxlpvtx",&isValid_pxlpvtx_,"isValid_pxlpvtx/I");
160	>	pvtxtree_->Branch("isFake_pxlpvtx",&isFake_pxlpvtx_,"isFake_pxlpvtx/I");
161	>	pvtxtree_->Branch("recx_pxlpvtx",&recx_pxlpvtx_,"recx_pxlpvtx[nrecPV_pxlpvtx]/D");
162	>	pvtxtree_->Branch("recy_pxlpvtx",&recy_pxlpvtx_,"recy_pxlpvtx[nrecPV_pxlpvtx]/D");
163	>	pvtxtree_->Branch("recz_pxlpvtx",&recz_pxlpvtx_,"recz_pxlpvtx[nrecPV_pxlpvtx]/D");
164	>	pvtxtree_->Branch("recx_err_pxlpvtx",&recx_err_pxlpvtx_,"recx_err_pxlpvtx[nrecPV_pxlpvtx]/D");
165	>	pvtxtree_->Branch("recy_err_pxlpvtx",&recy_err_pxlpvtx_,"recy_err_pxlpvtx[nrecPV_pxlpvtx]/D");
166	>	pvtxtree_->Branch("recz_err_pxlpvtx",&recz_err_pxlpvtx_,"recz_err_pxlpvtx[nrecPV_pxlpvtx]/D");
167	>
168	>	//Fill the variables in the twovtx pair (recvtx1, recvtx2)
169	>	pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
170	>	pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
171	>	pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
172	>	pvtxtree_->Branch("sumptsq1_spl_twovtx",&sumptsq1_spl_twovtx_,"sumptsq1_spl_twovtx[nrecPV_twovtx]/D");
173	>	pvtxtree_->Branch("sumptsq2_spl_twovtx",&sumptsq2_spl_twovtx_,"sumptsq2_spl_twovtx[nrecPV_twovtx]/D");
174	>	pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
175	>	pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
176	>	pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
177	>	pvtxtree_->Branch("deltaz_twovtx",&deltaz_twovtx_,"deltaz_twovtx[nrecPV_twovtx]/D");
178	>	pvtxtree_->Branch("errx_twovtx",&errx_twovtx_,"errx_twovtx[nrecPV_twovtx]/D");
179	>	pvtxtree_->Branch("erry_twovtx",&erry_twovtx_,"erry_twovtx[nrecPV_twovtx]/D");
180	>	pvtxtree_->Branch("errz_twovtx",&errz_twovtx_,"errz_twovtx[nrecPV_twovtx]/D");
181	>	pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
182	>	pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
183	>	pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
184	>
185	>	// MC variables
186	>	if(!realData_) {
187	>	pvtxtree_->Branch("nsimPV",&nsimPV_,"nsimPV/I");
188	>	pvtxtree_->Branch("nsimTrkPV",&nsimTrkPV_,"nsimTrkPV[nsimPV]/I");
189	>	pvtxtree_->Branch("simx",&simx_,"simx[nsimPV]/D");
190	>	pvtxtree_->Branch("simy",&simy_,"simy[nsimPV]/D");
191	>	pvtxtree_->Branch("simz",&simz_,"simz[nsimPV]/D");
192	>	pvtxtree_->Branch("simptsq",&simptsq_,"simptsq[nsimPV]/D");
193	>
194	>	// For pvtxs with all the tracks
195	>	pvtxtree_->Branch("nrecPV_mct",&nrecPV_mct_,"nrecPV_mct/I");
196	>	pvtxtree_->Branch("deltax_mct",&deltax_mct_,"deltax_mct[nrecPV_mct]/D");
197	>	pvtxtree_->Branch("deltay_mct",&deltay_mct_,"deltay_mct[nrecPV_mct]/D");
198	>	pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
199	>	pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
200	>	pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
201	>	pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
202	>	// For pvtxs with splittracks1
203	>	pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
204	>	pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
205	>	pvtxtree_->Branch("deltay_spl1_mct",&deltay_spl1_mct_,"deltay_spl1_mct[nrecPV_spl1_mct]/D");
206	>	pvtxtree_->Branch("deltaz_spl1_mct",&deltaz_spl1_mct_,"deltaz_spl1_mct[nrecPV_spl1_mct]/D");
207	>	pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
208	>	pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
209	>	pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
210	>	// For pvtxs with splittracks1
211	>	pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
212	>	pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
213	>	pvtxtree_->Branch("deltay_spl2_mct",&deltay_spl2_mct_,"deltay_spl2_mct[nrecPV_spl2_mct]/D");
214	>	pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
215	>	pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
216	>	pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
217	>	pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
218	>	}
219	>	}
220	>
221	>	setRootStyle();
222	>
223	>	// Create a root file for histograms
224	>	theFile = new TFile(histoFileName_.c_str(), "RECREATE");
225	>	// make diretories
226	>	theFile->mkdir("Summary");
227	>	theFile->cd();
228	>	theFile->mkdir("Others");
229	>	theFile->cd();
230	>
231	>	// Book Histograms:
232	>	h_pvtrk   = new PVHistograms();
233	>	h_pixvtx  = new PVHistograms();
234	>	h_misc    = new PVHistograms();
235	>	h_summary = new PVHistograms();
236	>	h_others  = new PVHistograms();
237	>
238	>	for(int i=0;i<3;i++) {
239	>	if(i == 0) h_pvtrk->Init("pvTrk");
240	>	else {
241	>	stringstream ss;
242	>	ss << i;
243	>	h_pvtrk->Init("pvTrk",ss.str(),"spl");
244	>	}
245	>	}
246	>	h_pixvtx->Init("pixVtx");
247	>	h_misc->Init("misc");
248
249	+	// Book MC only plots
250		if (!realData_) {
251	<	h_deltax       = fs->make<TH1D>("deltax","x-distance (Rec - Sim)",800,-0.04,0.04);
252	<	h_deltay       = fs->make<TH1D>("deltay","y-distance (Rec - Sim)",800,-0.04,0.04);
253	<	h_deltaz       = fs->make<TH1D>("deltaz","z-distance (Rec - Sim)",800,-0.04,0.04);
254	<	h_pullx        = fs->make<TH1D>("pullx","x-pull (Rec - Sim)",1000,-25.,25.);
255	<	h_pully        = fs->make<TH1D>("pully","y-pull (Rec - Sim)",1000,-25.,25.);
256	<	h_pullz        = fs->make<TH1D>("pullz","z-pull (Rec - Sim)",1000,-25.,25.);
257	<
258	<	// Summary of analysis:
259	<	if (analyze_) {
260	<	h_resx_nTrk    = fs->make<TH2D>("resx_ntrk","x-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
261	<	h_resx_nTrk->SetMarkerStyle(21);
262	<	h_resx_nTrk->SetMarkerColor(4);
263	<	h_resx_nTrk->GetXaxis()->SetTitle("Num of tracks");
264	<	h_resx_nTrk->GetYaxis()->SetTitle("#mum");
265	<	h_resy_nTrk    = fs->make<TH2D>("resy_ntrk","y-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
266	<	h_resy_nTrk->SetMarkerStyle(21);
267	<	h_resy_nTrk->SetMarkerColor(4);
268	<	h_resy_nTrk->GetXaxis()->SetTitle("Num of tracks");
110	<	h_resy_nTrk->GetYaxis()->SetTitle("#mum");
111	<	h_resz_nTrk    = fs->make<TH2D>("resz_ntrk","z-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
112	<	h_resz_nTrk->SetMarkerStyle(21);
113	<	h_resz_nTrk->SetMarkerColor(4);
114	<	h_resz_nTrk->GetXaxis()->SetTitle("Num of tracks");
115	<	h_resz_nTrk->GetYaxis()->SetTitle("#mum");
116	<	h_pullx_nTrk   = fs->make<TH2D>("pullx_ntrk","x-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
117	<	h_pullx_nTrk->SetMarkerStyle(21);
118	<	h_pullx_nTrk->SetMarkerColor(4);
119	<	h_pullx_nTrk->SetBit(TH1::kCanRebin);
120	<	h_pullx_nTrk->GetXaxis()->SetTitle("Num of tracks");
121	<	h_pully_nTrk   = fs->make<TH2D>("pully_ntrk","y-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
122	<	h_pully_nTrk->SetMarkerStyle(21);
123	<	h_pully_nTrk->SetMarkerColor(4);
124	<	h_pully_nTrk->SetBit(TH1::kCanRebin);
125	<	h_pully_nTrk->GetXaxis()->SetTitle("Num of tracks");
126	<	h_pullz_nTrk   = fs->make<TH2D>("pullz_ntrk","z-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
127	<	h_pullz_nTrk->SetMarkerStyle(21);
128	<	h_pullz_nTrk->SetMarkerColor(4);
129	<	h_pullz_nTrk->SetBit(TH1::kCanRebin);
130	<	h_pullz_nTrk->GetXaxis()->SetTitle("Num of tracks");
251	>	h_gen = new PVHistograms();
252	>	h_gen->Init("generator");
253	>	}
254	>	if (analyze_) h_ana = new PVHistograms();
255	>
256	>	//Book histograms sensitive to data/mc
257	>	for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
258	>	string suffix;
259	>	edm::LogInfo("Debug")<<"datamode = "<< *it<<endl;
260	>	switch(*it) {
261	>	case 0: suffix = "";
262	>	break;
263	>	case 1: suffix = "_spl1_mct";
264	>	break;
265	>	case 2: suffix = "_spl2_mct";
266	>	break;
267	>	case 3: suffix = "_mct";
268	>	break;
269		}
270	+	h_summary->Init("summary",suffix);
271	+
272	+	// Book residual, error and pull histograms for each nTrk bin
273		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
274		stringstream ss;
275		ss << ntrk;
276	<	string suffix = ss.str();
136	<	h["resx_"+suffix]  = fs->make<TH1D> (TString("deltax_"+suffix),"x-distance (Rec - Sim)",100,-0.02,0.02);
137	<	h["resx_"+suffix]->GetXaxis()->SetTitle("cm");
138	<	h["resy_"+suffix]  = fs->make<TH1D> (TString("deltay_"+suffix),"y-distance (Rec - Sim)",100,-0.02,0.02);
139	<	h["resy_"+suffix]->GetXaxis()->SetTitle("cm");
140	<	h["resz_"+suffix]  = fs->make<TH1D> (TString("deltaz_"+suffix),"z-distance (Rec - Sim)",100,-0.02,0.02);
141	<	h["resz_"+suffix]->GetXaxis()->SetTitle("cm");
142	<	h["pullx_"+suffix] = fs->make<TH1D> (TString("pullx_"+suffix),"x-pull (Rec - Sim)",100,-10.,10.);
143	<	h["pully_"+suffix] = fs->make<TH1D> (TString("pully_"+suffix),"y-pull (Rec - Sim)",100,-10.,10.);
144	<	h["pullz_"+suffix] = fs->make<TH1D> (TString("pullz_"+suffix),"z-pull (Rec - Sim)",100,-10.,10.);
145	<	suffix.clear();
276	>	h_others->Init("others",suffix,ss.str());
277		}
147	–	}
148	–	}
278
279	+	// Book residual and pull histograms only when analyzeOntheFly is enabled
280	+	if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
281	+	suffix.clear();
282	+	} // End of Book histograms sensitive to data/mc
283	+
284	+	}
285
286		PVStudy::~PVStudy()
287		{
288	<
289	<	// do anything here that needs to be done at desctruction time
290	<	// (e.g. close files, deallocate resources etc.)
288	>	theFile->cd();
289	>	theFile->cd("Summary");
290	>	h_pvtrk->Save();
291	>	h_pixvtx->Save();
292	>	h_gen->Save();
293	>	h_misc->Save();
294	>	h_summary->Save();
295	>	if (analyze_)
296	>	h_ana->Save();
297	>	theFile->cd();
298	>	theFile->cd("Others");
299	>	h_others->Save();
300
301	+	theFile->Close();
302		}
303
304	<
304	>	void PVStudy::setRootStyle() {
305	>	//
306	>	gROOT->SetStyle("Plain");
307	>	gStyle->SetFillColor(1);
308	>	gStyle->SetOptDate();
309	>	//   gStyle->SetOptStat(1111110);
310	>	//   gStyle->SetOptFit(0111);
311	>	//   gStyle->SetPadTickX(1);
312	>	//   gStyle->SetPadTickY(1);
313	>	gStyle->SetMarkerSize(0.5);
314	>	gStyle->SetMarkerStyle(8);
315	>	gStyle->SetGridStyle(3);
316	>	//gStyle->SetPadGridX(1);
317	>	//gStyle->SetPadGridY(1);
318	>	gStyle->SetPaperSize(TStyle::kA4);
319	>	gStyle->SetStatW(0.25);             // width of statistics box; default is 0.19
320	>	//   gStyle->SetStatH(0.10);             // height of statistics box; default is 0.1
321	>	gStyle->SetStatFormat("6.4g");      // leave default format for now
322	>	gStyle->SetTitleSize(0.055, "");    // size for pad title; default is 0.02
323	>	gStyle->SetLabelSize(0.03, "XYZ");  // size for axis labels; default is 0.04
324	>	gStyle->SetStatFontSize(0.08);      // size for stat. box
325	>	gStyle->SetTitleFont(32, "XYZ");    // times-bold-italic font (p. 153) for axes
326	>	gStyle->SetTitleFont(32, "");       // same for pad title
327	>	gStyle->SetLabelFont(32, "XYZ");    // same for axis labels
328	>	gStyle->SetStatFont(32);            // same for stat. box
329	>	gStyle->SetLabelOffset(0.006, "Y"); // default is 0.005
330	>	//
331	>	return;
332	>	}
333		//
334		// member functions
335		//
#	Line 165 | Line 338 | std::vector<PVStudy::simPrimaryVertex> P
338		std::vector<PVStudy::simPrimaryVertex> simpv;
339		const HepMC::GenEvent* evt=evtMC->GetEvent();
340		if (evt) {
341	<	if(verbose_) {
342	<	cout << "process id " << evt->signal_process_id()<<endl;
343	<	cout << "signal process vertex " << ( evt->signal_process_vertex() ?
344	<	evt->signal_process_vertex()->barcode() : 0 )   <<endl;
172	<	cout << "number of vertices " << evt->vertices_size() << endl;
173	<	}
341	>	edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
342	>	edm::LogInfo("SimPVs") << "[getSimPVs] signal process vertex " << ( evt->signal_process_vertex() ?
343	>	evt->signal_process_vertex()->barcode() : 0 ) <<endl;
344	>	edm::LogInfo("SimPVs") << "[getSimPVs] number of vertices " << evt->vertices_size() << endl;
345
346	<	int idx=0;
346	>	int idx=0; int npv=0;
347		for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
348	<	vitr != evt->vertices_end(); ++vitr )
349	<	{ // loop for vertex ...
350	<	HepMC::FourVector pos = (*vitr)->position();
351	<	//HepLorentzVector pos = (*vitr)->position();
352	<	if (pos.t()>0) { continue;} // ?
353	<
354	<	bool hasMotherVertex=false;
355	<	//      if(verbose_) cout << "mothers" << endl;
356	<	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
357	<	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
358	<	HepMC::GenVertex * mv=(*mother)->production_vertex();
359	<	if (mv) {hasMotherVertex=true;}
360	<	//      if(verbose_) {
361	<	//        cout << "\t";
191	<	//        (*mother)->print();
192	<	//      }
193	<	}
194	<
195	<	if(hasMotherVertex) {continue;}
196	<
197	<	// could be a new vertex, check  all primaries found so far to avoid multiple entries
198	<	const double mm2cm=0.1;
199	<	simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
200	<	simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
201	<	for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
202	<	v0!=simpv.end(); v0++){
203	<	if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
204	<	vp=&(*v0);
205	<	break;
348	>	vitr != evt->vertices_end(); ++vitr ) { // loop for vertex ...
349	>	HepMC::FourVector pos = (*vitr)->position();
350	>	//HepLorentzVector pos = (*vitr)->position();
351	>
352	>	// t component of PV:
353	>	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
354	>	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
355	>	//        edm::LogInfo("SimPVs") << "Status = " << (*mother)->status() << endl;
356	>	HepMC::GenVertex * mv=(*mother)->production_vertex();
357	>	if( ((*mother)->status() == 3) && (!mv)) {
358	>	//        edm::LogInfo("SimPVs") << "npv= " << npv << endl;
359	>	if (npv == 0) {
360	>	h_gen->Fill1d("genPart_cT", pos.t()); // mm
361	>	h_gen->Fill1d("genPart_T", pos.t()/299.792458); // ns
362		}
363	+	npv++;
364		}
365	+	}
366	+	//       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
367	+
368	+	bool hasMotherVertex=false;
369	+	if (verbose_) cout << "[getSimPVs] mothers of vertex[" << ++idx << "]: " << endl;
370	+	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
371	+	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
372	+	HepMC::GenVertex * mv=(*mother)->production_vertex();
373	+	// if (verbose_) cout << "Status = " << (*mother)->status() << endl;
374	+	if (mv) {
375	+	hasMotherVertex=true;
376	+	if(!verbose_) break; //if verbose_, print all particles of gen vertices
377	+	}
378	+	if(verbose_) {
379	+	cout << "\t";
380	+	(*mother)->print();
381	+	}
382	+	}
383
384	<	if(!vp){
385	<	// this is a new vertex
386	<	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
387	<	simpv.push_back(sv);
388	<	vp=&simpv.back();
389	<	}else{
390	<	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
391	<	}
392	<	vp->genVertex.push_back((*vitr)->barcode());
393	<	// collect final state descendants
394	<	for ( HepMC::GenVertex::particle_iterator
395	<	daughter  = (*vitr)->particles_begin(HepMC::descendants);
396	<	daughter != (*vitr)->particles_end(HepMC::descendants);
397	<	++daughter ) {
398	<	if (isFinalstateParticle(*daughter)){
399	<	if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
400	<	== vp->finalstateParticles.end()){
401	<	vp->finalstateParticles.push_back((*daughter)->barcode());
402	<	HepMC::FourVector m=(*daughter)->momentum();
403	<	// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
404	<	// but adding FourVectors seems not to be foreseen
405	<	vp->ptot.setPx(vp->ptot.px()+m.px());
406	<	vp->ptot.setPy(vp->ptot.py()+m.py());
407	<	vp->ptot.setPz(vp->ptot.pz()+m.pz());
408	<	vp->ptot.setE(vp->ptot.e()+m.e());
409	<	vp->ptsq+=(m.perp())*(m.perp());
410	<	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
411	<	vp->nGenTrk++;
412	<	}
384	>	if(hasMotherVertex) continue;
385	>
386	>	// could be a new vertex, check  all primaries found so far to avoid multiple entries
387	>	const double mm2cm=0.1;
388	>	simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
389	>	simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
390	>	for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
391	>	v0!=simpv.end(); v0++){
392	>	if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
393	>	vp=&(*v0);
394	>	break;
395	>	}
396	>	}
397	>
398	>	if(!vp){
399	>	// this is a new vertex
400	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
401	>	simpv.push_back(sv);
402	>	vp=&simpv.back();
403	>	}else{
404	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
405	>	}
406	>	vp->genVertex.push_back((*vitr)->barcode());
407	>	// collect final state descendants
408	>	for ( HepMC::GenVertex::particle_iterator daughter  = (*vitr)->particles_begin(HepMC::descendants);
409	>	daughter != (*vitr)->particles_end(HepMC::descendants);
410	>	++daughter ) {
411	>	if (isFinalstateParticle(*daughter)){
412	>	if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
413	>	== vp->finalstateParticles.end()){
414	>	vp->finalstateParticles.push_back((*daughter)->barcode());
415	>	HepMC::FourVector m=(*daughter)->momentum();
416	>	// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
417	>	// but adding FourVectors seems not to be foreseen
418	>	vp->ptot.setPx(vp->ptot.px()+m.px());
419	>	vp->ptot.setPy(vp->ptot.py()+m.py());
420	>	vp->ptot.setPz(vp->ptot.pz()+m.pz());
421	>	vp->ptot.setE(vp->ptot.e()+m.e());
422	>	vp->ptsq+=(m.perp())*(m.perp());
423	>	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
424	>	vp->nGenTrk++;
425		}
426		}
427		}
428	<	idx++;
429	<	}
428	>	}//loop MC vertices daughters
429	>	}//loop MC vertices
430		}
431		return simpv;
432		}
#	Line 272 | Line 459 | std::vector<PVStudy::simPrimaryVertex> P
459		ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
460		}
461		}
462	<	primary =  ( gv->position().t()==0);
462	>	//primary =  ( gv->position().t()==0 );
463	>	primary = true;
464	>	h_gen->Fill1d("genPart_cT", gv->position().t()); // mm
465	>	h_gen->Fill1d("genPart_T", gv->position().t()/299.792458); // ns
466	>
467		//resonance= ( gp->mother() && isResonance(gp->mother()));  // in CLHEP/HepMC days
468		// no more mother pointer in the improved HepMC GenParticle
469		resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
#	Line 319 | Line 510 | void
510		PVStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
511		{
512		using namespace edm;
513	<	using reco::TrackCollection;
513	>	using namespace reco;
514	>
515	>	edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
516	>	//=======================================================
517	>	// Initialize Root-tuple variables if needed
518	>	//=======================================================
519	>	//if(saventuple_)
520	>	SetVarToZero();
521	>
522	>	//=======================================================
523	>	// Track accessors
524	>	//=======================================================
525	>	//trackCollection
526	>	static const reco::TrackCollection s_empty_trackColl;
527	>	const reco::TrackCollection *trackColl = &s_empty_trackColl;
528	>	edm::Handle<reco::TrackCollection> trackCollectionHandle;
529	>	iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle);
530	>	if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
531	>	trackColl = trackCollectionHandle.product();
532	>	} else {
533	>	edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
534	>	}
535
536	<	ftree_->SetBranchAddress("residual",&fres);
537	<	ftree_->SetBranchAddress("error",&ferror);
538	<	ftree_->SetBranchAddress("nTrk",&fntrk);
536	>	//splitTrackCollection1
537	>	static const reco::TrackCollection s_empty_splitTrackColl1;
538	>	const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
539	>	edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
540	>	iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle);
541	>	if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
542	>	splitTrackColl1 = splitTrackCollection1Handle.product();
543	>	} else {
544	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
545	>	}
546
547	<	Handle<TrackCollection> tracks;
548	<	iEvent.getByLabel(tracksLabel_,tracks);
547	>	//splitTrackCollection2
548	>	static const reco::TrackCollection s_empty_splitTrackColl2;
549	>	const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
550	>	edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
551	>	iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle);
552	>	if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
553	>	splitTrackColl2 = splitTrackCollection2Handle.product();
554	>	} else {
555	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
556	>	}
557
331	–	Handle<reco::VertexCollection> recVtxs;
332	–	iEvent.getByLabel(vtxSample_, recVtxs);
558
559	+	//=======================================================
560	+	// Fill trackparameters of the input tracks to pvtx fitter
561	+	//=======================================================
562	+	edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
563	+	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype)
564	+	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
565	+	fillTrackHisto(trackColl, 0);
566	+	fillTrackHisto(splitTrackColl1, 1);
567	+	fillTrackHisto(splitTrackColl2, 2);
568	+	edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
569	+
570	+	//=======================================================
571	+	// PVTX accessors
572	+	//=======================================================
573	+	//vertexCollection
574	+	static const reco::VertexCollection s_empty_vertexColl;
575	+	const reco::VertexCollection *vertexColl = &s_empty_vertexColl;
576	+	edm::Handle<reco::VertexCollection> vertexCollectionHandle;
577	+	iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle);
578	+	if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
579	+	vertexColl = vertexCollectionHandle.product();
580	+	} else {
581	+	edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
582	+	}
583	+
584	+	//splitVertexCollection1
585	+	static const reco::VertexCollection s_empty_splitVertexColl1;
586	+	const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
587	+	edm::Handle<reco::VertexCollection> splitVertexCollection1Handle;
588	+	iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle);
589	+	if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
590	+	splitVertexColl1 = splitVertexCollection1Handle.product();
591	+	} else {
592	+	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
593	+	}
594	+
595	+	//splitVertexCollection2
596	+	static const reco::VertexCollection s_empty_splitVertexColl2;
597	+	const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
598	+	edm::Handle<reco::VertexCollection> splitVertexCollection2Handle;
599	+	iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle);
600	+	if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
601	+	splitVertexColl2 = splitVertexCollection2Handle.product();
602	+	} else {
603	+	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
604	+	}
605	+
606	+	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track and vertex collections"<<endl;
607	+
608	+	//=======================================================
609	+	// MC simvtx accessor
610	+	//=======================================================
611		if (!realData_) {
612	<	Handle<SimVertexContainer> simVtxs;
612	>	edm::Handle<SimVertexContainer> simVtxs;
613		iEvent.getByLabel( simG4_, simVtxs);
614
615	<	Handle<SimTrackContainer> simTrks;
615	>	edm::Handle<SimTrackContainer> simTrks;
616		iEvent.getByLabel( simG4_, simTrks);
617		}
618
619	<	if(verbose_) printRecVtxs(recVtxs);
620	<
619	>	//=======================================================
620	>	// GET PDT
621	>	//=======================================================
622		try{
623		iSetup.getData(pdt);
624		}catch(const Exception&){
625	<	cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
625	>	edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
626	>	}
627	>
628	>	//=======================================================
629	>	// GET pixelVertices
630	>	//=======================================================
631	>	static const reco::VertexCollection s_empty_pixelVertexColl;
632	>	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
633	>	edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
634	>	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
635	>	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
636	>	pixelVertexColl = pixelVertexCollectionHandle.product();
637	>	} else {
638	>	edm::LogInfo("Debug") << "[PVStudy] pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
639	>	}
640	>
641	>	//=======================================================
642	>	// Fill pixelVertices related histograms
643	>	//=======================================================
644	>	nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
645	>	if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
646	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
647	>	fillTrackHistoInPV(pixelVertexColl, 4, false, true);
648	>	h_pixvtx->Fill1d("dzErr_pxlpvtx", pixelVertexColl->begin()->zError());
649	>	// Get the dZ error of the tracks in the leading pixelVertexColl
650	>	for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
651	>	t!= (pixelVertexColl->begin())->tracks_end(); t++) {
652	>
653	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
654	>	//        h_pvtrk->Fill1d("trkPtPV", 0.);
655	>	}
656	>	else
657	>	h_pixvtx->Fill1d("trkdzErr_pxlpvtx", (**t).dzError());
658	>	}
659	>	// Fill the track->dz() in the PV difference from first pixelpvtx
660	>	for(reco::Vertex::trackRef_iterator t = (vertexColl->begin())->tracks_begin();
661	>	t!= (vertexColl->begin())->tracks_end(); t++) {
662	>	// illegal charge
663	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
664	>	//        h_pvtrk->Fill1d("trkPtPV", 0.);
665	>	}
666	>	else {
667	>	if(pixelVertexColl->size()>0) {
668	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz", (**t).dz() -  pixelVertexColl->begin()->z());
669	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_pxlpvtxdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
670	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_trkdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
671	>	h_pixvtx->Fill1d("trkdzErr_pvtx", (**t).dzError());
672	>	}
673	>	}
674	>	}
675	>	}
676	>
677	>	//=======================================================
678	>	// Fill number of reconstructed vertices
679	>	//=======================================================
680	>
681	>	edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
682	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
683	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
684	>	edm::LogInfo("Debug")<<"[PVStudy] Start filling pvtx parameters."<<endl;
685	>	if (verbose_) {
686	>	printRecVtxs(vertexCollectionHandle);
687	>	printRecVtxs(splitVertexCollection1Handle);
688	>	printRecVtxs(splitVertexCollection2Handle);
689		}
690	+	nrecPV_ = int (vertexColl->size());
691	+	nrecPV1_spl_ = int (splitVertexColl1->size());
692	+	nrecPV2_spl_ = int (splitVertexColl2->size());
693	+
694	+	h_pvtrk->Fill1d("nrecPV", nrecPV_);
695	+	h_pvtrk->Fill1d("nrecPV1_spl", nrecPV1_spl_);
696	+	h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
697	+	h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
698	+
699	+	//=================================================================
700	+	// Fill track parameter ntuple/hist for tracks used in recoVertices
701	+	//=================================================================
702	+
703	+	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
704	+	if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
705	+	fillTrackHistoInPV(vertexColl, 0, true, true);
706	+
707	+	if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
708	+	fillTrackHistoInPV(splitVertexColl1, 1, false, true);
709	+
710	+	if(splitVertexColl1->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
711	+	fillTrackHistoInPV(splitVertexColl2, 2, false, true);
712	+
713	+
714	+	//=======================================================
715	+	// Compare offlinePrimaryVertices with pixelVertices
716	+	//=======================================================
717	+	if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
718	+	&& (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {
719	+	h_pixvtx->Fill1d("nrecPV_minus_nrecPxlPV", double (nrecPV_ - nrecPV_pxlpvtx_));
720	+	// difference in reconstructed position of the leading pvtx
721	+	//edm::LogInfo("Debug")<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;
722	+	//edm::LogInfo("Debug")<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
723	+	h_pixvtx->Fill1d("recxPV_minus_recxPxlPV", recx_[0] - recx_pxlpvtx_[0]);
724	+	h_pixvtx->Fill1d("recyPV_minus_recyPxlPV", recy_[0] - recy_pxlpvtx_[0]);
725	+	h_pixvtx->Fill1d("reczPV_minus_reczPxlPV", recz_[0] - recz_pxlpvtx_[0]);
726	+	}
727	+
728
729	+
730	+	//==========================================================
731	+	// Fill secondary/primary min separations for multi-vertices
732	+	//==========================================================
733	+
734	+	if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {
735	+
736	+	double min_xsep = 9999.0;
737	+	double min_ysep = 9999.0;
738	+	double min_zsep = 9999.0;
739	+	double min_xsep_sign = 9999.0;
740	+	double min_ysep_sign = 9999.0;
741	+	double min_zsep_sign = 9999.0;
742	+	double min_ntrksep = 9999.0;
743	+	double min_sumpt2sep = 9999999.0;
744	+
745	+	edm::LogInfo("Debug")<<"[PVStudy] leading pvtx z = "<< vertexColl->begin()->z() <<endl;
746	+
747	+	// Looping through the secondary vertices to find the mininum separation to the primary
748	+	for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
749	+	v!=vertexColl->end(); ++v) {
750	+	if(v->isValid() && ! v->isFake()) {
751	+	// xsep
752	+	if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
753	+	min_xsep = fabs(v->x()- vertexColl->begin()->x());
754	+	// ysep
755	+	if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
756	+	min_ysep = fabs(v->y()- vertexColl->begin()->y());
757	+	// zsep
758	+	if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
759	+	min_zsep = fabs(v->z()- vertexColl->begin()->z());
760	+	// xsep_sign
761	+	double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
762	+	if(xsep_sign < min_xsep_sign)
763	+	min_xsep_sign = xsep_sign;
764	+	// ysep_sign
765	+	double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
766	+	if(ysep_sign < min_ysep_sign)
767	+	min_ysep_sign = ysep_sign;
768	+	// zsep_sign
769	+	double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
770	+	if(zsep_sign < min_zsep_sign)
771	+	min_zsep_sign = zsep_sign;
772	+	// ntrksep
773	+	if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
774	+	min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
775	+	// sumpt2sep
776	+	if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
777	+	min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
778	+	}
779	+	}
780	+	h_misc->Fill1d("min_xsep", min_xsep);
781	+	h_misc->Fill1d("min_ysep", min_ysep);
782	+	h_misc->Fill1d("min_zsep", min_zsep);
783	+	h_misc->Fill1d("min_xsep_sign", min_xsep_sign);
784	+	h_misc->Fill1d("min_ysep_sign", min_ysep_sign);
785	+	h_misc->Fill1d("min_zsep_sign", min_zsep_sign);
786	+	h_misc->Fill1d("min_ntrksep", min_ntrksep);
787	+	h_misc->Fill1d("min_sumpt2sep", min_sumpt2sep);
788	+
789	+	min_zsep_ = min_zsep;
790	+	min_ntrksep_ = min_ntrksep;
791	+	min_sumpt2sep_ = min_sumpt2sep;
792	+
793	+
794	+	} // End of  if(vertexColl->size()>1) {
795	+
796	+	edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
797	+
798	+	//=======================================================
799	+	//           PrimaryVertex Matching
800	+	// 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
801	+	// 2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
802	+	// Assume the first match is the primary vertex,
803	+	// since vertexColl are sorted in decreasing order of sum(pT)
804	+	//=======================================================
805	+
806	+	edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
807	+	reco::VertexCollection s_empty_matchedVertexColl1;
808	+	reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
809	+	matchedVertexColl1->reserve(splitVertexColl1->size());
810	+	reco::VertexCollection s_empty_matchedVertexColl2;
811	+	reco::VertexCollection *matchedVertexColl2 = &s_empty_matchedVertexColl2;
812	+	matchedVertexColl2->reserve(splitVertexColl2->size());
813	+
814	+	bool stopmatching = false;
815	+	for (size_t ivtx = 0; ivtx<splitVertexCollection1Handle->size() && !stopmatching; ++ivtx) {
816	+	reco::VertexRef recvtx1(splitVertexCollection1Handle, ivtx);
817	+	if( !(recvtx1->isValid()) || recvtx1->isFake()) break;
818	+	for (size_t jvtx = ivtx; jvtx < splitVertexCollection2Handle->size(); ++jvtx) {
819	+	//------------------------------------------------------------------------
820	+	//== If only considering matching the first vertex of the splitVertexColl
821	+	//------------------------------------------------------------------------
822	+	if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
823	+	reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
824	+	if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
825	+	edm::LogInfo("Debug")<<"[PVStudy] Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
826	+	edm::LogInfo("Debug")<<"[PVStudy] recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
827	+	edm::LogInfo("Debug")<<"[PVStudy] recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
828	+
829	+	if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
830	+	edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
831	+
832	+	int nTrkPV1 = recvtx1->tracksSize();
833	+	int nTrkPV2 = recvtx2->tracksSize();
834	+	double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
835	+	h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
836	+	h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
837	+	if(fabs(ntrkreldiff)<ntrkdiffcut_) {
838	+	edm::LogInfo("Debug")<<"[PVStudy] (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
839	+
840	+	matchedVertexColl1->push_back(*recvtx1);
841	+	matchedVertexColl2->push_back(*recvtx2);
842	+	stopmatching = true; // stop the matching when the first match is found!
843	+	break;
844	+	}
845	+	else
846	+	edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
847	+	}
848	+	else {
849	+	edm::LogInfo("Debug")<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
850	+	edm::LogInfo("Debug")<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
851	+	edm::LogInfo("Debug")<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
852	+	}
853	+	}
854	+	}
855	+	edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
856	+
857	+	//=======================================================
858	+	//  Analyze matchedVertexColls
859	+	//=======================================================
860	+	edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
861	+
862	+	// Compare the reconstructed vertex position and calculate resolution/pulls
863	+	if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
864	+	//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
865	+	fillTrackHistoInPV(matchedVertexColl1, 1, true, false);
866	+	fillTrackHistoInPV(matchedVertexColl2, 2, true, false);
867	+
868	+	reco::VertexCollection::const_iterator v1;
869	+	reco::VertexCollection::const_iterator v2;
870	+	nrecPV_twovtx_ = 0;
871	+	for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
872	+	v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
873	+	++v1, ++v2) {
874	+	fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
875	+	fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
876	+	fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
877	+	ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
878	+	ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
879	+	ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
880	+	fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
881	+
882	+	h_summary->Fill1d("deltax", fres_[0] );
883	+	h_summary->Fill1d("deltay", fres_[1] );
884	+	h_summary->Fill1d("deltaz", fres_[2] );
885	+	h_summary->Fill1d("pullx", fres_[0]/ferror_[0]);
886	+	h_summary->Fill1d("pully", fres_[1]/ferror_[1]);
887	+	h_summary->Fill1d("pullz", fres_[2]/ferror_[2]);
888	+	h_summary->Fill1d("errPVx", ferror_[0] );
889	+	h_summary->Fill1d("errPVy", ferror_[1] );
890	+	h_summary->Fill1d("errPVz", ferror_[2] );
891	+	pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
892	+	error(ferror_[0], ferror_[1], ferror_[2]),
893	+	fntrk_));
894	+	// Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
895	+	fillHisto(res(fres_[0], fres_[1], fres_[2]),
896	+	error(ferror_[0], ferror_[1], ferror_[2]),
897	+	fntrk_,0);
898	+	// Fill the ntuple variables
899	+	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
900	+	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
901	+	sumptsq1_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v1);
902	+	sumptsq2_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v2);
903	+	nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
904	+	deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
905	+	deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
906	+	deltaz_twovtx_[nrecPV_twovtx_] = fres_[2];
907	+	errx_twovtx_[nrecPV_twovtx_] = ferror_[0];
908	+	erry_twovtx_[nrecPV_twovtx_] = ferror_[1];
909	+	errz_twovtx_[nrecPV_twovtx_] = ferror_[2];
910	+	pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
911	+	pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
912	+	pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];
913	+	nrecPV_twovtx_++;
914	+	} // End of analyzing res/pull
915	+	} // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
916	+	else
917	+	edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
918	+
919	+	edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
920	+
921	+	//=======================================================
922	+	// Fill simulated vertices
923	+	//=======================================================
924		if (!realData_) {
925		bool MC=false;
926		Handle<HepMCProduct> evtMC;
927		iEvent.getByLabel("generator",evtMC);
928		if (!evtMC.isValid()) {
929		MC=false;
930	<	if(verbose_){
357	<	cout << "no HepMCProduct found"<< endl;
358	<	}
930	>	edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
931		} else {
932	<	if(verbose_){
361	<	cout << "generator HepMCProduct found"<< endl;
362	<	}
932	>	edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
933		MC=true;
934		}
935	<
935	>
936		if(MC){
937		// make a list of primary vertices:
938		std::vector<simPrimaryVertex> simpv;
939		simpv=getSimPVs(evtMC,"");
940		//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
941	<
942	<	int nsimtrk=0;
943	<	int nrectrk=0;
941	>	h_gen->Fill1d("nsimPV", simpv.size());
942	>	nsimPV_ = simpv.size();
943	>
944	>	int isimpv = 0;
945		for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
946	<	vsim!=simpv.end(); vsim++){
947	<	nsimtrk+=vsim->nGenTrk;
948	<	// look for a matching reconstructed vertex
949	<	vsim->recVtx=NULL;
950	<	for(reco::VertexCollection::const_iterator vrec=recVtxs->begin();
951	<	vrec!=recVtxs->end(); ++vrec){
952	<	nrectrk=vrec->tracksSize();
953	<	if(verbose_){
954	<	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
955	<	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
956	<	}
957	<	if ( matchVertex(*vsim,*vrec) ){
958	<	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
959	<	if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
960	<	|| (!vsim->recVtx) )
961	<	{
962	<	vsim->recVtx=&(*vrec);
963	<	}
964	<	}
394	<	}
395	<	if (vsim->recVtx){
396	<	if(verbose_) {
397	<	cout <<"primary matched " << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
398	<	}
399	<	fres[0] = vsim->recVtx->x()-vsim->x;
400	<	fres[1] = vsim->recVtx->y()-vsim->y;
401	<	fres[2] = vsim->recVtx->z()-vsim->z;
402	<	ferror[0] = vsim->recVtx->xError();
403	<	ferror[1] = vsim->recVtx->yError();
404	<	ferror[2] = vsim->recVtx->zError();
405	<	fntrk = nrectrk;
406	<	ftree_->Fill();
407	<
408	<	h_deltax->Fill( vsim->recVtx->x()-vsim->x );
409	<	h_deltay->Fill( vsim->recVtx->y()-vsim->y );
410	<	h_deltaz->Fill( vsim->recVtx->z()-vsim->z );
411	<	h_pullx->Fill( (vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError() );
412	<	h_pully->Fill( (vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError() );
413	<	h_pullz->Fill( (vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError() );
414	<	pvinfo.push_back(PVStudy::PVInfo(res(vsim->recVtx->x()-vsim->x,
415	<	vsim->recVtx->y()-vsim->y,
416	<	vsim->recVtx->z()-vsim->z),
417	<	error(vsim->recVtx->xError(),
418	<	vsim->recVtx->yError(),
419	<	vsim->recVtx->zError()),
420	<	nrectrk));
421	<	// Fill histo according to its track multiplicity
422	<	fillHisto(res(vsim->recVtx->x()-vsim->x,
423	<	vsim->recVtx->y()-vsim->y,
424	<	vsim->recVtx->z()-vsim->z),
425	<	pull((vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError(),
426	<	(vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError(),
427	<	(vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError()),
428	<	nrectrk);
429	<	}
430	<	else{  // no rec vertex found for this simvertex
431	<	if(verbose_) {
432	<	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
433	<	}
434	<	}
435	<	}
436	<	}
437	<	}// With MC
438	<
439	<	// Loop RecVtxs and find matched SimVtxs
440	<	for(reco::VertexCollection::const_iterator v=recVtxs->begin();
441	<	v!=recVtxs->end(); ++v){
442	<	try {
443	<	for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
444	<	t!=v->tracks_end(); t++) {
445	<	// illegal charge
446	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
447	<	//      h_trkPtPV->Fill(0.);
448	<	}
449	<	else {
450	<	h_trkPtPV->Fill((**t).pt());
451	<	}
452	<	}
453	<	}
454	<	catch (...) {
455	<	// exception thrown when trying to use linked track
456	<	//          h_trkPtPV->Fill(0.);
457	<	}
458	<
459	<	h_nTrkPV->Fill(v->tracksSize());
460	<
461	<	}
462	<
463	<	h_nTrk->Fill(tracks->size());
464	<
465	<	for(TrackCollection::const_iterator itTrack = tracks->begin();
466	<	itTrack != tracks->end();
467	<	++itTrack) {
468	<	int charge = 0;
469	<	charge = itTrack->charge();
470	<	h_trkPt->Fill(itTrack->pt());
946	>	vsim!=simpv.end(); vsim++, isimpv++){
947	>	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
948	>	simx_[isimpv] = vsim->x;
949	>	simy_[isimpv] = vsim->y;
950	>	simz_[isimpv] = vsim->y;
951	>	simptsq_[isimpv] = vsim->ptsq;
952	>
953	>	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
954	>	fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
955	>	fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
956	>	fillMCHisto(vsim, isimpv, vertexColl, 3);
957	>	}
958	>	} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
959	>	} // End of   if (!realData_) {
960	>
961	>	// Finally fill the ftree_
962	>	if(saventuple_) {
963	>	ftree_->Fill();
964	>	pvtxtree_->Fill();
965		}
472	–
966		}
967
475	–
968		// ------------ method called once each job just before starting event loop  ------------
969		void
970		PVStudy::beginJob()
#	Line 482 | Line 974 | PVStudy::beginJob()
974		// ------------ method called once each job just after ending the event loop  ------------
975		void
976		PVStudy::endJob() {
977	<	double sqt2 = sqrt(2.);
978	<	if (!realData_) {
979	<	if (verbose_) cout << "Analyzing PV info" << endl;
980	<	for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
981	<	if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
982	<	PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk);
983	<	if (analyze_) {
984	<	//      if ( ipv.res_.x() > 0 ) h_resx_nTrk->Fill(ntrk,ipv.res_.x()/sqt2, 1.);
985	<	//      if ( ipv.res_.y() > 0 ) h_resy_nTrk->Fill(ntrk,ipv.res_.y()/sqt2, 1.);
986	<	//      if ( ipv.res_.z() > 0 ) h_resz_nTrk->Fill(ntrk,ipv.res_.z()/sqt2, 1.);
987	<	if ( ipv.res_.x() > 0 ) h_resx_nTrk->Fill(ntrk,ipv.res_.x(), 1.);
988	<	if ( ipv.res_.y() > 0 ) h_resy_nTrk->Fill(ntrk,ipv.res_.y(), 1.);
989	<	if ( ipv.res_.z() > 0 ) h_resz_nTrk->Fill(ntrk,ipv.res_.z(), 1.);
990	<	if ( ipv.pull_.x() > 0 ) h_pullx_nTrk->Fill(ntrk,ipv.pull_.x(), 1.);
991	<	if ( ipv.pull_.y() > 0 ) h_pully_nTrk->Fill(ntrk,ipv.pull_.y(), 1.);
992	<	if ( ipv.pull_.z() > 0 ) h_pullz_nTrk->Fill(ntrk,ipv.pull_.z(), 1.);
977	>	edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl;
978	>	// Looping through the datamodes available
979	>	for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
980	>	string suffix;
981	>	edm::LogInfo("Analysis")<<"[endJob] datamode = "<< *it<<endl;
982	>	switch(*it) {
983	>	case 0: suffix = "";
984	>	break;
985	>	case 1: suffix = "_spl1_mct";
986	>	break;
987	>	case 2: suffix = "_spl2_mct";
988	>	break;
989	>	case 3: suffix = "_mct";
990	>	break;
991	>	}
992	>	suffix += "_nTrk";
993	>	if(analyze_) {
994	>	for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
995	>	edm::LogInfo("Analysis") << "[endJob] ntrk = " << ntrk << endl;
996	>	PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
997	>	if ( ipv.res_.x() > 0 ) h_ana->Fill2d("resx"+suffix, ntrk,ipv.res_.x());
998	>	if ( ipv.res_.y() > 0 ) h_ana->Fill2d("resy"+suffix, ntrk,ipv.res_.y());
999	>	if ( ipv.res_.z() > 0 ) h_ana->Fill2d("resz"+suffix, ntrk,ipv.res_.z());
1000	>	if ( ipv.pull_.x() > 0 ) h_ana->Fill2d("pullx"+suffix, ntrk,ipv.pull_.x());
1001	>	if ( ipv.pull_.y() > 0 ) h_ana->Fill2d("pully"+suffix, ntrk,ipv.pull_.y());
1002	>	if ( ipv.pull_.z() > 0 ) h_ana->Fill2d("pullz"+suffix, ntrk,ipv.pull_.z());
1003		}
1004		}
1005	+	suffix.clear();
1006	+	}
1007	+	if(saventuple_) {
1008	+	file_->cd();
1009	+	ftree_->Write();
1010	+	pvtxtree_->Write();
1011	+	file_->Close();
1012		}
504	–	file_->cd();
505	–	ftree_->Write();
506	–	file_->Close();
1013		}
1014
1015	+
1016	+
1017	+	// Match a recovertex with a simvertex
1018	+	// ? Should the cut be parameter dependent?
1019	+	// cut can be within n sigma of the vertex.zerror()
1020		bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
1021		const reco::Vertex & vrec)
1022		{
1023		return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1024		}
1025
1026	+	// Match two reconstructed vertices
1027	+	bool PVStudy::matchVertex( const reco::VertexRef & vrec1,
1028	+	const reco::VertexRef & vrec2,
1029	+	int zsigncut)
1030	+	{
1031	+	double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
1032	+	edm::LogInfo("Debug")<<"[matchVertex] The matching criteria in z is "<<cut<<endl;
1033	+	return (fabs(vrec1->z()-vrec2->z())<cut);
1034	+	}
1035	+
1036	+
1037		bool PVStudy::isResonance(const HepMC::GenParticle * p){
1038		double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
1039	<	if(verbose_) cout << "isResonance   " << p->pdg_id() << " " << ctau << endl;
1039	>	edm::LogInfo("Debug") << "[isResonance] isResonance   " << p->pdg_id() << " " << ctau << endl;
1040		return  ctau >0 && ctau <1e-6;
1041		}
1042
#	Line 583 | Line 1105 | void PVStudy::printSimTrks(const Handle<
1105		}
1106		}
1107
1108	<	void PVStudy::fillHisto(res r, pull p, int ntk){
1108	>	void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1109		stringstream ss;
1110		ss << ntk;
1111	<	string suffix = ss.str();
1111	>	string suffix;
1112	>	if(datamode == 0 )  suffix = "_" + ss.str();
1113	>	if(datamode == 1 )  suffix = "_spl1_mct_" + ss.str();
1114	>	if(datamode == 2 )  suffix = "_spl2_mct_" + ss.str();
1115	>	if(datamode == 3 )  suffix = "_mct_" + ss.str();
1116	>
1117		if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
1118		// Fill histograms of res and pull of ntk
1119	<	h["resx_"+suffix]->Fill(r.x());
1120	<	h["resy_"+suffix]->Fill(r.y());
1121	<	h["resz_"+suffix]->Fill(r.z());
1122	<	h["pullx_"+suffix]->Fill(p.x());
1123	<	h["pully_"+suffix]->Fill(p.y());
1124	<	h["pullz_"+suffix]->Fill(p.z());
1119	>	h_others->Fill1d("deltax"+suffix, r.x());
1120	>	h_others->Fill1d("deltay"+suffix, r.y());
1121	>	h_others->Fill1d("deltaz"+suffix, r.z());
1122	>	h_others->Fill1d("pullx"+suffix, r.x()/er.x());
1123	>	h_others->Fill1d("pully"+suffix, r.y()/er.y());
1124	>	h_others->Fill1d("pullz"+suffix, r.z()/er.z());
1125	>	h_others->Fill1d("errPVx"+suffix, er.x());
1126	>	h_others->Fill1d("errPVy"+suffix, er.y());
1127	>	h_others->Fill1d("errPVz"+suffix, er.z());
1128		}
1129
1130	<	PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk) {
1130	>	PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1131		map<int,double> data;
1132		data.clear();
1133	<	double fpar[2];
1133	>	double fpar[2] = {-999,-999};
1134		double cm2um = 10000;
1135		stringstream ss;
1136		ss << ntk;
1137	<	string suffix = ss.str();
1138	<
1137	>	string suffix = ss.str();
1138	>	if(datamode == 0 )  suffix = "_" + suffix;
1139	>	if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1140	>	if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1141	>	if(datamode == 3 )  suffix = "_mct_" + suffix;
1142	>
1143		if(analyze_) {
1144		for ( int i = 0; i < 6; ++i) {
611	–	cout << "Here" << endl;
1145		switch (i) {
1146		case 0:
1147	<	fit(h["resx_"+suffix], fpar);
1147	>	if (!h_others->ReadHisto1D("deltax"+suffix)) break;
1148	>	fit(h_others->ReadHisto1D("deltax"+suffix), fpar);
1149		data[i] = fpar[0]*cm2um;
1150		break;
1151		case 1:
1152	<	fit(h["resy_"+suffix], fpar);
1152	>	if (!h_others->ReadHisto1D("deltay"+suffix)) break;
1153	>	fit(h_others->ReadHisto1D("deltay"+suffix), fpar);
1154		data[i] = fpar[0]*cm2um;
1155		break;
1156		case 2:
1157	<	fit(h["resz_"+suffix], fpar);
1157	>	if (!h_others->ReadHisto1D("deltaz"+suffix)) break;
1158	>	fit(h_others->ReadHisto1D("deltaz"+suffix), fpar);
1159		data[i] = fpar[0]*cm2um;
1160		break;
1161		case 3:
1162	<	fit(h["pullx_"+suffix], fpar);
1162	>	if (!h_others->ReadHisto1D("pullx"+suffix)) break;
1163	>	fit(h_others->ReadHisto1D("pullx"+suffix), fpar);
1164		data[i] = fpar[0];
1165		break;
1166		case 4:
1167	<	fit(h["pully_"+suffix], fpar);
1167	>	if (!h_others->ReadHisto1D("pully"+suffix)) break;
1168	>	fit(h_others->ReadHisto1D("pully"+suffix), fpar);
1169		data[i] = fpar[0];
1170		break;
1171		case 5:
1172	<	fit(h["pullz_"+suffix], fpar);
1172	>	if (!h_others->ReadHisto1D("pullz"+suffix)) break;
1173	>	fit(h_others->ReadHisto1D("pullz"+suffix), fpar);
1174		data[i] = fpar[0];
1175		break;
1176		}
1177	<	if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1177	>	if (data[i] > 0) edm::LogInfo("Analysis") << "[Analysis] data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1178		}
1179		}
1180		else
#	Line 650 | Line 1189 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1189		ntk);
1190		}
1191
1192	<	void PVStudy::fit(TH1 *&hdist, double data[]){
1192	>	void PVStudy::fit(TH1 *hdist, double fitPars[]){
1193		TAxis *axis0 = hdist->GetXaxis();
1194		int nbin = axis0->GetLast();
1195		double nOF = hdist->GetBinContent(nbin+1);
1196		double nUF = hdist->GetBinContent(0);
1197	<	if ( verbose_ ){
1198	<	cout << "Last bin = " << nbin;
1199	<	cout << "; hdist: Overflow Entries = " << nOF;
1200	<	cout << "; Underflow Entries = " << nUF;
1201	<	cout << "; hdist->GetEntries() = " << hdist->GetEntries() << endl;
1202	<	}
1203	<	if (hdist->GetEntries() - nOF - nUF >= 5) { // FIXME: for reasonable Gaussian fit
1204	<	hdist->Fit("gaus","Q");
1205	<	TF1 *fgaus = hdist->GetFunction("gaus");
1206	<	if (verbose_) cout << "got function" << endl;
1207	<	data[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1197	>	//   double fitRange = axis0->GetBinUpEdge(nbin);
1198	>	double fitRange = axis0->GetXmax();
1199	>	double sigMax[2] = {0.,0.};
1200	>
1201	>	edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
1202	>	<< "; hdist: Overflow Entries = " << nOF
1203	>	<< "; Underflow Entries = " << nUF
1204	>	<< "; hdist->GetEntries() = " << hdist->GetEntries()
1205	>	<< "; fitting range = " << -fitRange << " to " << fitRange << endl;
1206	>
1207	>	if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit
1208	>	for (int bi = 0; bi < nbin; bi++) {
1209	>	if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) {
1210	>	sigMax[0] = axis0->GetBinLowEdge(bi);
1211	>	if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1212	>	}
1213	>	if ( (axis0->GetBinLowEdge(bi) >= 0) && (hdist->GetBinContent(bi) > 0) ) {
1214	>	sigMax[0] = axis0->GetBinUpEdge(bi);
1215	>	if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1216	>	}
1217	>	}
1218	>	//edm::LogInfo("Analysis") << "[Fit] Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1219	>
1220	>	TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1221	>	fgaus->SetParameter(1, 0.);
1222	>	fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
1223	>	fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
1224	>	fgaus->SetLineColor(4);
1225	>	hdist->Fit(fgaus,"Q");
1226	>	gPad->Update();
1227	>	TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
1228	>	s->SetOptStat(1111111);
1229	>	s->SetOptFit(0111);
1230	>	gPad->Update();
1231	>	//edm::LogInfo("Analysis") << "[Fit] got function" << endl;
1232	>	fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1233	>	}
1234	>	else
1235	>	fitPars[0] = -999;
1236	>	}
1237	>
1238	>	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype) {
1239	>	string suffix;
1240	>	suffix = ""; // for unsplit tracks
1241	>	if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1242	>	if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1243	>
1244	>	h_pvtrk->Fill1d("nTrk"+suffix, trackColl->size());
1245	>	for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1246	>	itTrack != trackColl->end();
1247	>	++itTrack) {
1248	>	h_pvtrk->Fill1d("trkPt"+suffix, itTrack->pt());
1249	>	h_pvtrk->Fill1d("trkDxy"+suffix, itTrack->dxy());
1250	>	h_pvtrk->Fill1d("trkDz"+suffix, itTrack->dz());
1251	>	h_pvtrk->Fill1d("trkEta"+suffix, itTrack->eta());
1252	>	h_pvtrk->Fill1d("trkPhi"+suffix, itTrack->phi());
1253	>	}
1254	>	}
1255	>
1256	>	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
1257	>	string suffix;
1258	>	suffix = ""; // for unsplit tracks
1259	>	if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1260	>	if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1261	>	int ivtx = 0;
1262	>	for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1263	>	v!=vertexColl->end(); ++v, ++ivtx) {
1264	>	if(fillNtuple) {
1265	>	if(datatype == 0) {
1266	>	nTrkPV_[ivtx] = v->tracksSize();
1267	>	sumptsq_[ivtx] = sumPtSquared(*v);
1268	>	isValid_[ivtx] = v->isValid();
1269	>	isFake_[ivtx] = v->isFake();
1270	>	recx_[ivtx] = v->x();
1271	>	recy_[ivtx] = v->y();
1272	>	recz_[ivtx] = v->z();
1273	>	recx_err_[ivtx] = v->xError();
1274	>	recy_err_[ivtx] = v->yError();
1275	>	recz_err_[ivtx] = v->zError();
1276	>	}
1277	>	if(datatype == 1) {
1278	>	nTrkPV1_spl_[ivtx] = v->tracksSize();
1279	>	sumptsq1_spl_[ivtx] = sumPtSquared(*v);
1280	>	isValid1_spl_[ivtx] = v->isValid();
1281	>	isFake1_spl_[ivtx] = v->isFake();
1282	>	recx1_spl_[ivtx] = v->x();
1283	>	recy1_spl_[ivtx] = v->y();
1284	>	recz1_spl_[ivtx] = v->z();
1285	>	recx1_err_spl_[ivtx] = v->xError();
1286	>	recy1_err_spl_[ivtx] = v->yError();
1287	>	recz1_err_spl_[ivtx] = v->zError();
1288	>
1289	>	}
1290	>	if(datatype == 2) {
1291	>	nTrkPV2_spl_[ivtx] = v->tracksSize();
1292	>	sumptsq2_spl_[ivtx] = sumPtSquared(*v);
1293	>	isValid2_spl_[ivtx] = v->isValid();
1294	>	isFake2_spl_[ivtx] = v->isFake();
1295	>	recx2_spl_[ivtx] = v->x();
1296	>	recy2_spl_[ivtx] = v->y();
1297	>	recz2_spl_[ivtx] = v->z();
1298	>	recx2_err_spl_[ivtx] = v->xError();
1299	>	recy2_err_spl_[ivtx] = v->yError();
1300	>	recz2_err_spl_[ivtx] = v->zError();
1301	>	}
1302	>	if(datatype == 4) { // for pixelVertices
1303	>	nTrkPV_pxlpvtx_[ivtx] = v->tracksSize();
1304	>	sumptsq_pxlpvtx_[ivtx] = sumPtSquared(*v);
1305	>	isValid_pxlpvtx_[ivtx] = v->isValid();
1306	>	isFake_pxlpvtx_[ivtx] = v->isFake();
1307	>	recx_pxlpvtx_[ivtx] = v->x();
1308	>	recy_pxlpvtx_[ivtx] = v->y();
1309	>	recz_pxlpvtx_[ivtx] = v->z();
1310	>	recx_err_pxlpvtx_[ivtx] = v->xError();
1311	>	recy_err_pxlpvtx_[ivtx] = v->yError();
1312	>	recz_err_pxlpvtx_[ivtx] = v->zError();
1313	>	}
1314	>	}
1315	>	}
1316	>	// For histogram only fill the leading pvtx parameters
1317	>	if(fillHisto) {
1318	>	h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1319	>	try {
1320	>	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1321	>	t!=vertexColl->begin()->tracks_end(); t++) {
1322	>	// illegal charge
1323	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1324	>	//      h_pvtrk->Fill1d("trkPtPV", 0.);
1325	>	}
1326	>	else {
1327	>	h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1328	>	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1329	>	h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1330	>	h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1331	>	h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1332	>	}
1333	>	}
1334	>	}
1335	>	catch (...) {
1336	>	// exception thrown when trying to use linked track
1337	>	//          h_pvtrk->Fill1d("trkPtPV", 0.);
1338	>	}
1339	>	}
1340	>
1341	>	}
1342	>
1343	>	void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
1344	>	{
1345	>	std::string suffix;
1346	>	// Set the vertexColl and histogram suffix according to datamode
1347	>	if (datatype == 1) {
1348	>	suffix = "_spl1_mct";
1349	>	nrecPV_spl1_mct_ = 0;
1350	>	}
1351	>	if (datatype == 2) {
1352	>	suffix = "_spl2_mct";
1353	>	nrecPV_spl2_mct_ = 0;
1354	>	}
1355	>	if (datatype == 3) {
1356	>	suffix = "_mct";
1357	>	nrecPV_mct_ = 0;
1358	>	}
1359	>
1360	>	//========================================================
1361	>	//  look for a matching reconstructed vertex in vertexColl
1362	>	//========================================================
1363	>
1364	>	for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1365	>	vrec!=vtxColl->end(); ++vrec){
1366	>	int nrectrk = vrec->tracksSize();
1367	>	vsim->recVtx=NULL;
1368	>
1369	>	edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1370	>	edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1371	>
1372	>	if ( matchVertex(*vsim,*vrec) && vrec->isValid() && !vrec->isFake() ) {
1373	>	vsim->recVtx=&(*vrec);
1374	>
1375	>	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1376	>	//if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
1377	>	//|| (!vsim->recVtx) )
1378	>	//vsim->recVtx=&(*vrec);
1379	>	//}
1380	>	edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1381	>
1382	>	double fres_mct[3];
1383	>	double ferror_mct[3];
1384	>
1385	>	fres_mct[0] = vsim->recVtx->x()-vsim->x;
1386	>	fres_mct[1] = vsim->recVtx->y()-vsim->y;
1387	>	fres_mct[2] = vsim->recVtx->z()-vsim->z;
1388	>	ferror_mct[0] = vsim->recVtx->xError();
1389	>	ferror_mct[1] = vsim->recVtx->yError();
1390	>	ferror_mct[2] = vsim->recVtx->zError();
1391	>
1392	>	h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1393	>	h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1394	>	h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1395	>	h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1396	>	h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1397	>	h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1398	>	h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1399	>	h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1400	>	h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1401	>	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1402	>	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1403	>	nrectrk));
1404	>	// Fill histo according to its track multiplicity
1405	>	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1406	>	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1407	>	nrectrk, datatype);
1408	>
1409	>	if(saventuple_) {
1410	>	//Fill the values for variables in ftree_
1411	>	if(datatype == 1) {
1412	>	nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1413	>	deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1414	>	deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1415	>	deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1416	>	pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
1417	>	pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1418	>	pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1419	>	nrecPV_spl1_mct_++;
1420	>	}
1421	>
1422	>	if(datatype == 2) {
1423	>	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1424	>	deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1425	>	deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1426	>	deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1427	>	pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
1428	>	pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1429	>	pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1430	>	nrecPV_spl2_mct_++;
1431	>	}
1432	>	if(datatype == 3) {
1433	>	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1434	>	deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1435	>	deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1436	>	deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1437	>	pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
1438	>	pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1439	>	pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1440	>	nrecPV_mct_++;
1441	>	}
1442	>	} // End of  if(saventuple_) {
1443	>	} //  if ( matchVertex(*vsim,*vrec) ) {
1444	>	else {  // no rec vertex found for this simvertex
1445	>	edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1446	>	}
1447	>	}
1448	>	}
1449	>
1450	>
1451	>	void PVStudy::SetVarToZero() {
1452	>	//Greg's variables
1453	>	fntrk_ = 0;
1454	>	//pvtx position (x,y,z) residual and error
1455	>	for(int i = 0; i<3;i++)    {
1456	>	fres_[i] = 0;
1457	>	ferror_[i] = 0;
1458	>	}
1459	>
1460	>	//Yanyan's variables
1461	>	// Number of reconstructed vertices
1462	>	nrecPV_ = 0;
1463	>	nrecPV1_spl_ = 0;
1464	>	nrecPV2_spl_ = 0;
1465	>	nrecPV_twovtx_ = 0;
1466	>	nsimPV_ = 0;
1467	>
1468	>	nrecPV_mct_ = 0;
1469	>	nrecPV_spl1_mct_ = 0;
1470	>	nrecPV_spl2_mct_ = 0;
1471	>
1472	>	// Mininum separation between the secondary pvtxes and leading pvtx
1473	>	min_zsep_ = 9999.0;
1474	>	min_ntrksep_ = 9999.0;
1475	>	min_sumpt2sep_ = -9999.0;
1476	>
1477	>
1478	>	for (int i = 0; i < nMaxPVs_; i++) {
1479	>	// recoVertices with all tracks
1480	>	nTrkPV_[i] = 0; // Number of tracks in the pvtx
1481	>	sumptsq_[i] = 0;
1482	>	isValid_[i] = -1;
1483	>	isFake_[i] = -1;
1484	>	recx_[i] = 0;
1485	>	recy_[i] = 0;
1486	>	recz_[i] = 0;
1487	>	recx_err_[i] = 0;
1488	>	recy_err_[i] = 0;
1489	>	recz_err_[i] = 0;
1490	>
1491	>	// recoVertices with splitTrack1
1492	>	nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx
1493	>	sumptsq1_spl_[i] = 0;
1494	>	isValid1_spl_[i] = -1;
1495	>	isFake1_spl_[i] = -1;
1496	>	recx1_spl_[i] = 0;
1497	>	recy1_spl_[i] = 0;
1498	>	recz1_spl_[i] = 0;
1499	>	recx1_err_spl_[i] = 0;
1500	>	recy1_err_spl_[i] = 0;
1501	>	recz1_err_spl_[i] = 0;
1502	>
1503	>	// recoVertices with splitTrack2
1504	>	nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx
1505	>	sumptsq2_spl_[i] = 0;
1506	>	isValid2_spl_[i] = -1;
1507	>	isFake2_spl_[i] = -1;
1508	>	recx2_spl_[i] = 0;
1509	>	recy2_spl_[i] = 0;
1510	>	recz2_spl_[i] = 0;
1511	>	recx2_err_spl_[i] = 0;
1512	>	recy2_err_spl_[i] = 0;
1513	>	recz2_err_spl_[i] = 0;
1514	>
1515	>	//pixelVertices
1516	>	nTrkPV_pxlpvtx_[i] = 0; // Number of tracks in the pvtx
1517	>	sumptsq_pxlpvtx_[i] = 0;
1518	>	isValid_pxlpvtx_[i] = -1;
1519	>	isFake_pxlpvtx_[i] = -1;
1520	>	recx_pxlpvtx_[i] = 0;
1521	>	recy_pxlpvtx_[i] = 0;
1522	>	recz_pxlpvtx_[i] = 0;
1523	>	recx_err_pxlpvtx_[i] = 0;
1524	>	recy_err_pxlpvtx_[i] = 0;
1525	>	recz_err_pxlpvtx_[i] = 0;
1526	>
1527	>	// matched two-vertices
1528	>	nTrkPV1_spl_twovtx_[i] = 0;
1529	>	nTrkPV2_spl_twovtx_[i] = 0;
1530	>	nTrkPV_twovtx_[i] = 0;
1531	>	sumptsq1_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1532	>	sumptsq2_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1533	>	deltax_twovtx_[i] = 0;
1534	>	deltay_twovtx_[i] = 0;
1535	>	deltaz_twovtx_[i] = 0;
1536	>	errx_twovtx_[i] = 0;
1537	>	erry_twovtx_[i] = 0;
1538	>	errz_twovtx_[i] = 0;
1539	>	pullx_twovtx_[i] = 0;
1540	>	pully_twovtx_[i] = 0;
1541	>	pullz_twovtx_[i] = 0;
1542	>
1543	>	//simpv
1544	>	nsimTrkPV_[i] = 0;
1545	>	simx_[i] = 0;
1546	>	simy_[i] = 0;
1547	>	simz_[i] = 0;
1548	>	simptsq_[i] = 0;
1549	>
1550	>	// residual and pulls with mc truth required
1551	>	deltax_mct_[i] = 0;
1552	>	deltay_mct_[i] = 0;
1553	>	deltaz_mct_[i] = 0;
1554	>	pullx_mct_[i] = 0;
1555	>	pully_mct_[i] = 0;
1556	>	pullz_mct_[i] = 0;
1557	>
1558	>	deltax_spl1_mct_[i] = 0;
1559	>	deltay_spl1_mct_[i] = 0;
1560	>	deltaz_spl1_mct_[i] = 0;
1561	>	pullx_spl1_mct_[i] = 0;
1562	>	pully_spl1_mct_[i] = 0;
1563	>	pullz_spl1_mct_[i] = 0;
1564	>
1565	>	deltax_spl2_mct_[i] = 0;
1566	>	deltay_spl2_mct_[i] = 0;
1567	>	deltaz_spl2_mct_[i] = 0;
1568	>	pullx_spl2_mct_[i] = 0;
1569	>	pully_spl2_mct_[i] = 0;
1570	>	pullz_spl2_mct_[i] = 0;
1571	>	}
1572	>
1573	>	}
1574	>
1575	>	double PVStudy::sumPtSquared(const reco::Vertex & v)  {
1576	>	double sum = 0.;
1577	>	double pT;
1578	>	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1579	>	pT = (**it).pt();
1580	>	sum += pT*pT;
1581		}
1582	+	return sum;
1583		}
1584
672	–	//define this as a plug-in
673	–	DEFINE_FWK_MODULE(PVStudy);

Diff Legend

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