Fork me on GitHub

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

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

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