Utils/src/RecoilTools.cc

#include "MitPhysics/Utils/interface/RecoilTools.h"
#include "MitPhysics/Utils/interface/JetTools.h"
#include <TFile.h>

ClassImp(mithep::RecoilTools)

using namespace mithep;

RecoilTools::RecoilTools(TString iJetLowPtMVAFile,TString iJetHighPtMVAFile,TString iCutFile) { 
  fJetIDMVA = new JetIDMVA();
  fJetIDMVA->Initialize( JetIDMVA::kMET,iJetLowPtMVAFile,iJetHighPtMVAFile,JetIDMVA::kBaseline,iCutFile);
}
//--------------------------------------------------------------------------------------------------
RecoilTools::~RecoilTools() { 
  delete fJetIDMVA;
}
//--------------------------------------------------------------------------------------------------
bool RecoilTools::filter(const PFJet *iJet,Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2) { 
  double pDEta1 = iJet->Eta() - iEta1;
  double pDPhi1 = fabs(iJet->Phi() - iPhi1); if(pDPhi1 > 2.*TMath::Pi()-pDPhi1) pDPhi1 = 2.*TMath::Pi()-pDPhi1;
  double pDR1   = sqrt(pDEta1*pDEta1 + pDPhi1*pDPhi1);
  if(pDR1 < 0.5) return false;
  double pDEta2 = iJet->Eta() - iEta2;
  double pDPhi2 = fabs(iJet->Phi() - iPhi2); if(pDPhi2 > 2.*TMath::Pi()-pDPhi2) pDPhi2 = 2.*TMath::Pi()-pDPhi2;
  double pDR2   = sqrt(pDEta2*pDEta2 + pDPhi2*pDPhi2);
  if(pDR2 < 0.5) return false;
  return true;
}

//--------------------------------------------------------------------------------------------------
Met RecoilTools::pfRecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
                          const PFCandidateCol *iCands) { 
  double lSumEt = 0;
  FourVectorM lVec(0,0,0,0);
  for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) { lVec -= iCands->At(i0)->Mom(); lSumEt += iCands->At(i0)->Et();}
  Met lPFMet(lVec.Px(),lVec.Py()); 
  lPFMet.SetSumEt(lSumEt);
  lPFMet.SetMex  (lPFMet.Mex()+iVisPt*cos(iVisPhi));  
  lPFMet.SetMey  (lPFMet.Mey()+iVisPt*sin(iVisPhi));
  lPFMet.SetSumEt(lPFMet.SumEt()-iVisSumEt);
  return lPFMet;
}
//--------------------------------------------------------------------------------------------------
Met RecoilTools::trackMet(const PFCandidateCol *iCands,const Vertex *iVertex,Double_t iDZCut) { 
  double trkMetx  = 0;
  double trkMety  = 0;
  double trkSumEt = 0; 
  for(UInt_t i0=0; i0<iCands->GetEntries(); i0++) {
    const PFCandidate *pfcand = iCands->At(i0);
    if( (pfcand->HasTrackerTrk() && (fabs(pfcand->TrackerTrk()->DzCorrected(*iVertex))< iDZCut)) ||
        (pfcand->HasGsfTrk()     && (fabs(pfcand->GsfTrk()->DzCorrected(*iVertex))    < iDZCut)) ) {
      trkMetx  -= pfcand->Px();
      trkMety  -= pfcand->Py();
      trkSumEt += pfcand->Pt();
    }
  }
  Met lMet(trkMetx,trkMety);
  lMet.SetSumEt(trkSumEt);
  return lMet;
}

//--------------------------------------------------------------------------------------------------
//Compute the recoil => here this requires the vector sum of the visible components
//VisPt    => Vector sum pT  of the visible non-recoiling components
//VisPhi   => Vector sum Phi of the visible non-recoiling components
//visSumEt => Vector sum Et  =f the visible non-recoiling components
Met RecoilTools::trackRecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
                             const PFCandidateCol *iCands,const Vertex *iVertex,double iDZCut) { 
  Met lTrkMet = trackMet(iCands,iVertex,iDZCut);
  lTrkMet.SetMex  (lTrkMet.Mex()+iVisPt*cos(iVisPhi));  
  lTrkMet.SetMey  (lTrkMet.Mey()+iVisPt*sin(iVisPhi));
  lTrkMet.SetSumEt(lTrkMet.SumEt()-iVisSumEt);
  return lTrkMet;
}
//--------------------------------------------------------------------------------------------------
void RecoilTools::addNeut(const PFJet *iJet,FourVectorM &iVec,Double_t &iSumEt, 
                          FactorizedJetCorrector *iJetCorrector,const PileupEnergyDensityCol *iPUEnergyDensity,
                          int iSign) { 
  FourVectorM lVec(0,0,0,0);
  double lPt = fJetIDMVA->correctedPt(iJet,iJetCorrector,iPUEnergyDensity);
  if(iJet->RawMom().Pt() < 10) lPt = TMath::Max(iJet->RawMom().Pt()-iJet->JetArea()*iPUEnergyDensity->At(0)->Rho(),0.);
  lPt *= (iJet->NeutralEmEnergy()/iJet->E() + iJet->NeutralHadronEnergy()/iJet->E());
  lVec.SetPt(lPt); lVec.SetEta(iJet->Eta()); lVec.SetPhi(iJet->Phi()); lVec.SetM(iJet->Mass());
  if(iSign > 0) iVec -= lVec;
  if(iSign < 0) iVec += lVec;
  iSumEt += lPt;
  //=== Above was a bug in the training
  //if(iSign > 0) iSumEt += lPt;
  //if(iSign < 0) iSumEt -= lPt;
}

//--------------------------------------------------------------------------------------------------
//Corrected Jets
void RecoilTools::addNeut(const PFJet *iJet,FourVectorM &iVec,Double_t &iSumEt,double iRho,int iSign) { 
  FourVectorM lVec(0,0,0,0);
  double lPt = iJet->Pt();
  if(iJet->RawMom().Pt() < 10) lPt = TMath::Max(iJet->RawMom().Pt()-iJet->JetArea()*iRho,0.);//to be fixed
  //if(iJet->RawMom().Pt() < 10) lPt = iJet->RawMom().Pt()*iJet->L1OffsetCorrectionScale();
  lPt *= (iJet->NeutralEmEnergy()/iJet->E() + iJet->NeutralHadronEnergy()/iJet->E());
  lVec.SetPt(lPt); lVec.SetEta(iJet->Eta()); lVec.SetPhi(iJet->Phi()); lVec.SetM(iJet->Mass());
  if(iSign > 0) iVec   -= lVec;
  if(iSign < 0) iVec   += lVec;
  iSumEt += lPt;
  //=== Above was a bug in the training
  //if(iSign > 0) iSumEt += lPt;
  //if(iSign < 0) iSumEt -= lPt;
}

//--------------------------------------------------------------------------------------------------
Met RecoilTools::NoPUMet( const PFJetCol       *iJets,FactorizedJetCorrector *iJetCorrector,
                          const PileupEnergyDensityCol *iPileupEnergyDensity,
                          const PFCandidateCol *iCands,const Vertex *iVertex,const VertexCol *iVertices,
                          Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,Double_t iDZCut) { 

  FourVectorM lVec        (0,0,0,0); double lSumEt          = 0; 
  for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) { 
    const PFCandidate *pPF = iCands->At(i0);
    const Track* pTrack = pPF->TrackerTrk();
    if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
    if(pTrack        ==  0                           ) continue;
     lSumEt   += pPF->Pt();  //=> another bug
     if(         !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
           (pPF->HasGsfTrk()     && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex))    <iDZCut)))) continue; 
    lVec     -= pPF->Mom();
  }
  int lNPass = 0;
  for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
    const PFJet *pJet = iJets->At(i0);
    if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2))                                       continue; //Quick cleaning==> if not done already
    if(!fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity)) continue;
    addNeut(pJet,lVec,lSumEt,iJetCorrector,iPileupEnergyDensity);
    lNPass++;
  }
  Met lMet(lVec.Px(),lVec.Py());
  lMet.SetSumEt( lSumEt);
  return lMet;
}
//--------------------------------------------------------------------------------------------------
//Corrected Jets
Met RecoilTools::NoPUMet( const PFJetCol       *iJets,
                          const PFCandidateCol *iCands,const Vertex *iVertex,const VertexCol *iVertices,Double_t iRho,
                          Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,Double_t iDZCut) { 

  FourVectorM lVec        (0,0,0,0); double lSumEt          = 0; 
  for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) { 
    const PFCandidate *pPF = iCands->At(i0);
    const Track* pTrack = pPF->TrackerTrk();
    if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
    if(pTrack        ==  0                           ) continue;
    lSumEt   += pPF->Pt();  //=> another bug
    if(  !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
           (pPF->HasGsfTrk()     && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex))    <iDZCut)))) continue; 
    lVec     -= pPF->Mom();
  }
  for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
    const PFJet *pJet = iJets->At(i0);
    if(!fJetIDMVA->pass(pJet,iVertex,iVertices))                          continue;
    if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2))                             continue; //Quick cleaning==> if not done already
    addNeut(pJet,lVec,lSumEt,iRho);
  }
  Met lMet(lVec.Px(),lVec.Py());
  lMet.SetSumEt(lSumEt);
  return lMet;
}
//--------------------------------------------------------------------------------------------------
Met RecoilTools::NoPURecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
                            const PFJetCol            *iJets,FactorizedJetCorrector *iJetCorrector,
                            const PileupEnergyDensityCol *iPileupEnergyDensity,
                            const PFCandidateCol   *iCands   ,const Vertex *iVertex,
                            const VertexCol *iVertices,
                            Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
                            Double_t iDZCut) { 
  
  Met lNoPUMet = NoPUMet(iJets,iJetCorrector,iPileupEnergyDensity,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
  lNoPUMet.SetMex  (lNoPUMet.Mex()+iVisPt*cos(iVisPhi));  
  lNoPUMet.SetMey  (lNoPUMet.Mey()+iVisPt*sin(iVisPhi));
  lNoPUMet.SetSumEt(lNoPUMet.SumEt()-iVisSumEt);
  return lNoPUMet;
}
//--------------------------------------------------------------------------------------------------
//Corrected Jets
Met RecoilTools::NoPURecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
                            const PFJetCol       *iJets,const PFCandidateCol *iCands,
                            const Vertex *iVertex,const VertexCol            *iVertices,Double_t iRho,
                            Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
                            Double_t iDZCut) { 
  
  Met lNoPUMet = NoPUMet(iJets,iCands,iVertex,iVertices,iRho,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
  lNoPUMet.SetMex  (lNoPUMet.Mex()+iVisPt*cos(iVisPhi));  
  lNoPUMet.SetMey  (lNoPUMet.Mey()+iVisPt*sin(iVisPhi));
  lNoPUMet.SetSumEt(lNoPUMet.SumEt()-iVisSumEt);
  return lNoPUMet;
}
//--------------------------------------------------------------------------------------------------
Met RecoilTools::PUCMet( const PFJetCol       *iJets,FactorizedJetCorrector *iJetCorrector,
                         const PileupEnergyDensityCol *iPileupEnergyDensity,
                         const PFCandidateCol *iCands,
                         const Vertex *iVertex,const VertexCol *iVertices,
                         Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
                         Double_t iDZCut) { 

  FourVectorM lVec        (0,0,0,0); double lSumEt          = 0; 
  for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) { 
    const PFCandidate *pPF = iCands->At(i0);
    const Track* pTrack = pPF->TrackerTrk();
    if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
    if(pTrack == 0                                   && 
       (pPF->PFType() == PFCandidate::eGamma         || 
        pPF->PFType() == PFCandidate::eEGammaHF      || 
        pPF->PFType() == PFCandidate::eNeutralHadron || 
        pPF->PFType() == PFCandidate::eHadronHF      ))
      {lVec -= pPF->Mom(); lSumEt += pPF->Pt();}
    if(pTrack        ==  0                           ) continue;
    if(  !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
           (pPF->HasGsfTrk()     && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex))    <iDZCut)))) continue; 
    lVec     -= pPF->Mom();
    lSumEt   += pPF->Pt();
  }
  for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
    const PFJet *pJet = iJets->At(i0);
    if(!JetTools::passPFLooseId(pJet))                                             continue;
    if(fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity)) continue;
    if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2))                                      continue; //Quick cleaning==> if not done already
    addNeut(pJet,lVec,lSumEt,iJetCorrector,iPileupEnergyDensity,-1);
  }
  Met lMet(lVec.Px(),lVec.Py());
  lMet.SetSumEt(lSumEt);
  return lMet;
}
//--------------------------------------------------------------------------------------------------
//Corrected jets
Met RecoilTools::PUCMet( const PFJetCol       *iJets,const PFCandidateCol *iCands,
                         const Vertex *iVertex,const VertexCol *iVertices,Double_t iRho,
                         Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
                         Double_t iDZCut) { 

  FourVectorM lVec        (0,0,0,0); double lSumEt          = 0; 
  for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) { 
    const PFCandidate *pPF = iCands->At(i0);
    const Track* pTrack = pPF->TrackerTrk();
    if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
    if(pTrack == 0                                   && 
       (pPF->PFType() == PFCandidate::eGamma         || 
        pPF->PFType() == PFCandidate::eEGammaHF      || 
        pPF->PFType() == PFCandidate::eNeutralHadron || 
        pPF->PFType() == PFCandidate::eHadronHF      ))
      {lVec -= pPF->Mom(); lSumEt += pPF->Pt();}
    if(pTrack        ==  0                           ) continue;
    if(  !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
           (pPF->HasGsfTrk()     && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex))    <iDZCut)))) continue; 
    lVec     -= pPF->Mom();
    lSumEt   += pPF->Pt();
  }
  for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
    const PFJet *pJet = iJets->At(i0);
    if(!JetTools::passPFLooseId(pJet))                                   continue;
    if(fJetIDMVA->pass(pJet,iVertex,iVertices))                          continue;
    if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2))                            continue; //Quick cleaning==> if not done already
    addNeut(pJet,lVec,lSumEt,iRho,-1);
  }
  Met lMet(lVec.Px(),lVec.Py());
  lMet.SetSumEt(lSumEt);
  return lMet;
}
//----> This MET is a bug need to fix it
//--------------------------------------------------------------------------------------------------
Met RecoilTools::PUCRecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
                           const PFJetCol       *iJets,FactorizedJetCorrector *iJetCorrector,
                           const PileupEnergyDensityCol *iPileupEnergyDensity,
                           const PFCandidateCol *iCands,const Vertex *iVertex,
                           const VertexCol *iVertices,
                           Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
                           Double_t iDZCut) { 
  Met lPUCMet = PUCMet(iJets,iJetCorrector,iPileupEnergyDensity,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
  lPUCMet.SetMex  (lPUCMet.Mex()+iVisPt*cos(iVisPhi));  
  lPUCMet.SetMey  (lPUCMet.Mey()+iVisPt*sin(iVisPhi));
  lPUCMet.SetSumEt(lPUCMet.SumEt()-iVisSumEt);
  return lPUCMet;
}
//--------------------------------------------------------------------------------------------------
//Corrected Jets
Met RecoilTools::PUCRecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
                           const PFJetCol       *iJets,
                           const PFCandidateCol *iCands,const Vertex *iVertex,
                           const VertexCol *iVertices,Double_t iRho,
                           Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
                           Double_t iDZCut) { 
  Met lPUCMet = PUCMet(iJets,iCands,iVertex,iVertices,iRho,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
  lPUCMet.SetMex  (lPUCMet.Mex()+iVisPt*cos(iVisPhi));  
  lPUCMet.SetMey  (lPUCMet.Mey()+iVisPt*sin(iVisPhi));
  lPUCMet.SetSumEt(lPUCMet.SumEt()-iVisSumEt);
  return lPUCMet;
}
//--------------------------------------------------------------------------------------------------
Met RecoilTools::PUMet( const PFJetCol       *iJets,FactorizedJetCorrector *iJetCorrector,
                        const PileupEnergyDensityCol *iPileupEnergyDensity,
                        const PFCandidateCol *iCands,const Vertex *iVertex,
                        const VertexCol      *iVertices,
                        Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
                        Double_t iDZCut) { 

  FourVectorM lVec        (0,0,0,0); double lSumEt          = 0; 
  for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) { 
    const PFCandidate *pPF = iCands->At(i0);
    const Track* pTrack = pPF->TrackerTrk();
    if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
    if(pTrack        ==  0                           ) continue;
    if(   ((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
           (pPF->HasGsfTrk()     && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex))    <iDZCut)))) continue; 
    lVec     -= pPF->Mom();
    lSumEt   += pPF->Pt();
  }
  for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
    const PFJet *pJet = iJets->At(i0);
    if(!JetTools::passPFLooseId(pJet))                                             continue;
    if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2))                                      continue; //Quick cleaning
    if(fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity)) continue;
    addNeut(pJet,lVec,lSumEt,iJetCorrector,iPileupEnergyDensity);
  }
  Met lMet(lVec.Px(),lVec.Py());
  lMet.SetSumEt(lSumEt);
  return lMet;
}
//--------------------------------------------------------------------------------------------------
//Corrected Jets
Met RecoilTools::PUMet( const PFJetCol       *iJets,
                        const PFCandidateCol *iCands,const Vertex *iVertex,const VertexCol *iVertices,Double_t iRho,
                        Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
                        Double_t iDZCut) { 

  FourVectorM lVec        (0,0,0,0); double lSumEt          = 0; 
  for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) { 
    const PFCandidate *pPF = iCands->At(i0);
    const Track* pTrack = pPF->TrackerTrk();
    if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
    if(pTrack        ==  0                           ) continue;
    if(   ((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) ||
           (pPF->HasGsfTrk()     && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex))    <iDZCut)))) continue; 
    lVec     -= pPF->Mom();
    lSumEt   += pPF->Pt();
  }
  for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
    const PFJet *pJet = iJets->At(i0);
    if(!JetTools::passPFLooseId(pJet))                                   continue;
    if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2))                            continue; //Quick cleaning
    if(fJetIDMVA->pass(pJet,iVertex,iVertices))                          continue;
    addNeut(pJet,lVec,lSumEt,iRho);
  }
  Met lMet(lVec.Px(),lVec.Py());
  lMet.SetSumEt(lSumEt);
  return lMet;
}
//--------------------------------------------------------------------------------------------------
Met RecoilTools::PURecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
                          const PFJetCol       *iJets,FactorizedJetCorrector *iJetCorrector,
                          const PileupEnergyDensityCol *iPileupEnergyDensity,
                          const PFCandidateCol *iCands,const Vertex *iVertex,
                          const VertexCol      *iVertices,
                          Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
                          Double_t iDZCut) { 
  Met lPUMet = PUMet(iJets,iJetCorrector,iPileupEnergyDensity,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
  lPUMet.SetMex  (lPUMet.Mex()+iVisPt*cos(iVisPhi));  
  lPUMet.SetMey  (lPUMet.Mey()+iVisPt*sin(iVisPhi));
  lPUMet.SetSumEt(lPUMet.SumEt()-iVisSumEt);
  return lPUMet;
}
//--------------------------------------------------------------------------------------------------
//Corrected Jets
Met RecoilTools::PURecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
                          const PFJetCol       *iJets,
                          const PFCandidateCol *iCands,const Vertex *iVertex,const VertexCol *iVertices,Double_t iRho,
                          Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
                          Double_t iDZCut) { 
  Met lPUMet = PUMet(iJets,iCands,iVertex,iVertices,iRho,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
  lPUMet.SetMex  (lPUMet.Mex()+iVisPt*cos(iVisPhi));  
  lPUMet.SetMey  (lPUMet.Mey()+iVisPt*sin(iVisPhi));
  lPUMet.SetSumEt(lPUMet.SumEt()-iVisSumEt);
  return lPUMet;
}
Revision:	1.11
Committed:	Tue May 1 17:19:56 2012 UTC (13 years ago) by pharris
Content type:	text/plain
Branch:	MAIN
Changes since 1.10:	+1 -1 lines
Log Message:	changed neutral fraction
#	User	Rev	Content
1	pharris	1.1	#include "MitPhysics/Utils/interface/RecoilTools.h"
2			#include "MitPhysics/Utils/interface/JetTools.h"
3			#include <TFile.h>
4
5			ClassImp(mithep::RecoilTools)
6
7			using namespace mithep;
8
9	pharris	1.3	RecoilTools::RecoilTools(TString iJetLowPtMVAFile,TString iJetHighPtMVAFile,TString iCutFile) {
10	pharris	1.1	fJetIDMVA = new JetIDMVA();
11	pharris	1.6	fJetIDMVA->Initialize( JetIDMVA::kMET,iJetLowPtMVAFile,iJetHighPtMVAFile,JetIDMVA::kBaseline,iCutFile);
12	pharris	1.1	}
13			//--------------------------------------------------------------------------------------------------
14			RecoilTools::~RecoilTools() {
15			delete fJetIDMVA;
16			}
17			//--------------------------------------------------------------------------------------------------
18			bool RecoilTools::filter(const PFJet *iJet,Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2) {
19			double pDEta1 = iJet->Eta() - iEta1;
20			double pDPhi1 = fabs(iJet->Phi() - iPhi1); if(pDPhi1 > 2.TMath::Pi()-pDPhi1) pDPhi1 = 2.TMath::Pi()-pDPhi1;
21			double pDR1 = sqrt(pDEta1pDEta1 + pDPhi1pDPhi1);
22			if(pDR1 < 0.5) return false;
23			double pDEta2 = iJet->Eta() - iEta2;
24			double pDPhi2 = fabs(iJet->Phi() - iPhi2); if(pDPhi2 > 2.TMath::Pi()-pDPhi2) pDPhi2 = 2.TMath::Pi()-pDPhi2;
25			double pDR2 = sqrt(pDEta2pDEta2 + pDPhi2pDPhi2);
26			if(pDR2 < 0.5) return false;
27			return true;
28			}
29
30			//--------------------------------------------------------------------------------------------------
31			Met RecoilTools::pfRecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
32	pharris	1.9	const PFCandidateCol *iCands) {
33			double lSumEt = 0;
34			FourVectorM lVec(0,0,0,0);
35	pharris	1.10	for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) { lVec -= iCands->At(i0)->Mom(); lSumEt += iCands->At(i0)->Et();}
36	pharris	1.9	Met lPFMet(lVec.Px(),lVec.Py());
37			lPFMet.SetSumEt(lSumEt);
38	pharris	1.1	lPFMet.SetMex (lPFMet.Mex()+iVisPt*cos(iVisPhi));
39			lPFMet.SetMey (lPFMet.Mey()+iVisPt*sin(iVisPhi));
40			lPFMet.SetSumEt(lPFMet.SumEt()-iVisSumEt);
41			return lPFMet;
42			}
43			//--------------------------------------------------------------------------------------------------
44			Met RecoilTools::trackMet(const PFCandidateCol iCands,const Vertex iVertex,Double_t iDZCut) {
45			double trkMetx = 0;
46			double trkMety = 0;
47			double trkSumEt = 0;
48	pharris	1.9	for(UInt_t i0=0; i0<iCands->GetEntries(); i0++) {
49			const PFCandidate *pfcand = iCands->At(i0);
50	pharris	1.8	if( (pfcand->HasTrackerTrk() && (fabs(pfcand->TrackerTrk()->DzCorrected(*iVertex))< iDZCut)) \|\|
51			(pfcand->HasGsfTrk() && (fabs(pfcand->GsfTrk()->DzCorrected(*iVertex)) < iDZCut)) ) {
52	pharris	1.1	trkMetx -= pfcand->Px();
53			trkMety -= pfcand->Py();
54			trkSumEt += pfcand->Pt();
55			}
56			}
57			Met lMet(trkMetx,trkMety);
58			lMet.SetSumEt(trkSumEt);
59			return lMet;
60			}
61
62			//--------------------------------------------------------------------------------------------------
63			//Compute the recoil => here this requires the vector sum of the visible components
64			//VisPt => Vector sum pT of the visible non-recoiling components
65			//VisPhi => Vector sum Phi of the visible non-recoiling components
66			//visSumEt => Vector sum Et =f the visible non-recoiling components
67			Met RecoilTools::trackRecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
68			const PFCandidateCol iCands,const Vertex iVertex,double iDZCut) {
69			Met lTrkMet = trackMet(iCands,iVertex,iDZCut);
70			lTrkMet.SetMex (lTrkMet.Mex()+iVisPt*cos(iVisPhi));
71			lTrkMet.SetMey (lTrkMet.Mey()+iVisPt*sin(iVisPhi));
72			lTrkMet.SetSumEt(lTrkMet.SumEt()-iVisSumEt);
73			return lTrkMet;
74			}
75			//--------------------------------------------------------------------------------------------------
76			void RecoilTools::addNeut(const PFJet *iJet,FourVectorM &iVec,Double_t &iSumEt,
77			FactorizedJetCorrector iJetCorrector,const PileupEnergyDensityCol iPUEnergyDensity,
78			int iSign) {
79			FourVectorM lVec(0,0,0,0);
80			double lPt = fJetIDMVA->correctedPt(iJet,iJetCorrector,iPUEnergyDensity);
81	pharris	1.6	if(iJet->RawMom().Pt() < 10) lPt = TMath::Max(iJet->RawMom().Pt()-iJet->JetArea()*iPUEnergyDensity->At(0)->Rho(),0.);
82	pharris	1.1	lPt *= (iJet->NeutralEmEnergy()/iJet->E() + iJet->NeutralHadronEnergy()/iJet->E());
83			lVec.SetPt(lPt); lVec.SetEta(iJet->Eta()); lVec.SetPhi(iJet->Phi()); lVec.SetM(iJet->Mass());
84			if(iSign > 0) iVec -= lVec;
85			if(iSign < 0) iVec += lVec;
86			iSumEt += lPt;
87	pharris	1.4	//=== Above was a bug in the training
88			//if(iSign > 0) iSumEt += lPt;
89			//if(iSign < 0) iSumEt -= lPt;
90	pharris	1.1	}
91	pharris	1.2
92			//--------------------------------------------------------------------------------------------------
93			//Corrected Jets
94	pharris	1.9	void RecoilTools::addNeut(const PFJet *iJet,FourVectorM &iVec,Double_t &iSumEt,double iRho,int iSign) {
95	pharris	1.2	FourVectorM lVec(0,0,0,0);
96			double lPt = iJet->Pt();
97	pharris	1.9	if(iJet->RawMom().Pt() < 10) lPt = TMath::Max(iJet->RawMom().Pt()-iJet->JetArea()*iRho,0.);//to be fixed
98			//if(iJet->RawMom().Pt() < 10) lPt = iJet->RawMom().Pt()*iJet->L1OffsetCorrectionScale();
99	pharris	1.11	lPt *= (iJet->NeutralEmEnergy()/iJet->E() + iJet->NeutralHadronEnergy()/iJet->E());
100	pharris	1.2	lVec.SetPt(lPt); lVec.SetEta(iJet->Eta()); lVec.SetPhi(iJet->Phi()); lVec.SetM(iJet->Mass());
101	pharris	1.4	if(iSign > 0) iVec -= lVec;
102			if(iSign < 0) iVec += lVec;
103	pharris	1.2	iSumEt += lPt;
104	pharris	1.4	//=== Above was a bug in the training
105			//if(iSign > 0) iSumEt += lPt;
106			//if(iSign < 0) iSumEt -= lPt;
107	pharris	1.2	}
108
109	pharris	1.1	//--------------------------------------------------------------------------------------------------
110			Met RecoilTools::NoPUMet( const PFJetCol iJets,FactorizedJetCorrector iJetCorrector,
111			const PileupEnergyDensityCol *iPileupEnergyDensity,
112	pharris	1.3	const PFCandidateCol iCands,const Vertex iVertex,const VertexCol *iVertices,
113	pharris	1.1	Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,Double_t iDZCut) {
114
115			FourVectorM lVec (0,0,0,0); double lSumEt = 0;
116			for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
117			const PFCandidate *pPF = iCands->At(i0);
118			const Track* pTrack = pPF->TrackerTrk();
119			if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
120			if(pTrack == 0 ) continue;
121	pharris	1.9	lSumEt += pPF->Pt(); //=> another bug
122			if( !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) \|\|
123	pharris	1.1	(pPF->HasGsfTrk() && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex)) <iDZCut)))) continue;
124			lVec -= pPF->Mom();
125			}
126	pharris	1.4	int lNPass = 0;
127	pharris	1.1	for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
128			const PFJet *pJet = iJets->At(i0);
129	pharris	1.4	if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2)) continue; //Quick cleaning==> if not done already
130	pharris	1.3	if(!fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity)) continue;
131	pharris	1.1	addNeut(pJet,lVec,lSumEt,iJetCorrector,iPileupEnergyDensity);
132	pharris	1.4	lNPass++;
133	pharris	1.1	}
134			Met lMet(lVec.Px(),lVec.Py());
135	pharris	1.4	lMet.SetSumEt( lSumEt);
136	pharris	1.1	return lMet;
137			}
138			//--------------------------------------------------------------------------------------------------
139	pharris	1.2	//Corrected Jets
140	pharris	1.3	Met RecoilTools::NoPUMet( const PFJetCol *iJets,
141	pharris	1.9	const PFCandidateCol iCands,const Vertex iVertex,const VertexCol *iVertices,Double_t iRho,
142	pharris	1.2	Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,Double_t iDZCut) {
143
144			FourVectorM lVec (0,0,0,0); double lSumEt = 0;
145			for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
146			const PFCandidate *pPF = iCands->At(i0);
147			const Track* pTrack = pPF->TrackerTrk();
148			if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
149			if(pTrack == 0 ) continue;
150	pharris	1.9	lSumEt += pPF->Pt(); //=> another bug
151	pharris	1.2	if( !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) \|\|
152			(pPF->HasGsfTrk() && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex)) <iDZCut)))) continue;
153			lVec -= pPF->Mom();
154			}
155			for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
156			const PFJet *pJet = iJets->At(i0);
157	pharris	1.3	if(!fJetIDMVA->pass(pJet,iVertex,iVertices)) continue;
158	pharris	1.2	if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2)) continue; //Quick cleaning==> if not done already
159	pharris	1.9	addNeut(pJet,lVec,lSumEt,iRho);
160	pharris	1.2	}
161			Met lMet(lVec.Px(),lVec.Py());
162			lMet.SetSumEt(lSumEt);
163			return lMet;
164			}
165			//--------------------------------------------------------------------------------------------------
166	pharris	1.1	Met RecoilTools::NoPURecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
167	pharris	1.3	const PFJetCol iJets,FactorizedJetCorrector iJetCorrector,
168	pharris	1.1	const PileupEnergyDensityCol *iPileupEnergyDensity,
169	pharris	1.3	const PFCandidateCol iCands ,const Vertex iVertex,
170			const VertexCol *iVertices,
171	pharris	1.1	Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
172			Double_t iDZCut) {
173
174	pharris	1.3	Met lNoPUMet = NoPUMet(iJets,iJetCorrector,iPileupEnergyDensity,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
175	pharris	1.1	lNoPUMet.SetMex (lNoPUMet.Mex()+iVisPt*cos(iVisPhi));
176			lNoPUMet.SetMey (lNoPUMet.Mey()+iVisPt*sin(iVisPhi));
177			lNoPUMet.SetSumEt(lNoPUMet.SumEt()-iVisSumEt);
178			return lNoPUMet;
179			}
180			//--------------------------------------------------------------------------------------------------
181	pharris	1.2	//Corrected Jets
182			Met RecoilTools::NoPURecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
183	pharris	1.3	const PFJetCol iJets,const PFCandidateCol iCands,
184	pharris	1.9	const Vertex iVertex,const VertexCol iVertices,Double_t iRho,
185	pharris	1.2	Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
186			Double_t iDZCut) {
187
188	pharris	1.9	Met lNoPUMet = NoPUMet(iJets,iCands,iVertex,iVertices,iRho,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
189	pharris	1.2	lNoPUMet.SetMex (lNoPUMet.Mex()+iVisPt*cos(iVisPhi));
190			lNoPUMet.SetMey (lNoPUMet.Mey()+iVisPt*sin(iVisPhi));
191			lNoPUMet.SetSumEt(lNoPUMet.SumEt()-iVisSumEt);
192			return lNoPUMet;
193			}
194			//--------------------------------------------------------------------------------------------------
195	pharris	1.1	Met RecoilTools::PUCMet( const PFJetCol iJets,FactorizedJetCorrector iJetCorrector,
196			const PileupEnergyDensityCol *iPileupEnergyDensity,
197	pharris	1.3	const PFCandidateCol *iCands,
198			const Vertex iVertex,const VertexCol iVertices,
199	pharris	1.1	Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
200			Double_t iDZCut) {
201
202			FourVectorM lVec (0,0,0,0); double lSumEt = 0;
203			for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
204			const PFCandidate *pPF = iCands->At(i0);
205			const Track* pTrack = pPF->TrackerTrk();
206			if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
207			if(pTrack == 0 &&
208			(pPF->PFType() == PFCandidate::eGamma \|\|
209			pPF->PFType() == PFCandidate::eEGammaHF \|\|
210			pPF->PFType() == PFCandidate::eNeutralHadron \|\|
211			pPF->PFType() == PFCandidate::eHadronHF ))
212			{lVec -= pPF->Mom(); lSumEt += pPF->Pt();}
213			if(pTrack == 0 ) continue;
214			if( !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) \|\|
215			(pPF->HasGsfTrk() && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex)) <iDZCut)))) continue;
216			lVec -= pPF->Mom();
217			lSumEt += pPF->Pt();
218			}
219			for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
220			const PFJet *pJet = iJets->At(i0);
221	pharris	1.3	if(!JetTools::passPFLooseId(pJet)) continue;
222			if(fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity)) continue;
223			if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2)) continue; //Quick cleaning==> if not done already
224	pharris	1.1	addNeut(pJet,lVec,lSumEt,iJetCorrector,iPileupEnergyDensity,-1);
225			}
226			Met lMet(lVec.Px(),lVec.Py());
227			lMet.SetSumEt(lSumEt);
228			return lMet;
229			}
230	pharris	1.2	//--------------------------------------------------------------------------------------------------
231			//Corrected jets
232	pharris	1.3	Met RecoilTools::PUCMet( const PFJetCol iJets,const PFCandidateCol iCands,
233	pharris	1.9	const Vertex iVertex,const VertexCol iVertices,Double_t iRho,
234	pharris	1.2	Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
235			Double_t iDZCut) {
236
237			FourVectorM lVec (0,0,0,0); double lSumEt = 0;
238			for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
239			const PFCandidate *pPF = iCands->At(i0);
240			const Track* pTrack = pPF->TrackerTrk();
241			if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
242			if(pTrack == 0 &&
243			(pPF->PFType() == PFCandidate::eGamma \|\|
244			pPF->PFType() == PFCandidate::eEGammaHF \|\|
245			pPF->PFType() == PFCandidate::eNeutralHadron \|\|
246			pPF->PFType() == PFCandidate::eHadronHF ))
247			{lVec -= pPF->Mom(); lSumEt += pPF->Pt();}
248			if(pTrack == 0 ) continue;
249			if( !((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) \|\|
250			(pPF->HasGsfTrk() && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex)) <iDZCut)))) continue;
251			lVec -= pPF->Mom();
252			lSumEt += pPF->Pt();
253			}
254			for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
255			const PFJet *pJet = iJets->At(i0);
256			if(!JetTools::passPFLooseId(pJet)) continue;
257	pharris	1.3	if(fJetIDMVA->pass(pJet,iVertex,iVertices)) continue;
258	pharris	1.2	if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2)) continue; //Quick cleaning==> if not done already
259	pharris	1.9	addNeut(pJet,lVec,lSumEt,iRho,-1);
260	pharris	1.2	}
261			Met lMet(lVec.Px(),lVec.Py());
262			lMet.SetSumEt(lSumEt);
263			return lMet;
264			}
265	pharris	1.1	//----> This MET is a bug need to fix it
266			//--------------------------------------------------------------------------------------------------
267			Met RecoilTools::PUCRecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
268			const PFJetCol iJets,FactorizedJetCorrector iJetCorrector,
269			const PileupEnergyDensityCol *iPileupEnergyDensity,
270			const PFCandidateCol iCands,const Vertex iVertex,
271	pharris	1.3	const VertexCol *iVertices,
272	pharris	1.1	Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
273			Double_t iDZCut) {
274	pharris	1.3	Met lPUCMet = PUCMet(iJets,iJetCorrector,iPileupEnergyDensity,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
275	pharris	1.1	lPUCMet.SetMex (lPUCMet.Mex()+iVisPt*cos(iVisPhi));
276			lPUCMet.SetMey (lPUCMet.Mey()+iVisPt*sin(iVisPhi));
277			lPUCMet.SetSumEt(lPUCMet.SumEt()-iVisSumEt);
278			return lPUCMet;
279			}
280			//--------------------------------------------------------------------------------------------------
281	pharris	1.2	//Corrected Jets
282			Met RecoilTools::PUCRecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
283			const PFJetCol *iJets,
284			const PFCandidateCol iCands,const Vertex iVertex,
285	pharris	1.9	const VertexCol *iVertices,Double_t iRho,
286	pharris	1.2	Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
287			Double_t iDZCut) {
288	pharris	1.9	Met lPUCMet = PUCMet(iJets,iCands,iVertex,iVertices,iRho,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
289	pharris	1.2	lPUCMet.SetMex (lPUCMet.Mex()+iVisPt*cos(iVisPhi));
290			lPUCMet.SetMey (lPUCMet.Mey()+iVisPt*sin(iVisPhi));
291			lPUCMet.SetSumEt(lPUCMet.SumEt()-iVisSumEt);
292			return lPUCMet;
293			}
294			//--------------------------------------------------------------------------------------------------
295	pharris	1.1	Met RecoilTools::PUMet( const PFJetCol iJets,FactorizedJetCorrector iJetCorrector,
296			const PileupEnergyDensityCol *iPileupEnergyDensity,
297			const PFCandidateCol iCands,const Vertex iVertex,
298	pharris	1.3	const VertexCol *iVertices,
299	pharris	1.1	Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
300			Double_t iDZCut) {
301
302			FourVectorM lVec (0,0,0,0); double lSumEt = 0;
303			for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
304			const PFCandidate *pPF = iCands->At(i0);
305			const Track* pTrack = pPF->TrackerTrk();
306			if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
307			if(pTrack == 0 ) continue;
308			if( ((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) \|\|
309			(pPF->HasGsfTrk() && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex)) <iDZCut)))) continue;
310			lVec -= pPF->Mom();
311			lSumEt += pPF->Pt();
312			}
313			for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
314			const PFJet *pJet = iJets->At(i0);
315	pharris	1.3	if(!JetTools::passPFLooseId(pJet)) continue;
316			if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2)) continue; //Quick cleaning
317			if(fJetIDMVA->pass(pJet,iVertex,iVertices,iJetCorrector,iPileupEnergyDensity)) continue;
318	pharris	1.1	addNeut(pJet,lVec,lSumEt,iJetCorrector,iPileupEnergyDensity);
319			}
320			Met lMet(lVec.Px(),lVec.Py());
321			lMet.SetSumEt(lSumEt);
322			return lMet;
323			}
324			//--------------------------------------------------------------------------------------------------
325	pharris	1.2	//Corrected Jets
326			Met RecoilTools::PUMet( const PFJetCol *iJets,
327	pharris	1.9	const PFCandidateCol iCands,const Vertex iVertex,const VertexCol *iVertices,Double_t iRho,
328	pharris	1.2	Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
329			Double_t iDZCut) {
330
331			FourVectorM lVec (0,0,0,0); double lSumEt = 0;
332			for(UInt_t i0 = 0; i0 < iCands->GetEntries(); i0++) {
333			const PFCandidate *pPF = iCands->At(i0);
334			const Track* pTrack = pPF->TrackerTrk();
335			if(pPF->GsfTrk()) pTrack = pPF->GsfTrk();
336			if(pTrack == 0 ) continue;
337			if( ((pPF->HasTrackerTrk() && (fabs(pPF->TrackerTrk()->DzCorrected(*iVertex))<iDZCut)) \|\|
338			(pPF->HasGsfTrk() && (fabs(pPF->GsfTrk()->DzCorrected(*iVertex)) <iDZCut)))) continue;
339			lVec -= pPF->Mom();
340			lSumEt += pPF->Pt();
341			}
342			for(UInt_t i0 = 0; i0 < iJets->GetEntries(); i0++) {
343			const PFJet *pJet = iJets->At(i0);
344			if(!JetTools::passPFLooseId(pJet)) continue;
345			if(!filter(pJet,iPhi1,iEta1,iPhi2,iEta2)) continue; //Quick cleaning
346	pharris	1.3	if(fJetIDMVA->pass(pJet,iVertex,iVertices)) continue;
347	pharris	1.9	addNeut(pJet,lVec,lSumEt,iRho);
348	pharris	1.2	}
349			Met lMet(lVec.Px(),lVec.Py());
350			lMet.SetSumEt(lSumEt);
351			return lMet;
352			}
353			//--------------------------------------------------------------------------------------------------
354	pharris	1.1	Met RecoilTools::PURecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
355			const PFJetCol iJets,FactorizedJetCorrector iJetCorrector,
356			const PileupEnergyDensityCol *iPileupEnergyDensity,
357			const PFCandidateCol iCands,const Vertex iVertex,
358	pharris	1.3	const VertexCol *iVertices,
359	pharris	1.1	Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
360			Double_t iDZCut) {
361	pharris	1.3	Met lPUMet = PUMet(iJets,iJetCorrector,iPileupEnergyDensity,iCands,iVertex,iVertices,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
362	pharris	1.1	lPUMet.SetMex (lPUMet.Mex()+iVisPt*cos(iVisPhi));
363			lPUMet.SetMey (lPUMet.Mey()+iVisPt*sin(iVisPhi));
364			lPUMet.SetSumEt(lPUMet.SumEt()-iVisSumEt);
365			return lPUMet;
366			}
367	pharris	1.2	//--------------------------------------------------------------------------------------------------
368			//Corrected Jets
369			Met RecoilTools::PURecoil(Double_t iVisPt,Double_t iVisPhi,Double_t iVisSumEt,
370			const PFJetCol *iJets,
371	pharris	1.9	const PFCandidateCol iCands,const Vertex iVertex,const VertexCol *iVertices,Double_t iRho,
372	pharris	1.2	Double_t iPhi1,Double_t iEta1,Double_t iPhi2,Double_t iEta2,
373			Double_t iDZCut) {
374	pharris	1.9	Met lPUMet = PUMet(iJets,iCands,iVertex,iVertices,iRho,iPhi1,iEta1,iPhi2,iEta2,iDZCut);
375	pharris	1.2	lPUMet.SetMex (lPUMet.Mex()+iVisPt*cos(iVisPhi));
376			lPUMet.SetMey (lPUMet.Mey()+iVisPt*sin(iVisPhi));
377			lPUMet.SetSumEt(lPUMet.SumEt()-iVisSumEt);
378			return lPUMet;
379			}