HLT/comp/AliL3ClusterFitter.cxx

   1 // @(#) $Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ALICE HLT Group
   5
   6 #include "AliL3StandardIncludes.h"
   7
   8 #include "AliL3Logging.h"
   9 #include "AliL3ClusterFitter.h"
  10 #include "AliL3FitUtilities.h"
  11 #include "AliL3DigitData.h"
  12 #include "AliL3ModelTrack.h"
  13 #include "AliL3TrackArray.h"
  14 #include "AliL3MemHandler.h"
  15 #include "AliL3HoughTrack.h"
  16 #include "AliL3SpacePointData.h"
  17 #include "AliL3Compress.h"
  18
  19 #if __GNUC__ == 3
  20 using namespace std;
  21 #endif
  22
  23 /** /class AliL3ClusterFitter
  24 //<pre>
  25 //_____________________________________________________________
  26 //
  27 //  AliL3ClusterFitter
  28 //
  29 </pre>
  30 */
  31
  32 ClassImp(AliL3ClusterFitter)
  33
  34 Int_t AliL3ClusterFitter::fBadFitError=0;
  35 Int_t AliL3ClusterFitter::fFitError=0;
  36 Int_t AliL3ClusterFitter::fResultError=0;
  37 Int_t AliL3ClusterFitter::fFitRangeError=0;
  38
  39 AliL3ClusterFitter::AliL3ClusterFitter()
  40 {
  41   plane=0;
  42   fNmaxOverlaps = 3;
  43   fChiSqMax[0]=fChiSqMax[1]=fChiSqMax[2]=12;
  44   fRowMin=-1;
  45   fRowMax=-1;
  46   fFitted=0;
  47   fFailed=0;
  48   fYInnerWidthFactor=1;
  49   fZInnerWidthFactor=1;
  50   fYOuterWidthFactor=1;
  51   fZOuterWidthFactor=1;
  52   fSeeds=0;
  53   fProcessTracks=0;
  54   fClusters=0;
  55   fNMaxClusters=0;
  56   fNClusters=0;
  57   fEvent=0;
  58 }
  59
  60 AliL3ClusterFitter::AliL3ClusterFitter(Char_t *path)
  61 {
  62   strcpy(fPath,path);
  63   plane=0;
  64   fNmaxOverlaps = 3;
  65   fChiSqMax[0]=fChiSqMax[1]=fChiSqMax[2]=12;
  66   fRowMin=-1;
  67   fRowMax=-1;
  68   fFitted=0;
  69   fFailed=0;
  70   fYInnerWidthFactor=1;
  71   fZInnerWidthFactor=1;
  72   fYOuterWidthFactor=1;
  73   fZOuterWidthFactor=1;
  74   fSeeds=0;
  75   fProcessTracks=0;
  76   fNMaxClusters=100000;
  77   fClusters=0;
  78   fNClusters=0;
  79   fEvent=0;
  80 }
  81
  82 AliL3ClusterFitter::~AliL3ClusterFitter()
  83 {
  84   if(fSeeds)
  85     delete fSeeds;
  86   if(fClusters)
  87     delete [] fClusters;
  88 }
  89
  90 void AliL3ClusterFitter::Init(Int_t slice,Int_t patch,Int_t *rowrange,AliL3TrackArray *tracks)
  91 {
  92   //Assuming tracklets found by the line transform
  93
  94   fSlice=slice;
  95   fPatch=patch;
  96
  97   if(rowrange[0] > AliL3Transform::GetLastRow(patch) || rowrange[1] < AliL3Transform::GetFirstRow(patch))
  98     cerr<<"AliL3ClusterFitter::Init : Wrong rows "<<rowrange[0]<<" "<<rowrange[1]<<endl;
  99   fRowMin=rowrange[0];
 100   fRowMax=rowrange[1];
 101
 102   if(fRowMin < 0)
 103     fRowMin = 0;
 104   if(fRowMax > AliL3Transform::GetLastRow(fPatch))
 105     fRowMax = AliL3Transform::GetLastRow(fPatch);
 106
 107   fFitted=fFailed=0;
 108
 109   Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
 110   Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
 111   Int_t bounds = ntimes*npads;
 112   if(fRow)
 113     delete [] fRow;
 114   fRow = new Digit[bounds];
 115   if(fTracks)
 116     delete fTracks;
 117
 118   fTracks = new AliL3TrackArray("AliL3ModelTrack");
 119
 120   for(Int_t i=0; i<tracks->GetNTracks(); i++)
 121     {
 122       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
 123       if(!track) continue;
 124       AliL3ModelTrack *mtrack = (AliL3ModelTrack*)fTracks->NextTrack();
 125       mtrack->Init(slice,patch);
 126       mtrack->SetTgl(track->GetTgl());
 127       mtrack->SetRowRange(rowrange[0],rowrange[1]);
 128       for(Int_t j=fRowMin; j<=fRowMax; j++)
 129         {
 130           Float_t hit[3];
 131           track->GetLineCrossingPoint(j,hit);
 132           hit[0] += AliL3Transform::Row2X(track->GetFirstRow());
 133           Float_t R = sqrt(hit[0]*hit[0] + hit[1]*hit[1]);
 134           hit[2] = R*track->GetTgl();
 135           Int_t se,ro;
 136           AliL3Transform::Slice2Sector(slice,j,se,ro);
 137           AliL3Transform::Local2Raw(hit,se,ro);
 138           if(hit[1]<0 || hit[1]>=AliL3Transform::GetNPads(j) || hit[2]<0 || hit[2]>=AliL3Transform::GetNTimeBins())
 139             {
 140               mtrack->SetPadHit(j,-1);
 141               mtrack->SetTimeHit(j,-1);
 142               continue;
 143             }
 144           mtrack->SetPadHit(j,hit[1]);
 145           mtrack->SetTimeHit(j,hit[2]);
 146           mtrack->SetCrossingAngleLUT(j,fabs(track->GetPsiLine() - AliL3Transform::Pi()/2));
 147           //if(mtrack->GetCrossingAngleLUT(j) > AliL3Transform::Deg2Rad(20))
 148           //  cout<<"Angle "<<mtrack->GetCrossingAngleLUT(j)<<" psiline "<<track->GetPsiLine()*180/3.1415<<endl;
 149           mtrack->CalculateClusterWidths(j);
 150         }
 151     }
 152   //  cout<<"Copied "<<fTracks->GetNTracks()<<" tracks "<<endl;
 153 }
 154
 155 void AliL3ClusterFitter::Init(Int_t slice,Int_t patch)
 156 {
 157   fSlice=slice;
 158   fPatch=patch;
 159
 160   fRowMin=AliL3Transform::GetFirstRow(patch);
 161   fRowMax=AliL3Transform::GetLastRow(patch);
 162
 163   fFitted=fFailed=0;
 164
 165   Int_t ntimes = AliL3Transform::GetNTimeBins()+1;
 166   Int_t npads = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(fPatch))+1;//Max num of pads.
 167   Int_t bounds = ntimes*npads;
 168   if(fRow)
 169     delete [] fRow;
 170   fRow = new Digit[bounds];
 171   if(fTracks)
 172     delete fTracks;
 173   fTracks = new AliL3TrackArray("AliL3ModelTrack");
 174
 175 }
 176
 177 void AliL3ClusterFitter::LoadLocalSegments()
 178 {
 179   Char_t filename[1024];
 180   sprintf(filename,"%s/hough/tracks_ho_%d_%d_%d.raw",fPath,fEvent,fSlice,fPatch);
 181   AliL3MemHandler mem;
 182   mem.SetBinaryInput(filename);
 183   mem.Binary2TrackArray(fTracks);
 184   mem.CloseBinaryInput();
 185
 186   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 187     {
 188       AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
 189       if(!track) continue;
 190
 191       track->CalculateHelix();
 192
 193       track->Init(fSlice,fPatch);
 194
 195       for(Int_t j=fRowMin; j<=fRowMax; j++)
 196         {
 197           //Calculate the crossing point between track and padrow
 198
 199           Float_t xyz_cross[3];
 200           if(!track->GetCrossingPoint(j,xyz_cross))
 201             continue;
 202
 203           Int_t sector,row;
 204           AliL3Transform::Slice2Sector(fSlice,j,sector,row);
 205           AliL3Transform::Local2Raw(xyz_cross,sector,row);
 206
 207           if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j) ||
 208              xyz_cross[2] < 0 || xyz_cross[2] >= AliL3Transform::GetNTimeBins()) //track goes out of range
 209             continue;
 210
 211           track->SetPadHit(j,xyz_cross[1]);
 212           track->SetTimeHit(j,xyz_cross[2]);
 213
 214           Float_t crossingangle = track->GetCrossingAngle(j);
 215           track->SetCrossingAngleLUT(j,crossingangle);
 216           track->CalculateClusterWidths(j);
 217           track->GetClusterModel(j)->fSlice = fSlice;
 218
 219         }
 220     }
 221 }
 222
 223 void AliL3ClusterFitter::LoadSeeds(Int_t *rowrange,Bool_t offline,Int_t eventnr,Float_t zvertex)
 224 {
 225   //Function assumes _global_ tracks written to a single file.
 226   cout<<"Loading the seeds"<<endl;
 227   Char_t fname[1024];
 228   fEvent = eventnr;
 229
 230   if(offline)
 231     sprintf(fname,"%s/offline/tracks_%d.raw",fPath,fEvent);
 232   else
 233     sprintf(fname,"%s/hough/tracks_%d.raw",fPath,fEvent);
 234
 235   cout<<"AliL3ClusterFitter::LoadSeeds : Loading input tracks from "<<fname<<endl;
 236
 237   AliL3MemHandler tfile;
 238   tfile.SetBinaryInput(fname);
 239
 240   if(fSeeds)
 241     delete fSeeds;
 242   fSeeds = new AliL3TrackArray("AliL3ModelTrack");
 243
 244   tfile.Binary2TrackArray(fSeeds);
 245   tfile.CloseBinaryInput();
 246
 247   //if(!offline)
 248   //fSeeds->QSort();
 249
 250   Int_t clustercount=0;
 251   for(Int_t i=0; i<fSeeds->GetNTracks(); i++)
 252     {
 253       AliL3ModelTrack *track = (AliL3ModelTrack*)fSeeds->GetCheckedTrack(i);
 254       if(!track) continue;
 255
 256       if(!offline)
 257         {
 258           if(i==0) cerr<<"AliL3ClusterFitter::LoadSeeds : Cutting on pt of 4 GeV!!"<<endl;
 259           if(track->GetPt() > 4.)
 260             {
 261               fSeeds->Remove(i);
 262               continue;
 263             }
 264         }
 265       clustercount += track->GetNHits();
 266       track->CalculateHelix();
 267
 268       Int_t nhits = track->GetNHits();
 269       UInt_t *hitids = track->GetHitNumbers();
 270
 271       Int_t origslice = (hitids[nhits-1]>>25)&0x7f;//Slice of innermost point
 272       /*
 273       if(i==0) cerr<<"Cluster fitter only in HALF TPC!!!"<<endl;
 274       if(origslice > 17)
 275         {
 276           fSeeds->Remove(i);
 277           continue;
 278         }
 279       */
 280       track->Init(origslice,-1);
 281       Int_t slice = origslice;
 282
 283       //for(Int_t j=rowrange[1]; j>=rowrange[0]; j--)
 284       for(Int_t j=rowrange[0]; j<=rowrange[1]; j++)
 285         {
 286
 287           //Calculate the crossing point between track and padrow
 288           Float_t angle = 0; //Perpendicular to padrow in local coordinates
 289           AliL3Transform::Local2GlobalAngle(&angle,slice);
 290           if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(j)))
 291             {
 292               //cerr<<"No crossing in slice "<<slice<<" padrow "<<j<<endl;
 293               continue;
 294               //track->Print();
 295               //exit(5);
 296             }
 297           Float_t xyz_cross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
 298           xyz_cross[2] += zvertex;
 299
 300           Int_t sector,row;
 301           AliL3Transform::Slice2Sector(slice,j,sector,row);
 302           AliL3Transform::Global2Raw(xyz_cross,sector,row);
 303           //cout<<"Examining slice "<<slice<<" row "<<j<<" pad "<<xyz_cross[1]<<" time "<<xyz_cross[2]<<endl;
 304           if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j)) //Track leaves the slice
 305             {
 306             newslice:
 307
 308               Int_t tslice=slice;
 309               Float_t lastcross=xyz_cross[1];
 310               if(xyz_cross[1] > 0)
 311                 {
 312                   if(slice == 17)
 313                     slice=0;
 314                   else if(slice == 35)
 315                     slice = 18;
 316                   else
 317                     slice += 1;
 318                 }
 319               else
 320                 {
 321                   if(slice == 0)
 322                     slice = 17;
 323                   else if(slice==18)
 324                     slice = 35;
 325                   else
 326                     slice -= 1;
 327                 }
 328               if(slice < 0 || slice>35)
 329                 {
 330                   cerr<<"Wrong slice "<<slice<<" on row "<<j<<endl;
 331                   exit(5);
 332                 }
 333               //cout<<"Track leaving, trying slice "<<slice<<endl;
 334               angle=0;
 335               AliL3Transform::Local2GlobalAngle(&angle,slice);
 336               if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(j)))
 337                 {
 338                   cerr<<"No crossing in slice "<<slice<<" padrow "<<j<<endl;
 339                   continue;
 340                   //track->Print();
 341                   //exit(5);
 342                 }
 343               xyz_cross[0] = track->GetPointX();
 344               xyz_cross[1] = track->GetPointY();
 345               xyz_cross[2] = track->GetPointZ();
 346               xyz_cross[2] += zvertex;
 347               Int_t sector,row;
 348               AliL3Transform::Slice2Sector(slice,j,sector,row);
 349               AliL3Transform::Global2Raw(xyz_cross,sector,row);
 350               if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j)) //track is in the borderline
 351                 {
 352                   if(xyz_cross[1] > 0 && lastcross > 0 || xyz_cross[1] < 0 && lastcross < 0)
 353                     goto newslice;
 354                   else
 355                     {
 356                       slice = tslice;//Track is on the border of two slices
 357                       continue;
 358                     }
 359                 }
 360             }
 361
 362           if(xyz_cross[2] < 0 || xyz_cross[2] >= AliL3Transform::GetNTimeBins())//track goes out of range
 363             continue;
 364
 365           if(xyz_cross[1] < 0 || xyz_cross[1] >= AliL3Transform::GetNPads(j))
 366             {
 367               cerr<<"Slice "<<slice<<" padrow "<<j<<" pad "<<xyz_cross[1]<<" time "<<xyz_cross[2]<<endl;
 368               track->Print();
 369               exit(5);
 370             }
 371
 372           track->SetPadHit(j,xyz_cross[1]);
 373           track->SetTimeHit(j,xyz_cross[2]);
 374           angle=0;
 375           AliL3Transform::Local2GlobalAngle(&angle,slice);
 376           Float_t crossingangle = track->GetCrossingAngle(j,slice);
 377           track->SetCrossingAngleLUT(j,crossingangle);
 378
 379           track->CalculateClusterWidths(j);
 380
 381           track->GetClusterModel(j)->fSlice = slice;
 382
 383         }
 384       memset(track->GetHitNumbers(),0,159*sizeof(UInt_t));//Reset the hitnumbers
 385       track->SetNHits(0);
 386     }
 387   fSeeds->Compress();
 388
 389   cout<<"Loaded "<<fSeeds->GetNTracks()<<" seeds and "<<clustercount<<" clusters"<<endl;
 390 }
 391
 392 void AliL3ClusterFitter::FindClusters()
 393 {
 394   if(!fTracks)
 395     {
 396       cerr<<"AliL3ClusterFitter::Process : No tracks"<<endl;
 397       return;
 398     }
 399   if(!fRowData)
 400     {
 401       cerr<<"AliL3ClusterFitter::Process : No data "<<endl;
 402       return;
 403     }
 404
 405   AliL3DigitRowData *rowPt = fRowData;
 406   AliL3DigitData *digPt=0;
 407
 408   Int_t pad,time;
 409   Short_t charge;
 410
 411   if(fRowMin < 0)
 412     {
 413       fRowMin = AliL3Transform::GetFirstRow(fPatch);
 414       fRowMax = AliL3Transform::GetLastRow(fPatch);
 415     }
 416   for(Int_t i=AliL3Transform::GetFirstRow(fPatch); i<=AliL3Transform::GetLastRow(fPatch); i++)
 417     {
 418       if((Int_t)rowPt->fRow < fRowMin)
 419         {
 420           AliL3MemHandler::UpdateRowPointer(rowPt);
 421           continue;
 422         }
 423       else if((Int_t)rowPt->fRow > fRowMax)
 424         break;
 425       else if((Int_t)rowPt->fRow != i)
 426         {
 427           cerr<<"AliL3ClusterFitter::FindClusters : Mismatching row numbering "<<i<<" "<<rowPt->fRow<<endl;
 428           exit(5);
 429         }
 430       fCurrentPadRow = i;
 431       memset((void*)fRow,0,(AliL3Transform::GetNTimeBins()+1)*(AliL3Transform::GetNPads(i)+1)*sizeof(Digit));
 432       digPt = (AliL3DigitData*)rowPt->fDigitData;
 433
 434       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 435         {
 436           pad = digPt[j].fPad;
 437           time = digPt[j].fTime;
 438           charge = digPt[j].fCharge;
 439           fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fCharge = charge;
 440           fRow[(AliL3Transform::GetNTimeBins()+1)*pad+time].fUsed = kFALSE;
 441           //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
 442         }
 443
 444       for(Int_t it=0; it<2; it++)
 445         {
 446           if(it==0)
 447             {
 448               fProcessTracks = fSeeds;
 449               fSeeding = kTRUE;
 450             }
 451           else
 452             {
 453               fProcessTracks = fTracks;
 454               fSeeding = kFALSE;
 455             }
 456           if(!fProcessTracks)
 457             continue;
 458
 459           for(Int_t k=0; k<fProcessTracks->GetNTracks(); k++)
 460             {
 461               AliL3ModelTrack *track = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(k);
 462               if(!track) continue;
 463
 464               if(fSeeding)
 465                 if(track->GetClusterModel(i)->fSlice != fSlice) continue;
 466
 467               if(track->GetPadHit(i) < 0 || track->GetPadHit(i) > AliL3Transform::GetNPads(i)-1 ||
 468                  track->GetTimeHit(i) < 0 || track->GetTimeHit(i) > AliL3Transform::GetNTimeBins()-1)
 469                 {
 470                   track->SetCluster(i,0,0,0,0,0,0);
 471                   continue;
 472                 }
 473
 474               if(CheckCluster(k) == kFALSE)
 475                 fFailed++;
 476             }
 477         }
 478       AliL3MemHandler::UpdateRowPointer(rowPt);
 479     }
 480
 481   fSeeding = kTRUE;
 482   AddClusters();
 483   fSeeding = kFALSE;
 484   AddClusters();
 485
 486   cout<<"Fitted "<<fFitted<<" clusters, failed "<<fFailed<<endl;
 487   cout<<"Distribution:"<<endl;
 488   cout<<"Bad fit "<<fBadFitError<<endl;
 489   cout<<"Fit error "<<fFitError<<endl;
 490   cout<<"Result error "<<fResultError<<endl;
 491   cout<<"Fit range error "<<fFitRangeError<<endl;
 492
 493 }
 494
 495 Bool_t AliL3ClusterFitter::CheckCluster(Int_t trackindex)
 496 {
 497   //Check if this is a single or overlapping cluster
 498
 499   AliL3ModelTrack *track = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(trackindex);
 500
 501   Int_t row = fCurrentPadRow;
 502
 503   if(track->IsSet(row)) //A fit has already be performed on this one
 504     return kTRUE;
 505
 506   //Define the cluster region of this hit:
 507   Int_t padr[2]={999,-1};
 508   Int_t timer[2]={999,-1};
 509
 510   if(!SetFitRange(track,padr,timer))
 511     {
 512       track->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
 513
 514       if(fDebug)
 515         cout<<"Failed to fit cluster at row "<<row<<" pad "<<(Int_t)rint(track->GetPadHit(row))<<" time "
 516             <<(Int_t)rint(track->GetTimeHit(row))<<" hitcharge "
 517             <<fRow[(AliL3Transform::GetNTimeBins()+1)*(Int_t)rint(track->GetPadHit(row))+(Int_t)rint(track->GetTimeHit(row))].fCharge<<endl;
 518       fFitRangeError++;
 519       return kFALSE;
 520     }
 521
 522   //Check if any other track contributes to this cluster:
 523
 524   for(Int_t t=trackindex+1; t<fProcessTracks->GetNTracks(); t++)
 525     {
 526       AliL3ModelTrack *tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(t);
 527       if(!tr) continue;
 528       if(fSeeding)
 529         if(tr->GetClusterModel(row)->fSlice != fSlice) continue;
 530
 531       if(tr->GetPadHit(row) > padr[0] && tr->GetPadHit(row) < padr[1] &&
 532          tr->GetTimeHit(row) > timer[0] && tr->GetTimeHit(row) < timer[1])
 533         {
 534           if(SetFitRange(tr,padr,timer))
 535             track->SetOverlap(row,t);
 536         }
 537       /*
 538       Int_t xyw = (Int_t)ceil(sqrt(tr->GetParSigmaY2(row))*GetYWidthFactor());
 539       Int_t zw = (Int_t)ceil(sqrt(tr->GetParSigmaZ2(row))*GetZWidthFactor());
 540       if(
 541          (tr->GetPadHit(row) - xyw > padr[0] && tr->GetPadHit(row) - xyw < padr[1] &&
 542           tr->GetTimeHit(row) - zw > timer[0] && tr->GetTimeHit(row) - zw < timer[1]) ||
 543
 544          (tr->GetPadHit(row) + xyw > padr[0] && tr->GetPadHit(row) + xyw < padr[1] &&
 545           tr->GetTimeHit(row) - zw > timer[0] && tr->GetTimeHit(row) - zw < timer[1]) ||
 546
 547          (tr->GetPadHit(row) - xyw > padr[0] && tr->GetPadHit(row) - xyw < padr[1] &&
 548           tr->GetTimeHit(row) + zw > timer[0] && tr->GetTimeHit(row) + zw < timer[1]) ||
 549
 550          (tr->GetPadHit(row) + xyw > padr[0] && tr->GetPadHit(row) + xyw < padr[1] &&
 551           tr->GetTimeHit(row) + zw > timer[0] && tr->GetTimeHit(row) + zw < timer[1])
 552          )
 553         {
 554
 555           if(SetFitRange(tr,padr,timer)) //Expand the cluster fit range
 556             {
 557               track->SetOverlap(row,t);    //Set overlap
 558             }
 559         }
 560       */
 561     }
 562
 563   if(fDebug)
 564     cout<<"Fitting cluster with "<<track->GetNOverlaps(fCurrentPadRow)<<" overlaps"<<endl;
 565
 566   FitClusters(track,padr,timer);
 567   //CalculateWeightedMean(track,padr,timer);
 568   return kTRUE;
 569 }
 570
 571 Bool_t AliL3ClusterFitter::SetFitRange(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
 572 {
 573   Int_t row = fCurrentPadRow;
 574   Int_t nt = AliL3Transform::GetNTimeBins()+1;
 575
 576   Int_t nsearchbins=0;
 577   if(row < AliL3Transform::GetNRowLow())
 578     nsearchbins=25;
 579   else if(row < AliL3Transform::GetNRowLow() + AliL3Transform::GetNRowUp1())
 580     nsearchbins=25;
 581   else
 582     nsearchbins=25;
 583
 584   /*
 585   Int_t padloop[49] = {0,0,0,-1,1,-1,1,-1,1,0,0,-1,1,-1,1 ,2,-2,2,-2,2,-2,2,-2,2,-2
 586                        ,0,1,2,3,3,3,3,3,3,3
 587                        ,2,1,0,-1,-2,-3
 588                        ,-3,-3,-3,-3,-3,-3,-1,-1};
 589   Int_t timeloop[49] = {0,1,-1,0,0,1,1,-1,-1,2,-2,2,2,-2,-2 ,0,0,1,1,-1,-1,2,2,-2,-2
 590                         ,-3,-3,-3,-3,-2,-1,0,1,2,3
 591                         ,3,3,3,3,3,3,2,1,0,-1,-2,-3,-3,-3};
 592   */
 593
 594   Int_t padloop[49] = {0,0,0,-1,1,-1,1,-1,1,0,0,-1,1,-1,1 ,2,-2,2,-2,2,-2,2,-2,2,-2
 595                        ,-3,3,-3,3,-3,3,-3,3,-3,3
 596                        ,0,0,-1,-1,1,1,-2,-2,2,2,-3,-3,3,3};
 597   Int_t timeloop[49] = {0,1,-1,0,0,1,1,-1,-1,2,-2,2,2,-2,-2 ,0,0,1,1,-1,-1,2,2,-2,-2
 598                         ,0,0,-1,-1,1,1,-2,-2,2,2
 599                         ,-3,3,-3,3,-3,3,-3,3,-3,3,-3,3,-3,3};
 600
 601   Int_t padhit = (Int_t)rint(track->GetPadHit(row));
 602   Int_t timehit = (Int_t)rint(track->GetTimeHit(row));
 603   Int_t padmax=-1;
 604   Int_t timemax=-1;
 605
 606   for(Int_t index=0; index<nsearchbins; index++)
 607     {
 608       if(IsMaximum(padhit + padloop[index],timehit + timeloop[index]))
 609         {
 610           padmax = padhit + padloop[index];
 611           timemax = timehit + timeloop[index];
 612           break;
 613         }
 614     }
 615
 616
 617   //Define the cluster region of this hit:
 618   //The region we look for, is centered at the local maxima
 619   //and expanded around using the parametrized cluster width
 620   //according to track parameters.
 621
 622   Int_t xyw = 3;
 623   Int_t zw =  5;
 624
 625   if(padmax>=0 && timemax>=0)
 626     {
 627       //Set the hit to the local maxima of the cluster.
 628       //Store the maxima in the cluster model structure,
 629       //-only temporary, it will be overwritten when calling SetCluster.
 630
 631       track->GetClusterModel(row)->fDPad = padmax;
 632       track->GetClusterModel(row)->fDTime = timemax;
 633
 634       Int_t i=padmax,j=timemax;
 635       for(Int_t pdir=-1; pdir<=1; pdir+=2)
 636         {
 637           i=padmax;
 638           while(abs(padmax-i) < xyw)
 639             {
 640               Bool_t chargeonpad=kFALSE;
 641               for(Int_t tdir=-1; tdir<=1; tdir+=2)
 642                 {
 643                   j=timemax;
 644                   while(abs(timemax-j) < zw)
 645                     {
 646                       if(i<0 || i>=AliL3Transform::GetNPads(row) || j<0 || j>=AliL3Transform::GetNTimeBins()) break;
 647                       if(fRow[nt*i+j].fCharge)
 648                         {
 649                           if(i < padrange[0]) padrange[0]=i;
 650                           if(i > padrange[1]) padrange[1]=i;
 651                           if(j < timerange[0]) timerange[0]=j;
 652                           if(j > timerange[1]) timerange[1]=j;
 653                           chargeonpad=kTRUE;
 654                         }
 655                       else
 656                         break;
 657                       j+=tdir;
 658                     }
 659                 }
 660               if(!chargeonpad)
 661                 break;
 662               i+=pdir;
 663             }
 664         }
 665       /*
 666         for(Int_t i=padmax-xyw; i<=padmax+xyw; i++)
 667         {
 668         for(Int_t j=timemax-zw; j<=timemax+zw; j++)
 669         {
 670         if(i<0 || i>=AliL3Transform::GetNPads(row) || j<0 || j>=AliL3Transform::GetNTimeBins()) continue;
 671         if(fRow[nt*i+j].fCharge)
 672         {
 673         if(i < padrange[0]) padrange[0]=i;
 674         if(i > padrange[1]) padrange[1]=i;
 675         if(j < timerange[0]) timerange[0]=j;
 676         if(j > timerange[1]) timerange[1]=j;
 677         }
 678         }
 679         }
 680       */
 681
 682       if(fDebug)
 683         cout<<"New padrange "<<padrange[0]<<" "<<padrange[1]<<" "<<" time "<<timerange[0]<<" "<<timerange[1]<<endl;
 684       return kTRUE;
 685     }
 686   return kFALSE;
 687 }
 688
 689 Bool_t AliL3ClusterFitter::IsMaximum(Int_t pad,Int_t time)
 690 {
 691   if(pad<0 || pad >= AliL3Transform::GetNPads(fCurrentPadRow) ||
 692      time<0 || time >= AliL3Transform::GetNTimeBins())
 693     return kFALSE;
 694   Int_t nt = AliL3Transform::GetNTimeBins()+1;
 695   if(fRow[nt*pad+time].fUsed == kTRUE) return kFALSE; //Peak has been assigned before
 696   Int_t charge = fRow[nt*pad+time].fCharge;
 697   if(charge == 1024 || charge==0) return kFALSE;
 698   //if(charge == 0) return kFALSE;
 699
 700   //fRow[nt*pad+time].fUsed = kTRUE;
 701   //return kTRUE;
 702
 703   if(charge < fRow[nt*(pad-1)+(time-1)].fCharge) return kFALSE;
 704   if(charge < fRow[nt*(pad)+(time-1)].fCharge) return kFALSE;
 705   if(charge < fRow[nt*(pad+1)+(time-1)].fCharge) return kFALSE;
 706   if(charge < fRow[nt*(pad-1)+(time)].fCharge) return kFALSE;
 707   if(charge < fRow[nt*(pad+1)+(time)].fCharge) return kFALSE;
 708   if(charge < fRow[nt*(pad-1)+(time+1)].fCharge) return kFALSE;
 709   if(charge < fRow[nt*(pad)+(time+1)].fCharge) return kFALSE;
 710   if(charge < fRow[nt*(pad+1)+(time+1)].fCharge) return kFALSE;
 711   fRow[nt*pad+time].fUsed = kTRUE;
 712   return kTRUE;
 713 }
 714
 715 void AliL3ClusterFitter::CalculateWeightedMean(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
 716 {
 717   Float_t sum=0;
 718   Int_t npads=0;
 719   Float_t pad=0,time=0;
 720   Int_t nt = AliL3Transform::GetNTimeBins()+1;
 721   for(Int_t i=padrange[0]; i<=padrange[1]; i++)
 722     {
 723       Int_t lsum=0;
 724       for(Int_t j=timerange[0]; j<=timerange[1]; j++)
 725         {
 726           lsum += fRow[nt*(i-1)+(j-1)].fCharge;
 727           time += j*fRow[nt*(i-1)+(j-1)].fCharge;
 728         }
 729       if(lsum)
 730         npads++;
 731       pad += i * lsum;
 732     }
 733   if(sum)
 734     {
 735       pad /= sum;
 736       time /= sum;
 737       track->SetCluster(fCurrentPadRow,pad,time,sum,0,0,npads);
 738     }
 739   else
 740     track->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
 741 }
 742
 743 void AliL3ClusterFitter::FitClusters(AliL3ModelTrack *track,Int_t *padrange,Int_t *timerange)
 744 {
 745   //Handle single and overlapping clusters
 746
 747   /*
 748   if( (Int_t)rint(track->GetClusterModel(fCurrentPadRow)->fDPad) <= 1 ||
 749       (Int_t)rint(track->GetClusterModel(fCurrentPadRow)->fDPad) >= AliL3Transform::GetNPads(fCurrentPadRow)-2 ||
 750       (Int_t)rint(track->GetClusterModel(fCurrentPadRow)->fDTime) <= 1 ||
 751       (Int_t)rint(track->GetClusterModel(fCurrentPadRow)->fDTime) >= AliL3Transform::GetNTimeBins()-2)
 752     {
 753       CalculateWeightedMean(track,padrange,timerange);
 754       return;
 755     }
 756   */
 757   Int_t size = FIT_PTS;
 758   Int_t max_tracks = FIT_MAXPAR/NUM_PARS;
 759   if(track->GetNOverlaps(fCurrentPadRow) > max_tracks)
 760     {
 761       cerr<<"AliL3ClusterFitter::FitOverlappingClusters : Too many overlapping tracks"<<endl;
 762       return;
 763     }
 764   Int_t *overlaps = track->GetOverlaps(fCurrentPadRow);
 765
 766   //Check if at least one cluster is not already fitted
 767   Bool_t all_fitted=kTRUE;
 768
 769   Int_t k=-1;
 770   while(k < track->GetNOverlaps(fCurrentPadRow))
 771     {
 772       AliL3ModelTrack *tr=0;
 773       if(k==-1)
 774         tr = track;
 775       else
 776         tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
 777       k++;
 778       if(!tr) continue;
 779       if(!tr->IsSet(fCurrentPadRow) && !tr->IsPresent(fCurrentPadRow))//cluster has not been set and is not present
 780         {
 781           all_fitted = kFALSE;
 782           break;
 783         }
 784     }
 785   if(all_fitted)
 786     {
 787       if(fDebug)
 788         cout<<"But all the clusters were already fitted on row "<<fCurrentPadRow<<endl;
 789       return;
 790     }
 791
 792   //Allocate fit parameters array; this is interface to the C code
 793   plane = new DPOINT[FIT_PTS];
 794   memset(plane,0,FIT_PTS*sizeof(DPOINT));
 795
 796   Double_t x[FIT_PTS],y[FIT_PTS],s[FIT_PTS];
 797
 798   //Fill the fit parameters:
 799   Double_t a[FIT_MAXPAR];
 800   Int_t lista[FIT_MAXPAR];
 801   Double_t dev[FIT_MAXPAR],chisq_f;
 802
 803   Int_t fit_pars=0;
 804
 805   Int_t n_overlaps=0;
 806   k=-1;
 807
 808   //Fill the overlapping tracks:
 809   while(k < track->GetNOverlaps(fCurrentPadRow))
 810     {
 811       AliL3ModelTrack *tr=0;
 812       if(k==-1)
 813         tr = track;
 814       else
 815         tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
 816       k++;
 817       if(!tr) continue;
 818
 819       if(tr->IsSet(fCurrentPadRow) && !tr->IsPresent(fCurrentPadRow)) continue;//Cluster fit failed before
 820
 821       //Use the local maxima as the input to the fitting routine.
 822       //The local maxima is temporary stored in the cluster model:
 823       Int_t hitpad = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDPad);
 824       Int_t hittime = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDTime);
 825       Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*hitpad + hittime].fCharge;
 826
 827       if(fDebug)
 828         cout<<"Fitting track cluster, pad "<<tr->GetPadHit(fCurrentPadRow)<<" time "
 829             <<tr->GetTimeHit(fCurrentPadRow)<<" charge "<<charge<<" at local maxima in pad "<<hitpad
 830             <<" time "<<hittime<<" xywidth "<<sqrt(tr->GetParSigmaY2(fCurrentPadRow))
 831             <<" zwidth "<<sqrt(tr->GetParSigmaZ2(fCurrentPadRow))<<endl;
 832
 833       if(charge==0)
 834         {
 835           cerr<<"Charge still zero!"<<endl;
 836           exit(5);
 837         }
 838
 839       a[n_overlaps*NUM_PARS+2] = hitpad;
 840       a[n_overlaps*NUM_PARS+4] = hittime;
 841
 842       if(!tr->IsSet(fCurrentPadRow)) //Cluster is not fitted before
 843         {
 844           a[n_overlaps*NUM_PARS+1] = charge;
 845           a[n_overlaps*NUM_PARS+3] = sqrt(tr->GetParSigmaY2(fCurrentPadRow)) * GetYWidthFactor();
 846           a[n_overlaps*NUM_PARS+5] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
 847           //a[n_overlaps*NUM_PARS+6] = sqrt(tr->GetParSigmaZ2(fCurrentPadRow)) * GetZWidthFactor();
 848           lista[n_overlaps*NUM_PARS + 1] = 1;
 849           lista[n_overlaps*NUM_PARS + 2] = 1;
 850           lista[n_overlaps*NUM_PARS + 3] = 0;
 851           lista[n_overlaps*NUM_PARS + 4] = 1;
 852           lista[n_overlaps*NUM_PARS + 5] = 0;
 853           //lista[n_overlaps*NUM_PARS + 6] = 0;
 854           fit_pars             += 3;          //<-------------------
 855         }
 856       else  //Cluster was fitted before
 857         {
 858           if(!tr->IsPresent(fCurrentPadRow))
 859             {
 860               cerr<<"AliL3ClusterFitter::FindClusters : Cluster not present; there is a bug here"<<endl;
 861               exit(5);
 862             }
 863           Int_t charge;
 864           Float_t xywidth,zwidth,pad,time;
 865           tr->GetPad(fCurrentPadRow,pad);
 866           tr->GetTime(fCurrentPadRow,time);
 867           tr->GetClusterCharge(fCurrentPadRow,charge);
 868           xywidth = sqrt(tr->GetParSigmaY2(fCurrentPadRow));
 869           zwidth = sqrt(tr->GetParSigmaZ2(fCurrentPadRow));
 870           if(fDebug)
 871             cout<<endl<<"Cluster had been fitted before, pad "<<pad<<" time "<<time<<" charge "<<charge<<" width "<<xywidth<<" "<<zwidth<<endl;
 872
 873           a[n_overlaps*NUM_PARS+2] = pad;
 874           a[n_overlaps*NUM_PARS+4] = time;
 875           a[n_overlaps*NUM_PARS+1] = charge;
 876           a[n_overlaps*NUM_PARS+3] = xywidth * GetYWidthFactor();
 877           a[n_overlaps*NUM_PARS+5] = zwidth * GetZWidthFactor();
 878           //a[n_overlaps*NUM_PARS+6] = zwidth * GetZWidthFactor();
 879
 880           lista[n_overlaps*NUM_PARS + 1] = 1;
 881           lista[n_overlaps*NUM_PARS + 2] = 0;
 882           lista[n_overlaps*NUM_PARS + 3] = 0;
 883           lista[n_overlaps*NUM_PARS + 4] = 0;
 884           lista[n_overlaps*NUM_PARS + 5] = 0;
 885           //lista[n_overlaps*NUM_PARS + 6] = 0;
 886           fit_pars             += 1;
 887         }
 888       n_overlaps++;
 889     }
 890
 891   if(n_overlaps==0) //No clusters here
 892     {
 893       delete [] plane;
 894       return;
 895     }
 896
 897   Int_t pad_num=0;
 898   Int_t time_num_max=0;
 899   Int_t ndata=0;
 900   Int_t tot_charge=0;
 901   if(fDebug)
 902     cout<<"Padrange "<<padrange[0]<<" "<<padrange[1]<<" timerange "<<timerange[0]<<" "<<timerange[1]<<endl;
 903   for(Int_t i=padrange[0]; i<=padrange[1]; i++)
 904     {
 905       Int_t max_charge = 0;
 906       Int_t time_num=0;
 907       for(Int_t j=timerange[0]; j<=timerange[1]; j++)
 908         {
 909           Int_t charge = fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fCharge;
 910
 911           if(charge <= 0) continue;
 912
 913           time_num++;
 914           if(charge > max_charge)
 915             {
 916               max_charge = charge;
 917               //time_num++;
 918             }
 919           if(fDebug)
 920             cout<<"Filling padrow "<<fCurrentPadRow<<" pad "<<i<<" time "<<j<<" charge "<<charge<<endl;
 921           tot_charge += charge;
 922           ndata++;
 923           if(ndata >= size)
 924             {
 925               cerr<<"Too many points; row "<<fCurrentPadRow<<" padrange "<<padrange[0]<<" "<<padrange[1]<<" timerange "
 926                   <<timerange[0]<<" "<<timerange[1]<<endl;
 927               exit(5);
 928             }
 929
 930           plane[ndata].u = (Double_t)i;
 931           plane[ndata].v = (Double_t)j;
 932           x[ndata]=ndata;
 933           y[ndata]=charge;
 934           s[ndata]= 1 + sqrt((Double_t)charge);
 935         }
 936       if(max_charge) //there was charge on this pad
 937         pad_num++;
 938       if(time_num_max < time_num)
 939         time_num_max = time_num;
 940     }
 941
 942   if(pad_num <= 1 || time_num_max <=1 || n_overlaps > fNmaxOverlaps || ndata <= fit_pars) //too few to do fit
 943     {
 944       SetClusterfitFalse(track);
 945       if(fDebug)
 946         cout<<"Too few digits or too many overlaps: "<<pad_num<<" "<<time_num_max<<" "<<n_overlaps<<" ndata "<<ndata<<" fit_pars "<<fit_pars<<endl;
 947       delete [] plane;
 948       return;
 949     }
 950
 951
 952   Int_t npars = n_overlaps * NUM_PARS;
 953   if(fDebug)
 954     cout<<"Number of overlapping clusters "<<n_overlaps<<endl;
 955   Int_t ret = lev_marq_fit( x, y, s, ndata, a, lista, dev, npars, &chisq_f, f2gauss5 );
 956
 957   if(ret<0)
 958     {
 959       SetClusterfitFalse(track);
 960       fFailed++;
 961       fFitError++;
 962       delete [] plane;
 963       return;
 964       //exit(5);
 965     }
 966
 967   chisq_f /= (ndata-fit_pars);
 968   if(fDebug)
 969     cout<<"Chisq "<<chisq_f<<endl;
 970
 971   Bool_t overlapping=kFALSE;
 972   if(track->GetNOverlaps(fCurrentPadRow) > 0)//There is a overlap
 973     overlapping=kTRUE;
 974
 975   k=-1;
 976   n_overlaps=0;
 977   while(k < track->GetNOverlaps(fCurrentPadRow))
 978     {
 979       AliL3ModelTrack *tr=0;
 980       if(k==-1)
 981         tr = track;
 982       else
 983         tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[k]);
 984       k++;
 985       if(!tr) continue;
 986       if(!tr->IsPresent(fCurrentPadRow))
 987         {
 988           if(tr->IsSet(fCurrentPadRow)) continue;//This cluster has been set before
 989
 990           Int_t lpatch;
 991           if(fCurrentPadRow < AliL3Transform::GetNRowLow())
 992             lpatch=0;
 993           else if(fCurrentPadRow < AliL3Transform::GetNRowLow() + AliL3Transform::GetNRowUp1())
 994             lpatch=1;
 995           else
 996             lpatch=2;
 997
 998           //if(chisq_f < fChiSqMax[(Int_t)overlapping])//cluster fit is good enough
 999           if(chisq_f < fChiSqMax[lpatch])//cluster fit is good enough
1000             {
1001               tot_charge = (Int_t)(2*AliL3Transform::Pi() * a[n_overlaps*NUM_PARS+1] * a[n_overlaps*NUM_PARS+3] * a[n_overlaps*NUM_PARS+5]);
1002               Float_t fpad = a[n_overlaps*NUM_PARS+2];
1003               Float_t ftime = a[n_overlaps*NUM_PARS+4];
1004               if(tot_charge < 0 || fpad < -1 || fpad > AliL3Transform::GetNPads(fCurrentPadRow) ||
1005                  ftime < -1 || ftime > AliL3Transform::GetNTimeBins())
1006                 {
1007                   if(fDebug)
1008                     cout<<"AliL3ClusterFitter::Fatal result(s) in fit; in slice "<<fSlice<<" row "<<fCurrentPadRow
1009                         <<"; pad "<<fpad<<" time "<<ftime<<" charge "<<tot_charge<<" xywidth "<<a[n_overlaps*NUM_PARS+3]
1010                         <<" zwidth "<<a[n_overlaps*NUM_PARS+5]<<" peakcharge "<<a[n_overlaps*NUM_PARS+1]<<endl;
1011                   tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
1012                   fFailed++;
1013                   fResultError++;
1014                   continue;
1015                 }
1016
1017               tr->SetCluster(fCurrentPadRow,fpad,ftime,tot_charge,0,0,pad_num);
1018               /*
1019                 tr->SetCluster(fCurrentPadRow,fpad,ftime,tot_charge,
1020                 pow(a[n_overlaps*NUM_PARS+3],2),
1021                 pow(a[n_overlaps*NUM_PARS+5],2),pad_num);
1022               */
1023               if(fDebug)
1024                 {
1025                   cout<<"Setting cluster in slice "<<fSlice<<" row "<<fCurrentPadRow<<" pad "<<a[n_overlaps*NUM_PARS+2]<<" time "<<a[n_overlaps*NUM_PARS+4]
1026                       <<" padwidth "<<a[n_overlaps*NUM_PARS+3]<<" timewidth "<<a[n_overlaps*NUM_PARS+5]
1027                       <<" charge "<<tot_charge<<endl;
1028                 }
1029               /*
1030               //Set the digits to used:
1031               for(Int_t i=padrange[0]; i<=padrange[1]; i++)
1032               for(Int_t j=timerange[0]; j<=timerange[1]; j++)
1033               fRow[(AliL3Transform::GetNTimeBins()+1)*i + j].fUsed = kTRUE;
1034               */
1035               fFitted++;
1036             }
1037           else //fit was too bad
1038             {
1039               if(fDebug)
1040                 cout<<"Cluster fit was too bad"<<endl;
1041               tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
1042               fBadFitError++;
1043               fFailed++;
1044             }
1045         }
1046       n_overlaps++;
1047     }
1048
1049   delete [] plane;
1050 }
1051
1052 void AliL3ClusterFitter::SetClusterfitFalse(AliL3ModelTrack *track)
1053 {
1054   //Cluster fit failed, so set the clusters to all the participating
1055   //tracks to zero.
1056
1057   Int_t i=-1;
1058   Int_t *overlaps = track->GetOverlaps(fCurrentPadRow);
1059   while(i < track->GetNOverlaps(fCurrentPadRow))
1060     {
1061       AliL3ModelTrack *tr=0;
1062       if(i==-1)
1063         tr = track;
1064       else
1065         tr = (AliL3ModelTrack*)fProcessTracks->GetCheckedTrack(overlaps[i]);
1066       i++;
1067       if(!tr) continue;
1068
1069       //Set the digit data to unused, so it can be fitted to another bastard:
1070       Int_t hitpad = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDPad);
1071       Int_t hittime = (Int_t)rint(tr->GetClusterModel(fCurrentPadRow)->fDTime);
1072       fRow[(AliL3Transform::GetNTimeBins()+1)*hitpad + hittime].fUsed = kFALSE;
1073
1074       tr->SetCluster(fCurrentPadRow,0,0,0,0,0,0);
1075     }
1076
1077 }
1078
1079
1080 void AliL3ClusterFitter::AddClusters()
1081 {
1082   if(!fClusters)
1083     {
1084       fClusters = new AliL3SpacePointData[fNMaxClusters];
1085       fNClusters=0;
1086     }
1087
1088   if(fDebug)
1089     cout<<"Writing cluster in slice "<<fSlice<<" patch "<<fPatch<<endl;
1090
1091   AliL3TrackArray *tracks=0;
1092   if(fSeeding==kTRUE)
1093     tracks = fSeeds;
1094   else
1095     tracks = fTracks;
1096
1097   if(!tracks)
1098     return;
1099
1100   for(Int_t i=0; i<tracks->GetNTracks(); i++)
1101     {
1102       AliL3ModelTrack *tr = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
1103       if(!tr) continue;
1104
1105       UInt_t *hitids = tr->GetHitNumbers();
1106       Int_t nhits = tr->GetNHits();
1107       for(Int_t i=fRowMax; i>=fRowMin; i--)
1108         {
1109           if(fSeeding)
1110             if(tr->GetClusterModel(i)->fSlice != fSlice) continue;
1111           if(!tr->IsPresent(i)) continue;
1112           fCurrentPadRow = i;
1113           Float_t pad,time,xywidth,zwidth;
1114           Int_t charge;
1115           tr->GetPad(i,pad);
1116           tr->GetTime(i,time);
1117           tr->GetClusterCharge(i,charge);
1118
1119           if(pad < -1 || pad >= AliL3Transform::GetNPads(i) ||
1120              time < -1 || time >= AliL3Transform::GetNTimeBins())
1121             {
1122               continue;
1123               cout<<"slice "<<fSlice<<" row "<<i<<" pad "<<pad<<" time "<<time<<endl;
1124               tr->Print();
1125               exit(5);
1126             }
1127
1128           tr->CalculateClusterWidths(i,kTRUE); //Parametrize errors
1129
1130           tr->GetSigmaY2(i,xywidth);
1131           tr->GetSigmaZ2(i,zwidth);
1132           Float_t xyz[3];
1133           Int_t sector,row;
1134           AliL3Transform::Slice2Sector(fSlice,i,sector,row);
1135
1136           AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
1137
1138           if(fNClusters >= fNMaxClusters)
1139             {
1140               cerr<<"AliL3ClusterFitter::AddClusters : Too many clusters "<<fNClusters<<endl;
1141               exit(5);
1142             }
1143
1144           fClusters[fNClusters].fX = xyz[0];
1145           fClusters[fNClusters].fY = xyz[1];
1146           fClusters[fNClusters].fZ = xyz[2];
1147           fClusters[fNClusters].fCharge = charge;
1148           fClusters[fNClusters].fPadRow = i;
1149           Int_t pa = AliL3Transform::GetPatch(i);
1150           if(xywidth==0 || zwidth==0)
1151             cerr<<"AliL3ClusterFitter::AddClusters : Cluster with zero width"<<endl;
1152           if(xywidth>0)
1153             fClusters[fNClusters].fSigmaY2 = xywidth*pow(AliL3Transform::GetPadPitchWidth(pa),2);
1154           else
1155             fClusters[fNClusters].fSigmaY2 = 1;
1156           if(zwidth>0)
1157             fClusters[fNClusters].fSigmaZ2 = zwidth*pow(AliL3Transform::GetZWidth(),2);
1158           else
1159             fClusters[fNClusters].fSigmaZ2 = 1;
1160           Int_t pat=fPatch;
1161           if(fPatch==-1)
1162             pat=0;
1163           fClusters[fNClusters].fID = fNClusters + ((fSlice&0x7f)<<25)+((pat&0x7)<<22);
1164
1165           if(nhits >= AliL3Transform::GetNRows())
1166             {
1167               cerr<<"AliL3ClusterFitter::AddClusters : Cluster counter of out range "<<nhits<<endl;
1168               exit(5);
1169             }
1170           hitids[nhits++] = fClusters[fNClusters].fID;
1171
1172 #ifdef do_mc
1173           Int_t trackID[3];
1174           Int_t fpad = (Int_t)rint(pad);
1175           Int_t ftime = (Int_t)rint(time);
1176           if(fpad < 0)
1177             fpad=0;
1178           if(fpad >= AliL3Transform::GetNPads(i))
1179             fpad = AliL3Transform::GetNPads(i)-1;
1180           if(ftime<0)
1181             ftime=0;
1182           if(ftime >= AliL3Transform::GetNTimeBins())
1183             ftime = AliL3Transform::GetNTimeBins()-1;
1184           GetTrackID(fpad,ftime,trackID);
1185           fClusters[fNClusters].fTrackID[0] = trackID[0];
1186           fClusters[fNClusters].fTrackID[1] = trackID[1];
1187           fClusters[fNClusters].fTrackID[2] = trackID[2];
1188 #endif
1189           //cout<<"Setting id "<<trackID[0]<<" on pad "<<pad<<" time "<<time<<" row "<<i<<endl;
1190           fNClusters++;
1191         }
1192
1193       //Copy back the number of assigned clusters
1194       tr->SetNHits(nhits);
1195
1196     }
1197 }
1198
1199 void AliL3ClusterFitter::WriteTracks(Int_t min_hits)
1200 {
1201   if(!fSeeds)
1202     return;
1203
1204   AliL3TrackArray *fakes = new AliL3TrackArray();
1205
1206   Int_t clustercount=0;
1207   for(Int_t i=0; i<fSeeds->GetNTracks(); i++)
1208     {
1209       AliL3ModelTrack *tr = (AliL3ModelTrack*)fSeeds->GetCheckedTrack(i);
1210       if(!tr) continue;
1211       if(tr->GetNHits() < min_hits)
1212         {
1213           fakes->AddLast(tr);
1214           fSeeds->Remove(i);
1215         }
1216       clustercount += tr->GetNHits();
1217     }
1218
1219   cout<<"Writing "<<clustercount<<" clusters"<<endl;
1220   fSeeds->Compress();
1221   AliL3MemHandler mem;
1222   Char_t filename[1024];
1223   sprintf(filename,"%s/fitter/tracks_%d.raw",fPath,fEvent);
1224   mem.SetBinaryOutput(filename);
1225   mem.TrackArray2Binary(fSeeds);
1226   mem.CloseBinaryOutput();
1227
1228   //Write the fake tracks to its own file
1229   mem.Free();
1230   sprintf(filename,"%s/fitter/tracks_fakes_%d.raw",fPath,fEvent);
1231   mem.SetBinaryOutput(filename);
1232   mem.TrackArray2Binary(fakes);
1233   mem.CloseBinaryOutput();
1234   delete fakes;
1235 }
1236
1237 void AliL3ClusterFitter::WriteClusters(Bool_t global)
1238 {
1239   AliL3MemHandler mem;
1240   if(fDebug)
1241     cout<<"Write "<<fNClusters<<" clusters to file"<<endl;
1242   Char_t filename[1024];
1243   sprintf(filename,"%s/fitter/points_%d_%d_%d.raw",fPath,fEvent,fSlice,fPatch);
1244   mem.SetBinaryOutput(filename);
1245   if(global == kTRUE)
1246     mem.Transform(fNClusters,fClusters,fSlice);
1247   mem.Memory2Binary(fNClusters,fClusters);
1248   mem.CloseBinaryOutput();
1249   mem.Free();
1250
1251   delete [] fClusters;
1252   fClusters=0;
1253   fNClusters=0;
1254 }