MarlinUtil  1.12.1
NNClusters.h
1 /* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
2 #ifndef NNClusters_h
3 #define NNClusters_h 1
4 
11 #include <list>
12 #include <vector>
13 
14 #include "lcio.h"
15 #include "EVENT/LCObject.h"
16 #include "EVENT/LCCollection.h"
17 #include "IMPL/ClusterImpl.h"
18 
19 #include "ClusterShapes.h"
20 #include "CLHEP/Vector/ThreeVector.h"
21 
22 // fix for transition from CLHEP 1.8 to 1.9
23 namespace CLHEP{}
24 using namespace CLHEP ;
25 
26 // forward declarations:
27 
28 template <class U>
30 
31 
47 template <class In, class Out, class Pred >
48 void cluster( In first, In last, Out result, Pred* pred ) {
49 
50  typedef typename In::value_type GenericHitPtr ;
51  typedef typename Pred::hit_type HitType ;
52 
53  typedef std::vector< GenericCluster<HitType >* > ClusterList ;
54 
55  ClusterList tmp ;
56  tmp.reserve( 1024 ) ;
57 
58 // int i(0),j(0) ;
59 
60  while( first != last ) {
61 
62 // j=i+1 ;
63 
64  for( In other = first+1 ; other != last ; other ++ ) {
65 
66 // std::cout << " in nested loop " << i << " , " << j << std::endl ;
67 
68  if( pred->mergeHits( (*first) , (*other) ) ) {
69 
70  if( (*first)->second == 0 && (*other)->second == 0 ) { // no cluster exists
71 
73 
74  cl->addHit( (*other) ) ;
75 
76  tmp.push_back( cl ) ;
77 
78  }
79  else if( (*first)->second != 0 && (*other)->second != 0 ) { // two clusters
80 
81 // if( (*first)->second == (*other)->second )
82 // std::cout << " Merging identical clusters !? " << std::endl ;
83 
84  if( (*first)->second != (*other)->second ) // this is a bug fix for old gcc 3.2 compiler !
85  (*first)->second->mergeClusters( (*other)->second ) ;
86 
87  } else { // one cluster exists
88 
89  if( (*first)->second != 0 ) {
90 
91  (*first)->second->addHit( (*other) ) ;
92 
93  } else {
94 
95  (*other)->second->addHit( (*first) ) ;
96  }
97  }
98 
99  } // dCut
100 // ++j ;
101  }
102 // ++i ;
103  ++first ;
104  }
105 
106  // remove empty clusters
107  // std::remove_copy_if( tmp.begin() , tmp.end() , result , &empty_list< GenericCluster<HitType > > ) ;
108 
109  for( typename ClusterList::iterator i = tmp.begin(); i != tmp.end() ; i++ ){
110 
111  if( (*i)->size() > 0 ) {
112 
113  result++ = *i ;
114  }
115  else { delete *i ; }
116  }
117 }
118 
127 template <class T>
128 class GenericHit : public std::pair< T*, GenericCluster<T>* >{
129 
130  typedef T value_type ;
131  typedef std::pair< T*, GenericCluster<T>* > Pair ;
132 
133 public:
134 
138  GenericHit(T* hit, int index0 = 0 ) : Index0( index0 ) {
139  Pair::first = hit ;
140  Pair::second = 0 ;
141  }
142 
145  GenericHit(T* hit , GenericCluster<T>* cl , int index0 = 0) : Index0( index0 ) {
146  Pair::first = hit ;
147  Pair::second = cl ;
148  }
149 
153  int Index0 ;
154 
155 protected:
157  GenericHit() ;
158 } ;
159 
160 
168 template <class T >
169 class GenericCluster : public std::list< GenericHit<T> * > {
170 
171 public :
172 
175  addHit( hit ) ;
176  }
177 
179  void addHit( GenericHit<T>* hit ) {
180 
181  hit->second = this ;
182  this->push_back( hit ) ;
183 
184  }
185 
188 
189  for( typename GenericCluster<T>::iterator it = cl->begin() ; it != cl->end() ; it++ ){
190  (*it)->second = this ;
191  }
192  this->merge( *cl ) ;
193  }
194 } ;
195 
196 
200 template <class T>
201 class GenericHitVec : public std::vector< GenericHit<T>* > {
202  typedef std::vector< GenericHit<T>* > Vector ;
203 public:
204  ~GenericHitVec() {
205  for( typename GenericHitVec::iterator i = Vector::begin() ; i != Vector::end() ; i++) delete *i ;
206  }
207 };
208 
209 
214 template <class T, class Pred>
215 void addToGenericHitVec(GenericHitVec<T>& v, LCCollection* col, Pred pred ){
216 
217  for( int i=0 ; i < col->getNumberOfElements() ; i++ ){
218 
219  T* hit = dynamic_cast<T*>( col->getElementAt( i) ) ;
220 
221  if( pred( hit ) ){
222 
223  v.push_back( new GenericHit<T>( hit ) ) ;
224  }
225  }
226 }
227 
232 template <class T, class Pred, class Order>
233 void addToGenericHitVec(GenericHitVec<T>& v, LCCollection* col, Pred pred , Order order ){
234 
235  for( int i=0 ; i < col->getNumberOfElements() ; i++ ){
236 
237  T* hit = dynamic_cast<T*>( col->getElementAt( i) ) ;
238 
239  if( pred( hit ) ){
240 
241  v.push_back( new GenericHit<T>( hit , order(hit) ) ) ;
242  }
243  }
244 }
245 
248 template <class T>
249 class GenericClusterVec : public std::list< GenericCluster<T>* > {
250  typedef std::list< GenericCluster<T>* > List ;
251 public:
252  ~GenericClusterVec() {
253  for( typename GenericClusterVec::iterator i = List::begin() ; i != List::end() ; i++) delete *i ;
254  }
255 };
256 
257 
262 template <class HitClass, typename PosType >
264 public:
265 
268  typedef HitClass hit_type ;
269 
271  NNDistance(float dCut) : _dCutSquared( dCut*dCut ) , _dCut(dCut) {}
272 
273 
276 
277  if( std::abs( h0->Index0 - h1->Index0 ) > 1 ) return false ;
278 
279  const PosType* pos0 = h0->first->getPosition() ;
280  const PosType* pos1 = h1->first->getPosition() ;
281 
282  return
283  ( pos0[0] - pos1[0] ) * ( pos0[0] - pos1[0] ) +
284  ( pos0[1] - pos1[1] ) * ( pos0[1] - pos1[1] ) +
285  ( pos0[2] - pos1[2] ) * ( pos0[2] - pos1[2] )
286  < _dCutSquared ;
287  }
288 
289 protected:
290  NNDistance() ;
291  float _dCutSquared ;
292  float _dCut ;
293 } ;
294 
295 
299 template <class T>
300 class EnergyCut{
301 public:
302  EnergyCut( double eCut ) : _eCut( eCut ) {}
303 
304  inline bool operator() (T* hit) { return hit->getEnergy() > _eCut ; }
305 
306 protected:
307 
308  EnergyCut() {} ;
309  double _eCut ;
310 } ;
311 
312 
317 template <class T, int N>
318 class ZIndex{
319 public:
321  ZIndex( float zmin , float zmax ) : _zmin( zmin ), _zmax( zmax ) {}
322 
323  inline int operator() (T* hit) {
324 
325  return (int) std::floor( ( hit->getPosition()[2] - _zmin ) / ( _zmax - _zmin ) * N ) ;
326  }
327 
328 protected:
329 
330  ZIndex() {} ;
331  float _zmin ;
332  float _zmax ;
333 } ;
334 
335 
339 template <class T>
340 struct LCIOCluster{
341 
342  inline lcio::Cluster* operator() (GenericCluster<T>* c) {
343 
344  ClusterImpl* clu = new ClusterImpl ;
345 
346  unsigned n = c->size() ;
347  unsigned i=0 ;
348 
349  float a[n], x[n], y[n], z[n] ;
350 
351  for( typename GenericCluster<T>::iterator hi = c->begin(); hi != c->end() ; hi++) {
352 
353  T* hit = (*hi)->first ;
354 
355  a[i] = hit->getEnergy() ;
356  x[i] = hit->getPosition()[0] ;
357  y[i] = hit->getPosition()[1] ;
358  z[i] = hit->getPosition()[2] ;
359 
360  clu->addHit( hit , a[i] ) ;
361 
362  ++i ;
363  }
364 
365  ClusterShapes cs( n,a,x,y,z) ;
366 
367  clu->setEnergy( cs.getTotalAmplitude() ) ;
368  clu->setPosition( cs.getCenterOfGravity() ) ;
369 
370  // direction of cluster's PCA
371  float* d = cs.getEigenVecInertia() ;
372 
373  Hep3Vector v( d[0], d[1], d[2] ) ;
374 
375  clu->setITheta( v.theta() ) ;
376  clu->setIPhi( v.phi() ) ;
377 
378  std::vector<float> param(5) ;
379 
380 // float* ev = cs.getEigenValInertia() ;
381 // param[0] = ev[0] ;
382 // param[1] = ev[1] ;
383 // param[2] = ev[2] ;
384 
385  param[0] = cs.getElipsoid_r1() ;
386  param[1] = cs.getElipsoid_r2() ;
387  param[2] = cs.getElipsoid_r3() ;
388  param[3] = cs.getElipsoid_vol() ;
389  param[4] = cs.getWidth() ;
390 
391  clu->setShape( param ) ;
392 
393  return clu ;
394 
395  }
396 
397 } ;
398 
399 
400 
401 #endif
402 
403 
HitClass hit_type
Required typedef for cluster algorithm.
Definition: NNClusters.h:268
NNDistance(float dCut)
C'tor takes merge distance.
Definition: NNClusters.h:271
float getElipsoid_r1()
largest spatial axis length of the ellipsoid derived by the inertia tensor (by their eigenvalues and ...
Definition: ClusterShapes.h:257
Simple predicate class for computing an index from N bins of the z-coordinate of LCObjects that have ...
Definition: NNClusters.h:318
Templated class for generic hit type objects that are to be clustered with an NN-like clustering algo...
Definition: NNClusters.h:128
void addHit(GenericHit< T > *hit)
Add a hit to this cluster - updates the hit's pointer to cluster.
Definition: NNClusters.h:179
float * getCenterOfGravity()
US spelling of getCentreOfGravity.
Definition: ClusterShapes.h:82
float * getEigenVecInertia()
array of the three main axes of inertia (9 entries) starting with the axis corresponding to the small...
Definition: ClusterShapes.cc:696
Templated class for generic clusters of GenericHits that are clustered with an NN-like clustering alg...
Definition: NNClusters.h:29
float getElipsoid_vol()
volume of the ellipsoid
Definition: ClusterShapes.h:276
float getTotalAmplitude()
returns the accumulated amplitude for the whole cluster (just the sum of the energy in all the entrie...
Definition: ClusterShapes.cc:665
Helper class that creates an lcio::Cluster from a generic cluster with hit types that have a getPosit...
Definition: NNClusters.h:340
GenericCluster(GenericHit< T > *hit)
C'tor that takes the first hit.
Definition: NNClusters.h:174
Helper vector of GenericCluster<T> taking care of memory management.
Definition: NNClusters.h:249
float getWidth()
'mean' width of the cluster perpendicular to the main principal axis, defined as: width := sqrt( 1/To...
Definition: ClusterShapes.cc:708
Simple predicate class for applying an energy cut to the objects of type T.
Definition: NNClusters.h:300
GenericHit(T *hit, int index0=0)
Default c'tor takes a pointer to the original hit type object.
Definition: NNClusters.h:138
GenericHit(T *hit, GenericCluster< T > *cl, int index0=0)
C'tor that also takes a pointer to the cluster this hit belongs to - in case seed hits/clusters are u...
Definition: NNClusters.h:145
Simple predicate class for NN clustering.
Definition: NNClusters.h:263
ZIndex(float zmin, float zmax)
C'tor takes zmin and zmax - NB index can be negative and larger than N.
Definition: NNClusters.h:321
int Index0
Index that can be used to code nearest neighbour bins, e.g.
Definition: NNClusters.h:153
Utility class to derive properties of clusters, such as centre of gravity, axes of inertia...
Definition: ClusterShapes.h:25
Helper vector of GenericHit<T> taking care of memory management, i.e.
Definition: NNClusters.h:201
void mergeClusters(GenericCluster< T > *cl)
Merges all hits from the other cluster cl into this cluster.
Definition: NNClusters.h:187
float getElipsoid_r3()
smallest spatial axis length of the ellipsoid derived by the inertia tensor (by their eigenvalues and...
Definition: ClusterShapes.h:271
float getElipsoid_r2()
medium spatial axis length of the ellipsoid derived by the inertia tensor (by their eigenvalues and e...
Definition: ClusterShapes.h:264
bool mergeHits(GenericHit< HitClass > *h0, GenericHit< HitClass > *h1)
Merge condition: true if distance is less than dCut given in the C'tor.
Definition: NNClusters.h:275