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.3 by yygao, Tue Apr 20 10:15:43 2010 UTC vs.
Revision 1.6 by yygao, Wed May 5 16:33:37 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");
106		zsigncut_                  = iConfig.getUntrackedParameter<double>("zsigncut");
107	+	ptcut_                     = iConfig.getUntrackedParameter<double>("ptcut");
108		analyze_                   = iConfig.getUntrackedParameter<bool>("analyze",false);
109		bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
110		reqCluster_               = iConfig.getUntrackedParameter<bool>("reqCluster",false);
#	Line 145 | Line 157 | PVEffAnalyzer::~PVEffAnalyzer()
157		//
158		// member functions
159		//
160	<	std::vector<PVEffAnalyzer::simPrimaryVertex> PVEffAnalyzer::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
160	>	std::vector<PVEffAnalyzer::simPrimaryVertex> PVEffAnalyzer::getSimPVs(const Handle<HepMCProduct> & evtMC, std::string suffix, double & ptcut_)
161		{
162		std::vector<PVEffAnalyzer::simPrimaryVertex> simpv;
163		const HepMC::GenEvent* evt=evtMC->GetEvent();
#	Line 232 | Line 244 | std::vector<PVEffAnalyzer::simPrimaryVer
244		vp->ptot.setPz(vp->ptot.pz()+m.pz());
245		vp->ptot.setE(vp->ptot.e()+m.e());
246		vp->ptsq+=(m.perp())*(m.perp());
247	<	if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
247	>	if ( (m.perp()>ptcut_) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
248		vp->nGenTrk++;
249		}
250		}
#	Line 243 | Line 255 | std::vector<PVEffAnalyzer::simPrimaryVer
255		return simpv;
256		}
257
246	–	std::vector<PVEffAnalyzer::simPrimaryVertex> PVEffAnalyzer::getSimPVs(const Handle<HepMCProduct> evtMC,
247	–	const Handle<SimVertexContainer> simVtxs,
248	–	const Handle<SimTrackContainer> simTrks)
249	–	{
250	–	// simvertices don't have enough information to decide,
251	–	// genVertices don't have the simulated coordinates  ( with VtxSmeared they might)
252	–	// go through simtracks to get the link between simulated and generated vertices
253	–	std::vector<PVEffAnalyzer::simPrimaryVertex> simpv;
254	–	int idx=0;
255	–	for(SimTrackContainer::const_iterator t=simTrks->begin();
256	–	t!=simTrks->end(); ++t){
257	–	if ( !(t->noVertex()) && !(t->type()==-99) ){
258	–	double ptsq=0;
259	–	bool primary=false;   // something coming directly from the primary vertex
260	–	bool resonance=false; // resonance
261	–	bool track=false;     // undecayed, charged particle
262	–	HepMC::GenParticle* gp=evtMC->GetEvent()->barcode_to_particle( (*t).genpartIndex() );
263	–	if (gp) {
264	–	HepMC::GenVertex * gv=gp->production_vertex();
265	–	if (gv) {
266	–	for ( HepMC::GenVertex::particle_iterator
267	–	daughter =  gv->particles_begin(HepMC::descendants);
268	–	daughter != gv->particles_end(HepMC::descendants);
269	–	++daughter ) {
270	–	if (isFinalstateParticle(*daughter)){
271	–	ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
272	–	}
273	–	}
274	–	//primary =  ( gv->position().t()==0 );
275	–	primary = true;
276	–	h_gen->Fill1d("genPart_cT", gv->position().t()); // mm
277	–	h_gen->Fill1d("genPart_T", gv->position().t()/299.792458); // ns
278	–
279	–	//resonance= ( gp->mother() && isResonance(gp->mother()));  // in CLHEP/HepMC days
280	–	// no more mother pointer in the improved HepMC GenParticle
281	–	resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
282	–	if (gp->status()==1){
283	–	//track=((pdt->particle(gp->pdg_id()))->charge() != 0);
284	–	track=not isCharged(gp);
285	–	}
286	–	}
287	–	}
288	–
289	–	const HepMC::FourVector & v=(*simVtxs)[t->vertIndex()].position();
290	–	//const HepLorentzVector & v=(*simVtxs)[t->vertIndex()].position();
291	–	if(primary or resonance){
292	–	{
293	–	// check all primaries found so far to avoid multiple entries
294	–	bool newVertex=true;
295	–	for(std::vector<PVEffAnalyzer::simPrimaryVertex>::iterator v0=simpv.begin();
296	–	v0!=simpv.end(); v0++){
297	–	if( (fabs(v0->x-v.x())<0.001) && (fabs(v0->y-v.y())<0.001) && (fabs(v0->z-v.z())<0.001) ){
298	–	if (track) {
299	–	v0->simTrackIndex.push_back(idx);
300	–	if (ptsq>(*v0).ptsq){(*v0).ptsq=ptsq;}
301	–	}
302	–	newVertex=false;
303	–	}
304	–	}
305	–	if(newVertex && !resonance){
306	–	PVEffAnalyzer::simPrimaryVertex anotherVertex(v.x(),v.y(),v.z());
307	–	if (track) anotherVertex.simTrackIndex.push_back(idx);
308	–	anotherVertex.ptsq=ptsq;
309	–	simpv.push_back(anotherVertex);
310	–	}
311	–	}//
312	–	}
313	–
314	–	}// simtrack has vertex and valid type
315	–	idx++;
316	–	}//simTrack loop
317	–	return simpv;
318	–	}
319	–
258		// ------------ method called to for each event  ------------
259		void
260		PVEffAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
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 331 | Line 272 | PVEffAnalyzer::analyze(const edm::Event&
272		edm::LogInfo("Debug")<<"[PVEffAnalyzer]"<<endl;
273
274		// ====== TrackCollection
275	<	static const reco::TrackCollection s_empty_trackColl;
335	<	const reco::TrackCollection *trackColl = &s_empty_trackColl;
336	<	edm::Handle<reco::TrackCollection> trackCollectionHandle;
275	>	edm::Handle<edm::View<reco::Track> > trackCollectionHandle;
276		iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle);
277	<	if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
339	<	trackColl = trackCollectionHandle.product();
340	<	} else {
341	<	edm::LogInfo("Debug") << "[PVEffAnalyzer] trackCollection cannot be found -> using empty collection of same type." <<endl;
342	<	}
277	>
278		// ====== splitTrackCollection1
279	<	static const reco::TrackCollection s_empty_splitTrackColl1;
345	<	const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
346	<	edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
279	>	edm::Handle<edm::View<reco::Track> > splitTrackCollection1Handle;
280		iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle);
281	<	if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
349	<	splitTrackColl1 = splitTrackCollection1Handle.product();
350	<	} else {
351	<	edm::LogInfo("Debug") << "[PVEffAnalyzer] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
352	<	}
281	>
282		// ====== splitTrackCollection2
283	<	static const reco::TrackCollection s_empty_splitTrackColl2;
355	<	const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
356	<	edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
283	>	edm::Handle<edm::View<reco::Track> > splitTrackCollection2Handle;
284		iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle);
358	–	if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
359	–	splitTrackColl2 = splitTrackCollection2Handle.product();
360	–	} else {
361	–	edm::LogInfo("Debug") << "[PVEffAnalyzer] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
362	–	}
285
286		// ======= PrimaryVertexCollection
287		static const reco::VertexCollection s_empty_vertexColl;
#	Line 435 | Line 357 | PVEffAnalyzer::analyze(const edm::Event&
357		TString suffix = "_mct";
358		// make a list of primary vertices:
359		std::vector<simPrimaryVertex> simpv;
360	<	simpv=getSimPVs(evtMC,"");
439	<	//   simpv=getSimPVs(evtMC, simVtxs, simTrks);
360	>	simpv=getSimPVs(evtMC,"",ptcut_);
361		h_gen->Fill1d("nsimPV", simpv.size());
362
363		int isimpv = 0;
#	Line 454 | 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));
381	<	}
382	<	} // End of Analyzing MC Efficiency
383	<
384	<	// Get the Builder for TrackClusters
385	<	edm::ESHandle<TransientTrackBuilder> theB;
386	<	iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB);
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	>	} // == end of getting the ntracks
423	>
424	>	// == Fill Efficiency Histogram
425	>	if(ntracks>1) {
426	>	h_summary->Fill1d(TString("eff_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("eff_numer_ntrack"+suffix), ntracks);
431	>	}
432	>
433	>	// == Fill FakeRate Histogram
434	>	if(!vertexColl->begin()->isFake ()) {
435	>	h_summary->Fill1d(TString("fakerate_denom_ntrack"+suffix), int(vertexColl->begin()->tracksSize()));
436	>	if (!isAssoVertex(vsim, *vertexColl->begin(), zsigncut_)  )
437	>	h_summary->Fill1d(TString("fakerate_numer_ntrack"+suffix), int(vertexColl->begin()->tracksSize()));
438	>	}
439	>	}  // end of MC analyzing
440	>	} // end of MC access
441
442	<	std::vector< std::vector<reco::TransientTrack> > clusters1 =  theTrackClusterizer.clusterize((*theB).build(splitTrackCollection1Handle));
443	<	std::vector< std::vector<reco::TransientTrack> > clusters2 =  theTrackClusterizer.clusterize((*theB).build(splitTrackCollection2Handle));
442	>	// == Start of Split-Method
443	>	// Event selections
444	>	if ( !isGoodSplitEvent( *vertexColl->begin())) return;
445	>	for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
446	>	t!=vertexColl->begin()->tracks_end(); t++) {
447	>	h_summary->Fill1d("trackweight",vertexColl->begin()->trackWeight(*t));
448	>	}
449	>
450	>	// == Start of FakeRate Analyzing
451	>
452	>	if(!splitVertexColl1->begin()->isFake()) {
453	>	h_summary->Fill1d("fakerate_denom_ntrack", int(splitVertexColl1->begin()->tracksSize()));
454	>	if ( !isGoodProbeVertex( *splitVertexColl1->begin(), *vertexColl->begin(), zsigncut_) )
455	>	h_summary->Fill1d("fakerate_numer_ntrack", int(splitVertexColl1->begin()->tracksSize()));
456	>	}
457	>
458	>	// == End of FakeRate Analyzing
459	>	// == Start of Efficiency Analyzing
460	>	std::vector< std::vector<reco::TransientTrack> > clusters1;
461	>	std::vector< std::vector<reco::TransientTrack> > clusters2;
462	>	if(splitTrackCollection1Handle.isValid())
463	>	clusters1 =  theTrackClusterizer.clusterize((*theB).build(splitTrackCollection1Handle));
464	>	if(splitTrackCollection2Handle.isValid())
465	>	clusters2 =  theTrackClusterizer.clusterize((*theB).build(splitTrackCollection2Handle));
466
467	<	bool withValidCluster_tag = false;
468	<	bool withValidCluster_probe = false;
472	<
473	<	for (vector< vector<reco::TransientTrack> >::const_iterator iclus = clusters1.begin(); iclus != clusters1.end(); iclus++) {
474	<	if((*iclus).size()>1)  withValidCluster_probe = true;
475	<	}
476	<	for (vector< vector<reco::TransientTrack> >::const_iterator iclus = clusters2.begin(); iclus != clusters2.end(); iclus++) {
477	<	if((*iclus).size()>1)  withValidCluster_tag = true;
478	<	}
467	>	// clustering requirement
468	>	if( !goodCluster (clusters1) || !goodCluster (clusters2) ) return;
469
470	<	// ======= Analyze efficiency with split method
471	<	if( reqCluster_ ) {
482	<	if ( !isGoodSplitEvent( *vertexColl->begin()) || !withValidCluster_tag || !withValidCluster_probe ) return;
483	<	}
484	<	else {
485	<	if ( !isGoodSplitEvent( *vertexColl->begin())) return;
486	<	}
487	<	if ( isGoodTagVertex( *splitVertexColl2->begin(), *vertexColl->begin(), zsigncut_) ) {
470	>	// TagVertex Requirement
471	>	if ( isGoodTagVertex( *(splitVertexColl2->begin()), *(vertexColl->begin()), zsigncut_) ) {
472		if(verbose_) {
473	<	std::cout<<"splitTrackColl1->size() = " << int(splitTrackColl1->size()) << std::endl;
474	<	std::cout<<"splitTrackColl2->size() = " << int(splitTrackColl2->size()) << std::endl;
473	>	std::cout<<"splitTrackCollection1Handle->size() = " << int(splitTrackCollection1Handle->size()) << std::endl;
474	>	std::cout<<"splitTrackCollection2Handle->size() = " << int(splitTrackCollection2Handle->size()) << std::endl;
475	>	}
476	>	int nGoodTracks = 0;
477	>	for(unsigned int i = 0; i < splitTrackCollection1Handle->size(); i++ ) {
478	>	edm::RefToBase<reco::Track> track(splitTrackCollection1Handle, i);
479	>	if(track->pt() > ptcut_) // && abs(track->eta()) < 2.5 )
480	>	nGoodTracks++;
481		}
482	<	h_summary->Fill1d("denom_ntrack", int(splitTrackColl1->size()));
483	<	if ( isGoodProbeVertex( *splitVertexColl1->begin(), *vertexColl->begin(), zsigncut_) )
484	<	h_summary->Fill1d("numer_ntrack", int(splitTrackColl1->size()) );
482	>
483	>	h_summary->Fill1d("eff_denom_ntrack", nGoodTracks);
484	>	// ProbeVertex Requirement
485	>	if ( isGoodProbeVertex( *splitVertexColl1->begin(), *vertexColl->begin(), zsigncut_) ) {
486	>	h_summary->Fill1d("eff_numer_ntrack", nGoodTracks);
487	>	h_summary->Fill1d("avgWeight_orgvtx_eff",avgWeight(*vertexColl->begin()) );
488	>	h_summary->Fill1d("avgWeight_tagvtx_eff",avgWeight(*splitVertexColl2->begin()) );
489	>	h_summary->Fill1d("avgWeight_probevtx_eff",avgWeight(*splitVertexColl1->begin()) );
490	>	}
491	>	else
492	>	{
493	>	if (nGoodTracks >= 2) {
494	>	h_summary->Fill1d("avgWeight_orgvtx_ineff",avgWeight(*vertexColl->begin()) );
495	>	h_summary->Fill1d("avgWeight_tagvtx_ineff",avgWeight(*splitVertexColl2->begin()) );
496	>	h_summary->Fill1d("avgWeight_probevtx_ineff",avgWeight(*splitVertexColl1->begin()) );
497	>
498	>	std::cout<<"**********Inefficiency Found************** " <<std::endl;
499	>	std::cout<<"==== Original Vertex: " <<std::endl;
500	>	//printRecVtx(*vertexColl->begin());
501	>	printRecVtxs(vertexCollectionHandle);
502	>	std::cout<<"==== Probe Vertex: " <<std::endl;
503	>	//printRecVtx(*splitVertexColl1->begin());
504	>	printRecVtxs(splitVertexCollection1Handle);
505	>	std::cout<<"==== Tag Vertex: " <<std::endl;
506	>	//printRecVtx(*splitVertexColl2->begin());
507	>	printRecVtxs(splitVertexCollection2Handle);
508	>	std::cout<<"==== Probe cluster information "<< std::endl;
509	>	printCluster(clusters1, beamSpot);
510	>	std::cout<<"==== Tag cluster information "<< std::endl;
511	>	printCluster(clusters2, beamSpot);
512	>	std::cout<<"**********Inefficiency Found************** " <<std::endl;
513	>	}
514	>	}
515		}
516	+	// == End of Efficiency Analyzing
517	+
518		}
519
520
#	Line 516 | Line 538 | PVEffAnalyzer::endJob() {
538		case 1: suffix = "_mct";
539		break;
540		}
541	<	if (analyze_)
542	<	MakeEff(h_summary->ReadHisto1D(TString("numer_ntrack"+suffix)), h_summary->ReadHisto1D(TString("denom_ntrack"+suffix)), h_summary->ReadHisto1D(TString("eff_ntrack"+suffix)), false, 1);
541	>	if (analyze_) {
542	>	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);
543	>	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);
544	>
545	>	}
546		}
547		}
548
#	Line 613 | Line 638 | void PVEffAnalyzer::printSimTrks(const H
638		// vtx.ndof = 2*Sum(trackWeight) - 3
639
640		bool PVEffAnalyzer::isGoodSplitEvent( const reco::Vertex & org_vtx)
616	–	//, const reco::Vertex & tag_vtx, const reco::Vertex & probe_vtx)
641		{
642	<	return (  org_vtx.tracksSize() >= 6 && org_vtx.ndof() >= 6 );
643	<	//      tag_vtx.tracksSize() > = 2 && probe_vtx.tracksSize() > = 2      ) ;
642	>	if (  org_vtx.tracksSize() >= 6 && org_vtx.ndof() >= 6 && avgWeight(org_vtx) > 0.75  ) return true;
643	>	//if (  org_vtx.tracksSize() >= 6 && org_vtx.ndof() >= 6 ) return true;
644	>	else
645	>	return false;
646		}
647
648		// Comparing the tag vertex with the unsplit vertex position in z
#	Line 624 | Line 650 | bool PVEffAnalyzer::isGoodTagVertex( con
650		{
651		if(tag_vtx.isValid () && !tag_vtx.isFake()) {
652		float zsign = TMath::Abs(tag_vtx.z() - org_vtx.z() ) / TMath::Max(tag_vtx.zError(), org_vtx.zError() );
653	+	//float zsign = TMath::Abs(tag_vtx.z() - org_vtx.z() ) / sqrt(pow(tag_vtx.zError(),2)+pow(org_vtx.zError(),2) );
654		return (zsign < zsign_cut) ;
655		}
656		else
#	Line 635 | Line 662 | bool PVEffAnalyzer::isGoodProbeVertex( c
662		{
663		if(probe_vtx.isValid () && !probe_vtx.isFake()) {
664		float zsign = TMath::Abs(probe_vtx.z() - org_vtx.z() ) / TMath::Max(probe_vtx.zError(), org_vtx.zError() );
665	+	//float zsign = TMath::Abs(probe_vtx.z() - org_vtx.z() ) / sqrt(pow(probe_vtx.zError(),2)+pow(org_vtx.zError(),2) );
666		return (zsign < zsign_cut) ;
667		}
668		else
#	Line 652 | Line 680 | bool PVEffAnalyzer::isAssoVertex(const P
680		}
681
682
683	+
684		void PVEffAnalyzer::MakeEff(TH1D* numer, TH1D* denom, TH1D* eff,  //TGraphAsymmErrors* & gr_eff,
685		const bool rebin, const Float_t n ) {
686		eff->Divide( numer, denom, 1,1,"B" );
#	Line 686 | Line 715 | void PVEffAnalyzer::MakeEff(TH1D* numer,
715		eff->SetBinError(i, error);
716		}
717		}
718	+
719	+	bool PVEffAnalyzer::goodCluster(std::vector< std::vector<reco::TransientTrack> > & clusters)
720	+	{
721	+	if( clusters.size() == 1 && ( clusters.begin()->size() > 1 )) return true;
722	+	else
723	+	return false;
724	+	}
725	+
726	+	void PVEffAnalyzer::printCluster(std::vector< std::vector<reco::TransientTrack> > & clusters, const Point & beamSpot)
727	+	{
728	+	//std::cout<<"theTrackClusterizer.zSeparation()"<<theTrackClusterizer.zSeparation()<<std::endl;
729	+	cout << "#clusters " << std::setw(3) << clusters.size() << endl;
730	+	int idx_cluster = 0;
731	+	for (vector< vector<reco::TransientTrack> >::const_iterator iclus = clusters.begin(); iclus != clusters.end(); iclus++, idx_cluster++) {
732	+	if((*iclus).size() == 0) continue;
733	+	cout << "=== Cluster # " << std::setw(3) << idx_cluster << std::setw(3) << " nTracks" << std::setw(3) << (*iclus).size() <<endl;
734	+	for (int idx_track = 0; idx_track < (*iclus).size() ; idx_track++) {
735	+	reco::Track t = (*iclus).at(idx_track).track();
736	+	cout << "== Trk #  " << std::setw(4) << idx_track
737	+	<< " nHit     " << std::setw(4) << t.hitPattern().numberOfValidHits() <<endl
738	+	<< " pt       " << std::setw(9) << std::fixed << std::setprecision(4) << t.pt()
739	+	<< " dpt      " << std::setw(9) << t.ptError() << endl
740	+	<< " dxy(BS)  " << std::setw(9) << t.dxy(beamSpot)
741	+	<< " sigmaDxy " << std::setw(9) << t.dxyError() << endl
742	+	<< " dz(BS)   " << std::setw(9) << t.dz(beamSpot)
743	+	<< " sigmaDz  " << std::setw(9) << t.dzError() << endl
744	+	<< " xPCA     " << std::setw(9) << t.vertex().x()
745	+	<< " yPCA     " << std::setw(9) << t.vertex().y()
746	+	<< " zPCA     " << std::setw(9) << t.vertex().z()
747	+	<< endl;
748	+	}
749	+	}
750	+
751	+	}
752	+	double PVEffAnalyzer::avgWeight(const reco::Vertex & vtx)
753	+	{
754	+	if(!vtx.isValid() || vtx.isFake() )
755	+	return 0;
756	+	else
757	+	return (vtx.ndof()+3.0)/2.0/double(vtx.tracksSize());
758	+	}

Diff Legend

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