Update in the interface for HLT and for visualization
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 26 May 2008 13:02:22 +0000 (13:02 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 26 May 2008 13:02:22 +0000 (13:02 +0000)
TRD/AliTRDtrackV1.cxx
TRD/AliTRDtrackV1.h
TRD/AliTRDtrackerV1.cxx
TRD/AliTRDtrackerV1.h

index dacf653..8d5f2a1 100644 (file)
@@ -37,56 +37,99 @@ ClassImp(AliTRDtrackV1)
 ///////////////////////////////////////////////////////////////////////////////
 
 //_______________________________________________________________
-AliTRDtrackV1::AliTRDtrackV1() 
-  :AliTRDtrack()
+AliTRDtrackV1::AliTRDtrackV1() : AliKalmanTrack()
+  ,fPIDquality(0)
+  ,fDE(0.)
+  ,fBackupTrack(0x0)
 {
   //
   // Default constructor
   //
-       
-       for(int ip=0; ip<6; ip++){
-               fTrackletIndex[ip] = -1;
-               fTracklet[ip].Reset();
-       }
+  for(int i =0; i<3; i++) fBudget[i] = 0.;
+  
+  Float_t pid = 1./AliPID::kSPECIES;
+  for(int is =0; is<AliPID::kSPECIES; is++) fPID[is] = pid;
+
+  for(int ip=0; ip<kNplane; ip++){
+    fTrackletIndex[ip] = 0xffff;
+    fTracklet[ip]      = 0x0;
+  }
 }
 
 //_______________________________________________________________
-AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) 
-  :AliTRDtrack(t)
+AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref) : AliKalmanTrack(ref)
+  ,fPIDquality(ref.fPIDquality)
+  ,fDE(ref.fDE)
+  ,fBackupTrack(0x0)
 {
   //
-  // Constructor from AliESDtrack
+  // Copy constructor
   //
 
-       //AliInfo(Form("alpha %f", GetAlpha()));
-       t.GetTRDtracklets(&fTrackletIndex[0]);
-       for(int ip=0; ip<6; ip++) fTracklet[ip].Reset();
+  for(int ip=0; ip<kNplane; ip++){ 
+    fTrackletIndex[ip] = ref.fTrackletIndex[ip];
+    fTracklet[ip]      = ref.fTracklet[ip];
+  }
+
+  for (Int_t i = 0; i < 3;i++) fBudget[i]      = ref.fBudget[i];
+
+       for(Int_t is = 0; is<AliPID::kSPECIES; is++) fPID[is] = ref.fPID[is];
+  
+  AliKalmanTrack::SetNumberOfClusters(ref.GetNumberOfClusters());
 }
 
 //_______________________________________________________________
-AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref) 
-  :AliTRDtrack(ref)
+AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) : AliKalmanTrack()
+  ,fPIDquality(0)
+  ,fDE(0.)
+  ,fBackupTrack(0x0)
 {
   //
-  // Copy constructor
+  // Constructor from AliESDtrack
   //
 
-       for(int ip=0; ip<6; ip++){ 
-               fTrackletIndex[ip] = ref.fTrackletIndex[ip];
-               fTracklet[ip]      = ref.fTracklet[ip];
-       }
+  SetLabel(t.GetLabel());
+  SetChi2(0.0);
+  SetMass(t.GetMass());
+  AliKalmanTrack::SetNumberOfClusters(t.GetTRDncls()); 
+  Int_t ti[kNplane]; t.GetTRDtracklets(&ti[0]);
+  for(int ip=0; ip<kNplane; ip++){ 
+    fTrackletIndex[ip] = ti[ip] < 0 ? 0xffff : ti[ip];
+    fTracklet[ip]      = 0x0;
+  }
+  for(int i =0; i<3; i++) fBudget[i] = 0.;
+  
+  Float_t pid = 1./AliPID::kSPECIES;
+  for(int is =0; is<AliPID::kSPECIES; is++) fPID[is] = pid;
+
+  const AliExternalTrackParam *par = &t;
+  if (t.GetStatus() & AliESDtrack::kTRDbackup) { 
+    par = t.GetOuterParam();
+    if (!par) {
+      AliError("No backup info!"); 
+      par = &t;
+    }
+  }
+  Set(par->GetX() 
+     ,par->GetAlpha()
+     ,par->GetParameter()
+     ,par->GetCovariance());
+
+  if(t.GetStatus() & AliESDtrack::kTIME) {
+    StartTimeIntegral();
+    Double_t times[10]; 
+    t.GetIntegratedTimes(times); 
+    SetIntegratedTimes(times);
+    SetIntegratedLength(t.GetIntegratedLength());
+  }
 }
 
 //_______________________________________________________________
-// AliTRDtrackV1::~AliTRDtrackV1()
-// {
-//     
-// }
-       
-//_______________________________________________________________
 AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Double_t cov[15]
-             , Double_t x, Double_t alpha)
-  :AliTRDtrack()
+             , Double_t x, Double_t alpha) : AliKalmanTrack()
+  ,fPIDquality(0)
+  ,fDE(0.)
+  ,fBackupTrack(0x0)
 {
   //
   // The stand alone tracking constructor
@@ -117,83 +160,79 @@ AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Do
                     , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
   
        Set(x,alpha,pp,cc);
-  Double_t s = GetSnp();
-  Double_t t = GetTgl();
-
-       Int_t ncl = 0, nclPlane; AliTRDcluster *c = 0x0;
-       for(int iplane=0; iplane<6; iplane++){
-               fTrackletIndex[iplane] = -1;
-               fTracklet[iplane] = trklts[iplane];
-               nclPlane = 0;
-               for(int ic = 0; ic<AliTRDseed::knTimebins; ic++){
-                       if(!fTracklet[iplane].IsUsable(ic)) continue;
-                       if(!(c = fTracklet[iplane].GetClusters(ic))) continue;
-                       
-                       fIndex[ncl]      = fTracklet[iplane].GetIndexes(ic);
-                       Double_t q = TMath::Abs(c->GetQ());
-                       fClusters[ncl]   = c;
-                       // temporary !!!
-                       // This is not correctly. Has to be updated in the AliTRDtrackerV1::FollowBackProlonagation()
-                       fdQdl[ncl]       = q * (s*s < 1.) ? TMath::Sqrt((1-s*s)/(1+t*t)) : 1.;
-                       ncl++;
-                       nclPlane++;
-               }
-               //printf("%d N clusters plane %d [%d %d].\n", iplane, nclPlane, fTracklet[iplane].GetN2(), trklts[iplane].GetN());
+       for(int iplane=0; iplane<kNplane; iplane++){
+               fTrackletIndex[iplane] = 0xffff;
+               fTracklet[iplane] = &trklts[iplane];
        }
-       //printf("N clusters in AliTRDtrackV1 %d.\n", ncl);
-       SetNumberOfClusters(/*ncl*/);
+       SetNumberOfClusters();
 }
 
 //_______________________________________________________________
-Bool_t AliTRDtrackV1::CookLabel(Float_t wrong)
+AliTRDtrackV1::~AliTRDtrackV1()
 {
-       // set MC label for this tracklet
+  if(fBackupTrack) {
+    delete fBackupTrack;
+  }
+  fBackupTrack = 0x0;
 
+  if(TestBit(1)){
+    for(Int_t ip=0; ip<kNplane; ip++){
+      if(fTracklet[ip]) delete fTracklet[ip];
+      fTracklet[ip] = 0x0;
+    }
+  }
+}
+       
+//_______________________________________________________________
+Bool_t AliTRDtrackV1::CookLabel(Float_t wrong)
+{
+  // set MC label for this track
+  
   Int_t s[kMAXCLUSTERSPERTRACK][2];
   for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
     s[i][0] = -1;
     s[i][1] =  0;
   }
-
+  
   Bool_t labelAdded;
-       Int_t label;
-       AliTRDcluster *c    = 0x0;
-  for (Int_t ip = 0; ip < AliESDtrack::kTRDnPlanes; ip++) {
-               if(fTrackletIndex[ip] < 0) continue;
-               for (Int_t ic = 0; ic < AliTRDseed::knTimebins; ic++) {
-                       if(!(c = fTracklet[ip].GetClusters(ic))) continue;
-                       for (Int_t k = 0; k < 3; k++) { 
-                               label      = c->GetLabel(k);
-                               labelAdded = kFALSE; 
-                               Int_t j = 0;
-                               if (label >= 0) {
-                                       while ((!labelAdded) && (j < kMAXCLUSTERSPERTRACK)) {
-                                               if ((s[j][0] == label) || 
-                                                               (s[j][1] ==     0)) {
-                                                       s[j][0] = label; 
-                                                       s[j][1]++; 
-                                                       labelAdded = kTRUE;
-                                               }
-                                               j++;
-                                       }
-                               }
-                       }
-               }
-       }
-
+  Int_t label;
+  AliTRDcluster *c    = 0x0;
+  for (Int_t ip = 0; ip < kNplane; ip++) {
+    if(fTrackletIndex[ip] == 0xffff) continue;
+    for (Int_t ic = 0; ic < AliTRDseed::knTimebins; ic++) {
+      if(!(c = fTracklet[ip]->GetClusters(ic))) continue;
+      for (Int_t k = 0; k < 3; k++) { 
+        label      = c->GetLabel(k);
+        labelAdded = kFALSE; 
+        Int_t j = 0;
+        if (label >= 0) {
+          while ((!labelAdded) && (j < kMAXCLUSTERSPERTRACK)) {
+            if ((s[j][0] == label) || 
+                (s[j][1] ==     0)) {
+              s[j][0] = label; 
+              s[j][1]++; 
+              labelAdded = kTRUE;
+            }
+            j++;
+          }
+        }
+      }
+    }
+  }
+  
   Int_t max = 0;
   label = -123456789;
   for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
     if (s[i][1] <= max) continue;
-               max   = s[i][1]; 
-               label = s[i][0];
+    max   = s[i][1]; 
+    label = s[i][0];
   }
-
+  
   if ((1. - Float_t(max)/GetNumberOfClusters()) > wrong) label = -label;
-
+  
   SetLabel(label); 
-       
-       return kTRUE;
+  
+  return kTRUE;
 }
 
 //_______________________________________________________________
@@ -202,66 +241,85 @@ Bool_t AliTRDtrackV1::CookPID()
   //
   // Cook the PID information
   //
-
-// CookdEdx();  // truncated mean ... do we still need it ?
-
-// CookdEdxTimBin(seed->GetID());
-       
-  // Sets the a priori probabilities
+  
+  // Reset the a priori probabilities
+  Double_t pid = 1. / AliPID::kSPECIES;
   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
-    fPID[ispec] = 1.0 / AliPID::kSPECIES;      
+    fPID[ispec] = pid; 
+  }
+  fPIDquality = 0;
+  
+  // steer PID calculation @ tracklet level
+  Double_t *prob = 0x0;
+  for(int ip=0; ip<kNplane; ip++){
+    if(fTrackletIndex[ip] == 0xffff) continue;
+    if(!fTracklet[ip]->IsOK()) continue;
+    if(!(prob = fTracklet[ip]->GetProbability())) return kFALSE;
+    
+    Int_t nspec = 0; // quality check of tracklet dEdx
+    for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
+      if(prob[ispec] < 0.) continue;
+      fPID[ispec] *= prob[ispec];
+      nspec++;
+    }
+    if(!nspec) continue;
+    
+    fPIDquality++;
   }
-       
-       // steer PID calculation @ tracklet level
-       Double_t *prob = 0x0;
-       fPIDquality = 0;
-       for(int itrklt=0; itrklt<AliESDtrack::kTRDnPlanes; itrklt++){
-    //for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) fdEdxPlane[itrklt][iSlice] = -1.;
-
-               if(fTrackletIndex[itrklt]<0) continue;
-               if(!fTracklet[itrklt].IsOK()) continue;
-               if(!(prob = fTracklet[itrklt].GetProbability())) return kFALSE;
-               
-               Int_t nspec = 0; // quality check of tracklet dEdx
-               for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
-                       if(prob[ispec] < 0.) continue;
-                       fPID[ispec] *= prob[ispec];
-                       nspec++;
-               }
-               if(!nspec) continue;
-               
-               fPIDquality++;
-       }
   
   // no tracklet found for PID calculations
   if(!fPIDquality) return kTRUE;
-
-       // slot for PID calculation @ track level
-       
-
+  
+  // slot for PID calculation @ track level
+  
+  
   // normalize probabilities
   Double_t probTotal = 0.0;
   for (Int_t is = 0; is < AliPID::kSPECIES; is++) probTotal += fPID[is];
-
-
+  
+  
   if (probTotal <= 0.0) {
     AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference data.");
     return kFALSE;
   }
-
+  
   for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) fPID[iSpecies] /= probTotal;
-
+  
   return kTRUE;
 }
 
-//_______________________________________________________________
-Float_t AliTRDtrackV1::GetMomentum(Int_t plane) const
+//_____________________________________________________________________________
+Double_t AliTRDtrackV1::GetBz() const 
 {
   //
-  // Get the momentum at a given plane
+  // Returns Bz component of the magnetic field (kG)
   //
 
-       return plane >=0 && plane < 6 && fTrackletIndex[plane] >= 0 ? fTracklet[plane].GetMomentum() : -1.;
+  if (AliTracker::UniformField()) return AliTracker::GetBz();
+
+  Double_t r[3]; 
+  GetXYZ(r);
+  return AliTracker::GetBz(r);
+}
+
+//_______________________________________________________________
+Int_t  AliTRDtrackV1::GetClusterIndex(Int_t id) const
+{
+  Int_t n = 0;
+  for(Int_t ip=0; ip<kNplane; ip++){
+    if(!fTracklet[ip]) continue;
+    if(n+fTracklet[ip]->GetN() < id){ 
+      n+=fTracklet[ip]->GetN();
+      continue;
+    }
+    for(Int_t ic=34; ic>=0; ic--){
+      if(!fTracklet[ip]->GetClusters(ic)) continue;
+      n++;
+      if(n<id) continue;
+      return fTracklet[ip]->GetIndexes(ic);
+    }
+  }
+  return -1;
 }
 
 //_______________________________________________________________
@@ -270,121 +328,276 @@ Double_t AliTRDtrackV1::GetPredictedChi2(const AliTRDseedV1 *trklt) const
   //
   // Get the predicted chi2
   //
-
+  
   Double_t x      = trklt->GetX0();
   Double_t p[2]   = { trklt->GetYat(x)
                     , trklt->GetZat(x) };
   Double_t cov[3];
-       trklt->GetCovAt(x, cov);
-
+  trklt->GetCovAt(x, cov);
+  
   return AliExternalTrackParam::GetPredictedChi2(p, cov);
+}
 
+       
+//_____________________________________________________________________________
+void AliTRDtrackV1::MakeBackupTrack()
+{
+  //
+  // Creates a backup track
+  //
+
+  if(fBackupTrack) {
+    fBackupTrack->~AliTRDtrackV1();
+    new(fBackupTrack) AliTRDtrackV1((AliTRDtrackV1&)(*this));
+  }
+  fBackupTrack = new AliTRDtrackV1((AliTRDtrackV1&)(*this));
 }
 
-//_______________________________________________________________
-Bool_t AliTRDtrackV1::IsOwner() const
+//_____________________________________________________________________________
+Int_t AliTRDtrackV1::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
+{
+  //
+  // Find a prolongation at given x
+  // Return 0 if it does not exist
+  //  
+
+  Double_t bz = GetBz();
+  if (!AliExternalTrackParam::GetYAt(xk,bz,y)) return 0;
+  if (!AliExternalTrackParam::GetZAt(xk,bz,z)) return 0;
+
+  return 1;  
+
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDtrackV1::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
 {
   //
-  // Check whether track owns the tracklets
+  // Propagates this track to a reference plane defined by "xk" [cm] 
+  // correcting for the mean crossed material.
   //
+  // "xx0"  - thickness/rad.length [units of the radiation length] 
+  // "xrho" - thickness*density    [g/cm^2] 
+  // 
 
-       for (Int_t ip = 0; ip < AliESDtrack::kTRDnPlanes; ip++) {
-               if(fTrackletIndex[ip] < 0) continue;
-               if(!fTracklet[ip].IsOwner()) return kFALSE;
-       }
-       return kTRUE;
+  if (xk == GetX()) {
+    return kTRUE;
+  }
+
+  Double_t oldX = GetX();
+  Double_t oldY = GetY();
+  Double_t oldZ = GetZ();
+
+  Double_t bz   = GetBz();
+
+  if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
+    return kFALSE;
+  }
+
+  Double_t x = GetX();
+  Double_t y = GetY();
+  Double_t z = GetZ();
+
+  if (oldX < xk) {
+    xrho = -xrho;
+    if (IsStartedTimeIntegral()) {
+      Double_t l2  = TMath::Sqrt((x-oldX)*(x-oldX) 
+                               + (y-oldY)*(y-oldY) 
+                               + (z-oldZ)*(z-oldZ));
+      Double_t crv = AliExternalTrackParam::GetC(bz);
+      if (TMath::Abs(l2*crv) > 0.0001) {
+        // Make correction for curvature if neccesary
+        l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX) 
+                             + (y-oldY)*(y-oldY));
+        l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
+        l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ));
+      }
+      AddTimeStep(l2);
+    }
+  }
+
+  if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0,xrho,GetMass())) { 
+    return kFALSE;
+  }
+
+  {
+
+    // Energy losses
+    Double_t p2    = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt());
+    Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
+    if ((beta2 < 1.0e-10) || 
+        ((5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0)) {
+      return kFALSE;
+    }
+
+    Double_t dE    = 0.153e-3 / beta2 
+                   * (TMath::Log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
+                   * xrho;
+    fBudget[0] += xrho;
+
+    /*
+    // Suspicious part - think about it ?
+    Double_t kinE =  TMath::Sqrt(p2);
+    if (dE > 0.8*kinE) dE = 0.8 * kinE;  //      
+    if (dE < 0)        dE = 0.0;         // Not valid region for Bethe bloch 
+    */
+    fDE += dE;
+
+    /*
+    // Suspicious ! I.B.
+    Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE));   // Energy loss fluctuation 
+    Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
+    fCcc += sigmac2;
+    fCee += fX*fX * sigmac2;  
+    */
+
+  }
+
+  return kTRUE;
 }
-       
+
 //_____________________________________________________________________________
-void AliTRDtrackV1::MakeBackupTrack()
+Int_t   AliTRDtrackV1::PropagateToR(Double_t r,Double_t step)
 {
   //
-  // Creates a backup track
+  // Propagate track to the radial position
+  // Rotation always connected to the last track position
   //
 
-  if (fBackupTrack) {
-    fBackupTrack->~AliTRDtrack();
-    new(fBackupTrack) AliTRDtrack((AliTRDtrack&)(*this));
+  Double_t xyz0[3];
+  Double_t xyz1[3];
+  Double_t y;
+  Double_t z; 
+
+  Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
+  // Direction +-
+  Double_t dir    = (radius > r) ? -1.0 : 1.0;   
+
+  for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
+
+    GetXYZ(xyz0);      
+    Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
+    Rotate(alpha,kTRUE);
+    GetXYZ(xyz0);      
+    GetProlongation(x,y,z);
+    xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha); 
+    xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
+    xyz1[2] = z;
+    Double_t param[7];
+    AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
+    if (param[1] <= 0) {
+      param[1] = 100000000;
+    }
+    PropagateTo(x,param[1],param[0]*param[4]);
+
+  } 
+
+  GetXYZ(xyz0);        
+  Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
+  Rotate(alpha,kTRUE);
+  GetXYZ(xyz0);        
+  GetProlongation(r,y,z);
+  xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha); 
+  xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
+  xyz1[2] = z;
+  Double_t param[7];
+  AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
+
+  if (param[1] <= 0) {
+    param[1] = 100000000;
   }
-  fBackupTrack = new AliTRDtrack((AliTRDtrack&)(*this));
+  PropagateTo(r,param[1],param[0]*param[4]);
+
+  return 0;
+
 }
 
+//_____________________________________________________________________________
+Bool_t AliTRDtrackV1::Rotate(Double_t alpha, Bool_t absolute)
+{
+  //
+  // Rotates track parameters in R*phi plane
+  // if absolute rotation alpha is in global system
+  // otherwise alpha rotation is relative to the current rotation angle
+  //  
+
+  if (absolute) alpha -= GetAlpha();
+  //else fNRotate++;
+
+  return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
+}
 
 //___________________________________________________________
 void AliTRDtrackV1::SetNumberOfClusters() 
 {
 // Calculate the number of clusters attached to this track
        
-       Int_t ncls = 0;
-       for(int ip=0; ip<6; ip++){
-               if(fTrackletIndex[ip] >= 0) ncls += fTracklet[ip].GetN();
-       }
-       AliKalmanTrack::SetNumberOfClusters(ncls);      
+  Int_t ncls = 0;
+  for(int ip=0; ip<kNplane; ip++){
+    if(fTracklet[ip] && fTrackletIndex[ip] != 0xffff) ncls += fTracklet[ip]->GetN();
+  }
+  AliKalmanTrack::SetNumberOfClusters(ncls);   
 }
 
        
 //_______________________________________________________________
-void AliTRDtrackV1::SetOwner(Bool_t own)
+void AliTRDtrackV1::SetOwner()
 {
   //
   // Toggle ownership of tracklets
   //
 
-       for (Int_t ip = 0; ip < AliESDtrack::kTRDnPlanes; ip++) {
-               if(fTrackletIndex[ip] < 0) continue;
-               //AliInfo(Form("p[%d] index[%d]", ip, fTrackletIndex[ip]));
-               fTracklet[ip].SetOwner(own);
-       }
+  if(TestBit(1)) return;
+  for (Int_t ip = 0; ip < kNplane; ip++) {
+    if(fTrackletIndex[ip] == 0xffff) continue;
+    fTracklet[ip] = new AliTRDseedV1(*fTracklet[ip]);
+    fTracklet[ip]->SetOwner();
+  }
+  SetBit(1);
 }
 
 //_______________________________________________________________
-void AliTRDtrackV1::SetTracklet(AliTRDseedV1 *trklt, Int_t plane, Int_t index)
+void AliTRDtrackV1::SetTracklet(AliTRDseedV1 *trklt, Int_t index)
 {
   //
   // Set the tracklets
   //
-
-       if(plane < 0 || plane >= AliESDtrack::kTRDnPlanes) return;
-       fTracklet[plane]      = (*trklt);
-       fTrackletIndex[plane] = index;
+  Int_t plane = trklt->GetPlane();
+  fTracklet[plane]      = trklt;
+  fTrackletIndex[plane] = index;
 }
 
 //_______________________________________________________________
 Bool_t  AliTRDtrackV1::Update(AliTRDseedV1 *trklt, Double_t chisq)
 {
   //
-  // Update track parameters
+  // Update track and tracklet parameters 
   //
-
+  
   Double_t x      = trklt->GetX0();
   Double_t p[2]   = { trklt->GetYat(x)
                     , trklt->GetZat(x) };
   Double_t cov[3];
-       trklt->GetCovAt(x, cov);
-       
-//     Print();
-//     AliInfo(Form("cov[%f %f %f]", cov[0], cov[1], cov[2]));
-
+  trklt->GetCovAt(x, cov);
+  
+  
   if(!AliExternalTrackParam::Update(p, cov)) return kFALSE;
-       //Print();
-
-       AliTRDcluster *c = 0x0;
-       Int_t ic = 0; while(!(c = trklt->GetClusters(ic))) ic++;
-       AliTracker::FillResiduals(this, p, cov, c->GetVolumeId());
-
-
+  
+  AliTRDcluster *c = 0x0;
+  Int_t ic = 0; while(!(c = trklt->GetClusters(ic))) ic++;
+  AliTracker::FillResiduals(this, p, cov, c->GetVolumeId());
+  
+  
   // Register info to track
-//   Int_t n      = GetNumberOfClusters();
-//   fIndex[n]    = index;
-//   fClusters[n] = c;
-  SetNumberOfClusters(/*GetNumberOfClusters()+trklt->GetN()*/);
+  SetNumberOfClusters();
   SetChi2(GetChi2() + chisq);
-       
-       // update tracklet
-       trklt->SetMomentum(GetP());
-       trklt->SetSnp(GetSnp());
-       trklt->SetTgl(GetTgl());
-       return kTRUE;
+  
+  // update tracklet
+  trklt->SetMomentum(GetP());
+  trklt->SetSnp(GetSnp());
+  trklt->SetTgl(GetTgl());
+  return kTRUE;
 }
 
 //_______________________________________________________________
@@ -395,14 +608,12 @@ void AliTRDtrackV1::UpdateESDtrack(AliESDtrack *track)
   //
 
   track->SetNumberOfTRDslices(kNslice);
-       
+
   for (Int_t ip = 0; ip < kNplane; ip++) {
-      if(fTrackletIndex[ip] < 0) continue;
-      fTracklet[ip].CookdEdx(kNslice);
-      Float_t *dedx = fTracklet[ip].GetdEdx();
-      for (Int_t js = 0; js < kNslice; js++) { 
-          track->SetTRDslice(dedx[js], ip, js);
-      }
+    if(fTrackletIndex[ip] == 0xffff) continue;
+    fTracklet[ip]->CookdEdx(kNslice);
+    Float_t *dedx = fTracklet[ip]->GetdEdx();
+    for (Int_t js = 0; js < kNslice; js++) track->SetTRDslice(dedx[js], ip, js);
   }
 
   // copy PID to ESD
index dec649c..e9248d8 100644 (file)
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
-#ifndef ALITRDTRACK_H
-#include "AliTRDtrack.h"
+#ifndef ALIKALMANTRACK_H
+#include "AliKalmanTrack.h"
 #endif
 
-class AliTRDseedV1;
-class AliESDtrack;
+#ifndef ALIESDTRACK_H
+#include "AliESDtrack.h"
+#endif
+
+#ifndef ALITRDSEEDV1_H
+#include "AliTRDseedV1.h"
+#endif
 
-class AliTRDtrackV1 : public AliTRDtrack
+class AliTRDcluster;
+class AliESDtrack;
+class AliTRDtrackV1 : public AliKalmanTrack
 {
+public:
+  enum { kMAXCLUSTERSPERTRACK = 210 };
+    
+  enum { kNdet      = 540
+        , kNstacks   =  90
+        , kNplane    =   AliESDtrack::kTRDnPlanes
+        , kNcham     =   5
+        , kNsect     =  18
+        , kNslice    =   3
+        , kNMLPslice =   8 };
+  
+  enum AliTRDPIDMethod {
+          kNN = 0
+        , kLQ = 1 };
 
- public:
-       AliTRDtrackV1();
-       AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Double_t cov[15], Double_t x, Double_t alpha);
-       AliTRDtrackV1(const AliESDtrack &ref);
-       AliTRDtrackV1(const AliTRDtrackV1 &ref);
-       AliTRDtrackV1 &operator=(const AliTRDtrackV1 &ref) { *(new(this) AliTRDtrackV1(ref));
-                                                             return *this; }
-
-       Bool_t         CookPID();
-       Bool_t         CookLabel(Float_t wrong);
-       Float_t        GetMomentum(Int_t plane) const;
-       Double_t       GetPredictedChi2(const AliTRDseedV1 *tracklet) const;
-        Double_t       GetPredictedChi2(const AliCluster* /*c*/) const                   { return 0.0; }
-       const AliTRDseedV1* GetTracklet(Int_t plane) const {return plane >=0 && plane <6 ? &fTracklet[plane] : 0x0;}
-       Int_t*         GetTrackletIndexes() {return &fTrackletIndex[0];}
-       
-       Bool_t         IsOwner() const;
-       
+  AliTRDtrackV1();
+  AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Double_t cov[15], Double_t x, Double_t alpha);
+  AliTRDtrackV1(const AliESDtrack &ref);
+  AliTRDtrackV1(const AliTRDtrackV1 &ref);
+  virtual ~AliTRDtrackV1();
+  AliTRDtrackV1 &operator=(const AliTRDtrackV1 &ref) { *(new(this) AliTRDtrackV1(ref)); return *this; }
+  
+  Bool_t         CookPID();
+  Bool_t         CookLabel(Float_t wrong);
+  AliTRDtrackV1* GetBackupTrack() const {return fBackupTrack;}
+  Double_t       GetBudget(Int_t i) const { return fBudget[i];}
+  Double_t       GetC() const { return AliExternalTrackParam::GetC(GetBz());}
+  Int_t          GetClusterIndex(Int_t id) const;
+  Float_t        GetEdep() const {return fDE;}
+  inline Float_t GetMomentum(Int_t plane) const;
+  inline Int_t   GetNCross();
+  Double_t       GetPIDsignal() const {return 0.;}
+  Double_t       GetPID(Int_t is) const { return (is >=0 && is < AliPID::kSPECIES) ? fPID[is] : -1.;}
+  Double_t       GetPredictedChi2(const AliTRDseedV1 *tracklet) const;
+  Double_t       GetPredictedChi2(const AliCluster* /*c*/) const                   { return 0.0; }
+  Int_t          GetProlongation(Double_t xk, Double_t &y, Double_t &z);
+  AliTRDseedV1*  GetTracklet(Int_t plane) const {return plane >=0 && plane <kNplane ? fTracklet[plane] : 0x0;}
+  Int_t          GetTrackletIndex(Int_t plane) const          { return (plane>=0 && plane<kNplane) ? fTrackletIndex[plane] : -1;} 
+  UShort_t*      GetTrackletIndexes() {return &fTrackletIndex[0];}
+  
+  Bool_t         IsOwner() const {return TestBit(1);};
+  Bool_t         IsStopped() const {return TestBit(2);};
+  
   void           MakeBackupTrack();
-       void           SetNumberOfClusters();
-       void           SetOwner(Bool_t own = kTRUE);
-       void           SetTracklet(AliTRDseedV1 *trklt, Int_t plane, Int_t index);
-       Bool_t         Update(AliTRDseedV1 *tracklet, Double_t chi2);
-        Bool_t         Update(const AliTRDcluster *c, Double_t chi2, Int_t index, Double_t h01) 
-                                                           { return AliTRDtrack::Update(c,chi2,index,h01); };
-        Bool_t         Update(const AliCluster *, Double_t, Int_t)                        { return kFALSE; };
-       void           UpdateESDtrack(AliESDtrack *t);
+  Bool_t         PropagateTo(Double_t xr, Double_t x0 = 8.72, Double_t rho = 5.86e-3);
+  Int_t          PropagateToR(Double_t xr, Double_t step);
+  Bool_t         Rotate(Double_t angle, Bool_t absolute = kFALSE);
+  void           SetBudget(Int_t i, Double_t b) {if(i>=0 && i<3) fBudget[i] = b;}
+  void           SetNumberOfClusters();
+  void           SetOwner();
+  void           SetStopped(Bool_t stop) {SetBit(2, stop);}
+  void           SetTracklet(AliTRDseedV1 *trklt,  Int_t index);
+  inline Float_t StatusForTOF();
+  Bool_t         Update(AliTRDseedV1 *tracklet, Double_t chi2);
+  //Bool_t         Update(const AliTRDcluster *c, Double_t chi2, Int_t index, Double_t h01){ return AliTRDtrack::Update(c,chi2,index,h01); };
+  Bool_t         Update(const AliCluster *, Double_t, Int_t)                        { return kFALSE; };
+  void           UpdateESDtrack(AliESDtrack *t);
+
+protected:
+  Double_t       GetBz() const;
+
+private:
+  UChar_t      fPIDquality;           //  No of planes used for PID calculation        
+  UShort_t     fTrackletIndex[kNplane];//  Tracklets index in the tracker list
+  Double32_t   fPID[AliPID::kSPECIES];//  PID probabilities
+  Double32_t   fBudget[3];            //  Integrated material budget
+  Double32_t   fDE;                   //  Integrated delta energy
+  AliTRDseedV1 *fTracklet[kNplane];   //  Tracklets array defining the track
+  AliTRDtrackV1 *fBackupTrack;        // Backup track
 
-       ClassDef(AliTRDtrackV1, 1)         // development TRD track
 
+  ClassDef(AliTRDtrackV1, 2)          // new TRD track
 };
 
+//____________________________________________________
+inline Float_t AliTRDtrackV1::GetMomentum(Int_t plane) const
+{
+  return plane >=0 && plane < kNplane && fTrackletIndex[plane] != 0xff ? fTracklet[plane]->GetMomentum() : -1.;
+}
+
+//____________________________________________________
+inline Int_t AliTRDtrackV1::GetNCross()
+{
+  Int_t ncross = 0;
+  for(Int_t ip=0; ip<kNplane; ip++){
+    if(!fTracklet[ip]) continue;
+    ncross += fTracklet[ip]->IsRowCross();
+  }
+  return ncross;
+}
+
+//____________________________________________________________________________
+inline Float_t AliTRDtrackV1::StatusForTOF()
+{
+  //
+  // Defines the status of the TOF extrapolation
+  //
+
+  if(!fTracklet[5]) return 0.;
+
+  // Definition of res ????
+  Float_t res = /*(0.2 + 0.8 * (fN / (fNExpected + 5.0))) **/ (0.4 + 0.6 * fTracklet[5]->GetN() / 20.0);
+  res *= (0.25 + 0.8 * 40.0 / (40.0 + fBudget[2]));
+  return res;
+}
 
 #endif
 
index fb7e1ca..c3409f6 100644 (file)
@@ -162,7 +162,7 @@ Bool_t AliTRDtrackerV1::GetTrackPoint(Int_t index, AliTrackPoint &p) const
 {
        //AliInfo(Form("Asking for tracklet %d", index));
        
-       if(index<0) return kFALSE;
+       if(index<0 || index == 0xffff) return kFALSE;
        AliTRDseedV1 *tracklet = 0x0; 
        if(!(tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index))) return kFALSE;
        
@@ -281,10 +281,10 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
                if ((status & AliESDtrack::kTRDout) != 0) continue;
        
                // Do the back prolongation
-               Int_t   lbl         = seed->GetLabel();
                new(&track) AliTRDtrackV1(*seed);
                //track->Print();
-               track.SetSeedLabel(lbl);
+               //Int_t   lbl         = seed->GetLabel();
+               //track.SetSeedLabel(lbl);
                seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup); // Make backup
                Float_t p4          = track.GetC();
                if((expectedClr = FollowBackProlongation(track))){
@@ -314,8 +314,8 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
                        Int_t foundClr = track.GetNumberOfClusters();
                        if (foundClr >= foundMin) {
                                //AliInfo(Form("Making backup track ncls [%d]...", foundClr));
-                               track.CookdEdx();
-                               track.CookdEdxTimBin(seed->GetID());
+                               //track.CookdEdx();
+                               //track.CookdEdxTimBin(seed->GetID());
                                track.CookLabel(1. - fgkLabelFraction);
                                if(track.GetBackupTrack()) UseClusters(track.GetBackupTrack());
                                
@@ -355,7 +355,7 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
                }
                
                // Propagation to the TOF (I.Belikov)
-               if (track.GetStop() == kFALSE) {
+               if (track.IsStopped() == kFALSE) {
                        Double_t xtof  = 371.0;
                        Double_t xTOF0 = 370.0;
                
@@ -523,7 +523,7 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
                if(x > t.GetX()+fgkMaxStep) continue;
 
                // append tracklet to track
-               t.SetTracklet(tracklet, iplane, index);
+               t.SetTracklet(tracklet, index);
                
                if (x < (t.GetX()-fgkMaxStep) && !PropagateToX(t, x+fgkMaxStep, fgkMaxStep)) break;
                if (!AdjustSector(&t)) break;
@@ -561,7 +561,7 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
                for(int iplane=0; iplane<6; iplane++){
                        AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
                        if(!tracklet) continue;
-                       t.SetTracklet(tracklet, iplane, index);
+                       t.SetTracklet(tracklet, index);
                }
 
                Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
@@ -608,12 +608,15 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
        Double_t clength = AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick();
        AliTRDtrackingChamber *chamber = 0x0;
        
+  AliTRDseedV1 tracklet, *ptrTracklet = 0x0;
+
        // Loop through the TRD planes
        for (Int_t iplane = 0; iplane < AliTRDgeometry::Nplan(); iplane++) {
                // BUILD TRACKLET IF NOT ALREADY BUILT
                Double_t x = 0., y, z, alpha;
-               AliTRDseedV1 tracklet(*t.GetTracklet(iplane));
-               if(!tracklet.IsOK()){
+               ptrTracklet  = t.GetTracklet(iplane);
+               if(!ptrTracklet){
+                 new(&tracklet) AliTRDseedV1(iplane);
                        alpha = t.GetAlpha();
                        Int_t sector = Int_t(alpha/AliTRDgeometry::GetAlpha() + (alpha>0. ? 0 : AliTRDgeometry::kNsect));
 
@@ -641,7 +644,7 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
                                tracklet.SetPlane(iplane);
                                tracklet.SetX0(x);
                                if(!tracklet.Init(&t)){
-                                       t.SetStop(kTRUE);
+                                       t.SetStopped(kTRUE);
                                        return nClustersExpected;
                                }
                                if(!tracklet.AttachClustersIter(chamber, 1000.)) continue;
@@ -667,8 +670,8 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
                if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) break;
                
                // load tracklet to the tracker and the track
-               Int_t index = SetTracklet(&tracklet);
-               t.SetTracklet(&tracklet, iplane, index);
+               ptrTracklet = SetTracklet(&tracklet);
+               t.SetTracklet(ptrTracklet, fTracklets->GetEntriesFast()-1);
        
        
                // Calculate the mean material budget along the path inside the chamber
@@ -676,7 +679,7 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
                Double_t xyz0[3]; // entry point 
                t.GetXYZ(xyz0);
                alpha = t.GetAlpha();
-               x = tracklet.GetX0();
+               x = ptrTracklet->GetX0();
                if (!t.GetProlongation(x, y, z)) break;
                Double_t xyz1[3]; // exit point
                xyz1[0] =  x * TMath::Cos(alpha) - y * TMath::Sin(alpha); 
@@ -691,19 +694,18 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
                // Propagate and update track
                t.PropagateTo(x, xx0, xrho);
                if (!AdjustSector(&t)) break;
-               Double_t maxChi2 = t.GetPredictedChi2(&tracklet);
-               if (maxChi2<1e+10 && t.Update(&tracklet, maxChi2)){ 
-                       nClustersExpected += tracklet.GetN();
-                       t.SetTracklet(&tracklet, iplane, index);
-                       UpdateTracklet(&tracklet, index);
+               Double_t maxChi2 = t.GetPredictedChi2(ptrTracklet);
+               if (maxChi2<1e+10 && t.Update(ptrTracklet, maxChi2)){ 
+                       nClustersExpected += ptrTracklet->GetN();
+                       //t.SetTracklet(&tracklet, index);
                }
                // Reset material budget if 2 consecutive gold
-               if(iplane>0 && tracklet.GetN() + t.GetTracklet(iplane-1)->GetN() > 20) t.SetBudget(2, 0.);
+               if(iplane>0 && t.GetTracklet(iplane-1) && ptrTracklet->GetN() + t.GetTracklet(iplane-1)->GetN() > 20) t.SetBudget(2, 0.);
 
                // Make backup of the track until is gold
                // TO DO update quality check of the track.
                // consider comparison with fTimeBinsRange
-               Float_t ratio0 = tracklet.GetN() / Float_t(fgNTimeBins);
+               Float_t ratio0 = ptrTracklet->GetN() / Float_t(fgNTimeBins);
                //Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);    
                //printf("tracklet.GetChi2() %f     [< 18.0]\n", tracklet.GetChi2()); 
                //printf("ratio0    %f              [>   0.8]\n", ratio0);
@@ -1062,6 +1064,172 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro
        return chi2track;
 }
 
+
+//_________________________________________________________________________
+Double_t AliTRDtrackerV1::FitRiemanTilt(AliTRDtrackV1 *track, AliTRDseedV1 *tracklets, Bool_t sigError, Int_t np, AliTrackPoint *points)
+{
+       //
+       // Performs a Riemann fit taking tilting pad correction into account
+       // The equation of a Riemann circle, where the y position is substituted by the 
+       // measured y-position taking pad tilting into account, has to be transformed
+       // into a 4-dimensional hyperplane equation
+       // Riemann circle: (x-x0)^2 + (y-y0)^2 -R^2 = 0
+       // Measured y-Position: ymeas = y - tan(phiT)(zc - zt)
+       //          zc: center of the pad row
+       //          zt: z-position of the track
+       // The z-position of the track is assumed to be linear dependent on the x-position
+       // Transformed equation: a + b * u + c * t + d * v  + e * w - 2 * (ymeas + tan(phiT) * zc) * t = 0
+       // Transformation:       u = 2 * x * t
+       //                       v = 2 * tan(phiT) * t
+       //                       w = 2 * tan(phiT) * (x - xref) * t
+       //                       t = 1 / (x^2 + ymeas^2)
+       // Parameters:           a = -1/y0
+       //                       b = x0/y0
+       //                       c = (R^2 -x0^2 - y0^2)/y0
+       //                       d = offset
+       //                       e = dz/dx
+       // If the offset respectively the slope in z-position is impossible, the parameters are fixed using 
+       // results from the simple riemann fit. Afterwards the fit is redone.
+       // The curvature is calculated according to the formula:
+       //                       curv = a/(1 + b^2 + c*a) = 1/R
+       //
+       // Paramters:   - Array of tracklets (connected to the track candidate)
+       //              - Flag selecting the error definition
+       // Output:      - Chi2 values of the track (in Parameter list)
+       //
+       TLinearFitter *fitter = GetTiltedRiemanFitter();
+       fitter->StoreData(kTRUE);
+       fitter->ClearPoints();
+       AliTRDLeastSquare zfitter;
+       AliTRDcluster *cl = 0x0;
+
+  AliTRDseedV1 work[kNPlanes], *tracklet = 0x0;
+  if(!tracklets){
+    for(Int_t ipl = 0; ipl < kNPlanes; ipl++){
+      if(!(tracklet = track->GetTracklet(ipl))) continue;
+      if(!tracklet->IsOK()) continue;
+      new(&work[ipl]) AliTRDseedV1(*tracklet);
+    }
+    tracklets = &work[0];
+  }
+
+       Double_t xref = CalculateReferenceX(tracklets);
+       Double_t x, y, z, t, tilt, dx, w, we;
+       Double_t uvt[4];
+       Int_t nPoints = 0;
+       // Containers for Least-square fitter
+       for(Int_t ipl = 0; ipl < kNPlanes; ipl++){
+               if(!tracklets[ipl].IsOK()) continue;
+               for(Int_t itb = 0; itb < fgNTimeBins; itb++){
+                       if(!(cl = tracklets[ipl].GetClusters(itb))) continue;
+                       if (!tracklets[ipl].IsUsable(itb)) continue;
+                       x = cl->GetX();
+                       y = cl->GetY();
+                       z = cl->GetZ();
+                       tilt = tracklets[ipl].GetTilt();
+                       dx = x - xref;
+                       // Transformation
+                       t = 1./(x*x + y*y);
+                       uvt[0] = 2. * x * t;
+                       uvt[1] = t;
+                       uvt[2] = 2. * tilt * t;
+                       uvt[3] = 2. * tilt * dx * t;
+                       w = 2. * (y + tilt*z) * t;
+                       // error definition changes for the different calls
+                       we = 2. * t;
+                       we *= sigError ? tracklets[ipl].GetSigmaY() : 0.2;
+                       fitter->AddPoint(uvt, w, we);
+                       zfitter.AddPoint(&x, z, static_cast<Double_t>(TMath::Sqrt(cl->GetSigmaZ2())));
+                       nPoints++;
+               }
+       }
+       fitter->Eval();
+       Double_t z0    = fitter->GetParameter(3);
+       Double_t dzdx  = fitter->GetParameter(4);
+
+
+  // Linear fitter  - not possible to make boundaries
+  // Do not accept non possible z and dzdx combinations
+  Bool_t accept = kTRUE;
+  Double_t zref = 0.0;
+  for (Int_t iLayer = 0; iLayer < kNPlanes; iLayer++) {
+    if(!tracklets[iLayer].IsOK()) continue;
+    zref = z0 + dzdx * (tracklets[iLayer].GetX0() - xref);
+    if (TMath::Abs(tracklets[iLayer].GetZProb() - zref) > tracklets[iLayer].GetPadLength() * 0.5 + 1.0) 
+      accept = kFALSE;
+  }
+  if (!accept) {
+       zfitter.Eval();
+    Double_t dzmf      = zfitter.GetFunctionParameter(1);
+    Double_t zmf       = zfitter.GetFunctionValue(&xref);
+    fitter->FixParameter(3, zmf);
+    fitter->FixParameter(4, dzmf);
+    fitter->Eval();
+    fitter->ReleaseParameter(3);
+    fitter->ReleaseParameter(4);
+    z0   = fitter->GetParameter(3); // = zmf ?
+    dzdx = fitter->GetParameter(4); // = dzmf ?
+  }
+
+  // Calculate Curvature
+  Double_t a    =  fitter->GetParameter(0);
+  Double_t b    =  fitter->GetParameter(1);
+  Double_t c    =  fitter->GetParameter(2);
+  Double_t y0   = 1. / a;
+  Double_t x0   = -b * y0;
+  Double_t R    = TMath::Sqrt(y0*y0 + x0*x0 - c*y0);
+  Double_t C    =  1.0 + b*b - c*a;
+  if (C > 0.0) C  =  a / TMath::Sqrt(C);
+
+  // Calculate chi2 of the fit 
+  Double_t chi2 = fitter->GetChisquare()/Double_t(nPoints);
+
+  // Update the tracklets
+  if(!track){
+    for(Int_t ip = 0; ip < kNPlanes; ip++) {
+      x = tracklets[ip].GetX0();
+      Double_t tmp = TMath::Sqrt(R*R-(x-x0)*(x-x0));  
+
+      // y:     R^2 = (x - x0)^2 + (y - y0)^2
+      //     =>   y = y0 +/- Sqrt(R^2 - (x - x0)^2)
+      tracklets[ip].SetYref(0, y0 - (y0>0.?1.:-1)*tmp);
+      //     => dy/dx = (x - x0)/Sqrt(R^2 - (x - x0)^2) 
+      tracklets[ip].SetYref(1, (x - x0) / tmp);
+      tracklets[ip].SetZref(0, z0 + dzdx * (x - xref));
+      tracklets[ip].SetZref(1, dzdx);
+      tracklets[ip].SetC(C);
+      tracklets[ip].SetChi2(chi2);
+    }
+  }
+
+  //update track points array
+  if(np && points){
+    Float_t xyz[3];
+    for(int ip=0; ip<np; ip++){
+      points[ip].GetXYZ(xyz);
+      xyz[1] = y0 - (y0>0.?1.:-1)*TMath::Sqrt(R*R-(xyz[0]-x0)*(xyz[0]-x0));
+      xyz[2] = z0 + dzdx * (xyz[0] - xref);
+      points[ip].SetXYZ(xyz);
+    }
+  }
+  
+  if(AliTRDReconstructor::StreamLevel() >=5){
+    TTreeSRedirector &cstreamer = *fgDebugStreamer;
+    Int_t eventNumber                  = AliTRDtrackerDebug::GetEventNumber();
+    Int_t candidateNumber      = AliTRDtrackerDebug::GetCandidateNumber();
+    Double_t chi2z = CalculateChi2Z(tracklets, z0, dzdx, xref);
+    cstreamer << "FitRiemanTilt"
+        << "EventNumber="                      << eventNumber
+        << "CandidateNumber="  << candidateNumber
+        << "xref="                                             << xref
+        << "Chi2Z="                                            << chi2z
+        << "\n";
+  }
+  return chi2;
+}
+
+
+
 //_________________________________________________________________________
 Float_t AliTRDtrackerV1::CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset, Double_t slope, Double_t xref)
 {
@@ -1302,13 +1470,13 @@ AliTRDseedV1* AliTRDtrackerV1::GetTracklet(AliTRDtrackV1 *track, Int_t p, Int_t
        // Detailed description
        //
        idx = track->GetTrackletIndex(p);
-       AliTRDseedV1 *tracklet = idx<0 ? 0x0 : (AliTRDseedV1*)fTracklets->UncheckedAt(idx);
+       AliTRDseedV1 *tracklet = (idx==0xffff) ? 0x0 : (AliTRDseedV1*)fTracklets->UncheckedAt(idx);
 
        return tracklet;
 }
 
 //____________________________________________________________________
-Int_t AliTRDtrackerV1::SetTracklet(AliTRDseedV1 *tracklet)
+AliTRDseedV1* AliTRDtrackerV1::SetTracklet(AliTRDseedV1 *tracklet)
 {
        // Add this tracklet to the list of tracklets stored in the tracker
        //
@@ -1327,28 +1495,9 @@ Int_t AliTRDtrackerV1::SetTracklet(AliTRDseedV1 *tracklet)
                fTracklets->SetOwner(kTRUE);
        }
        Int_t nentries = fTracklets->GetEntriesFast();
-       new ((*fTracklets)[nentries]) AliTRDseedV1(*tracklet);
-       return nentries;
+       return new ((*fTracklets)[nentries]) AliTRDseedV1(*tracklet);
 }
 
-//____________________________________________________________________
-Bool_t AliTRDtrackerV1::UpdateTracklet(AliTRDseedV1 *tracklet, Int_t index){
-       //
-       // Update Tracklet in the tracklet container
-       // 
-       // Parameters:
-       //  - the tracklet information
-       //  - the index of the tracklet which has to be updated
-       //
-       // Output:
-       //  - True if successfull
-       //  - false if the container doesn't exist or the index is out of range
-       //
-       if(!fTracklets || index >= fTracklets->GetEntriesFast()) return kFALSE;
-       
-       new((*fTracklets)[index]) AliTRDseedV1(*tracklet);
-       return kTRUE;
-}
 
 //____________________________________________________________________
 Int_t AliTRDtrackerV1::Clusters2TracksSM(Int_t sector, AliESDEvent *esd)
@@ -2184,8 +2333,8 @@ AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params)
                delete track;
                track = 0x0;
        } else {
-               track->CookdEdx();
-               track->CookdEdxTimBin(-1);
+               //track->CookdEdx();
+               //track->CookdEdxTimBin(-1);
                track->CookLabel(.9);
                // computes PID for track
                track->CookPID();
@@ -2639,6 +2788,20 @@ AliCluster* AliTRDtrackerV1::GetCluster(Int_t idx) const
 }
 
 //____________________________________________________________________
+AliTRDseedV1* AliTRDtrackerV1::GetTracklet(Int_t idx) const
+{
+       Int_t ntrklt = fTracklets->GetEntriesFast();
+       return idx >= 0 || idx < ntrklt ? (AliTRDseedV1*)fTracklets->UncheckedAt(idx) : 0x0;
+}
+
+//____________________________________________________________________
+AliKalmanTrack* AliTRDtrackerV1::GetTrack(Int_t idx) const
+{
+       Int_t ntrk = fTracks->GetEntriesFast();
+       return idx >= 0 || idx < ntrk ? (AliKalmanTrack*)fTracks->UncheckedAt(idx) : 0x0;
+}
+
+//____________________________________________________________________
 Float_t AliTRDtrackerV1::CalculateReferenceX(AliTRDseedV1 *tracklets){
        //
        // Calculates the reference x-position for the tilted Rieman fit defined as middle
index 95120f9..d872fd6 100644 (file)
@@ -49,80 +49,87 @@ class AliTrackPoint;
 class AliTRDtrackerV1 : public AliTracker
 {
 public:
-       enum{
-               kMaxLayersPerSector   = 1000
-               , kMaxTimeBinIndex    = 216
-               , kTrackingSectors    = 18
-               , kNTimeBins          = 35
-               , kNPlanes            = 6
-               , kNSeedPlanes        = 4
-               , kMaxTracksStack     = 100
-               , kNConfigs           = 15
-       };
-       AliTRDtrackerV1();
-       virtual ~AliTRDtrackerV1();
+  enum{
+    kMaxLayersPerSector   = 1000
+    , kMaxTimeBinIndex    = 216
+    , kTrackingSectors    = 18
+    , kNTimeBins          = 35
+    , kNPlanes            = 6
+    , kNSeedPlanes        = 4
+    , kMaxTracksStack     = 100
+    , kNConfigs           = 15
+  };
+  AliTRDtrackerV1();
+  virtual ~AliTRDtrackerV1();
   
   //temporary
   AliTRDtrackingSector* GetTrackingSector(Int_t sec) {return &fTrSec[sec];}
   
-       Int_t          Clusters2Tracks(AliESDEvent *esd);
+  Int_t           Clusters2Tracks(AliESDEvent *esd);
   static TTreeSRedirector* DebugStreamer() {return fgDebugStreamer;}
-  static void SetNTimeBins(Int_t nTimeBins){fgNTimeBins = nTimeBins; }
-  AliCluster*    GetCluster(Int_t index) const;
-       static void    GetExtrapolationConfig(Int_t iconfig, Int_t planes[2]);
-       static const Int_t   GetNTimeBins() {return fgNTimeBins;}
-       static void    GetSeedingConfig(Int_t iconfig, Int_t planes[4]);
-       static TLinearFitter *GetTiltedRiemanFitter();
-       static TLinearFitter *GetTiltedRiemanFitterConstraint();
-       static AliRieman *GetRiemanFitter();
-       static void    FitRieman(AliTRDcluster **clusters, Double_t chi2[2]);
-       static Float_t FitRieman(AliTRDseedV1 *tracklets, Double_t *chi2, Int_t *planes = 0x0);
-       static Float_t FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Double_t zVertex);
-       static Float_t FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigError);
-       
-       Int_t          FollowBackProlongation(AliTRDtrackV1 &t);
-       Int_t          FollowProlongation(AliTRDtrackV1 &t);
-  Int_t          LoadClusters(TTree *cTree);
-       Int_t          PropagateBack(AliESDEvent *event);
-  Int_t          ReadClusters(TClonesArray* &array, TTree *in) const;
-       Int_t          RefitInward(AliESDEvent *event);
-       void           UnloadClusters();
+  AliCluster*     GetCluster(Int_t index) const;
+  AliTRDseedV1*   GetTracklet(Int_t index) const;
+  AliKalmanTrack* GetTrack(Int_t index) const;
+  TClonesArray*   GetListOfClusters() {return fClusters;}
+  TClonesArray*   GetListOfTracklets() {return fTracklets;}
+  TClonesArray*   GetListOfTracks() {return fTracks;}
+  static const Int_t     GetNTimeBins() {return fgNTimeBins;}
+  static void     GetExtrapolationConfig(Int_t iconfig, Int_t planes[2]);
+  static void     GetSeedingConfig(Int_t iconfig, Int_t planes[4]);
+  static TLinearFitter*  GetTiltedRiemanFitter();
+  static TLinearFitter*  GetTiltedRiemanFitterConstraint();
+  static AliRieman*      GetRiemanFitter();
+  static void     FitRieman(AliTRDcluster **clusters, Double_t chi2[2]);
+  static Float_t  FitRieman(AliTRDseedV1 *tracklets, Double_t *chi2, Int_t *planes = 0x0);
+  static Float_t  FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Double_t zVertex);
+  static Float_t  FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigError);
+  
+       static Double_t FitRiemanTilt(AliTRDtrackV1 *trk, AliTRDseedV1 *tracklets = 0x0, Bool_t err=0, Int_t np = 0, AliTrackPoint *points = 0x0);
+       //static Double_t FitKalman(AliTRDtrackV1 *trk, AliTRDseedV1 *tracklets = 0x0, Bool_t up=0, Int_t np = 0, AliTrackPoint *points = 0x0);
+
+  Int_t           FollowBackProlongation(AliTRDtrackV1 &t);
+  Int_t           FollowProlongation(AliTRDtrackV1 &t);
+  Int_t           LoadClusters(TTree *cTree);
+  Int_t           PropagateBack(AliESDEvent *event);
+  Int_t           ReadClusters(TClonesArray* &array, TTree *in) const;
+  Int_t           RefitInward(AliESDEvent *event);
+  static void     SetNTimeBins(Int_t nTimeBins){fgNTimeBins = nTimeBins; }
+  void            UnloadClusters();
   
-  static Int_t   Freq(Int_t n, const Int_t *inlist, Int_t *outlist, Bool_t down); // to be removed 
+  static Int_t    Freq(Int_t n, const Int_t *inlist, Int_t *outlist, Bool_t down); // to be removed 
 
-       class AliTRDLeastSquare{
-       public:
-                                       AliTRDLeastSquare();
-                                       ~AliTRDLeastSquare(){};
-               
-               void            AddPoint(Double_t *x, Double_t y, Double_t sigmaY);
-               void            RemovePoint(Double_t *x, Double_t y, Double_t sigmaY);
-               void            Eval();
-               
-               Double_t        GetFunctionParameter(Int_t ParNumber) const {return fParams[ParNumber];}
-               Double_t        GetFunctionValue(Double_t *xpos) const;
-               void            GetCovarianceMatrix(Double_t *storage) const;
-       private:
-                                       AliTRDLeastSquare(const AliTRDLeastSquare &);
-               AliTRDLeastSquare& operator=(const AliTRDLeastSquare &);
-               Double_t        fParams[2];                                             // Fitparameter 
-               Double_t        fCovarianceMatrix[3];                   // Covariance Matrix
-               Double_t        fSums[6];                                               // Sums
-       };
+  class AliTRDLeastSquare{
+  public:
+    AliTRDLeastSquare();
+    ~AliTRDLeastSquare(){};
+    
+    void               AddPoint(Double_t *x, Double_t y, Double_t sigmaY);
+    void               RemovePoint(Double_t *x, Double_t y, Double_t sigmaY);
+    void               Eval();
+    
+    Double_t   GetFunctionParameter(Int_t ParNumber) const {return fParams[ParNumber];}
+    Double_t   GetFunctionValue(Double_t *xpos) const;
+    void               GetCovarianceMatrix(Double_t *storage) const;
+  private:
+          AliTRDLeastSquare(const AliTRDLeastSquare &);
+    AliTRDLeastSquare& operator=(const AliTRDLeastSquare &);
+    Double_t   fParams[2];                                             // Fitparameter 
+    Double_t   fCovarianceMatrix[3];                   // Covariance Matrix
+    Double_t   fSums[6];                                               // Sums
+  };
 
 protected:
-  Bool_t         AdjustSector(AliTRDtrackV1 *track); 
-       Double_t       BuildSeedingConfigs(AliTRDtrackingChamber **stack, Int_t *configs);
-       static Float_t CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset, Double_t slope, Double_t xref);
-       Int_t          Clusters2TracksSM(Int_t sector, AliESDEvent *esd);
-       Int_t          Clusters2TracksStack(AliTRDtrackingChamber **stack, TClonesArray *esdTrackList);
-       AliTRDseedV1*  GetTracklet(AliTRDtrackV1 *trk, Int_t plane, Int_t &idx);
-       Bool_t         GetTrackPoint(Int_t index, AliTrackPoint &p) const;      
-       Int_t          MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *sseed, Int_t *ipar);
-       AliTRDtrackV1* MakeTrack(AliTRDseedV1 *seeds, Double_t *params);
-  Int_t          PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t maxStep);
-       Int_t          SetTracklet(AliTRDseedV1 *tracklet);
-       Bool_t         UpdateTracklet(AliTRDseedV1 *tracklet, Int_t index);
+  static Bool_t  AdjustSector(AliTRDtrackV1 *track); 
+  Double_t       BuildSeedingConfigs(AliTRDtrackingChamber **stack, Int_t *configs);
+  static Float_t CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset, Double_t slope, Double_t xref);
+  Int_t          Clusters2TracksSM(Int_t sector, AliESDEvent *esd);
+  Int_t          Clusters2TracksStack(AliTRDtrackingChamber **stack, TClonesArray *esdTrackList);
+  AliTRDseedV1*  GetTracklet(AliTRDtrackV1 *trk, Int_t plane, Int_t &idx);
+  Bool_t         GetTrackPoint(Int_t index, AliTrackPoint &p) const;   
+  Int_t          MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *sseed, Int_t *ipar);
+  AliTRDtrackV1* MakeTrack(AliTRDseedV1 *seeds, Double_t *params);
+  static Int_t   PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t maxStep);
+  AliTRDseedV1*  SetTracklet(AliTRDseedV1 *tracklet);
 
 private:
        AliTRDtrackerV1(const AliTRDtrackerV1 &tracker);