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