ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/MitPhysics/Utils/src/JetTools.cc

Comparing UserCode/MitPhysics/Utils/src/JetTools.cc (file contents):
Revision 1.6 by mzanetti, Tue Jan 18 14:30:20 2011 UTC vs.
Revision 1.21 by pharris, Thu Apr 26 09:50:37 2012 UTC

#	Line 1 | Line 1
1		#include "MitPhysics/Utils/interface/JetTools.h"
2	+	#include <algorithm>
3	+	#include <vector>
4
5		ClassImp(mithep::JetTools)
6
7		using namespace mithep;
8	+
9
10		JetTools::JetTools()
11		{
#	Line 27 | Line 30 | Double_t JetTools::NJettiness(const Part
30		for(int j=0;j<int(jets->GetEntries());j++){
31		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
32		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-particles->At(i)->Eta()))
33	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),particles->At(i)->Phi()))));
33	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),particles->At(i)->Phi())))));
34		}
35		fval = fval + fvalpart;
36		}
#	Line 49 | Line 52 | Double_t JetTools::NJettiness(const PFCa
52		for(int j=0;j<int(jets->GetEntries());j++){
53		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
54		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-pfCandidates->At(i)->Eta()))
55	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi()))));
55	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi())))));
56		}
57		fval = fval + fvalpart;
58		}
#	Line 71 | Line 74 | Double_t JetTools::NJettiness(const Trac
74		for(int j=0;j<int(jets->GetEntries());j++){
75		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
76		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-tracks->At(i)->Eta()))
77	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),tracks->At(i)->Phi()))));
77	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),tracks->At(i)->Phi())))));
78		}
79		fval = fval + fvalpart;
80		}
#	Line 93 | Line 96 | Double_t JetTools::NJettiness(const JetO
96		for(int j=0;j<int(jets->GetEntries());j++){
97		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
98		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-jetsS->At(i)->Eta()))
99	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),jetsS->At(i)->Phi()))));
99	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),jetsS->At(i)->Phi())))));
100		}
101		fval = fval + fvalpart;
102		}
#	Line 115 | Line 118 | Double_t JetTools::NJettiness(const Calo
118		for(int j=0;j<int(jets->GetEntries());j++){
119		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
120		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-calos->At(i)->Eta()))
121	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),calos->At(i)->Phi()))));
121	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),calos->At(i)->Phi())))));
122		}
123		fval = fval + fvalpart;
124		}
#	Line 288 | Line 291 | Double_t JetTools::MtHiggs(const Particl
291		- met->Px()*dilepton->Px() - met->Py()*dilepton->Py()
292		- leptons->At(0)->Px()*leptons->At(1)->Px() - leptons->At(0)->Py()*leptons->At(1)->Py());
293		}
294	+	else if(nsel == 7){ // Use of M mass and mnu == 0
295	+	double deltaPhiDileptonMet = fabs(MathUtils::DeltaPhi(dilepton->Phi(),
296	+	met->Phi()));
297	+	mtHiggs = 2.0*dilepton->Pt()*met->Pt()*(1.0 - cos(deltaPhiDileptonMet));
298	+	}
299
300		if(nsel >= 0 && nsel <= 4){
301		mtHiggs = mll*mll + mnu*mnu + 2.0*(enell*enenn - enex*enex - eney*eney);
302		}
303	+
304		if(mtHiggs <= 0) mtHiggs = 0.0;
305		else             mtHiggs = TMath::Sqrt(mtHiggs);
306
#	Line 300 | Line 309 | Double_t JetTools::MtHiggs(const Particl
309		return mtHiggs;
310		}
311
312	<	void JetTools::Alpha(Double_t AlphaVar[2], const TrackCol *tracks, Jet *jet, const VertexCol *vertices, Double_t  delta_z, Double_t delta_cone){
313	<	AlphaVar[0] = -1.0;
314	<	AlphaVar[1] = -1.0;
306	<	if(tracks->GetEntries() <= 0) return;
312	>	Double_t JetTools::Beta(const TrackCol *tracks, Jet *jet, const Vertex *vertex, Double_t  delta_z, Double_t delta_cone){
313	>
314	>	if(tracks->GetEntries() <= 0) return 1.0;
315
316		double Pt_jets_X = 0. ;
317		double Pt_jets_Y = 0. ;
#	Line 314 | Line 322 | void JetTools::Alpha(Double_t AlphaVar[2
322		if(MathUtils::DeltaR(tracks->At(i)->Mom(),jet->Mom()) < delta_cone){
323		Pt_jets_X_tot += tracks->At(i)->Px();
324		Pt_jets_Y_tot += tracks->At(i)->Py();
325	<	double pDz = TMath::Abs(tracks->At(i)->DzCorrected(*vertices->At(0)));
325	>	double pDz = TMath::Abs(tracks->At(i)->DzCorrected(*vertex));
326		if(pDz < delta_z){
327		Pt_jets_X += tracks->At(i)->Px();
328		Pt_jets_Y += tracks->At(i)->Py();
#	Line 322 | Line 330 | void JetTools::Alpha(Double_t AlphaVar[2
330		}
331		}
332
333	<	if(jet->Pt() > 0)
334	<	AlphaVar[0] = sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / jet->Pt();
335	<	if(Pt_jets_X_tot > 0 || Pt_jets_Y_tot > 0)
336	<	AlphaVar[1] = sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot);
333	>	if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
334	>	return sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot);
335	>
336	>	return 1.0;
337	>	}
338	>
339	>
340	>	Double_t JetTools::Beta(const PFJet *jet, const Vertex *vertex, Double_t  delta_z){
341	>	double Pt_jets= 0. ;
342	>	double Pt_jetsTot = 0. ;
343	>
344	>	for(UInt_t i=0;i<jet->NPFCands();i++){
345	>	if(jet->PFCand(i)->TrackerTrk()){
346	>	Pt_jetsTot += jet->PFCand(i)->TrackerTrk()->Pt();
347	>	double pDz = TMath::Abs(jet->PFCand(i)->TrackerTrk()->DzCorrected(*vertex));
348	>	if(pDz < delta_z){
349	>	Pt_jets += jet->PFCand(i)->TrackerTrk()->Pt();
350	>	}
351	>	}
352	>	}
353	>
354	>	Double_t beta = 0.;
355	>	if (Pt_jetsTot > 0)
356	>	beta = Pt_jets/Pt_jetsTot;
357	>
358	>	return beta;
359	>	}
360	>
361	>	Double_t JetTools::Beta2(const PFJet *jet, const Vertex *vertex, Double_t  delta_z){
362	>	double Pt_jets= 0. ;
363	>	double Pt_jetsTot = 0. ;
364	>
365	>	for(UInt_t i=0;i<jet->NPFCands();i++){
366	>	if(jet->PFCand(i)->BestTrk()){
367	>	Pt_jetsTot += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
368	>	double pDz = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertex));
369	>	if(pDz < delta_z){
370	>	Pt_jets += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
371	>	}
372	>	}
373	>	}
374	>
375	>	Double_t beta = 1.0;
376	>	if (Pt_jetsTot > 0)
377	>	beta = Pt_jets/Pt_jetsTot;
378	>	return beta;
379	>	}
380	>
381	>
382	>	Bool_t  JetTools::PassBetaVertexAssociationCut(const PFJet *jet, const Vertex *referenceVertex, const VertexCol *vertices, Double_t delta_z) {
383	>
384	>	Bool_t passBetaCut = kTRUE;
385	>	if(vertices->GetEntries() > 0) {
386	>	Double_t Beta = JetTools::Beta(jet, referenceVertex, 0.2);
387	>	Double_t Beta_other = 0.0;
388	>	for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
389	>	if (referenceVertex == vertices->At(nv)) continue;
390	>	Double_t BetaAux = JetTools::Beta(jet, vertices->At(nv), 0.2);
391	>	if(BetaAux > Beta_other) Beta_other = BetaAux;
392	>	}
393	>	if(Beta_other > Beta) passBetaCut = kFALSE;
394	>	}
395	>
396	>	return passBetaCut;
397	>
398	>	}
399	>
400	>	Bool_t  JetTools::PassBeta2VertexAssociationCut(const PFJet *jet, const Vertex *referenceVertex, const VertexCol *vertices, Double_t delta_z) {
401	>
402	>	Bool_t passBetaCut = kTRUE;
403	>	if(vertices->GetEntries() > 0) {
404	>	Double_t Beta = JetTools::Beta2(jet, referenceVertex, 0.2);
405	>	Double_t Beta_other = 0.0;
406	>	for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
407	>	if (referenceVertex == vertices->At(nv)) continue;
408	>	Double_t BetaAux = JetTools::Beta2(jet, vertices->At(nv), 0.2);
409	>	if(BetaAux > Beta_other) Beta_other = BetaAux;
410	>	}
411	>	if(Beta_other > Beta) passBetaCut = kFALSE;
412	>	}
413	>
414	>	return passBetaCut;
415	>
416	>	}
417	>
418	>
419	>	Int_t JetTools::MaxBetaVertexIndex(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){
420	>
421	>	Int_t vertexIndex = -1;
422	>	double beta = -0.1;
423	>	for (UInt_t v=0; v < vertices->GetEntries(); v++){
424	>	Double_t betaTmp = JetTools::Beta(jet, vertices->At(v), delta_z);
425	>	if (betaTmp > beta) {
426	>	beta = betaTmp;
427	>	vertexIndex = v;
428	>	}
429	>	}
430	>	return vertexIndex;
431	>
432	>	}
433	>
434	>	Int_t JetTools::MaxBeta2VertexIndex(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){
435	>
436	>	Int_t vertexIndex = -1;
437	>	double beta = -0.1;
438	>	for (UInt_t v=0; v < vertices->GetEntries(); v++){
439	>	Double_t betaTmp = JetTools::Beta2(jet, vertices->At(v), delta_z);
440	>	if (betaTmp > beta) {
441	>	beta = betaTmp;
442	>	vertexIndex = v;
443	>	}
444	>	}
445	>	return vertexIndex;
446	>
447		}
448
449	+
450	+	Int_t JetTools::JetToPVAssociation(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){
451	+
452	+	std::vector<float> verticesPt2(vertices->GetEntries());
453	+	for(UInt_t i=0;i<jet->NPFCands();i++){
454	+	if(jet->PFCand(i)->BestTrk()){
455	+	double minDZ = delta_z;
456	+	int trackVertexIndex = -1;
457	+	for (UInt_t v=0; v < vertices->GetEntries(); v++){
458	+	if (minDZ > TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)))) {
459	+	minDZ = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)));
460	+	trackVertexIndex = v;
461	+	}
462	+	}
463	+	if (trackVertexIndex < 0) continue;
464	+	verticesPt2[trackVertexIndex]+= jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
465	+	}
466	+	}
467	+
468	+	Int_t vertexIndex = 0;
469	+	float pt2Max = 0;
470	+	for (uint i=0; i < verticesPt2.size(); ++i){
471	+	if (pt2Max < verticesPt2[i]) {
472	+	pt2Max = verticesPt2[i];
473	+	vertexIndex = i;
474	+	}
475	+	}
476	+	return vertexIndex;
477	+	}
478	+	const PFCandidate* JetTools::leadCand(const PFJet *iJet,int iPFType,bool i2nd) {
479	+	int lCount = 0;
480	+	const PFCandidate *lCand = 0;
481	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
482	+	lCand = iJet->PFCand(i0);
483	+	if(iPFType != -1 && lCand->PFType() != iPFType) continue;
484	+	if(lCount == 0 && !i2nd) break;
485	+	if(lCount >  0)          break;
486	+	lCount++;
487	+	}
488	+	return lCand;
489	+	}
490	+	Double_t JetTools::impactParameter(const PFJet *iJet,const Vertex *iVertex,bool iDZ) {
491	+	double lDZCorr = -1000;
492	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
493	+	const PFCandidate *pCand = iJet->PFCand(i0);
494	+	if(pCand->Trk() == 0) continue;
495	+	//if(pCand->Pt() < 1.) continue; => previous iterations
496	+	if(iDZ)  lDZCorr = fabs(pCand->TrackerTrk()->DzCorrected(*iVertex));
497	+	if(!iDZ) lDZCorr = fabs(pCand->TrackerTrk()->D0Corrected(*iVertex));
498	+	break;
499	+	}
500	+	return lDZCorr;
501	+	}
502	+	Double_t JetTools::dRMean(const PFJet *iJet,int iPFType) {
503	+	double lDRMean = 0;
504	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
505	+	const PFCandidate *pCand = iJet->PFCand(i0);
506	+	if(iPFType != -1 && pCand->PFType() != iPFType) continue;
507	+	double pDR = MathUtils::DeltaR(iJet->Mom(),pCand->Mom());
508	+	lDRMean    += pDR*(pCand->Pt())/iJet->RawMom().Pt();
509	+	}
510	+	return lDRMean;
511	+	}
512	+	Double_t JetTools::frac(const PFJet *iJet,Double_t iDRMax,Double_t iDRMin,Int_t iPFType) {
513	+	double lFrac = 0;
514	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
515	+	const PFCandidate *pCand = iJet->PFCand(i0);
516	+	if(iPFType != -1 && pCand->PFType() != iPFType) continue;
517	+	Double_t pDR = MathUtils::DeltaR(iJet->Mom(),pCand->Mom());
518	+	if(pDR > iDRMax) continue;
519	+	if(pDR < iDRMax-0.1) continue;
520	+	lFrac += pCand->Pt()/iJet->Pt();
521	+	}
522	+	return lFrac;
523	+	}
524	+	Double_t JetTools::betaStar(const PFJet *iJet,const Vertex *iVertex,const VertexCol* iVertices,Double_t iDZCut) {
525	+	Double_t lTotal = 0;
526	+	Double_t lPileup = 0;
527	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
528	+	const PFCandidate* pPF   = iJet->PFCand(i0);
529	+	const Track* pTrack      = pPF->TrackerTrk();
530	+	//if(pPF->GsfTrk()) pTrack = pPF->GsfTrk(); ==> not used in CMSSW
531	+	if(pTrack == 0) continue;
532	+	lTotal += pPF->Pt();
533	+	double pDZPV  = fabs(pTrack->DzCorrected(*iVertex));
534	+	double pDZMin = pDZPV;
535	+	for(unsigned int i1 = 0; i1 < iVertices->GetEntries(); i1++) {
536	+	const Vertex *pV = iVertices->At(i1);
537	+	if(pV->Ndof() < 4 ||
538	+	(pV->Position() - iVertex->Position()).R() > 0.02 ) continue;
539	+	pDZMin = TMath::Min(pDZMin,fabs(pTrack->DzCorrected(*pV)));
540	+	}
541	+	if(pDZPV > 0.2 && pDZMin < 0.2) lPileup += pTrack->Pt();
542	+	}
543	+	if(lTotal == 0) lTotal = 1;
544	+	return lPileup/(lTotal);
545	+	}
546	+	Bool_t  JetTools::passPFLooseId(const PFJet *iJet) {
547	+	if(iJet->E()                              == 0)       return false;
548	+	if(iJet->NeutralHadronEnergy()/iJet->E()  >  0.99)    return false;
549	+	if(iJet->NeutralEmEnergy()/iJet->E()      >  0.99)    return false;
550	+	if(iJet->NConstituents()                  <  2)       return false;
551	+	if(iJet->ChargedHadronEnergy()/iJet->E()  <= 0     && fabs(iJet->Eta()) < 2.4 ) return false;
552	+	if(iJet->ChargedEmEnergy()/iJet->E()      >  0.99  && fabs(iJet->Eta()) < 2.4 ) return false;
553	+	if(iJet->ChargedMultiplicity()            < 1      && fabs(iJet->Eta()) < 2.4 ) return false;
554	+	return true;
555	+	}
556	+
557	+	/*
558	+	double JetTools::genFrac(const PFJet *iJet) {
559	+	double lTrueFrac = 0;
560	+	for(UInt_t i0 = 0; i0 < fParticles->GetEntries(); i0++) {
561	+	const MCParticle *p = fParticles->At(i0);
562	+	if(p->Status() != 1) continue;
563	+	double pDEta = iJet->Eta() - p->Eta();
564	+	double pDPhi = fabs(iJet->Phi()-p->Phi()); if(pDPhi > 2.*TMath::Pi() - pDPhi) pDPhi =  2.*TMath::Pi() - pDPhi;
565	+	double pDR   = sqrt(pDEta*pDEta + pDPhi*pDPhi);
566	+	if(pDR > 0.5) continue;
567	+	lTrueFrac += p->Pt();
568	+	}
569	+	lTrueFrac/=iJet->Pt();
570	+	return lTrueFrac;
571	+	}
572	+	*/

Diff Legend

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