3 // Author: Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
4 //*-- Copyright © ALICE HLT Group
5 /** /class AliL3HoughTransformerLUT
7 //_____________________________________________________________
8 // AliL3HoughTransformerLUT
10 // Hough transformation class using Look-UP-Tables
15 #include "AliL3StandardIncludes.h"
17 #include "AliL3Logging.h"
18 #include "AliL3HoughTransformerLUT.h"
19 #include "AliL3Transform.h"
20 #include "AliL3MemHandler.h"
21 #include "AliL3DigitData.h"
22 #include "AliL3Histogram.h"
23 #include "AliL3HistogramAdaptive.h"
29 ClassImp(AliL3HoughTransformerLUT)
31 AliL3HoughTransformerLUT::AliL3HoughTransformerLUT() : AliL3HoughBaseTransformer()
65 AliL3HoughTransformerLUT::AliL3HoughTransformerLUT(Int_t slice,Int_t patch,Int_t nEtaSegments) : AliL3HoughBaseTransformer(slice,patch,nEtaSegments)
97 Init(slice,patch,nEtaSegments);
100 AliL3HoughTransformerLUT::~AliL3HoughTransformerLUT()
113 delete[] fLUTEtaReal;
118 void AliL3HoughTransformerLUT::DeleteHistograms()
120 // deletes all histograms
122 delete[] fLUT2sinphi0;
123 delete[] fLUT2cosphi0;
134 for(Int_t i=0; i<fNEtas; i++)
136 if(!fParamSpace[i]) continue;
137 delete fParamSpace[i];
139 delete [] fParamSpace;
143 void AliL3HoughTransformerLUT::Init(Int_t slice,Int_t patch,Int_t nEtaSegments,Int_t /*nSeqs*/)
146 AliL3HoughBaseTransformer::Init(slice,patch,nEtaSegments);
148 //delete old LUT tables
156 delete[] fLUTEtaReal;
161 fMinRow=AliL3Transform::GetFirstRow(patch);;
162 fMaxRow=AliL3Transform::GetLastRow(patch);;
163 fNRows=AliL3Transform::GetNRows(patch);;
165 if(fNEtas!=GetNEtaSegments()) {
166 cerr << "AliL3HoughTransformerLUT::Init -> Error: Number of Etas must be equal!" << endl;
170 AliL3Transform::Slice2Sector(slice,fMinRow,fSector,fSectorRow);
171 fZSign = slice < 18 ? 1:-1; //see AliL3Transformer
173 fZLengthPlusOff=AliL3Transform::GetZLength()+AliL3Transform::GetZOffset();
174 fTimeWidth=AliL3Transform::GetZWidth();
177 if(fSector<AliL3Transform::GetNSectorLow())
178 fPadPitch=AliL3Transform::GetPadPitchWidthLow();
180 fPadPitch=AliL3Transform::GetPadPitchWidthUp();
182 Float_t etaMax=GetEtaMax();
183 Float_t etaMin=GetEtaMin();
184 fEtaSlice=(etaMax-etaMin)/nEtaSegments;
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++){
190 fLUTX[rr]=Float_t(AliL3Transform::Row2X(rr+fMinRow));
191 fLUTY[rr]=Float_t(0.5*(AliL3Transform::GetNPads(rr+fMinRow)-1)*fPadPitch);
192 //cout << rr << ": " << (Float_t)fLUTX[rr] << " " << (Float_t)fLUTY[rr] << endl;
195 //lookup tables for rz2s <=> etas
196 //will only be used ifdef FULLLUT
197 fLUTEta=new Float_t[fNEtas];
198 fLUTEtaReal=new Float_t[fNEtas];
199 for(Int_t rr=0;rr<fNEtas;rr++){
200 Float_t eta=etaMin+(rr+1)*fEtaSlice;
201 fLUTEta[rr]=CalcRoverZ2(eta);
202 fLUTEtaReal[rr]=eta-0.5*fEtaSlice;
203 //cout << rr << ": " << eta << " " << fLUTEtaReal[rr] << " " << GetEta(rr,fSlice) << " - " << fLUTEta[rr] << endl;
207 void AliL3HoughTransformerLUT::CreateHistograms(Float_t ptmin,Float_t ptmax,Float_t ptres,
208 Int_t nybin,Float_t psi)
211 //_Only_ to be used in case of the adaptive histograms!
212 //phimax is given in radians!!
216 cerr<<"AliL3HoughTransformerLUT::CreateHistograms: Error in ptrange "<<ptmin<<" "<<ptmax<<endl;
221 cerr<<"AliL3HoughTransformerLUT::CreateHistograms: Wrong psi-angle "<<psi<<endl;
225 fParamSpace = new AliL3Histogram*[fNEtas];
226 Char_t histname[256];
227 for(Int_t i=0; i<GetNEtaSegments(); i++)
229 sprintf(histname,"paramspace_%d",i);
230 fParamSpace[i] = new AliL3HistogramAdaptive(histname,ptmin,ptmax,ptres,nybin,-psi,psi);
233 //create lookup table for sin and cos
235 fLUTphi0=new Float_t[fNPhi0];
236 fLUT2sinphi0=new Float_t[fNPhi0];
237 fLUT2cosphi0=new Float_t[fNPhi0];
238 fLUTKappa=new Float_t[fNPhi0];
239 AliL3Histogram *hist=fParamSpace[0];
241 for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
243 Float_t phi0 = hist->GetBinCenterY(b);
245 fLUT2sinphi0[i]=2.*sin(phi0);
246 fLUT2cosphi0[i]=2.*cos(phi0);
249 //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
253 void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Float_t ptMin,Int_t nybin,Float_t phimin,Float_t phimax)
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
260 //ptMin = mimium Pt of track (corresponding to maximum kappa)
261 //phi_min = mimimum phi0 (degrees)
262 //phi_max = maximum phi0 (degrees)
264 Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/ptMin;
265 //Double_t torad = AliL3Transform::Pi()/180;
266 CreateHistograms(nxbin,-1.*x,x,nybin,phimin/**torad*/,phimax/**torad*/);
269 void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,Int_t nybin,Float_t ymin,Float_t ymax)
271 fParamSpace = new AliL3Histogram*[fNEtas];
273 Char_t histname[256];
274 for(Int_t i=0; i<fNEtas; i++)
276 sprintf(histname,"paramspace_%d",i);
277 fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
280 //create lookup table for sin and cos
282 fLUTphi0=new Float_t[fNPhi0];
283 fLUT2sinphi0=new Float_t[fNPhi0];
284 fLUT2cosphi0=new Float_t[fNPhi0];
285 fLUTKappa=new Float_t[fNPhi0];
287 AliL3Histogram *hist=fParamSpace[0];
289 for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
291 Float_t phi0 = hist->GetBinCenterY(b);
293 fLUT2sinphi0[i]=2.*sin(phi0);
294 fLUT2cosphi0[i]=2.*cos(phi0);
296 //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
301 void AliL3HoughTransformerLUT::Reset()
303 //Reset all the histograms
307 LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
308 <<"No histograms to reset"<<ENDLOG;
312 for(Int_t i=0; i<fNEtas; i++)
313 fParamSpace[i]->Reset();
316 Int_t AliL3HoughTransformerLUT::GetEtaIndex(Double_t eta)
318 //Return the histogram index of the corresponding eta.
321 /* try to imitate a circuit -> should
322 go into the VHDL implementation of transformer */
323 Float_t rz2=CalcRoverZ2(eta);
324 return FindIndex(rz2);
325 #else /* optimize for speed on the computer */
326 Double_t index = (eta-GetEtaMin())/fEtaSlice;
331 AliL3Histogram *AliL3HoughTransformerLUT::GetHistogram(Int_t etaIndex)
334 if(!fParamSpace || etaIndex >= fNEtas || etaIndex < 0)
336 if(!fParamSpace[etaIndex])
338 return fParamSpace[etaIndex];
341 Double_t AliL3HoughTransformerLUT::GetEta(Int_t etaIndex,Int_t slice) const
344 if(etaIndex >= fNEtas || etaIndex < 0){
345 LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::GetEta","Index") << "Index out of range."<<ENDLOG;
349 LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::GetEta","Index") << "Given slice does not match internal slice."<<ENDLOG;
353 return(fLUTEtaReal[etaIndex]);
356 void AliL3HoughTransformerLUT::TransformCircle()
358 //Transform the input data with a circle HT.
359 //The function loops over all the data, and transforms each pixel with the equation:
361 //kappa = 2/R^2 * (y*cos(phi0) - x*sin(phi0) )
363 //where R^2 = x^2 + y^2
365 //After a day of testing it is proven that this h
367 AliL3DigitRowData *tempPt = GetDataPointer();
370 LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
371 <<"No input data "<<ENDLOG;
375 Int_t lowch=GetLowerThreshold();
376 Int_t highch=GetUpperThreshold();
378 //Loop over the padrows:
379 for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
381 //Get the data on this padrow:
382 AliL3DigitData *digPt = tempPt->fDigitData;
383 if(i != (Int_t)tempPt->fRow)
385 LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
386 <<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
393 //store x for this row
395 //accumulate charge per etaslice
400 //Loop over the data on this padrow:
401 for(UInt_t j=0; j<tempPt->fNDigit; j++)
403 UShort_t charge = digPt[j].fCharge;
406 if((Int_t)charge <= lowch)
408 if ((Int_t)charge > highch)
411 UChar_t pad = digPt[j].fPad;
412 UShort_t time = digPt[j].fTime;
414 if(fLastPad!=pad){ //only update if necessary
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
425 //calculate Y for the new pad
433 Float_t z = CalcZ(time);
436 Float_t rz2 = 1 + r2/z2;
437 Int_t etaIndex = FindIndex(rz2,fLastIndex);
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));
442 Int_t etaIndex = GetEtaIndex(eta);
444 if(etaIndex < 0 || etaIndex >= fNEtas){
445 LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<etaIndex<<"."<<ENDLOG;
450 if(fLastIndex==-1) fLastIndex=etaIndex;
453 if(fLastIndex!=etaIndex){ //enter old values first
454 AddCurveToHistogram(etaIndex);
461 //in case there is charge, store it before moving to the next row
462 if(fAccCharge>0) AddCurveToHistogram(-1);
464 //Move the data pointer to the next padrow:
465 AliL3MemHandler::UpdateRowPointer(tempPt);
469 void AliL3HoughTransformerLUT::Print()
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;
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;
492 cout << "fLUTEta " << fNEtas << endl;
493 for(Int_t i=0;i<fNEtas;i++) cout << "fLUTEta[" << i << "]=" << fLUTEta[i] << endl;
494 cout << "fLUTEtaReal " << fNEtas << endl;
495 for(Int_t i=0;i<fNEtas;i++) cout << "fLUTEtaReal[" << i << "]=" << fLUTEtaReal[i] << endl;
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;