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.14 by jengbou, Fri Dec 4 20:48:04 2009 UTC vs.
Revision 1.18 by jengbou, Fri Jan 29 21:36:22 2010 UTC

#	Line 12 | Line 12 | Implementation:
12		*/
13		//
14		// Original Author:  "Geng-yuan Jeng/UC Riverside"
15	<	//                   "Yanyan Gao/Fermilab ygao@fnal.gov"
15	>	//                    "Yanyan Gao/Fermilab ygao@fnal.gov"
16		//         Created:  Thu Aug 20 11:55:40 CDT 2009
17		// $Id$
18		//
#	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	+
55
56		//root
57		#include <TROOT.h>
#	Line 60 | Line 62 | Implementation:
62		#include <TPad.h>
63
64		using namespace std;
65	+	typedef math::XYZTLorentzVectorF LorentzVector;
66	+	typedef math::XYZPoint Point;
67
68		PVStudy::PVStudy(const edm::ParameterSet& iConfig)
69		{
#	Line 72 | Line 76 | PVStudy::PVStudy(const edm::ParameterSet
76		splitTrackCollection2Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
77		vertexCollectionTag_       = iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection");
78		splitVertexCollection1Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
79	<	splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
76	<	pixelVertexCollectionTag_  = iConfig.getParameter<edm::InputTag>("pixelVertexCollectionTag");
79	>	splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
80		verbose_                   = iConfig.getUntrackedParameter<bool>("verbose",false);
81		realData_                  = iConfig.getUntrackedParameter<bool>("realData",false);
82		analyze_                   = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
83		saventuple_                = iConfig.getUntrackedParameter<bool>("saventuple",false);
84		outputfilename_            = iConfig.getUntrackedParameter<string>("OutputFileName");
85		ntrkdiffcut_               = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
83	–	zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
86		nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
87		nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
88	+	zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
89	+	bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
90	+	//   applyEvtClean_             = iConfig.getUntrackedParameter<bool>("applyEvtClean",false);
91		histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
92	<
92	>
93		//now do what ever initialization is needed
94		pvinfo.clear();
95
91	–
96		// Specify the data mode vector
97		if(realData_) datamode.push_back(0);
98		else {
#	Line 97 | Line 101 | PVStudy::PVStudy(const edm::ParameterSet
101		datamode.push_back(2);
102		datamode.push_back(3);
103		}
104	<
101	<	// Create ntuple files if needed
104	>	// Create ntuple files if needed
105		if(saventuple_) {
106		file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
107		ftree_ = new TTree("mytree","mytree");
#	Line 217 | Line 220 | PVStudy::PVStudy(const edm::ParameterSet
220		pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
221		}
222		}
223	<
223	>
224		setRootStyle();
225
226	+	//========================================
227	+	// Booking histograms
228	+	//========================================
229	+
230		// Create a root file for histograms
231		theFile = new TFile(histoFileName_.c_str(), "RECREATE");
232		// make diretories
#	Line 227 | Line 234 | PVStudy::PVStudy(const edm::ParameterSet
234		theFile->cd();
235		theFile->mkdir("Others");
236		theFile->cd();
237	+	if (analyze_) {
238	+	theFile->mkdir("Results");
239	+	theFile->cd();
240	+	}
241
242		// Book Histograms:
243		h_pvtrk   = new PVHistograms();
#	Line 240 | Line 251 | PVStudy::PVStudy(const edm::ParameterSet
251		else {
252		stringstream ss;
253		ss << i;
254	<	h_pvtrk->Init("pvTrk",ss.str(),"spl");
254	>	h_pvtrk->Init("pvTrk", ss.str(),"spl");
255		}
256		}
257		h_pixvtx->Init("pixVtx");
#	Line 280 | Line 291 | PVStudy::PVStudy(const edm::ParameterSet
291		if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
292		suffix.clear();
293		} // End of Book histograms sensitive to data/mc
294	+
295
296		}
297
298		PVStudy::~PVStudy()
299		{
300	+
301	+	// do anything here that needs to be done at desctruction time
302	+	// (e.g. close files, deallocate resources etc.)
303		theFile->cd();
304		theFile->cd("Summary");
305		h_pvtrk->Save();
306		h_pixvtx->Save();
292	–	h_gen->Save();
307		h_misc->Save();
308		h_summary->Save();
309	<	if (analyze_)
309	>	if (!realData_)
310	>	h_gen->Save();
311	>	if (analyze_) {
312	>	theFile->cd();
313	>	theFile->cd("Results");
314		h_ana->Save();
315	+	}
316		theFile->cd();
317		theFile->cd("Others");
318		h_others->Save();
#	Line 335 | Line 354 | void PVStudy::setRootStyle() {
354		//
355		std::vector<PVStudy::simPrimaryVertex> PVStudy::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
356		{
357	<	std::vector<PVStudy::simPrimaryVertex> simpv;
357	>	std::vector<PVStudy::simPrimaryVertex> simpv;
358		const HepMC::GenEvent* evt=evtMC->GetEvent();
359		if (evt) {
360		edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
#	Line 517 | Line 536 | PVStudy::analyze(const edm::Event& iEven
536		// Initialize Root-tuple variables if needed
537		//=======================================================
538		//if(saventuple_)
539	<	SetVarToZero();
539	>	SetVarToZero();
540
541		//=======================================================
542		// Track accessors
#	Line 532 | Line 551 | PVStudy::analyze(const edm::Event& iEven
551		} else {
552		edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
553		}
535	–
554		//splitTrackCollection1
555		static const reco::TrackCollection s_empty_splitTrackColl1;
556		const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
#	Line 543 | Line 561 | PVStudy::analyze(const edm::Event& iEven
561		} else {
562		edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
563		}
546	–
564		//splitTrackCollection2
565		static const reco::TrackCollection s_empty_splitTrackColl2;
566		const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
#	Line 552 | Line 569 | PVStudy::analyze(const edm::Event& iEven
569		if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
570		splitTrackColl2 = splitTrackCollection2Handle.product();
571		} else {
572	<	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
572	>	edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
573		}
574	<
558	<
559	<	//=======================================================
560	<	// Fill trackparameters of the input tracks to pvtx fitter
561	<	//=======================================================
562	<	edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
563	<	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype)
564	<	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
565	<	fillTrackHisto(trackColl, 0);
566	<	fillTrackHisto(splitTrackColl1, 1);
567	<	fillTrackHisto(splitTrackColl2, 2);
568	<	edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
569	<
574	>
575		//=======================================================
576		// PVTX accessors
577		//=======================================================
#	Line 578 | Line 583 | PVStudy::analyze(const edm::Event& iEven
583		if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
584		vertexColl = vertexCollectionHandle.product();
585		} else {
586	<	edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
586	>	edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
587		}
583	–
588		//splitVertexCollection1
589		static const reco::VertexCollection s_empty_splitVertexColl1;
590		const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
#	Line 590 | Line 594 | PVStudy::analyze(const edm::Event& iEven
594		splitVertexColl1 = splitVertexCollection1Handle.product();
595		} else {
596		edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
597	<	}
594	<
597	>	}
598		//splitVertexCollection2
599		static const reco::VertexCollection s_empty_splitVertexColl2;
600		const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
#	Line 602 | Line 605 | PVStudy::analyze(const edm::Event& iEven
605		} else {
606		edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
607		}
608	+
609	+
610	+	//=======================================================
611	+	// GET pixelVertices
612	+	//=======================================================
613	+	static const reco::VertexCollection s_empty_pixelVertexColl;
614	+	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
615	+	edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
616	+	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
617	+	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
618	+	pixelVertexColl = pixelVertexCollectionHandle.product();
619	+	} else {
620	+	edm::LogInfo("Debug") << "[PVStudy] pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
621	+	}
622	+
623	+	//=======================================================
624	+	// BeamSpot accessors
625	+	//=======================================================
626	+
627	+	edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
628	+	iEvent.getByLabel(bsSrc,recoBeamSpotHandle);
629	+	reco::BeamSpot bs = *recoBeamSpotHandle;
630	+	const Point beamSpot = recoBeamSpotHandle.isValid() ? Point(recoBeamSpotHandle->x0(), recoBeamSpotHandle->y0(), recoBeamSpotHandle->z0()) : Point(0, 0, 0);
631
632	<	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track and vertex collections"<<endl;
632	>	edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex and pixelvertices collections"<<endl;
633
634		//=======================================================
635		// MC simvtx accessor
#	Line 624 | Line 650 | PVStudy::analyze(const edm::Event& iEven
650		}catch(const Exception&){
651		edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
652		}
653	<
653	>
654		//=======================================================
655	<	// GET pixelVertices
655	>	// Apply event cleaning for firstcoll data and MC
656		//=======================================================
657	<	static const reco::VertexCollection s_empty_pixelVertexColl;
658	<	const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
659	<	edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
660	<	iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
661	<	if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
662	<	pixelVertexColl = pixelVertexCollectionHandle.product();
663	<	} else {
664	<	edm::LogInfo("Debug") << "[PVStudy] pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
657	>
658	>	// =====Loop over the trackColl to tet the fraction of HighPurity tracks
659	>	int nTracks = 0;
660	>	int nHighPurityTracks=0;
661	>	double fHighPurity=0;
662	>
663	>	for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
664	>	TrackRef tkref(trackCollectionHandle,i);
665	>	if(tkref->quality(reco::TrackBase::highPurity)) ++nHighPurityTracks;
666	>	}
667	>	if(nTracks>0)
668	>	fHighPurity =  double(nHighPurityTracks)/double(nTracks);
669	>	if(verbose_)
670	>	cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;
671	>
672	>	if( (fHighPurity<0.2 && nTracks>10) || vertexColl->begin()->isFake()) {
673	>	//  glob_runno_ = iEvent.id().run();
674	>	//  glob_evtno_ = iEvent.id().event();
675	>	//  glob_ls_   = iEvent.luminosityBlock();
676	>	//  glob_bx_  = iEvent.bunchCrossing();
677	>	return;
678		}
679
680		//=======================================================
681	+	// Fill trackparameters of the input tracks to pvtx fitter
682	+	//=======================================================
683	+	edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
684	+	//fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs)
685	+	// datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
686	+	fillTrackHisto(trackColl, 0, beamSpot);
687	+	fillTrackHisto(splitTrackColl1, 1, beamSpot);
688	+	fillTrackHisto(splitTrackColl2, 2, beamSpot);
689	+	edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
690	+
691	+
692	+	//=======================================================
693		// Fill pixelVertices related histograms
694		//=======================================================
695		nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
696		if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
697	<	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
698	<	fillTrackHistoInPV(pixelVertexColl, 4, false, true);
697	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
698	>	fillTrackHistoInPV(pixelVertexColl, 4, false, true, beamSpot);
699		h_pixvtx->Fill1d("dzErr_pxlpvtx", pixelVertexColl->begin()->zError());
700		// Get the dZ error of the tracks in the leading pixelVertexColl
701		for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
#	Line 696 | Line 747 | PVStudy::analyze(const edm::Event& iEven
747		h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
748		h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
749
750	<	//=================================================================
750	>
751	>	//=================================================================
752		// Fill track parameter ntuple/hist for tracks used in recoVertices
753		//=================================================================
754
755	<	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
755	>	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
756		if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
757	<	fillTrackHistoInPV(vertexColl, 0, true, true);
757	>	fillTrackHistoInPV(vertexColl, 0, true, true, beamSpot);
758
759		if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
760	<	fillTrackHistoInPV(splitVertexColl1, 1, false, true);
709	<
710	<	if(splitVertexColl1->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
711	<	fillTrackHistoInPV(splitVertexColl2, 2, false, true);
760	>	fillTrackHistoInPV(splitVertexColl1, 1, false, true, beamSpot);
761
762	+	if(splitVertexColl2->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
763	+	fillTrackHistoInPV(splitVertexColl2, 2, false, true, beamSpot);
764
765		//=======================================================
766		// Compare offlinePrimaryVertices with pixelVertices
#	Line 725 | Line 776 | PVStudy::analyze(const edm::Event& iEven
776		h_pixvtx->Fill1d("reczPV_minus_reczPxlPV", recz_[0] - recz_pxlpvtx_[0]);
777		}
778
779	<
780	<
730	<	//==========================================================
779	>
780	>	//==========================================================
781		// Fill secondary/primary min separations for multi-vertices
782		//==========================================================
783
#	Line 794 | Line 844 | PVStudy::analyze(const edm::Event& iEven
844		} // End of  if(vertexColl->size()>1) {
845
846		edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
847	+
848
849	+
850		//=======================================================
851		//           PrimaryVertex Matching
852		// 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
#	Line 802 | Line 854 | PVStudy::analyze(const edm::Event& iEven
854		// Assume the first match is the primary vertex,
855		// since vertexColl are sorted in decreasing order of sum(pT)
856		//=======================================================
857	<
857	>
858		edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
859		reco::VertexCollection s_empty_matchedVertexColl1;
860		reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
#	Line 825 | Line 877 | PVStudy::analyze(const edm::Event& iEven
877		edm::LogInfo("Debug")<<"[PVStudy] Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
878		edm::LogInfo("Debug")<<"[PVStudy] recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
879		edm::LogInfo("Debug")<<"[PVStudy] recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
880	<
880	>
881	>	double twovtxsig = (recvtx1->z()-recvtx2->z())/max(recvtx1->zError(), recvtx2->zError());
882	>	h_misc->Fill1d("twovtxzsign", twovtxsig);
883		if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
884		edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
885
#	Line 862 | Line 916 | PVStudy::analyze(const edm::Event& iEven
916		// Compare the reconstructed vertex position and calculate resolution/pulls
917		if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
918		//fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
919	<	fillTrackHistoInPV(matchedVertexColl1, 1, true, false);
920	<	fillTrackHistoInPV(matchedVertexColl2, 2, true, false);
919	>	fillTrackHistoInPV(matchedVertexColl1, 1, true, false, beamSpot);
920	>	fillTrackHistoInPV(matchedVertexColl2, 2, true, false, beamSpot);
921
922		reco::VertexCollection::const_iterator v1;
923		reco::VertexCollection::const_iterator v2;
924		nrecPV_twovtx_ = 0;
925		for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
926		v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
927	<	++v1, ++v2) {
927	>	++v1, ++v2) {
928	>	h_misc->Fill1d("ndofPVDiff", v1->ndof() - v2->ndof());
929	>	h_misc->Fill1d("ndofPVRelDiff", (v1->ndof()-v2->ndof())/(v1->ndof()+v2->ndof()));
930		fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
931		fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
932		fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
#	Line 912 | Line 968 | PVStudy::analyze(const edm::Event& iEven
968		pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];
969		nrecPV_twovtx_++;
970		} // End of analyzing res/pull
971	+
972	+	//=======================================================
973	+	// Fill simulated vertices
974	+	//=======================================================
975	+	if (!realData_) {
976	+	bool MC=false;
977	+	Handle<HepMCProduct> evtMC;
978	+	iEvent.getByLabel("generator",evtMC);
979	+	if (!evtMC.isValid()) {
980	+	MC=false;
981	+	edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
982	+	} else {
983	+	edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
984	+	MC=true;
985	+	}
986	+
987	+	if(MC){
988	+	// make a list of primary vertices:
989	+	std::vector<simPrimaryVertex> simpv;
990	+	simpv=getSimPVs(evtMC,"");
991	+	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
992	+	h_gen->Fill1d("nsimPV", simpv.size());
993	+	nsimPV_ = simpv.size();
994	+
995	+	int isimpv = 0;
996	+	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
997	+	vsim!=simpv.end(); vsim++, isimpv++){
998	+	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
999	+	simx_[isimpv] = vsim->x;
1000	+	simy_[isimpv] = vsim->y;
1001	+	simz_[isimpv] = vsim->y;
1002	+	simptsq_[isimpv] = vsim->ptsq;
1003	+
1004	+	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
1005	+	fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
1006	+	fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
1007	+	fillMCHisto(vsim, isimpv, vertexColl, 3);
1008	+	}
1009	+	} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
1010	+	} // End of   if (!realData_) {
1011	+
1012		} // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
1013		else
1014		edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
1015
1016		edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
1017
1018	<	//=======================================================
922	<	// Fill simulated vertices
923	<	//=======================================================
924	<	if (!realData_) {
925	<	bool MC=false;
926	<	Handle<HepMCProduct> evtMC;
927	<	iEvent.getByLabel("generator",evtMC);
928	<	if (!evtMC.isValid()) {
929	<	MC=false;
930	<	edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
931	<	} else {
932	<	edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
933	<	MC=true;
934	<	}
935	<
936	<	if(MC){
937	<	// make a list of primary vertices:
938	<	std::vector<simPrimaryVertex> simpv;
939	<	simpv=getSimPVs(evtMC,"");
940	<	//     simpv=getSimPVs(evtMC, simVtxs, simTrks);
941	<	h_gen->Fill1d("nsimPV", simpv.size());
942	<	nsimPV_ = simpv.size();
943	<
944	<	int isimpv = 0;
945	<	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
946	<	vsim!=simpv.end(); vsim++, isimpv++){
947	<	nsimTrkPV_[isimpv]  =vsim->nGenTrk;
948	<	simx_[isimpv] = vsim->x;
949	<	simy_[isimpv] = vsim->y;
950	<	simz_[isimpv] = vsim->y;
951	<	simptsq_[isimpv] = vsim->ptsq;
952	<
953	<	//fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
954	<	fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
955	<	fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
956	<	fillMCHisto(vsim, isimpv, vertexColl, 3);
957	<	}
958	<	} // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
959	<	} // End of   if (!realData_) {
1018	>
1019
1020		// Finally fill the ftree_
1021		if(saventuple_) {
1022		ftree_->Fill();
1023		pvtxtree_->Fill();
1024		}
1025	+
1026		}
1027
1028	+
1029	+
1030		// ------------ method called once each job just before starting event loop  ------------
1031		void
1032		PVStudy::beginJob()
#	Line 1010 | Line 1072 | PVStudy::endJob() {
1072		pvtxtree_->Write();
1073		file_->Close();
1074		}
1075	<	}
1014	<
1075	>
1076
1077	+	}
1078
1079		// Match a recovertex with a simvertex
1080		// ? Should the cut be parameter dependent?
#	Line 1020 | Line 1082 | PVStudy::endJob() {
1082		bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
1083		const reco::Vertex & vrec)
1084		{
1085	<	return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1085	>	if(vrec.isFake() || !vrec.isValid()) return false;
1086	>	else
1087	>	return (fabs(vsim.z-vrec.z())<0.0500); // =500um
1088		}
1089
1090		// Match two reconstructed vertices
#	Line 1127 | Line 1191 | void PVStudy::fillHisto(res r, error er,
1191		h_others->Fill1d("errPVz"+suffix, er.z());
1192		}
1193
1194	+
1195		PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
1196		map<int,double> data;
1197		data.clear();
#	Line 1189 | Line 1254 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1254		ntk);
1255		}
1256
1257	+
1258		void PVStudy::fit(TH1 *hdist, double fitPars[]){
1259		TAxis *axis0 = hdist->GetXaxis();
1260		int nbin = axis0->GetLast();
1261		double nOF = hdist->GetBinContent(nbin+1);
1262		double nUF = hdist->GetBinContent(0);
1263		//   double fitRange = axis0->GetBinUpEdge(nbin);
1264	<	double fitRange = axis0->GetXmax();
1264	>	// double fitRange = axis0->GetXmax();
1265	>	double fitRange = 2*hdist->GetRMS();
1266		double sigMax[2] = {0.,0.};
1267
1268		edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
#	Line 1219 | Line 1286 | void PVStudy::fit(TH1 *hdist, double fit
1286
1287		TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
1288		fgaus->SetParameter(1, 0.);
1289	<	fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
1289	>	fgaus->SetParLimits(1, -fitRange/10., fitRange/10.);
1290		fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
1291		fgaus->SetLineColor(4);
1292	<	hdist->Fit(fgaus,"Q");
1292	>	hdist->Fit(fgaus,"QLRM");
1293		gPad->Update();
1294		TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
1295		s->SetOptStat(1111111);
#	Line 1235 | Line 1302 | void PVStudy::fit(TH1 *hdist, double fit
1302		fitPars[0] = -999;
1303		}
1304
1305	<	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype) {
1305	>
1306	>	void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs) {
1307		string suffix;
1308		suffix = ""; // for unsplit tracks
1309		if(datatype == 1) suffix = "1_spl"; // for splittracks 1
#	Line 1247 | Line 1315 | void PVStudy::fillTrackHisto(const reco:
1315		++itTrack) {
1316		h_pvtrk->Fill1d("trkPt"+suffix, itTrack->pt());
1317		h_pvtrk->Fill1d("trkDxy"+suffix, itTrack->dxy());
1318	+	h_pvtrk->Fill1d("trkDxyCorr"+suffix, itTrack->dxy(bs));
1319		h_pvtrk->Fill1d("trkDz"+suffix, itTrack->dz());
1320		h_pvtrk->Fill1d("trkEta"+suffix, itTrack->eta());
1321		h_pvtrk->Fill1d("trkPhi"+suffix, itTrack->phi());
1322		}
1323		}
1324
1325	<	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
1325	>	void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
1326		string suffix;
1327		suffix = ""; // for unsplit tracks
1328		if(datatype == 1) suffix = "1_spl"; // for splittracks 1
#	Line 1261 | Line 1330 | void PVStudy::fillTrackHistoInPV(const r
1330		int ivtx = 0;
1331		for(reco::VertexCollection::const_iterator v=vertexColl->begin();
1332		v!=vertexColl->end(); ++v, ++ivtx) {
1333	+	if(v->isFake()) continue;
1334		if(fillNtuple) {
1335		if(datatype == 0) {
1336		nTrkPV_[ivtx] = v->tracksSize();
#	Line 1316 | Line 1386 | void PVStudy::fillTrackHistoInPV(const r
1386		// For histogram only fill the leading pvtx parameters
1387		if(fillHisto) {
1388		h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
1389	+	h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
1390	+	int nHWTrkPV_ = 0;
1391	+
1392		try {
1393		for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
1394		t!=vertexColl->begin()->tracks_end(); t++) {
#	Line 1325 | Line 1398 | void PVStudy::fillTrackHistoInPV(const r
1398		}
1399		else {
1400		h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1401	<	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1401	>	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
1402	>	h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
1403		h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
1404		h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
1405		h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
1406	+
1407	+	if(vertexColl->begin()->trackWeight(*t) > 0.5)
1408	+	nHWTrkPV_++;
1409		}
1410		}
1411		}
#	Line 1336 | Line 1413 | void PVStudy::fillTrackHistoInPV(const r
1413		// exception thrown when trying to use linked track
1414		//          h_pvtrk->Fill1d("trkPtPV", 0.);
1415		}
1416	+	h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
1417		}
1418	+
1419
1420		}
1421
#	Line 1369 | Line 1448 | void PVStudy::fillMCHisto(std::vector<si
1448		edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
1449		edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
1450
1451	<	if ( matchVertex(*vsim,*vrec) && vrec->isValid() && !vrec->isFake() ) {
1451	>	if ( matchVertex(*vsim,*vrec)  ) {
1452		vsim->recVtx=&(*vrec);
1453
1454		// if the matching critera are fulfilled, accept the rec-vertex that is closest in z
#	Line 1447 | Line 1526 | void PVStudy::fillMCHisto(std::vector<si
1526		}
1527		}
1528
1450	–
1529		void PVStudy::SetVarToZero() {
1530		//Greg's variables
1531		fntrk_ = 0;
#	Line 1569 | Line 1647 | void PVStudy::SetVarToZero() {
1647		pully_spl2_mct_[i] = 0;
1648		pullz_spl2_mct_[i] = 0;
1649		}
1572	–
1650		}
1651
1652		double PVStudy::sumPtSquared(const reco::Vertex & v)  {
#	Line 1582 | Line 1659 | double PVStudy::sumPtSquared(const reco:
1659		return sum;
1660		}
1661
1662	+
1663	+

Diff Legend

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