]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliExternalTrackParam.cxx
Additional protection
[u/mrichter/AliRoot.git] / STEER / AliExternalTrackParam.cxx
index 8b94b8211cb074ffc28ea023d6940d7dd7fd2aab..d91056f91f003d703c20b81f42f77b7a523a4d2e 100644 (file)
@@ -27,7 +27,7 @@
 ///////////////////////////////////////////////////////////////////////////////
 #include "AliExternalTrackParam.h"
 #include "AliKalmanTrack.h"
-#include "AliTracker.h"
+#include "AliESDVertex.h"
 
 
 ClassImp(AliExternalTrackParam)
@@ -89,10 +89,17 @@ Double_t AliExternalTrackParam::GetP() const {
   // This function returns the track momentum
   // Results for (nearly) straight tracks are meaningless !
   //---------------------------------------------------------------------
-  if (TMath::Abs(fP[4])<=0) return 0;
+  if (TMath::Abs(fP[4])<=kAlmost0) return kVeryBig;
   return TMath::Sqrt(1.+ fP[3]*fP[3])/TMath::Abs(fP[4]);
 }
 
+Double_t AliExternalTrackParam::Get1P() const {
+  //---------------------------------------------------------------------
+  // This function returns the 1/(track momentum)
+  //---------------------------------------------------------------------
+  return TMath::Abs(fP[4])/TMath::Sqrt(1.+ fP[3]*fP[3]);
+}
+
 //_______________________________________________________________________
 Double_t AliExternalTrackParam::GetD(Double_t x,Double_t y,Double_t b) const {
   //------------------------------------------------------------------
@@ -100,6 +107,7 @@ Double_t AliExternalTrackParam::GetD(Double_t x,Double_t y,Double_t b) const {
   // with respect to a point with global coordinates (x,y)
   // in the magnetic field "b" (kG)
   //------------------------------------------------------------------
+  if (TMath::Abs(b) < kAlmost0Field) return GetLinearD(x,y);
   Double_t rp4=kB2C*b*fP[4];
 
   Double_t xt=fX, yt=fP[0];
@@ -163,7 +171,7 @@ CorrectForMaterial(Double_t d,  Double_t x0, Double_t mass) {
   }
 
   //Energy losses************************
-  if (x0!=0.) {
+  if (x0!=0. && beta2<1) {
      d*=x0;
      Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d;
      if (beta2/(1-beta2)>3.5*3.5)
@@ -204,6 +212,11 @@ Bool_t AliExternalTrackParam::Rotate(Double_t alpha) {
   fP0= -x*sa + fP0*ca;
   fP2=  sf*ca - cf*sa;
 
+  if (TMath::Abs(cf)<kAlmost0) {
+    AliError(Form("Too small cosine value %f",cf)); 
+    cf = kAlmost0;
+  } 
+
   Double_t rr=(ca+sf/cf*sa);  
 
   fC00 *= (ca*ca);
@@ -223,8 +236,12 @@ Bool_t AliExternalTrackParam::PropagateTo(Double_t xk, Double_t b) {
   //----------------------------------------------------------------
   // Propagate this track to the plane X=xk (cm) in the field "b" (kG)
   //----------------------------------------------------------------
-  Double_t crv=kB2C*b*fP[4];
   Double_t dx=xk-fX;
+  if (TMath::Abs(dx)<=kAlmost0)  return kTRUE;
+
+  Double_t crv=kB2C*b*fP[4];
+  if (TMath::Abs(b) < kAlmost0Field) crv=0.;
+
   Double_t f1=fP[2], f2=f1 + crv*dx;
   if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
   if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
@@ -241,7 +258,7 @@ Bool_t AliExternalTrackParam::PropagateTo(Double_t xk, Double_t b) {
 
   fX=xk;
   fP0 += dx*(f1+f2)/(r1+r2);
-  fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
+  fP1 += dx*(r2 + f2*(f1+f2)/(r1+r2))*fP3;  // Many thanks to P.Hristov !
   fP2 += dx*crv;
 
   //f = F - 1
@@ -527,6 +544,45 @@ PropagateToDCA(AliExternalTrackParam *p, Double_t b) {
 
 
 
+
+Bool_t AliExternalTrackParam::PropagateToDCA(const AliESDVertex *vtx, Double_t b, Double_t maxd){
+  //
+  // Try to relate this track to the vertex "vtx", 
+  // if the (rough) transverse impact parameter is not bigger then "maxd". 
+  //            Magnetic field is "b" (kG).
+  //
+  // a) The track gets extapolated to the DCA to the vertex.
+  // b) The impact parameters and their covariance matrix are calculated.
+  //
+  //    In the case of success, the returned value is kTRUE
+  //    (otherwise, it's kFALSE)
+  //  
+  Double_t alpha=GetAlpha();
+  Double_t sn=TMath::Sin(alpha), cs=TMath::Cos(alpha);
+  Double_t x=GetX(), y=GetParameter()[0], snp=GetParameter()[2];
+  Double_t xv= vtx->GetXv()*cs + vtx->GetYv()*sn;
+  Double_t yv=-vtx->GetXv()*sn + vtx->GetYv()*cs;
+  x-=xv; y-=yv;
+
+  //Estimate the impact parameter neglecting the track curvature
+  Double_t d=TMath::Abs(x*snp - y*TMath::Sqrt(1.- snp*snp));
+  if (d > maxd) return kFALSE; 
+
+  //Propagate to the DCA
+  Double_t crv=0.299792458e-3*b*GetParameter()[4];
+  Double_t tgfv=-(crv*x - snp)/(crv*y + TMath::Sqrt(1.-snp*snp));
+  sn=tgfv/TMath::Sqrt(1.+ tgfv*tgfv); cs=TMath::Sqrt(1.- sn*sn);
+
+  x = xv*cs + yv*sn;
+  yv=-xv*sn + yv*cs; xv=x;
+
+  if (!Propagate(alpha+TMath::ASin(sn),xv,b)) return kFALSE;
+  return kTRUE;
+}
+
+
+
+
 Bool_t Local2GlobalMomentum(Double_t p[3],Double_t alpha) {
   //----------------------------------------------------------------
   // This function performs local->global transformation of the
@@ -542,7 +598,7 @@ Bool_t Local2GlobalMomentum(Double_t p[3],Double_t alpha) {
   //    p[2] = pz
   // Results for (nearly) straight tracks are meaningless !
   //----------------------------------------------------------------
-  if (TMath::Abs(p[0])<=0)        return kFALSE;
+  if (TMath::Abs(p[0])<=kAlmost0) return kFALSE;
   if (TMath::Abs(p[1])> kAlmost1) return kFALSE;
 
   Double_t pt=1./TMath::Abs(p[0]);
@@ -603,7 +659,7 @@ Bool_t AliExternalTrackParam::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
   //
   // Results for (nearly) straight tracks are meaningless !
   //---------------------------------------------------------------------
-  if (TMath::Abs(fP[4])<=0) {
+  if (TMath::Abs(fP[4])<=kAlmost0) {
      for (Int_t i=0; i<21; i++) cv[i]=0.;
      return kFALSE;
   }
@@ -665,8 +721,11 @@ AliExternalTrackParam::GetXYZAt(Double_t x, Double_t b, Double_t *r) const {
   // the radial position "x" (cm) in the magnetic field "b" (kG)
   //---------------------------------------------------------------------
   Double_t dx=x-fX;
+  if(TMath::Abs(dx)<=kAlmost0) return GetXYZ(r);
+
   Double_t f1=fP[2], f2=f1 + dx*fP[4]*b*kB2C;
 
+  if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
   if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
   
   Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
@@ -676,7 +735,6 @@ AliExternalTrackParam::GetXYZAt(Double_t x, Double_t b, Double_t *r) const {
   return Local2GlobalPosition(r,fAlpha);
 }
 
-
 //_____________________________________________________________________________
 void AliExternalTrackParam::Print(Option_t* /*option*/) const
 {
@@ -695,13 +753,16 @@ void AliExternalTrackParam::Print(Option_t* /*option*/) const
 }
 
 
-Bool_t AliExternalTrackParam::PropagateTo(Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t rotateTo){
+Bool_t AliExternalTrackParam::PropagateTo(Double_t xToGo, Double_t b, Double_t mass, Double_t maxStep, Bool_t rotateTo){
   //----------------------------------------------------------------
-  // Propagate this track to the plane X=xk (cm) 
-  // correction for unhomogenity of the magnetic field and the
+  //
+  // Very expensive function !  Don't abuse it !
+  //
+  // Propagates this track to the plane X=xk (cm) 
+  // in the magnetic field "b" (kG),
   // the correction for the material is included
   //
-  //  Require acces to magnetic field and geomanager
+  //  Requires acces to geomanager
   //
   // mass     - mass used in propagation - used for energy loss correction
   // maxStep  - maximal step for propagation
@@ -711,18 +772,21 @@ Bool_t AliExternalTrackParam::PropagateTo(Double_t xToGo, Double_t mass, Double_
   Double_t dir      = (xpos<xToGo) ? 1.:-1.;
   //
   while ( (xToGo-xpos)*dir > kEpsilon){
+    if (TMath::Abs(fP[2]) >= kAlmost1) return kFALSE;
     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
     Double_t x    = xpos+step;
     Double_t xyz0[3],xyz1[3],param[7];
     GetXYZ(xyz0);   //starting global position
-    Float_t  pos0[3] = {xyz0[0],xyz0[1],xyz0[2]};
-    Double_t magZ = AliTracker::GetBz(pos0);
-    if (!GetXYZAt(x,magZ,xyz1)) return kFALSE;   // no prolongation
+    if (!GetXYZAt(x,b,xyz1)) return kFALSE;   // no prolongation
     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);       
-    if (!PropagateTo(x,magZ))  return kFALSE;
-    Double_t distance = param[4];
-    if (!CorrectForMaterial(distance,param[1],param[0],mass)) return kFALSE;
+    if (!PropagateTo(x,b))  return kFALSE;
+
+    Double_t rho=param[0],x0=param[1],distance=param[4];
+    Double_t d=distance*rho/x0;
+
+    if (!CorrectForMaterial(d,x0,mass)) return kFALSE;
     if (rotateTo){
+      if (TMath::Abs(fP[2]) >= kAlmost1) return kFALSE;
       GetXYZ(xyz0);   // global position
       Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
       if (!Rotate(alphan)) return kFALSE;
@@ -732,37 +796,6 @@ Bool_t AliExternalTrackParam::PropagateTo(Double_t xToGo, Double_t mass, Double_
   return kTRUE;
 }
 
-//_____________________________________________________________________________
-Bool_t AliExternalTrackParam::CorrectForMaterial(Double_t d, Double_t x0, Double_t rho, Double_t mass)
-{
-  //
-  // Take into account material effects assuming:
-  // x0  - mean rad length
-  // rho - mean density
 
-  //
-  // multiple scattering
-  //
-  if (mass<=0) {
-    AliError("Non-positive mass");
-    return kFALSE;
-  }
-  Double_t p2=(1.+ fP[3]*fP[3])/(fP[4]*fP[4]);
-  Double_t beta2=p2/(p2 + mass*mass);
-  Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
-  //
-  fC[5] += theta2*(1.- fP[2]*fP[2])*(1. + fP[3]*fP[3]);
-  fC[9] += theta2*(1. + fP[3]*fP[3])*(1. + fP[3]*fP[3]);
-  fC[13] += theta2*fP[3]*fP[4]*(1. + fP[3]*fP[3]);
-  fC[14] += theta2*fP[3]*fP[4]*fP[3]*fP[4];
-  //
-  Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2+1e-10)) - beta2)*d*rho;  
-  fP[4] *=(1.- TMath::Sqrt(p2+mass*mass)/p2*dE);
-  //
-  Double_t sigmade = 0.02*TMath::Sqrt(TMath::Abs(dE));   // energy loss fluctuation 
-  Double_t sigmac2 = sigmade*sigmade*fP[4]*fP[4]*(p2+mass*mass)/(p2*p2);
-  fC[14] += sigmac2;
-  return kTRUE;
-}