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

Diff Legend

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