Fork me on GitHub

source: svn/trunk/src/BFieldProp.cc@ 252

Last change on this file since 252 was 248, checked in by Xavier Rouby, 16 years ago

new TrackPropagation::bfield function
File size: 16.8 KB
RevLine 
[53]1/*
2 * ---- Delphes ----
3 * A Fast Simulator for general purpose LHC detector
4 * S. Ovyn ~~~~ severine.ovyn@uclouvain.be
5 *
6 * Center for Particle Physics and Phenomenology (CP3)
7 * Universite Catholique de Louvain (UCL)
8 * Louvain-la-Neuve, Belgium
9 * */
10
[219]11#include "BFieldProp.h"
[53]12#include<cmath>
13using namespace std;
14
15
16//------------------------------------------------------------------------------
17
[219]18TrackPropagation::TrackPropagation(){
19 DET = new RESOLution();
20 init();
21}
[53]22
[219]23TrackPropagation::TrackPropagation(const string& DetDatacard){
24 DET = new RESOLution();
25 DET->ReadDataCard(DetDatacard);
26 init();
27}
[53]28
[219]29TrackPropagation::TrackPropagation(const RESOLution* DetDatacard){
30 DET= new RESOLution(*DetDatacard);
31 init();
32}
33
34TrackPropagation::TrackPropagation(const TrackPropagation & tp){
35 MAXITERATION = tp.MAXITERATION;
36 DET = new RESOLution(*(tp.DET));
37 R_max = tp.R_max; z_max = tp.z_max;
38 B_x = tp.B_x; B_y = tp.B_y; B_z = tp.B_z;
39 q = tp.q; phi_0 = tp.phi_0;
40 gammam= tp.gammam; omega = tp.omega;
41 r = tp.r; rr = tp.rr;
42 x_c = tp.x_c; y_c = tp.y_c;
43 R_c = tp.R_c; Phi_c = tp.Phi_c;
44 t = tp.t; t_z = tp.t_z; t_T = tp.t_T;
45 x_t = tp.x_t; y_t = tp.y_t; z_t = tp.z_t;
46 R_t = tp.R_t; Phi_t = tp.Phi_t;
47 Theta_t=tp.Theta_t; Eta_t = tp.Eta_t;
48 Px_t = tp.Px_t; Py_t = tp.Py_t; Pz_t = tp.Pz_t;
49 PT_t = tp.PT_t; p_t = tp.p_t; E_t = tp.E_t;
50 loop_overflow_counter = tp.loop_overflow_counter;
51}
52
53TrackPropagation& TrackPropagation::operator=(const TrackPropagation & tp) {
54 if(this==&tp) return *this;
55 MAXITERATION = tp.MAXITERATION;
56 DET = new RESOLution(*(tp.DET));
57 R_max = tp.R_max; z_max = tp.z_max;
58 B_x = tp.B_x; B_y = tp.B_y; B_z = tp.B_z;
59 q = tp.q; phi_0 = tp.phi_0;
60 gammam= tp.gammam; omega = tp.omega;
61 r = tp.r; rr = tp.rr;
62 x_c = tp.x_c; y_c = tp.y_c;
63 R_c = tp.R_c; Phi_c = tp.Phi_c;
64 t = tp.t; t_z = tp.t_z; t_T = tp.t_T;
65 x_t = tp.x_t; y_t = tp.y_t; z_t = tp.z_t;
66 R_t = tp.R_t; Phi_t = tp.Phi_t;
67 Theta_t=tp.Theta_t; Eta_t = tp.Eta_t;
68 Px_t = tp.Px_t; Py_t = tp.Py_t; Pz_t = tp.Pz_t;
69 PT_t = tp.PT_t; p_t = tp.p_t; E_t = tp.E_t;
70 loop_overflow_counter = tp.loop_overflow_counter;
71 return *this;
72}
73
74void TrackPropagation::init() {
75 MAXITERATION = 10000;
76 q= UNDEFINED; phi_0= UNDEFINED; gammam= UNDEFINED; omega=UNDEFINED; r=UNDEFINED;
77 x_c=UNDEFINED; y_c=UNDEFINED; R_c=UNDEFINED; Phi_c=UNDEFINED;
78 rr=UNDEFINED; t=UNDEFINED; t_z=UNDEFINED; t_T=UNDEFINED;
79 x_t=UNDEFINED; y_t=UNDEFINED; z_t=UNDEFINED;
80 R_t=UNDEFINED; Phi_t=UNDEFINED; Theta_t=UNDEFINED; Eta_t=UNDEFINED;
81 Px_t=UNDEFINED; Py_t=UNDEFINED; Pz_t=UNDEFINED; PT_t=UNDEFINED; p_t=UNDEFINED; E_t=UNDEFINED;
82
83 // DET has been initialised in the constructors
84 // magnetic field parameters
[193]85 R_max = DET->TRACK_radius;
86 z_max = DET->TRACK_length/2.;
87 B_x = DET->TRACK_bfield_x;
88 B_y = DET->TRACK_bfield_y;
89 B_z = DET->TRACK_bfield_z;
90
91 loop_overflow_counter=0;
[53]92}
93
[219]94
95
[193]96void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &momentum) {
[53]97
[193]98 q = Charge(Part->PID);
99 if(q==0) return;
100
101 if(R_max ==0) { cout << "ERROR: magnetic field has no lateral extention\n"; return;}
102 if(z_max==0) { cout << "ERROR: magnetic field has no longitudinal extention\n"; return;}
103
[199]104 if (B_x== 0 && B_y== 0) { // faster if only B_z
[193]105 if (B_z==0) return; // nothing to do
106
107 // initial conditions:
108 // p_X0 = Part->Px, p_Y0 = Part->Py, p_Z0 = Part->Pz, p_T0 = Part->PT;
109 // X_0 = Part->X, Y_0 = Part->Y, Z_0 = Part->Z;
110
111 // 1. initial transverse momentum p_{T0} : Part->PT
112 // initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0)
113 // relativistic gamma : gamma = E/mc² ; gammam = gamma \times m
114 // giration frequency \omega = q/(gamma m) B_z
115 // helix radius r = p_T0 / (omega gamma m)
116 phi_0 = -atan2(Part->Px,Part->Py);
117 gammam = Part->E; // here c==1
118 //cout << "gammam" << gammam << "\t gamma" << gammam/Part->M << endl;
119 omega = q * B_z /gammam;
120 r = Part->PT / (omega * gammam);
121
122 // 2. Helix parameters : center coordinates in transverse plane
123 // x_c = x_0 - r*cos(phi_0) and y_c = y_0 - r*sin(phi_0)
124 // R_c = \sqrt{x_c² + y_c²} and \Phi_c = atan{y_c/x_c}
[199]125 x_c = Part->X - r*cos(phi_0); /// TEST !!
[193]126 y_c = Part->Y - r*sin(phi_0);
127 R_c = sqrt(pow(x_c,2.) + pow(y_c,2.) );
128 Phi_c = atan2(y_c,x_c);
129
130 // 3. time evaluation t = min(t_T, t_z)
131 // t_T : time to exit from the sides
[199]132 // t_T= [ Phi_c - phi_0 + atan( (R_max^2 - (R_c^2 + r^2))/(2rR_c) ) ]/omega
[193]133 // t_z : time to exit from the front or the back
[199]134 // t_z = gamma * m /p_z0 \times (-z_0 + z_max * sign(p_z0))
[193]135 rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable
136 t_T=0;
[219]137 int sign_pz= (Part->Pz >0) ? 1 : -1;
138 t_z = gammam / Part->Pz * (-Part->Z + z_max*sign_pz ) ;
[193]139 if ( fabs(R_c - r) > R_max || R_c + r < R_max ) t = t_z;
140 else {
141 t_T = (Phi_c - phi_0 + atan2( (R_max + rr)*(R_max - rr) , 2*r*R_c ) ) / omega;
142 t = min(t_T,t_z);
143 }
144
145 // 4. position in terms of x(t), y(t), z(t)
146 // x(t) = x_c + r cos (omega t + phi_0)
147 // y(t) = y_c + r sin (omega t + phi_0)
148 // z(t) = z_0 + (p_Z0/gammam) t
149 x_t = x_c + r * cos(omega * t + phi_0);
150 y_t = y_c + r * sin(omega * t + phi_0);
151 z_t = Part->Z + Part->Pz / gammam * t;
152
153 // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t)
154 // R(t) = sqrt(x(t)² + y(t)²)
155 // Phi(t) = atan(y(t)/x(t))
156 // Theta(t) = atan(R(t)/z(t))
157 // Eta(t) = -ln tan (Theta(t)/2)
158 R_t = sqrt( pow(x_t,2.) + pow(y_t,2.) );
159 Phi_t = atan2( y_t, x_t);
[219]160 if(R_t>0) {
[199]161 Theta_t = acos( z_t / sqrt(z_t*z_t+ R_t*R_t));
162 Eta_t = - log(tan(Theta_t/2.));
163 } else{
164 Theta_t=0; Eta_t = 9999;
165 }
[219]166
[199]167 Px_t = - Part->PT * sin(omega*t + phi_0);
168 Py_t = Part->PT * cos(omega*t + phi_0);
169 Pz_t = Part->Pz;
170 PT_t = sqrt(Px_t*Px_t + Py_t*Py_t);
171 p_t = sqrt(PT_t*PT_t + Pz_t*Pz_t);
172 E_t=sqrt(Part->M*Part->M +p_t);
[219]173 //if(p_t != fabs(Pz_t) ) Eta_t = log( (p_t+Pz_t)/(p_t-Pz_t) )/2.;
174 //if(p_t>0) Theta_t = acos(Pz_t/p_t);
[199]175 momentum.SetPxPyPzE(Px_t,Py_t,Pz_t,E_t);
[193]176
[199]177// test zone ---
178/*
179 cout << cos(atan(R_t/z_t)) << "\t" << cos(Theta_t) << "\t" << cos(momentum.Theta()) << "\t" << Pz_t/temp_p << endl;
180 double Eta_t1 = log( (E+Pz_t)/(E-Pz_t) )/2.;
181 double Eta_t2 = log( (temp_p+Pz_t)/(temp_p-Pz_t) )/2.;
182 if(0 && fabs(Eta_t -Eta_t2)>1e-310) {
183 cout << "ERROR-BUG: Eta_t != Eta_t2\n";
184 cout << "Eta_t= " << Eta_t << "\t Eta_t1= " << Eta_t1 << "\t Eta_t2= " << Eta_t2 << endl;
[193]185 }
186
[199]187 double R_t2 = sqrt( pow(R_c,2.) + pow(r,2.) + 2*r*R_c*cos(phi_0 + omega*t - Phi_c) ); // cross-check
188 if(fabs(R_t - R_t2) > 1e-7)
189 cout << "ERROR-BUG: R_t != R_t2: R_t=" << R_t << " R_t2=" << R_t2 << " R_t - R_t2 =" << R_t - R_t2 << endl;
190 if( fabs(E - gammam) > 1e-3 ) {
[193]191 cout << "ERROR-BUG: energy is not conserved in src/BFieldProp.cc\n";
192 cout << "E - momentum.E() = " << fabs(E - momentum.E()) << " gammam - E " << fabs(gammam -E) << endl; }
193 if( fabs(PT_t - Part->PT) > 1e-10 ) {
[199]194 cout << "ERROR-BUG: PT is not conversed in src/BFieldProp.cc. ";
[193]195 cout << "(at " << 100*(PT_t - Part->PT) << "%)\n";
196 }
197 if(momentum.Pz() != Pz_t)
198 cout << "ERROR-BUG: Pz is not conserved in src/BFieldProp.cc\n";
199
[199]200 double temp_p0=sqrt(Part->PT*Part->PT + Part->Pz*Part->Pz);
201 if(fabs( (temp_p-temp_p0)*(temp_p+temp_p0) )>1e-10 ) {
202 cout << "ERROR-BUG: momentum |vec{p}| is not conserved in src/BFieldProp.cc\n";
203 cout << temp_p << "\t" << temp_p0 << endl;
204 }
205
206 // if x_c == y_c ==0 (set it by hand!), easy cross-check
207 //cout << "tan(phi_p)= " << momentum.Py()/momentum.Px() << "\t -1/tan(phi_x)= " << -x_t/y_t << endl;
208*/
209
[193]210 } else { // if B_x or B_y are non zero: longer computation
211
[53]212 float Xvertex1 = Part->X;
213 float Yvertex1 = Part->Y;
214 float Zvertex1 = Part->Z;
215
[193]216 //out of tracking coverage?
217 if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > R_max){return;}
218 if(fabs(Zvertex1) > z_max){return;}
[53]219
[193]220 double px = Part->Px / 0.003;
221 double py = Part->Py / 0.003;
222 double pz = Part->Pz / 0.003;
223 double pt = Part->PT / 0.003; // sqrt(px*px+py*py);
[199]224 double p = sqrt(pz*pz + pt*pt); //sqrt(px*px+py*py+pz*pz);
[59]225
[193]226 double M = Part->M;
227 double vx = px/M;
228 double vy = py/M;
229 double vz = pz/M;
230 double qm = q/M;
[53]231
[193]232 double ax = qm*(B_z*vy - B_y*vz);
233 double ay = qm*(B_x*vz - B_z*vx);
234 double az = qm*(B_y*vx - B_x*vy);
235 double dt = 1/p;
236 if(pt<266 && vz < 0.0012) dt = fabs(0.001/vz); // ?????
[59]237
[193]238 double xold=Xvertex1; double x=xold;
239 double yold=Yvertex1; double y=yold;
240 double zold=Zvertex1; double z=zold;
[59]241
[193]242 double VTold = pt/M; //=sqrt(vx*vx+vy*vy);
[59]243
[193]244 unsigned int k = 0;
245 double VTratio=0;
246 double R_max2 = R_max*R_max;
247 double r2=0; // will be x*x+y*y
[100]248
[193]249 while(k < MAXITERATION){
[59]250 k++;
251
252 vx += ax*dt;
253 vy += ay*dt;
254 vz += az*dt;
255
[193]256 VTratio = VTold/sqrt(vx*vx+vy*vy);
[59]257 vx *= VTratio;
258 vy *= VTratio;
259
[193]260 ax = qm*(B_z*vy - B_y*vz);
261 ay = qm*(B_x*vz - B_z*vx);
262 az = qm*(B_y*vx - B_x*vy);
[59]263
264 x += vx*dt;
265 y += vy*dt;
266 z += vz*dt;
[193]267 r2 = x*x + y*y;
[59]268
[193]269 if( r2 > R_max2 ){
270 x /= r2/R_max2;
271 y /= r2/R_max2;
272 break;
273 }
274 if( fabs(z)>z_max)break;
[59]275
276 xold = x;
277 yold = y;
278 zold = z;
[193]279 } // while loop
280
281 if(k == MAXITERATION) loop_overflow_counter++;
282 //cout << "too short loop in " << loop_overflow_counter << " cases" << endl;
[59]283
[193]284 if(x!=0 && y!=0 && z!=0) {
285 float Theta = atan2(sqrt(r2),z);
286 double eta = -log(tan(Theta/2.));
[53]287 double phi = atan2(y,x);
[193]288 momentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
289 }
290
291 } // if b_x or b_y non zero
[53]292}
[248]293
294
295
296void TrackPropagation::bfield(const TRootGenParticle *Part, float& etacalo, float& phicalo) {
297
298 // initialisation, valid for z_max==0, R_max==0 and q==0
299 etacalo = Part->Eta;
300 phicalo = -atan2(Part->Px,Part->Py);
301
302 q = Charge(Part->PID);
303 if(q==0) return;
304 if(R_max ==0) { cout << "ERROR: magnetic field has no lateral extention\n"; return;}
305 if(z_max==0) { cout << "ERROR: magnetic field has no longitudinal extention\n"; return;}
306
307 if (B_x== 0 && B_y== 0) { // faster if only B_z
308 if (B_z==0) return; // nothing to do
309
310 // initial conditions:
311 // p_X0 = Part->Px, p_Y0 = Part->Py, p_Z0 = Part->Pz, p_T0 = Part->PT;
312 // X_0 = Part->X, Y_0 = Part->Y, Z_0 = Part->Z;
313
314 // 1. initial transverse momentum p_{T0} : Part->PT
315 // initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0)
316 // relativistic gamma : gamma = E/mc² ; gammam = gamma \times m
317 // giration frequency \omega = q/(gamma m) B_z
318 // helix radius r = p_T0 / (omega gamma m)
319 phi_0 = -atan2(Part->Px,Part->Py);
320 gammam = Part->E; // here c==1
321 //cout << "gammam" << gammam << "\t gamma" << gammam/Part->M << endl;
322 omega = q * B_z /gammam;
323 r = Part->PT / (omega * gammam);
324
325 // 2. Helix parameters : center coordinates in transverse plane
326 // x_c = x_0 - r*cos(phi_0) and y_c = y_0 - r*sin(phi_0)
327 // R_c = \sqrt{x_c² + y_c²} and \Phi_c = atan{y_c/x_c}
328 x_c = Part->X - r*cos(phi_0); /// TEST !!
329 y_c = Part->Y - r*sin(phi_0);
330 R_c = sqrt(pow(x_c,2.) + pow(y_c,2.) );
331 Phi_c = atan2(y_c,x_c);
332
333 // 3. time evaluation t = min(t_T, t_z)
334 // t_T : time to exit from the sides
335 // t_T= [ Phi_c - phi_0 + atan( (R_max^2 - (R_c^2 + r^2))/(2rR_c) ) ]/omega
336 // t_z : time to exit from the front or the back
337 // t_z = gamma * m /p_z0 \times (-z_0 + z_max * sign(p_z0))
338 rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable
339 t_T=0;
340 int sign_pz= (Part->Pz >0) ? 1 : -1;
341 t_z = gammam / Part->Pz * (-Part->Z + z_max*sign_pz ) ;
342 if ( fabs(R_c - r) > R_max || R_c + r < R_max ) t = t_z;
343 else {
344 t_T = (Phi_c - phi_0 + atan2( (R_max + rr)*(R_max - rr) , 2*r*R_c ) ) / omega;
345 t = min(t_T,t_z);
346 }
347
348 // 4. position in terms of x(t), y(t), z(t)
349 // x(t) = x_c + r cos (omega t + phi_0)
350 // y(t) = y_c + r sin (omega t + phi_0)
351 // z(t) = z_0 + (p_Z0/gammam) t
352 x_t = x_c + r * cos(omega * t + phi_0);
353 y_t = y_c + r * sin(omega * t + phi_0);
354 z_t = Part->Z + Part->Pz / gammam * t;
355
356 // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t)
357 // R(t) = sqrt(x(t)² + y(t)²)
358 // Phi(t) = atan(y(t)/x(t))
359 // Theta(t) = atan(R(t)/z(t))
360 // Eta(t) = -ln tan (Theta(t)/2)
361 R_t = sqrt( pow(x_t,2.) + pow(y_t,2.) );
362 Phi_t = atan2( y_t, x_t);
363 if(R_t>0) {
364 Theta_t = acos( z_t / sqrt(z_t*z_t+ R_t*R_t));
365 Eta_t = - log(tan(Theta_t/2.));
366 } else{
367 Theta_t=0; Eta_t = 9999;
368 }
369/* Not needed here. but these formulae are correct -------
370 Px_t = - Part->PT * sin(omega*t + phi_0);
371 Py_t = Part->PT * cos(omega*t + phi_0);
372 Pz_t = Part->Pz;
373 PT_t = sqrt(Px_t*Px_t + Py_t*Py_t);
374 p_t = sqrt(PT_t*PT_t + Pz_t*Pz_t);
375 E_t=sqrt(Part->M*Part->M +p_t);
376 //if(p_t != fabs(Pz_t) ) Eta_t = log( (p_t+Pz_t)/(p_t-Pz_t) )/2.;
377 //if(p_t>0) Theta_t = acos(Pz_t/p_t);
378 momentum.SetPxPyPzE(Px_t,Py_t,Pz_t,E_t);
379*/
380 etacalo = Eta_t;
381 phicalo = Phi_t;
382 return;
383// test zone ---
384/*
385 cout << cos(atan(R_t/z_t)) << "\t" << cos(Theta_t) << "\t" << cos(momentum.Theta()) << "\t" << Pz_t/temp_p << endl;
386 double Eta_t1 = log( (E+Pz_t)/(E-Pz_t) )/2.;
387 double Eta_t2 = log( (temp_p+Pz_t)/(temp_p-Pz_t) )/2.;
388 if(0 && fabs(Eta_t -Eta_t2)>1e-310) {
389 cout << "ERROR-BUG: Eta_t != Eta_t2\n";
390 cout << "Eta_t= " << Eta_t << "\t Eta_t1= " << Eta_t1 << "\t Eta_t2= " << Eta_t2 << endl;
391 }
392
393 double R_t2 = sqrt( pow(R_c,2.) + pow(r,2.) + 2*r*R_c*cos(phi_0 + omega*t - Phi_c) ); // cross-check
394 if(fabs(R_t - R_t2) > 1e-7)
395 cout << "ERROR-BUG: R_t != R_t2: R_t=" << R_t << " R_t2=" << R_t2 << " R_t - R_t2 =" << R_t - R_t2 << endl;
396 if( fabs(E - gammam) > 1e-3 ) {
397 cout << "ERROR-BUG: energy is not conserved in src/BFieldProp.cc\n";
398 cout << "E - momentum.E() = " << fabs(E - momentum.E()) << " gammam - E " << fabs(gammam -E) << endl; }
399 if( fabs(PT_t - Part->PT) > 1e-10 ) {
400 cout << "ERROR-BUG: PT is not conversed in src/BFieldProp.cc. ";
401 cout << "(at " << 100*(PT_t - Part->PT) << "%)\n";
402 }
403 if(momentum.Pz() != Pz_t)
404 cout << "ERROR-BUG: Pz is not conserved in src/BFieldProp.cc\n";
405
406 double temp_p0=sqrt(Part->PT*Part->PT + Part->Pz*Part->Pz);
407 if(fabs( (temp_p-temp_p0)*(temp_p+temp_p0) )>1e-10 ) {
408 cout << "ERROR-BUG: momentum |vec{p}| is not conserved in src/BFieldProp.cc\n";
409 cout << temp_p << "\t" << temp_p0 << endl;
410 }
411
412 // if x_c == y_c ==0 (set it by hand!), easy cross-check
413 //cout << "tan(phi_p)= " << momentum.Py()/momentum.Px() << "\t -1/tan(phi_x)= " << -x_t/y_t << endl;
414*/
415
416 } else { // if B_x or B_y are non zero: longer computation
417
418 float Xvertex1 = Part->X;
419 float Yvertex1 = Part->Y;
420 float Zvertex1 = Part->Z;
421
422 //out of tracking coverage?
423 if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > R_max){return;}
424 if(fabs(Zvertex1) > z_max){return;}
425
426 double px = Part->Px / 0.003;
427 double py = Part->Py / 0.003;
428 double pz = Part->Pz / 0.003;
429 double pt = Part->PT / 0.003; // sqrt(px*px+py*py);
430 double p = sqrt(pz*pz + pt*pt); //sqrt(px*px+py*py+pz*pz);
431
432 double M = Part->M;
433 double vx = px/M;
434 double vy = py/M;
435 double vz = pz/M;
436 double qm = q/M;
437
438 double ax = qm*(B_z*vy - B_y*vz);
439 double ay = qm*(B_x*vz - B_z*vx);
440 double az = qm*(B_y*vx - B_x*vy);
441 double dt = 1/p;
442 if(pt<266 && vz < 0.0012) dt = fabs(0.001/vz); // ?????
443
444 double xold=Xvertex1; double x=xold;
445 double yold=Yvertex1; double y=yold;
446 double zold=Zvertex1; double z=zold;
447
448 double VTold = pt/M; //=sqrt(vx*vx+vy*vy);
449
450 unsigned int k = 0;
451 double VTratio=0;
452 double R_max2 = R_max*R_max;
453 double r2=0; // will be x*x+y*y
454
455 while(k < MAXITERATION){
456 k++;
457
458 vx += ax*dt;
459 vy += ay*dt;
460 vz += az*dt;
461
462 VTratio = VTold/sqrt(vx*vx+vy*vy);
463 vx *= VTratio;
464 vy *= VTratio;
465
466 ax = qm*(B_z*vy - B_y*vz);
467 ay = qm*(B_x*vz - B_z*vx);
468 az = qm*(B_y*vx - B_x*vy);
469
470 x += vx*dt;
471 y += vy*dt;
472 z += vz*dt;
473 r2 = x*x + y*y;
474
475 if( r2 > R_max2 ){
476 x /= r2/R_max2;
477 y /= r2/R_max2;
478 break;
479 }
480 if( fabs(z)>z_max)break;
481
482 xold = x;
483 yold = y;
484 zold = z;
485 } // while loop
486
487 if(k == MAXITERATION) loop_overflow_counter++;
488 //cout << "too short loop in " << loop_overflow_counter << " cases" << endl;
489 float Theta=0;
490 if(x!=0 && y!=0 && z!=0) {
491 Theta = atan2(sqrt(r2),z);
492 etacalo = -log(tan(Theta/2.));
493 phicalo = atan2(y,x);
494 //momentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
495 }
496 } // if b_x or b_y non zero
497}
Note: See TracBrowser for help on using the repository browser.