Changes done for new aliroot version. Faster calculation through saving LUT for kappa.
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformerVhdl.cxx
CommitLineData
b46b53c1 1//$Id$
2
3// Author: Constantin Loizides <mailto:loizides@fi.uib.no>
6173606e 4//*-- Copyright&Copy CL
5
e6cdf93b 6#include "AliL3StandardIncludes.h"
b46b53c1 7
45f0bc53 8#include "AliL3RootTypes.h"
b46b53c1 9#include "AliL3Logging.h"
e6cdf93b 10#include "AliL3MemHandler.h"
b46b53c1 11#include "AliL3Transform.h"
12#include "AliL3DigitData.h"
b46b53c1 13#include "AliL3HoughTransformerVhdl.h"
18659a6b 14#include "AliL3FFloat.h"
b46b53c1 15
e6cdf93b 16#if GCCVERSION == 3
17using namespace std;
18#endif
19
6173606e 20/** \class AliL3HoughTransformerVhdl
21// <pre>
b46b53c1 22//_____________________________________________________________
23// AliL3HoughTransformerVhdl
24//
6173606e 25// Hough transformation class for VHDL comparism.
b46b53c1 26//
6173606e 27//</pre>
28*/
b46b53c1 29
30ClassImp(AliL3HoughTransformerVhdl)
31
6173606e 32
b46b53c1 33AliL3HoughTransformerVhdl::AliL3HoughTransformerVhdl()
45f0bc53 34 : AliL3HoughTransformerLUT()
e6cdf93b 35
b46b53c1 36{
45f0bc53 37 fEpsilon=0;
38 fSinEpsilon=0;
39 fCosEpsilon=1;
40 fIts=0;
b46b53c1 41}
42
45f0bc53 43AliL3HoughTransformerVhdl::AliL3HoughTransformerVhdl(Int_t slice,Int_t patch,Int_t n_eta_segments,Int_t n_its)
44 : AliL3HoughTransformerLUT(slice,patch,n_eta_segments)
b46b53c1 45{
45f0bc53 46 fEpsilon=0;
47 fSinEpsilon=0;
48 fCosEpsilon=1;
49 fIts=n_its;
6173606e 50}
51
52AliL3HoughTransformerVhdl::~AliL3HoughTransformerVhdl()
53{
6173606e 54
6173606e 55}
56
57void AliL3HoughTransformerVhdl::CreateHistograms(Int_t nxbin,Double_t pt_min,Int_t nybin,Double_t phimin,Double_t phimax)
58{
45f0bc53 59 AliL3HoughTransformerLUT::CreateHistograms(nxbin,pt_min,nybin,phimin,phimax);
6173606e 60}
61
62void AliL3HoughTransformerVhdl::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
63 Int_t nybin,Double_t ymin,Double_t ymax)
64{
45f0bc53 65 AliL3HoughTransformerLUT::CreateHistograms(nxbin,xmin,xmax,nybin,ymin,ymax);
6173606e 66
45f0bc53 67 fEpsilon=(ymax-ymin)/nybin;
68 fSinEpsilon=sin(fEpsilon);
69 fCosEpsilon=cos(fEpsilon);
6173606e 70
18659a6b 71 fNxbin=nxbin;
72 fXmin=xmin;
73 fXmax=xmax;
74 fNybin=nybin;
75 fYmin=ymin;
76 fYmax=ymax;
45f0bc53 77 // cout << fEpsilon << " - " << (xmax-xmin)/nxbin << endl;
6173606e 78}
79
45f0bc53 80void AliL3HoughTransformerVhdl::Init(Int_t slice,Int_t patch,Int_t n_eta_segments,Int_t n_its)
6173606e 81{
45f0bc53 82 AliL3HoughTransformerLUT::Init(slice,patch,n_eta_segments);
6173606e 83}
84
85void AliL3HoughTransformerVhdl::TransformCircle()
86{
87 //Transform the input data with a circle HT.
45f0bc53 88 //The function loops over all the data, and transforms each pixel with the equation:
6173606e 89 //
45f0bc53 90 //kappa = lastkappa +- epsilon * lastkappaprime
b46b53c1 91 //
45f0bc53 92 //kappaprime = lastkappaprime -+ epsilon * lastkappa
b46b53c1 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
18659a6b 98
b46b53c1 99 AliL3DigitRowData *tempPt = GetDataPointer();
100 if(!tempPt)
101 {
102 LOG(AliL3Log::kError,"AliL3HoughTransformerVhdl::TransformCircle","Data")
103 <<"No input data "<<ENDLOG;
104 return;
105 }
45f0bc53 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
b46b53c1 113
114 //Loop over the padrows:
45f0bc53 115 for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
b46b53c1 116 {
117 //Get the data on this padrow:
118 AliL3DigitData *digPt = tempPt->fDigitData;
119 if(i != (Int_t)tempPt->fRow)
120 {
45f0bc53 121 LOG(AliL3Log::kError,"AliL3HoughTransformerVhdl::TransformCircle","Data")
b490510d 122 <<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
b46b53c1 123 continue;
124 }
6173606e 125
18659a6b 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
b46b53c1 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;
45f0bc53 138
139 //check threshold
140 if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
141 continue;
6173606e 142
b46b53c1 143 UChar_t pad = digPt[j].fPad;
144 UShort_t time = digPt[j].fTime;
6173606e 145
18659a6b 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
45f0bc53 231 Float_t z = CalcZ(time);
6173606e 232
45f0bc53 233 //find eta slice
234 Float_t rz2 = 1 + (x*x+y*y)/(z*z);
235 Int_t eta_index = FindIndex(rz2);
6173606e 236
45f0bc53 237 if(eta_index < 0 || eta_index >= fNEtas){
6173606e 238 //LOG(AliL3Log::kWarning,"AliL3HoughTransformerVhdl::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<eta_index<<"."<<ENDLOG;
b46b53c1 239 continue;
6173606e 240 }
45f0bc53 241
b46b53c1 242 //Get the correct histogrampointer:
243 AliL3Histogram *hist = fParamSpace[eta_index];
45f0bc53 244 if(!hist){
245 //LOG(AliL3Log::kWarning,"AliL3HoughTransformerVhdl::TransformCircle","Histograms")<<"Error getting histogram in index "<<eta_index<<"."<<ENDLOG;
246 continue;
247 }
b46b53c1 248
45f0bc53 249 for(Int_t b=0; b<fNPhi0; b++){
45f0bc53 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;
b46b53c1 257 }
45f0bc53 258
259 rel_error+=err;
260 counter++;
261#endif
18659a6b 262 //hist->Fill(kappa,fLUTphi0[b],charge);
45f0bc53 263
264 //cout << kappa << " " << fLUTphi0[b] << " " << charge << endl;
265 }
b46b53c1 266 }
b46b53c1 267 //Move the data pointer to the next padrow:
268 AliL3MemHandler::UpdateRowPointer(tempPt);
269 }
b46b53c1 270
45f0bc53 271#ifdef use_error
272 cout <<"Max Error: " << max_error << endl;
273 cout <<"Rel Error average: " << rel_error/counter << endl;
6173606e 274#endif
18659a6b 275
45f0bc53 276}
277
278void AliL3HoughTransformerVhdl::Print()
279{
280 AliL3HoughTransformerLUT::Print();
281
282 cout << "fEpsilon: " << fEpsilon << endl;
283 cout << "fIts: " << fIts << endl;
284}
18659a6b 285
286void 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}