TCosmicFitter Class Reference

A class to fit a TTrackBase object to a helix. More...

#include <TCosmicFitter.h>

Inheritance diagram for TCosmicFitter:

Public Member Functions
	TCosmicFitter (const std::string &name)
	Constructor.
virtual	~TCosmicFitter ()
	Destructor.
void	dump (const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
	dumps debug information.
int	fit (TTrackBase &) const
int	fit (TTrackBase &, float t0Offset) const
int	fitWithCathode (TTrackBase &, float t0Offset=0., float windowSize=0.6, int SysCorr=0)
	TCosmicFitter (const std::string &name)
	Constructor.
virtual	~TCosmicFitter ()
	Destructor.
void	dump (const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
	dumps debug information.
int	fit (TTrackBase &) const
int	fit (TTrackBase &, float t0Offset) const
int	fitWithCathode (TTrackBase &, float t0Offset=0., float windowSize=0.6, int SysCorr=0)
	TCosmicFitter (const std::string &name)
	Constructor.
virtual	~TCosmicFitter ()
	Destructor.
void	dump (const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
	dumps debug information.
int	fit (TTrackBase &) const
int	fit (TTrackBase &, float t0Offset) const
int	fitWithCathode (TTrackBase &, float t0Offset=0., float windowSize=0.6, int SysCorr=0)
Public Member Functions inherited from TMFitter
	TMFitter (const std::string &name)
	Constructor.
virtual	~TMFitter ()
	Destructor.
const std::string &	name (void) const
	returns name.
void	dump (const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
	TMFitter (const std::string &name)
	Constructor.
virtual	~TMFitter ()
	Destructor.
const std::string &	name (void) const
	returns name.
void	dump (const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
	TMFitter (const std::string &name)
	Constructor.
virtual	~TMFitter ()
	Destructor.
const std::string &	name (void) const
	returns name.
void	dump (const std::string &message=std::string(""), const std::string &prefix=std::string("")) const

Additional Inherited Members
Protected Member Functions inherited from TMFitter
void	fitDone (TTrackBase &) const
	sets the fitted flag. (Bad implementation)
void	fitDone (TTrackBase &) const
	sets the fitted flag. (Bad implementation)
void	fitDone (TTrackBase &) const
	sets the fitted flag. (Bad implementation)

Detailed Description

A class to fit a TTrackBase object to a helix.

Definition at line 44 of file InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TCosmicFitter.h.

Constructor & Destructor Documentation

TCosmicFitter() [1/3]

TCosmicFitter::TCosmicFitter

(

const std::string &

name

)

Constructor.

Definition at line 78 of file TCosmicFitter.cxx.

    : TMFitter( name ), m_pmgnIMF( nullptr ) {}

~TCosmicFitter() [1/3]

TCosmicFitter::~TCosmicFitter ( )

virtual

Destructor.

Definition at line 81 of file TCosmicFitter.cxx.

{}

TCosmicFitter() [2/3]

TCosmicFitter::TCosmicFitter

(

const std::string &

name

)

Constructor.

~TCosmicFitter() [2/3]

virtual TCosmicFitter::~TCosmicFitter ( )

virtual

Destructor.

TCosmicFitter() [3/3]

TCosmicFitter::TCosmicFitter

(

const std::string &

name

)

Constructor.

~TCosmicFitter() [3/3]

virtual TCosmicFitter::~TCosmicFitter ( )

virtual

Destructor.

Member Function Documentation

dump() [1/3]

void TCosmicFitter::dump	(	const std::string &	message = std::string(""),
		const std::string &	prefix = std::string("") ) const

dumps debug information.

dump() [2/3]

void TCosmicFitter::dump	(	const std::string &	message = std::string(""),
		const std::string &	prefix = std::string("") ) const

dumps debug information.

dump() [3/3]

void TCosmicFitter::dump	(	const std::string &	message = std::string(""),
		const std::string &	prefix = std::string("") ) const

dumps debug information.

fit() [1/6]

int TCosmicFitter::fit ( TTrackBase & b ) const

virtual

Implements TMFitter.

Definition at line 83 of file TCosmicFitter.cxx.

                                            {
 
  //  double t0_bes;
  //  TrkReco::t0(t0_bes);
  //  const double t0_bes = TrkReco::t0();
 
  //...Already fitted ?...
  if ( b.fitted() ) return TFitAlreadyFitted;
 
  int err = fit( b, 0. );
  if ( !err ) b._fitted = true;
 
  return err;
}

Referenced by fit().

fit() [2/6]

int TCosmicFitter::fit ( TTrackBase & ) const

virtual

Implements TMFitter.

fit() [3/6]

int TCosmicFitter::fit ( TTrackBase & ) const

virtual

Implements TMFitter.

fit() [4/6]

int TCosmicFitter::fit	(	TTrackBase &	b,
		float	t0Offset ) const

Definition at line 98 of file TCosmicFitter.cxx.

                                                            {
 
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    TCosmicFitter::fit ..." << std::endl;
#endif
 
  IMdcCalibFunSvc* l_mdcCalFunSvc;
  Gaudi::svcLocator()->service( "MdcCalibFunSvc", l_mdcCalFunSvc );
 
  //...Check # of hits...
  if ( b.links().length() < 5 ) return TFitErrorFewHits;
  unsigned nValid = 0;
  for ( unsigned i = 0; i < b.links().length(); i++ )
  {
    unsigned state = b.links()[i]->hit()->state();
    if ( state & WireHitInvalidForFit ) continue;
    if ( state & WireHitFittingValid ) ++nValid;
  }
  if ( nValid < 5 ) return TFitErrorFewHits;
 
  //...Type check...
  //    if (b.type() != Track) return TFitUnavailable;
  if ( b.objectType() != Track ) return TFitUnavailable;
  TTrack& t = (TTrack&)b;
 
  // read t0 from TDS-------------------------------------//
  double       _t0_bes = -1.;
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "TCosmicFitter" );
 
  IDataProviderSvc* eventSvc = NULL;
  Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
 
  SmartDataPtr<RecEsTimeCol> aevtimeCol( eventSvc, "/Event/Recon/RecEsTimeCol" );
  if ( aevtimeCol )
  {
    RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
    _t0_bes                         = ( *iter_evt )->getTest();
    //      cout<<"  tev: "<<setw(4)<<_t0_bes<<endl;
  }
  else { log << MSG::WARNING << "Could not find RecEsTimeCol" << endmsg; }
  //    if (_t0_bes < 7.) _t0_bes = 7.;
  //----------------------------------------------------//
 
  //...Setup...
  Vector a( 5 ), da( 5 );
  a = t.helix().a();
  Vector    dxda( 5 );
  Vector    dyda( 5 );
  Vector    dzda( 5 );
  Vector    dDda( 5 );
  Vector    dchi2da( 5 );
  SymMatrix d2chi2d2a( 5, 0 );
  double    chi2;
  // double chi2Old = 10e99;
  double       chi2Old     = DBL_MAX;
  const double convergence = Convergence;
  bool         allAxial    = true;
  Matrix       e( 3, 3 );
  Vector       f( 3 );
  int          err    = 0;
  double       factor = 1.0; // jtanaka0715
 
  Vector maxDouble( 5 );
  for ( unsigned i = 0; i < 5; i++ ) maxDouble[i] = ( FLT_MAX );
 
  //...Fitting loop...
  unsigned nTrial = 0;
  while ( nTrial < NTrailMax )
  {
 
    //...Set up...
    chi2 = 0.;
    for ( unsigned j = 0; j < 5; j++ ) dchi2da[j] = 0.;
    d2chi2d2a = SymMatrix( 5, 0 );
 
#if CAL_TOF_HELIX
    // cal the time pass through the chamber
    const float Rad   = 81; // cm
    float       dRho  = t.helix().a()[0];
    float       Lproj = sqrt( Rad * Rad - dRho * dRho );
    float       tlmd  = t.helix().a()[4];
    float       fct   = sqrt( 1. + tlmd * tlmd );
#endif
 
    //...Loop with hits...
    unsigned i = 0;
    while ( TMLink* l = t.links()[i++] )
    {
      const TMDCWireHit& h = *l->hit();
 
      //...Check state...
      if ( h.state() & WireHitInvalidForFit ) continue;
      if ( !( h.state() & WireHitFittingValid ) ) continue;
 
      //...Check wire...
      if ( !nTrial )
        if ( h.wire()->stereo() ) allAxial = false;
 
      //...Cal. closest points...
      int doSagCorrection = 0;
      //        if(nTrial && !allAxial ) doSagCorrection = 1; //Liuqg
      t.approach( *l, doSagCorrection );
      double            dPhi    = l->dPhi();
      const HepPoint3D& onTrack = l->positionOnTrack();
      const HepPoint3D& onWire  = l->positionOnWire();
 
#ifdef TRKRECO_DEBUG_DETAIL
//      std::cout << "    in fit " << onTrack << ", " << onWire;
//      h.dump();
#endif
 
      //...Obtain drift distance and its error...
      unsigned leftRight = WireHitRight;
      if ( onWire.cross( onTrack ).z() < 0. ) leftRight = WireHitLeft;
      double distance  = h.drift( leftRight );
      double eDistance = h.dDrift( leftRight );
      //...
      if ( nTrial && !allAxial )
      {
        int side = leftRight;
 
        double Tfly = _t0_bes / 220. * ( 110. - onWire.y() );
#if CAL_TOF_HELIX
        float Rad_wire  = sqrt( x[0] * x[0] + x[1] * x[1] ); // cm
        float L_wire    = sqrt( Rad_wire * Rad_wire - dRho * dRho );
        float flyLength = Lproj - L_wire;
        if ( x[1] < 0 ) flyLength = Lproj + L_wire;
        Tfly = flyLength / 29.98 * sqrt( 1.0 + ( 0.105 / 0.5 ) * ( 0.105 / 0.5 ) ) *
               fct; // approxiate... cm/ns
#endif
        float time = l->hit()->reccdc()->tdc - Tfly;
 
        int    wire    = h.wire()->localId();
        int    layerId = h.wire()->layerId();
        float  dist    = distance;
        float  edist   = eDistance;
        double T0      = l_mdcCalFunSvc->getT0( layerId, wire );
        double rawadc  = 0.;
        rawadc         = l->hit()->reccdc()->adc;
 
        double timeWalk  = l_mdcCalFunSvc->getTimeWalk( layerId, rawadc );
        double drifttime = time - T0 - timeWalk;
        dist  = l_mdcCalFunSvc->driftTimeToDist( drifttime, layerId, wire, side, 0.0 );
        edist = l_mdcCalFunSvc->getSigma( layerId, side, dist, 0.0 );
        dist  = dist / 10.0; // mm->cm
        edist = edist / 10.0;
        //         if(layerId<20) edist = 9999.;
 
        // zsl         calcdc_driftdist_(&doPropDelay,&wire,&side,p,x,&time,&dist,&edist);
 
        distance  = (double)dist;
        eDistance = (double)edist;
 
        l->drift( distance, 0 ); // update
        l->drift( distance, 1 );
        l->dDrift( eDistance, 0 );
        l->dDrift( eDistance, 1 );
        l->tof( Tfly );
      }
      double eDistance2 = eDistance * eDistance;
 
      //...Residual...
      HepVector3D v         = onTrack - onWire;
      double      vmag      = v.mag();
      double      dDistance = vmag - distance;
 
      //...dxda...
      this->dxda( *l, t.helix(), dPhi, dxda, dyda, dzda, doSagCorrection );
 
      //...Chi2 related...
      // dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag;
      HepVector3D vw = h.wire()->direction();
      dDda = ( vmag > 0. ) ? ( ( v.x() * ( 1. - vw.x() * vw.x() ) - v.y() * vw.x() * vw.y() -
                                 v.z() * vw.x() * vw.z() ) *
                                   dxda +
                               ( v.y() * ( 1. - vw.y() * vw.y() ) - v.z() * vw.y() * vw.z() -
                                 v.x() * vw.y() * vw.x() ) *
                                   dyda +
                               ( v.z() * ( 1. - vw.z() * vw.z() ) - v.x() * vw.z() * vw.x() -
                                 v.y() * vw.z() * vw.y() ) *
                                   dzda ) /
                                 vmag
                           : Vector( 5, 0 );
      if ( vmag <= 0.0 )
      {
        std::cout << "    in fit " << onTrack << ", " << onWire;
        h.dump();
      }
      dchi2da += ( dDistance / eDistance2 ) * dDda;
      d2chi2d2a += vT_times_v( dDda ) / eDistance2;
      double pChi2 = dDistance * dDistance / eDistance2;
      chi2 += pChi2;
 
      //...Store results...
      l->update( onTrack, onWire, leftRight, pChi2 );
    }
 
    //...Check condition...
    double change = chi2Old - chi2;
    if ( fabs( change ) < convergence ) break;
    // if (change < 0.) break;
    // jtanaka -- from traffs -- Ozaki-san added this part to traffs.
    if ( change < 0. )
    {
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout << "chi2Old, chi2=" << chi2Old << " " << chi2 << std::endl;
#endif
      // change to the old value.
      a += factor * da;
      //      a[2] = 0.000000001;
      t._helix->a( a );
 
      chi2 = 0.;
      for ( unsigned j = 0; j < 5; j++ ) dchi2da[j] = 0.;
      d2chi2d2a = SymMatrix( 5, 0 );
 
      //...Loop with hits...
      unsigned i = 0;
      while ( TMLink* l = t.links()[i++] )
      {
        const TMDCWireHit& h = *l->hit();
 
        //...Check state...
        if ( h.state() & WireHitInvalidForFit ) continue;
        if ( !( h.state() & WireHitFittingValid ) ) continue;
 
        //...Check wire...
        if ( !nTrial )
          if ( h.wire()->stereo() ) allAxial = false;
 
        //...Cal. closest points...
        int doSagCorrection = 0;
        //            if( nTrial  && !allAxial ) doSagCorrection = 1;  //Liuqg
        t.approach( *l, doSagCorrection );
        double            dPhi    = l->dPhi();
        const HepPoint3D& onTrack = l->positionOnTrack();
        const HepPoint3D& onWire  = l->positionOnWire();
 
#ifdef TRKRECO_DEBUG_DETAIL
        // std::cout << "    in fit in case of change < 0. " << onTrack << ", " << onWire;
        // h.dump();
#endif
 
        //...Obtain drift distance and its error...
        unsigned leftRight = WireHitRight;
        if ( onWire.cross( onTrack ).z() < 0. ) leftRight = WireHitLeft;
        double distance  = h.drift( leftRight );
        double eDistance = h.dDrift( leftRight );
        if ( nTrial && !allAxial )
        {
          int side = leftRight;
 
          double Tfly = _t0_bes / 220. * ( 110. - onWire.y() );
#if CAL_TOF_HELIX
          float Rad_wire  = sqrt( x[0] * x[0] + x[1] * x[1] ); // cm
          float L_wire    = sqrt( Rad_wire * Rad_wire - dRho * dRho );
          float flyLength = Lproj - L_wire;
          if ( x[1] < 0 ) flyLength = Lproj + L_wire;
          Tfly = flyLength / 29.98 * sqrt( 1.0 + ( 0.105 / 10 ) * ( 0.105 / 10 ) ) *
                 fct; // approxiate... cm/ns
#endif
 
          float time = l->hit()->reccdc()->tdc - Tfly;
 
          int   wire    = h.wire()->localId();
          int   layerId = h.wire()->layerId();
          float dist    = distance;
          float edist   = eDistance;
          //               int doPropDelay = 1; //
 
          double T0        = l_mdcCalFunSvc->getT0( layerId, wire );
          double rawadc    = 0.;
          rawadc           = l->hit()->reccdc()->adc;
          double timeWalk  = l_mdcCalFunSvc->getTimeWalk( layerId, rawadc );
          double drifttime = time - T0 - timeWalk;
          dist  = l_mdcCalFunSvc->driftTimeToDist( drifttime, layerId, wire, side, 0.0 );
          edist = l_mdcCalFunSvc->getSigma( layerId, side, dist, 0.0 );
          dist  = dist / 10.0; // mm->cm
          edist = edist / 10.0;
          //           if (layerId<20) edist = 9999.;
 
          // zsl calcdc_driftdist_(&doPropDelay,&wire,&side,p,x,&time,&dist,&edist);
 
          distance  = (double)dist;
          eDistance = (double)edist;
 
          l->drift( distance, 0 ); // update
          l->drift( distance, 1 );
          l->dDrift( eDistance, 0 );
          l->dDrift( eDistance, 1 );
          l->tof( Tfly );
        }
        double eDistance2 = eDistance * eDistance;
 
        //...Residual...
        HepVector3D v         = onTrack - onWire;
        double      vmag      = v.mag();
        double      dDistance = vmag - distance;
 
        //...dxda...
        this->dxda( *l, t.helix(), dPhi, dxda, dyda, dzda, doSagCorrection );
 
        //...Chi2 related...
        // dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag;
        HepVector3D vw = h.wire()->direction();
        dDda = ( vmag > 0. ) ? ( ( v.x() * ( 1. - vw.x() * vw.x() ) - v.y() * vw.x() * vw.y() -
                                   v.z() * vw.x() * vw.z() ) *
                                     dxda +
                                 ( v.y() * ( 1. - vw.y() * vw.y() ) - v.z() * vw.y() * vw.z() -
                                   v.x() * vw.y() * vw.x() ) *
                                     dyda +
                                 ( v.z() * ( 1. - vw.z() * vw.z() ) - v.x() * vw.z() * vw.x() -
                                   v.y() * vw.z() * vw.y() ) *
                                     dzda ) /
                                   vmag
                             : Vector( 5, 0 );
        if ( vmag <= 0.0 )
        {
          std::cout << "    in fit " << onTrack << ", " << onWire;
          h.dump();
        }
        dchi2da += ( dDistance / eDistance2 ) * dDda;
        d2chi2d2a += vT_times_v( dDda ) / eDistance2;
        double pChi2 = dDistance * dDistance / eDistance2;
        chi2 += pChi2;
 
        //...Store results...
        l->update( onTrack, onWire, leftRight, pChi2 );
      }
      // break;
      factor *= 0.75;
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout << "factor = " << factor << std::endl;
      std::cout << "chi2 = " << chi2 << std::endl;
#endif
      if ( factor < 0.01 ) break;
    }
 
    chi2Old = chi2;
 
    //...Cal. helix parameters for next loop...
    if ( allAxial )
    {
      f     = dchi2da.sub( 1, 3 );
      e     = d2chi2d2a.sub( 1, 3 );
      f     = solve( e, f );
      da[0] = f[0];
      da[1] = f[1];
      da[2] = f[2];
      da[3] = 0.;
      da[4] = 0.;
    }
    else { da = solve( d2chi2d2a, dchi2da ); }
 
#ifdef TRKRECO_DEBUG_DETAIL
    // std::cout << "    fit " << nTrial << " : " << da << std::endl;
    // std::cout << "        d2chi " << d2chi2d2a << std::endl;
    // std::cout << "        dchi2 " << dchi2da << std::endl;
#endif
 
    a -= factor * da;
    //  a[2] = 0.000000001;
    t._helix->a( a );
    ++nTrial;
  }
 
  //...Cal. error matrix...
  SymMatrix Ea( 5, 0 );
  unsigned  dim;
  if ( allAxial )
  {
    dim = 3;
    SymMatrix Eb( 3, 0 ), Ec( 3, 0 );
    Eb       = d2chi2d2a.sub( 1, 3 );
    Ec       = Eb.inverse( err );
    Ea[0][0] = Ec[0][0];
    Ea[0][1] = Ec[0][1];
    Ea[0][2] = Ec[0][2];
    Ea[1][1] = Ec[1][1];
    Ea[1][2] = Ec[1][2];
    Ea[2][2] = Ec[2][2];
  }
  else
  {
    dim = 5;
    Ea  = d2chi2d2a.inverse( err );
  }
 
  //...Store information...
  if ( !err )
  {
    t._helix->a( a );
    t._helix->Ea( Ea );
  }
  else { err = TFitFailed; }
 
  t._ndf  = nValid - dim;
  t._chi2 = chi2;
  // t._fitted = true;
 
  return err;
}

fit() [5/6]

int TCosmicFitter::fit	(	TTrackBase &	,
		float	t0Offset ) const

fit() [6/6]

int TCosmicFitter::fit	(	TTrackBase &	,
		float	t0Offset ) const

fitWithCathode() [1/3]

int TCosmicFitter::fitWithCathode	(	TTrackBase &	,
		float	t0Offset = 0.,
		float	windowSize = 0.6,
		int	SysCorr = 0 )

fitWithCathode() [2/3]

int TCosmicFitter::fitWithCathode	(	TTrackBase &	,
		float	t0Offset = 0.,
		float	windowSize = 0.6,
		int	SysCorr = 0 )

fitWithCathode() [3/3]

int TCosmicFitter::fitWithCathode	(	TTrackBase &	,
		float	t0Offset = 0.,
		float	windowSize = 0.6,
		int	SysCorr = 0 )

---

The documentation for this class was generated from the following files:

• InstallArea/x86_64-el9-gcc13-dbg/include/TrkReco/TCosmicFitter.h
• InstallArea/x86_64-el9-gcc13-opt/include/TrkReco/TCosmicFitter.h
• Reconstruction/TrkReco/include/TrkReco/TCosmicFitter.h
• TCosmicFitter.cxx