Fork me on GitHub

Changeset 434 in svn


Ignore:
Timestamp:
Jun 17, 2009, 10:34:40 PM (15 years ago)
Author:
Xavier Rouby
Message:

new Bfield bug free

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/BFieldProp.cc

    r380 r434  
    109109  // magnetic field parameters
    110110   R_max = DET->TRACK_radius/100.; //[m]
    111    z_max = DET->TRACK_length/200.; //[m]
     111   z_max = DET->TRACK_length/100.; //[m]
    112112   B_x = DET->TRACK_bfield_x*tesla;
    113113   B_y = DET->TRACK_bfield_y*tesla;
    … …  
    132132  if (!DET->FLAG_bfield ) return;
    133133
    134   double M;
     134  double M; // GeV/c²
    135135  //int q1  = ChargeVal(Part->PID) *eplus; // in units of 'e'
    136136  if(Part->M < -999) { // unitialised!
    137137     PdgParticle pdg_part = DET->PDGtable[Part->PID];
    138138     q  = pdg_part.charge() *eplus; // in units of 'e'
    139      M  = pdg_part.mass(); // GeV
     139     M  = pdg_part.mass(); // GeV/c²
    140140  } else  { q = Part->Charge; M = Part->M; }
    141141
    … …  
    154154  if (B_x== 0 && B_y== 0) { // faster if only B_z
    155155    if (B_z==0) return; // nothing to do
    156 
    157156
    158157     //in test mode, just run once
    … …  
    168167    //     giration frequency \omega = q/(gamma m) B_z
    169168    //     helix radius r = p_T0 / (omega gamma m)
    170     double Px = Part->Px; // [GeV/c]
     169    double Px = Part->Px;                    // [GeV/c]
    171170    double Py = Part->Py;
    172171    double Pz = Part->Pz;
    173172    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
    317189      x_c = X + r*sin(phi_0);
    318190      y_c = Y - r*cos(phi_0);
    … …  
    321193      Phi = Phi_c;
    322194      if(x_c<0) Phi += pi;
    323       //if(Phi<0)Phi = 2*pi+Phi; // ne pas le mettre pour que le test1 fonctionne
    324195
    325196    // 3. time evaluation t = min(t_T, t_z)
    326197    //    t_T : time to exit from the sides
    327198    //    t_z : time to exit from the front or the back
    328       rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable
    329       t_T=0; //[ns]
     199      rr = sqrt( pow(R_c,2.) + pow(r,2.) );     // temp variable [m]
     200      t_T=0;                                    //[ns]
    330201      int sign_pz= (Pz >0) ? 1 : -1;
    331202      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 );
    333204      if( t_z <0) cout << "ERROR: t_z <0 !" << endl;
    334205
    … …  
    365236                double t_Tb = min(t4,min(t5,t6));
    366237                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";}
    372238                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                 }
    377239         }
    378240      }
    379 //r = PT / (omega * gammam );
    380241
    381242    // 4. position in terms of x(t), y(t), z(t)
    382243      x_t = x_c + r * sin(omega * t - phi_0);
    383244      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;
    385246
    386247    // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t)
    … …  
    392253      }
    393254      else {
    394                 Theta_t=0; Eta_t = UNDEFINED;
     255                Theta_t=0;
     256                Eta_t = UNDEFINED;
    395257      }
    396 } // method2
    397258
    398259     if(test) {
    399260        cout << endl << endl;
    400         cout << "method " << method << "----------------\n";
    401261        cout << "x0,y0,z0= " << X << ", " << Y << ", " << Z << endl;
    402262        cout << "px0,py0,pz0= " << Px << ", " << Py << ", " << Pz << endl;
    … …  
    414274     
    415275/*      Not needed here. but these formulae are correct -------
    416         // method1
     276        // method1 (removed)
    417277        Px_t = - PT * sin(omega*t + phi_0);
    418278        Py_t =   PT * cos(omega*t + phi_0);
    … …  
    431291*/
    432292
    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;
    436293        Part->EtaCalo = Eta_t;
    437294        Part->PhiCalo = Phi_t;
Note: See TracChangeset for help on using the changeset viewer.