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)
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
virtual int	fit (TTrackBase &) const =0

Additional Inherited Members
Protected Member Functions inherited from TMFitter
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 45 of file TCosmicFitter.h.

Constructor & Destructor Documentation

TCosmicFitter()

TCosmicFitter::TCosmicFitter

(

const std::string &

name

)

Constructor.

Definition at line 84 of file TCosmicFitter.cxx.

                                                   : TMFitter(name) {
//jialk
   StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF); 
  if(scmgn!=StatusCode::SUCCESS) { 
    std::cout<< "Unable to open Magnetic field service"<<std::endl;
  }
 
}

~TCosmicFitter()

TCosmicFitter::~TCosmicFitter ( )

virtual

Destructor.

Definition at line 93 of file TCosmicFitter.cxx.

                              {
}

Member Function Documentation

dump()

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

dumps debug information.

fit() [1/2]

int TCosmicFitter::fit ( TTrackBase &  b ) const

virtual

Implements TMFitter.

Definition at line 97 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 TBuilderCosmic::buildStereo(), TBuilder0::fit(), and fit().

fit() [2/2]

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

Definition at line 113 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" << endreq;
    }
//    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;
}

fitWithCathode()

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:

• TCosmicFitter.h
• TCosmicFitter.cxx