GsfEleCalPFClusterIsoCut.cc

Go to the documentation of this file.

#include "PhysicsTools/SelectorUtils/interface/CutApplicatorWithEventContentBase.h"
#include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
#include "RecoEgamma/EgammaTools/interface/EffectiveAreas.h"


class GsfEleCalPFClusterIsoCut : public CutApplicatorWithEventContentBase {
public:

  enum IsoType {UNDEF=-1, ISO_ECAL=0, ISO_HCAL=1};

  GsfEleCalPFClusterIsoCut(const edm::ParameterSet& c);
  
  result_type operator()(const reco::GsfElectronPtr&) const final;

  void setConsumes(edm::ConsumesCollector&) final;
  void getEventContent(const edm::EventBase&) final;

  double value(const reco::CandidatePtr& cand) const final;

  CandidateType candidateType() const final { 
    return ELECTRON; 
  }

private:
  // Cut values
  const float isoCutEBLowPt_,isoCutEBHighPt_,isoCutEELowPt_,isoCutEEHighPt_;
  // Configuration
  const int     isoType_;
  const float   ptCutOff_;
  const float   barrelCutOff_;
  bool          isRelativeIso_;
  // Effective area constants
  EffectiveAreas effectiveAreas_;
  // The rho
  edm::Handle< double > rhoHandle_;

  constexpr static char rhoString_     [] = "rho";
};

constexpr char GsfEleCalPFClusterIsoCut::rhoString_[];

DEFINE_EDM_PLUGIN(CutApplicatorFactory,
          GsfEleCalPFClusterIsoCut,
          "GsfEleCalPFClusterIsoCut");

GsfEleCalPFClusterIsoCut::GsfEleCalPFClusterIsoCut(const edm::ParameterSet& c) :
  CutApplicatorWithEventContentBase(c),
  isoCutEBLowPt_(c.getParameter<double>("isoCutEBLowPt")),
  isoCutEBHighPt_(c.getParameter<double>("isoCutEBHighPt")),
  isoCutEELowPt_(c.getParameter<double>("isoCutEELowPt")),
  isoCutEEHighPt_(c.getParameter<double>("isoCutEEHighPt")),
  isoType_(c.getParameter<int>("isoType")),
  ptCutOff_(c.getParameter<double>("ptCutOff")),
   barrelCutOff_(c.getParameter<double>("barrelCutOff")),
  isRelativeIso_(c.getParameter<bool>("isRelativeIso")),
  effectiveAreas_( (c.getParameter<edm::FileInPath>("effAreasConfigFile")).fullPath())
{
  
  edm::InputTag rhoTag = c.getParameter<edm::InputTag>("rho");
  contentTags_.emplace(rhoString_,rhoTag);

}

void GsfEleCalPFClusterIsoCut::setConsumes(edm::ConsumesCollector& cc) {

  auto rho = cc.consumes<double>(contentTags_[rhoString_]);
  contentTokens_.emplace(rhoString_, rho);
}

void GsfEleCalPFClusterIsoCut::getEventContent(const edm::EventBase& ev) {  

  ev.getByLabel(contentTags_[rhoString_],rhoHandle_);
}

CutApplicatorBase::result_type 
GsfEleCalPFClusterIsoCut::
operator()(const reco::GsfElectronPtr& cand) const{  

  // Establish the cut value
  double absEta = std::abs(cand->superCluster()->eta());

  const pat::Electron* elPat = dynamic_cast<const pat::Electron*>(cand.get());
  if( !elPat ) {
    throw cms::Exception("ERROR: this VID selection is meant to be run on miniAOD/PAT only")
      << std::endl << "Change input format to PAT/miniAOD or contact Egamma experts" 
      << std::endl << std::endl;
  }

  const float isoCut =
    ( cand->pt() < ptCutOff_ ?
      ( absEta < barrelCutOff_ ? isoCutEBLowPt_ : isoCutEELowPt_ ) 
      :
      ( absEta < barrelCutOff_ ? isoCutEBHighPt_ : isoCutEEHighPt_ ) );

  const float  eA = effectiveAreas_.getEffectiveArea( absEta );
  const float rho = rhoHandle_.isValid() ? (float)(*rhoHandle_) : 0; // std::max likes float arguments

  float isoValue = -999;
  if( isoType_ == ISO_ECAL ){
    isoValue = elPat->ecalPFClusterIso();
  }else if( isoType_ == ISO_HCAL ){
    isoValue = elPat->hcalPFClusterIso();
  }else{
    throw cms::Exception("ERROR: unknown type requested for PF cluster isolation.")
      << std::endl << "Check VID configuration." << std::endl;
  }
  float isoValueCorr = std::max(0.0f, isoValue - rho*eA);
  
  // Apply the cut and return the result
  // Scale by pT if the relative isolation is requested but avoid division by 0
  return isoValueCorr < isoCut*(isRelativeIso_ ? cand->pt() : 1.);
}

double GsfEleCalPFClusterIsoCut::value(const reco::CandidatePtr& cand) const {

  reco::GsfElectronPtr ele(cand);
  // Establish the cut value
  double absEta = std::abs(ele->superCluster()->eta());

  const pat::Electron *elPat = dynamic_cast<const pat::Electron*>(ele.get());
  if( !elPat ){
    throw cms::Exception("ERROR: this VID selection is meant to be run on miniAOD/PAT only")
      << std::endl << "Change input format to PAT/miniAOD or contact Egamma experts" 
      << std::endl << std::endl;
  }

  const float  eA = effectiveAreas_.getEffectiveArea( absEta );
  const float rho = rhoHandle_.isValid() ? (float)(*rhoHandle_) : 0; // std::max likes float arguments

  float isoValue = -999;
  if( isoType_ == ISO_ECAL ){
    isoValue = elPat->ecalPFClusterIso();
  }else if( isoType_ == ISO_HCAL ){
    isoValue = elPat->hcalPFClusterIso();
  }else{
    throw cms::Exception("ERROR: unknown type requested for PF cluster isolation.")
      << std::endl << "Check VID configuration." << std::endl;
  } 
  float isoValueCorr = std::max(0.0f, isoValue - rho*eA);
  
  // Divide by pT if the relative isolation is requested
  if( isRelativeIso_ )
    isoValueCorr /= ele->pt();

  // Apply the cut and return the result
  return isoValueCorr;
}