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.33 by pharris, Thu Mar 21 17:50:13 2013 UTC

#	Line 1 | Line 1
1		#include "MitPhysics/Utils/interface/JetTools.h"
2	+	#include <algorithm>
3	+	#include <vector>
4	+	#include "TMatrixDSym.h"
5	+	#include "TMatrixDSymEigen.h"
6
7		ClassImp(mithep::JetTools)
8
9		using namespace mithep;
10	+
11
12		JetTools::JetTools()
13		{
#	Line 27 | Line 32 | Double_t JetTools::NJettiness(const Part
32		for(int j=0;j<int(jets->GetEntries());j++){
33		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
34		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-particles->At(i)->Eta()))
35	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),particles->At(i)->Phi()))));
35	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),particles->At(i)->Phi())))));
36		}
37		fval = fval + fvalpart;
38		}
#	Line 49 | Line 54 | Double_t JetTools::NJettiness(const PFCa
54		for(int j=0;j<int(jets->GetEntries());j++){
55		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
56		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-pfCandidates->At(i)->Eta()))
57	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi()))));
57	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi())))));
58		}
59		fval = fval + fvalpart;
60		}
#	Line 71 | Line 76 | Double_t JetTools::NJettiness(const Trac
76		for(int j=0;j<int(jets->GetEntries());j++){
77		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
78		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-tracks->At(i)->Eta()))
79	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),tracks->At(i)->Phi()))));
79	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),tracks->At(i)->Phi())))));
80		}
81		fval = fval + fvalpart;
82		}
#	Line 93 | Line 98 | Double_t JetTools::NJettiness(const JetO
98		for(int j=0;j<int(jets->GetEntries());j++){
99		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
100		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-jetsS->At(i)->Eta()))
101	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),jetsS->At(i)->Phi()))));
101	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),jetsS->At(i)->Phi())))));
102		}
103		fval = fval + fvalpart;
104		}
#	Line 115 | Line 120 | Double_t JetTools::NJettiness(const Calo
120		for(int j=0;j<int(jets->GetEntries());j++){
121		fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
122		(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-calos->At(i)->Eta()))
123	<	- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),calos->At(i)->Phi()))));
123	>	- 2 * TMath::Cos(fabs(MathUtils::DeltaPhi(jets->At(j)->Phi(),calos->At(i)->Phi())))));
124		}
125		fval = fval + fvalpart;
126		}
#	Line 288 | Line 293 | Double_t JetTools::MtHiggs(const Particl
293		- met->Px()*dilepton->Px() - met->Py()*dilepton->Py()
294		- leptons->At(0)->Px()*leptons->At(1)->Px() - leptons->At(0)->Py()*leptons->At(1)->Py());
295		}
296	+	else if(nsel == 7){ // Use of M mass and mnu == 0
297	+	double deltaPhiDileptonMet = fabs(MathUtils::DeltaPhi(dilepton->Phi(),
298	+	met->Phi()));
299	+	mtHiggs = 2.0*dilepton->Pt()*met->Pt()*(1.0 - cos(deltaPhiDileptonMet));
300	+	}
301
302		if(nsel >= 0 && nsel <= 4){
303		mtHiggs = mll*mll + mnu*mnu + 2.0*(enell*enenn - enex*enex - eney*eney);
304		}
305	+
306		if(mtHiggs <= 0) mtHiggs = 0.0;
307		else             mtHiggs = TMath::Sqrt(mtHiggs);
308
#	Line 300 | Line 311 | Double_t JetTools::MtHiggs(const Particl
311		return mtHiggs;
312		}
313
314	<	void JetTools::Alpha(Double_t AlphaVar[2], const TrackCol *tracks, Jet *jet, const VertexCol *vertices, Double_t  delta_z, Double_t delta_cone){
315	<	AlphaVar[0] = -1.0;
316	<	AlphaVar[1] = -1.0;
306	<	if(tracks->GetEntries() <= 0) return;
314	>	Double_t JetTools::Beta(const TrackCol *tracks, Jet *jet, const Vertex *vertex, Double_t  delta_z, Double_t delta_cone){
315	>
316	>	if(tracks->GetEntries() <= 0) return 1.0;
317
318		double Pt_jets_X = 0. ;
319		double Pt_jets_Y = 0. ;
#	Line 314 | Line 324 | void JetTools::Alpha(Double_t AlphaVar[2
324		if(MathUtils::DeltaR(tracks->At(i)->Mom(),jet->Mom()) < delta_cone){
325		Pt_jets_X_tot += tracks->At(i)->Px();
326		Pt_jets_Y_tot += tracks->At(i)->Py();
327	<	double pDz = TMath::Abs(tracks->At(i)->DzCorrected(*vertices->At(0)));
327	>	double pDz = TMath::Abs(tracks->At(i)->DzCorrected(*vertex));
328		if(pDz < delta_z){
329		Pt_jets_X += tracks->At(i)->Px();
330		Pt_jets_Y += tracks->At(i)->Py();
#	Line 322 | Line 332 | void JetTools::Alpha(Double_t AlphaVar[2
332		}
333		}
334
335	<	if(jet->Pt() > 0)
336	<	AlphaVar[0] = sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / jet->Pt();
337	<	if(Pt_jets_X_tot > 0 || Pt_jets_Y_tot > 0)
338	<	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);
335	>	if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
336	>	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);
337	>
338	>	return 1.0;
339	>	}
340	>
341	>
342	>	Double_t JetTools::Beta(const PFJet *jet, const Vertex *vertex, Double_t  delta_z){
343	>	double Pt_jets= 0. ;
344	>	double Pt_jetsTot = 0. ;
345	>
346	>	for(UInt_t i=0;i<jet->NPFCands();i++){
347	>	if(jet->PFCand(i)->TrackerTrk()){
348	>	Pt_jetsTot += jet->PFCand(i)->TrackerTrk()->Pt();
349	>	double pDz = TMath::Abs(jet->PFCand(i)->TrackerTrk()->DzCorrected(*vertex));
350	>	if(pDz < delta_z){
351	>	Pt_jets += jet->PFCand(i)->TrackerTrk()->Pt();
352	>	}
353	>	}
354	>	}
355	>
356	>	Double_t beta = 0.;
357	>	if (Pt_jetsTot > 0)
358	>	beta = Pt_jets/Pt_jetsTot;
359	>
360	>	return beta;
361	>	}
362	>
363	>	Double_t JetTools::Beta2(const PFJet *jet, const Vertex *vertex, Double_t  delta_z){
364	>	double Pt_jets= 0. ;
365	>	double Pt_jetsTot = 0. ;
366	>
367	>	for(UInt_t i=0;i<jet->NPFCands();i++){
368	>	if(jet->PFCand(i)->BestTrk()){
369	>	Pt_jetsTot += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
370	>	double pDz = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertex));
371	>	if(pDz < delta_z){
372	>	Pt_jets += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
373	>	}
374	>	}
375	>	}
376	>
377	>	Double_t beta = 1.0;
378	>	if (Pt_jetsTot > 0)
379	>	beta = Pt_jets/Pt_jetsTot;
380	>	return beta;
381	>	}
382	>
383	>
384	>	Bool_t  JetTools::PassBetaVertexAssociationCut(const PFJet *jet, const Vertex *referenceVertex, const VertexCol *vertices, Double_t delta_z) {
385	>
386	>	Bool_t passBetaCut = kTRUE;
387	>	if(vertices->GetEntries() > 0) {
388	>	Double_t Beta = JetTools::Beta(jet, referenceVertex, 0.2);
389	>	Double_t Beta_other = 0.0;
390	>	for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
391	>	if (referenceVertex == vertices->At(nv)) continue;
392	>	Double_t BetaAux = JetTools::Beta(jet, vertices->At(nv), 0.2);
393	>	if(BetaAux > Beta_other) Beta_other = BetaAux;
394	>	}
395	>	if(Beta_other > Beta) passBetaCut = kFALSE;
396	>	}
397	>
398	>	return passBetaCut;
399	>
400	>	}
401	>
402	>	Bool_t  JetTools::PassBeta2VertexAssociationCut(const PFJet *jet, const Vertex *referenceVertex, const VertexCol *vertices, Double_t delta_z) {
403	>
404	>	Bool_t passBetaCut = kTRUE;
405	>	if(vertices->GetEntries() > 0) {
406	>	Double_t Beta = JetTools::Beta2(jet, referenceVertex, 0.2);
407	>	Double_t Beta_other = 0.0;
408	>	for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
409	>	if (referenceVertex == vertices->At(nv)) continue;
410	>	Double_t BetaAux = JetTools::Beta2(jet, vertices->At(nv), 0.2);
411	>	if(BetaAux > Beta_other) Beta_other = BetaAux;
412	>	}
413	>	if(Beta_other > Beta) passBetaCut = kFALSE;
414	>	}
415	>
416	>	return passBetaCut;
417	>
418	>	}
419	>
420	>
421	>	Int_t JetTools::MaxBetaVertexIndex(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){
422	>
423	>	Int_t vertexIndex = -1;
424	>	double beta = -0.1;
425	>	for (UInt_t v=0; v < vertices->GetEntries(); v++){
426	>	Double_t betaTmp = JetTools::Beta(jet, vertices->At(v), delta_z);
427	>	if (betaTmp > beta) {
428	>	beta = betaTmp;
429	>	vertexIndex = v;
430	>	}
431	>	}
432	>	return vertexIndex;
433	>
434		}
435
436	+	Int_t JetTools::MaxBeta2VertexIndex(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){
437	+
438	+	Int_t vertexIndex = -1;
439	+	double beta = -0.1;
440	+	for (UInt_t v=0; v < vertices->GetEntries(); v++){
441	+	Double_t betaTmp = JetTools::Beta2(jet, vertices->At(v), delta_z);
442	+	if (betaTmp > beta) {
443	+	beta = betaTmp;
444	+	vertexIndex = v;
445	+	}
446	+	}
447	+	return vertexIndex;
448	+
449	+	}
450	+
451	+
452	+	Int_t JetTools::JetToPVAssociation(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){
453	+
454	+	std::vector<float> verticesPt2(vertices->GetEntries());
455	+	for(UInt_t i=0;i<jet->NPFCands();i++){
456	+	if(jet->PFCand(i)->BestTrk()){
457	+	double minDZ = delta_z;
458	+	int trackVertexIndex = -1;
459	+	for (UInt_t v=0; v < vertices->GetEntries(); v++){
460	+	if (minDZ > TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)))) {
461	+	minDZ = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)));
462	+	trackVertexIndex = v;
463	+	}
464	+	}
465	+	if (trackVertexIndex < 0) continue;
466	+	verticesPt2[trackVertexIndex]+= jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
467	+	}
468	+	}
469	+
470	+	Int_t vertexIndex = 0;
471	+	float pt2Max = 0;
472	+	for (uint i=0; i < verticesPt2.size(); ++i){
473	+	if (pt2Max < verticesPt2[i]) {
474	+	pt2Max = verticesPt2[i];
475	+	vertexIndex = i;
476	+	}
477	+	}
478	+	return vertexIndex;
479	+	}
480	+	const PFCandidate* JetTools::leadCand(const PFJet *iJet,int iPFType,bool i2nd) {
481	+	int lCount = 0;
482	+	const PFCandidate *lCand = 0;
483	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
484	+	lCand = iJet->PFCand(i0);
485	+	if(iPFType != -1 && lCand->PFType() != iPFType) continue;
486	+	if(lCount == 0 && !i2nd) break;
487	+	if(lCount >  0)          break;
488	+	lCount++;
489	+	}
490	+	return lCand;
491	+	}
492	+	Double_t JetTools::impactParameter(const PFJet *iJet,const Vertex *iVertex,bool iDZ) {
493	+	double lDZCorr = -1000;
494	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
495	+	const PFCandidate *pCand = iJet->PFCand(i0);
496	+	if(pCand->TrackerTrk() == 0) continue;
497	+	//if(pCand->Pt() < 1.) continue; => previous iterations
498	+	if(iDZ)  lDZCorr = fabs(pCand->TrackerTrk()->DzCorrected(*iVertex));
499	+	if(!iDZ) lDZCorr = fabs(pCand->TrackerTrk()->D0Corrected(*iVertex));
500	+	break;
501	+	}
502	+	return lDZCorr;
503	+	}
504	+	Double_t JetTools::dRMean(const PFJet *iJet,int iPFType) {
505	+	double lDRMean = 0;
506	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
507	+	const PFCandidate *pCand = iJet->PFCand(i0);
508	+	if(iPFType != -1 && pCand->PFType() != iPFType) continue;
509	+	double pDR = MathUtils::DeltaR(iJet->Mom(),pCand->Mom());
510	+	lDRMean    += pDR*(pCand->Pt())/iJet->RawMom().Pt();
511	+	}
512	+	return lDRMean;
513	+	}
514	+	Double_t JetTools::dR2Mean(const PFJet *iJet,int iPFType) {
515	+	double lDR2Mean = 0;
516	+	double lSumPt2 = 0;
517	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
518	+	const PFCandidate *pCand = iJet->PFCand(i0);
519	+	if(iPFType != -1 && pCand->PFType() != iPFType) continue;
520	+	lSumPt2   += pCand->Pt()*pCand->Pt();
521	+	double pDR = MathUtils::DeltaR(iJet->Mom(),pCand->Mom());
522	+	lDR2Mean    += pDR*pDR*(pCand->Pt()*pCand->Pt());
523	+	}
524	+	lDR2Mean/=lSumPt2;
525	+	return lDR2Mean;
526	+	}
527	+	Double_t JetTools::frac(const PFJet *iJet,Double_t iDRMax,Double_t iDRMin,Int_t iPFType) {
528	+	double lFrac = 0;
529	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
530	+	const PFCandidate *pCand = iJet->PFCand(i0);
531	+	if(iPFType != -1 && pCand->PFType() != iPFType) continue;
532	+	Double_t pDR = MathUtils::DeltaR(iJet->Mom(),pCand->Mom());
533	+	if(pDR > iDRMax) continue;
534	+	if(pDR < iDRMax-0.1) continue;
535	+	lFrac += pCand->Pt()/iJet->RawMom().Pt();
536	+	}
537	+	return lFrac;
538	+	}
539	+	Double_t JetTools::betaStar(const PFJet *iJet,const Vertex *iVertex,const VertexCol* iVertices,Double_t iDZCut) {
540	+	Double_t lTotal = 0;
541	+	Double_t lPileup = 0;
542	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
543	+	const PFCandidate* pPF   = iJet->PFCand(i0);
544	+	const Track* pTrack      = pPF->TrackerTrk();
545	+	//if(pPF->GsfTrk()) pTrack = pPF->GsfTrk(); ==> not used in CMSSW
546	+	if(pTrack == 0) continue;
547	+	lTotal += pTrack->Pt();
548	+	double pDZPV  = fabs(pTrack->DzCorrected(*iVertex));
549	+	double pDZMin = pDZPV;
550	+	for(unsigned int i1 = 0; i1 < iVertices->GetEntries(); i1++) {
551	+	const Vertex *pV = iVertices->At(i1);
552	+	if(pV->Ndof() < 4 ||
553	+	(pV->Position() - iVertex->Position()).R() < 0.02 ) continue;
554	+	pDZMin = TMath::Min(pDZMin,fabs(pTrack->DzCorrected(*pV)));
555	+	}
556	+	if(pDZPV > 0.2 && pDZMin < 0.2) lPileup += pTrack->Pt();
557	+	}
558	+	if(lTotal == 0) lTotal = 1;
559	+	return lPileup/(lTotal);
560	+	}
561	+	Double_t JetTools::betaStarClassic(const PFJet *iJet,const Vertex *iVertex,const VertexCol* iVertices) {
562	+	Double_t lTotal = 0;
563	+	Double_t lPileup = 0;
564	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
565	+	const PFCandidate* pPF   = iJet->PFCand(i0);
566	+	const Track* pTrack      = pPF->TrackerTrk();
567	+	//if(pPF->GsfTrk()) pTrack = pPF->GsfTrk(); ==> not used in CMSSW
568	+	if(pTrack == 0) continue;
569	+	lTotal += pTrack->Pt();
570	+	bool isPV    = iVertex->HasTrack(pPF->TrackerTrk());
571	+	bool isOtherV = false;
572	+	for(unsigned int i1 = 0; i1 < iVertices->GetEntries(); i1++) {
573	+	const Vertex *pV = iVertices->At(i1);
574	+	if(isOtherV || isPV) continue;
575	+	if(pV->Ndof() < 4 ||
576	+	(pV->Position() - iVertex->Position()).R() < 0.02 ) continue;
577	+	isOtherV    = pV->HasTrack(pPF->TrackerTrk());
578	+	}
579	+	if(!isPV && isOtherV) lPileup += pTrack->Pt();
580	+	}
581	+	if(lTotal == 0) lTotal = 1;
582	+	return lPileup/(lTotal);
583	+	}
584	+	Bool_t  JetTools::passPFLooseId(const PFJet *iJet) {
585	+	if(iJet->RawMom().E()                              == 0)       return false;
586	+	if(iJet->NeutralHadronEnergy()/iJet->RawMom().E()  >  0.99)    return false;
587	+	if(iJet->NeutralEmEnergy()/iJet->RawMom().E()      >  0.99)    return false;
588	+	if(iJet->NConstituents()                           <  2)       return false;
589	+	if(iJet->ChargedHadronEnergy()/iJet->RawMom().E()  <= 0     && fabs(iJet->Eta()) < 2.4 ) return false;
590	+	if(iJet->ChargedEmEnergy()/iJet->RawMom().E()      >  0.99  && fabs(iJet->Eta()) < 2.4 ) return false;
591	+	if(iJet->ChargedMultiplicity()                     < 1      && fabs(iJet->Eta()) < 2.4 ) return false;
592	+	//if(fabs(iJet->Eta())                               > 4.99) return false;
593	+	return true;
594	+	}
595	+	//Jet Width Variables
596	+	double JetTools::W(const PFJet *iJet,int iPFType,int iType) {
597	+	double lPtD    = 0;
598	+	double lSumPt  = 0;
599	+	double lSumPt2 = 0;
600	+	TMatrixDSym lCovMatrix(2); lCovMatrix = 0.;
601	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
602	+	const PFCandidate *pCand = iJet->PFCand(i0);
603	+	if(iPFType != -1 && pCand->PFType() != iPFType) continue;
604	+	double pDEta = iJet->Eta() - pCand->Eta();
605	+	double pDPhi = fabs(iJet->Phi()-pCand->Phi()); if(pDPhi > 2.*TMath::Pi() - pDPhi) pDPhi =  2.*TMath::Pi() - pDPhi;
606	+	lCovMatrix(0,0) += pCand->Pt()*pCand->Pt()*pDEta*pDEta;
607	+	lCovMatrix(0,1) += pCand->Pt()*pCand->Pt()*pDEta*pDPhi;
608	+	lCovMatrix(1,1) += pCand->Pt()*pCand->Pt()*pDPhi*pDPhi;
609	+	lPtD            += pCand->Pt()*pCand->Pt();
610	+	lSumPt          += pCand->Pt();
611	+	lSumPt2         += pCand->Pt()*pCand->Pt();
612	+	}
613	+	lCovMatrix(0,0) /= lSumPt2;
614	+	lCovMatrix(0,1) /= lSumPt2;
615	+	lCovMatrix(1,1) /= lSumPt2;
616	+	lCovMatrix(1,0)  = lCovMatrix(0,1);
617	+	lPtD             = sqrt(lPtD);
618	+	lPtD            /= lSumPt;
619	+	double lEtaW     = sqrt(lCovMatrix(0,0));
620	+	double lPhiW     = sqrt(lCovMatrix(1,1));
621	+	double lJetW     = 0.5*(lEtaW+lPhiW);
622	+	TVectorD lEigVals(2); lEigVals = TMatrixDSymEigen(lCovMatrix).GetEigenValues();
623	+	double lMajW     = sqrt(fabs(lEigVals(0)));
624	+	double lMinW     = sqrt(fabs(lEigVals(1)));
625	+
626	+	if(iType == 1) return lMajW;
627	+	if(iType == 2) return lMinW;
628	+	if(iType == 3) return lEtaW;
629	+	if(iType == 4) return lPhiW;
630	+	if(iType == 5) return lJetW;
631	+	return lPtD; //ptRMS
632	+	}
633	+	/*
634	+	double JetTools::genFrac(const PFJet *iJet) {
635	+	double lTrueFrac = 0;
636	+	for(UInt_t i0 = 0; i0 < fParticles->GetEntries(); i0++) {
637	+	const MCParticle *p = fParticles->At(i0);
638	+	if(p->Status() != 1) continue;
639	+	double pDEta = iJet->Eta() - p->Eta();
640	+	double pDPhi = fabs(iJet->Phi()-p->Phi()); if(pDPhi > 2.*TMath::Pi() - pDPhi) pDPhi =  2.*TMath::Pi() - pDPhi;
641	+	double pDR   = sqrt(pDEta*pDEta + pDPhi*pDPhi);
642	+	if(pDR > 0.5) continue;
643	+	lTrueFrac += p->Pt();
644	+	}
645	+	lTrueFrac/=iJet->Pt();
646	+	return lTrueFrac;
647	+	}
648	+	*/
649	+	/*
650	+	double* JetTools::subStructure(const PFJet *iJet) {
651	+	for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
652	+	const PFCandidate *pCand = iJet->PFCand(i0);
653	+	//Fast Jet
654	+	FJparticles.push_back( fastjet::PseudoJet( pCand->Px(),pCand->Py(),pCand->Pz(),pCand->Energy() ) );
655	+	}
656	+	//ReCluster
657	+	fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, 0.5);
658	+	int activeAreaRepeats = 1;
659	+	double ghostArea = 0.01;
660	+	double ghostEtaMax = 5.0;
661	+	fastjet::ActiveAreaSpec fjActiveArea(ghostEtaMax,activeAreaRepeats,ghostArea);
662	+	fjActiveArea.set_fj2_placement(true);
663	+
664	+	fastjet::AreaDefinition fjAreaDefinition(fastjet::active_area_explicit_ghosts, fjActiveArea );
665	+	fastjet::ClusterSequenceArea thisClustering(FJparticles, jetDef, fjAreaDefinition);
666	+	std::vector<fastjet::PseudoJet> out_jets = sorted_by_pt(thisClustering.inclusive_jets(0.0));
667	+	//std::cout << "===>  Size " << FJparticles.size() << " -- " << out_jets.size() << std::endl;
668	+	// Adding substructure
669	+	// define groomers
670	+	fastjet::Filter trimmer( fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.03)));
671	+	fastjet::Filter filter( fastjet::Filter(fastjet::JetDefinition(fastjet::cambridge_algorithm, 0.3), fastjet::SelectorNHardest(3)));
672	+	fastjet::Pruner pruner(fastjet::cambridge_algorithm, 0.1, 0.5);
673	+
674	+	std::vector<fastjet::Transformer const *> transformers;
675	+	transformers.push_back(&trimmer);
676	+	transformers.push_back(&filter);
677	+	transformers.push_back(&pruner);
678	+
679	+	// define n-subjettiness
680	+	NsubParameters paraNsub = NsubParameters(1.0, 0.8);
681	+	Nsubjettiness routine(nsub_kt_axes, paraNsub);
682	+	for (unsigned i0 = 0; i0 < out_jets.size(); i0++) {
683	+	int i1 = -1;
684	+	if(out_jets.at(i0).pt()  < 1.) continue;
685	+	//std::cout << " ---> " << out_jets.at(i0).pt() << std::endl;
686	+	for ( std::vector<fastjet::Transformer const *>::const_iterator itransf = transformers.begin(), itransfEnd = transformers.end(); itransf != itransfEnd; ++itransf ) {
687	+	i1++;
688	+	fastjet::PseudoJet transformedJet = out_jets.at(i0);
689	+	transformedJet = (**itransf)(transformedJet);
690	+	TLorentzVector jetcorr(transformedJet.px() * jec,transformedJet.py() * jec,transformedJet.pz() * jec,transformedJet.e() * jec);
691	+	if(i1 == 0) internalId_.trimmass_   = jetcorr.M();
692	+	if(i1 == 0) internalId_.trimarea_   = transformedJet.area();
693	+	if(i1 == 1) internalId_.filtermass_ = jetcorr.M();
694	+	if(i1 == 1) internalId_.filterarea_ = transformedJet.area();
695	+	if(i1 == 2) internalId_.prunedmass_ = jetcorr.M();
696	+	if(i1 == 2) internalId_.prunedarea_ = transformedJet.area();
697	+	if (transformedJet.constituents().size() > 1 && i1 == 2 ) {
698	+	std::vector<fastjet::PseudoJet> subjets = transformedJet.associated_cluster_sequence()->exclusive_subjets(transformedJet,2);
699	+	internalId_.nsubjets_     = subjets.size();
700	+	internalId_.massdrop_     = subjets.at(0).m()/transformedJet.m();
701	+	internalId_.massdropcorr_ = subjets.at(0).m()/internalId_.prunedmass_;
702	+	}
703	+	}
704	+	internalId_.tau1_ = routine.getTau(1, out_jets.at(i0).constituents());
705	+	internalId_.tau2_ = routine.getTau(2, out_jets.at(i0).constituents());
706	+	internalId_.tau3_ = routine.getTau(3, out_jets.at(i0).constituents());
707	+	internalId_.tau4_ = routine.getTau(4, out_jets.at(i0).constituents());
708	+	internalId_.tau2tau1_ = internalId_.tau2_/internalId_.tau1_;
709	+	}
710	+	}
711	+	}
712	+	*/

Diff Legend

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