New stand alone tracking by Alexadru and Martin
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Nov 2007 09:21:46 +0000 (09:21 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 29 Nov 2007 09:21:46 +0000 (09:21 +0000)
21 files changed:
TRD/AliTRDReconstructor.cxx
TRD/AliTRDdigitizer.cxx
TRD/AliTRDpropagationLayer.cxx [new file with mode: 0644]
TRD/AliTRDpropagationLayer.h [new file with mode: 0644]
TRD/AliTRDrecoParam.cxx [new file with mode: 0644]
TRD/AliTRDrecoParam.h [new file with mode: 0644]
TRD/AliTRDseed.cxx
TRD/AliTRDseed.h
TRD/AliTRDseedV1.cxx [new file with mode: 0644]
TRD/AliTRDseedV1.h [new file with mode: 0644]
TRD/AliTRDstackLayer.cxx [new file with mode: 0644]
TRD/AliTRDstackLayer.h [new file with mode: 0644]
TRD/AliTRDtrack.cxx
TRD/AliTRDtracker.cxx
TRD/AliTRDtracker.h
TRD/AliTRDtrackerFitter.cxx [new file with mode: 0644]
TRD/AliTRDtrackerFitter.h [new file with mode: 0644]
TRD/AliTRDtrackerV1.cxx [new file with mode: 0644]
TRD/AliTRDtrackerV1.h [new file with mode: 0644]
TRD/TRDrecLinkDef.h
TRD/libTRDrec.pkg

index df406e5..030725a 100644 (file)
@@ -36,6 +36,8 @@
 #include "AliTRDgtuTrack.h"
 #include "AliTRDrawData.h"
 #include "AliTRDdigitsManager.h"
+#include "AliTRDtrackerV1.h"
+#include "AliTRDrecoParam.h"
 
 ClassImp(AliTRDReconstructor)
 
@@ -90,6 +92,7 @@ void AliTRDReconstructor::Reconstruct(TTree *digitsTree
   //
   // Reconstruct clusters
   //
+
   AliInfo("Reconstruct TRD clusters from Digits [Digit TTree -> Cluster TTree]");
 
   AliTRDclusterizer clusterer("clusterer","TRD clusterizer");
@@ -106,7 +109,9 @@ AliTracker *AliTRDReconstructor::CreateTracker() const
   // Create a TRD tracker
   //
 
-  return new AliTRDtracker(NULL);
+  //return new AliTRDtracker(NULL);
+
+  return new AliTRDtrackerV1(NULL, AliTRDrecoParam::GetLowFluxParam());
 
 }
 
index 564cca3..e9c9e26 100644 (file)
@@ -1291,7 +1291,7 @@ Bool_t AliTRDdigitizer::MakeDigits()
 
             for (iTime = 0; iTime < nTimeTotal; iTime++) {   
               // Store the amplitude of the digit if above threshold
-              // if (outADC[iTime] > (simParam->GetADCbaseline() + simParam->GetADCthreshold())) {
+              //if (outADC[iTime] > (simParam->GetADCbaseline() + simParam->GetADCthreshold())) {
               if (outADC[iTime] != 0 ) {   // Now this is enough because there is ZS in raw simulator
                 nDigits++;
                 digits->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
@@ -1501,7 +1501,7 @@ Bool_t AliTRDdigitizer::ConvertSDigits()
 
           for (iTime = 0; iTime < nTimeTotal; iTime++) {
             // Store the amplitude of the digit if above threshold
-            // if (outADC[iTime] > (adcBaseline + adcThreshold)) {
+            //if (outADC[iTime] > (adcBaseline + adcThreshold)) {
            if (outADC[iTime] != 0) {  // now this is ok because there is ZS in raw simulation
               digitsOut->SetDataUnchecked(iRow,iCol,iTime,((Int_t) outADC[iTime]));
              // Copy the dictionary
diff --git a/TRD/AliTRDpropagationLayer.cxx b/TRD/AliTRDpropagationLayer.cxx
new file mode 100644 (file)
index 0000000..93d088f
--- /dev/null
@@ -0,0 +1,353 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+//  The TRD propagation layer                                             //
+//                                                                        //
+//  Authors:                                                              //
+//    Marian Ivanov <M.Ivanov@gsi.de>                                     //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                     //
+//    Markus Fasel <M.Fasel@gsi.de>                                       //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#include "string.h"
+
+#include "TMath.h"
+
+#include "AliTRDpropagationLayer.h"
+//#include "AliTRDtracker.h"
+#include "AliTRDcluster.h"
+#include "AliTRDgeometry.h"
+
+//_____________________________________________________________________________
+AliTRDpropagationLayer::AliTRDpropagationLayer()
+  :TObject()
+  ,fN(0)
+  ,fSec(0)
+  ,fClusters(NULL)
+  ,fIndex(NULL)
+  ,fX(0.)
+  ,fdX(0.)
+  ,fRho(0.)
+  ,fX0(0.)
+  ,fTimeBinIndex(0)
+  ,fPlane(0)
+  ,fYmax(0)
+  ,fYmaxSensitive(0)
+  ,fHole(kFALSE)
+  ,fHoleZc(0)
+  ,fHoleZmax(0)
+  ,fHoleYc(0)
+  ,fHoleYmax(0)
+  ,fHoleRho(0)
+  ,fHoleX0(0)
+{
+  //
+  // Default constructor
+  //
+
+}
+
+//_____________________________________________________________________________
+AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
+                                             , Double_t radLength, Int_t tbIndex, Int_t plane)
+  :TObject()
+  ,fN(0)
+  ,fSec(0)
+  ,fClusters(NULL)
+  ,fIndex(NULL)
+  ,fX(x)
+  ,fdX(dx)
+  ,fRho(rho)
+  ,fX0(radLength)
+  ,fTimeBinIndex(tbIndex)
+  ,fPlane(plane)
+  ,fYmax(0)
+  ,fYmaxSensitive(0)
+  ,fHole(kFALSE)
+  ,fHoleZc(0)
+  ,fHoleZmax(0)
+  ,fHoleYc(0)
+  ,fHoleYmax(0)
+  ,fHoleRho(0)
+  ,fHoleX0(0)
+{ 
+  //
+  // AliTRDpropagationLayer constructor
+  //
+
+  for (Int_t i = 0; i < (Int_t)kZones; i++) {
+    fZc[i]   = 0; 
+    fZmax[i] = 0;
+  }
+
+  if (fTimeBinIndex >= 0) { 
+    fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
+    fIndex    = new UInt_t[kMaxClusterPerTimeBin];
+  }
+
+  for (Int_t i = 0; i < 5; i++) {
+    fIsHole[i] = kFALSE;
+  }
+
+}
+
+//_____________________________________________________________________________
+AliTRDpropagationLayer::AliTRDpropagationLayer(const AliTRDpropagationLayer &p)
+  :TObject((TObject&)p)
+  ,fN(p.fN)
+  ,fSec(p.fSec)
+  ,fClusters(0x0)
+  ,fIndex(0x0)
+  ,fX(p.fX)
+  ,fdX(p.fdX)
+  ,fRho(p.fRho)
+  ,fX0(p.fX0)
+  ,fTimeBinIndex(p.fTimeBinIndex)
+  ,fPlane(p.fPlane)
+  ,fYmax(p.fYmax)
+  ,fYmaxSensitive(p.fYmaxSensitive)
+  ,fHole(p.fHole)
+  ,fHoleZc(p.fHoleZc)
+  ,fHoleZmax(p.fHoleZmax)
+  ,fHoleYc(p.fHoleYc)
+  ,fHoleYmax(p.fHoleYmax)
+  ,fHoleRho(p.fHoleRho)
+  ,fHoleX0(p.fHoleX0)
+{
+  //
+  // AliTRDpropagationLayer copy constructor
+  //
+
+  for (Int_t i = 0; i < (Int_t)kZones; i++) {
+    fZc[i]   = p.fZc[i]; 
+    fZmax[i] = p.fZmax[i];
+               fIsHole[i] = p.fIsHole[i];
+               fZmaxSensitive[i] = p.fZmaxSensitive[i];  
+       }
+
+       // Make a deep copy of the Clusters array and the Index array unless they are needed in class AliTRDstackLayer
+       Int_t arrsize = (Int_t)kMaxClusterPerTimeBin;
+        if (fTimeBinIndex >= 0) { 
+    fClusters = new AliTRDcluster*[arrsize];
+    fIndex    = new UInt_t[arrsize];
+  }
+       memset(fIndex, 0, sizeof(UInt_t)*arrsize);
+       memset(fClusters, 0, sizeof(AliTRDcluster *)*arrsize);
+       for(Int_t i = 0; i < arrsize; i++){
+               fClusters[i] = p.fClusters[i];
+               fIndex[i] = p.fIndex[i];
+       }
+}
+//_____________________________________________________________________________
+AliTRDpropagationLayer::~AliTRDpropagationLayer()
+{
+  //
+  // Destructor
+  //
+
+  if (fTimeBinIndex >= 0) { 
+    delete[] fClusters;
+    delete[] fIndex;
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDpropagationLayer::Copy(TObject &o) const 
+{
+  //
+  // Copy function
+  //
+
+  AliTRDpropagationLayer &p = (AliTRDpropagationLayer &)o; 
+  p.fN   = fN;
+  p.fSec = fSec;
+  p.fX = fX;
+  p.fdX = fdX;
+  p.fRho = fRho;
+  p.fX0  = fX0;
+  p.fTimeBinIndex = fTimeBinIndex;
+  p.fPlane  = fPlane;
+  p.fYmax = fYmax;
+  p.fYmaxSensitive  = fYmaxSensitive;
+  p.fHole = fHole;
+  p.fHoleZc = fHoleZc;
+  p.fHoleZmax = fHoleZmax;
+  p.fHoleYc = fHoleYc;
+  p.fHoleYmax = fHoleYmax;
+  p.fHoleRho = fHoleRho;
+  p.fHoleX0 = fHoleX0;
+
+  for (Int_t i = 0; i < (Int_t)kZones; i++) {
+    p.fZc[i]   = fZc[i]; 
+    p.fZmax[i] = fZmax[i];
+    p.fIsHole[i] = fIsHole[i];
+    p.fZmaxSensitive[i] = fZmaxSensitive[i];  
+  }
+       
+  // Make a deep copy of the Clusters array and the Index array
+  // unless they are needed in class AliTRDstackLayer
+  if (fTimeBinIndex >= 0) { 
+    if (!p.fClusters) 
+      p.fClusters = new AliTRDcluster*[(Int_t)kMaxClusterPerTimeBin];
+    if (!p.fIndex)
+      p.fIndex    = new UInt_t[(Int_t)kMaxClusterPerTimeBin];
+  }
+  for (Int_t i = 0; i < (Int_t)kMaxClusterPerTimeBin; i++){
+    //overwrite
+    p.fClusters[i] = fClusters[i];
+    p.fIndex[i] = fIndex[i];
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDpropagationLayer::SetZ(Double_t *center, Double_t *w, Double_t *wsensitive )
+{
+  //
+  // Set centers and the width of sectors
+  //
+
+  for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
+    fZc[icham]            = center[icham];  
+    fZmax[icham]          = w[icham];
+    fZmaxSensitive[icham] = wsensitive[icham];
+  }  
+
+}
+
+//_____________________________________________________________________________
+void AliTRDpropagationLayer::SetHoles(Bool_t *holes)
+{
+  //
+  // Set centers and the width of sectors
+  //
+
+  fHole = kFALSE;
+
+  for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
+    fIsHole[icham] = holes[icham]; 
+    if (holes[icham]) {
+      fHole = kTRUE;
+    }
+  }  
+
+}
+
+//_____________________________________________________________________________
+void AliTRDpropagationLayer::InsertCluster(AliTRDcluster *c, UInt_t index) 
+{
+  //
+  // Insert cluster in cluster array.
+  // Clusters are sorted according to Y coordinate.  
+  //
+
+  if (fTimeBinIndex < 0) { 
+    //AliWarning("Attempt to insert cluster into non-sensitive time bin!\n");
+    return;
+  }
+
+  if (fN == (Int_t) kMaxClusterPerTimeBin) {
+    //AliWarning("Too many clusters !\n"); 
+    return;
+  }
+
+  if (fN == 0) {
+    fIndex[0]       = index; 
+    fClusters[fN++] = c; 
+    return;
+  }
+
+  Int_t i = Find(c->GetY());
+  memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
+  memmove(fIndex   +i+1,fIndex   +i,(fN-i)*sizeof(UInt_t)); 
+  fIndex[i]    = index; 
+  fClusters[i] = c; 
+  fN++;
+
+}  
+
+//_____________________________________________________________________________
+Int_t AliTRDpropagationLayer::Find(Float_t y) const
+{
+  //
+  // Returns index of the cluster nearest in Y    
+  //
+
+  if (fN <= 0) {
+    return 0;
+  }
+  if (y <= fClusters[0]->GetY()) {
+    return 0;
+  }
+  if (y >  fClusters[fN-1]->GetY()) {
+    return fN;
+  }
+
+  Int_t b = 0;
+  Int_t e = fN - 1;
+  Int_t m = (b + e) / 2;
+
+  for ( ; b < e; m = (b + e) / 2) {
+    if (y > fClusters[m]->GetY()) {
+      b = m + 1;
+    }
+    else {
+      e = m;
+    }
+  }
+
+  return m;
+
+}    
+
+//_____________________________________________________________________________
+Int_t AliTRDpropagationLayer::FindNearestCluster(Float_t y, Float_t z
+                                               , Float_t maxroad
+                                               , Float_t maxroadz) const 
+{
+  //
+  // Returns index of the cluster nearest to the given y,z
+  //
+
+  Int_t   index   = -1;
+  Int_t   maxn    = fN;
+  Float_t mindist = maxroad;                   
+
+  for (Int_t i = Find(y-maxroad); i < maxn; i++) {
+    AliTRDcluster *c = (AliTRDcluster *) (fClusters[i]);
+    Float_t ycl = c->GetY();
+    if (ycl > (y + maxroad)) {
+      break;
+    }
+    if (TMath::Abs(c->GetZ() - z) > maxroadz) {
+      continue;
+    }
+    if (TMath::Abs(ycl - y)       < mindist) {
+      mindist = TMath::Abs(ycl - y);
+      index   = fIndex[i];
+    }
+  }                                            
+
+  return index;
+
+}             
+
diff --git a/TRD/AliTRDpropagationLayer.h b/TRD/AliTRDpropagationLayer.h
new file mode 100644 (file)
index 0000000..36bfe94
--- /dev/null
@@ -0,0 +1,106 @@
+#ifndef ALITRDPROPAGATIONLAYER_H
+#define ALITRDPROPAGATIONLAYER_H   
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */ 
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+//  The TRD propagation layer                                             //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ROOT_TObject
+#include "TObject.h"
+#endif
+
+class AliTRDcluster;
+
+class AliTRDpropagationLayer : public TObject
+{      
+
+ public:       
+
+  enum { kMaxClusterPerTimeBin = 2300
+       , kZones                = 5    };
+       
+  AliTRDpropagationLayer();
+  AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
+                      , Double_t x0, Int_t tbIndex, Int_t plane);
+  AliTRDpropagationLayer(const AliTRDpropagationLayer &/*p*/);
+  ~AliTRDpropagationLayer();
+  AliTRDpropagationLayer &operator=(const AliTRDpropagationLayer &/*p*/) { return *this; }
+
+  operator Int_t() const                                        { return fN;                    }
+  AliTRDcluster *operator[](Int_t i)                            { return fClusters[i];          }
+  virtual void Copy(TObject &o) const;
+       
+  void         SetZmax(Int_t cham, Double_t center, Double_t w) { fZc[cham]      = center;
+                                                                  fZmax[cham]    = w;           }
+  void         SetYmax(Double_t w, Double_t wsensitive)         { fYmax          = w;
+                                                                  fYmaxSensitive = wsensitive;  }
+
+  void         SetZ(Double_t* center, Double_t *w, Double_t *wsensitive);
+  void         SetHoles(Bool_t* holes);
+  //void         SetHole(Double_t Zmax, Double_t Ymax
+  //                  , Double_t rho = 1.29e-3, Double_t x0 = 36.66
+  //                   , Double_t Yc = 0.0, Double_t Zc = 0.0);
+  void         SetSector(const Int_t sec)                       { fSec           = sec;         }
+  void         SetX(Double_t x)                                 { fX             = x;           }
+       
+  Double_t     GetYmax() const                                  { return fYmax;                 }
+  Double_t     GetZmax(Int_t c) const                           { return fZmax[c];              }
+  Double_t     GetZc(Int_t c) const                             { return fZc[c];                }
+  UInt_t       GetIndex(Int_t i) const                          { return fIndex[i];             }
+  Double_t     GetX() const                                     { return fX;                    }
+  Double_t     GetdX() const                                    { return fdX;                   }
+  Int_t        GetTimeBinIndex() const                          { return fTimeBinIndex;         }
+  Int_t        GetPlane() const                                 { return fPlane;                }
+  Int_t        GetSector() const                                { return fSec;                  }
+
+  Bool_t       IsHole(Int_t zone) const                         { return fIsHole[zone];         }
+  Bool_t       IsSensitive() const                              { return (fTimeBinIndex >= 0) ? kTRUE : kFALSE;}
+
+  void         Clear(const Option_t * /*o*/)                    { ; } 
+  void         Clear()                                          { for (Int_t i = 0; i < fN; i++) 
+                                                                    fClusters[i] = NULL;
+                                                                 fN = 0;                       }
+       
+  void         InsertCluster(AliTRDcluster *c, UInt_t index);
+  Int_t        Find(Float_t y) const;
+  Int_t        FindNearestCluster(Float_t y, Float_t z, Float_t maxroad, Float_t maxroadz) const;
+
+ protected:
+
+       Int_t         fN;                     // Number of clusters
+       Int_t         fSec;                   // Sector mumber
+       AliTRDcluster **fClusters;            //[fN] Array of pointers to clusters
+       UInt_t        *fIndex;                //[fN] Array of cluster indexes
+       Double_t      fX;                     // X coordinate of the middle plane
+       Double_t      fdX;                    // Radial thickness of the time bin
+       Double_t      fRho;                   // Default density of the material
+       Double_t      fX0;                    // Default radiation length
+       Int_t         fTimeBinIndex;          // Plane * F(local_tb)
+       Int_t         fPlane;                 // Plane number
+       Double_t      fZc[kZones];            // Z position of the center for 5 active areas
+       Double_t      fZmax[kZones];          // Half of active area length in Z
+       Double_t      fZmaxSensitive[kZones]; // Sensitive area for detection Z
+       Bool_t        fIsHole[kZones];        // Is hole in given sector
+       Double_t      fYmax;                  // Half of active area length in Y
+       Double_t      fYmaxSensitive;         // Half of active area length in Y
+       Bool_t        fHole;                  // kTRUE if there is a hole in the layer
+       Double_t      fHoleZc;                // Z of the center of the hole
+       Double_t      fHoleZmax;              // Half of the hole length in Z
+       Double_t      fHoleYc;                // Y of the center of the hole
+       Double_t      fHoleYmax;              // Half of the hole length in Y
+       Double_t      fHoleRho;               // Density of the gas in the hole
+       Double_t      fHoleX0;                // Radiation length of the gas in the hole
+
+       ClassDef(AliTRDpropagationLayer, 1)   // Propagation layer for TRD tracker
+
+};
+
+#endif 
+
diff --git a/TRD/AliTRDrecoParam.cxx b/TRD/AliTRDrecoParam.cxx
new file mode 100644 (file)
index 0000000..fff4343
--- /dev/null
@@ -0,0 +1,75 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Parameter class for the TRD reconstruction                               //
+//                                                                           //
+//  Authors:                                                                 //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                        //
+//    Markus Fasel <M.Fasel@gsi.de>                                          //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDrecoParam.h"
+
+ClassImp(AliTRDrecoParam)
+
+//______________________________________________________________
+AliTRDrecoParam::AliTRDrecoParam()
+  :AliDetectorRecoParam()
+  ,fkMaxTheta(1.0)
+  ,fkMaxPhi(2.0)
+  ,fkRoad0y(6.0)
+  ,fkRoad0z(8.5) 
+  ,fkRoad1y(2.0)
+  ,fkRoad1z(20.0)      
+  ,fkRoad2y(3.0)
+  ,fkRoad2z(20.0)
+  ,fkPlaneQualityThreshold(5.0)// 4.2? under Investigation
+  ,fkFindable(.333)
+  ,fkChi2Z(14./*12.5*/)
+  ,fkChi2Y(.25)
+  ,fkTrackLikelihood(-15.)
+{
+  //
+  // Default constructor
+  //
+
+}
+
+//______________________________________________________________
+AliTRDrecoParam *AliTRDrecoParam::GetLowFluxParam()
+{
+  //
+  // Parameters for the low flux environment
+  //
+
+  return new AliTRDrecoParam();
+
+}
+
+//______________________________________________________________
+AliTRDrecoParam *AliTRDrecoParam::GetHighFluxParam()
+{
+  //
+  // Parameters for the high flux environment
+  //
+
+  return new AliTRDrecoParam();
+
+}
diff --git a/TRD/AliTRDrecoParam.h b/TRD/AliTRDrecoParam.h
new file mode 100644 (file)
index 0000000..96b6bee
--- /dev/null
@@ -0,0 +1,71 @@
+#ifndef ALITRDRECOPARAM_H
+#define ALITRDRECOPARAM_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+//  Parameter class for the TRD reconstruction                            //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALIDETECTORRECOPARAM_H
+#include "AliDetectorRecoParam.h"
+#endif
+
+class AliTRDrecoParam : public AliDetectorRecoParam
+{
+
+  public:
+
+       AliTRDrecoParam();
+       ~AliTRDrecoParam() { }
+
+       Double_t GetChi2Y() const                 { return fkChi2Y;    }
+       Double_t GetChi2Z() const                 { return fkChi2Z;    }
+       Double_t GetFindableClusters() const      { return fkFindable; }
+       Double_t GetMaxTheta() const              { return fkMaxTheta; }
+       Double_t GetMaxPhi() const                { return fkMaxPhi;   }
+
+       Double_t GetRoad0y() const                { return fkRoad0y;   }
+       Double_t GetRoad0z() const                { return fkRoad0z;   }
+
+       Double_t GetRoad1y() const                { return fkRoad1y;   }
+       Double_t GetRoad1z() const                { return fkRoad1z;   }
+
+       Double_t GetRoad2y() const                { return fkRoad2y;   }
+       Double_t GetRoad2z() const                { return fkRoad2z;   }
+
+       Double_t GetPlaneQualityThreshold() const { return fkPlaneQualityThreshold; }
+
+       Double_t GetTrackLikelihood() const       { return fkTrackLikelihood;       }
+       
+       static   AliTRDrecoParam *GetLowFluxParam();
+        static   AliTRDrecoParam *GetHighFluxParam();
+
+ private:
+
+       Double_t fkMaxTheta;              // Maximum theta
+       Double_t fkMaxPhi;                // Maximum phi
+
+       Double_t fkRoad0y;                // Road for middle cluster
+       Double_t fkRoad0z;                // Road for middle cluster
+
+       Double_t fkRoad1y;                // Road in y for seeded cluster
+       Double_t fkRoad1z;                // Road in z for seeded cluster
+
+       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
+       Double_t fkChi2Y;                 // Max chi2 on the y direction for seeding clusters Rieman fit
+       Double_t fkTrackLikelihood;       // Track likelihood for tracklets Rieman fit
+
+       ClassDef(AliTRDrecoParam, 1)      // Reconstruction parameters for TRD detector
+
+};
+#endif
index 54b5505..f7c5357 100644 (file)
@@ -27,6 +27,7 @@
 #include "AliMathBase.h"
 
 #include "AliTRDseed.h"
+#include "AliTRDcalibDB.h"
 #include "AliTRDcluster.h"
 #include "AliTRDtracker.h"
 
@@ -57,7 +58,7 @@ AliTRDseed::AliTRDseed()
   // Default constructor  
   //
 
-  for (Int_t i = 0; i < 25; i++) {
+  for (Int_t i = 0; i < knTimebins; i++) {
     fX[i]        = 0;   // x position
     fY[i]        = 0;   // y position
     fZ[i]        = 0;   // z position
@@ -103,7 +104,7 @@ AliTRDseed::AliTRDseed(const AliTRDseed &s)
   // Copy constructor  
   //
 
-  for (Int_t i = 0; i < 25; i++) {
+  for (Int_t i = 0; i < knTimebins; i++) {
     fX[i]        = s.fX[i];        // x position
     fY[i]        = s.fY[i];        // y position
     fZ[i]        = s.fZ[i];        // z position
@@ -125,18 +126,61 @@ AliTRDseed::AliTRDseed(const AliTRDseed &s)
 }
 
 //_____________________________________________________________________________
+void AliTRDseed::Copy(TObject &o) const {
+       AliTRDseed &seed = (AliTRDseed &)o;
+       seed.fTilt = fTilt;
+  seed.fPadLength = fPadLength;
+  seed.fX0 = fX0;
+  seed.fSigmaY = fSigmaY;
+  seed.fSigmaY2 = fSigmaY2;
+  seed.fMeanz = fMeanz;
+  seed.fZProb = fZProb;
+  seed.fN = fN;
+  seed.fN2 = fN2;
+  seed.fNUsed = fNUsed;
+  seed.fFreq = fFreq;
+  seed.fNChange = fNChange;
+  seed.fMPads = fMPads;
+  seed.fC = fC;
+  seed.fCC = fCC;
+  seed.fChi2 = fChi2;
+  seed.fChi2Z = fChi2Z;
+       for (Int_t i = 0; i < knTimebins; i++) {
+    seed.fX[i]        = fX[i];
+    seed.fY[i]        = fY[i]; 
+    seed.fZ[i]        = fZ[i]; 
+    seed.fIndexes[i]  = fIndexes[i]; 
+    seed.fClusters[i] = fClusters[i]; 
+    seed.fUsable[i]   = fUsable[i]; 
+  }
+
+  for (Int_t i = 0; i < 2; i++) {
+    seed.fYref[i]     = fYref[i];
+    seed.fZref[i]     = fZref[i];
+    seed.fYfit[i]     = fYfit[i];
+    seed.fYfitR[i]    = fYfitR[i]; 
+    seed.fZfit[i]     = fZfit[i]; 
+    seed.fZfitR[i]    = fZfitR[i];  
+    seed.fLabels[i]   = fLabels[i]; 
+  }
+
+  TObject::Copy(seed);
+
+}
+
+//_____________________________________________________________________________
 void AliTRDseed::Reset()
 {
   //
   // Reset seed
   //
 
-  for (Int_t i = 0; i < 25; i++) {
+  for (Int_t i = 0; i < knTimebins; i++) {
     fX[i]        = 0;  // X position
     fY[i]        = 0;  // Y position
     fZ[i]        = 0;  // Z position
     fIndexes[i]  = 0;  // Indexes
-    fClusters[i] = 0;  // Clusters
+    fClusters[i] = 0x0;  // Clusters
     fUsable[i]   = kFALSE;    
   }
 
@@ -168,11 +212,14 @@ 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 < 25; i++) {
+  for (Int_t i = 0; i < nTimeBins+1; i++) {
     if (!fClusters[i]) continue;
     for (Int_t ilab = 0; ilab < 3; ilab++) {
       if (fClusters[i]->GetLabel(ilab) >= 0) {
@@ -198,7 +245,10 @@ void AliTRDseed::UseClusters()
   // Use clusters
   //
 
-  for (Int_t i = 0; i < 25; i++) {
+  AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+  Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+  for (Int_t i = 0; i < nTimeBins+1; i++) {
     if (!fClusters[i]) continue;
     if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
   }
@@ -213,10 +263,16 @@ void AliTRDseed::Update()
   //
 
   const Float_t kRatio  = 0.8;
-  const Int_t   kClmin  = 6;
+  const Int_t   kClmin  = 5;
   const Float_t kmaxtan = 2;
 
-  if (TMath::Abs(fYref[1]) > kmaxtan) return;              // Track inclined too much
+  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);
+               return;              // Track inclined too much
+       }
 
   Float_t  sigmaexp  = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
   Float_t  ycrosscor = fPadLength * fTilt * 0.5;           // Y correction for crossing 
@@ -230,17 +286,18 @@ void AliTRDseed::Update()
   Double_t sumwz;
   Double_t sumwxz;
 
-  Int_t    zints[25];                    // Histograming of the z coordinate 
+       // Buffering: Leave it constant fot Performance issues
+  Int_t    zints[knTimebins];            // Histograming of the z coordinate 
                                          // Get 1 and second max probable coodinates in z
-  Int_t    zouts[50];       
-  Float_t  allowedz[25];                 // Allowed z for given time bin
-  Float_t  yres[25];                     // Residuals from reference
+  Int_t    zouts[2*knTimebins];       
+  Float_t  allowedz[knTimebins];         // Allowed z for given time bin
+  Float_t  yres[knTimebins];             // Residuals from reference
   Float_t  anglecor = fTilt * fZref[1];  // Correction to the angle
   
   
   fN  = 0; 
   fN2 = 0;
-  for (Int_t i = 0; i < 25; i++) {
+  for (Int_t i = 0; i < nTimeBins; i++) {
     yres[i] = 10000.0;
     if (!fClusters[i]) continue;
     yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i];   // Residual y
@@ -248,11 +305,17 @@ void AliTRDseed::Update()
     fN++;    
   }
 
-  if (fN < kClmin) return;
+  if (fN < kClmin){
+               //printf("Exit fN < kClmin: fN = %d\n", fN);
+               return; 
+       }
   Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
   fZProb   = zouts[0];
   if (nz <= 1) zouts[3] = 0;
-  if (zouts[1] + zouts[3] < kClmin) return;
+  if (zouts[1] + zouts[3] < kClmin) {
+               //printf("Exit zouts[1] = %d, zouts[3] = %d\n",zouts[1],zouts[3]);
+               return;
+       }
   
   // Z distance bigger than pad - length
   if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) {
@@ -261,7 +324,7 @@ void AliTRDseed::Update()
   
   Int_t  breaktime = -1;
   Bool_t mbefore   = kFALSE;
-  Int_t  cumul[25][2];
+  Int_t  cumul[knTimebins][2];
   Int_t  counts[2] = { 0, 0 };
   
   if (zouts[3] >= 3) {
@@ -271,22 +334,22 @@ void AliTRDseed::Update()
     // with maximal numebr of accepted clusters
     //
     fNChange = 1;
-    for (Int_t i = 0; i < 25; i++) {
+    for (Int_t i = 0; i < nTimeBins; 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 < 24; i++) {
-      Int_t after  = cumul[24][0] - cumul[i][0];
+    for (Int_t i = 0; i < nTimeBins; i++) {
+      Int_t after  = cumul[nTimeBins][0] - cumul[i][0];
       Int_t before = cumul[i][1];
       if (after + before > maxcount) { 
        maxcount  = after + before; 
        breaktime = i;
        mbefore   = kFALSE;
       }
-      after  = cumul[24][1] - cumul[i][1];
+      after  = cumul[nTimeBins-1][1] - cumul[i][1];
       before = cumul[i][0];
       if (after + before > maxcount) { 
        maxcount  = after + before; 
@@ -299,18 +362,18 @@ void AliTRDseed::Update()
 
   }
 
-  for (Int_t i = 0; i < 25; i++) {
+  for (Int_t i = 0; i < nTimeBins+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[24]) && (fZref[1] < 0)) || 
-      ((allowedz[0] < allowedz[24]) && (fZref[1] > 0))) {
+  if (((allowedz[0] > allowedz[nTimeBins]) && (fZref[1] < 0)) || 
+      ((allowedz[0] < allowedz[nTimeBins]) && (fZref[1] > 0))) {
     //
     // Tracklet z-direction not in correspondance with track z direction 
     //
     fNChange = 0;
-    for (Int_t i = 0; i < 25; i++) {
+    for (Int_t i = 0; i < nTimeBins+1; i++) {
       allowedz[i] = zouts[0];  // Only longest taken
     } 
   }
@@ -319,7 +382,7 @@ void AliTRDseed::Update()
     //
     // Cross pad -row tracklet  - take the step change into account
     //
-    for (Int_t i = 0; i < 25; i++) {
+    for (Int_t i = 0; i < nTimeBins+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
@@ -330,20 +393,21 @@ void AliTRDseed::Update()
     }
   }
   
-  Double_t yres2[25];
+  Double_t yres2[knTimebins];
   Double_t mean;
   Double_t sigma;
-  for (Int_t i = 0; i < 25; i++) {
+  for (Int_t i = 0; i < nTimeBins+1; i++) {
     if (!fClusters[i]) continue;
     if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
     yres2[fN2] = yres[i];
     fN2++;
   }
   if (fN2 < kClmin) {
+               //printf("Exit fN2 < kClmin: fN2 = %d\n", fN2);
     fN2 = 0;
     return;
   }
-  AliMathBase::EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*kRatio-2));
+  AliMathBase::EvaluateUni(fN2,yres2,mean,sigma, Int_t(fN2*kRatio-2.));
   if (sigma < sigmaexp * 0.8) {
     sigma = sigmaexp;
   }
@@ -362,7 +426,7 @@ void AliTRDseed::Update()
   fMeanz = 0;
   fMPads = 0;
 
-  for (Int_t i = 0; i < 25; i++) {
+  for (Int_t i = 0; i < nTimeBins+1; i++) {
 
     fUsable[i] = kFALSE;
     if (!fClusters[i]) continue;
@@ -387,6 +451,7 @@ void AliTRDseed::Update()
   }
 
   if (fN2 < kClmin){
+               //printf("Exit fN2 < kClmin(2): fN2 = %d\n",fN2);
     fN2 = 0;
     return;
   }
@@ -403,7 +468,7 @@ void AliTRDseed::Update()
   fYfitR[1]    = (sumw   * sumwxy - sumwx * sumwy)  / det;
   
   fSigmaY2 = 0;
-  for (Int_t i = 0; i < 25; i++) {    
+  for (Int_t i = 0; i < nTimeBins+1; i++) {    
     if (!fUsable[i]) continue;
     Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
     fSigmaY2 += delta*delta;
@@ -430,8 +495,11 @@ void AliTRDseed::UpdateUsed()
   // Update used seed
   //
 
+  AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+  Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
   fNUsed = 0;
-  for (Int_t i = 0; i < 25; i++) {
+  for (Int_t i = 0; i < nTimeBins; i++) {
     if (!fClusters[i]) {
       continue;
     }
@@ -452,6 +520,10 @@ Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
   // Fitting with tilting pads - kz not fixed
   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;
@@ -462,7 +534,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 < 25; itime++) {
+    for (Int_t itime = 0; itime < nTimeBins+1; itime++) {
 
       if (!cseed[iLayer].fUsable[itime]) continue;
       // x relative to the midle chamber
@@ -506,12 +578,9 @@ Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
   //       
   Bool_t acceptablez = kTRUE;
   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
-    if (cseed[iLayer].IsOK()) {
-      Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
-      if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) {
-       acceptablez = kFALSE;
-      }
-    }
+    if (!cseed[iLayer].IsOK()) continue;
+    Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
+    if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) acceptablez = kFALSE;
   }
   if (!acceptablez) {
     Double_t zmf  = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0);
@@ -532,8 +601,8 @@ Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
   params[2] =  fitterT2.GetParameter(2);           
   Double_t curvature =  1.0 + params[1] * params[1] - params[2] * params[0];
 
+       
   for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
-
     Double_t  x  = cseed[iLayer].fX0;
     Double_t  y  = 0;
     Double_t  dy = 0;
@@ -554,9 +623,9 @@ Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
     if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) {
       Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1); 
       if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) {
-       Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0));
-       if (params[0] < 0) res *= -1.0;
-       dy = res;
+                               Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0));
+                               if (params[0] < 0) res *= -1.0;
+                               dy = res;
       }
     }
     z  = rpolz0 + rpolz1 * (x - xref2);
index f9af380..11e6e69 100644 (file)
@@ -19,6 +19,7 @@ class AliTRDcluster;
 class AliTRDseed : public TObject {
 
  public:
+       enum { knTimebins = 35 };
 
   AliTRDseed(); 
   AliTRDseed(const AliTRDseed &s);
@@ -27,7 +28,6 @@ class AliTRDseed : public TObject {
   AliTRDseed      &operator=(const AliTRDseed &s)           { *(new(this) AliTRDseed(s)); 
                                                               return *this;          }
 
-
   static  Float_t  FitRiemanTilt(AliTRDseed *seed, Bool_t error);
           void     UseClusters();
           void     Update();
@@ -35,7 +35,7 @@ class AliTRDseed : public TObject {
           void     UpdateUsed();
           void     Reset();
 
-          Bool_t   IsOK() const                             { return fN2 > 8;        }
+          Bool_t   IsOK() const                             { return fN2 > 4;        }
           Bool_t   IsUsable(Int_t i) const                  { return fUsable[i];     }
  
           Float_t  GetTilt() const                          { return fTilt;          }
@@ -93,17 +93,19 @@ class AliTRDseed : public TObject {
           void     SetChi2(Float_t chi2)                    { fChi2        = chi2;   }
           void     SetChi2Z(Float_t chi2z)                  { fChi2Z       = chi2z;  }
 
- private:
+ protected:
 
+          void     Copy(TObject &o) const;
+          
           Float_t  fTilt;               //  Tilting angle
           Float_t  fPadLength;          //  Pad length
           Float_t  fX0;                 //  X0 position
-          Float_t  fX[25];              //! X position
-          Float_t  fY[25];              //! Y position
-          Float_t  fZ[25];              //! Z position
-          Int_t    fIndexes[25];        //! Indexes
-          AliTRDcluster *fClusters[25]; //! Clusters
-          Bool_t   fUsable[25];         //! Indication  - usable cluster
+          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
diff --git a/TRD/AliTRDseedV1.cxx b/TRD/AliTRDseedV1.cxx
new file mode 100644 (file)
index 0000000..039826e
--- /dev/null
@@ -0,0 +1,683 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+//  The TRD track seed                                                    //
+//                                                                        //
+//  Authors:                                                              //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                     //
+//    Markus Fasel <M.Fasel@gsi.de>                                       //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#include "TMath.h"
+#include "TLinearFitter.h"
+
+#include "AliLog.h"
+#include "AliMathBase.h"
+
+#include "AliTRDseedV1.h"
+#include "AliTRDcluster.h"
+#include "AliTRDcalibDB.h"
+#include "AliTRDstackLayer.h"
+#include "AliTRDrecoParam.h"
+
+#define SEED_DEBUG
+
+ClassImp(AliTRDseedV1)
+
+//____________________________________________________________________
+AliTRDseedV1::AliTRDseedV1(Int_t layer, AliTRDrecoParam *p) 
+  :AliTRDseed()
+  ,fLayer(layer)
+  ,fTimeBins(0)
+  ,fOwner(kFALSE)
+  ,fRecoParam(p)
+{
+  //
+  // Constructor
+  //
+
+       //AliInfo("");
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       fTimeBins = cal->GetNumberOfTimeBins();
+
+}
+
+//____________________________________________________________________
+AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref, Bool_t owner)
+  :AliTRDseed((AliTRDseed&)ref)
+  ,fLayer(ref.fLayer)
+  ,fTimeBins(ref.fTimeBins)
+  ,fOwner(kFALSE)
+  ,fRecoParam(ref.fRecoParam)
+{
+  //
+  // Copy Constructor performing a deep copy
+  //
+
+       //AliInfo("");
+
+       if(owner){
+               for(int ic=0; ic<fTimeBins; ic++){
+                       if(!fClusters[ic]) continue;
+                       fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
+               }
+               fOwner = kTRUE;
+       }
+
+}
+
+//____________________________________________________________________
+AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
+{
+  //
+  // Assignment Operator using the copy function
+  //
+
+       //AliInfo("");
+       if(this != &ref){
+               ref.Copy(*this);
+       }
+       return *this;
+
+}
+
+//____________________________________________________________________
+AliTRDseedV1::~AliTRDseedV1()
+{
+  //
+  // Destructor. The RecoParam object belongs to the underlying tracker.
+  //
+
+       //AliInfo(Form("fOwner[%s]", fOwner?"YES":"NO"));
+
+       if(fOwner) delete [] fClusters;
+}
+
+//____________________________________________________________________
+void AliTRDseedV1::Copy(TObject &ref) const
+{
+  //
+  // Copy function
+  //
+
+       //AliInfo("");
+       AliTRDseedV1 &target = (AliTRDseedV1 &)ref; 
+       
+       target.fLayer     = fLayer;
+       target.fTimeBins  = fTimeBins;
+       target.fRecoParam = fRecoParam;
+       AliTRDseed::Copy(target);
+}
+
+//____________________________________________________________________
+Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
+{
+  //
+  // Returns a quality measurement of the current seed
+  //
+
+       Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
+       return .5 * (18.0 - fN2)
+               + 10.* TMath::Abs(fYfit[1] - fYref[1])
+               + 5.* TMath::Abs(fYfit[0] - fYref[0] + zcorr)
+               + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength;
+}
+
+//____________________________________________________________________
+Bool_t AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
+                                       , Float_t quality
+                                       , Bool_t kZcorr
+                                       , AliTRDcluster *c)
+{
+  //
+  // Iterative process to register clusters to the seed.
+  // In iteration 0 we try only one pad-row and if quality not
+  // sufficient we try 2 pad-rows (about 5% of tracks cross 2 pad-rows)
+  //
+       
+       if(!fRecoParam){
+               AliError("Seed can not be used without a valid RecoParam.");
+               return kFALSE;
+       }
+       
+       Float_t  tquality;
+       Double_t kroady = fRecoParam->GetRoad1y();
+       Double_t kroadz = fPadLength * .5 + 1.;
+       
+       // initialize configuration parameters
+       Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
+       Int_t   niter = kZcorr ? 1 : 2;
+       
+       Double_t yexp, zexp;
+       Int_t ncl = 0;
+       // start seed update
+       for (Int_t iter = 0; iter < niter; iter++) {
+       //AliInfo(Form("iter = %i", iter));
+               ncl = 0;
+               for (Int_t iTime = 0; iTime < fTimeBins; iTime++) {
+                       // define searching configuration
+                       Double_t dxlayer = layer[iTime].GetX() - fX0;
+                       if(c){
+                               zexp = c->GetZ();
+                               //Try 2 pad-rows in second iteration
+                               if (iter > 0) {
+                                       zexp = fZref[0] + fZref[1] * dxlayer - zcorr;
+                                       if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
+                                       if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
+                               }
+                       } else zexp = fZref[0];
+                       yexp  = fYref[0] + fYref[1] * dxlayer - zcorr;
+                       // get  cluster
+//                     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);
+//                     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("Index = %i\n",index);
+                       if (index < 0) continue;
+                       
+                       // Register cluster
+                       AliTRDcluster *cl = (AliTRDcluster*) layer[iTime].GetCluster(index);
+                       
+                       //printf("Cluster %i(0x%x): x = %3.3f, y = %3.3f, z = %3.3f\n", index, cl, cl->GetX(), cl->GetY(), cl->GetZ());
+
+                       Int_t GlobalIndex = layer[iTime].GetGlobalIndex(index);
+                       fIndexes[iTime]  = GlobalIndex;
+                       fClusters[iTime] = cl;
+                       fX[iTime]        = dxlayer;
+                       fY[iTime]        = cl->GetY();
+                       fZ[iTime]        = cl->GetZ();
+       
+                       // Debugging
+                       ncl++;
+               }
+
+#ifdef SEED_DEBUG
+//             Int_t nclusters = 0;
+//             Float_t fD[iter] = 0.;
+//             for(int ic=0; ic<fTimeBins+1; ic++){
+//                     AliTRDcluster *ci = fClusters[ic];
+//                     if(!ci) continue;
+//                     for(int jc=ic+1; jc<fTimeBins+1; jc++){
+//                             AliTRDcluster *cj = fClusters[jc];
+//                             if(!cj) continue;
+//                             fD[iter] += TMath::Sqrt((ci->GetY()-cj->GetY())*(ci->GetY()-cj->GetY())+
+//                             (ci->GetZ()-cj->GetZ())*(ci->GetZ()-cj->GetZ()));
+//                             nclusters++;
+//                     }
+//             }
+//             if(nclusters) fD[iter] /= float(nclusters);
+#endif
+
+               AliTRDseed::Update();
+
+               if(IsOK()){
+                       tquality = GetQuality(kZcorr);
+                       if(tquality < quality) break;
+                       else quality = tquality;
+               }
+               kroadz *= 2.;
+       } // Loop: iter
+       if (!IsOK()) return kFALSE;
+
+       CookLabels();
+       UpdateUsed();
+       return kTRUE;   
+}
+
+//____________________________________________________________________
+Bool_t AliTRDseedV1::AttachClustersProj(AliTRDstackLayer *layer
+                                       , Float_t /*quality*/
+                                       , Bool_t kZcorr
+                                       , AliTRDcluster *c)
+{
+  //
+  // Projective algorithm to attach clusters to seeding tracklets
+  //
+  // Parameters
+  //
+  // Output
+  //
+  // Detailed description
+  // 1. Collapse x coordinate for the full detector plane
+  // 2. truncated mean on y (r-phi) direction
+  // 3. purge clusters
+  // 4. truncated mean on z direction
+  // 5. purge clusters
+  // 6. fit tracklet
+  //   
+
+       if(!fRecoParam){
+               AliError("Seed can not be used without a valid RecoParam.");
+               return kFALSE;
+       }
+
+       const Int_t knTimeBins = 35;
+       const Int_t kClusterCandidates = 2 * knTimeBins;
+       
+       //define roads
+       Double_t kroady = fRecoParam->GetRoad1y();
+       Double_t kroadz = fPadLength * 1.5 + 1.;
+       // correction to y for the tilting angle
+       Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
+
+       // working variables
+       AliTRDcluster *clusters[kClusterCandidates];
+       Double_t cond[4], yexp[knTimeBins], zexp[knTimeBins],
+               yres[kClusterCandidates], zres[kClusterCandidates];
+       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++) {
+               fX[iTime] = layer[iTime].GetX() - fX0;
+               zexp[iTime] = fZref[0] + fZref[1] * fX[iTime];
+               yexp[iTime] = fYref[0] + fYref[1] * fX[iTime] - zcorr;
+               
+               // build condition and process clusters
+               cond[0] = yexp[iTime] - kroady; cond[1] = yexp[iTime] + kroady;
+               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]);
+                       clusters[nYclusters] = c;
+                       yres[nYclusters++] = c->GetY() - yexp[iTime];
+                       if(nYclusters >= kClusterCandidates) {
+                               AliWarning(Form("Cluster candidates reached limit %d. Some may be lost.", kClusterCandidates));
+                               kEXIT = kTRUE;
+                               break;
+                       }
+               }
+               tboundary[iTime] = nYclusters;
+               if(kEXIT) break;
+       }
+       
+       // Evaluate truncated mean on the y direction
+       Double_t mean, sigma;
+       AliMathBase::EvaluateUni(nYclusters, yres, mean, sigma, Int_t(nYclusters*.8)-2);
+       //purge cluster candidates
+       Int_t nZclusters = 0;
+       for(Int_t ic = 0; ic<nYclusters; ic++){
+               if(yres[ic] - mean > 4. * sigma){
+                       clusters[ic] = 0x0;
+                       continue;
+               }
+               zres[nZclusters++] = clusters[ic]->GetZ() - zexp[clusters[ic]->GetLocalTimeBin()];
+       }
+       
+       // Evaluate truncated mean on the z direction
+       AliMathBase::EvaluateUni(nZclusters, zres, mean, sigma, Int_t(nZclusters*.8)-2);
+       //purge cluster candidates
+       for(Int_t ic = 0; ic<nZclusters; ic++){
+               if(zres[ic] - mean > 4. * sigma){
+                       clusters[ic] = 0x0;
+                       continue;
+               }
+       }
+
+       
+       // Select only one cluster/TimeBin
+       Int_t lastCluster = 0;
+       fN2 = 0;
+       for (Int_t iTime = 0; iTime < fTimeBins; iTime++) {
+               ncl = tboundary[iTime] - lastCluster;
+               if(!ncl) continue;
+               if(ncl == 1){
+                       c = clusters[lastCluster];
+               } else if(ncl > 1){
+                       Float_t dold = 9999.; Int_t iptr = lastCluster;
+                       for(int ic=lastCluster; ic<tboundary[iTime]; ic++){
+                               if(!clusters[ic]) continue;
+                               Float_t y = yexp[iTime] - clusters[ic]->GetY();
+                               Float_t z = zexp[iTime] - clusters[ic]->GetZ();
+                               Float_t d = y * y + z * z;
+                               if(d > dold) continue;
+                               dold = d;
+                               iptr = ic;
+                       }
+                       c = clusters[iptr];
+               }
+               //Int_t GlobalIndex = layer[iTime].GetGlobalIndex(index);
+               //fIndexes[iTime]  = GlobalIndex;
+               fClusters[iTime] = c;
+               fY[iTime]        = c->GetY();
+               fZ[iTime]        = c->GetZ();
+               lastCluster = tboundary[iTime];
+               fN2++;
+       }
+       
+       // number of minimum numbers of clusters expected for the tracklet
+       Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fTimeBins);
+  if (fN2 < kClmin){
+               AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
+    fN2 = 0;
+    return kFALSE;
+  }
+       AliTRDseed::Update();
+       
+//     // fit tracklet and update clusters
+//     if(!FitTracklet()) return kFALSE;
+//     UpdateUsed();
+       return kTRUE;
+}
+
+//____________________________________________________________________
+Bool_t AliTRDseedV1::FitTracklet()
+{
+  //
+  // Linear fit of the tracklet
+  //
+  // Parameters :
+  //
+  // Output :
+  //  True if successful
+  //
+  // Detailed description
+  // 2. Check if tracklet crosses pad row boundary
+  // 1. Calculate residuals in the y (r-phi) direction
+  // 3. Do a Least Square Fit to the data
+  //
+
+       //Float_t  sigmaexp  = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
+  Float_t  ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing
+  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];//
+       
+       fN = 0;
+       for (Int_t iTime = 0; iTime < fTimeBins; iTime++) {
+    if (!fClusters[iTime]) continue;
+    yres[iTime] = fY[iTime] - fYref[0] - (fYref[1] + anglecor) * fX[iTime];
+               zint[fN++] = Int_t(fZ[iTime]);
+       }
+
+       // calculate pad row boundary crosses
+       Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fTimeBins);
+       Int_t nz = AliMathBase::Freq(fN, zint, zout, kFALSE);
+  fZProb   = zout[0];
+  if(nz <= 1) zout[3] = 0;
+  if(zout[1] + zout[3] < kClmin) {
+               AliWarning(Form("Not enough clusters to fit the cross boundary tracklet %d [%d].", zout[1]+zout[3], kClmin));
+               return kFALSE;
+       }
+  // Z distance bigger than pad - length
+  if (TMath::Abs(zout[0]-zout[2]) > fPadLength) zout[3]=0;
+
+  Double_t sumw   = 0., 
+               sumwx  = 0.,
+               sumwx2 = 0.,
+               sumwy  = 0.,
+               sumwxy = 0.,
+               sumwz  = 0.,
+               sumwxz = 0.;
+       Int_t npads;
+  fMPads = 0;
+       fMeanz = 0.;
+       for(int iTime=0; iTime<fTimeBins; iTime++){
+    fUsable[iTime] = kFALSE;
+    if (!fClusters[iTime]) continue;
+               npads = fClusters[iTime]->GetNPads();
+
+               fUsable[iTime] = kTRUE;
+    fN2++;
+    fMPads += npads;
+    Float_t weight = 1.0;
+    if(npads > 5) weight = 0.2;
+    else if(npads > 4) weight = 0.5;
+    sumw   += weight; 
+    sumwx  += fX[iTime] * weight;
+    sumwx2 += fX[iTime] * fX[iTime] * weight;
+    sumwy  += weight * yres[iTime];
+    sumwxy += weight * yres[iTime] * fX[iTime];
+    sumwz  += weight * fZ[iTime];    
+    sumwxz += weight * fZ[iTime] * fX[iTime];
+       }
+  if (fN2 < kClmin){
+               AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
+    fN2 = 0;
+    return kFALSE;
+  }
+  fMeanz = sumwz / sumw;
+       fNChange = 0;
+
+       // Tracklet on boundary
+  Float_t correction = 0;
+  if (fNChange > 0) {
+    if (fMeanz < fZProb) correction =  ycrosscor;
+    if (fMeanz > fZProb) correction = -ycrosscor;
+  }
+
+  Double_t det = sumw * sumwx2 - sumwx * sumwx;
+  fYfitR[0]    = (sumwx2 * sumwy  - sumwx * sumwxy) / det;
+  fYfitR[1]    = (sumw   * sumwxy - sumwx * sumwy)  / det;
+  
+  fSigmaY2 = 0;
+  for (Int_t i = 0; i < fTimeBins+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));
+  
+  fZfitR[0]  = (sumwx2 * sumwz  - sumwx * sumwxz) / det;
+  fZfitR[1]  = (sumw   * sumwxz - sumwx * sumwz)  / det;
+  fZfit[0]   = (sumwx2 * sumwz  - sumwx * sumwxz) / det;
+  fZfit[1]   = (sumw   * sumwxz - sumwx * sumwz)  / det;
+  fYfitR[0] += fYref[0] + correction;
+  fYfitR[1] += fYref[1];
+  fYfit[0]   = fYfitR[0];
+  fYfit[1]   = fYfitR[1];
+
+       return kTRUE;
+}
+
+//_____________________________________________________________________________
+Float_t AliTRDseedV1::FitRiemanTilt(AliTRDseedV1 *cseed, Bool_t terror)
+{
+  //
+  // Fit the Rieman tilt
+  //
+
+  // Fitting with tilting pads - kz not fixed
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       Int_t nTimeBins = cal->GetNumberOfTimeBins();
+  TLinearFitter fitterT2(4,"hyp4");  
+  fitterT2.StoreData(kTRUE);
+  Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
+  
+  Int_t npointsT = 0;
+  fitterT2.ClearPoints();
+
+  for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+//             printf("\nLayer %d\n", iLayer);
+//     cseed[iLayer].Print();
+               if (!cseed[iLayer].IsOK()) continue;
+    Double_t tilt = cseed[iLayer].fTilt;
+
+    for (Int_t itime = 0; itime < nTimeBins+1; itime++) {
+//                     printf("\ttime %d\n", itime);
+      if (!cseed[iLayer].fUsable[itime]) continue;
+      // x relative to the midle chamber
+      Double_t x = cseed[iLayer].fX[itime] + cseed[iLayer].fX0 - xref2;  
+      Double_t y = cseed[iLayer].fY[itime];
+      Double_t z = cseed[iLayer].fZ[itime];
+
+      //
+      // Tilted rieman
+      //
+      Double_t uvt[6];
+      Double_t x2 = cseed[iLayer].fX[itime] + cseed[iLayer].fX0;      // Global x
+      Double_t t  = 1.0 / (x2*x2 + y*y);
+      uvt[1]  = t;
+      uvt[0]  = 2.0 * x2   * uvt[1];
+      uvt[2]  = 2.0 * tilt * uvt[1];
+      uvt[3]  = 2.0 * tilt *uvt[1] * x;              
+      uvt[4]  = 2.0 * (y + tilt * z) * uvt[1];
+      
+      Double_t error = 2.0 * uvt[1];
+      if (terror) {
+        error *= cseed[iLayer].fSigmaY;
+      }
+      else {
+        error *= 0.2; //Default error
+      }
+//                     printf("\tadd point :\n");
+//                     for(int i=0; i<5; i++) printf("%f ", uvt[i]);
+//                     printf("\n");
+      fitterT2.AddPoint(uvt,uvt[4],error);
+      npointsT++;
+
+    }
+
+  }
+  fitterT2.Eval();
+  Double_t rpolz0 = fitterT2.GetParameter(3);
+  Double_t rpolz1 = fitterT2.GetParameter(4);      
+
+  //
+  // Linear fitter  - not possible to make boundaries
+  // non accept non possible z and dzdx combination
+  //       
+  Bool_t acceptablez = kTRUE;
+  for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+    if (cseed[iLayer].IsOK()) {
+      Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].fX0 - xref2);
+      if (TMath::Abs(cseed[iLayer].fZProb - zT2) > cseed[iLayer].fPadLength * 0.5 + 1.0) {
+       acceptablez = kFALSE;
+      }
+    }
+  }
+  if (!acceptablez) {
+    Double_t zmf  = cseed[2].fZref[0] + cseed[2].fZref[1] * (xref2 - cseed[2].fX0);
+    Double_t dzmf = (cseed[2].fZref[1] + cseed[3].fZref[1]) * 0.5;
+    fitterT2.FixParameter(3,zmf);
+    fitterT2.FixParameter(4,dzmf);
+    fitterT2.Eval();
+    fitterT2.ReleaseParameter(3);
+    fitterT2.ReleaseParameter(4);
+    rpolz0 = fitterT2.GetParameter(3);
+    rpolz1 = fitterT2.GetParameter(4);
+  }
+  
+  Double_t chi2TR = fitterT2.GetChisquare() / Float_t(npointsT);  
+  Double_t params[3];
+  params[0] =  fitterT2.GetParameter(0);
+  params[1] =  fitterT2.GetParameter(1);
+  params[2] =  fitterT2.GetParameter(2);           
+  Double_t curvature =  1.0 + params[1] * params[1] - params[2] * params[0];
+
+  for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+
+    Double_t  x  = cseed[iLayer].fX0;
+    Double_t  y  = 0;
+    Double_t  dy = 0;
+    Double_t  z  = 0;
+    Double_t  dz = 0;
+
+    // y
+    Double_t res2 = (x * params[0] + params[1]);
+    res2 *= res2;
+    res2  = 1.0 - params[2]*params[0] + params[1]*params[1] - res2;
+    if (res2 >= 0) {
+      res2 = TMath::Sqrt(res2);
+      y    = (1.0 - res2) / params[0];
+    }
+
+    //dy
+    Double_t x0 = -params[1] / params[0];
+    if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) {
+      Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1); 
+      if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) {
+       Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0));
+       if (params[0] < 0) res *= -1.0;
+       dy = res;
+      }
+    }
+    z  = rpolz0 + rpolz1 * (x - xref2);
+    dz = rpolz1;
+    cseed[iLayer].fYref[0] = y;
+    cseed[iLayer].fYref[1] = dy;
+    cseed[iLayer].fZref[0] = z;
+    cseed[iLayer].fZref[1] = dz;
+    cseed[iLayer].fC       = curvature;
+    
+  }
+
+  return chi2TR;
+
+}
+
+//___________________________________________________________________
+void AliTRDseedV1::Print()
+{
+  //
+  // Printing the seedstatus
+  //
+
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       Int_t nTimeBins = cal->GetNumberOfTimeBins();
+       
+       printf("Seed status :\n");
+       printf("  fTilt      = %f\n", fTilt);
+       printf("  fPadLength = %f\n", fPadLength);
+       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"
+                , ic
+                , fX[ic]
+                , fY[ic]
+                , fZ[ic]
+                , fIndexes[ic]
+                , ((Int_t) fClusters[ic])
+                , isUsable);
+        }
+
+       printf("  fYref[0] =%f fYref[1] =%f\n", fYref[0], fYref[1]);
+       printf("  fZref[0] =%f fZref[1] =%f\n", fZref[0], fZref[1]);
+       printf("  fYfit[0] =%f fYfit[1] =%f\n", fYfit[0], fYfit[1]);
+       printf("  fYfitR[0]=%f fYfitR[1]=%f\n", fYfitR[0], fYfitR[1]);
+       printf("  fZfit[0] =%f fZfit[1] =%f\n", fZfit[0], fZfit[1]);
+       printf("  fZfitR[0]=%f fZfitR[1]=%f\n", fZfitR[0], fZfitR[1]);
+       printf("  fSigmaY =%f\n", fSigmaY);
+       printf("  fSigmaY2=%f\n", fSigmaY2);            
+       printf("  fMeanz  =%f\n", fMeanz);
+       printf("  fZProb  =%f\n", fZProb);
+       printf("  fLabels[0]=%d fLabels[1]=%d\n", fLabels[0], fLabels[1]);
+       printf("  fN      =%d\n", fN);
+       printf("  fN2     =%d (>8 isOK)\n",fN2);
+       printf("  fNUsed  =%d\n", fNUsed);
+       printf("  fFreq   =%d\n", fFreq);
+       printf("  fNChange=%d\n",  fNChange);
+       printf("  fMPads  =%f\n", fMPads);
+       
+       printf("  fC      =%f\n", fC);        
+       printf("  fCC     =%f\n",fCC);      
+       printf("  fChi2   =%f\n", fChi2);  
+       printf("  fChi2Z  =%f\n", fChi2Z);
+
+}
diff --git a/TRD/AliTRDseedV1.h b/TRD/AliTRDseedV1.h
new file mode 100644 (file)
index 0000000..9b6b6e2
--- /dev/null
@@ -0,0 +1,100 @@
+#ifndef ALITRDSEEDV1_H
+#define ALITRDSEEDV1_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+//  The TRD track seed                                                    //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALITRDSEED_H
+#include "AliTRDseed.h"
+#endif
+
+#ifndef ALIRIEMAN_H
+#include "AliRieman.h"
+#endif
+
+class TTreeSRedirector;
+
+class AliRieman;
+
+class AliTRDstackLayer;
+class AliTRDcluster;
+class AliTRDrecoParam;
+
+class AliTRDseedV1 : public AliTRDseed
+{
+
+ public:
+
+       AliTRDseedV1(Int_t layer = -1, AliTRDrecoParam *p=0x0);
+       ~AliTRDseedV1();
+       AliTRDseedV1(const AliTRDseedV1 &ref, Bool_t owner=kFALSE);
+       AliTRDseedV1& operator=(const AliTRDseedV1 &);
+
+       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);
+       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();
+
+              void    SetLayer(Int_t l)                      { fLayer     = l;   }
+              void    SetNTimeBins(Int_t nTB)                { fTimeBins  = nTB; }
+              void    SetRecoParam(AliTRDrecoParam *p)       { fRecoParam = p;   }
+
+ protected:
+
+       void Copy(TObject &ref) const;
+
+ 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 
+
+       ClassDef(AliTRDseedV1, 1)    //  New TRD seed 
+
+};
+
+//____________________________________________________________
+inline Float_t AliTRDseedV1::GetChi2Z(const Float_t z) const
+{
+       Float_t z1  = (z == 0.) ? fMeanz : z;
+       Float_t chi = fZref[0] - z1;
+       return chi*chi;
+}
+
+//____________________________________________________________
+inline Float_t AliTRDseedV1::GetChi2Y(const Float_t y) const
+{
+       Float_t y1  = (y == 0.) ? fYfitR[0] : y;
+       Float_t chi = fYref[0] - y1;
+       return chi*chi;
+}
+
+//____________________________________________________________
+inline void AliTRDseedV1::Update(const AliRieman *rieman)
+{
+       fZref[0] = rieman->GetZat(fX0);
+       fZref[1] = rieman->GetDZat(fX0);
+       fYref[0] = rieman->GetYat(fX0);
+       fYref[1] = rieman->GetDYat(fX0);
+}
+
+#endif
+
diff --git a/TRD/AliTRDstackLayer.cxx b/TRD/AliTRDstackLayer.cxx
new file mode 100644 (file)
index 0000000..65863ed
--- /dev/null
@@ -0,0 +1,550 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Organization of clusters at the level of 1 TRD chamber.                  //
+//  The data structure is used for tracking at the stack level.              //
+//                                                                           //
+//  Functionalities:                                                         //
+//   1. cluster organization and sorting                                     //
+//   2. fast data navigation                                                 //
+//                                                                           //
+//  Authors:                                                                 //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                        //
+//    Markus Fasel <M.Fasel@gsi.de>                                          //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include <TObject.h>
+#include <TROOT.h>
+#include <TMath.h>
+#include <TStopwatch.h>
+#include <TTreeStream.h>
+
+#include "AliLog.h"
+#include "AliTRDcluster.h"
+
+#include "AliTRDstackLayer.h"
+#include "AliTRDrecoParam.h"
+#include "AliTRDReconstructor.h"
+
+#define DEBUG
+
+ClassImp(AliTRDstackLayer)
+
+//_____________________________________________________________________________
+AliTRDstackLayer::AliTRDstackLayer(Double_t z0, Double_t zLength
+                                 , UChar_t stackNr, AliTRDrecoParam *p)
+  :AliTRDpropagationLayer()
+  ,fOwner(kFALSE)
+  ,fStackNr(stackNr)
+  ,fNRows(kMaxRows)
+  ,fZ0(z0)
+  ,fZLength(zLength)
+  ,fRecoParam(p)
+  ,fDebugStream(0x0)
+{
+  //
+  // Default constructor (Only provided to use AliTRDstackLayer with arrays)
+  //
+
+       for(int i=0; i<kMaxRows; i++) fPositions[i] = 0;
+}
+
+//_____________________________________________________________________________
+AliTRDstackLayer::AliTRDstackLayer(const AliTRDpropagationLayer &layer, Double_t
+z0, Double_t zLength, UChar_t stackNr, AliTRDrecoParam *p):
+       AliTRDpropagationLayer(layer)
+       ,fOwner(kFALSE)
+       ,fStackNr(stackNr)
+       ,fNRows(kMaxRows)
+       ,fZ0(z0)
+       ,fZLength(zLength)
+       ,fRecoParam(p)
+       ,fDebugStream(0x0)
+{
+// Standard constructor.
+// Initialize also the underlying AliTRDpropagationLayer using the copy constructor.
+
+       for(int i=0; i<kMaxRows; i++) fPositions[i] = 0;
+}
+
+//_____________________________________________________________________________
+AliTRDstackLayer::AliTRDstackLayer(const AliTRDpropagationLayer &layer):
+       AliTRDpropagationLayer(layer)
+       ,fOwner(kFALSE)
+       ,fStackNr(0)
+       ,fNRows(kMaxRows)
+       ,fZ0(0)
+       ,fZLength(0)
+       ,fRecoParam(0x0)
+       ,fDebugStream(0x0)
+{
+// Standard constructor using only AliTRDpropagationLayer.
+       
+       for(int i=0; i<kMaxRows; i++) fPositions[i] = 0;
+}
+
+//_____________________________________________________________________________
+AliTRDstackLayer::AliTRDstackLayer(const AliTRDstackLayer &layer):
+       AliTRDpropagationLayer(layer)
+       ,fOwner(layer.fOwner)
+       ,fStackNr(layer.fStackNr)
+       ,fNRows(layer.fNRows)
+       ,fZ0(layer.fZ0)
+       ,fZLength(layer.fZLength)
+       ,fRecoParam(layer.fRecoParam)
+       ,fDebugStream(layer.fDebugStream)
+{
+// Copy Constructor (performs a deep copy)
+       
+       for(Int_t i = 0; i < kMaxRows; i++) fPositions[i] = layer.fPositions[i];
+//     BuildIndices();
+}
+
+//_____________________________________________________________________________
+AliTRDstackLayer &AliTRDstackLayer::operator=(const AliTRDpropagationLayer &layer)
+{
+// Assignment operator from an AliTRDpropagationLayer
+
+       if (this != &layer) layer.Copy(*this);
+       return *this;
+}
+
+//_____________________________________________________________________________
+AliTRDstackLayer &AliTRDstackLayer::operator=(const AliTRDstackLayer &layer)
+{
+// Assignment operator
+
+       if (this != &layer) layer.Copy(*this);
+  return *this;
+}
+
+//_____________________________________________________________________________
+void AliTRDstackLayer::Copy(TObject &o) const
+{
+// Copy method. Performs a deep copy of all data from this object to object o.
+       
+       AliTRDstackLayer &layer = (AliTRDstackLayer &)o;
+       layer.fZ0          = fZ0;
+       layer.fOwner       = kFALSE;
+       layer.fNRows       = fNRows;
+       layer.fZLength     = fZLength;
+       layer.fStackNr     = fStackNr;
+       layer.fDebugStream = fDebugStream;
+       layer.fRecoParam   = fRecoParam;
+       
+       AliTRDpropagationLayer::Copy(layer); // copies everything into layer
+       for(UChar_t i = 0; i < kMaxRows; i++) layer.fPositions[i] = 0;
+//     layer.BuildIndices();
+}
+
+//_____________________________________________________________________________
+AliTRDstackLayer::~AliTRDstackLayer()
+{
+// Destructor
+       if(fOwner) for(int ic=0; ic<fN; ic++) delete fClusters[ic];
+}
+
+//_____________________________________________________________________________
+void AliTRDstackLayer::SetRange(const Float_t z0, const Float_t zLength)
+{
+// Sets the range in z-direction
+//
+// Parameters:
+//   z0      : starting position of layer in the z direction
+//   zLength : length of layer in the z direction 
+
+       fZ0 = (z0 <= z0 + zLength) ? z0 : z0 + zLength;
+       fZLength = TMath::Abs(zLength);
+}
+
+//_____________________________________________________________________________
+void AliTRDstackLayer::BuildIndices(Int_t iter)
+{
+// Rearrangement of the clusters belonging to the propagation layer for the stack.
+//
+// Detailed description
+//
+// The array indices of all clusters in one PropagationLayer are stored in
+// array. The array is divided into several bins.
+// The clusters are sorted in increasing order of their y coordinate.
+//
+// Sorting algorithm: TreeSearch
+//
+
+       if(!fN) return;
+
+       // Select clusters that belong to the Stack
+       Int_t nClStack = 0;                                     // Internal counter
+       for(Int_t i = 0; i < fN; i++){
+               Double_t zval = fClusters[i]->GetZ();
+               if(zval < fZ0 || zval > fZ0 + fZLength || fClusters[i]->IsUsed()){
+                       fClusters[i] = 0x0;
+                       fIndex[i] = 9999;
+               } else nClStack++;
+       }
+       if(nClStack > kMaxClustersLayer) AliWarning(Form("Number of clusters in stack %d exceed buffer size %d", nClStack, kMaxClustersLayer));
+       
+       // Nothing in this Stack 
+       if(!nClStack){
+               delete fClusters;
+               delete fIndex;
+               fClusters = 0x0;
+               fIndex = 0x0;
+               fN = 0;
+               memset(fPositions, 0, sizeof(UChar_t) * 16);
+               return;
+       }
+       
+       // Make a copy
+       AliTRDcluster *helpCL[kMaxClustersLayer];
+       Int_t helpInd[kMaxClustersLayer];
+       nClStack = 0;
+       for(Int_t i = 0; i < TMath::Min(fN, kMaxClustersLayer); i++){
+               if(!fClusters[i]) continue;
+               helpCL[nClStack]  = fClusters[i];
+               helpInd[nClStack] = fIndex[i];
+               fClusters[i]      = 0x0;
+               fIndex[i]         = 9999;
+               nClStack++;
+       }
+       
+       // do clusters arrangement
+       fN =  nClStack;
+       nClStack = 0;
+       memset(fPositions,0,sizeof(UChar_t)*16);                                // Reset Positions array
+       for(Int_t i = 0; i < fN; i++){
+               // boundarie check
+               AliTRDcluster *cl = helpCL[i];
+               Double_t zval = cl->GetZ();
+               UChar_t TreeIndex = (UChar_t)(TMath::Abs(fZ0 - zval)/fZLength * fNRows);
+               if(TreeIndex > fNRows - 1) TreeIndex = fNRows - 1;
+               // Insert Leaf
+               Int_t pos = FindYPosition(cl->GetY(), TreeIndex, i);
+               if(pos == -1){          // zbin is empty;
+                       Int_t upper = (TreeIndex == fNRows - 1) ? nClStack : fPositions[TreeIndex + 1];
+                       memmove(fClusters + upper + 1, fClusters + upper, (sizeof(AliTRDcluster *))*(nClStack-upper));
+                       memmove(fIndex + upper + 1, fIndex + upper, (sizeof(UInt_t))*(nClStack-upper));
+                       fClusters[upper] = cl;
+                       fIndex[upper] = helpInd[i]; 
+                       // Move All pointer one position back
+                       for(UChar_t j = TreeIndex + 1; j < fNRows; j++) fPositions[j]++;
+                       nClStack++;
+               } else {                // zbin not empty
+                       memmove(fClusters + pos + 2, fClusters + pos+1, (sizeof(AliTRDcluster *))*(nClStack-(pos+1)));
+                       memmove(fIndex + pos + 2, fIndex + pos+1, (sizeof(UInt_t))*(nClStack-(pos+1)));
+                       fClusters[pos + 1] = cl;        //fIndex[i];
+                       fIndex[pos + 1] = helpInd[i];
+                       // Move All pointer one position back
+                       for(UChar_t j = TreeIndex + 1; j < fNRows; j++) fPositions[j]++;        
+                       nClStack++;
+               }
+
+               
+               // Debug Streaming
+#ifdef DEBUG
+               if(fDebugStream && AliTRDReconstructor::StreamLevel() >= 3){
+                       TTreeSRedirector &cstream = *fDebugStream;
+                       cstream << "BuildIndices"
+                       << "StackNr="  << fStackNr
+                       << "SectorNr=" << fSec
+                       << "Iter="     << iter
+                       << "C.="       << cl
+                       << "TreeIdx="  << TreeIndex
+                       << "\n";
+               }
+#endif 
+       }
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDstackLayer::FindYPosition(Double_t y, UChar_t z, Int_t nClusters) const
+{
+//
+// Tree search Algorithm to find the nearest left cluster for a given
+// y-position in a certain z-bin (in fact AVL-tree). 
+// Making use of the fact that clusters are sorted in y-direction.
+//
+// Parameters:
+//   y : y position of the reference point in tracking coordinates
+//   z : z reference bin.
+//   nClusters : 
+//
+// Output :
+// Index of the nearest left cluster in the StackLayer indexing (-1 if no clusters are found)
+//
+
+       Int_t start = fPositions[z];            // starting Position of the bin
+       Int_t upper = (Int_t)((z != fNRows - 1) ? fPositions[z+1] : nClusters); // ending Position of the bin 
+       Int_t end = upper - 1; // ending Position of the bin 
+       if(end < start) return -1; // Bin is empty
+       Int_t middle = static_cast<Int_t>((start + end)/2);
+       // 1st Part: climb down the tree: get the next cluster BEFORE ypos
+       while(start + 1 < end){
+               if(y >= fClusters[middle]->GetY()) start = middle;
+               else end = middle;
+               middle = static_cast<Int_t>((start + end)/2);
+       }
+       if(y > fClusters[end]->GetY()) return end;
+       return start;
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDstackLayer::FindNearestYCluster(Double_t y, UChar_t z) const
+{
+//
+// Tree search Algorithm to find the nearest cluster for a given
+// y-position in a certain z-bin (in fact AVL-tree). 
+// Making use of the fact that clusters are sorted in y-direction.
+//
+// Parameters:
+//   y : y position of the reference point in tracking coordinates
+//   z : z reference bin.
+//
+// Output 
+// Index of the nearest cluster in the StackLayer indexing (-1 if no clusters are found)
+//
+
+       Int_t position = FindYPosition(y, z, fN);
+       if(position == -1) return position; // bin empty
+       // FindYPosition always returns the left Neighbor. We don't know if the left or the right Neighbor is nearest
+       // to the Reference y-position, so test both
+       Int_t upper = (Int_t)((z < fNRows-1) ? fPositions[z+1] : fN); // ending Position of the bin
+       if((position + 1) < (upper)){
+               if(TMath::Abs(y - fClusters[position + 1]->GetY()) < TMath::Abs(y - fClusters[position]->GetY())) return position + 1;
+               else return position;
+       }
+       return position;
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDstackLayer::SearchNearestCluster(Double_t y, Double_t z, Double_t maxroady, Double_t maxroadz) const
+{
+//
+// Finds the nearest cluster from a given point in a defined range.
+// Distance is determined in a 2D space by the 2-Norm.
+//
+// Parameters :
+//   y : y position of the reference point in tracking coordinates
+//   z : z reference bin.
+//   maxroady : maximum searching distance in y direction
+//   maxroadz : maximum searching distance in z direction
+//
+// Output 
+// Index of the nearest cluster in the StackLayer indexing (-1 if no cluster is found).
+// Cluster can be accessed with the operator[] or GetCluster(Int_t index)
+//
+// Detail description
+//
+// The following steps are perfomed:
+// 1. Get the expected z bins inside maxroadz.
+// 2. For each z bin find nearest y cluster.
+// 3. Select best candidate
+//
+       Int_t   index   = -1;
+       // initial minimal distance will be represented as ellipse: semi-major = z-direction
+       // later 2-Norm will be used  
+//     Float_t nExcentricity = TMath::Sqrt(maxroadz*maxroadz - maxroad*maxroad)/maxroadz;
+       Float_t mindist = maxroadz;
+       
+       // not very nice but unfortunately neccessarry: we have ho check the neighbors in both directions (+ and -) too. How 
+       // much neighbors depends on the Quotient maxroadz/fZLength   
+       UChar_t maxRows = 3;
+       UChar_t zpos[kMaxRows];
+  // Float_t mindist = TMath::Sqrt(maxroad*maxroad + maxroadz*maxroadz);
+//     UChar_t myZbin = FindTreePosition(z, fZ0 + fZLength/2, fZLength/4, 8, 8, kFALSE);
+       UChar_t myZbin = (UChar_t)(TMath::Abs(fZ0 - z)/fZLength * fNRows);
+       if(z < fZ0) myZbin = 0;
+       if(z > fZ0 + fZLength) myZbin = fNRows - 1;
+
+       UChar_t nNeighbors = 0;
+       for(UChar_t i = 0; i < maxRows; i++){
+               if((myZbin - 1 + i) < 0) continue;
+               if((myZbin - 1 + i) > fNRows - 1) break;
+               zpos[nNeighbors] = myZbin - 1 + i;
+               nNeighbors++;
+       }
+       Float_t ycl = 0, zcl = 0;
+       for(UChar_t neighbor = 0; neighbor < nNeighbors; neighbor++){   // Always test the neighbors too
+               Int_t pos = FindNearestYCluster(y,zpos[neighbor]);
+               if(pos == -1) continue; // No cluster in bin
+               AliTRDcluster *c = (AliTRDcluster *) (fClusters[pos]);
+               if(c->IsUsed()) continue;               // we are only interested in unused clusters
+               ycl = c->GetY();
+               // Too far away in y-direction (Prearrangement)
+               if (TMath::Abs(ycl - y) > maxroady) continue;
+
+               zcl = c->GetZ();
+               // Too far away in z-Direction
+               // (Prearrangement since we have not so many bins to test)
+               if (TMath::Abs(zcl - z) > maxroadz) continue;
+               
+               Float_t dist;           // distance defined as 2-Norm   
+               // if we havent found a Particle that is in the ellipse around (y,z) with maxroad as semi-minor and
+               // maxroadz as semi-major, we take the radius of the ellipse concerning the cluster as mindist, later we 
+               // take the 2-Norm when we found a cluster inside the ellipse (The value 10000 is taken because it is surely
+               // large enough to be usable as an indicator whether we have found a nearer cluster or not)
+//             if(mindist > 10000.){
+//                     Float_t phi = ((zcl - z) == 0) ? TMath::Pi()/2 : TMath::ATan((ycl - y)/(zcl - z));
+//                     mindist = maxroad/TMath::Sqrt(1 - nExcentricity*nExcentricity * (TMath::Cos(phi))*(TMath::Cos(phi)));
+//             }
+               dist = TMath::Max(TMath::Abs(y-ycl),TMath::Abs(z-zcl)); // infinity Norm
+//             dist = TMath::Sqrt((ycl - y)*(ycl - y) + (zcl - z)*(zcl - z));
+               if((Int_t)(dist * 100000) < (Int_t)(mindist * 100000)){
+               //if((dist = TMath::Sqrt((ycl - y)*(ycl - y) + (zcl - z)*(zcl - z))) < mindist){
+                       mindist = dist;
+                       index   = pos;
+               }       
+       }
+       // This is the Array Position in fIndex2D of the Nearest cluster: if a
+       // cluster is called, then the function has to retrieve the Information
+       // which is Stored in the Array called, the function
+       return index;
+}
+
+//_____________________________________________________________________________
+void AliTRDstackLayer::BuildCond(AliTRDcluster *cl, Double_t *cond, UChar_t Layer, Double_t theta, Double_t phi)
+{
+// Helper function to calculate the area where to expect a cluster in THIS
+// layer. 
+//
+// Parameters :
+//   cl    : 
+//   cond  :
+//   Layer : 
+//   theta : 
+//   phi   : 
+//
+// Detail description
+//
+// Helper function to calculate the area where to expect a cluster in THIS
+// layer. by using the information of a former cluster in another layer
+// and the angle in theta- and phi-direction between layer 0 and layer 3.
+// If the layer is zero, initial conditions are calculated. Otherwise a
+// linear interpolation is performed. 
+//Begin_Html
+//<img src="gif/build_cond.gif">
+//End_Html
+//
+
+       if(!fRecoParam){
+               AliError("Reconstruction parameters not initialized.");
+               return;
+       }
+       
+       if(Layer == 0){
+               cond[0] = cl->GetY();                   // center: y-Direction
+               cond[1] = cl->GetZ();                   // center: z-Direction
+               cond[2] = fRecoParam->GetMaxPhi()   * (cl->GetX() - GetX()) + 1.0;                      // deviation: y-Direction
+               cond[3] = fRecoParam->GetMaxTheta() * (cl->GetX() - GetX()) + 1.0;                      // deviation: z-Direction
+       } else {
+               cond[0] = cl->GetY() + phi   * (GetX() - cl->GetX()); 
+               cond[1] = cl->GetZ() + theta * (GetX() - cl->GetX());
+               cond[2] = fRecoParam->GetRoad0y() + phi;
+               cond[3] = fRecoParam->GetRoad0z();
+       }
+}
+
+//_____________________________________________________________________________
+void AliTRDstackLayer::GetClusters(Double_t *cond, Int_t *index, Int_t& ncl, Int_t BufferSize)
+{
+// Finds all clusters situated in this layer inside a rectangle  given by the center an ranges.
+//
+// Parameters :
+//   cond  :
+//   index : 
+//   ncl :
+//   BufferSize   :
+//
+// Output :
+//
+// Detail description
+//
+// Function returs an array containing the indices in the stacklayer of
+// the clusters found an  the number of found clusters in the stacklayer
+
+       ncl = 0;
+       memset(index, 0, BufferSize*sizeof(Int_t));
+       if(fN == 0) return;
+       
+       //Boundary checks
+       Double_t zvals[2];
+       zvals[0] = ((cond[1] - cond[3]) < fZ0) ? fZ0 : (cond[1] - cond[3]);
+       zvals[1] = ((cond[1] + cond[3]) < fZ0 + fZLength) ? (cond[1] + cond[3]) : fZ0 + fZLength;
+
+       UChar_t zlo = (fZ0>zvals[0]) ? 0 : (UChar_t)(TMath::Abs(fZ0 - zvals[0])/fZLength * fNRows);
+       UChar_t zhi = (fZ0+fZLength<zvals[1]) ? fNRows - 1 : (UChar_t)(TMath::Abs(fZ0 - zvals[1])/fZLength * fNRows);
+
+       //Preordering in Direction z saves a lot of loops (boundary checked)
+       for(UChar_t z = zlo; z <= zhi; z++){
+               UInt_t upper = (z != fNRows-1) ? fPositions[z+1] : fN;
+               for(Int_t y = fPositions[z]; y < (Int_t)upper; y++){
+                       if(ncl == BufferSize){
+                               AliWarning("Buffer size riched. Some clusters may be lost.");
+                               return; //Buffer filled
+                       }
+                       if(fClusters[y]->GetY() > (cond[0] + cond[2])) break;                   // Abbortion conditions!!!
+                       if(fClusters[y]->GetY() < (cond[0] - cond[2])) continue;        // Too small
+                       if(((Int_t)((fClusters[y]->GetZ())*1000) < (Int_t)(zvals[0]*1000)) || ((Int_t)((fClusters[y]->GetZ())*1000) > (Int_t)(zvals[1]*1000))) continue;
+                       index[ncl] = y;
+                       ncl++;
+               }
+       }
+       if(ncl>fN) AliError(Form("Clusters found %d > %d (clusters in layer)", ncl, fN));
+}
+
+//_____________________________________________________________________________
+AliTRDcluster *AliTRDstackLayer::GetNearestCluster(Double_t *cond)
+{
+// Function returning a pointer to the nearest cluster (nullpointer if not successfull).
+//
+// Parameters :
+//   cond  :
+//
+// Output :
+//   pointer to the nearest cluster (nullpointer if not successfull).
+// 
+// Detail description
+//
+// returns a pointer to the nearest cluster (nullpointer if not
+// successfull) by the help of the method FindNearestCluster
+       
+       
+       Double_t maxroad  = fRecoParam->GetRoad2y();
+       Double_t maxroadz = fRecoParam->GetRoad2z();
+       
+       Int_t index = SearchNearestCluster(cond[0],cond[1],maxroad,maxroadz);
+       AliTRDcluster *returnCluster = 0x0;
+       if(index != -1) returnCluster = (AliTRDcluster *) fClusters[index];
+       return returnCluster;
+}
+
+//_____________________________________________________________________________
+void AliTRDstackLayer::PrintClusters() const
+{
+// Prints the position of each cluster in the stacklayer on the stdout
+//
+       printf("fDebugStream = %p\n", fDebugStream);
+       printf("fRecoParam   = %p\n", fRecoParam);
+       
+       for(Int_t i = 0; i < fN; i++){
+               printf("AliTRDstackLayer: index=%i, Cluster: X = %3.3f, Y = %3.3f, Z = %3.3f\n", i,  fClusters[i]->GetX(),fClusters[i]->GetY(),fClusters[i]->GetZ());
+               if(fClusters[i]->IsUsed()) printf("cluster allready used. rejected in search algorithm\n");
+       }
+}
diff --git a/TRD/AliTRDstackLayer.h b/TRD/AliTRDstackLayer.h
new file mode 100644 (file)
index 0000000..1717d9f
--- /dev/null
@@ -0,0 +1,87 @@
+#ifndef ALITRDSTACKLAYER_H
+#define ALITRDSTACKLAYER_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+//  A TRD layer in a single stack                                         //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALITRDPROPAGATIONLAYER_H
+#include "AliTRDpropagationLayer.h"
+#endif
+
+#ifndef ALITRDCLUSTER_H
+#include "AliTRDcluster.h"
+#endif
+
+class AliTRDrecoParam;
+
+class AliTRDstackLayer : public AliTRDpropagationLayer
+{
+  public:
+       enum{
+               kMaxClustersLayer = 150,
+               kMaxRows = 16
+       };
+
+       AliTRDstackLayer(Double_t z0 = 0., Double_t zLength = 0., UChar_t stackNr = 0
+                       , AliTRDrecoParam *p=0x0);
+       AliTRDstackLayer(const AliTRDpropagationLayer &layer, Double_t z0
+                       , Double_t zLength, UChar_t stackNr, AliTRDrecoParam *p = 0x0);
+       AliTRDstackLayer(const AliTRDpropagationLayer &layer);
+       AliTRDstackLayer(const AliTRDstackLayer &layer);
+       ~AliTRDstackLayer();
+       AliTRDstackLayer   &operator=(const AliTRDpropagationLayer &myLayer);
+       AliTRDstackLayer   &operator=(const AliTRDstackLayer &myLayer);
+       AliTRDcluster      *operator[](const Int_t i) {
+               return ((i < fN) && (i >= 0)) ? fClusters[i] : 0x0;
+       }
+
+       void           BuildIndices(Int_t iter = 0);
+       void           BuildCond(AliTRDcluster *cl, Double_t *cond, UChar_t Layer, Double_t theta=0., Double_t phi=0.);
+       AliTRDcluster* GetCluster(Int_t index) const {return index < fN ? fClusters[index] : 0x0;}
+       Int_t          GetGlobalIndex(const Int_t index) const {return ((index < fN) && (index >= 0)) ? fIndex[index] : 0; }
+       void           GetClusters(Double_t *cond, Int_t *index, Int_t& ncl, Int_t BufferSize = kMaxClustersLayer);
+       AliTRDcluster* GetNearestCluster(Double_t *cond);
+
+       Double_t       GetZ0()                           const { return fZ0; }
+       Double_t       GetDZ0()                          const { return fZLength; }
+       Int_t          GetNClusters()                    const { return fN; }
+       UInt_t         GetStackNr()                      const { return fStackNr; }
+       void           PrintClusters()                   const;
+       Int_t          SearchNearestCluster(const Double_t y, const Double_t z, const Double_t Roady, const Double_t Roadz) const;
+       void           SetRange(Float_t z0, Float_t zLength);
+       void           SetNRows(const Int_t nRows){ fNRows = nRows; }
+       void           SetStackNr(const UInt_t stackNr){ fStackNr = stackNr; }
+       void           SetOwner(Bool_t own = kTRUE) {fOwner = own;}
+       void           SetClustersArray(AliTRDcluster **cl, Int_t nClusters){fClusters = cl; fN = nClusters;}
+       void           SetIndexArray(UInt_t *indexArray){fIndex = indexArray;}
+       void           SetDebugStream(TTreeSRedirector *debug) {fDebugStream = debug;}
+       void           SetRecoParam(AliTRDrecoParam *p) {fRecoParam = p;}
+
+private:
+       void           Copy(TObject &o) const;
+       Int_t          FindYPosition(Double_t y, UChar_t z, Int_t nClusters) const;
+       Int_t          FindNearestYCluster(Double_t y, UChar_t z) const;
+
+private:
+       Bool_t            fOwner;               //  owner of the clusters
+       UChar_t           fStackNr;             //  stack number in supermodule
+       UChar_t           fNRows;               //  number of pad rows in the chamber
+       UChar_t           fPositions[kMaxRows]; //  starting index of clusters in pad row 
+       Double_t          fZ0;                  //  starting position of the layer in Z direction
+       Double_t          fZLength;             //  length of the layer in Z direction
+       AliTRDrecoParam  *fRecoParam;           //! reconstruction parameters
+       TTreeSRedirector *fDebugStream;         //! debug streamer
+       
+       ClassDef(AliTRDstackLayer, 1)           //  stack propagation layer
+
+};
+#endif // ALITRDSTACKLAYER_H_
+
index 9ea0edb..475083b 100644 (file)
@@ -15,9 +15,6 @@
 
 /* $Id$ */
 
-//#include <Riostream.h>
-
-//#include <TMath.h>
 #include <TVector2.h>
 
 #include "AliTracker.h"
@@ -26,7 +23,6 @@
 #include "AliTRDgeometry.h" 
 #include "AliTRDcluster.h" 
 #include "AliTRDtrack.h"
-//#include "AliTRDtracklet.h"
 #include "AliTRDcalibDB.h"
 #include "Cal/AliTRDCalPID.h"
 
@@ -506,8 +502,8 @@ void AliTRDtrack::CookdEdxTimBin(const Int_t/* tid*/)
     }
 
     tb     = cluster->GetLocalTimeBin();
-    if ((tb == 0) || (tb >= ntb)) {
-      AliWarning(Form("time bin 0 or > %d in cluster %d", ntb, iClus));
+    if ((tb < 0) || (tb >= ntb)) {
+      AliWarning(Form("time bin < 0 or > %d in cluster %d", ntb, iClus));
       AliInfo(Form("dQ/dl %f", fdQdl[iClus]));
       continue;
     }
index 758d31e..033743f 100644 (file)
@@ -399,65 +399,6 @@ Int_t  AliTRDtracker::GetLastPlane(AliTRDtrack *track)
   return lastplane;
 
 }
-
-//_____________________________________________________________________________
-Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *event)
-{
-  //
-  // Finds tracks within the TRD. The ESD event is expected to contain seeds 
-  // at the outer part of the TRD. The seeds
-  // are found within the TRD if fAddTRDseeds is TRUE. 
-  // The tracks are propagated to the innermost time bin 
-  // of the TRD and the ESD event is updated
-  //
-
-  Int_t   timeBins = fTrSec[0]->GetNumberOfTimeBins();
-  Float_t foundMin = fgkMinClustersInTrack * timeBins; 
-  Int_t   nseed    = 0;
-  Int_t   found    = 0;
-  //Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
-
-  Int_t n = event->GetNumberOfTracks();
-  for (Int_t i = 0; i < n; i++) {
-
-    AliESDtrack *seed = event->GetTrack(i);
-    ULong_t status = seed->GetStatus();
-    if ((status & AliESDtrack::kTRDout) == 0) {
-      continue;
-    }
-    if ((status & AliESDtrack::kTRDin)  != 0) {
-      continue;
-    }
-    nseed++;
-    
-    AliTRDtrack *seed2 = new AliTRDtrack(*seed);
-    //seed2->ResetCovariance(); 
-    AliTRDtrack *pt    = new AliTRDtrack(*seed2,seed2->GetAlpha());
-    AliTRDtrack &t     = *pt; 
-    FollowProlongation(t); 
-    if (t.GetNumberOfClusters() >= foundMin) {
-      UseClusters(&t);
-      CookLabel(pt,1 - fgkLabelFraction);
-      //t.CookdEdx();
-    }
-    found++;
-
-    Double_t xTPC = 250.0;
-    if (PropagateToX(t,xTPC,fgkMaxStep)) {
-      seed->UpdateTrackParams(pt, AliESDtrack::kTRDin);
-    }  
-    delete seed2;
-    delete pt;
-
-  }
-
-  AliInfo(Form("Number of loaded seeds: %d",nseed));
-  AliInfo(Form("Number of found tracks from loaded seeds: %d",found));
-  AliInfo(Form("Total number of found tracks: %d",found));
-    
-  return 0;    
-
-}     
      
 //_____________________________________________________________________________
 Int_t AliTRDtracker::PropagateBack(AliESDEvent *event) 
@@ -471,259 +412,244 @@ Int_t AliTRDtracker::PropagateBack(AliESDEvent *event)
   // by the TPC tracker.   
   //  
 
-  Int_t   found    = 0;     // number of tracks found
-  Float_t foundMin = 20.0;
-
-  // Sort tracks according to covariance of local Y and Z
-  Int_t    nSeed   = event->GetNumberOfTracks();
-  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];      
-  }
-  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]);
-    fHBackfit->Fill(0);   // All seeds
-
-    // Check the seed status
-    ULong_t status = seed->GetStatus();
-    if ((status & AliESDtrack::kTPCout) == 0) {
-      fHBackfit->Fill(1); // TPC outer edge reached
-      continue;
-    }
-    if ((status & AliESDtrack::kTRDout) != 0) {
-      fHBackfit->Fill(2); // TRD outer edge reached (does this happen ?)
-      continue;
-    }
-
-    // Do the back prolongation
-    Int_t   lbl         = seed->GetLabel();
-    AliTRDtrack *track  = new AliTRDtrack(*seed);
-    track->SetSeedLabel(lbl);
-    seed->UpdateTrackParams(track,AliESDtrack::kTRDbackup); // Make backup
-    fNseeds++;
-    Float_t p4          = track->GetC();
-    Int_t   expectedClr = FollowBackProlongation(*track);
-    fHBackfit->Fill(3);   // Back prolongation done
-    fHX->Fill(track->GetX());
-    // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-    if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) || 
-        (track->Pt()                                     > 0.8)) {
-      
-      fHBackfit->Fill(4);
-      
-      // 
-      // Make backup for back propagation 
-      //
-      
-      Int_t foundClr = track->GetNumberOfClusters();
-      if (foundClr >= foundMin) {
-       track->CookdEdx(); 
-       track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
-       CookLabel(track,1 - fgkLabelFraction);
-       if (track->GetBackupTrack()) {
-          UseClusters(track->GetBackupTrack());
-       }
-
-        // 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;
+       Int_t   found    = 0;     // number of tracks found
+       Float_t foundMin = 20.0;
        
-        // Full gold track
-       if (track->GetChi2() / track->GetNumberOfClusters() < 5) {  
-         //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
-         if (track->GetBackupTrack()) {
-            seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
-           
-         }
-         isGold = kTRUE;
-         //fHBackfit->Fill()
-       }
-
-        // 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;
-       }
-
-       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;
-         }
+       Int_t    nSeed   = event->GetNumberOfTracks();
+       if(!nSeed){
+               // run stand alone tracking
+               if (AliTRDReconstructor::SeedingOn()) Clusters2Tracks(event);
+               return 0;
        }
-
-       if ((track->StatusForTOF() > 0) &&
-            (track->GetNCross() == 0) && 
-            (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected())  > 0.4)) {
-         //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
+       
+       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];
        }
-
-      }
-    }
-    /**/
-
-    /**/
-    // 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";
-      }
-    }
-    /**/
-
-    // Propagation to the TOF (I.Belikov)    
-    if (track->GetStop() == kFALSE) {
-      fHBackfit->Fill(5);
-
-      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) {
+       // Sort tracks according to covariance of local Y and Z
+       TMath::Sort(nSeed,quality,index,kFALSE);
        
-       fHBackfit->Fill(6);
-       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) {
+       // Backpropagate all seeds
+       for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
        
-       fHBackfit->Fill(7);
-       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())) {
-         fHBackfit->Fill(8);
-         delete track;
-         continue;
-       }
-      } 
-      else if (y < -ymax) {
-       if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
-         fHBackfit->Fill(9);
-         delete track;
-         continue;
-       }
-      }
-         
-      if (track->PropagateTo(xtof)) {
-       seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
-       fHBackfit->Fill(10);
-
-        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) {
-         fHBackfit->Fill(11);
-          found++;
+               // Get the seeds in sorted sequence
+               AliESDtrack *seed = event->GetTrack(index[iSeed]);
+               fHBackfit->Fill(0);   // All seeds
+       
+               // Check the seed status
+               ULong_t status = seed->GetStatus();
+               if ((status & AliESDtrack::kTPCout) == 0) {
+                       fHBackfit->Fill(1); // TPC outer edge reached
+                       continue;
+               }
+               if ((status & AliESDtrack::kTRDout) != 0) {
+                       fHBackfit->Fill(2); // TRD outer edge reached (does this happen ?)
+                       continue;
+               }
+       
+               // Do the back prolongation
+               Int_t   lbl         = seed->GetLabel();
+               AliTRDtrack *track  = new AliTRDtrack(*seed);
+               track->SetSeedLabel(lbl);
+               seed->UpdateTrackParams(track,AliESDtrack::kTRDbackup); // Make backup
+               fNseeds++;
+               Float_t p4          = track->GetC();
+               Int_t   expectedClr = FollowBackProlongation(*track);
+               fHBackfit->Fill(3);   // Back prolongation done
+               fHX->Fill(track->GetX());
+       
+               if ((TMath::Abs(track->GetC() - p4) / TMath::Abs(p4) < 0.2) ||
+                               (track->Pt()                                     > 0.8)) {
+                       
+                       fHBackfit->Fill(4);
+                       
+                       //
+                       // Make backup for back propagation
+                       //
+                       
+                       Int_t foundClr = track->GetNumberOfClusters();
+                       if (foundClr >= foundMin) {
+                               track->CookdEdx();
+                               track->CookdEdxTimBin(seed->GetID()); // A.Bercuci 25.07.07
+                               CookLabel(track,1 - fgkLabelFraction);
+                               if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack());
+                               
+                               
+                               // 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) {
+                                       //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
+                                       if (track->GetBackupTrack()) {
+                                               seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup);
+                                       }
+                                       isGold = kTRUE;
+                                       //fHBackfit->Fill()
+                               }
+       
+                               // 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;
+                               }
+                               
+                               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;
+                                       }
+                               }
+       
+                               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";
+                       }
+               }
+               /**/
+       
+               // Propagation to the TOF (I.Belikov)
+               if (track->GetStop() == kFALSE) {
+                       fHBackfit->Fill(5);
+       
+                       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) {
+                               fHBackfit->Fill(6);
+                               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) {
+                               fHBackfit->Fill(7);
+                               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())) {
+                                       fHBackfit->Fill(8);
+                                       delete track;
+                                       continue;
+                               }
+                       }
+                       else if (y < -ymax) {
+                               if (!track->Rotate(-AliTRDgeometry::GetAlpha())) {
+                                       fHBackfit->Fill(9);
+                                       delete track;
+                                       continue;
+                               }
+                       }
+                                       
+                       if (track->PropagateTo(xtof)) {
+                               seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
+                               fHBackfit->Fill(10);
+       
+                               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) {
+                                       fHBackfit->Fill(11);
+                                                               found++;
+                               }
+                       }
+               }
+               else {                  
+                       fHBackfit->Fill(12);
+                       
+                       if ((track->GetNumberOfClusters() >              15) &&
+                                       (track->GetNumberOfClusters() > 0.5*expectedClr)) {
+                               seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
+                               fHBackfit->Fill(13);
+       
+                               //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++;
+                       }
+               }
+       
+               seed->SetTRDQuality(track->StatusForTOF());
+               seed->SetTRDBudget(track->GetBudget(0));
+       
+               fHBackfit->Fill(14);
+               delete track;
        }
-      }
-
-    }
-    else {
-      
-      fHBackfit->Fill(12);
-      
-      if ((track->GetNumberOfClusters() >              15) &&
-          (track->GetNumberOfClusters() > 0.5*expectedClr)) {
        
-       seed->UpdateTrackParams(track,AliESDtrack::kTRDout);
-       fHBackfit->Fill(13);
-
-       //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++;
-      }
-
-    }
-
-    seed->SetTRDQuality(track->StatusForTOF());    
-    seed->SetTRDBudget(track->GetBudget(0));    
-  
-    fHBackfit->Fill(14);
-    delete track;
-  }
-
-  AliInfo(Form("Number of seeds: %d",fNseeds));
-  AliInfo(Form("Number of back propagated TRD tracks: %d",found));
-  
-  // New seeding
-  if (AliTRDReconstructor::SeedingOn()) {
-    MakeSeedsMI(3,5,event);
-  }
-
-  fSeeds->Clear(); 
-  fNseeds = 0;
-
-  delete [] index;
-  delete [] quality;
-  
-  SaveLogHists();
+       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;
+       
+       SaveLogHists();
   
   return 0;
-
 }
 
 //_____________________________________________________________________________
@@ -900,7 +826,7 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t)
     // Get material budget
     Double_t xyz0[3];
     Double_t xyz1[3];
-    Double_t param[7];
+    Double_t param[7]; 
     Double_t x;
     Double_t y;
     Double_t z;
@@ -1127,13 +1053,11 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t)
       fHClSearch->Fill(hb+6);
       continue;
     }
-    //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
     //
     // Propagate and update track
     //
     for (Int_t itime = GetTimeBinsPerPlane()-1; itime >= 0; itime--) {
-
       Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
       expectedNumberOfClusters++;       
       t.SetNExpected(t.GetNExpected() + 1);
@@ -1147,50 +1071,48 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t)
       x = timeBin.GetX();
 
       if (timeBin) {   
-
-       if (clusters[ilayer] > 0) {
-         index = clusters[ilayer];
-         cl    = (AliTRDcluster *)GetCluster(index);
-         //Double_t h01 = GetTiltFactor(cl);   // I.B's fix
-          //maxChi2=t.GetPredictedChi2(cl,h01); //
-       }
+                               if (clusters[ilayer] > 0) {
+                                       index = clusters[ilayer];
+                                       cl    = (AliTRDcluster *)GetCluster(index);
+                                       //Double_t h01 = GetTiltFactor(cl);   // I.B's fix
+                                                               //maxChi2=t.GetPredictedChi2(cl,h01); //
+                               }
        
         if (cl) {
-
-         //if (cl->GetNPads() < 5) 
-         Double_t dxsample = timeBin.GetdX();
-         t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample)); 
-          Double_t h01      = GetTiltFactor(cl);
-         Int_t    det      = cl->GetDetector();    
-         Int_t    plane    = fGeom->GetPlane(det);
-         if (t.GetX() > 345.0) {
-           t.SetNLast(t.GetNLast() + 1);
-           t.SetChi2Last(t.GetChi2Last() + maxChi2);
-         }
-         Double_t xcluster = cl->GetX();
-         t.PropagateTo(xcluster,xx0,xrho);
-         maxChi2 = t.GetPredictedChi2(cl,h01);
-
-         if (maxChi2<1e+10)
-           if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
-             if (!t.Update(cl,maxChi2,index,h01)) {
-               // ????
-             }
-           } // else SetCluster(cl, GetNumberOfClusters()-1); // A.Bercuci 25.07.07
+                                       //if (cl->GetNPads() < 5)
+                                       Double_t dxsample = timeBin.GetdX();
+                                       t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
+                                                               Double_t h01      = GetTiltFactor(cl);
+                                       Int_t    det      = cl->GetDetector();
+                                       Int_t    plane    = fGeom->GetPlane(det);
+                                       if (t.GetX() > 345.0) {
+                                               t.SetNLast(t.GetNLast() + 1);
+                                               t.SetChi2Last(t.GetChi2Last() + maxChi2);
+                                       }
+                                       Double_t xcluster = cl->GetX();
+                                       t.PropagateTo(xcluster,xx0,xrho);
+                                       maxChi2 = t.GetPredictedChi2(cl,h01);
+
+                                       if (maxChi2<1e+10)
+                                               if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
+                                                       if (!t.Update(cl,maxChi2,index,h01)) {
+                                       // ????
+                                                       }
+                                               } // else SetCluster(cl, GetNumberOfClusters()-1); // A.Bercuci 25.07.07
   
 
           if (calibra->GetMITracking()) {
             calibra->UpdateHistograms(cl,&t);
           }
 
-         // Reset material budget if 2 consecutive gold
-         if (plane > 0) { 
-           if ((t.GetTracklets(plane).GetN() + t.GetTracklets(plane-1).GetN()) > 20) {
-             t.SetBudget(2,0.0);
-           }
-         }
-
-       }
+                                       // Reset material budget if 2 consecutive gold
+                                       if (plane > 0) {
+                                               if ((t.GetTracklets(plane).GetN() + t.GetTracklets(plane-1).GetN()) > 20) {
+                                                       t.SetBudget(2,0.0);
+                                               }
+                                       }
+                       
+                               }
 
       }
 
@@ -1290,7 +1212,8 @@ Int_t AliTRDtracker::LoadClusters(TTree *cTree)
   // differs from that of TRD sectors
   //
 
-  if (ReadClusters(fClusters,cTree)) {
+       
+  if (ReadClusters(fClusters, cTree)) {
     AliError("Problem with reading the clusters !");
     return 1;
   }
@@ -1367,7 +1290,7 @@ void AliTRDtracker::UnloadClusters()
 }
 
 //_____________________________________________________________________________
-void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESDEvent *esd)
+Int_t AliTRDtracker::Clusters2Tracks(AliESDEvent *esd)
 {
   //
   // Creates  seeds using clusters between  position inner plane  and outer plane 
@@ -2531,6 +2454,7 @@ void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESDEvent *e
 
   delete [] pseed;
 
+       return 0;
 }
           
 //_____________________________________________________________________________
@@ -2555,7 +2479,7 @@ Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *clusterTree) const
   // Loop through all entries in the tree
   Int_t nEntries   = (Int_t) clusterTree->GetEntries();
   Int_t nbytes     = 0;
-  AliTRDcluster *c = 0;
+  AliTRDcluster *c = 0x0;
   for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {    
     
     // Import the tree
@@ -2563,17 +2487,15 @@ Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *clusterTree) const
     
     // Get the number of points in the detector
     Int_t nCluster = clusterArray->GetEntriesFast();  
-    
     // Loop through all TRD digits
     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
-      c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster);
-      AliTRDcluster *co = c;
-      array->AddLast(co);
+      if(!(c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster))) continue;
+      array->AddLast(c);
+                       //printf("Add cluster 0x%x.\n", c);
       clusterArray->RemoveAt(iCluster); 
     }
 
   }
-
   delete clusterArray;
 
   return 0;
@@ -2781,68 +2703,6 @@ Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
 
 }
 
-//_____________________________________________________________________________
-AliTRDtracker::AliTRDpropagationLayer
-             ::AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
-                                    , Double_t radLength, Int_t tbIndex, Int_t plane)
-  :fN(0)
-  ,fSec(0)
-  ,fClusters(NULL)
-  ,fIndex(NULL)
-  ,fX(x)
-  ,fdX(dx)
-  ,fRho(rho)
-  ,fX0(radLength)
-  ,fTimeBinIndex(tbIndex)
-  ,fPlane(plane)
-  ,fYmax(0)
-  ,fYmaxSensitive(0)
-  ,fHole(kFALSE)
-  ,fHoleZc(0)
-  ,fHoleZmax(0)
-  ,fHoleYc(0)
-  ,fHoleYmax(0)
-  ,fHoleRho(0)
-  ,fHoleX0(0)
-{ 
-  //
-  // AliTRDpropagationLayer constructor
-  //
-
-  for (Int_t i = 0; i < (Int_t) kZones; i++) {
-    fZc[i]   = 0; 
-    fZmax[i] = 0;
-  }
-
-  if (fTimeBinIndex >= 0) { 
-    fClusters = new AliTRDcluster*[kMaxClusterPerTimeBin];
-    fIndex    = new UInt_t[kMaxClusterPerTimeBin];
-  }
-
-  for (Int_t i = 0; i < 5; i++) {
-    fIsHole[i] = kFALSE;
-  }
-
-}
-
-//_____________________________________________________________________________
-void AliTRDtracker::AliTRDpropagationLayer
-                  ::SetHole(Double_t Zmax, Double_t Ymax, Double_t rho
-                          , Double_t radLength, Double_t Yc, Double_t Zc) 
-{
-  //
-  // Sets hole in the layer 
-  //
-
-  fHole     = kTRUE;
-  fHoleZc   = Zc;
-  fHoleZmax = Zmax;
-  fHoleYc   = Yc;
-  fHoleYmax = Ymax;
-  fHoleRho  = rho;
-  fHoleX0   = radLength;
-
-}
 
 //_____________________________________________________________________________
 AliTRDtracker::AliTRDtrackingSector
@@ -3160,140 +3020,6 @@ Int_t AliTRDtracker::AliTRDtrackingSector
 }             
 
 //_____________________________________________________________________________
-void AliTRDtracker::AliTRDpropagationLayer
-                  ::SetZ(Double_t *center, Double_t *w, Double_t *wsensitive )
-{
-  //
-  // set centers and the width of sectors
-  //
-
-  for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
-    fZc[icham]            = center[icham];  
-    fZmax[icham]          = w[icham];
-    fZmaxSensitive[icham] = wsensitive[icham];
-  }  
-
-}
-
-//_____________________________________________________________________________
-void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
-{
-  //
-  // set centers and the width of sectors
-  //
-
-  fHole = kFALSE;
-
-  for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
-    fIsHole[icham] = holes[icham]; 
-    if (holes[icham]) {
-      fHole = kTRUE;
-    }
-  }  
-
-}
-
-//_____________________________________________________________________________
-void AliTRDtracker::AliTRDpropagationLayer
-                  ::InsertCluster(AliTRDcluster *c, UInt_t index) 
-{
-  //
-  // Insert cluster in cluster array.
-  // Clusters are sorted according to Y coordinate.  
-  //
-
-  if (fTimeBinIndex < 0) { 
-    //AliWarning("Attempt to insert cluster into non-sensitive time bin!\n");
-    return;
-  }
-
-  if (fN == (Int_t) kMaxClusterPerTimeBin) {
-    //AliWarning("Too many clusters !\n"); 
-    return;
-  }
-
-  if (fN == 0) {
-    fIndex[0]       = index; 
-    fClusters[fN++] = c; 
-    return;
-  }
-
-  Int_t i = Find(c->GetY());
-  memmove(fClusters+i+1,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
-  memmove(fIndex   +i+1,fIndex   +i,(fN-i)*sizeof(UInt_t)); 
-  fIndex[i]    = index; 
-  fClusters[i] = c; 
-  fN++;
-
-}  
-
-//_____________________________________________________________________________
-Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const 
-{
-  //
-  // Returns index of the cluster nearest in Y    
-  //
-
-  if (fN <= 0) {
-    return 0;
-  }
-  if (y <= fClusters[0]->GetY()) {
-    return 0;
-  }
-  if (y >  fClusters[fN-1]->GetY()) {
-    return fN;
-  }
-
-  Int_t b = 0;
-  Int_t e = fN - 1;
-  Int_t m = (b + e) / 2;
-
-  for ( ; b < e; m = (b + e) / 2) {
-    if (y > fClusters[m]->GetY()) {
-      b = m + 1;
-    }
-    else {
-      e = m;
-    }
-  }
-
-  return m;
-
-}    
-
-//_____________________________________________________________________________
-Int_t AliTRDtracker::AliTRDpropagationLayer
-                   ::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad
-                                      , Float_t maxroadz) const 
-{
-  //
-  // Returns index of the cluster nearest to the given y,z
-  //
-
-  Int_t   index   = -1;
-  Int_t   maxn    = fN;
-  Float_t mindist = maxroad;                   
-
-  for (Int_t i = Find(y-maxroad); i < maxn; i++) {
-    AliTRDcluster *c = (AliTRDcluster *) (fClusters[i]);
-    Float_t ycl = c->GetY();
-    if (ycl > (y + maxroad)) {
-      break;
-    }
-    if (TMath::Abs(c->GetZ() - z) > maxroadz) {
-      continue;
-    }
-    if (TMath::Abs(ycl - y)       < mindist) {
-      mindist = TMath::Abs(ycl - y);
-      index   = fIndex[i];
-    }
-  }                                            
-
-  return index;
-
-}             
-
-//_____________________________________________________________________________
 Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c) 
 {
   //
@@ -3313,7 +3039,6 @@ Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c)
 
 }
 
-
 //_____________________________________________________________________________
 Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
                                 , AliTRDtrack *track
@@ -4042,9 +3767,24 @@ Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist
 AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params)
 {
   //
-  // Register a seed
+  // Build a TRD track out of tracklet candidates
+  //
+  // Parameters :
+  //   seeds  : array of tracklets
+  //   params : track parameters (see MakeSeeds() function body for a detailed description)
+  //
+  // Output :
+  //   The TRD track.
+  //
+  // Detailed description
+  //
+  // 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];
@@ -4059,13 +3799,15 @@ AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params)
   AliTRDcluster *cl = 0;
 
   for (Int_t ilayer = 0; ilayer < 6; ilayer++) {
-    if (seeds[ilayer].IsOK()) {
-      for (Int_t itime = 22; itime > 0; itime--) {
-       if (seeds[ilayer].GetIndexes(itime) > 0) {
-         index = seeds[ilayer].GetIndexes(itime);
-         cl    = seeds[ilayer].GetClusters(itime);
-         break;
-       }
+               
+               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);
+                                       //printf("l[%d] index[%d] tb[%d] cptr[%p]\n", ilayer, index, itime, cl);
+                                       break;
+                               }
       }
     }
     if (index > 0) {
index 0067e56..288f59c 100644 (file)
 //                                                                        //
 ////////////////////////////////////////////////////////////////////////////
 
+#ifndef ROOT_TObjArray
 #include "TObjArray.h" 
+#endif
 
+#ifndef ALITRACKER_H
 #include "AliTracker.h" 
+#endif
+
+#ifndef ALITRDPROPAGATIONLAYER_H
+#include "AliTRDpropagationLayer.h"
+#endif
 
 class TFile;
 class TTree;
@@ -30,6 +38,7 @@ class AliTRDtracklet;
 class AliTRDcluster;
 class AliTRDseed;
 class AliESDEvent;
+class AliTRDpropagationLayer;
 
 ///////////////////////////////////////////////////////////////////////////////
 //                                                                           //
@@ -47,8 +56,6 @@ class AliTRDtracker : public AliTracker {
 
   enum { kMaxLayersPerSector   = 1000
        , kMaxTimeBinIndex      = 216
-       , kMaxClusterPerTimeBin = 2300
-       , kZones                = 5
        , kTrackingSectors      = 18   };
   
   AliTRDtracker();
@@ -93,91 +100,17 @@ class AliTRDtracker : public AliTracker {
   Int_t            ReadClusters(TObjArray *array, TTree *in) const;
   AliTRDcluster   *GetCluster(AliTRDtrack *track, Int_t plane, Int_t timebin, UInt_t &index);
   Int_t            FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack *track
-                              , Int_t *clusters, AliTRDtracklet &tracklet);
+                              , Int_t *clusters, AliTRDtracklet &tracklet);  
   
  protected:
-  
-  //__________________________________________________________________________________________________
-  class AliTRDpropagationLayer {
-    
-   public: 
-    
-    AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
-                        , Double_t x0, Int_t tbIndex, Int_t plane); 
-    AliTRDpropagationLayer(const AliTRDpropagationLayer &/*p*/);
-    ~AliTRDpropagationLayer() {  if (fTimeBinIndex >= 0) { 
-                                   delete[] fClusters; 
-                                  delete[] fIndex; 
-                                 } 
-                               }
-    AliTRDpropagationLayer &operator=(const AliTRDpropagationLayer &/*p*/) 
-                                                              { return *this; }
-
-    operator Int_t() const                                    { return fN;    }
-    AliTRDcluster  *operator[](Int_t i)                       { return fClusters[i];          }
-    
-    void     SetZmax(Int_t cham, Double_t center, Double_t w) { fZc[cham]      = center; 
-                                                                fZmax[cham]    = w;           }
-    void     SetYmax(Double_t w, Double_t wsensitive)         { fYmax          = w;
-                                                                fYmaxSensitive = wsensitive;  }
 
-    void     SetZ(Double_t* center, Double_t *w, Double_t *wsensitive);
-    void     SetHoles(Bool_t* holes);
-    void     SetHole(Double_t Zmax, Double_t Ymax
-                  , Double_t rho = 1.29e-3, Double_t x0 = 36.66
-                  , Double_t Yc = 0.0, Double_t Zc = 0.0);
-    
-    Double_t GetYmax() const           { return fYmax;                  }
-    Double_t GetZmax(Int_t c) const    { return fZmax[c];               }
-    Double_t GetZc(Int_t c) const      { return fZc[c];                 }
-    UInt_t   GetIndex(Int_t i) const   { return fIndex[i];              }  
-    Double_t GetX() const              { return fX;                     }
-    Double_t GetdX() const             { return fdX;                    }
-    Int_t    GetTimeBinIndex() const   { return fTimeBinIndex;          }     
-    Int_t    GetPlane() const          { return fPlane;                 }
-    Bool_t   IsHole(Int_t zone) const  { return fIsHole[zone];          }              
-    Bool_t   IsSensitive() const       { return (fTimeBinIndex >= 0) ? kTRUE : kFALSE;} 
-    
-    void     Clear() { 
-      for (Int_t i = 0; i < fN; i++) fClusters[i] = NULL;
-      fN = 0;
-    }
-    
-    void     InsertCluster(AliTRDcluster *c, UInt_t index);
-    Int_t    Find(Float_t y) const; 
-    Int_t    FindNearestCluster(Float_t y, Float_t z, Float_t maxroad, Float_t maxroadz) const;
-    void     SetX(Double_t x)          { fX = x; }
-    
-   private:     
-    
-    Int_t                     fN;                            // This is fN
-    Int_t                     fSec;                          // Sector mumber
-    AliTRDcluster           **fClusters;                     // Array of pointers to clusters
-    UInt_t                   *fIndex;                        // Array of cluster indexes
-    Double_t                  fX;                            // X coordinate of the middle plane
-    Double_t                  fdX;                           // Radial thickness of the time bin
-    Double_t                  fRho;                          // Default density of the material 
-    Double_t                  fX0;                           // Default radiation length 
-    Int_t                     fTimeBinIndex;                 // Plane * F(local_tb)  
-    Int_t                     fPlane;                        // Plane number
-    Double_t                  fZc[kZones];                   // Z position of the center for 5 active areas
-    Double_t                  fZmax[kZones];                 // Half of active area length in Z
-    Double_t                  fZmaxSensitive[kZones];        // Sensitive area for detection Z     
-    Bool_t                    fIsHole[kZones];               // Is hole in given sector       
-    Double_t                  fYmax;                         // Half of active area length in Y
-    Double_t                  fYmaxSensitive;                // Half of active area length in Y
-    
-    Bool_t                    fHole;                         // kTRUE if there is a hole in the layer
-    Double_t                  fHoleZc;                       // Z of the center of the hole 
-    Double_t                  fHoleZmax;                     // Half of the hole length in Z
-    Double_t                  fHoleYc;                       // Y of the center of the hole 
-    Double_t                  fHoleYmax;                     // Half of the hole length in Y 
-    Double_t                  fHoleRho;                      // Density of the gas in the hole 
-    Double_t                  fHoleX0;                       // Radiation length of the gas in the hole 
-    
-  };
+  Bool_t           AdjustSector(AliTRDtrack *track); 
+  AliTRDtrack     *RegisterSeed(AliTRDseed *seeds, Double_t *params);
+  Int_t            FollowBackProlongation(AliTRDtrack &t);
+  //void             MakeSeedsMI(Int_t inner, Int_t outer, AliESDEvent *esd = 0);
   
+ protected:
+
   //__________________________________________________________________________________________________
   class AliTRDtrackingSector {
     
@@ -197,6 +130,7 @@ class AliTRDtracker : public AliTracker {
     Int_t    GetLayerNumber(Int_t tb) const        { return fTimeBinIndex[tb];   }
     Double_t GetX(Int_t pl) const                  { return fLayers[pl]->GetX(); }
     AliTRDpropagationLayer* GetLayer(Int_t i)      { return fLayers[i];          }
+    Int_t    GetSector() const {return fGeomSector;}   
     
     void     MapTimeBinLayers();
     Int_t    Find(Double_t x) const; 
@@ -213,6 +147,8 @@ class AliTRDtracker : public AliTracker {
     
   };
   
+ protected:
+
   AliTRDgeometry          *fGeom;                          // Pointer to TRD geometry
   AliTRDtrackingSector    *fTrSec[kTrackingSectors];       // Array of tracking sectors;    
   Int_t                    fNclusters;                     // Number of clusters in TRD 
@@ -246,23 +182,21 @@ class AliTRDtracker : public AliTracker {
   TH2D                    *fHMinD;                         // QA histogram
   TH1D                    *fHDeltaX;                       // QA histogram
   TH1D                    *fHXCl;                          // QA histogram
-  
-  Bool_t   AdjustSector(AliTRDtrack *track); 
-  
+
  private:
-  
-  AliTRDtrack *RegisterSeed(AliTRDseed *seeds, Double_t *params);
-  void     MakeSeedsMI(Int_t inner, Int_t outer, AliESDEvent *esd = 0);
-  
-  Int_t    FollowBackProlongation(AliTRDtrack &t);
-  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;
-  
+
+  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
   
-  ClassDef(AliTRDtracker,3)                                // TRD tracker
+  ClassDef(AliTRDtracker, 4)                               // TRD tracker
     
 };
+
+
 #endif 
diff --git a/TRD/AliTRDtrackerFitter.cxx b/TRD/AliTRDtrackerFitter.cxx
new file mode 100644 (file)
index 0000000..8f0943e
--- /dev/null
@@ -0,0 +1,402 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Track fitter                                                             //
+//                                                                           //
+//  Authors:                                                                 //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                        //
+//    Markus Fasel <M.Fasel@gsi.de>                                          //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include <TROOT.h>
+#include <TLinearFitter.h>
+#include <TMath.h>
+#include <TTreeStream.h>
+
+#include "AliLog.h"
+
+#include "AliTRDseed.h"
+#include "AliTRDcalibDB.h"
+#include "AliTRDcluster.h"
+#include "AliTRDReconstructor.h"
+#include "AliTRDseedV1.h"
+#include "AliTRDtrackerFitter.h"
+
+#define DEBUG
+
+ClassImp(AliTRDtrackerFitter)
+
+//_________________________________________________________________________
+AliTRDtrackerFitter::AliTRDtrackerFitter()
+  :TObject()
+  ,fFitterTC(0x0)
+  ,fFitterT2(0x0)
+  ,fRieman1(0x0)
+  ,fRieman2(0x0)
+  ,fChi2TR(0)
+  ,fChi2TC(0)
+  ,fCR(0)
+  ,fCC(0)
+  ,fDca(0)
+  ,fDzmf(0)
+  ,fZmf(0)
+  ,fNlayers(0)
+  ,fDebugStream(0x0)
+{
+  //
+  // Constructor
+  //
+
+       fRieman1   = new AliRieman(1000);
+       fRieman2  = new AliRieman(1000);
+       fFitterTC = new TLinearFitter(2,"hyp2");
+       fFitterT2 = new TLinearFitter(4,"hyp4");
+       fFitterTC->StoreData(kTRUE);
+       fFitterT2->StoreData(kTRUE);
+}
+
+//_________________________________________________________________________
+AliTRDtrackerFitter::AliTRDtrackerFitter(const AliTRDtrackerFitter &f)
+  :TObject(f)
+  ,fFitterTC(0x0)
+  ,fFitterT2(0x0)
+  ,fRieman1(0x0)
+  ,fRieman2(0x0)
+  ,fChi2TR(f.fChi2TR)
+  ,fChi2TC(f.fChi2TC)
+  ,fCR(f.fCR)
+  ,fCC(f.fCC)
+  ,fDca(f.fDca)
+  ,fDzmf(f.fDzmf)
+  ,fZmf(f.fZmf)
+  ,fNlayers(f.fNlayers)
+  ,fDebugStream(0x0)
+{
+  //
+  // Copy Constructor (performs a deep copy) 
+  //
+
+       fRieman1  = new AliRieman(*f.fRieman1);
+       fRieman2  = new AliRieman(*f.fRieman2);
+       fFitterTC = new TLinearFitter(*f.fFitterTC);
+       fFitterT2 = new TLinearFitter(*f.fFitterT2);
+       fFitterTC->StoreData(kTRUE);
+       fFitterT2->StoreData(kTRUE);
+}
+
+//_________________________________________________________________________
+AliTRDtrackerFitter::~AliTRDtrackerFitter()
+{
+  //
+  // Destructor
+  //
+
+       delete fRieman1;
+       delete fRieman2;
+       delete fFitterTC;
+       delete fFitterT2;
+}
+
+//_________________________________________________________________________
+AliTRDtrackerFitter &AliTRDtrackerFitter::operator=(const AliTRDtrackerFitter& fitter)
+{
+  //
+  // Assignment operator using the copy Method of this class
+  //
+       if(this != &fitter){
+               fitter.Copy(*this);
+       }
+       return *this;
+}
+
+//_________________________________________________________________________
+void AliTRDtrackerFitter::Copy(TObject &f) const 
+{
+  //
+  // Copy method. 
+  // Performs a deep copy. Mainly used in the Assignment operator
+  //
+
+  AliTRDtrackerFitter &tf = (AliTRDtrackerFitter &) f;
+
+       if(tf.fFitterTC)
+               delete tf.fFitterTC;
+       tf.fFitterTC = new TLinearFitter(*fFitterTC); 
+       if(tf.fFitterT2)
+               delete tf.fFitterT2;
+       tf.fFitterT2 = new TLinearFitter(*fFitterT2);
+       if(tf.fRieman1)
+               delete tf.fRieman1;
+       tf.fRieman1 = new AliRieman(*fRieman1);
+       if(tf.fRieman2)
+               delete fRieman2;
+       tf.fRieman2 = new AliRieman(*fRieman2);
+       tf.fChi2TR = fChi2TR;
+       tf.fChi2TC = fChi2TC;
+       tf.fCR = fCR;
+       tf.fCC = fCC;
+       tf.fDca =fDca;
+       tf.fDzmf = fDzmf;
+       tf.fZmf =fZmf;
+       tf.fNlayers = fNlayers;
+       tf.fDebugStream = 0x0;
+       fFitterTC->StoreData(kTRUE);
+       fFitterT2->StoreData(kTRUE);
+}
+
+//_________________________________________________________________________
+void AliTRDtrackerFitter::FitRieman(AliTRDcluster **cl, Int_t nLayers)
+{
+  //
+  // Performs a Rieman Fit including 4 clusters (one in each layer)
+  //
+       fRieman1->Reset();
+       for(Int_t i = 0; i < nLayers; i++)
+               fRieman1->AddPoint(cl[i]->GetX(), cl[i]->GetY(), cl[i]->GetZ(), 1, 10);
+       fRieman1->Update();
+}
+
+//_________________________________________________________________________
+void AliTRDtrackerFitter::FitRieman(AliTRDseedV1 *cseed, Int_t *planes)
+{
+  //
+  // Performs a Rieman fit using four respectively six seeds
+  // 2 times used: 1st with 4 parameters and Later in the full track
+  // fitting Part with 6 parameters
+  //
+
+       Int_t lplanes[] = {0, 1, 2, 3, 4, 5}, *tplanes;
+       Int_t nplanes;
+       if(planes){
+               nplanes = 4;
+               tplanes = planes;
+       } else {
+               nplanes = 6;
+               tplanes = &lplanes[0];
+       }
+       fRieman1->Reset();
+       for(Int_t iLayer = 0; iLayer < nplanes; iLayer++){
+               if(!cseed[tplanes[iLayer]].IsOK()) continue;
+               fRieman1->AddPoint(cseed[tplanes[iLayer]].GetX0(), cseed[tplanes[iLayer]].GetYfitR(0), cseed[tplanes[iLayer]].GetZProb(), 1, 10);
+       }
+       fRieman1->Update();
+}
+
+//_________________________________________________________________________
+Double_t AliTRDtrackerFitter::FitHyperplane(AliTRDseedV1 *cseed
+                                          , Double_t chi2ZF
+                                          , Double_t zFitter)
+{
+  //
+  // performs a hyperplane fit
+  // 
+  // Two linear fitters:  one in which kz is fixed and one in which kz is a free parameter
+  // checking if values are acceptable and improves the fitresults
+  // Calculates curvature parameters and chisquares
+  // Return: likelihood value as a measurement for the seedquality
+  //
+
+  //const Int_t kMaxTimebins = 100; // Buffer
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       Int_t nTimeBins = cal->GetNumberOfTimeBins();
+       Int_t npointsT = 0;
+       fFitterTC->ClearPoints();
+       fFitterT2->ClearPoints();
+       
+       Double_t fHl[6];
+       Double_t xref2 = (cseed[2].GetX0() + cseed[3].GetX0())/2;
+       fRieman2->Reset();
+       for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+               if (!cseed[iLayer].IsOK()) continue;
+               fHl[iLayer] = cseed[iLayer].GetTilt();
+
+               for (Int_t itime = 0; itime < nTimeBins; itime++) {
+                       if (!cseed[iLayer].IsUsable(itime)) continue;
+
+                       // X relative to the middle chamber
+                       Double_t x  = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0() - xref2;
+                       Double_t y  = cseed[iLayer].GetY(itime);
+                       Double_t z  = cseed[iLayer].GetZ(itime);
+                       // ExB correction to the correction
+                       // Tilted fRieman1
+                       Double_t uvt[6];
+                       // Global x
+                       Double_t x2 = cseed[iLayer].GetX(itime) + cseed[iLayer].GetX0();
+                       Double_t t  = 1.0 / (x2*x2 + y*y);
+                       uvt[1] = t;                 // t
+                       uvt[0] = 2.0 * x2 * uvt[1]; // u
+                       uvt[2] = 2.0 * fHl[iLayer] * uvt[1];
+                       uvt[3] = 2.0 * fHl[iLayer] * x * uvt[1];
+                       uvt[4] = 2.0 * (y + fHl[iLayer]*z) * uvt[1];
+                       Double_t error = 2.0 * 0.2 * uvt[1];
+                       fFitterT2->AddPoint(uvt,uvt[4],error);
+
+                       //Constrained fRieman1
+                       z = cseed[iLayer].GetZ(itime);
+                       uvt[0] = 2.0 * x2 * t; // u
+                       uvt[1] = 2.0 * fHl[iLayer] * x2 * uvt[1];
+                       uvt[2] = 2.0 * (y + fHl[iLayer] * (z - zFitter)) * t;   // parameter that is most propably wrong
+                       fFitterTC->AddPoint(uvt,uvt[2],error);
+                       fRieman2->AddPoint(x2,y,z,1,10);
+                       npointsT++;
+               }
+       } // Loop: iLayer
+
+       fRieman2->Update();
+       fFitterTC->Eval();
+       fFitterT2->Eval();
+       Double_t rpolz0 = fFitterT2->GetParameter(3);
+       Double_t rpolz1 = fFitterT2->GetParameter(4);
+
+       // Linear fitter  - not possible to make boundaries
+       // Do not accept non possible z and dzdx combinations
+       Bool_t acceptablez = kTRUE;
+       for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+               if(!cseed[iLayer].IsOK()) continue;
+               Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].GetX0() - xref2);
+               if (TMath::Abs(cseed[iLayer].GetZProb() - zT2) > cseed[iLayer].GetPadLength() * 0.5 + 1.0) acceptablez = kFALSE;
+       }
+       if (!acceptablez) {
+               fFitterT2->FixParameter(3,fRieman1->GetZat(xref2));
+               fFitterT2->FixParameter(4,fRieman1->GetDZat(xref2));
+               fFitterT2->Eval();
+               fFitterT2->ReleaseParameter(3);
+               fFitterT2->ReleaseParameter(4);
+               rpolz0 = fFitterT2->GetParameter(3);
+               rpolz1 = fFitterT2->GetParameter(4);
+       }
+
+       fChi2TR = fFitterT2->GetChisquare() / Float_t(npointsT);
+       fChi2TC = fFitterTC->GetChisquare() / Float_t(npointsT);
+       Double_t polz1c = fFitterTC->GetParameter(2);
+       Double_t polz0c = polz1c * xref2;
+       Double_t aC     =  fFitterTC->GetParameter(0);
+       Double_t bC     =  fFitterTC->GetParameter(1);
+       fCC     =  aC / TMath::Sqrt(bC * bC + 1.0); // Curvature
+       Double_t aR     =  fFitterT2->GetParameter(0);
+       Double_t bR     =  fFitterT2->GetParameter(1);
+       Double_t dR     =  fFitterT2->GetParameter(2);
+       fCR     =  1.0 + bR*bR - dR*aR;
+       fDca    =  0.0;
+       if (fCR > 0.0) {
+               fDca = -dR / (TMath::Sqrt(1.0 + bR*bR - dR*aR) + TMath::Sqrt(1.0 + bR*bR));
+               fCR  =  aR / TMath::Sqrt(fCR);
+       }
+
+       
+       Double_t chi2ZT2 = 0.0;
+       Double_t chi2ZTC = 0.0;
+       for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+               if(!cseed[iLayer].IsOK()) continue;
+               Double_t zT2 = rpolz0 + rpolz1 * (cseed[iLayer].GetX0() - xref2);
+               Double_t zTC = polz0c + polz1c * (cseed[iLayer].GetX0() - xref2);
+               chi2ZT2 += TMath::Abs(cseed[iLayer].GetMeanz() - zT2);
+               chi2ZTC += TMath::Abs(cseed[iLayer].GetMeanz() - zTC);
+       }
+       chi2ZT2 /= TMath::Max((fNlayers - 3.0),1.0);
+       chi2ZTC /= TMath::Max((fNlayers - 3.0),1.0);
+        
+       
+       AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
+       Float_t sumdaf = 0.0;
+       for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+               if(!cseed[iLayer].IsOK()) continue;
+               sumdaf += TMath::Abs((cseed[iLayer].GetYfit(1) - cseed[iLayer].GetYref(1))/ cseed[iLayer].GetSigmaY2());
+       }
+       sumdaf /= Float_t (fNlayers - 2.0);
+
+       // Cook Likelihoods
+       Double_t likezf     = TMath::Exp(-chi2ZF * 0.14);
+       Double_t likechi2C  = TMath::Exp(-fChi2TC * 0.677);
+       Double_t likechi2TR = TMath::Exp(-fChi2TR * 0.78);
+       Double_t likeaf     = TMath::Exp(-sumdaf * 3.23);
+       Double_t likef = likezf * likechi2TR * likeaf;
+
+#ifdef DEBUG
+       if(fDebugStream && AliTRDReconstructor::StreamLevel() >= 2){
+               TTreeSRedirector &treeStreamer = *fDebugStream;
+               treeStreamer << "FitHyperplane"
+                       << "seed0.="     << &cseed[0]
+                       << "seed1.="     << &cseed[1]
+                       << "seed2.="     << &cseed[2]
+                       << "seed3.="     << &cseed[3]
+                       << "seed4.="     << &cseed[4]
+                       << "seed5.="     << &cseed[5]
+                       << "chi2TR="     << fChi2TR
+                       << "chi2TC="     << fChi2TC
+                       << "chi2ZT2="    << chi2ZT2
+                       << "chi2ZTC="    << chi2ZTC
+                       << "CC="         << fCC
+                       << "CR="         << fCR
+                       << "DR="         << dR
+                       << "DCA="        << fDca
+                       << "Polz0="      << polz0c
+                       << "Polz1="      << polz1c
+                       << "RPolz0="     << rpolz0
+                       << "RPolz1="     << rpolz1
+                       << "Likechi2C="  << likechi2C
+                       << "Likechi2TR=" << likechi2TR
+                       << "Likezf="     << likezf
+                       << "Likeaf="     << likeaf
+                       << "LikeF="      << likef
+                       << "Rieman1.="   << fRieman1
+                       << "Rieman2.="   << fRieman2
+                       << "\n";
+       }
+#endif
+
+       return likef;
+}
+
+//_________________________________________________________________________
+void AliTRDtrackerFitter::GetHyperplaneFitChi2(Double_t *chisquares) const
+{
+  //
+  // Getter method returns the chisquares of the hyperplane fit
+  //
+
+       chisquares[0] = fChi2TR;
+       chisquares[1] = fChi2TC;
+}
+
+//_________________________________________________________________________
+void AliTRDtrackerFitter::GetHyperplaneFitResults(Double_t *params) const
+{
+  //
+  // Getter Method returning the Curvature parameters of the hyperplane fit
+  //
+
+       params[0] = fCC;
+       params[1] = fCR;
+       params[2] = fDca;
+}
+
+//_________________________________________________________________________
+void AliTRDtrackerFitter::Reset()
+{
+  //
+  // Resets the object
+  //
+
+       fRieman1->Reset();
+       fRieman2->Reset();
+       //fFitterTC->Clear();
+       //fFitterT2->Clear();
+}
diff --git a/TRD/AliTRDtrackerFitter.h b/TRD/AliTRDtrackerFitter.h
new file mode 100644 (file)
index 0000000..a423655
--- /dev/null
@@ -0,0 +1,80 @@
+#ifndef ALITRDTRACKERFITTER_H
+#define ALITRDTRACKERFITTER_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Track fitter                                                             //
+//                                                                           //
+//  Authors:                                                                 //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                        //
+//    Markus Fasel <M.Fasel@gsi.de>                                          //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef AliRIEMAN_H
+#include <AliRieman.h>
+#endif
+
+#ifndef ROOT_Rtypes
+#include "Rtypes.h"
+#endif
+
+// implementations in the future
+class TObject;
+class TLinearFitter;
+class TTreeSRedirector;
+
+class AliTRDcluster;
+class AliTRDseedV1;
+
+class AliTRDtrackerFitter : public TObject {
+ public:
+
+       AliTRDtrackerFitter();
+       AliTRDtrackerFitter(const AliTRDtrackerFitter &);
+       AliTRDtrackerFitter& operator=(const AliTRDtrackerFitter &);
+       virtual ~AliTRDtrackerFitter();
+
+       void            Copy(TObject &f) const;
+       AliRieman      *GetRiemanFitter() const { return fRieman1; }
+       void            GetHyperplaneFitChi2(Double_t *chisquares) const;
+       void            GetHyperplaneFitResults(Double_t *params) const;
+       Double_t        GetRiemanCurvature() const { return fRieman1->GetC(); }
+
+       void            FitRieman(AliTRDcluster **cl, Int_t nLayers);
+       void            FitRieman(AliTRDseedV1 *cseed, Int_t *planes=0x0);
+       Double_t        FitHyperplane(AliTRDseedV1 *cseed, Double_t chi2ZF, Double_t zval);
+       void            Reset();
+       void            SetDebugStream(TTreeSRedirector *debug){fDebugStream = debug;}
+       void            SetLayers(const Int_t nlayers) {fNlayers = nlayers;}
+
+ private:
+
+       // Linear fitters in planes
+       TLinearFitter    *fFitterTC;       // Fitting with tilting pads - kz fixed - kz= Z/x, + vertex const
+       TLinearFitter    *fFitterT2;       // Fitting with tilting pads - kz not fixed
+       AliRieman        *fRieman1;        // Rieman fitter
+       AliRieman        *fRieman2;        // Rieman fitter
+       // Linear Fitter chisquares
+       Double_t          fChi2TR;         // Chisquared
+       Double_t          fChi2TC;         // Chisquared
+       // Hyperplane Fit results
+       Double_t          fCR;             // Track radius
+       Double_t          fCC;             // Track curvature
+       Double_t          fDca;            // DCA       
+       // Helping Parameters
+       Double_t          fDzmf;           // Something
+       Double_t          fZmf;            // Something else
+       Int_t             fNlayers;        // Some layers
+
+       // Debug
+       TTreeSRedirector *fDebugStream;    //! debugger
+
+       ClassDef(AliTRDtrackerFitter,1)    // TRD track fitter
+};
+#endif // ALITRDTRACKERFITTER_H
diff --git a/TRD/AliTRDtrackerV1.cxx b/TRD/AliTRDtrackerV1.cxx
new file mode 100644 (file)
index 0000000..82bb77a
--- /dev/null
@@ -0,0 +1,1865 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Track finder                                                             //
+//                                                                           //
+//  Authors:                                                                 //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                        //
+//    Markus Fasel <M.Fasel@gsi.de>                                          //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include <Riostream.h>
+#include <stdio.h>
+#include <string.h>
+
+#include <TBranch.h>
+#include <TFile.h>
+#include <TGraph.h>
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TLinearFitter.h>
+#include <TObjArray.h> 
+#include <TROOT.h>
+#include <TTree.h>  
+#include <TClonesArray.h>
+#include <TRandom.h>
+#include <TTreeStream.h>
+
+#include "AliLog.h"
+#include "AliESDEvent.h"
+#include "AliAlignObj.h"
+#include "AliRieman.h"
+#include "AliTrackPointArray.h"
+
+#include "AliTRDtracker.h"
+#include "AliTRDtrackerV1.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDpadPlane.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDcluster.h" 
+#include "AliTRDtrack.h"
+#include "AliTRDseed.h"
+#include "AliTRDcalibDB.h"
+#include "AliTRDCommonParam.h"
+#include "AliTRDReconstructor.h"
+#include "AliTRDCalibraFillHisto.h"
+#include "AliTRDtrackerFitter.h"
+#include "AliTRDstackLayer.h"
+#include "AliTRDrecoParam.h"
+#include "AliTRDseedV1.h"
+
+#define DEBUG
+
+ClassImp(AliTRDtrackerV1)
+Double_t AliTRDtrackerV1::fTopologicQA[kNConfigs] = {
+               0.1112, 0.1112, 0.1112, 0.0786, 0.0786,
+               0.0786, 0.0786, 0.0579, 0.0579, 0.0474,
+               0.0474, 0.0408, 0.0335, 0.0335, 0.0335
+};
+
+//____________________________________________________________________
+AliTRDtrackerV1::AliTRDtrackerV1(AliTRDrecoParam *p) 
+  :AliTRDtracker()
+  ,fSieveSeeding(0)
+  ,fRecoParam(p)
+  ,fFitter(0x0)
+  ,fDebugStreamer(0x0)
+{
+  //
+  // Default constructor. Nothing is initialized.
+  //
+
+}
+
+//____________________________________________________________________
+AliTRDtrackerV1::AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p) 
+  :AliTRDtracker(in)
+  ,fSieveSeeding(0)
+  ,fRecoParam(p)
+  ,fFitter(0x0)
+  ,fDebugStreamer(0x0)
+{
+  //
+  // Standard constructor.
+  // Setting of the geometry file, debug output (if enabled)
+  // and initilize fitter helper.
+  //
+
+       fFitter = new AliTRDtrackerFitter();
+
+#ifdef DEBUG
+       fDebugStreamer    = new TTreeSRedirector("TRDdebug.root");
+       fFitter->SetDebugStream(fDebugStreamer);
+#endif
+
+}
+  
+//____________________________________________________________________
+AliTRDtrackerV1::~AliTRDtrackerV1()
+{ 
+  //
+  // Destructor
+  //
+
+       if(fDebugStreamer) delete fDebugStreamer;
+       if(fFitter) delete fFitter;
+       if(fRecoParam) delete fRecoParam;
+
+}
+
+//____________________________________________________________________
+Int_t AliTRDtrackerV1::Clusters2Tracks(AliESDEvent *esd)
+{
+  //
+  // Steering stand alone tracking for full TRD detector
+  //
+  // Parameters :
+  //   esd     : The ESD event. On output it contains 
+  //             the ESD tracks found in TRD.
+  //
+  // Output :
+  //   Number of tracks found in the TRD detector.
+  // 
+  // Detailed description
+  // 1. Launch individual SM trackers. 
+  //    See AliTRDtrackerV1::Clusters2TracksSM() for details.
+  //
+
+       if(!fRecoParam){
+               AliError("Reconstruction configuration not initialized. Call first AliTRDtrackerV1::SetRecoParam().");
+               return 0;
+       }
+       
+       //AliInfo("Start Track Finder ...");
+       Int_t ntracks = 0;
+       for(int ism=0; ism<AliTRDtracker::kTrackingSectors; ism++){
+                       //AliInfo(Form("Processing supermodule %i ...", ism));
+                       ntracks += Clusters2TracksSM(fTrSec[ism], esd);
+       }
+       AliInfo(Form("Found %d TRD tracks.", ntracks));
+       return ntracks;
+}
+
+//____________________________________________________________________
+Int_t AliTRDtrackerV1::Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *sector
+                                       , AliESDEvent *esd)
+{
+  //
+  // Steer tracking for one SM.
+  //
+  // Parameters :
+  //   sector  : Array of (SM) propagation layers containing clusters
+  //   esd     : The current ESD event. On output it contains the also
+  //             the ESD (TRD) tracks found in this SM. 
+  //
+  // Output :
+  //   Number of tracks found in this TRD supermodule.
+  // 
+  // Detailed description
+  //
+  // 1. Unpack AliTRDpropagationLayers objects for each stack.
+  // 2. Launch stack tracking. 
+  //    See AliTRDtrackerV1::Clusters2TracksStack() for details.
+  // 3. Pack results in the ESD event.
+  //
+       
+       AliTRDpadPlane *pp = 0x0;
+       
+       // allocate space for esd tracks in this SM
+       TClonesArray esdTrackList("AliESDtrack", 2*kMaxTracksStack);
+       esdTrackList.SetOwner();
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       Int_t nTimeBins = cal->GetNumberOfTimeBins();
+       const Int_t kFindable = Int_t(fRecoParam->GetFindableClusters()*6.*nTimeBins);
+       
+       Int_t ntracks = 0;
+       Int_t nClStack = 0;
+       for(int istack = 0; istack<AliTRDpropagationLayer::kZones; istack++){
+               AliTRDstackLayer stackLayer[kNPlanes*kNTimeBins];
+               
+               nClStack = 0;
+               //AliInfo(Form("Processing stack %i ...",istack));
+               //AliInfo("Building stack propagation layers ...");
+               for(int ilayer=0; ilayer<kNPlanes*nTimeBins; ilayer++){
+                       pp = fGeom->GetPadPlane((Int_t)(ilayer/nTimeBins), istack);
+                       Double_t stacklength = (pp->GetNrows() - 2) * pp->GetLengthIPad() 
+                                             + 2 * pp->GetLengthOPad() + 2 * pp->GetLengthRim();
+                       //Debug
+                       Double_t z0  = fGeom->GetRow0((Int_t)(ilayer/nTimeBins),istack,0);
+                       const AliTRDpropagationLayer smLayer(*(sector->GetLayer(ilayer)));
+                       stackLayer[ilayer] = smLayer;
+#ifdef DEBUG
+                       stackLayer[ilayer].SetDebugStream(fDebugStreamer);
+#endif                 
+                       stackLayer[ilayer].SetRange(z0 - stacklength, stacklength);
+                       stackLayer[ilayer].SetSector(sector->GetSector());
+                       stackLayer[ilayer].SetStackNr(istack);
+                       stackLayer[ilayer].SetNRows(pp->GetNrows());
+                       stackLayer[ilayer].SetRecoParam(fRecoParam);
+                       stackLayer[ilayer].BuildIndices();
+                       nClStack += stackLayer[ilayer].GetNClusters();
+               }
+               //AliInfo(Form("Finish building stack propagation layers. nClusters %d.", nClStack));
+               if(nClStack < kFindable) continue;
+               ntracks += Clusters2TracksStack(&stackLayer[0], &esdTrackList);
+       }
+       //AliInfo(Form("Found %d tracks in SM", ntracks));
+       
+       for(int itrack=0; itrack<ntracks; itrack++) 
+          esd->AddTrack((AliESDtrack*)esdTrackList[itrack]);
+
+       return ntracks;
+}
+
+//____________________________________________________________________
+Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
+                                          , TClonesArray *esdTrackList)
+{
+  //
+  // Make tracks in one TRD stack.
+  //
+  // Parameters :
+  //   layer  : Array of stack propagation layers containing clusters
+  //   esdTrackList  : Array of ESD tracks found by the stand alone tracker. 
+  //                   On exit the tracks found in this stack are appended.
+  //
+  // Output :
+  //   Number of tracks found in this stack.
+  // 
+  // Detailed description
+  //
+  // 1. Find the 3 most useful seeding chambers. See BuildSeedingConfigs() for details.
+  // 2. Steer AliTRDtrackerV1::MakeSeeds() for 3 seeding layer configurations. 
+  //    See AliTRDtrackerV1::MakeSeeds() for more details.
+  // 3. Arrange track candidates in decreasing order of their quality
+  // 4. Classify tracks in 5 categories according to:
+  //    a) number of layers crossed
+  //    b) track quality 
+  // 5. Sign clusters by tracks in decreasing order of track quality
+  // 6. Build AliTRDtrack out of seeding tracklets
+  // 7. Cook MC label
+  // 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
+
+       //Double_t alpha = AliTRDgeometry::GetAlpha();
+       //Double_t shift = .5 * alpha;
+       Int_t configs[kNConfigs];
+       
+       // Build initial seeding configurations
+       Double_t quality = BuildSeedingConfigs(layer, configs);
+#ifdef DEBUG
+               if(AliTRDReconstructor::StreamLevel() > 1) 
+                  AliInfo(Form("Plane config %d %d %d Quality %f"
+                              , configs[0], configs[1], configs[2], quality));
+#endif
+
+       // Initialize contors
+       Int_t ntracks,      // number of TRD track candidates
+             ntracks1,     // number of registered TRD tracks/iter
+             ntracks2 = 0; // number of all registered TRD tracks in stack
+       fSieveSeeding = 0;
+       do{
+               // Loop over seeding configurations
+               ntracks = 0; ntracks1 = 0;
+               for (Int_t iconf = 0; iconf<3; iconf++) {
+                       pars[0] = configs[iconf];
+                       pars[1] = layer->GetStackNr();
+                       pars[2] = ntracks;
+                       ntracks = MakeSeeds(layer, &sseed[6*ntracks], pars);
+                       if(ntracks == kMaxTracksStack) break;
+               }
+#ifdef DEBUG
+               if(AliTRDReconstructor::StreamLevel() > 1) 
+                  AliInfo(Form("Candidate TRD tracks %d in stack %d.", ntracks, pars[1]));
+#endif         
+               if(!ntracks) break;
+               
+               // Sort the seeds according to their quality
+               Int_t sort[kMaxTracksStack];
+               TMath::Sort(ntracks, fTrackQuality, sort, kTRUE);
+       
+               // Initialize number of tracks so far and logic switches
+               Int_t ntracks0 = esdTrackList->GetEntriesFast();
+               Bool_t signedTrack[kMaxTracksStack];
+               Bool_t fakeTrack[kMaxTracksStack];
+               for (Int_t i=0; i<ntracks; i++){
+                       signedTrack[i] = kFALSE;
+                       fakeTrack[i] = kFALSE;
+               }
+               //AliInfo("Selecting track candidates ...");
+               
+               // Sieve clusters in decreasing order of track quality
+               Double_t trackParams[7];
+//             AliTRDseedV1 *lseed = 0x0;
+               Int_t jSieve = 0, candidates;
+               do{
+                       //AliInfo(Form("\t\tITER = %i ", jSieve));
+
+                       // Check track candidates
+                       candidates = 0;
+                       for (Int_t itrack = 0; itrack < ntracks; itrack++) {
+                               Int_t trackIndex = sort[itrack];
+                               if (signedTrack[trackIndex] || fakeTrack[trackIndex]) continue;
+       
+                               
+                               // Calculate track parameters from tracklets seeds
+                               Int_t labelsall[1000];
+                               Int_t nlabelsall = 0;
+                               Int_t naccepted  = 0;
+                               Int_t ncl        = 0;
+                               Int_t nused      = 0;
+                               Int_t nlayers    = 0;
+                               Int_t findable   = 0;
+                               for (Int_t jLayer = 0; jLayer < kNPlanes; jLayer++) {
+                                       Int_t jseed = kNPlanes*trackIndex+jLayer;
+                                       if (TMath::Abs(sseed[jseed].GetYref(0) / sseed[jseed].GetX0()) < 0.15) 
+                                          findable++;
+       
+                                       if(!sseed[jseed].IsOK()) continue;
+                                       sseed[jseed].UpdateUsed();
+                                       ncl   += sseed[jseed].GetN2();
+                                       nused += sseed[jseed].GetNUsed();
+                                       nlayers++;
+       
+                                       // Cooking label
+                                       for (Int_t itime = 0; itime < nTimeBins; itime++) {
+                                               if(!sseed[jseed].IsUsable(itime)) continue;
+                                               naccepted++;
+                                               Int_t tindex = 0, ilab = 0;
+                                               while(ilab<3 && (tindex = sseed[jseed].GetClusters(itime)->GetLabel(ilab)) >= 0){
+                                                       labelsall[nlabelsall++] = tindex;
+                                                       ilab++;
+                                               }
+                                       }
+                               }
+                               // Filter duplicated tracks
+                               if (nused > 30){
+                                       //printf("Skip nused %d\n", nused);
+                                       fakeTrack[trackIndex] = kTRUE;
+                                       continue;
+                               }
+                               if (Float_t(nused)/ncl >= .25){
+                                       //printf("Skip nused/ncl >= .25\n");
+                                       fakeTrack[trackIndex] = kTRUE;
+                                       continue;
+                               }
+                               
+                               // Classify tracks
+                               Bool_t skip = kFALSE;
+                               switch(jSieve){
+                               case 0:
+                                       if(nlayers < 6) {skip = kTRUE; break;}
+                                       if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
+                                       break;
+       
+                               case 1:
+                                       if(nlayers < findable){skip = kTRUE; break;}
+                                       if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -4.){skip = kTRUE; break;}
+                                       break;
+       
+                               case 2:
+                                       if ((nlayers == findable) || (nlayers == 6)) { skip = kTRUE; break;}
+                                       if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -6.0){skip = kTRUE; break;}
+                                       break;
+       
+                               case 3:
+                                       if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;}
+                                       break;
+       
+                               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;}
+                                       break;
+                               }
+                               if(skip){
+                                       candidates++;
+                                       //printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused);
+                                       continue;
+                               }
+                               signedTrack[trackIndex] = kTRUE;
+                                               
+
+                               // Build track label - what happens if measured data ???
+                               Int_t labels[1000];
+                               Int_t outlab[1000];
+                               Int_t nlab = 0;
+                               for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+                                       Int_t jseed = kNPlanes*trackIndex+iLayer;
+                                       if(!sseed[jseed].IsOK()) continue;
+                                       for(int ilab=0; ilab<2; ilab++){
+                                               if(sseed[jseed].GetLabels(ilab) < 0) continue;
+                                               labels[nlab] = sseed[jseed].GetLabels(ilab);
+                                               nlab++;
+                                       }
+                               }
+                               Freq(nlab,labels,outlab,kFALSE);
+                               Int_t   label     = outlab[0];
+                               Int_t   frequency = outlab[1];
+                               Freq(nlabelsall,labelsall,outlab,kFALSE);
+                               Int_t   label1    = outlab[0];
+                               Int_t   label2    = outlab[2];
+                               Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted);
+       
+                               
+                               // Sign clusters
+                               AliTRDcluster *cl = 0x0; Int_t clusterIndex = -1;
+                               for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
+                                       Int_t jseed = kNPlanes*trackIndex+jLayer;
+                                       if(!sseed[jseed].IsOK()) continue;
+                                       if(TMath::Abs(sseed[jseed].GetYfit(1) - sseed[jseed].GetYfit(1)) >= .2) continue; // check this condition with Marian
+                                       sseed[jseed].UseClusters();
+                                       if(!cl){
+                                               Int_t ic = 0;
+                                               while(!(cl = sseed[jseed].GetClusters(ic))) ic++;
+                                               clusterIndex =  sseed[jseed].GetIndexes(ic);
+                                       }
+                               }
+                               if(!cl) continue;
+
+                               
+                               // Build track parameters
+                               AliTRDseedV1 *lseed =&sseed[trackIndex*6];
+                               Int_t idx = 0;
+                               while(idx<3 && !lseed->IsOK()) {
+                                       idx++;
+                                       lseed++;
+                               }
+                               Double_t cR = lseed->GetC();
+                               trackParams[1] = lseed->GetYref(0);
+                               trackParams[2] = lseed->GetZref(0);
+                               trackParams[3] = lseed->GetX0() * cR - TMath::Sin(TMath::ATan(lseed->GetYref(1)));
+                               trackParams[4] = lseed->GetZref(1) / TMath::Sqrt(1. + lseed->GetYref(1) * lseed->GetYref(1));
+                               trackParams[5] = cR;
+                               trackParams[0] = lseed->GetX0();
+                               trackParams[6] = layer[0].GetSector();/* *alpha+shift;  // Supermodule*/
+
+#ifdef DEBUG
+                               if(AliTRDReconstructor::StreamLevel() > 1) printf("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1]);
+                               
+                               if(AliTRDReconstructor::StreamLevel() >= 1){
+                                       Int_t sector = layer[0].GetSector();
+                                       Int_t nclusters = 0;
+                                       AliTRDseedV1 *dseed[6];
+                                       for(int is=0; is<6; is++){
+                                               dseed[is] = new AliTRDseedV1(sseed[trackIndex*6+is], kTRUE);
+                                               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 = *fDebugStreamer;
+                                       cstreamer << "Clusters2TracksStack"
+                                               << "Iter="      << fSieveSeeding
+                                               << "Like="      << fTrackQuality[trackIndex]
+                                               << "S0.="       << dseed[0]
+                                               << "S1.="       << dseed[1]
+                                               << "S2.="       << dseed[2]
+                                               << "S3.="       << dseed[3]
+                                               << "S4.="       << dseed[4]
+                                               << "S5.="       << dseed[5]
+                                               << "p0=" << trackParams[0]
+                                               << "p1=" << trackParams[1]
+                                               << "p2=" << trackParams[2]
+                                               << "p3=" << trackParams[3]
+                                               << "p4=" << trackParams[4]
+                                               << "p5=" << trackParams[5]
+                                               << "p6=" << trackParams[6]
+                                               << "Sector="    << sector
+                                               << "Stack="     << pars[1]
+                                               << "Label="     << label
+                                               << "Label1="    << label1
+                                               << "Label2="    << label2
+                                               << "FakeRatio=" << fakeratio
+                                               << "Freq="      << frequency
+                                               << "Ncl="       << ncl
+                                               << "NLayers="   << nlayers
+                                               << "Findable="  << findable
+                                               << "NUsed="     << nused
+                                               << "\n";
+                                       //???for(int is=0; is<6; is++) delete dseed[is];
+                               }
+#endif
+                       
+                               AliTRDtrack *track = AliTRDtrackerV1::RegisterSeed(&sseed[trackIndex*kNPlanes], trackParams);
+                               if(!track){
+                                       AliWarning("Fail to build a TRD Track.");
+                                       continue;
+                               }
+                               AliESDtrack esdTrack;
+                               esdTrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
+                               esdTrack.SetLabel(track->GetLabel());
+                               new ((*esdTrackList)[ntracks0++]) AliESDtrack(esdTrack);
+                               ntracks1++;
+                       }
+
+                       jSieve++;
+               } while(jSieve<5 && candidates); // end track candidates sieve
+               if(!ntracks1) break;
+
+               // increment counters
+               ntracks2 += ntracks1;
+               fSieveSeeding++;
+
+               // Rebuild plane configurations and indices taking only unused clusters into account
+               quality = BuildSeedingConfigs(layer, configs);
+               //if(quality < fRecoParam->GetPlaneQualityThreshold()) break;
+               
+               for(Int_t il = 0; il < kNPlanes * nTimeBins; 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));
+#endif
+       } while(fSieveSeeding<10); // end stack clusters sieve
+       
+
+
+       //AliInfo(Form("Registered TRD tracks %d in stack %d.", ntracks2, pars[1]));
+
+       return ntracks2;
+}
+
+//___________________________________________________________________
+Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDstackLayer *layers
+                                            , Int_t *configs)
+{
+  //
+  // Assign probabilities to chambers according to their
+  // capability of producing seeds.
+  // 
+  // Parameters :
+  //
+  //   layers : Array of stack propagation layers for all 6 chambers in one stack
+  //   configs : On exit array of configuration indexes (see GetSeedingConfig()
+  // for details) in the decreasing order of their seeding probabilities. 
+  //
+  // Output :
+  //
+  //  Return top configuration quality 
+  //
+  // Detailed description:
+  //
+  // To each chamber seeding configuration (see GetSeedingConfig() for
+  // the list of all configurations) one defines 2 quality factors:
+  //  - an apriori topological quality (see GetSeedingConfig() for details) and
+  //  - a data quality based on the uniformity of the distribution of
+  //    clusters over the x range (time bins population). See CookChamberQA() for details.
+  // 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]);
+               //printf("chamberQA[%d] = %f\n", iplane, chamberQA[iplane]);
+       }
+
+       Double_t tconfig[kNConfigs];
+       Int_t planes[4];
+       for(int iconf=0; iconf<kNConfigs; iconf++){
+               GetSeedingConfig(iconf, planes);
+               tconfig[iconf] = fTopologicQA[iconf];
+               for(int iplane=0; iplane<4; iplane++) tconfig[iconf] *= chamberQA[planes[iplane]]; 
+       }
+       
+       TMath::Sort(kNConfigs, tconfig, configs, kTRUE);
+       return tconfig[configs[0]];
+}
+
+//____________________________________________________________________
+Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
+                               , AliTRDseedV1 *sseed
+                               , Int_t *ipar)
+{
+  //
+  // Make tracklet seeds in the TRD stack.
+  //
+  // Parameters :
+  //   layers : Array of stack propagation layers containing clusters
+  //   sseed  : Array of empty tracklet seeds. On exit they are filled.
+  //   ipar   : Control parameters:
+  //       ipar[0] -> seeding chambers configuration
+  //       ipar[1] -> stack index
+  //       ipar[2] -> number of track candidates found so far
+  //
+  // Output :
+  //   Number of tracks candidates found.
+  // 
+  // Detailed description
+  //
+  // The following steps are performed:
+  // 1. Select seeding layers from seeding chambers
+  // 2. Select seeding clusters from the seeding AliTRDpropagationLayerStack.
+  //   The clusters are taken from layer 3, layer 0, layer 1 and layer 2, in
+  //   this order. The parameters controling the range of accepted clusters in
+  //   layer 0, 1, and 2 are defined in AliTRDstackLayer::BuildCond().
+  // 3. Helix fit of the cluster set. (see AliTRDtrackerFitter::FitRieman(AliTRDcluster**))
+  // 4. Initialize seeding tracklets in the seeding chambers.
+  // 5. Filter 0.
+  //   Chi2 in the Y direction less than threshold ... (1./(3. - sLayer))
+  //   Chi2 in the Z direction less than threshold ... (1./(3. - sLayer))
+  // 6. Attach clusters to seeding tracklets and find linear approximation of
+  //   the tracklet (see AliTRDseedV1::AttachClustersIter()). The number of used
+  //   clusters used by current seeds should not exceed ... (25).
+  // 7. Filter 1.
+  //   All 4 seeding tracklets should be correctly constructed (see
+  //   AliTRDseedV1::AttachClustersIter())
+  // 8. Helix fit of the seeding tracklets
+  // 9. Filter 2.
+  //   Likelihood calculation of the fit. (See AliTRDtrackerV1::CookLikelihood() for details)
+  // 10. Extrapolation of the helix fit to the other 2 chambers:
+  //    a) Initialization of extrapolation tracklet with fit parameters
+  //    b) Helix fit of tracklets
+  //    c) Attach clusters and linear interpolation to extrapolated tracklets
+  //    d) Helix fit of tracklets
+  // 11. Improve seeding tracklets quality by reassigning clusters.
+  //      See AliTRDtrackerV1::ImproveSeedQuality() for details.
+  // 12. Helix fit of all 6 seeding tracklets and chi2 calculation
+  // 13. Hyperplane fit and track quality calculation. See AliTRDtrackerFitter::FitHyperplane() for details.
+  // 14. Cooking labels for tracklets. Should be done only for MC
+  // 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
+       Int_t index[AliTRDstackLayer::kMaxClustersLayer], jndex[AliTRDstackLayer::kMaxClustersLayer];
+       // chi2 storage
+       // chi2[0] = tracklet chi2 on the Z direction
+       // chi2[1] = tracklet chi2 on the R direction
+       Double_t chi2[4];
+
+
+       // this should be data member of AliTRDtrack
+       Double_t seedQuality[kMaxTracksStack];
+       
+       // unpack control parameters
+       Int_t config  = ipar[0];
+       Int_t istack  = ipar[1];
+       Int_t ntracks = ipar[2];
+       Int_t planes[kNSeedPlanes]; GetSeedingConfig(config, planes);   
+#ifdef DEBUG
+               if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks));
+#endif
+       
+       // Init chambers 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();
+       }
+
+       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]);
+
+
+       // Start finding seeds
+       Int_t icl = 0;
+       while((c[3] = (*layer[3])[icl++])){
+               if(!c[3]) continue;
+               layer[0]->BuildCond(c[3], cond0, 0);
+               layer[0]->GetClusters(cond0, index, ncl);
+               Int_t jcl = 0;
+               while(jcl<ncl) {
+                       c[0] = (*layer[0])[index[jcl++]];
+                       if(!c[0]) continue;
+                       Double_t dx    = c[3]->GetX() - c[0]->GetX();
+                       Double_t theta = (c[3]->GetZ() - c[0]->GetZ())/dx;
+                       Double_t phi   = (c[3]->GetY() - c[0]->GetY())/dx;
+                       layer[1]->BuildCond(c[0], cond1, 1, theta, phi);
+                       layer[1]->GetClusters(cond1, jndex, mcl);
+
+                       Int_t kcl = 0;
+                       while(kcl<mcl) {
+                               c[1] = (*layer[1])[jndex[kcl++]];
+                               if(!c[1]) continue;
+                               layer[2]->BuildCond(c[1], cond2, 2, theta, phi);
+                               c[2] = layer[2]->GetNearestCluster(cond2);
+                               if(!c[2]) continue;
+                               
+                               //AliInfo("Seeding clusters found. Building seeds ...");
+                               //for(Int_t i = 0; i < kNSeedPlanes; i++) printf("%i. coordinates: x = %3.3f, y = %3.3f, z = %3.3f\n", i, c[i]->GetX(), c[i]->GetY(), c[i]->GetZ());
+                               for (Int_t il = 0; il < 6; il++) cseed[il].Reset();
+
+                               fFitter->Reset();
+
+                               fFitter->FitRieman(c, kNSeedPlanes);
+
+                               chi2[0] = 0.; chi2[1] = 0.;
+                               AliTRDseedV1 *tseed = 0x0;
+                               for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
+                                       tseed = &cseed[planes[iLayer]];
+                                       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());
+                                       chi2[0] += tseed->GetChi2Z(c[iLayer]->GetZ());
+                                       chi2[1] += tseed->GetChi2Y(c[iLayer]->GetY());
+                               }
+
+                               Bool_t isFake = kFALSE;
+                               if (c[0]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
+                               if (c[1]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
+                               if (c[2]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
+#ifdef DEBUG
+                               if(AliTRDReconstructor::StreamLevel() >= 2){
+                                       Float_t yref[4], ycluster[4];
+                                       for(int il=0; il<4; il++){
+                                               tseed = &cseed[planes[il]];
+                                               yref[il] = tseed->GetYref(0);
+                                               ycluster[il] = c[il]->GetY();
+                                       }
+                                       Float_t threshold = .5;//1./(3. - sLayer);
+                                       Int_t ll = c[3]->GetLabel(0);
+                                       TTreeSRedirector &cs0 = *fDebugStreamer;
+                                                       cs0 << "MakeSeeds0"
+                                                       <<"isFake=" << isFake
+                                                       <<"label=" << ll
+                                                       <<"threshold=" << threshold
+                                                       <<"chi2=" << chi2[1]
+                                                       <<"yref0="<<yref[0]
+                                                       <<"yref1="<<yref[1]
+                                                       <<"yref2="<<yref[2]
+                                                       <<"yref3="<<yref[3]
+                                                       <<"ycluster0="<<ycluster[0]
+                                                       <<"ycluster1="<<ycluster[1]
+                                                       <<"ycluster2="<<ycluster[2]
+                                                       <<"ycluster3="<<ycluster[3]
+                                                       <<"\n";
+                               }
+#endif
+
+                               if(chi2[0] > fRecoParam->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){
+                                       //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0]));
+                                       continue;
+                               }
+                               if(chi2[1] > fRecoParam->GetChi2Y()/*1./(3. - sLayer)*//*iter*/){
+                                       //AliInfo(Form("Failed chi2 filter on chi2Y [%f].", chi2[1]));
+                                       continue;
+                               }
+                               //AliInfo("Passed chi2 filter.");
+
+#ifdef DEBUG
+                               if(AliTRDReconstructor::StreamLevel() >= 2){
+                                       Float_t minmax[2] = { -100.0,  100.0 };
+                                       for (Int_t iLayer = 0; iLayer < 4; iLayer++) {
+                                               Float_t max = c[iLayer]->GetZ() + cseed[planes[iLayer]].GetPadLength() * 0.5 + 1.0 - cseed[planes[iLayer]].GetZref(0);
+                                               if (max < minmax[1]) minmax[1] = max;
+                                               Float_t min = c[iLayer]->GetZ()-cseed[planes[iLayer]].GetPadLength() * 0.5 - 1.0 - cseed[planes[iLayer]].GetZref(0);
+                                               if (min > minmax[0]) minmax[0] = min;
+                                       }
+                                       Double_t xpos[4];
+                                       for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX();
+                                       TTreeSRedirector &cstreamer = *fDebugStreamer;
+                                                       cstreamer << "MakeSeeds1"
+                                               << "isFake=" << isFake
+                                               << "config="   << config
+                                               << "Cl0.="   << c[0]
+                                               << "Cl1.="   << c[1]
+                                               << "Cl2.="   << c[2]
+                                               << "Cl3.="   << c[3]
+                                               << "X0="     << xpos[0] //layer[sLayer]->GetX()
+                                               << "X1="     << xpos[1] //layer[sLayer + 1]->GetX()
+                                               << "X2="     << xpos[2] //layer[sLayer + 2]->GetX()
+                                               << "X3="     << xpos[3] //layer[sLayer + 3]->GetX()
+                                               << "Y2exp="  << cond2[0]
+                                               << "Z2exp="  << cond2[1]
+                                               << "Chi2R="  << chi2[0]
+                                               << "Chi2Z="  << chi2[1]
+                                               << "Seed0.=" << &cseed[planes[0]]
+                                               << "Seed1.=" << &cseed[planes[1]]
+                                               << "Seed2.=" << &cseed[planes[2]]
+                                               << "Seed3.=" << &cseed[planes[3]]
+                                               << "Zmin="   << minmax[0]
+                                               << "Zmax="   << minmax[1]
+                                               << "\n" ;
+                               }               
+#endif
+                               // try attaching clusters to tracklets
+                               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();
+                                       if(nUsedCl > 25) break;
+                                       nlayers++;
+                               }
+                               if(nlayers < kNSeedPlanes){ 
+                                       //AliInfo("Failed updating all seeds.");
+                                       continue;
+                               }
+                               // fit tracklets and cook likelihood
+                               chi2[0] = 0.; chi2[1] = 0.;
+                               fFitter->FitRieman(&cseed[0], &planes[0]);
+                               AliRieman *rim = fFitter->GetRiemanFitter();
+                               for(int iLayer=0; iLayer<4; iLayer++){
+                                       cseed[planes[iLayer]].Update(rim);
+                                       chi2[0] += (Float_t)cseed[planes[iLayer]].GetChi2Z();
+                                       chi2[1] += cseed[planes[iLayer]].GetChi2Y();
+                               }
+                               Double_t chi2r = chi2[1], chi2z = chi2[0];
+                               Double_t like = CookLikelihood(&cseed[0], planes, chi2); // to be checked
+                               if (TMath::Log(1.E-9 + like) < fRecoParam->GetTrackLikelihood()){
+                                       //AliInfo(Form("Failed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
+                                       continue;
+                               }
+                               //AliInfo(Form("Passed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
+
+
+                               // book preliminary results
+                               seedQuality[ntracks] = like;
+                               fSeedLayer[ntracks]  = config;/*sLayer;*/
+
+                               // attach clusters to the extrapolation seeds
+                               Int_t lextrap[2];
+                               GetExtrapolationConfig(config, lextrap);
+                               Int_t nusedf   = 0; // debug value
+                               for(int iLayer=0; iLayer<2; iLayer++){
+                                       Int_t jLayer = lextrap[iLayer];
+                                       
+                                       // prepare extrapolated seed
+                                       cseed[jLayer].Reset();
+                                       cseed[jLayer].SetRecoParam(fRecoParam);
+                                       cseed[jLayer].SetLayer(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
+
+                                       // 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;
+                                       cseed[jLayer] = tseed;
+                                       nusedf += cseed[jLayer].GetNUsed(); // debug value
+                                       AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
+                               }
+                               //AliInfo("Extrapolation done.");
+
+                               ImproveSeedQuality(layers, cseed);
+                               //AliInfo("Improve seed quality done.");
+
+                               nlayers   = 0;
+                               Int_t nclusters = 0;
+                               Int_t findable  = 0;
+                               for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+                                       if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) findable++;
+                                       if (!cseed[iLayer].IsOK()) continue;
+                                       nclusters += cseed[iLayer].GetN2();
+                                       nlayers++;
+                               }
+                               if (nlayers < 3){ 
+                                       //AliInfo("Failed quality check on seeds.");
+                                       continue;
+                               }
+
+                               // fit full track and cook likelihoods
+                               fFitter->FitRieman(&cseed[0]);
+                               Double_t chi2ZF = 0., chi2RF = 0.;
+                               for(int ilayer=0; ilayer<6; ilayer++){
+                                       cseed[ilayer].Update(fFitter->GetRiemanFitter());
+                                       if (!cseed[ilayer].IsOK()) continue;
+                                       //tchi2 = cseed[ilayer].GetChi2Z();
+                                       //printf("layer %d chi2 %e\n", ilayer, tchi2);
+                                       chi2ZF += cseed[ilayer].GetChi2Z();
+                                       chi2RF += cseed[ilayer].GetChi2Y();
+                               }
+                               chi2ZF /= TMath::Max((nlayers - 3.), 1.);
+                               chi2RF /= TMath::Max((nlayers - 3.), 1.);
+
+                               // do the final track fitting
+                               fFitter->SetLayers(nlayers);
+#ifdef DEBUG
+                               fFitter->SetDebugStream(fDebugStreamer);
+#endif
+                               fTrackQuality[ntracks] = fFitter->FitHyperplane(&cseed[0], chi2ZF, GetZ());
+                               Double_t param[3];
+                               Double_t chi2[2];
+                               fFitter->GetHyperplaneFitResults(param);
+                               fFitter->GetHyperplaneFitChi2(chi2);
+                               //AliInfo("Hyperplane fit done\n");
+
+                               // finalize tracklets
+                               Int_t labels[12];
+                               Int_t outlab[24];
+                               Int_t nlab = 0;
+                               for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+                                       if (!cseed[iLayer].IsOK()) continue;
+
+                                       if (cseed[iLayer].GetLabels(0) >= 0) {
+                                               labels[nlab] = cseed[iLayer].GetLabels(0);
+                                               nlab++;
+                                       }
+
+                                       if (cseed[iLayer].GetLabels(1) >= 0) {
+                                               labels[nlab] = cseed[iLayer].GetLabels(1);
+                                               nlab++;
+                                       }
+                               }
+                               Freq(nlab,labels,outlab,kFALSE);
+                               Int_t label     = outlab[0];
+                               Int_t frequency = outlab[1];
+                               for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
+                                       cseed[iLayer].SetFreq(frequency);
+                                       cseed[iLayer].SetC(param[1]/*cR*/);
+                                       cseed[iLayer].SetCC(param[0]/*cC*/);
+                                       cseed[iLayer].SetChi2(chi2[0]);
+                                       cseed[iLayer].SetChi2Z(chi2ZF);
+                               }
+           
+#ifdef DEBUG
+                               if(AliTRDReconstructor::StreamLevel() >= 2){
+                                       Double_t curv = (fFitter->GetRiemanFitter())->GetC();
+                                       TTreeSRedirector &cstreamer = *fDebugStreamer;
+                                       cstreamer << "MakeSeeds2"
+                                               << "C="       << curv
+                                               << "Chi2R="   << chi2r
+                                               << "Chi2Z="   << chi2z
+                                               << "Chi2TR="  << chi2[0]
+                                               << "Chi2TC="  << chi2[1]
+                                               << "Chi2RF="  << chi2RF
+                                               << "Chi2ZF="  << chi2ZF
+                                               << "Ncl="     << nclusters
+                                               << "Nlayers=" << nlayers
+                                               << "NUsedS="  << nUsedCl
+                                               << "NUsed="   << nusedf
+                                               << "Findable" << findable
+                                               << "Like="    << like
+                                               << "S0.="     << &cseed[0]
+                                               << "S1.="     << &cseed[1]
+                                               << "S2.="     << &cseed[2]
+                                               << "S3.="     << &cseed[3]
+                                               << "S4.="     << &cseed[4]
+                                               << "S5.="     << &cseed[5]
+                                               << "Label="   << label
+                                               << "Freq="    << frequency
+                                               << "\n";
+                               }
+#endif
+                               
+                               ntracks++;
+                               if(ntracks == kMaxTracksStack){
+                                       AliWarning(Form("Number of seeds reached maximum allowed (%d) in stack.", kMaxTracksStack));
+                                       return ntracks;
+                               }
+                               cseed += 6;
+                       }
+               }
+       }
+       for(int isl=0; isl<4; isl++) delete layer[isl];
+       
+       return ntracks;
+}
+
+//_____________________________________________________________________________
+AliTRDtrack *AliTRDtrackerV1::RegisterSeed(AliTRDseedV1 *seeds, Double_t *params)
+{
+  //
+  // Build a TRD track out of tracklet candidates
+  //
+  // Parameters :
+  //   seeds  : array of tracklets
+  //   params : track parameters (see MakeSeeds() function body for a detailed description)
+  //
+  // Output :
+  //   The TRD track.
+  //
+  // Detailed description
+  //
+  // 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];
+
+  c[ 0] = 0.2;
+  c[ 1] = 0.0; c[ 2] = 2.0;
+  c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02;
+  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
+       track->PropagateTo(params[0]-5.0);
+  track->ResetCovariance(1);
+  Int_t rc = FollowBackProlongation(*track);
+  if (rc < 30) {
+    delete track;
+    track = 0;
+  }
+  else {
+    track->CookdEdx();
+    track->CookdEdxTimBin(-1);
+    CookLabel(track,0.9);
+  }
+
+  return track;
+}
+
+//____________________________________________________________________
+void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers
+                                       , AliTRDseedV1 *cseed)
+{
+  //
+  // Sort tracklets according to "quality" and try to "improve" the first 4 worst
+  //
+  // Parameters :
+  //  layers : Array of propagation layers for a stack/supermodule
+  //  cseed  : Array of 6 seeding tracklets which has to be improved
+  // 
+  // Output :
+  //   cssed : Improved seeds
+  // 
+  // Detailed description
+  //
+  // Iterative procedure in which new clusters are searched for each
+  // tracklet seed such that the seed quality (see AliTRDseed::GetQuality())
+  // can be maximized. If some optimization is found the old seeds are replaced.
+  //
+       
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       Int_t nTimeBins = cal->GetNumberOfTimeBins();
+       
+       // make a local working copy
+       AliTRDseedV1 bseed[6];
+       for (Int_t jLayer = 0; jLayer < 6; jLayer++) bseed[jLayer] = cseed[jLayer];
+
+
+       Float_t lastquality = 10000.0;
+       Float_t lastchi2    = 10000.0;
+       Float_t chi2        =  1000.0;
+
+       for (Int_t iter = 0; iter < 4; iter++) {
+               Float_t sumquality = 0.0;
+               Float_t squality[6];
+               Int_t   sortindexes[6];
+
+               for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
+                       squality[jLayer]  = bseed[jLayer].IsOK() ? bseed[jLayer].GetQuality(kTRUE) : -1.;
+                       sumquality +=squality[jLayer];
+               }
+               if ((sumquality >= lastquality) || (chi2       >     lastchi2)) break;
+
+               
+               lastquality = sumquality;
+               lastchi2    = chi2;
+               if (iter > 0) for (Int_t jLayer = 0; jLayer < 6; jLayer++) cseed[jLayer] = bseed[jLayer];
+
+               
+               TMath::Sort(6, squality, sortindexes, kFALSE);
+               for (Int_t jLayer = 5; jLayer > 1; jLayer--) {
+                       Int_t bLayer = sortindexes[jLayer];
+                       bseed[bLayer].AttachClustersIter(&layers[bLayer*nTimeBins], squality[bLayer], kTRUE);
+               }
+
+               chi2 = AliTRDseedV1::FitRiemanTilt(bseed,kTRUE);
+       } // Loop: iter
+}
+
+//____________________________________________________________________
+Double_t  AliTRDtrackerV1::CookPlaneQA(AliTRDstackLayer *layers)
+{
+  //
+  // Calculate plane quality for seeding.
+  // 
+  //
+  // Parameters :
+  //   layers : Array of propagation layers for this plane.
+  //
+  // Output :
+  //   plane quality factor for seeding
+  // 
+  // Detailed description
+  //
+  // The quality of the plane for seeding is higher if:
+  //  1. the average timebin population is closer to an integer number
+  //  2. the distribution of clusters/timebin is closer to a uniform distribution.
+  //    - the slope of the first derivative of a parabolic fit is small or
+  //    - the slope of a linear fit is small
+  //
+
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+//     Double_t x;
+//     TLinearFitter fitter(1, "pol1");
+//     fitter.ClearPoints();
+       Int_t ncl = 0;
+       Int_t nused = 0;
+       Int_t nClLayer;
+       for(int itb=0; itb<nTimeBins; itb++){
+               //x = layer[itb].GetX();
+               //printf("x[%d] = %f nCls %d\n", itb, x, layer[itb].GetNClusters());
+               //if(!layer[itb].GetNClusters()) continue;
+               //fitter.AddPoint(&x, layer[itb].GetNClusters(), 1.);
+               nClLayer = layers[itb].GetNClusters();
+               ncl += nClLayer;
+               for(Int_t incl = 0; incl < nClLayer; incl++)
+                       if((layers[itb].GetCluster(incl))->IsUsed()) nused++;
+       }
+       
+       // calculate the deviation of the mean number of clusters from the
+       // closest integer values
+       Float_t ncl_med = float(ncl-nused)/nTimeBins;
+       Int_t ncli = Int_t(ncl_med);
+       Float_t ncl_dev = TMath::Abs(ncl_med - TMath::Max(ncli, 1));
+       ncl_dev -= (ncl_dev>.5) && ncli ? .5 : 0.; 
+       /*Double_t quality = */ return TMath::Exp(-2.*ncl_dev);
+       
+//     // get slope of the derivative
+//     if(!fitter.Eval()) return quality;
+//     fitter.PrintResults(3);
+//     Double_t a = fitter.GetParameter(1);
+// 
+//     printf("ncl_dev(%f)  a(%f)\n", ncl_dev, a);
+//     return quality*TMath::Exp(-a);
+}
+
+//____________________________________________________________________
+Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed
+                                       , Int_t planes[4]
+                                       , Double_t *chi2)
+{
+  //
+  // Calculate the probability of this track candidate.
+  //
+  // Parameters :
+  //   cseeds : array of candidate tracklets
+  //   planes : array of seeding planes (see seeding configuration)
+  //   chi2   : chi2 values (on the Z and Y direction) from the rieman fit of the track.
+  //
+  // Output :
+  //   likelihood value
+  // 
+  // Detailed description
+  //
+  // The track quality is estimated based on the following 4 criteria:
+  //  1. precision of the rieman fit on the Y direction (likea)
+  //  2. chi2 on the Y direction (likechi2y)
+  //  3. chi2 on the Z direction (likechi2z)
+  //  4. number of attached clusters compared to a reference value 
+  //     (see AliTRDrecoParam::fkFindable) (likeN)
+  //
+  // The distributions for each type of probabilities are given below as of
+  // (date). They have to be checked to assure consistency of estimation.
+  //
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       Int_t nTimeBins = cal->GetNumberOfTimeBins();
+       // ratio of the total number of clusters/track which are expected to be found by the tracker.
+       Float_t fgFindable = fRecoParam->GetFindableClusters();
+
+       
+       Int_t nclusters = 0;
+       Double_t sumda = 0.;
+       for(UChar_t ilayer = 0; ilayer < 4; ilayer++){
+               Int_t jlayer = planes[ilayer];
+               nclusters += cseed[jlayer].GetN2();
+               sumda += TMath::Abs(cseed[jlayer].GetYfitR(1) - cseed[jlayer].GetYref(1));
+       }
+       Double_t likea     = TMath::Exp(-sumda*10.6);
+       Double_t likechi2y  = 0.0000000001;
+       if (chi2[1] < 0.5) likechi2y += TMath::Exp(-TMath::Sqrt(chi2[1]) * 7.73);
+       Double_t likechi2z = TMath::Exp(-chi2[0] * 0.088) / TMath::Exp(-chi2[0] * 0.019);
+       Int_t enc = Int_t(fgFindable*4.*nTimeBins);     // Expected Number Of Clusters, normally 72
+       Double_t likeN     = TMath::Exp(-(enc - nclusters) * 0.19);
+       
+       Double_t like      = likea * likechi2y * likechi2z * likeN;
+
+#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 = *fDebugStreamer;
+               cstreamer << "CookLikelihood"
+                       << "sumda="     << sumda
+                       << "chi0="      << chi2[0]
+                       << "chi1="      << chi2[1]
+                       << "likea="     << likea
+                       << "likechi2y=" << likechi2y
+                       << "likechi2z=" << likechi2z
+                       << "nclusters=" << nclusters
+                       << "likeN="     << likeN
+                       << "like="      << like
+                       << "\n";
+       }
+#endif
+
+       return like;
+}
+
+//___________________________________________________________________
+void AliTRDtrackerV1::GetMeanCLStack(AliTRDstackLayer *layers
+                                   , Int_t *planes
+                                   , Double_t *params)
+{
+  //
+  // Determines the Mean number of clusters per layer.
+  // Needed to determine good Seeding Layers
+  //
+  // Parameters:
+  //    - Array of AliTRDstackLayers
+  //    - Container for the params
+  //
+  // Detailed description
+  //
+  // Two Iterations:
+  // In the first Iteration the mean is calculted using all layers.
+  // After this, all layers outside the 1-sigma-region are rejected.
+  // Then the mean value and the standard-deviation are calculted a second
+  // time in order to select all layers in the 1-sigma-region as good-candidates.
+  //
+
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       Int_t nTimeBins = cal->GetNumberOfTimeBins();
+       
+       Float_t mean = 0, stdev = 0;
+       Double_t ncl[kNTimeBins*kNSeedPlanes], mcl[kNTimeBins*kNSeedPlanes];
+       Int_t position = 0;
+       memset(ncl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
+       memset(mcl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes);
+       Int_t nused = 0;
+       for(Int_t ipl = 0; ipl < kNSeedPlanes; ipl++){
+               for(Int_t ils = 0; ils < nTimeBins; ils++){
+                       position = planes[ipl]*nTimeBins + ils;
+                       ncl[ipl * nTimeBins + ils] = layers[position].GetNClusters();
+                       nused = 0;
+                       for(Int_t icl = 0; icl < ncl[ipl * nTimeBins + ils]; icl++)
+                               if((layers[position].GetCluster(icl))->IsUsed()) nused++;
+                       ncl[ipl * nTimeBins + ils] -= nused;
+               }
+       }
+       // Declaration of quartils:
+       //Double_t qvals[3] = {0.0, 0.0, 0.0};
+       //Double_t qprop[3] = {0.16667, 0.5, 0.83333};
+       // Iterations
+       Int_t counter;
+       Double_t *array;
+       Int_t *limit;
+       Int_t nLayers = nTimeBins * kNSeedPlanes;
+       for(Int_t iter = 0; iter < 2; iter++){
+               array = (iter == 0) ? &ncl[0] : &mcl[0];
+               limit = (iter == 0) ? &nLayers : &counter;
+               counter = 0;
+               if(iter == 1){
+                       for(Int_t i = 0; i < nTimeBins *kNSeedPlanes; i++){
+                               if((ncl[i] >  mean + stdev) || (ncl[i] <  mean - stdev)) continue; // Outside 1-sigma region
+//                             if((ncl[i] >  qvals[2]) || (ncl[i] <  qvals[0])) continue; // Outside 1-sigma region
+                               if(ncl[i] == 0) continue;                                                // 0-Layers also rejected
+                               mcl[counter] = ncl[i];
+                               counter++;
+                       }
+               }
+               if(*limit == 0) break;
+               printf("Limit = %d\n", *limit);
+               //using quartils instead of mean and RMS 
+//             TMath::Quantiles(*limit,3,array,qvals,qprop,kFALSE);
+               mean = TMath::Median(*limit, array, 0x0);
+               stdev  = TMath::RMS(*limit, array);
+       }
+//     printf("Quantiles: 0.16667 = %3.3f, 0.5 = %3.3f, 0.83333 = %3.3f\n", qvals[0],qvals[1],qvals[2]);
+//     memcpy(params,qvals,sizeof(Double_t)*3);
+       params[1] = (Double_t)TMath::Nint(mean);
+       params[0] = (Double_t)TMath::Nint(mean - stdev);
+       params[2] = (Double_t)TMath::Nint(mean + stdev);
+
+}
+
+//___________________________________________________________________
+Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers
+                                      , Double_t *params)
+{
+  //
+  // Algorithm to find optimal seeding layer
+  // Layers inside one sigma region (given by Quantiles) are sorted
+  // according to their difference.
+  // All layers outside are sorted according t
+  //
+  // Parameters:
+  //     - Array of AliTRDstackLayers (in the current plane !!!)
+  //     - Container for the Indices of the seeding Layer candidates
+  //
+  // Output:
+  //     - Number of Layers inside the 1-sigma-region
+  //
+  // The optimal seeding layer should contain the mean number of
+  // custers in the layers in one chamber.
+  //
+
+       //printf("Params: %3.3f, %3.3f, %3.3f\n", params[0], params[1], params[2]);
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       const Int_t kMaxClustersLayer = AliTRDstackLayer::kMaxClustersLayer;
+       Int_t nTimeBins = cal->GetNumberOfTimeBins();
+       Int_t ncl[kNTimeBins], indices[kNTimeBins], bins[kMaxClustersLayer];
+       memset(ncl, 0, sizeof(Int_t)*kNTimeBins);
+       memset(indices, 0, sizeof(Int_t)*kNTimeBins);
+       memset(bins, 0, sizeof(Int_t)*kMaxClustersLayer);
+       Int_t nused = 0;
+       for(Int_t ils = 0; ils < nTimeBins; ils++){
+               ncl[ils] = layers[ils].GetNClusters();
+               nused = 0;
+               for(Int_t icl = 0; icl < ncl[ils]; icl++)
+                       if((layers[ils].GetCluster(icl))->IsUsed()) nused++;
+               ncl[ils] -= nused;
+       }
+       
+       Float_t mean = params[1];
+       for(Int_t ils = 0; ils < nTimeBins; ils++){
+               memmove(indices + bins[ncl[ils]+1] + 1, indices + bins[ncl[ils]+1], sizeof(Int_t)*(nTimeBins - ils));
+               indices[bins[ncl[ils]+1]] = ils;
+               for(Int_t i = ncl[ils]+1; i < kMaxClustersLayer; i++)
+                       bins[i]++;
+       }
+       
+       //for(Int_t i = 0; i < nTimeBins; i++) printf("Bin %d = %d\n", i, bins[i]);
+       Int_t sbin = -1;
+       Int_t nElements;
+       Int_t position = 0;
+       TRandom *r = new TRandom();
+       Int_t iter = 0;
+       while(1){
+               while(sbin < (Int_t)params[0] || sbin > (Int_t)params[2]){
+                       // Randomly selecting one bin
+                       sbin = (Int_t)r->Poisson(mean);
+               }
+               printf("Bin = %d\n",sbin);
+               //Randomly selecting one Layer in the bin
+               nElements = bins[sbin + 1] - bins[sbin];
+               printf("nElements = %d\n", nElements);
+               if(iter == 5){
+                       position = (Int_t)(gRandom->Rndm()*(nTimeBins-1));
+                       break;
+               }
+               else if(nElements==0){
+                       iter++;
+                       continue;
+               }
+               position = (Int_t)(gRandom->Rndm()*(nElements-1)) + bins[sbin];
+               break;
+       }
+       delete r;
+       return indices[position];
+}
+
+//____________________________________________________________________
+AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDstackLayer *layers
+                                                 , AliTRDseedV1/*AliRieman*/ *reference)
+{
+  //
+  // Finds a seeding Cluster for the extrapolation chamber.
+  //
+  // The seeding cluster should be as close as possible to the assumed
+  // track which is represented by a Rieman fit.
+  // Therefore the selecting criterion is the minimum distance between
+  // the best fitting cluster and the Reference which is derived from
+  // the AliTRDseed. Because all layers are assumed to be equally good
+  // a linear search is performed.
+  //
+  // Imput parameters: - layers: array of AliTRDstackLayers (in one chamber!!!)
+  //                   - sfit: the reference
+  //
+  // Output:           - the best seeding cluster
+  //
+
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       Int_t nTimeBins = cal->GetNumberOfTimeBins();
+       
+       // distances as squared distances
+       Int_t index = 0;
+       Float_t ypos = 0.0, zpos = 0.0, distance = 0.0, nearest_distance =100000.0; 
+       ypos = reference->GetYref(0);
+       zpos = reference->GetZref(0);
+       AliTRDcluster *currentBest = 0x0, *temp = 0x0;
+       for(Int_t ils = 0; ils < nTimeBins; ils++){
+               // Reference positions
+//             ypos = reference->GetYat(layers[ils].GetX());
+//             zpos = reference->GetZat(layers[ils].GetX());
+               index = layers[ils].SearchNearestCluster(ypos, zpos, fRecoParam->GetRoad2y(), fRecoParam->GetRoad2z());
+               if(index == -1) continue;
+               temp = layers[ils].GetCluster(index);
+               if(!temp) continue;
+               distance = (temp->GetY() - ypos) * (temp->GetY() - ypos) + (temp->GetZ() - zpos) * (temp->GetZ() - zpos);
+               if(distance < nearest_distance){
+                       nearest_distance = distance;
+                       currentBest = temp;
+               }
+       }
+       return currentBest;
+}
+
+//____________________________________________________________________
+AliTRDstackLayer *AliTRDtrackerV1::MakeSeedingLayer(AliTRDstackLayer *layers
+                                                  , Int_t plane)
+{
+  //
+  // Creates a seeding layer
+  //
+       
+       // constants
+       const Int_t kMaxRows = 16;
+       const Int_t kMaxCols = 144;
+       const Int_t kMaxPads = 2304;
+       
+       // Get the calculation
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       Int_t nTimeBins = cal->GetNumberOfTimeBins();
+       
+       // Get the geometrical data of the chamber
+       AliTRDpadPlane *pp = fGeom->GetPadPlane(plane, layers[0].GetStackNr());
+       Int_t nCols = pp->GetNcols();
+       Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd());
+       Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd());
+       Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd());
+       Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd());
+       Int_t nRows = pp->GetNrows();
+       Float_t binlength = (ymax - ymin)/nCols; 
+       //AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength));
+       
+       // Fill the histogram
+       Int_t arrpos;
+       Float_t ypos;
+       Int_t irow, nClusters;  
+       Int_t *histogram[kMaxRows];                                                                                     // 2D-Histogram
+       Int_t hvals[kMaxPads];  memset(hvals, 0, sizeof(Int_t)*kMaxPads);       
+       Float_t *sigmas[kMaxRows];
+       Float_t svals[kMaxPads];        memset(svals, 0, sizeof(Float_t)*kMaxPads);     
+       AliTRDcluster *c = 0x0;
+       for(Int_t irs = 0; irs < kMaxRows; irs++){
+               histogram[irs] = &hvals[irs*kMaxCols];
+               sigmas[irs] = &svals[irs*kMaxCols];
+       }
+       for(Int_t iTime = 0; iTime < nTimeBins; iTime++){
+               nClusters = layers[iTime].GetNClusters();
+               for(Int_t incl = 0; incl < nClusters; incl++){
+                       c = layers[iTime].GetCluster(incl);     
+                       ypos = c->GetY();
+                       if(ypos > ymax && ypos < ymin) continue;
+                       irow = pp->GetPadRowNumber(c->GetZ());                          // Zbin
+                       if(irow < 0)continue;
+                       arrpos = static_cast<Int_t>((ypos - ymin)/binlength);
+                       if(ypos == ymax) arrpos = nCols - 1;
+                       histogram[irow][arrpos]++;
+                       sigmas[irow][arrpos] += c->GetSigmaZ2();
+               }
+       }
+       
+// Now I have everything in the histogram, do the selection
+//     printf("Starting the analysis\n");
+       //Int_t nPads = nCols * nRows;
+       // This is what we are interested in: The center of gravity of the best candidates
+       Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads);
+       Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads);
+       Float_t *cogy[kMaxRows];
+       Float_t *cogz[kMaxRows];
+       // Lookup-Table storing coordinates according ti the bins
+       Float_t yLengths[kMaxCols];
+       Float_t zLengths[kMaxRows];
+       for(Int_t icnt = 0; icnt < nCols; icnt++){
+               yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2;
+       }
+       for(Int_t icnt = 0; icnt < nRows; icnt++){
+               zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2;
+       }
+
+       // A bitfield is used to mask the pads as usable
+       Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads];
+       for(UChar_t icount = 0; icount < nRows; icount++){
+               cogy[icount] = &cogyvals[icount*kMaxCols];
+               cogz[icount] = &cogzvals[icount*kMaxCols];
+       }
+       // In this array the array position of the best candidates will be stored
+       Int_t cand[kMaxTracksStack];
+       Float_t sigcands[kMaxTracksStack];
+       
+       // helper variables
+       Int_t indices[kMaxPads]; memset(indices, 0, sizeof(Int_t)*kMaxPads);
+       Int_t nCandidates = 0;
+       Float_t norm, cogv;
+       // histogram filled -> Select best bins
+       TMath::Sort(kMaxPads, hvals, indices);                  // bins storing a 0 should not matter
+       // Set Threshold
+       Int_t maximum = hvals[indices[0]];      // best
+       Int_t threshold = static_cast<UChar_t>(maximum * fRecoParam->GetFindableClusters());
+       Int_t col, row, lower, lower1, upper, upper1;
+       for(Int_t ib = 0; ib < kMaxPads; ib++){
+               if(nCandidates >= kMaxTracksStack){
+                       AliWarning(Form("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, kMaxTracksStack));
+                       break;
+               }
+               // Positions
+               row = indices[ib]/nCols;
+               col = indices[ib]%nCols;
+               // here will be the threshold condition:
+               if((mask[col] & (1 << row)) != 0) continue;             // Pad is masked: continue
+               if(histogram[row][col] < TMath::Max(threshold, 1)){     // of course at least one cluster is needed
+                       break;                  // number of clusters below threshold: break;
+               } 
+               // passing: Mark the neighbors
+               lower  = TMath::Max(col - 1, 0); upper  = TMath::Min(col + 2, nCols);
+               lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols);
+               for(Int_t ic = lower; ic < upper; ++ic)
+                       for(Int_t ir = lower1; ir < upper1; ++ir){
+                               if(ic == col && ir == row) continue;
+                               mask[ic] |= (1 << ir);
+                       }
+               // Storing the position in an array
+               // testing for neigboring
+               cogv = 0;
+               norm = 0;
+               lower = TMath::Max(col - 1,0);
+               upper = TMath::Min(col + 2, nCols);
+               for(Int_t inb = lower; inb < upper; ++inb){
+                       cogv += yLengths[inb] * histogram[row][inb];
+                       norm += histogram[row][inb];
+               }
+               cogy[row][col] = cogv / norm;
+               cogv = 0; norm = 0;
+               lower = TMath::Max(row - 1, 0);
+               upper = TMath::Min(row + 2, nRows);
+               for(Int_t inb = lower; inb < upper; ++inb){
+                       cogv += zLengths[inb] * histogram[inb][col];
+                       norm += histogram[inb][col];
+               }
+               cogz[row][col] = cogv /  norm;
+               // passed the filter
+               cand[nCandidates] = row*kMaxCols + col; // store the position of a passig candidate into an Array
+               sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption
+               // Analysis output
+               nCandidates++;
+       }
+       AliTRDstackLayer *fakeLayer = new AliTRDstackLayer(layers[0].GetZ0(), layers[0].GetDZ0(), layers[0].GetStackNr());
+       fakeLayer->SetX((TMath::Abs(layers[nTimeBins-1].GetX() + layers[0].GetX()))/2);
+       fakeLayer->SetSector(layers[0].GetSector());
+       AliTRDcluster **fakeClusters = 0x0;
+       UInt_t *fakeIndices = 0x0;
+       if(nCandidates){
+               fakeClusters = new AliTRDcluster*[nCandidates];
+               fakeIndices = new UInt_t[nCandidates];
+               UInt_t fakeIndex = 0;
+               for(Int_t ican = 0; ican < nCandidates; ican++){
+                       fakeClusters[ican] = new AliTRDcluster();
+                       fakeClusters[ican]->SetX(fakeLayer->GetX());
+                       fakeClusters[ican]->SetY(cogyvals[cand[ican]]);
+                       fakeClusters[ican]->SetZ(cogzvals[cand[ican]]);
+                       fakeClusters[ican]->SetSigmaZ2(sigcands[ican]);
+                       fakeIndices[ican] = fakeIndex++;// fantasy number
+               }
+       }
+       fakeLayer->SetRecoParam(fRecoParam);
+       fakeLayer->SetClustersArray(fakeClusters, nCandidates);
+       fakeLayer->SetIndexArray(fakeIndices);
+       fakeLayer->SetNRows(nRows);
+       fakeLayer->BuildIndices();
+       //fakeLayer->PrintClusters();
+       
+#ifdef DEBUG
+       if(AliTRDReconstructor::StreamLevel() >= 3){
+               TMatrixT<double> 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 = *fDebugStreamer;
+               cstreamer << "MakeSeedingLayer"
+                       << "Iteration="  << fSieveSeeding
+                       << "plane="      << plane
+                       << "ymin="       << ymin
+                       << "ymax="       << ymax
+                       << "zmin="       << zmin
+                       << "zmax="       << zmax
+                       << "L.="         << fakeLayer
+                       << "Histogram.=" << &hist
+                       << "\n";
+       }
+#endif
+       return fakeLayer;
+}
+
+//____________________________________________________________________
+void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig, Int_t planes[4])
+{
+  //
+  // Map seeding configurations to detector planes.
+  //
+  // Parameters :
+  //   iconfig : configuration index
+  //   planes  : member planes of this configuration. On input empty.
+  //
+  // Output :
+  //   planes : contains the planes which are defining the configuration
+  // 
+  // Detailed description
+  //
+  // Here is the list of seeding planes configurations together with
+  // their topological classification:
+  //
+  //  0 - 5432 TQ 0
+  //  1 - 4321 TQ 0
+  //  2 - 3210 TQ 0
+  //  3 - 5321 TQ 1
+  //  4 - 4210 TQ 1
+  //  5 - 5431 TQ 1
+  //  6 - 4320 TQ 1
+  //  7 - 5430 TQ 2
+  //  8 - 5210 TQ 2
+  //  9 - 5421 TQ 3
+  // 10 - 4310 TQ 3
+  // 11 - 5410 TQ 4
+  // 12 - 5420 TQ 5
+  // 13 - 5320 TQ 5
+  // 14 - 5310 TQ 5
+  //
+  // The topologic quality is modeled as follows:
+  // 1. The general model is define by the equation:
+  //  p(conf) = exp(-conf/2)
+  // 2. According to the topologic classification, configurations from the same
+  //    class are assigned the agerage value over the model values.
+  // 3. Quality values are normalized.
+  // 
+  // The topologic quality distribution as function of configuration is given below:
+  //Begin_Html
+  // <img src="gif/topologicQA.gif">
+  //End_Html
+  //
+
+       switch(iconfig){
+       case 0: // 5432 TQ 0
+               planes[0] = 2;
+               planes[1] = 3;
+               planes[2] = 4;
+               planes[3] = 5;
+               break;
+       case 1: // 4321 TQ 0
+               planes[0] = 1;
+               planes[1] = 2;
+               planes[2] = 3;
+               planes[3] = 4;
+               break;
+       case 2: // 3210 TQ 0
+               planes[0] = 0;
+               planes[1] = 1;
+               planes[2] = 2;
+               planes[3] = 3;
+               break;
+       case 3: // 5321 TQ 1
+               planes[0] = 1;
+               planes[1] = 2;
+               planes[2] = 3;
+               planes[3] = 5;
+               break;
+       case 4: // 4210 TQ 1
+               planes[0] = 0;
+               planes[1] = 1;
+               planes[2] = 2;
+               planes[3] = 4;
+               break;
+       case 5: // 5431 TQ 1
+               planes[0] = 1;
+               planes[1] = 3;
+               planes[2] = 4;
+               planes[3] = 5;
+               break;
+       case 6: // 4320 TQ 1
+               planes[0] = 0;
+               planes[1] = 2;
+               planes[2] = 3;
+               planes[3] = 4;
+               break;
+       case 7: // 5430 TQ 2
+               planes[0] = 0;
+               planes[1] = 3;
+               planes[2] = 4;
+               planes[3] = 5;
+               break;
+       case 8: // 5210 TQ 2
+               planes[0] = 0;
+               planes[1] = 1;
+               planes[2] = 2;
+               planes[3] = 5;
+               break;
+       case 9: // 5421 TQ 3
+               planes[0] = 1;
+               planes[1] = 2;
+               planes[2] = 4;
+               planes[3] = 5;
+               break;
+       case 10: // 4310 TQ 3
+               planes[0] = 0;
+               planes[1] = 1;
+               planes[2] = 3;
+               planes[3] = 4;
+               break;
+       case 11: // 5410 TQ 4
+               planes[0] = 0;
+               planes[1] = 1;
+               planes[2] = 4;
+               planes[3] = 5;
+               break;
+       case 12: // 5420 TQ 5
+               planes[0] = 0;
+               planes[1] = 2;
+               planes[2] = 4;
+               planes[3] = 5;
+               break;
+       case 13: // 5320 TQ 5
+               planes[0] = 0;
+               planes[1] = 2;
+               planes[2] = 3;
+               planes[3] = 5;
+               break;
+       case 14: // 5310 TQ 5
+               planes[0] = 0;
+               planes[1] = 1;
+               planes[2] = 3;
+               planes[3] = 5;
+               break;
+       }
+}
+
+//____________________________________________________________________
+void AliTRDtrackerV1::GetExtrapolationConfig(Int_t iconfig, Int_t planes[2])
+{
+  //
+  // Returns the extrapolation planes for a seeding configuration.
+  //
+  // Parameters :
+  //   iconfig : configuration index
+  //   planes  : planes which are not in this configuration. On input empty.
+  //
+  // Output :
+  //   planes : contains the planes which are not in the configuration
+  // 
+  // Detailed description
+  //
+
+       switch(iconfig){
+       case 0: // 5432 TQ 0
+               planes[0] = 1;
+               planes[1] = 0;
+               break;
+       case 1: // 4321 TQ 0
+               planes[0] = 5;
+               planes[1] = 0;
+               break;
+       case 2: // 3210 TQ 0
+               planes[0] = 4;
+               planes[1] = 5;
+               break;
+       case 3: // 5321 TQ 1
+               planes[0] = 4;
+               planes[1] = 0;
+               break;
+       case 4: // 4210 TQ 1
+               planes[0] = 5;
+               planes[1] = 3;
+               break;
+       case 5: // 5431 TQ 1
+               planes[0] = 2;
+               planes[1] = 0;
+               break;
+       case 6: // 4320 TQ 1
+               planes[0] = 5;
+               planes[1] = 1;
+               break;
+       case 7: // 5430 TQ 2
+               planes[0] = 2;
+               planes[1] = 1;
+               break;
+       case 8: // 5210 TQ 2
+               planes[0] = 4;
+               planes[1] = 3;
+               break;
+       case 9: // 5421 TQ 3
+               planes[0] = 3;
+               planes[1] = 0;
+               break;
+       case 10: // 4310 TQ 3
+               planes[0] = 5;
+               planes[1] = 2;
+               break;
+       case 11: // 5410 TQ 4
+               planes[0] = 3;
+               planes[1] = 2;
+               break;
+       case 12: // 5420 TQ 5
+               planes[0] = 3;
+               planes[1] = 1;
+               break;
+       case 13: // 5320 TQ 5
+               planes[0] = 4;
+               planes[1] = 1;
+               break;
+       case 14: // 5310 TQ 5
+               planes[0] = 4;
+               planes[1] = 2;
+               break;
+       }
+}
diff --git a/TRD/AliTRDtrackerV1.h b/TRD/AliTRDtrackerV1.h
new file mode 100644 (file)
index 0000000..6852d74
--- /dev/null
@@ -0,0 +1,93 @@
+#ifndef ALITRDTRACKERV1_H
+#define ALITRDTRACKERV1_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */ 
+
+/* $Id$ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+//  The TRD tracker                                                       //
+//                                                                        //
+//  Authors:                                                              //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                     //
+//    Markus Fasel <M.Fasel@gsi.de>                                       //
+//                                                                        //
+/////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALITRDTRACKER_H
+#include "AliTRDtracker.h"
+#endif
+
+#define DEBUG
+
+/**************************************************************************
+ * Class Status see source file                                           *
+ **************************************************************************/
+class TFile;
+class TTreeSRedirector;
+class TClonesArray;
+
+class AliRieman;
+class AliESDEvent;
+
+class AliTRDseedV1;
+class AliTRDstackLayer;
+class AliTRDtrackerFitter;
+class AliTRDrecoParam;
+
+class AliTRDtrackerV1 : public AliTRDtracker
+{
+
+ public:
+       enum{
+               kNTimeBins = 35,
+               kNPlanes = 6,
+               kNSeedPlanes = 4,
+               kMaxTracksStack = 100,
+               kNConfigs = 15
+       };
+       AliTRDtrackerV1(AliTRDrecoParam *p = 0x0);
+       AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p);
+       ~AliTRDtrackerV1();
+  
+       Int_t          Clusters2Tracks(AliESDEvent *esd);
+       void           GetSeedingConfig(Int_t iconfig, Int_t planes[4]);
+       void           GetExtrapolationConfig(Int_t iconfig, Int_t planes[2]);
+       void           SetRecoParam(AliTRDrecoParam *p){fRecoParam = p;}
+       
+ protected:
+
+       Double_t       BuildSeedingConfigs(AliTRDstackLayer *layer, Int_t *configs);
+       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);
+       Int_t          GetSeedingLayers(AliTRDstackLayer *layers, Double_t *params);
+       void           GetMeanCLStack(AliTRDstackLayer *layers, Int_t *planes, Double_t *params);
+       AliTRDcluster *FindSeedingCluster(AliTRDstackLayer *layers, AliTRDseedV1/*AliRieman*/ *sfit);
+       void           ImproveSeedQuality(AliTRDstackLayer *layer, AliTRDseedV1 *cseed);
+       Int_t          MakeSeeds(AliTRDstackLayer *layers, AliTRDseedV1 *sseed, Int_t *ipar);
+       AliTRDstackLayer *MakeSeedingLayer(AliTRDstackLayer *layers, Int_t Plane);
+       AliTRDtrack*   RegisterSeed(AliTRDseedV1 *seeds, Double_t *params);
+
+ private:
+
+       AliTRDtrackerV1(const AliTRDtrackerV1 &tracker);
+       AliTRDtrackerV1 &operator=(const AliTRDtrackerV1 &tracker);
+
+ private:
+
+       static Double_t      fTopologicQA[kNConfigs];         //  Topologic quality
+       Double_t             fTrackQuality[kMaxTracksStack];  //  Track quality 
+       Int_t                fSeedLayer[kMaxTracksStack];     //  Seed layer
+       Int_t                fSieveSeeding;                   //! Seeding iterator
+       AliTRDrecoParam     *fRecoParam;                      //  Reconstruction parameters
+       AliTRDtrackerFitter *fFitter;                         //! Fitter class of the tracker
+       TTreeSRedirector    *fDebugStreamer;                  //! Debug stream of the tracker
+
+       ClassDef(AliTRDtrackerV1, 1)                          //  Stand alone tracker development class
+
+};
+#endif
index 7051cbe..c61c39a 100644 (file)
@@ -20,6 +20,7 @@
 #pragma link C++ class  AliTRDtrack+;
 #pragma link C++ class  AliTRDtracklet+;
 #pragma link C++ class  AliTRDtracker+;
+#pragma link C++ class  AliTRDpropagationLayer+;
 #pragma link C++ class  AliTRDseed+;
 
 #pragma link C++ class  AliTRDpidESD+;
 
 #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+;
+
+
 #endif
index 75ea308..10e7a7b 100644 (file)
@@ -10,9 +10,15 @@ SRCS= AliTRDcluster.cxx \
       AliTRDpidESD.cxx \
       AliTRDRecParam.cxx \
       AliTRDReconstructor.cxx \
-      AliTRDtrackingAnalysis.cxx
+      AliTRDtrackingAnalysis.cxx \
+      AliTRDrecoParam.cxx \
+      AliTRDseedV1.cxx \
+      AliTRDtrackerV1.cxx \
+      AliTRDpropagationLayer.cxx \
+      AliTRDstackLayer.cxx \
+      AliTRDtrackerFitter.cxx
 
-HDRS= $(SRCS:.cxx=.h)                
+HDRS= $(SRCS:.cxx=.h)
 
 DHDR= TRDrecLinkDef.h