1c218e523697d27fbe46577f783674fe8b4c72c2
[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 }
51
52 AliL3HoughTransformerLUT::AliL3HoughTransformerLUT(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
53 {
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;
74
75   Init(slice,patch,n_eta_segments);
76 }
77
78 AliL3HoughTransformerLUT::~AliL3HoughTransformerLUT()
79 {
80   DeleteHistograms();
81
82 #ifdef do_mc
83   if(fTrackID)
84     {
85       for(Int_t i=0; i<fNEtas; i++)
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   }
99
100   if(fNEtas){
101     delete[] fLUTEta;
102     fNEtas=0;
103   }
104 }
105
106 void 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){
119     for(Int_t i=0; i<fNEtas; i++)
120       {
121         if(!fParamSpace[i]) continue;
122         delete fParamSpace[i];
123       }
124     delete [] fParamSpace;
125   }
126 }
127
128 void 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;
148   if(fNEtas!=GetNEtaSegments()) {
149     cerr << "AliL3HoughTransformerLUT::Init -> Error: Number of Etas must be equal!" << endl;
150     exit(1);
151   }
152
153   AliL3Transform::Slice2Sector(slice,fMinRow,fSector,fSectorRow);
154   fZSign = slice < 18 ? 1:-1;
155   fSlice = slice;
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++){
181     Float_t eta=etamin_+(rr+1)*fEtaSlice;
182     fLUTEta[rr]=CalcRoverZ2(eta);
183     //cout << rr << ": " << eta << " " << fLUTEta[rr] << " " << GetEta(rr,fSlice) << endl;
184   }
185 }
186
187 void 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
203 void 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);
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++)
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;
237     fLUT2sinphi0[i]=2.*sin(phi0);
238     fLUT2cosphi0[i]=2.*cos(phi0);
239     //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
240   }  
241 }
242
243 void 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   
254   for(Int_t i=0; i<fNEtas; i++)
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);
261       for(Int_t i=0; i<fNEtas; i++)
262         memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
263     }
264 #endif
265 }
266
267 Int_t AliL3HoughTransformerLUT::GetEtaIndex(Double_t eta)
268 {
269   //Return the histogram index of the corresponding eta. 
270   
271   Float_t rz2=CalcRoverZ2(eta);
272   return FindIndex(rz2);
273 }
274
275 Int_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;
283 }
284
285 inline AliL3Histogram *AliL3HoughTransformerLUT::GetHistogram(Int_t eta_index)
286 {
287   if(!fParamSpace || eta_index >= fNEtas || eta_index < 0)
288     return 0;
289   if(!fParamSpace[eta_index])
290     return 0;
291   return fParamSpace[eta_index];
292 }
293
294 Float_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
302 Float_t AliL3HoughTransformerLUT::CalcEta(Float_t roverz2)
303 {
304   Float_t rz=sqrt(roverz2);
305   if(fZSign<0) rz=-rz;
306   Float_t ret=(1+rz)/(rz-1);
307   ret=0.5*log(ret);
308   return ret;
309 }
310
311 Float_t AliL3HoughTransformerLUT::CalcX(Int_t row)
312 {
313   return fLUTX[row];
314 }
315
316 Float_t AliL3HoughTransformerLUT::CalcY(Int_t pad,Int_t row)
317 {
318   return pad*fPadPitch-fLUTY[row];
319 }
320
321 Float_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
329 Double_t AliL3HoughTransformerLUT::GetEta(Int_t eta_index,Int_t slice)
330 {
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);
341 }
342
343 void AliL3HoughTransformerLUT::TransformCircle()
344 {
345   //Transform the input data with a circle HT.
346   //The function loops over all the data, and transforms each pixel with the equation:
347   // 
348   //kappa = 2/R^2 * (y*cos(phi0) - x*sin(phi0) )
349   //
350   //where R^2 = x^2 + y^2
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
356   AliL3DigitRowData *tempPt = GetDataPointer();
357   if(!tempPt)
358     {
359       LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
360         <<"No input data "<<ENDLOG;
361       return;
362     }
363   
364   //Loop over the padrows:
365   for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
366     {
367       //Get the data on this padrow:
368       AliL3DigitData *digPt = tempPt->fDigitData;
369       if(i != (Int_t)tempPt->fRow)
370         {
371           LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
372             <<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
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;
380
381           //check threshold
382           if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
383             continue;
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;
398             continue;
399           }       
400
401           //Get the correct histogrampointer:
402           AliL3Histogram *hist = fParamSpace[eta_index];
403           if(!hist){
404             //LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"Error getting histogram in index "<<eta_index<<"."<<ENDLOG;
405             continue;
406           }
407
408           //Fill the histogram along the phirange
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
416 #ifdef do_mcc
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               }
433 #endif
434           }
435         }
436       //Move the data pointer to the next padrow:
437       AliL3MemHandler::UpdateRowPointer(tempPt);
438     }
439 }
440
441 Int_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   
449 #ifdef do_mcc
450   if(eta_index < 0 || eta_index > fNEtas)
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
475 void 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 }