Minor bugfix
[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   fSlice=0;
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   fLUTKappa=0;
51   
52   fLastPad=0;
53   fLastIndex=0;
54 }
55
56 AliL3HoughTransformerLUT::AliL3HoughTransformerLUT(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
57 {
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;
78   fLUTKappa=0;
79
80   fLastPad=0;
81   fLastIndex=0;
82
83   Init(slice,patch,n_eta_segments);
84 }
85
86 AliL3HoughTransformerLUT::~AliL3HoughTransformerLUT()
87 {
88   DeleteHistograms();
89
90 #ifdef do_mc
91 #if 0
92   if(fTrackID)
93     {
94       for(Int_t i=0; i<fNEtas; i++)
95         {
96           if(!fTrackID[i]) continue;
97           delete fTrackID[i];
98         }
99       delete [] fTrackID;
100     }
101 #endif
102 #endif
103
104   if(fNRows){
105     delete[] fLUTX;
106     delete[] fLUTY;
107     fNRows=0;
108   }
109
110   if(fNEtas){
111     delete[] fLUTEta;
112     fNEtas=0;
113   }
114 }
115
116 void AliL3HoughTransformerLUT::DeleteHistograms()
117 {
118   if(fNPhi0){
119     delete[] fLUT2sinphi0;
120     delete[] fLUT2cosphi0;
121     delete[] fLUTphi0;
122     delete[] fLUTKappa;
123     fNPhi0=0;
124     fLUTphi0=0;
125     fLUT2sinphi0=0;
126     fLUT2cosphi0=0;
127     fLUTKappa=0;
128   }
129
130   if(!fParamSpace){
131     for(Int_t i=0; i<fNEtas; i++)
132       {
133         if(!fParamSpace[i]) continue;
134         delete fParamSpace[i];
135       }
136     delete [] fParamSpace;
137   }
138 }
139
140 void 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;
160   if(fNEtas!=GetNEtaSegments()) {
161     cerr << "AliL3HoughTransformerLUT::Init -> Error: Number of Etas must be equal!" << endl;
162     exit(1);
163   }
164
165   AliL3Transform::Slice2Sector(slice,fMinRow,fSector,fSectorRow);
166   fZSign = slice < 18 ? 1:-1;
167   fSlice = slice;
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++){
193     Float_t eta=etamin_+(rr+1)*fEtaSlice;
194     fLUTEta[rr]=CalcRoverZ2(eta);
195     //cout << rr << ": " << eta << " " << fLUTEta[rr] << " " << GetEta(rr,fSlice) << endl;
196   }
197 }
198
199 void 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
215 void 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     {
229       cerr << "Did not program changes in AliL3TransformerLUT for DO_MC!" << endl;
230       exit(1);
231 #if 0
232       AliL3Histogram *hist = fParamSpace[0];
233       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
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++)
237         fTrackID[i] = new TrackIndex[ncells];
238 #endif
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];
248   fLUTKappa=new Float_t[fNPhi0];
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;
254     fLUT2sinphi0[i]=2.*sin(phi0);
255     fLUT2cosphi0[i]=2.*cos(phi0);
256     fLUTKappa[i]=0.;
257     //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
258   }  
259 }
260
261 void 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   
272   for(Int_t i=0; i<fNEtas; i++)
273     fParamSpace[i]->Reset();
274 #ifdef do_mc
275   if(fDoMC)
276     {
277       cerr << "Did not program changes in AliL3TransformerLUT for DO_MC!" << endl;
278       exit(1);
279 #if 0
280       AliL3Histogram *hist = fParamSpace[0];
281       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
282       for(Int_t i=0; i<fNEtas; i++)
283         memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
284 #endif
285     }
286 #endif
287 }
288
289 Int_t AliL3HoughTransformerLUT::GetEtaIndex(Double_t eta)
290 {
291   //Return the histogram index of the corresponding eta. 
292   
293   Float_t rz2=CalcRoverZ2(eta);
294   return FindIndex(rz2);
295 }
296
297 Int_t AliL3HoughTransformerLUT::FindIndex(Float_t rz2, Int_t start)
298 {
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     }
311     index++;
312   }
313   //cout << start << " - " << index << ": " << rz2 << " " << fLUTEta[index] << endl;
314   return index;
315 }
316
317 inline AliL3Histogram *AliL3HoughTransformerLUT::GetHistogram(Int_t eta_index)
318 {
319   if(!fParamSpace || eta_index >= fNEtas || eta_index < 0)
320     return 0;
321   if(!fParamSpace[eta_index])
322     return 0;
323   return fParamSpace[eta_index];
324 }
325
326 Float_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
334 Float_t AliL3HoughTransformerLUT::CalcEta(Float_t roverz2)
335 {
336   Float_t rz=sqrt(roverz2);
337   if(fZSign<0) rz=-rz;
338   Float_t ret=(1+rz)/(rz-1);
339   ret=0.5*log(ret);
340   return ret;
341 }
342
343 Float_t AliL3HoughTransformerLUT::CalcX(Int_t row)
344 {
345   return fLUTX[row];
346 }
347
348 Float_t AliL3HoughTransformerLUT::CalcY(Int_t pad,Int_t row)
349 {
350   return pad*fPadPitch-fLUTY[row];
351 }
352
353 Float_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
361 Double_t AliL3HoughTransformerLUT::GetEta(Int_t eta_index,Int_t slice)
362 {
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);
373 }
374
375 void AliL3HoughTransformerLUT::TransformCircle()
376 {
377   //Transform the input data with a circle HT.
378   //The function loops over all the data, and transforms each pixel with the equation:
379   // 
380   //kappa = 2/R^2 * (y*cos(phi0) - x*sin(phi0) )
381   //
382   //where R^2 = x^2 + y^2
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
388   AliL3DigitRowData *tempPt = GetDataPointer();
389   if(!tempPt)
390     {
391       LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
392         <<"No input data "<<ENDLOG;
393       return;
394     }
395   
396   //Loop over the padrows:
397   for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
398     {
399       //Get the data on this padrow:
400       AliL3DigitData *digPt = tempPt->fDigitData;
401       if(i != (Int_t)tempPt->fRow)
402         {
403           LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
404             <<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
405           continue;
406         }
407
408       //calculate x for this row
409       Float_t x = CalcX(row);
410       Float_t x2=x*x;
411       Float_t r2=0;
412
413       //start a new row
414       fLastPad=-1;
415
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;
420
421           //check threshold
422           if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
423             continue;
424
425           UChar_t pad = digPt[j].fPad;
426           UShort_t time = digPt[j].fTime;
427
428           if(fLastPad!=pad){ //only update if necessary
429             fLastIndex=fNEtas-1;   
430   
431             //calculate hough for this pad
432             Float_t y = CalcY(pad,row);
433             Float_t y2 = y*y;
434             r2 = x2 + y2;
435             Float_t R2 = 1. / r2;
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
442           //find eta slice
443           Float_t z = CalcZ(time);
444           Float_t z2=z*z;
445           Float_t rz2 = 1 + r2/z2;
446           Int_t eta_index = FindIndex(rz2,fLastIndex);
447           //cout << row << " " << (int)pad << " " << (int)time << " " << eta_index <<" " <<fLastIndex<< endl;
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;
450             continue;
451           }       
452           //Get the correct histogrampointer:
453           AliL3Histogram *hist = fParamSpace[eta_index];
454           if(!hist){
455             //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"Error getting histogram in index "<<eta_index<<"."<<ENDLOG;
456             continue;
457           }
458           fLastIndex=eta_index;
459
460           //Fill the histogram along the phirange
461           for(Int_t b=0; b<fNPhi0; b++){
462             hist->Fill(fLUTKappa[b],fLUTphi0[b],charge);
463             //hist->Fill(fLUTKappa[b],fLUTphi0[b],1);
464             //cout << kappa << " " << fLUTphi0[b] << " " << charge << endl;
465
466 #ifdef do_mcc
467             if(fDoMC)
468               {
469                 cerr << "Did not program changes in AliL3TransformerLUT for DO_MC!" << endl;
470                 exit(1);
471 #if 0
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                   }
485 #endif
486               }
487 #endif
488           }
489         }
490       
491       //Move the data pointer to the next padrow:
492       AliL3MemHandler::UpdateRowPointer(tempPt);
493     }
494 }
495
496 Int_t AliL3HoughTransformerLUT::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
497 {
498   if(!fDoMC)
499     {
500       cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
501       return -1;
502     }
503   
504 #ifdef do_mcc
505   cerr << "Did not program changes in AliL3TransformerLUT for DO_MC!" << endl;
506   exit(1);
507 #if 0
508   if(eta_index < 0 || eta_index > fNEtas)
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;
528 #endif
529 #endif
530   cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
531   return -1;
532 }
533
534 void 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 }