Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of...
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 May 2001 17:07:58 +0000 (17:07 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 May 2001 17:07:58 +0000 (17:07 +0000)
TRD/AliTRD.cxx
TRD/AliTRD.h
TRD/AliTRDclusterizerV0.cxx
TRD/AliTRDclusterizerV1.cxx
TRD/AliTRDdigitizer.h
TRD/AliTRDtrack.cxx
TRD/AliTRDtrack.h
TRD/AliTRDtracker.cxx
TRD/AliTRDtracker.h
TRD/AliTRDtrackingSector.cxx
TRD/AliTRDtrackingSector.h

index efa195a..fd840b6 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.29  2001/05/21 16:45:47  hristov
+Last minute changes (C.Blume)
+
 Revision 1.28  2001/05/16 14:57:27  alibrary
 New files for folders and Stack
 
@@ -267,7 +270,7 @@ AliTRD::~AliTRD()
 
 //_____________________________________________________________________________
 void AliTRD::AddCluster(Float_t *pos, Int_t *digits, Int_t det, Float_t amp
-                       , Int_t *tracks, Int_t iType)
+                       , Int_t *tracks, Float_t sigmaY2, Int_t iType)
 {
   //
   // Add a cluster for the TRD
@@ -295,7 +298,7 @@ void AliTRD::AddCluster(Float_t *pos, Int_t *digits, Int_t det, Float_t amp
   c->SetY(- (col0 + padCol * colSize));
   c->SetZ(   row0 + padRow * rowSize);
   
-  c->SetSigmaY2(0.05 * 0.05);
+  c->SetSigmaY2((sigmaY2 + 1./12.) * colSize*colSize);   
   c->SetSigmaZ2(rowSize * rowSize / 12.);
 
   switch (iType) {
index 519cf43..36c1f6f 100644 (file)
@@ -36,7 +36,8 @@ class AliTRD : public AliDetector {
   virtual void       AddHit(Int_t track, Int_t det, Float_t *hits, Int_t q);
   virtual void       AddDigit(Int_t *digits, Int_t *amp);    
   virtual void       AddCluster(Float_t *pos, Int_t *digits
-                              , Int_t det, Float_t amp, Int_t *tracks, Int_t iType);
+                              , Int_t det, Float_t amp, Int_t *tracks
+                              , Float_t sigmaY2, Int_t iType);
   virtual void       BuildGeometry();
   virtual void       Copy(TObject &trd);
   virtual void       CreateGeometry();
index 535bf4a..19f6a70 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.8  2001/05/07 08:06:44  cblume
+Speedup of the code. Create only AliTRDcluster
+
 Revision 1.7  2001/02/14 18:22:26  cblume
 Change in the geometry of the padplane
 
@@ -317,7 +320,7 @@ Bool_t AliTRDclusterizerV0::MakeClusters()
           Int_t detector  = recPoint1->GetDetector();
           Int_t digits[3] = {0};
          Int_t tr[9] = {-1}; 
-          fTRD->AddCluster(smear,digits,detector,0.0,tr,0);
+          fTRD->AddCluster(smear,digits,detector,0.0,tr,0,0);
 
        }
 
index d05df73..6c79379 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.12  2001/05/21 17:42:58  hristov
+Constant casted to avoid the ambiguity
+
 Revision 1.11  2001/05/21 16:45:47  hristov
 Last minute changes (C.Blume)
 
@@ -71,6 +74,7 @@ Add new TRD classes
 #include <TF1.h>
 #include <TTree.h>
 #include <TH1.h>
+#include <TFile.h>
 
 #include "AliRun.h"
 
@@ -245,6 +249,17 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
   // Get the geometry
   AliTRDgeometry *geo = fTRD->GetGeometry();
 
+  Float_t timeBinSize = geo->GetTimeBinSize();
+  // Half of ampl.region
+  const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.; 
+
+  AliTRDdigitizer *digitizer = (AliTRDdigitizer*) fInputFile->Get("digitizer");
+  printf("AliTRDclusterizerV1::MakeCluster -- ");
+  printf("Got digitizer\n");
+  Float_t omegaTau = digitizer->GetOmegaTau();
+  printf("AliTRDclusterizerV1::MakeCluster -- ");
+  printf("OmegaTau = %f \n",omegaTau);
   printf("AliTRDclusterizerV1::MakeCluster -- ");
   printf("Start creating clusters.\n");
 
@@ -347,9 +362,9 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
           for ( col = 2;  col <  nColMax;    col++) {
             for (time = 0; time < nTimeTotal; time++) {
 
-              Int_t signalL = digits->GetDataUnchecked(row,col  ,time);
-              Int_t signalM = digits->GetDataUnchecked(row,col-1,time);
-              Int_t signalR = digits->GetDataUnchecked(row,col-2,time);
+              Int_t signalL = TMath::Abs(digits->GetDataUnchecked(row,col  ,time));
+              Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time));
+              Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
  
              // Look for the maximum
               if (signalM >= maxThresh) {
@@ -521,6 +536,18 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
 
                }
 
+                Float_t clusterSigmaY2 = (clusterSignal[2] + clusterSignal[0]) / clusterCharge 
+                                       - (clusterPads[1]-col-0.5) * (clusterPads[1]-col-0.5);
+
+                // Correct for ExB displacement
+                if (digitizer->GetExB()) { 
+                  Int_t   local_time_bin = (Int_t) clusterPads[2];
+                  Float_t driftLength    = local_time_bin * timeBinSize + kAmWidth;
+                  Float_t colSize        = geo->GetColPadSize(iplan);
+                  Float_t deltaY         = omegaTau*driftLength/colSize;
+                  clusterPads[1]         = clusterPads[1] - deltaY;
+                }
+                                       
                 if (fVerbose) {
                   printf("-----------------------------------------------------------\n");
                   printf("Create cluster no. %d\n",nClusters);
@@ -549,6 +576,7 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
                                 ,idet
                                 ,clusterCharge
                                 ,clusterTracks
+                               ,clusterSigmaY2
                                 ,iType);
 
               }
@@ -617,8 +645,8 @@ Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
                      / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
 
     // Set cluster charge ratio
-    Float_t ampLeft  = padSignal[1];
-    Float_t ampRight = padSignal[3];
+    Float_t ampLeft  = padSignal[1] / PadResponse(0 - maxLeft );
+    Float_t ampRight = padSignal[3] / PadResponse(0 - maxRight);
 
     // Apply pad response to parameters
     newLeftSignal[0]  = ampLeft  * PadResponse(-1 - maxLeft);
index dd825eb..9fdc87b 100644 (file)
@@ -70,6 +70,7 @@ class AliTRDdigitizer : public TNamed {
           Float_t      GetDiffusionT() const            { return fDiffusionT;         };
           Float_t      GetDiffusionL() const            { return fDiffusionL;         };
           Float_t      GetElAttachProp() const          { return fElAttachProp;       };
+          Int_t        GetExB() const                   { return fExBOn;              };
           Float_t      GetOmegaTau() const              { return fOmegaTau;           };
           Float_t      GetDriftVelocity() const         { return fDriftVelocity;      };
           Float_t      GetPadCoupling() const           { return fPadCoupling;        };
index a2b5f4c..c9f11d1 100644 (file)
@@ -15,9 +15,6 @@
 
 /*
 $Log$
-Revision 1.5  2001/02/05 14:49:32  hristov
-Compare() declared const (R.Brun)
-
 Revision 1.4  2000/12/08 16:07:02  cblume
 Update of the tracking by Sergei
 
@@ -46,8 +43,8 @@ ClassImp(AliTRDtrack)
 
 //_____________________________________________________________________________
 
-AliTRDtrack::AliTRDtrack(UInt_t index, const Double_t xx[5],
-const Double_t cc[15], Double_t xref, Double_t alpha) {
+AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index, 
+const Double_t xx[5], const Double_t cc[15], Double_t xref, Double_t alpha) {
   //-----------------------------------------------------------------
   // This is the main track constructor.
   //-----------------------------------------------------------------
@@ -67,11 +64,12 @@ const Double_t cc[15], Double_t xref, Double_t alpha) {
   fCty=cc[10]; fCtz=cc[11]; fCtc=cc[12]; fCte=cc[13]; fCtt=cc[14];
 
   fN=0;
-  fIndex[fN++]=index;
-
-  fLhPion     = 0.0;
-  fLhElectron = 0.0;
+  fIndex[fN]=index;
 
+  Float_t q = c->GetQ();
+  Double_t s = fX*fC - fE, t=fT;
+  q *= TMath::Sqrt((1-s*s)/(1+t*t)); 
+  fdQdl[fN++] = q;
 }                              
            
 //_____________________________________________________________________________
@@ -97,11 +95,10 @@ AliTRDtrack::AliTRDtrack(const AliTRDtrack& t) {
   fCty=t.fCty;  fCtz=t.fCtz;  fCtc=t.fCtc;  fCte=t.fCte;  fCtt=t.fCtt;
 
   fN=t.fN;
-  for (Int_t i=0; i<fN; i++) fIndex[i]=t.fIndex[i];
-
-  fLhPion     = t.fLhPion;
-  fLhElectron = t.fLhElectron;
-
+  for (Int_t i=0; i<fN; i++) {
+    fIndex[i]=t.fIndex[i];
+    fdQdl[i]=t.fdQdl[i];
+  }
 }                                                       
 
 //_____________________________________________________________________________
@@ -130,6 +127,37 @@ Int_t AliTRDtrack::Compare(const TObject *o) const {
   return 0;
 }                
 
+//_____________________________________________________________________________
+void AliTRDtrack::CookdEdx(Double_t low, Double_t up) {
+  //-----------------------------------------------------------------
+  // Calculates dE/dX within the "low" and "up" cuts.
+  //-----------------------------------------------------------------
+  Int_t i;
+  Int_t nc=GetNclusters();
+
+  Float_t sorted[200];
+  for (i=0; i<200; i++) sorted[i]=fdQdl[i];
+
+  Int_t swap; 
+  do {
+    swap=0;
+    for (i=0; i<nc-1; i++) {
+      if (sorted[i]<=sorted[i+1]) continue;
+      Float_t tmp=sorted[i];
+      sorted[i]=sorted[i+1]; sorted[i+1]=tmp;
+      swap++;
+    }
+  } while (swap);
+
+  Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
+  Float_t dedx=0;
+  for (i=nl; i<=nu; i++) dedx += sorted[i];
+  dedx /= (nu-nl+1);
+  SetdEdx(dedx);
+}                     
+
+
+
 //_____________________________________________________________________________
 Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
 {
@@ -211,25 +239,12 @@ Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho,Double_t pm)
 }     
 
 
-//_____________________________________________________________________________
-void AliTRDtrack::PropagateToVertex(Double_t x0,Double_t rho,Double_t pm)
-{
-  // This function propagates tracks to the "vertex".
-
-  Double_t c=fC*fX - fE;
-  Double_t tgf=-fE/(fC*fY + sqrt(1-c*c));
-  Double_t snf=tgf/sqrt(1.+ tgf*tgf);
-  Double_t xv=(fE+snf)/fC;
-  PropagateTo(xv,x0,rho,pm); 
-}          
-
-
 //_____________________________________________________________________________
 void AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index)
 {
   // Assignes found cluster to the track and updates track information
 
-  Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2()*12;
+  Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
   r00+=fCyy; r01+=fCzy; r11+=fCzz;
   Double_t det=r00*r11 - r01*r01;
   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
@@ -272,7 +287,13 @@ void AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index)
 
   fCtt-=k40*c04+k41*c14;
 
-  fIndex[fN++]=index;
+  fIndex[fN]=index;
+
+  Float_t q = c->GetQ();
+  Double_t s = fX*fC - fE, t=fT;
+  q *= TMath::Sqrt((1-s*s)/(1+t*t)); 
+  fdQdl[fN++] = q;
+
   fChi2 += chisq;   
 
   //  cerr<<"in update: fIndex["<<fN<<"] = "<<index<<endl;
@@ -339,12 +360,11 @@ Int_t AliTRDtrack::Rotate(Double_t alpha)
 }                         
 
 
-
-
 //_____________________________________________________________________________
 Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c) const
 {
-  Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2()*12;
+  /*
+  Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
   r00+=fCyy; r01+=fCzy; r11+=fCzz;
 
   Double_t det=r00*r11 - r01*r01;
@@ -357,6 +377,13 @@ Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c) const
   Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
 
   return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;  
+  */
+
+  Double_t dy=c->GetY() - fY;
+  Double_t r00=c->GetSigmaY2();
+
+  return (dy*dy)/r00;
+
 }            
 
 
@@ -410,10 +437,7 @@ void AliTRDtrack::Streamer(TBuffer &R__b)
       R__b >> fCtt;
       R__b >> fN;
       for (Int_t i=0; i<fN; i++) R__b >> fIndex[i];
-      if (R__v > 1) {
-        R__b >> fLhElectron;
-        R__b >> fLhPion;
-      }
+      for (Int_t i=0; i<fN; i++) R__b >> fdQdl[i];
    } else {                                
       R__b.WriteVersion(AliTRDtrack::IsA());
       TObject::Streamer(R__b);
@@ -444,8 +468,7 @@ void AliTRDtrack::Streamer(TBuffer &R__b)
       R__b << fCtt;
       R__b << fN;
       for (Int_t i=0; i<fN; i++) R__b << fIndex[i];
-      R__b << fLhElectron;
-      R__b << fLhPion;
+      for (Int_t i=0; i<fN; i++) R__b << fdQdl[i];
    }
 }                                                          
 
index 540bf8f..2dea83a 100644 (file)
@@ -15,18 +15,21 @@ class AliTRDtrack : public TObject {
 public:
 
    AliTRDtrack() { fN=0;}
-   AliTRDtrack(UInt_t index, const Double_t xx[5],
+   AliTRDtrack(const AliTRDcluster *c, UInt_t index, const Double_t xx[5],
                const Double_t cc[15], Double_t xr, Double_t alpha);  
    AliTRDtrack(const AliTRDtrack& t);    
 
-   Int_t  Compare(const TObject *o) const;
+   Int_t    Compare(const TObject *o) const;
+   void     CookdEdx(Double_t low=0.05, Double_t up=0.70);   
 
    Double_t GetAlpha() const {return fAlpha;}
    Double_t GetC()     const {return fC;}
    Int_t    GetClusterIndex(Int_t i) const {return fIndex[i];}    
+   Float_t  GetClusterdQdl(Int_t i) const {return fdQdl[i];}    
    Double_t GetChi2()  const {return fChi2;}
+   Double_t GetZchi2()  const {return fZchi2;}
    void     GetCovariance(Double_t cov[15]) const;  
-   Double_t GetdEdX()  const {return fdEdx;}
+   Double_t GetdEdx()  const {return fdEdx;}
    Double_t GetEta()   const {return fE;}
    Int_t    GetLabel() const {return fLab;}
    Int_t    GetNclusters() const {return fN;}
@@ -34,8 +37,7 @@ public:
      return TMath::Abs(GetPt())*sqrt(1.+GetTgl()*GetTgl());
    }
    Double_t GetPredictedChi2(const AliTRDcluster*) const ;
-   //Double_t GetPt()    const {return 0.3*0.2/GetC()/100;}
-   Double_t GetPt()    const {return 0.3*0.4/GetC()/100;}
+   Double_t GetPt()    const {return 0.299792458*0.4/GetC()/100;}
    void     GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const ;
    Double_t GetSigmaC2()   const {return fCcc;}
    Double_t GetSigmaTgl2() const {return fCtt;}
@@ -45,16 +47,13 @@ public:
    Double_t GetX()   const {return fX;}
    Double_t GetY()   const {return fY;} // returns running Y
    Double_t GetZ()   const {return fZ;}
-
-   Float_t  GetLikelihoodPion()     const { return fLhPion;     };
    Float_t  GetLikelihoodElectron() const { return fLhElectron; };
-
    Bool_t   IsSortable() const {return kTRUE;}
 
    Int_t    PropagateTo(Double_t xr,
                    Double_t x0=8.72,Double_t rho=5.86e-3,Double_t pm=0.139);
-   void     PropagateToVertex(
-                   Double_t x0=36.66,Double_t rho=1.2e-3,Double_t pm=0.139);
    Int_t    Rotate(Double_t angle);
 
    void     SetLabel(Int_t lab=0) {fLab=lab;}
@@ -62,13 +61,13 @@ public:
 
    void     Update(const AliTRDcluster* c, Double_t chi2, UInt_t i);
 
-   void     SetLikelihoodPion(Float_t l)     { fLhPion     = l; };
-   void     SetLikelihoodElectron(Float_t l) { fLhElectron = l; };
+   void     SetLikelihoodElectron(Float_t l) { fLhElectron = l; };   
 
 protected:
 
    Int_t    fLab;         // track label  
-   Double_t fChi2;        // total chi2 value for the track  
+   Double_t fChi2;        // total chi2 value for R*phi measurements
+   Double_t fZchi2;       // total chi2 value for Z measurements  
    Float_t  fdEdx;        // dE/dx 
 
    Double_t fAlpha;       // rotation angle
@@ -88,11 +87,12 @@ protected:
 
    Short_t fN;             // number of clusters associated with the track
    UInt_t  fIndex[200];    // global indexes of these clusters  
+   Float_t fdQdl[200];     // cluster amplitudes corrected for track angles    
+                          
+   Float_t fLhElectron;    // Likelihood to be an electron    
 
-   Float_t fLhElectron;    // Likelihood to be an electron
-   Float_t fLhPion;        // Likelihood to be a pion             
+   ClassDef(AliTRDtrack,2) // TRD reconstructed tracks
 
-   ClassDef(AliTRDtrack,2)  // TRD reconstructed tracks
 };                     
 
 
index 85e686c..c423db0 100644 (file)
                                                       
 /*
 $Log$
-Revision 1.10  2001/05/07 08:08:05  cblume
-Update of TRD code
-
-Revision 1.9  2001/02/14 18:22:26  cblume
-Change in the geometry of the padplane
-
 Revision 1.8  2000/12/20 13:00:44  cblume
 Modifications for the HP-compiler
 
@@ -70,25 +64,28 @@ Add the tracking code
 
 ClassImp(AliTRDtracker) 
 
-  const  Int_t     AliTRDtracker::fSeedGap            = 35;  
-  const  Int_t     AliTRDtracker::fSeedStep           = 5;   
+  const  Float_t     AliTRDtracker::fSeedDepth          = 0.5; 
+  const  Float_t     AliTRDtracker::fSeedStep           = 0.05;   
+  const  Float_t     AliTRDtracker::fSeedGap            = 0.25;  
+
+  const  Float_t     AliTRDtracker::fMaxSeedDeltaZ12    = 40.;  
+  const  Float_t     AliTRDtracker::fMaxSeedDeltaZ      = 25.;  
+  const  Float_t     AliTRDtracker::fMaxSeedC           = 0.0052; 
+  const  Float_t     AliTRDtracker::fMaxSeedTan         = 1.2;  
+  const  Float_t     AliTRDtracker::fMaxSeedVertexZ     = 150.; 
 
+  const  Double_t    AliTRDtracker::fSeedErrorSY        = 0.2;
+  const  Double_t    AliTRDtracker::fSeedErrorSY3       = 2.5;
+  const  Double_t    AliTRDtracker::fSeedErrorSZ        = 0.1;
 
-  const  Float_t   AliTRDtracker::fMinClustersInTrack = 0.5;  
-  const  Float_t   AliTRDtracker::fMinClustersInSeed  = 0.5;  
-  const  Float_t   AliTRDtracker::fSeedDepth          = 0.5; 
-  const  Float_t   AliTRDtracker::fSkipDepth          = 0.2;
-  const  Float_t   AliTRDtracker::fMaxSeedDeltaZ      = 30.;  
-  const  Float_t   AliTRDtracker::fMaxSeedC           = 0.01; 
-  const  Float_t   AliTRDtracker::fMaxSeedTan         = 1.2;  
-  const  Float_t   AliTRDtracker::fMaxSeedVertexZ     = 200.; 
-  const  Float_t   AliTRDtracker::fLabelFraction      = 0.5;  
-  const  Float_t   AliTRDtracker::fWideRoad           = 30.;
+  const  Float_t     AliTRDtracker::fMinClustersInSeed  = 0.7;  
 
-  const  Double_t  AliTRDtracker::fMaxChi2            = 12.; 
-  const  Double_t  AliTRDtracker::fSeedErrorSY        = 0.1;
-  const  Double_t  AliTRDtracker::fSeedErrorSY3       = 2.5;
-  const  Double_t  AliTRDtracker::fSeedErrorSZ        = 0.1;
+  const  Float_t     AliTRDtracker::fMinClustersInTrack = 0.5;  
+  const  Float_t     AliTRDtracker::fSkipDepth          = 0.05;
+  const  Float_t     AliTRDtracker::fLabelFraction      = 0.5;  
+  const  Float_t     AliTRDtracker::fWideRoad           = 20.;
+
+  const  Double_t    AliTRDtracker::fMaxChi2            = 24.; 
 
 //____________________________________________________________________
 AliTRDtracker::AliTRDtracker()
@@ -107,6 +104,7 @@ AliTRDtracker::AliTRDtracker()
   fNtracks   = 0;
   fTracks    = NULL;
 
+  fSY2corr = 0.025;
 }   
 
 //____________________________________________________________________
@@ -123,6 +121,7 @@ AliTRDtracker::AliTRDtracker(const Text_t* name, const Text_t* title)
   fNtracks   = 0;
   fTracks    = new TObjArray(10000);
 
+  fSY2corr = 0.025;
 }   
 
 //___________________________________________________________________
@@ -135,18 +134,42 @@ AliTRDtracker::~AliTRDtracker()
 }   
 
 //___________________________________________________________________
-void AliTRDtracker::Clusters2Tracks()
+void AliTRDtracker::Clusters2Tracks(TH1F *hs, TH1F *hd)
 {
   Int_t inner, outer;
   Int_t fTotalNofTimeBins = fGeom->GetTimeMax() * AliTRDgeometry::Nplan();
-  Int_t nSteps = (Int_t) (fTotalNofTimeBins * fSeedDepth)/fSeedStep;
+  Int_t nSteps = (Int_t) (fSeedDepth / fSeedStep);
+  Int_t gap = (Int_t) (fTotalNofTimeBins * fSeedGap);
+  Int_t step = (Int_t) (fTotalNofTimeBins * fSeedStep);
+
+
+  //  nSteps = 1;
+
+  if (!fClusters) return; 
+
+  AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect];
+  SetUpSectors(fTrSec);
+
+  // make a first turn looking for seed ends in the same (n,n) 
+  // and in the adjacent sectors (n,n+1)
 
   for(Int_t i=0; i<nSteps; i++) {
     printf("step %d out of %d \n", i+1, nSteps);
-    outer=fTotalNofTimeBins-1-i*fSeedStep; inner=outer-fSeedGap;
-    MakeSeeds(inner,outer);
-    FindTracks();
+    outer=fTotalNofTimeBins-1-i*step; inner=outer-gap;
+    MakeSeeds(inner,outer, fTrSec, 1, hs, hd);
+    FindTracks(fTrSec, hs, hd);
   } 
+
+  // make a second turn looking for seed ends in next-to-adjacent 
+  // sectors (n,n+2)
+
+  for(Int_t i=0; i<nSteps; i++) {
+    printf("step %d out of %d \n", i+1, nSteps);
+    outer=fTotalNofTimeBins-1-i*step; inner=outer-gap;
+    MakeSeeds(inner, outer, fTrSec, 2, hs, hd);
+    FindTracks(fTrSec,hs,hd);
+  } 
+
 }          
 
 //_____________________________________________________________________
@@ -154,7 +177,7 @@ Double_t AliTRDtracker::ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt)
 {
   // Parametrised "expected" error of the cluster reconstruction in Y 
 
-  Double_t s = 0.2;    
+  Double_t s = 0.08 * 0.08;    
   return s;
 }
 
@@ -163,17 +186,14 @@ Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t r, Double_t tgl)
 {
   // Parametrised "expected" error of the cluster reconstruction in Z 
 
-  // Take an example pad size for the time being, check back
-  // again with Sergei
-  Double_t s, pad = fGeom->GetRowPadSize(0,2,0);
-  s = pad * pad /12.;  
+  Double_t s = 6 * 6 /12.;  
   return s;
 }                  
 
 //_____________________________________________________________________
 Double_t f1trd(Double_t x1,Double_t y1,
-                      Double_t x2,Double_t y2,
-                      Double_t x3,Double_t y3)
+              Double_t x2,Double_t y2,
+              Double_t x3,Double_t y3)
 {
   // Initial approximation of the track curvature
   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
@@ -191,8 +211,8 @@ Double_t f1trd(Double_t x1,Double_t y1,
 
 //_____________________________________________________________________
 Double_t f2trd(Double_t x1,Double_t y1,
-                      Double_t x2,Double_t y2,
-                      Double_t x3,Double_t y3)
+              Double_t x2,Double_t y2,
+              Double_t x3,Double_t y3)
 {
   // Initial approximation of the track curvature times center of curvature
   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
@@ -210,8 +230,8 @@ Double_t f2trd(Double_t x1,Double_t y1,
 
 //_____________________________________________________________________
 Double_t f3trd(Double_t x1,Double_t y1,
-                      Double_t x2,Double_t y2,
-                      Double_t z1,Double_t z2)
+              Double_t x2,Double_t y2,
+              Double_t z1,Double_t z2)
 {
   // Initial approximation of the tangent of the track dip angle
   // Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
@@ -223,13 +243,20 @@ Double_t f3trd(Double_t x1,Double_t y1,
 //___________________________________________________________________
 
 Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec,
-                            Int_t s, Int_t rf)
+                            Int_t s, Int_t rf, Int_t matched_index, 
+                                     TH1F *hs, TH1F *hd)
 {
   // Starting from current position on track=t this function tries 
   // to extrapolate the track up to timeBin=rf and to confirm prolongation
   // if a close cluster is found. *sec is a pointer to allocated
   // array of sectors, in which the initial sector has index=s. 
 
+
+  //  TH1F *hsame = hs;     
+  //  TH1F *hdiff = hd;   
+
+  //  Bool_t good_match;
+
   const Int_t TIME_BINS_TO_SKIP=Int_t(fSkipDepth*sec->GetNtimeBins());
   Int_t try_again=TIME_BINS_TO_SKIP;
 
@@ -237,22 +264,29 @@ Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec,
 
   Int_t ns=Int_t(2*TMath::Pi()/alpha+0.5);
 
+  Double_t x0, rho;
+
   for (Int_t nr=sec->GetTimeBinNumber(t.GetX())-1; nr>=rf; nr--) {
-    Double_t x=sec->GetX(nr), ymax=sec->GetMaxY(nr);
-    if (!t.PropagateTo(x)) {
+
+    Double_t x=sec->GetX(nr);
+    Double_t ymax=x*TMath::Tan(0.5*alpha);
+
+    rho = 0.00295; x0 = 11.09;  // TEC
+    if(sec->TECframe(nr,t.GetY(),t.GetZ())) { 
+      rho = 1.7; x0 = 33.0;     // G10 frame 
+    } 
+    if(TMath::Abs(x - t.GetX()) > 3) { 
+      rho = 0.0559; x0 = 55.6;  // radiator
+    }
+    if (!t.PropagateTo(x,x0,rho,0.139)) {
       cerr<<"Can't propagate to x = "<<x<<endl;
       return 0;
     }
 
-    AliTRDcluster *cl=0;
-    UInt_t index=0;
-
-    Double_t max_chi2=fMaxChi2;
-
     AliTRDtimeBin& time_bin=sec[s][nr];
     Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
     Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
-    Double_t road=5.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
+    Double_t road=25.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
 
     if (road>fWideRoad) {
       if (t.GetNclusters() > 4) {
@@ -261,48 +295,94 @@ Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec,
       return 0;
     }       
 
+    AliTRDcluster *cl=0;
+    UInt_t index=0;
+    //    Int_t ncl = 0;
+
+    Double_t max_chi2=fMaxChi2;
+
     if (time_bin) {
+
       for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
         AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
+
+       //      good_match = kFALSE;
+       //      if((c->GetTrackIndex(0) == matched_index) ||
+       //   (c->GetTrackIndex(1) == matched_index) ||
+       //   (c->GetTrackIndex(2) == matched_index)) good_match = kTRUE;
+       //        if(hsame) hsame->Fill(TMath::Abs(c->GetY()-y)/road);
+       //        if(hdiff) hdiff->Fill(road);
+
         if (c->GetY() > y+road) break;
         if (c->IsUsed() > 0) continue;
 
-       if((c->GetZ()-z)*(c->GetZ()-z) > 25. + sz2) continue;
+       //      if(good_match) hsame->Fill(TMath::Abs(c->GetZ()-z));
+       //      else hdiff->Fill(TMath::Abs(c->GetZ()-z));
+
+       //      if(!good_match) continue;
+
+       if((c->GetZ()-z)*(c->GetZ()-z) > 3 * 12 * sz2) continue;
 
         Double_t chi2=t.GetPredictedChi2(c);
 
+       //      if((c->GetTrackIndex(0) == matched_index) ||
+       //         (c->GetTrackIndex(1) == matched_index) ||
+       //         (c->GetTrackIndex(2) == matched_index))
+       //        hdiff->Fill(chi2);
+
+       //      ncl++;
+
         if (chi2 > max_chi2) continue;
         max_chi2=chi2;
         cl=c;
         index=time_bin.GetIndex(i);
       }   
+      
+      if(!cl) {
+
+       for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
+         AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
+
+         if (c->GetY() > y+road) break;
+         if (c->IsUsed() > 0) continue;          
+         if((c->GetZ()-z)*(c->GetZ()-z) > 3.25 * 12 * sz2) continue;
+         
+         Double_t chi2=t.GetPredictedChi2(c);
+         
+         //      ncl++;
+
+         if (chi2 > max_chi2) continue;
+         max_chi2=chi2;
+         cl=c;
+         index=time_bin.GetIndex(i);
+       }   
+      }
+      
     }
+
     if (cl) {
 
-      //      Float_t l=sec->GetPitch();
-      //      t.SetSampledEdx(cl->fQ/l,Int_t(t));  
-     
       t.Update(cl,max_chi2,index);
 
       try_again=TIME_BINS_TO_SKIP;
     } else {
       if (try_again==0) break;
       if (y > ymax) {
-       cerr<<"y > ymax: "<<y<<" > "<<ymax<<endl;
+       //      cerr<<"y > ymax: "<<y<<" > "<<ymax<<endl;
          s = (s+1) % ns;
          if (!t.Rotate(alpha)) {
           cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
           return 0;
         }
       } else if (y <-ymax) {
-       cerr<<"y < -ymax: "<<y<<" < "<<-ymax<<endl;
+       //      cerr<<"y < -ymax: "<<y<<" < "<<-ymax<<endl;
          s = (s-1+ns) % ns;
          if (!t.Rotate(-alpha)) {
           cerr<<"Failed to rotate, alpha = "<<alpha<<endl;
           return 0;
         }
       }
-      try_again--;
+      if(!sec->TECframe(nr,y,z)) try_again--;
     }
   }
 
@@ -314,7 +394,7 @@ Int_t AliTRDtracker::FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec,
 //_____________________________________________________________________________
 void AliTRDtracker::GetEvent(const Char_t *hitfile, const Char_t *clusterfile)
 {
-  // Opens a ROOT-file with TRD-clusters and reads in the cluster-tree
+  // Opens a ROOT-file with TRD-clusters and reads the cluster-tree in
 
   ReadClusters(fClusters, clusterfile);
 
@@ -324,7 +404,7 @@ void AliTRDtracker::GetEvent(const Char_t *hitfile, const Char_t *clusterfile)
   if (!fInputFile) {
     printf("AliTRDtracker::Open -- ");
     printf("Open the ALIROOT-file %s.\n",hitfile);
-    fInputFile = new TFile(hitfile);
+    fInputFile = new TFile(hitfile,"UPDATE");
   }
   else {
     printf("AliTRDtracker::Open -- ");
@@ -343,6 +423,7 @@ void AliTRDtracker::GetEvent(const Char_t *hitfile, const Char_t *clusterfile)
     printf("Could not find AliRun object.\n");
   }
 
+  /*  
   // Import the Trees for the event nEvent in the file
   Int_t nparticles = gAlice->GetEvent(fEvent);
   cerr<<"nparticles = "<<nparticles<<endl;
@@ -351,6 +432,7 @@ void AliTRDtracker::GetEvent(const Char_t *hitfile, const Char_t *clusterfile)
     printf("AliTRDtracker::GetEvent -- ");
     printf("No entries in the trees for event %d.\n",fEvent);
   }
+  */  
 
   AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD");
   fGeom = TRD->GetGeometry();
@@ -369,38 +451,43 @@ void AliTRDtracker::SetUpSectors(AliTRDtrackingSector *sec)
 
   //  Sort clusters into AliTRDtimeBin's within AliTRDtrackSector's 
 
-  cerr<<"MakeSeeds: sorting clusters"<<endl;
+  cerr<<"SetUpSectors: sorting clusters"<<endl;
               
   Int_t ncl=fClusters->GetEntriesFast();
   UInt_t index;
   while (ncl--) {
-    printf("\r %d left",ncl); 
+    printf("\r %d left  ",ncl); 
     AliTRDcluster *c=(AliTRDcluster*)fClusters->UncheckedAt(ncl);
     Int_t detector=c->GetDetector(), local_time_bin=c->GetLocalTimeBin();
     Int_t sector=fGeom->GetSector(detector);
 
-    Int_t tracking_sector=sector;
-    if(sector > 0) tracking_sector = AliTRDgeometry::kNsect - sector;
+    Int_t tracking_sector = AliTRDgeometry::kNsect - sector - 1;
 
     Int_t tb=sec[sector].GetTimeBin(detector,local_time_bin); 
     index=ncl;
     sec[tracking_sector][tb].InsertCluster(c,index);
+
   }    
   printf("\r\n");
 }
 
 
 //_____________________________________________________________________________
-void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer)
+void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, 
+                              AliTRDtrackingSector* fTrSec, Int_t turn,
+                             TH1F *hs, TH1F *hd)
 {
   // Creates track seeds using clusters in timeBins=i1,i2
 
   Int_t i2 = inner, i1 = outer; 
+  Int_t ti[3], to[3];
+  Int_t nprim = 85600/2;
 
-  if (!fClusters) return; 
 
-  AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect];
-  SetUpSectors(fTrSec);
+  TH1F *hsame = hs;
+  TH1F *hdiff = hd;   
+  Bool_t match = false;
+  Int_t matched_index;
 
   // find seeds
 
@@ -410,51 +497,114 @@ void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer)
   Double_t alpha=AliTRDgeometry::GetAlpha();
   Double_t shift=AliTRDgeometry::GetAlpha()/2.;
   Double_t cs=cos(alpha), sn=sin(alpha);  
+  Double_t cs2=cos(2.*alpha), sn2=sin(2.*alpha);  
 
   Double_t x1 =fTrSec[0].GetX(i1);
   Double_t xx2=fTrSec[0].GetX(i2);  
 
+
+  printf("\n");
+
+  if((turn != 1)&&(turn != 2)) {
+    printf("*** Error in MakeSeeds: unexpected turn = %d  \n", turn);
+    return;
+  }
+
+
   for (Int_t ns=0; ns<max_sec; ns++) {
 
+    printf("\n MakeSeeds: sector %d \n", ns); 
+
+    Int_t nl2=fTrSec[(ns-2+max_sec)%max_sec][i2]; 
     Int_t nl=fTrSec[(ns-1+max_sec)%max_sec][i2];
     Int_t nm=fTrSec[ns][i2];
     Int_t nu=fTrSec[(ns+1)%max_sec][i2];
+    Int_t nu2=fTrSec[(ns+2)%max_sec][i2]; 
 
     AliTRDtimeBin& r1=fTrSec[ns][i1];
 
     for (Int_t is=0; is < r1; is++) {
       Double_t y1=r1[is]->GetY(), z1=r1[is]->GetZ();
+      for(Int_t ii=0; ii<3; ii++) to[ii] = r1[is]->GetTrackIndex(ii); 
 
-      for (Int_t js=0; js < nl+nm+nu; js++) {
+      for (Int_t js=0; js < nl2+nl+nm+nu+nu2; js++) {
+         
        const AliTRDcluster *cl;
-        Double_t x2,   y2,   z2;
-        Double_t x3=0.,y3=0.;  
+       Double_t x2,   y2,   z2;
+       Double_t x3=0., y3=0.;  
 
-        if (js<nl) {
+       if (js<nl2) {
+         if(turn != 2) continue;
+         AliTRDtimeBin& r2=fTrSec[(ns-2+max_sec)%max_sec][i2];
+         cl=r2[js];
+         y2=cl->GetY(); z2=cl->GetZ();
+         for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
+
+         x2= xx2*cs2+y2*sn2;
+         y2=-xx2*sn2+y2*cs2;
+       }        
+       else if (js<nl2+nl) {
+         if(turn != 1) continue;
          AliTRDtimeBin& r2=fTrSec[(ns-1+max_sec)%max_sec][i2];
-         cl=r2[js]; 
+         cl=r2[js-nl2];
          y2=cl->GetY(); z2=cl->GetZ();
+         for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
 
-          x2= xx2*cs+y2*sn;
-          y2=-xx2*sn+y2*cs;          
+         x2= xx2*cs+y2*sn;
+         y2=-xx2*sn+y2*cs;
 
-        } else
-          if (js<nl+nm) {
-           AliTRDtimeBin& r2=fTrSec[ns][i2];
-           cl=r2[js-nl];
-           x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
-          } else {
-           AliTRDtimeBin& r2=fTrSec[(ns+1)%max_sec][i2];
-            cl=r2[js-nl-nm];
-           y2=cl->GetY(); z2=cl->GetZ();
+       }
+       else if (js<nl2+nl+nm) {
+         if(turn != 1) continue;
+         AliTRDtimeBin& r2=fTrSec[ns][i2];
+         cl=r2[js-nl2-nl];
+         x2=xx2; y2=cl->GetY(); z2=cl->GetZ();
+         for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
+       }
+       else if (js<nl2+nl+nm+nu) {
+         if(turn != 1) continue;
+         AliTRDtimeBin& r2=fTrSec[(ns+1)%max_sec][i2];
+         cl=r2[js-nl2-nl-nm];
+         y2=cl->GetY(); z2=cl->GetZ();
+         for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
 
-            x2=xx2*cs-y2*sn;
-            y2=xx2*sn+y2*cs;   
+         x2=xx2*cs-y2*sn;
+         y2=xx2*sn+y2*cs;
 
-          }
+       }
+       else {
+         if(turn != 2) continue;
+         AliTRDtimeBin& r2=fTrSec[(ns+2)%max_sec][i2];
+         cl=r2[js-nl2-nl-nm-nu];
+         y2=cl->GetY(); z2=cl->GetZ();
+         for(Int_t ii=0; ii<3; ii++) ti[ii] = cl->GetTrackIndex(ii);
+         
+         x2=xx2*cs2-y2*sn2;
+         y2=xx2*sn2+y2*cs2;
+       }         
+       
 
-        Double_t zz=z1 - z1/x1*(x1-x2);
+       match = false;
+       matched_index = -1;
+       for (Int_t ii=0; ii<3; ii++) {
+         // cerr<<"ti["<<ii<<"] = "<<ti[ii]<<"; to["<<ii<<"] = "<<to[ii]<<endl;
+         if(ti[ii] < 0) continue;
+         if(ti[ii] >= nprim) continue;
+         for (Int_t kk=0; kk<3; kk++) {
+           if(to[kk] < 0) continue;
+           if(to[kk] >= nprim) continue;
+           if(ti[ii] == to[kk]) {
+             //cerr<<"ti["<<ii<<"] = "<<ti[ii]<<" = "<<to[kk]<<" = to["<<kk<<"]"<<endl;
+             matched_index = ti[ii];
+             match = true;
+           }
+         }
+       }                 
        
+       if(TMath::Abs(z1-z2) > fMaxSeedDeltaZ12) continue;
+
+        Double_t zz=z1 - z1/x1*(x1-x2);
+
         if (TMath::Abs(zz-z2)>fMaxSeedDeltaZ) continue;   
 
         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
@@ -464,7 +614,7 @@ void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer)
         x[1]=z1;
         x[2]=f1trd(x1,y1,x2,y2,x3,y3);
 
-        if (TMath::Abs(x[2]) >= fMaxSeedC) continue;
+        if (TMath::Abs(x[2]) > fMaxSeedC) continue;
 
         x[3]=f2trd(x1,y1,x2,y2,x3,y3);
 
@@ -476,6 +626,7 @@ void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer)
 
         Double_t a=asin(x[3]);
         Double_t zv=z1 - x[4]/x[2]*(a+asin(x[2]*x1-x[3]));
+
         if (TMath::Abs(zv)>fMaxSeedVertexZ) continue;    
 
         Double_t sy1=r1[is]->GetSigmaY2(), sz1=r1[is]->GetSigmaZ2();
@@ -503,15 +654,21 @@ void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer)
         c[14]=f40*sy1*f40+f41*sz1*f41+f42*sy2*f42+f43*sz2*f43;   
 
         UInt_t index=r1.GetIndex(is);
-        AliTRDtrack *track=new AliTRDtrack(index, x, c, x1, ns*alpha+shift); 
+       
+        AliTRDtrack *track=new AliTRDtrack(r1[is],index,x,c,x1,ns*alpha+shift); 
+
+        Int_t rc=FindProlongation(*track,fTrSec,ns,i2,matched_index,hsame,hdiff);
+
+       //      if (match) hsame->Fill((Float_t) track->GetNclusters());
+       //      else hdiff->Fill((Float_t) track->GetNclusters());  
+       //      delete track;
+       //      continue;
 
-        Int_t rc=FindProlongation(*track,fTrSec,ns,i2);
-       
         if ((rc < 0) || 
             (track->GetNclusters() < (i1-i2)*fMinClustersInSeed)) delete track;
         else { 
          fSeeds->AddLast(track); fNseeds++; 
-         cerr<<"found seed "<<fNseeds<<endl;
+         printf("\r found seed %d  ", fNseeds);
        }
       }
     }
@@ -520,15 +677,67 @@ void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer)
   fSeeds->Sort();
 }          
 
-//___________________________________________________________________
-void AliTRDtracker::FindTracks(
+//_____________________________________________________________________________
+void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename
 {
-  if (!fClusters) return; 
+  //
+  // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0) 
+  // from the file. The names of the cluster tree and branches 
+  // should match the ones used in AliTRDclusterizer::WriteClusters()
+  //
 
-  AliTRDtrackingSector fTrSec[AliTRDgeometry::kNsect];
-  SetUpSectors(fTrSec);
+  TDirectory *savedir=gDirectory; 
+
+  TFile *file = TFile::Open(filename);
+  if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;} 
+
+  TTree *ClusterTree = (TTree*)file->Get("ClusterTree");
 
-  // find tracks
+  TObjArray *ClusterArray = new TObjArray(400); 
+  ClusterTree->GetBranch("TRDcluster")->SetAddress(&ClusterArray); 
+  
+  Int_t nEntries = (Int_t) ClusterTree->GetEntries();
+  printf("found %d entries in %s.\n",nEntries,ClusterTree->GetName());
+
+  // Loop through all entries in the tree
+  Int_t nbytes;
+  AliTRDcluster *c = 0;
+
+  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();  
+    printf("Read %d clusters from entry %d \n", nCluster, iEntry);
+
+    // Loop through all TRD digits
+    for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
+      c = (AliTRDcluster*)ClusterArray->UncheckedAt(iCluster);
+      AliTRDcluster *co = new AliTRDcluster(*c);
+      co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
+      array->AddLast(co);
+      delete ClusterArray->RemoveAt(iCluster); 
+    }
+  }
+
+  file->Close();                   
+  delete ClusterArray;
+  savedir->cd(); 
+  
+}
+
+//___________________________________________________________________
+void AliTRDtracker::FindTracks(AliTRDtrackingSector* fTrSec, TH1F *hs, TH1F *hd) 
+{
+  //
+  // Finds tracks in TRD
+  //
+
+  TH1F *hsame = hs;
+  TH1F *hdiff = hd;   
 
   Int_t num_of_time_bins = fTrSec[0].GetNtimeBins(); 
   Int_t nseed=fSeeds->GetEntriesFast();
@@ -546,11 +755,14 @@ void AliTRDtracker::FindTracks()
     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
     Int_t ns=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;
 
-    if (FindProlongation(t,fTrSec,ns)) {
+    Int_t label = GetTrackLabel(t);
+
+    if (FindProlongation(t,fTrSec,ns,0,label,hsame,hdiff)) {
       cerr<<"No of clusters in the track = "<<t.GetNclusters()<<endl; 
       if (t.GetNclusters() >= Int_t(fMinClustersInTrack*num_of_time_bins)) {
        Int_t label = GetTrackLabel(t);
        t.SetLabel(label);
+       t.CookdEdx();
        UseClusters(t);
 
         AliTRDtrack *pt = new AliTRDtrack(t);
@@ -560,6 +772,7 @@ void AliTRDtracker::FindTracks()
       }                         
     }     
     delete fSeeds->RemoveAt(i);  
+    fNseeds--;
   }            
 }
 
@@ -638,7 +851,17 @@ Int_t AliTRDtracker::WriteTracks(const Char_t *filename) {
 
   TDirectory *savedir=gDirectory;   
 
-  TFile *out=TFile::Open(filename,"RECREATE");
+  //TFile *out=TFile::Open(filename,"RECREATE");
+  TFile *out = (TFile*) gROOT->GetListOfFiles()->FindObject(filename);
+  if (!out) {
+    printf("AliTRDtracker::Open -- ");
+    printf("Open the ALIROOT-file %s.\n",filename);
+    out = new TFile(filename,"RECREATE");
+  }
+  else {
+    printf("AliTRDtracker::Open -- ");
+    printf("%s is already open.\n",filename);
+  }
 
   TTree tracktree("TreeT","Tree with TRD tracks");
 
@@ -660,68 +883,6 @@ Int_t AliTRDtracker::WriteTracks(const Char_t *filename) {
   savedir->cd();  
 
   cerr<<"WriteTracks: done"<<endl;
-  return 0;
-}
-
-//_____________________________________________________________________________
-void AliTRDtracker::ReadClusters(TObjArray *array, const Char_t *filename, 
-Int_t option) 
-{
-  //
-  // Reads AliTRDclusters (option >= 0) or AliTRDrecPoints (option < 0) 
-  // from the file. The names of the cluster tree and branches 
-  // should match the ones used in AliTRDclusterizer::WriteClusters()
-  //
-
-  TDirectory *savedir=gDirectory; 
-
-  TFile *file = TFile::Open(filename);
-  if (!file->IsOpen()) {printf("Can't open file %s !\n",filename); return;} 
-
-  TTree *tree = (TTree*)file->Get("ClusterTree");
-  Int_t nentr = (Int_t) tree->GetEntries();
-  printf("found %d entries in %s.\n",nentr,tree->GetName());
-
-  TBranch *branch;
-  TObjArray *ioArray = new TObjArray(400);
-
-  if( option < 0 ) {
-    //branch = tree->GetBranch("RecPoints");
-    // changed CBL
-    branch = tree->GetBranch("TRDrecPoints");
-
-    for (Int_t i=0; i<nentr; i++) {
-      branch->SetAddress(&ioArray);
-      tree->GetEvent(i);
-      Int_t npoints = ioArray->GetEntriesFast();
-      printf("Read %d rec. points from entry %d \n", npoints, i);
-
-      for(Int_t j=0; j<npoints; j++) {
-       AliTRDrecPoint *p=(AliTRDrecPoint*)ioArray->UncheckedAt(j);
-       array->AddLast(p);
-       ioArray->RemoveAt(j);
-      }
-    }
-  }
-  else {
-    branch = tree->GetBranch("TRDcluster");
-
-    for (Int_t i=0; i<nentr; i++) {
-      branch->SetAddress(&ioArray);
-      tree->GetEvent(i);
-      Int_t npoints = ioArray->GetEntriesFast();
-      printf("Read %d clusters from entry %d \n", npoints, i);
-
-      for(Int_t j=0; j<npoints; j++) {
-       AliTRDcluster *c=(AliTRDcluster*)ioArray->UncheckedAt(j);
-       array->AddLast(c);
-       ioArray->RemoveAt(j);
-      }
-    }
-  }
-
-  file->Close();                   
-  savedir->cd(); 
-  
+  return 1;
 }
 
index e0b7d88..bc67eb3 100644 (file)
@@ -5,6 +5,7 @@
  * See cxx source for full Copyright notice                               */ 
 
 #include <TNamed.h>
+#include <TH1.h>   
 
 class TFile;
 class TParticle;
@@ -12,10 +13,9 @@ class TParticlePDG;
 class TObjArray;
 
 class AliTRDgeometry;
-class AliTRDtrackingSector;
 class AliTRDtrack;
 class AliTRDmcTrack;
-
+class AliTRDtrackingSector;
 
 class AliTRDtracker : public TNamed { 
 
@@ -25,20 +25,41 @@ class AliTRDtracker : public TNamed {
   AliTRDtracker(const Text_t* name, const Text_t* title);
   ~AliTRDtracker(); 
 
-  virtual void  Clusters2Tracks(); 
+  virtual void  Clusters2Tracks(TH1F *hs, TH1F *hd); 
   Double_t      ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt);
   Double_t      ExpectedSigmaZ2(Double_t r, Double_t tgl);
   Int_t         FindProlongation(AliTRDtrack& t, AliTRDtrackingSector *sec,
-                                 Int_t s, Int_t rf=0);
+                              Int_t s, Int_t rf=0, Int_t matched_index = -1,
+                                TH1F *hs=0, TH1F *hd=0);
   void          GetEvent(const Char_t *hitfile, const Char_t *clusterfile);
   void          SetUpSectors(AliTRDtrackingSector *sec);
-  virtual void  MakeSeeds(Int_t inner, Int_t outer);
-  virtual void  FindTracks();
+  virtual void  MakeSeeds(Int_t inner, Int_t outer, AliTRDtrackingSector *sec,
+                         Int_t turn, TH1F *hs, TH1F *hd);
+  virtual void  FindTracks(AliTRDtrackingSector *sec, TH1F *hs, TH1F *hd);
   virtual void  UseClusters(AliTRDtrack t);
   virtual Int_t GetTrackLabel(AliTRDtrack t);
   Int_t         WriteTracks(const Char_t *filename); 
-  void          ReadClusters(TObjArray *array, const Char_t *filename, 
-                             Int_t option = 1);
+  void          ReadClusters(TObjArray *array, const Char_t *filename);
+
+  Float_t  GetSeedGap()       const {return fSeedGap;}   
+  Float_t  GetSeedStep()      const {return fSeedStep;}
+  Float_t  GetSeedDepth()     const {return fSeedDepth;}
+  Float_t  GetSkipDepth()     const {return fSkipDepth;}
+  Double_t GetMaxChi2()       const {return fMaxChi2;}
+  Float_t  GetMaxSeedC()      const {return fMaxSeedC;}
+  Float_t  GetMaxSeedTan()    const {return fMaxSeedTan;}
+  Double_t GetSeedErrorSY()   const {return fSeedErrorSY;}
+  Double_t GetSeedErrorSY3()  const {return fSeedErrorSY3;}
+  Double_t GetSeedErrorSZ()   const {return fSeedErrorSZ;}
+  Float_t  GetLabelFraction() const {return fLabelFraction;}
+  Float_t  GetWideRoad()      const {return fWideRoad;}
+
+  Float_t  GetMinClustersInTrack() const {return fMinClustersInTrack;}
+  Float_t  GetMinClustersInSeed()  const {return fMinClustersInSeed;} 
+  Float_t  GetMaxSeedDeltaZ()      const {return fMaxSeedDeltaZ;}
+  Float_t  GetMaxSeedVertexZ()     const {return fMaxSeedVertexZ;}
+
+  void     SetSY2corr(Float_t w)    {fSY2corr = w;}
 
  protected:
 
@@ -55,10 +76,14 @@ class AliTRDtracker : public TNamed {
   Int_t            fNtracks;          // Number of reconstructed tracks 
   TObjArray        *fTracks;          // List of reconstructed tracks   
 
-  static const Int_t    fSeedGap;  // Distance between inner and outer
-                                      // time bin in seeding
+  Float_t          fSY2corr;          // Correction coefficient for
+                                      // cluster SigmaY2 
+
+  static const Float_t  fSeedGap;   // Distance between inner and outer
+                                    // time bin in seeding 
+                                   // (fraction of all time bins) 
   
-  static const Int_t    fSeedStep;    // Step in iterations
+  static const Float_t  fSeedStep;    // Step in iterations
   static const Float_t         fSeedDepth;   // Fraction of TRD allocated for seeding
   static const Float_t  fSkipDepth;   // Fraction of TRD which can be skipped
                                       // in track prolongation            
@@ -66,7 +91,8 @@ class AliTRDtracker : public TNamed {
        
   static const Float_t  fMinClustersInTrack; // min fraction of clusters in track
   static const Float_t  fMinClustersInSeed;  // min fraction of clusters in seed
-  static const Float_t  fMaxSeedDeltaZ;  // max dZ in MakeSeeds
+  static const Float_t  fMaxSeedDeltaZ;   // max dZ in MakeSeeds
+  static const Float_t  fMaxSeedDeltaZ12; // max abs(z1-z2) in MakeSeeds
   static const Float_t  fMaxSeedC;       // max initial curvature in MakeSeeds
   static const Float_t  fMaxSeedTan;     // max initial Tangens(lambda) in MakeSeeds
   static const Float_t  fMaxSeedVertexZ; // max vertex Z in MakeSeeds
index 157d00a..68a5171 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.5  2000/12/08 16:07:02  cblume
+Update of the tracking by Sergei
+
 Revision 1.4  2000/10/16 01:16:53  cblume
 Changed timebin 0 to be the one closest to the readout
 
@@ -85,21 +88,20 @@ void AliTRDtrackingSector::SetUp()
 
 //______________________________________________________
 
-Double_t AliTRDtrackingSector::GetX(Int_t l) const
+Double_t AliTRDtrackingSector::GetX(Int_t tb) const
 {
-  if( (l<0) || (l>fN-1)) {
+  if( (tb<0) || (tb>fN-1)) {
     fprintf(stderr,"AliTRDtrackingSector::GetX: TimeBin index is out of range !\n");
     return -99999.;
   }
   else { 
     
     Int_t tb_per_plane = fN/AliTRDgeometry::Nplan();
-    Int_t plane = l/tb_per_plane;
-    Int_t time_slice = l%(Int_t(AliTRDgeometry::DrThick()
-                               /fTimeBinSize) + 1);
+    Int_t local_tb = tb_per_plane - tb%tb_per_plane - 1;
 
+    Int_t plane = tb/tb_per_plane;
     Float_t t0 = fGeom->GetTime0(plane);
-    Double_t x = t0 - time_slice * fTimeBinSize;
+    Double_t x = t0 - (local_tb + 0.5) * fTimeBinSize;
 
     return x;
   }
@@ -107,43 +109,31 @@ Double_t AliTRDtrackingSector::GetX(Int_t l) const
 
 //______________________________________________________
 
-Double_t AliTRDtrackingSector::GetMaxY(Int_t l) const 
-{ 
-
-  if((l<(fN-1)) && (l>-1)) { 
-    Int_t tb_per_plane = fN/AliTRDgeometry::Nplan();
-    Int_t plane = l/tb_per_plane;
-    return fGeom->GetChamberWidth(plane); 
-  }
-  else {
-    fprintf(stderr,
-    "AliTRDtrackingSector::GetMaxY: TimeBin index is out of range !\n");
-    if(l<0) return fGeom->GetChamberWidth(0);
-    else return fGeom->GetChamberWidth(AliTRDgeometry::Nplan()-1);
-  }
-}
-
-//______________________________________________________
-
 Int_t AliTRDtrackingSector::GetTimeBinNumber(Double_t x) const
 {
   Float_t r_out = fGeom->GetTime0(AliTRDgeometry::Nplan()-1); 
   Float_t r_in = fGeom->GetTime0(0) - AliTRDgeometry::DrThick();
 
+
   if(x >= r_out) return fN-1;
   if(x <= r_in) return 0;
 
-  Float_t gap = fGeom->GetTime0(1) - fGeom->GetTime0(0);
-
-  Int_t plane = Int_t((x - r_in + fTimeBinSize/2)/gap);
-
-  Int_t local_tb = Int_t((fGeom->GetTime0(plane)-x)/fTimeBinSize + 0.5);
-
-
+  Int_t plane;
+  for (plane = AliTRDgeometry::Nplan()-1; plane >= 0; plane--) {
+    if(x > (fGeom->GetTime0(plane) - AliTRDgeometry::DrThick())) break;
+  }  
   Int_t tb_per_plane = fN/AliTRDgeometry::Nplan();
+  Int_t local_tb = Int_t((fGeom->GetTime0(plane)-x)/fTimeBinSize);
 
-  Int_t time_bin = plane * (Int_t(AliTRDgeometry::DrThick()/fTimeBinSize) + 1)
-                   + (tb_per_plane - 1 - local_tb);
+  if((local_tb < 0) || (local_tb >= tb_per_plane)) {
+    printf("AliTRDtrackingSector::GetTimeBinNumber: \n");
+    printf("local time bin %d is out of bounds [0..%d]: x = %f \n",
+          local_tb,tb_per_plane-1,x);
+    return -1;
+  }
+      
+  Int_t time_bin = (plane + 1) * tb_per_plane - 1 - local_tb;
 
   return time_bin;
 }
@@ -156,9 +146,49 @@ Int_t AliTRDtrackingSector::GetTimeBin(Int_t det, Int_t local_tb) const
 
   Int_t tb_per_plane = fN/AliTRDgeometry::Nplan();
 
-  Int_t time_bin = plane * (Int_t(AliTRDgeometry::DrThick()/fTimeBinSize) + 1)
-                   + (tb_per_plane - 1 - local_tb);
+  Int_t time_bin = (plane + 1) * tb_per_plane - 1 - local_tb;
 
   return time_bin;
 }
 
+
+//______________________________________________________
+
+Bool_t AliTRDtrackingSector::TECframe(Int_t tb, Double_t y, Double_t z) const
+{
+// 
+// Returns <true> if point defined by <x(tb),y,z> is within 
+// the TEC G10 frame, otherwise returns <false>  
+//  
+
+  if((tb > (fN-1)) || (tb < 0)) return kFALSE; 
+
+  Int_t tb_per_plane = fN/AliTRDgeometry::Nplan();
+  Int_t plane = tb/tb_per_plane;
+  
+  Double_t x = GetX(tb);
+  y = TMath::Abs(y);
+
+  if((y > fGeom->GetChamberWidth(plane)/2.) &&
+     (y < x*TMath::Tan(0.5*AliTRDgeometry::GetAlpha()))) return kTRUE; 
+
+  Double_t zmin, zmax;
+  Float_t  fRowPadSize, fRow0;
+  Int_t    nPadRows;
+
+  for(Int_t iCha = 1; iCha < AliTRDgeometry::Ncham(); iCha++) {
+
+    fRow0 = fGeom->GetRow0(plane,iCha-1,0);
+    fRowPadSize = fGeom->GetRowPadSize(plane,iCha-1,0);
+    nPadRows = fGeom->GetRowMax(plane,iCha-1,0);
+    zmin = fRow0 - fRowPadSize/2 + fRowPadSize * nPadRows;
+
+    fRow0 = fGeom->GetRow0(plane,iCha,0);
+    fRowPadSize = fGeom->GetRowPadSize(plane,iCha,0);
+    zmax = fRow0 - fRowPadSize/2;
+
+    if((z > zmin) && (z < zmax)) return kTRUE;     
+  }
+
+  return kFALSE;
+}
index e79f549..7e38676 100644 (file)
@@ -23,13 +23,10 @@ public:
  
   AliTRDtimeBin& operator[](Int_t i);
   Int_t GetNtimeBins() const { return fN; }
-  Double_t GetX(Int_t l) const;
-  Double_t GetMaxY(Int_t l) const; 
-  //Double_t GetAlpha() const { return 2*TMath::Pi()/kNsect; } 
-  Int_t GetTimeBinNumber(Double_t x) const;
-  Int_t GetTimeBin(Int_t det, Int_t local_tb) const;
-  Float_t GetPitch() const {return fTimeBinSize;}   
-
+  Double_t GetX(Int_t tb) const;
+  Int_t   GetTimeBinNumber(Double_t x) const;
+  Int_t   GetTimeBin(Int_t det, Int_t local_tb) const;
+  Bool_t  TECframe(Int_t tb, Double_t y, Double_t z) const;
 
 protected:
 
@@ -37,10 +34,10 @@ protected:
   AliTRDgeometry          *fGeom;       // Pointer to TRD geometry
   AliTRDtimeBin           *fTimeBin;    // Pointer to array of AliTRDtimeBin
   Float_t                  fTimeBinSize;  // Time bin size in cm  
-
-
+                                             
   ClassDef(AliTRDtrackingSector,1)  // Provides tools to address clusters which lay within one sector
 
 }; 
 
+
 #endif