StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
THelix3d.cxx
1 #include <stdlib.h>
2 #include <math.h>
3 #include "TError.h"
4 #include "TArrayD.h"
5 #include "TCanvas.h"
6 #include "TGraph.h"
7 #include "TSystem.h"
8 #include "TMath.h"
9 #include "TCernLib.h"
10 #include "TVector3.h"
11 #include "TVectorD.h"
12 #include "TRandom.h"
13 #include "TRandom2.h"
14 #include "THelix3d.h"
15 #include "TRandomVector.h"
16 #include "StMatrixD.hh"
17 #include "TH1.h"
18 #include <cassert>
19 
20 
21 static double EmxSign(int n,const double *a);
22 
23 #if 0
24 //_____________________________________________________________________________
25 inline static void spheric(const double D[3]
26  ,double &cosL,double &sinL
27  ,double &cosP,double &sinP)
28 {
29  sinL = D[2];
30  cosL = (1.-sinL)*(1+sinL);
31  cosL = (cosL<=0)? 0.:sqrt(cosL);
32  if (cosL <1e-6) {//track along Z-axis
33  cosP = 1; sinP = 0; cosL = 1e-6;
34  } else {
35  cosP = D[0]/cosL; sinP = D[1]/cosL;
36  }
37 }
38 //_____________________________________________________________________________
39 inline static void satlit(const double &cosL,const double &sinL
40  ,const double &cosP,const double &sinP
41  ,double U[3][3])
42 {
43  U[0][0]= cosL*cosP; U[0][1] = cosL*sinP;U[0][2]= sinL;
44  U[1][0]=- sinP; U[1][1] = cosP; U[1][2]=0;
45  U[2][0]=-sinL*cosP; U[2][1] = -sinL*sinP; U[2][2]=cosL;
46 }
47 #endif
48 //_____________________________________________________________________________
49 //#define dot(a,b) (a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
50 
51 inline static double dot(const double *a,const double *b)
52 {return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];}
53 
54 
55 //_____________________________________________________________________________
56 inline static double nor(double *a)
57 {
58  double N = dot(a,a); if (fabs(N-1)<1e-7) return N;
59  if (N<=0) return N;
60  N = sqrt(N);
61  a[0]/=N; a[1]/=N; a[2]/=N;
62  return N;
63 }
64 //_____________________________________________________________________________
65 inline static void lin(const double a[3],double f, const double b[3], double c[3])
66 {
67  c[0] = a[0]+f*b[0];c[1] = a[1]+f*b[1];c[2] = a[2]+f*b[2];
68 }
69 //_____________________________________________________________________________
70 inline static void cop(const double a[3],double c[3])
71 { c[0] = a[0]; c[1] = a[1];c[2] = a[2]; }
72 //_____________________________________________________________________________
73 inline static void cro(const double a[3],const double b[3], double c[3])
74 {
75  c[0] = a[1]*b[2]-b[1]*a[2];
76  c[1] = a[2]*b[0]-b[2]*a[0];
77  c[2] = a[0]*b[1]-b[0]*a[1];
78 }
79 //_____________________________________________________________________________
80 //_____________________________________________________________________________
81 //_____________________________________________________________________________
82 //_____________________________________________________________________________
83 //_____________________________________________________________________________
84 ClassImp(THelix3d)
85 //_____________________________________________________________________________
86 // in THDer3d_t direction of track represented as
87 /*
88  mTkDir[*][2]*cos(alfa)*cos(beta) +
89  mTkDir[*][1]*cos(alfa)*sin(beta) +
90  mTkDir[*][0]*sin(alfa)
91 
92 Or when alfa & beta are small
93  mTkDir[*][2] +
94  mTkDir[*][1]*beta +
95  mTkDir[*][0]*alfa
96 */
97 //_____________________________________________________________________________
98 void TkDir_t::Backward()
99 {
100 // U = 0.0233245 0.828218 -0.55992
101 // Uinv= -0.0233245 -0.828218 0.55992
102 // V = -0.30932 0.538569 0.78375
103 // Vinv= -0.30932 0.538569 0.78375
104 // kKU is immediately after kKT
105  TCL::vscale(mTkd[kKT],-1.,mTkd[kKT],3*2);
106 }
107 //_____________________________________________________________________________
108 int TkDir_t::Same(const TkDir_t &as) const
109 {
110 static const double eps = 0.01;
111  for (int i=0;i<3;i++) {
112  for (int j=0;j<3; j++) {
113  if (fabs(mTkd[i][j]-as[i][j])>eps) return 0;
114  } }
115  return 1;
116 }
117 //_____________________________________________________________________________
118 void THDer3d_t::Clear()
119 {
120  mLen=0;
121  memset(mDer[0] ,0,sizeof(mDer ));
122  for (int i=0;i< 5;i++) { mDer[i][i] = 1;}
123 
124  for (int jk=0; jk<2;jk++) { mTkDir[jk].Clear();}
125 }
126 //_____________________________________________________________________________
127 THDer3d_t &THDer3d_t::operator*=(const THDer3d_t &by)
128 {
129  Mtx55D_t res;
130  TCL::mxmpy(by.mDer[0],mDer[0],res[0],5,5,5);
131  TCL::ucopy(res[0],mDer[0],5*5);
132  TCL::ucopy(by.mTkDir[1][0],mTkDir[1][0],3*3);
133  mLen+=by.mLen;
134  return *this;
135 }
136 //_____________________________________________________________________________
137 void THDer3d_t::Backward()
138 {
139 // kU = U*dX = -1
140 // kV = V*dX = +1
141 // kT -1
142 // kFita = kU*dT = 1
143 // kLama = kV*dT = -1
144 // kPinv = -1
145 // kU kV kFita kLama kPinv
146 // kU 1 -1 -1 1 1
147 // kV -1 1 1 -1 -1
148 // kFita -1 1 1 -1 -1
149 // kLama 1 -1 -1 1 1
150 // kPinv 1 -1 -1 1 1
151 
152  mTkDir[0].Backward();
153  mTkDir[1].Backward();
154 // U, V, Fita, Lama, Pinv... = 0.001 0.002 0.003 0.004 0.001
155 // Um,Vm,FitaM, LamaM.PinvM.. = -0.001 0.002 0.003 -0.004 -0.001
156 // kU-,kV+,kFita+,kLama-,kPinv-
157 
158 static const int idx[]= {kU,kV , kU,kFita
159  ,kV,kU , kV,kLama , kV,kPinv
160  ,kFita,kU, kFita,kLama, kFita,kPinv
161  ,kLama,kV, kLama,kFita
162  ,kPinv,kV, kPinv,kFita, -1,-1};
163 
164  for (int i=0;idx[i]>=0; i+=2) {mDer[idx[i]][idx[i+1]]*=-1.;}
165 }
166 //_____________________________________________________________________________
167 THelix3d::THelix3d(int charge,const double *xyz,const double *mom,const double *Mag)
168 {
169 // Generalization of THelixTrack_ for arbitrary direction of mag field
170 // V.Perevoztchikov
171 //
172  memset(fBeg3d,0,fEnd3d-fBeg3d);
173  Set(charge,xyz,mom,Mag);
174 }
175 //_____________________________________________________________________________
176 THelix3d::THelix3d()
177 {
178  memset(fBeg3d,0,fEnd3d-fBeg3d);
179 }
180 //_____________________________________________________________________________
181 THelix3d::~THelix3d()
182 {
183  delete fEmx3d;
184  delete fDer3d;
185 }
186 //_____________________________________________________________________________
187 void THelix3d::Clear(const char*)
188 {
189  memset(fBeg3d,0,fMed3d-fBeg3d);
190 
191  if (fDer3d) fDer3d->Clear();
192  if (fEmx3d) fEmx3d->Clear();
193 }
194 //_____________________________________________________________________________
195 void THelix3d::Set(int charge,const double *xyz,const double *mom,const double *mag)
196 {
197  if (charge) fCharge = charge;
198  if (xyz) memcpy(fX3d,xyz,sizeof(fX3d));
199  if (mom) memcpy(fP3d,mom,sizeof(fP3d));
200  if (mag) memcpy(fH3d,mag,sizeof(fH3d[0])*3);
201  Build();
202 }
203 //_____________________________________________________________________________
204 void THelix3d::Build()
205 {
206  fMom = sqrt(TCL::vdot(fP3d,fP3d,3));
207  TCL::vscale(fP3d,1./fMom,fD3d,3);
208  fH3d[3] = sqrt(TCL::vdot(fH3d,fH3d,3));
209  MakeLocal(fH3d,fD3d,fLoc);
210  ToLocal();
211  fPinv = -fCharge/fMom;
212  double pt = sqrt(TCL::vdot(fP,fP,2))*fMom;
213  fRho = -fCharge*fH3d[3]/pt;
214  THelixTrack_::Build();
215  if (!fEmx3d) return;
216  TkDir_t tkDir;
217  DoTkDir(tkDir);
218  fEmx3d->Update(tkDir);
219  fDer3d->Clear();
220  fDer3d->SetTkDir(0,tkDir);
221  fDer3d->SetTkDir(1,tkDir);
222 }
223 //_____________________________________________________________________________
224 THelix3d::THelix3d(const THelix3d *from): THelixTrack_(from)
225 {
226  memcpy(fBeg3d,from->fBeg3d,fEnd3d-fBeg3d);
227  fEmx3d=0; fEmx=0;
228 }
229 //_____________________________________________________________________________
230 THelix3d::THelix3d(const THelix3d &from): THelixTrack_(&from)
231 {
232  memcpy(fBeg3d,from.fBeg3d,fEnd3d-fBeg3d);
233  if (fEmx3d) {fEmx3d = new THEmx3d_t; *fEmx3d = *from.fEmx3d;}
234  if (fDer3d) {fDer3d = new THDer3d_t; *fDer3d = *from.fDer3d;}
235 }
236 //_____________________________________________________________________________
237 void THelix3d::ToGlobal(const double locX[3],const double locD[3]
238  , double gloX[3], double gloD[3]
239  , double gloP[3])const
240 // Convert local variables(THelixTrack_) into global (THelix3d)
241 
242 {
243  double myX[3];
244  TCL::vmatr (locX,fLoc[0],myX ,3,3);
245  TCL::vadd (myX ,fX3d ,gloX,3 );
246  TCL::vmatr (locD,fLoc[0],gloD,3,3);
247  TCL::vscale(gloD,fMom ,gloP,3 );
248 }
249 //_____________________________________________________________________________
250 void THelix3d::GetdDdL(double dDdL[3]) const
251 {
252  double dDdLloc[3] = {-fP[1]*fRho*fCosL,fP[0]*fRho*fCosL,0};
253  TCL::vmatr (dDdLloc,fLoc[0],dDdL,3,3);
254 
255 }//_____________________________________________________________________________
256 void THelix3d::ToGlobal()
257 {
258  ToGlobal(fX,fP,fX3d,fD3d,fP3d);
259 }
260 //_____________________________________________________________________________
261 void THelix3d::ToLocal()
262 {
263  TCL::vzero(fX,3);
264  TCL::vmatl(fLoc[0],fD3d,fP,3,3);
265  memcpy(fPpre,fP,sizeof(fPpre));
266  fCosL = sqrt(fabs((1.-fP[2])*(1+fP[2])));
267 }
268 //_____________________________________________________________________________
269 void THelix3d::Backward()
270 {
271  TCL::vscale(fP3d ,-1.,fP3d ,3);
272  TCL::vscale(fD3d ,-1.,fD3d ,3);
273  fCharge = - fCharge;
274  fPinv = - fPinv;
275  THelixTrack_::Backward();
276  if (fEmx3d) fEmx3d->Backward();
277  if (fDer3d) fDer3d->Backward();
278  ToLocal();
279 }
280 //_____________________________________________________________________________
281 double THelix3d::Move(double step,double F[5][5])
282 {
283  if (fabs(step) <1e-11) return step;
284  memcpy(fX3dPre,fX3d,sizeof(fX3dPre));
285  memcpy(fP3dPre,fP3d,sizeof(fP3dPre));
286  memcpy(fD3dPre,fD3d,sizeof(fD3dPre));
287  if (!IsDerOn() && !F) { //No derivatives needed
288  ToLocal();
289  fLen = THelixTrack_::Move(step);
290  ToGlobal();
291  } else {
292  ToLocal();
293  fLen = THelixTrack_::Move(step);
294  ToGlobal();
295  MakeMtx();
296  if (F) memcpy(F[0],fDer3d[0],sizeof(double)*5*5);
297  if (fEmx3d) fEmx3d->Move(fDer3d);
298 
299  }
300 
301  return fLen;
302 }
303 //_____________________________________________________________________________
304 double THelix3d::Eval(double step, double xyz[3], double mom[3])
305 {
306  fLen = step;
307  if (fabs(step)<1e-11) {
308  if (xyz) memcpy(xyz,fX3d,sizeof(fX3d));
309  if (mom) memcpy(mom,fP3d,sizeof(fP3d));
310  } else {
311  ToLocal();
312  double myXE[3],myDE[3],myX3dE[3],myD3dE[3],myP3dE[3];
313  THelixTrack_::Eval(step,myXE,myDE);
314  ToGlobal(myXE,myDE,myX3dE,myD3dE,myP3dE);
315  if (xyz) memcpy(xyz,myX3dE,sizeof(myX3dE));
316  if (mom) memcpy(mom,myP3dE,sizeof(myP3dE));
317  }
318  return step;
319 }
320 //_____________________________________________________________________________
321 double THelix3d::Path(const double point[3],double xyz[3], double mom[3])
322 {
323  double myPoint[6],myXyz[3],myMom[3];
324  ToLocal();
325  TCL::vsub(point,fX3d,myPoint+3,3);
326  TCL::vmatl(fLoc[0],myPoint+3,myPoint,3,3);
327  double s = THelixTrack_::Path(myPoint,myXyz,myMom);
328  if (xyz) {
329  TCL::vmatr(myXyz,fLoc[0],xyz,3,3);
330  TCL::vadd(xyz,fX3d,xyz,3);
331  }
332  if (mom) {
333  TCL::vmatr(mom,fLoc[0],myMom,3,3);
334  TCL::vscale(mom,fMom,mom,3);
335  }
336  return s;
337 }
338 //______________________________________________________________________________
339 double THelix3d::Dca(const double point[3],double *dcaErr)
340 {
341 if(point || dcaErr){};
342 assert(!"THelix3d::Dca Not implemented");
343 return 0;
344 }
345 //______________________________________________________________________________
346 double THelix3d::Path(double x,double y)
347 {
348 static const double kEps=1e-5,kRef=1e-4,kBigEps=1e-1;
349  double keepStep = 1e11;
350  THelix3d TH(this);
351  const double *X = TH.Pos();
352  const double *D = TH.Dir();
353  double maxStep = GetPeriod()/6.28;
354  maxStep/= sqrt(1. - D[2]*D[2])+1e-11;
355  fLen = 0;
356  double qaLeft,posLeft=0,qaRite=0,posRite=0,qaMine,qaStep,dx[2];
357  for (int iter=0; iter<=20; iter++) {
358  dx[0] = x-X[0];dx[1] = y-X[1];
359  qaMine = dx[0]*D[0] + dx[1]*D[1];
360  qaMine /=(1.-D[2])*(1.+D[2]);
361  if (fabs(qaMine) < kEps || fabs(qaMine) < kRef*fabs(fLen)) goto END;;
362  if (!iter) {qaLeft = qaMine; qaStep = qaMine; }
363  else if (qaMine*qaLeft>0) {qaLeft = qaMine; qaStep*= 2 ; }
364  else {qaRite = qaMine; posLeft=-qaStep; break ;}
365  if (fabs(qaStep)>maxStep) qaStep = (qaStep<0)? -maxStep:maxStep;
366  TH.Move(qaStep) ;
367  fLen += qaStep;
368 // assert(fabs(fLen)<4*maxStep);
369  if (fabs(fLen)>4*maxStep) return 1e11;
370  }
371  assert(qaRite);
372  for (int iter=0; iter<=20; iter++) {
373  double denom = qaRite-qaLeft;
374  if (iter&1 || fabs(denom) < 1e-4) { //Half step
375  qaStep = 0.5*(posRite+posLeft);
376  } else { //Linear approach
377  qaStep = -(qaLeft*posRite-qaRite*posLeft)/denom;
378  }
379  TH.Move(qaStep) ;
380  fLen += qaStep;
381  dx[0] = x-X[0];dx[1] = y-X[1];
382  qaMine = dx[0]*D[0]+dx[1]*D[1];
383  qaMine /=(1.-D[2])*(1.+D[2]);
384  if (fabs(qaMine) < kEps || fabs(qaMine) < kRef*fabs(fLen)) goto END;;
385  if (qaLeft*qaMine>0) { qaLeft = qaMine; posLeft=0; posRite-=qaStep;}
386  else { qaRite = qaMine; posRite=0; posLeft-=qaStep;}
387  assert(fabs(posRite-posLeft)<=keepStep);
388  keepStep= fabs(posRite-posLeft);
389  }
390 
391 END: return (fabs(qaMine)<kBigEps)? fLen:1e11;
392  }
393 //______________________________________________________________________________
394 double THelix3d::Dca(double x,double y, double *dcaErr)
395 {
396 if (x || y){}
397 return 0;
398 assert(!"THelix3d::Dca(x,y,dcaError) n oy implemented");
399 }
400 
401 //static double mydPdP0[3][3];
402 //static double mydXdP0[3][3];
403 
404 //______________________________________________________________________________
405 void THelix3d::MakeMtx()
406 {
407 static const double Zero[3]= {0,0,1}; //Start of coordinate in local sys
408  if (!fDer3d) fDer3d = new THDer3d_t;
409  auto &der = *fDer3d;
410  der.Clear();
411  TkDir_t &TkDir = der.mTkDir[0];
412  TkDir_t &TkDirE = der.mTkDir[1];
413 
414  assert(fLen!=0);
415  double Dhlx[5][5],to3d[5][5],to3di[5][5],tmp[5*5];
416  THelixTrack_::MakeMtx(fLen,Dhlx);
417  THelixTrack_ preHlx(Zero,fPpre,fRho);
418  ConvertErrs( &preHlx, fH3d[3],0, 0,to3di);
419  ConvertErrs((const THelixTrack_ *)this, fH3d[3],0,to3d, 0);
420 
421  TCL::mxmpy(to3d[0],Dhlx[0],tmp,5,5,5);
422  TCL::mxmpy(tmp,to3di[0],der[0],5,5,5);
423 
424  MakeTkDir(fD3dPre,fH3d,TkDir );
425  MakeTkDir(fD3d ,fH3d,TkDirE);
426 
427 
428  der.mLen = fLen;
429 }
430 
431 //______________________________________________________________________________
432 THEmx3d_t *THelix3d::SetEmx(THEmx3d_t *emx)
433 {
434  SetDerOn();
435  if (!fEmx3d) { fEmx3d = new THEmx3d_t;}
436  else { fEmx3d->Clear() ;}
437  if (!fDer3d) { fDer3d = new THDer3d_t;}
438  else { fDer3d->Clear() ;}
439  if (emx) { *fEmx3d = *emx;};
440  TkDir_t &tkDir = fEmx3d->TkDir();
441  DoTkDir(tkDir);
442  fDer3d->Clear();
443  fDer3d->SetTkDir(0,tkDir);
444  fDer3d->SetTkDir(1,tkDir);
445  return fEmx3d;
446 }
447 //______________________________________________________________________________
448 THEmx3d_t *THelix3d::SetEmx(const double G[15])
449 {
450  SetDerOn();
451  if (!fEmx3d) { fEmx3d = new THEmx3d_t;}
452  else { fEmx3d->Clear() ;}
453  fEmx3d->Set(G);
454  if (!fDer3d) { fDer3d = new THDer3d_t;}
455  else { fDer3d->Clear() ;}
456  TkDir_t &tkDir = fEmx3d->TkDir();
457  DoTkDir(tkDir);
458  fDer3d->Clear();
459  fDer3d->SetTkDir(0,tkDir);
460  fDer3d->SetTkDir(1,tkDir);
461  return fEmx3d;
462 }
463 //______________________________________________________________________________
464 void THelix3d::MakeTkDir( const double T[3],const double H[3],TkDir_t &tkDir)
465 {
466 // H mag field direction
467 // T track direction
468 // Zaxis along T,Yaxis is normal to H and T, X axis (X*H) >=0
469 
470  cop(T,tkDir[kKT]); //Zaxis along T
471  nor(tkDir[kKT]);
472  cop(H,tkDir[kKV]);
473  double N = nor(tkDir[kKV]);
474  if (N<1e-11) {
475  tkDir[kKV][0]=0; tkDir[kKV][1]=0; tkDir[kKV][2]=1; N=1;
476  }
477  double HT = dot(tkDir[kKV],tkDir[kKT]);
478  lin(tkDir[kKV],-HT,tkDir[kKT],tkDir[kKV]);
479  N = nor(tkDir[kKV]); //Y axis = H - (H*T)*T,then norm
480  assert(N>1e-11);
481  cro(tkDir[kKV],tkDir[kKT],tkDir[kKU]); //Yaxis normal to X,Z
482  nor(tkDir[kKU]);
483 
484 #if 1
485  for (int i=0;i< 3;i++) {
486  for (int k=0;k<=i;k++) {
487  double dot = TCL::vdot(tkDir[i],tkDir[k],3);
488  if (i==k) dot--;
489  assert(fabs(dot)<1e-4);
490  }}
491 
492  for (int i=0;i<3;i++) {
493  int j=(i+1)%3;int k=(j+1)%3;
494  double qwe =(TVector3(tkDir[i]).Cross(TVector3(tkDir[j]))).Dot(TVector3(tkDir[k]));
495  assert(fabs(qwe-1)<1e-4);
496  }
497 
498 #endif
499 }
500 //______________________________________________________________________________
501 void THelix3d::DoTkDir(TkDir_t &tkDir)
502 {
503  const double *t =fD3d;
504  MakeTkDir(t,fH3d,tkDir);
505  GetdDdL(tkDir[kKdDdL]);
506 }
507 //______________________________________________________________________________
508 void THelix3d::MakeLocal( const double H[3],const double T[3],double tkDir[3][3])
509 {
510 // H mag field direction
511 // T track direction
512 // Zaxis along T,Yaxis is normal to H and T, X axis (X*H) >=0
513 
514  cop(H,tkDir[2]); //Zaxis along Z
515  nor( tkDir[2]);
516  cop(T,tkDir[1]);
517  nor( tkDir[1]);
518  double HT = dot(tkDir[1],tkDir[2]);
519  lin(tkDir[1],-HT,tkDir[2],tkDir[1]);
520  nor(tkDir[1]); //Y axis = H - (H*T)*T,then norm
521  cro(tkDir[1],tkDir[2],tkDir[0]); //Xaxis normal to Y,Z
522  nor(tkDir[0]);
523 }
524 //_____________________________________________________________________________
525 void THEmx3d_t::Clear()
526 {
527  mLen=0;mTimes[0]=0;mTimes[1]=0;
528  memset(&mUU,0,sizeof(mUU)*15);
529  for (int i=0,li=0;i< 5;li+=++i) { (&mUU)[li+i] = 1;}
530  mTkDir.Clear();
531 }
532 //______________________________________________________________________________
533 void THEmx3d_t::Set(double const err[15],TkDir_t *tkDir)
534 {
535  TCL::ucopy(err,&mUU,15);
536 assert(tkDir || fabs(mTkDir[kKU][0])+fabs(mTkDir[kKU][1])+fabs(mTkDir[kKU][2])>0.1);
537  if (!tkDir) return;
538  TCL::ucopy((*tkDir)[0],mTkDir[0],3*3);
539 }
540 //______________________________________________________________________________
541 void THEmx3d_t::Set(double eUU,double eUV,double eVV)
542 {
543  mUU=eUU;mUV=eUV;mVV=eVV;
544 assert(fabs(mTkDir[kKU][0])+fabs(mTkDir[kKU][1])+fabs(mTkDir[kKU][2])>0.1);
545 }
546 //______________________________________________________________________________
547 void THEmx3d_t::Add(double Theta2,double Orth2,double PinvRr)
548 {
549  mTimes[1]++;
550  mFF+= Theta2;
551  mLL+= Theta2;
552  mUU+= Orth2;
553  mVV+= Orth2;
554  mPP+= PinvRr;
555 }
556 //______________________________________________________________________________
557 void THEmx3d_t::Move(const THDer3d_t *der)
558 {
559  assert(mUU>0);
560  assert(mTkDir.Same(der->TkDir(0)));
561  mTimes[0]++;
562  mLen +=der->Len();
563  double oErr[15];
564  memcpy(oErr,*this,sizeof(oErr));
565  TCL::trasat((*der)[0],oErr,*this,5,5);
566  TCL::ucopy(der->TkDir(1)[0],mTkDir[0],3*3);
567 }
568 //_____________________________________________________________________________
569 void THEmx3d_t::Backward()
570 {
571 // kU kV kFita kLama kPinv
572 // kU 1 -1 -1 1 1
573 // kV -1 1 1 -1 -1
574 // kFita -1 1 1 -1 -1
575 // kLama 1 -1 -1 1 1
576 // kPinv 1 -1 -1 1 1
577 // kU-,kV+,kFita+,kLama-,kPinv-
578  mTkDir.Backward();
579 
580  mUV*=-1; mUF*=-1;
581  mVL*=-1; mVP*=-1;
582  mFL*=-1; mFP*=-1;
583 }
584 //_____________________________________________________________________________
585 void THEmx3d_t::Update(const TkDir_t &tkDir)
586 {
587 // U1U0*fita0 + U1V0*lama0 = fita1
588 // V1U0*fita0 + V1V0*lama0 = lama1
589 //
590 // U1U0*u0 +U1V0*v0 = u1
591 // V1U0*u0 +V1V0*v0 = V1
592 
593  double dots[3][3];
594  enum {kU,kV,kFita,kLama,kPinv};
595 // for (int i=0;i<3;i++) {
596 // for (int j=0;j<3;j++) {
597 // dots[i][j] = dot(tkDir[i],mTkDir[j]);
598 // } }
599  for (int i=kU;i<=kV;i++) {
600  for (int j=kU;j<=kV;j++) {
601  dots[i][j] = dot(tkDir[i+kKU],mTkDir[j+kKU]);
602  } }
603 
604  double T1U0 = dot(tkDir[kKT],mTkDir[kKU]);
605  double T1V0 = dot(tkDir[kKT],mTkDir[kKV]);
606  double UdDdL = dot(tkDir[kKU],tkDir[kKdDdL]);
607  double VdDdL = dot(tkDir[kKV],tkDir[kKdDdL]);
608 
609 
610  double T[5][5]={{0}};
611  T[kU][kU] = dots[kU][kU];
612  T[kU][kV] = dots[kU][kV];
613  T[kV][kU] = dots[kV][kU];
614  T[kV][kV] = dots[kV][kV];
615  T[kFita][kFita] = dots[kU][kU];
616  T[kFita][kLama] = dots[kU][kV];
617  T[kLama][kFita] = dots[kV][kU];
618  T[kLama][kLama] = dots[kV][kV];
619  T[kPinv][kPinv] = 1;
620  T[kFita][kU] = -(UdDdL)*T1U0;
621  T[kFita][kV] = -(UdDdL)*T1V0;
622  T[kLama][kU] = -(VdDdL)*T1U0;
623  T[kLama][kV] = -(VdDdL)*T1V0;
624 
625  double preRR[5*(5+1)/2];
626  memcpy(preRR,&mUU,sizeof(preRR));
627  TCL::trasat(T[0],preRR,&mUU,5,5);
628  assert(mPP>0);
629  memcpy(mTkDir[0],tkDir[0],sizeof(mTkDir));
630 
631 }
632 //______________________________________________________________________________
633 double THEmx3d_t::Sign() const
634 {
635  const double *E = &mUU;
636  return EmxSign(5,E);
637 }
638 //______________________________________________________________________________
639 double THEmx3d_t::Trace() const {return mUU+mVV+mFF+mLL+mPP;}
640 //______________________________________________________________________________
641 void THEmx3d_t::Print(const char *tit) const
642 {
643 static const char *N="UVXFLP";
644  if (!tit) tit = "";
645  printf("THEmx3d_t::::Print(%s) ==\n",tit);
646  const double *e = &mUU;
647  for (int i=0,li=0;i< 5;li+=++i) {
648  printf("%c ",N[i]);
649  for (int j=0;j<=i;j++) {
650  printf("%g\t",e[li+j]);}
651  printf("\n");
652  }
653 }
654 //____________________________________________________________
655 double mySign(const double *a,int n)
656 {
657  enum {kMaxDim=5,kMaxSize=(kMaxDim*(kMaxDim+1))/2 };
658  double ans=3e33;
659  double *aa = (double *)a;
660  double save = aa[0]; if (!save) aa[0] = 1;
661  double B[kMaxSize];
662  double *myb = B;
663  if (n>kMaxDim) myb = new double[(n*(n+1))/2];
664  double *b = myb;
665  // trchlu.F -- translated by f2c (version 19970219).
666  //
667  //see original documentation of CERNLIB package F112
668 
669  /* Local variables */
670  int ipiv, kpiv, i__, j;
671  double r__, dc;
672  int id, kd;
673  double sum;
674 
675 
676  /* CERN PROGLIB# F112 TRCHLU .VERSION KERNFOR 4.16 870601 */
677  /* ORIG. 18/12/74 W.HART */
678 
679 
680  /* Parameter adjuTments */
681  --b; --a;
682 
683  /* Function Body */
684  ipiv = 0;
685 
686  i__ = 0;
687 
688  do {
689  ++i__;
690  ipiv += i__;
691  kpiv = ipiv;
692  r__ = a[ipiv];
693 
694  for (j = i__; j <= n; ++j) {
695  sum = 0.;
696  if (i__ == 1) goto L40;
697  if (r__ == 0.) goto L42;
698  id = ipiv - i__ + 1;
699  kd = kpiv - i__ + 1;
700 
701  do {
702  sum += b[kd] * b[id];
703  ++kd; ++id;
704  } while (id < ipiv);
705 
706 L40:
707  sum = a[kpiv] - sum;
708 L42:
709  if (j != i__) b[kpiv] = sum * r__;
710  else {
711  if (sum<ans) ans = sum;
712  if (sum<=0.) goto RETN;
713  dc = sqrt(sum);
714  b[kpiv] = dc;
715  if (r__ > 0.) r__ = (double)1. / dc;
716  }
717  kpiv += j;
718  }
719 
720  } while (i__ < n);
721 
722 RETN: aa[0]=save;
723  if (myb!=B) delete [] myb;
724  return ans;
725 } /* trchlu_ */
726 //______________________________________________________________________________
727 //______________________________________________________________________________
728 static double myDif(TVectorD &a,TVectorD &b, int &idx)
729 {
730 
731  int n = a.GetNrows();
732  double dif=0; idx = -1;
733  for (int i=0;i<n;i++) {
734  double mydif = fabs(a[i])+fabs(b[i]);
735  if (mydif<1e-4) continue;
736  mydif = fabs(a[i]-b[i])/mydif;
737  if (mydif <= dif) continue;
738  dif = mydif; idx = i;
739  }
740  return dif;
741 }
742 //______________________________________________________________________________
743 void THelix3d::Test()
744 {
745  double PtGev = 1.,Rho = 1./100, Hz = PtGev*Rho ;
746  double ZER[3]= {0,0,0};
747  double HZ[3] = {0,0,Hz};
748  double XZ[3] = {100,100,100};
749  double HH[3],XX[3],PP[3],PZ[3] ;
750 
751  TVector3 vHZ(HZ);
752  vHZ.RotateX(1.); vHZ.RotateY(1.); vHZ.GetXYZ(HH);
753 
754  TVector3 vXZ(XZ);
755  vXZ.RotateX(1.); vXZ.RotateY(1.); vXZ.GetXYZ(XX);
756 
757  TVector3 vPZ(PtGev,0,4./3*PtGev);
758  /*vPZ.RotateZ(1.);*/ vPZ.GetXYZ(PZ);
759  vPZ.RotateX(1.); vPZ.RotateY(1.); vPZ.GetXYZ(PP);
760 
761  for (int charge=1;charge >=-1;charge-=2) {
762 
763  THelixTrack_ TH(ZER,PZ,-Rho*charge);
764  double s1 = TH.Path(XZ);
765  THelix3d T3(charge,ZER,PP,HH);
766  double s2 = T3.Path(XX);
767  printf ("Charge = %d S1,S2 = %g %g \n",charge,s1,s2);
768  }
769 
770 }
771 //______________________________________________________________________________
772 void THelix3d::Test2()
773 {
774 printf("THelix3d::Test2() Empty\n");
775 }
776 //______________________________________________________________________________
777 void THelix3d::TestDer2()
778 {
779 // Test for derivatives + changing direction
780 
781  double Ptot = 5, Pt=3;
782 //double curv = 1./100,Hz = Pt*curv;
783  double curv = 1./30,Hz = Pt*curv;
784  TVector3 XbegV = TVector3(10,20,30);
785  TVector3 PbegV = TVector3(Pt,0,Pt/3*4);
786  double HZ[3] = {0,0,Hz};
787 
788 
789  double L = 100.;
790 
791  double Mag[3];
792  TVector3 magV(HZ); magV.GetXYZ(Mag );
793 #if 1
794  double r1=gRandom->Rndm(),r2=gRandom->Rndm();
795  magV.RotateX(r1); magV.RotateY(r2); magV.GetXYZ(Mag );
796  PbegV.RotateX(r1); PbegV.RotateY(r2);
797  XbegV.RotateX(r1); XbegV.RotateY(r2);
798 #endif
799 
800 // XbegV.Print("XbegV");
801 // PbegV.Print("PbegV");
802 // magV.Print("magV");
803 
804  printf("===============PbegV:= \n"); PbegV.Print("");
805 
806  TVector3 XendV,PendV;
807  const THDer3d_t *myDer = 0;
808  for (int ichar=-1; ichar<=1; ichar +=2) {
809  THelix3d TH0(ichar,(double*)&XbegV[0],(double*)&PbegV[0],Mag);
810  TH0.SetDerOn();
811  TH0.Move(L);
812  TH0.Backward(); //At the end TH0 helix inverted
813  myDer = TH0.Der(); //Inverted derivatives
814  TH0.Eval(0,&XendV[0],&PendV[0]); // Got X & P at the end
815  auto DendV = PendV.Unit();
816 
817  TVector3 UbegV(myDer->TkDir(0)[kKU]); //Phi vector
818  TVector3 VbegV(myDer->TkDir(0)[kKV]); //Lam vector
819  TVector3 TbegV(myDer->TkDir(0)[kKT]); //Along track vector
820 
821  TVector3 UendV(myDer->TkDir(1)[kKU]); //Phi vector
822  TVector3 VendV(myDer->TkDir(1)[kKV]); //Lam vector
823  TVector3 TendV(myDer->TkDir(1)[kKT]); //Along track vector
824  assert((DendV-TendV).Mag()<1e-3);
825 
826 
827  TVectorD der(5);
828  static const char* dp[5] = {"U","V","Fita","Lama","Pinv"};
829  for (int J=0;J<5;J++) {
830 
831  printf("\nTest dFit/d%s\n\n",dp[J]);
832  for (int i=0;i<5;i++) {der[i] = (*myDer)[i][J];}
833 
834 
835  THelix3d TH1(ichar,(double*)&XbegV[0],(double*)&PbegV[0],Mag);
836  TH1.Backward();
837  TVector3 P1begV(TH1.Mom());
838  TVector3 X1begV(TH1.Pos());
839  TkDir_t myTkDir;
840  THelix3d::MakeTkDir((double*)&P1begV[0],Mag,myTkDir);
841  double eps = (TVectorD(9,myTkDir[0])-TVectorD(9,myDer->TkDir(0)[0])).NormInf();
842  assert(eps<1e-3);
843 
844 
845 //?? double delta = 1e-4;
846  double delta = 1e-6;
847  TVector3 X1endV,P1endV;
848  assert(TH0.Pinv()*ichar>0);
849  assert(TH1.Pinv()*ichar>0);
850  switch(J) {
851  case kU:{
852  X1begV += UbegV*delta;
853  break;
854  }
855  case kV: {
856  X1begV += VbegV*delta;
857  break;
858  }
859  case kPinv: {
860  double Pinv1 = TH1.Pinv()+delta;
861  double Ptot1 = fabs(1./Pinv1);
862  P1begV.SetMag(Ptot1);
863  break;
864  }
865  case kLama: { //lamda variation (alfa)
866  P1begV+=VbegV*(Ptot*delta); P1begV.SetMag(Ptot);
867  break;
868  }
869 
870  case kFita: { //Phi variation (beta)
871  P1begV+=UbegV*(Ptot*delta); P1begV.SetMag(Ptot);
872  break;
873  }
874  default: assert(0 && "Wrong J");
875  }
876  TH1.Set(0,(double*)&X1begV[0],(double*)&P1begV[0],0);
877  assert(TH1.Pinv()*ichar>0);
878 
879  TH1.Move(-L);
880 
881  TH1.Eval(0,&X1endV[0],&P1endV[0]);
882  auto D1endV = P1endV.Unit();
883  double dL = -(X1endV-XendV).Dot(PendV)/(PendV.Dot(D1endV));
884  TH1.Move(dL);
885  TH1.Eval(0,&X1endV[0],&P1endV[0]);
886 
887  TVectorD dif(5);
888  TVector3 difX = X1endV-XendV;
889  TVector3 difD = (P1endV.Unit()-PendV.Unit());
890  dif[kU] = difX.Dot(UendV);
891  dif[kV] = difX.Dot(VendV);
892  dif[kPinv] = TH1.Pinv()-TH0.Pinv();
893  dif[kFita] = difD.Dot(UendV);
894  dif[kLama] = difD.Dot(VendV);
895 
896 
897  assert(fabs(dif.NormInf())<1.);
898  dif*=1./delta;
899 
900  int idx=-1;
901  double maxa = myDif(dif,der,idx);
902 
903  if (maxa>1e-1) printf("########################[%d]= %g\n",idx,maxa);
904  der.Print("ANALITIC");
905  dif.Print("NUMERIC ");
906  }
907  } //end charge loop
908 }
909 //______________________________________________________________________________
910 void THelix3d::TestDer()
911 {
912 // Test for derivatives + changing direction
913 
914  double Pt=1;
915  double curv = 1./30,Hz = Pt*curv;
916  TVector3 XbegV = TVector3(0,0,0);
917  TVector3 PbegV = TVector3(Pt,0,Pt/3*4);
918  double Ptot = PbegV.Mag();
919  double cosBeg = Pt/Ptot;
920 
921  double HZ[3] = {0,0,Hz};
922 
923 
924  double L = 100.;
925  L = M_PI/curv/cosBeg;
926 
927 
928 
929  double Mag[3];
930  TVector3 magV(HZ); magV.GetXYZ(Mag );
931 #if 1
932  double r1=gRandom->Rndm(),r2=gRandom->Rndm();
933  magV.RotateX(r1); magV.RotateY(r2); magV.GetXYZ(Mag );
934  PbegV.RotateX(r1); PbegV.RotateY(r2);
935  XbegV.RotateX(r1); XbegV.RotateY(r2);
936 #endif
937 
938 
939  printf("===============PbegV:= \n"); PbegV.Print("");
940 
941  TVector3 XendV,PendV;
942  const THDer3d_t *myDer = 0;
943  for (int ichar=-1; ichar<=1; ichar +=2) {
944  THelix3d TH0(ichar,(double*)&XbegV[0],(double*)&PbegV[0],Mag);
945  TH0.SetDerOn();
946  TH0.Move(L);
947  //?? TH0.Backward(); //At the end TH0 helix inverted
948  myDer = TH0.Der(); //Inverted derivatives
949  TH0.Eval(0,&XendV[0],&PendV[0]); // Got X & P at the end
950  auto DendV = PendV.Unit();
951 
952  TVector3 UbegV(myDer->TkDir(0)[kKU]); //Phi vector
953  TVector3 VbegV(myDer->TkDir(0)[kKV]); //Lam vector
954  TVector3 TbegV(myDer->TkDir(0)[kKT]); //Along track vector
955 
956  TVector3 UendV(myDer->TkDir(1)[kKU]); //Phi vector
957  TVector3 VendV(myDer->TkDir(1)[kKV]); //Lam vector
958  TVector3 TendV(myDer->TkDir(1)[kKT]); //Along track vector
959  assert((DendV-TendV).Mag()<1e-3);
960 
961 
962  TVectorD der(5);
963  static const char* dp[5] = {"U","V","Fita","Lama","Pinv"};
964  for (int J=0;J<5;J++) {
965 
966  printf("\nTest dFit/d%s\n",dp[J]);
967  printf("==============\n");
968  for (int i=0;i<5;i++) {der[i] = (*myDer)[i][J];}
969 
970 
971  THelix3d TH1(ichar,(double*)&XbegV[0],(double*)&PbegV[0],Mag);
972 //?? TH1.Backward();
973  TVector3 P1begV(TH1.Mom());
974  TVector3 X1begV(TH1.Pos());
975  TkDir_t myTkDir;
976  THelix3d::MakeTkDir((double*)&P1begV[0],Mag,myTkDir);
977  double eps = (TVectorD(9,myTkDir[0])-TVectorD(9,myDer->TkDir(0)[0])).NormInf();
978  assert(eps<1e-3);
979 
980 
981  double delta = 1e-6;
982  TVector3 X1endV,P1endV;
983  switch(J) {
984  case kU:{
985  X1begV += UbegV*delta;
986  break;
987  }
988  case kV: {
989  X1begV += VbegV*delta;
990  break;
991  }
992  case kPinv: {
993  double Pinv1 = TH1.Pinv()+delta;
994  double Ptot1 = fabs(1./Pinv1);
995  P1begV.SetMag(Ptot1);
996  break;
997  }
998  case kLama: { //lamda variation (alfa)
999  P1begV+=VbegV*(Ptot*delta); P1begV.SetMag(Ptot);
1000  break;
1001  }
1002 
1003  case kFita: { //Phi variation (beta)
1004  P1begV+=UbegV*(Ptot*delta); P1begV.SetMag(Ptot);
1005  break;
1006  }
1007  default: assert(0 && "Wrong J");
1008  }
1009  TH1.Set(0,(double*)&X1begV[0],(double*)&P1begV[0],0);
1010  TH1.Move(L);
1011  TH1.Eval(0,&X1endV[0],&P1endV[0]);
1012  auto D1endV = P1endV.Unit();
1013  double dL = -(X1endV-XendV).Dot(PendV)/(PendV.Dot(D1endV));
1014  TH1.Move(dL);
1015  TH1.Eval(0,&X1endV[0],&P1endV[0]);
1016 
1017  TVectorD dif(5);
1018  TVector3 difX = X1endV-XendV;
1019  TVector3 difD = (P1endV.Unit()-PendV.Unit());
1020  dif[kU] = difX.Dot(UendV);
1021  dif[kV] = difX.Dot(VendV);
1022  dif[kPinv] = TH1.Pinv()-TH0.Pinv();
1023  dif[kFita] = difD.Dot(UendV);
1024  dif[kLama] = difD.Dot(VendV);
1025 
1026 
1027  assert(fabs(dif.NormInf())<1.);
1028  dif*=1./delta;
1029 
1030  int idx=-1;
1031  double maxa = myDif(dif,der,idx);
1032 
1033  if (maxa>1e-1) printf("########################[%d]= %g\n",idx,maxa);
1034  der.Print("ANALITIC");
1035  dif.Print("NUMERIC ");
1036  }
1037  } //end charge loop
1038 }
1039 //______________________________________________________________________________
1040 void THelix3d::TestErr(int charge)
1041 {
1042  double PtGev = 3.,Curv = 1./100, Hz = PtGev*Curv;
1043  double HZ[3] = {0,0,Hz};
1044  double Pbeg[3],P1beg[3],Xbeg[3]={0},X1beg[3];
1045  double Pend[3],P1end[3],Xend[3] ,X1end[3];
1046 
1047  double L = 300.;
1048 
1049 // Canvas + Histograms
1050 enum {kNCanvs=10};
1051 static TCanvas *myCanvas[kNCanvs] = {0};
1052 static int myDiv[] = {2,5,5,5,5,0};
1053  for (int ic = 0;myDiv[ic];ic++) {
1054  if (myDiv[ic]<0) continue;
1055  TString ts("C"); ts+=ic;
1056  if (!myCanvas[ic]) myCanvas[ic] = new TCanvas(ts.Data(),ts.Data(),600,800);
1057  myCanvas[ic]->Clear();myCanvas[ic]->Divide(1,myDiv[ic]);
1058  }
1059 static TH1F *hXi2[2] = {0};
1060  for (int jk=0;jk<2;jk++) {
1061  TString ts("Xi2_"); ts += jk;
1062  delete hXi2[jk];
1063  hXi2[jk]= new TH1F(ts.Data(),ts.Data(),50,0,0);
1064  myCanvas[0]->cd(jk+1); hXi2[jk]->Draw();
1065  }
1066 // Now test of DCA
1067 static TH1F * hh[5]={0};
1068  for (int ih=0;ih<myDiv[1];ih++) {
1069 const char * tit[]={"U","V","Fita","Lama","Pinv",0};
1070  delete hh[ih];
1071  hh[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1072  myCanvas[1]->cd(ih+1); hh[ih]->Draw();
1073  }
1074 // Now test of TRandomVector
1075 static TH1F * hrt[5]={0};
1076  for (int ih=0;ih<myDiv[1];ih++) {
1077 const char * tit[]={"U0","V0","Pinv0","Fita0","Lama0",0};
1078  delete hrt[ih];
1079  hrt[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1080  myCanvas[2]->cd(ih+1); hrt[ih]->Draw();
1081  }
1082 // Now test of Corr
1083 static TH1F * hcr[10]={0};
1084 {
1085 static const char *tit[]={"UV","UF","VF","UL","VL","FL","UP","VP","FP","LP"};
1086  for (int ih=0;ih<10;ih++) {
1087  delete hcr[ih];
1088  hcr[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1089  if (ih<5) myCanvas[3]->cd(ih+1 );
1090  else myCanvas[4]->cd(ih+1-5);
1091  hcr[ih]->Draw();
1092  }}
1093 
1094 
1095  TVector3 PbegV(PtGev,0,PtGev*3/4);
1096  double Ptot = PbegV.Mag();
1097  PbegV.RotateZ(gRandom->Rndm()*3.1415*2); PbegV.GetXYZ(Pbeg);
1098  TVector3 XbegV(Xbeg);
1099 
1100  double Mag[3]={HZ[0],HZ[1],HZ[2]};
1101 #if 11
1102  TVector3 magV(HZ);
1103  double r1 = gRandom->Rndm(),r2 = gRandom->Rndm();
1104  magV.RotateX(r1); magV.RotateY(r2); magV.GetXYZ(Mag);
1105  PbegV.RotateX(r1); PbegV.RotateY(r2); PbegV.GetXYZ(Pbeg);
1106 #endif
1107  TVectorD dia(5);
1108  dia[kU]= 0.1; dia[kV]= 0.2; dia[kPinv]= 0.1/Ptot; dia[kFita]= 1./360; dia[kLama]= 2./360;
1109  for(int jk=0;jk<5;jk++) {dia[jk]*=dia[jk];}
1110 
1111  TRandomVector RV(dia);
1112  auto &EMX = RV.GetMtx();
1113  auto &val = RV.GetLam();
1114  dia.Print("DIA");
1115  val.Print("VAL");
1116 
1117 
1118  double Gbeg[15],GbegD[5];
1119  for (int i=0,li=0;i< 5;li+=++i) {
1120  GbegD[i] = sqrt(EMX[i][i]);
1121  for (int j=0;j<=i;j++) {
1122  Gbeg[li+j] = EMX[i][j];
1123  } }
1124  assert( mySign(Gbeg,5)>0);
1125 
1126 
1127  double GbegI[15];
1128  TCL::trsinv(Gbeg,GbegI,5);
1129  assert( mySign(GbegI,5)>0);
1130 
1131  THelix3d TH0(charge,Xbeg,Pbeg,Mag);
1132  THEmx3d_t *emx = new THEmx3d_t(Gbeg);
1133  TH0.SetEmx(emx);
1134  TH0.Move(L);
1135  TH0.Eval(0,Xend,Pend);
1136  auto *myEmx = TH0.Emx();
1137  const double *Gend = *myEmx;
1138  assert(mySign(Gend,5)>0);
1139  double GendD[5];
1140  for (int i=0,li=0;i< 5;li+=++i) {
1141  GendD[i] = sqrt(Gend[li+i]);
1142  }
1143  double GendI[15];
1144  TCL::trsinv(Gend,GendI,5);
1145 
1146  assert(mySign(GendI,5)>0);
1147 
1148  auto &tkDir = TH0.TkDir(0);
1149  auto &tkDirE = TH0.TkDir(1);
1150  TVector3 PendV(Pend),XendV(Xend);
1151  TVectorD UendV(5);
1152  for (int j=0;j<3;j++) { UendV[j+2]=Pend[j];}
1153 
1154  double myG[15]={0};
1155  int nIter = 1000000;
1156  for (int iter=0;iter<nIter;iter++) {
1157  TVectorD delta = RV.Gaus();
1158  double chi2;
1159  TCL::trasat(delta.GetMatrixArray(),GbegI,&chi2,1,5);
1160  assert(chi2>0);
1161  hXi2[0]->Fill(chi2);
1162 
1163  for (int ih=0;ih<myDiv[1];ih++) { hrt[ih]->Fill(delta[ih]/GbegD[ih]);}
1164  TVector3 P1begV = PbegV;
1165  TVector3 X1begV = XbegV;
1166 
1167  X1begV += TVector3(tkDir[kKU])*delta[kU];
1168  X1begV += TVector3(tkDir[kKV])*delta[kV];
1169  double Pinv = -charge/Ptot;
1170  double Pinv1 = Pinv + delta[kPinv];
1171  double Ptot1 = fabs(1./Pinv1);
1172  P1begV.SetMag(Ptot1);
1173  P1begV += TVector3(tkDir[kKV])*(delta[kLama]*Ptot1);P1begV.SetMag(Ptot1);
1174  P1begV += TVector3(tkDir[kKU])*(delta[kFita]*Ptot1);P1begV.SetMag(Ptot1);
1175 
1176  X1begV.GetXYZ(X1beg);
1177  P1begV.GetXYZ(P1beg);
1178  THelix3d TH1(charge,X1beg,P1beg,Mag);
1179  double s = TH1.Path(Xend);
1180  TH1.Move(s);
1181  TH1.Eval(0,X1end,P1end);
1182  TVector3 X1endV(X1end),P1endV(P1end);
1183  TVectorD U1endV(5);
1184  TVector3 difX = X1endV-XendV;
1185  TVector3 difD = (P1endV.Unit()-PendV.Unit());
1186 
1187  U1endV[kU] = difX.Dot(TVector3(tkDirE[kKU]));
1188  U1endV[kV] = difX.Dot(TVector3(tkDirE[kKV]));
1189  U1endV[kPinv] = (Pinv1-Pinv);
1190  U1endV[kLama] = difD.Dot(TVector3(tkDirE[kKV]));
1191  U1endV[kFita] = difD.Dot(TVector3(tkDirE[kKU]));
1192 
1193  for (int ih=0;ih<5;ih++) { hh[ih]->Fill(U1endV[ih]/GendD[ih]);}
1194  TCL::trasat(U1endV.GetMatrixArray(),GendI,&chi2,1,5);
1195  hXi2[1]->Fill(chi2);
1196 
1197  for (int i=0,li=0;i<5;li+=++i) {
1198  for (int j=0;j<=i;j++) {
1199  myG[li+j]+=U1endV[i]*U1endV[j];
1200  } }
1201 
1202  int jk=0;
1203  for (int i=0,li=0;i<5;li+=++i) {
1204  for (int j=0;j<i;j++) {
1205  double d = U1endV[i]*U1endV[j] - Gend[li+j];
1206  d/= GendD[i]*GendD[j];
1207  hcr[jk++]->Fill(d);
1208  } }
1209  }
1210 
1211  TCL::vscale(myG,1./nIter,myG,15);
1212 
1213  printf("Numerical vs Analitical Error matrices:\n");
1214  for (int i=0,li=0;i<5;li+=++i) {
1215  for (int j=0;j<=i;j++) {
1216  printf("\t%g ",myG[li+j]);
1217  }
1218  printf("\n");
1219  for (int j=0;j<=i;j++) {
1220  printf("\t%g ",Gend[li+j]);
1221  }
1222  printf("\n");
1223  }
1224 
1225  for (int i=0,li=0;i<5;li+=++i) {
1226  for (int j=0;j<=i;j++) {
1227  printf("\t%g ",Gend[li+j]/(GendD[i]*GendD[j]));
1228  }
1229  printf("\n");
1230  }
1231 
1232  for (int i=0;myCanvas[i];i++) {
1233  if (!myCanvas[i]) continue;
1234  myCanvas[i]->Modified();myCanvas[i]->Update();}
1235 
1236  while(!gSystem->ProcessEvents()){gSystem->Sleep(200);};
1237 }
1238 //______________________________________________________________________________
1239 void THelix3d::TestErr2(int charge)
1240 {
1241  double PtGev = 1.,Curv = 1./30, Hz = PtGev*Curv;
1242  double HZ[3] = {0,0,Hz};
1243  double Pbeg[3],P1beg[3],Xbeg[3]={0},X1beg[3];
1244  double Pend[3],P1end[3],Xend[3] ,X1end[3];
1245 
1246 // Canvas + Histograms
1247 enum {kNCanvs=10};
1248 static TCanvas *myCanvas[kNCanvs] = {0};
1249 static int myDiv[] = {2,5,5,5,5,0};
1250  for (int ic = 0;myDiv[ic];ic++) {
1251  if (myDiv[ic]<0) continue;
1252  TString ts("C"); ts+=ic;
1253  if (!myCanvas[ic]) myCanvas[ic] = new TCanvas(ts.Data(),ts.Data(),600,800);
1254  myCanvas[ic]->Clear();myCanvas[ic]->Divide(1,myDiv[ic]);
1255  }
1256 static TH1F *hXi2[2] = {0};
1257  for (int jk=0;jk<2;jk++) {
1258  TString ts("Xi2_"); ts += jk;
1259  delete hXi2[jk];
1260  hXi2[jk]= new TH1F(ts.Data(),ts.Data(),50,0,0);
1261  myCanvas[0]->cd(jk+1); hXi2[jk]->Draw();
1262  }
1263 // Now test of DCA
1264 static TH1F * hh[5]={0};
1265  for (int ih=0;ih<myDiv[1];ih++) {
1266 const char * tit[]={"U","V","Pinv","Fita","Lama",0};
1267  delete hh[ih];
1268  hh[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1269  myCanvas[1]->cd(ih+1); hh[ih]->Draw();
1270  }
1271 // Now test of TRandomVector
1272 static TH1F * hrt[5]={0};
1273  for (int ih=0;ih<myDiv[1];ih++) {
1274 const char * tit[]={"U0","V0","Pinv0","Fita0","Lama0",0};
1275  delete hrt[ih];
1276  hrt[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1277  myCanvas[2]->cd(ih+1); hrt[ih]->Draw();
1278  }
1279 // Now test of Corr
1280 static TH1F * hcr[10]={0};
1281 {
1282 static const char *tit[]={"UV","UP","VP","UL","VL","PL","UF","VF","PF","LF"};
1283  for (int ih=0;ih<10;ih++) {
1284  delete hcr[ih];
1285  hcr[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1286  if (ih<5) myCanvas[3]->cd(ih+1 );
1287  else myCanvas[4]->cd(ih+1-5);
1288  hcr[ih]->Draw();
1289  }}
1290 //============================================================================
1291 
1292  TVector3 PbegV(PtGev,0,PtGev*3);
1293  double Ptot = PbegV.Mag();
1294  PbegV.RotateZ(gRandom->Rndm()*3.1415*2); PbegV.GetXYZ(Pbeg);
1295  TVector3 XbegV(Xbeg);
1296 
1297  double L = 100.;
1298 
1299  double Mag[3]={HZ[0],HZ[1],HZ[2]};
1300 #if 11
1301  TVector3 magV(HZ);
1302  double r1 = gRandom->Rndm(),r2 = gRandom->Rndm();
1303  magV.RotateX(r1); magV.RotateY(r2); magV.GetXYZ(Mag);
1304  PbegV.RotateX(r1); PbegV.RotateY(r2); PbegV.GetXYZ(Pbeg);
1305 #endif
1306  TVectorD dia(5);
1307  dia[kU]= 0.1; dia[kV]= 0.2; dia[kPinv]= 0.1/Ptot; dia[kFita]= 1./360; dia[kLama]= 2./360;
1308 
1309  TRandomVector RV(dia);
1310  auto &EMX = RV.GetMtx();
1311  auto &val = RV.GetLam();
1312  dia.Print("DIA");
1313  val.Print("VAL");
1314 
1315 
1316  double Gbeg[15],GbegD[5];
1317  for (int i=0,li=0;i< 5;li+=++i) {
1318  GbegD[i] = sqrt(EMX[i][i]);
1319  for (int j=0;j<=i;j++) {
1320  Gbeg[li+j] = EMX[i][j];
1321  } }
1322  assert( mySign(Gbeg,5)>0);
1323 
1324 
1325 
1326  THelix3d TH0(charge,Xbeg,Pbeg,Mag);
1327 //==================================
1328  THEmx3d_t *emx = new THEmx3d_t(Gbeg);
1329 
1330  THelix3d TH0i(&TH0);
1331  TH0i.SetEmx(emx);
1332  TH0i.Backward();
1333  double GbegI[15];
1334  memcpy(Gbeg,*TH0i.Emx(),sizeof(Gbeg));
1335  TRandomVector RVi(5,Gbeg);
1336  TCL::trsinv(Gbeg,GbegI,5);
1337  assert( mySign(GbegI,5)>0);
1338 
1339  TH0.SetEmx(emx);
1340  TH0.Move(L);
1341  TH0.Backward();
1342 
1343  double tkDir[3][3],tkDirE[3][3];
1344  memcpy(tkDir[0] ,TH0.TkDir(0)[0],sizeof(tkDir ));
1345  memcpy(tkDirE[0],TH0.TkDir(1)[0],sizeof(tkDirE));
1346 
1347  TH0.Eval(0,Xend,Pend);
1348  auto *myEmx = TH0.Emx();
1349  const double *Gend = *myEmx;
1350  assert(mySign(Gend,5)>0);
1351  double GendD[5];
1352  for (int i=0,li=0;i< 5;li+=++i) {
1353  GendD[i] = sqrt(Gend[li+i]);
1354  }
1355 
1356  double GendI[15];
1357  TCL::trsinv(Gend,GendI,5);
1358 
1359  TVector3 PendV(Pend),XendV(Xend);
1360  TVectorD UendV(5);
1361  for (int j=0;j<3;j++) { UendV[j+2]=Pend[j];}
1362 
1363  double myG[15]={0};
1364  int nIter = 1000000;
1365  for (int iter=0;iter<nIter;iter++) {
1366  TVectorD delta = RVi.Gaus();
1367  double chi2;
1368  TCL::trasat(delta.GetMatrixArray(),GbegI,&chi2,1,5);
1369  assert(chi2>0);
1370  hXi2[0]->Fill(chi2);
1371 
1372  for (int ih=0;ih<myDiv[1];ih++) { hrt[ih]->Fill(delta[ih]/GbegD[ih]);}
1373  TVector3 P1begV(TH0i.Mom());
1374  TVector3 X1begV(TH0i.Pos());
1375 
1376  X1begV += TVector3(tkDir[kKU])*delta[kU];
1377  X1begV += TVector3(tkDir[kKV])*delta[kV];
1378  double Pinv = TH0i.Pinv();
1379  double Pinv1 = Pinv + delta[kPinv];
1380  double Ptot1 = fabs(1./Pinv1);
1381  P1begV.SetMag(Ptot1);
1382  P1begV += TVector3(tkDir[kKV])*(delta[kLama]*Ptot1);P1begV.SetMag(Ptot1);
1383  P1begV += TVector3(tkDir[kKU])*(delta[kFita]*Ptot1);P1begV.SetMag(Ptot1);
1384 
1385  X1begV.GetXYZ(X1beg);
1386  P1begV.GetXYZ(P1beg);
1387  THelix3d TH1(TH0i.Charge(),X1beg,P1beg,Mag);
1388  TH1.Move(-L);
1389  TH1.Eval(0,X1end,P1end);
1390  TVector3 X1endV(X1end),P1endV(P1end);
1391  TVectorD U1endV(5);
1392  TVector3 difX = X1endV-XendV;
1393  TVector3 difD = (P1endV.Unit()-PendV.Unit());
1394 
1395  U1endV[kU] = difX.Dot(TVector3(tkDirE[kKU]));
1396  U1endV[kV] = difX.Dot(TVector3(tkDirE[kKV]));
1397  U1endV[kFita] = difD.Dot(TVector3(tkDirE[kKU]));
1398  U1endV[kLama] = difD.Dot(TVector3(tkDirE[kKV]));
1399  U1endV[kPinv] = (Pinv1-Pinv);
1400 
1401  for (int ih=0;ih<5;ih++) { hh[ih]->Fill(U1endV[ih]/GendD[ih]);}
1402  TCL::trasat(U1endV.GetMatrixArray(),GendI,&chi2,1,5);
1403  hXi2[1]->Fill(chi2);
1404 
1405  for (int i=0,li=0;i<5;li+=++i) {
1406  for (int j=0;j<=i;j++) {
1407  myG[li+j]+=U1endV[i]*U1endV[j];
1408  } }
1409 
1410  int jk=0;
1411  for (int i=0,li=0;i<5;li+=++i) {
1412  for (int j=0;j<i;j++) {
1413  double d = U1endV[i]*U1endV[j] - Gend[li+j];
1414  d/= GendD[i]*GendD[j];
1415  hcr[jk++]->Fill(d);
1416  } }
1417  }
1418 
1419  TCL::vscale(myG,1./nIter,myG,15);
1420 
1421  printf("Numerical vs Analitical Error matrices:\n");
1422  for (int i=0,li=0;i<5;li+=++i) {
1423  for (int j=0;j<=i;j++) {
1424  printf("\t%g ",myG[li+j]);
1425  }
1426  printf("\n");
1427  for (int j=0;j<=i;j++) {
1428  printf("\t%g ",Gend[li+j]);
1429  }
1430  printf("\n");
1431  }
1432 
1433 
1434  for (int i=0;myCanvas[i];i++) {
1435  if (!myCanvas[i]) continue;
1436  myCanvas[i]->Modified();myCanvas[i]->Update();}
1437 
1438  while(!gSystem->ProcessEvents()){gSystem->Sleep(200);};
1439 }
1440 //______________________________________________________________________________
1441 void THelix3d::TestErr3(int charge)
1442 {
1443 // double PtGev = 1.,Curv = 1./100, Hz = PtGev*Curv;
1444  double PtGev = 1.,Curv = 1./30, Hz = PtGev*Curv;
1445  double HZ[3] = {0,0,Hz};
1446  double Pbeg[3],P1beg[3],Xbeg[3]={0},X1beg[3];
1447  double Pend[3],P1end[3],Xend[3] ,X1end[3];
1448  double myTkRot = 0.1;
1449 
1450 // Canvas + Histograms
1451 enum {kNCanvs=10};
1452 static TCanvas *myCanvas[kNCanvs] = {0};
1453 static int myDiv[] = {2,5,5,5,5,0};
1454  for (int ic = 0;myDiv[ic];ic++) {
1455  if (myDiv[ic]<0) continue;
1456  TString ts("C"); ts+=ic;
1457  if (!myCanvas[ic]) myCanvas[ic] = new TCanvas(ts.Data(),ts.Data(),600,800);
1458  myCanvas[ic]->Clear();myCanvas[ic]->Divide(1,myDiv[ic]);
1459  }
1460 static TH1F *hXi2[2] = {0};
1461  for (int jk=0;jk<2;jk++) {
1462  TString ts("Xi2_"); ts += jk;
1463  delete hXi2[jk];
1464  hXi2[jk]= new TH1F(ts.Data(),ts.Data(),50,0,0);
1465  myCanvas[0]->cd(jk+1); hXi2[jk]->Draw();
1466  }
1467 // Now test of DCA
1468 static TH1F * hh[5]={0};
1469  for (int ih=0;ih<myDiv[1];ih++) {
1470 const char * tit[]={"U","V","Fita","Lama","Pinv",0};
1471  delete hh[ih];
1472  hh[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1473  myCanvas[1]->cd(ih+1); hh[ih]->Draw();
1474  }
1475 // Now test of TRandomVector
1476 static TH1F * hrt[5]={0};
1477  for (int ih=0;ih<myDiv[1];ih++) {
1478 const char * tit[]={"U0","V0","Pinv0","Fita0","Lama0",0};
1479  delete hrt[ih];
1480  hrt[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1481  myCanvas[2]->cd(ih+1); hrt[ih]->Draw();
1482  }
1483 // Now test of Corr
1484 static TH1F * hcr[10]={0};
1485 {
1486 static const char *tit[]={"UV","UF","VF","UL","VL","FL","UP","VP","FP","LP"};
1487  for (int ih=0;ih<10;ih++) {
1488  delete hcr[ih];
1489  hcr[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1490  if (ih<5) myCanvas[3]->cd(ih+1 );
1491  else myCanvas[4]->cd(ih+1-5);
1492  hcr[ih]->Draw();
1493  }}
1494 
1495 
1496  TVector3 PbegV(PtGev,0,PtGev*3);
1497  double Ptot = PbegV.Mag();
1498  PbegV.RotateZ(gRandom->Rndm()*3.1415*2); PbegV.GetXYZ(Pbeg);
1499  TVector3 XbegV(Xbeg);
1500 
1501  double Mag[3]={HZ[0],HZ[1],HZ[2]};
1502 #if 11
1503  TVector3 magV(HZ);
1504  double r1 = gRandom->Rndm(),r2 = gRandom->Rndm();
1505  magV.RotateX(r1); magV.RotateY(r2); magV.GetXYZ(Mag);
1506  PbegV.RotateX(r1); PbegV.RotateY(r2); PbegV.GetXYZ(Pbeg);
1507 #endif
1508  TVectorD dia(5);
1509  dia[kU]= 0.1; dia[kV]= 0.2; dia[kPinv]= 0.1/Ptot; dia[kFita]= 1./360; dia[kLama]= 2./360;
1510 
1511  TRandomVector RV(dia);
1512  auto &EMX = RV.GetMtx();
1513  auto &val = RV.GetLam();
1514  dia.Print("DIA");
1515  val.Print("VAL");
1516 
1517 
1518  double Gbeg[15],GbegD[5];
1519  for (int i=0,li=0;i< 5;li+=++i) {
1520  GbegD[i] = sqrt(EMX[i][i]);
1521  for (int j=0;j<=i;j++) {
1522  Gbeg[li+j] = EMX[i][j];
1523  } }
1524  assert( mySign(Gbeg,5)>0);
1525 
1526 
1527  double GbegI[15];
1528  TCL::trsinv(Gbeg,GbegI,5);
1529  assert( mySign(GbegI,5)>0);
1530 
1531 // initial state
1532  THelix3d TH0(charge,Xbeg,Pbeg,Mag);
1533  TH0.SetEmx(Gbeg);
1534  THEmx3d_t emx = *TH0.Emx();
1535  TkDir_t &tkDir = emx.TkDir();
1536 
1537 // now forget Pbeg, it was only temporary. Choose another Pbeg but close to
1538 // previous. But errors and TkDir will belong to new Pbeg. For this case
1539 // tkDir[kkT] != new Pbeg
1540 
1541 
1542 
1543  TH0.Eval(0,Xbeg,Pbeg);
1544  PbegV.RotateX(myTkRot);
1545  PbegV.RotateY(myTkRot);
1546  PbegV.GetXYZ(Pbeg);
1547  TH0.Set(0,0,Pbeg);
1548  TkDir_t tkDirE = TH0.TkDir(0);
1549 // Pend == Pbeg
1550  TVector3 PendV(Pbeg),XendV(Xbeg);
1551  PendV.GetXYZ(Pend); XendV.GetXYZ(Xend);
1552 
1553 
1554 
1555  auto *emxE = TH0.Emx();
1556  const double *Gend = *emxE;
1557 
1558 
1559  assert(mySign(Gend,5)>0);
1560  double GendD[5];
1561  for (int i=0,li=0;i< 5;li+=++i) {
1562  GendD[i] = sqrt(Gend[li+i]);
1563  }
1564  double GendI[15];
1565  TCL::trsinv(Gend,GendI,5);
1566 
1567  assert(mySign(GendI,5)>0);
1568  TVectorD UendV(5);
1569 
1570  for (int j=0;j<3;j++) { UendV[j+2]=Pend[j];}
1571 
1572  double myG[15]={0};
1573  int nIter = 1000000;
1574  for (int iter=0;iter<nIter;iter++) {
1575  TVectorD delta = RV.Gaus();
1576  double chi2;
1577  TCL::trasat(delta.GetMatrixArray(),GbegI,&chi2,1,5);
1578  assert(chi2>0);
1579  assert(chi2<100);
1580  hXi2[0]->Fill(chi2);
1581 
1582  for (int ih=0;ih<myDiv[1];ih++) { hrt[ih]->Fill(delta[ih]/GbegD[ih]);}
1583  TVector3 P1begV = PbegV;
1584  TVector3 D1begV = P1begV.Unit();
1585  TVector3 X1begV = XbegV;
1586  double cosL = D1begV.Dot(tkDir.V());
1587  double cosLF= D1begV.Dot(tkDir.T());
1588 
1589  cosL=1;cosLF=1;
1590 
1591  X1begV += TVector3(tkDir.U())*delta[kU];
1592  X1begV += TVector3(tkDir.V())*delta[kV];
1593  double Pinv = -charge/Ptot;
1594  double Pinv1 = Pinv + delta[kPinv];
1595  double Ptot1 = fabs(1./Pinv1);
1596  P1begV.SetMag(Ptot1);
1597  P1begV += TVector3(tkDir[kKV])*(delta[kLama]*Ptot1*cosL );P1begV.SetMag(Ptot1);
1598  P1begV += TVector3(tkDir[kKU])*(delta[kFita]*Ptot1*cosLF);P1begV.SetMag(Ptot1);
1599 
1600  X1begV.GetXYZ(X1beg);
1601  P1begV.GetXYZ(P1beg);
1602  THelix3d TH1(charge,X1beg,P1beg,Mag);
1603  double aa = (X1begV-XbegV).Dot(tkDirE.T());
1604  double bb = P1begV.Unit().Dot(tkDirE.T());
1605  double tau = -aa/bb;
1606  assert(fabs(tau)<1);
1607  TH1.Move(tau);
1608  TH1.Eval(0,X1end,P1end);
1609  TVector3 X1endV(X1end),P1endV(P1end);
1610  TVectorD U1endV(5);
1611  TVector3 difX = X1endV-XendV;
1612  TVector3 difD = (P1endV.Unit()-PendV.Unit());
1613 
1614  U1endV[kU] = difX.Dot(tkDirE[kKU]);
1615  U1endV[kV] = difX.Dot(tkDirE[kKV]);
1616  U1endV[kPinv] = (Pinv1-Pinv);
1617  U1endV[kLama] = difD.Dot(tkDirE[kKV]);
1618  U1endV[kFita] = difD.Dot(tkDirE[kKU]);
1619 
1620  for (int ih=0;ih<5;ih++) { hh[ih]->Fill(U1endV[ih]/GendD[ih]);}
1621  TCL::trasat(U1endV.GetMatrixArray(),GendI,&chi2,1,5);
1622  hXi2[1]->Fill(chi2);
1623 
1624  for (int i=0,li=0;i<5;li+=++i) {
1625  for (int j=0;j<=i;j++) {
1626  myG[li+j]+=U1endV[i]*U1endV[j];
1627  } }
1628 
1629  int jk=0;
1630  for (int i=0,li=0;i<5;li+=++i) {
1631  for (int j=0;j<i;j++) {
1632  double d = U1endV[i]*U1endV[j] - Gend[li+j];
1633  d/= GendD[i]*GendD[j];
1634  hcr[jk++]->Fill(d);
1635  } }
1636  }
1637 
1638  TCL::vscale(myG,1./nIter,myG,15);
1639 
1640  printf("Numerical vs Analitical Error matrices:\n");
1641  for (int i=0,li=0;i<5;li+=++i) {
1642  for (int j=0;j<=i;j++) {
1643  printf("\t%g ",myG[li+j]);
1644  }
1645  printf("\n");
1646  for (int j=0;j<=i;j++) {
1647  printf("\t%g ",Gend[li+j]);
1648  }
1649  printf("\n");
1650  }
1651 
1652  for (int i=0,li=0;i<5;li+=++i) {
1653  for (int j=0;j<=i;j++) {
1654  printf("\t%g ",Gend[li+j]/(GendD[i]*GendD[j]));
1655  }
1656  printf("\n");
1657  }
1658 
1659  for (int i=0;myCanvas[i];i++) {
1660  if (!myCanvas[i]) continue;
1661  myCanvas[i]->Modified();myCanvas[i]->Update();}
1662 
1663  while(!gSystem->ProcessEvents()){gSystem->Sleep(200);};
1664 }
1665 //______________________________________________________________________________
1666 void THelix3d::ConvertErrs(const THelixTrack_ *he, const double Bzp
1667  , double G[15],double D[5][5],double Di[5][5])
1668 {
1669 // in THelixTrack_:
1670 // Let Dx,Dy,Dz direction of track
1671 // dH: along vector (-Dy,Dx,0)
1672 // dZ: along Zaxis
1673 // dA: delta azimuth angle (in X,Y plane;
1674 // dC: delta of curvature;
1675 // dL = dLambda, angle between track and X,Y plane
1676 // B: direction of mag field and it is along Z
1677 // T = (cosLam*cosPhi,coaLam*sinPhi,sinLam) == track dir
1678 // H = (-sinPhi,cosPhi,0)
1679 // Z = ( 0, 0,1)
1680 // C = curvature
1681 // THelixTrack_ space point modification = H*h + Z*z
1682 //
1683 // // in THelix3d:
1684 // V: vector orthogonal T and in plane (T,B) in our case B==Z
1685 //
1686 // V ~= Z-(Z*T)*T = Z - sinLam*T
1687 // V ~= (-sinLam*coaLam*cosPhi,-sinLam*coaLam*sinPhi,1-sinLam**2)
1688 // V ~= (-sinLam*cosPhi,-sinLam*sinPhi,cosLam)*cosLam
1689 // After normalization:
1690 //
1691 //
1692 // V = (-sinLam*cosPhi,-sinLam*sinPhi,cosLam)
1693 // U = (- sinPhi, cosPhi, 0)
1694 //
1695 // We see that H == U
1696 //
1697 // THelixTrack_: space point modification = H*h + Z*z
1698 // THelix3d: space point modification = U*u + V*v
1699 //
1700 //
1701 // U*h + Z*z +T*t = U*u + V*v
1702 // t is moving along track from crossing with THelixTrack_ plane
1703 // orthogonal (Tx,Ty,0)
1704 // to THelix3d plane orthogonal (Tx,Ty,Tz)
1705 // Mult by T
1706 // (T*Z)z + t = 0
1707 // t = - (T*Z)*z
1708 //
1709 // U*h + Z*z -T*(T*Z)*z = U*u + V*v
1710 // U*h + (Z-T*(T*Z))*z = U*u + V*v
1711 //
1712 // From the above (Z-T*(T*Z)) == V*cosLam
1713 // U*h + V*cosLam*z = U*u + V*v
1714 //
1715 // Hence h == u and cosLam*z = v
1716 //
1717 // du/dh = 1
1718 // dv/dz = cosLam
1719 //
1720 
1721 
1722 // ======================================================================
1723 // Now angles:
1724 //
1725 // U*h + Z*z +T*t = U*u + V*v
1726 // T = Ex*cosLam*cosPhi+Ey*cosLam*sinPhi+Ez*sinLam
1727 //
1728 // dT = (-Ex*cosLam*sinPhi+Ey*cosLam*cosPhi )*dPhi
1729 // + (-Ex*sinLam*cosPhi-Ey*sinLam*sinPhi +Ez*CosLam)*dLam
1730 //
1731 // dT = (-Ex*sinPhi+Ey*cosPhi)*cosLam *dPhi
1732 // + ((-Ex*cosPhi-Ey*sinPhi)*sinLam +Ez*CosLam)*dLam
1733 // = U*fita +V*Lama
1734 //
1735 // dT = U*cosLam*dPhi+ V*dLam
1736 // fita = cosLam*dPhi
1737 // lama = dLam
1738 
1739 
1740 // dPhi = dPhi0 - (T*Z)*z*cosLam*Rho
1741 // dPhi = dPhi0 - sinLam*cosLam*Rho*z
1742 //
1743 // dT = U*cosLam*(dPhi-(T*Z)*z*cosLam*Rho) + V*dLam
1744 // = U*fita + V*Lama
1745 //
1746 // dFita/dPhi = cosLam
1747 // dFita/dz = - sinLam*cosLam*cosLam*Rho
1748 // dLamadLam = 1
1749 //
1750 //
1751 // Now Pinv
1752 // Pinv = -fCharge/fMom;
1753 // Rho = -fCharge*B[2]/pt;
1754 // -fCharge/pt = Rho/B[2]
1755 // -fCharge/(pt/cosLam) = Rho/B[2]*cosLam
1756 // -fCharge/(P) = Rho/B[2]*cosLam
1757 // Pinv = Rho*cosLam/B[2]
1758 //
1759 // dPinv/dRho = cosLam/B[2]
1760 // dPinv/dLam =-sinLam*Rho/B[2]
1761 //
1762 //=======================================================================
1763 // In total
1764 //=======================================================================
1765 // du = dH
1766 // dv = cosLam*dz
1767 // dFita = cosLam*dPhi - sinLam*cosLam*cosLam*Rho*dz
1768 // dLama = dLam
1769 // dPinv = cosLam/B[2]*dRho-sinLam*Rho/B[2]*dLam
1770 //=======================================================================
1771 // du/dh = 1
1772 // dv/dz = cosLam
1773 // dFita/dPhi = cosLam
1774 // dFita/dz = - sinLam*cosLam *cosLam*Rho
1775 // dLamadLam = 1
1776 // dPinv/dRho = cosLam/B[2]
1777 // dPinv/dLam =-sinLam*Rho/B[2]
1778 //
1779 // dH = du
1780 // dPhi = dFita/cosLam + sinLam*Rho*dv
1781 // dRho = dPinv/cosLam*B[2] +tanLam*Rho*dLama
1782 // dz = dv/cosLam
1783 // dLam = dLama
1784 //
1785 // ======================================================================
1786 static const double kMinBz = 0.001494399 * 1e-3;
1787  const double *T = he->Dir();
1788  double Rho = he->GetRho();
1789  double sinLam = T[2], cosLam = sqrt((1-sinLam)*(1+sinLam));
1790 //double cosPhi = T[0]/cosLam, sinPhi = T[1]/cosLam;
1791 
1792 
1793  double Bz = (fabs(Bzp) < kMinBz)? kMinBz:Bzp;
1794  double dUdH = 1;
1795  double dVdZ = cosLam;
1796  double dFita_dPhi = cosLam;
1797  double dFita_dZ = -sinLam*cosLam*cosLam*Rho;
1798  double dLama_dLam = 1;
1799  double dPinv_dRho = cosLam/Bz;
1800  double dPinv_dLam = -Rho*sinLam/Bz;
1801 
1802 
1803  double Mtx[5][5] ={
1804  // H Phi Curv Z Lam
1805  //----------------------------------------------------------
1806  /*U*/ {dUdH, 0, 0, 0, 0},
1807  /*V*/ {0, 0, 0, dVdZ, 0},
1808  /*Fita*/ {0, dFita_dPhi, 0,dFita_dZ, 0},
1809  /*Lama*/ {0, 0, 0, 0,dLama_dLam},
1810  /*Pinv*/ {0, 0,dPinv_dRho, 0,dPinv_dLam}};
1811  //-----------------------------------------------------------
1812  if ( D) memcpy(D[0],Mtx[0],sizeof(Mtx));
1813 
1814  if (Di) {
1815  double dHdU = 1;
1816  double dPhi_dFita = 1./cosLam;
1817  double dPhi_dV = sinLam*Rho;
1818  double dRho_dPinv = Bz/cosLam;
1819  double dRho_dLama = sinLam/cosLam*Rho;
1820  double dZdV = 1./cosLam;
1821  double dLam_dLama = 1;
1822  double Xtm[5][5] ={
1823 
1824  // U V Fita Lama Pinv
1825  //----------------------------------------------------------------------
1826  /*H*/ {dHdU, 0, 0, 0, 0},
1827  /*Phi*/ {0, dPhi_dV, dPhi_dFita, 0, 0},
1828  /*Curv*/ {0, 0, 0, dRho_dLama, dRho_dPinv},
1829  /*Z*/ {0, dZdV, 0, 0, 0},
1830  /*Lam*/ {0, 0, 0, dLam_dLama, 0}};
1831  //----------------------------------------------------------------------
1832 
1833  double qwe[5][5];
1834  TCL::mxmpy(Mtx[0],Xtm[0],qwe[0],5,5,5);
1835  for (int i=0;i<5;i++) {
1836  for (int j=0;j<5;j++) {
1837  double ttt = qwe[i][j]; if (i==j) ttt--;
1838  assert(fabs(ttt)<1e-8);
1839  } }
1840  memcpy(Di[0],Xtm[0],sizeof(Mtx));
1841  }
1842 
1843  if (!G) return;
1844 
1845  double hemx[15];
1846  TCL::ucopy(*he->Emx(),hemx,15);
1847  TCL::trasat(Mtx[0],hemx,G,5,5);
1848  for (int i=0,li=0;i<5;li+=++i) {
1849  assert(hemx[li+i]>0);
1850  assert( G[li+i]>0);
1851  }
1852 }
1853 //______________________________________________________________________________
1854 void THelix3d::TestConvertErrs()
1855 {
1856  double PtGev = 1.,Curv = 1./100, Mag[3] = {0,0,PtGev*Curv};
1857  double Pbeg[3],Xbeg[3]={0};
1858 
1859  TVector3 PbegV(PtGev,0,PtGev/3*4);
1860  int icharge = (-Curv*PtGev/Mag[2]<0)? -1:1;
1861 #if 0
1862  PbegV.RotateZ(gRandom->Rndm()*3.1415);
1863 #endif
1864  PbegV.GetXYZ(Pbeg);
1865  TVector3 XbegV(Xbeg);
1866 
1867  double L = 100.;
1868  TVectorD dia(5);
1869  dia[0]= 0.1; dia[1]= 3./360; dia[2]= Curv*0.01; dia[3]= 0.2; dia[4]= 2./360;
1870  for (int i=0;i<5;i++) {dia[i]*=dia[i];}
1871 
1872  TRandomVector RV(dia);
1873  auto &EMX = RV.GetMtx();
1874  auto &val = RV.GetLam();
1875  dia.Print("DIA");
1876  val.Print("VAL");
1877 
1878 
1879  double Gbeg[15],Gbeg1[15],Gend[15];
1880  for (int i=0,li=0;i< 5;li+=++i) {
1881  for (int j=0;j<=i;j++) {
1882  Gbeg[li+j] = EMX[i][j];
1883  } }
1884 
1885  THelixTrack_ TH0( Xbeg,Pbeg,Curv);
1886  THelix3d TH1(icharge ,Xbeg,Pbeg,Mag );
1887  assert (fabs(TH1.GetRho()-Curv)<1e-5);
1888  TH0.SetEmx(Gbeg);
1889 
1890  ConvertErrs(&TH0,Mag[2],Gbeg1);
1891 
1892  TVectorD Gbeg0V(15,Gbeg ); Gbeg0V.Print("Gbeg ");
1893  TVectorD Gbeg1V(15,Gbeg1 ); Gbeg1V.Print("Gbeg1");
1894 
1895  auto *emx = TH1.SetEmx(Gbeg1);
1896  emx->Print("Converted1");
1897 
1898  TH0.Move(L);
1899  TH1.Move(L);
1900  TVector3 vX0(TH0.Pos());
1901  TVector3 vX1(TH1.Pos());
1902  double eps = (vX1-vX0).Mag();
1903  assert(eps<1e-6);
1904  TVector3 vD0(TH0.Dir());
1905  TVector3 vD1(TH1.Dir());
1906  eps = (vD1-vD0).Mag();
1907  assert(eps<1e-6);
1908 
1909  ConvertErrs(&TH0,Mag[2],Gend);
1910 //const double *Gend0 = *TH0.Emx();
1911  const double *Gend1 = *TH1.Emx();
1912  TVectorD Gend0V(15,Gend ); Gend0V.Print("Gend ");
1913  TVectorD Gend1V(15,Gend1 ); Gend1V.Print("Gend1");
1914 
1915  double GendD[5],Gend1D[5];
1916  const char *tit = "UVFLP";
1917  for (int i=0,li=0;i< 5;li+=++i) {
1918  GendD [i] = sqrt(Gend [li+i]);
1919  Gend1D[i] = sqrt(Gend1[li+i]);
1920  double dif = 2*(GendD[i]-Gend1D[i])/(GendD[i]+Gend1D[i]);
1921  printf("[%c][%c] %g %g dif = %g\n",tit[i],tit[i],Gend [li+i],Gend1[li+i],dif);
1922  for (int j=0;j<i;j++) {
1923  double cor0 = Gend [li+j]/(GendD [i]*GendD [j]);
1924  double cor1 = Gend1[li+j]/(Gend1D[i]*Gend1D[j]);
1925  dif = (cor0-cor1);
1926  printf("[%c][%c] %g %g dif = %g\n",tit[i],tit[j],cor0,cor1,dif);
1927  }}
1928 
1929 }
1930 //______________________________________________________________________________
1931 void THelix3d::TestConvertDers()
1932 {
1933 // Test for derivatives + changing direction
1934 enum {kkH,kkPhi,kkCur,kkZ,kkLam};
1935 
1936  double Pt=1;
1937  double Curv = 1./100,Hz = Pt*Curv;
1938  TVector3 XbegV = TVector3(0,0,0);
1939 //?? TVector3 PbegV = TVector3(Pt,0,Pt/3*4);
1940  TVector3 PbegV = TVector3(Pt,0,0.01);
1941  TVector3 PbegV1(PbegV);
1942  TVector3 DbegV = PbegV.Unit();
1943  double HZ[3] = {0,0,Hz};
1944  int icharge = (-Curv*Pt/HZ[2]<0)? -1:1;
1945 
1946  double L = 100.,dL;
1947 
1948  double Dhlx[5][5],Dh3d[5][5],ConvEnd[5][5],ConvBeg[5][5];
1949  double Der1[5][5],NumDer1[5][5]; //dH3d/dHlx
1950  double Der2[5][5],NumDer2[5][5]; //dH3d/dHlx
1951  double NumDer3[5][5];
1952 
1953  THelixTrack_ THbeg((double*)&XbegV[0],(double*)&DbegV[0],Curv);
1954  double CosBeg = THbeg.GetCos();
1955  printf(" ==================== CosBeg = %g ==========================\n",CosBeg);
1956 
1957 //?? L = 0.5*CosBeg/Curv*M_PI; ////????
1958 
1959 
1960  THelixTrack_ THend((double*)&XbegV[0],(double*)&DbegV[0],Curv);
1961  THend.Move(L,Dhlx);
1962  double CosEnd = THend.GetCos();
1963  TVector3 XendV(THend.Pos());
1964  TVector3 DendV(THend.Dir());
1965 
1966  THelix3d T3beg(icharge,(double*)&XbegV[0],(double*)&PbegV[0],HZ);
1967  THelix3d T3end(icharge,(double*)&XbegV[0],(double*)&PbegV[0],HZ);
1968  double Pinv = T3beg.Pinv();
1969  T3end.SetDerOn();
1970  T3end.Move(L,Dh3d);
1971 
1972 // Conv = d(3d)/d(Hlx)
1973  T3beg.ConvertErrs(&THbeg,HZ[2],0,ConvBeg);
1974  T3end.ConvertErrs(&THend,HZ[2],0,ConvEnd);
1975 
1976  TCL::mxmpy(ConvEnd[0],Dhlx[0],Der1[0],5,5,5);
1977  TCL::mxmpy(Dh3d[0],ConvBeg[0],Der2[0],5,5,5);
1978  auto myDer = T3end.Der(); //derivatives
1979 
1980  TVector3 HbegV(myDer->TkDir(0)[kKU]); //H vector
1981  TVector3 UbegV(myDer->TkDir(0)[kKU]); //U vector
1982  TVector3 VbegV(myDer->TkDir(0)[kKV]); //V vector
1983  TVector3 ZbegV(0.,0.,1.); //zVECTOR
1984  TVector3 TbegV(myDer->TkDir(0)[kKT]); //Along track vector
1985 
1986  TVector3 UendV(myDer->TkDir(1)[kKU]); //U Phi vector
1987  TVector3 VendV(myDer->TkDir(1)[kKV]); //V Lam vector
1988  TVector3 TendV(myDer->TkDir(1)[kKT]); //Along track vector
1989 
1990  TVectorD der(5);
1991  for (int J=0;J<5;J++) {
1992  TVector3 X1begV(XbegV),D1begV(DbegV);
1993  double Curv1 = Curv;
1994  double delta = 1e-4;
1995  switch(J) {
1996  case kkH:{
1997  X1begV += HbegV*delta;
1998  break;
1999  }
2000  case kkPhi: { //Phi variation (beta)
2001  D1begV+=HbegV*(delta*CosBeg); D1begV.SetMag(1.);
2002  break;
2003  }
2004  case kkCur: {
2005  Curv1+=delta;
2006  break;
2007  }
2008  case kkZ: {
2009  X1begV += ZbegV*delta;
2010  break;
2011  }
2012  case kkLam: { //lamda variation (alfa)
2013  D1begV+=VbegV*delta; D1begV.SetMag(1.);
2014  break;
2015  }
2016 
2017  default: assert(0 && "Wrong J");
2018  }
2019 // Now separate THelx3d derivative
2020  TVector3 X2begV(XbegV),D2begV(DbegV),P2begV;
2021  double P2invBeg=Pinv;
2022  double delta2 = delta;
2023  switch(J) {
2024  case kU:{
2025  X2begV += UbegV*delta2;
2026  break;
2027  }
2028  case kFita: { //Fita variation (beta)
2029  D2begV+=UbegV*(delta2); D2begV.SetMag(1.);
2030 
2031  break;
2032  }
2033  case kPinv: {
2034  P2invBeg+= delta2;
2035  break;
2036  }
2037  case kV: {
2038  X2begV += VbegV*delta2;
2039  break;
2040  }
2041  case kLama: { //lama variation
2042  delta2 = 3.14/180;
2043  D2begV+=VbegV*delta2; D1begV.SetMag(1.);
2044  break;
2045  }
2046 
2047  default: assert(0 && "Wrong J");
2048  }
2049 
2050 
2051 
2052 // By THelixTrack_
2053  THelixTrack_ TH1((double*)&X1begV[0],(double*)&D1begV[0],Curv1);
2054  double CosBeg1 = TH1.GetCos();
2055  TH1.Move(L);
2056  TVector3 X1endV(TH1.Pos());
2057  TVector3 D1endV(TH1.Dir());
2058  dL = -(X1endV-XendV).Dot(DendV)/(DendV.Dot(D1endV));
2059  TH1.Move(dL);
2060  double CosEnd1 = TH1.GetCos();
2061  X1endV = TVector3(TH1.Pos());
2062  D1endV = TVector3(TH1.Dir());
2063 
2064 
2065  TVectorD dif(5);
2066  TVector3 difX = X1endV-XendV;
2067  TVector3 difD = (D1endV-DendV);
2068  dif[kU] = difX.Dot(UendV);
2069  dif[kV] = difX.Dot(VendV);
2070  dif[kPinv] = -icharge*(Curv1*CosEnd1-Curv*CosEnd)/Hz;
2071  dif[kFita] = difD.Dot(UendV);
2072  dif[kLama] = difD.Dot(VendV);
2073  assert(fabs(dif.NormInf())<1.);
2074  dif*=1./delta;
2075  for (int i=0;i<5;i++) { NumDer1[i][J] = dif[i];}
2076 
2077 // By THelix3d
2078  double P1tot1 = Hz/(Curv1*CosBeg1);
2079  auto P1begV = D1begV*=P1tot1;
2080  THelix3d T31(icharge,(double*)&X1begV[0],(double*)&P1begV[0],HZ);
2081  T31.Move(L);
2082  TVector3 P1endV;
2083  T31.Eval(0,&X1endV[0],&P1endV[0]);
2084  D1endV = P1endV.Unit();
2085  dL = -(X1endV-XendV).Dot(DendV)/(DendV.Dot(D1endV));
2086  TH1.Move(dL);
2087  X1endV = TVector3(T31.Pos());
2088  D1endV = TVector3(T31.Dir());
2089  double PinvEnd = T31.Pinv();
2090  difX = X1endV-XendV;
2091  difD = (D1endV-DendV);
2092  dif[kU] = difX.Dot(UendV);
2093  dif[kV] = difX.Dot(VendV);
2094  dif[kPinv] = PinvEnd-Pinv;
2095  dif[kFita] = difD.Dot(UendV);
2096  dif[kLama] = difD.Dot(VendV);
2097  assert(fabs(dif.NormInf())<1.);
2098  dif*=1./delta;
2099  for (int i=0;i<5;i++) { NumDer2[i][J] = dif[i];}
2100 
2101 
2102 // By THelixTrack_ but 3d derivative
2103  double CosBeg2 = sqrt((1.-D2begV.CosTheta())*(1.+D2begV.CosTheta()));
2104  double Curv2beg = -icharge*P2invBeg/CosBeg2*Hz;
2105 #if 1
2106  THelixTrack_ TH2((double*)&X2begV[0],(double*)&D2begV[0],Curv2beg);
2107 #endif
2108 #if 0
2109  P2begV = D2begV*(1/fabs(P2invBeg));
2110  THelix3d TH2(icharge,(double*)&X2begV[0],(double*)&P2begV[0],HZ);
2111 #endif
2112 
2113  TH2.Move(L);
2114 
2115 
2116 
2117  double CosEnd2 = TH2.GetCos();
2118  double CurvEnd2 = TH2.GetRho();
2119  TVector3 X2endV(TH2.Pos());
2120  TVector3 D2endV(TH2.Dir());
2121  difX = X2endV-XendV;
2122  difD = (D2endV-DendV);
2123  dif[kU] = difX.Dot(UendV);
2124  dif[kV] = difX.Dot(VendV);
2125  dif[kPinv] = -icharge*(CurvEnd2*CosEnd2-Curv*CosEnd)/Hz;
2126  dif[kFita] = difD.Dot(UendV);
2127  dif[kLama] = difD.Dot(VendV);
2128 //??? assert(fabs(dif.NormInf())<1.);
2129  dif*=1./delta2;
2130  for (int i=0;i<5;i++) { NumDer3[i][J] = dif[i];}
2131 
2132 
2133 
2134  }
2135 
2136 static const char* LabHlx[]={"H ","Phi","Cur","Z ","Lam"};
2137 static const char* LabH3d[]={"U ","V ","Fit","Ama","Pin"};
2138  for (int ihlx=0;ihlx<5;ihlx++) {
2139  for (int ih3d=0;ih3d<5;ih3d++) {
2140  printf("[%s][%s] = %g == %g == %g == %g\n",LabH3d[ih3d],LabHlx[ihlx]
2141  ,Der1[ih3d][ihlx]
2142  ,Der2[ih3d][ihlx]
2143  ,NumDer1[ih3d][ihlx]
2144  ,NumDer2[ih3d][ihlx]);
2145  }}
2146 printf("\n\n");
2147  for (int ih3d=0;ih3d<5;ih3d++) {
2148  for (int jh3d=0;jh3d<5;jh3d++) {
2149  printf("[%s][%s] = %g == %g \n",LabH3d[ih3d],LabH3d[jh3d]
2150  ,Dh3d[ih3d][jh3d]
2151  ,NumDer3[ih3d][jh3d]);
2152  }}
2153 
2154 }
2155 //____________________________________________________________
2156 double EmxSign(int n,const double *a)
2157 {
2158  double ans=3e33;
2159  double buf[55];
2160  double *B = (n<=10) ? buf : new double[n];
2161  double *b = B;
2162  // trchlu.F -- translated by f2c (version 19970219).
2163  //
2164  //see original documentation of CERNLIB package F112
2165 
2166  /* Local variables */
2167  int ipiv, kpiv, i__, j;
2168  double r__, dc;
2169  int id, kd;
2170  double sum;
2171 
2172 
2173  /* CERN PROGLIB# F112 TRCHLU .VERSION KERNFOR 4.16 870601 */
2174  /* ORIG. 18/12/74 W.HART */
2175 
2176 
2177  /* Parameter adjuTments */
2178  --b; --a;
2179 
2180  /* Function Body */
2181  ipiv = 0;
2182 
2183  i__ = 0;
2184 
2185  do {
2186  ++i__;
2187  ipiv += i__;
2188  kpiv = ipiv;
2189  r__ = a[ipiv];
2190 
2191  for (j = i__; j <= n; ++j) {
2192  sum = 0.;
2193  if (i__ == 1) goto L40;
2194  if (r__ == 0.) goto L42;
2195  id = ipiv - i__ + 1;
2196  kd = kpiv - i__ + 1;
2197 
2198  do {
2199  sum += b[kd] * b[id];
2200  ++kd; ++id;
2201  } while (id < ipiv);
2202 
2203 L40:
2204  sum = a[kpiv] - sum;
2205 L42:
2206  if (j != i__) b[kpiv] = sum * r__;
2207  else {
2208  if (sum<ans) ans = sum;
2209  if (sum<0.) goto RETN;
2210  dc = sqrt(sum);
2211  b[kpiv] = dc;
2212  if (r__ > 0.) r__ = (double)1. / dc;
2213  }
2214  kpiv += j;
2215  }
2216 
2217  } while (i__ < n);
2218 
2219 RETN: if (B!=buf) delete B;
2220  return ans;
2221 } /* trchlu_ */
THelixTrack_::Backward
void Backward()
Change direction.
Definition: THelixTrack_.cxx:338
THelix3d::Eval
double Eval(double step, double xyz[3]=0, double mom[3]=0)
Evaluate params with given step along helix.
Definition: THelix3d.cxx:304
TkDir_t
Definition: THelix3d.h:15
THEmx3d_t
Definition: THelix3d.h:62
THelixTrack_::Path
double Path(double stmax, const double *surf, int nsurf, double *x=0, double *dir=0, int nearest=0) const
Definition: THelixTrack_.cxx:443
TCL::trsinv
static float * trsinv(const float *g, float *gi, int n)
Definition: TCernLib.cxx:1160
THelixTrack_::Move
double Move(double step)
Move along helix.
Definition: THelixTrack_.cxx:421
THelix3d::Pos
const double * Pos() const
distance and DCAxy and DCAz to given space point (with error matrix)
Definition: THelix3d.h:161
THelix3d::Backward
void Backward()
Change direction.
Definition: THelix3d.cxx:269
THelix3d::Move
double Move(double step, double F[5][5]=0)
Move along helix.
Definition: THelix3d.cxx:281
TCL::trasat
static float * trasat(const float *a, const float *s, float *r, int m, int n)
Definition: TCernLib.cxx:540
THelixTrack_
Definition: THelixTrack_.h:214
THelixTrack_::MakeMtx
void MakeMtx(double step, double F[5][5])
Definition: THelixTrack_.cxx:388
THelixTrack_::Eval
double Eval(double step, double *xyz, double *dir=0, double *rho=0) const
Evaluate params with given step along helix.
Definition: THelixTrack_.cxx:1030
THelix3d::Path
double Path(const double point[3], double xyz[3]=0, double mom[3]=0)
Distance to nearest point to given space point.
Definition: THelix3d.cxx:321
TDataSet::Path
virtual TString Path() const
return the full path of this data set
Definition: TDataSet.cxx:626
THelix3d::Dca
double Dca(const double point[3], double *dcaErr=0)
DCA to given space point (with error matrix)
Definition: THelix3d.cxx:339
THDer3d_t
Definition: THelix3d.h:40
Generated by   doxygen 1.8.5