Changes done for new aliroot version. Faster calculation through saving LUT for kappa.
[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(fTrackID)
92     {
93       for(Int_t i=0; i<fNEtas; i++)
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   }
107
108   if(fNEtas){
109     delete[] fLUTEta;
110     fNEtas=0;
111   }
112 }
113
114 void AliL3HoughTransformerLUT::DeleteHistograms()
115 {
116   if(fNPhi0){
117     delete[] fLUT2sinphi0;
118     delete[] fLUT2cosphi0;
119     delete[] fLUTphi0;
120     delete[] fLUTKappa;
121     fNPhi0=0;
122     fLUTphi0=0;
123     fLUT2sinphi0=0;
124     fLUT2cosphi0=0;
125     fLUTKappa=0;
126   }
127
128   if(!fParamSpace){
129     for(Int_t i=0; i<fNEtas; i++)
130       {
131         if(!fParamSpace[i]) continue;
132         delete fParamSpace[i];
133       }
134     delete [] fParamSpace;
135   }
136 }
137
138 void 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;
158   if(fNEtas!=GetNEtaSegments()) {
159     cerr << "AliL3HoughTransformerLUT::Init -> Error: Number of Etas must be equal!" << endl;
160     exit(1);
161   }
162
163   AliL3Transform::Slice2Sector(slice,fMinRow,fSector,fSectorRow);
164   fZSign = slice < 18 ? 1:-1;
165   fSlice = slice;
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++){
191     Float_t eta=etamin_+(rr+1)*fEtaSlice;
192     fLUTEta[rr]=CalcRoverZ2(eta);
193     //cout << rr << ": " << eta << " " << fLUTEta[rr] << " " << GetEta(rr,fSlice) << endl;
194   }
195 }
196
197 void 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
213 void 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);
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++)
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];
242   fLUTKappa=new Float_t[fNPhi0];
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;
248     fLUT2sinphi0[i]=2.*sin(phi0);
249     fLUT2cosphi0[i]=2.*cos(phi0);
250     fLUTKappa[i]=0.;
251     //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
252   }  
253 }
254
255 void 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   
266   for(Int_t i=0; i<fNEtas; i++)
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);
273       for(Int_t i=0; i<fNEtas; i++)
274         memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
275     }
276 #endif
277 }
278
279 Int_t AliL3HoughTransformerLUT::GetEtaIndex(Double_t eta)
280 {
281   //Return the histogram index of the corresponding eta. 
282   
283   Float_t rz2=CalcRoverZ2(eta);
284   return FindIndex(rz2);
285 }
286
287 Int_t AliL3HoughTransformerLUT::FindIndex(Float_t rz2, Int_t start)
288 {
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     }
301   }
302   //cout << start << " - " << index << ": " << rz2 << " " << fLUTEta[index] << endl;
303   return index;
304 }
305
306 inline AliL3Histogram *AliL3HoughTransformerLUT::GetHistogram(Int_t eta_index)
307 {
308   if(!fParamSpace || eta_index >= fNEtas || eta_index < 0)
309     return 0;
310   if(!fParamSpace[eta_index])
311     return 0;
312   return fParamSpace[eta_index];
313 }
314
315 Float_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
323 Float_t AliL3HoughTransformerLUT::CalcEta(Float_t roverz2)
324 {
325   Float_t rz=sqrt(roverz2);
326   if(fZSign<0) rz=-rz;
327   Float_t ret=(1+rz)/(rz-1);
328   ret=0.5*log(ret);
329   return ret;
330 }
331
332 Float_t AliL3HoughTransformerLUT::CalcX(Int_t row)
333 {
334   return fLUTX[row];
335 }
336
337 Float_t AliL3HoughTransformerLUT::CalcY(Int_t pad,Int_t row)
338 {
339   return pad*fPadPitch-fLUTY[row];
340 }
341
342 Float_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
350 Double_t AliL3HoughTransformerLUT::GetEta(Int_t eta_index,Int_t slice)
351 {
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);
362 }
363
364 void AliL3HoughTransformerLUT::TransformCircle()
365 {
366   //Transform the input data with a circle HT.
367   //The function loops over all the data, and transforms each pixel with the equation:
368   // 
369   //kappa = 2/R^2 * (y*cos(phi0) - x*sin(phi0) )
370   //
371   //where R^2 = x^2 + y^2
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
377   AliL3DigitRowData *tempPt = GetDataPointer();
378   if(!tempPt)
379     {
380       LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
381         <<"No input data "<<ENDLOG;
382       return;
383     }
384   
385   //Loop over the padrows:
386   for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
387     {
388       //Get the data on this padrow:
389       AliL3DigitData *digPt = tempPt->fDigitData;
390       if(i != (Int_t)tempPt->fRow)
391         {
392           LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
393             <<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
394           continue;
395         }
396
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
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;
408
409           //check threshold
410           if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
411             continue;
412
413           UChar_t pad = digPt[j].fPad;
414           UShort_t time = digPt[j].fTime;
415
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
428           Float_t z = CalcZ(time);
429
430           //find eta slice
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;
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;
437             continue;
438           }       
439
440           //Get the correct histogrampointer:
441           AliL3Histogram *hist = fParamSpace[eta_index];
442           if(!hist){
443             //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"Error getting histogram in index "<<eta_index<<"."<<ENDLOG;
444             continue;
445           }
446
447           //Fill the histogram along the phirange
448
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];
452             hist->Fill(kappa,fLUTphi0[b],charge);
453             //cout << kappa << " " << fLUTphi0[b] << " " << charge << endl;
454
455 #ifdef do_mcc
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               }
472 #endif
473           }
474         }
475       //Move the data pointer to the next padrow:
476       AliL3MemHandler::UpdateRowPointer(tempPt);
477     }
478 }
479
480 Int_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   
488 #ifdef do_mcc
489   if(eta_index < 0 || eta_index > fNEtas)
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
514 void 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 }