]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/AliExternalTrackParam.cxx
First V0 MC Analysis from H.Ricaud
[u/mrichter/AliRoot.git] / STEER / AliExternalTrackParam.cxx
index 651ae7e80b162350b95a8bd41ace96f2c3a6994f..6f3906916a0294ff6e2459ccad8a1ee8cb2f8223 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(),
+  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())<kAlmost0){
+    xyz2[0]=x;
+    if (!param2->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;
+}