#include <HGCalImagingAlgo.h>

Classes
struct	Hexel

Public Types
typedef math::XYZPoint	Point
	point in the space More...

enum	VerbosityLevel { pDEBUG = 0, pWARNING = 1, pINFO = 2, pERROR = 3 }

Public Member Functions
void	computeThreshold ()

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

void	getEventSetup (const edm::EventSetup &es)

	HGCalImagingAlgo ()

	HGCalImagingAlgo (const std::vector< double > &vecDeltas_in, double kappa_in, double ecut_in, reco::CaloCluster::AlgoId algoId_in, bool dependSensor_in, const std::vector< double > &dEdXweights_in, const std::vector< double > &thicknessCorrection_in, const std::vector< double > &fcPerMip_in, double fcPerEle_in, const std::vector< double > &nonAgedNoises_in, double noiseMip_in, VerbosityLevel the_verbosity=pERROR)

	HGCalImagingAlgo (const std::vector< double > &vecDeltas_in, double kappa_in, double ecut_in, double showerSigma, reco::CaloCluster::AlgoId algoId_in, bool dependSensor_in, const std::vector< double > &dEdXweights_in, const std::vector< double > &thicknessCorrection_in, const std::vector< double > &fcPerMip_in, double fcPerEle_in, const std::vector< double > &nonAgedNoises_in, double noiseMip_in, VerbosityLevel the_verbosity=pERROR)

void	makeClusters ()

void	populate (const HGCRecHitCollection &hits)

void	reset ()

void	setVerbosity (VerbosityLevel the_verbosity)

virtual	~HGCalImagingAlgo ()

Static Public Attributes
static const unsigned int	maxlayer = 52

Private Types
typedef KDTreeNodeInfoT< Hexel, 2 >	KDNode

typedef KDTreeLinkerAlgo< Hexel, 2 >	KDTree

Private Member Functions
double	calculateDistanceToHigher (std::vector< KDNode > &)

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

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

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

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

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

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

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

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

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)

Private Attributes
reco::CaloCluster::AlgoId	algoId

unsigned int	cluster_offset

std::vector< reco::BasicCluster >	clusters_v

std::vector< std::vector< KDNode > >	current_v

std::vector< double >	dEdXweights

bool	dependSensor

double	ecut

double	fcPerEle

std::vector< double >	fcPerMip

bool	initialized

double	kappa

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

hgcal::RecHitTools	rhtools_

double	sigma2

std::vector< double >	thicknessCorrection

std::vector< std::vector< double > >	thresholds

std::vector< std::vector< double > >	v_sigmaNoise

std::vector< double >	vecDeltas

VerbosityLevel	verbosity

Static Private Attributes
static const unsigned int	lastLayerEE = 28

static const unsigned int	lastLayerFH = 40

static const unsigned int	maxNumberOfWafersPerLayer = 796

Detailed Description

Definition at line 44 of file HGCalImagingAlgo.h.

Member Typedef Documentation

typedef KDTreeNodeInfoT<Hexel,2> HGCalImagingAlgo::KDNode

private

Definition at line 260 of file HGCalImagingAlgo.h.

typedef KDTreeLinkerAlgo<Hexel,2> HGCalImagingAlgo::KDTree

private

Definition at line 259 of file HGCalImagingAlgo.h.

typedef math::XYZPoint HGCalImagingAlgo::Point

point in the space

Definition at line 160 of file HGCalImagingAlgo.h.

Member Enumeration Documentation

enum HGCalImagingAlgo::VerbosityLevel

Enumerator
pDEBUG
pWARNING
pINFO
pERROR

Definition at line 50 of file HGCalImagingAlgo.h.

50 { pDEBUG = 0, pWARNING = 1, pINFO = 2, pERROR = 3 };

HGCalImagingAlgo::pDEBUG

Definition: HGCalImagingAlgo.h:50

HGCalImagingAlgo::pERROR

Definition: HGCalImagingAlgo.h:50

HGCalImagingAlgo::pWARNING

Definition: HGCalImagingAlgo.h:50

HGCalImagingAlgo::pINFO

Definition: HGCalImagingAlgo.h:50

Constructor & Destructor Documentation

HGCalImagingAlgo::HGCalImagingAlgo ( )

inline

Definition at line 52 of file HGCalImagingAlgo.h.

                    : vecDeltas(), kappa(1.), ecut(0.), cluster_offset(0),
         sigma2(1.0),
         algoId(reco::CaloCluster::undefined),
         verbosity(pERROR),initialized(false){
 }

HGCalImagingAlgo::HGCalImagingAlgo	(	const std::vector< double > &	vecDeltas_in,
		double	kappa_in,
		double	ecut_in,
		reco::CaloCluster::AlgoId	algoId_in,
		bool	dependSensor_in,
		const std::vector< double > &	dEdXweights_in,
		const std::vector< double > &	thicknessCorrection_in,
		const std::vector< double > &	fcPerMip_in,
		double	fcPerEle_in,
		const std::vector< double > &	nonAgedNoises_in,
		double	noiseMip_in,
		VerbosityLevel	the_verbosity = `pERROR`
	)

inline

Definition at line 58 of file HGCalImagingAlgo.h.

References maxlayer, and maxpos.

                                                         :
         vecDeltas(vecDeltas_in), kappa(kappa_in),
         ecut(ecut_in),
         cluster_offset(0),
         sigma2(1.0),
         algoId(algoId_in),
         dependSensor(dependSensor_in),
         dEdXweights(dEdXweights_in),
         thicknessCorrection(thicknessCorrection_in),
         fcPerMip(fcPerMip_in),
         fcPerEle(fcPerEle_in),
         nonAgedNoises(nonAgedNoises_in),
         noiseMip(noiseMip_in),
         verbosity(the_verbosity),
         initialized(false),
         points(2*(maxlayer+1)),
         minpos(2*(maxlayer+1),{
                 {0.0f,0.0f}
         }),
         maxpos(2*(maxlayer+1),{ {0.0f,0.0f} })
 {
 }

HGCalImagingAlgo::HGCalImagingAlgo	(	const std::vector< double > &	vecDeltas_in,
		double	kappa_in,
		double	ecut_in,
		double	showerSigma,
		reco::CaloCluster::AlgoId	algoId_in,
		bool	dependSensor_in,
		const std::vector< double > &	dEdXweights_in,
		const std::vector< double > &	thicknessCorrection_in,
		const std::vector< double > &	fcPerMip_in,
		double	fcPerEle_in,
		const std::vector< double > &	nonAgedNoises_in,
		double	noiseMip_in,
		VerbosityLevel	the_verbosity = `pERROR`
	)

inline

Definition at line 90 of file HGCalImagingAlgo.h.

References maxlayer, and maxpos.

                                                         : vecDeltas(vecDeltas_in), kappa(kappa_in),
         ecut(ecut_in),
         cluster_offset(0),
         sigma2(std::pow(showerSigma,2.0)),
         algoId(algoId_in),
         dependSensor(dependSensor_in),
         dEdXweights(dEdXweights_in),
         thicknessCorrection(thicknessCorrection_in),
         fcPerMip(fcPerMip_in),
         fcPerEle(fcPerEle_in),
         nonAgedNoises(nonAgedNoises_in),
         noiseMip(noiseMip_in),
         verbosity(the_verbosity),
         initialized(false),
         points(2*(maxlayer+1)),
     minpos(2*(maxlayer+1),{
                 {0.0f,0.0f}
         }),
     maxpos(2*(maxlayer+1),{ {0.0f,0.0f} })
 {
 }

virtual HGCalImagingAlgo::~HGCalImagingAlgo ( )

inlinevirtual

Definition at line 122 of file HGCalImagingAlgo.h.

123 {

124 }

Member Function Documentation

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

private

Definition at line 235 of file HGCalImagingAlgo.cc.

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

Referenced by distance(), and makeClusters().

                                                                        {
 
 
   //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];
       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 455 of file HGCalImagingAlgo.cc.

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

Referenced by distance(), getClusters(), and shareEnergy().

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

private

Definition at line 209 of file HGCalImagingAlgo.cc.

References data, distance(), runEdmFileComparison::found, mps_fire::i, lastLayerEE, lastLayerFH, KDTreeLinkerAlgo< DATA, DIM >::search(), and vecDeltas.

Referenced by distance(), and makeClusters().

                                                                                                          {
 
   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;
         if(nd[i].data.rho > maxdensity) maxdensity = nd[i].data.rho;
       }
     } // end loop found
   } // end loop nodes
   return maxdensity;
 }

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

private

Definition at line 166 of file HGCalImagingAlgo.cc.

References data, mps_fire::i, x, y, and z.

Referenced by distance(), and getClusters().

                                                                     {
   float total_weight = 0.;
   float x = 0.;
   float y = 0.;
   float z = 0.;
   unsigned int v_size = v.size();
   unsigned int maxEnergyIndex = 0;
   float maxEnergyValue = 0;
   bool haloOnlyCluster = true;
 
   // loop over hits in cluster candidate building up weight for
   // energy-weighted position calculation and determining the maximum
   // energy hit in case this is a halo-only cluster
   for (unsigned int i = 0; i < v_size; i++){
     if(!v[i].data.isHalo){
       haloOnlyCluster = false;
       total_weight += v[i].data.weight;
       x += v[i].data.x*v[i].data.weight;
       y += v[i].data.y*v[i].data.weight;
       z += v[i].data.z*v[i].data.weight;
     }
     else {
       if (v[i].data.weight > maxEnergyValue) {
         maxEnergyValue = v[i].data.weight;
         maxEnergyIndex = i;
       }
     }
   }
 
   if (!haloOnlyCluster) {
     if (total_weight != 0) {
       return math::XYZPoint( x/total_weight,
              y/total_weight,
              z/total_weight );
     }
   }
   else if (v_size > 0) {
     // return position of hit with maximum energy
     return math::XYZPoint(v[maxEnergyIndex].data.x, v[maxEnergyIndex].data.y, 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 439 of file HGCalImagingAlgo.cc.

References mps_fire::i, mps_fire::result, mps_merge::weight, x, y, and z.

Referenced by distance(), getClusters(), and shareEnergy().

                                                                      {
   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);
   double norm_inv = 1.0/norm;
   result *= norm_inv;
   return result;
 }

void HGCalImagingAlgo::computeThreshold ( )

Definition at line 560 of file HGCalImagingAlgo.cc.

References dEdXweights, ecut, fcPerEle, fcPerMip, DetId::Forward, hgcal::RecHitTools::getGeometry(), hgcal::RecHitTools::getLayerWithOffset(), hgcal::RecHitTools::getSiThickness(), CaloGeometry::getValidDetIds(), hgcal::RecHitTools::getWafer(), HGCEE, HGCHEF, initialized, lastLayerFH, maxlayer, maxNumberOfWafersPerLayer, noiseMip, nonAgedNoises, rhtools_, thicknessCorrection, thresholds, v_sigmaNoise, and hgcal::RecHitTools::zside().

Referenced by populate(), and reset().

                                         {
 
   if(initialized) return; // only need to calculate thresholds once
   const std::vector<DetId>& listee(rhtools_.getGeometry()->getValidDetIds(DetId::Forward,ForwardSubdetector::HGCEE));
   const std::vector<DetId>& listfh(rhtools_.getGeometry()->getValidDetIds(DetId::Forward,ForwardSubdetector::HGCHEF));
 
   std::vector<double> dummy;
   dummy.resize(maxNumberOfWafersPerLayer, 0);
   thresholds.resize(maxlayer, dummy);
   v_sigmaNoise.resize(maxlayer, dummy);
   int previouswafer=-999;
 
   for(unsigned icalo=0;icalo<2;++icalo)
     {
       const std::vector<DetId>& listDetId( icalo==0 ? listee : listfh);
 
       for(auto& detid: listDetId)
     {
       int wafer = rhtools_.getWafer(detid);
       if(wafer==previouswafer) continue;
       previouswafer = wafer;
       // no need to do it twice
       if(rhtools_.zside(detid)<0) continue;
       int layer = rhtools_.getLayerWithOffset(detid);
       float thickness = rhtools_.getSiThickness(detid);
       int thickIndex = -1;
       if( thickness>99. && thickness<101.) thickIndex=0;
       else if( thickness>199. && thickness<201. ) thickIndex=1;
       else if( thickness>299. && thickness<301. ) thickIndex=2;
       else assert( thickIndex>0 && "ERROR - silicon thickness has a nonsensical value" );
       float sigmaNoise = 0.001 * fcPerEle * nonAgedNoises[thickIndex] * dEdXweights[layer] / (fcPerMip[thickIndex] * thicknessCorrection[thickIndex]);
       thresholds[layer-1][wafer]=sigmaNoise*ecut;
       v_sigmaNoise[layer-1][wafer] = sigmaNoise;
     }
     }
 
   // now BH, much faster
   for ( unsigned ilayer=lastLayerFH+1;ilayer<=maxlayer;++ilayer)
     {
       float sigmaNoise = 0.001 * noiseMip * dEdXweights[ilayer];
       std::vector<double> bhDummy_thresholds;
       std::vector<double> bhDummy_sigmaNoise;
       bhDummy_thresholds.push_back(sigmaNoise*ecut);
       bhDummy_sigmaNoise.push_back(sigmaNoise);
       thresholds[ilayer-1]=bhDummy_thresholds;
       v_sigmaNoise[ilayer-1]=bhDummy_sigmaNoise;
     }
   initialized=true;
 }

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

inlineprivate

Definition at line 289 of file HGCalImagingAlgo.h.

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

Referenced by calculateLocalDensity(), findAndAssignClusters(), and findLocalMaximaInCluster().

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

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

inlineprivate

Definition at line 284 of file HGCalImagingAlgo.h.

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

Referenced by calculateDistanceToHigher(), and distance().

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

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

private

Definition at line 286 of file HGCalImagingAlgo.cc.

References KDTreeLinkerAlgo< DATA, DIM >::build(), KDTreeLinkerAlgo< DATA, DIM >::clear(), cluster_offset, gather_cfg::cout, current_v, data, dependSensor, distance(), runEdmFileComparison::found, mps_fire::i, kappa, lastLayerEE, lastLayerFH, pINFO, KDTreeLinkerAlgo< DATA, DIM >::search(), sort_by_delta(), sorted_indices(), vecDeltas, and verbosity.

Referenced by distance(), and makeClusters().

                                                                                                                                            {
 
   //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 clusterIndex = 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 = clusterIndex;
     if (verbosity < pINFO)
       {
         std::cout << "Adding new cluster with index " << clusterIndex+cluster_offset << std::endl;
         std::cout << "Cluster center is hit " << ds[i] << std::endl;
       }
     clusterIndex++;
   }
 
   //at this point clusterIndex is equal to the number of cluster centers - if it is zero we are
   //done
   if(clusterIndex==0) return clusterIndex;
 
   //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 "<< clusterIndex << std::endl;
     }
   current_v.resize(cluster_offset+clusterIndex);
 
   //assign points closer than dc to other clusters to border region
   //and find critical border density
   std::vector<double> rho_b(clusterIndex,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;
       current_v[ci+cluster_offset].push_back(nd[i]);
       if (verbosity < pINFO)
       {
         std::cout << "Pushing hit " << i << " into cluster with index " << ci+cluster_offset << std::endl;
         std::cout << "Size now " << current_v[ci+cluster_offset].size() << std::endl;
       }
     }
   }
 
   //prepare the offset for the next layer if there is one
   if (verbosity < pINFO)
     {
       std::cout << "moving cluster offset by " << clusterIndex << std::endl;
     }
   cluster_offset += clusterIndex;
   return clusterIndex;
 }

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

private

Definition at line 411 of file HGCalImagingAlgo.cc.

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

Referenced by distance(), and getClusters().

                                                                                                {
   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( distance(cluster[i].data,cluster[j].data) < delta_c && i != j) {
     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 )

Definition at line 92 of file HGCalImagingAlgo.cc.

References algoId, calculateEnergyWithFraction(), calculatePosition(), calculatePositionWithFraction(), clusters_v, gather_cfg::cout, current_v, data, reco::CaloID::DET_HGCAL_ENDCAP, MillePedeFileConverter_cfg::e, findLocalMaximaInCluster(), dedxEstimators_cff::fraction, mps_fire::i, pINFO, position, shareEnergy(), and verbosity.

Referenced by setVerbosity().

                                                                        {
 
   reco::CaloID caloID = reco::CaloID::DET_HGCAL_ENDCAP;
   std::vector< std::pair<DetId, float> > thisCluster;
   for (unsigned int i = 0; i < current_v.size(); ++i){
     double energy = 0;
     Point position;
 
     if( doSharing ) {
 
       std::vector<unsigned> seeds = findLocalMaximaInCluster(current_v[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(current_v[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(current_v[i],fractions[isub]);
     Point position = calculatePositionWithFraction(current_v[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(current_v[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.push_back(reco::BasicCluster(energy, position, caloID, thisCluster,
                         algoId));
     thisCluster.clear();
       }
     }else{
       position = calculatePosition(current_v[i]); // energy-weighted position
     //   std::vector< KDNode >::iterator it;
       for (auto& it: current_v[i])
       {
         energy += it.data.isHalo ? 0. : it.data.weight;
       // use fraction to store whether this is a Halo hit or not
       thisCluster.emplace_back(std::pair<DetId, float>(it.data.detid,(it.data.isHalo?0.:1.)));
     };
       if (verbosity < pINFO)
     {
       std::cout << "******** NEW CLUSTER (HGCIA) ********" << std::endl;
       std::cout << "Index          " << i                   << std::endl;
       std::cout << "No. of cells = " << current_v[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.push_back(reco::BasicCluster(energy, position, caloID, thisCluster,
                           algoId));
       thisCluster.clear();
     }
   }
   return clusters_v;
 }

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

inline

Definition at line 138 of file HGCalImagingAlgo.h.

References hgcal::RecHitTools::getEventSetup(), and rhtools_.

                                            {
         rhtools_.getEventSetup(es);
 }

void HGCalImagingAlgo::makeClusters ( )

Definition at line 69 of file HGCalImagingAlgo.cc.

References calculateDistanceToHigher(), calculateLocalDensity(), findAndAssignClusters(), mps_fire::i, maxlayer, maxpos, minpos, and points.

Referenced by setVerbosity().

 {
 
   // used for speedy search
   std::vector<KDTree> hit_kdtree(2*(maxlayer+1));
 
   //assign all hits in each layer to a cluster core or halo
   for (unsigned int i = 0; i <= 2*maxlayer+1; ++i) {
     KDTreeBox bounds(minpos[i][0],maxpos[i][0],
              minpos[i][1],maxpos[i][1]);
 
     hit_kdtree[i].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[i], 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[i],maxdensity,bounds,actualLayer);
   }
   //make the cluster vector
 }

void HGCalImagingAlgo::populate ( const HGCRecHitCollection & hits )

Definition at line 13 of file HGCalImagingAlgo.cc.

References computeThreshold(), dependSensor, CaloRecHit::detid(), ecut, CaloRecHit::energy(), hgcal::RecHitTools::getLayerWithOffset(), hgcal::RecHitTools::getPosition(), hgcal::RecHitTools::getSiThickness(), hgcal::RecHitTools::getWafer(), mps_fire::i, createfilelist::int, hgcal::RecHitTools::isHalfCell(), lastLayerFH, hpstanc_transforms::max, maxlayer, maxpos, min(), minpos, eostools::move(), points, position, rhtools_, edm::SortedCollection< T, SORT >::size(), thresholds, v_sigmaNoise, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and hgcal::RecHitTools::zside().

Referenced by setVerbosity().

                                                               {
   //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.;
     // set sigmaNoise default value 1 to use kappa value directly in case of sensor-independent thresholds
     float sigmaNoise = 1.;
     if(dependSensor){
       if (layer<= lastLayerFH) // only EE and FH have silicon sensors
         thickness = rhtools_.getSiThickness(detid);
       double storedThreshold=thresholds[layer-1][layer<=lastLayerFH ? rhtools_.getWafer(detid) : 0];
       sigmaNoise = v_sigmaNoise[layer-1][layer<=lastLayerFH ? rhtools_.getWafer(detid) : 0];
 
       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( std::move( 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 )

inline

Definition at line 142 of file HGCalImagingAlgo.h.

References cluster_offset, clusters_v, computeThreshold(), current_v, mps_fire::i, maxpos, minpos, points, and edm::swap().

             {
         current_v.clear();
         clusters_v.clear();
         cluster_offset = 0;
         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::setVerbosity ( VerbosityLevel the_verbosity )

inline

Definition at line 126 of file HGCalImagingAlgo.h.

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

 {
         verbosity = the_verbosity;
 }

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

private

Definition at line 464 of file HGCalImagingAlgo.cc.

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

Referenced by distance(), and getClusters().

                                                                {
   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 double toleranceScaling = std::pow(std::max(1.0,seeds.size()-1.0),2.0);
   std::vector<Point> prevCentroids;
   std::vector<double> frac(seeds.size()), dist2(seeds.size());
   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 < seeds.size(); ++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 < seeds.size(); ++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(seeds.size());
     double diff2 = 0.0;
     for( unsigned i = 0; i < seeds.size(); ++i ) {
       centroids[i] = calculatePositionWithFraction(incluster,outclusters[i]);
       energies[i]  = calculateEnergyWithFraction(incluster,outclusters[i]);
       // calculate convergence parameters
       const double delta2 = (prevCentroids[i]-centroids[i]).perp2();
       if( delta2 > diff2 ) diff2 = delta2;
     }
     //update convergance parameter outside loop
     diff = std::sqrt(diff2);
   }
 }

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

inlineprivate

Definition at line 266 of file HGCalImagingAlgo.h.

References begin, end, training_settings::idx, and findQualityFiles::v.

Referenced by findAndAssignClusters().

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