Changeset 434 in svn for trunk/src/BFieldProp.cc

Timestamp:

Jun 17, 2009, 10:34:40 PM (15 years ago)

Author:

Xavier Rouby

Message:

new Bfield bug free

File:

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

trunk/src/BFieldProp.cc

@@ -109 +109 @@
   // magnetic field parameters
    R_max = DET->TRACK_radius/100.; //[m]
-   z_max = DET->TRACK_length/200.; //[m]
+   z_max = DET->TRACK_length/100.; //[m]
    B_x = DET->TRACK_bfield_x*tesla;
    B_y = DET->TRACK_bfield_y*tesla;
@@ -132 +132 @@
   if (!DET->FLAG_bfield ) return;
 
-  double M;
+  double M; // GeV/c²
   //int q1  = ChargeVal(Part->PID) *eplus; // in units of 'e'
   if(Part->M < -999) { // unitialised!
      PdgParticle pdg_part = DET->PDGtable[Part->PID];
      q  = pdg_part.charge() *eplus; // in units of 'e'
-     M  = pdg_part.mass(); // GeV
+     M  = pdg_part.mass(); // GeV/c²
   } else  { q = Part->Charge; M = Part->M; }
 
@@ -154 +154 @@
   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
@@ -168 +167 @@
     //     giration frequency \omega = q/(gamma m) B_z
     //     helix radius r = p_T0 / (omega gamma m)
-    double Px = Part->Px; // [GeV/c]
+    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 = 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)
-    //    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; //[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 = 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 + 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;
-      }
-
-      //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]
-
+    double E  = Part->E;                     // [GeV]
+    double Phi = UNDEFINED;
+
+    float c_light = 2.99792458E+8;
+    gammam = E*1E9/(c_light*c_light);               // gammam in [eV/c²]
+    omega = q * B_z / (gammam);                     // omega is here in [ 89875518 / s]  
+    //cout << "omega*gammam = B_z in BFieldProp.cc: " << fabs(omega*gammam) - B_z  << endl;
+    r = PT / (omega * gammam)  *1E9/c_light;        // in [m]  
+
+    // test mode ?
+    bool test=false; if(test) loop_overflow_counter++;
+
+    double delta= UNDEFINED;
+    phi_0 = atan2(Py,Px); // [rad] in [-pi ; pi ]
+
+    // 2. helix axis coordinates
       x_c = X + r*sin(phi_0);
       y_c = Y - r*cos(phi_0);
@@ -321 +193 @@
       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]
+      rr = sqrt( pow(R_c,2.) + pow(r,2.) );     // temp variable [m]
+      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 );
+      else t_z = gammam / (Pz*1E9/c_light) * (-Z + z_max*sign_pz );
       if( t_z <0) cout << "ERROR: t_z <0 !" << endl;
 
@@ -365 +236 @@
                 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;
+      z_t = Z + Pz*1E9/c_light / gammam * t;
 
     // 5. position in terms of Theta(t), Phi(t), R(t), Eta(t)
@@ -392 +253 @@
       }
       else {
-                Theta_t=0; Eta_t = UNDEFINED;
+                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;
@@ -414 +274 @@
       
 /*      Not needed here. but these formulae are correct -------
-        // method1
+        // method1 (removed)
         Px_t = - PT * sin(omega*t + phi_0); 
         Py_t =   PT * cos(omega*t + phi_0);
@@ -431 +291 @@
 */
 
-//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;