Fork me on GitHub

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

Last change on this file since 262 was 260, checked in by severine ovyn, 16 years ago

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