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