PFPhotonIsolationCalculator.cc

Go to the documentation of this file.

#include "RecoEgamma/PhotonIdentification/interface/PFPhotonIsolationCalculator.h"
#include <cmath>
#include "DataFormats/Math/interface/deltaR.h"



#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
#include "DataFormats/Common/interface/RefToPtr.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/IPTools/interface/IPTools.h"



void PFPhotonIsolationCalculator::setup(const edm::ParameterSet& conf) {
 

 iParticleType_ = conf.getParameter<int>("particleType");
  if (  iParticleType_ ==1 ) { 

    fConeSize_     = conf.getParameter<double>("coneDR");
    iNumberOfRings_ = conf.getParameter<int>("numberOfRings");
    fRingSize_     = conf.getParameter<double>("ringSize");
    //
    bApplyVeto_    = conf.getParameter<bool>("applyVeto");
    bApplyPFPUVeto_    = conf.getParameter<bool>("applyPFPUVeto");
    bApplyDzDxyVeto_    = conf.getParameter<bool>("applyDzDxyVeto");
    bApplyMissHitPhVeto_    = conf.getParameter<bool>("applyMissHitPhVeto");
    bDeltaRVetoBarrel_    = conf.getParameter<bool>("deltaRVetoBarrel");
    bDeltaRVetoEndcap_    = conf.getParameter<bool>("deltaRVetoEndcap");
    bRectangleVetoBarrel_    = conf.getParameter<bool>("rectangleVetoBarrel");
    bRectangleVetoEndcap_    = conf.getParameter<bool>("rectangleVetoEndcap");
    bUseCrystalSize_    = conf.getParameter<bool>("useCrystalSize");
    //
    fDeltaRVetoBarrelPhotons_    = conf.getParameter<double>("deltaRVetoBarrelPhotons");   
    fDeltaRVetoBarrelNeutrals_    = conf.getParameter<double>("deltaRVetoBarrelNeutrals");   
    fDeltaRVetoBarrelCharged_    = conf.getParameter<double>("deltaRVetoBarrelCharged");   
    fDeltaRVetoEndcapPhotons_    = conf.getParameter<double>("deltaRVetoEndcapPhotons");   
    fDeltaRVetoEndcapNeutrals_    = conf.getParameter<double>("deltaRVetoEndcapNeutrals");   
    fDeltaRVetoEndcapCharged_    = conf.getParameter<double>("deltaRVetoEndcapCharged");   
    fNumberOfCrystalEndcapPhotons_ = conf.getParameter<double>("numberOfCrystalEndcapPhotons");  
    //
    fRectangleDeltaPhiVetoBarrelPhotons_    = conf.getParameter<double>("rectangleDeltaPhiVetoBarrelPhotons");
    fRectangleDeltaPhiVetoBarrelNeutrals_   = conf.getParameter<double>("rectangleDeltaPhiVetoBarrelNeutrals");
    fRectangleDeltaPhiVetoBarrelCharged_    = conf.getParameter<double>("rectangleDeltaPhiVetoBarrelCharged");
    fRectangleDeltaPhiVetoEndcapPhotons_    = conf.getParameter<double>("rectangleDeltaPhiVetoEndcapPhotons");
    fRectangleDeltaPhiVetoEndcapNeutrals_   = conf.getParameter<double>("rectangleDeltaPhiVetoEndcapNeutrals");
    fRectangleDeltaPhiVetoEndcapCharged_    = conf.getParameter<double>("rectangleDeltaPhiVetoEndcapCharged");
    //
    fRectangleDeltaEtaVetoBarrelPhotons_    = conf.getParameter<double>("rectangleDeltaEtaVetoBarrelPhotons");
    fRectangleDeltaEtaVetoBarrelNeutrals_   = conf.getParameter<double>("rectangleDeltaEtaVetoBarrelNeutrals");
    fRectangleDeltaEtaVetoBarrelCharged_    = conf.getParameter<double>("rectangleDeltaEtaVetoBarrelCharged");
    fRectangleDeltaEtaVetoEndcapPhotons_    = conf.getParameter<double>("rectangleDeltaEtaVetoEndcapPhotons");
    fRectangleDeltaEtaVetoEndcapNeutrals_   = conf.getParameter<double>("rectangleDeltaEtaVetoEndcapNeutrals");
    fRectangleDeltaEtaVetoEndcapCharged_    = conf.getParameter<double>("rectangleDeltaEtaVetoEndcapCharged");
    //
    bCheckClosestZVertex_    = conf.getParameter<bool>("checkClosestZVertex");
    initializeRings(iNumberOfRings_, fConeSize_);
    
  }  
 



}



//--------------------------------------------------------------------------------------------------

PFPhotonIsolationCalculator::PFPhotonIsolationCalculator() {
  // Constructor.
}



//--------------------------------------------------------------------------------------------------
PFPhotonIsolationCalculator::~PFPhotonIsolationCalculator()
{

}


void PFPhotonIsolationCalculator::calculate(const reco::Photon* aPho, 
                        const edm::Handle<reco::PFCandidateCollection> pfCandidateHandle,
                        //                      reco::VertexRef vtx,
                        edm::Handle< reco::VertexCollection >& vertices,
                        const edm::Event& e , const edm::EventSetup& es,
                        reco::Photon::PflowIsolationVariables& phoisol03) {

  reco::VertexRef vtx(vertices, 0);
  this->fGetIsolation(&*aPho,  pfCandidateHandle, vtx, vertices);
  phoisol03.chargedHadronIso = this->getIsolationCharged() ;
  phoisol03.neutralHadronIso = this->getIsolationNeutral();
  phoisol03.photonIso        = this->getIsolationPhoton();


  
 }




//--------------------------------------------------------------------------------------------------
void PFPhotonIsolationCalculator::initializeRings(int iNumberOfRings, float fRingSize){
 
  fIsolationInRings_.clear();
  for(int isoBin =0;isoBin<iNumberOfRings;isoBin++){
    float fTemp = 0.0;
    fIsolationInRings_.push_back(fTemp);
    
    float fTempPhoton = 0.0;
    fIsolationInRingsPhoton_.push_back(fTempPhoton);

    float fTempNeutral = 0.0;
    fIsolationInRingsNeutral_.push_back(fTempNeutral);

    float fTempCharged = 0.0;
    fIsolationInRingsCharged_.push_back(fTempCharged);

    float fTempChargedAll = 0.0;
    fIsolationInRingsChargedAll_.push_back(fTempChargedAll);

  }

  fConeSize_ = fRingSize * (float)iNumberOfRings;

}
  
 


//--------------------------------------------------------------------------------------------------
float PFPhotonIsolationCalculator::fGetIsolation(const reco::Photon * photon, const  edm::Handle<reco::PFCandidateCollection> pfCandidateHandle ,reco::VertexRef vtx, edm::Handle< reco::VertexCollection > vertices) {
 
  fGetIsolationInRings( photon, pfCandidateHandle, vtx, vertices);
  fIsolation_ = fIsolationInRings_[0];
  return fIsolation_;
}



//--------------------------------------------------------------------------------------------------
std::vector<float > PFPhotonIsolationCalculator::fGetIsolationInRings(const reco::Photon * photon, const  edm::Handle<reco::PFCandidateCollection> pfCandidateHandle, reco::VertexRef vtx, edm::Handle< reco::VertexCollection > vertices) {


  int isoBin;
  
  for(isoBin =0;isoBin<iNumberOfRings_;isoBin++){
    fIsolationInRings_[isoBin]=0.;
    fIsolationInRingsPhoton_[isoBin]=0.;
    fIsolationInRingsNeutral_[isoBin]=0.;
    fIsolationInRingsCharged_[isoBin]=0.;
    fIsolationInRingsChargedAll_[isoBin]=0.;
  }
  
  iMissHits_ = 0;

  refSC = photon->superCluster();
  pivotInBarrel = std::fabs((refSC->position().eta()))<1.479;

  for(unsigned iPF=0; iPF< pfCandidateHandle->size(); iPF++) {

    reco::PFCandidateRef pfParticle(reco::PFCandidateRef(pfCandidateHandle, iPF));

    if (pfParticle->superClusterRef().isNonnull() &&
        photon->superCluster().isNonnull() &&
        pfParticle->superClusterRef() == photon->superCluster()) 
      continue;
    

  
    if(pfParticle->pdgId()==22){
      // Set the vertex of reco::Photon to the first PV
      math::XYZVector direction = math::XYZVector(photon->superCluster()->x() - pfParticle->vx(),
                          photon->superCluster()->y() - pfParticle->vy(),
                                                  photon->superCluster()->z() - pfParticle->vz());
      
      fEta_ = direction.Eta();
      fPhi_ = direction.Phi();
      fVx_ = pfParticle->vx();
      fVy_ = pfParticle->vy();
      fVz_ = pfParticle->vz();

      float fDeltaR =isPhotonParticleVetoed(pfParticle); 
      if( fDeltaR >=0.){
        isoBin = (int)(fDeltaR/fRingSize_);
        fIsolationInRingsPhoton_[isoBin] = fIsolationInRingsPhoton_[isoBin] + pfParticle->pt();
      }
      
    }else if(std::abs(pfParticle->pdgId())==130){
      // Set the vertex of reco::Photon to the first PV
      math::XYZVector direction = math::XYZVector(photon->superCluster()->x() - pfParticle->vx(),
                                                  photon->superCluster()->y() - pfParticle->vy(),
                                                  photon->superCluster()->z() - pfParticle->vz());

      fEta_ = direction.Eta();
      fPhi_ = direction.Phi();
      fVx_ = pfParticle->vx();
      fVy_ = pfParticle->vy();
      fVz_ = pfParticle->vz();
      float fDeltaR =  isNeutralParticleVetoed( pfParticle); 
      if( fDeltaR>=0.){
    isoBin = (int)(fDeltaR/fRingSize_);
        fIsolationInRingsNeutral_[isoBin] = fIsolationInRingsNeutral_[isoBin] + pfParticle->pt();
      }

      //}else if(abs(pfParticle.pdgId()) == 11 ||abs(pfParticle.pdgId()) == 13 || abs(pfParticle.pdgId()) == 211){
    }else if(std::abs(pfParticle->pdgId()) == 211){
      // Set the vertex of reco::Photon to the first PV
      math::XYZVector direction = math::XYZVector(photon->superCluster()->x() - (*vtx).x(),
                                                  photon->superCluster()->y() - (*vtx).y(),
                                                  photon->superCluster()->z() - (*vtx).z());

      fEta_ = direction.Eta();
      fPhi_ = direction.Phi();
      fVx_ = (*vtx).x();
      fVy_ = (*vtx).y();
      fVz_ = (*vtx).z();
      float fDeltaR = isChargedParticleVetoed( pfParticle, vtx, vertices);
      if( fDeltaR >=0.){
        isoBin = (int)(fDeltaR/fRingSize_);
        fIsolationInRingsCharged_[isoBin] = fIsolationInRingsCharged_[isoBin] + pfParticle->pt();
      }

    }
  }

 
  for(int isoBin =0;isoBin<iNumberOfRings_;isoBin++){
    fIsolationInRings_[isoBin]= fIsolationInRingsPhoton_[isoBin]+ fIsolationInRingsNeutral_[isoBin] + fIsolationInRingsCharged_[isoBin];
  }
  
  return fIsolationInRings_;
}




//--------------------------------------------------------------------------------------------------
float PFPhotonIsolationCalculator::isPhotonParticleVetoed( const reco::PFCandidate* pfIsoCand ){
  
  
  float fDeltaR = deltaR(fEta_,fPhi_,pfIsoCand->eta(),pfIsoCand->phi());

  if(fDeltaR > fConeSize_)
    return -999.;
  
  float fDeltaPhi = deltaPhi(fPhi_,pfIsoCand->phi());
  float fDeltaEta = fEta_-pfIsoCand->eta();

  if(!bApplyVeto_)
    return fDeltaR;
 
  //NOTE: get the direction for the EB/EE transition region from the deposit just to be in synch with the isoDep
  // this will be changed in the future
  
  if(bApplyMissHitPhVeto_) {
    if(iMissHits_ > 0)
      if(pfIsoCand->mva_nothing_gamma() > 0.99) {
        if(pfIsoCand->superClusterRef().isNonnull() && refSC.isNonnull()) {
         if(pfIsoCand->superClusterRef() == refSC)
         return -999.;
        }
      }
  }

  if(pivotInBarrel){
    if(bDeltaRVetoBarrel_){
      if(fDeltaR < fDeltaRVetoBarrelPhotons_)
        return -999.;
    }
    
    if(bRectangleVetoBarrel_){
      if(std::abs(fDeltaEta) < fRectangleDeltaEtaVetoBarrelPhotons_ && std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoBarrelPhotons_ ){
        return -999.;
      }
    }
  }else{
    if (bUseCrystalSize_ == true) {
      fDeltaRVetoEndcapPhotons_ = 0.00864*std::fabs(refSC->position().z()/sqrt(refSC->position().perp2()))*fNumberOfCrystalEndcapPhotons_;
    }

    if(bDeltaRVetoEndcap_){
      if(fDeltaR < fDeltaRVetoEndcapPhotons_)
        return -999.;
    }
    if(bRectangleVetoEndcap_){
      if(std::abs(fDeltaEta) < fRectangleDeltaEtaVetoEndcapPhotons_ && std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoEndcapPhotons_ ){
         return -999.;
      }
    }
  }

  return fDeltaR;
}



//--------------------------------------------------------------------------------------------------
float PFPhotonIsolationCalculator::isPhotonParticleVetoed( const reco::PFCandidateRef pfIsoCand ){
  
  
  float fDeltaR = deltaR(fEta_,fPhi_,pfIsoCand->eta(),pfIsoCand->phi());

  if(fDeltaR > fConeSize_)
    return -999.;
  
  float fDeltaPhi = deltaPhi(fPhi_,pfIsoCand->phi());
  float fDeltaEta = fEta_-pfIsoCand->eta();

  if(!bApplyVeto_)
    return fDeltaR;
 
  //NOTE: get the direction for the EB/EE transition region from the deposit just to be in synch with the isoDep
  // this will be changed in the future
  
  if(bApplyMissHitPhVeto_) {
    if(iMissHits_ > 0)
      if(pfIsoCand->mva_nothing_gamma() > 0.99) {
        if(pfIsoCand->superClusterRef().isNonnull() && refSC.isNonnull()) {
         if(pfIsoCand->superClusterRef() == refSC)
         return -999.;
        }
      }
  }

  if(pivotInBarrel){
    if(bDeltaRVetoBarrel_){
      if(fDeltaR < fDeltaRVetoBarrelPhotons_)
        return -999.;
    }
    
    if(bRectangleVetoBarrel_){
      if(std::abs(fDeltaEta) < fRectangleDeltaEtaVetoBarrelPhotons_ && std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoBarrelPhotons_){
        return -999.;
      }
    }
  }else{
    if (bUseCrystalSize_ == true) {
      fDeltaRVetoEndcapPhotons_ = 0.00864*std::fabs(refSC->position().z()/sqrt(refSC->position().perp2()))*fNumberOfCrystalEndcapPhotons_;
    }

    if(bDeltaRVetoEndcap_){
      if(fDeltaR < fDeltaRVetoEndcapPhotons_)
        return -999.;
    }
    if(bRectangleVetoEndcap_){
      if(std::abs(fDeltaEta) < fRectangleDeltaEtaVetoEndcapPhotons_ && std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoEndcapPhotons_){
         return -999.;
      }
    }
  }

  return fDeltaR;
}


//--------------------------------------------------------------------------------------------------
float PFPhotonIsolationCalculator::isNeutralParticleVetoed( const reco::PFCandidate* pfIsoCand ){

  float fDeltaR = deltaR(fEta_,fPhi_,pfIsoCand->eta(),pfIsoCand->phi());
  
  if(fDeltaR > fConeSize_)
    return -999;
  
  float fDeltaPhi = deltaPhi(fPhi_,pfIsoCand->phi());
  float fDeltaEta = fEta_-pfIsoCand->eta();

  if(!bApplyVeto_)
    return fDeltaR;

  //NOTE: get the direction for the EB/EE transition region from the deposit just to be in synch with the isoDep
  // this will be changed in the future
  if(pivotInBarrel){
    if(!bDeltaRVetoBarrel_ && !bRectangleVetoBarrel_){
      return fDeltaR;
    }
    
    if(bDeltaRVetoBarrel_ ){
        if(fDeltaR < fDeltaRVetoBarrelNeutrals_)
         return -999.;
      }
      if(bRectangleVetoBarrel_){
        if(std::abs(fDeltaEta) < fRectangleDeltaEtaVetoBarrelNeutrals_ && std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoBarrelNeutrals_){
         return -999.;
        }
      }
      
    }else{
     if(!bDeltaRVetoEndcap_  &&!  bRectangleVetoEndcap_){
       return fDeltaR;
     }
      if(bDeltaRVetoEndcap_){
        if(fDeltaR < fDeltaRVetoEndcapNeutrals_)
         return -999.;
      }
      if(bRectangleVetoEndcap_){
        if(std::abs(fDeltaEta) < fRectangleDeltaEtaVetoEndcapNeutrals_ && std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoEndcapNeutrals_){
         return -999.;
        }
      }
  }

  return fDeltaR;
}


//--------------------------------------------------------------------------------------------------
float PFPhotonIsolationCalculator::isNeutralParticleVetoed( const reco::PFCandidateRef pfIsoCand ){

  float fDeltaR = deltaR(fEta_,fPhi_,pfIsoCand->eta(),pfIsoCand->phi());
  
  if(fDeltaR > fConeSize_)
    return -999;
  
  float fDeltaPhi = deltaPhi(fPhi_,pfIsoCand->phi());
  float fDeltaEta = fEta_-pfIsoCand->eta();

  if(!bApplyVeto_)
    return fDeltaR;

  //NOTE: get the direction for the EB/EE transition region from the deposit just to be in synch with the isoDep
  // this will be changed in the future
  if(pivotInBarrel){
    if(!bDeltaRVetoBarrel_  && !bRectangleVetoBarrel_ ){
      return fDeltaR;
    }
    
    if(bDeltaRVetoBarrel_){
        if(fDeltaR < fDeltaRVetoBarrelNeutrals_)
         return -999.;
      }
      if(bRectangleVetoBarrel_){
        if(std::abs(fDeltaEta) < fRectangleDeltaEtaVetoBarrelNeutrals_ && std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoBarrelNeutrals_){
         return -999.;
        }
      }
      
    }else{
     if(!bDeltaRVetoEndcap_ && !bRectangleVetoEndcap_ ){
       return fDeltaR;
     }
      if(bDeltaRVetoEndcap_ ){
        if(fDeltaR < fDeltaRVetoEndcapNeutrals_ )
         return -999.;
      }
      if(bRectangleVetoEndcap_){
        if(std::abs(fDeltaEta) < fRectangleDeltaEtaVetoEndcapNeutrals_ && std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoEndcapNeutrals_){
         return -999.;
        }
      }
  }

  return fDeltaR;
}




//----------------------------------------------------------------------------------------------------
float PFPhotonIsolationCalculator::isChargedParticleVetoed(const reco::PFCandidate* pfIsoCand, edm::Handle< reco::VertexCollection > vertices ){
  //need code to handle special conditions
  
  return -999;
}

//-----------------------------------------------------------------------------------------------------
float PFPhotonIsolationCalculator::isChargedParticleVetoed(const reco::PFCandidate* pfIsoCand,reco::VertexRef vtxMain, edm::Handle< reco::VertexCollection > vertices ){
  
  reco::VertexRef vtx = chargedHadronVertex(vertices, *pfIsoCand );
  if(vtx.isNull())
    return -999.;
  
// float fVtxMainX = (*vtxMain).x();
// float fVtxMainY = (*vtxMain).y();
  float fVtxMainZ = (*vtxMain).z();

  if(bApplyPFPUVeto_) {
    if(vtx != vtxMain)
      return -999.;
  }
    

  if(bApplyDzDxyVeto_) {
    if(iParticleType_==kPhoton){
      
      float dz = std::fabs( pfIsoCand->trackRef()->dz( (*vtxMain).position() ) );
      if (dz > 0.2)
        return -999.;
        
      double dxy = pfIsoCand->trackRef()->dxy( (*vtxMain).position() );
      if (std::fabs(dxy) > 0.1)
        return -999.;
      
      /*
float dz = fabs(vtx->z() - fVtxMainZ);
if (dz > 1.)
        return -999.;
double dxy = ( -(vtx->x() - fVtxMainX)*pfIsoCand->py() + (vtx->y() - fVtxMainY)*pfIsoCand->px()) / pfIsoCand->pt();
if(fabs(dxy) > 0.2)
        return -999.;
*/
    }else{
      
      
      float dz = std::fabs(vtx->z() - fVtxMainZ);
      if (dz > 1.)
        return -999.;
      
      double dxy = ( -(vtx->x() - fVx_)*pfIsoCand->py() + (vtx->y() - fVy_)*pfIsoCand->px()) / pfIsoCand->pt();
      if(std::fabs(dxy) > 0.1)
        return -999.;
    }
  }

  float fDeltaR = deltaR(pfIsoCand->eta(),pfIsoCand->phi(),fEta_,fPhi_);

  if(fDeltaR > fConeSize_)
    return -999.;

  float fDeltaPhi = deltaPhi(fPhi_,pfIsoCand->phi());
  float fDeltaEta = fEta_-pfIsoCand->eta();
  

// std::abscout << " charged hadron: DR " << fDeltaR
//          << " pt " << pfIsoCand->pt() << " eta,phi " << pfIsoCand->eta() << ", " << pfIsoCand->phi()
//          << " fVtxMainZ " << (*vtxMain).z() << " cand z " << vtx->z() << std::endl;
  

  if(!bApplyVeto_)
    return fDeltaR;
  
  //NOTE: get the direction for the EB/EE transition region from the deposit just to be in synch with the isoDep
  // this will be changed in the future
  if(pivotInBarrel){
    if(!bDeltaRVetoBarrel_ &&!bRectangleVetoBarrel_ ){
      return fDeltaR;
    }
    
    if(bDeltaRVetoBarrel_){
        if(fDeltaR < fDeltaRVetoBarrelCharged_)
         return -999.;
      }
      if(bRectangleVetoBarrel_){
        if(std::abs(fDeltaEta) < fRectangleDeltaEtaVetoBarrelCharged_ && std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoBarrelCharged_){
         return -999.;
        }
      }
      
    }else{
     if(!bDeltaRVetoEndcap_  && !bRectangleVetoEndcap_){
       return fDeltaR;
     }
      if(bDeltaRVetoEndcap_){
        if(fDeltaR < fDeltaRVetoEndcapCharged_)
         return -999.;
      }
      if(bRectangleVetoEndcap_){
        if(std::abs(fDeltaEta) < fRectangleDeltaEtaVetoEndcapCharged_ && std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoEndcapCharged_){
         return -999.;
        }
      }
  }
                
  


  return fDeltaR;
}


//-----------------------------------------------------------------------------------------------------
float PFPhotonIsolationCalculator::isChargedParticleVetoed(const reco::PFCandidateRef pfIsoCand,reco::VertexRef vtxMain, edm::Handle< reco::VertexCollection > vertices ){
  
  reco::VertexRef vtx = chargedHadronVertex(vertices, *pfIsoCand );
  if(vtx.isNull())
    return -999.;
  
// float fVtxMainX = (*vtxMain).x();
// float fVtxMainY = (*vtxMain).y();
  float fVtxMainZ = (*vtxMain).z();

  if(bApplyPFPUVeto_) {
    if(vtx != vtxMain)
      return -999.;
  }
    

  if(bApplyDzDxyVeto_) {
    if(iParticleType_==kPhoton){
      
      float dz = std::fabs( pfIsoCand->trackRef()->dz( (*vtxMain).position() ) );
      if (dz > 0.2)
        return -999.;
        
      double dxy = pfIsoCand->trackRef()->dxy( (*vtxMain).position() );
      if (std::fabs(dxy) > 0.1)
        return -999.;
      
      /*
float dz = fabs(vtx->z() - fVtxMainZ);
if (dz > 1.)
        return -999.;
double dxy = ( -(vtx->x() - fVtxMainX)*pfIsoCand->py() + (vtx->y() - fVtxMainY)*pfIsoCand->px()) / pfIsoCand->pt();
if(fabs(dxy) > 0.2)
        return -999.;
*/
    }else{
      
      
      float dz = std::fabs(vtx->z() - fVtxMainZ);
      if (dz > 1.)
        return -999.;
      
      double dxy = ( -(vtx->x() - fVx_)*pfIsoCand->py() + (vtx->y() - fVy_)*pfIsoCand->px()) / pfIsoCand->pt();
      if(std::fabs(dxy) > 0.1)
        return -999.;
    }
  }

  float fDeltaR = deltaR(pfIsoCand->eta(),pfIsoCand->phi(),fEta_,fPhi_);

  if(fDeltaR > fConeSize_)
    return -999.;

  float fDeltaPhi = deltaPhi(fPhi_,pfIsoCand->phi());
  float fDeltaEta = fEta_-pfIsoCand->eta();
  

// std::cout << " charged hadron: DR " << fDeltaR
//          << " pt " << pfIsoCand->pt() << " eta,phi " << pfIsoCand->eta() << ", " << pfIsoCand->phi()
//          << " fVtxMainZ " << (*vtxMain).z() << " cand z " << vtx->z() << std::endl;
  

  if(!bApplyVeto_)
    return fDeltaR;
  
  //NOTE: get the direction for the EB/EE transition region from the deposit just to be in synch with the isoDep
  // this will be changed in the future
  if(pivotInBarrel){
    if(!bDeltaRVetoBarrel_ && !bRectangleVetoBarrel_) {
      return fDeltaR;
    }
    
    if(bDeltaRVetoBarrel_){
        if(fDeltaR < fDeltaRVetoBarrelCharged_)
         return -999.;
      }
      if(bRectangleVetoBarrel_){
        if(std::abs(fDeltaEta) < fRectangleDeltaEtaVetoBarrelCharged_ && std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoBarrelCharged_){
         return -999.;
        }
      }
      
    }else{
     if(!bDeltaRVetoEndcap_ &&!bRectangleVetoEndcap_ ){
       return fDeltaR;
     }
      if(bDeltaRVetoEndcap_){
        if(fDeltaR < fDeltaRVetoEndcapCharged_)
         return -999.;
      }
      if(bRectangleVetoEndcap_){
        if(std::abs(fDeltaEta) < fRectangleDeltaEtaVetoEndcapCharged_ && std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoEndcapCharged_){
         return -999.;
        }
      }
  }
                
  


  return fDeltaR;
}



//--------------------------------------------------------------------------------------------------
reco::VertexRef PFPhotonIsolationCalculator::chargedHadronVertex( edm::Handle< reco::VertexCollection > verticesColl, const reco::PFCandidate& pfcand ){

  //code copied from Florian's PFNoPU class (corrected removing the double loop....)
    
  auto const & track = pfcand.trackRef();

  size_t iVertex = 0;
  unsigned int index=0;
  unsigned int nFoundVertex = 0;

  float bestweight=0;
  
  const reco::VertexCollection& vertices = *(verticesColl.product());

  for( auto const & vtx :  vertices) {
    float w = vtx.trackWeight(track); // 0 if does not belong here
    //select the vertex for which the track has the highest weight
    if (w > bestweight){ // should we break here?
        bestweight=w;
    iVertex=index;
    nFoundVertex++;
    }
    ++index; 
  }
 
 
  
  if (nFoundVertex>0){
    if (nFoundVertex!=1)
      edm::LogWarning("TrackOnTwoVertex")<<"a track is shared by at least two verteces. Used to be an assert";
    return reco::VertexRef( verticesColl, iVertex);
  }
  // no vertex found with this track.

  // optional: as a secondary solution, associate the closest vertex in z
  if (  bCheckClosestZVertex_ ) {

    double dzmin = 10000.;
    double ztrack = pfcand.vertex().z();
    bool foundVertex = false;
    index = 0;
    for( reco::VertexCollection::const_iterator iv=vertices.begin(); iv!=vertices.end(); ++iv, ++index) {

      double dz = std::fabs(ztrack - iv->z());
      if(dz<dzmin) {
        dzmin = dz;
        iVertex = index;
        foundVertex = true;
      }
    }

    if( foundVertex )
      return reco::VertexRef( verticesColl, iVertex);
  
  }
   
  return reco::VertexRef( );
}



int PFPhotonIsolationCalculator::matchPFObject(const reco::Photon* photon, const reco::PFCandidateCollection* Candidates ){
  
  Int_t iMatch = -1;

  int i=0;
  for(reco::PFCandidateCollection::const_iterator iPF=Candidates->begin();iPF !=Candidates->end();iPF++){
    const reco::PFCandidate& pfParticle = (*iPF);
    if((((pfParticle.pdgId()==22 ) || std::abs(pfParticle.pdgId())==11) )){
     
      if(pfParticle.superClusterRef()==photon->superCluster())
        iMatch= i;
     
    }
    
    i++;
  }
  
/*
if(iMatch == -1){
i=0;
float fPt = -1;
for(reco::PFCandidateCollection::const_iterator iPF=Candidates->begin();iPF !=Candidates->end();iPF++){
const reco::PFCandidate& pfParticle = (*iPF);
if((((pfParticle.pdgId()==22 ) || TMath::Abs(pfParticle.pdgId())==11) )){
        if(pfParticle.pt()>fPt){
         fDeltaR = deltaR(pfParticle.eta(),pfParticle.phi(),photon->eta(),photon->phi());
         if(fDeltaR<0.1){
         iMatch = i;
         fPt = pfParticle.pt();
         }
        }
}
i++;
}
}
*/

  return iMatch;

}





int PFPhotonIsolationCalculator::matchPFObject(const reco::GsfElectron* electron, const reco::PFCandidateCollection* Candidates ){
  
  Int_t iMatch = -1;

  int i=0;
  for(reco::PFCandidateCollection::const_iterator iPF=Candidates->begin();iPF !=Candidates->end();iPF++){
    const reco::PFCandidate& pfParticle = (*iPF);
    if((((pfParticle.pdgId()==22 ) || std::abs(pfParticle.pdgId())==11) )){
     
      if(pfParticle.superClusterRef()==electron->superCluster())
        iMatch= i;
     
    }
    
    i++;
  }
  
  if(iMatch == -1){
    i=0;
    float fPt = -1;
    for(reco::PFCandidateCollection::const_iterator iPF=Candidates->begin();iPF !=Candidates->end();iPF++){
      const reco::PFCandidate& pfParticle = (*iPF);
      if((((pfParticle.pdgId()==22 ) || std::abs(pfParticle.pdgId())==11) )){
        if(pfParticle.pt()>fPt){
         float fDeltaR = deltaR(pfParticle.eta(),pfParticle.phi(),electron->eta(),electron->phi());
         if(fDeltaR<0.1){
         iMatch = i;
         fPt = pfParticle.pt();
         }
        }
      }
      i++;
    }
  }
  
  return iMatch;

}