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