]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HLT/comp/AliL3Modeller.cxx
New version of TOF tracker which uses TOF clusters as an input (A. De Caro)
[u/mrichter/AliRoot.git] / HLT / comp / AliL3Modeller.cxx
index 1d068273cc92d478ef23f9d47fa0549ec7f6d52c..fd6f0a0a48520ef52a506a1a58bbae3def712a6f 100644 (file)
@@ -1,59 +1,71 @@
-//$Id$
+// @(#) $Id$
 
 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
-//*-- Copyright &copy ASV
+//*-- Copyright &copy ALICE HLT Group
+//_____________________________________________________________
+// AliL3Modeller
+//
+// Class for modeling TPC data.
+// 
+// This performs the cluster finding, based on track parameters.
+// Basically it propagates the tracks to all padrows, and looks 
+// for a corresponding cluster. For the moment only cog is calculated,
+// and no deconvolution is done. 
+
 
 #include "AliL3StandardIncludes.h"
 
+#include "AliL3Logging.h"
 #include "AliL3Modeller.h"
 #include "AliL3MemHandler.h"
-#ifdef use_aliroot
-#include "AliL3FileHandler.h"
-#endif
 #include "AliL3TrackArray.h"
 #include "AliL3ModelTrack.h"
 #include "AliL3DigitData.h"
 #include "AliL3Transform.h"
+#include "AliL3SpacePointData.h"
 
-#if GCCVERSION == 3
-using namespace std;
+#ifdef use_aliroot
+#include "AliL3FileHandler.h"
 #endif
 
-//_____________________________________________________________
-// AliL3Modeller
-//
-// Class for modeling TPC data.
-// 
-// This performs the cluster finding, based on track parameters.
-// Basically it propagates the tracks to all padrows, and looks 
-// for a corresponding cluster. For the moment only cog is calculated,
-// and no deconvolution is done. 
+#if __GNUC__ >= 3
+using namespace std;
+#endif
 
 ClassImp(AliL3Modeller)
 
 AliL3Modeller::AliL3Modeller()
 {
+  // default constructor
   fMemHandler=0;
   fTracks=0;
+  fRow=0;
   fTrackThreshold=0;
-  fPadOverlap=0;
   SetOverlap();
   SetTrackThreshold();
+  SetSearchRange();
+  SetMaxClusterRange(0,0);
+  fDebug=kFALSE;
 }
 
 
 AliL3Modeller::~AliL3Modeller()
 {
+  // destructor
   if(fMemHandler)
     delete fMemHandler;
   if(fTracks)
     delete fTracks;
+  if(fRow)
+    delete [] fRow;
 }
 
 void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path,Bool_t houghtracks,Bool_t binary)
 {
+  // Initialization
   fSlice = slice;
   fPatch = patch;
+  fHoughTracks=houghtracks;
 
   sprintf(fPath,"%s",path);
   
@@ -61,8 +73,11 @@ void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path,
   
   Char_t fname[100];
   AliL3MemHandler *file = new AliL3MemHandler();
-  //sprintf(fname,"%s/tracks_tr_%d_0.raw",trackdata,fSlice); //output tracks from the tracker (no merging)
-  sprintf(fname,"%s/tracks_ho_%d.raw",trackdata,fSlice);
+  if(!houghtracks)
+    sprintf(fname,"%s/tracks_tr_%d_0.raw",trackdata,fSlice); //output tracks from the tracker (no merging)
+  else 
+    sprintf(fname,"%s/tracks_ho_%d.raw",trackdata,fSlice);
+  //sprintf(fname,"%s/tracks_ho_%d_%d.raw",trackdata,fSlice,fPatch);
   if(!file->SetBinaryInput(fname))
     {
       cerr<<"AliL3Modeller::Init : Error opening trackfile: "<<fname<<endl;
@@ -72,11 +87,6 @@ void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path,
   file->CloseBinaryInput();
   delete file;
   
-  if(houghtracks)
-    cout<<"AliL3Modeller is assuming local hough tracksegments!"<<endl;
-  else
-    cout<<"AliL3Modeller is assuming global tracks!"<<endl;
-
   if(!houghtracks)
     fTracks->QSort();
   
@@ -93,9 +103,11 @@ void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path,
       track->CalculateHelix();
     }    
   
-  CalculateCrossingPoints();
+  Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
+  Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
+  Int_t bounds = ntimes*npads;
+  fRow = new Digit[bounds];
   
-  CheckForOverlaps();
   
   UInt_t ndigits=0;
   AliL3DigitRowData *digits=0;
@@ -106,7 +118,7 @@ void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path,
     {
       sprintf(fname,"%s/digitfile.root",fPath);
       fMemHandler->SetAliInput(fname);
-      digits = fMemHandler->AliDigits2Memory(ndigits);
+      digits = fMemHandler->AliAltroDigits2Memory(ndigits);
     }
   else
     {
@@ -143,7 +155,9 @@ void AliL3Modeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path,
 
 void AliL3Modeller::FindClusters()
 {
-  
+  // Finds clusters
+  if(fDebug)
+    cout<<"AliL3Modeller::FindClusters : Processing slice "<<fSlice<<" patch "<<fPatch<<endl;
   if(!fTracks)
     {
       cerr<<"AliL3Modeller::Process : No tracks"<<endl;
@@ -158,20 +172,20 @@ void AliL3Modeller::FindClusters()
   AliL3DigitRowData *rowPt = fRowData;
   AliL3DigitData *digPt=0;
 
-  Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
-  Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
-  Int_t bounds = ntimes*npads;
-  Digit *row = new Digit[bounds];
-  
-  Int_t seq_charge;
-  Int_t pad,time,index;
+  Int_t pad,time;
   Short_t charge;
   Cluster cluster;
-
+  ClusterRegion region[200];
+  
   for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
     {
+      if(i != (Int_t)rowPt->fRow)
+       {
+         cerr<<"AliL3Modeller::FindClusters : Mismatching rownumbering "<<i<<" "<<rowPt->fRow<<endl;
+         return;
+       }
       fCurrentPadRow = i;
-      memset((void*)row,0,ntimes*npads*sizeof(Digit));
+      memset((void*)fRow,0,(AliL3Transform::GetNTimeBins()+1)*(AliL3Transform::GetNPads(i)+1)*sizeof(Digit));
       digPt = (AliL3DigitData*)rowPt->fDigitData;
       //cout<<"Loading row "<<i<<" with "<<(Int_t)rowPt->fNDigit<<" digits"<<endl;
       for(UInt_t j=0; j<rowPt->fNDigit; j++)
@@ -179,8 +193,8 @@ void AliL3Modeller::FindClusters()
          pad = digPt[j].fPad;
          time = digPt[j].fTime;
          charge = digPt[j].fCharge;
-         row[ntimes*pad+time].fCharge = charge;
-         row[ntimes*pad+time].fUsed = kFALSE;
+         fRow[(AliL3Transform::GetNTimeBins()+1)*pad + time].fCharge = charge;
+         fRow[(AliL3Transform::GetNTimeBins()+1)*pad + time].fUsed = kFALSE;
          //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
        }
       
@@ -189,83 +203,40 @@ void AliL3Modeller::FindClusters()
          AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(k);
          if(!track) continue;
          
-         if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0 || track->GetOverlap(i)>=0)
+         if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0 || track->GetNOverlaps(i)>0)//track->GetOverlap(i)>=0)
            {
-             track->SetCluster(i,0,0,0,0,0,0); //The track has left the patch.
+             //cout<<"Track "<<k<<" is empty on row "<<i<<" "<<track->GetPadHit(i)<<" "<<track->GetTimeHit(i)<<endl;
+             track->SetCluster(i,0,0,0,0,0,0); //The track has left the patch, or it is overlapping
              continue;
            }
          
-         Int_t hitpad = (Int_t)rint(track->GetPadHit(i));
-         Int_t hittime = (Int_t)rint(track->GetTimeHit(i));
-         //cout<<"Checking track on row "<<i<<" with pad "<<hitpad<<" time "<<hittime<<endl;
-         pad = hitpad;
-         time = hittime;
-         Int_t padsign=-1;
-         Int_t timesign=-1;
+         Int_t minpad,mintime,maxpad,maxtime;
+         minpad = mintime = 999;
+         maxpad = maxtime = 0;
          
          memset(&cluster,0,sizeof(Cluster));
+         LocateCluster(track,region,minpad,maxpad);//,mintime,maxtime);
+         if(maxpad - minpad + 1 > fMaxPads ||  // maxtime - mintime + 1 > fMaxTimebins ||
+            maxpad - minpad < 1)               //  || maxtime - mintime < 1)
+           {
+             //cout<<"Cluster not found on row "<<i<<" maxpad "<<maxpad<<" minpad "<<minpad<<" maxtime "<<maxtime<<" mintime "<<mintime
+             //  <<" padhit "<<track->GetPadHit(i)<<" timehit "<<track->GetTimeHit(i)<<endl;
+               
+             track->SetCluster(i,0,0,0,0,0,0);
+             continue;
+           }
          
          Int_t npads=0;
-         while(1)//Process this padrow
+         for(pad=minpad; pad<=maxpad; pad++)
            {
-             if(pad < 0 || pad >= AliL3Transform::GetNPads(i)) 
-               {
-                 //cout<<"Pad = "<<pad<<" on row "<<i<<endl;
-                 FillCluster(track,&cluster,i,npads);
-                 break;
-               }
-             seq_charge=0;
-             timesign=-1;
-             time = hittime;
-             
-             while(1) //Process sequence on this pad:
+             Int_t ntimes=0;
+             for(time=region[pad].fMintime; time<=region[pad].fMaxtime; time++)
                {
-                 if(time < 0) break;
-                 index = ntimes*pad + time;
-                 if(index < 0 || index >= bounds)
-                   {
-                     cerr<<"AliL3Modeller::FindClusters : Index out of range : "<<index
-                       <<" on row "<<i<<" pad "<<pad<<" time "<<time<<endl;
-                     break;
-                   }
-                 
-                 charge = row[index].fCharge;
-                 if(charge==0 && timesign==-1) //zero charge on this timebin, perform checks:
-                   {
-                     if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap) //No charge found on this pad, look further.
-                       {
-                         time--;
-                         continue;
-                       }
-                     else //Boundary reached, or we have found one end of the sequence,->start looking in the other time direction
-                       {
-                         time = hittime+1;
-                         timesign=1;
-                         continue;
-                       }
-                   }
-                 else if(charge==0 && timesign==1)//zero charge on this timebin, perform checks:
-                   {
-                     if(seq_charge==0 && abs(time-hittime) <= fTimeOverlap)//No charge found on this pad, look further
-                       {
-                         time++;
-                         continue;
-                       }
-                     else //Boundary reached, or we have found the other end of the sequence, stop looking on this pad.
-                       {
-                         //if(fCurrentPadRow==31)
-                         //   cerr<<"Breaking off at pad "<<pad<<" and time "<<time<<endl;
-                         break;
-                       }
-                   }
-                 
-                 if(row[ntimes*pad+time].fUsed==kTRUE) //Don't use digits several times. This leads to mult. rec.tracks.
-                   {
-                     time += timesign;
-                     continue;
-                   }
-                 
-                 seq_charge += charge;
+                 charge = fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge;
+                 if(!charge) continue;
+                 if(fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed == kTRUE)
+                   continue;
+                 ntimes++;
                  
                  //Update the cluster parameters with this timebin
                  cluster.fTime += time*charge;
@@ -273,66 +244,19 @@ void AliL3Modeller::FindClusters()
                  cluster.fCharge += charge;
                  cluster.fSigmaY2 += pad*pad*charge;
                  cluster.fSigmaZ2 += time*time*charge;
-                 
-                 row[ntimes*pad+time].fUsed = kTRUE;
-                 time += timesign;
-               }
-             
-             
-             if(seq_charge)//There was something on this pad, so keep looking on the neighbouring pad
-               {
-                 pad += padsign;
-                 npads++;
-               }
-             else //Nothing more on this pad, goto next pad
-               {
-                 if(padsign==-1) 
-                   {
-                     if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap && pad > 0)
-                       {
-                         pad--; //In this case, we haven't found anything yet, 
-                       }        //so we will try to expand our search within the natural boundaries.
-                     else
-                       {
-                         pad=hitpad+1; 
-                         padsign=1; 
-                       }
-                     continue;
-                   }
-                 
-                 else if(padsign==1)
-                   {
-                     if(cluster.fCharge==0 && abs(pad-hitpad) <= fPadOverlap && pad < AliL3Transform::GetNPads(i)-2)
-                       {
-                         pad++;     //In this case, we haven't found anything yet, 
-                         continue;  //so we will try to expand our search within the natural boundaries.
-                       }
-                     else //We are out of range, or cluster if finished.
-                       {
-                         //if(fCurrentPadRow==31)
-                         //cout<<"Out of range; charge "<<cluster.fCharge<<" paddiff "<<abs(pad-hitpad)<<endl;
-                         FillCluster(track,&cluster,i,npads);
-                         break;
-                       }
-                   }
-                 else //Nothing more in this cluster
-                   {
-                     //if(fCurrentPadRow==31)
-                     //cout<<"Filling final cluster"<<endl;
-                     FillCluster(track,&cluster,i,npads);
-                     break;
-                   } 
+                 fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed = kTRUE;
                }
+             if(ntimes)
+               npads++;
            }
-         //cout<<"done"<<endl;
+         FillCluster(track,&cluster,i,npads);
        }
-      FillZeros(rowPt,row);
+      FillZeros(rowPt);
       fMemHandler->UpdateRowPointer(rowPt);
     }
-  delete [] row;
-  cout<<"done processing"<<endl;
-  
+  //cout<<"done processing"<<endl;
   
+
   //Debug:
   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
     {
@@ -344,8 +268,165 @@ void AliL3Modeller::FindClusters()
   
 }
 
+
+void AliL3Modeller::LocateCluster(AliL3ModelTrack *track,ClusterRegion *region,Int_t &padmin,Int_t &padmax)
+{
+  //Set the cluster range
+  //This method searches for _all_ nonzeros timebins which are neigbours.
+  //This makes it rather impractical when dealing with high occupancy,
+  //because then you might have very large "cluster" areas from low
+  //pt electrons/noise. 
+  
+  Int_t row=fCurrentPadRow,charge,prtmin=0,prtmax=999;
+  Int_t hitpad = (Int_t)rint(track->GetPadHit(row));
+  Int_t hittime = (Int_t)rint(track->GetTimeHit(row));
+  Int_t tmin = hittime;
+  Int_t tmax = tmin;
+    
+  Int_t clustercharge=0;
+  Int_t pad=hitpad;
+  Bool_t pm = kTRUE;
+  Int_t npads=0,middlemax=tmax,middlemin=tmin;
+  while(1)
+    {
+      Bool_t padpr=kFALSE;
+      Int_t time = hittime;
+      Bool_t tm = kTRUE;
+      if(pad < 0)
+       {
+         padmin = 0;
+         pad = hitpad+1;
+         pm = kFALSE;
+         prtmin = middlemin;
+         prtmax = middlemax;
+         continue;
+       }
+      else if(pad >= AliL3Transform::GetNPads(row))
+       {
+         padmax = AliL3Transform::GetNPads(row)-1;
+         break;
+       }
+      
+      tmin = 999;
+      tmax = 0;
+      //if(row==0)
+      //cout<<"Starting to look in pad "<<pad<<" time "<<time<<endl;
+      while(1)
+       {
+         if(time < 0)
+           {
+             time = hittime+1;
+             tm = kFALSE;
+           }
+         else if(time >= AliL3Transform::GetNTimeBins())
+           {
+             //timemax = AliL3Transform::GetNTimeBins()-1;
+             break;
+           }
+         charge = fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge;
+         //if(row==0)
+         //cout<<"charge "<<charge<<" at pad "<<pad<<" time "<<time<<endl;
+         if(charge>0)
+           {
+             clustercharge+=charge;
+             padpr = kTRUE;
+             if(time < tmin)
+               tmin = time;
+             if(time > tmax)
+               tmax = time;
+             if(tm)
+               time--;
+             else
+               time++;
+           }
+         else
+           {
+             if(tm)
+               {
+                 //if(abs(time - hittime) < fTimeSearch && padpr == kFALSE)//Keep looking
+                 if(time > prtmin && npads!=0)
+                   time--;
+                 else
+                   {
+                     time = hittime+1;
+                     tm=kFALSE;
+                   }
+               }
+             //else if(abs(time-hittime) < fTimeSearch && padpr == kFALSE)//Keep looking
+             else if(time < prtmax && npads != 0)
+               time++;
+             else
+               break;
+           }
+       }
+      if(npads==0)
+       {
+         middlemax = tmax;
+         middlemin = tmin;
+       }
+      //      if(row==0)
+      //cout<<"tmax "<<tmax<<" tmin "<<tmin<<" prtmin "<<prtmin<<" ptrmax "<<prtmax<<endl;
+      
+      if(padpr && tmax >= prtmin && tmin <= prtmax)//Sequence is overlapping with the previous
+       {
+         //if(row==0)
+         //cout<<"Incrementing pad "<<endl;
+         npads++;
+         
+         region[pad].fMintime=tmin;
+         region[pad].fMaxtime=tmax;
+         
+         /*
+         if(tmin < timemin)
+           timemin=tmin;
+         if(tmax > timemax)
+           timemax=tmax;
+         */
+         if(pad < padmin)
+           padmin = pad;
+         if(pad > padmax)
+           padmax = pad;
+         if(pm)
+           pad--;
+         else
+           pad++;
+         
+         prtmin = tmin;
+         prtmax = tmax;
+       }
+      else
+       {
+         if(pm)
+           {
+             if(abs(pad-hitpad)<fPadSearch && clustercharge == 0)
+               pad--;
+             else
+               {
+                 //if(row==0)
+                 //cout<<"Setting new pad "<<hitpad+1<<endl;
+                 pad = hitpad+1;
+                 pm = kFALSE;
+                 prtmin = middlemin;
+                 prtmax = middlemax;
+                 continue;
+               }
+           }
+         else 
+           {
+             if(abs(pad-hitpad)<fPadSearch && clustercharge==0)
+               pad++;
+             else
+               break;
+           }
+       }
+    }
+  
+}
+
+
 void AliL3Modeller::FillCluster(AliL3ModelTrack *track,Cluster *cluster,Int_t row,Int_t npads)
 {
+  // Fill clusters
   if(cluster->fCharge==0)
     {
       track->SetCluster(row,0,0,0,0,0,0);
@@ -357,20 +438,34 @@ void AliL3Modeller::FillCluster(AliL3ModelTrack *track,Cluster *cluster,Int_t ro
   Float_t sigmaY2,sigmaZ2;
   CalcClusterWidth(cluster,sigmaY2,sigmaZ2);
   track->SetCluster(row,fpad,ftime,fcharge,sigmaY2,sigmaZ2,npads);
+#ifdef do_mc
+  Int_t trackID[3];
+  GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
+  track->SetClusterLabel(row,trackID);
+#endif
 }
 
-void AliL3Modeller::FillZeros(AliL3DigitRowData *rowPt,Digit *row)
+
+
+void AliL3Modeller::FillZeros(AliL3DigitRowData *rowPt,Bool_t reversesign)
 {
   //Fill zero where data has been used.
-  
-  Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
+
   AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
   for(UInt_t j=0; j<rowPt->fNDigit; j++)
     {
       Int_t pad = digPt[j].fPad;
       Int_t time = digPt[j].fTime;
-      if(row[ntimes*pad+time].fUsed==kTRUE)
-       digPt[j].fCharge = 0;
+      if(fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed==kTRUE)
+       {
+         if(reversesign)
+           {
+             if(digPt[j].fCharge < 1024)
+               digPt[j].fCharge += 1024;
+           }
+         else
+           digPt[j].fCharge = 0;
+       }
     }
 }
 
@@ -381,7 +476,7 @@ void AliL3Modeller::WriteRemaining()
   AliL3DigitRowData *rowPt;
   rowPt = (AliL3DigitRowData*)fRowData;
   Int_t digitcount=0;
-  Int_t ndigits[(AliL3Transform::GetNRows(fPatch))];
+  Int_t *ndigits=new Int_t[(AliL3Transform::GetNRows(fPatch))];
   for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
     {
       AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
@@ -437,17 +532,43 @@ void AliL3Modeller::WriteRemaining()
   Char_t fname[100];
   AliL3MemHandler *mem = new AliL3MemHandler();
   sprintf(fname,"%s/comp/remains_%d_%d.raw",fPath,fSlice,fPatch);
+  mem->Init(fSlice,fPatch);
   mem->SetBinaryOutput(fname);
   mem->Memory2CompBinary((UInt_t)AliL3Transform::GetNRows(fPatch),(AliL3DigitRowData*)data);
   mem->CloseBinaryOutput();
   delete mem;
   delete [] data;
+  delete [] ndigits;
 }
 
+void AliL3Modeller::RemoveBadTracks()
+{
+  //Remove tracsk which should not be included in the compression scheme.
+
+  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
+    {
+      AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
+      if(!track) continue;
+
+      if(track->GetPt() < 0.08)
+       {
+         fTracks->Remove(i);
+         continue;
+       }
+
+      if(!fHoughTracks)
+       if(track->GetNHits() < fTrackThreshold)
+         fTracks->Remove(i);
+    }
+  fTracks->Compress();
+  
+}
 
 void AliL3Modeller::CalculateCrossingPoints()
 {
-  cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
+  // calculates crossing points
+  if(fDebug)
+    cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
   if(!fTracks)
     {
       cerr<<"AliL3Modeller::CalculateCrossingPoints(): No tracks"<<endl;
@@ -465,12 +586,12 @@ void AliL3Modeller::CalculateCrossingPoints()
 
          if(!track->GetCrossingPoint(i,hit)) 
            {
-             cerr<<"AliL3Modeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
-               "First point "<<track->GetFirstPointX()<<
-               " nhits "<<track->GetNHits()<<endl;//" tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
-               //"Center "<<track->GetCenterX()<<" "<<track->GetCenterY()<<endl<<endl<<
-               //"--------"<<endl;
-             fTracks->Remove(j);
+             //cerr<<"AliL3Modeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
+             //        " pt "<<track->GetPt()<<
+             //        " tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
+             //fTracks->Remove(j);
+             track->SetPadHit(i,-1);
+             track->SetTimeHit(i,-1);
              continue;
            }
          //cout<<"X "<<hit[0]<<" Y "<<hit[1]<<" Z "<<hit[2]<<" tgl "<<track->GetTgl()<<endl;
@@ -485,10 +606,13 @@ void AliL3Modeller::CalculateCrossingPoints()
              track->SetTimeHit(i,-1);
              continue;
            }
-         
-         
+
          track->SetPadHit(i,hit[1]);
          track->SetTimeHit(i,hit[2]);
+         track->CalculateClusterWidths(i);
+
+         Double_t beta = track->GetCrossingAngle(i);
+         track->SetCrossingAngleLUT(i,beta);
          
          //if(hit[1]<0 || hit[2]>445)
          //if(hit[2]<0 || hit[2]>445)
@@ -497,48 +621,62 @@ void AliL3Modeller::CalculateCrossingPoints()
        }
     }
   fTracks->Compress();
-  cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
+  if(fDebug)
+    cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
 }
 
-void AliL3Modeller::CheckForOverlaps()
+void AliL3Modeller::CheckForOverlaps(Float_t dangle,Int_t *rowrange)
 {
   //Flag the tracks that overlap
   
-  cout<<"Checking for overlaps...";
+  if(fDebug)
+    cout<<"Checking for overlaps on "<<fTracks->GetNTracks()<<endl;
   Int_t counter=0;
-  for(Int_t i=0; i<fTracks->GetNTracks(); i++)
+  
+  for(Int_t k=AliL3Transform::GetFirstRow(fPatch); k<=AliL3Transform::GetLastRow(fPatch); k++)
     {
-      AliL3ModelTrack *track1 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
-      if(!track1) continue;
-      for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
+      if(rowrange)
        {
-         AliL3ModelTrack *track2 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
-         if(!track2) continue;
-         for(Int_t k=AliL3Transform::GetFirstRow(fPatch); k<=AliL3Transform::GetLastRow(fPatch); k++)
+         if(k < rowrange[0]) continue;
+         if(k > rowrange[1]) break;
+       }
+      for(Int_t i=0; i<fTracks->GetNTracks(); i++)
+       {
+         AliL3ModelTrack *track1 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
+         if(!track1) continue;
+         if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0) continue;
+         
+         for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
            {
-             if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0 ||
-                track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0)
-               continue;
-             
-             if(track1->GetOverlap(k)>=0 || track2->GetOverlap(k)>=0) continue;
+             AliL3ModelTrack *track2 = (AliL3ModelTrack*)fTracks->GetCheckedTrack(j);
+             if(!track2) continue;
+             if(track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0) continue;
              
              if(abs((Int_t)rint(track1->GetPadHit(k))-(Int_t)rint(track2->GetPadHit(k))) <= fPadOverlap &&
                 abs((Int_t)rint(track1->GetTimeHit(k))-(Int_t)rint(track2->GetTimeHit(k))) <= fTimeOverlap)
                {
-                 track2->SetOverlap(k,i);
-                 //track1->SetOverlap(k,j);
+                 if(dangle>0 && fabs(track1->GetCrossingAngleLUT(k) - track2->GetCrossingAngleLUT(k)) < dangle)
+                   fTracks->Remove(j);
+                 
+                 //cout<<"row "<<k<<" "<<i<<" "<<j<<" "<<track1->GetPadHit(k)<<" "<<track2->GetPadHit(k)<<" "<<fabs(track1->GetCrossingAngleLUT(k) - track2->GetCrossingAngleLUT(k))<<endl;
+
+                 else
+                   track1->SetOverlap(k,j);
                  counter++;
                }
            }
        }
     }
-  cout<<"found "<<counter<<" done"<<endl;
+  fTracks->Compress();
+  if(fDebug)
+    cout<<"and there are "<<fTracks->GetNTracks()<<" track left"<<endl;
+  //cout<<"found "<<counter<<" done"<<endl;
 }
 
 
 void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
 {
-  
+  // calculates cluster's width
   Float_t padw,timew;
   
   padw = AliL3Transform::GetPadPitchWidth(fPatch);
@@ -567,7 +705,8 @@ void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigma
   
 
   
-  /*Constants added by offline
+  /*
+    Constants added by offline
     if(s2 != 0)
     {
     sigmaZ2 = sigmaZ2*0.169;
@@ -576,3 +715,44 @@ void AliL3Modeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigma
     }
   */
 }
+
+#ifdef do_mc
+void AliL3Modeller::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
+{
+  // Gets track ID
+  AliL3DigitRowData *rowPt = (AliL3DigitRowData*)fRowData;
+  
+  trackID[0]=trackID[1]=trackID[2]=-2;
+  
+  for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
+    {
+      if(rowPt->fRow < (UInt_t)fCurrentPadRow)
+       {
+         AliL3MemHandler::UpdateRowPointer(rowPt);
+         continue;
+       }
+      AliL3DigitData *digPt = (AliL3DigitData*)rowPt->fDigitData;
+      for(UInt_t j=0; j<rowPt->fNDigit; j++)
+       {
+         Int_t cpad = digPt[j].fPad;
+         Int_t ctime = digPt[j].fTime;
+         if(cpad != pad) continue;
+         if(ctime != time) continue;
+         //if(cpad != pad && ctime != ctime) continue;
+         //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
+         trackID[0] = digPt[j].fTrackID[0];
+         trackID[1] = digPt[j].fTrackID[1];
+         trackID[2] = digPt[j].fTrackID[2];
+         break;
+         //cout<<"Reading trackID "<<trackID[0]<<endl;
+       }
+      break;
+    }
+#else
+  void AliL3Modeller::GetTrackID(Int_t /*pad*/,Int_t /*time*/,Int_t */*trackID*/)
+{
+  // Does nothing if do_mc undefined
+  return;
+#endif
+}
+