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