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.2 by jengbou, Thu Sep 3 19:29:38 2009 UTC vs.
Revision 1.8 by jengbou, Fri Oct 16 00:11:30 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 50 | Line 51 | Implementation:
51		#include <SimDataFormats/Track/interface/SimTrackContainer.h>
52
53		//root
54	<	#include "TROOT.h"
55	<	#include "TF1.h"
56	<	#include "TString.h"
57	<	#include "TStyle.h"
58	<
54	>	#include <TROOT.h>
55	>	#include <TF1.h>
56	>	#include <TString.h>
57	>	#include <TStyle.h>
58	>	#include <TPaveStats.h>
59	>	#include <TPad.h>
60
61		using namespace std;
62
63	<	PVStudy::PVStudy(const edm::ParameterSet& iConfig):
62	<	nTrkMin_(10),nTrkMax_(150)
63	>	PVStudy::PVStudy(const edm::ParameterSet& iConfig)
64		{
65	<	simG4_           = iConfig.getParameter<edm::InputTag>( "simG4" );
66	<	tracksLabel_     = iConfig.getUntrackedParameter<edm::InputTag>("tracks");
67	<	vtxSample_       = iConfig.getUntrackedParameter<edm::InputTag>("vtxSample");
68	<	verbose_         = iConfig.getUntrackedParameter<bool>("verbose", false);
69	<	realData_        = iConfig.getUntrackedParameter<bool>("realData",false);
70	<	analyze_         = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
65	>	//=======================================================================
66	>	// Get configuration for TrackTupleMaker
67	>	//=======================================================================
68	>	simG4_                     = iConfig.getParameter<edm::InputTag>( "simG4" );
69	>	trackCollectionTag_        = iConfig.getUntrackedParameter<edm::InputTag>("trackCollection");
70	>	splitTrackCollection1Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection1");
71	>	splitTrackCollection2Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
72	>	vertexCollectionTag_       = iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection");
73	>	splitVertexCollection1Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
74	>	splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
75	>	verbose_                   = iConfig.getUntrackedParameter<bool>("verbose",false);
76	>	realData_                  = iConfig.getUntrackedParameter<bool>("realData",false);
77	>	analyze_                   = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
78	>	saventuple_                = iConfig.getUntrackedParameter<bool>("saventuple",false);
79	>	outputfilename_            = iConfig.getUntrackedParameter<string>("OutputFileName");
80	>	ntrkdiffcut_               = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
81	>	nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
82	>	nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
83
84		//now do what ever initialization is needed
85	<	pvanainfo.clear();
85	>	pvinfo.clear();
86
74	–	edm::Service<TFileService> fs;
75	–	//   TFileDirectory subDir = fs->mkdir( "Summary" );
76	–	//   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
87
88	<	h_nTrk           = fs->make<TH1F>("nTrk","Num of total tracks",300,0,300);
89	<	h_nTrkPV         = fs->make<TH1F>("nTrkPV","Num of Tracks from PV",300,0,300);
90	<	h_trkPt          = fs->make<TH1D>("trkPt","Pt of all tracks",100,0,100);
91	<	h_trkPtPV        = fs->make<TH1D>("trkPtPV","Pt dist of tracks from PV",100,0,100);
88	>	// Specify the data mode vector
89	>	if(realData_) datamode.push_back(0);
90	>	else {
91	>	datamode.push_back(0);
92	>	datamode.push_back(1);
93	>	datamode.push_back(2);
94	>	datamode.push_back(3);
95	>	}
96
97	<	if (!realData_) {
98	<	h_deltax       = fs->make<TH1D>("deltax","x-distance (Rec - Sim)",800,-0.04,0.04);
99	<	h_deltay       = fs->make<TH1D>("deltay","y-distance (Rec - Sim)",800,-0.04,0.04);
100	<	h_deltaz       = fs->make<TH1D>("deltaz","z-distance (Rec - Sim)",800,-0.04,0.04);
101	<	h_pullx        = fs->make<TH1D>("pullx","x-pull (Rec - Sim)",1000,-25.,25.);
102	<	h_pully        = fs->make<TH1D>("pully","y-pull (Rec - Sim)",1000,-25.,25.);
103	<	h_pullz        = fs->make<TH1D>("pullz","z-pull (Rec - Sim)",1000,-25.,25.);
97	>	// Create ntuple files if needed
98	>	if(saventuple_) {
99	>	file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
100	>	ftree_ = new TTree("mytree","mytree");
101	>	ftree_->AutoSave();
102	>	ftree_->Branch("residual",&fres_,"fres_[3]/D");
103	>	ftree_->Branch("error",&ferror_,"ferror_[3]/D");
104	>	ftree_->Branch("nTrkPV",&fntrk_,"fntrk_/I");
105	>	ftree_->Branch("nrecPV",&nrecPV_,"nrecPV_/I");
106	>	ftree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl_/I");
107	>	ftree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl_/I");
108	>
109	>	// If it is MC, get the resA, resB
110	>	if(!realData_) {
111	>	ftree_->Branch("residual_spl1_mct",&fres_spl1_mct_,"fres_spl1_mct_[3]/D");
112	>	ftree_->Branch("error_spl1_mct",&ferror_spl1_mct_,"ferror_spl1_mct_[3]/D");
113	>	ftree_->Branch("nTrkPV_spl1_mct",&fntrk_spl1_mct_,"fntrk_spl1_mct_/I");
114	>	ftree_->Branch("residual_spl2_mct",&fres_spl2_mct_,"fres_spl2_mct_[3]/D");
115	>	ftree_->Branch("error_spl2_mct",&ferror_spl2_mct_,"ferror_spl2_mct_[3]/D");
116	>	ftree_->Branch("nTrkPV_spl2_mct",&fntrk_spl2_mct_,"fntrk_spl2_mct_/I");
117	>	ftree_->Branch("residual_mct",&fres_mct_,"fres_mct_[3]/D");
118	>	ftree_->Branch("error_mct",&ferror_mct_,"ferror_mct_[3]/D");
119	>	ftree_->Branch("nTrkPV_mct",&fntrk_mct_,"fntrk_mct_/I");
120	>	}
121	>	}
122
123	<	// Summary of analysis:
124	<	if (analyze_) {
125	<	h_resx_nTrk    = fs->make<TH2D>("resx_ntrk","x-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
126	<	h_resx_nTrk->SetMarkerStyle(21);
127	<	h_resx_nTrk->GetYaxis()->SetTitle("#mum");
128	<	h_resy_nTrk    = fs->make<TH2D>("resy_ntrk","y-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
129	<	h_resy_nTrk->SetMarkerStyle(21);
130	<	h_resy_nTrk->GetYaxis()->SetTitle("#mum");
131	<	h_resz_nTrk    = fs->make<TH2D>("resz_ntrk","z-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
132	<	h_resz_nTrk->SetMarkerStyle(21);
133	<	h_resz_nTrk->GetYaxis()->SetTitle("#mum");
134	<	h_pullx_nTrk   = fs->make<TH2D>("pullx_ntrk","x-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
135	<	h_pullx_nTrk->SetMarkerStyle(21);
136	<	h_pullx_nTrk->SetBit(TH1::kCanRebin);
137	<	h_pullx_nTrk->GetYaxis()->SetTitle("cm");
138	<	h_pully_nTrk   = fs->make<TH2D>("pully_ntrk","y-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
139	<	h_pully_nTrk->SetMarkerStyle(21);
140	<	h_pully_nTrk->SetBit(TH1::kCanRebin);
141	<	h_pully_nTrk->GetYaxis()->SetTitle("cm");
142	<	h_pullz_nTrk   = fs->make<TH2D>("pullz_ntrk","z-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
143	<	h_pullz_nTrk->SetMarkerStyle(21);
144	<	h_pullz_nTrk->SetBit(TH1::kCanRebin);
145	<	h_pullz_nTrk->GetYaxis()->SetTitle("cm");
123	>	edm::Service<TFileService> fs;
124	>	TFileDirectory subDir = fs->mkdir( "Summary" );
125	>	TFileDirectory subDir1 = fs->mkdir( "Others" );
126	>	//   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
127	>
128	>	setRootStyle();
129	>
130	>	//Book histograms for track and pvtx parameters for each splitted pvtx
131	>	//i=0 : x (as in unsplit track collection)
132	>	//i=1 : x1_spl_
133	>	//i=2 : x2_spl_
134	>
135	>	for(int i=0;i<3;i++)
136	>	{
137	>	string suffix;
138	>	if(i == 0) suffix = "";
139	>	else {
140	>	stringstream ss;
141	>	ss << i;
142	>	suffix =ss.str()+"_spl";
143	>	}
144	>	h["nTrk"+suffix]   = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
145	>	h["trkPt"+suffix]  = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
146	>	h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#Eta of rec tracks "+suffix),100,-3,3);
147	>	h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
148	>	h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-5,5);
149	>	h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
150	>
151	>	h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
152	>	h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
153	>	h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#Eta of rec tracks in PV"+suffix),100,-3,3);
154	>	h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
155	>	h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);
156	>	h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
157	>	h["nrecPV"+suffix]   = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);
158	>	suffix.clear();
159		}
160	+	h["nrecPVDiff"]     = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
161	+	h["nTrkPVDiff"]     = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
162	+	h["nTrkPVRelDiff"]  = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
163	+
164	+	//Book histograms sensitive to data/mc
165	+	for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
166	+	string suffix;
167	+	if(verbose_)  cout<<"datamode = "<< *it<<endl;
168	+	switch(*it) {
169	+	case 0: suffix = "";
170	+	break;
171	+	case 1: suffix = "_spl1_mct";
172	+	break;
173	+	case 2: suffix = "_spl2_mct";
174	+	break;
175	+	case 3: suffix = "_mct";
176	+	break;
177	+	}
178	+	h["deltax"+suffix] = subDir.make<TH1D>(TString("deltax"+suffix), TString("x-residual pvtx"+suffix),800,-0.04,0.04);
179	+	h["deltay"+suffix] = subDir.make<TH1D>(TString("deltay"+suffix), TString("y-residual pvtx"+suffix),800,-0.04,0.04);
180	+	h["deltaz"+suffix] = subDir.make<TH1D>(TString("deltaz"+suffix), TString("z-residual pvtx"+suffix),800,-0.04,0.04);
181	+	h["pullx"+suffix]  = subDir.make<TH1D>(TString("pullx"+suffix), TString("x-pull pvtx"+suffix),800,-10,10);
182	+	h["pully"+suffix]  = subDir.make<TH1D>(TString("pully"+suffix), TString("y-pull pvtx"+suffix),800,-10,10);
183	+	h["pullz"+suffix]  = subDir.make<TH1D>(TString("pullz"+suffix), TString("z-pull pvtx"+suffix),800,-10,10);
184	+	h["errPVx"+suffix] = subDir.make<TH1D>(TString("errPVx"+suffix), TString("X"+suffix+" vertex error"),100,0.,0.02);
185	+	h["errPVy"+suffix] = subDir.make<TH1D>(TString("errPVy"+suffix), TString("Y"+suffix+" vertex error"),100,0.,0.02);
186	+	h["errPVz"+suffix] = subDir.make<TH1D>(TString("errPVz"+suffix), TString("Z"+suffix+" vertex error"),100,0.,0.02);
187	+
188	+	// Book residual, error and pull histograms for each nTrk bin
189		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
190		stringstream ss;
191		ss << ntrk;
192	<	string suffix = ss.str();
193	<	h["resx_"+suffix]  = fs->make<TH1D> (TString("deltax_"+suffix),"x-distance (Rec - Sim)",400,-0.02,0.02);
194	<	h["resx_"+suffix]->GetXaxis()->SetTitle("cm");
195	<	h["resy_"+suffix]  = fs->make<TH1D> (TString("deltay_"+suffix),"y-distance (Rec - Sim)",400,-0.02,0.02);
196	<	h["resy_"+suffix]->GetXaxis()->SetTitle("cm");
197	<	h["resz_"+suffix]  = fs->make<TH1D> (TString("deltaz_"+suffix),"z-distance (Rec - Sim)",400,-0.02,0.02);
198	<	h["resz_"+suffix]->GetXaxis()->SetTitle("cm");
199	<	h["pullx_"+suffix] = fs->make<TH1D> (TString("pullx_"+suffix),"x-pull (Rec - Sim)",200,-10.,10.);
200	<	h["pullx_"+suffix]->GetXaxis()->SetTitle("cm");
201	<	h["pully_"+suffix] = fs->make<TH1D> (TString("pully_"+suffix),"y-pull (Rec - Sim)",200,-10.,10.);
202	<	h["pully_"+suffix]->GetXaxis()->SetTitle("cm");
203	<	h["pullz_"+suffix] = fs->make<TH1D> (TString("pullz_"+suffix),"z-pull (Rec - Sim)",200,-10.,10.);
204	<	h["pullz_"+suffix]->GetXaxis()->SetTitle("cm");
205	<	suffix.clear();
206	<	}
192	>	string suffix2 = suffix+"_"+ss.str();
193	>	h["deltax"+suffix2]  = subDir1.make<TH1D>(TString("deltax"+suffix2), TString("x-residual of pvtx"+suffix),100,-0.02,0.02);
194	>	h["deltax"+suffix2]->GetXaxis()->SetTitle("cm");
195	>	h["deltay"+suffix2]  = subDir1.make<TH1D>(TString("deltay"+suffix2), TString("y-residual of pvtx"+suffix),100,-0.02,0.02);
196	>	h["deltay"+suffix2]->GetXaxis()->SetTitle("cm");
197	>	h["deltaz"+suffix2]  = subDir1.make<TH1D>(TString("deltaz"+suffix2), TString("z-residual of pvtx"+suffix),100,-0.02,0.02);
198	>	h["deltaz"+suffix2]->GetXaxis()->SetTitle("cm");
199	>	h["pullx"+suffix2] = subDir1.make<TH1D>(TString("pullx"+suffix2), TString("x-pull of pvtx"+suffix),100,-10.,10.);
200	>	h["pully"+suffix2] = subDir1.make<TH1D>(TString("pully"+suffix2), TString("y-pull of pvtx"+suffix),100,-10.,10.);
201	>	h["pullz"+suffix2] = subDir1.make<TH1D>(TString("pullz"+suffix2), TString("z-pull of pvtx"+suffix),100,-10.,10.);
202	>	h["errPVx"+suffix2] = subDir1.make<TH1D>(TString("errPVx"+suffix2), TString("X"+suffix+" vertex error"),100,0.,0.02);
203	>	h["errPVy"+suffix2] = subDir1.make<TH1D>(TString("errPVy"+suffix2), TString("Y"+suffix+" vertex error"),100,0.,0.02);
204	>	h["errPVz"+suffix2] = subDir1.make<TH1D>(TString("errPVz"+suffix2), TString("Z"+suffix+" vertex error"),100,0.,0.02);
205	>	suffix2.clear();
206	>	} // End of   for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
207	>
208	>	// Book residual and pull histograms only when analyzeOntheFly is enabled
209	>	if(analyze_) {
210	>	h2["resx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resx"+suffix+"_nTrk"), TString("x-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
211	>	h2["resx"+suffix+"_nTrk"]->SetMarkerStyle(21);
212	>	h2["resx"+suffix+"_nTrk"]->SetMarkerColor(4);
213	>	h2["resx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
214	>	h2["resx"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
215	>	h2["resy"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resy"+suffix+"_nTrk"), TString("y-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
216	>	h2["resy"+suffix+"_nTrk"]->SetMarkerStyle(21);
217	>	h2["resy"+suffix+"_nTrk"]->SetMarkerColor(4);
218	>	h2["resy"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
219	>	h2["resy"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
220	>	h2["resz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resz"+suffix+"_nTrk"), TString("z-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
221	>	h2["resz"+suffix+"_nTrk"]->SetMarkerStyle(21);
222	>	h2["resz"+suffix+"_nTrk"]->SetMarkerColor(4);
223	>	h2["resz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
224	>	h2["resz"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
225	>	h2["pullx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullx"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
226	>	h2["pullx"+suffix+"_nTrk"]->SetMarkerStyle(21);
227	>	h2["pullx"+suffix+"_nTrk"]->SetMarkerColor(4);
228	>	h2["pullx"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
229	>	h2["pullx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
230	>	h2["pully"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pully"+suffix+"_nTrk"), TString("y-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
231	>	h2["pully"+suffix+"_nTrk"]->SetMarkerStyle(21);
232	>	h2["pully"+suffix+"_nTrk"]->SetMarkerColor(4);
233	>	h2["pully"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
234	>	h2["pully"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
235	>	h2["pullz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullz"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
236	>	h2["pullz"+suffix+"_nTrk"]->SetMarkerStyle(21);
237	>	h2["pullz"+suffix+"_nTrk"]->SetMarkerColor(4);
238	>	h2["pullz"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
239	>	h2["pullz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
240	>	} // End of   if(analyze_) {
241	>	suffix.clear();
242	>	} // End of Book histograms sensitive to data/mc
243	>	// Book MC only plots
244	>	if (!realData_) {
245	>	h["genPart_T"]      = subDir.make<TH1D>("genPart_T","t component of gen particles",300,-0.5,0.5);
246	>	h["genPart_T"]->GetXaxis()->SetTitle("t (nanosecond)");
247	>	h["genPart_cT"]     = subDir.make<TH1D>("genPart_cT","ct component of gen particles",300,-150.,150.);
248	>	h["genPart_cT"]->GetXaxis()->SetTitle("ct (mm)");
249	>	h["nsimPV"]         = subDir.make<TH1D>("nsimPV","Num of sim PV",51,-0.5,50.5);
250		}
251		}
252
136	–
253		PVStudy::~PVStudy()
254		{
255
#	Line 142 | Line 258 | PVStudy::~PVStudy()
258
259		}
260
261	<
261	>	void PVStudy::setRootStyle() {
262	>	//
263	>	gROOT->SetStyle("Plain");
264	>	gStyle->SetFillColor(1);
265	>	gStyle->SetOptDate();
266	>	//   gStyle->SetOptStat(1111110);
267	>	//   gStyle->SetOptFit(0111);
268	>	//   gStyle->SetPadTickX(1);
269	>	//   gStyle->SetPadTickY(1);
270	>	gStyle->SetMarkerSize(0.5);
271	>	gStyle->SetMarkerStyle(8);
272	>	gStyle->SetGridStyle(3);
273	>	//gStyle->SetPadGridX(1);
274	>	//gStyle->SetPadGridY(1);
275	>	gStyle->SetPaperSize(TStyle::kA4);
276	>	gStyle->SetStatW(0.25);             // width of statistics box; default is 0.19
277	>	//   gStyle->SetStatH(0.10);             // height of statistics box; default is 0.1
278	>	gStyle->SetStatFormat("6.4g");      // leave default format for now
279	>	gStyle->SetTitleSize(0.055, "");    // size for pad title; default is 0.02
280	>	gStyle->SetLabelSize(0.03, "XYZ");  // size for axis labels; default is 0.04
281	>	gStyle->SetStatFontSize(0.08);      // size for stat. box
282	>	gStyle->SetTitleFont(32, "XYZ");    // times-bold-italic font (p. 153) for axes
283	>	gStyle->SetTitleFont(32, "");       // same for pad title
284	>	gStyle->SetLabelFont(32, "XYZ");    // same for axis labels
285	>	gStyle->SetStatFont(32);            // same for stat. box
286	>	gStyle->SetLabelOffset(0.006, "Y"); // default is 0.005
287	>	//
288	>	return;
289	>	}
290		//
291		// member functions
292		//
#	Line 158 | Line 302 | std::vector<PVStudy::simPrimaryVertex> P
302		cout << "number of vertices " << evt->vertices_size() << endl;
303		}
304
305	<	int idx=0;
305	>	int idx=0; int npv=0;
306		for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
307	<	vitr != evt->vertices_end(); ++vitr )
308	<	{ // loop for vertex ...
309	<	HepMC::FourVector pos = (*vitr)->position();
310	<	//HepLorentzVector pos = (*vitr)->position();
311	<	if (pos.t()>0) { continue;} // ?
307	>	vitr != evt->vertices_end(); ++vitr ) { // loop for vertex ...
308	>	HepMC::FourVector pos = (*vitr)->position();
309	>	//HepLorentzVector pos = (*vitr)->position();
310	>
311	>	// t component of PV:
312	>	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
313	>	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
314	>	//        std::cout << "Status = " << (*mother)->status() << std::endl;
315	>	HepMC::GenVertex * mv=(*mother)->production_vertex();
316	>	if( ((*mother)->status() == 3) && (!mv)) {
317	>	//        std::cout << "npv= " << npv << std::endl;
318	>	if (npv == 0) {
319	>	h["genPart_cT"]->Fill(pos.t()); // mm
320	>	h["genPart_T"]->Fill(pos.t()/299.792458); // ns
321	>	}
322	>	npv++;
323	>	}
324	>	}
325	>	//       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
326
327	<	bool hasMotherVertex=false;
328	<	//      if(verbose_) cout << "mothers" << endl;
329	<	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
330	<	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
331	<	HepMC::GenVertex * mv=(*mother)->production_vertex();
332	<	if (mv) {hasMotherVertex=true;}
333	<	//      if(verbose_) {
334	<	//        cout << "\t";
335	<	//        (*mother)->print();
336	<	//      }
327	>	bool hasMotherVertex=false;
328	>	if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl;
329	>	for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
330	>	mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
331	>	HepMC::GenVertex * mv=(*mother)->production_vertex();
332	>	//        std::cout << "Status = " << (*mother)->status() << std::endl;
333	>	if (mv) {
334	>	hasMotherVertex=true;
335	>	if(!verbose_) break; //if verbose_, print all particles of gen vertices
336	>	}
337	>	if(verbose_) {
338	>	cout << "\t";
339	>	(*mother)->print();
340		}
341	+	}
342
343	<	if(hasMotherVertex) {continue;}
343	>	if(hasMotherVertex) continue;
344
345	<	// could be a new vertex, check  all primaries found so far to avoid multiple entries
346	<	const double mm2cm=0.1;
347	<	simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
348	<	simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
349	<	for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
350	<	v0!=simpv.end(); v0++){
351	<	if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
352	<	vp=&(*v0);
353	<	break;
192	<	}
345	>	// could be a new vertex, check  all primaries found so far to avoid multiple entries
346	>	const double mm2cm=0.1;
347	>	simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
348	>	simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
349	>	for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
350	>	v0!=simpv.end(); v0++){
351	>	if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
352	>	vp=&(*v0);
353	>	break;
354		}
355	+	}
356
357	<	if(!vp){
358	<	// this is a new vertex
359	<	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
360	<	simpv.push_back(sv);
361	<	vp=&simpv.back();
362	<	}else{
363	<	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
364	<	}
365	<	vp->genVertex.push_back((*vitr)->barcode());
366	<	// collect final state descendants
367	<	for ( HepMC::GenVertex::particle_iterator
368	<	daughter  = (*vitr)->particles_begin(HepMC::descendants);
369	<	daughter != (*vitr)->particles_end(HepMC::descendants);
370	<	++daughter ) {
371	<	if (isFinalstateParticle(*daughter)){
372	<	if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
373	<	== vp->finalstateParticles.end()){
374	<	vp->finalstateParticles.push_back((*daughter)->barcode());
375	<	HepMC::FourVector m=(*daughter)->momentum();
376	<	// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
377	<	// but adding FourVectors seems not to be foreseen
378	<	vp->ptot.setPx(vp->ptot.px()+m.px());
379	<	vp->ptot.setPy(vp->ptot.py()+m.py());
380	<	vp->ptot.setPz(vp->ptot.pz()+m.pz());
381	<	vp->ptot.setE(vp->ptot.e()+m.e());
382	<	vp->ptsq+=(m.perp())*(m.perp());
383	<	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
222	<	vp->nGenTrk++;
223	<	}
357	>	if(!vp){
358	>	// this is a new vertex
359	>	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
360	>	simpv.push_back(sv);
361	>	vp=&simpv.back();
362	>	}else{
363	>	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
364	>	}
365	>	vp->genVertex.push_back((*vitr)->barcode());
366	>	// collect final state descendants
367	>	for ( HepMC::GenVertex::particle_iterator daughter  = (*vitr)->particles_begin(HepMC::descendants);
368	>	daughter != (*vitr)->particles_end(HepMC::descendants);
369	>	++daughter ) {
370	>	if (isFinalstateParticle(*daughter)){
371	>	if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
372	>	== vp->finalstateParticles.end()){
373	>	vp->finalstateParticles.push_back((*daughter)->barcode());
374	>	HepMC::FourVector m=(*daughter)->momentum();
375	>	// the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
376	>	// but adding FourVectors seems not to be foreseen
377	>	vp->ptot.setPx(vp->ptot.px()+m.px());
378	>	vp->ptot.setPy(vp->ptot.py()+m.py());
379	>	vp->ptot.setPz(vp->ptot.pz()+m.pz());
380	>	vp->ptot.setE(vp->ptot.e()+m.e());
381	>	vp->ptsq+=(m.perp())*(m.perp());
382	>	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
383	>	vp->nGenTrk++;
384		}
385		}
386		}
227	–	idx++;
387		}
388	+	}
389		}
390		return simpv;
391		}
#	Line 305 | Line 465 | void
465		PVStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
466		{
467		using namespace edm;
468	<
309	<	using reco::TrackCollection;
468	>	using namespace reco;
469
470	<	Handle<TrackCollection> tracks;
471	<	iEvent.getByLabel(tracksLabel_,tracks);
472	<
473	<	Handle<reco::VertexCollection> recVtxs;
474	<	iEvent.getByLabel(vtxSample_, recVtxs);
470	>	if (verbose_)
471	>	cout<<"PVStudy::analyze():"<<endl;
472	>	//=======================================================
473	>	// Initialize Root-tuple variables if needed
474	>	//=======================================================
475	>	//if(saventuple_)
476	>	SetVarToZero();
477	>
478	>	//=======================================================
479	>	// Track accessors
480	>	//=======================================================
481	>	//trackCollection
482	>	static const reco::TrackCollection s_empty_trackColl;
483	>	const reco::TrackCollection *trackColl = &s_empty_trackColl;
484	>	edm::Handle<reco::TrackCollection> trackCollectionHandle;
485	>	iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle);
486	>	if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
487	>	trackColl = trackCollectionHandle.product();
488	>	} else {
489	>	cout << "trackCollection cannot be found -> using empty collection of same type." <<endl;
490	>	}
491	>	//splitTrackCollection1
492	>	static const reco::TrackCollection s_empty_splitTrackColl1;
493	>	const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
494	>	edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
495	>	iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle);
496	>	if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
497	>	splitTrackColl1 = splitTrackCollection1Handle.product();
498	>	} else {
499	>	cout << "splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
500	>	}
501	>	//splitTrackCollection2
502	>	static const reco::TrackCollection s_empty_splitTrackColl2;
503	>	const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
504	>	edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
505	>	iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle);
506	>	if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
507	>	splitTrackColl2 = splitTrackCollection2Handle.product();
508	>	} else {
509	>	cout << "splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
510	>	}
511	>
512	>	//=======================================================
513	>	// PVTX accessors
514	>	//=======================================================
515	>	//vertexCollection
516	>	static const reco::VertexCollection s_empty_vertexColl;
517	>	const reco::VertexCollection *vertexColl = &s_empty_vertexColl;
518	>	edm::Handle<reco::VertexCollection> vertexCollectionHandle;
519	>	iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle);
520	>	if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
521	>	vertexColl = vertexCollectionHandle.product();
522	>	} else {
523	>	cout << "vertexCollection cannot be found -> using empty collection of same type." <<endl;
524	>	}
525	>	//splitVertexCollection1
526	>	static const reco::VertexCollection s_empty_splitVertexColl1;
527	>	const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
528	>	edm::Handle<reco::VertexCollection> splitVertexCollection1Handle;
529	>	iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle);
530	>	if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
531	>	splitVertexColl1 = splitVertexCollection1Handle.product();
532	>	} else {
533	>	cout << "splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
534	>	}
535	>	//splitVertexCollection2
536	>	static const reco::VertexCollection s_empty_splitVertexColl2;
537	>	const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
538	>	edm::Handle<reco::VertexCollection> splitVertexCollection2Handle;
539	>	iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle);
540	>	if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
541	>	splitVertexColl2 = splitVertexCollection2Handle.product();
542	>	} else {
543	>	cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
544	>	}
545
546	+	if(verbose_) cout<<"Done accessing the track and vertex collections"<<endl;
547	+
548	+	//=======================================================
549	+	// MC simvtx accessor
550	+	//=======================================================
551		if (!realData_) {
552	<	Handle<SimVertexContainer> simVtxs;
552	>	edm::Handle<SimVertexContainer> simVtxs;
553		iEvent.getByLabel( simG4_, simVtxs);
554
555	<	Handle<SimTrackContainer> simTrks;
555	>	edm::Handle<SimTrackContainer> simTrks;
556		iEvent.getByLabel( simG4_, simTrks);
557		}
558
559	<	if(verbose_) printRecVtxs(recVtxs);
560	<
559	>	//=======================================================
560	>	// GET PDT
561	>	//=======================================================
562		try{
563		iSetup.getData(pdt);
564		}catch(const Exception&){
565		cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
566		}
567
568	+	//=======================================================
569	+	// Fill trackparameters of the input tracks to pvtx fitter
570	+	//=======================================================
571	+
572	+	// GeneralTracks
573	+	h["nTrk"]->Fill(trackColl->size());
574	+	for(TrackCollection::const_iterator itTrack = trackColl->begin();
575	+	itTrack != trackColl->end();
576	+	++itTrack) {
577	+	h["trkPt"]->Fill(itTrack->pt());
578	+	h["trkDxy"]->Fill(itTrack->dxy());
579	+	h["trkDz"]->Fill(itTrack->dz());
580	+	h["trkEta"]->Fill(itTrack->eta());
581	+	h["trkPhi"]->Fill(itTrack->phi());
582	+	}
583	+
584	+	//SplittedTracks1
585	+	h["nTrk1_spl"]->Fill(splitTrackColl1->size());
586	+	for(TrackCollection::const_iterator itTrack = splitTrackColl1->begin();
587	+	itTrack != splitTrackColl1->end();
588	+	++itTrack) {
589	+	h["trkPt1_spl"]->Fill(itTrack->pt());
590	+	h["trkDxy1_spl"]->Fill(itTrack->dxy());
591	+	h["trkDz1_spl"]->Fill(itTrack->dz());
592	+	h["trkEta1_spl"]->Fill(itTrack->eta());
593	+	h["trkPhi1_spl"]->Fill(itTrack->phi());
594	+	}
595	+	//SplittedTracks2
596	+	h["nTrk2_spl"]->Fill(splitTrackColl2->size());
597	+	for(TrackCollection::const_iterator itTrack = splitTrackColl2->begin();
598	+	itTrack != splitTrackColl2->end();
599	+	++itTrack) {
600	+	h["trkPt2_spl"]->Fill(itTrack->pt());
601	+	h["trkDxy2_spl"]->Fill(itTrack->dxy());
602	+	h["trkDz2_spl"]->Fill(itTrack->dz());
603	+	h["trkEta2_spl"]->Fill(itTrack->eta());
604	+	h["trkPhi2_spl"]->Fill(itTrack->phi());
605	+	}
606	+
607	+	if(verbose_)
608	+	cout<<"Done filling track parameters of the input tracks to pvtx fitter."<<endl;
609	+
610	+	//=======================================================
611	+	// Fill reconstructed vertices
612	+	//=======================================================
613	+	if(verbose_)  {
614	+	cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
615	+	printRecVtxs(vertexCollectionHandle);
616	+	cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
617	+	printRecVtxs(splitVertexCollection1Handle);
618	+	cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
619	+	printRecVtxs(splitVertexCollection2Handle);
620	+	}
621	+
622	+	//Fill the tree variables
623	+	nrecPV_ = vertexColl->size();
624	+	nrecPV1_spl_ = splitVertexColl1->size();
625	+	nrecPV2_spl_ = splitVertexColl2->size();
626	+
627	+	h["nrecPV"]->Fill(nrecPV_);
628	+	h["nrecPV1_spl"]->Fill(nrecPV1_spl_);
629	+	h["nrecPV2_spl"]->Fill(nrecPV2_spl_);
630	+	h["nrecPVDiff"]->Fill(double(nrecPV1_spl_)-double(nrecPV2_spl_));
631	+
632	+	// Fill the track paramters for the vertexColl
633	+	for(reco::VertexCollection::const_iterator v=vertexColl->begin();
634	+	v!=vertexColl->end(); ++v){
635	+	h["nTrkPV"]->Fill(v->tracksSize());
636	+	try {
637	+	for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
638	+	t!=v->tracks_end(); t++) {
639	+	// illegal charge
640	+	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
641	+	//      h["trkPtPV"]->Fill(0.);
642	+	}
643	+	else {
644	+	h["trkPtPV"]->Fill((**t).pt());
645	+	h["trkDxyPV"]->Fill((**t).dxy());
646	+	h["trkDzPV"]->Fill((**t).dz());
647	+	h["trkEtaPV"]->Fill((**t).eta());
648	+	h["trkPhiPV"]->Fill((**t).phi());
649	+	}
650	+	}
651	+	}
652	+	catch (...) {
653	+	// exception thrown when trying to use linked track
654	+	//          h["trkPtPV"]->Fill(0.);
655	+	}
656	+	}
657	+
658	+
659	+	// == ignore multi-vertices (for now)
660	+	//take single pvt collection and compare
661	+	if(splitVertexColl1->size() == 1 && splitVertexColl2->size() == 1 ) {
662	+	const reco::Vertex recvtx1 = *(splitVertexColl1->begin());
663	+	const reco::Vertex recvtx2 = *(splitVertexColl2->begin());
664	+	if (matchVertex(recvtx1, recvtx2))  {
665	+	int nTrkPV1 = recvtx1.tracksSize();
666	+	int nTrkPV2 = recvtx2.tracksSize();
667	+	h["nTrkPV1_spl"]->Fill(nTrkPV1);
668	+	h["nTrkPV2_spl"]->Fill(nTrkPV2);
669	+	double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
670	+	h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
671	+	h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
672	+
673	+	// Fill the track paramters for the splitVertexColl1
674	+	try {
675	+	for(reco::Vertex::trackRef_iterator t = recvtx1.tracks_begin();
676	+	t!=recvtx1.tracks_end(); t++) {
677	+	// illegal charge
678	+	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
679	+	//    h["trkPtPV"]->Fill(0.);
680	+	}
681	+	else {
682	+	h["trkPtPV1_spl"]->Fill((**t).pt());
683	+	h["trkDxyPV1_spl"]->Fill((**t).dxy());
684	+	h["trkDzPV1_spl"]->Fill((**t).dz());
685	+	h["trkEtaPV1_spl"]->Fill((**t).eta());
686	+	h["trkPhiPV1_spl"]->Fill((**t).phi());
687	+	}
688	+	}
689	+	}
690	+	catch (...) {
691	+	// exception thrown when trying to use linked track
692	+	//        h["trkPtPV"]->Fill(0.);
693	+	}
694	+
695	+	// Fill the track paramters for the splitVertexColl2
696	+	try {
697	+	for(reco::Vertex::trackRef_iterator t = recvtx2.tracks_begin();
698	+	t!=recvtx2.tracks_end(); t++) {
699	+	// illegal charge
700	+	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
701	+	//    h["trkPtPV"]->Fill(0.);
702	+	}
703	+	else {
704	+	h["trkPtPV2_spl"]->Fill((**t).pt());
705	+	h["trkDxyPV2_spl"]->Fill((**t).dxy());
706	+	h["trkDzPV2_spl"]->Fill((**t).dz());
707	+	h["trkEtaPV2_spl"]->Fill((**t).eta());
708	+	h["trkPhiPV2_spl"]->Fill((**t).phi());
709	+	}
710	+	}
711	+	}
712	+	catch (...) {
713	+	// exception thrown when trying to use linked track
714	+	//        h["trkPtPV"]->Fill(0.);
715	+	}
716	+
717	+	if(verbose_)
718	+	cout<<"Done filling track parameters in reconstruced vertices."<<endl;
719	+
720	+	if(fabs(ntrkreldiff)<ntrkdiffcut_) {
721	+	fres_[0] = (recvtx1.x()-recvtx2.x())/sqrt(2.0);
722	+	fres_[1] = (recvtx1.y()-recvtx2.y())/sqrt(2.0);
723	+	fres_[2] = (recvtx1.z()-recvtx2.z())/sqrt(2.0);
724	+	ferror_[0] = sqrt(pow(recvtx1.xError(),2)+pow(recvtx2.xError(),2))/sqrt(2);
725	+	ferror_[1] = sqrt(pow(recvtx1.yError(),2)+pow(recvtx2.yError(),2))/sqrt(2);
726	+	ferror_[2] = sqrt(pow(recvtx1.zError(),2)+pow(recvtx2.zError(),2))/sqrt(2);
727	+	fntrk_ =  (nTrkPV1+nTrkPV2)/2;
728	+
729	+	h["deltax"]->Fill( fres_[0] );
730	+	h["deltay"]->Fill( fres_[1] );
731	+	h["deltaz"]->Fill( fres_[2] );
732	+	h["pullx"]->Fill( fres_[0]/ferror_[0]);
733	+	h["pully"]->Fill( fres_[1]/ferror_[1]);
734	+	h["pullz"]->Fill( fres_[2]/ferror_[2]);
735	+	h["errPVx"]->Fill( ferror_[0] );
736	+	h["errPVy"]->Fill( ferror_[1] );
737	+	h["errPVz"]->Fill( ferror_[2] );
738	+	pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
739	+	error(ferror_[0], ferror_[1], ferror_[2]),
740	+	fntrk_));
741	+	// Fill histo according to its track multiplicity, set datamode = 0
742	+	fillHisto(res(fres_[0], fres_[1], fres_[2]),
743	+	error(ferror_[0], ferror_[1], ferror_[2]),
744	+	fntrk_,0);
745	+
746	+	if (saventuple_) ftree_->Fill();
747	+	} // End of  if(fabs(ntrkreldiff)<ntrkdiffcut_) {
748	+	} // End of  if (matchVertex(recvtx1, recvtx2))  {
749	+	else   {
750	+	cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
751	+	cout<<"recvtx1.z() = " << recvtx1.z() <<"; recvtx2.z() = " <<  recvtx2.z() <<endl;
752	+	}
753	+	} // End of fill histograms for two-vertices Method
754	+
755	+	if(verbose_)
756	+	cout<<"Done filling res/pull with two-vertices method"<<endl;
757	+
758	+	//=======================================================
759	+	// Fill simulated vertices
760	+	//=======================================================
761		if (!realData_) {
762		bool MC=false;
763		Handle<HepMCProduct> evtMC;
#	Line 345 | Line 773 | PVStudy::analyze(const edm::Event& iEven
773		}
774		MC=true;
775		}
776	<
776	>
777		if(MC){
778		// make a list of primary vertices:
779		std::vector<simPrimaryVertex> simpv;
780		simpv=getSimPVs(evtMC,"");
781		//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
782	<
782	>	h["nsimPV"]->Fill(simpv.size());
783		int nsimtrk=0;
784		int nrectrk=0;
785	+
786		for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
787		vsim!=simpv.end(); vsim++){
788		nsimtrk+=vsim->nGenTrk;
789	<	// look for a matching reconstructed vertex
790	<	vsim->recVtx=NULL;
791	<	for(reco::VertexCollection::const_iterator vrec=recVtxs->begin();
792	<	vrec!=recVtxs->end(); ++vrec){
793	<	nrectrk=vrec->tracksSize();
794	<	if(verbose_){
795	<	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
796	<	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
797	<	}
798	<	if ( matchVertex(*vsim,*vrec) ){
799	<	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
800	<	if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
801	<	|| (!vsim->recVtx) )
802	<	{
789	>	// Loop over the datamodes, for MC, fill only 1, 2, 3
790	>	for (int i=1;i<=3;i++) {
791	>	static const reco::VertexCollection s_empty_vtxColl;
792	>	const reco::VertexCollection *vtxColl = &s_empty_vtxColl;
793	>	std::string suffix;
794	>	// Set the vertexColl and histogram suffix according to datamode
795	>	if (i == 1) {
796	>	vtxColl = splitVertexColl1;
797	>	suffix = "_spl1_mct";
798	>	}
799	>	if (i == 2) {
800	>	vtxColl =  splitVertexColl2;
801	>	suffix = "_spl2_mct";
802	>	}
803	>	if (i == 3) {
804	>	vtxColl =  vertexColl;
805	>	suffix = "_mct";
806	>	}
807	>	//========================================================
808	>	//  look for a matching reconstructed vertex in vertexColl
809	>	//========================================================
810	>	vsim->recVtx=NULL;
811	>	for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
812	>	vrec!=vtxColl->end(); ++vrec){
813	>	nrectrk=vrec->tracksSize();
814	>	if(verbose_){
815	>	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
816	>	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
817	>	}
818	>	if ( matchVertex(*vsim,*vrec) ){
819	>	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
820	>	if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
821	>	|| (!vsim->recVtx) )
822		vsim->recVtx=&(*vrec);
823	+	}// End of finding the matched reco_vertex
824	+
825	+	if (vsim->recVtx){
826	+	if(verbose_) {
827	+	cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
828	+	}
829	+	double fres_mct[3];
830	+	double ferror_mct[3];
831	+	int fntrk_mct;
832	+
833	+	fres_mct[0] = vsim->recVtx->x()-vsim->x;
834	+	fres_mct[1] = vsim->recVtx->y()-vsim->y;
835	+	fres_mct[2] = vsim->recVtx->z()-vsim->z;
836	+	ferror_mct[0] = vsim->recVtx->xError();
837	+	ferror_mct[1] = vsim->recVtx->yError();
838	+	ferror_mct[2] = vsim->recVtx->zError();
839	+	fntrk_mct = nrectrk;
840	+
841	+	//Fill the values for variables in ftree_
842	+	if(i == 1) {
843	+	for(int j = 0;j<3;j++)
844	+	fres_spl1_mct_[j] = fres_mct[j];
845	+	fntrk_spl1_mct_ =  fntrk_mct;
846		}
847	+	if(i == 2) {
848	+	for(int j = 0;j<3;j++)
849	+	fres_spl2_mct_[j] = fres_mct[j];
850	+	fntrk_spl2_mct_ =  fntrk_mct;
851	+	}
852	+	if(i == 3) {
853	+	for(int j=0;j<3;j++)
854	+	fres_mct_[j] = fres_mct[j];
855	+	fntrk_mct_ =  fntrk_mct;
856	+	}
857	+
858	+	h["deltax"+suffix]->Fill( fres_mct[0] );
859	+	h["deltay"+suffix]->Fill( fres_mct[1] );
860	+	h["deltaz"+suffix]->Fill( fres_mct[2] );
861	+	h["pullx"+suffix]->Fill( fres_mct[0]/ferror_mct[0] );
862	+	h["pully"+suffix]->Fill( fres_mct[1]/ferror_mct[1] );
863	+	h["pullz"+suffix]->Fill( fres_mct[2]/ferror_mct[2] );
864	+	h["errPVx"+suffix]->Fill( ferror_mct[0] );
865	+	h["errPVy"+suffix]->Fill( ferror_mct[1] );
866	+	h["errPVz"+suffix]->Fill( ferror_mct[2] );
867	+	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
868	+	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
869	+	fntrk_mct));
870	+	// Fill histo according to its track multiplicity
871	+	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
872	+	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
873	+	fntrk_mct, i);
874	+	if (saventuple_) ftree_->Fill();
875	+	suffix.clear();
876	+	}
877	+	else {  // no rec vertex found for this simvertex
878	+	if(verbose_)
879	+	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
880	+	}
881		}
882	<	}
883	<	if (vsim->recVtx){
884	<	if(verbose_) {
885	<	cout <<"primary matched " << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
381	<	}
382	<	h_deltax->Fill( vsim->recVtx->x()-vsim->x );
383	<	h_deltay->Fill( vsim->recVtx->y()-vsim->y );
384	<	h_deltaz->Fill( vsim->recVtx->z()-vsim->z );
385	<	h_pullx->Fill( (vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError() );
386	<	h_pully->Fill( (vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError() );
387	<	h_pullz->Fill( (vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError() );
388	<	pvanainfo.push_back(PVStudy::PVAnaInfo( vsim->recVtx->x()-vsim->x,
389	<	vsim->recVtx->y()-vsim->y,
390	<	vsim->recVtx->z()-vsim->z,
391	<	(vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError(),
392	<	(vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError(),
393	<	(vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError(),
394	<	nrectrk )
395	<	);
396	<	}
397	<	else{  // no rec vertex found for this simvertex
398	<	if(verbose_) {
399	<	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
400	<	}
401	<	}
402	<	}
403	<	}
404	<	}// With MC
405	<
406	<	// Loop RecVtxs and find matched SimVtxs
407	<	for(reco::VertexCollection::const_iterator v=recVtxs->begin();
408	<	v!=recVtxs->end(); ++v){
409	<	try {
410	<	for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
411	<	t!=v->tracks_end(); t++) {
412	<	// illegal charge
413	<	if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
414	<	//      h_trkPtPV->Fill(0.);
415	<	}
416	<	else {
417	<	h_trkPtPV->Fill((**t).pt());
418	<	}
419	<	}
420	<	}
421	<	catch (...) {
422	<	// exception thrown when trying to use linked track
423	<	//          h_trkPtPV->Fill(0.);
424	<	}
425	<
426	<	h_nTrkPV->Fill(v->tracksSize());
427	<
428	<	}
429	<
430	<	h_nTrk->Fill(tracks->size());
431	<
432	<	for(TrackCollection::const_iterator itTrack = tracks->begin();
433	<	itTrack != tracks->end();
434	<	++itTrack) {
435	<	int charge = 0;
436	<	charge = itTrack->charge();
437	<	h_trkPt->Fill(itTrack->pt());
438	<	}
439	<
882	>	} // === End of Looping through vertexColl
883	>	}
884	>	} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
885	>	} // End of Filling MC variables
886		}
887
442	–
888		// ------------ method called once each job just before starting event loop  ------------
889		void
890		PVStudy::beginJob()
#	Line 449 | Line 894 | PVStudy::beginJob()
894		// ------------ method called once each job just after ending the event loop  ------------
895		void
896		PVStudy::endJob() {
897	<	double sqt2 = sqrt(2.);
898	<	if (!realData_) {
899	<	if (verbose_) cout << "Analyzing PV info" << endl;
900	<	for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
901	<	if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
902	<	PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk);
903	<	if (analyze_) {
904	<	if ( ipv.resx > 0 ) h_resx_nTrk->Fill(ntrk,ipv.resx/sqt2, 1.);
905	<	if ( ipv.resy > 0 ) h_resy_nTrk->Fill(ntrk,ipv.resy/sqt2, 1.);
906	<	if ( ipv.resz > 0 ) h_resz_nTrk->Fill(ntrk,ipv.resz/sqt2, 1.);
907	<	if ( ipv.pullx > 0 ) h_pullx_nTrk->Fill(ntrk,ipv.pullx, 1.);
908	<	if ( ipv.pully > 0 ) h_pully_nTrk->Fill(ntrk,ipv.pully, 1.);
909	<	if ( ipv.pullz > 0 ) h_pullz_nTrk->Fill(ntrk,ipv.pullz, 1.);
897	>	if (verbose_) cout << "Analyzing PV info" << endl;
898	>	// Looping through the datamodes available
899	>	for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
900	>	string suffix;
901	>	if(verbose_)  cout<<"datamode = "<< *it<<endl;
902	>	switch(*it) {
903	>	case 0: suffix = "";
904	>	break;
905	>	case 1: suffix = "_spl1_mct";
906	>	break;
907	>	case 2: suffix = "_spl2_mct";
908	>	break;
909	>	case 3: suffix = "_mct";
910	>	break;
911	>	}
912	>	if (analyze_) {
913	>	for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
914	>	PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
915	>	if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
916	>	if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
917	>	if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
918	>	if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
919	>	if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
920	>	if ( ipv.pull_.y() > 0 ) h2["pully"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
921	>	if ( ipv.pull_.z() > 0 ) h2["pullz"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.z(), 1.);
922		}
923		}
924	+	suffix.clear();
925	+	}
926	+	if(saventuple_) {
927	+	file_->cd();
928	+	ftree_->Write();
929	+	file_->Close();
930		}
931		}
932
933	+
934	+
935	+	// Match a recovertex with a simvertex
936	+	// ? Should the cut be parameter dependent?
937	+	// cut can be within n sigma of the vertex.zerror()
938		bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
939		const reco::Vertex & vrec)
940		{
941		return (fabs(vsim.z-vrec.z())<0.0500); // =500um
942		}
943
944	+	// Match two reconstructed vertices
945	+	bool PVStudy::matchVertex( const reco::Vertex & vrec1,
946	+	const reco::Vertex & vrec2)
947	+	{
948	+	return (fabs(vrec1.z()-vrec2.z())<0.0500); // =500um
949	+	}
950	+
951	+
952		bool PVStudy::isResonance(const HepMC::GenParticle * p){
953		double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
954		if(verbose_) cout << "isResonance   " << p->pdg_id() << " " << ctau << endl;
#	Line 544 | Line 1020 | void PVStudy::printSimTrks(const Handle<
1020		}
1021		}
1022
1023	<	PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk) {
1023	>	void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1024	>	stringstream ss;
1025	>	ss << ntk;
1026	>	string suffix = ss.str();
1027	>	if(datamode == 0 )  suffix = "_" + suffix;
1028	>	if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1029	>	if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1030	>	if(datamode == 3 )  suffix = "_mct_" + suffix;
1031	>
1032	>	if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
1033	>	// Fill histograms of res and pull of ntk
1034	>	h["deltax"+suffix]->Fill(r.x());
1035	>	h["deltay"+suffix]->Fill(r.y());
1036	>	h["deltaz"+suffix]->Fill(r.z());
1037	>	h["pullx"+suffix]->Fill(r.x()/er.x());
1038	>	h["pully"+suffix]->Fill(r.y()/er.y());
1039	>	h["pullz"+suffix]->Fill(r.z()/er.z());
1040	>	h["errPVx"+suffix]->Fill(er.x());
1041	>	h["errPVy"+suffix]->Fill(er.y());
1042	>	h["errPVz"+suffix]->Fill(er.z());
1043	>	}
1044	>
1045	>	PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1046		map<int,double> data;
1047		data.clear();
1048	<	double xmin=0.,xmax=0.;
551	<	double fpar[2];
1048	>	double fpar[2] = {-999,-999};
1049		double cm2um = 10000;
1050		stringstream ss;
1051		ss << ntk;
1052	<	string suffix = ss.str();
1053	<
1054	<	for ( int i = 0; i < 6; ++i) {
1055	<	switch(i){
1056	<	case 0:
1057	<	case 1:
1058	<	case 2:
1059	<	if (verbose_) cout << "Analyzing Resolution:" << endl;
1060	<	xmin = -0.01; xmax = 0.01;
1061	<	break;
1062	<	default:
1063	<	if (verbose_) cout << "Analyzing Pull:" << endl;
1064	<	xmin = -10.; xmax = 10;
1065	<	break;
1066	<	}
1067	<	TH1D *hdist = new TH1D("hdist","distribution to be fitted", 1000,xmin,xmax);
1068	<
1069	<	for (std::vector<PVStudy::PVAnaInfo>::const_iterator ipvinfo = pvanainfo.begin();
1070	<	ipvinfo != pvanainfo.end(); ++ipvinfo){
1071	<	if ( verbose_ && i==0 ) {
1072	<	cout << "Num of tracks of the PV = " << ipvinfo->nTrk;
1073	<	cout << "; delta x = " << ipvinfo->resx << endl;
1074	<	}
1075	<	if (ipvinfo->nTrk == ntk) {
1076	<	switch(i) {
1077	<	case 0: if (verbose_) cout << "Matched resx" << endl;
1078	<	hdist->Fill(ipvinfo->resx);
1079	<	h["resx_"+suffix]->Fill(ipvinfo->resx);
1080	<	break;
1081	<	case 1: if (verbose_) cout << "Matched resy" << endl;
1082	<	hdist->Fill(ipvinfo->resy);
1083	<	h["resy_"+suffix]->Fill(ipvinfo->resy);
1084	<	break;
588	<	case 2: if (verbose_) cout << "Matched resz" << endl;
589	<	hdist->Fill(ipvinfo->resz);
590	<	h["resz_"+suffix]->Fill(ipvinfo->resz);
591	<	break;
592	<	case 3: if (verbose_) cout << "Matched pullx" << endl;
593	<	hdist->Fill(ipvinfo->pullx);
594	<	h["pullx_"+suffix]->Fill(ipvinfo->pullx);
595	<	break;
596	<	case 4: if (verbose_) cout << "Matched pully" << endl;
597	<	hdist->Fill(ipvinfo->pully);
598	<	h["pully_"+suffix]->Fill(ipvinfo->pully);
599	<	break;
600	<	case 5: if (verbose_) cout << "Matched pullz" << endl;
601	<	hdist->Fill(ipvinfo->pullz);
602	<	h["pullz_"+suffix]->Fill(ipvinfo->pullz);
603	<	break;
604	<	}
605	<	}
606	<	else {
607	<	//      switch(i) {
608	<	//      case 0: if (verbose_) cout << "bad resx" << endl; break;
609	<	//      case 1: if (verbose_) cout << "bad resy" << endl; break;
610	<	//      case 2: if (verbose_) cout << "bad resz" << endl; break;
611	<	//      case 3: if (verbose_) cout << "bad pullx" << endl; break;
612	<	//      case 4: if (verbose_) cout << "bad pully" << endl; break;
613	<	//      case 5: if (verbose_) cout << "bad pullz" << endl; break;
614	<	//      }
1052	>	string suffix = ss.str();
1053	>	if(datamode == 0 )  suffix = "_" + suffix;
1054	>	if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1055	>	if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1056	>	if(datamode == 3 )  suffix = "_mct_" + suffix;
1057	>
1058	>	if(analyze_) {
1059	>	for ( int i = 0; i < 6; ++i) {
1060	>	switch (i) {
1061	>	case 0:
1062	>	fit(h["deltax"+suffix], fpar);
1063	>	data[i] = fpar[0]*cm2um;
1064	>	break;
1065	>	case 1:
1066	>	fit(h["deltay"+suffix], fpar);
1067	>	data[i] = fpar[0]*cm2um;
1068	>	break;
1069	>	case 2:
1070	>	fit(h["deltaz"+suffix], fpar);
1071	>	data[i] = fpar[0]*cm2um;
1072	>	break;
1073	>	case 3:
1074	>	fit(h["pullx"+suffix], fpar);
1075	>	data[i] = fpar[0];
1076	>	break;
1077	>	case 4:
1078	>	fit(h["pully"+suffix], fpar);
1079	>	data[i] = fpar[0];
1080	>	break;
1081	>	case 5:
1082	>	fit(h["pullz"+suffix], fpar);
1083	>	data[i] = fpar[0];
1084	>	break;
1085		}
1086	+	if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1087		}
1088	<	if(analyze_) {
1089	<	TAxis *axis0 = hdist->GetXaxis();
1090	<	int nbin = axis0->GetLast();
1091	<	double nOF = hdist->GetBinContent(nbin+1);
1092	<	double nUF = hdist->GetBinContent(0);
1093	<	if ( verbose_ && i==0 ){
1094	<	cout << "Last bin = " << nbin;
1095	<	cout << "; hdist: Overflow Entries = " << nOF;
1096	<	cout << "; Underflow Entries = " << nUF;
1097	<	cout << "; hdist->GetEntries() = " << hdist->GetEntries() << endl;
1098	<	}
1099	<	if (hdist->GetEntries() - nOF - nUF >= 5) { // FIXME: for reasonable Gaussian fit
1100	<	if(i<3){
1101	<	hdist->Fit("gaus","Q0");
1102	<	TF1 *fgaus = hdist->GetFunction("gaus");
1103	<	if (verbose_) cout << "got function" << endl;
1104	<	fpar[0] = ((fgaus->GetParameter(1))?fgaus->GetParameter(1):-999);
1105	<	fpar[1] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1106	<	delete fgaus;
1107	<	}
1108	<	else {
1109	<	hdist->Fit("gausn","Q0");
1110	<	TF1 *fgaus = hdist->GetFunction("gausn");
1111	<	if (verbose_) cout << "got function" << endl;
1112	<	fpar[0] = ((fgaus->GetParameter(1))?fgaus->GetParameter(1):-999);
1113	<	fpar[1] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1114	<	delete fgaus;
1115	<	}
1116	<	if (verbose_) cout << "fpar[2] = (" << fpar[0] << "," << fpar[1] << ")" << endl;
1117	<	switch (i) {
1118	<	case 0:
1119	<	h["resx_"+suffix]->Fit("gaus","Q");
1120	<	data[i] = fpar[1]*cm2um;
1121	<	break;
651	<	case 1:
652	<	h["resy_"+suffix]->Fit("gaus","Q");
653	<	data[i] = fpar[1]*cm2um;
654	<	break;
655	<	case 2:
656	<	h["resz_"+suffix]->Fit("gaus","Q");
657	<	data[i] = fpar[1]*cm2um;
658	<	break;
659	<	case 3:
660	<	h["pullx_"+suffix]->Fit("gausn","Q");
661	<	data[i] = fpar[1];
662	<	break;
663	<	case 4:
664	<	h["pully_"+suffix]->Fit("gausn","Q");
665	<	data[i] = fpar[1];
666	<	break;
667	<	case 5:
668	<	h["pullz_"+suffix]->Fit("gausn","Q");
669	<	data[i] = fpar[1];
670	<	break;
671	<	}
672	<	if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" cm") << endl;
1088	>	}
1089	>	else
1090	>	for ( int i = 0; i < 6; ++i) {
1091	>	data[i] = -999;
1092	>	}
1093	>
1094	>	return PVStudy::PVAnaInfo(res(data[0], data[1], data[2]),
1095	>	error(0.,0.,0.),
1096	>	pull(data[3], data[4], data[5]),
1097	>	error(0.,0.,0.),
1098	>	ntk);
1099	>	}
1100	>
1101	>	void PVStudy::fit(TH1 *&hdist, double fitPars[]){
1102	>	TAxis *axis0 = hdist->GetXaxis();
1103	>	int nbin = axis0->GetLast();
1104	>	double nOF = hdist->GetBinContent(nbin+1);
1105	>	double nUF = hdist->GetBinContent(0);
1106	>	//   double fitRange = axis0->GetBinUpEdge(nbin);
1107	>	double fitRange = axis0->GetXmax();
1108	>	double sigMax[2] = {0.,0.};
1109	>
1110	>	if ( verbose_ ){
1111	>	cout << "Last bin = " << nbin;
1112	>	cout << "; hdist: Overflow Entries = " << nOF;
1113	>	cout << "; Underflow Entries = " << nUF;
1114	>	cout << "; hdist->GetEntries() = " << hdist->GetEntries();
1115	>	cout << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1116	>	}
1117	>	if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit
1118	>	for (int bi = 0; bi < nbin; bi++) {
1119	>	if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) {
1120	>	sigMax[0] = axis0->GetBinLowEdge(bi);
1121	>	if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1122		}
1123	<	else {
1124	<	if ( verbose_) cout << "Empty Histo! No Fit!!!" << endl;
1125	<	data[i] = -999;
1123	>	if ( (axis0->GetBinLowEdge(bi) >= 0) && (hdist->GetBinContent(bi) > 0) ) {
1124	>	sigMax[0] = axis0->GetBinUpEdge(bi);
1125	>	if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1126		}
1127		}
1128	<	else
1129	<	data[i] = -999;
1130	<	delete hdist;
1128	>	if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1129	>
1130	>	TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1131	>	fgaus->SetParameter(1, 0.);
1132	>	fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
1133	>	fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
1134	>	fgaus->SetLineColor(4);
1135	>	hdist->Fit(fgaus,"Q");
1136	>	gPad->Update();
1137	>	TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
1138	>	s->SetOptStat(1111111);
1139	>	s->SetOptFit(0111);
1140	>	gPad->Update();
1141	>	if (verbose_) cout << "got function" << endl;
1142	>	fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1143		}
1144	+	else
1145	+	fitPars[0] = -999;
1146	+	}
1147
1148	<	return PVStudy::PVAnaInfo(data[0], data[1], data[2],
1149	<	data[3], data[4], data[5],
1150	<	ntk);
1148	>	void PVStudy::SetVarToZero() {
1149	>
1150	>	// number of reconstructed vertices
1151	>	nrecPV_ = 0;
1152	>	nrecPV1_spl_ = 0;
1153	>	nrecPV2_spl_ = 0;
1154	>	// number of tracks in the vertex
1155	>	fntrk_ = 0;
1156	>	fntrk_spl1_mct_ = 0;
1157	>	fntrk_spl2_mct_ = 0;
1158	>	fntrk_mct_ = 0;
1159	>	//pvtx position (x,y,z) residual and error
1160	>	for(int i = 0; i<3;i++)
1161	>	{
1162	>	fres_[i] = 0;
1163	>	ferror_[i] = 0;
1164	>	fres_spl1_mct_[i] = 0;
1165	>	ferror_spl1_mct_[i] = 0;
1166	>	fres_spl2_mct_[i] = 0;
1167	>	ferror_spl2_mct_[i] = 0;
1168	>	fres_mct_[i] = 0;
1169	>	ferror_mct_[i] = 0;
1170	>	}
1171		}
1172
1173	<	//define this as a plug-in
690	<	DEFINE_FWK_MODULE(PVStudy);
1173	>

Diff Legend

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