EgammaHcalIsolation.cc

Go to the documentation of this file.

//*****************************************************************************
// File:      EgammaHcalIsolation.cc
// ----------------------------------------------------------------------------
// OrigAuth:  Matthias Mozer
// Institute: IIHE-VUB
//=============================================================================
//*****************************************************************************
//C++ includes
#include <vector>
#include <functional>

//ROOT includes
#include <Math/VectorUtil.h>

//CMSSW includes
#include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHcalIsolation.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "RecoCaloTools/Selectors/interface/CaloConeSelector.h"
#include "DataFormats/RecoCandidate/interface/RecoCandidate.h"

using namespace std;

double scaleToE(const double& eta) { return 1.0; }
double scaleToEt(const double& eta) { return sin(2*atan(exp(-eta))); }

EgammaHcalIsolation::EgammaHcalIsolation (
        double extRadius,
        double intRadius,
        double eLowB,
        double eLowE,
        double etLowB,
        double etLowE,
        edm::ESHandle<CaloGeometry> theCaloGeom ,
        HBHERecHitMetaCollection*  mhbhe
) :
  extRadius_(extRadius),
  intRadius_(intRadius),
  eLowB_(eLowB),
  eLowE_(eLowE),
  etLowB_(etLowB),
  etLowE_(etLowE),
  theCaloGeom_(theCaloGeom) ,  
  mhbhe_(mhbhe)
{
    //set up the geometry and selector
    const CaloGeometry* caloGeom = theCaloGeom_.product();
    doubleConeSel_ = new CaloDualConeSelector (intRadius_ ,extRadius_, caloGeom, DetId::Hcal);
}

EgammaHcalIsolation::~EgammaHcalIsolation ()
{
    delete doubleConeSel_;
}

double EgammaHcalIsolation::getHcalSum(const GlobalPoint &pclu, const HcalDepth &depth,
                                       double(*scale)(const double&) ) const
{
    double sum = 0.;
    if (mhbhe_) 
    {
        //Compute the HCAL energy behind ECAL
        double eta;
        std::auto_ptr<CaloRecHitMetaCollectionV> chosen = doubleConeSel_->select(pclu,*mhbhe_);
        CaloRecHitMetaCollectionV::const_iterator i;
        for (i = chosen->begin () ; i!= chosen->end () ; ++i) 
        {
            eta = theCaloGeom_.product()->getPosition(i->detid()).eta();
            HcalDetId hcalDetId(i->detid());
            if(hcalDetId.subdet() == HcalBarrel &&              //Is it in the barrel?
               i->energy() > eLowB_ &&                          //Does it pass the min energy?
               i->energy()*scaleToEt(eta) > etLowB_ &&          //Does it pass the min et?
               (depth == AllDepths || depth == Depth1)) {                    //Are we asking for the first depth?
                    sum += i->energy() * scale(eta);
            }
            if(hcalDetId.subdet() == HcalEndcap &&              //Is it in the endcap?
               i->energy() > eLowE_ &&                          //Does it pass the min energy?
               i->energy()*scaleToEt(eta) > etLowE_ ) {         //Does it pass the min et?
               switch(depth) {                                  //Which depth?
                    case AllDepths: sum += i->energy() * scale(eta); break;
                    case Depth1: sum += (isDepth2(i->detid())) ? 0 : i->energy() * scale(eta); break;
                    case Depth2: sum += (isDepth2(i->detid())) ? i->energy() * scale(eta) : 0; break;
               }
            }
        }
    } 
    return sum ;
}

bool EgammaHcalIsolation::isDepth2(const DetId& detId) const {
    
    if( (HcalDetId(detId).depth()==2 && HcalDetId(detId).ietaAbs()>=18 && HcalDetId(detId).ietaAbs()<27) ||
        (HcalDetId(detId).depth()==3 && HcalDetId(detId).ietaAbs()==27) ) {

        return true;

   } else {
        return false;
   }
}