StRoot: StarRoot/THelixTrack.h Source File

00001 #ifndef THELIXTRACK_H
00002 #define THELIXTRACK_H
00003 #include "TObject.h"
00004 #include "TArrayD.h"
00005 #include "TPolinom.h"
00006 class TCEmx_t
00007 { 
00008 public:
00009 const double *Arr() const       { return &mHH;}
00010       double *Arr()     { return &mHH;}
00011 const double &operator[](int idx) const         { return (&mHH)[idx];}
00012       double &operator[](int idx)               { return (&mHH)[idx];}
00013       TCEmx_t &operator*=(double f)             
00014               { for (int i=0;i<6;i++) {Arr()[i]*=f;} return *this;}
00015 void Clear()                    { memset(this,0,sizeof(*this));}
00016      TCEmx_t()                  { Clear();}
00017 void Set(const double *err);    
00018 void Move(double const F[3][3]);
00019 double Sign() const;
00020 void Backward();
00021 public:
00022 
00023 double
00024 mHH,
00025 mHA, mAA,
00026 mHC, mAC, mCC;
00027 };
00028 
00029 class THEmx_t
00030 { 
00031 public:
00032      THEmx_t()                  { Clear();}
00033 const double *Arr() const       { return &mHH;}
00034       double *Arr()     { return &mHH;}
00035 const double &operator[](int idx) const         { return (&mHH)[idx];}
00036       double &operator[](int idx)               { return (&mHH)[idx];}
00037       THEmx_t &operator*=(double f)             
00038               { for (int i=0;i<15;i++) {Arr()[i]*=f;} return *this;}
00039 void Clear()                    { memset(this,0,sizeof(*this));}
00040 void Set(const double *err);
00041 void Set(const double *errxy,const double *errz);
00042 void Move(double const F[5][5]);
00043 void Backward();
00044 void Print(const char *tit=0) const;
00045 double Sign() const;
00046 public:
00047 //  dA: delta azimuth angle; dH: error along ort to dir and Z axis
00048 //  dC: error of curvature;  dZ == dZ; dL = dLambda
00049 double
00050 mHH,
00051 mHA, mAA,
00052 mHC, mAC, mCC,
00053 mHZ, mAZ, mCZ, mZZ,
00054 mHL, mAL, mCL, mZL, mLL;
00055 };
00056 
00057 
00058 class TCircle: public TObject
00059 {
00060 friend class THelixTrack;
00061 friend class TCircleFitter;
00062 public:
00063  TCircle();
00064  TCircle(const double *x,const double *dir,double rho);
00065  TCircle(const TCircle& fr);
00066  TCircle(const TCircle* fr);    //special ctr without errors
00067 ~TCircle();
00068 TCircle &operator=(const TCircle& fr);
00069  void Set(const double *x=0,const double *dir=0,const double rho=0);
00070 virtual void  Clear(const char *opt="");
00071 const double* Pos() const       {return fX;  } 
00072       double* Pos()             {return fX;  } 
00073 const double* Dir() const       {return fD;  } 
00074       double  Rho() const       {return fRho;}
00075       double& Rho()             {return fRho;}
00076       void    Nor(double *norVec) const; 
00077       void    SetEmx(const double *err=0);
00078 const TCEmx_t *Emx() const      {return fEmx;} 
00079       TCEmx_t *Emx()            {return fEmx;} 
00080   void GetCenter(double center[2]) const;
00081 double Path(const double pnt[2]) const;
00082 double Path(const double pnt[2], const double exy[3]) const;
00083 double Path(const TCircle &tc,double *s2=0) const;
00084 double Move(double step);
00085 void   Rot(double angle);
00086 void   Rot(double cosa,double sina);
00087 void   Backward();
00088 double Eval(double step,double *xy,double *dir=0) const;
00089 void   Show(int nPts,const double *Pts,int pstep=2);
00090 virtual void   Print(const char* chopt = "") const;
00091 void   SetStrait(int strait=1)          {SetBit(1,strait) ;}
00092 int    IsStrait()                       {return TestBit(1);}
00093 
00094 //      static funs
00095 static void Test2();
00096 static void Test3();
00097 static void Test4();
00098 static void TestMtx();
00099 
00100 private:
00101 void  MoveErrs(double l);
00102 void  MakeMtx (double l,double F[3][3]);
00103 
00104 protected:
00105 double fX[2];
00106 double fD[2];
00107 double fRho;
00108 TCEmx_t *fEmx; //let h = fX[1]*fD[0], a=atan2(fD[1],fD[0]),c=fRho
00109                 // hh,
00110                 // ah,aa,
00111                 // ch,ca,cc
00112 ClassDef(TCircle,0)
00113 };
00114 class TCircleFitterAux;
00115 class TCircleFitterAux
00116 {
00117   public:
00118   static int dSize() {return sizeof(TCircleFitterAux)/sizeof(double);}
00119   public:
00120   double x,y,z;         //x,y,z of measured point
00121   double exy[3];        //err matrix(xx,xy,yy) of x,y
00122   double ezz;           //error of z
00123   double wt;            //calculated weight
00124 
00125 };
00126 class TCircleFitter: public TCircle
00127 {
00128 public:
00129        TCircleFitter();
00130 int    Size() const                     {return fN;}
00131 int    Used() const                     {return fNuse;}
00132 void   Add (double x,double y,const double *errs=0); 
00133 void   Add (double x,double y,double z); 
00134 void   AddErr(const double *errs,double errz=0); 
00135 void   AddErr(double errh,double errz=0); 
00136 void   AddZ(double z,double err2z=0);
00137 double Fit();   
00138 void   MakeErrs();
00139 double FixAt(const double vals[5],int flag); 
00140 void   Skip(int idx); 
00141 double GetZ0() const                    {return fZ0    ;}
00142 double GetTanL() const                  {return fTanL  ;}
00143 void   SetCase(int kase=0)              {fCase=kase    ;}
00144 int    GetCase() const                  {return fKase  ;}
00145 double Chi2()   const                   {return fChi2  ;}
00146 int    Ndf()    const                   {return fNdf   ;}
00147 double Chi2Z () const                   {return fChi2Z ;}
00148 void   SetNdf(int ndf);                 
00149 double EvalChi2();
00150 void   Clear(const char *opt ="");
00151 void   Print(const char* chopt = "") const;
00152 const double *GetX(int i=0) const;
00153       double *GetX(int i=0);
00154 TCircleFitterAux* GetAux(int i) const;
00155 
00156 static void Test(int iTest=0);
00157 static void TestCorr(int kode=0);
00158 private:
00159 double f();
00160 double df(int i);
00161 double d2f(int i,int j);
00162 double Rho2();
00163 double dRho2(int i);
00164 double d2Rho2(int i,int j);
00165 double F();
00166 double dF(int i);
00167 double d2F(int i,int j);
00168 
00169 private:
00170 TArrayD fArr;
00171 char   fBeg[1];
00172 int    fN;
00173 int    fNuse;
00174 int    fCase;   //default case 0=automatic,1=small angle,2=Chernov/Ososkov
00175 int    fKase;   //used case
00176 int    fBack;
00177 TCircleFitterAux* fAux;
00178 double fCos,fSin;       //Direction of local coordinate
00179 double fNor[2];
00180 double fPol[6];
00181 double fXgravity;
00182 double fYgravity;
00183 double fXx;
00184 double fYy;
00185 double fXy;
00186 double fXrr;
00187 double fYrr;
00188 double fRrrr;
00189 double fRr;
00190 double fWtot;
00191 double fRadius2;
00192 double fXd, fYd, fG1;           //fXd,fYd coordianate of circle center in local sys 
00193 double fNx, fNy;                //Direction from circle center to point on circle in local sys 
00194 double fXCenter,fYCenter;       //coordianate of circle center in globa
00195 double fCov[6],fA,fB,fC,fH;
00196 double fR;                      //radius of curvatur
00197 double fRd;                     //distance to center in local system
00198 double fCorrR,fCorrB;
00199 double fChi2;
00200 int    fNdf;
00201 double fZ0,fTanL;
00202 double fChi2Z;
00203 char   fEnd[1];
00204 ClassDef(TCircleFitter,0)
00205 };
00206 
00207 
00208 class THelixTrack : public TObject 
00209 {
00210 public:
00211 
00212         THelixTrack();
00213         THelixTrack(const double *xyz,const double *dir,double rho,double drho=0);
00214         THelixTrack(const THelixTrack &from);
00215         THelixTrack(const THelixTrack *from);   //Special ctr without errs
00216 virtual ~THelixTrack();
00217 THelixTrack &operator=(const THelixTrack &from);
00218         void Set   (const double *xyz,const double *dir,double rho,double drho=0);
00219         void Set   (double rho,double drho=0);
00220         void SetEmx(const double*  err2xy,const double*  err2z);
00221         void SetEmx(const double*  err=0);
00222     THEmx_t *Emx() const                        {return fEmx;}
00223         void StiEmx(double emx[21]) const;
00224         void GetSpot(const double axis[3][3],double emx[3]) const;
00225         void Fill  (TCircle &circ) const;
00227         void Backward();
00229         double Move(double step);
00230         double Move(double step,double F[5][5]);
00232         double Eval(double step, double *xyz, double *dir,double &rho) const;
00233         double Step(double step, double *xyz, double *dir,double &rho) const
00234        {return Eval(       step,         xyz,         dir,        rho);}
00236         void   Get (double *xyz, double *dir,double &rho) const {Step(0.,xyz,dir,rho);}
00237         double Eval(double step, double *xyz, double *dir=0) const;
00238         double Step(double step, double *xyz, double *dir=0) const
00239        {return Eval(       step,         xyz,         dir  );}
00240         void   Get (double *xyz, double *dir=0) const {Step(0.,xyz,dir);}
00246         double Step(double stmax, const double *surf, int nsurf
00247                    ,double *x=0, double *dir=0, int nearest=0) const;
00248         double Path(double stmax, const double *surf, int nsurf
00249                    ,double *x=0, double *dir=0, int nearest=0) const
00250                    {return Step(stmax,surf,nsurf,x,dir,nearest);}
00251 
00253         double Step(const double point[3],double *xyz=0, double *dir=0) const;
00254         double Path(const double point[3],double *xyz=0, double *dir=0) const 
00255                    {return Step(point,xyz,dir);}
00257         double Dca(const double point[3],double *dcaErr=0) const;
00258 
00260         double Path(double x,double y) const ;
00262         double Dca(double x,double y,double *dcaErr=0) const ;
00263 
00267 
00268         double Path(const THelixTrack &hlx,double *s2=0) const ;
00269 
00274         double PathX(const THelixTrack &hlx,double *s2=0
00275                     ,double *dist=0, double *xyz=0) const;
00276 
00278         double Dca(const double point[3]
00279                   ,double &dcaXY,double &dcaZ,double dcaEmx[3],int kind=3) const;
00280 
00281         const double *GetXYZ()  const {return fX;}
00282         const double *Pos()     const {return fX;}
00283               double *Pos()           {return fX;}
00284         const double *GetDir()  const {return fP;}
00285         const double *Dir()     const {return fP;}
00286               double *Dir()           {return fP;}
00287         double GetRho()         const {return fRho ;}
00288         double GetDRho()        const {return fDRho ;}
00289         double GetCos()         const {return fCosL;}
00290         double GetSin()         const {return fP[2];}
00291         double GetTan()         const {return fP[2]/fCosL;}
00292         double GetPeriod() const ;
00293         void Rot(double angle);
00294         void Rot(double cosa,double sina);
00295 
00296         void Show(double len, const THelixTrack *other=0) const;
00297         void Print(Option_t *opt="") const;
00298 //      statics
00299 static  void InvertMtx(double derivs[5][5]);
00300 static  void Test1();
00301 static  void Test2();
00302 static  void Test3();
00303 static  void Test4();
00304 static  void Test5();
00305 static  void TestMtx();
00306 static  void TestDer();
00307 private:
00310         void MakeMtx(double step,double F[5][5]);
00311 protected:
00312         double Step(double stmin,double stmax, const double *surf, int nsurf
00313                    ,double *x=0, double *dir=0,int nearest=0) const;
00314         double StepHZ(const double *surf, int nsurf
00315                    ,double *x=0, double *dir=0,int nearest=0) const;
00316         void Build();
00317         char   fBeg[1];  
00318         double fX[3];
00319         double fP[3];
00320         double fRho;
00321         double fDRho;
00322         double fCosL;
00323         THEmx_t *fEmx;
00324         char fEnd[1];
00325 ClassDef(THelixTrack,0)
00326 };
00327 class THelixFitter: public THelixTrack
00328 {
00329 public:
00330        THelixFitter();
00331       ~THelixFitter();
00332 int    Size() const                     {return fCircleFitter.Size();}
00333 int    Used() const                     {return fCircleFitter.Used();}
00334 void   Add (double x,double y,double z); 
00335 void   AddErr(const double *err2xy,double err2z); 
00336 void   AddErr(double errhh,double errzz); 
00337 double Fit();   
00338 void   MakeErrs();
00339 double FixAt(const double vals[5],int flag=1); 
00340 void   Skip(int idx); 
00341 void   SetCase(int kase=0)              {fCircleFitter.SetCase(kase);}
00342 int    GetCase() const                  {return fCircleFitter.GetCase();}
00343 double Chi2() const                     {return fChi2;}
00344 int    Ndf()  const                     {return fCircleFitter.Ndf()+fPoli1Fitter.Ndf();}
00345 double Chi2XY () const                  {return fCircleFitter.Chi2();}
00346 double Chi2SZ () const                  {return fPoli1Fitter.Chi2() ;}
00347 int    NdfXY ()  const                  {return fCircleFitter.Ndf();}
00348 int    NdfSZ ()  const                  {return fPoli1Fitter.Ndf() ;}
00349 TCircleFitterAux* GetAux(int i) const   {return fCircleFitter.GetAux(i);}
00350 double EvalChi2();
00351 void   Clear(const char *opt ="");
00352 void   Print(const char* chopt = "") const;
00353 void   Show() const;
00354 
00355 static void Test(int kase=0);
00356 private:
00357 void Update(int kase);
00358 private:
00359 TCircleFitter fCircleFitter;
00360 TPoliFitter   fPoli1Fitter;
00361 double fChi2;
00362 ClassDef(THelixFitter,0)
00363 };
00364 #endif // THELIXTRACK_H