Skip to content

Commit

Permalink
Merge pull request #44513 from felicepantaleo/fasterCLUE2D
Browse files Browse the repository at this point in the history 
ticlv5: Optimize CLUE2D
  • Loading branch information
@cmsbuild
cmsbuild committed Apr 1, 2024
2 parents 1b3f760 + 14c46b8 commit a471799
Show file tree
Hide file tree
Showing 3 changed files with 69 additions and 9 deletions.
54 changes: 52 additions & 2 deletions RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,6 +61,9 @@ void HGCalCLUEAlgoT<T, STRATEGY>::populate(const HGCRecHitCollection& hits) {
const GlobalPoint position(rhtools_.getPosition(detid));
int offset = ((rhtools_.zside(detid) + 1) >> 1) * maxlayer_;
int layer = layerOnSide + offset;
// setting the layer position only once per layer
if (cells_[layer].layerDim3 == std::numeric_limits<float>::infinity())
cells_[layer].layerDim3 = position.z();

cells_[layer].detid.emplace_back(detid);
if constexpr (std::is_same_v<STRATEGY, HGCalScintillatorStrategy>) {
Expand Down Expand Up @@ -158,17 +161,64 @@ std::vector<reco::BasicCluster> HGCalCLUEAlgoT<T, STRATEGY>::getClusters(bool) {
std::vector<std::pair<DetId, float>> thisCluster;

for (auto& cl : cellsIdInCluster) {
math::XYZPoint position = math::XYZPoint(0.f, 0.f, 0.f);
float maxEnergyValue = std::numeric_limits<float>::min();
int maxEnergyCellIndex = -1;
DetId maxEnergyDetId;
float energy = 0.f;
int seedDetId = -1;

float x = 0.f;
float y = 0.f;
float z = cellsOnLayer.layerDim3;
// TODO Felice: maybe use the seed for the position calculation
for (auto cellIdx : cl) {
energy += cellsOnLayer.weight[cellIdx];
if (cellsOnLayer.weight[cellIdx] > maxEnergyValue) {
maxEnergyValue = cellsOnLayer.weight[cellIdx];
maxEnergyCellIndex = cellIdx;
maxEnergyDetId = cellsOnLayer.detid[cellIdx];
}
thisCluster.emplace_back(cellsOnLayer.detid[cellIdx], 1.f);
if (cellsOnLayer.isSeed[cellIdx]) {
seedDetId = cellsOnLayer.detid[cellIdx];
}
}

float total_weight_log = 0.f;
float total_weight = energy;

if constexpr (std::is_same_v<STRATEGY, HGCalSiliconStrategy>) {
auto thick = rhtools_.getSiThickIndex(maxEnergyDetId);
for (auto cellIdx : cl) {
const float d1 = cellsOnLayer.dim1[cellIdx] - cellsOnLayer.dim1[maxEnergyCellIndex];
const float d2 = cellsOnLayer.dim2[cellIdx] - cellsOnLayer.dim2[maxEnergyCellIndex];
if ((d1 * d1 + d2 * d2) < positionDeltaRho2_) {
float Wi = std::max(thresholdW0_[thick] + std::log(cellsOnLayer.weight[cellIdx] / energy), 0.);
x += cellsOnLayer.dim1[cellIdx] * Wi;
y += cellsOnLayer.dim2[cellIdx] * Wi;
total_weight_log += Wi;
}
}
} else {
for (auto cellIdx : cl) {
auto position = rhtools_.getPosition(cellsOnLayer.detid[cellIdx]);
x += position.x() * cellsOnLayer.weight[cellIdx];
y += position.y() * cellsOnLayer.weight[cellIdx];
}
}

if constexpr (std::is_same_v<STRATEGY, HGCalSiliconStrategy>) {
total_weight = total_weight_log;
}

if (total_weight != 0.) {
float inv_tot_weight = 1.f / total_weight;
x *= inv_tot_weight;
y *= inv_tot_weight;
} else {
x = y = 0.f;
}
math::XYZPoint position = math::XYZPoint(x, y, z);

auto globalClusterIndex = cellsOnLayer.clusterIndex[cl[0]] + firstClusterIdx;

clusters_v_[globalClusterIndex] =
Expand Down
7 changes: 6 additions & 1 deletion RecoLocalCalo/HGCalRecProducers/plugins/HGCalCLUEAlgo.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,9 @@
#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"

// C/C++ headers
#include <set>
#include <string>
#include <vector>
#include <limits>

template <typename TILE, typename STRATEGY>
class HGCalCLUEAlgoT : public HGCalClusteringAlgoBase {
Expand All @@ -45,6 +45,8 @@ class HGCalCLUEAlgoT : public HGCalClusteringAlgoBase {
fcPerEle_(ps.getParameter<double>("fcPerEle")),
nonAgedNoises_(ps.getParameter<std::vector<double>>("noises")),
noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("noise_MIP")),
thresholdW0_(ps.getParameter<std::vector<double>>("thresholdW0")),
positionDeltaRho2_(ps.getParameter<double>("positionDeltaRho2")),
use2x2_(ps.getParameter<bool>("use2x2")),
initialized_(false) {}

Expand Down Expand Up @@ -137,6 +139,8 @@ class HGCalCLUEAlgoT : public HGCalClusteringAlgoBase {
double noiseMip_;
std::vector<std::vector<double>> thresholds_;
std::vector<std::vector<double>> v_sigmaNoise_;
std::vector<double> thresholdW0_;
double positionDeltaRho2_;

bool use2x2_;

Expand All @@ -159,6 +163,7 @@ class HGCalCLUEAlgoT : public HGCalClusteringAlgoBase {
std::vector<float> sigmaNoise;
std::vector<std::vector<int>> followers;
std::vector<bool> isSeed;
float layerDim3 = std::numeric_limits<float>::infinity();

void clear() {
detid.clear();
Expand Down
17 changes: 11 additions & 6 deletions RecoLocalCalo/HGCalRecProducers/plugins/HGCalLayerClusterProducer.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,7 @@ class HGCalLayerClusterProducer : public edm::stream::EDProducer<> {
double positionDeltaRho2_;
hgcal::RecHitTools rhtools_;
edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
const bool calculatePositionInAlgo_;

/**
* @brief Sets algoId accordingly to the detector type
Expand Down Expand Up @@ -107,7 +108,8 @@ HGCalLayerClusterProducer::HGCalLayerClusterProducer(const edm::ParameterSet& ps
detector_(ps.getParameter<std::string>("detector")), // one of EE, FH, BH, HFNose
timeClname_(ps.getParameter<std::string>("timeClname")),
hitsTime_(ps.getParameter<unsigned int>("nHitsTime")),
caloGeomToken_(consumesCollector().esConsumes<CaloGeometry, CaloGeometryRecord>()) {
caloGeomToken_(consumesCollector().esConsumes<CaloGeometry, CaloGeometryRecord>()),
calculatePositionInAlgo_(ps.getParameter<bool>("calculatePositionInAlgo")) {
setAlgoId(); //sets algo id according to detector type
hits_token_ = consumes<HGCRecHitCollection>(ps.getParameter<edm::InputTag>("recHits"));

Expand Down Expand Up @@ -139,6 +141,7 @@ void HGCalLayerClusterProducer::fillDescriptions(edm::ConfigurationDescriptions&
desc.add<edm::InputTag>("recHits", edm::InputTag("HGCalRecHit", "HGCEERecHits"));
desc.add<std::string>("timeClname", "timeLayerCluster");
desc.add<unsigned int>("nHitsTime", 3);
desc.add<bool>("calculatePositionInAlgo", true);
descriptions.add("hgcalLayerClusters", desc);
}

Expand Down Expand Up @@ -212,10 +215,10 @@ std::pair<float, float> HGCalLayerClusterProducer::calculateTime(

float rhTimeE = rechit->timeError();
//check on timeError to exclude scintillator
if (rhTimeE < 0.)
if (rhTimeE < 0.f)
continue;
timeClhits.push_back(rechit->time());
timeErrorClhits.push_back(1. / (rhTimeE * rhTimeE));
timeErrorClhits.push_back(1.f / (rhTimeE * rhTimeE));
}
hgcalsimclustertime::ComputeClusterTime timeEstimator;
timeCl = timeEstimator.fixSizeHighestDensity(timeClhits, timeErrorClhits, hitsTime_);
Expand Down Expand Up @@ -249,12 +252,14 @@ void HGCalLayerClusterProducer::produce(edm::Event& evt, const edm::EventSetup&
times.reserve(clusters->size());

for (unsigned i = 0; i < clusters->size(); ++i) {
const reco::CaloCluster& sCl = (*clusters)[i];
(*clusters)[i].setPosition(calculatePosition(hitmap, sCl.hitsAndFractions()));
reco::CaloCluster& sCl = (*clusters)[i];
if (!calculatePositionInAlgo_) {
sCl.setPosition(calculatePosition(hitmap, sCl.hitsAndFractions()));
}
if (detector_ != "BH") {
times.push_back(calculateTime(hitmap, sCl.hitsAndFractions(), sCl.size()));
} else {
times.push_back(std::pair<float, float>(-99., -1.));
times.push_back(std::pair<float, float>(-99.f, -1.f));
}
}

Expand Down

0 comments on commit a471799

Please sign in to comment.