Merged Cvetans RowHoughTransformer, Anders latest developments in comp
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformerRow.cxx
1 // @(#) $Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ALICE HLT Group
5
6 #include "AliL3StandardIncludes.h"
7
8 #include "AliL3Logging.h"
9 #include "AliL3MemHandler.h"
10 #include "AliL3Transform.h"
11 #include "AliL3DigitData.h"
12 #include "AliL3HistogramAdaptive.h"
13 #include "AliL3HoughTransformerRow.h"
14
15 #if __GNUC__ == 3
16 using namespace std;
17 #endif
18
19 //_____________________________________________________________
20 // AliL3HoughTransformerRow
21 //
22 // TPC rows Hough transformation class
23 //
24
25 ClassImp(AliL3HoughTransformerRow)
26
27 UChar_t **AliL3HoughTransformerRow::fRowCount = 0;
28 UChar_t **AliL3HoughTransformerRow::fGapCount = 0;
29 UChar_t **AliL3HoughTransformerRow::fCurrentRowCount = 0;
30 #ifdef do_mc
31 TrackIndex **AliL3HoughTransformerRow::fTrackID = 0;
32 #endif
33
34
35 AliL3HoughTransformerRow::AliL3HoughTransformerRow()
36 {
37   //Default constructor
38   fParamSpace = 0;
39   fDoMC = kFALSE;;
40
41   fLUT2sinphi0up=0;
42   fLUT2sinphi0low=0;
43   fLUT2cosphi0up=0;
44   fLUT2cosphi0low=0;
45
46   fLUTforwardZ=0;
47   fLUTforwardZ2=0;
48   fLUTbackwardZ=0;
49   fLUTbackwardZ2=0;
50 }
51
52 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)
53 {
54   //Normal constructor
55   fParamSpace = 0;
56   fDoMC = kFALSE;
57   fDoMC=kFALSE;
58 #ifdef do_mc
59   fDoMC = kTRUE;
60 #endif
61
62   fLUT2sinphi0up=0;
63   fLUT2sinphi0low=0;
64   fLUT2cosphi0up=0;
65   fLUT2cosphi0low=0;
66
67   fLUTforwardZ=0;
68   fLUTforwardZ2=0;
69   fLUTbackwardZ=0;
70   fLUTbackwardZ2=0;
71 }
72
73 AliL3HoughTransformerRow::~AliL3HoughTransformerRow()
74 {
75   DeleteHistograms();
76 #ifdef do_mc
77   if(fTrackID)
78     {
79       for(Int_t i=0; i<GetNEtaSegments(); i++)
80         {
81           if(!fTrackID[i]) continue;
82           delete fTrackID[i];
83         }
84       delete [] fTrackID;
85       fTrackID = 0;
86     }
87 #endif
88
89   if(fRowCount)
90     {
91       for(Int_t i=0; i<GetNEtaSegments(); i++)
92         {
93           if(!fRowCount[i]) continue;
94           delete [] fRowCount[i];
95         }
96       delete [] fRowCount;
97       fRowCount = 0;
98     }
99   if(fGapCount)
100     {
101       for(Int_t i=0; i<GetNEtaSegments(); i++)
102         {
103           if(!fGapCount[i]) continue;
104           delete [] fGapCount[i];
105         }
106       delete [] fGapCount;
107       fGapCount = 0;
108     }
109   if(fCurrentRowCount)
110     {
111       for(Int_t i=0; i<GetNEtaSegments(); i++)
112         {
113           if(fCurrentRowCount[i])
114             delete [] fCurrentRowCount[i];
115         }
116       delete [] fCurrentRowCount;
117       fCurrentRowCount = 0;
118     }
119 }
120
121 void AliL3HoughTransformerRow::DeleteHistograms()
122 {
123   delete[] fLUT2sinphi0up;
124   delete[] fLUT2cosphi0up;
125   delete[] fLUT2sinphi0low;
126   delete[] fLUT2cosphi0low;
127
128   fLUT2sinphi0up=0;
129   fLUT2cosphi0up=0;
130   fLUT2sinphi0low=0;
131   fLUT2cosphi0low=0;
132
133   delete[] fLUTforwardZ;
134   delete[] fLUTforwardZ2;
135   delete[] fLUTbackwardZ;
136   delete[] fLUTbackwardZ2;
137
138   fLUTforwardZ=0;
139   fLUTforwardZ2=0;
140   fLUTbackwardZ=0;
141   fLUTbackwardZ2=0;
142
143   if(!fParamSpace)
144     return;
145   for(Int_t i=0; i<GetNEtaSegments(); i++)
146     {
147       if(!fParamSpace[i]) continue;
148       delete fParamSpace[i];
149     }
150   delete [] fParamSpace;
151 }
152
153 void AliL3HoughTransformerRow::CreateHistograms(Int_t nxbin,Float_t pt_min,
154                                                 Int_t nybin,Float_t phimin,Float_t phimax)
155 {
156   //Create the histograms (parameter space).
157   //These are 2D histograms, span by kappa (curvature of track) 
158   //and phi0 (emission angle with x-axis).
159   //The arguments give the range and binning; 
160   //nxbin = #bins in kappa
161   //nybin = #bins in phi0
162   //pt_min = mimium Pt of track (corresponding to maximum kappa)
163   //phi_min = mimimum phi0 
164   //phi_max = maximum phi0 
165     
166   Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
167   //Double_t torad = AliL3Transform::Pi()/180;
168   CreateHistograms(nxbin,-1.*x,x,nybin,phimin/**torad*/,phimax/**torad*/);
169 }
170
171 void AliL3HoughTransformerRow::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
172                                                 Int_t nybin,Float_t ymin,Float_t ymax)
173 {
174   
175   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
176   
177   Char_t histname[256];
178   for(Int_t i=0; i<GetNEtaSegments(); i++)
179     {
180       sprintf(histname,"paramspace_%d",i);
181       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
182     }
183   
184 #ifdef do_mc
185   if(fDoMC)
186     {
187       AliL3Histogram *hist = fParamSpace[0];
188       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
189       if(!fTrackID)
190         {
191           cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
192           fTrackID = new TrackIndex*[GetNEtaSegments()];
193           for(Int_t i=0; i<GetNEtaSegments(); i++)
194             fTrackID[i] = new TrackIndex[ncells];
195         }
196     }
197 #endif
198   AliL3Histogram *hist = fParamSpace[0];
199   Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
200   if(!fRowCount)
201     {
202       cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fRowCount"<<endl;
203       fRowCount = new UChar_t*[GetNEtaSegments()];
204       for(Int_t i=0; i<GetNEtaSegments(); i++)
205         fRowCount[i] = new UChar_t[ncells];
206     }
207   if(!fGapCount)
208     {
209       cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fGapCount"<<endl;
210       fGapCount = new UChar_t*[GetNEtaSegments()];
211       for(Int_t i=0; i<GetNEtaSegments(); i++)
212         fGapCount[i] = new UChar_t[ncells];
213     }
214   if(!fCurrentRowCount)
215     {
216       cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fCurrentRowCount"<<endl;
217       fCurrentRowCount = new UChar_t*[GetNEtaSegments()];
218       for(Int_t i=0; i<GetNEtaSegments(); i++)
219         fCurrentRowCount[i] = new UChar_t[ncells];
220     }
221
222   //create lookup table for sin and cos
223   Int_t nbins = hist->GetNbinsY()+2;
224   fLUT2sinphi0up=new Float_t[nbins];
225   fLUT2cosphi0up=new Float_t[nbins];
226   fLUT2sinphi0low=new Float_t[nbins];
227   fLUT2cosphi0low=new Float_t[nbins];
228   Float_t hist_bin=hist->GetBinWidthY()/2.0;
229   for(Int_t i=hist->GetFirstYbin(); i<=hist->GetLastYbin(); i++){
230     Float_t phi0=hist->GetBinCenterY(i);
231     fLUT2sinphi0low[i]=2.*sin(phi0+hist_bin);
232     fLUT2cosphi0low[i]=2.*cos(phi0+hist_bin);
233     fLUT2sinphi0up[i]=2.*sin(phi0-hist_bin);
234     fLUT2cosphi0up[i]=2.*cos(phi0-hist_bin);
235   }  
236
237   //create lookup table for z of the digits
238   Int_t ntimebins = AliL3Transform::GetNTimeBins();
239   fLUTforwardZ = new Float_t[ntimebins];
240   fLUTforwardZ2 = new Float_t[ntimebins];
241   fLUTbackwardZ = new Float_t[ntimebins];
242   fLUTbackwardZ2 = new Float_t[ntimebins];
243   for(Int_t i=0; i<ntimebins; i++){
244     Float_t z;
245     z=AliL3Transform::GetZFast(0,i,GetZVertex());
246     fLUTforwardZ[i]=z;
247     fLUTforwardZ2[i]=z*z;
248     z=AliL3Transform::GetZFast(18,i,GetZVertex());
249     fLUTbackwardZ[i]=z;
250     fLUTbackwardZ2[i]=z*z;
251   }
252
253 }
254
255 void AliL3HoughTransformerRow::Reset()
256 {
257   //Reset all the histograms
258
259   if(!fParamSpace)
260     {
261       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
262         <<"No histograms to reset"<<ENDLOG;
263       return;
264     }
265   
266   for(Int_t i=0; i<GetNEtaSegments(); i++)
267     fParamSpace[i]->Reset();
268   
269 #ifdef do_mc
270   if(fDoMC)
271     {
272       AliL3Histogram *hist = fParamSpace[0];
273       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
274       for(Int_t i=0; i<GetNEtaSegments(); i++)
275         memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
276     }
277 #endif
278   AliL3Histogram *hist = fParamSpace[0];
279   Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
280   for(Int_t i=0; i<GetNEtaSegments(); i++)
281     {
282       memset(fRowCount[i],0,ncells*sizeof(UChar_t));
283       memset(fGapCount[i],1,ncells*sizeof(UChar_t));
284       memset(fCurrentRowCount[i],160,ncells*sizeof(UChar_t));
285     }
286 }
287
288 Int_t AliL3HoughTransformerRow::GetEtaIndex(Double_t eta)
289 {
290   //Return the histogram index of the corresponding eta. 
291
292   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
293   Double_t index = (eta-GetEtaMin())/etaslice;
294   return (Int_t)index;
295 }
296
297 inline AliL3Histogram *AliL3HoughTransformerRow::GetHistogram(Int_t eta_index)
298 {
299   if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
300     return 0;
301   if(!fParamSpace[eta_index])
302     return 0;
303   return fParamSpace[eta_index];
304 }
305
306 Double_t AliL3HoughTransformerRow::GetEta(Int_t eta_index,Int_t slice)
307 {
308   Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
309   Double_t eta=0;
310   eta=(Double_t)((eta_index+0.5)*eta_slice);
311   return eta;
312 }
313
314 struct EtaRow {
315   UChar_t start_pad;
316   UChar_t end_pad;
317   Bool_t found;
318   Float_t start_y;
319 #ifdef do_mc
320   Int_t mc_labels[MaxTrack];
321 #endif
322 };
323
324 void AliL3HoughTransformerRow::TransformCircle()
325 {
326   Int_t n_eta_segments = GetNEtaSegments();
327   Double_t eta_min = GetEtaMin();
328   Double_t eta_slice = (GetEtaMax() - eta_min)/n_eta_segments;
329
330   Int_t lower_threshold = GetLowerThreshold();
331
332   //Assumes that all the etaslice histos are the same!
333   AliL3Histogram *h = fParamSpace[0];
334   Int_t first_bin = h->GetFirstYbin();
335   Int_t last_bin = h->GetLastYbin();
336   Float_t x_min = h->GetXmin();
337   Float_t x_max = h->GetXmax();
338   Float_t x_bin = (x_max-x_min)/h->GetNbinsX();
339   Int_t first_binx = h->GetFirstXbin()-1;
340   Int_t last_binx = h->GetLastXbin()+1;
341   Int_t nbinx = h->GetNbinsX()+2;
342
343   UChar_t last_pad;
344   Int_t last_eta_index;
345   EtaRow *eta_clust = new EtaRow[n_eta_segments];
346
347   AliL3DigitRowData *tempPt = GetDataPointer();
348   if(!tempPt)
349     {
350       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
351         <<"No input data "<<ENDLOG;
352       return;
353     }
354
355   Int_t ipatch = GetPatch();
356   Double_t pad_pitch = AliL3Transform::GetPadPitchWidth(ipatch);
357   Int_t islice = GetSlice();
358   Float_t *fLUTz;
359   Float_t *fLUTz2;
360   if(islice < 18) {
361     fLUTz = fLUTforwardZ;
362     fLUTz2 = fLUTforwardZ2;
363   }
364   else {
365     fLUTz = fLUTbackwardZ;
366     fLUTz2 = fLUTbackwardZ2;
367   }
368
369   Int_t npads_in_patch = AliL3Transform::GetNPads(AliL3Transform::GetLastRow(ipatch))+2;
370   Bool_t **IsPrevRow = new Bool_t*[n_eta_segments];
371   Int_t nrows_in_patch = AliL3Transform::GetLastRow(ipatch)-AliL3Transform::GetFirstRow(ipatch)+2;
372   for(Int_t i=0; i<n_eta_segments; i++) {
373     IsPrevRow[i] = new Bool_t[npads_in_patch*nrows_in_patch];
374     memset(IsPrevRow[i],1,npads_in_patch*nrows_in_patch*sizeof(Bool_t));
375   }
376
377   //Loop over the padrows:
378   for(UChar_t i=AliL3Transform::GetFirstRow(ipatch); i<=AliL3Transform::GetLastRow(ipatch); i++)
379     {
380       Int_t npads = AliL3Transform::GetNPads((Int_t)i)-1;
381
382       last_pad = 255;
383       //Flush eta clusters array
384       memset(eta_clust,0,n_eta_segments*sizeof(EtaRow));  
385
386       Float_t x = AliL3Transform::Row2X((Int_t)i);
387       Float_t x2 = x*x;
388       Float_t y,r2;
389
390       Int_t lrow = (i-AliL3Transform::GetFirstRow(ipatch))*npads_in_patch;
391
392       //Get the data on this padrow:
393       AliL3DigitData *digPt = tempPt->fDigitData;
394       if((Int_t)i != (Int_t)tempPt->fRow)
395         {
396           cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<(Int_t)i<<" "<<(Int_t)tempPt->fRow<<endl;
397           continue;
398         }
399       //      cout<<" Starting row "<<i<<endl;
400       //Loop over the data on this padrow:
401       for(UInt_t j=0; j<tempPt->fNDigit; j++)
402         {
403           UShort_t charge = digPt[j].fCharge;
404           if((Int_t)charge <= lower_threshold)
405             continue;
406           UChar_t pad = digPt[j].fPad;
407           UShort_t time = digPt[j].fTime;
408           if(pad != last_pad)
409             {
410               y = (pad-0.5*npads)*pad_pitch;
411               Float_t y2 = y*y;
412               r2 = x2 + y2;
413               last_eta_index = -1;
414             }
415
416           //Transform data to local cartesian coordinates:
417           Float_t z = fLUTz[(Int_t)time];
418           Float_t z2 = fLUTz2[(Int_t)time];
419           //Calculate the eta:
420           Double_t r = sqrt(r2+z2);
421           Double_t eta = 0.5 * log((r+z)/(r-z));
422           //Get the corresponding index, which determines which histogram to fill:
423           Int_t eta_index = (Int_t)((eta-eta_min)/eta_slice);
424
425 #ifndef do_mc
426           if(eta_index == last_eta_index) continue;
427 #endif
428           //      cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<eta_index<<endl;
429           
430           if(eta_index < 0 || eta_index >= n_eta_segments)
431             continue;
432
433           if(!eta_clust[eta_index].found)
434             {
435               eta_clust[eta_index].start_pad = pad;
436               eta_clust[eta_index].end_pad = pad;
437               eta_clust[eta_index].start_y = y - pad_pitch/2.0;
438               eta_clust[eta_index].found = 1;
439 #ifdef do_mc
440               if(fDoMC)
441                 {
442                   for(Int_t t=0; t<3; t++)
443                     {
444                       Int_t label = digPt[j].fTrackID[t];
445                       if(label < 0) break;
446                       UInt_t c;
447                       for(c=0; c<(MaxTrack-1); c++)
448                         if(eta_clust[eta_index].mc_labels[c] == label || eta_clust[eta_index].mc_labels[c] == 0)
449                           break;
450                       //                      if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
451                       eta_clust[eta_index].mc_labels[c] = label;
452                     }
453                 }
454 #endif
455               continue;
456             }
457           else
458             {
459               if(pad <= (eta_clust[eta_index].end_pad + 1))
460                 {
461                   eta_clust[eta_index].end_pad = pad;
462 #ifdef do_mc
463                   if(fDoMC)
464                     {
465                       for(Int_t t=0; t<3; t++)
466                         {
467                           Int_t label = digPt[j].fTrackID[t];
468                           if(label < 0) break;
469                           UInt_t c;
470                           for(c=0; c<(MaxTrack-1); c++)
471                             if(eta_clust[eta_index].mc_labels[c] == label || eta_clust[eta_index].mc_labels[c] == 0)
472                               break;
473                           //                      if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
474                           eta_clust[eta_index].mc_labels[c] = label;
475                         }
476                     }
477 #endif
478                 }
479               else
480                 {
481
482                   Bool_t fill_cluster = kFALSE;
483                   Bool_t *IsPrevRow2 = &IsPrevRow[eta_index][lrow+eta_clust[eta_index].start_pad];
484                   for(Int_t ipad = 0; ipad <= (eta_clust[eta_index].end_pad - eta_clust[eta_index].start_pad + 2); ipad++)
485                     {
486                       if(*IsPrevRow2)
487                         {
488                           fill_cluster = kTRUE;
489                           break;
490                         }
491                       IsPrevRow2++;
492                     }
493
494                   //              if(eta_index == 51)
495                   //                cout<<" Cluster to fill:"<<" row "<<(Int_t)i<<" pads from "<<(Int_t)eta_clust[eta_index].start_pad<<"("<<eta_clust[eta_index].start_y<<") to "<<(Int_t)eta_clust[eta_index].end_pad<<"("<<(eta_clust[eta_index].end_pad-0.5*(npads-1))*pad_pitch<<") fill "<<fill_cluster<<endl;
496
497                   if(fill_cluster) {
498                   //              cout<<" Cluster found at etaslice "<<eta_index<<" from pad "<<(Int_t)eta_clust[eta_index].start_pad<<" to pad "<<(Int_t)eta_clust[eta_index].end_pad<<endl;
499                   //              Bool_t fill_next_row = kFALSE;
500
501                   UChar_t *nrows = fRowCount[eta_index];
502                   UChar_t *ngaps = fGapCount[eta_index];
503                   UChar_t *currentrow = fCurrentRowCount[eta_index];
504
505                   //Do the transformation:
506                   Float_t start_y = eta_clust[eta_index].start_y;
507                   if(eta_clust[eta_index].start_pad == 0)
508                     start_y -= 3.0*pad_pitch;
509                   Float_t R1 = x2 + start_y*start_y;
510                   Float_t x_over_R1 = x/R1;
511                   Float_t start_y_over_R1 = start_y/R1;
512                   Float_t end_y = (eta_clust[eta_index].end_pad-0.5*(npads-1))*pad_pitch;
513                   if(eta_clust[eta_index].end_pad == npads)
514                     end_y += 3.0*pad_pitch;
515                   Float_t R2 = x2 + end_y*end_y; 
516                   Float_t x_over_R2 = x/R2;
517                   Float_t end_y_over_R2 = end_y/R2;
518
519                   //Fill the histogram along the phirange
520                   for(Int_t b=first_bin; b<=last_bin; b++)
521                     {
522                       Float_t kappa1 = start_y_over_R1*fLUT2cosphi0low[b]-x_over_R1*fLUT2sinphi0low[b];
523                       if(kappa1>x_max) continue;
524                       Float_t kappa2 = end_y_over_R2*fLUT2cosphi0up[b]-x_over_R2*fLUT2sinphi0up[b];
525                       if(kappa2<x_min) break;
526
527                       Int_t binx1 = 1 + (Int_t)((kappa1-x_min)/x_bin);
528                       if(binx1<first_binx) binx1 = first_binx;
529                       Int_t binx2 = 1 + (Int_t)((kappa2-x_min)/x_bin);
530                       if(binx2>last_binx) binx2 = last_binx;
531                       //                      if(eta_index == 51)
532                       //                        cout<<"     phi bin "<<b<<" kappa1 "<<kappa1<<" kappa2 "<<kappa2<<" from bin "<<binx1<<" to "<<binx2<<endl;
533                       Int_t temp_bin = b*nbinx + binx1;
534                       UChar_t *nrows2 = nrows + temp_bin;
535                       UChar_t *ngaps2 = ngaps + temp_bin;
536                       UChar_t *currentrow2 = currentrow + temp_bin;
537                       for(Int_t bin=binx1;bin<=binx2;bin++)
538                         {
539                           //Assumes threshold about 32
540                           if(*ngaps2 < 3) {
541                             if(i != *currentrow2)
542                               {
543                                 (*nrows2)++;
544                                 if(i > (*currentrow2 + 1))
545                                   (*ngaps2)++;
546                                 *currentrow2=i;
547                               }
548 #ifdef do_mc
549                             if(fDoMC)
550                               {
551                                 for(UInt_t t=0;t<(MaxTrack-1); t++)
552                                   {
553                                     Int_t label = eta_clust[eta_index].mc_labels[t];
554                                     if(label == 0) break;
555                                     UInt_t c;
556                                     Int_t temp_bin2 = b*nbinx + bin;
557                                     for(c=0; c<(MaxTrack-1); c++)
558                                       if(fTrackID[eta_index][temp_bin2].fLabel[c] == label || fTrackID[eta_index][temp_bin2].fNHits[c] == 0)
559                                         break;
560                                     //                              if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
561                                     fTrackID[eta_index][temp_bin2].fLabel[c] = label;
562                                     if(fTrackID[eta_index][temp_bin2].fCurrentRow[c] != i) {
563                                       fTrackID[eta_index][temp_bin2].fNHits[c]++;
564                                       fTrackID[eta_index][temp_bin2].fCurrentRow[c] = i;
565                                     }
566                                   }
567                               }
568 #endif
569                           }
570                           nrows2++;
571                           ngaps2++;
572                           currentrow2++;
573                         }
574                     }
575                   //End of the transformation
576
577                   /*              if(!fill_next_row) {
578                     Bool_t *IsCurrentRow = &IsPrevRow[eta_index][lrow+npads_in_patch];
579                     memset(&IsCurrentRow[eta_clust[eta_index].start_pad + 1],0,eta_clust[eta_index].end_pad - eta_clust[eta_index].start_pad + 1);
580                     } */
581                   }
582                   else {
583                     Bool_t *IsCurrentRow = &IsPrevRow[eta_index][lrow+npads_in_patch+eta_clust[eta_index].start_pad+1];
584                     memset(IsCurrentRow,0,eta_clust[eta_index].end_pad - eta_clust[eta_index].start_pad + 1);
585                   }
586
587                   Bool_t *IsCurrentRow = &IsPrevRow[eta_index][lrow+npads_in_patch+eta_clust[eta_index].end_pad+2];
588                   memset(IsCurrentRow,0,pad - eta_clust[eta_index].end_pad - 1);
589
590                   eta_clust[eta_index].start_pad = pad;
591                   eta_clust[eta_index].end_pad = pad;
592                   eta_clust[eta_index].start_y = y - pad_pitch/2.0;
593
594 #ifdef do_mc
595                   if(fDoMC)
596                     {
597                       memset(eta_clust[eta_index].mc_labels,0,MaxTrack);
598                       for(Int_t t=0; t<3; t++)
599                         {
600                           Int_t label = digPt[j].fTrackID[t];
601                           if(label < 0) break;
602                           UInt_t c;
603                           for(c=0; c<(MaxTrack-1); c++)
604                             if(eta_clust[eta_index].mc_labels[c] == label || eta_clust[eta_index].mc_labels[c] == 0)
605                               break;
606                           //                      if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
607                           eta_clust[eta_index].mc_labels[c] = label;
608                         }
609                     }
610 #endif
611                 }
612             }
613           last_pad = pad;
614           last_eta_index = eta_index;
615         }
616       //Write remaining clusters
617       for(Int_t eta_index = 0;eta_index < n_eta_segments;eta_index++)
618         {
619           //Check for empty row
620           if((eta_clust[eta_index].start_pad == 0) && (eta_clust[eta_index].end_pad == 0)) continue;
621
622           UChar_t *nrows = fRowCount[eta_index];
623           UChar_t *ngaps = fGapCount[eta_index];
624           UChar_t *currentrow = fCurrentRowCount[eta_index];
625
626           //Do the transformation:
627           Float_t start_y = eta_clust[eta_index].start_y;
628           if(eta_clust[eta_index].start_pad == 0)
629             start_y -= 3.0*pad_pitch;
630           Float_t R1 = x2 + start_y*start_y; 
631           Float_t x_over_R1 = x/R1;
632           Float_t start_y_over_R1 = start_y/R1;
633           Float_t end_y = (eta_clust[eta_index].end_pad-0.5*(npads-1))*pad_pitch;
634           if(eta_clust[eta_index].end_pad == npads)
635             end_y += 3.0*pad_pitch;
636           Float_t R2 = x2 + end_y*end_y; 
637           Float_t x_over_R2 = x/R2;
638           Float_t end_y_over_R2 = end_y/R2;
639
640           //Fill the histogram along the phirange
641           for(Int_t b=first_bin; b<=last_bin; b++)
642             {
643               Float_t kappa1 = start_y_over_R1*fLUT2cosphi0low[b]-x_over_R1*fLUT2sinphi0low[b];
644               if(kappa1>x_max) continue;
645               Float_t kappa2 = end_y_over_R2*fLUT2cosphi0up[b]-x_over_R2*fLUT2sinphi0up[b];
646               if(kappa2<x_min) break;
647
648               Int_t binx1 = 1 + (Int_t)((kappa1-x_min)/x_bin);
649               if(binx1<first_binx) binx1 = first_binx;
650               Int_t binx2 = 1 + (Int_t)((kappa2-x_min)/x_bin);
651               if(binx2>last_binx) binx2 = last_binx;
652
653               Int_t temp_bin = b*nbinx + binx1;
654               UChar_t *nrows2 = nrows + temp_bin;
655               UChar_t *ngaps2 = ngaps + temp_bin;
656               UChar_t *currentrow2 = currentrow + temp_bin;
657               for(Int_t bin=binx1;bin<=binx2;bin++)
658                 {
659                   //Assumes threshold about 32
660                   if(*ngaps2 < 3) {
661                     if(i != *currentrow2)
662                       {
663                         (*nrows2)++;
664                         if(i > (*currentrow2 + 1))
665                           (*ngaps2)++;
666                         *currentrow2=i;
667                       }
668 #ifdef do_mc
669                     if(fDoMC)
670                       {
671                         for(UInt_t t=0;t<(MaxTrack-1); t++)
672                           {
673                             Int_t label = eta_clust[eta_index].mc_labels[t];
674                             if(label == 0) break;
675                             UInt_t c;
676                             Int_t temp_bin2 = b*nbinx + bin;
677                             for(c=0; c<(MaxTrack-1); c++)
678                               if(fTrackID[eta_index][temp_bin2].fLabel[c] == label || fTrackID[eta_index][temp_bin2].fNHits[c] == 0)
679                                 break;
680                             //                      if(c == MaxTrack-1) cout<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
681                             fTrackID[eta_index][temp_bin2].fLabel[c] = label;
682                             if(fTrackID[eta_index][temp_bin2].fCurrentRow[c] != i) {
683                               fTrackID[eta_index][temp_bin2].fNHits[c]++;
684                               fTrackID[eta_index][temp_bin2].fCurrentRow[c] = i;
685                             }
686                           }
687                       }
688 #endif
689                   }
690                   nrows2++;
691                   ngaps2++;
692                   currentrow2++;
693                 }
694             }
695           //End of the transformation
696
697           if(eta_clust[eta_index].end_pad < npads)
698             {
699               Bool_t *IsCurrentRow = &IsPrevRow[eta_index][lrow+npads_in_patch+eta_clust[eta_index].end_pad+2];
700               memset(IsCurrentRow,0,npads - eta_clust[eta_index].end_pad);
701             }
702
703         }
704
705       //Move the data pointer to the next padrow:
706       AliL3MemHandler::UpdateRowPointer(tempPt);
707     }
708
709
710   for(Int_t i=0; i<n_eta_segments; i++)
711     {
712       delete [] IsPrevRow[i];
713     }
714   delete [] IsPrevRow;
715
716   delete [] eta_clust;
717 }
718
719 Int_t AliL3HoughTransformerRow::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
720 {
721   if(!fDoMC)
722     {
723       cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
724       return -1;
725     }
726   
727 #ifdef do_mc
728   if(eta_index < 0 || eta_index > GetNEtaSegments())
729     {
730       cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl;
731       return -1;
732     }
733   AliL3Histogram *hist = fParamSpace[eta_index];
734   Int_t bin = hist->FindBin(kappa,psi);
735   Int_t label=-1;
736   Int_t max=0;
737   for(UInt_t i=0; i<(MaxTrack-1); i++)
738     {
739       Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
740       if(nhits == 0) break;
741       if(nhits > max)
742         {
743           max = nhits;
744           label = fTrackID[eta_index][bin].fLabel[i];
745         }
746     }
747   Int_t label2=-1;
748   Int_t max2=0;
749   for(UInt_t i=0; i<(MaxTrack-1); i++)
750     {
751       Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
752       if(nhits == 0) break;
753       if(nhits > max2)
754         {
755           if(fTrackID[eta_index][bin].fLabel[i]!=label) {
756             max2 = nhits;
757             label2 = fTrackID[eta_index][bin].fLabel[i];
758           }
759         }
760     }
761   //  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;
762   return label;
763 #endif
764   cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
765   return -1;
766 }
767
768 UChar_t *AliL3HoughTransformerRow::GetRowCount(Int_t eta_index)
769 {
770   return fRowCount[eta_index];
771 }
772
773 UChar_t *AliL3HoughTransformerRow::GetGapCount(Int_t eta_index)
774 {
775   return fGapCount[eta_index];
776 }