X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliExternalTrackParam.cxx;h=6f3906916a0294ff6e2459ccad8a1ee8cb2f8223;hb=04cd9bcd360d53ac0d1772d679c603d51201ed54;hp=651ae7e80b162350b95a8bd41ace96f2c3a6994f;hpb=59c0deeafca58240227ef1193d01cc7beb7cd5c5;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliExternalTrackParam.cxx b/STEER/AliExternalTrackParam.cxx index 651ae7e80b1..6f3906916a0 100644 --- a/STEER/AliExternalTrackParam.cxx +++ b/STEER/AliExternalTrackParam.cxx @@ -25,15 +25,18 @@ // are implemented. // Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch // /////////////////////////////////////////////////////////////////////////////// +#include #include "AliExternalTrackParam.h" #include "AliESDVertex.h" #include "AliLog.h" ClassImp(AliExternalTrackParam) +Double32_t AliExternalTrackParam::fgMostProbablePt=kMostProbablePt; + //_____________________________________________________________________________ AliExternalTrackParam::AliExternalTrackParam() : - TObject(), + AliVParticle(), fX(0), fAlpha(0) { @@ -46,7 +49,7 @@ AliExternalTrackParam::AliExternalTrackParam() : //_____________________________________________________________________________ AliExternalTrackParam::AliExternalTrackParam(const AliExternalTrackParam &track): - TObject(track), + AliVParticle(track), fX(track.fX), fAlpha(track.fAlpha) { @@ -57,11 +60,30 @@ AliExternalTrackParam::AliExternalTrackParam(const AliExternalTrackParam &track) for (Int_t i = 0; i < 15; i++) fC[i] = track.fC[i]; } +//_____________________________________________________________________________ +AliExternalTrackParam& AliExternalTrackParam::operator=(const AliExternalTrackParam &trkPar) +{ + // + // assignment operator + // + + if (this!=&trkPar) { + AliVParticle::operator=(trkPar); + fX = trkPar.fX; + fAlpha = trkPar.fAlpha; + + 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]; + } + + return *this; +} + //_____________________________________________________________________________ AliExternalTrackParam::AliExternalTrackParam(Double_t x, Double_t alpha, const Double_t param[5], const Double_t covar[15]) : - TObject(), + AliVParticle(), fX(x), fAlpha(alpha) { @@ -178,9 +200,69 @@ 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)) { + //------------------------------------------------------------------ + // 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]; + + //Apply angle correction, if requested + if(anglecorr) { + Double_t angle=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2)); + xOverX0 *=angle; + xTimesRho *=angle; + } + + Double_t p=GetP(); + Double_t p2=p*p; + Double_t beta2=p2/(p2 + mass*mass); + + //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) @@ -219,8 +301,8 @@ Bool_t AliExternalTrackParam::CorrectForMaterial 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)); } @@ -381,7 +463,7 @@ Double_t p[3], Double_t bz) const { //+++++++++++++++++++++++++++++++++++++++++ GetXYZ(x); - if (TMath::Abs(Get1Pt()) < kAlmost0 || TMath::Abs(bz) < kAlmost0Field ){ //straight-line tracks + if (OneOverPt() < kAlmost0 || TMath::Abs(bz) < kAlmost0Field ){ //straight-line tracks Double_t unit[3]; GetDirection(unit); x[0]+=unit[0]*len; x[1]+=unit[1]*len; @@ -457,6 +539,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 @@ -716,54 +906,6 @@ Bool_t AliExternalTrackParam::PropagateToDCA(const AliESDVertex *vtx, Double_t b } - - -Bool_t Local2GlobalMomentum(Double_t p[3],Double_t alpha) { - //---------------------------------------------------------------- - // This function performs local->global transformation of the - // track momentum. - // When called, the arguments are: - // p[0] = 1/pt of the track; - // p[1] = sine of local azim. angle of the track momentum; - // p[2] = tangent of the track momentum dip angle; - // alpha - rotation angle. - // The result is returned as: - // p[0] = px - // p[1] = py - // p[2] = pz - // Results for (nearly) straight tracks are meaningless ! - //---------------------------------------------------------------- - if (TMath::Abs(p[0])<=kAlmost0) return kFALSE; - if (TMath::Abs(p[1])> kAlmost1) return kFALSE; - - Double_t pt=1./TMath::Abs(p[0]); - Double_t cs=TMath::Cos(alpha), sn=TMath::Sin(alpha); - Double_t r=TMath::Sqrt(1 - p[1]*p[1]); - p[0]=pt*(r*cs - p[1]*sn); p[1]=pt*(p[1]*cs + r*sn); p[2]=pt*p[2]; - - return kTRUE; -} - -Bool_t Local2GlobalPosition(Double_t r[3],Double_t alpha) { - //---------------------------------------------------------------- - // This function performs local->global transformation of the - // track position. - // When called, the arguments are: - // r[0] = local x - // r[1] = local y - // r[2] = local z - // alpha - rotation angle. - // The result is returned as: - // r[0] = global x - // r[1] = global y - // r[2] = global z - //---------------------------------------------------------------- - Double_t cs=TMath::Cos(alpha), sn=TMath::Sin(alpha), x=r[0]; - r[0]=x*cs - r[1]*sn; r[1]=x*sn + r[1]*cs; - - return kTRUE; -} - void AliExternalTrackParam::GetDirection(Double_t d[3]) const { //---------------------------------------------------------------- // This function returns a unit vector along the track direction @@ -778,7 +920,7 @@ void AliExternalTrackParam::GetDirection(Double_t d[3]) const { d[2]=fP[3]/norm; } -Bool_t AliExternalTrackParam::GetPxPyPz(Double_t *p) const { +Bool_t AliExternalTrackParam::GetPxPyPz(Double_t p[3]) const { //--------------------------------------------------------------------- // This function returns the global track momentum components // Results for (nearly) straight tracks are meaningless ! @@ -787,6 +929,127 @@ Bool_t AliExternalTrackParam::GetPxPyPz(Double_t *p) const { return Local2GlobalMomentum(p,fAlpha); } +Double_t AliExternalTrackParam::Px() const { + //--------------------------------------------------------------------- + // Returns x-component of momentum + // Result for (nearly) straight tracks is meaningless ! + //--------------------------------------------------------------------- + + Double_t p[3]={kVeryBig,kVeryBig,kVeryBig}; + GetPxPyPz(p); + + return p[0]; +} + +Double_t AliExternalTrackParam::Py() const { + //--------------------------------------------------------------------- + // Returns y-component of momentum + // Result for (nearly) straight tracks is meaningless ! + //--------------------------------------------------------------------- + + Double_t p[3]={kVeryBig,kVeryBig,kVeryBig}; + GetPxPyPz(p); + + 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 + //--------------------------------------------------------------------- + + Double_t r[3]={0.,0.,0.}; + GetXYZ(r); + + return r[0]; +} + +Double_t AliExternalTrackParam::Yv() const { + //--------------------------------------------------------------------- + // Returns y-component of first track point + //--------------------------------------------------------------------- + + Double_t r[3]={0.,0.,0.}; + GetXYZ(r); + + 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 + + return 0.5*TMath::Pi() - TMath::ATan(fP[3]); +} + +Double_t AliExternalTrackParam::Phi() const { + //--------------------------------------------------------------------- + // Returns the azimuthal angle of momentum + // 0 <= phi < 2*pi + //--------------------------------------------------------------------- + + Double_t phi=TMath::ASin(fP[2]) + fAlpha; + if (phi<0.) phi+=2.*TMath::Pi(); + else if (phi>=2.*TMath::Pi()) phi-=2.*TMath::Pi(); + + return phi; +} + +Double_t AliExternalTrackParam::M() const { + // return particle mass + + // No mass information available so far. + // Redifine in derived class! + + return -999.; +} + +Double_t AliExternalTrackParam::E() const { + // return particle energy + + // No PID information available so far. + // Redifine in derived class! + + return -999.; +} + +Double_t AliExternalTrackParam::Eta() const { + // return pseudorapidity + + return -TMath::Log(TMath::Tan(0.5 * Theta())); +} + +Double_t AliExternalTrackParam::Y() const { + // return rapidity + + // No PID information available so far. + // Redifine in derived class! + + return -999.; +} + Bool_t AliExternalTrackParam::GetXYZ(Double_t *r) const { //--------------------------------------------------------------------- // This function returns the global track position @@ -825,6 +1088,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; @@ -951,3 +1218,43 @@ Double_t AliExternalTrackParam::GetSnpAt(Double_t x,Double_t b) const { Double_t res = fP[2]+dx*crv; return res; } + +Bool_t AliExternalTrackParam::GetDistance(AliExternalTrackParam *param2, Double_t x, Double_t dist[3], Double_t bz){ + //------------------------------------------------------------------------ + // Get the distance between two tracks at the local position x + // working in the local frame of this track. + // Origin : Marian.Ivanov@cern.ch + //----------------------------------------------------------------------- + Double_t xyz[3]; + Double_t xyz2[3]; + xyz[0]=x; + if (!GetYAt(x,bz,xyz[1])) return kFALSE; + if (!GetZAt(x,bz,xyz[2])) return kFALSE; + // + // + if (TMath::Abs(GetAlpha()-param2->GetAlpha())GetYAt(x,bz,xyz2[1])) return kFALSE; + if (!param2->GetZAt(x,bz,xyz2[2])) return kFALSE; + }else{ + // + Double_t xyz1[3]; + Double_t dfi = param2->GetAlpha()-GetAlpha(); + Double_t ca = TMath::Cos(dfi), sa = TMath::Sin(dfi); + xyz2[0] = xyz[0]*ca+xyz[1]*sa; + xyz2[1] = -xyz[0]*sa+xyz[1]*ca; + // + xyz1[0]=xyz2[0]; + if (!param2->GetYAt(xyz2[0],bz,xyz1[1])) return kFALSE; + if (!param2->GetZAt(xyz2[0],bz,xyz1[2])) return kFALSE; + // + xyz2[0] = xyz1[0]*ca-xyz1[1]*sa; + xyz2[1] = +xyz1[0]*sa+xyz1[1]*ca; + xyz2[2] = xyz1[2]; + } + dist[0] = xyz[0]-xyz2[0]; + dist[1] = xyz[1]-xyz2[1]; + dist[2] = xyz[2]-xyz2[2]; + + return kTRUE; +}