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.17 by yygao, Tue Jan 26 12:06:49 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	+	// 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		{
#	Line 72 | Line 90 | PVStudy::PVStudy(const edm::ParameterSet
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");
76	<	pixelVertexCollectionTag_  = iConfig.getParameter<edm::InputTag>("pixelVertexCollectionTag");
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");
83	–	zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
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	<
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 97 | Line 116 | PVStudy::PVStudy(const edm::ParameterSet
116		datamode.push_back(2);
117		datamode.push_back(3);
118		}
119	<
101	<	// 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 217 | 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		}
238	<
238	>
239		setRootStyle();
240
241	+	//========================================
242	+	// Booking histograms
243	+	//========================================
244	+
245		// Create a root file for histograms
246		theFile = new TFile(histoFileName_.c_str(), "RECREATE");
247		// make diretories
#	Line 227 | Line 249 | PVStudy::PVStudy(const edm::ParameterSet
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();
#	Line 240 | Line 266 | PVStudy::PVStudy(const edm::ParameterSet
266		else {
267		stringstream ss;
268		ss << i;
269	<	h_pvtrk->Init("pvTrk",ss.str(),"spl");
269	>	h_pvtrk->Init("pvTrk", ss.str(),"spl");
270		}
271		}
272		h_pixvtx->Init("pixVtx");
#	Line 280 | Line 306 | PVStudy::PVStudy(const edm::ParameterSet
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
311		}
312
313		PVStudy::~PVStudy()
314		{
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();
292	–	h_gen->Save();
322		h_misc->Save();
323		h_summary->Save();
324	<	if (analyze_)
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();
#	Line 335 | 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		edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
#	Line 517 | Line 551 | PVStudy::analyze(const edm::Event& iEven
551		// Initialize Root-tuple variables if needed
552		//=======================================================
553		//if(saventuple_)
554	<	SetVarToZero();
554	>	SetVarToZero();
555
556		//=======================================================
557		// Track accessors
#	Line 532 | Line 566 | PVStudy::analyze(const edm::Event& iEven
566		} else {
567		edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
568		}
535	–
569		//splitTrackCollection1
570		static const reco::TrackCollection s_empty_splitTrackColl1;
571		const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
#	Line 543 | Line 576 | PVStudy::analyze(const edm::Event& iEven
576		} else {
577		edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
578		}
546	–
579		//splitTrackCollection2
580		static const reco::TrackCollection s_empty_splitTrackColl2;
581		const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
#	Line 552 | Line 584 | PVStudy::analyze(const edm::Event& iEven
584		if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
585		splitTrackColl2 = splitTrackCollection2Handle.product();
586		} else {
587	<	edm::LogInfo("Debug") << "[PVStudy] 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	<
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	<
589	>
590		//=======================================================
591		// PVTX accessors
592		//=======================================================
#	Line 578 | Line 598 | PVStudy::analyze(const edm::Event& iEven
598		if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
599		vertexColl = vertexCollectionHandle.product();
600		} else {
601	<	edm::LogInfo("Debug") << "[PVStudy] 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		}
583	–
603		//splitVertexCollection1
604		static const reco::VertexCollection s_empty_splitVertexColl1;
605		const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
#	Line 590 | Line 609 | PVStudy::analyze(const edm::Event& iEven
609		splitVertexColl1 = splitVertexCollection1Handle.product();
610		} else {
611		edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
612	<	}
594	<
612	>	}
613		//splitVertexCollection2
614		static const reco::VertexCollection s_empty_splitVertexColl2;
615		const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
#	Line 602 | Line 620 | PVStudy::analyze(const edm::Event& iEven
620		} else {
621		edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
622		}
623	+
624
625	<	edm::LogInfo("Debug")<<"[PVStudy] 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 624 | Line 665 | PVStudy::analyze(const edm::Event& iEven
665		}catch(const Exception&){
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	<	edm::LogInfo("Debug") << "[PVStudy] 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	<	//fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple) {
746	<	fillTrackHistoInPV(pixelVertexColl, 4, false, true);
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();
#	Line 696 | Line 795 | PVStudy::analyze(const edm::Event& iEven
795		h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
796		h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
797
798	<	//=================================================================
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);
709	<
710	<	if(splitVertexColl1->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
711	<	fillTrackHistoInPV(splitVertexColl2, 2, false, true);
808	>	fillTrackHistoInPV(splitVertexColl1, 1, false, true, beamSpot);
809
810	+	if(splitVertexColl2->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
811	+	fillTrackHistoInPV(splitVertexColl2, 2, false, true, beamSpot);
812
813		//=======================================================
814		// Compare offlinePrimaryVertices with pixelVertices
#	Line 725 | Line 824 | PVStudy::analyze(const edm::Event& iEven
824		h_pixvtx->Fill1d("reczPV_minus_reczPxlPV", recz_[0] - recz_pxlpvtx_[0]);
825		}
826
827	<
828	<
730	<	//==========================================================
827	>
828	>	//==========================================================
829		// Fill secondary/primary min separations for multi-vertices
830		//==========================================================
831
#	Line 794 | Line 892 | PVStudy::analyze(const edm::Event& iEven
892		} // End of  if(vertexColl->size()>1) {
893
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 802 | 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	<
905	>
906		edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
907		reco::VertexCollection s_empty_matchedVertexColl1;
908		reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
#	Line 825 | Line 925 | PVStudy::analyze(const edm::Event& iEven
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	<
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		edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
933
#	Line 862 | Line 964 | PVStudy::analyze(const edm::Event& iEven
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 912 | 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		edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
1063
1064		edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
1065
1066	<	//=======================================================
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_) {
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 1010 | Line 1120 | PVStudy::endJob() {
1120		pvtxtree_->Write();
1121		file_->Close();
1122		}
1123	<	}
1014	<
1123	>
1124
1125	+	}
1126
1127		// Match a recovertex with a simvertex
1128		// ? Should the cut be parameter dependent?
#	Line 1020 | 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 1127 | Line 1239 | void PVStudy::fillHisto(res r, error er,
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 1189 | Line 1302 | PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo
1302		ntk);
1303		}
1304
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		edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
#	Line 1219 | Line 1334 | void PVStudy::fit(TH1 *hdist, double fit
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);
#	Line 1235 | Line 1350 | void PVStudy::fit(TH1 *hdist, double fit
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
#	Line 1247 | Line 1363 | void PVStudy::fillTrackHisto(const reco:
1363		++itTrack) {
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 1261 | 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 1316 | Line 1434 | void PVStudy::fillTrackHistoInPV(const r
1434		// For histogram only fill the leading pvtx parameters
1435		if(fillHisto) {
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++) {
#	Line 1325 | Line 1446 | void PVStudy::fillTrackHistoInPV(const r
1446		}
1447		else {
1448		h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
1449	<	h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
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		}
#	Line 1336 | Line 1461 | void PVStudy::fillTrackHistoInPV(const r
1461		// exception thrown when trying to use linked track
1462		//          h_pvtrk->Fill1d("trkPtPV", 0.);
1463		}
1464	+	h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
1465		}
1466	+
1467
1468		}
1469
#	Line 1369 | Line 1496 | void PVStudy::fillMCHisto(std::vector<si
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) && vrec->isValid() && !vrec->isFake() ) {
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 1447 | Line 1574 | void PVStudy::fillMCHisto(std::vector<si
1574		}
1575		}
1576
1450	–
1577		void PVStudy::SetVarToZero() {
1578		//Greg's variables
1579		fntrk_ = 0;
#	Line 1569 | Line 1695 | void PVStudy::SetVarToZero() {
1695		pully_spl2_mct_[i] = 0;
1696		pullz_spl2_mct_[i] = 0;
1697		}
1572	–
1698		}
1699
1700		double PVStudy::sumPtSquared(const reco::Vertex & v)  {
#	Line 1582 | Line 1707 | double PVStudy::sumPtSquared(const reco:
1707		return sum;
1708		}
1709
1710	+
1711	+

Diff Legend

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