Changes done for new aliroot version. Faster calculation through saving LUT for kappa.
[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;
3bcf971b 36 fSlice=0;
d5d754e5 37 fSectorRow=0;
38 fZSign=0;
39 fZLengthPlusOff=0;
40 fTimeWidth=0;
41 fPadPitch=0;
42 fEtaSlice=0;
43
44 fLUTX=0;
45 fLUTY=0;
46 fLUTEta=0;
47 fLUTphi0=0;
48 fLUT2sinphi0=0;
49 fLUT2cosphi0=0;
18659a6b 50 fLUTKappa=0;
51
52 fLastPad=0;
53 fLastIndex=0;
d5d754e5 54}
55
56AliL3HoughTransformerLUT::AliL3HoughTransformerLUT(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
57{
3bcf971b 58 fMinRow=0;
59 fMaxRow=0;
60
61 fNRows=0;
62 fNEtas=0;
63 fNPhi0=0;
64 fSector=0;
65 fSectorRow=0;
66 fZSign=0;
67 fZLengthPlusOff=0;
68 fTimeWidth=0;
69 fPadPitch=0;
70 fEtaSlice=0;
71
72 fLUTX=0;
73 fLUTY=0;
74 fLUTEta=0;
75 fLUTphi0=0;
76 fLUT2sinphi0=0;
77 fLUT2cosphi0=0;
18659a6b 78 fLUTKappa=0;
79
80 fLastPad=0;
81 fLastIndex=0;
d5d754e5 82
83 Init(slice,patch,n_eta_segments);
84}
85
86AliL3HoughTransformerLUT::~AliL3HoughTransformerLUT()
87{
88 DeleteHistograms();
89
90#ifdef do_mc
91 if(fTrackID)
92 {
3bcf971b 93 for(Int_t i=0; i<fNEtas; i++)
d5d754e5 94 {
95 if(!fTrackID[i]) continue;
96 delete fTrackID[i];
97 }
98 delete [] fTrackID;
99 }
100#endif
101
102 if(fNRows){
103 delete[] fLUTX;
104 delete[] fLUTY;
105 fNRows=0;
106 }
3bcf971b 107
d5d754e5 108 if(fNEtas){
109 delete[] fLUTEta;
110 fNEtas=0;
111 }
112}
113
114void AliL3HoughTransformerLUT::DeleteHistograms()
115{
116 if(fNPhi0){
117 delete[] fLUT2sinphi0;
118 delete[] fLUT2cosphi0;
119 delete[] fLUTphi0;
18659a6b 120 delete[] fLUTKappa;
d5d754e5 121 fNPhi0=0;
122 fLUTphi0=0;
123 fLUT2sinphi0=0;
124 fLUT2cosphi0=0;
18659a6b 125 fLUTKappa=0;
d5d754e5 126 }
127
128 if(!fParamSpace){
3bcf971b 129 for(Int_t i=0; i<fNEtas; i++)
d5d754e5 130 {
131 if(!fParamSpace[i]) continue;
132 delete fParamSpace[i];
133 }
134 delete [] fParamSpace;
135 }
136}
137
138void AliL3HoughTransformerLUT::Init(Int_t slice,Int_t patch,Int_t n_eta_segments)
139{
140 AliL3HoughBaseTransformer::Init(slice,patch,n_eta_segments);
141
142 //delete old LUT tables
143 if(fNRows){
144 fNRows=0;
145 delete[] fLUTX;
146 delete[] fLUTY;
147 }
148 if(fNEtas){
149 delete[] fLUTEta;
150 fNEtas=0;
151 }
152
153 //member values
154 fMinRow=AliL3Transform::GetFirstRow(patch);;
155 fMaxRow=AliL3Transform::GetLastRow(patch);;
156 fNRows=AliL3Transform::GetNRows(patch);;
157 fNEtas=n_eta_segments;
3bcf971b 158 if(fNEtas!=GetNEtaSegments()) {
159 cerr << "AliL3HoughTransformerLUT::Init -> Error: Number of Etas must be equal!" << endl;
160 exit(1);
161 }
d5d754e5 162
163 AliL3Transform::Slice2Sector(slice,fMinRow,fSector,fSectorRow);
164 fZSign = slice < 18 ? 1:-1;
3bcf971b 165 fSlice = slice;
d5d754e5 166 fZLengthPlusOff=AliL3Transform::GetZLength()+AliL3Transform::GetZOffset();
167 fTimeWidth=AliL3Transform::GetZWidth();
168
169 fPadPitch=0;
170 if(fSector<AliL3Transform::GetNSectorLow())
171 fPadPitch=AliL3Transform::GetPadPitchWidthLow();
172 else
173 fPadPitch=AliL3Transform::GetPadPitchWidthUp();
174
175 Float_t etamax_=GetEtaMax();
176 Float_t etamin_=GetEtaMin();
177 fEtaSlice=(etamax_-etamin_)/n_eta_segments;
178
179 //lookup tables for X and Y
180 fLUTX=new Float_t[fNRows];
181 fLUTY=new Float_t[fNRows];
182 for(Int_t rr=0;rr<fNRows;rr++){
183 fLUTX[rr]=Float_t(AliL3Transform::Row2X(rr+fMinRow));
184 fLUTY[rr]=Float_t(0.5*(AliL3Transform::GetNPads(rr+fMinRow)-1)*fPadPitch);
185 //cout << rr << ": " << (Float_t)fLUTX[rr] << " " << (Float_t)fLUTY[rr] << endl;
186 }
187
188 //lookup tables for rz2s <=> etas
189 fLUTEta=new Float_t[fNEtas];
190 for(Int_t rr=0;rr<fNEtas;rr++){
3bcf971b 191 Float_t eta=etamin_+(rr+1)*fEtaSlice;
192 fLUTEta[rr]=CalcRoverZ2(eta);
193 //cout << rr << ": " << eta << " " << fLUTEta[rr] << " " << GetEta(rr,fSlice) << endl;
d5d754e5 194 }
195}
196
197void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Double_t pt_min,Int_t nybin,Double_t phimin,Double_t phimax)
198{
199 //Create the histograms (parameter space).
200 //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
201 //The arguments give the range and binning;
202 //nxbin = #bins in kappa
203 //nybin = #bins in phi0
204 //pt_min = mimium Pt of track (corresponding to maximum kappa)
205 //phi_min = mimimum phi0 (degrees)
206 //phi_max = maximum phi0 (degrees)
207
208 Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
209 Double_t torad = AliL3Transform::Pi()/180;
210 CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
211}
212
213void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,Int_t nybin,Double_t ymin,Double_t ymax)
214{
215 fParamSpace = new AliL3Histogram*[fNEtas];
216
217 Char_t histname[256];
218 for(Int_t i=0; i<fNEtas; i++)
219 {
220 sprintf(histname,"paramspace_%d",i);
221 fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
222 }
223
224#ifdef do_mc
225 if(fDoMC)
226 {
227 AliL3Histogram *hist = fParamSpace[0];
228 Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
3bcf971b 229 cout<<"Allocating "<<fNEtas*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
230 fTrackID = new TrackIndex*[fNEtas];
231 for(Int_t i=0; i<fNEtas; i++)
d5d754e5 232 fTrackID[i] = new TrackIndex[ncells];
233 }
234#endif
235
236 //create lookup table for sin and cos
237 fNPhi0=nybin+1;
238
239 fLUTphi0=new Float_t[fNPhi0];
240 fLUT2sinphi0=new Float_t[fNPhi0];
241 fLUT2cosphi0=new Float_t[fNPhi0];
18659a6b 242 fLUTKappa=new Float_t[fNPhi0];
d5d754e5 243 Float_t diff=(ymax-ymin)/nybin;
244 Float_t phi0=ymin-0.5*diff;
245 for(Int_t i=0; i<fNPhi0; i++){
246 phi0+=diff;
247 fLUTphi0[i]=phi0;
3bcf971b 248 fLUT2sinphi0[i]=2.*sin(phi0);
249 fLUT2cosphi0[i]=2.*cos(phi0);
18659a6b 250 fLUTKappa[i]=0.;
3bcf971b 251 //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
d5d754e5 252 }
253}
254
255void AliL3HoughTransformerLUT::Reset()
256{
257 //Reset all the histograms
258
259 if(!fParamSpace)
260 {
261 LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
262 <<"No histograms to reset"<<ENDLOG;
263 return;
264 }
265
3bcf971b 266 for(Int_t i=0; i<fNEtas; i++)
d5d754e5 267 fParamSpace[i]->Reset();
268#ifdef do_mc
269 if(fDoMC)
270 {
271 AliL3Histogram *hist = fParamSpace[0];
272 Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
3bcf971b 273 for(Int_t i=0; i<fNEtas; i++)
d5d754e5 274 memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
275 }
276#endif
277}
278
d5d754e5 279Int_t AliL3HoughTransformerLUT::GetEtaIndex(Double_t eta)
280{
281 //Return the histogram index of the corresponding eta.
3bcf971b 282
283 Float_t rz2=CalcRoverZ2(eta);
284 return FindIndex(rz2);
285}
286
18659a6b 287Int_t AliL3HoughTransformerLUT::FindIndex(Float_t rz2, Int_t start)
3bcf971b 288{
18659a6b 289 //could improve search through devide and conquere strategy
290
291 Int_t index=start;
292 if(index==-100){
293 index=0;
294 while((index<fNEtas)&&(rz2<=fLUTEta[index])){
295 index++;
296 }
297 } else {
298 while((index>=0)&&(rz2>fLUTEta[index])){
299 index--;
300 }
3bcf971b 301 }
18659a6b 302 //cout << start << " - " << index << ": " << rz2 << " " << fLUTEta[index] << endl;
3bcf971b 303 return index;
d5d754e5 304}
305
306inline AliL3Histogram *AliL3HoughTransformerLUT::GetHistogram(Int_t eta_index)
307{
3bcf971b 308 if(!fParamSpace || eta_index >= fNEtas || eta_index < 0)
d5d754e5 309 return 0;
310 if(!fParamSpace[eta_index])
311 return 0;
312 return fParamSpace[eta_index];
313}
314
315Float_t AliL3HoughTransformerLUT::CalcRoverZ2(Float_t eta)
316{
317 Float_t e=exp(2*eta);
318 Float_t ret=(e+1)/(e-1);
319 ret*=ret;
320 return ret;
321}
322
323Float_t AliL3HoughTransformerLUT::CalcEta(Float_t roverz2)
324{
325 Float_t rz=sqrt(roverz2);
3bcf971b 326 if(fZSign<0) rz=-rz;
d5d754e5 327 Float_t ret=(1+rz)/(rz-1);
328 ret=0.5*log(ret);
329 return ret;
330}
331
332Float_t AliL3HoughTransformerLUT::CalcX(Int_t row)
333{
334 return fLUTX[row];
335}
336
337Float_t AliL3HoughTransformerLUT::CalcY(Int_t pad,Int_t row)
338{
339 return pad*fPadPitch-fLUTY[row];
340}
341
342Float_t AliL3HoughTransformerLUT::CalcZ(Int_t time)
343{
344 Float_t ret=time*fTimeWidth;
345 if(fZSign>0) ret=fZLengthPlusOff-ret;
346 else ret=ret-fZLengthPlusOff;
347 return ret;
348}
349
350Double_t AliL3HoughTransformerLUT::GetEta(Int_t eta_index,Int_t slice)
351{
3bcf971b 352 if(eta_index >= fNEtas || eta_index < 0){
353 //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::GetEta","Index") << "Index out of range."<<ENDLOG;
354 return 0.;
355 }
356 if(slice != fSlice){
357 //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::GetEta","Index") << "Given slice does not match internal slice."<<ENDLOG;
358 return 0.;
359 }
360
361 return (CalcEta(fLUTEta[eta_index])-0.5*fEtaSlice);
d5d754e5 362}
363
364void AliL3HoughTransformerLUT::TransformCircle()
365{
d5d754e5 366 //Transform the input data with a circle HT.
3bcf971b 367 //The function loops over all the data, and transforms each pixel with the equation:
d5d754e5 368 //
3bcf971b 369 //kappa = 2/R^2 * (y*cos(phi0) - x*sin(phi0) )
d5d754e5 370 //
3bcf971b 371 //where R^2 = x^2 + y^2
d5d754e5 372 //
373 //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
374 //which histogram in which the pixel should be transformed, the eta-value is calcluated
375 //and the proper histogram index is found by GetEtaIndex(eta).
376
d5d754e5 377 AliL3DigitRowData *tempPt = GetDataPointer();
378 if(!tempPt)
379 {
3bcf971b 380 LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
d5d754e5 381 <<"No input data "<<ENDLOG;
382 return;
383 }
384
385 //Loop over the padrows:
3bcf971b 386 for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
d5d754e5 387 {
388 //Get the data on this padrow:
389 AliL3DigitData *digPt = tempPt->fDigitData;
390 if(i != (Int_t)tempPt->fRow)
391 {
3bcf971b 392 LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
b490510d 393 <<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
d5d754e5 394 continue;
395 }
396
18659a6b 397 Float_t x = CalcX(row);
398 Float_t x2=x*x;
399 Float_t y=0,y2=0;
400 Float_t r2=0;
401 Float_t R2=0;
402 fLastPad=-1;
403
d5d754e5 404 //Loop over the data on this padrow:
405 for(UInt_t j=0; j<tempPt->fNDigit; j++)
406 {
407 UShort_t charge = digPt[j].fCharge;
3bcf971b 408
409 //check threshold
d5d754e5 410 if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
411 continue;
3bcf971b 412
413 UChar_t pad = digPt[j].fPad;
414 UShort_t time = digPt[j].fTime;
415
18659a6b 416 if(fLastPad!=pad){ //only update if necessary
417 fLastIndex=fNEtas-1;
418 y = CalcY(pad,row);
419 y2 = y*y;
420 r2 = x2 + y2;
421 R2 = 1. / r2;
422 for(Int_t b=0; b<fNPhi0; b++)
423 fLUTKappa[b]=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
424
425 fLastPad=pad;
426 }
427
3bcf971b 428 Float_t z = CalcZ(time);
429
430 //find eta slice
18659a6b 431 Float_t rz2 = 1 + r2/(z*z);
432 Int_t eta_index = FindIndex(rz2,fLastIndex);
433 //cout << row << " " << (int)pad << " " << (int)time << " " << eta_index <<" " <<fLastIndex<< endl;
434 fLastIndex=eta_index;
3bcf971b 435 if(eta_index < 0 || eta_index >= fNEtas){
436 //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<eta_index<<"."<<ENDLOG;
d5d754e5 437 continue;
3bcf971b 438 }
439
d5d754e5 440 //Get the correct histogrampointer:
441 AliL3Histogram *hist = fParamSpace[eta_index];
3bcf971b 442 if(!hist){
443 //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"Error getting histogram in index "<<eta_index<<"."<<ENDLOG;
444 continue;
445 }
446
d5d754e5 447 //Fill the histogram along the phirange
3bcf971b 448
18659a6b 449 for(Int_t b=0; b<fNPhi0; b++){
450 //Float_t kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
451 Float_t kappa=fLUTKappa[b];
3bcf971b 452 hist->Fill(kappa,fLUTphi0[b],charge);
453 //cout << kappa << " " << fLUTphi0[b] << " " << charge << endl;
454
b490510d 455#ifdef do_mcc
3bcf971b 456 if(fDoMC)
457 {
458 Int_t bin = hist->FindBin(kappa,phi0);
459 for(Int_t t=0; t<3; t++)
460 {
461 Int_t label = digPt[j].fTrackID[t];
462 if(label < 0) break;
463 UInt_t c;
464 for(c=0; c<MaxTrack; c++)
465 if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0)
466 break;
467 if(c == MaxTrack-1) LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","MCData") <<"Array reached maximum!! "<<c<<endl;
468 fTrackID[eta_index][bin].fLabel[c] = label;
469 fTrackID[eta_index][bin].fNHits[c]++;
470 }
471 }
d5d754e5 472#endif
3bcf971b 473 }
d5d754e5 474 }
d5d754e5 475 //Move the data pointer to the next padrow:
476 AliL3MemHandler::UpdateRowPointer(tempPt);
477 }
d5d754e5 478}
479
480Int_t AliL3HoughTransformerLUT::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
481{
482 if(!fDoMC)
483 {
484 cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
485 return -1;
486 }
487
b490510d 488#ifdef do_mcc
3bcf971b 489 if(eta_index < 0 || eta_index > fNEtas)
d5d754e5 490 {
491 cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl;
492 return -1;
493 }
494 AliL3Histogram *hist = fParamSpace[eta_index];
495 Int_t bin = hist->FindBin(kappa,psi);
496 Int_t label=-1;
497 Int_t max=0;
498 for(UInt_t i=0; i<MaxTrack; i++)
499 {
500 Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
501 if(nhits == 0) break;
502 if(nhits > max)
503 {
504 max = nhits;
505 label = fTrackID[eta_index][bin].fLabel[i];
506 }
507 }
508 return label;
509#endif
510 cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
511 return -1;
512}
513
514void AliL3HoughTransformerLUT::Print()
515{
516 cout << "fSlice: " << GetSlice() << endl;
517 cout << "fPatch: " << GetPatch() << endl;
518 cout << "fSector: " << fSector << endl;
519 cout << "fSectorRow: " << fSectorRow << endl;
520 cout << "fMinRow: " << fMinRow << endl;
521 cout << "fMaxRow: " << fMaxRow << endl;
522 cout << "fNRows: " << fNRows << endl;
523 cout << "fNEtas: " << fNEtas << endl;
524 cout << "fNPhi0: " << fNPhi0 << endl;
525 cout << "fZSign: " << fZSign << endl;
526 cout << "fZLengthPlusOff: " << fZLengthPlusOff << endl;
527 cout << "fPadPitch: " << fPadPitch << endl;
528 cout << "fTimeWidth: " << fTimeWidth << endl;
529
530 if(!fNRows) return;
531 cout << "fLUTX " << fNRows << endl;
532 for(Int_t i=0;i<fNRows;i++) cout << "fLUTX[" << i << "]=" << (Float_t)fLUTX[i] << endl;
533 cout << "fLUTY " << fNRows << endl;
534 for(Int_t i=0;i<fNRows;i++) cout << "fLUTY[" << i << "]=" << fLUTY[i] << endl;
535 if(!fNEtas) return;
536 cout << "fLUTEta " << fNEtas << endl;
537 for(Int_t i=0;i<fNEtas;i++) cout << "fLUTEta[" << i << "]=" << fLUTEta[i] << endl;
538 if(!fNPhi0) return;
539 cout << "fLUTphi0 " << fNPhi0 << endl;
540 for(Int_t i=0;i<fNPhi0;i++) cout << "fLUTPhi0[" << i << "]=" << fLUTphi0[i] << endl;
541 cout << "fLUT2sinphi0 " << fNPhi0 << endl;
542 for(Int_t i=0;i<fNPhi0;i++) cout << "fLUT2sinphi0[" << i << "]=" << fLUT2sinphi0[i] << endl;
543 cout << "fLUT2cosphi0 " << fNPhi0 << endl;
544 for(Int_t i=0;i<fNPhi0;i++) cout << "fLUT2cosphi0[" << i << "]=" << fLUT2cosphi0[i] << endl;
545}