HGCalImagingAlgo Class Reference

#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
std::vector< reco::BasicCluster >	getClusters (bool)
void	getEventSetup (const edm::EventSetup &es)
	HGCalImagingAlgo ()
	HGCalImagingAlgo (float delta_c_in, double kappa_in, double ecut_in, reco::CaloCluster::AlgoId algoId_in, VerbosityLevel the_verbosity=pERROR)
	HGCalImagingAlgo (float delta_c_in, double kappa_in, double ecut_in, double showerSigma, reco::CaloCluster::AlgoId algoId_in, VerbosityLevel the_verbosity=pERROR)
void	makeClusters (const HGCRecHitCollection &hits)
void	reset ()
void	setVerbosity (VerbosityLevel the_verbosity)
virtual	~HGCalImagingAlgo ()

Private Types
typedef KDTreeNodeInfoT< Hexel, 2 >	KDNode
typedef KDTreeLinkerAlgo< Hexel, 2 >	KDTree

Private Member Functions
double	calculateDistanceToHigher (std::vector< KDNode > &, KDTree &)
double	calculateEnergyWithFraction (const std::vector< KDNode > &, const std::vector< double > &)
double	calculateLocalDensity (std::vector< KDNode > &, KDTree &)
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 &)
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
float	delta_c
double	ecut
double	kappa
std::vector< std::vector< Hexel > >	points
hgcal::RecHitTools	rhtools_
double	sigma2
VerbosityLevel	verbosity

Static Private Attributes
static const unsigned int	maxlayer = 52

Detailed Description

Definition at line 42 of file HGCalImagingAlgo.h.

Member Typedef Documentation

typedef KDTreeNodeInfoT<Hexel,2> HGCalImagingAlgo::KDNode

private

Definition at line 179 of file HGCalImagingAlgo.h.

typedef KDTreeLinkerAlgo<Hexel,2> HGCalImagingAlgo::KDTree

private

Definition at line 178 of file HGCalImagingAlgo.h.

typedef math::XYZPoint HGCalImagingAlgo::Point

point in the space

Definition at line 102 of file HGCalImagingAlgo.h.

Member Enumeration Documentation

enum HGCalImagingAlgo::VerbosityLevel

Enumerator
pDEBUG
pWARNING
pINFO
pERROR

Definition at line 48 of file HGCalImagingAlgo.h.

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

Constructor & Destructor Documentation

HGCalImagingAlgo::HGCalImagingAlgo ( )

inline

Definition at line 50 of file HGCalImagingAlgo.h.

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

HGCalImagingAlgo::HGCalImagingAlgo ( float  delta_c_in, double  kappa_in, double  ecut_in, reco::CaloCluster::AlgoId  algoId_in, VerbosityLevel  the_verbosity = pERROR  )

inline

Definition at line 56 of file HGCalImagingAlgo.h.

                                                  : delta_c(delta_c_in), kappa(kappa_in), 
                                ecut(ecut_in),    
                                cluster_offset(0),
                                sigma2(1.0),
                                algoId(algoId_in),
                                verbosity(the_verbosity){
  }

HGCalImagingAlgo::HGCalImagingAlgo ( float  delta_c_in, double  kappa_in, double  ecut_in, double  showerSigma, reco::CaloCluster::AlgoId  algoId_in, VerbosityLevel  the_verbosity = pERROR  )

inline

Definition at line 67 of file HGCalImagingAlgo.h.

                                                  : delta_c(delta_c_in), kappa(kappa_in), 
                                ecut(ecut_in),    
                                cluster_offset(0),
                                sigma2(std::pow(showerSigma,2.0)),
                                algoId(algoId_in),
                                verbosity(the_verbosity){
  }

virtual HGCalImagingAlgo::~HGCalImagingAlgo ( )

inlinevirtual

Definition at line 79 of file HGCalImagingAlgo.h.

    {
    }

Member Function Documentation

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

private

Definition at line 191 of file HGCalImagingAlgo.cc.

References data, distance2(), i, j, 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.size()>0) 
    maxdensity = nd[rs[0]].data.rho;
  else
    return maxdensity; // there are no hits
  double dist2 = 2500.0;
  //start by setting delta for the highest density hit to 
  //the most distant hit - this is a convention

  for(unsigned int j = 0; j < nd.size(); j++){
    double tmp = distance2(nd[rs[0]].data, nd[j].data);
    dist2 = tmp > dist2 ? tmp : dist2;
  }
  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;
  
  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 399 of file HGCalImagingAlgo.cc.

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

private

Definition at line 174 of file HGCalImagingAlgo.cc.

References data, delta_c, distance(), runEdmFileComparison::found, i, j, and KDTreeLinkerAlgo< DATA, DIM >::search().

Referenced by distance(), and makeClusters().

                                                                               {
  double maxdensity = 0.;
  for(unsigned int i = 0; i < nd.size(); ++i){
    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);
    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;
      }
    }
  }
  return maxdensity;
}

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

private

Definition at line 155 of file HGCalImagingAlgo.cc.

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

Referenced by distance(), and getClusters().

                                                                    {
  float total_weight = 0.;
  float x = 0.;
  float y = 0.;
  float z = 0.;
  for (unsigned int i = 0; i < v.size(); i++){
    if(!v[i].data.isHalo){
      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;
    }
  }
  
  return math::XYZPoint( x/total_weight, 
             y/total_weight, 
             z/total_weight );
} 

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

private

Definition at line 383 of file HGCalImagingAlgo.cc.

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

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

inlineprivate

Definition at line 202 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 197 of file HGCalImagingAlgo.h.

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

private

Definition at line 241 of file HGCalImagingAlgo.cc.

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

Referenced by distance(), and makeClusters().

                                                                                                                 {

  //this is called once per layer...
  //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;

  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


  for(unsigned int i =0; i < nd.size(); ++i){

    //    std::cout << " delta " << lp[ds[i]].delta << " rho " << lp[ds[i]].rho << std::endl;
    if(nd[ds[i]].data.delta < delta_c) break; // no more cluster centers to be looked at 
    if(nd[ds[i]].data.rho < maxdensity/kappa  /* || lp[ds[i]].rho<0.001*/) continue; 
    //skip this as a potential cluster center because it fails the density cut

    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 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){
      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);
    
      for(unsigned int j = 1; j < found.size(); j++){
    //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){
        //this is not an isolated hit
        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;
  }

  //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 && nd[i].data.rho < rho_b[ci])
      nd[i].data.isHalo = true;
    if(nd[i].data.clusterIndex!=-1){ 
      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 355 of file HGCalImagingAlgo.cc.

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

Referenced by distance(), and getClusters().

                                                                                               {
  std::vector<unsigned> result;
  std::vector<bool> seed(cluster.size(),true);
 
  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) {
    //std::cout << "hit-to-hit distance = " << distance(cluster[i],cluster[j]) << std::endl;
    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) {
      //std::cout << "seed at " << i << " with energy " << cluster[i].weight << std::endl;
      result.push_back(i);
    }
  }

  //std::cout << "Found " << result.size() << " sub-clusters in input cluster of length: " << cluster.size() << std::endl;

  return result;
}

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

(

bool

doSharing

)

Definition at line 79 of file HGCalImagingAlgo.cc.

References algoId, begin, calculateEnergyWithFraction(), calculatePosition(), calculatePositionWithFraction(), clusters_v, gather_cfg::cout, current_v, data, reco::CaloID::DET_HGCAL_ENDCAP, MillePedeFileConverter_cfg::e, findLocalMaximaInCluster(), dedxEstimators_cff::fraction, 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]);
      //std::cout << " sharing found " << seeds.size() << " sub-cluster seeds in cluster " << i << std::endl;
      
      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]);

    //std::cout << "Fractions*Energies: ";
    for( unsigned ihit = 0; ihit < fractions[isub].size(); ++ihit ) {
      const double fraction = fractions[isub][ihit];
      if( fraction > 1e-7 ) {
        //std::cout << fraction << "*" << current_v[i][ihit].weight << " ";
        effective_hits += fraction;
        thisCluster.emplace_back(current_v[i][ihit].data.detid,fraction);
      }
    }
    //std::cout << std::endl;
    
    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]);    
      std::vector< KDNode >::iterator it;
      for (it = current_v[i].begin(); it != current_v[i].end(); it++)
    {
      energy += (*it).data.isHalo ? 0. : (*it).data.weight;
      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 94 of file HGCalImagingAlgo.h.

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

{ rhtools_.getEventSetup(es); }

void HGCalImagingAlgo::makeClusters

(

const HGCRecHitCollection &

hits

)

Definition at line 16 of file HGCalImagingAlgo.cc.

References calculateDistanceToHigher(), calculateLocalDensity(), gather_cfg::cout, delta_c, CaloRecHit::detid(), ecut, CaloRecHit::energy(), findAndAssignClusters(), hgcal::RecHitTools::getLayerWithOffset(), hgcal::RecHitTools::getPosition(), i, createfilelist::int, hgcal::RecHitTools::isHalfCell(), kappa, hpstanc_transforms::max, maxlayer, min(), eostools::move(), pINFO, points, position, rhtools_, edm::SortedCollection< T, SORT >::size(), findQualityFiles::size, verbosity, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and HcalDigiParam_cfi::zside.

Referenced by setVerbosity().

{
  //used for speedy search 
  std::vector<std::vector<KDNode> >points(2*(maxlayer+1));
  std::vector<KDTree> hit_kdtree(2*(maxlayer+1));

  std::vector<std::array<float,2> > minpos(2*(maxlayer+1),{ {0.0f,0.0f} }),maxpos(2*(maxlayer+1),{ {0.0f,0.0f} });

  //  std::vector<std::vector<Hexel> > points(2*(maxlayer+1)); //a vector of vectors of hexels, one for each layer
  //@@EM todo: the number of layers should be obtained programmatically - the range is 1-n instead of 0-n-1...


  if (verbosity < pINFO)
    {
      std::cout << "-------------------------------------------------------------" << std::endl;
      std::cout << "HGC Imaging algorithm invoked for " << std::endl;
      std::cout << "delta_c " << delta_c << " kappa " << kappa;
      //      if( doSharing ) std::cout << " showerSigma " << std::sqrt(sigma2);
      std::cout << std::endl;
    }

  //loop over all hits and create the Hexel structure, skip energies below ecut
  for (unsigned int i=0;i<hits.size();++i) {
    const HGCRecHit& hgrh = hits[i];
    if(hgrh.energy() < ecut) continue; 
    DetId detid = hgrh.detid();

    int layer = rhtools_.getLayerWithOffset(detid)+int(HGCalDetId(detid).zside()>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,&rhtools_),position.x(),position.y());
    if(points[layer].size()==0){
      minpos[layer][0] = position.x(); minpos[layer][1] = position.y();
      maxpos[layer][0] = position.x(); maxpos[layer][1] = position.y();
    }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]);
      
    }
  }

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

    double maxdensity = calculateLocalDensity(points[i],hit_kdtree[i]);
    // std::cout << "layer " << i << " max density " << maxdensity 
    //        << " total hits " << points[i].size() << std::endl;
    calculateDistanceToHigher(points[i],hit_kdtree[i]);
    findAndAssignClusters(points[i],hit_kdtree[i],maxdensity,bounds);
    //    std::cout << "found " << nclusters << " clusters" << std::endl;
  }
  //make the cluster vector
}

void HGCalImagingAlgo::reset ( void  )

inline

Definition at line 96 of file HGCalImagingAlgo.h.

References cluster_offset, clusters_v, and current_v.

              {
    current_v.clear();
    clusters_v.clear();
    cluster_offset = 0;
  }

void HGCalImagingAlgo::setVerbosity ( VerbosityLevel  the_verbosity )

inline

Definition at line 83 of file HGCalImagingAlgo.h.

References getClusters(), hfClusterShapes_cfi::hits, makeClusters(), 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 408 of file HGCalImagingAlgo.cc.

References calculateEnergyWithFraction(), calculatePositionWithFraction(), data, mps_update::diff, MillePedeFileConverter_cfg::e, JetChargeProducer_cfi::exp, cropTnPTrees::frac, dedxEstimators_cff::fraction, i, j, 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;
  }

  //std::cout << "saving seeds" << std::endl;

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

  //std::cout << "seeds saved" << std::endl;

  // initialize clusters to be shared
  // centroids start off at seed positions
  // seeds always have fraction 1.0, to stabilize fit
  //std::cout << "initializing fit" << std::endl;
  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; 
      } 
    }
  }

  //return;

  //std::cout << "fit initialized" << std::endl;

  // 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 fraction(0.0), fracTot(0.0), d2(0.0);
      for( unsigned j = 0; j < seeds.size(); ++j ) {
    fraction = 0.0;
    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::cout << " iteration = " << iter << " diff = " << diff << std::endl;
  } 
}

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

inlineprivate

Definition at line 185 of file HGCalImagingAlgo.h.

References i, training_settings::idx, and findQualityFiles::v.

Referenced by findAndAssignClusters().

                                                             {
    std::vector<size_t> idx(v.size());
    for (size_t i = 0; i != idx.size(); ++i) idx[i] = i;
    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

reco::CaloCluster::AlgoId HGCalImagingAlgo::algoId

private

Definition at line 128 of file HGCalImagingAlgo.h.

Referenced by getClusters().

unsigned int HGCalImagingAlgo::cluster_offset

private

Definition at line 117 of file HGCalImagingAlgo.h.

Referenced by findAndAssignClusters(), and reset().

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

private

Definition at line 123 of file HGCalImagingAlgo.h.

Referenced by getClusters(), and reset().

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

private

Definition at line 183 of file HGCalImagingAlgo.h.

Referenced by findAndAssignClusters(), getClusters(), and reset().

float HGCalImagingAlgo::delta_c

private

Definition at line 110 of file HGCalImagingAlgo.h.

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

double HGCalImagingAlgo::ecut

private

Definition at line 114 of file HGCalImagingAlgo.h.

Referenced by makeClusters().

double HGCalImagingAlgo::kappa

private

Definition at line 111 of file HGCalImagingAlgo.h.

Referenced by findAndAssignClusters(), and makeClusters().

const unsigned int HGCalImagingAlgo::maxlayer = 52

staticprivate

Definition at line 107 of file HGCalImagingAlgo.h.

Referenced by makeClusters().

std::vector<std::vector<Hexel> > HGCalImagingAlgo::points

private

Definition at line 193 of file HGCalImagingAlgo.h.

Referenced by makeClusters().

hgcal::RecHitTools HGCalImagingAlgo::rhtools_

private

Definition at line 125 of file HGCalImagingAlgo.h.

Referenced by getEventSetup(), and makeClusters().

double HGCalImagingAlgo::sigma2

private

Definition at line 120 of file HGCalImagingAlgo.h.

Referenced by shareEnergy().

VerbosityLevel HGCalImagingAlgo::verbosity

private

Definition at line 131 of file HGCalImagingAlgo.h.

Referenced by findAndAssignClusters(), getClusters(), makeClusters(), and setVerbosity().