Little changes to make g++ version 3.2 compile the hough library.
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformer.cxx
1 //$Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ASV 
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 GCCVERSION == 3
16 using namespace std;
17 #endif
18
19 //_____________________________________________________________
20 // AliL3HoughTransformer
21 //
22 // Hough transformation class
23 //
24
25 ClassImp(AliL3HoughTransformer)
26
27 AliL3HoughTransformer::AliL3HoughTransformer()
28 {
29   //Default constructor
30   fParamSpace = 0;
31   fDoMC = kFALSE;
32 #ifdef do_mc
33   fTrackID = 0;
34 #endif
35 }
36
37 AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
38 {
39   //Normal constructor
40   fParamSpace = 0;
41   fDoMC = kFALSE;
42 #ifdef do_mc
43   fTrackID = 0;
44 #endif
45 }
46
47 AliL3HoughTransformer::~AliL3HoughTransformer()
48 {
49   DeleteHistograms();
50 #ifdef do_mc
51   if(fTrackID)
52     {
53       for(Int_t i=0; i<GetNEtaSegments(); i++)
54         {
55           if(!fTrackID[i]) continue;
56           delete fTrackID[i];
57         }
58       delete [] fTrackID;
59     }
60 #endif
61 }
62
63 //void AliL3HoughTransformer::Init(Int_t slice=0,Int_t patch=0,Int_t n_eta_segments=100){}
64
65 void AliL3HoughTransformer::DeleteHistograms()
66 {
67   if(!fParamSpace)
68     return;
69   for(Int_t i=0; i<GetNEtaSegments(); i++)
70     {
71       if(!fParamSpace[i]) continue;
72       delete fParamSpace[i];
73     }
74   delete [] fParamSpace;
75 }
76
77 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min,
78                                              Int_t nybin,Double_t phimin,Double_t phimax)
79 {
80   //Create the histograms (parameter space).
81   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
82   //The arguments give the range and binning; 
83   //nxbin = #bins in kappa
84   //nybin = #bins in phi0
85   //pt_min = mimium Pt of track (corresponding to maximum kappa)
86   //phi_min = mimimum phi0 (degrees)
87   //phi_max = maximum phi0 (degrees)
88     
89   Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
90   Double_t torad = AliL3Transform::Pi()/180;
91   CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
92 }
93
94 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
95                                              Int_t nybin,Double_t ymin,Double_t ymax)
96 {
97   
98   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
99   
100   Char_t histname[256];
101   for(Int_t i=0; i<GetNEtaSegments(); i++)
102     {
103       sprintf(histname,"paramspace_%d",i);
104       //fParamSpace[i] = new AliL3HistogramAdaptive(histname,0.1,1,0.05,64,ymin,ymax);
105       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
106     }
107   
108 #ifdef do_mc
109   if(fDoMC)
110     {
111       AliL3Histogram *hist = fParamSpace[0];
112       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
113       cout<<"Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
114       fTrackID = new TrackIndex*[GetNEtaSegments()];
115       for(Int_t i=0; i<GetNEtaSegments(); i++)
116         fTrackID[i] = new TrackIndex[ncells];
117     }
118 #endif
119 }
120
121 void AliL3HoughTransformer::Reset()
122 {
123   //Reset all the histograms
124
125   if(!fParamSpace)
126     {
127       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
128         <<"No histograms to reset"<<ENDLOG;
129       return;
130     }
131   
132   for(Int_t i=0; i<GetNEtaSegments(); i++)
133     fParamSpace[i]->Reset();
134 #ifdef do_mc
135   if(fDoMC)
136     {
137       AliL3Histogram *hist = fParamSpace[0];
138       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
139       for(Int_t i=0; i<GetNEtaSegments(); i++)
140         memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
141     }
142 #endif
143 }
144
145 Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta)
146 {
147   //Return the histogram index of the corresponding eta. 
148
149   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
150   Double_t index = (eta-GetEtaMin())/etaslice;
151   return (Int_t)index;
152 }
153
154 inline AliL3Histogram *AliL3HoughTransformer::GetHistogram(Int_t eta_index)
155 {
156   if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
157     return 0;
158   if(!fParamSpace[eta_index])
159     return 0;
160   return fParamSpace[eta_index];
161 }
162
163 Double_t AliL3HoughTransformer::GetEta(Int_t eta_index,Int_t slice)
164 {
165   Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
166   Double_t eta=(Double_t)((eta_index+0.5)*eta_slice);
167   if(slice>17) eta*=-1;
168   return eta;
169 }
170
171 void AliL3HoughTransformer::TransformCircle()
172 {
173   //Transform the input data with a circle HT.
174   //The function loops over all the data, and transforms each pixel with the equations:
175   // 
176   //kappa = 2/R*sin(phi - phi0)
177   //
178   //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
179   //
180   //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
181   //which histogram in which the pixel should be transformed, the eta-value is calcluated
182   //and the proper histogram index is found by GetEtaIndex(eta).
183
184
185   AliL3DigitRowData *tempPt = GetDataPointer();
186   if(!tempPt)
187     {
188       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
189         <<"No input data "<<ENDLOG;
190       return;
191     }
192   
193   //Loop over the padrows:
194   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
195     {
196       //Get the data on this padrow:
197       AliL3DigitData *digPt = tempPt->fDigitData;
198       if(i != (Int_t)tempPt->fRow)
199         {
200           cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl;
201           continue;
202         }
203
204       //Loop over the data on this padrow:
205       for(UInt_t j=0; j<tempPt->fNDigit; j++)
206         {
207           UShort_t charge = digPt[j].fCharge;
208           UChar_t pad = digPt[j].fPad;
209           UShort_t time = digPt[j].fTime;
210           if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
211             continue;
212           Int_t sector,row;
213           Float_t xyz[3];
214           
215           //Transform data to local cartesian coordinates:
216           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
217           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
218           
219           //Calculate the eta:
220           Double_t eta = AliL3Transform::GetEta(xyz);
221           
222           //Get the corresponding index, which determines which histogram to fill:
223           Int_t eta_index = GetEtaIndex(eta);
224           if(eta_index < 0 || eta_index >= GetNEtaSegments())
225             continue;
226           
227           //Get the correct histogrampointer:
228           AliL3Histogram *hist = fParamSpace[eta_index];
229           if(!hist)
230             {
231               printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index);
232               continue;
233             }
234
235           //Do the transformation:
236           Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]); 
237           Float_t phi = AliL3Transform::GetPhi(xyz);
238           
239           //Fill the histogram along the phirange
240           for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
241             {
242               Float_t phi0 = hist->GetBinCenterY(b);
243               Float_t kappa = 2*sin(phi - phi0)/R;
244               hist->Fill(kappa,phi0,charge);
245 #ifdef do_mc
246               if(fDoMC)
247                 {
248                   Int_t bin = hist->FindBin(kappa,phi0);
249                   for(Int_t t=0; t<3; t++)
250                     {
251                       Int_t label = digPt[j].fTrackID[t];
252                       if(label < 0) break;
253                       UInt_t c;
254                       for(c=0; c<MaxTrack; c++)
255                         if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0)
256                           break;
257                       if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
258                       fTrackID[eta_index][bin].fLabel[c] = label;
259                       fTrackID[eta_index][bin].fNHits[c]++;
260                     }
261                 }
262 #endif
263             }
264         }
265       
266       //Move the data pointer to the next padrow:
267       AliL3MemHandler::UpdateRowPointer(tempPt);
268     }
269 }
270
271 void AliL3HoughTransformer::TransformCircleC(Int_t row_range)
272 {
273   //Circle transform, using combinations of every 2 points lying
274   //on different padrows and within the same etaslice.
275   
276   AliL3DigitRowData *tempPt = GetDataPointer();
277   if(!tempPt)
278     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
279       <<"No input data "<<ENDLOG;
280  
281   Int_t counter=0;
282   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
283     {
284       counter += tempPt->fNDigit;
285       AliL3MemHandler::UpdateRowPointer(tempPt);
286     }
287   
288   struct Digit {
289     Int_t row;
290     Double_t r;
291     Double_t phi;
292     Int_t eta_index;
293     Int_t charge;
294     Int_t trackID[3];
295   };
296   
297   Digit *digits = new Digit[counter];
298   cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
299   
300   Int_t total_digits=counter;
301   Int_t sector,row,tot_charge,pad,time,charge;
302   Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
303   Float_t xyz[3];
304   
305   counter=0;
306   tempPt = GetDataPointer();
307   
308   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
309     {
310       AliL3DigitData *digPt = tempPt->fDigitData;
311       for(UInt_t di=0; di<tempPt->fNDigit; di++)
312         {
313           charge = digPt[di].fCharge;
314           pad = digPt[di].fPad;
315           time = digPt[di].fTime;
316           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
317           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
318           eta = AliL3Transform::GetEta(xyz);
319           digits[counter].row = i;
320           digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
321           digits[counter].phi = atan2(xyz[1],xyz[0]);
322           digits[counter].eta_index = GetEtaIndex(eta);
323           digits[counter].charge = charge;
324 #ifdef do_mc
325           if(fDoMC)
326             {
327               digits[counter].trackID[0] = digPt[di].fTrackID[0];
328               digits[counter].trackID[1] = digPt[di].fTrackID[1];
329               digits[counter].trackID[2] = digPt[di].fTrackID[2];
330             }
331 #endif
332           counter++;
333         }
334       AliL3MemHandler::UpdateRowPointer(tempPt);
335     }
336   
337   for(Int_t i=0; i<total_digits; i++)
338     {
339       if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue;
340       Int_t ind = digits[i].eta_index;
341       
342       for(Int_t j=i+1; j<total_digits; j++)
343         {
344           if(digits[i].row == digits[j].row) continue;
345           if(digits[i].eta_index != digits[j].eta_index) continue;
346           if(digits[i].row + row_range < digits[j].row) break;
347           
348           //Get the correct histogrampointer:
349           AliL3Histogram *hist = fParamSpace[ind];
350           if(!hist)
351             {
352               printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",ind);
353               continue;
354             }
355           
356           r1 = digits[i].r;
357           phi1 = digits[i].phi;
358           r2 = digits[j].r;
359           phi2 = digits[j].phi;
360           phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
361           kappa = 2*sin(phi2-phi_0)/r2;
362           tot_charge = digits[i].charge + digits[j].charge;
363           hist->Fill(kappa,phi_0,tot_charge);
364 #ifdef do_mc
365           if(fDoMC)
366             {
367               Int_t bin = hist->FindBin(kappa,phi_0);
368               for(Int_t l=0; l<3; l++)
369                 {
370                   for(Int_t m=0; m<3; m++)
371                     {
372                       if(digits[i].trackID[l] == digits[j].trackID[m])
373                         {
374                           Int_t label = digits[i].trackID[l];
375                           if(label < 0) continue;
376                           UInt_t c;
377                           for(c=0; c<MaxTrack; c++)
378                             if(fTrackID[ind][bin].fLabel[c] == label || fTrackID[ind][bin].fNHits[c] == 0)
379                               break;
380                           if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircleC : Array reached maximum!! "<<c<<endl;
381                           fTrackID[ind][bin].fLabel[c] = label;
382                           fTrackID[ind][bin].fNHits[c]++;
383                         }
384                     }
385                 }
386             }
387 #endif
388         }
389     }
390   delete [] digits;
391 }
392
393 void AliL3HoughTransformer::TransformLine()
394 {
395   //Do a line transform on the data.
396
397   
398   AliL3DigitRowData *tempPt = GetDataPointer();
399   if(!tempPt)
400     {
401       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data")
402         <<"No input data "<<ENDLOG;
403       return;
404     }
405     
406   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
407     {
408       AliL3DigitData *digPt = tempPt->fDigitData;
409       if(i != (Int_t)tempPt->fRow)
410         {
411           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
412           continue;
413         }
414       for(UInt_t j=0; j<tempPt->fNDigit; j++)
415         {
416           UShort_t charge = digPt[j].fCharge;
417           UChar_t pad = digPt[j].fPad;
418           UShort_t time = digPt[j].fTime;
419           if(charge < GetLowerThreshold())
420             continue;
421           Int_t sector,row;
422           Float_t xyz[3];
423           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
424           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
425           Float_t eta = AliL3Transform::GetEta(xyz);
426           Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
427           if(eta_index < 0 || eta_index >= GetNEtaSegments())
428             continue;
429           
430           //Get the correct histogram:
431           AliL3Histogram *hist = fParamSpace[eta_index];
432           if(!hist)
433             {
434               printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
435               continue;
436             }
437           for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
438             {
439               Double_t theta = hist->GetBinCenterX(xbin);
440               Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
441               hist->Fill(theta,rho,charge);
442             }
443         }
444       AliL3MemHandler::UpdateRowPointer(tempPt);
445     }
446   
447 }
448
449 Int_t AliL3HoughTransformer::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
450 {
451   if(!fDoMC)
452     {
453       cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
454       return -1;
455     }
456   
457 #ifdef do_mc
458   if(eta_index < 0 || eta_index > GetNEtaSegments())
459     {
460       cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl;
461       return -1;
462     }
463   AliL3Histogram *hist = fParamSpace[eta_index];
464   Int_t bin = hist->FindBin(kappa,psi);
465   Int_t label=-1;
466   Int_t max=0;
467   for(UInt_t i=0; i<MaxTrack; i++)
468     {
469       Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
470       if(nhits == 0) break;
471       if(nhits > max)
472         {
473           max = nhits;
474           label = fTrackID[eta_index][bin].fLabel[i];
475         }
476     }
477   return label;
478 #endif
479   cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
480   return -1;
481 }
482