#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

	HGCalImagingAlgo (const edm::ParameterSet &ps)

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)

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

void	setVerbosity (VerbosityLevel the_verbosity)

virtual	~HGCalClusteringAlgoBase ()

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

Private Types
typedef KDTreeNodeInfoT< Hexel, 2 >	KDNode

typedef KDTreeLinkerAlgo< Hexel, 2 >	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 &, 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
Static Public Attributes inherited from HGCalClusteringAlgoBase
static const unsigned int	lastLayerEE = 28

static const unsigned int	lastLayerFH = 40

static const unsigned int	maxlayer = 52

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 31 of file HGCalImagingAlgo.h.

Member Typedef Documentation

typedef KDTreeNodeInfoT<Hexel,2> HGCalImagingAlgo::KDNode

private

Definition at line 209 of file HGCalImagingAlgo.h.

typedef KDTreeLinkerAlgo<Hexel,2> HGCalImagingAlgo::KDTree

private

Definition at line 208 of file HGCalImagingAlgo.h.

typedef math::XYZPoint HGCalImagingAlgo::Point

point in the space

Definition at line 122 of file HGCalImagingAlgo.h.

Constructor & Destructor Documentation

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

inline

Definition at line 35 of file HGCalImagingAlgo.h.

References HGCalClusteringAlgoBase::maxlayer, and maxpos_.

   : HGCalClusteringAlgoBase(
       (HGCalClusteringAlgoBase::VerbosityLevel)ps.getUntrackedParameter<unsigned int>("verbosity",3),
       reco::CaloCluster::undefined),
      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>("value")),
      initialized_(false),
      points_(2*(maxlayer+1)),
      minpos_(2*(maxlayer+1),{ {0.0f,0.0f} }),
      maxpos_(2*(maxlayer+1),{ {0.0f,0.0f} }) {}

HGCalImagingAlgo::~HGCalImagingAlgo ( )

inlineoverride

Definition at line 57 of file HGCalImagingAlgo.h.

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

57 {}

Member Function Documentation

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

private

Definition at line 299 of file HGCalImagingAlgo.cc.

References data, mps_fire::i, hgcal_clustering::sorted_indices(), mathSSE::sqrt(), and tmp.

Referenced by distance().

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

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

private

Definition at line 541 of file HGCalImagingAlgo.cc.

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

Referenced by distance().

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

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

private

Definition at line 266 of file HGCalImagingAlgo.cc.

References data, SoftLeptonByDistance_cfi::distance, runEdmFileComparison::found, mps_fire::i, SiStripPI::max, and KDTreeLinkerAlgo< DATA >::search().

Referenced by distance().

                                                                                {
 
   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 <= lastLayerFH)
     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<KDNode> 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].data) < delta_c) {
         nd[i].data.rho += found[j].data.weight;
         maxdensity = std::max(maxdensity, nd[i].data.rho);
       }
     } // end loop found
   }   // end loop nodes
   return maxdensity;
 }

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

private

Definition at line 198 of file HGCalImagingAlgo.cc.

References data, mps_fire::i, training_settings::idx, cmsBatch::log, and SiStripPI::max.

Referenced by distance().

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

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

private

Definition at line 526 of file HGCalImagingAlgo.cc.

References mps_fire::i, mps_fire::result, and mps_merge::weight.

Referenced by distance().

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

void HGCalImagingAlgo::computeThreshold ( )

Definition at line 650 of file HGCalImagingAlgo.cc.

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

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

inlineprivate

Definition at line 237 of file HGCalImagingAlgo.h.

References calculateDistanceToHigher(), calculateEnergyWithFraction(), calculateLocalDensity(), calculatePosition(), calculatePositionWithFraction(), distance2(), findAndAssignClusters(), findLocalMaximaInCluster(), setDensity(), shareEnergy(), and mathSSE::sqrt().

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

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

inlineprivate

Definition at line 232 of file HGCalImagingAlgo.h.

References PVValHelper::dx, PVValHelper::dy, HGCalImagingAlgo::Hexel::x, and HGCalImagingAlgo::Hexel::y.

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

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

inlinestatic

Definition at line 89 of file HGCalImagingAlgo.h.

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

                                                                    {
   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<double>("value", 0 );
   iDesc.add<edm::ParameterSetDescription>("noiseMip", descNestedNoiseMIP);
 }

int HGCalImagingAlgo::findAndAssignClusters	(	std::vector< KDNode > &	,
		KDTree &	,
		double	,
		KDTreeBox &	,
		const unsigned int	,
		std::vector< std::vector< KDNode > > &
	)		const

private

Definition at line 353 of file HGCalImagingAlgo.cc.

References KDTreeLinkerAlgo< DATA >::build(), KDTreeLinkerAlgo< DATA >::clear(), gather_cfg::cout, data, SoftLeptonByDistance_cfi::distance, runEdmFileComparison::found, mps_fire::i, KDTreeLinkerAlgo< DATA >::search(), and hgcal_clustering::sorted_indices().

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;
   float delta_c; // critical distance
   if (layer <= lastLayerEE)
     delta_c = vecDeltas_[0];
   else if (layer <= lastLayerFH)
     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<KDNode> 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].data.clusterIndex != -1) {
           float dist = distance(found[j].data, nd[i].data);
           if (dist < delta_c && found[j].data.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].data.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;
 }

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

private

Definition at line 498 of file HGCalImagingAlgo.cc.

References data, SoftLeptonByDistance_cfi::distance, MillePedeFileConverter_cfg::e, mps_fire::i, mps_fire::result, and SurveyInfoScenario_cff::seed.

Referenced by distance().

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

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

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 114 of file HGCalImagingAlgo.cc.

References AlignmentPI::calculatePosition(), gather_cfg::cout, data, reco::CaloID::DET_HGCAL_ENDCAP, MillePedeFileConverter_cfg::e, randomXiThetaGunProducer_cfi::energy, dedxEstimators_cff::fraction, mps_fire::i, hgcal_clustering::max_index(), and position.

Referenced by ~HGCalImagingAlgo().

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

Density HGCalImagingAlgo::getDensity ( )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 693 of file HGCalImagingAlgo.cc.

Referenced by reset().

                                     {
   return density_;
 }

void HGCalImagingAlgo::makeClusters ( )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 86 of file HGCalImagingAlgo.cc.

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

Referenced by ~HGCalImagingAlgo().

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

void HGCalImagingAlgo::populate ( const HGCRecHitCollection & hits )

overridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 17 of file HGCalImagingAlgo.cc.

References CaloRecHit::detid(), CaloRecHit::energy(), mps_fire::i, createfilelist::int, SiStripPI::max, min(), position, and edm::SortedCollection< T, SORT >::size().

Referenced by ~HGCalImagingAlgo().

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

void HGCalImagingAlgo::reset ( void )

inlineoverridevirtual

Implements HGCalClusteringAlgoBase.

Definition at line 70 of file HGCalImagingAlgo.h.

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

                       {
         clusters_v_.clear();
         layerClustersPerLayer_.clear();
         for( auto& it: points_)
         {
                 it.clear();
                 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.;
         }
 }

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

private

Definition at line 684 of file HGCalImagingAlgo.cc.

References mps_fire::i.

Referenced by distance().

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

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

private

Definition at line 550 of file HGCalImagingAlgo.cc.

References data, diffTreeTool::diff, MillePedeFileConverter_cfg::e, JetChargeProducer_cfi::exp, cropTnPTrees::frac, dedxEstimators_cff::fraction, mps_fire::i, SiStripPI::max, particleFlowClusterECALUncorrected_cfi::minFracTot, eostools::move(), perp2(), funct::pow(), mathSSE::sqrt(), particleFlowClusterECALUncorrected_cfi::stoppingTolerance, HGCalImagingAlgo::Hexel::x, HGCalImagingAlgo::Hexel::y, and HGCalImagingAlgo::Hexel::z.

Referenced by distance().

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

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

inlineprivate

Definition at line 214 of file HGCalImagingAlgo.h.

References begin, end, training_settings::idx, 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

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

private

Definition at line 141 of file HGCalImagingAlgo.h.

Referenced by reset().

std::vector<double> HGCalImagingAlgo::dEdXweights_

private

Definition at line 148 of file HGCalImagingAlgo.h.

Density HGCalImagingAlgo::density_

private

Definition at line 144 of file HGCalImagingAlgo.h.

bool HGCalImagingAlgo::dependSensor_

private

Definition at line 147 of file HGCalImagingAlgo.h.

double HGCalImagingAlgo::ecut_

private

Definition at line 135 of file HGCalImagingAlgo.h.

double HGCalImagingAlgo::fcPerEle_

private

Definition at line 151 of file HGCalImagingAlgo.h.

std::vector<double> HGCalImagingAlgo::fcPerMip_

private

Definition at line 150 of file HGCalImagingAlgo.h.

bool HGCalImagingAlgo::initialized_

private

Definition at line 158 of file HGCalImagingAlgo.h.

double HGCalImagingAlgo::kappa_

private

Definition at line 132 of file HGCalImagingAlgo.h.

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

private

Definition at line 212 of file HGCalImagingAlgo.h.

Referenced by reset().

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

private

Definition at line 228 of file HGCalImagingAlgo.h.

Referenced by HGCalImagingAlgo(), and reset().

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

private

Definition at line 227 of file HGCalImagingAlgo.h.

Referenced by reset().

double HGCalImagingAlgo::noiseMip_

private

Definition at line 153 of file HGCalImagingAlgo.h.

std::vector<double> HGCalImagingAlgo::nonAgedNoises_

private

Definition at line 152 of file HGCalImagingAlgo.h.

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

private

Definition at line 224 of file HGCalImagingAlgo.h.

Referenced by reset().

std::vector<double> HGCalImagingAlgo::positionDeltaRho_c_

private

Definition at line 128 of file HGCalImagingAlgo.h.

double HGCalImagingAlgo::sigma2_

private

Definition at line 138 of file HGCalImagingAlgo.h.

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

private

Definition at line 155 of file HGCalImagingAlgo.h.

std::vector<double> HGCalImagingAlgo::thicknessCorrection_

private

Definition at line 149 of file HGCalImagingAlgo.h.

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

private

Definition at line 154 of file HGCalImagingAlgo.h.

std::vector<double> HGCalImagingAlgo::thresholdW0_

private

Definition at line 127 of file HGCalImagingAlgo.h.

std::vector<double> HGCalImagingAlgo::vecDeltas_

private

Definition at line 131 of file HGCalImagingAlgo.h.

Classes

Public Types

Public Member Functions

Static Public Member Functions

Private Types

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation