Fork me on GitHub

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

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

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