StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
THelixTrack.h
1 #ifndef THELIXTRACK_H
2 #define THELIXTRACK_H
3 #include "TObject.h"
4 #include "TArrayD.h"
5 #include "TPolinom.h"
6 //..............................................................................
7 class TCEmx_t
8 {
9 public:
10 const double *Arr() const { return &mHH;}
11  double *Arr() { return &mHH;}
12 const double &operator[](int idx) const { return (&mHH)[idx];}
13  double &operator[](int idx) { return (&mHH)[idx];}
14  TCEmx_t &operator*=(double f)
15  { for (int i=0;i<6;i++) {Arr()[i]*=f;} return *this;}
16 void Clear() { memset(this,0,sizeof(*this));}
17  TCEmx_t() { Clear();}
18 void Set(const double *err);
19 void Move(double const F[3][3]);
20 double Sign() const;
21 void Backward();
22 public:
23 
24 double
25 mHH,
26 mHA, mAA,
27 mHC, mAC, mCC;
28 };
29 
30 //..............................................................................
31 class THEmx_t
32 {
33 public:
34  THEmx_t() { Clear();}
35 operator const double* () { return &mHH;}
36 operator double* () { return &mHH;}
37 const double *Arr() const { return &mHH;}
38  double *Arr() { return &mHH;}
39 const double &operator[](int idx) const { return (&mHH)[idx];}
40  double &operator[](int idx) { return (&mHH)[idx];}
41  THEmx_t &operator*=(double f)
42  { for (int i=0;i<15;i++) {Arr()[i]*=f;} return *this;}
43 void Clear() { memset(this,0,sizeof(*this));}
44 void Set(const double *err);
45 void Set(const double *errxy,const double *errz);
46 void Move(double const F[5][5]);
47 void Backward();
48 void Print(const char *tit=0) const;
49 double Sign() const;
50 double MaxCorr() const;
51 public:
52 // dA: delta azimuth angle; dH: error along ort to dir and Z axis
53 // dC: error of curvature; dZ == dZ; dL = dLambda
54 double
55 mHH,
56 mHA, mAA,
57 mHC, mAC, mCC,
58 mHZ, mAZ, mCZ, mZZ,
59 mHL, mAL, mCL, mZL, mLL;
60 };
61 
62 
63 //..............................................................................
64 class TCircle: public TObject
65 {
66 friend class THelixTrack;
67 friend class TCircleFitter;
68 public:
69  TCircle();
70  TCircle(const double *x,const double *dir,double rho);
71  TCircle(const TCircle& fr);
72  TCircle(const TCircle* fr); //special ctr without errors
73 ~TCircle();
74 TCircle &operator=(const TCircle& fr);
75  void Set(const double *x=0,const double *dir=0,const double rho=0);
76 virtual void Clear(const char *opt="");
77 const double* Pos() const {return fX; }
78  double* Pos() {return fX; }
79 const double* Dir() const {return fD; }
80  double Rho() const {return fRho;}
81  double& Rho() {return fRho;}
82  void Nor(double *norVec) const;
83  void SetEmx(const double *err=0);
84 const TCEmx_t *Emx() const {return fEmx;}
85  TCEmx_t *Emx() {return fEmx;}
86  void GetCenter(double center[2]) const;
87 double Path(const double pnt[2]) const;
88 double Path(const double pnt[2], const double exy[3]) const;
89 double Path(const TCircle &tc,double *s2=0) const;
90 double Move(double step);
91 void Rot(double angle);
92 void Rot(double cosa,double sina);
93 void Backward();
94 double Eval(double step,double *xy,double *dir=0) const;
95 void Show(int nPts,const double *Pts,int pstep=2) const;
96 virtual void Print(const char* chopt = "") const;
97 void SetStrait(int strait=1) {SetBit(1,strait) ;}
98 int IsStrait() {return TestBit(1);}
99 
100 // static funs
101 static void Test2();
102 static void Test3();
103 static void Test4();
104 static void TestMtx();
105 
106 private:
107 void MoveErrs(double l);
108 void MakeMtx (double l,double F[3][3]);
109 
110 protected:
111 double fX[2];
112 double fD[2];
113 double fRho;
114 TCEmx_t *fEmx; //let h = fX[1]*fD[0], a=atan2(fD[1],fD[0]),c=fRho
115  // hh,
116  // ah,aa,
117  // ch,ca,cc
118 ClassDef(TCircle,0)
119 };
120 //..............................................................................
122 {
123  public:
124  static int dSize() {return sizeof(TCircleFitterAux)/sizeof(double);}
125  public:
126  double x,y,z; //x,y,z of measured point
127  double exy[3]; //err matrix(xx,xy,yy) of x,y
128  double ezz; //error of z
129  double wt; //calculated weight
130 
131 };
132 //..............................................................................
133 class TCircleFitter: public TCircle
134 {
135 public:
136  TCircleFitter();
137 int Size() const {return fN;}
138 int Used() const {return fNuse;}
139 void Add (double x,double y,const double *errs=0);
140 void Add (double x,double y,double z);
141 void AddErr(const double *errs,double errz=0);
142 void AddErr(double errh,double errz=0);
143 void AddZ(double z,double err2z=0);
144 double Fit();
145 void MakeErrs();
146 double FixAt(const double vals[5],int flag);
147 void Skip(int idx);
148 double GetZ0() const {return fZ0 ;}
149 double GetTanL() const {return fTanL ;}
150 void SetCase(int kase=0) {fCase=kase ;}
151 int GetCase() const {return fKase ;}
152 double Chi2() const {return fChi2 ;}
153 int Ndf() const {return fNdf ;}
154 double Chi2Z () const {return fChi2Z ;}
155 void SetNdf(int ndf);
156 double EvalChi2();
157 void Clear(const char *opt ="");
158 void Print(const char* chopt = "") const;
159 const double *GetX(int i=0) const;
160  double *GetX(int i=0);
161 void Show() const;
162 TCircleFitterAux* GetAux(int i) const;
163 
164 static void Test(int iTest=0);
165 static void TestCorr(int kode=0);
166 private:
167 double f();
168 double df(int i);
169 double d2f(int i,int j);
170 double Rho2();
171 double dRho2(int i);
172 double d2Rho2(int i,int j);
173 double F();
174 double dF(int i);
175 double d2F(int i,int j);
176 
177 private:
178 TArrayD fArr;
179 char fBeg[1];
180 int fN;
181 int fNuse;
182 int fCase; //default case 0=automatic,1=small angle,2=Chernov/Ososkov
183 int fKase; //used case
184 int fBack;
185 TCircleFitterAux* fAux;
186 double fCos,fSin; //Direction of local coordinate
187 double fNor[2];
188 double fPol[6];
189 double fXgravity;
190 double fYgravity;
191 double fXx;
192 double fXy;
193 double fYy;
194 double fXrr;
195 double fYrr;
196 double fRrrr;
197 double fRr;
198 double fWtot;
199 double fRadius2;
200 double fXd, fYd, fG1; //fXd,fYd coordianate of circle center in local sys
201 double fNx, fNy; //Direction from circle center to point on circle in local sys
202 double fXCenter,fYCenter; //coordianate of circle center in globa
203 double fCov[6],fA,fB,fC,fH;
204 double fR; //radius of curvatur
205 double fRd; //distance to center in local system
206 double fCorrR,fCorrB;
207 double fChi2;
208 int fNdf;
209 double fZ0,fTanL;
210 double fChi2Z;
211 char fEnd[1];
212 ClassDef(TCircleFitter,0)
213 };
214 
215 
216 //..............................................................................
217 class THelixTrack : public TObject
218 {
219 public:
220 
221  THelixTrack();
222  THelixTrack(const double *xyz,const double *dir,double rho,double drho=0);
223  THelixTrack(const THelixTrack &from);
224  THelixTrack(const THelixTrack *from); //Special ctr without errs
225 virtual ~THelixTrack();
226 THelixTrack &operator=(const THelixTrack &from);
227  void Set (const double *xyz,const double *dir,double rho,double drho=0);
228  void Set (double rho,double drho=0);
229  void SetEmx(const double* err2xy,const double* err2z);
230  void SetEmx(const double* err=0);
231  THEmx_t *Emx() const {return fEmx;}
232  void StiEmx(double emx[21]) const;
233  void GetSpot(const double axis[3][3],double emx[3]) const;
234  void Fill (TCircle &circ) const;
236  void Backward();
238  double Move(double step);
239  double Move(double step,double F[5][5]);
241  double Eval(double step, double *xyz, double *dir,double &rho) const;
242  double Step(double step, double *xyz, double *dir,double &rho) const
243  {return Eval( step, xyz, dir, rho);}
245  void Get (double *xyz, double *dir,double &rho) const {Step(0.,xyz,dir,rho);}
246  double Eval(double step, double *xyz, double *dir=0) const;
247  double Step(double step, double *xyz, double *dir=0) const
248  {return Eval( step, xyz, dir );}
249  void Get (double *xyz, double *dir=0) const {Step(0.,xyz,dir);}
255  double Step(double stmax, const double *surf, int nsurf
256  ,double *x=0, double *dir=0, int nearest=0) const;
257  double Path(double stmax, const double *surf, int nsurf
258  ,double *x=0, double *dir=0, int nearest=0) const
259  {return Step(stmax,surf,nsurf,x,dir,nearest);}
260 
262  double Step(const double point[3],double *xyz=0, double *dir=0) const;
263  double Path(const double point[3],double *xyz=0, double *dir=0) const
264  {return Step(point,xyz,dir);}
266  double Dca(const double point[3],double *dcaErr=0) const;
267 
269  double Path(double x,double y) const ;
271  double Dca(double x,double y,double *dcaErr=0) const ;
272 
276 
277  double Path(const THelixTrack &hlx,double *s2=0) const ;
278 
283  double PathX(const THelixTrack &hlx,double *s2=0
284  ,double *dist=0, double *xyz=0) const;
285 
287  double Dca(const double point[3]
288  ,double &dcaXY,double &dcaZ,double dcaEmx[3],int kind=3) const;
289 
290  const double *GetXYZ() const {return fX;}
291  const double *Pos() const {return fX;}
292  double *Pos() {return fX;}
293  const double *GetDir() const {return fP;}
294  const double *Dir() const {return fP;}
295  double *Dir() {return fP;}
296  double GetRho() const {return fRho ;}
297  double GetDRho() const {return fDRho ;}
298  double GetCos() const {return fCosL;}
299  double GetSin() const {return fP[2];}
300  double GetTan() const {return fP[2]/fCosL;}
301  double GetPeriod() const ;
302  void Rot(double angle);
303  void Rot(double cosa,double sina);
304 
305  void Show(double len, const THelixTrack *other=0) const;
306  void Print(Option_t *opt="") const;
307 // statics
308 static void InvertMtx(double derivs[5][5]);
309 static void Test1();
310 static void Test2();
311 static void Test3();
312 static void Test4();
313 static void Test5();
314 static void TestMtx();
315 static void TestDer();
316 static void TestTwoHlx();
317 private:
320  void MakeMtx(double step,double F[5][5]);
321 protected:
322  double Step(double stmin,double stmax, const double *surf, int nsurf
323  ,double *x=0, double *dir=0,int nearest=0) const;
324  double StepHZ(const double *surf, int nsurf
325  ,double *x=0, double *dir=0,int nearest=0) const;
326  void Build();
327  char fBeg[1];
328  double fX[3];
329  double fP[3];
330  double fRho;
331  double fDRho;
332  double fCosL;
333  THEmx_t *fEmx;
334  char fEnd[1];
335 ClassDef(THelixTrack,0)
336 };
338 {
339 public:
340  THelixFitter();
341  ~THelixFitter();
342 int Size() const {return fCircleFitter.Size();}
343 int Used() const {return fCircleFitter.Used();}
344 void Add (double x,double y,double z);
345 void AddErr(const double *err2xy,double err2z);
346 void AddErr(double errhh,double errzz);
347 double Fit();
348 void MakeErrs();
349 double FixAt(const double vals[5],int flag=1);
350 void Skip(int idx);
351 void SetCase(int kase=0) {fCircleFitter.SetCase(kase);}
352 int GetCase() const {return fCircleFitter.GetCase();}
353 double Chi2() const {return fChi2;}
354 int Ndf() const {return fCircleFitter.Ndf()+fPoli1Fitter.Ndf();}
355 double Chi2XY () const {return fCircleFitter.Chi2();}
356 double Chi2SZ () const {return fPoli1Fitter.Chi2() ;}
357 int NdfXY () const {return fCircleFitter.Ndf() ;}
358 int NdfSZ () const {return fPoli1Fitter.Ndf() ;}
359 TCircleFitterAux* GetAux(int i) const {return fCircleFitter.GetAux(i);}
360 double EvalChi2();
361 void Clear(const char *opt ="");
362 void Print(const char* chopt = "") const;
363 void Show() const;
364 
365 static void Test(int kase=0);
366 private:
367 void Update(int kase);
368 private:
369 TCircleFitter fCircleFitter;
370 TPoliFitter fPoli1Fitter;
371 double fChi2;
372 ClassDef(THelixFitter,0)
373 };
374 
375 //..............................................................................
376 #include <vector>
377 class THelixKFitterAux { public: double x[3],e[6],wt,xi2; };
378 typedef std::vector<THelixKFitterAux> THelixKFitterAuxV;
379 
381 {
382 public:
383  THelixKFitter() {fFitingShow=0;Clear();}
384  ~THelixKFitter(){;}
385 void Add (const double x[3]);
386 void AddErr(const double err[6]);
387 double Fit();
388 double Chi2() const {return fChi2 ;}
389 int Ndf() const {return 2*fAux.size()-5 ;}
390 int Size() const {return fAux.size() ;}
391 const THelixKFitterAux* GetAux(int i) const {return &fAux[i] ;}
392 void Clear(const char * ="") {fAux.clear();fChi2=0 ;}
393 void Print(const char* chopt = "") const;
394 void Show() const;
395 void SetFitingShow() {fFitingShow = new std::vector<double>;}
396 static void Test(int nev=10000);
397 private:
398 std::vector<double> *fFitingShow;
399 THelixKFitterAuxV fAux;
400 double fChi2;
401 ClassDef(THelixKFitter,0)
402 };
403 #endif // THELIXTRACK_H
double Eval(double step, double *xyz, double *dir, double &rho) const
Evaluate params with given step along helix.
double PathX(const THelixTrack &hlx, double *s2=0, double *dist=0, double *xyz=0) const
void Get(double *xyz, double *dir, double &rho) const
Get current parameters.
Definition: THelixTrack.h:245
void Backward()
Change direction.
double Move(double step)
Move along helix.
Definition: AgUStep.h:71
void GetSpot(const double axis[3][3], double emx[3]) const
double FixAt(const double vals[5], int flag)
double Dca(const double point[3], double *dcaErr=0) const
DCA to given space point (with error matrix)