L3 becomes HLT
[u/mrichter/AliRoot.git] / HLT / hough / AliHLTHoughTransformerLUT.cxx
CommitLineData
3e87ef69 1// @(#) $Id$
d5d754e5 2
3// Author: Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
3e87ef69 4//*-- Copyright &copy ALICE HLT Group
4aa41877 5/** /class AliHLTHoughTransformerLUT
f5271e68 6//<pre>
7//_____________________________________________________________
4aa41877 8// AliHLTHoughTransformerLUT
f5271e68 9//
10// Hough transformation class using Look-UP-Tables
11//
12//</pre>
13*/
d5d754e5 14
4aa41877 15#include "AliHLTStandardIncludes.h"
d5d754e5 16
4aa41877 17#include "AliHLTLogging.h"
18#include "AliHLTHoughTransformerLUT.h"
19#include "AliHLTTransform.h"
20#include "AliHLTMemHandler.h"
21#include "AliHLTDigitData.h"
22#include "AliHLTHistogram.h"
23#include "AliHLTHistogramAdaptive.h"
d5d754e5 24
c6747af6 25#if __GNUC__ >= 3
d5d754e5 26using namespace std;
27#endif
28
4aa41877 29ClassImp(AliHLTHoughTransformerLUT)
d5d754e5 30
4aa41877 31AliHLTHoughTransformerLUT::AliHLTHoughTransformerLUT() : AliHLTHoughBaseTransformer()
d5d754e5 32{
f5271e68 33 // default contructor
1f1942b8 34 fParamSpace=0;
d5d754e5 35 fMinRow=0;
36 fMaxRow=0;
37
38 fNRows=0;
39 fNEtas=0;
40 fNPhi0=0;
41 fSector=0;
3bcf971b 42 fSlice=0;
d5d754e5 43 fSectorRow=0;
44 fZSign=0;
45 fZLengthPlusOff=0;
46 fTimeWidth=0;
47 fPadPitch=0;
48 fEtaSlice=0;
49
50 fLUTX=0;
51 fLUTY=0;
52 fLUTEta=0;
3e87ef69 53 fLUTEtaReal=0;
d5d754e5 54 fLUTphi0=0;
55 fLUT2sinphi0=0;
56 fLUT2cosphi0=0;
18659a6b 57 fLUTKappa=0;
700db259 58
18659a6b 59 fLastPad=0;
60 fLastIndex=0;
3e87ef69 61 fAccCharge=0;
1f1942b8 62 fX=fY=0.;
d5d754e5 63}
64
4aa41877 65AliHLTHoughTransformerLUT::AliHLTHoughTransformerLUT(Int_t slice,Int_t patch,Int_t nEtaSegments) : AliHLTHoughBaseTransformer(slice,patch,nEtaSegments)
d5d754e5 66{
f5271e68 67 // constructor
1f1942b8 68 fParamSpace=0;
3bcf971b 69 fMinRow=0;
70 fMaxRow=0;
71
72 fNRows=0;
73 fNEtas=0;
74 fNPhi0=0;
75 fSector=0;
76 fSectorRow=0;
77 fZSign=0;
78 fZLengthPlusOff=0;
79 fTimeWidth=0;
80 fPadPitch=0;
81 fEtaSlice=0;
82
83 fLUTX=0;
84 fLUTY=0;
85 fLUTEta=0;
3e87ef69 86 fLUTEtaReal=0;
3bcf971b 87 fLUTphi0=0;
88 fLUT2sinphi0=0;
89 fLUT2cosphi0=0;
18659a6b 90 fLUTKappa=0;
91
92 fLastPad=0;
93 fLastIndex=0;
3e87ef69 94 fAccCharge=0;
1f1942b8 95 fX=fY=0.;
d5d754e5 96
f5271e68 97 Init(slice,patch,nEtaSegments);
d5d754e5 98}
99
4aa41877 100AliHLTHoughTransformerLUT::~AliHLTHoughTransformerLUT()
d5d754e5 101{
f5271e68 102 // destructor
d5d754e5 103 DeleteHistograms();
104
d5d754e5 105 if(fNRows){
106 delete[] fLUTX;
107 delete[] fLUTY;
108 fNRows=0;
109 }
3bcf971b 110
d5d754e5 111 if(fNEtas){
112 delete[] fLUTEta;
3e87ef69 113 delete[] fLUTEtaReal;
d5d754e5 114 fNEtas=0;
115 }
116}
117
4aa41877 118void AliHLTHoughTransformerLUT::DeleteHistograms()
d5d754e5 119{
f5271e68 120 // deletes all histograms
d5d754e5 121 if(fNPhi0){
122 delete[] fLUT2sinphi0;
123 delete[] fLUT2cosphi0;
124 delete[] fLUTphi0;
18659a6b 125 delete[] fLUTKappa;
d5d754e5 126 fNPhi0=0;
127 fLUTphi0=0;
128 fLUT2sinphi0=0;
129 fLUT2cosphi0=0;
18659a6b 130 fLUTKappa=0;
d5d754e5 131 }
132
133 if(!fParamSpace){
3bcf971b 134 for(Int_t i=0; i<fNEtas; i++)
d5d754e5 135 {
136 if(!fParamSpace[i]) continue;
137 delete fParamSpace[i];
138 }
139 delete [] fParamSpace;
140 }
141}
142
4aa41877 143void AliHLTHoughTransformerLUT::Init(Int_t slice,Int_t patch,Int_t nEtaSegments,Int_t /*nSeqs*/)
d5d754e5 144{
f5271e68 145 // Initialization
4aa41877 146 AliHLTHoughBaseTransformer::Init(slice,patch,nEtaSegments);
d5d754e5 147
148 //delete old LUT tables
149 if(fNRows){
150 fNRows=0;
151 delete[] fLUTX;
152 delete[] fLUTY;
153 }
154 if(fNEtas){
155 delete[] fLUTEta;
3e87ef69 156 delete[] fLUTEtaReal;
d5d754e5 157 fNEtas=0;
158 }
159
160 //member values
4aa41877 161 fMinRow=AliHLTTransform::GetFirstRow(patch);;
162 fMaxRow=AliHLTTransform::GetLastRow(patch);;
163 fNRows=AliHLTTransform::GetNRows(patch);;
f5271e68 164 fNEtas=nEtaSegments;
3bcf971b 165 if(fNEtas!=GetNEtaSegments()) {
4aa41877 166 cerr << "AliHLTHoughTransformerLUT::Init -> Error: Number of Etas must be equal!" << endl;
3bcf971b 167 exit(1);
168 }
d5d754e5 169
4aa41877 170 AliHLTTransform::Slice2Sector(slice,fMinRow,fSector,fSectorRow);
171 fZSign = slice < 18 ? 1:-1; //see AliHLTTransformer
3bcf971b 172 fSlice = slice;
4aa41877 173 fZLengthPlusOff=AliHLTTransform::GetZLength()+AliHLTTransform::GetZOffset();
174 fTimeWidth=AliHLTTransform::GetZWidth();
d5d754e5 175
176 fPadPitch=0;
4aa41877 177 if(fSector<AliHLTTransform::GetNSectorLow())
178 fPadPitch=AliHLTTransform::GetPadPitchWidthLow();
d5d754e5 179 else
4aa41877 180 fPadPitch=AliHLTTransform::GetPadPitchWidthUp();
d5d754e5 181
f5271e68 182 Float_t etaMax=GetEtaMax();
183 Float_t etaMin=GetEtaMin();
184 fEtaSlice=(etaMax-etaMin)/nEtaSegments;
d5d754e5 185
186 //lookup tables for X and Y
187 fLUTX=new Float_t[fNRows];
188 fLUTY=new Float_t[fNRows];
189 for(Int_t rr=0;rr<fNRows;rr++){
4aa41877 190 fLUTX[rr]=Float_t(AliHLTTransform::Row2X(rr+fMinRow));
191 fLUTY[rr]=Float_t(0.5*(AliHLTTransform::GetNPads(rr+fMinRow)-1)*fPadPitch);
d5d754e5 192 //cout << rr << ": " << (Float_t)fLUTX[rr] << " " << (Float_t)fLUTY[rr] << endl;
193 }
194
195 //lookup tables for rz2s <=> etas
1f1942b8 196 //will only be used ifdef FULLLUT
d5d754e5 197 fLUTEta=new Float_t[fNEtas];
3e87ef69 198 fLUTEtaReal=new Float_t[fNEtas];
d5d754e5 199 for(Int_t rr=0;rr<fNEtas;rr++){
f5271e68 200 Float_t eta=etaMin+(rr+1)*fEtaSlice;
3bcf971b 201 fLUTEta[rr]=CalcRoverZ2(eta);
3e87ef69 202 fLUTEtaReal[rr]=eta-0.5*fEtaSlice;
203 //cout << rr << ": " << eta << " " << fLUTEtaReal[rr] << " " << GetEta(rr,fSlice) << " - " << fLUTEta[rr] << endl;
d5d754e5 204 }
205}
206
4aa41877 207void AliHLTHoughTransformerLUT::CreateHistograms(Float_t ptmin,Float_t ptmax,Float_t ptres,
1f1942b8 208 Int_t nybin,Float_t psi)
209{
210 //Create histograms.
211 //_Only_ to be used in case of the adaptive histograms!
212 //phimax is given in radians!!
213
214 if(ptmin > ptmax)
215 {
4aa41877 216 cerr<<"AliHLTHoughTransformerLUT::CreateHistograms: Error in ptrange "<<ptmin<<" "<<ptmax<<endl;
1f1942b8 217 return;
218 }
219 if(psi < 0)
220 {
4aa41877 221 cerr<<"AliHLTHoughTransformerLUT::CreateHistograms: Wrong psi-angle "<<psi<<endl;
1f1942b8 222 return;
223 }
224
4aa41877 225 fParamSpace = new AliHLTHistogram*[fNEtas];
1f1942b8 226 Char_t histname[256];
227 for(Int_t i=0; i<GetNEtaSegments(); i++)
228 {
229 sprintf(histname,"paramspace_%d",i);
4aa41877 230 fParamSpace[i] = new AliHLTHistogramAdaptive(histname,ptmin,ptmax,ptres,nybin,-psi,psi);
1f1942b8 231 }
232
233 //create lookup table for sin and cos
234 fNPhi0=nybin;
235 fLUTphi0=new Float_t[fNPhi0];
236 fLUT2sinphi0=new Float_t[fNPhi0];
237 fLUT2cosphi0=new Float_t[fNPhi0];
238 fLUTKappa=new Float_t[fNPhi0];
4aa41877 239 AliHLTHistogram *hist=fParamSpace[0];
1f1942b8 240 Int_t i=0;
241 for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
242 {
243 Float_t phi0 = hist->GetBinCenterY(b);
244 fLUTphi0[i]=phi0;
245 fLUT2sinphi0[i]=2.*sin(phi0);
246 fLUT2cosphi0[i]=2.*cos(phi0);
247 fLUTKappa[i]=0.;
248 i++;
249 //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
250 }
251}
252
4aa41877 253void AliHLTHoughTransformerLUT::CreateHistograms(Int_t nxbin,Float_t ptMin,Int_t nybin,Float_t phimin,Float_t phimax)
d5d754e5 254{
255 //Create the histograms (parameter space).
256 //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
257 //The arguments give the range and binning;
258 //nxbin = #bins in kappa
259 //nybin = #bins in phi0
f5271e68 260 //ptMin = mimium Pt of track (corresponding to maximum kappa)
d5d754e5 261 //phi_min = mimimum phi0 (degrees)
262 //phi_max = maximum phi0 (degrees)
263
4aa41877 264 Double_t x = AliHLTTransform::GetBFact()*AliHLTTransform::GetBField()/ptMin;
265 //Double_t torad = AliHLTTransform::Pi()/180;
1f1942b8 266 CreateHistograms(nxbin,-1.*x,x,nybin,phimin/**torad*/,phimax/**torad*/);
d5d754e5 267}
268
4aa41877 269void AliHLTHoughTransformerLUT::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,Int_t nybin,Float_t ymin,Float_t ymax)
d5d754e5 270{
4aa41877 271 fParamSpace = new AliHLTHistogram*[fNEtas];
d5d754e5 272
273 Char_t histname[256];
274 for(Int_t i=0; i<fNEtas; i++)
275 {
276 sprintf(histname,"paramspace_%d",i);
4aa41877 277 fParamSpace[i] = new AliHLTHistogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
d5d754e5 278 }
279
d5d754e5 280 //create lookup table for sin and cos
1f1942b8 281 fNPhi0=nybin;
d5d754e5 282 fLUTphi0=new Float_t[fNPhi0];
283 fLUT2sinphi0=new Float_t[fNPhi0];
284 fLUT2cosphi0=new Float_t[fNPhi0];
18659a6b 285 fLUTKappa=new Float_t[fNPhi0];
1f1942b8 286
4aa41877 287 AliHLTHistogram *hist=fParamSpace[0];
1f1942b8 288 Int_t i=0;
289 for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
290 {
291 Float_t phi0 = hist->GetBinCenterY(b);
292 fLUTphi0[i]=phi0;
293 fLUT2sinphi0[i]=2.*sin(phi0);
294 fLUT2cosphi0[i]=2.*cos(phi0);
295 fLUTKappa[i]=0.;
296 //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
297 i++;
298 }
d5d754e5 299}
300
4aa41877 301void AliHLTHoughTransformerLUT::Reset()
d5d754e5 302{
303 //Reset all the histograms
304
305 if(!fParamSpace)
306 {
4aa41877 307 LOG(AliHLTLog::kWarning,"AliHLTHoughTransformer::Reset","Histograms")
d5d754e5 308 <<"No histograms to reset"<<ENDLOG;
309 return;
310 }
311
3bcf971b 312 for(Int_t i=0; i<fNEtas; i++)
d5d754e5 313 fParamSpace[i]->Reset();
d5d754e5 314}
315
4aa41877 316Int_t AliHLTHoughTransformerLUT::GetEtaIndex(Double_t eta) const
d5d754e5 317{
318 //Return the histogram index of the corresponding eta.
3bcf971b 319
3e87ef69 320#ifdef FULLLUT
321 /* try to imitate a circuit -> should
322 go into the VHDL implementation of transformer */
3bcf971b 323 Float_t rz2=CalcRoverZ2(eta);
324 return FindIndex(rz2);
3e87ef69 325#else /* optimize for speed on the computer */
326 Double_t index = (eta-GetEtaMin())/fEtaSlice;
327 return (Int_t)index;
328#endif
3bcf971b 329}
330
4aa41877 331AliHLTHistogram *AliHLTHoughTransformerLUT::GetHistogram(Int_t etaIndex)
d5d754e5 332{
f5271e68 333 // gets hitogram
334 if(!fParamSpace || etaIndex >= fNEtas || etaIndex < 0)
d5d754e5 335 return 0;
f5271e68 336 if(!fParamSpace[etaIndex])
d5d754e5 337 return 0;
f5271e68 338 return fParamSpace[etaIndex];
d5d754e5 339}
340
4aa41877 341Double_t AliHLTHoughTransformerLUT::GetEta(Int_t etaIndex,Int_t slice) const
d5d754e5 342{
f5271e68 343 // gets eta
344 if(etaIndex >= fNEtas || etaIndex < 0){
4aa41877 345 LOG(AliHLTLog::kWarning,"AliHLTHoughTransformerLUT::GetEta","Index") << "Index out of range."<<ENDLOG;
3bcf971b 346 return 0.;
347 }
348 if(slice != fSlice){
4aa41877 349 LOG(AliHLTLog::kWarning,"AliHLTHoughTransformerLUT::GetEta","Index") << "Given slice does not match internal slice."<<ENDLOG;
3bcf971b 350 return 0.;
351 }
352
f5271e68 353 return(fLUTEtaReal[etaIndex]);
3e87ef69 354}
355
4aa41877 356void AliHLTHoughTransformerLUT::TransformCircle()
d5d754e5 357{
d5d754e5 358 //Transform the input data with a circle HT.
3bcf971b 359 //The function loops over all the data, and transforms each pixel with the equation:
d5d754e5 360 //
3bcf971b 361 //kappa = 2/R^2 * (y*cos(phi0) - x*sin(phi0) )
d5d754e5 362 //
3bcf971b 363 //where R^2 = x^2 + y^2
d5d754e5 364 //
1f1942b8 365 //After a day of testing it is proven that this h
d5d754e5 366
4aa41877 367 AliHLTDigitRowData *tempPt = GetDataPointer();
d5d754e5 368 if(!tempPt)
369 {
4aa41877 370 LOG(AliHLTLog::kError,"AliHLTHoughTransformerLUT::TransformCircle","Data")
d5d754e5 371 <<"No input data "<<ENDLOG;
372 return;
373 }
1f1942b8 374
375 Int_t lowch=GetLowerThreshold();
376 Int_t highch=GetUpperThreshold();
377
d5d754e5 378 //Loop over the padrows:
3bcf971b 379 for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
d5d754e5 380 {
381 //Get the data on this padrow:
4aa41877 382 AliHLTDigitData *digPt = tempPt->fDigitData;
d5d754e5 383 if(i != (Int_t)tempPt->fRow)
384 {
4aa41877 385 LOG(AliHLTLog::kError,"AliHLTHoughTransformerLUT::TransformCircle","Data")
386 <<"AliHLTHoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
d5d754e5 387 continue;
388 }
389
efb8e06f 390
391 //start a new row
18659a6b 392 fLastPad=-1;
1f1942b8 393 //store x for this row
394 fX = CalcX(row);
395 //accumulate charge per etaslice
3e87ef69 396 fAccCharge=0;
1f1942b8 397 //last histogram
398 fLastIndex=-1;
3e87ef69 399
d5d754e5 400 //Loop over the data on this padrow:
401 for(UInt_t j=0; j<tempPt->fNDigit; j++)
402 {
403 UShort_t charge = digPt[j].fCharge;
3bcf971b 404
405 //check threshold
1f1942b8 406 if((Int_t)charge <= lowch)
d5d754e5 407 continue;
1f1942b8 408 if ((Int_t)charge > highch)
409 charge=highch;
3bcf971b 410
411 UChar_t pad = digPt[j].fPad;
412 UShort_t time = digPt[j].fTime;
413
18659a6b 414 if(fLastPad!=pad){ //only update if necessary
3e87ef69 415
1f1942b8 416 //if there is accumalated charge,
417 //put it in the old histogram
418 //using the old X and Y values
419 if(fAccCharge>0) AddCurveToHistogram(-1); //fLastIndex will be set to -1
3e87ef69 420
421#ifdef FULLLUT
422 fLastIndex=fNEtas-1;
423#endif
424
1f1942b8 425 //calculate Y for the new pad
426 fY = CalcY(pad,row);
427 //remember new pad
18659a6b 428 fLastPad=pad;
429 }
430
1f1942b8 431#ifdef FULLLUT
3bcf971b 432 //find eta slice
efb8e06f 433 Float_t z = CalcZ(time);
3e87ef69 434
1f1942b8 435 Float_t z2=z*z;
efb8e06f 436 Float_t rz2 = 1 + r2/z2;
f5271e68 437 Int_t etaIndex = FindIndex(rz2,fLastIndex);
3e87ef69 438#else
1f1942b8 439 Float_t z = CalcZ(time);
440 Double_t r = sqrt(fX*fX+fY*fY+z*z);
441 Double_t eta = 0.5 * log((r+z)/(r-z));
f5271e68 442 Int_t etaIndex = GetEtaIndex(eta);
3e87ef69 443#endif
f5271e68 444 if(etaIndex < 0 || etaIndex >= fNEtas){
4aa41877 445 LOG(AliHLTLog::kWarning,"AliHLTHoughTransformerLUT::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<etaIndex<<"."<<ENDLOG;
d5d754e5 446 continue;
3bcf971b 447 }
3e87ef69 448
449#ifndef FULLLUT
f5271e68 450 if(fLastIndex==-1) fLastIndex=etaIndex;
d5d754e5 451#endif
3e87ef69 452
f5271e68 453 if(fLastIndex!=etaIndex){ //enter old values first
454 AddCurveToHistogram(etaIndex);
3bcf971b 455 }
3e87ef69 456
457 fAccCharge+=charge;
458 //fAccCharge+=1;
d5d754e5 459 }
3e87ef69 460
461 //in case there is charge, store it before moving to the next row
462 if(fAccCharge>0) AddCurveToHistogram(-1);
efb8e06f 463
d5d754e5 464 //Move the data pointer to the next padrow:
4aa41877 465 AliHLTMemHandler::UpdateRowPointer(tempPt);
d5d754e5 466 }
d5d754e5 467}
468
4aa41877 469void AliHLTHoughTransformerLUT::Print()
d5d754e5 470{
f5271e68 471 // debug printout
d5d754e5 472 cout << "fSlice: " << GetSlice() << endl;
473 cout << "fPatch: " << GetPatch() << endl;
474 cout << "fSector: " << fSector << endl;
475 cout << "fSectorRow: " << fSectorRow << endl;
476 cout << "fMinRow: " << fMinRow << endl;
477 cout << "fMaxRow: " << fMaxRow << endl;
478 cout << "fNRows: " << fNRows << endl;
479 cout << "fNEtas: " << fNEtas << endl;
480 cout << "fNPhi0: " << fNPhi0 << endl;
481 cout << "fZSign: " << fZSign << endl;
482 cout << "fZLengthPlusOff: " << fZLengthPlusOff << endl;
483 cout << "fPadPitch: " << fPadPitch << endl;
484 cout << "fTimeWidth: " << fTimeWidth << endl;
485
486 if(!fNRows) return;
487 cout << "fLUTX " << fNRows << endl;
488 for(Int_t i=0;i<fNRows;i++) cout << "fLUTX[" << i << "]=" << (Float_t)fLUTX[i] << endl;
489 cout << "fLUTY " << fNRows << endl;
490 for(Int_t i=0;i<fNRows;i++) cout << "fLUTY[" << i << "]=" << fLUTY[i] << endl;
491 if(!fNEtas) return;
492 cout << "fLUTEta " << fNEtas << endl;
493 for(Int_t i=0;i<fNEtas;i++) cout << "fLUTEta[" << i << "]=" << fLUTEta[i] << endl;
3e87ef69 494 cout << "fLUTEtaReal " << fNEtas << endl;
495 for(Int_t i=0;i<fNEtas;i++) cout << "fLUTEtaReal[" << i << "]=" << fLUTEtaReal[i] << endl;
d5d754e5 496 if(!fNPhi0) return;
497 cout << "fLUTphi0 " << fNPhi0 << endl;
498 for(Int_t i=0;i<fNPhi0;i++) cout << "fLUTPhi0[" << i << "]=" << fLUTphi0[i] << endl;
499 cout << "fLUT2sinphi0 " << fNPhi0 << endl;
500 for(Int_t i=0;i<fNPhi0;i++) cout << "fLUT2sinphi0[" << i << "]=" << fLUT2sinphi0[i] << endl;
501 cout << "fLUT2cosphi0 " << fNPhi0 << endl;
502 for(Int_t i=0;i<fNPhi0;i++) cout << "fLUT2cosphi0[" << i << "]=" << fLUT2cosphi0[i] << endl;
503}