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.6 by yygao, Wed May 5 16:33:37 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		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	<
111	<
110	>	reqCluster_               = iConfig.getUntrackedParameter<bool>("reqCluster",false);
111	>
112		// Specify the data mode vector
113		if(realData_) datamode.push_back(0);
114		else {
#	Line 120 | Line 144 | PVEffAnalyzer::PVEffAnalyzer(const edm::
144
145		PVEffAnalyzer::~PVEffAnalyzer()
146		{
123	–
147		// do anything here that needs to be done at desctruction time
148		// (e.g. close files, deallocate resources etc.)
149		theFile->cd();
#	Line 134 | 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 221 | 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 232 | Line 255 | std::vector<PVEffAnalyzer::simPrimaryVer
255		return simpv;
256		}
257
235	–	std::vector<PVEffAnalyzer::simPrimaryVertex> PVEffAnalyzer::getSimPVs(const Handle<HepMCProduct> evtMC,
236	–	const Handle<SimVertexContainer> simVtxs,
237	–	const Handle<SimTrackContainer> simTrks)
238	–	{
239	–	// simvertices don't have enough information to decide,
240	–	// genVertices don't have the simulated coordinates  ( with VtxSmeared they might)
241	–	// go through simtracks to get the link between simulated and generated vertices
242	–	std::vector<PVEffAnalyzer::simPrimaryVertex> simpv;
243	–	int idx=0;
244	–	for(SimTrackContainer::const_iterator t=simTrks->begin();
245	–	t!=simTrks->end(); ++t){
246	–	if ( !(t->noVertex()) && !(t->type()==-99) ){
247	–	double ptsq=0;
248	–	bool primary=false;   // something coming directly from the primary vertex
249	–	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	–
278	–	const HepMC::FourVector & v=(*simVtxs)[t->vertIndex()].position();
279	–	//const HepLorentzVector & v=(*simVtxs)[t->vertIndex()].position();
280	–	if(primary or resonance){
281	–	{
282	–	// check all primaries found so far to avoid multiple entries
283	–	bool newVertex=true;
284	–	for(std::vector<PVEffAnalyzer::simPrimaryVertex>::iterator v0=simpv.begin();
285	–	v0!=simpv.end(); v0++){
286	–	if( (fabs(v0->x-v.x())<0.001) && (fabs(v0->y-v.y())<0.001) && (fabs(v0->z-v.z())<0.001) ){
287	–	if (track) {
288	–	v0->simTrackIndex.push_back(idx);
289	–	if (ptsq>(*v0).ptsq){(*v0).ptsq=ptsq;}
290	–	}
291	–	newVertex=false;
292	–	}
293	–	}
294	–	if(newVertex && !resonance){
295	–	PVEffAnalyzer::simPrimaryVertex anotherVertex(v.x(),v.y(),v.z());
296	–	if (track) anotherVertex.simTrackIndex.push_back(idx);
297	–	anotherVertex.ptsq=ptsq;
298	–	simpv.push_back(anotherVertex);
299	–	}
300	–	}//
301	–	}
302	–
303	–	}// simtrack has vertex and valid type
304	–	idx++;
305	–	}//simTrack loop
306	–	return simpv;
307	–	}
308	–
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 320 | Line 272 | PVEffAnalyzer::analyze(const edm::Event&
272		edm::LogInfo("Debug")<<"[PVEffAnalyzer]"<<endl;
273
274		// ====== TrackCollection
275	<	static const reco::TrackCollection s_empty_trackColl;
324	<	const reco::TrackCollection *trackColl = &s_empty_trackColl;
325	<	edm::Handle<reco::TrackCollection> trackCollectionHandle;
275	>	edm::Handle<edm::View<reco::Track> > trackCollectionHandle;
276		iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle);
277	<	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	<	}
277	>
278		// ====== splitTrackCollection1
279	<	static const reco::TrackCollection s_empty_splitTrackColl1;
334	<	const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
335	<	edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
279	>	edm::Handle<edm::View<reco::Track> > splitTrackCollection1Handle;
280		iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle);
281	<	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	<	}
281	>
282		// ====== splitTrackCollection2
283	<	static const reco::TrackCollection s_empty_splitTrackColl2;
344	<	const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
345	<	edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
283	>	edm::Handle<edm::View<reco::Track> > splitTrackCollection2Handle;
284		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	–	}
285
286		// ======= PrimaryVertexCollection
287		static const reco::VertexCollection s_empty_vertexColl;
#	Line 424 | 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,"");
428	<	//   simpv=getSimPVs(evtMC, simVtxs, simTrks);
360	>	simpv=getSimPVs(evtMC,"",ptcut_);
361		h_gen->Fill1d("nsimPV", simpv.size());
362
363		int isimpv = 0;
#	Line 443 | 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
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	>	// == 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	<	// ======= Analyze efficiency with split method
468	<	if ( !isGoodSplitEvent( *vertexColl->begin()) ) return;
469	<	if ( isGoodTagVertex( *splitVertexColl2->begin(), *vertexColl->begin(), zsigncut_) ) {
467	>	// clustering requirement
468	>	if( !goodCluster (clusters1) || !goodCluster (clusters2) ) return;
469	>
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 483 | Line 538 | PVEffAnalyzer::endJob() {
538		case 1: suffix = "_mct";
539		break;
540		}
541	<	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 579 | Line 638 | void PVEffAnalyzer::printSimTrks(const H
638		// vtx.ndof = 2*Sum(trackWeight) - 3
639
640		bool PVEffAnalyzer::isGoodSplitEvent( const reco::Vertex & org_vtx)
582	–	//, 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 590 | 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 601 | 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 618 | 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 648 | Line 711 | void PVEffAnalyzer::MakeEff(TH1D* numer,
711
712		for (int i = 1; i<eff->GetNbinsX()+1; i++) {
713		if(numer->GetBinContent(i) == 0 || denom->GetBinContent(i) == 0 ) continue;
714	<	float error = eff->GetBinContent(i)*sqrt(pow(numer->GetBinError(i)/numer->GetBinContent(i),2)
652	<	+pow(denom->GetBinError(i)/denom->GetBinContent(i),2));
714	>	float error = sqrt(eff->GetBinContent(i)*(1-eff->GetBinContent(i))/denom->GetBinContent(i)); // Binominal-Error
715		eff->SetBinError(i, error);
716		}
717		}
718
719	<	/*
720	<	// Copied from Validation/RecoTrack
721	<	void PVEffAnalyzer::setUpVectors(bool isRealData, int & nTrkMin_, int & nTrkMax_ ) {
722	<	std::vector<int>   totTAGv_ntracks, totPROBEv_ntracks;
723	<	if(!isRealData)
724	<	std::vector<int>   totSIMv_ntracks,totASSv_ntracks;
725	<	for (int k=nTrkMin_;k<nTrkMax_+1;k++) {
726	<	totTAGv_ntracks.push_back(0);
727	<	totPROBEv_ntracks.push_back(0);
728	<	if(!realData) {
729	<	totSIMv_ntracks.push_back(0);
730	<	totASSv_ntracks.push_back(0);
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
672	–	totTAG_ntracks.push_back(totTAGv_ntracks);
673	–	totPROBE_ntracks.push_back(totPROBEv_ntracks);
674	–	if(!realData) {
675	–	totSIM_ntracks.push_back(totSIMv_ntracks);
676	–	totASS_ntracks.push_back(totASSv_ntracks);
677	–	}
751		}
752	<
753	<	*/
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