Fork me on GitHub

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

Last change on this file since 343 was 294, checked in by Xavier Rouby, 15 years ago

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