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