Skip to content

Commit

Permalink
adding HCAL rechits to miniAOD
Browse files Browse the repository at this point in the history
  • Loading branch information
Sam-Harper committed Jun 8, 2018
1 parent 30c7b03 commit 2e888d5
Show file tree
Hide file tree
Showing 10 changed files with 203 additions and 104 deletions.
Expand Up @@ -32,7 +32,7 @@
'keep recoCaloClusters_reducedEgamma_*_*',
'keep EcalRecHitsSorted_reducedEgamma_*_*',
'keep recoGsfTracks_reducedEgamma_*_*',

'keep HBHERecHitsSorted_reducedEgamma_*_*',
'drop *_*_caloTowers_*',
'drop *_*_pfCandidates_*',
'drop *_*_genJets_*',
Expand Down
67 changes: 67 additions & 0 deletions RecoEgamma/EgammaIsolationAlgos/interface/EGHcalRecHitSelector.h
@@ -0,0 +1,67 @@
#ifndef RecoEgamam_EgammaIsolationAlgos_EGHcalRecHitSelector_h
#define RecoEgamam_EgammaIsolationAlgos_EGHcalRecHitSelector_h

#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"

#include "Geometry/CaloTopology/interface/CaloTowerConstituentsMap.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"

#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"

class EGHcalRecHitSelector {
public:
explicit EGHcalRecHitSelector(const edm::ParameterSet& config);
~EGHcalRecHitSelector(){}

template<typename CollType>
void addDetIds(const reco::SuperCluster& superClus,const HBHERecHitCollection& recHits,
CollType& detIdsToStore)const;

void setup(const edm::EventSetup& iSetup){iSetup.get<CaloGeometryRecord>().get(towerMap_);}

static edm::ParameterSetDescription makePSetDescription();

private:
static int calDIPhi(int iPhi1,int iPhi2);
static int calDIEta(int iEta1,int iEta2);

int maxDIEta_;
int maxDIPhi_;
float minEnergyHCAL_;
edm::ESHandle<CaloTowerConstituentsMap> towerMap_;

};

template <typename CollType>
void EGHcalRecHitSelector::
addDetIds(const reco::SuperCluster& superClus,const HBHERecHitCollection& recHits,
CollType& detIdsToStore)const
{
DetId seedId = superClus.seed()->seed();
if(seedId.det() != DetId::Ecal && seedId.det() != DetId::Forward) {
edm::LogError("EgammaIsoHcalDetIdCollectionProducerError") << "Somehow the supercluster has a seed which is not ECAL, something is badly wrong";
}
//so we are using CaloTowers to get the iEta/iPhi of the HCAL rec hit behind the seed cluster, this might go funny on tower 28 but shouldnt matter there

if( seedId.det() == DetId::Forward ) return;

CaloTowerDetId towerId(towerMap_->towerOf(seedId));
int seedHcalIEta = towerId.ieta();
int seedHcalIPhi = towerId.iphi();

for(auto& recHit : recHits){
int dIEtaAbs = std::abs(calDIEta(seedHcalIEta,recHit.id().ieta()));
int dIPhiAbs = std::abs(calDIPhi(seedHcalIPhi,recHit.id().iphi()));

if(dIEtaAbs<=maxDIEta_ && dIPhiAbs<=maxDIPhi_ && recHit.energy()>minEnergyHCAL_) detIdsToStore.insert(detIdsToStore.end(),recHit.id().rawId());
}
}


#endif
Expand Up @@ -6,7 +6,8 @@
#include "DataFormats/DetId/interface/DetIdCollection.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"

EgammaIsoHcalDetIdCollectionProducer::EgammaIsoHcalDetIdCollectionProducer(const edm::ParameterSet& iConfig)
EgammaIsoHcalDetIdCollectionProducer::EgammaIsoHcalDetIdCollectionProducer(const edm::ParameterSet& iConfig):
hcalHitSelector_(iConfig.getParameter<edm::ParameterSet>("hitSelection"))
{

recHitsToken_ =
Expand All @@ -27,31 +28,36 @@ EgammaIsoHcalDetIdCollectionProducer::EgammaIsoHcalDetIdCollectionProducer(const

interestingDetIdCollection_ = iConfig.getParameter<std::string>("interestingDetIdCollection");

maxDIEta_ = iConfig.getParameter<int>("maxDIEta");
maxDIPhi_ = iConfig.getParameter<int>("maxDIPhi");

minEnergyHCAL_= iConfig.getParameter<double>("minEnergyHCAL");

//register your products
produces< DetIdCollection > (interestingDetIdCollection_) ;
}


void EgammaIsoHcalDetIdCollectionProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions)
{
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("recHitsLabel",edm::InputTag("hbhereco"));
desc.add<edm::InputTag>("elesLabel",edm::InputTag("gedGsfElectrons"));
desc.add<edm::InputTag>("phosLabel",edm::InputTag("gedPhotons"));
desc.add<edm::InputTag>("superClustersLabel",edm::InputTag("particleFlowEGamma"));
desc.add<double>("minSCEt",20);
desc.add<double>("minEleEt",20);
desc.add<double>("minPhoEt",20);
desc.add<std::string>("interestingDetIdCollection","");
desc.add<edm::ParameterSetDescription>("hitSelection",EGHcalRecHitSelector::makePSetDescription());
descriptions.add(("interestingGedEgammaIsoHCALDetId"), desc);
}


void EgammaIsoHcalDetIdCollectionProducer::beginRun (edm::Run const& run, const edm::EventSetup & iSetup)
{
iSetup.get<CaloGeometryRecord>().get(towerMap_);
// std::cout <<" got geom "<<towerMap_.isValid()<<" "<<&(*towerMap_)<<std::endl;
hcalHitSelector_.setup(iSetup);
}

// ------------ method called to produce the data ------------
void
EgammaIsoHcalDetIdCollectionProducer::produce (edm::Event& iEvent,
const edm::EventSetup& iSetup)
{

// take BasicClusters
edm::Handle<reco::SuperClusterCollection> superClusters;
iEvent.getByToken(superClustersToken_, superClusters);

Expand All @@ -64,27 +70,31 @@ EgammaIsoHcalDetIdCollectionProducer::produce (edm::Event& iEvent,
edm::Handle<HBHERecHitCollection> recHits;
iEvent.getByToken(recHitsToken_,recHits);

//Create empty output collections
std::vector<DetId> indexToStore;
indexToStore.reserve(100);

if(eles.isValid() && recHits.isValid()){
for(auto& ele : *eles){

float scEt = ele.superCluster()->energy()*std::sin(ele.superCluster()->position().theta());
if(scEt > minEleEt_ || ele.et()> minEleEt_) addDetIds(*ele.superCluster(),*recHits,indexToStore);
if(scEt > minEleEt_ || ele.et()> minEleEt_){
hcalHitSelector_.addDetIds(*ele.superCluster(),*recHits,indexToStore);
}
}
}
if(phos.isValid() && recHits.isValid()){
for(auto& pho : *phos){
float scEt = pho.superCluster()->energy()*std::sin(pho.superCluster()->position().theta());
if(scEt > minPhoEt_ || pho.et()> minPhoEt_) addDetIds(*pho.superCluster(),*recHits,indexToStore);
if(scEt > minPhoEt_ || pho.et()> minPhoEt_){
hcalHitSelector_.addDetIds(*pho.superCluster(),*recHits,indexToStore);
}
}
}
if(superClusters.isValid() && recHits.isValid()){
for(auto& sc : *superClusters){
float scEt = sc.energy()*std::sin(sc.position().theta());
if(scEt > minSCEt_) addDetIds(sc,*recHits,indexToStore);
if(scEt > minSCEt_){
hcalHitSelector_.addDetIds(sc,*recHits,indexToStore);
}
}
}

Expand All @@ -97,53 +107,3 @@ EgammaIsoHcalDetIdCollectionProducer::produce (edm::Event& iEvent,
iEvent.put(std::move(detIdCollection), interestingDetIdCollection_ );

}

//some nasty hardcoded badness
int calDIEta(int iEta1,int iEta2)
{

int dEta = iEta1-iEta2;
if(iEta1*iEta2<0) {//-ve to +ve transistion and no crystal at zero
if(dEta<0) dEta++;
else dEta--;
}
return dEta;
}

//some nasty hardcoded badness
int calDIPhi(int iPhi1,int iPhi2)
{

int dPhi = iPhi1-iPhi2;

if(dPhi>72/2) dPhi-=72;
else if(dPhi<-72/2) dPhi+=72;

return dPhi;

}


void
EgammaIsoHcalDetIdCollectionProducer::addDetIds(const reco::SuperCluster& superClus,const HBHERecHitCollection& recHits,std::vector<DetId>& detIdsToStore)
{
DetId seedId = superClus.seed()->seed();
if(seedId.det() != DetId::Ecal && seedId.det() != DetId::Forward) {
edm::LogError("EgammaIsoHcalDetIdCollectionProducerError") << "Somehow the supercluster has a seed which is not ECAL, something is badly wrong";
}
//so we are using CaloTowers to get the iEta/iPhi of the HCAL rec hit behind the seed cluster, this might go funny on tower 28 but shouldnt matter there

if( seedId.det() == DetId::Forward ) return;

CaloTowerDetId towerId(towerMap_->towerOf(seedId));
int seedHcalIEta = towerId.ieta();
int seedHcalIPhi = towerId.iphi();

for(auto& recHit : recHits){
int dIEtaAbs = std::abs(calDIEta(seedHcalIEta,recHit.id().ieta()));
int dIPhiAbs = std::abs(calDIPhi(seedHcalIPhi,recHit.id().iphi()));

if(dIEtaAbs<=maxDIEta_ && dIPhiAbs<=maxDIPhi_ && recHit.energy()>minEnergyHCAL_) detIdsToStore.push_back(recHit.id().rawId());
}

}
Expand Up @@ -37,21 +37,20 @@ Modified from the ECAL version "InterestingDetIdCollectionProducer" to be HCAL
#include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
#include "DataFormats/EgammaCandidates/interface/Photon.h"

#include "Geometry/CaloTopology/interface/CaloTowerConstituentsMap.h"

#include "RecoEgamma/EgammaIsolationAlgos/interface/EGHcalRecHitSelector.h"


class EgammaIsoHcalDetIdCollectionProducer : public edm::stream::EDProducer<> {
public:
//! ctor
explicit EgammaIsoHcalDetIdCollectionProducer(const edm::ParameterSet&);
static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

void beginRun (edm::Run const&, const edm::EventSetup&) final;
//! producer
void produce(edm::Event &, const edm::EventSetup&) override;

private:
void addDetIds(const reco::SuperCluster& superClus,const HBHERecHitCollection& recHits,std::vector<DetId>& detIdsToStore);


private:
// ----------member data ---------------------------
edm::EDGetTokenT<HBHERecHitCollection> recHitsToken_;
edm::EDGetTokenT<reco::SuperClusterCollection> superClustersToken_;
Expand All @@ -64,12 +63,7 @@ class EgammaIsoHcalDetIdCollectionProducer : public edm::stream::EDProducer<> {
float minEleEt_;
float minPhoEt_;

int maxDIEta_;
int maxDIPhi_;

float minEnergyHCAL_;

edm::ESHandle<CaloTowerConstituentsMap> towerMap_;
EGHcalRecHitSelector hcalHitSelector_;

};

Expand Down
Expand Up @@ -68,18 +68,21 @@
interestingGamIsoDetIdEE.outerRadius = 0.6
interestingGamIsoDetIdEE.innerRadius = 0.0

import RecoEgamma.EgammaIsolationAlgos.interestingEgammaIsoHCALDetIdModule_cff
interestingGedEgammaIsoHCALDetId = RecoEgamma.EgammaIsolationAlgos.interestingEgammaIsoHCALDetIdModule_cff.interestingEgammaIsoHCALDetId.clone()
import RecoEgamma.EgammaIsolationAlgos.interestingGedEgammaIsoHCALDetId_cfi
interestingGedEgammaIsoHCALDetId = RecoEgamma.EgammaIsolationAlgos.interestingGedEgammaIsoHCALDetId_cfi.interestingGedEgammaIsoHCALDetId.clone()
interestingEgammaIsoHCALSel = cms.PSet(
maxDIEta=cms.int32(5),
maxDIPhi=cms.int32(5),
minEnergyHCAL = cms.double(0.8),
)
interestingGedEgammaIsoHCALDetId.recHitsLabel=cms.InputTag("hbhereco")
interestingGedEgammaIsoHCALDetId.elesLabel=cms.InputTag("gedGsfElectrons")
interestingGedEgammaIsoHCALDetId.phosLabel=cms.InputTag("gedPhotons")
interestingGedEgammaIsoHCALDetId.superClustersLabel=cms.InputTag("particleFlowEGamma")
interestingGedEgammaIsoHCALDetId.minSCEt=cms.double(20)
interestingGedEgammaIsoHCALDetId.minEleEt=cms.double(20)
interestingGedEgammaIsoHCALDetId.minPhoEt=cms.double(20)
interestingGedEgammaIsoHCALDetId.maxDIEta=cms.int32(5)
interestingGedEgammaIsoHCALDetId.maxDIPhi=cms.int32(5)
interestingGedEgammaIsoHCALDetId.minEnergyHCAL = cms.double(0.8)
interestingGedEgammaIsoHCALDetId.hitSelection=interestingEgammaIsoHCALSel

## OOT Photons
interestingOotEgammaIsoHCALDetId = interestingGedEgammaIsoHCALDetId.clone()
Expand Down

This file was deleted.

37 changes: 37 additions & 0 deletions RecoEgamma/EgammaIsolationAlgos/src/EGHcalRecHitSelector.cc
@@ -0,0 +1,37 @@
#include "RecoEgamma/EgammaIsolationAlgos/interface/EGHcalRecHitSelector.h"

EGHcalRecHitSelector::EGHcalRecHitSelector(const edm::ParameterSet& config):
maxDIEta_(config.getParameter<int>("maxDIEta")),
maxDIPhi_(config.getParameter<int>("maxDIPhi")),
minEnergyHCAL_(config.getParameter<double>("minEnergyHCAL"))
{

}

edm::ParameterSetDescription EGHcalRecHitSelector::makePSetDescription()
{
edm::ParameterSetDescription desc;
desc.add<int>("maxDIEta",5);
desc.add<int>("maxDIPhi",5);
desc.add<double>("minEnergyHCAL",0.8);
return desc;
}

int EGHcalRecHitSelector::calDIEta(int iEta1,int iEta2)
{
int dEta = iEta1-iEta2;
if(iEta1*iEta2<0) {//-ve to +ve transistion and no crystal at zero
if(dEta<0) dEta++;
else dEta--;
}
return dEta;
}


int EGHcalRecHitSelector::calDIPhi(int iPhi1,int iPhi2)
{
int dPhi = iPhi1-iPhi2;
if(dPhi>72/2) dPhi-=72;
else if(dPhi<-72/2) dPhi+=72;
return dPhi;
}

0 comments on commit 2e888d5

Please sign in to comment.