FairRoot
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FairRKPropagator.cxx
Go to the documentation of this file.
1 /********************************************************************************
2  * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH *
3  * *
4  * This software is distributed under the terms of the *
5  * GNU Lesser General Public Licence (LGPL) version 3, *
6  * copied verbatim in the file "LICENSE" *
7  ********************************************************************************/
8 #include "FairRKPropagator.h"
9 
10 #include "FairField.h" // for FairField
11 #include "FairTrackParH.h" // for FairTrackParH
12 #include "FairTrackParP.h" // for FairTrackParP
13 #include "TDatabasePDG.h" // for TDatabasePDG
14 #include "TMath.h" // for Sqrt
15 #include "TMathBase.h" // for Abs
16 #include "TVector3.h" // for TVector3, operator-, operator*
17 
18 #include <TGenericClassInfo.h> // for TGenericClassInfo
19 #include <TParticlePDG.h> // for TParticlePDG
20 #include <fairlogger/Logger.h> // for LOG
21 #include <math.h> // for sqrt
22 #include <stdio.h> // for printf
23 #include <stdlib.h> // for abs
24 
25 ClassImp(FairRKPropagator);
26 
27 //______________________________________________________________________________
28 FairRKPropagator::FairRKPropagator(FairField* field)
29  : FairPropagator("FairRKPropagator", "Runge-Kutta propagator")
30  , fMaxStep(10.0)
31  , fMagField(field)
32  , fPropagationFlag(NONE)
33  , fDefPlaneV0()
34  , fDefPlaneV1()
35  , fDefPlaneV2()
36  , fPCAPropagationType(0)
37  , fPCAPropagationDir(1)
38  , fPCAPropagationPar(nullptr)
39 
40 {
41  // fMaxStep=10.0;
42 }
43 //______________________________________________________________________________
44 FairRKPropagator::~FairRKPropagator()
45 {
46  // Destructor.
47 }
48 
49 double FairRKPropagator::GetChargeFromPDG(int pdg)
50 {
51  TDatabasePDG* pdgDB = TDatabasePDG::Instance();
52  return pdgDB->GetParticle(pdg)->Charge() / 3.; // Charge(): charge in units of |e|/3
53 }
54 
55 bool FairRKPropagator::Propagate(FairTrackParH* TStart, FairTrackParH* TEnd, int PDG)
56 {
57  if (fPropagationFlag == TOPLANE) {
58  float x1[3] = {(float)TStart->GetX(), (float)TStart->GetY(), (float)TStart->GetZ()};
59  float p1[3] = {(float)TStart->GetPx(), (float)TStart->GetPy(), (float)TStart->GetPz()};
60  float x2[3] = {0., 0., 0.};
61  float p2[3] = {0., 0., 0.};
62  bool ret = Propagate(x1, p1, x2, p2, PDG);
63  if (!ret)
64  return false;
65  TEnd->SetX(x2[0]);
66  TEnd->SetY(x2[1]);
67  TEnd->SetZ(x2[2]);
68  TEnd->SetPx(p2[0]);
69  TEnd->SetPy(p2[1]);
70  TEnd->SetPz(p2[2]);
71  TEnd->SetQp(GetChargeFromPDG(PDG) / sqrt(p2[0] * p2[0] + p2[1] * p2[1] + p2[2] * p2[2]));
72  TEnd->SetDX(0.);
73  TEnd->SetDY(0.);
74  TEnd->SetDZ(0.);
75  TEnd->SetDPx(0.);
76  TEnd->SetDPy(0.);
77  TEnd->SetDPz(0.);
78  return true;
79  }
80  LOG(warning) << "FairRKPropagator::Propagate not implemented yet or plane to propagate not set";
81  return false;
82 }
83 
84 bool FairRKPropagator::Propagate(FairTrackParP* TStart, FairTrackParH* TEnd, int PDG)
85 {
86  if (fPropagationFlag == TOPLANE) {
87  float x1[3] = {(float)TStart->GetX(), (float)TStart->GetY(), (float)TStart->GetZ()};
88  float p1[3] = {(float)TStart->GetPx(), (float)TStart->GetPy(), (float)TStart->GetPz()};
89  float x2[3] = {0., 0., 0.};
90  float p2[3] = {0., 0., 0.};
91  bool ret = Propagate(x1, p1, x2, p2, PDG);
92  if (!ret)
93  return false;
94  TEnd->SetX(x2[0]);
95  TEnd->SetY(x2[1]);
96  TEnd->SetZ(x2[2]);
97  TEnd->SetPx(p2[0]);
98  TEnd->SetPy(p2[1]);
99  TEnd->SetPz(p2[2]);
100  TEnd->SetQp(GetChargeFromPDG(PDG) / sqrt(p2[0] * p2[0] + p2[1] * p2[1] + p2[2] * p2[2]));
101  TEnd->SetDX(0.);
102  TEnd->SetDY(0.);
103  TEnd->SetDZ(0.);
104  TEnd->SetDPx(0.);
105  TEnd->SetDPy(0.);
106  TEnd->SetDPz(0.);
107  return true;
108  }
109  LOG(warning) << "FairRKPropagator::Propagate not implemented yet or plane to propagate not set";
110  return false;
111 }
112 
113 bool FairRKPropagator::Propagate(FairTrackParP* TStart, FairTrackParP* TEnd, int PDG)
114 {
115  if (fPropagationFlag == TOPLANE) {
116  float x1[3] = {(float)TStart->GetX(), (float)TStart->GetY(), (float)TStart->GetZ()};
117  float p1[3] = {(float)TStart->GetPx(), (float)TStart->GetPy(), (float)TStart->GetPz()};
118  float x2[3] = {0., 0., 0.};
119  float p2[3] = {0., 0., 0.};
120  bool ret = Propagate(x1, p1, x2, p2, PDG);
121  if (!ret)
122  return false;
123  TEnd->SetX(x2[0]);
124  TEnd->SetY(x2[1]);
125  TEnd->SetZ(x2[2]);
126  TEnd->SetPx(p2[0]);
127  TEnd->SetPy(p2[1]);
128  TEnd->SetPz(p2[2]);
129  TEnd->SetQp(GetChargeFromPDG(PDG) / sqrt(p2[0] * p2[0] + p2[1] * p2[1] + p2[2] * p2[2]));
130  TEnd->SetDX(0.);
131  TEnd->SetDY(0.);
132  TEnd->SetDZ(0.);
133  TEnd->SetDPx(0.);
134  TEnd->SetDPy(0.);
135  TEnd->SetDPz(0.);
136  return true;
137  }
138  LOG(warning) << "FairRKPropagator::Propagate not implemented yet or plane to propagate not set";
139  return false;
140 }
141 
142 bool FairRKPropagator::Propagate(FairTrackParH* TStart, FairTrackParP* TEnd, int PDG)
143 {
144  if (fPropagationFlag == TOPLANE) {
145  float x1[3] = {(float)TStart->GetX(), (float)TStart->GetY(), (float)TStart->GetZ()};
146  float p1[3] = {(float)TStart->GetPx(), (float)TStart->GetPy(), (float)TStart->GetPz()};
147  float x2[3] = {0., 0., 0.};
148  float p2[3] = {0., 0., 0.};
149  bool ret = Propagate(x1, p1, x2, p2, PDG);
150  if (!ret)
151  return false;
152  TEnd->SetX(x2[0]);
153  TEnd->SetY(x2[1]);
154  TEnd->SetZ(x2[2]);
155  TEnd->SetPx(p2[0]);
156  TEnd->SetPy(p2[1]);
157  TEnd->SetPz(p2[2]);
158  TEnd->SetQp(GetChargeFromPDG(PDG) / sqrt(p2[0] * p2[0] + p2[1] * p2[1] + p2[2] * p2[2]));
159  TEnd->SetDX(0.);
160  TEnd->SetDY(0.);
161  TEnd->SetDZ(0.);
162  TEnd->SetDPx(0.);
163  TEnd->SetDPy(0.);
164  TEnd->SetDPz(0.);
165  return true;
166  }
167  LOG(warning) << "FairRKPropagator::Propagate not implemented yet";
168  return false;
169 }
170 
171 bool FairRKPropagator::Propagate(float* x1, float* p1, float* x2, float* p2, int PDG)
172 {
173  double charge = GetChargeFromPDG(PDG);
174  double momIn = TMath::Sqrt(p1[0] * p1[0] + p1[1] * p1[1] + p1[2] * p1[2]);
175  if (momIn == 0.)
176  return false;
177  double vecIn[7] = {x1[0], x1[1], x1[2], p1[0] / momIn, p1[1] / momIn, p1[2] / momIn, momIn};
178  double vecOut[7];
179  double vec1[3] = {fDefPlaneV1.X(), fDefPlaneV1.Y(), fDefPlaneV1.Z()};
180  double vec2[3] = {fDefPlaneV2.X(), fDefPlaneV2.Y(), fDefPlaneV2.Z()};
181  double vec3[3] = {fDefPlaneV0.X(), fDefPlaneV0.Y(), fDefPlaneV0.Z()};
182  PropagateToPlane(charge, vecIn, vec1, vec2, vec3, vecOut);
183  x2[0] = vecOut[0];
184  x2[1] = vecOut[1];
185  x2[2] = vecOut[2];
186  p2[0] = vecOut[3] * vecOut[6];
187  p2[1] = vecOut[4] * vecOut[6];
188  p2[2] = vecOut[5] * vecOut[6];
189  return true;
190 }
191 
192 bool FairRKPropagator::SetDestinationPlane(const TVector3& v0, const TVector3& v1, const TVector3& v2)
193 {
194  fDefPlaneV0 = v0;
195  fDefPlaneV1 = v1;
196  fDefPlaneV2 = v2;
197  fPropagationFlag = TOPLANE;
198  return true;
199 }
200 
201 bool FairRKPropagator::SetOriginPlane([[gnu::unused]] const TVector3& v0, [[gnu::unused]] const TVector3& v1)
202 {
203  LOG(warning) << "FairRKPropagator::SetLengthToPropagateTo not implemented yet";
204  return false;
205 }
206 
207 bool FairRKPropagator::SetDestinationVolume([[gnu::unused]] std::string volName,
208  [[gnu::unused]] int copyNo,
209  [[gnu::unused]] int option)
210 {
211  LOG(warning) << "FairRKPropagator::SetDestinationVolume not implemented yet";
212  return false;
213 }
214 
215 bool FairRKPropagator::SetDestinationLength([[gnu::unused]] float length)
216 {
217  LOG(warning) << "FairRKPropagator::SetDestinationLength not implemented yet";
218  return false;
219 }
220 
221 //______________________________________________________________________________
222 bool FairRKPropagator::SetPCAPropagation(int pca, int dir, FairTrackParP* par)
223 {
224  if (abs(dir) != 1) {
225  LOG(warning) << "Set PCA dir to +1 or -1.";
226  return false;
227  }
228  fPCAPropagationType = pca;
229  fPCAPropagationDir = dir;
230  fPCAPropagationPar = par;
231  return true;
232 }
233 
234 //______________________________________________________________________________
235 PCAOutputStruct FairRKPropagator::FindPCA(int PCA,
236  int PDGCode,
237  TVector3 Point,
238  TVector3 Wire1,
239  TVector3 Wire2,
240  [[gnu::unused]] double MaxDistance)
241 {
242  PCAOutputStruct pcastruct = PCAOutputStruct();
243 
244  if (PCA != 1 && PCA != 2) {
245  LOG(info) << "FairRKPropagator::FindPCA implemented for point (pca=1) and wire (pca=2) only";
246  return pcastruct;
247  }
248  double charge = GetChargeFromPDG(PDGCode);
249  double momIn = fPCAPropagationDir * // if set to -1, it will back propagate
250  TMath::Sqrt(fPCAPropagationPar->GetPx() * fPCAPropagationPar->GetPx()
251  + fPCAPropagationPar->GetPy() * fPCAPropagationPar->GetPy()
252  + fPCAPropagationPar->GetPz() * fPCAPropagationPar->GetPz());
253  if (momIn == 0.)
254  return pcastruct;
255  double vecIn[7] = {fPCAPropagationPar->GetX(),
256  fPCAPropagationPar->GetY(),
257  fPCAPropagationPar->GetZ(),
258  fPCAPropagationPar->GetPx() / momIn,
259  fPCAPropagationPar->GetPy() / momIn,
260  fPCAPropagationPar->GetPz() / momIn,
261  momIn};
262 
263  double diff;
264  if (PCA == 1)
265  diff = sqrt((vecIn[0] - Point.X()) * (vecIn[0] - Point.X()) + (vecIn[1] - Point.Y()) * (vecIn[1] - Point.Y())
266  + (vecIn[2] - Point.Z()) * (vecIn[2] - Point.Z()));
267  else // if ( PCA == 2 )
268  diff = CalculatePointToWireDistance(
269  TVector3(fPCAPropagationPar->GetX(), fPCAPropagationPar->GetY(), fPCAPropagationPar->GetZ()),
270  Wire1,
271  Wire2,
272  pcastruct.OnWirePCA);
273 
274  fMaxStep = diff / 25;
275  double res_old = diff;
276  double res = 100.0;
277  double vecOut[7];
278  double vecOutT[7];
279  for (int i = 0; i < 7; i++) {
280  vecOut[i] = 0;
281  vecOutT[i] = 0;
282  }
283 
284  int nIter = 0;
285  pcastruct.TrackLength = 0.;
286 
287  do {
288  double stepLength = Step(charge, vecIn, vecOut);
289  double newDiff;
290  if (PCA == 1)
291  newDiff = sqrt((vecOut[0] - Point.X()) * (vecOut[0] - Point.X())
292  + (vecOut[1] - Point.Y()) * (vecOut[1] - Point.Y())
293  + (vecOut[2] - Point.Z()) * (vecOut[2] - Point.Z()));
294  else // if ( PCA == 2 )
295  newDiff = CalculatePointToWireDistance(
296  TVector3(vecOut[0], vecOut[1], vecOut[2]), Wire1, Wire2, pcastruct.OnWirePCA);
297 
298  res = newDiff / diff;
299  if (TMath::Abs(res) < 0.01 || res > res_old) {
300  break;
301  } else {
302  for (int i = 0; i < 7; i++) {
303  vecOutT[i] = vecOut[i];
304  vecIn[i] = vecOut[i];
305  }
306  res_old = res;
307  pcastruct.TrackLength += stepLength;
308  }
309  if (nIter++ > 1000) {
310  break;
311  }
312  } while (1);
313  if (res > res_old)
314  for (int k = 0; k < 7; k++) {
315  vecOut[k] = vecOutT[k];
316  }
317 
318  pcastruct.OnTrackPCA.SetX(vecOut[0]);
319  pcastruct.OnTrackPCA.SetY(vecOut[1]);
320  pcastruct.OnTrackPCA.SetZ(vecOut[2]);
321 
322  if (PCA == 1)
323  pcastruct.Distance =
324  sqrt((vecOut[0] - Point.X()) * (vecOut[0] - Point.X()) + (vecOut[1] - Point.Y()) * (vecOut[1] - Point.Y())
325  + (vecOut[2] - Point.Z()) * (vecOut[2] - Point.Z()));
326  else // if ( pcastruct.PCA == 2 )
327  pcastruct.Distance =
328  CalculatePointToWireDistance(TVector3(vecOut[0], vecOut[1], vecOut[2]), Wire1, Wire2, pcastruct.OnWirePCA);
329  pcastruct.PCAStatusFlag = 0;
330  return pcastruct;
331 }
332 //______________________________________________________________________________
333 
334 //______________________________________________________________________________
335 double FairRKPropagator::CalculatePointToWireDistance(TVector3 point, TVector3 wire1, TVector3 wire2, TVector3& vwi)
336 {
337  TVector3 ab = wire2 - wire1;
338  TVector3 av = point - wire1;
339 
340  if (av.Dot(ab) <= 0.0) { // Point is lagging behind start of the segment, so perpendicular distance is not viable
341  vwi = wire1;
342  return av.Mag(); // Use distance to start of segment instead
343  }
344 
345  TVector3 bv = point - wire2;
346 
347  if (bv.Dot(ab) >= 0.0) { // Point is advanced past the end of the segment, so perpendicular distance is not viable
348  vwi = wire2;
349  return bv.Mag(); // Use distance to end of the segment instead
350  }
351 
352  vwi = (ab.Dot(av) / ab.Mag() / ab.Mag()) * ab;
353  vwi += wire1;
354  return (ab.Cross(av)).Mag() / ab.Mag(); // Perpendicular distance of point to segment
355 }
356 //______________________________________________________________________________
357 
358 //______________________________________________________________________________
359 void FairRKPropagator::PropagateToPlane(double Charge,
360  double* vecRKIn,
361  double* vec1,
362  double* vec2,
363  double* vec3,
364  double* vecOut)
365 {
370  double Norm[3];
371  double Mag;
372  double dist[3];
373  double distance[3];
374  double vecRKoutT[7];
375 
376  for (int i = 0; i < 7; i++) {
377  vecRKoutT[i] = 0;
378  }
379 
380  Norm[0] = vec1[1] * vec2[2] - vec2[2] * vec2[1]; // a2b3 − a3b2,
381  Norm[1] = vec1[2] * vec2[0] - vec1[0] * vec2[2]; // a3b1 − a1b3;
382  Norm[2] = vec1[0] * vec2[1] - vec1[1] * vec2[0];
383 
384  Mag = TMath::Sqrt(Norm[0] * Norm[0] + Norm[1] * Norm[1] + Norm[2] * Norm[2]);
385 
386  // printf(" Mag = %f \n ", Mag);
387 
388  Norm[0] = Norm[0] / Mag;
389  Norm[1] = Norm[1] / Mag;
390  Norm[2] = Norm[2] / Mag;
391  // printf(" after normalization : Normal = %f %f %f \n ", Norm[0],Norm[1],Norm[2]);
392 
393  dist[0] = vecRKIn[0] - vec3[0];
394  dist[1] = vecRKIn[1] - vec3[1];
395  dist[2] = vecRKIn[2] - vec3[2];
396 
397  distance[0] = Norm[0] * dist[0];
398  distance[1] = Norm[1] * dist[1];
399  distance[2] = Norm[2] * dist[2];
400  // printf(" distance = %f %f %f \n ", distance[0],distance[1],distance[2]);
401  double diff = TMath::Abs(distance[0] + distance[1] + distance[2]);
402  fMaxStep = diff;
403  double res = 100.0;
404  double res_old = 100.0;
405 
406  double vecRKOut[7];
407  // for (int i=0; i< 7; i++) {vecRKOut[i]=0;}
408  int nIter = 0;
409 
410  // printf("I am in CPU code %f %f %f res= %f diff = %f \n ", vecRKIn[0], vecRKIn[1],vecRKIn[2], res, diff);
411 
412  do {
413  Step(Charge, vecRKIn, vecRKOut);
414  dist[0] = (vecRKOut[0] - vec3[0]) * Norm[0];
415  dist[1] = (vecRKOut[1] - vec3[1]) * Norm[1];
416  dist[2] = (vecRKOut[2] - vec3[2]) * Norm[2];
417  fMaxStep = TMath::Sqrt(dist[0] * dist[0] + dist[1] * dist[1] + dist[2] * dist[2]);
418  res = TMath::Abs(fMaxStep / diff);
419  // printf("After %i step %f %f %f res = %f \n", nIter ,vecRKOut[0], vecRKOut[1],vecRKOut[2] , res);
420  if (res < 0.001 || res > res_old) {
421  break;
422  } else {
423  for (int i = 0; i < 7; i++) {
424  vecRKIn[i] = vecRKOut[i];
425  vecRKoutT[i] = vecRKOut[i];
426  }
427  res_old = res;
428  }
429  if (nIter++ > 1000) {
430  break;
431  }
432  } while (1);
433  // printf("The results is %f %f %f , no of iter %i \n", vecRKOut[0],vecRKOut[1],vecRKOut[2], nIter);
434  // printf("\n");
435  for (int i = 0; i < 7; i++) {
436  if (res > res_old) {
437  vecOut[i] = vecRKoutT[i];
438  } else {
439  vecOut[i] = vecRKOut[i];
440  }
441  }
442 }
443 //______________________________________________________________________________
444 void FairRKPropagator::Propagate(double Charge, double* vecRKIn, double* Pos)
445 {
446  double diff = Pos[2] - vecRKIn[2];
447  fMaxStep = diff / 25;
448  double res_old = diff;
449  double res = 100.0;
450  double vecRKOut[7];
451  double vecRKOutT[7];
452  for (int i = 0; i < 7; i++) {
453  vecRKOut[i] = 0;
454  vecRKOutT[i] = 0;
455  }
456 
457  int nIter = 0;
458  do {
459  Step(Charge, vecRKIn, vecRKOut);
460  res = (vecRKOut[2] - Pos[2]) / diff;
461  if (TMath::Abs(res) < 0.01 || res > res_old) {
462  break;
463  } else {
464  for (int i = 0; i < 7; i++) {
465  vecRKOutT[i] = vecRKOut[i];
466  vecRKIn[i] = vecRKOut[i];
467  }
468  }
469  if (nIter++ > 1000) {
470  break;
471  }
472  } while (1);
473  if (res > res_old)
474  for (int k = 0; k < 7; k++) {
475  vecRKOut[k] = vecRKOutT[k];
476  }
477  for (int k = 0; k < 3; k++) {
478  printf(" vecRKOut[%i] =%f ", k, vecRKOut[k]);
479  }
480  printf("\n");
481 }
482 //______________________________________________________________________________
483 double FairRKPropagator::Step(double Charge, double* vecRKIn, double* vecOut)
484 {
485 
486  double vecRKOut[7];
487  for (int i = 0; i < 7; i++) {
488  vecRKOut[i] = 0;
489  }
490  // for (int i=0; i< 7; i++) printf( "vectRKIn(%i)=%f \n",i ,vecRKIn[i]);
491  // printf(" ---------------------------------------------------------------- \n");
492  double stepLength = OneStepRungeKutta(Charge, fMaxStep, vecRKIn, vecRKOut);
493  // printf(" now at x=%f y=%f z=%f \n", vecRKOut[0],vecRKOut[1],vecRKOut[2]);
494  vecOut[0] = vecRKOut[0];
495  vecOut[1] = vecRKOut[1];
496  vecOut[2] = vecRKOut[2];
497  vecOut[6] = vecRKOut[6];
498  vecOut[3] = vecRKOut[3];
499  vecOut[4] = vecRKOut[4];
500  vecOut[5] = vecRKOut[5];
501 
502  return stepLength;
503 }
504 
505 //______________________________________________________________________________
506 double FairRKPropagator::OneStepRungeKutta(double charge, double step, double* vect, double* vout)
507 {
508 
509  // Wrapper to step with method RungeKutta.
510 
532 
533  double f[4];
534  double xyzt[3];
535  double secxs[4], secys[4], seczs[4]; // hxp[3];
536  // double /*g1 , g2, g3, g4, g5, g6,*/;
537  // double /*f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost*/;
538 
539  double track_length = 0.;
540 
541  double maxit = 1992;
542  double maxcut = 11;
543 
544  const double hmin = 1e-4; // !!! MT ADD, should be member
545  const double kdlt = 1e-6; // !!! MT CHANGE from 1e-4, should be member
546  const double kdlt32 = kdlt / 32.;
547  const double kthird = 1. / 3.;
548  const double khalf = 0.5;
549  const double kec = 2.99792458e-4;
550  const double kpisqua = 9.86960440109;
551  /* const int kix = 0;
552  const int kiy = 1;
553  const int kiz = 2;
554  const int kipx = 3;
555  const int kipy = 4;
556  const int kipz = 5;
557  */
558  // *.
559  // *. ------------------------------------------------------------------
560  // *.
561  // * this constant is for units cm,gev/c and kgauss
562  // *
563  int iter = 0;
564  int ncut = 0;
565  for (int j = 0; j < 7; j++) {
566  vout[j] = vect[j];
567  }
568 
569  double pinv = kec * charge / vect[6];
570  double tl = 0.;
571  double h = step;
572 
573  do {
574  double rest = step - tl;
575  // 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);
576  if (TMath::Abs(h) > TMath::Abs(rest)) {
577  h = rest;
578  }
579 
580  fMagField->GetFieldValue(vout, f);
581 
582  // * start of integration
583  double x = vout[0];
584  double y = vout[1];
585  double z = vout[2];
586  double a = vout[3];
587  double b = vout[4];
588  double c = vout[5];
589 
590  double h2 = khalf * h;
591  double h4 = khalf * h2;
592  double ph = pinv * h;
593  double ph2 = khalf * ph;
594 
595  // printf(" ------------------------------------------- h2 = %f\n",h2);
596 
597  secxs[0] = (b * f[2] - c * f[1]) * ph2;
598  secys[0] = (c * f[0] - a * f[2]) * ph2;
599  seczs[0] = (a * f[1] - b * f[0]) * ph2;
600  double ang2 = (secxs[0] * secxs[0] + secys[0] * secys[0] + seczs[0] * seczs[0]);
601  if (ang2 > kpisqua) {
602  break;
603  }
604 
605  double dxt = h2 * a + h4 * secxs[0];
606  double dyt = h2 * b + h4 * secys[0];
607  double dzt = h2 * c + h4 * seczs[0];
608  double xt = x + dxt;
609  double yt = y + dyt;
610  double zt = z + dzt;
611  // printf(" Position 1 at xt=%f yt=%f zt=%f \n", xt, yt, zt);
612  // printf(" differance dxt=%f dyt=%f dzt=%f \n", dxt, dyt, dzt);
613  // * second intermediate point
614  double est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
615  if (est > h) {
616  if (ncut++ > maxcut) {
617  break;
618  }
619  h *= khalf;
620  continue;
621  }
622 
623  xyzt[0] = xt;
624  xyzt[1] = yt;
625  xyzt[2] = zt;
626 
627  fMagField->GetFieldValue(xyzt, f);
628 
629  double at = a + secxs[0];
630  double bt = b + secys[0];
631  double ct = c + seczs[0];
632 
633  secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
634  secys[1] = (ct * f[0] - at * f[2]) * ph2;
635  seczs[1] = (at * f[1] - bt * f[0]) * ph2;
636  at = a + secxs[1];
637  bt = b + secys[1];
638  ct = c + seczs[1];
639  secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
640  secys[2] = (ct * f[0] - at * f[2]) * ph2;
641  seczs[2] = (at * f[1] - bt * f[0]) * ph2;
642  dxt = h * (a + secxs[2]);
643  dyt = h * (b + secys[2]);
644  dzt = h * (c + seczs[2]);
645  xt = x + dxt;
646  yt = y + dyt;
647  zt = z + dzt;
648  at = a + 2. * secxs[2];
649  bt = b + 2. * secys[2];
650  ct = c + 2. * seczs[2];
651  // printf(" Position 2 at xt=%f yt=%f zt=%f \n", xt, yt, zt);
652 
653  est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
654  if (est > 2. * TMath::Abs(h)) {
655  if (ncut++ > maxcut) {
656  break;
657  }
658  h *= khalf;
659  continue;
660  }
661 
662  xyzt[0] = xt;
663  xyzt[1] = yt;
664  xyzt[2] = zt;
665 
666  fMagField->GetFieldValue(xyzt, f);
667 
668  z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
669  y = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
670  x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
671  // printf(" Position 3 at x=%f y=%f z=%f \n", x, y, z);
672  secxs[3] = (bt * f[2] - ct * f[1]) * ph2;
673  secys[3] = (ct * f[0] - at * f[2]) * ph2;
674  seczs[3] = (at * f[1] - bt * f[0]) * ph2;
675  a = a + (secxs[0] + secxs[3] + 2. * (secxs[1] + secxs[2])) * kthird;
676  b = b + (secys[0] + secys[3] + 2. * (secys[1] + secys[2])) * kthird;
677  c = c + (seczs[0] + seczs[3] + 2. * (seczs[1] + seczs[2])) * kthird;
678 
679  est = TMath::Abs(secxs[0] + secxs[3] - (secxs[1] + secxs[2]))
680  + TMath::Abs(secys[0] + secys[3] - (secys[1] + secys[2]))
681  + TMath::Abs(seczs[0] + seczs[3] - (seczs[1] + seczs[2]));
682 
683  if (est > kdlt && TMath::Abs(h) > hmin) {
684  if (ncut++ > maxcut) {
685  break;
686  }
687  h *= khalf;
688  continue;
689  }
690 
691  ncut = 0;
692  // * if too many iterations, go to helix
693  if (iter++ > maxit) {
694  break;
695  }
696 
697  tl += h;
698  if (est < kdlt32) {
699  h *= 2.;
700  }
701  double cba = 1. / TMath::Sqrt(a * a + b * b + c * c);
702  vout[0] = x;
703  vout[1] = y;
704  vout[2] = z;
705  vout[3] = cba * a;
706  vout[4] = cba * b;
707  vout[5] = cba * c;
708 
709  rest = step - tl;
710  track_length = tl;
711  // printf(" Position 4 at x=%f y=%f z=%f Step = %f \n", x, y, z, step );
712 
713  if (step < 0.) {
714  rest = -rest;
715  }
716  if (rest < 1.e-5 * TMath::Abs(step)) {
717  return track_length;
718  }
719 
720  } while (1);
721 
722  return track_length;
723 
724  // angle too big, use helix
725  /*
726  f1 = f[0];
727  f2 = f[1];
728  f3 = f[2];
729  f4 = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
730  rho = -f4*pinv;
731  tet = rho * step;
732 
733  hnorm = 1./f4;
734  f1 = f1*hnorm;
735  f2 = f2*hnorm;
736  f3 = f3*hnorm;
737 
738  hxp[0] = f2*vect[kipz] - f3*vect[kipy];
739  hxp[1] = f3*vect[kipx] - f1*vect[kipz];
740  hxp[2] = f1*vect[kipy] - f2*vect[kipx];
741 
742  hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];
743 
744  rho1 = 1./rho;
745  sint = TMath::Sin(tet);
746  cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);
747 
748  g1 = sint*rho1;
749  g2 = cost*rho1;
750  g3 = (tet-sint) * hp*rho1;
751  g4 = -cost;
752  g5 = sint;
753  g6 = cost * hp;
754 
755  vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
756  vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
757  vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
758 
759  vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
760  vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
761  vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
762  */
763 }
FairTrackPar::GetPz
virtual Double_t GetPz() const
Definition: FairTrackPar.h:62
FairRKPropagator::SetDestinationPlane
virtual bool SetDestinationPlane(const TVector3 &v0, const TVector3 &v1, const TVector3 &v2)
Definition: FairRKPropagator.cxx:192
FairRKPropagator::FindPCA
virtual PCAOutputStruct FindPCA(int PCA, int PDGCode, TVector3 Point, TVector3 Wire1, TVector3 Wire2, double MaxDistance)
Definition: FairRKPropagator.cxx:235
FairTrackPar::SetDPy
virtual void SetDPy(Double_t dpy)
Definition: FairTrackPar.h:85
FairRKPropagator::Step
double Step(double Charge, double *vecRKIn, double *vecOut)
Definition: FairRKPropagator.cxx:483
FairTrackPar::SetDPz
virtual void SetDPz(Double_t dpz)
Definition: FairTrackPar.h:86
FairRKPropagator::SetDestinationVolume
virtual bool SetDestinationVolume(std::string volName, int copyNo, int option)
Definition: FairRKPropagator.cxx:207
FairTrackPar::SetPy
virtual void SetPy(Double_t py)
Definition: FairTrackPar.h:81
FairTrackParH::GetY
Double_t GetY()
Definition: FairTrackParH.cxx:668
ClassImp
ClassImp(FairEventBuilder)
FairRKPropagator::Propagate
void Propagate(double Charge, double *vecRKIn, double *Pos)
Definition: FairRKPropagator.cxx:444
FairTrackPar::SetY
virtual void SetY(Double_t y)
Definition: FairTrackPar.h:73
FairTrackPar::GetZ
virtual Double_t GetZ()
Definition: FairTrackPar.h:51
FairRKPropagator::SetPCAPropagation
virtual bool SetPCAPropagation(int pca, int dir=1, FairTrackParP *par=nullptr)
Definition: FairRKPropagator.cxx:222
FairTrackParH::GetX
Double_t GetX()
Definition: FairTrackParH.cxx:660
FairTrackPar::SetZ
virtual void SetZ(Double_t z)
Definition: FairTrackPar.h:74
FairTrackPar::GetPx
virtual Double_t GetPx() const
Definition: FairTrackPar.h:60
FairTrackParH::GetZ
Double_t GetZ()
Definition: FairTrackParH.cxx:676
PCAOutputStruct::OnTrackPCA
TVector3 OnTrackPCA
Definition: FairPropagator.h:50
FairTrackPar::SetDPx
virtual void SetDPx(Double_t dpx)
Definition: FairTrackPar.h:84
PCAOutputStruct::OnWirePCA
TVector3 OnWirePCA
Definition: FairPropagator.h:51
FairRKPropagator::OneStepRungeKutta
double OneStepRungeKutta(double charge, double step, double *vect, double *vout)
Definition: FairRKPropagator.cxx:506
par
FairMQExParamsParOne * par
Definition: fill_parameters.C:11
FairTrackPar::SetPz
virtual void SetPz(Double_t pz)
Definition: FairTrackPar.h:82
FairTrackPar::SetPx
virtual void SetPx(Double_t px)
Definition: FairTrackPar.h:80
FairField::GetFieldValue
virtual void GetFieldValue(const Double_t point[3], Double_t *bField)
Definition: FairField.cxx:44
FairTrackPar::GetPy
virtual Double_t GetPy() const
Definition: FairTrackPar.h:61
FairTrackPar::SetX
virtual void SetX(Double_t x)
Definition: FairTrackPar.h:72
FairTrackPar::GetY
virtual Double_t GetY()
Definition: FairTrackPar.h:50
FairTrackPar::SetDY
virtual void SetDY(Double_t dy)
Definition: FairTrackPar.h:77
FairRKPropagator::SetDestinationLength
virtual bool SetDestinationLength(float length)
Definition: FairRKPropagator.cxx:215
FairTrackPar::GetX
virtual Double_t GetX()
Definition: FairTrackPar.h:49
FairTrackPar::SetQp
void SetQp(Double_t qp)
Definition: FairTrackPar.h:88
FairTrackPar::SetDX
virtual void SetDX(Double_t dx)
Definition: FairTrackPar.h:76
FairRKPropagator::SetOriginPlane
virtual bool SetOriginPlane(const TVector3 &v0, const TVector3 &v1)
Definition: FairRKPropagator.cxx:201
FairTrackPar::SetDZ
virtual void SetDZ(Double_t dz)
Definition: FairTrackPar.h:78
FairRKPropagator::~FairRKPropagator
virtual ~FairRKPropagator()
Definition: FairRKPropagator.cxx:44
FairRKPropagator::PropagateToPlane
void PropagateToPlane(double Charge, double *vecRKIn, double *vec1, double *vec2, double *vec3, double *vecOut)
Definition: FairRKPropagator.cxx:359