ATLAS Offline Software
VKgrkuta.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "TrkVKalVrtCore/VKgrkuta.h"
6 #include "TrkVKalVrtCore/CommonPars.h"
7 #include "TrkVKalVrtCore/TrkVKalVrtCore.h"
8 #include "TrkVKalVrtCore/VKalVrtBMag.h"
9 #include <cmath>
10 #include <iostream>
11 
12 namespace Trk {
13 
14 void vkgrkuta_(const double charge, const double step, double *vect, double *vout, VKalVrtControlBase* CONTROL)
15 {
16  double equiv_2[3], equiv_5[3];
17  long int iter, ncut, j;
18  double cost, pinv, rest, sint, a, b, c, f[4], hst;
19  double hnorm, secxs[4], secys[4], seczs[4], f1, f2, h2, f3, h4, f4;
20  double g1, g2, g3, g4, g5, g6, at, bt, ct, ph, hp, fx, fy, fz, tl;
21  double ph2, cba, rho, est, tet, hxp[3], dxt, dyt, dzt, ang2, rho1;
22 
23 
24 
25 #define xyz (equiv_2)
26 #define x (equiv_2)
27 #define y (equiv_2 + 1)
28 #define z (equiv_2 + 2)
29 #define xyzt (equiv_5)
30 #define xt (equiv_5)
31 #define yt (equiv_5 + 1)
32 #define zt (equiv_5 + 2)
33 
34 
35 /* *******************************************************************/
36 /* */
37 /* Runge-Kutta method for tracking a particle through a magnetic */
38 /* . field. Uses Nystroem algorithm (See Handbook Nat. Bur. of */
39 /* . Standards, procedure 25.5.20) */
40 /* . */
41 /* . Input parameters */
42 /* . CHARGE Particle charge */
43 /* . STEP Step size */
44 /* . VECT Initial co-ords,direction cosines,momentum */
45 /* . Output parameters */
46 /* . VOUT Output co-ords,direction cosines,momentum */
47 /* . User routine called */
48 /* CALL GUFLD(X,F) */
49 /* Authors R.Brun, M.Hansroul ********* */
50 /* V.Perevoztchikov (CUT STEP implementation) */
51 /* Taken from GEANT3.21 */
52 /* Modified for VKalVrt V.Kostyukhin */
53 /*********************************************************************/
54 /* This constant is for units CM,GEV/C and KGAUSS */
55 
56  /* Parameter adjustments */
57  --vout;
58  --vect;
59 
60  iter = 0;
61  ncut = 0;
62  for (j = 1; j <= 7; ++j) {vout[j] = vect[j];}
63  pinv = (charge) * 2.9979251e-2 / vect[7]; // New for MM, MEV/C and KGAUSS
64  tl = 0.;
65  hst = step;
66 
67 //std::cout <<" Now in grkuta="<<vect[1]<<", "<<vect[2]<<", "<<vect[3]<<'\n';
68 
69 L20:
70  rest = step - tl;
71  if (std::abs(hst) > std::abs(rest)) hst = rest;
72 /* ***** CALL GUFLD(VOUT,F) */
73  Trk::vkalMagFld::getMagFld( vout[1], vout[2], vout[3], fx, fy, fz, CONTROL);
74  f[0] = fx*10.; //VK Convert returned field in Tesla into kGauss for old code
75  f[1] = fy*10.;
76  f[2] = fz*10.;
77 
78 /* Start of integration */
79 
80  *x = vout[1];
81  *y = vout[2];
82  *z = vout[3];
83  a = vout[4];
84  b = vout[5];
85  c = vout[6];
86 
87  h2 = hst * .5;
88  h4 = h2 * .5;
89  ph = pinv * hst;
90  ph2 = ph * .5;
91  secxs[0] = (b * f[2] - c * f[1]) * ph2;
92  secys[0] = (c * f[0] - a * f[2]) * ph2;
93  seczs[0] = (a * f[1] - b * f[0]) * ph2;
94  ang2 = secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0];
95  if (ang2 > 9.86960440109) goto L40;
96  dxt = h2 * a + h4 * secxs[0];
97  dyt = h2 * b + h4 * secys[0];
98  dzt = h2 * c + h4 * seczs[0];
99  *xt = *x + dxt;
100  *yt = *y + dyt;
101  *zt = *z + dzt;
102 
103 /* Second intermediate point */
104 
105  est = std::abs(dxt) + std::abs(dyt) + std::abs(dzt);
106  if (est > hst) goto L30;
107 
108 
109 /* ***** CALL GUFLD(XYZT,F) */
110  Trk::vkalMagFld::getMagFld( xyzt[0], xyzt[1], xyzt[2], fx, fy, fz, CONTROL);
111  f[0] = fx*10.; //VK Convert returned field in Tesla into kGauss for old code
112  f[1] = fy*10.;
113  f[2] = fz*10.;
114  at = a + secxs[0];
115  bt = b + secys[0];
116  ct = c + seczs[0];
117 
118  secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
119  secys[1] = (ct * f[0] - at * f[2]) * ph2;
120  seczs[1] = (at * f[1] - bt * f[0]) * ph2;
121  at = a + secxs[1];
122  bt = b + secys[1];
123  ct = c + seczs[1];
124  secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
125  secys[2] = (ct * f[0] - at * f[2]) * ph2;
126  seczs[2] = (at * f[1] - bt * f[0]) * ph2;
127  dxt = hst * (a + secxs[2]);
128  dyt = hst * (b + secys[2]);
129  dzt = hst * (c + seczs[2]);
130  *xt = *x + dxt;
131  *yt = *y + dyt;
132  *zt = *z + dzt;
133  at = a + secxs[2] * 2.;
134  bt = b + secys[2] * 2.;
135  ct = c + seczs[2] * 2.;
136 
137  est = std::abs(dxt) + std::abs(dyt) + std::abs(dzt);
138  if (est > std::abs(hst) * 2.) goto L30;
139 /* ***** CALL GUFLD(XYZT,F) */
140  Trk::vkalMagFld::getMagFld( xyzt[0], xyzt[1], xyzt[2], fx, fy, fz, CONTROL);
141  f[0] = fx*10.; //VK Convert returned field in Tesla into kGauss for old code
142  f[1] = fy*10.;
143  f[2] = fz*10.;
144 
145  *z += (c + (seczs[0] + seczs[1] + seczs[2]) /3.) * hst;
146  *y += (b + (secys[0] + secys[1] + secys[2]) /3.) * hst;
147  *x += (a + (secxs[0] + secxs[1] + secxs[2]) /3.) * hst;
148 
149  secxs[3] = (bt * f[2] - ct * f[1]) * ph2;
150  secys[3] = (ct * f[0] - at * f[2]) * ph2;
151  seczs[3] = (at * f[1] - bt * f[0]) * ph2;
152  a += (secxs[0] + secxs[3] + (secxs[1] + secxs[2]) * 2.) /3.;
153  b += (secys[0] + secys[3] + (secys[1] + secys[2]) * 2.) /3.;
154  c += (seczs[0] + seczs[3] + (seczs[1] + seczs[2]) * 2.) /3.;
155 
156  est = std::abs(secxs[0] + secxs[3] - (secxs[1] + secxs[2]))
157  + std::abs(secys[0] + secys[3] - (secys[1] + secys[2]))
158  + std::abs(seczs[0] + seczs[3] - (seczs[1] + seczs[2]));
159 
160  if (est > 1e-4 && std::abs(hst) > 1e-4) goto L30;
161  ++iter;
162  ncut = 0;
163 
164 
165 /* If too many iterations, go to HELIX */
166  if (iter > 1992) goto L40;
167 
168  tl += hst;
169  if (est < 3.125e-6)
170  hst *= 2.;
171  cba = 1. / sqrt(a * a + b * b + c * c);
172  vout[1] = *x;
173  vout[2] = *y;
174  vout[3] = *z;
175  vout[4] = cba * a;
176  vout[5] = cba * b;
177  vout[6] = cba * c;
178  rest = step - tl;
179  if (step < 0.) rest = -rest;
180  if (rest > std::abs(step) * 1e-5) goto L20;
181 
182  return;
183 
184 /* * CUT STEP */
185 L30:
186  ++ncut;
187 /* If too many cuts , go to HELIX */
188  if (ncut > 11) {
189  goto L40;
190  }
191  hst *= .5;
192  goto L20;
193 
194 /* * ANGLE TOO BIG, USE HELIX */
195 L40:
196  f1 = f[0];
197  f2 = f[1];
198  f3 = f[2];
199  f4 = sqrt(f1*f1 + f2*f2 + f3*f3);
200  rho = -f4 * pinv;
201  tet = rho * step;
202  if (tet != 0.) {
203  hnorm = 1. / f4;
204  f1 *= hnorm;
205  f2 *= hnorm;
206  f3 *= hnorm;
207 
208  hxp[0] = f2 * vect[6] - f3 * vect[5];
209  hxp[1] = f3 * vect[4] - f1 * vect[6];
210  hxp[2] = f1 * vect[5] - f2 * vect[4];
211  hp = f1 * vect[4] + f2 * vect[5] + f3 * vect[6];
212 
213  rho1 = 1. / rho;
214  sint = sin(tet);
215  cost = sin(tet/2.) * sin(tet/2.) * 2.;
216 
217  g1 = sint * rho1;
218  g2 = cost * rho1;
219  g3 = (tet - sint) * hp * rho1;
220  g4 = -cost;
221  g5 = sint;
222  g6 = cost * hp;
223  vout[1] = vect[1] + (g1 * vect[4] + g2 * hxp[0] + g3 * f1);
224  vout[2] = vect[2] + (g1 * vect[5] + g2 * hxp[1] + g3 * f2);
225  vout[3] = vect[3] + (g1 * vect[6] + g2 * hxp[2] + g3 * f3);
226  vout[4] = vect[4] + (g4 * vect[4] + g5 * hxp[0] + g6 * f1);
227  vout[5] = vect[5] + (g4 * vect[5] + g5 * hxp[1] + g6 * f2);
228  vout[6] = vect[6] + (g4 * vect[6] + g5 * hxp[2] + g6 * f3);
229 
230  } else {
231  vout[1] = vect[1] + step * vect[4];
232  vout[2] = vect[2] + step * vect[5];
233  vout[3] = vect[3] + step * vect[6];
234 
235  }
236 
237  }
238 #undef xyz
239 #undef zt
240 #undef yt
241 #undef xt
242 #undef z
243 #undef y
244 #undef x
245 #undef xyzt
246 
247 } /* End of namespace */
248 
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
Trk::y
@ y
Definition: ParamDefs.h:56
CommonPars.h
cost
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)
Definition: hcg.cxx:921
yt
#define yt
Trk::z
@ z
global position (cartesian)
Definition: ParamDefs.h:57
Trk::vkalMagFld::getMagFld
static void getMagFld(const double, const double, const double, double &, double &, double &, VKalVrtControlBase *)
Definition: VKalVrtBMag.cxx:24
VKalVrtBMag.h
read_hist_ntuple.f2
f2
Definition: read_hist_ntuple.py:20
xAOD::EgammaParameters::f3
@ f3
fraction of energy reconstructed in 3rd sampling
Definition: EgammaEnums.h:54
xt
#define xt
fitman.g1
g1
Definition: fitman.py:619
zt
#define zt
TrkVKalVrtCore.h
mc.g4
g4
Definition: mc.PhPy8EG_A14NNPDF23_DY_VLQ_example.py:25
hist_file_dump.f
f
Definition: hist_file_dump.py:135
calibdata.ct
ct
Definition: calibdata.py:418
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
fitman.g2
g2
Definition: fitman.py:624
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
charge
double charge(const T &p)
Definition: AtlasPID.h:756
VKgrkuta.h
xyzt
#define xyzt
a
TList * a
Definition: liststreamerinfos.cxx:10
LArCellBinning.step
step
Definition: LArCellBinning.py:158
dqt_zlumi_alleff_HIST.tl
tl
Definition: dqt_zlumi_alleff_HIST.py:73
Trk::VKalVrtControlBase
Definition: TrkVKalVrtCore.h:26
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
Trk::x
@ x
Definition: ParamDefs.h:55
python.compressB64.c
def c
Definition: compressB64.py:93
Trk::vkgrkuta_
void vkgrkuta_(const double charge, const double step, double *vect, double *vout, VKalVrtControlBase *CONTROL)
Definition: VKgrkuta.cxx:14
fitman.rho
rho
Definition: fitman.py:532
read_hist_ntuple.f1
f1
Definition: read_hist_ntuple.py:4
Generated on Tue Dec 24 2024 21:22:36 for ATLAS Offline Software by   doxygen 1.8.18