TCosmicFitter Class Reference

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

#include <TCosmicFitter.h>

Inheritance diagram for TCosmicFitter:

List of all members.

Public Member Functions
void	dump (const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
	dumps debug information.
void	dump (const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
	dumps debug information.
int	fit (TTrackBase &, float t0Offset) const
int	fit (TTrackBase &) const
int	fit (TTrackBase &, float t0Offset) const
int	fit (TTrackBase &) const
int	fitWithCathode (TTrackBase &, float t0Offset=0., float windowSize=0.6, int SysCorr=0)
int	fitWithCathode (TTrackBase &, float t0Offset=0., float windowSize=0.6, int SysCorr=0)
const std::string &	name (void) const
	returns name.
const std::string &	name (void) const
	returns name.
	TCosmicFitter (const std::string &name)
	Constructor.
	TCosmicFitter (const std::string &name)
	Constructor.
virtual	~TCosmicFitter ()
	Destructor.
virtual	~TCosmicFitter ()
	Destructor.
Protected Member Functions
void	fitDone (TTrackBase &) const
	sets the fitted flag. (Bad implementation)
void	fitDone (TTrackBase &) const
	sets the fitted flag. (Bad implementation)
Private Member Functions
int	dxda (const TMLink &link, const Helix &helix, double dPhi, HepVector &dxda, HepVector &dyda, HepVector &dzda, int doSagCorrection) const
	calculates dXda. 'link' and 'dPhi' are inputs. Others are outputs.
int	dxda (const TMLink &link, const Helix &helix, double dPhi, HepVector &dxda, HepVector &dyda, HepVector &dzda, int doSagCorrection) const
	calculates dXda. 'link' and 'dPhi' are inputs. Others are outputs.
Private Attributes
IMagneticFieldSvc *	m_pmgnIMF
IMagneticFieldSvc *	m_pmgnIMF

---

Detailed Description

A class to fit a TTrackBase object to a helix.

---

Constructor & Destructor Documentation

TCosmicFitter::TCosmicFitter (  const std::string &  name  )

Constructor. : 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 (   )  [virtual]

Destructor. { }

TCosmicFitter::TCosmicFitter (  const std::string &  name  )

Constructor.

virtual TCosmicFitter::~TCosmicFitter (   )  [virtual]

Destructor.

---

Member Function Documentation

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

dumps debug information. Reimplemented from TMFitter.

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

dumps debug information. Reimplemented from TMFitter.

int TCosmicFitter::dxda (  const TMLink &  link, const Helix &  helix, double  dPhi, HepVector &  dxda, HepVector &  dyda, HepVector &  dzda, int  doSagCorrection )  const [private]

calculates dXda. 'link' and 'dPhi' are inputs. Others are outputs.

int TCosmicFitter::dxda (  const TMLink &  link, const Helix &  helix, double  dPhi, HepVector &  dxda, HepVector &  dyda, HepVector &  dzda, int  doSagCorrection )  const [private]

calculates dXda. 'link' and 'dPhi' are inputs. Others are outputs.

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

int TCosmicFitter::fit (  TTrackBase &   )  const [virtual]

Implements TMFitter.

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

{ #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; dist = l_mdcCalFunSvc->rawTimeNoTOFToDist(time, 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; // dist = l_mdcCalFunSvc->rawTimeNoTOFToDist(time, 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; }

int TCosmicFitter::fit (  TTrackBase &   )  const [virtual]

Implements TMFitter. { // 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; }

void TMFitter::fitDone (  TTrackBase &   )  const [protected, inherited]

sets the fitted flag. (Bad implementation)

void TMFitter::fitDone (  TTrackBase &   )  const [protected, inherited]

sets the fitted flag. (Bad implementation) { t._fitted = true; }

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

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

const std::string& TMFitter::name (  void   )  const [inherited]

returns name.

const std::string & TMFitter::name (  void   )  const [inline, inherited]

returns name. { return _name; }

---

Member Data Documentation

IMagneticFieldSvc* TCosmicFitter::m_pmgnIMF [private]

IMagneticFieldSvc* TCosmicFitter::m_pmgnIMF [private]

---

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

• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/InstallArea/include/TrkReco/TrkReco/TCosmicFitter.h
• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Reconstruction/TrkReco/TrkReco-00-06-12/TrkReco/TCosmicFitter.h
• /data0/mdestefa/bes3soft/Boss_6.5.5/dist/6.5.5/Reconstruction/TrkReco/TrkReco-00-06-12/src/TCosmicFitter.cxx

---