Changeset 248 in svn
- Timestamp:
- Feb 5, 2009, 6:10:28 PM (16 years ago)
- Location:
- trunk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/interface/BFieldProp.h
r223 r248 30 30 ~TrackPropagation() {delete DET;}; 31 31 void init(); // for constructors 32 32 33 // Propagation and bfield are very similar. At the end, after code cleaning, 34 // onle bfield will remain in this class 33 35 void Propagation(const TRootGenParticle *Part,TLorentzVector &genMomentum); 36 void bfield(const TRootGenParticle *Part, float& etacalo, float& phicalo); 34 37 35 38 -
trunk/src/BFieldProp.cc
r219 r248 72 72 } 73 73 74 75 76 74 void TrackPropagation::init() { 77 75 MAXITERATION = 10000; … … 293 291 } // if b_x or b_y non zero 294 292 } 293 294 295 296 void TrackPropagation::bfield(const TRootGenParticle *Part, float& etacalo, float& phicalo) { 297 298 // initialisation, valid for z_max==0, R_max==0 and q==0 299 etacalo = Part->Eta; 300 phicalo = -atan2(Part->Px,Part->Py); 301 302 q = Charge(Part->PID); 303 if(q==0) return; 304 if(R_max ==0) { cout << "ERROR: magnetic field has no lateral extention\n"; return;} 305 if(z_max==0) { cout << "ERROR: magnetic field has no longitudinal extention\n"; return;} 306 307 if (B_x== 0 && B_y== 0) { // faster if only B_z 308 if (B_z==0) return; // nothing to do 309 310 // initial conditions: 311 // p_X0 = Part->Px, p_Y0 = Part->Py, p_Z0 = Part->Pz, p_T0 = Part->PT; 312 // X_0 = Part->X, Y_0 = Part->Y, Z_0 = Part->Z; 313 314 // 1. initial transverse momentum p_{T0} : Part->PT 315 // initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0) 316 // relativistic gamma : gamma = E/mc² ; gammam = gamma \times m 317 // giration frequency \omega = q/(gamma m) B_z 318 // helix radius r = p_T0 / (omega gamma m) 319 phi_0 = -atan2(Part->Px,Part->Py); 320 gammam = Part->E; // here c==1 321 //cout << "gammam" << gammam << "\t gamma" << gammam/Part->M << endl; 322 omega = q * B_z /gammam; 323 r = Part->PT / (omega * gammam); 324 325 // 2. Helix parameters : center coordinates in transverse plane 326 // x_c = x_0 - r*cos(phi_0) and y_c = y_0 - r*sin(phi_0) 327 // R_c = \sqrt{x_c² + y_c²} and \Phi_c = atan{y_c/x_c} 328 x_c = Part->X - r*cos(phi_0); /// TEST !! 329 y_c = Part->Y - r*sin(phi_0); 330 R_c = sqrt(pow(x_c,2.) + pow(y_c,2.) ); 331 Phi_c = atan2(y_c,x_c); 332 333 // 3. time evaluation t = min(t_T, t_z) 334 // t_T : time to exit from the sides 335 // t_T= [ Phi_c - phi_0 + atan( (R_max^2 - (R_c^2 + r^2))/(2rR_c) ) ]/omega 336 // t_z : time to exit from the front or the back 337 // t_z = gamma * m /p_z0 \times (-z_0 + z_max * sign(p_z0)) 338 rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable 339 t_T=0; 340 int sign_pz= (Part->Pz >0) ? 1 : -1; 341 t_z = gammam / Part->Pz * (-Part->Z + z_max*sign_pz ) ; 342 if ( fabs(R_c - r) > R_max || R_c + r < R_max ) t = t_z; 343 else { 344 t_T = (Phi_c - phi_0 + atan2( (R_max + rr)*(R_max - rr) , 2*r*R_c ) ) / omega; 345 t = min(t_T,t_z); 346 } 347 348 // 4. position in terms of x(t), y(t), z(t) 349 // x(t) = x_c + r cos (omega t + phi_0) 350 // y(t) = y_c + r sin (omega t + phi_0) 351 // z(t) = z_0 + (p_Z0/gammam) t 352 x_t = x_c + r * cos(omega * t + phi_0); 353 y_t = y_c + r * sin(omega * t + phi_0); 354 z_t = Part->Z + Part->Pz / gammam * t; 355 356 // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t) 357 // R(t) = sqrt(x(t)² + y(t)²) 358 // Phi(t) = atan(y(t)/x(t)) 359 // Theta(t) = atan(R(t)/z(t)) 360 // Eta(t) = -ln tan (Theta(t)/2) 361 R_t = sqrt( pow(x_t,2.) + pow(y_t,2.) ); 362 Phi_t = atan2( y_t, x_t); 363 if(R_t>0) { 364 Theta_t = acos( z_t / sqrt(z_t*z_t+ R_t*R_t)); 365 Eta_t = - log(tan(Theta_t/2.)); 366 } else{ 367 Theta_t=0; Eta_t = 9999; 368 } 369 /* Not needed here. but these formulae are correct ------- 370 Px_t = - Part->PT * sin(omega*t + phi_0); 371 Py_t = Part->PT * cos(omega*t + phi_0); 372 Pz_t = Part->Pz; 373 PT_t = sqrt(Px_t*Px_t + Py_t*Py_t); 374 p_t = sqrt(PT_t*PT_t + Pz_t*Pz_t); 375 E_t=sqrt(Part->M*Part->M +p_t); 376 //if(p_t != fabs(Pz_t) ) Eta_t = log( (p_t+Pz_t)/(p_t-Pz_t) )/2.; 377 //if(p_t>0) Theta_t = acos(Pz_t/p_t); 378 momentum.SetPxPyPzE(Px_t,Py_t,Pz_t,E_t); 379 */ 380 etacalo = Eta_t; 381 phicalo = Phi_t; 382 return; 383 // test zone --- 384 /* 385 cout << cos(atan(R_t/z_t)) << "\t" << cos(Theta_t) << "\t" << cos(momentum.Theta()) << "\t" << Pz_t/temp_p << endl; 386 double Eta_t1 = log( (E+Pz_t)/(E-Pz_t) )/2.; 387 double Eta_t2 = log( (temp_p+Pz_t)/(temp_p-Pz_t) )/2.; 388 if(0 && fabs(Eta_t -Eta_t2)>1e-310) { 389 cout << "ERROR-BUG: Eta_t != Eta_t2\n"; 390 cout << "Eta_t= " << Eta_t << "\t Eta_t1= " << Eta_t1 << "\t Eta_t2= " << Eta_t2 << endl; 391 } 392 393 double R_t2 = sqrt( pow(R_c,2.) + pow(r,2.) + 2*r*R_c*cos(phi_0 + omega*t - Phi_c) ); // cross-check 394 if(fabs(R_t - R_t2) > 1e-7) 395 cout << "ERROR-BUG: R_t != R_t2: R_t=" << R_t << " R_t2=" << R_t2 << " R_t - R_t2 =" << R_t - R_t2 << endl; 396 if( fabs(E - gammam) > 1e-3 ) { 397 cout << "ERROR-BUG: energy is not conserved in src/BFieldProp.cc\n"; 398 cout << "E - momentum.E() = " << fabs(E - momentum.E()) << " gammam - E " << fabs(gammam -E) << endl; } 399 if( fabs(PT_t - Part->PT) > 1e-10 ) { 400 cout << "ERROR-BUG: PT is not conversed in src/BFieldProp.cc. "; 401 cout << "(at " << 100*(PT_t - Part->PT) << "%)\n"; 402 } 403 if(momentum.Pz() != Pz_t) 404 cout << "ERROR-BUG: Pz is not conserved in src/BFieldProp.cc\n"; 405 406 double temp_p0=sqrt(Part->PT*Part->PT + Part->Pz*Part->Pz); 407 if(fabs( (temp_p-temp_p0)*(temp_p+temp_p0) )>1e-10 ) { 408 cout << "ERROR-BUG: momentum |vec{p}| is not conserved in src/BFieldProp.cc\n"; 409 cout << temp_p << "\t" << temp_p0 << endl; 410 } 411 412 // if x_c == y_c ==0 (set it by hand!), easy cross-check 413 //cout << "tan(phi_p)= " << momentum.Py()/momentum.Px() << "\t -1/tan(phi_x)= " << -x_t/y_t << endl; 414 */ 415 416 } else { // if B_x or B_y are non zero: longer computation 417 418 float Xvertex1 = Part->X; 419 float Yvertex1 = Part->Y; 420 float Zvertex1 = Part->Z; 421 422 //out of tracking coverage? 423 if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > R_max){return;} 424 if(fabs(Zvertex1) > z_max){return;} 425 426 double px = Part->Px / 0.003; 427 double py = Part->Py / 0.003; 428 double pz = Part->Pz / 0.003; 429 double pt = Part->PT / 0.003; // sqrt(px*px+py*py); 430 double p = sqrt(pz*pz + pt*pt); //sqrt(px*px+py*py+pz*pz); 431 432 double M = Part->M; 433 double vx = px/M; 434 double vy = py/M; 435 double vz = pz/M; 436 double qm = q/M; 437 438 double ax = qm*(B_z*vy - B_y*vz); 439 double ay = qm*(B_x*vz - B_z*vx); 440 double az = qm*(B_y*vx - B_x*vy); 441 double dt = 1/p; 442 if(pt<266 && vz < 0.0012) dt = fabs(0.001/vz); // ????? 443 444 double xold=Xvertex1; double x=xold; 445 double yold=Yvertex1; double y=yold; 446 double zold=Zvertex1; double z=zold; 447 448 double VTold = pt/M; //=sqrt(vx*vx+vy*vy); 449 450 unsigned int k = 0; 451 double VTratio=0; 452 double R_max2 = R_max*R_max; 453 double r2=0; // will be x*x+y*y 454 455 while(k < MAXITERATION){ 456 k++; 457 458 vx += ax*dt; 459 vy += ay*dt; 460 vz += az*dt; 461 462 VTratio = VTold/sqrt(vx*vx+vy*vy); 463 vx *= VTratio; 464 vy *= VTratio; 465 466 ax = qm*(B_z*vy - B_y*vz); 467 ay = qm*(B_x*vz - B_z*vx); 468 az = qm*(B_y*vx - B_x*vy); 469 470 x += vx*dt; 471 y += vy*dt; 472 z += vz*dt; 473 r2 = x*x + y*y; 474 475 if( r2 > R_max2 ){ 476 x /= r2/R_max2; 477 y /= r2/R_max2; 478 break; 479 } 480 if( fabs(z)>z_max)break; 481 482 xold = x; 483 yold = y; 484 zold = z; 485 } // while loop 486 487 if(k == MAXITERATION) loop_overflow_counter++; 488 //cout << "too short loop in " << loop_overflow_counter << " cases" << endl; 489 float Theta=0; 490 if(x!=0 && y!=0 && z!=0) { 491 Theta = atan2(sqrt(r2),z); 492 etacalo = -log(tan(Theta/2.)); 493 phicalo = atan2(y,x); 494 //momentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E); 495 } 496 } // if b_x or b_y non zero 497 }
Note:
See TracChangeset
for help on using the changeset viewer.