457200b5f5ac678bf910ac0d68297c2819e051c6
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformerRow.cxx
1 // @(#) $Id$
2
3
4 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
5 //*-- Copyright &copy ALICE HLT Group
6
7 #include "AliL3StandardIncludes.h"
8
9 #include "AliL3Logging.h"
10 #include "AliL3MemHandler.h"
11 #include "AliL3Transform.h"
12 #include "AliL3DigitData.h"
13 #include "AliL3HistogramAdaptive.h"
14 #include "AliL3HoughTrack.h"
15 #include "AliL3HoughTransformerRow.h"
16
17 #if GCCVERSION == 3
18 using namespace std;
19 #endif
20
21 //_____________________________________________________________
22 // AliL3HoughTransformerRow
23 //
24 // TPC rows Hough transformation class
25 //
26
27 ClassImp(AliL3HoughTransformerRow)
28
29 UChar_t **AliL3HoughTransformerRow::fgRowCount = 0;
30 UChar_t **AliL3HoughTransformerRow::fgGapCount = 0;
31 UChar_t **AliL3HoughTransformerRow::fgCurrentRowCount = 0;
32 #ifdef do_mc
33 AliL3TrackIndex **AliL3HoughTransformerRow::fgTrackID = 0;
34 #endif
35 UChar_t *AliL3HoughTransformerRow::fgTrackNRows = 0;
36 UChar_t *AliL3HoughTransformerRow::fgTrackFirstRow = 0;
37 UChar_t *AliL3HoughTransformerRow::fgTrackLastRow = 0;
38
39 Float_t AliL3HoughTransformerRow::fgBeta1 = AliL3Transform::Row2X(79)/pow(AliL3Transform::Row2X(79),2);
40 Float_t AliL3HoughTransformerRow::fgBeta2 = (AliL3Transform::Row2X(158)+6.0)/pow((AliL3Transform::Row2X(158)+6.0),2);
41
42 AliL3HoughTransformerRow::AliL3HoughTransformerRow()
43 {
44   //Default constructor
45   fParamSpace = 0;
46   fDoMC = kFALSE;;
47
48   fLUTforwardZ=0;
49   fLUTforwardZ2=0;
50   fLUTbackwardZ=0;
51   fLUTbackwardZ2=0;
52 }
53
54 AliL3HoughTransformerRow::AliL3HoughTransformerRow(Int_t slice,Int_t patch,Int_t netasegments,Bool_t /*DoMC*/,Float_t zvertex) : AliL3HoughBaseTransformer(slice,patch,netasegments,zvertex)
55 {
56   //Normal constructor
57   fParamSpace = 0;
58   fDoMC = kFALSE;
59   fDoMC=kFALSE;
60 #ifdef do_mc
61   fDoMC = kTRUE;
62 #endif
63
64   fLUTforwardZ=0;
65   fLUTforwardZ2=0;
66   fLUTbackwardZ=0;
67   fLUTbackwardZ2=0;
68 }
69
70 AliL3HoughTransformerRow::~AliL3HoughTransformerRow()
71 {
72   //Destructor
73   DeleteHistograms();
74 #ifdef do_mc
75   if(fgTrackID)
76     {
77       for(Int_t i=0; i<GetNEtaSegments(); i++)
78         {
79           if(!fgTrackID[i]) continue;
80           delete fgTrackID[i];
81         }
82       delete [] fgTrackID;
83       fgTrackID = 0;
84     }
85 #endif
86
87   if(fgRowCount)
88     {
89       for(Int_t i=0; i<GetNEtaSegments(); i++)
90         {
91           if(!fgRowCount[i]) continue;
92           delete [] fgRowCount[i];
93         }
94       delete [] fgRowCount;
95       fgRowCount = 0;
96     }
97   if(fgGapCount)
98     {
99       for(Int_t i=0; i<GetNEtaSegments(); i++)
100         {
101           if(!fgGapCount[i]) continue;
102           delete [] fgGapCount[i];
103         }
104       delete [] fgGapCount;
105       fgGapCount = 0;
106     }
107   if(fgCurrentRowCount)
108     {
109       for(Int_t i=0; i<GetNEtaSegments(); i++)
110         {
111           if(fgCurrentRowCount[i])
112             delete [] fgCurrentRowCount[i];
113         }
114       delete [] fgCurrentRowCount;
115       fgCurrentRowCount = 0;
116     }
117   if(fgTrackNRows)
118     {
119       delete [] fgTrackNRows;
120       fgTrackNRows = 0;
121     }
122   if(fgTrackFirstRow)
123     {
124       delete [] fgTrackFirstRow;
125       fgTrackFirstRow = 0;
126     }
127   if(fgTrackLastRow)
128     {
129       delete [] fgTrackLastRow;
130       fgTrackLastRow = 0;
131     }
132 }
133
134 void AliL3HoughTransformerRow::DeleteHistograms()
135 {
136   // Clean up
137   delete[] fLUTforwardZ;
138   delete[] fLUTforwardZ2;
139   delete[] fLUTbackwardZ;
140   delete[] fLUTbackwardZ2;
141
142   fLUTforwardZ=0;
143   fLUTforwardZ2=0;
144   fLUTbackwardZ=0;
145   fLUTbackwardZ2=0;
146
147   if(!fParamSpace)
148     return;
149   for(Int_t i=0; i<GetNEtaSegments(); i++)
150     {
151       if(!fParamSpace[i]) continue;
152       delete fParamSpace[i];
153     }
154   delete [] fParamSpace;
155 }
156
157 void AliL3HoughTransformerRow::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
158                                                 Int_t nybin,Float_t ymin,Float_t ymax)
159 {
160   //Create the histograms (parameter space)  
161   //nxbin = #bins in X
162   //nybin = #bins in Y
163   //xmin xmax ymin ymax = histogram limits in X and Y
164   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
165   
166   Char_t histname[256];
167   for(Int_t i=0; i<GetNEtaSegments(); i++)
168     {
169       sprintf(histname,"paramspace_%d",i);
170       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
171     }
172   
173 #ifdef do_mc
174   if(fDoMC)
175     {
176       AliL3Histogram *hist = fParamSpace[0];
177       Int_t ncellsx = (hist->GetNbinsX()+3)/2;
178       Int_t ncellsy = (hist->GetNbinsY()+3)/2;
179       Int_t ncells = ncellsx*ncellsy;
180       if(!fgTrackID)
181         {
182           LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
183             <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(AliL3TrackIndex)<<" bytes to fgTrackID"<<ENDLOG;
184           fgTrackID = new AliL3TrackIndex*[GetNEtaSegments()];
185           for(Int_t i=0; i<GetNEtaSegments(); i++)
186             fgTrackID[i] = new AliL3TrackIndex[ncells];
187         }
188     }
189 #endif
190   AliL3Histogram *hist = fParamSpace[0];
191   Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
192   if(!fgRowCount)
193     {
194       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
195         <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fgRowCount"<<ENDLOG;
196       fgRowCount = new UChar_t*[GetNEtaSegments()];
197       for(Int_t i=0; i<GetNEtaSegments(); i++)
198         fgRowCount[i] = new UChar_t[ncells];
199     }
200   if(!fgGapCount)
201     {
202       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
203         <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fgGapCount"<<ENDLOG;
204       fgGapCount = new UChar_t*[GetNEtaSegments()];
205       for(Int_t i=0; i<GetNEtaSegments(); i++)
206         fgGapCount[i] = new UChar_t[ncells];
207     }
208   if(!fgCurrentRowCount)
209     {
210       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
211         <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fgCurrentRowCount"<<ENDLOG;
212       fgCurrentRowCount = new UChar_t*[GetNEtaSegments()];
213       for(Int_t i=0; i<GetNEtaSegments(); i++)
214         fgCurrentRowCount[i] = new UChar_t[ncells];
215     }
216
217   if(!fgTrackNRows)
218     {
219       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
220         <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fgTrackNRows"<<ENDLOG;
221       fgTrackNRows = new UChar_t[ncells];
222       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
223         <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fgTrackFirstRow"<<ENDLOG;
224       fgTrackFirstRow = new UChar_t[ncells];
225       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
226         <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fgTrackLastRow"<<ENDLOG;
227       fgTrackLastRow = new UChar_t[ncells];
228
229       AliL3HoughTrack track;
230       Int_t xmin = hist->GetFirstXbin();
231       Int_t xmax = hist->GetLastXbin();
232       Int_t ymin = hist->GetFirstYbin();
233       Int_t ymax = hist->GetLastYbin();
234       Int_t nxbins = hist->GetNbinsX()+2;
235       for(Int_t ybin=ymin; ybin<=ymax; ybin++)
236         {
237           for(Int_t xbin=xmin; xbin<=xmax; xbin++)
238             {
239               //cvetan: we get strange warning on gcc-2.95
240               //warning: large integer implicitly truncated to unsigned type
241               fgTrackNRows[xbin + ybin*nxbins] = 99999;
242               for(Int_t deltay = -1; deltay <= 1; deltay += 2) {
243                 for(Int_t deltax = -1; deltax <= 1; deltax += 2) {
244                 
245                   Float_t xtrack = hist->GetPreciseBinCenterX((Float_t)xbin+0.5*(Float_t)deltax);
246                   Float_t ytrack = hist->GetPreciseBinCenterY((Float_t)ybin+0.5*(Float_t)deltay);
247
248                   Float_t psi = atan((xtrack-ytrack)/(fgBeta1-fgBeta2));
249                   Float_t kappa = 2.0*(xtrack*cos(psi)-fgBeta1*sin(psi));
250                   track.SetTrackParameters(kappa,psi,1);
251                   Bool_t firstrow = kFALSE;
252                   UInt_t maxfirstrow = 0;
253                   UInt_t maxlastrow = 0;
254                   UInt_t curfirstrow = 0;
255                   UInt_t curlastrow = 0;
256                   for(Int_t j=AliL3Transform::GetFirstRow(0); j<=AliL3Transform::GetLastRow(5); j++)
257                     {
258                       Float_t hit[3];
259                       if(!track.GetCrossingPoint(j,hit)) continue;
260                       AliL3Transform::LocHLT2Raw(hit,0,j);
261                       if(hit[1]>=0 && hit[1]<AliL3Transform::GetNPads(j))
262                         {
263                           if(!firstrow) {
264                             curfirstrow = j;
265                             firstrow = kTRUE;
266                           }
267                           curlastrow = j;
268                         }
269                       else {
270                         if(firstrow) {
271                           firstrow = kFALSE;
272                           if((curlastrow-curfirstrow) >= (maxlastrow-maxfirstrow)) {
273                             maxfirstrow = curfirstrow;
274                             maxlastrow = curlastrow;
275                           }
276                         }
277                       }
278                     }
279                   if((curlastrow-curfirstrow) >= (maxlastrow-maxfirstrow)) {
280                     maxfirstrow = curfirstrow;
281                     maxlastrow = curlastrow;
282                   }
283                   if((maxlastrow-maxfirstrow) < fgTrackNRows[xbin + ybin*nxbins]) {
284                     fgTrackNRows[xbin + ybin*nxbins] = maxlastrow-maxfirstrow;
285                     fgTrackFirstRow[xbin + ybin*nxbins] = maxfirstrow;
286                     fgTrackLastRow[xbin + ybin*nxbins] = maxlastrow;
287                   }
288                 }
289               }
290               //              cout<<" fgTrackNRows "<<xbin<<" "<<ybin<<" "<<(Int_t)fgTrackNRows[xbin + ybin*nxbins]<<" "<<endl;
291             }
292         }
293     }
294
295   //create lookup table for z of the digits
296   Int_t ntimebins = AliL3Transform::GetNTimeBins();
297   fLUTforwardZ = new Float_t[ntimebins];
298   fLUTforwardZ2 = new Float_t[ntimebins];
299   fLUTbackwardZ = new Float_t[ntimebins];
300   fLUTbackwardZ2 = new Float_t[ntimebins];
301   for(Int_t i=0; i<ntimebins; i++){
302     Float_t z;
303     z=AliL3Transform::GetZFast(0,i,GetZVertex());
304     fLUTforwardZ[i]=z;
305     fLUTforwardZ2[i]=z*z;
306     z=AliL3Transform::GetZFast(18,i,GetZVertex());
307     fLUTbackwardZ[i]=z;
308     fLUTbackwardZ2[i]=z*z;
309   }
310
311 }
312
313 void AliL3HoughTransformerRow::Reset()
314 {
315   //Reset all the histograms. Should be done when processing new slice
316
317   if(!fParamSpace)
318     {
319       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
320         <<"No histograms to reset"<<ENDLOG;
321       return;
322     }
323   
324   for(Int_t i=0; i<GetNEtaSegments(); i++)
325     fParamSpace[i]->Reset();
326   
327 #ifdef do_mc
328   if(fDoMC)
329     {
330       AliL3Histogram *hist = fParamSpace[0];
331       Int_t ncellsx = (hist->GetNbinsX()+3)/2;
332       Int_t ncellsy = (hist->GetNbinsY()+3)/2;
333       Int_t ncells = ncellsx*ncellsy;
334       for(Int_t i=0; i<GetNEtaSegments(); i++)
335         memset(fgTrackID[i],0,ncells*sizeof(AliL3TrackIndex));
336     }
337 #endif
338   AliL3Histogram *hist = fParamSpace[0];
339   Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
340   for(Int_t i=0; i<GetNEtaSegments(); i++)
341     {
342       memset(fgRowCount[i],0,ncells*sizeof(UChar_t));
343       memset(fgGapCount[i],1,ncells*sizeof(UChar_t));
344       //      memset(fgCurrentRowCount[i],160,ncells*sizeof(UChar_t));
345       memcpy(fgCurrentRowCount[i],fgTrackFirstRow,ncells*sizeof(UChar_t));
346     }
347 }
348
349 Int_t AliL3HoughTransformerRow::GetEtaIndex(Double_t eta) const
350 {
351   //Return the histogram index of the corresponding eta. 
352
353   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
354   Double_t index = (eta-GetEtaMin())/etaslice;
355   return (Int_t)index;
356 }
357
358 inline AliL3Histogram *AliL3HoughTransformerRow::GetHistogram(Int_t etaindex)
359 {
360   // Return a pointer to the histogram which contains etaindex eta slice
361   if(!fParamSpace || etaindex >= GetNEtaSegments() || etaindex < 0)
362     return 0;
363   if(!fParamSpace[etaindex])
364     return 0;
365   return fParamSpace[etaindex];
366 }
367
368 Double_t AliL3HoughTransformerRow::GetEta(Int_t etaindex,Int_t /*slice*/) const
369 {
370   // Return eta calculated in the middle of the eta slice
371   Double_t etaslice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
372   Double_t eta=0;
373   eta=(Double_t)((etaindex+0.5)*etaslice);
374   return eta;
375 }
376
377 struct AliL3EtaRow {
378   UChar_t fStartPad; //First pad in the cluster
379   UChar_t fEndPad; //Last pad in the cluster
380   Bool_t fIsFound; //Is the cluster already found
381   Float_t fStartY; //Y position of the first pad in the cluster
382 #ifdef do_mc
383   Int_t fMcLabels[MaxTrack]; //Array to store mc labels inside cluster
384 #endif
385 };
386
387 void AliL3HoughTransformerRow::TransformCircle()
388 {
389   //Do the Hough Transform
390   Float_t beta1 = fgBeta1;
391   Float_t beta2 = fgBeta2;
392   Float_t beta1minusbeta2 = fgBeta1 - fgBeta2;
393
394   Int_t netasegments = GetNEtaSegments();
395   Double_t etamin = GetEtaMin();
396   Double_t etaslice = (GetEtaMax() - etamin)/netasegments;
397
398   Int_t lowerthreshold = GetLowerThreshold();
399
400   //Assumes that all the etaslice histos are the same!
401   AliL3Histogram *h = fParamSpace[0];
402   Float_t ymin = h->GetYmin();
403   //Float_t y_max = h->GetYmax();
404   Float_t histbin = h->GetBinWidthY();
405   Int_t firstbin = h->GetFirstYbin();
406   Int_t lastbin = h->GetLastYbin();
407   Float_t xmin = h->GetXmin();
408   Float_t xmax = h->GetXmax();
409   Float_t xbin = (xmax-xmin)/h->GetNbinsX();
410   Int_t firstbinx = h->GetFirstXbin()-1;
411   Int_t lastbinx = h->GetLastXbin()+1;
412   Int_t nbinx = h->GetNbinsX()+2;
413
414   UChar_t lastpad;
415   Int_t lastetaindex;
416   AliL3EtaRow *etaclust = new AliL3EtaRow[netasegments];
417
418   AliL3DigitRowData *tempPt = GetDataPointer();
419   if(!tempPt)
420     {
421       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
422         <<"No input data "<<ENDLOG;
423       return;
424     }
425
426   Int_t ipatch = GetPatch();
427   Double_t padpitch = AliL3Transform::GetPadPitchWidth(ipatch);
428   Int_t islice = GetSlice();
429   Float_t *fLUTz;
430   Float_t *fLUTz2;
431   if(islice < 18) {
432     fLUTz = fLUTforwardZ;
433     fLUTz2 = fLUTforwardZ2;
434   }
435   else {
436     fLUTz = fLUTbackwardZ;
437     fLUTz2 = fLUTbackwardZ2;
438   }
439
440   //Loop over the padrows:
441   for(UChar_t i=AliL3Transform::GetFirstRow(ipatch); i<=AliL3Transform::GetLastRow(ipatch); i++)
442     {
443       Int_t npads = AliL3Transform::GetNPads((Int_t)i)-1;
444
445       lastpad = 255;
446       //Flush eta clusters array
447       memset(etaclust,0,netasegments*sizeof(AliL3EtaRow));  
448
449       Float_t x = AliL3Transform::Row2X((Int_t)i);
450       Float_t x2 = x*x;
451       Float_t y=0,r2=0;
452
453       //Get the data on this padrow:
454       AliL3DigitData *digPt = tempPt->fDigitData;
455       if((Int_t)i != (Int_t)tempPt->fRow)
456         {
457           LOG(AliL3Log::kError,"AliL3HoughTransformerRow::TransformCircle","Data")
458             <<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<(Int_t)i<<" "<<(Int_t)tempPt->fRow<<ENDLOG;
459           continue;
460         }
461       //      cout<<" Starting row "<<i<<endl;
462       //Loop over the data on this padrow:
463       for(UInt_t j=0; j<tempPt->fNDigit; j++)
464         {
465           UShort_t charge = digPt[j].fCharge;
466           if((Int_t)charge <= lowerthreshold)
467             continue;
468           UChar_t pad = digPt[j].fPad;
469           UShort_t time = digPt[j].fTime;
470
471           if(pad != lastpad)
472             {
473               y = (pad-0.5*npads)*padpitch;
474               Float_t y2 = y*y;
475               r2 = x2 + y2;
476               lastetaindex = -1;
477             }
478
479           //Transform data to local cartesian coordinates:
480           Float_t z = fLUTz[(Int_t)time];
481           Float_t z2 = fLUTz2[(Int_t)time];
482           //Calculate the eta:
483           Double_t r = sqrt(r2+z2);
484           Double_t eta = 0.5 * log((r+z)/(r-z));
485           //Get the corresponding index, which determines which histogram to fill:
486           Int_t etaindex = (Int_t)((eta-etamin)/etaslice);
487
488 #ifndef do_mc
489           if(etaindex == lastetaindex) continue;
490 #endif
491           //      cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<etaindex<<endl;
492           
493           if(etaindex < 0 || etaindex >= netasegments)
494             continue;
495
496           if(!etaclust[etaindex].fIsFound)
497             {
498               etaclust[etaindex].fStartPad = pad;
499               etaclust[etaindex].fEndPad = pad;
500               etaclust[etaindex].fStartY = y - padpitch/2.0;
501               etaclust[etaindex].fIsFound = 1;
502 #ifdef do_mc
503               if(fDoMC)
504                 {
505                   for(Int_t t=0; t<3; t++)
506                     {
507                       Int_t label = digPt[j].fTrackID[t];
508                       if(label < 0) break;
509                       UInt_t c;
510                       for(c=0; c<(MaxTrack-1); c++)
511                         if(etaclust[etaindex].fMcLabels[c] == label || etaclust[etaindex].fMcLabels[c] == 0)
512                           break;
513                       //                      if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
514                       etaclust[etaindex].fMcLabels[c] = label;
515                     }
516                 }
517 #endif
518               continue;
519             }
520           else
521             {
522               if(pad <= (etaclust[etaindex].fEndPad + 1))
523                 {
524                   etaclust[etaindex].fEndPad = pad;
525 #ifdef do_mc
526                   if(fDoMC)
527                     {
528                       for(Int_t t=0; t<3; t++)
529                         {
530                           Int_t label = digPt[j].fTrackID[t];
531                           if(label < 0) break;
532                           UInt_t c;
533                           for(c=0; c<(MaxTrack-1); c++)
534                             if(etaclust[etaindex].fMcLabels[c] == label || etaclust[etaindex].fMcLabels[c] == 0)
535                               break;
536                           //                      if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
537                           etaclust[etaindex].fMcLabels[c] = label;
538                         }
539                     }
540 #endif
541                 }
542               else
543                 {
544                   Bool_t fillcluster = kTRUE;
545                   if(fillcluster) {
546
547                   UChar_t *nrows = fgRowCount[etaindex];
548                   UChar_t *ngaps = fgGapCount[etaindex];
549                   UChar_t *currentrow = fgCurrentRowCount[etaindex];
550                   UChar_t *lastrow = fgTrackLastRow;
551
552                   //Do the transformation:
553                   Float_t starty = etaclust[etaindex].fStartY;
554                   if(etaclust[etaindex].fStartPad == 0)
555                     starty -= 0.0*padpitch;
556                   Float_t r1 = x2 + starty*starty;
557                   Float_t xoverr1 = x/r1;
558                   Float_t startyoverr1 = starty/r1;
559                   Float_t endy = (etaclust[etaindex].fEndPad-0.5*(npads-1))*padpitch;
560                   if(etaclust[etaindex].fEndPad == npads)
561                     endy += 0.0*padpitch;
562                   Float_t r2 = x2 + endy*endy; 
563                   Float_t xoverr2 = x/r2;
564                   Float_t endyoverr2 = endy/r2;
565                   Float_t a1 = beta1minusbeta2/(xoverr1-beta2);
566                   Float_t b1 = (xoverr1-beta1)/(xoverr1-beta2);
567                   Float_t a2 = beta1minusbeta2/(xoverr2-beta2);
568                   Float_t b2 = (xoverr2-beta1)/(xoverr2-beta2);
569
570                   Float_t alpha1 = (a1*startyoverr1+b1*ymin-xmin)/xbin;
571                   Float_t deltaalpha1 = b1*histbin/xbin;
572                   if(b1<0)
573                     alpha1 += deltaalpha1;
574                   Float_t alpha2 = (a2*endyoverr2+b2*ymin-xmin)/xbin;
575                   Float_t deltaalpha2 = b2*histbin/xbin;
576                   if(b2>=0)
577                     alpha2 += deltaalpha2;
578
579                   //Fill the histogram along the phirange
580                   for(Int_t b=firstbin; b<=lastbin; b++, alpha1 += deltaalpha1, alpha2 += deltaalpha2)
581                     {
582                       Int_t binx1 = 1 + (Int_t)alpha1;
583                       if(binx1>lastbinx) continue;
584                       if(binx1<firstbinx) binx1 = firstbinx;
585                       Int_t binx2 = 1 + (Int_t)alpha2;
586                       if(binx2<firstbinx) continue;
587                       if(binx2>lastbinx) binx2 = lastbinx;
588 #ifdef do_mc
589                       if(binx2<binx1) {
590                         LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::TransformCircle()","")
591                           <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<x<<" "<<starty<<" "<<endy<<ENDLOG;
592                       }
593 #endif
594                       Int_t tempbin = b*nbinx;
595                       UChar_t *nrows2 = nrows + tempbin;
596                       UChar_t *ngaps2 = ngaps + tempbin;
597                       UChar_t *currentrow2 = currentrow + tempbin;
598                       UChar_t *lastrow2 = lastrow + tempbin;
599                       for(Int_t bin=binx1;bin<=binx2;bin++)
600                         {
601                           if(ngaps2[bin] < MAX_N_GAPS) {
602                             if(i > (currentrow2[bin] + MAX_GAP_SIZE) && i < lastrow2[bin]) {
603                               ngaps2[bin] = MAX_N_GAPS;
604                               continue;
605                             }
606                             if(i > currentrow2[bin] && i < lastrow2[bin])
607                               {
608                                 nrows2[bin]++;
609                                 if(i > (currentrow2[bin] + 1))
610                                   ngaps2[bin]++;
611                                 currentrow2[bin]=i;
612                               }
613 #ifdef do_mc
614                             if(fDoMC)
615                               {
616                                 for(UInt_t t=0;t<(MaxTrack-1); t++)
617                                   {
618                                     Int_t label = etaclust[etaindex].fMcLabels[t];
619                                     if(label == 0) break;
620                                     UInt_t c;
621                                     Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2)) + (Int_t)(bin/2);
622                                     for(c=0; c<(MaxTrack-1); c++)
623                                       if(fgTrackID[etaindex][tempbin2].fLabel[c] == label || fgTrackID[etaindex][tempbin2].fNHits[c] == 0)
624                                         break;
625                                     //                              if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
626                                     fgTrackID[etaindex][tempbin2].fLabel[c] = label;
627                                     if(fgTrackID[etaindex][tempbin2].fCurrentRow[c] != i) {
628                                       fgTrackID[etaindex][tempbin2].fNHits[c]++;
629                                       fgTrackID[etaindex][tempbin2].fCurrentRow[c] = i;
630                                     }
631                                   }
632                               }
633 #endif
634                           }
635                         }
636                     }
637                   //End of the transformation
638
639                   }
640
641                   etaclust[etaindex].fStartPad = pad;
642                   etaclust[etaindex].fEndPad = pad;
643                   etaclust[etaindex].fStartY = y - padpitch/2.0;
644
645 #ifdef do_mc
646                   if(fDoMC)
647                     {
648                       memset(etaclust[etaindex].fMcLabels,0,MaxTrack);
649                       for(Int_t t=0; t<3; t++)
650                         {
651                           Int_t label = digPt[j].fTrackID[t];
652                           if(label < 0) break;
653                           UInt_t c;
654                           for(c=0; c<(MaxTrack-1); c++)
655                             if(etaclust[etaindex].fMcLabels[c] == label || etaclust[etaindex].fMcLabels[c] == 0)
656                               break;
657                           //                      if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
658                           etaclust[etaindex].fMcLabels[c] = label;
659                         }
660                     }
661 #endif
662                 }
663             }
664           lastpad = pad;
665           lastetaindex = etaindex;
666         }
667       //Write remaining clusters
668       for(Int_t etaindex = 0;etaindex < netasegments;etaindex++)
669         {
670           //Check for empty row
671           if((etaclust[etaindex].fStartPad == 0) && (etaclust[etaindex].fEndPad == 0)) continue;
672
673           UChar_t *nrows = fgRowCount[etaindex];
674           UChar_t *ngaps = fgGapCount[etaindex];
675           UChar_t *currentrow = fgCurrentRowCount[etaindex];
676           UChar_t *lastrow = fgTrackLastRow;
677
678           //Do the transformation:
679           Float_t starty = etaclust[etaindex].fStartY;
680           if(etaclust[etaindex].fStartPad == 0)
681             starty -= 0.0*padpitch;
682           Float_t r1 = x2 + starty*starty; 
683           Float_t xoverr1 = x/r1;
684           Float_t startyoverr1 = starty/r1;
685           Float_t endy = (etaclust[etaindex].fEndPad-0.5*(npads-1))*padpitch;
686           if(etaclust[etaindex].fEndPad == npads)
687             endy += 0.0*padpitch;
688           Float_t r2 = x2 + endy*endy; 
689           Float_t xoverr2 = x/r2;
690           Float_t endyoverr2 = endy/r2;
691           Float_t a1 = beta1minusbeta2/(xoverr1-beta2);
692           Float_t b1 = (xoverr1-beta1)/(xoverr1-beta2);
693           Float_t a2 = beta1minusbeta2/(xoverr2-beta2);
694           Float_t b2 = (xoverr2-beta1)/(xoverr2-beta2);
695
696           Float_t alpha1 = (a1*startyoverr1+b1*ymin-xmin)/xbin;
697           Float_t deltaalpha1 = b1*histbin/xbin;
698           if(b1<0)
699             alpha1 += deltaalpha1;
700           Float_t alpha2 = (a2*endyoverr2+b2*ymin-xmin)/xbin;
701           Float_t deltaalpha2 = b2*histbin/xbin;
702           if(b2>=0)
703             alpha2 += deltaalpha2;
704
705           //Fill the histogram along the phirange
706           for(Int_t b=firstbin; b<=lastbin; b++, alpha1 += deltaalpha1, alpha2 += deltaalpha2)
707             {
708               Int_t binx1 = 1 + (Int_t)alpha1;
709               if(binx1>lastbinx) continue;
710               if(binx1<firstbinx) binx1 = firstbinx;
711               Int_t binx2 = 1 + (Int_t)alpha2;
712               if(binx2<firstbinx) continue;
713               if(binx2>lastbinx) binx2 = lastbinx;
714 #ifdef do_mc
715               if(binx2<binx1) {
716                 LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::TransformCircle()","")
717                   <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<x<<" "<<starty<<" "<<endy<<ENDLOG;
718               }
719 #endif
720               Int_t tempbin = b*nbinx;
721               UChar_t *nrows2 = nrows + tempbin;
722               UChar_t *ngaps2 = ngaps + tempbin;
723               UChar_t *currentrow2 = currentrow + tempbin;
724               UChar_t *lastrow2 = lastrow + tempbin;
725               for(Int_t bin=binx1;bin<=binx2;bin++)
726                 {
727                   if(ngaps2[bin] < MAX_N_GAPS) {
728                     if(i > (currentrow2[bin] + MAX_GAP_SIZE) && i < lastrow2[bin]) {
729                       ngaps2[bin] = MAX_N_GAPS;
730                       continue;
731                     }
732                     if(i > currentrow2[bin] && i < lastrow2[bin])
733                       {
734                         nrows2[bin]++;
735                         if(i > (currentrow2[bin] + 1))
736                           ngaps2[bin]++;
737                         currentrow2[bin]=i;
738                       }
739 #ifdef do_mc
740                     if(fDoMC)
741                       {
742                         for(UInt_t t=0;t<(MaxTrack-1); t++)
743                           {
744                             Int_t label = etaclust[etaindex].fMcLabels[t];
745                             if(label == 0) break;
746                             UInt_t c;
747                             Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2)) + (Int_t)(bin/2);
748                             for(c=0; c<(MaxTrack-1); c++)
749                               if(fgTrackID[etaindex][tempbin2].fLabel[c] == label || fgTrackID[etaindex][tempbin2].fNHits[c] == 0)
750                                 break;
751                             //                      if(c == MaxTrack-1) cout<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
752                             fgTrackID[etaindex][tempbin2].fLabel[c] = label;
753                             if(fgTrackID[etaindex][tempbin2].fCurrentRow[c] != i) {
754                               fgTrackID[etaindex][tempbin2].fNHits[c]++;
755                               fgTrackID[etaindex][tempbin2].fCurrentRow[c] = i;
756                             }
757                           }
758                       }
759 #endif
760                   }
761                 }
762             }
763           //End of the transformation
764
765         }
766
767       //Move the data pointer to the next padrow:
768       AliL3MemHandler::UpdateRowPointer(tempPt);
769     }
770
771   delete [] etaclust;
772 }
773
774 Int_t AliL3HoughTransformerRow::GetTrackID(Int_t etaindex,Double_t alpha1,Double_t alpha2) const
775 {
776   // Returns the MC label for a given peak found in the Hough space
777   if(!fDoMC)
778     {
779       LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID","Data")
780         <<"Flag switched off"<<ENDLOG;
781       return -1;
782     }
783   
784 #ifdef do_mc
785   if(etaindex < 0 || etaindex > GetNEtaSegments())
786     {
787       LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID","Data")
788         <<"Wrong etaindex "<<etaindex<<ENDLOG;
789       return -1;
790     }
791   AliL3Histogram *hist = fParamSpace[etaindex];
792   Int_t bin = hist->FindLabelBin(alpha1,alpha2);
793   if(bin==-1) {
794     LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID()","")
795       <<"Track candidate outside Hough space boundaries: "<<alpha1<<" "<<alpha2<<ENDLOG;
796     return -1;
797   }
798   Int_t label=-1;
799   Int_t max=0;
800   for(UInt_t i=0; i<(MaxTrack-1); i++)
801     {
802       Int_t nhits=fgTrackID[etaindex][bin].fNHits[i];
803       if(nhits == 0) break;
804       if(nhits > max)
805         {
806           max = nhits;
807           label = fgTrackID[etaindex][bin].fLabel[i];
808         }
809     }
810   Int_t label2=-1;
811   Int_t max2=0;
812   for(UInt_t i=0; i<(MaxTrack-1); i++)
813     {
814       Int_t nhits=fgTrackID[etaindex][bin].fNHits[i];
815       if(nhits == 0) break;
816       if(nhits > max2)
817         {
818           if(fgTrackID[etaindex][bin].fLabel[i]!=label) {
819             max2 = nhits;
820             label2 = fgTrackID[etaindex][bin].fLabel[i];
821           }
822         }
823     }
824   if(max2 !=0 ) {
825     LOG(AliL3Log::kDebug,"AliL3HoughTransformerRow::GetTrackID()","")
826       <<" TrackID"<<" label "<<label<<" max "<<max<<" label2 "<<label2<<" max2 "<<max2<<" "<<(Float_t)max2/(Float_t)max<<" "<<fgTrackID[etaindex][bin].fLabel[MaxTrack-1]<<" "<<(Int_t)fgTrackID[etaindex][bin].fNHits[MaxTrack-1]<<ENDLOG;
827   }
828   return label;
829 #endif
830   LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID()","")
831     <<"Compile with do_mc flag!"<<ENDLOG;
832   return -1;
833 }
834
835 UChar_t *AliL3HoughTransformerRow::GetRowCount(Int_t etaindex)
836 {
837   return fgRowCount[etaindex];
838 }
839
840 UChar_t *AliL3HoughTransformerRow::GetGapCount(Int_t etaindex)
841 {
842   return fgGapCount[etaindex];
843 }
844
845 UChar_t *AliL3HoughTransformerRow::GetCurrentRowCount(Int_t etaindex)
846 {
847   return fgCurrentRowCount[etaindex];
848 }
849
850 UChar_t *AliL3HoughTransformerRow::GetTrackNRows()
851 {
852   return fgTrackNRows;
853 }
854
855
856 UChar_t *AliL3HoughTransformerRow::GetTrackFirstRow()
857 {
858   return fgTrackFirstRow;
859 }
860
861 UChar_t *AliL3HoughTransformerRow::GetTrackLastRow()
862 {
863   return fgTrackLastRow;
864 }