CMS 3D CMS Logo

List of all members | Classes | Public Types | Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes
HGCalCLUEAlgoT< TILE, STRATEGY > Class Template Reference

#include <HGCalCLUEAlgo.h>

Inheritance diagram for HGCalCLUEAlgoT< TILE, STRATEGY >:
HGCalClusteringAlgoBase

Classes

struct  CellsOnLayer
 

Public Types

typedef math::XYZPoint Point
 point in the space More...
 
- Public Types inherited from HGCalClusteringAlgoBase
enum  VerbosityLevel { pDEBUG = 0, pWARNING = 1, pINFO = 2, pERROR = 3 }
 

Public Member Functions

void computeThreshold ()
 
std::vector< reco::BasicClustergetClusters (bool) override
 
void getEventSetupPerAlgorithm (const edm::EventSetup &es) override
 
 HGCalCLUEAlgoT (const edm::ParameterSet &ps)
 
void makeClusters () override
 
void populate (const HGCRecHitCollection &hits) override
 
void reset () override
 
 ~HGCalCLUEAlgoT () override
 
- Public Member Functions inherited from HGCalClusteringAlgoBase
virtual hgcal_clustering::Density getDensity ()
 
void getEventSetup (const edm::EventSetup &es, hgcal::RecHitTools rhtools)
 
 HGCalClusteringAlgoBase (VerbosityLevel v, reco::CaloCluster::AlgoId algo)
 
void setAlgoId (reco::CaloCluster::AlgoId algo, bool isNose=false)
 
void setVerbosity (VerbosityLevel the_verbosity)
 
virtual ~HGCalClusteringAlgoBase ()
 

Static Public Member Functions

static void fillPSetDescription (edm::ParameterSetDescription &iDesc)
 

Private Member Functions

void calculateDistanceToHigher (const TILE &lt, const unsigned int layerId, float delta)
 
void calculateLocalDensity (const TILE &lt, const unsigned int layerId, float delta)
 
void calculateLocalDensity (const TILE &lt, const unsigned int layerId, float delta, HGCalSiliconStrategy strategy)
 
void calculateLocalDensity (const TILE &lt, const unsigned int layerId, float delta, HGCalScintillatorStrategy strategy)
 
float distance (const TILE &lt, int cell1, int cell2, int layerId) const
 
float distance2 (const TILE &lt, int cell1, int cell2, int layerId) const
 
int findAndAssignClusters (const unsigned int layerId, float delta)
 
void prepareDataStructures (const unsigned int layerId)
 

Private Attributes

std::vector< CellsOnLayercells_
 
std::vector< double > dEdXweights_
 
int deltasi_index_regemfac_
 
bool dependSensor_
 
double ecut_
 
double fcPerEle_
 
std::vector< double > fcPerMip_
 
bool initialized_
 
double kappa_
 
unsigned maxNumberOfThickIndices_
 
double noiseMip_
 
std::vector< double > nonAgedNoises_
 
std::vector< int > numberOfClustersPerLayer_
 
float outlierDeltaFactor_ = 2.f
 
double positionDeltaRho2_
 
double sciThicknessCorrection_
 
std::vector< double > thicknessCorrection_
 
std::vector< std::vector< double > > thresholds_
 
std::vector< double > thresholdW0_
 
bool use2x2_
 
std::vector< std::vector< double > > v_sigmaNoise_
 
std::vector< double > vecDeltas_
 

Additional Inherited Members

- Public Attributes inherited from HGCalClusteringAlgoBase
unsigned int firstLayerBH_
 
bool isNose_
 
unsigned int lastLayerEE_
 
unsigned int lastLayerFH_
 
unsigned int maxlayer_
 
int scintMaxIphi_
 
- Protected Attributes inherited from HGCalClusteringAlgoBase
reco::CaloCluster::AlgoId algoId_
 
edm::ESGetToken< CaloGeometry, CaloGeometryRecordcaloGeomToken_
 
std::vector< reco::BasicClusterclusters_v_
 
hgcal::RecHitTools rhtools_
 
VerbosityLevel verbosity_
 

Detailed Description

template<typename TILE, typename STRATEGY>
class HGCalCLUEAlgoT< TILE, STRATEGY >

Definition at line 35 of file HGCalCLUEAlgo.h.

Member Typedef Documentation

◆ Point

template<typename TILE , typename STRATEGY >
typedef math::XYZPoint HGCalCLUEAlgoT< TILE, STRATEGY >::Point

point in the space

Definition at line 128 of file HGCalCLUEAlgo.h.

Constructor & Destructor Documentation

◆ HGCalCLUEAlgoT()

template<typename TILE , typename STRATEGY >
HGCalCLUEAlgoT< TILE, STRATEGY >::HGCalCLUEAlgoT ( const edm::ParameterSet ps)
inline

Definition at line 37 of file HGCalCLUEAlgo.h.

References edm::ParameterSet::getParameter(), and AlCaHLTBitMon_QueryRunRegistry::string.

39  (HGCalClusteringAlgoBase::VerbosityLevel)ps.getUntrackedParameter<unsigned int>("verbosity", 3),
41  vecDeltas_(ps.getParameter<std::vector<double>>("deltac")),
42  kappa_(ps.getParameter<double>("kappa")),
43  ecut_(ps.getParameter<double>("ecut")),
44  dependSensor_(ps.getParameter<bool>("dependSensor")),
45  dEdXweights_(ps.getParameter<std::vector<double>>("dEdXweights")),
46  thicknessCorrection_(ps.getParameter<std::vector<double>>("thicknessCorrection")),
47  sciThicknessCorrection_(ps.getParameter<double>("sciThicknessCorrection")),
48  deltasi_index_regemfac_(ps.getParameter<int>("deltasi_index_regemfac")),
49  maxNumberOfThickIndices_(ps.getParameter<unsigned>("maxNumberOfThickIndices")),
50  fcPerMip_(ps.getParameter<std::vector<double>>("fcPerMip")),
51  fcPerEle_(ps.getParameter<double>("fcPerEle")),
52  nonAgedNoises_(ps.getParameter<std::vector<double>>("noises")),
53  noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("noise_MIP")),
54  thresholdW0_(ps.getParameter<std::vector<double>>("thresholdW0")),
55  positionDeltaRho2_(ps.getParameter<double>("positionDeltaRho2")),
56  use2x2_(ps.getParameter<bool>("use2x2")),
57  initialized_(false) {
58 #if DEBUG_CLUSTERS_ALPAKA
59  moduleType_ = ps.getParameter<std::string>("type");
60 #endif
61  }
std::vector< double > nonAgedNoises_
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
std::vector< double > vecDeltas_
int deltasi_index_regemfac_
T getUntrackedParameter(std::string const &, T const &) const
HGCalClusteringAlgoBase(VerbosityLevel v, reco::CaloCluster::AlgoId algo)
double positionDeltaRho2_
double sciThicknessCorrection_
std::vector< double > thicknessCorrection_
std::vector< double > fcPerMip_
std::vector< double > dEdXweights_
unsigned maxNumberOfThickIndices_
std::vector< double > thresholdW0_

◆ ~HGCalCLUEAlgoT()

template<typename TILE , typename STRATEGY >
HGCalCLUEAlgoT< TILE, STRATEGY >::~HGCalCLUEAlgoT ( )
inlineoverride

Definition at line 63 of file HGCalCLUEAlgo.h.

63 {}

Member Function Documentation

◆ calculateDistanceToHigher()

template<typename TILE , typename STRATEGY >
void HGCalCLUEAlgoT< T, STRATEGY >::calculateDistanceToHigher ( const TILE &  lt,
const unsigned int  layerId,
float  delta 
)
private

Definition at line 347 of file HGCalCLUEAlgo.cc.

References ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), dumpMFGeometry_cfg::delta, mps_fire::i, dqmiolumiharvest::j, LogDebug, SiStripPI::max, allConversions_cfi::maxDelta, FastTimerService_cff::range, and mathSSE::sqrt().

347  {
348  auto& cellsOnLayer = cells_[layerId];
349  unsigned int numberOfCells = cellsOnLayer.detid.size();
350 
351  for (unsigned int i = 0; i < numberOfCells; i++) {
352  // initialize delta and nearest higher for i
354  float i_delta = maxDelta;
355  int i_nearestHigher = -1;
356  float rho_max = 0.f;
358  std::array<int, 4> search_box = lt.searchBox(cellsOnLayer.dim1[i] - range,
359  cellsOnLayer.dim1[i] + range,
360  cellsOnLayer.dim2[i] - range,
361  cellsOnLayer.dim2[i] + range);
362  // loop over all bins in the search box
363  for (int dim1Bin = search_box[0]; dim1Bin < search_box[1] + 1; ++dim1Bin) {
364  for (int dim2Bin = search_box[2]; dim2Bin < search_box[3] + 1; ++dim2Bin) {
365  // get the id of this bin
366  size_t binId = lt.getGlobalBinByBin(dim1Bin, dim2Bin);
367  if constexpr (std::is_same_v<STRATEGY, HGCalScintillatorStrategy>)
368  binId = lt.getGlobalBinByBin(dim1Bin, (dim2Bin % T::type::nRows));
369  // get the size of this bin
370  size_t binSize = lt[binId].size();
371 
372  // loop over all hits in this bin
373  for (unsigned int j = 0; j < binSize; j++) {
374  unsigned int otherId = lt[binId][j];
375  float dist = distance2(lt, i, otherId, layerId);
376  bool foundHigher =
377  (cellsOnLayer.rho[otherId] > cellsOnLayer.rho[i]) ||
378  (cellsOnLayer.rho[otherId] == cellsOnLayer.rho[i] && cellsOnLayer.detid[otherId] > cellsOnLayer.detid[i]);
379  if (!foundHigher) {
380  continue;
381  }
382  if ((dist < i_delta) || ((dist == i_delta) && (cellsOnLayer.rho[otherId] > rho_max)) ||
383  ((dist == i_delta) && (cellsOnLayer.rho[otherId] == rho_max) &&
384  (cellsOnLayer.detid[otherId] > cellsOnLayer.detid[i]))) {
385  rho_max = cellsOnLayer.rho[otherId];
386  i_delta = dist;
387  i_nearestHigher = otherId;
388  }
389  }
390  }
391  }
392  bool foundNearestHigherInSearchBox = (i_delta != maxDelta);
393  if (foundNearestHigherInSearchBox) {
394  cellsOnLayer.delta[i] = std::sqrt(i_delta);
395  cellsOnLayer.nearestHigher[i] = i_nearestHigher;
396  } else {
397  // otherwise delta is guaranteed to be larger outlierDeltaFactor_*delta_c
398  // we can safely maximize delta to be maxDelta
399  cellsOnLayer.delta[i] = maxDelta;
400  cellsOnLayer.nearestHigher[i] = -1;
401  }
402 
403  LogDebug("HGCalCLUEAlgo") << "Debugging calculateDistanceToHigher: \n"
404  << " cell: " << i << " eta: " << cellsOnLayer.dim1[i] << " phi: " << cellsOnLayer.dim2[i]
405  << " energy: " << cellsOnLayer.weight[i] << " density: " << cellsOnLayer.rho[i]
406  << " nearest higher: " << cellsOnLayer.nearestHigher[i]
407  << " distance: " << cellsOnLayer.delta[i] << "\n";
408  }
409 }
float outlierDeltaFactor_
T sqrt(T t)
Definition: SSEVec.h:23
std::vector< CellsOnLayer > cells_
float distance2(const TILE &lt, int cell1, int cell2, int layerId) const
#define LogDebug(id)

◆ calculateLocalDensity() [1/3]

template<typename TILE , typename STRATEGY >
void HGCalCLUEAlgoT< TILE, STRATEGY >::calculateLocalDensity ( const TILE &  lt,
const unsigned int  layerId,
float  delta 
)
private

◆ calculateLocalDensity() [2/3]

template<typename TILE , typename STRATEGY >
void HGCalCLUEAlgoT< TILE, STRATEGY >::calculateLocalDensity ( const TILE &  lt,
const unsigned int  layerId,
float  delta,
HGCalSiliconStrategy  strategy 
)
private

◆ calculateLocalDensity() [3/3]

template<typename TILE , typename STRATEGY >
void HGCalCLUEAlgoT< TILE, STRATEGY >::calculateLocalDensity ( const TILE &  lt,
const unsigned int  layerId,
float  delta,
HGCalScintillatorStrategy  strategy 
)
private

◆ computeThreshold()

template<typename T , typename STRATEGY >
void HGCalCLUEAlgoT< T, STRATEGY >::computeThreshold ( )

Definition at line 457 of file HGCalCLUEAlgo.cc.

References LogDebug.

457  {
458  // To support the TDR geometry and also the post-TDR one (v9 onwards), we
459  // need to change the logic of the vectors containing signal to noise and
460  // thresholds. The first 3 indices will keep on addressing the different
461  // thicknesses of the Silicon detectors in CE_E , the next 3 indices will address
462  // the thicknesses of the Silicon detectors in CE_H, while the last one, number 6 (the
463  // seventh) will address the Scintillators. This change will support both
464  // geometries at the same time.
465 
466  if (initialized_)
467  return; // only need to calculate thresholds once
468 
469  initialized_ = true;
470 
471  std::vector<double> dummy;
472 
473  dummy.resize(maxNumberOfThickIndices_ + !isNose_, 0); // +1 to accomodate for the Scintillators
474  thresholds_.resize(maxlayer_, dummy);
475  v_sigmaNoise_.resize(maxlayer_, dummy);
476 
477  for (unsigned ilayer = 1; ilayer <= maxlayer_; ++ilayer) {
478  for (unsigned ithick = 0; ithick < maxNumberOfThickIndices_; ++ithick) {
479  float sigmaNoise = 0.001f * fcPerEle_ * nonAgedNoises_[ithick] * dEdXweights_[ilayer] /
480  (fcPerMip_[ithick] * thicknessCorrection_[ithick]);
481  thresholds_[ilayer - 1][ithick] = sigmaNoise * ecut_;
482  v_sigmaNoise_[ilayer - 1][ithick] = sigmaNoise;
483  LogDebug("HGCalCLUEAlgo") << "ilayer: " << ilayer << " nonAgedNoises: " << nonAgedNoises_[ithick]
484  << " fcPerEle: " << fcPerEle_ << " fcPerMip: " << fcPerMip_[ithick]
485  << " noiseMip: " << fcPerEle_ * nonAgedNoises_[ithick] / fcPerMip_[ithick]
486  << " sigmaNoise: " << sigmaNoise << "\n";
487  }
488 
489  if (!isNose_) {
490  float scintillators_sigmaNoise = 0.001f * noiseMip_ * dEdXweights_[ilayer] / sciThicknessCorrection_;
491  thresholds_[ilayer - 1][maxNumberOfThickIndices_] = ecut_ * scintillators_sigmaNoise;
492  v_sigmaNoise_[ilayer - 1][maxNumberOfThickIndices_] = scintillators_sigmaNoise;
493  LogDebug("HGCalCLUEAlgo") << "ilayer: " << ilayer << " noiseMip: " << noiseMip_
494  << " scintillators_sigmaNoise: " << scintillators_sigmaNoise << "\n";
495  }
496  }
497 }
std::vector< double > nonAgedNoises_
double sciThicknessCorrection_
std::vector< double > thicknessCorrection_
std::vector< double > fcPerMip_
std::vector< double > dEdXweights_
std::vector< std::vector< double > > v_sigmaNoise_
unsigned maxNumberOfThickIndices_
std::vector< std::vector< double > > thresholds_
#define LogDebug(id)

◆ distance()

template<typename TILE , typename STRATEGY >
float HGCalCLUEAlgoT< TILE, STRATEGY >::distance ( const TILE &  lt,
int  cell1,
int  cell2,
int  layerId 
) const
inlineprivate

Definition at line 222 of file HGCalCLUEAlgo.h.

References hgcalTopologyTester_cfi::cell1, hgcalTopologyTester_cfi::cell2, HGCalCLUEAlgoT< TILE, STRATEGY >::cells_, and mathSSE::sqrt().

222  { // 2-d distance on the layer (x-y)
223  return std::sqrt(lt.distance2(cells_[layerId].dim1[cell1],
224  cells_[layerId].dim2[cell1],
225  cells_[layerId].dim1[cell2],
226  cells_[layerId].dim2[cell2]));
227  }
T sqrt(T t)
Definition: SSEVec.h:23
std::vector< CellsOnLayer > cells_

◆ distance2()

template<typename TILE , typename STRATEGY >
float HGCalCLUEAlgoT< TILE, STRATEGY >::distance2 ( const TILE &  lt,
int  cell1,
int  cell2,
int  layerId 
) const
inlineprivate

Definition at line 215 of file HGCalCLUEAlgo.h.

References hgcalTopologyTester_cfi::cell1, hgcalTopologyTester_cfi::cell2, and HGCalCLUEAlgoT< TILE, STRATEGY >::cells_.

215  { // 2-d distance on the layer (x-y)
216  return (lt.distance2(cells_[layerId].dim1[cell1],
217  cells_[layerId].dim2[cell1],
218  cells_[layerId].dim1[cell2],
219  cells_[layerId].dim2[cell2]));
220  }
std::vector< CellsOnLayer > cells_

◆ fillPSetDescription()

template<typename TILE , typename STRATEGY >
static void HGCalCLUEAlgoT< TILE, STRATEGY >::fillPSetDescription ( edm::ParameterSetDescription iDesc)
inlinestatic

Definition at line 92 of file HGCalCLUEAlgo.h.

References edm::ParameterSetDescription::add(), edm::ParameterSetDescription::addUntracked(), and AlCaHLTBitMon_QueryRunRegistry::string.

92  {
93  iDesc.add<std::vector<double>>("thresholdW0", {2.9, 2.9, 2.9});
94  iDesc.add<double>("positionDeltaRho2", 1.69);
95  iDesc.add<std::vector<double>>("deltac",
96  {
97  1.3,
98  1.3,
99  1.3,
100  0.0315, // for scintillator
101  });
102  iDesc.add<bool>("dependSensor", true);
103  iDesc.add<double>("ecut", 3.0);
104  iDesc.add<double>("kappa", 9.0);
105  iDesc.addUntracked<unsigned int>("verbosity", 3);
106  iDesc.add<std::vector<double>>("dEdXweights", {});
107  iDesc.add<std::vector<double>>("thicknessCorrection", {});
108  iDesc.add<double>("sciThicknessCorrection", 0.9);
109  iDesc.add<int>("deltasi_index_regemfac", 3);
110  iDesc.add<unsigned>("maxNumberOfThickIndices", 6);
111  iDesc.add<std::vector<double>>("fcPerMip", {});
112  iDesc.add<double>("fcPerEle", 0.0);
113  iDesc.add<std::vector<double>>("noises", {});
114  edm::ParameterSetDescription descNestedNoiseMIP;
115  descNestedNoiseMIP.add<bool>("scaleByDose", false);
116  descNestedNoiseMIP.add<unsigned int>("scaleByDoseAlgo", 0);
117  descNestedNoiseMIP.add<double>("scaleByDoseFactor", 1.);
118  descNestedNoiseMIP.add<std::string>("doseMap", "");
119  descNestedNoiseMIP.add<std::string>("sipmMap", "");
120  descNestedNoiseMIP.add<double>("referenceIdark", -1);
121  descNestedNoiseMIP.add<double>("referenceXtalk", -1);
122  descNestedNoiseMIP.add<double>("noise_MIP", 1. / 100.);
123  iDesc.add<edm::ParameterSetDescription>("noiseMip", descNestedNoiseMIP);
124  iDesc.add<bool>("use2x2", true); // use 2x2 or 3x3 scenario for scint density calculation
125  }
ParameterDescriptionBase * addUntracked(U const &iLabel, T const &value)
ParameterDescriptionBase * add(U const &iLabel, T const &value)

◆ findAndAssignClusters()

template<typename T , typename STRATEGY >
int HGCalCLUEAlgoT< T, STRATEGY >::findAndAssignClusters ( const unsigned int  layerId,
float  delta 
)
private

Definition at line 412 of file HGCalCLUEAlgo.cc.

References dumpMFGeometry_cfg::delta, mps_fire::i, and dqmiolumiharvest::j.

412  {
413  // this is called once per layer and endcap...
414  // so when filling the cluster temporary vector of Hexels we resize each time
415  // by the number of clusters found. This is always equal to the number of
416  // cluster centers...
417  unsigned int nClustersOnLayer = 0;
418  auto& cellsOnLayer = cells_[layerId];
419  unsigned int numberOfCells = cellsOnLayer.detid.size();
420  std::vector<int> localStack;
421  // find cluster seeds and outlier
422  for (unsigned int i = 0; i < numberOfCells; i++) {
423  float rho_c = kappa_ * cellsOnLayer.sigmaNoise[i];
424  // initialize clusterIndex
425  cellsOnLayer.clusterIndex[i] = -1;
426  bool isSeed = (cellsOnLayer.delta[i] > delta) && (cellsOnLayer.rho[i] >= rho_c);
427  bool isOutlier = (cellsOnLayer.delta[i] > outlierDeltaFactor_ * delta) && (cellsOnLayer.rho[i] < rho_c);
428  if (isSeed) {
429  cellsOnLayer.clusterIndex[i] = nClustersOnLayer;
430  cellsOnLayer.isSeed[i] = true;
431  nClustersOnLayer++;
432  localStack.push_back(i);
433 
434  } else if (!isOutlier) {
435  cellsOnLayer.followers[cellsOnLayer.nearestHigher[i]].push_back(i);
436  }
437  }
438 
439  // need to pass clusterIndex to their followers
440  while (!localStack.empty()) {
441  int endStack = localStack.back();
442  auto& thisSeed = cellsOnLayer.followers[endStack];
443  localStack.pop_back();
444 
445  // loop over followers
446  for (int j : thisSeed) {
447  // pass id to a follower
448  cellsOnLayer.clusterIndex[j] = cellsOnLayer.clusterIndex[endStack];
449  // push this follower to localStack
450  localStack.push_back(j);
451  }
452  }
453  return nClustersOnLayer;
454 }
float outlierDeltaFactor_
std::vector< CellsOnLayer > cells_

◆ getClusters()

template<typename T , typename STRATEGY >
std::vector< reco::BasicCluster > HGCalCLUEAlgoT< T, STRATEGY >::getClusters ( bool  )
overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 137 of file HGCalCLUEAlgo.cc.

References haddnano::cl, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), d1, reco::CaloID::DET_HGCAL_ENDCAP, hcalRecHitTable_cff::energy, mps_fire::i, dqm-mbProfile::log, SiStripPI::max, SiStripPI::min, unpackBuffers-CaloStage1::offsets, position, and x.

137  {
138  std::vector<int> offsets(numberOfClustersPerLayer_.size(), 0);
139 
140  int maxClustersOnLayer = numberOfClustersPerLayer_[0];
141 
142  for (unsigned layerId = 1; layerId < offsets.size(); ++layerId) {
143  offsets[layerId] = offsets[layerId - 1] + numberOfClustersPerLayer_[layerId - 1];
144 
145  maxClustersOnLayer = std::max(maxClustersOnLayer, numberOfClustersPerLayer_[layerId]);
146  }
147 
148  auto totalNumberOfClusters = offsets.back() + numberOfClustersPerLayer_.back();
149  clusters_v_.resize(totalNumberOfClusters);
150  std::vector<std::vector<int>> cellsIdInCluster;
151  cellsIdInCluster.reserve(maxClustersOnLayer);
152 
153  for (unsigned int layerId = 0; layerId < 2 * maxlayer_ + 2; ++layerId) {
154  cellsIdInCluster.resize(numberOfClustersPerLayer_[layerId]);
155  auto& cellsOnLayer = cells_[layerId];
156  unsigned int numberOfCells = cellsOnLayer.detid.size();
157  auto firstClusterIdx = offsets[layerId];
158 
159  for (unsigned int i = 0; i < numberOfCells; ++i) {
160  auto clusterIndex = cellsOnLayer.clusterIndex[i];
161  if (clusterIndex != -1)
162  cellsIdInCluster[clusterIndex].push_back(i);
163  }
164 
165  std::vector<std::pair<DetId, float>> thisCluster;
166 
167  for (auto& cl : cellsIdInCluster) {
168  float maxEnergyValue = std::numeric_limits<float>::min();
169  int maxEnergyCellIndex = -1;
170  DetId maxEnergyDetId;
171  float energy = 0.f;
172  int seedDetId = -1;
173  float x = 0.f;
174  float y = 0.f;
175  float z = cellsOnLayer.layerDim3;
176  // TODO Felice: maybe use the seed for the position calculation
177  for (auto cellIdx : cl) {
178  energy += cellsOnLayer.weight[cellIdx];
179  if (cellsOnLayer.weight[cellIdx] > maxEnergyValue) {
180  maxEnergyValue = cellsOnLayer.weight[cellIdx];
181  maxEnergyCellIndex = cellIdx;
182  maxEnergyDetId = cellsOnLayer.detid[cellIdx];
183  }
184  thisCluster.emplace_back(cellsOnLayer.detid[cellIdx], 1.f);
185  if (cellsOnLayer.isSeed[cellIdx]) {
186  seedDetId = cellsOnLayer.detid[cellIdx];
187  }
188  }
189 
190  float total_weight_log = 0.f;
191  float total_weight = energy;
192 
193  if constexpr (std::is_same_v<STRATEGY, HGCalSiliconStrategy>) {
194  auto thick = rhtools_.getSiThickIndex(maxEnergyDetId);
195  for (auto cellIdx : cl) {
196  const float d1 = cellsOnLayer.dim1[cellIdx] - cellsOnLayer.dim1[maxEnergyCellIndex];
197  const float d2 = cellsOnLayer.dim2[cellIdx] - cellsOnLayer.dim2[maxEnergyCellIndex];
198  if ((d1 * d1 + d2 * d2) < positionDeltaRho2_) {
199  float Wi = std::max(thresholdW0_[thick] + std::log(cellsOnLayer.weight[cellIdx] / energy), 0.);
200  x += cellsOnLayer.dim1[cellIdx] * Wi;
201  y += cellsOnLayer.dim2[cellIdx] * Wi;
202  total_weight_log += Wi;
203  }
204  }
205  } else {
206  for (auto cellIdx : cl) {
207  auto position = rhtools_.getPosition(cellsOnLayer.detid[cellIdx]);
208  x += position.x() * cellsOnLayer.weight[cellIdx];
209  y += position.y() * cellsOnLayer.weight[cellIdx];
210  }
211  }
212 
213  if constexpr (std::is_same_v<STRATEGY, HGCalSiliconStrategy>) {
214  total_weight = total_weight_log;
215  }
216 
217  if (total_weight != 0.) {
218  float inv_tot_weight = 1.f / total_weight;
219  x *= inv_tot_weight;
220  y *= inv_tot_weight;
221  } else {
222  x = cellsOnLayer.dim1[maxEnergyCellIndex];
223  y = cellsOnLayer.dim2[maxEnergyCellIndex];
224  }
226 
227  auto globalClusterIndex = cellsOnLayer.clusterIndex[cl[0]] + firstClusterIdx;
228 
229  clusters_v_[globalClusterIndex] =
231  clusters_v_[globalClusterIndex].setSeed(seedDetId);
232  thisCluster.clear();
233  }
234 
235  cellsIdInCluster.clear();
236  }
237  return clusters_v_;
238 }
CaloCluster BasicCluster
std::vector< reco::BasicCluster > clusters_v_
double positionDeltaRho2_
GlobalPoint getPosition(const DetId &id) const
Definition: RecHitTools.cc:140
reco::CaloCluster::AlgoId algoId_
Definition: DetId.h:17
std::vector< int > numberOfClustersPerLayer_
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
std::vector< double > thresholdW0_
std::vector< CellsOnLayer > cells_
static int position[264][3]
Definition: ReadPGInfo.cc:289
static constexpr float d1
int getSiThickIndex(const DetId &) const
Definition: RecHitTools.cc:216

◆ getEventSetupPerAlgorithm()

template<typename T , typename STRATEGY >
void HGCalCLUEAlgoT< T, STRATEGY >::getEventSetupPerAlgorithm ( const edm::EventSetup es)
overridevirtual

Reimplemented from HGCalClusteringAlgoBase.

Definition at line 20 of file HGCalCLUEAlgo.cc.

20  {
21  cells_.clear();
23  cells_.resize(2 * (maxlayer_ + 1));
24  numberOfClustersPerLayer_.resize(2 * (maxlayer_ + 1), 0);
25 }
std::vector< int > numberOfClustersPerLayer_
std::vector< CellsOnLayer > cells_

◆ makeClusters()

template<typename T , typename STRATEGY >
void HGCalCLUEAlgoT< T, STRATEGY >::makeClusters ( )
overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 98 of file HGCalCLUEAlgo.cc.

References ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), dumpMFGeometry_cfg::delta, hitfit::delta_r(), hgcalUtils::DumpLegacySoA::dumpInfos(), mps_fire::i, and LogDebug.

98  {
99  // assign all hits in each layer to a cluster core
100  tbb::this_task_arena::isolate([&] {
101  tbb::parallel_for(size_t(0), size_t(2 * maxlayer_ + 2), [&](size_t i) {
103  T lt;
104  lt.clear();
105  lt.fill(cells_[i].dim1, cells_[i].dim2);
106 
107  float delta;
108  if constexpr (std::is_same_v<STRATEGY, HGCalSiliconStrategy>) {
109  // maximum search distance (critical distance) for local density calculation
110  float delta_c;
111  if (i % maxlayer_ < lastLayerEE_)
112  delta_c = vecDeltas_[0];
113  else if (i % maxlayer_ < (firstLayerBH_ - 1))
114  delta_c = vecDeltas_[1];
115  else
116  delta_c = vecDeltas_[2];
117  delta = delta_c;
118  } else {
119  float delta_r = vecDeltas_[3];
120  delta = delta_r;
121  }
122  LogDebug("HGCalCLUEAlgo") << "maxlayer: " << maxlayer_ << " lastLayerEE: " << lastLayerEE_
123  << " firstLayerBH: " << firstLayerBH_ << "\n";
124 
128  });
129  });
130 #if DEBUG_CLUSTERS_ALPAKA
131  hgcalUtils::DumpLegacySoA dumperLegacySoA;
132  dumperLegacySoA.dumpInfos(cells_, moduleType_);
133 #endif
134 }
std::vector< double > vecDeltas_
void calculateDistanceToHigher(const TILE &lt, const unsigned int layerId, float delta)
void dumpInfos(const T &cells, const std::string &moduleType) const
double delta_r(const Fourvec &a, const Fourvec &b)
Find the distance between two four-vectors in the two-dimensional space .
Definition: fourvec.cc:238
std::vector< int > numberOfClustersPerLayer_
void calculateLocalDensity(const TILE &lt, const unsigned int layerId, float delta)
std::vector< CellsOnLayer > cells_
void prepareDataStructures(const unsigned int layerId)
long double T
int findAndAssignClusters(const unsigned int layerId, float delta)
#define LogDebug(id)

◆ populate()

template<typename T , typename STRATEGY >
void HGCalCLUEAlgoT< T, STRATEGY >::populate ( const HGCRecHitCollection hits)
overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 28 of file HGCalCLUEAlgo.cc.

References ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), DetId::det(), CaloRecHit::detid(), CaloRecHit::energy(), DetId::HGCalHSi, HGCHEF, hfClusterShapes_cfi::hits, mps_fire::i, infinity, hltrates_dqm_sourceclient-live_cfg::offset, position, and DetId::subdetId().

28  {
29  // loop over all hits and create the Hexel structure, skip energies below ecut
30  if (dependSensor_) {
31  // for each layer and wafer calculate the thresholds (sigmaNoise and energy)
32  // once
34  }
35 
36  for (unsigned int i = 0; i < hits.size(); ++i) {
37  const HGCRecHit& hgrh = hits[i];
38  DetId detid = hgrh.detid();
39  unsigned int layerOnSide = (rhtools_.getLayerWithOffset(detid) - 1);
40 
41  // set sigmaNoise default value 1 to use kappa value directly in case of
42  // sensor-independent thresholds
43  float sigmaNoise = 1.f;
44  if (dependSensor_) {
45  int thickness_index = rhtools_.getSiThickIndex(detid);
46  if (thickness_index == -1)
47  thickness_index = maxNumberOfThickIndices_;
48 
49  double storedThreshold = thresholds_[layerOnSide][thickness_index];
50  if (detid.det() == DetId::HGCalHSi || detid.subdetId() == HGCHEF) {
51  storedThreshold = thresholds_[layerOnSide][thickness_index + deltasi_index_regemfac_];
52  }
53  sigmaNoise = v_sigmaNoise_[layerOnSide][thickness_index];
54 
55  if (hgrh.energy() < storedThreshold)
56  continue; // this sets the ZS threshold at ecut times the sigma noise
57  // for the sensor
58  }
59  if (!dependSensor_ && hgrh.energy() < ecut_)
60  continue;
62  int offset = ((rhtools_.zside(detid) + 1) >> 1) * maxlayer_;
63  int layer = layerOnSide + offset;
64  // setting the layer position only once per layer
66  cells_[layer].layerDim3 = position.z();
67 
68  cells_[layer].detid.emplace_back(detid);
69  if constexpr (std::is_same_v<STRATEGY, HGCalScintillatorStrategy>) {
70  cells_[layer].dim1.emplace_back(position.eta());
71  cells_[layer].dim2.emplace_back(position.phi());
72  } // else, isSilicon == true and eta phi values will not be used
73  else {
74  cells_[layer].dim1.emplace_back(position.x());
75  cells_[layer].dim2.emplace_back(position.y());
76  }
77  cells_[layer].weight.emplace_back(hgrh.energy());
78  cells_[layer].sigmaNoise.emplace_back(sigmaNoise);
79  }
80 }
void computeThreshold()
constexpr const DetId & detid() const
Definition: CaloRecHit.h:33
constexpr Detector det() const
get the detector field from this detid
Definition: DetId.h:46
int zside(const DetId &id) const
Definition: RecHitTools.cc:174
int deltasi_index_regemfac_
constexpr float energy() const
Definition: CaloRecHit.h:29
GlobalPoint getPosition(const DetId &id) const
Definition: RecHitTools.cc:140
const double infinity
constexpr int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:48
Definition: DetId.h:17
std::vector< std::vector< double > > v_sigmaNoise_
unsigned maxNumberOfThickIndices_
std::vector< CellsOnLayer > cells_
static int position[264][3]
Definition: ReadPGInfo.cc:289
std::vector< std::vector< double > > thresholds_
int getSiThickIndex(const DetId &) const
Definition: RecHitTools.cc:216
unsigned int getLayerWithOffset(const DetId &) const
Definition: RecHitTools.cc:381

◆ prepareDataStructures()

template<typename T , typename STRATEGY >
void HGCalCLUEAlgoT< T, STRATEGY >::prepareDataStructures ( const unsigned int  layerId)
private

Definition at line 83 of file HGCalCLUEAlgo.cc.

References f, and MainPageGenerator::l.

83  {
84  auto cellsSize = cells_[l].detid.size();
85  cells_[l].rho.resize(cellsSize, 0.f);
86  cells_[l].delta.resize(cellsSize, 9999999);
87  cells_[l].nearestHigher.resize(cellsSize, -1);
88  cells_[l].clusterIndex.resize(cellsSize, -1);
89  cells_[l].followers.resize(cellsSize);
90  cells_[l].isSeed.resize(cellsSize, false);
91 }
double f[11][100]
std::vector< CellsOnLayer > cells_

◆ reset()

template<typename TILE , typename STRATEGY >
void HGCalCLUEAlgoT< TILE, STRATEGY >::reset ( void  )
inlineoverridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 77 of file HGCalCLUEAlgo.h.

References hgcalTBTopologyTester_cfi::cells, HGCalCLUEAlgoT< TILE, STRATEGY >::cells_, haddnano::cl, HGCalClusteringAlgoBase::clusters_v_, and HGCalCLUEAlgoT< TILE, STRATEGY >::numberOfClustersPerLayer_.

77  {
78  clusters_v_.clear();
79  clusters_v_.shrink_to_fit();
80  for (auto& cl : numberOfClustersPerLayer_) {
81  cl = 0;
82  }
83 
84  for (auto& cells : cells_) {
85  cells.clear();
86  cells.shrink_to_fit();
87  }
88  }
std::vector< reco::BasicCluster > clusters_v_
std::vector< int > numberOfClustersPerLayer_
std::vector< CellsOnLayer > cells_

Member Data Documentation

◆ cells_

template<typename TILE , typename STRATEGY >
std::vector<CellsOnLayer> HGCalCLUEAlgoT< TILE, STRATEGY >::cells_
private

◆ dEdXweights_

template<typename TILE , typename STRATEGY >
std::vector<double> HGCalCLUEAlgoT< TILE, STRATEGY >::dEdXweights_
private

Definition at line 141 of file HGCalCLUEAlgo.h.

◆ deltasi_index_regemfac_

template<typename TILE , typename STRATEGY >
int HGCalCLUEAlgoT< TILE, STRATEGY >::deltasi_index_regemfac_
private

Definition at line 144 of file HGCalCLUEAlgo.h.

◆ dependSensor_

template<typename TILE , typename STRATEGY >
bool HGCalCLUEAlgoT< TILE, STRATEGY >::dependSensor_
private

Definition at line 140 of file HGCalCLUEAlgo.h.

◆ ecut_

template<typename TILE , typename STRATEGY >
double HGCalCLUEAlgoT< TILE, STRATEGY >::ecut_
private

Definition at line 136 of file HGCalCLUEAlgo.h.

◆ fcPerEle_

template<typename TILE , typename STRATEGY >
double HGCalCLUEAlgoT< TILE, STRATEGY >::fcPerEle_
private

Definition at line 147 of file HGCalCLUEAlgo.h.

◆ fcPerMip_

template<typename TILE , typename STRATEGY >
std::vector<double> HGCalCLUEAlgoT< TILE, STRATEGY >::fcPerMip_
private

Definition at line 146 of file HGCalCLUEAlgo.h.

◆ initialized_

template<typename TILE , typename STRATEGY >
bool HGCalCLUEAlgoT< TILE, STRATEGY >::initialized_
private

Definition at line 158 of file HGCalCLUEAlgo.h.

◆ kappa_

template<typename TILE , typename STRATEGY >
double HGCalCLUEAlgoT< TILE, STRATEGY >::kappa_
private

Definition at line 133 of file HGCalCLUEAlgo.h.

◆ maxNumberOfThickIndices_

template<typename TILE , typename STRATEGY >
unsigned HGCalCLUEAlgoT< TILE, STRATEGY >::maxNumberOfThickIndices_
private

Definition at line 145 of file HGCalCLUEAlgo.h.

◆ noiseMip_

template<typename TILE , typename STRATEGY >
double HGCalCLUEAlgoT< TILE, STRATEGY >::noiseMip_
private

Definition at line 149 of file HGCalCLUEAlgo.h.

◆ nonAgedNoises_

template<typename TILE , typename STRATEGY >
std::vector<double> HGCalCLUEAlgoT< TILE, STRATEGY >::nonAgedNoises_
private

Definition at line 148 of file HGCalCLUEAlgo.h.

◆ numberOfClustersPerLayer_

template<typename TILE , typename STRATEGY >
std::vector<int> HGCalCLUEAlgoT< TILE, STRATEGY >::numberOfClustersPerLayer_
private

Definition at line 209 of file HGCalCLUEAlgo.h.

Referenced by HGCalCLUEAlgoT< TILE, STRATEGY >::reset().

◆ outlierDeltaFactor_

template<typename TILE , typename STRATEGY >
float HGCalCLUEAlgoT< TILE, STRATEGY >::outlierDeltaFactor_ = 2.f
private

Definition at line 160 of file HGCalCLUEAlgo.h.

◆ positionDeltaRho2_

template<typename TILE , typename STRATEGY >
double HGCalCLUEAlgoT< TILE, STRATEGY >::positionDeltaRho2_
private

Definition at line 153 of file HGCalCLUEAlgo.h.

◆ sciThicknessCorrection_

template<typename TILE , typename STRATEGY >
double HGCalCLUEAlgoT< TILE, STRATEGY >::sciThicknessCorrection_
private

Definition at line 143 of file HGCalCLUEAlgo.h.

◆ thicknessCorrection_

template<typename TILE , typename STRATEGY >
std::vector<double> HGCalCLUEAlgoT< TILE, STRATEGY >::thicknessCorrection_
private

Definition at line 142 of file HGCalCLUEAlgo.h.

◆ thresholds_

template<typename TILE , typename STRATEGY >
std::vector<std::vector<double> > HGCalCLUEAlgoT< TILE, STRATEGY >::thresholds_
private

Definition at line 150 of file HGCalCLUEAlgo.h.

◆ thresholdW0_

template<typename TILE , typename STRATEGY >
std::vector<double> HGCalCLUEAlgoT< TILE, STRATEGY >::thresholdW0_
private

Definition at line 152 of file HGCalCLUEAlgo.h.

◆ use2x2_

template<typename TILE , typename STRATEGY >
bool HGCalCLUEAlgoT< TILE, STRATEGY >::use2x2_
private

Definition at line 155 of file HGCalCLUEAlgo.h.

◆ v_sigmaNoise_

template<typename TILE , typename STRATEGY >
std::vector<std::vector<double> > HGCalCLUEAlgoT< TILE, STRATEGY >::v_sigmaNoise_
private

Definition at line 151 of file HGCalCLUEAlgo.h.

◆ vecDeltas_

template<typename TILE , typename STRATEGY >
std::vector<double> HGCalCLUEAlgoT< TILE, STRATEGY >::vecDeltas_
private

Definition at line 132 of file HGCalCLUEAlgo.h.