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