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