Fork me on GitHub

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

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

new algorithm inserted but not yet active
File size: 7.7 KB
Line 
1/*
2 * ---- Delphes ----
3 * A Fast Simulator for general purpose LHC detector
4 * S. Ovyn ~~~~ severine.ovyn@uclouvain.be
5 *
6 * Center for Particle Physics and Phenomenology (CP3)
7 * Universite Catholique de Louvain (UCL)
8 * Louvain-la-Neuve, Belgium
9 * */
10
11#include "interface/BFieldProp.h"
12#include<cmath>
13#include "TMath.h"
14using namespace std;
15
16
17//------------------------------------------------------------------------------
18
19TrackPropagation::TrackPropagation(const string DetDatacard):
20 MAXITERATION(10000), q(-9999.), phi_0(-9999.), gammam(-9999.), omega(-9999.), r(-9999.),
21 x_c(-9999.), y_c(-9999.), R_c(-9999.), Phi_c(-9999.),
22 rr(-9999.), t(-9999.), t_z(-9999.), t_T(-9999.),
23 x_t(-9999.), y_t(-9999.), z_t(-9999.),
24 R_t(-9999.), Phi_t(-9999.), Theta_t(-9999.), Eta_t(-9999.) {
25
26 // if(DetDatacard="") { DET = new RESOLution(); }
27 // else DET = new RESOLution(DetDatacard);
28 DET = new RESOLution();
29 DET->ReadDataCard(DetDatacard);
30
31 // magnetic field parameters
32 R_max = DET->TRACK_radius;
33 z_max = DET->TRACK_length/2.;
34 B_x = DET->TRACK_bfield_x;
35 B_y = DET->TRACK_bfield_y;
36 B_z = DET->TRACK_bfield_z;
37
38 loop_overflow_counter=0;
39}
40
41void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &momentum) {
42
43 q = Charge(Part->PID);
44 if(q==0) return;
45
46 if(R_max ==0) { cout << "ERROR: magnetic field has no lateral extention\n"; return;}
47 if(z_max==0) { cout << "ERROR: magnetic field has no longitudinal extention\n"; return;}
48
49 if (0 && B_x== 0 && B_y== 0) { // faster if only B_z
50 if (B_z==0) return; // nothing to do
51
52 // initial conditions:
53 // p_X0 = Part->Px, p_Y0 = Part->Py, p_Z0 = Part->Pz, p_T0 = Part->PT;
54 // X_0 = Part->X, Y_0 = Part->Y, Z_0 = Part->Z;
55
56 // 1. initial transverse momentum p_{T0} : Part->PT
57 // initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0)
58 // relativistic gamma : gamma = E/mc² ; gammam = gamma \times m
59 // giration frequency \omega = q/(gamma m) B_z
60 // helix radius r = p_T0 / (omega gamma m)
61 phi_0 = -atan2(Part->Px,Part->Py);
62 gammam = Part->E; // here c==1
63 //cout << "gammam" << gammam << "\t gamma" << gammam/Part->M << endl;
64 omega = q * B_z /gammam;
65 r = Part->PT / (omega * gammam);
66
67 // 2. Helix parameters : center coordinates in transverse plane
68 // x_c = x_0 - r*cos(phi_0) and y_c = y_0 - r*sin(phi_0)
69 // R_c = \sqrt{x_c² + y_c²} and \Phi_c = atan{y_c/x_c}
70 x_c = Part->X - r*cos(phi_0);
71 y_c = Part->Y - r*sin(phi_0);
72 R_c = sqrt(pow(x_c,2.) + pow(y_c,2.) );
73 Phi_c = atan2(y_c,x_c);
74
75 // 3. time evaluation t = min(t_T, t_z)
76 // t_T : time to exit from the sides
77 // t_z : time to exit from the front or the back
78 rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable
79 t_T=0;
80 t_z = gammam / Part->Pz * (-Part->Z + DET-> TRACK_length/2. * TMath::Sign((Float_t)1.,(Float_t)Part->Pz) ) ;
81 if ( fabs(R_c - r) > R_max || R_c + r < R_max ) t = t_z;
82 else {
83 t_T = (Phi_c - phi_0 + atan2( (R_max + rr)*(R_max - rr) , 2*r*R_c ) ) / omega;
84 t = min(t_T,t_z);
85 }
86
87 // 4. position in terms of x(t), y(t), z(t)
88 // x(t) = x_c + r cos (omega t + phi_0)
89 // y(t) = y_c + r sin (omega t + phi_0)
90 // z(t) = z_0 + (p_Z0/gammam) t
91 x_t = x_c + r * cos(omega * t + phi_0);
92 y_t = y_c + r * sin(omega * t + phi_0);
93 z_t = Part->Z + Part->Pz / gammam * t;
94
95 // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t)
96 // R(t) = sqrt(x(t)² + y(t)²)
97 // Phi(t) = atan(y(t)/x(t))
98 // Theta(t) = atan(R(t)/z(t))
99 // Eta(t) = -ln tan (Theta(t)/2)
100 R_t = sqrt( pow(x_t,2.) + pow(y_t,2.) );
101 double R_t2 = sqrt( pow(R_c,2.) + pow(r,2.) + 2*r*R_c*cos(phi_0 + omega*t - Phi_c) ); // cross-check
102 if(fabs(R_t - R_t2) > 1e-11)
103 cout << "ERROR-BUG: R_t != R_t2: R_t=" << R_t << " R_t2=" << R_t2 << " R_t - R_t2 =" << R_t - R_t2 << endl;
104 Phi_t = atan2( y_t, x_t);
105 Theta_t = atan2( R_t , z_t);
106 Eta_t = - log(tan(Theta_t/2.));
107
108 double Px_t = - Part->PT * sin(omega*t + phi_0);
109 double Py_t = Part->PT * cos(omega*t + phi_0);
110 double Pz_t = Part->Pz;
111 double PT_t = sqrt(Px_t*Px_t + Py_t*Py_t);
112 double E=sqrt(Part->M*Part->M + PT_t*PT_t + Pz_t*Pz_t); // cross-check
113 momentum.SetPxPyPzE(Px_t,Py_t,Pz_t,E);
114 TLorentzVector mom_temp;
115 mom_temp.SetPtEtaPhiM(PT_t,Eta_t,Phi_t+PI,Part->M);
116
117 if (momentum != mom_temp) {
118// cout << "momentum et mom_temp different\n";
119/* if (momentum.Px() != mom_temp.Px() ) cout << " px: " << momentum.Px() << " & " << mom_temp.Px() << endl;
120 if (momentum.Py() != mom_temp.Py() ) cout << " py: " << momentum.Py() << " & " << mom_temp.Py() << endl;
121 if (momentum.Pz() != mom_temp.Pz() ) cout << " pz: " << momentum.Pz() << " & " << mom_temp.Pz() << endl;
122 if (momentum.Pt() != mom_temp.Pt() ) cout << " pt: " << momentum.Pt() << " & " << mom_temp.Pt() << endl;*/
123// cout << "eta: " << momentum.Eta() << " " << mom_temp.Eta() << endl;
124// cout << "the: " << momentum.Theta() << " " << mom_temp.Theta() << endl;
125// cout << "phi: " << momentum.Phi() << " " << mom_temp.Phi() << endl;
126 }
127 //cout << Part->PT << " or " << momentum.Pt() << " or " << mom_temp.Pt() << endl;
128
129 if( fabs(E - gammam) > 1e-4 ) {
130 cout << "ERROR-BUG: energy is not conserved in src/BFieldProp.cc\n";
131 cout << "E - momentum.E() = " << fabs(E - momentum.E()) << " gammam - E " << fabs(gammam -E) << endl; }
132 if( fabs(PT_t - Part->PT) > 1e-10 ) {
133 cout << "ERROR-BUG: PT is not conserved in src/BFieldProp.cc. ";
134 cout << "(at " << 100*(PT_t - Part->PT) << "%)\n";
135 }
136 if(momentum.Pz() != Pz_t)
137 cout << "ERROR-BUG: Pz is not conserved in src/BFieldProp.cc\n";
138
139 //momentum.SetPtEtaPhiE(Part->PT,Eta_t,Phi_t,Part->E);
140
141 } else { // if B_x or B_y are non zero: longer computation
142
143 float Xvertex1 = Part->X;
144 float Yvertex1 = Part->Y;
145 float Zvertex1 = Part->Z;
146
147 //out of tracking coverage?
148 if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > R_max){return;}
149 if(fabs(Zvertex1) > z_max){return;}
150
151 double px = Part->Px / 0.003;
152 double py = Part->Py / 0.003;
153 double pz = Part->Pz / 0.003;
154 double pt = Part->PT / 0.003; // sqrt(px*px+py*py);
155 double p = sqrt(px*px+py*py+pz*pz);
156
157 double M = Part->M;
158 double vx = px/M;
159 double vy = py/M;
160 double vz = pz/M;
161 double qm = q/M;
162
163 double ax = qm*(B_z*vy - B_y*vz);
164 double ay = qm*(B_x*vz - B_z*vx);
165 double az = qm*(B_y*vx - B_x*vy);
166 double dt = 1/p;
167 if(pt<266 && vz < 0.0012) dt = fabs(0.001/vz); // ?????
168
169 double xold=Xvertex1; double x=xold;
170 double yold=Yvertex1; double y=yold;
171 double zold=Zvertex1; double z=zold;
172
173 double VTold = pt/M; //=sqrt(vx*vx+vy*vy);
174
175 unsigned int k = 0;
176 double VTratio=0;
177 double R_max2 = R_max*R_max;
178 double r2=0; // will be x*x+y*y
179
180 while(k < MAXITERATION){
181 k++;
182
183 vx += ax*dt;
184 vy += ay*dt;
185 vz += az*dt;
186
187 VTratio = VTold/sqrt(vx*vx+vy*vy);
188 vx *= VTratio;
189 vy *= VTratio;
190
191 ax = qm*(B_z*vy - B_y*vz);
192 ay = qm*(B_x*vz - B_z*vx);
193 az = qm*(B_y*vx - B_x*vy);
194
195 x += vx*dt;
196 y += vy*dt;
197 z += vz*dt;
198 r2 = x*x + y*y;
199
200 if( r2 > R_max2 ){
201 x /= r2/R_max2;
202 y /= r2/R_max2;
203 break;
204 }
205 if( fabs(z)>z_max)break;
206
207 xold = x;
208 yold = y;
209 zold = z;
210 } // while loop
211
212 if(k == MAXITERATION) loop_overflow_counter++;
213 //cout << "too short loop in " << loop_overflow_counter << " cases" << endl;
214
215 if(x!=0 && y!=0 && z!=0) {
216 float Theta = atan2(sqrt(r2),z);
217 double eta = -log(tan(Theta/2.));
218 double phi = atan2(y,x);
219 momentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
220 }
221
222 } // if b_x or b_y non zero
223}
Note: See TracBrowser for help on using the repository browser.