de7faf66e3d2ea0e50468c84a6dd90406159368e
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformerVhdl.cxx
1 //$Id$
2
3 // Author: Constantin Loizides <mailto:loizides@fi.uib.no>
4 //*-- Copyright&Copy CL
5
6 #include "AliL3StandardIncludes.h"
7
8 #include "AliL3RootTypes.h"
9 #include "AliL3Logging.h"
10 #include "AliL3MemHandler.h"
11 #include "AliL3Transform.h"
12 #include "AliL3DigitData.h"
13 #include "AliL3HoughTransformerVhdl.h"
14 #include "AliL3FFloat.h"
15
16 #if GCCVERSION == 3
17 using namespace std;
18 #endif
19
20 /** \class AliL3HoughTransformerVhdl
21 // <pre>
22 //_____________________________________________________________
23 // AliL3HoughTransformerVhdl
24 //
25 // Hough transformation class for VHDL comparism.
26 //
27 //</pre>
28 */
29
30 ClassImp(AliL3HoughTransformerVhdl)
31
32
33 AliL3HoughTransformerVhdl::AliL3HoughTransformerVhdl()
34   : AliL3HoughTransformerLUT()
35
36 {
37   fEpsilon=0;
38   fSinEpsilon=0;
39   fCosEpsilon=1;
40   fIts=0;
41 }
42
43 AliL3HoughTransformerVhdl::AliL3HoughTransformerVhdl(Int_t slice,Int_t patch,Int_t n_eta_segments,Int_t n_its) 
44   : AliL3HoughTransformerLUT(slice,patch,n_eta_segments)
45 {
46   fEpsilon=0;
47   fSinEpsilon=0;
48   fCosEpsilon=1;
49   fIts=n_its;
50 }
51
52 AliL3HoughTransformerVhdl::~AliL3HoughTransformerVhdl()
53 {
54
55 }
56
57 void AliL3HoughTransformerVhdl::CreateHistograms(Int_t nxbin,Double_t pt_min,Int_t nybin,Double_t phimin,Double_t phimax)
58 {
59   AliL3HoughTransformerLUT::CreateHistograms(nxbin,pt_min,nybin,phimin,phimax);
60 }
61
62 void AliL3HoughTransformerVhdl::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
63                                                  Int_t nybin,Double_t ymin,Double_t ymax)
64 {
65   AliL3HoughTransformerLUT::CreateHistograms(nxbin,xmin,xmax,nybin,ymin,ymax);
66
67   fEpsilon=(ymax-ymin)/nybin;
68   fSinEpsilon=sin(fEpsilon);
69   fCosEpsilon=cos(fEpsilon);
70
71   fNxbin=nxbin;
72   fXmin=xmin;
73   fXmax=xmax;
74   fNybin=nybin;
75   fYmin=ymin;
76   fYmax=ymax;
77   //  cout << fEpsilon << " - " << (xmax-xmin)/nxbin << endl;
78 }
79
80 void AliL3HoughTransformerVhdl::Init(Int_t slice,Int_t patch,Int_t n_eta_segments,Int_t n_its)
81 {
82   AliL3HoughTransformerLUT::Init(slice,patch,n_eta_segments);
83 }
84
85 void AliL3HoughTransformerVhdl::TransformCircle()
86 {
87   //Transform the input data with a circle HT.
88   //The function loops over all the data, and transforms each pixel with the equation:
89   // 
90   //kappa      = lastkappa +- epsilon * lastkappaprime
91   //
92   //kappaprime = lastkappaprime -+ epsilon * lastkappa 
93   //
94   //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
95   //which histogram in which the pixel should be transformed, the eta-value is calcluated
96   //and the proper histogram index is found by GetEtaIndex(eta).
97
98
99   AliL3DigitRowData *tempPt = GetDataPointer();
100   if(!tempPt)
101     {
102       LOG(AliL3Log::kError,"AliL3HoughTransformerVhdl::TransformCircle","Data")
103         <<"No input data "<<ENDLOG;
104       return;
105     }
106
107 //#define use_error
108 #ifdef use_error
109       Float_t max_error=0;
110       Float_t rel_error=0;
111       Int_t counter=0;
112 #endif
113   
114   //Loop over the padrows:
115   for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
116     {
117       //Get the data on this padrow:
118       AliL3DigitData *digPt = tempPt->fDigitData;
119       if(i != (Int_t)tempPt->fRow)
120         {
121           LOG(AliL3Log::kError,"AliL3HoughTransformerVhdl::TransformCircle","Data")
122             <<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
123           continue;
124         }
125
126       Float_t x = CalcX(row);
127       Float_t x2=x*x;
128       Float_t y=0,y2=0;
129       Float_t r2=0;
130       Float_t R2=0;
131       fLastPad=-1;
132
133
134       //Loop over the data on this padrow:
135       for(UInt_t j=0; j<tempPt->fNDigit; j++)
136         {
137           UShort_t charge = digPt[j].fCharge;
138
139           //check threshold
140           if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
141             continue;
142
143           UChar_t pad = digPt[j].fPad;
144           UShort_t time = digPt[j].fTime;
145
146           if(fLastPad!=pad){ //only update if necessary
147             fLastIndex=fNEtas-1;   
148             y = CalcY(pad,row);
149             y2 = y*y;
150             r2 = x2 + y2;
151             R2 = 1. / r2;
152             for(Int_t b=0; b<fNPhi0; b++)
153               fLUTKappa[b]=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
154
155             //Fill the histogram along the phirange
156             Float_t kappa=0;
157             Float_t kappaprime=0;
158             Float_t lastkappa=0;
159             Float_t lastkappaprime=0;
160             Int_t its=0;
161             Float_t R2=1/(x*x+y*y);
162             Float_t A=2*y*R2;
163             Float_t B=-2*x*R2;
164
165             //starting point
166             Float_t phi = fLUTphi0[fNPhi0/2];
167             kappa=A*cos(phi)+B*sin(phi);
168             kappaprime=B*cos(phi)-A*sin(phi);
169             its=fIts;
170             lastkappa=kappa;
171             lastkappaprime=kappaprime;
172             //    hist->Fill(kappa,phi,charge);
173
174             for(Int_t b=fNPhi0/2+1; b<fNPhi0; b++){ 
175               Float_t exact_kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
176               
177               phi=fLUTphi0[b];
178
179               if(its==0) { //initialize or re-adjust values
180                 kappa=A*cos(phi)+B*sin(phi); //equals exact_kappa!
181                 kappaprime=B*cos(phi)-A*sin(phi);
182                 its=fIts;
183               } else {
184                 //kappa=fCosEpsilon*lastkappa+fSinEpsilon*lastkappaprime;
185                 //kappaprime=fCosEpsilon*lastkappaprime-fSinEpsilon*lastkappa;
186                 kappa=lastkappa+fEpsilon*lastkappaprime;
187                 kappaprime=lastkappaprime-fEpsilon*lastkappa;
188                 its--;
189               }
190
191               lastkappa=kappa;
192               lastkappaprime=kappaprime;
193
194               //            hist->Fill(kappa,phi,charge);
195             } // end positive running values
196
197             phi = fLUTphi0[fNPhi0/2];
198             kappa=A*cos(phi)+B*sin(phi);
199             kappaprime=B*cos(phi)-A*sin(phi);
200             its=fIts;
201             lastkappa=kappa;
202             lastkappaprime=kappaprime;
203             //hist->Fill(kappa,fLUTphi0[b],charge);
204
205             for(Int_t b=fNPhi0/2-1; b>=0; b--){ 
206               Float_t exact_kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
207
208               Float_t phi = fLUTphi0[b];
209
210               if(its==0) { //initialize or re-adjust values
211                 kappa=A*cos(phi)+B*sin(phi); //equals exact_kappa!
212                 kappaprime=B*cos(phi)-A*sin(phi);
213                 its=fIts;
214               } else {
215                 //kappa=fCosEpsilon*lastkappa-fSinEpsilon*lastkappaprime;
216                 //kappaprime=fCosEpsilon*lastkappaprime+fSinEpsilon*lastkappa;
217                 kappa=lastkappa-fEpsilon*lastkappaprime;
218                 kappaprime=lastkappaprime+fEpsilon*lastkappa;
219                 its--;
220               }
221
222               lastkappa=kappa;
223               lastkappaprime=kappaprime;
224               
225               //            hist->Fill(kappa,phi,charge);
226             }
227
228             fLastPad=pad;
229           }
230
231           Float_t z = CalcZ(time);
232
233           //find eta slice
234           Float_t rz2 = 1 + (x*x+y*y)/(z*z);
235           Int_t eta_index = FindIndex(rz2);
236
237           if(eta_index < 0 || eta_index >= fNEtas){
238             //LOG(AliL3Log::kWarning,"AliL3HoughTransformerVhdl::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<eta_index<<"."<<ENDLOG;
239             continue;
240           }       
241
242           //Get the correct histogrampointer:
243           AliL3Histogram *hist = fParamSpace[eta_index];
244           if(!hist){
245             //LOG(AliL3Log::kWarning,"AliL3HoughTransformerVhdl::TransformCircle","Histograms")<<"Error getting histogram in index "<<eta_index<<"."<<ENDLOG;
246             continue;
247           }
248
249           for(Int_t b=0; b<fNPhi0; b++){ 
250
251 #ifdef use_error
252             Float_t err=fabs((exact_kappa-kappa)/exact_kappa);
253
254             if(err>max_error) {
255               //cout << max_error << " - " << err << " " << kappa << " " << exact_kappa << " " << kappa/exact_kappa<< endl;
256               max_error=err;
257             }
258
259             rel_error+=err;
260             counter++;
261 #endif
262             //hist->Fill(kappa,fLUTphi0[b],charge);
263
264             //cout << kappa << " " << fLUTphi0[b] << " " << charge << endl;
265           }
266         }
267       //Move the data pointer to the next padrow:
268       AliL3MemHandler::UpdateRowPointer(tempPt);
269     }
270
271 #ifdef use_error
272   cout <<"Max Error: " << max_error << endl;
273   cout <<"Rel Error average: " << rel_error/counter << endl;
274 #endif
275
276 }
277
278 void AliL3HoughTransformerVhdl::Print()
279 {
280   AliL3HoughTransformerLUT::Print();
281
282   cout << "fEpsilon: " << fEpsilon << endl;
283   cout << "fIts: " << fIts << endl;
284 }
285
286 void AliL3HoughTransformerVhdl::PrintVhdl()
287 {
288   cout << "fSlice := " << GetSlice() << ";" << endl;
289   cout << "fPatch := " << GetPatch() << ";" << endl;
290   //  cout << "fSector := " << fSector << ";" << endl;
291   //  cout << "fSectorRow := " << fSectorRow << ";" << endl;
292   //  cout << "fMinRow := " << fMinRow << ";" << endl;
293   //  cout << "fMaxRow := " << fMaxRow << ";" << endl;
294   //  cout << "fNRows := " << fNRows << ";" << endl;
295   //  cout << "fNEtas := " << fNEtas << ";" << endl;
296   //  cout << "fNPhi0 := " << fNPhi0 << ";" << endl;
297   cout << "fZSign := " << fZSign << ";" << endl;
298   cout << "fZLengthPlusOff := " << fZLengthPlusOff << ";" << endl;
299   cout << "fPadPitch := " << fPadPitch << ";" << endl;
300   cout << "fTimeWidth := " << fTimeWidth << ";" << endl;
301   
302   if(!fNRows) return;
303   cout << "fNLUTX :=" << fNRows << ";" << endl;
304   cout << "fLUTX := ( ";
305   for(Int_t i=0;i<fNRows-1;i++) cout << fLUTX[i] << ", ";
306   cout << fLUTX[fNRows-1] << " );\n" << endl; 
307   cout << "fNLUTY := " << fNRows << ";" << endl;
308   cout << "fLUTY := ( ";
309   for(Int_t i=0;i<fNRows-1;i++) cout << fLUTY[i] << ", ";
310   cout << fLUTY[fNRows-1] << " );\n" << endl;
311  
312   if(!fNEtas) return;
313   cout << "fNLUTEta := " << fNEtas << ";" << endl;
314   cout << "fLUTEta := ("; 
315   for(Int_t i=0;i<fNEtas-1;i++) cout << fLUTEta[i] << ", ";
316   cout << fLUTEta[fNEtas-1] << " );\n" << endl; 
317
318   if(!fNPhi0) return;
319   cout << "fNxbin := " << fNxbin << endl;
320   cout << "fxmin := " << fXmin << endl;
321   cout << "fxmax := " << fXmax << endl;
322   cout << "fNybin := " << fNybin << endl;
323   cout << "fymin := " << fYmin << endl;
324   cout << "fymax := " << fYmax << endl;
325
326   cout << "\nfNLUTphi0 := " << fNPhi0 << ";" << endl;
327   cout << "fLUTphi0 := (";
328   for(Int_t i=0;i<fNPhi0-1;i++) cout << fLUTphi0[i] << ", ";
329   cout << fLUTphi0[fNPhi0-1] << " );\n" << endl; 
330
331   cout << "fLUT2sinphi0 := (";
332   for(Int_t i=0;i<fNPhi0-1;i++) cout << fLUT2sinphi0[i] << ", ";
333   cout << fLUT2sinphi0[fNPhi0-1] << " );\n" << endl; 
334
335   cout << "fLUT2cosphi0 := (";
336   for(Int_t i=0;i<fNPhi0;i++) cout << fLUT2cosphi0[i] << ", ";
337   cout << fLUT2cosphi0[fNPhi0-1] << " );\n" << endl; 
338 }