// 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(),
return -d;
}
-Bool_t AliExternalTrackParam::
-CorrectForMaterial(Double_t d, Double_t x0, Double_t mass) {
+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)
//Energy losses************************
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)
- dE=0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2)*d;
-
- fP4*=(1.- TMath::Sqrt(p2 + mass*mass)/p2*dE);
+ 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 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;
}
+Double_t ApproximateBetheBloch(Double_t beta2) {
+ //------------------------------------------------------------------
+ // This is an approximation of the Bethe-Bloch formula with
+ // the density effect taken into account at beta*gamma > 3.5
+ // (the approximation is reasonable only for solid materials)
+ //------------------------------------------------------------------
+ if (beta2/(1-beta2)>3.5*3.5)
+ return 0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2);
+
+ return 0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2);
+}
+
Bool_t AliExternalTrackParam::Rotate(Double_t alpha) {
//------------------------------------------------------------------
// Transform this track to the local coord. system rotated
return kTRUE;
}
+void AliExternalTrackParam::Propagate(Double_t len, Double_t x[3],
+Double_t p[3], Double_t bz) const {
+ //+++++++++++++++++++++++++++++++++++++++++
+ // Origin: K. Shileev (Kirill.Shileev@cern.ch)
+ // Extrapolate track along simple helix in magnetic field
+ // Arguments: len -distance alogn helix, [cm]
+ // bz - mag field, [kGaus]
+ // Returns: x and p contain extrapolated positon and momentum
+ // The momentum returned for straight-line tracks is meaningless !
+ //+++++++++++++++++++++++++++++++++++++++++
+ GetXYZ(x);
+
+ 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;
+ x[2]+=unit[2]*len;
+
+ p[0]=unit[0]/kAlmost0;
+ p[1]=unit[1]/kAlmost0;
+ p[2]=unit[2]/kAlmost0;
+ } else {
+ GetPxPyPz(p);
+ Double_t pp=GetP();
+ Double_t a = -kB2C*bz*GetSign();
+ Double_t rho = a/pp;
+ x[0] += p[0]*TMath::Sin(rho*len)/a - p[1]*(1-TMath::Cos(rho*len))/a;
+ x[1] += p[1]*TMath::Sin(rho*len)/a + p[0]*(1-TMath::Cos(rho*len))/a;
+ x[2] += p[2]*len/pp;
+
+ Double_t p0=p[0];
+ p[0] = p0 *TMath::Cos(rho*len) - p[1]*TMath::Sin(rho*len);
+ p[1] = p[1]*TMath::Cos(rho*len) + p0 *TMath::Sin(rho*len);
+ }
+}
+
+Bool_t AliExternalTrackParam::Intersect(Double_t pnt[3], Double_t norm[3],
+Double_t bz) const {
+ //+++++++++++++++++++++++++++++++++++++++++
+ // Origin: K. Shileev (Kirill.Shileev@cern.ch)
+ // Finds point of intersection (if exists) of the helix with the plane.
+ // Stores result in fX and fP.
+ // Arguments: planePoint,planeNorm - the plane defined by any plane's point
+ // and vector, normal to the plane
+ // Returns: kTrue if helix intersects the plane, kFALSE otherwise.
+ //+++++++++++++++++++++++++++++++++++++++++
+ Double_t x0[3]; GetXYZ(x0); //get track position in MARS
+
+ //estimates initial helix length up to plane
+ Double_t s=
+ (pnt[0]-x0[0])*norm[0] + (pnt[1]-x0[1])*norm[1] + (pnt[2]-x0[2])*norm[2];
+ Double_t dist=99999,distPrev=dist;
+ Double_t x[3],p[3];
+ while(TMath::Abs(dist)>0.00001){
+ //calculates helix at the distance s from x0 ALONG the helix
+ Propagate(s,x,p,bz);
+
+ //distance between current helix position and plane
+ dist=(x[0]-pnt[0])*norm[0]+(x[1]-pnt[1])*norm[1]+(x[2]-pnt[2])*norm[2];
+
+ if(TMath::Abs(dist) >= TMath::Abs(distPrev)) {return kFALSE;}
+ distPrev=dist;
+ s-=dist;
+ }
+ //on exit pnt is intersection point,norm is track vector at that point,
+ //all in MARS
+ for (Int_t i=0; i<3; i++) {pnt[i]=x[i]; norm[i]=p[i];}
+ return kTRUE;
+}
+
Double_t
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
//----------------------------------------------------------------
Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
Double_t snp=fP[2];
- Double_t csp =TMath::Sqrt(1.- snp*snp);
+ Double_t csp =TMath::Sqrt((1.- snp)*(1.+snp));
Double_t norm=TMath::Sqrt(1.+ fP[3]*fP[3]);
d[0]=(csp*cs - snp*sn)/norm;
d[1]=(snp*cs + csp*sn)/norm;
}
Double_t pt=1./TMath::Abs(fP[4]);
Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
- Double_t r=TMath::Sqrt(1-fP[2]*fP[2]);
+ Double_t r=TMath::Sqrt((1.-fP[2])*(1.+fP[2]));
Double_t m00=-sn, m10=cs;
Double_t m23=-pt*(sn + fP[2]*cs/r), m43=-pt*pt*(r*cs - fP[2]*sn);
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;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff