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.10 by yygao, Sun Nov 8 02:56:57 2009 UTC vs.
Revision 1.17 by yygao, Tue Jan 26 12:06:49 2010 UTC

#	Line 25 | Line 25 | Implementation:
25		#include <vector>
26		#include <iostream>
27		#include <sstream>
28	–	#include <algorithm>
28
29		// user include files
30		#include "FWCore/Framework/interface/Frameworkfwd.h"
#	Line 50 | Line 49 | Implementation:
49		#include <SimDataFormats/Vertex/interface/SimVertexContainer.h>
50		#include <SimDataFormats/Track/interface/SimTrack.h>
51		#include <SimDataFormats/Track/interface/SimTrackContainer.h>
52	+	// BeamSpot
53	+	#include "DataFormats/BeamSpot/interface/BeamSpot.h"
54	+	// Trigger
55	+	#include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
56	+	#include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
57	+	#include "CondFormats/L1TObjects/interface/L1GtPrescaleFactors.h"
58	+	#include "CondFormats/DataRecord/interface/L1GtPrescaleFactorsAlgoTrigRcd.h"
59	+	#include "CondFormats/DataRecord/interface/L1GtPrescaleFactorsTechTrigRcd.h"
60	+	#include "CondFormats/L1TObjects/interface/L1GtTriggerMask.h"
61	+	#include "CondFormats/DataRecord/interface/L1GtTriggerMaskAlgoTrigRcd.h"
62	+	#include "CondFormats/DataRecord/interface/L1GtTriggerMaskTechTrigRcd.h"
63	+	#include "CondFormats/DataRecord/interface/L1GtTriggerMaskVetoAlgoTrigRcd.h"
64	+	#include "CondFormats/DataRecord/interface/L1GtTriggerMaskVetoTechTrigRcd.h"
65	+	#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
66	+	#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
67	+	#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
68	+
69
70		//root
71		#include <TROOT.h>
#	Line 60 | Line 76 | Implementation:
76		#include <TPad.h>
77
78		using namespace std;
79	+	typedef math::XYZTLorentzVectorF LorentzVector;
80	+	typedef math::XYZPoint Point;
81
82		PVStudy::PVStudy(const edm::ParameterSet& iConfig)
83		{
84		//=======================================================================
85		// Get configuration for TrackTupleMaker
86		//=======================================================================
87	<	simG4_           = iConfig.getParameter<edm::InputTag>( "simG4" );
88	<	trackCollectionTag_     = iConfig.getUntrackedParameter<edm::InputTag>("trackCollection");
89	<	splitTrackCollection1Tag_     = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection1");
90	<	splitTrackCollection2Tag_     = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
87	>	simG4_                     = iConfig.getParameter<edm::InputTag>( "simG4" );
88	>	trackCollectionTag_        = iConfig.getUntrackedParameter<edm::InputTag>("trackCollection");
89	>	splitTrackCollection1Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection1");
90	>	splitTrackCollection2Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
91		vertexCollectionTag_       = iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection");
92	<	splitVertexCollection1Tag_       = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
93	<	splitVertexCollection2Tag_       = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
94	<	pixelVertexCollectionTag_      = iConfig.getParameter<edm::InputTag>("pixelVertexCollectionTag");
95	<	verbose_         = iConfig.getUntrackedParameter<bool>("verbose",false);
96	<	realData_        = iConfig.getUntrackedParameter<bool>("realData",false);
97	<	analyze_         = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
98	<	saventuple_         = iConfig.getUntrackedParameter<bool>("saventuple",false);
99	<	outputfilename_  = iConfig.getUntrackedParameter<string>("OutputFileName");
100	<	ntrkdiffcut_         = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
101	<	zsigncut_         = iConfig.getUntrackedParameter<int>("zsigncut");
102	<	nTrkMin_ =  iConfig.getUntrackedParameter<int>("nTrkMin");
103	<	nTrkMax_ =  iConfig.getUntrackedParameter<int>("nTrkMax");
104	<
92	>	splitVertexCollection1Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
93	>	splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
94	>	verbose_                   = iConfig.getUntrackedParameter<bool>("verbose",false);
95	>	realData_                  = iConfig.getUntrackedParameter<bool>("realData",false);
96	>	analyze_                   = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
97	>	saventuple_                = iConfig.getUntrackedParameter<bool>("saventuple",false);
98	>	outputfilename_            = iConfig.getUntrackedParameter<string>("OutputFileName");
99	>	ntrkdiffcut_               = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
100	>	nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
101	>	nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
102	>	zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
103	>	bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
104	>	applyEvtClean_            = iConfig.getUntrackedParameter<bool>("applyEvtClean",false);
105	>	histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
106	>
107		//now do what ever initialization is needed
108		pvinfo.clear();
109
110	<
110	>
111		// Specify the data mode vector
112		if(realData_) datamode.push_back(0);
113		else {
#	Line 96 | Line 116 | PVStudy::PVStudy(const edm::ParameterSet
116		datamode.push_back(2);
117		datamode.push_back(3);
118		}
119	<
100	<	// Create ntuple files if needed
119	>	// Create ntuple files if needed
120		if(saventuple_) {
121		file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
122		ftree_ = new TTree("mytree","mytree");
#	Line 120 | Line 139 | PVStudy::PVStudy(const edm::ParameterSet
139		pvtxtree_->Branch("recx_err",&recx_err_,"recx_err[nrecPV]/D");
140		pvtxtree_->Branch("recy_err",&recy_err_,"recy_err[nrecPV]/D");
141		pvtxtree_->Branch("recz_err",&recz_err_,"recz_err[nrecPV]/D");
142	<
142	>
143		pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
144		pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
145		pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
#	Line 150 | Line 169 | PVStudy::PVStudy(const edm::ParameterSet
169		pvtxtree_->Branch("recx2_err_spl",&recx2_err_spl_,"recx2_err_spl[nrecPV2_spl]/D");
170		pvtxtree_->Branch("recy2_err_spl",&recy2_err_spl_,"recy2_err_spl[nrecPV2_spl]/D");
171		pvtxtree_->Branch("recz2_err_spl",&recz2_err_spl_,"recz2_err_spl[nrecPV2_spl]/D");
172	<
172	>
173	>	//pixeVertices
174	>	pvtxtree_->Branch("nrecPV_pxlpvtx",&nrecPV_pxlpvtx_,"nrecPV_pxlpvtx/I");
175	>	pvtxtree_->Branch("nTrkPV_pxlpvtx",&nTrkPV_pxlpvtx_,"nTrkPV_pxlpvtx[nrecPV_pxlpvtx]/I");
176	>	pvtxtree_->Branch("sumptsq_pxlpvtx",&sumptsq_pxlpvtx_,"sumptsq_pxlpvtx[nrecPV_pxlpvtx]/D");
177	>	pvtxtree_->Branch("isValid_pxlpvtx",&isValid_pxlpvtx_,"isValid_pxlpvtx/I");
178	>	pvtxtree_->Branch("isFake_pxlpvtx",&isFake_pxlpvtx_,"isFake_pxlpvtx/I");
179	>	pvtxtree_->Branch("recx_pxlpvtx",&recx_pxlpvtx_,"recx_pxlpvtx[nrecPV_pxlpvtx]/D");
180	>	pvtxtree_->Branch("recy_pxlpvtx",&recy_pxlpvtx_,"recy_pxlpvtx[nrecPV_pxlpvtx]/D");
181	>	pvtxtree_->Branch("recz_pxlpvtx",&recz_pxlpvtx_,"recz_pxlpvtx[nrecPV_pxlpvtx]/D");
182	>	pvtxtree_->Branch("recx_err_pxlpvtx",&recx_err_pxlpvtx_,"recx_err_pxlpvtx[nrecPV_pxlpvtx]/D");
183	>	pvtxtree_->Branch("recy_err_pxlpvtx",&recy_err_pxlpvtx_,"recy_err_pxlpvtx[nrecPV_pxlpvtx]/D");
184	>	pvtxtree_->Branch("recz_err_pxlpvtx",&recz_err_pxlpvtx_,"recz_err_pxlpvtx[nrecPV_pxlpvtx]/D");
185	>
186		//Fill the variables in the twovtx pair (recvtx1, recvtx2)
187		pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
188		pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
#	Line 203 | Line 235 | PVStudy::PVStudy(const edm::ParameterSet
235		pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
236		}
237		}
206	–
207	–	edm::Service<TFileService> fs;
208	–	TFileDirectory subDir = fs->mkdir( "Summary" );
209	–	TFileDirectory subDir1 = fs->mkdir( "Others" );
210	–	//   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
211	–
212	–	setRootStyle();
238
239	<	//Book histograms for track and pvtx parameters for each splitted pvtx
215	<	//i=0 : x (as in unsplit track collection)
216	<	//i=1 : x1_spl_
217	<	//i=2 : x2_spl_
239	>	setRootStyle();
240
241	<
242	<	for(int i=0;i<3;i++)
243	<	{
244	<	string suffix;
245	<	if(i == 0) suffix = "";
246	<	else {
247	<	stringstream ss;
248	<	ss << i;
249	<	suffix =ss.str()+"_spl";
250	<	}
251	<	h["nTrk"+suffix]   = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
252	<	h["trkPt"+suffix]  = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
253	<	h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#eta of rec tracks "+suffix),100,-3,3);
254	<	h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
255	<	h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-0.5,0.5);
256	<	h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
257	<
258	<	h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
259	<	h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
260	<	h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#eta of rec tracks in PV"+suffix),100,-3,3);
261	<	h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
262	<	h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);
263	<	h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);
264	<	h["nrecPV"+suffix]   = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);
265	<	suffix.clear();
241	>	//========================================
242	>	// Booking histograms
243	>	//========================================
244	>
245	>	// Create a root file for histograms
246	>	theFile = new TFile(histoFileName_.c_str(), "RECREATE");
247	>	// make diretories
248	>	theFile->mkdir("Summary");
249	>	theFile->cd();
250	>	theFile->mkdir("Others");
251	>	theFile->cd();
252	>	if (analyze_) {
253	>	theFile->mkdir("Results");
254	>	theFile->cd();
255	>	}
256	>
257	>	// Book Histograms:
258	>	h_pvtrk   = new PVHistograms();
259	>	h_pixvtx  = new PVHistograms();
260	>	h_misc    = new PVHistograms();
261	>	h_summary = new PVHistograms();
262	>	h_others  = new PVHistograms();
263	>
264	>	for(int i=0;i<3;i++) {
265	>	if(i == 0) h_pvtrk->Init("pvTrk");
266	>	else {
267	>	stringstream ss;
268	>	ss << i;
269	>	h_pvtrk->Init("pvTrk", ss.str(),"spl");
270		}
271	<	h["nrecPVDiff"]     = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
272	<	h["nTrkPVDiff"]     = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
273	<	h["nTrkPVRelDiff"]  = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
274	<
275	<	// Histograms on comparing the multi-vertices
276	<	// Difference in reconstructed vtx position
277	<	h["min_xsep"]   = subDir.make<TH1D>("min_xsep", "min x diff of primary and secondary pvtx",300,0,0.1);
278	<	h["min_xsep_sign"]   = subDir.make<TH1D>("min_xsep_sign", "min x diff in signf of primary and secondary pvtx",300,0,5);
279	<	h["min_ysep"]   = subDir.make<TH1D>("min_ysep", "min y diff of primary and secondary pvtx",300,0,0.1);
280	<	h["min_ysep_sign"]   = subDir.make<TH1D>("min_ysep_sign", "min y diff in signf of primary and secondary pvtx",300,0,5);
255	<	h["min_zsep"]   = subDir.make<TH1D>("min_zsep", "min z diff of primary and secondary pvtx",300,0,5);
256	<	h["min_zsep_sign"]   = subDir.make<TH1D>("min_zsep_sign", "min z diff in signf of primary and secondary pvtx",300,0,200);
257	<	// Difference in reconstructed vtx position
258	<	h["min_ntrksep"]   = subDir.make<TH1D>("min_ntrksep", "min nTrk diff of primary and secondary pvtx",201,-50.5,150.5);
259	<	h["min_sumpt2sep"]   = subDir.make<TH1D>("min_sumpt2sep", "min sumpt2 diff of primary and secondary pvtx",300,0,10000);
271	>	}
272	>	h_pixvtx->Init("pixVtx");
273	>	h_misc->Init("misc");
274	>
275	>	// Book MC only plots
276	>	if (!realData_) {
277	>	h_gen = new PVHistograms();
278	>	h_gen->Init("generator");
279	>	}
280	>	if (analyze_) h_ana = new PVHistograms();
281
282		//Book histograms sensitive to data/mc
283		for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
284		string suffix;
285	<	if(verbose_)  cout<<"datamode = "<< *it<<endl;
285	>	edm::LogInfo("Debug")<<"datamode = "<< *it<<endl;
286		switch(*it) {
287		case 0: suffix = "";
288		break;
#	Line 272 | Line 293 | PVStudy::PVStudy(const edm::ParameterSet
293		case 3: suffix = "_mct";
294		break;
295		}
296	<	h["deltax"+suffix] = subDir.make<TH1D>(TString("deltax"+suffix), TString("x-residual pvtx"+suffix),800,-0.04,0.04);
276	<	h["deltay"+suffix] = subDir.make<TH1D>(TString("deltay"+suffix), TString("y-residual pvtx"+suffix),800,-0.04,0.04);
277	<	h["deltaz"+suffix] = subDir.make<TH1D>(TString("deltaz"+suffix), TString("z-residual pvtx"+suffix),800,-0.04,0.04);
278	<	h["pullx"+suffix]  = subDir.make<TH1D>(TString("pullx"+suffix), TString("x-pull pvtx"+suffix),800,-10,10);
279	<	h["pully"+suffix]  = subDir.make<TH1D>(TString("pully"+suffix), TString("y-pull pvtx"+suffix),800,-10,10);
280	<	h["pullz"+suffix]  = subDir.make<TH1D>(TString("pullz"+suffix), TString("z-pull pvtx"+suffix),800,-10,10);
281	<	h["errPVx"+suffix] = subDir.make<TH1D>(TString("errPVx"+suffix), TString("X"+suffix+" vertex error"),100,0.,0.02);
282	<	h["errPVy"+suffix] = subDir.make<TH1D>(TString("errPVy"+suffix), TString("Y"+suffix+" vertex error"),100,0.,0.02);
283	<	h["errPVz"+suffix] = subDir.make<TH1D>(TString("errPVz"+suffix), TString("Z"+suffix+" vertex error"),100,0.,0.02);
296	>	h_summary->Init("summary",suffix);
297
298		// Book residual, error and pull histograms for each nTrk bin
299		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
300		stringstream ss;
301		ss << ntrk;
302	<	string suffix2 = suffix+"_"+ss.str();
303	<	h["deltax"+suffix2]  = subDir1.make<TH1D>(TString("deltax"+suffix2), TString("x-residual of pvtx"+suffix),100,-0.02,0.02);
304	<	h["deltax"+suffix2]->GetXaxis()->SetTitle("cm");
292	<	h["deltay"+suffix2]  = subDir1.make<TH1D>(TString("deltay"+suffix2), TString("y-residual of pvtx"+suffix),100,-0.02,0.02);
293	<	h["deltay"+suffix2]->GetXaxis()->SetTitle("cm");
294	<	h["deltaz"+suffix2]  = subDir1.make<TH1D>(TString("deltaz"+suffix2), TString("z-residual of pvtx"+suffix),100,-0.02,0.02);
295	<	h["deltaz"+suffix2]->GetXaxis()->SetTitle("cm");
296	<	h["pullx"+suffix2] = subDir1.make<TH1D>(TString("pullx"+suffix2), TString("x-pull of pvtx"+suffix),100,-10.,10.);
297	<	h["pully"+suffix2] = subDir1.make<TH1D>(TString("pully"+suffix2), TString("y-pull of pvtx"+suffix),100,-10.,10.);
298	<	h["pullz"+suffix2] = subDir1.make<TH1D>(TString("pullz"+suffix2), TString("z-pull of pvtx"+suffix),100,-10.,10.);
299	<	h["errPVx"+suffix2] = subDir1.make<TH1D>(TString("errPVx"+suffix2), TString("X"+suffix+" vertex error"),100,0.,0.02);
300	<	h["errPVy"+suffix2] = subDir1.make<TH1D>(TString("errPVy"+suffix2), TString("Y"+suffix+" vertex error"),100,0.,0.02);
301	<	h["errPVz"+suffix2] = subDir1.make<TH1D>(TString("errPVz"+suffix2), TString("Z"+suffix+" vertex error"),100,0.,0.02);
302	<	suffix2.clear();
303	<	} // End of   for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
304	<
302	>	h_others->Init("others",suffix,ss.str());
303	>	}
304	>
305		// Book residual and pull histograms only when analyzeOntheFly is enabled
306	<	if(analyze_) {
307	<	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);
308	<	h2["resx"+suffix+"_nTrk"]->SetMarkerStyle(21);
309	<	h2["resx"+suffix+"_nTrk"]->SetMarkerColor(4);
310	<	h2["resx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
311	<	h2["resx"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
312	<	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);
313	<	h2["resy"+suffix+"_nTrk"]->SetMarkerStyle(21);
314	<	h2["resy"+suffix+"_nTrk"]->SetMarkerColor(4);
315	<	h2["resy"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
316	<	h2["resy"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
317	<	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);
318	<	h2["resz"+suffix+"_nTrk"]->SetMarkerStyle(21);
319	<	h2["resz"+suffix+"_nTrk"]->SetMarkerColor(4);
320	<	h2["resz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
321	<	h2["resz"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
322	<	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.);
323	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerStyle(21);
324	<	h2["pullx"+suffix+"_nTrk"]->SetMarkerColor(4);
325	<	h2["pullx"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
326	<	h2["pullx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
327	<	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.);
328	<	h2["pully"+suffix+"_nTrk"]->SetMarkerStyle(21);
329	<	h2["pully"+suffix+"_nTrk"]->SetMarkerColor(4);
330	<	h2["pully"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
331	<	h2["pully"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
332	<	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.);
333	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerStyle(21);
334	<	h2["pullz"+suffix+"_nTrk"]->SetMarkerColor(4);
335	<	h2["pullz"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
336	<	h2["pullz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
337	<	} // End of   if(analyze_) {
306	>	if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
307		suffix.clear();
308		} // End of Book histograms sensitive to data/mc
309
310	<	// Book histograms about pixelVertices
342	<	h["trkdz_pxlpvtxdz"]   = subDir.make<TH1D>("trkdz_pxlpvtxdz", "(Track dz - pixelpvtx dz) in cm",300,-0.5,0.5);
343	<	h["trkdz_pxlpvtxdz_pxlpvtxdzerr"]   = subDir.make<TH1D>("trkdz_pxlpvtxdz_pxlpvtxdzerr", "|Track dz - pixelpvtx dz| / pxlpvtxdzErr",300,0,100);
344	<	h["trkdz_pxlpvtxdz_trkdzerr"]   = subDir.make<TH1D>("trkdz_pxlpvtxdz_trkdzerr", "|Track dz - pixelpvtx dz| / trkdzErr",300,0,50);
345	<	h["trkdzErr_pxlpvtx"]   = subDir.make<TH1D>("trkdzErr_pxlpvtxdz", "Track dzErr of leading pixelpvtx ",300,0,0.5);
346	<	h["trkdzErr_pvtx"]   = subDir.make<TH1D>("trkdzErr_pvtx", "Track dzErr of the leading pvtx ",300,0,0.5);
347	<	h["dzErr_pxlpvtx"]   = subDir.make<TH1D>("dzErr_pxlpvtx", "zError of the leading pvtx ",300,0,0.5);
348	<
349	<	// Book MC only plots
350	<	if (!realData_) {
351	<	h["genPart_T"]      = subDir.make<TH1D>("genPart_T","t component of gen particles",300,-0.5,0.5);
352	<	h["genPart_T"]->GetXaxis()->SetTitle("t (nanosecond)");
353	<	h["genPart_cT"]     = subDir.make<TH1D>("genPart_cT","ct component of gen particles",300,-150.,150.);
354	<	h["genPart_cT"]->GetXaxis()->SetTitle("ct (mm)");
355	<	h["nsimPV"]         = subDir.make<TH1D>("nsimPV","Num of sim PV",51,-0.5,50.5);
356	<	}
310	>
311		}
312
313		PVStudy::~PVStudy()
#	Line 361 | Line 315 | PVStudy::~PVStudy()
315
316		// do anything here that needs to be done at desctruction time
317		// (e.g. close files, deallocate resources etc.)
318	+	theFile->cd();
319	+	theFile->cd("Summary");
320	+	h_pvtrk->Save();
321	+	h_pixvtx->Save();
322	+	h_misc->Save();
323	+	h_summary->Save();
324	+	if (!realData_)
325	+	h_gen->Save();
326	+	if (analyze_) {
327	+	theFile->cd();
328	+	theFile->cd("Results");
329	+	h_ana->Save();
330	+	}
331	+	theFile->cd();
332	+	theFile->cd("Others");
333	+	h_others->Save();
334
335	+	theFile->Close();
336		}
337
338		void PVStudy::setRootStyle() {
#	Line 398 | Line 369 | void PVStudy::setRootStyle() {
369		//
370		std::vector<PVStudy::simPrimaryVertex> PVStudy::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
371		{
372	<	std::vector<PVStudy::simPrimaryVertex> simpv;
372	>	std::vector<PVStudy::simPrimaryVertex> simpv;
373		const HepMC::GenEvent* evt=evtMC->GetEvent();
374		if (evt) {
375	<	if(verbose_) {
376	<	cout << "process id " << evt->signal_process_id()<<endl;
377	<	cout << "signal process vertex " << ( evt->signal_process_vertex() ?
378	<	evt->signal_process_vertex()->barcode() : 0 )   <<endl;
408	<	cout << "number of vertices " << evt->vertices_size() << endl;
409	<	}
375	>	edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
376	>	edm::LogInfo("SimPVs") << "[getSimPVs] signal process vertex " << ( evt->signal_process_vertex() ?
377	>	evt->signal_process_vertex()->barcode() : 0 ) <<endl;
378	>	edm::LogInfo("SimPVs") << "[getSimPVs] number of vertices " << evt->vertices_size() << endl;
379
380		int idx=0; int npv=0;
381		for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
#	Line 417 | Line 386 | std::vector<PVStudy::simPrimaryVertex> P
386		// t component of PV:
387		for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
388		mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
389	<	//        std::cout << "Status = " << (*mother)->status() << std::endl;
389	>	//        edm::LogInfo("SimPVs") << "Status = " << (*mother)->status() << endl;
390		HepMC::GenVertex * mv=(*mother)->production_vertex();
391		if( ((*mother)->status() == 3) && (!mv)) {
392	<	//        std::cout << "npv= " << npv << std::endl;
392	>	//        edm::LogInfo("SimPVs") << "npv= " << npv << endl;
393		if (npv == 0) {
394	<	h["genPart_cT"]->Fill(pos.t()); // mm
395	<	h["genPart_T"]->Fill(pos.t()/299.792458); // ns
394	>	h_gen->Fill1d("genPart_cT", pos.t()); // mm
395	>	h_gen->Fill1d("genPart_T", pos.t()/299.792458); // ns
396		}
397		npv++;
398		}
#	Line 431 | Line 400 | std::vector<PVStudy::simPrimaryVertex> P
400		//       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
401
402		bool hasMotherVertex=false;
403	<	if(verbose_) cout << "mothers of vertex[" << ++idx << "]: " << endl;
403	>	if (verbose_) cout << "[getSimPVs] mothers of vertex[" << ++idx << "]: " << endl;
404		for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
405		mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
406		HepMC::GenVertex * mv=(*mother)->production_vertex();
407	<	//        std::cout << "Status = " << (*mother)->status() << std::endl;
407	>	// if (verbose_) cout << "Status = " << (*mother)->status() << endl;
408		if (mv) {
409		hasMotherVertex=true;
410		if(!verbose_) break; //if verbose_, print all particles of gen vertices
#	Line 462 | Line 431 | std::vector<PVStudy::simPrimaryVertex> P
431
432		if(!vp){
433		// this is a new vertex
434	<	if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
434	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
435		simpv.push_back(sv);
436		vp=&simpv.back();
437		}else{
438	<	if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
438	>	edm::LogInfo("SimPVs") << "[getSimPVs] this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
439		}
440		vp->genVertex.push_back((*vitr)->barcode());
441		// collect final state descendants
#	Line 490 | Line 459 | std::vector<PVStudy::simPrimaryVertex> P
459		}
460		}
461		}
462	<	}
463	<	}
462	>	}//loop MC vertices daughters
463	>	}//loop MC vertices
464		}
465		return simpv;
466		}
#	Line 524 | Line 493 | std::vector<PVStudy::simPrimaryVertex> P
493		ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
494		}
495		}
496	<	primary =  ( gv->position().t()==0);
496	>	//primary =  ( gv->position().t()==0 );
497	>	primary = true;
498	>	h_gen->Fill1d("genPart_cT", gv->position().t()); // mm
499	>	h_gen->Fill1d("genPart_T", gv->position().t()/299.792458); // ns
500	>
501		//resonance= ( gp->mother() && isResonance(gp->mother()));  // in CLHEP/HepMC days
502		// no more mother pointer in the improved HepMC GenParticle
503		resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
#	Line 573 | Line 546 | PVStudy::analyze(const edm::Event& iEven
546		using namespace edm;
547		using namespace reco;
548
549	<	if (verbose_)
577	<	cout<<"PVStudy::analyze():"<<endl;
549	>	edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
550		//=======================================================
551		// Initialize Root-tuple variables if needed
552		//=======================================================
553		//if(saventuple_)
554	<	SetVarToZero();
554	>	SetVarToZero();
555
556		//=======================================================
557		// Track accessors
#	Line 592 | Line 564 | PVStudy::analyze(const edm::Event& iEven
564		if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
565		trackColl = trackCollectionHandle.product();
566		} else {
567	<	cout << "trackCollection cannot be found -> using empty collection of same type." <<endl;
567	>	edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
568		}
597	–
569		//splitTrackCollection1
570		static const reco::TrackCollection s_empty_splitTrackColl1;
571		const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
#	Line 603 | Line 574 | PVStudy::analyze(const edm::Event& iEven
574		if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
575		splitTrackColl1 = splitTrackCollection1Handle.product();
576		} else {
577	<	cout << "splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
577	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
578		}
608	–
579		//splitTrackCollection2
580		static const reco::TrackCollection s_empty_splitTrackColl2;
581		const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
#	Line 614 | Line 584 | PVStudy::analyze(const edm::Event& iEven
584		if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
585		splitTrackColl2 = splitTrackCollection2Handle.product();
586		} else {
587	<	cout << "splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
587	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
588		}
589	<
620	<
621	<	//=======================================================
622	<	// Fill trackparameters of the input tracks to pvtx fitter
623	<	//=======================================================
624	<	if(verbose_)
625	<	cout<<"Start filling track parameters of the input tracks to pvtx fitter."<<endl;
626	<	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype)
627	<	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
628	<	fillTrackHisto(trackColl, 0);
629	<	fillTrackHisto(splitTrackColl1, 1);
630	<	fillTrackHisto(splitTrackColl2, 2);
631	<	if(verbose_)
632	<	cout<<"End filling track parameters of the input tracks to pvtx fitter."<<endl;
633	<
589	>
590		//=======================================================
591		// PVTX accessors
592		//=======================================================
#	Line 642 | Line 598 | PVStudy::analyze(const edm::Event& iEven
598		if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
599		vertexColl = vertexCollectionHandle.product();
600		} else {
601	<	cout << "vertexCollection cannot be found -> using empty collection of same type." <<endl;
601	>	edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
602		}
647	–
603		//splitVertexCollection1
604		static const reco::VertexCollection s_empty_splitVertexColl1;
605		const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
#	Line 653 | Line 608 | PVStudy::analyze(const edm::Event& iEven
608		if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
609		splitVertexColl1 = splitVertexCollection1Handle.product();
610		} else {
611	<	cout << "splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
612	<	}
658	<
611	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
612	>	}
613		//splitVertexCollection2
614		static const reco::VertexCollection s_empty_splitVertexColl2;
615		const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
#	Line 664 | Line 618 | PVStudy::analyze(const edm::Event& iEven
618		if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
619		splitVertexColl2 = splitVertexCollection2Handle.product();
620		} else {
621	<	cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
621	>	edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
622		}
623	+
624
625	<	if(verbose_) cout<<"End accessing the track and vertex collections"<<endl;
625	>	//=======================================================
626	>	// GET pixelVertices
627	>	//=======================================================
628	>	static const reco::VertexCollection s_empty_pixelVertexColl;
629	>	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
630	>	edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
631	>	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
632	>	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
633	>	pixelVertexColl = pixelVertexCollectionHandle.product();
634	>	} else {
635	>	edm::LogInfo("Debug") << "[PVStudy] pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
636	>	}
637	>
638	>	//=======================================================
639	>	// BeamSpot accessors
640	>	//=======================================================
641	>
642	>	edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
643	>	iEvent.getByLabel(bsSrc,recoBeamSpotHandle);
644	>	reco::BeamSpot bs = *recoBeamSpotHandle;
645	>	const Point beamSpot = recoBeamSpotHandle.isValid() ? Point(recoBeamSpotHandle->x0(), recoBeamSpotHandle->y0(), recoBeamSpotHandle->z0()) : Point(0, 0, 0);
646	>
647	>	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex and pixelvertices collections"<<endl;
648
649		//=======================================================
650		// MC simvtx accessor
#	Line 686 | Line 663 | PVStudy::analyze(const edm::Event& iEven
663		try{
664		iSetup.getData(pdt);
665		}catch(const Exception&){
666	<	cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
666	>	edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
667		}
668	<
668	>
669		//=======================================================
670	<	// GET pixelVertices
670	>	// Apply event cleaning for firstcoll data and MC
671		//=======================================================
672	<	static const reco::VertexCollection s_empty_pixelVertexColl;
673	<	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
674	<	edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
675	<	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
676	<	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
677	<	pixelVertexColl = pixelVertexCollectionHandle.product();
678	<	} else {
679	<	cout << "pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
680	<	}
672	>	if(applyEvtClean_) {
673	>	// =====Select on the trigger bits
674	>	edm::ESHandle<L1GtTriggerMenu> menuRcd;
675	>	iSetup.get<L1GtTriggerMenuRcd>().get(menuRcd) ;
676	>	const L1GtTriggerMenu* menu = menuRcd.product();
677	>	edm::Handle< L1GlobalTriggerReadoutRecord > gtRecord;
678	>	iEvent.getByLabel( edm::InputTag("gtDigis"), gtRecord);
679	>
680	>	bool isTechBit0 = false;
681	>	bool isTechBit40 = false;
682	>	bool isBeamHalo = false;
683	>
684	>	TechnicalTriggerWord tw = gtRecord->technicalTriggerWord();
685	>	if ( ! tw.empty() ) {
686	>	// loop over dec. bit to get total rate (no overlap)
687	>	for ( int i = 0; i < 64; ++i ) {
688	>	if ( tw[i] ) {
689	>	//cout<<"technical number "<<i<<"  "<<endl;
690	>	if (i == 0) isTechBit0  = true;
691	>	if (i < 40 && i > 35) isBeamHalo = true;   // The beamHalo bits are 36-39
692	>	if (i == 40 || i == 41) isTechBit40 = true;
693	>	}
694	>	}
695	>	} // =====End select on the trigger bits
696	>	// =====Loop over the trackColl to tet the fraction of HighPurity tracks
697	>	int nTracks = 0;
698	>	int nHighPurityTracks=0;
699	>	double fHighPurity=0;
700	>
701	>	for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
702	>	TrackRef tkref(trackCollectionHandle,i);
703	>	if( (tkref->qualityMask() & 4 ) == 4) ++nHighPurityTracks;
704	>	}
705	>	if(nTracks>0)
706	>	fHighPurity =  double(nHighPurityTracks)/double(nTracks);
707	>	if(verbose_)
708	>	cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;
709	>
710	>	// Apply the event selection for MC events
711	>	if(realData_) {
712	>	if(!isTechBit40 || isBeamHalo || (fHighPurity<0.2&&nTracks>10) || vertexColl->begin()->isFake()) {
713	>	glob_runno_ = iEvent.id().run();
714	>	glob_evtno_ = iEvent.id().event();
715	>	glob_ls_   = iEvent.luminosityBlock();
716	>	glob_bx_  = iEvent.bunchCrossing();
717	>	return;
718	>	}
719	>	else {
720	>	// TechBit 0, beamHalo are not simulted in MC
721	>	if ( !isTechBit0 || !isTechBit40 || isBeamHalo || (fHighPurity<0.2&&nTracks>10) || vertexColl->begin()->isFake()) return;
722	>	}
723	>	}
724	>	} // End of Apply event cleaning cuts for firstcoll data
725	>
726	>
727	>
728	>	//=======================================================
729	>	// Fill trackparameters of the input tracks to pvtx fitter
730	>	//=======================================================
731	>	edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
732	>	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs)
733	>	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
734	>	fillTrackHisto(trackColl, 0, beamSpot);
735	>	fillTrackHisto(splitTrackColl1, 1, beamSpot);
736	>	fillTrackHisto(splitTrackColl2, 2, beamSpot);
737	>	edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
738	>
739
740		//=======================================================
741		// Fill pixelVertices related histograms
742		//=======================================================
743	+	nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
744		if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
745	<	h["dzErr_pxlpvtx"]->Fill( pixelVertexColl->begin()->zError());
745	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
746	>	fillTrackHistoInPV(pixelVertexColl, 4, false, true, beamSpot);
747	>	h_pixvtx->Fill1d("dzErr_pxlpvtx", pixelVertexColl->begin()->zError());
748		// Get the dZ error of the tracks in the leading pixelVertexColl
749		for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
750		t!= (pixelVertexColl->begin())->tracks_end(); t++) {
751
752		if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
753	<	//        h["trkPtPV"]->Fill(0.);
753	>	//        h_pvtrk->Fill1d("trkPtPV", 0.);
754		}
755		else
756	<	h["trkdzErr_pxlpvtx"]->Fill((**t).dzError());
756	>	h_pixvtx->Fill1d("trkdzErr_pxlpvtx", (**t).dzError());
757		}
758		// Fill the track->dz() in the PV difference from first pixelpvtx
759		for(reco::Vertex::trackRef_iterator t = (vertexColl->begin())->tracks_begin();
760		t!= (vertexColl->begin())->tracks_end(); t++) {
761		// illegal charge
762		if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
763	<	//        h["trkPtPV"]->Fill(0.);
763	>	//        h_pvtrk->Fill1d("trkPtPV", 0.);
764		}
765		else {
766		if(pixelVertexColl->size()>0) {
767	<	h["trkdz_pxlpvtxdz"]->Fill((**t).dz() -  pixelVertexColl->begin()->z());
768	<	h["trkdz_pxlpvtxdz_pxlpvtxdzerr"]->Fill(fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
769	<	h["trkdz_pxlpvtxdz_trkdzerr"]->Fill(fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
770	<	h["trkdzErr_pvtx"]->Fill((**t).dzError());
767	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz", (**t).dz() -  pixelVertexColl->begin()->z());
768	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_pxlpvtxdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
769	>	h_pixvtx->Fill1d("trkdz_pxlpvtxdz_trkdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
770	>	h_pixvtx->Fill1d("trkdzErr_pvtx", (**t).dzError());
771		}
772		}
773		}
774		}
737	–
775
776		//=======================================================
777	<	// Fill number of reconstructed vertices ootb
777	>	// Fill number of reconstructed vertices
778		//=======================================================
779
780	<	if(verbose_)  {
781	<	cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
780	>	edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
781	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
782	>	edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
783	>	edm::LogInfo("Debug")<<"[PVStudy] Start filling pvtx parameters."<<endl;
784	>	if (verbose_) {
785		printRecVtxs(vertexCollectionHandle);
746	–	cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
786		printRecVtxs(splitVertexCollection1Handle);
748	–	cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
787		printRecVtxs(splitVertexCollection2Handle);
750	–	cout<<"Start filling pvtx parameters."<<endl;
788		}
752	–
789		nrecPV_ = int (vertexColl->size());
790		nrecPV1_spl_ = int (splitVertexColl1->size());
791		nrecPV2_spl_ = int (splitVertexColl2->size());
792
793	<	h["nrecPV"]->Fill(nrecPV_);
794	<	h["nrecPV1_spl"]->Fill(nrecPV1_spl_);
795	<	h["nrecPV2_spl"]->Fill(nrecPV2_spl_);
796	<	h["nrecPVDiff"]->Fill(double(nrecPV1_spl_)-double(nrecPV2_spl_));
797	<
798	<	//=================================================================
793	>	h_pvtrk->Fill1d("nrecPV", nrecPV_);
794	>	h_pvtrk->Fill1d("nrecPV1_spl", nrecPV1_spl_);
795	>	h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
796	>	h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
797	>
798	>
799	>	//=================================================================
800		// Fill track parameter ntuple/hist for tracks used in recoVertices
801		//=================================================================
802
803	<	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
803	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
804		if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
805	<	fillTrackHistoInPV(vertexColl, 0, true, true);
805	>	fillTrackHistoInPV(vertexColl, 0, true, true, beamSpot);
806
807		if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
808	<	fillTrackHistoInPV(splitVertexColl1, 1, false, true);
808	>	fillTrackHistoInPV(splitVertexColl1, 1, false, true, beamSpot);
809
810	<	if(splitVertexColl1->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
811	<	fillTrackHistoInPV(splitVertexColl2, 2, false, true);
810	>	if(splitVertexColl2->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
811	>	fillTrackHistoInPV(splitVertexColl2, 2, false, true, beamSpot);
812
813	<	//==========================================================
813	>	//=======================================================
814	>	// Compare offlinePrimaryVertices with pixelVertices
815	>	//=======================================================
816	>	if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
817	>	&& (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {
818	>	h_pixvtx->Fill1d("nrecPV_minus_nrecPxlPV", double (nrecPV_ - nrecPV_pxlpvtx_));
819	>	// difference in reconstructed position of the leading pvtx
820	>	//edm::LogInfo("Debug")<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;
821	>	//edm::LogInfo("Debug")<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
822	>	h_pixvtx->Fill1d("recxPV_minus_recxPxlPV", recx_[0] - recx_pxlpvtx_[0]);
823	>	h_pixvtx->Fill1d("recyPV_minus_recyPxlPV", recy_[0] - recy_pxlpvtx_[0]);
824	>	h_pixvtx->Fill1d("reczPV_minus_reczPxlPV", recz_[0] - recz_pxlpvtx_[0]);
825	>	}
826	>
827	>
828	>	//==========================================================
829		// Fill secondary/primary min separations for multi-vertices
830		//==========================================================
831
#	Line 788 | Line 840 | PVStudy::analyze(const edm::Event& iEven
840		double min_ntrksep = 9999.0;
841		double min_sumpt2sep = 9999999.0;
842
843	<	if(verbose_) cout<<"leading pvtx z = "<< vertexColl->begin()->z() <<endl;
843	>	edm::LogInfo("Debug")<<"[PVStudy] leading pvtx z = "<< vertexColl->begin()->z() <<endl;
844
845		// Looping through the secondary vertices to find the mininum separation to the primary
846		for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
#	Line 823 | Line 875 | PVStudy::analyze(const edm::Event& iEven
875		min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
876		}
877		}
878	<	h["min_xsep"]->Fill(min_xsep);
879	<	h["min_ysep"]->Fill(min_ysep);
880	<	h["min_zsep"]->Fill(min_zsep);
881	<	h["min_xsep_sign"]->Fill(min_xsep_sign);
882	<	h["min_ysep_sign"]->Fill(min_ysep_sign);
883	<	h["min_zsep_sign"]->Fill(min_zsep_sign);
884	<	h["min_ntrksep"]->Fill(min_ntrksep);
885	<	h["min_sumpt2sep"]->Fill(min_sumpt2sep);
878	>	h_misc->Fill1d("min_xsep", min_xsep);
879	>	h_misc->Fill1d("min_ysep", min_ysep);
880	>	h_misc->Fill1d("min_zsep", min_zsep);
881	>	h_misc->Fill1d("min_xsep_sign", min_xsep_sign);
882	>	h_misc->Fill1d("min_ysep_sign", min_ysep_sign);
883	>	h_misc->Fill1d("min_zsep_sign", min_zsep_sign);
884	>	h_misc->Fill1d("min_ntrksep", min_ntrksep);
885	>	h_misc->Fill1d("min_sumpt2sep", min_sumpt2sep);
886
887		min_zsep_ = min_zsep;
888		min_ntrksep_ = min_ntrksep;
#	Line 839 | Line 891 | PVStudy::analyze(const edm::Event& iEven
891
892		} // End of  if(vertexColl->size()>1) {
893
894	<	if(verbose_)
895	<	cout<<"End filling pvtx parameters."<<endl;
894	>	edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
895	>
896
897	+
898		//=======================================================
899		//           PrimaryVertex Matching
900		// 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
#	Line 849 | Line 902 | PVStudy::analyze(const edm::Event& iEven
902		// Assume the first match is the primary vertex,
903		// since vertexColl are sorted in decreasing order of sum(pT)
904		//=======================================================
905	<
906	<	if(verbose_)   cout<<"PVStudy::analyze(): matching pvtxs "<<endl;
905	>
906	>	edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
907		reco::VertexCollection s_empty_matchedVertexColl1;
908		reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
909		matchedVertexColl1->reserve(splitVertexColl1->size());
#	Line 869 | Line 922 | PVStudy::analyze(const edm::Event& iEven
922		if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
923		reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
924		if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
925	<	if(verbose_) {
926	<	cout<<"Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
927	<	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
928	<	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
929	<	}
925	>	edm::LogInfo("Debug")<<"[PVStudy] Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
926	>	edm::LogInfo("Debug")<<"[PVStudy] recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
927	>	edm::LogInfo("Debug")<<"[PVStudy] recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
928	>
929	>	double twovtxsig = (recvtx1->z()-recvtx2->z())/max(recvtx1->zError(), recvtx2->zError());
930	>	h_misc->Fill1d("twovtxzsign", twovtxsig);
931		if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
932	<	if(verbose_) {
933	<	cout<<"The two splitted vertices match in Z. "<<endl;
880	<	}
932	>	edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
933	>
934		int nTrkPV1 = recvtx1->tracksSize();
935		int nTrkPV2 = recvtx2->tracksSize();
936		double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
937	<	h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
938	<	h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
937	>	h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
938	>	h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
939		if(fabs(ntrkreldiff)<ntrkdiffcut_) {
940	<	if(verbose_) {
941	<	cout<<"(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
889	<	}
940	>	edm::LogInfo("Debug")<<"[PVStudy] (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
941	>
942		matchedVertexColl1->push_back(*recvtx1);
943		matchedVertexColl2->push_back(*recvtx2);
944		stopmatching = true; // stop the matching when the first match is found!
945		break;
946		}
947		else
948	<	cout<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
948	>	edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
949		}
950		else {
951	<	cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
952	<	cout<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
953	<	cout<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
951	>	edm::LogInfo("Debug")<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
952	>	edm::LogInfo("Debug")<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
953	>	edm::LogInfo("Debug")<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
954		}
955		}
956		}
957	<	if(verbose_)     cout<<"PVStudy::analyze(): End matching pvtxs"<<endl;
957	>	edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
958
959		//=======================================================
960		//  Analyze matchedVertexColls
961		//=======================================================
962	<	if(verbose_) {
911	<	cout<<"Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
912	<	}
962	>	edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
963
964		// Compare the reconstructed vertex position and calculate resolution/pulls
965		if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
966		//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
967	<	fillTrackHistoInPV(matchedVertexColl1, 1, true, false);
968	<	fillTrackHistoInPV(matchedVertexColl2, 2, true, false);
967	>	fillTrackHistoInPV(matchedVertexColl1, 1, true, false, beamSpot);
968	>	fillTrackHistoInPV(matchedVertexColl2, 2, true, false, beamSpot);
969
970		reco::VertexCollection::const_iterator v1;
971		reco::VertexCollection::const_iterator v2;
972		nrecPV_twovtx_ = 0;
973		for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
974		v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
975	<	++v1, ++v2) {
975	>	++v1, ++v2) {
976	>	h_misc->Fill1d("ndofPVDiff", v1->ndof() - v2->ndof());
977	>	h_misc->Fill1d("ndofPVRelDiff", (v1->ndof()-v2->ndof())/(v1->ndof()+v2->ndof()));
978		fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
979		fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
980		fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
#	Line 931 | Line 983 | PVStudy::analyze(const edm::Event& iEven
983		ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
984		fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
985
986	<	h["deltax"]->Fill( fres_[0] );
987	<	h["deltay"]->Fill( fres_[1] );
988	<	h["deltaz"]->Fill( fres_[2] );
989	<	h["pullx"]->Fill( fres_[0]/ferror_[0]);
990	<	h["pully"]->Fill( fres_[1]/ferror_[1]);
991	<	h["pullz"]->Fill( fres_[2]/ferror_[2]);
992	<	h["errPVx"]->Fill( ferror_[0] );
993	<	h["errPVy"]->Fill( ferror_[1] );
994	<	h["errPVz"]->Fill( ferror_[2] );
986	>	h_summary->Fill1d("deltax", fres_[0] );
987	>	h_summary->Fill1d("deltay", fres_[1] );
988	>	h_summary->Fill1d("deltaz", fres_[2] );
989	>	h_summary->Fill1d("pullx", fres_[0]/ferror_[0]);
990	>	h_summary->Fill1d("pully", fres_[1]/ferror_[1]);
991	>	h_summary->Fill1d("pullz", fres_[2]/ferror_[2]);
992	>	h_summary->Fill1d("errPVx", ferror_[0] );
993	>	h_summary->Fill1d("errPVy", ferror_[1] );
994	>	h_summary->Fill1d("errPVz", ferror_[2] );
995		pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
996		error(ferror_[0], ferror_[1], ferror_[2]),
997		fntrk_));
#	Line 964 | Line 1016 | PVStudy::analyze(const edm::Event& iEven
1016		pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];
1017		nrecPV_twovtx_++;
1018		} // End of analyzing res/pull
1019	+
1020	+	//=======================================================
1021	+	// Fill simulated vertices
1022	+	//=======================================================
1023	+	if (!realData_) {
1024	+	bool MC=false;
1025	+	Handle<HepMCProduct> evtMC;
1026	+	iEvent.getByLabel("generator",evtMC);
1027	+	if (!evtMC.isValid()) {
1028	+	MC=false;
1029	+	edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
1030	+	} else {
1031	+	edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
1032	+	MC=true;
1033	+	}
1034	+
1035	+	if(MC){
1036	+	// make a list of primary vertices:
1037	+	std::vector<simPrimaryVertex> simpv;
1038	+	simpv=getSimPVs(evtMC,"");
1039	+	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
1040	+	h_gen->Fill1d("nsimPV", simpv.size());
1041	+	nsimPV_ = simpv.size();
1042	+
1043	+	int isimpv = 0;
1044	+	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
1045	+	vsim!=simpv.end(); vsim++, isimpv++){
1046	+	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
1047	+	simx_[isimpv] = vsim->x;
1048	+	simy_[isimpv] = vsim->y;
1049	+	simz_[isimpv] = vsim->y;
1050	+	simptsq_[isimpv] = vsim->ptsq;
1051	+
1052	+	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
1053	+	fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
1054	+	fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
1055	+	fillMCHisto(vsim, isimpv, vertexColl, 3);
1056	+	}
1057	+	} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
1058	+	} // End of   if (!realData_) {
1059	+
1060		} // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
1061		else
1062	<	cout<<"WARNING: Cannot find matching pvtxs"<<endl;
1062	>	edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
1063
1064	<	if(verbose_)     cout<<"End analyzing the matchedVertexColl"<<endl;
1064	>	edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
1065
1066	<	//=======================================================
974	<	// Fill simulated vertices
975	<	//=======================================================
976	<	if (!realData_) {
977	<	bool MC=false;
978	<	Handle<HepMCProduct> evtMC;
979	<	iEvent.getByLabel("generator",evtMC);
980	<	if (!evtMC.isValid()) {
981	<	MC=false;
982	<	if(verbose_){
983	<	cout << "no HepMCProduct found"<< endl;
984	<	}
985	<	} else {
986	<	if(verbose_){
987	<	cout << "generator HepMCProduct found"<< endl;
988	<	}
989	<	MC=true;
990	<	}
991	<
992	<	if(MC){
993	<	// make a list of primary vertices:
994	<	std::vector<simPrimaryVertex> simpv;
995	<	simpv=getSimPVs(evtMC,"");
996	<	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
997	<	h["nsimPV"]->Fill(simpv.size());
998	<	nsimPV_ = simpv.size();
999	<
1000	<	int isimpv = 0;
1001	<	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
1002	<	vsim!=simpv.end(); vsim++, isimpv++){
1003	<	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
1004	<	simx_[isimpv] = vsim->x;
1005	<	simy_[isimpv] = vsim->y;
1006	<	simz_[isimpv] = vsim->y;
1007	<	simptsq_[isimpv] = vsim->ptsq;
1008	<
1009	<	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
1010	<	fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
1011	<	fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
1012	<	fillMCHisto(vsim, isimpv, vertexColl, 3);
1013	<	}
1014	<	} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
1015	<	} // End of   if (!realData_) {
1066	>
1067
1068		// Finally fill the ftree_
1069		if(saventuple_) {
1070		ftree_->Fill();
1071		pvtxtree_->Fill();
1072		}
1073	+
1074		}
1075
1076	+
1077	+
1078		// ------------ method called once each job just before starting event loop  ------------
1079		void
1080		PVStudy::beginJob()
#	Line 1030 | Line 1084 | PVStudy::beginJob()
1084		// ------------ method called once each job just after ending the event loop  ------------
1085		void
1086		PVStudy::endJob() {
1087	<	if (verbose_) cout << "Analyzing PV info" << endl;
1087	>	edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl;
1088		// Looping through the datamodes available
1089		for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
1090		string suffix;
1091	<	if(verbose_)  cout<<"datamode = "<< *it<<endl;
1091	>	edm::LogInfo("Analysis")<<"[endJob] datamode = "<< *it<<endl;
1092		switch(*it) {
1093		case 0: suffix = "";
1094		break;
#	Line 1045 | Line 1099 | PVStudy::endJob() {
1099		case 3: suffix = "_mct";
1100		break;
1101		}
1102	+	suffix += "_nTrk";
1103		if(analyze_) {
1104		for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
1105	+	edm::LogInfo("Analysis") << "[endJob] ntrk = " << ntrk << endl;
1106		PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
1107	<	if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
1108	<	if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
1109	<	if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
1110	<	if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
1111	<	if ( ipv.pull_.y() > 0 ) h2["pully"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
1112	<	if ( ipv.pull_.z() > 0 ) h2["pullz"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.z(), 1.);
1107	>	if ( ipv.res_.x() > 0 ) h_ana->Fill2d("resx"+suffix, ntrk,ipv.res_.x());
1108	>	if ( ipv.res_.y() > 0 ) h_ana->Fill2d("resy"+suffix, ntrk,ipv.res_.y());
1109	>	if ( ipv.res_.z() > 0 ) h_ana->Fill2d("resz"+suffix, ntrk,ipv.res_.z());
1110	>	if ( ipv.pull_.x() > 0 ) h_ana->Fill2d("pullx"+suffix, ntrk,ipv.pull_.x());
1111	>	if ( ipv.pull_.y() > 0 ) h_ana->Fill2d("pully"+suffix, ntrk,ipv.pull_.y());
1112	>	if ( ipv.pull_.z() > 0 ) h_ana->Fill2d("pullz"+suffix, ntrk,ipv.pull_.z());
1113		}
1114	<	}
1114	>	}
1115		suffix.clear();
1116		}
1117		if(saventuple_) {
#	Line 1064 | Line 1120 | PVStudy::endJob() {
1120		pvtxtree_->Write();
1121		file_->Close();
1122		}
1123	<	}
1068	<
1123	>
1124
1125	+	}
1126
1127		// Match a recovertex with a simvertex
1128		// ? Should the cut be parameter dependent?
#	Line 1074 | Line 1130 | PVStudy::endJob() {
1130		bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
1131		const reco::Vertex & vrec)
1132		{
1133	<	return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1133	>	if(vrec.isFake() || !vrec.isValid()) return false;
1134	>	else
1135	>	return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1136		}
1137
1138		// Match two reconstructed vertices
#	Line 1083 | Line 1141 | bool PVStudy::matchVertex( const reco::V
1141		int zsigncut)
1142		{
1143		double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
1144	<	if(verbose_)
1087	<	cout<<"The matching criteria in z is "<<cut<<endl;
1144	>	edm::LogInfo("Debug")<<"[matchVertex] The matching criteria in z is "<<cut<<endl;
1145		return (fabs(vrec1->z()-vrec2->z())<cut);
1146		}
1147
1148
1149		bool PVStudy::isResonance(const HepMC::GenParticle * p){
1150		double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
1151	<	if(verbose_) cout << "isResonance   " << p->pdg_id() << " " << ctau << endl;
1151	>	edm::LogInfo("Debug") << "[isResonance] isResonance   " << p->pdg_id() << " " << ctau << endl;
1152		return  ctau >0 && ctau <1e-6;
1153		}
1154
#	Line 1163 | Line 1220 | void PVStudy::printSimTrks(const Handle<
1220		void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1221		stringstream ss;
1222		ss << ntk;
1223	<	string suffix = ss.str();
1224	<	if(datamode == 0 )  suffix = "_" + suffix;
1225	<	if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
1226	<	if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
1227	<	if(datamode == 3 )  suffix = "_mct_" + suffix;
1223	>	string suffix;
1224	>	if(datamode == 0 )  suffix = "_" + ss.str();
1225	>	if(datamode == 1 )  suffix = "_spl1_mct_" + ss.str();
1226	>	if(datamode == 2 )  suffix = "_spl2_mct_" + ss.str();
1227	>	if(datamode == 3 )  suffix = "_mct_" + ss.str();
1228
1229		if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
1230		// Fill histograms of res and pull of ntk
1231	<	h["deltax"+suffix]->Fill(r.x());
1232	<	h["deltay"+suffix]->Fill(r.y());
1233	<	h["deltaz"+suffix]->Fill(r.z());
1234	<	h["pullx"+suffix]->Fill(r.x()/er.x());
1235	<	h["pully"+suffix]->Fill(r.y()/er.y());
1236	<	h["pullz"+suffix]->Fill(r.z()/er.z());
1237	<	h["errPVx"+suffix]->Fill(er.x());
1238	<	h["errPVy"+suffix]->Fill(er.y());
1239	<	h["errPVz"+suffix]->Fill(er.z());
1231	>	h_others->Fill1d("deltax"+suffix, r.x());
1232	>	h_others->Fill1d("deltay"+suffix, r.y());
1233	>	h_others->Fill1d("deltaz"+suffix, r.z());
1234	>	h_others->Fill1d("pullx"+suffix, r.x()/er.x());
1235	>	h_others->Fill1d("pully"+suffix, r.y()/er.y());
1236	>	h_others->Fill1d("pullz"+suffix, r.z()/er.z());
1237	>	h_others->Fill1d("errPVx"+suffix, er.x());
1238	>	h_others->Fill1d("errPVy"+suffix, er.y());
1239	>	h_others->Fill1d("errPVz"+suffix, er.z());
1240		}
1241
1242	+
1243		PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1244		map<int,double> data;
1245		data.clear();
#	Line 1199 | Line 1257 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1257		for ( int i = 0; i < 6; ++i) {
1258		switch (i) {
1259		case 0:
1260	<	fit(h["deltax"+suffix], fpar);
1260	>	if (!h_others->ReadHisto1D("deltax"+suffix)) break;
1261	>	fit(h_others->ReadHisto1D("deltax"+suffix), fpar);
1262		data[i] = fpar[0]*cm2um;
1263		break;
1264		case 1:
1265	<	fit(h["deltay"+suffix], fpar);
1265	>	if (!h_others->ReadHisto1D("deltay"+suffix)) break;
1266	>	fit(h_others->ReadHisto1D("deltay"+suffix), fpar);
1267		data[i] = fpar[0]*cm2um;
1268		break;
1269		case 2:
1270	<	fit(h["deltaz"+suffix], fpar);
1270	>	if (!h_others->ReadHisto1D("deltaz"+suffix)) break;
1271	>	fit(h_others->ReadHisto1D("deltaz"+suffix), fpar);
1272		data[i] = fpar[0]*cm2um;
1273		break;
1274		case 3:
1275	<	fit(h["pullx"+suffix], fpar);
1275	>	if (!h_others->ReadHisto1D("pullx"+suffix)) break;
1276	>	fit(h_others->ReadHisto1D("pullx"+suffix), fpar);
1277		data[i] = fpar[0];
1278		break;
1279		case 4:
1280	<	fit(h["pully"+suffix], fpar);
1280	>	if (!h_others->ReadHisto1D("pully"+suffix)) break;
1281	>	fit(h_others->ReadHisto1D("pully"+suffix), fpar);
1282		data[i] = fpar[0];
1283		break;
1284		case 5:
1285	<	fit(h["pullz"+suffix], fpar);
1285	>	if (!h_others->ReadHisto1D("pullz"+suffix)) break;
1286	>	fit(h_others->ReadHisto1D("pullz"+suffix), fpar);
1287		data[i] = fpar[0];
1288		break;
1289		}
1290	<	if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1290	>	if (data[i] > 0) edm::LogInfo("Analysis") << "[Analysis] data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
1291		}
1292		}
1293		else
#	Line 1238 | Line 1302 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1302		ntk);
1303		}
1304
1305	<	void PVStudy::fit(TH1 *&hdist, double fitPars[]){
1305	>
1306	>	void PVStudy::fit(TH1 *hdist, double fitPars[]){
1307		TAxis *axis0 = hdist->GetXaxis();
1308		int nbin = axis0->GetLast();
1309		double nOF = hdist->GetBinContent(nbin+1);
1310		double nUF = hdist->GetBinContent(0);
1311		//   double fitRange = axis0->GetBinUpEdge(nbin);
1312	<	double fitRange = axis0->GetXmax();
1312	>	// double fitRange = axis0->GetXmax();
1313	>	double fitRange = 2*hdist->GetRMS();
1314		double sigMax[2] = {0.,0.};
1315
1316	<	if ( verbose_ ){
1317	<	cout << "Last bin = " << nbin;
1318	<	cout << "; hdist: Overflow Entries = " << nOF;
1319	<	cout << "; Underflow Entries = " << nUF;
1320	<	cout << "; hdist->GetEntries() = " << hdist->GetEntries();
1321	<	cout << "; fitting range = " << -fitRange << " to " << fitRange << endl;
1256	<	}
1316	>	edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
1317	>	<< "; hdist: Overflow Entries = " << nOF
1318	>	<< "; Underflow Entries = " << nUF
1319	>	<< "; hdist->GetEntries() = " << hdist->GetEntries()
1320	>	<< "; fitting range = " << -fitRange << " to " << fitRange << endl;
1321	>
1322		if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit
1323		for (int bi = 0; bi < nbin; bi++) {
1324		if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) {
#	Line 1265 | Line 1330 | void PVStudy::fit(TH1 *&hdist, double fi
1330		if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
1331		}
1332		}
1333	<	if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1333	>	//edm::LogInfo("Analysis") << "[Fit] Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
1334
1335		TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1336		fgaus->SetParameter(1, 0.);
1337	<	fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
1337	>	fgaus->SetParLimits(1, -fitRange/10., fitRange/10.);
1338		fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
1339		fgaus->SetLineColor(4);
1340	<	hdist->Fit(fgaus,"Q");
1340	>	hdist->Fit(fgaus,"QLRM");
1341		gPad->Update();
1342		TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
1343		s->SetOptStat(1111111);
1344		s->SetOptFit(0111);
1345		gPad->Update();
1346	<	if (verbose_) cout << "got function" << endl;
1346	>	//edm::LogInfo("Analysis") << "[Fit] got function" << endl;
1347		fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
1348		}
1349		else
1350		fitPars[0] = -999;
1351		}
1352
1353	<	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype) {
1353	>
1354	>	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs) {
1355		string suffix;
1356		suffix = ""; // for unsplit tracks
1357		if(datatype == 1) suffix = "1_spl"; // for splittracks 1
1358		if(datatype == 2) suffix = "2_spl"; // for splittracks 2
1359
1360	<	h["nTrk"+suffix]->Fill(trackColl->size());
1360	>	h_pvtrk->Fill1d("nTrk"+suffix, trackColl->size());
1361		for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
1362		itTrack != trackColl->end();
1363		++itTrack) {
1364	<	h["trkPt"+suffix]->Fill(itTrack->pt());
1365	<	h["trkDxy"+suffix]->Fill(itTrack->dxy());
1366	<	h["trkDz"+suffix]->Fill(itTrack->dz());
1367	<	h["trkEta"+suffix]->Fill(itTrack->eta());
1368	<	h["trkPhi"+suffix]->Fill(itTrack->phi());
1364	>	h_pvtrk->Fill1d("trkPt"+suffix, itTrack->pt());
1365	>	h_pvtrk->Fill1d("trkDxy"+suffix, itTrack->dxy());
1366	>	h_pvtrk->Fill1d("trkDxyCorr"+suffix, itTrack->dxy(bs));
1367	>	h_pvtrk->Fill1d("trkDz"+suffix, itTrack->dz());
1368	>	h_pvtrk->Fill1d("trkEta"+suffix, itTrack->eta());
1369	>	h_pvtrk->Fill1d("trkPhi"+suffix, itTrack->phi());
1370		}
1371		}
1372
1373	<	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
1373	>	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
1374		string suffix;
1375		suffix = ""; // for unsplit tracks
1376		if(datatype == 1) suffix = "1_spl"; // for splittracks 1
#	Line 1311 | Line 1378 | void PVStudy::fillTrackHistoInPV(const r
1378		int ivtx = 0;
1379		for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1380		v!=vertexColl->end(); ++v, ++ivtx) {
1381	+	if(v->isFake()) continue;
1382		if(fillNtuple) {
1383		if(datatype == 0) {
1384		nTrkPV_[ivtx] = v->tracksSize();
#	Line 1323 | Line 1391 | void PVStudy::fillTrackHistoInPV(const r
1391		recx_err_[ivtx] = v->xError();
1392		recy_err_[ivtx] = v->yError();
1393		recz_err_[ivtx] = v->zError();
1326	–
1327	–
1394		}
1395		if(datatype == 1) {
1396		nTrkPV1_spl_[ivtx] = v->tracksSize();
#	Line 1351 | Line 1417 | void PVStudy::fillTrackHistoInPV(const r
1417		recy2_err_spl_[ivtx] = v->yError();
1418		recz2_err_spl_[ivtx] = v->zError();
1419		}
1420	+	if(datatype == 4) { // for pixelVertices
1421	+	nTrkPV_pxlpvtx_[ivtx] = v->tracksSize();
1422	+	sumptsq_pxlpvtx_[ivtx] = sumPtSquared(*v);
1423	+	isValid_pxlpvtx_[ivtx] = v->isValid();
1424	+	isFake_pxlpvtx_[ivtx] = v->isFake();
1425	+	recx_pxlpvtx_[ivtx] = v->x();
1426	+	recy_pxlpvtx_[ivtx] = v->y();
1427	+	recz_pxlpvtx_[ivtx] = v->z();
1428	+	recx_err_pxlpvtx_[ivtx] = v->xError();
1429	+	recy_err_pxlpvtx_[ivtx] = v->yError();
1430	+	recz_err_pxlpvtx_[ivtx] = v->zError();
1431	+	}
1432		}
1433		}
1434		// For histogram only fill the leading pvtx parameters
1435		if(fillHisto) {
1436	<	h["nTrkPV"+suffix]->Fill(vertexColl->begin()->tracksSize());
1436	>	h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1437	>	h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1438	>	int nHWTrkPV_ = 0;
1439	>
1440		try {
1441		for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1442		t!=vertexColl->begin()->tracks_end(); t++) {
1443		// illegal charge
1444		if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
1445	<	//      h["trkPtPV"]->Fill(0.);
1445	>	//      h_pvtrk->Fill1d("trkPtPV", 0.);
1446		}
1447		else {
1448	<	h["trkPtPV"+suffix]->Fill((**t).pt());
1449	<	h["trkDxyPV"+suffix]->Fill((**t).dxy());
1450	<	h["trkDzPV"+suffix]->Fill((**t).dz());
1451	<	h["trkEtaPV"+suffix]->Fill((**t).eta());
1452	<	h["trkPhiPV"+suffix]->Fill((**t).phi());
1448	>	h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1449	>	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1450	>	h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
1451	>	h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1452	>	h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1453	>	h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1454	>
1455	>	if(vertexColl->begin()->trackWeight(*t) > 0.5)
1456	>	nHWTrkPV_++;
1457		}
1458		}
1459		}
1460		catch (...) {
1461		// exception thrown when trying to use linked track
1462	<	//          h["trkPtPV"]->Fill(0.);
1462	>	//          h_pvtrk->Fill1d("trkPtPV", 0.);
1463		}
1464	+	h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
1465		}
1466	+
1467
1468		}
1469
#	Line 1406 | Line 1493 | void PVStudy::fillMCHisto(std::vector<si
1493		int nrectrk = vrec->tracksSize();
1494		vsim->recVtx=NULL;
1495
1496	<	if(verbose_){
1497	<	cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1498	<	cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1499	<	}
1413	<	if ( matchVertex(*vsim,*vrec) ) {
1496	>	edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1497	>	edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1498	>
1499	>	if ( matchVertex(*vsim,*vrec)  ) {
1500		vsim->recVtx=&(*vrec);
1501
1502		// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
#	Line 1418 | Line 1504 | void PVStudy::fillMCHisto(std::vector<si
1504		//|| (!vsim->recVtx) )
1505		//vsim->recVtx=&(*vrec);
1506		//}
1507	<	if(verbose_)
1422	<	cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1507	>	edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1508
1509		double fres_mct[3];
1510		double ferror_mct[3];
#	Line 1431 | Line 1516 | void PVStudy::fillMCHisto(std::vector<si
1516		ferror_mct[1] = vsim->recVtx->yError();
1517		ferror_mct[2] = vsim->recVtx->zError();
1518
1519	<	h["deltax"+suffix]->Fill( fres_mct[0] );
1520	<	h["deltay"+suffix]->Fill( fres_mct[1] );
1521	<	h["deltaz"+suffix]->Fill( fres_mct[2] );
1522	<	h["pullx"+suffix]->Fill( fres_mct[0]/ferror_mct[0] );
1523	<	h["pully"+suffix]->Fill( fres_mct[1]/ferror_mct[1] );
1524	<	h["pullz"+suffix]->Fill( fres_mct[2]/ferror_mct[2] );
1525	<	h["errPVx"+suffix]->Fill( ferror_mct[0] );
1526	<	h["errPVy"+suffix]->Fill( ferror_mct[1] );
1527	<	h["errPVz"+suffix]->Fill( ferror_mct[2] );
1519	>	h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1520	>	h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1521	>	h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1522	>	h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1523	>	h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1524	>	h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1525	>	h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1526	>	h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1527	>	h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1528		pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1529		error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1530		nrectrk));
#	Line 1484 | Line 1569 | void PVStudy::fillMCHisto(std::vector<si
1569		} // End of  if(saventuple_) {
1570		} //  if ( matchVertex(*vsim,*vrec) ) {
1571		else {  // no rec vertex found for this simvertex
1572	<	if(verbose_)
1488	<	cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1572	>	edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1573		}
1574		}
1575		}
1576
1493	–
1577		void PVStudy::SetVarToZero() {
1578		//Greg's variables
1579		fntrk_ = 0;
1580		//pvtx position (x,y,z) residual and error
1581		for(int i = 0; i<3;i++)    {
1582	<	fres_[i] = 0;
1583	<	ferror_[i] = 0;
1582	>	fres_[i] = 0;
1583	>	ferror_[i] = 0;
1584		}
1585
1586		//Yanyan's variables
#	Line 1531 | Line 1614 | void PVStudy::SetVarToZero() {
1614		recy_err_[i] = 0;
1615		recz_err_[i] = 0;
1616
1534	–
1535	–
1617		// recoVertices with splitTrack1
1618		nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx
1619		sumptsq1_spl_[i] = 0;
#	Line 1545 | Line 1626 | void PVStudy::SetVarToZero() {
1626		recy1_err_spl_[i] = 0;
1627		recz1_err_spl_[i] = 0;
1628
1548	–
1629		// recoVertices with splitTrack2
1630		nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx
1631		sumptsq2_spl_[i] = 0;
#	Line 1558 | Line 1638 | void PVStudy::SetVarToZero() {
1638		recy2_err_spl_[i] = 0;
1639		recz2_err_spl_[i] = 0;
1640
1641	+	//pixelVertices
1642	+	nTrkPV_pxlpvtx_[i] = 0; // Number of tracks in the pvtx
1643	+	sumptsq_pxlpvtx_[i] = 0;
1644	+	isValid_pxlpvtx_[i] = -1;
1645	+	isFake_pxlpvtx_[i] = -1;
1646	+	recx_pxlpvtx_[i] = 0;
1647	+	recy_pxlpvtx_[i] = 0;
1648	+	recz_pxlpvtx_[i] = 0;
1649	+	recx_err_pxlpvtx_[i] = 0;
1650	+	recy_err_pxlpvtx_[i] = 0;
1651	+	recz_err_pxlpvtx_[i] = 0;
1652	+
1653		// matched two-vertices
1654		nTrkPV1_spl_twovtx_[i] = 0;
1655		nTrkPV2_spl_twovtx_[i] = 0;
#	Line 1603 | Line 1695 | void PVStudy::SetVarToZero() {
1695		pully_spl2_mct_[i] = 0;
1696		pullz_spl2_mct_[i] = 0;
1697		}
1606	–
1698		}
1699
1700		double PVStudy::sumPtSquared(const reco::Vertex & v)  {
#	Line 1616 | Line 1707 | double PVStudy::sumPtSquared(const reco:
1707		return sum;
1708		}
1709
1710	+
1711	+

Diff Legend

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