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.4 by yygao, Tue Apr 20 12:39:48 2010 UTC vs.
Revision 1.5 by yygao, Fri Apr 30 12:19:52 2010 UTC

#	Line 41 | 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 50 | 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	<
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"
#	Line 67 | 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;
#	Line 75 | Line 83 | typedef math::XYZTLorentzVectorF Lorentz
83		typedef math::XYZPoint Point;
84
85		PVEffAnalyzer::PVEffAnalyzer(const edm::ParameterSet& iConfig)
86	<	:theTrackClusterizer(iConfig.getParameter<edm::ParameterSet>("TkClusParameters"))
86	>	:theTrackClusterizer(iConfig.getParameter<edm::ParameterSet>("TkClusParameters")),
87	>	theTrackFilter(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters"))
88		{
89		//=======================================================================
90		// Get configuration for TrackTupleMaker
#	Line 89 | 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		nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
105		nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
#	Line 250 | Line 261 | PVEffAnalyzer::analyze(const edm::Event&
261		{
262		using namespace edm;
263		using namespace reco;
264	<
264	>
265	>	edm::ESHandle<TransientTrackBuilder> theB;
266	>	iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB);
267	>
268		//========================================================================
269		// Step 0: Prepare root variables and get information from the Event
270		//========================================================================
#	Line 258 | Line 272 | PVEffAnalyzer::analyze(const edm::Event&
272		edm::LogInfo("Debug")<<"[PVEffAnalyzer]"<<endl;
273
274		// ====== TrackCollection
275	<	static const reco::TrackCollection s_empty_trackColl;
262	<	const reco::TrackCollection *trackColl = &s_empty_trackColl;
263	<	edm::Handle<reco::TrackCollection> trackCollectionHandle;
275	>	edm::Handle<edm::View<reco::Track> > trackCollectionHandle;
276		iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle);
277	<	if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
266	<	trackColl = trackCollectionHandle.product();
267	<	} else {
268	<	edm::LogInfo("Debug") << "[PVEffAnalyzer] trackCollection cannot be found -> using empty collection of same type." <<endl;
269	<	}
277	>
278		// ====== splitTrackCollection1
279	<	static const reco::TrackCollection s_empty_splitTrackColl1;
272	<	const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
273	<	edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
279	>	edm::Handle<edm::View<reco::Track> > splitTrackCollection1Handle;
280		iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle);
281	<	if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
276	<	splitTrackColl1 = splitTrackCollection1Handle.product();
277	<	} else {
278	<	edm::LogInfo("Debug") << "[PVEffAnalyzer] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
279	<	}
281	>
282		// ====== splitTrackCollection2
283	<	static const reco::TrackCollection s_empty_splitTrackColl2;
282	<	const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
283	<	edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
283	>	edm::Handle<edm::View<reco::Track> > splitTrackCollection2Handle;
284		iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle);
285	–	if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
286	–	splitTrackColl2 = splitTrackCollection2Handle.product();
287	–	} else {
288	–	edm::LogInfo("Debug") << "[PVEffAnalyzer] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
289	–	}
285
286		// ======= PrimaryVertexCollection
287		static const reco::VertexCollection s_empty_vertexColl;
#	Line 380 | Line 375 | PVEffAnalyzer::analyze(const edm::Event&
375
376		// Just analyze the first vertex in the collection
377		simPrimaryVertex vsim = *simpv.begin();
378	<	h_summary->Fill1d(TString("denom_ntrack"+suffix), int(vsim.nGenTrk));
379	<	if ( isAssoVertex(vsim, *vertexColl->begin(), zsigncut_)  )
380	<	h_summary->Fill1d(TString("numer_ntrack"+suffix), int(vsim.nGenTrk));
378	>	int ntracks(0);
379	>	if( useTP_ )
380	>	ntracks = vsim.nGenTrk;
381	>	else {
382	>	reco::RecoToSimCollection recSimColl;
383	>	if (useAssociator_) {
384	>	// get TrackingParticle
385	>	edm::Handle<TrackingParticleCollection>  TPCollectionHandle ;
386	>	iEvent.getByLabel("mergedtruth","MergedTrackTruth",TPCollectionHandle);
387	>	// get associators
388	>	edm::ESHandle<TrackAssociatorBase> theAssociator;
389	>	iSetup.get<TrackAssociatorRecord>().get("TrackAssociatorByHitsRecoDenom",theAssociator);
390	>	recSimColl = (theAssociator.product())->associateRecoToSim( trackCollectionHandle,
391	>	TPCollectionHandle,
392	>	& iEvent);
393	>	}
394	>
395	>	// == use the nubmer of tracks that are filtered and clusterized
396	>	// copy code from CMSSW/RecoVertex/PrimaryVertexProducer/src/PrimaryVertexProducerAlgorithm.cc
397	>	std::vector<reco::TransientTrack> seltks;
398	>	for(unsigned int i=0; i<trackCollectionHandle->size(); ++i){
399	>	edm::RefToBase<reco::Track> track(trackCollectionHandle, i);
400	>	if(useAssociator_) {
401	>	bool isSimMatched(false);
402	>	std::vector<std::pair<TrackingParticleRef, double> > tp;
403	>	if(recSimColl.find(track) != recSimColl.end()){
404	>	tp = recSimColl[track];
405	>	if (tp.size()!=0)
406	>	isSimMatched = true;
407	>	}
408	>	if(!isSimMatched) continue;
409	>	}
410	>	TransientTrack t_track = (*theB).build(*track);
411	>	t_track.setBeamSpot(bs);
412	>	if (theTrackFilter(t_track)) seltks.push_back(t_track);
413	>	}
414	>	std::vector< std::vector<reco::TransientTrack> > clusters =  theTrackClusterizer.clusterize(seltks);
415	>	if (clusters.size() == 1 && (clusters.begin()->size()) >1 ) {
416	>	for(unsigned idx_tt = 0; idx_tt < clusters.begin()->size(); idx_tt++) {
417	>	reco::Track track = (*(clusters.begin())).at(idx_tt).track();
418	>	if(track.pt()>ptcut_) // && TMath::Abs(track.eta())<2.5) // count only good tracks
419	>	ntracks++;
420	>	}
421	>	}
422	>	}
423	>
424	>	// == end of getting the ntracks
425	>	if(ntracks>1) {
426	>	h_summary->Fill1d(TString("denom_ntrack"+suffix), ntracks);
427	>	if(!vertexColl->begin()->isFake())
428	>	h_summary->Fill1d(TString("deltazSign"+suffix), vsim.z - vertexColl->begin()->z());
429	>	if ( isAssoVertex(vsim, *vertexColl->begin(), zsigncut_)  )
430	>	h_summary->Fill1d(TString("numer_ntrack"+suffix), ntracks);
431	>	}
432		}
433		} // End of Analyzing MC Efficiency
434
435	<	// Get the Builder for TrackClusters
436	<	edm::ESHandle<TransientTrackBuilder> theB;
437	<	iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB);
438	<
439	<	std::vector< std::vector<reco::TransientTrack> > clusters1 =  theTrackClusterizer.clusterize((*theB).build(splitTrackCollection1Handle));
440	<	std::vector< std::vector<reco::TransientTrack> > clusters2 =  theTrackClusterizer.clusterize((*theB).build(splitTrackCollection2Handle));
435	>	// == Start of Split-Method
436	>	// Event selections
437	>	if ( !isGoodSplitEvent( *vertexColl->begin())) return;
438	>	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
439	>	t!=vertexColl->begin()->tracks_end(); t++) {
440	>	h_summary->Fill1d("trackweight",vertexColl->begin()->trackWeight(*t));
441	>	}
442	>	std::vector< std::vector<reco::TransientTrack> > clusters1;
443	>	std::vector< std::vector<reco::TransientTrack> > clusters2;
444	>	if(splitTrackCollection1Handle.isValid())
445	>	clusters1 =  theTrackClusterizer.clusterize((*theB).build(splitTrackCollection1Handle));
446	>	if(splitTrackCollection2Handle.isValid())
447	>	clusters2 =  theTrackClusterizer.clusterize((*theB).build(splitTrackCollection2Handle));
448
449	<	bool withValidCluster_tag = false;
450	<	bool withValidCluster_probe = false;
449	>	// clustering requirement
450	>	if( !goodCluster (clusters1) || !goodCluster (clusters2) ) return;
451
452	<
453	<	for (vector< vector<reco::TransientTrack> >::const_iterator iclus = clusters1.begin(); iclus != clusters1.end(); iclus++) {
401	<	if((*iclus).size()>1)  withValidCluster_probe = true;
402	<	}
403	<	for (vector< vector<reco::TransientTrack> >::const_iterator iclus = clusters2.begin(); iclus != clusters2.end(); iclus++) {
404	<	if((*iclus).size()>1)  withValidCluster_tag = true;
405	<	}
406	<
407	<	// ======= Analyze efficiency with split method
408	<	if( reqCluster_ ) {
409	<	if ( !isGoodSplitEvent( *vertexColl->begin()) || !withValidCluster_tag || !withValidCluster_probe ) return;
410	<	}
411	<	else {
412	<	if ( !isGoodSplitEvent( *vertexColl->begin())) return;
413	<	}
414	<	if ( isGoodTagVertex( *splitVertexColl2->begin(), *vertexColl->begin(), zsigncut_) ) {
452	>	// TagVertex Requirement
453	>	if ( isGoodTagVertex( *(splitVertexColl2->begin()), *(vertexColl->begin()), zsigncut_) ) {
454		if(verbose_) {
455	<	std::cout<<"splitTrackColl1->size() = " << int(splitTrackColl1->size()) << std::endl;
456	<	std::cout<<"splitTrackColl2->size() = " << int(splitTrackColl2->size()) << std::endl;
455	>	std::cout<<"splitTrackCollection1Handle->size() = " << int(splitTrackCollection1Handle->size()) << std::endl;
456	>	std::cout<<"splitTrackCollection2Handle->size() = " << int(splitTrackCollection2Handle->size()) << std::endl;
457		}
419	–	// the same cut as in the sim->nGenTrk
458		int nGoodTracks = 0;
459	<	for(reco::TrackCollection::const_iterator track = splitTrackColl1->begin();
460	<	track!= splitTrackColl1->end(); ++track) {
461	<	if(track->pt() > ptcut_ && abs(track->eta()) < 2.5 )
459	>	for(unsigned int i = 0; i < splitTrackCollection1Handle->size(); i++ ) {
460	>	edm::RefToBase<reco::Track> track(splitTrackCollection1Handle, i);
461	>	if(track->pt() > ptcut_) // && abs(track->eta()) < 2.5 )
462		nGoodTracks++;
463		}
464
465		h_summary->Fill1d("denom_ntrack", nGoodTracks);
466	<	if ( isGoodProbeVertex( *splitVertexColl1->begin(), *vertexColl->begin(), zsigncut_) )
466	>	// ProbeVertex Requirement
467	>	if ( isGoodProbeVertex( *splitVertexColl1->begin(), *vertexColl->begin(), zsigncut_) ) {
468		h_summary->Fill1d("numer_ntrack", nGoodTracks);
469	+	h_summary->Fill1d("avgWeight_orgvtx_eff",avgWeight(*vertexColl->begin()) );
470	+	h_summary->Fill1d("avgWeight_tagvtx_eff",avgWeight(*splitVertexColl2->begin()) );
471	+	h_summary->Fill1d("avgWeight_probevtx_eff",avgWeight(*splitVertexColl1->begin()) );
472	+	}
473	+	else
474	+	{
475	+	if (nGoodTracks >= 2) {
476	+	h_summary->Fill1d("avgWeight_orgvtx_ineff",avgWeight(*vertexColl->begin()) );
477	+	h_summary->Fill1d("avgWeight_tagvtx_ineff",avgWeight(*splitVertexColl2->begin()) );
478	+	h_summary->Fill1d("avgWeight_probevtx_ineff",avgWeight(*splitVertexColl1->begin()) );
479	+
480	+	std::cout<<"**********Inefficiency Found************** " <<std::endl;
481	+	std::cout<<"==== Original Vertex: " <<std::endl;
482	+	//printRecVtx(*vertexColl->begin());
483	+	printRecVtxs(vertexCollectionHandle);
484	+	std::cout<<"==== Probe Vertex: " <<std::endl;
485	+	//printRecVtx(*splitVertexColl1->begin());
486	+	printRecVtxs(splitVertexCollection1Handle);
487	+	std::cout<<"==== Tag Vertex: " <<std::endl;
488	+	//printRecVtx(*splitVertexColl2->begin());
489	+	printRecVtxs(splitVertexCollection2Handle);
490	+	std::cout<<"==== Probe cluster information "<< std::endl;
491	+	printCluster(clusters1, beamSpot);
492	+	std::cout<<"==== Tag cluster information "<< std::endl;
493	+	printCluster(clusters2, beamSpot);
494	+	std::cout<<"**********Inefficiency Found************** " <<std::endl;
495	+	}
496	+	}
497		}
498		}
499
#	Line 548 | Line 615 | void PVEffAnalyzer::printSimTrks(const H
615		// vtx.ndof = 2*Sum(trackWeight) - 3
616
617		bool PVEffAnalyzer::isGoodSplitEvent( const reco::Vertex & org_vtx)
551	–	//, const reco::Vertex & tag_vtx, const reco::Vertex & probe_vtx)
618		{
619	<	return (  org_vtx.tracksSize() >= 6 && org_vtx.ndof() >= 6 );
620	<	//      tag_vtx.tracksSize() > = 2 && probe_vtx.tracksSize() > = 2      ) ;
619	>	if (  org_vtx.tracksSize() >= 6 && org_vtx.ndof() >= 6 && avgWeight(org_vtx) > 0.75  ) return true;
620	>	//if (  org_vtx.tracksSize() >= 6 && org_vtx.ndof() >= 6 ) return true;
621	>	else
622	>	return false;
623		}
624
625		// Comparing the tag vertex with the unsplit vertex position in z
#	Line 559 | Line 627 | bool PVEffAnalyzer::isGoodTagVertex( con
627		{
628		if(tag_vtx.isValid () && !tag_vtx.isFake()) {
629		float zsign = TMath::Abs(tag_vtx.z() - org_vtx.z() ) / TMath::Max(tag_vtx.zError(), org_vtx.zError() );
630	+	//float zsign = TMath::Abs(tag_vtx.z() - org_vtx.z() ) / sqrt(pow(tag_vtx.zError(),2)+pow(org_vtx.zError(),2) );
631		return (zsign < zsign_cut) ;
632		}
633		else
#	Line 570 | Line 639 | bool PVEffAnalyzer::isGoodProbeVertex( c
639		{
640		if(probe_vtx.isValid () && !probe_vtx.isFake()) {
641		float zsign = TMath::Abs(probe_vtx.z() - org_vtx.z() ) / TMath::Max(probe_vtx.zError(), org_vtx.zError() );
642	+	//float zsign = TMath::Abs(probe_vtx.z() - org_vtx.z() ) / sqrt(pow(probe_vtx.zError(),2)+pow(org_vtx.zError(),2) );
643		return (zsign < zsign_cut) ;
644		}
645		else
#	Line 587 | Line 657 | bool PVEffAnalyzer::isAssoVertex(const P
657		}
658
659
660	+
661		void PVEffAnalyzer::MakeEff(TH1D* numer, TH1D* denom, TH1D* eff,  //TGraphAsymmErrors* & gr_eff,
662		const bool rebin, const Float_t n ) {
663		eff->Divide( numer, denom, 1,1,"B" );
#	Line 622 | Line 693 | void PVEffAnalyzer::MakeEff(TH1D* numer,
693		}
694		}
695
696	+	bool PVEffAnalyzer::goodCluster(std::vector< std::vector<reco::TransientTrack> > & clusters)
697	+	{
698	+	if( clusters.size() == 1 && ( clusters.begin()->size() > 1 )) return true;
699	+	else
700	+	return false;
701	+	}
702	+
703	+	void PVEffAnalyzer::printCluster(std::vector< std::vector<reco::TransientTrack> > & clusters, const Point & beamSpot)
704	+	{
705	+	//std::cout<<"theTrackClusterizer.zSeparation()"<<theTrackClusterizer.zSeparation()<<std::endl;
706	+	cout << "#clusters " << std::setw(3) << clusters.size() << endl;
707	+	int idx_cluster = 0;
708	+	for (vector< vector<reco::TransientTrack> >::const_iterator iclus = clusters.begin(); iclus != clusters.end(); iclus++, idx_cluster++) {
709	+	if((*iclus).size() == 0) continue;
710	+	cout << "=== Cluster # " << std::setw(3) << idx_cluster << std::setw(3) << " nTracks" << std::setw(3) << (*iclus).size() <<endl;
711	+	for (int idx_track = 0; idx_track < (*iclus).size() ; idx_track++) {
712	+	reco::Track t = (*iclus).at(idx_track).track();
713	+	cout << "== Trk #  " << std::setw(4) << idx_track
714	+	<< " nHit     " << std::setw(4) << t.hitPattern().numberOfValidHits() <<endl
715	+	<< " pt       " << std::setw(9) << std::fixed << std::setprecision(4) << t.pt()
716	+	<< " dpt      " << std::setw(9) << t.ptError() << endl
717	+	<< " dxy(BS)  " << std::setw(9) << t.dxy(beamSpot)
718	+	<< " sigmaDxy " << std::setw(9) << t.dxyError() << endl
719	+	<< " dz(BS)   " << std::setw(9) << t.dz(beamSpot)
720	+	<< " sigmaDz  " << std::setw(9) << t.dzError() << endl
721	+	<< " xPCA     " << std::setw(9) << t.vertex().x()
722	+	<< " yPCA     " << std::setw(9) << t.vertex().y()
723	+	<< " zPCA     " << std::setw(9) << t.vertex().z()
724	+	<< endl;
725	+	}
726	+	}
727	+
728	+	}
729	+	double PVEffAnalyzer::avgWeight(const reco::Vertex & vtx)
730	+	{
731	+	if(!vtx.isValid() || vtx.isFake() )
732	+	return 0;
733	+	else
734	+	return (vtx.ndof()+3.0)/2.0/double(vtx.tracksSize());
735	+	}

Diff Legend

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