Fork me on GitHub

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

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

new Bfield bug free
File size: 14.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 "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/100.; //[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; // GeV/c²
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/c²
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 //in test mode, just run once
158 if (loop_overflow_counter) return;
159
160 // initial conditions:
161 // p_X0 = Part->Px, p_Y0 = Part->Py, p_Z0 = Part->Pz, p_T0 = Part->PT;
162 // X_0 = Part->X, Y_0 = Part->Y, Z_0 = Part->Z;
163
164 // 1. initial transverse momentum p_{T0} : Part->PT
165 // initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0)
166 // relativistic gamma : gamma = E/mc² ; gammam = gamma \times m
167 // giration frequency \omega = q/(gamma m) B_z
168 // helix radius r = p_T0 / (omega gamma m)
169 double Px = Part->Px; // [GeV/c]
170 double Py = Part->Py;
171 double Pz = Part->Pz;
172 double PT = Part->PT;
173 double E = Part->E; // [GeV]
174 double Phi = UNDEFINED;
175
176 float c_light = 2.99792458E+8;
177 gammam = E*1E9/(c_light*c_light); // gammam in [eV/c²]
178 omega = q * B_z / (gammam); // omega is here in [ 89875518 / s]
179 //cout << "omega*gammam = B_z in BFieldProp.cc: " << fabs(omega*gammam) - B_z << endl;
180 r = PT / (omega * gammam) *1E9/c_light; // in [m]
181
182 // test mode ?
183 bool test=false; if(test) loop_overflow_counter++;
184
185 double delta= UNDEFINED;
186 phi_0 = atan2(Py,Px); // [rad] in [-pi ; pi ]
187
188 // 2. helix axis coordinates
189 x_c = X + r*sin(phi_0);
190 y_c = Y - r*cos(phi_0);
191 R_c = sqrt( pow(x_c,2.) + pow(y_c,2.) );
192 Phi_c = atan2(y_c,x_c);
193 Phi = Phi_c;
194 if(x_c<0) Phi += pi;
195
196 // 3. time evaluation t = min(t_T, t_z)
197 // t_T : time to exit from the sides
198 // t_z : time to exit from the front or the back
199 rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable [m]
200 t_T=0; //[ns]
201 int sign_pz= (Pz >0) ? 1 : -1;
202 if(Pz==0) t_z = 1E99;
203 else t_z = gammam / (Pz*1E9/c_light) * (-Z + z_max*sign_pz );
204 if( t_z <0) cout << "ERROR: t_z <0 !" << endl;
205
206 if ( fabs(R_c - fabs(r)) > R_max || R_c + fabs(r) < R_max ) t = t_z;
207 else {
208 if(r==0) cout << "r ==0 !" << endl;
209 if(R_c==0) cout << "R_c ==0 !" << endl;
210 if(r==0|| R_c ==0) t_T=1E99;
211 else {
212 double asinrho = asin( (R_max + rr)*(R_max - rr) / (2*fabs(r)*R_c) );
213 delta = phi_0 - Phi;
214 if(delta<-pi) delta += 2*pi;
215 if(delta> pi) delta -= 2*pi;
216 double t1 = (delta + asinrho) / omega;
217 double t2 = (delta + pi - asinrho) / omega;
218 double t3 = (delta + pi + asinrho) / omega;
219 double t4 = (delta - asinrho) / omega;
220 double t5 = (delta - pi - asinrho) / omega;
221 double t6 = (delta - pi + asinrho) / omega;
222
223 if(test) {
224 cout << "t4 = " << t4 << "\t t5 = " << t5 << "\t t_6=" << t6 << "\t t_3 = " << t3 << endl;
225 cout << "t1 = " << t1 << "\t t2 = " << t2 << "\t t_T=" << t_T << "\t t_z = " << t_z << endl;
226 cout << "delta= " << delta << endl;
227 }
228 if(t1<0)t1=1E99;
229 if(t2<0)t2=1E99;
230 if(t3<0)t3=1E99;
231 if(t4<0)t4=1E99;
232 if(t5<0)t5=1E99;
233 if(t6<0)t6=1E99;
234
235 double t_Ta = min(t1,min(t2,t3));
236 double t_Tb = min(t4,min(t5,t6));
237 t_T = min(t_Ta,t_Tb);
238 t = min(t_T,t_z);
239 }
240 }
241
242 // 4. position in terms of x(t), y(t), z(t)
243 x_t = x_c + r * sin(omega * t - phi_0);
244 y_t = y_c + r * cos(omega * t - phi_0);
245 z_t = Z + Pz*1E9/c_light / gammam * t;
246
247 // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t)
248 R_t = sqrt( pow(x_t,2.) + pow(y_t,2.) );
249 Phi_t = atan2( y_t, x_t);
250 if(R_t>0) {
251 Theta_t = acos( z_t / sqrt(z_t*z_t+ R_t*R_t));
252 Eta_t = - log(tan(Theta_t/2.));
253 }
254 else {
255 Theta_t=0;
256 Eta_t = UNDEFINED;
257 }
258
259 if(test) {
260 cout << endl << endl;
261 cout << "x0,y0,z0= " << X << ", " << Y << ", " << Z << endl;
262 cout << "px0,py0,pz0= " << Px << ", " << Py << ", " << Pz << endl;
263 cout << "r = " << r << "R_max = " << R_max << "\t phi_0=" << phi_0 << endl;
264 cout << "gammam= " << gammam << "\t omega=" << omega << "\t PT = " << PT << endl;
265 cout << "x_c = " << x_c << "\t y_c = " << y_c << "\t R_c = " << R_c << "\t Phi = " << Phi << endl;
266 cout << "omega t = " << omega*t << "\t";
267 cout << "cos(omega t -phi0)= " << cos(omega*t-phi_0) << "\t sin(omega t -phi0)= " << sin(omega*t-phi_0) << endl;
268 cout << "t_T = " << t_T << "\t t_z = " << t_z << "\t r = " << r << endl;
269 cout << "x_t = " << x_t << "\t y_t = " << y_t << "\t z_t = " << z_t << endl;
270 cout << "R_t = " << R_t << "\t Phi_t = " << Phi_t << "\t";
271 cout << "Theta_t = " << Theta_t << "\t Eta_t = " << Eta_t << endl;
272 }
273
274
275/* Not needed here. but these formulae are correct -------
276 // method1 (removed)
277 Px_t = - PT * sin(omega*t + phi_0);
278 Py_t = PT * cos(omega*t + phi_0);
279
280 // method2
281 Px_t = PT * cos(phi_0 - omega*t);
282 Py_t = PT * sin(phi_0 - omega*t);
283
284 Pz_t = Pz;
285 PT_t = sqrt(Px_t*Px_t + Py_t*Py_t);
286 p_t = sqrt(PT_t*PT_t + Pz_t*Pz_t);
287 E_t=sqrt(M*M +p_t*p_t);
288 //if(p_t != fabs(Pz_t) ) Eta_t = log( (p_t+Pz_t)/(p_t-Pz_t) )/2.;
289 //if(p_t>0) Theta_t = acos(Pz_t/p_t)>;
290 momentum.SetPxPyPzE(Px_t,Py_t,Pz_t,E_t);
291*/
292
293 Part->EtaCalo = Eta_t;
294 Part->PhiCalo = Phi_t;
295
296// test zone ---
297/*
298
299 cout << "r = " << r << " et " << fabs(PT/(q*B_z)) << endl;
300 cout << cos(atan(R_t/z_t)) << "\t" << cos(Theta_t) << "\t" << cos(momentum.Theta()) << "\t" << Pz_t/temp_p << endl;
301 double Eta_t1 = log( (E+Pz_t)/(E-Pz_t) )/2.;
302 double Eta_t2 = log( (temp_p+Pz_t)/(temp_p-Pz_t) )/2.;
303 if(0 && fabs(Eta_t -Eta_t2)>1e-310) {
304 cout << "ERROR-BUG: Eta_t != Eta_t2\n";
305 cout << "Eta_t= " << Eta_t << "\t Eta_t1= " << Eta_t1 << "\t Eta_t2= " << Eta_t2 << endl;
306 }
307
308 double R_t2 = sqrt( pow(R_c,2.) + pow(r,2.) + 2*r*R_c*cos(phi_0 + omega*t - Phi_c) ); // cross-check
309 if(fabs(R_t - R_t2) > 1e-7)
310 cout << "ERROR-BUG: R_t != R_t2: R_t=" << R_t << " R_t2=" << R_t2 << " R_t - R_t2 =" << R_t - R_t2 << endl;
311 if( fabs(E - gammam) > 1e-3 ) {
312 cout << "ERROR-BUG: energy is not conserved in src/BFieldProp.cc\n";
313 cout << "E - momentum.E() = " << fabs(E - momentum.E()) << " gammam - E " << fabs(gammam -E) << endl; }
314 if( fabs(PT_t - Part->PT) > 1e-10 ) {
315 cout << "ERROR-BUG: PT is not conversed in src/BFieldProp.cc. ";
316 cout << "(at " << 100*(PT_t - Part->PT) << "%)\n";
317 }
318 if(momentum.Pz() != Pz_t)
319 cout << "ERROR-BUG: Pz is not conserved in src/BFieldProp.cc\n";
320
321 double temp_p0=sqrt(Part->PT*Part->PT + Part->Pz*Part->Pz);
322 if(fabs( (temp_p-temp_p0)*(temp_p+temp_p0) )>1e-10 ) {
323 cout << "ERROR-BUG: momentum |vec{p}| is not conserved in src/BFieldProp.cc\n";
324 cout << temp_p << "\t" << temp_p0 << endl;
325 }
326
327 // if x_c == y_c ==0 (set it by hand!), easy cross-check
328 //cout << "tan(phi_p)= " << momentum.Py()/momentum.Px() << "\t -1/tan(phi_x)= " << -x_t/y_t << endl;
329*/
330 return;
331
332 } else { // if B_x or B_y are non zero: longer computation
333//cout << "bfield de loic\n";
334 float Xvertex1 = Part->X;
335 float Yvertex1 = Part->Y;
336 float Zvertex1 = Part->Z;
337
338 double px = Part->Px / 0.003;
339 double py = Part->Py / 0.003;
340 double pz = Part->Pz / 0.003;
341 double pt = Part->PT / 0.003; // sqrt(px*px+py*py);
342 double p = sqrt(pz*pz + pt*pt); //sqrt(px*px+py*py+pz*pz);
343
344 //double M = Part->M; // see above
345 double vx = px/M;
346 double vy = py/M;
347 double vz = pz/M;
348 double qm = q/M;
349
350//double v = sqrt(vx*vx + vy*vy + vz*vz)/3E8;
351//cout << "v = " << v;
352//double gamma = 1./sqrt(1-v*v);
353//cout << "gamma = " << gamma << endl;
354
355 double ax = qm*(B_z*vy - B_y*vz);
356 double ay = qm*(B_x*vz - B_z*vx);
357 double az = qm*(B_y*vx - B_x*vy);
358 double dt = 1/p;
359 if(pt<266 && vz < 0.0012) dt = fabs(0.001/vz); // ?????
360
361 double xold=Xvertex1; double x=xold;
362 double yold=Yvertex1; double y=yold;
363 double zold=Zvertex1; double z=zold;
364
365 double VTold = pt/M; //=sqrt(vx*vx+vy*vy);
366
367 unsigned int k = 0;
368 double VTratio=0;
369 double R_max2 = R_max*R_max;
370 double r2=0; // will be x*x+y*y
371
372 while(k < MAXITERATION){
373 k++;
374
375 vx += ax*dt;
376 vy += ay*dt;
377 vz += az*dt;
378
379 VTratio = VTold/sqrt(vx*vx+vy*vy);
380 vx *= VTratio;
381 vy *= VTratio;
382
383 ax = qm*(B_z*vy - B_y*vz);
384 ay = qm*(B_x*vz - B_z*vx);
385 az = qm*(B_y*vx - B_x*vy);
386
387 x += vx*dt;
388 y += vy*dt;
389 z += vz*dt;
390 r2 = x*x + y*y;
391
392 if( r2 > R_max2 ){
393 x /= r2/R_max2;
394 y /= r2/R_max2;
395 break;
396 }
397 if( fabs(z)>z_max)break;
398
399 xold = x;
400 yold = y;
401 zold = z;
402 } // while loop
403
404 if(k == MAXITERATION) loop_overflow_counter++;
405 //cout << "too short loop in " << loop_overflow_counter << " cases" << endl;
406 float Theta=0;
407 if(x!=0 && y!=0 && z!=0) {
408 Theta = atan2(sqrt(r2),z);
409 Part->EtaCalo = -log(tan(Theta/2.));
410 Part->PhiCalo = atan2(y,x);
411 //momentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
412 }
413 } // if b_x or b_y non zero
414}
Note: See TracBrowser for help on using the repository browser.