Bug corrected.
[u/mrichter/AliRoot.git] / STEER / AliExternalTrackParam.cxx
index a6f0f67..295830b 100644 (file)
@@ -25,6 +25,9 @@
 // are implemented.
 // Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch                            //
 ///////////////////////////////////////////////////////////////////////////////
+#include <cassert>
+
+#include <TVectorD.h>
 #include <TMatrixDSym.h>
 #include <TPolyMarker3D.h>
 #include <TVector3.h>
@@ -62,6 +65,7 @@ AliExternalTrackParam::AliExternalTrackParam(const AliExternalTrackParam &track)
   //
   for (Int_t i = 0; i < 5; i++) fP[i] = track.fP[i];
   for (Int_t i = 0; i < 15; i++) fC[i] = track.fC[i];
+  CheckCovariance();
 }
 
 //_____________________________________________________________________________
@@ -78,6 +82,7 @@ AliExternalTrackParam& AliExternalTrackParam::operator=(const AliExternalTrackPa
 
     for (Int_t i = 0; i < 5; i++) fP[i] = trkPar.fP[i];
     for (Int_t i = 0; i < 15; i++) fC[i] = trkPar.fC[i];
+    CheckCovariance();
   }
 
   return *this;
@@ -96,6 +101,7 @@ AliExternalTrackParam::AliExternalTrackParam(Double_t x, Double_t alpha,
   //
   for (Int_t i = 0; i < 5; i++)  fP[i] = param[i];
   for (Int_t i = 0; i < 15; i++) fC[i] = covar[i];
+  CheckCovariance();
 }
 
 //_____________________________________________________________________________
@@ -234,6 +240,8 @@ void AliExternalTrackParam::Set(Double_t xyz[3],Double_t pxpypz[3],
   fC[13] = b1/b3-b2*fC[8]/b3;
   fC[9 ] = TMath::Abs((cv[20]-fC[14]*(m45*m45)-fC[13]*2.*m35*m45)/(m35*m35));
 
+  CheckCovariance();
+
   return;
 }
 
@@ -258,6 +266,7 @@ void AliExternalTrackParam::AddCovariance(const Double_t c[15]) {
     fC[3] +=c[3];  fC[4] +=c[4];  fC[5] +=c[5];
     fC[6] +=c[6];  fC[7] +=c[7];  fC[8] +=c[8];  fC[9] +=c[9];
     fC[10]+=c[10]; fC[11]+=c[11]; fC[12]+=c[12]; fC[13]+=c[13]; fC[14]+=c[14];
+    CheckCovariance();
 }
 
 
@@ -345,14 +354,17 @@ 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, Bool_t anglecorr, 
- Double_t (*Bethe)(Double_t)) {
+Bool_t AliExternalTrackParam::CorrectForMeanMaterialdEdx
+(Double_t xOverX0,  Double_t xTimesRho, Double_t mass, 
+ Double_t dEdx,
+ Bool_t anglecorr) {
   //------------------------------------------------------------------
   // 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).
+  // "dEdx" - mean enery loss (GeV/(g/cm^2)
+  // "anglecorr" - switch for the angular correction
   //------------------------------------------------------------------
   Double_t &fP2=fP[2];
   Double_t &fP3=fP[3];
@@ -392,10 +404,12 @@ Bool_t AliExternalTrackParam::CorrectForMeanMaterial
   //Calculating the energy loss corrections************************
   Double_t cP4=1.;
   if ((xTimesRho != 0.) && (beta2 < 1.)) {
-     Double_t dE=Bethe(p/mass)*xTimesRho;
+     Double_t dE=dEdx*xTimesRho;
      Double_t e=TMath::Sqrt(p2 + mass*mass);
      if ( TMath::Abs(dE) > 0.3*e ) return kFALSE; //30% energy loss is too much!
-     cP4 = (1.- e/p2*dE);
+     //cP4 = (1.- e/p2*dE);
+     if ( (1.+ dE/p2*(dE + 2*e)) < 0. ) return kFALSE;
+     cP4 = 1./TMath::Sqrt(1.+ dE/p2*(dE + 2*e));  //A precise formula by Ruben !
      if (TMath::Abs(fP4*cP4)>100.) return kFALSE; //Do not track below 10 MeV/c
 
 
@@ -413,9 +427,62 @@ Bool_t AliExternalTrackParam::CorrectForMeanMaterial
   fC44 += cC44;
   fP4  *= cP4;
 
+  CheckCovariance();
+
   return kTRUE;
 }
 
+Bool_t AliExternalTrackParam::CorrectForMeanMaterial
+(Double_t xOverX0,  Double_t xTimesRho, Double_t mass, 
+ Bool_t anglecorr,
+ 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).
+  // "anglecorr" - switch for the angular correction
+  // "Bethe" - function calculating the energy loss (GeV/(g/cm^2)) 
+  //------------------------------------------------------------------
+  
+  Double_t bg=GetP()/mass;
+  Double_t dEdx=Bethe(bg);
+
+  return CorrectForMeanMaterialdEdx(xOverX0,xTimesRho,mass,dEdx,anglecorr);
+}
+
+Bool_t AliExternalTrackParam::CorrectForMeanMaterialZA
+(Double_t xOverX0, Double_t xTimesRho, Double_t mass,
+ Double_t zOverA,
+ Double_t density,
+ Double_t exEnergy,
+ Double_t jp1,
+ Double_t jp2,
+ Bool_t anglecorr) {
+  //------------------------------------------------------------------
+  // This function corrects the track parameters for the crossed material
+  // using the full Geant-like Bethe-Bloch formula parameterization
+  // "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).
+  // "density"  - mean density (g/cm^3)
+  // "zOverA"   - mean Z/A
+  // "exEnergy" - mean exitation energy (GeV)
+  // "jp1"      - density effect first junction point
+  // "jp2"      - density effect second junction point
+  // "anglecorr" - switch for the angular correction
+  //
+  //  The default values of the parameters are for silicon 
+  //
+  //------------------------------------------------------------------
+
+  Double_t bg=GetP()/mass;
+  Double_t dEdx=BetheBlochGeant(bg,density,jp1,jp2,exEnergy,zOverA);
+
+  return CorrectForMeanMaterialdEdx(xOverX0,xTimesRho,mass,dEdx,anglecorr);
+}
+
+
 
 Bool_t AliExternalTrackParam::CorrectForMaterial
 (Double_t d,  Double_t x0, Double_t mass, Double_t (*Bethe)(Double_t)) {
@@ -428,58 +495,9 @@ Bool_t AliExternalTrackParam::CorrectForMaterial
   // "x0"   - the radiation length (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);
-  d*=TMath::Sqrt((1.+ fP3*fP3)/((1.-fP2)*(1.+fP2)));
-
-  //Multiple scattering******************
-  Double_t cC22 = 0.;
-  Double_t cC33 = 0.;
-  Double_t cC43 = 0.;
-  Double_t cC44 = 0.;
-  if (d!=0) {
-     Double_t theta2=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(d);
-     //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*TMath::Abs(d)*9.36*2.33;
-     if(theta2>TMath::Pi()*TMath::Pi()) return kFALSE;
-     cC22 = theta2*(1.-fP2)*(1.+fP2)*(1. + fP3*fP3);
-     cC33 = theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
-     cC43 = theta2*fP3*fP4*(1. + fP3*fP3);
-     cC44 = theta2*fP3*fP4*fP3*fP4;
-  }
 
-  //Energy losses************************
-  Double_t cP4=1.;
-  if (x0!=0. && beta2<1) {
-     d*=x0;
-     Double_t dE=Bethe(p/mass)*d;
-     Double_t e=TMath::Sqrt(p2 + mass*mass);
-     if ( TMath::Abs(dE) > 0.3*e ) return kFALSE; //30% energy loss is too much!
-     cP4 = (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)); 
-     cC44 += ((sigmadE*e/p2*fP4)*(sigmadE*e/p2*fP4)); 
-  }
+  return CorrectForMeanMaterial(d,x0*d,mass,kTRUE,Bethe);
 
-  fC22 += cC22;
-  fC33 += cC33;
-  fC43 += cC43;
-  fC44 += cC44;
-  fP4  *= cP4;
-
-  return kTRUE;
 }
 
 Double_t AliExternalTrackParam::BetheBlochAleph(Double_t bg,
@@ -639,6 +657,8 @@ Bool_t AliExternalTrackParam::Rotate(Double_t alpha) {
   fC40 *= ca;
   fC42 *= rr;
 
+  CheckCovariance();
+
   return kTRUE;
 }
 
@@ -652,7 +672,8 @@ Bool_t AliExternalTrackParam::PropagateTo(Double_t xk, Double_t b) {
   Double_t crv=GetC(b);
   if (TMath::Abs(b) < kAlmost0Field) crv=0.;
 
-  Double_t f1=fP[2], f2=f1 + crv*dx;
+  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;
 
@@ -667,9 +688,23 @@ Bool_t AliExternalTrackParam::PropagateTo(Double_t xk, Double_t b) {
   Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
 
   fX=xk;
-  fP0 += dx*(f1+f2)/(r1+r2);
-  fP1 += dx*(r2 + f2*(f1+f2)/(r1+r2))*fP3;  // Many thanks to P.Hristov !
-  fP2 += dx*crv;
+  double dy2dx = (f1+f2)/(r1+r2);
+  fP0 += dx*dy2dx;
+  if (TMath::Abs(x2r)<0.05) {
+    fP1 += dx*(r2 + f2*dy2dx)*fP3;  // Many thanks to P.Hristov !
+    fP2 += x2r;
+  }
+  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
+    fP1 += rot/crv*fP3;
+    fP2  = TMath::Sin(rot + TMath::ASin(fP2));
+  }
 
   //f = F - 1
    
@@ -711,6 +746,8 @@ Bool_t AliExternalTrackParam::PropagateTo(Double_t xk, Double_t b) {
   fC32 += b32;
   fC42 += b42;
 
+  CheckCovariance();
+
   return kTRUE;
 }
 
@@ -1107,6 +1144,8 @@ Bool_t AliExternalTrackParam::Update(Double_t p[2], Double_t cov[3]) {
 
   fC44-=k40*c04+k41*c14; 
 
+  CheckCovariance();
+
   return kTRUE;
 }
 
@@ -1387,7 +1426,6 @@ Double_t b[3], Double_t maxd, Double_t dz[2], Double_t covar[3]) {
   return kTRUE;
 }
 
-
 void AliExternalTrackParam::GetDirection(Double_t d[3]) const {
   //----------------------------------------------------------------
   // This function returns a unit vector along the track direction
@@ -1435,18 +1473,6 @@ Double_t AliExternalTrackParam::Py() const {
   return p[1];
 }
 
-Double_t AliExternalTrackParam::Pz() const {
-  //---------------------------------------------------------------------
-  // Returns z-component of momentum
-  // Result for (nearly) straight tracks is meaningless !
-  //---------------------------------------------------------------------
-
-  Double_t p[3]={kVeryBig,kVeryBig,kVeryBig};
-  GetPxPyPz(p);
-
-  return p[2];
-}
-
 Double_t AliExternalTrackParam::Xv() const {
   //---------------------------------------------------------------------
   // Returns x-component of first track point
@@ -1469,17 +1495,6 @@ Double_t AliExternalTrackParam::Yv() const {
   return r[1];
 }
 
-Double_t AliExternalTrackParam::Zv() const {
-  //---------------------------------------------------------------------
-  // Returns z-component of first track point
-  //---------------------------------------------------------------------
-
-  Double_t r[3]={0.,0.,0.};
-  GetXYZ(r);
-
-  return r[2];
-}
-
 Double_t AliExternalTrackParam::Theta() const {
   // return theta angle of momentum
 
@@ -1814,7 +1829,7 @@ void AliExternalTrackParam::g3helx3(Double_t qfield,
   Double_t tet = rho*step;
 
   Double_t tsint, sintt, sint, cos1t; 
-  if (TMath::Abs(tet) > 0.15) {
+  if (TMath::Abs(tet) > 0.03) {
      sint  = TMath::Sin(tet);
      sintt = sint/tet;
      tsint = (tet - sint)/tet;
@@ -1852,11 +1867,13 @@ Bool_t AliExternalTrackParam::PropagateToBxByBz(Double_t xk, const Double_t b[3]
 
   Double_t dx=xk-fX;
   if (TMath::Abs(dx)<=kAlmost0)  return kTRUE;
+  if (TMath::Abs(fP[4])<=kAlmost0) return kFALSE;
 
   Double_t crv=GetC(b[2]);
   if (TMath::Abs(b[2]) < kAlmost0Field) crv=0.;
 
-  Double_t f1=fP[2], f2=f1 + crv*dx;
+  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;
 
@@ -1911,12 +1928,14 @@ Bool_t AliExternalTrackParam::PropagateToBxByBz(Double_t xk, const Double_t b[3]
   fC32 += b32;
   fC42 += b42;
 
+  CheckCovariance();
   
   // Appoximate step length
-  Double_t step=dx*TMath::Abs(r2 + f2*(f1+f2)/(r1+r2));
+  double dy2dx = (f1+f2)/(r1+r2);
+  Double_t step = (TMath::Abs(x2r)<0.05) ? dx*TMath::Abs(r2 + f2*dy2dx)  // chord
+    : 2.*TMath::ASin(0.5*dx*TMath::Sqrt(1.+dy2dx*dy2dx)*crv)/crv;        // arc
   step *= TMath::Sqrt(1.+ GetTgl()*GetTgl());
 
-
   // Get the track's (x,y,z) and (px,py,pz) in the Global System
   Double_t r[3]; GetXYZ(r);
   Double_t p[3]; GetPxPyPz(p);
@@ -2076,3 +2095,74 @@ Bool_t AliExternalTrackParam::Translate(Double_t *vTrasl,Double_t *covV){
 
   return kTRUE;
  }
+
+void AliExternalTrackParam::CheckCovariance() {
+
+  // This function forces the diagonal elements of the covariance matrix to be positive.
+  // In case the diagonal element is bigger than the maximal allowed value, it is set to
+  // the limit and the off-diagonal elements that correspond to it are set to zero.
+
+    fC[0] = TMath::Abs(fC[0]);
+    if (fC[0]>kC0max) {
+      fC[0] = kC0max;
+      fC[1] = 0;
+      fC[3] = 0;
+      fC[6] = 0;
+      fC[10] = 0;
+    }
+    fC[2] = TMath::Abs(fC[2]);
+    if (fC[2]>kC2max) {
+      fC[2] = kC2max;
+      fC[1] = 0;
+      fC[4] = 0;
+      fC[7] = 0;
+      fC[11] = 0;
+    }
+    fC[5] = TMath::Abs(fC[5]);
+    if (fC[5]>kC5max) {
+      fC[5] = kC5max;
+      fC[3] = 0;
+      fC[4] = 0;
+      fC[8] = 0;
+      fC[12] = 0;
+    }
+    fC[9] = TMath::Abs(fC[9]);
+    if (fC[9]>kC9max) {
+      fC[9] = kC9max;
+      fC[6] = 0;
+      fC[7] = 0;
+      fC[8] = 0;
+      fC[13] = 0;
+    }
+    fC[14] = TMath::Abs(fC[14]);
+    if (fC[14]>kC14max) {
+      fC[14] = kC14max;
+      fC[10] = 0;
+      fC[11] = 0;
+      fC[12] = 0;
+      fC[13] = 0;
+    }
+    
+    // The part below is used for tests and normally is commented out    
+//     TMatrixDSym m(5);
+//     TVectorD eig(5);
+    
+//     m(0,0)=fC[0];
+//     m(1,0)=fC[1];  m(1,1)=fC[2];
+//     m(2,0)=fC[3];  m(2,1)=fC[4];  m(2,2)=fC[5];
+//     m(3,0)=fC[6];  m(3,1)=fC[7];  m(3,2)=fC[8];  m(3,3)=fC[9];
+//     m(4,0)=fC[10]; m(4,1)=fC[11]; m(4,2)=fC[12]; m(4,3)=fC[13]; m(4,4)=fC[14];
+    
+//     m(0,1)=m(1,0);
+//     m(0,2)=m(2,0); m(1,2)=m(2,1);
+//     m(0,3)=m(3,0); m(1,3)=m(3,1); m(2,3)=m(3,2);
+//     m(0,4)=m(4,0); m(1,4)=m(4,1); m(2,4)=m(4,2); m(3,4)=m(4,3);
+//     m.EigenVectors(eig);
+
+//     //    assert(eig(0)>=0 && eig(1)>=0 && eig(2)>=0 && eig(3)>=0 && eig(4)>=0);
+//     if (!(eig(0)>=0 && eig(1)>=0 && eig(2)>=0 && eig(3)>=0 && eig(4)>=0)) {
+//       AliWarning("Negative eigenvalues of the covariance matrix!");
+//       this->Print();
+//       eig.Print();
+//     }
+}