- Timestamp:
- Jun 17, 2009, 10:34:40 PM (15 years ago)
- File:
-
- 1 edited
-
trunk/src/BFieldProp.cc (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/BFieldProp.cc
r380 r434 109 109 // magnetic field parameters 110 110 R_max = DET->TRACK_radius/100.; //[m] 111 z_max = DET->TRACK_length/ 200.; //[m]111 z_max = DET->TRACK_length/100.; //[m] 112 112 B_x = DET->TRACK_bfield_x*tesla; 113 113 B_y = DET->TRACK_bfield_y*tesla; … … 132 132 if (!DET->FLAG_bfield ) return; 133 133 134 double M; 134 double M; // GeV/c² 135 135 //int q1 = ChargeVal(Part->PID) *eplus; // in units of 'e' 136 136 if(Part->M < -999) { // unitialised! 137 137 PdgParticle pdg_part = DET->PDGtable[Part->PID]; 138 138 q = pdg_part.charge() *eplus; // in units of 'e' 139 M = pdg_part.mass(); // GeV 139 M = pdg_part.mass(); // GeV/c² 140 140 } else { q = Part->Charge; M = Part->M; } 141 141 … … 154 154 if (B_x== 0 && B_y== 0) { // faster if only B_z 155 155 if (B_z==0) return; // nothing to do 156 157 156 158 157 //in test mode, just run once … … 168 167 // giration frequency \omega = q/(gamma m) B_z 169 168 // helix radius r = p_T0 / (omega gamma m) 170 double Px = Part->Px; // [GeV/c]169 double Px = Part->Px; // [GeV/c] 171 170 double Py = Part->Py; 172 171 double Pz = Part->Pz; 173 172 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 242 double delta= UNDEFINED; 243 244 if (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 312 else { 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 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 317 189 x_c = X + r*sin(phi_0); 318 190 y_c = Y - r*cos(phi_0); … … 321 193 Phi = Phi_c; 322 194 if(x_c<0) Phi += pi; 323 //if(Phi<0)Phi = 2*pi+Phi; // ne pas le mettre pour que le test1 fonctionne324 195 325 196 // 3. time evaluation t = min(t_T, t_z) 326 197 // t_T : time to exit from the sides 327 198 // t_z : time to exit from the front or the back 328 rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable329 t_T=0; //[ns]199 rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable [m] 200 t_T=0; //[ns] 330 201 int sign_pz= (Pz >0) ? 1 : -1; 331 202 if(Pz==0) t_z = 1E99; 332 else t_z = gammam / (Pz*1E -9* 2.99792458E+8) * (-Z + z_max*sign_pz );203 else t_z = gammam / (Pz*1E9/c_light) * (-Z + z_max*sign_pz ); 333 204 if( t_z <0) cout << "ERROR: t_z <0 !" << endl; 334 205 … … 365 236 double t_Tb = min(t4,min(t5,t6)); 366 237 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 238 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 239 } 378 240 } 379 //r = PT / (omega * gammam );380 241 381 242 // 4. position in terms of x(t), y(t), z(t) 382 243 x_t = x_c + r * sin(omega * t - phi_0); 383 244 y_t = y_c + r * cos(omega * t - phi_0); 384 z_t = Z + Pz*1E -9* 2.99792458E+8/ gammam * t;245 z_t = Z + Pz*1E9/c_light / gammam * t; 385 246 386 247 // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t) … … 392 253 } 393 254 else { 394 Theta_t=0; Eta_t = UNDEFINED; 255 Theta_t=0; 256 Eta_t = UNDEFINED; 395 257 } 396 } // method2397 258 398 259 if(test) { 399 260 cout << endl << endl; 400 cout << "method " << method << "----------------\n";401 261 cout << "x0,y0,z0= " << X << ", " << Y << ", " << Z << endl; 402 262 cout << "px0,py0,pz0= " << Px << ", " << Py << ", " << Pz << endl; … … 414 274 415 275 /* Not needed here. but these formulae are correct ------- 416 // method1 276 // method1 (removed) 417 277 Px_t = - PT * sin(omega*t + phi_0); 418 278 Py_t = PT * cos(omega*t + phi_0); … … 431 291 */ 432 292 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 293 Part->EtaCalo = Eta_t; 437 294 Part->PhiCalo = Phi_t;
Note:
See TracChangeset
for help on using the changeset viewer.
