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

Timestamp:

Feb 2, 2009, 12:33:21 PM (16 years ago)

Author:

Xavier Rouby

Message:

JetUtils.cc

File:

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

trunk/src/BFieldProp.cc

@@ -9 +9 @@
  *             */
 
-#include "interface/BFieldProp.h"
+#include "BFieldProp.h"
 #include<cmath>
-#include "TMath.h"
 using namespace std;
 
@@ -17 +16 @@
 //------------------------------------------------------------------------------
 
-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.),
- Px_t(-9999), Py_t(-9999), Pz_t(-9999), PT_t(-9999), p_t(-9999), E_t(-9999) {
-
- // if(DetDatacard="") { DET = new RESOLution(); }
- // else DET = new RESOLution(DetDatacard);
- DET = new RESOLution();
- DET->ReadDataCard(DetDatacard);
-
- // magnetic field parameters
+TrackPropagation::TrackPropagation(){
+   DET = new RESOLution();
+   init();
+}
+
+TrackPropagation::TrackPropagation(const string& DetDatacard){
+   DET = new RESOLution();
+   DET->ReadDataCard(DetDatacard);
+   init();      
+}
+
+TrackPropagation::TrackPropagation(const RESOLution* DetDatacard){ 
+   DET= new RESOLution(*DetDatacard);
+   init();
+}
+
+TrackPropagation::TrackPropagation(const TrackPropagation & tp){
+   MAXITERATION = tp.MAXITERATION;
+   DET   = new RESOLution(*(tp.DET));
+   R_max = tp.R_max;    z_max = tp.z_max;
+   B_x   = tp.B_x;      B_y   = tp.B_y;         B_z   = tp.B_z;
+   q     = tp.q;        phi_0 = tp.phi_0;       
+   gammam= tp.gammam;   omega = tp.omega;
+   r     = tp.r;        rr    = tp.rr;
+   x_c   = tp.x_c;      y_c   = tp.y_c; 
+   R_c   = tp.R_c;      Phi_c = tp.Phi_c;
+   t     = tp.t;        t_z   = tp.t_z;         t_T   = tp.t_T; 
+   x_t   = tp.x_t;      y_t   = tp.y_t;         z_t   = tp.z_t;
+   R_t   = tp.R_t;      Phi_t = tp.Phi_t;       
+   Theta_t=tp.Theta_t;  Eta_t = tp.Eta_t;
+   Px_t  = tp.Px_t;     Py_t  = tp.Py_t;        Pz_t  = tp.Pz_t;
+   PT_t  = tp.PT_t;     p_t   = tp.p_t;         E_t   = tp.E_t;
+   loop_overflow_counter = tp.loop_overflow_counter;
+}
+
+TrackPropagation& TrackPropagation::operator=(const TrackPropagation & tp) {
+   if(this==&tp) return *this;
+   MAXITERATION = tp.MAXITERATION;
+   DET   = new RESOLution(*(tp.DET));
+   R_max = tp.R_max;    z_max = tp.z_max;
+   B_x   = tp.B_x;      B_y   = tp.B_y;         B_z   = tp.B_z;
+   q     = tp.q;        phi_0 = tp.phi_0;
+   gammam= tp.gammam;   omega = tp.omega;
+   r     = tp.r;        rr    = tp.rr;
+   x_c   = tp.x_c;      y_c   = tp.y_c;
+   R_c   = tp.R_c;      Phi_c = tp.Phi_c;
+   t     = tp.t;        t_z   = tp.t_z;         t_T   = tp.t_T;
+   x_t   = tp.x_t;      y_t   = tp.y_t;         z_t   = tp.z_t;
+   R_t   = tp.R_t;      Phi_t = tp.Phi_t;
+   Theta_t=tp.Theta_t;  Eta_t = tp.Eta_t;
+   Px_t  = tp.Px_t;     Py_t  = tp.Py_t;        Pz_t  = tp.Pz_t;
+   PT_t  = tp.PT_t;     p_t   = tp.p_t;         E_t   = tp.E_t;
+   loop_overflow_counter = tp.loop_overflow_counter;
+   return *this;
+}
+
+
+
+void TrackPropagation::init() {
+ MAXITERATION = 10000; 
+ q= UNDEFINED;  phi_0= UNDEFINED; gammam= UNDEFINED; omega=UNDEFINED; r=UNDEFINED;
+ x_c=UNDEFINED; y_c=UNDEFINED;    R_c=UNDEFINED;     Phi_c=UNDEFINED;
+ rr=UNDEFINED;  t=UNDEFINED;      t_z=UNDEFINED;     t_T=UNDEFINED;
+ x_t=UNDEFINED; y_t=UNDEFINED;    z_t=UNDEFINED;
+ R_t=UNDEFINED; Phi_t=UNDEFINED;  Theta_t=UNDEFINED; Eta_t=UNDEFINED;
+ Px_t=UNDEFINED; Py_t=UNDEFINED;  Pz_t=UNDEFINED;    PT_t=UNDEFINED; p_t=UNDEFINED; E_t=UNDEFINED;
+
+  // DET has been initialised in the constructors
+  // magnetic field parameters
    R_max = DET->TRACK_radius;
    z_max = DET->TRACK_length/2.;
@@ -39 +93 @@
    loop_overflow_counter=0;
 }
+
+
 
 void TrackPropagation::Propagation(const TRootGenParticle *Part,TLorentzVector &momentum) {
@@ -81 +137 @@
       rr = sqrt( pow(R_c,2.) + pow(r,2.) ); // temp variable
       t_T=0;
-      t_z = gammam / Part->Pz * (-Part->Z + z_max* TMath::Sign((Float_t)1.,(Float_t)Part->Pz) ) ;
+      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 {
@@ -103 +160 @@
       R_t   = sqrt( pow(x_t,2.) + pow(y_t,2.)  );
       Phi_t = atan2( y_t, x_t);
-/*      if(R_t>0) {
+      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.));
@@ -109 +166 @@
                 Theta_t=0; Eta_t = 9999;
       }
-*/
+
         Px_t = - Part->PT * sin(omega*t + phi_0);
         Py_t =   Part->PT * cos(omega*t + phi_0);
@@ -116 +173 @@
         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);
+        //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);