HLT/hough/AliL3HoughTransformerRow.cxx

   1 // @(#) $Id$
   2
   3
   4 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   5 //*-- Copyright &copy ALICE HLT Group
   6
   7 #include "AliL3StandardIncludes.h"
   8
   9 #include "AliL3Logging.h"
  10 #include "AliL3MemHandler.h"
  11 #include "AliL3Transform.h"
  12 #include "AliL3DigitData.h"
  13 #include "AliL3HistogramAdaptive.h"
  14 #include "AliL3HoughTrack.h"
  15 #include "AliL3HoughTransformerRow.h"
  16
  17 #if GCCVERSION == 3
  18 using namespace std;
  19 #endif
  20
  21 //_____________________________________________________________
  22 // AliL3HoughTransformerRow
  23 //
  24 // TPC rows Hough transformation class
  25 //
  26
  27 ClassImp(AliL3HoughTransformerRow)
  28
  29 Float_t AliL3HoughTransformerRow::fgBeta1 = 1.0/AliL3Transform::Row2X(79);
  30 Float_t AliL3HoughTransformerRow::fgBeta2 = 1.0/(AliL3Transform::Row2X(158)*(1.0+tan(AliL3Transform::Pi()*10/180)*tan(AliL3Transform::Pi()*10/180)));
  31
  32 AliL3HoughTransformerRow::AliL3HoughTransformerRow()
  33 {
  34   //Default constructor
  35   fParamSpace = 0;
  36
  37   fGapCount = 0;
  38   fCurrentRowCount = 0;
  39 #ifdef do_mc
  40   fTrackID = 0;
  41 #endif
  42   fTrackNRows = 0;
  43   fTrackFirstRow = 0;
  44   fTrackLastRow = 0;
  45   fInitialGapCount = 0;
  46
  47   fPrevBin = 0;
  48   fNextBin = 0;
  49   fNextRow = 0;
  50
  51   fStartPadParams = 0;
  52   fEndPadParams = 0;
  53   fLUTr2 = 0;
  54   fLUTforwardZ = 0;
  55   fLUTforwardZ2 = 0;
  56   fLUTbackwardZ = 0;
  57   fLUTbackwardZ2 = 0;
  58
  59 }
  60
  61 AliL3HoughTransformerRow::AliL3HoughTransformerRow(Int_t slice,Int_t patch,Int_t netasegments,Bool_t /*DoMC*/,Float_t zvertex) : AliL3HoughBaseTransformer(slice,patch,netasegments,zvertex)
  62 {
  63   //Normal constructor
  64   fParamSpace = 0;
  65
  66   fGapCount = 0;
  67   fCurrentRowCount = 0;
  68 #ifdef do_mc
  69   fTrackID = 0;
  70 #endif
  71
  72   fTrackNRows = 0;
  73   fTrackFirstRow = 0;
  74   fTrackLastRow = 0;
  75   fInitialGapCount = 0;
  76
  77   fPrevBin = 0;
  78   fNextBin = 0;
  79   fNextRow = 0;
  80
  81   fStartPadParams = 0;
  82   fEndPadParams = 0;
  83   fLUTr2 = 0;
  84   fLUTforwardZ = 0;
  85   fLUTforwardZ2 = 0;
  86   fLUTbackwardZ = 0;
  87   fLUTbackwardZ2 = 0;
  88
  89 }
  90
  91 AliL3HoughTransformerRow::~AliL3HoughTransformerRow()
  92 {
  93   //Destructor
  94   DeleteHistograms();
  95   if(fLastTransformer) return;
  96 #ifdef do_mc
  97   if(fTrackID)
  98     {
  99       for(Int_t i=0; i<GetNEtaSegments(); i++)
 100         {
 101           if(!fTrackID[i]) continue;
 102           delete fTrackID[i];
 103         }
 104       delete [] fTrackID;
 105       fTrackID = 0;
 106     }
 107 #endif
 108
 109   if(fGapCount)
 110     {
 111       for(Int_t i=0; i<GetNEtaSegments(); i++)
 112         {
 113           if(!fGapCount[i]) continue;
 114           delete [] fGapCount[i];
 115         }
 116       delete [] fGapCount;
 117       fGapCount = 0;
 118     }
 119   if(fCurrentRowCount)
 120     {
 121       for(Int_t i=0; i<GetNEtaSegments(); i++)
 122         {
 123           if(fCurrentRowCount[i])
 124             delete [] fCurrentRowCount[i];
 125         }
 126       delete [] fCurrentRowCount;
 127       fCurrentRowCount = 0;
 128     }
 129   if(fTrackNRows)
 130     {
 131       delete [] fTrackNRows;
 132       fTrackNRows = 0;
 133     }
 134   if(fTrackFirstRow)
 135     {
 136       delete [] fTrackFirstRow;
 137       fTrackFirstRow = 0;
 138     }
 139   if(fTrackLastRow)
 140     {
 141       delete [] fTrackLastRow;
 142       fTrackLastRow = 0;
 143     }
 144   if(fInitialGapCount)
 145     {
 146       delete [] fInitialGapCount;
 147       fInitialGapCount = 0;
 148     }
 149   if(fPrevBin)
 150     {
 151       for(Int_t i=0; i<GetNEtaSegments(); i++)
 152         {
 153           if(!fPrevBin[i]) continue;
 154           delete [] fPrevBin[i];
 155         }
 156       delete [] fPrevBin;
 157       fPrevBin = 0;
 158     }
 159   if(fNextBin)
 160     {
 161       for(Int_t i=0; i<GetNEtaSegments(); i++)
 162         {
 163           if(!fNextBin[i]) continue;
 164           delete [] fNextBin[i];
 165         }
 166       delete [] fNextBin;
 167       fNextBin = 0;
 168     }
 169   if(fNextRow)
 170     {
 171       for(Int_t i=0; i<GetNEtaSegments(); i++)
 172         {
 173           if(!fNextRow[i]) continue;
 174           delete [] fNextRow[i];
 175         }
 176       delete [] fNextRow;
 177       fNextRow = 0;
 178     }
 179   if(fStartPadParams)
 180     {
 181       for(Int_t i = AliL3Transform::GetFirstRow(0); i<=AliL3Transform::GetLastRow(5); i++)
 182         {
 183           if(!fStartPadParams[i]) continue;
 184           delete [] fStartPadParams[i];
 185         }
 186       delete [] fStartPadParams;
 187       fStartPadParams = 0;
 188     }
 189   if(fEndPadParams)
 190     {
 191       for(Int_t i = AliL3Transform::GetFirstRow(0); i<=AliL3Transform::GetLastRow(5); i++)
 192         {
 193           if(!fEndPadParams[i]) continue;
 194           delete [] fEndPadParams[i];
 195         }
 196       delete [] fEndPadParams;
 197       fEndPadParams = 0;
 198     }
 199   if(fLUTr2)
 200     {
 201       for(Int_t i = AliL3Transform::GetFirstRow(0); i<=AliL3Transform::GetLastRow(5); i++)
 202         {
 203           if(!fLUTr2[i]) continue;
 204           delete [] fLUTr2[i];
 205         }
 206       delete [] fLUTr2;
 207       fLUTr2 = 0;
 208     }
 209   if(fLUTforwardZ)
 210     {
 211       delete[] fLUTforwardZ;
 212       fLUTforwardZ=0;
 213     }
 214   if(fLUTforwardZ2)
 215     {
 216       delete[] fLUTforwardZ2;
 217       fLUTforwardZ2=0;
 218     }
 219   if(fLUTbackwardZ)
 220     {
 221       delete[] fLUTbackwardZ;
 222       fLUTbackwardZ=0;
 223     }
 224   if(fLUTbackwardZ2)
 225     {
 226       delete[] fLUTbackwardZ2;
 227       fLUTbackwardZ2=0;
 228     }
 229 }
 230
 231 void AliL3HoughTransformerRow::DeleteHistograms()
 232 {
 233   // Clean up
 234   if(!fParamSpace)
 235     return;
 236   for(Int_t i=0; i<GetNEtaSegments(); i++)
 237     {
 238       if(!fParamSpace[i]) continue;
 239       delete fParamSpace[i];
 240     }
 241   delete [] fParamSpace;
 242 }
 243
 244 struct AliL3TrackLength {
 245   Bool_t fIsFilled;
 246   UInt_t fFirstRow;
 247   UInt_t fLastRow;
 248   Float_t fTrackPt;
 249 };
 250
 251 void AliL3HoughTransformerRow::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
 252                                                 Int_t nybin,Float_t ymin,Float_t ymax)
 253 {
 254   //Create the histograms (parameter space)
 255   //nxbin = #bins in X
 256   //nybin = #bins in Y
 257   //xmin xmax ymin ymax = histogram limits in X and Y
 258   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
 259
 260   Char_t histname[256];
 261   for(Int_t i=0; i<GetNEtaSegments(); i++)
 262     {
 263       sprintf(histname,"paramspace_%d",i);
 264       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
 265     }
 266
 267   if(fLastTransformer) {
 268
 269     fGapCount = ((AliL3HoughTransformerRow *)fLastTransformer)->fGapCount;
 270     fCurrentRowCount = ((AliL3HoughTransformerRow *)fLastTransformer)->fCurrentRowCount;
 271 #ifdef do_mc
 272     fTrackID = ((AliL3HoughTransformerRow *)fLastTransformer)->fTrackID;
 273 #endif
 274     fTrackNRows = ((AliL3HoughTransformerRow *)fLastTransformer)->fTrackNRows;
 275     fTrackFirstRow = ((AliL3HoughTransformerRow *)fLastTransformer)->fTrackFirstRow;
 276     fTrackLastRow = ((AliL3HoughTransformerRow *)fLastTransformer)->fTrackLastRow;
 277     fInitialGapCount = ((AliL3HoughTransformerRow *)fLastTransformer)->fInitialGapCount;
 278
 279     fPrevBin = ((AliL3HoughTransformerRow *)fLastTransformer)->fPrevBin;
 280     fNextBin = ((AliL3HoughTransformerRow *)fLastTransformer)->fNextBin;
 281     fNextRow = ((AliL3HoughTransformerRow *)fLastTransformer)->fNextRow;
 282
 283     fStartPadParams = ((AliL3HoughTransformerRow *)fLastTransformer)->fStartPadParams;
 284     fEndPadParams = ((AliL3HoughTransformerRow *)fLastTransformer)->fEndPadParams;
 285     fLUTr2 = ((AliL3HoughTransformerRow *)fLastTransformer)->fLUTr2;
 286     fLUTforwardZ = ((AliL3HoughTransformerRow *)fLastTransformer)->fLUTforwardZ;
 287     fLUTforwardZ2 = ((AliL3HoughTransformerRow *)fLastTransformer)->fLUTforwardZ2;
 288     fLUTbackwardZ = ((AliL3HoughTransformerRow *)fLastTransformer)->fLUTbackwardZ;
 289     fLUTbackwardZ2 = ((AliL3HoughTransformerRow *)fLastTransformer)->fLUTbackwardZ2;
 290
 291     return;
 292   }
 293 #ifdef do_mc
 294   {
 295     AliL3Histogram *hist = fParamSpace[0];
 296     Int_t ncellsx = (hist->GetNbinsX()+3)/2;
 297     Int_t ncellsy = (hist->GetNbinsY()+3)/2;
 298     Int_t ncells = ncellsx*ncellsy;
 299     if(!fTrackID)
 300       {
 301         LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 302           <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(AliL3TrackIndex)<<" bytes to fTrackID"<<ENDLOG;
 303         fTrackID = new AliL3TrackIndex*[GetNEtaSegments()];
 304         for(Int_t i=0; i<GetNEtaSegments(); i++)
 305           fTrackID[i] = new AliL3TrackIndex[ncells];
 306       }
 307   }
 308 #endif
 309   AliL3Histogram *hist = fParamSpace[0];
 310   Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
 311   if(!fGapCount)
 312     {
 313       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 314         <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fGapCount"<<ENDLOG;
 315       fGapCount = new UChar_t*[GetNEtaSegments()];
 316       for(Int_t i=0; i<GetNEtaSegments(); i++)
 317         fGapCount[i] = new UChar_t[ncells];
 318     }
 319   if(!fCurrentRowCount)
 320     {
 321       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 322         <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fCurrentRowCount"<<ENDLOG;
 323       fCurrentRowCount = new UChar_t*[GetNEtaSegments()];
 324       for(Int_t i=0; i<GetNEtaSegments(); i++)
 325         fCurrentRowCount[i] = new UChar_t[ncells];
 326     }
 327   if(!fPrevBin)
 328     {
 329       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 330         <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fPrevBin"<<ENDLOG;
 331       fPrevBin = new UChar_t*[GetNEtaSegments()];
 332       for(Int_t i=0; i<GetNEtaSegments(); i++)
 333         fPrevBin[i] = new UChar_t[ncells];
 334     }
 335   if(!fNextBin)
 336     {
 337       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 338         <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fNextBin"<<ENDLOG;
 339       fNextBin = new UChar_t*[GetNEtaSegments()];
 340       for(Int_t i=0; i<GetNEtaSegments(); i++)
 341         fNextBin[i] = new UChar_t[ncells];
 342     }
 343   Int_t ncellsy = hist->GetNbinsY()+2;
 344   if(!fNextRow)
 345     {
 346       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 347         <<"Transformer: Allocating "<<GetNEtaSegments()*ncellsy*sizeof(UChar_t)<<" bytes to fNextRow"<<ENDLOG;
 348       fNextRow = new UChar_t*[GetNEtaSegments()];
 349       for(Int_t i=0; i<GetNEtaSegments(); i++)
 350         fNextRow[i] = new UChar_t[ncellsy];
 351     }
 352
 353   if(!fTrackNRows)
 354     {
 355       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 356         <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fTrackNRows"<<ENDLOG;
 357       fTrackNRows = new UChar_t[ncells];
 358       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 359         <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fTrackFirstRow"<<ENDLOG;
 360       fTrackFirstRow = new UChar_t[ncells];
 361       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 362         <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fTrackLastRow"<<ENDLOG;
 363       fTrackLastRow = new UChar_t[ncells];
 364       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 365         <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fInitialGapCount"<<ENDLOG;
 366       fInitialGapCount = new UChar_t[ncells];
 367
 368
 369       AliL3HoughTrack track;
 370       Int_t xmin = hist->GetFirstXbin();
 371       Int_t xmax = hist->GetLastXbin();
 372       Int_t xmiddle = (hist->GetNbinsX()+1)/2;
 373       Int_t ymin = hist->GetFirstYbin();
 374       Int_t ymax = hist->GetLastYbin();
 375       Int_t nxbins = hist->GetNbinsX()+2;
 376       Int_t nxgrid = (hist->GetNbinsX()+3)/2+1;
 377       Int_t nygrid = hist->GetNbinsY()+3;
 378
 379       AliL3TrackLength *tracklength = new AliL3TrackLength[nxgrid*nygrid];
 380       memset(tracklength,0,nxgrid*nygrid*sizeof(AliL3TrackLength));
 381
 382       for(Int_t ybin=ymin-1; ybin<=(ymax+1); ybin++)
 383         {
 384           for(Int_t xbin=xmin-1; xbin<=xmiddle; xbin++)
 385             {
 386               fTrackNRows[xbin + ybin*nxbins] = 255;
 387               for(Int_t deltay = 0; deltay <= 1; deltay++) {
 388                 for(Int_t deltax = 0; deltax <= 1; deltax++) {
 389
 390                   AliL3TrackLength *curtracklength = &tracklength[(xbin + deltax) + (ybin + deltay)*nxgrid];
 391                   UInt_t maxfirstrow = 0;
 392                   UInt_t maxlastrow = 0;
 393                   Float_t maxtrackpt = 0;
 394                   if(curtracklength->fIsFilled) {
 395                     maxfirstrow = curtracklength->fFirstRow;
 396                     maxlastrow = curtracklength->fLastRow;
 397                     maxtrackpt = curtracklength->fTrackPt;
 398                   }
 399                   else {
 400                     Float_t xtrack = hist->GetPreciseBinCenterX((Float_t)xbin+0.5*(Float_t)(2*deltax-1));
 401                     Float_t ytrack = hist->GetPreciseBinCenterY((Float_t)ybin+0.5*(Float_t)(2*deltay-1));
 402
 403                     Float_t psi = atan((xtrack-ytrack)/(fgBeta1-fgBeta2));
 404                     Float_t kappa = 2.0*(xtrack*cos(psi)-fgBeta1*sin(psi));
 405                     track.SetTrackParameters(kappa,psi,1);
 406                     Bool_t firstrow = kFALSE;
 407                     UInt_t curfirstrow = 0;
 408                     UInt_t curlastrow = 0;
 409
 410                     Double_t centerx = track.GetCenterX();
 411                     Double_t centery = track.GetCenterY();
 412                     Double_t radius = track.GetRadius();
 413
 414                     for(Int_t j=AliL3Transform::GetFirstRow(0); j<=AliL3Transform::GetLastRow(5); j++)
 415                       {
 416                         Float_t hit[3];
 417                         //                    if(!track.GetCrossingPoint(j,hit)) continue;
 418                         hit[0] = AliL3Transform::Row2X(j);
 419                         Double_t aa = (hit[0] - centerx)*(hit[0] - centerx);
 420                         Double_t r2 = radius*radius;
 421                         if(aa > r2)
 422                           continue;
 423
 424                         Double_t aa2 = sqrt(r2 - aa);
 425                         Double_t y1 = centery + aa2;
 426                         Double_t y2 = centery - aa2;
 427                         hit[1] = y1;
 428                         if(fabs(y2) < fabs(y1)) hit[1] = y2;
 429
 430                         hit[2] = 0;
 431
 432                         AliL3Transform::LocHLT2Raw(hit,0,j);
 433                         hit[1] += 0.5;
 434                         if(hit[1]>=0 && hit[1]<AliL3Transform::GetNPads(j))
 435                           {
 436                             if(!firstrow) {
 437                               curfirstrow = j;
 438                               firstrow = kTRUE;
 439                             }
 440                             curlastrow = j;
 441                           }
 442                         else {
 443                           if(firstrow) {
 444                             firstrow = kFALSE;
 445                             if((curlastrow-curfirstrow) >= (maxlastrow-maxfirstrow)) {
 446                               maxfirstrow = curfirstrow;
 447                               maxlastrow = curlastrow;
 448                             }
 449                           }
 450                         }
 451                       }
 452                     if((curlastrow-curfirstrow) >= (maxlastrow-maxfirstrow)) {
 453                       maxfirstrow = curfirstrow;
 454                       maxlastrow = curlastrow;
 455                     }
 456
 457                     maxtrackpt = track.GetPt();
 458
 459                     curtracklength->fIsFilled = kTRUE;
 460                     curtracklength->fFirstRow = maxfirstrow;
 461                     curtracklength->fLastRow = maxlastrow;
 462                     curtracklength->fTrackPt = maxtrackpt;
 463
 464                   }
 465                   if((maxlastrow-maxfirstrow) < fTrackNRows[xbin + ybin*nxbins]) {
 466                     fTrackNRows[xbin + ybin*nxbins] = maxlastrow-maxfirstrow;
 467                     fInitialGapCount[xbin + ybin*nxbins] = 1;
 468                     if((maxlastrow-maxfirstrow+1)<MIN_TRACK_LENGTH)
 469                       fInitialGapCount[xbin + ybin*nxbins] = MAX_N_GAPS+1;
 470                     if(maxtrackpt < 0.9*0.1*AliL3Transform::GetSolenoidField())
 471                       fInitialGapCount[xbin + ybin*nxbins] = MAX_N_GAPS;
 472                     fTrackFirstRow[xbin + ybin*nxbins] = maxfirstrow;
 473                     fTrackLastRow[xbin + ybin*nxbins] = maxlastrow;
 474
 475                     Int_t xmirror = xmax - xbin + 1;
 476                     Int_t ymirror = ymax - ybin + 1;
 477                     fTrackNRows[xmirror + ymirror*nxbins] = fTrackNRows[xbin + ybin*nxbins];
 478                     fInitialGapCount[xmirror + ymirror*nxbins] = fInitialGapCount[xbin + ybin*nxbins];
 479                     fTrackFirstRow[xmirror + ymirror*nxbins] = fTrackFirstRow[xbin + ybin*nxbins];
 480                     fTrackLastRow[xmirror + ymirror*nxbins] = fTrackLastRow[xbin + ybin*nxbins];
 481                   }
 482                 }
 483               }
 484               //              cout<<" fTrackNRows "<<(Int_t)fInitialGapCount[xbin + ybin*nxbins]<<" "<<xbin<<" "<<ybin<<" "<<(Int_t)fTrackNRows[xbin + ybin*nxbins]<<" "<<(Int_t)fTrackFirstRow[xbin + ybin*nxbins]<<" "<<(Int_t)fTrackLastRow[xbin + ybin*nxbins]<<" "<<endl;
 485             }
 486         }
 487       delete tracklength;
 488     }
 489
 490   if(!fStartPadParams)
 491     {
 492       Int_t nrows = AliL3Transform::GetLastRow(5) - AliL3Transform::GetFirstRow(0) + 1;
 493       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 494         <<"Transformer: Allocating about "<<nrows*100*sizeof(AliL3PadHoughParams)<<" bytes to fStartPadParams"<<ENDLOG;
 495       fStartPadParams = new AliL3PadHoughParams*[nrows];
 496       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 497         <<"Transformer: Allocating about "<<nrows*100*sizeof(AliL3PadHoughParams)<<" bytes to fEndPadParams"<<ENDLOG;
 498       fEndPadParams = new AliL3PadHoughParams*[nrows];
 499       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 500         <<"Transformer: Allocating about "<<nrows*100*sizeof(Float_t)<<" bytes to fLUTr2"<<ENDLOG;
 501       fLUTr2 = new Float_t*[nrows];
 502
 503       Float_t beta1 = fgBeta1;
 504       Float_t beta2 = fgBeta2;
 505       Float_t beta1minusbeta2 = fgBeta1 - fgBeta2;
 506       Float_t ymin = hist->GetYmin();
 507       Float_t histbin = hist->GetBinWidthY();
 508       Float_t xmin = hist->GetXmin();
 509       Float_t xmax = hist->GetXmax();
 510       Float_t xbin = (xmax-xmin)/hist->GetNbinsX();
 511       Int_t firstbinx = hist->GetFirstXbin();
 512       Int_t lastbinx = hist->GetLastXbin();
 513       Int_t nbinx = hist->GetNbinsX()+2;
 514       Int_t firstbin = hist->GetFirstYbin();
 515       Int_t lastbin = hist->GetLastYbin();
 516       for(Int_t i=AliL3Transform::GetFirstRow(0); i<=AliL3Transform::GetLastRow(5); i++)
 517         {
 518           Int_t npads = AliL3Transform::GetNPads(i);
 519           Int_t ipatch = AliL3Transform::GetPatch(i);
 520           Double_t padpitch = AliL3Transform::GetPadPitchWidth(ipatch);
 521           Float_t x = AliL3Transform::Row2X(i);
 522           Float_t x2 = x*x;
 523
 524           fStartPadParams[i] = new AliL3PadHoughParams[npads];
 525           fEndPadParams[i] = new AliL3PadHoughParams[npads];
 526           fLUTr2[i] = new Float_t[npads];
 527           for(Int_t pad=0; pad<npads; pad++)
 528             {
 529               Float_t y = (pad-0.5*(npads-1))*padpitch;
 530               fLUTr2[i][pad] = x2 + y*y;
 531               Float_t starty = (pad-0.5*npads)*padpitch;
 532               Float_t r1 = x2 + starty*starty;
 533               Float_t xoverr1 = x/r1;
 534               Float_t startyoverr1 = starty/r1;
 535               Float_t endy = (pad-0.5*(npads-2))*padpitch;
 536               Float_t r2 = x2 + endy*endy;
 537               Float_t xoverr2 = x/r2;
 538               Float_t endyoverr2 = endy/r2;
 539               Float_t a1 = beta1minusbeta2/(xoverr1-beta2);
 540               Float_t b1 = (xoverr1-beta1)/(xoverr1-beta2);
 541               Float_t a2 = beta1minusbeta2/(xoverr2-beta2);
 542               Float_t b2 = (xoverr2-beta1)/(xoverr2-beta2);
 543
 544               Float_t alpha1 = (a1*startyoverr1+b1*ymin-xmin)/xbin;
 545               Float_t deltaalpha1 = b1*histbin/xbin;
 546               if(b1<0)
 547                 alpha1 += deltaalpha1;
 548               Float_t alpha2 = (a2*endyoverr2+b2*ymin-xmin)/xbin;
 549               Float_t deltaalpha2 = b2*histbin/xbin;
 550               if(b2>=0)
 551                 alpha2 += deltaalpha2;
 552
 553               fStartPadParams[i][pad].fAlpha = alpha1;
 554               fStartPadParams[i][pad].fDeltaAlpha = deltaalpha1;
 555               fEndPadParams[i][pad].fAlpha = alpha2;
 556               fEndPadParams[i][pad].fDeltaAlpha = deltaalpha2;
 557
 558               //Find the first and last bin rows to be filled
 559               Bool_t binfound1 = kFALSE;
 560               Bool_t binfound2 = kFALSE;
 561               Int_t firstbin1 = lastbin;
 562               Int_t firstbin2 = lastbin;
 563               Int_t lastbin1 = firstbin;
 564               Int_t lastbin2 = firstbin;
 565               for(Int_t b=firstbin; b<=lastbin; b++, alpha1 += deltaalpha1, alpha2 += deltaalpha2)
 566                 {
 567                   Int_t binx1 = 1 + (Int_t)alpha1;
 568                   if(binx1<=lastbinx) {
 569                     UChar_t initialgapcount;
 570                     if(binx1>=firstbinx)
 571                       initialgapcount = fInitialGapCount[binx1 + b*nbinx];
 572                     else
 573                       initialgapcount = fInitialGapCount[firstbinx + b*nbinx];
 574                     if(initialgapcount != MAX_N_GAPS) {
 575                       if(!binfound1) {
 576                         firstbin1 = b;
 577                         binfound1 = kTRUE;
 578                       }
 579                       lastbin1 = b;
 580                     }
 581                   }
 582                   Int_t binx2 = 1 + (Int_t)alpha2;
 583                   if(binx2>=firstbin) {
 584                     UChar_t initialgapcount;
 585                     if(binx2<=lastbinx)
 586                       initialgapcount = fInitialGapCount[binx2 + b*nbinx];
 587                     else
 588                       initialgapcount = fInitialGapCount[lastbinx + b*nbinx];
 589                     if(initialgapcount != MAX_N_GAPS) {
 590                       if(!binfound2) {
 591                         firstbin2 = b;
 592                         binfound2 = kTRUE;
 593                       }
 594                       lastbin2 = b;
 595                     }
 596                   }
 597                 }
 598               fStartPadParams[i][pad].fFirstBin=firstbin1;
 599               fStartPadParams[i][pad].fLastBin=lastbin1;
 600               fEndPadParams[i][pad].fFirstBin=firstbin2;
 601               fEndPadParams[i][pad].fLastBin=lastbin2;
 602             }
 603         }
 604     }
 605
 606   //create lookup table for z of the digits
 607   if(!fLUTforwardZ)
 608     {
 609       Int_t ntimebins = AliL3Transform::GetNTimeBins();
 610       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 611         <<"Transformer: Allocating "<<ntimebins*sizeof(Float_t)<<" bytes to fLUTforwardZ"<<ENDLOG;
 612       fLUTforwardZ = new Float_t[ntimebins];
 613       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 614         <<"Transformer: Allocating "<<ntimebins*sizeof(Float_t)<<" bytes to fLUTforwardZ2"<<ENDLOG;
 615       fLUTforwardZ2 = new Float_t[ntimebins];
 616       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 617         <<"Transformer: Allocating "<<ntimebins*sizeof(Float_t)<<" bytes to fLUTbackwardZ"<<ENDLOG;
 618       fLUTbackwardZ = new Float_t[ntimebins];
 619       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 620         <<"Transformer: Allocating "<<ntimebins*sizeof(Float_t)<<" bytes to fLUTbackwardZ2"<<ENDLOG;
 621       fLUTbackwardZ2 = new Float_t[ntimebins];
 622       for(Int_t i=0; i<ntimebins; i++){
 623         Float_t z;
 624         z=AliL3Transform::GetZFast(0,i,GetZVertex());
 625         fLUTforwardZ[i]=z;
 626         fLUTforwardZ2[i]=z*z;
 627         z=AliL3Transform::GetZFast(18,i,GetZVertex());
 628         fLUTbackwardZ[i]=z;
 629         fLUTbackwardZ2[i]=z*z;
 630       }
 631     }
 632 }
 633
 634 void AliL3HoughTransformerRow::Reset()
 635 {
 636   //Reset all the histograms. Should be done when processing new slice
 637
 638   if(!fParamSpace)
 639     {
 640       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
 641         <<"No histograms to reset"<<ENDLOG;
 642       return;
 643     }
 644
 645   for(Int_t i=0; i<GetNEtaSegments(); i++)
 646     fParamSpace[i]->Reset();
 647
 648 #ifdef do_mc
 649   {
 650     AliL3Histogram *hist = fParamSpace[0];
 651     Int_t ncellsx = (hist->GetNbinsX()+3)/2;
 652     Int_t ncellsy = (hist->GetNbinsY()+3)/2;
 653     Int_t ncells = ncellsx*ncellsy;
 654     for(Int_t i=0; i<GetNEtaSegments(); i++)
 655       memset(fTrackID[i],0,ncells*sizeof(AliL3TrackIndex));
 656   }
 657 #endif
 658   AliL3Histogram *hist = fParamSpace[0];
 659   Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
 660   for(Int_t i=0; i<GetNEtaSegments(); i++)
 661     {
 662       memcpy(fGapCount[i],fInitialGapCount,ncells*sizeof(UChar_t));
 663       memcpy(fCurrentRowCount[i],fTrackFirstRow,ncells*sizeof(UChar_t));
 664     }
 665 }
 666
 667 Int_t AliL3HoughTransformerRow::GetEtaIndex(Double_t eta) const
 668 {
 669   //Return the histogram index of the corresponding eta.
 670
 671   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
 672   Double_t index = (eta-GetEtaMin())/etaslice;
 673   return (Int_t)index;
 674 }
 675
 676 inline AliL3Histogram *AliL3HoughTransformerRow::GetHistogram(Int_t etaindex)
 677 {
 678   // Return a pointer to the histogram which contains etaindex eta slice
 679   if(!fParamSpace || etaindex >= GetNEtaSegments() || etaindex < 0)
 680     return 0;
 681   if(!fParamSpace[etaindex])
 682     return 0;
 683   return fParamSpace[etaindex];
 684 }
 685
 686 Double_t AliL3HoughTransformerRow::GetEta(Int_t etaindex,Int_t /*slice*/) const
 687 {
 688   // Return eta calculated in the middle of the eta slice
 689   Double_t etaslice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
 690   Double_t eta=0;
 691   eta=(Double_t)((etaindex+0.5)*etaslice);
 692   return eta;
 693 }
 694
 695 void AliL3HoughTransformerRow::TransformCircle()
 696 {
 697   if(GetDataPointer())
 698     TransformCircleFromDigitArray();
 699   else if(fTPCRawStream)
 700     TransformCircleFromRawStream();
 701 }
 702 void AliL3HoughTransformerRow::TransformCircleFromDigitArray()
 703 {
 704   //Do the Hough Transform
 705
 706   Int_t netasegments = GetNEtaSegments();
 707   Double_t etamin = GetEtaMin();
 708   Double_t etaslice = (GetEtaMax() - etamin)/netasegments;
 709
 710   Int_t lowerthreshold = GetLowerThreshold();
 711
 712   //Assumes that all the etaslice histos are the same!
 713   AliL3Histogram *h = fParamSpace[0];
 714   Int_t firstbiny = h->GetFirstYbin();
 715   Int_t firstbinx = h->GetFirstXbin();
 716   Int_t lastbinx = h->GetLastXbin();
 717   Int_t nbinx = h->GetNbinsX()+2;
 718
 719   UChar_t lastpad;
 720   Int_t lastetaindex;
 721   AliL3EtaRow *etaclust = new AliL3EtaRow[netasegments];
 722
 723   AliL3DigitRowData *tempPt = GetDataPointer();
 724   if(!tempPt)
 725     {
 726       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
 727         <<"No input data "<<ENDLOG;
 728       return;
 729     }
 730
 731   Int_t ipatch = GetPatch();
 732   Int_t ilastpatch = GetLastPatch();
 733   Int_t islice = GetSlice();
 734   Float_t *lutz;
 735   Float_t *lutz2;
 736   if(islice < 18) {
 737     lutz = fLUTforwardZ;
 738     lutz2 = fLUTforwardZ2;
 739   }
 740   else {
 741     lutz = fLUTbackwardZ;
 742     lutz2 = fLUTbackwardZ2;
 743   }
 744
 745   //Loop over the padrows:
 746   for(UChar_t i=AliL3Transform::GetFirstRow(ipatch); i<=AliL3Transform::GetLastRow(ipatch); i++)
 747     {
 748       lastpad = 255;
 749       //Flush eta clusters array
 750       memset(etaclust,0,netasegments*sizeof(AliL3EtaRow));
 751
 752       Float_t x = AliL3Transform::Row2X((Int_t)i);
 753       Float_t x2 = x*x;
 754       Float_t radius2=0;
 755
 756       //Get the data on this padrow:
 757       AliL3DigitData *digPt = tempPt->fDigitData;
 758       if((Int_t)i != (Int_t)tempPt->fRow)
 759         {
 760           LOG(AliL3Log::kError,"AliL3HoughTransformerRow::TransformCircle","Data")
 761             <<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<(Int_t)i<<" "<<(Int_t)tempPt->fRow<<ENDLOG;
 762           continue;
 763         }
 764       //      cout<<" Starting row "<<i<<endl;
 765       //Loop over the data on this padrow:
 766       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 767         {
 768           UShort_t charge = digPt[j].fCharge;
 769           if((Int_t)charge <= lowerthreshold)
 770             continue;
 771           UChar_t pad = digPt[j].fPad;
 772           UShort_t time = digPt[j].fTime;
 773
 774           if(pad != lastpad)
 775             {
 776               radius2 = fLUTr2[(Int_t)i][(Int_t)pad];
 777               lastetaindex = -1;
 778             }
 779
 780           //Transform data to local cartesian coordinates:
 781           Float_t z = lutz[(Int_t)time];
 782           Float_t z2 = lutz2[(Int_t)time];
 783           // Acceptance cut : to be verified
 784           //      if(radius2<0.72406166*z2) continue;
 785           if(x2<0.8464*z2) continue;
 786           //Calculate the eta:
 787           Double_t r = sqrt(radius2+z2);
 788           Double_t eta = 0.5 * log((r+z)/(r-z));
 789           //Get the corresponding index, which determines which histogram to fill:
 790           Int_t etaindex = (Int_t)((eta-etamin)/etaslice);
 791
 792 #ifndef do_mc
 793           if(etaindex == lastetaindex) continue;
 794 #endif
 795           //      cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<etaindex<<endl;
 796
 797           if(etaindex < 0 || etaindex >= netasegments)
 798             continue;
 799
 800           if(!etaclust[etaindex].fIsFound)
 801             {
 802               etaclust[etaindex].fStartPad = pad;
 803               etaclust[etaindex].fEndPad = pad;
 804               etaclust[etaindex].fIsFound = 1;
 805 #ifdef do_mc
 806               FillClusterMCLabels(digPt[j],etaclust[etaindex]);
 807 #endif
 808               continue;
 809             }
 810           else
 811             {
 812               if(pad <= (etaclust[etaindex].fEndPad + 1))
 813                 {
 814                   etaclust[etaindex].fEndPad = pad;
 815 #ifdef do_mc
 816                   FillClusterMCLabels(digPt[j],etaclust[etaindex]);
 817 #endif
 818                 }
 819               else
 820                 {
 821                   FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
 822
 823                   etaclust[etaindex].fStartPad = pad;
 824                   etaclust[etaindex].fEndPad = pad;
 825
 826 #ifdef do_mc
 827                   memset(etaclust[etaindex].fMcLabels,0,MaxTrack);
 828                   FillClusterMCLabels(digPt[j],etaclust[etaindex]);
 829 #endif
 830                 }
 831             }
 832           lastpad = pad;
 833           lastetaindex = etaindex;
 834         }
 835       //Write remaining clusters
 836       for(Int_t etaindex = 0;etaindex < netasegments;etaindex++)
 837         {
 838           //Check for empty row
 839           if((etaclust[etaindex].fStartPad == 0) && (etaclust[etaindex].fEndPad == 0)) continue;
 840
 841           FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
 842         }
 843
 844       //Move the data pointer to the next padrow:
 845       AliL3MemHandler::UpdateRowPointer(tempPt);
 846     }
 847
 848   delete [] etaclust;
 849 }
 850
 851 void AliL3HoughTransformerRow::TransformCircleFromRawStream()
 852 {
 853   //Do the Hough Transform
 854
 855   Int_t netasegments = GetNEtaSegments();
 856   Double_t etamin = GetEtaMin();
 857   Double_t etaslice = (GetEtaMax() - etamin)/netasegments;
 858
 859   Int_t lowerthreshold = GetLowerThreshold();
 860
 861   //Assumes that all the etaslice histos are the same!
 862   AliL3Histogram *h = fParamSpace[0];
 863   Int_t firstbiny = h->GetFirstYbin();
 864   Int_t firstbinx = h->GetFirstXbin();
 865   Int_t lastbinx = h->GetLastXbin();
 866   Int_t nbinx = h->GetNbinsX()+2;
 867
 868   Int_t lastetaindex;
 869   AliL3EtaRow *etaclust = new AliL3EtaRow[netasegments];
 870
 871   if(!fTPCRawStream)
 872     {
 873       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
 874         <<"No input data "<<ENDLOG;
 875       return;
 876     }
 877
 878   Int_t ipatch = GetPatch();
 879   UChar_t rowmin = AliL3Transform::GetFirstRowOnDDL(ipatch);
 880   UChar_t rowmax = AliL3Transform::GetLastRowOnDDL(ipatch);
 881   //  Int_t ntimebins = AliL3Transform::GetNTimeBins();
 882   Int_t ilastpatch = GetLastPatch();
 883   Int_t islice = GetSlice();
 884   Float_t *lutz;
 885   Float_t *lutz2;
 886   if(islice < 18) {
 887     lutz = fLUTforwardZ;
 888     lutz2 = fLUTforwardZ2;
 889   }
 890   else {
 891     lutz = fLUTbackwardZ;
 892     lutz2 = fLUTbackwardZ2;
 893   }
 894
 895   //Flush eta clusters array
 896   memset(etaclust,0,netasegments*sizeof(AliL3EtaRow));
 897
 898   UChar_t i=0;
 899   Int_t npads=0;
 900   Float_t x=0;
 901   Float_t x2=0;
 902   Float_t radius2=0;
 903   UChar_t pad=0;
 904
 905   //Loop over the rawdata:
 906   while (fTPCRawStream->Next()) {
 907
 908     if(fTPCRawStream->IsNewSector() || fTPCRawStream->IsNewRow()) {
 909
 910       //Write remaining clusters
 911       for(Int_t etaindex = 0;etaindex < netasegments;etaindex++)
 912         {
 913           //Check for empty row
 914           if((etaclust[etaindex].fStartPad == 0) && (etaclust[etaindex].fEndPad == 0)) continue;
 915
 916           FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
 917         }
 918
 919       Int_t sector=fTPCRawStream->GetSector();
 920       Int_t row=fTPCRawStream->GetRow();
 921       Int_t slice,srow;
 922       AliL3Transform::Sector2Slice(slice,srow,sector,row);
 923       if(slice!=islice){
 924         LOG(AliL3Log::kError,"AliL3DDLDataFileHandler::DDLDigits2Memory","Slice")
 925           <<AliL3Log::kDec<<"Found slice "<<slice<<", expected "<<islice<<ENDLOG;
 926         continue;
 927       }
 928
 929       i=(UChar_t)srow;
 930       npads = AliL3Transform::GetNPads(srow)-1;
 931
 932       //Flush eta clusters array
 933       memset(etaclust,0,netasegments*sizeof(AliL3EtaRow));
 934
 935       x = AliL3Transform::Row2X((Int_t)i);
 936       x2 = x*x;
 937       radius2=0;
 938
 939     }
 940
 941     if((i<rowmin)||(i>rowmax))continue;
 942
 943     //    cout<<" Starting row "<<(UInt_t)i<<endl;
 944     //Loop over the data on this padrow:
 945     if(fTPCRawStream->IsNewRow() || fTPCRawStream->IsNewPad()) {
 946       pad=fTPCRawStream->GetPad();
 947       /*
 948       if((pad<0)||(pad>=(npads+1))){
 949         LOG(AliL3Log::kError,"AliL3DDLDataFileHandler::DDLDigits2Memory","Pad")
 950           <<AliL3Log::kDec<<"Pad value out of bounds "<<pad<<" "
 951           <<npads+1<<ENDLOG;
 952         continue;
 953       }
 954       */
 955       radius2 = fLUTr2[(Int_t)i][(Int_t)pad];
 956       lastetaindex = -1;
 957     }
 958
 959     UShort_t time=fTPCRawStream->GetTime();
 960     /*
 961     if((time<0)||(time>=ntimebins)){
 962       LOG(AliL3Log::kError,"AliL3DDLDataFileHandler::DDLDigits2Memory","Time")
 963         <<AliL3Log::kDec<<"Time out of bounds "<<time<<" "
 964         <<AliL3Transform::GetNTimeBins()<<ENDLOG;
 965       continue;
 966     }
 967     */
 968
 969     if(fTPCRawStream->GetSignal() <= lowerthreshold)
 970       continue;
 971
 972     //Transform data to local cartesian coordinates:
 973     Float_t z = lutz[(Int_t)time];
 974     Float_t z2 = lutz2[(Int_t)time];
 975     //    if(radius2<0.72406166*z2) continue;
 976     if(x2<0.8464*z2) continue;
 977     //Calculate the eta:
 978     Double_t r = sqrt(radius2+z2);
 979     Double_t eta = 0.5 * log((r+z)/(r-z));
 980     //Get the corresponding index, which determines which histogram to fill:
 981     Int_t etaindex = (Int_t)((eta-etamin)/etaslice);
 982
 983 #ifndef do_mc
 984     if(etaindex == lastetaindex) continue;
 985 #endif
 986     //    cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<etaindex<<endl;
 987
 988     if(etaindex < 0 || etaindex >= netasegments)
 989       continue;
 990
 991     if(!etaclust[etaindex].fIsFound)
 992       {
 993         etaclust[etaindex].fStartPad = pad;
 994         etaclust[etaindex].fEndPad = pad;
 995         etaclust[etaindex].fIsFound = 1;
 996         continue;
 997       }
 998     else
 999       {
1000         if(pad <= (etaclust[etaindex].fEndPad + 1))
1001           {
1002             etaclust[etaindex].fEndPad = pad;
1003           }
1004         else
1005           {
1006             FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
1007
1008             etaclust[etaindex].fStartPad = pad;
1009             etaclust[etaindex].fEndPad = pad;
1010
1011           }
1012       }
1013     lastetaindex = etaindex;
1014   }
1015
1016   //Write remaining clusters
1017   for(Int_t etaindex = 0;etaindex < netasegments;etaindex++)
1018     {
1019       //Check for empty row
1020       if((etaclust[etaindex].fStartPad == 0) && (etaclust[etaindex].fEndPad == 0)) continue;
1021
1022       FillCluster(i,etaindex,etaclust,ilastpatch,firstbinx,lastbinx,nbinx,firstbiny);
1023     }
1024
1025   delete [] etaclust;
1026 }
1027
1028 Int_t AliL3HoughTransformerRow::GetTrackID(Int_t etaindex,Double_t alpha1,Double_t alpha2) const
1029 {
1030   // Returns the MC label for a given peak found in the Hough space
1031 #ifndef do_mc
1032   LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID","Data")
1033     <<"Flag switched off"<<ENDLOG;
1034   return -1;
1035 #else
1036   if(etaindex < 0 || etaindex > GetNEtaSegments())
1037     {
1038       LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID","Data")
1039         <<"Wrong etaindex "<<etaindex<<ENDLOG;
1040       return -1;
1041     }
1042   AliL3Histogram *hist = fParamSpace[etaindex];
1043   Int_t bin = hist->FindLabelBin(alpha1,alpha2);
1044   if(bin==-1) {
1045     LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID()","")
1046       <<"Track candidate outside Hough space boundaries: "<<alpha1<<" "<<alpha2<<ENDLOG;
1047     return -1;
1048   }
1049   Int_t label=-1;
1050   Int_t max=0;
1051   for(UInt_t i=0; i<(MaxTrack-1); i++)
1052     {
1053       Int_t nhits=fTrackID[etaindex][bin].fNHits[i];
1054       if(nhits == 0) break;
1055       if(nhits > max)
1056         {
1057           max = nhits;
1058           label = fTrackID[etaindex][bin].fLabel[i];
1059         }
1060     }
1061   Int_t label2=-1;
1062   Int_t max2=0;
1063   for(UInt_t i=0; i<(MaxTrack-1); i++)
1064     {
1065       Int_t nhits=fTrackID[etaindex][bin].fNHits[i];
1066       if(nhits == 0) break;
1067       if(nhits > max2)
1068         {
1069           if(fTrackID[etaindex][bin].fLabel[i]!=label) {
1070             max2 = nhits;
1071             label2 = fTrackID[etaindex][bin].fLabel[i];
1072           }
1073         }
1074     }
1075   if(max2 !=0 ) {
1076     LOG(AliL3Log::kDebug,"AliL3HoughTransformerRow::GetTrackID()","")
1077       <<" TrackID"<<" label "<<label<<" max "<<max<<" label2 "<<label2<<" max2 "<<max2<<" "<<(Float_t)max2/(Float_t)max<<" "<<fTrackID[etaindex][bin].fLabel[MaxTrack-1]<<" "<<(Int_t)fTrackID[etaindex][bin].fNHits[MaxTrack-1]<<ENDLOG;
1078   }
1079   return label;
1080 #endif
1081   LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID()","")
1082     <<"Compile with do_mc flag!"<<ENDLOG;
1083   return -1;
1084 }
1085
1086 UChar_t *AliL3HoughTransformerRow::GetGapCount(Int_t etaindex)
1087 {
1088   return fGapCount[etaindex];
1089 }
1090
1091 UChar_t *AliL3HoughTransformerRow::GetCurrentRowCount(Int_t etaindex)
1092 {
1093   return fCurrentRowCount[etaindex];
1094 }
1095
1096 UChar_t *AliL3HoughTransformerRow::GetTrackNRows()
1097 {
1098   return fTrackNRows;
1099 }
1100
1101
1102 UChar_t *AliL3HoughTransformerRow::GetTrackFirstRow()
1103 {
1104   return fTrackFirstRow;
1105 }
1106
1107 UChar_t *AliL3HoughTransformerRow::GetTrackLastRow()
1108 {
1109   return fTrackLastRow;
1110 }
1111
1112 UChar_t *AliL3HoughTransformerRow::GetPrevBin(Int_t etaindex)
1113 {
1114   return fPrevBin[etaindex];
1115 }
1116
1117 UChar_t *AliL3HoughTransformerRow::GetNextBin(Int_t etaindex)
1118 {
1119   return fNextBin[etaindex];
1120 }
1121
1122 UChar_t *AliL3HoughTransformerRow::GetNextRow(Int_t etaindex)
1123 {
1124   return fNextRow[etaindex];
1125 }
1126
1127 inline void AliL3HoughTransformerRow::FillClusterRow(UChar_t i,Int_t binx1,Int_t binx2,UChar_t *ngaps2,UChar_t *currentrow2,UChar_t *lastrow2
1128 #ifdef do_mc
1129                            ,AliL3EtaRow etaclust,AliL3TrackIndex *trackid
1130 #endif
1131                            )
1132 {
1133   for(Int_t bin=binx1;bin<=binx2;bin++)
1134     {
1135       if(ngaps2[bin] < MAX_N_GAPS) {
1136         if(i < lastrow2[bin] && i > currentrow2[bin]) {
1137           ngaps2[bin] += (i-currentrow2[bin]-1);
1138           currentrow2[bin]=i;
1139         }
1140 #ifdef do_mc
1141         if(i < lastrow2[bin] && i >= currentrow2[bin]) {
1142           for(UInt_t t=0;t<(MaxTrack-1); t++)
1143             {
1144               Int_t label = etaclust.fMcLabels[t];
1145               if(label == 0) break;
1146               UInt_t c;
1147               Int_t tempbin2 = (Int_t)(bin/2);
1148               for(c=0; c<(MaxTrack-1); c++)
1149                 if(trackid[tempbin2].fLabel[c] == label || trackid[tempbin2].fNHits[c] == 0)
1150                   break;
1151               trackid[tempbin2].fLabel[c] = label;
1152               if(trackid[tempbin2].fCurrentRow[c] != i) {
1153                 trackid[tempbin2].fNHits[c]++;
1154                 trackid[tempbin2].fCurrentRow[c] = i;
1155               }
1156             }
1157         }
1158 #endif
1159       }
1160     }
1161
1162 }
1163
1164 inline void AliL3HoughTransformerRow::FillCluster(UChar_t i,Int_t etaindex,AliL3EtaRow *etaclust,Int_t ilastpatch,Int_t firstbinx,Int_t lastbinx,Int_t nbinx,Int_t firstbiny)
1165 {
1166   UChar_t *ngaps = fGapCount[etaindex];
1167   UChar_t *currentrow = fCurrentRowCount[etaindex];
1168   UChar_t *lastrow = fTrackLastRow;
1169   UChar_t *prevbin = fPrevBin[etaindex];
1170   UChar_t *nextbin = fNextBin[etaindex];
1171   UChar_t *nextrow = fNextRow[etaindex];
1172 #ifdef do_mc
1173   AliL3TrackIndex *trackid = fTrackID[etaindex];
1174 #endif
1175
1176   //Do the transformation:
1177   AliL3PadHoughParams *startparams = &fStartPadParams[(Int_t)i][etaclust[etaindex].fStartPad];
1178   AliL3PadHoughParams *endparams = &fEndPadParams[(Int_t)i][etaclust[etaindex].fEndPad];
1179
1180   Float_t alpha1 = startparams->fAlpha;
1181   Float_t deltaalpha1 = startparams->fDeltaAlpha;
1182   Float_t alpha2 = endparams->fAlpha;
1183   Float_t deltaalpha2 = endparams->fDeltaAlpha;
1184
1185   Int_t firstbin1 = startparams->fFirstBin;
1186   Int_t firstbin2 = endparams->fFirstBin;
1187   Int_t firstbin = firstbin1;
1188   if(firstbin>firstbin2) firstbin = firstbin2;
1189
1190   Int_t lastbin1 = startparams->fLastBin;
1191   Int_t lastbin2 = endparams->fLastBin;
1192   Int_t lastbin = lastbin1;
1193   if(lastbin<lastbin2) lastbin = lastbin2;
1194
1195   alpha1 += (firstbin-firstbiny)*deltaalpha1;
1196   alpha2 += (firstbin-firstbiny)*deltaalpha2;
1197
1198   //Fill the histogram along the alpha2 range
1199   if(ilastpatch == -1) {
1200     for(Int_t b=firstbin; b<=lastbin; b++, alpha1 += deltaalpha1, alpha2 += deltaalpha2)
1201       {
1202         Int_t binx1 = 1 + (Int_t)alpha1;
1203         if(binx1>lastbinx) continue;
1204         if(binx1<firstbinx) binx1 = firstbinx;
1205         Int_t binx2 = 1 + (Int_t)alpha2;
1206         if(binx2<firstbinx) continue;
1207         if(binx2>lastbinx) binx2 = lastbinx;
1208 #ifdef do_mc
1209         if(binx2<binx1) {
1210           LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::TransformCircle()","")
1211             <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<etaclust[etaindex].fStartPad<<" "<<etaclust[etaindex].fEndPad<<ENDLOG;
1212         }
1213 #endif
1214         Int_t tempbin = b*nbinx;
1215         UChar_t *ngaps2 = ngaps + tempbin;
1216         UChar_t *currentrow2 = currentrow + tempbin;
1217         UChar_t *lastrow2 = lastrow + tempbin;
1218 #ifdef do_mc
1219         Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2));
1220         AliL3TrackIndex *trackid2 = trackid + tempbin2;
1221 #endif
1222         FillClusterRow(i,binx1,binx2,ngaps2,currentrow2,lastrow2
1223 #ifdef do_mc
1224                        ,etaclust[etaindex],trackid2
1225 #endif
1226                        );
1227       }
1228   }
1229   else {
1230     for(Int_t b=firstbin; b<=lastbin; b++, alpha1 += deltaalpha1, alpha2 += deltaalpha2)
1231       {
1232         Int_t binx1 = 1 + (Int_t)alpha1;
1233         if(binx1>lastbinx) continue;
1234         if(binx1<firstbinx) binx1 = firstbinx;
1235         Int_t binx2 = 1 + (Int_t)alpha2;
1236         if(binx2<firstbinx) continue;
1237         if(binx2>lastbinx) binx2 = lastbinx;
1238 #ifdef do_mc
1239         if(binx2<binx1) {
1240           LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::TransformCircle()","")
1241             <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<etaclust[etaindex].fStartPad<<" "<<etaclust[etaindex].fEndPad<<ENDLOG;
1242         }
1243 #endif
1244         if(nextrow[b] > (b + 1)) {
1245           Int_t deltab = (nextrow[b] - b - 1);
1246           b += deltab;
1247           alpha1 += deltaalpha1*deltab;
1248           alpha2 += deltaalpha2*deltab;
1249           continue;
1250         }
1251         Int_t tempbin = b*nbinx;
1252         binx1 = (UInt_t)nextbin[binx1 + tempbin];
1253         binx2 = (UInt_t)prevbin[binx2 + tempbin];
1254         if(binx2<binx1) continue;
1255         UChar_t *ngaps2 = ngaps + tempbin;
1256         UChar_t *currentrow2 = currentrow + tempbin;
1257         UChar_t *lastrow2 = lastrow + tempbin;
1258 #ifdef do_mc
1259         Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2));
1260         AliL3TrackIndex *trackid2 = trackid + tempbin2;
1261 #endif
1262         FillClusterRow(i,binx1,binx2,ngaps2,currentrow2,lastrow2
1263 #ifdef do_mc
1264                        ,etaclust[etaindex],trackid2
1265 #endif
1266                        );
1267       }
1268   }
1269 }
1270
1271 #ifdef do_mc
1272 inline void AliL3HoughTransformerRow::FillClusterMCLabels(AliL3DigitData digpt,AliL3EtaRow etaclust)
1273 {
1274   for(Int_t t=0; t<3; t++)
1275     {
1276       Int_t label = digpt.fTrackID[t];
1277       if(label < 0) break;
1278       UInt_t c;
1279       for(c=0; c<(MaxTrack-1); c++)
1280         if(etaclust.fMcLabels[c] == label || etaclust.fMcLabels[c] == 0)
1281           break;
1282
1283       etaclust.fMcLabels[c] = label;
1284     }
1285 }
1286 #endif