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