| 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 |
|
|
| 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++) { |
| 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->Pt(); |
| 535 |
> |
lFrac += pCand->Pt()/iJet->RawMom().Pt(); |
| 536 |
|
} |
| 537 |
|
return lFrac; |
| 538 |
|
} |
| 544 |
|
const Track* pTrack = pPF->TrackerTrk(); |
| 545 |
|
//if(pPF->GsfTrk()) pTrack = pPF->GsfTrk(); ==> not used in CMSSW |
| 546 |
|
if(pTrack == 0) continue; |
| 547 |
< |
lTotal += pPF->Pt(); |
| 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; |
| 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(); |
| 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->E() == 0) return false; |
| 586 |
< |
if(iJet->NeutralHadronEnergy()/iJet->E() > 0.99) return false; |
| 587 |
< |
if(iJet->NeutralEmEnergy()/iJet->E() > 0.99) return false; |
| 588 |
< |
if(iJet->NConstituents() < 2) return false; |
| 589 |
< |
if(iJet->ChargedHadronEnergy()/iJet->E() <= 0 && fabs(iJet->Eta()) < 2.4 ) return false; |
| 590 |
< |
if(iJet->ChargedEmEnergy()/iJet->E() > 0.99 && fabs(iJet->Eta()) < 2.4 ) return false; |
| 591 |
< |
if(iJet->ChargedMultiplicity() < 1 && fabs(iJet->Eta()) < 2.4 ) return false; |
| 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 |
< |
|
| 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; |
| 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 |
+ |
*/ |
