4955c20fb6da81cf0bd3f0c279fd75d5bd05c44c
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformer.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 "AliL3HoughTransformer.h"
10 #include "AliL3MemHandler.h"
11 #include "AliL3Transform.h"
12 #include "AliL3DigitData.h"
13 #include "AliL3HistogramAdaptive.h"
14
15 #if __GNUC__ == 3
16 using namespace std;
17 #endif
18
19 /** \class AliL3HoughTransformer
20 <pre>
21 //_____________________________________________________________
22 // AliL3HoughTransformer
23 //
24 // Hough transformation class
25 //
26 </pre>
27 */
28
29 ClassImp(AliL3HoughTransformer)
30
31 AliL3HoughTransformer::AliL3HoughTransformer()
32 {
33   //Default constructor
34   fParamSpace = 0;
35   fDoMC = kFALSE;;
36   fEtaOverlap=kFALSE;
37 #ifdef do_mc
38   fTrackID = 0;
39 #endif
40 }
41
42 AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t netasegments,Bool_t DoEtaOverlap,Bool_t /*DoMC*/) : AliL3HoughBaseTransformer(slice,patch,netasegments)
43 {
44   //Normal constructor
45   fParamSpace = 0;
46   fDoMC = kFALSE;
47   fEtaOverlap = DoEtaOverlap;
48   fDoMC=kFALSE;
49 #ifdef do_mc
50   fTrackID = 0;
51 #endif
52 }
53
54 AliL3HoughTransformer::~AliL3HoughTransformer()
55 {
56   // Dtor
57   DeleteHistograms();
58 #ifdef do_mc
59   if(fTrackID)
60     {
61       for(Int_t i=0; i<GetNEtaSegments(); i++)
62         {
63           if(!fTrackID[i]) continue;
64           delete fTrackID[i];
65         }
66       delete [] fTrackID;
67     }
68 #endif
69 }
70
71 void AliL3HoughTransformer::DeleteHistograms()
72 {
73   // Clean up
74   if(!fParamSpace)
75     return;
76   for(Int_t i=0; i<GetNEtaSegments(); i++)
77     {
78       if(!fParamSpace[i]) continue;
79       delete fParamSpace[i];
80     }
81   delete [] fParamSpace;
82 }
83
84 void AliL3HoughTransformer::CreateHistograms(Float_t ptmin,Float_t ptmax,Float_t ptres,
85                                              Int_t nybin,Float_t psi)
86 {
87   //Create histograms.
88   //_Only_ to be used in case of the adaptive histograms!
89   //phimax is given in radians!!
90   
91   if(ptmin > ptmax)
92     {
93       cerr<<"AliL3HoughTransformer::CreateHistograms: Error in ptrange "<<ptmin<<" "<<ptmax<<endl;
94       return;
95     }
96   if(psi < 0)
97     {
98       cerr<<"AliL3HoughTransformer::CreateHistograms: Wrong psi-angle "<<psi<<endl;
99       return;
100     }
101   
102   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
103   Char_t histname[256];
104   Int_t i;
105   for(i=0; i<GetNEtaSegments(); i++)
106     {
107       sprintf(histname,"paramspace_%d",i);
108       fParamSpace[i] = new AliL3HistogramAdaptive(histname,ptmin,ptmax,ptres,nybin,-psi,psi);
109     }
110 }
111
112 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Float_t ptmin,
113                                              Int_t nybin,Float_t phimin,Float_t phimax)
114 {
115   //Create the histograms (parameter space).
116   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
117   //The arguments give the range and binning; 
118   //nxbin = #bins in kappa
119   //nybin = #bins in phi0
120   //ptmin = mimium Pt of track (corresponding to maximum kappa)
121   //phimin = mimimum phi0 
122   //phimax = maximum phi0 
123     
124   Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/ptmin;
125   //Double_t torad = AliL3Transform::Pi()/180;
126   
127   CreateHistograms(nxbin,-1.*x,x,nybin,phimin/**torad*/,phimax/**torad*/);
128 }
129
130 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
131                                              Int_t nybin,Float_t ymin,Float_t ymax)
132 {
133   //Create the histograms (parameter space).
134   //nxbin = #bins in X
135   //nybin = #bins in Y
136   //xmin xmax ymin ymax = histogram limits in X and Y
137   
138   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
139   
140   Char_t histname[256];
141   for(Int_t i=0; i<GetNEtaSegments(); i++)
142     {
143       sprintf(histname,"paramspace_%d",i);
144       //fParamSpace[i] = new AliL3HistogramAdaptive(histname,0.5,1.5,0.05,nybin,ymin,ymax);
145       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
146     }
147   
148 #ifdef do_mc
149   if(fDoMC)
150     {
151       AliL3Histogram *hist = fParamSpace[0];
152       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
153       cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(AliL3TrackIndex)<<" bytes to fTrackID"<<endl;
154       fTrackID = new AliL3TrackIndex*[GetNEtaSegments()];
155       for(Int_t i=0; i<GetNEtaSegments(); i++)
156         fTrackID[i] = new AliL3TrackIndex[ncells];
157     }
158 #endif
159 }
160
161 void AliL3HoughTransformer::Reset()
162 {
163   //Reset all the histograms
164
165   if(!fParamSpace)
166     {
167       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
168         <<"No histograms to reset"<<ENDLOG;
169       return;
170     }
171   
172   for(Int_t i=0; i<GetNEtaSegments(); i++)
173     fParamSpace[i]->Reset();
174   
175 #ifdef do_mc
176   if(fDoMC)
177     {
178       AliL3Histogram *hist = fParamSpace[0];
179       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
180       for(Int_t i=0; i<GetNEtaSegments(); i++)
181         memset(fTrackID[i],0,ncells*sizeof(AliL3TrackIndex));
182     }
183 #endif
184 }
185
186 Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta) const
187 {
188   //Return the histogram index of the corresponding eta. 
189
190   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
191   Double_t index = (eta-GetEtaMin())/etaslice;
192   return (Int_t)index;
193 }
194
195 void AliL3HoughTransformer::GetEtaIndexes(Double_t eta,Int_t *indexes) const
196 {
197   //Return histogram indexes in case of overlapping etaslices.
198   
199   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
200   Int_t index = (Int_t)((eta-GetEtaMin())/etaslice);
201   if(index%2 == 0)
202     {
203       indexes[0] = index;
204       indexes[1] = index - 1;
205     }
206   else
207     {
208       indexes[0] = index - 1;
209       indexes[1] = index;
210     }
211 }
212
213 AliL3Histogram *AliL3HoughTransformer::GetHistogram(Int_t etaindex)
214 {
215   // Return a pointer to the histogram which contains etaindex eta slice
216   if(!fParamSpace || etaindex >= GetNEtaSegments() || etaindex < 0)
217     return 0;
218   if(!fParamSpace[etaindex])
219     return 0;
220   return fParamSpace[etaindex];
221 }
222
223 Double_t AliL3HoughTransformer::GetEta(Int_t etaindex,Int_t /*slice*/) const
224 {
225   // Return eta calculated in the middle of the eta slice
226   Double_t etaslice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
227   Double_t eta=0;
228   if(fEtaOverlap)
229     {
230       Int_t index = etaindex + 1;
231       eta=(Double_t)((index)*etaslice);
232     }
233   else
234     eta=(Double_t)((etaindex+0.5)*etaslice);
235   return eta;
236 }
237
238 void AliL3HoughTransformer::TransformCircle()
239 {
240   //Transform the input data with a circle HT.
241   //The function loops over all the data, and transforms each pixel with the equations:
242   // 
243   //kappa = 2/r*sin(phi - phi0)
244   //
245   //where r = sqrt(x*x +y*y), and phi = arctan(y/x)
246   //
247   //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
248   //which histogram in which the pixel should be transformed, the eta-value is calculated
249   //and the proper histogram index is found by GetEtaIndex(eta).
250
251
252   AliL3DigitRowData *tempPt = GetDataPointer();
253   if(!tempPt)
254     {
255       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
256         <<"No input data "<<ENDLOG;
257       return;
258     }
259   
260   //Loop over the padrows:
261   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
262     {
263       //Get the data on this padrow:
264       AliL3DigitData *digPt = tempPt->fDigitData;
265       if(i != (Int_t)tempPt->fRow)
266         {
267           cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl;
268           continue;
269         }
270       
271       //Loop over the data on this padrow:
272       for(UInt_t j=0; j<tempPt->fNDigit; j++)
273         {
274           UShort_t charge = digPt[j].fCharge;
275           UChar_t pad = digPt[j].fPad;
276           UShort_t time = digPt[j].fTime;
277           if((Int_t)charge <= GetLowerThreshold())
278             continue;
279           
280           if((Int_t)charge > GetUpperThreshold())
281             charge = GetUpperThreshold();
282           
283           Int_t sector,row;
284           Float_t xyz[3];
285
286           //Transform data to local cartesian coordinates:
287           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
288           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
289                   
290           //Calculate the eta:
291           Double_t eta = AliL3Transform::GetEta(xyz);
292           
293           //Get the corresponding index, which determines which histogram to fill:
294           Int_t etaindex = GetEtaIndex(eta);
295                   
296           if(etaindex < 0 || etaindex >= GetNEtaSegments())
297             continue;
298           
299           //Get the correct histogrampointer:
300           AliL3Histogram *hist = fParamSpace[etaindex];
301           if(!hist)
302             {
303               cerr<<"AliL3HoughTransformer::TransformCircle : Error getting histogram in index "<<etaindex<<endl;
304               continue;
305             }
306           
307           //Do the transformation:
308           Float_t r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); 
309           Float_t phi = AliL3Transform::GetPhi(xyz);
310           
311           
312           //Fill the histogram along the phirange
313           for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
314             {
315               Float_t phi0 = hist->GetBinCenterY(b);
316               Float_t kappa = 2*sin(phi - phi0)/r;
317               //hist->Fill(kappa,phi0,(int)rint(log((Float_t)charge)));
318               hist->Fill(kappa,phi0,charge);
319               //hist->Fill(kappa,phi0,1);
320 #ifdef do_mc
321               if(fDoMC)
322                 {
323                   Int_t bin = hist->FindBin(kappa,phi0);
324                   for(Int_t t=0; t<3; t++)
325                     {
326                       Int_t label = digPt[j].fTrackID[t];
327                       if(label < 0) break;
328                       UInt_t c;
329                       for(c=0; c<MaxTrack; c++)
330                         if(fTrackID[etaindex][bin].fLabel[c] == label || fTrackID[etaindex][bin].fNHits[c] == 0)
331                           break;
332                       if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
333                       fTrackID[etaindex][bin].fLabel[c] = label;
334                       fTrackID[etaindex][bin].fNHits[c]++;
335                     }
336                 }
337 #endif
338             }
339         }
340       
341       //Move the data pointer to the next padrow:
342       AliL3MemHandler::UpdateRowPointer(tempPt);
343     }
344 }
345
346 struct AliL3Digit {
347   Int_t fRow; // Digit padrow
348   Double_t fR; // Digit radius in local coordinate system
349   Double_t fPhi; // Digit Phi angle in local coordinate system
350   Int_t fCharge; // Digit charge
351   AliL3Digit *fNext; // Next digit
352 };
353
354 struct AliL3EtaContainer {
355   AliL3Digit *fFirst; //First digit
356   AliL3Digit *fLast; //Last digit
357 };
358
359 void AliL3HoughTransformer::TransformCircleC(Int_t *rowrange,Int_t every)
360 {
361   //Circle transform, using combinations of every 2 points lying
362   //on different padrows and within the same etaslice.
363   
364   AliL3DigitRowData *tempPt = GetDataPointer();
365   if(!tempPt)
366     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
367       <<"No input data "<<ENDLOG;
368   
369   Int_t minrow = AliL3Transform::GetFirstRow(GetPatch());
370   Int_t maxrow = AliL3Transform::GetLastRow(GetPatch());
371   if(rowrange)
372     {
373       minrow = rowrange[0];
374       maxrow = rowrange[1];
375       if(minrow < AliL3Transform::GetFirstRow(GetPatch()) || minrow >= AliL3Transform::GetLastRow(GetPatch()))
376         minrow = AliL3Transform::GetFirstRow(GetPatch());
377       if(maxrow < AliL3Transform::GetFirstRow(GetPatch()) || maxrow >= AliL3Transform::GetLastRow(GetPatch()))
378         maxrow = AliL3Transform::GetLastRow(GetPatch());
379       if(minrow > maxrow || maxrow==minrow)
380         {
381           cerr<<"AliL3HoughTransformer::TransformCircleC : Bad row range "<<minrow<<" "<<maxrow<<endl;
382           return;
383         }
384     }
385   else
386     {
387       minrow = AliL3Transform::GetFirstRow(GetPatch());
388       maxrow = AliL3Transform::GetLastRow(GetPatch());
389     }
390       
391   Int_t counter=0;
392   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
393     {
394       counter += tempPt->fNDigit;
395       AliL3MemHandler::UpdateRowPointer(tempPt);
396     }
397   
398   Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1);
399   AliL3EtaContainer *etaPt = new AliL3EtaContainer[bound];
400   memset(etaPt,0,bound*sizeof(AliL3EtaContainer));  
401   
402   AliL3Digit *digits = new AliL3Digit[counter];
403   cout<<"Allocating "<<counter*sizeof(AliL3Digit)<<" bytes to digitsarray"<<endl;
404   memset(digits,0,counter*sizeof(AliL3Digit));
405
406   Int_t sector,row,totcharge,pad,time,charge;
407   Double_t r1,r2,phi1,phi2,eta,kappa,phi0;
408   Float_t xyz[3];
409   
410   counter=0;
411   tempPt = GetDataPointer();
412   
413   cout<<"Calculating digits in patch "<<GetPatch()<<endl;
414   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
415     {
416       AliL3DigitData *digPt = tempPt->fDigitData;
417       for(UInt_t di=0; di<tempPt->fNDigit; di++)
418         {
419           charge = digPt[di].fCharge;
420           pad = digPt[di].fPad;
421           time = digPt[di].fTime;
422           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
423           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
424           eta = AliL3Transform::GetEta(xyz);
425           
426           digits[counter].fRow = i;
427           digits[counter].fR = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
428           digits[counter].fPhi = atan2(xyz[1],xyz[0]);
429           digits[counter].fCharge = charge;
430
431           if(!fEtaOverlap)
432             {
433               Int_t etaindex = GetEtaIndex(eta);
434               
435               Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex;
436               
437               if(index > 0 && index < bound) 
438                 {
439                   if(etaPt[index].fFirst == 0)
440                     etaPt[index].fFirst = &digits[counter];
441                   else
442                     (etaPt[index].fLast)->fNext = &digits[counter];
443                   etaPt[index].fLast = &digits[counter];
444                 }
445             }
446           else
447             {
448               Int_t etaindex[2];
449               GetEtaIndexes(eta,etaindex);
450               Int_t index[2];
451               index[0] = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex[0];
452               index[1] = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex[1];
453               if(index[0] == index[1])
454                 {
455                   cerr<<"Same etaindexes "<<index[0]<<" "<<index[1]<<endl;
456                   exit(5);
457                 }
458               
459               Int_t ind = index[0];
460               if(ind > 0 && ind < bound)
461                 {
462                   if(etaPt[ind].fFirst == 0)
463                     etaPt[ind].fFirst = &digits[counter];
464                   else
465                     (etaPt[ind].fLast)->fNext = &digits[counter];
466                   etaPt[ind].fLast = &digits[counter];
467                 }
468               
469               ind = index[1];
470               if(ind > 0 && ind < bound)
471                 {
472                   if(etaPt[ind].fFirst == 0)
473                     etaPt[ind].fFirst = &digits[counter];
474                   else
475                     (etaPt[ind].fLast)->fNext = &digits[counter];
476                   etaPt[ind].fLast = &digits[counter];
477                 }
478             }
479
480           counter++;
481         }
482       AliL3MemHandler::UpdateRowPointer(tempPt);
483     }
484   
485   cout<<"Doing the combinatorics"<<endl;
486   
487   AliL3Digit *dPt1,*dPt2;
488   
489   for(Int_t e=0; e<GetNEtaSegments(); e++)
490     {
491       for(Int_t i=minrow; i<=maxrow; i+=every)
492         {
493           Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e;
494           
495           for(dPt1 = (AliL3Digit*)etaPt[index1].fFirst; dPt1 != 0; dPt1 = (AliL3Digit*)dPt1->fNext)
496             {
497               for(Int_t j=i+every; j<=maxrow; j+=every)
498                 {
499                   Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e;
500                   
501                   for(dPt2 = (AliL3Digit*)etaPt[index2].fFirst; dPt2 != 0; dPt2 = (AliL3Digit*)dPt2->fNext)
502                     {
503                       if(dPt1->fRow == dPt2->fRow)
504                         {
505                           cerr<<"same row; indexes "<<index1<<" "<<index2<<endl;
506                           exit(5);
507                         }
508                       
509                       //Get the correct histogrampointer:
510                       AliL3Histogram *hist = fParamSpace[e];
511                       if(!hist)
512                         {
513                           printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",i);
514                           continue;
515                         }
516                       
517                       //Do the transform:
518                       r1 = dPt1->fR;
519                       phi1 = dPt1->fPhi;
520                       r2 = dPt2->fR;
521                       phi2 = dPt2->fPhi;
522                       phi0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
523                       kappa = 2*sin(phi2-phi0)/r2;
524                       totcharge = dPt1->fCharge + dPt2->fCharge;
525                       hist->Fill(kappa,phi0,totcharge);
526                       
527                     }
528                 }
529             }
530         }
531     }
532
533   cout<<"done"<<endl;
534   delete [] etaPt;
535   delete [] digits;
536
537 }
538
539 void AliL3HoughTransformer::TransformLine(Int_t *rowrange,Float_t *phirange)
540 {
541   //Do a line transform on the data.
542
543
544   AliL3DigitRowData *tempPt = GetDataPointer();
545   if(!tempPt)
546     {
547       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data")
548         <<"No input data "<<ENDLOG;
549       return;
550     }
551   
552   Int_t minrow = AliL3Transform::GetFirstRow(GetPatch());
553   Int_t maxrow = AliL3Transform::GetLastRow(GetPatch());
554   if(rowrange)
555     {
556       minrow = rowrange[0];
557       maxrow = rowrange[1];
558       if(minrow < AliL3Transform::GetFirstRow(GetPatch()) || minrow >= AliL3Transform::GetLastRow(GetPatch()))
559         minrow = AliL3Transform::GetFirstRow(GetPatch());
560       if(maxrow < AliL3Transform::GetFirstRow(GetPatch()) || maxrow >= AliL3Transform::GetLastRow(GetPatch()))
561         maxrow = AliL3Transform::GetLastRow(GetPatch());
562       if(minrow > maxrow || maxrow==minrow)
563         {
564           cerr<<"AliL3HoughTransformer::TransformCircleC : Bad row range "<<minrow<<" "<<maxrow<<endl;
565           return;
566         }
567     }
568
569   for(Int_t i=minrow; i<=maxrow; i++)
570     {
571       AliL3DigitData *digPt = tempPt->fDigitData;
572       if(i != (Int_t)tempPt->fRow)
573         {
574           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
575           continue;
576         }
577       for(UInt_t j=0; j<tempPt->fNDigit; j++)
578         {
579           UShort_t charge = digPt[j].fCharge;
580           UChar_t pad = digPt[j].fPad;
581           UShort_t time = digPt[j].fTime;
582           if(charge < GetLowerThreshold())
583             continue;
584           Int_t sector,row;
585           Float_t xyz[3];
586           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
587           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
588           
589           if(phirange)
590             {
591               Float_t phi = AliL3Transform::GetPhi(xyz);
592               if(phi < phirange[0] || phi > phirange[1])
593                 continue;
594             }
595           Float_t eta = AliL3Transform::GetEta(xyz);
596           Int_t etaindex = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
597           if(etaindex < 0 || etaindex >= GetNEtaSegments())
598             continue;
599           
600           xyz[0] = xyz[0] - AliL3Transform::Row2X(minrow);
601           
602           //Get the correct histogram:
603           AliL3Histogram *hist = fParamSpace[etaindex];
604           if(!hist)
605             {
606               printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",etaindex);
607               continue;
608             }
609           for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
610             {
611               Double_t theta = hist->GetBinCenterX(xbin);
612               Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
613               hist->Fill(theta,rho,charge);
614             }
615         }
616       AliL3MemHandler::UpdateRowPointer(tempPt);
617     }
618   
619 }
620
621 struct AliL3LDigit {
622   Int_t fRow; // Digit rowpad
623   Int_t fCharge; // Digit charge
624   Float_t fY; // Y position of the digit in the local coor system
625   AliL3LDigit *fNext; // Next digit
626 };
627 struct AliL3LEtaContainer {
628   AliL3LDigit *fFirst; //First digit
629   AliL3LDigit *fLast; //Last digit
630 };
631 void AliL3HoughTransformer::TransformLineC(Int_t *rowrange,Float_t *phirange)
632 {
633   //Circle transform ??
634   AliL3DigitRowData *tempPt = GetDataPointer();
635   if(!tempPt)
636     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
637       <<"No input data "<<ENDLOG;
638
639       
640   Int_t counter=0;
641   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
642     {
643       counter += tempPt->fNDigit;
644       AliL3MemHandler::UpdateRowPointer(tempPt);
645     }
646   
647   Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1);
648   AliL3LEtaContainer *etaPt = new AliL3LEtaContainer[bound];
649   memset(etaPt,0,bound*sizeof(AliL3LEtaContainer));  
650   
651   AliL3LDigit *digits = new AliL3LDigit[counter];
652   cout<<"Allocating "<<counter*sizeof(AliL3LDigit)<<" bytes to digitsarray"<<endl;
653   memset(digits,0,counter*sizeof(AliL3LDigit));
654
655   Int_t sector,row;
656   Float_t xyz[3];
657   
658   counter=0;
659   tempPt = GetDataPointer();
660   
661   cout<<"Calculating digits in patch "<<GetPatch()<<endl;
662   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
663     {
664       AliL3DigitData *digPt = tempPt->fDigitData;
665       for(UInt_t di=0; di<tempPt->fNDigit; di++)
666         {
667           Int_t charge = digPt[di].fCharge;
668           Int_t pad = digPt[di].fPad;
669           Int_t time = digPt[di].fTime;
670           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
671           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
672           Double_t eta = AliL3Transform::GetEta(xyz);
673           
674           Float_t phi = atan2(xyz[1],xyz[0]);
675           if(phi < phirange[0] || phi > phirange[1]) continue;
676           
677           digits[counter].fRow = i;
678           digits[counter].fY = xyz[1];
679           digits[counter].fCharge = charge;
680           
681           Int_t etaindex = GetEtaIndex(eta);
682           Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + etaindex;
683           
684           if(index > 0 && index < bound) 
685             {
686               if(etaPt[index].fFirst == 0)
687                 etaPt[index].fFirst = &digits[counter];
688               else
689                 (etaPt[index].fLast)->fNext = &digits[counter];
690               etaPt[index].fLast = &digits[counter];
691             }
692           counter++;
693         }
694       AliL3MemHandler::UpdateRowPointer(tempPt);
695     }
696   
697   cout<<"Doing the combinatorics"<<endl;
698   
699   AliL3LDigit *dPt1,*dPt2;
700   
701   for(Int_t e=0; e<GetNEtaSegments(); e++)
702     {
703       for(Int_t i=rowrange[0]; i<=rowrange[1]; i++)
704         {
705           Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e;
706           
707           for(dPt1 = (AliL3LDigit*)etaPt[index1].fFirst; dPt1 != 0; dPt1 = (AliL3LDigit*)dPt1->fNext)
708             {
709               for(Int_t j=i+1; j<=rowrange[1]; j++)
710                 {
711                   Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e;
712                   
713                   for(dPt2 = (AliL3LDigit*)etaPt[index2].fFirst; dPt2 != 0; dPt2 = (AliL3LDigit*)dPt2->fNext)
714                     {
715                       if(dPt1->fRow == dPt2->fRow)
716                         {
717                           cerr<<"same row; indexes "<<index1<<" "<<index2<<endl;
718                           exit(5);
719                         }
720                       
721                       //Get the correct histogrampointer:
722                       AliL3Histogram *hist = fParamSpace[e];
723                       if(!hist)
724                         {
725                           printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",i);
726                           continue;
727                         }
728                       
729                       //Do the transform:
730                       float x1 = AliL3Transform::Row2X(dPt1->fRow) - AliL3Transform::Row2X(rowrange[0]);
731                       float x2 = AliL3Transform::Row2X(dPt2->fRow) - AliL3Transform::Row2X(rowrange[0]);
732                       float y1 = dPt1->fY;
733                       float y2 = dPt2->fY;
734                       float theta = atan2(x2-x1,y1-y2);
735                       float rho = x1*cos(theta)+y1*sin(theta);
736                       hist->Fill(theta,rho,1);//dPt1->charge+dPt2->charge);
737                     }
738                 }
739             }
740         }
741     }
742
743   cout<<"done"<<endl;
744   delete [] etaPt;
745   delete [] digits;
746 }
747
748 #ifdef do_mc
749 Int_t AliL3HoughTransformer::GetTrackID(Int_t etaindex,Double_t kappa,Double_t psi) const
750 {
751   // Returns the MC label for a given peak found in the Hough space
752   if(!fDoMC)
753     {
754       cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
755       return -1;
756     }
757   
758   if(etaindex < 0 || etaindex > GetNEtaSegments())
759     {
760       cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<etaindex<<endl;
761       return -1;
762     }
763   AliL3Histogram *hist = fParamSpace[etaindex];
764   Int_t bin = hist->FindBin(kappa,psi);
765   Int_t label=-1;
766   Int_t max=0;
767   for(UInt_t i=0; i<MaxTrack; i++)
768     {
769       Int_t nhits=fTrackID[etaindex][bin].fNHits[i];
770       if(nhits == 0) break;
771       if(nhits > max)
772         {
773           max = nhits;
774           label = fTrackID[etaindex][bin].fLabel[i];
775         }
776     }
777   //nhits = max;
778   return label;
779 #else
780   Int_t AliL3HoughTransformer::GetTrackID(Int_t /*etaindex*/,Double_t /*kappa*/,Double_t /*psi*/) const
781 {
782   // Returns the MC label for a given peak found in the Hough space
783   if(!fDoMC)
784     {
785       cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
786       return -1;
787     }
788   
789   cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
790   return -1;
791 #endif
792 }
793