Update of tracking code (alignment/calibration awareness)
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Feb 2008 12:54:44 +0000 (12:54 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Feb 2008 12:54:44 +0000 (12:54 +0000)
18 files changed:
TRD/AliTRDReconstructor.cxx
TRD/AliTRDReconstructor.h
TRD/AliTRDchamberTimeBin.cxx [new file with mode: 0644]
TRD/AliTRDchamberTimeBin.h [new file with mode: 0644]
TRD/AliTRDcluster.cxx
TRD/AliTRDrecoParam.cxx
TRD/AliTRDseedV1.cxx
TRD/AliTRDseedV1.h
TRD/AliTRDtrackV1.cxx
TRD/AliTRDtrackV1.h
TRD/AliTRDtrackerDebug.cxx [new file with mode: 0644]
TRD/AliTRDtrackerDebug.h [new file with mode: 0644]
TRD/AliTRDtrackerV1.cxx
TRD/AliTRDtrackerV1.h
TRD/AliTRDtrackingChamber.cxx [new file with mode: 0644]
TRD/AliTRDtrackingChamber.h [new file with mode: 0644]
TRD/AliTRDtrackingSector.cxx [new file with mode: 0644]
TRD/AliTRDtrackingSector.h [new file with mode: 0644]

index 030725a..1b2d6a6 100644 (file)
@@ -43,6 +43,15 @@ ClassImp(AliTRDReconstructor)
 
 Bool_t AliTRDReconstructor::fgkSeedingOn  = kFALSE;
 Int_t  AliTRDReconstructor::fgStreamLevel = 0;      // Stream (debug) level
+AliTRDrecoParam* AliTRDReconstructor::fgRecoParam = 0x0; 
+
+
+//_____________________________________________________________________________
+AliTRDReconstructor::~AliTRDReconstructor()                  { 
+       if(fgRecoParam) delete fgRecoParam;
+}
+
+
 
 //_____________________________________________________________________________
 void AliTRDReconstructor::ConvertDigits(AliRawReader *rawReader
@@ -111,7 +120,9 @@ AliTracker *AliTRDReconstructor::CreateTracker() const
 
   //return new AliTRDtracker(NULL);
 
-  return new AliTRDtrackerV1(NULL, AliTRDrecoParam::GetLowFluxParam());
+       // TODO move it to rec.C. check TPC
+       fgRecoParam = AliTRDrecoParam::GetLowFluxParam();
+  return new AliTRDtrackerV1();
 
 }
 
index 517d987..eb00bb2 100644 (file)
 #include "AliReconstructor.h"
 
 class AliRawReader;
-
-class AliTRDReconstructor: public AliReconstructor {
+class AliTRDrecoParam;
+class AliTRDReconstructor: public AliReconstructor 
+{
 
  public:
 
-  AliTRDReconstructor():AliReconstructor()                       { };
-  virtual ~AliTRDReconstructor()                                 { };
-
-  //virtual Bool_t      HasDigitConversion() const                 { return kTRUE; };
-  virtual Bool_t      HasDigitConversion() const                 { return kFALSE; };
-  virtual void        ConvertDigits(AliRawReader *rawReader, TTree *digitsTree) const;
+  AliTRDReconstructor():AliReconstructor()                              { };
+  AliTRDReconstructor(const AliTRDReconstructor &r):AliReconstructor(r) { };
+  virtual ~AliTRDReconstructor();
+  AliTRDReconstructor& operator = (const AliTRDReconstructor& /*r*/) 
+                                                                 { return *this;            };
 
-  virtual void        Reconstruct(AliRawReader *rawReader, TTree *clusterTree) const;
-  virtual void        Reconstruct(TTree *digitsTree, TTree *clusterTree) const;
+  virtual Bool_t           HasDigitConversion() const            { return kFALSE;           };
+  virtual void             ConvertDigits(AliRawReader *rawReader, TTree *digitsTree) const;
 
-  virtual AliTracker *CreateTracker() const;
+  virtual void             Reconstruct(AliRawReader *rawReader, TTree *clusterTree) const;
+  virtual void             Reconstruct(TTree *digitsTree, TTree *clusterTree) const;
+  static  AliTRDrecoParam *RecoParam()                           { return fgRecoParam;      }
+  virtual AliTracker      *CreateTracker() const;
 
-  virtual void        FillESD(AliRawReader */*rawReader*/, TTree *clusterTree, AliESDEvent *esd) const
-  {FillESD((TTree*)NULL,clusterTree,esd);}
-  virtual void        FillESD(TTree *digitsTree, TTree *clusterTree, AliESDEvent *esd) const;
+  virtual void             FillESD(AliRawReader */*rawReader*/, TTree *clusterTree, AliESDEvent *esd) const
+                                                                 { FillESD((TTree * )NULL
+                                                                 , clusterTree
+                                                                 , esd);                    }
+  virtual void             FillESD(TTree *digitsTree, TTree *clusterTree, AliESDEvent *esd) const;
 
-  static  void        SetSeedingOn(Bool_t seeding)               { fgkSeedingOn  = seeding; }  
-  static  void        SetStreamLevel(Int_t level)                { fgStreamLevel = level;   }
+  static  void             SetSeedingOn(Bool_t seeding)          { fgkSeedingOn  = seeding; }  
+  static  void             SetStreamLevel(Int_t level)           { fgStreamLevel = level;   }
 
-  static  Bool_t      SeedingOn()                                { return fgkSeedingOn;     }
-  static  Int_t       StreamLevel()                              { return fgStreamLevel;    }
+  static  Bool_t           SeedingOn()                           { return fgkSeedingOn;     }
+  static  Int_t            StreamLevel()                         { return fgStreamLevel;    }
 
  private:
 
-  static  Bool_t   fgkSeedingOn;  //  Set flag for seeding during reconstruction
-  static  Int_t    fgStreamLevel; //  Flag for streaming
+  static  Bool_t           fgkSeedingOn;  //  Set flag for seeding during reconstruction
+  static  Int_t            fgStreamLevel; //  Flag for streaming
+  static  AliTRDrecoParam *fgRecoParam;   //  Reconstruction parameters
 
-  ClassDef(AliTRDReconstructor,0) //  Class for the TRD reconstruction
+  ClassDef(AliTRDReconstructor,0)         //  Class for the TRD reconstruction
 
 };
 
diff --git a/TRD/AliTRDchamberTimeBin.cxx b/TRD/AliTRDchamberTimeBin.cxx
new file mode 100644 (file)
index 0000000..fb0e0eb
--- /dev/null
@@ -0,0 +1,661 @@
+/**************************************************************************
+ * 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: AliTRDchamberTimeBin.cxx 23313 2008-01-11 14:56:43Z cblume $ */
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  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 "AliTRDchamberTimeBin.h"
+#include "AliTRDrecoParam.h"
+#include "AliTRDReconstructor.h"
+#include "AliTRDtrackerV1.h"
+
+
+ClassImp(AliTRDchamberTimeBin)
+
+//_____________________________________________________________________________
+AliTRDchamberTimeBin::AliTRDchamberTimeBin(Int_t plane, Int_t stack, Int_t sector, Double_t z0, Double_t zLength)
+  :TObject()
+  ,fOwner(kFALSE)
+  ,fPlane(plane)
+  ,fStack(stack)
+  ,fSector(sector)
+  ,fNRows(kMaxRows)
+  ,fN(0)
+  ,fX(0.)
+  ,fZ0(z0)
+  ,fZLength(zLength)
+{
+  //
+  // Default constructor (Only provided to use AliTRDchamberTimeBin with arrays)
+  //
+
+       for(int i=0; i<kMaxRows; i++) fPositions[i] = 0xff;
+       for(int ic=0; ic<kMaxClustersLayer; ic++){
+               fClusters[ic] = 0x0;
+               fIndex[ic]    = 0xffff;
+       }
+}
+
+// //_____________________________________________________________________________
+// AliTRDchamberTimeBin::AliTRDchamberTimeBin(const AliTRDpropagationLayer &layer, Double_t
+// z0, Double_t zLength, UChar_t stackNr):
+//     TObject()
+//     ,fOwner(kFALSE)
+//   ,fPlane(0xff)
+//   ,fStack(0xff)
+//   ,fSector(0xff)
+//   ,fNRows(kMaxRows)
+//   ,fN(0)
+//   ,fX(0.)
+//     ,fZ0(z0)
+//     ,fZLength(zLength)
+// {
+// // Standard constructor.
+// // Initialize also the underlying AliTRDpropagationLayer using the copy constructor.
+// 
+//     SetT0(layer.IsT0());
+//     for(int i=0; i<kMaxRows; i++) fPositions[i] = 0xff;
+//     for(int ic=0; ic<kMaxClustersLayer; ic++){
+//             fClusters[ic] = 0x0;
+//             fIndex[ic]    = 0xffff;
+//     }
+// }
+// 
+// //_____________________________________________________________________________
+// AliTRDchamberTimeBin::AliTRDchamberTimeBin(const AliTRDpropagationLayer &layer):
+//     TObject()
+//     ,fOwner(kFALSE)
+//   ,fPlane(0xff)
+//   ,fStack(0xff)
+//   ,fSector(0xff)
+//   ,fNRows(kMaxRows)
+//   ,fN(0)
+//   ,fX(0.)
+//     ,fZ0(0)
+//     ,fZLength(0)
+// {
+// // Standard constructor using only AliTRDpropagationLayer.
+//     
+//     SetT0(layer.IsT0());
+//     for(int i=0; i<kMaxRows; i++) fPositions[i] = 0xff;
+//     for(int ic=0; ic<kMaxClustersLayer; ic++){
+//             fClusters[ic] = 0x0;
+//             fIndex[ic]    = 0xffff;
+//     }
+// }
+// //_____________________________________________________________________________
+// AliTRDchamberTimeBin &AliTRDchamberTimeBin::operator=(const AliTRDpropagationLayer &layer)
+// {
+// // Assignment operator from an AliTRDpropagationLayer
+// 
+//     if (this != &layer) layer.Copy(*this);
+//     return *this;
+// }
+// 
+
+//_____________________________________________________________________________
+AliTRDchamberTimeBin::AliTRDchamberTimeBin(const AliTRDchamberTimeBin &layer):
+       TObject()
+       ,fOwner(layer.fOwner)
+  ,fPlane(layer.fPlane)
+  ,fStack(layer.fStack)
+  ,fSector(layer.fSector)
+       ,fNRows(layer.fNRows)
+  ,fN(layer.fN)
+  ,fX(layer.fX)
+       ,fZ0(layer.fZ0)
+       ,fZLength(layer.fZLength)
+{
+// Copy Constructor (performs a deep copy)
+       
+       SetT0(layer.IsT0());
+       for(int i=0; i<kMaxRows; i++) fPositions[i] = layer.fPositions[i];
+       memcpy(&fClusters[0], &layer.fClusters[0], kMaxClustersLayer*sizeof(UChar_t));
+       memcpy(&fIndex[0], &layer.fIndex[0], kMaxClustersLayer*sizeof(UInt_t));
+
+
+//     BuildIndices();
+}
+
+//_____________________________________________________________________________
+AliTRDchamberTimeBin &AliTRDchamberTimeBin::operator=(const AliTRDchamberTimeBin &layer)
+{
+// Assignment operator
+
+       if (this != &layer) layer.Copy(*this);
+  return *this;
+}
+
+//_____________________________________________________________________________
+void AliTRDchamberTimeBin::Copy(TObject &o) const
+{
+// Copy method. Performs a deep copy of all data from this object to object o.
+       
+       AliTRDchamberTimeBin &layer = (AliTRDchamberTimeBin &)o;
+       layer.fOwner       = kFALSE;
+       layer.fPlane       = fPlane;
+       layer.fStack       = fStack;
+       layer.fSector      = fSector;
+       layer.fNRows       = fNRows;
+       layer.fN           = fN;
+       layer.fX           = fX;
+       layer.fZ0          = fZ0;
+       layer.fZLength     = fZLength;
+       layer.SetT0(IsT0());
+       
+       for(int i = 0; i < kMaxRows; i++) layer.fPositions[i] = 0;
+       
+       for(int i=0; i<kMaxRows; i++) layer.fPositions[i] = fPositions[i];
+       memcpy(&layer.fClusters[0], &fClusters[0], kMaxClustersLayer*sizeof(UChar_t));
+       memcpy(&layer.fIndex[0], &fIndex[0], kMaxClustersLayer*sizeof(UInt_t));
+       
+       TObject::Copy(layer); // copies everything into layer
+       
+//     layer.BuildIndices();
+}
+
+//_____________________________________________________________________________
+AliTRDchamberTimeBin::~AliTRDchamberTimeBin()
+{
+// Destructor
+       if(fOwner) for(int ic=0; ic<fN; ic++) delete fClusters[ic];
+}
+
+//_____________________________________________________________________________
+void AliTRDchamberTimeBin::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 AliTRDchamberTimeBin::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) kMaxClustersLayer) {
+    //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++;
+
+}  
+
+
+//_____________________________________________________________________________
+void AliTRDchamberTimeBin::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++){
+               if(fClusters[i]->IsUsed()){
+                       fClusters[i] = 0x0;
+                       fIndex[i] = 0xffff;
+               } else nClStack++;
+       }
+       if(nClStack > kMaxClustersLayer) AliWarning(Form("Number of clusters in stack %d exceed buffer size %d. Truncating.", nClStack, kMaxClustersLayer));
+               
+       // Nothing in this time bin. Reset indexes 
+       if(!nClStack){
+               fN = 0;
+               memset(&fPositions[0], 0xff, sizeof(UChar_t) * kMaxRows);
+               memset(&fClusters[0], 0x0, sizeof(AliTRDcluster*) * kMaxClustersLayer);
+               memset(&fIndex[0], 0xffff, sizeof(UInt_t) * kMaxClustersLayer);
+               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]         = 0xffff;
+               nClStack++;
+       }
+       
+       // do clusters arrangement
+       fX = 0.;
+       fN =  nClStack;
+       nClStack = 0;
+       // Reset Positions array
+       memset(fPositions, 0, sizeof(UChar_t)*kMaxRows);
+       for(Int_t i = 0; i < fN; i++){
+               // boundary check
+               AliTRDcluster *cl = helpCL[i];
+               UChar_t rowIndex = cl->GetPadRow();
+               // Insert Leaf
+               Int_t pos = FindYPosition(cl->GetY(), rowIndex, i);
+               if(pos == -1){          // zbin is empty;
+                       Int_t upper = (rowIndex == fNRows - 1) ? nClStack : fPositions[rowIndex + 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 = rowIndex + 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 = rowIndex + 1; j < fNRows; j++) fPositions[j]++; 
+                       nClStack++;
+               }
+
+               // calculate mean x
+               fX += cl->GetX(); 
+               
+               // Debug Streaming
+               if(AliTRDtrackerV1::DebugStreamer() && AliTRDReconstructor::StreamLevel() >= 3){
+                       TTreeSRedirector &cstream = *AliTRDtrackerV1::DebugStreamer();
+                       cstream << "BuildIndices"
+                       << "Plane="    << fPlane
+                       << "Stack="    << fStack
+                       << "Sector="   << fSector
+                       << "Iter="     << iter
+                       << "C.="       << cl
+                       << "rowIndex=" << rowIndex
+                       << "\n";
+               }
+       }
+
+//     AliInfo("Positions");
+//     for(int ir=0; ir<fNRows; ir++) printf("pos[%d] %d\n", ir, fPositions[ir]);
+
+       fX /= fN;
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDchamberTimeBin::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 AliTRDchamberTimeBin::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 AliTRDchamberTimeBin::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 AliTRDchamberTimeBin::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 = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - z)/fZLength * fNRows);
+       if(z < fZ0) myZbin = fNRows - 1;
+       if(z > fZ0 + fZLength) myZbin = 0;
+       //printf("\n%f < %f < %f [%d]\n", fZ0, z, fZ0 + fZLength, myZbin);
+       //for(int ic=0; ic<fN; ic++) printf("%d z = %f row %d\n", ic, fClusters[ic]->GetZ(), fClusters[ic]->GetPadRow());
+
+       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){ 
+                       //printf("y[%f] ycl[%f] roady[%f]\n", y, ycl, 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 AliTRDchamberTimeBin::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(!AliTRDReconstructor::RecoParam()){
+               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] = AliTRDReconstructor::RecoParam()->GetMaxPhi()   * (cl->GetX() - GetX()) + 1.0;                        // deviation: y-Direction
+               cond[3] = AliTRDReconstructor::RecoParam()->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] = AliTRDReconstructor::RecoParam()->GetRoad0y() + phi;
+               cond[3] = AliTRDReconstructor::RecoParam()->GetRoad0z();
+       }
+}
+
+//_____________________________________________________________________________
+void AliTRDchamberTimeBin::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 - 1.E-3;
+
+       UChar_t zhi = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - zvals[0])/fZLength * fNRows);
+       UChar_t zlo = fNRows - 1 - (UChar_t)(TMath::Abs(fZ0 - zvals[1])/fZLength * fNRows);
+       
+       
+//     AliInfo(Form("yc[%f] zc[%f] dy[%f] dz[%f]", cond[0], cond[1], cond[2], cond[3]));
+//     PrintClusters();
+//     AliInfo(Form("zlo[%f] zhi[%f]", zvals[0], zvals[1]));
+//     AliInfo(Form("zlo[%d] zhi[%d]", zlo, zhi));
+       
+       //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;
+               //AliInfo(Form("z[%d] y [%d %d]", z, fPositions[z], upper));
+               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))){/*printf("exit z\n"); TODO*/ continue;}
+                       index[ncl] = y;
+                       ncl++;
+               }
+       }
+       if(ncl>fN) AliError(Form("Clusters found %d > %d (clusters in layer)", ncl, fN));
+}
+
+//_____________________________________________________________________________
+AliTRDcluster *AliTRDchamberTimeBin::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  = AliTRDReconstructor::RecoParam()->GetRoad2y();
+       Double_t maxroadz = AliTRDReconstructor::RecoParam()->GetRoad2z();
+       
+       Int_t index = SearchNearestCluster(cond[0],cond[1],maxroad,maxroadz);
+       AliTRDcluster *returnCluster = 0x0;
+       if(index != -1) returnCluster = (AliTRDcluster *) fClusters[index];
+       return returnCluster;
+}
+
+//_____________________________________________________________________________
+void AliTRDchamberTimeBin::PrintClusters() const
+{
+// Prints the position of each cluster in the stacklayer on the stdout
+//
+       printf("\nnRows = %d\n", fNRows);
+       printf("Z0 = %f\n", fZ0);
+       printf("Z1 = %f\n", fZ0+fZLength);
+       printf("clusters in AliTRDchamberTimeBin %d\n", fN);
+       for(Int_t i = 0; i < fN; i++){
+               printf("AliTRDchamberTimeBin: index=%i, Cluster: X = %3.3f [%d] Y = %3.3f [%d] Z = %3.3f [%d]\n", i,  fClusters[i]->GetX(), fClusters[i]->GetLocalTimeBin(), fClusters[i]->GetY(), fClusters[i]->GetPadCol(), fClusters[i]->GetZ(), fClusters[i]->GetPadRow());
+               if(fClusters[i]->IsUsed()) printf("cluster allready used. rejected in search algorithm\n");
+       }
+}
diff --git a/TRD/AliTRDchamberTimeBin.h b/TRD/AliTRDchamberTimeBin.h
new file mode 100644 (file)
index 0000000..c1f3a30
--- /dev/null
@@ -0,0 +1,106 @@
+#ifndef ALITRDCHAMBERTIMEBIN_H
+#define ALITRDCHAMBERTIMEBIN_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliTRDchamberTimeBin.h 22646 2007-11-29 18:13:40Z cblume $ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+//  A TRD layer in a single stack                                         //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ROOT_TObject
+#include "TObject.h"
+#endif
+
+#ifndef ALITRDCLUSTER_H
+#include "AliTRDcluster.h"
+#endif
+
+class AliTRDchamberTimeBin : public TObject
+{
+public:
+       enum{
+               kMaxClustersLayer = 150
+               ,kMaxRows = 16
+       };
+
+       AliTRDchamberTimeBin(Int_t plane=-1, Int_t stack=-1, Int_t sector=-1, Double_t z0=-1., Double_t zLength=-1.);
+       //AliTRDchamberTimeBin(const AliTRDpropagationLayer &layer, Double_t z0, Double_t zLength, UChar_t stackNr);
+       //AliTRDchamberTimeBin(const AliTRDpropagationLayer &layer);
+       AliTRDchamberTimeBin(const AliTRDchamberTimeBin &layer);
+       ~AliTRDchamberTimeBin();
+       //AliTRDchamberTimeBin   &operator=(const AliTRDpropagationLayer &myLayer);
+  operator Int_t() const                                        { return fN;                    }
+       AliTRDchamberTimeBin   &operator=(const AliTRDchamberTimeBin &myLayer);
+       AliTRDcluster      *operator[](const Int_t i) const {
+               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.);
+  inline void    Clear(const Option_t *opt = 0x0);
+  AliTRDcluster* GetCluster(Int_t index) const {return index < fN && index >= 0 ? 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       GetX()                            const {
+       return fX;      }
+       Double_t       GetZ0()                           const { return fZ0;     }
+       Double_t       GetDZ0()                          const { return fZLength;}
+       Int_t          GetNClusters()                    const { return fN; }
+       Int_t          GetPlane()                        const { return fPlane;  }
+       Int_t          GetStack()                        const { return fStack;  }
+       Int_t          GetSector()                       const { return fSector; }
+       void           InsertCluster(AliTRDcluster *c, UInt_t index);
+
+       Bool_t         IsT0() const {return TestBit(1);}
+       
+       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           SetPlane(Int_t plane){ fPlane = plane; }
+       void           SetStack(Int_t stack){ fStack = stack; }
+       void           SetSector(Int_t sector){ fSector = sector; }
+       void           SetOwner(Bool_t own = kTRUE) {fOwner = own;}
+       void           SetT0(Bool_t set=kTRUE) {SetBit(1, set);}
+       void           SetX(Double_t x) {fX = x;}
+private:
+       void           Copy(TObject &o) const;
+       Int_t          Find(Float_t y) 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
+       Char_t        fPlane;               // Plane number
+       Char_t        fStack;               //  stack number in supermodule
+       Char_t        fSector;              // Sector mumber
+       Char_t        fNRows;               //  number of pad rows in the chamber
+       UChar_t       fPositions[kMaxRows]; //  starting index of clusters in pad row 
+       Int_t         fN;                   // number of clusters
+       AliTRDcluster *fClusters[kMaxClustersLayer];            //Array of pointers to clusters
+       UInt_t        fIndex[kMaxClustersLayer];                //Array of cluster indexes
+       Double_t      fX;                   //  radial position of tb
+       
+       // obsolete !!
+       Double_t      fZ0;                  //  starting position of the layer in Z direction
+       Double_t      fZLength;             //  length of the layer in Z direction
+       
+       ClassDef(AliTRDchamberTimeBin, 2)   //  tracking propagation layer for one time bin in chamber
+
+};
+
+inline 
+void AliTRDchamberTimeBin::Clear(const Option_t *) 
+{ 
+       for (Int_t i = 0; i < fN; i++) fClusters[i] = NULL;
+       fN = 0; 
+}
+
+#endif // ALITRDCHAMBERTIMEBIN_H
+
index 6a4f1f2..0cc0a22 100644 (file)
@@ -94,6 +94,7 @@ AliTRDcluster::AliTRDcluster(const AliTRDcluster &c)
   // Copy constructor 
   //
 
+  SetBit(1, c.IsInChamber());
   SetLabel(c.GetLabel(0),0);
   SetLabel(c.GetLabel(1),1);
   SetLabel(c.GetLabel(2),2);
index e771511..4ce768d 100644 (file)
@@ -44,7 +44,7 @@ AliTRDrecoParam::AliTRDrecoParam()
   ,fkRoad2z(20.0)
   ,fkPlaneQualityThreshold(5.0)// 4.2? under Investigation
   ,fkFindable(.333)
-  ,fkChi2Z(14./*12.5*/)
+  ,fkChi2Z(30./*14.*//*12.5*/)
   ,fkChi2Y(.25)
   ,fkTrackLikelihood(-15.)
 {
index 020f4d6..ff89e04 100644 (file)
@@ -27,6 +27,8 @@
 
 #include "TMath.h"
 #include "TLinearFitter.h"
+#include "TClonesArray.h" // tmp
+#include <TTreeStream.h>
 
 #include "AliLog.h"
 #include "AliMathBase.h"
 #include "AliTRDcluster.h"
 #include "AliTRDtrack.h"
 #include "AliTRDcalibDB.h"
-#include "AliTRDstackLayer.h"
+#include "AliTRDchamberTimeBin.h"
+#include "AliTRDtrackingChamber.h"
+#include "AliTRDtrackerV1.h"
+#include "AliTRDReconstructor.h"
 #include "AliTRDrecoParam.h"
 #include "AliTRDgeometry.h"
 #include "Cal/AliTRDCalPID.h"
 
-#define SEED_DEBUG
-
 ClassImp(AliTRDseedV1)
 
 //____________________________________________________________________
-AliTRDseedV1::AliTRDseedV1(Int_t layer, AliTRDrecoParam *p) 
+AliTRDseedV1::AliTRDseedV1(Int_t plane) 
   :AliTRDseed()
-  ,fPlane(layer)
+  ,fPlane(plane)
   ,fOwner(kFALSE)
   ,fMom(0.)
   ,fSnp(0.)
   ,fTgl(0.)
   ,fdX(0.)
-  ,fRecoParam(p)
 {
   //
   // Constructor
@@ -71,7 +73,6 @@ AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
   ,fSnp(ref.fSnp)
   ,fTgl(ref.fTgl)
   ,fdX(ref.fdX)
-  ,fRecoParam(ref.fRecoParam)
 {
   //
   // Copy Constructor performing a deep copy
@@ -132,7 +133,6 @@ void AliTRDseedV1::Copy(TObject &ref) const
        target.fSnp           = fSnp;
        target.fTgl           = fTgl;
        target.fdX            = fdX;
-       target.fRecoParam     = fRecoParam;
        
        for(int islice=0; islice < knSlices; islice++) target.fdEdx[islice] = fdEdx[islice];
        for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) target.fProb[ispec] = fProb[ispec];
@@ -254,7 +254,7 @@ Double_t* AliTRDseedV1::GetProbability()
   }
 
   // Retrieve the CDB container class with the parametric detector response
-  const AliTRDCalPID *pd = calibration->GetPIDObject(fRecoParam->GetPIDMethod());
+  const AliTRDCalPID *pd = calibration->GetPIDObject(AliTRDReconstructor::RecoParam()->GetPIDMethod());
   if (!pd) {
     AliError("No access to AliTRDCalPID object");
     return 0x0;
@@ -265,7 +265,7 @@ Double_t* AliTRDseedV1::GetProbability()
   /// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
   
   //calculate dE/dx
-  CookdEdx(fRecoParam->GetNdEdxSlices());
+  CookdEdx(AliTRDReconstructor::RecoParam()->GetNdEdxSlices());
   
   // Sets the a priori probabilities
   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
@@ -283,9 +283,10 @@ Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
   //
 
        Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
-       return .5 * (18.0 - fN2)
+       return 
+                 .5 * TMath::Abs(18.0 - fN2)
                + 10.* TMath::Abs(fYfit[1] - fYref[1])
-               + 5.* TMath::Abs(fYfit[0] - fYref[0] + zcorr)
+               + 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
                + 2. * TMath::Abs(fMeanz - fZref[0]) / fPadLength;
 }
 
@@ -294,12 +295,14 @@ void AliTRDseedV1::GetCovAt(Double_t /*x*/, Double_t *cov) const
 {
 // Computes covariance in the y-z plane at radial point x
 
-       const Float_t k0= .2; // to be checked in FindClusters
-       Double_t sy20   = k0*TMath::Tan(fYfit[1]); sy20 *= sy20;
-       
-       Double_t sy2    = fSigmaY2*fSigmaY2 + sy20;
+       Int_t ic = 0; while (!fClusters[ic]) ic++; 
+  AliTRDcalibDB *fCalib = AliTRDcalibDB::Instance();
+       Double_t exB         = fCalib->GetOmegaTau(fCalib->GetVdriftAverage(fClusters[ic]->GetDetector()), -AliTracker::GetBz()*0.1);
+
+       Double_t sy2    = fSigmaY2*fSigmaY2 + .2*(fYfit[1]-exB)*(fYfit[1]-exB);
        Double_t sz2    = fPadLength/12.;
 
+
        //printf("Yfit[1] %f sy20 %f SigmaY2 %f\n", fYfit[1], sy20, fSigmaY2);
 
        cov[0] = sy2;
@@ -332,26 +335,51 @@ void AliTRDseedV1::SetOwner(Bool_t own)
 }
 
 //____________________________________________________________________
-Bool_t AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
-                                       , Float_t quality
-                                       , Bool_t kZcorr
-                                       , AliTRDcluster *c)
+Bool_t AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, 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)
   //
+       // debug level 7
+       //
        
-       if(!fRecoParam){
+       if(!AliTRDReconstructor::RecoParam()){
                AliError("Seed can not be used without a valid RecoParam.");
                return kFALSE;
        }
 
-       //AliInfo(Form("TimeBins = %d TimeBinsRange = %d", fgTimeBins, fTimeBinsRange));
+       AliTRDchamberTimeBin *layer = 0x0;
+       if(AliTRDReconstructor::StreamLevel()>=7 && c){
+               TClonesArray clusters("AliTRDcluster", 24);
+               clusters.SetOwner(kTRUE);
+               AliTRDcluster *cc = 0x0;
+               Int_t det=-1, ncl, ncls = 0;
+               for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
+                       if(!(layer = chamber->GetTB(iTime))) continue;
+                       if(!(ncl = Int_t(*layer))) continue;
+                       for(int ic=0; ic<ncl; ic++){ 
+                               cc = (*layer)[ic];
+                               det = cc->GetDetector();
+                               new(clusters[ncls++]) AliTRDcluster(*cc);
+                       }
+               }
+               AliInfo(Form("N clusters[%d] = %d", fPlane, ncls));
+               
+               Int_t ref = c ? 1 : 0;
+               TTreeSRedirector &cstreamer = *AliTRDtrackerV1::DebugStreamer();
+               cstreamer << "AttachClustersIter"
+                       << "det="        << det 
+                       << "ref="        << ref 
+                       << "clusters.="  << &clusters
+                       << "tracklet.="  << this
+                       << "cl.="        << c
+                       << "\n";        
+       }
 
        Float_t  tquality;
-       Double_t kroady = fRecoParam->GetRoad1y();
+       Double_t kroady = AliTRDReconstructor::RecoParam()->GetRoad1y();
        Double_t kroadz = fPadLength * .5 + 1.;
        
        // initialize configuration parameters
@@ -362,11 +390,13 @@ Bool_t    AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
        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 < fgTimeBins; iTime++) {
+                       if(!(layer = chamber->GetTB(iTime))) continue;
+                       if(!Int_t(*layer)) continue;
+                       
                        // define searching configuration
-                       Double_t dxlayer = layer[iTime].GetX() - fX0;
+                       Double_t dxlayer = layer->GetX() - fX0;
                        if(c){
                                zexp = c->GetZ();
                                //Try 2 pad-rows in second iteration
@@ -375,23 +405,23 @@ Bool_t    AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
                                        if (zexp > c->GetZ()) zexp = c->GetZ() + fPadLength*0.5;
                                        if (zexp < c->GetZ()) zexp = c->GetZ() - fPadLength*0.5;
                                }
-                       } else zexp = fZref[0];
+                       } else zexp = fZref[0] + (kZcorr ? fZref[1] * dxlayer : 0.);
                        yexp  = fYref[0] + fYref[1] * dxlayer - zcorr;
                        
                        // Get and register cluster
-                       Int_t    index = layer[iTime].SearchNearestCluster(yexp, zexp, kroady, kroadz);
+                       Int_t    index = layer->SearchNearestCluster(yexp, zexp, kroady, kroadz);
                        if (index < 0) continue;
-                       AliTRDcluster *cl = (AliTRDcluster*) layer[iTime].GetCluster(index);
+                       AliTRDcluster *cl = (*layer)[index];
                        
-                       Int_t globalIndex = layer[iTime].GetGlobalIndex(index);
-                       fIndexes[iTime]  = globalIndex;
+                       fIndexes[iTime]  = layer->GetGlobalIndex(index);
                        fClusters[iTime] = cl;
                        fY[iTime]        = cl->GetY();
                        fZ[iTime]        = cl->GetZ();
                        ncl++;
                }
+       if(AliTRDReconstructor::StreamLevel()>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fPlane, ncl));
                
-               if(ncl){        
+               if(ncl>1){      
                        // calculate length of the time bin (calibration aware)
                        Int_t irp = 0; Float_t x[2]; Int_t tb[2];
                        for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
@@ -405,7 +435,8 @@ Bool_t      AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
        
                        // update X0 from the clusters (calibration/alignment aware)
                        for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
-                               if(!layer[iTime].IsT0()) continue;
+                               if(!(layer = chamber->GetTB(iTime))) continue;
+                               if(!layer->IsT0()) continue;
                                if(fClusters[iTime]){ 
                                        fX0 = fClusters[iTime]->GetX();
                                        break;
@@ -429,6 +460,7 @@ Bool_t      AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
                        
                        AliTRDseed::Update();
                }
+       if(AliTRDReconstructor::StreamLevel()>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fPlane, fN2));
                
                if(IsOK()){
                        tquality = GetQuality(kZcorr);
@@ -445,7 +477,7 @@ Bool_t      AliTRDseedV1::AttachClustersIter(AliTRDstackLayer *layer
 }
 
 //____________________________________________________________________
-Bool_t AliTRDseedV1::AttachClusters(AliTRDstackLayer *layer
+Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber
                                        ,Bool_t kZcorr)
 {
   //
@@ -464,7 +496,7 @@ Bool_t      AliTRDseedV1::AttachClusters(AliTRDstackLayer *layer
   // 6. fit tracklet
   //   
 
-       if(!fRecoParam){
+       if(!AliTRDReconstructor::RecoParam()){
                AliError("Seed can not be used without a valid RecoParam.");
                return kFALSE;
        }
@@ -472,7 +504,7 @@ Bool_t      AliTRDseedV1::AttachClusters(AliTRDstackLayer *layer
        const Int_t kClusterCandidates = 2 * knTimebins;
        
        //define roads
-       Double_t kroady = fRecoParam->GetRoad1y();
+       Double_t kroady = AliTRDReconstructor::RecoParam()->GetRoad1y();
        Double_t kroadz = fPadLength * 1.5 + 1.;
        // correction to y for the tilting angle
        Float_t zcorr = kZcorr ? fTilt * (fZProb - fZref[0]) : 0.;
@@ -484,18 +516,22 @@ Bool_t    AliTRDseedV1::AttachClusters(AliTRDstackLayer *layer
        Int_t ncl, *index = 0x0, tboundary[knTimebins];
        
        // Do cluster projection
+       AliTRDchamberTimeBin *layer = 0x0;
        Int_t nYclusters = 0; Bool_t kEXIT = kFALSE;
        for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
-               fX[iTime] = layer[iTime].GetX() - fX0;
+               if(!(layer = chamber->GetTB(iTime))) continue;
+               if(!Int_t(*layer)) continue;
+               
+               fX[iTime] = layer->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);
+               layer->GetClusters(cond, index, ncl);
                for(Int_t ic = 0; ic<ncl; ic++){
-                       AliTRDcluster *c = layer[iTime].GetCluster(index[ic]);
+                       AliTRDcluster *c = layer->GetCluster(index[ic]);
                        clusters[nYclusters] = c;
                        yres[nYclusters++] = c->GetY() - yexp[iTime];
                        if(nYclusters >= kClusterCandidates) {
@@ -511,7 +547,7 @@ Bool_t      AliTRDseedV1::AttachClusters(AliTRDstackLayer *layer
        // 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
+       // purge cluster candidates
        Int_t nZclusters = 0;
        for(Int_t ic = 0; ic<nYclusters; ic++){
                if(yres[ic] - mean > 4. * sigma){
@@ -523,7 +559,7 @@ Bool_t      AliTRDseedV1::AttachClusters(AliTRDstackLayer *layer
        
        // Evaluate truncated mean on the z direction
        AliMathBase::EvaluateUni(nZclusters, zres, mean, sigma, Int_t(nZclusters*.8)-2);
-       //purge cluster candidates
+       // purge cluster candidates
        for(Int_t ic = 0; ic<nZclusters; ic++){
                if(zres[ic] - mean > 4. * sigma){
                        clusters[ic] = 0x0;
@@ -538,11 +574,9 @@ Bool_t     AliTRDseedV1::AttachClusters(AliTRDstackLayer *layer
        for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
                ncl = tboundary[iTime] - lastCluster;
                if(!ncl) continue;
-               AliTRDcluster *c = 0x0;
-               if(ncl == 1){
-                       c = clusters[lastCluster];
-               } else if(ncl > 1){
-                       Float_t dold = 9999.; Int_t iptr = lastCluster;
+               Int_t iptr = lastCluster;
+               if(ncl > 1){
+                       Float_t dold = 9999.;
                        for(int ic=lastCluster; ic<tboundary[iTime]; ic++){
                                if(!clusters[ic]) continue;
                                Float_t y = yexp[iTime] - clusters[ic]->GetY();
@@ -552,19 +586,17 @@ Bool_t    AliTRDseedV1::AttachClusters(AliTRDstackLayer *layer
                                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();
+               fIndexes[iTime]  = chamber->GetTB(iTime)->GetGlobalIndex(iptr);
+               fClusters[iTime] = clusters[iptr];
+               fY[iTime]        = clusters[iptr]->GetY();
+               fZ[iTime]        = clusters[iptr]->GetZ();
                lastCluster      = tboundary[iTime];
                fN2++;
        }
        
        // number of minimum numbers of clusters expected for the tracklet
-       Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fgTimeBins);
+       Int_t kClmin = Int_t(AliTRDReconstructor::RecoParam()->GetFindableClusters()*fgTimeBins);
   if (fN2 < kClmin){
                AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
     fN2 = 0;
@@ -623,7 +655,7 @@ Bool_t AliTRDseedV1::Fit()
        }
 
        // calculate pad row boundary crosses
-       Int_t kClmin = Int_t(fRecoParam->GetFindableClusters()*fgTimeBins);
+       Int_t kClmin = Int_t(AliTRDReconstructor::RecoParam()->GetFindableClusters()*fgTimeBins);
        Int_t nz = AliMathBase::Freq(fN, zint, zout, kFALSE);
   fZProb   = zout[0];
   if(nz <= 1) zout[3] = 0;
@@ -704,137 +736,6 @@ Bool_t AliTRDseedV1::Fit()
        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];
-      error *= terror ? cseed[iLayer].fSigmaY : .2;
-
-//                     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()
index c8055ab..bb0847c 100644 (file)
@@ -27,11 +27,9 @@ class TTreeSRedirector;
 
 class AliRieman;
 
-class AliTRDstackLayer;
+class AliTRDtrackingChamber;
 class AliTRDcluster;
-class AliTRDrecoParam;
 class AliTRDtrack;
-
 class AliTRDseedV1 : public AliTRDseed
 {
 
@@ -41,23 +39,22 @@ class AliTRDseedV1 : public AliTRDseed
          knSlices = 10
        };
 
-       AliTRDseedV1(Int_t layer = -1, AliTRDrecoParam *p=0x0);
+       AliTRDseedV1(Int_t plane = -1);
        ~AliTRDseedV1();
        AliTRDseedV1(const AliTRDseedV1 &ref);
        AliTRDseedV1& operator=(const AliTRDseedV1 &ref);
 
-       Bool_t  AttachClustersIter(AliTRDstackLayer *layer, Float_t quality, Bool_t kZcorr = kFALSE
+       Bool_t  AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t quality, Bool_t kZcorr = kFALSE
                                  , AliTRDcluster *c=0x0);
-       Bool_t  AttachClusters(AliTRDstackLayer *layer, Bool_t kZcorr = kFALSE);
+       Bool_t  AttachClusters(AliTRDtrackingChamber *chamber, Bool_t kZcorr = kFALSE);
        void    CookdEdx(Int_t nslices);
-       static  Float_t FitRiemanTilt(AliTRDseedV1 * cseed, Bool_t terror);
        Bool_t  Fit();
 
               void      Init(AliTRDtrack *track);
        inline void      Init(const AliRieman *fit);
        
-       inline Float_t   GetChi2Z(const Float_t z = 0.) const;
-       inline Float_t   GetChi2Y(const Float_t y = 0.) const;
+       inline Float_t   GetChi2Z(const Float_t z = 999.) const;
+       inline Float_t   GetChi2Y(const Float_t y = 999.) const;
               void      GetCovAt(Double_t x, Double_t *cov) const;
               Float_t*  GetdEdx() {return &fdEdx[0];}
               Float_t   GetdQdl(Int_t ic) const;
@@ -78,7 +75,6 @@ class AliTRDseedV1 : public AliTRDseed
               void      SetMomentum(Double_t mom) {fMom = mom;}
               void      SetOwner(Bool_t own = kTRUE);
               void      SetPlane(Int_t p)                      { fPlane     = p;   }
-              void      SetRecoParam(AliTRDrecoParam *p)       { fRecoParam = p;   }
               void      SetSnp(Double_t snp) {fSnp = snp;}
               void      SetTgl(Double_t tgl) {fTgl = tgl;}
 
@@ -96,7 +92,6 @@ class AliTRDseedV1 : public AliTRDseed
        Float_t          fdX;                     // length of time bin
        Float_t          fdEdx[knSlices];         //  dE/dx measurements for tracklet
        Double_t         fProb[AliPID::kSPECIES]; //  PID probabilities
-       AliTRDrecoParam *fRecoParam;              //! Local copy of the reco params 
 
        ClassDef(AliTRDseedV1, 1)                 //  New TRD seed 
 
@@ -105,7 +100,7 @@ class AliTRDseedV1 : public AliTRDseed
 //____________________________________________________________
 inline Float_t AliTRDseedV1::GetChi2Z(const Float_t z) const
 {
-       Float_t z1  = (z == 0.) ? fMeanz : z;
+       Float_t z1  = (z == 999.) ? fMeanz : z;
        Float_t chi = fZref[0] - z1;
        return chi*chi;
 }
@@ -113,7 +108,7 @@ inline Float_t AliTRDseedV1::GetChi2Z(const Float_t z) const
 //____________________________________________________________
 inline Float_t AliTRDseedV1::GetChi2Y(const Float_t y) const
 {
-       Float_t y1  = (y == 0.) ? fYfitR[0] : y;
+       Float_t y1  = (y == 999.) ? fYfitR[0] : y;
        Float_t chi = fYref[0] - y1;
        return chi*chi;
 }
index c9a4542..771e57a 100644 (file)
 
 /* $Id$ */
 
+#include "AliESDtrack.h"
+
 #include "AliTRDtrackV1.h"
 #include "AliTRDcluster.h"
 #include "AliTRDcalibDB.h"
 #include "AliTRDrecoParam.h"
 
-#include "AliESDtrack.h"
-
 ClassImp(AliTRDtrackV1)
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -38,7 +38,6 @@ ClassImp(AliTRDtrackV1)
 //_______________________________________________________________
 AliTRDtrackV1::AliTRDtrackV1() 
   :AliTRDtrack()
-  ,fRecoParam(0x0)
 {
   //
   // Default constructor
@@ -53,7 +52,6 @@ AliTRDtrackV1::AliTRDtrackV1()
 //_______________________________________________________________
 AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) 
   :AliTRDtrack(t)
-  ,fRecoParam(0x0)
 {
   //
   // Constructor from AliESDtrack
@@ -67,7 +65,6 @@ AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t)
 //_______________________________________________________________
 AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref) 
   :AliTRDtrack(ref)
-  ,fRecoParam(ref.fRecoParam)
 {
   //
   // Copy constructor
@@ -86,11 +83,9 @@ AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref)
 // }
        
 //_______________________________________________________________
-AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5]
-                           , const Double_t cov[15]
-                           , Double_t x, Double_t alpha)
+AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Double_t cov[15]
+             , Double_t x, Double_t alpha)
   :AliTRDtrack()
-  ,fRecoParam(0x0)
 {
   //
   // The stand alone tracking constructor
@@ -145,7 +140,7 @@ AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5]
                //printf("%d N clusters plane %d [%d %d].\n", iplane, nclPlane, fTracklet[iplane].GetN2(), trklts[iplane].GetN());
        }
        //printf("N clusters in AliTRDtrackV1 %d.\n", ncl);
-       SetNumberOfClusters(ncl);
+       SetNumberOfClusters(/*ncl*/);
 }
 
 //_______________________________________________________________
@@ -247,6 +242,19 @@ Bool_t AliTRDtrackV1::IsOwner() const
        return kTRUE;
 }
        
+//___________________________________________________________
+void AliTRDtrackV1::SetNumberOfClusters() 
+{
+// Calculate the number of clusters attached to this track
+       
+       Int_t ncls = 0;
+       for(int ip=0; ip<6; ip++){
+               if(fTrackletIndex[ip] >= 0) ncls += fTracklet[ip].GetN();
+       }
+       AliKalmanTrack::SetNumberOfClusters(ncls);      
+}
+
+       
 //_______________________________________________________________
 void AliTRDtrackV1::SetOwner(Bool_t own)
 {
@@ -268,7 +276,7 @@ void AliTRDtrackV1::SetTracklet(AliTRDseedV1 *trklt, Int_t plane, Int_t index)
   // Set the tracklets
   //
 
-       if(plane < 0 || plane >=6) return;
+       if(plane < 0 || plane >= AliESDtrack::kNPlane) return;
        fTracklet[plane]      = (*trklt);
        fTrackletIndex[plane] = index;
 }
@@ -296,7 +304,7 @@ Bool_t  AliTRDtrackV1::Update(AliTRDseedV1 *trklt, Double_t chisq)
 //   Int_t n      = GetNumberOfClusters();
 //   fIndex[n]    = index;
 //   fClusters[n] = c;
-  SetNumberOfClusters(GetNumberOfClusters()+trklt->GetN());
+  SetNumberOfClusters(/*GetNumberOfClusters()+trklt->GetN()*/);
   SetChi2(GetChi2() + chisq);
        
        // update tracklet
index 90b1736..521b93b 100644 (file)
@@ -31,7 +31,6 @@ class AliTRDtrackV1 : public AliTRDtrack
 
        Bool_t         CookPID();
        Float_t        GetMomentum(Int_t plane) const;
-       inline Int_t   GetNumberOfClusters() const;
        Double_t       GetPredictedChi2(const AliTRDseedV1 *tracklet) const;
         Double_t       GetPredictedChi2(const AliCluster* /*c*/) const                   { return 0.0; }
        const AliTRDseedV1* GetTracklet(Int_t plane) const {return plane >=0 && plane <6 ? &fTracklet[plane] : 0x0;}
@@ -39,40 +38,19 @@ class AliTRDtrackV1 : public AliTRDtrack
        
        Bool_t         IsOwner() const;
        
+       void           SetNumberOfClusters();
        void           SetOwner(Bool_t own = kTRUE);
        void           SetTracklet(AliTRDseedV1 *trklt, Int_t plane, Int_t index);
        Bool_t         Update(AliTRDseedV1 *tracklet, Double_t chi2);
-        Bool_t         Update(const AliTRDcluster*, Double_t, Int_t, Double_t) { return kFALSE; };
-        Bool_t         Update(const AliCluster *, Double_t, Int_t)             { return kFALSE; };
+        Bool_t         Update(const AliTRDcluster *c, Double_t chi2, Int_t index, Double_t h01) 
+                                                           { return AliTRDtrack::Update(c,chi2,index,h01); };
+        Bool_t         Update(const AliCluster *, Double_t, Int_t)                        { return kFALSE; };
        void           UpdateESDtrack(AliESDtrack *t);
 
- protected:
-       AliTRDrecoParam *fRecoParam;       // reconstruction parameters
-
        ClassDef(AliTRDtrackV1, 1)         // development TRD track
 
 };
 
-//___________________________________________________________
-inline Int_t AliTRDtrackV1::GetNumberOfClusters() const
-{
-/*     Int_t ntrklts = GetNumberOfTracklets();
-       printf("AliTRDtrackV1::GetNumberOfClusters() %d\n", ntrklts);
-       return ntrklts;*/
-       Int_t ncls = 0;
-       for(int ip=0; ip<6; ip++)
-               if(fTrackletIndex[ip] >= 0) ncls+=fTracklet[ip].GetN();
-               
-       return ncls;
-}
-
-// //___________________________________________________________
-// inline Int_t AliTRDtrackV1::GetNumberOfTracklets() const
-// {
-//     Int_t ntrklt = 0;
-//     for(int ip=0; ip<6; ip++) if(fTrackletIndex[ip] >= 0) ntrklt++;
-//     return ntrklt;
-// }
 
 #endif
 
diff --git a/TRD/AliTRDtrackerDebug.cxx b/TRD/AliTRDtrackerDebug.cxx
new file mode 100644 (file)
index 0000000..bc7a53f
--- /dev/null
@@ -0,0 +1,311 @@
+/**************************************************************************
+ * 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: AliTRDtrackerDebug.cxx 23810 2008-02-08 09:00:27Z hristov $ */
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Tracker debug streamer                                                   //
+//                                                                           //
+//  Authors:                                                                 //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                        //
+//    Markus Fasel <M.Fasel@gsi.de>                                          //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDtrackerDebug.h"
+
+#include "TFile.h"
+#include "TTree.h"
+#include "TTreeStream.h"
+#include "TLinearFitter.h"
+
+#include "AliLog.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDtrackV1.h"
+#include "AliTRDseedV1.h"
+#include "AliTRDseed.h"
+#include "AliTRDcluster.h"
+
+ClassImp(AliTRDtrackerDebug)
+
+//____________________________________________________
+AliTRDtrackerDebug::AliTRDtrackerDebug() : AliTRDtrackerV1()
+       ,fOutputStreamer(0x0)
+       ,fTree(0x0)
+       ,fTracklet(0x0)
+       ,fTrack(0x0)
+       ,fNClusters(0)
+       ,fAlpha(0.)
+{
+        //
+       // Default constructor
+       //
+       fOutputStreamer = new TTreeSRedirector("TRD.Debug.root");
+}
+
+//____________________________________________________
+AliTRDtrackerDebug::~AliTRDtrackerDebug()
+{
+       // destructor
+       
+       delete fOutputStreamer;
+}
+
+
+//____________________________________________________
+void AliTRDtrackerDebug::Draw(const Option_t *)
+{
+// steer draw function
+}
+
+
+//____________________________________________________
+Bool_t AliTRDtrackerDebug::Init()
+{
+// steer linking data for various debug streams        
+       
+       fTrack = new AliTRDtrackV1();
+       fTree->SetBranchAddress("ncl", &fNClusters);
+       fTree->SetBranchAddress("track.", &fTrack);
+       return kTRUE;
+}
+
+//____________________________________________________
+Bool_t AliTRDtrackerDebug::Open(const char *method)
+{
+       // Connect to the tracker debug file
+       
+       TDirectory *savedir = gDirectory; 
+       TFile::Open("TRD.TrackerDebugger.root");
+       fTree = (TTree*)gFile->Get(method);
+       if(!fTree){
+               AliInfo(Form("Can not find debug stream for the %s method.\n", method));
+               savedir->cd();
+               return kFALSE;
+       }
+       savedir->cd();
+       return kTRUE;
+}
+
+//____________________________________________________
+Int_t AliTRDtrackerDebug::Process()
+{
+// steer debug process threads
+       
+       for(int it = 0; it<fTree->GetEntries(); it++){
+               if(!fTree->GetEntry(it)) continue;
+               if(!fNClusters) continue;
+               fAlpha = fTrack->GetAlpha();
+               //printf("Processing track %d [%d] ...\n", it, fNClusters);
+               ResidualsTrackletsTrack();
+
+               const AliTRDseedV1 *tracklet = 0x0;
+               for(int ip = 5; ip>=0; ip--){
+                       if(!(tracklet = fTrack->GetTracklet(ip))) continue;
+                       if(!tracklet->GetN()) continue;
+                       
+                       ResidualsClustersTrack(tracklet);
+                       ResidualsClustersTracklet(tracklet);
+                       ResidualsClustersParametrisation(tracklet);
+               }
+       }
+       return kTRUE;
+}
+
+
+//____________________________________________________
+void AliTRDtrackerDebug::ResidualsClustersTrack(const AliTRDseedV1 *tracklet)
+{
+// Calculate averange distances from clusters to the TRD track 
+       
+       Double_t x[3]; 
+       AliTRDcluster *c = 0x0;
+       for(int ic=0; ic<35/*AliTRDseed:knTimebins*/; ic++){
+               if(!(c = tracklet->GetClusters(ic))) continue;
+               Double_t xc = c->GetX(), yc = c->GetY(), zc = c->GetZ();
+
+               // propagate track to cluster 
+               PropagateToX(*fTrack, xc, 2.); 
+               fTrack->GetXYZ(x);
+               
+               // transform to local tracking coordinates
+               //Double_t xg =  x[0] * TMath::Cos(fAlpha) + x[1] * TMath::Sin(fAlpha); 
+               Double_t yg = -x[0] * TMath::Sin(fAlpha) + x[1] * TMath::Cos(fAlpha);
+
+               // apply tilt pad correction
+               yc+= (zc - x[2]) * tracklet->GetTilt();
+               
+               Double_t dy = yc-yg;
+
+               TTreeSRedirector &cstreamer = *fOutputStreamer;
+               cstreamer << "ResidualsClustersTrack"
+                       << "c.="   << c
+                       << "dy="   << dy
+                       << "\n";
+       }
+}
+
+//____________________________________________________
+void AliTRDtrackerDebug::ResidualsClustersTracklet(const AliTRDseedV1 *tracklet) const
+{
+// Calculates distances from clusters to tracklets
+       
+       Double_t x0 = tracklet->GetX0(), 
+                y0 = tracklet->GetYfit(0), 
+                ys = tracklet->GetYfit(1);
+                //z0 = tracklet->GetZfit(0), 
+                //zs = tracklet->GetZfit(1);
+       
+       AliTRDcluster *c = 0x0;
+       for(int ic=0; ic<35/*AliTRDseed:knTimebins*/; ic++){
+               if(!(c = tracklet->GetClusters(ic))) continue;
+               Double_t xc = c->GetX(), yc = c->GetY()/*, zc = c->GetZ()*/;
+               Double_t dy = yc- (y0-(x0-xc)*ys);
+
+               //To draw  use : 
+    //ResidualsClustersTracklet->Draw("TMath::Abs(10.*dy):TMath::ATan(ys)*TMath::RadToDeg()>>h(20, -40, 40)", "", "prof");
+               TTreeSRedirector &cstreamer = *fOutputStreamer;
+               cstreamer << "ResidualsClustersTracklet"
+                       << "c.="   << c
+                       << "ys="   << ys
+                       << "dy="   << dy
+                       << "\n";
+       }
+}
+
+//____________________________________________________
+void AliTRDtrackerDebug::ResidualsClustersParametrisation(const AliTRDseedV1 *tracklet) const
+{
+// Calculates distances from clusters to Rieman fit.
+       
+       // store cluster positions
+       Double_t x0 = tracklet->GetX0();
+       AliTRDcluster *c = 0x0;
+       
+       Double_t x[2]; Int_t ncl, mcl, jc;
+       TLinearFitter fitter(3, "hyp2");
+       for(int ic=0; ic<35/*AliTRDseed:knTimebins*/; ic++){
+               if(!(c = tracklet->GetClusters(ic))) continue;
+               Double_t xc = c->GetX(), yc = c->GetY()/*, zc = c->GetZ()*/;
+               
+               jc = ic; ncl = 0; mcl=0; fitter.ClearPoints();
+               while(ncl<6){
+                       // update index
+                       mcl++;
+                       jc = ic + ((mcl&1)?-1:1)*(mcl>>1);
+
+                       if(jc<0 || jc>=35) continue;
+                       if(!(c = tracklet->GetClusters(jc))) continue;
+
+                       x[0] = c->GetX()-x0;
+                       x[1] = x[0]*x[0];
+                       fitter.AddPoint(x, c->GetY(), c->GetSigmaY2());
+                       ncl++;
+               }
+               fitter.Eval();
+               Double_t dy = yc - fitter.GetParameter(0) -fitter.GetParameter(1) * (xc-x0) - fitter.GetParameter(2)* (xc-x0)*(xc-x0); 
+       
+               TTreeSRedirector &cstreamer = *fOutputStreamer;
+               cstreamer << "ResidualsClustersParametrisation"
+                       << "dy="   << dy
+                       << "\n";
+       }
+}
+
+
+//____________________________________________________
+void AliTRDtrackerDebug::ResidualsTrackletsTrack() const
+{
+// Calculates distances from tracklets to the TRD track.
+       
+       if(fTrack->GetNumberOfTracklets() < 6) return;
+
+       // build a working copy of the tracklets attached to the track 
+       // and initialize working variables fX, fY and fZ
+       AliTRDcluster *c = 0x0;
+       AliTRDseedV1 tracklet[6] = {0x0, 0x0, 0x0, 0x0, 0x0, 0x0};
+       const AliTRDseedV1 *ctracklet = 0x0;
+       for(int ip = 0; ip<6; ip++){
+               if(!(ctracklet = fTrack->GetTracklet(ip))) continue;
+               tracklet[ip] = (*ctracklet); 
+               Double_t x0 = tracklet[ip].GetX0();
+               for(int ic=0; ic<35/*AliTRDseed:knTimebins*/; ic++){
+                       if(!(c = tracklet[ip].GetClusters(ic))) continue;
+                       Double_t xc = c->GetX(), yc = c->GetY(), zc = c->GetZ();
+                       tracklet[ip].SetX(ic, xc-x0);
+                       tracklet[ip].SetY(ic, yc);
+                       tracklet[ip].SetZ(ic, zc);
+               }
+       }
+       
+       // Do a Rieman fit (with tilt correction) for all tracklets 
+       // except the one which is tested. 
+       // (Based on AliTRDseed::IsOK() return false)
+       for(int ip=0; ip<6; ip++){
+               // reset tracklet to be tested
+               Double_t x0 = tracklet[ip].GetX0();
+               new(&tracklet[ip]) AliTRDseedV1();
+               tracklet[ip].SetX0(x0);
+
+               // fit Rieman with tilt correction
+               AliTRDseedV1::FitRiemanTilt(&tracklet[0], kTRUE);
+
+               // make a copy of the fit result
+               Double_t 
+                       y0   = tracklet[ip].GetYref(0),
+                       dydx = tracklet[ip].GetYref(1),
+                       z0   = tracklet[ip].GetZref(0),
+                       dzdx = tracklet[ip].GetZref(1);
+
+               // restore tracklet
+               tracklet[ip] = (*fTrack->GetTracklet(ip)); 
+               for(int ic=0; ic<35/*AliTRDseed:knTimebins*/; ic++){
+                       if(!(c = tracklet[ip].GetClusters(ic))) continue;
+                       Double_t xc = c->GetX(), yc = c->GetY(), zc = c->GetZ();
+                       tracklet[ip].SetX(ic, xc-x0);
+                       tracklet[ip].SetY(ic, yc);
+                       tracklet[ip].SetZ(ic, zc);
+               }               
+               
+               // fit clusters
+               AliTRDseedV1 ts(tracklet[ip]);
+               ts.SetYref(0, y0); ts.SetYref(1, dydx);
+               ts.SetZref(0, z0); ts.SetZref(1, dzdx);
+               ts.Update();
+
+               // save results for plotting
+               Int_t plane   = tracklet[ip].GetPlane();
+               Double_t dy   = tracklet[ip].GetYfit(0) - ts.GetYfit(0);
+               Double_t tgy  = tracklet[ip].GetYfit(1);
+               Double_t dtgy = (tracklet[ip].GetYfit(1) - ts.GetYfit(1))/(1. + tracklet[ip].GetYfit(1) * ts.GetYfit(1));
+               Double_t dz   = tracklet[ip].GetZfit(0) - ts.GetZfit(0);
+               Double_t tgz  = tracklet[ip].GetZfit(1);
+               Double_t dtgz = (tracklet[ip].GetZfit(1) - ts.GetZfit(1))/(1. + tracklet[ip].GetZfit(1) * ts.GetZfit(1));
+               TTreeSRedirector &cstreamer = *fOutputStreamer;
+               cstreamer << "ResidualsTrackletsTrack"
+                       << "ref.="   << &tracklet[ip]
+                       << "fit.="   << &ts
+                       << "plane="  << plane
+                       << "dy="     << dy
+                       << "tgy="    << tgy
+                       << "dtgy="   << dtgy
+                       << "dz="     << dz
+                       << "tgz="    << tgz
+                       << "dtgz="   << dtgz
+                       << "\n";
+       }
+}
+
diff --git a/TRD/AliTRDtrackerDebug.h b/TRD/AliTRDtrackerDebug.h
new file mode 100644 (file)
index 0000000..8f6ac60
--- /dev/null
@@ -0,0 +1,62 @@
+#ifndef ALITRDTRACKERDEBUG_H
+#define ALITRDTRACKERDEBUG_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliTRDtrackerDebug.h 22646 2007-11-29 18:13:40Z cblume $ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+//  Reader for the TRD tracker debug streamer                             // 
+//                                                                        // 
+//  Authors:                                                              //
+//                                                                        //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                     //
+//                                                                        // 
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALITRDTRACKERV1_H
+#include "AliTRDtrackerV1.h"
+#endif
+
+class TTree;
+class TTreeSRedirector;
+class AliTRDtrackV1;
+class AliTRDseedV1;
+
+class AliTRDtrackerDebug : public AliTRDtrackerV1
+{
+public:
+       AliTRDtrackerDebug();
+       ~AliTRDtrackerDebug();
+
+       void                            Draw(const Option_t *);
+
+       Bool_t      Init();
+       Bool_t      Open(const char *method);
+       Int_t       Process();
+
+       void        ResidualsClustersTrack(const AliTRDseedV1 *tracklet);
+       void        ResidualsClustersTracklet(const AliTRDseedV1 *tracklet) const;
+       void        ResidualsClustersParametrisation(const AliTRDseedV1 *tracklet) const;
+       void        ResidualsTrackletsTrack() const;
+
+
+private:
+       AliTRDtrackerDebug(const AliTRDtrackerDebug &);
+       AliTRDtrackerDebug& operator=(const AliTRDtrackerDebug &);
+
+
+  TTreeSRedirector *fOutputStreamer;                 //!Output streamer
+       TTree            *fTree;       // debug tree
+       AliTRDseedV1     *fTracklet;   // current tracklet
+       AliTRDtrackV1    *fTrack;      // current TRD track
+       Int_t            fNClusters;   // N clusters for current track
+       Float_t          fAlpha;       // sector
+
+       ClassDef(AliTRDtrackerDebug, 1) // debug suite of the TRD tracker
+};
+
+#endif
+
index 26b4f0c..c2fe6d8 100644 (file)
@@ -36,7 +36,6 @@
 #include <TH1D.h>
 #include <TH2D.h>
 #include <TLinearFitter.h>
-#include <TObjArray.h> 
 #include <TROOT.h>
 #include <TTree.h>  
 #include <TClonesArray.h>
@@ -49,8 +48,8 @@
 #include "AliRieman.h"
 #include "AliTrackPointArray.h"
 
-#include "AliTRDtracker.h"
 #include "AliTRDtrackerV1.h"
+#include "AliTRDtrackingChamber.h"
 #include "AliTRDgeometry.h"
 #include "AliTRDpadPlane.h"
 #include "AliTRDgeometry.h"
 #include "AliTRDCommonParam.h"
 #include "AliTRDReconstructor.h"
 #include "AliTRDCalibraFillHisto.h"
-#include "AliTRDtrackerFitter.h"
-#include "AliTRDstackLayer.h"
+#include "AliTRDchamberTimeBin.h"
 #include "AliTRDrecoParam.h"
 #include "AliTRDseedV1.h"
 #include "AliTRDtrackV1.h"
 #include "Cal/AliTRDCalDet.h"
 
-#define DEBUG
 
 ClassImp(AliTRDtrackerV1)
+
+
+const  Float_t  AliTRDtrackerV1::fgkMinClustersInTrack =  0.5;  //
+const  Float_t  AliTRDtrackerV1::fgkLabelFraction      =  0.8;  //
+const  Double_t AliTRDtrackerV1::fgkMaxChi2            = 12.0;  //
+const  Double_t AliTRDtrackerV1::fgkMaxSnp             =  0.95; // Maximum local sine of the azimuthal angle
+const  Double_t AliTRDtrackerV1::fgkMaxStep            =  2.0;  // Maximal step size in propagation 
 Double_t AliTRDtrackerV1::fgTopologicQA[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
 };
+TTreeSRedirector *AliTRDtrackerV1::fgDebugStreamer = 0x0;
+AliRieman* AliTRDtrackerV1::fgRieman = 0x0;
+TLinearFitter* AliTRDtrackerV1::fgTiltedRieman = 0x0;
+TLinearFitter* AliTRDtrackerV1::fgTiltedRiemanConstrained = 0x0;
 
 //____________________________________________________________________
-AliTRDtrackerV1::AliTRDtrackerV1(AliTRDrecoParam *p) 
-  :AliTRDtracker()
-  ,fSieveSeeding(0)
+AliTRDtrackerV1::AliTRDtrackerV1() 
+  :AliTracker()
+  ,fGeom(new AliTRDgeometry())
+  ,fClusters(0x0)
   ,fTracklets(0x0)
-  ,fRecoParam(p)
-  ,fFitter(0x0)
-{
-  //
-  // Default constructor. Nothing is initialized.
-  //
-
-}
-
-//____________________________________________________________________
-AliTRDtrackerV1::AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p) 
-  :AliTRDtracker(in)
+  ,fTracks(0x0)
+  ,fTimeBinsPerPlane(0)
   ,fSieveSeeding(0)
-  ,fTracklets(0x0)
-  ,fRecoParam(p)
-  ,fFitter(0x0)
 {
   //
-  // Standard constructor.
-  // Setting of the geometry file, debug output (if enabled)
-  // and initilize fitter helper.
-  //
-
-       fFitter = new AliTRDtrackerFitter();
-
-#ifdef DEBUG
-       fFitter->SetDebugStream(fDebugStreamer);
-#endif
+  // Default constructor.
+  // 
+  if (!AliTRDcalibDB::Instance()) {
+    AliFatal("Could not get calibration object");
+  }
+  fTimeBinsPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
 
-}
+       for (Int_t isector = 0; isector < AliTRDgeometry::kNsect; isector++) new(&fTrSec[isector]) AliTRDtrackingSector(fGeom, isector, fTimeBinsPerPlane);
   
+  if(AliTRDReconstructor::StreamLevel() > 1){
+               TDirectory *savedir = gDirectory; 
+               fgDebugStreamer    = new TTreeSRedirector("TRD.TrackerDebug.root");
+       savedir->cd();
+       }
+}
+
 //____________________________________________________________________
 AliTRDtrackerV1::~AliTRDtrackerV1()
 { 
   //
   // Destructor
   //
-
-       if(fFitter) delete fFitter;
-       if(fRecoParam) delete fRecoParam;
+       
+       if(fgDebugStreamer) delete fgDebugStreamer;
+       if(fgRieman) delete fgRieman;
+       if(fgTiltedRieman) delete fgTiltedRieman;
+       if(fgTiltedRiemanConstrained) delete fgTiltedRiemanConstrained;
+       if(fTracks) {fTracks->Delete(); delete fTracks;}
        if(fTracklets) {fTracklets->Delete(); delete fTracklets;}
+       if(fClusters) {fClusters->Delete(); delete fClusters;}
+       if(fGeom) delete fGeom;
 }
 
 //____________________________________________________________________
@@ -143,34 +147,96 @@ Int_t AliTRDtrackerV1::Clusters2Tracks(AliESDEvent *esd)
   //    See AliTRDtrackerV1::Clusters2TracksSM() for details.
   //
 
-       if(!fRecoParam){
-               AliError("Reconstruction configuration not initialized. Call first AliTRDtrackerV1::SetRecoParam().");
+       if(!AliTRDReconstructor::RecoParam()){
+               AliError("Reconstruction configuration not initialized. Call first AliTRDReconstructor::SetRecoParam().");
                return 0;
        }
        
        //AliInfo("Start Track Finder ...");
        Int_t ntracks = 0;
-       for(int ism=0; ism<AliTRDtracker::kTrackingSectors; ism++){
+       for(int ism=0; ism<AliTRDgeometry::kNsect; ism++){
+//     for(int ism=1; ism<2; ism++){
                        //AliInfo(Form("Processing supermodule %i ...", ism));
-                       ntracks += Clusters2TracksSM(fTrSec[ism], esd);
+                       ntracks += Clusters2TracksSM(ism, esd);
        }
-       AliInfo(Form("Found %d TRD tracks.", ntracks));
+  AliInfo(Form("Number of found tracks : %d", ntracks));
        return ntracks;
 }
 
 
 //_____________________________________________________________________________
-Bool_t AliTRDtrackerV1::GetTrackPoint(Int_t index, AliTrackPoint &/*p*/) const
+Bool_t AliTRDtrackerV1::GetTrackPoint(Int_t index, AliTrackPoint &p) const
 {
        //AliInfo(Form("Asking for tracklet %d", index));
        
        if(index<0) return kFALSE;
-       //AliTRDseedV1 *tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index);
-       // etc
+       AliTRDseedV1 *tracklet = 0x0; 
+       if(!(tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index))) return kFALSE;
+       
+       // get detector for this tracklet
+       AliTRDcluster *cl = 0x0;
+       Int_t ic = 0; do; while(!(cl = tracklet->GetClusters(ic++)));    
+       Int_t  idet     = cl->GetDetector();
+               
+       Double_t local[3];
+       local[0] = tracklet->GetX0(); 
+       local[1] = tracklet->GetYfit(0);
+       local[2] = tracklet->GetZfit(0);
+       Double_t global[3];
+       fGeom->RotateBack(idet, local, global);
+       p.SetXYZ(global[0],global[1],global[2]);
+       
+       
+       // setting volume id
+       AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
+       switch (fGeom->GetPlane(idet)) {
+       case 0:
+               iLayer = AliGeomManager::kTRD1;
+               break;
+       case 1:
+               iLayer = AliGeomManager::kTRD2;
+               break;
+       case 2:
+               iLayer = AliGeomManager::kTRD3;
+               break;
+       case 3:
+               iLayer = AliGeomManager::kTRD4;
+               break;
+       case 4:
+               iLayer = AliGeomManager::kTRD5;
+               break;
+       case 5:
+               iLayer = AliGeomManager::kTRD6;
+               break;
+       };
+       Int_t    modId = fGeom->GetSector(idet) * fGeom->Ncham() + fGeom->GetChamber(idet);
+       UShort_t volid = AliGeomManager::LayerToVolUID(iLayer, modId);
+       p.SetVolumeID(volid);
+               
        return kTRUE;
 }
 
+//____________________________________________________________________
+TLinearFitter* AliTRDtrackerV1::GetTiltedRiemanFitter()
+{
+       if(!fgTiltedRieman) fgTiltedRieman = new TLinearFitter(4, "hyp4");
+       return fgTiltedRieman;
+}
 
+//____________________________________________________________________
+TLinearFitter* AliTRDtrackerV1::GetTiltedRiemanFitterConstraint()
+{
+       if(!fgTiltedRiemanConstrained) fgTiltedRiemanConstrained = new TLinearFitter(2, "hyp2");
+       return fgTiltedRiemanConstrained;
+}
+       
+//____________________________________________________________________ 
+AliRieman* AliTRDtrackerV1::GetRiemanFitter()
+{
+       if(!fgRieman) fgRieman = new AliRieman(AliTRDtrackingChamber::kNTimeBins * AliTRDgeometry::kNplan);
+       return fgRieman;
+}
+       
 //_____________________________________________________________________________
 Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) 
 {
@@ -222,7 +288,6 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
                //track->Print();
                track->SetSeedLabel(lbl);
                seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); // Make backup
-               fNseeds++;
                Float_t p4          = track->GetC();
                Int_t   expectedClr = FollowBackProlongation(*track);
                //AliInfo(Form("\nTRACK %d Clusters %d [%d] in chi2 %f", index[iSeed], expectedClr, track->GetNumberOfClusters(), track->GetChi2()));
@@ -294,7 +359,7 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
        
                /**/
                // Debug part of tracking
-/*             TTreeSRedirector &cstream = *fDebugStreamer;
+/*             TTreeSRedirector &cstream = *fgDebugStreamer;
                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()) {
@@ -397,25 +462,13 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
        
                seed->SetTRDQuality(track->StatusForTOF());
                seed->SetTRDBudget(track->GetBudget(0));
-               
-//             if ((seed->GetStatus()&AliESDtrack::kTRDin)!=0 ) printf("TRDin ");      
-//             if ((seed->GetStatus()&AliESDtrack::kTRDbackup)!=0 ) printf("TRDbackup ");      
-//             if ((seed->GetStatus()&AliESDtrack::kTRDStop)!=0 ) printf("TRDstop ");  
-//             if ((seed->GetStatus()&AliESDtrack::kTRDout)!=0 ) printf("TRDout ");    
-//             printf("\n");
                delete track;
-
-               //seed->GetExternalCovariance(cov);
-               //AliInfo(Form("track %d cov[%f %f] 2", index[iSeed], cov[0], cov[2]));
        }
        
 
-       AliInfo(Form("Number of seeds: %d",fNseeds));
-       AliInfo(Form("Number of back propagated TRD tracks: %d",found));
-               
-       //fSeeds->Clear(); 
-       fNseeds = 0;
-       
+       AliInfo(Form("Number of seeds: %d", nSeed));
+       AliInfo(Form("Number of back propagated TRD tracks: %d", found));
+                       
        delete [] index;
        delete [] quality;
        
@@ -526,53 +579,39 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
   
        //AliInfo("");
        Int_t    nClustersExpected = 0;
-       Float_t clength = AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick();
        Int_t lastplane = 5; //GetLastPlane(&t);
        for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
-               //AliInfo(Form("plane %d", iplane));
-    Int_t row1 = GetGlobalTimeBin(0, iplane, 0); // to be modified to the true time bin in the geometrical acceptance
-               //AliInfo(Form("row1 %d", row1));
-
-    // Propagate track close to the plane if neccessary
-               AliTRDpropagationLayer *layer = fTrSec[0]->GetLayer(row1);
-               Double_t currentx  = layer->GetX();
-    if (currentx < (-fgkMaxStep + t.GetX())) 
-      if (!PropagateToX(t, currentx+fgkMaxStep, fgkMaxStep)) break;
-
+    Int_t   index   = 0;
+    AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
+               if(!tracklet) continue;
+               if(!tracklet->IsOK()) AliWarning("tracklet not OK");
+               
+               t.SetTracklet(tracklet, iplane, index);
+               
+               Double_t x  = tracklet->GetX0();
+    if (x < (t.GetX()-fgkMaxStep) && !PropagateToX(t, x+fgkMaxStep, fgkMaxStep)) break;
     if (!AdjustSector(&t)) break;
      
-               Int_t row0    = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
-               //AliInfo(Form("row0 %d", row0));
-
     // Start global position
     Double_t xyz0[3];
     t.GetXYZ(xyz0);
 
                // End global position
-    Double_t x = fTrSec[0]->GetLayer(row0)->GetX(), y, z;
+    Double_t alpha = t.GetAlpha(), y, z;
     if (!t.GetProlongation(x,y,z)) break;    
     Double_t xyz1[3];
-    xyz1[0] =  x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
-    xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
+    xyz1[0] =  x * TMath::Cos(alpha) - y * TMath::Sin(alpha);
+    xyz1[1] =  x * TMath::Sin(alpha) + y * TMath::Cos(alpha);
     xyz1[2] =  z;
-               
+                               
     // Get material budget
     Double_t param[7];
-    AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
+    AliTracker::MeanMaterialBudget(xyz0, xyz1, param);
     Double_t xrho= param[0]*param[4];
     Double_t xx0 = param[1]; // Get mean propagation parameters
 
-    // Propagate and update
-    //Int_t sector = t.GetSector();
-    Int_t   index   = 0;
-               //AliInfo(Form("sector %d", sector));
-    AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
-         //AliInfo(Form("tracklet %p @ %d", tracklet, index));
-               if(!tracklet) continue;
-               //AliInfo(Form("reco %p", tracklet->GetRecoParam()));
-               t.SetTracklet(tracklet, iplane, index);
-               
-               t.PropagateTo(tracklet->GetX0() - clength, xx0, xrho);
+    // Propagate and update            
+               t.PropagateTo(x, xx0, xrho);
          if (!AdjustSector(&t)) break;
          
     Double_t maxChi2 = t.GetPredictedChi2(tracklet);
@@ -581,7 +620,6 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
        }
   }
 
-#ifdef DEBUG
        if(AliTRDReconstructor::StreamLevel() > 1){
                Int_t index;
                for(int iplane=0; iplane<6; iplane++){
@@ -590,13 +628,12 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
                        t.SetTracklet(tracklet, iplane, index);
                }
 
-               TTreeSRedirector &cstreamer = *fDebugStreamer;
+               TTreeSRedirector &cstreamer = *fgDebugStreamer;
                cstreamer << "FollowProlongation"
                        << "ncl="      << nClustersExpected
                        << "track.="   << &t
                        << "\n";
        }
-#endif
 
   return nClustersExpected;
 
@@ -630,121 +667,80 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
 //
 
        Int_t nClustersExpected = 0;
+  Double_t clength = AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick();
+  AliTRDtrackingChamber *chamber = 0x0;
+  
   // Loop through the TRD planes
   for (Int_t iplane = 0; iplane < AliTRDgeometry::Nplan(); iplane++) {
-               //AliInfo(Form("Processing plane %d ...", iplane));
-               // Get the global time bin for the first local time bin of the given plane
-    Int_t row0 = GetGlobalTimeBin(0, iplane, fTimeBinsPerPlane-1);
-
-               // Retrive first propagation layer in the chamber
-               AliTRDpropagationLayer *layer = fTrSec[0]->GetLayer(row0);
+               // BUILD TRACKLET IF NOT ALREADY BUILT
+               Double_t x = 0., y, z, alpha;
+    AliTRDseedV1 tracklet(*t.GetTracklet(iplane));
+               if(!tracklet.IsOK()){
+               alpha = t.GetAlpha();
+               Int_t sector = Int_t(alpha/AliTRDgeometry::GetAlpha() + (alpha>0. ? 0 : AliTRDgeometry::kNsect));
 
-               // Get the X coordinates of the propagation layer for the first time bin
-    Double_t currentx = layer->GetX(); // what if X is not defined ???
-    if (currentx < t.GetX()) continue;
-    
-               // Get the global time bin for the last local time bin of the given plane
-               Int_t row1 = GetGlobalTimeBin(0, iplane, 0);
+               if(!fTrSec[sector].GetNChambers()) continue;
+               
+               if((x = fTrSec[sector].GetX(iplane)) < 1.) continue;
+               
+                       if (!t.GetProlongation(x, y, z)) break;
+                       Int_t stack = fGeom->GetChamber(z, iplane);
+                       Int_t nCandidates = stack >= 0 ? 1 : 2;
+                       z -= stack >= 0 ? 0. : 4.; 
+                       
+                       for(int icham=0; icham<nCandidates; icham++, z+=8){
+                               if((stack = fGeom->GetChamber(z, iplane)) < 0) continue;
+                       
+                               if(!(chamber = fTrSec[sector].GetChamber(stack, iplane))) continue;
+                       
+                               if(chamber->GetNClusters() < fTimeBinsPerPlane*AliTRDReconstructor::AliTRDReconstructor::RecoParam()->GetFindableClusters()) continue;
+                       
+                               x = chamber->GetX();
+                       
+                               AliTRDpadPlane *pp = fGeom->GetPadPlane(iplane, stack);
+                               tracklet.SetTilt(TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle()));
+                               tracklet.SetPadLength(pp->GetLengthIPad());
+                               tracklet.SetPlane(iplane);
+                               tracklet.SetX0(x);
+                               tracklet.Init(&t);
+                               if(!tracklet.AttachClustersIter(chamber, 1000.)) continue;
+                               tracklet.Init(&t);
+               
+               if(tracklet.GetN() < fTimeBinsPerPlane * AliTRDReconstructor::RecoParam()->GetFindableClusters()) continue;
+                       
+                               break;
+                       }
+               }
+    if(!tracklet.IsOK()){
+                       if(x < 1.) continue; //temporary
+                       if(!PropagateToX(t, x-fgkMaxStep, fgkMaxStep)) break;
+                       if(!AdjustSector(&t)) break;
+                       if(TMath::Abs(t.GetSnp()) > fgkMaxSnp) break;
+       continue;
+    }
     
                // Propagate closer to the current chamber if neccessary 
-    if (currentx > (fgkMaxStep + t.GetX()) && !PropagateToX(t, currentx-fgkMaxStep, fgkMaxStep)) break;
-   
-               // Rotate track to adjacent sector if neccessary
+    x -= clength;
+    if (x > (fgkMaxStep + t.GetX()) && !PropagateToX(t, x-fgkMaxStep, fgkMaxStep)) break;
     if (!AdjustSector(&t)) break;
-               Int_t sector =  Int_t(TMath::Abs(t.GetAlpha()/AliTRDgeometry::GetAlpha()));
-               if(t.GetAlpha() < 0) sector = AliTRDgeometry::Nsect() - sector-1;
-               
-               //AliInfo(Form("sector %d [%f]", sector, t.GetAlpha()));
-
-               // Check whether azimuthal angle is getting too large
     if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) break;
-   
-    //Calculate global entry and exit positions of the track in chamber (only track prolongation)
-    Double_t xyz0[3]; // entry point 
-               t.GetXYZ(xyz0);   
-               //printf("Entry global x[%7.3f] y[%7.3f] z[%7.3f]\n", xyz0[0], xyz0[1], xyz0[2]);
-                               
-               // Get local Y and Z at the X-position of the end of the chamber
-               Double_t x0 = fTrSec[sector]->GetLayer(row1)->GetX(), y, z;
-               if (!t.GetProlongation(x0, y, z)) break;
-               //printf("Exit  local x[%7.3f] y[%7.3f] z[%7.3f]\n", x0, y, z);
-
-               Double_t xyz1[3]; // exit point
-               xyz1[0] =  x0 * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha()); 
-    xyz1[1] = +x0 * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
-    xyz1[2] =  z;
-
-               //printf("Exit  global x[%7.3f] y[%7.3f] z[%7.3f]\n", xyz1[0], xyz1[1], xyz1[2]);
-               // Find tracklet along the path inside the chamber
-    AliTRDseedV1 tracklet(*t.GetTracklet(iplane));
-               // if the track is not already build (e.g. stand alone tracker) we build it now.
-               if(!tracklet.GetN()){ // a better check has to be implemented TODO!!!!!!!
-               
-                       //AliInfo(Form("Building tracklet for plane %d ...", iplane));
-                       // check if we are inside detection volume
-                       Int_t ichmb = fGeom->GetChamber(xyz0[2], iplane); 
-                       if(ichmb<0) ichmb = fGeom->GetChamber(xyz1[2], iplane);
-                       if(ichmb<0){
-                               // here we should decide what to do with the track. The space between the pads in 2 chambers is 4cm+. Is it making sense to continue building the tracklet here TODO????
-                               AliWarning(Form("Track prolongated in the interspace between TRD detectors in plane %d. Skip plane. To be fixed !", iplane));
-                               continue;
-                       }
-       
-                       // temporary until the functionalities of AliTRDpropagationLayer and AliTRDstackLayer are merged TODO
-                       AliTRDpadPlane *pp = fGeom->GetPadPlane(iplane, ichmb);
-                       Int_t nrows = pp->GetNrows();
-                       Double_t stacklength = pp->GetRow0ROC() - pp->GetRowEndROC();/*(nrows - 2) * pp->GetLengthIPad()        + 2 * pp->GetLengthOPad() + (nrows - 1) * pp->GetRowSpacing();*/
-                       Double_t z0  = fGeom->GetRow0(iplane, ichmb, 0);
-
-                       Int_t nClustersChmb = 0;
-                       AliTRDstackLayer stackLayer[35];
-                       for(int itb=0; itb<fTimeBinsPerPlane; itb++){
-                               const AliTRDpropagationLayer ksmLayer(*(fTrSec[sector]->GetLayer(row1 - itb)));
-                               stackLayer[itb] = ksmLayer;
-#ifdef DEBUG
-                               stackLayer[itb].SetDebugStream(fDebugStreamer);
-#endif                 
-                               stackLayer[itb].SetRange(z0 - stacklength, stacklength);
-                               stackLayer[itb].SetSector(sector);
-                               stackLayer[itb].SetStackNr(ichmb);
-                               stackLayer[itb].SetNRows(nrows);
-                               stackLayer[itb].SetRecoParam(fRecoParam);
-                               stackLayer[itb].BuildIndices();
-                               nClustersChmb += stackLayer[itb].GetNClusters();
-                       }
-                       if(!nClustersChmb) continue;
-                       //AliInfo(Form("Detector p[%d] c[%d]. Building tracklet from %d clusters ... ", iplane, ichmb, nClustersChmb));
-
-                       tracklet.SetRecoParam(fRecoParam);
-                       tracklet.SetTilt(TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle()));
-                       tracklet.SetPadLength(pp->GetLengthIPad());
-                       tracklet.SetPlane(iplane);
-                       //Int_t tbRange   = fTimeBinsPerPlane; //Int_t(AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght() * AliTRDCommonParam::Instance()->GetSamplingFrequency()/AliTRDcalibDB::Instance()->GetVdriftDet()->GetValue(det));
-                       //printf("%d hl[%f] pl[%f] tb[%d]\n", il, hL[il], padlength[il], tbRange[il]);
-                       //tracklet.SetNTimeBinsRange(tbRange);
-                       tracklet.SetX0(x0);
-                       tracklet.Init(&t);
-                       if(!tracklet.AttachClustersIter(stackLayer, 1000.)) continue;
-                       tracklet.Init(&t);
-
-                       //if(!tracklet.AttachClusters(stackLayer, kTRUE)) continue;
-                       //if(!tracklet.Fit()) continue;
-               }
-               Int_t ncl = tracklet.GetN();
-               //AliInfo(Form("N clusters %d", ncl));
-               
-               // Discard tracklet if bad quality.
-               //Check if this condition is not already checked during building of the tracklet
-    if(ncl < fTimeBinsPerPlane * fRecoParam->GetFindableClusters()){
-                       //AliInfo(Form("Discard tracklet for %d nclusters", ncl));
-                       continue;
-               }
                
                // load tracklet to the tracker and the track
                Int_t index = SetTracklet(&tracklet);
                t.SetTracklet(&tracklet, iplane, index);
-               
+   
+   
                // Calculate the mean material budget along the path inside the chamber
+    //Calculate global entry and exit positions of the track in chamber (only track prolongation)
+    Double_t xyz0[3]; // entry point 
+               t.GetXYZ(xyz0);
+               alpha = t.GetAlpha();
+               x = tracklet.GetX0();
+               if (!t.GetProlongation(x, y, z)) break;
+               Double_t xyz1[3]; // exit point
+               xyz1[0] =  x * TMath::Cos(alpha) - y * TMath::Sin(alpha); 
+    xyz1[1] = +x * TMath::Sin(alpha) + y * TMath::Cos(alpha);
+    xyz1[2] =  z;
     Double_t param[7];
                AliTracker::MeanMaterialBudget(xyz0, xyz1, param);      
     // The mean propagation parameters
@@ -752,19 +748,19 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
     Double_t xx0  = param[1]; // radiation length
                
                // Propagate and update track
-               t.PropagateTo(tracklet.GetX0(), xx0, xrho);
+               t.PropagateTo(x, xx0, xrho);
          if (!AdjustSector(&t)) break;
                Double_t maxChi2 = t.GetPredictedChi2(&tracklet);
                if (maxChi2<1e+10 && t.Update(&tracklet, maxChi2)){ 
-                       nClustersExpected += ncl;
+                       nClustersExpected += tracklet.GetN();
                }
                // Reset material budget if 2 consecutive gold
-               if(iplane>0 && ncl + t.GetTracklet(iplane-1)->GetN() > 20) t.SetBudget(2, 0.);
+               if(iplane>0 && tracklet.GetN() + t.GetTracklet(iplane-1)->GetN() > 20) t.SetBudget(2, 0.);
 
                // Make backup of the track until is gold
                // TO DO update quality check of the track.
                // consider comparison with fTimeBinsRange
-               Float_t ratio0 = ncl / Float_t(fTimeBinsPerPlane);
+               Float_t ratio0 = tracklet.GetN() / Float_t(fTimeBinsPerPlane);
                //Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);    
     //printf("tracklet.GetChi2() %f     [< 18.0]\n", tracklet.GetChi2()); 
                //printf("ratio0    %f              [>   0.8]\n", ratio0);
@@ -784,19 +780,391 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
                
        } // end planes loop
 
-#ifdef DEBUG
        if(AliTRDReconstructor::StreamLevel() > 1){
-               TTreeSRedirector &cstreamer = *fDebugStreamer;
+               TTreeSRedirector &cstreamer = *fgDebugStreamer;
                cstreamer << "FollowBackProlongation"
                        << "ncl="      << nClustersExpected
                        << "track.="   << &t
                        << "\n";
        }
-#endif
        
        return nClustersExpected;
 }
 
+//_________________________________________________________________________
+Float_t AliTRDtrackerV1::FitRieman(AliTRDseedV1 *tracklets, Double_t *chi2, Int_t *planes){
+//
+// Fits a Riemann-circle to the given points without tilting pad correction.
+// The fit is performed using an instance of the class AliRieman (equations 
+// and transformations see documentation of this class)
+// Afterwards all the tracklets are Updated
+//
+// Parameters: - Array of tracklets (AliTRDseedV1)
+//             - Storage for the chi2 values (beginning with direction z)  
+//             - Seeding configuration
+// Output:     - The curvature
+//
+  AliRieman *fitter = AliTRDtrackerV1::GetRiemanFitter();
+       fitter->Reset();
+  Int_t allplanes[] = {0, 1, 2, 3, 4, 5};
+  Int_t *ppl = &allplanes[0];
+       Int_t maxLayers = 6;
+  if(planes){
+    maxLayers = 4;
+    ppl = planes;
+  }
+  for(Int_t il = 0; il < maxLayers; il++){
+               if(!tracklets[ppl[il]].IsOK()) continue;
+    fitter->AddPoint(tracklets[ppl[il]].GetX0(), tracklets[ppl[il]].GetYfitR(0), tracklets[ppl[il]].GetZProb(),1,10);
+  }
+  fitter->Update();
+  // Set the reference position of the fit and calculate the chi2 values
+  memset(chi2, 0, sizeof(Double_t) * 2);
+       for(Int_t il = 0; il < maxLayers; il++){
+               // Reference positions
+               tracklets[ppl[il]].Init(fitter);
+               
+               // chi2
+               if((!tracklets[ppl[il]].IsOK()) && (!planes)) continue;
+               chi2[0] += tracklets[ppl[il]].GetChi2Z();
+               chi2[1] += tracklets[ppl[il]].GetChi2Y();
+       }
+       return fitter->GetC();
+}
+
+//_________________________________________________________________________
+void AliTRDtrackerV1::FitRieman(AliTRDcluster **seedcl, Double_t chi2[2])
+{
+//
+// Performs a Riemann helix fit using the seedclusters as spacepoints
+// Afterwards the chi2 values are calculated and the seeds are updated
+//
+// Parameters: - The four seedclusters
+//             - The tracklet array (AliTRDseedV1)
+//             - The seeding configuration
+//             - Chi2 array
+//
+// debug level 2
+//
+       AliRieman *fitter = AliTRDtrackerV1::GetRiemanFitter();
+       fitter->Reset();
+       for(Int_t i = 0; i < 4; i++)
+               fitter->AddPoint(seedcl[i]->GetX(), seedcl[i]->GetY(), seedcl[i]->GetZ(), 1, 10);
+       fitter->Update();
+       
+       
+       // Update the seed and calculated the chi2 value
+       chi2[0] = 0; chi2[1] = 0;
+       for(Int_t ipl = 0; ipl < kNSeedPlanes; ipl++){
+               // chi2
+               chi2[0] += (seedcl[ipl]->GetZ() - fitter->GetZat(seedcl[ipl]->GetX())) * (seedcl[ipl]->GetZ() - fitter->GetZat(seedcl[ipl]->GetX()));
+               chi2[1] += (seedcl[ipl]->GetY() - fitter->GetYat(seedcl[ipl]->GetX())) * (seedcl[ipl]->GetY() - fitter->GetYat(seedcl[ipl]->GetX()));
+       }       
+}
+
+
+//_________________________________________________________________________
+Float_t AliTRDtrackerV1::FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Double_t zVertex)
+{
+//
+// Fits a helix to the clusters. Pad tilting is considered. As constraint it is 
+// assumed that the vertex position is set to 0.
+// This method is very usefull for high-pt particles
+// Basis for the fit: (x - x0)^2 + (y - y0)^2 - R^2 = 0
+//      x0, y0: Center of the circle
+// Measured y-position: ymeas = y - tan(phiT)(zc - zt)
+//      zc: center of the pad row
+// Equation which has to be fitted (after transformation):
+// a + b * u + e * v + 2*(ymeas + tan(phiT)(z - zVertex))*t = 0
+// Transformation:
+// t = 1/(x^2 + y^2)
+// u = 2 * x * t
+// v = 2 * x * tan(phiT) * t
+// Parameters in the equation: 
+//    a = -1/y0, b = x0/y0, e = dz/dx
+//
+// The Curvature is calculated by the following equation:
+//               - curv = a/Sqrt(b^2 + 1) = 1/R
+// Parameters:   - the 6 tracklets
+//               - the Vertex constraint
+// Output:       - the Chi2 value of the track
+//
+// debug level 5
+//
+
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       Int_t nTimeBins = cal->GetNumberOfTimeBins();
+       
+       TLinearFitter *fitter = GetTiltedRiemanFitterConstraint();
+       fitter->StoreData(kTRUE);
+       fitter->ClearPoints();
+
+       Float_t x, y, z, w, t, error, tilt;
+       Double_t uvt[2];
+       Int_t nPoints = 0;
+  for(Int_t ipl = 0; ipl < AliTRDgeometry::kNplan; ipl++){
+               if(!tracklets[ipl].IsOK()) continue;
+               for(Int_t itb = 0; itb < nTimeBins; itb++){
+                       if(!tracklets[ipl].IsUsable(itb)) continue;
+                       x = tracklets[ipl].GetX(itb) + tracklets[ipl].GetX0();
+                       y = tracklets[ipl].GetY(itb);
+                       z = tracklets[ipl].GetZ(itb);
+                       tilt = tracklets[ipl].GetTilt();
+                       // Transformation
+                       t = 1/(x * x + y * y);
+                       uvt[0] = 2 * x* t;      
+                       uvt[1] =  2.0 * tilt * x * t;
+                       w = 2.0 * (y + tilt * (z - zVertex)) * t;
+                       error = 2 * 0.2 * t;
+                       fitter->AddPoint(uvt, w, error);
+                       nPoints++;
+               }
+       }
+       fitter->Eval();
+
+       // Calculate curvature
+       Double_t a = fitter->GetParameter(0);
+       Double_t b = fitter->GetParameter(0);
+       Double_t curvature = a/TMath::Sqrt(b*b + 1);
+
+       Float_t chi2track = fitter->GetChisquare()/Double_t(nPoints);
+       for(Int_t ip = 0; ip < AliTRDtrackerV1::kNPlanes; ip++)
+               tracklets[ip].SetCC(curvature);
+
+       if(AliTRDReconstructor::StreamLevel() >= 5){
+               //Linear Model on z-direction
+         Double_t xref = (tracklets[2].GetX0() + tracklets[3].GetX0())/2;              // Relative to the middle of the stack
+               Double_t slope = fitter->GetParameter(2);
+               Double_t zref = slope * xref;
+               Float_t chi2Z = CalculateChi2Z(tracklets, zref, slope);
+               TTreeSRedirector &treeStreamer = *fgDebugStreamer;
+               treeStreamer << "FitTiltedRiemanConstraint"
+                       << "Curvature=" << curvature
+                       << "Chi2Track=" << chi2track
+                       << "Chi2Z="                     << chi2Z
+                       << "zref="                      << zref
+                       << "\n";
+       }
+       return chi2track;
+}
+
+//_________________________________________________________________________
+Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigError)
+{
+//
+// Performs a Riemann fit taking tilting pad correction into account
+// The equation of a Riemann circle, where the y position is substituted by the 
+// measured y-position taking pad tilting into account, has to be transformed
+// into a 4-dimensional hyperplane equation
+// Riemann circle: (x-x0)^2 + (y-y0)^2 -R^2 = 0
+// Measured y-Position: ymeas = y - tan(phiT)(zc - zt)
+//          zc: center of the pad row
+//          zt: z-position of the track
+// The z-position of the track is assumed to be linear dependent on the x-position
+// Transformed equation: a + b * u + c * t + d * v  + e * w - 2 * (ymeas + tan(phiT) * zc) * t = 0
+// Transformation:       u = 2 * x * t
+//                       v = 2 * tan(phiT) * t
+//                       w = 2 * tan(phiT) * (x - xref) * t
+//                       t = 1 / (x^2 + ymeas^2)
+// Parameters:           a = -1/y0
+//                       b = x0/y0
+//                       c = (R^2 -x0^2 - y0^2)/y0
+//                       d = offset
+//                       e = dz/dx
+// If the offset respectively the slope in z-position is impossible, the parameters are fixed using 
+// results from the simple riemann fit. Afterwards the fit is redone.
+// The curvature is calculated according to the formula:
+//                       curv = a/(1 + b^2 + c*a) = 1/R
+//
+// Paramters:   - Array of tracklets (connected to the track candidate)
+//              - Flag selecting the error definition
+// Output:      - Chi2 value of the track
+//
+// debug level 5
+//
+
+       AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
+       Int_t nTimeBins = cal->GetNumberOfTimeBins();
+
+       TLinearFitter *fitter = GetTiltedRiemanFitter();
+  fitter->StoreData(kTRUE);
+       fitter->ClearPoints();
+       
+       // Calculate the reference position:
+       Int_t nDistances = 0;
+       Float_t meanDistance = 0.;
+       Int_t startIndex = 5;
+       for(Int_t il =5; il > 0; il--){
+               if(tracklets[il].IsOK() && tracklets[il -1].IsOK()){
+                       meanDistance += tracklets[il].GetX0() - tracklets[il -1].GetX0();
+                       nDistances++;
+               }
+               if(tracklets[il].IsOK()) startIndex = il;
+       }
+       meanDistance /= nDistances;
+       if(tracklets[0].IsOK()) startIndex = 0;
+       Double_t xref = tracklets[startIndex].GetX0() + (2.5 - startIndex) * meanDistance - 0.5 * (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
+
+       Float_t x, y, z, t, tilt, xdelta, rhs, error;
+       Float_t dzMean = 0;     Int_t dzcounter = 0;    // A reference z and a reference slope is used if the fitresults in z-direction are not acceptable
+       Double_t uvt[4];
+       Int_t nPoints = 0;
+       for(Int_t ipl = 0; ipl < AliTRDgeometry::kNplan; ipl++){
+               if(!tracklets[ipl].IsOK()) continue;
+               dzMean += tracklets[ipl].GetZfitR(1);
+               dzcounter++;
+               for(Int_t itb = 0; itb < nTimeBins; itb++){
+                       if (!tracklets[ipl].IsUsable(itb)) continue;
+                       x = tracklets[ipl].GetX(itb) + tracklets[ipl].GetX0();
+                       y = tracklets[ipl].GetY(itb);
+                       z = tracklets[ipl].GetZ(itb);
+                       tilt = tracklets[ipl].GetTilt();
+                       xdelta = x - xref;
+                       // Transformation
+                       t = 1/(x*x + y*y);
+                       uvt[0] = 2.0 * x * t;
+                       uvt[1] = t;
+                       uvt[2] = 2.0 * tilt * t;
+                       uvt[3] = 2.0 * tilt * xdelta * t;
+                       rhs = 2.0 * (y + tilt*z) * t;
+                       // error definition changes for the different calls
+                       error = 2.0 * t;
+                       error *= sigError ? tracklets[ipl].GetSigmaY() : 0.2;
+                       fitter->AddPoint(uvt, rhs, error);
+                       nPoints++;
+               }
+       }
+       
+       fitter->Eval();
+
+       Double_t offset = fitter->GetParameter(3);
+       Double_t slope  = fitter->GetParameter(4);
+
+       // Linear fitter  - not possible to make boundaries
+       // Do not accept non possible z and dzdx combinations
+       Bool_t acceptablez = kTRUE;
+       Double_t zref = 0.0;
+       for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNplan; iLayer++) {
+               if(!tracklets[iLayer].IsOK()) continue;
+               zref = offset + slope * (tracklets[iLayer].GetX0() - xref);
+               if (TMath::Abs(tracklets[iLayer].GetZProb() - zref) > tracklets[iLayer].GetPadLength() * 0.5 + 1.0) 
+                       acceptablez = kFALSE;
+       }
+       if (!acceptablez) {
+               dzMean /= dzcounter;
+               Double_t zmf  = tracklets[startIndex].GetZfitR(0) + dzMean * (xref - tracklets[startIndex].GetX0()); // Z-Position of the track at the middle of a stack assuming a linear dependence on x (approximation)
+               fgTiltedRieman->FixParameter(3, zmf);
+               fgTiltedRieman->FixParameter(4, dzMean);
+               fitter->Eval();
+               fitter->ReleaseParameter(3);
+               fitter->ReleaseParameter(4);
+               offset = fitter->GetParameter(3);
+               slope = fitter->GetParameter(4);
+       }
+
+       // Calculate Curvarture
+       Double_t a     =  fitter->GetParameter(0);
+       Double_t b     =  fitter->GetParameter(1);
+       Double_t c     =  fitter->GetParameter(2);
+       Double_t curvature =  1.0 + b*b - c*a;
+       Double_t dca  =  0.0;                                                                                                                   // Distance to closest approach
+       if (curvature > 0.0) {
+               dca = -c / (TMath::Sqrt(1.0 + b*b - c*a) + TMath::Sqrt(1.0 + b*b));
+               curvature  =  a / TMath::Sqrt(curvature);
+       }
+
+       Double_t chi2track = fitter->GetChisquare()/Double_t(nPoints);
+
+       // Update the tracklets
+       Double_t dy, dz;
+       for(Int_t iLayer = 0; iLayer < AliTRDtrackerV1::kNPlanes; iLayer++) {
+
+               x  = tracklets[iLayer].GetX0();
+               y  = 0;
+               z  = 0;
+               dy = 0;
+               dz = 0;
+
+               // y:     R^2 = (x - x0)^2 + (y - y0)^2
+               //     =>   y = y0 +/- Sqrt(R^2 - (x - x0)^2)
+               //          R = Sqrt() = 1/Curvature
+               //     =>   y = y0 +/- Sqrt(1/Curvature^2 - (x - x0)^2)  
+               Double_t res = (x * a + b);                                                             // = (x - x0)/y0
+               res *= res;
+               res  = 1.0 - c * a + b * b - res;                                       // = (R^2 - (x - x0)^2)/y0^2
+               if (res >= 0) {
+                       res = TMath::Sqrt(res);
+                       y    = (1.0 - res) / a;
+               }
+
+               // dy:      R^2 = (x - x0)^2 + (y - y0)^2
+               //     =>     y = +/- Sqrt(R^2 - (x - x0)^2) + y0
+               //     => dy/dx = (x - x0)/Sqrt(R^2 - (x - x0)^2) 
+               // Curvature: cr = 1/R = a/Sqrt(1 + b^2 - c*a)
+               //     => dy/dx =  (x - x0)/(1/(cr^2) - (x - x0)^2) 
+               Double_t x0 = -b / a;
+               if (-c * a + b * b + 1 > 0) {
+                       if (1.0/(curvature * curvature) - (x - x0) * (x - x0) > 0.0) {
+                               Double_t yderiv = (x - x0) / TMath::Sqrt(1.0/(curvature * curvature) - (x - x0) * (x - x0));
+                               if (a < 0) yderiv *= -1.0;
+                               dy = yderiv;
+                       }
+               }
+               z  = offset + slope * (x - xref);
+               dz = slope;
+               tracklets[iLayer].SetYref(0, y);
+               tracklets[iLayer].SetYref(1, dy);
+               tracklets[iLayer].SetZref(0, z);
+               tracklets[iLayer].SetZref(1, dz);
+               tracklets[iLayer].SetC(curvature);
+               tracklets[iLayer].SetChi2(chi2track);
+  }
+
+
+       if(AliTRDReconstructor::StreamLevel() >= 5){
+               Double_t chi2Z = CalculateChi2Z(tracklets, offset, slope);
+               TTreeSRedirector &treeStreamer = *fgDebugStreamer;
+               treeStreamer << "FitTiltedRieman"
+                       << "error="         << sigError
+                       << "Curvature="     << curvature
+                       << "Chi2track="     << chi2track
+                       << "Chi2Z="         << chi2Z
+                       << "D="             << c
+                       << "DCA="           << dca
+                       << "Offset="        << offset
+                       << "Slope="         << slope
+                       << "\n";
+       }
+
+       return chi2track;
+}
+
+//_________________________________________________________________________
+Float_t AliTRDtrackerV1::CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset, Double_t slope)
+{
+//
+// Calculates the chi2-value of the track in z-Direction including tilting pad correction.
+// A linear dependence on the x-value serves as a model.
+// The parameters are related to the tilted Riemann fit.
+// Parameters: - Array of tracklets (AliTRDseedV1) related to the track candidate
+//             - the offset for the reference x
+//             - the slope
+// Output:     - The Chi2 value of the track in z-Direction
+//
+       Double_t xref = .5 * (tracklets[2].GetX0() + tracklets[3].GetX0());
+       Float_t chi2Z = 0, nLayers = 0;
+       for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNplan; iLayer++) {
+               if(!tracklets[iLayer].IsOK()) continue;
+               Double_t z = offset + slope * (tracklets[iLayer].GetX0() - xref);
+               chi2Z += TMath::Abs(tracklets[iLayer].GetMeanz() - z);
+               nLayers++;
+       }
+       chi2Z /= TMath::Max((nLayers - 3.0),1.0);
+       return chi2Z;
+}
+
+
+
 //_____________________________________________________________________________
 Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t maxStep)
 {
@@ -844,7 +1212,7 @@ Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t m
 
     // Calculate the mean material budget between start and
     // end point of this prolongation step
-    AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
+    AliTracker::MeanMaterialBudget(xyz0, xyz1, param);
 
     // Propagate the track to the X-position after the next step
     if (!t.PropagateTo(x,param[1],param[0]*param[4])) {
@@ -863,6 +1231,91 @@ Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t m
 
 }
 
+
+//_____________________________________________________________________________
+Int_t AliTRDtrackerV1::ReadClusters(TClonesArray* &array, TTree *clusterTree) const
+{
+  //
+  // Reads AliTRDclusters from the file. 
+  // The names of the cluster tree and branches 
+  // should match the ones used in AliTRDclusterizer::WriteClusters()
+  //
+
+  Int_t nsize = Int_t(clusterTree->GetTotBytes() / (sizeof(AliTRDcluster))); 
+  TObjArray *clusterArray = new TObjArray(nsize+1000); 
+  
+  TBranch *branch = clusterTree->GetBranch("TRDcluster");
+  if (!branch) {
+    AliError("Can't get the branch !");
+    return 1;
+  }
+  branch->SetAddress(&clusterArray); 
+  
+  if(!fClusters){ 
+       array = new TClonesArray("AliTRDcluster", nsize);
+       array->SetOwner(kTRUE);
+  }
+  
+  // Loop through all entries in the tree
+  Int_t nEntries   = (Int_t) clusterTree->GetEntries();
+  Int_t nbytes     = 0;
+  Int_t ncl        = 0;
+  AliTRDcluster *c = 0x0;
+  for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) {
+    // Import the tree
+    nbytes += clusterTree->GetEvent(iEntry);  
+    
+    // Get the number of points in the detector
+    Int_t nCluster = clusterArray->GetEntriesFast();  
+    for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
+      if(!(c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster))) continue;
+      new((*fClusters)[ncl++]) AliTRDcluster(*c);
+      clusterArray->RemoveAt(iCluster); 
+    }
+
+  }
+  delete clusterArray;
+
+  return 0;
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDtrackerV1::LoadClusters(TTree *cTree)
+{
+  //
+  // Fills clusters into TRD tracking_sectors 
+  // Note that the numbering scheme for the TRD tracking_sectors 
+  // differs from that of TRD sectors
+  //
+
+       
+  if (ReadClusters(fClusters, cTree)) {
+    AliError("Problem with reading the clusters !");
+    return 1;
+  }
+  Int_t ncl  = fClusters->GetEntriesFast(), nin = 0;
+  Int_t icl = ncl;
+  while (icl--) {
+    AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(icl);
+               if(c->IsInChamber()) nin++;
+    Int_t detector       = c->GetDetector();
+    Int_t sector         = fGeom->GetSector(detector);
+    Int_t stack          = fGeom->GetChamber(detector);
+    Int_t plane          = fGeom->GetPlane(detector);
+               
+               fTrSec[sector].GetChamber(stack, plane, kTRUE)->InsertCluster(c, icl);
+  }
+  AliInfo(Form("Clusters %d in %6.2f %%", ncl, 100.*float(nin)/ncl));
+       
+       for(int isector =0; isector<AliTRDgeometry::kNsect; isector++){ 
+               if(!fTrSec[isector].GetNChambers()) continue;
+               fTrSec[isector].Init();
+  }
+  
+  return 0;
+}
+
+
 //____________________________________________________________________
 void AliTRDtrackerV1::UnloadClusters() 
 { 
@@ -870,38 +1323,41 @@ void AliTRDtrackerV1::UnloadClusters()
   // Clears the arrays of clusters and tracks. Resets sectors and timebins 
   //
 
-  Int_t i;
-  Int_t nentr;
+       if(fTracks) fTracks->Delete(); 
+  if(fTracklets) fTracklets->Delete();
+  if(fClusters) fClusters->Delete();
 
-  nentr = fClusters->GetEntriesFast();
-  for (i = 0; i < nentr; i++) {
-    delete fClusters->RemoveAt(i);
-  }
-  fNclusters = 0;
-  
-       if(fTracklets){
-       for (i = 0; i < fTracklets->GetEntriesFast(); i++) delete fTracklets->RemoveAt(i);
-       }
+  for (int i = 0; i < AliTRDgeometry::kNsect; i++) fTrSec[i].Clear();
 
-  nentr = fSeeds->GetEntriesFast();
-  for (i = 0; i < nentr; i++) {
-    delete fSeeds->RemoveAt(i);
-  }
+}
 
-  nentr = fTracks->GetEntriesFast();
-       for (i = 0; i < nentr; i++) {
-    delete fTracks->RemoveAt(i);
-  }
+//_____________________________________________________________________________
+Bool_t AliTRDtrackerV1::AdjustSector(AliTRDtrackV1 *track) 
+{
+  //
+  // Rotates the track when necessary
+  //
+
+  Double_t alpha = AliTRDgeometry::GetAlpha(); 
+  Double_t y     = track->GetY();
+  Double_t ymax  = track->GetX()*TMath::Tan(0.5*alpha);
 
-  Int_t nsec = AliTRDgeometry::kNsect;
-  for (i = 0; i < nsec; i++) {    
-    for(Int_t pl = 0; pl < fTrSec[i]->GetNumberOfLayers(); pl++) {
-      fTrSec[i]->GetLayer(pl)->Clear();
+  if      (y >  ymax) {
+    if (!track->Rotate( alpha)) {
+      return kFALSE;
     }
-  }
+  } 
+  else if (y < -ymax) {
+    if (!track->Rotate(-alpha)) {
+      return kFALSE;   
+    }
+  } 
+
+  return kTRUE;
 
 }
 
+
 //____________________________________________________________________
 AliTRDseedV1* AliTRDtrackerV1::GetTracklet(AliTRDtrackV1 *track, Int_t p, Int_t &idx)
 {
@@ -917,23 +1373,7 @@ AliTRDseedV1* AliTRDtrackerV1::GetTracklet(AliTRDtrackV1 *track, Int_t p, Int_t
 // Detailed description
 //
        idx = track->GetTrackletIndex(p);
-       //AliInfo(Form("looking for tracklet in plane %d idx %d [%d]", p, idx, track->GetTrackletIndex(p)));
        AliTRDseedV1 *tracklet = idx<0 ? 0x0 : (AliTRDseedV1*)fTracklets->UncheckedAt(idx);
-       //AliInfo(Form("found 0x%x @ %d", tracklet, idx));
-
-//   Int_t *index = track->GetTrackletIndexes();
-//   for (UInt_t i = 0; i < 6; i++) AliInfo(Form("index[%d] = %d", i, index[i]));
-// 
-//     for (UInt_t i = 0; i < 6/*kMaxTimeBinIndex*/; i++) {
-//     if (index[i] < 0) continue;
-// 
-//             tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index[i]);
-//             if(!tracklet) break;
-// 
-//             if(tracklet->GetPlane() != p) continue;
-// 
-//             idx = index[i];
-//     }
 
        return tracklet;
 }
@@ -963,8 +1403,7 @@ Int_t AliTRDtrackerV1::SetTracklet(AliTRDseedV1 *tracklet)
 }
 
 //____________________________________________________________________
-Int_t AliTRDtrackerV1::Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *sector
-                                       , AliESDEvent *esd)
+Int_t AliTRDtrackerV1::Clusters2TracksSM(Int_t sector, AliESDEvent *esd)
 {
   //
   // Steer tracking for one SM.
@@ -985,57 +1424,36 @@ Int_t AliTRDtrackerV1::Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *se
   // 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 ksmLayer(*(sector->GetLayer(ilayer)));
-                       stackLayer[ilayer] = ksmLayer;
-#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();
+       Int_t nTracks   = 0;
+       Int_t nChambers = 0;
+       AliTRDtrackingChamber **stack = 0x0, *chamber = 0x0;
+       for(int istack = 0; istack<AliTRDgeometry::kNcham; istack++){
+               if(!(stack = fTrSec[sector].GetStack(istack))) continue;
+               nChambers = 0;
+               for(int iplane=0; iplane<AliTRDgeometry::kNplan; iplane++){
+                       if(!(chamber = stack[iplane])) continue;
+                       if(chamber->GetNClusters() < fTimeBinsPerPlane * AliTRDReconstructor::RecoParam()->GetFindableClusters()) continue;
+                       nChambers++;
+                       //AliInfo(Form("sector %d stack %d plane %d clusters %d", sector, istack, iplane, chamber->GetNClusters()));
                }
-               //AliInfo(Form("Finish building stack propagation layers. nClusters %d.", nClStack));
-               if(nClStack < kFindable) continue;
-               ntracks += Clusters2TracksStack(&stackLayer[0], &esdTrackList);
+               if(nChambers < 4) continue;
+               //AliInfo(Form("Doing stack %d", istack));
+               nTracks += Clusters2TracksStack(stack, &esdTrackList);
        }
-       //AliInfo(Form("Found %d tracks in SM", ntracks));
+       //AliInfo(Form("Found %d tracks in SM %d [%d]\n", nTracks, sector, esd->GetNumberOfTracks()));
        
-       for(int itrack=0; itrack<ntracks; itrack++) 
+       for(int itrack=0; itrack<nTracks; itrack++) 
           esd->AddTrack((AliESDtrack*)esdTrackList[itrack]);
 
-       return ntracks;
+       return nTracks;
 }
 
 //____________________________________________________________________
-Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
-                                          , TClonesArray *esdTrackList)
+Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClonesArray *esdTrackList)
 {
   //
   // Make tracks in one TRD stack.
@@ -1063,6 +1481,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
   // 8. Build ESD track and register it to the output list
   //
 
+       AliTRDtrackingChamber *chamber = 0x0;
        AliTRDseedV1 sseed[kMaxTracksStack*6]; // to be initialized
        Int_t pars[4]; // MakeSeeds parameters
 
@@ -1071,13 +1490,12 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
        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
-
+       Double_t quality = BuildSeedingConfigs(stack, configs);
+       if(AliTRDReconstructor::StreamLevel() > 1){
+       AliInfo(Form("Plane config %d %d %d Quality %f"
+    , configs[0], configs[1], configs[2], quality));
+       }
+       
        // Initialize contors
        Int_t ntracks,      // number of TRD track candidates
              ntracks1,     // number of registered TRD tracks/iter
@@ -1088,14 +1506,12 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                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);
+                       pars[1] = ntracks;
+                       ntracks = MakeSeeds(stack, &sseed[6*ntracks], pars);
                        if(ntracks == kMaxTracksStack) break;
                }
-#ifdef DEBUG
-               if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Candidate TRD tracks %d in stack %d iteration %d.", ntracks, pars[1], fSieveSeeding));
-#endif         
+               if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Candidate TRD tracks %d in iteration %d.", ntracks, fSieveSeeding));
+               
                if(!ntracks) break;
                
                // Sort the seeds according to their quality
@@ -1136,10 +1552,9 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                                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;
+                                       if (TMath::Abs(sseed[jseed].GetYref(0) / sseed[jseed].GetX0()) < 0.15) findable++;
+       
                                        sseed[jseed].UpdateUsed();
                                        ncl   += sseed[jseed].GetN2();
                                        nused += sseed[jseed].GetNUsed();
@@ -1158,12 +1573,12 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                                }
                                // Filter duplicated tracks
                                if (nused > 30){
-                                       printf("Skip %d nused %d\n", trackIndex, nused);
+                                       //printf("Skip %d nused %d\n", trackIndex, nused);
                                        fakeTrack[trackIndex] = kTRUE;
                                        continue;
                                }
                                if (Float_t(nused)/ncl >= .25){
-                                       printf("Skip %d nused/ncl >= .25\n", trackIndex);
+                                       //printf("Skip %d nused/ncl >= .25\n", trackIndex);
                                        fakeTrack[trackIndex] = kTRUE;
                                        continue;
                                }
@@ -1197,7 +1612,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                                }
                                if(skip){
                                        candidates++;
-                                       printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused);
+                                       //printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused);
                                        continue;
                                }
                                signedTrack[trackIndex] = kTRUE;
@@ -1255,24 +1670,22 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                                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*/
+                               Int_t ich = 0; while(!(chamber = stack[ich])) ich++;
+                               trackParams[6] = fGeom->GetSector(chamber->GetDetector());/* *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();
+                               if(AliTRDReconstructor::StreamLevel() > 1){
+                                       AliInfo(Form("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1]));
+                                       
                                        Int_t nclusters = 0;
                                        AliTRDseedV1 *dseed[6];
                                        for(int is=0; is<6; is++){
                                                dseed[is] = new AliTRDseedV1(sseed[trackIndex*6+is]);
                                                dseed[is]->SetOwner();
                                                nclusters += sseed[is].GetN2();
-                                               //for(int ic=0; ic<30; ic++) if(sseed[trackIndex*6+is].GetClusters(ic)) printf("l[%d] tb[%d] cptr[%p]\n", is, ic, sseed[trackIndex*6+is].GetClusters(ic));
                                        }
                                        //Int_t eventNrInFile = esd->GetEventNumberInFile();
                                        //AliInfo(Form("Number of clusters %d.", nclusters));
-                                       TTreeSRedirector &cstreamer = *fDebugStreamer;
+                                       TTreeSRedirector &cstreamer = *fgDebugStreamer;
                                        cstreamer << "Clusters2TracksStack"
                                                << "Iter="      << fSieveSeeding
                                                << "Like="      << fTrackQuality[trackIndex]
@@ -1289,8 +1702,6 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                                                << "p4=" << trackParams[4]
                                                << "p5=" << trackParams[5]
                                                << "p6=" << trackParams[6]
-                                               << "Sector="    << sector
-                                               << "Stack="     << pars[1]
                                                << "Label="     << label
                                                << "Label1="    << label1
                                                << "Label2="    << label2
@@ -1301,16 +1712,14 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                                                << "Findable="  << findable
                                                << "NUsed="     << nused
                                                << "\n";
-                                       //???for(int is=0; is<6; is++) delete dseed[is];
                                }
-#endif
                        
-                               AliTRDtrackV1 *track = AliTRDtrackerV1::MakeTrack(&sseed[trackIndex*kNPlanes], trackParams);
+                               AliTRDtrackV1 *track = MakeTrack(&sseed[trackIndex*kNPlanes], trackParams);
                                if(!track){
-                                       AliWarning("Fail to build a TRD Track.");
+                                       //AliWarning("Fail to build a TRD Track.");
                                        continue;
                                }
-                               AliInfo("End of MakeTrack()");
+                               //AliInfo("End of MakeTrack()");
                                AliESDtrack esdTrack;
                                esdTrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
                                esdTrack.SetLabel(track->GetLabel());
@@ -1327,14 +1736,17 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
                fSieveSeeding++;
 
                // Rebuild plane configurations and indices taking only unused clusters into account
-               quality = BuildSeedingConfigs(layer, configs);
-               //if(quality < fRecoParam->GetPlaneQualityThreshold()) break;
+               quality = BuildSeedingConfigs(stack, configs);
+               if(quality < 1.E-7) break; //AliTRDReconstructor::RecoParam()->GetPlaneQualityThreshold()) break;
                
-               for(Int_t il = 0; il < kNPlanes * fTimeBinsPerPlane; il++) layer[il].BuildIndices(fSieveSeeding);
+               for(Int_t ip = 0; ip < kNPlanes; ip++){ 
+                       if(!(chamber = stack[ip])) continue;
+                       chamber->Build(fGeom);//Indices(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
+               if(AliTRDReconstructor::StreamLevel() > 1){ 
+                       AliInfo(Form("Sieve level %d Plane config %d %d %d Quality %f", fSieveSeeding, configs[0], configs[1], configs[2], quality));
+               }
        } while(fSieveSeeding<10); // end stack clusters sieve
        
 
@@ -1345,8 +1757,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer
 }
 
 //___________________________________________________________________
-Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDstackLayer *layers
-                                            , Int_t *configs)
+Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDtrackingChamber **stack, Int_t *configs)
 {
   //
   // Assign probabilities to chambers according to their
@@ -1372,10 +1783,11 @@ Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDstackLayer *layers
   // The overall chamber quality is given by the product of this 2 contributions.
   // 
 
-       Double_t chamberQA[kNPlanes];
+       Double_t chamberQ[kNPlanes];
+       AliTRDtrackingChamber *chamber = 0x0;
        for(int iplane=0; iplane<kNPlanes; iplane++){
-               chamberQA[iplane] = MakeSeedingPlanes(&layers[iplane*fTimeBinsPerPlane]);
-               //printf("chamberQA[%d] = %f\n", iplane, chamberQA[iplane]);
+               if(!(chamber = stack[iplane])) continue;
+               chamberQ[iplane] = (chamber = stack[iplane]) ?  chamber->GetQuality(fTimeBinsPerPlane) : 0.;
        }
 
        Double_t tconfig[kNConfigs];
@@ -1383,17 +1795,19 @@ Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDstackLayer *layers
        for(int iconf=0; iconf<kNConfigs; iconf++){
                GetSeedingConfig(iconf, planes);
                tconfig[iconf] = fgTopologicQA[iconf];
-               for(int iplane=0; iplane<4; iplane++) tconfig[iconf] *= chamberQA[planes[iplane]]; 
+               for(int iplane=0; iplane<4; iplane++) tconfig[iconf] *= chamberQ[planes[iplane]]; 
        }
        
        TMath::Sort(kNConfigs, tconfig, configs, kTRUE);
+//     AliInfo(Form("q[%d] = %f", configs[0], tconfig[configs[0]]));
+//     AliInfo(Form("q[%d] = %f", configs[1], tconfig[configs[1]]));
+//     AliInfo(Form("q[%d] = %f", configs[2], tconfig[configs[2]]));
+       
        return tconfig[configs[0]];
 }
 
 //____________________________________________________________________
-Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
-                               , AliTRDseedV1 *sseed
-                               , Int_t *ipar)
+Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *sseed, Int_t *ipar)
 {
   //
   // Make tracklet seeds in the TRD stack.
@@ -1416,7 +1830,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
   // 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().
+  //   layer 0, 1, and 2 are defined in AliTRDchamberTimeBin::BuildCond().
   // 3. Helix fit of the cluster set. (see AliTRDtrackerFitter::FitRieman(AliTRDcluster**))
   // 4. Initialize seeding tracklets in the seeding chambers.
   // 5. Filter 0.
@@ -1444,10 +1858,11 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
   // 15. Register seeds.
   //
 
+       AliTRDtrackingChamber *chamber = 0x0;
        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];
+       Int_t index[AliTRDchamberTimeBin::kMaxClustersLayer], jndex[AliTRDchamberTimeBin::kMaxClustersLayer];
        // chi2 storage
        // chi2[0] = tracklet chi2 on the Z direction
        // chi2[1] = tracklet chi2 on the R direction
@@ -1459,39 +1874,44 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
        
        // unpack control parameters
        Int_t config  = ipar[0];
-       Int_t istack  = ipar[1];
-       Int_t ntracks = ipar[2];
+       Int_t ntracks = ipar[1];
        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
-       Int_t det/*, tbRange[6]*/; // time bins inside the detector geometry
+       Int_t ic = 0; while(!(chamber = stack[ic])) ic++;
+       Int_t istack = fGeom->GetChamber(chamber->GetDetector());
        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] = pp->GetLengthIPad();
-               det           = il; // to be fixed !!!!!
-               //tbRange[il]   = fTimeBinsPerPlane; //Int_t(AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght() * AliTRDCommonParam::Instance()->GetSamplingFrequency()/AliTRDcalibDB::Instance()->GetVdriftDet()->GetValue(det));
-               //printf("%d hl[%f] pl[%f] tb[%d]\n", il, hL[il], padlength[il], tbRange[il]);
+       AliTRDpadPlane *pp = 0x0;
+       for(int iplane=0; iplane<kNPlanes; iplane++){
+               pp                = fGeom->GetPadPlane(iplane, istack);
+               hL[iplane]        = TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle());
+               padlength[iplane] = pp->GetLengthIPad();
+       }
+       
+       if(AliTRDReconstructor::StreamLevel() > 1){
+               AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks));
        }
 
-       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] * fTimeBinsPerPlane], planes[isl]);
-
-
+       Int_t nlayers = 0;
+       AliTRDchamberTimeBin *layer[] = {0x0, 0x0, 0x0, 0x0};
+       for(int isl=0; isl<kNSeedPlanes; isl++){ 
+               if(!(chamber = stack[planes[isl]])) continue;
+               if(!(layer[isl] = chamber->GetSeedingLayer(fGeom))) continue;
+               nlayers++;
+               //AliInfo(Form("seeding plane %d clusters %d", planes[isl], Int_t(*layer[isl])));
+       }
+       if(nlayers < 4) return 0;
+       
+       
        // Start finding seeds
+       Double_t cond0[4], cond1[4], cond2[4];
        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);
+               //printf("Found c[3] candidates 0 %d\n", ncl);
                Int_t jcl = 0;
                while(jcl<ncl) {
                        c[0] = (*layer[0])[index[jcl++]];
@@ -1501,6 +1921,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                        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);
+                       //printf("Found c[0] candidates 1 %d\n", mcl);
 
                        Int_t kcl = 0;
                        while(kcl<mcl) {
@@ -1508,78 +1929,63 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                if(!c[1]) continue;
                                layer[2]->BuildCond(c[1], cond2, 2, theta, phi);
                                c[2] = layer[2]->GetNearestCluster(cond2);
+                               //printf("Found c[1] candidate 2 %p\n", c[2]);
                                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());
+//                             AliInfo("Seeding clusters found. Building seeds ...");
+//                             for(Int_t i = 0; i < kNSeedPlanes; i++) printf("%i. coordinates: x = %6.3f, y = %6.3f, z = %6.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);
+                               FitRieman(c, chi2);
 
-                               chi2[0] = 0.; chi2[1] = 0.;
                                AliTRDseedV1 *tseed = 0x0;
                                for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
                                        Int_t jLayer = planes[iLayer];
                                        tseed = &cseed[jLayer];
-                                       tseed->SetRecoParam(fRecoParam);
                                        tseed->SetPlane(jLayer);
                                        tseed->SetTilt(hL[jLayer]);
                                        tseed->SetPadLength(padlength[jLayer]);
-                                       //tseed->SetNTimeBinsRange(tbRange[jLayer]);
-                                       tseed->SetX0(layer[iLayer]->GetX());//layers[jLayer*fTimeBinsPerPlane].GetX());
-
-                                       tseed->Init(fFitter->GetRiemanFitter());
-                                       // temporary until new AttachClusters()
-                                       tseed->SetX0(layers[(jLayer+1)*fTimeBinsPerPlane-1].GetX());
-                                       chi2[0] += tseed->GetChi2Z(c[iLayer]->GetZ());
-                                       chi2[1] += tseed->GetChi2Y(c[iLayer]->GetY());
+                                       tseed->SetX0(stack[jLayer]->GetX());
+                                       tseed->Init(GetRiemanFitter());
                                }
 
                                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);
+                                       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;
+                                       
+                                       Float_t yref[4];
+                                       for(int il=0; il<4; il++) yref[il] = cseed[planes[il]].GetYref(0);
                                        Int_t ll = c[3]->GetLabel(0);
-                                       TTreeSRedirector &cs0 = *fDebugStreamer;
+                                       TTreeSRedirector &cs0 = *fgDebugStreamer;
                                                        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]
+                                                       <<"chi2z=" << chi2[0]
+                                                       <<"chi2y=" << chi2[1]
+                                                       <<"yref0=" << yref[0]
+                                                       <<"yref1=" << yref[1]
+                                                       <<"yref2=" << yref[2]
+                                                       <<"yref3=" << yref[3]
+                                                       <<"c0.="   << c[0]
+                                                       <<"c1.="   << c[1]
+                                                       <<"c2.="   << c[2]
+                                                       <<"c3.="   << c[3]
                                                        <<"\n";
                                }
-#endif
 
-                               if(chi2[0] > fRecoParam->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){
+                               if(chi2[0] > AliTRDReconstructor::RecoParam()->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*/){
+                               if(chi2[1] > AliTRDReconstructor::RecoParam()->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++) {
@@ -1590,7 +1996,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                        }
                                        Double_t xpos[4];
                                        for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX();
-                                       TTreeSRedirector &cstreamer = *fDebugStreamer;
+                                       TTreeSRedirector &cstreamer = *fgDebugStreamer;
                                                        cstreamer << "MakeSeeds1"
                                                << "isFake=" << isFake
                                                << "config="   << config
@@ -1614,33 +2020,25 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                                << "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++){
                                        Int_t jLayer = planes[iLayer];
-                                       if(!cseed[jLayer].AttachClustersIter(&layers[jLayer*fTimeBinsPerPlane], 5., kFALSE, c[iLayer])) continue;
+                                       if(!cseed[jLayer].AttachClustersIter(stack[jLayer], 5., kFALSE, c[iLayer])) continue;
                                        nUsedCl += cseed[jLayer].GetNUsed();
                                        if(nUsedCl > 25) break;
                                        nlayers++;
                                }
                                if(nlayers < kNSeedPlanes){ 
-                                       //AliInfo("Failed updating all seeds.");
+                                       //AliInfo(Form("Failed updating all seeds %d [%d].", nlayers, kNSeedPlanes));
                                        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]].Init(rim);
-                                       chi2[0] += (Float_t)cseed[planes[iLayer]].GetChi2Z();
-                                       chi2[1] += cseed[planes[iLayer]].GetChi2Y();
-                               }
-                               Double_t chi2r = chi2[1], chi2z = chi2[0];
+                               FitRieman(&cseed[0], chi2, &planes[0]);
                                Double_t like = CookLikelihood(&cseed[0], planes, chi2); // to be checked
-                               if (TMath::Log(1.E-9 + like) < fRecoParam->GetTrackLikelihood()){
+                               if (TMath::Log(1.E-9 + like) < AliTRDReconstructor::RecoParam()->GetTrackLikelihood()){
                                        //AliInfo(Form("Failed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
                                        continue;
                                }
@@ -1657,72 +2055,41 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                Int_t nusedf   = 0; // debug value
                                for(int iLayer=0; iLayer<2; iLayer++){
                                        Int_t jLayer = lextrap[iLayer];
-                                       
+                                       if(!(chamber = stack[jLayer])) continue;
+                                               
                                        // prepare extrapolated seed
                                        cseed[jLayer].Reset();
-                                       cseed[jLayer].SetRecoParam(fRecoParam);
                                        cseed[jLayer].SetPlane(jLayer);
                                        cseed[jLayer].SetTilt(hL[jLayer]);
-                                       cseed[jLayer].SetX0(layers[(jLayer +1) * fTimeBinsPerPlane-1].GetX());
+                                       cseed[jLayer].SetX0(chamber->GetX());
                                        cseed[jLayer].SetPadLength(padlength[jLayer]);
-                                       //cseed[jLayer].SetNTimeBinsRange(tbRange[jLayer]);
-                                       cseed[jLayer].Init(rim);
-//                                     AliTRDcluster *cd = FindSeedingCluster(&layers[jLayer*fTimeBinsPerPlane], &cseed[jLayer]);
-//                                     if(cd == 0x0) continue;
 
                                        // fit extrapolated seed
-                                       AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
+                                       FitTiltedRieman(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*fTimeBinsPerPlane], 1000.)) continue;
+                                       if(!tseed.AttachClustersIter(chamber, 1000.)) continue;
                                        cseed[jLayer] = tseed;
                                        nusedf += cseed[jLayer].GetNUsed(); // debug value
-                                       AliTRDseedV1::FitRiemanTilt(cseed, kTRUE);
                                }
+                               FitTiltedRieman(cseed, kTRUE);
                                //AliInfo("Extrapolation done.");
 
-                               ImproveSeedQuality(layers, cseed);
+                               if(ImproveSeedQuality(stack, cseed) < 4) continue;
                                //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].Init(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);
+                               Double_t curv = FitRieman(&cseed[0], chi2);
+                               Double_t chi2ZF = chi2[0] / TMath::Max((nlayers - 3.), 1.);
+                               Double_t chi2RF = chi2[1] / TMath::Max((nlayers - 3.), 1.);
+
+                               // do the final track fitting (Once with vertex constraint and once without vertex constraint)
+                               Double_t chi2Vals[3];
+                               chi2Vals[0] = FitTiltedRieman(&cseed[0], kFALSE);
+                               chi2Vals[1] = FitTiltedRiemanConstraint(&cseed[0], GetZ());
+                               chi2Vals[2] = chi2ZF;
+                               fTrackQuality[ntracks] = CalculateTrackLikelihood(&cseed[0], &chi2Vals[0]);
                                //AliInfo("Hyperplane fit done\n");
 
                                // finalize tracklets
@@ -1747,29 +2114,20 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                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;
+                                       TTreeSRedirector &cstreamer = *fgDebugStreamer;
                                        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]
@@ -1781,11 +2139,11 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDstackLayer *layers
                                                << "Freq="    << frequency
                                                << "\n";
                                }
-#endif
                                
                                ntracks++;
                                if(ntracks == kMaxTracksStack){
                                        AliWarning(Form("Number of seeds reached maximum allowed (%d) in stack.", kMaxTracksStack));
+                                       for(int isl=0; isl<4; isl++) delete layer[isl];
                                        return ntracks;
                                }
                                cseed += 6;
@@ -1829,7 +2187,7 @@ AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params)
        track->PropagateTo(params[0]-5.0);
   track->ResetCovariance(1);
   Int_t nc = FollowBackProlongation(*track);
-       AliInfo(Form("N clusters for track %d", nc));
+       //AliInfo(Form("N clusters for track %d", nc));
        if (nc < 30) {
     delete track;
     track = 0x0;
@@ -1849,8 +2207,7 @@ void AliTRDtrackerV1::CookLabel(AliKalmanTrack */*pt*/, Float_t /*wrong*/) const
 }
 
 //____________________________________________________________________
-void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers
-                                       , AliTRDseedV1 *cseed)
+Int_t AliTRDtrackerV1::ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDseedV1 *cseed)
 {
   //
   // Sort tracklets according to "quality" and try to "improve" the first 4 worst
@@ -1869,14 +2226,12 @@ void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers
   // 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
+       AliTRDtrackingChamber *chamber = 0x0;
        AliTRDseedV1 bseed[6];
+       Int_t nLayers = 0;
        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;
@@ -1888,88 +2243,80 @@ void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers
 
                for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
                        squality[jLayer]  = bseed[jLayer].IsOK() ? bseed[jLayer].GetQuality(kTRUE) : -1.;
-                       sumquality +=squality[jLayer];
+                       sumquality += squality[jLayer];
                }
                if ((sumquality >= lastquality) || (chi2       >     lastchi2)) break;
 
-               
+               nLayers = 0;
                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);
+                       if(!(chamber = stack[bLayer])) continue;
+                       bseed[bLayer].AttachClustersIter(chamber, squality[bLayer], kTRUE);
+                       if(bseed[bLayer].IsOK()) nLayers++;
                }
 
-               chi2 = AliTRDseedV1::FitRiemanTilt(bseed,kTRUE);
+               chi2 = FitTiltedRieman(bseed, kTRUE);
        } // Loop: iter
+       
+       // we are sure that at least 2 tracklets are OK !
+       return nLayers+2;
 }
 
-//____________________________________________________________________
-Double_t  AliTRDtrackerV1::MakeSeedingPlanes(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 AliTRDtrackerV1::CalculateTrackLikelihood(AliTRDseedV1 *tracklets, Double_t *chi2){
+//
+// Calculates the Track Likelihood value. This parameter serves as main quality criterion for 
+// the track selection
+// The likelihood value containes:
+//    - The chi2 values from the both fitters and the chi2 values in z-direction from a linear fit
+//    - The Sum of the Parameter  |slope_ref - slope_fit|/Sigma of the tracklets
+// For all Parameters an exponential dependency is used
+//
+// Parameters: - Array of tracklets (AliTRDseedV1) related to the track candidate
+//             - Array of chi2 values: 
+//                 * Non-Constrained Tilted Riemann fit
+//                 * Vertex-Constrained Tilted Riemann fit
+//                 * z-Direction from Linear fit
+// Output:     - The calculated track likelihood
+//
+// debug level 2
+//
 
-//     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++;
+       Double_t sumdaf = 0, nLayers = 0;
+       for (Int_t iLayer = 0; iLayer < kNPlanes; iLayer++) {
+               if(!tracklets[iLayer].IsOK()) continue;
+               sumdaf += TMath::Abs((tracklets[iLayer].GetYfit(1) - tracklets[iLayer].GetYref(1))/ tracklets[iLayer].GetSigmaY2());
+               nLayers++;
        }
+       sumdaf /= Float_t (nLayers - 2.0);
        
-       // calculate the deviation of the mean number of clusters from the
-       // closest integer values
-       Float_t nclMed = float(ncl-nused)/nTimeBins;
-       Int_t ncli = Int_t(nclMed);
-       Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
-       nclDev -= (nclDev>.5) && ncli ? 0. : 1.; 
-       /*Double_t quality = */ return TMath::Exp(2.*nclDev);
-
-//     // 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 likeChi2Z  = TMath::Exp(-chi2[2] * 0.14);                      // Chi2Z 
+       Double_t likeChi2TC = TMath::Exp(-chi2[1] * 0.677);                     // Constrained Tilted Riemann
+       Double_t likeChi2TR = TMath::Exp(-chi2[0] * 0.78);                      // Non-constrained Tilted Riemann
+       Double_t likeAF     = TMath::Exp(-sumdaf * 3.23);
+       Double_t trackLikelihood     = likeChi2Z * likeChi2TR * likeAF;
+
+       if(AliTRDReconstructor::StreamLevel() >= 2){
+               TTreeSRedirector &cstreamer = *fgDebugStreamer;
+               cstreamer << "CalculateTrackLikelihood0"
+                       << "LikeChi2Z="         << likeChi2Z
+                       << "LikeChi2TR="        << likeChi2TR
+                       << "LikeChi2TC="        << likeChi2TC
+                       << "LikeAF="                    << likeAF
+                       << "TrackLikelihood=" << trackLikelihood
+                       << "\n";
+       }
+
+       return trackLikelihood;
 }
 
 //____________________________________________________________________
-Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed
-                                       , Int_t planes[4]
+Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4]
                                        , Double_t *chi2)
 {
   //
@@ -1999,7 +2346,7 @@ Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed
        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();
+       Float_t fgFindable = AliTRDReconstructor::RecoParam()->GetFindableClusters();
 
        
        Int_t nclusters = 0;
@@ -2018,10 +2365,9 @@ Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed
        
        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;
+               TTreeSRedirector &cstreamer = *fgDebugStreamer;
                cstreamer << "CookLikelihood"
                        << "sumda="     << sumda
                        << "chi0="      << chi2[0]
@@ -2034,22 +2380,20 @@ Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed
                        << "like="      << like
                        << "\n";
        }
-#endif
 
        return like;
 }
 
+
 //___________________________________________________________________
-void AliTRDtrackerV1::GetMeanCLStack(AliTRDstackLayer *layers
-                                   , Int_t *planes
-                                   , Double_t *params)
+void AliTRDtrackerV1::GetMeanCLStack(AliTRDtrackingChamber *chamber, Int_t *planes, Double_t *params)
 {
   //
   // Determines the Mean number of clusters per layer.
   // Needed to determine good Seeding Layers
   //
   // Parameters:
-  //    - Array of AliTRDstackLayers
+  //    - Array of AliTRDchamberTimeBins
   //    - Container for the params
   //
   // Detailed description
@@ -2061,23 +2405,21 @@ void AliTRDtrackerV1::GetMeanCLStack(AliTRDstackLayer *layers
   // 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;
+       AliTRDchamberTimeBin *layers = chamber->GetTB(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();
+               for(Int_t ils = 0; ils < fTimeBinsPerPlane; ils++){
+                       position = planes[ipl]*fTimeBinsPerPlane + ils;
+                       ncl[ipl * fTimeBinsPerPlane + ils] = layers[position].GetNClusters();
                        nused = 0;
-                       for(Int_t icl = 0; icl < ncl[ipl * nTimeBins + ils]; icl++)
+                       for(Int_t icl = 0; icl < ncl[ipl * fTimeBinsPerPlane + ils]; icl++)
                                if((layers[position].GetCluster(icl))->IsUsed()) nused++;
-                       ncl[ipl * nTimeBins + ils] -= nused;
+                       ncl[ipl * fTimeBinsPerPlane + ils] -= nused;
                }
        }
        // Declaration of quartils:
@@ -2087,13 +2429,13 @@ void AliTRDtrackerV1::GetMeanCLStack(AliTRDstackLayer *layers
        Int_t counter;
        Double_t *array;
        Int_t *limit;
-       Int_t nLayers = nTimeBins * kNSeedPlanes;
+       Int_t nLayers = fTimeBinsPerPlane * 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++){
+                       for(Int_t i = 0; i < fTimeBinsPerPlane *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
@@ -2117,8 +2459,7 @@ void AliTRDtrackerV1::GetMeanCLStack(AliTRDstackLayer *layers
 }
 
 //___________________________________________________________________
-Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers
-                                      , Double_t *params)
+Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDtrackingChamber *chamber, Double_t *params)
 {
   //
   // Algorithm to find optimal seeding layer
@@ -2127,7 +2468,7 @@ Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers
   // All layers outside are sorted according t
   //
   // Parameters:
-  //     - Array of AliTRDstackLayers (in the current plane !!!)
+  //     - Array of AliTRDchamberTimeBins (in the current plane !!!)
   //     - Container for the Indices of the seeding Layer candidates
   //
   // Output:
@@ -2138,15 +2479,15 @@ Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers
   //
 
        //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();
+       const Int_t kMaxClustersLayer = AliTRDchamberTimeBin::kMaxClustersLayer;
        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);
+       
+       AliTRDchamberTimeBin *layers = chamber->GetTB(0);
        Int_t nused = 0;
-       for(Int_t ils = 0; ils < nTimeBins; ils++){
+       for(Int_t ils = 0; ils < fTimeBinsPerPlane; ils++){
                ncl[ils] = layers[ils].GetNClusters();
                nused = 0;
                for(Int_t icl = 0; icl < ncl[ils]; icl++)
@@ -2155,8 +2496,8 @@ Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers
        }
        
        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));
+       for(Int_t ils = 0; ils < fTimeBinsPerPlane; ils++){
+               memmove(indices + bins[ncl[ils]+1] + 1, indices + bins[ncl[ils]+1], sizeof(Int_t)*(fTimeBinsPerPlane - ils));
                indices[bins[ncl[ils]+1]] = ils;
                for(Int_t i = ncl[ils]+1; i < kMaxClustersLayer; i++)
                        bins[i]++;
@@ -2178,7 +2519,7 @@ Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers
                nElements = bins[sbin + 1] - bins[sbin];
                printf("nElements = %d\n", nElements);
                if(iter == 5){
-                       position = (Int_t)(gRandom->Rndm()*(nTimeBins-1));
+                       position = (Int_t)(gRandom->Rndm()*(fTimeBinsPerPlane-1));
                        break;
                }
                else if(nElements==0){
@@ -2193,8 +2534,7 @@ Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers
 }
 
 //____________________________________________________________________
-AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDstackLayer *layers
-                                                 , AliTRDseedV1/*AliRieman*/ *reference)
+AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDtrackingChamber *chamber, AliTRDseedV1* reference) const
 {
   //
   // Finds a seeding Cluster for the extrapolation chamber.
@@ -2206,14 +2546,12 @@ AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDstackLayer *layers
   // 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!!!)
+  // Imput parameters: - layers: array of AliTRDchamberTimeBins (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;
@@ -2221,11 +2559,12 @@ AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDstackLayer *layers
        ypos = reference->GetYref(0);
        zpos = reference->GetZref(0);
        AliTRDcluster *currentBest = 0x0, *temp = 0x0;
-       for(Int_t ils = 0; ils < nTimeBins; ils++){
+       AliTRDchamberTimeBin *layers = chamber->GetTB(0);
+       for(Int_t ils = 0; ils < fTimeBinsPerPlane; 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());
+               index = layers[ils].SearchNearestCluster(ypos, zpos, AliTRDReconstructor::RecoParam()->GetRoad2y(), AliTRDReconstructor::RecoParam()->GetRoad2z());
                if(index == -1) continue;
                temp = layers[ils].GetCluster(index);
                if(!temp) continue;
@@ -2238,192 +2577,6 @@ AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDstackLayer *layers
        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){
-               TMatrixD hist(nRows, nCols);
-               for(Int_t i = 0; i < nRows; i++)
-                       for(Int_t j = 0; j < nCols; j++)
-                               hist(i,j) = histogram[i][j];
-               TTreeSRedirector &cstreamer = *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])
@@ -2645,3 +2798,70 @@ void AliTRDtrackerV1::GetExtrapolationConfig(Int_t iconfig, Int_t planes[2])
                break;
        }
 }
+
+//____________________________________________________________________
+AliCluster* AliTRDtrackerV1::GetCluster(Int_t idx) const
+{
+       Int_t ncls = fClusters->GetEntriesFast();
+       return idx >= 0 || idx < ncls ? (AliCluster*)fClusters->UncheckedAt(idx) : 0x0;
+}
+
+
+//_____________________________________________________________________________
+Int_t AliTRDtrackerV1::Freq(Int_t n, const Int_t *inlist
+                        , Int_t *outlist, Bool_t down)
+{    
+  //
+  // Sort eleements according occurancy 
+  // The size of output array has is 2*n 
+  //
+
+  if (n <= 0) {
+    return 0;
+  }
+
+  Int_t *sindexS = new Int_t[n];   // Temporary array for sorting
+  Int_t *sindexF = new Int_t[2*n];   
+  for (Int_t i = 0; i < n; i++) {
+    sindexF[i] = 0;
+  }
+
+  TMath::Sort(n,inlist,sindexS,down); 
+  Int_t last     = inlist[sindexS[0]];
+  Int_t val      = last;
+  sindexF[0]     = 1;
+  sindexF[0+n]   = last;
+  Int_t countPos = 0;
+
+  // Find frequency
+  for (Int_t i = 1; i < n; i++) {
+    val = inlist[sindexS[i]];
+    if (last == val) {
+      sindexF[countPos]++;
+    }
+    else {      
+      countPos++;
+      sindexF[countPos+n] = val;
+      sindexF[countPos]++;
+      last                = val;
+    }
+  }
+  if (last == val) {
+    countPos++;
+  }
+
+  // Sort according frequency
+  TMath::Sort(countPos,sindexF,sindexS,kTRUE);
+
+  for (Int_t i = 0; i < countPos; i++) {
+    outlist[2*i  ] = sindexF[sindexS[i]+n];
+    outlist[2*i+1] = sindexF[sindexS[i]];
+  }
+
+  delete [] sindexS;
+  delete [] sindexF;
+  
+  return countPos;
+
+}
index 26cb1ce..aa7e361 100644 (file)
@@ -1,5 +1,6 @@
 #ifndef ALITRDTRACKERV1_H
 #define ALITRDTRACKERV1_H
+
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */ 
 
 //  The TRD tracker                                                       //
 //                                                                        //
 //  Authors:                                                              //
+//    Marian Ivanov <M.Ivanov@gsi.de>                                     //
 //    Alex Bercuci <A.Bercuci@gsi.de>                                     //
+//    Jouri Belikov <J.Belikov@cern.ch>                                   //
 //    Markus Fasel <M.Fasel@gsi.de>                                       //
 //                                                                        //
-/////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////
 
-#ifndef ALITRDTRACKER_H
-#include "AliTRDtracker.h"
+#ifndef ALITRACKER_H
+#include "AliTracker.h"
 #endif
 
-#define DEBUG
+#ifndef ALITRDTRACKINGSECTOR_H
+#include "AliTRDtrackingSector.h"
+#endif
 
 /**************************************************************************
  * Class Status see source file                                           *
 class TFile;
 class TTreeSRedirector;
 class TClonesArray;
+class TLinearFitter;
 
 class AliRieman;
 class AliESDEvent;
+class AliCluster;
 
+class AliTRDcluster;
 class AliTRDseedV1;
-class AliTRDstackLayer;
+class AliTRDtrackingChamber;
+class AliTRDchamberTimeBin;
 class AliTRDtrackerFitter;
-class AliTRDrecoParam;
 class AliTRDtrackV1;
 class AliTrackPoint;
-class AliTRDtrackerV1 : public AliTRDtracker
+class AliTRDtrackerV1 : public AliTracker
 {
-
- public:
+public:
        enum{
-               kNTimeBins = 35,
-               kNPlanes = 6,
-               kNSeedPlanes = 4,
-               kMaxTracksStack = 100,
-               kNConfigs = 15
+               kMaxLayersPerSector   = 1000
+               , kMaxTimeBinIndex    = 216
+               , kTrackingSectors    = 18
+               , kNTimeBins          = 35
+               , kNPlanes            = 6
+               , kNSeedPlanes        = 4
+               , kMaxTracksStack     = 100
+               , kNConfigs           = 15
        };
-       AliTRDtrackerV1(AliTRDrecoParam *p = 0x0);
-       AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p);
+       AliTRDtrackerV1();
        virtual ~AliTRDtrackerV1();
   
+  //temporary
+  AliTRDtrackingSector* GetTrackingSector(Int_t sec) {return &fTrSec[sec];}
+  
        Int_t          Clusters2Tracks(AliESDEvent *esd);
+  static TTreeSRedirector* DebugStreamer() {return fgDebugStreamer;}
+  AliCluster*    GetCluster(Int_t index) const;
        static void    GetExtrapolationConfig(Int_t iconfig, Int_t planes[2]);
        static void    GetSeedingConfig(Int_t iconfig, Int_t planes[4]);
+       static TLinearFitter *GetTiltedRiemanFitter();
+       static TLinearFitter *GetTiltedRiemanFitterConstraint();
+       static AliRieman *GetRiemanFitter();
+       static void    FitRieman(AliTRDcluster **clusters, Double_t chi2[2]);
+       static Float_t FitRieman(AliTRDseedV1 *tracklets, Double_t *chi2, Int_t *planes = 0x0);
+       static Float_t FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Double_t zVertex);
+       static Float_t FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigError);
+       
        Int_t          FollowBackProlongation(AliTRDtrackV1 &t);
        Int_t          FollowProlongation(AliTRDtrackV1 &t);
+  Int_t          LoadClusters(TTree *cTree);
        Int_t          PropagateBack(AliESDEvent *event);
+  Int_t          ReadClusters(TClonesArray* &array, TTree *in) const;
        Int_t          RefitInward(AliESDEvent *event);
-       void           SetRecoParam(AliTRDrecoParam *p){fRecoParam = p;}
        void           UnloadClusters();
+  
+  static Int_t   Freq(Int_t n, const Int_t *inlist, Int_t *outlist, Bool_t down); // to be removed 
 
- protected:
-       Double_t       BuildSeedingConfigs(AliTRDstackLayer *layer, Int_t *configs);
-       Int_t          Clusters2TracksSM( AliTRDtracker::AliTRDtrackingSector *sector, AliESDEvent *esd);
-       Int_t          Clusters2TracksStack(AliTRDstackLayer *layer, TClonesArray *esdTrackList);
+protected:
+  Bool_t         AdjustSector(AliTRDtrackV1 *track); 
+       Double_t       BuildSeedingConfigs(AliTRDtrackingChamber **stack, Int_t *configs);
+       static Float_t CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset, Double_t slope);
+       Int_t          Clusters2TracksSM(Int_t sector, AliESDEvent *esd);
+       Int_t          Clusters2TracksStack(AliTRDtrackingChamber **stack, TClonesArray *esdTrackList);
        void           CookLabel(AliKalmanTrack *pt, Float_t wrong) const;
-       Int_t          GetSeedingLayers(AliTRDstackLayer *layers, Double_t *params);
-       void           GetMeanCLStack(AliTRDstackLayer *layers, Int_t *planes, Double_t *params);
+       AliTRDcluster* FindSeedingCluster(AliTRDtrackingChamber *chamber, AliTRDseedV1 *sfit) const;
+       Int_t          GetSeedingLayers(AliTRDtrackingChamber *chamber, Double_t *params);
+       void           GetMeanCLStack(AliTRDtrackingChamber *chamber, Int_t *planes, Double_t *params);
        AliTRDseedV1*  GetTracklet(AliTRDtrackV1 *trk, Int_t plane, Int_t &idx);
-       virtual Bool_t GetTrackPoint(Int_t index, AliTrackPoint &p) const;
-       AliTRDcluster *FindSeedingCluster(AliTRDstackLayer *layers, AliTRDseedV1/*AliRieman*/ *sfit);
-       
-       Double_t       MakeSeedingPlanes(AliTRDstackLayer *layer);
-       AliTRDstackLayer *MakeSeedingLayer(AliTRDstackLayer *layers, Int_t Plane);
-       Int_t          MakeSeeds(AliTRDstackLayer *layers, AliTRDseedV1 *sseed, Int_t *ipar);
-       AliTRDtrackV1*   MakeTrack(AliTRDseedV1 *seeds, Double_t *params);
+       Bool_t         GetTrackPoint(Int_t index, AliTrackPoint &p) const;      
+       Int_t          MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *sseed, Int_t *ipar);
+       AliTRDtrackV1* MakeTrack(AliTRDseedV1 *seeds, Double_t *params);
   Int_t          PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t maxStep);
        Int_t          SetTracklet(AliTRDseedV1 *tracklet);
 
@@ -85,19 +109,37 @@ private:
        AliTRDtrackerV1(const AliTRDtrackerV1 &tracker);
        AliTRDtrackerV1 &operator=(const AliTRDtrackerV1 &tracker);
        Double_t       CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4], Double_t *chi2);
-       void           ImproveSeedQuality(AliTRDstackLayer *layer, AliTRDseedV1 *cseed);
+       Double_t       CalculateTrackLikelihood(AliTRDseedV1 *tracklets, Double_t *chi2);
+       Int_t          ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDseedV1 *tracklet);
 
 private:
+  AliTRDgeometry          *fGeom;                          // Pointer to TRD geometry
+  AliTRDtrackingSector    fTrSec[kTrackingSectors];       // Array of tracking sectors;    
 
+       TClonesArray        *fClusters;                       // List of clusters
+       TClonesArray        *fTracklets;                      // List of tracklets
+       TClonesArray        *fTracks;                         // List of tracks
+  Int_t                    fTimeBinsPerPlane;              // Timebins per plane in track prolongation 
+  
+  // should go to the recoParam
+  static const Double_t    fgkMaxChi2;                     // Max increment in track chi2 
+  static const Float_t     fgkMinClustersInTrack;          // Min number of clusters in track
+  static const Float_t     fgkLabelFraction;               // Min fraction of same label
+  static const Double_t    fgkMaxSnp;                      // Maximal snp for tracking
+  static const Double_t    fgkMaxStep;                     // Maximal step for tracking  
+       
+       // stand alone tracking
        static Double_t      fgTopologicQA[kNConfigs];         //  Topologic quality
        Double_t             fTrackQuality[kMaxTracksStack];  //  Track quality 
        Int_t                fSeedLayer[kMaxTracksStack];     //  Seed layer
        Int_t                fSieveSeeding;                   //! Seeding iterator
-       TClonesArray        *fTracklets;                      // List of tracklets for all sectors
-       AliTRDrecoParam     *fRecoParam;                      //  Reconstruction parameters
-       AliTRDtrackerFitter *fFitter;                         //! Fitter class of the tracker
+       
+       static TLinearFitter *fgTiltedRieman;                   //  Fitter for the tilted Rieman fit without vertex constriant
+       static TLinearFitter *fgTiltedRiemanConstrained;        //  Fitter for the tilted Rieman fit with vertex constraint     
+       static AliRieman     *fgRieman;                             //  Fitter for the untilted Rieman fit
+  static TTreeSRedirector *fgDebugStreamer;                 //!Debug streamer
 
-       ClassDef(AliTRDtrackerV1, 1)                          //  Stand alone tracker development class
+       ClassDef(AliTRDtrackerV1, 2)                          //  TRD tracker development class
 
 };
 #endif
diff --git a/TRD/AliTRDtrackingChamber.cxx b/TRD/AliTRDtrackingChamber.cxx
new file mode 100644 (file)
index 0000000..f3ed5d6
--- /dev/null
@@ -0,0 +1,350 @@
+/**************************************************************************
+ * 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: AliTRDtrackingDebug.cxx 23810 2008-02-08 09:00:27Z hristov $ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+//  Tracking in one chamber                                               //
+//                                                                        //
+//  Authors:                                                              //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                     //
+//    Markus Fasel <M.Fasel@gsi.de>                                       //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDtrackingChamber.h"
+
+#include "TMath.h"
+#include "TMatrixTBase.h"
+#include <TTreeStream.h>
+
+#include "AliTRDReconstructor.h"
+#include "AliTRDrecoParam.h"
+#include "AliTRDtrackerV1.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDpadPlane.h"
+#include "AliTRDcalibDB.h"
+
+ClassImp(AliTRDtrackingChamber)
+
+//_______________________________________________________
+AliTRDtrackingChamber::AliTRDtrackingChamber(Int_t det) :
+       fDetector(det)
+       ,fX0(0.)
+{}  
+
+//_______________________________________________________
+void AliTRDtrackingChamber::Clear(const Option_t *opt)
+{
+       for(Int_t itb=0; itb<kNTimeBins; itb++) fTB[itb].Clear(opt);
+}
+
+//_______________________________________________________
+void AliTRDtrackingChamber::InsertCluster(AliTRDcluster *c, Int_t index)
+{
+       fTB[c->GetLocalTimeBin()].InsertCluster(c, index);
+}
+
+//_______________________________________________________
+Bool_t AliTRDtrackingChamber::Build(AliTRDgeometry *geo)
+{
+// Init chamber and all time bins (AliTRDchamberTimeBin)
+// Calculates radial position of the chamber based on 
+// radial positions of the time bins (calibration/alignment aware)
+//
+       Int_t stack = geo->GetChamber(fDetector);
+       Int_t plane = geo->GetPlane(fDetector);
+       AliTRDpadPlane *pp = geo->GetPadPlane(plane, stack);
+       Double_t zl = pp->GetRow0ROC() - pp->GetRowEndROC();
+       Double_t z0 = geo->GetRow0(plane, stack, 0) - zl;
+       Int_t nrows = pp->GetNrows();
+       
+       Int_t index[50], jtb = 0;
+       for(Int_t itb=0; itb<kNTimeBins; itb++){ 
+               if(!fTB[itb]) continue;
+               fTB[itb].SetRange(z0, zl);
+               fTB[itb].SetNRows(nrows);
+               fTB[itb].BuildIndices();
+               index[jtb++] = itb;
+       }       
+       if(jtb<2) return kFALSE;
+       
+       
+       // ESTIMATE POSITION OF PAD PLANE FOR THIS CHAMBER
+       Double_t x0 = fTB[index[0]].GetX();
+       Double_t x1 = fTB[index[1]].GetX();
+       Double_t dx = (x0 - x1)/(index[1] - index[0]); 
+       fX0 = x0 + dx*(index[0] - (Int_t)AliTRDcalibDB::Instance()->GetT0Average(fDetector));   
+       return kTRUE;
+}
+       
+//_______________________________________________________      
+Int_t AliTRDtrackingChamber::GetNClusters() const
+{
+// Returns number of clusters in chamber
+//
+       Int_t n = 0;
+       for(Int_t itb=0; itb<kNTimeBins; itb++){ 
+               n += Int_t(fTB[itb]);
+       }
+       return n;       
+}      
+
+//_______________________________________________________
+Double_t AliTRDtrackingChamber::GetQuality(Int_t timeBins)
+{
+  //
+  // Calculate chamber 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
+  //
+
+       Int_t ncl   = 0;
+       Int_t nused = 0;
+       Int_t nClLayer;
+       for(int itb=0; itb<kNTimeBins; itb++){
+               if(!(nClLayer = fTB[itb].GetNClusters())) continue;
+               ncl += nClLayer;
+               for(Int_t incl = 0; incl < nClLayer; incl++){
+                       if((fTB[itb].GetCluster(incl))->IsUsed()) nused++;
+               }
+       }
+       
+       // calculate the deviation of the mean number of clusters from the
+       // closest integer values
+       Float_t nclMed = float(ncl-nused)/timeBins;
+       Int_t ncli = Int_t(nclMed);
+       Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1));
+       nclDev -= (nclDev>.5) && ncli ? 1. : 0.;
+       return TMath::Exp(-5.*TMath::Abs(nclDev));
+
+//     // 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);
+
+}
+
+
+//_______________________________________________________
+AliTRDchamberTimeBin *AliTRDtrackingChamber::GetSeedingLayer(AliTRDgeometry *geo)
+{
+  //
+  // Creates a seeding layer
+  //
+       
+       // constants
+       const Int_t kMaxRows = 16;
+       const Int_t kMaxCols = 144;
+       const Int_t kMaxPads = 2304;
+               
+       // Get the geometrical data of the chamber
+       Int_t plane = geo->GetPlane(fDetector);
+       Int_t stack = geo->GetChamber(fDetector);
+       Int_t sector= geo->GetSector(fDetector);
+       AliTRDpadPlane *pp = geo->GetPadPlane(plane, stack);
+       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());
+       Float_t z0 = -1., zl = -1.;
+       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 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 < kNTimeBins; iTime++){
+               if(!(nClusters = fTB[iTime].GetNClusters())) continue;
+               z0 = fTB[iTime].GetZ0();
+               zl = fTB[iTime].GetDZ0();
+               for(Int_t incl = 0; incl < nClusters; incl++){
+                       c = fTB[iTime].GetCluster(incl);        
+                       histogram[c->GetPadRow()][c->GetPadCol()]++;
+                       sigmas[c->GetPadRow()][c->GetPadCol()] += c->GetSigmaZ2();
+               }
+       }
+       
+// Now I have everything in the histogram, do the selection
+       //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 to 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[AliTRDtrackerV1::kMaxTracksStack];
+       Float_t sigcands[AliTRDtrackerV1::kMaxTracksStack];
+       
+       // helper variables
+       Int_t indices[kMaxPads]; memset(indices, -1, sizeof(Int_t)*kMaxPads);
+       Int_t nCandidates = 0;
+       Float_t norm, cogv;
+       // histogram filled -> Select best bins
+       Int_t nPads = nCols * nRows;
+       TMath::Sort(nPads, hvals, indices);                     // bins storing a 0 should not matter
+       // Set Threshold
+       Int_t maximum = hvals[indices[0]];      // best
+       Int_t threshold = Int_t(maximum * AliTRDReconstructor::RecoParam()->GetFindableClusters());
+       Int_t col, row, lower, lower1, upper, upper1;
+       for(Int_t ib = 0; ib < nPads; ib++){
+               if(nCandidates >= AliTRDtrackerV1::kMaxTracksStack){
+                       printf("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, AliTRDtrackerV1::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] = Float_t(cogv) /  norm;
+               // passed the filter
+               cand[nCandidates] = row*nCols + 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++;
+       }
+       if(!nCandidates) return 0x0;
+       
+       Float_t pos[3], sig[2];
+       Short_t signal[7]; memset(&signal[0], 0, 7*sizeof(Short_t));
+       AliTRDchamberTimeBin *fakeLayer = new AliTRDchamberTimeBin(plane, stack, sector, z0, zl);
+       AliTRDcluster *cluster = 0x0;
+       if(nCandidates){
+               UInt_t fakeIndex = 0;
+               for(Int_t ican = 0; ican < nCandidates; ican++){
+                       row = cand[ican] / nCols;
+                       col = cand[ican] % nCols;
+                       //temporary
+                       Int_t n = 0; Double_t x = 0., y = 0., z = 0.;
+                       for(int itb=0; itb<kNTimeBins; itb++){
+                               if(!(nClusters = fTB[itb].GetNClusters())) continue;
+                               for(Int_t incl = 0; incl < nClusters; incl++){
+                                       c = fTB[itb].GetCluster(incl);  
+                                       if(c->GetPadRow() != row) continue;
+                                       if(TMath::Abs(c->GetPadCol() - col) > 2) continue;
+                                       x += c->GetX();
+                                       y += c->GetY();
+                                       z += c->GetZ();
+                                       n++;
+                               }
+                       }
+                       pos[0] = x/n;
+                       pos[1] = y/n;
+                       pos[2] = z/n;
+                       sig[0] = .02;
+                       sig[1] = sigcands[ican];
+                       cluster = new AliTRDcluster(fDetector, 0., pos, sig, 0x0, 3, signal, col, row, 0, 0, 0., 0);
+                       fakeLayer->InsertCluster(cluster, fakeIndex++);
+               }
+       }
+       fakeLayer->SetNRows(nRows);
+       fakeLayer->SetOwner();
+       fakeLayer->BuildIndices();
+       //fakeLayer->PrintClusters();
+       
+       if(AliTRDReconstructor::StreamLevel() >= 3){
+               //TMatrixD hist(nRows, nCols);
+               //for(Int_t i = 0; i < nRows; i++)
+               //      for(Int_t j = 0; j < nCols; j++)
+               //              hist(i,j) = histogram[i][j];
+               TTreeSRedirector &cstreamer = *AliTRDtrackerV1::DebugStreamer();
+               cstreamer << "GetSeedingLayer"
+               << "plane="      << plane
+               << "ymin="       << ymin
+               << "ymax="       << ymax
+               << "zmin="       << zmin
+               << "zmax="       << zmax
+               << "L.="         << fakeLayer
+               //<< "Histogram.=" << &hist
+               << "\n";
+       }
+       
+       return fakeLayer;
+}
+
diff --git a/TRD/AliTRDtrackingChamber.h b/TRD/AliTRDtrackingChamber.h
new file mode 100644 (file)
index 0000000..be945fd
--- /dev/null
@@ -0,0 +1,56 @@
+#ifndef ALITRDTRACKINGCHAMBER_H
+#define ALITRDTRACKINGCHAMBER_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliTRDtrackingChamber.h 22646 2007-11-29 18:13:40Z cblume $ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// Data container for one TRD chamber                                     // 
+//                                                                        // 
+// Authors:                                                               //
+//                                                                        //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                     //
+//                                                                        // 
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ALITRDCHAMBERTIMEBIN_H
+#include "AliTRDchamberTimeBin.h"
+#endif
+
+class AliTRDgeometry;
+
+class AliTRDtrackingChamber
+{
+
+public:
+       enum{
+               kNTimeBins = 35
+       };
+       AliTRDtrackingChamber(Int_t det);
+       virtual ~AliTRDtrackingChamber(){}
+       
+       Bool_t   Build(AliTRDgeometry *geo);
+  void     Clear(const Option_t *opt = 0x0);
+       Int_t    GetDetector() const {return fDetector;}
+       Int_t    GetNClusters() const;
+       Double_t GetQuality(Int_t ntb);
+       AliTRDchamberTimeBin *GetSeedingLayer(AliTRDgeometry *geo);
+       Float_t  GetX()        const {return fX0;}
+       AliTRDchamberTimeBin* GetTB(int tb) {return tb >= 0 && tb < kNTimeBins ? &fTB[tb] : 0x0;}
+       void     InsertCluster(AliTRDcluster *c, Int_t index);
+       
+
+private:
+       Int_t         fDetector;  // detector number
+       Float_t       fX0;        // approximate position of the pad plane
+       
+       AliTRDchamberTimeBin fTB[kNTimeBins];    // time bins 
+       
+       
+       ClassDef(AliTRDtrackingChamber, 1)  // TRD tracker container for one chamber
+};
+
+#endif  // ALITRDTRACKINGCHAMBER_H
diff --git a/TRD/AliTRDtrackingSector.cxx b/TRD/AliTRDtrackingSector.cxx
new file mode 100644 (file)
index 0000000..ccec101
--- /dev/null
@@ -0,0 +1,189 @@
+/**************************************************************************
+ * 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: AliTRDtrackingSector.cxx 23810 2008-02-08 09:00:27Z hristov $ */
+
+///////////////////////////////////////////////////////////////////////////////
+//                                                                           //
+//  Tracking data container for one sector                                   //
+//                                                                           //
+//  Authors:                                                                 //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                        //
+//    Markus Fasel <M.Fasel@gsi.de>                                          //
+//                                                                           //
+///////////////////////////////////////////////////////////////////////////////
+
+#include "AliTRDtrackingSector.h"
+#include "AliTRDcalibDB.h"
+#include "AliTRDCommonParam.h"
+#include "AliTRDgeometry.h"
+#include "AliTRDpadPlane.h"
+#include "AliTRDtrackingChamber.h"
+
+ClassImp(AliTRDtrackingSector)
+
+//_____________________________________________________________________________
+AliTRDtrackingSector::AliTRDtrackingSector()
+  :fTimeBinsPerPlane(0)
+  ,fSector(-1)
+  ,fN(0)
+  ,fGeom(0x0)
+{
+       // Default constructor
+       
+       for(int ic=0; ic<kNChambersSector; ic++){
+               fChamber[ic] = 0x0;
+               fIndex[ic]   = -1;
+       }
+       for(int ip=0; ip<AliTRDgeometry::kNplan; ip++) fX0[ip] = 0.;
+}
+
+//_____________________________________________________________________________
+AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry *geo, Int_t gs, Int_t tb)
+  :fTimeBinsPerPlane(tb)
+  ,fSector(gs)
+  ,fN(0)
+  ,fGeom(geo)
+{
+  //
+  // AliTRDtrackingSector Constructor
+  //
+
+       for(int ic=0; ic<kNChambersSector; ic++){
+               fChamber[ic] = 0x0;
+               fIndex[ic]   = -1;
+       }
+       for(int ip=0; ip<AliTRDgeometry::kNplan; ip++) fX0[ip] = 0.;
+}
+
+//_____________________________________________________________________________
+AliTRDtrackingSector::AliTRDtrackingSector(const AliTRDtrackingSector &/*t*/)
+  :fTimeBinsPerPlane(0)
+  ,fSector(-1)
+  ,fN(0)
+  ,fGeom(0x0)
+{
+  //
+  // Copy constructor
+  //
+
+}
+
+//_____________________________________________________________________________
+AliTRDtrackingSector::~AliTRDtrackingSector()
+{
+  //
+  // Destructor
+  //
+
+}
+               
+//_____________________________________________________________________________
+void AliTRDtrackingSector::Init()
+{              
+//     Steer building of tracking chambers and build tracking sector.
+//     Propagate radial position information (calibration/alignment aware) from chambers to sector level
+//
+       
+       AliTRDtrackingChamber *tc = 0x0; int ic = 0; 
+       while((tc = fChamber[ic++])) tc->Build(fGeom);
+               
+       Int_t np;
+       for(int ip=0; ip<AliTRDgeometry::kNplan; ip++){
+               fX0[ip] = 0.; np = 0;
+               for(int is=0; is<AliTRDgeometry::kNcham; is++){
+                       Int_t idx = is*AliTRDgeometry::kNplan + ip;
+                       if(fIndex[idx]<0) continue;
+                       tc = GetChamber(fIndex[idx]);
+                       fX0[ip] += tc->GetX(); np++; 
+               }
+               if(!np){
+                       //printf("Could not estimate radial position  of plane %d in sector %d.\n", ip, fSector);
+                       continue;
+               }
+               fX0[ip] /= Float_t(np);
+       }
+}
+//_____________________________________________________________________________
+void AliTRDtrackingSector::Clear(const Option_t *opt)
+{
+// Reset counters and steer chamber clear
+
+       for(Int_t ich=0; ich<fN; ich++){ 
+               fChamber[ich]->Clear(opt);
+               delete fChamber[ich]; fChamber[ich] = 0x0;   // I would avoid
+       }       
+       for(Int_t ich=0; ich<kNChambersSector; ich++) fIndex[ich] = -1;
+       fN = 0;
+}
+
+//_____________________________________________________________________________
+AliTRDtrackingChamber* AliTRDtrackingSector::GetChamber(Int_t stack, Int_t plane, Bool_t build)
+{
+// Return chamber at position (stack, plane) in current 
+// sector or build a new one if it is not already created
+       
+       Int_t ch = stack*AliTRDgeometry::kNplan + plane;
+       if(fIndex[ch] >= 0) return fChamber[Int_t(fIndex[ch])];
+       else if(!build) return 0x0;
+       
+       // CHAMBER HAS TO BE BUILD
+       Int_t rch = ch;do rch--; while(rch>=0 && fIndex[rch]<0);
+       fIndex[ch] = rch >=0 ? fIndex[rch]+1 : 0; 
+       fN++;
+       
+       memmove(&fChamber[Int_t(fIndex[ch])+1], &fChamber[Int_t(fIndex[ch])], (kNChambersSector-fIndex[ch]-1)*sizeof(void*));
+       for(Int_t ic = ch+1; ic<kNChambersSector; ic++) fIndex[ic] += fIndex[ic] >= 0 ? 1 : 0;
+       
+       return fChamber[Int_t(fIndex[ch])] = new AliTRDtrackingChamber(AliTRDgeometry::GetDetector(plane, stack, fSector));
+}
+
+//_____________________________________________________________________________
+AliTRDtrackingChamber** AliTRDtrackingSector::GetStack(Int_t stack)
+{
+// Return chamber at position (stack, plane) in current 
+// sector or build a new one if it is not already created
+       
+       if(stack<0 || stack>=AliTRDgeometry::kNcham) return 0x0;
+       
+       Int_t ich, n = 0;
+       for(int ip=0; ip<AliTRDgeometry::kNplan; ip++){
+               ich = stack*AliTRDgeometry::kNplan + ip;
+               if(fIndex[ich] < 0) fStack[ip] = 0x0; 
+               else{
+                       fStack[ip] = fChamber[Int_t(fIndex[ich])];
+                       n++;
+               }
+       }
+       
+       return n ? &fStack[0] : 0x0;
+}
+
+//_____________________________________________________________________________
+void AliTRDtrackingSector::Print(Option_t *)
+{
+// Dump info about this tracking sector and the tracking chamber within
+// 
+
+       printf("\tSector %2d\n", fSector);
+       for(int ip=0; ip<6; ip++){
+               for(int is =0; is<5; is++){
+                       Int_t ch = is*AliTRDgeometry::kNplan + ip;
+                       printf("%2d[%2d] ", fIndex[ch], fIndex[ch]>=0 ? fChamber[Int_t(fIndex[ch])]->GetNClusters() : 0);
+               }
+               printf("\n");
+       }
+
+}
diff --git a/TRD/AliTRDtrackingSector.h b/TRD/AliTRDtrackingSector.h
new file mode 100644 (file)
index 0000000..9f37e8b
--- /dev/null
@@ -0,0 +1,76 @@
+#ifndef ALITRDTRACKINGSECTOR_H
+#define ALITRDTRACKINGSECTOR_H
+
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliTRDtrackingSector.h 22646 2007-11-29 18:13:40Z cblume $ */
+
+////////////////////////////////////////////////////////////////////////////
+//                                                                        //
+// Data container for one TRD supermodule                                 // 
+//                                                                        // 
+// Authors:                                                               //
+//                                                                        //
+//    Marian Ivanov <M.Ivanov@gsi.de>                                     //
+//    Alex Bercuci <A.Bercuci@gsi.de>                                     //
+//                                                                        //
+////////////////////////////////////////////////////////////////////////////
+
+#ifndef ROOT_Rtypes
+#include "Rtypes.h"
+#endif
+
+class AliTRDtrackingChamber;
+class AliTRDgeometry;
+
+class AliTRDtrackingSector 
+{
+
+public:
+       enum{
+               kNChambersSector   = 30
+               , kNplane          = 6
+               , kTrackingSectors = 18
+       };
+       
+       AliTRDtrackingSector();
+       AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs, Int_t tbs);
+       AliTRDtrackingSector(const AliTRDtrackingSector &/*t*/);
+       virtual ~AliTRDtrackingSector();
+       
+       AliTRDtrackingSector &operator=(const AliTRDtrackingSector &/*t*/) { return *this; }
+       
+  void     Clear(const Option_t *opt = 0x0);
+       Int_t    GetNChambers() const             { return fN; }
+       Double_t GetX(Int_t pl) const                  { return pl >=0 && pl < kNplane ? fX0[pl] : 0.; }
+       AliTRDtrackingChamber* GetChamber(Int_t i) const  { return i>=0 && i < fN ? fChamber[i] : 0x0;  }
+       AliTRDtrackingChamber* GetChamber(Int_t stack, Int_t plane, Bool_t build = kFALSE);
+       AliTRDtrackingChamber** GetStack(Int_t stack);
+       Int_t    GetSector() const {return fSector;}    
+       void     Init();
+       
+       void     SetGeometry(AliTRDgeometry *geo) {fGeom = geo;}
+        
+       // temporary ... some of this functions are obsolete
+       void Print(Option_t *opt = 0x0);
+//     void     MapTimeBinLayers();
+//     Int_t    Find(Double_t x) const; 
+//     void     InsertLayer(AliTRDchamberTimeBin *pl);
+//     Int_t    CookTimeBinIndex(Int_t plane, Int_t localTB) const;     
+       
+private:
+       UChar_t        fTimeBinsPerPlane; // no. of time bins
+       Char_t         fSector;           // Sector# in AliTRDgeometry
+       UChar_t        fN;                // Total number of chambers allocated
+       Char_t         fIndex[kNChambersSector]; // indexes of allocated chambers
+       Float_t        fX0[kNplane];      // average position of pad plane for each plane
+       AliTRDgeometry *fGeom;            // Geometry
+       AliTRDtrackingChamber *fChamber[kNChambersSector];// chambers   
+       AliTRDtrackingChamber *fStack[kNplane]; //! temporary holding one stack
+       ClassDef(AliTRDtrackingSector, 1) // TRD tracker container for one sector
+};
+
+
+#endif // ALITRDTRACKINGSECTOR_H
+