Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding HBHE rechits to miniAOD #23540

Merged
merged 2 commits into from Jun 22, 2018
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
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
72 changes: 72 additions & 0 deletions RecoEgamma/EgammaIsolationAlgos/interface/EGHcalRecHitSelector.h
@@ -0,0 +1,72 @@
#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);
float getMinEnergyHCAL_(HcalDetId id)const;

int maxDIEta_;
int maxDIPhi_;
float minEnergyHB_;
float minEnergyHEDepth1_;
float minEnergyHEDefault_;

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";
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe we should better return here

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

return;
}
//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()>getMinEnergyHCAL_(recHit.id())) 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,23 @@
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),
minEnergyHB = cms.double(0.8),
minEnergyHEDepth1 = cms.double(0.1),
minEnergyHEDefault = cms.double(0.2),
)
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.

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

#include<limits>

EGHcalRecHitSelector::EGHcalRecHitSelector(const edm::ParameterSet& config):
maxDIEta_(config.getParameter<int>("maxDIEta")),
maxDIPhi_(config.getParameter<int>("maxDIPhi")),
minEnergyHB_(config.getParameter<double>("minEnergyHB")),
minEnergyHEDepth1_(config.getParameter<double>("minEnergyHEDepth1")),
minEnergyHEDefault_(config.getParameter<double>("minEnergyHEDefault"))
{

}

edm::ParameterSetDescription EGHcalRecHitSelector::makePSetDescription()
{
edm::ParameterSetDescription desc;
desc.add<int>("maxDIEta",5);
desc.add<int>("maxDIPhi",5);
desc.add<double>("minEnergyHB",0.8);
desc.add<double>("minEnergyHEDepth1",0.1);
desc.add<double>("minEnergyHEDefault",0.2);
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;
}

float EGHcalRecHitSelector::getMinEnergyHCAL_(HcalDetId id)const
{
if(id.subdetId()==HcalBarrel) return minEnergyHB_;
else if(id.subdetId()==HcalEndcap){
if(id.depth()==1) return minEnergyHEDepth1_;
else return minEnergyHEDefault_;
}else return std::numeric_limits<float>::max();

}