ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/Jeng/PVStudy/plugins/PVEffAnalyzer.cc

Comparing UserCode/Jeng/PVStudy/plugins/PVEffAnalyzer.cc (file contents):
Revision 1.1 by yygao, Mon Apr 19 16:07:17 2010 UTC vs.
Revision 1.8 by yygao, Fri Jun 4 17:46:09 2010 UTC

#	Line 1 | Line 1
1	+
2		// -*- C++ -*-
3		//
4		// Package:    PVEffAnalyzer
#	Line 40 | Line 41 | Implementation:
41		#include "FWCore/ServiceRegistry/interface/Service.h"
42		#include "CommonTools/UtilAlgos/interface/TFileService.h"
43		#include "UserCode/PVStudy/interface/PVEffAnalyzer.h"
44	<	//
44	>	#include "DataFormats/Common/interface/RefToBaseVector.h"
45	>	#include "DataFormats/Common/interface/RefToBase.h"
46	>	// Vertex
47		#include "DataFormats/VertexReco/interface/Vertex.h"
48		#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
49
#	Line 49 | Line 52 | Implementation:
52		#include <SimDataFormats/Vertex/interface/SimVertexContainer.h>
53		#include <SimDataFormats/Track/interface/SimTrack.h>
54		#include <SimDataFormats/Track/interface/SimTrackContainer.h>
55	+	// TrackingParticle Associators
56	+	#include "SimTracker/TrackAssociation/interface/TrackAssociatorByHits.h"
57	+	#include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
58	+	#include "SimTracker/Records/interface/TrackAssociatorRecord.h"
59		// BeamSpot
60		#include "DataFormats/BeamSpot/interface/BeamSpot.h"
61
62	+	// For the clusters
63	+	#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
64	+	#include "TrackingTools/Records/interface/TransientTrackRecord.h"
65	+	#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
66	+	#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
67	+
68		//root
69		#include <TROOT.h>
70		#include <TF1.h>
#	Line 60 | Line 73 | Implementation:
73		#include <TPaveStats.h>
74		#include <TPad.h>
75
76	+
77	+
78		using namespace std;
79	+	using namespace edm;
80	+	using namespace reco;
81	+
82		typedef math::XYZTLorentzVectorF LorentzVector;
83		typedef math::XYZPoint Point;
84
85		PVEffAnalyzer::PVEffAnalyzer(const edm::ParameterSet& iConfig)
86	+	:theTrackClusterizer(iConfig.getParameter<edm::ParameterSet>("TkClusParameters")),
87	+	theTrackFilter(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters"))
88		{
89		//=======================================================================
90		// Get configuration for TrackTupleMaker
#	Line 78 | Line 98 | PVEffAnalyzer::PVEffAnalyzer(const edm::
98		splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
99		verbose_                   = iConfig.getUntrackedParameter<bool>("verbose",false);
100		realData_                  = iConfig.getUntrackedParameter<bool>("realData",false);
101	+	useTP_                     = iConfig.getUntrackedParameter<bool>("useTP",false);
102	+	useAssociator_             = iConfig.getUntrackedParameter<bool>("useAssociator",false);
103		histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
104	+	outputfilename_            = iConfig.getUntrackedParameter<string>("OutputFileName");
105		nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
106		nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
107		zsigncut_                  = iConfig.getUntrackedParameter<double>("zsigncut");
108	+	ptcut_                     = iConfig.getUntrackedParameter<double>("ptcut");
109	+	analyze_                   = iConfig.getUntrackedParameter<bool>("analyze",false);
110	+	saventuple_                   = iConfig.getUntrackedParameter<bool>("saventuple",false);
111		bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
112	<
113	<
112	>	reqCluster_               = iConfig.getUntrackedParameter<bool>("reqCluster",false);
113	>
114		// Specify the data mode vector
115		if(realData_) datamode.push_back(0);
116		else {
117		datamode.push_back(0);
118		datamode.push_back(1);
119		}
120	+
121	+	if(saventuple_) {
122	+	file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
123	+	// pvtxtree_ analyzing the pvtxs ootb
124	+	pvtxtree_ = new TTree("pvtxtree","pvtxtree");
125	+	if(!realData_) {
126	+	pvtxtree_->Branch("nsimPV",&nsimPV_,"nsimPV/I");
127	+	pvtxtree_->Branch("nsimTrkPV",&nsimTrkPV_,"nsimTrkPV[nsimPV]/I");
128	+	pvtxtree_->Branch("simx",&simx_,"simx[nsimPV]/D");
129	+	pvtxtree_->Branch("simy",&simy_,"simy[nsimPV]/D");
130	+	pvtxtree_->Branch("simz",&simz_,"simz[nsimPV]/D");
131	+	pvtxtree_->Branch("simptsq",&simptsq_,"simptsq[nsimPV]/D");
132	+	pvtxtree_->Branch("simDeltaZ",&simDeltaZ_,"simDeltaZ[nsimPV]/D");
133	+	}
134	+	}
135
136		theFile = new TFile(histoFileName_.c_str(), "RECREATE");
137		theFile->mkdir("Summary");
138		theFile->cd();
139
140	+
141		// Book MC only plots
142		if (!realData_) {
143		h_gen = new PVEffHistograms();
#	Line 120 | Line 162 | PVEffAnalyzer::PVEffAnalyzer(const edm::
162
163		PVEffAnalyzer::~PVEffAnalyzer()
164		{
123	–
165		// do anything here that needs to be done at desctruction time
166		// (e.g. close files, deallocate resources etc.)
167		theFile->cd();
#	Line 134 | Line 175 | PVEffAnalyzer::~PVEffAnalyzer()
175		//
176		// member functions
177		//
178	<	std::vector<PVEffAnalyzer::simPrimaryVertex> PVEffAnalyzer::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
178	>	std::vector<PVEffAnalyzer::simPrimaryVertex> PVEffAnalyzer::getSimPVs(const Handle<HepMCProduct> & evtMC, std::string suffix, double & ptcut_)
179		{
180		std::vector<PVEffAnalyzer::simPrimaryVertex> simpv;
181		const HepMC::GenEvent* evt=evtMC->GetEvent();
#	Line 221 | Line 262 | std::vector<PVEffAnalyzer::simPrimaryVer
262		vp->ptot.setPz(vp->ptot.pz()+m.pz());
263		vp->ptot.setE(vp->ptot.e()+m.e());
264		vp->ptsq+=(m.perp())*(m.perp());
265	<	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
265	>	if ( (m.perp()>ptcut_) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
266		vp->nGenTrk++;
267		}
268		}
#	Line 232 | Line 273 | std::vector<PVEffAnalyzer::simPrimaryVer
273		return simpv;
274		}
275
276	<	std::vector<PVEffAnalyzer::simPrimaryVertex> PVEffAnalyzer::getSimPVs(const Handle<HepMCProduct> evtMC,
277	<	const Handle<SimVertexContainer> simVtxs,
278	<	const Handle<SimTrackContainer> simTrks)
279	<	{
280	<	// simvertices don't have enough information to decide,
281	<	// genVertices don't have the simulated coordinates  ( with VtxSmeared they might)
282	<	// go through simtracks to get the link between simulated and generated vertices
283	<	std::vector<PVEffAnalyzer::simPrimaryVertex> simpv;
284	<	int idx=0;
285	<	for(SimTrackContainer::const_iterator t=simTrks->begin();
286	<	t!=simTrks->end(); ++t){
287	<	if ( !(t->noVertex()) && !(t->type()==-99) ){
288	<	double ptsq=0;
289	<	bool primary=false;   // something coming directly from the primary vertex
290	<	bool resonance=false; // resonance
250	<	bool track=false;     // undecayed, charged particle
251	<	HepMC::GenParticle* gp=evtMC->GetEvent()->barcode_to_particle( (*t).genpartIndex() );
252	<	if (gp) {
253	<	HepMC::GenVertex * gv=gp->production_vertex();
254	<	if (gv) {
255	<	for ( HepMC::GenVertex::particle_iterator
256	<	daughter =  gv->particles_begin(HepMC::descendants);
257	<	daughter != gv->particles_end(HepMC::descendants);
258	<	++daughter ) {
259	<	if (isFinalstateParticle(*daughter)){
260	<	ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
261	<	}
262	<	}
263	<	//primary =  ( gv->position().t()==0 );
264	<	primary = true;
265	<	h_gen->Fill1d("genPart_cT", gv->position().t()); // mm
266	<	h_gen->Fill1d("genPart_T", gv->position().t()/299.792458); // ns
267	<
268	<	//resonance= ( gp->mother() && isResonance(gp->mother()));  // in CLHEP/HepMC days
269	<	// no more mother pointer in the improved HepMC GenParticle
270	<	resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
271	<	if (gp->status()==1){
272	<	//track=((pdt->particle(gp->pdg_id()))->charge() != 0);
273	<	track=not isCharged(gp);
274	<	}
275	<	}
276	>
277	>	// get sim pv from TrackingParticles/TrackingVertex
278	>	std::vector<PVEffAnalyzer::simPrimaryVertex> PVEffAnalyzer::getSimPVs( const edm::Handle<TrackingVertexCollection> & tVC )
279	>	{
280	>	std::vector<simPrimaryVertex> simpv;
281	>	//std::cout <<"number of vertices " << tVC->size() << std::endl;
282	>
283	>	if(verbose_){std::cout << "getSimPVs TrackingVertexCollection " << std::endl;}
284	>
285	>	for (TrackingVertexCollection::const_iterator v = tVC -> begin(); v != tVC -> end(); ++v) {
286	>
287	>	if(verbose_){
288	>	std::cout << (v->eventId()).event() << v -> position() << v->g4Vertices().size() <<" "  <<v->genVertices().size() <<  "   t=" <<v->position().t()*1.e12 <<"    ==0:" <<(v->position().t()>0) << std::endl;
289	>	for( TrackingVertex::g4v_iterator gv=v->g4Vertices_begin(); gv!=v->g4Vertices_end(); gv++){
290	>	std::cout << *gv << std::endl;
291		}
292	+	std::cout << "----------" << std::endl;
293	+	}
294	+
295	+	//    bool hasMotherVertex=false;
296	+	if ((unsigned int) v->eventId().event()<simpv.size()) continue;
297	+	//if (v->position().t()>0) { continue;} // skip secondary vertices (obsolete + doesn't work)
298	+	if (fabs(v->position().z())>1000) continue;  // skip funny junk vertices
299	+
300	+	// could be a new vertex, check  all primaries found so far to avoid multiple entries
301	+	const double mm=1.0; // for tracking vertices
302	+	simPrimaryVertex sv(v->position().x()*mm,v->position().y()*mm,v->position().z()*mm);
303	+
304	+	for (TrackingParticleRefVector::iterator iTrack = v->daughterTracks_begin(); iTrack != v->daughterTracks_end(); ++iTrack){
305	+	//cout <<((**iTrack).eventId()).event() << " ";  // an iterator of Refs, dereference twice. Cool eyh?
306	+	sv.eventId=(**iTrack).eventId();
307	+	sv.ptsq += pow((**iTrack).pt(),2);
308	+	sv.nGenTrk++;
309	+	}
310
311	<	const HepMC::FourVector & v=(*simVtxs)[t->vertIndex()].position();
312	<	//const HepLorentzVector & v=(*simVtxs)[t->vertIndex()].position();
313	<	if(primary or resonance){
314	<	{
315	<	// check all primaries found so far to avoid multiple entries
316	<	bool newVertex=true;
317	<	for(std::vector<PVEffAnalyzer::simPrimaryVertex>::iterator v0=simpv.begin();
318	<	v0!=simpv.end(); v0++){
319	<	if( (fabs(v0->x-v.x())<0.001) && (fabs(v0->y-v.y())<0.001) && (fabs(v0->z-v.z())<0.001) ){
320	<	if (track) {
321	<	v0->simTrackIndex.push_back(idx);
322	<	if (ptsq>(*v0).ptsq){(*v0).ptsq=ptsq;}
323	<	}
324	<	newVertex=false;
325	<	}
326	<	}
327	<	if(newVertex && !resonance){
328	<	PVEffAnalyzer::simPrimaryVertex anotherVertex(v.x(),v.y(),v.z());
329	<	if (track) anotherVertex.simTrackIndex.push_back(idx);
330	<	anotherVertex.ptsq=ptsq;
331	<	simpv.push_back(anotherVertex);
332	<	}
333	<	}//
311	>	simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
312	>	for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
313	>	v0!=simpv.end(); v0++){
314	>	if( (sv.eventId==v0->eventId) && (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
315	>	vp=&(*v0);
316	>	break;
317	>	}
318	>	}
319	>	if(!vp){
320	>	// this is a new vertex
321	>	if(verbose_){std::cout << "this is a new vertex " << sv.eventId.event() << "   "  << sv.x << " " << sv.y << " " << sv.z <<std::endl;}
322	>	simpv.push_back(sv);
323	>	vp=&simpv.back();
324	>	}else{
325	>	if(verbose_){std::cout << "this is not a new vertex"  << sv.x << " " << sv.y << " " << sv.z <<std::endl;}
326	>	}
327	>
328	>
329	>	// Loop over daughter track(s)
330	>	if(verbose_){
331	>	for (TrackingVertex::tp_iterator iTP = v -> daughterTracks_begin(); iTP != v -> daughterTracks_end(); ++iTP) {
332	>	std::cout << "  Daughter momentum:      " << (*(*iTP)).momentum();
333	>	std::cout << "  Daughter type     " << (*(*iTP)).pdgId();
334	>	std::cout << std::endl;
335		}
336	<
337	<	}// simtrack has vertex and valid type
338	<	idx++;
339	<	}//simTrack loop
336	>	}
337	>	}
338	>	if(verbose_){
339	>	cout << "------- PVEffAnalyzer simPVs from TrackingVertices -------" <<  endl;
340	>	for(std::vector<simPrimaryVertex>::iterator v0=simpv.begin();
341	>	v0!=simpv.end(); v0++){
342	>	cout << "z=" << v0->z << "  event=" << v0->eventId.event() << endl;
343	>	}
344	>	cout << "-----------------------------------------------" << endl;
345	>	}
346		return simpv;
347		}
348
349	+
350	+
351	+
352		// ------------ method called to for each event  ------------
353		void
354		PVEffAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
355		{
356		using namespace edm;
357		using namespace reco;
358	<
358	>
359	>	edm::ESHandle<TransientTrackBuilder> theB;
360	>	iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB);
361	>
362		//========================================================================
363		// Step 0: Prepare root variables and get information from the Event
364		//========================================================================
#	Line 320 | Line 366 | PVEffAnalyzer::analyze(const edm::Event&
366		edm::LogInfo("Debug")<<"[PVEffAnalyzer]"<<endl;
367
368		// ====== TrackCollection
369	<	static const reco::TrackCollection s_empty_trackColl;
324	<	const reco::TrackCollection *trackColl = &s_empty_trackColl;
325	<	edm::Handle<reco::TrackCollection> trackCollectionHandle;
369	>	edm::Handle<edm::View<reco::Track> > trackCollectionHandle;
370		iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle);
371	<	if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
328	<	trackColl = trackCollectionHandle.product();
329	<	} else {
330	<	edm::LogInfo("Debug") << "[PVEffAnalyzer] trackCollection cannot be found -> using empty collection of same type." <<endl;
331	<	}
371	>
372		// ====== splitTrackCollection1
373	<	static const reco::TrackCollection s_empty_splitTrackColl1;
334	<	const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
335	<	edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
373	>	edm::Handle<edm::View<reco::Track> > splitTrackCollection1Handle;
374		iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle);
375	<	if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
338	<	splitTrackColl1 = splitTrackCollection1Handle.product();
339	<	} else {
340	<	edm::LogInfo("Debug") << "[PVEffAnalyzer] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
341	<	}
375	>
376		// ====== splitTrackCollection2
377	<	static const reco::TrackCollection s_empty_splitTrackColl2;
344	<	const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
345	<	edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
377	>	edm::Handle<edm::View<reco::Track> > splitTrackCollection2Handle;
378		iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle);
347	–	if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
348	–	splitTrackColl2 = splitTrackCollection2Handle.product();
349	–	} else {
350	–	edm::LogInfo("Debug") << "[PVEffAnalyzer] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
351	–	}
379
380		// ======= PrimaryVertexCollection
381		static const reco::VertexCollection s_empty_vertexColl;
#	Line 410 | Line 437 | PVEffAnalyzer::analyze(const edm::Event&
437
438		// ======= Analyze MC efficiency
439		if (!realData_) {
440	<	bool MC=false;
440	>
441	>	bool MC=false; // the flag to see if the MC level information is complete
442	>
443		Handle<HepMCProduct> evtMC;
444		iEvent.getByLabel("generator",evtMC);
445	<	if (!evtMC.isValid()) {
446	<	MC=false;
447	<	edm::LogInfo("Debug") << "[PVEffAnalyzer] no HepMCProduct found"<< endl;
448	<	} else {
449	<	edm::LogInfo("Debug") << "[PVEffAnalyzer] generator HepMCProduct found"<< endl;
445	>
446	>	edm::Handle<TrackingParticleCollection>  TPCollectionH ;
447	>	bool gotTP=iEvent.getByLabel("mergedtruth","MergedTrackTruth",TPCollectionH);
448	>
449	>	edm::Handle<TrackingVertexCollection>    TVCollectionH ;
450	>	bool gotTV=iEvent.getByLabel("mergedtruth","MergedTrackTruth",TVCollectionH);
451	>
452	>	std::vector<simPrimaryVertex> simpv;
453	>
454	>	if(gotTV){
455		MC=true;
456	+	if(verbose_){
457	+	cout << "Found Tracking Vertices " << endl;
458	+	}
459	+	simpv=getSimPVs(TVCollectionH);
460		}
461	+	else if (evtMC.isValid()) {
462	+	MC=true;
463	+	edm::LogInfo("Debug") << "[PVEffAnalyzer] generator HepMCProduct found"<< endl;
464	+	} else {
465	+	edm::LogInfo("Debug") << "[PVEffAnalyzer] no HepMCProduct found"<< endl;
466	+	MC=false;
467	+	}
468	+
469	+
470		if(MC){
471		TString suffix = "_mct";
472		// make a list of primary vertices:
426	–	std::vector<simPrimaryVertex> simpv;
427	–	simpv=getSimPVs(evtMC,"");
428	–	//   simpv=getSimPVs(evtMC, simVtxs, simTrks);
473		h_gen->Fill1d("nsimPV", simpv.size());
474
475		int isimpv = 0;
476	<
477	<	for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
478	<	vsim!=simpv.end(); vsim++, isimpv++){
479	<	//nsimTrkPV_[isimpv]  =vsim->nGenTrk;
480	<	//simx_[isimpv] = vsim->x;
481	<	//simy_[isimpv] = vsim->y;
482	<	//simz_[isimpv] = vsim->y;
483	<	//simptsq_[isimpv] = vsim->ptsq;
484	<	if(verbose_ && simpv.size() > 1 )
485	<	std::cout<<"Simulated Vertex # " << isimpv << ": ptsq = " << vsim->ptsq<<std::endl;
476	>	nsimPV_ = simpv.size();
477	>
478	>	for(int isimpv = 0; isimpv<int(simpv.size());isimpv++) {
479	>	simPrimaryVertex vsim = simpv.at(isimpv);
480	>	nsimTrkPV_[isimpv]  =vsim.nGenTrk;
481	>	simx_[isimpv] = vsim.x;
482	>	simy_[isimpv] = vsim.y;
483	>	simz_[isimpv] = vsim.z;
484	>	simptsq_[isimpv] = vsim.ptsq;
485	>	if(simpv.size() > 1 )
486	>	if(verbose_)
487	>	std::cout<<"Simulated Vertex # " << isimpv << ": ptsq = " << vsim.ptsq<<std::endl;
488	>	if(isimpv >=1 )
489	>	simDeltaZ_[isimpv-1] = vsim.z-simz_[0];
490		}
491
492		// Just analyze the first vertex in the collection
493		simPrimaryVertex vsim = *simpv.begin();
494	<	h_summary->Fill1d(TString("denom_ntrack"+suffix), int(vsim.nGenTrk));
495	<	if ( isAssoVertex(vsim, *vertexColl->begin(), zsigncut_)  )
496	<	h_summary->Fill1d(TString("numer_ntrack"+suffix), int(vsim.nGenTrk));
497	<	}
498	<	} // End of Analyzing MC Efficiency
499	<
500	<	// ======= Analyze efficiency with split method
501	<	if ( !isGoodSplitEvent( *vertexColl->begin()) ) return;
502	<	if ( isGoodTagVertex( *splitVertexColl2->begin(), *vertexColl->begin(), zsigncut_) ) {
494	>	int ntracks(0);
495	>	if( useTP_ )
496	>	ntracks = vsim.nGenTrk;
497	>	else {
498	>	reco::RecoToSimCollection recSimColl;
499	>	if (useAssociator_) {
500	>	// get TrackingParticle
501	>	edm::Handle<TrackingParticleCollection>  TPCollectionHandle ;
502	>	iEvent.getByLabel("mergedtruth","MergedTrackTruth",TPCollectionHandle);
503	>	// get associators
504	>	edm::ESHandle<TrackAssociatorBase> theAssociator;
505	>	iSetup.get<TrackAssociatorRecord>().get("TrackAssociatorByHitsRecoDenom",theAssociator);
506	>	recSimColl = (theAssociator.product())->associateRecoToSim( trackCollectionHandle,
507	>	TPCollectionHandle,
508	>	& iEvent);
509	>	}
510	>
511	>	// == use the nubmer of tracks that are filtered and clusterized
512	>	// copy code from CMSSW/RecoVertex/PrimaryVertexProducer/src/PrimaryVertexProducerAlgorithm.cc
513	>	std::vector<reco::TransientTrack> seltks;
514	>	for(unsigned int i=0; i<trackCollectionHandle->size(); ++i){
515	>	edm::RefToBase<reco::Track> track(trackCollectionHandle, i);
516	>	if(useAssociator_) {
517	>	bool isSimMatched(false);
518	>	std::vector<std::pair<TrackingParticleRef, double> > tp;
519	>	if(recSimColl.find(track) != recSimColl.end()){
520	>	tp = recSimColl[track];
521	>	if (tp.size()!=0)
522	>	isSimMatched = true;
523	>	}
524	>	if(!isSimMatched) continue;
525	>	}
526	>	TransientTrack t_track = (*theB).build(*track);
527	>	t_track.setBeamSpot(bs);
528	>	if (theTrackFilter(t_track)) seltks.push_back(t_track);
529	>	}
530	>	std::vector< std::vector<reco::TransientTrack> > clusters =  theTrackClusterizer.clusterize(seltks);
531	>	if (clusters.size() == 1 && (clusters.begin()->size()) >1 ) {
532	>	for(unsigned idx_tt = 0; idx_tt < clusters.begin()->size(); idx_tt++) {
533	>	reco::Track track = (*(clusters.begin())).at(idx_tt).track();
534	>	if(track.pt()>ptcut_) // && TMath::Abs(track.eta())<2.5) // count only good tracks
535	>	ntracks++;
536	>	}
537	>	}
538	>	} // == end of getting the ntracks
539	>
540	>	// == Fill Efficiency Histogram
541	>	if(ntracks>1) {
542	>	h_summary->Fill1d(TString("eff_denom_ntrack"+suffix), ntracks);
543	>	if(!vertexColl->begin()->isFake())
544	>	h_summary->Fill1d(TString("deltazSign"+suffix), vsim.z - vertexColl->begin()->z());
545	>	if ( isAssoVertex(vsim, *vertexColl->begin(), zsigncut_)  )
546	>	h_summary->Fill1d(TString("eff_numer_ntrack"+suffix), ntracks);
547	>	}
548	>
549	>	// == Fill FakeRate Histogram
550	>	if(!vertexColl->begin()->isFake ()) {
551	>	h_summary->Fill1d(TString("fakerate_denom_ntrack"+suffix), int(vertexColl->begin()->tracksSize()));
552	>	if (!isAssoVertex(vsim, *vertexColl->begin(), zsigncut_)  )
553	>	h_summary->Fill1d(TString("fakerate_numer_ntrack"+suffix), int(vertexColl->begin()->tracksSize()));
554	>	}
555	>	}  // end of MC analyzing
556	>	} // end of MC access
557	>
558	>	// == Start of Split-Method
559	>	// Event selections
560	>	if ( !isGoodSplitEvent( *vertexColl->begin())) return;
561	>	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
562	>	t!=vertexColl->begin()->tracks_end(); t++) {
563	>	h_summary->Fill1d("trackweight",vertexColl->begin()->trackWeight(*t));
564	>	}
565	>
566	>	// == Start of FakeRate Analyzing
567	>
568	>	if(!splitVertexColl1->begin()->isFake()) {
569	>	h_summary->Fill1d("fakerate_denom_ntrack", int(splitVertexColl1->begin()->tracksSize()));
570	>	if ( !isGoodProbeVertex( *splitVertexColl1->begin(), *vertexColl->begin(), zsigncut_) )
571	>	h_summary->Fill1d("fakerate_numer_ntrack", int(splitVertexColl1->begin()->tracksSize()));
572	>	}
573	>
574	>	// == End of FakeRate Analyzing
575	>	// == Start of Efficiency Analyzing
576	>	std::vector< std::vector<reco::TransientTrack> > clusters1;
577	>	std::vector< std::vector<reco::TransientTrack> > clusters2;
578	>	if(splitTrackCollection1Handle.isValid())
579	>	clusters1 =  theTrackClusterizer.clusterize((*theB).build(splitTrackCollection1Handle));
580	>	if(splitTrackCollection2Handle.isValid())
581	>	clusters2 =  theTrackClusterizer.clusterize((*theB).build(splitTrackCollection2Handle));
582	>
583	>	// clustering requirement
584	>	if( !goodCluster (clusters1) || !goodCluster (clusters2) ) return;
585	>
586	>	// TagVertex Requirement
587	>	if ( isGoodTagVertex( *(splitVertexColl2->begin()), *(vertexColl->begin()), zsigncut_) ) {
588		if(verbose_) {
589	<	std::cout<<"splitTrackColl1->size() = " << int(splitTrackColl1->size()) << std::endl;
590	<	std::cout<<"splitTrackColl2->size() = " << int(splitTrackColl2->size()) << std::endl;
589	>	std::cout<<"splitTrackCollection1Handle->size() = " << int(splitTrackCollection1Handle->size()) << std::endl;
590	>	std::cout<<"splitTrackCollection2Handle->size() = " << int(splitTrackCollection2Handle->size()) << std::endl;
591	>	}
592	>	int nGoodTracks = 0;
593	>	for(unsigned int i = 0; i < splitTrackCollection1Handle->size(); i++ ) {
594	>	edm::RefToBase<reco::Track> track(splitTrackCollection1Handle, i);
595	>	if(track->pt() > ptcut_) // && abs(track->eta()) < 2.5 )
596	>	nGoodTracks++;
597		}
598	<	h_summary->Fill1d("denom_ntrack", int(splitTrackColl1->size()));
599	<	if ( isGoodProbeVertex( *splitVertexColl1->begin(), *vertexColl->begin(), zsigncut_) )
600	<	h_summary->Fill1d("numer_ntrack", int(splitTrackColl1->size()) );
598	>
599	>	h_summary->Fill1d("eff_denom_ntrack", nGoodTracks);
600	>	// ProbeVertex Requirement
601	>	if ( isGoodProbeVertex( *splitVertexColl1->begin(), *vertexColl->begin(), zsigncut_) ) {
602	>	h_summary->Fill1d("eff_numer_ntrack", nGoodTracks);
603	>	h_summary->Fill1d("avgWeight_orgvtx_eff",avgWeight(*vertexColl->begin()) );
604	>	h_summary->Fill1d("avgWeight_tagvtx_eff",avgWeight(*splitVertexColl2->begin()) );
605	>	h_summary->Fill1d("avgWeight_probevtx_eff",avgWeight(*splitVertexColl1->begin()) );
606	>	}
607	>	else
608	>	{
609	>	if (nGoodTracks >= 2) {
610	>	h_summary->Fill1d("avgWeight_orgvtx_ineff",avgWeight(*vertexColl->begin()) );
611	>	h_summary->Fill1d("avgWeight_tagvtx_ineff",avgWeight(*splitVertexColl2->begin()) );
612	>	h_summary->Fill1d("avgWeight_probevtx_ineff",avgWeight(*splitVertexColl1->begin()) );
613	>
614	>	std::cout<<"**********Inefficiency Found************** " <<std::endl;
615	>	std::cout<<"==== Original Vertex: " <<std::endl;
616	>	//printRecVtx(*vertexColl->begin());
617	>	printRecVtxs(vertexCollectionHandle);
618	>	std::cout<<"==== Probe Vertex: " <<std::endl;
619	>	//printRecVtx(*splitVertexColl1->begin());
620	>	printRecVtxs(splitVertexCollection1Handle);
621	>	std::cout<<"==== Tag Vertex: " <<std::endl;
622	>	//printRecVtx(*splitVertexColl2->begin());
623	>	printRecVtxs(splitVertexCollection2Handle);
624	>	std::cout<<"==== Probe cluster information "<< std::endl;
625	>	printCluster(clusters1, beamSpot);
626	>	std::cout<<"==== Tag cluster information "<< std::endl;
627	>	printCluster(clusters2, beamSpot);
628	>	std::cout<<"**********Inefficiency Found************** " <<std::endl;
629	>	}
630	>	}
631		}
632	+	// == End of Efficiency Analyzing
633	+	// == Save Ntuple
634	+	if(saventuple_) {
635	+	pvtxtree_->Fill();
636	+	}
637	+
638		}
639
640
#	Line 473 | Line 648 | PVEffAnalyzer::beginJob()
648		void
649		PVEffAnalyzer::endJob() {
650		edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl;
651	<
651	>
652	>	if(saventuple_) {
653	>	file_->cd();
654	>	pvtxtree_->Write();
655	>	file_->Close();
656	>	}
657	>
658	>
659		for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
660		string suffix;
661		edm::LogInfo("Debug")<<"datamode = "<< *it<<endl;
#	Line 483 | Line 665 | PVEffAnalyzer::endJob() {
665		case 1: suffix = "_mct";
666		break;
667		}
668	<	MakeEff(h_summary->ReadHisto1D(TString("numer_ntrack"+suffix)), h_summary->ReadHisto1D(TString("denom_ntrack"+suffix)), h_summary->ReadHisto1D(TString("eff_ntrack"+suffix)), false, 1);
668	>	if (analyze_) {
669	>	MakeEff(h_summary->ReadHisto1D(TString("eff_numer_ntrack"+suffix)), h_summary->ReadHisto1D(TString("eff_denom_ntrack"+suffix)), h_summary->ReadHisto1D(TString("eff_ntrack"+suffix)), false, 1);
670	>	MakeEff(h_summary->ReadHisto1D(TString("fakerate_numer_ntrack"+suffix)), h_summary->ReadHisto1D(TString("fakerate_denom_ntrack"+suffix)), h_summary->ReadHisto1D(TString("fakerate_ntrack"+suffix)), false, 1);
671	>	}
672		}
673		}
674
#	Line 579 | Line 764 | void PVEffAnalyzer::printSimTrks(const H
764		// vtx.ndof = 2*Sum(trackWeight) - 3
765
766		bool PVEffAnalyzer::isGoodSplitEvent( const reco::Vertex & org_vtx)
582	–	//, const reco::Vertex & tag_vtx, const reco::Vertex & probe_vtx)
767		{
768	<	return (  org_vtx.tracksSize() >= 6 && org_vtx.ndof() >= 6 );
769	<	//      tag_vtx.tracksSize() > = 2 && probe_vtx.tracksSize() > = 2      ) ;
768	>	if (  org_vtx.tracksSize() >= 6 && org_vtx.ndof() >= 6 && avgWeight(org_vtx) > 0.75  ) return true;
769	>	//if (  org_vtx.tracksSize() >= 6 && org_vtx.ndof() >= 6 ) return true;
770	>	else
771	>	return false;
772		}
773
774		// Comparing the tag vertex with the unsplit vertex position in z
#	Line 590 | Line 776 | bool PVEffAnalyzer::isGoodTagVertex( con
776		{
777		if(tag_vtx.isValid () && !tag_vtx.isFake()) {
778		float zsign = TMath::Abs(tag_vtx.z() - org_vtx.z() ) / TMath::Max(tag_vtx.zError(), org_vtx.zError() );
779	+	//float zsign = TMath::Abs(tag_vtx.z() - org_vtx.z() ) / sqrt(pow(tag_vtx.zError(),2)+pow(org_vtx.zError(),2) );
780		return (zsign < zsign_cut) ;
781		}
782		else
#	Line 601 | Line 788 | bool PVEffAnalyzer::isGoodProbeVertex( c
788		{
789		if(probe_vtx.isValid () && !probe_vtx.isFake()) {
790		float zsign = TMath::Abs(probe_vtx.z() - org_vtx.z() ) / TMath::Max(probe_vtx.zError(), org_vtx.zError() );
791	+	//float zsign = TMath::Abs(probe_vtx.z() - org_vtx.z() ) / sqrt(pow(probe_vtx.zError(),2)+pow(org_vtx.zError(),2) );
792		return (zsign < zsign_cut) ;
793		}
794		else
#	Line 618 | Line 806 | bool PVEffAnalyzer::isAssoVertex(const P
806		}
807
808
809	+
810		void PVEffAnalyzer::MakeEff(TH1D* numer, TH1D* denom, TH1D* eff,  //TGraphAsymmErrors* & gr_eff,
811		const bool rebin, const Float_t n ) {
812		eff->Divide( numer, denom, 1,1,"B" );
#	Line 648 | Line 837 | void PVEffAnalyzer::MakeEff(TH1D* numer,
837
838		for (int i = 1; i<eff->GetNbinsX()+1; i++) {
839		if(numer->GetBinContent(i) == 0 || denom->GetBinContent(i) == 0 ) continue;
840	<	float error = eff->GetBinContent(i)*sqrt(pow(numer->GetBinError(i)/numer->GetBinContent(i),2)
652	<	+pow(denom->GetBinError(i)/denom->GetBinContent(i),2));
840	>	float error = sqrt(eff->GetBinContent(i)*(1-eff->GetBinContent(i))/denom->GetBinContent(i)); // Binominal-Error
841		eff->SetBinError(i, error);
842		}
843		}
844
845	<	/*
846	<	// Copied from Validation/RecoTrack
847	<	void PVEffAnalyzer::setUpVectors(bool isRealData, int & nTrkMin_, int & nTrkMax_ ) {
848	<	std::vector<int>   totTAGv_ntracks, totPROBEv_ntracks;
849	<	if(!isRealData)
850	<	std::vector<int>   totSIMv_ntracks,totASSv_ntracks;
851	<	for (int k=nTrkMin_;k<nTrkMax_+1;k++) {
852	<	totTAGv_ntracks.push_back(0);
853	<	totPROBEv_ntracks.push_back(0);
854	<	if(!realData) {
855	<	totSIMv_ntracks.push_back(0);
856	<	totASSv_ntracks.push_back(0);
845	>	bool PVEffAnalyzer::goodCluster(std::vector< std::vector<reco::TransientTrack> > & clusters)
846	>	{
847	>	if( clusters.size() == 1 && ( clusters.begin()->size() > 1 )) return true;
848	>	else
849	>	return false;
850	>	}
851	>
852	>	void PVEffAnalyzer::printCluster(std::vector< std::vector<reco::TransientTrack> > & clusters, const Point & beamSpot)
853	>	{
854	>	//std::cout<<"theTrackClusterizer.zSeparation()"<<theTrackClusterizer.zSeparation()<<std::endl;
855	>	cout << "#clusters " << std::setw(3) << clusters.size() << endl;
856	>	int idx_cluster = 0;
857	>	for (vector< vector<reco::TransientTrack> >::const_iterator iclus = clusters.begin(); iclus != clusters.end(); iclus++, idx_cluster++) {
858	>	if((*iclus).size() == 0) continue;
859	>	cout << "=== Cluster # " << std::setw(3) << idx_cluster << std::setw(3) << " nTracks" << std::setw(3) << (*iclus).size() <<endl;
860	>	for (int idx_track = 0; idx_track < int((*iclus).size()) ; idx_track++) {
861	>	reco::Track t = (*iclus).at(idx_track).track();
862	>	cout << "== Trk #  " << std::setw(4) << idx_track
863	>	<< " nHit     " << std::setw(4) << t.hitPattern().numberOfValidHits() <<endl
864	>	<< " pt       " << std::setw(9) << std::fixed << std::setprecision(4) << t.pt()
865	>	<< " dpt      " << std::setw(9) << t.ptError() << endl
866	>	<< " dxy(BS)  " << std::setw(9) << t.dxy(beamSpot)
867	>	<< " sigmaDxy " << std::setw(9) << t.dxyError() << endl
868	>	<< " dz(BS)   " << std::setw(9) << t.dz(beamSpot)
869	>	<< " sigmaDz  " << std::setw(9) << t.dzError() << endl
870	>	<< " xPCA     " << std::setw(9) << t.vertex().x()
871	>	<< " yPCA     " << std::setw(9) << t.vertex().y()
872	>	<< " zPCA     " << std::setw(9) << t.vertex().z()
873	>	<< endl;
874		}
875		}
876
877	<	totTAG_ntracks.push_back(totTAGv_ntracks);
878	<	totPROBE_ntracks.push_back(totPROBEv_ntracks);
879	<	if(!realData) {
880	<	totSIM_ntracks.push_back(totSIMv_ntracks);
881	<	totASS_ntracks.push_back(totASSv_ntracks);
882	<	}
877	>	}
878	>	double PVEffAnalyzer::avgWeight(const reco::Vertex & vtx)
879	>	{
880	>	if(!vtx.isValid() || vtx.isFake() )
881	>	return 0;
882	>	else
883	>	return (vtx.ndof()+3.0)/2.0/double(vtx.tracksSize());
884		}
885
886	<	*/
886	>
887	>	void PVEffAnalyzer::SetVarToZero() {
888	>	nsimPV_ = 0;
889	>	//simpv
890	>	for(int i = 0; i < nMaxPVs_; i++) {
891	>	nsimTrkPV_[i] = 0;
892	>	simx_[i] = 0;
893	>	simy_[i] = 0;
894	>	simz_[i] = 0;
895	>	simptsq_[i] = 0;
896	>	if(i>=1)
897	>	simDeltaZ_[i-1] = 0;
898	>	}
899	>	}

Diff Legend

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