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