Fork me on GitHub

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

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

new header in all files
File size: 15.0 KB
Line 
1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
6** \------------- arXiv:0903.2225v1 ------------/ **
7** **
8** **
9** This package uses: **
10** ------------------ **
11** ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86 **
12** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
13** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
14** FROG: [hep-ex/0901.2718v1] **
15** HepMC: Comput. Phys. Commun.134 (2001) 41 **
16** **
17** ------------------------------------------------------------------ **
18** **
19** Main authors: **
20** ------------- **
21** **
22** Severine Ovyn Xavier Rouby **
23** severine.ovyn@uclouvain.be xavier.rouby@cern **
24** **
25** Center for Particle Physics and Phenomenology (CP3) **
26** Universite catholique de Louvain (UCL) **
27** Louvain-la-Neuve, Belgium **
28** **
29** Copyright (C) 2008-2009, **
30** All rights reserved. **
31** **
32***********************************************************************/
33
34#include "BFieldProp.h"
35#include "PdgParticle.h"
36#include "SystemOfUnits.h"
37#include "PhysicalConstants.h"
38#include<cmath>
39using namespace std;
40
41
42//------------------------------------------------------------------------------
43extern const float UNDEFINED;
44
45TrackPropagation::TrackPropagation(){
46 DET = new RESOLution();
47 init();
48}
49
50TrackPropagation::TrackPropagation(const string& DetDatacard){
51 DET = new RESOLution();
52 DET->ReadDataCard(DetDatacard);
53 init();
54}
55
56TrackPropagation::TrackPropagation(const RESOLution* DetDatacard){
57 DET= new RESOLution(*DetDatacard);
58 init();
59}
60
61TrackPropagation::TrackPropagation(const TrackPropagation & tp){
62 MAXITERATION = tp.MAXITERATION;
63 DET = new RESOLution(*(tp.DET));
64 R_max = tp.R_max; z_max = tp.z_max;
65 B_x = tp.B_x; B_y = tp.B_y; B_z = tp.B_z;
66 q = tp.q; phi_0 = tp.phi_0;
67 gammam= tp.gammam; omega = tp.omega;
68 r = tp.r; rr = tp.rr;
69 x_c = tp.x_c; y_c = tp.y_c;
70 R_c = tp.R_c; Phi_c = tp.Phi_c;
71 t = tp.t; t_z = tp.t_z; t_T = tp.t_T;
72 x_t = tp.x_t; y_t = tp.y_t; z_t = tp.z_t;
73 R_t = tp.R_t; Phi_t = tp.Phi_t;
74 Theta_t=tp.Theta_t; Eta_t = tp.Eta_t;
75 Px_t = tp.Px_t; Py_t = tp.Py_t; Pz_t = tp.Pz_t;
76 PT_t = tp.PT_t; p_t = tp.p_t; E_t = tp.E_t;
77 loop_overflow_counter = tp.loop_overflow_counter;
78}
79
80TrackPropagation& TrackPropagation::operator=(const TrackPropagation & tp) {
81 if(this==&tp) return *this;
82 MAXITERATION = tp.MAXITERATION;
83 DET = new RESOLution(*(tp.DET));
84 R_max = tp.R_max; z_max = tp.z_max;
85 B_x = tp.B_x; B_y = tp.B_y; B_z = tp.B_z;
86 q = tp.q; phi_0 = tp.phi_0;
87 gammam= tp.gammam; omega = tp.omega;
88 r = tp.r; rr = tp.rr;
89 x_c = tp.x_c; y_c = tp.y_c;
90 R_c = tp.R_c; Phi_c = tp.Phi_c;
91 t = tp.t; t_z = tp.t_z; t_T = tp.t_T;
92 x_t = tp.x_t; y_t = tp.y_t; z_t = tp.z_t;
93 R_t = tp.R_t; Phi_t = tp.Phi_t;
94 Theta_t=tp.Theta_t; Eta_t = tp.Eta_t;
95 Px_t = tp.Px_t; Py_t = tp.Py_t; Pz_t = tp.Pz_t;
96 PT_t = tp.PT_t; p_t = tp.p_t; E_t = tp.E_t;
97 loop_overflow_counter = tp.loop_overflow_counter;
98 return *this;
99}
100
101void TrackPropagation::init() {
102 MAXITERATION = 10000;
103 q= UNDEFINED; phi_0= UNDEFINED; gammam= UNDEFINED; omega=UNDEFINED; r=UNDEFINED;
104 x_c=UNDEFINED; y_c=UNDEFINED; R_c=UNDEFINED; Phi_c=UNDEFINED;
105 rr=UNDEFINED; t=UNDEFINED; t_z=UNDEFINED; t_T=UNDEFINED;
106 x_t=UNDEFINED; y_t=UNDEFINED; z_t=UNDEFINED;
107 R_t=UNDEFINED; Phi_t=UNDEFINED; Theta_t=UNDEFINED; Eta_t=UNDEFINED;
108 Px_t=UNDEFINED; Py_t=UNDEFINED; Pz_t=UNDEFINED; PT_t=UNDEFINED; p_t=UNDEFINED; E_t=UNDEFINED;
109
110 // DET has been initialised in the constructors
111 // magnetic field parameters
112 R_max = DET->TRACK_radius/100.; //[m]
113 z_max = DET->TRACK_length/100.; //[m]
114 B_x = DET->TRACK_bfield_x*tesla;
115 B_y = DET->TRACK_bfield_y*tesla;
116 B_z = DET->TRACK_bfield_z;
117
118 loop_overflow_counter=0;
119}
120
121
122
123
124
125
126void TrackPropagation::bfield(TRootGenParticle *Part) {
127
128
129 // initialisation, valid for z_max==0, R_max==0 and q==0
130 Part->EtaCalo = Part->Eta;
131 Part->PhiCalo = Part->Phi;//-atan2(Part->Px,Part->Py);
132
133 // trivial cases
134 if (!DET->FLAG_bfield ) return;
135
136 double M; // GeV/c²
137 //int q1 = ChargeVal(Part->PID) *eplus; // in units of 'e'
138 if(Part->M < -999) { // unitialised!
139 PdgParticle pdg_part = DET->PDGtable[Part->PID];
140 q = pdg_part.charge() *eplus; // in units of 'e'
141 M = pdg_part.mass(); // GeV/c²
142 } else { q = Part->Charge; M = Part->M; }
143
144 if(q==0) return;
145 if(R_max==0) { cout << "ERROR: magnetic field has no lateral extention\n"; return;}
146 if(z_max==0) { cout << "ERROR: magnetic field has no longitudinal extention\n"; return;}
147
148 double X = Part->X/1000.;//[m]
149 double Y = Part->Y/1000.;//[m]
150 double Z = Part->Z/1000.;//[m]
151
152 // out of tracking coverage?
153 if(sqrt(X*X+Y*Y) > R_max){return;}
154 if(fabs(Z) > z_max){return;}
155
156 if (B_x== 0 && B_y== 0) { // faster if only B_z
157 if (B_z==0) return; // nothing to do
158
159 //in test mode, just run once
160 if (loop_overflow_counter) return;
161
162 // initial conditions:
163 // p_X0 = Part->Px, p_Y0 = Part->Py, p_Z0 = Part->Pz, p_T0 = Part->PT;
164 // X_0 = Part->X, Y_0 = Part->Y, Z_0 = Part->Z;
165
166 // 1. initial transverse momentum p_{T0} : Part->PT
167 // initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0)
168 // relativistic gamma : gamma = E/mc² ; gammam = gamma \times m
169 // giration frequency \omega = q/(gamma m) B_z
170 // helix radius r = p_T0 / (omega gamma m)
171 double Px = Part->Px; // [GeV/c]
172 double Py = Part->Py;
173 double Pz = Part->Pz;
174 double PT = Part->PT;
175 double E = Part->E; // [GeV]
176 double Phi = UNDEFINED;
177
178 float c_light = 2.99792458E+8;
179 gammam = E*1E9/(c_light*c_light); // gammam in [eV/c²]
180 omega = q * B_z / (gammam); // omega is here in [ 89875518 / s]
181 //cout << "omega*gammam = B_z in BFieldProp.cc: " << fabs(omega*gammam) - B_z << endl;
182 r = PT / (omega * gammam) *1E9/c_light; // in [m]
183
184 // test mode ?
185 bool test=false; if(test) loop_overflow_counter++;
186
187 double delta= UNDEFINED;
188 phi_0 = atan2(Py,Px); // [rad] in [-pi ; pi ]
189
190 // 2. helix axis coordinates
191 x_c = X + r*sin(phi_0);
192 y_c = Y - r*cos(phi_0);
193 R_c = sqrt( pow(x_c,2.) + pow(y_c,2.) );
194 Phi_c = atan2(y_c,x_c);
195 Phi = Phi_c;
196 if(x_c<0) Phi += pi;
197
198 // 3. time evaluation t = min(t_T, t_z)
199 // t_T : time to exit from the sides
200 // t_z : time to exit from the front or the back
201 rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable [m]
202 t_T=0; //[ns]
203 int sign_pz= (Pz >0) ? 1 : -1;
204 if(Pz==0) t_z = 1E99;
205 else t_z = gammam / (Pz*1E9/c_light) * (-Z + z_max*sign_pz );
206 if( t_z <0) cout << "ERROR: t_z <0 !" << endl;
207
208 if ( fabs(R_c - fabs(r)) > R_max || R_c + fabs(r) < R_max ) t = t_z;
209 else {
210 if(r==0) cout << "r ==0 !" << endl;
211 if(R_c==0) cout << "R_c ==0 !" << endl;
212 if(r==0|| R_c ==0) t_T=1E99;
213 else {
214 double asinrho = asin( (R_max + rr)*(R_max - rr) / (2*fabs(r)*R_c) );
215 delta = phi_0 - Phi;
216 if(delta<-pi) delta += 2*pi;
217 if(delta> pi) delta -= 2*pi;
218 double t1 = (delta + asinrho) / omega;
219 double t2 = (delta + pi - asinrho) / omega;
220 double t3 = (delta + pi + asinrho) / omega;
221 double t4 = (delta - asinrho) / omega;
222 double t5 = (delta - pi - asinrho) / omega;
223 double t6 = (delta - pi + asinrho) / omega;
224
225 if(test) {
226 cout << "t4 = " << t4 << "\t t5 = " << t5 << "\t t_6=" << t6 << "\t t_3 = " << t3 << endl;
227 cout << "t1 = " << t1 << "\t t2 = " << t2 << "\t t_T=" << t_T << "\t t_z = " << t_z << endl;
228 cout << "delta= " << delta << endl;
229 }
230 if(t1<0)t1=1E99;
231 if(t2<0)t2=1E99;
232 if(t3<0)t3=1E99;
233 if(t4<0)t4=1E99;
234 if(t5<0)t5=1E99;
235 if(t6<0)t6=1E99;
236
237 double t_Ta = min(t1,min(t2,t3));
238 double t_Tb = min(t4,min(t5,t6));
239 t_T = min(t_Ta,t_Tb);
240 t = min(t_T,t_z);
241 }
242 }
243
244 // 4. position in terms of x(t), y(t), z(t)
245 x_t = x_c + r * sin(omega * t - phi_0);
246 y_t = y_c + r * cos(omega * t - phi_0);
247 z_t = Z + Pz*1E9/c_light / gammam * t;
248
249 // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t)
250 R_t = sqrt( pow(x_t,2.) + pow(y_t,2.) );
251 Phi_t = atan2( y_t, x_t);
252 if(R_t>0) {
253 Theta_t = acos( z_t / sqrt(z_t*z_t+ R_t*R_t));
254 Eta_t = - log(tan(Theta_t/2.));
255 }
256 else {
257 Theta_t=0;
258 Eta_t = UNDEFINED;
259 }
260
261 if(test) {
262 cout << endl << endl;
263 cout << "x0,y0,z0= " << X << ", " << Y << ", " << Z << endl;
264 cout << "px0,py0,pz0= " << Px << ", " << Py << ", " << Pz << endl;
265 cout << "r = " << r << "R_max = " << R_max << "\t phi_0=" << phi_0 << endl;
266 cout << "gammam= " << gammam << "\t omega=" << omega << "\t PT = " << PT << endl;
267 cout << "x_c = " << x_c << "\t y_c = " << y_c << "\t R_c = " << R_c << "\t Phi = " << Phi << endl;
268 cout << "omega t = " << omega*t << "\t";
269 cout << "cos(omega t -phi0)= " << cos(omega*t-phi_0) << "\t sin(omega t -phi0)= " << sin(omega*t-phi_0) << endl;
270 cout << "t_T = " << t_T << "\t t_z = " << t_z << "\t r = " << r << endl;
271 cout << "x_t = " << x_t << "\t y_t = " << y_t << "\t z_t = " << z_t << endl;
272 cout << "R_t = " << R_t << "\t Phi_t = " << Phi_t << "\t";
273 cout << "Theta_t = " << Theta_t << "\t Eta_t = " << Eta_t << endl;
274 }
275
276
277/* Not needed here. but these formulae are correct -------
278 // method1 (removed)
279 Px_t = - PT * sin(omega*t + phi_0);
280 Py_t = PT * cos(omega*t + phi_0);
281
282 // method2
283 Px_t = PT * cos(phi_0 - omega*t);
284 Py_t = PT * sin(phi_0 - omega*t);
285
286 Pz_t = Pz;
287 PT_t = sqrt(Px_t*Px_t + Py_t*Py_t);
288 p_t = sqrt(PT_t*PT_t + Pz_t*Pz_t);
289 E_t=sqrt(M*M +p_t*p_t);
290 //if(p_t != fabs(Pz_t) ) Eta_t = log( (p_t+Pz_t)/(p_t-Pz_t) )/2.;
291 //if(p_t>0) Theta_t = acos(Pz_t/p_t)>;
292 momentum.SetPxPyPzE(Px_t,Py_t,Pz_t,E_t);
293*/
294
295 Part->EtaCalo = Eta_t;
296 Part->PhiCalo = Phi_t;
297
298// test zone ---
299/*
300
301 cout << "r = " << r << " et " << fabs(PT/(q*B_z)) << endl;
302 cout << cos(atan(R_t/z_t)) << "\t" << cos(Theta_t) << "\t" << cos(momentum.Theta()) << "\t" << Pz_t/temp_p << endl;
303 double Eta_t1 = log( (E+Pz_t)/(E-Pz_t) )/2.;
304 double Eta_t2 = log( (temp_p+Pz_t)/(temp_p-Pz_t) )/2.;
305 if(0 && fabs(Eta_t -Eta_t2)>1e-310) {
306 cout << "ERROR-BUG: Eta_t != Eta_t2\n";
307 cout << "Eta_t= " << Eta_t << "\t Eta_t1= " << Eta_t1 << "\t Eta_t2= " << Eta_t2 << endl;
308 }
309
310 double R_t2 = sqrt( pow(R_c,2.) + pow(r,2.) + 2*r*R_c*cos(phi_0 + omega*t - Phi_c) ); // cross-check
311 if(fabs(R_t - R_t2) > 1e-7)
312 cout << "ERROR-BUG: R_t != R_t2: R_t=" << R_t << " R_t2=" << R_t2 << " R_t - R_t2 =" << R_t - R_t2 << endl;
313 if( fabs(E - gammam) > 1e-3 ) {
314 cout << "ERROR-BUG: energy is not conserved in src/BFieldProp.cc\n";
315 cout << "E - momentum.E() = " << fabs(E - momentum.E()) << " gammam - E " << fabs(gammam -E) << endl; }
316 if( fabs(PT_t - Part->PT) > 1e-10 ) {
317 cout << "ERROR-BUG: PT is not conversed in src/BFieldProp.cc. ";
318 cout << "(at " << 100*(PT_t - Part->PT) << "%)\n";
319 }
320 if(momentum.Pz() != Pz_t)
321 cout << "ERROR-BUG: Pz is not conserved in src/BFieldProp.cc\n";
322
323 double temp_p0=sqrt(Part->PT*Part->PT + Part->Pz*Part->Pz);
324 if(fabs( (temp_p-temp_p0)*(temp_p+temp_p0) )>1e-10 ) {
325 cout << "ERROR-BUG: momentum |vec{p}| is not conserved in src/BFieldProp.cc\n";
326 cout << temp_p << "\t" << temp_p0 << endl;
327 }
328
329 // if x_c == y_c ==0 (set it by hand!), easy cross-check
330 //cout << "tan(phi_p)= " << momentum.Py()/momentum.Px() << "\t -1/tan(phi_x)= " << -x_t/y_t << endl;
331*/
332 return;
333
334 } else { // if B_x or B_y are non zero: longer computation
335//cout << "bfield de loic\n";
336 float Xvertex1 = Part->X;
337 float Yvertex1 = Part->Y;
338 float Zvertex1 = Part->Z;
339
340 double px = Part->Px / 0.003;
341 double py = Part->Py / 0.003;
342 double pz = Part->Pz / 0.003;
343 double pt = Part->PT / 0.003; // sqrt(px*px+py*py);
344 double p = sqrt(pz*pz + pt*pt); //sqrt(px*px+py*py+pz*pz);
345
346 //double M = Part->M; // see above
347 double vx = px/M;
348 double vy = py/M;
349 double vz = pz/M;
350 double qm = q/M;
351
352//double v = sqrt(vx*vx + vy*vy + vz*vz)/3E8;
353//cout << "v = " << v;
354//double gamma = 1./sqrt(1-v*v);
355//cout << "gamma = " << gamma << endl;
356
357 double ax = qm*(B_z*vy - B_y*vz);
358 double ay = qm*(B_x*vz - B_z*vx);
359 double az = qm*(B_y*vx - B_x*vy);
360 double dt = 1/p;
361 if(pt<266 && vz < 0.0012) dt = fabs(0.001/vz); // ?????
362
363 double xold=Xvertex1; double x=xold;
364 double yold=Yvertex1; double y=yold;
365 double zold=Zvertex1; double z=zold;
366
367 double VTold = pt/M; //=sqrt(vx*vx+vy*vy);
368
369 unsigned int k = 0;
370 double VTratio=0;
371 double R_max2 = R_max*R_max;
372 double r2=0; // will be x*x+y*y
373
374 while(k < MAXITERATION){
375 k++;
376
377 vx += ax*dt;
378 vy += ay*dt;
379 vz += az*dt;
380
381 VTratio = VTold/sqrt(vx*vx+vy*vy);
382 vx *= VTratio;
383 vy *= VTratio;
384
385 ax = qm*(B_z*vy - B_y*vz);
386 ay = qm*(B_x*vz - B_z*vx);
387 az = qm*(B_y*vx - B_x*vy);
388
389 x += vx*dt;
390 y += vy*dt;
391 z += vz*dt;
392 r2 = x*x + y*y;
393
394 if( r2 > R_max2 ){
395 x /= r2/R_max2;
396 y /= r2/R_max2;
397 break;
398 }
399 if( fabs(z)>z_max)break;
400
401 xold = x;
402 yold = y;
403 zold = z;
404 } // while loop
405
406 if(k == MAXITERATION) loop_overflow_counter++;
407 //cout << "too short loop in " << loop_overflow_counter << " cases" << endl;
408 float Theta=0;
409 if(x!=0 && y!=0 && z!=0) {
410 Theta = atan2(sqrt(r2),z);
411 Part->EtaCalo = -log(tan(Theta/2.));
412 Part->PhiCalo = atan2(y,x);
413 //momentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
414 }
415 } // if b_x or b_y non zero
416}
Note: See TracBrowser for help on using the repository browser.