Fast Hough transformer using extensivle LUT for geometry and cos/sin functions.
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformerLUT.cxx
CommitLineData
d5d754e5 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
16using namespace std;
17#endif
18
19//_____________________________________________________________
20// AliL3HoughTransformerLUT
21//
22// Hough transformation class using Look-UP-Tables
23//
24
25ClassImp(AliL3HoughTransformerLUT)
26
27AliL3HoughTransformerLUT::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
51AliL3HoughTransformerLUT::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
58AliL3HoughTransformerLUT::~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
85void 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
107void 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
160void 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
176void 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
216void 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
241Int_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
249inline 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
258Float_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
266Float_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
275Float_t AliL3HoughTransformerLUT::CalcX(Int_t row)
276{
277 return fLUTX[row];
278}
279
280Float_t AliL3HoughTransformerLUT::CalcY(Int_t pad,Int_t row)
281{
282 return pad*fPadPitch-fLUTY[row];
283}
284
285Float_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
293Double_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
301void 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
403Int_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
437void 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}