Added method to get local X at arbitrary R, moving in selected direction

diff --git a/STEER/STEERBase/AliExternalTrackParam.cxx b/STEER/STEERBase/AliExternalTrackParam.cxx
index 0899a6fb0c76773ea2997692cd458c254471b7f1..f300ef3aafe865264215f2a79d17679c1b497327 100644 (file)
--- a/STEER/STEERBase/AliExternalTrackParam.cxx
+++ b/STEER/STEERBase/AliExternalTrackParam.cxx
@@ -2300,3 +2300,121 @@ Bool_t AliExternalTrackParam::ConstrainToVertex(const AliVVertex* vtx, Double_t
 
   return kTRUE;
 }
+
+//___________________________________________________________________________________________
+Bool_t AliExternalTrackParam::GetXatLabR(Double_t r,Double_t &x, Double_t bz, Int_t dir) const
+{
+  // Get local X of the track position estimated at the radius lab radius r. 
+  // The track curvature is accounted exactly
+  //
+  // The flag "dir" can be used to remove the ambiguity of which intersection to take (out of 2 possible)
+  // 0  - take the intersection closest to the current track position
+  // >0 - go along the track (increasing fX)
+  // <0 - go backward (decreasing fX)
+  //
+  const Double_t &fy=fP[0], &sn = fP[2];
+  //
+  double crv = GetC(bz);
+  if (TMath::Abs(crv)<=kAlmost0) { // this is a straight track
+    if (TMath::Abs(sn)>=kAlmost1) { // || to Y axis
+      double det = (r-fX)*(r+fX);
+      if (det<0) return kFALSE;     // does not reach raduis r
+      x = fX;
+      if (dir==0) return kTRUE;
+      det = TMath::Sqrt(det);
+      if (dir>0) {                       // along the track direction
+       if (sn>0) {if (fy>det)  return kFALSE;} // track is along Y axis and above the circle
+       else      {if (fy<-det) return kFALSE;} // track is against Y axis amd belo the circle
+      }
+      else if(dir>0) {                                    // agains track direction
+       if (sn>0) {if (fy<-det) return kFALSE;} // track is along Y axis
+        else if (fy>det)  return kFALSE;        // track is against Y axis
+      }
+    }
+    else if (TMath::Abs(sn)<=kAlmost0) { // || to X axis
+      double det = (r-fy)*(r+fy);
+      if (det<0) return kFALSE;     // does not reach raduis r
+      det = TMath::Sqrt(det);
+      if (!dir) {
+       x = fX>0  ? det : -det;    // choose the solution requiring the smalest step
+       return kTRUE;
+      }
+      else if (dir>0) {                    // along the track direction
+       if      (fX > det) return kFALSE;  // current point is in on the right from the circle
+       else if (fX <-det) x = -det;       // on the left
+       else               x =  det;       // within the circle
+      }
+      else {                               // against the track direction
+       if      (fX <-det) return kFALSE;  
+       else if (fX > det) x =  det;
+       else               x = -det;
+      }
+    }
+    else {                                 // general case of straight line
+      double cs = TMath::Sqrt((1-sn)*(1+sn));
+      double xsyc = fX*sn-fy*cs;
+      double det = (r-xsyc)*(r+xsyc);
+      if (det<0) return kFALSE;    // does not reach raduis r
+      det = TMath::Sqrt(det);
+      double xcys = fX*cs+fy*sn;
+      double t = -xcys;
+      if (dir==0) t += t>0 ? -det:det;  // chose the solution requiring the smalest step
+      else if (dir>0) {                 // go in increasing fX direction. ( t+-det > 0)
+       if (t>=-det) t += -det;         // take minimal step giving t>0
+       else return kFALSE;             // both solutions have negative t
+      }
+      else {                            // go in increasing fX direction. (t+-det < 0)
+       if (t<det) t -= det;            // take minimal step giving t<0
+       else return kFALSE;             // both solutions have positive t
+      }
+      x = fX + cs*t;
+    }
+  }
+  else {                                 // helix
+    // get center of the track circle
+    double tR = 1./crv;   // track radius (for the moment signed)
+    double cs = TMath::Sqrt((1-sn)*(1+sn));
+    double x0 = fX - sn*tR;
+    double y0 = fy + cs*tR;
+    double r0 = TMath::Sqrt(x0*x0+y0*y0);
+    //    printf("Xc:%+e Yc:%+e\n",x0,y0);
+    //
+    if (r0<=kAlmost0) return kFALSE;            // the track is concentric to circle
+    tR = TMath::Abs(tR);
+    double tR2r0 = tR/r0;
+    double g = 0.5*(r*r/(r0*tR) - tR2r0 - 1./tR2r0);
+    double det = (1.-g)*(1.+g);
+    if (det<0) return kFALSE;         // does not reach raduis r
+    det = TMath::Sqrt(det);
+    //
+    // the intersection happens in 2 points: {x0+tR*C,y0+tR*S} 
+    // with C=f*c0+-|s0|*det and S=f*s0-+c0 sign(s0)*det
+    // where s0 and c0 make direction for the circle center (=x0/r0 and y0/r0)
+    //
+    double tmp = 1.+g*tR2r0;
+    x = x0*tmp; 
+    double y = y0*tmp;
+    if (TMath::Abs(y0)>kAlmost0) { // when y0==0 the x,y is unique
+      double dfx = tR2r0*TMath::Abs(y0)*det;
+      double dfy = tR2r0*x0*TMath::Sign(det,y0);
+      if (dir==0) {                    // chose the one which corresponds to smallest step 
+       double delta = (x-fX)*dfx-(y-fy)*dfy; // the choice of + in C will lead to smaller step if delta<0
+       if (delta<0) x += dfx;
+       else         x -= dfx;
+      }
+      else if (dir>0) {  // along track direction: x must be > fX
+       x -= dfx; // try the smallest step (dfx is positive)
+       if (x<fX && (x+=dfx+dfx)<fX) return kFALSE;
+      }
+      else { // backward: x must be < fX
+       x += dfx; // try the smallest step (dfx is positive)
+       if (x>fX && (x-=dfx+dfx)>fX) return kFALSE;
+      }
+    }
+    else { // special case: track touching the circle just in 1 point
+      if ( (dir>0&&x<fX) || (dir<0&&x>fX) ) return kFALSE; 
+    }
+  }
+  //
+  return kTRUE;
+}

diff --git a/STEER/STEERBase/AliExternalTrackParam.h b/STEER/STEERBase/AliExternalTrackParam.h
index 8b8f5f171b029a8444f01b83c204cc0c0925eca3..fb0e98b49c00aa38d7fa4d41d55e97224f22d774 100644 (file)
--- a/STEER/STEERBase/AliExternalTrackParam.h
+++ b/STEER/STEERBase/AliExternalTrackParam.h
@@ -216,6 +216,7 @@ class AliExternalTrackParam: public AliVTrack {
   Bool_t GetZAt(Double_t x,  Double_t b,  Double_t &z) const;
   void Print(Option_t* option = "") const;
   Double_t GetSnpAt(Double_t x,Double_t b) const;
+  Bool_t GetXatLabR(Double_t r,Double_t &x, Double_t bz, Int_t dir=0) const;
 
   //Deprecated
   Bool_t CorrectForMaterial(Double_t d, Double_t x0, Double_t mass,