move from Bz to BxByBz in track propagation

[u/mrichter/AliRoot.git] / STEER / AliHelix.cxx

diff --git a/STEER/AliHelix.cxx b/STEER/AliHelix.cxx
index 1810a0f68f3b412210a284bc4d7e5df39fee44c5..ce4e5aa9c27ba28be7e9631593ee5b4326e20079 100644 (file)
--- a/STEER/AliHelix.cxx
+++ b/STEER/AliHelix.cxx
@@ -23,7 +23,7 @@
 
 #include "AliHelix.h"
 #include "AliKalmanTrack.h"
-#include "AliExternalTrackParam.h"
+#include "AliTracker.h"
 #include "TMath.h"
 ClassImp(AliHelix)
 
@@ -54,13 +54,20 @@ AliHelix::AliHelix(const AliKalmanTrack &t)
   alpha=t.GetAlpha();
   //
   //circle parameters
-  fHelix[4]=fHelix[4]/t.GetConvConst();    // C
+  //PH Sometimes fP4 and fHelix[4] are very big and the calculation
+  //PH of the Sqrt cannot be done. To be investigated...
+  fHelix[4]=fHelix[4]/(-1000/0.299792458/AliTracker::GetBz());    // C
   cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
 
   Double_t xc, yc, rc;
   rc  =  1/fHelix[4];
   xc  =  x-fHelix[2]*rc;
-  yc  =  fHelix[0]+TMath::Sqrt(1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4])/fHelix[4];
+  Double_t dummy = 1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4];
+  if (dummy<0) {
+    AliError(Form("The argument of the Sqrt is %f => set to 0\n",dummy));
+    dummy = 0;
+  }
+  yc  =  fHelix[0]+TMath::Sqrt(dummy)/fHelix[4];
   
   fHelix[6] = xc*cs - yc*sn;
   fHelix[7] = xc*sn + yc*cs;
@@ -88,17 +95,24 @@ AliHelix::AliHelix(const AliExternalTrackParam &t)
   Double_t alpha,x,cs,sn;
   const Double_t *param =t.GetParameter(); 
   for (Int_t i=0;i<5;i++) fHelix[i]=param[i]; 
-  x = t.X();
-  alpha=t.Alpha();
+  x = t.GetX();
+  alpha=t.GetAlpha();
   //
   //circle parameters
-  fHelix[4]=fHelix[4]/AliKalmanTrack::GetConvConst();    // C
+  //PH Sometimes fP4 and fHelix[4] are very big and the calculation
+  //PH of the Sqrt cannot be done. To be investigated...
+  fHelix[4]=fHelix[4]/(-1000/0.299792458/AliTracker::GetBz());    // C
   cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
 
   Double_t xc, yc, rc;
   rc  =  1/fHelix[4];
   xc  =  x-fHelix[2]*rc;
-  yc  =  fHelix[0]+TMath::Sqrt(1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4])/fHelix[4];
+  Double_t dummy = 1-(x-xc)*(x-xc)*fHelix[4]*fHelix[4];
+  if (dummy<0) {
+    AliError(Form("The argument of the Sqrt is %f => set to 0\n",dummy));
+    dummy = 0;
+  }
+  yc  =  fHelix[0]+TMath::Sqrt(dummy)/fHelix[4];
   
   fHelix[6] = xc*cs - yc*sn;
   fHelix[7] = xc*sn + yc*cs;
@@ -123,7 +137,7 @@ AliHelix::AliHelix(Double_t x[3], Double_t p[3], Double_t charge, Double_t conve
   //
   Double_t pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
   if (TMath::Abs(conversion)<0.00000001) 
-    conversion = AliKalmanTrack::GetConvConst();
+    conversion = -1000/0.299792458/AliTracker::GetBz();
   //
   //  
   fHelix[4] = charge/(conversion*pt); // C
@@ -161,7 +175,7 @@ void  AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion, D
   // return  momentum at given phase
   Double_t x[3],g[3],gg[3];
   Evaluate(phase,x,g,gg);
-  if (TMath::Abs(conversion)<0.0001) conversion = AliKalmanTrack::GetConvConst();
+  if (TMath::Abs(conversion)<0.0001) conversion = -1000/0.299792458/AliTracker::GetBz();
   Double_t mt = TMath::Sqrt(g[0]*g[0]+g[1]*g[1]);
   p[0] = fHelix[8]*g[0]/(mt*conversion);
   p[1] = fHelix[8]*g[1]/(mt*conversion);
@@ -341,6 +355,11 @@ Int_t    AliHelix::GetRPHIintersections(AliHelix &h, Double_t phase[2][2], Doubl
   //  Double_t * c1 = &fHelix[6];
   //Double_t * c2 = &(h.fHelix[6]);
   //  Double_t  c1[3] = {fHelix[5],fHelix[0],fHelix[8]};
+
+  // PH initiaziation in case of return
+  phase[0][0]=phase[0][1]=phase[1][0]=phase[1][1]=0;
+  ri[0]=ri[1]=1000000;
+
   Double_t  c1[3] = {0,0,fHelix[8]};
   Double_t  c2[3] = {h.fHelix[5]-fHelix[5],h.fHelix[0]-fHelix[0],h.fHelix[8]};
 
@@ -355,9 +374,9 @@ Int_t    AliHelix::GetRPHIintersections(AliHelix &h, Double_t phase[2][2], Doubl
     x0[0] = (d+c1[2]-c2[2])*c2[0]/(2*d)+ fHelix[5];
     y0[0] = (d+c1[2]-c2[2])*c2[1]/(2*d)+ fHelix[0];
     //    return 0;
-    phase[0][0] = GetPhase(x0[0],y0[0]);
-    phase[0][1] = h.GetPhase(x0[0],y0[0]);
-    ri[0] = x0[0]*x0[0]+y0[0]*y0[0];
+    phase[1][0] = phase[0][0] = GetPhase(x0[0],y0[0]);
+    phase[1][1] = phase[0][1] = h.GetPhase(x0[0],y0[0]);
+    ri[1] = ri[0] = x0[0]*x0[0]+y0[0]*y0[0];
     return 1;
   }
   if ( (d+c2[2])<c1[2]){
@@ -365,12 +384,12 @@ Int_t    AliHelix::GetRPHIintersections(AliHelix &h, Double_t phase[2][2], Doubl
     //
     Double_t xx = c2[0]+ c2[0]*c2[2]/d+ fHelix[5];
     Double_t yy = c2[1]+ c2[1]*c2[2]/d+ fHelix[0]; 
-    phase[0][1] = h.GetPhase(xx,yy);
+    phase[1][1] = phase[0][1] = h.GetPhase(xx,yy);
     //
     Double_t xx2 = c2[0]*c1[2]/d+ fHelix[5];
     Double_t yy2 = c2[1]*c1[2]/d+ fHelix[0]; 
-    phase[0][0] = GetPhase(xx2,yy2);
-    ri[0] = xx*xx+yy*yy;
+    phase[1][0] = phase[0][0] = GetPhase(xx2,yy2);
+    ri[1] = ri[0] = xx*xx+yy*yy;
     return 1;
   }
 
@@ -379,12 +398,12 @@ Int_t    AliHelix::GetRPHIintersections(AliHelix &h, Double_t phase[2][2], Doubl
     //
     Double_t xx = -c2[0]*c1[2]/d+ fHelix[5];
     Double_t yy = -c2[1]*c1[2]/d+ fHelix[0]; 
-    phase[0][1] = GetPhase(xx,yy);
+    phase[1][1] = phase[0][1] = GetPhase(xx,yy);
     //
     Double_t xx2 = c2[0]- c2[0]*c2[2]/d+ fHelix[5];
     Double_t yy2 = c2[1]- c2[1]*c2[2]/d+ fHelix[0]; 
-    phase[0][0] = h.GetPhase(xx2,yy2);
-    ri[0] = xx*xx+yy*yy;
+    phase[1][0] = phase[0][0] = h.GetPhase(xx2,yy2);
+    ri[1] = ri[0] = xx*xx+yy*yy;
     return 1;
   }