HLT/hough/AliL3HoughTransformerRow.cxx

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