]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
New tracker by Markus and Alexandru
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 9 Jan 2008 13:57:03 +0000 (13:57 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 9 Jan 2008 13:57:03 +0000 (13:57 +0000)
16 files changed:
TRD/AliTRDgeometry.cxx
TRD/AliTRDrecoParam.cxx
TRD/AliTRDrecoParam.h
TRD/AliTRDseed.cxx
TRD/AliTRDseed.h
TRD/AliTRDseedV1.cxx
TRD/AliTRDseedV1.h
TRD/AliTRDtrack.cxx
TRD/AliTRDtrack.h
TRD/AliTRDtrackV1.cxx
TRD/AliTRDtrackV1.h
TRD/AliTRDtracker.h
TRD/AliTRDtrackerV1.cxx
TRD/AliTRDtrackerV1.h
TRD/TRDrecLinkDef.h
TRD/libTRDrec.pkg

index c47c0b00bfe73ae5affbd1e9c929d22de6cdf151..889b5af2d3c4b4447e783d74cf5272df3df7d2a5 100644 (file)
@@ -2420,33 +2420,28 @@ Int_t AliTRDgeometry::GetChamber(Double_t z, Int_t plane)
   //
   // Reconstruct the chamber number from the z position and plane number
   //
-  // The return function has to be protected for positiveness !!
-  //
+  // The return function has to be protected for pozitiveness !!
 
-  if ((plane <         0) || 
-      (plane >= fgkNplan)) {
-    return -1;
-  }
+       if(plane<0 || plane>=fgkNplan) return -1;
        
-  Int_t    ichmb = fgkNcham;
-  Double_t zmin;
-  Double_t zmax;
-
-  do {
-    ichmb--;
-    if (ichmb < 0) break;
-    AliTRDpadPlane *pp = GetPadPlane(plane,ichmb);
-    zmax = pp->GetRow0();
-    Int_t nrows = pp->GetNrows();
-    zmin = zmax - 2*pp->GetLengthOPad() 
-                - (nrows-2)*pp->GetLengthIPad() 
-                - (nrows-1)*pp->GetRowSpacing(); 
-  } while((z < zmin) || (z > zmax));
+       Int_t ichmb = fgkNcham;
+       Double_t zmin, zmax;
+  //printf("Looking for z[%7.3f] in plane %d\n", z, plane);
+  do{
+       ichmb--;
+       if(ichmb<0) break;
+       AliTRDpadPlane *pp = GetPadPlane(plane, ichmb);
+    zmax  = pp->GetRow0();
+    // why don't we have the function AliTRDpadPlane::GetLength() ???
+               Int_t nrows = pp->GetNrows();
+               zmin = zmax - 2*pp->GetLengthOPad() - (nrows-2)*pp->GetLengthIPad() - (nrows-1)*pp->GetRowSpacing(); // pp->GetLength();
+       //printf("%d %7.3f %7.3f\n", ichmb, zmin, zmax);
+  } while(z < zmin || z > zmax);
   
   return ichmb;
-
 }
 
+
 //_____________________________________________________________________________
 Int_t AliTRDgeometry::GetSector(Int_t d) const
 {
index fff434364b167ed0eb18afcc9b3da2e8383ea3d9..e771511afee89afd09782287892b7625582671f8 100644 (file)
@@ -32,6 +32,8 @@ ClassImp(AliTRDrecoParam)
 //______________________________________________________________
 AliTRDrecoParam::AliTRDrecoParam()
   :AliDetectorRecoParam()
+  ,fkClusterSharing(0)
+  ,fkPIDMethod(1) // LQ PID
   ,fkMaxTheta(1.0)
   ,fkMaxPhi(2.0)
   ,fkRoad0y(6.0)
index 96b6beeb1af5efe8ee5bdcb71387f5e31529297d..323b50bdad74b08d39d33bcb44e958102203d294 100644 (file)
 
 class AliTRDrecoParam : public AliDetectorRecoParam
 {
-
   public:
-
+       enum{
+         kNNslices = 8,
+         kLQslices = 3
+       };
+       
        AliTRDrecoParam();
        ~AliTRDrecoParam() { }
 
        Double_t GetChi2Y() const                 { return fkChi2Y;    }
        Double_t GetChi2Z() const                 { return fkChi2Z;    }
+       Bool_t   GetClusterSharing() const        { return fkClusterSharing;}
        Double_t GetFindableClusters() const      { return fkFindable; }
        Double_t GetMaxTheta() const              { return fkMaxTheta; }
        Double_t GetMaxPhi() const                { return fkMaxPhi;   }
-
+       Int_t    GetNdEdxSlices() const           { return fkPIDMethod ? kNNslices : kLQslices;}
+       Int_t    GetPIDMethod() const             { return fkPIDMethod;}
        Double_t GetRoad0y() const                { return fkRoad0y;   }
        Double_t GetRoad0z() const                { return fkRoad0z;   }
 
@@ -44,9 +49,13 @@ class AliTRDrecoParam : public AliDetectorRecoParam
        
        static   AliTRDrecoParam *GetLowFluxParam();
         static   AliTRDrecoParam *GetHighFluxParam();
-
+       void     SetClusterSharing(Bool_t share = kTRUE) {fkClusterSharing = share;}
+       void     SetPIDMethod(Int_t pid=1)        { fkPIDMethod = pid ? 1 : 0;}
+       
  private:
 
+       Bool_t   fkClusterSharing;        // Toggle cluster sharing
+       Int_t    fkPIDMethod;             // PID method selector 0(LQ) 1(NN)
        Double_t fkMaxTheta;              // Maximum theta
        Double_t fkMaxPhi;                // Maximum phi
 
@@ -58,7 +67,7 @@ class AliTRDrecoParam : public AliDetectorRecoParam
 
        Double_t fkRoad2y;                // Road in y for extrapolated cluster
        Double_t fkRoad2z;                // Road in z for extrapolated cluster
-
+       
        Double_t fkPlaneQualityThreshold; // Quality threshold
        Double_t fkFindable;              // Ratio of clusters from a track in one chamber which are at minimum supposed to be found.
        Double_t fkChi2Z;                 // Max chi2 on the z direction for seeding clusters fit
index 7d1129f1814592d1c3932893531b5096946926a4..c5a368ad587ba3a17062a77af1c109513b8299e9 100644 (file)
 
 ClassImp(AliTRDseed)
 
+Int_t AliTRDseed::fgTimeBins = 0;
+
 //_____________________________________________________________________________
 AliTRDseed::AliTRDseed() 
   :TObject()
+  ,fTimeBinsRange(0)
+  ,fTimeBin0(0)
   ,fTilt(0)
   ,fPadLength(0)
   ,fX0(0)
@@ -82,6 +86,8 @@ AliTRDseed::AliTRDseed()
 //_____________________________________________________________________________
 AliTRDseed::AliTRDseed(const AliTRDseed &s)
   :TObject(s)
+  ,fTimeBinsRange(s.fTimeBinsRange)
+  ,fTimeBin0(s.fTimeBin0)
   ,fTilt(s.fTilt)
   ,fPadLength(s.fPadLength)
   ,fX0(s.fX0)
@@ -126,14 +132,15 @@ AliTRDseed::AliTRDseed(const AliTRDseed &s)
 }
 
 //_____________________________________________________________________________
-void AliTRDseed::Copy(TObject &o) const 
+void AliTRDseed::Copy(TObject &o) const
 {
-  //
-  // Copy function
-  //
+       //printf("AliTRDseed::Copy()\n");
 
-  AliTRDseed &seed = (AliTRDseed &)o;
-  seed.fTilt = fTilt;
+       AliTRDseed &seed = (AliTRDseed &)o;
+  
+       seed.fTimeBinsRange = fTimeBinsRange;
+       seed.fTimeBin0 = fTimeBin0;
+       seed.fTilt = fTilt;
   seed.fPadLength = fPadLength;
   seed.fX0 = fX0;
   seed.fSigmaY = fSigmaY;
@@ -150,8 +157,7 @@ void AliTRDseed::Copy(TObject &o) const
   seed.fCC = fCC;
   seed.fChi2 = fChi2;
   seed.fChi2Z = fChi2Z;
-
-  for (Int_t i = 0; i < knTimebins; i++) {
+       for (Int_t i = 0; i < knTimebins; i++) {
     seed.fX[i]        = fX[i];
     seed.fY[i]        = fY[i]; 
     seed.fZ[i]        = fZ[i]; 
@@ -218,14 +224,11 @@ void AliTRDseed::CookLabels()
   // Cook 2 labels for seed
   //
 
-  AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-  Int_t nTimeBins = cal->GetNumberOfTimeBins();
-
   Int_t labels[200];
   Int_t out[200];
   Int_t nlab = 0;
 
-  for (Int_t i = 0; i < nTimeBins+1; i++) {
+  for (Int_t i = 0; i < fgTimeBins+1; i++) {
     if (!fClusters[i]) continue;
     for (Int_t ilab = 0; ilab < 3; ilab++) {
       if (fClusters[i]->GetLabel(ilab) >= 0) {
@@ -251,10 +254,7 @@ void AliTRDseed::UseClusters()
   // Use clusters
   //
 
-  AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-  Int_t nTimeBins = cal->GetNumberOfTimeBins();
-
-  for (Int_t i = 0; i < nTimeBins+1; i++) {
+  for (Int_t i = 0; i < fgTimeBins+1; i++) {
     if (!fClusters[i]) continue;
     if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
   }
@@ -268,12 +268,15 @@ void AliTRDseed::Update()
   // Update the seed.
   //
 
+
+       // linear fit on the y direction
+       // dy|x = (yc|x - dz|x*tg(tilt)) - (y0 + dy/dx|x * x )
+       // dz|x = zc|x - (z0 + dz/dx|x) 
+
   const Float_t kRatio  = 0.8;
   const Int_t   kClmin  = 5;
   const Float_t kmaxtan = 2;
 
-  AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-  Int_t nTimeBins = cal->GetNumberOfTimeBins();
 
   if (TMath::Abs(fYref[1]) > kmaxtan){
                //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
@@ -303,7 +306,7 @@ void AliTRDseed::Update()
   
   fN  = 0; 
   fN2 = 0;
-  for (Int_t i = 0; i < nTimeBins; i++) {
+  for (Int_t i = 0; i < fTimeBinsRange; i++) {
     yres[i] = 10000.0;
     if (!fClusters[i]) continue;
     yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i];   // Residual y
@@ -340,22 +343,22 @@ void AliTRDseed::Update()
     // with maximal numebr of accepted clusters
     //
     fNChange = 1;
-    for (Int_t i = 0; i < nTimeBins; i++) {
+    for (Int_t i = 0; i < fTimeBinsRange; i++) {
       cumul[i][0] = counts[0];
       cumul[i][1] = counts[1];
       if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
       if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
     }
     Int_t  maxcount = 0;
-    for (Int_t i = 0; i < nTimeBins; i++) {
-      Int_t after  = cumul[nTimeBins][0] - cumul[i][0];
+    for (Int_t i = 0; i < fTimeBinsRange; i++) {
+      Int_t after  = cumul[fTimeBinsRange][0] - cumul[i][0];
       Int_t before = cumul[i][1];
       if (after + before > maxcount) { 
        maxcount  = after + before; 
        breaktime = i;
        mbefore   = kFALSE;
       }
-      after  = cumul[nTimeBins-1][1] - cumul[i][1];
+      after  = cumul[fTimeBinsRange-1][1] - cumul[i][1];
       before = cumul[i][0];
       if (after + before > maxcount) { 
        maxcount  = after + before; 
@@ -368,18 +371,18 @@ void AliTRDseed::Update()
 
   }
 
-  for (Int_t i = 0; i < nTimeBins+1; i++) {
+  for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
     if (i >  breaktime) allowedz[i] =   mbefore  ? zouts[2] : zouts[0];
     if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
   }  
 
-  if (((allowedz[0] > allowedz[nTimeBins]) && (fZref[1] < 0)) || 
-      ((allowedz[0] < allowedz[nTimeBins]) && (fZref[1] > 0))) {
+  if (((allowedz[0] > allowedz[fTimeBinsRange]) && (fZref[1] < 0)) ||
+      ((allowedz[0] < allowedz[fTimeBinsRange]) && (fZref[1] > 0))) {
     //
     // Tracklet z-direction not in correspondance with track z direction 
     //
     fNChange = 0;
-    for (Int_t i = 0; i < nTimeBins+1; i++) {
+    for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
       allowedz[i] = zouts[0];  // Only longest taken
     } 
   }
@@ -388,7 +391,7 @@ void AliTRDseed::Update()
     //
     // Cross pad -row tracklet  - take the step change into account
     //
-    for (Int_t i = 0; i < nTimeBins+1; i++) {
+    for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
       if (!fClusters[i]) continue; 
       if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
       yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i];   // Residual y
@@ -402,7 +405,7 @@ void AliTRDseed::Update()
   Double_t yres2[knTimebins];
   Double_t mean;
   Double_t sigma;
-  for (Int_t i = 0; i < nTimeBins+1; i++) {
+  for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
     if (!fClusters[i]) continue;
     if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
     yres2[fN2] = yres[i];
@@ -432,12 +435,12 @@ void AliTRDseed::Update()
   fMeanz = 0;
   fMPads = 0;
 
-  for (Int_t i = 0; i < nTimeBins+1; i++) {
+  for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
 
     fUsable[i] = kFALSE;
     if (!fClusters[i]) continue;
-    if (TMath::Abs(fZ[i] - allowedz[i]) > 2)  continue;
-    if (TMath::Abs(yres[i] - mean) > 4.0 * sigma) continue;
+    if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = 0x0; continue;}
+    if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = 0x0;  continue;}
     fUsable[i] = kTRUE;
     fN2++;
     fMPads += fClusters[i]->GetNPads();
@@ -445,7 +448,10 @@ void AliTRDseed::Update()
     if (fClusters[i]->GetNPads() > 4) weight = 0.5;
     if (fClusters[i]->GetNPads() > 5) weight = 0.2;
    
+       
     Double_t x = fX[i];
+    //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
+    
     sumw   += weight; 
     sumwx  += x * weight; 
     sumwx2 += x*x * weight;
@@ -474,12 +480,14 @@ void AliTRDseed::Update()
   fYfitR[1]    = (sumw   * sumwxy - sumwx * sumwy)  / det;
   
   fSigmaY2 = 0;
-  for (Int_t i = 0; i < nTimeBins+1; i++) {    
+  for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
     if (!fUsable[i]) continue;
     Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
     fSigmaY2 += delta*delta;
   }
   fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
+       // TEMPORARY UNTIL covariance properly calculated
+       fSigmaY2 = TMath::Max(fSigmaY2, Float_t(.1));
   
   fZfitR[0]  = (sumwx2 * sumwz  - sumwx * sumwxz) / det;
   fZfitR[1]  = (sumw   * sumwxz - sumwx * sumwz)  / det;
@@ -489,7 +497,9 @@ void AliTRDseed::Update()
   fYfitR[1] += fYref[1];
   fYfit[0]   = fYfitR[0];
   fYfit[1]   = fYfitR[1];
-    
+
+       //printf("y0 = %7.3f tgy = %7.3f z0 = %7.3f tgz = %7.3f \n", fYfitR[0], fYfitR[1], fZfitR[0], fZfitR[1]);
+
   UpdateUsed();
 
 }
@@ -501,17 +511,11 @@ void AliTRDseed::UpdateUsed()
   // Update used seed
   //
 
-  AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-  Int_t nTimeBins = cal->GetNumberOfTimeBins();
-
   fNUsed = 0;
-  for (Int_t i = 0; i < nTimeBins; i++) {
-    if (!fClusters[i]) {
-      continue;
-    }
-    if ((fClusters[i]->IsUsed())) {
-      fNUsed++;
-    }
+  for (Int_t i = 0; i < fgTimeBins; i++) {
+    if (!fClusters[i]) continue;
+               if(!fUsable[i]) continue;   
+    if ((fClusters[i]->IsUsed())) fNUsed++;
   }
 
 }
@@ -527,9 +531,6 @@ Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
   TLinearFitter fitterT2(4,"hyp4");  
   fitterT2.StoreData(kTRUE);
        
-  AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-  Int_t nTimeBins = cal->GetNumberOfTimeBins();
-
   Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
   
   Int_t npointsT = 0;
@@ -540,7 +541,7 @@ Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
     if (!cseed[iLayer].IsOK()) continue;
     Double_t tilt = cseed[iLayer].fTilt;
 
-    for (Int_t itime = 0; itime < nTimeBins+1; itime++) {
+    for (Int_t itime = 0; itime < fgTimeBins+1; itime++) {
 
       if (!cseed[iLayer].fUsable[itime]) continue;
       // x relative to the midle chamber
index 11e6e69bb2975923ebdcbfd5ee42e409ef9446d6..ff2bff417c060b3a42e73550a3960d2b305483d6 100644 (file)
@@ -1,6 +1,5 @@
 #ifndef ALITRDSEED_H
-#define ALITRDSEED_H   
-
+#define ALITRDSEED_H
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */ 
 
@@ -23,7 +22,7 @@ class AliTRDseed : public TObject {
 
   AliTRDseed(); 
   AliTRDseed(const AliTRDseed &s);
-  ~AliTRDseed() {};                 
+  ~AliTRDseed() {};
 
   AliTRDseed      &operator=(const AliTRDseed &s)           { *(new(this) AliTRDseed(s)); 
                                                               return *this;          }
@@ -64,7 +63,10 @@ class AliTRDseed : public TObject {
           Float_t  GetCC() const                            { return fCC;            }
           Float_t  GetChi2() const                          { return fChi2;          }
           Float_t  GetChi2Z() const                         { return fChi2Z;         }
+          Int_t    GetNTimeBins() const                   { return fgTimeBins; }
+          Int_t    GetNTimeBinsRange() const              { return fTimeBinsRange; }
 
+                               
           void     SetTilt(Float_t tilt)                    { fTilt        = tilt;   }
           void     SetPadLength(Float_t len)                { fPadLength   = len;    }
           void     SetX0(Float_t x0)                        { fX0          = x0;     }
@@ -92,44 +94,49 @@ class AliTRDseed : public TObject {
           void     SetCC(Float_t cc)                        { fCC          = cc;     }
           void     SetChi2(Float_t chi2)                    { fChi2        = chi2;   }
           void     SetChi2Z(Float_t chi2z)                  { fChi2Z       = chi2z;  }
+  static  void     SetNTimeBins(Int_t nTB)                { fgTimeBins  = nTB; }
+          void     SetNTimeBinsRange(Int_t nTB)           { fTimeBinsRange  = nTB; }
 
  protected:
 
           void     Copy(TObject &o) const;
           
-          Float_t  fTilt;               //  Tilting angle
-          Float_t  fPadLength;          //  Pad length
-          Float_t  fX0;                 //  X0 position
-          Float_t  fX[knTimebins];      //! X position
-          Float_t  fY[knTimebins];      //! Y position
-          Float_t  fZ[knTimebins];      //! Z position
-          Int_t    fIndexes[knTimebins];//! Indexes
-          AliTRDcluster *fClusters[knTimebins]; // Clusters
-          Bool_t   fUsable[knTimebins]; //! Indication  - usable cluster
-          Float_t  fYref[2];            //  Reference y
-          Float_t  fZref[2];            //  Reference z
-          Float_t  fYfit[2];            //  Y fit position +derivation
-          Float_t  fYfitR[2];           //  Y fit position +derivation
-          Float_t  fZfit[2];            //  Z fit position
-          Float_t  fZfitR[2];           //  Z fit position
-          Float_t  fSigmaY;             //  "Robust" sigma in Y - constant fit
-          Float_t  fSigmaY2;            //  "Robust" sigma in Y - line fit
-          Float_t  fMeanz;              //  Mean vaue of z
-          Float_t  fZProb;              //  Max probbable z
-          Int_t    fLabels[2];          //  Labels
-          Int_t    fN;                  //  Number of associated clusters
-          Int_t    fN2;                 //  Number of not crossed
-          Int_t    fNUsed;              //  Number of used clusters
-          Int_t    fFreq;               //  Frequency
-          Int_t    fNChange;            //  Change z counter
-          Float_t  fMPads;              //  Mean number of pads per cluster
-
-          Float_t  fC;                  //  Curvature
-          Float_t  fCC;                 //  Curvature with constrain
-          Float_t  fChi2;               //  Global chi2
-          Float_t  fChi2Z;              //  Global chi2
-
-  ClassDef(AliTRDseed,1)                //  Seed for a local TRD track
+  static  Int_t          fgTimeBins;            //  Local copy of the total number of TB
+          Int_t          fTimeBinsRange;        //  Number of time bins in the geometrical range of the detector
+          Int_t          fTimeBin0;             //  Time bin corresponding to t0
+          Float_t        fTilt;                 //  Tilting angle
+          Float_t        fPadLength;            //  Pad length
+          Float_t        fX0;                   //  X0 position
+          Float_t        fX[knTimebins];        //! X position
+          Float_t        fY[knTimebins];        //! Y position
+          Float_t        fZ[knTimebins];        //! Z position
+          Int_t          fIndexes[knTimebins];  //! Indexes
+          AliTRDcluster *fClusters[knTimebins]; //  Clusters
+          Bool_t         fUsable[knTimebins];   //  Indication  - usable cluster
+          Float_t        fYref[2];              //  Reference y
+          Float_t        fZref[2];              //  Reference z
+          Float_t        fYfit[2];              //  Y fit position +derivation
+          Float_t        fYfitR[2];             //  Y fit position +derivation
+          Float_t        fZfit[2];              //  Z fit position
+          Float_t        fZfitR[2];             //  Z fit position
+          Float_t        fSigmaY;               //  "Robust" sigma in Y - constant fit
+          Float_t        fSigmaY2;              //  "Robust" sigma in Y - line fit
+          Float_t        fMeanz;                //  Mean vaue of z
+          Float_t        fZProb;                //  Max probbable z
+          Int_t          fLabels[2];            //  Labels
+          Int_t          fN;                    //  Number of associated clusters
+          Int_t          fN2;                   //  Number of not crossed
+          Int_t          fNUsed;                //  Number of used clusters
+          Int_t          fFreq;                 //  Frequency
+          Int_t          fNChange;              //  Change z counter
+          Float_t        fMPads;                //  Mean number of pads per cluster
+
+          Float_t        fC;                    //  Curvature
+          Float_t        fCC;                   //  Curvature with constrain
+          Float_t        fChi2;                 //  Global chi2
+          Float_t        fChi2Z;                //  Global chi2
+
+  ClassDef(AliTRDseed,1)                        //  Seed for a local TRD track
 
 };
 
index 2cbcca783698d570513751e98a55e6a93698f906..450f8ffca360dd5aa3dcc0365ca8a0ea1c7a338b 100644 (file)
 
 #include "AliTRDseedV1.h"
 #include "AliTRDcluster.h"
+#include "AliTRDtrack.h"
 #include "AliTRDcalibDB.h"
 #include "AliTRDstackLayer.h"
 #include "AliTRDrecoParam.h"
+#include "AliTRDgeometry.h"
+#include "Cal/AliTRDCalPID.h"
 
 #define SEED_DEBUG
 
@@ -44,27 +47,28 @@ ClassImp(AliTRDseedV1)
 //____________________________________________________________________
 AliTRDseedV1::AliTRDseedV1(Int_t layer, AliTRDrecoParam *p) 
   :AliTRDseed()
-  ,fLayer(layer)
-  ,fTimeBins(0)
+  ,fPlane(layer)
   ,fOwner(kFALSE)
+  ,fMom(0.)
   ,fRecoParam(p)
 {
   //
   // Constructor
   //
-
-       //AliInfo("");
-       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-       fTimeBins = cal->GetNumberOfTimeBins();
-
+       for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = 0.;
+       for(int itb=0; itb < knTimebins; itb++){
+               fdQdl[itb]  = 0.;
+               fdQ[itb]    = 0.;
+       }
+       for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec]  = -1.;
 }
 
 //____________________________________________________________________
-AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref, Bool_t owner)
+AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
   :AliTRDseed((AliTRDseed&)ref)
-  ,fLayer(ref.fLayer)
-  ,fTimeBins(ref.fTimeBins)
+  ,fPlane(ref.fPlane)
   ,fOwner(kFALSE)
+  ,fMom(ref.fMom)
   ,fRecoParam(ref.fRecoParam)
 {
   //
@@ -72,17 +76,16 @@ AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref, Bool_t owner)
   //
 
        //AliInfo("");
-
-       if(owner){
-               for(int ic=0; ic<fTimeBins; ic++){
-                       if(!fClusters[ic]) continue;
-                       fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
-               }
-               fOwner = kTRUE;
+       if(ref.fOwner) SetOwner();
+       for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = ref.fdEdx[islice];
+       for(int itb=0; itb < knTimebins; itb++){ 
+               fdQdl[itb] = ref.fdQdl[itb];
+               fdQ[itb]   = ref.fdQ[itb];
        }
-
+       for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = ref.fProb[ispec];
 }
 
+
 //____________________________________________________________________
 AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
 {
@@ -107,7 +110,13 @@ AliTRDseedV1::~AliTRDseedV1()
 
        //AliInfo(Form("fOwner[%s]", fOwner?"YES":"NO"));
 
-       if(fOwner) delete [] fClusters;
+       if(fOwner) 
+               for(int itb=0; itb<knTimebins; itb++){
+                       if(!fClusters[itb]) continue; 
+                       //AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
+                       delete fClusters[itb];
+                       fClusters[itb] = 0x0;
+               }
 }
 
 //____________________________________________________________________
@@ -120,12 +129,136 @@ void AliTRDseedV1::Copy(TObject &ref) const
        //AliInfo("");
        AliTRDseedV1 &target = (AliTRDseedV1 &)ref; 
        
-       target.fLayer     = fLayer;
-       target.fTimeBins  = fTimeBins;
-       target.fRecoParam = fRecoParam;
+       target.fPlane         = fPlane;
+       target.fMom           = fMom;
+       target.fRecoParam     = fRecoParam;
+       
+       for(int islice=0; islice < knSlices; islice++) target.fdEdx[islice] = fdEdx[islice];
+       for(int itb=0; itb < knTimebins; itb++){
+               target.fdQdl[itb] = fdQdl[itb];
+               target.fdQ[itb]   = fdQ[itb];
+       }
+       for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) target.fProb[ispec] = fProb[ispec];
+       
        AliTRDseed::Copy(target);
 }
 
+
+//____________________________________________________________
+void AliTRDseedV1::Init(AliTRDtrack *track)
+{
+// Initialize this tracklet using the track information
+//
+// Parameters:
+//   track - the TRD track used to initialize the tracklet
+// 
+// Detailed description
+// The function sets the starting point and direction of the
+// tracklet according to the information from the TRD track.
+// 
+// Caution
+// The TRD track has to be propagated to the beginning of the
+// chamber where the tracklet will be constructed
+//
+
+       Double_t y, z; 
+       track->GetProlongation(fX0, y, z);
+       fYref[0] = y;
+       fYref[1] = track->GetSnp() < 0. ? track->GetTgl() : -track->GetTgl();
+       fZref[0] = z;
+       // TO DO 
+       // tilting pad correction !!
+       fZref[1] = 0.; // TMath::Tan(track->Theta());
+
+       //printf("Tracklet ref x[%7.3f] y[%7.3f] z[%7.3f], snp[%f] tgl[%f]\n", fX0, fYref[0], fZref[0], track->GetSnp(), track->GetTgl());
+}
+
+//____________________________________________________________________
+Double_t* AliTRDseedV1::GetProbability()
+{      
+// Fill probability array for tracklet from the DB.
+//
+// Parameters
+//
+// Output
+//   returns pointer to the probability array and 0x0 if missing DB access 
+//
+// Detailed description
+
+       
+       // retrive calibration db
+  AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
+  if (!calibration) {
+    AliError("No access to calibration data");
+    return 0x0;
+  }
+
+  // Retrieve the CDB container class with the parametric detector response
+  const AliTRDCalPID *pd = calibration->GetPIDObject(fRecoParam->GetPIDMethod());
+  if (!pd) {
+    AliError("No access to AliTRDCalPID object");
+    return 0x0;
+  }
+       
+       // calculate tracklet length TO DO
+  Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
+  /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
+  
+  //calculate dE/dx
+  CookdEdx(fRecoParam->GetNdEdxSlices());
+  
+  // Sets the a priori probabilities
+  for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
+    fProb[ispec] = pd->GetProbability(ispec, fMom, &fdEdx[0], length, fPlane); 
+  }
+
+       return &fProb[0];
+}
+
+//____________________________________________________________________
+void AliTRDseedV1::CookdEdx(Int_t nslices)
+{
+// Calculates average dE/dx for all slices and store them in the internal array fdEdx. 
+//
+// Parameters:
+//  nslices : number of slices for which dE/dx should be calculated
+// Output:
+//
+// Detailed description
+// Calculates average dE/dx for all slices. Depending on the PID methode 
+// the number of slices can be 3 (LQ) or 8(NN). The calculation is based 
+// on previously calculated quantities dQ/dl of each cluster. The 
+// following effects are included in the calculation:
+// 1. calibration values for t0 and vdrift
+// 2. cluster sharing (optional see AliTRDrecoParam::SetClusterSharing())
+//
+
+       Int_t nclusters[knSlices];
+       for(int i=0; i<knSlices; i++){ 
+               fdEdx[i]     = 0.;
+               nclusters[i] = 0;
+       }
+       
+       AliTRDcluster *cluster = 0x0;
+       for(int ic=0; ic<fgTimeBins; ic++){
+               if(!(cluster = fClusters[ic])) continue;
+               Int_t tb = cluster->GetLocalTimeBin();
+               
+               // consider calibration effects
+               if(tb < fTimeBin0 || tb >= fTimeBin0+fTimeBinsRange) continue;
+       
+               // consider cluster sharing ... TO DO
+               //if(fRecoParam->GetClusterSharing() && cluster->GetSharing()) continue;
+               
+               Int_t slice = (tb-fTimeBin0)*nslices/fTimeBinsRange;
+               fdEdx[slice]   += fdQdl[ic];
+               nclusters[slice]++;
+       } // End of loop over clusters
+
+       // calculate mean charge per slice
+       for(int is=0; is<nslices; is++) if(nclusters[is]) fdEdx[is] /= nclusters[is];
+}
+
 //____________________________________________________________________
 Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
 {
@@ -140,6 +273,53 @@ Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
                + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength;
 }
 
+//____________________________________________________________________
+void AliTRDseedV1::GetCovAt(Double_t /*x*/, Double_t *cov) const
+{
+// Computes covariance in the y-z plane at radial point x
+
+       const Float_t k0= .2; // to be checked in FindClusters
+       Double_t sy20   = k0*TMath::Tan(fYfit[1]); sy20 *= sy20;
+       
+       Double_t sy2    = fSigmaY2*fSigmaY2 + sy20;
+       Double_t sz2    = fPadLength/12.;
+
+       //printf("Yfit[1] %f sy20 %f SigmaY2 %f\n", fYfit[1], sy20, fSigmaY2);
+
+       cov[0] = sy2;
+       cov[1] = fTilt*(sy2-sz2);
+       cov[2] = sz2;
+}
+
+//____________________________________________________________________
+void AliTRDseedV1::SetdQdl(Double_t length)
+{
+       for(int ic=0; ic<fgTimeBins; ic++) fdQdl[ic] = fdQ[ic] *length;
+}
+
+//____________________________________________________________________
+void AliTRDseedV1::SetOwner(Bool_t own)
+{
+       //AliInfo(Form("own [%s] fOwner[%s]", own?"YES":"NO", fOwner?"YES":"NO"));
+       
+       if(own){
+               for(int ic=0; ic<knTimebins; ic++){
+                       if(!fClusters[ic]) continue;
+                       fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
+               }
+               fOwner = kTRUE;
+       } else {
+               if(fOwner){
+                       for(int ic=0; ic<knTimebins; ic++){
+                               if(!fClusters[ic]) continue;
+                               delete fClusters[ic];
+                               //fClusters[ic] = tracker->GetClusters(index) TODO
+                       }
+               }
+               fOwner = kFALSE;
+       }
+}
+
 //____________________________________________________________________
 Bool_t AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
                                        , Float_t quality
@@ -156,7 +336,9 @@ Bool_t      AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
                AliError("Seed can not be used without a valid RecoParam.");
                return kFALSE;
        }
-       
+
+       //AliInfo(Form("TimeBins = %d TimeBinsRange = %d", fgTimeBins, fTimeBinsRange));
+
        Float_t  tquality;
        Double_t kroady = fRecoParam->GetRoad1y();
        Double_t kroadz = fPadLength * .5 + 1.;
@@ -171,7 +353,7 @@ Bool_t      AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
        for (Int_t iter = 0; iter < niter; iter++) {
        //AliInfo(Form("iter = %i", iter));
                ncl = 0;
-               for (Int_t iTime = 0; iTime < fTimeBins; iTime++) {
+               for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
                        // define searching configuration
                        Double_t dxlayer = layer[iTime].GetX() - fX0;
                        if(c){
@@ -188,11 +370,14 @@ Bool_t    AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
 //                     printf("xexp = %3.3f ,yexp = %3.3f, zexp = %3.3f\n",layer[iTime].GetX(),yexp,zexp);
 //                     printf("layer[%i].GetNClusters() = %i\n", iTime, layer[iTime].GetNClusters());
                        Int_t    index = layer[iTime].SearchNearestCluster(yexp, zexp, kroady, kroadz);
+
+//                     printf("%d[%d] x[%7.3f | %7.3f] y[%7.3f] z[%7.3f]\n", iTime, layer[iTime].GetNClusters(), dxlayer, layer[iTime].GetX(), yexp, zexp);
 //                     for(Int_t iclk = 0; iclk < layer[iTime].GetNClusters(); iclk++){
 //                             AliTRDcluster *testcl = layer[iTime].GetCluster(iclk);
-//                             printf("Cluster %i: x = %3.3f, y = %3.3f, z = %3.3f\n",iclk,testcl->GetX(), testcl->GetY(), testcl->GetZ());
+//                             printf("Cluster %i: %d x = %7.3f, y = %7.3f, z = %7.3f\n", iclk, testcl->GetLocalTimeBin(), testcl->GetX(), testcl->GetY(), testcl->GetZ());
 //                     }
 //                     printf("Index = %i\n",index);
+
                        if (index < 0) continue;
                        
                        // Register cluster
@@ -206,7 +391,8 @@ Bool_t      AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
                        fX[iTime]        = dxlayer;
                        fY[iTime]        = cl->GetY();
                        fZ[iTime]        = cl->GetZ();
-       
+                       fdQ[iTime]       = cl->GetQ()/layer[iTime].GetdX();
+                       
                        // Debugging
                        ncl++;
                }
@@ -214,10 +400,10 @@ Bool_t    AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
 #ifdef SEED_DEBUG
 //             Int_t nclusters = 0;
 //             Float_t fD[iter] = 0.;
-//             for(int ic=0; ic<fTimeBins+1; ic++){
+//             for(int ic=0; ic<fgTimeBins+1; ic++){
 //                     AliTRDcluster *ci = fClusters[ic];
 //                     if(!ci) continue;
-//                     for(int jc=ic+1; jc<fTimeBins+1; jc++){
+//                     for(int jc=ic+1; jc<fgTimeBins+1; jc++){
 //                             AliTRDcluster *cj = fClusters[jc];
 //                             if(!cj) continue;
 //                             fD[iter] += TMath::Sqrt((ci->GetY()-cj->GetY())*(ci->GetY()-cj->GetY())+
@@ -245,10 +431,8 @@ Bool_t     AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
 }
 
 //____________________________________________________________________
-Bool_t AliTRDseedV1::AttachClustersProj(AliTRDstackLayer *layer
-                                       , Float_t /*quality*/
-                                       , Bool_t kZcorr
-                                       , AliTRDcluster *c)
+Bool_t AliTRDseedV1::AttachClusters(AliTRDstackLayer *layer
+                                       ,Bool_t kZcorr)
 {
   //
   // Projective algorithm to attach clusters to seeding tracklets
@@ -271,8 +455,7 @@ Bool_t      AliTRDseedV1::AttachClustersProj(AliTRDstackLayer *layer
                return kFALSE;
        }
 
-       const Int_t knTimeBins = 35;
-       const Int_t kClusterCandidates = 2 * knTimeBins;
+       const Int_t kClusterCandidates = 2 * knTimebins;
        
        //define roads
        Double_t kroady = fRecoParam->GetRoad1y();
@@ -282,13 +465,13 @@ Bool_t    AliTRDseedV1::AttachClustersProj(AliTRDstackLayer *layer
 
        // working variables
        AliTRDcluster *clusters[kClusterCandidates];
-       Double_t cond[4], yexp[knTimeBins], zexp[knTimeBins],
+       Double_t cond[4], yexp[knTimebins], zexp[knTimebins],
                yres[kClusterCandidates], zres[kClusterCandidates];
-       Int_t ncl, *index = 0x0, tboundary[knTimeBins];
+       Int_t ncl, *index = 0x0, tboundary[knTimebins];
        
        // Do cluster projection
        Int_t nYclusters = 0; Bool_t kEXIT = kFALSE;
-       for (Int_t iTime = 0; iTime < fTimeBins; iTime++) {
+       for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
                fX[iTime] = layer[iTime].GetX() - fX0;
                zexp[iTime] = fZref[0] + fZref[1] * fX[iTime];
                yexp[iTime] = fYref[0] + fYref[1] * fX[iTime] - zcorr;
@@ -298,7 +481,7 @@ Bool_t      AliTRDseedV1::AttachClustersProj(AliTRDstackLayer *layer
                cond[2] = zexp[iTime] - kroadz; cond[3] = zexp[iTime] + kroadz;
                layer[iTime].GetClusters(cond, index, ncl);
                for(Int_t ic = 0; ic<ncl; ic++){
-                       c = layer[iTime].GetCluster(index[ic]);
+                       AliTRDcluster *c = layer[iTime].GetCluster(index[ic]);
                        clusters[nYclusters] = c;
                        yres[nYclusters++] = c->GetY() - yexp[iTime];
                        if(nYclusters >= kClusterCandidates) {
@@ -338,9 +521,10 @@ Bool_t     AliTRDseedV1::AttachClustersProj(AliTRDstackLayer *layer
        // Select only one cluster/TimeBin
        Int_t lastCluster = 0;
        fN2 = 0;
-       for (Int_t iTime = 0; iTime < fTimeBins; iTime++) {
+       for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
                ncl = tboundary[iTime] - lastCluster;
                if(!ncl) continue;
+               AliTRDcluster *c = 0x0;
                if(ncl == 1){
                        c = clusters[lastCluster];
                } else if(ncl > 1){
@@ -356,32 +540,42 @@ Bool_t    AliTRDseedV1::AttachClustersProj(AliTRDstackLayer *layer
                        }
                        c = clusters[iptr];
                }
-               //Int_t globalIndex = layer[iTime].GetGlobalIndex(index);
-               //fIndexes[iTime]  = globalIndex;
+               //Int_t GlobalIndex = layer[iTime].GetGlobalIndex(index);
+               //fIndexes[iTime]  = GlobalIndex;
                fClusters[iTime] = c;
                fY[iTime]        = c->GetY();
                fZ[iTime]        = c->GetZ();
-               lastCluster = tboundary[iTime];
+               fdQ[iTime]       = c->GetQ()/layer[iTime].GetdX();
+               lastCluster      = tboundary[iTime];
                fN2++;
        }
        
        // number of minimum numbers of clusters expected for the tracklet
-       Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fTimeBins);
+       Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fgTimeBins);
   if (fN2 < kClmin){
                AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
     fN2 = 0;
     return kFALSE;
   }
-       AliTRDseed::Update();
+
+       // update used clusters
+       fNUsed = 0;
+       for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
+               if(!fClusters[iTime]) continue;
+               if((fClusters[iTime]->IsUsed())) fNUsed++;
+       }
+
+  if (fN2-fNUsed < kClmin){
+               AliWarning(Form("Too many clusters already in use %d (from %d).", fNUsed, fN2));
+    fN2 = 0;
+    return kFALSE;
+  }
        
-//     // fit tracklet and update clusters
-//     if(!FitTracklet()) return kFALSE;
-//     UpdateUsed();
        return kTRUE;
 }
 
 //____________________________________________________________________
-Bool_t AliTRDseedV1::FitTracklet()
+Bool_t AliTRDseedV1::Fit()
 {
   //
   // Linear fit of the tracklet
@@ -402,20 +596,20 @@ Bool_t AliTRDseedV1::FitTracklet()
   Float_t  anglecor = fTilt * fZref[1];  // Correction to the angle
 
        // calculate residuals
-       const Int_t knTimeBins = 35;
-       Float_t yres[knTimeBins]; // y (r-phi) residuals
-       Int_t zint[knTimeBins],   // Histograming of the z coordinate 
-             zout[2*knTimeBins];//
+       Float_t yres[knTimebins]; // y (r-phi) residuals
+       Int_t zint[knTimebins],   // Histograming of the z coordinate
+             zout[2*knTimebins];//
        
        fN = 0;
-       for (Int_t iTime = 0; iTime < fTimeBins; iTime++) {
+       for (Int_t iTime = 0; iTime < fTimeBinsRange; iTime++) {
     if (!fClusters[iTime]) continue;
     yres[iTime] = fY[iTime] - fYref[0] - (fYref[1] + anglecor) * fX[iTime];
-               zint[fN++] = Int_t(fZ[iTime]);
+               zint[fN] = Int_t(fZ[iTime]);
+               fN++;
        }
 
        // calculate pad row boundary crosses
-       Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fTimeBins);
+       Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fTimeBinsRange);
        Int_t nz = AliMathBase::Freq(fN, zint, zout, kFALSE);
   fZProb   = zout[0];
   if(nz <= 1) zout[3] = 0;
@@ -437,7 +631,8 @@ Bool_t AliTRDseedV1::FitTracklet()
        Int_t npads;
   fMPads = 0;
        fMeanz = 0.;
-       for(int iTime=0; iTime<fTimeBins; iTime++){
+       // we will use only the clusters which are in the detector range
+       for(int iTime=0; iTime<fTimeBinsRange; iTime++){
     fUsable[iTime] = kFALSE;
     if (!fClusters[iTime]) continue;
                npads = fClusters[iTime]->GetNPads();
@@ -476,7 +671,7 @@ Bool_t AliTRDseedV1::FitTracklet()
   fYfitR[1]    = (sumw   * sumwxy - sumwx * sumwy)  / det;
   
   fSigmaY2 = 0;
-  for (Int_t i = 0; i < fTimeBins+1; i++) {
+  for (Int_t i = 0; i < fTimeBinsRange+1; i++) {
     if (!fUsable[i]) continue;
     Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
     fSigmaY2 += delta*delta;
@@ -647,13 +842,13 @@ void AliTRDseedV1::Print()
        printf("  fX0        = %f\n", fX0);
        for(int ic=0; ic<nTimeBins; ic++) {
           const Char_t *isUsable = fUsable[ic]?"Yes":"No";
-         printf("  %d X[%f] Y[%f] Z[%f] Indexes[%d] clusters[%#x] usable[%s]\n"
+         printf("  %d X[%f] Y[%f] Z[%f] Indexes[%d] clusters[%p] usable[%s]\n"
                 , ic
                 , fX[ic]
                 , fY[ic]
                 , fZ[ic]
                 , fIndexes[ic]
-                , ((void *) fClusters[ic])
+                , ((void*) fClusters[ic])
                 , isUsable);
         }
 
index 8a3b1958cfd48ae9afd93ce24d6ea66327f4fe73..0e428d5633bd53ee8686e74cbaa811e54150f803 100644 (file)
 #include "AliTRDseed.h"
 #endif
 
+#ifndef ALIPID_H
+#include "AliPID.h"
+#endif
+
 #ifndef ALIRIEMAN_H
 #include "AliRieman.h"
 #endif
@@ -26,35 +30,54 @@ class AliRieman;
 class AliTRDstackLayer;
 class AliTRDcluster;
 class AliTRDrecoParam;
+class AliTRDtrack;
 
 class AliTRDseedV1 : public AliTRDseed
 {
 
- public:
+  public:
+
+       enum {
+         knSlices = 10
+       };
 
        AliTRDseedV1(Int_t layer = -1, AliTRDrecoParam *p=0x0);
        ~AliTRDseedV1();
-       AliTRDseedV1(const AliTRDseedV1 &ref, Bool_t owner=kFALSE);
+       AliTRDseedV1(const AliTRDseedV1 &ref);
        AliTRDseedV1& operator=(const AliTRDseedV1 &ref);
 
-       Bool_t  AttachClustersIter(AliTRDstackLayer *layer, Float_t quality, Bool_t kZcorr = kFALSE, AliTRDcluster *c=0x0);
-       Bool_t  AttachClustersProj(AliTRDstackLayer *layer, Float_t quality, Bool_t kZcorr = kFALSE, AliTRDcluster *c=0x0);
+       Bool_t  AttachClustersIter(AliTRDstackLayer *layer, Float_t quality, Bool_t kZcorr = kFALSE
+                                 , AliTRDcluster *c=0x0);
+       Bool_t  AttachClusters(AliTRDstackLayer *layer, Bool_t kZcorr = kFALSE);
+       void    CookdEdx(Int_t nslices);
        static  Float_t FitRiemanTilt(AliTRDseedV1 * cseed, Bool_t terror);
-       Bool_t  FitTracklet();
-       
-//             Bool_t  AttachClusters(Double_t *dx, Float_t quality, Bool_t kZcorr=kFALSE, AliTRDcluster *c=0x0);
-       inline Float_t GetChi2Z(const Float_t z = 0.) const;
-       inline Float_t GetChi2Y(const Float_t y = 0.) const;
-              Float_t GetQuality(Bool_t kZcorr) const;
-              Int_t   GetLayer() const                       { return fLayer;    }
-
-       inline void    Update(const AliRieman *rieman);
-               void    Print(Option_t * /*o*/) const          { }
-              void    Print();
+       Bool_t  Fit();
 
-              void    SetLayer(Int_t l)                      { fLayer     = l;   }
-              void    SetNTimeBins(Int_t nTB)                { fTimeBins  = nTB; }
-              void    SetRecoParam(AliTRDrecoParam *p)       { fRecoParam = p;   }
+              void      Init(AliTRDtrack *track);
+       inline void      Init(const AliRieman *fit);
+       
+       inline Float_t   GetChi2Z(const Float_t z = 0.) const;
+       inline Float_t   GetChi2Y(const Float_t y = 0.) const;
+              void      GetCovAt(Double_t x, Double_t *cov) const;
+              Float_t*  GetdEdx() {return &fdEdx[0];}
+              Double_t* GetdQdl() {return &fdQdl[0];}
+              Double_t  GetMomentum() const {return fMom;}
+              Int_t     GetN() const {return fN2;}
+              Float_t   GetQuality(Bool_t kZcorr) const;
+              Int_t     GetPlane() const                       { return fPlane;    }
+              Double_t* GetProbability();
+              Double_t  GetYat(Double_t x) const {return fYfitR[0] + fYfitR[1] * (x - fX0);}
+              Double_t  GetZat(Double_t x) const {return fZfitR[0] + fZfitR[1] * (x - fX0);}
+                                
+              Bool_t    IsOwner() const {return fOwner;}
+         void      Print(Option_t * /*o*/) const          { }
+              void      Print();
+              
+              void      SetdQdl(Double_t length);
+              void      SetMomentum(Double_t mom) {fMom = mom;}
+              void      SetOwner(Bool_t own = kTRUE);
+              void      SetPlane(Int_t p)                      { fPlane     = p;   }
+              void      SetRecoParam(AliTRDrecoParam *p)       { fRecoParam = p;   }
 
  protected:
 
@@ -62,12 +85,16 @@ class AliTRDseedV1 : public AliTRDseed
 
  private:
 
-       Int_t            fLayer;     //  layer for this seed
-       Int_t            fTimeBins;  //  local copy of the DB info
-       Bool_t           fOwner;     //  owner of the clusters
-       AliTRDrecoParam *fRecoParam; //! local copy of the reco params 
+       Int_t            fPlane;                  //  TRD plane
+       Bool_t           fOwner;                  //  Toggle ownership of clusters
+       Float_t          fMom;                    //  Momentum estimate for tracklet [GeV/c]
+       Float_t          fdEdx[knSlices];         //  dE/dx measurements for tracklet
+       Double_t         fdQdl[knTimebins];       //  dQ/dl for all clusters attached to tracklet 
+       Double_t         fdQ[knTimebins];         //! dQ for all clusters attached to tracklet TODO migrate to AliTRDcluster
+        Double_t         fProb[AliPID::kSPECIES]; //  PID probabilities
+       AliTRDrecoParam *fRecoParam;              //! Local copy of the reco params 
 
-       ClassDef(AliTRDseedV1, 1)    //  New TRD seed 
+       ClassDef(AliTRDseedV1, 1)                 //  New TRD seed 
 
 };
 
@@ -88,7 +115,7 @@ inline Float_t AliTRDseedV1::GetChi2Y(const Float_t y) const
 }
 
 //____________________________________________________________
-inline void AliTRDseedV1::Update(const AliRieman *rieman)
+inline void AliTRDseedV1::Init(const AliRieman *rieman)
 {
        fZref[0] = rieman->GetZat(fX0);
        fZref[1] = rieman->GetDZat(fX0);
index 8570e29461079fa5c1d1f1b1ea37267d21540bc5..8371049c4cc82766eaa8333f0b1637d26985e368 100644 (file)
@@ -41,6 +41,7 @@ AliTRDtrack::AliTRDtrack()
   ,fSeedLab(-1)
   ,fdEdx(0)
   ,fDE(0)
+  ,fPIDquality(0)
   ,fClusterOwner(kFALSE)
   ,fPIDmethod(kLQ)
   ,fStopped(kFALSE)
@@ -90,6 +91,7 @@ AliTRDtrack::AliTRDtrack(AliTRDcluster *c, Int_t index
   ,fSeedLab(-1)
   ,fdEdx(0)
   ,fDE(0)
+  ,fPIDquality(0)
   ,fClusterOwner(kFALSE)
   ,fPIDmethod(kLQ)
   ,fStopped(kFALSE)
@@ -170,6 +172,7 @@ AliTRDtrack::AliTRDtrack(const AliTRDtrack &t/*, const Bool_t owner*/)
   ,fSeedLab(t.GetSeedLabel())
   ,fdEdx(t.fdEdx)
   ,fDE(t.fDE)
+  ,fPIDquality(t.fPIDquality)
   ,fClusterOwner(kTRUE)
   ,fPIDmethod(t.fPIDmethod)
   ,fStopped(t.fStopped)
@@ -225,7 +228,8 @@ AliTRDtrack::AliTRDtrack(const AliTRDtrack &t/*, const Bool_t owner*/)
     fBudget[i]      = t.fBudget[i];
   }
 
-}                                
+       for(Int_t ispec = 0; ispec<AliPID::kSPECIES; ispec++) fPID[ispec] = t.fPID[ispec];
+}
 
 //_____________________________________________________________________________
 AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t /*alpha*/) 
@@ -233,6 +237,7 @@ AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t /*alpha*/)
   ,fSeedLab(-1)
   ,fdEdx(t.GetPIDsignal())
   ,fDE(0)
+  ,fPIDquality(0)
   ,fClusterOwner(kFALSE)
   ,fPIDmethod(kLQ)
   ,fStopped(kFALSE)
@@ -277,7 +282,7 @@ AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t /*alpha*/)
     fBudget[i]      = 0;
   }
 
-}              
+}
 
 //_____________________________________________________________________________
 AliTRDtrack::AliTRDtrack(const AliESDtrack &t)
@@ -285,6 +290,7 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack &t)
   ,fSeedLab(-1)
   ,fdEdx(t.GetTRDsignal())
   ,fDE(0)
+  ,fPIDquality(0)
   ,fClusterOwner(kFALSE)
   ,fPIDmethod(kLQ)
   ,fStopped(kFALSE)
@@ -876,8 +882,6 @@ Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index
     return kFALSE;
   }
 
-  AliTracker::FillResiduals(this,p,cov,c->GetVolumeId());
-
   // Register cluster to track
   Int_t n      = GetNumberOfClusters();
   fIndex[n]    = index;
index c705f28bf9eb5d781d1df501cdaac937fb68b0fb..8feeb17031787eaafde649d3fb2ec4b2c27eb213 100644 (file)
 
 #include "AliTRDtracklet.h"
 
+#ifndef ALITRDSEEDV1_H
+#include "AliTRDseedV1.h"
+#endif
+
 class AliESDtrack;
 class AliTrackReference;
 class AliTRDcluster;
-
 class AliTRDtrack : public AliKalmanTrack {
 
  public:
@@ -44,7 +47,7 @@ class AliTRDtrack : public AliKalmanTrack {
                    , Double_t xr, Double_t alpha);
          AliTRDtrack(const AliTRDtrack &t/*, const Bool_t owner = kTRUE*/);    
          AliTRDtrack(const AliESDtrack &t);
       ~AliTRDtrack();
virtual ~AliTRDtrack();
          AliTRDtrack(const AliKalmanTrack &t, Double_t alpha); 
 
          void            ResetClusters()                              { SetChi2(0.0); 
@@ -87,6 +90,13 @@ class AliTRDtrack : public AliKalmanTrack {
          Float_t         GetBudget(Int_t i) const                     { return fBudget[i];                   }
          Float_t         GetChi2Last() const                          { return fChi2Last;                    }
          AliTRDtrack    *GetBackupTrack()                             { return fBackupTrack;                 }
+         // dummy to bridge the function in AliTRDtrackV1
+        //Int_t          GetNumberOfClusters() const                   { printf("AliTRDtrack::GetNumberOfClusters()\n"); 
+         //                                                               return AliKalmanTrack::GetNumberOfClusters();   }
+ inline virtual Int_t    GetNumberOfTracklets() const;
+ virtual Int_t           GetTrackletIndex(Int_t plane) const          { return plane>=0 && plane<6 
+                                                                             ? fTrackletIndex[plane] : -1;   } 
+
 
          void            SetdEdx(Double_t dedx)                       { fdEdx                      = dedx;   }
          void            SetStop(Bool_t stop)                         { fStopped                   = stop;   }
@@ -163,8 +173,8 @@ class AliTRDtrack : public AliKalmanTrack {
          Float_t  fDE;                                //  Integrated delta energy
          Float_t  fdEdxPlane[kNplane][kNslice];       //  dE/dx from all 6 planes in 3 slices each
          Int_t    fTimBinPlane[kNplane];              //  Time bin of Max cluster from all 6 planes
+         UChar_t  fPIDquality;                        //  No of planes used for PID calculation        
          Double_t fPID[AliPID::kSPECIES];             //  PID probabilities
-
         Float_t  fMom[kNplane];                      //  Track momentum at chamber entrance
         Float_t  fSnp[kNplane];                      //  Track direction
         Float_t  fTgl[kNplane];                      //  Track direction
@@ -188,9 +198,21 @@ class AliTRDtrack : public AliKalmanTrack {
   AliTRDtracklet  fTracklets[6];                      //  Tracklets
          Float_t  fBudget[3];                         //  Integrated material budget
   AliTRDtrack    *fBackupTrack;                       //! Backup track
+       
+        Int_t    fTrackletIndex[6];                  //  Tracklets index in the tracker list
+  AliTRDseedV1    fTracklet[6];                       //  Tracklets array defining the track
 
   ClassDef(AliTRDtrack,9)                             //  TRD reconstructed tracks
 
 };                     
 
+//___________________________________________________________
+inline Int_t AliTRDtrack::GetNumberOfTracklets() const
+{
+       Int_t ntrklt = 0;
+       for(int ip=0; ip<6; ip++) if(fTrackletIndex[ip] >= 0) ntrklt++;
+       return ntrklt;
+}
+
+
 #endif   
index 93d0b04a0ff1972e5fc39299f6703f0657396739..38f5fe63689503d780c151a58ced2d7d0f42684b 100644 (file)
 
 /* $Id$ */
 
-///////////////////////////////////////////////////////////////////////////////
-//                                                                           //
-//  TRD track                                                                //
-//                                                                           //
-//  Authors:                                                                 //
-//    Alex Bercuci <A.Bercuci@gsi.de>                                        //
-//    Markus Fasel <M.Fasel@gsi.de>                                          //
-//                                                                           //
-///////////////////////////////////////////////////////////////////////////////
-
 #include "AliTRDtrackV1.h"
 #include "AliTRDcluster.h"
 #include "AliTRDcalibDB.h"
 
 ClassImp(AliTRDtrackV1)
 
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Represents a reconstructed TRD track                                     //
+//  Local TRD Kalman track                                                   //
+//                                                                           //
+//  Authors:                                                                 //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                        //
+//    Markus Fasel <M.Fasel@gsi.de>                                          //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
 
 //_______________________________________________________________
-AliTRDtrackV1::AliTRDtrackV1() :
-       AliTRDtrack()
-       ,fRecoParam(0x0)
-{      
+AliTRDtrackV1::AliTRDtrackV1() 
+  :AliTRDtrack()
+  ,fRecoParam(0x0)
+{
   //
   // Default constructor
   //
-
+       
        for(int ip=0; ip<6; ip++){
                fTrackletIndex[ip] = -1;
                fTracklet[ip].Reset();
@@ -51,12 +51,12 @@ AliTRDtrackV1::AliTRDtrackV1() :
 }
 
 //_______________________________________________________________
-AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) :
-       AliTRDtrack(t)
-       ,fRecoParam(0x0)
+AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) 
+  :AliTRDtrack(t)
+  ,fRecoParam(0x0)
 {
   //
-  // Standard constructor
+  // Constructor from AliESDtrack
   //
 
        //AliInfo(Form("alpha %f", GetAlpha()));
@@ -65,14 +65,18 @@ AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) :
 }
 
 //_______________________________________________________________
-AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &t) :
-       AliTRDtrack(t)
-       ,fRecoParam(0x0)
+AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref) 
+  :AliTRDtrack(ref)
+  ,fRecoParam(ref.fRecoParam)
 {
   //
   // Copy constructor
   //
 
+       for(int ip=0; ip<6; ip++){ 
+               fTrackletIndex[ip] = ref.fTrackletIndex[ip];
+               fTracklet[ip]      = ref.fTracklet[ip];
+       }
 }
 
 //_______________________________________________________________
@@ -80,21 +84,23 @@ AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &t) :
 // {
 //     
 // }
-
        
 //_______________________________________________________________
-AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Double_t cov[15], Double_t x, Double_t alpha) :
-       AliTRDtrack()
-       ,fRecoParam(0x0)
+AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5]
+                           , const Double_t cov[15]
+                           , Double_t x, Double_t alpha)
+  :AliTRDtrack()
+  ,fRecoParam(0x0)
 {
-// The stand alone tracking constructor
-// TEMPORARY !!!!!!!!!!!
-// to check :
-// 1. covariance matrix
-// 2. dQdl calculation
-
+  //
+  // The stand alone tracking constructor
+  // TEMPORARY !!!!!!!!!!!
+  // to check :
+  // 1. covariance matrix
+  // 2. dQdl calculation
+  //
 
-       Double_t cnv   = 1.0 / (GetBz() * kB2C);
+  Double_t cnv   = 1.0 / (GetBz() * kB2C);
 
   Double_t pp[5] = { p[0]    
                    , p[1]
@@ -145,41 +151,75 @@ AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Do
 //_______________________________________________________________
 Bool_t AliTRDtrackV1::CookPID()
 {
+  //
+  // Cook the PID information
+  //
+
 // CookdEdx();  // truncated mean ... do we still need it ?
 
 // CookdEdxTimBin(seed->GetID());
-       for(int itrklt=0; itrklt<6; itrklt++){
-    for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) fdEdxPlane[itrklt][iSlice] = -1.;
+       
+  // Sets the a priori probabilities
+  for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
+    fPID[ispec] = 1.0 / AliPID::kSPECIES;      
+  }
+       
+       // steer PID calculation @ tracklet level
+       Double_t *prob = 0x0;
+       fPIDquality = 0;
+       for(int itrklt=0; itrklt<AliESDtrack::kNPlane; itrklt++){
+    //for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) fdEdxPlane[itrklt][iSlice] = -1.;
 
                if(fTrackletIndex[itrklt]<0) continue;
-               fTracklet[itrklt].CookdEdx(fdEdxPlane[itrklt]);
+               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
+       
+
+  // normalize probabilities
+  Double_t probTotal = 0.0;
+  for (Int_t is = 0; is < AliPID::kSPECIES; is++) probTotal += fPID[is];
 
-       // retrive calibration db
-  AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
-  if (!calibration) {
-    AliError("No access to calibration data");
-    return kFALSE;
-  }
 
-  // Retrieve the CDB container class with the parametric detector response
-  const AliTRDCalPID *pd = calibration->GetPIDObject(0/*fRecoParam->GetPIDMethod()*/);
-  if (!pd) {
-    AliError("No access to AliTRDCalPID object");
+  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;
   }
 
-// CookPID(pidQ);
-
+  for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) fPID[iSpecies] /= probTotal;
 
+  return kTRUE;
 }
 
+//_______________________________________________________________
+Float_t AliTRDtrackV1::GetMomentum(Int_t plane) const
+{
+  //
+  // Get the momentum at a given plane
+  //
+
+       return plane >=0 && plane < 6 && fTrackletIndex[plane] >= 0 ? fTracklet[plane].GetMomentum() : -1.;
+}
 
 //_______________________________________________________________
 Double_t AliTRDtrackV1::GetPredictedChi2(const AliTRDseedV1 *trklt) const
 {
   //
-  // Returns the predicted chi2
+  // Get the predicted chi2
   //
 
   Double_t x      = trklt->GetX0();
@@ -192,11 +232,39 @@ Double_t AliTRDtrackV1::GetPredictedChi2(const AliTRDseedV1 *trklt) const
 
 }
 
+//_______________________________________________________________
+Bool_t AliTRDtrackV1::IsOwner() const
+{
+  //
+  // Check whether track owns the tracklets
+  //
+
+       for (Int_t ip = 0; ip < AliESDtrack::kNPlane; ip++) {
+               if(fTrackletIndex[ip] < 0) continue;
+               if(!fTracklet[ip].IsOwner()) return kFALSE;
+       }
+       return kTRUE;
+}
+       
+//_______________________________________________________________
+void AliTRDtrackV1::SetOwner(Bool_t own)
+{
+  //
+  // Toggle ownership of tracklets
+  //
+
+       for (Int_t ip = 0; ip < AliESDtrack::kNPlane; ip++) {
+               if(fTrackletIndex[ip] < 0) continue;
+               //AliInfo(Form("p[%d] index[%d]", ip, fTrackletIndex[ip]));
+               fTracklet[ip].SetOwner(own);
+       }
+}
+
 //_______________________________________________________________
 void AliTRDtrackV1::SetTracklet(AliTRDseedV1 *trklt, Int_t plane, Int_t index)
 {
   //
-  // Sets the tracklets
+  // Set the tracklets
   //
 
        if(plane < 0 || plane >=6) return;
@@ -205,13 +273,13 @@ void AliTRDtrackV1::SetTracklet(AliTRDseedV1 *trklt, Int_t plane, Int_t index)
 }
 
 //_______________________________________________________________
-Bool_t  AliTRDtrackV1::Update(const  AliTRDseedV1 *trklt, Double_t chisq)
+Bool_t  AliTRDtrackV1::Update(AliTRDseedV1 *trklt, Double_t chisq)
 {
   //
-  // Update the track
+  // Update track parameters
   //
 
-       Double_t x      = trklt->GetX0();
+  Double_t x      = trklt->GetX0();
   Double_t p[2]   = { trklt->GetYat(x)
                     , trklt->GetZat(x) };
   Double_t cov[3];
@@ -229,21 +297,32 @@ Bool_t  AliTRDtrackV1::Update(const  AliTRDseedV1 *trklt, Double_t chisq)
 //   fClusters[n] = c;
   SetNumberOfClusters(GetNumberOfClusters()+trklt->GetN());
   SetChi2(GetChi2() + chisq);
-
+       
+       // update tracklet
+       trklt->SetMomentum(GetP());
+  Double_t s = GetSnp(), t = GetTgl();
+       trklt->SetdQdl(TMath::Sqrt((1.0 - s*s) / (1.0 + t*t)));
        return kTRUE;
 }
 
 //_______________________________________________________________
-void AliTRDtrackV1::UpdateESDdEdx(AliESDtrack *track)
+void AliTRDtrackV1::UpdateESDtrack(AliESDtrack *track)
 {
   //
-  // Update the dedx info
+  // Update the ESD track
   //
-
-       for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
-               for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
-                       track->SetTRDsignals(fdEdxPlane[i][j], i, j);
-               }
-               track->SetTRDTimBin(fTimBinPlane[i], i);
+       
+       // copy dEdx to ESD
+       Float_t *dedx = 0x0;
+       for (Int_t ip = 0; ip < AliESDtrack::kNPlane; ip++) {
+               if(fTrackletIndex[ip] < 0) continue;
+               fTracklet[ip].CookdEdx(AliESDtrack::kNSlice);
+               dedx = fTracklet[ip].GetdEdx();
+               for (Int_t js = 0; js < AliESDtrack::kNSlice; js++) track->SetTRDsignals(dedx[js], ip, js);
+               //track->SetTRDTimBin(fTimBinPlane[i], i);
        }
+
+       // copy PID to ESD
+       track->SetTRDpid(fPID);
+       track->SetTRDpidQuality(fPIDquality);
 }
index 52a003306b6152b9785fe455e34984e7c8a7df9b..90b173699e5ccea85dc3e6be37dee53ddb82e01c 100644 (file)
@@ -5,15 +5,11 @@
 
 /* $Id$ */
 
-////////////////////////////////////////////////////////////////////////////
-//                                                                        //
-//  The TRD track                                                         //
-//                                                                        //
-//  Authors:                                                              //
-//    Alex Bercuci <A.Bercuci@gsi.de>                                     //
-//    Markus Fasel <M.Fasel@gsi.de>                                       //
-//                                                                        //
-////////////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Represents a reconstructed TRD track                                     //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
 
 #ifndef ALITRDTRACK_H
 #include "AliTRDtrack.h"
 
 class AliTRDseedV1;
 class AliESDtrack;
-class AliCluster;
 
 class AliTRDtrackV1 : public AliTRDtrack
 {
-public:
+
+ 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*/){ return *this; }
+       AliTRDtrackV1(const AliTRDtrackV1 &ref);
+       AliTRDtrackV1 &operator=(const AliTRDtrackV1 &ref) { *(new(this) AliTRDtrackV1(ref));
+                                                             return *this; }
 
-       Bool_t        CookPID();
-       inline Int_t  GetNumberOfClusters() const;
-       Double_t      GetPredictedChi2(const AliTRDseedV1 *tracklet) const;
-        Double_t      GetPredictedChi2(const AliCluster* /*c*/) const                   { return 0.0; }
+       Bool_t         CookPID();
+       Float_t        GetMomentum(Int_t plane) const;
+       inline Int_t   GetNumberOfClusters() 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];}
-       void          SetTracklet(AliTRDseedV1 *trklt, Int_t plane, Int_t index);
-       Bool_t        Update(const  AliTRDseedV1 *tracklet, Double_t chi2);
-virtual Bool_t        Update(const AliCluster*, Double_t, Int_t) { return kTRUE; }
-       void          UpdateESDdEdx(AliESDtrack *t);
+       Int_t*         GetTrackletIndexes() {return &fTrackletIndex[0];}
+       
+       Bool_t         IsOwner() const;
+       
+       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*, Double_t, Int_t, Double_t) { return kFALSE; };
+        Bool_t         Update(const AliCluster *, Double_t, Int_t)             { return kFALSE; };
+       void           UpdateESDtrack(AliESDtrack *t);
 
-protected:
+ protected:
        AliTRDrecoParam *fRecoParam;       // reconstruction parameters
+
        ClassDef(AliTRDtrackV1, 1)         // development TRD track
+
 };
 
 //___________________________________________________________
index 288f59c1203065eb9976f13ee2e9cad7ba53a4bb..77d3f7eb7dc7dee24fec70df103dbdb0af79118f 100644 (file)
@@ -73,7 +73,7 @@ class AliTRDtracker : public AliTracker {
   Float_t          GetMinClustersInTrack() const  { return fgkMinClustersInTrack; }
   Int_t            GetLastPlane(AliTRDtrack *track);
   Double_t         GetTiltFactor(const AliTRDcluster *c);
-  Bool_t           GetTrackPoint(Int_t index, AliTrackPoint& p) const;
+  virtual Bool_t   GetTrackPoint(Int_t index, AliTrackPoint& p) const;
   Double_t         GetX(Int_t sec, Int_t plane, Int_t localTB) const;
   Double_t         GetX(Int_t sec, Int_t pl) const
                                                   { return fTrSec[sec]->GetLayer(pl)->GetX();            }
@@ -183,14 +183,13 @@ class AliTRDtracker : public AliTracker {
   TH1D                    *fHDeltaX;                       // QA histogram
   TH1D                    *fHXCl;                          // QA histogram
 
- private:
+ protected:
 
   Int_t            FollowProlongation(AliTRDtrack &t);
   Int_t            PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxStep);
   Double_t         ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) const;
   Double_t         ExpectedSigmaZ2(Double_t r, Double_t tgl) const;
 
- private:  
  
   TTreeSRedirector        *fDebugStreamer;                 //!Debug streamer
   
index 37441a904ae4dd616c2bee75b49dbe03bc6a8370..a976537b6aedd8fc7671e29cd7d18fcd8eae6ca7 100644 (file)
@@ -64,6 +64,8 @@
 #include "AliTRDstackLayer.h"
 #include "AliTRDrecoParam.h"
 #include "AliTRDseedV1.h"
+#include "AliTRDtrackV1.h"
+#include "Cal/AliTRDCalDet.h"
 
 #define DEBUG
 
@@ -78,9 +80,9 @@ Double_t AliTRDtrackerV1::fgTopologicQA[kNConfigs] = {
 AliTRDtrackerV1::AliTRDtrackerV1(AliTRDrecoParam *p) 
   :AliTRDtracker()
   ,fSieveSeeding(0)
+  ,fTracklets(0x0)
   ,fRecoParam(p)
   ,fFitter(0x0)
-  ,fDebugStreamerV1(0x0)
 {
   //
   // Default constructor. Nothing is initialized.
@@ -92,9 +94,9 @@ AliTRDtrackerV1::AliTRDtrackerV1(AliTRDrecoParam *p)
 AliTRDtrackerV1::AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p) 
   :AliTRDtracker(in)
   ,fSieveSeeding(0)
+  ,fTracklets(0x0)
   ,fRecoParam(p)
   ,fFitter(0x0)
-  ,fDebugStreamerV1(0x0)
 {
   //
   // Standard constructor.
@@ -105,8 +107,7 @@ AliTRDtrackerV1::AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p)
        fFitter = new AliTRDtrackerFitter();
 
 #ifdef DEBUG
-       fDebugStreamerV1 = new TTreeSRedirector("TRDdebug.root");
-       fFitter->SetDebugStream(fDebugStreamerV1);
+       fFitter->SetDebugStream(fDebugStreamer);
 #endif
 
 }
@@ -118,10 +119,9 @@ AliTRDtrackerV1::~AliTRDtrackerV1()
   // Destructor
   //
 
-       if(fDebugStreamerV1) delete fDebugStreamerV1;
        if(fFitter) delete fFitter;
        if(fRecoParam) delete fRecoParam;
-
+       if(fTracklets) {fTracklets->Delete(); delete fTracklets;}
 }
 
 //____________________________________________________________________
@@ -157,6 +157,734 @@ Int_t AliTRDtrackerV1::Clusters2Tracks(AliESDEvent *esd)
        return ntracks;
 }
 
+
+//_____________________________________________________________________________
+Bool_t AliTRDtrackerV1::GetTrackPoint(Int_t /*index*/, AliTrackPoint &/*p*/) const
+{
+       //AliInfo(Form("Asking for tracklet %d", index));
+       
+       if(index<0) return kFALSE;
+       //AliTRDseedV1 *tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index);
+       // etc
+       return kTRUE;
+}
+
+
+//_____________________________________________________________________________
+Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) 
+{
+  //
+  // Gets seeds from ESD event. The seeds are AliTPCtrack's found and
+  // backpropagated by the TPC tracker. Each seed is first propagated 
+  // to the TRD, and then its prolongation is searched in the TRD.
+  // If sufficiently long continuation of the track is found in the TRD
+  // the track is updated, otherwise it's stored as originaly defined 
+  // by the TPC tracker.   
+  //  
+
+       Int_t   found    = 0;     // number of tracks found
+       Float_t foundMin = 20.0;
+       
+       AliTRDseed::SetNTimeBins(fTimeBinsPerPlane);
+       Int_t    nSeed   = event->GetNumberOfTracks();
+       if(!nSeed){
+               // run stand alone tracking
+               if (AliTRDReconstructor::SeedingOn()) Clusters2Tracks(event);
+               return 0;
+       }
+       
+       Float_t *quality = new Float_t[nSeed];
+       Int_t   *index   = new Int_t[nSeed];
+       for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
+               AliESDtrack *seed = event->GetTrack(iSeed);
+               Double_t covariance[15];
+               seed->GetExternalCovariance(covariance);
+               quality[iSeed] = covariance[0] + covariance[2];
+       }
+       // Sort tracks according to covariance of local Y and Z
+       TMath::Sort(nSeed,quality,index,kFALSE);
+       
+       // Backpropagate all seeds
+       for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
+       
+               // Get the seeds in sorted sequence
+               AliESDtrack *seed = event->GetTrack(index[iSeed]);
+       
+               // Check the seed status
+               ULong_t status = seed->GetStatus();
+               if ((status & AliESDtrack::kTPCout) == 0) continue;
+               if ((status & AliESDtrack::kTRDout) != 0) continue;
+       
+               // Do the back prolongation
+               Int_t   lbl         = seed->GetLabel();
+               AliTRDtrackV1 *track  = new AliTRDtrackV1(*seed);
+               //track->Print();
+               track->SetSeedLabel(lbl);
+               seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); // Make backup
+               fNseeds++;
+               Float_t p4          = track->GetC();
+               Int_t   expectedClr = FollowBackProlongation(*track);
+               //AliInfo(Form("\nTRACK %d Clusters %d [%d] in chi2 %f", index[iSeed], expectedClr, track->GetNumberOfClusters(), track->GetChi2()));
+               //track->Print();
+
+               //Double_t cov[15];
+               //seed->GetExternalCovariance(cov);
+               //AliInfo(Form("track %d cov[%f %f] 0", index[iSeed], cov[0], cov[2]));
+
+               if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
+                               (track->Pt() > 0.8)) {
+                       //
+                       // Make backup for back propagation
+                       //
+                       Int_t foundClr = track->GetNumberOfClusters();
+                       if (foundClr >= foundMin) {
+                               //AliInfo(Form("Making backup track ncls [%d]...", foundClr));
+                               track->CookdEdx();
+                               track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
+                               CookLabel(track,1 - fgkLabelFraction);
+                               if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
+                               
+                               
+               //seed->GetExternalCovariance(cov);
+               //AliInfo(Form("track %d cov[%f %f] 0 test", index[iSeed], cov[0], cov[2]));
+
+                               // Sign only gold tracks
+                               if (track->GetChi2() / track->GetNumberOfClusters() < 4) {
+                                       if ((seed->GetKinkIndex(0)      ==   0) &&
+                                                       (track->Pt()                <  1.5)) UseClusters(track);
+                               }
+                               Bool_t isGold = kFALSE;
+       
+                               // Full gold track
+                               if (track->GetChi2() / track->GetNumberOfClusters() < 5) {
+                                       if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
+
+                                       isGold = kTRUE;
+                               }
+               //seed->GetExternalCovariance(cov);
+               //AliInfo(Form("track %d cov[%f %f] 00", index[iSeed], cov[0], cov[2]));
+       
+                               // Almost gold track
+                               if ((!isGold)  && (track->GetNCross() == 0) &&
+                                               (track->GetChi2() / track->GetNumberOfClusters()  < 7)) {
+                                       //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
+                                       if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
+                                       
+                                       isGold = kTRUE;
+                               }
+               //seed->GetExternalCovariance(cov);
+               //AliInfo(Form("track %d cov[%f %f] 01", index[iSeed], cov[0], cov[2]));
+                               
+                               if ((!isGold) && (track->GetBackupTrack())) {
+                                       if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) && ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) {
+                                               seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
+                                               isGold = kTRUE;
+                                       }
+                               }
+               //seed->GetExternalCovariance(cov);
+               //AliInfo(Form("track %d cov[%f %f] 02", index[iSeed], cov[0], cov[2]));
+       
+                               //if ((track->StatusForTOF() > 0) && (track->GetNCross() == 0) && (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected())  > 0.4)) {
+                                       //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
+                               //}
+                       }
+               }
+               /**/
+       
+               /**/
+               // Debug part of tracking
+/*             TTreeSRedirector &cstream = *fDebugStreamer;
+               Int_t eventNrInFile = event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
+               if (AliTRDReconstructor::StreamLevel() > 0) {
+                       if (track->GetBackupTrack()) {
+                               cstream << "Tracks"
+                               << "EventNrInFile="  << eventNrInFile
+                               << "ESD.="     << seed
+                               << "trd.="     << track
+                               << "trdback.=" << track->GetBackupTrack()
+                               << "\n";
+                       }
+                       else {
+                               cstream << "Tracks"
+                               << "EventNrInFile="  << eventNrInFile
+                               << "ESD.="     << seed
+                               << "trd.="     << track
+                               << "trdback.=" << track
+                               << "\n";
+                       }
+               }*/
+               /**/
+       
+               //seed->GetExternalCovariance(cov);
+               //AliInfo(Form("track %d cov[%f %f] 1", index[iSeed], cov[0], cov[2]));
+
+               // Propagation to the TOF (I.Belikov)
+               if (track->GetStop() == kFALSE) {
+                       //AliInfo("Track not stopped in TRD ...");
+                       Double_t xtof  = 371.0;
+                       Double_t xTOF0 = 370.0;
+               
+                       Double_t c2    = track->GetSnp() + track->GetC() * (xtof - track->GetX());
+                       if (TMath::Abs(c2) >= 0.99) {
+                               delete track;
+                               continue;
+                       }
+                       
+                       PropagateToX(*track,xTOF0,fgkMaxStep);
+       
+                       // Energy losses taken to the account - check one more time
+                       c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX());
+                       if (TMath::Abs(c2) >= 0.99) {
+                               delete track;
+                               continue;
+                       }
+                       
+                       //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
+                       //      fHBackfit->Fill(7);
+                       //delete track;
+                       //      continue;
+                       //}
+       
+                       Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
+                       Double_t y;
+                       track->GetYAt(xtof,GetBz(),y);
+                       if (y >  ymax) {
+                               if (!track->Rotate( AliTRDgeometry::GetAlpha())) {
+                                       delete track;
+                                       continue;
+                               }
+                       }else if (y < -ymax) {
+                               if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
+                                       delete track;
+                                       continue;
+                               }
+                       }
+                                       
+                       if (track->PropagateTo(xtof)) {
+                               //AliInfo("set kTRDout");
+                               seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
+       
+                               for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
+                                       for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
+                                               seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
+                                       }
+                                       seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
+                               }
+                               //seed->SetTRDtrack(new AliTRDtrack(*track));
+                               if (track->GetNumberOfClusters() > foundMin) found++;
+                       }
+               } else {
+                       //AliInfo("Track stopped in TRD ...");
+                       
+                       if ((track->GetNumberOfClusters() >              15) &&
+                                       (track->GetNumberOfClusters() > 0.5*expectedClr)) {
+                               seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
+       
+                               //seed->SetStatus(AliESDtrack::kTRDStop);
+                               for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
+                                       for (Int_t j = 0; j <AliESDtrack::kNSlice; j++) {
+                                               seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
+                                       }
+                                       seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
+                               }
+                               //seed->SetTRDtrack(new AliTRDtrack(*track));
+                               found++;
+                       }
+               }
+
+               //if (((t->GetStatus()&AliESDtrack::kTRDout)!=0 )
+       
+               seed->SetTRDQuality(track->StatusForTOF());
+               seed->SetTRDBudget(track->GetBudget(0));
+               
+//             if ((seed->GetStatus()&AliESDtrack::kTRDin)!=0 ) printf("TRDin ");      
+//             if ((seed->GetStatus()&AliESDtrack::kTRDbackup)!=0 ) printf("TRDbackup ");      
+//             if ((seed->GetStatus()&AliESDtrack::kTRDStop)!=0 ) printf("TRDstop ");  
+//             if ((seed->GetStatus()&AliESDtrack::kTRDout)!=0 ) printf("TRDout ");    
+//             printf("\n");
+               delete track;
+
+               //seed->GetExternalCovariance(cov);
+               //AliInfo(Form("track %d cov[%f %f] 2", index[iSeed], cov[0], cov[2]));
+       }
+       
+
+       AliInfo(Form("Number of seeds: %d",fNseeds));
+       AliInfo(Form("Number of back propagated TRD tracks: %d",found));
+               
+       //fSeeds->Clear(); 
+       fNseeds = 0;
+       
+       delete [] index;
+       delete [] quality;
+       
+  return 0;
+}
+
+
+//____________________________________________________________________
+Int_t AliTRDtrackerV1::RefitInward(AliESDEvent *event)
+{
+  //
+  // Refits tracks within the TRD. The ESD event is expected to contain seeds 
+  // at the outer part of the TRD. 
+  // The tracks are propagated to the innermost time bin 
+  // of the TRD and the ESD event is updated
+  // Origin: Thomas KUHR (Thomas.Kuhr@cern.ch)
+  //
+
+  Int_t   nseed    = 0; // contor for loaded seeds
+  Int_t   found    = 0; // contor for updated TRD tracks
+  AliTRDtrackV1 track;
+  for (Int_t itrack = 0; itrack < event->GetNumberOfTracks(); itrack++) {
+    AliESDtrack *seed = event->GetTrack(itrack);
+               new(&track) AliTRDtrackV1(*seed);
+
+    if (track.GetX() < 270.0) {
+      seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup);
+      //AliInfo(Form("Remove for X = %7.3f [270.]\n", track.GetX()));
+                       continue;
+    }
+
+    ULong_t status = seed->GetStatus();
+    if((status & AliESDtrack::kTRDout) == 0) continue;
+    if((status & AliESDtrack::kTRDin)  != 0) continue;
+    nseed++; 
+
+    track.ResetCovariance(50.0);
+
+               // do the propagation and processing
+    FollowProlongation(track);
+    // computes PID for track
+    track.CookPID();
+    // update calibration references using this track
+               //track.Calibrate();
+
+               // Prolongate to TPC
+    Double_t xTPC = 250.0;
+    if (PropagateToX(track, xTPC, fgkMaxStep)) { //  -with update
+      seed->UpdateTrackParams(&track, AliESDtrack::kTRDrefit);
+      track.UpdateESDtrack(seed);
+       // Add TRD track to ESDfriendTrack
+                       if (AliTRDReconstructor::StreamLevel() > 0 /*&& quality TODO*/){ 
+                               AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track);
+                               calibTrack->SetOwner();
+                               seed->AddCalibObject(calibTrack);
+                       }
+                       found++;
+    } else {  // - without update
+      AliTRDtrackV1 tt(*seed);
+      if (PropagateToX(tt, xTPC, fgkMaxStep)) seed->UpdateTrackParams(&tt, AliESDtrack::kTRDrefit);
+    }
+  }
+  AliInfo(Form("Number of loaded seeds: %d",nseed));
+  AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
+  
+       return 0;
+}
+
+
+//____________________________________________________________________
+Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
+{
+// Extrapolates the TRD track in the TPC direction.
+//
+// Parameters
+//   t : the TRD track which has to be extrapolated
+// 
+// Output
+//   number of clusters attached to the track
+//
+// Detailed description
+//
+// Starting from current radial position of track <t> this function
+// extrapolates the track through the 6 TRD layers. The following steps
+// are being performed for each plane:
+// 1. prepare track:
+//   a. get plane limits in the local x direction
+//   b. check crossing sectors 
+//   c. check track inclination
+// 2. search tracklet in the tracker list (see GetTracklet() for details)
+// 3. evaluate material budget using the geo manager
+// 4. propagate and update track using the tracklet information.
+//
+// Debug level 2
+//
+  
+       //AliInfo("");
+       Int_t    nClustersExpected = 0;
+       Int_t lastplane = 5; //GetLastPlane(&t);
+       for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
+               //AliInfo(Form("plane %d", iplane));
+    Int_t row1 = GetGlobalTimeBin(0, iplane, 0); // to be modified to the true time bin in the geometrical acceptance
+               //AliInfo(Form("row1 %d", row1));
+
+    // Propagate track close to the plane if neccessary
+               AliTRDpropagationLayer *layer = fTrSec[0]->GetLayer(row1);
+               Double_t currentx  = layer->GetX();
+    if (currentx < (-fgkMaxStep + t.GetX())) 
+      if (!PropagateToX(t, currentx+fgkMaxStep, fgkMaxStep)) break;
+
+    if (!AdjustSector(&t)) break;
+     
+               Int_t row0    = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
+               //AliInfo(Form("row0 %d", row0));
+
+    // Start global position
+    Double_t xyz0[3];
+    t.GetXYZ(xyz0);
+
+               // End global position
+    Double_t x = fTrSec[0]->GetLayer(row0)->GetX(), y, z;
+    if (!t.GetProlongation(x,y,z)) break;    
+    Double_t xyz1[3];
+    xyz1[0] =  x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
+    xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
+    xyz1[2] =  z;
+               
+    // Get material budget
+    Double_t param[7];
+    AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
+    Double_t xrho= param[0]*param[4];
+    Double_t xx0 = param[1]; // Get mean propagation parameters
+
+    // Propagate and update
+    //Int_t sector = t.GetSector();
+    Int_t   index   = 0;
+               //AliInfo(Form("sector %d", sector));
+    AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
+         //AliInfo(Form("tracklet 0x%x @ %d", tracklet, index));
+               if(!tracklet) continue;
+
+               t.PropagateTo(tracklet->GetX0(), xx0, xrho); // not correct
+         if (!AdjustSector(&t)) break;
+         
+    Double_t maxChi2 = t.GetPredictedChi2(tracklet);
+         if (maxChi2 < 1e+10 && t.Update(tracklet, maxChi2)){ 
+               nClustersExpected += tracklet->GetN();
+       }
+  }
+
+#ifdef DEBUG
+       if(AliTRDReconstructor::StreamLevel() > 1){
+               Int_t index;
+               for(int iplane=0; iplane<6; iplane++){
+                       AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
+                       if(!tracklet) continue;
+                       t.SetTracklet(tracklet, iplane, index);
+               }
+
+               TTreeSRedirector &cstreamer = *fDebugStreamer;
+               cstreamer << "FollowProlongation"
+                       << "ncl="      << nClustersExpected
+                       << "track.="   << &t
+                       << "\n";
+       }
+#endif
+
+  return nClustersExpected;
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
+{
+// Extrapolates the TRD track in the TOF direction.
+//
+// Parameters
+//   t : the TRD track which has to be extrapolated
+// 
+// Output
+//   number of clusters attached to the track
+//
+// Detailed description
+//
+// Starting from current radial position of track <t> this function
+// extrapolates the track through the 6 TRD layers. The following steps
+// are being performed for each plane:
+// 1. prepare track:
+//   a. get plane limits in the local x direction
+//   b. check crossing sectors 
+//   c. check track inclination
+// 2. build tracklet (see AliTRDseed::AttachClusters() for details)
+// 3. evaluate material budget using the geo manager
+// 4. propagate and update track using the tracklet information.
+//
+// Debug level 2
+//
+
+       Int_t nClustersExpected = 0;
+
+  // Loop through the TRD planes
+  for (Int_t iplane = 0; iplane < AliTRDgeometry::Nplan(); iplane++) {
+               //AliInfo(Form("Processing plane %d ...", iplane));
+               // Get the global time bin for the first local time bin of the given plane
+    Int_t row0 = GetGlobalTimeBin(0, iplane, fTimeBinsPerPlane-1);
+
+               // Retrive first propagation layer in the chamber
+               AliTRDpropagationLayer *layer = fTrSec[0]->GetLayer(row0);
+
+               // Get the X coordinates of the propagation layer for the first time bin
+    Double_t currentx = layer->GetX(); // what if X is not defined ???
+    if (currentx < t.GetX()) continue;
+    
+               // Get the global time bin for the last local time bin of the given plane
+               Int_t row1 = GetGlobalTimeBin(0, iplane, 0);
+    
+               // Propagate closer to the current chamber if neccessary 
+    if (currentx > (fgkMaxStep + t.GetX()) && !PropagateToX(t, currentx-fgkMaxStep, fgkMaxStep)) break;
+   
+               // Rotate track to adjacent sector if neccessary
+    if (!AdjustSector(&t)) break;
+               Int_t sector =  Int_t(TMath::Abs(t.GetAlpha()/AliTRDgeometry::GetAlpha()));
+               if(t.GetAlpha() < 0) sector = AliTRDgeometry::Nsect() - sector-1;
+               
+               //AliInfo(Form("sector %d [%f]", sector, t.GetAlpha()));
+
+               // Check whether azimuthal angle is getting too large
+    if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) break;
+   
+    //Calculate global entry and exit positions of the track in chamber (only track prolongation)
+    Double_t xyz0[3]; // entry point 
+               t.GetXYZ(xyz0);   
+               //printf("Entry global x[%7.3f] y[%7.3f] z[%7.3f]\n", xyz0[0], xyz0[1], xyz0[2]);
+                               
+               // Get local Y and Z at the X-position of the end of the chamber
+               Double_t x0 = fTrSec[sector]->GetLayer(row1)->GetX(), y, z;
+               if (!t.GetProlongation(x0, y, z)) break;
+               //printf("Exit  local x[%7.3f] y[%7.3f] z[%7.3f]\n", x0, y, z);
+
+               Double_t xyz1[3]; // exit point
+               xyz1[0] =  x0 * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha()); 
+    xyz1[1] = +x0 * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
+    xyz1[2] =  z;
+
+               //printf("Exit  global x[%7.3f] y[%7.3f] z[%7.3f]\n", xyz1[0], xyz1[1], xyz1[2]);
+               // Find tracklet along the path inside the chamber
+    AliTRDseedV1 tracklet(*t.GetTracklet(iplane));
+               // if the track is not already build (e.g. stand alone tracker) we build it now.
+               if(!tracklet.GetN()){ // a better check has to be implemented TODO!!!!!!!
+               
+                       //AliInfo(Form("Building tracklet for plane %d ...", iplane));
+                       // check if we are inside detection volume
+                       Int_t ichmb = fGeom->GetChamber(xyz0[2], iplane); 
+                       if(ichmb<0) ichmb = fGeom->GetChamber(xyz1[2], iplane);
+                       if(ichmb<0){
+                               // here we should decide what to do with the track. The space between the pads in 2 chambers is 4cm+. Is it making sense to continue building the tracklet here TODO????
+                               AliWarning(Form("Track prolongated in the interspace between TRD detectors in plane %d. Skip plane. To be fixed !", iplane));
+                               continue;
+                       }
+       
+                       // temporary until the functionalities of AliTRDpropagationLayer and AliTRDstackLayer are merged TODO
+                       AliTRDpadPlane *pp = fGeom->GetPadPlane(iplane, ichmb);
+                       Int_t nrows = pp->GetNrows();
+                       Double_t stacklength = pp->GetRow0ROC() - pp->GetRowEndROC();/*(nrows - 2) * pp->GetLengthIPad()        + 2 * pp->GetLengthOPad() + (nrows - 1) * pp->GetRowSpacing();*/
+                       Double_t z0  = fGeom->GetRow0(iplane, ichmb, 0);
+
+                       Int_t nClustersChmb = 0;
+                       AliTRDstackLayer stackLayer[35];
+                       for(int itb=0; itb<fTimeBinsPerPlane; itb++){
+                               const AliTRDpropagationLayer ksmLayer(*(fTrSec[sector]->GetLayer(row1 - itb)));
+                               stackLayer[itb] = ksmLayer;
+#ifdef DEBUG
+                               stackLayer[itb].SetDebugStream(fDebugStreamer);
+#endif                 
+                               stackLayer[itb].SetRange(z0 - stacklength, stacklength);
+                               stackLayer[itb].SetSector(sector);
+                               stackLayer[itb].SetStackNr(ichmb);
+                               stackLayer[itb].SetNRows(nrows);
+                               stackLayer[itb].SetRecoParam(fRecoParam);
+                               stackLayer[itb].BuildIndices();
+                               nClustersChmb += stackLayer[itb].GetNClusters();
+                       }
+                       //AliInfo(Form("Detector p[%d] c[%d]. Building tracklet from %d clusters ... ", iplane, ichmb, nClustersChmb));
+
+                       tracklet.SetRecoParam(fRecoParam);
+                       tracklet.SetTilt(TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle()));
+                       tracklet.SetPadLength(pp->GetLengthIPad());
+                       tracklet.SetPlane(iplane);
+                       Int_t tbRange   = fTimeBinsPerPlane; //Int_t(AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght() * AliTRDCommonParam::Instance()->GetSamplingFrequency()/AliTRDcalibDB::Instance()->GetVdriftDet()->GetValue(det));
+                       //printf("%d hl[%f] pl[%f] tb[%d]\n", il, hL[il], padlength[il], tbRange[il]);
+                       tracklet.SetNTimeBinsRange(tbRange);
+                       tracklet.SetX0(x0);
+                       tracklet.Init(&t);
+                       if(!tracklet.AttachClustersIter(stackLayer, 1000.)) continue;
+
+                       //if(!tracklet.AttachClusters(stackLayer, kTRUE)) continue;
+                       //if(!tracklet.Fit()) continue;
+               }
+               Int_t ncl = tracklet.GetN();
+               //AliInfo(Form("N clusters %d", ncl));
+               
+               // Discard tracklet if bad quality.
+               //Check if this condition is not already checked during building of the tracklet
+    if(ncl < fTimeBinsPerPlane * fRecoParam->GetFindableClusters()){
+                       //AliInfo(Form("Discard tracklet for %d nclusters", ncl));
+                       continue;
+               }
+               
+               // load tracklet to the tracker and the track
+               Int_t index = SetTracklet(&tracklet);
+               t.SetTracklet(&tracklet, iplane, index);
+               
+               // Calculate the mean material budget along the path inside the chamber
+    Double_t param[7];
+               AliTracker::MeanMaterialBudget(xyz0, xyz1, param);      
+    // The mean propagation parameters
+    Double_t xrho = param[0]*param[4]; // density*length
+    Double_t xx0  = param[1]; // radiation length
+               
+               // Propagate and update track
+               t.PropagateTo(tracklet.GetX0(), xx0, xrho);
+         if (!AdjustSector(&t)) break;
+               Double_t maxChi2 = t.GetPredictedChi2(&tracklet);
+               if (maxChi2<1e+10 && t.Update(&tracklet, maxChi2)){ 
+                       nClustersExpected += ncl;
+               }
+               // Reset material budget if 2 consecutive gold
+               if(iplane>0 && ncl + 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 = ncl / Float_t(fTimeBinsPerPlane);
+               //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);
+               //printf("ratio1     %f             [>   0.6]\n", ratio1); 
+               //printf("ratio0+ratio1 %f          [>   1.5]\n", ratio0+ratio1); 
+               //printf("t.GetNCross()  %d         [==    0]\n", t.GetNCross()); 
+               //printf("TMath::Abs(t.GetSnp()) %f [<  0.85]\n", TMath::Abs(t.GetSnp()));
+               //printf("t.GetNumberOfClusters() %d [>    20]\n", t.GetNumberOfClusters());
+    
+               if (//(tracklet.GetChi2()      <  18.0) && TO DO check with FindClusters and move it to AliTRDseed::Update 
+        (ratio0                  >   0.8) && 
+        //(ratio1                  >   0.6) && 
+        //(ratio0+ratio1           >   1.5) && 
+        (t.GetNCross()           ==    0) && 
+        (TMath::Abs(t.GetSnp())  <  0.85) &&
+        (t.GetNumberOfClusters() >    20)) t.MakeBackupTrack();
+               
+       } // end planes loop
+
+#ifdef DEBUG
+       if(AliTRDReconstructor::StreamLevel() > 1){
+               TTreeSRedirector &cstreamer = *fDebugStreamer;
+               cstreamer << "FollowBackProlongation"
+                       << "ncl="      << nClustersExpected
+                       << "track.="   << &t
+                       << "\n";
+       }
+#endif
+       
+       return nClustersExpected;
+}
+
+//____________________________________________________________________
+void AliTRDtrackerV1::UnloadClusters() 
+{ 
+  //
+  // Clears the arrays of clusters and tracks. Resets sectors and timebins 
+  //
+
+  Int_t i;
+  Int_t nentr;
+
+  nentr = fClusters->GetEntriesFast();
+       //AliInfo(Form("clearing %d clusters", nentr));
+  for (i = 0; i < nentr; i++) {
+    delete fClusters->RemoveAt(i);
+  }
+  fNclusters = 0;
+
+  nentr = fTracklets->GetEntriesFast();
+       //AliInfo(Form("clearing %d tracklets", nentr));
+  for (i = 0; i < nentr; i++) {
+    delete fTracklets->RemoveAt(i);
+  }
+
+  nentr = fSeeds->GetEntriesFast();
+       //AliInfo(Form("clearing %d seeds", nentr));
+  for (i = 0; i < nentr; i++) {
+    delete fSeeds->RemoveAt(i);
+  }
+
+  nentr = fTracks->GetEntriesFast();
+  //AliInfo(Form("clearing %d tracks", nentr));
+       for (i = 0; i < nentr; i++) {
+    delete fTracks->RemoveAt(i);
+  }
+
+  Int_t nsec = AliTRDgeometry::kNsect;
+  for (i = 0; i < nsec; i++) {    
+    for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
+      fTrSec[i]->GetLayer(pl)->Clear();
+    }
+  }
+
+}
+
+//____________________________________________________________________
+AliTRDseedV1* AliTRDtrackerV1::GetTracklet(AliTRDtrackV1 *track, Int_t p, Int_t &idx)
+{
+// Find tracklet for TRD track <track>
+// Parameters
+// - track
+// - sector
+// - plane
+// - index
+// Output
+// tracklet
+// index
+// Detailed description
+//
+       idx = track->GetTrackletIndex(p);
+       //AliInfo(Form("looking for tracklet in plane %d idx %d [%d]", p, idx, track->GetTrackletIndex(p)));
+       AliTRDseedV1 *tracklet = idx<0 ? 0x0 : (AliTRDseedV1*)fTracklets->UncheckedAt(idx);
+       //AliInfo(Form("found 0x%x @ %d", tracklet, idx));
+
+//   Int_t *index = track->GetTrackletIndexes();
+//   for (UInt_t i = 0; i < 6; i++) AliInfo(Form("index[%d] = %d", i, index[i]));
+// 
+//     for (UInt_t i = 0; i < 6/*kMaxTimeBinIndex*/; i++) {
+//     if (index[i] < 0) continue;
+// 
+//             tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index[i]);
+//             if(!tracklet) break;
+// 
+//             if(tracklet->GetPlane() != p) continue;
+// 
+//             idx = index[i];
+//     }
+
+       return tracklet;
+}
+
+//____________________________________________________________________
+Int_t AliTRDtrackerV1::SetTracklet(AliTRDseedV1 */*tracklet*/)
+{
+// Add this tracklet to the list of tracklets stored in the tracker
+//
+// Parameters
+//   - tracklet : pointer to the tracklet to be added to the list
+//
+// Output
+//   - the index of the new tracklet in the tracker tracklets list
+//
+// Detailed description
+// Build the tracklets list if it is not yet created (late initialization)
+// and adds the new tracklet to the list.
+//
+       if(!fTracklets){
+               fTracklets = new TClonesArray("AliTRDseedV1", AliTRDgeometry::Nsect()*kMaxTracksStack);
+               fTracklets->SetOwner(kTRUE);
+       }
+       Int_t nentries = fTracklets->GetEntriesFast();
+       //AliTRDseedV1 *t = new ((*fTracklets)[nentries]) AliTRDseedV1(*tracklet);
+       //AliInfo(Form("0x%x @ %d", t, nentries));
+       return nentries;
+}
+
 //____________________________________________________________________
 Int_t AliTRDtrackerV1::Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *sector
                                        , AliESDEvent *esd)
@@ -203,10 +931,10 @@ Int_t AliTRDtrackerV1::Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *se
                                              + 2 * pp->GetLengthOPad() + 2 * pp->GetLengthRim();
                        //Debug
                        Double_t z0  = fGeom->GetRow0((Int_t)(ilayer/nTimeBins),istack,0);
-                       const AliTRDpropagationLayer kSMlayer(*(sector->GetLayer(ilayer)));
-                       stackLayer[ilayer] = kSMlayer;
+                       const AliTRDpropagationLayer ksmLayer(*(sector->GetLayer(ilayer)));
+                       stackLayer[ilayer] = ksmLayer;
 #ifdef DEBUG
-                       stackLayer[ilayer].SetDebugStream(fDebugStreamerV1);
+                       stackLayer[ilayer].SetDebugStream(fDebugStreamer);
 #endif                 
                        stackLayer[ilayer].SetRange(z0 - stacklength, stacklength);
                        stackLayer[ilayer].SetSector(sector->GetSector());
@@ -258,10 +986,8 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
   // 8. Build ESD track and register it to the output list
   //
 
-       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-       Int_t nTimeBins = cal->GetNumberOfTimeBins();
        AliTRDseedV1 sseed[kMaxTracksStack*6]; // to be initialized
-       Int_t pars[3]; // MakeSeeds parameters
+       Int_t pars[4]; // MakeSeeds parameters
 
        //Double_t alpha = AliTRDgeometry::GetAlpha();
        //Double_t shift = .5 * alpha;
@@ -291,8 +1017,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                        if(ntracks == kMaxTracksStack) break;
                }
 #ifdef DEBUG
-               if(AliTRDReconstructor::StreamLevel() > 1) 
-                  AliInfo(Form("Candidate TRD tracks %d in stack %d.", ntracks, pars[1]));
+               if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Candidate TRD tracks %d in stack %d iteration %d.", ntracks, pars[1], fSieveSeeding));
 #endif         
                if(!ntracks) break;
                
@@ -344,7 +1069,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                                        nlayers++;
        
                                        // Cooking label
-                                       for (Int_t itime = 0; itime < nTimeBins; itime++) {
+                                       for (Int_t itime = 0; itime < fTimeBinsPerPlane; itime++) {
                                                if(!sseed[jseed].IsUsable(itime)) continue;
                                                naccepted++;
                                                Int_t tindex = 0, ilab = 0;
@@ -356,12 +1081,12 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                                }
                                // Filter duplicated tracks
                                if (nused > 30){
-                                       //printf("Skip nused %d\n", nused);
+                                       printf("Skip %d nused %d\n", trackIndex, nused);
                                        fakeTrack[trackIndex] = kTRUE;
                                        continue;
                                }
                                if (Float_t(nused)/ncl >= .25){
-                                       //printf("Skip nused/ncl >= .25\n");
+                                       printf("Skip %d nused/ncl >= .25\n", trackIndex);
                                        fakeTrack[trackIndex] = kTRUE;
                                        continue;
                                }
@@ -390,12 +1115,12 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
        
                                case 4:
                                        if (nlayers == 3){skip = kTRUE; break;}
-                                       if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) - nused/(nlayers-3.0) < -15.0){skip = kTRUE; break;}
+                                       //if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) - nused/(nlayers-3.0) < -15.0){skip = kTRUE; break;}
                                        break;
                                }
                                if(skip){
                                        candidates++;
-                                       //printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused);
+                                       printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused);
                                        continue;
                                }
                                signedTrack[trackIndex] = kTRUE;
@@ -463,13 +1188,14 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                                        Int_t nclusters = 0;
                                        AliTRDseedV1 *dseed[6];
                                        for(int is=0; is<6; is++){
-                                               dseed[is] = new AliTRDseedV1(sseed[trackIndex*6+is], kTRUE);
+                                               dseed[is] = new AliTRDseedV1(sseed[trackIndex*6+is]);
+                                               dseed[is]->SetOwner();
                                                nclusters += sseed[is].GetN2();
                                                //for(int ic=0; ic<30; ic++) if(sseed[trackIndex*6+is].GetClusters(ic)) printf("l[%d] tb[%d] cptr[%p]\n", is, ic, sseed[trackIndex*6+is].GetClusters(ic));
                                        }
                                        //Int_t eventNrInFile = esd->GetEventNumberInFile();
                                        //AliInfo(Form("Number of clusters %d.", nclusters));
-                                       TTreeSRedirector &cstreamer = *fDebugStreamerV1;
+                                       TTreeSRedirector &cstreamer = *fDebugStreamer;
                                        cstreamer << "Clusters2TracksStack"
                                                << "Iter="      << fSieveSeeding
                                                << "Like="      << fTrackQuality[trackIndex]
@@ -502,11 +1228,12 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                                }
 #endif
                        
-                               AliTRDtrack *track = AliTRDtrackerV1::RegisterSeed(&sseed[trackIndex*kNPlanes], trackParams);
+                               AliTRDtrackV1 *track = AliTRDtrackerV1::MakeTrack(&sseed[trackIndex*kNPlanes], trackParams);
                                if(!track){
                                        AliWarning("Fail to build a TRD Track.");
                                        continue;
                                }
+                               AliInfo("End of MakeTrack()");
                                AliESDtrack esdTrack;
                                esdTrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
                                esdTrack.SetLabel(track->GetLabel());
@@ -526,7 +1253,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                quality = BuildSeedingConfigs(layer, configs);
                //if(quality < fRecoParam->GetPlaneQualityThreshold()) break;
                
-               for(Int_t il = 0; il < kNPlanes * nTimeBins; il++) layer[il].BuildIndices(fSieveSeeding);
+               for(Int_t il = 0; il < kNPlanes * fTimeBinsPerPlane; il++) layer[il].BuildIndices(fSieveSeeding);
 
 #ifdef DEBUG
                                if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Sieve level %d Plane config %d %d %d Quality %f", fSieveSeeding, configs[0], configs[1], configs[2], quality));
@@ -568,12 +1295,9 @@ Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDstackLayer *layers
   // The overall chamber quality is given by the product of this 2 contributions.
   // 
 
-       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-       Int_t nTimeBins = cal->GetNumberOfTimeBins();
-
        Double_t chamberQA[kNPlanes];
        for(int iplane=0; iplane<kNPlanes; iplane++){
-               chamberQA[iplane] = CookPlaneQA(&layers[iplane*nTimeBins]);
+               chamberQA[iplane] = MakeSeedingPlanes(&layers[iplane*fTimeBinsPerPlane]);
                //printf("chamberQA[%d] = %f\n", iplane, chamberQA[iplane]);
        }
 
@@ -643,8 +1367,6 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
   // 15. Register seeds.
   //
 
-       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-       Int_t nTimeBins = cal->GetNumberOfTimeBins();
        AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0}; // initilize seeding clusters
        AliTRDseedV1 *cseed = &sseed[0]; // initialize tracklets for first track
        Int_t ncl, mcl; // working variable for looping over clusters
@@ -668,19 +1390,23 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
 #endif
        
        // Init chambers geometry
+       Int_t det, tbRange[6]; // time bins inside the detector geometry
        Double_t hL[kNPlanes];       // Tilting angle
        Float_t padlength[kNPlanes]; // pad lenghts
        AliTRDpadPlane *pp;
        for(int il=0; il<kNPlanes; il++){
                pp = fGeom->GetPadPlane(il, istack); // istack has to be imported
                hL[il]        = TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle());
-               padlength[il] = 10.; //pp->GetLengthIPad();
+               padlength[il] = pp->GetLengthIPad();
+               det           = il; // to be fixed !!!!!
+               tbRange[il]   = fTimeBinsPerPlane; //Int_t(AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght() * AliTRDCommonParam::Instance()->GetSamplingFrequency()/AliTRDcalibDB::Instance()->GetVdriftDet()->GetValue(det));
+               //printf("%d hl[%f] pl[%f] tb[%d]\n", il, hL[il], padlength[il], tbRange[il]);
        }
 
        Double_t cond0[4], cond1[4], cond2[4];
        // make seeding layers (to be moved in Clusters2TracksStack)
        AliTRDstackLayer *layer[] = {0x0, 0x0, 0x0, 0x0};
-       for(int isl=0; isl<kNSeedPlanes; isl++) layer[isl] = MakeSeedingLayer(&layers[planes[isl] * nTimeBins], planes[isl]);
+       for(int isl=0; isl<kNSeedPlanes; isl++) layer[isl] = MakeSeedingLayer(&layers[planes[isl] * fTimeBinsPerPlane], planes[isl]);
 
 
        // Start finding seeds
@@ -718,14 +1444,18 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                chi2[0] = 0.; chi2[1] = 0.;
                                AliTRDseedV1 *tseed = 0x0;
                                for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
-                                       tseed = &cseed[planes[iLayer]];
+                                       Int_t jLayer = planes[iLayer];
+                                       tseed = &cseed[jLayer];
                                        tseed->SetRecoParam(fRecoParam);
-                                       tseed->SetLayer(planes[iLayer]);
-                                       tseed->SetTilt(hL[planes[iLayer]]);
-                                       tseed->SetPadLength(TMath::Sqrt(c[iLayer]->GetSigmaZ2()*12));
-                                       tseed->SetX0(layer[iLayer]->GetX());
-
-                                       tseed->Update(fFitter->GetRiemanFitter());
+                                       tseed->SetPlane(jLayer);
+                                       tseed->SetTilt(hL[jLayer]);
+                                       tseed->SetPadLength(padlength[jLayer]);
+                                       tseed->SetNTimeBinsRange(tbRange[jLayer]);
+                                       tseed->SetX0(layer[iLayer]->GetX());//layers[jLayer*fTimeBinsPerPlane].GetX());
+
+                                       tseed->Init(fFitter->GetRiemanFitter());
+                                       // temporary until new AttachClusters()
+                                       tseed->SetX0(layers[(jLayer+1)*fTimeBinsPerPlane-1].GetX());
                                        chi2[0] += tseed->GetChi2Z(c[iLayer]->GetZ());
                                        chi2[1] += tseed->GetChi2Y(c[iLayer]->GetY());
                                }
@@ -744,7 +1474,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                        }
                                        Float_t threshold = .5;//1./(3. - sLayer);
                                        Int_t ll = c[3]->GetLabel(0);
-                                       TTreeSRedirector &cs0 = *fDebugStreamerV1;
+                                       TTreeSRedirector &cs0 = *fDebugStreamer;
                                                        cs0 << "MakeSeeds0"
                                                        <<"isFake=" << isFake
                                                        <<"label=" << ll
@@ -783,7 +1513,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                        }
                                        Double_t xpos[4];
                                        for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX();
-                                       TTreeSRedirector &cstreamer = *fDebugStreamerV1;
+                                       TTreeSRedirector &cstreamer = *fDebugStreamer;
                                                        cstreamer << "MakeSeeds1"
                                                << "isFake=" << isFake
                                                << "config="   << config
@@ -812,10 +1542,9 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                Int_t nUsedCl = 0;
                                Int_t nlayers = 0;
                                for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
-                                       AliTRDseedV1 tseed = cseed[planes[iLayer]];
-                                       if(!tseed.AttachClustersIter(&layers[planes[iLayer]*nTimeBins], 5., kFALSE, c[iLayer])) continue;
-                                       cseed[planes[iLayer]] = tseed;
-                                       nUsedCl += cseed[planes[iLayer]].GetNUsed();
+                                       Int_t jLayer = planes[iLayer];
+                                       if(!cseed[jLayer].AttachClustersIter(&layers[jLayer*fTimeBinsPerPlane], 5., kFALSE, c[iLayer])) continue;
+                                       nUsedCl += cseed[jLayer].GetNUsed();
                                        if(nUsedCl > 25) break;
                                        nlayers++;
                                }
@@ -828,7 +1557,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                fFitter->FitRieman(&cseed[0], &planes[0]);
                                AliRieman *rim = fFitter->GetRiemanFitter();
                                for(int iLayer=0; iLayer<4; iLayer++){
-                                       cseed[planes[iLayer]].Update(rim);
+                                       cseed[planes[iLayer]].Init(rim);
                                        chi2[0] += (Float_t)cseed[planes[iLayer]].GetChi2Z();
                                        chi2[1] += cseed[planes[iLayer]].GetChi2Y();
                                }
@@ -855,24 +1584,21 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                        // prepare extrapolated seed
                                        cseed[jLayer].Reset();
                                        cseed[jLayer].SetRecoParam(fRecoParam);
-                                       cseed[jLayer].SetLayer(jLayer);
+                                       cseed[jLayer].SetPlane(jLayer);
                                        cseed[jLayer].SetTilt(hL[jLayer]);
-                                       cseed[jLayer].SetX0(layers[jLayer * nTimeBins + (nTimeBins/2)].GetX());  // ????????
-                                       //cseed[jLayer].SetPadLength(??????????);
-                                       cseed[jLayer].Update(rim);
-                                       
-                                       AliTRDcluster *cd = FindSeedingCluster(&layers[jLayer*nTimeBins], &cseed[jLayer]);
-                                       if(cd == 0x0) continue;
-//                                     if(!cd) continue;
-                                       cseed[jLayer].SetPadLength(TMath::Sqrt(cd->GetSigmaZ2() * 12.));
-                                       cseed[jLayer].SetX0(cd->GetX());        // reference defined by a seedingLayer which is defined by the x-coordinate of the layers inside
+                                       cseed[jLayer].SetX0(layers[(jLayer +1) * fTimeBinsPerPlane-1].GetX());
+                                       cseed[jLayer].SetPadLength(padlength[jLayer]);
+                                       cseed[jLayer].SetNTimeBinsRange(tbRange[jLayer]);
+                                       cseed[jLayer].Init(rim);
+//                                     AliTRDcluster *cd = FindSeedingCluster(&layers[jLayer*fTimeBinsPerPlane], &cseed[jLayer]);
+//                                     if(cd == 0x0) continue;
 
                                        // fit extrapolated seed
                                        AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
                                        if ((jLayer == 0) && !(cseed[1].IsOK())) continue;
                                        if ((jLayer == 5) && !(cseed[4].IsOK())) continue;
                                        AliTRDseedV1 tseed = cseed[jLayer];
-                                       if(!tseed.AttachClustersIter(&layers[jLayer*nTimeBins], 1000.)) continue;
+                                       if(!tseed.AttachClustersIter(&layers[jLayer*fTimeBinsPerPlane], 1000.)) continue;
                                        cseed[jLayer] = tseed;
                                        nusedf += cseed[jLayer].GetNUsed(); // debug value
                                        AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
@@ -900,7 +1626,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                fFitter->FitRieman(&cseed[0]);
                                Double_t chi2ZF = 0., chi2RF = 0.;
                                for(int ilayer=0; ilayer<6; ilayer++){
-                                       cseed[ilayer].Update(fFitter->GetRiemanFitter());
+                                       cseed[ilayer].Init(fFitter->GetRiemanFitter());
                                        if (!cseed[ilayer].IsOK()) continue;
                                        //tchi2 = cseed[ilayer].GetChi2Z();
                                        //printf("layer %d chi2 %e\n", ilayer, tchi2);
@@ -913,7 +1639,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                // do the final track fitting
                                fFitter->SetLayers(nlayers);
 #ifdef DEBUG
-                               fFitter->SetDebugStream(fDebugStreamerV1);
+                               fFitter->SetDebugStream(fDebugStreamer);
 #endif
                                fTrackQuality[ntracks] = fFitter->FitHyperplane(&cseed[0], chi2ZF, GetZ());
                                Double_t param[3];
@@ -953,7 +1679,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
 #ifdef DEBUG
                                if(AliTRDReconstructor::StreamLevel() >= 2){
                                        Double_t curv = (fFitter->GetRiemanFitter())->GetC();
-                                       TTreeSRedirector &cstreamer = *fDebugStreamerV1;
+                                       TTreeSRedirector &cstreamer = *fDebugStreamer;
                                        cstreamer << "MakeSeeds2"
                                                << "C="       << curv
                                                << "Chi2R="   << chi2r
@@ -995,7 +1721,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
 }
 
 //_____________________________________________________________________________
-AliTRDtrack *AliTRDtrackerV1::RegisterSeed(AliTRDseedV1 *seeds, Double_t *params)
+AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params)
 {
   //
   // Build a TRD track out of tracklet candidates
@@ -1012,9 +1738,6 @@ AliTRDtrack *AliTRDtrackerV1::RegisterSeed(AliTRDseedV1 *seeds, Double_t *params
   // To be discussed with Marian !!
   //
 
-       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-       Int_t nTimeBins = cal->GetNumberOfTimeBins();
-       
   Double_t alpha = AliTRDgeometry::GetAlpha();
   Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
   Double_t c[15];
@@ -1025,47 +1748,29 @@ AliTRDtrack *AliTRDtrackerV1::RegisterSeed(AliTRDseedV1 *seeds, Double_t *params
   c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0;  c[ 9] = 0.1;
   c[10] = 0.0; c[11] = 0.0; c[12] = 0.0;  c[13] = 0.0; c[14] = params[5]*params[5]*0.01;
 
-  Int_t index = 0;
-  AliTRDcluster *cl = 0;
-
-  for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
-    if (seeds[ilayer].IsOK()) {
-      for (Int_t itime = nTimeBins - 1; itime > 0; itime--) {
-       if (seeds[ilayer].GetIndexes(itime) > 0) {
-         index = seeds[ilayer].GetIndexes(itime);
-         cl    = seeds[ilayer].GetClusters(itime);
-         break;
-       }
-      }
-    }
-    if (index > 0) {
-      break;
-    }
-  }
-  if (cl == 0) return 0;
-  AliTRDtrack *track = new AliTRDtrack(cl
-                                      ,index
-                                      ,&params[1]
-                                      ,c
-                                      ,params[0]
-                                      ,params[6]*alpha+shift);
-       // SetCluster(cl, 0); // A. Bercuci
+  AliTRDtrackV1 *track = new AliTRDtrackV1(seeds, &params[1], c, params[0], params[6]*alpha+shift);
        track->PropagateTo(params[0]-5.0);
   track->ResetCovariance(1);
-  Int_t rc = FollowBackProlongation(*track);
-  if (rc < 30) {
+  Int_t nc = FollowBackProlongation(*track);
+       AliInfo(Form("N clusters for track %d", nc));
+       if (nc < 30) {
     delete track;
-    track = 0;
-  }
-  else {
-    track->CookdEdx();
-    track->CookdEdxTimBin(-1);
-    CookLabel(track,0.9);
+    track = 0x0;
+  } else {
+//     track->CookdEdx();
+//     track->CookdEdxTimBin(-1);
+//     CookLabel(track, 0.9);
   }
 
   return track;
 }
 
+//____________________________________________________________________
+void AliTRDtrackerV1::CookLabel(AliKalmanTrack */*pt*/, Float_t /*wrong*/) const
+{
+       // to be implemented, preferably at the level of TRD tracklet. !!!!!!!
+}
+
 //____________________________________________________________________
 void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers
                                        , AliTRDseedV1 *cseed)
@@ -1127,7 +1832,7 @@ void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers
 }
 
 //____________________________________________________________________
-Double_t  AliTRDtrackerV1::CookPlaneQA(AliTRDstackLayer *layers)
+Double_t  AliTRDtrackerV1::MakeSeedingPlanes(AliTRDstackLayer *layers)
 {
   //
   // Calculate plane quality for seeding.
@@ -1181,7 +1886,7 @@ Double_t  AliTRDtrackerV1::CookPlaneQA(AliTRDstackLayer *layers)
 //     fitter.PrintResults(3);
 //     Double_t a = fitter.GetParameter(1);
 // 
-//     printf("nclDev(%f)  a(%f)\n", nclDev, a);
+//     printf("ncl_dev(%f)  a(%f)\n", ncl_dev, a);
 //     return quality*TMath::Exp(-a);
 }
 
@@ -1239,7 +1944,7 @@ Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed
 #ifdef DEBUG
        //AliInfo(Form("sumda(%f) chi2[0](%f) chi2[1](%f) likea(%f) likechi2y(%f) likechi2z(%f) nclusters(%d) likeN(%f)", sumda, chi2[0], chi2[1], likea, likechi2y, likechi2z, nclusters, likeN));
        if(AliTRDReconstructor::StreamLevel() >= 2){
-               TTreeSRedirector &cstreamer = *fDebugStreamerV1;
+               TTreeSRedirector &cstreamer = *fDebugStreamer;
                cstreamer << "CookLikelihood"
                        << "sumda="     << sumda
                        << "chi0="      << chi2[0]
@@ -1623,11 +2328,11 @@ AliTRDstackLayer *AliTRDtrackerV1::MakeSeedingLayer(AliTRDstackLayer *layers
        
 #ifdef DEBUG
        if(AliTRDReconstructor::StreamLevel() >= 3){
-               TMatrixT<double> hist(nRows, nCols);
+               TMatrixD hist(nRows, nCols);
                for(Int_t i = 0; i < nRows; i++)
                        for(Int_t j = 0; j < nCols; j++)
                                hist(i,j) = histogram[i][j];
-               TTreeSRedirector &cstreamer = *fDebugStreamerV1;
+               TTreeSRedirector &cstreamer = *fDebugStreamer;
                cstreamer << "MakeSeedingLayer"
                        << "Iteration="  << fSieveSeeding
                        << "plane="      << plane
@@ -1644,8 +2349,7 @@ AliTRDstackLayer *AliTRDtrackerV1::MakeSeedingLayer(AliTRDstackLayer *layers
 }
 
 //____________________________________________________________________
-void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig
-                                     , Int_t planes[4]) const
+void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig, Int_t planes[4])
 {
   //
   // Map seeding configurations to detector planes.
@@ -1786,8 +2490,7 @@ void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig
 }
 
 //____________________________________________________________________
-void AliTRDtrackerV1::GetExtrapolationConfig(Int_t iconfig
-                                           , Int_t planes[2]) const
+void AliTRDtrackerV1::GetExtrapolationConfig(Int_t iconfig, Int_t planes[2])
 {
   //
   // Returns the extrapolation planes for a seeding configuration.
index 9c69df97d8e67c34751227a4fdce92bdbb49de5b..e28b9ea53111769d43f6d7d3e91988676375466f 100644 (file)
@@ -36,7 +36,8 @@ class AliTRDseedV1;
 class AliTRDstackLayer;
 class AliTRDtrackerFitter;
 class AliTRDrecoParam;
-
+class AliTRDtrackV1;
+class AliTrackPoint;
 class AliTRDtrackerV1 : public AliTRDtracker
 {
 
@@ -50,42 +51,50 @@ class AliTRDtrackerV1 : public AliTRDtracker
        };
        AliTRDtrackerV1(AliTRDrecoParam *p = 0x0);
        AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p);
-       ~AliTRDtrackerV1();
+       virtual ~AliTRDtrackerV1();
   
        Int_t          Clusters2Tracks(AliESDEvent *esd);
-       void           GetSeedingConfig(Int_t iconfig, Int_t planes[4]) const;
-       void           GetExtrapolationConfig(Int_t iconfig, Int_t planes[2]) const;
+       static void    GetExtrapolationConfig(Int_t iconfig, Int_t planes[2]);
+       static void    GetSeedingConfig(Int_t iconfig, Int_t planes[4]);
+       Int_t          FollowBackProlongation(AliTRDtrackV1 &t);
+       Int_t          FollowProlongation(AliTRDtrackV1 &t);
+       Int_t          PropagateBack(AliESDEvent *event);
+       Int_t          RefitInward(AliESDEvent *event);
        void           SetRecoParam(AliTRDrecoParam *p){fRecoParam = p;}
-       
- protected:
+       void           UnloadClusters();
 
+ protected:
        Double_t       BuildSeedingConfigs(AliTRDstackLayer *layer, Int_t *configs);
-       Int_t          Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *sector, AliESDEvent *esd);
+       Int_t          Clusters2TracksSM( AliTRDtracker::AliTRDtrackingSector *sector, AliESDEvent *esd);
        Int_t          Clusters2TracksStack(AliTRDstackLayer *layer, TClonesArray *esdTrackList);
-       Double_t       CookPlaneQA(AliTRDstackLayer *layer);
-       Double_t       CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4], Double_t *chi2);
+       void           CookLabel(AliKalmanTrack *pt, Float_t wrong) const;
        Int_t          GetSeedingLayers(AliTRDstackLayer *layers, Double_t *params);
        void           GetMeanCLStack(AliTRDstackLayer *layers, Int_t *planes, Double_t *params);
+       AliTRDseedV1*  GetTracklet(AliTRDtrackV1 *trk, Int_t plane, Int_t &idx);
+       virtual Bool_t GetTrackPoint(Int_t index, AliTrackPoint &p) const;
        AliTRDcluster *FindSeedingCluster(AliTRDstackLayer *layers, AliTRDseedV1/*AliRieman*/ *sfit);
-       void           ImproveSeedQuality(AliTRDstackLayer *layer, AliTRDseedV1 *cseed);
-       Int_t          MakeSeeds(AliTRDstackLayer *layers, AliTRDseedV1 *sseed, Int_t *ipar);
+       
+       Double_t       MakeSeedingPlanes(AliTRDstackLayer *layer);
        AliTRDstackLayer *MakeSeedingLayer(AliTRDstackLayer *layers, Int_t Plane);
-       AliTRDtrack*   RegisterSeed(AliTRDseedV1 *seeds, Double_t *params);
-
- private:
+       Int_t          MakeSeeds(AliTRDstackLayer *layers, AliTRDseedV1 *sseed, Int_t *ipar);
+       AliTRDtrackV1*   MakeTrack(AliTRDseedV1 *seeds, Double_t *params);
+       Int_t          SetTracklet(AliTRDseedV1 *tracklet);
 
+private:
        AliTRDtrackerV1(const AliTRDtrackerV1 &tracker);
        AliTRDtrackerV1 &operator=(const AliTRDtrackerV1 &tracker);
+       Double_t       CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4], Double_t *chi2);
+       void           ImproveSeedQuality(AliTRDstackLayer *layer, AliTRDseedV1 *cseed);
 
- private:
+private:
 
-       static Double_t      fgTopologicQA[kNConfigs];        //  Topologic quality
+       static Double_t      fgTopologicQA[kNConfigs];         //  Topologic quality
        Double_t             fTrackQuality[kMaxTracksStack];  //  Track quality 
        Int_t                fSeedLayer[kMaxTracksStack];     //  Seed layer
        Int_t                fSieveSeeding;                   //! Seeding iterator
+       TClonesArray        *fTracklets;                      // List of tracklets for all sectors
        AliTRDrecoParam     *fRecoParam;                      //  Reconstruction parameters
        AliTRDtrackerFitter *fFitter;                         //! Fitter class of the tracker
-       TTreeSRedirector    *fDebugStreamerV1;                //! Debug stream of the tracker
 
        ClassDef(AliTRDtrackerV1, 1)                          //  Stand alone tracker development class
 
index 29dab8abab65ecea341722549f960eb0b8ce490a..19120f9af497a19cceb231cab8defecd1d2dfcf9 100644 (file)
 
 #pragma link C++ class  AliTRDcluster+;
 #pragma link C++ class  AliTRDclusterMI+;
-
 #pragma link C++ class  AliTRDclusterizer+;
-
 #pragma link C++ class  AliTRDclusterCorrection+;
-
 #pragma link C++ class  AliTRDtransform+;
 
 #pragma link C++ class  AliTRDtrack+;
 #pragma link C++ class  AliTRDReconstructor+;
 #pragma link C++ class  AliTRDRecParam+;
 
-#pragma link C++ class  AliTRDtrackingAnalysis+;
-
-#pragma link C++ class AliTRDrecoParam+;
-#pragma link C++ class AliTRDseedV1+;
-#pragma link C++ class AliTRDtrackerV1+;
-#pragma link C++ class AliTRDstackLayer+;
-#pragma link C++ class AliTRDtrackerFitter+;
-#pragma link C++ class AliTRDQADataMakerRec+;
+#pragma link C++ class  AliTRDrecoParam+;
+#pragma link C++ class  AliTRDseedV1+;
+#pragma link C++ class  AliTRDtrackV1+;
+#pragma link C++ class  AliTRDtrackerV1+;
+#pragma link C++ class  AliTRDstackLayer+;
+#pragma link C++ class  AliTRDtrackerFitter+;
 
+#pragma link C++ class  AliTRDQADataMakerRec+;
+#pragma link C++ class  AliTRDtrackingAnalysis+;
 
 #endif
index 5830ad77f9efbe33982bb1d5ac38928ed85c75af..09f546da1bf3d0fbdfbf5f51d8035b1572f27573 100644 (file)
@@ -13,6 +13,7 @@ SRCS= AliTRDcluster.cxx \
       AliTRDtrackingAnalysis.cxx \
       AliTRDrecoParam.cxx \
       AliTRDseedV1.cxx \
+      AliTRDtrackV1.cxx \
       AliTRDtrackerV1.cxx \
       AliTRDpropagationLayer.cxx \
       AliTRDstackLayer.cxx \