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.14 by ceballos, Wed Apr 20 10:06:15 2011 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 30 | 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 52 | 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 74 | 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 96 | 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 118 | 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 292 | Line 294 | Double_t JetTools::MtHiggs(const Particl
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 = MathUtils::DeltaPhi(dilepton->Phi(),
298 <                                                     met->Phi());
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  
# Line 340 | Line 342 | Double_t JetTools::Beta(const TrackCol *
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 <
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    }
355  
356 <  Double_t beta = 1.0;
356 >  Double_t beta = 0.;
357    if (Pt_jetsTot > 0)
358      beta = Pt_jets/Pt_jetsTot;
359  
# Line 475 | Line 477 | Int_t JetTools::JetToPVAssociation(const
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 + */

Diff Legend

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