test
CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Classes | Public Types | Public Member Functions | Private Types | Private Member Functions | Private Attributes | Static Private Attributes
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::BasicClustergetClusters (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)
 
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::BasicClusterclusters_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

Definition at line 179 of file HGCalImagingAlgo.h.

Definition at line 178 of file HGCalImagingAlgo.h.

point in the space

Definition at line 102 of file HGCalImagingAlgo.h.

Member Enumeration Documentation

Enumerator
pDEBUG 
pWARNING 
pINFO 
pERROR 

Definition at line 48 of file HGCalImagingAlgo.h.

Constructor & Destructor Documentation

HGCalImagingAlgo::HGCalImagingAlgo ( )
inline
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.

59  : delta_c(delta_c_in), kappa(kappa_in),
60  ecut(ecut_in),
61  cluster_offset(0),
62  sigma2(1.0),
63  algoId(algoId_in),
64  verbosity(the_verbosity){
65  }
reco::CaloCluster::AlgoId algoId
unsigned int cluster_offset
VerbosityLevel 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.

71  : delta_c(delta_c_in), kappa(kappa_in),
72  ecut(ecut_in),
73  cluster_offset(0),
74  sigma2(std::pow(showerSigma,2.0)),
75  algoId(algoId_in),
76  verbosity(the_verbosity){
77  }
reco::CaloCluster::AlgoId algoId
unsigned int cluster_offset
VerbosityLevel verbosity
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
virtual HGCalImagingAlgo::~HGCalImagingAlgo ( )
inlinevirtual

Definition at line 79 of file HGCalImagingAlgo.h.

80  {
81  }

Member Function Documentation

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

Definition at line 198 of file HGCalImagingAlgo.cc.

References data, distance(), i, j, sorted_indices(), and tmp.

Referenced by makeClusters().

198  {
199 
200 
201  //sort vector of Hexels by decreasing local density
202  std::vector<size_t> rs = sorted_indices(nd);
203 
204  double maxdensity = 0.0;
205  int nearestHigher = -1;
206 
207 
208  if(rs.size()>0)
209  maxdensity = nd[rs[0]].data.rho;
210  else
211  return maxdensity; // there are no hits
212  double dist = 50.0;
213  //start by setting delta for the highest density hit to
214  //the most distant hit - this is a convention
215 
216  for(unsigned int j = 0; j < nd.size(); j++){
217  double tmp = distance(nd[rs[0]].data, nd[j].data);
218  dist = tmp > dist ? tmp : dist;
219  }
220  nd[rs[0]].data.delta = dist;
221  nd[rs[0]].data.nearestHigher = nearestHigher;
222 
223  //now we save the largest distance as a starting point
224 
225  double max_dist = dist;
226 
227  for(unsigned int oi = 1; oi < nd.size(); ++oi){ // start from second-highest density
228  dist = max_dist;
229  unsigned int i = rs[oi];
230  // we only need to check up to oi since hits
231  // are ordered by decreasing density
232  // and all points coming BEFORE oi are guaranteed to have higher rho
233  // and the ones AFTER to have lower rho
234  for(unsigned int oj = 0; oj < oi; oj++){
235  unsigned int j = rs[oj];
236  double tmp = distance(nd[i].data, nd[j].data);
237  if(tmp <= dist){ //this "<=" instead of "<" addresses the (rare) case when there are only two hits
238  dist = tmp;
239  nearestHigher = j;
240  }
241  }
242  nd[i].data.delta = dist;
243  nd[i].data.nearestHigher = nearestHigher; //this uses the original unsorted hitlist
244  }
245  return maxdensity;
246 }
int i
Definition: DBlmapReader.cc:9
double distance(const Hexel &pt1, const Hexel &pt2)
int j
Definition: DBlmapReader.cc:9
std::vector< size_t > sorted_indices(const std::vector< T > &v)
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
double HGCalImagingAlgo::calculateEnergyWithFraction ( const std::vector< KDNode > &  hits,
const std::vector< double > &  fractions 
)
private

Definition at line 406 of file HGCalImagingAlgo.cc.

References i, and mps_fire::result.

Referenced by getClusters(), and shareEnergy().

407  {
408  double result = 0.0;
409  for( unsigned i = 0 ; i < hits.size(); ++i ) {
410  result += fractions[i]*hits[i].data.weight;
411  }
412  return result;
413 }
int i
Definition: DBlmapReader.cc:9
tuple result
Definition: mps_fire.py:84
double HGCalImagingAlgo::calculateLocalDensity ( std::vector< KDNode > &  nd,
KDTree lp 
)
private

Definition at line 181 of file HGCalImagingAlgo.cc.

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

Referenced by makeClusters().

181  {
182  double maxdensity = 0.;
183  for(unsigned int i = 0; i < nd.size(); ++i){
184  KDTreeBox search_box(nd[i].dims[0]-delta_c,nd[i].dims[0]+delta_c,
185  nd[i].dims[1]-delta_c,nd[i].dims[1]+delta_c);
186  std::vector<KDNode> found;
187  lp.search(search_box,found);
188  for(unsigned int j = 0; j < found.size(); j++){
189  if(distance(nd[i].data,found[j].data) < delta_c){
190  nd[i].data.rho += found[j].data.weight;
191  if(nd[i].data.rho > maxdensity) maxdensity = nd[i].data.rho;
192  }
193  }
194  }
195  return maxdensity;
196 }
int i
Definition: DBlmapReader.cc:9
double distance(const Hexel &pt1, const Hexel &pt2)
int j
Definition: DBlmapReader.cc:9
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
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 getClusters().

155  {
156  float total_weight = 0.;
157  float x = 0.;
158  float y = 0.;
159  float z = 0.;
160  for (unsigned int i = 0; i < v.size(); i++){
161  if(!v[i].data.isHalo){
162  total_weight += v[i].data.weight;
163  x += v[i].data.x*v[i].data.weight;
164  y += v[i].data.y*v[i].data.weight;
165  z += v[i].data.z*v[i].data.weight;
166  }
167  }
168 
169  return math::XYZPoint( x/total_weight,
170  y/total_weight,
171  z/total_weight );
172 }
int i
Definition: DBlmapReader.cc:9
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
math::XYZPoint HGCalImagingAlgo::calculatePositionWithFraction ( const std::vector< KDNode > &  hits,
const std::vector< double > &  fractions 
)
private

Definition at line 390 of file HGCalImagingAlgo.cc.

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

Referenced by getClusters(), and shareEnergy().

391  {
392  double norm(0.0), x(0.0), y(0.0), z(0.0);
393  for( unsigned i = 0; i < hits.size(); ++i ) {
394  const double weight = fractions[i]*hits[i].data.weight;
395  norm += weight;
396  x += weight*hits[i].data.x;
397  y += weight*hits[i].data.y;
398  z += weight*hits[i].data.z;
399  }
401  double norm_inv = 1.0/norm;
402  result *= norm_inv;
403  return result;
404 }
int i
Definition: DBlmapReader.cc:9
tuple result
Definition: mps_fire.py:84
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
int weight
Definition: histoStyle.py:50
double HGCalImagingAlgo::distance ( const Hexel pt1,
const Hexel pt2 
)
private

Definition at line 174 of file HGCalImagingAlgo.cc.

References mathSSE::sqrt(), HGCalImagingAlgo::Hexel::x, and HGCalImagingAlgo::Hexel::y.

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

174  {
175  const double dx = pt1.x - pt2.x;
176  const double dy = pt1.y - pt2.y;
177  return std::sqrt(dx*dx + dy*dy);
178 }
T sqrt(T t)
Definition: SSEVec.h:18
int HGCalImagingAlgo::findAndAssignClusters ( std::vector< KDNode > &  nd,
KDTree lp,
double  maxdensity,
KDTreeBox bounds 
)
private

Definition at line 248 of file HGCalImagingAlgo.cc.

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

Referenced by makeClusters().

248  {
249 
250  //this is called once per layer...
251  //so when filling the cluster temporary vector of Hexels we resize each time by the number
252  //of clusters found. This is always equal to the number of cluster centers...
253 
254  unsigned int clusterIndex = 0;
255 
256  std::vector<size_t> rs = sorted_indices(nd); // indices sorted by decreasing rho
257  std::vector<size_t> ds = sort_by_delta(nd); // sort in decreasing distance to higher
258 
259 
260  for(unsigned int i =0; i < nd.size(); ++i){
261 
262  // std::cout << " delta " << lp[ds[i]].delta << " rho " << lp[ds[i]].rho << std::endl;
263  if(nd[ds[i]].data.delta < delta_c) break; // no more cluster centers to be looked at
264  if(nd[ds[i]].data.rho < maxdensity/kappa /* || lp[ds[i]].rho<0.001*/) continue;
265  //skip this as a potential cluster center because it fails the density cut
266 
267  nd[ds[i]].data.clusterIndex = clusterIndex;
268  if (verbosity < pINFO)
269  {
270  std::cout << "Adding new cluster with index " << clusterIndex+cluster_offset << std::endl;
271  std::cout << "Cluster center is hit " << ds[i] << std::endl;
272  }
273  clusterIndex++;
274  }
275 
276  //at this point clusterIndex is equal to the number of cluster centers - if it is zero we are
277  //done
278  if(clusterIndex==0) return clusterIndex;
279 
280  //assign to clusters, using the nearestHigher set from previous step (always set except
281  // for top density hit that is skipped...
282  for(unsigned int oi =1; oi < nd.size(); ++oi){
283  unsigned int i = rs[oi];
284  int ci = nd[i].data.clusterIndex;
285  if(ci == -1){
286  nd[i].data.clusterIndex = nd[nd[i].data.nearestHigher].data.clusterIndex;
287  }
288  }
289 
290  //make room in the temporary cluster vector for the additional clusterIndex clusters
291  // from this layer
292  if (verbosity < pINFO)
293  {
294  std::cout << "resizing cluster vector by "<< clusterIndex << std::endl;
295  }
296  current_v.resize(cluster_offset+clusterIndex);
297 
298  //assign points closer than dc to other clusters to border region
299  //and find critical border density
300  std::vector<double> rho_b(clusterIndex,0.);
301  lp.clear();
302  lp.build(nd,bounds);
303  //now loop on all hits again :( and check: if there are hits from another cluster within d_c -> flag as border hit
304  for(unsigned int i = 0; i < nd.size(); ++i){
305  int ci = nd[i].data.clusterIndex;
306  bool flag_isolated = true;
307  if(ci != -1){
308  KDTreeBox search_box(nd[i].dims[0]-delta_c,nd[i].dims[0]+delta_c,
309  nd[i].dims[1]-delta_c,nd[i].dims[1]+delta_c);
310  std::vector<KDNode> found;
311  lp.search(search_box,found);
312 
313  for(unsigned int j = 1; j < found.size(); j++){
314  //check if the hit is not within d_c of another cluster
315  if(found[j].data.clusterIndex!=-1){
316  float dist = distance(found[j].data,nd[i].data);
317  if(dist < delta_c && found[j].data.clusterIndex!=ci){
318  //in which case we assign it to the border
319  nd[i].data.isBorder = true;
320  break;
321  }
322  //because we are using two different containers, we have to make sure that we don't unflag the
323  // hit when it finds *itself* closer than delta_c
324  if(dist < delta_c && dist != 0. && found[j].data.clusterIndex==ci){
325  //this is not an isolated hit
326  flag_isolated = false;
327  }
328  }
329  }
330  if(flag_isolated) nd[i].data.isBorder = true; //the hit is more than delta_c from any of its brethren
331  }
332  //check if this border hit has density larger than the current rho_b and update
333  if(nd[i].data.isBorder && rho_b[ci] < nd[i].data.rho)
334  rho_b[ci] = nd[i].data.rho;
335  }
336 
337  //flag points in cluster with density < rho_b as halo points, then fill the cluster vector
338  for(unsigned int i = 0; i < nd.size(); ++i){
339  int ci = nd[i].data.clusterIndex;
340  if(ci!=-1 && nd[i].data.rho < rho_b[ci])
341  nd[i].data.isHalo = true;
342  if(nd[i].data.clusterIndex!=-1){
343  current_v[ci+cluster_offset].push_back(nd[i]);
344  if (verbosity < pINFO)
345  {
346  std::cout << "Pushing hit " << i << " into cluster with index " << ci+cluster_offset << std::endl;
347  std::cout << "Size now " << current_v[ci+cluster_offset].size() << std::endl;
348  }
349  }
350  }
351 
352  //prepare the offset for the next layer if there is one
353  if (verbosity < pINFO)
354  {
355  std::cout << "moving cluster offset by " << clusterIndex << std::endl;
356  }
357  cluster_offset += clusterIndex;
358  return clusterIndex;
359 }
int i
Definition: DBlmapReader.cc:9
std::vector< std::vector< KDNode > > current_v
double distance(const Hexel &pt1, const Hexel &pt2)
std::vector< size_t > sort_by_delta(const std::vector< KDNode > &v)
unsigned int cluster_offset
VerbosityLevel verbosity
int j
Definition: DBlmapReader.cc:9
std::vector< size_t > sorted_indices(const std::vector< T > &v)
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
tuple cout
Definition: gather_cfg.py:145
std::vector< unsigned > HGCalImagingAlgo::findLocalMaximaInCluster ( const std::vector< KDNode > &  cluster)
private

Definition at line 362 of file HGCalImagingAlgo.cc.

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

Referenced by getClusters().

362  {
363  std::vector<unsigned> result;
364  std::vector<bool> seed(cluster.size(),true);
365 
366  for( unsigned i = 0; i < cluster.size(); ++i ) {
367  for( unsigned j = 0; j < cluster.size(); ++j ) {
368  if( distance(cluster[i].data,cluster[j].data) < delta_c && i != j) {
369  //std::cout << "hit-to-hit distance = " << distance(cluster[i],cluster[j]) << std::endl;
370  if( cluster[i].data.weight < cluster[j].data.weight ) {
371  seed[i] = false;
372  break;
373  }
374  }
375  }
376  }
377 
378  for( unsigned i = 0 ; i < cluster.size(); ++i ) {
379  if( seed[i] && cluster[i].data.weight > 5e-4) {
380  //std::cout << "seed at " << i << " with energy " << cluster[i].weight << std::endl;
381  result.push_back(i);
382  }
383  }
384 
385  //std::cout << "Found " << result.size() << " sub-clusters in input cluster of length: " << cluster.size() << std::endl;
386 
387  return result;
388 }
int i
Definition: DBlmapReader.cc:9
double distance(const Hexel &pt1, const Hexel &pt2)
tuple result
Definition: mps_fire.py:84
int j
Definition: DBlmapReader.cc:9
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
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, alignCSCRings::e, relval_parameters_module::energy, findLocalMaximaInCluster(), HLT_25ns10e33_v2_cff::fraction, i, pINFO, position, shareEnergy(), and verbosity.

79  {
80 
82  std::vector< std::pair<DetId, float> > thisCluster;
83  for (unsigned int i = 0; i < current_v.size(); ++i){
84  double energy = 0;
86 
87  if( doSharing ) {
88 
89  std::vector<unsigned> seeds = findLocalMaximaInCluster(current_v[i]);
90  //std::cout << " sharing found " << seeds.size() << " sub-cluster seeds in cluster " << i << std::endl;
91 
92  std::vector<std::vector<double> > fractions;
93  // first pass can have noise it in
94  shareEnergy(current_v[i],seeds,fractions);
95 
96  // reset and run second pass after vetoing seeds
97  // that result in trivial clusters (less than 2 effective cells)
98 
99 
100  for( unsigned isub = 0; isub < fractions.size(); ++isub ) {
101  double effective_hits = 0.0;
102  double energy = calculateEnergyWithFraction(current_v[i],fractions[isub]);
103  Point position = calculatePositionWithFraction(current_v[i],fractions[isub]);
104 
105  //std::cout << "Fractions*Energies: ";
106  for( unsigned ihit = 0; ihit < fractions[isub].size(); ++ihit ) {
107  const double fraction = fractions[isub][ihit];
108  if( fraction > 1e-7 ) {
109  //std::cout << fraction << "*" << current_v[i][ihit].weight << " ";
110  effective_hits += fraction;
111  thisCluster.emplace_back(current_v[i][ihit].data.detid,fraction);
112  }
113  }
114  //std::cout << std::endl;
115 
116  if (verbosity < pINFO)
117  {
118  std::cout << "\t******** NEW CLUSTER (SHARING) ********" << std::endl;
119  std::cout << "\tEff. No. of cells = " << effective_hits << std::endl;
120  std::cout << "\t Energy = " << energy << std::endl;
121  std::cout << "\t Phi = " << position.phi() << std::endl;
122  std::cout << "\t Eta = " << position.eta() << std::endl;
123  std::cout << "\t*****************************" << std::endl;
124  }
125  clusters_v.push_back(reco::BasicCluster(energy, position, caloID, thisCluster,
126  algoId));
127  thisCluster.clear();
128  }
129  }else{
130  position = calculatePosition(current_v[i]);
131  std::vector< KDNode >::iterator it;
132  for (it = current_v[i].begin(); it != current_v[i].end(); it++)
133  {
134  energy += (*it).data.isHalo ? 0. : (*it).data.weight;
135  thisCluster.emplace_back(std::pair<DetId, float>((*it).data.detid,((*it).data.isHalo?0.:1.)));
136  };
137  if (verbosity < pINFO)
138  {
139  std::cout << "******** NEW CLUSTER (HGCIA) ********" << std::endl;
140  std::cout << "Index " << i << std::endl;
141  std::cout << "No. of cells = " << current_v[i].size() << std::endl;
142  std::cout << " Energy = " << energy << std::endl;
143  std::cout << " Phi = " << position.phi() << std::endl;
144  std::cout << " Eta = " << position.eta() << std::endl;
145  std::cout << "*****************************" << std::endl;
146  }
147  clusters_v.push_back(reco::BasicCluster(energy, position, caloID, thisCluster,
148  algoId));
149  thisCluster.clear();
150  }
151  }
152  return clusters_v;
153 }
double calculateEnergyWithFraction(const std::vector< KDNode > &, const std::vector< double > &)
int i
Definition: DBlmapReader.cc:9
std::vector< std::vector< KDNode > > current_v
math::XYZPoint calculatePositionWithFraction(const std::vector< KDNode > &, const std::vector< double > &)
reco::CaloCluster::AlgoId algoId
void shareEnergy(const std::vector< KDNode > &, const std::vector< unsigned > &, std::vector< std::vector< double > > &)
math::XYZPoint Point
VerbosityLevel verbosity
std::vector< unsigned > findLocalMaximaInCluster(const std::vector< KDNode > &)
math::XYZPoint calculatePosition(std::vector< KDNode > &)
#define begin
Definition: vmac.h:30
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
static int position[264][3]
Definition: ReadPGInfo.cc:509
tuple cout
Definition: gather_cfg.py:145
std::vector< reco::BasicCluster > clusters_v
void HGCalImagingAlgo::getEventSetup ( const edm::EventSetup es)
inline

Definition at line 94 of file HGCalImagingAlgo.h.

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

94 { rhtools_.getEventSetup(es); }
hgcal::RecHitTools rhtools_
void getEventSetup(const edm::EventSetup &)
Definition: RecHitTools.cc:66
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, hgcal::RecHitTools::isHalfCell(), kappa, bookConverter::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 ecaldqm::zside().

17 {
18  //used for speedy search
19  std::vector<std::vector<KDNode> >points(2*(maxlayer+1));
20  std::vector<KDTree> hit_kdtree(2*(maxlayer+1));
21 
22  std::vector<std::array<float,2> > minpos(2*(maxlayer+1),{ {0.0f,0.0f} }),maxpos(2*(maxlayer+1),{ {0.0f,0.0f} });
23 
24  // std::vector<std::vector<Hexel> > points(2*(maxlayer+1)); //a vector of vectors of hexels, one for each layer
25  //@@EM todo: the number of layers should be obtained programmatically - the range is 1-n instead of 0-n-1...
26 
27 
28  if (verbosity < pINFO)
29  {
30  std::cout << "-------------------------------------------------------------" << std::endl;
31  std::cout << "HGC Imaging algorithm invoked for " << std::endl;
32  std::cout << "delta_c " << delta_c << " kappa " << kappa;
33  // if( doSharing ) std::cout << " showerSigma " << std::sqrt(sigma2);
34  std::cout << std::endl;
35  }
36 
37  //loop over all hits and create the Hexel structure, skip energies below ecut
38  for (unsigned int i=0;i<hits.size();++i) {
39  const HGCRecHit& hgrh = hits[i];
40  if(hgrh.energy() < ecut) continue;
41  DetId detid = hgrh.detid();
42 
43  int layer = rhtools_.getLayerWithOffset(detid)+int(HGCalDetId(detid).zside()>0)*(maxlayer+1);
44 
45  // determine whether this is a half-hexagon
46  bool isHalf = rhtools_.isHalfCell(detid);
47  const GlobalPoint position( std::move( rhtools_.getPosition( detid ) ) );
48  //here's were the KDNode is passed its dims arguments - note that these are *copied* from the Hexel
49  points[layer].emplace_back(Hexel(hgrh,detid,isHalf,&rhtools_),position.x(),position.y());
50  if(points[layer].size()==0){
51  minpos[layer][0] = position.x(); minpos[layer][1] = position.y();
52  maxpos[layer][0] = position.x(); maxpos[layer][1] = position.y();
53  }else{
54  minpos[layer][0] = std::min((float)position.x(),minpos[layer][0]);
55  minpos[layer][1] = std::min((float)position.y(),minpos[layer][1]);
56  maxpos[layer][0] = std::max((float)position.x(),maxpos[layer][0]);
57  maxpos[layer][1] = std::max((float)position.y(),maxpos[layer][1]);
58 
59  }
60  }
61 
62  //assign all hits in each layer to a cluster core or halo
63  for (unsigned int i = 0; i <= 2*maxlayer+1; ++i) {
64  KDTreeBox bounds(minpos[i][0],maxpos[i][0],
65  minpos[i][1],maxpos[i][1]);
66 
67  hit_kdtree[i].build(points[i],bounds);
68 
69  double maxdensity = calculateLocalDensity(points[i],hit_kdtree[i]);
70  // std::cout << "layer " << i << " max density " << maxdensity
71  // << " total hits " << points[i].size() << std::endl;
72  calculateDistanceToHigher(points[i],hit_kdtree[i]);
73  findAndAssignClusters(points[i],hit_kdtree[i],maxdensity,bounds);
74  // std::cout << "found " << nclusters << " clusters" << std::endl;
75  }
76  //make the cluster vector
77 }
double calculateDistanceToHigher(std::vector< KDNode > &, KDTree &)
bool isHalfCell(const DetId &) const
Definition: RecHitTools.cc:158
hgcal::RecHitTools rhtools_
int i
Definition: DBlmapReader.cc:9
static const unsigned int maxlayer
const DetId & detid() const
Definition: CaloRecHit.h:21
double calculateLocalDensity(std::vector< KDNode > &, KDTree &)
int zside(DetId const &)
float energy() const
Definition: CaloRecHit.h:17
def move
Definition: eostools.py:510
std::vector< std::vector< Hexel > > points
VerbosityLevel verbosity
T min(T a, T b)
Definition: MathUtil.h:58
Definition: DetId.h:18
unsigned int getLayerWithOffset(const DetId &) const
Definition: RecHitTools.cc:126
GlobalPoint getPosition(const DetId &id) const
Definition: RecHitTools.cc:73
size_type size() const
static int position[264][3]
Definition: ReadPGInfo.cc:509
tuple cout
Definition: gather_cfg.py:145
int findAndAssignClusters(std::vector< KDNode > &, KDTree &, double, KDTreeBox &)
tuple size
Write out results.
void HGCalImagingAlgo::reset ( void  )
inline

Definition at line 96 of file HGCalImagingAlgo.h.

References cluster_offset, clusters_v, and current_v.

96  {
97  current_v.clear();
98  clusters_v.clear();
99  cluster_offset = 0;
100  }
std::vector< std::vector< KDNode > > current_v
unsigned int cluster_offset
std::vector< reco::BasicCluster > clusters_v
void HGCalImagingAlgo::setVerbosity ( VerbosityLevel  the_verbosity)
inline

Definition at line 83 of file HGCalImagingAlgo.h.

References verbosity.

84  {
85  verbosity = the_verbosity;
86  }
VerbosityLevel verbosity
void HGCalImagingAlgo::shareEnergy ( const std::vector< KDNode > &  incluster,
const std::vector< unsigned > &  seeds,
std::vector< std::vector< double > > &  outclusters 
)
private

Definition at line 415 of file HGCalImagingAlgo.cc.

References calculateEnergyWithFraction(), calculatePositionWithFraction(), data, mps_update::diff, alignCSCRings::e, create_public_lumi_plots::exp, cropTnPTrees::frac, HLT_25ns10e33_v2_cff::fraction, i, j, bookConverter::max, eostools::move(), perp2(), funct::pow(), sigma2, mathSSE::sqrt(), HGCalImagingAlgo::Hexel::x, HGCalImagingAlgo::Hexel::y, and HGCalImagingAlgo::Hexel::z.

Referenced by getClusters().

417  {
418  std::vector<bool> isaseed(incluster.size(),false);
419  outclusters.clear();
420  outclusters.resize(seeds.size());
421  std::vector<Point> centroids(seeds.size());
422  std::vector<double> energies(seeds.size());
423 
424  if( seeds.size() == 1 ) { // short circuit the case of a lone cluster
425  outclusters[0].clear();
426  outclusters[0].resize(incluster.size(),1.0);
427  return;
428  }
429 
430  //std::cout << "saving seeds" << std::endl;
431 
432  // create quick seed lookup
433  for( unsigned i = 0; i < seeds.size(); ++i ) {
434  isaseed[seeds[i]] = true;
435  }
436 
437  //std::cout << "seeds saved" << std::endl;
438 
439  // initialize clusters to be shared
440  // centroids start off at seed positions
441  // seeds always have fraction 1.0, to stabilize fit
442  //std::cout << "initializing fit" << std::endl;
443  for( unsigned i = 0; i < seeds.size(); ++i ) {
444  outclusters[i].resize(incluster.size(),0.0);
445  for( unsigned j = 0; j < incluster.size(); ++j ) {
446  if( j == seeds[i] ) {
447  outclusters[i][j] = 1.0;
448  centroids[i] = math::XYZPoint(incluster[j].data.x,incluster[j].data.y,incluster[j].data.z);
449  energies[i] = incluster[j].data.weight;
450  }
451  }
452  }
453 
454  //return;
455 
456  //std::cout << "fit initialized" << std::endl;
457 
458  // run the fit while we are less than max iterations, and clusters are still moving
459  const double minFracTot = 1e-20;
460  unsigned iter = 0;
461  const unsigned iterMax = 50;
463  const double stoppingTolerance = 1e-8;
464  const double toleranceScaling = std::pow(std::max(1.0,seeds.size()-1.0),2.0);
465  std::vector<Point> prevCentroids;
466  std::vector<double> frac(seeds.size()), dist2(seeds.size());
467  while( iter++ < iterMax && diff > stoppingTolerance*toleranceScaling ) {
468  for( unsigned i = 0; i < incluster.size(); ++i ) {
469  const Hexel& ihit = incluster[i].data;
470  double fraction(0.0), fracTot(0.0), d2(0.0);
471  for( unsigned j = 0; j < seeds.size(); ++j ) {
472  fraction = 0.0;
473  d2 = ( std::pow(ihit.x - centroids[j].x(),2.0) +
474  std::pow(ihit.y - centroids[j].y(),2.0) +
475  std::pow(ihit.z - centroids[j].z(),2.0) )/sigma2;
476  dist2[j] = d2;
477  // now we set the fractions up based on hit type
478  if( i == seeds[j] ) { // this cluster's seed
479  fraction = 1.0;
480  } else if( isaseed[i] ) {
481  fraction = 0.0;
482  } else {
483  fraction = energies[j]*std::exp( -0.5*d2 );
484  }
485  fracTot += fraction;
486  frac[j] = fraction;
487  }
488  // now that we have calculated all fractions for all hits
489  // assign the new fractions
490  for( unsigned j = 0; j < seeds.size(); ++j ) {
491  if( fracTot > minFracTot ||
492  ( i == seeds[j] && fracTot > 0.0 ) ) {
493  outclusters[j][i] = frac[j]/fracTot;
494  } else {
495  outclusters[j][i] = 0.0;
496  }
497  }
498  }
499 
500  // save previous centroids
501  prevCentroids = std::move(centroids);
502  // finally update the position of the centroids from the last iteration
503  centroids.resize(seeds.size());
504  double diff2 = 0.0;
505  for( unsigned i = 0; i < seeds.size(); ++i ) {
506  centroids[i] = calculatePositionWithFraction(incluster,outclusters[i]);
507  energies[i] = calculateEnergyWithFraction(incluster,outclusters[i]);
508  // calculate convergence parameters
509  const double delta2 = (prevCentroids[i]-centroids[i]).perp2();
510  if( delta2 > diff2 ) diff2 = delta2;
511  }
512  //update convergance parameter outside loop
513  diff = std::sqrt(diff2);
514  //std::cout << " iteration = " << iter << " diff = " << diff << std::endl;
515  }
516 }
double calculateEnergyWithFraction(const std::vector< KDNode > &, const std::vector< double > &)
int i
Definition: DBlmapReader.cc:9
math::XYZPoint calculatePositionWithFraction(const std::vector< KDNode > &, const std::vector< double > &)
list diff
Definition: mps_update.py:85
T sqrt(T t)
Definition: SSEVec.h:18
def move
Definition: eostools.py:510
int j
Definition: DBlmapReader.cc:9
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
T perp2() const
Squared magnitude of transverse component.
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
std::vector<size_t> HGCalImagingAlgo::sort_by_delta ( const std::vector< KDNode > &  v)
inlineprivate

Definition at line 185 of file HGCalImagingAlgo.h.

References i, and findQualityFiles::v.

Referenced by findAndAssignClusters().

185  {
186  std::vector<size_t> idx(v.size());
187  for (size_t i = 0; i != idx.size(); ++i) idx[i] = i;
188  sort(idx.begin(), idx.end(),
189  [&v](size_t i1, size_t i2) {return v[i1].data.delta > v[i2].data.delta;});
190  return idx;
191  }
int i
Definition: DBlmapReader.cc:9

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
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().