| 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 |
|
|
| 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; |
| 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; |
| 592 |
> |
//if(fabs(iJet->Eta()) > 4.99) return false; |
| 593 |
|
return true; |
| 594 |
|
} |
| 595 |
< |
|
| 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; |
