X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=STEER%2FAliKalmanTrack.cxx;h=92df4af77cfd24124241509e795ccc97cfd820f3;hp=b01f2d5b649c448e7c44785e343e2342b921e47b;hb=96c33f8f31e546eab077ada93d9c3fc663d36ce4;hpb=f37d970d39405521c25e9753ee52970a06cd051b

diff --git a/STEER/AliKalmanTrack.cxx b/STEER/AliKalmanTrack.cxx
index b01f2d5b649..92df4af77cf 100644
--- a/STEER/AliKalmanTrack.cxx
+++ b/STEER/AliKalmanTrack.cxx
@@ -20,242 +20,62 @@
 //   that is the base for AliTPCtrack, AliITStrackV2 and AliTRDtrack
 //        Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
 //-------------------------------------------------------------------------
+#include <TGeoManager.h>
 
 #include "AliKalmanTrack.h"
-#include "AliLog.h"
-#include "AliPDG.h"
-#include "TPDGCode.h"
-#include "TDatabasePDG.h"
 
 ClassImp(AliKalmanTrack)
 
-Double_t AliKalmanTrack::fgConvConst;
-
 //_______________________________________________________________________
-AliKalmanTrack::AliKalmanTrack():
-  fLab(-3141593),
+  AliKalmanTrack::AliKalmanTrack():AliExternalTrackParam(),
+  fFakeRatio(0),
   fChi2(0),
-  fMass(0.13957),
-  fN(0)
+  fMass(AliPID::ParticleMass(AliPID::kPion)),
+  fLab(-3141593),
+  fN(0),
+  fStartTimeIntegral(kFALSE),
+  fIntegratedLength(0)
 {
   //
   // Default constructor
   //
-    if (fgConvConst==0) {
-      AliFatal("The magnetic field has not been set!");
-    }
-    
-    fStartTimeIntegral = kFALSE;
-    fIntegratedLength = 0;
-    for(Int_t i=0; i<5; i++) fIntegratedTime[i] = 0;
+
+  for(Int_t i=0; i<AliPID::kSPECIES; i++) fIntegratedTime[i] = 0;
 }
 
-//_______________________________________________________________________
 AliKalmanTrack::AliKalmanTrack(const AliKalmanTrack &t):
-  TObject(t),
-  fLab(t.fLab),
+  AliExternalTrackParam(t),
   fFakeRatio(t.fFakeRatio),
   fChi2(t.fChi2),
   fMass(t.fMass),
-  fN(t.fN)
+  fLab(t.fLab),
+  fN(t.fN),
+  fStartTimeIntegral(t.fStartTimeIntegral),
+  fIntegratedLength(t.fIntegratedLength)
 {
   //
   // Copy constructor
   //
-  if (fgConvConst==0) {
-    AliFatal("The magnetic field has not been set!");
-  }
-
-  fStartTimeIntegral = t.fStartTimeIntegral;
-  fIntegratedLength = t.fIntegratedLength;
   
-  for (Int_t i=0; i<5; i++) 
-    fIntegratedTime[i] = t.fIntegratedTime[i];
-}
-
-//_______________________________________________________________________
-Double_t AliKalmanTrack::GetX() const
-{
-  // Returns the X coordinate of the current track position
-  AliWarning("Method must be overloaded !");
-  return 0.;
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::GetdEdx() const
-{
-  // Returns the dE/dx of the track
-  AliWarning("Method must be overloaded !");
-  return 0.;
-}
-
-//_______________________________________________________________________
-Double_t AliKalmanTrack::GetY() const
-{
-  // Returns the Y coordinate of the current track position
-  Double_t par[5];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  return par[0];
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::GetZ() const
-{
-  // Returns the Z coordinate of the current track position
-  Double_t par[5];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  return par[1];
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::GetSnp() const
-{
-  // Returns the Sin(phi), where phi is the angle between the transverse
-  // momentum (in xOy plane) and the X axis
-  Double_t par[5];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  return par[2];
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::GetTgl() const
-{
-  // Returns the Tan(lambda), where lambda is the dip angle (between 
-  // the bending plane (xOy) and the momentum of the track
-  Double_t par[5];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  return par[3];
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::Get1Pt() const
-{
-  // Returns 1/pT
-  Double_t par[5];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  return par[4];
-}
-
-//_______________________________________________________________________
-Double_t AliKalmanTrack::Phi() const
-{
-// return global phi of track
-
-  Double_t par[5];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  if (par[2] >  1.) par[2] =  1.;
-  if (par[2] < -1.) par[2] = -1.;
-  Double_t phi = TMath::ASin(par[2]) + GetAlpha();
-  while (phi < 0) phi += TMath::TwoPi();
-  while (phi > TMath::TwoPi()) phi -= TMath::TwoPi();
-  return phi;
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::SigmaPhi() const
-{
-// return error of global phi of track
-
-  Double_t par[5];
-  Double_t cov[15];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  GetExternalCovariance(cov);
-  return TMath::Sqrt(TMath::Abs(cov[5] / (1. - par[2]*par[2])));
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::Theta() const
-{
-// return global theta of track
-
-  Double_t par[5];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  return TMath::Pi()/2. - TMath::ATan(par[3]);
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::SigmaTheta() const
-{
-// return error of global theta of track
-
-  Double_t par[5];
-  Double_t cov[15];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  GetExternalCovariance(cov);
-  return TMath::Sqrt(TMath::Abs(cov[5])) / (1. + par[3]*par[3]);
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::Eta() const
-{
-// return global eta of track
-
-  return -TMath::Log(TMath::Tan(Theta()/2.));
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::Px() const
-{
-// return x component of track momentum
-
-  Double_t par[5];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  Double_t phi = TMath::ASin(par[2]) + GetAlpha();
-  return TMath::Cos(phi) / TMath::Abs(par[4]);
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::Py() const
-{
-// return y component of track momentum
-
-  Double_t par[5];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  Double_t phi = TMath::ASin(par[2]) + GetAlpha();
-  return TMath::Sin(phi) / TMath::Abs(par[4]);
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::Pz() const
-{
-// return z component of track momentum
-
-  Double_t par[5];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  return par[3] / TMath::Abs(par[4]);
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::Pt() const
-{
-// return transverse component of track momentum
-
-  Double_t par[5];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  return 1. / TMath::Abs(par[4]);
-}
-//_______________________________________________________________________
-Double_t AliKalmanTrack::SigmaPt() const
-{
-// return error of transverse component of track momentum
-
-  Double_t par[5];
-  Double_t cov[15];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  GetExternalCovariance(cov);
-  return TMath::Sqrt(cov[14]) / TMath::Abs(par[4]);
+  for (Int_t i=0; i<AliPID::kSPECIES; i++)
+      fIntegratedTime[i] = t.fIntegratedTime[i];
+}
+
+AliKalmanTrack& AliKalmanTrack::operator=(const AliKalmanTrack&o){
+  if(this!=&o){
+    AliExternalTrackParam::operator=(o);
+    fLab = o.fLab;
+    fFakeRatio = o.fFakeRatio;
+    fChi2 = o.fChi2;
+    fMass = o.fMass;
+    fN = o.fN;
+    fStartTimeIntegral = o.fStartTimeIntegral;
+    for(Int_t i = 0;i<AliPID::kSPECIES;++i)fIntegratedTime[i] = o.fIntegratedTime[i];
+    fIntegratedLength = o.fIntegratedLength;
+  }
+  return *this;
 }
-//_______________________________________________________________________
-Double_t AliKalmanTrack::P() const
-{
-// return total track momentum
 
-  Double_t par[5];
-  Double_t localX = GetX();
-  GetExternalParameters(localX, par);
-  return 1. / TMath::Abs(par[4] * TMath::Cos(TMath::ATan(par[3])));
-}
 //_______________________________________________________________________
 void AliKalmanTrack::StartTimeIntegral() 
 {
@@ -270,9 +90,10 @@ void AliKalmanTrack::StartTimeIntegral()
   //  AliWarning("Reseting Recorded Time.");
 
   fStartTimeIntegral = kTRUE;
-  for(Int_t i=0; i<fgkTypes; i++) fIntegratedTime[i] = 0;  
+  for(Int_t i=0; i<AliPID::kSPECIES; i++) fIntegratedTime[i] = 0;  
   fIntegratedLength = 0;
 }
+
 //_______________________________________________________________________
 void AliKalmanTrack:: AddTimeStep(Double_t length) 
 {
@@ -302,9 +123,6 @@ void AliKalmanTrack:: AddTimeStep(Double_t length)
   
   fIntegratedLength += length;
 
-  static Int_t pdgCode[fgkTypes]  = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
-  TDatabasePDG *db = TDatabasePDG::Instance();
-
   Double_t xr, param[5];
   Double_t pt, tgl;
   
@@ -316,9 +134,9 @@ void AliKalmanTrack:: AddTimeStep(Double_t length)
 
   if (length > 100) return;
 
-  for (Int_t i=0; i<fgkTypes; i++) {
+  for (Int_t i=0; i<AliPID::kSPECIES; i++) {
     
-    Double_t mass = db->GetParticle(pdgCode[i])->Mass();
+    Double_t mass = AliPID::ParticleMass(i);
     Double_t correction = TMath::Sqrt( pt*pt * (1 + tgl*tgl) + mass * mass ) / p;
     Double_t time = length * correction / kcc;
 
@@ -327,7 +145,6 @@ void AliKalmanTrack:: AddTimeStep(Double_t length)
 }
 
 //_______________________________________________________________________
-
 Double_t AliKalmanTrack::GetIntegratedTime(Int_t pdg) const 
 {
   // Sylwester Radomski, GSI
@@ -346,194 +163,18 @@ Double_t AliKalmanTrack::GetIntegratedTime(Int_t pdg) const
     return 0.;
   }
 
-  static Int_t pdgCode[fgkTypes] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
-
-  for (Int_t i=0; i<fgkTypes; i++)
-    if (pdgCode[i] == TMath::Abs(pdg)) return fIntegratedTime[i];
+  for (Int_t i=0; i<AliPID::kSPECIES; i++)
+    if (AliPID::ParticleCode(i) == TMath::Abs(pdg)) return fIntegratedTime[i];
 
   AliWarning(Form("Particle type [%d] not found", pdg));
   return 0;
 }
 
 void AliKalmanTrack::GetIntegratedTimes(Double_t *times) const {
-  for (Int_t i=0; i<fgkTypes; i++) times[i]=fIntegratedTime[i];
+  for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fIntegratedTime[i];
 }
 
 void AliKalmanTrack::SetIntegratedTimes(const Double_t *times) {
-  for (Int_t i=0; i<fgkTypes; i++) fIntegratedTime[i]=times[i];
-}
-
-//_______________________________________________________________________
-
-void AliKalmanTrack::PrintTime() const
-{
-  // Sylwester Radomski, GSI
-  // S.Radomski@gsi.de
-  //
-  // For testing
-  // Prints time for all hypothesis
-  //
-
-  static Int_t pdgCode[fgkTypes] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
-
-  for (Int_t i=0; i<fgkTypes; i++)
-    printf("%d: %.2f  ", pdgCode[i], fIntegratedTime[i]);
-  printf("\n");  
-}
-
-static void External2Helix(const AliKalmanTrack *t, Double_t helix[6]) { 
-  //--------------------------------------------------------------------
-  // External track parameters -> helix parameters 
-  //--------------------------------------------------------------------
-  Double_t alpha,x,cs,sn;
-  t->GetExternalParameters(x,helix); alpha=t->GetAlpha();
-
-  cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
-  helix[5]=x*cs - helix[0]*sn;            // x0
-  helix[0]=x*sn + helix[0]*cs;            // y0
-//helix[1]=                               // z0
-  helix[2]=TMath::ASin(helix[2]) + alpha; // phi0
-//helix[3]=                               // tgl
-  helix[4]=helix[4]/t->GetConvConst();    // C
+  for (Int_t i=0; i<AliPID::kSPECIES; i++) fIntegratedTime[i]=times[i];
 }
 
-static void Evaluate(const Double_t *h, Double_t t,
-                     Double_t r[3],  //radius vector
-                     Double_t g[3],  //first defivatives
-                     Double_t gg[3]) //second derivatives
-{
-  //--------------------------------------------------------------------
-  // Calculate position of a point on a track and some derivatives
-  //--------------------------------------------------------------------
-  Double_t phase=h[4]*t+h[2];
-  Double_t sn=TMath::Sin(phase), cs=TMath::Cos(phase);
-
-  r[0] = h[5] + (sn - h[6])/h[4];
-  r[1] = h[0] - (cs - h[7])/h[4];  
-  r[2] = h[1] + h[3]*t;
-
-  g[0] = cs; g[1]=sn; g[2]=h[3];
-  
-  gg[0]=-h[4]*sn; gg[1]=h[4]*cs; gg[2]=0.;
-}
-
-Double_t AliKalmanTrack::
-GetDCA(const AliKalmanTrack *p, Double_t &xthis, Double_t &xp) const {
-  //------------------------------------------------------------
-  // Returns the (weighed !) distance of closest approach between 
-  // this track and the track passed as the argument.
-  // Other returned values:
-  //   xthis, xt - coordinates of tracks' reference planes at the DCA 
-  //-----------------------------------------------------------
-  Double_t dy2=GetSigmaY2() + p->GetSigmaY2();
-  Double_t dz2=GetSigmaZ2() + p->GetSigmaZ2();
-  Double_t dx2=dy2; 
-
-  //dx2=dy2=dz2=1.;
-
-  Double_t p1[8]; External2Helix(this,p1);
-  p1[6]=TMath::Sin(p1[2]); p1[7]=TMath::Cos(p1[2]);
-  Double_t p2[8]; External2Helix(p,p2);
-  p2[6]=TMath::Sin(p2[2]); p2[7]=TMath::Cos(p2[2]);
-
-
-  Double_t r1[3],g1[3],gg1[3]; Double_t t1=0.;
-  Evaluate(p1,t1,r1,g1,gg1);
-  Double_t r2[3],g2[3],gg2[3]; Double_t t2=0.;
-  Evaluate(p2,t2,r2,g2,gg2);
-
-  Double_t dx=r2[0]-r1[0], dy=r2[1]-r1[1], dz=r2[2]-r1[2];
-  Double_t dm=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
-
-  Int_t max=27;
-  while (max--) {
-     Double_t gt1=-(dx*g1[0]/dx2 + dy*g1[1]/dy2 + dz*g1[2]/dz2);
-     Double_t gt2=+(dx*g2[0]/dx2 + dy*g2[1]/dy2 + dz*g2[2]/dz2);
-     Double_t h11=(g1[0]*g1[0] - dx*gg1[0])/dx2 + 
-                  (g1[1]*g1[1] - dy*gg1[1])/dy2 +
-                  (g1[2]*g1[2] - dz*gg1[2])/dz2;
-     Double_t h22=(g2[0]*g2[0] + dx*gg2[0])/dx2 + 
-                  (g2[1]*g2[1] + dy*gg2[1])/dy2 +
-                  (g2[2]*g2[2] + dz*gg2[2])/dz2;
-     Double_t h12=-(g1[0]*g2[0]/dx2 + g1[1]*g2[1]/dy2 + g1[2]*g2[2]/dz2);
-
-     Double_t det=h11*h22-h12*h12;
-
-     Double_t dt1,dt2;
-     if (TMath::Abs(det)<1.e-33) {
-        //(quasi)singular Hessian
-        dt1=-gt1; dt2=-gt2;
-     } else {
-        dt1=-(gt1*h22 - gt2*h12)/det; 
-        dt2=-(h11*gt2 - h12*gt1)/det;
-     }
-
-     if ((dt1*gt1+dt2*gt2)>0) {dt1=-dt1; dt2=-dt2;}
-
-     //check delta(phase1) ?
-     //check delta(phase2) ?
-
-     if (TMath::Abs(dt1)/(TMath::Abs(t1)+1.e-3) < 1.e-4)
-     if (TMath::Abs(dt2)/(TMath::Abs(t2)+1.e-3) < 1.e-4) {
-        if ((gt1*gt1+gt2*gt2) > 1.e-4/dy2/dy2) 
-	  AliWarning(" stopped at not a stationary point !");
-        Double_t lmb=h11+h22; lmb=lmb-TMath::Sqrt(lmb*lmb-4*det);
-        if (lmb < 0.) 
-	  AliWarning(" stopped at not a minimum !");
-        break;
-     }
-
-     Double_t dd=dm;
-     for (Int_t div=1 ; ; div*=2) {
-        Evaluate(p1,t1+dt1,r1,g1,gg1);
-        Evaluate(p2,t2+dt2,r2,g2,gg2);
-        dx=r2[0]-r1[0]; dy=r2[1]-r1[1]; dz=r2[2]-r1[2];
-        dd=dx*dx/dx2 + dy*dy/dy2 + dz*dz/dz2;
-	if (dd<dm) break;
-        dt1*=0.5; dt2*=0.5;
-        if (div>512) {
-           AliWarning(" overshoot !"); break;
-        }   
-     }
-     dm=dd;
-
-     t1+=dt1;
-     t2+=dt2;
-
-  }
-
-  if (max<=0) AliWarning(" too many iterations !");
-
-  Double_t cs=TMath::Cos(GetAlpha());
-  Double_t sn=TMath::Sin(GetAlpha());
-  xthis=r1[0]*cs + r1[1]*sn;
-
-  cs=TMath::Cos(p->GetAlpha());
-  sn=TMath::Sin(p->GetAlpha());
-  xp=r2[0]*cs + r2[1]*sn;
-
-  return TMath::Sqrt(dm*TMath::Sqrt(dy2*dz2));
-}
-
-Double_t AliKalmanTrack::
-PropagateToDCA(AliKalmanTrack *p, Double_t d, Double_t x0) {
-  //--------------------------------------------------------------
-  // Propagates this track and the argument track to the position of the
-  // distance of closest approach. 
-  // Returns the (weighed !) distance of closest approach.
-  //--------------------------------------------------------------
-  Double_t xthis,xp;
-  Double_t dca=GetDCA(p,xthis,xp);
-
-  if (!PropagateTo(xthis,d,x0)) {
-    //AliWarning(" propagation failed !");
-    return 1e+33;
-  }  
-
-  if (!p->PropagateTo(xp,d,x0)) {
-    //AliWarning(" propagation failed !";
-    return 1e+33;
-  }  
-
-  return dca;
-}