HGCalImagingAlgo Class Reference

#include <HGCalImagingAlgo.h>

Inheritance diagram for HGCalImagingAlgo:

Classes
struct	Hexel

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
	HGCalImagingAlgo (const edm::ParameterSet &ps, edm::ConsumesCollector iC)
void	makeClusters () override
void	populate (const HGCRecHitCollection &hits) override
void	reset () override
	~HGCalImagingAlgo () 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 Types
typedef KDTreeNodeInfo< Hexel >	KDNode
typedef KDTreeLinkerAlgo< Hexel >	KDTree

Private Member Functions
double	calculateDistanceToHigher (std::vector< KDNode > &) const
double	calculateEnergyWithFraction (const std::vector< KDNode > &, const std::vector< double > &)
double	calculateLocalDensity (std::vector< KDNode > &, KDTree &, const unsigned int) const
math::XYZPoint	calculatePosition (std::vector< KDNode > &) const
math::XYZPoint	calculatePositionWithFraction (const std::vector< KDNode > &, const std::vector< double > &)
double	distance (const Hexel &pt1, const Hexel &pt2) const
double	distance2 (const Hexel &pt1, const Hexel &pt2) const
int	findAndAssignClusters (std::vector< KDNode > &, KDTree &, double, KDTreeBox< 2 > &, const unsigned int, std::vector< std::vector< KDNode >> &) const
std::vector< unsigned >	findLocalMaximaInCluster (const std::vector< KDNode > &)
void	setDensity (const std::vector< KDNode > &nd)
void	shareEnergy (const std::vector< KDNode > &, const std::vector< unsigned > &, std::vector< std::vector< double >> &)
std::vector< size_t >	sort_by_delta (const std::vector< KDNode > &v) const

Private Attributes
std::vector< reco::BasicCluster >	clusters_v_
std::vector< double >	dEdXweights_
Density	density_
bool	dependSensor_
double	ecut_
double	fcPerEle_
std::vector< double >	fcPerMip_
bool	initialized_
double	kappa_
std::vector< std::vector< std::vector< KDNode > > >	layerClustersPerLayer_
std::vector< std::array< float, 2 > >	maxpos_
std::vector< std::array< float, 2 > >	minpos_
double	noiseMip_
std::vector< double >	nonAgedNoises_
std::vector< std::vector< KDNode > >	points_
std::vector< double >	positionDeltaRho_c_
double	sigma2_
std::vector< std::vector< double > >	sigmaNoise_
std::vector< double >	thicknessCorrection_
std::vector< std::vector< double > >	thresholds_
std::vector< double >	thresholdW0_
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

Definition at line 31 of file HGCalImagingAlgo.h.

Member Typedef Documentation

KDNode

typedef KDTreeNodeInfo<Hexel> HGCalImagingAlgo::KDNode

private

Definition at line 221 of file HGCalImagingAlgo.h.

KDTree

typedef KDTreeLinkerAlgo<Hexel> HGCalImagingAlgo::KDTree

private

Definition at line 220 of file HGCalImagingAlgo.h.

Point

typedef math::XYZPoint HGCalImagingAlgo::Point

point in the space

Definition at line 120 of file HGCalImagingAlgo.h.

Constructor & Destructor Documentation

HGCalImagingAlgo()

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

inline

Definition at line 33 of file HGCalImagingAlgo.h.

      : HGCalClusteringAlgoBase(
            (HGCalClusteringAlgoBase::VerbosityLevel)ps.getUntrackedParameter<unsigned int>("verbosity", 3),
            reco::CaloCluster::undefined,
            iC),
        thresholdW0_(ps.getParameter<std::vector<double>>("thresholdW0")),
        positionDeltaRho_c_(ps.getParameter<std::vector<double>>("positionDeltaRho_c")),
        vecDeltas_(ps.getParameter<std::vector<double>>("deltac")),
        kappa_(ps.getParameter<double>("kappa")),
        ecut_(ps.getParameter<double>("ecut")),
        sigma2_(1.0),
        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")),
        initialized_(false) {}

~HGCalImagingAlgo()

HGCalImagingAlgo::~HGCalImagingAlgo ( )

inlineoverride

Definition at line 53 of file HGCalImagingAlgo.h.

{}

Member Function Documentation

calculateDistanceToHigher()

double HGCalImagingAlgo::calculateDistanceToHigher ( std::vector< KDNode > &  nd ) const

private

Definition at line 291 of file HGCalImagingAlgo.cc.

References data, cms::cuda::for(), mps_fire::i, dqmiolumiharvest::j, hgcal_clustering::sorted_indices(), mathSSE::sqrt(), and createJobs::tmp.

                                                                              {
  // sort vector of Hexels by decreasing local density
  std::vector<size_t> rs = sorted_indices(nd);

  double maxdensity = 0.0;
  int nearestHigher = -1;

  if (!rs.empty())
    maxdensity = nd[rs[0]].data.rho;
  else
    return maxdensity;  // there are no hits
  double dist2 = 0.;
  // start by setting delta for the highest density hit to
  // the most distant hit - this is a convention

  for (auto &j : nd) {
    double tmp = distance2(nd[rs[0]].data, j.data);
    if (tmp > dist2)
      dist2 = tmp;
  }
  nd[rs[0]].data.delta = std::sqrt(dist2);
  nd[rs[0]].data.nearestHigher = nearestHigher;

  // now we save the largest distance as a starting point
  const double max_dist2 = dist2;
  const unsigned int nd_size = nd.size();

  for (unsigned int oi = 1; oi < nd_size; ++oi) {  // start from second-highest density
    dist2 = max_dist2;
    unsigned int i = rs[oi];
    // we only need to check up to oi since hits
    // are ordered by decreasing density
    // and all points coming BEFORE oi are guaranteed to have higher rho
    // and the ones AFTER to have lower rho
    for (unsigned int oj = 0; oj < oi; ++oj) {
      unsigned int j = rs[oj];
      // Limit the search box
      if ((nd[i].data.x - nd[j].data.x) * (nd[i].data.x - nd[j].data.x) > dist2)
        continue;
      if ((nd[i].data.y - nd[j].data.y) * (nd[i].data.y - nd[j].data.y) > dist2)
        continue;
      double tmp = distance2(nd[i].data, nd[j].data);
      if (tmp <= dist2) {  // this "<=" instead of "<" addresses the (rare) case
                           // when there are only two hits
        dist2 = tmp;
        nearestHigher = j;
      }
    }
    nd[i].data.delta = std::sqrt(dist2);
    nd[i].data.nearestHigher = nearestHigher;  // this uses the original unsorted hitlist
  }
  return maxdensity;
}

calculateEnergyWithFraction()

double HGCalImagingAlgo::calculateEnergyWithFraction ( const std::vector< KDNode > &  hits, const std::vector< double > &  fractions  )

private

Definition at line 522 of file HGCalImagingAlgo.cc.

References hfClusterShapes_cfi::hits, mps_fire::i, and mps_fire::result.

                                                                                         {
  double result = 0.0;
  for (unsigned i = 0; i < hits.size(); ++i) {
    result += fractions[i] * hits[i].data.weight;
  }
  return result;
}

calculateLocalDensity()

double HGCalImagingAlgo::calculateLocalDensity ( std::vector< KDNode > &  nd, KDTree &  lp, const unsigned int  layer  ) const

private

Definition at line 262 of file HGCalImagingAlgo.cc.

References data, HLT_2022v15_cff::distance, newFWLiteAna::found, mps_fire::i, dqmiolumiharvest::j, SiStripPI::max, and KDTreeLinkerAlgo< DATA, DIM >::search().

                                                                                                                {
  double maxdensity = 0.;
  float delta_c;  // maximum search distance (critical distance) for local
                  // density calculation
  if (layer <= lastLayerEE_)
    delta_c = vecDeltas_[0];
  else if (layer < firstLayerBH_)
    delta_c = vecDeltas_[1];
  else
    delta_c = vecDeltas_[2];

  // for each node calculate local density rho and store it
  for (unsigned int i = 0; i < nd.size(); ++i) {
    // speec up search by looking within +/- delta_c window only
    KDTreeBox search_box(
        nd[i].dims[0] - delta_c, nd[i].dims[0] + delta_c, nd[i].dims[1] - delta_c, nd[i].dims[1] + delta_c);
    std::vector<Hexel> found;
    lp.search(search_box, found);
    const unsigned int found_size = found.size();
    for (unsigned int j = 0; j < found_size; j++) {
      if (distance(nd[i].data, found[j]) < delta_c) {
        nd[i].data.rho += found[j].weight;
        maxdensity = std::max(maxdensity, nd[i].data.rho);
      }
    }  // end loop found
  }    // end loop nodes
  return maxdensity;
}

calculatePosition()

math::XYZPoint HGCalImagingAlgo::calculatePosition ( std::vector< KDNode > &  v ) const

private

Definition at line 196 of file HGCalImagingAlgo.cc.

References data, mps_fire::i, heavyIonCSV_trainingSettings::idx, dqm-mbProfile::log, SiStripPI::max, findQualityFiles::v, and x.

                                                                           {
  float total_weight = 0.f;
  float x = 0.f;
  float y = 0.f;

  unsigned int v_size = v.size();
  unsigned int maxEnergyIndex = 0;
  float maxEnergyValue = 0;

  // loop over hits in cluster candidate
  // determining the maximum energy hit
  for (unsigned int i = 0; i < v_size; i++) {
    if (v[i].data.weight > maxEnergyValue) {
      maxEnergyValue = v[i].data.weight;
      maxEnergyIndex = i;
    }
  }

  // Si cell or Scintillator. Used to set approach and parameters
  int thick = rhtools_.getSiThickIndex(v[maxEnergyIndex].data.detid);

  // for hits within positionDeltaRho_c_ from maximum energy hit
  // build up weight for energy-weighted position
  // and save corresponding hits indices
  std::vector<unsigned int> innerIndices;
  for (unsigned int i = 0; i < v_size; i++) {
    if (thick == -1 || distance2(v[i].data, v[maxEnergyIndex].data) < positionDeltaRho_c_[thick]) {
      innerIndices.push_back(i);

      float rhEnergy = v[i].data.weight;
      total_weight += rhEnergy;
      // just fill x, y for scintillator
      // for Si it is overwritten later anyway
      if (thick == -1) {
        x += v[i].data.x * rhEnergy;
        y += v[i].data.y * rhEnergy;
      }
    }
  }
  // just loop on reduced vector of interesting indices
  // to compute log weighting
  if (thick != -1 && total_weight != 0.) {  // Silicon case
    float total_weight_log = 0.f;
    float x_log = 0.f;
    float y_log = 0.f;
    for (auto idx : innerIndices) {
      float rhEnergy = v[idx].data.weight;
      if (rhEnergy == 0.)
        continue;
      float Wi = std::max(thresholdW0_[thick] + log(rhEnergy / total_weight), 0.);
      x_log += v[idx].data.x * Wi;
      y_log += v[idx].data.y * Wi;
      total_weight_log += Wi;
    }
    total_weight = total_weight_log;
    x = x_log;
    y = y_log;
  }

  if (total_weight != 0.) {
    auto inv_tot_weight = 1. / total_weight;
    return math::XYZPoint(x * inv_tot_weight, y * inv_tot_weight, v[maxEnergyIndex].data.z);
  }
  return math::XYZPoint(0, 0, 0);
}

calculatePositionWithFraction()

math::XYZPoint HGCalImagingAlgo::calculatePositionWithFraction ( const std::vector< KDNode > &  hits, const std::vector< double > &  fractions  )

private

Definition at line 507 of file HGCalImagingAlgo.cc.

References hfClusterShapes_cfi::hits, mps_fire::i, mps_fire::result, mps_merge::weight, and x.

                                                                                                   {
  double norm(0.0), x(0.0), y(0.0), z(0.0);
  for (unsigned i = 0; i < hits.size(); ++i) {
    const double weight = fractions[i] * hits[i].data.weight;
    norm += weight;
    x += weight * hits[i].data.x;
    y += weight * hits[i].data.y;
    z += weight * hits[i].data.z;
  }
  math::XYZPoint result(x, y, z);
  result /= norm;
  return result;
}

computeThreshold()

void HGCalImagingAlgo::computeThreshold

(

)

Definition at line 628 of file HGCalImagingAlgo.cc.

References hgcalLayerClusters_cfi::maxNumberOfThickIndices.

                                        {
  // 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);
  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_;
      sigmaNoise_[ilayer - 1][ithick] = sigmaNoise;
    }
    float scintillators_sigmaNoise = 0.001f * noiseMip_ * dEdXweights_[ilayer];
    thresholds_[ilayer - 1][maxNumberOfThickIndices] = ecut_ * scintillators_sigmaNoise;
    sigmaNoise_[ilayer - 1][maxNumberOfThickIndices] = scintillators_sigmaNoise;
  }
}

distance()

double HGCalImagingAlgo::distance ( const Hexel &  pt1, const Hexel &  pt2  ) const

inlineprivate

Definition at line 244 of file HGCalImagingAlgo.h.

References distance2(), HLT_2022v15_cff::pt1, HLT_2022v15_cff::pt2, and mathSSE::sqrt().

                                                                   {  //2-d distance on the layer (x-y)
    return std::sqrt(distance2(pt1, pt2));
  }

distance2()

double HGCalImagingAlgo::distance2 ( const Hexel &  pt1, const Hexel &  pt2  ) const

inlineprivate

Definition at line 239 of file HGCalImagingAlgo.h.

References PVValHelper::dx, PVValHelper::dy, HLT_2022v15_cff::pt1, and HLT_2022v15_cff::pt2.

Referenced by distance().

                                                                    {  //distance squared
    const double dx = pt1.x - pt2.x;
    const double dy = pt1.y - pt2.y;
    return (dx * dx + dy * dy);
  }                                                                   //distance squaredq

fillPSetDescription()

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

inlinestatic

Definition at line 89 of file HGCalImagingAlgo.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>>("positionDeltaRho_c", {1.3, 1.3, 1.3});
    iDesc.add<std::vector<double>>("deltac",
                                   {
                                       2.0,
                                       2.0,
                                       5.0,
                                   });
    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<double>("scaleByDoseFactor", 1.);
    iDesc.add<edm::ParameterSetDescription>("scaleByDose", descNestedNoiseMIP);
    descNestedNoiseMIP.add<std::string>("doseMap", "");
    iDesc.add<edm::ParameterSetDescription>("doseMap", descNestedNoiseMIP);
    descNestedNoiseMIP.add<double>("noise_MIP", 1. / 100.);
    iDesc.add<edm::ParameterSetDescription>("noiseMip", descNestedNoiseMIP);
  }

findAndAssignClusters()

int HGCalImagingAlgo::findAndAssignClusters ( std::vector< KDNode > &  nd, KDTree &  lp, double  maxdensity, KDTreeBox< 2 > &  bounds, const unsigned int  layer, std::vector< std::vector< KDNode >> &  clustersOnLayer  ) const

private

Definition at line 344 of file HGCalImagingAlgo.cc.

References KDTreeLinkerAlgo< DATA, DIM >::build(), KDTreeLinkerAlgo< DATA, DIM >::clear(), gather_cfg::cout, data, HLT_2022v15_cff::distance, newFWLiteAna::found, mps_fire::i, dqmiolumiharvest::j, KDTreeLinkerAlgo< DATA, DIM >::search(), and hgcal_clustering::sorted_indices().

                                                                                                 {
  // 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;
  float delta_c;  // critical distance
  if (layer <= lastLayerEE_)
    delta_c = vecDeltas_[0];
  else if (layer < firstLayerBH_)
    delta_c = vecDeltas_[1];
  else
    delta_c = vecDeltas_[2];

  std::vector<size_t> rs = sorted_indices(nd);  // indices sorted by decreasing rho
  std::vector<size_t> ds = sort_by_delta(nd);   // sort in decreasing distance to higher

  const unsigned int nd_size = nd.size();
  for (unsigned int i = 0; i < nd_size; ++i) {
    if (nd[ds[i]].data.delta < delta_c)
      break;  // no more cluster centers to be looked at
    if (dependSensor_) {
      float rho_c = kappa_ * nd[ds[i]].data.sigmaNoise;
      if (nd[ds[i]].data.rho < rho_c)
        continue;  // set equal to kappa times noise threshold

    } else if (nd[ds[i]].data.rho * kappa_ < maxdensity)
      continue;

    nd[ds[i]].data.clusterIndex = nClustersOnLayer;
    if (verbosity_ < pINFO) {
      std::cout << "Adding new cluster with index " << nClustersOnLayer << std::endl;
      std::cout << "Cluster center is hit " << ds[i] << std::endl;
    }
    nClustersOnLayer++;
  }

  // at this point nClustersOnLayer is equal to the number of cluster centers -
  // if it is zero we are  done
  if (nClustersOnLayer == 0)
    return nClustersOnLayer;

  // assign remaining points to clusters, using the nearestHigher set from
  // previous step (always set except
  // for top density hit that is skipped...)
  for (unsigned int oi = 1; oi < nd_size; ++oi) {
    unsigned int i = rs[oi];
    int ci = nd[i].data.clusterIndex;
    if (ci == -1) {  // clusterIndex is initialised with -1 if not yet used in cluster
      nd[i].data.clusterIndex = nd[nd[i].data.nearestHigher].data.clusterIndex;
    }
  }

  // make room in the temporary cluster vector for the additional clusterIndex
  // clusters
  // from this layer
  if (verbosity_ < pINFO) {
    std::cout << "resizing cluster vector by " << nClustersOnLayer << std::endl;
  }
  clustersOnLayer.resize(nClustersOnLayer);

  // assign points closer than dc to other clusters to border region
  // and find critical border density
  std::vector<double> rho_b(nClustersOnLayer, 0.);
  lp.clear();
  lp.build(nd, bounds);
  // now loop on all hits again :( and check: if there are hits from another
  // cluster within d_c -> flag as border hit
  for (unsigned int i = 0; i < nd_size; ++i) {
    int ci = nd[i].data.clusterIndex;
    bool flag_isolated = true;
    if (ci != -1) {
      KDTreeBox search_box(
          nd[i].dims[0] - delta_c, nd[i].dims[0] + delta_c, nd[i].dims[1] - delta_c, nd[i].dims[1] + delta_c);
      std::vector<Hexel> found;
      lp.search(search_box, found);

      const unsigned int found_size = found.size();
      for (unsigned int j = 0; j < found_size; j++) {  // start from 0 here instead of 1
        // check if the hit is not within d_c of another cluster
        if (found[j].clusterIndex != -1) {
          float dist = distance(found[j], nd[i].data);
          if (dist < delta_c && found[j].clusterIndex != ci) {
            // in which case we assign it to the border
            nd[i].data.isBorder = true;
            break;
          }
          // because we are using two different containers, we have to make sure
          // that we don't unflag the
          // hit when it finds *itself* closer than delta_c
          if (dist < delta_c && dist != 0. && found[j].clusterIndex == ci) {
            // in this case it is not an isolated hit
            // the dist!=0 is because the hit being looked at is also inside the
            // search box and at dist==0
            flag_isolated = false;
          }
        }
      }
      if (flag_isolated)
        nd[i].data.isBorder = true;  // the hit is more than delta_c from any of its brethren
    }
    // check if this border hit has density larger than the current rho_b and
    // update
    if (nd[i].data.isBorder && rho_b[ci] < nd[i].data.rho)
      rho_b[ci] = nd[i].data.rho;
  }  // end loop all hits

  // flag points in cluster with density < rho_b as halo points, then fill the
  // cluster vector
  for (unsigned int i = 0; i < nd_size; ++i) {
    int ci = nd[i].data.clusterIndex;
    if (ci != -1) {
      if (nd[i].data.rho <= rho_b[ci])
        nd[i].data.isHalo = true;
      clustersOnLayer[ci].push_back(nd[i]);
      if (verbosity_ < pINFO) {
        std::cout << "Pushing hit " << i << " into cluster with index " << ci << std::endl;
      }
    }
  }

  // prepare the offset for the next layer if there is one
  if (verbosity_ < pINFO) {
    std::cout << "moving cluster offset by " << nClustersOnLayer << std::endl;
  }
  return nClustersOnLayer;
}

findLocalMaximaInCluster()

std::vector< unsigned > HGCalImagingAlgo::findLocalMaximaInCluster ( const std::vector< KDNode > &  cluster )

private

Definition at line 479 of file HGCalImagingAlgo.cc.

References data, HLT_2022v15_cff::distance, MillePedeFileConverter_cfg::e, mps_fire::i, dqmiolumiharvest::j, mps_fire::result, and fileCollector::seed.

                                                                                               {
  std::vector<unsigned> result;
  std::vector<bool> seed(cluster.size(), true);
  float delta_c = 2.;

  for (unsigned i = 0; i < cluster.size(); ++i) {
    for (unsigned j = 0; j < cluster.size(); ++j) {
      if (i != j and distance(cluster[i].data, cluster[j].data) < delta_c) {
        if (cluster[i].data.weight < cluster[j].data.weight) {
          seed[i] = false;
          break;
        }
      }
    }
  }

  for (unsigned i = 0; i < cluster.size(); ++i) {
    if (seed[i] && cluster[i].data.weight > 5e-4) {
      // seed at i with energy cluster[i].weight
      result.push_back(i);
    }
  }

  // Found result.size() sub-clusters in input cluster of length cluster.size()

  return result;
}

getClusters()

std::vector< reco::BasicCluster > HGCalImagingAlgo::getClusters ( bool  doSharing )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 119 of file HGCalImagingAlgo.cc.

References AlignmentPI::calculatePosition(), gather_cfg::cout, data, reco::CaloID::DET_HGCAL_ENDCAP, hgcalLayerClusters_cfi::doSharing, MillePedeFileConverter_cfg::e, HCALHighEnergyHPDFilter_cfi::energy, HLT_2022v15_cff::fraction, mps_fire::i, hgcal_clustering::max_index(), position, and DetachedQuadStep_cff::seeds.

                                                                        {
  reco::CaloID caloID = reco::CaloID::DET_HGCAL_ENDCAP;
  std::vector<std::pair<DetId, float>> thisCluster;
  for (auto &clsOnLayer : layerClustersPerLayer_) {
    for (unsigned int i = 0; i < clsOnLayer.size(); ++i) {
      double energy = 0;
      Point position;

      //Will save the maximum density hit of the cluster
      size_t rsmax = max_index(clsOnLayer[i]);

      if (doSharing) {
        std::vector<unsigned> seeds = findLocalMaximaInCluster(clsOnLayer[i]);
        // sharing found seeds.size() sub-cluster seeds in cluster i

        std::vector<std::vector<double>> fractions;
        // first pass can have noise it in
        shareEnergy(clsOnLayer[i], seeds, fractions);

        // reset and run second pass after vetoing seeds
        // that result in trivial clusters (less than 2 effective cells)

        for (unsigned isub = 0; isub < fractions.size(); ++isub) {
          double effective_hits = 0.0;
          double energy = calculateEnergyWithFraction(clsOnLayer[i], fractions[isub]);
          Point position = calculatePositionWithFraction(clsOnLayer[i], fractions[isub]);

          for (unsigned ihit = 0; ihit < fractions[isub].size(); ++ihit) {
            const double fraction = fractions[isub][ihit];
            if (fraction > 1e-7) {
              effective_hits += fraction;
              thisCluster.emplace_back(clsOnLayer[i][ihit].data.detid, fraction);
            }
          }

          if (verbosity_ < pINFO) {
            std::cout << "\t******** NEW CLUSTER (SHARING) ********" << std::endl;
            std::cout << "\tEff. No. of cells = " << effective_hits << std::endl;
            std::cout << "\t     Energy       = " << energy << std::endl;
            std::cout << "\t     Phi          = " << position.phi() << std::endl;
            std::cout << "\t     Eta          = " << position.eta() << std::endl;
            std::cout << "\t*****************************" << std::endl;
          }
          clusters_v_.emplace_back(energy, position, caloID, thisCluster, algoId_);
          if (!clusters_v_.empty()) {
            clusters_v_.back().setSeed(clsOnLayer[i][rsmax].data.detid);
          }
          thisCluster.clear();
        }
      } else {
        position = calculatePosition(clsOnLayer[i]);  // energy-weighted position
        //   std::vector< KDNode >::iterator it;
        for (auto &it : clsOnLayer[i]) {
          energy += it.data.isHalo ? 0. : it.data.weight;
          // use fraction to store whether this is a Halo hit or not
          thisCluster.emplace_back(it.data.detid, (it.data.isHalo ? 0.f : 1.f));
        }
        if (verbosity_ < pINFO) {
          std::cout << "******** NEW CLUSTER (HGCIA) ********" << std::endl;
          std::cout << "Index          " << i << std::endl;
          std::cout << "No. of cells = " << clsOnLayer[i].size() << std::endl;
          std::cout << "     Energy     = " << energy << std::endl;
          std::cout << "     Phi        = " << position.phi() << std::endl;
          std::cout << "     Eta        = " << position.eta() << std::endl;
          std::cout << "*****************************" << std::endl;
        }
        clusters_v_.emplace_back(energy, position, caloID, thisCluster, algoId_);
        if (!clusters_v_.empty()) {
          clusters_v_.back().setSeed(clsOnLayer[i][rsmax].data.detid);
        }
        thisCluster.clear();
      }
    }
  }
  return clusters_v_;
}

getDensity()

Density HGCalImagingAlgo::getDensity ( )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 668 of file HGCalImagingAlgo.cc.

{ return density_; }

getEventSetupPerAlgorithm()

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

overridevirtual

Reimplemented from HGCalClusteringAlgoBase.

Definition at line 17 of file HGCalImagingAlgo.cc.

                                                                        {
  points_.clear();
  minpos_.clear();
  maxpos_.clear();
  points_.resize(2 * (maxlayer_ + 1));
  minpos_.resize(2 * (maxlayer_ + 1), {{0.0f, 0.0f}});
  maxpos_.resize(2 * (maxlayer_ + 1), {{0.0f, 0.0f}});
}

makeClusters()

void HGCalImagingAlgo::makeClusters ( )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 93 of file HGCalImagingAlgo.cc.

References KDTreeLinkerAlgo< DATA, DIM >::build(), mps_fire::i, and AlCaHLTBitMon_ParallelJobs::p.

                                    {
  layerClustersPerLayer_.resize(2 * maxlayer_ + 2);
  // assign all hits in each layer to a cluster core or halo
  tbb::this_task_arena::isolate([&] {
    tbb::parallel_for(size_t(0), size_t(2 * maxlayer_ + 2), [&](size_t i) {
      KDTreeBox bounds(minpos_[i][0], maxpos_[i][0], minpos_[i][1], maxpos_[i][1]);
      KDTree hit_kdtree;
      hit_kdtree.build(points_[i], bounds);

      unsigned int actualLayer =
          i > maxlayer_ ? (i - (maxlayer_ + 1)) : i;  // maps back from index used for KD trees to actual layer

      double maxdensity = calculateLocalDensity(points_[i], hit_kdtree, actualLayer);  // also stores rho (energy
                                                                                       // density) for each point (node)
      // calculate distance to nearest point with higher density storing
      // distance (delta) and point's index
      calculateDistanceToHigher(points_[i]);
      findAndAssignClusters(points_[i], hit_kdtree, maxdensity, bounds, actualLayer, layerClustersPerLayer_[i]);
    });
  });
  //Now that we have the density per point we can store it
  for (auto const &p : points_) {
    setDensity(p);
  }
}

populate()

void HGCalImagingAlgo::populate ( const HGCRecHitCollection &  hits )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 26 of file HGCalImagingAlgo.cc.

References CaloRecHit::detid(), CaloRecHit::energy(), hfClusterShapes_cfi::hits, mps_fire::i, createfilelist::int, SiStripPI::max, SiStripPI::min, position, and Calorimetry_cff::thickness.

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

  std::vector<bool> firstHit(2 * (maxlayer_ + 1), true);
  for (unsigned int i = 0; i < hits.size(); ++i) {
    const HGCRecHit &hgrh = hits[i];
    DetId detid = hgrh.detid();
    unsigned int layer = rhtools_.getLayerWithOffset(detid);
    float thickness = 0.f;
    // set sigmaNoise default value 1 to use kappa value directly in case of
    // sensor-independent thresholds
    float sigmaNoise = 1.f;
    if (dependSensor_) {
      thickness = rhtools_.getSiThickness(detid);
      int thickness_index = rhtools_.getSiThickIndex(detid);
      if (thickness_index == -1)
        thickness_index = 3;
      double storedThreshold = thresholds_[layer - 1][thickness_index];
      sigmaNoise = sigmaNoise_[layer - 1][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;

    // map layers from positive endcap (z) to layer + maxlayer+1 to prevent
    // mixing up hits from different sides
    layer += int(rhtools_.zside(detid) > 0) * (maxlayer_ + 1);

    // determine whether this is a half-hexagon
    bool isHalf = rhtools_.isHalfCell(detid);
    const GlobalPoint position(rhtools_.getPosition(detid));

    // here's were the KDNode is passed its dims arguments - note that these are
    // *copied* from the Hexel
    points_[layer].emplace_back(
        Hexel(hgrh, detid, isHalf, sigmaNoise, thickness, &rhtools_), position.x(), position.y());

    // for each layer, store the minimum and maximum x and y coordinates for the
    // KDTreeBox boundaries
    if (firstHit[layer]) {
      minpos_[layer][0] = position.x();
      minpos_[layer][1] = position.y();
      maxpos_[layer][0] = position.x();
      maxpos_[layer][1] = position.y();
      firstHit[layer] = false;
    } else {
      minpos_[layer][0] = std::min((float)position.x(), minpos_[layer][0]);
      minpos_[layer][1] = std::min((float)position.y(), minpos_[layer][1]);
      maxpos_[layer][0] = std::max((float)position.x(), maxpos_[layer][0]);
      maxpos_[layer][1] = std::max((float)position.y(), maxpos_[layer][1]);
    }

  }  // end loop hits
}

reset()

void HGCalImagingAlgo::reset ( void  )

inlineoverridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 67 of file HGCalImagingAlgo.h.

References clusters_v_, mps_fire::i, layerClustersPerLayer_, maxpos_, minpos_, points_, and edm::swap().

                        {
    clusters_v_.clear();
    clusters_v_.shrink_to_fit();
    layerClustersPerLayer_.clear();
    layerClustersPerLayer_.shrink_to_fit();
    for (auto &it : points_) {
      it.clear();
      it.shrink_to_fit();
      std::vector<KDNode>().swap(it);
    }
    for (unsigned int i = 0; i < minpos_.size(); i++) {
      minpos_[i][0] = 0.;
      minpos_[i][1] = 0.;
      maxpos_[i][0] = 0.;
      maxpos_[i][1] = 0.;
    }
  }

setDensity()

void HGCalImagingAlgo::setDensity ( const std::vector< KDNode > &  nd )

private

Definition at line 660 of file HGCalImagingAlgo.cc.

References mps_fire::i.

                                                             {
  // for each node calculate local density rho and store it
  for (auto &i : nd) {
    density_[i.data.detid] = i.data.rho;
  }  // end loop nodes
}

shareEnergy()

void HGCalImagingAlgo::shareEnergy ( const std::vector< KDNode > &  incluster, const std::vector< unsigned > &  seeds, std::vector< std::vector< double >> &  outclusters  )

private

Definition at line 531 of file HGCalImagingAlgo.cc.

References data, change_name::diff, MillePedeFileConverter_cfg::e, JetChargeProducer_cfi::exp, DivergingColor::frac, HLT_2022v15_cff::fraction, mps_fire::i, dqmiolumiharvest::j, SiStripPI::max, HLT_2022v15_cff::minFracTot, eostools::move(), perp2(), funct::pow(), DetachedQuadStep_cff::seeds, mathSSE::sqrt(), HLT_2022v15_cff::stoppingTolerance, HGCalImagingAlgo::Hexel::x, HGCalImagingAlgo::Hexel::y, and HGCalImagingAlgo::Hexel::z.

                                                                              {
  std::vector<bool> isaseed(incluster.size(), false);
  outclusters.clear();
  outclusters.resize(seeds.size());
  std::vector<Point> centroids(seeds.size());
  std::vector<double> energies(seeds.size());

  if (seeds.size() == 1) {  // short circuit the case of a lone cluster
    outclusters[0].clear();
    outclusters[0].resize(incluster.size(), 1.0);
    return;
  }

  // saving seeds

  // create quick seed lookup
  for (unsigned i = 0; i < seeds.size(); ++i) {
    isaseed[seeds[i]] = true;
  }

  // initialize clusters to be shared
  // centroids start off at seed positions
  // seeds always have fraction 1.0, to stabilize fit
  // initializing fit
  for (unsigned i = 0; i < seeds.size(); ++i) {
    outclusters[i].resize(incluster.size(), 0.0);
    for (unsigned j = 0; j < incluster.size(); ++j) {
      if (j == seeds[i]) {
        outclusters[i][j] = 1.0;
        centroids[i] = math::XYZPoint(incluster[j].data.x, incluster[j].data.y, incluster[j].data.z);
        energies[i] = incluster[j].data.weight;
      }
    }
  }

  // run the fit while we are less than max iterations, and clusters are still
  // moving
  const double minFracTot = 1e-20;
  unsigned iter = 0;
  const unsigned iterMax = 50;
  double diff = std::numeric_limits<double>::max();
  const double stoppingTolerance = 1e-8;
  const auto numberOfSeeds = seeds.size();
  auto toleranceScaling = numberOfSeeds > 2 ? (numberOfSeeds - 1) * (numberOfSeeds - 1) : 1;
  std::vector<Point> prevCentroids;
  std::vector<double> frac(numberOfSeeds), dist2(numberOfSeeds);
  while (iter++ < iterMax && diff > stoppingTolerance * toleranceScaling) {
    for (unsigned i = 0; i < incluster.size(); ++i) {
      const Hexel &ihit = incluster[i].data;
      double fracTot(0.0);
      for (unsigned j = 0; j < numberOfSeeds; ++j) {
        double fraction = 0.0;
        double d2 = (std::pow(ihit.x - centroids[j].x(), 2.0) + std::pow(ihit.y - centroids[j].y(), 2.0) +
                     std::pow(ihit.z - centroids[j].z(), 2.0)) /
                    sigma2_;
        dist2[j] = d2;
        // now we set the fractions up based on hit type
        if (i == seeds[j]) {  // this cluster's seed
          fraction = 1.0;
        } else if (isaseed[i]) {
          fraction = 0.0;
        } else {
          fraction = energies[j] * std::exp(-0.5 * d2);
        }
        fracTot += fraction;
        frac[j] = fraction;
      }
      // now that we have calculated all fractions for all hits
      // assign the new fractions
      for (unsigned j = 0; j < numberOfSeeds; ++j) {
        if (fracTot > minFracTot || (i == seeds[j] && fracTot > 0.0)) {
          outclusters[j][i] = frac[j] / fracTot;
        } else {
          outclusters[j][i] = 0.0;
        }
      }
    }

    // save previous centroids
    prevCentroids = std::move(centroids);
    // finally update the position of the centroids from the last iteration
    centroids.resize(numberOfSeeds);
    double diff2 = 0.0;
    for (unsigned i = 0; i < numberOfSeeds; ++i) {
      centroids[i] = calculatePositionWithFraction(incluster, outclusters[i]);
      energies[i] = calculateEnergyWithFraction(incluster, outclusters[i]);
      // calculate convergence parameters
      const double delta2 = (prevCentroids[i] - centroids[i]).perp2();
      diff2 = std::max(delta2, diff2);
    }
    // update convergance parameter outside loop
    diff = std::sqrt(diff2);
  }
}

sort_by_delta()

std::vector<size_t> HGCalImagingAlgo::sort_by_delta ( const std::vector< KDNode > &  v ) const

inlineprivate

Definition at line 225 of file HGCalImagingAlgo.h.

References testProducerWithPsetDescEmpty_cfi::i1, testProducerWithPsetDescEmpty_cfi::i2, heavyIonCSV_trainingSettings::idx, jetUpdater_cfi::sort, and findQualityFiles::v.

                                                                    {
    std::vector<size_t> idx(v.size());
    std::iota(std::begin(idx), std::end(idx), 0);
    sort(idx.begin(), idx.end(), [&v](size_t i1, size_t i2) { return v[i1].data.delta > v[i2].data.delta; });
    return idx;
  }

Member Data Documentation

clusters_v_

std::vector<reco::BasicCluster> HGCalImagingAlgo::clusters_v_

private

Definition at line 138 of file HGCalImagingAlgo.h.

Referenced by reset().

dEdXweights_

std::vector<double> HGCalImagingAlgo::dEdXweights_

private

Definition at line 145 of file HGCalImagingAlgo.h.

density_

Density HGCalImagingAlgo::density_

private

Definition at line 141 of file HGCalImagingAlgo.h.

dependSensor_

bool HGCalImagingAlgo::dependSensor_

private

Definition at line 144 of file HGCalImagingAlgo.h.

ecut_

double HGCalImagingAlgo::ecut_

private

Definition at line 132 of file HGCalImagingAlgo.h.

fcPerEle_

double HGCalImagingAlgo::fcPerEle_

private

Definition at line 148 of file HGCalImagingAlgo.h.

fcPerMip_

std::vector<double> HGCalImagingAlgo::fcPerMip_

private

Definition at line 147 of file HGCalImagingAlgo.h.

initialized_

bool HGCalImagingAlgo::initialized_

private

Definition at line 155 of file HGCalImagingAlgo.h.

kappa_

double HGCalImagingAlgo::kappa_

private

Definition at line 129 of file HGCalImagingAlgo.h.

layerClustersPerLayer_

std::vector<std::vector<std::vector<KDNode> > > HGCalImagingAlgo::layerClustersPerLayer_

private

Definition at line 223 of file HGCalImagingAlgo.h.

Referenced by reset().

maxpos_

std::vector<std::array<float, 2> > HGCalImagingAlgo::maxpos_

private

Definition at line 236 of file HGCalImagingAlgo.h.

Referenced by reset().

minpos_

std::vector<std::array<float, 2> > HGCalImagingAlgo::minpos_

private

Definition at line 235 of file HGCalImagingAlgo.h.

Referenced by reset().

noiseMip_

double HGCalImagingAlgo::noiseMip_

private

Definition at line 150 of file HGCalImagingAlgo.h.

nonAgedNoises_

std::vector<double> HGCalImagingAlgo::nonAgedNoises_

private

Definition at line 149 of file HGCalImagingAlgo.h.

points_

std::vector<std::vector<KDNode> > HGCalImagingAlgo::points_

private

Definition at line 232 of file HGCalImagingAlgo.h.

Referenced by reset().

positionDeltaRho_c_

std::vector<double> HGCalImagingAlgo::positionDeltaRho_c_

private

Definition at line 125 of file HGCalImagingAlgo.h.

sigma2_

double HGCalImagingAlgo::sigma2_

private

Definition at line 135 of file HGCalImagingAlgo.h.

sigmaNoise_

std::vector<std::vector<double> > HGCalImagingAlgo::sigmaNoise_

private

Definition at line 152 of file HGCalImagingAlgo.h.

thicknessCorrection_

std::vector<double> HGCalImagingAlgo::thicknessCorrection_

private

Definition at line 146 of file HGCalImagingAlgo.h.

thresholds_

std::vector<std::vector<double> > HGCalImagingAlgo::thresholds_

private

Definition at line 151 of file HGCalImagingAlgo.h.

thresholdW0_

std::vector<double> HGCalImagingAlgo::thresholdW0_

private

Definition at line 124 of file HGCalImagingAlgo.h.

vecDeltas_

std::vector<double> HGCalImagingAlgo::vecDeltas_

private

Definition at line 128 of file HGCalImagingAlgo.h.