THistogram.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: THistogram.cxx,v 1.8 2010/03/31 09:58:59 liucy Exp $
//-----------------------------------------------------------------------------
// Filename : THistogram.h
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class for a histogram used in tracking.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#include <stdio.h>
#include <stdlib.h>
 
#include "TrkReco/TCircle.h"
#include "TrkReco/THistogram.h"
 
THistogram::THistogram( unsigned nBins ) : _nBins( nBins ) {
  _binSize = 2. * M_PI / (float)_nBins;
  _bins    = (unsigned*)malloc( _nBins * sizeof( unsigned ) );
  _masks   = (bool*)malloc( _nBins * sizeof( bool ) );
  _links   = (AList<TMLink>**)malloc( _nBins * sizeof( AList<TMLink>* ) );
  for ( unsigned i = 0; i < _nBins; i++ )
  {
    _bins[i]  = 0;
    _masks[i] = false;
    _links[i] = new AList<TMLink>;
  }
}
 
THistogram::~THistogram() {
  free( _bins );
  free( _masks );
  for ( unsigned i = 0; i < _nBins; i++ ) delete _links[i];
  free( _links );
}
 
void THistogram::dump( const std::string& msg, const std::string& pre ) const {
  std::cout << pre;
  std::cout << "THistogram dump:#bins=" << _nBins << std::endl;
  unsigned nLoops = _nBins / 15 + 1;
  unsigned n0     = 0;
  unsigned n1     = 0;
  for ( unsigned i = 0; i < nLoops; i++ )
  {
    for ( unsigned j = 0; j < 15; j++ )
    {
      if ( n0 == _nBins ) break;
      printf( "%4d", n0 );
      ++n0;
    }
    std::cout << std::endl;
    for ( unsigned j = 0; j < 15; j++ )
    {
      if ( n1 == _nBins ) break;
      if ( !_masks[n1] ) printf( "%4d", _bins[n1] );
      else printf( "-%3d", _bins[n1] );
      ++n1;
    }
    std::cout << std::endl;
  }
 
  if ( msg.find( "detail" ) != std::string::npos )
  {
    for ( unsigned i = 0; i < _nBins; i++ )
    {
      std::cout << "bin " << i << " : ";
      for ( unsigned j = 0; j < _links[i]->length(); j++ )
      { std::cout << ( *_links[i] )[j]->wire()->name() << ","; }
      std::cout << std::endl;
    }
  }
 
  return;
}
 
void THistogram::fillX( const AList<TMLink>& links ) {
  _all            = (AList<TMLink>&)links;
  unsigned nLinks = links.length();
  double   offset = _binSize / 4.;
  for ( unsigned i = 0; i < nLinks; i++ )
  {
    TMLink*           l   = links[i];
    const HepPoint3D& p   = l->position();
    unsigned          pos = (unsigned)floor( ( p.x() + offset ) / _binSize );
 
    //...Why is this needed?...
    pos %= _nBins;
 
    ++_bins[pos];
    _links[pos]->append( l );
  }
}
 
void THistogram::fillY( const AList<TMLink>& links ) {
  _all            = (AList<TMLink>&)links;
  unsigned nLinks = links.length();
  for ( unsigned i = 0; i < nLinks; i++ )
  {
    TMLink*           l   = links[i];
    const HepPoint3D& p   = l->position();
    unsigned          pos = (unsigned)floor( p.y() / _binSize );
 
    //...Why is this needed?...
    pos %= _nBins;
 
    ++_bins[pos];
    _links[pos]->append( l );
  }
}
 
void THistogram::fillPhi( const AList<TMLink>& links ) {
  _all            = (AList<TMLink>&)links;
  unsigned nLinks = links.length();
  double   offset = _binSize / 4.;
  for ( unsigned i = 0; i < nLinks; i++ )
  {
    TMLink*           l   = links[i];
    const HepPoint3D& p   = l->position();
    float             phi = atan2( p.y(), p.x() ) + M_PI;
    //  std::cout<<"atan "<<atan2(-1., -1.)<<std::endl;
    unsigned pos = (unsigned)floor( ( phi + offset ) / _binSize );
 
    //  std::cout <<"layer "<<l->wire()->layerId()<<" cell "<<l->wire()->localId()
    //        <<" x "<<p.x()<<" y "<<p.y()<<" pos "<<pos<<" xypos hit
    //"<<l->hit()->xyPosition()<<" xypos wire "<<l->wire()->xyPosition()<<" drift
    //"<<l->hit()->drift(0)<< "+-" <<l->hit()->dDrift(0)<<std::endl;    std::cout<<"binsize
    //"<<_binSize<<" offset "<<offset<<" phi "<<phi<<std::endl;
    //...Why is this needed?...
    pos %= _nBins;
 
    ++_bins[pos];
    _links[pos]->append( l );
  }
}
 
void THistogram::remove( const AList<TMLink>& links ) {
  for ( unsigned i = 0; i < _nBins; i++ )
  {
    _links[i]->remove( links );
    _bins[i] = _links[i]->length();
  }
  _all.remove( links );
}
 
AList<TMLink> THistogram::contents( unsigned center, unsigned width ) const {
  AList<TMLink> links;
  for ( int i = -(int)width; i <= (int)width; i++ )
  { links.append( *bin( (int)center + i ) ); }
  return links;
}
 
AList<TMLink> THistogram::contents( int start, int end ) const {
  AList<TMLink> links;
  for ( int i = start; i <= end; i++ ) links.append( *bin( i ) );
  return links;
}
 
AList<TSegment0> THistogram::clusters0( void ) const {
  AList<TSegment0> list;
 
  //...Serach for empty bin...
  unsigned begin = 0;
  while ( _bins[begin] > 0 ) begin++;
  if ( begin == _nBins ) return list;
 
  //...Start searching...
  unsigned loop = 0;
  while ( loop < _nBins )
  {
    ++loop;
    unsigned id = ( begin + loop ) % _nBins;
    if ( _bins[id] )
    {
      unsigned   size = 0;
      TSegment0* c    = new TSegment0();
      while ( _bins[id] )
      {
        if ( _bins[id] ) ++size;
        c->append( *_links[id] );
        ++loop;
        id = ( begin + loop ) % _nBins;
        if ( loop == _nBins ) break;
      }
      list.append( c );
    }
  }
  return list;
}
 
AList<TSegment0> THistogram::segments0( void ) const {
 
  //...Obtain raw clusters...
  AList<TSegment0> list = clusters0();
  unsigned         n    = list.length();
  if ( n == 0 ) return list;
 
  //...Examine each cluster...
  AList<TSegment0> splitted;
  for ( unsigned i = 0; i < n; i++ )
  {
    TSegment0* c = list[i];
 
    AList<TSegment0> newClusters = c->split();
    if ( newClusters.length() == 0 )
    {
      c->solveDualHits();
      continue;
    }
 
    list.append( newClusters );
    splitted.append( c );
#ifdef TRKRECO_DEBUG_DETAIL
    c->dump( "hits", "    " );
    std::cout << "    ... splitted as" << std::endl;
    for ( unsigned j = 0; j < newClusters.length(); j++ )
    {
      std::cout << "    " << j << " : ";
      newClusters[j]->dump( "hits" );
    }
#endif
  }
  list.remove( splitted );
  HepAListDeleteAll( splitted );
 
  return list;
}
 
AList<TSegment> THistogram::clusters( void ) const {
  AList<TSegment> list;
  //    std::cout<<"enter clusters"<<std::endl;
  //...Serach for empty bin...
  unsigned begin = 0;
  while ( _bins[begin] > 0 ) begin++;
  if ( begin == _nBins ) return list;
 
  //...Start searching...
  unsigned loop = 0;
  while ( loop < _nBins )
  {
    ++loop;
    unsigned id = ( begin + loop ) % _nBins;
    if ( _bins[id] )
    {
      unsigned  size = 0;
      TSegment* c    = new TSegment();
      while ( _bins[id] )
      {
        if ( _bins[id] ) ++size;
        c->append( *_links[id] );
        ++loop;
        id = ( begin + loop ) % _nBins;
        if ( loop == _nBins ) break;
      }
      list.append( c );
    }
  }
  return list;
}
 
AList<TSegment> THistogram::segments( void ) const {
  // yuany
  /*
  for (unsigned i = 0; i < _nBins; i++) {
      std::cout<<"i "<<i<<" bins[i] "<<_bins[i]<<std::endl;
    }
  */
  //...Obtain raw clusters...
  AList<TSegment> list = clusters();
  unsigned        n    = list.length();
  if ( n == 0 ) return list;
 
  //...Examine each cluster...
  AList<TSegment> splitted;
  for ( unsigned i = 0; i < n; i++ )
  {
    TSegment* c = list[i];
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    base segment : ";
    c->dump( "hits" );
#endif
 
    AList<TSegment> newClusters = c->split();
    if ( newClusters.length() == 0 )
    {
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout << "    ... Solving dual hits" << std::endl;
#endif
      c->solveDualHits();
      continue;
    }
 
    list.append( newClusters );
    splitted.append( c );
#ifdef TRKRECO_DEBUG_DETAIL
    c->dump( "hits", "    " );
    std::cout << "    ... splitted as" << std::endl;
    for ( unsigned j = 0; j < newClusters.length(); j++ )
    {
      std::cout << "    " << j << " : ";
      newClusters[j]->dump( "hits" );
    }
#endif
  }
  list.remove( splitted );
  HepAListDeleteAll( splitted );
 
  // yuany
  //    n = list.length();
  //    for (unsigned i = 0; i < n; i++) {
  //    TSegment * c = list[i];
  //    std::cout << "    base segment : ";
  //    c->dump("hits");
  //    }
 
  return list;
}