ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitPhysics/Utils/src/JetTools.cc
(Generate patch)

Comparing UserCode/MitPhysics/Utils/src/JetTools.cc (file contents):
Revision 1.19 by pharris, Fri Apr 13 14:24:00 2012 UTC vs.
Revision 1.32 by pharris, Tue Sep 25 15:39:15 2012 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  
# Line 342 | Line 344 | Double_t JetTools::Beta(const PFJet *jet
344    double Pt_jetsTot = 0. ;
345    
346    for(UInt_t i=0;i<jet->NPFCands();i++){
347 <    if(jet->PFCand(i)->BestTrk()){
348 <      Pt_jetsTot += jet->PFCand(i)->BestTrk()->Pt();
349 <      double pDz = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertex));
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)->BestTrk()->Pt();
351 >        Pt_jets += jet->PFCand(i)->TrackerTrk()->Pt();
352        }
353      }
354    }
# Line 491 | Line 493 | Double_t JetTools::impactParameter(const
493    double lDZCorr = -1000;
494    for(UInt_t i0 = 0; i0 < iJet->NPFCands(); i0++) {
495      const PFCandidate *pCand = iJet->PFCand(i0);
496 <    if(pCand->BestTrk() == 0) continue;
496 >    if(pCand->TrackerTrk() == 0) continue;
497      //if(pCand->Pt() < 1.) continue; => previous iterations
498 <    if(iDZ)  lDZCorr = fabs(pCand->BestTrk()->DzCorrected(*iVertex));
499 <    if(!iDZ) lDZCorr = fabs(pCand->BestTrk()->D0Corrected(*iVertex));
498 >    if(iDZ)  lDZCorr = fabs(pCand->TrackerTrk()->DzCorrected(*iVertex));
499 >    if(!iDZ) lDZCorr = fabs(pCand->TrackerTrk()->D0Corrected(*iVertex));
500      break;
501    }
502    return lDZCorr;
# Line 509 | Line 511 | Double_t JetTools::dRMean(const PFJet *i
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++) {
# Line 517 | Line 532 | Double_t JetTools::frac(const PFJet *iJe
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   }
# Line 527 | Line 542 | Double_t JetTools::betaStar(const PFJet
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();
545 >    //if(pPF->GsfTrk()) pTrack = pPF->GsfTrk(); ==> not used in CMSSW
546      if(pTrack == 0) continue;
547 <
548 <    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(fabs(pTrack->DzCorrected(*pV))             < iDZCut  &&
553 <         (pV->Position() - iVertex->Position()).R() > 0.02 )
554 <        { lPileup  +=pPF->Pt();}
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->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;

Diff Legend

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