A major update in the tracking code is available in the TRUNK. This
authorabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Mar 2009 15:44:34 +0000 (15:44 +0000)
committerabercuci <abercuci@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 19 Mar 2009 15:44:34 +0000 (15:44 +0000)
update is mainly triggered to create the infrastructure for
recalculating the cluster position during tracking but it is also
answering some comments/suggestions from the users. Here is a list of
the features which are now available:

- use of GeoManager in the tracking. This is needed to correctly
  initialize the radial position of the following :
      - anode wire plane
      - chamber entrance

- configurable fiducial area. for the moment I have a 5mm exclusion area around the pad plane limits

- proper definition of ESD status bits, in particular the TRDin and TRDStop

- history of tracking per TRD track. One would be able now
  to query what happened to the track in each layer on the
  way from TPC to TOF. At the level of AliTRDtrackV1 one
  can access this information via the function
   inline UChar_t GetStatusTRD(Int_t ly=-1) const;
  There are also general (static) functions to be used in the future with the ESD track:
   inline static Bool_t IsTrackError(ETRDtrackError error, UInt_t status);
   inline static Bool_t IsLayerError(ETRDlayerError error, Int_t layer, UInt_t status);

TRD/AliTRDgeometry.cxx
TRD/AliTRDgeometry.h
TRD/AliTRDseedV1.cxx
TRD/AliTRDseedV1.h
TRD/AliTRDtrackV1.cxx
TRD/AliTRDtrackV1.h
TRD/AliTRDtrackerV1.cxx
TRD/AliTRDtrackerV1.h
TRD/qaRec/AliTRDpidChecker.cxx

index 0c2df28..ad7ee57 100644 (file)
@@ -2881,5 +2881,39 @@ Bool_t AliTRDgeometry::IsHole(Int_t /*la*/, Int_t st, Int_t se) const
     return kTRUE; 
   }
   return kFALSE;
+}
 
+//_____________________________________________________________________________
+Bool_t AliTRDgeometry::IsOnBoundary(Int_t det, Float_t y, Float_t z, Float_t eps) const
+{
+  Int_t ly = GetLayer(det);
+  if ((ly <          0) || 
+      (ly >= fgkNlayer)) return kTRUE;
+       
+  Int_t stk = GetStack(det);
+  if ((stk <          0) || 
+      (stk >= fgkNstack)) return kTRUE;
+
+  AliTRDpadPlane *pp = (AliTRDpadPlane*) fPadPlaneArray->At(GetDetectorSec(ly, stk));
+  if(!pp) return kTRUE;
+
+  Double_t max  = pp->GetRow0();
+  Int_t n = pp->GetNrows();
+  Double_t min = max - 2 * pp->GetLengthOPad() 
+                 - (n-2) * pp->GetLengthIPad() 
+                 - (n-1) * pp->GetRowSpacing();
+  if(z < min+eps || z > max-eps){ 
+    //printf("z : min[%7.2f (%7.2f)] %7.2f max[(%7.2f) %7.2f]\n", min, min+eps, z, max-eps, max);
+    return kTRUE;
+  }
+  min  = pp->GetCol0();
+  n = pp->GetNcols();
+  max = min +2 * pp->GetWidthOPad() 
+       + (n-2) * pp->GetWidthIPad() 
+       + (n-1) * pp->GetColSpacing();
+  if(y < min+eps || y > max-eps){ 
+    //printf("y : min[%7.2f (%7.2f)] %7.2f max[(%7.2f) %7.2f]\n", min, min+eps, y, max-eps, max);
+    return kTRUE;
+  }
+  return kFALSE;
 }
index e80a171..af7bcfd 100644 (file)
@@ -39,7 +39,7 @@ class AliTRDgeometry : public AliGeometry {
   virtual Int_t    IsVersion()                                            { return 1;               }
   virtual Bool_t   Impact(const TParticle* ) const                        { return kTRUE;           }
   virtual Bool_t   IsHole(Int_t la, Int_t st, Int_t se) const;
-
+  virtual Bool_t   IsOnBoundary(Int_t det, Float_t y, Float_t z, Float_t eps = .5) const;
   virtual Bool_t   RotateBack(Int_t det, Double_t *loc, Double_t *glb) const;
 
           Bool_t   ChamberInGeometry(Int_t det);
index ec20f45..32a9026 100644 (file)
@@ -710,6 +710,17 @@ void AliTRDseedV1::SetOwner()
   SetBit(kOwner);
 }
 
+//____________________________________________________________
+void AliTRDseedV1::SetPadPlane(AliTRDpadPlane *p)
+{
+// Shortcut method to initialize pad geometry.
+  if(!p) return;
+  SetTilt(TMath::Tan(TMath::DegToRad()*p->GetTiltingAngle()));
+  SetPadLength(p->GetLengthIPad());
+  SetPadWidth(p->GetWidthIPad());
+}
+
+
 // //____________________________________________________________________
 // Bool_t      AliTRDseedV1::AttachClustersIter(AliTRDtrackingChamber *chamber, Float_t quality, Bool_t kZcorr, AliTRDcluster *c)
 // {
@@ -1025,22 +1036,6 @@ Bool_t   AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *chamber, Bool_t tilt)
   Int_t dtb = tb[1] - tb[0];
   fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
 
-  // update X0 from the clusters (calibration/alignment aware) TODO remove dependence on x0 !!
-  for (Int_t it = 0; it < AliTRDtrackerV1::GetNTimeBins(); it++) {
-    if(!(layer = chamber->GetTB(it))) continue;
-    if(!layer->IsT0()) continue;
-    if(fClusters[it]){ 
-      fX0 = fClusters[it]->GetX();
-      break;
-    } else { // we have to infere the position of the anode wire from the other clusters
-      for (Int_t jt = it+1; jt < AliTRDtrackerV1::GetNTimeBins(); jt++) {
-        if(!fClusters[jt]) continue;
-        fX0 = fClusters[jt]->GetX() + fdX * (jt - it);
-        break;
-      }
-    }
-  }    
-
   return kTRUE;
 }
 
@@ -1564,7 +1559,7 @@ void AliTRDseedV1::Print(Option_t *o) const
   // Printing the seedstatus
   //
 
-  AliInfo(Form("Det[%3d] Pad[L[%5.2f] W[%5.2f] Tilt[%+6.2f]]", fDet, GetPadLength(), GetPadWidth(), GetTilt()));
+  AliInfo(Form("Det[%3d] X0[%7.2f] Pad[L[%5.2f] W[%5.2f] Tilt[%+6.2f]]", fDet, fX0, GetPadLength(), GetPadWidth(), GetTilt()));
   AliInfo(Form("N[%2d] Nused[%2d] Nshared[%2d] [%d]", GetN(), GetNUsed(), GetNShared(), fN));
 
   Double_t cov[3], x=GetX();
index 4d43f2e..98091d5 100644 (file)
@@ -42,6 +42,7 @@ class AliRieman;
 class AliTRDtrackingChamber;
 class AliTRDtrackV1;
 class AliTRDReconstructor;
+class AliTRDpadPlane;
 class AliTRDseedV1 : public AliTRDtrackletBase
 {
 public:
@@ -156,6 +157,7 @@ public:
   void      SetStandAlone(Bool_t st) { SetBit(kStandAlone, st); }
   void      SetMomentum(Double_t mom){ fMom = mom;}
   void      SetOwner();
+  void      SetPadPlane(AliTRDpadPlane *p);
   void      SetPadLength(Float_t l)  { fPad[0] = l;}
   void      SetPadWidth(Float_t w)   { fPad[1] = w;}
   void      SetTilt(Float_t tilt)    { fPad[2] = tilt; }
index bdf9f45..97ec135 100644 (file)
@@ -40,7 +40,7 @@ ClassImp(AliTRDtrackV1)
 
 //_______________________________________________________________
 AliTRDtrackV1::AliTRDtrackV1() : AliKalmanTrack()
-  ,fPIDquality(0)
+  ,fStatus(0)
   ,fDE(0.)
   ,fReconstructor(0x0)
   ,fBackupTrack(0x0)
@@ -65,7 +65,7 @@ AliTRDtrackV1::AliTRDtrackV1() : AliKalmanTrack()
 
 //_______________________________________________________________
 AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref) : AliKalmanTrack(ref)
-  ,fPIDquality(ref.fPIDquality)
+  ,fStatus(ref.fStatus)
   ,fDE(ref.fDE)
   ,fReconstructor(ref.fReconstructor)
   ,fBackupTrack(0x0)
@@ -94,7 +94,7 @@ AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref) : AliKalmanTrack(ref)
 
 //_______________________________________________________________
 AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) : AliKalmanTrack()
-  ,fPIDquality(0)
+  ,fStatus(0)
   ,fDE(0.)
   ,fReconstructor(0x0)
   ,fBackupTrack(0x0)
@@ -144,7 +144,7 @@ AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) : AliKalmanTrack()
 //_______________________________________________________________
 AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Double_t cov[15]
              , Double_t x, Double_t alpha) : AliKalmanTrack()
-  ,fPIDquality(0)
+  ,fStatus(0)
   ,fDE(0.)
   ,fReconstructor(0x0)
   ,fBackupTrack(0x0)
@@ -283,48 +283,78 @@ Bool_t AliTRDtrackV1::CookPID()
   
   // Reset the a priori probabilities
   Double_t pid = 1. / AliPID::kSPECIES;
-  for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
-    fPID[ispec] = pid; 
+  for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fPID[ispec] = pid; 
+
+  UChar_t fPIDquality = SetNumberOfTrackletsPID(kTRUE);
+  // no tracklet found for PID calculations
+  if(!fPIDquality) return kFALSE;
+  
+  // slot for PID calculation @ track level
+  
+  // normalize probabilities
+  Double_t probTotal = 0.0;
+  for (Int_t is = 0; is < AliPID::kSPECIES; is++) probTotal += fPID[is];
+  
+  
+  if (probTotal <= 0.0) {
+    AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference data.");
+    return kFALSE;
   }
-  fPIDquality = 0;
   
-  // steer PID calculation @ tracklet level
+  for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) fPID[iSpecies] /= probTotal;
+  
+  return kTRUE;
+}
+
+//___________________________________________________________
+UChar_t AliTRDtrackV1::GetNumberOfTrackletsPID() const
+{
+// Retrieve number of tracklets used for PID calculation. 
+
+  UChar_t fPIDquality = 0;
   Float_t *prob = 0x0;
   for(int ip=0; ip<kNplane; ip++){
     if(fTrackletIndex[ip] == 0xffff) continue;
     if(!fTracklet[ip]->IsOK()) continue;
-    if(!(prob = fTracklet[ip]->GetProbability(kTRUE))) return kFALSE;
+    if(!(prob = fTracklet[ip]->GetProbability(kFALSE))) return 0;
     
     Int_t nspec = 0; // quality check of tracklet dEdx
     for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
       if(prob[ispec] < 0.) continue;
-      fPID[ispec] *= prob[ispec];
       nspec++;
     }
     if(!nspec) continue;
     
     fPIDquality++;
   }
+  return fPIDquality;
+}
+
+//___________________________________________________________
+UChar_t AliTRDtrackV1::SetNumberOfTrackletsPID(Bool_t recalc)
+{
+// Retrieve number of tracklets used for PID calculation. // Recalculated PID at tracklet level by quering the PID DB.
+
+  UChar_t fPIDquality = 0;
   
-  // no tracklet found for PID calculations
-  if(!fPIDquality) return kTRUE;
-  
-  // slot for PID calculation @ track level
-  
-  
-  // normalize probabilities
-  Double_t probTotal = 0.0;
-  for (Int_t is = 0; is < AliPID::kSPECIES; is++) probTotal += fPID[is];
-  
-  
-  if (probTotal <= 0.0) {
-    AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference data.");
-    return kFALSE;
+  // steer PID calculation @ tracklet level
+  Float_t *prob = 0x0;
+  for(int ip=0; ip<kNplane; ip++){
+    if(fTrackletIndex[ip] == 0xffff) continue;
+    if(!fTracklet[ip]->IsOK()) continue;
+    if(!(prob = fTracklet[ip]->GetProbability(recalc))) return 0;
+    
+    Int_t nspec = 0; // quality check of tracklet dEdx
+    for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
+      if(prob[ispec] < 0.) continue;
+      fPID[ispec] *= prob[ispec];
+      nspec++;
+    }
+    if(!nspec) continue;
+    
+    fPIDquality++;
   }
-  
-  for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) fPID[iSpecies] /= probTotal;
-  
-  return kTRUE;
+  return fPIDquality;
 }
 
 //_______________________________________________________________
@@ -389,13 +419,19 @@ Double_t AliTRDtrackV1::GetPredictedChi2(const AliTRDseedV1 *trklt) const
 }
 
 //_______________________________________________________________
+Int_t AliTRDtrackV1::GetSector() const
+{
+  return Int_t(GetAlpha()/AliTRDgeometry::GetAlpha() + (GetAlpha()>0. ? 0 : AliTRDgeometry::kNsector));
+}
+
+//_______________________________________________________________
 Bool_t AliTRDtrackV1::IsEqual(const TObject *o) const
 {
   if (!o) return kFALSE;
   const AliTRDtrackV1 *inTrack = dynamic_cast<const AliTRDtrackV1*>(o);
   if (!inTrack) return kFALSE;
   
-  if ( fPIDquality != inTrack->GetPIDquality() ) return kFALSE;
+  //if ( fPIDquality != inTrack->GetPIDquality() ) return kFALSE;
   
   for(Int_t i = 0; i < AliPID::kSPECIES; i++){
     if ( fPID[i] != inTrack->GetPID(i) ) return kFALSE;
@@ -629,11 +665,11 @@ Int_t   AliTRDtrackV1::PropagateToR(Double_t r,Double_t step)
 //_____________________________________________________________________________
 void AliTRDtrackV1::Print(Option_t *o) const
 {
-  AliInfo(Form("PID q[%d] [%4.1f %4.1f %4.1f %4.1f %4.1f]", fPIDquality, 1.E2*fPID[0], 1.E2*fPID[1], 1.E2*fPID[2], 1.E2*fPID[3], 1.E2*fPID[4]));
+  AliInfo(Form("PID [%4.1f %4.1f %4.1f %4.1f %4.1f]", 1.E2*fPID[0], 1.E2*fPID[1], 1.E2*fPID[2], 1.E2*fPID[3], 1.E2*fPID[4]));
   AliInfo(Form("Material[%5.2f %5.2f %5.2f]", fBudget[0], fBudget[1], fBudget[2]));
 
   AliInfo(Form("x[%7.2f] t[%7.4f] alpha[%f] mass[%f]", GetX(), GetIntegratedLength(), GetAlpha(), fMass));
-  AliInfo(Form("Ntr[%1d] Ncl[%3d] lab[%3d]", GetNumberOfTracklets(), fN, fLab));
+  AliInfo(Form("Ntr[%1d] NtrPID[%1d] Ncl[%3d] lab[%3d]", GetNumberOfTracklets(), GetNumberOfTrackletsPID(), fN, fLab));
 
   if(strcmp(o, "a")!=0) return;
   printf("|X| = (");
@@ -800,10 +836,11 @@ void AliTRDtrackV1::UpdateESDtrack(AliESDtrack *track)
     for (Int_t js = 0; js < nslices; js++, dedx++) track->SetTRDslice(*dedx, ip, js);
   }
 
-  if(!fPIDquality) track->SetTRDntracklets(n);
+  UChar_t nPID = GetNumberOfTrackletsPID();
+  if(!nPID) track->SetTRDntracklets(n);
   else {
     track->SetTRDpid(fPID);
-    n |= (fPIDquality<<3);
+    n |= (nPID<<3);
     track->SetTRDntracklets(n);
   }
 }
index 72e2eeb..c44ac1f 100644 (file)
@@ -29,7 +29,7 @@ class AliTRDReconstructor;
 class AliTRDtrackV1 : public AliKalmanTrack
 {
 public:
-  enum ETRDtrackV1Size { 
+  enum ETRDtrackSize { 
     kNdet      = AliTRDgeometry::kNdet
    ,kNstacks   = AliTRDgeometry::kNstack*AliTRDgeometry::kNsector
    ,kNplane    = AliTRDgeometry::kNlayer
@@ -41,10 +41,31 @@ public:
   };
   
   // bits from 0-13 are reserved by ROOT (see TObject.h)
-  enum ETRDtrackV1Status {
-    kOwner   = BIT(14)
-   ,kStopped = BIT(15) 
-   ,kKink    = BIT(16) 
+  enum ETRDtrackStatus {
+    kOwner     = BIT(14)
+   ,kStopped   = BIT(15) 
+   ,kKink      = BIT(16) 
+  };
+
+  // propagation/update error codes (up to 4 bits)
+  enum ETRDtrackError {
+    kProlongation = 1
+   ,kPropagation
+   ,kAdjustSector
+   ,kSnp
+   ,kTrackletInit
+   ,kUpdate
+  };
+
+  // data/clusters/tracklet error codes (up to 4 bits/layer)
+  enum ETRDlayerError {
+    kGeometry = 1
+   ,kBoundary
+   ,kNoClusters
+   ,kNoAttach
+   ,kNoClustersTracklet
+   ,kNoFit
+   ,kChi2
   };
 
   AliTRDtrackV1();
@@ -66,10 +87,12 @@ public:
   inline Int_t   GetNumberOfTracklets() const;
   Double_t       GetPIDsignal() const   { return 0.;}
   Double_t       GetPID(Int_t is) const { return (is >=0 && is < AliPID::kSPECIES) ? fPID[is] : -1.;}
-  UChar_t        GetPIDquality() const  { return fPIDquality;}
+  UChar_t        GetNumberOfTrackletsPID() const;
   Double_t       GetPredictedChi2(const AliTRDseedV1 *tracklet) const;
   Double_t       GetPredictedChi2(const AliCluster* /*c*/) const                   { return 0.0; }
   Int_t          GetProlongation(Double_t xk, Double_t &y, Double_t &z);
+  inline UChar_t GetStatusTRD(Int_t ly=-1) const;
+  Int_t          GetSector() const;
   AliTRDseedV1*  GetTracklet(Int_t plane) const {return plane >=0 && plane <kNplane ? fTracklet[plane] : 0x0;}
   Int_t          GetTrackletIndex(Int_t plane) const          { return (plane>=0 && plane<kNplane) ? fTrackletIndex[plane] : -1;}
   AliExternalTrackParam*
@@ -83,7 +106,9 @@ public:
   Bool_t         IsOwner() const   { return TestBit(kOwner);};
   Bool_t         IsStopped() const { return TestBit(kStopped);};
   Bool_t         IsElectron() const;
-  
+  inline static Bool_t IsTrackError(ETRDtrackError error, UInt_t status);
+  inline static Bool_t IsLayerError(ETRDlayerError error, Int_t layer, UInt_t status);
+
   void           MakeBackupTrack();
   void           Print(Option_t *o="") const;
 
@@ -94,9 +119,11 @@ public:
   void           SetEdep(Double32_t inDE){fDE = inDE;};
   void           SetKink(Bool_t k)        { SetBit(kKink, k);}
   void           SetNumberOfClusters();
+  UChar_t        SetNumberOfTrackletsPID(Bool_t recalc);
   void           SetOwner();
   void           SetPID(Short_t is, Double_t inPID){if (is >=0 && is < AliPID::kSPECIES) fPID[is]=inPID;};
-  void           SetPIDquality(UChar_t inPIDquality){fPIDquality = inPIDquality;};
+  void           SetPIDquality(UChar_t /*inPIDquality*/){/*fPIDquality = inPIDquality*/;};
+  inline void    SetStatus(UChar_t stat, Int_t ly=-1);
   void           SetStopped(Bool_t stop) {SetBit(kStopped, stop);}
   void           SetTracklet(AliTRDseedV1 *trklt,  Int_t index);
   void           SetTrackLow();
@@ -110,18 +137,18 @@ public:
   void           UpdateESDtrack(AliESDtrack *t);
 
 private:
-  UChar_t      fPIDquality;           //  No of planes used for PID calculation        
+  UInt_t       fStatus;                //  Bit map for the status of propagation
   UShort_t     fTrackletIndex[kNplane];//  Tracklets index in the tracker list
-  Double32_t   fPID[AliPID::kSPECIES];//  PID probabilities
-  Double32_t   fBudget[3];            //  Integrated material budget
-  Double32_t   fDE;                   //  Integrated delta energy
+  Double32_t   fPID[AliPID::kSPECIES]; //  PID probabilities
+  Double32_t   fBudget[3];             //  Integrated material budget
+  Double32_t   fDE;                    //  Integrated delta energy
   const AliTRDReconstructor *fReconstructor;//! reconstructor link 
-  AliTRDseedV1 *fTracklet[kNplane];   //  Tracklets array defining the track
-  AliTRDtrackV1 *fBackupTrack;        // Backup track
-  AliExternalTrackParam *fTrackLow;   // parameters of the track which enter TRD from below (TPC) 
+  AliTRDtrackV1 *fBackupTrack;         //! Backup track
+  AliTRDseedV1 *fTracklet[kNplane];    //  Tracklets array defining the track
+  AliExternalTrackParam *fTrackLow;    // parameters of the track which enter TRD from below (TPC) 
   AliExternalTrackParam *fTrackHigh;  // parameters of the track which enter TRD from above (HMPID, PHOS) 
 
-  ClassDef(AliTRDtrackV1, 4)          // new TRD track
+  ClassDef(AliTRDtrackV1, 5)          // TRD track - tracklet based
 };
 
 //____________________________________________________
@@ -152,6 +179,25 @@ inline Int_t AliTRDtrackV1::GetNumberOfTracklets() const
   return n;
 }
 
+//____________________________________________________
+inline UChar_t AliTRDtrackV1::GetStatusTRD(Int_t ly) const
+{
+  if(ly<kNplane) return (fStatus>>((ly+1)*4))&0xf;
+  return -1;
+}
+
+//____________________________________________________
+inline Bool_t AliTRDtrackV1::IsTrackError(ETRDtrackError error, UInt_t status)
+{
+  return (status&0xf)==UChar_t(error);
+}
+
+//____________________________________________________
+inline Bool_t AliTRDtrackV1::IsLayerError(ETRDlayerError error, Int_t ly, UInt_t status)
+{
+  if(ly>=kNplane || ly<0) return kFALSE;
+  return ((status>>((ly+1)*4))&0xf) == UChar_t(error);
+}
 
 //____________________________________________________
 inline void AliTRDtrackV1::SetReconstructor(const AliTRDReconstructor *rec)
@@ -163,11 +209,18 @@ inline void AliTRDtrackV1::SetReconstructor(const AliTRDReconstructor *rec)
   fReconstructor = rec;
 }
 
+//____________________________________________________
+inline void AliTRDtrackV1::SetStatus(UChar_t status, Int_t ly)
+{
+  if(ly<kNplane) fStatus|=((status&0xf)<<((ly+1)*4));
+  return;
+}
+
 
 //____________________________________________________________________________
 inline Float_t AliTRDtrackV1::StatusForTOF()
 {
-  //
+  // OBSOLETE
   // Defines the status of the TOF extrapolation
   //
 
index d32fee1..7d5e949 100644 (file)
@@ -35,6 +35,8 @@
 #include <TTree.h>  
 #include <TClonesArray.h>
 #include <TTreeStream.h>
+#include <TGeoMatrix.h>
+#include <TGeoManager.h>
 
 #include "AliLog.h"
 #include "AliMathBase.h"
@@ -82,7 +84,7 @@ TLinearFitter* AliTRDtrackerV1::fgTiltedRiemanConstrained = 0x0;
 AliTRDtrackerV1::AliTRDtrackerV1(AliTRDReconstructor *rec) 
   :AliTracker()
   ,fReconstructor(0x0)
-  ,fGeom(new AliTRDgeometry())
+  ,fGeom(0x0)
   ,fClusters(0x0)
   ,fTracklets(0x0)
   ,fTracks(0x0)
@@ -91,19 +93,41 @@ AliTRDtrackerV1::AliTRDtrackerV1(AliTRDReconstructor *rec)
   //
   // Default constructor.
   // 
+  
+  SetReconstructor(rec); // initialize reconstructor
+
+  // initialize geometry
+  if(!AliGeomManager::GetGeometry()){
+    AliFatal("Could not get geometry.");
+  }
+  fGeom = new AliTRDgeometry();
+  fGeom->CreateClusterMatrixArray();
+  TGeoHMatrix *matrix = 0x0;
+  Double_t loc[] = {0., 0., 0.};
+  Double_t glb[] = {0., 0., 0.};
+  for(Int_t ily=kNPlanes; ily--;){
+    if(!(matrix = fGeom->GetClusterMatrix(AliTRDgeometry::GetDetector(ily, 2, 0)))){
+      AliFatal(Form("Could not get matrix for chamber @ 0 2 %d.", ily));
+      continue;
+    }
+    matrix->LocalToMaster(loc, glb);
+    fR[ily] = glb[0]+ AliTRDgeometry::AnodePos()-.5*AliTRDgeometry::AmThick() - AliTRDgeometry::DrThick();
+  }
+
+  // initialize calibration values
   AliTRDcalibDB *trd = 0x0;
   if (!(trd = AliTRDcalibDB::Instance())) {
-    AliFatal("Could not get calibration object");
+    AliFatal("Could not get calibration.");
   }
-
   if(!fgNTimeBins) fgNTimeBins = trd->GetNumberOfTimeBins();
 
+  // initialize cluster containers
   for (Int_t isector = 0; isector < AliTRDgeometry::kNsector; isector++) new(&fTrSec[isector]) AliTRDtrackingSector(fGeom, isector);
   
-  for(Int_t isl =0; isl<kNSeedPlanes; isl++) fSeedTB[isl] = 0x0;
-
-  // Initialize debug stream
-  if(rec) SetReconstructor(rec);
+  // initialize arrays
+  memset(fTrackQuality, 0, kMaxTracksStack*sizeof(Double_t));
+  memset(fSeedLayer, 0, kMaxTracksStack*sizeof(Int_t));
+  memset(fSeedTB, 0, kNSeedPlanes*sizeof(AliTRDchamberTimeBin*));
 }
 
 //____________________________________________________________________
@@ -248,61 +272,87 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
   AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
   if (!calibra) AliInfo("Could not get Calibra instance\n");
   
-  Int_t   found    = 0;     // number of tracks found
+  // Define scalers
+  Int_t nFound   = 0, // number of tracks found
+        nSeeds   = 0, // total number of ESD seeds
+        nTRDseeds= 0, // number of seeds in the TRD acceptance
+        nTPCseeds= 0; // number of TPC seeds
   Float_t foundMin = 20.0;
   
   Float_t *quality = 0x0;
   Int_t   *index   = 0x0;
-  Int_t    nSeed   = event->GetNumberOfTracks();
-  if(nSeed){  
-    quality = new Float_t[nSeed];
-    index   = new Int_t[nSeed];
-    for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
+  nSeeds   = event->GetNumberOfTracks();
+  // Sort tracks according to quality 
+  // (covariance in the yz plane)
+  if(nSeeds){  
+    quality = new Float_t[nSeeds];
+    index   = new Int_t[nSeeds];
+    for (Int_t iSeed = nSeeds; iSeed--;) {
       AliESDtrack *seed = event->GetTrack(iSeed);
       Double_t covariance[15];
       seed->GetExternalCovariance(covariance);
       quality[iSeed] = covariance[0] + covariance[2];
     }
-    // Sort tracks according to covariance of local Y and Z
-    TMath::Sort(nSeed, quality, index,kFALSE);
+    TMath::Sort(nSeeds, quality, index,kFALSE);
   }
   
-  // Backpropagate all seeds
+  // Propagate all seeds
   Int_t   expectedClr;
   AliTRDtrackV1 track;
-  for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
+  for (Int_t iSeed = 0; iSeed < nSeeds; iSeed++) {
   
     // Get the seeds in sorted sequence
     AliESDtrack *seed = event->GetTrack(index[iSeed]);
+    Float_t p4  = seed->GetC(seed->GetBz());
   
     // Check the seed status
     ULong_t status = seed->GetStatus();
     if ((status & AliESDtrack::kTPCout) == 0) continue;
     if ((status & AliESDtrack::kTRDout) != 0) continue;
-  
-    // Do the back prolongation
+
+    // Propagate to the entrance in the TRD mother volume
     new(&track) AliTRDtrackV1(*seed);
-    track.SetReconstructor(fReconstructor);
-    track.SetKink(Bool_t(seed->GetKinkIndex(0)));
-    //Int_t   lbl         = seed->GetLabel();
-    //track.SetSeedLabel(lbl);
+    if(AliTRDgeometry::GetXtrdBeg() > (fgkMaxStep + track.GetX()) && !PropagateToX(track, AliTRDgeometry::GetXtrdBeg(), fgkMaxStep)){ 
+      seed->UpdateTrackParams(&track, AliESDtrack::kTRDStop);
+      continue;
+    }    
+    if(!AdjustSector(&track)){
+      seed->UpdateTrackParams(&track, AliESDtrack::kTRDStop);
+      continue;
+    }
+    if(TMath::Abs(track.GetSnp()) > fgkMaxSnp) {
+      seed->UpdateTrackParams(&track, AliESDtrack::kTRDStop);
+      continue;
+    }
+
+    nTPCseeds++;
 
-    // Make backup and mark entrance in the TRD
-    seed->UpdateTrackParams(&track, AliESDtrack::kTRDin);
+    // store track status at TRD entrance
     seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup);
-    Float_t p4  = track.GetC(track.GetBz());
-    expectedClr = FollowBackProlongation(track);
 
-    if (expectedClr<0) continue; // Back prolongation failed
+    // prepare track and do propagation in the TRD
+    track.SetReconstructor(fReconstructor);
+    track.SetKink(Bool_t(seed->GetKinkIndex(0)));
+    expectedClr = FollowBackProlongation(track);
+    // check if track entered the TRD fiducial volume
+    if(track.GetTrackLow()){ 
+      seed->UpdateTrackParams(&track, AliESDtrack::kTRDin);
+      nTRDseeds++;
+    }
+    // check if track was stopped in the TRD
+    if (expectedClr<0){      
+      seed->UpdateTrackParams(&track, AliESDtrack::kTRDStop);
+      continue;
+    }
 
     if(expectedClr){
-      found++;  
+      nFound++;  
       // computes PID for track
       track.CookPID();
       // update calibration references using this track
       if(calibra->GetHisto2d()) calibra->UpdateHistogramsV1(&track);
       // save calibration object
-      if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0 /*&& quality TODO*/){ 
+      if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0){ 
         AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track);
         calibTrack->SetOwner();
         seed->AddCalibObject(calibTrack);
@@ -349,63 +399,46 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event)
             isGold = kTRUE;
           }
         }
-  
-        //if ((track->StatusForTOF() > 0) && (track->GetNCross() == 0) && (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected())  > 0.4)) {
-        //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
-        //}
       }
     }
     
-    // Propagation to the TOF (I.Belikov)
-    if (track.IsStopped() == kFALSE) {
-      Double_t xtof  = 371.0;
-      Double_t xTOF0 = 370.0;
-    
-      Double_t c2    = track.GetSnp() + track.GetC(track.GetBz()) * (xtof - track.GetX());
-      if (TMath::Abs(c2) >= 0.99) continue;
-      
-      if (!PropagateToX(track, xTOF0, fgkMaxStep)) continue;
-  
-      // Energy losses taken to the account - check one more time
-      c2 = track.GetSnp() + track.GetC(track.GetBz()) * (xtof - track.GetX());
-      if (TMath::Abs(c2) >= 0.99) continue;
-      
-      //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) {
-      //       fHBackfit->Fill(7);
-      //delete track;
-      //       continue;
-      //}
-  
-      Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha());
-      Double_t y;
-      track.GetYAt(xtof,GetBz(),y);
-      if (y >  ymax) {
-        if (!track.Rotate( AliTRDgeometry::GetAlpha())) continue;      
-      }else if (y < -ymax) {
-        if (!track.Rotate(-AliTRDgeometry::GetAlpha())) continue;
-      }
-          
-      if (track.PropagateTo(xtof)) {
-        seed->UpdateTrackParams(&track, AliESDtrack::kTRDout);
-        track.UpdateESDtrack(seed);
-      }
-    } else {                   
-      if ((track.GetNumberOfClusters() > 15) && (track.GetNumberOfClusters() > 0.5*expectedClr)) {
-        seed->UpdateTrackParams(&track, AliESDtrack::kTRDout);
-  
-        track.UpdateESDtrack(seed);
+    // Propagation to the TOF
+    if(!(seed->GetStatus()&AliESDtrack::kTRDStop)) {
+      Int_t sm = track.GetSector();
+      // default value in case we have problems with the geometry.
+      Double_t xtof  = 371.; 
+      //Calculate radial position of the beginning of the TOF
+      //mother volume. In order to avoid mixing of the TRD 
+      //and TOF modules some hard values are needed. This are:
+      //1. The path to the TOF module.
+      //2. The width of the TOF (29.05 cm)
+      //(with the help of Annalisa de Caro Mar-17-2009)
+      if(gGeoManager){
+        gGeoManager->cd(Form("/ALIC_1/B077_1/BSEGMO%d_1/BTOF%d_1", sm, sm));
+        TGeoHMatrix *m = 0x0;
+        Double_t loc[]={0., 0., -.5*29.05}, glob[3];
+        
+        if((m=gGeoManager->GetCurrentMatrix())){
+          m->LocalToMaster(loc, glob);
+          xtof = TMath::Sqrt(glob[0]*glob[0]+glob[1]*glob[1]);
+        }
       }
+      if(xtof > (fgkMaxStep + track.GetX()) && !PropagateToX(track, xtof, fgkMaxStep)) continue;
+      if(!AdjustSector(&track)) continue;
+      if(TMath::Abs(track.GetSnp()) > fgkMaxSnp) continue;
+
+      seed->UpdateTrackParams(&track, AliESDtrack::kTRDout);
+      // TODO obsolete - delete
+      seed->SetTRDQuality(track.StatusForTOF()); 
     }
-  
-    seed->SetTRDQuality(track.StatusForTOF());
     seed->SetTRDBudget(track.GetBudget(0));
   }
   if(index) delete [] index;
   if(quality) delete [] quality;
   
 
-  AliInfo(Form("Number of TPC seeds: %d", nSeed));
-  AliInfo(Form("Number of back propagated TRD tracks: %d", found));
+  AliInfo(Form("Number of TPC seeds: %d (%d)", nTRDseeds, nTPCseeds));
+  AliInfo(Form("Number of propagated TRD tracks: %d", nFound));
       
   // run stand alone tracking
   if (fReconstructor->IsSeeding()) Clusters2Tracks(event);
@@ -576,10 +609,12 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
 
     Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
     TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker);
+    AliTRDtrackV1 track(t);
+    track.SetOwner();
     cstreamer << "FollowProlongation"
         << "EventNumber="      << eventNumber
         << "ncl="                                      << nClustersExpected
-        //<< "track.="                 << &t
+        << "track.="                   << &track
         << "\n";
   }
 
@@ -614,8 +649,8 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
   // Debug level 2
   //
 
-  Int_t nClustersExpected = 0;
-  Double_t clength = .5*AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick();
+  Int_t n = 0;
+  Double_t driftLength = .5*AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick();
   AliTRDtrackingChamber *chamber = 0x0;
   
   AliTRDseedV1 tracklet, *ptrTracklet = 0x0;
@@ -628,136 +663,201 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
   } 
   Bool_t kStoreIn = kTRUE;
 
-
   // Loop through the TRD layers
-  for (Int_t ilayer = 0; ilayer < kNPlanes; ilayer++) {
-    // BUILD TRACKLET IF NOT ALREADY BUILT
-    Double_t x = 0., x0, y, z, alpha;
-    ptrTracklet  = tracklets[ilayer];
-    if(!ptrTracklet){
-      ptrTracklet = new(&tracklet) AliTRDseedV1(ilayer);
-      ptrTracklet->SetReconstructor(fReconstructor);
-      ptrTracklet->SetKink(t.IsKink());
-      alpha = t.GetAlpha();
-      Int_t sector = Int_t(alpha/AliTRDgeometry::GetAlpha() + (alpha>0. ? 0 : AliTRDgeometry::kNsector));
+  TGeoHMatrix *matrix = 0x0;
+  Double_t x, y, z;
+  for (Int_t ily=0, sm=-1, stk=-1, det=-1; ily < AliTRDgeometry::kNlayer; ily++) {
+    // rough estimate of the entry point
+    if (!t.GetProlongation(fR[ily], y, z)){
+      n=-1; 
+      t.SetStatus(AliTRDtrackV1::kProlongation);
+      break;
+    }
 
-      if(!fTrSec[sector].GetNChambers()) continue;
-      
-      if((x = fTrSec[sector].GetX(ilayer)) < 1.) continue;
-    
-      // Propagate closer to the current layer
-      x0 = x - 1.5*clength;
-      if (x0 > (fgkMaxStep + t.GetX()) && !PropagateToX(t, x0-fgkMaxStep, fgkMaxStep)) return -1/*nClustersExpected*/;
-      if (!AdjustSector(&t)) return -1/*nClustersExpected*/;
-      if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) return -1/*nClustersExpected*/;
-
-      if (!t.GetProlongation(x, y, z)) return -1/*nClustersExpected*/;
-      Int_t stack = fGeom->GetStack(z, ilayer);
-      Int_t nCandidates = stack >= 0 ? 1 : 2;
-      z -= stack >= 0 ? 0. : 4.; 
-      
-      for(int icham=0; icham<nCandidates; icham++, z+=8){
-        if((stack = fGeom->GetStack(z, ilayer)) < 0) continue;
-      
-        if(!(chamber = fTrSec[sector].GetChamber(stack, ilayer))) continue;
-      
-        if(chamber->GetNClusters() < fgNTimeBins*fReconstructor->GetRecoParam() ->GetFindableClusters()) continue;
-      
-        x = chamber->GetX();
-      
-        AliTRDpadPlane *pp = fGeom->GetPadPlane(ilayer, stack);
-        tracklet.SetTilt(TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle()));
-        tracklet.SetPadLength(pp->GetLengthIPad());
-        tracklet.SetPadWidth(pp->GetWidthIPad());
-        tracklet.SetDetector(chamber->GetDetector());
-        tracklet.SetX0(x);
-        tracklet.UpDate(&t);
-//         if(!tracklet.Init(&t)){
-//           t.SetStopped(kTRUE);
-//           return nClustersExpected;
-//         }
-        if(!tracklet.AttachClusters(chamber, kTRUE)) continue;
-        //if(!tracklet.AttachClustersIter(chamber, 1000.)) continue;
-        //tracklet.Init(&t);
-        
-        if(tracklet.GetN() < fgNTimeBins*fReconstructor->GetRecoParam() ->GetFindableClusters()) continue;
-      
+    // find sector / stack / detector
+    sm = t.GetSector();
+    // TODO cross check with y value !
+    stk = fGeom->GetStack(z, ily);
+    det = stk>=0 ? AliTRDgeometry::GetDetector(ily, stk, sm) : -1;
+    matrix = det>=0 ? fGeom->GetClusterMatrix(det) : 0x0;
+
+    // check if supermodule/chamber is installed
+    if( !fGeom->GetSMstatus(sm) ||
+        stk<0. ||
+        fGeom->IsHole(ily, stk, sm) ||
+        !matrix ){ 
+      // propagate to the default radial position
+      if(fR[ily] > (fgkMaxStep + t.GetX()) && !PropagateToX(t, fR[ily], fgkMaxStep)){
+        n=-1; 
+        t.SetStatus(AliTRDtrackV1::kPropagation);
         break;
       }
-      ptrTracklet->UpdateUsed();
-    }
-    if(!ptrTracklet->IsOK()){
-      ptrTracklet->Reset();
-      if(x < 1.) continue; //temporary
-      if(!PropagateToX(t, x-fgkMaxStep, fgkMaxStep)) return -1/*nClustersExpected*/;
-      if(!AdjustSector(&t)) return -1/*nClustersExpected*/;
-      if(TMath::Abs(t.GetSnp()) > fgkMaxSnp) return -1/*nClustersExpected*/;
+      if(!AdjustSector(&t)){
+        n=-1; 
+        t.SetStatus(AliTRDtrackV1::kAdjustSector);
+        break;
+      }
+      if(TMath::Abs(t.GetSnp()) > fgkMaxSnp){
+        n=-1; 
+        t.SetStatus(AliTRDtrackV1::kSnp);
+        break;
+      }
+      t.SetStatus(AliTRDtrackV1::kGeometry, ily);
       continue;
     }
-    
-    // Propagate closer to the current chamber if neccessary 
-    x -= clength;
-    if (x > (fgkMaxStep + t.GetX()) && !PropagateToX(t, x-fgkMaxStep, fgkMaxStep)) return -1/*nClustersExpected*/;
-    if (!AdjustSector(&t)) return -1/*nClustersExpected*/;
-    if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) return -1/*nClustersExpected*/;
-    
-    // load tracklet to the tracker and the track
-    ptrTracklet->UseClusters();
-    ptrTracklet->Fit(kFALSE); // no tilt correction
-    ptrTracklet = SetTracklet(ptrTracklet);
-    t.SetTracklet(ptrTracklet, fTracklets->GetEntriesFast()-1);
-  
-  
-    // 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 = ptrTracklet->GetX(); //GetX0();
-    if (!t.GetProlongation(x, y, z)) return -1/*nClustersExpected*/;
-    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];
-    if(AliTracker::MeanMaterialBudget(xyz0, xyz1, param)<=0.) return -1;       
-    // The mean propagation parameters
-    Double_t xrho = param[0]*param[4]; // density*length
-    Double_t xx0  = param[1]; // radiation length
-    
-    // Propagate and update track
-    if (!t.PropagateTo(x, xx0, xrho)) return -1/*nClustersExpected*/;
-    if (!AdjustSector(&t)) return -1/*nClustersExpected*/;
 
+    // retrieve rotation matrix for the current chamber
+    Double_t loc[] = {AliTRDgeometry::AnodePos()- driftLength, 0., 0.};
+    Double_t glb[] = {0., 0., 0.};
+    matrix->LocalToMaster(loc, glb);
+
+    // Propagate to the radial distance of the current layer
+    x = glb[0] - fgkMaxStep;
+    if(x > (fgkMaxStep + t.GetX()) && !PropagateToX(t, x, fgkMaxStep)){
+      n=-1; 
+      t.SetStatus(AliTRDtrackV1::kPropagation);
+      break;
+    }
+    if(!AdjustSector(&t)){
+      n=-1; 
+      t.SetStatus(AliTRDtrackV1::kAdjustSector);
+      break;
+    }
+    if(TMath::Abs(t.GetSnp()) > fgkMaxSnp) {
+      n=-1; 
+      t.SetStatus(AliTRDtrackV1::kSnp);
+      break;
+    }
+    Bool_t RECALCULATE = kFALSE;
+    if(sm != t.GetSector()){
+      sm = t.GetSector(); 
+      RECALCULATE = kTRUE;
+    }
+    if(stk != fGeom->GetStack(z, ily)){
+      stk = fGeom->GetStack(z, ily);
+      RECALCULATE = kTRUE;
+    }
+    if(RECALCULATE){
+      det = AliTRDgeometry::GetDetector(ily, stk, sm);
+      if(!(matrix = fGeom->GetClusterMatrix(det))){ 
+        t.SetStatus(AliTRDtrackV1::kGeometry, ily);
+        continue;
+      }
+      matrix->LocalToMaster(loc, glb);
+      x = glb[0] - fgkMaxStep;
+    }
+
+    // check if track is well inside fiducial volume 
+    if (!t.GetProlongation(x+fgkMaxStep, y, z)) {
+      n=-1; 
+      t.SetStatus(AliTRDtrackV1::kProlongation);
+      break;
+    }
+    if(fGeom->IsOnBoundary(det, y, z, .5)){ 
+      t.SetStatus(AliTRDtrackV1::kBoundary, ily);
+      continue;
+    }
+    // mark track as entering the FIDUCIAL volume of TRD
     if(kStoreIn){
       t.SetTrackLow(); 
       kStoreIn = kFALSE;
     }
-    Double_t maxChi2 = t.GetPredictedChi2(ptrTracklet);
-    if (!t.Update(ptrTracklet, maxChi2)) return -1/*nClustersExpected*/;
-    ptrTracklet->UpDate(&t);
 
-    if (maxChi2<1e+10) { 
-      nClustersExpected += ptrTracklet->GetN();
-      //t.SetTracklet(&tracklet, index);
+    ptrTracklet  = tracklets[ily];
+    if(!ptrTracklet){ // BUILD TRACKLET
+      // check data in supermodule
+      if(!fTrSec[sm].GetNChambers()){ 
+        t.SetStatus(AliTRDtrackV1::kNoClusters, ily);
+        continue;
+      }
+      if(fTrSec[sm].GetX(ily) < 1.){ 
+        t.SetStatus(AliTRDtrackV1::kNoClusters, ily);
+        continue;
+      }
+      
+      // check data in chamber
+      if(!(chamber = fTrSec[sm].GetChamber(stk, ily))){ 
+        t.SetStatus(AliTRDtrackV1::kNoClusters, ily);
+        continue;
+      }
+      if(chamber->GetNClusters() < fgNTimeBins*fReconstructor->GetRecoParam() ->GetFindableClusters()){ 
+        t.SetStatus(AliTRDtrackV1::kNoClusters, ily);
+        continue;
+      }      
+      // build tracklet
+      ptrTracklet = new(&tracklet) AliTRDseedV1(det);
+      ptrTracklet->SetReconstructor(fReconstructor);
+      ptrTracklet->SetKink(t.IsKink());
+      ptrTracklet->SetPadPlane(fGeom->GetPadPlane(ily, stk));
+      ptrTracklet->SetX0(glb[0]+driftLength);
+      if(!tracklet.Init(&t)){
+        n=-1; 
+        t.SetStatus(AliTRDtrackV1::kTrackletInit);
+        break;
+      }
+      if(!tracklet.AttachClusters(chamber, kTRUE)){   
+        t.SetStatus(AliTRDtrackV1::kNoAttach, ily);
+        continue;
+      }
+      if(tracklet.GetN() < fgNTimeBins*fReconstructor->GetRecoParam() ->GetFindableClusters()){
+        t.SetStatus(AliTRDtrackV1::kNoClustersTracklet, ily);
+        continue;
+      }
+      ptrTracklet->UpdateUsed();
+    }
+
+    // propagate track to the radial position of the tracklet
+    ptrTracklet->UseClusters(); // TODO ? do we need this here ?
+    // fit tracklet no tilt correction
+    if(!ptrTracklet->Fit(kFALSE)){
+      t.SetStatus(AliTRDtrackV1::kNoFit, ily);
+      continue;
+    } 
+    x = ptrTracklet->GetX(); //GetX0();
+    if(x > (fgkMaxStep + t.GetX()) && !PropagateToX(t, x, fgkMaxStep)) {
+      n=-1; 
+      t.SetStatus(AliTRDtrackV1::kPropagation);
+      break;
+    }
+    if(!AdjustSector(&t)) {
+      n=-1; 
+      t.SetStatus(AliTRDtrackV1::kAdjustSector);
+      break;
     }
+    if(TMath::Abs(t.GetSnp()) > fgkMaxSnp) {
+      n=-1; 
+      t.SetStatus(AliTRDtrackV1::kSnp);
+      break;
+    }
+  
+    // update Kalman with the TRD measurement
+    Double_t chi2 = t.GetPredictedChi2(ptrTracklet);
+    if(chi2>1e+10){
+      t.SetStatus(AliTRDtrackV1::kChi2, ily);
+      continue; 
+    }
+    if(!t.Update(ptrTracklet, chi2)) {
+      n=-1; 
+      t.SetStatus(AliTRDtrackV1::kUpdate);
+      break;
+    }
+
+    // load tracklet to the tracker
+    ptrTracklet->UpDate(&t);
+    ptrTracklet = SetTracklet(ptrTracklet);
+    t.SetTracklet(ptrTracklet, fTracklets->GetEntriesFast()-1);
+    n += ptrTracklet->GetN();
+
     // Reset material budget if 2 consecutive gold
-    if(ilayer>0 && t.GetTracklet(ilayer-1) && ptrTracklet->GetN() + t.GetTracklet(ilayer-1)->GetN() > 20) t.SetBudget(2, 0.);
+//     if(ilayer>0 && t.GetTracklet(ilayer-1) && ptrTracklet->GetN() + t.GetTracklet(ilayer-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 = ptrTracklet->GetN() / Float_t(fgNTimeBins);
     //Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1);       
-    //printf("tracklet.GetChi2() %f     [< 18.0]\n", tracklet.GetChi2()); 
-    //printf("ratio0    %f              [>   0.8]\n", ratio0);
-    //printf("ratio1     %f             [>   0.6]\n", ratio1); 
-    //printf("ratio0+ratio1 %f          [>   1.5]\n", ratio0+ratio1); 
-    //printf("t.GetNCross()  %d         [==    0]\n", t.GetNCross()); 
-    //printf("TMath::Abs(t.GetSnp()) %f [<  0.85]\n", TMath::Abs(t.GetSnp()));
-    //printf("t.GetNumberOfClusters() %d [>    20]\n", t.GetNumberOfClusters());
     
-    if (//(tracklet.GetChi2()      <  18.0) && TO DO check with FindClusters and move it to AliTRDseed::Update 
+    if( (chi2                    <  18.0) &&  
         (ratio0                  >   0.8) && 
         //(ratio1                  >   0.6) && 
         //(ratio0+ratio1           >   1.5) && 
@@ -767,20 +867,22 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
       t.MakeBackupTrack();
     }
   } // end layers loop
+  //printf("clusters[%d] chi2[%f] x[%f] status[%d ", n, t.GetChi2(), t.GetX(), t.GetStatusTRD());
+  //for(int i=0; i<6; i++) printf("%d ", t.GetStatusTRD(i)); printf("]\n");
 
   if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){
     TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker);
     Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
-    //AliTRDtrackV1 *debugTrack = new AliTRDtrackV1(t);
-    //debugTrack->SetOwner();
+    AliTRDtrackV1 track(t);
+    track.SetOwner();
     cstreamer << "FollowBackProlongation"
-        << "EventNumber="                      << eventNumber
-        << "ncl="                                                      << nClustersExpected
-        //<< "track.="                                 << debugTrack
+        << "EventNumber=" << eventNumber
+        << "ncl="         << n
+        << "track.="      << &track
         << "\n";
   }
   
-  return nClustersExpected;
+  return n;
 }
 
 //_________________________________________________________________________
@@ -1913,7 +2015,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon
 
   // Build initial seeding configurations
   Double_t quality = BuildSeedingConfigs(stack, configs);
-  if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){
+  if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 10){
     AliInfo(Form("Plane config %d %d %d Quality %f"
     , configs[0], configs[1], configs[2], quality));
   }
@@ -1941,7 +2043,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon
       ntracks = MakeSeeds(stack, &sseed[6*ntracks], pars);
       if(ntracks == kMaxTracksStack) break;
     }
-    if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1) AliInfo(Form("Candidate TRD tracks %d in iteration %d.", ntracks, fSieveSeeding));
+    if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 10) AliInfo(Form("Candidate TRD tracks %d in iteration %d.", ntracks, fSieveSeeding));
     
     if(!ntracks) break;
     
@@ -2053,7 +2155,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon
         trackParams[6] = fGeom->GetSector(chamber->GetDetector());/* *alpha+shift;     // Supermodule*/
 
         if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){
-          AliInfo(Form("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1]));
+          //AliInfo(Form("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1]));
 
           AliTRDseedV1 *dseed[6];
           for(Int_t iseed = AliTRDgeometry::kNlayer; iseed--;) dseed[iseed] = new AliTRDseedV1(lseed[iseed]);
@@ -2130,7 +2232,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon
       chamber->Build(fGeom, cal);//Indices(fSieveSeeding);
     }
 
-    if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){ 
+    if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 10){ 
       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
@@ -2169,15 +2271,15 @@ Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDtrackingChamber **stack, Int
   // The overall chamber quality is given by the product of this 2 contributions.
   // 
 
-  Double_t chamberQ[kNPlanes];
+  Double_t chamberQ[kNPlanes];memset(chamberQ, 0, kNPlanes*sizeof(Double_t));
   AliTRDtrackingChamber *chamber = 0x0;
   for(int iplane=0; iplane<kNPlanes; iplane++){
     if(!(chamber = stack[iplane])) continue;
     chamberQ[iplane] = (chamber = stack[iplane]) ?  chamber->GetQuality() : 0.;
   }
 
-  Double_t tconfig[kNConfigs];
-  Int_t planes[4];
+  Double_t tconfig[kNConfigs];memset(tconfig, 0, kNConfigs*sizeof(Double_t));
+  Int_t planes[] = {0, 0, 0, 0};
   for(int iconf=0; iconf<kNConfigs; iconf++){
     GetSeedingConfig(iconf, planes);
     tconfig[iconf] = fgTopologicQA[iconf];
@@ -2254,10 +2356,6 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss
   // chi2[1] = tracklet chi2 on the R direction
   Double_t chi2[4];
 
-       // Default positions for the anode wire in all 6 Layers in case of a stack with missing clusters
-       // Positions taken using cosmic data taken with SM3 after rebuild
-  Double_t x_def[kNPlanes] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2};
-
   // this should be data member of AliTRDtrack
   Double_t seedQuality[kMaxTracksStack];
   
@@ -2281,7 +2379,21 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss
     padwidth[iplane] = pp->GetWidthIPad();
   }
   
-  if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){
+  // Init anode wire position for chambers
+  Double_t x0[kNPlanes] = {0., 0., 0., 0., 0., 0.},       // anode wire position
+           driftLength = .5*AliTRDgeometry::AmThick() - AliTRDgeometry::DrThick(); // drift length
+  TGeoHMatrix *matrix = 0x0;
+  Double_t loc[] = {AliTRDgeometry::AnodePos(), 0., 0.};
+  Double_t glb[] = {0., 0., 0.};
+  AliTRDtrackingChamber **cIter = &stack[0];
+  for(int iLayer=kNPlanes; iLayer--; cIter++){
+    if(!(*cIter)) continue;
+    if(!(matrix = fGeom->GetClusterMatrix((*cIter)->GetDetector()))) continue;
+    matrix->LocalToMaster(loc, glb);
+    x0[iLayer] = glb[0];
+  }
+
+  if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 10){
     AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks));
   }
 
@@ -2293,7 +2405,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss
     if(!chamber->GetSeedingLayer(fSeedTB[isl], fGeom, fReconstructor)) continue;
     nlayers++;
   }
-  if(nlayers < 4) return ntracks;
+  if(nlayers < kNSeedPlanes) return ntracks;
   
   
   // Start finding seeds
@@ -2334,12 +2446,13 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss
         AliTRDseedV1 *tseed = &cseed[0];
         AliTRDtrackingChamber **cIter = &stack[0];
         for(int iLayer=0; iLayer<kNPlanes; iLayer++, tseed++, cIter++){
-          tseed->SetDetector((*cIter) ? (*cIter)->GetDetector() : -1);
+          Int_t det = (*cIter) ? (*cIter)->GetDetector() : -1;
+          tseed->SetDetector(det);
           tseed->SetTilt(hL[iLayer]);
           tseed->SetPadLength(padlength[iLayer]);
           tseed->SetPadWidth(padwidth[iLayer]);
           tseed->SetReconstructor(fReconstructor);
-          tseed->SetX0((*cIter) ? (*cIter)->GetX() : x_def[iLayer]);
+          tseed->SetX0(det<0 ? fR[iLayer]+driftLength : x0[iLayer]);
           tseed->Init(GetRiemanFitter());
           tseed->SetStandAlone(kTRUE);
         }
index 49a9e87..6acad41 100644 (file)
@@ -3,7 +3,6 @@
 
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */ 
-
 /* $Id$ */
 
 ////////////////////////////////////////////////////////////////////////////
@@ -73,10 +72,10 @@ public:
   AliCluster*     GetCluster(Int_t index) const;
   AliTRDseedV1*   GetTracklet(Int_t index) const;
   AliKalmanTrack* GetTrack(Int_t index) const;
-  TClonesArray*   GetListOfClusters() const {return fClusters;}
-  TClonesArray*   GetListOfTracklets() const {return fTracklets;}
-  TClonesArray*   GetListOfTracks() const {return fTracks;}
-  static Int_t     GetNTimeBins() {return fgNTimeBins;}
+  TClonesArray*   GetListOfClusters() const  { return fClusters;}
+  TClonesArray*   GetListOfTracklets() const { return fTracklets;}
+  TClonesArray*   GetListOfTracks() const    { return fTracks;}
+  static Int_t    GetNTimeBins()             { return fgNTimeBins;}
   static void     GetExtrapolationConfig(Int_t iconfig, Int_t planes[2]);
   static void     GetSeedingConfig(Int_t iconfig, Int_t planes[4]);
   static TLinearFitter*  GetTiltedRiemanFitter();
@@ -91,7 +90,7 @@ public:
        static Double_t FitLine(const AliTRDtrackV1 *trk, AliTRDseedV1 *tracklets = 0x0, Bool_t err=0, Int_t np = 0, AliTrackPoint *points = 0x0);
        static Double_t FitKalman(AliTRDtrackV1 *trk, AliTRDseedV1 *tracklets = 0x0, Bool_t up=0, Int_t np = 0, AliTrackPoint *points = 0x0);
 
-  Bool_t          IsClustersOwner() const {return TestBit(kOwner);}
+  Bool_t          IsClustersOwner() const    { return TestBit(kOwner);}
   void            SetClustersOwner(Bool_t own=kTRUE) {SetBit(kOwner, own); if(!own) fClusters = 0x0;}
 
   Int_t           FollowBackProlongation(AliTRDtrackV1 &t);
@@ -106,7 +105,6 @@ public:
   void            SetReconstructor(const AliTRDReconstructor *rec){ fReconstructor = rec; }
   void            UnloadClusters();
 //   void            UseClusters(const AliKalmanTrack *t, Int_t from = 0) const;
-//   static Int_t    Freq(Int_t n, const Int_t *inlist, Int_t *outlist, Bool_t down); // to be removed 
 
   class AliTRDLeastSquare{
   public:
@@ -137,6 +135,7 @@ protected:
   Int_t          Clusters2TracksStack(AliTRDtrackingChamber **stack, TClonesArray *esdTrackList);
   AliTRDseedV1*  GetTracklet(AliTRDtrackV1 *trk, Int_t plane, Int_t &idx);
   Bool_t         GetTrackPoint(Int_t index, AliTrackPoint &p) const;   
+  Float_t        GetR4Layer(Int_t ly) const { return fR[ly];}
   Int_t          MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *sseed, Int_t *ipar);
   AliTRDtrackV1* MakeTrack(AliTRDseedV1 *seeds, Double_t *params);
   AliTRDtrackV1* SetTrack(AliTRDtrackV1 *track);
@@ -154,20 +153,20 @@ private:
   Float_t     GetChi2Z(AliTRDseedV1 *tracklets) const;
 
 private:
-  const AliTRDReconstructor *fReconstructor;                 // reconstructor manager
-  AliTRDgeometry      *fGeom;                          // Pointer to TRD geometry
-  AliTRDtrackingSector fTrSec[kTrackingSectors];       // Array of tracking sectors;    
-
+  const AliTRDReconstructor *fReconstructor;            // reconstructor manager
+  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
-  
+  Float_t              fR[kNPlanes];                    //! rough radial position of each TRD layer
+
   // 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  
+  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
@@ -181,7 +180,7 @@ private:
        static TLinearFitter *fgTiltedRiemanConstrained;      //  Fitter for the tilted Rieman fit with vertex constraint       
        static AliRieman     *fgRieman;                       //  Fitter for the untilted Rieman fit
        
-       ClassDef(AliTRDtrackerV1, 2)                          //  TRD tracker development class
+       ClassDef(AliTRDtrackerV1, 3)                          //  TRD tracker - tracklet based tracking
 
 };
 #endif
index 03625d1..f2500b5 100644 (file)
@@ -277,7 +277,7 @@ TH1 *AliTRDpidChecker::PlotLQ(const AliTRDtrackV1 *track)
 
   fReconstructor -> SetOption("!nn");
   cTrack.CookPID();
-  if(cTrack.GetPIDquality() < fMinNTracklets) return 0x0;
+  if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return 0x0;
 
   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
   hPIDLQ -> Fill(FindBin(species, momentum), cTrack.GetPID(AliPID::kElectron));
@@ -336,7 +336,7 @@ TH1 *AliTRDpidChecker::PlotNN(const AliTRDtrackV1 *track)
 
   fReconstructor -> SetOption("nn");
   cTrack.CookPID();
-  if(cTrack.GetPIDquality() < fMinNTracklets) return 0x0;
+  if(cTrack.GetNumberOfTrackletsPID() < fMinNTracklets) return 0x0;
 
   Int_t species = AliTRDpidUtil::Pdg2Pid(pdg);
   hPIDNN -> Fill(FindBin(species, momentum), cTrack.GetPID(AliPID::kElectron));
@@ -367,7 +367,7 @@ TH1 *AliTRDpidChecker::PlotESD(const AliTRDtrackV1 *track)
   if(!(status&AliESDtrack::kTRDin)) return 0x0;
 
   if(!CheckTrackQuality(fTrack)) return 0x0;
-  if(fTrack->GetPIDquality() < fMinNTracklets) return 0x0;
+  if(fTrack->GetNumberOfTrackletsPID() < fMinNTracklets) return 0x0;
   
   if(!(fEfficiency = dynamic_cast<TObjArray *>(fContainer->At(kEfficiency)))){
     AliWarning("No Histogram defined.");