FTSegment.cxx

Go to the documentation of this file.

#include "FTSegment.h"
 
#include "FTLayer.h"
#include "FTList.h"
#include "FTSuperLayer.h"
#include "FTTrack.h"
#include "FTWire.h"
#include "MdcParameter.h"
 
#include <math.h>
 
FTSegment::FTSegment( FTSuperLayer* super, FTList<FTWire*>& hits )
    : _wireHits( hits )
    , _superLayer( super )
    , _innerBoundHits( 0 )
    , _outerBoundHits( 0 )
    , _firstCallUpdate3D( true )
    , _sList( 0 )
    , _zList( 0 ) {
  StatusCode scmgn = Gaudi::svcLocator()->service( "MagneticFieldSvc", m_pmgnIMF );
  if ( scmgn != StatusCode::SUCCESS )
  { std::cout << "Unable to open Magnetic field service" << std::endl; }
}
 
void FTSegment::connect_outer( const FTList<FTWire*>& outerHits,
                               const FTList<FTWire*>& outerBound ) {
  _wireHits.insert( _wireHits.end(), outerHits.begin(), outerHits.end() );
  _outerBoundHits.clear();
  _outerBoundHits.insert( _outerBoundHits.end(), outerBound.begin(), outerBound.end() );
}
 
void FTSegment::connect_inner( const FTList<FTWire*>& innerHits,
                               const FTList<FTWire*>& innerBound ) {
  _wireHits.insert( _wireHits.end(), innerHits.begin(), innerHits.end() );
  _innerBoundHits.clear();
  _innerBoundHits.insert( _innerBoundHits.end(), innerBound.begin(), innerBound.end() );
}
 
void FTSegment::connect_outer( FTWire* h ) {
  _wireHits.push_back( h );
  _outerBoundHits.clear();
  _outerBoundHits.push_back( h );
}
 
void FTSegment::connect_inner( FTWire* h ) {
  _wireHits.push_back( h );
  _innerBoundHits.clear();
  _innerBoundHits.push_back( h );
}
 
FTSuperLayer*          FTSegment::superLayer() const { return _superLayer; }
FTList<FTWire*>&       FTSegment::wireHits() { return _wireHits; }
const FTList<FTWire*>& FTSegment::innerBoundHits() { return _innerBoundHits; }
const FTList<FTWire*>& FTSegment::outerBoundHits() { return _outerBoundHits; }
float                  FTSegment::outgoingX() const { return _outgoing_x; }
float                  FTSegment::outgoingY() const { return _outgoing_y; }
float                  FTSegment::incomingX() const { return _incoming_x; }
float                  FTSegment::incomingY() const { return _incoming_y; }
float                  FTSegment::r() const { return _r; }
 
float FTSegment::s() const {
  if ( _sList.size() ) return _sList[0];
  else
  {
    std::cout << "FTSegment::s() _sList.size() = 0 ! this = " << this << std::endl;
    return 0;
  }
}
 
float FTSegment::z() const {
  if ( _zList.size() ) return _zList[0];
  else
  {
    std::cout << "FTSegment::s() _zList.size() = 0 ! this = " << this << std::endl;
    return 0;
  }
}
 
float FTSegment::kappa() const { return _kappa; }
 
FTTrack* FTSegment::track() const {
  if ( _trackList.size() ) return _trackList.front();
  else return nullptr;
}
 
FTTrack* FTSegment::track( FTTrack* src ) {
  if ( _trackList.size() ) _trackList.front() = src;
  else _trackList.push_back( src );
  return _trackList.front();
}
 
float FTSegment::incomingPhi() const {
  if ( _innerBoundHits.size() > 1 )
  {
    auto  it      = _innerBoundHits.begin();
    float phi_1st = ( *it )->phi();
    it++;
    float phi_2nd = ( *it )->phi();
    float phi     = 0.5 * ( phi_1st + phi_2nd );
    if ( fabs( phi_1st - phi_2nd ) > M_PI )
    { return ( phi > M_PI ) ? ( phi - M_PI ) : ( phi + M_PI ); }
    return phi;
  }
  else return _innerBoundHits.front()->phi();
}
 
float FTSegment::outgoingPhi() const {
  if ( _outerBoundHits.size() > 1 )
  {
    auto  it      = _outerBoundHits.begin();
    float phi_1st = ( *it )->phi();
    it++;
    float phi_2nd = ( *it )->phi();
    float phi     = 0.5 * ( phi_1st + phi_2nd );
    if ( fabs( phi_1st - phi_2nd ) > M_PI )
    { return ( phi > M_PI ) ? ( phi - M_PI ) : ( phi + M_PI ); }
    return phi;
  }
  else return _outerBoundHits.front()->phi();
}
 
float FTSegment::attr_const() const {
  float x_len = _outgoing_x - _incoming_x;
  float y_len = _outgoing_y - _incoming_y;
  float phi_l = outgoingPhi() - incomingPhi();
  return sin( phi_l ) / sqrt( x_len * x_len + y_len * y_len );
}
 
MdcParameter* FTSegment::param = MdcParameter::instance();
 
int FTSegment::examine() {
  int stat = -1;
 
  array<FTList<FTWire*>, 4> wireHitsByLocalLayer;
 
  for ( auto& hptr : _wireHits )
  { wireHitsByLocalLayer[hptr->layer()->localLayerId()].push_back( hptr ); }
 
  float length_phi   = 0;
  int   innerBoundId = 0;
  int   outerBoundId = 0;
  int   AlreadyFound = 0;
 
  const int nLayer = _superLayer->nLayer();
 
  // to look for the innerBound and outBand of the hits
  for ( int i = 0; i != nLayer; i++ )
  {
    int N_wireHits = wireHitsByLocalLayer[i].size();
    if ( N_wireHits )
    {
      float localPhi_max = 0.;
      float localPhi_min = 7.;
 
      for ( int j = 0; j != N_wireHits; j++ )
      {
        float localPhi = wireHitsByLocalLayer[i][j]->phi();
        if ( localPhi > localPhi_max ) localPhi_max = localPhi;
        if ( localPhi < localPhi_min ) localPhi_min = localPhi;
      }
 
      float local_length_phi = localPhi_max - localPhi_min;
      if ( local_length_phi > M_PI ) local_length_phi = 2 * M_PI - local_length_phi;
      if ( local_length_phi > length_phi ) length_phi = local_length_phi;
      outerBoundId = i;
      AlreadyFound = 1;
    }
    else if ( !AlreadyFound ) { innerBoundId++; }
    else { break; }
  }
 
  // length_phi needed to decide more carefully
  if ( length_phi > FTSuperLayer::_maxDphi[_superLayer->superLayerId()] * 2 * M_PI / 180 )
  {           // by X.-R. Lu
    stat = 3; // to be divided
  }
  else
  {
    if ( innerBoundId > 0 )
    {
      stat = 2; // outer short
    }
    else
    {
      if ( outerBoundId < nLayer - 1 )
      {
        stat = 1; // inner short
      }
      else
      {
        stat = 0; // long simple
      }
    }
  }
  _innerBoundHits.insert( _innerBoundHits.end(), wireHitsByLocalLayer[innerBoundId].begin(),
                          wireHitsByLocalLayer[innerBoundId].end() );
  _outerBoundHits.insert( _outerBoundHits.end(), wireHitsByLocalLayer[outerBoundId].begin(),
                          wireHitsByLocalLayer[outerBoundId].end() );
  return stat;
}
 
void FTSegment::update() {
 
  updateInOutXY();
 
  float in_r  = _innerBoundHits.front()->layer()->r();
  float out_r = _outerBoundHits.front()->layer()->r();
  float sbX   = _incoming_x - _outgoing_x;
  float sbY   = _incoming_y - _outgoing_y;
 
  _kappa = 2. * ( -1. / 2.99792458 / m_pmgnIMF->getReferField() ) *
           ( _outgoing_x * _incoming_y - _outgoing_y * _incoming_x ) /
           ( in_r * out_r * sqrt( sbX * sbX + sbY * sbY ) );
 
  _r = 0.5 * ( in_r + out_r );
}
 
int FTSegment::update3D( FTTrack* track ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc ).ignore();
 
  MsgStream log( msgSvc, "FTFinder" );
 
  if ( _firstCallUpdate3D )
  {
    updateInOutXY();
    _firstCallUpdate3D = false;
  }
 
  const Lpav& la = track->lpav();
  double      d;
  double      inS;
 
  const FTLayer* bound = _innerBoundHits.front()->layer();
 
  log << MSG::DEBUG << "started inS calculation with layer limit " << bound->limit() << endmsg;
  if ( !( la.sd( (double)bound->r(), (double)_incoming_x, (double)_incoming_y,
                 (double)bound->limit(), // multiplied with an factor due to bad rphi fit
                 inS, d ) ) )
    return 0;
 
  double inZ = bound->z( d );
  double outS;
 
  bound = _outerBoundHits.front()->layer();
  log << MSG::DEBUG << "started outS calculation with layer limit " << bound->limit()
      << endmsg;
  if ( !( la.sd( (double)bound->r(), (double)_outgoing_x, (double)_outgoing_y,
                 (double)bound->limit(), // multiplied with an factor due to bad rphi fit
                 outS, d ) ) )
    return 0;
 
  double outZ = bound->z( d );
  _sList.push_back( inS + outS );
  _zList.push_back( inZ + outZ );
  _trackList.push_back( track );
  log << MSG::DEBUG << "-----> coupled! DONE! " << endmsg;
 
  return 1;
}
 
void FTSegment::linkStereoSegments() {
  float    min_D_z  = 9999.;
  float    S        = 0.;
  float    Z        = 0.;
  FTTrack* selected = nullptr;
 
  int j = 0;
  for ( auto t : _trackList )
  {
    float z_tmp = _zList[j];
    float s_tmp = _sList[j];
    float D_z   = fabs( t->d_z( s_tmp, z_tmp ) );
    if ( D_z < min_D_z )
    {
      selected = t;
      min_D_z  = D_z;
      S        = s_tmp;
      Z        = z_tmp;
    }
    j++;
  }
  if ( selected ) { selected->append_stereo( this, S, Z ); }
}
 
#ifndef OnlineMode
void FTSegment::printout( void ) {
 
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
 
  MsgStream log( msgSvc, "FTSegment" );
  log << MSG::DEBUG << "pintout of one segment:" << endmsg;
 
  for ( auto h : _wireHits )
  {
    log << MSG::DEBUG << "hit: layer " << h->layer()->layerId() << " phi " << h->phi()
        << endmsg;
  }
}
#endif
 
void FTSegment::updateInOutXY() {
  int n = _innerBoundHits.size();
  int m = _outerBoundHits.size();
 
  if ( n > 1 )
  {
    // std::list doesn't support random access, so we have to use an iterator
    auto  it = _innerBoundHits.begin();
    float x0 = ( *it )->x(), y0 = ( *it )->y();
    it++;
    float x1 = ( *it )->x(), y1 = ( *it )->y();
 
    _incoming_x = 0.5 * ( x0 + x1 );
    _incoming_y = 0.5 * ( y0 + y1 );
  }
  else
  {
    _incoming_x = _innerBoundHits.front()->x();
    _incoming_y = _innerBoundHits.front()->y();
  }
 
  if ( m > 1 )
  {
    // std::list doesn't support random access, so we have to use an iterator
    auto  it = _outerBoundHits.begin();
    float x0 = ( *it )->x(), y0 = ( *it )->y();
    it++;
    float x1 = ( *it )->x(), y1 = ( *it )->y();
 
    _outgoing_x = 0.5 * ( x0 + x1 );
    _outgoing_y = 0.5 * ( y0 + y1 );
  }
  else
  {
    _outgoing_x = _outerBoundHits.front()->x();
    _outgoing_y = _outerBoundHits.front()->y();
  }
}