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