// are implemented.
// Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch //
///////////////////////////////////////////////////////////////////////////////
+#include <TMatrixDSym.h>
#include "AliExternalTrackParam.h"
-#include "AliKalmanTrack.h"
#include "AliESDVertex.h"
-
+#include "AliLog.h"
ClassImp(AliExternalTrackParam)
+Double32_t AliExternalTrackParam::fgMostProbablePt=kMostProbablePt;
+
//_____________________________________________________________________________
AliExternalTrackParam::AliExternalTrackParam() :
+ TObject(),
fX(0),
fAlpha(0)
{
for (Int_t i = 0; i < 15; i++) fC[i] = 0;
}
+//_____________________________________________________________________________
+AliExternalTrackParam::AliExternalTrackParam(const AliExternalTrackParam &track):
+ TObject(track),
+ fX(track.fX),
+ fAlpha(track.fAlpha)
+{
+ //
+ // copy constructor
+ //
+ 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];
+}
+
//_____________________________________________________________________________
AliExternalTrackParam::AliExternalTrackParam(Double_t x, Double_t alpha,
const Double_t param[5],
const Double_t covar[15]) :
+ TObject(),
fX(x),
fAlpha(alpha)
{
}
//_____________________________________________________________________________
-AliExternalTrackParam::AliExternalTrackParam(const AliKalmanTrack& track) :
- fX(0),
- fAlpha(track.GetAlpha())
-{
- //
- //
- track.GetExternalParameters(fX,fP);
- track.GetExternalCovariance(fC);
-}
-
-//_____________________________________________________________________________
-void AliExternalTrackParam::Set(const AliKalmanTrack& track) {
+void AliExternalTrackParam::Set(Double_t x, Double_t alpha,
+ const Double_t p[5], const Double_t cov[15]) {
//
+ // Sets the parameters
//
- fAlpha=track.GetAlpha();
- track.GetExternalParameters(fX,fP);
- track.GetExternalCovariance(fC);
+ fX=x;
+ fAlpha=alpha;
+ for (Int_t i = 0; i < 5; i++) fP[i] = p[i];
+ for (Int_t i = 0; i < 15; i++) fC[i] = cov[i];
}
//_____________________________________________________________________________
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)
Double_t &fC43=fC[13];
Double_t &fC44=fC[14];
- Double_t p2=(1.+ fP3*fP3)/(fP4*fP4);
+ Double_t p=GetP();
+ Double_t p2=p*p;
Double_t beta2=p2/(p2 + mass*mass);
d*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
//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
return kTRUE;
}
+void AliExternalTrackParam::GetDirection(Double_t d[3]) const {
+ //----------------------------------------------------------------
+ // This function returns a unit vector along the track direction
+ // in the global coordinate system.
+ //----------------------------------------------------------------
+ Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
+ Double_t snp=fP[2];
+ 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;
+ d[2]=fP[3]/norm;
+}
+
Bool_t AliExternalTrackParam::GetPxPyPz(Double_t *p) const {
//---------------------------------------------------------------------
// This function returns the global track momentum components
}
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;
return kTRUE;
}
+Bool_t
+AliExternalTrackParam::GetZAt(Double_t x, Double_t b, Double_t &z) const {
+ //---------------------------------------------------------------------
+ // This function returns the local Z-coordinate of the intersection
+ // point between this track and the reference plane "x" (cm).
+ // Magnetic field "b" (kG)
+ //---------------------------------------------------------------------
+ Double_t dx=x-fX;
+ if(TMath::Abs(dx)<=kAlmost0) {z=fP[1]; return kTRUE;}
+
+ 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=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
+ z = fP[1] + dx*(r2 + f2*(f1+f2)/(r1+r2))*fP[3]; // Many thanks to P.Hristov !
+ return kTRUE;
+}
+
Bool_t
AliExternalTrackParam::GetXYZAt(Double_t x, Double_t b, Double_t *r) const {
//---------------------------------------------------------------------
Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
r[0] = x;
r[1] = fP[0] + dx*(f1+f2)/(r1+r2);
- r[2] = fP[1] + dx*(r2 + f2*(f1+f2)/(r1+r2))*fP[3];//Many thanks to P.Hristov !
+ r[2] = fP[1] + dx*(f1+f2)/(f1*r2 + f2*r1)*fP[3];
return Local2GlobalPosition(r,fAlpha);
}
Double_t res = fP[2]+dx*crv;
return res;
}
-
-Bool_t AliExternalTrackParam::PropagateTo(Double_t xToGo, Double_t b, Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp){
- //----------------------------------------------------------------
- //
- // 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
- //
- // Requires acces to geomanager
- //
- // mass - mass used in propagation - used for energy loss correction
- // maxStep - maximal step for propagation
- //----------------------------------------------------------------
- const Double_t kEpsilon = 0.00001;
- Double_t xpos = GetX();
- Double_t dir = (xpos<xToGo) ? 1.:-1.;
- //
- while ( (xToGo-xpos)*dir > kEpsilon){
- 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
- if (!GetXYZAt(x,b,xyz1)) return kFALSE; // no prolongation
- xyz1[2]+=kEpsilon; // waiting for bug correction in geo
- AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
- if (TMath::Abs(GetSnpAt(x,b)) >= maxSnp) 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]) >= maxSnp) return kFALSE;
- GetXYZ(xyz0); // global position
- Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
- //
- Double_t ca=TMath::Cos(alphan-fAlpha), sa=TMath::Sin(alphan-fAlpha);
- Double_t sf=fP[2], cf=TMath::Sqrt(1.- fP[2]*fP[2]);
- Double_t sinNew = sf*ca - cf*sa;
- if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
- if (!Rotate(alphan)) return kFALSE;
- }
- xpos = GetX();
- }
- return kTRUE;
-}
-
-
-
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff