HGCalIsoCalculator.cc

Go to the documentation of this file.

/*
 * HGCalIsoCalculator.cc
 *
 *  Created on: 13 Oct 2017
 *      Author: jkiesele, ncsmith
 */

#include "DataFormats/Math/interface/deltaR.h"
#include "RecoEgamma/EgammaTools//interface/HGCalIsoCalculator.h"

HGCalIsoCalculator::HGCalIsoCalculator():dr2_(0.15*0.15),mindr2_(0),rechittools_(nullptr),debug_(false),nlayers_(30){
    setNRings(5);
}

HGCalIsoCalculator::~HGCalIsoCalculator(){
}

void HGCalIsoCalculator::setRecHits(
        edm::Handle<HGCRecHitCollection> hitsEE,
        edm::Handle<HGCRecHitCollection> hitsFH,
        edm::Handle<HGCRecHitCollection> hitsBH
     ){
    recHitsEE_=hitsEE;
    recHitsFH_=hitsFH;
    recHitsBH_=hitsBH;

    if(!rechittools_)
        throw cms::Exception("HGCalIsoCalculator::produceHGCalIso: rechittools not set");
    
    hitEtaPhiCache_.clear();
    hitEtaPhiCache_.reserve(recHitsEE_->size()+recHitsFH_->size()+recHitsBH_->size());

    // Since HGCal is not projective and the rechits don't cache any
    // eta,phi, we make our own here
    auto makeEtaPhiPair = [this](const auto& hit) {
        const GlobalPoint position = rechittools_->getPosition(hit.id());
        float eta=rechittools_->getEta(position, 0); //assume vertex at z=0
        float phi=rechittools_->getPhi(position);
        return std::make_pair(eta, phi);
    };

    for(const auto& hit : *recHitsEE_) hitEtaPhiCache_.push_back(makeEtaPhiPair(hit));
    for(const auto& hit : *recHitsFH_) hitEtaPhiCache_.push_back(makeEtaPhiPair(hit));
    for(const auto& hit : *recHitsBH_) hitEtaPhiCache_.push_back(makeEtaPhiPair(hit));
}

void HGCalIsoCalculator::produceHGCalIso(const reco::CaloClusterPtr & seed){

    if(!rechittools_)
        throw cms::Exception("HGCalIsoCalculator::produceHGCalIso: rechittools not set");

    for(auto& r:isoringdeposits_)
        r=0;

    // make local temporaries to pass to the lambda
    // avoids recomputing every iteration
    const float seedEta=seed->eta();
    const float seedPhi=seed->phi();
    const std::vector<std::pair<DetId, float>> & seedHitsAndFractions = seed->hitsAndFractions();

    auto checkAndFill = [this, &seedEta, &seedPhi, &seedHitsAndFractions](const HGCRecHit& hit, std::pair<float,float> etaphiVal) {
        float deltar2=reco::deltaR2(etaphiVal.first,etaphiVal.second,seedEta,seedPhi);
        if(deltar2>dr2_ || deltar2<mindr2_) return;

        unsigned int layer=rechittools_->getLayerWithOffset(hit.id());
        if(layer>=nlayers_) return;

        //do not consider hits associated to the photon cluster
        if( std::none_of(seedHitsAndFractions.begin(), seedHitsAndFractions.end(), [&hit](const auto& seedhit){return hit.id()==seedhit.first;}) ){
            const unsigned int ring=ringasso_.at(layer);
            isoringdeposits_.at(ring)+=hit.energy();
        }
    };

    // The cache order is EE,FH,BH, so we should loop over them the same here
    auto itEtaPhiCache = hitEtaPhiCache_.cbegin();
    for(const auto& hit : *recHitsEE_){
      checkAndFill(hit, *itEtaPhiCache);
      itEtaPhiCache++;
    }
    for(const auto& hit : *recHitsFH_){
      checkAndFill(hit, *itEtaPhiCache);
      itEtaPhiCache++;
    }
    for(const auto& hit : *recHitsBH_){
      checkAndFill(hit, *itEtaPhiCache);
      itEtaPhiCache++;
    }
}

void HGCalIsoCalculator::setNRings(const size_t nrings){
    if(nrings>nlayers_)
        throw std::logic_error("PhotonHGCalIsoCalculator::setNRings: max number of rings reached");

    ringasso_.clear();
    isoringdeposits_.clear();
    unsigned int separator=nlayers_/nrings;
    size_t counter=0;
    for(size_t i=0;i<nlayers_+1;i++){
        ringasso_.push_back(counter);
        //the last ring might be larger.
        if(i && !(i%separator) && (int)counter<(int)nrings-1){
            counter++;
        }
    }
    isoringdeposits_.resize(nrings,0);
}

const float HGCalIsoCalculator::getIso(const unsigned int ring)const{
    if(ring>=isoringdeposits_.size())
        throw cms::Exception("HGCalIsoCalculator::getIso: ring index out of range");
    return isoringdeposits_[ring];
}