FTSuperLayer.cxx

Go to the documentation of this file.

#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/StatusCode.h"
#include <iostream>
 
#include "MdcGeomSvc/IMdcGeomSvc.h"
#include "MdcGeomSvc/MdcGeoLayer.h"
 
#include "FTGeom.h"
#include "FTLayer.h"
#include "FTSegment.h"
#include "FTSuperLayer.h"
#include "FTWire.h"
 
// ############################################################# //
 
FTSuperLayer::FTSuperLayer( const int superLayerID, const int firstLayerID,
                            const int firstWireID, const int Nlayer, const int NWire )
    : _superLayerId( superLayerID )
    , _complicated_segments( 0 )
    , _wireHits( 0 )
    , _singleHits( 0 )
    , _segments( 0 )
    , _firstLayerId( firstLayerID )
    , _firstWireId( firstWireID )
    , _Nlayer( Nlayer )
    , _Nwire( NWire ) {}
 
const FTList<FTWire*>& FTSuperLayer::wireHits() { return _wireHits; }
const int              FTSuperLayer::nWire() const { return _Nwire; }
const int              FTSuperLayer::nLayer() const { return _Nlayer; }
const int              FTSuperLayer::localMaxId() const { return ( _Nwire - 1 ); }
const int              FTSuperLayer::layerMaxId() const { return ( _Nlayer - 1 ); }
FTList<FTSegment*>&    FTSuperLayer::segments() { return _segments; }
const int              FTSuperLayer::superLayerId() const { return _superLayerId; }
FTList<FTSegment*>&    FTSuperLayer::complecated_segments() { return _complicated_segments; }
void                   FTSuperLayer::appendHit( FTWire* h ) { _wireHits.push_back( h ); }
 
// #undef inline
 
// ############################################################# //
 
MdcParameter* FTSuperLayer::param = MdcParameter::instance();
 
const std::array<float, 11> FTSuperLayer::_maxDphi = { 11.9, 9.72, 6.26, 4.86, 3.81, 2.37,
                                                       2.08, 1.76, 1.53, 2.33, 1.88 };
 
void FTSuperLayer::clear() {
  for ( auto w : _wireHits ) { w->state( FTWireHitInvalid ); }
  _wireHits.clear();
 
  _singleHits.clear();
 
  // for ( auto s : _segments ) { delete s; }
  _segments.clear();
 
  for ( auto s : _complicated_segments ) { delete s; }
  _complicated_segments.clear();
}
 
void FTSuperLayer::mkSegmentList() {
  clustering();
 
  FTList<FTSegment*> inner_short;
  FTList<FTSegment*> outer_short;
  FTList<FTSegment*> mid_short;
 
  for ( auto it = _segments.begin(); it != _segments.end(); )
  {
    FTSegment* s = *it;
 
    switch ( s->examine() )
    {
 
    case 0: // long simple
      it++;
      break;
 
    case 1: // inner short
      inner_short.push_back( s );
      it = _segments.erase( it ); // "it" updated
      break;
 
    case 2: // outer short
      outer_short.push_back( s );
      it = _segments.erase( it ); // "it" updated
      break;
 
    case 3: // to be divided
      _complicated_segments.push_back( s );
      it = _segments.erase( it ); // "it" updated
      break;
 
    default:
      mid_short.push_back( s );
      it = _segments.erase( it ); // "it" updated
      break;
    }
 
    // make sure that "it" has been updated
  }
 
  connect_short_segments( inner_short, outer_short );
  connect_singleHit( inner_short );
  connect_singleHit( outer_short );
  connect_singleHit( mid_short );
 
  for ( auto i : _segments )
  {
    i->update();
 
#ifndef OnlineMode
    i->printout();
#endif
  }
}
 
void FTSuperLayer::reduce_noise( FTList<FTList<float>>& Estime ) {
  if ( _superLayerId > 9 ) return;
  if ( _segments.size() > 30 ) return;
 
  float tdc[30][4] = { 0 };
  float r[4]       = { 0 };
  float phi[4]     = { 0 };
  int   wireId[4]  = { 0 };
 
  IMdcGeomSvc* mdcGeomSvc;
  StatusCode   sc = Gaudi::svcLocator()->service( "MdcGeomSvc", mdcGeomSvc );
 
  float T1 = 0.5 * 0.1 * ( mdcGeomSvc->Layer( 4 * _superLayerId + 0 )->PCSiz() ) / 0.004;
  float T2 = 0.5 * 0.1 * ( mdcGeomSvc->Layer( 4 * _superLayerId + 1 )->PCSiz() ) / 0.004;
  float T3 = 0.5 * 0.1 * ( mdcGeomSvc->Layer( 4 * _superLayerId + 2 )->PCSiz() ) / 0.004;
  float T4 = 0.5 * 0.1 * ( mdcGeomSvc->Layer( 4 * _superLayerId + 3 )->PCSiz() ) / 0.004;
 
  auto seg_it = _segments.begin();
  for ( int i = 0; seg_it != _segments.end(); seg_it++, i++ )
  {
    float chi2 = -99;
    float t0c  = 0;
 
    FTList<FTWire*>& hits = ( *seg_it )->wireHits();
 
    for ( auto h : hits )
    {
 
      const int layerId = h->layer()->layerId();
 
      wireId[layerId % 4] = h->localId();
      phi[layerId % 4]    = h->phi();
 
      tdc[i][layerId % 4] = h->time();
      r[layerId % 4]      = h->layer()->r();
    }
 
    if ( tdc[i][0] != 0 && tdc[i][1] != 0 && tdc[i][2] != 0 && tdc[i][3] != 0 &&
         ( wireId[0] == wireId[1] && wireId[0] == wireId[2] && wireId[0] == wireId[3] ||
           wireId[0] == wireId[1] - 1 && wireId[0] == wireId[2] && wireId[1] == wireId[3] ) )
    {
      float r0 = r[0] - r[1] - ( r[2] - r[1] ) / 2;
      float r1 = -( r[2] - r[1] ) / 2;
      float r2 = ( r[2] - r[1] ) / 2;
      float r3 = r[3] - r[2] + ( r[2] - r[1] ) / 2;
 
      float t_i  = tdc[i][0] + tdc[i][2];
      float t_j  = tdc[i][1] + tdc[i][3];
      float l_j  = T2 + T4;
      float r_i  = r0 + r2;
      float r_j  = r1 + r3;
      float r_2k = r0 * r0 + r1 * r1 + r2 * r2 + r3 * r3;
      float rt_i = r0 * tdc[i][0] + r2 * tdc[i][2];
      float rt_j = r1 * tdc[i][1] + r3 * tdc[i][3];
      float rl_j = r1 * T2 + r3 * T4;
 
      float deno = 4 * r_2k - 2 * ( r_i * r_i + r_j * r_j );
 
      if ( deno != 0 )
      {
        float Pa = ( 4 * ( rt_i - rt_j + rl_j ) - ( t_i + t_j - l_j ) * ( r_i - r_j ) -
                     ( t_i - t_j + l_j ) * ( r_i + r_j ) ) /
                   deno;
        float Pb  = 0.25 * ( t_i - t_j + l_j - ( r_i + r_j ) * Pa );
        float Ang = fabs( 90. - fabs( atan( Pa ) * 180. / 3.141593 ) );
 
        t0c = -0.25 * ( ( r_i - r_j ) * Pa - t_i - t_j + l_j );
 
        float pa1 = ( tdc[i][0] - tdc[i][2] ) / ( r0 - r2 );
        float pa2 = ( T4 - T2 - tdc[i][3] + tdc[i][1] ) / ( r3 - r1 );
 
        chi2 =
            ( tdc[i][0] - t0c - r0 * Pa - Pb ) * ( tdc[i][0] - t0c - r0 * Pa - Pb ) +
            ( T2 - tdc[i][1] + t0c - r1 * Pa - Pb ) * ( T2 - tdc[i][1] + t0c - r1 * Pa - Pb ) +
            ( tdc[i][2] - t0c - r2 * Pa - Pb ) * ( tdc[i][2] - t0c - r2 * Pa - Pb ) +
            ( T4 - tdc[i][3] + t0c - r3 * Pa - Pb ) * ( T4 - tdc[i][3] + t0c - r3 * Pa - Pb );
      }
    }
 
    if ( chi2 < ( param->_chi2_segfit ) && _superLayerId > 4 && chi2 > -1 )
    { Estime[_superLayerId].push_back( t0c ); }
 
    // continue; // FIXME
    // _complicated_segments.push_back( _segments[i] );
    // for ( int k = 0; k < 4; k++ ) { tdc[i][k] = 0; }
    // n = _segments.remove( i );
  }
}
 
// cluster hits to segments
void FTSuperLayer::clustering() {
  //                  +---+---+
  //                  | 4 | 5 |
  //                +-+-+---+-+-+    r
  //  Neighbor ID   | 2 | * | 3 |    ^
  //                +-+-+---+-+-+    |
  //                  | 0 | 1 |      +--> -phi
  //                  +---+---+
  //
  if ( !_wireHits.size() ) return;
 
  for ( auto h : _wireHits ) h->chk_left_and_right();
 
  for ( auto hptr : _wireHits )
  {
    if ( hptr->stateAND( FTWireHitAppendedOrInvalid ) ) continue;
 
    FTList<FTWire*> hits;
    hits.push_back( hptr );
    hptr->stateOR( FTWireHitAppended );
 
    for ( int j = 0; j < hits.size(); j++ )
    { // loop over hits in a cluster
      FTWire*            h       = hits[j];
      const unsigned int checked = h->state();
      // const unsigned int * mptr = _neighborsMask;
      auto neighbors = h->getNeighbors();
      // check 6 neighbors
      for ( auto n : neighbors )
      {
        if ( n->stateAND( FTWireHitAppendedOrInvalid ) ) continue;
        hits.push_back( n );
        n->stateOR( FTWireHitAppended );
      }
    }
 
    if ( hits.size() != 1 ) { _segments.push_back( new FTSegment( this, hits ) ); }
    else { _singleHits.push_back( hptr ); }
  }
}
 
// to match two phi-near short segments and connect them
void FTSuperLayer::connect_short_segments( FTList<FTSegment*>& inner_short,
                                           FTList<FTSegment*>& outer_short ) {
  for ( auto inner_it = inner_short.begin();
        inner_it != inner_short.end(); ) // no iner_it++ since we need to erase something
  {
    auto inner = *inner_it;
 
    const FTWire*  in_outerBoundHit = inner->outerBoundHits().front();
    const FTLayer* in_outerBound    = in_outerBoundHit->layer();
    float          in_outerPhi      = inner->outgoingPhi();
 
    auto  min_dphi_iter = outer_short.end();
    float min_dphi      = 2 * M_PI;
    int   dLayer_cache  = -1;
    float out_innerPhi, D_phi;
 
    for ( auto outer_it = outer_short.begin(); outer_it != outer_short.end(); outer_it++ )
    {
      const int out_innerLayerId =
          ( *outer_it )->innerBoundHits().front()->layer()->localLayerId();
      int dLayer_cache_tmp = out_innerLayerId - in_outerBound->localLayerId();
 
      if ( dLayer_cache_tmp < 1 ) continue;
 
      out_innerPhi = ( *outer_it )->incomingPhi();
      D_phi        = fabs( in_outerPhi - out_innerPhi );
      if ( D_phi > M_PI ) D_phi = 2 * M_PI - D_phi;
      if ( D_phi < min_dphi )
      {
        min_dphi      = D_phi;
        min_dphi_iter = outer_it;
        dLayer_cache  = dLayer_cache_tmp;
      }
    }
 
    if ( min_dphi_iter == outer_short.end() )
    {
      inner_it++;
      continue;
    }
 
    switch ( dLayer_cache )
    {
    case 1:
      if ( min_dphi > _maxDphi[_superLayerId] * M_PI / 180. )
      {
        inner_it++;
        continue;
      }
      break;
    case 2:
      if ( min_dphi > _maxDphi[_superLayerId] * M_PI / 120. )
      {
        inner_it++;
        continue;
      }
      break;
    case 3:
      if ( min_dphi > _maxDphi[_superLayerId] * M_PI / 80. )
      {
        inner_it++;
        continue;
      }
      break;
    default: inner_it++; continue;
    }
 
    inner->connect_outer( ( *min_dphi_iter )->wireHits(),
                          ( *min_dphi_iter )->outerBoundHits() );
 
    outer_short.erase( min_dphi_iter );
    inner_it = inner_short.erase( inner_it );
    _segments.push_back( inner );
  }
}
 
void FTSuperLayer::connect_singleHit( FTList<FTSegment*>& short_sgmnts ) {
  int minLength = 0;
  if ( superLayerId() == 2 || superLayerId() == 3 || superLayerId() == 4 ||
       superLayerId() == 9 || superLayerId() == 10 )
    minLength = 1;
 
  // for ( int i = 0; i < short_sgmnts.size(); i++ )
  for ( FTSegment* s : short_sgmnts )
  {
    const FTWire*  outerBoundHit = s->outerBoundHits().front();
    const FTWire*  innerBoundHit = s->innerBoundHits().front();
    const FTLayer* outerBound    = outerBoundHit->layer();
    const FTLayer* innerBound    = innerBoundHit->layer();
    float          outerPhi      = s->outgoingPhi();
    float          innerPhi      = s->incomingPhi();
 
    auto  min_dphi_iter = _singleHits.end();
    float min_dphi      = 2 * M_PI;
    int   dLayer_cache  = -1;
 
    bool inner_flag_cache     = false;
    bool inner_flag_cache_set = false;
 
    for ( auto h_iter = _singleHits.begin(); h_iter != _singleHits.end(); h_iter++ )
    {
      const FTWire* h                = *h_iter;
      const int     hLocalLayerId    = h->layer()->localLayerId();
      int           dLayer_cache_tmp = -1;
      bool          inner_flag_cache_tmp;
      float         D_phi;
 
      if ( ( dLayer_cache_tmp = innerBound->localLayerId() - hLocalLayerId ) > 0 )
      {
        D_phi                = fabs( innerPhi - h->phi() );
        inner_flag_cache_tmp = true;
      }
      else if ( ( dLayer_cache_tmp = hLocalLayerId - outerBound->localLayerId() ) > 0 )
      {
        D_phi                = fabs( outerPhi - h->phi() );
        inner_flag_cache_tmp = false;
      }
      else { continue; }
      if ( D_phi > M_PI ) D_phi = 2 * M_PI - D_phi;
 
      if ( D_phi < min_dphi )
      {
        min_dphi             = D_phi;
        min_dphi_iter        = h_iter;
        dLayer_cache         = dLayer_cache_tmp;
        inner_flag_cache     = inner_flag_cache_tmp;
        inner_flag_cache_set = true;
      }
    }
 
    int not_append = 1;
    if ( min_dphi_iter != _singleHits.end() )
    {
      switch ( dLayer_cache )
      {
      case 1:
        if ( min_dphi < _maxDphi[_superLayerId] * M_PI / 180. ) not_append = 0;
        break;
      case 2:
        if ( min_dphi < _maxDphi[_superLayerId] * M_PI / 120. ) not_append = 0;
        break;
      case 3:
        if ( min_dphi < _maxDphi[_superLayerId] * M_PI / 80. ) not_append = 0;
        break;
      default: break;
      }
    }
 
    if ( !not_append )
    {
      if ( !inner_flag_cache_set )
      {
        std::cerr << "FTSuperLayer::connect_singleHit: inner_flag_cache is not set"
                  << std::endl;
        abort();
      }
 
      if ( inner_flag_cache ) { s->connect_inner( *min_dphi_iter ); }
      else { s->connect_outer( *min_dphi_iter ); }
      _segments.push_back( s );
      _singleHits.erase( min_dphi_iter );
    }
    else
    {
      if ( outerBound->localLayerId() - innerBound->localLayerId() > minLength )
      { _segments.push_back( s ); }
      else
      {
        if ( !_superLayerId ) _complicated_segments.push_back( s );
        else delete s;
      }
    }
  }
}
 
void FTSuperLayer::reAppendSalvage() {
  for ( auto it = _segments.begin(); it != _segments.end(); )
  {
    if ( ( *it )->track() ) it++;
    else
    {
      _complicated_segments.push_back( *it );
      it = _segments.erase( it );
    }
  }
}
 
#ifdef FZISAN_DEBUG
void FTSuperLayer::showBinary( unsigned src ) {
  unsigned mask = 512;
  while ( mask >>= 1 ) { unsigned bit = ( mask & src ) ? 1 : 0; }
}
#endif