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