FTTrack.cxx

Go to the documentation of this file.

#include "FTTrack.h"
#include "FTFinder.h"
#include "FTSegment.h"
#include "FTTrack.h"
#include "FTWire.h"
#include "MdcParameter.h"
 
#include <algorithm>
#include <iterator>
#include <stdlib.h>
#include <tuple>
 
#ifndef OnlineMode
#  include "AIDA/IAxis.h"
#  include "AIDA/IHistogram1D.h"
#  include "AIDA/IHistogram2D.h"
#  include "AIDA/IHistogram3D.h"
#  include "AIDA/IHistogramFactory.h"
 
extern AIDA::IHistogram1D* g_sigmaxy;
extern AIDA::IHistogram1D* g_sigmaz;
extern AIDA::IHistogram1D* g_chi2xy;
extern AIDA::IHistogram1D* g_chi2sz;
 
#endif
 
// ##################################################### //
 
// #ifdef FTTrack_NO_INLINE
// #  define inline
// #else
// #  undef inline
// #endif
 
FTTrack::FTTrack( FTList<FTSegment*>& axial_segments, float kappa, float chi2_kappa )
    : _la( nullptr )
    , _za( nullptr )
    , _helix( nullptr )
    , _axial_segments( axial_segments )
    , _kappa( kappa )
    , _chi2_kappa( chi2_kappa ) {
  StatusCode scmgn = Gaudi::svcLocator()->service( "MagneticFieldSvc", m_pmgnIMF );
  if ( scmgn != StatusCode::SUCCESS )
  { std::cout << "Unable to open MagneticField service" << std::endl; }
}
 
FTTrack::~FTTrack() {
  // for ( auto i : _axial_segments ) delete i; // FIXME: this is a memory leak
  // delete _la;
  // delete _za;
  // delete _helix;
}
 
const FTList<FTSegment*>& FTTrack::axial_segments() { return _axial_segments; }
const FTList<FTSegment*>& FTTrack::stereo_segments() { return _stereo_segments; }
float                     FTTrack::kappa_tmp() const { return _kappa; }
float                     FTTrack::chi2_kappa_tmp() const { return _chi2_kappa; }
const Lpav&               FTTrack::lpav() const { return *_la; }
const zav&                FTTrack::Zav() const { return *_za; }
Helix*                    FTTrack::helix() const { return _helix; }
float                     FTTrack::SigmaZ( float z ) { return _SigmaZ += z; }
float                     FTTrack::SigmaS( float s ) { return _SigmaS += s; }
float                     FTTrack::SigmaSZ( float sz ) { return _SigmaSZ += sz; }
float                     FTTrack::SigmaSS( float ss ) { return _SigmaSS += ss; }
void                      FTTrack::setFTFinder( FTFinder* finder ) { m_ftFinder = finder; }
 
void FTTrack::append_stereo_cache( FTSegment* src ) {
  _stereo_segments_cache.push_back( src );
}
 
void FTTrack::append_stereo( FTSegment* src, float s = 0, float z = 0 ) {
  _stereo_segments.push_back( src );
  _SigmaS += s;
  _SigmaZ += z;
  _SigmaSZ += s * z;
  _SigmaSS += s * s;
}
 
float FTTrack::d_z( float s, float z ) const {
  int n = _stereo_segments.size();
  if ( !n ) return 0;
  float Dz = ( n == 1 ) ? z - s * _SigmaZ / _SigmaS
                        : z - ( s * ( _SigmaSZ - _SigmaS * _SigmaZ ) + _SigmaSS * _SigmaZ -
                                _SigmaS * _SigmaSZ ) /
                                  ( _SigmaSS - _SigmaS * _SigmaS );
  return ( Dz < 10 * ( 7 - n ) ) ? Dz : 9999.;
}
 
void FTTrack::updateSZ() {
  int n = _stereo_segments_cache.size();
  if ( n == 1 )
  {
    FTSegment* s     = _stereo_segments_cache[0];
    float      s_tmp = s->s();
    float      z_tmp = s->z();
    if ( d_z( s_tmp, z_tmp ) < 9998. )
    {
      _SigmaS += s_tmp;
      _SigmaZ += z_tmp;
      _SigmaSS += s_tmp * s_tmp;
      _SigmaSZ += s_tmp * z_tmp;
      _stereo_segments.push_back( s );
    }
    // _stereo_segments_cache.clear();
  }
  else if ( n )
  {
    _stereo_segments_by_superLayer.push_back( _stereo_segments_cache );
    // _stereo_segments_cache = FTList<FTSegment*>( 3 );
  }
  _stereo_segments_cache.clear();
}
 
// #undef inline
 
// ##################################################### //
 
MdcParameter* FTTrack::param = MdcParameter::instance();
 
int FTTrack::r_phiFit() {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc ).ignore();
 
  MsgStream log( msgSvc, "FTFinder" );
 
  // static const float alpha(333.564095);
  if ( fabs( _kappa ) > 1.2 / param->_minPt ) return 0;
 
  _la   = new Lpav;
  int n = _axial_segments.size();
  for ( int i = 0; i != n; i++ )
  {
    FTList<FTWire*>& hits = _axial_segments[i]->wireHits();
 
    for ( auto it = hits.begin();
          it != hits.end(); ) // erase operation occurs, so no "it++" here
    {
      FTWire& h = **it;
      if ( h.stateAND( FTWireFittingInvalid ) )
      {
        it = hits.erase( it );
        continue;
      }
 
      double par = h.distance() / ( 0.25 * h.layer()->csize() );
      _la->add_point( (double)h.x(), (double)h.y(), exp( -par * par ) );
 
      it++;
    }
  }
 
  double    chi2 = _la->fit();
  HepVector a    = _la->Hpar( m_ftFinder->pivot );
 
  log << MSG::DEBUG << " chi2/_la cut(1): " << chi2 << " / " << _la->nc() << " a1("
      << param->_minDr << "): " << a( 1 ) << " a3(" << 1.05 / param->_minPt << "):" << a( 3 )
      << endmsg;
 
  if ( chi2 / _la->nc() > 1. ) return 0;
  if ( fabs( a( 3 ) ) > ( 1.05 / param->_minPt ) ) return 0;
  if ( fabs( a( 1 ) ) > param->_minDr ) return 0;
 
  _la->clear();
 
  log << MSG::DEBUG << " passed chi2/_la, a(3) and a(1) cut" << endmsg;
 
  int layer0 = 0;
  for ( int i = 0; i != n; i++ )
  {
    FTList<FTWire*>& hits = _axial_segments[i]->wireHits();
    if ( !_axial_segments[i]->superLayer()->superLayerId() ) layer0 = 1;
    int m = hits.size();
 
    for ( auto it = hits.begin(); it != hits.end(); it++ )
    {
      FTWire&     h     = **it;
      const float x     = h.x();
      const float y     = h.y();
      double      d0    = _la->d( (double)x, (double)y );
      float       delta = h.distance() / h.layer()->r();
      if ( fabs( d0 ) > 0.7 * h.layer()->csize() ) continue; // remove bad hits
      if ( d0 > 0 )
      { // left or right
        d0 = -d0;
      }
      else { delta = -delta; }
      _la->add_point( x - delta * y, y + delta * x, 1 );
    }
  }
 
  if ( layer0 == 0 )
  { // salvage hits from complecated segments in layer0
    FTList<FTSegment*>& salvage = m_ftFinder->superLayer( 0 )->complecated_segments();
    HepVector           center  = _la->center();
    const float         xc      = center( 1 );
    const float         yc      = center( 2 );
    const float         rc      = a( 1 ) + ( -1. / 2.99792458 / m_pmgnIMF->getReferField() ) /
                                  a( 3 ); // rho+drho(signed)
    int nn = salvage.size();
 
    for ( int i = 0; i != nn; i++ )
    {
      int              salvaged = 0;
      FTList<FTWire*>& hits     = salvage[i]->wireHits();
 
      for ( auto it = hits.begin(); it != hits.end(); it++ )
      {
        FTWire& h = **it;
        float   x = h.x();
        float   y = h.y();
        float   r = h.layer()->r();
        if ( ( y * xc - x * yc ) / ( r * rc ) < 0.707 ) break;
        double d0    = _la->d( (double)x, (double)y );
        float  delta = h.distance() / r;
        if ( fabs( d0 ) > 0.7 * h.layer()->csize() ) continue; // remove bad hits
        if ( d0 > 0 )
        { // left or right
          h.stateOR( FTWireHitRight );
          d0 = -d0;
        }
        else
        {
          h.stateOR( FTWireHitLeft );
          delta = -delta;
        }
        _la->add_point( x - delta * y, y + delta * x, 1 );
        salvaged = 1;
      }
 
      if ( salvaged )
      {
        _axial_segments.push_back( salvage[i] );
        break;
      }
    }
  }
  chi2 = _la->fit(); // refit using drift distance
 
  _stereo_segments.clear();
  _stereo_segments_cache.clear();
  _stereo_segments_by_superLayer.clear();
 
  _za      = new zav;
  _SigmaS  = 0.;
  _SigmaZ  = 0.;
  _SigmaSS = 0.;
  _SigmaSZ = 0.;
  return 1;
}
 
int FTTrack::r_phiReFit( float vx, float vy, int vtx_flag ) {
  if ( vtx_flag ) _la->fit( vx, vy, 20 * _la->chisq() / _la->nc() );
 
  _la->clear();
 
  int n = _axial_segments.size();
  for ( int i = 0; i != n; i++ )
  {
    FTList<FTWire*>& hits = _axial_segments[i]->wireHits();
 
    for ( auto it = hits.begin(); it != hits.end(); it++ )
    {
      FTWire&     h        = **it;
      const float x        = h.x();
      const float y        = h.y();
      double      d0       = _la->d( (double)x, (double)y );
      float       cellsize = h.layer()->csize();
      if ( m_ftFinder->evtTiming )
      {
        float time = h.time() + m_ftFinder->evtTiming;
        if ( time < -100 ) continue; // discard the bad TDC time
        h.distance( m_ftFinder->t2x( h.layer(), time ) );
        // h.distance(time*40*0.0001);
      }
      float delta = h.distance() / h.layer()->r();
      if ( fabs( d0 ) > 0.5 * cellsize ) continue; // remove bad hits
      if ( d0 > 0 )
      { // left or right
        h.stateOR( FTWireHitRight );
        d0 = -d0;
      }
      else
      {
        h.stateOR( FTWireHitLeft );
        delta = -delta;
      }
      h.setChi2( h.distance() + d0 );
 
#ifndef OnlineMode
      // fill histogram
      g_sigmaxy->fill( h.distance() + d0, 1.0 );
#endif
 
      _la->add_point( x - delta * y, y + delta * x, 1 );
    }
  }
 
  HepVector b    = _la->Hpar( m_ftFinder->pivot );
  double    chi2 = _la->fit(); // refit using drift distance
 
#ifndef OnlineMode
  g_chi2xy->fill( chi2 / _la->nc(), 1.0 );
#endif
 
  return 1;
}
 
int FTTrack::r_phi2Fit( float vx, float vy, int vtx_flag ) {
  if ( vtx_flag ) _la->fit( vx, vy, 20 * _la->chisq() / _la->nc() );
  int n = _axial_segments.size();
  _la->clear();
 
  int k = 0, l = 0;
 
  float temp = 0;
  for ( int i = 0; i != n; i++ )
  {
    FTList<FTWire*>& hits = _axial_segments[i]->wireHits();
 
    for ( auto it = hits.begin(); it != hits.end(); it++ )
    {
      k++;
      FTWire& h = **it;
 
      const float x        = h.x();
      const float y        = h.y();
      double      d0       = _la->d( (double)x, (double)y );
      float       cellsize = h.layer()->csize();
      if ( m_ftFinder->evtTiming )
      {
        float time = h.time() + m_ftFinder->evtTiming;
        if ( time < -100 ) continue;
        h.distance( m_ftFinder->t2x( h.layer(), time ) );
      }
      float delta = h.distance() / h.layer()->r();
      if ( fabs( d0 ) > 0.5 * cellsize ) continue; // remove bad hits
      if ( h.distance() < 0.1 || h.distance() > 0.6 ) continue;
 
      if ( d0 > 0 )
      { // left or right
        h.stateOR( FTWireHitRight );
        d0 = -d0;
      }
      else
      {
        h.stateOR( FTWireHitLeft );
        delta = -delta;
      }
 
      _la->add_point( x - delta * y, y + delta * x, 1 );
 
      if ( temp < fabs( h.distance() + d0 ) )
      {
        temp = fabs( h.distance() + d0 );
        l    = k;
      }
    }
  }
  HepVector b    = _la->Hpar( m_ftFinder->pivot );
  double    chi2 = _la->fit(); // refit using drift distance
 
  if ( k > 5 ) { return l; }
  else { return -99; }
}
 
int FTTrack::r_phi3Fit( int l, float vx, float vy, int vtx_flag ) {
  if ( vtx_flag ) _la->fit( vx, vy, 20 * _la->chisq() / _la->nc() );
  int n = _axial_segments.size();
  _la->clear();
  int k = 0;
  for ( int i = 0; i ^ n; i++ )
  {
    FTList<FTWire*>& hits = _axial_segments[i]->wireHits();
 
    for ( auto it = hits.begin();
          it != hits.end(); ) // erase operation occurs, so no "it++" here
    {
      k++;
      FTWire& h = **it;
 
      const float x        = h.x();
      const float y        = h.y();
      double      d0       = _la->d( (double)x, (double)y );
      float       cellsize = h.layer()->csize();
      if ( m_ftFinder->evtTiming )
      {
        float time = h.time() + m_ftFinder->evtTiming;
        if ( time < -100 ) continue;
        h.distance( m_ftFinder->t2x( h.layer(), time ) );
        // h.distance(time*40*0.0001);
      }
      float delta = h.distance() / h.layer()->r();
      if ( fabs( d0 ) > 0.5 * cellsize ) continue; // remove bad hits
      if ( d0 > 0 )
      { // left or right
        h.stateOR( FTWireHitRight );
        d0 = -d0;
      }
      else
      {
        h.stateOR( FTWireHitLeft );
        delta = -delta;
      }
      if ( k == l )
      {
        it = hits.erase( it );
        continue;
      }
      _la->add_point( x - delta * y, y + delta * x, 1 );
 
      it++;
    }
  }
  HepVector b    = _la->Hpar( m_ftFinder->pivot );
  double    chi2 = _la->fit(); // refit using drift distance
  return 1;
}
 
int FTTrack::r_phi4Fit( float vx, float vy, int vtx_flag ) {
  if ( vtx_flag ) _la->fit( vx, vy, 20 * _la->chisq() / _la->nc() );
 
  _la->clear();
 
  int n = _axial_segments.size();
  for ( int i = 0; i != n; i++ )
  {
    FTList<FTWire*>& hits = _axial_segments[i]->wireHits();
 
    for ( auto it = hits.begin(); it != hits.end(); it++ )
    {
      FTWire* h = *it;
 
      const float x        = h->x();
      const float y        = h->y();
      double      d0       = _la->d( (double)x, (double)y );
      float       cellsize = h->layer()->csize();
      if ( m_ftFinder->evtTiming )
      {
        float time = h->time() + m_ftFinder->evtTiming;
        if ( time < -100 ) continue;
        h->distance( m_ftFinder->t2x( h->layer(), time ) );
      }
      float delta = h->distance() / h->layer()->r();
      if ( fabs( d0 ) > 0.5 * cellsize ) continue; // remove bad hits
      if ( d0 > 0 )
      { // left or right
        h->stateOR( FTWireHitRight );
        d0 = -d0;
      }
      else
      {
        h->stateOR( FTWireHitLeft );
        delta = -delta;
      }
 
      _la->add_point( x - delta * y, y + delta * x, 1 );
      h->setChi2( h->distance() + d0 );
    }
  }
  HepVector b    = _la->Hpar( m_ftFinder->pivot );
  double    chi2 = _la->fit(); // refit using drift distance
 
  return 1;
}
 
int FTTrack::linkStereoSegments() {
  // delete _stereo_segments_cache;
  int n = _stereo_segments_by_superLayer.size();
  for ( int i = 0; i != n; i++ )
  {
    FTList<FTSegment*>& segments = _stereo_segments_by_superLayer[i];
 
    int        m        = segments.size();
    float      min_D_z  = 9998.;
    float      S        = 0.;
    float      Z        = 0.;
    FTSegment* selected = nullptr;
    for ( int j = 0; j != m; j++ )
    {
      FTSegment* s     = segments[j];
      float      s_tmp = s->s();
      float      z_tmp = s->z();
      double     D_z   = fabs( d_z( s_tmp, z_tmp ) );
      if ( D_z < min_D_z )
      {
        selected = s;
        min_D_z  = D_z;
        S        = s_tmp;
        Z        = z_tmp;
      }
    }
    if ( selected )
    {
      _stereo_segments.push_back( selected );
      _SigmaS += S;
      _SigmaZ += Z;
      _SigmaSZ += S * Z;
      _SigmaSS += S * S;
    }
  }
 
  // for ( auto i : _stereo_segments_by_superLayer )
  // {
  //   for ( auto j : i ) delete j;
  // }
 
  return 0;
}
 
int FTTrack::s_zFit() {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc ).ignore();
 
  MsgStream log( msgSvc, "FTFinder" );
 
  HepVector a = _la->Hpar( m_ftFinder->pivot );
 
  int n = _stereo_segments.size();
  log << MSG::DEBUG << "number of stereo segments: " << n << endmsg;
  if ( n < ( param->_nseg ) )
  {
    a( 4 ) = -9999.;
    a( 5 ) = -9999.;
    _helix = new Helix( m_ftFinder->pivot, a );
    log << MSG::DEBUG << "cut by _nseg" << param->_nseg << endmsg;
    return 0;
  }
 
  FTList<double>  zList;
  FTList<double>  sList;
  FTList<FTWire*> hList;
  for ( int i = 0; i != n; i++ )
  {
    FTList<FTWire*>& hits = _stereo_segments[i]->wireHits();
 
    for ( auto it = hits.begin(); it != hits.end(); it++ )
    {
      FTWire& h = **it;
 
      double z;
      if ( !( h.z( *_la, z ) ) ) continue;
 
      double s        = _la->s( h.layer()->r() );
      float  cellsize = h.layer()->csize();
      log << MSG::DEBUG << "cellsize: " << cellsize << endmsg;
 
      if ( m_ftFinder->evtTiming )
      {
        float time = h.time() + m_ftFinder->evtTiming;
        if ( time < -100 ) continue; // discard the bad TDC time
        double distance = m_ftFinder->t2x( h.layer(), time );
        h.distance( distance );
        log << MSG::DEBUG << "m_ftFinder->evtTiming: " << m_ftFinder->evtTiming
            << " TDC time: " << h.time() << " distance: " << distance << endmsg;
      }
 
      double par = h.distance() / ( 0.25 * cellsize );
      _za->add( s, z, exp( -par * par ) );
      sList.push_back( s );
      zList.push_back( z );
      hList.push_back( &h );
    }
  }
 
  if ( hList.size() < ( param->_nlength ) )
  {
    a( 4 ) = -9999.;
    a( 5 ) = -9999.;
    _helix = new Helix( m_ftFinder->pivot, a );
    log << MSG::DEBUG << "cut by _nlength: " << param->_nlength << endmsg;
    return 0;
  }
 
  double chi2 = _za->calculate();
  _za->clear();
 
  // ----------------------------------
  for ( auto it = std::make_tuple( zList.begin(), sList.begin(), hList.begin() );
        it != std::make_tuple( zList.end(), sList.end(), hList.end() ); )
  {
    auto itz = std::get<0>( it );
    auto its = std::get<1>( it );
    auto ith = std::get<2>( it );
 
    double d          = _za->d( *its, *itz );
    float  z_distance = ( *ith )->distance_z();
 
    if ( fabs( fabs( d ) - z_distance ) > ( param->_z_cut1 ) )
    {
      log << MSG::DEBUG << "cut by _z_cut1: " << param->_z_cut1 << endmsg;
 
      itz = zList.erase( itz );
      its = sList.erase( its );
      ith = hList.erase( ith );
      it  = std::make_tuple( itz, its, ith );
      continue;
    }
 
    _za->add( *its, ( d > 0 ) ? *itz - z_distance : *itz + z_distance, 1. );
 
    it = std::make_tuple( ++itz, ++its, ++ith );
  }
 
  if ( _za->nc() < ( param->_nc ) )
  {
    a( 4 ) = -9999.;
    a( 5 ) = -9999.;
    _helix = new Helix( m_ftFinder->pivot, a );
    log << MSG::DEBUG << "cut by _nc: " << param->_nc << endmsg;
    return 0;
  }
 
  chi2 = _za->calculate();
  _za->clear();
 
  for ( auto it = std::make_tuple( zList.begin(), sList.begin(), hList.begin() );
        it != std::make_tuple( zList.end(), sList.end(), hList.end() );
        it = std::make_tuple( ++std::get<0>( it ), ++std::get<1>( it ), ++std::get<2>( it ) ) )
  {
    auto itz = std::get<0>( it );
    auto its = std::get<1>( it );
    auto ith = std::get<2>( it );
 
    double d          = _za->d( *its, *itz );
    float  z_distance = ( *ith )->distance_z();
    ( *ith )->setChi2( fabs( d ) - z_distance );
 
#ifndef OnlineMode
    // fill ntuple
    g_sigmaz->fill( fabs( d ) - z_distance, 1.0 );
#endif
 
    if ( fabs( fabs( d ) - z_distance ) > ( param->_z_cut2 ) ) continue;
    _za->add( *its, ( d > 0 ) ? *itz - z_distance : *itz + z_distance, 1. );
  }
 
  if ( _za->nc() < ( param->_nc ) )
  {
    a( 4 ) = -9999.;
    a( 5 ) = -9999.;
    _helix = new Helix( m_ftFinder->pivot, a );
    log << MSG::DEBUG << "cut by _nc" << param->_nc << endmsg;
    return 0;
  }
 
  chi2 = _za->calculate();
#ifndef OnlineMode
  g_chi2sz->fill( chi2 / _za->nc(), 1.0 );
#endif
 
  a( 4 ) = _za->b(); // dz
  a( 5 ) = _za->a(); // tanLambda
  _helix = new Helix( m_ftFinder->pivot, a );
  return 1;
}
 
int FTTrack::get_nhits() {
  int nhits = 0;
  for ( auto i : _axial_segments ) nhits += i->wireHits().size();
  for ( auto i : _stereo_segments ) nhits += i->wireHits().size();
  return nhits;
}
 
#ifndef OnlineMode
void FTTrack::printout() {
  for ( auto i : _axial_segments ) i->printout();
  for ( auto i : _stereo_segments ) i->printout();
}
#endif