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.20 by yygao, Thu Feb 18 10:59:48 2010 UTC vs.
Revision 1.23 by yygao, Mon May 31 14:19:45 2010 UTC

#	Line 38 | Line 38 | Implementation:
38		#include "DataFormats/TrackReco/interface/Track.h"
39		#include "DataFormats/TrackReco/interface/TrackFwd.h"
40		#include "FWCore/ServiceRegistry/interface/Service.h"
41	<	#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
41	>	#include "CommonTools/UtilAlgos/interface/TFileService.h"
42		#include "UserCode/PVStudy/interface/PVStudy.h"
43		//
44		#include "DataFormats/VertexReco/interface/Vertex.h"
#	Line 86 | Line 86 | PVStudy::PVStudy(const edm::ParameterSet
86		nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
87		nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
88		zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
89	+	useHWTrk_                  = iConfig.getUntrackedParameter<bool>("useHWTrk",false);
90	+	avgTrkPtInPVMin_                  = iConfig.getUntrackedParameter<double>("avgTrkPtInPVMin");
91	+	avgTrkPtInPVMax_                  = iConfig.getUntrackedParameter<double>("avgTrkPtInPVMax");
92		bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
93
94		//now do what ever initialization is needed
#	Line 129 | Line 132 | PVStudy::PVStudy(const edm::ParameterSet
132		//Fill the variables in the twovtx pair (recvtx1, recvtx2)
133		// Information related to the analyzing the two-vertex method
134
132	–
133	–
135		pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
136		pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
137		pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
#	Line 219 | Line 220 | PVStudy::PVStudy(const edm::ParameterSet
220		pvtxtree_->Branch("ndofPV_spl2_mct",&ndofPV_spl2_mct_,"ndofPV_spl2_mct[nrecPV_spl2_mct]/D");
221		pvtxtree_->Branch("normchi2PV_spl2_mct",&normchi2PV_spl2_mct_,"normchi2PV_spl2_mct[nrecPV_spl2_mct]/D");
222		pvtxtree_->Branch("avgPtPV_spl2_mct",&avgPtPV_spl2_mct_,"avgPtPV_spl2_mct[nrecPV_spl2_mct]/D");
222	–
223		}
224		}
225
#	Line 643 | Line 643 | PVStudy::analyze(const edm::Event& iEven
643		// Step 1:  Apply event cleaning for data and MC
644		//          WARNING: event selection cut are hard coded!!
645		//========================================================================
646	–	// =====Loop over the trackColl to get the fraction of HighPurity tracks
647	–	int nTracks = 0;
648	–	int nHighPurityTracks=0;
649	–	double fHighPurity=0;
650	–
651	–	for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
652	–	TrackRef tkref(trackCollectionHandle,i);
653	–	if(tkref->quality(reco::TrackBase::highPurity)) ++nHighPurityTracks;
654	–	}
655	–	if(nTracks>0)
656	–	fHighPurity =  double(nHighPurityTracks)/double(nTracks);
657	–	if(verbose_)
658	–	cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;
659	–
660	–	if( (fHighPurity<0.2 && nTracks>10) || vertexColl->begin()->isFake()) return;
646
647		glob_runno_ = iEvent.id().run();
648		glob_evtno_ = iEvent.id().event();
649		glob_ls_   = iEvent.luminosityBlock();
650		glob_bx_  = iEvent.bunchCrossing();
666	–	return;
667	–
651
652		//========================================================================
653		// Step 2: Fill histograms for the splitting consistency checks
#	Line 723 | Line 706 | PVStudy::analyze(const edm::Event& iEven
706		for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
707		v!=vertexColl->end(); ++v) {
708		if(v->isValid() && ! v->isFake()) {
726	–	// xsep
709		if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
710		min_xsep = fabs(v->x()- vertexColl->begin()->x());
729	–	// ysep
711		if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
712		min_ysep = fabs(v->y()- vertexColl->begin()->y());
732	–	// zsep
713		if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
714		min_zsep = fabs(v->z()- vertexColl->begin()->z());
735	–	// xsep_sign
715		double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
716		if(xsep_sign < min_xsep_sign)
717		min_xsep_sign = xsep_sign;
739	–	// ysep_sign
718		double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
719		if(ysep_sign < min_ysep_sign)
720		min_ysep_sign = ysep_sign;
743	–	// zsep_sign
721		double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
722		if(zsep_sign < min_zsep_sign)
723		min_zsep_sign = zsep_sign;
747	–	// ntrksep
724		if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
725		min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
750	–	// sumpt2sep
726		if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
727		min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
728		}
#	Line 803 | Line 778 | PVStudy::analyze(const edm::Event& iEven
778		if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
779		edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
780
781	<	int nTrkPV1 = recvtx1->tracksSize();
782	<	int nTrkPV2 = recvtx2->tracksSize();
781	>	int nTrkPV1, nTrkPV2;
782	>	if(useHWTrk_) {
783	>	nTrkPV1 = nHWTrkRecVtx(*recvtx1);
784	>	nTrkPV2 = nHWTrkRecVtx(*recvtx2);
785	>	}
786	>	else {
787	>	nTrkPV1 = recvtx1->tracksSize();
788	>	nTrkPV2 = recvtx2->tracksSize();
789	>	}
790		double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
791		h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
792		h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
#	Line 849 | Line 831 | PVStudy::analyze(const edm::Event& iEven
831		simy_[isimpv] = vsim->y;
832		simz_[isimpv] = vsim->y;
833		simptsq_[isimpv] = vsim->ptsq;
834	<	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
853	<	fillMCHisto(vsim, isimpv, vertexColl, 3);
834	>	fillMCHisto(vsim, isimpv, vertexColl, 3, avgTrkPtInPVMin_, avgTrkPtInPVMax_);
835		}
836		}
837
#	Line 865 | Line 846 | PVStudy::analyze(const edm::Event& iEven
846		int isimpv = 0;
847		for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
848		vsim!=simpv.end(); vsim++, isimpv++){
849	<	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
850	<	fillMCHisto(vsim, isimpv, matchedVertexColl1, 1);
870	<	fillMCHisto(vsim, isimpv, matchedVertexColl2, 2);
849	>	fillMCHisto(vsim, isimpv, matchedVertexColl1, 1, avgTrkPtInPVMin_, avgTrkPtInPVMax_);
850	>	fillMCHisto(vsim, isimpv, matchedVertexColl2, 2, avgTrkPtInPVMin_, avgTrkPtInPVMax_);
851		}
852		}
853
#	Line 882 | Line 862 | PVStudy::analyze(const edm::Event& iEven
862		nrecPV_twovtx_ = 0;
863		for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
864		v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
865	<	++v1, ++v2) {
865	>	++v1, ++v2) {
866	>
867	>	// ==================================================================
868	>	// With the option to fill the histograms at a given average pT range
869	>	// ==================================================================
870	>	if (avgPtRecVtx(*v1) < avgTrkPtInPVMin_ || avgPtRecVtx(*v1) > avgTrkPtInPVMax_ ) continue;
871	>	if (avgPtRecVtx(*v2) < avgTrkPtInPVMin_ || avgPtRecVtx(*v2) > avgTrkPtInPVMax_ ) continue;
872	>
873		h_misc->Fill1d("ndofPVDiff", v1->ndof() - v2->ndof());
874		h_misc->Fill1d("ndofPVRelDiff", (v1->ndof()-v2->ndof())/(v1->ndof()+v2->ndof()));
875		fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
#	Line 891 | Line 878 | PVStudy::analyze(const edm::Event& iEven
878		ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
879		ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
880		ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
881	<	fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
881	>
882	>	int nTrkPV1, nTrkPV2;
883	>	if(useHWTrk_) {
884	>	nTrkPV1 = nHWTrkRecVtx(*v1);
885	>	nTrkPV2 = nHWTrkRecVtx(*v2);
886	>	}
887	>	else {
888	>	nTrkPV1 = v1->tracksSize();
889	>	nTrkPV2 = v2->tracksSize();
890	>	}
891	>
892	>	fntrk_ = (nTrkPV1 + nTrkPV2)/2;
893
894		h_summary->Fill1d("deltax", fres_[0] );
895		h_summary->Fill1d("deltay", fres_[1] );
#	Line 906 | Line 904 | PVStudy::analyze(const edm::Event& iEven
904		error(ferror_[0], ferror_[1], ferror_[2]),
905		fntrk_));
906		// Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
907	+
908		fillHisto(res(fres_[0], fres_[1], fres_[2]),
909		error(ferror_[0], ferror_[1], ferror_[2]),
910		fntrk_,0);
#	Line 924 | Line 923 | PVStudy::analyze(const edm::Event& iEven
923
924
925		//SplittedVertex
926	<
928	<	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
926	>	nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = nTrkPV1;
927		ndofPV1_spl_twovtx_[nrecPV_twovtx_] = v1->ndof();
928		normchi2PV1_spl_twovtx_[nrecPV_twovtx_] = v1->normalizedChi2();
929		avgPtPV1_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v1);
930		errx1_spl_twovtx_[nrecPV_twovtx_] = v1->xError();
931		erry1_spl_twovtx_[nrecPV_twovtx_] = v1->yError();
932		errz1_spl_twovtx_[nrecPV_twovtx_] = v1->zError();
933	<
934	<	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
933	>
934	>	nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = nTrkPV2;
935		ndofPV2_spl_twovtx_[nrecPV_twovtx_] = v2->ndof();
936		normchi2PV2_spl_twovtx_[nrecPV_twovtx_] = v2->normalizedChi2();
937		avgPtPV2_spl_twovtx_[nrecPV_twovtx_] = avgPtRecVtx(*v2);
#	Line 942 | Line 940 | PVStudy::analyze(const edm::Event& iEven
940		errz2_spl_twovtx_[nrecPV_twovtx_] = v2->zError();
941
942		nrecPV_twovtx_++;
945	–
946	–
943
944		// Print some information of the two tracks events
945		if(verbose_ && fntrk_ < 4) {
#	Line 1129 | Line 1125 | void PVStudy::printSimTrks(const Handle<
1125		}
1126		}
1127
1128	<	void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1128	>	void PVStudy::fillHisto(res r, error er, int ntk, int datamode) {
1129		stringstream ss;
1130		ss << ntk;
1131		string suffix;
#	Line 1292 | Line 1288 | void PVStudy::fillTrackHistoInPV(const r
1288		v!=vertexColl->end(); ++v, ++ivtx) {
1289		if(!v->isFake()) {
1290		h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1291	<	h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1292	<	int nHWTrkPV_ = 0;
1291	>	h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1292	>	h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkRecVtx(*v));
1293		try {
1294		for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1295		t!=vertexColl->begin()->tracks_end(); t++) {
#	Line 1308 | Line 1304 | void PVStudy::fillTrackHistoInPV(const r
1304		h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1305		h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1306		h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1311	–	if(vertexColl->begin()->trackWeight(*t) > 0.5)
1312	–	nHWTrkPV_++;
1307		}
1308		}
1309		}
#	Line 1317 | Line 1311 | void PVStudy::fillTrackHistoInPV(const r
1311		// exception thrown when trying to use linked track
1312		//        h_pvtrk->Fill1d("trkPtPV", 0.);
1313		}
1320	–	h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
1314		}
1315		}
1316		}
1317
1318	<	void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
1318	>	void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl,
1319	>	int datatype, double avgTrkPtInPVMin_, double avgTrkPtInPVMax_)
1320		{
1321		std::string suffix;
1322		// Set the vertexColl and histogram suffix according to datamode
#	Line 1338 | Line 1332 | void PVStudy::fillMCHisto(std::vector<si
1332		suffix = "_mct";
1333		nrecPV_mct_ = 0;
1334		}
1335	<
1335	>
1336		//========================================================
1337	<	//  look for a matching reconstructed vertex in vertexColl
1338	<	//========================================================
1339	<
1337	>	//  For each simulated vertex, look for a match in the vertexColl
1338	>	//  If more than 1 recVtx is found, use the one with closest in Z
1339	>	//========================================================
1340	>
1341	>	// === look for a matching reconstructed vertex in vertexColl
1342		for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
1343		vrec!=vtxColl->end(); ++vrec){
1344		vsim->recVtx=NULL;
1345		edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1346		edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1347	<
1347	>
1348		if ( matchVertex(*vsim,*vrec)  ) {
1349		vsim->recVtx=&(*vrec);
1350	<
1355	<	// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1350	>	//if the matching critera are fulfilled, accept the rec-vertex that is closest in z
1351		//if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
1352		//|| (!vsim->recVtx) )
1353		//vsim->recVtx=&(*vrec);
1359	–	//}
1360	–	edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1361	–
1362	–	double fres_mct[3];
1363	–	double ferror_mct[3];
1364	–	int nrectrk =  int(vsim->recVtx->tracksSize());
1365	–	double ndofPV = vsim->recVtx->ndof();
1366	–	double normchi2PV =  vsim->recVtx->normalizedChi2();
1367	–	double avgPtPV = avgPtRecVtx(*(vsim->recVtx));
1368	–
1369	–	fres_mct[0] = vsim->recVtx->x()-vsim->x;
1370	–	fres_mct[1] = vsim->recVtx->y()-vsim->y;
1371	–	fres_mct[2] = vsim->recVtx->z()-vsim->z;
1372	–	ferror_mct[0] = vsim->recVtx->xError();
1373	–	ferror_mct[1] = vsim->recVtx->yError();
1374	–	ferror_mct[2] = vsim->recVtx->zError();
1375	–
1376	–	h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1377	–	h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1378	–	h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1379	–	h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1380	–	h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1381	–	h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1382	–	h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1383	–	h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1384	–	h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1385	–	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1386	–	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1387	–	nrectrk));
1388	–	// Fill histo according to its track multiplicity
1389	–	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1390	–	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1391	–	nrectrk, datatype);
1392	–
1393	–	if(saventuple_) {
1394	–	//Fill the values for variables in pvtxtree_
1395	–	if(datatype == 1) {
1396	–	nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1397	–	ndofPV_spl1_mct_[nrecPV_spl1_mct_] =  ndofPV;
1398	–	normchi2PV_spl1_mct_[nrecPV_spl1_mct_] =  normchi2PV;
1399	–	avgPtPV_spl1_mct_[nrecPV_spl1_mct_] = avgPtPV;
1400	–	deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1401	–	deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[1];
1402	–	deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2];
1403	–	pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
1404	–	pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1405	–	pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1406	–	errx_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[0];
1407	–	erry_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[1];
1408	–	errz_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[2];
1409	–	nrecPV_spl1_mct_++;
1410	–	}
1411	–	if(datatype == 2) {
1412	–	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1413	–	ndofPV_spl2_mct_[nrecPV_spl2_mct_] =  ndofPV;
1414	–	normchi2PV_spl2_mct_[nrecPV_spl2_mct_] =  normchi2PV;
1415	–	avgPtPV_spl2_mct_[nrecPV_spl2_mct_] = avgPtPV;
1416	–	deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1417	–	deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[1];
1418	–	deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2];
1419	–	pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
1420	–	pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1421	–	pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1422	–	errx_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[0];
1423	–	erry_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[1];
1424	–	errz_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[2];
1425	–	nrecPV_spl2_mct_++;
1426	–	}
1427	–	if(datatype == 3) {
1428	–	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1429	–	ndofPV_mct_[nrecPV_mct_] =  ndofPV;
1430	–	normchi2PV_mct_[nrecPV_mct_] =  normchi2PV;
1431	–	avgPtPV_mct_[nrecPV_mct_] = avgPtPV;
1432	–	deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1433	–	deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1434	–	deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1435	–	pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
1436	–	pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1437	–	pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1438	–	errx_mct_[nrecPV_mct_] =  ferror_mct[0];
1439	–	erry_mct_[nrecPV_mct_] =  ferror_mct[1];
1440	–	errz_mct_[nrecPV_mct_] =  ferror_mct[2];
1441	–	nrecPV_mct_++;
1442	–	}
1443	–	} // End of  if(saventuple_) {
1444	–	} //  if ( matchVertex(*vsim,*vrec) ) {
1445	–	else {  // no rec vertex found for this simvertex
1446	–	edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1354		}
1355		}
1356	+
1357	+	// === If match found fill the residual and pulls
1358	+	if(vsim->recVtx) {
1359	+	edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1360	+	double avgPtPV = avgPtRecVtx(*(vsim->recVtx));
1361	+	if(avgPtPV>avgTrkPtInPVMax_ || avgPtPV<avgTrkPtInPVMin_) return;
1362	+
1363	+	// if study the resolution/pull for all tracks in the PVtx
1364	+	int nrectrk;
1365	+	if(useHWTrk_)
1366	+	nrectrk=nHWTrkRecVtx(*(vsim->recVtx));
1367	+	else
1368	+	nrectrk =  int(vsim->recVtx->tracksSize());
1369	+
1370	+	double fres_mct[3];
1371	+	double ferror_mct[3];
1372	+
1373	+	double ndofPV = vsim->recVtx->ndof();
1374	+	double normchi2PV =  vsim->recVtx->normalizedChi2();
1375	+
1376	+	fres_mct[0] = vsim->recVtx->x()-vsim->x;
1377	+	fres_mct[1] = vsim->recVtx->y()-vsim->y;
1378	+	fres_mct[2] = vsim->recVtx->z()-vsim->z;
1379	+	ferror_mct[0] = vsim->recVtx->xError();
1380	+	ferror_mct[1] = vsim->recVtx->yError();
1381	+	ferror_mct[2] = vsim->recVtx->zError();
1382	+
1383	+	h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
1384	+	h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
1385	+	h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
1386	+	h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
1387	+	h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
1388	+	h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
1389	+	h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
1390	+	h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
1391	+	h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
1392	+	pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1393	+	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1394	+	nrectrk));
1395	+	// Fill histo according to its track multiplicity
1396	+	fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
1397	+	error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
1398	+	nrectrk, datatype);
1399	+
1400	+	if(saventuple_) {
1401	+	//Fill the values for variables in pvtxtree_
1402	+	if(datatype == 1) {
1403	+	nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
1404	+	ndofPV_spl1_mct_[nrecPV_spl1_mct_] =  ndofPV;
1405	+	normchi2PV_spl1_mct_[nrecPV_spl1_mct_] =  normchi2PV;
1406	+	avgPtPV_spl1_mct_[nrecPV_spl1_mct_] = avgPtPV;
1407	+	deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
1408	+	deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[1];
1409	+	deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2];
1410	+	pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
1411	+	pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
1412	+	pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
1413	+	errx_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[0];
1414	+	erry_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[1];
1415	+	errz_spl1_mct_[nrecPV_spl1_mct_] =  ferror_mct[2];
1416	+	nrecPV_spl1_mct_++;
1417	+	}
1418	+	if(datatype == 2) {
1419	+	nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
1420	+	ndofPV_spl2_mct_[nrecPV_spl2_mct_] =  ndofPV;
1421	+	normchi2PV_spl2_mct_[nrecPV_spl2_mct_] =  normchi2PV;
1422	+	avgPtPV_spl2_mct_[nrecPV_spl2_mct_] = avgPtPV;
1423	+	deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
1424	+	deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[1];
1425	+	deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2];
1426	+	pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
1427	+	pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
1428	+	pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
1429	+	errx_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[0];
1430	+	erry_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[1];
1431	+	errz_spl2_mct_[nrecPV_spl2_mct_] =  ferror_mct[2];
1432	+	nrecPV_spl2_mct_++;
1433	+	}
1434	+	if(datatype == 3) {
1435	+	nTrkPV_mct_[nrecPV_mct_] = nrectrk;
1436	+	ndofPV_mct_[nrecPV_mct_] =  ndofPV;
1437	+	normchi2PV_mct_[nrecPV_mct_] =  normchi2PV;
1438	+	avgPtPV_mct_[nrecPV_mct_] = avgPtPV;
1439	+	deltax_mct_[nrecPV_mct_] =  fres_mct[0];
1440	+	deltay_mct_[nrecPV_mct_] =  fres_mct[0];
1441	+	deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
1442	+	pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
1443	+	pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
1444	+	pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
1445	+	errx_mct_[nrecPV_mct_] =  ferror_mct[0];
1446	+	erry_mct_[nrecPV_mct_] =  ferror_mct[1];
1447	+	errz_mct_[nrecPV_mct_] =  ferror_mct[2];
1448	+	nrecPV_mct_++;
1449	+	}
1450	+	} // End of  if(saventuple_) {
1451	+	} //  if ( matchVertex(*vsim,*vrec) ) {
1452	+	else {  // no rec vertex found for this simvertex
1453	+	edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
1454	+	}
1455		}
1456
1457	+
1458		void PVStudy::SetVarToZero() {
1459		fntrk_ = 0;
1460		//pvtx position (x,y,z) residual and error
#	Line 1577 | Line 1584 | double PVStudy::avgPtRecVtx(const reco::
1584
1585		if(v.isFake() || !v.isValid() || v.tracksSize()==0 ) return 0;
1586		else {
1587	+	int nHWTrk = 0;
1588		double sumpT = 0.;
1589		for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1590	<	sumpT +=  (**it).pt();
1590	>	if(v.trackWeight(*it) > 0.5 ) {
1591	>	sumpT +=  (**it).pt();
1592	>	nHWTrk++;
1593	>	}
1594		}
1595	<	return sumpT/double(v.tracksSize());
1595	>	if(nHWTrk > 0)
1596	>	return sumpT/double(nHWTrk);
1597	>	else
1598	>	return 0;
1599		}
1600		}
1601
1602	<
1603	<
1602	>	int PVStudy::nHWTrkRecVtx(const reco::Vertex & v)  {
1603	>	int nHWTrkPV = 0;
1604	>	for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
1605	>	if(v.trackWeight(*it) > 0.5)
1606	>	nHWTrkPV++;
1607	>	}
1608	>	return nHWTrkPV;
1609	>	}

Diff Legend

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