THelixFitter Class Reference

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

#include <THelixFitter.h>

Inheritance diagram for THelixFitter:

List of all members.

Public Member Functions
	THelixFitter (const std::string &name)
	Constructor.
virtual	~THelixFitter ()
	Destructor.
void	dump (const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
	dumps debug information.
bool	fit2D (void) const
	sets/returns 2D flag.
bool	fit2D (bool)
bool	freeT0 (void) const
	sets/returns free T0 flag.
bool	freeT0 (bool)
unsigned	corrections (void) const
	sets/returns correctin flag.
unsigned	corrections (unsigned)
bool	sag (void) const
	sets/returns sag correction flag.
bool	sag (bool)
bool	propagation (void) const
	sets/returns propagation-delay correction flag.
bool	propagation (bool)
bool	tof (void) const
	sets/returns propagation-delay correction flag.
bool	tof (bool)
bool	tanl (void) const
	sets/returns tanLambda correction flag.
bool	tanl (bool)
double	preChi2 (void) const
	returns sum of chi2 before fit.
double	chi2 (void) const
	returns sum of chi2 aftter fit.
IMagneticFieldSvc *	getMagneticFieldPointer (void) const
int	fit (TTrackBase &) const
int	fit (TTrackBase &, double *pre_chi2, double *fitted_chi2) const
int	fit (TTrackBase &, float t0Offset, double *pre_chi2=NULL, double *fitted_chi2=NULL) const
int	fit (TTrackBase &, float &tev, float &tev_err, double *pre_chi2=NULL, double *fitted_chi2=NULL) const
const std::string &	name (void) const
	returns name.
Protected Member Functions
void	fitDone (TTrackBase &) const
	sets the fitted flag. (Bad implementation)
Private Member Functions
int	main (TTrackBase &, float t0Offset, double *pre_chi2=NULL, double *fitted_chi2=NULL) const
	main routine for fixed T0.
int	main (TTrackBase &, float &tev, float &tev_err, double *pre_chi2=NULL, double *fitted_chi2=NULL) const
	main routine for free T0.
void	drift (const TTrack &, TMLink &, float t0Offset, double &distance, double &itsError) const
	calculates drift distance and its error.
void	drift (const TTrack &, TMLink &, float t0Offset, double &distance, double &itsError, double &ddda) const
	calculates drift distance and its error for free T0 case.
int	dxda (const TMLink &link, const Helix &helix, double dPhi, HepVector &dxda, HepVector &dyda, HepVector &dzda) const
	calculates dXda. 'link' and 'dPhi' are inputs. Others are outputs.
Private Attributes
bool	_fit2D
bool	_freeT0
unsigned	_corrections
bool	_sag
int	_propagation
bool	_tof
bool	_tanl
double	_driftTime
double	_pre_chi2
double	_fitted_chi2
IMagneticFieldSvc *	m_pmgnIMF

---

Detailed Description

A class to fit a TTrackBase object to a helix.

Definition at line 48 of file THelixFitter.h.

---

Constructor & Destructor Documentation

THelixFitter::THelixFitter

(

const std::string &

name

)

Constructor.

Definition at line 122 of file THelixFitter.cxx.

: TMFitter(name),
  _fit2D(false),
  _freeT0(false),
  _sag(false),
  _propagation(false),
  _tof(false),
  _tanl(false),
  _pre_chi2(0.),
  _fitted_chi2(0.),
  m_pmgnIMF(0) {

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

}

THelixFitter::~THelixFitter

(

)

[virtual]

Destructor.

Definition at line 144 of file THelixFitter.cxx.

                            {
}

---

Member Function Documentation

double THelixFitter::chi2

(

void

)

const [inline]

returns sum of chi2 aftter fit.

Definition at line 305 of file THelixFitter.h.

References _fitted_chi2.

Referenced by main().

                             {
    return _fitted_chi2;
}

unsigned THelixFitter::corrections

(

unsigned

)

[inline]

Definition at line 293 of file THelixFitter.h.

References _corrections.

                                    {
    return _corrections = a;
}

unsigned THelixFitter::corrections

(

void

)

const [inline]

sets/returns correctin flag.

Definition at line 287 of file THelixFitter.h.

References _corrections.

                                    {
    return _corrections;
}

void THelixFitter::drift	(	const TTrack &	,
		TMLink &	,
		float	t0Offset,
		double &	distance,
		double &	itsError,
		double &	ddda
	)			const [private]

calculates drift distance and its error for free T0 case.

by wangjk, propagation correction

readout in the west of MDC

readout in the east of MDC

Definition at line 1973 of file THelixFitter.cxx.

References _propagation, _tof, MdcRec_wirhit::adc, TMDCWire::backwardPosition(), TMDCWireHit::dDrift(), check_raw_filter::dist, TMLink::dPhi(), TMDCWireHit::drift(), IMdcCalibFunSvc::driftTimeToDist(), TMDCWire::forwardPosition(), IMdcCalibFunSvc::getSigma(), IMdcCalibFunSvc::getTimeWalk(), TMLink::hit(), TMDCWire::layerId(), TMDCWire::localId(), msgSvc(), TMLink::positionOnTrack(), TMLink::positionOnWire(), TMDCWireHit::reccdc(), s, TMLink::setDriftTime(), t(), MdcRec_wirhit::tdc, tof(), Bes_Common::WARNING, TMDCWireHit::wire(), WireHitLeft, and WireHitRight.

                                         {
//=========================================
// be cautious here
    const TMDCWireHit & h = * l.hit();
    const HepPoint3D & onTrack = l.positionOnTrack();
    const HepPoint3D & onWire = l.positionOnWire();
    unsigned leftRight = WireHitRight;
    if (onWire.cross(onTrack).z() < 0.) leftRight = WireHitLeft;

    //...No correction...
    if ((tev == 0.) && (! _propagation) && (! _tof)) {
        distance = h.drift(leftRight);
        err = h.dDrift(leftRight);
        return;
    }

    //...TOF correction...
    float tof = 0.;
    if (_tof) {
        int imass = 3;
        float tl = t.helix().a()[4];
        float f = sqrt(1. + tl * tl);
        float s = fabs(t.helix().curv()) * fabs(l.dPhi()) * f;
        float p = f / fabs(t.helix().a()[2]);
//zsl   calcdc_tof2_(& imass, & p, & s, & tof);
        float Mass[5] = {0.000511,0.105,0.140,0.493677,0.98327};
        tof = s/29.98*sqrt(1.0+(Mass[imass-1]/p)*(Mass[imass-1]/p));
        //cout<<"tof:"<<tof<<endl;
    }

    //...T0 and propagation corrections...
    int wire = h.wire()->localId();
    int layerId = h.wire()->layerId();
//    int wire = h.wire()->id();
    int side = leftRight;


    const double zhit = onWire.z();
    const double vinner = 220.0e8; // unit is cm/s
    const double vouter = 240.0e8;
    
    double tprop = 0.;
    const double vprop = (layerId<8) ? vinner : vouter;
    if (0 == layerId%2){
      tprop = fabs((zhit - h.wire()->backwardPosition()[2]))/vprop; 
    }else{ 
      tprop = fabs(( zhit - h.wire()->forwardPosition()[2]))/vprop;
    }


//    if (side==0) side = -1;
//    HepVector3D tp = t.helix().momentum(l.dPhi());
//    float p[3] = {tp.x(),tp.y(),tp.z()};
//    float x[3] = {onWire.x(), onWire.y(), onWire.z()};
    float time = h.reccdc()->tdc + tev - tof - 1.0e9*tprop;
    //    float dist_p;
    //    float dist_m;
//    float dist;
//    float edist;

//  cout<<"tdc: "<<h.reccdc()->tdc<<"  tev: "<<tev<<"  tof: "<<tof<<endl;


    double EsT0 = -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();
      EsT0 = (*iter_evt)->getTest();
    }else{
      log << MSG::WARNING << "Could not find RecEsTimeCol" << endreq;
    }

    double rawTime = 0.;
    rawTime = h.reccdc()->tdc;
    double rawadc = 0.;
    rawadc =h.reccdc()->adc;
//    double timewalk = CalibFunSvc_->getTimeWalk(layerid, 0.0);
//    double timeDrift = rawTime - tof - tprop - EsTo - toWalk; 


//----cal drift----//
    IMdcCalibFunSvc* l_mdcCalFunSvc;
    Gaudi::svcLocator() -> service("MdcCalibFunSvc", l_mdcCalFunSvc);

    double timeWalk = l_mdcCalFunSvc->getTimeWalk(layerId, rawadc);
    double drifttime = rawTime - tof - 1.0e9*tprop - EsT0 - timeWalk;

    l.setDriftTime(drifttime);


    double dist;
    double edist;
    int prop = _propagation;

    //    //calculate derivative w.r.t. time in blute force way; need update 
    //    //in future to speed up
    //    float time_p = time + 0.1;
    //    calcdc_driftdist_(& prop,
    //                        & wire,
    //                        & side,
    //                        p,
    //                        x,
    //                        & time_p,
    //                        & dist_p,
    //                        & edist);
    //
    //        float time_m = time - 0.1;
    //        calcdc_driftdist_(& prop,
    //                        & wire,                 
    //                      & side,
    //                        p,
    //                        x,
    //                        & time_m,
    //                        & dist_m,
    //                        & edist);
    //    dddt = (dist_p - dist_m)*5.;
    //               std::cout << "side=" << side << std::endl;
    //               std::cout << "dddt=" << dddt << std::endl;

//    float dist2[2];
//    float sigma_d2[2];
//    float deriv2[2];
    float time_tmp = time;

    //----cal drift----//
//    IMdcCalibFunSvc* l_mdcCalFunSvc;
    Gaudi::svcLocator() -> service("MdcCalibFunSvc", l_mdcCalFunSvc);
  //  dist = l_mdcCalFunSvc->rawTimeNoTOFToDist(time_tmp, layerId, wire, side, 0.0);
    dist =  l_mdcCalFunSvc->driftTimeToDist(time_tmp,layerId, wire, side,0.0);//by liucy 2010/05/12
    edist = l_mdcCalFunSvc->getSigma(layerId, side, dist, 0.0);
    dist = dist/10.0; //mm->cm
    edist = edist/10.0;

    distance = dist;
    err = edist;

//    cout<<"dist: "<<dist<<" edist: "<<edist<<endl;

    float time_p = time - 0.1;
    float time_m = time + 0.1;
//    float dist_p = l_mdcCalFunSvc->rawTimeNoTOFToDist(time_p, layerId, wire, side, 0.0);
//    float dist_m = l_mdcCalFunSvc->rawTimeNoTOFToDist(time_m, layerId, wire, side, 0.0);
//    dist_p = dist_p/10.;
//    dist_m = dist_m/10.;
//    dddt = (dist_m - dist_p)/0.2;
//    dddt = 0;
//    dddt = 0.004;

//    float dist2[2] = {dist, dist};
//    float sigma_d2[2] = {edist, edist};
//    float deriv2[2] = {dddt_tmp, dddt_tmp}; // cm/ns
//    float deriv2[2] = {0.00, 0.00};
//    float deriv2[2] = {0.004, 0.004};
//cout<<"dddt_tmp: "<<dddt_tmp<<endl;

//zsl    calcdc_driftdist3_(& prop,
//                     & wire,                
//                     p,
//                     x,
//                     & time_tmp,
//                     dist2,
//                     sigma_d2, 
//                     deriv2);
    //n.b. input and output time are slightly different because of prop. 
    //delay corr. in driftdist3.
    //    calcdc_driftdist_(& prop,
    //                        & wire,
    //                        & side,
    //                        p,
    //                        x,
    //                        & time,
    //                        & dist,
    //                        & edist);

/*    if (side==-1) {
      //            std::cout << " " << std::endl;
      //            std::cout << dist << " " << dist2[0] << " " <<dist2[1] << std::endl;
      //            std::cout << edist << " " << sigma_d2[0] << " " << sigma_d2[1] << std::endl;
      //            std::cout << dddt  << " " << 0.001*deriv2[0] << std::endl;
      dist  = dist2[0];
      edist = sigma_d2[0];
      //dddt = 0.001*deriv2[0];
      dddt = deriv2[0];
    } else if(side== 1) {
      //            std::cout << " " << std::endl;
      //            std::cout << dist << " " << dist2[1] << " " <<dist2[0] << std::endl;
      //            std::cout << edist << " " << sigma_d2[1] << " " << sigma_d2[0] << std::endl;
      //            std::cout << dddt  << " " << 0.001*deriv2[1] << std::endl;
      dist  = dist2[1];
      edist = sigma_d2[1];
      //dddt = 0.001*deriv2[1];
      dddt = deriv2[1];
    }  

    distance = (double) dist;
//    std::cout << "time,distance,dddt="<<time<<" "<<distance<<"  "<<dddt<<std::endl;
    err = (double) edist;
*/
    return;
}

void THelixFitter::drift	(	const TTrack &	,
		TMLink &	,
		float	t0Offset,
		double &	distance,
		double &	itsError
	)			const [private]

calculates drift distance and its error.

by wangjk, propagation correction

Definition at line 1253 of file THelixFitter.cxx.

References _propagation, _tof, MdcRec_wirhit::adc, alpha, cos(), TMDCWireHit::dDrift(), check_raw_filter::dist, TMLink::dPhi(), TMDCWireHit::drift(), IMdcCalibFunSvc::driftTimeToDist(), IMagneticFieldSvc::getReferField(), IMdcCalibFunSvc::getSigma(), IMdcCalibFunSvc::getT0(), IMdcCalibFunSvc::getTimeWalk(), IMdcCalibFunSvc::getTprop(), TMLink::hit(), TMDCWire::layerId(), TMDCWire::localId(), m_pmgnIMF, msgSvc(), pi, TMLink::positionOnTrack(), TMLink::positionOnWire(), TMDCWireHit::reccdc(), s, TMLink::setDriftTime(), sin(), t(), tanl(), MdcRec_wirhit::tdc, tof(), Bes_Common::WARNING, TMDCWireHit::wire(), WireHitLeft, and WireHitRight.

Referenced by main().

                                        {
  //be cautious here
    const TMDCWireHit & h = * l.hit();
    const HepPoint3D & onTrack = l.positionOnTrack();
    const HepPoint3D & onWire = l.positionOnWire();
    unsigned leftRight = WireHitRight;
    //yzhang delete TEMP
    if (onWire.cross(onTrack).z() < 0.) leftRight = WireHitLeft;
    int charge = (t.helix().a()[2]>0) ? 1: -1;//track charge//yzhang 2012-05-03 TEMP
    //m_pmgnIMF NULL set mutable 2012-05-04 
    if(NULL == m_pmgnIMF) {
      Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
      if(NULL == m_pmgnIMF) {
        std::cout<< __FILE__<<" "<<__LINE__<<" Unable to open Magnetic field service"<<std::endl;
      }
    }
    const double Bz = -1000*m_pmgnIMF->getReferField();
#ifdef TRKRECO_DEBUG
    std::cout << " orignal ambig "<<leftRight<<" charge "<< charge <<" Bz "<<Bz<<std::endl;
#endif
    //if(charge * Bz <0){
    //  leftRight = (onWire.cross(onTrack).z() < 0.)
    //    ? WireHitLeft : WireHitRight;
    //}else{
    //  leftRight = (onWire.cross(onTrack).z() < 0.)
    //    ? WireHitRight : WireHitLeft;
    //}
#ifdef TRKRECO_DEBUG
    std::cout << " new ambig "<<leftRight<<std::endl;
#endif
    //zhangy

    //...No correction...
    if ((t0Offset == 0.) && (! _propagation) && (! _tof)) {
        distance = h.drift(leftRight);
        err = h.dDrift(leftRight);
        return;
    }

//    float x[3]={ onWire.x(), onWire.y(), onWire.z()};
//    float Tcenter = 0;

    //...TOF correction...
    float tof = 0.;
    if (_tof) {
        int imass = 3;
        float tl = t.helix().a()[4];
        float f = sqrt(1. + tl * tl);
        float s = fabs(t.helix().curv()) * fabs(l.dPhi()) * f;
        float p = f / fabs(t.helix().a()[2]);
//zsl   calcdc_tof2_(& imass, & p, & s, & tof);
        float Mass[5] = {0.000511,0.105,0.140,0.493677,0.98327};
        tof = s/29.98*sqrt(1.0+(Mass[imass-1]/p)*(Mass[imass-1]/p));

/*      float x[3]={ onWire.x(), onWire.y(), onWire.z()};
      float Rad = 81; // cm
      float dRho = t.helix().a()[0];
      float Phi0 = t.helix().a()[1];
      float a2 = (dRho*cos(Phi0))*(dRho*cos(Phi0));
      float b2 = (dRho*sin(Phi0))*(dRho*sin(Phi0));
      float Lproj = sqrt(Rad*Rad - a2 - b2);
      Tcenter = Lproj/29.98; //approxiate... cm/ns
      tof = s/29.98;  //cosmic
      if (x[1]>0) tof = -tof;
      */
    }

    //...T0 and propagation corrections...
    int wire = h.wire()->localId();
    int layerId = h.wire()->layerId();
    int side = leftRight;
//    if (side == 0) side = -1;
//    HepVector3D tp = t.helix().momentum(l.dPhi());
//    float p[3] = {tp.x(), tp.y(), tp.z()};
//    float x[3] = {onWire.x(), onWire.y(), onWire.z()};
//    float time = h.reccdc()->tdc + t0Offset - tof;
//    float dist;
//    float edist;
    double dist;
    double edist;
    int prop = _propagation;

    const double x0 = t.helix().pivot().x();
    const double y0 = t.helix().pivot().y();
    Hep3Vector pivot_helix(x0,y0,0);
    const double dr    =  t.helix().a()[0];
    const double phi0  =  t.helix().a()[1];
    const double kappa =  t.helix().a()[2];
    const double dz    =  t.helix().a()[3];
    const double tanl  =  t.helix().a()[4];

    //yzhang 2012-05-03 delete
    //const double alpha = 333.564095;
    //yzhang add
    const double alpha= 333.564095/(-1000*(m_pmgnIMF->getReferField()));
    //zhangy

    const double cox = x0 + dr*cos(phi0) + alpha*cos(phi0)/kappa;
    const double coy = y0 + dr*sin(phi0) + alpha*sin(phi0)/kappa;


    unsigned nHits = t.links().length();
    unsigned nStereos = 0;
    unsigned firstLyr = 44;
    unsigned lastLyr = 0;


    HepPoint3D ontrack = l.positionOnTrack();
    HepPoint3D onwire = l.positionOnWire();
    HepPoint3D dir(ontrack.y()-coy,cox-ontrack.x(),0);           
    double pos_phi=onwire.phi();
    double dir_phi=dir.phi();
    while(pos_phi>pi){pos_phi-=pi;}
    while(pos_phi<0){pos_phi+=pi;}
    while(dir_phi>pi){dir_phi-=pi;}
    while(dir_phi<0){dir_phi+=pi;}
    double entrangle=dir_phi-pos_phi;
    while(entrangle>pi/2)entrangle-=pi;
    while(entrangle<(-1)*pi/2)entrangle+=pi; 
    double zhit = onWire.z();
    
#ifdef TRKRECO_DEBUG 
    std::cout<<" onWire: "<<onWire<<std::endl;
    std::cout<<" zhit: "<<zhit<<std::endl;
#endif

    const double vinner = 220.0e8; // unit is cm/s
    const double vouter = 240.0e8; 
    double vprop = 300.0e9;
//    double tprop = 0.;
    
    
    if(layerId<8) {
      vprop = vinner;
    }
    else {
      vprop = vouter;
    }
  
    
    double EsT0 = -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();
      EsT0 = (*iter_evt)->getTest();
    }else{
      log << MSG::WARNING << "Could not find RecEsTimeCol" << endreq;
    }
   
    double rawTime = 0.;
    rawTime = h.reccdc()->tdc;
    double rawadc = 0.;
    rawadc =h.reccdc()->adc;
    //    double timewalk = CalibFunSvc_->getTimeWalk(layerid, 0.0);
    //    double timeDrift = rawTime - tof - tprop - EsTo - toWalk; 
    

    //----cal drift----//
    IMdcCalibFunSvc* l_mdcCalFunSvc;
    Gaudi::svcLocator() -> service("MdcCalibFunSvc", l_mdcCalFunSvc);

    double tprop = l_mdcCalFunSvc->getTprop(layerId,zhit*10.);

    double timeWalk = l_mdcCalFunSvc->getTimeWalk(layerId, rawadc);

    double T0 = l_mdcCalFunSvc->getT0(layerId,wire);
    double drifttime = rawTime - tof - tprop - timeWalk -T0;
    l.setDriftTime(drifttime);

#ifdef TRKRECO_DEBUG
    std::cout<<"timewalk is : "<< timeWalk << std::endl ;
    std::cout<<"T0 is : "<< T0 << std::endl ;
    std::cout<<"EsT0 is : "<< EsT0 << std::endl ;
    std::cout<<"tprop is : "<< tprop << std::endl ;
    std::cout<<"tof is : "<< tof << std::endl ;
    std::cout<<"rawTime is : "<< rawTime << std::endl ;
    std::cout<<"driftTime is : "<< drifttime << std::endl ;
#endif
//    dist = l_mdcCalFunSvc->rawTimeNoTOFToDist(h.reccdc()->tdc - tof - tprop-timeWalk, layerId, wire, side, entrangle); //default entranceangle is 0
    dist = l_mdcCalFunSvc->driftTimeToDist(drifttime,layerId, wire, side, entrangle);//by liucy 2010/05/12
    edist = l_mdcCalFunSvc->getSigma(layerId, side, dist, entrangle,tanl,zhit*10,rawadc);
    dist = dist/10.0; //mm->cm
    edist = edist/10.0;

//zsl    calcdc_driftdist_(& prop,
//                    & wire,
//                    & side,
//                    p,
//                    x,
//                    & time,
//                    & dist,
//                    & edist);

//      dist = (h.reccdc()->tdc-tof) * 40./10000;

//    distance = (double) dist;
//    err = (double) edist;
    distance = dist;
    err = edist;

    return;
}

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

dumps debug information.

Reimplemented from TMFitter.

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

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

Referenced by main().

int THelixFitter::fit	(	TTrackBase &	,
		float &	tev,
		float &	tev_err,
		double *	pre_chi2 = NULL,
		double *	fitted_chi2 = NULL
	)			const [inline]

Definition at line 267 of file THelixFitter.h.

References TTrackBase::_fitted, and main().

                                                               {
    a._fitted = false;
    return main(a, tev, tev_err, pre_chi2, fitted_chi2);
}

int THelixFitter::fit	(	TTrackBase &	,
		float	t0Offset,
		double *	pre_chi2 = NULL,
		double *	fitted_chi2 = NULL
	)			const [inline]

Definition at line 254 of file THelixFitter.h.

References TTrackBase::_fitted, _freeT0, and main().

                                                               {
    a._fitted = false;
    if (! _freeT0) return main(a, t0Offset, pre_chi2, fitted_chi2);
    else {
        float tev = t0Offset;
        float tevError;
        return main(a, tev, tevError, pre_chi2, fitted_chi2);
    }
}

int THelixFitter::fit	(	TTrackBase &	,
		double *	pre_chi2,
		double *	fitted_chi2
	)			const [inline]

Definition at line 242 of file THelixFitter.h.

References _freeT0, and main().

                                                               {
    if (! _freeT0) return main(a, 0., pre_chi2, fitted_chi2);
    else {
        float tev = 0.;
        float tevError;
        return main(a, tev, tevError, pre_chi2, fitted_chi2);
    }
}

int THelixFitter::fit

(

TTrackBase &

)

const [inline, virtual]

Implements TMFitter.

Definition at line 231 of file THelixFitter.h.

References _freeT0, and main().

Referenced by TBuilder::build(), TBuilder::buildRphi(), TBuilderCurl::buildStereo(), TBuilder::buildStereo(), TBuilderCurl::buildStereoMC(), TPerfectFinder::doit(), TTrackManager::mask(), TTrackManager::maskCurl(), TTrackManager::maskNormal(), TTrackManager::merge(), TTrackManager::refit(), TBuilder::salvage(), TCurlFinder::salvage3DTrack(), TTrackManager::T0(), and TTrackManager::T0Fit().

                                      {
    if (! _freeT0) return main(a, 0.);
    else {
        float tev = 0.;
        float tevError;
        return main(a, tev, tevError);
    }
}

bool THelixFitter::fit2D

(

bool

)

[inline]

Definition at line 175 of file THelixFitter.h.

References _fit2D.

                          {
    return _fit2D = a;
}

bool THelixFitter::fit2D

(

void

)

const [inline]

sets/returns 2D flag.

Definition at line 169 of file THelixFitter.h.

References _fit2D.

Referenced by TTrackManager::determineT0().

                              {
    return _fit2D;
}

void TMFitter::fitDone

(

TTrackBase &

)

const [protected, inherited]

sets the fitted flag. (Bad implementation)

Definition at line 24 of file TMFitter.cxx.

References t().

Referenced by TRobustLineFitter::fit(), TLineFitter::fit(), and TCircleFitter::fit().

                                      {
    t._fitted = true;
}

bool THelixFitter::freeT0

(

bool

)

[inline]

Definition at line 281 of file THelixFitter.h.

References _freeT0.

                           {
    return _freeT0 = a;
}

bool THelixFitter::freeT0

(

void

)

const [inline]

sets/returns free T0 flag.

Definition at line 275 of file THelixFitter.h.

References _freeT0.

Referenced by TTrackManager::fittingFlag(), and TBuilder::TBuilder().

                               {
    return _freeT0;
}

IMagneticFieldSvc* THelixFitter::getMagneticFieldPointer

(

void

)

const [inline]

Definition at line 98 of file THelixFitter.h.

References m_pmgnIMF.

Referenced by TBuilder::build().

{return m_pmgnIMF;}

int THelixFitter::main	(	TTrackBase &	,
		float &	tev,
		float &	tev_err,
		double *	pre_chi2 = NULL,
		double *	fitted_chi2 = NULL
	)			const [private]

main routine for free T0.

Definition at line 1471 of file THelixFitter.cxx.

References _fit2D, _fitted_chi2, _pre_chi2, _sag, AxialHits(), DBL_MAX, Helix::direction(), drift(), dxda(), ganga-rec::j, NStereoHits(), NTrailMax, TTrackBase::objectType(), t(), TFitAlreadyFitted, TFitErrorFewHits, TFitFailed, TFitUnavailable, Track, WireHitLeft, and WireHitRight.

                                                                {
//=====================================================================
    //...Initialize
    _pre_chi2 = _fitted_chi2 = 0.;
    if(pre_chi2)*pre_chi2 = _pre_chi2;
    if(fitted_chi2)*fitted_chi2 = _fitted_chi2;

    //...Type check...
    if (b.objectType() != Track) return TFitUnavailable;
    TTrack & t = (TTrack &) b;

    //...Already fitted ?...
    if (t.fitted()) return TFitAlreadyFitted;

    //...Count # of hits...
    AList<TMLink> cores = t.cores();
    if (_fit2D) cores = AxialHits(cores);
    unsigned nCores = cores.length();
    unsigned nStereoCores = NStereoHits(cores);

    //...2D or 3D...
    bool fitBy2D = _fit2D;
    if (! fitBy2D) fitBy2D = (! nStereoCores);

    //...Check # of hits...
    if (! fitBy2D) {
        if ((nStereoCores < 2) || (nCores - nStereoCores < 3))
            return TFitErrorFewHits;
    }
    else {
        if (nCores < 3) return TFitErrorFewHits;
    }

    //...Setup...
    Vector a(6), da(6);
    Vector a_5dim(5);
    for (unsigned j = 0; j < 5; j++) a[j] = t.helix().a()[j];
    a[5] = tev;
    Vector dxda(5);
    Vector dyda(5);
    Vector dzda(5);
    Vector dDda(6);
    Vector dDda_5dim(5);
    Vector dchi2da(6);
    SymMatrix d2chi2d2a(6, 0);
    static const SymMatrix zero6(6, 0);
    double chi2;
    double chi2Old = DBL_MAX;
//    const double convergence = Convergence;
    const double convergence = 1.0e-4;   //Liuqg
    bool allAxial = true;
    Matrix e(3, 3);
    Vector f(3);
    int err = 0;
    double factor = 1.0;//jtanaka0715

    //...Fitting loop...
    unsigned nTrial = 0;
    while (nTrial < NTrailMax) {

        //...Set up...
        chi2 = 0.;
        for (unsigned j = 0; j < 6; j++) dchi2da[j] = 0.;
        d2chi2d2a = zero6;

        //...Loop with hits...
        unsigned i = 0;
        while (TMLink * l = cores[i++]) {
            const TMDCWireHit & h = * l->hit();

            //...Cal. closest points...
            t.approach(* l, _sag);
            double dPhi = l->dPhi();
            const HepPoint3D & onTrack = l->positionOnTrack();
            const HepPoint3D & onWire = l->positionOnWire();
            unsigned leftRight = (onWire.cross(onTrack).z() < 0.) ? WireHitLeft : WireHitRight;

            //...Obtain drift distance and its error...
            double distance;
            double eDistance;
            double dddt;
            drift(t, * l, tev, distance, eDistance, dddt);
//          l->drift(distance,0);
//          l->drift(distance,1);
//          l->dDrift(distance,0);
//          l->dDrift(distance,1);
           

            double inv_eDistance2 = 1./(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);

            //...Chi2 related...
            // dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag;
            double vw[3] = { h.wire()->direction().x(),
                             h.wire()->direction().y(),
                             h.wire()->direction().z() };
            double vwxy = vw[0]*vw[1];
            double vwyz = vw[1]*vw[2];
            double vwzx = vw[2]*vw[0];
            dDda_5dim = (vmag > 0.)
                ? ((v.x() * (1. - vw[0] * vw[0]) -
                    v.y() * vwxy - v.z() * vwzx)
                   * dxda + 
                   (v.y() * (1. - vw[1] * vw[1]) -
                    v.z() * vwyz - v.x() * vwxy)
                   * dyda + 
                   (v.z() * (1. - vw[2] * vw[2]) -
                    v.x() * vwzx - v.y() * vwyz)
                   * dzda) / vmag
                : Vector(5,0);
            if (vmag <= 0.0) {
                std::cout << "    in fit " << onTrack << ", " << onWire;
                h.dump();
            }
            //      for (unsigned j = 0; j < 5; j++) dDda[j] = dDda_5dim[j];
            dDda[0] = dDda_5dim[0];
            dDda[1] = dDda_5dim[1];
            dDda[2] = dDda_5dim[2];
            dDda[3] = dDda_5dim[3];
            dDda[4] = dDda_5dim[4];
            dDda[5] = -dddt;

            dchi2da += (dDistance * inv_eDistance2) * dDda;
            d2chi2d2a += vT_times_v(dDda) * inv_eDistance2;
            double pChi2 = dDistance * dDistance * inv_eDistance2;
            chi2 += pChi2;

            //...Store results...
            l->update(onTrack, onWire, leftRight, pChi2);
        }

        //...Save chi2 information...
        if(nTrial == 0){
          _pre_chi2 = chi2;
          _fitted_chi2 = chi2;
        }else _fitted_chi2 = chi2;

        //...Check condition...
        double change = chi2Old - chi2;
        if (fabs(change) < convergence) break;
        //temp
                factor = 1.0;
        //temp
        if (change < 0.) {
//          a += factor * da;
//          t._helix->a(a);
//          break;
            factor = 0.5;
        }

        chi2Old = chi2;

        //...Cal. helix parameters for next loop...
        if (fitBy2D) {
            f = dchi2da.sub(1, 4);
            e = d2chi2d2a.sub(1, 4);
            f = solve(e, f);
            da[0] = f[0];
            da[1] = f[1];
            da[2] = f[2];
            da[3] = f[3];
            da[4] = 0.;
            da[5] = 0.;
        }
        else {
            da = solve(d2chi2d2a, dchi2da);
        }

        a -= factor * da;

        //      for (unsigned j = 0; j < 5; j++) a_5dim[j] = a[j];
        a_5dim[0] = a[0];
        a_5dim[1] = a[1];
        a_5dim[2] = a[2];
        a_5dim[3] = a[3];
        a_5dim[4] = a[4];
        t._helix->a(a_5dim);
        tev = a[5];
        //temp
        //      if(nTrial == 0) std::cout << "initial chi2=" <<chi2 << std::endl;
        //temp
        ++nTrial;

#ifdef TRKRECO_DEBUG
        std::cout << "fit " << nTrial-1<< " : " << chi2 << " : " << change << std::endl;
#endif
    }

    //...Cal. error matrix...
    SymMatrix Ea(6, 0);
    unsigned dim;
    if (fitBy2D) {
        dim = 4;
        SymMatrix Eb(4, 0), Ec(4, 0);
        Eb = d2chi2d2a.sub(1, 4);
        Ec = Eb.inverse(err);
        Ea[0][0] = Ec[0][0];
        Ea[0][1] = Ec[0][1];
        Ea[0][2] = Ec[0][2];
        Ea[0][3] = Ec[0][3];
        Ea[1][1] = Ec[1][1];
        Ea[1][2] = Ec[1][2];
        Ea[1][3] = Ec[1][3];
        Ea[2][2] = Ec[2][2];
        Ea[2][3] = Ec[2][3];
        Ea[3][3] = Ec[3][3];
    }
    else {
        dim = 6;
        Ea = d2chi2d2a.inverse(err);
        //        std::cout << "err flg=" << err << std::endl;
    }


    //...Store information...
    if (! err) {
        for (unsigned j = 0; j < 5; j++) a_5dim[j] = a[j];
        SymMatrix Ea_5dim(5, 0);
        Ea_5dim = Ea.sub(1, 5);
        t._helix->a(a_5dim);
        t._helix->Ea(Ea_5dim);
        tev = a[5];
        tev_err = sqrt(Ea[5][5]);
        //temp
        //      std::cout << "nTrial=" << nTrial << std::endl;
        //      std::cout << "chi2="   << chi2   << std::endl;
        //      std::cout << "pt:" << fabs(1./a_5dim[2])<< std::endl;
        //      std::cout << "tev,tev_err="<<tev<<" "<<tev_err<<std::endl;
        //temp

        t._fitted = true;
    }
    else {
        err = TFitFailed;
    }

    t._ndf = nCores - dim;
    t._chi2 = chi2;

    if(pre_chi2)*pre_chi2 = _pre_chi2;
    if(fitted_chi2)*fitted_chi2 = _fitted_chi2;

    return err;
}

int THelixFitter::main	(	TTrackBase &	,
		float	t0Offset,
		double *	pre_chi2 = NULL,
		double *	fitted_chi2 = NULL
	)			const [private]

main routine for fixed T0.

Definition at line 150 of file THelixFitter.cxx.

References _fit2D, _fitted_chi2, _pre_chi2, _sag, AxialHits(), chi2(), Convergence, DBL_MAX, TMDCWire::direction(), drift(), TMDCWireHit::dump(), dxda(), showlog::err, genRecEmupikp::i, ganga-rec::j, NStereoHits(), NTrailMax, TTrackBase::objectType(), TMDCWireHit::state(), t(), TFitAlreadyFitted, TFitErrorFewHits, TFitFailed, TFitUnavailable, Track, TrkRecoHelixFitterChisqMax, v, TMDCWireHit::wire(), WireHitLeft, WireHitPatternLeft, WireHitPatternRight, and WireHitRight.

Referenced by fit().

                                                                {
#ifdef TRKRECO_DEBUG
/*
        if (first) {
        first = false;
        extern BelleTupleManager * BASF_Histogram;
        BelleTupleManager * m = BASF_Histogram;
        _nCall[0] = m->histogram("HF nCall all", 1, 0., 1.);
        _nCall[1] = m->histogram("HF nCall conf f2d l0", 1, 0., 1.);
        _nCall[2] = m->histogram("HF nCall conf f3d l0", 1, 0., 1.);
        _nCall[3] = m->histogram("HF nCall conf f2d l1", 1, 0., 1.);
        _nCall[4] = m->histogram("HF nCall conf f3d l1", 1, 0., 1.);
        _nCall[5] = m->histogram("HF nCall conf s2d", 1, 0., 1.);
        _nCall[6] = m->histogram("HF nCall conf s3d", 1, 0., 1.);
        _nCall[7] = m->histogram("HF nCall other", 1, 0., 1.);
        _nTrial[0] = m->histogram("HF nTrial all", 100, 0., 100.);
        _nTrial[1] = m->histogram("HF nTrial conf f2d l0", 100, 0., 100.);
        _nTrial[2] = m->histogram("HF nTrial conf f3d l0", 100, 0., 100.);
        _nTrial[3] = m->histogram("HF nTrial conf f2d l1", 100, 0., 100.);
        _nTrial[4] = m->histogram("HF nTrial conf f3d l1", 100, 0., 100.);
        _nTrial[5] = m->histogram("HF nTrial conf s2d", 100, 0., 100.);
        _nTrial[6] = m->histogram("HF nTrial conf s3d", 100, 0., 100.);
        _nTrial[7] = m->histogram("HF nTrial other", 100, 0., 100.);
        _pull[0][0][0] = m->histogram("HF pull ax true all",
                                      100, 0., 5000.);
        _pull[0][0][2] = m->histogram("HF pull ax true conf f3d l0",
                                      100, 0., 5000.);
        _pull[0][0][4] = m->histogram("HF pull ax true conf f3d l1",
                                      100, 0., 5000.);
        _pull[0][0][6] = m->histogram("HF pull ax true conf s3d",
                                      100, 0., 5000.);
        _pull[0][0][7] = m->histogram("HF pull ax true other",
                                      100, 0., 5000.);
        _pull[1][0][0] = m->histogram("HF pull st true all",
                                      100, 0., 5000.);
        _pull[1][0][2] = m->histogram("HF pull st true conf f3d l0", 
                                      100, 0., 5000.);
        _pull[1][0][4] = m->histogram("HF pull st true conf f3d l1",
                                      100, 0., 5000.);
        _pull[1][0][6] = m->histogram("HF pull st true conf s3d",
                                      100, 0., 5000.);
        _pull[1][0][7] = m->histogram("HF pull st true other",
                                      100, 0., 5000.);
        _pull[0][1][0] = m->histogram("HF pull ax wrong all",
                                      100, 0., 5000.);
        _pull[0][1][2] = m->histogram("HF pull ax wrong conf f3d l0",
                                      100, 0., 5000.);
        _pull[0][1][4] = m->histogram("HF pull ax wrong conf f3d l1",
                                      100, 0., 5000.);
        _pull[0][1][6] = m->histogram("HF pull ax wrong conf s3d",
                                      100, 0., 5000.);
        _pull[0][1][7] = m->histogram("HF pull ax wrong other",
                                      100, 0., 5000.);
        _pull[1][1][0] = m->histogram("HF pull st wrong all",
                                      100, 0., 5000.);
        _pull[1][1][2] = m->histogram("HF pull st wrong conf f3d l0",
                                      100, 0., 5000.);
        _pull[1][1][4] = m->histogram("HF pull st wrong conf f3d l1",
                                      100, 0., 5000.);
        _pull[1][1][6] = m->histogram("HF pull st wrong conf s3d",
                                      100, 0., 5000.);
        _pull[1][1][7] = m->histogram("HF pull st wrong other",
                                      100, 0., 5000.);
        _nTrialPositive = m->histogram("HF nTrial +", 100, 0., 100.);
        _nTrialNegative = m->histogram("HF nTrial -", 100, 0., 100.);
    }
    _nCall[0]->accumulate(.5);
    if (TConformalFinder::_stage == ConformalOutside)
        _nCall[7]->accumulate(.5);
    else if (TConformalFinder::_stage == ConformalFast2DLevel0)
        _nCall[1]->accumulate(.5);
    else if (TConformalFinder::_stage == ConformalFast3DLevel0)
        _nCall[2]->accumulate(.5);
    else if (TConformalFinder::_stage == ConformalFast2DLevel1)
        _nCall[3]->accumulate(.5);
    else if (TConformalFinder::_stage == ConformalFast3DLevel1)
        _nCall[4]->accumulate(.5);
    else if (TConformalFinder::_stage == ConformalSlow2D)
        _nCall[5]->accumulate(.5);
    else if (TConformalFinder::_stage == ConformalSlow3D)
        _nCall[6]->accumulate(.5);
    bool posi = true;
    const TTrackHEP & hep = Links2HEP(b.links());
*/
#endif
    //...Initialize
    _pre_chi2 = _fitted_chi2 = 0.;
    if(pre_chi2)*pre_chi2 = _pre_chi2;
    if(fitted_chi2)*fitted_chi2 = _fitted_chi2;

    //...Type check...
    if (b.objectType() != Track) return TFitUnavailable;
    TTrack & t = (TTrack &) b;

    //...Already fitted ?...
    if (t.fitted()) return TFitAlreadyFitted;

    //...Count # of hits...
    AList<TMLink> cores = t.cores();
#ifdef TRKRECO_DEBUG
    cout<<__FILE__<<" cores in helix = "<<cores.length()<<endl;
#endif
    if (_fit2D) cores = AxialHits(cores);
    unsigned nCores = cores.length();
    unsigned nStereoCores = NStereoHits(cores);

    //...2D or 3D...
    bool fitBy2D = _fit2D;
    if (! fitBy2D) fitBy2D = (! nStereoCores);

    //...Check # of hits...
    if (! fitBy2D) {
        if ((nStereoCores < 2) || (nCores - nStereoCores < 3))
            return TFitErrorFewHits;
    }
    else {
        if (nCores < 3) return TFitErrorFewHits;
    }

    //...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);
    static const SymMatrix zero5(5, 0);
    double chi2;
    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

    //...For bad hit rejection...(by JT, 2001/04/12)...
    int flagBad = 0;
    if (TrkRecoHelixFitterChisqMax != 0.)
        flagBad = 1;
    AList<TMLink> initBadWires; 
    unsigned nInitBadWires = 0;
    Vector initBadDchi2da(5);
    SymMatrix initBadD2chi2d2a(5, 0);
    for (unsigned j = 0; j < 5; ++j) initBadDchi2da[j] = 0.;
    double initBadChi2 = 0.;

    //...Fitting loop...
    unsigned nTrial = 0;
    while (nTrial < NTrailMax) {

        //...Set up...
        chi2 = 0.;
        for (unsigned j = 0; j < 5; j++) dchi2da[j] = 0.;
        d2chi2d2a = zero5;

    //yuany
#ifdef TRKRECO_DEBUG
        cout<<"helix fitter a5 "<<t.helix().a()<< " cores.length "<<cores.length()<<endl;
#endif
        //...Loop with hits...
        unsigned i = 0;
        while (TMLink * l = cores[i++]) {
            const TMDCWireHit & h = * l->hit();
    //yuany
#ifdef TRKRECO_DEBUG
            cout<<"trial "<<nTrial<<" wire name "<<l->wire()->name()<<"    L "<<(h.state()&WireHitPatternLeft)<<" R "<<(h.state()&WireHitPatternRight)<<" link "<<l->leftRight()<<endl;
#endif
            //...Cal. closest points...
            t.approach(* l, _sag);
            double dPhi = l->dPhi();
            const HepPoint3D & onTrack = l->positionOnTrack();
            const HepPoint3D & onWire = l->positionOnWire();
            unsigned leftRight = (onWire.cross(onTrack).z() < 0.)
                ? WireHitLeft : WireHitRight;

            //...Obtain drift distance and its error...
            double distance;
            double eDistance;
            drift(t, * l, t0Offset, distance, eDistance);
            l->drift(distance,0);
            l->drift(distance,1);
            l->dDrift(eDistance,0);
            l->dDrift(eDistance,1);


            double inv_eDistance2 = 1. / (eDistance * eDistance);
            
            //...Residual...
            HepVector3D v = onTrack - onWire;
            double vmag = v.mag();
#ifdef TRKRECO_DEBUG
            std::cout<<"THelixFitter::eDistance-----"<<eDistance<<endl;
            cout<<"  vmag = "<<vmag<<"   distance = "<<distance<<"  eDistance = "<<eDistance<<endl;
#endif
            double dDistance = vmag - distance;
            
            //...dxda...
            this->dxda(* l, t.helix(), dPhi, dxda, dyda, dzda);

            //...Chi2 related...
            // dDda = (v.x() * dxda + v.y() * dyda + v.z() * dzda) / vmag;
            // Vector3 vw = h.wire()->direction();
            double vw[3] = { h.wire()->direction().x(),
                             h.wire()->direction().y(),
                             h.wire()->direction().z() };
            double vwxy = vw[0]*vw[1];
            double vwyz = vw[1]*vw[2];
            double vwzx = vw[2]*vw[0];
            dDda = (vmag > 0.)
                ? ((v.x() * (1. - vw[0] * vw[0]) -
                    v.y() * vwxy - v.z() * vwzx)
                   * dxda + 
                   (v.y() * (1. - vw[1] * vw[1]) -
                    v.z() * vwyz - v.x() * vwxy)
                   * dyda + 
                   (v.z() * (1. - vw[2] * vw[2]) -
                    v.x() * vwzx - v.y() * vwyz)
                   * dzda) / vmag
                : Vector(5,0);
            if (vmag <= 0.0) {
                std::cout << "    in fit " << onTrack << ", " << onWire;
                h.dump();
            }
            double pChi2 = dDistance * dDistance * inv_eDistance2;
#ifdef TRKRECO_DEBUG
            std::cout<<"THelixFitter::dDistance"<<dDistance<<" inv_eDistance2 "<<inv_eDistance2<<endl; 
            cout<<"Liuqg check ... .. pChi2 = "<<pChi2<<endl;
#endif

            //...Bad hit rejection...
            if (flagBad && nTrial == 0) {
                if (pChi2 > TrkRecoHelixFitterChisqMax) {
                    initBadWires.append(l);
                    initBadDchi2da += (dDistance * inv_eDistance2) * dDda;
                    initBadD2chi2d2a += vT_times_v(dDda) * inv_eDistance2;
                    initBadChi2 += pChi2;
                }
            }
            else {
                dchi2da += (dDistance * inv_eDistance2) * dDda;
                d2chi2d2a += vT_times_v(dDda) * inv_eDistance2;
                chi2 += pChi2;      

                //...Store results...
                l->update(onTrack, onWire, leftRight, pChi2);
            }

#ifdef TRKRECO_DEBUG
//          if ((! fitBy2D) && (nTrial == 0)) {
//              unsigned as = 0;
//              if (l->hit()->wire()->stereo()) as = 1;
//              unsigned mt = 0;
//              if (& hep != l->hit()->mc()->hep()) mt = 1;

//              _pull[as][mt][0]->accumulate(pChi2);
//              if (TConformalFinder::_stage == ConformalOutside)
//                  _pull[as][mt][7]->accumulate(pChi2);
//              else if (TConformalFinder::_stage == ConformalFast3DLevel0)
//                  _pull[as][mt][2]->accumulate(pChi2);
//              else if (TConformalFinder::_stage == ConformalFast3DLevel1)
//                  _pull[as][mt][4]->accumulate(pChi2);
//              else if (TConformalFinder::_stage == ConformalSlow3D)
//                  _pull[as][mt][6]->accumulate(pChi2);
//          }
#endif
        }

        //...Bad hit rejection...
        if (flagBad && nTrial == 0) {
            if ((initBadWires.length() == 1 || initBadWires.length() == 2) &&
                nCores >= 20 &&
                chi2 / (double)(nCores - initBadWires.length()) < 10.) {
                cores.remove(initBadWires);
                nInitBadWires = initBadWires.length();
            }
            else if (initBadWires.length() != 0) {
                dchi2da += initBadDchi2da;
                d2chi2d2a += initBadD2chi2d2a;
                chi2 += initBadChi2;
            }
        }

        //...Save chi2 information...
        if (nTrial == 0) {
            _pre_chi2 = chi2;
            _fitted_chi2 = chi2;
        }
        else {
            _fitted_chi2 = chi2;
        }

        //...Check condition...
        double change = chi2Old - chi2;
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout<<" chi2 change  "<<change <<" convergence  "<<convergence <<" break?  "<<(fabs(change) < convergence == true)<<std::endl;
#endif
        if (fabs(change) < convergence) break;
        if (change < 0.) {
//          a += factor * da;
//          t._helix->a(a);
//          break;
            factor = 0.5;
        }
        chi2Old = chi2;

        //...Cal. helix parameters for next loop...
        if (fitBy2D) {
            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);
        }

        a -= factor * da;
        t._helix->a(a);
        ++nTrial;

        // jtanaka 001008
        //if( fabs(a[3]) > 200. ){
        // yiwasaki 001010
        if( fabs(a[3]) > 1000. ){
          // stop "fit" and return error.
          //std::cout << "Stop Fit... " << a << std::endl;
          err = 1;
          break;
        }
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "            fit " << nTrial-1<< " : " << chi2 << " : "
                  << change << std::endl;
#endif
    }

#ifdef TRKRECO_DEBUG
/*
    _nTrial[0]->accumulate(float(nTrial) + .5);
    if (TConformalFinder::_stage == ConformalOutside)
        _nTrial[7]->accumulate(float(nTrial) + .5);
    else if (TConformalFinder::_stage == ConformalFast2DLevel0)
        _nTrial[1]->accumulate(float(nTrial) + .5);
    else if (TConformalFinder::_stage == ConformalFast3DLevel0)
        _nTrial[2]->accumulate(float(nTrial) + .5);
    else if (TConformalFinder::_stage == ConformalFast2DLevel1)
        _nTrial[3]->accumulate(float(nTrial) + .5);
    else if (TConformalFinder::_stage == ConformalFast3DLevel1)
        _nTrial[4]->accumulate(float(nTrial) + .5);
    else if (TConformalFinder::_stage == ConformalSlow2D)
        _nTrial[5]->accumulate(float(nTrial) + .5);
    else if (TConformalFinder::_stage == ConformalSlow3D)
        _nTrial[6]->accumulate(float(nTrial) + .5);

    if (posi) _nTrialPositive->accumulate((float) nTrial + .5);
    else _nTrialNegative->accumulate((float) nTrial + .5);
*/
#endif

    //...Cal. error matrix...
    SymMatrix Ea(5, 0);
    unsigned dim;
    if (fitBy2D) {
        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);
        t._fitted = true;
    }
    else {
        err = TFitFailed;
    }

    t._charge = copysign(1., a[2]);    
    t._ndf = nCores - dim -nInitBadWires;
    t._chi2 = chi2;

    //...Treatment for bad wires...
    if (nInitBadWires) {
        for  (unsigned i = 0; i < nInitBadWires; i++) {
            TMLink * l = initBadWires[i];
            t.approach(* l, _sag);
            double dPhi = l->dPhi();
            const HepPoint3D & onTrack = l->positionOnTrack();
            const HepPoint3D & onWire = l->positionOnWire();
            HepVector3D v = onTrack - onWire;
            double vmag = v.mag();
            unsigned leftRight = (onWire.cross(onTrack).z() < 0.)
                ? WireHitLeft : WireHitRight;
            double distance;
            double eDistance;
            drift(t, * l, t0Offset, distance, eDistance);
            l->drift(distance,0);
            l->drift(distance,1);
            l->dDrift(eDistance,0);
            l->dDrift(eDistance,1);


            double inv_eDistance2 = 1. / (eDistance * eDistance);
            double dDistance = vmag - distance;
            double pChi2 = dDistance * dDistance * inv_eDistance2;
            l->update(onTrack, onWire, leftRight, pChi2);
        }
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "            HelixFitter : # of rejected hits="
                  << nInitBadWires << endl;
#endif
    }

    if (pre_chi2) * pre_chi2 = _pre_chi2;
    if (fitted_chi2) * fitted_chi2 = _fitted_chi2;

    return err;
}

const std::string & TMFitter::name

(

void

)

const [inline, inherited]

returns name.

Definition at line 73 of file TMFitter.h.

References TMFitter::_name.

                         {
    return _name;
}

double THelixFitter::preChi2

(

void

)

const [inline]

returns sum of chi2 before fit.

Definition at line 299 of file THelixFitter.h.

References _pre_chi2.

                                {
    return _pre_chi2;
}

bool THelixFitter::propagation

(

bool

)

[inline]

Definition at line 199 of file THelixFitter.h.

References _propagation, and propagation().

                                {
    if (a) _propagation = 1;
    else   _propagation = 0;
    return propagation();
}

bool THelixFitter::propagation

(

void

)

const [inline]

sets/returns propagation-delay correction flag.

Definition at line 193 of file THelixFitter.h.

References _propagation.

Referenced by TBuilderCurl::buildStereo(), TTrackManager::determineT0(), TTrackManager::fittingFlag(), propagation(), TBuilderCurl::setParam(), TBuilder::TBuilder(), and TBuilderCurl::TBuilderCurl().

                                    {
    return (bool) _propagation;
}

bool THelixFitter::sag

(

bool

)

[inline]

Definition at line 187 of file THelixFitter.h.

References _sag.

                        {
    return _sag = a;
}

bool THelixFitter::sag

(

void

)

const [inline]

sets/returns sag correction flag.

Definition at line 181 of file THelixFitter.h.

References _sag.

Referenced by TBuilderCurl::buildStereo(), TTrackManager::fittingFlag(), TBuilderCurl::setParam(), TBuilder::TBuilder(), and TBuilderCurl::TBuilderCurl().

                            {
    return _sag;
}

bool THelixFitter::tanl

(

bool

)

[inline]

Definition at line 225 of file THelixFitter.h.

References _tanl.

                         {
    return _tanl = a;
}

bool THelixFitter::tanl

(

void

)

const [inline]

sets/returns tanLambda correction flag.

Definition at line 219 of file THelixFitter.h.

References _tanl.

Referenced by drift().

                             {
    return _tanl;
}

bool THelixFitter::tof

(

bool

)

[inline]

Definition at line 213 of file THelixFitter.h.

References _tof.

                        {
    return _tof = a;
}

bool THelixFitter::tof

(

void

)

const [inline]

sets/returns propagation-delay correction flag.

Definition at line 207 of file THelixFitter.h.

References _tof.

Referenced by TBuilderCurl::buildStereo(), TTrackManager::determineT0(), drift(), TTrackManager::fittingFlag(), TBuilderCurl::setParam(), TBuilder::TBuilder(), and TBuilderCurl::TBuilderCurl().

                            {
    return _tof;
}

---

Member Data Documentation

unsigned THelixFitter::_corrections [private]

Definition at line 139 of file THelixFitter.h.

Referenced by corrections().

double THelixFitter::_driftTime [private]

Definition at line 147 of file THelixFitter.h.

bool THelixFitter::_fit2D [private]

Definition at line 137 of file THelixFitter.h.

Referenced by fit2D(), and main().

double THelixFitter::_fitted_chi2 [mutable, private]

Definition at line 149 of file THelixFitter.h.

Referenced by chi2(), and main().

bool THelixFitter::_freeT0 [private]

Definition at line 138 of file THelixFitter.h.

Referenced by fit(), and freeT0().

double THelixFitter::_pre_chi2 [mutable, private]

Definition at line 148 of file THelixFitter.h.

Referenced by main(), and preChi2().

int THelixFitter::_propagation [private]

Definition at line 142 of file THelixFitter.h.

Referenced by drift(), and propagation().

bool THelixFitter::_sag [private]

Definition at line 141 of file THelixFitter.h.

Referenced by main(), and sag().

bool THelixFitter::_tanl [private]

Definition at line 144 of file THelixFitter.h.

Referenced by tanl().

bool THelixFitter::_tof [private]

Definition at line 143 of file THelixFitter.h.

Referenced by drift(), and tof().

IMagneticFieldSvc* THelixFitter::m_pmgnIMF [mutable, private]

Definition at line 151 of file THelixFitter.h.

Referenced by drift(), and getMagneticFieldPointer().

---