ATLAS Offline Software
Loading...
Searching...
No Matches
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
12namespace Trk {
13
14void 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
69L20:
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 */
185L30:
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 */
195L40:
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
charge
double charge(const T &p)
Definition AtlasPID.h:997
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
xt
#define xt
yt
#define yt
zt
#define zt
xyzt
#define xyzt
Trk::vkalMagFld::getMagFld
static void getMagFld(const double, const double, const double, double &, double &, double &, const VKalVrtControlBase *)
Definition VKalVrtBMag.cxx:24
cost
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)
Definition hcg.cxx:922
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition FakeTrackBuilder.h:9
Trk::x
@ x
Definition ParamDefs.h:55
Trk::z
@ z
global position (cartesian)
Definition ParamDefs.h:57
Trk::y
@ y
Definition ParamDefs.h:56
Trk::vkgrkuta_
void vkgrkuta_(const double charge, const double step, double *vect, double *vout, VKalVrtControlBase *CONTROL)
Definition VKgrkuta.cxx:14
Generated on for ATLAS Offline Software by doxygen 1.14.0