Changeset 193 in svn

Timestamp:

Jan 26, 2009, 3:43:42 PM (16 years ago)

Author:

Xavier Rouby

Message:

new algorithm inserted but not yet active

Location:

trunk

Files:

• interface/BFieldProp.h (modified) (2 diffs)
• src/BFieldProp.cc (modified) (3 diffs)

trunk/interface/BFieldProp.h

@@ -12 +12 @@
  *             */
 
-#include <vector>
 #include "TLorentzVector.h"
-
 #include "interface/SmearUtil.h"
-
 #include "Utilities/ExRootAnalysis/interface/BlockClasses.h"
 #include "Utilities/ExRootAnalysis/interface/TSimpleArray.h"
@@ -26 +23 @@
  public:
   // Constructor
-    TrackPropagation(string DetDatacard);
+    TrackPropagation(const string DetDatacard);
  
     void Propagation(const TRootGenParticle *Part,TLorentzVector &genMomentum);
 
-    int MAXITERATION;
-    int MINSEGLENGTH;
 
  private:
+   unsigned int MAXITERATION;
+   RESOLution *DET;
 
-   RESOLution *DET;  
+
+   /// radial/longitudinal extensions of magnetic field volume
+   double R_max, z_max;
+   /// magnetic field components
+   double B_x, B_y, B_z;
+   /// particle charge
+   double q;
+   /// initial coordinate
+   double phi_0;
+   /// relativistic gamma \times m
+   double gammam;
+   /// giration frequency and radius
+   double omega, r;
+   /// center of the helix in the transverse plane
+   double x_c, y_c, R_c, Phi_c;
+   // variable for an intermediate computing
+   double rr;
+   /// times for exiting the magnetic field volume
+   double t, t_z, t_T;
+   /// coordinates of the exit point
+   double x_t, y_t, z_t, R_t, Phi_t, Theta_t, Eta_t;
+
+   unsigned int loop_overflow_counter;
 };
   

trunk/src/BFieldProp.cc

@@ -10 +10 @@
 
 #include "interface/BFieldProp.h"
-#include "TRandom.h"
-
-#include <iostream>
-#include <sstream>
-#include <fstream>
-#include <iomanip>
-
 #include<cmath>
-
-
+#include "TMath.h"
 using namespace std;
 
@@ -25 +17 @@
 //------------------------------------------------------------------------------
 
-TrackPropagation::TrackPropagation(string DetDatacard) {
-
+TrackPropagation::TrackPropagation(const string DetDatacard):
+ MAXITERATION(10000), q(-9999.), phi_0(-9999.), gammam(-9999.), omega(-9999.), r(-9999.),
+ x_c(-9999.), y_c(-9999.), R_c(-9999.), Phi_c(-9999.),
+ 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.) {
+
+ // if(DetDatacard="") { DET = new RESOLution(); }
+ // else DET = new RESOLution(DetDatacard);
  DET = new RESOLution();
  DET->ReadDataCard(DetDatacard);
- MAXITERATION = 10000;
- MINSEGLENGTH = 70;
-
+
+ // 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;
+   B_z = DET->TRACK_bfield_z;
+
+   loop_overflow_counter=0;
 }
 
-void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &genMomentum)
-{
-
-  double q  = Charge(Part->PID);
-  if(q==0)return;
-  
+void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &momentum) {
+
+  q  = Charge(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 (0 && 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 = 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); 
+      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_z : time to exit from the front or the back
+      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) ) ;
+      if ( fabs(R_c - r) > R_max || R_c + r  < R_max ) t = t_z;
+      else {
+         t_T = (Phi_c - phi_0 + atan2( (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.)  );
+      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) 
+                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 ) {
+                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 << "(at " << 100*(PT_t - Part->PT) << "%)\n";
+                }
+        if(momentum.Pz() != Pz_t)
+                cout << "ERROR-BUG: Pz is not conserved in src/BFieldProp.cc\n";
+
+      //momentum.SetPtEtaPhiE(Part->PT,Eta_t,Phi_t,Part->E);
+    
+  } else { // if B_x or B_y are non zero: longer computation
+
   float Xvertex1 = Part->X;
   float Yvertex1 = Part->Y;
   float Zvertex1 = Part->Z;
   
-  //out of trackibg coverage?
-  if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > DET->TRACK_radius){return;}
-  if(fabs(Zvertex1) > DET->TRACK_length/2.){return;}
+  //out of tracking coverage?
+  if(sqrt(Xvertex1*Xvertex1+Yvertex1*Yvertex1) > R_max){return;}
+  if(fabs(Zvertex1) > z_max){return;}
   
-  float Px = Part->Px;
-  float Py = Part->Py;
-  float Pz = Part->Pz;
-    
-  double px = Px / 0.003;
-  double py = Py / 0.003;
-  double pz = Pz / 0.003;
-  double pt = sqrt(px*px+py*py);
+  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(px*px+py*py+pz*pz);
  
-  if(q!=0){
-     double M  = Part->M;
-     double vx = px/M;
-     double vy = py/M;
-     double vz = pz/M;
-
-     double Bx = DET->TRACK_bfield_x;
-     double By = DET->TRACK_bfield_y;
-     double Bz = DET->TRACK_bfield_z;
-
-     double ax =  (q/M)*(Bz*vy - By*vz);
-     double ay =  (q/M)*(Bx*vz - Bz*vx);
-     double az =  (q/M)*(By*vx - Bx*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 = sqrt(vx*vx+vy*vy);
- 
-     int k = 0;
-     
-     double Radius=DET->TRACK_radius;
-     double Length=DET->TRACK_length/2.;
-
-     while(k < MAXITERATION){
+  double M  = Part->M;
+  double vx = px/M;
+  double vy = py/M;
+  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++;
 
@@ -92 +185 @@
         vz += az*dt;
 
-        double VTratio = VTold/sqrt(vx*vx+vy*vy);
+        VTratio = VTold/sqrt(vx*vx+vy*vy);
         vx *= VTratio;
         vy *= VTratio;
 
-        ax = (q/M)*(Bz*vy - By*vz);
-        ay = (q/M)*(Bx*vz - Bz*vx);
-        az = (q/M)*(By*vx - Bx*vy);
+        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;
-
-       if( (x*x+y*y) > Radius*Radius ){ x /= (x*x+y*y)/(Radius*Radius); y /= (x*x+y*y)/(Radius*Radius); break;}
-       if( fabs(z)>Length)break;
+        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;
-     }
- 
-      if(x!=0 && y!=0 && z!=0)
-        {
-          double eta;
-          float Theta = atan2(sqrt(x*x+y*y),z);
-          eta  = -log(tan(Theta/2));
+  } // 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);
-          genMomentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
-        }
-    }
+          momentum.SetPtEtaPhiE(Part->PT,eta,phi,Part->E);
+  }
+    
+ } // if b_x or b_y non zero
 }