0aac048497f53ce2b052af340ce8e0a45f647784
[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
15 #if GCCVERSION == 3
16 using namespace std;
17 #endif
18
19 /** \class AliL3HoughTransformerVhdl
20 // <pre>
21 //_____________________________________________________________
22 // AliL3HoughTransformerVhdl
23 //
24 // Hough transformation class for VHDL comparism.
25 //
26 //</pre>
27 */
28
29 ClassImp(AliL3HoughTransformerVhdl)
30
31
32 AliL3HoughTransformerVhdl::AliL3HoughTransformerVhdl()
33   : AliL3HoughTransformerLUT()
34
35 {
36   fEpsilon=0;
37   fSinEpsilon=0;
38   fCosEpsilon=1;
39   fIts=0;
40 }
41
42 AliL3HoughTransformerVhdl::AliL3HoughTransformerVhdl(Int_t slice,Int_t patch,Int_t n_eta_segments,Int_t n_its) 
43   : AliL3HoughTransformerLUT(slice,patch,n_eta_segments)
44 {
45   fEpsilon=0;
46   fSinEpsilon=0;
47   fCosEpsilon=1;
48   fIts=n_its;
49 }
50
51 AliL3HoughTransformerVhdl::~AliL3HoughTransformerVhdl()
52 {
53
54 }
55
56 void AliL3HoughTransformerVhdl::CreateHistograms(Int_t nxbin,Double_t pt_min,Int_t nybin,Double_t phimin,Double_t phimax)
57 {
58   AliL3HoughTransformerLUT::CreateHistograms(nxbin,pt_min,nybin,phimin,phimax);
59 }
60
61 void AliL3HoughTransformerVhdl::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
62                                                  Int_t nybin,Double_t ymin,Double_t ymax)
63 {
64   AliL3HoughTransformerLUT::CreateHistograms(nxbin,xmin,xmax,nybin,ymin,ymax);
65
66   fEpsilon=(ymax-ymin)/nybin;
67   fSinEpsilon=sin(fEpsilon);
68   fCosEpsilon=cos(fEpsilon);
69
70   //  cout << fEpsilon << " - " << (xmax-xmin)/nxbin << endl;
71 }
72
73 void AliL3HoughTransformerVhdl::Init(Int_t slice,Int_t patch,Int_t n_eta_segments,Int_t n_its)
74 {
75   AliL3HoughTransformerLUT::Init(slice,patch,n_eta_segments);
76 }
77
78 void AliL3HoughTransformerVhdl::TransformCircle()
79 {
80   //Transform the input data with a circle HT.
81   //The function loops over all the data, and transforms each pixel with the equation:
82   // 
83   //kappa      = lastkappa +- epsilon * lastkappaprime
84   //
85   //kappaprime = lastkappaprime -+ epsilon * lastkappa 
86   //
87   //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
88   //which histogram in which the pixel should be transformed, the eta-value is calcluated
89   //and the proper histogram index is found by GetEtaIndex(eta).
90
91   AliL3DigitRowData *tempPt = GetDataPointer();
92   if(!tempPt)
93     {
94       LOG(AliL3Log::kError,"AliL3HoughTransformerVhdl::TransformCircle","Data")
95         <<"No input data "<<ENDLOG;
96       return;
97     }
98
99 //#define use_error
100 #ifdef use_error
101       Float_t max_error=0;
102       Float_t rel_error=0;
103       Int_t counter=0;
104 #endif
105   
106   //Loop over the padrows:
107   for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
108     {
109       //Get the data on this padrow:
110       AliL3DigitData *digPt = tempPt->fDigitData;
111       if(i != (Int_t)tempPt->fRow)
112         {
113           LOG(AliL3Log::kError,"AliL3HoughTransformerVhdl::TransformCircle","Data")
114             <<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
115           continue;
116         }
117
118       //Loop over the data on this padrow:
119       for(UInt_t j=0; j<tempPt->fNDigit; j++)
120         {
121           UShort_t charge = digPt[j].fCharge;
122
123           //check threshold
124           if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
125             continue;
126
127           UChar_t pad = digPt[j].fPad;
128           UShort_t time = digPt[j].fTime;
129
130           Float_t x = CalcX(row);
131           Float_t y = CalcY(pad,row);
132           Float_t z = CalcZ(time);
133
134           //find eta slice
135           Float_t rz2 = 1 + (x*x+y*y)/(z*z);
136           Int_t eta_index = FindIndex(rz2);
137
138           if(eta_index < 0 || eta_index >= fNEtas){
139             //LOG(AliL3Log::kWarning,"AliL3HoughTransformerVhdl::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<eta_index<<"."<<ENDLOG;
140             continue;
141           }       
142
143           //Get the correct histogrampointer:
144           AliL3Histogram *hist = fParamSpace[eta_index];
145           if(!hist){
146             //LOG(AliL3Log::kWarning,"AliL3HoughTransformerVhdl::TransformCircle","Histograms")<<"Error getting histogram in index "<<eta_index<<"."<<ENDLOG;
147             continue;
148           }
149
150           //Fill the histogram along the phirange
151           Float_t kappa=0;
152           Float_t kappaprime=0;
153           Float_t lastkappa=0;
154           Float_t lastkappaprime=0;
155           Int_t its=0;
156           Float_t R2=1/(x*x+y*y);
157           Float_t A=2*y*R2;
158           Float_t B=-2*x*R2;
159
160 #if 1
161           //starting point
162           Float_t phi = fLUTphi0[fNPhi0/2];
163           kappa=A*cos(phi)+B*sin(phi);
164           kappaprime=B*cos(phi)-A*sin(phi);
165           its=fIts;
166           lastkappa=kappa;
167           lastkappaprime=kappaprime;
168           hist->Fill(kappa,phi,charge);
169
170           for(Int_t b=fNPhi0/2+1; b<fNPhi0; b++){ 
171             Float_t exact_kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
172
173             phi=fLUTphi0[b];
174
175             if(its==0) { //initialize or re-adjust values
176               kappa=A*cos(phi)+B*sin(phi); //equals exact_kappa!
177               kappaprime=B*cos(phi)-A*sin(phi);
178               its=fIts;
179             } else {
180               //kappa=fCosEpsilon*lastkappa+fSinEpsilon*lastkappaprime;
181               //kappaprime=fCosEpsilon*lastkappaprime-fSinEpsilon*lastkappa;
182               kappa=lastkappa+fEpsilon*lastkappaprime;
183               kappaprime=lastkappaprime-fEpsilon*lastkappa;
184               its--;
185             }
186
187             lastkappa=kappa;
188             lastkappaprime=kappaprime;
189
190             hist->Fill(kappa,phi,charge);
191           } // end positive running values
192
193           phi = fLUTphi0[fNPhi0/2];
194           kappa=A*cos(phi)+B*sin(phi);
195           kappaprime=B*cos(phi)-A*sin(phi);
196           its=fIts;
197           lastkappa=kappa;
198           lastkappaprime=kappaprime;
199           //hist->Fill(kappa,fLUTphi0[b],charge);
200
201           for(Int_t b=fNPhi0/2-1; b>=0; b--){ 
202             Float_t exact_kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
203
204             Float_t phi = fLUTphi0[b];
205
206             if(its==0) { //initialize or re-adjust values
207               kappa=A*cos(phi)+B*sin(phi); //equals exact_kappa!
208               kappaprime=B*cos(phi)-A*sin(phi);
209               its=fIts;
210             } else {
211               //kappa=fCosEpsilon*lastkappa-fSinEpsilon*lastkappaprime;
212               //kappaprime=fCosEpsilon*lastkappaprime+fSinEpsilon*lastkappa;
213               kappa=lastkappa-fEpsilon*lastkappaprime;
214               kappaprime=lastkappaprime+fEpsilon*lastkappa;
215               its--;
216             }
217
218             lastkappa=kappa;
219             lastkappaprime=kappaprime;
220
221             hist->Fill(kappa,phi,charge);
222           }
223
224 #else
225           for(Int_t b=0; b<fNPhi0; b++){ 
226             Float_t exact_kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
227
228             Float_t phi = fLUTphi0[b];
229
230             if(its==0) { //initialize or re-adjust values
231               kappa=A*cos(phi)+B*sin(phi); //equals exact_kappa!
232               kappaprime=B*cos(phi)-A*sin(phi);
233               its=fIts;
234             } else {
235               //kappa=fCosEpsilon*lastkappa+fSinEpsilon*lastkappaprime;
236               //kappaprime=fCosEpsilon*lastkappaprime-fSinEpsilon*lastkappa;
237               kappa=lastkappa+fEpsilon*lastkappaprime;
238               kappaprime=lastkappaprime-fEpsilon*lastkappa;
239               its--;
240             }
241
242 #ifdef use_error
243             Float_t err=fabs((exact_kappa-kappa)/exact_kappa);
244
245             if(err>max_error) {
246               //cout << max_error << " - " << err << " " << kappa << " " << exact_kappa << " " << kappa/exact_kappa<< endl;
247               max_error=err;
248             }
249
250             rel_error+=err;
251             counter++;
252 #endif
253
254             lastkappa=kappa;
255             lastkappaprime=kappaprime;
256
257             hist->Fill(kappa,fLUTphi0[b],charge);
258
259             //cout << kappa << " " << fLUTphi0[b] << " " << charge << endl;
260           }
261 #endif /*test*/
262
263         }
264       //Move the data pointer to the next padrow:
265       AliL3MemHandler::UpdateRowPointer(tempPt);
266     }
267
268 #ifdef use_error
269   cout <<"Max Error: " << max_error << endl;
270   cout <<"Rel Error average: " << rel_error/counter << endl;
271 #endif
272 }
273
274 void AliL3HoughTransformerVhdl::Print()
275 {
276   AliL3HoughTransformerLUT::Print();
277
278   cout << "fEpsilon: " << fEpsilon << endl;
279   cout << "fIts: " << fIts << endl;
280 }