]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliExternalTrackParam.cxx
Updates.
[u/mrichter/AliRoot.git] / STEER / AliExternalTrackParam.cxx
index 2f4eeda5cc61390424b90b689f5fb946d686ff7e..79a27dfaabbe43270f93912dea5a1689bc6c3ca1 100644 (file)
 // are implemented.
 // Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch                            //
 ///////////////////////////////////////////////////////////////////////////////
+#include <TMatrixDSym.h>
 #include "AliExternalTrackParam.h"
 #include "AliESDVertex.h"
 #include "AliLog.h"
 
 ClassImp(AliExternalTrackParam)
 
+Double32_t AliExternalTrackParam::fgMostProbablePt=kMostProbablePt;
 //_____________________________________________________________________________
 AliExternalTrackParam::AliExternalTrackParam() :
   TObject(),
@@ -178,9 +181,63 @@ Double_t AliExternalTrackParam::GetLinearD(Double_t xv,Double_t yv) const {
   return -d;
 }
 
+Bool_t AliExternalTrackParam::CorrectForMeanMaterial
+(Double_t xOverX0,  Double_t xTimesRho, Double_t mass, 
+Double_t (*Bethe)(Double_t)) {
+  //------------------------------------------------------------------
+  // This function corrects the track parameters for the crossed material.
+  // "xOverX0"   - X/X0, the thickness in units of the radiation length.
+  // "xTimesRho" - is the product length*density (g/cm^2). 
+  // "mass" - the mass of this particle (GeV/c^2).
+  //------------------------------------------------------------------
+  Double_t &fP2=fP[2];
+  Double_t &fP3=fP[3];
+  Double_t &fP4=fP[4];
+
+  Double_t &fC22=fC[5];
+  Double_t &fC33=fC[9];
+  Double_t &fC43=fC[13];
+  Double_t &fC44=fC[14];
+
+  Double_t p=GetP();
+  Double_t p2=p*p;
+  Double_t beta2=p2/(p2 + mass*mass);
+  xOverX0*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
+
+  //Multiple scattering******************
+  if (xOverX0 != 0) {
+     Double_t theta2=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(xOverX0);
+     //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*TMath::Abs(d)*9.36*2.33;
+     fC22 += theta2*(1.- fP2*fP2)*(1. + fP3*fP3);
+     fC33 += theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
+     fC43 += theta2*fP3*fP4*(1. + fP3*fP3);
+     fC44 += theta2*fP3*fP4*fP3*fP4;
+  }
+
+  //Energy losses************************
+  if ((xTimesRho != 0.) && (beta2 < 1.)) {
+     Double_t dE=Bethe(beta2)*xTimesRho;
+     Double_t e=TMath::Sqrt(p2 + mass*mass);
+     if ( TMath::Abs(dE) > 0.3*e ) return kFALSE; //30% energy loss is too much!
+     fP4*=(1.- e/p2*dE);
+
+     // Approximate energy loss fluctuation (M.Ivanov)
+     const Double_t knst=0.07; // To be tuned.  
+     Double_t sigmadE=knst*TMath::Sqrt(TMath::Abs(dE)); 
+     fC44+=((sigmadE*e/p2*fP4)*(sigmadE*e/p2*fP4)); 
+  }
+
+  return kTRUE;
+}
+
+
 Bool_t AliExternalTrackParam::CorrectForMaterial
 (Double_t d,  Double_t x0, Double_t mass, Double_t (*Bethe)(Double_t)) {
   //------------------------------------------------------------------
+  //                    Deprecated function !   
+  //       Better use CorrectForMeanMaterial instead of it.
+  //
   // This function corrects the track parameters for the crossed material
   // "d"    - the thickness (fraction of the radiation length)
   // "x0"   - the radiation length (g/cm^2) 
@@ -215,11 +272,12 @@ Bool_t AliExternalTrackParam::CorrectForMaterial
      d*=x0;
      Double_t dE=Bethe(beta2)*d;
      Double_t e=TMath::Sqrt(p2 + mass*mass);
+     if ( TMath::Abs(dE) > 0.3*e ) return kFALSE; //30% energy loss is too much!
      fP4*=(1.- e/p2*dE);
 
      // Approximate energy loss fluctuation (M.Ivanov)
-     const Double_t cnst=0.07; // To be tuned.  
-     Double_t sigmadE=cnst*TMath::Sqrt(TMath::Abs(dE)); 
+     const Double_t knst=0.07; // To be tuned.  
+     Double_t sigmadE=knst*TMath::Sqrt(TMath::Abs(dE)); 
      fC44+=((sigmadE*e/p2*fP4)*(sigmadE*e/p2*fP4)); 
  
   }
@@ -380,7 +438,7 @@ Double_t p[3], Double_t bz) const {
   //+++++++++++++++++++++++++++++++++++++++++    
   GetXYZ(x);
     
-  if (TMath::Abs(Get1Pt()) < kAlmost0){ //straight-line tracks
+  if (TMath::Abs(Get1Pt()) < kAlmost0 || TMath::Abs(bz) < kAlmost0Field ){ //straight-line tracks
      Double_t unit[3]; GetDirection(unit);
      x[0]+=unit[0]*len;   
      x[1]+=unit[1]*len;   
@@ -456,6 +514,114 @@ AliExternalTrackParam::GetPredictedChi2(Double_t p[2],Double_t cov[3]) const {
   return (d*szz*d - 2*d*sdz*z + z*sdd*z)/det;
 }
 
+Double_t AliExternalTrackParam::
+GetPredictedChi2(Double_t p[3],Double_t covyz[3],Double_t covxyz[3]) const {
+  //----------------------------------------------------------------
+  // Estimate the chi2 of the 3D space point "p" and
+  // the full covariance matrix "covyz" and "covxyz"
+  //
+  // Cov(x,x) ... :   covxyz[0]
+  // Cov(y,x) ... :   covxyz[1]  covyz[0]
+  // Cov(z,x) ... :   covxyz[2]  covyz[1]  covyz[2]
+  //----------------------------------------------------------------
+
+  Double_t res[3] = {
+    GetX() - p[0],
+    GetY() - p[1],
+    GetZ() - p[2]
+  };
+
+  Double_t f=GetSnp();
+  if (TMath::Abs(f) >= kAlmost1) return kVeryBig;
+  Double_t r=TMath::Sqrt(1.- f*f);
+  Double_t a=f/r, b=GetTgl()/r;
+
+  Double_t s2=333.*333.;  //something reasonably big (cm^2)
+  TMatrixDSym v(3);
+  v(0,0)=  s2;  v(0,1)=  a*s2;                 v(0,2)=  b*s2;;
+  v(1,0)=a*s2;  v(1,1)=a*a*s2 + GetSigmaY2();  v(1,2)=a*b*s2 + GetSigmaZY();
+  v(2,0)=b*s2;  v(2,1)=a*b*s2 + GetSigmaZY();  v(2,2)=b*b*s2 + GetSigmaZ2();
+
+  v(0,0)+=covxyz[0]; v(0,1)+=covxyz[1]; v(0,2)+=covxyz[2];
+  v(1,0)+=covxyz[1]; v(1,1)+=covyz[0];  v(1,2)+=covyz[1];
+  v(2,0)+=covxyz[2]; v(2,1)+=covyz[1];  v(2,2)+=covyz[2];
+
+  v.Invert();
+  if (!v.IsValid()) return kVeryBig;
+
+  Double_t chi2=0.;
+  for (Int_t i = 0; i < 3; i++)
+    for (Int_t j = 0; j < 3; j++) chi2 += res[i]*res[j]*v(i,j);
+
+  return chi2;  
+
+
+}
+
+Bool_t AliExternalTrackParam::
+PropagateTo(Double_t p[3],Double_t covyz[3],Double_t covxyz[3],Double_t bz) {
+  //----------------------------------------------------------------
+  // Propagate this track to the plane 
+  // the 3D space point "p" (with the covariance matrix "covyz" and "covxyz")
+  // belongs to.
+  // The magnetic field is "bz" (kG)
+  //
+  // The track curvature and the change of the covariance matrix
+  // of the track parameters are negleted !
+  // (So the "step" should be small compared with 1/curvature)
+  //----------------------------------------------------------------
+
+  Double_t f=GetSnp();
+  if (TMath::Abs(f) >= kAlmost1) return kFALSE;
+  Double_t r=TMath::Sqrt(1.- f*f);
+  Double_t a=f/r, b=GetTgl()/r;
+
+  Double_t s2=333.*333.;  //something reasonably big (cm^2)
+  TMatrixDSym tV(3);
+  tV(0,0)=  s2;  tV(0,1)=  a*s2;  tV(0,2)=  b*s2;
+  tV(1,0)=a*s2;  tV(1,1)=a*a*s2;  tV(1,2)=a*b*s2;
+  tV(2,0)=b*s2;  tV(2,1)=a*b*s2;  tV(2,2)=b*b*s2;
+
+  TMatrixDSym pV(3);
+  pV(0,0)=covxyz[0]; pV(0,1)=covxyz[1]; pV(0,2)=covxyz[2];
+  pV(1,0)=covxyz[1]; pV(1,1)=covyz[0];  pV(1,2)=covyz[1];
+  pV(2,0)=covxyz[2]; pV(2,1)=covyz[1];  pV(2,2)=covyz[2];
+
+  TMatrixDSym tpV(tV);
+  tpV+=pV;
+  tpV.Invert();
+  if (!tpV.IsValid()) return kFALSE;
+
+  TMatrixDSym pW(3),tW(3);
+  for (Int_t i=0; i<3; i++)
+    for (Int_t j=0; j<3; j++) {
+      pW(i,j)=tW(i,j)=0.;
+      for (Int_t k=0; k<3; k++) {
+       pW(i,j) += tV(i,k)*tpV(k,j);
+       tW(i,j) += pV(i,k)*tpV(k,j);
+      }
+    }
+
+  Double_t t[3] = {GetX(), GetY(), GetZ()};
+
+  Double_t x=0.;
+  for (Int_t i=0; i<3; i++) x += (tW(0,i)*t[i] + pW(0,i)*p[i]);  
+  Double_t crv=GetC(bz);
+  if (TMath::Abs(b) < kAlmost0Field) crv=0.;
+  f += crv*(x-fX);
+  if (TMath::Abs(f) >= kAlmost1) return kFALSE;
+  fX=x;  
+
+  fP[0]=0.;
+  for (Int_t i=0; i<3; i++) fP[0] += (tW(1,i)*t[i] + pW(1,i)*p[i]);  
+  fP[1]=0.;
+  for (Int_t i=0; i<3; i++) fP[1] += (tW(2,i)*t[i] + pW(2,i)*p[i]);  
+
+  return kTRUE;  
+}
+
 Bool_t AliExternalTrackParam::Update(Double_t p[2], Double_t cov[3]) {
   //------------------------------------------------------------------
   // Update the track parameters with the space point "p" having
@@ -824,6 +990,10 @@ Bool_t AliExternalTrackParam::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
   Double_t m24= pt*(cs - fP[2]*sn/r), m44=-pt*pt*(r*sn + fP[2]*cs);
   Double_t m35=pt, m45=-pt*pt*fP[3];
 
+  m43*=GetSign();
+  m44*=GetSign();
+  m45*=GetSign();
+
   cv[0 ] = fC[0]*m00*m00;
   cv[1 ] = fC[0]*m00*m10; 
   cv[2 ] = fC[0]*m10*m10;