]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/hough/AliL3HoughTransformerVhdl.cxx
Merged Cvetans RowHoughTransformer, Anders latest developments in comp
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformerVhdl.cxx
CommitLineData
3e87ef69 1// @(#) $Id$
b46b53c1 2
3// Author: Constantin Loizides <mailto:loizides@fi.uib.no>
3e87ef69 4//*-- Copyright &copy ALICE HLT Group
6173606e 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
0bd0c1ef 16#if __GNUC__ == 3
e6cdf93b 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
b2a02bce 57void AliL3HoughTransformerVhdl::CreateHistograms(Int_t nxbin,Float_t pt_min,Int_t nybin,Float_t phimin,Float_t phimax)
6173606e 58{
45f0bc53 59 AliL3HoughTransformerLUT::CreateHistograms(nxbin,pt_min,nybin,phimin,phimax);
6173606e 60}
61
b2a02bce 62void AliL3HoughTransformerVhdl::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
63 Int_t nybin,Float_t ymin,Float_t ymax)
6173606e 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;
4ec77aa8 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++){
3e87ef69 175 //Float_t exact_kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
18659a6b 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;
3e87ef69 193 //hist->Fill(kappa,phi,charge);
18659a6b 194 } // end positive running values
195
196 phi = fLUTphi0[fNPhi0/2];
197 kappa=A*cos(phi)+B*sin(phi);
198 kappaprime=B*cos(phi)-A*sin(phi);
199 its=fIts;
200 lastkappa=kappa;
201 lastkappaprime=kappaprime;
202 //hist->Fill(kappa,fLUTphi0[b],charge);
203
204 for(Int_t b=fNPhi0/2-1; b>=0; b--){
3e87ef69 205 //Float_t exact_kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
18659a6b 206
207 Float_t phi = fLUTphi0[b];
208
209 if(its==0) { //initialize or re-adjust values
210 kappa=A*cos(phi)+B*sin(phi); //equals exact_kappa!
211 kappaprime=B*cos(phi)-A*sin(phi);
212 its=fIts;
213 } else {
214 //kappa=fCosEpsilon*lastkappa-fSinEpsilon*lastkappaprime;
215 //kappaprime=fCosEpsilon*lastkappaprime+fSinEpsilon*lastkappa;
216 kappa=lastkappa-fEpsilon*lastkappaprime;
217 kappaprime=lastkappaprime+fEpsilon*lastkappa;
218 its--;
219 }
220
221 lastkappa=kappa;
222 lastkappaprime=kappaprime;
3e87ef69 223 //hist->Fill(kappa,phi,charge);
18659a6b 224 }
225
226 fLastPad=pad;
227 }
228
45f0bc53 229 Float_t z = CalcZ(time);
6173606e 230
45f0bc53 231 //find eta slice
232 Float_t rz2 = 1 + (x*x+y*y)/(z*z);
233 Int_t eta_index = FindIndex(rz2);
6173606e 234
45f0bc53 235 if(eta_index < 0 || eta_index >= fNEtas){
6173606e 236 //LOG(AliL3Log::kWarning,"AliL3HoughTransformerVhdl::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<eta_index<<"."<<ENDLOG;
b46b53c1 237 continue;
6173606e 238 }
45f0bc53 239
b46b53c1 240 //Get the correct histogrampointer:
241 AliL3Histogram *hist = fParamSpace[eta_index];
45f0bc53 242 if(!hist){
243 //LOG(AliL3Log::kWarning,"AliL3HoughTransformerVhdl::TransformCircle","Histograms")<<"Error getting histogram in index "<<eta_index<<"."<<ENDLOG;
244 continue;
245 }
b46b53c1 246
45f0bc53 247 for(Int_t b=0; b<fNPhi0; b++){
45f0bc53 248
249#ifdef use_error
250 Float_t err=fabs((exact_kappa-kappa)/exact_kappa);
251
252 if(err>max_error) {
253 //cout << max_error << " - " << err << " " << kappa << " " << exact_kappa << " " << kappa/exact_kappa<< endl;
254 max_error=err;
b46b53c1 255 }
45f0bc53 256
257 rel_error+=err;
258 counter++;
259#endif
18659a6b 260 //hist->Fill(kappa,fLUTphi0[b],charge);
45f0bc53 261
262 //cout << kappa << " " << fLUTphi0[b] << " " << charge << endl;
263 }
b46b53c1 264 }
b46b53c1 265 //Move the data pointer to the next padrow:
266 AliL3MemHandler::UpdateRowPointer(tempPt);
267 }
b46b53c1 268
45f0bc53 269#ifdef use_error
270 cout <<"Max Error: " << max_error << endl;
271 cout <<"Rel Error average: " << rel_error/counter << endl;
6173606e 272#endif
18659a6b 273
45f0bc53 274}
275
276void AliL3HoughTransformerVhdl::Print()
277{
278 AliL3HoughTransformerLUT::Print();
279
280 cout << "fEpsilon: " << fEpsilon << endl;
281 cout << "fIts: " << fIts << endl;
282}
18659a6b 283
284void AliL3HoughTransformerVhdl::PrintVhdl()
285{
286 cout << "fSlice := " << GetSlice() << ";" << endl;
287 cout << "fPatch := " << GetPatch() << ";" << endl;
4ec77aa8 288 //cout << "fSector := " << fSector << ";" << endl;
289 //cout << "fSectorRow := " << fSectorRow << ";" << endl;
290 //cout << "fMinRow := " << fMinRow << ";" << endl;
291 //cout << "fMaxRow := " << fMaxRow << ";" << endl;
292 //cout << "fNRows := " << fNRows << ";" << endl;
293 //cout << "fNEtas := " << fNEtas << ";" << endl;
294 //cout << "fNPhi0 := " << fNPhi0 << ";" << endl;
18659a6b 295 cout << "fZSign := " << fZSign << ";" << endl;
296 cout << "fZLengthPlusOff := " << fZLengthPlusOff << ";" << endl;
297 cout << "fPadPitch := " << fPadPitch << ";" << endl;
298 cout << "fTimeWidth := " << fTimeWidth << ";" << endl;
299
300 if(!fNRows) return;
301 cout << "fNLUTX :=" << fNRows << ";" << endl;
302 cout << "fLUTX := ( ";
303 for(Int_t i=0;i<fNRows-1;i++) cout << fLUTX[i] << ", ";
304 cout << fLUTX[fNRows-1] << " );\n" << endl;
305 cout << "fNLUTY := " << fNRows << ";" << endl;
306 cout << "fLUTY := ( ";
307 for(Int_t i=0;i<fNRows-1;i++) cout << fLUTY[i] << ", ";
308 cout << fLUTY[fNRows-1] << " );\n" << endl;
309
310 if(!fNEtas) return;
311 cout << "fNLUTEta := " << fNEtas << ";" << endl;
312 cout << "fLUTEta := (";
313 for(Int_t i=0;i<fNEtas-1;i++) cout << fLUTEta[i] << ", ";
314 cout << fLUTEta[fNEtas-1] << " );\n" << endl;
315
316 if(!fNPhi0) return;
4ec77aa8 317 cout << "-- Kappa -- " << endl;
18659a6b 318 cout << "fNxbin := " << fNxbin << endl;
319 cout << "fxmin := " << fXmin << endl;
320 cout << "fxmax := " << fXmax << endl;
4ec77aa8 321 Float_t kdiff=(fXmax-fXmin)/fNxbin;
322 Float_t k=fXmin;
323 cout << "fKappa := ";
324 for(Int_t i=0;i<fNxbin;i++){
325 k+=kdiff;
326 cout << k << ", ";
327 }
328 cout << k+kdiff << " );\n" << endl;
329
330 cout << "-- Phi --" << endl;
18659a6b 331 cout << "fNybin := " << fNybin << endl;
332 cout << "fymin := " << fYmin << endl;
333 cout << "fymax := " << fYmax << endl;
334
335 cout << "\nfNLUTphi0 := " << fNPhi0 << ";" << endl;
336 cout << "fLUTphi0 := (";
337 for(Int_t i=0;i<fNPhi0-1;i++) cout << fLUTphi0[i] << ", ";
338 cout << fLUTphi0[fNPhi0-1] << " );\n" << endl;
339
4ec77aa8 340 cout << "\nfNLUT2sinphi0 := " << fNPhi0 << ";" << endl;
18659a6b 341 cout << "fLUT2sinphi0 := (";
342 for(Int_t i=0;i<fNPhi0-1;i++) cout << fLUT2sinphi0[i] << ", ";
343 cout << fLUT2sinphi0[fNPhi0-1] << " );\n" << endl;
344
4ec77aa8 345 cout << "\nfNLUT2cosphi0 := " << fNPhi0 << ";" << endl;
18659a6b 346 cout << "fLUT2cosphi0 := (";
3e87ef69 347 for(Int_t i=0;i<fNPhi0-1;i++) cout << fLUT2cosphi0[i] << ", ";
18659a6b 348 cout << fLUT2cosphi0[fNPhi0-1] << " );\n" << endl;
4ec77aa8 349
350 //cout << "\nfEpsilon := " << fEpsilon << endl;
351 //cout << "\nfSinEpsilon := " << fSinEpsilon << endl;
352 //cout << "\nfCosEpsilon := " << fCosEpsilon << endl;
18659a6b 353}