#include <HGCalCLUEAlgo.h>

Inheritance diagram for HGCalCLUEAlgoT< TILE >:

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::BasicCluster >	getClusters (bool) override

Density	getDensity () override

void	getEventSetupPerAlgorithm (const edm::EventSetup &es) override

	HGCalCLUEAlgoT (const edm::ParameterSet &ps, edm::ConsumesCollector iC)

void	makeClusters () override

void	populate (const HGCRecHitCollection &hits) override

void	reset () override

	~HGCalCLUEAlgoT () override

Public Member Functions inherited from HGCalClusteringAlgoBase
void	getEventSetup (const edm::EventSetup &es)

	HGCalClusteringAlgoBase (VerbosityLevel v, reco::CaloCluster::AlgoId algo, edm::ConsumesCollector iC)

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_c, float delta_r)

void	calculateLocalDensity (const TILE &lt, const unsigned int layerId, float delta_c, float delta_r)

math::XYZPoint	calculatePosition (const std::vector< int > &v, const unsigned int layerId) const

float	distance (int cell1, int cell2, int layerId, bool isEtaPhi) const

float	distance2 (int cell1, int cell2, int layerId, bool isEtaPhi) const

int	findAndAssignClusters (const unsigned int layerId, float delta_c, float delta_r)

void	prepareDataStructures (const unsigned int layerId)

void	setDensity (const unsigned int layerId)

Private Attributes
std::vector< CellsOnLayer >	cells_

std::vector< double >	dEdXweights_

int	deltasi_index_regemfac_

Density	density_

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

const 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, CaloGeometryRecord >	caloGeomToken_

std::vector< reco::BasicCluster >	clusters_v_

hgcal::RecHitTools	rhtools_

VerbosityLevel	verbosity_

Detailed Description

template<typename TILE>
class HGCalCLUEAlgoT< TILE >

Definition at line 30 of file HGCalCLUEAlgo.h.

Member Typedef Documentation

◆ Point

template<typename TILE >

typedef math::XYZPoint HGCalCLUEAlgoT< TILE >::Point

point in the space

Definition at line 120 of file HGCalCLUEAlgo.h.

Constructor & Destructor Documentation

◆ HGCalCLUEAlgoT()

template<typename TILE >

HGCalCLUEAlgoT< TILE >::HGCalCLUEAlgoT	(	const edm::ParameterSet &	ps,
		edm::ConsumesCollector	iC
	)

inline

Definition at line 32 of file HGCalCLUEAlgo.h.

       : HGCalClusteringAlgoBase(
             (HGCalClusteringAlgoBase::VerbosityLevel)ps.getUntrackedParameter<unsigned int>("verbosity", 3),
             reco::CaloCluster::undefined,
             iC),
         thresholdW0_(ps.getParameter<std::vector<double>>("thresholdW0")),
         positionDeltaRho2_(ps.getParameter<double>("positionDeltaRho2")),
         vecDeltas_(ps.getParameter<std::vector<double>>("deltac")),
         kappa_(ps.getParameter<double>("kappa")),
         ecut_(ps.getParameter<double>("ecut")),
         dependSensor_(ps.getParameter<bool>("dependSensor")),
         dEdXweights_(ps.getParameter<std::vector<double>>("dEdXweights")),
         thicknessCorrection_(ps.getParameter<std::vector<double>>("thicknessCorrection")),
         sciThicknessCorrection_(ps.getParameter<double>("sciThicknessCorrection")),
         deltasi_index_regemfac_(ps.getParameter<int>("deltasi_index_regemfac")),
         maxNumberOfThickIndices_(ps.getParameter<unsigned>("maxNumberOfThickIndices")),
         fcPerMip_(ps.getParameter<std::vector<double>>("fcPerMip")),
         fcPerEle_(ps.getParameter<double>("fcPerEle")),
         nonAgedNoises_(ps.getParameter<std::vector<double>>("noises")),
         noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("noise_MIP")),
         use2x2_(ps.getParameter<bool>("use2x2")),
         initialized_(false) {}

◆ ~HGCalCLUEAlgoT()

template<typename TILE >

HGCalCLUEAlgoT< TILE >::~HGCalCLUEAlgoT ( )

inlineoverride

Definition at line 55 of file HGCalCLUEAlgo.h.

55 {}

Member Function Documentation

◆ calculateDistanceToHigher()

template<typename TILE >

void HGCalCLUEAlgoT< T >::calculateDistanceToHigher	(	const TILE &	lt,
		const unsigned int	layerId,
		float	delta_c,
		float	delta_r
	)

private

Definition at line 345 of file HGCalCLUEAlgo.cc.

                                                                  {
   auto& cellsOnLayer = cells_[layerId];
   unsigned int numberOfCells = cellsOnLayer.detid.size();
   bool isOnlySi(false);
   if (rhtools_.isOnlySilicon(layerId))
     isOnlySi = true;
  
   for (unsigned int i = 0; i < numberOfCells; i++) {
     bool isSi = isOnlySi || cellsOnLayer.isSi[i];
     // initialize delta and nearest higher for i
     float maxDelta = std::numeric_limits<float>::max();
     float i_delta = maxDelta;
     int i_nearestHigher = -1;
     if (isSi) {
       float delta = delta_c;
       // get search box for ith hit
       // guarantee to cover a range "outlierDeltaFactor_*delta_c"
       auto range = outlierDeltaFactor_ * delta;
       std::array<int, 4> search_box = lt.searchBox(
           cellsOnLayer.x[i] - range, cellsOnLayer.x[i] + range, cellsOnLayer.y[i] - range, cellsOnLayer.y[i] + range);
       // loop over all bins in the search box
       for (int xBin = search_box[0]; xBin < search_box[1] + 1; ++xBin) {
         for (int yBin = search_box[2]; yBin < search_box[3] + 1; ++yBin) {
           // get the id of this bin
           size_t binId = lt.getGlobalBinByBin(xBin, yBin);
           // get the size of this bin
           size_t binSize = lt[binId].size();
  
           // loop over all hits in this bin
           for (unsigned int j = 0; j < binSize; j++) {
             unsigned int otherId = lt[binId][j];
             bool otherSi = isOnlySi || cellsOnLayer.isSi[otherId];
             if (otherSi) {  //silicon cells cannot talk to scintillator cells
               float dist = distance(i, otherId, layerId, false);
               bool foundHigher = (cellsOnLayer.rho[otherId] > cellsOnLayer.rho[i]) ||
                                  (cellsOnLayer.rho[otherId] == cellsOnLayer.rho[i] &&
                                   cellsOnLayer.detid[otherId] > cellsOnLayer.detid[i]);
               // if dist == i_delta, then last comer being the nearest higher
               if (foundHigher && dist <= i_delta) {
                 // update i_delta
                 i_delta = dist;
                 // update i_nearestHigher
                 i_nearestHigher = otherId;
               }
             }
           }
         }
       }
  
       bool foundNearestHigherInSearchBox = (i_delta != maxDelta);
       if (foundNearestHigherInSearchBox) {
         cellsOnLayer.delta[i] = i_delta;
         cellsOnLayer.nearestHigher[i] = i_nearestHigher;
       } else {
         // otherwise delta is guaranteed to be larger outlierDeltaFactor_*delta_c
         // we can safely maximize delta to be maxDelta
         cellsOnLayer.delta[i] = maxDelta;
         cellsOnLayer.nearestHigher[i] = -1;
       }
     } else {
       //similar to silicon
       float delta = delta_r;
       auto range = outlierDeltaFactor_ * delta;
       std::array<int, 4> search_box = lt.searchBoxEtaPhi(cellsOnLayer.eta[i] - range,
                                                          cellsOnLayer.eta[i] + range,
                                                          cellsOnLayer.phi[i] - range,
                                                          cellsOnLayer.phi[i] + range);
       // loop over all bins in the search box
       for (int xBin = search_box[0]; xBin < search_box[1] + 1; ++xBin) {
         for (int yBin = search_box[2]; yBin < search_box[3] + 1; ++yBin) {
           // get the id of this bin
           size_t binId = lt.getGlobalBinByBinEtaPhi(xBin, yBin);
           // get the size of this bin
           size_t binSize = lt[binId].size();
  
           // loop over all hits in this bin
           for (unsigned int j = 0; j < binSize; j++) {
             unsigned int otherId = lt[binId][j];
             if (!cellsOnLayer.isSi[otherId]) {  //scintillator cells cannot talk to silicon cells
               float dist = distance(i, otherId, layerId, true);
               bool foundHigher = (cellsOnLayer.rho[otherId] > cellsOnLayer.rho[i]) ||
                                  (cellsOnLayer.rho[otherId] == cellsOnLayer.rho[i] &&
                                   cellsOnLayer.detid[otherId] > cellsOnLayer.detid[i]);
               // if dist == i_delta, then last comer being the nearest higher
               if (foundHigher && dist <= i_delta) {
                 // update i_delta
                 i_delta = dist;
                 // update i_nearestHigher
                 i_nearestHigher = otherId;
               }
             }
           }
         }
       }
  
       bool foundNearestHigherInSearchBox = (i_delta != maxDelta);
       if (foundNearestHigherInSearchBox) {
         cellsOnLayer.delta[i] = i_delta;
         cellsOnLayer.nearestHigher[i] = i_nearestHigher;
       } else {
         // otherwise delta is guaranteed to be larger outlierDeltaFactor_*delta_c
         // we can safely maximize delta to be maxDelta
         cellsOnLayer.delta[i] = maxDelta;
         cellsOnLayer.nearestHigher[i] = -1;
       }
     }
     LogDebug("HGCalCLUEAlgo") << "Debugging calculateDistanceToHigher: \n"
                               << "  cell: " << i << " isSilicon: " << cellsOnLayer.isSi[i]
                               << " eta: " << cellsOnLayer.eta[i] << " phi: " << cellsOnLayer.phi[i]
                               << " energy: " << cellsOnLayer.weight[i] << " density: " << cellsOnLayer.rho[i]
                               << " nearest higher: " << cellsOnLayer.nearestHigher[i]
                               << " distance: " << cellsOnLayer.delta[i] << "\n";
   }
 }

References dumpMFGeometry_cfg::delta, hitfit::delta_r(), HLT_FULL_cff::distance, mps_fire::i, dqmiolumiharvest::j, LogDebug, SiStripPI::max, allConversions_cfi::maxDelta, FastTimerService_cff::range, photonAnalyzer_cfi::xBin, and photonAnalyzer_cfi::yBin.

◆ calculateLocalDensity()

template<typename TILE >

void HGCalCLUEAlgoT< T >::calculateLocalDensity	(	const TILE &	lt,
		const unsigned int	layerId,
		float	delta_c,
		float	delta_r
	)

private

Definition at line 247 of file HGCalCLUEAlgo.cc.

                                                                                                                    {
   auto& cellsOnLayer = cells_[layerId];
   unsigned int numberOfCells = cellsOnLayer.detid.size();
   bool isOnlySi(false);
   if (rhtools_.isOnlySilicon(layerId))
     isOnlySi = true;
  
   for (unsigned int i = 0; i < numberOfCells; i++) {
     bool isSi = isOnlySi || cellsOnLayer.isSi[i];
     if (isSi) {
       float delta = delta_c;
       std::array<int, 4> search_box = lt.searchBox(
           cellsOnLayer.x[i] - delta, cellsOnLayer.x[i] + delta, cellsOnLayer.y[i] - delta, cellsOnLayer.y[i] + delta);
  
       for (int xBin = search_box[0]; xBin < search_box[1] + 1; ++xBin) {
         for (int yBin = search_box[2]; yBin < search_box[3] + 1; ++yBin) {
           int binId = lt.getGlobalBinByBin(xBin, yBin);
           size_t binSize = lt[binId].size();
  
           for (unsigned int j = 0; j < binSize; j++) {
             unsigned int otherId = lt[binId][j];
             bool otherSi = isOnlySi || cellsOnLayer.isSi[otherId];
             if (otherSi) {  //silicon cells cannot talk to scintillator cells
               if (distance(i, otherId, layerId, false) < delta) {
                 cellsOnLayer.rho[i] += (i == otherId ? 1.f : 0.5f) * cellsOnLayer.weight[otherId];
               }
             }
           }
         }
       }
     } else {
       float delta = delta_r;
       std::array<int, 4> search_box = lt.searchBoxEtaPhi(cellsOnLayer.eta[i] - delta,
                                                          cellsOnLayer.eta[i] + delta,
                                                          cellsOnLayer.phi[i] - delta,
                                                          cellsOnLayer.phi[i] + delta);
       cellsOnLayer.rho[i] += cellsOnLayer.weight[i];
       float northeast(0), northwest(0), southeast(0), southwest(0), all(0);
  
       for (int etaBin = search_box[0]; etaBin < search_box[1] + 1; ++etaBin) {
         for (int phiBin = search_box[2]; phiBin < search_box[3] + 1; ++phiBin) {
           int binId = lt.getGlobalBinByBinEtaPhi(etaBin, phiBin);
           size_t binSize = lt[binId].size();
  
           for (unsigned int j = 0; j < binSize; j++) {
             unsigned int otherId = lt[binId][j];
             if (!cellsOnLayer.isSi[otherId]) {  //scintillator cells cannot talk to silicon cells
               if (distance(i, otherId, layerId, true) < delta) {
                 int iPhi = HGCScintillatorDetId(cellsOnLayer.detid[i]).iphi();
                 int otherIPhi = HGCScintillatorDetId(cellsOnLayer.detid[otherId]).iphi();
                 int iEta = HGCScintillatorDetId(cellsOnLayer.detid[i]).ieta();
                 int otherIEta = HGCScintillatorDetId(cellsOnLayer.detid[otherId]).ieta();
                 int dIPhi = otherIPhi - iPhi;
                 dIPhi += abs(dIPhi) < 2 ? 0
                          : dIPhi < 0    ? scintMaxIphi_
                                         : -scintMaxIphi_;  // cells with iPhi=288 and iPhi=1 should be neiboring cells
                 int dIEta = otherIEta - iEta;
                 LogDebug("HGCalCLUEAlgo") << "  Debugging calculateLocalDensity for Scintillator: \n"
                                           << "    cell: " << otherId << " energy: " << cellsOnLayer.weight[otherId]
                                           << " otherIPhi: " << otherIPhi << " iPhi: " << iPhi
                                           << " otherIEta: " << otherIEta << " iEta: " << iEta << "\n";
  
                 if (otherId != i) {
                   auto neighborCellContribution = 0.5f * cellsOnLayer.weight[otherId];
                   all += neighborCellContribution;
                   if (dIPhi >= 0 && dIEta >= 0)
                     northeast += neighborCellContribution;
                   if (dIPhi <= 0 && dIEta >= 0)
                     southeast += neighborCellContribution;
                   if (dIPhi >= 0 && dIEta <= 0)
                     northwest += neighborCellContribution;
                   if (dIPhi <= 0 && dIEta <= 0)
                     southwest += neighborCellContribution;
                 }
                 LogDebug("HGCalCLUEAlgo") << "  Debugging calculateLocalDensity for Scintillator: \n"
                                           << "    northeast: " << northeast << " southeast: " << southeast
                                           << " northwest: " << northwest << " southwest: " << southwest << "\n";
               }
             }
           }
         }
       }
       float neighborsval = (std::max(northeast, northwest) > std::max(southeast, southwest))
                                ? std::max(northeast, northwest)
                                : std::max(southeast, southwest);
       if (use2x2_)
         cellsOnLayer.rho[i] += neighborsval;
       else
         cellsOnLayer.rho[i] += all;
     }
     LogDebug("HGCalCLUEAlgo") << "Debugging calculateLocalDensity: \n"
                               << "  cell: " << i << " isSilicon: " << cellsOnLayer.isSi[i]
                               << " eta: " << cellsOnLayer.eta[i] << " phi: " << cellsOnLayer.phi[i]
                               << " energy: " << cellsOnLayer.weight[i] << " density: " << cellsOnLayer.rho[i] << "\n";
   }
 }

References funct::abs(), python.cmstools::all(), dumpMFGeometry_cfg::delta, hitfit::delta_r(), HLT_FULL_cff::distance, etaBin(), mps_fire::i, HGCScintillatorDetId::ieta(), L1TowerCalibrationProducer_cfi::iEta, HGCScintillatorDetId::iphi(), dqmiolumiharvest::j, LogDebug, SiStripPI::max, BeamMonitor_cff::phiBin, photonAnalyzer_cfi::xBin, and photonAnalyzer_cfi::yBin.

◆ calculatePosition()

template<typename T >

math::XYZPoint HGCalCLUEAlgoT< T >::calculatePosition	(	const std::vector< int > &	v,
		const unsigned int	layerId
	)		const

private

Definition at line 187 of file HGCalCLUEAlgo.cc.

                                                                                                            {
   float total_weight = 0.f;
   float x = 0.f;
   float y = 0.f;
  
   unsigned int maxEnergyIndex = 0;
   float maxEnergyValue = 0.f;
  
   auto& cellsOnLayer = cells_[layerId];
  
   // loop over hits in cluster candidate
   // determining the maximum energy hit
   for (auto i : v) {
     total_weight += cellsOnLayer.weight[i];
     if (cellsOnLayer.weight[i] > maxEnergyValue) {
       maxEnergyValue = cellsOnLayer.weight[i];
       maxEnergyIndex = i;
     }
   }
  
   // Si cell or Scintillator. Used to set approach and parameters
   auto thick = rhtools_.getSiThickIndex(cellsOnLayer.detid[maxEnergyIndex]);
   bool isSiliconCell = (thick != -1);
  
   // TODO: this is recomputing everything twice and overwriting the position with log weighting position
   if (isSiliconCell) {
     float total_weight_log = 0.f;
     float x_log = 0.f;
     float y_log = 0.f;
     for (auto i : v) {
       //for silicon only just use 1+6 cells = 1.3cm for all thicknesses
       if (distance2(i, maxEnergyIndex, layerId, false) > positionDeltaRho2_)
         continue;
       float rhEnergy = cellsOnLayer.weight[i];
       float Wi = std::max(thresholdW0_[thick] + std::log(rhEnergy / total_weight), 0.);
       x_log += cellsOnLayer.x[i] * Wi;
       y_log += cellsOnLayer.y[i] * Wi;
       total_weight_log += Wi;
     }
  
     total_weight = total_weight_log;
     x = x_log;
     y = y_log;
   } else {
     for (auto i : v) {
       float rhEnergy = cellsOnLayer.weight[i];
  
       x += cellsOnLayer.x[i] * rhEnergy;
       y += cellsOnLayer.y[i] * rhEnergy;
     }
   }
   if (total_weight != 0.) {
     float inv_tot_weight = 1.f / total_weight;
     return math::XYZPoint(
         x * inv_tot_weight, y * inv_tot_weight, rhtools_.getPosition(cellsOnLayer.detid[maxEnergyIndex]).z());
   } else
     return math::XYZPoint(0.f, 0.f, 0.f);
 }

References f, mps_fire::i, dqm-mbProfile::log, SiStripPI::max, and findQualityFiles::v.

◆ computeThreshold()

template<typename T >

void HGCalCLUEAlgoT< T >::computeThreshold ( )

Definition at line 512 of file HGCalCLUEAlgo.cc.

                                          {
   // To support the TDR geometry and also the post-TDR one (v9 onwards), we
   // need to change the logic of the vectors containing signal to noise and
   // thresholds. The first 3 indices will keep on addressing the different
   // thicknesses of the Silicon detectors in CE_E , the next 3 indices will address
   // the thicknesses of the Silicon detectors in CE_H, while the last one, number 6 (the
   // seventh) will address the Scintillators. This change will support both
   // geometries at the same time.
  
   if (initialized_)
     return;  // only need to calculate thresholds once
  
   initialized_ = true;
  
   std::vector<double> dummy;
  
   dummy.resize(maxNumberOfThickIndices_ + !isNose_, 0);  // +1 to accomodate for the Scintillators
   thresholds_.resize(maxlayer_, dummy);
   v_sigmaNoise_.resize(maxlayer_, dummy);
  
   for (unsigned ilayer = 1; ilayer <= maxlayer_; ++ilayer) {
     for (unsigned ithick = 0; ithick < maxNumberOfThickIndices_; ++ithick) {
       float sigmaNoise = 0.001f * fcPerEle_ * nonAgedNoises_[ithick] * dEdXweights_[ilayer] /
                          (fcPerMip_[ithick] * thicknessCorrection_[ithick]);
       thresholds_[ilayer - 1][ithick] = sigmaNoise * ecut_;
       v_sigmaNoise_[ilayer - 1][ithick] = sigmaNoise;
       LogDebug("HGCalCLUEAlgo") << "ilayer: " << ilayer << " nonAgedNoises: " << nonAgedNoises_[ithick]
                                 << " fcPerEle: " << fcPerEle_ << " fcPerMip: " << fcPerMip_[ithick]
                                 << " noiseMip: " << fcPerEle_ * nonAgedNoises_[ithick] / fcPerMip_[ithick]
                                 << " sigmaNoise: " << sigmaNoise << "\n";
     }
  
     if (!isNose_) {
       float scintillators_sigmaNoise = 0.001f * noiseMip_ * dEdXweights_[ilayer] / sciThicknessCorrection_;
       thresholds_[ilayer - 1][maxNumberOfThickIndices_] = ecut_ * scintillators_sigmaNoise;
       v_sigmaNoise_[ilayer - 1][maxNumberOfThickIndices_] = scintillators_sigmaNoise;
       LogDebug("HGCalCLUEAlgo") << "ilayer: " << ilayer << " noiseMip: " << noiseMip_
                                 << " scintillators_sigmaNoise: " << scintillators_sigmaNoise << "\n";
     }
   }
 }

References LogDebug.

◆ distance()

template<typename TILE >

float HGCalCLUEAlgoT< TILE >::distance	(	int	cell1,
		int	cell2,
		int	layerId,
		bool	isEtaPhi
	)		const

inlineprivate

Definition at line 228 of file HGCalCLUEAlgo.h.

                                                                                 {  // 2-d distance on the layer (x-y)
     return std::sqrt(distance2(cell1, cell2, layerId, isEtaPhi));
   }

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

◆ distance2()

template<typename TILE >

float HGCalCLUEAlgoT< TILE >::distance2	(	int	cell1,
		int	cell2,
		int	layerId,
		bool	isEtaPhi
	)		const

inlineprivate

Definition at line 216 of file HGCalCLUEAlgo.h.

                                                                                  {  // distance squared
     if (isEtaPhi) {
       const float dphi = reco::deltaPhi(cells_[layerId].phi[cell1], cells_[layerId].phi[cell2]);
       const float deta = cells_[layerId].eta[cell1] - cells_[layerId].eta[cell2];
       return (deta * deta + dphi * dphi);
     } else {
       const float dx = cells_[layerId].x[cell1] - cells_[layerId].x[cell2];
       const float dy = cells_[layerId].y[cell1] - cells_[layerId].y[cell2];
       return (dx * dx + dy * dy);
     }
   }

References hgcalTopologyTester_cfi::cell1, hgcalTopologyTester_cfi::cell2, HGCalCLUEAlgoT< TILE >::cells_, reco::deltaPhi(), PVValHelper::dx, PVValHelper::dy, and phi.

Referenced by HGCalCLUEAlgoT< TILE >::distance().

◆ fillPSetDescription()

template<typename TILE >

static void HGCalCLUEAlgoT< TILE >::fillPSetDescription ( edm::ParameterSetDescription & iDesc )

inlinestatic

Definition at line 87 of file HGCalCLUEAlgo.h.

                                                                      {
     iDesc.add<std::vector<double>>("thresholdW0", {2.9, 2.9, 2.9});
     iDesc.add<double>("positionDeltaRho2", 1.69);
     iDesc.add<std::vector<double>>("deltac",
                                    {
                                        1.3,
                                        1.3,
                                        5.0,
                                        0.0315,  // for scintillator
                                    });
     iDesc.add<bool>("dependSensor", true);
     iDesc.add<double>("ecut", 3.0);
     iDesc.add<double>("kappa", 9.0);
     iDesc.addUntracked<unsigned int>("verbosity", 3);
     iDesc.add<std::vector<double>>("dEdXweights", {});
     iDesc.add<std::vector<double>>("thicknessCorrection", {});
     iDesc.add<double>("sciThicknessCorrection", 0.9);
     iDesc.add<int>("deltasi_index_regemfac", 3);
     iDesc.add<unsigned>("maxNumberOfThickIndices", 6);
     iDesc.add<std::vector<double>>("fcPerMip", {});
     iDesc.add<double>("fcPerEle", 0.0);
     iDesc.add<std::vector<double>>("noises", {});
     edm::ParameterSetDescription descNestedNoiseMIP;
     descNestedNoiseMIP.add<bool>("scaleByDose", false);
     descNestedNoiseMIP.add<unsigned int>("scaleByDoseAlgo", 0);
     descNestedNoiseMIP.add<double>("scaleByDoseFactor", 1.);
     descNestedNoiseMIP.add<std::string>("doseMap", "");
     descNestedNoiseMIP.add<double>("noise_MIP", 1. / 100.);
     iDesc.add<edm::ParameterSetDescription>("noiseMip", descNestedNoiseMIP);
     iDesc.add<bool>("use2x2", true);  // use 2x2 or 3x3 scenario for scint density calculation
   }

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

◆ findAndAssignClusters()

template<typename T >

int HGCalCLUEAlgoT< T >::findAndAssignClusters	(	const unsigned int	layerId,
		float	delta_c,
		float	delta_r
	)

private

Definition at line 464 of file HGCalCLUEAlgo.cc.

                                                                                                      {
   // this is called once per layer and endcap...
   // so when filling the cluster temporary vector of Hexels we resize each time
   // by the number  of clusters found. This is always equal to the number of
   // cluster centers...
   unsigned int nClustersOnLayer = 0;
   auto& cellsOnLayer = cells_[layerId];
   unsigned int numberOfCells = cellsOnLayer.detid.size();
   std::vector<int> localStack;
   // find cluster seeds and outlier
   for (unsigned int i = 0; i < numberOfCells; i++) {
     float rho_c = kappa_ * cellsOnLayer.sigmaNoise[i];
     bool isSi = rhtools_.isOnlySilicon(layerId) || cellsOnLayer.isSi[i];
     float delta = isSi ? delta_c : delta_r;
  
     // initialize clusterIndex
     cellsOnLayer.clusterIndex[i] = -1;
     bool isSeed = (cellsOnLayer.delta[i] > delta) && (cellsOnLayer.rho[i] >= rho_c);
     bool isOutlier = (cellsOnLayer.delta[i] > outlierDeltaFactor_ * delta) && (cellsOnLayer.rho[i] < rho_c);
     if (isSeed) {
       cellsOnLayer.clusterIndex[i] = nClustersOnLayer;
       cellsOnLayer.isSeed[i] = true;
       nClustersOnLayer++;
       localStack.push_back(i);
  
     } else if (!isOutlier) {
       cellsOnLayer.followers[cellsOnLayer.nearestHigher[i]].push_back(i);
     }
   }
  
   // need to pass clusterIndex to their followers
   while (!localStack.empty()) {
     int endStack = localStack.back();
     auto& thisSeed = cellsOnLayer.followers[endStack];
     localStack.pop_back();
  
     // loop over followers
     for (int j : thisSeed) {
       // pass id to a follower
       cellsOnLayer.clusterIndex[j] = cellsOnLayer.clusterIndex[endStack];
       // push this follower to localStack
       localStack.push_back(j);
     }
   }
   return nClustersOnLayer;
 }

References dumpMFGeometry_cfg::delta, hitfit::delta_r(), mps_fire::i, and dqmiolumiharvest::j.

◆ getClusters()

template<typename T >

std::vector< reco::BasicCluster > HGCalCLUEAlgoT< T >::getClusters ( bool )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 131 of file HGCalCLUEAlgo.cc.

                                                                {
   std::vector<int> offsets(numberOfClustersPerLayer_.size(), 0);
  
   int maxClustersOnLayer = numberOfClustersPerLayer_[0];
  
   for (unsigned layerId = 1; layerId < offsets.size(); ++layerId) {
     offsets[layerId] = offsets[layerId - 1] + numberOfClustersPerLayer_[layerId - 1];
  
     maxClustersOnLayer = std::max(maxClustersOnLayer, numberOfClustersPerLayer_[layerId]);
   }
  
   auto totalNumberOfClusters = offsets.back() + numberOfClustersPerLayer_.back();
   clusters_v_.resize(totalNumberOfClusters);
   std::vector<std::vector<int>> cellsIdInCluster;
   cellsIdInCluster.reserve(maxClustersOnLayer);
  
   for (unsigned int layerId = 0; layerId < 2 * maxlayer_ + 2; ++layerId) {
     cellsIdInCluster.resize(numberOfClustersPerLayer_[layerId]);
     auto& cellsOnLayer = cells_[layerId];
     unsigned int numberOfCells = cellsOnLayer.detid.size();
     auto firstClusterIdx = offsets[layerId];
  
     for (unsigned int i = 0; i < numberOfCells; ++i) {
       auto clusterIndex = cellsOnLayer.clusterIndex[i];
       if (clusterIndex != -1)
         cellsIdInCluster[clusterIndex].push_back(i);
     }
  
     std::vector<std::pair<DetId, float>> thisCluster;
  
     for (auto& cl : cellsIdInCluster) {
       auto position = calculatePosition(cl, layerId);
       float energy = 0.f;
       int seedDetId = -1;
  
       for (auto cellIdx : cl) {
         energy += cellsOnLayer.weight[cellIdx];
         thisCluster.emplace_back(cellsOnLayer.detid[cellIdx], 1.f);
         if (cellsOnLayer.isSeed[cellIdx]) {
           seedDetId = cellsOnLayer.detid[cellIdx];
         }
       }
       auto globalClusterIndex = cellsOnLayer.clusterIndex[cl[0]] + firstClusterIdx;
  
       clusters_v_[globalClusterIndex] =
           reco::BasicCluster(energy, position, reco::CaloID::DET_HGCAL_ENDCAP, thisCluster, algoId_);
       clusters_v_[globalClusterIndex].setSeed(seedDetId);
       thisCluster.clear();
     }
  
     cellsIdInCluster.clear();
   }
   return clusters_v_;
 }

References AlignmentPI::calculatePosition(), GetRecoTauVFromDQM_MC_cff::cl, reco::CaloID::DET_HGCAL_ENDCAP, HCALHighEnergyHPDFilter_cfi::energy, mps_fire::i, SiStripPI::max, unpackBuffers-CaloStage1::offsets, and position.

◆ getDensity()

template<typename T >

Density HGCalCLUEAlgoT< T >::getDensity ( )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 563 of file HGCalCLUEAlgo.cc.

                                       {
   return density_;
 }

◆ getEventSetupPerAlgorithm()

template<typename T >

void HGCalCLUEAlgoT< T >::getEventSetupPerAlgorithm ( const edm::EventSetup & es )

overridevirtual

Reimplemented from HGCalClusteringAlgoBase.

Definition at line 19 of file HGCalCLUEAlgo.cc.

                                                                          {
   cells_.clear();
   numberOfClustersPerLayer_.clear();
   cells_.resize(2 * (maxlayer_ + 1));
   numberOfClustersPerLayer_.resize(2 * (maxlayer_ + 1), 0);
 }

◆ makeClusters()

template<typename T >

void HGCalCLUEAlgoT< T >::makeClusters ( )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 99 of file HGCalCLUEAlgo.cc.

                                      {
   // assign all hits in each layer to a cluster core
   tbb::this_task_arena::isolate([&] {
     tbb::parallel_for(size_t(0), size_t(2 * maxlayer_ + 2), [&](size_t i) {
       prepareDataStructures(i);
       T lt;
       lt.clear();
       lt.fill(cells_[i].x, cells_[i].y, cells_[i].eta, cells_[i].phi, cells_[i].isSi);
       float delta_c;  // maximum search distance (critical distance) for local
                       // density calculation
       if (i % maxlayer_ < lastLayerEE_)
         delta_c = vecDeltas_[0];
       else if (i % maxlayer_ < (firstLayerBH_ - 1))
         delta_c = vecDeltas_[1];
       else
         delta_c = vecDeltas_[2];
       float delta_r = vecDeltas_[3];
       LogDebug("HGCalCLUEAlgo") << "maxlayer: " << maxlayer_ << " lastLayerEE: " << lastLayerEE_
                                 << " firstLayerBH: " << firstLayerBH_ << "\n";
  
       calculateLocalDensity(lt, i, delta_c, delta_r);
       calculateDistanceToHigher(lt, i, delta_c, delta_r);
       numberOfClustersPerLayer_[i] = findAndAssignClusters(i, delta_c, delta_r);
     });
   });
   //Now that we have the density per point we can store it
   for (unsigned int i = 0; i < 2 * maxlayer_ + 2; ++i) {
     setDensity(i);
   }
 }

References hitfit::delta_r(), PVValHelper::eta, mps_fire::i, and LogDebug.

◆ populate()

template<typename T >

void HGCalCLUEAlgoT< T >::populate ( const HGCRecHitCollection & hits )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 27 of file HGCalCLUEAlgo.cc.

                                                                 {
   // loop over all hits and create the Hexel structure, skip energies below ecut
  
   if (dependSensor_) {
     // for each layer and wafer calculate the thresholds (sigmaNoise and energy)
     // once
     computeThreshold();
   }
  
   for (unsigned int i = 0; i < hits.size(); ++i) {
     const HGCRecHit& hgrh = hits[i];
     DetId detid = hgrh.detid();
     unsigned int layerOnSide = (rhtools_.getLayerWithOffset(detid) - 1);
  
     // set sigmaNoise default value 1 to use kappa value directly in case of
     // sensor-independent thresholds
     float sigmaNoise = 1.f;
     if (dependSensor_) {
       int thickness_index = rhtools_.getSiThickIndex(detid);
       if (thickness_index == -1)
         thickness_index = maxNumberOfThickIndices_;
  
       double storedThreshold = thresholds_[layerOnSide][thickness_index];
       if (detid.det() == DetId::HGCalHSi || detid.subdetId() == HGCHEF) {
         storedThreshold = thresholds_[layerOnSide][thickness_index + deltasi_index_regemfac_];
       }
       sigmaNoise = v_sigmaNoise_[layerOnSide][thickness_index];
  
       if (hgrh.energy() < storedThreshold)
         continue;  // this sets the ZS threshold at ecut times the sigma noise
                    // for the sensor
     }
     if (!dependSensor_ && hgrh.energy() < ecut_)
       continue;
     const GlobalPoint position(rhtools_.getPosition(detid));
     int offset = ((rhtools_.zside(detid) + 1) >> 1) * maxlayer_;
     int layer = layerOnSide + offset;
  
     cells_[layer].detid.emplace_back(detid);
     cells_[layer].x.emplace_back(position.x());
     cells_[layer].y.emplace_back(position.y());
     if (!rhtools_.isOnlySilicon(layer)) {
       cells_[layer].isSi.emplace_back(rhtools_.isSilicon(detid));
       cells_[layer].eta.emplace_back(position.eta());
       cells_[layer].phi.emplace_back(position.phi());
     }  // else, isSilicon == true and eta phi values will not be used
     cells_[layer].weight.emplace_back(hgrh.energy());
     cells_[layer].sigmaNoise.emplace_back(sigmaNoise);
   }
 }

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

◆ prepareDataStructures()

template<typename T >

void HGCalCLUEAlgoT< T >::prepareDataStructures ( const unsigned int layerId )

private

Definition at line 79 of file HGCalCLUEAlgo.cc.

                                                             {
   auto cellsSize = cells_[l].detid.size();
   cells_[l].rho.resize(cellsSize, 0.f);
   cells_[l].delta.resize(cellsSize, 9999999);
   cells_[l].nearestHigher.resize(cellsSize, -1);
   cells_[l].clusterIndex.resize(cellsSize, -1);
   cells_[l].followers.resize(cellsSize);
   cells_[l].isSeed.resize(cellsSize, false);
   if (rhtools_.isOnlySilicon(l)) {
     cells_[l].isSi.resize(cellsSize, true);
     cells_[l].eta.resize(cellsSize, 0.f);
     cells_[l].phi.resize(cellsSize, 0.f);
   }
 }

References f, and cmsLHEtoEOSManager::l.

◆ reset()

template<typename TILE >

void HGCalCLUEAlgoT< TILE >::reset ( void )

inlineoverridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 69 of file HGCalCLUEAlgo.h.

                         {
     clusters_v_.clear();
     clusters_v_.shrink_to_fit();
     for (auto& cl : numberOfClustersPerLayer_) {
       cl = 0;
     }
  
     for (auto& cells : cells_) {
       cells.clear();
       cells.shrink_to_fit();
     }
     density_.clear();
   }

References cells, HGCalCLUEAlgoT< TILE >::cells_, GetRecoTauVFromDQM_MC_cff::cl, HGCalClusteringAlgoBase::clusters_v_, HGCalCLUEAlgoT< TILE >::density_, and HGCalCLUEAlgoT< TILE >::numberOfClustersPerLayer_.

◆ setDensity()

template<typename T >

void HGCalCLUEAlgoT< T >::setDensity ( const unsigned int layerId )

private

Definition at line 555 of file HGCalCLUEAlgo.cc.

                                                              {
   auto& cellsOnLayer = cells_[layerId];
   unsigned int numberOfCells = cellsOnLayer.detid.size();
   for (unsigned int i = 0; i < numberOfCells; ++i)
     density_[cellsOnLayer.detid[i]] = cellsOnLayer.rho[i];
 }

References mps_fire::i.

Member Data Documentation

◆ cells_

template<typename TILE >

std::vector<CellsOnLayer> HGCalCLUEAlgoT< TILE >::cells_

private

Definition at line 212 of file HGCalCLUEAlgo.h.

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

◆ dEdXweights_

template<typename TILE >

std::vector<double> HGCalCLUEAlgoT< TILE >::dEdXweights_

private

Definition at line 140 of file HGCalCLUEAlgo.h.

◆ deltasi_index_regemfac_

template<typename TILE >

int HGCalCLUEAlgoT< TILE >::deltasi_index_regemfac_

private

Definition at line 143 of file HGCalCLUEAlgo.h.

◆ density_

template<typename TILE >

Density HGCalCLUEAlgoT< TILE >::density_

private

Definition at line 135 of file HGCalCLUEAlgo.h.

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

◆ dependSensor_

template<typename TILE >

bool HGCalCLUEAlgoT< TILE >::dependSensor_

private

Definition at line 139 of file HGCalCLUEAlgo.h.

◆ ecut_

template<typename TILE >

double HGCalCLUEAlgoT< TILE >::ecut_

private

Definition at line 132 of file HGCalCLUEAlgo.h.

◆ fcPerEle_

template<typename TILE >

double HGCalCLUEAlgoT< TILE >::fcPerEle_

private

Definition at line 146 of file HGCalCLUEAlgo.h.

◆ fcPerMip_

template<typename TILE >

std::vector<double> HGCalCLUEAlgoT< TILE >::fcPerMip_

private

Definition at line 145 of file HGCalCLUEAlgo.h.

◆ initialized_

template<typename TILE >

bool HGCalCLUEAlgoT< TILE >::initialized_

private

Definition at line 155 of file HGCalCLUEAlgo.h.

◆ kappa_

template<typename TILE >

double HGCalCLUEAlgoT< TILE >::kappa_

private

Definition at line 129 of file HGCalCLUEAlgo.h.

◆ maxNumberOfThickIndices_

template<typename TILE >

unsigned HGCalCLUEAlgoT< TILE >::maxNumberOfThickIndices_

private

Definition at line 144 of file HGCalCLUEAlgo.h.

◆ noiseMip_

template<typename TILE >

double HGCalCLUEAlgoT< TILE >::noiseMip_

private

Definition at line 148 of file HGCalCLUEAlgo.h.

◆ nonAgedNoises_

template<typename TILE >

std::vector<double> HGCalCLUEAlgoT< TILE >::nonAgedNoises_

private

Definition at line 147 of file HGCalCLUEAlgo.h.

◆ numberOfClustersPerLayer_

template<typename TILE >

std::vector<int> HGCalCLUEAlgoT< TILE >::numberOfClustersPerLayer_

private

Definition at line 214 of file HGCalCLUEAlgo.h.

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

◆ outlierDeltaFactor_

template<typename TILE >

float HGCalCLUEAlgoT< TILE >::outlierDeltaFactor_ = 2.f

private

Definition at line 157 of file HGCalCLUEAlgo.h.

◆ positionDeltaRho2_

template<typename TILE >

const double HGCalCLUEAlgoT< TILE >::positionDeltaRho2_

private

Definition at line 125 of file HGCalCLUEAlgo.h.

◆ sciThicknessCorrection_

template<typename TILE >

double HGCalCLUEAlgoT< TILE >::sciThicknessCorrection_

private

Definition at line 142 of file HGCalCLUEAlgo.h.

◆ thicknessCorrection_

template<typename TILE >

std::vector<double> HGCalCLUEAlgoT< TILE >::thicknessCorrection_

private

Definition at line 141 of file HGCalCLUEAlgo.h.

◆ thresholds_

template<typename TILE >

std::vector<std::vector<double> > HGCalCLUEAlgoT< TILE >::thresholds_

private

Definition at line 149 of file HGCalCLUEAlgo.h.

◆ thresholdW0_

template<typename TILE >

std::vector<double> HGCalCLUEAlgoT< TILE >::thresholdW0_

private

Definition at line 124 of file HGCalCLUEAlgo.h.

◆ use2x2_

template<typename TILE >

bool HGCalCLUEAlgoT< TILE >::use2x2_

private

Definition at line 152 of file HGCalCLUEAlgo.h.

◆ v_sigmaNoise_

template<typename TILE >

std::vector<std::vector<double> > HGCalCLUEAlgoT< TILE >::v_sigmaNoise_

private

Definition at line 150 of file HGCalCLUEAlgo.h.

◆ vecDeltas_

template<typename TILE >

std::vector<double> HGCalCLUEAlgoT< TILE >::vecDeltas_

private

Definition at line 128 of file HGCalCLUEAlgo.h.

Classes

Public Types

Public Member Functions

Static Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

template<typename TILE> class HGCalCLUEAlgoT< TILE >

Member Typedef Documentation

◆ Point

Constructor & Destructor Documentation

◆ HGCalCLUEAlgoT()

◆ ~HGCalCLUEAlgoT()

Member Function Documentation

◆ calculateDistanceToHigher()

◆ calculateLocalDensity()

◆ calculatePosition()

◆ computeThreshold()

◆ distance()

◆ distance2()

◆ fillPSetDescription()

◆ findAndAssignClusters()

◆ getClusters()

◆ getDensity()

◆ getEventSetupPerAlgorithm()

◆ makeClusters()

◆ populate()

◆ prepareDataStructures()

◆ reset()

◆ setDensity()

Member Data Documentation

◆ cells_

◆ dEdXweights_

◆ deltasi_index_regemfac_

◆ density_

◆ dependSensor_

◆ ecut_

◆ fcPerEle_

◆ fcPerMip_

◆ initialized_

◆ kappa_

◆ maxNumberOfThickIndices_

◆ noiseMip_

◆ nonAgedNoises_

◆ numberOfClustersPerLayer_

◆ outlierDeltaFactor_

◆ positionDeltaRho2_

◆ sciThicknessCorrection_

◆ thicknessCorrection_

◆ thresholds_

◆ thresholdW0_

◆ use2x2_

◆ v_sigmaNoise_

◆ vecDeltas_

template<typename TILE>
class HGCalCLUEAlgoT< TILE >