Changeset 294 in svn

Timestamp:

Mar 9, 2009, 12:30:27 AM (16 years ago)

Author:

Xavier Rouby

Message:

bug free bfield?

File:

• trunk/src/BFieldProp.cc (modified) (8 diffs)

trunk/src/BFieldProp.cc

@@ -31 +31 @@
 
 #include "BFieldProp.h"
+#include "SystemOfUnits.h"
+#include "PhysicalConstants.h"
 #include<cmath>
 using namespace std;
@@ -105 +107 @@
   // DET has been initialised in the constructors
   // magnetic field parameters
-   R_max = DET->TRACK_radius;
-   z_max = DET->TRACK_length/2.;
-   B_x = DET->TRACK_bfield_x;
-   B_y = DET->TRACK_bfield_y;
+   R_max = DET->TRACK_radius/100.; //[m]
+   z_max = DET->TRACK_length/200.; //[m]
+   B_x = DET->TRACK_bfield_x*tesla;
+   B_y = DET->TRACK_bfield_y*tesla;
    B_z = DET->TRACK_bfield_z;
 
@@ -116 +118 @@
 
 
-void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &momentum) {
-
-  q  = ChargeVal(Part->PID);
+
+
+
+void TrackPropagation::bfield(TRootGenParticle *Part) {
+
+
+  // initialisation, valid for z_max==0, R_max==0 and q==0
+  Part->EtaCalo = Part->Eta;
+  Part->PhiCalo = Part->Phi;//-atan2(Part->Px,Part->Py);
+
+  // trivial cases
+  if (!DET->FLAG_bfield ) return;
+
+  q  = ChargeVal(Part->PID) *eplus; // in units of 'e'
   if(q==0) return;
-
-  if(R_max ==0) { cout << "ERROR: magnetic field has no lateral extention\n"; return;}
+  if(R_max==0) { cout << "ERROR: magnetic field has no lateral extention\n"; return;}
   if(z_max==0) { cout << "ERROR: magnetic field has no longitudinal extention\n"; return;}
+
+  double X  = Part->X/1000.;//[m] 
+  double Y  = Part->Y/1000.;//[m]
+  double Z  = Part->Z/1000.;//[m]
+
+  // out of tracking coverage?
+  if(sqrt(X*X+Y*Y) > R_max){return;}
+  if(fabs(Z) > z_max){return;}
 
   if (B_x== 0 && B_y== 0) { // faster if only B_z
     if (B_z==0) return; // nothing to do
 
+
+     //in test mode, just run once
+     if (loop_overflow_counter) return;
+  
     // initial conditions:
       // p_X0 = Part->Px, p_Y0 = Part->Py, p_Z0 = Part->Pz, p_T0 = Part->PT;
@@ -136 +160 @@
     //     giration frequency \omega = q/(gamma m) B_z
     //     helix radius r = p_T0 / (omega gamma m)
-      phi_0 = -atan2(Part->Px,Part->Py);
-      gammam = Part->E; // here c==1
-      //cout << "gammam" << gammam << "\t gamma" << gammam/Part->M << endl;
-      omega = q * B_z /gammam;
-      //r = Part->PT / (omega * gammam);
-      r = fabs(Part->PT / (omega * gammam));
+    double Px = Part->Px; // [GeV/c]
+    double Py = Part->Py;
+    double Pz = Part->Pz;
+    double PT = Part->PT;
+    double E  = Part->E;              // [GeV]
+    double M  = Part->M;              // [GeV]/c²
+      double Phi = UNDEFINED;
+
+      unsigned int method =2; 
+      gammam = E * 1E-9 ;                   // gammam in [eV]
+      omega = q * e_SI * B_z *  2.99792458E+8* 2.99792458E+8 / gammam;  // omega is here in [rad/s]  
+      r = PT * 1E-9 *  2.99792458E+8 / (omega * gammam );             // in [m]  m2 /( s)
+
+      // test mode
+      bool test=false; if(test) loop_overflow_counter++;
+
+    // test -- point 1) // tests faciles: changer le signe de q et de Py
+    if(test && false) {
+      q = -e_SI; X = Y = Z = Px = Pz = M = 0; 
+      Py= R_max * (q*B_z); 
+      PT = sqrt(Px*Px + Py*Py + Pz*Pz); 
+      E = PT* 2.99792458E+8; gammam = PT/ 2.99792458E+8;
+      omega = q/e_SI * 2.99792458E+8/R_max; // omega has a sign!
+      r = PT / (omega * gammam );           // r has a sign!
+    }
+    // test -- point 2)
+    if(test && false) { 
+      q = e_SI; X = R_max/2.; Y = Z = Px = Pz = M = 0; 
+      Py= -R_max * (q*B_z);
+      PT = sqrt(Px*Px + Py*Py + Pz*Pz);
+      E = PT* 2.99792458E+8; gammam = PT/ 2.99792458E+8;
+      omega = q/e_SI * 2.99792458E+8/R_max;
+      r = PT / (omega * gammam );
+    }
+    // test -- point 3)
+    if(test && false) {
+      q = -e_SI; X = 0; Y= -R_max/2.; Z = Px = Pz = M = 0;
+      Py= R_max * (q*B_z);
+      PT = sqrt(Px*Px + Py*Py + Pz*Pz);
+      E = PT* 2.99792458E+8; gammam = PT/ 2.99792458E+8;
+      omega = q/e_SI * 2.99792458E+8/R_max;
+      r = PT / (omega * gammam );
+    }
+    // test -- point 4)
+    if(test && false) {
+      q = -e_SI; X = R_max/2.; Y = Z = Py = Pz = M = 0;
+      Px= R_max * (q*B_z);
+      PT = sqrt(Px*Px + Py*Py + Pz*Pz);
+      E = PT* 2.99792458E+8; gammam = PT/ 2.99792458E+8;
+      omega = q/e_SI * 2.99792458E+8/R_max;
+      r = PT / (omega * gammam );
+    }
+    // test -- point 5)
+    if(test && false) {
+      q = e_SI; X = -R_max/2.; Y = Z = Px = Pz = M = 0;
+      Py= -R_max * (q*B_z);
+      PT = sqrt(Px*Px + Py*Py + Pz*Pz);
+      E = PT* 2.99792458E+8; gammam = PT/ 2.99792458E+8;
+      omega = q/e_SI * 2.99792458E+8/R_max;
+      r = PT / (omega * gammam );
+    }
+    // test -- point 6)
+    if(test && true) {
+      q = e_SI; Y = -R_max/2.; X = Z = Py = Pz = M = 0;
+      Px= -R_max * (q*B_z);
+      PT = sqrt(Px*Px + Py*Py + Pz*Pz);
+      E = PT* 2.99792458E+8; gammam = PT/ 2.99792458E+8;
+      omega = q/e_SI * 2.99792458E+8/R_max;
+      r = PT / (omega * gammam );
+    }
+
+
+double delta= UNDEFINED;
+
+if (method==1) {
+
+      phi_0 = -atan2(Px,Py);                // [rad]
+      //if(phi_0<0)phi_0 = 2*pi+phi_0;        // [rad], in [0 - 2 pi]
 
     // 2.  Helix parameters : center coordinates in transverse plane
     //     x_c = x_0 - r*cos(phi_0) and y_c = y_0 - r*sin(phi_0)
     //     R_c = \sqrt{x_c² + y_c²} and \Phi_c = atan{y_c/x_c}
-      x_c = Part->X - r*cos(phi_0);  /// TEST !!
-      y_c = Part->Y - r*sin(phi_0); 
+      x_c = X - r*cos(phi_0);
+      y_c = Y - r*sin(phi_0);
       R_c = sqrt(pow(x_c,2.) + pow(y_c,2.) );
       Phi_c = atan2(y_c,x_c);
+      //if(Phi_c<0)Phi_c = 2*pi+Phi_c;
+      Phi = Phi_c;
+      //r=fabs(r); 
 
     // 3. time evaluation t = min(t_T, t_z)
@@ -157 +256 @@
     //    t_z = gamma * m /p_z0 \times (-z_0 + z_max * sign(p_z0))
       rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable
-      t_T=0;
-      int sign_pz= (Part->Pz >0) ? 1 : -1;
-      t_z = gammam / Part->Pz * (-Part->Z + z_max*sign_pz ) ;
+      t_T=0; //[ns]
+      int sign_pz= (Pz >0) ? 1 : -1;
+      if(Pz==0) t_z = 1E99;
+      else t_z = gammam / (Pz*1E-9* 2.99792458E+8) * (-Z + z_max*sign_pz );
+      if( t_z <0) cout << "ERROR: t_z <0 !" << endl;
+
       if ( fabs(R_c - r) > R_max || R_c + r  < R_max ) t = t_z;
       else {
          if(r==0|| R_c ==0) t_T=1E99;
-         //else t_T = (Phi_c - phi_0 + atan2( (R_max + rr)*(R_max - rr) , 2*r*R_c ) ) / omega; 
-         else t_T = (Phi_c - phi_0 + acos( (R_max + rr)*(R_max - rr) / (2*r*R_c) ) ) / omega; 
-         t = min(t_T,t_z);
+         else { 
+                //t_T = fabs((Phi_c - phi_0 - acos( (R_max + rr)*(R_max - rr) / (2*r*R_c) ))/omega) ;
+                double A = fabs(acos( (R_max + rr)*(R_max - rr) / (2*r*R_c) )) ;
+                double a = phi_0 + A;
+       // validé si acos >0 , mais quid si acos < 0?
+
+                double tm = (Phi_c - a) / omega;
+                double tp = (-Phi_c + a) / omega;
+                //cout << "t- = " << tm << "\t t+ = " << tp << endl;
+                if(tm<0) t_T = tp;
+                else if(tp<0) t_T=tm;
+                else t_T = min(tm,tp);
+
+                //cout << "t_T = " << t_T << "\t T_z = " << t_z << "\t r = " << r << endl;
+                t = min(t_T,t_z); // t is output here in [ns], which is compatible with omega
+         }
       }
 
     // 4. position in terms of x(t), y(t), z(t)
-    //    x(t) = x_c + r cos (omega t + phi_0)
-    //    y(t) = y_c + r sin (omega t + phi_0)
-    //    z(t) = z_0 + (p_Z0/gammam) t
       x_t = x_c + r * cos(omega * t + phi_0); 
       y_t = y_c + r * sin(omega * t + phi_0); 
-      z_t = Part->Z + Part->Pz / gammam * t;
+      z_t = Z + Pz*1E-9* 2.99792458E+8 / gammam * t;
 
     // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t)
-    //    R(t) = sqrt(x(t)² + y(t)²)
-    //    Phi(t) = atan(y(t)/x(t))
-    //    Theta(t) = atan(R(t)/z(t))
-    //    Eta(t) = -ln tan (Theta(t)/2)
       R_t   = sqrt( pow(x_t,2.) + pow(y_t,2.)  );
       Phi_t = atan2( y_t, x_t);
@@ -187 +295 @@
               Eta_t = - log(tan(Theta_t/2.));
       } else{
-                Theta_t=0; Eta_t = 9999;
+                Theta_t=0; Eta_t = UNDEFINED;
       }
 
-        Px_t = - Part->PT * sin(omega*t + phi_0);
-        Py_t =   Part->PT * cos(omega*t + phi_0);
-        Pz_t =   Part->Pz;
+      //if(phi_0 <0 || phi_0 > 2*pi) cout <<"ERROR: phi_0 out of range: [0 ; 2pi] " << phi_0 << endl;
+      //if(Phi_c <0 || Phi_c > 2*pi) cout <<"ERROR: Phi_c out of range: [0 ; 2pi] " << Phi_c << endl;
+
+} // method 1
+else {
+
+      phi_0 = atan2(Py,Px); // [rad] in [-pi ; pi ]
+      //if(phi_0<0)phi_0 = 2*pi+phi_0;        // [rad], in [0 - 2 pi]
+
+      x_c = X + r*sin(phi_0);
+      y_c = Y - r*cos(phi_0);
+      R_c = sqrt( pow(x_c,2.) + pow(y_c,2.) );
+      Phi_c = atan2(y_c,x_c);
+      Phi = Phi_c;
+      if(x_c<0) Phi += pi;
+      //if(Phi<0)Phi = 2*pi+Phi; // ne pas le mettre pour que le test1 fonctionne
+
+    // 3. time evaluation t = min(t_T, t_z)
+    //    t_T : time to exit from the sides
+    //    t_z : time to exit from the front or the back
+      rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable
+      t_T=0; //[ns]
+      int sign_pz= (Pz >0) ? 1 : -1;
+      if(Pz==0) t_z = 1E99;
+      else t_z = gammam / (Pz*1E-9* 2.99792458E+8) * (-Z + z_max*sign_pz );
+      if( t_z <0) cout << "ERROR: t_z <0 !" << endl;
+
+      if ( fabs(R_c - fabs(r)) > R_max || R_c + fabs(r)  < R_max ) t = t_z;
+      else {
+         if(r==0) cout << "r ==0 !" << endl;
+         if(R_c==0) cout << "R_c ==0 !" << endl;
+         if(r==0|| R_c ==0) t_T=1E99;
+         else { 
+                double asinrho = asin( (R_max + rr)*(R_max - rr) / (2*fabs(r)*R_c)  );
+                delta = phi_0 - Phi;
+                if(delta<-pi) delta += 2*pi;
+                if(delta> pi) delta -= 2*pi;
+                double t1 = (delta + asinrho) / omega;
+                double t2 = (delta + pi - asinrho) / omega;
+                double t3 = (delta + pi + asinrho) / omega;
+                double t4 = (delta - asinrho) / omega;
+                double t5 = (delta - pi - asinrho) / omega;
+                double t6 = (delta - pi + asinrho) / omega;
+
+                if(test) {
+                        cout << "t4 = " << t4 << "\t t5 = " << t5 << "\t t_6=" << t6 << "\t t_3 = " << t3 << endl;
+                        cout << "t1 = " << t1 << "\t t2 = " << t2 << "\t t_T=" << t_T << "\t t_z = " << t_z << endl;
+                        cout << "delta= " << delta << endl;
+                }
+                if(t1<0)t1=1E99;
+                if(t2<0)t2=1E99;
+                if(t3<0)t3=1E99;
+                if(t4<0)t4=1E99;
+                if(t5<0)t5=1E99;
+                if(t6<0)t6=1E99;
+
+                double t_Ta = min(t1,min(t2,t3));
+                double t_Tb = min(t4,min(t5,t6));
+                t_T = min(t_Ta,t_Tb);
+
+                //if(t1<0) t_T = t2;
+                //else if(t2<0) t_T = t1;
+                //else {t_T = min(t1,t2); /*cout << "**";*/}
+                //if (t1<0 && t2<0) {t_T = fabs(min(t1,t2)); cout << "\tbad!!\n";}
+                t = min(t_T,t_z); 
+                if(test) {
+                  // cout << "t4 = " << t4 << "\t t5 = " << t5 << "\t t_6=" << t6 << "\t t_3 = " << t3 << endl;
+                  // cout << "t1 = " << t1 << "\t t2 = " << t2 << "\t t_T=" << t_T << "\t t_z = " << t_z << endl;
+                }
+         }
+      }
+//r = PT / (omega * gammam );
+
+    // 4. position in terms of x(t), y(t), z(t)
+      x_t = x_c + r * sin(omega * t - phi_0);
+      y_t = y_c + r * cos(omega * t - phi_0);
+      z_t = Z + Pz*1E-9* 2.99792458E+8 / gammam * t;
+
+    // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t)
+      R_t   = sqrt( pow(x_t,2.) + pow(y_t,2.)  );
+      Phi_t = atan2( y_t, x_t);
+      if(R_t>0) {
+                Theta_t = acos( z_t / sqrt(z_t*z_t+ R_t*R_t));
+                Eta_t = - log(tan(Theta_t/2.));
+      }
+      else {
+                Theta_t=0; Eta_t = UNDEFINED;
+      }
+} // method2
+
+     if(test) {
+        cout << endl << endl;
+        cout << "method " << method << "----------------\n";
+        cout << "x0,y0,z0= " << X << ", " << Y << ", " << Z << endl;
+        cout << "px0,py0,pz0= " << Px << ", " << Py << ", " << Pz << endl;
+        cout << "r = " << r << "R_max = " << R_max << "\t phi_0=" << phi_0 << endl;
+        cout << "gammam= " << gammam << "\t omega=" << omega << "\t PT = " << PT << endl;
+        cout << "x_c = " << x_c << "\t y_c = " << y_c << "\t R_c = " << R_c << "\t Phi = " << Phi << endl;
+        cout << "omega t = " << omega*t << "\t";
+        cout << "cos(omega t -phi0)= " << cos(omega*t-phi_0) << "\t sin(omega t -phi0)= " << sin(omega*t-phi_0) << endl;
+        cout << "t_T = " << t_T << "\t t_z = " << t_z << "\t r = " << r << endl;
+        cout << "x_t = " << x_t << "\t y_t = " << y_t << "\t z_t = " << z_t << endl;
+        cout << "R_t = " << R_t << "\t Phi_t = " << Phi_t << "\t";
+        cout << "Theta_t = " << Theta_t << "\t Eta_t = " << Eta_t << endl;
+     }
+
+      
+/*      Not needed here. but these formulae are correct -------
+        // method1
+        Px_t = - PT * sin(omega*t + phi_0); 
+        Py_t =   PT * cos(omega*t + phi_0);
+
+        // method2
+        Px_t =   PT * cos(phi_0 - omega*t);  
+        Py_t =   PT * sin(phi_0 - omega*t);
+
+        Pz_t =   Pz;
         PT_t =   sqrt(Px_t*Px_t + Py_t*Py_t);
         p_t = sqrt(PT_t*PT_t + Pz_t*Pz_t);
-        E_t=sqrt(Part->M*Part->M +p_t);
+        E_t=sqrt(M*M +p_t*p_t);
         //if(p_t != fabs(Pz_t) ) Eta_t = log( (p_t+Pz_t)/(p_t-Pz_t) )/2.;
-        //if(p_t>0) Theta_t = acos(Pz_t/p_t);
+        //if(p_t>0) Theta_t = acos(Pz_t/p_t)>;
         momentum.SetPxPyPzE(Px_t,Py_t,Pz_t,E_t);
+*/
+
+//cout << "R_c = " << R_c << " r " << r << " rr = " << rr << " R_t" << R_t << endl;
+//cout << "x_t = " << x_t << " x_c = " << x_c <<  "\ty_t = " << y_t << " y_c = " << y_c << " z_t " << z_t << endl;
+//cout << "Eta = " << Part->Eta << " Eta_t " << Eta_t << "Phi = " << Part->Phi << " Phi_t= " << Phi_t << "-----" << endl;
+        Part->EtaCalo = Eta_t;
+        Part->PhiCalo = Phi_t;
 
 // test zone ---
 /*
+
+        cout << "r = " << r << " et " << fabs(PT/(q*B_z)) << endl;
         cout << cos(atan(R_t/z_t)) << "\t" << cos(Theta_t) << "\t" << cos(momentum.Theta()) << "\t" << Pz_t/temp_p << endl;
         double Eta_t1 = log( (E+Pz_t)/(E-Pz_t) )/2.;
@@ -232 +463 @@
    //cout << "tan(phi_p)= " << momentum.Py()/momentum.Px() << "\t -1/tan(phi_x)= " << -x_t/y_t << endl;
 */
+        return;
 
   } else { // if B_x or B_y are non zero: longer computation
-
+//cout << "bfield de loic\n";
   float Xvertex1 = Part->X;
   float Yvertex1 = Part->Y;
   float Zvertex1 = Part->Z;
-  
-  //out of tracking coverage?
-  if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > R_max){return;}
-  if(fabs(Zvertex1) > z_max){return;}
   
   double px = Part->Px / 0.003;
@@ -254 +482 @@
   double vz = pz/M;
   double qm = q/M;
-
-  double ax =  qm*(B_z*vy - B_y*vz);
-  double ay =  qm*(B_x*vz - B_z*vx);
-  double az =  qm*(B_y*vx - B_x*vy);
-  double dt = 1/p;
-  if(pt<266 && vz < 0.0012)       dt = fabs(0.001/vz); // ?????
  
-  double xold=Xvertex1;         double x=xold;
-  double yold=Yvertex1;         double y=yold;
-  double zold=Zvertex1;         double z=zold;
-
-  double VTold = pt/M; //=sqrt(vx*vx+vy*vy);
- 
-  unsigned int k = 0;
-  double VTratio=0;
-  double R_max2 = R_max*R_max;
-  double r2=0; // will be x*x+y*y
-
-  while(k < MAXITERATION){
-        k++;
-
-        vx += ax*dt;
-        vy += ay*dt;
-        vz += az*dt;
-
-        VTratio = VTold/sqrt(vx*vx+vy*vy);
-        vx *= VTratio;
-        vy *= VTratio;
-
-        ax = qm*(B_z*vy - B_y*vz);
-        ay = qm*(B_x*vz - B_z*vx);
-        az = qm*(B_y*vx - B_x*vy);
-
-        x  += vx*dt;
-        y  += vy*dt;
-        z  += vz*dt;
-        r2  = x*x + y*y;
-
-       if( r2 > R_max2 ){ 
-                x /= r2/R_max2; 
-                y /= r2/R_max2; 
-                break;
-       }
-       if( fabs(z)>z_max)break;
-
-       xold = x;
-       yold = y;
-       zold = z;
-  } // while loop
-
-  if(k == MAXITERATION) loop_overflow_counter++;
-  //cout << "too short loop in " << loop_overflow_counter << " cases" << endl;
- 
-  if(x!=0 && y!=0 && z!=0) {
-          float Theta = atan2(sqrt(r2),z);
-          double eta  = -log(tan(Theta/2.));
-          double phi = atan2(y,x);
-          momentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
-  }
-    
- } // if b_x or b_y non zero
-}
-
-
-
-void TrackPropagation::bfield(TRootGenParticle *Part) {
-
-  // initialisation, valid for z_max==0, R_max==0 and q==0
-  Part->EtaCalo = Part->Eta;
-  Part->PhiCalo = Part->Phi;//-atan2(Part->Px,Part->Py);
-
-  if (!DET->FLAG_bfield ) return;
-
-  q  = ChargeVal(Part->PID);
-  if(q==0) return;
-  if(R_max ==0) { cout << "ERROR: magnetic field has no lateral extention\n"; return;}
-  if(z_max==0) { cout << "ERROR: magnetic field has no longitudinal extention\n"; return;}
-
-  if (B_x== 0 && B_y== 0) { // faster if only B_z
-    if (B_z==0) return; // nothing to do
-
-    // initial conditions:
-      // p_X0 = Part->Px, p_Y0 = Part->Py, p_Z0 = Part->Pz, p_T0 = Part->PT;
-      // X_0 = Part->X, Y_0 = Part->Y, Z_0 = Part->Z;
-
-    // 1.  initial transverse momentum p_{T0} : Part->PT
-    //     initial transverse momentum direction \phi_0 = -atan(p_X0/p_Y0) 
-    //     relativistic gamma : gamma = E/mc² ; gammam = gamma \times m
-    //     giration frequency \omega = q/(gamma m) B_z
-    //     helix radius r = p_T0 / (omega gamma m)
-      phi_0 = -atan2(Part->Px,Part->Py);
-      gammam = Part->E; // here c==1
-      //cout << "gammam" << gammam << "\t gamma" << gammam/Part->M << endl;
-      omega = q * B_z /gammam;
-      r = fabs(Part->PT / (omega * gammam));
-
-    // 2.  Helix parameters : center coordinates in transverse plane
-    //     x_c = x_0 - r*cos(phi_0) and y_c = y_0 - r*sin(phi_0)
-    //     R_c = \sqrt{x_c² + y_c²} and \Phi_c = atan{y_c/x_c}
-      x_c = Part->X - r*cos(phi_0);  /// TEST !!
-      y_c = Part->Y - r*sin(phi_0); 
-      R_c = sqrt(pow(x_c,2.) + pow(y_c,2.) );
-      Phi_c = atan2(y_c,x_c);
-
-    // 3. time evaluation t = min(t_T, t_z)
-    //    t_T : time to exit from the sides
-    //    t_T= [ Phi_c - phi_0 + acos( (R_max^2 - (R_c^2 + r^2))/(2rR_c) )  ]/omega
-    //    t_z : time to exit from the front or the back
-    //    t_z = gamma * m /p_z0 \times (-z_0 + z_max * sign(p_z0))
-      rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable
-      t_T=0;
-      int sign_pz= (Part->Pz >0) ? 1 : -1;
-      t_z = gammam / Part->Pz * (-Part->Z + z_max*sign_pz ) ;
-      if ( fabs(R_c - r) > R_max || R_c + r  < R_max ) t = t_z;
-      else {
-         if(r==0|| R_c ==0) t_T=1E99;
-         //else t_T = (Phi_c - phi_0 + atan2( (R_max + rr)*(R_max - rr) , 2*r*R_c ) ) / omega;
-         else t_T = (Phi_c - phi_0 + acos( (R_max + rr)*(R_max - rr) / (2*r*R_c) ) ) / omega;
-         t = min(t_T,t_z);
-      }
-
-    // 4. position in terms of x(t), y(t), z(t)
-    //    x(t) = x_c + r cos (omega t + phi_0)
-    //    y(t) = y_c + r sin (omega t + phi_0)
-    //    z(t) = z_0 + (p_Z0/gammam) t
-      x_t = x_c + r * cos(omega * t + phi_0); 
-      y_t = y_c + r * sin(omega * t + phi_0); 
-      z_t = Part->Z + Part->Pz / gammam * t;
-
-    // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t)
-    //    R(t) = sqrt(x(t)² + y(t)²)
-    //    Phi(t) = atan(y(t)/x(t))
-    //    Theta(t) = atan(R(t)/z(t))
-    //    Eta(t) = -ln tan (Theta(t)/2)
-      R_t   = sqrt( pow(x_t,2.) + pow(y_t,2.)  );
-      Phi_t = atan2( y_t, x_t);
-      if(R_t>0) {
-              Theta_t = acos( z_t / sqrt(z_t*z_t+ R_t*R_t));
-              Eta_t = - log(tan(Theta_t/2.));
-      } else{
-                Theta_t=0; Eta_t = UNDEFINED;
-      }
-/*      Not needed here. but these formulae are correct -------
-        Px_t = - Part->PT * sin(omega*t + phi_0);
-        Py_t =   Part->PT * cos(omega*t + phi_0);
-        Pz_t =   Part->Pz;
-        PT_t =   sqrt(Px_t*Px_t + Py_t*Py_t);
-        p_t = sqrt(PT_t*PT_t + Pz_t*Pz_t);
-        E_t=sqrt(Part->M*Part->M +p_t);
-        //if(p_t != fabs(Pz_t) ) Eta_t = log( (p_t+Pz_t)/(p_t-Pz_t) )/2.;
-        //if(p_t>0) Theta_t = acos(Pz_t/p_t);
-        momentum.SetPxPyPzE(Px_t,Py_t,Pz_t,E_t);
-*/
-        Part->EtaCalo = Eta_t;
-        Part->PhiCalo = Phi_t;
-        return;
-// test zone ---
-/*
-        cout << cos(atan(R_t/z_t)) << "\t" << cos(Theta_t) << "\t" << cos(momentum.Theta()) << "\t" << Pz_t/temp_p << endl;
-        double Eta_t1 = log( (E+Pz_t)/(E-Pz_t) )/2.;
-        double Eta_t2 = log( (temp_p+Pz_t)/(temp_p-Pz_t) )/2.;
-        if(0 &&  fabs(Eta_t -Eta_t2)>1e-310) {
-                cout << "ERROR-BUG: Eta_t != Eta_t2\n";
-                cout << "Eta_t= " << Eta_t << "\t Eta_t1= " << Eta_t1 << "\t Eta_t2= " << Eta_t2 << endl;
-        }
-
-        double R_t2  = sqrt( pow(R_c,2.) + pow(r,2.) + 2*r*R_c*cos(phi_0 + omega*t - Phi_c) ); // cross-check
-        if(fabs(R_t - R_t2) > 1e-7) 
-                cout << "ERROR-BUG: R_t != R_t2:  R_t=" << R_t << " R_t2=" << R_t2 << " R_t - R_t2 =" << R_t - R_t2 << endl;
-        if( fabs(E - gammam) > 1e-3 ) {
-                cout << "ERROR-BUG: energy is not conserved in src/BFieldProp.cc\n";
-                cout << "E - momentum.E() = " <<  fabs(E - momentum.E()) <<  " gammam - E " << fabs(gammam -E) << endl; }
-        if( fabs(PT_t - Part->PT) > 1e-10 ) {
-                cout << "ERROR-BUG: PT is not conversed in src/BFieldProp.cc. ";
-                cout << "(at " << 100*(PT_t - Part->PT) << "%)\n";
-                }
-        if(momentum.Pz() != Pz_t)
-                cout << "ERROR-BUG: Pz is not conserved in src/BFieldProp.cc\n";
-
-        double temp_p0=sqrt(Part->PT*Part->PT + Part->Pz*Part->Pz);
-        if(fabs( (temp_p-temp_p0)*(temp_p+temp_p0) )>1e-10  ) {
-                cout << "ERROR-BUG: momentum |vec{p}| is not conserved in src/BFieldProp.cc\n";
-                cout << temp_p << "\t" << temp_p0 << endl;
-        }
-
-   // if x_c == y_c ==0 (set it by hand!), easy cross-check
-   //cout << "tan(phi_p)= " << momentum.Py()/momentum.Px() << "\t -1/tan(phi_x)= " << -x_t/y_t << endl;
-*/
-
-  } else { // if B_x or B_y are non zero: longer computation
-//cout << "bfield de loic\n";
-  float Xvertex1 = Part->X;
-  float Yvertex1 = Part->Y;
-  float Zvertex1 = Part->Z;
-  
-  //out of tracking coverage?
-  if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > R_max){return;}
-  if(fabs(Zvertex1) > z_max){return;}
-  
-  double px = Part->Px / 0.003;
-  double py = Part->Py / 0.003;
-  double pz = Part->Pz / 0.003;
-  double pt = Part->PT / 0.003; // sqrt(px*px+py*py);
-  double p  = sqrt(pz*pz + pt*pt); //sqrt(px*px+py*py+pz*pz);
- 
-  double M  = Part->M;
-  double vx = px/M;
-  double vy = py/M;
-  double vz = pz/M;
-  double qm = q/M;
+//double v = sqrt(vx*vx + vy*vy + vz*vz)/3E8;
+//cout << "v = " << v;
+//double gamma = 1./sqrt(1-v*v);
+//cout << "gamma = " << gamma << endl;
 
   double ax =  qm*(B_z*vy - B_y*vz);