EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FairRKPropagator.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FairRKPropagator.cxx
1 #include "FairRKPropagator.h"
2 
3 
4 #include "TMath.h"
5 #include "TVector3.h"
6 
7 ClassImp(FairRKPropagator);
8 
9 //______________________________________________________________________________
10 FairRKPropagator::FairRKPropagator(FairField* field)
11  : TObject(),
12  fMaxStep(10.0),
13  fMagField (field)
14 {
15  // fMaxStep=10.0;
16 }
17 //______________________________________________________________________________
18 FairRKPropagator::~FairRKPropagator()
19 {
20  // Destructor.
21 }
22 //______________________________________________________________________________
23 void FairRKPropagator::PropagatToPlane(Double_t Charge, Double_t* vecRKIn, Double_t* vec1, Double_t* vec2, Double_t* vec3, Double_t* vecOut)
24 {
29  Double_t Norm[3];
30  Double_t Mag;
31  Double_t dist[3];
32  Double_t distance[3];
33  Double_t vecRKoutT[7];
34 
35  Norm[0]=vec1[1]*vec2[2] - vec2[2] * vec2[1]; // a2b3 − a3b2,
36  Norm[1]=vec1[2]*vec2[0] - vec1[0] * vec2[2]; // a3b1 − a1b3;
37  Norm[2]=vec1[0]*vec2[1] - vec1[1] * vec2[0];
38 
39  Mag=TMath::Sqrt(Norm[0]*Norm[0]+Norm[1]*Norm[1]+Norm[2]*Norm[2]);
40 
41  // printf(" Mag = %f \n ", Mag);
42 
43  Norm[0]=Norm[0]/Mag;
44  Norm[1]=Norm[1]/Mag;
45  Norm[2]=Norm[2]/Mag;
46  // printf(" after normalization : Normal = %f %f %f \n ", Norm[0],Norm[1],Norm[2]);
47 
48  dist[0]=vecRKIn[0]-vec3[0];
49  dist[1]=vecRKIn[1]-vec3[1];
50  dist[2]=vecRKIn[2]-vec3[2];
51 
52  distance[0]=Norm[0]*dist[0];
53  distance[1]=Norm[1]*dist[1];
54  distance[2]=Norm[2]*dist[2];
55 // printf(" distance = %f %f %f \n ", distance[0],distance[1],distance[2]);
56  Double_t diff = TMath::Abs(distance[0]+distance[1]+distance[2]);
57  fMaxStep = diff;
58  Double_t res = 100.0;
59  Double_t res_old = 100.0;
60 
61  Double_t vecRKOut[7];
62  // for (Int_t i=0; i< 7; i++) {vecRKOut[i]=0;}
63  Int_t nIter=0;
64 
65  // printf("I am in CPU code %f %f %f res= %f diff = %f \n ", vecRKIn[0], vecRKIn[1],vecRKIn[2], res, diff);
66 
67  do {
68  Step(Charge,vecRKIn,vecRKOut);
69  dist[0]=(vecRKOut[0]-vec3[0])*Norm[0];
70  dist[1]=(vecRKOut[1]-vec3[1])*Norm[1];
71  dist[2]=(vecRKOut[2]-vec3[2])*Norm[2];
72  fMaxStep=TMath::Sqrt(dist[0]*dist[0]+dist[1]*dist[1]+dist[2]*dist[2]);
73  res=TMath::Abs(fMaxStep/diff);
74  // printf("After %i step %f %f %f res = %f \n", nIter ,vecRKOut[0], vecRKOut[1],vecRKOut[2] , res);
75  if( res< 0.001 || res >res_old ) {
76  break;
77  } else {
78  for (Int_t i=0; i< 3; i++) {
79  vecRKIn[i]=vecRKOut[i];
80  vecRKoutT[i]=vecRKOut[i];
81  }
82  res_old=res;
83  }
84  if(nIter++>1000) { break; }
85  } while(1);
86 // printf("The results is %f %f %f , no of iter %i \n", vecRKOut[0],vecRKOut[1],vecRKOut[2], nIter);
87 // printf("\n");
88  for (Int_t i=0; i< 3; i++) {
89  if (res > res_old) { vecOut[i]=vecRKoutT[i]; }
90  else { vecOut[i]=vecRKOut[i]; }
91  }
92 }
93 //______________________________________________________________________________
94 void FairRKPropagator::Propagat(Double_t Charge, Double_t* vecRKIn, Double_t* Pos)
95 {
96  Double_t diff = Pos[2] - vecRKIn[2];
97  fMaxStep = diff/25;
98  Double_t res_old= diff;
99  Double_t res = 100.0;
100  Double_t vecRKOut[7];
101  Double_t vecRKOutT[7];
102  for (Int_t i=0; i< 7; i++) {vecRKOut[i]=0; vecRKOutT[i]=0;}
103 
104  Int_t nIter=0;
105  do {
106  Step(Charge,vecRKIn,vecRKOut);
107  res=(vecRKOut[2]-Pos[2])/diff;
108  if( TMath::Abs(res)< 0.01 || res >res_old ) {
109  break;
110  } else {
111  for (Int_t i=0; i< 3; i++) {
112  vecRKOutT[i]=vecRKOut[i];
113  vecRKIn[i]=vecRKOut[i];
114  }
115  }
116  if(nIter++>1000) { break; }
117  } while(1);
118  if (res > res_old) for (Int_t k=0; k< 7; k++) { vecRKOut[k]=vecRKOutT[k]; }
119  for (Int_t k=0; k< 3; k++) { printf(" vecRKOut[%i] =%f ", k, vecRKOut[k] ); }
120  printf("\n");
121 }
122 //______________________________________________________________________________
123 void FairRKPropagator::Step(Double_t Charge, Double_t* vecRKIn, Double_t* vecOut)
124 {
125 
126  Double_t vecRKOut[7];
127  for (Int_t i=0; i< 7; i++) { vecRKOut[i]=0; }
128 // for (Int_t i=0; i< 7; i++) printf( "vectRKIn(%i)=%f \n",i ,vecRKIn[i]);
129 // printf(" ---------------------------------------------------------------- \n");
130  OneStepRungeKutta(Charge,fMaxStep, vecRKIn, vecRKOut);
131  // printf(" now at x=%f y=%f z=%f \n", vecRKOut[0],vecRKOut[1],vecRKOut[2]);
132  vecOut[0] = vecRKOut[0];
133  vecOut[1] = vecRKOut[1];
134  vecOut[2] = vecRKOut[2];
135  vecOut[6] = vecRKOut[6];
136  vecOut[3] = vecRKOut[3]/vecOut[6];
137  vecOut[4] = vecRKOut[4]/vecOut[6];
138  vecOut[5] = vecRKOut[5]/vecOut[6];
139 
140 }
141 
142 //______________________________________________________________________________
143 void FairRKPropagator::OneStepRungeKutta(Double_t charge, Double_t step,
144  Double_t* vect, Double_t* vout)
145 {
146 
147  // Wrapper to step with method RungeKutta.
148 
170 
171  Double_t h2, h4, f[4];
172  Double_t xyzt[3], a=0, b=0, c=0, ph,ph2;
173  Double_t secxs[4],secys[4],seczs[4]; //hxp[3];
174  Double_t /*g1 , g2, g3, g4, g5, g6,*/ ang2, dxt, dyt, dzt;
175  Double_t est, at, bt, ct, cba;
176 // Double_t /*f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost*/;
177 
178  Double_t x=0;
179  Double_t y=0;
180  Double_t z=0;
181 
182  Double_t xt;
183  Double_t yt;
184  Double_t zt;
185 
186 // Double_t maxit = 1992;
187  Double_t maxit = 10;
188  Double_t maxcut = 11;
189 
190  const Double_t hmin = 1e-4; // !!! MT ADD, should be member
191  const Double_t kdlt = 1e-3; // !!! MT CHANGE from 1e-4, should be member
192  const Double_t kdlt32 = kdlt/32.;
193  const Double_t kthird = 1./3.;
194  const Double_t khalf = 0.5;
195  const Double_t kec = 2.9979251e-3;
196  const Double_t kpisqua = 9.86960440109;
197  /* const Int_t kix = 0;
198  const Int_t kiy = 1;
199  const Int_t kiz = 2;
200  const Int_t kipx = 3;
201  const Int_t kipy = 4;
202  const Int_t kipz = 5;
203  */
204  // *.
205  // *. ------------------------------------------------------------------
206  // *.
207  // * this constant is for units cm,gev/c and kgauss
208  // *
209  Int_t iter = 0;
210  Int_t ncut = 0;
211  for(Int_t j = 0; j < 7; j++) {
212  vout[j] = vect[j];
213  }
214 
215  Double_t pinv = kec * charge / vect[6];
216  Double_t tl = 0.;
217  Double_t h = step;
218  Double_t rest;
219 
220  do {
221  rest = step - tl;
222  printf(" Step no. %i x=%f y=%f z=%f px/p = %f py/p =%f pz/p= %f \n", iter, x,y,z,a,b,c);
223  if (TMath::Abs(h) > TMath::Abs(rest)) {
224  h = rest;
225  }
226 
227 
228  fMagField->GetFieldValue( vout, f);
229  f[0] = -1.0*f[0];
230  f[1] = -1.0*f[1];
231  f[2] = -1.0*f[2];
232 
233 
234 // printf(" Field Values x=%f y=%f z=%f \n", f[0],f[1],f[2]);
235 // f[0] = -fH.fB.fX;
236 // f[1] = -fH.fB.fY;
237 // f[2] = -fH.fB.fZ;
238 
239  // * start of integration
240  x = vout[0];
241  y = vout[1];
242  z = vout[2];
243  a = vout[3];
244  b = vout[4];
245  c = vout[5];
246 
247  h2 = khalf * h;
248  h4 = khalf * h2;
249  ph = pinv * h;
250  ph2 = khalf * ph;
251 
252 // printf(" ------------------------------------------- h2 = %f\n",h2);
253 
254  secxs[0] = (b * f[2] - c * f[1]) * ph2;
255  secys[0] = (c * f[0] - a * f[2]) * ph2;
256  seczs[0] = (a * f[1] - b * f[0]) * ph2;
257  ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]);
258  if (ang2 > kpisqua) { break; }
259 
260  dxt = h2 * a + h4 * secxs[0];
261  dyt = h2 * b + h4 * secys[0];
262  dzt = h2 * c + h4 * seczs[0];
263  xt = x + dxt;
264  yt = y + dyt;
265  zt = z + dzt;
266 // printf(" Position 1 at xt=%f yt=%f zt=%f \n", xt, yt, zt);
267 // printf(" differance dxt=%f dyt=%f dzt=%f \n", dxt, dyt, dzt);
268  // * second intermediate point
269  est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
270  if (est > h) {
271  if (ncut++ > maxcut) { break; }
272  h *= khalf;
273  continue;
274  }
275 
276  xyzt[0] = xt;
277  xyzt[1] = yt;
278  xyzt[2] = zt;
279 
280  fMagField->GetFieldValue( xyzt, f);
281 
282  f[0] = -f[0];
283  f[1] = -f[1];
284  f[2] = -f[2];
285 
286 
287 // printf(" Field Values at x=%f y=%f z=%f , Bx=%f By=%f Bz=%f \n", xyzt[0], xyzt[1], xyzt[2] ,f[0],f[1],f[2]);
288 // fH.fB = fMagFieldObj->GetField(xt, yt, zt);
289 // f[0] = -fH.fB.fX;
290 // f[1] = -fH.fB.fY;
291 // f[2] = -fH.fB.fZ;
292 
293 
294 
295 
296  at = a + secxs[0];
297  bt = b + secys[0];
298  ct = c + seczs[0];
299 
300  secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
301  secys[1] = (ct * f[0] - at * f[2]) * ph2;
302  seczs[1] = (at * f[1] - bt * f[0]) * ph2;
303  at = a + secxs[1];
304  bt = b + secys[1];
305  ct = c + seczs[1];
306  secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
307  secys[2] = (ct * f[0] - at * f[2]) * ph2;
308  seczs[2] = (at * f[1] - bt * f[0]) * ph2;
309  dxt = h * (a + secxs[2]);
310  dyt = h * (b + secys[2]);
311  dzt = h * (c + seczs[2]);
312  xt = x + dxt;
313  yt = y + dyt;
314  zt = z + dzt;
315  at = a + 2.*secxs[2];
316  bt = b + 2.*secys[2];
317  ct = c + 2.*seczs[2];
318  // printf(" Position 2 at xt=%f yt=%f zt=%f \n", xt, yt, zt);
319 
320  est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt);
321  if (est > 2.*TMath::Abs(h)) {
322  if (ncut++ > maxcut) { break; }
323  h *= khalf;
324  continue;
325  }
326 
327  xyzt[0] = xt;
328  xyzt[1] = yt;
329  xyzt[2] = zt;
330 
331  fMagField->GetFieldValue( xyzt, f);
332  f[0] = -1.0*f[0];
333  f[1] = -1.0*f[1];
334  f[2] = -1.0*f[2];
335 
336  // fH.fB = fMagFieldObj->GetField(xt, yt, zt);
337  // f[0] = -fH.fB.fX;
338  // f[1] = -fH.fB.fY;
339  // f[2] = -fH.fB.fZ;
340 
341  z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
342  y = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
343  x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
344  // printf(" Position 3 at x=%f y=%f z=%f \n", x, y, z);
345  secxs[3] = (bt*f[2] - ct*f[1])* ph2;
346  secys[3] = (ct*f[0] - at*f[2])* ph2;
347  seczs[3] = (at*f[1] - bt*f[0])* ph2;
348  a = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
349  b = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
350  c = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;
351 
352  est = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
353  + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
354  + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));
355 
356  if (est > kdlt && TMath::Abs(h) > hmin) {
357  if (ncut++ > maxcut) { break; }
358  h *= khalf;
359  continue;
360  }
361 
362  ncut = 0;
363  // * if too many iterations, go to helix
364  if (iter++ > maxit) { break; }
365 
366  tl += h;
367  if (est < kdlt32) {
368  h *= 2.;
369  }
370  cba = 1./ TMath::Sqrt(a*a + b*b + c*c);
371  vout[0] = x;
372  vout[1] = y;
373  vout[2] = z;
374  vout[3] = cba*a;
375  vout[4] = cba*b;
376  vout[5] = cba*c;
377 
378 
379  rest = step - tl;
380  // printf(" Position 4 at x=%f y=%f z=%f Step = %f \n", x, y, z, step );
381 
382  if (step < 0.) { rest = -rest; }
383  if (rest < 1.e-5*TMath::Abs(step)) { return; }
384 
385  } while(1);
386 
387  // angle too big, use helix
388  /*
389  f1 = f[0];
390  f2 = f[1];
391  f3 = f[2];
392  f4 = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
393  rho = -f4*pinv;
394  tet = rho * step;
395 
396  hnorm = 1./f4;
397  f1 = f1*hnorm;
398  f2 = f2*hnorm;
399  f3 = f3*hnorm;
400 
401  hxp[0] = f2*vect[kipz] - f3*vect[kipy];
402  hxp[1] = f3*vect[kipx] - f1*vect[kipz];
403  hxp[2] = f1*vect[kipy] - f2*vect[kipx];
404 
405  hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];
406 
407  rho1 = 1./rho;
408  sint = TMath::Sin(tet);
409  cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);
410 
411  g1 = sint*rho1;
412  g2 = cost*rho1;
413  g3 = (tet-sint) * hp*rho1;
414  g4 = -cost;
415  g5 = sint;
416  g6 = cost * hp;
417 
418  vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
419  vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
420  vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
421 
422  vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
423  vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
424  vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
425  */
426 
427 }