#include <HGCalCLUEAlgo.h>

Inheritance diagram for HGCalCLUEAlgo:

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

	HGCalCLUEAlgo (const edm::ParameterSet &ps)

void	makeClusters () override

void	populate (const HGCRecHitCollection &hits) override

void	reset () override

	~HGCalCLUEAlgo () override

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

	HGCalClusteringAlgoBase (VerbosityLevel v, reco::CaloCluster::AlgoId algo)

void	setAlgoId (reco::CaloCluster::AlgoId algo)

void	setVerbosity (VerbosityLevel the_verbosity)

virtual	~HGCalClusteringAlgoBase ()

Static Public Member Functions
static void	fillPSetDescription (edm::ParameterSetDescription &iDesc)

Private Member Functions
void	calculateDistanceToHigher (const HGCalLayerTiles &lt, const unsigned int layerId, float delta_c, float delta_r)

void	calculateLocalDensity (const HGCalLayerTiles &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_

Density	density_

bool	dependSensor_

double	ecut_

double	fcPerEle_

std::vector< double >	fcPerMip_

bool	initialized_

double	kappa_

double	noiseMip_

std::vector< double >	nonAgedNoises_

std::vector< int >	numberOfClustersPerLayer_

float	outlierDeltaFactor_ = 2.f

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_

unsigned int	lastLayerEE_

unsigned int	lastLayerFH_

unsigned int	maxlayer_

int	scintMaxIphi_

Protected Attributes inherited from HGCalClusteringAlgoBase
reco::CaloCluster::AlgoId	algoId_

std::vector< reco::BasicCluster >	clusters_v_

hgcal::RecHitTools	rhtools_

VerbosityLevel	verbosity_

Detailed Description

Definition at line 29 of file HGCalCLUEAlgo.h.

Member Typedef Documentation

typedef math::XYZPoint HGCalCLUEAlgo::Point

point in the space

Definition at line 107 of file HGCalCLUEAlgo.h.

Constructor & Destructor Documentation

HGCalCLUEAlgo::HGCalCLUEAlgo ( const edm::ParameterSet & ps )

inline

Definition at line 31 of file HGCalCLUEAlgo.h.

       : HGCalClusteringAlgoBase(
             (HGCalClusteringAlgoBase::VerbosityLevel)ps.getUntrackedParameter<unsigned int>("verbosity", 3),
             reco::CaloCluster::undefined),
         thresholdW0_(ps.getParameter<std::vector<double>>("thresholdW0")),
         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")),
         fcPerMip_(ps.getParameter<std::vector<double>>("fcPerMip")),
         fcPerEle_(ps.getParameter<double>("fcPerEle")),
         nonAgedNoises_(ps.getParameter<edm::ParameterSet>("noises").getParameter<std::vector<double>>("values")),
         noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("noise_MIP")),
         use2x2_(ps.getParameter<bool>("use2x2")),
         initialized_(false) {}

HGCalCLUEAlgo::~HGCalCLUEAlgo ( )

inlineoverride

Definition at line 49 of file HGCalCLUEAlgo.h.

References getClusters(), getEventSetupPerAlgorithm(), hfClusterShapes_cfi::hits, makeClusters(), and populate().

49 {}

Member Function Documentation

void HGCalCLUEAlgo::calculateDistanceToHigher	(	const HGCalLayerTiles &	lt,
		const unsigned int	layerId,
		float	delta_c,
		float	delta_r
	)

private

Definition at line 334 of file HGCalCLUEAlgo.cc.

References dumpMFGeometry_cfg::delta, hitfit::delta_r(), HLT_2018_cff::distance, HGCalLayerTiles::getGlobalBinByBin(), HGCalLayerTiles::getGlobalBinByBinEtaPhi(), mps_fire::i, dqmiolumiharvest::j, LogDebug, SiStripPI::max, allConversions_cfi::maxDelta, FastTimerService_cff::range, HGCalLayerTiles::searchBox(), HGCalLayerTiles::searchBoxEtaPhi(), photonAnalyzer_cfi::xBin, and photonAnalyzer_cfi::yBin.

Referenced by distance().

                                                              {
   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";
   }
 }

void HGCalCLUEAlgo::calculateLocalDensity	(	const HGCalLayerTiles &	lt,
		const unsigned int	layerId,
		float	delta_c,
		float	delta_r
	)

private

Definition at line 233 of file HGCalCLUEAlgo.cc.

References funct::abs(), Vispa.Plugins.EdmBrowser.EdmDataAccessor::all(), dumpMFGeometry_cfg::delta, hitfit::delta_r(), HLT_2018_cff::distance, muonRecoAnalyzer_cfi::etaBin, HGCalLayerTiles::getGlobalBinByBin(), HGCalLayerTiles::getGlobalBinByBinEtaPhi(), mps_fire::i, HGCScintillatorDetId::ieta(), HGCScintillatorDetId::iphi(), dqmiolumiharvest::j, LogDebug, SiStripPI::max, BeamMonitor_cff::phiBin, HGCalLayerTiles::searchBox(), HGCalLayerTiles::searchBoxEtaPhi(), photonAnalyzer_cfi::xBin, and photonAnalyzer_cfi::yBin.

Referenced by distance().

                                                          {
   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";
   }
 }

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

private

Definition at line 177 of file HGCalCLUEAlgo.cc.

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

Referenced by distance().

                                                                                                        {
   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) {
       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);
 }

void HGCalCLUEAlgo::computeThreshold ( )

Definition at line 499 of file HGCalCLUEAlgo.cc.

References LogDebug.

Referenced by reset().

                                      {
   // 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, while the last one, number 3 (the
   // fourth) 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;
   const unsigned maxNumberOfThickIndices = 3;
   dummy.resize(maxNumberOfThickIndices + 1, 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";
     }
     float scintillators_sigmaNoise = 0.001f * noiseMip_ * dEdXweights_[ilayer];
     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";
   }
 }

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

inlineprivate

Definition at line 194 of file HGCalCLUEAlgo.h.

References calculateDistanceToHigher(), calculateLocalDensity(), calculatePosition(), hitfit::delta_r(), distance2(), findAndAssignClusters(), prepareDataStructures(), setDensity(), mathSSE::sqrt(), and findQualityFiles::v.

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

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

inlineprivate

Definition at line 182 of file HGCalCLUEAlgo.h.

References hgcalTopologyTester_cfi::cell1, hgcalTopologyTester_cfi::cell2, reco::deltaPhi(), PVValHelper::dx, PVValHelper::dy, and HGCalCLUEAlgo::CellsOnLayer::phi.

Referenced by distance().

                                                                                  {  // 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);
     }
   }

static void HGCalCLUEAlgo::fillPSetDescription ( edm::ParameterSetDescription & iDesc )

inlinestatic

Definition at line 77 of file HGCalCLUEAlgo.h.

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

                                                                      {
     iDesc.add<std::vector<double>>("thresholdW0", {2.9, 2.9, 2.9});
     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<std::vector<double>>("fcPerMip", {});
     iDesc.add<double>("fcPerEle", 0.0);
     edm::ParameterSetDescription descNestedNoises;
     descNestedNoises.add<std::vector<double>>("values", {});
     iDesc.add<edm::ParameterSetDescription>("noises", descNestedNoises);
     edm::ParameterSetDescription descNestedNoiseMIP;
     descNestedNoiseMIP.add<bool>("scaleByDose", false);
     descNestedNoiseMIP.add<unsigned int>("scaleByDoseAlgo", 0);
     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
   }

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

private

Definition at line 452 of file HGCalCLUEAlgo.cc.

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

Referenced by distance().

                                                                                                  {
   // 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;
 }

std::vector< reco::BasicCluster > HGCalCLUEAlgo::getClusters ( bool )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 122 of file HGCalCLUEAlgo.cc.

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.

Referenced by ~HGCalCLUEAlgo().

                                                            {
   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_;
 }

Density HGCalCLUEAlgo::getDensity ( )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 544 of file HGCalCLUEAlgo.cc.

Referenced by reset().

544 { return density_; }

HGCalCLUEAlgo::density_

Density density_

Definition: HGCalCLUEAlgo.h:121

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

overridevirtual

Reimplemented from HGCalClusteringAlgoBase.

Definition at line 18 of file HGCalCLUEAlgo.cc.

Referenced by ~HGCalCLUEAlgo().

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

void HGCalCLUEAlgo::makeClusters ( )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 91 of file HGCalCLUEAlgo.cc.

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

Referenced by ~HGCalCLUEAlgo().

                                  {
   // 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);
       HGCalLayerTiles 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);
   }
 }

void HGCalCLUEAlgo::populate ( const HGCRecHitCollection & hits )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 25 of file HGCalCLUEAlgo.cc.

References CaloRecHit::detid(), CaloRecHit::energy(), mps_fire::i, hltrates_dqm_sourceclient-live_cfg::offset, position, and edm::SortedCollection< T, SORT >::size().

Referenced by ~HGCalCLUEAlgo().

                                                             {
   // 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 = 3;
       double storedThreshold = thresholds_[layerOnSide][thickness_index];
       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);
   }
 }

void HGCalCLUEAlgo::prepareDataStructures ( const unsigned int layerId )

private

Definition at line 72 of file HGCalCLUEAlgo.cc.

References f, and cmsLHEtoEOSManager::l.

Referenced by distance().

                                                         {
   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);
   }
 }

void HGCalCLUEAlgo::reset ( void )

inlineoverridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 63 of file HGCalCLUEAlgo.h.

References postprocess-scan-build::cells, cells_, GetRecoTauVFromDQM_MC_cff::cl, HGCalClusteringAlgoBase::clusters_v_, computeThreshold(), getDensity(), and numberOfClustersPerLayer_.

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

void HGCalCLUEAlgo::setDensity ( const unsigned int layerId )

private

Definition at line 537 of file HGCalCLUEAlgo.cc.

References mps_fire::i.

Referenced by distance().

                                                          {
   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];
 }

Member Data Documentation

std::vector<CellsOnLayer> HGCalCLUEAlgo::cells_

private

Definition at line 178 of file HGCalCLUEAlgo.h.

Referenced by reset().

std::vector<double> HGCalCLUEAlgo::dEdXweights_

private

Definition at line 126 of file HGCalCLUEAlgo.h.

Density HGCalCLUEAlgo::density_

private

Definition at line 121 of file HGCalCLUEAlgo.h.

bool HGCalCLUEAlgo::dependSensor_

private

Definition at line 125 of file HGCalCLUEAlgo.h.

double HGCalCLUEAlgo::ecut_

private

Definition at line 118 of file HGCalCLUEAlgo.h.

double HGCalCLUEAlgo::fcPerEle_

private

Definition at line 129 of file HGCalCLUEAlgo.h.

std::vector<double> HGCalCLUEAlgo::fcPerMip_

private

Definition at line 128 of file HGCalCLUEAlgo.h.

bool HGCalCLUEAlgo::initialized_

private

Definition at line 138 of file HGCalCLUEAlgo.h.

double HGCalCLUEAlgo::kappa_

private

Definition at line 115 of file HGCalCLUEAlgo.h.

double HGCalCLUEAlgo::noiseMip_

private

Definition at line 131 of file HGCalCLUEAlgo.h.

std::vector<double> HGCalCLUEAlgo::nonAgedNoises_

private

Definition at line 130 of file HGCalCLUEAlgo.h.

std::vector<int> HGCalCLUEAlgo::numberOfClustersPerLayer_

private

Definition at line 180 of file HGCalCLUEAlgo.h.

Referenced by reset().

float HGCalCLUEAlgo::outlierDeltaFactor_ = 2.f

private

Definition at line 140 of file HGCalCLUEAlgo.h.

std::vector<double> HGCalCLUEAlgo::thicknessCorrection_

private

Definition at line 127 of file HGCalCLUEAlgo.h.

std::vector<std::vector<double> > HGCalCLUEAlgo::thresholds_

private

Definition at line 132 of file HGCalCLUEAlgo.h.

std::vector<double> HGCalCLUEAlgo::thresholdW0_

private

Definition at line 111 of file HGCalCLUEAlgo.h.

bool HGCalCLUEAlgo::use2x2_

private

Definition at line 135 of file HGCalCLUEAlgo.h.

std::vector<std::vector<double> > HGCalCLUEAlgo::v_sigmaNoise_

private

Definition at line 133 of file HGCalCLUEAlgo.h.

std::vector<double> HGCalCLUEAlgo::vecDeltas_

private

Definition at line 114 of file HGCalCLUEAlgo.h.

Classes

Public Types

Public Member Functions

Static Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation