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.9 by yygao, Thu Oct 29 02:05:54 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	+	//=======================================================================
67	+	// Get configuration for TrackTupleMaker
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);
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	>	verbose_         = iConfig.getUntrackedParameter<bool>("verbose",false);
77		realData_        = iConfig.getUntrackedParameter<bool>("realData",false);
78		analyze_         = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
79	<	outputfilename_  = iConfig.getUntrackedParameter<string>("OutputFileName");
79	>	saventuple_         = iConfig.getUntrackedParameter<bool>("saventuple",false);
80	>	outputfilename_  = iConfig.getUntrackedParameter<string>("OutputFileName");
81	>	ntrkdiffcut_         = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
82	>	zsigncut_         = iConfig.getUntrackedParameter<int>("zsigncut");
83	>	nTrkMin_ =  iConfig.getUntrackedParameter<int>("nTrkMin");
84	>	nTrkMax_ =  iConfig.getUntrackedParameter<int>("nTrkMax");
85
86		//now do what ever initialization is needed
87		pvinfo.clear();
88
75	–	file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
76	–	ftree_ = new TTree("mytree","mytree");
77	–	ftree_->AutoSave();
78	–	ftree_->Branch("residual",&fres,"fres[3]/D");
79	–	ftree_->Branch("error",&ferror,"ferror[3]/D");
80	–	ftree_->Branch("nTrk",&fntrk,"fntrk/I");
89
90	<	edm::Service<TFileService> fs;
91	<	//   TFileDirectory subDir = fs->mkdir( "Summary" );
92	<	//   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
90	>	// Specify the data mode vector
91	>	if(realData_) datamode.push_back(0);
92	>	else {
93	>	datamode.push_back(0);
94	>	datamode.push_back(1);
95	>	datamode.push_back(2);
96	>	datamode.push_back(3);
97	>	}
98
99	<	h_nTrk           = fs->make<TH1F>("nTrk","Num of total tracks",300,0,300);
100	<	h_nTrkPV         = fs->make<TH1F>("nTrkPV","Num of Tracks from PV",300,0,300);
101	<	h_trkPt          = fs->make<TH1D>("trkPt","Pt of all tracks",100,0,100);
102	<	h_trkPtPV        = fs->make<TH1D>("trkPtPV","Pt dist of tracks from PV",100,0,100);
99	>	// Create ntuple files if needed
100	>	if(saventuple_) {
101	>	file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
102	>	ftree_ = new TTree("mytree","mytree");
103	>	ftree_->AutoSave();
104	>	ftree_->Branch("residual",&fres_,"fres_[3]/D");
105	>	ftree_->Branch("error",&ferror_,"ferror_[3]/D");
106	>	ftree_->Branch("nTrkPV",&fntrk_,"fntrk_/I");
107	>
108	>	// pvtxtree_ analyzing the pvtxs ootb
109	>	pvtxtree_ = new TTree("pvtxtree","pvtxtree");
110	>	//pvtx using all tracks
111	>	pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
112	>	pvtxtree_->Branch("nTrkPV",&nTrkPV_,"nTrkPV[nrecPV]/I");
113	>	pvtxtree_->Branch("recx",&recx_,"recx[nrecPV]/D");
114	>	pvtxtree_->Branch("recy",&recy_,"recy[nrecPV]/D");
115	>	pvtxtree_->Branch("recz",&recz_,"recz[nrecPV]/D");
116	>	pvtxtree_->Branch("recx_err",&recx_err_,"recx_err[nrecPV]/D");
117	>	pvtxtree_->Branch("recy_err",&recy_err_,"recy_err[nrecPV]/D");
118	>	pvtxtree_->Branch("recz_err",&recz_err_,"recz_err[nrecPV]/D");
119	>	//pvtx using splittracks1
120	>	pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
121	>	pvtxtree_->Branch("nTrkPV1_spl",&nTrkPV1_spl_,"nTrkPV1_spl[nrecPV1_spl]/I");
122	>	pvtxtree_->Branch("recx1_spl",&recx1_spl_,"recx1_spl[nrecPV1_spl]/D");
123	>	pvtxtree_->Branch("recy1_spl",&recy1_spl_,"recy1_spl[nrecPV1_spl]/D");
124	>	pvtxtree_->Branch("recz1_spl",&recz1_spl_,"recz1_spl[nrecPV1_spl]/D");
125	>	pvtxtree_->Branch("recx1_err_spl",&recx1_err_spl_,"recx1_err_spl[nrecPV1_spl]/D");
126	>	pvtxtree_->Branch("recy1_err_spl",&recy1_err_spl_,"recy1_err_spl[nrecPV1_spl]/D");
127	>	pvtxtree_->Branch("recz1_err_spl",&recz1_err_spl_,"recz1_err_spl[nrecPV1_spl]/D");
128	>	//pvtx using splittracks2
129	>	pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
130	>	pvtxtree_->Branch("nTrkPV2_spl",&nTrkPV2_spl_,"nTrkPV2_spl[nrecPV2_spl]/I");
131	>	pvtxtree_->Branch("recx2_spl",&recx2_spl_,"recx2_spl[nrecPV2_spl]/D");
132	>	pvtxtree_->Branch("recy2_spl",&recy2_spl_,"recy2_spl[nrecPV2_spl]/D");
133	>	pvtxtree_->Branch("recz2_spl",&recz2_spl_,"recz2_spl[nrecPV2_spl]/D");
134	>	pvtxtree_->Branch("recx2_err_spl",&recx2_err_spl_,"recx2_err_spl[nrecPV2_spl]/D");
135	>	pvtxtree_->Branch("recy2_err_spl",&recy2_err_spl_,"recy2_err_spl[nrecPV2_spl]/D");
136	>	pvtxtree_->Branch("recz2_err_spl",&recz2_err_spl_,"recz2_err_spl[nrecPV2_spl]/D");
137	>	//Fill the variables in the twovtx pair (recvtx1, recvtx2)
138	>	pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
139	>	pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
140	>	pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
141	>	pvtxtree_->Branch("sumptsq1_spl_twovtx",&sumptsq1_spl_twovtx_,"sumptsq1_spl_twovtx[nrecPV_twovtx]/D");
142	>	pvtxtree_->Branch("sumptsq2_spl_twovtx",&sumptsq2_spl_twovtx_,"sumptsq2_spl_twovtx[nrecPV_twovtx]/D");
143	>	pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
144	>	pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
145	>	pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
146	>	pvtxtree_->Branch("deltaz_twovtx",&deltaz_twovtx_,"deltaz_twovtx[nrecPV_twovtx]/D");
147	>	pvtxtree_->Branch("errx_twovtx",&errx_twovtx_,"errx_twovtx[nrecPV_twovtx]/D");
148	>	pvtxtree_->Branch("erry_twovtx",&erry_twovtx_,"erry_twovtx[nrecPV_twovtx]/D");
149	>	pvtxtree_->Branch("errz_twovtx",&errz_twovtx_,"errz_twovtx[nrecPV_twovtx]/D");
150	>	pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
151	>	pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
152	>	pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
153	>
154	>	// MC variables
155	>	if(!realData_) {
156	>	pvtxtree_->Branch("nsimPV",&nsimPV_,"nsimPV/I");
157	>	pvtxtree_->Branch("nsimTrkPV",&nsimTrkPV_,"nsimTrkPV[nsimPV]/I");
158	>	pvtxtree_->Branch("simx",&simx_,"simx[nsimPV]/D");
159	>	pvtxtree_->Branch("simy",&simy_,"simy[nsimPV]/D");
160	>	pvtxtree_->Branch("simz",&simz_,"simz[nsimPV]/D");
161	>	pvtxtree_->Branch("simptsq",&simptsq_,"simptsq[nsimPV]/D");
162	>
163	>	// For pvtxs with all the tracks
164	>	pvtxtree_->Branch("nrecPV_mct",&nrecPV_mct_,"nrecPV_mct/I");
165	>	pvtxtree_->Branch("deltax_mct",&deltax_mct_,"deltax_mct[nrecPV_mct]/D");
166	>	pvtxtree_->Branch("deltay_mct",&deltay_mct_,"deltay_mct[nrecPV_mct]/D");
167	>	pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
168	>	pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
169	>	pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
170	>	pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
171	>	// For pvtxs with splittracks1
172	>	pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
173	>	pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
174	>	pvtxtree_->Branch("deltay_spl1_mct",&deltay_spl1_mct_,"deltay_spl1_mct[nrecPV_spl1_mct]/D");
175	>	pvtxtree_->Branch("deltaz_spl1_mct",&deltaz_spl1_mct_,"deltaz_spl1_mct[nrecPV_spl1_mct]/D");
176	>	pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
177	>	pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
178	>	pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
179	>	// For pvtxs with splittracks1
180	>	pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
181	>	pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
182	>	pvtxtree_->Branch("deltay_spl2_mct",&deltay_spl2_mct_,"deltay_spl2_mct[nrecPV_spl2_mct]/D");
183	>	pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
184	>	pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
185	>	pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
186	>	pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
187	>	}
188	>	}
189
190	<	if (!realData_) {
191	<	h_deltax       = fs->make<TH1D>("deltax","x-distance (Rec - Sim)",800,-0.04,0.04);
192	<	h_deltay       = fs->make<TH1D>("deltay","y-distance (Rec - Sim)",800,-0.04,0.04);
193	<	h_deltaz       = fs->make<TH1D>("deltaz","z-distance (Rec - Sim)",800,-0.04,0.04);
194	<	h_pullx        = fs->make<TH1D>("pullx","x-pull (Rec - Sim)",1000,-25.,25.);
195	<	h_pully        = fs->make<TH1D>("pully","y-pull (Rec - Sim)",1000,-25.,25.);
196	<	h_pullz        = fs->make<TH1D>("pullz","z-pull (Rec - Sim)",1000,-25.,25.);
197	<
198	<	// Summary of analysis:
199	<	if (analyze_) {
200	<	h_resx_nTrk    = fs->make<TH2D>("resx_ntrk","x-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
201	<	h_resx_nTrk->SetMarkerStyle(21);
202	<	h_resx_nTrk->SetMarkerColor(4);
203	<	h_resx_nTrk->GetXaxis()->SetTitle("Num of tracks");
204	<	h_resx_nTrk->GetYaxis()->SetTitle("#mum");
205	<	h_resy_nTrk    = fs->make<TH2D>("resy_ntrk","y-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
206	<	h_resy_nTrk->SetMarkerStyle(21);
207	<	h_resy_nTrk->SetMarkerColor(4);
208	<	h_resy_nTrk->GetXaxis()->SetTitle("Num of tracks");
209	<	h_resy_nTrk->GetYaxis()->SetTitle("#mum");
210	<	h_resz_nTrk    = fs->make<TH2D>("resz_ntrk","z-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
211	<	h_resz_nTrk->SetMarkerStyle(21);
212	<	h_resz_nTrk->SetMarkerColor(4);
213	<	h_resz_nTrk->GetXaxis()->SetTitle("Num of tracks");
214	<	h_resz_nTrk->GetYaxis()->SetTitle("#mum");
215	<	h_pullx_nTrk   = fs->make<TH2D>("pullx_ntrk","x-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
216	<	h_pullx_nTrk->SetMarkerStyle(21);
217	<	h_pullx_nTrk->SetMarkerColor(4);
218	<	h_pullx_nTrk->SetBit(TH1::kCanRebin);
219	<	h_pullx_nTrk->GetXaxis()->SetTitle("Num of tracks");
220	<	h_pully_nTrk   = fs->make<TH2D>("pully_ntrk","y-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
221	<	h_pully_nTrk->SetMarkerStyle(21);
222	<	h_pully_nTrk->SetMarkerColor(4);
223	<	h_pully_nTrk->SetBit(TH1::kCanRebin);
224	<	h_pully_nTrk->GetXaxis()->SetTitle("Num of tracks");
225	<	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");
190	>	edm::Service<TFileService> fs;
191	>	TFileDirectory subDir = fs->mkdir( "Summary" );
192	>	TFileDirectory subDir1 = fs->mkdir( "Others" );
193	>	//   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
194	>
195	>	setRootStyle();
196	>
197	>	//Book histograms for track and pvtx parameters for each splitted pvtx
198	>	//i=0 : x (as in unsplit track collection)
199	>	//i=1 : x1_spl_
200	>	//i=2 : x2_spl_
201	>
202	>	for(int i=0;i<3;i++)
203	>	{
204	>	string suffix;
205	>	if(i == 0) suffix = "";
206	>	else {
207	>	stringstream ss;
208	>	ss << i;
209	>	suffix =ss.str()+"_spl";
210	>	}
211	>	h["nTrk"+suffix]   = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
212	>	h["trkPt"+suffix]  = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
213	>	h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#eta of rec tracks "+suffix),100,-3,3);
214	>	h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
215	>	h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-5,5);
216	>	h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
217	>
218	>	h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
219	>	h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
220	>	h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#eta of rec tracks in PV"+suffix),100,-3,3);
221	>	h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
222	>	h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);
223	>	h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
224	>	h["nrecPV"+suffix]   = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);
225	>	suffix.clear();
226		}
227	+	h["nrecPVDiff"]     = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
228	+	h["nTrkPVDiff"]     = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
229	+	h["nTrkPVRelDiff"]  = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
230	+
231	+	//Book histograms sensitive to data/mc
232	+	for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
233	+	string suffix;
234	+	if(verbose_)  cout<<"datamode = "<< *it<<endl;
235	+	switch(*it) {
236	+	case 0: suffix = "";
237	+	break;
238	+	case 1: suffix = "_spl1_mct";
239	+	break;
240	+	case 2: suffix = "_spl2_mct";
241	+	break;
242	+	case 3: suffix = "_mct";
243	+	break;
244	+	}
245	+	h["deltax"+suffix] = subDir.make<TH1D>(TString("deltax"+suffix), TString("x-residual pvtx"+suffix),800,-0.04,0.04);
246	+	h["deltay"+suffix] = subDir.make<TH1D>(TString("deltay"+suffix), TString("y-residual pvtx"+suffix),800,-0.04,0.04);
247	+	h["deltaz"+suffix] = subDir.make<TH1D>(TString("deltaz"+suffix), TString("z-residual pvtx"+suffix),800,-0.04,0.04);
248	+	h["pullx"+suffix]  = subDir.make<TH1D>(TString("pullx"+suffix), TString("x-pull pvtx"+suffix),800,-10,10);
249	+	h["pully"+suffix]  = subDir.make<TH1D>(TString("pully"+suffix), TString("y-pull pvtx"+suffix),800,-10,10);
250	+	h["pullz"+suffix]  = subDir.make<TH1D>(TString("pullz"+suffix), TString("z-pull pvtx"+suffix),800,-10,10);
251	+	h["errPVx"+suffix] = subDir.make<TH1D>(TString("errPVx"+suffix), TString("X"+suffix+" vertex error"),100,0.,0.02);
252	+	h["errPVy"+suffix] = subDir.make<TH1D>(TString("errPVy"+suffix), TString("Y"+suffix+" vertex error"),100,0.,0.02);
253	+	h["errPVz"+suffix] = subDir.make<TH1D>(TString("errPVz"+suffix), TString("Z"+suffix+" vertex error"),100,0.,0.02);
254	+
255	+	// Book residual, error and pull histograms for each nTrk bin
256		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
257		stringstream ss;
258		ss << ntrk;
259	<	string suffix = ss.str();
260	<	h["resx_"+suffix]  = fs->make<TH1D> (TString("deltax_"+suffix),"x-distance (Rec - Sim)",100,-0.02,0.02);
261	<	h["resx_"+suffix]->GetXaxis()->SetTitle("cm");
262	<	h["resy_"+suffix]  = fs->make<TH1D> (TString("deltay_"+suffix),"y-distance (Rec - Sim)",100,-0.02,0.02);
263	<	h["resy_"+suffix]->GetXaxis()->SetTitle("cm");
264	<	h["resz_"+suffix]  = fs->make<TH1D> (TString("deltaz_"+suffix),"z-distance (Rec - Sim)",100,-0.02,0.02);
265	<	h["resz_"+suffix]->GetXaxis()->SetTitle("cm");
266	<	h["pullx_"+suffix] = fs->make<TH1D> (TString("pullx_"+suffix),"x-pull (Rec - Sim)",100,-10.,10.);
267	<	h["pully_"+suffix] = fs->make<TH1D> (TString("pully_"+suffix),"y-pull (Rec - Sim)",100,-10.,10.);
268	<	h["pullz_"+suffix] = fs->make<TH1D> (TString("pullz_"+suffix),"z-pull (Rec - Sim)",100,-10.,10.);
269	<	suffix.clear();
270	<	}
259	>	string suffix2 = suffix+"_"+ss.str();
260	>	h["deltax"+suffix2]  = subDir1.make<TH1D>(TString("deltax"+suffix2), TString("x-residual of pvtx"+suffix),100,-0.02,0.02);
261	>	h["deltax"+suffix2]->GetXaxis()->SetTitle("cm");
262	>	h["deltay"+suffix2]  = subDir1.make<TH1D>(TString("deltay"+suffix2), TString("y-residual of pvtx"+suffix),100,-0.02,0.02);
263	>	h["deltay"+suffix2]->GetXaxis()->SetTitle("cm");
264	>	h["deltaz"+suffix2]  = subDir1.make<TH1D>(TString("deltaz"+suffix2), TString("z-residual of pvtx"+suffix),100,-0.02,0.02);
265	>	h["deltaz"+suffix2]->GetXaxis()->SetTitle("cm");
266	>	h["pullx"+suffix2] = subDir1.make<TH1D>(TString("pullx"+suffix2), TString("x-pull of pvtx"+suffix),100,-10.,10.);
267	>	h["pully"+suffix2] = subDir1.make<TH1D>(TString("pully"+suffix2), TString("y-pull of pvtx"+suffix),100,-10.,10.);
268	>	h["pullz"+suffix2] = subDir1.make<TH1D>(TString("pullz"+suffix2), TString("z-pull of pvtx"+suffix),100,-10.,10.);
269	>	h["errPVx"+suffix2] = subDir1.make<TH1D>(TString("errPVx"+suffix2), TString("X"+suffix+" vertex error"),100,0.,0.02);
270	>	h["errPVy"+suffix2] = subDir1.make<TH1D>(TString("errPVy"+suffix2), TString("Y"+suffix+" vertex error"),100,0.,0.02);
271	>	h["errPVz"+suffix2] = subDir1.make<TH1D>(TString("errPVz"+suffix2), TString("Z"+suffix+" vertex error"),100,0.,0.02);
272	>	suffix2.clear();
273	>	} // End of   for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
274	>
275	>	// Book residual and pull histograms only when analyzeOntheFly is enabled
276	>	if(analyze_) {
277	>	h2["resx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resx"+suffix+"_nTrk"), TString("x-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
278	>	h2["resx"+suffix+"_nTrk"]->SetMarkerStyle(21);
279	>	h2["resx"+suffix+"_nTrk"]->SetMarkerColor(4);
280	>	h2["resx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
281	>	h2["resx"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
282	>	h2["resy"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resy"+suffix+"_nTrk"), TString("y-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
283	>	h2["resy"+suffix+"_nTrk"]->SetMarkerStyle(21);
284	>	h2["resy"+suffix+"_nTrk"]->SetMarkerColor(4);
285	>	h2["resy"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
286	>	h2["resy"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
287	>	h2["resz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resz"+suffix+"_nTrk"), TString("z-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
288	>	h2["resz"+suffix+"_nTrk"]->SetMarkerStyle(21);
289	>	h2["resz"+suffix+"_nTrk"]->SetMarkerColor(4);
290	>	h2["resz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
291	>	h2["resz"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
292	>	h2["pullx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullx"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
293	>	h2["pullx"+suffix+"_nTrk"]->SetMarkerStyle(21);
294	>	h2["pullx"+suffix+"_nTrk"]->SetMarkerColor(4);
295	>	h2["pullx"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
296	>	h2["pullx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
297	>	h2["pully"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pully"+suffix+"_nTrk"), TString("y-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
298	>	h2["pully"+suffix+"_nTrk"]->SetMarkerStyle(21);
299	>	h2["pully"+suffix+"_nTrk"]->SetMarkerColor(4);
300	>	h2["pully"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
301	>	h2["pully"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
302	>	h2["pullz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullz"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
303	>	h2["pullz"+suffix+"_nTrk"]->SetMarkerStyle(21);
304	>	h2["pullz"+suffix+"_nTrk"]->SetMarkerColor(4);
305	>	h2["pullz"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
306	>	h2["pullz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
307	>	} // End of   if(analyze_) {
308	>	suffix.clear();
309	>	} // End of Book histograms sensitive to data/mc
310	>	// Book MC only plots
311	>	if (!realData_) {
312	>	h["genPart_T"]      = subDir.make<TH1D>("genPart_T","t component of gen particles",300,-0.5,0.5);
313	>	h["genPart_T"]->GetXaxis()->SetTitle("t (nanosecond)");
314	>	h["genPart_cT"]     = subDir.make<TH1D>("genPart_cT","ct component of gen particles",300,-150.,150.);
315	>	h["genPart_cT"]->GetXaxis()->SetTitle("ct (mm)");
316	>	h["nsimPV"]         = subDir.make<TH1D>("nsimPV","Num of sim PV",51,-0.5,50.5);
317		}
318		}
319
150	–
320		PVStudy::~PVStudy()
321		{
322
#	Line 156 | Line 325 | PVStudy::~PVStudy()
325
326		}
327
328	<
328	>	void PVStudy::setRootStyle() {
329	>	//
330	>	gROOT->SetStyle("Plain");
331	>	gStyle->SetFillColor(1);
332	>	gStyle->SetOptDate();
333	>	//   gStyle->SetOptStat(1111110);
334	>	//   gStyle->SetOptFit(0111);
335	>	//   gStyle->SetPadTickX(1);
336	>	//   gStyle->SetPadTickY(1);
337	>	gStyle->SetMarkerSize(0.5);
338	>	gStyle->SetMarkerStyle(8);
339	>	gStyle->SetGridStyle(3);
340	>	//gStyle->SetPadGridX(1);
341	>	//gStyle->SetPadGridY(1);
342	>	gStyle->SetPaperSize(TStyle::kA4);
343	>	gStyle->SetStatW(0.25);             // width of statistics box; default is 0.19
344	>	//   gStyle->SetStatH(0.10);             // height of statistics box; default is 0.1
345	>	gStyle->SetStatFormat("6.4g");      // leave default format for now
346	>	gStyle->SetTitleSize(0.055, "");    // size for pad title; default is 0.02
347	>	gStyle->SetLabelSize(0.03, "XYZ");  // size for axis labels; default is 0.04
348	>	gStyle->SetStatFontSize(0.08);      // size for stat. box
349	>	gStyle->SetTitleFont(32, "XYZ");    // times-bold-italic font (p. 153) for axes
350	>	gStyle->SetTitleFont(32, "");       // same for pad title
351	>	gStyle->SetLabelFont(32, "XYZ");    // same for axis labels
352	>	gStyle->SetStatFont(32);            // same for stat. box
353	>	gStyle->SetLabelOffset(0.006, "Y"); // default is 0.005
354	>	//
355	>	return;
356	>	}
357		//
358		// member functions
359		//
#	Line 172 | Line 369 | std::vector<PVStudy::simPrimaryVertex> P
369		cout << "number of vertices " << evt->vertices_size() << endl;
370		}
371
372	<	int idx=0;
372	>	int idx=0; int npv=0;
373		for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
374	<	vitr != evt->vertices_end(); ++vitr )
375	<	{ // loop for vertex ...
376	<	HepMC::FourVector pos = (*vitr)->position();
377	<	//HepLorentzVector pos = (*vitr)->position();
378	<	if (pos.t()>0) { continue;} // ?
379	<
380	<	bool hasMotherVertex=false;
381	<	//      if(verbose_) cout << "mothers" << endl;
382	<	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
383	<	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
384	<	HepMC::GenVertex * mv=(*mother)->production_vertex();
385	<	if (mv) {hasMotherVertex=true;}
386	<	//      if(verbose_) {
387	<	//        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;
374	>	vitr != evt->vertices_end(); ++vitr ) { // loop for vertex ...
375	>	HepMC::FourVector pos = (*vitr)->position();
376	>	//HepLorentzVector pos = (*vitr)->position();
377	>
378	>	// t component of PV:
379	>	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
380	>	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
381	>	//        std::cout << "Status = " << (*mother)->status() << std::endl;
382	>	HepMC::GenVertex * mv=(*mother)->production_vertex();
383	>	if( ((*mother)->status() == 3) && (!mv)) {
384	>	//        std::cout << "npv= " << npv << std::endl;
385	>	if (npv == 0) {
386	>	h["genPart_cT"]->Fill(pos.t()); // mm
387	>	h["genPart_T"]->Fill(pos.t()/299.792458); // ns
388		}
389	+	npv++;
390	+	}
391	+	}
392	+	//       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
393	+
394	+	bool hasMotherVertex=false;
395	+	if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl;
396	+	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
397	+	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
398	+	HepMC::GenVertex * mv=(*mother)->production_vertex();
399	+	//        std::cout << "Status = " << (*mother)->status() << std::endl;
400	+	if (mv) {
401	+	hasMotherVertex=true;
402	+	if(!verbose_) break; //if verbose_, print all particles of gen vertices
403	+	}
404	+	if(verbose_) {
405	+	cout << "\t";
406	+	(*mother)->print();
407		}
408	+	}
409
410	<	if(!vp){
411	<	// this is a new vertex
412	<	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
413	<	simpv.push_back(sv);
414	<	vp=&simpv.back();
415	<	}else{
416	<	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
417	<	}
418	<	vp->genVertex.push_back((*vitr)->barcode());
419	<	// collect final state descendants
420	<	for ( HepMC::GenVertex::particle_iterator
421	<	daughter  = (*vitr)->particles_begin(HepMC::descendants);
422	<	daughter != (*vitr)->particles_end(HepMC::descendants);
423	<	++daughter ) {
424	<	if (isFinalstateParticle(*daughter)){
425	<	if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
426	<	== vp->finalstateParticles.end()){
427	<	vp->finalstateParticles.push_back((*daughter)->barcode());
428	<	HepMC::FourVector m=(*daughter)->momentum();
429	<	// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
430	<	// but adding FourVectors seems not to be foreseen
431	<	vp->ptot.setPx(vp->ptot.px()+m.px());
432	<	vp->ptot.setPy(vp->ptot.py()+m.py());
433	<	vp->ptot.setPz(vp->ptot.pz()+m.pz());
434	<	vp->ptot.setE(vp->ptot.e()+m.e());
435	<	vp->ptsq+=(m.perp())*(m.perp());
436	<	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
437	<	vp->nGenTrk++;
438	<	}
410	>	if(hasMotherVertex) continue;
411	>
412	>	// could be a new vertex, check  all primaries found so far to avoid multiple entries
413	>	const double mm2cm=0.1;
414	>	simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
415	>	simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
416	>	for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
417	>	v0!=simpv.end(); v0++){
418	>	if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
419	>	vp=&(*v0);
420	>	break;
421	>	}
422	>	}
423	>
424	>	if(!vp){
425	>	// this is a new vertex
426	>	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
427	>	simpv.push_back(sv);
428	>	vp=&simpv.back();
429	>	}else{
430	>	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
431	>	}
432	>	vp->genVertex.push_back((*vitr)->barcode());
433	>	// collect final state descendants
434	>	for ( HepMC::GenVertex::particle_iterator daughter  = (*vitr)->particles_begin(HepMC::descendants);
435	>	daughter != (*vitr)->particles_end(HepMC::descendants);
436	>	++daughter ) {
437	>	if (isFinalstateParticle(*daughter)){
438	>	if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
439	>	== vp->finalstateParticles.end()){
440	>	vp->finalstateParticles.push_back((*daughter)->barcode());
441	>	HepMC::FourVector m=(*daughter)->momentum();
442	>	// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
443	>	// but adding FourVectors seems not to be foreseen
444	>	vp->ptot.setPx(vp->ptot.px()+m.px());
445	>	vp->ptot.setPy(vp->ptot.py()+m.py());
446	>	vp->ptot.setPz(vp->ptot.pz()+m.pz());
447	>	vp->ptot.setE(vp->ptot.e()+m.e());
448	>	vp->ptsq+=(m.perp())*(m.perp());
449	>	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
450	>	vp->nGenTrk++;
451		}
452		}
453		}
241	–	idx++;
454		}
455	+	}
456		}
457		return simpv;
458		}
#	Line 319 | Line 532 | void
532		PVStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
533		{
534		using namespace edm;
535	<	using reco::TrackCollection;
536	<
537	<	ftree_->SetBranchAddress("residual",&fres);
538	<	ftree_->SetBranchAddress("error",&ferror);
539	<	ftree_->SetBranchAddress("nTrk",&fntrk);
540	<
541	<	Handle<TrackCollection> tracks;
542	<	iEvent.getByLabel(tracksLabel_,tracks);
543	<
544	<	Handle<reco::VertexCollection> recVtxs;
545	<	iEvent.getByLabel(vtxSample_, recVtxs);
535	>	using namespace reco;
536	>
537	>	if (verbose_)
538	>	cout<<"PVStudy::analyze():"<<endl;
539	>	//=======================================================
540	>	// Initialize Root-tuple variables if needed
541	>	//=======================================================
542	>	//if(saventuple_)
543	>	SetVarToZero();
544	>
545	>	//=======================================================
546	>	// Track accessors
547	>	//=======================================================
548	>	//trackCollection
549	>	static const reco::TrackCollection s_empty_trackColl;
550	>	const reco::TrackCollection *trackColl = &s_empty_trackColl;
551	>	edm::Handle<reco::TrackCollection> trackCollectionHandle;
552	>	iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle);
553	>	if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
554	>	trackColl = trackCollectionHandle.product();
555	>	} else {
556	>	cout << "trackCollection cannot be found -> using empty collection of same type." <<endl;
557	>	}
558	>	//splitTrackCollection1
559	>	static const reco::TrackCollection s_empty_splitTrackColl1;
560	>	const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
561	>	edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
562	>	iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle);
563	>	if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
564	>	splitTrackColl1 = splitTrackCollection1Handle.product();
565	>	} else {
566	>	cout << "splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
567	>	}
568	>	//splitTrackCollection2
569	>	static const reco::TrackCollection s_empty_splitTrackColl2;
570	>	const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
571	>	edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
572	>	iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle);
573	>	if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
574	>	splitTrackColl2 = splitTrackCollection2Handle.product();
575	>	} else {
576	>	cout << "splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
577	>	}
578	>
579	>	//=======================================================
580	>	// PVTX accessors
581	>	//=======================================================
582	>	//vertexCollection
583	>	static const reco::VertexCollection s_empty_vertexColl;
584	>	const reco::VertexCollection *vertexColl = &s_empty_vertexColl;
585	>	edm::Handle<reco::VertexCollection> vertexCollectionHandle;
586	>	iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle);
587	>	if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
588	>	vertexColl = vertexCollectionHandle.product();
589	>	} else {
590	>	cout << "vertexCollection cannot be found -> using empty collection of same type." <<endl;
591	>	}
592	>	//splitVertexCollection1
593	>	static const reco::VertexCollection s_empty_splitVertexColl1;
594	>	const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
595	>	edm::Handle<reco::VertexCollection> splitVertexCollection1Handle;
596	>	iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle);
597	>	if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
598	>	splitVertexColl1 = splitVertexCollection1Handle.product();
599	>	} else {
600	>	cout << "splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
601	>	}
602	>	//splitVertexCollection2
603	>	static const reco::VertexCollection s_empty_splitVertexColl2;
604	>	const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
605	>	edm::Handle<reco::VertexCollection> splitVertexCollection2Handle;
606	>	iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle);
607	>	if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
608	>	splitVertexColl2 = splitVertexCollection2Handle.product();
609	>	} else {
610	>	cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
611	>	}
612
613	+	if(verbose_) cout<<"End accessing the track and vertex collections"<<endl;
614	+
615	+	//=======================================================
616	+	// MC simvtx accessor
617	+	//=======================================================
618		if (!realData_) {
619	<	Handle<SimVertexContainer> simVtxs;
619	>	edm::Handle<SimVertexContainer> simVtxs;
620		iEvent.getByLabel( simG4_, simVtxs);
621
622	<	Handle<SimTrackContainer> simTrks;
622	>	edm::Handle<SimTrackContainer> simTrks;
623		iEvent.getByLabel( simG4_, simTrks);
624		}
625
626	<	if(verbose_) printRecVtxs(recVtxs);
627	<
626	>	//=======================================================
627	>	// GET PDT
628	>	//=======================================================
629		try{
630		iSetup.getData(pdt);
631		}catch(const Exception&){
632		cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
633		}
634
635	+	//=======================================================
636	+	// Fill trackparameters of the input tracks to pvtx fitter
637	+	//=======================================================
638	+	if(verbose_)
639	+	cout<<"Start filling track parameters of the input tracks to pvtx fitter."<<endl;
640	+	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype)
641	+	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
642	+	fillTrackHisto(trackColl, 0);
643	+	fillTrackHisto(splitTrackColl1, 1);
644	+	fillTrackHisto(splitTrackColl2, 2);
645	+	if(verbose_)
646	+	cout<<"End filling track parameters of the input tracks to pvtx fitter."<<endl;
647	+
648	+	//=======================================================
649	+	// Fill reconstructed vertices ootb
650	+	//=======================================================
651	+	if(verbose_)  {
652	+	cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
653	+	printRecVtxs(vertexCollectionHandle);
654	+	cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
655	+	printRecVtxs(splitVertexCollection1Handle);
656	+	cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
657	+	printRecVtxs(splitVertexCollection2Handle);
658	+	cout<<"Start filling pvtx parameters."<<endl;
659	+	}
660	+
661	+	nrecPV_ = int (vertexColl->size());
662	+	nrecPV1_spl_ = int (splitVertexColl1->size());
663	+	nrecPV2_spl_ = int (splitVertexColl2->size());
664	+
665	+	h["nrecPV"]->Fill(nrecPV_);
666	+	h["nrecPV1_spl"]->Fill(nrecPV1_spl_);
667	+	h["nrecPV2_spl"]->Fill(nrecPV2_spl_);
668	+	h["nrecPVDiff"]->Fill(double(nrecPV1_spl_)-double(nrecPV2_spl_));
669	+
670	+	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
671	+	fillTrackHistoInPV(vertexColl, 0, true, true);
672	+	fillTrackHistoInPV(splitVertexColl1, 1, false, true);
673	+	fillTrackHistoInPV(splitVertexColl2, 2, false, true);
674	+	if(verbose_)
675	+	cout<<"End filling pvtx parameters."<<endl;
676	+
677	+	//=======================================================
678	+	//           PrimaryVertex Matching
679	+	// 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
680	+	// 2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
681	+	// Assume the first match is the primary vertex,
682	+	// since vertexColl are sorted in decreasing order of sum(pT)
683	+	//=======================================================
684	+
685	+	if(verbose_)   cout<<"PVStudy::analyze(): matching pvtxs "<<endl;
686	+	reco::VertexCollection s_empty_matchedVertexColl1;
687	+	reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
688	+	matchedVertexColl1->reserve(splitVertexColl1->size());
689	+	reco::VertexCollection s_empty_matchedVertexColl2;
690	+	reco::VertexCollection *matchedVertexColl2 = &s_empty_matchedVertexColl2;
691	+	matchedVertexColl2->reserve(splitVertexColl2->size());
692	+
693	+	for (size_t ivtx = 0; ivtx<splitVertexCollection1Handle->size(); ++ivtx) {
694	+	reco::VertexRef recvtx1(splitVertexCollection1Handle, ivtx);
695	+	for (size_t jvtx = ivtx; jvtx < splitVertexCollection2Handle->size(); ++jvtx) {
696	+	reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
697	+	if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
698	+	if(verbose_) {
699	+	cout<<"The two splitted vertices match in Z: "<<endl;
700	+	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
701	+	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
702	+	}
703	+	int nTrkPV1 = recvtx1->tracksSize();
704	+	int nTrkPV2 = recvtx2->tracksSize();
705	+	double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
706	+	h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
707	+	h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
708	+	if(fabs(ntrkreldiff)<ntrkdiffcut_) {
709	+	if(verbose_) {
710	+	cout<<"(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
711	+	}
712	+	matchedVertexColl1->push_back(*recvtx1);
713	+	matchedVertexColl2->push_back(*recvtx2);
714	+	continue;
715	+	}
716	+	else
717	+	cout<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
718	+	}
719	+	else {
720	+	cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
721	+	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
722	+	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
723	+	}
724	+	}
725	+	}
726	+	if(verbose_)     cout<<"PVStudy::analyze(): End matching pvtxs"<<endl;
727	+
728	+	//=======================================================
729	+	//  Analyze matchedVertexColls
730	+	//=======================================================
731	+	if(verbose_) {
732	+	cout<<"Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
733	+	}
734	+
735	+
736	+	// Compare the reconstructed vertex position and calculate resolution/pulls
737	+	if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
738	+	//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
739	+	fillTrackHistoInPV(matchedVertexColl1, 1, true, false);
740	+	fillTrackHistoInPV(matchedVertexColl2, 2, true, false);
741	+
742	+	reco::VertexCollection::const_iterator v1;
743	+	reco::VertexCollection::const_iterator v2;
744	+	nrecPV_twovtx_ = 0;
745	+	for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
746	+	v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
747	+	++v1, ++v2) {
748	+	fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
749	+	fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
750	+	fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
751	+	ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
752	+	ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
753	+	ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
754	+	fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
755	+
756	+	h["deltax"]->Fill( fres_[0] );
757	+	h["deltay"]->Fill( fres_[1] );
758	+	h["deltaz"]->Fill( fres_[2] );
759	+	h["pullx"]->Fill( fres_[0]/ferror_[0]);
760	+	h["pully"]->Fill( fres_[1]/ferror_[1]);
761	+	h["pullz"]->Fill( fres_[2]/ferror_[2]);
762	+	h["errPVx"]->Fill( ferror_[0] );
763	+	h["errPVy"]->Fill( ferror_[1] );
764	+	h["errPVz"]->Fill( ferror_[2] );
765	+	pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
766	+	error(ferror_[0], ferror_[1], ferror_[2]),
767	+	fntrk_));
768	+	// Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
769	+	fillHisto(res(fres_[0], fres_[1], fres_[2]),
770	+	error(ferror_[0], ferror_[1], ferror_[2]),
771	+	fntrk_,0);
772	+	// Fill the ntuple variables
773	+	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
774	+	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
775	+	sumptsq1_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v1);
776	+	sumptsq2_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v2);
777	+	nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
778	+	deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
779	+	deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
780	+	deltaz_twovtx_[nrecPV_twovtx_] = fres_[2];
781	+	errx_twovtx_[nrecPV_twovtx_] = ferror_[0];
782	+	erry_twovtx_[nrecPV_twovtx_] = ferror_[1];
783	+	errz_twovtx_[nrecPV_twovtx_] = ferror_[2];
784	+	pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
785	+	pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
786	+	pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];
787	+	nrecPV_twovtx_++;
788	+	} // End of analyzing res/pull
789	+	} // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
790	+	else
791	+	cout<<"WARNING: Cannot find matching pvtxs"<<endl;
792	+
793	+	if(verbose_)     cout<<"End analyzing the matchedVertexColl"<<endl;
794	+
795	+	//=======================================================
796	+	// Fill simulated vertices
797	+	//=======================================================
798		if (!realData_) {
799		bool MC=false;
800		Handle<HepMCProduct> evtMC;
#	Line 362 | Line 810 | PVStudy::analyze(const edm::Event& iEven
810		}
811		MC=true;
812		}
813	<
813	>
814		if(MC){
815		// make a list of primary vertices:
816		std::vector<simPrimaryVertex> simpv;
817		simpv=getSimPVs(evtMC,"");
818		//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
819	<
820	<	int nsimtrk=0;
821	<	int nrectrk=0;
819	>	h["nsimPV"]->Fill(simpv.size());
820	>	nsimPV_ = simpv.size();
821	>
822	>	int isimpv = 0;
823		for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
824	<	vsim!=simpv.end(); vsim++){
825	<	nsimtrk+=vsim->nGenTrk;
826	<	// look for a matching reconstructed vertex
827	<	vsim->recVtx=NULL;
828	<	for(reco::VertexCollection::const_iterator vrec=recVtxs->begin();
829	<	vrec!=recVtxs->end(); ++vrec){
830	<	nrectrk=vrec->tracksSize();
831	<	if(verbose_){
832	<	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
833	<	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
834	<	}
835	<	if ( matchVertex(*vsim,*vrec) ){
836	<	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
837	<	if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
838	<	|| (!vsim->recVtx) )
839	<	{
840	<	vsim->recVtx=&(*vrec);
841	<	}
842	<	}
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());
824	>	vsim!=simpv.end(); vsim++, isimpv++){
825	>	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
826	>	simx_[isimpv] = vsim->x;
827	>	simy_[isimpv] = vsim->y;
828	>	simz_[isimpv] = vsim->y;
829	>	simptsq_[isimpv] = vsim->ptsq;
830	>
831	>	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
832	>	fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
833	>	fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
834	>	fillMCHisto(vsim, isimpv, vertexColl, 3);
835	>	}
836	>	} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
837	>	} // End of   if (!realData_) {
838	>
839	>	// Finally fill the ftree_
840	>	if(saventuple_) {
841	>	ftree_->Fill();
842	>	pvtxtree_->Fill();
843		}
472	–
844		}
845
475	–
846		// ------------ method called once each job just before starting event loop  ------------
847		void
848		PVStudy::beginJob()
#	Line 482 | Line 852 | PVStudy::beginJob()
852		// ------------ method called once each job just after ending the event loop  ------------
853		void
854		PVStudy::endJob() {
855	<	double sqt2 = sqrt(2.);
856	<	if (!realData_) {
857	<	if (verbose_) cout << "Analyzing PV info" << endl;
858	<	for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
859	<	if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
860	<	PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk);
861	<	if (analyze_) {
862	<	//      if ( ipv.res_.x() > 0 ) h_resx_nTrk->Fill(ntrk,ipv.res_.x()/sqt2, 1.);
863	<	//      if ( ipv.res_.y() > 0 ) h_resy_nTrk->Fill(ntrk,ipv.res_.y()/sqt2, 1.);
864	<	//      if ( ipv.res_.z() > 0 ) h_resz_nTrk->Fill(ntrk,ipv.res_.z()/sqt2, 1.);
865	<	if ( ipv.res_.x() > 0 ) h_resx_nTrk->Fill(ntrk,ipv.res_.x(), 1.);
866	<	if ( ipv.res_.y() > 0 ) h_resy_nTrk->Fill(ntrk,ipv.res_.y(), 1.);
867	<	if ( ipv.res_.z() > 0 ) h_resz_nTrk->Fill(ntrk,ipv.res_.z(), 1.);
868	<	if ( ipv.pull_.x() > 0 ) h_pullx_nTrk->Fill(ntrk,ipv.pull_.x(), 1.);
499	<	if ( ipv.pull_.y() > 0 ) h_pully_nTrk->Fill(ntrk,ipv.pull_.y(), 1.);
500	<	if ( ipv.pull_.z() > 0 ) h_pullz_nTrk->Fill(ntrk,ipv.pull_.z(), 1.);
501	<	}
855	>	if (verbose_) cout << "Analyzing PV info" << endl;
856	>	// Looping through the datamodes available
857	>	for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
858	>	string suffix;
859	>	if(verbose_)  cout<<"datamode = "<< *it<<endl;
860	>	switch(*it) {
861	>	case 0: suffix = "";
862	>	break;
863	>	case 1: suffix = "_spl1_mct";
864	>	break;
865	>	case 2: suffix = "_spl2_mct";
866	>	break;
867	>	case 3: suffix = "_mct";
868	>	break;
869		}
870	+	if(analyze_) {
871	+	for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
872	+	PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
873	+	if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
874	+	if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
875	+	if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
876	+	if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
877	+	if ( ipv.pull_.y() > 0 ) h2["pully"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
878	+	if ( ipv.pull_.z() > 0 ) h2["pullz"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.z(), 1.);
879	+	}
880	+	}
881	+	suffix.clear();
882	+	}
883	+	if(saventuple_) {
884	+	file_->cd();
885	+	ftree_->Write();
886	+	pvtxtree_->Write();
887	+	file_->Close();
888		}
504	–	file_->cd();
505	–	ftree_->Write();
506	–	file_->Close();
889		}
890
891	+
892	+
893	+	// Match a recovertex with a simvertex
894	+	// ? Should the cut be parameter dependent?
895	+	// cut can be within n sigma of the vertex.zerror()
896		bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
897		const reco::Vertex & vrec)
898		{
899		return (fabs(vsim.z-vrec.z())<0.0500); // =500um
900		}
901
902	+	// Match two reconstructed vertices
903	+	bool PVStudy::matchVertex( const reco::VertexRef & vrec1,
904	+	const reco::VertexRef & vrec2,
905	+	int zsigncut)
906	+	{
907	+	double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
908	+	if(verbose_)
909	+	cout<<"The matching criteria in z is "<<cut<<endl;
910	+	return (fabs(vrec1->z()-vrec2->z())<cut);
911	+	}
912	+
913	+
914		bool PVStudy::isResonance(const HepMC::GenParticle * p){
915		double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
916		if(verbose_) cout << "isResonance   " << p->pdg_id() << " " << ctau << endl;
#	Line 583 | Line 982 | void PVStudy::printSimTrks(const Handle<
982		}
983		}
984
985	<	void PVStudy::fillHisto(res r, pull p, int ntk){
985	>	void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
986		stringstream ss;
987		ss << ntk;
988		string suffix = ss.str();
989	+	if(datamode == 0 )  suffix = "_" + suffix;
990	+	if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
991	+	if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
992	+	if(datamode == 3 )  suffix = "_mct_" + suffix;
993	+
994		if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
995		// Fill histograms of res and pull of ntk
996	<	h["resx_"+suffix]->Fill(r.x());
997	<	h["resy_"+suffix]->Fill(r.y());
998	<	h["resz_"+suffix]->Fill(r.z());
999	<	h["pullx_"+suffix]->Fill(p.x());
1000	<	h["pully_"+suffix]->Fill(p.y());
1001	<	h["pullz_"+suffix]->Fill(p.z());
996	>	h["deltax"+suffix]->Fill(r.x());
997	>	h["deltay"+suffix]->Fill(r.y());
998	>	h["deltaz"+suffix]->Fill(r.z());
999	>	h["pullx"+suffix]->Fill(r.x()/er.x());
1000	>	h["pully"+suffix]->Fill(r.y()/er.y());
1001	>	h["pullz"+suffix]->Fill(r.z()/er.z());
1002	>	h["errPVx"+suffix]->Fill(er.x());
1003	>	h["errPVy"+suffix]->Fill(er.y());
1004	>	h["errPVz"+suffix]->Fill(er.z());
1005		}
1006
1007	<	PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk) {
1007	>	PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1008		map<int,double> data;
1009		data.clear();
1010	<	double fpar[2];
1010	>	double fpar[2] = {-999,-999};
1011		double cm2um = 10000;
1012		stringstream ss;
1013		ss << ntk;
1014	<	string suffix = ss.str();
1015	<
1014	>	string suffix = ss.str();
1015	>	if(datamode == 0 )  suffix = "_" + suffix;
1016	>	if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1017	>	if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1018	>	if(datamode == 3 )  suffix = "_mct_" + suffix;
1019	>
1020		if(analyze_) {
1021		for ( int i = 0; i < 6; ++i) {
611	–	cout << "Here" << endl;
1022		switch (i) {
1023		case 0:
1024	<	fit(h["resx_"+suffix], fpar);
1024	>	fit(h["deltax"+suffix], fpar);
1025		data[i] = fpar[0]*cm2um;
1026		break;
1027		case 1:
1028	<	fit(h["resy_"+suffix], fpar);
1028	>	fit(h["deltay"+suffix], fpar);
1029		data[i] = fpar[0]*cm2um;
1030		break;
1031		case 2:
1032	<	fit(h["resz_"+suffix], fpar);
1032	>	fit(h["deltaz"+suffix], fpar);
1033		data[i] = fpar[0]*cm2um;
1034		break;
1035		case 3:
1036	<	fit(h["pullx_"+suffix], fpar);
1036	>	fit(h["pullx"+suffix], fpar);
1037		data[i] = fpar[0];
1038		break;
1039		case 4:
1040	<	fit(h["pully_"+suffix], fpar);
1040	>	fit(h["pully"+suffix], fpar);
1041		data[i] = fpar[0];
1042		break;
1043		case 5:
1044	<	fit(h["pullz_"+suffix], fpar);
1044	>	fit(h["pullz"+suffix], fpar);
1045		data[i] = fpar[0];
1046		break;
1047		}
#	Line 650 | Line 1060 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1060		ntk);
1061		}
1062
1063	<	void PVStudy::fit(TH1 *&hdist, double data[]){
1063	>	void PVStudy::fit(TH1 *&hdist, double fitPars[]){
1064		TAxis *axis0 = hdist->GetXaxis();
1065		int nbin = axis0->GetLast();
1066		double nOF = hdist->GetBinContent(nbin+1);
1067		double nUF = hdist->GetBinContent(0);
1068	+	//   double fitRange = axis0->GetBinUpEdge(nbin);
1069	+	double fitRange = axis0->GetXmax();
1070	+	double sigMax[2] = {0.,0.};
1071	+
1072		if ( verbose_ ){
1073		cout << "Last bin = " << nbin;
1074		cout << "; hdist: Overflow Entries = " << nOF;
1075		cout << "; Underflow Entries = " << nUF;
1076	<	cout << "; hdist->GetEntries() = " << hdist->GetEntries() << endl;
1076	>	cout << "; hdist->GetEntries() = " << hdist->GetEntries();
1077	>	cout << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1078		}
1079	<	if (hdist->GetEntries() - nOF - nUF >= 5) { // FIXME: for reasonable Gaussian fit
1080	<	hdist->Fit("gaus","Q");
1081	<	TF1 *fgaus = hdist->GetFunction("gaus");
1079	>	if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit
1080	>	for (int bi = 0; bi < nbin; bi++) {
1081	>	if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) {
1082	>	sigMax[0] = axis0->GetBinLowEdge(bi);
1083	>	if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1084	>	}
1085	>	if ( (axis0->GetBinLowEdge(bi) >= 0) && (hdist->GetBinContent(bi) > 0) ) {
1086	>	sigMax[0] = axis0->GetBinUpEdge(bi);
1087	>	if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1088	>	}
1089	>	}
1090	>	if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1091	>
1092	>	TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1093	>	fgaus->SetParameter(1, 0.);
1094	>	fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
1095	>	fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
1096	>	fgaus->SetLineColor(4);
1097	>	hdist->Fit(fgaus,"Q");
1098	>	gPad->Update();
1099	>	TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
1100	>	s->SetOptStat(1111111);
1101	>	s->SetOptFit(0111);
1102	>	gPad->Update();
1103		if (verbose_) cout << "got function" << endl;
1104	<	data[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1104	>	fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1105	>	}
1106	>	else
1107	>	fitPars[0] = -999;
1108	>	}
1109	>
1110	>	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype) {
1111	>	string suffix;
1112	>	suffix = ""; // for unsplit tracks
1113	>	if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1114	>	if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1115	>
1116	>	h["nTrk"+suffix]->Fill(trackColl->size());
1117	>	for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1118	>	itTrack != trackColl->end();
1119	>	++itTrack) {
1120	>	h["trkPt"+suffix]->Fill(itTrack->pt());
1121	>	h["trkDxy"+suffix]->Fill(itTrack->dxy());
1122	>	h["trkDz"+suffix]->Fill(itTrack->dz());
1123	>	h["trkEta"+suffix]->Fill(itTrack->eta());
1124	>	h["trkPhi"+suffix]->Fill(itTrack->phi());
1125	>	}
1126	>	}
1127	>
1128	>	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
1129	>	string suffix;
1130	>	suffix = ""; // for unsplit tracks
1131	>	if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1132	>	if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1133	>	int ivtx = 0;
1134	>	for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1135	>	v!=vertexColl->end(); ++v, ++ivtx) {
1136	>	if(fillNtuple) {
1137	>	if(datatype == 0) {
1138	>	nTrkPV_[ivtx] = v->tracksSize();
1139	>	recx_[ivtx] = v->x();
1140	>	recy_[ivtx] = v->y();
1141	>	recz_[ivtx] = v->z();
1142	>	recx_err_[ivtx] = v->xError();
1143	>	recy_err_[ivtx] = v->yError();
1144	>	recz_err_[ivtx] = v->zError();
1145	>	}
1146	>	if(datatype == 1) {
1147	>	nTrkPV1_spl_[ivtx] = v->tracksSize();
1148	>	recx1_spl_[ivtx] = v->x();
1149	>	recy1_spl_[ivtx] = v->y();
1150	>	recz1_spl_[ivtx] = v->z();
1151	>	recx1_err_spl_[ivtx] = v->xError();
1152	>	recy1_err_spl_[ivtx] = v->yError();
1153	>	recz1_err_spl_[ivtx] = v->zError();
1154	>	}
1155	>	if(datatype == 2) {
1156	>	nTrkPV2_spl_[ivtx] = v->tracksSize();
1157	>	recx2_spl_[ivtx] = v->x();
1158	>	recy2_spl_[ivtx] = v->y();
1159	>	recz2_spl_[ivtx] = v->z();
1160	>	recx2_err_spl_[ivtx] = v->xError();
1161	>	recy2_err_spl_[ivtx] = v->yError();
1162	>	recz2_err_spl_[ivtx] = v->zError();
1163	>	}
1164	>	}
1165	>
1166	>	if(fillHisto) {
1167	>	h["nTrkPV"+suffix]->Fill(v->tracksSize());
1168	>	try {
1169	>	for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
1170	>	t!=v->tracks_end(); t++) {
1171	>	// illegal charge
1172	>	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1173	>	//    h["trkPtPV"]->Fill(0.);
1174	>	}
1175	>	else {
1176	>	h["trkPtPV"+suffix]->Fill((**t).pt());
1177	>	h["trkDxyPV"+suffix]->Fill((**t).dxy());
1178	>	h["trkDzPV"+suffix]->Fill((**t).dz());
1179	>	h["trkEtaPV"+suffix]->Fill((**t).eta());
1180	>	h["trkPhiPV"+suffix]->Fill((**t).phi());
1181	>	}
1182	>	}
1183	>	}
1184	>	catch (...) {
1185	>	// exception thrown when trying to use linked track
1186	>	//        h["trkPtPV"]->Fill(0.);
1187	>	}
1188	>	}
1189	>	}
1190	>	}
1191	>
1192	>	void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
1193	>	{
1194	>	std::string suffix;
1195	>	// Set the vertexColl and histogram suffix according to datamode
1196	>	if (datatype == 1) {
1197	>	suffix = "_spl1_mct";
1198	>	nrecPV_spl1_mct_ = 0;
1199	>	}
1200	>	if (datatype == 2) {
1201	>	suffix = "_spl2_mct";
1202	>	nrecPV_spl2_mct_ = 0;
1203	>	}
1204	>	if (datatype == 3) {
1205	>	suffix = "_mct";
1206	>	nrecPV_mct_ = 0;
1207	>	}
1208	>
1209	>	//========================================================
1210	>	//  look for a matching reconstructed vertex in vertexColl
1211	>	//========================================================
1212	>
1213	>	for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1214	>	vrec!=vtxColl->end(); ++vrec){
1215	>	int nrectrk = vrec->tracksSize();
1216	>	vsim->recVtx=NULL;
1217	>
1218	>	if(verbose_){
1219	>	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1220	>	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1221	>	}
1222	>	if ( matchVertex(*vsim,*vrec) ) {
1223	>	vsim->recVtx=&(*vrec);
1224	>
1225	>	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1226	>	//if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
1227	>	//|| (!vsim->recVtx) )
1228	>	//vsim->recVtx=&(*vrec);
1229	>	//}
1230	>	if(verbose_)
1231	>	cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1232	>
1233	>	double fres_mct[3];
1234	>	double ferror_mct[3];
1235	>
1236	>	fres_mct[0] = vsim->recVtx->x()-vsim->x;
1237	>	fres_mct[1] = vsim->recVtx->y()-vsim->y;
1238	>	fres_mct[2] = vsim->recVtx->z()-vsim->z;
1239	>	ferror_mct[0] = vsim->recVtx->xError();
1240	>	ferror_mct[1] = vsim->recVtx->yError();
1241	>	ferror_mct[2] = vsim->recVtx->zError();
1242	>
1243	>	h["deltax"+suffix]->Fill( fres_mct[0] );
1244	>	h["deltay"+suffix]->Fill( fres_mct[1] );
1245	>	h["deltaz"+suffix]->Fill( fres_mct[2] );
1246	>	h["pullx"+suffix]->Fill( fres_mct[0]/ferror_mct[0] );
1247	>	h["pully"+suffix]->Fill( fres_mct[1]/ferror_mct[1] );
1248	>	h["pullz"+suffix]->Fill( fres_mct[2]/ferror_mct[2] );
1249	>	h["errPVx"+suffix]->Fill( ferror_mct[0] );
1250	>	h["errPVy"+suffix]->Fill( ferror_mct[1] );
1251	>	h["errPVz"+suffix]->Fill( ferror_mct[2] );
1252	>	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1253	>	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1254	>	nrectrk));
1255	>	// Fill histo according to its track multiplicity
1256	>	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1257	>	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1258	>	nrectrk, datatype);
1259	>
1260	>	if(saventuple_) {
1261	>	//Fill the values for variables in ftree_
1262	>	if(datatype == 1) {
1263	>	nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1264	>	deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1265	>	deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1266	>	deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1267	>	pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
1268	>	pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1269	>	pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1270	>	nrecPV_spl1_mct_++;
1271	>	}
1272	>
1273	>	if(datatype == 2) {
1274	>	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1275	>	deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1276	>	deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1277	>	deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1278	>	pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
1279	>	pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1280	>	pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1281	>	nrecPV_spl2_mct_++;
1282	>	}
1283	>	if(datatype == 3) {
1284	>	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1285	>	deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1286	>	deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1287	>	deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1288	>	pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
1289	>	pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1290	>	pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1291	>	nrecPV_mct_++;
1292	>	}
1293	>	} // End of  if(saventuple_) {
1294	>	} //  if ( matchVertex(*vsim,*vrec) ) {
1295	>	else {  // no rec vertex found for this simvertex
1296	>	if(verbose_)
1297	>	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1298	>	}
1299		}
1300		}
1301
1302	<	//define this as a plug-in
1303	<	DEFINE_FWK_MODULE(PVStudy);
1302	>
1303	>	void PVStudy::SetVarToZero() {
1304	>	//Greg's variables
1305	>	fntrk_ = 0;
1306	>	//pvtx position (x,y,z) residual and error
1307	>	for(int i = 0; i<3;i++)    {
1308	>	fres_[i] = 0;
1309	>	ferror_[i] = 0;
1310	>	}
1311	>
1312	>	//Yanyan's variables
1313	>	// Number of reconstructed vertices  // for GeneralTracks
1314	>	nrecPV_ = 0;
1315	>	nrecPV1_spl_ = 0;
1316	>	nrecPV2_spl_ = 0;
1317	>	nrecPV_twovtx_ = 0;
1318	>	nsimPV_ = 0;
1319	>	nrecPV_mct_ = 0;
1320	>	nrecPV_spl1_mct_ = 0;
1321	>	nrecPV_spl2_mct_ = 0;
1322	>
1323	>	for (int i = 0; i < nMaxPVs_; i++) {
1324	>	nTrkPV_[i] = 0; // Number of tracks in the pvtx
1325	>	recx_[i] = 0;
1326	>	recy_[i] = 0;
1327	>	recz_[i] = 0;
1328	>	recx_err_[i] = 0;
1329	>	recy_err_[i] = 0;
1330	>	recz_err_[i] = 0;
1331	>
1332	>	nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx
1333	>	recx1_spl_[i] = 0;
1334	>	recy1_spl_[i] = 0;
1335	>	recz1_spl_[i] = 0;
1336	>	recx1_err_spl_[i] = 0;
1337	>	recy1_err_spl_[i] = 0;
1338	>	recz1_err_spl_[i] = 0;
1339	>
1340	>	nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx
1341	>
1342	>	recx2_spl_[i] = 0;
1343	>	recy2_spl_[i] = 0;
1344	>	recz2_spl_[i] = 0;
1345	>	recx2_err_spl_[i] = 0;
1346	>	recy2_err_spl_[i] = 0;
1347	>	recz2_err_spl_[i] = 0;
1348	>
1349	>	nTrkPV1_spl_twovtx_[i] = 0;
1350	>	nTrkPV2_spl_twovtx_[i] = 0;
1351	>	nTrkPV_twovtx_[i] = 0;
1352	>	sumptsq1_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1353	>	sumptsq2_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
1354	>	deltax_twovtx_[i] = 0;
1355	>	deltay_twovtx_[i] = 0;
1356	>	deltaz_twovtx_[i] = 0;
1357	>	errx_twovtx_[i] = 0;
1358	>	erry_twovtx_[i] = 0;
1359	>	errz_twovtx_[i] = 0;
1360	>	pullx_twovtx_[i] = 0;
1361	>	pully_twovtx_[i] = 0;
1362	>	pullz_twovtx_[i] = 0;
1363	>
1364	>	nsimTrkPV_[i] = 0;
1365	>	simx_[i] = 0;
1366	>	simy_[i] = 0;
1367	>	simz_[i] = 0;
1368	>	simptsq_[i] = 0;
1369	>
1370	>	deltax_mct_[i] = 0;
1371	>	deltay_mct_[i] = 0;
1372	>	deltaz_mct_[i] = 0;
1373	>	pullx_mct_[i] = 0;
1374	>	pully_mct_[i] = 0;
1375	>	pullz_mct_[i] = 0;
1376	>
1377	>	deltax_spl1_mct_[i] = 0;
1378	>	deltay_spl1_mct_[i] = 0;
1379	>	deltaz_spl1_mct_[i] = 0;
1380	>	pullx_spl1_mct_[i] = 0;
1381	>	pully_spl1_mct_[i] = 0;
1382	>	pullz_spl1_mct_[i] = 0;
1383	>
1384	>	deltax_spl2_mct_[i] = 0;
1385	>	deltay_spl2_mct_[i] = 0;
1386	>	deltaz_spl2_mct_[i] = 0;
1387	>	pullx_spl2_mct_[i] = 0;
1388	>	pully_spl2_mct_[i] = 0;
1389	>	pullz_spl2_mct_[i] = 0;
1390	>	}
1391	>
1392	>	}
1393	>
1394	>	double PVStudy::sumPtSquared(const reco::Vertex & v)  {
1395	>	double sum = 0.;
1396	>	double pT;
1397	>	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1398	>	pT = (**it).pt();
1399	>	sum += pT*pT;
1400	>	}
1401	>	return sum;
1402	>	}

Diff Legend

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