LeptonSelection/src/IsolationSelection.cc

#include  <math.h>

#include "IsolationSelection.h"
#include "IsolationSelectionDefs.h"

#include "MathUtils.h"
#include "MuonTools.h"
#include "MuonIDMVA.h"
#include "ElectronTools.h"
#include "ElectronIDMVA.h"

using namespace mithep;

mithep::MuonIDMVA     * muIsoMVA;
mithep::MuonTools       muT;
mithep::ElectronIDMVA * eleIsoMVA;
mithep::ElectronTools   eleT;

// global hack to sync
double gChargedIso;
double gGammaIso;
double gNeutralIso;

extern vector<bool> PFnoPUflag;

//--------------------------------------------------------------------------------------------------
Float_t computePFMuonIso(const mithep::Muon *muon, 
                         const mithep::Vertex * vtx,
                         const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                         const Double_t dRMax)
//--------------------------------------------------------------------------------------------------
{
  const Double_t dRMin    = 0;
  const Double_t neuPtMin = 1.0;
  const Double_t dzMax    = 0.1;
    
  Double_t zLepton = (muon->BestTrk()) ? muon->BestTrk()->DzCorrected(*vtx) : 0.0;
  
  Float_t iso=0;
  for(UInt_t ipf=0; ipf<fPFCandidates->GetEntries(); ipf++) {
    const PFCandidate *pfcand = fPFCandidates->At(ipf);
    
    if(!pfcand->HasTrk() && (pfcand->Pt()<=neuPtMin)) continue;  // pT cut on neutral particles
    
    // exclude THE muon
    if(pfcand->TrackerTrk() && muon->TrackerTrk() && (pfcand->TrackerTrk()==muon->TrackerTrk())) continue;
    
    // dz cut
    Double_t dz = (pfcand->BestTrk()) ? fabs(pfcand->BestTrk()->DzCorrected(*vtx) - zLepton) : 0;
    if(dz >= dzMax) continue;
    
    // check iso cone
    Double_t dr = MathUtils::DeltaR(muon->Mom(), pfcand->Mom());
    if(dr<dRMax && dr>=dRMin)
      iso += pfcand->Pt(); 
  }
  
  return iso;
}

//--------------------------------------------------------------------------------------------------
Float_t computePFEleIso(const mithep::Electron *electron, 
                        const mithep::Vertex * fVertex,
                        const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                        const Double_t dRMax)
//--------------------------------------------------------------------------------------------------
{
  const Double_t dRMin    = 0;
  const Double_t neuPtMin = 1.0;
  const Double_t dzMax    = 0.1;
    
  Double_t zLepton = (electron->BestTrk()) ? electron->BestTrk()->DzCorrected(*fVertex) : 0.0;
  
  Float_t iso=0;
  for(UInt_t ipf=0; ipf<fPFCandidates->GetEntries(); ipf++) {
    const PFCandidate *pfcand = (PFCandidate*)(fPFCandidates->At(ipf));
    
    if(!pfcand->HasTrk() && (pfcand->Pt()<=neuPtMin)) continue;  // pT cut on neutral particles
    
    // dz cut
    Double_t dz = (pfcand->BestTrk()) ? fabs(pfcand->BestTrk()->DzCorrected(*fVertex) - zLepton) : 0;
    if(dz >= dzMax) continue;
    
    // remove THE electron
    if(pfcand->TrackerTrk() && electron->TrackerTrk() && (pfcand->TrackerTrk()==electron->TrackerTrk())) continue;
    if(pfcand->GsfTrk()     && electron->GsfTrk()     && (pfcand->GsfTrk()==electron->GsfTrk()))         continue;
    
    // check iso cone
    Double_t dr = MathUtils::DeltaR(electron->Mom(), pfcand->Mom());
    if(dr<dRMax && dr>=dRMin) {
      // eta-strip veto for photons
      if((pfcand->PFType() == PFCandidate::eGamma) && fabs(electron->Eta() - pfcand->Eta()) < 0.025) continue;
      
      // Inner cone (one tower = dR < 0.07) veto for non-photon neutrals
      if(!pfcand->HasTrk() && (pfcand->PFType() == PFCandidate::eNeutralHadron) && 
         (MathUtils::DeltaR(electron->Mom(), pfcand->Mom()) < 0.07)) continue;
      
      iso += pfcand->Pt();
    }
  }
  
  return iso;
};

//--------------------------------------------------------------------------------------------------
bool pairwiseIsoSelection( ControlFlags &ctrl, 
                           vector<SimpleLepton> &lepvec, 
                           float rho ) 
//--------------------------------------------------------------------------------------------------
{

  bool passiso=true;

  for( int i=0; i<lepvec.size(); i++ ) 
    { 
      
      if( !(lepvec[i].is4l) ) continue;
      
      float effArea_ecal_i, effArea_hcal_i;
      if( lepvec[i].isEB ) { 
        if( lepvec[i].type == 11 ) {
          effArea_ecal_i = 0.101;
          effArea_hcal_i = 0.021;
        } else { 
          effArea_ecal_i = 0.074;
          effArea_hcal_i = 0.022;
        }
      } else {
        if( lepvec[i].type == 11 ) {
          effArea_ecal_i = 0.046;
          effArea_hcal_i = 0.040;
        } else {
          effArea_ecal_i = 0.045;
          effArea_hcal_i = 0.030;
        }
      }
      
      float isoEcal_corr_i = lepvec[i].isoEcal - (effArea_ecal_i*rho);
      float isoHcal_corr_i = lepvec[i].isoHcal - (effArea_hcal_i*rho);

      for( int j=i+1; j<lepvec.size(); j++ ) 
        { 

          if( !(lepvec[j].is4l) ) continue;

          float effArea_ecal_j, effArea_hcal_j;
          if( lepvec[j].isEB ) { 
            if( lepvec[j].type == 11 ) {
              effArea_ecal_j = 0.101;
              effArea_hcal_j = 0.021;
            } else { 
              effArea_ecal_j = 0.074;
              effArea_hcal_j = 0.022;
            }
          } else {
            if( lepvec[j].type == 11 ) {
              effArea_ecal_j = 0.046;
              effArea_hcal_j = 0.040;
            } else {
              effArea_ecal_j = 0.045;
              effArea_hcal_j = 0.030;
            }
          }

          float isoEcal_corr_j = lepvec[j].isoEcal - (effArea_ecal_j*rho);
          float isoHcal_corr_j = lepvec[j].isoHcal - (effArea_hcal_j*rho);
          float RIso_i = (lepvec[i].isoTrk+isoEcal_corr_i+isoHcal_corr_i)/lepvec[i].vec.Pt();
          float RIso_j = (lepvec[j].isoTrk+isoEcal_corr_j+isoHcal_corr_j)/lepvec[j].vec.Pt();       
          float comboIso = RIso_i + RIso_j; 
          
          if( comboIso > 0.35 ) { 
            if( ctrl.debug ) cout << "combo failing for indices: " << i << "," << j << endl;
            passiso = false;
            return passiso;
          }
        }
    }
  
  return passiso;
}

//--------------------------------------------------------------------------------------------------
SelectionStatus muonIsoSelection(ControlFlags &ctrl, 
                                 const mithep::Muon * mu, 
                                 const mithep::Vertex * vtx, 
                                 const mithep::Array<mithep::PFCandidate> * fPFCandidateCol   ) 
//--------------------------------------------------------------------------------------------------
{
  float reliso = computePFMuonIso(mu,vtx,fPFCandidateCol,0.3)/mu->Pt();
  bool isEB = (fabs(mu->Eta()) < 1.479 ? 1 : 0 );   
  bool failiso = false; 
  if( isEB && mu->Pt() > 20 && reliso > PFISO_MU_LOOSE_EB_HIGHPT ) {  
    failiso = true;
  }
  if( isEB && mu->Pt() < 20 && reliso > PFISO_MU_LOOSE_EB_LOWPT ) { 
    failiso = true;
  }
  if( !(isEB) && mu->Pt() > 20 && reliso > PFISO_MU_LOOSE_EE_HIGHPT ) { 
    failiso = true;
  }
  if( !(isEB) && mu->Pt() < 20 && reliso > PFISO_MU_LOOSE_EE_LOWPT ) { 
    failiso = true;
  }

  SelectionStatus status;
  if( !failiso ) status.setStatus(SelectionStatus::LOOSEISO);
  if( !failiso ) status.setStatus(SelectionStatus::TIGHTISO);
  return status;

};

//--------------------------------------------------------------------------------------------------
SelectionStatus electronIsoSelection(ControlFlags &ctrl, 
                                     const mithep::Electron * ele, 
                                     const mithep::Vertex *fVertex,
                                     const mithep::Array<mithep::PFCandidate> * fPFCandidates) 
//--------------------------------------------------------------------------------------------------
{ 

  bool failiso=false;

  float reliso = computePFEleIso(ele,fVertex,fPFCandidates,0.4)/ele->Pt();

  if( ele->IsEB() && ele->Pt() > 20 && reliso > PFISO_ELE_LOOSE_EB_HIGHPT ) {
    failiso = true;
  }
  if( ele->IsEB() && ele->Pt() < 20 && reliso > PFISO_ELE_LOOSE_EB_LOWPT ) {
    failiso = true;
  }
  if( !(ele->IsEB()) && ele->Pt() > 20 && reliso > PFISO_ELE_LOOSE_EE_HIGHPT ) {
    failiso = true;
  }
  if( !(ele->IsEB()) && ele->Pt() < 20 && reliso > PFISO_ELE_LOOSE_EE_LOWPT ) {
    failiso = true;
  }

  SelectionStatus status;
  if( !failiso ) { 
    status.orStatus(SelectionStatus::LOOSEISO);
    status.orStatus(SelectionStatus::TIGHTISO);
  }
  if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
  return status;

}


bool noIso(ControlFlags &, vector<SimpleLepton> &, float rho) { 

        return true;
}


//--------------------------------------------------------------------------------------------------
SelectionStatus muonIsoMVASelection(ControlFlags &ctrl, 
                                    const mithep::Muon * mu, 
                                    const mithep::Vertex * vtx,
                                    const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                                    const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity, 
                                    mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
                                    vector<const mithep::Muon*> muonsToVeto,
                                    vector<const mithep::Electron*> electronsToVeto) 
//--------------------------------------------------------------------------------------------------
{ 

  if( ctrl.debug ) { 
    cout << "muonIsoMVASelection :: muons to veto " << endl;
    for( int i=0; i<muonsToVeto.size(); i++ ) { 
      const mithep::Muon * vmu = muonsToVeto[i];
      cout << "\tpt: " << vmu->Pt() 
           << "\teta: " << vmu->Eta() 
           << "\tphi: " << vmu->Phi() 
           << endl;
    }
    cout << "muonIsoMVASelection :: electrson to veto " << endl;
    for( int i=0; i<electronsToVeto.size(); i++ ) { 
      const mithep::Electron * vel = electronsToVeto[i];
      cout << "\tpt: " << vel->Pt() 
           << "\teta: " << vel->Eta() 
           << "\tphi: " << vel->Phi() 
           << endl;
    }
  }
  bool failiso=false;

  //
  // tmp iso rings
  //
  Double_t tmpChargedIso_DR0p0To0p1  = 0;
  Double_t tmpChargedIso_DR0p1To0p2  = 0;
  Double_t tmpChargedIso_DR0p2To0p3  = 0;
  Double_t tmpChargedIso_DR0p3To0p4  = 0;
  Double_t tmpChargedIso_DR0p4To0p5  = 0;
  Double_t tmpChargedIso_DR0p5To0p7  = 0;

  Double_t tmpGammaIso_DR0p0To0p1  = 0;
  Double_t tmpGammaIso_DR0p1To0p2  = 0;
  Double_t tmpGammaIso_DR0p2To0p3  = 0;
  Double_t tmpGammaIso_DR0p3To0p4  = 0;
  Double_t tmpGammaIso_DR0p4To0p5  = 0;
  Double_t tmpGammaIso_DR0p5To0p7  = 0;

  Double_t tmpNeutralHadronIso_DR0p0To0p1  = 0;
  Double_t tmpNeutralHadronIso_DR0p1To0p2  = 0;
  Double_t tmpNeutralHadronIso_DR0p2To0p3  = 0;
  Double_t tmpNeutralHadronIso_DR0p3To0p4  = 0;
  Double_t tmpNeutralHadronIso_DR0p4To0p5  = 0;
  Double_t tmpNeutralHadronIso_DR0p5To0p7  = 0;

        
  //
  // final rings for the MVA
  //
  Double_t fChargedIso_DR0p0To0p1;
  Double_t fChargedIso_DR0p1To0p2;
  Double_t fChargedIso_DR0p2To0p3;
  Double_t fChargedIso_DR0p3To0p4;
  Double_t fChargedIso_DR0p4To0p5;
  Double_t fChargedIso_DR0p5To0p7;

  Double_t fGammaIso_DR0p0To0p1;
  Double_t fGammaIso_DR0p1To0p2;
  Double_t fGammaIso_DR0p2To0p3;
  Double_t fGammaIso_DR0p3To0p4;
  Double_t fGammaIso_DR0p4To0p5;
  Double_t fGammaIso_DR0p5To0p7;

  Double_t fNeutralHadronIso_DR0p0To0p1;
  Double_t fNeutralHadronIso_DR0p1To0p2;
  Double_t fNeutralHadronIso_DR0p2To0p3;
  Double_t fNeutralHadronIso_DR0p3To0p4;
  Double_t fNeutralHadronIso_DR0p4To0p5;
  Double_t fNeutralHadronIso_DR0p5To0p7;


  //
  //Loop over PF Candidates
  //
  for(int k=0; k<fPFCandidates->GetEntries(); ++k) {

    if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack

    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);

    Double_t deta = (mu->Eta() - pf->Eta());
    Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
    Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
    if (dr > 1.0) continue;

    if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;

    //
    // Lepton Footprint Removal
    //
    Bool_t IsLeptonFootprint = kFALSE;
    if (dr < 1.0) {

      //
      // Check for electrons
      //
      for (Int_t q=0; q < electronsToVeto.size(); ++q) {
        const mithep::Electron *tmpele = electronsToVeto[q];
        // 4l electron
        if( pf->HasTrackerTrk() ) { 
          if( pf->TrackerTrk() == tmpele->TrackerTrk() )
            IsLeptonFootprint = kTRUE;
        }
        if( pf->HasGsfTrk() ) { 
          if( pf->GsfTrk() == tmpele->GsfTrk() )
            IsLeptonFootprint = kTRUE;
        }
        // PF charged
        if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) >= 1.479 
            && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) 
          IsLeptonFootprint = kTRUE;
        // PF gamma
        if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) >= 1.479 
            && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) 
          IsLeptonFootprint = kTRUE;
      } // loop over electrons 
      
      /* KH - commented for sync
      //
      // Check for muons
      //
      for (Int_t q=0; q < muonsToVeto.size(); ++q) {
        const mithep::Muon *tmpmu = muonsToVeto[q];
        // 4l muon
        if( pf->HasTrackerTrk() ) { 
          if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
            IsLeptonFootprint = kTRUE;
        }
        // PF charged
        if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) 
          IsLeptonFootprint = kTRUE;
      } // loop over muons
      */

    if (IsLeptonFootprint) 
      continue;

    //
    // Charged Iso Rings
    //
    if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {

      if( dr < 0.01 ) continue; // only for muon iso mva?
      if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;

//       if( pf->HasTrackerTrk() ) { 
//      if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
//      if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " " 
//                            << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
//                            << dr << endl;
//       }
//       if( pf->HasGsfTrk() ) { 
//      if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
//      if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " " 
//                            << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
//                            << dr << endl;
//       }

      // Footprint Veto
      if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
      if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
      if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
      if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
      if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
      if (dr >= 0.5 && dr < 0.7) tmpChargedIso_DR0p5To0p7 += pf->Pt();
    }

    //
    // Gamma Iso Rings
    //
    else if (abs(pf->PFType()) == PFCandidate::eGamma) {
      if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
      if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
      if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
      if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
      if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
      if (dr >= 0.5 && dr < 0.7) tmpGammaIso_DR0p5To0p7 += pf->Pt();
    }

    //
    // Other Neutral Iso Rings 
    //
    else {
      if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
      if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
      if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
      if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
      if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
      if (dr >= 0.5 && dr < 0.7) tmpNeutralHadronIso_DR0p5To0p7 += pf->Pt();
    }

    }

  }

  fChargedIso_DR0p0To0p1   = fmin((tmpChargedIso_DR0p0To0p1)/mu->Pt(), 2.5);
  fChargedIso_DR0p1To0p2   = fmin((tmpChargedIso_DR0p1To0p2)/mu->Pt(), 2.5);
  fChargedIso_DR0p2To0p3   = fmin((tmpChargedIso_DR0p2To0p3)/mu->Pt(), 2.5);
  fChargedIso_DR0p3To0p4   = fmin((tmpChargedIso_DR0p3To0p4)/mu->Pt(), 2.5);
  fChargedIso_DR0p4To0p5   = fmin((tmpChargedIso_DR0p4To0p5)/mu->Pt(), 2.5);


  double rho = 0;
  if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho()))) 
    rho = fPUEnergyDensity->At(0)->Rho();

//   if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta()))) 
//     rho = fPUEnergyDensity->At(0)->RhoLowEta();
  
  // WARNING!!!!  
  // hardcode for sync ...
  EffectiveAreaVersion = muT.kMuEAData2011;
  // WARNING!!!!  

 
  fGammaIso_DR0p0To0p1 = fmax(fmin((tmpGammaIso_DR0p0To0p1
                                  -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p0To0p1,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                 ,2.5)
                             ,0.0);
  fGammaIso_DR0p1To0p2 = fmax(fmin((tmpGammaIso_DR0p1To0p2
                                  -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p1To0p2,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                 ,2.5)
                             ,0.0);
  fGammaIso_DR0p2To0p3 = fmax(fmin((tmpGammaIso_DR0p2To0p3 
                                  -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p2To0p3,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                 ,2.5)
                             ,0.0);
  fGammaIso_DR0p3To0p4 = fmax(fmin((tmpGammaIso_DR0p3To0p4 
                                  -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p3To0p4,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                 ,2.5) 
                             ,0.0);
  fGammaIso_DR0p4To0p5 = fmax(fmin((tmpGammaIso_DR0p4To0p5 
                                  -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p4To0p5,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                 ,2.5) 
                             ,0.0);


  fNeutralHadronIso_DR0p0To0p1 = fmax(fmin((tmpNeutralHadronIso_DR0p0To0p1
                                          -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p0To0p1,
                                                                 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                         , 2.5)
                                     , 0.0);
  fNeutralHadronIso_DR0p1To0p2 = fmax(fmin((tmpNeutralHadronIso_DR0p1To0p2
                                            -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p1To0p2, 
                                                                   mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                           , 2.5)
                                       , 0.0);
  fNeutralHadronIso_DR0p2To0p3 = fmax(fmin((tmpNeutralHadronIso_DR0p2To0p3
                                          -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p2To0p3, 
                                                                 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                         , 2.5)
                                     , 0.0);
  fNeutralHadronIso_DR0p3To0p4 = fmax(fmin((tmpNeutralHadronIso_DR0p3To0p4 
                                          -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p3To0p4, 
                                                                 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
                                         , 2.5)
                                     , 0.0);
  fNeutralHadronIso_DR0p4To0p5 = fmax(fmin((tmpNeutralHadronIso_DR0p4To0p5 
                                          -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p4To0p5, 
                                                                 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
                                         , 2.5)
                                     , 0.0);


  double mvaval = muIsoMVA->MVAValue_IsoRings( mu->Pt(),
                                               mu->Eta(),
                                               mu->IsGlobalMuon(),
                                               mu->IsTrackerMuon(),
                                               fChargedIso_DR0p0To0p1,
                                               fChargedIso_DR0p1To0p2,
                                               fChargedIso_DR0p2To0p3,
                                               fChargedIso_DR0p3To0p4,
                                               fChargedIso_DR0p4To0p5,
                                               fGammaIso_DR0p0To0p1,
                                               fGammaIso_DR0p1To0p2,
                                               fGammaIso_DR0p2To0p3,
                                               fGammaIso_DR0p3To0p4,
                                               fGammaIso_DR0p4To0p5,
                                               fNeutralHadronIso_DR0p0To0p1,
                                               fNeutralHadronIso_DR0p1To0p2,
                                               fNeutralHadronIso_DR0p2To0p3,
                                               fNeutralHadronIso_DR0p3To0p4,
                                               fNeutralHadronIso_DR0p4To0p5,
                                               ctrl.debug);

  SelectionStatus status;
  bool pass;

  pass = false;
  if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
      && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_OPT_BIN0)   pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
           && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_OPT_BIN1)  pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
           && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_OPT_BIN2)  pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
           && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_OPT_BIN3)  pass = true;
  else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_OPT_BIN4)  pass = true;
  else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_OPT_BIN5)  pass = true;
  if( pass ) status.orStatus(SelectionStatus::LOOSEISO);

  /*
  pass = false;
  if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
      && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN0)   pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
           && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN1)  pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
           && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN2)  pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
           && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN3)  pass = true;
  else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN4)  pass = true;
  else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN5)  pass = true;
  if( pass ) status.orStatus(SelectionStatus::TIGHTISO);
  */

  //  pass &= (fChargedIso_DR0p0To0p1 + fChargedIso_DR0p1To0p2 + fChargedIso_DR0p2To0p3 < 0.7);

  status.isoMVA = mvaval;

  if(ctrl.debug)  {
    cout << "returning status : " << hex << status.getStatus() << dec << endl;
    cout << "MVAVAL : " << status.isoMVA << endl;
  }
  return status;

}


//--------------------------------------------------------------------------------------------------
SelectionStatus muonIsoMVASelection(ControlFlags &ctrl, 
                                    const mithep::Muon * mu, 
                                    const mithep::Vertex * vtx,
                                    const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                                    float rho,
                                    //const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity, 
                                    mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
                                    vector<const mithep::Muon*> muonsToVeto,
                                    vector<const mithep::Electron*> electronsToVeto) 
//--------------------------------------------------------------------------------------------------
// hacked version
{ 

  if( ctrl.debug ) { 
    cout << "muonIsoMVASelection :: muons to veto " << endl;
    for( int i=0; i<muonsToVeto.size(); i++ ) { 
      const mithep::Muon * vmu = muonsToVeto[i];
      cout << "\tpt: " << vmu->Pt() 
           << "\teta: " << vmu->Eta() 
           << "\tphi: " << vmu->Phi() 
           << endl;
    }
    cout << "muonIsoMVASelection :: electrson to veto " << endl;
    for( int i=0; i<electronsToVeto.size(); i++ ) { 
      const mithep::Electron * vel = electronsToVeto[i];
      cout << "\tpt: " << vel->Pt() 
           << "\teta: " << vel->Eta() 
           << "\tphi: " << vel->Phi() 
           << endl;
    }
  }
  bool failiso=false;

  //
  // tmp iso rings
  //
  Double_t tmpChargedIso_DR0p0To0p1  = 0;
  Double_t tmpChargedIso_DR0p1To0p2  = 0;
  Double_t tmpChargedIso_DR0p2To0p3  = 0;
  Double_t tmpChargedIso_DR0p3To0p4  = 0;
  Double_t tmpChargedIso_DR0p4To0p5  = 0;
  Double_t tmpChargedIso_DR0p5To0p7  = 0;

  Double_t tmpGammaIso_DR0p0To0p1  = 0;
  Double_t tmpGammaIso_DR0p1To0p2  = 0;
  Double_t tmpGammaIso_DR0p2To0p3  = 0;
  Double_t tmpGammaIso_DR0p3To0p4  = 0;
  Double_t tmpGammaIso_DR0p4To0p5  = 0;
  Double_t tmpGammaIso_DR0p5To0p7  = 0;

  Double_t tmpNeutralHadronIso_DR0p0To0p1  = 0;
  Double_t tmpNeutralHadronIso_DR0p1To0p2  = 0;
  Double_t tmpNeutralHadronIso_DR0p2To0p3  = 0;
  Double_t tmpNeutralHadronIso_DR0p3To0p4  = 0;
  Double_t tmpNeutralHadronIso_DR0p4To0p5  = 0;
  Double_t tmpNeutralHadronIso_DR0p5To0p7  = 0;

        
  //
  // final rings for the MVA
  //
  Double_t fChargedIso_DR0p0To0p1;
  Double_t fChargedIso_DR0p1To0p2;
  Double_t fChargedIso_DR0p2To0p3;
  Double_t fChargedIso_DR0p3To0p4;
  Double_t fChargedIso_DR0p4To0p5;
  Double_t fChargedIso_DR0p5To0p7;

  Double_t fGammaIso_DR0p0To0p1;
  Double_t fGammaIso_DR0p1To0p2;
  Double_t fGammaIso_DR0p2To0p3;
  Double_t fGammaIso_DR0p3To0p4;
  Double_t fGammaIso_DR0p4To0p5;
  Double_t fGammaIso_DR0p5To0p7;

  Double_t fNeutralHadronIso_DR0p0To0p1;
  Double_t fNeutralHadronIso_DR0p1To0p2;
  Double_t fNeutralHadronIso_DR0p2To0p3;
  Double_t fNeutralHadronIso_DR0p3To0p4;
  Double_t fNeutralHadronIso_DR0p4To0p5;
  Double_t fNeutralHadronIso_DR0p5To0p7;


  //
  //Loop over PF Candidates
  //
  for(int k=0; k<fPFCandidates->GetEntries(); ++k) {

    if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack

    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);

    Double_t deta = (mu->Eta() - pf->Eta());
    Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
    Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
    if (dr > 1.0) continue;

    if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;

    //
    // Lepton Footprint Removal
    //
    Bool_t IsLeptonFootprint = kFALSE;
    if (dr < 1.0) {

      //
      // Check for electrons
      //
      for (Int_t q=0; q < electronsToVeto.size(); ++q) {
        const mithep::Electron *tmpele = electronsToVeto[q];
        // 4l electron
        if( pf->HasTrackerTrk() ) { 
          if( pf->TrackerTrk() == tmpele->TrackerTrk() )
            IsLeptonFootprint = kTRUE;
        }
        if( pf->HasGsfTrk() ) { 
          if( pf->GsfTrk() == tmpele->GsfTrk() )
            IsLeptonFootprint = kTRUE;
        }
        // PF charged
        if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) >= 1.479 
            && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) 
          IsLeptonFootprint = kTRUE;
        // PF gamma
        if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) >= 1.479 
            && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) 
          IsLeptonFootprint = kTRUE;
      } // loop over electrons 
      
      /* KH - commented for sync
      //
      // Check for muons
      //
      for (Int_t q=0; q < muonsToVeto.size(); ++q) {
        const mithep::Muon *tmpmu = muonsToVeto[q];
        // 4l muon
        if( pf->HasTrackerTrk() ) { 
          if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
            IsLeptonFootprint = kTRUE;
        }
        // PF charged
        if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) 
          IsLeptonFootprint = kTRUE;
      } // loop over muons
      */

    if (IsLeptonFootprint) 
      continue;

    //
    // Charged Iso Rings
    //
    if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {

      if( dr < 0.01 ) continue; // only for muon iso mva?
      if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;

//       if( pf->HasTrackerTrk() ) { 
//      if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
//      if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " " 
//                            << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
//                            << dr << endl;
//       }
//       if( pf->HasGsfTrk() ) { 
//      if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
//      if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " " 
//                            << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
//                            << dr << endl;
//       }

      // Footprint Veto
      if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
      if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
      if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
      if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
      if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();
      if (dr >= 0.5 && dr < 0.7) tmpChargedIso_DR0p5To0p7 += pf->Pt();
    }

    //
    // Gamma Iso Rings
    //
    else if (abs(pf->PFType()) == PFCandidate::eGamma) {
      if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
      if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
      if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
      if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
      if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
      if (dr >= 0.5 && dr < 0.7) tmpGammaIso_DR0p5To0p7 += pf->Pt();
    }

    //
    // Other Neutral Iso Rings 
    //
    else {
      if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
      if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
      if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
      if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
      if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
      if (dr >= 0.5 && dr < 0.7) tmpNeutralHadronIso_DR0p5To0p7 += pf->Pt();
    }

    }

  }

  fChargedIso_DR0p0To0p1   = fmin((tmpChargedIso_DR0p0To0p1)/mu->Pt(), 2.5);
  fChargedIso_DR0p1To0p2   = fmin((tmpChargedIso_DR0p1To0p2)/mu->Pt(), 2.5);
  fChargedIso_DR0p2To0p3   = fmin((tmpChargedIso_DR0p2To0p3)/mu->Pt(), 2.5);
  fChargedIso_DR0p3To0p4   = fmin((tmpChargedIso_DR0p3To0p4)/mu->Pt(), 2.5);
  fChargedIso_DR0p4To0p5   = fmin((tmpChargedIso_DR0p4To0p5)/mu->Pt(), 2.5);


//   double rho = 0;
//   if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho()))) 
//     rho = fPUEnergyDensity->At(0)->Rho();
//   if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta()))) 
//     rho = fPUEnergyDensity->At(0)->RhoLowEta();
  
  // WARNING!!!!  
  // hardcode for sync ...
  EffectiveAreaVersion = muT.kMuEAData2011;
  // WARNING!!!!  

 
  fGammaIso_DR0p0To0p1 = fmax(fmin((tmpGammaIso_DR0p0To0p1
                                  -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p0To0p1,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                 ,2.5)
                             ,0.0);
  fGammaIso_DR0p1To0p2 = fmax(fmin((tmpGammaIso_DR0p1To0p2
                                  -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p1To0p2,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                 ,2.5)
                             ,0.0);
  fGammaIso_DR0p2To0p3 = fmax(fmin((tmpGammaIso_DR0p2To0p3 
                                  -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p2To0p3,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                 ,2.5)
                             ,0.0);
  fGammaIso_DR0p3To0p4 = fmax(fmin((tmpGammaIso_DR0p3To0p4 
                                  -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p3To0p4,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                 ,2.5) 
                             ,0.0);
  fGammaIso_DR0p4To0p5 = fmax(fmin((tmpGammaIso_DR0p4To0p5 
                                  -rho*muT.MuonEffectiveArea(muT.kMuGammaIsoDR0p4To0p5,mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                 ,2.5) 
                             ,0.0);


  fNeutralHadronIso_DR0p0To0p1 = fmax(fmin((tmpNeutralHadronIso_DR0p0To0p1
                                          -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p0To0p1,
                                                                 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                         , 2.5)
                                     , 0.0);
  fNeutralHadronIso_DR0p1To0p2 = fmax(fmin((tmpNeutralHadronIso_DR0p1To0p2
                                            -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p1To0p2, 
                                                                   mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                           , 2.5)
                                       , 0.0);
  fNeutralHadronIso_DR0p2To0p3 = fmax(fmin((tmpNeutralHadronIso_DR0p2To0p3
                                          -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p2To0p3, 
                                                                 mu->Eta(),EffectiveAreaVersion))/mu->Pt()
                                         , 2.5)
                                     , 0.0);
  fNeutralHadronIso_DR0p3To0p4 = fmax(fmin((tmpNeutralHadronIso_DR0p3To0p4 
                                          -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p3To0p4, 
                                                                 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
                                         , 2.5)
                                     , 0.0);
  fNeutralHadronIso_DR0p4To0p5 = fmax(fmin((tmpNeutralHadronIso_DR0p4To0p5 
                                          -rho*muT.MuonEffectiveArea(muT.kMuNeutralHadronIsoDR0p4To0p5, 
                                                                 mu->Eta(), EffectiveAreaVersion))/mu->Pt()
                                         , 2.5)
                                     , 0.0);


  double mvaval = muIsoMVA->MVAValue_IsoRings( mu->Pt(),
                                               mu->Eta(),
                                               mu->IsGlobalMuon(),
                                               mu->IsTrackerMuon(),
                                               fChargedIso_DR0p0To0p1,
                                               fChargedIso_DR0p1To0p2,
                                               fChargedIso_DR0p2To0p3,
                                               fChargedIso_DR0p3To0p4,
                                               fChargedIso_DR0p4To0p5,
                                               fGammaIso_DR0p0To0p1,
                                               fGammaIso_DR0p1To0p2,
                                               fGammaIso_DR0p2To0p3,
                                               fGammaIso_DR0p3To0p4,
                                               fGammaIso_DR0p4To0p5,
                                               fNeutralHadronIso_DR0p0To0p1,
                                               fNeutralHadronIso_DR0p1To0p2,
                                               fNeutralHadronIso_DR0p2To0p3,
                                               fNeutralHadronIso_DR0p3To0p4,
                                               fNeutralHadronIso_DR0p4To0p5,
                                               ctrl.debug);

  SelectionStatus status;
  bool pass;

  pass = false;
  if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
      && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN0)   pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
           && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN1)  pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
           && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN2)  pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
           && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN3)  pass = true;
  else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN4)  pass = true;
  else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_LOOSE_FORPFID_CUT_BIN5)  pass = true;
  if( pass ) status.orStatus(SelectionStatus::LOOSEISO);

  /*
  pass = false;
  if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
      && fabs(mu->Eta()) <= 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN0)   pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon()
           && fabs(mu->Eta()) <= 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN1)  pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
           && fabs(mu->Eta()) > 1.5 && mu->Pt() <= 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN2)  pass = true;
  else if( mu->IsGlobalMuon() && mu->IsTrackerMuon() 
           && fabs(mu->Eta()) > 1.5 && mu->Pt() > 10 && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN3)  pass = true;
  else if( !(mu->IsGlobalMuon()) && mu->IsTrackerMuon() && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN4)  pass = true;
  else if( mu->IsGlobalMuon() && !(mu->IsTrackerMuon()) && mvaval >= MUON_ISOMVA_TIGHT_FORPFID_CUT_BIN5)  pass = true;
  if( pass ) status.orStatus(SelectionStatus::TIGHTISO);
  */

  //  pass &= (fChargedIso_DR0p0To0p1 + fChargedIso_DR0p1To0p2 + fChargedIso_DR0p2To0p3 < 0.7);

  status.isoMVA = mvaval;

  if(ctrl.debug)  {
    cout << "returning status : " << hex << status.getStatus() << dec << endl;
    cout << "MVAVAL : " << status.isoMVA << endl;
  }
  return status;

}


//--------------------------------------------------------------------------------------------------
void initMuonIsoMVA() { 
//--------------------------------------------------------------------------------------------------
  muIsoMVA = new mithep::MuonIDMVA();
  vector<string> weightFiles;
  weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_barrel_lowpt.weights.xml");
  weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_barrel_highpt.weights.xml");
  weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_endcap_lowpt.weights.xml");
  weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_endcap_highpt.weights.xml");
  weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_tracker.weights.xml");
  weightFiles.push_back("./data/MuonIsoMVAWeights/MuonIsoMVA_BDTG_V0_global.weights.xml");
  muIsoMVA->Initialize( "MuonIsoMVA",
                        mithep::MuonIDMVA::kIsoRingsV0,
                        kTRUE, weightFiles);
}


//--------------------------------------------------------------------------------------------------
double  muonPFIso04(ControlFlags &ctrl, 
                    const mithep::Muon * mu, 
                    const mithep::Vertex * vtx,
                    const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                    const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity, 
                    mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
                    vector<const mithep::Muon*> muonsToVeto,
                    vector<const mithep::Electron*> electronsToVeto) 
//--------------------------------------------------------------------------------------------------
{ 

  extern double gChargedIso;  
  extern double  gGammaIso;      
  extern double  gNeutralIso;
  
  if( ctrl.debug ) { 
    cout << "muonIsoMVASelection :: muons to veto " << endl;
    for( int i=0; i<muonsToVeto.size(); i++ ) { 
      const mithep::Muon * vmu = muonsToVeto[i];
      cout << "\tpt: " << vmu->Pt() 
           << "\teta: " << vmu->Eta() 
           << "\tphi: " << vmu->Phi() 
           << endl;
    }
    cout << "muonIsoMVASelection :: electrson to veto " << endl;
    for( int i=0; i<electronsToVeto.size(); i++ ) { 
      const mithep::Electron * vel = electronsToVeto[i];
      cout << "\tpt: " << vel->Pt() 
           << "\teta: " << vel->Eta() 
           << "\tphi: " << vel->Phi() 
           << endl;
    }
  }

  //
  // final iso 
  //
  Double_t fChargedIso  = 0.0;
  Double_t fGammaIso  = 0.0;
  Double_t fNeutralHadronIso  = 0.0;

  //
  //Loop over PF Candidates
  //
  for(int k=0; k<fPFCandidates->GetEntries(); ++k) {

    if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack
    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);

    Double_t deta = (mu->Eta() - pf->Eta());
    Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
    Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
    if (dr > 0.4) continue;

    if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;

    //
    // Lepton Footprint Removal
    //
    Bool_t IsLeptonFootprint = kFALSE;
    if (dr < 1.0) {

      //
      // Check for electrons
      //
      for (Int_t q=0; q < electronsToVeto.size(); ++q) {
        const mithep::Electron *tmpele = electronsToVeto[q];
        // 4l electron
        if( pf->HasTrackerTrk() ) { 
          if( pf->TrackerTrk() == tmpele->TrackerTrk() )
            IsLeptonFootprint = kTRUE;
        }
        if( pf->HasGsfTrk() ) { 
          if( pf->GsfTrk() == tmpele->GsfTrk() )
            IsLeptonFootprint = kTRUE;
        }
        // PF charged
        if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479 
            && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
          if( ctrl.debug) cout << "\tcharged trk, dR (" 
                               << mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta())
                               << " matches 4L ele ..." << endl;
          IsLeptonFootprint = kTRUE;
        }
        // PF gamma
        if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479 
            && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) 
          IsLeptonFootprint = kTRUE;
      } // loop over electrons 

      /* KH - comment for sync      
      //
      // Check for muons
      //
      for (Int_t q=0; q < muonsToVeto.size(); ++q) {
        const mithep::Muon *tmpmu = muonsToVeto[q];
        // 4l muon
        if( pf->HasTrackerTrk() ) { 
          if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
            IsLeptonFootprint = kTRUE;
        }
        // PF charged
        if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) 
          IsLeptonFootprint = kTRUE;
      } // loop over muons
      */

    if (IsLeptonFootprint) 
      continue;

    //
    // Charged Iso 
    //
    if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {

      //if( dr < 0.01 ) continue; // only for muon iso mva?
      if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;


//       if( pf->HasTrackerTrk() ) { 
//      if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
//      if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " " 
//                            << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
//                            << dr << endl;
//       }
//       if( pf->HasGsfTrk() ) { 
//      if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
//      if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " " 
//                            << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
//                            << dr << endl;
//       }


      fChargedIso += pf->Pt();
    }

    //
    // Gamma Iso 
    //
    else if (abs(pf->PFType()) == PFCandidate::eGamma) {
      // KH, add to sync 
      if( pf->Pt() > 0.5 ) 
      fGammaIso += pf->Pt();
    }

    //
    // Other Neutrals
    //
    else {
      // KH, add to sync 
      if( pf->Pt() > 0.5 ) 
        fNeutralHadronIso += pf->Pt();
    }
    
    }
    
  }

  double rho=0;
  if( (EffectiveAreaVersion == mithep::MuonTools::kMuEAFall11MC) ||
      (EffectiveAreaVersion == mithep::MuonTools::kMuEAData2011) ) { 
    if (!(isnan(fPUEnergyDensity->At(0)->RhoKt6PFJetsForIso25()) || 
          isinf(fPUEnergyDensity->At(0)->RhoKt6PFJetsForIso25()))) 
      rho = fPUEnergyDensity->At(0)->RhoKt6PFJetsForIso25();
    //rho = fPUEnergyDensity->At(0)->Rho();
    // !!!!!!!!!!!!! TMP HACK FOR SYNC !!!!!!!!!!!!!!!!!!!!!
    EffectiveAreaVersion  = mithep::MuonTools::kMuEAData2011;
    // !!!!!!!!!!!!! TMP HACK FOR SYNC !!!!!!!!!!!!!!!!!!!!!
  } else { 
    if (!(isnan(fPUEnergyDensity->At(0)->RhoKt6PFJetsCentralNeutral()) || 
          isinf(fPUEnergyDensity->At(0)->RhoKt6PFJetsCentralNeutral()))) 
      rho = fPUEnergyDensity->At(0)->RhoKt6PFJetsCentralNeutral();
    // !!!!!!!!!!!!! TMP HACK FOR SYNC !!!!!!!!!!!!!!!!!!!!!
    EffectiveAreaVersion  = mithep::MuonTools::kMuEAData2012;
    // !!!!!!!!!!!!! TMP HACK FOR SYNC !!!!!!!!!!!!!!!!!!!!!
  }
  if(ctrl.debug) cout << "rho: " << rho << endl;

  double pfIso = fChargedIso + fmax(0.0,(fGammaIso + fNeutralHadronIso 
                                        -rho*muT.MuonEffectiveArea(muT.kMuGammaAndNeutralHadronIso04,
                                                                   mu->Eta(),EffectiveAreaVersion)));
  gChargedIso = fChargedIso;
  gGammaIso   = fGammaIso;
  gNeutralIso = fNeutralHadronIso;
  
  return pfIso;
}


//--------------------------------------------------------------------------------------------------
// hacked version
double  muonPFIso04(ControlFlags &ctrl, 
                    const mithep::Muon * mu, 
                    const mithep::Vertex * vtx,
                    const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                    float rho,
                    mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
                    vector<const mithep::Muon*> muonsToVeto,
                    vector<const mithep::Electron*> electronsToVeto) 
//--------------------------------------------------------------------------------------------------
{ 

  extern double gChargedIso;  
  extern double  gGammaIso;      
  extern double  gNeutralIso;
  
  if( ctrl.debug ) { 
    cout << "muonIsoMVASelection :: muons to veto " << endl;
    for( int i=0; i<muonsToVeto.size(); i++ ) { 
      const mithep::Muon * vmu = muonsToVeto[i];
      cout << "\tpt: " << vmu->Pt() 
           << "\teta: " << vmu->Eta() 
           << "\tphi: " << vmu->Phi() 
           << endl;
    }
    cout << "muonIsoMVASelection :: electrson to veto " << endl;
    for( int i=0; i<electronsToVeto.size(); i++ ) { 
      const mithep::Electron * vel = electronsToVeto[i];
      cout << "\tpt: " << vel->Pt() 
           << "\teta: " << vel->Eta() 
           << "\tphi: " << vel->Phi() 
           << endl;
    }
  }

  //
  // final iso 
  //
  Double_t fChargedIso  = 0.0;
  Double_t fGammaIso  = 0.0;
  Double_t fNeutralHadronIso  = 0.0;

  //
  //Loop over PF Candidates
  //
  for(int k=0; k<fPFCandidates->GetEntries(); ++k) {

    if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack
    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);

    Double_t deta = (mu->Eta() - pf->Eta());
    Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(mu->Phi()),Double_t(pf->Phi()));
    Double_t dr = mithep::MathUtils::DeltaR(mu->Phi(),mu->Eta(), pf->Phi(), pf->Eta());
    if (dr > 0.4) continue;

    if (pf->HasTrackerTrk() && (pf->TrackerTrk() == mu->TrackerTrk()) ) continue;

    //
    // Lepton Footprint Removal
    //
    Bool_t IsLeptonFootprint = kFALSE;
    if (dr < 1.0) {

      //
      // Check for electrons
      //
      for (Int_t q=0; q < electronsToVeto.size(); ++q) {
        const mithep::Electron *tmpele = electronsToVeto[q];
        // 4l electron
        if( pf->HasTrackerTrk() ) { 
          if( pf->TrackerTrk() == tmpele->TrackerTrk() )
            IsLeptonFootprint = kTRUE;
        }
        if( pf->HasGsfTrk() ) { 
          if( pf->GsfTrk() == tmpele->GsfTrk() )
            IsLeptonFootprint = kTRUE;
        }
        // PF charged
        if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479 
            && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) 
          IsLeptonFootprint = kTRUE;
        // PF gamma
        if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479 
            && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) 
          IsLeptonFootprint = kTRUE;
      } // loop over electrons 

      /* KH - comment for sync      
      //
      // Check for muons
      //
      for (Int_t q=0; q < muonsToVeto.size(); ++q) {
        const mithep::Muon *tmpmu = muonsToVeto[q];
        // 4l muon
        if( pf->HasTrackerTrk() ) { 
          if( pf->TrackerTrk() == tmpmu->TrackerTrk() )
            IsLeptonFootprint = kTRUE;
        }
        // PF charged
        if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) 
          IsLeptonFootprint = kTRUE;
      } // loop over muons
      */

    if (IsLeptonFootprint) 
      continue;

    //
    // Charged Iso 
    //
    if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {

      //if( dr < 0.01 ) continue; // only for muon iso mva?
      if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;


//       if( pf->HasTrackerTrk() ) { 
//      if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
//      if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " " 
//                            << abs(pf->TrackerTrk()->DzCorrected(vtx)) << " "
//                            << dr << endl;
//       }
//       if( pf->HasGsfTrk() ) { 
//      if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
//      if( ctrl.debug ) cout << "charged:: " << pf->PFType() << " " << pf->Pt() << " " 
//                            << abs(pf->GsfTrk()->DzCorrected(vtx)) << " "
//                            << dr << endl;
//       }


      fChargedIso += pf->Pt();
    }

    //
    // Gamma Iso 
    //
    else if (abs(pf->PFType()) == PFCandidate::eGamma) {
      // KH, add to sync 
      if( pf->Pt() > 0.5 ) 
      fGammaIso += pf->Pt();
    }

    //
    // Other Neutrals
    //
    else {
      // KH, add to sync 
      if( pf->Pt() > 0.5 ) 
        fNeutralHadronIso += pf->Pt();
    }
    
    }
    
  }
  
//   double rho = 0;
//   if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho()))) 
//     rho = fPUEnergyDensity->At(0)->Rho();

  // WARNING!!!!  
  // hardcode for sync ...
  EffectiveAreaVersion = muT.kMuEAData2011;
  // WARNING!!!!  


  double pfIso = fChargedIso + fmax(0.0,(fGammaIso + fNeutralHadronIso 
                                        -rho*muT.MuonEffectiveArea(muT.kMuGammaAndNeutralHadronIso04,
                                                                   mu->Eta(),EffectiveAreaVersion)));
  gChargedIso = fChargedIso;
  gGammaIso   = fGammaIso;
  gNeutralIso = fNeutralHadronIso;
  
  return pfIso;
}


//--------------------------------------------------------------------------------------------------
SelectionStatus muonReferenceIsoSelection(ControlFlags &ctrl, 
                                          const mithep::Muon * mu, 
                                          const mithep::Vertex * vtx,
                                          const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                                          const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity, 
                                          mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
                                          vector<const mithep::Muon*> muonsToVeto,
                                          vector<const mithep::Electron*> electronsToVeto) 
//--------------------------------------------------------------------------------------------------
{
  
  SelectionStatus status;

  double pfIso = muonPFIso04( ctrl, mu, vtx, fPFCandidates, fPUEnergyDensity, 
                              EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
  //  cout << "--------------> setting muon isoPF04 to" << pfIso << endl;
  status.isoPF04 = pfIso;
  status.chisoPF04 = gChargedIso;
  status.gaisoPF04 = gGammaIso;
  status.neisoPF04 = gNeutralIso;

  bool pass = false;
  if( (pfIso/mu->Pt()) < MUON_REFERENCE_PFISO_CUT ) pass = true;
  
  if( pass ) {
    status.orStatus(SelectionStatus::LOOSEISO);
    status.orStatus(SelectionStatus::TIGHTISO);
  }
  if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
  return status;
  
}


//--------------------------------------------------------------------------------------------------
// hacked version
SelectionStatus muonReferenceIsoSelection(ControlFlags &ctrl, 
                                          const mithep::Muon * mu, 
                                          const mithep::Vertex * vtx,
                                          const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                                          float rho,
                                          mithep::MuonTools::EMuonEffectiveAreaTarget EffectiveAreaVersion,
                                          vector<const mithep::Muon*> muonsToVeto,
                                          vector<const mithep::Electron*> electronsToVeto) 
//--------------------------------------------------------------------------------------------------
{
  
  SelectionStatus status;
  
  double pfIso = muonPFIso04( ctrl, mu, vtx, fPFCandidates, rho, 
                              EffectiveAreaVersion, muonsToVeto ,electronsToVeto );

  status.isoPF04 = pfIso;
  status.chisoPF04 = gChargedIso;
  status.gaisoPF04 = gGammaIso;
  status.neisoPF04 = gNeutralIso;

  bool pass = false;
  if( (pfIso/mu->Pt()) < MUON_REFERENCE_PFISO_CUT ) pass = true;
  
  if( pass ) {
    status.orStatus(SelectionStatus::LOOSEISO);
    status.orStatus(SelectionStatus::TIGHTISO);
  }
  if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
  return status;
  
}


//--------------------------------------------------------------------------------------------------
SelectionStatus electronIsoMVASelection(ControlFlags &ctrl, 
                                        const mithep::Electron * ele, 
                                        const mithep::Vertex * vtx,
                                        const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                                        const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity, 
                                        mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
                                        vector<const mithep::Muon*> muonsToVeto,
                                        vector<const mithep::Electron*> electronsToVeto) 
//--------------------------------------------------------------------------------------------------
{ 

  if( ctrl.debug ) { 
    cout << "electronIsoMVASelection :: muons to veto " << endl;
    for( int i=0; i<muonsToVeto.size(); i++ ) { 
      const mithep::Muon * vmu = muonsToVeto[i];
      cout << "\tpt: " << vmu->Pt() 
           << "\teta: " << vmu->Eta() 
           << "\tphi: " << vmu->Phi() 
           << endl;
    }
    cout << "electronIsoMVASelection :: electrson to veto " << endl;
    for( int i=0; i<electronsToVeto.size(); i++ ) { 
      const mithep::Electron * vel = electronsToVeto[i];
      cout << "\tpt: " << vel->Pt() 
           << "\teta: " << vel->Eta() 
           << "\tphi: " << vel->Phi() 
           << "\ttrk: " << vel->TrackerTrk()
           << endl;
    }
  }

  bool failiso=false;

  //
  // tmp iso rings
  //
  Double_t tmpChargedIso_DR0p0To0p1  = 0;
  Double_t tmpChargedIso_DR0p1To0p2  = 0;
  Double_t tmpChargedIso_DR0p2To0p3  = 0;
  Double_t tmpChargedIso_DR0p3To0p4  = 0;
  Double_t tmpChargedIso_DR0p4To0p5  = 0;

  Double_t tmpGammaIso_DR0p0To0p1  = 0;
  Double_t tmpGammaIso_DR0p1To0p2  = 0;
  Double_t tmpGammaIso_DR0p2To0p3  = 0;
  Double_t tmpGammaIso_DR0p3To0p4  = 0;
  Double_t tmpGammaIso_DR0p4To0p5  = 0;


  Double_t tmpNeutralHadronIso_DR0p0To0p1  = 0;
  Double_t tmpNeutralHadronIso_DR0p1To0p2  = 0;
  Double_t tmpNeutralHadronIso_DR0p2To0p3  = 0;
  Double_t tmpNeutralHadronIso_DR0p3To0p4  = 0;
  Double_t tmpNeutralHadronIso_DR0p4To0p5  = 0;

        
  //
  // final rings for the MVA
  //
  Double_t fChargedIso_DR0p0To0p1;
  Double_t fChargedIso_DR0p1To0p2;
  Double_t fChargedIso_DR0p2To0p3;
  Double_t fChargedIso_DR0p3To0p4;
  Double_t fChargedIso_DR0p4To0p5;

  Double_t fGammaIso_DR0p0To0p1;
  Double_t fGammaIso_DR0p1To0p2;
  Double_t fGammaIso_DR0p2To0p3;
  Double_t fGammaIso_DR0p3To0p4;
  Double_t fGammaIso_DR0p4To0p5;

  Double_t fNeutralHadronIso_DR0p0To0p1;
  Double_t fNeutralHadronIso_DR0p1To0p2;
  Double_t fNeutralHadronIso_DR0p2To0p3;
  Double_t fNeutralHadronIso_DR0p3To0p4;
  Double_t fNeutralHadronIso_DR0p4To0p5;


  //
  //Loop over PF Candidates
  //
  for(int k=0; k<fPFCandidates->GetEntries(); ++k) {

    if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack

    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
    Double_t deta = (ele->Eta() - pf->Eta());
    Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
    Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
    if (dr > 1.0) continue;

    if(ctrl.debug) { 
      cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt();
      if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(*vtx);
      cout << endl;
    }


    if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
         (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) continue;
    

    //
    // Lepton Footprint Removal
    //
    Bool_t IsLeptonFootprint = kFALSE;
    if (dr < 1.0) {


      //
      // Check for electrons
      //

      for (Int_t q=0; q < electronsToVeto.size(); ++q) {
        const mithep::Electron *tmpele = electronsToVeto[q];
        double tmpdr = mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta());

        // 4l electron
        if( pf->HasTrackerTrk()  ) { 
          if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
            if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
            IsLeptonFootprint = kTRUE;
          }
        }
        if( pf->HasGsfTrk()  ) { 
          if( pf->GsfTrk() == tmpele->GsfTrk() ) {
            if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
            IsLeptonFootprint = kTRUE;
          }
        }
        // PF charged
        if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) >= 1.479 && tmpdr < 0.015) {
          if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
          IsLeptonFootprint = kTRUE;
        }
        // PF gamma
        if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) >= 1.479 
            && tmpdr < 0.08) {
          if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
          IsLeptonFootprint = kTRUE;
        }
      } // loop over electrons 


      /* KH - comment for sync            
      //
      // Check for muons
      //
      for (Int_t q=0; q < muonsToVeto.size(); ++q) {
        const mithep::Muon *tmpmu = muonsToVeto[q];
        // 4l muon
        if( pf->HasTrackerTrk() ) { 
          if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
            if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
            IsLeptonFootprint = kTRUE;
          }
        }
        // PF charged
        if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
          if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
          IsLeptonFootprint = kTRUE;
        }
      } // loop over muons
      */

    if (IsLeptonFootprint) 
      continue;

    //
    // Charged Iso Rings
    //
    if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {

//       if( pf->HasGsfTrk() ) {
//       if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
//       } else if( pf->HasTrackerTrk() ){
//      if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
//       }

      // Veto any PFmuon, or PFEle
      if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;

      // Footprint Veto
      if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;

      if( ctrl.debug) cout << "charged:: pt: " << pf->Pt() 
                           << "\ttype: " << pf->PFType() 
                           << "\ttrk: " << pf->TrackerTrk() << endl;

      if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
      if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
      if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
      if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
      if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();

    }

    //
    // Gamma Iso Rings
    //
    else if (abs(pf->PFType()) == PFCandidate::eGamma) {

      if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.08) continue;

      if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " " 
                           << dr << endl;

      if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
      if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
      if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
      if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
      if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
    }

    //
    // Other Neutral Iso Rings 
    //
    else {
      if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " " 
                           << dr << endl;
      if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
      if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
      if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
      if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
      if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
    }

    }

  }

  fChargedIso_DR0p0To0p1   = fmin((tmpChargedIso_DR0p0To0p1)/ele->Pt(), 2.5);
  fChargedIso_DR0p1To0p2   = fmin((tmpChargedIso_DR0p1To0p2)/ele->Pt(), 2.5);
  fChargedIso_DR0p2To0p3   = fmin((tmpChargedIso_DR0p2To0p3)/ele->Pt(), 2.5);
  fChargedIso_DR0p3To0p4   = fmin((tmpChargedIso_DR0p3To0p4)/ele->Pt(), 2.5);
  fChargedIso_DR0p4To0p5   = fmin((tmpChargedIso_DR0p4To0p5)/ele->Pt(), 2.5);

  if(ctrl.debug) {
    cout << "fChargedIso_DR0p0To0p1 : " << fChargedIso_DR0p0To0p1  << endl;
    cout << "fChargedIso_DR0p1To0p2 : " << fChargedIso_DR0p1To0p2  << endl; 
    cout << "fChargedIso_DR0p2To0p3 : " << fChargedIso_DR0p2To0p3  << endl;
    cout << "fChargedIso_DR0p3To0p4 : " << fChargedIso_DR0p3To0p4  << endl;
    cout << "fChargedIso_DR0p4To0p5 : " << fChargedIso_DR0p4To0p5  << endl;
  }


  double rho = 0;
  if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho()))) 
    rho = fPUEnergyDensity->At(0)->Rho();
  //   if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta()))) 
  //     rho = fPUEnergyDensity->At(0)->RhoLowEta();
  
  // WARNING!!!!  
  // hardcode for sync ...
  EffectiveAreaVersion = eleT.kEleEAData2011;
  // WARNING!!!!  

  if( ctrl.debug) { 
    cout << "RHO: " << rho << endl;
    cout << "eta: " << ele->SCluster()->Eta() << endl;
    cout << "target: " << EffectiveAreaVersion << endl;
    cout << "effA 0-1: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
                                                       ele->SCluster()->Eta(),
                                                       EffectiveAreaVersion) 
         << endl;
    cout << "effA 1-2: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2,
                                                       ele->SCluster()->Eta(),
                                                       EffectiveAreaVersion) 
         << endl;
    cout << "effA 2-3: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3,
                                                       ele->SCluster()->Eta(),
                                                       EffectiveAreaVersion) 
         << endl;
    cout << "effA 3-4: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4,
                                                       ele->SCluster()->Eta(),
                                                       EffectiveAreaVersion) 
         << endl;
  }

  fGammaIso_DR0p0To0p1 = fmax(fmin((tmpGammaIso_DR0p0To0p1 
                                  -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p0To0p1,
                                                              ele->SCluster()->Eta(),
                                                              EffectiveAreaVersion))/ele->Pt()
                                 ,2.5)
                             ,0.0);
  fGammaIso_DR0p1To0p2 = fmax(fmin((tmpGammaIso_DR0p1To0p2
                                  -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p1To0p2,
                                                              ele->SCluster()->Eta(),
                                                              EffectiveAreaVersion))/ele->Pt()
                                 ,2.5)
                             ,0.0);
  fGammaIso_DR0p2To0p3 = fmax(fmin((tmpGammaIso_DR0p2To0p3
                                  -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p2To0p3,
                                                              ele->SCluster()->Eta()
                                                              ,EffectiveAreaVersion))/ele->Pt()
                                 ,2.5)
                             ,0.0);
  fGammaIso_DR0p3To0p4 = fmax(fmin((tmpGammaIso_DR0p3To0p4 
                                  -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p3To0p4,
                                                              ele->SCluster()->Eta(),
                                                              EffectiveAreaVersion))/ele->Pt()
                                 ,2.5) 
                             ,0.0);
  fGammaIso_DR0p4To0p5 = fmax(fmin((tmpGammaIso_DR0p4To0p5 
                                  -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p4To0p5,
                                                              ele->SCluster()->Eta(),
                                                              EffectiveAreaVersion))/ele->Pt()
                                 ,2.5) 
                             ,0.0);


  if( ctrl.debug) {
    cout << "fGammaIso_DR0p0To0p1: " << fGammaIso_DR0p0To0p1 << endl;
    cout << "fGammaIso_DR0p1To0p2: " << fGammaIso_DR0p1To0p2 << endl;
    cout << "fGammaIso_DR0p2To0p3: " << fGammaIso_DR0p2To0p3 << endl;
    cout << "fGammaIso_DR0p3To0p4: " << fGammaIso_DR0p3To0p4 << endl;
    cout << "fGammaIso_DR0p4To0p5: " << fGammaIso_DR0p4To0p5 << endl;
  }

  fNeutralHadronIso_DR0p0To0p1 = fmax(fmin((tmpNeutralHadronIso_DR0p0To0p1
                                          -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
                                                                 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
                                         , 2.5)
                                     , 0.0);
  fNeutralHadronIso_DR0p1To0p2 = fmax(fmin((tmpNeutralHadronIso_DR0p1To0p2
                                            -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2, 
                                                                   ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
                                           , 2.5)
                                       , 0.0);
  fNeutralHadronIso_DR0p2To0p3 = fmax(fmin((tmpNeutralHadronIso_DR0p2To0p3
                                          -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3, 
                                                                 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
                                         , 2.5)
                                     , 0.0);
  fNeutralHadronIso_DR0p3To0p4 = fmax(fmin((tmpNeutralHadronIso_DR0p3To0p4 
                                          -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4, 
                                                                 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
                                         , 2.5)
                                     , 0.0);
  fNeutralHadronIso_DR0p4To0p5 = fmax(fmin((tmpNeutralHadronIso_DR0p4To0p5 
                                          -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p4To0p5, 
                                                                 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
                                         , 2.5)
                                     , 0.0);

  if( ctrl.debug) {
    cout << "fNeutralHadronIso_DR0p0To0p1: " << fNeutralHadronIso_DR0p0To0p1 << endl;
    cout << "fNeutralHadronIso_DR0p1To0p2: " << fNeutralHadronIso_DR0p1To0p2 << endl;
    cout << "fNeutralHadronIso_DR0p2To0p3: " << fNeutralHadronIso_DR0p2To0p3 << endl;
    cout << "fNeutralHadronIso_DR0p3To0p4: " << fNeutralHadronIso_DR0p3To0p4 << endl;
    cout << "fNeutralHadronIso_DR0p4To0p5: " << fNeutralHadronIso_DR0p4To0p5 << endl;
  }

  double mvaval = eleIsoMVA->MVAValue_IsoRings( ele->Pt(),
                                                ele->SCluster()->Eta(),
                                                fChargedIso_DR0p0To0p1,
                                                fChargedIso_DR0p1To0p2,
                                                fChargedIso_DR0p2To0p3,
                                                fChargedIso_DR0p3To0p4,
                                                fChargedIso_DR0p4To0p5,
                                                fGammaIso_DR0p0To0p1,
                                                fGammaIso_DR0p1To0p2,
                                                fGammaIso_DR0p2To0p3,
                                                fGammaIso_DR0p3To0p4,
                                                fGammaIso_DR0p4To0p5,
                                                fNeutralHadronIso_DR0p0To0p1,
                                                fNeutralHadronIso_DR0p1To0p2,
                                                fNeutralHadronIso_DR0p2To0p3,
                                                fNeutralHadronIso_DR0p3To0p4,
                                                fNeutralHadronIso_DR0p4To0p5,
                                                ctrl.debug);

  SelectionStatus status;
  status.isoMVA = mvaval;
  bool pass = false;

  Int_t subdet = 0;
  if (fabs(ele->SCluster()->Eta()) < 0.8) subdet = 0;
  else if (fabs(ele->SCluster()->Eta()) < 1.479) subdet = 1;
  else subdet = 2;

  Int_t ptBin = 0;
  if (ele->Pt() >= 10.0) ptBin = 1;
  
  Int_t MVABin = -1;
  if (subdet == 0 && ptBin == 0) MVABin = 0;
  if (subdet == 1 && ptBin == 0) MVABin = 1;
  if (subdet == 2 && ptBin == 0) MVABin = 2;
  if (subdet == 0 && ptBin == 1) MVABin = 3;
  if (subdet == 1 && ptBin == 1) MVABin = 4;
  if (subdet == 2 && ptBin == 1) MVABin = 5;

  pass = false;
  if( MVABin == 0 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_OPT_BIN0 ) pass = true;
  if( MVABin == 1 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_OPT_BIN1 ) pass = true;
  if( MVABin == 2 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_OPT_BIN2 ) pass = true;
  if( MVABin == 3 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_OPT_BIN3 ) pass = true;
  if( MVABin == 4 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_OPT_BIN4 ) pass = true;
  if( MVABin == 5 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_OPT_BIN5 ) pass = true;
  //  pass &= (fChargedIso_DR0p0To0p1 + fChargedIso_DR0p1To0p2 + fChargedIso_DR0p2To0p3 < 0.7);
  if( pass ) status.orStatus(SelectionStatus::LOOSEISO);

//   pass = false;
//   if( MVABin == 0 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN0 ) pass = true;
//   if( MVABin == 1 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN1 ) pass = true;
//   if( MVABin == 2 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN2 ) pass = true;
//   if( MVABin == 3 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN3 ) pass = true;
//   if( MVABin == 4 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN4 ) pass = true;
//   if( MVABin == 5 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN5 ) pass = true;
//   if( pass ) status.orStatus(SelectionStatus::TIGHTISO);

  if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
  return status;
  
}


//--------------------------------------------------------------------------------------------------
SelectionStatus electronIsoMVASelection(ControlFlags &ctrl, 
                                        const mithep::Electron * ele, 
                                        const mithep::Vertex * vtx,
                                        const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                                        float rho,
                                        //const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity, 
                                        mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
                                        vector<const mithep::Muon*> muonsToVeto,
                                        vector<const mithep::Electron*> electronsToVeto) 
//--------------------------------------------------------------------------------------------------
// hacked version
{ 
  if( ctrl.debug ) { 
    cout << "================> hacked ele Iso MVA <======================" << endl; 
  }

  if( ctrl.debug ) { 
    cout << "electronIsoMVASelection :: muons to veto " << endl;
    for( int i=0; i<muonsToVeto.size(); i++ ) { 
      const mithep::Muon * vmu = muonsToVeto[i];
      cout << "\tpt: " << vmu->Pt() 
           << "\teta: " << vmu->Eta() 
           << "\tphi: " << vmu->Phi() 
           << endl;
    }
    cout << "electronIsoMVASelection :: electrson to veto " << endl;
    for( int i=0; i<electronsToVeto.size(); i++ ) { 
      const mithep::Electron * vel = electronsToVeto[i];
      cout << "\tpt: " << vel->Pt() 
           << "\teta: " << vel->Eta() 
           << "\tphi: " << vel->Phi() 
           << "\ttrk: " << vel->TrackerTrk()
           << endl;
    }
  }

  bool failiso=false;

  //
  // tmp iso rings
  //
  Double_t tmpChargedIso_DR0p0To0p1  = 0;
  Double_t tmpChargedIso_DR0p1To0p2  = 0;
  Double_t tmpChargedIso_DR0p2To0p3  = 0;
  Double_t tmpChargedIso_DR0p3To0p4  = 0;
  Double_t tmpChargedIso_DR0p4To0p5  = 0;

  Double_t tmpGammaIso_DR0p0To0p1  = 0;
  Double_t tmpGammaIso_DR0p1To0p2  = 0;
  Double_t tmpGammaIso_DR0p2To0p3  = 0;
  Double_t tmpGammaIso_DR0p3To0p4  = 0;
  Double_t tmpGammaIso_DR0p4To0p5  = 0;


  Double_t tmpNeutralHadronIso_DR0p0To0p1  = 0;
  Double_t tmpNeutralHadronIso_DR0p1To0p2  = 0;
  Double_t tmpNeutralHadronIso_DR0p2To0p3  = 0;
  Double_t tmpNeutralHadronIso_DR0p3To0p4  = 0;
  Double_t tmpNeutralHadronIso_DR0p4To0p5  = 0;

        
  //
  // final rings for the MVA
  //
  Double_t fChargedIso_DR0p0To0p1;
  Double_t fChargedIso_DR0p1To0p2;
  Double_t fChargedIso_DR0p2To0p3;
  Double_t fChargedIso_DR0p3To0p4;
  Double_t fChargedIso_DR0p4To0p5;

  Double_t fGammaIso_DR0p0To0p1;
  Double_t fGammaIso_DR0p1To0p2;
  Double_t fGammaIso_DR0p2To0p3;
  Double_t fGammaIso_DR0p3To0p4;
  Double_t fGammaIso_DR0p4To0p5;

  Double_t fNeutralHadronIso_DR0p0To0p1;
  Double_t fNeutralHadronIso_DR0p1To0p2;
  Double_t fNeutralHadronIso_DR0p2To0p3;
  Double_t fNeutralHadronIso_DR0p3To0p4;
  Double_t fNeutralHadronIso_DR0p4To0p5;


  //
  //Loop over PF Candidates
  //
  for(int k=0; k<fPFCandidates->GetEntries(); ++k) {

    if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack

    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
    Double_t deta = (ele->Eta() - pf->Eta());
    Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
    Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());
    if (dr > 1.0) continue;

    if(ctrl.debug) { 
      cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt();
      if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(*vtx);
      cout << endl;
    }


    if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
         (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) continue;
    

    //
    // Lepton Footprint Removal
    //
    Bool_t IsLeptonFootprint = kFALSE;
    if (dr < 1.0) {


      //
      // Check for electrons
      //

      for (Int_t q=0; q < electronsToVeto.size(); ++q) {
        const mithep::Electron *tmpele = electronsToVeto[q];
        double tmpdr = mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta());

        // 4l electron
        if( pf->HasTrackerTrk()  ) { 
          if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
            if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
            IsLeptonFootprint = kTRUE;
          }
        }
        if( pf->HasGsfTrk()  ) { 
          if( pf->GsfTrk() == tmpele->GsfTrk() ) {
            if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
            IsLeptonFootprint = kTRUE;
          }
        }
        // PF charged
        if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) >= 1.479 && tmpdr < 0.015) {
          if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
          IsLeptonFootprint = kTRUE;
        }
        // PF gamma
        if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) >= 1.479 
            && tmpdr < 0.08) {
          if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
          IsLeptonFootprint = kTRUE;
        }
      } // loop over electrons 


      /* KH - comment for sync            
      //
      // Check for muons
      //
      for (Int_t q=0; q < muonsToVeto.size(); ++q) {
        const mithep::Muon *tmpmu = muonsToVeto[q];
        // 4l muon
        if( pf->HasTrackerTrk() ) { 
          if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
            if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
            IsLeptonFootprint = kTRUE;
          }
        }
        // PF charged
        if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
          if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
          IsLeptonFootprint = kTRUE;
        }
      } // loop over muons
      */

    if (IsLeptonFootprint) 
      continue;

    //
    // Charged Iso Rings
    //
    if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {

//       if( pf->HasGsfTrk() ) {
//       if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;
//       } else if( pf->HasTrackerTrk() ){
//      if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
//       }

      // Veto any PFmuon, or PFEle
      if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) continue;

      // Footprint Veto
      if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;

      if( ctrl.debug) cout << "charged:: pt: " << pf->Pt() 
                           << "\ttype: " << pf->PFType() 
                           << "\ttrk: " << pf->TrackerTrk() << endl;

      if (dr < 0.1) tmpChargedIso_DR0p0To0p1 += pf->Pt();
      if (dr >= 0.1 && dr < 0.2) tmpChargedIso_DR0p1To0p2 += pf->Pt();
      if (dr >= 0.2 && dr < 0.3) tmpChargedIso_DR0p2To0p3 += pf->Pt();
      if (dr >= 0.3 && dr < 0.4) tmpChargedIso_DR0p3To0p4 += pf->Pt();
      if (dr >= 0.4 && dr < 0.5) tmpChargedIso_DR0p4To0p5 += pf->Pt();

    }

    //
    // Gamma Iso Rings
    //
    else if (abs(pf->PFType()) == PFCandidate::eGamma) {

      if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.08) continue;

      if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " " 
                           << dr << endl;

      if (dr < 0.1) tmpGammaIso_DR0p0To0p1 += pf->Pt();
      if (dr >= 0.1 && dr < 0.2) tmpGammaIso_DR0p1To0p2 += pf->Pt();
      if (dr >= 0.2 && dr < 0.3) tmpGammaIso_DR0p2To0p3 += pf->Pt();
      if (dr >= 0.3 && dr < 0.4) tmpGammaIso_DR0p3To0p4 += pf->Pt();
      if (dr >= 0.4 && dr < 0.5) tmpGammaIso_DR0p4To0p5 += pf->Pt();
    }

    //
    // Other Neutral Iso Rings 
    //
    else {
      if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " " 
                           << dr << endl;
      if (dr < 0.1) tmpNeutralHadronIso_DR0p0To0p1 += pf->Pt();
      if (dr >= 0.1 && dr < 0.2) tmpNeutralHadronIso_DR0p1To0p2 += pf->Pt();
      if (dr >= 0.2 && dr < 0.3) tmpNeutralHadronIso_DR0p2To0p3 += pf->Pt();
      if (dr >= 0.3 && dr < 0.4) tmpNeutralHadronIso_DR0p3To0p4 += pf->Pt();
      if (dr >= 0.4 && dr < 0.5) tmpNeutralHadronIso_DR0p4To0p5 += pf->Pt();
    }

    }

  }

  fChargedIso_DR0p0To0p1   = fmin((tmpChargedIso_DR0p0To0p1)/ele->Pt(), 2.5);
  fChargedIso_DR0p1To0p2   = fmin((tmpChargedIso_DR0p1To0p2)/ele->Pt(), 2.5);
  fChargedIso_DR0p2To0p3   = fmin((tmpChargedIso_DR0p2To0p3)/ele->Pt(), 2.5);
  fChargedIso_DR0p3To0p4   = fmin((tmpChargedIso_DR0p3To0p4)/ele->Pt(), 2.5);
  fChargedIso_DR0p4To0p5   = fmin((tmpChargedIso_DR0p4To0p5)/ele->Pt(), 2.5);

  if(ctrl.debug) {
    cout << "fChargedIso_DR0p0To0p1 : " << fChargedIso_DR0p0To0p1  << endl;
    cout << "fChargedIso_DR0p1To0p2 : " << fChargedIso_DR0p1To0p2  << endl; 
    cout << "fChargedIso_DR0p2To0p3 : " << fChargedIso_DR0p2To0p3  << endl;
    cout << "fChargedIso_DR0p3To0p4 : " << fChargedIso_DR0p3To0p4  << endl;
    cout << "fChargedIso_DR0p4To0p5 : " << fChargedIso_DR0p4To0p5  << endl;
  }


  //  rho=0;
  //  double rho = 0;
  //   if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho()))) 
  //     rho = fPUEnergyDensity->At(0)->Rho();
  //   if (!(isnan(fPUEnergyDensity->At(0)->RhoLowEta()) || isinf(fPUEnergyDensity->At(0)->RhoLowEta()))) 
  //     rho = fPUEnergyDensity->At(0)->RhoLowEta();
  
  // WARNING!!!!  
  // hardcode for sync ...
  EffectiveAreaVersion = eleT.kEleEAData2011;
  // WARNING!!!!  

  if( ctrl.debug) { 
    cout << "RHO: " << rho << endl;
    cout << "eta: " << ele->SCluster()->Eta() << endl;
    cout << "target: " << EffectiveAreaVersion << endl;
    cout << "effA 0-1: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
                                                       ele->SCluster()->Eta(),
                                                       EffectiveAreaVersion) 
         << endl;
    cout << "effA 1-2: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2,
                                                       ele->SCluster()->Eta(),
                                                       EffectiveAreaVersion) 
         << endl;
    cout << "effA 2-3: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3,
                                                       ele->SCluster()->Eta(),
                                                       EffectiveAreaVersion) 
         << endl;
    cout << "effA 3-4: " << eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4,
                                                       ele->SCluster()->Eta(),
                                                       EffectiveAreaVersion) 
         << endl;
  }

  fGammaIso_DR0p0To0p1 = fmax(fmin((tmpGammaIso_DR0p0To0p1 
                                  -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p0To0p1,
                                                              ele->SCluster()->Eta(),
                                                              EffectiveAreaVersion))/ele->Pt()
                                 ,2.5)
                             ,0.0);
  fGammaIso_DR0p1To0p2 = fmax(fmin((tmpGammaIso_DR0p1To0p2
                                  -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p1To0p2,
                                                              ele->SCluster()->Eta(),
                                                              EffectiveAreaVersion))/ele->Pt()
                                 ,2.5)
                             ,0.0);
  fGammaIso_DR0p2To0p3 = fmax(fmin((tmpGammaIso_DR0p2To0p3
                                  -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p2To0p3,
                                                              ele->SCluster()->Eta()
                                                              ,EffectiveAreaVersion))/ele->Pt()
                                 ,2.5)
                             ,0.0);
  fGammaIso_DR0p3To0p4 = fmax(fmin((tmpGammaIso_DR0p3To0p4 
                                  -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p3To0p4,
                                                              ele->SCluster()->Eta(),
                                                              EffectiveAreaVersion))/ele->Pt()
                                 ,2.5) 
                             ,0.0);
  fGammaIso_DR0p4To0p5 = fmax(fmin((tmpGammaIso_DR0p4To0p5 
                                  -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaIsoDR0p4To0p5,
                                                              ele->SCluster()->Eta(),
                                                              EffectiveAreaVersion))/ele->Pt()
                                 ,2.5) 
                             ,0.0);


  if( ctrl.debug) {
    cout << "fGammaIso_DR0p0To0p1: " << fGammaIso_DR0p0To0p1 << endl;
    cout << "fGammaIso_DR0p1To0p2: " << fGammaIso_DR0p1To0p2 << endl;
    cout << "fGammaIso_DR0p2To0p3: " << fGammaIso_DR0p2To0p3 << endl;
    cout << "fGammaIso_DR0p3To0p4: " << fGammaIso_DR0p3To0p4 << endl;
    cout << "fGammaIso_DR0p4To0p5: " << fGammaIso_DR0p4To0p5 << endl;
  }

  fNeutralHadronIso_DR0p0To0p1 = fmax(fmin((tmpNeutralHadronIso_DR0p0To0p1
                                          -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p0To0p1,
                                                                 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
                                         , 2.5)
                                     , 0.0);
  fNeutralHadronIso_DR0p1To0p2 = fmax(fmin((tmpNeutralHadronIso_DR0p1To0p2
                                            -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p1To0p2, 
                                                                   ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
                                           , 2.5)
                                       , 0.0);
  fNeutralHadronIso_DR0p2To0p3 = fmax(fmin((tmpNeutralHadronIso_DR0p2To0p3
                                          -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p2To0p3, 
                                                                 ele->SCluster()->Eta(),EffectiveAreaVersion))/ele->Pt()
                                         , 2.5)
                                     , 0.0);
  fNeutralHadronIso_DR0p3To0p4 = fmax(fmin((tmpNeutralHadronIso_DR0p3To0p4 
                                          -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p3To0p4, 
                                                                 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
                                         , 2.5)
                                     , 0.0);
  fNeutralHadronIso_DR0p4To0p5 = fmax(fmin((tmpNeutralHadronIso_DR0p4To0p5 
                                          -rho*eleT.ElectronEffectiveArea(eleT.kEleNeutralHadronIsoDR0p4To0p5, 
                                                                 ele->SCluster()->Eta(), EffectiveAreaVersion))/ele->Pt()
                                         , 2.5)
                                     , 0.0);

  if( ctrl.debug) {
    cout << "fNeutralHadronIso_DR0p0To0p1: " << fNeutralHadronIso_DR0p0To0p1 << endl;
    cout << "fNeutralHadronIso_DR0p1To0p2: " << fNeutralHadronIso_DR0p1To0p2 << endl;
    cout << "fNeutralHadronIso_DR0p2To0p3: " << fNeutralHadronIso_DR0p2To0p3 << endl;
    cout << "fNeutralHadronIso_DR0p3To0p4: " << fNeutralHadronIso_DR0p3To0p4 << endl;
    cout << "fNeutralHadronIso_DR0p4To0p5: " << fNeutralHadronIso_DR0p4To0p5 << endl;
  }

  double mvaval = eleIsoMVA->MVAValue_IsoRings( ele->Pt(),
                                                ele->SCluster()->Eta(),
                                                fChargedIso_DR0p0To0p1,
                                                fChargedIso_DR0p1To0p2,
                                                fChargedIso_DR0p2To0p3,
                                                fChargedIso_DR0p3To0p4,
                                                fChargedIso_DR0p4To0p5,
                                                fGammaIso_DR0p0To0p1,
                                                fGammaIso_DR0p1To0p2,
                                                fGammaIso_DR0p2To0p3,
                                                fGammaIso_DR0p3To0p4,
                                                fGammaIso_DR0p4To0p5,
                                                fNeutralHadronIso_DR0p0To0p1,
                                                fNeutralHadronIso_DR0p1To0p2,
                                                fNeutralHadronIso_DR0p2To0p3,
                                                fNeutralHadronIso_DR0p3To0p4,
                                                fNeutralHadronIso_DR0p4To0p5,
                                                ctrl.debug);

  SelectionStatus status;
  status.isoMVA = mvaval;
  bool pass = false;

  Int_t subdet = 0;
  if (fabs(ele->SCluster()->Eta()) < 0.8) subdet = 0;
  else if (fabs(ele->SCluster()->Eta()) < 1.479) subdet = 1;
  else subdet = 2;

  Int_t ptBin = 0;
  if (ele->Pt() >= 10.0) ptBin = 1;
  
  Int_t MVABin = -1;
  if (subdet == 0 && ptBin == 0) MVABin = 0;
  if (subdet == 1 && ptBin == 0) MVABin = 1;
  if (subdet == 2 && ptBin == 0) MVABin = 2;
  if (subdet == 0 && ptBin == 1) MVABin = 3;
  if (subdet == 1 && ptBin == 1) MVABin = 4;
  if (subdet == 2 && ptBin == 1) MVABin = 5;

  pass = false;
  if( MVABin == 0 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN0 ) pass = true;
  if( MVABin == 1 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN1 ) pass = true;
  if( MVABin == 2 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN2 ) pass = true;
  if( MVABin == 3 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN3 ) pass = true;
  if( MVABin == 4 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN4 ) pass = true;
  if( MVABin == 5 && mvaval > ELECTRON_LOOSE_ISOMVA_CUT_BIN5 ) pass = true;
  if( pass ) status.orStatus(SelectionStatus::LOOSEISO);

//   pass = false;
//   if( MVABin == 0 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN0 ) pass = true;
//   if( MVABin == 1 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN1 ) pass = true;
//   if( MVABin == 2 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN2 ) pass = true;
//   if( MVABin == 3 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN3 ) pass = true;
//   if( MVABin == 4 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN4 ) pass = true;
//   if( MVABin == 5 && mvaval > ELECTRON_TIGHT_ISOMVA_CUT_BIN5 ) pass = true;
//   if( pass ) status.orStatus(SelectionStatus::TIGHTISO);

  if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
  return status;
  
}


//--------------------------------------------------------------------------------------------------
void initElectronIsoMVA() { 
//--------------------------------------------------------------------------------------------------
  eleIsoMVA = new mithep::ElectronIDMVA();
  vector<string> weightFiles;
  weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_BarrelPt5To10.weights.xml");
  weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_EndcapPt5To10.weights.xml");
  weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_BarrelPt10ToInf.weights.xml");
  weightFiles.push_back("../MitPhysics/data/ElectronMVAWeights/ElectronIso_BDTG_V0_EndcapPt10ToInf.weights.xml");
  eleIsoMVA->Initialize( "ElectronIsoMVA",
                        mithep::ElectronIDMVA::kIsoRingsV0,
                        kTRUE, weightFiles);
}


//--------------------------------------------------------------------------------------------------
float electronPFIso04(ControlFlags &ctrl, 
                      const mithep::Electron * ele, 
                      const mithep::Vertex * vtx,
                      const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                      const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity, 
                      mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
                      vector<const mithep::Muon*> muonsToVeto,
                      vector<const mithep::Electron*> electronsToVeto) 
//--------------------------------------------------------------------------------------------------
{ 
  /*
  if( ctrl.debug ) { 
    cout << "electronIsoMVASelection :: muons to veto " << endl;
    for( int i=0; i<muonsToVeto.size(); i++ ) { 
      const mithep::Muon * vmu = muonsToVeto[i];
      cout << "\tpt: " << vmu->Pt() 
           << "\teta: " << vmu->Eta() 
           << "\tphi: " << vmu->Phi() 
           << endl;
    }
    cout << "electronIsoMVASelection :: electrons to veto " << endl;
    for( int i=0; i<electronsToVeto.size(); i++ ) { 
      const mithep::Electron * vel = electronsToVeto[i];
      cout << "\tpt: " << vel->Pt() 
           << "\teta: " << vel->Eta() 
           << "\tphi: " << vel->Phi() 
           << "\ttrk: " << vel->TrackerTrk()
           << endl;
    }
  }
  */

  //
  // final iso
  //
  Double_t fChargedIso = 0.0;
  Double_t fGammaIso = 0.0;
  Double_t fNeutralHadronIso = 0.0;


  //
  //Loop over PF Candidates
  //
  for(int k=0; k<fPFCandidates->GetEntries(); ++k) {


    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
    Double_t deta = (ele->Eta() - pf->Eta());
    Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
    Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());

    if (dr > 0.4) continue;
    if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack

    if(ctrl.debug) { 
      cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt() << "\tdR: " << dr;
      if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(*vtx)
                                     << "\ttrk: " << pf->HasTrackerTrk() 
                                     << "\tgsf: " << pf->HasGsfTrk(); 
      
      cout << endl;
    }


    //
    // sync : I don't think theyre doing this ...
    //
    //     if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
    //   (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) { 
    //       if( ctrl.debug ) cout << "\tskipping, matches to the electron ..."  << endl;
    //       continue;
    //     }


    //
    // Lepton Footprint Removal
    //
    Bool_t IsLeptonFootprint = kFALSE;
    if (dr < 1.0) {

      //
      // Check for electrons
      //
      for (Int_t q=0; q < electronsToVeto.size(); ++q) {
        const mithep::Electron *tmpele = electronsToVeto[q];
        /*
        // 4l electron
        if( pf->HasTrackerTrk()  ) { 
          if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
            if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
            IsLeptonFootprint = kTRUE;
          }
        }
        if( pf->HasGsfTrk()  ) { 
          if( pf->GsfTrk() == tmpele->GsfTrk() ) {
            if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
            IsLeptonFootprint = kTRUE;
          }
        }
        */
        // PF charged
        if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479 
            && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
          if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
          IsLeptonFootprint = kTRUE;
        }
        // PF gamma
        if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479 
            && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) {
          if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
          IsLeptonFootprint = kTRUE;
        }
      } // loop over electrons 

      /* KH - comment for sync            
      //
      // Check for muons
      //
      for (Int_t q=0; q < muonsToVeto.size(); ++q) {
        const mithep::Muon *tmpmu = muonsToVeto[q];
        // 4l muon
        if( pf->HasTrackerTrk() ) { 
          if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
            if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
            IsLeptonFootprint = kTRUE;
          }
        }
        // PF charged
        if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
          if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
          IsLeptonFootprint = kTRUE;
        }
      } // loop over muons
      */

    if (IsLeptonFootprint) 
      continue;

    //
    // Charged Iso
    //
    if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {

//       if( pf->HasTrackerTrk() )
//      if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
//       if( pf->HasGsfTrk() )
//      if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;

      // Veto any PFmuon, or PFEle
      if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) {
        if( ctrl.debug ) cout << "\t skipping, pf is and ele or mu .." <<endl;
        continue;
      }

      // Footprint Veto
      if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;

      if( ctrl.debug) cout << "charged:: pt: " << pf->Pt() 
                           << "\ttype: " << pf->PFType() 
                           << "\ttrk: " << pf->TrackerTrk() << endl;

      fChargedIso += pf->Pt();
    }

    //
    // Gamma Iso 
    //
    else if (abs(pf->PFType()) == PFCandidate::eGamma) {

      if (fabs(ele->SCluster()->Eta()) > 1.479) {
        if (mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta()) < 0.08) continue;
      }
      if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " " 
                           << dr << endl;
      // KH, add to sync 
      //      if( pf->Pt() > 0.5 ) 
        fGammaIso += pf->Pt();
    }

    //
    // Neutral Iso 
    //
    else {
      if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " " 
                           << dr << endl;
      // KH, add to sync 
      //      if( pf->Pt() > 0.5 ) 
        fNeutralHadronIso += pf->Pt();
    }

    }

  }


  double rho=0;
  if( (EffectiveAreaVersion == mithep::ElectronTools::kEleEAFall11MC) ||
      (EffectiveAreaVersion == mithep::ElectronTools::kEleEAData2011) ) { 
    if (!(isnan(fPUEnergyDensity->At(0)->RhoKt6PFJetsForIso25()) || 
          isinf(fPUEnergyDensity->At(0)->RhoKt6PFJetsForIso25()))) 
      rho = fPUEnergyDensity->At(0)->RhoKt6PFJetsForIso25();
    // !!!!!!!!!!!!! TMP HACK FOR SYNC !!!!!!!!!!!!!!!!!!!!!
    EffectiveAreaVersion  = mithep::ElectronTools::kEleEAData2011;
    // !!!!!!!!!!!!! TMP HACK FOR SYNC !!!!!!!!!!!!!!!!!!!!!
  } else { 
    if (!(isnan(fPUEnergyDensity->At(0)->RhoKt6PFJets()) || 
          isinf(fPUEnergyDensity->At(0)->RhoKt6PFJets()))) 
      rho = fPUEnergyDensity->At(0)->RhoKt6PFJets();
    // !!!!!!!!!!!!! TMP HACK FOR SYNC !!!!!!!!!!!!!!!!!!!!!
    EffectiveAreaVersion  = mithep::ElectronTools::kEleEAData2012;
    // !!!!!!!!!!!!! TMP HACK FOR SYNC !!!!!!!!!!!!!!!!!!!!!
  }
  if(ctrl.debug) cout << "rho: " << rho << endl;

  double pfIso = fChargedIso + fmax(0.0,(fGammaIso + fNeutralHadronIso 
                                        -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaAndNeutralHadronIso04,
                                                                   ele->Eta(),EffectiveAreaVersion)));


  gChargedIso = fChargedIso;
  gGammaIso = fGammaIso;
  gNeutralIso = fNeutralHadronIso;  
  return pfIso;
}


//--------------------------------------------------------------------------------------------------
// hacked version
float electronPFIso04(ControlFlags &ctrl, 
                      const mithep::Electron * ele, 
                      const mithep::Vertex * vtx,
                      const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                      float rho,
                      mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
                      vector<const mithep::Muon*> muonsToVeto,
                      vector<const mithep::Electron*> electronsToVeto) 
//--------------------------------------------------------------------------------------------------
{ 

  if( ctrl.debug ) { 
    cout << "electronIsoMVASelection :: muons to veto " << endl;
    for( int i=0; i<muonsToVeto.size(); i++ ) { 
      const mithep::Muon * vmu = muonsToVeto[i];
      cout << "\tpt: " << vmu->Pt() 
           << "\teta: " << vmu->Eta() 
           << "\tphi: " << vmu->Phi() 
           << endl;
    }
    cout << "electronIsoMVASelection :: electrons to veto " << endl;
    for( int i=0; i<electronsToVeto.size(); i++ ) { 
      const mithep::Electron * vel = electronsToVeto[i];
      cout << "\tpt: " << vel->Pt() 
           << "\teta: " << vel->Eta() 
           << "\tphi: " << vel->Phi() 
           << "\ttrk: " << vel->TrackerTrk()
           << endl;
    }
  }


  //
  // final iso
  //
  Double_t fChargedIso = 0.0;
  Double_t fGammaIso = 0.0;
  Double_t fNeutralHadronIso = 0.0;


  //
  //Loop over PF Candidates
  //
  for(int k=0; k<fPFCandidates->GetEntries(); ++k) {


    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
    Double_t deta = (ele->Eta() - pf->Eta());
    Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(ele->Phi()),Double_t(pf->Phi()));
    Double_t dr = mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta());

    if (dr > 0.4) continue;
    if( !(PFnoPUflag[k]) ) continue; // my PF no PU hack

    if(ctrl.debug) { 
      cout << "pf :: type: " << pf->PFType() << "\tpt: " << pf->Pt() << "\tdR: " << dr;
      if( pf->HasTrackerTrk() ) cout << "\tdZ: " << pf->TrackerTrk()->DzCorrected(*vtx)
                                     << "\ttrk: " << pf->HasTrackerTrk() 
                                     << "\tgsf: " << pf->HasGsfTrk(); 
      
      cout << endl;
    }


    //
    // sync : I don't think theyre doing this ...
    //
    //     if ( (pf->HasTrackerTrk() && (pf->TrackerTrk() == ele->TrackerTrk())) ||
    //   (pf->HasGsfTrk() && (pf->GsfTrk() == ele->GsfTrk()))) { 
    //       if( ctrl.debug ) cout << "\tskipping, matches to the electron ..."  << endl;
    //       continue;
    //     }


    //
    // Lepton Footprint Removal
    //
    Bool_t IsLeptonFootprint = kFALSE;
    if (dr < 1.0) {

      //
      // Check for electrons
      //
      for (Int_t q=0; q < electronsToVeto.size(); ++q) {
        const mithep::Electron *tmpele = electronsToVeto[q];
        /*
        // 4l electron
        if( pf->HasTrackerTrk()  ) { 
          if( pf->TrackerTrk() == tmpele->TrackerTrk() ) {
            if( ctrl.debug) cout << "\tcharged tktrk, matches 4L ele ..." << endl;
            IsLeptonFootprint = kTRUE;
          }
        }
        if( pf->HasGsfTrk()  ) { 
          if( pf->GsfTrk() == tmpele->GsfTrk() ) {
            if( ctrl.debug) cout << "\tcharged gsftrk, matches 4L ele ..." << endl;
            IsLeptonFootprint = kTRUE;
          }
        }
        */
        // PF charged
        if (pf->Charge() != 0 && fabs(tmpele->SCluster()->Eta()) > 1.479 
            && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.015) {
          if( ctrl.debug) cout << "\tcharged trk, dR matches 4L ele ..." << endl;
          IsLeptonFootprint = kTRUE;
        }
        // PF gamma
        if (abs(pf->PFType()) == PFCandidate::eGamma && fabs(tmpele->SCluster()->Eta()) > 1.479 
            && mithep::MathUtils::DeltaR(tmpele->Phi(),tmpele->Eta(), pf->Phi(), pf->Eta()) < 0.08) {
          if( ctrl.debug) cout << "\tPF gamma, matches 4L ele ..." << endl;
          IsLeptonFootprint = kTRUE;
        }
      } // loop over electrons 

      /* KH - comment for sync            
      //
      // Check for muons
      //
      for (Int_t q=0; q < muonsToVeto.size(); ++q) {
        const mithep::Muon *tmpmu = muonsToVeto[q];
        // 4l muon
        if( pf->HasTrackerTrk() ) { 
          if (pf->TrackerTrk() == tmpmu->TrackerTrk() ){
            if( ctrl.debug) cout << "\tmatches 4L mu ..." << endl;
            IsLeptonFootprint = kTRUE;
          }
        }
        // PF charged
        if (pf->Charge() != 0 && mithep::MathUtils::DeltaR(tmpmu->Phi(),tmpmu->Eta(), pf->Phi(), pf->Eta()) < 0.01) {
          if( ctrl.debug) cout << "\tcharged trk, dR matches 4L mu ..." << endl;
          IsLeptonFootprint = kTRUE;
        }
      } // loop over muons
      */

    if (IsLeptonFootprint) 
      continue;

    //
    // Charged Iso
    //
    if (pf->Charge() != 0 && (pf->HasTrackerTrk()||pf->HasGsfTrk()) ) {

//       if( pf->HasTrackerTrk() )
//      if (abs(pf->TrackerTrk()->DzCorrected(vtx)) > 0.2) continue;
//       if( pf->HasGsfTrk() )
//      if (abs(pf->GsfTrk()->DzCorrected(vtx)) > 0.2) continue;

      // Veto any PFmuon, or PFEle
      if (abs(pf->PFType()) == PFCandidate::eElectron || abs(pf->PFType()) == PFCandidate::eMuon) {
        if( ctrl.debug ) cout << "\t skipping, pf is and ele or mu .." <<endl;
        continue;
      }

      // Footprint Veto
      if (fabs(ele->SCluster()->Eta()) > 1.479 && dr < 0.015) continue;

      if( ctrl.debug) cout << "charged:: pt: " << pf->Pt() 
                           << "\ttype: " << pf->PFType() 
                           << "\ttrk: " << pf->TrackerTrk() << endl;

      fChargedIso += pf->Pt();
    }

    //
    // Gamma Iso 
    //
    else if (abs(pf->PFType()) == PFCandidate::eGamma) {

      if (fabs(ele->SCluster()->Eta()) > 1.479) {
        if (mithep::MathUtils::DeltaR(ele->Phi(),ele->Eta(), pf->Phi(), pf->Eta()) < 0.08) continue;
      }
      if( ctrl.debug) cout << "gamma:: " << pf->Pt() << " " 
                           << dr << endl;
      // KH, add to sync 
      //      if( pf->Pt() > 0.5 ) 
        fGammaIso += pf->Pt();
    }

    //
    // Neutral Iso 
    //
    else {
      if( ctrl.debug) cout << "neutral:: " << pf->Pt() << " " 
                           << dr << endl;
      // KH, add to sync 
      //      if( pf->Pt() > 0.5 ) 
        fNeutralHadronIso += pf->Pt();
    }

    }

  }

//   double rho = 0;
//   if (!(isnan(fPUEnergyDensity->At(0)->Rho()) || isinf(fPUEnergyDensity->At(0)->Rho()))) 
//     rho = fPUEnergyDensity->At(0)->Rho();

  // WARNING!!!!  
  // hardcode for sync ...
  EffectiveAreaVersion = eleT.kEleEAData2011;
  // WARNING!!!!  


  double pfIso = fChargedIso + fmax(0.0,(fGammaIso + fNeutralHadronIso 
                                        -rho*eleT.ElectronEffectiveArea(eleT.kEleGammaAndNeutralHadronIso04,
                                                                   ele->Eta(),EffectiveAreaVersion)));


  gChargedIso = fChargedIso;
  gGammaIso = fGammaIso;
  gNeutralIso = fNeutralHadronIso;  
  return pfIso;
}


//--------------------------------------------------------------------------------------------------
SelectionStatus electronReferenceIsoSelection(ControlFlags &ctrl, 
                                              const mithep::Electron * ele, 
                                              const mithep::Vertex * vtx,
                                              const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                                              const mithep::Array<mithep::PileupEnergyDensity> * fPUEnergyDensity, 
                                              mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
                                              vector<const mithep::Muon*> muonsToVeto,
                                              vector<const mithep::Electron*> electronsToVeto) 
//--------------------------------------------------------------------------------------------------
{

  SelectionStatus status;

  double pfIso = electronPFIso04( ctrl, ele, vtx, fPFCandidates, fPUEnergyDensity, 
                                  EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
  //  cout << "--------------> setting electron isoPF04 to " << pfIso << endl;
  status.isoPF04 = pfIso;
  status.chisoPF04 = gChargedIso;
  status.gaisoPF04 = gGammaIso;
  status.neisoPF04 = gNeutralIso;

  bool pass = false;
  if( (pfIso/ele->Pt()) < ELECTRON_REFERENCE_PFISO_CUT ) pass = true;

  if( pass ) {
    status.orStatus(SelectionStatus::LOOSEISO);
    status.orStatus(SelectionStatus::TIGHTISO);
  }
  if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
  return status;

}


//--------------------------------------------------------------------------------------------------
// hacked version
SelectionStatus electronReferenceIsoSelection(ControlFlags &ctrl, 
                                              const mithep::Electron * ele, 
                                              const mithep::Vertex * vtx,
                                              const mithep::Array<mithep::PFCandidate> * fPFCandidates,
                                              float rho,
                                              mithep::ElectronTools::EElectronEffectiveAreaTarget EffectiveAreaVersion,
                                              vector<const mithep::Muon*> muonsToVeto,
                                              vector<const mithep::Electron*> electronsToVeto) 
//--------------------------------------------------------------------------------------------------
{

  SelectionStatus status;

  double pfIso = electronPFIso04( ctrl, ele, vtx, fPFCandidates, rho, 
                                  EffectiveAreaVersion, muonsToVeto ,electronsToVeto );
  status.isoPF04 = pfIso;
  status.chisoPF04 = gChargedIso;
  status.gaisoPF04 = gGammaIso;
  status.neisoPF04 = gNeutralIso;
  bool pass = false;
  if( (pfIso/ele->Pt()) < ELECTRON_REFERENCE_PFISO_CUT ) pass = true;

  if( pass ) {
    status.orStatus(SelectionStatus::LOOSEISO);
    status.orStatus(SelectionStatus::TIGHTISO);
  }
  if(ctrl.debug) cout << "returning status : " << hex << status.getStatus() << dec << endl;
  return status;

}


//--------------------------------------------------------------------------------------------------
double  dbetaCorrectedIsoDr03(ControlFlags & ctrl, 
                              const mithep::PFCandidate * photon, 
                              const mithep::Muon * lepton, 
                              const mithep::Array<mithep::PFCandidate> * fPFCandidates)
//--------------------------------------------------------------------------------------------------
{ 

  //
  // final iso 
  //
  Double_t fChargedIso  = 0.0;
  Double_t fGammaIso  = 0.0;
  Double_t fNeutralHadronIso  = 0.0;
  Double_t fpfPU  = 0.0;

  //
  // Loop over PF Candidates
  //
  for(int k=0; k<fPFCandidates->GetEntries(); ++k) {

    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
    
    Double_t deta = (photon->Eta() - pf->Eta());
    Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(photon->Phi()),Double_t(pf->Phi()));
    Double_t dr = mithep::MathUtils::DeltaR(photon->Phi(),photon->Eta(), pf->Phi(), pf->Eta());
    if (dr > 0.3) continue;

    if( !(PFnoPUflag[k]) && pf->Charge() != 0 ) {
      if( pf->Pt() >= 0.2 && dr > 0.01 )
        fpfPU += pf->Pt();
      continue;
    }
    
    //
    // skip this photon
    //
    if( abs(pf->PFType()) == mithep::PFCandidate::eGamma && 
        pf->Et() == photon->Et() ) continue;
    
      
    //
    // Charged Iso 
    //
    if (pf->Charge() != 0 ) {
      if( dr > 0.01 && pf->Pt() >= 0.2 && 
          !(pf->TrackerTrk() == lepton->TrackerTrk()) )
        fChargedIso += pf->Pt();
    }
    
    //
    // Gamma Iso 
    //
    else if (abs(pf->PFType()) == mithep::PFCandidate::eGamma) {
      if( pf->Pt() > 0.5 && dr > 0.01) 
        fGammaIso += pf->Pt();
    }
    
    //
    // Other Neutrals
    //
    else {
      if( pf->Pt() > 0.5 && dr > 0.01) 
        fNeutralHadronIso += pf->Pt();
    }
    
  }
  
  if( ctrl.debug ) { 
    cout << "photon dbetaIso :: " << endl;
    cout << "\tfChargedIso: " << fChargedIso
         << "\tfGammaIso: " << fGammaIso
         << "\tfNeutralHadronIso: " << fNeutralHadronIso
         << "\tfpfPU: " << fpfPU 
         << endl;
  }
  double pfIso = fChargedIso + fGammaIso + fNeutralHadronIso - 0.5*fpfPU; 
  return pfIso/photon->Pt();
}


//--------------------------------------------------------------------------------------------------
double  dbetaCorrectedIsoDr03(ControlFlags & ctrl, 
                              const mithep::PFCandidate * photon, 
                              const mithep::Electron * lepton, 
                              const mithep::Array<mithep::PFCandidate> * fPFCandidates)
//--------------------------------------------------------------------------------------------------
{ 

  //
  // final iso 
  //
  Double_t fChargedIso  = 0.0;
  Double_t fGammaIso  = 0.0;
  Double_t fNeutralHadronIso  = 0.0;
  Double_t fpfPU  = 0.0;

  //
  // Loop over PF Candidates
  //
  for(int k=0; k<fPFCandidates->GetEntries(); ++k) {

    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
    
    Double_t deta = (photon->Eta() - pf->Eta());
    Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(photon->Phi()),Double_t(pf->Phi()));
    Double_t dr = mithep::MathUtils::DeltaR(photon->Phi(),photon->Eta(), pf->Phi(), pf->Eta());
    if (dr > 0.3) continue;

    if( !(PFnoPUflag[k]) && pf->Charge() != 0 ) {
      if( pf->Pt() >= 0.2 && dr > 0.01 )
        fpfPU += pf->Pt();
      continue;
    }
    
    //
    // skip this photon
    //
    if( abs(pf->PFType()) == mithep::PFCandidate::eGamma && 
        pf->Et() == photon->Et() ) continue;
    
      
    //
    // Charged Iso 
    //
    if (pf->Charge() != 0 ) {
      if( dr > 0.01 && pf->Pt() >= 0.2 && 
          !(pf->TrackerTrk() == lepton->TrackerTrk()) )
        fChargedIso += pf->Pt();
    }
    
    //
    // Gamma Iso 
    //
    else if (abs(pf->PFType()) == mithep::PFCandidate::eGamma) {
      if( pf->Pt() > 0.5 && dr > 0.01) 
        fGammaIso += pf->Pt();
    }
    
    //
    // Other Neutrals
    //
    else {
      if( pf->Pt() > 0.5 && dr > 0.01) 
        fNeutralHadronIso += pf->Pt();
    }
    
  }
  
  if( ctrl.debug ) { 
    cout << "photon dbetaIso :: " << endl;
    cout << "\tfChargedIso: " << fChargedIso
         << "\tfGammaIso: " << fGammaIso
         << "\tfNeutralHadronIso: " << fNeutralHadronIso
         << "\tfpfPU: " << fpfPU 
         << endl;
  }
  double pfIso = fChargedIso + fGammaIso + fNeutralHadronIso - 0.5*fpfPU; 
  return pfIso/photon->Pt();
}


//--------------------------------------------------------------------------------------------------
double  betaCorrectedIsoDr03(ControlFlags & ctrl, 
                             const mithep::PFCandidate * photon, 
                             const mithep::Muon * lepton, 
                             const mithep::Array<mithep::PFCandidate> * fPFCandidates)
//--------------------------------------------------------------------------------------------------
{ 

  //
  // final iso 
  //
  Double_t fChargedIso  = 0.0;
  Double_t fGammaIso  = 0.0;
  Double_t fNeutralHadronIso  = 0.0;
  Double_t fpfPU  = 0.0;

  //
  // Loop over PF Candidates
  //
  for(int k=0; k<fPFCandidates->GetEntries(); ++k) {

    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
    
    Double_t deta = (photon->Eta() - pf->Eta());
    Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(photon->Phi()),Double_t(pf->Phi()));
    Double_t dr = mithep::MathUtils::DeltaR(photon->Phi(),photon->Eta(), pf->Phi(), pf->Eta());
    if (dr > 0.3) continue;

    if( !(PFnoPUflag[k]) && pf->Charge() != 0 ) {
      if( pf->Pt() >= 0.2 && dr > 0.01 )
        fpfPU += pf->Pt();
      continue;
    }
    
    //
    // skip this photon
    //
    if( abs(pf->PFType()) == mithep::PFCandidate::eGamma && 
        pf->Et() == photon->Et() ) continue;
    
      
    //
    // Charged Iso 
    //
    if (pf->Charge() != 0 ) {
      if( dr > 0.01 && pf->Pt() >= 0.2 && 
          !(pf->TrackerTrk() == lepton->TrackerTrk()) )
        fChargedIso += pf->Pt();
    }
    
    //
    // Gamma Iso 
    //
    else if (abs(pf->PFType()) == mithep::PFCandidate::eGamma) {
      if( pf->Pt() > 0.5 && dr > 0.01) 
        fGammaIso += pf->Pt();
    }
    
    //
    // Other Neutrals
    //
    else {
      if( pf->Pt() > 0.5 && dr > 0.01) 
        fNeutralHadronIso += pf->Pt();
    }
    
  }
  
  if( ctrl.debug ) { 
    cout << "photon dbetaIso :: " << endl;
    cout << "\tfChargedIso: " << fChargedIso
         << "\tfGammaIso: " << fGammaIso
         << "\tfNeutralHadronIso: " << fNeutralHadronIso
         << "\tfpfPU: " << fpfPU 
         << endl;
  }
  double pfIso = fChargedIso + fGammaIso + fNeutralHadronIso + fpfPU; 
  return pfIso/photon->Pt();
}


//--------------------------------------------------------------------------------------------------
double  betaCorrectedIsoDr03(ControlFlags & ctrl, 
                              const mithep::PFCandidate * photon, 
                              const mithep::Electron * lepton, 
                              const mithep::Array<mithep::PFCandidate> * fPFCandidates)
//--------------------------------------------------------------------------------------------------
{ 

  //
  // final iso 
  //
  Double_t fChargedIso  = 0.0;
  Double_t fGammaIso  = 0.0;
  Double_t fNeutralHadronIso  = 0.0;
  Double_t fpfPU  = 0.0;

  //
  // Loop over PF Candidates
  //
  for(int k=0; k<fPFCandidates->GetEntries(); ++k) {

    const mithep::PFCandidate *pf = (mithep::PFCandidate*)((*fPFCandidates)[k]);
    
    Double_t deta = (photon->Eta() - pf->Eta());
    Double_t dphi = mithep::MathUtils::DeltaPhi(Double_t(photon->Phi()),Double_t(pf->Phi()));
    Double_t dr = mithep::MathUtils::DeltaR(photon->Phi(),photon->Eta(), pf->Phi(), pf->Eta());
    if (dr > 0.3) continue;

    if( !(PFnoPUflag[k]) && pf->Charge() != 0 ) {
      if( pf->Pt() >= 0.2 && dr > 0.01 )
        fpfPU += pf->Pt();
      continue;
    }
    
    //
    // skip this photon
    //
    if( abs(pf->PFType()) == mithep::PFCandidate::eGamma && 
        pf->Et() == photon->Et() ) continue;
    
      
    //
    // Charged Iso 
    //
    if (pf->Charge() != 0 ) {
      if( dr > 0.01 && pf->Pt() >= 0.2 && 
          !(pf->TrackerTrk() == lepton->TrackerTrk()) )
        fChargedIso += pf->Pt();
    }
    
    //
    // Gamma Iso 
    //
    else if (abs(pf->PFType()) == mithep::PFCandidate::eGamma) {
      if( pf->Pt() > 0.5 && dr > 0.01) 
        fGammaIso += pf->Pt();
    }
    
    //
    // Other Neutrals
    //
    else {
      if( pf->Pt() > 0.5 && dr > 0.01) 
        fNeutralHadronIso += pf->Pt();
    }
    
  }
  
  if( ctrl.debug ) { 
    cout << "photon dbetaIso :: " << endl;
    cout << "\tfChargedIso: " << fChargedIso
         << "\tfGammaIso: " << fGammaIso
         << "\tfNeutralHadronIso: " << fNeutralHadronIso
         << "\tfpfPU: " << fpfPU 
         << endl;
  }
  double pfIso = fChargedIso + fGammaIso + fNeutralHadronIso + fpfPU; 
  return pfIso/photon->Pt();
}


Revision:	1.29
Committed:	Sun Jun 3 15:51:51 2012 UTC (12 years, 11 months ago) by khahn
Content type:	text/plain
Branch:	MAIN
Changes since 1.28:	+332 -0 lines
Log Message:	FSR syncing...