StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
TRungeKutta.cxx
1 #include <stdlib.h>
2 #include <math.h>
3 #include "TError.h"
4 #include "TCanvas.h"
5 #include "TGraph.h"
6 #include "TSystem.h"
7 #include "TMath.h"
8 #include "TCernLib.h"
9 #include "TVector3.h"
10 #include "TVectorD.h"
11 #include "TRandom.h"
12 #include "TRandomVector.h"
13 #include "TH1.h"
14 #include "TCernLib.h"
15 #include "TRungeKutta.h"
16 
17 TRKuttaMag *TRungeKutta::fgMag=0;
18 
19 
20 //_____________________________________________________________________________
21 //_____________________________________________________________________________
22 ClassImp(TRungeKutta)
23 //_____________________________________________________________________________
24 void TRungeKutta::TRKutaPars_t::Print(const char *tit) const
25 {
26  if (tit && *tit) {printf("TRKutaPars_t::Print(%s)\n",tit);}
27  printf("X= %g %g %g D= %g %g %g Pinv = %g\n"
28  ,X[0],X[1],X[2],D[0],D[1],D[2],Pinv);
29 }
30 //_____________________________________________________________________________
31 void TRungeKutta::Print(const char *tit) const
32 {
33  if (tit && *tit) {printf("TRungeKutta::Print(%s)\n",tit);}
34  fInp.Print("");
35 }
36 
37 //_____________________________________________________________________________
38 TRungeKutta::TRungeKutta(int charge,const double *xyz,const double *mom,TRKuttaMag *Mag)
39 {
40 // Generalization of RungeKutta algorithm with derivative and error propagationf
41 // V.Perevoztchikov
42 //
43  memset(fBeg,0,fEnd-fBeg);
44  SetMag(Mag);
45  Set(charge,xyz,mom);
46 }
47 //_____________________________________________________________________________
48 TRungeKutta::TRungeKutta()
49 {
50  memset(fBeg,0,fEnd-fBeg);
51  SetMag(0);
52 }
53 //_____________________________________________________________________________
54 TRungeKutta::~TRungeKutta()
55 {
56  delete fEmx3d;
57  delete fDer3d;
58 }
59 //_____________________________________________________________________________
60 TRungeKutta::TRungeKutta(const THelixTrack_ &from,TRKuttaMag* mag)
61 {
62  memset(fBeg,0,fEnd-fBeg);
63  SetMag(mag);
64  double B[3],P[3];
65  (*fGUFLD)(from.Pos(),B);
66  assert(fabs(B[0])+fabs(B[1])<1e-3*fabs(B[2]));
67  double rho = from.GetRho();
68  int iq = ((rho>0)==(B[2]>0)) ? -1:1;
69  TCL::vscale(from.Dir(),B[2]/rho,P,3);
70  Set(iq,from.Pos(),P);
71 }
72 //_____________________________________________________________________________
73 void TRungeKutta::SetMag(TRKuttaMag* mag)
74 {
75  fGUFLD = (mag)? mag:fgMag;
76  fBInp[2] = 1e11;
77  assert(fGUFLD);
78 }
79 //_____________________________________________________________________________
80 void TRungeKutta::ZetMag(const double *mag)
81 {
82  memcpy(fBInp,mag,sizeof(fBInp));
83  memcpy(fBOut,mag,sizeof(fBOut));
84 }
85 //_____________________________________________________________________________
86 void TRungeKutta::Set(int charge,const double *xyz,const double *mom)
87 {
88  if (charge) fCharge = charge;
89  memcpy(fInp.X,xyz,sizeof(fInp.X));
90  memcpy(fInp.P,mom,sizeof(fInp.P));
91  fInp.Pinv = -fCharge/sqrt(TCL::vdot(fInp.P,fInp.P,3));
92  TCL::vscale(fInp.P,fabs(fInp.Pinv),fInp.D,3);
93  if (fBInp[2]<1e3) return;
94  (*fGUFLD)(fInp.X,fBInp);
95 
96 }
97 //_____________________________________________________________________________
98 TRungeKutta::TRungeKutta(const TRungeKutta *from)
99 {
100  *this = *from; fEmx3d=0,fDer3d=0;
101 }
102 //_____________________________________________________________________________
103 TRungeKutta::TRungeKutta(const TRungeKutta &from)
104 {
105  *this = from;
106  if (fEmx3d) {fEmx3d = new THEmx3d_t; *fEmx3d = *from.fEmx3d;}
107  if (fDer3d) {fDer3d = new THDer3d_t; *fDer3d = *from.fDer3d;}
108 }
109 //_____________________________________________________________________________
110 void TRungeKutta::Backward()
111 {
112  TCL::vscale(fInp.P,-1.,fInp.P,3);
113  TCL::vscale(fInp.D,-1.,fInp.D,3);
114  fCharge = - fCharge;
115  fInp.Pinv = -fInp.Pinv;
116  if (fEmx3d) fEmx3d->Backward();
117  if (fDer3d) fDer3d->Backward();
118 
119 }
120 //_____________________________________________________________________________
121 double TRungeKutta::Move(double step)
122 {
123  fLen = step;
124  if (fabs(fLen)<=kMicron) return 0;
125  grkuta(fCharge ,fLen,fInp,fOut) ;
126  fOut.Update();
127  Move();
128  return fLen;
129 }
130 //______________________________________________________________________________
131 double TRungeKutta::GetCurv() const //Full curvature in current point
132 {
133  double H = sqrt(fBInp[0]*fBInp[0]+fBInp[1]*fBInp[1]+fBInp[2]*fBInp[2]);
134  return fInp.Pinv*H;
135 }
136 //______________________________________________________________________________
137 int TRungeKutta::GetDir() const //
138 {
139  double qwe = TCL::vdot(Pos(),Dir(),3);
140  return (qwe<0)? 0:1;
141 }
142 //______________________________________________________________________________
143 double TRungeKutta::Eval(double step, double xyz[3], double mom[3]) const
144 {
145  fLen = step;
146  if (fabs(step)<=0) {
147  if (xyz) memcpy(xyz,fInp.X,sizeof(fInp.X));
148  if (mom) memcpy(mom,fInp.P,sizeof(fInp.P));
149  return fLen;
150  }
151  grkuta(fCharge ,fLen,fInp,fOut) ;
152  fOut.Update();
153  if (xyz) memcpy(xyz,fOut.X,sizeof(fOut.X));
154  if (mom) memcpy(mom,fOut.P,sizeof(fOut.P));
155  return fLen;
156 }
157 //______________________________________________________________________________
158 double TRungeKutta::Move()
159 {
160  if (IsDerOn()) {
161  MakeMtx();
162  }
163  fInp = fOut;
164  return fLen;
165 }
166 //______________________________________________________________________________
167 double TRungeKutta::Path(const double point[3], int idir,int ndca) const
168 {
169 static const double kEps=1e-7,kRef=1e-6;
170 assert(ndca==2 || ndca==3);
171  double dx[3],gate[4]={-1e11,1e11};
172  int myDir = idir, fail = 13,myState=0;
173  fLen = 0; fOut = fInp;
174  for (int iter=0; iter<=100; iter++) {
175  double radius = 1./(fabs(GetCurv())+1e-11);
176  double maxStep = (radius<100)? radius:100;
177 
178  auto inp = fOut;
179 assert(fLen>=gate[0] && fLen<=gate[1]);
180  TCL::vsub(point,inp.X,dx,ndca);
181  double dis = sqrt(TCL::vdot(dx,dx,ndca));
182  if (maxStep>dis) maxStep=dis;
183  double tau = TCL::vdot(dx,inp.D,ndca);
184  if (ndca==2) tau/=sqrt((1.-inp.D[2])*(1.+inp.D[2]));
185  double space = fabs(inp.X[0])+fabs(inp.X[1]);
186  if (fabs(tau)<kEps || fabs(tau)<kRef*space) {
187  fLen+=tau; fail=0; break;
188  }
189 
190 
191  if (myDir>0) {
192  if (tau<0) { tau= maxStep;} else { myDir = 0;}
193  } else if (myDir<0){
194  if (tau>0) { tau=-maxStep;} else { myDir = 0;}
195  }
196  if (fabs(tau)>maxStep) tau = (tau<0)? -maxStep:maxStep;
197  if (tau<0) { myState|=2;
198  if (gate[1]>fLen) {gate[1]=fLen;gate[3]=-tau;}
199  }
200  else { myState|=1;
201  if (gate[0]<fLen) {gate[0]=fLen;gate[2] = tau;}
202  }
203  if (myState==3) {
204  tau = ((gate[0]-fLen)*gate[3]+(gate[1]-fLen)*gate[2])/(gate[3]+gate[2]);
205  int myLen = fLen+tau;
206  if (iter>20 || myLen<=gate[0] || myLen>=gate[1]) tau = (gate[0]-fLen+gate[1]-fLen)/2;
207  } else {
208  tau = (tau<0)? -maxStep:maxStep;
209  }
210  grkuta(fCharge ,tau,inp,fOut) ;
211  fLen += tau;
212 assert(fLen>=gate[0] && fLen<=gate[1]);
213  }
214  if (fail) return 2e22;
215 
216  fOut.Update();
217  return fLen;
218 }
219 //______________________________________________________________________________
220 double TRungeKutta::Path(double x,double y, int idir) const
221 {
222  double gde[2];
223  gde[0] = x;gde[1] = y;
224  return Path(gde,idir,2);
225 }
226 #if 0
227 //______________________________________________________________________________
228 double TRungeKutta::Path(const double point[3], int idir) const
229 {
230 static const double kEps=1e-5,kRef=1e-5,kBigEps=1e-1;
231  double qaLeft,posLeft=0,qaRite=0,posRite=0,qaMine,qaStep,dx[3];
232  double radius = 1./(fabs(GetCurv())+1e-11);
233  radius/= sqrt(1. - fOut.D[2]*fOut.D[2])+1e-11;
234  double maxStep = (radius<100)? radius:100;
235  double minStep = 0.1;
236 
237  fLen = 0; fOut = fInp;
238  for (int iter=0; iter<=20; iter++) {
239  auto inp = fOut;
240  TCL::vsub(inp.X,point,dx,3);
241  qaMine = -TCL::vdot(dx,inp.D,3);
242  if (!iter) {qaLeft = qaMine;}
243  else if (qaMine*qaLeft>0) {qaLeft = qaMine;}
244  else {qaRite = qaMine; posLeft=-qaStep; break ;}
245  qaStep = fabs(qaMine); if (iter) qaStep*=2;
246  if (qaStep > maxStep) qaStep = maxStep;
247  if (qaStep < minStep) qaStep = minStep;
248  if (idir) { qaStep*=idir;} else if (qaMine<0) {qaStep*=-1;};
249 
250  grkuta(fCharge ,qaStep,inp,fOut) ;
251  fLen += qaStep;
252  if (fabs(fLen)>6*radius) return 1e11;
253  assert(fabs(fLen)<6*radius);
254  }
255  if (!qaRite) return 2e22;
256  assert(qaRite);
257  for (int iter=0; iter<=20; iter++) {
258  auto inp = fOut;
259  double denom = qaRite-qaLeft;
260  if ((iter&3)==3 || fabs(denom) < 1e-4) { //Half step
261  qaStep = 0.5*(posRite+posLeft);
262  } else { //Linear approach
263  qaStep = -(qaLeft*posRite-qaRite*posLeft)/denom;
264  }
265  grkuta(fCharge ,qaStep,inp,fOut) ;
266  fLen += qaStep;
267  TCL::vsub(fOut.X,point,dx,3);
268  qaMine = -TCL::vdot(dx,fOut.D,3);
269 
270  if (qaLeft*qaMine>0) { qaLeft = qaMine; posLeft=0; posRite-=qaStep;}
271  else { qaRite = qaMine; posRite=0; posLeft-=qaStep;}
272  if (fabs(qaMine)>kBigEps) continue;
273  double dis2 = TCL::vdot(dx,dx,3);;
274  if (fabs(qaMine) < kEps || qaMine*qaMine < kRef*kRef*dis2) break;
275  }
276 
277  fOut.Update();
278  return (fabs(qaMine)<kBigEps)? fLen:1e11;
279 }
280 #endif
281 //______________________________________________________________________________
282 double TRungeKutta::Dca(const double point[3],double *dcaErr) const
283 {
284 if(point || dcaErr){};
285 assert(!"TRungeKutta::Dca Not implemented");
286 return 0;
287 }
288 #if 0
289 //______________________________________________________________________________
290 double TRungeKutta::Path(double x,double y, int idir) const
291 {
292 static const double kEps=1e-5,kRef=1e-5,kBigEps=1e-1;
293  double qaLeft,posLeft=0,qaRite=0,posRite=0,qaMine,qaStep,dx[2];
294  double minStep = 0.1;
295  double radius = 1./(fabs(GetCurv())+1e-11);
296  radius/= sqrt(1. - fOut.D[2]*fOut.D[2])+1e-11;
297  double maxStep = (radius<100)? radius:100;
298 
299  fLen = 0; fOut = fInp;
300  for (int iter=0; iter<=20; iter++) {
301  auto inp = fOut;
302  dx[0] = inp.X[0]-x;dx[1] = inp.X[1]-y;
303  qaMine = -(dx[0]*inp.D[0] + dx[1]*inp.D[1]);
304  qaMine /=(1.-inp.D[2])*(1.+inp.D[2]);
305  if (!iter) {qaLeft = qaMine;}
306  else if (qaMine*qaLeft>0) {qaLeft = qaMine;}
307  else {qaRite = qaMine; posLeft=-qaStep; break ;}
308  qaStep = fabs(qaMine); if (iter) qaStep*=2;
309  if (qaStep > maxStep) qaStep = maxStep;
310  if (qaStep < minStep) qaStep = minStep;
311  if (idir) { qaStep*=idir;} else if (qaMine<0) {qaStep*=-1;};
312 
313 
314  grkuta(fCharge ,qaStep,inp,fOut) ;
315  fLen += qaStep;
316  if (fabs(fLen)>6*radius) return 1e11;
317  assert(fabs(fLen)<6*radius);
318  if (fabs(fLen)>6*radius) return 2e22;
319  }
320  if(!qaRite) return 3e33;;
321  assert(qaRite);
322  for (int iter=0; iter<=20; iter++) {
323  auto inp = fOut;
324  double denom = qaRite-qaLeft;
325  if ((iter&3)==3 || fabs(denom) < 1e-4) { //Half step
326  qaStep = 0.5*(posRite+posLeft);
327  } else { //Linear approach
328  qaStep = -(qaLeft*posRite-qaRite*posLeft)/denom;
329  }
330  grkuta(fCharge ,qaStep,inp,fOut) ;
331  fLen += qaStep;
332  dx[0] = fOut.X[0]-x;dx[1] = fOut.X[1]-y;
333  qaMine = -(dx[0]*fOut.D[0]+dx[1]*fOut.D[1]);
334  qaMine /= (1.-fOut.D[2])*(1.+fOut.D[2]);
335 
336  if (qaLeft*qaMine>0) { qaLeft = qaMine; posLeft=0; posRite-=qaStep;}
337  else { qaRite = qaMine; posRite=0; posLeft-=qaStep;}
338  if (fabs(qaMine)>kBigEps) continue;
339  double dis2 = dx[0]*dx[0]+dx[1]*dx[1];
340  if (fabs(qaMine) < kEps || qaMine*qaMine< kRef*kRef*dis2) break;
341  }
342 
343  fOut.Update();
344  return (fabs(qaMine)<kBigEps)? fLen:1e11;
345 }
346 #endif
347 //______________________________________________________________________________
348 double TRungeKutta::Dca(double x,double y, double *dcaErr) const
349 {
350 if (x || y){}
351 return 0;
352 assert(!"TRungeKutta::Dca(x,y,dcaError) n oy implemented");
353 }
354 //______________________________________________________________________________
355 void TRungeKutta::MakeMtx()
356 {
357  if (fabs(fLen) < kMicron) return;
358  if (!fDer3d) fDer3d = new THDer3d_t;
359 
360 // Move from Beg to End
361  THelix3d TH0(fCharge,fInp.X,fInp.P,fBInp);
362  TH0.SetDerOn();
363  TH0.SetEmx(fEmx3d);
364  TH0.Move(fLen);
365  auto *myDer = TH0.Der();
366  *fDer3d = *myDer;
367  THEmx3d_t Emx1 = *(TH0.Emx());
368 
369 // Move from End to Beg and then back to End
370  TH0.Set(fCharge,fOut.X,fOut.P,fBOut);
371  TH0.SetDerOn(0);
372  TH0.Move(-fLen);
373  TH0.SetEmx(fEmx3d);
374  TH0.SetDerOn();
375  TH0.Move(fLen);
376  THEmx3d_t &Emx2 = *(TH0.Emx());
377  Emx1.Update(Emx2.TkDir());
378 
379  TCL::vlinco(Emx1,0.5,Emx2,0.5,(*fEmx3d),15);
380  fEmx3d->TkDir() = Emx1.TkDir();
381 
382 }
383 //______________________________________________________________________________
384 void TRungeKutta::SetEmx(const double err[15])
385 {
386  SetDerOn();
387  if (!fEmx3d) {fEmx3d = new THEmx3d_t;}
388  THelix3d::MakeTkDir(fInp.D,fBInp,fEmx3d->TkDir());
389  if (err) fEmx3d->Set(err);
390 }
391 //______________________________________________________________________________
392 void TRungeKutta::SetEmx(const THEmx3d_t *err)
393 {
394  SetDerOn();
395  if (!fEmx3d) {fEmx3d = new THEmx3d_t;}
396  if (err) *fEmx3d = *err;
397  THelix3d::MakeTkDir(fInp.D,fBInp,fEmx3d->TkDir());
398 }
399 //______________________________________________________________________________
400 TRungeKutta &TRungeKutta::operator=(const TRungeKutta &from)
401 {
402  memcpy(fBeg,from.fBeg,fEnd-fBeg+1);
403  if (fEmx3d) {
404  fEmx3d = new THEmx3d_t;
405  *fEmx3d = *from.fEmx3d;
406  }
407  if (fDer3d) {
408  fDer3d = new THDer3d_t;
409  *fDer3d = *from.fDer3d;
410  }
411  return *this;
412 }
413 //______________________________________________________________________________
414 void TRungeKutta::grkuta(double CHARGE,double STEP
415  ,const TRKutaPars_t &VECT,TRKutaPars_t &VOUT) const
416 {
417 /*
418 * Revision 1.1.1.1 1997/11/03 15:30:47 atlascvs
419 * Importing CERNLIB version 08.21.
420 *
421 * Revision 1.1.1.1 1995/10/24 10:21:42 cernlib
422 * Geant
423 *
424 *
425 #include "geant321/pilot.h"
426 *CMZ : 3.21/02 29/03/94 15.41.23 by S.Giani
427 *-- Author :
428  SUBROUTINE GRKUTA (CHARGE,STEP,VECT,VOUT)
429 C.
430 C. ******************************************************************
431 C. * *
432 C. * Runge-Kutta method for tracking a particle through a magnetic *
433 C. * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of *
434 C. * Standards, procedure 25.5.20) *
435 C. * *
436 C. * Input parameters *
437 C. * CHARGE Particle charge *
438 C. * STEP Step size *
439 C. * VECT Initial co-ords,direction cosines,momentum *
440 C. * Output parameters *
441 C. * VOUT Output co-ords,direction cosines,momentum *
442 C. * User routine called *
443 C. * CALL GUFLD(X,F) *
444 C. * *
445 C. * ==>Called by : <USER>, GUSWIM *
446 C. * Authors R.Brun, M.Hansroul ********* *
447 C. * V.Perevoztchikov (CUT STEP implementation) *
448 C. * *
449 C. * *
450 C. ******************************************************************
451 C.
452 */
453  double F[4];
454  double XYZT[3], XYZ[3];
455  double SECXS[4],SECYS[4],SECZS[4],HXP[3];
456 // EQUIVALENCE (X,XYZ(1)),(Y,XYZ(2)),(Z,XYZ(3)),
457 // + (XT,XYZT(1)),(YT,XYZT(2)),(ZT,XYZT(3))
458 //*
459  double& X = XYZ[0]; double& Y = XYZ[1]; double& Z = XYZ[2];
460  double& XT = XYZT[0]; double& YT = XYZT[1]; double& ZT = XYZT[2];
461 
462  enum {MAXIT = 1992, MAXCUT = 11};
463  static const double EC=2.9979251e-4,DLT=1e-6,DLT32=DLT/32;;
464  static const double ZERO=0, ONE=1, TWO=2, THREE=3;
465  static const double THIRD=ONE/THREE, HALF=ONE/TWO;
466  static const double PISQUA=.986960440109e+01;
467 // PARAMETER (IX=1,IY=2,IZ=3,IPX=4,IPY=5,IPZ=6)
468 /*.
469 *. ------------------------------------------------------------------
470 *.
471 * This constant is for units CM,GEV/C and KGAUSS
472 */
473  double A,B,C,AT,BT,CT;
474  double dA,dB,dC,EZT;
475  double PINV,TL,H,H2,H4,REST,F1,F2,F3,F4,PH,PH2,DXT,DYT,DZT ;
476  double EST,CBA,RHO,TET,HNORM,RHO1,SINT,COST;
477  double G1,G2,G3,G4,G5,G6,HP,ANG2;
478  if (EC||ZERO||EZT){}
479  fNStps = 0;
480  int ITER = 0;
481  int NCUT = 0;
482  VOUT = VECT;
483  //VP PINV = EC * CHARGE / VECT.Mom;
484  PINV = -VECT.Pinv; //the opposite sign to STAR definition
485  TL = 0.;
486  H = STEP;
487 //*
488 //*
489  L20: REST = STEP-TL;
490  if (fabs(H) > fabs(REST)) H = REST;
491  (*fGUFLD)(VOUT.X,F);
492 
493 //*
494 //* Start of integration
495 //*
496  X = VOUT.X[0];
497  Y = VOUT.X[1];
498  Z = VOUT.X[2];
499  A = VOUT.D[0];
500  B = VOUT.D[1];
501  C = VOUT.D[2];
502 //*
503  H2 = HALF * H;
504  H4 = HALF * H2;
505  PH = PINV * H;
506  PH2 = HALF * PH;
507  SECXS[0] = (B * F[2] - C * F[1]) * PH2;
508  SECYS[0] = (C * F[0] - A * F[2]) * PH2;
509  SECZS[0] = (A * F[1] - B * F[0]) * PH2;
510  ANG2 = (SECXS[0]*SECXS[0] + SECYS[0]*SECYS[0] + SECZS[0]*SECZS[0]);
511  if (ANG2 > PISQUA) goto L40;
512  DXT = H2 * A + H4 * SECXS[0];
513  DYT = H2 * B + H4 * SECYS[0];
514  DZT = H2 * C + H4 * SECZS[0];
515  XT = X + DXT;
516  YT = Y + DYT;
517  ZT = Z + DZT;
518 //*
519 //* Second intermediate point
520 //*
521  EST = fabs(DXT)+fabs(DYT)+fabs(DZT);
522  if (EST > fabs(H)) goto L30;
523 
524  (*fGUFLD)(XYZT,F);
525  AT = A + SECXS[0];
526  BT = B + SECYS[0];
527  CT = C + SECZS[0];
528 //*
529  SECXS[1] = (BT * F[2] - CT * F[1]) * PH2;
530  SECYS[1] = (CT * F[0] - AT * F[2]) * PH2;
531  SECZS[1] = (AT * F[1] - BT * F[0]) * PH2;
532  AT = A + SECXS[1];
533  BT = B + SECYS[1];
534  CT = C + SECZS[1];
535  SECXS[2] = (BT * F[2] - CT * F[1]) * PH2;
536  SECYS[2] = (CT * F[0] - AT * F[2]) * PH2;
537  SECZS[2] = (AT * F[1] - BT * F[0]) * PH2;
538  DXT = H * (A + SECXS[2]);
539  DYT = H * (B + SECYS[2]);
540  DZT = H * (C + SECZS[2]);
541  XT = X + DXT;
542  YT = Y + DYT;
543  ZT = Z + DZT;
544  AT = A + TWO*SECXS[2];
545  BT = B + TWO*SECYS[2];
546  CT = C + TWO*SECZS[2];
547 //*
548  EST = fabs(DXT)+fabs(DYT)+fabs(DZT);
549  if (EST > 2.*fabs(H)) goto L30;
550 
551  (*fGUFLD)(XYZT,F);
552  memcpy(fBOut,F,sizeof(fBOut));
553 //*
554  Z = Z + (C + (SECZS[0] + SECZS[1] + SECZS[2]) * THIRD) * H;
555  Y = Y + (B + (SECYS[0] + SECYS[1] + SECYS[2]) * THIRD) * H;
556  X = X + (A + (SECXS[0] + SECXS[1] + SECXS[2]) * THIRD) * H;
557 //*
558  SECXS[3] = (BT*F[2] - CT*F[1])* PH2;
559  SECYS[3] = (CT*F[0] - AT*F[2])* PH2;
560  SECZS[3] = (AT*F[1] - BT*F[0])* PH2;
561  dA = (SECXS[0]+SECXS[3]+TWO * (SECXS[1]+SECXS[2])) * THIRD;
562  dB = (SECYS[0]+SECYS[3]+TWO * (SECYS[1]+SECYS[2])) * THIRD;
563  dC = (SECZS[0]+SECZS[3]+TWO * (SECZS[1]+SECZS[2])) * THIRD;
564  A+=dA; B+=dB; C+=dC; EZT = fabs(dA)+fabs(dB) + fabs(dC);
565 //*
566  EST = fabs(SECXS[0]+SECXS[3] - (SECXS[1]+SECXS[2]))
567  + fabs(SECYS[0]+SECYS[3] - (SECYS[1]+SECYS[2]))
568  + fabs(SECZS[0]+SECZS[3] - (SECZS[1]+SECZS[2]));
569 //*
570  if (/*EST > 1e-2*EZT &&*/ EST > DLT && fabs(H) > 1.E-2) goto L30;
571  ITER = ITER + 1;
572  NCUT = 0;
573 //* If too many iterations, go to HELIX
574  if (ITER > MAXIT) goto L40;
575 //*;
576  TL = TL + H;
577  if (EST < (DLT32)) {//then;
578  H = H*TWO;
579  }//endif
580  CBA = ONE/ sqrt(A*A + B*B + C*C);
581 // iqwe++;
582 // printf("%d UUUUUUUUUUU dX = %g %g %g\n",iqwe, X-VOUT.X[0],Y-VOUT.X[1],Z-VOUT.X[2]);
583 
584 
585  VOUT.X[0] = X;
586  VOUT.X[1] = Y;
587  VOUT.X[2] = Z;
588  VOUT.D[0] = CBA*A;
589  VOUT.D[1] = CBA*B;
590  VOUT.D[2] = CBA*C;
591  REST = STEP - TL;
592  if (STEP < 0.) REST = -REST;
593 
594  if (REST > 1.E-5*fabs(STEP)) goto L20;
595  if (REST > 1.E-5) goto L20; //VP
596 //*
597  goto L999;
598 //*
599 //** CUT STEP
600  L30: NCUT = NCUT + 1;
601 //* If too many cuts , go to HELIX
602  if (NCUT > MAXCUT) goto L40;
603  H = H*HALF;
604  goto L20;
605 //*
606 //** ANGLE TOO BIG, USE HELIX
607  L40: F1 = F[0];
608  F2 = F[1];
609  F3 = F[2];
610  F4 = sqrt(F1*F1+F2*F2+F3*F3);
611  RHO = -F4*PINV;
612  TET = RHO * STEP;
613  if(fabs(TET) > 1e-6) {//then;
614  HNORM = ONE/F4;
615  F1 = F1*HNORM;
616  F2 = F2*HNORM;
617  F3 = F3*HNORM;
618 //*
619  HXP[0] = F2*VECT.D[2] - F3*VECT.D[1];
620  HXP[1] = F3*VECT.D[0] - F1*VECT.D[2];
621  HXP[2] = F1*VECT.D[1] - F2*VECT.D[0];
622 
623  HP = F1*VECT.D[0] + F2*VECT.D[1] + F3*VECT.D[2];
624 //*
625  RHO1 = ONE/RHO;
626  SINT = sin(TET);
627  COST = TWO*pow(sin(HALF*TET),2);
628 //*
629  G1 = SINT*RHO1;
630  G2 = COST*RHO1;
631  G3 = (TET-SINT) * HP*RHO1;
632  G4 = -COST;
633  G5 = SINT;
634  G6 = COST * HP;
635 
636  VOUT.X[0] = VECT.X[0] + (G1*VECT.D[0] + G2*HXP[0] + G3*F1);
637  VOUT.X[1] = VECT.X[1] + (G1*VECT.D[1] + G2*HXP[1] + G3*F2);
638  VOUT.X[0] = VECT.X[2] + (G1*VECT.D[2] + G2*HXP[2] + G3*F3);
639 
640  VOUT.D[0] = VECT.D[0] + (G4*VECT.D[0] + G5*HXP[0] + G6*F1);
641  VOUT.D[1] = VECT.D[1] + (G4*VECT.D[1] + G5*HXP[1] + G6*F2);
642  VOUT.D[2] = VECT.D[2] + (G4*VECT.D[2] + G5*HXP[2] + G6*F3);
643 //*
644  } else {
645  VOUT.X[0] = VECT.X[0] + STEP*VECT.D[0];
646  VOUT.X[1] = VECT.X[1] + STEP*VECT.D[1];
647  VOUT.X[0] = VECT.X[2] + STEP*VECT.D[2];
648 //*
649  }//endif
650 //*
651  L999:;
652  fNStps = ITER;
653 
654 }
655 //______________________________________________________________________________
656 //______________________________________________________________________________
657 //______________________________________________________________________________
658 //______________________________________________________________________________
659 class TestMag: public TRKuttaMag
660 {
661  public:
662  TestMag (double mag[3]){memcpy (mMag,mag,sizeof(mMag));}
663  // methods
664  virtual void operator()(const double X[3], double B[3])
665  { memcpy(B,mMag,sizeof(mMag));};
666  protected:
667  // data members
668  double mMag[3];
669 };
670 
671 //______________________________________________________________________________
672 void TRungeKutta::Test(int flag)
673 {
674 // Flag == 3: 3 Dca test
675 // Flag == 2: 2 Dca test
676 
677  assert(flag==2 || flag==3);
678  double PtGev = 3.,Rho = 1./100, Hz = PtGev*Rho ;
679  double ZER[3]= {0,0,0};
680  double HZ[3] = {0,0,Hz};
681  double XZ[3] = {100,100,100};
682  double PZ[3] = {PtGev,0,PtGev*4./3};
683  double HH[3],XX[3],PP[3] ;
684 
685  TVector3 vPZ(PZ); vPZ.RotateZ(gRandom->Rndm()); vPZ.GetXYZ(PZ);
686  TCL::ucopy(HZ,HH,3);
687  TCL::ucopy(XZ,XX,3);
688  TCL::ucopy(PZ,PP,3);
689 
690  if (flag==3) {
691  double r1 = gRandom->Rndm(),r2 = gRandom->Rndm();
692  TVector3 vHZ(HZ);
693  vHZ.RotateX(r1); vHZ.RotateY(r2); vHZ.GetXYZ(HH);
694 
695  TVector3 vXZ(XZ);
696  vXZ.RotateX(r1); vXZ.RotateY(r2); vXZ.GetXYZ(XX);
697  TVector3 vPZ(PZ);
698  vPZ.RotateX(r1); vPZ.RotateY(r2); vPZ.GetXYZ(PP);
699  }
700  TestMag myMag(HH);
701  double s1,s2;
702  for (int charge=1;charge >=-1;charge-=2) {
703 
704  THelix3d TH(charge,ZER,PP,HH);
705  TRungeKutta T3(charge,ZER,PP,&myMag);
706  if (flag==3) {
707  s1 = TH.Path(XX);
708  s2 = T3.Path(XX);
709  } else {
710  s1 = TH.Path(XX[0],XX[1]);
711  s2 = T3.Path(XX[0],XX[1]);
712  }
713  TH.Move(s1); T3.Move(s2);
714  double dif[3];
715  TCL::vsub(TH.Pos(),XX,dif,3);
716  double ThTest = TCL::vdot(TH.Dir(),dif,flag);
717  TCL::vsub(T3.Pos(),XX,dif,3);
718  double TrTest = TCL::vdot(T3.Dir(),dif,flag);
719  printf ("Charge = %d ThTest,TrTest = %g/%g %g/%g\n",charge,ThTest,s1,TrTest,s2);
720  }
721 
722 }
723 //______________________________________________________________________________
724 class Test2Mag: public TRKuttaMag
725 {
726  public:
727  Test2Mag (double mag[3]){memcpy (mMag,mag,sizeof(mMag));}
728  // methods
729  virtual void operator()(const double X[3], double B[3])
730  {
731  double fak = 1./(1. + (X[0]*X[0]+X[1]*X[1]+X[2]*X[2]/100)*1e-5);
732  TCL::vscale(mMag,fak,B,3);};
733  int birth() { return 1946;}
734 
735  protected:
736  // data members
737  double mMag[3];
738 };
739 //______________________________________________________________________________
740 void TRungeKutta::Test2()
741 {
742  double PtGev = 3,Rho = 1./50, Hz = PtGev*Rho ;
743  double ZER[3]= {0,0,0};
744  double HZ[3] = {0,0,Hz};
745  double PZ[3] = {PtGev,0,PtGev/3*4};
746  TVector3 vPZ(PZ); vPZ.RotateZ(gRandom->Rndm()); vPZ.GetXYZ(PZ);
747 
748  Test2Mag myMag(HZ);
749 
750  for (int charge=1;charge >=-1;charge-=2) {
751  double myD[3]={PZ[0],PZ[1],PZ[2]};
752  double myX[3]={0,0,0};
753  double myH[3],s1=0;
754  double step = 1.;
755  int nSteps = 100./step+1e-6;
756  for (int is=0;is<nSteps;is++) {
757  myMag(myX,myH);
758  double myRho = -charge*myH[2]/PtGev;
759  THelixTrack_ TH(myX,myD,myRho);
760  TH.Move(step);
761  TCL::ucopy(TH.Pos(),myX,3);
762  TCL::ucopy(TH.Dir(),myD,3);
763  s1+=step;
764  }
765 
766  TRungeKutta T3(charge,ZER,PZ,&myMag);
767  double s2 = T3.Path(myX);
768  T3.Move();
769  const double *pos = T3.Pos();
770  printf("Helix.X = %g %g %g\n",myX[0],myX[1],myX[2]);
771  printf("RKuta.X = %g %g %g\n",pos[0],pos[1],pos[2]);
772  printf ("Charge = %d S1,S2 = %g %g \n",charge,s1,s2);
773  }
774 
775 }
776 //______________________________________________________________________________
777 void TRungeKutta::Test3()
778 {
779 // enum {kX=0,kY=1,kZ=2,kPx=0,kPy=1,kPz=2,kP0x=0,kP0y=1,kP0z=2};
780  double PtGev = 3.,Curv = 1./100, Hz = PtGev*Curv;
781  double HZ[3] = {0,0,Hz};
782  double Pbeg[3],Xbeg[3]={0},Pend[3];
783 
784  double X[200],Y[200],Z[200],A[200];
785 
786  TVector3 PbegV(PtGev,0,PtGev/3*4);
787  PbegV.RotateZ(gRandom->Rndm()*3.1415*2); PbegV.GetXYZ(Pbeg);
788  TVector3 XbegV(Xbeg);
789 
790  double Mag[3];
791  TVector3 magV(HZ);
792  double r1 =gRandom->Rndm(),r2 =gRandom->Rndm();
793  magV.RotateX(r1); magV.RotateY(r2); magV.GetXYZ(Mag);
794  Test2Mag myMag(Mag);
795  double stp = 3.1415/180/10;
796  int nJk = 100;
797  double myMax = -999, myMin=999;
798  double Plst[3];
799  for (int jk=-1;jk<nJk;jk++) {
800  PbegV.GetXYZ(Pbeg);
801  TRungeKutta TR0(1,Xbeg,Pbeg,&myMag);
802  TR0.Move(100.);
803  TCL::ucopy(TR0.Mom(),Pend,3);
804  if (jk>=0) {
805  A[jk] = jk*stp;
806  X[jk] = (Pend[0]-Plst[0])/stp;
807  Y[jk] = (Pend[1]-Plst[1])/stp;
808  Z[jk] = (Pend[2]-Plst[2])/stp;
809  if (myMax<X[jk]) myMax=X[jk];
810  if (myMax<Y[jk]) myMax=Y[jk];
811  if (myMax<Z[jk]) myMax=Z[jk];
812  if (myMin>X[jk]) myMin=X[jk];
813  if (myMin>Y[jk]) myMin=Y[jk];
814  if (myMin>Z[jk]) myMin=Z[jk];
815  }
816  PbegV.RotateY(stp);
817  TCL::ucopy(Pend,Plst,3);
818 
819  }
820 static TGraph *Gx=0,*Gy=0,*Gz=0;
821 static TCanvas *myCanvas=0;
822  delete Gx; delete Gy; delete Gz; delete myCanvas;
823  Gx = new TGraph(nJk,A,X); Gx->SetLineColor(kRed );Gx->SetMarkerColor(kRed );
824  Gy = new TGraph(nJk,A,Y); Gy->SetLineColor(kBlue );Gy->SetMarkerColor(kBlue );
825  Gz = new TGraph(nJk,A,Z); Gz->SetLineColor(kGreen);Gz->SetMarkerColor(kGreen);
826 
827 
828 
829  Gx->SetMaximum(myMax); Gx->SetMinimum(myMin);
830  Gy->SetMaximum(myMax); Gy->SetMinimum(myMin);
831  Gz->SetMaximum(myMax); Gz->SetMinimum(myMin);
832 
833  myCanvas = new TCanvas("TRungeKutta__Test3","",600,800);
834  Gx->Draw("AP");
835  Gy->Draw("same P");
836  Gz->Draw("same P");
837 
838  myCanvas->Modified();
839  myCanvas->Update();
840  while(!gSystem->ProcessEvents()){gSystem->Sleep(200);};
841 
842 }
843 //______________________________________________________________________________
844 void TRungeKutta::TestBak()
845 {
846 // enum {kX=0,kY=1,kZ=2,kPx=0,kPy=1,kPz=2,kP0x=0,kP0y=1,kP0z=2};
847  double PtGev = 3.,Curv = 1./100, Hz = PtGev*Curv;
848  TVector3 PbegV(PtGev,0,PtGev/3*4);
849  PbegV.RotateZ(gRandom->Rndm()*3.1415*2);
850  TVector3 XbegV(10,20,30),X1endV,P1endV;
851 
852  double L = 10.;
853 
854  TVector3 magV(0.,0.,Hz);
855  double r1 =gRandom->Rndm(),r2 =gRandom->Rndm();
856  magV.RotateX(r1); magV.RotateY(r2);
857  PbegV.RotateX(r1); PbegV.RotateY(r2);
858  Test2Mag myMag((double*)&magV[0]);
859 // TestMag myMag((double*)&magV[0]);
860 
861  for (int ichar=-1; ichar<=1; ichar +=2) {
862  auto DbegV = PbegV.Unit();
863  TRungeKutta TR0(ichar,(double*)&XbegV[0],(double*)&PbegV[0],&myMag);
864  TR0.Print("Normal");
865  TVector3 P0endV,X0endV,D0endV;
866  TVector3 P1endV,X1endV,D1endV;
867  TR0.SetDerOn();
868  TR0.Move(L);
869  TR0.Eval(0,&X0endV[0],&P0endV[0]);
870  D0endV = P0endV.Unit();
871  const THDer3d_t &totDer0 =*TR0.Der();
872 
873  TRungeKutta TR1(ichar,(double*)&XbegV[0],(double*)&PbegV[0],&myMag);
874  TR1.SetDerOn();
875  TR1.Backward();
876  TR1.Print("Negativ");
877  TR1.Move(-L);
878  TR1.Backward();
879  TR1.Eval(0,&X1endV[0],&P1endV[0]);
880  D1endV = P1endV.Unit();
881  const THDer3d_t &totDer1 =*TR1.Der();
882  TVector3 dif = X1endV-X0endV;
883  for (int i=0;i<3;i++) {
884  printf("X[%d] = %g \t%g \t%g\n",i,X0endV[i],X1endV[i],dif[i]);
885  }
886  dif = D1endV-D0endV;
887  for (int i=0;i<3;i++) {
888  printf("D[%d] = %g \t%g \t%g\n",i,D0endV[i],D1endV[i],dif[i]);
889  }
890 static const char* dp[5] = {"U","V","Fita","Lama","Pinv"};
891 
892  for (int i=0;i<5;i++) {
893  for (int j=0;j<5;j++) {
894  printf("Der[%s][%s] = %g \t%g \t%g\n",dp[i],dp[j]
895  ,totDer0[i][j],totDer1[i][j]
896  ,totDer1[i][j]-totDer0[i][j]);
897  } }
898  printf("\n");
899  }
900 }
901 
902 
903 
904 //______________________________________________________________________________
905 void TRungeKutta::TestDer()
906 {
907 // enum {kX=0,kY=1,kZ=2,kPx=0,kPy=1,kPz=2,kP0x=0,kP0y=1,kP0z=2};
908  double PtGev = 3.,Curv = 1./100, Hz = PtGev*Curv;
909  double tkSys[2][4][3]={{{0}}};
910  TVector3 PbegV(PtGev,0,PtGev/3*4);
911  PbegV.RotateZ(gRandom->Rndm()*3.1415*2);
912  TVector3 XbegV(10,20,30),X1endV,P1endV;
913 
914  double L = 100.;
915 
916  TVector3 magV(0.,0.,Hz);
917  double r1 =gRandom->Rndm(),r2 =gRandom->Rndm();
918  magV.RotateX(r1); magV.RotateY(r2);
919 // PbegV.RotateX(r1); PbegV.RotateY(r2);
920  Test2Mag myMag((double*)&magV[0]);
921 // TestMag myMag((double*)&magV[0]);
922 
923 // XbegV.Print("XbegV");
924 // PbegV.Print("PbegV");
925 // magV.Print("magV");
926  double totDif = 0;
927  for (int ichar=-1; ichar<=1; ichar +=2) {
928  double Pinv = -ichar/PbegV.Mag();
929  auto DbegV = PbegV.Unit();
930  TRungeKutta TR0(ichar,(double*)&XbegV[0],(double*)&PbegV[0],&myMag);
931  TVector3 PendV,XendV,DendV;
932  TR0.SetDerOn();
933  TR0.Move(L);
934  TCL::ucopy(TR0.Dir(),tkSys[1][kKT],3);
935  TR0.Eval(0,&XendV[0],&PendV[0]);
936  DendV = PendV.Unit();
937  const THDer3d_t &totDer =*TR0.Der();
938  TCL::ucopy(totDer.TkDir(0)[0],tkSys[0][0],2*4*3);
939 
940 
941  TVectorD der(5);
942 static const char* dp[5] = {"U","V","Fita","Lama","Pinv"};
943  for (int J=0;J<5;J++) {
944 
945  printf("\nTest d/d%s\n\n",dp[J]);
946  for (int i=0;i<5;i++) {der[i] = totDer[i][J];}
947  double delta = 1e-4;
948 
949  double P1inv = Pinv;
950  TVector3 D1begV = DbegV;
951  TVector3 X1begV = XbegV;
952  if (J<2) { X1begV += TVector3(tkSys[0][kKU+J ])*delta;}
953  else if (J<4) { D1begV += TVector3(tkSys[0][kKU+J-2])*delta;}
954  else { P1inv += delta ;}
955 
956  auto P1begV = D1begV; P1begV*=(1./fabs(P1inv));
957  TRungeKutta TR1(ichar,(double*)&X1begV[0],(double*)&P1begV[0],&myMag);
958  TR1.Eval(L,&X1endV[0],&P1endV[0]);
959  auto D1endV = P1endV.Unit();
960  TVectorD U1endV(5);
961  TVector3 difX = X1endV-XendV;
962  TVector3 difD = D1endV-DendV;
963  for (int j=0;j<5;j++) {
964  if (j<2) { U1endV[j] = difX.Dot(TVector3(tkSys[1][kKU+j ]));}
965  else if (j<4) { U1endV[j] = difD.Dot(TVector3(tkSys[1][kKU+j-2]));}
966  else { U1endV[4] = -ichar/P1endV.Mag()-Pinv ;}
967  }
968 
969  TVectorD dif = U1endV;
970  dif*=1./delta;
971  double pct=0; int ipct=-1;
972  for (int i=0;i<5;i++) {
973  if (fabs(dif[i])<1e-2) dif[i] = 0;
974  if (fabs(der[i])<1e-2) der[i] = 0;
975  double bott = 0.5*(fabs(dif[i])+fabs(der[i]));
976  if (bott<1e-3) continue;
977  double myPct = (dif[i]-der[i])/bott;
978  if (fabs(myPct)<fabs(pct)) continue;
979  pct = myPct; ipct = i;
980  }
981  if (pct >0.2) printf("################################ [%d] = %d%% \n",ipct,int(pct*100));
982  dif.Print("NUMERIC ");
983  der.Print("ANALITIC");
984  totDif+= (dif-der).Norm1();
985  }
986  } //end charge loop
987  printf ("========> TotalDiff = %g\n",totDif);
988 
989 
990 }
991 //______________________________________________________________________________
992 static double myDif(TVectorD &a,TVectorD &b, int &idx)
993 {
994 
995  int n = a.GetNrows();
996  double dif=0; idx = -1;
997  for (int i=0;i<n;i++) {
998  double mydif = fabs(a[i])+fabs(b[i]);
999  if (mydif<1e-4) continue;
1000  mydif = fabs(a[i]-b[i])/mydif;
1001  if (mydif <= dif) continue;
1002  dif = mydif; idx = i;
1003  }
1004  return dif;
1005 }
1006 //______________________________________________________________________________
1007 void TRungeKutta::TestDer2()
1008 {
1009 // Test for derivatives + changing direction
1010 
1011  double Ptot = 5, Pt=3;
1012 //double curv = 1./100,Hz = Pt*curv;
1013  double curv = 1./30,Hz = Pt*curv;
1014  TVector3 XbegV = TVector3(10,20,30);
1015  TVector3 PbegV = TVector3(Pt,0,Pt/3*4);
1016  double HZ[3] = {0,0,Hz};
1017 
1018 
1019  double L = 100.;
1020 
1021  double Mag[3];
1022  TVector3 magV(HZ); magV.GetXYZ(Mag );
1023 #if 1
1024  double r1=gRandom->Rndm(),r2=gRandom->Rndm();
1025  magV.RotateX(r1); magV.RotateY(r2); magV.GetXYZ(Mag );
1026  PbegV.RotateX(r1); PbegV.RotateY(r2);
1027  XbegV.RotateX(r1); XbegV.RotateY(r2);
1028 #endif
1029 // Test2Mag myMag((double*)&magV[0]);
1030  TestMag myMag((double*)&magV[0]);
1031 
1032 // XbegV.Print("XbegV");
1033 // PbegV.Print("PbegV");
1034 // magV.Print("magV");
1035 
1036  printf("===============PbegV:= \n"); PbegV.Print("");
1037 
1038  TVector3 XendV,PendV;
1039  const THDer3d_t *myDer = 0;
1040  for (int ichar=-1; ichar<=1; ichar +=2) {
1041  TRungeKutta TR0(ichar,(double*)&XbegV[0],(double*)&PbegV[0],&myMag);
1042  TR0.SetDerOn();
1043  TR0.Move(L);
1044  TR0.Backward(); //At the end TR0 helix inverted
1045  myDer = TR0.Der(); //Inverted derivatives
1046  TR0.Eval(0,&XendV[0],&PendV[0]); // Got X & P at the end
1047  auto DendV = PendV.Unit();
1048 
1049  TVector3 UbegV(myDer->TkDir(0)[kKU]); //Phi vector
1050  TVector3 VbegV(myDer->TkDir(0)[kKV]); //Lam vector
1051  TVector3 TbegV(myDer->TkDir(0)[kKT]); //Along track vector
1052 
1053  TVector3 UendV(myDer->TkDir(1)[kKU]); //Phi vector
1054  TVector3 VendV(myDer->TkDir(1)[kKV]); //Lam vector
1055  TVector3 TendV(myDer->TkDir(1)[kKT]); //Along track vector
1056  assert((DendV-TendV).Mag()<1e-3);
1057 
1058 
1059  TVectorD der(5);
1060  static const char* dp[5] = {"U","V","Fita","Lama","Pinv"};
1061  for (int J=0;J<5;J++) {
1062 
1063  printf("\nTest dFit/d%s\n\n",dp[J]);
1064  for (int i=0;i<5;i++) {der[i] = (*myDer)[i][J];}
1065 
1066 
1067  TRungeKutta TR1(ichar,(double*)&XbegV[0],(double*)&PbegV[0],&myMag);
1068  TR1.Backward();
1069  TVector3 P1begV(TR1.Mom());
1070  TVector3 X1begV(TR1.Pos());
1071  TkDir_t myTkDir;
1072  THelix3d::MakeTkDir((double*)&P1begV[0],Mag,myTkDir);
1073  double eps = (TVectorD(9,myTkDir[0])-TVectorD(9,myDer->TkDir(0)[0])).NormInf();
1074  assert(eps<1e-3);
1075 
1076 
1077 //?? double delta = 1e-4;
1078  double delta = 1e-6;
1079  TVector3 X1endV,P1endV;
1080  assert(TR0.Pinv()*ichar>0);
1081  assert(TR1.Pinv()*ichar>0);
1082  switch(J) {
1083  case kU:{
1084  X1begV += UbegV*delta;
1085  break;
1086  }
1087  case kV: {
1088  X1begV += VbegV*delta;
1089  break;
1090  }
1091  case kPinv: {
1092  double Pinv1 = TR1.Pinv()+delta;
1093  double Ptot1 = fabs(1./Pinv1);
1094  P1begV.SetMag(Ptot1);
1095  break;
1096  }
1097  case kLama: { //lamda variation (alfa)
1098  P1begV+=VbegV*(Ptot*delta); P1begV.SetMag(Ptot);
1099  break;
1100  }
1101 
1102  case kFita: { //Phi variation (beta)
1103  P1begV+=UbegV*(Ptot*delta); P1begV.SetMag(Ptot);
1104  break;
1105  }
1106  default: assert(0 && "Wrong J");
1107  }
1108  TR1.Set(0,(double*)&X1begV[0],(double*)&P1begV[0]);
1109  assert(TR1.Pinv()*ichar>0);
1110  TR1.Eval(-L,&X1endV[0],&P1endV[0]);
1111  TVectorD dif(5);
1112  TVector3 difX = X1endV-XendV;
1113  TVector3 difD = (P1endV.Unit()-PendV.Unit());
1114  dif[kU] = difX.Dot(UendV);
1115  dif[kV] = difX.Dot(VendV);
1116  dif[kPinv] = TR1.Pinv()-TR0.Pinv();
1117  dif[kFita] = difD.Dot(UendV);
1118  dif[kLama] = difD.Dot(VendV);
1119 
1120 
1121  assert(fabs(dif.NormInf())<1.);
1122  dif*=1./delta;
1123 
1124  int idx=-1;
1125  double maxa = myDif(dif,der,idx);
1126 
1127  if (maxa>1e-1) printf("########################[%d]= %g\n",idx,maxa);
1128  der.Print("ANALITIC");
1129  dif.Print("NUMERIC ");
1130  }
1131  } //end charge loop
1132 }
1133 //______________________________________________________________________________
1134 void TRungeKutta::TestErr(int charge)
1135 {
1136  double PtGev = 3.,Curv = 1./100, Hz = PtGev*Curv;
1137  double HZ[3] = {0,0,Hz};
1138  double Pbeg[3],P1beg[3],Xbeg[3]={0},X1beg[3];
1139  double Pend[3],P1end[3],Xend[3] ,X1end[3];
1140  double L = 300.; int nL = 1;
1141 
1142 // Canvas + Histograms
1143 enum {kNCanvs=10};
1144 static TCanvas *myCanvas[kNCanvs] = {0};
1145 static int myDiv[] = {2,5,5,5,5,0};
1146  for (int ic = 0;myDiv[ic];ic++) {
1147  if (myDiv[ic]<0) continue;
1148  TString ts("C"); ts+=ic;
1149  if (!myCanvas[ic]) myCanvas[ic] = new TCanvas(ts.Data(),ts.Data(),600,800);
1150  myCanvas[ic]->Clear();myCanvas[ic]->Divide(1,myDiv[ic]);
1151  }
1152 static TH1F *hXi2[2] = {0};
1153  for (int jk=0;jk<2;jk++) {
1154  TString ts("Xi2_"); ts += jk;
1155  delete hXi2[jk];
1156  hXi2[jk]= new TH1F(ts.Data(),ts.Data(),50,0,50);
1157  myCanvas[0]->cd(jk+1); hXi2[jk]->Draw();
1158  }
1159 // Now test of DCA
1160 static TH1F * hh[5]={0};
1161  for (int ih=0;ih<myDiv[1];ih++) {
1162 static const char *tit[5]={"U","V","Fita","Lama","Pinv"};
1163  delete hh[ih];
1164  hh[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1165  myCanvas[1]->cd(ih+1); hh[ih]->Draw();
1166  }
1167 // Now test of TRandomVector
1168 static TH1F * hrt[5]={0};
1169  for (int ih=0;ih<myDiv[1];ih++) {
1170 static const char *tit[5]={"rU","rV","rFita","rLama","rPinv"};
1171  delete hrt[ih];
1172  hrt[ih] = new TH1F(tit[ih],tit[ih],100,-10,+10);
1173  myCanvas[2]->cd(ih+1); hrt[ih]->Draw();
1174  }
1175 // Now test of Corr
1176 static TH1F * hcr[10]={0};
1177 {
1178 static const char *tit[]={"UV","UF","VF","UL","VL","FL","UP","VP","FP","LP"};
1179  for (int ih=0;ih<10;ih++) {
1180  delete hcr[ih];
1181  hcr[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1182  if (ih<5) myCanvas[3]->cd(ih+1 );
1183  else myCanvas[4]->cd(ih+1-5);
1184  hcr[ih]->Draw();
1185  }}
1186 //========================================================================
1187 
1188  TVector3 PbegV(PtGev,0,PtGev/3*4);
1189  double Ptot = PbegV.Mag();
1190  PbegV.RotateZ(gRandom->Rndm()*3.1415*2); PbegV.GetXYZ(Pbeg);
1191  TVector3 XbegV(Xbeg);
1192 
1193  double Mag[3]={HZ[0],HZ[1],HZ[2]};
1194  TVector3 magV(HZ);
1195  double r1 = gRandom->Rndm(),r2 = gRandom->Rndm();
1196  magV.RotateX(r1); magV.RotateY(r2); magV.GetXYZ(Mag);
1197  PbegV.RotateX(r1); PbegV.RotateY(r2); PbegV.GetXYZ(Pbeg);
1198  // Test2Mag myMag(Mag);
1199  TestMag myMag(Mag);
1200  TVectorD dia(5);
1201  dia[kU]= 0.1; dia[kV]= 0.2; dia[kPinv]= 0.1/Ptot, dia[kFita]= 1./360; dia[kLama]= 2./360;
1202  for(int jk=0;jk<5;jk++) {dia[jk]*=dia[jk];}
1203  TRandomVector RV(dia);
1204  auto &EMX = RV.GetMtx();
1205  auto &val = RV.GetLam();
1206  dia.Print("DIA");
1207  val.Print("VAL");
1208 
1209 
1210  double Gbeg[15],GbegD[5];
1211  for (int i=0,li=0;i< 5;li+=++i) {
1212  GbegD[i] = sqrt(EMX[i][i]);
1213  for (int j=0;j<=i;j++) {
1214  Gbeg[li+j] = EMX[i][j];
1215  } }
1216 
1217 
1218  double GbegI[15];
1219  TCL::trsinv(Gbeg,GbegI,5);
1220 
1221  TRungeKutta TH0(charge,Xbeg,Pbeg,&myMag);
1222  TH0.SetEmx(Gbeg);
1223 
1224  for (int split=0;split<nL;split++) {
1225  TH0.Move(L/nL);
1226  }
1227  double tkDir[3][3],tkDirE[3][3];
1228  memcpy(tkDir[0] ,TH0.TkDir(0)[0],sizeof(tkDir ));
1229  memcpy(tkDirE[0],TH0.TkDir(1)[0],sizeof(tkDirE));
1230 
1231  TH0.Eval(0,Xend,Pend);
1232  auto *myEmx = TH0.Emx();
1233  const double *Gend = *myEmx;
1234  double GendD[5];
1235  for (int i=0,li=0;i< 5;li+=++i) {
1236  GendD[i] = sqrt(Gend[li+i]);
1237  }
1238  double GendI[15];
1239  TCL::trsinv(Gend,GendI,5);
1240 
1241 
1242 // auto &tkDir = TH0.TkDir(0);
1243 // auto &tkDirE = TH0.TkDir(1);
1244  TVector3 PendV(Pend),XendV(Xend);
1245  TVectorD UendV(5);
1246  for (int j=0;j<3;j++) { UendV[j+2]=Pend[j];}
1247 
1248  double myG[15]={0};
1249  int nIter = 1000000;
1250  for (int iter=0;iter<nIter;iter++) {
1251  TVectorD delta = RV.Gaus();
1252  double chi2;
1253  TCL::trasat(delta.GetMatrixArray(),GbegI,&chi2,1,5);
1254  assert(chi2>0);
1255  hXi2[0]->Fill(chi2);
1256 
1257  for (int ih=0;ih<myDiv[1];ih++) { hrt[ih]->Fill(delta[ih]/GbegD[ih]);}
1258  double Pinv = -charge/PbegV.Mag(),P1inv = Pinv;
1259  TVector3 P1begV = PbegV,D1begV = P1begV.Unit();
1260  TVector3 X1begV = XbegV;
1261  for (int i=0;i<5;i++) {
1262  if (i<2) { X1begV += TVector3(tkDir[kKU+i ])*delta[i];}
1263  else if (i<4) { D1begV += TVector3(tkDir[kKU+i-2])*delta[i];}
1264  else { P1inv += delta[kPinv]; ;}
1265  }
1266  X1begV.GetXYZ(X1beg);
1267  P1begV = D1begV*(1./fabs(P1inv));
1268  P1begV.GetXYZ(P1beg);
1269  TRungeKutta TH1(charge,X1beg,P1beg,&myMag);
1270 
1272 
1273  double s = TH1.Path(Xend);
1274  TH1.Move(s);
1275  TH1.Eval(0,X1end,P1end);
1276  TVector3 X1endV(X1end),P1endV(P1end);
1277  TVectorD U1endV(5);
1278  TVector3 difX = X1endV-XendV;
1279  TVector3 difD = P1endV.Unit()-PendV.Unit();
1280  for (int j=0;j<5;j++) {
1281  if (j<2) { U1endV[j] = difX.Dot(TVector3(tkDirE[kKU+j ]));}
1282  else if (j<kPinv) { U1endV[j] = difD.Dot(TVector3(tkDirE[kKU+j-2]));}
1283  else { U1endV[4] = -charge/P1endV.Mag()-Pinv ;}
1284  }
1285 
1286  for (int ih=0;ih<5;ih++) { hh[ih]->Fill(U1endV[ih]/GendD[ih]);}
1287  TCL::trasat(U1endV.GetMatrixArray(),GendI,&chi2,1,5);
1288  hXi2[1]->Fill(chi2);
1289 
1290  for (int i=0,li=0;i<5;li+=++i) {
1291  for (int j=0;j<=i;j++) {
1292  myG[li+j]+=U1endV[i]*U1endV[j];
1293  } }
1294 
1295  int jk=0;
1296  for (int i=0,li=0;i<5;li+=++i) {
1297  for (int j=0;j<i;j++) {
1298  double d = U1endV[i]*U1endV[j] - Gend[li+j];
1299  d/= GendD[i]*GendD[j];
1300  hcr[jk++]->Fill(d);
1301  } }
1302 
1303 
1304  }
1305 
1306 
1307 
1308  TCL::vscale(myG,1./nIter,myG,15);
1309 
1310  printf("Numerical vs Analitical Error matrices:\n");
1311  for (int i=0,li=0;i<5;li+=++i) {
1312  for (int j=0;j<=i;j++) {
1313  printf("\t%g ",myG[li+j]);
1314  }
1315  printf("\n");
1316  for (int j=0;j<=i;j++) {
1317  printf("\t%g ",Gend[li+j]);
1318  }
1319  printf("\n");
1320  }
1321 
1322  for (int i=0,li=0;i<5;li+=++i) {
1323  for (int j=0;j<=i;j++) {
1324  printf("\t%g ",Gend[li+j]/(GendD[i]*GendD[j]));
1325  }
1326  printf("\n");
1327  for (int j=0;j<=i;j++) {
1328  printf("\t%g ",myG[li+j]/ (GendD[i]*GendD[j]));
1329  }
1330  printf("\n");
1331  }
1332 
1333 
1334 
1335  for (int i=0;myCanvas[i];i++) {
1336  if (!myCanvas[i]) continue;
1337  myCanvas[i]->Modified();myCanvas[i]->Update();}
1338 
1339  while(!gSystem->ProcessEvents()){gSystem->Sleep(200);};
1340 }
1341 
1342 //______________________________________________________________________________
1343 void TRungeKutta::TestErr2(int charge)
1344 {
1345  double PtGev = 1.,Curv = 1./100, Hz = PtGev*Curv;
1346  double HZ[3] = {0,0,Hz};
1347  double Pbeg[3],P1beg[3],Xbeg[3]={0},X1beg[3];
1348  double Pend[3],P1end[3],Xend[3] ,X1end[3];
1349 
1350 // Canvas + Histograms
1351 enum {kNCanvs=10};
1352 static TCanvas *myCanvas[kNCanvs] = {0};
1353 static int myDiv[] = {2,5,5,5,5,0};
1354  for (int ic = 0;myDiv[ic];ic++) {
1355  if (myDiv[ic]<0) continue;
1356  TString ts("C"); ts+=ic;
1357  if (!myCanvas[ic]) myCanvas[ic] = new TCanvas(ts.Data(),ts.Data(),600,800);
1358  myCanvas[ic]->Clear();myCanvas[ic]->Divide(1,myDiv[ic]);
1359  }
1360 static TH1F *hXi2[2] = {0};
1361  for (int jk=0;jk<2;jk++) {
1362  TString ts("Xi2_"); ts += jk;
1363  delete hXi2[jk];
1364  hXi2[jk]= new TH1F(ts.Data(),ts.Data(),50,0,0);
1365  myCanvas[0]->cd(jk+1); hXi2[jk]->Draw();
1366  }
1367 // Now test of DCA
1368 static TH1F * hh[5]={0};
1369  for (int ih=0;ih<myDiv[1];ih++) {
1370 const char * tit[]={"U","V","Pinv","Fita","Lama",0};
1371  delete hh[ih];
1372  hh[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1373  myCanvas[1]->cd(ih+1); hh[ih]->Draw();
1374  }
1375 // Now test of TRandomVector
1376 static TH1F * hrt[5]={0};
1377  for (int ih=0;ih<myDiv[1];ih++) {
1378 const char * tit[]={"U0","V0","Pinv0","Fita0","Lama0",0};
1379  delete hrt[ih];
1380  hrt[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1381  myCanvas[2]->cd(ih+1); hrt[ih]->Draw();
1382  }
1383 // Now test of Corr
1384 static TH1F * hcr[10]={0};
1385 {
1386 static const char *tit[]={"UV","UP","VP","UL","VL","PL","UF","VF","PF","LF"};
1387  for (int ih=0;ih<10;ih++) {
1388  delete hcr[ih];
1389  hcr[ih] = new TH1F(tit[ih],tit[ih],100,0,0);
1390  if (ih<5) myCanvas[3]->cd(ih+1 );
1391  else myCanvas[4]->cd(ih+1-5);
1392  hcr[ih]->Draw();
1393  }}
1394 //============================================================================
1395 
1396  TVector3 PbegV(PtGev,0,PtGev*3);
1397  double Ptot = PbegV.Mag();
1398  PbegV.RotateZ(gRandom->Rndm()*3.1415*2); PbegV.GetXYZ(Pbeg);
1399  TVector3 XbegV(Xbeg);
1400 
1401  double L = 100.;
1402 
1403  double Mag[3]={HZ[0],HZ[1],HZ[2]};
1404 #if 11
1405  TVector3 magV(HZ);
1406  double r1 = gRandom->Rndm(),r2 = gRandom->Rndm();
1407  magV.RotateX(r1); magV.RotateY(r2); magV.GetXYZ(Mag);
1408  PbegV.RotateX(r1); PbegV.RotateY(r2); PbegV.GetXYZ(Pbeg);
1409 #endif
1410  Test2Mag myMag(Mag);
1411 
1412 
1413 
1414  TVectorD dia(5);
1415  dia[kU]= 0.1; dia[kV]= 0.2; dia[kPinv]= 0.1/Ptot; dia[kFita]= 1./360; dia[kLama]= 2./360;
1416 
1417  TRandomVector RV(dia);
1418  auto &EMX = RV.GetMtx();
1419  auto &val = RV.GetLam();
1420  dia.Print("DIA");
1421  val.Print("VAL");
1422 
1423 
1424  double Gbeg[15],GbegD[5];
1425  for (int i=0,li=0;i< 5;li+=++i) {
1426  GbegD[i] = sqrt(EMX[i][i]);
1427  for (int j=0;j<=i;j++) {
1428  Gbeg[li+j] = EMX[i][j];
1429  } }
1430 
1431 
1432 
1433 
1434  TRungeKutta TH0(charge,Xbeg,Pbeg,&myMag);
1435 //==================================
1436  THEmx3d_t *emx = new THEmx3d_t(Gbeg);
1437 
1438  TRungeKutta TH0i(&TH0);
1439  TH0i.SetEmx(emx);
1440  TH0i.Backward();
1441  double GbegI[15];
1442  memcpy(Gbeg,*TH0i.Emx(),sizeof(Gbeg));
1443  TRandomVector RVi(5,Gbeg);
1444  TCL::trsinv(Gbeg,GbegI,5);
1445 
1446  TH0.SetEmx(emx);
1447  TH0.Move(L);
1448  TH0.Backward();
1449 
1450  double tkDir[3][3],tkDirE[3][3];
1451  memcpy(tkDir[0] ,TH0.TkDir(0)[0],sizeof(tkDir ));
1452  memcpy(tkDirE[0],TH0.TkDir(1)[0],sizeof(tkDirE));
1453 
1454  TH0.Eval(0,Xend,Pend);
1455  auto *myEmx = TH0.Emx();
1456  const double *Gend = *myEmx;
1457  double GendD[5];
1458  for (int i=0,li=0;i< 5;li+=++i) {
1459  GendD[i] = sqrt(Gend[li+i]);
1460  }
1461 
1462  double GendI[15];
1463  TCL::trsinv(Gend,GendI,5);
1464 
1465  TVector3 PendV(Pend),XendV(Xend);
1466  TVectorD UendV(5);
1467  for (int j=0;j<3;j++) { UendV[j+2]=Pend[j];}
1468 
1469  double myG[15]={0};
1470  int nIter = 1000000;
1471  for (int iter=0;iter<nIter;iter++) {
1472  TVectorD delta = RVi.Gaus();
1473  double chi2;
1474  TCL::trasat(delta.GetMatrixArray(),GbegI,&chi2,1,5);
1475  assert(chi2>0);
1476  hXi2[0]->Fill(chi2);
1477 
1478  for (int ih=0;ih<myDiv[1];ih++) { hrt[ih]->Fill(delta[ih]/GbegD[ih]);}
1479  TVector3 P1begV(TH0i.Mom());
1480  TVector3 X1begV(TH0i.Pos());
1481 
1482  X1begV += TVector3(tkDir[kKU])*delta[kU];
1483  X1begV += TVector3(tkDir[kKV])*delta[kV];
1484  double Pinv = TH0i.Pinv();
1485  double Pinv1 = Pinv + delta[kPinv];
1486  double Ptot1 = fabs(1./Pinv1);
1487  P1begV.SetMag(Ptot1);
1488  P1begV += TVector3(tkDir[kKV])*(delta[kLama]*Ptot1);P1begV.SetMag(Ptot1);
1489  P1begV += TVector3(tkDir[kKU])*(delta[kFita]*Ptot1);P1begV.SetMag(Ptot1);
1490 
1491  X1begV.GetXYZ(X1beg);
1492  P1begV.GetXYZ(P1beg);
1493  TRungeKutta TH1(TH0i.Charge(),X1beg,P1beg,&myMag);
1494  TH1.Move(-L);
1495  TH1.Eval(0,X1end,P1end);
1496  TVector3 X1endV(X1end),P1endV(P1end);
1497  TVectorD U1endV(5);
1498  TVector3 difX = X1endV-XendV;
1499  TVector3 difD = (P1endV.Unit()-PendV.Unit());
1500 
1501  U1endV[kU] = difX.Dot(TVector3(tkDirE[kKU]));
1502  U1endV[kV] = difX.Dot(TVector3(tkDirE[kKV]));
1503  U1endV[kFita] = difD.Dot(TVector3(tkDirE[kKU]));
1504  U1endV[kLama] = difD.Dot(TVector3(tkDirE[kKV]));
1505  U1endV[kPinv] = (Pinv1-Pinv);
1506 
1507  for (int ih=0;ih<5;ih++) { hh[ih]->Fill(U1endV[ih]/GendD[ih]);}
1508  TCL::trasat(U1endV.GetMatrixArray(),GendI,&chi2,1,5);
1509  hXi2[1]->Fill(chi2);
1510 
1511  for (int i=0,li=0;i<5;li+=++i) {
1512  for (int j=0;j<=i;j++) {
1513  myG[li+j]+=U1endV[i]*U1endV[j];
1514  } }
1515 
1516  int jk=0;
1517  for (int i=0,li=0;i<5;li+=++i) {
1518  for (int j=0;j<i;j++) {
1519  double d = U1endV[i]*U1endV[j] - Gend[li+j];
1520  d/= GendD[i]*GendD[j];
1521  hcr[jk++]->Fill(d);
1522  } }
1523  }
1524 
1525  TCL::vscale(myG,1./nIter,myG,15);
1526 
1527  printf("Numerical vs Analitical Error matrices:\n");
1528  for (int i=0,li=0;i<5;li+=++i) {
1529  for (int j=0;j<=i;j++) {
1530  printf("\t%g ",myG[li+j]);
1531  }
1532  printf("\n");
1533  for (int j=0;j<=i;j++) {
1534  printf("\t%g ",Gend[li+j]);
1535  }
1536  printf("\n");
1537  }
1538 
1539 
1540  for (int i=0;myCanvas[i];i++) {
1541  if (!myCanvas[i]) continue;
1542  myCanvas[i]->Modified();myCanvas[i]->Update();}
1543 
1544  while(!gSystem->ProcessEvents()){gSystem->Sleep(200);};
1545 }
1546 //______________________________________________________________________________
1547 //______________________________________________________________________________
1548 void TRungeKutta::TestSign()
1549 {
1550 
1551  double PtGev = 3.,Rho = 1./100, Hz = PtGev*Rho ;
1552  double ZER[3]= {0,0,0};
1553  double HZ[3] = {0,0,Hz};
1554  double PZ[3] = {PtGev,0,PtGev*4./3};
1555 
1556  for (int charge = -1 ;charge<=1;charge +=2 ) {
1557  for (double L =-100 ;L<=100 ;L +=200) {
1558  for (int ih =-1 ;ih<=1 ;ih +=2 ) {
1559  HZ[2] = fabs(HZ[2])*ih;
1560  TestMag myMag(HZ);
1561  double curv = -(HZ[2]*charge)/PtGev;
1562 
1563  THelixTrack_ th(ZER,PZ,curv);
1564  THelix3d t3(charge,ZER,PZ,HZ);
1565  TRungeKutta tr(charge,ZER,PZ,&myMag);
1566 
1567  th.Move(L);
1568  t3.Move(L);
1569  tr.Move(L);
1570 
1571  printf("charge = %d \tL = %g , THelixTreack \tX=%g \tY=%g \tZ=%g \tDx=%g \tDy=%g \tDz=%g\n"
1572  , charge,L,th.Pos()[0],th.Pos()[1],th.Pos()[2]
1573  , th.Dir()[0],th.Dir()[1],th.Dir()[2]);
1574  printf("charge = %d \tL = %g , THelix3d \tX=%g \tY=%g \tZ=%g \tDx=%g \tDy=%g \tDz=%g\n"
1575  , charge,L,t3.Pos()[0],t3.Pos()[1],t3.Pos()[2]
1576  , t3.Dir()[0],t3.Dir()[1],t3.Dir()[2]);
1577  printf("charge = %d \tL = %g , TRungeKutta \tX=%g \tY=%g \tZ=%g \tDx=%g \tDy=%g \tDz=%g\n"
1578  , charge,L,tr.Pos()[0],tr.Pos()[1],tr.Pos()[2]
1579  , tr.Dir()[0],tr.Dir()[1],tr.Dir()[2]);
1580  printf("\n");
1581  } } }
1582 }
TRungeKutta::Move
double Move()
Move to point founded before.
Definition: TRungeKutta.cxx:158
TRungeKutta::Eval
double Eval(double step, double xyz[3]=0, double mom[3]=0) const
Evaluate params with given step along helix.
Definition: TRungeKutta.cxx:143
TRungeKutta::Dca
double Dca(double x, double y, double *dcaErr=0) const
DCA to given 2dim point (with error matrix)
Definition: TRungeKutta.cxx:348
TkDir_t
Definition: THelix3d.h:15
THEmx3d_t
Definition: THelix3d.h:62
TRungeKutta::TestErr
static void TestErr(int charge=1)
Definition: TRungeKutta.cxx:1134
TRungeKutta::Backward
void Backward()
Change direction.
Definition: TRungeKutta.cxx:110
TCL::trsinv
static float * trsinv(const float *g, float *gi, int n)
Definition: TCernLib.cxx:1160
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
TRungeKutta::Move
double Move(double step)
Move along track.
Definition: TRungeKutta.cxx:121
THDer3d_t
Definition: THelix3d.h:40
TRungeKutta::Path
double Path(const double point[3], int idir=0, int ndca=3) const
Definition: TRungeKutta.cxx:167
Generated by   doxygen 1.8.5