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

Timestamp:

Jan 28, 2009, 12:03:03 AM (16 years ago)

Author:

Xavier Rouby

Message:

fast bfield is working

File:

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

trunk/src/BFieldProp.cc

@@ -22 +22 @@
  rr(-9999.), t(-9999.), t_z(-9999.), t_T(-9999.),
  x_t(-9999.), y_t(-9999.), z_t(-9999.), 
- R_t(-9999.), Phi_t(-9999.), Theta_t(-9999.), Eta_t(-9999.) {
+ R_t(-9999.), Phi_t(-9999.), Theta_t(-9999.), Eta_t(-9999.),
+ Px_t(-9999), Py_t(-9999), Pz_t(-9999), PT_t(-9999), p_t(-9999), E_t(-9999) {
 
  // if(DetDatacard="") { DET = new RESOLution(); }
@@ -47 +48 @@
   if(z_max==0) { cout << "ERROR: magnetic field has no longitudinal extention\n"; return;}
 
-  if (0 && B_x== 0 && B_y== 0) { // faster if only B_z
+  if (B_x== 0 && B_y== 0) { // faster if only B_z
     if (B_z==0) return; // nothing to do
 
@@ -68 +69 @@
     //     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); 
+      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.) );
@@ -75 +76 @@
     // 3. time evaluation t = min(t_T, t_z)
     //    t_T : time to exit from the sides
+    //    t_T= [ Phi_c - phi_0 + atan( (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;
-      t_z = gammam / Part->Pz * (-Part->Z + DET-> TRACK_length/2. * TMath::Sign((Float_t)1.,(Float_t)Part->Pz) ) ;
+      t_z = gammam / Part->Pz * (-Part->Z + z_max* TMath::Sign((Float_t)1.,(Float_t)Part->Pz) ) ;
       if ( fabs(R_c - r) > R_max || R_c + r  < R_max ) t = t_z;
       else {
@@ -99 +102 @@
     //    Eta(t) = -ln tan (Theta(t)/2)
       R_t   = sqrt( pow(x_t,2.) + pow(y_t,2.)  );
-      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-11) 
+      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 = 9999;
+      }
+*/
+        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);
+
+// 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;
-      Phi_t = atan2( y_t, x_t);
-      Theta_t = atan2( R_t , z_t);
-      Eta_t = - log(tan(Theta_t/2.));
-
-        double Px_t = - Part->PT * sin(omega*t + phi_0);
-        double Py_t =   Part->PT * cos(omega*t + phi_0);
-        double Pz_t =   Part->Pz;
-        double PT_t =   sqrt(Px_t*Px_t + Py_t*Py_t);
-        double E=sqrt(Part->M*Part->M + PT_t*PT_t + Pz_t*Pz_t); // cross-check
-        momentum.SetPxPyPzE(Px_t,Py_t,Pz_t,E);
-        TLorentzVector mom_temp;
-        mom_temp.SetPtEtaPhiM(PT_t,Eta_t,Phi_t+PI,Part->M);
-
-        if (momentum != mom_temp) { 
-//              cout << "momentum et mom_temp different\n";
-/*              if (momentum.Px() != mom_temp.Px() ) cout << "  px: " << momentum.Px() << " & " << mom_temp.Px() << endl;
-                if (momentum.Py() != mom_temp.Py() ) cout << "  py: " << momentum.Py() << " & " << mom_temp.Py() << endl;
-                if (momentum.Pz() != mom_temp.Pz() ) cout << "  pz: " << momentum.Pz() << " & " << mom_temp.Pz() << endl;
-                if (momentum.Pt() != mom_temp.Pt() ) cout << "  pt: " << momentum.Pt() << " & " << mom_temp.Pt() << endl;*/
-//              cout << "eta: " << momentum.Eta() << " " << mom_temp.Eta() << endl;
-//              cout << "the: " << momentum.Theta() << " " << mom_temp.Theta() << endl;
-//              cout << "phi: " << momentum.Phi() << " " << mom_temp.Phi() << endl;
-        }
-        //cout << Part->PT << " or " << momentum.Pt() << " or " << mom_temp.Pt() << endl;
-
-        if( fabs(E - gammam) > 1e-4 ) {
+        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 conserved in src/BFieldProp.cc. ";
+                cout << "ERROR-BUG: PT is not conversed in src/BFieldProp.cc. ";
                 cout << "(at " << 100*(PT_t - Part->PT) << "%)\n";
                 }
@@ -137 +143 @@
                 cout << "ERROR-BUG: Pz is not conserved in src/BFieldProp.cc\n";
 
-      //momentum.SetPtEtaPhiE(Part->PT,Eta_t,Phi_t,Part->E);
-    
+        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
 
@@ -153 +167 @@
   double pz = Part->Pz / 0.003;
   double pt = Part->PT / 0.003; // sqrt(px*px+py*py);
-  double p  = sqrt(px*px+py*py+pz*pz);
+  double p  = sqrt(pz*pz + pt*pt); //sqrt(px*px+py*py+pz*pz);
  
   double M  = Part->M;