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