Skip to content

Commit

Permalink
Merge pull request #29049 from hatakeyamak/PFRecHitNavigatorForHcalWi…
Browse files Browse the repository at this point in the history
…thDenseId

PFRecHit navigator for HCAL with dense index
  • Loading branch information
cmsbuild committed Mar 9, 2020
2 parents 7f96e94 + 29da198 commit 761c1de
Show file tree
Hide file tree
Showing 13 changed files with 262 additions and 951 deletions.
31 changes: 31 additions & 0 deletions HLTrigger/Configuration/python/customizeHLTforCMSSW.py
Expand Up @@ -164,10 +164,41 @@ def customiseFor2017DtUnpacking(process):

return process

def customiseFor29049(process) :

listHltPFRecHitHBHE=['hltParticleFlowRecHitHBHE',
'hltParticleFlowRecHitHBHEForEgamma',
'hltParticleFlowRecHitHBHEForEgammaUnseeded',
'hltParticleFlowRecHitHBHEForMuons',
'hltParticleFlowRecHitHBHERegForMuons']
for att in listHltPFRecHitHBHE:
if hasattr(process,att):
prod = getattr(process, att)
pset_navi = prod.navigator
if hasattr(pset_navi, "sigmaCut"): delattr(pset_navi,'sigmaCut')
if hasattr(pset_navi, "timeResolutionCalc"): delattr(pset_navi,'timeResolutionCalc')
pset_navi.name = cms.string("PFRecHitHCALDenseIdNavigator")
pset_navi.hcalEnums = cms.vint32(1,2)

listHltPFRecHitHF=['hltParticleFlowRecHitHF',
'hltParticleFlowRecHitHFForEgammaUnseeded']
for att in listHltPFRecHitHF:
if hasattr(process,att):
prod = getattr(process, att)
pset_navi = prod.navigator
if hasattr(pset_navi, "barrel"): delattr(pset_navi,'barrel')
if hasattr(pset_navi, "endcap"): delattr(pset_navi,'endcap')
pset_navi.name = cms.string("PFRecHitHCALDenseIdNavigator")
pset_navi.hcalEnums = cms.vint32(4)

return process

# CMSSW version specific customizations
def customizeHLTforCMSSW(process, menuType="GRun"):

# add call to action function in proper order: newest last!
# process = customiseFor12718(process)

process = customiseFor29049(process)

return process
Expand Up @@ -59,13 +59,13 @@ class HGCRecHitNavigator : public PFRecHitNavigatorBase {
}
}

void beginEvent(const edm::EventSetup& iSetup) override {
void init(const edm::EventSetup& iSetup) override {
if (nullptr != eeNav_)
eeNav_->beginEvent(iSetup);
eeNav_->init(iSetup);
if (nullptr != hefNav_)
hefNav_->beginEvent(iSetup);
hefNav_->init(iSetup);
if (nullptr != hebNav_)
hebNav_->beginEvent(iSetup);
hebNav_->init(iSetup);
}

void associateNeighbours(reco::PFRecHit& hit,
Expand Down
Expand Up @@ -26,7 +26,7 @@ class PFECALHashNavigator : public PFRecHitNavigatorBase {
neighbourmapcalculated_ = false;
}

void beginEvent(const edm::EventSetup& iSetup) override {
void init(const edm::EventSetup& iSetup) override {
edm::ESHandle<CaloGeometry> geoHandle;
iSetup.get<CaloGeometryRecord>().get(geoHandle);

Expand Down
192 changes: 192 additions & 0 deletions RecoParticleFlow/PFClusterProducer/interface/PFHCALDenseIdNavigator.h
@@ -0,0 +1,192 @@
#ifndef RecoParticleFlow_PFClusterProducer_PFHCALDenseIdNavigator_h
#define RecoParticleFlow_PFClusterProducer_PFHCALDenseIdNavigator_h

#include "FWCore/Framework/interface/ESWatcher.h"

#include "RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"

#include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/EcalDetId/interface/ESDetId.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"

#include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
#include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
#include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
#include "Geometry/CaloTopology/interface/HcalTopology.h"

#include "Geometry/CaloTopology/interface/CaloTowerTopology.h"
#include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"

template <typename DET, typename TOPO, bool ownsTopo = true>
class PFHCALDenseIdNavigator : public PFRecHitNavigatorBase {
public:
~PFHCALDenseIdNavigator() override {
if (!ownsTopo) {
topology_.release();
}
}

PFHCALDenseIdNavigator(const edm::ParameterSet& iConfig) {
vhcalEnum_ = iConfig.getParameter<std::vector<int>>("hcalEnums");
}

void init(const edm::EventSetup& iSetup) override {
bool check = theRecNumberWatcher_.check(iSetup);
if (!check)
return;

edm::ESHandle<HcalTopology> hcalTopology;
iSetup.get<HcalRecNumberingRecord>().get(hcalTopology);
topology_.release();
topology_.reset(hcalTopology.product());

// Fill a vector of valid denseid's
edm::ESHandle<CaloGeometry> hGeom;
iSetup.get<CaloGeometryRecord>().get(hGeom);
const CaloGeometry& caloGeom = *hGeom;

std::vector<DetId> vecHcal;
std::vector<unsigned int> vDenseIdHcal;
neighboursHcal_.clear();
for (auto hcalSubdet : vhcalEnum_) {
std::vector<DetId> vecDetIds(caloGeom.getValidDetIds(DetId::Hcal, hcalSubdet));
vecHcal.insert(vecHcal.end(), vecDetIds.begin(), vecDetIds.end());
}
for (auto hDetId : vecHcal) {
vDenseIdHcal.push_back(topology_.get()->detId2denseId(hDetId));
}
std::sort(vDenseIdHcal.begin(), vDenseIdHcal.end());

// Fill a vector of cell neighbours
denseIdHcalMax_ = *max_element(vDenseIdHcal.begin(), vDenseIdHcal.end());
denseIdHcalMin_ = *min_element(vDenseIdHcal.begin(), vDenseIdHcal.end());
neighboursHcal_.resize(denseIdHcalMax_ - denseIdHcalMin_ + 1);

for (auto denseid : vDenseIdHcal) {
DetId N(0);
DetId E(0);
DetId S(0);
DetId W(0);
DetId NW(0);
DetId NE(0);
DetId SW(0);
DetId SE(0);
std::vector<DetId> neighbours(9, DetId(0));

// the centre
unsigned denseid_c = denseid;
DetId detid_c = topology_.get()->denseId2detId(denseid_c);
CaloNavigator<DET> navigator(detid_c, topology_.get());

// Using enum in Geometry/CaloTopology/interface/CaloDirection.h
// Order: CENTER(NONE),SOUTH,SOUTHEAST,SOUTHWEST,EAST,WEST,NORTHEAST,NORTHWEST,NORTH
neighbours.at(NONE) = detid_c;

navigator.home();
N = navigator.north();
neighbours.at(NORTH) = N;
if (N != DetId(0)) {
NE = navigator.east();
} else {
navigator.home();
E = navigator.east();
NE = navigator.north();
}
neighbours.at(NORTHEAST) = NE;

navigator.home();
S = navigator.south();
neighbours.at(SOUTH) = S;
if (S != DetId(0)) {
SW = navigator.west();
} else {
navigator.home();
W = navigator.west();
SW = navigator.south();
}
neighbours.at(SOUTHWEST) = SW;

navigator.home();
E = navigator.east();
neighbours.at(EAST) = E;
if (E != DetId(0)) {
SE = navigator.south();
} else {
navigator.home();
S = navigator.south();
SE = navigator.east();
}
neighbours.at(SOUTHEAST) = SE;

navigator.home();
W = navigator.west();
neighbours.at(WEST) = W;
if (W != DetId(0)) {
NW = navigator.north();
} else {
navigator.home();
N = navigator.north();
NW = navigator.west();
}
neighbours.at(NORTHWEST) = NW;

unsigned index = getIdx(denseid_c);
neighboursHcal_[index] = neighbours;
}
}

void associateNeighbours(reco::PFRecHit& hit,
std::unique_ptr<reco::PFRecHitCollection>& hits,
edm::RefProd<reco::PFRecHitCollection>& refProd) override {
DetId detid(hit.detId());
unsigned denseid = topology_.get()->detId2denseId(detid);

std::vector<DetId> neighbours(9, DetId(0));

if (denseid < denseIdHcalMin_ || denseid > denseIdHcalMax_) {
edm::LogWarning("PFRecHitHCALCachedNavigator") << " DenseId for this cell is out of the range." << std::endl;
} else if (!validNeighbours(denseid)) {
edm::LogWarning("PFRecHitHCALCachedNavigator")
<< " DenseId for this cell does not have the neighbour information." << std::endl;
} else {
unsigned index = getIdx(denseid);
neighbours = neighboursHcal_.at(index);
}

associateNeighbour(neighbours.at(NORTH), hit, hits, refProd, 0, 1, 0); // N
associateNeighbour(neighbours.at(NORTHEAST), hit, hits, refProd, 1, 1, 0); // NE
associateNeighbour(neighbours.at(SOUTH), hit, hits, refProd, 0, -1, 0); // S
associateNeighbour(neighbours.at(SOUTHWEST), hit, hits, refProd, -1, -1, 0); // SW
associateNeighbour(neighbours.at(EAST), hit, hits, refProd, 1, 0, 0); // E
associateNeighbour(neighbours.at(SOUTHEAST), hit, hits, refProd, 1, -1, 0); // SE
associateNeighbour(neighbours.at(WEST), hit, hits, refProd, -1, 0, 0); // W
associateNeighbour(neighbours.at(NORTHWEST), hit, hits, refProd, -1, 1, 0); // NW
}

bool validNeighbours(const unsigned int denseid) const {
bool ok = true;
unsigned index = getIdx(denseid);
if (neighboursHcal_.at(index).size() != 9)
ok = false; // the neighbour vector size should be 3x3
return ok;
}

unsigned int getIdx(const unsigned int denseid) const {
unsigned index = denseid - denseIdHcalMin_;
return index;
}

protected:
edm::ESWatcher<HcalRecNumberingRecord> theRecNumberWatcher_;
std::unique_ptr<const TOPO> topology_;
std::vector<int> vhcalEnum_;
std::vector<std::vector<DetId>> neighboursHcal_;
unsigned int denseIdHcalMax_;
unsigned int denseIdHcalMin_;
};

#endif
Expand Up @@ -14,9 +14,9 @@ class PFRecHitDualNavigator : public PFRecHitNavigatorBase {
endcapNav_ = new endcap(iConfig.getParameter<edm::ParameterSet>("endcap"));
}

void beginEvent(const edm::EventSetup& iSetup) override {
barrelNav_->beginEvent(iSetup);
endcapNav_->beginEvent(iSetup);
void init(const edm::EventSetup& iSetup) override {
barrelNav_->init(iSetup);
endcapNav_->init(iSetup);
}

void associateNeighbours(reco::PFRecHit& hit,
Expand Down
Expand Up @@ -32,7 +32,7 @@ class PFRecHitNavigatorBase {

virtual ~PFRecHitNavigatorBase() = default;

virtual void beginEvent(const edm::EventSetup&) = 0;
virtual void init(const edm::EventSetup&) = 0;
virtual void associateNeighbours(reco::PFRecHit&,
std::unique_ptr<reco::PFRecHitCollection>&,
edm::RefProd<reco::PFRecHitCollection>&) = 0;
Expand Down

0 comments on commit 761c1de

Please sign in to comment.