NNClusters.h

/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
#ifndef NNClusters_h
#define NNClusters_h 1

#include <list>
#include <vector>

#include "lcio.h"
#include "EVENT/LCObject.h"
#include "EVENT/LCCollection.h"
#include "IMPL/ClusterImpl.h"

#include "ClusterShapes.h"
#include "CLHEP/Vector/ThreeVector.h"

// fix for transition from CLHEP 1.8 to 1.9
namespace CLHEP{}
using namespace CLHEP ;

// forward declarations:

template <class U>
class GenericCluster ;


template <class In, class Out, class Pred > 
void cluster( In first, In last, Out result, Pred* pred ) {

  typedef typename In::value_type GenericHitPtr ;
  typedef typename Pred::hit_type HitType ;

  typedef std::vector< GenericCluster<HitType >* >  ClusterList ;

  ClusterList tmp ; 
  tmp.reserve( 1024 ) ;
  
//   int i(0),j(0) ;

  while( first != last ) {

//     j=i+1 ;

    for( In other = first+1 ; other != last  ; other ++ ) {
      
//       std::cout << "  in nested loop " << i << " , " << j << std::endl ;

      if( pred->mergeHits( (*first) , (*other) ) ) {
        
        if( (*first)->second == 0 && (*other)->second == 0 ) {  // no cluster exists
          
          GenericCluster<HitType >* cl = new GenericCluster<HitType >( (*first) ) ;

          cl->addHit( (*other) ) ;
          
          tmp.push_back( cl ) ;

        }
        else if( (*first)->second != 0 && (*other)->second != 0 ) { // two clusters
          
//        if(  (*first)->second == (*other)->second ) 
//          std::cout << " Merging identical clusters !? " << std::endl ;

          if(  (*first)->second != (*other)->second )  // this is a bug fix for old gcc 3.2 compiler !
            (*first)->second->mergeClusters( (*other)->second ) ;
          
        } else {  // one cluster exists
          
          if( (*first)->second != 0 ) {
            
            (*first)->second->addHit( (*other)  ) ;
            
          } else {                           
            
            (*other)->second->addHit( (*first)  ) ;
          }
        }
        
      } // dCut 
//       ++j ;
    }
//     ++i ;
    ++first ;
  }

  // remove empty clusters 
  //   std::remove_copy_if( tmp.begin() , tmp.end() , result , &empty_list< GenericCluster<HitType > > ) ;

  for( typename ClusterList::iterator i = tmp.begin(); i !=  tmp.end() ; i++ ){

    if( (*i)->size() > 0 ) {

      result++ = *i ;
    } 
    else {  delete *i ; }
  }
}

template <class T>
class GenericHit : public  std::pair< T*, GenericCluster<T>* >{

  typedef T value_type ;
  typedef std::pair< T*, GenericCluster<T>* > Pair ;

public:

  GenericHit(T* hit, int index0 = 0 ) : Index0( index0 )   {  
    Pair::first = hit ; 
    Pair::second = 0 ;
  }

  GenericHit(T* hit ,  GenericCluster<T>* cl , int index0 = 0) : Index0( index0 ) {
    Pair::first =  hit ;
    Pair::second = cl ;
  }
  
  int Index0 ;
  
protected:
  GenericHit() ;
} ;


template <class T >
class GenericCluster : public std::list< GenericHit<T> * > {

public :

  GenericCluster( GenericHit<T>* hit) {
    addHit( hit ) ;
  }
  
  void addHit( GenericHit<T>* hit ) {
    
    hit->second = this ;
    this->push_back( hit ) ;

  }

  void mergeClusters( GenericCluster<T>* cl ) {
    
    for( typename GenericCluster<T>::iterator it = cl->begin() ; it !=  cl->end() ; it++ ){
      (*it)->second = this  ;
    }
    this->merge( *cl ) ;
  }
} ;


template <class T> 
class GenericHitVec : public std::vector< GenericHit<T>* > {
  typedef std::vector< GenericHit<T>* > Vector ;
public:
  ~GenericHitVec() {
    for( typename GenericHitVec::iterator i = Vector::begin() ; i != Vector::end() ; i++) delete *i ;
  }
};


template <class T, class Pred> 
void addToGenericHitVec(GenericHitVec<T>& v, LCCollection* col, Pred pred ){

  for( int i=0 ; i < col->getNumberOfElements() ; i++ ){

    T* hit = dynamic_cast<T*>( col->getElementAt( i) ) ;

    if( pred( hit ) ){

      v.push_back( new GenericHit<T>( hit ) ) ;
    }
  }
} 

template <class T, class Pred, class Order> 
void addToGenericHitVec(GenericHitVec<T>& v, LCCollection* col, Pred pred , Order order ){

  for( int i=0 ; i < col->getNumberOfElements() ; i++ ){

    T* hit = dynamic_cast<T*>( col->getElementAt( i) ) ;

    if( pred( hit ) ){

      v.push_back( new GenericHit<T>( hit , order(hit) ) ) ;
    }
  }
} 

template <class T> 
class GenericClusterVec : public std::list< GenericCluster<T>* > {
  typedef std::list< GenericCluster<T>* > List ;
public:
  ~GenericClusterVec() {
    for( typename GenericClusterVec::iterator i = List::begin() ; i != List::end() ; i++) delete *i ;
  }
};


template <class HitClass, typename PosType > 
class NNDistance{
public:

  typedef HitClass hit_type ;

  NNDistance(float dCut) : _dCutSquared( dCut*dCut ) , _dCut(dCut)  {} 


  inline bool mergeHits( GenericHit<HitClass>* h0, GenericHit<HitClass>* h1){
    
    if( std::abs( h0->Index0 - h1->Index0 ) > 1 ) return false ;

    const PosType* pos0 =  h0->first->getPosition() ;
    const PosType* pos1 =  h1->first->getPosition() ;
    
    return 
      ( pos0[0] - pos1[0] ) * ( pos0[0] - pos1[0] ) +
      ( pos0[1] - pos1[1] ) * ( pos0[1] - pos1[1] ) +
      ( pos0[2] - pos1[2] ) * ( pos0[2] - pos1[2] )   
      < _dCutSquared ;
  }
  
protected:
  NNDistance() ;
  float _dCutSquared ;
  float _dCut ;
} ;


template <class T>
class EnergyCut{
public:
  EnergyCut( double eCut ) : _eCut( eCut ) {}  

  inline bool operator() (T* hit) {  return hit->getEnergy() > _eCut ; }

protected:

  EnergyCut() {} ;
  double _eCut ;
} ;


template <class T, int N>
class ZIndex{
public:
  ZIndex( float zmin , float zmax ) : _zmin( zmin ), _zmax( zmax ) {}  

  inline int operator() (T* hit) {  

    return (int) std::floor( ( hit->getPosition()[2] - _zmin ) / ( _zmax - _zmin ) * N ) ; 
  }

protected:

  ZIndex() {} ;
  float _zmin ;
  float _zmax ;
} ;


template <class T>
struct LCIOCluster{

  inline lcio::Cluster* operator() (GenericCluster<T>* c) {  
    
    ClusterImpl* clu = new ClusterImpl ;
    
    unsigned n = c->size() ;
    unsigned i=0 ;

    float a[n], x[n], y[n], z[n] ;

    for( typename GenericCluster<T>::iterator hi = c->begin(); hi != c->end() ; hi++) {
      
      T* hit = (*hi)->first ;  

      a[i] = hit->getEnergy() ;
      x[i] = hit->getPosition()[0] ;
      y[i] = hit->getPosition()[1] ;
      z[i] = hit->getPosition()[2] ;

      clu->addHit( hit , a[i] ) ;
      
      ++i ;
    }
    
    ClusterShapes cs( n,a,x,y,z) ;

    clu->setEnergy( cs.getTotalAmplitude()  ) ;
    clu->setPosition( cs.getCenterOfGravity() ) ;

    // direction of cluster's PCA
    float* d = cs.getEigenVecInertia() ;

    Hep3Vector v( d[0], d[1], d[2] ) ;
    
    clu->setITheta( v.theta() )  ;
    clu->setIPhi(   v.phi() ) ;
  
    std::vector<float> param(5) ;

//     float* ev = cs.getEigenValInertia() ;
//     param[0] = ev[0] ;
//     param[1] = ev[1] ;
//     param[2] = ev[2] ;

    param[0] = cs.getElipsoid_r1() ;
    param[1] = cs.getElipsoid_r2() ;
    param[2] = cs.getElipsoid_r3() ;
    param[3] = cs.getElipsoid_vol() ;
    param[4] = cs.getWidth() ;

    clu->setShape( param ) ;

    return clu ;

  }

} ;



#endif