Improve precision of AliExternalTrackParam:: GetZat, GetXYZat

Suppress unnecessary checks in AliAODTrack

diff --git a/STEER/AOD/AliAODTrack.cxx b/STEER/AOD/AliAODTrack.cxx
index a33d5baf584da800732e65c09f067e25a93ff86e..269d908e6ef9d65a5521cdae13949c8590bec3c5 100644 (file)
--- a/STEER/AOD/AliAODTrack.cxx
+++ b/STEER/AOD/AliAODTrack.cxx
@@ -900,7 +900,6 @@ Bool_t AliAODTrack::GetXYZAt(Double_t x, Double_t b, Double_t *r) const
   //conversion of track parameter representation is
   //based on the implementation of AliExternalTrackParam::Set(...)
   //maybe some of this code can be moved to AliVTrack to avoid code duplication
-  const double kSafe = 1e-5;
   Double_t alpha=0.0;
   Double_t radPos2 = fPosition[0]*fPosition[0]+fPosition[1]*fPosition[1];  
   Double_t radMax  = 45.; // approximately ITS outer radius
@@ -912,27 +911,10 @@ Bool_t AliAODTrack::GetXYZAt(Double_t x, Double_t b, Double_t *r) const
      TMath::DegToRad()*(20*((((Int_t)(phiPos*TMath::RadToDeg()))/20))+10);
   }
   //
-  Double_t cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
-  // protection:  avoid alpha being too close to 0 or +-pi/2
-  if (TMath::Abs(sn)<kSafe) {
-    alpha = kSafe;
-    cs=TMath::Cos(alpha);
-    sn=TMath::Sin(alpha);
-  }
-  else if (cs<kSafe) {
-    alpha -= TMath::Sign(kSafe, alpha);
-    cs=TMath::Cos(alpha);
-    sn=TMath::Sin(alpha);    
-  }
-  
   // Get the vertex of origin and the momentum
   TVector3 ver(fPosition[0],fPosition[1],fPosition[2]);
   TVector3 mom(Px(),Py(),Pz());
   //
-  // avoid momenta along axis
-  if (TMath::Abs(mom[0])<kSafe) mom[0] = TMath::Sign(kSafe*TMath::Abs(mom[1]), mom[0]);
-  if (TMath::Abs(mom[1])<kSafe) mom[1] = TMath::Sign(kSafe*TMath::Abs(mom[0]), mom[1]);
-
   // Rotate to the local coordinate system
   ver.RotateZ(-alpha);
   mom.RotateZ(-alpha);

diff --git a/STEER/STEERBase/AliExternalTrackParam.cxx b/STEER/STEERBase/AliExternalTrackParam.cxx
index 18416deb551f247a7ebbd2e74efbdfe4a7398c78..28fd45bdd3f97e01013b7a2615eca26ed5fe9037 100644 (file)
--- a/STEER/STEERBase/AliExternalTrackParam.cxx
+++ b/STEER/STEERBase/AliExternalTrackParam.cxx
@@ -1810,14 +1810,28 @@ AliExternalTrackParam::GetZAt(Double_t x, Double_t b, Double_t &z) const {
   //---------------------------------------------------------------------
   Double_t dx=x-fX;
   if(TMath::Abs(dx)<=kAlmost0) {z=fP[1]; return kTRUE;}
-
-  Double_t f1=fP[2], f2=f1 + dx*GetC(b);
+  Double_t crv=GetC(b);
+  Double_t x2r = crv*dx;
+  Double_t f1=fP[2], f2=f1 + x2r;
 
   if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
   if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
   
   Double_t r1=sqrt((1.-f1)*(1.+f1)), r2=sqrt((1.-f2)*(1.+f2));
-  z = fP[1] + dx*(r2 + f2*(f1+f2)/(r1+r2))*fP[3]; // Many thanks to P.Hristov !
+  double dy2dx = (f1+f2)/(r1+r2);
+  if (TMath::Abs(x2r)<0.05) {
+    z = fP[1] + dx*(r2 + f2*dy2dx)*fP[3]; // Many thanks to P.Hristov !    
+  }
+  else {
+    // for small dx/R the linear apporximation of the arc by the segment is OK,
+    // but at large dx/R the error is very large and leads to incorrect Z propagation
+    // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
+    // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
+    // Similarly, the rotation angle in linear in dx only for dx<<R
+    double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx);   // distance from old position to new one
+    double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
+    z = fP[1] + rot/crv*fP[3];    
+  }
   return kTRUE;
 }
 
@@ -1829,17 +1843,29 @@ AliExternalTrackParam::GetXYZAt(Double_t x, Double_t b, Double_t *r) const {
   //---------------------------------------------------------------------
   Double_t dx=x-fX;
   if(TMath::Abs(dx)<=kAlmost0) return GetXYZ(r);
-
-  Double_t f1=fP[2], f2=f1 + dx*GetC(b);
-
+  Double_t crv=GetC(b);
+  Double_t x2r = crv*dx;
+  Double_t f1=fP[2], f2=f1 + dx*crv;
   if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
   if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
   
   Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
+  double dy2dx = (f1+f2)/(r1+r2);
   r[0] = x;
-  r[1] = fP[0] + dx*(f1+f2)/(r1+r2);
-  r[2] = fP[1] + dx*(r2 + f2*(f1+f2)/(r1+r2))*fP[3];//Thanks to Andrea & Peter
-
+  r[1] = fP[0] + dx*dy2dx;
+  if (TMath::Abs(x2r)<0.05) {
+    r[2] = fP[1] + dx*(r2 + f2*dy2dx)*fP[3];//Thanks to Andrea & Peter
+  }
+  else {
+    // for small dx/R the linear apporximation of the arc by the segment is OK,
+    // but at large dx/R the error is very large and leads to incorrect Z propagation
+    // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
+    // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
+    // Similarly, the rotation angle in linear in dx only for dx<<R
+    double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx);   // distance from old position to new one
+    double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
+    r[2] = fP[1] + rot/crv*fP[3];
+  }
   return Local2GlobalPosition(r,fAlpha);
 }