//$Id$
// Author: Constantin Loizides <mailto:loizides@fi.uib.no>
-//*-- Copyright © CL and ASV
+//*-- Copyright&Copy CL
-#include "AliL3MemHandler.h"
+#include "AliL3StandardIncludes.h"
+
+#include "AliL3RootTypes.h"
#include "AliL3Logging.h"
+#include "AliL3MemHandler.h"
#include "AliL3Transform.h"
#include "AliL3DigitData.h"
-#include "AliL3Histogram.h"
#include "AliL3HoughTransformerVhdl.h"
+#include "AliL3FFloat.h"
+
+#if GCCVERSION == 3
+using namespace std;
+#endif
+/** \class AliL3HoughTransformerVhdl
+// <pre>
//_____________________________________________________________
// AliL3HoughTransformerVhdl
//
-// Hough transformation class
+// Hough transformation class for VHDL comparism.
//
+//</pre>
+*/
ClassImp(AliL3HoughTransformerVhdl)
+
AliL3HoughTransformerVhdl::AliL3HoughTransformerVhdl()
+ : AliL3HoughTransformerLUT()
+
{
- //Default constructor
- fParamSpace = 0;
+ fEpsilon=0;
+ fSinEpsilon=0;
+ fCosEpsilon=1;
+ fIts=0;
}
-AliL3HoughTransformerVhdl::AliL3HoughTransformerVhdl(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
+AliL3HoughTransformerVhdl::AliL3HoughTransformerVhdl(Int_t slice,Int_t patch,Int_t n_eta_segments,Int_t n_its)
+ : AliL3HoughTransformerLUT(slice,patch,n_eta_segments)
{
- //Normal constructor
- fParamSpace = 0;
+ fEpsilon=0;
+ fSinEpsilon=0;
+ fCosEpsilon=1;
+ fIts=n_its;
}
AliL3HoughTransformerVhdl::~AliL3HoughTransformerVhdl()
{
- DeleteHistograms();
-}
-void AliL3HoughTransformerVhdl::DeleteHistograms()
-{
- if(!fParamSpace)
- return;
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- if(!fParamSpace[i]) continue;
- delete fParamSpace[i];
- }
- delete [] fParamSpace;
}
-void AliL3HoughTransformerVhdl::CreateHistograms(Int_t nxbin,Double_t pt_min,
- Int_t nybin,Double_t phimin,Double_t phimax)
+void AliL3HoughTransformerVhdl::CreateHistograms(Int_t nxbin,Double_t pt_min,Int_t nybin,Double_t phimin,Double_t phimax)
{
- //Create the histograms (parameter space).
- //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
- //The arguments give the range and binning;
- //nxbin = #bins in kappa
- //nybin = #bins in phi0
- //pt_min = mimium Pt of track (corresponding to maximum kappa)
- //phi_min = mimimum phi0 (degrees)
- //phi_max = maximum phi0 (degrees)
-
- Double_t bfact = 0.0029980;
- Double_t bfield = AliL3Transform::GetBField();
- Double_t x = bfact*bfield/pt_min;
- Double_t torad = AliL3Transform::Pi()/180;
- CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
+ AliL3HoughTransformerLUT::CreateHistograms(nxbin,pt_min,nybin,phimin,phimax);
}
void AliL3HoughTransformerVhdl::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
- Int_t nybin,Double_t ymin,Double_t ymax)
+ Int_t nybin,Double_t ymin,Double_t ymax)
{
-
- fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
-
- Char_t histname[256];
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- {
- sprintf(histname,"paramspace_%d",i);
- fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
- }
-}
+ AliL3HoughTransformerLUT::CreateHistograms(nxbin,xmin,xmax,nybin,ymin,ymax);
-void AliL3HoughTransformerVhdl::Reset()
-{
- //Reset all the histograms
+ fEpsilon=(ymax-ymin)/nybin;
+ fSinEpsilon=sin(fEpsilon);
+ fCosEpsilon=cos(fEpsilon);
- if(!fParamSpace)
- {
- LOG(AliL3Log::kWarning,"AliL3HoughTransformerVhdl::Reset","Histograms")
- <<"No histograms to reset"<<ENDLOG;
- return;
- }
-
- for(Int_t i=0; i<GetNEtaSegments(); i++)
- fParamSpace[i]->Reset();
+ fNxbin=nxbin;
+ fXmin=xmin;
+ fXmax=xmax;
+ fNybin=nybin;
+ fYmin=ymin;
+ fYmax=ymax;
+ // cout << fEpsilon << " - " << (xmax-xmin)/nxbin << endl;
}
-
-Int_t AliL3HoughTransformerVhdl::GetEtaIndex(Double_t eta)
+void AliL3HoughTransformerVhdl::Init(Int_t slice,Int_t patch,Int_t n_eta_segments,Int_t n_its)
{
- //Return the histogram index of the corresponding eta.
-
- Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
- Double_t index = (eta-GetEtaMin())/etaslice;
- return (Int_t)index;
+ AliL3HoughTransformerLUT::Init(slice,patch,n_eta_segments);
}
void AliL3HoughTransformerVhdl::TransformCircle()
{
//Transform the input data with a circle HT.
- //The function loops over all the data, and transforms each pixel with the equations:
+ //The function loops over all the data, and transforms each pixel with the equation:
//
- //kappa = 2/R*sin(phi - phi0)
+ //kappa = lastkappa +- epsilon * lastkappaprime
//
- //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
+ //kappaprime = lastkappaprime -+ epsilon * lastkappa
//
//Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
//which histogram in which the pixel should be transformed, the eta-value is calcluated
<<"No input data "<<ENDLOG;
return;
}
+
+//#define use_error
+#ifdef use_error
+ Float_t max_error=0;
+ Float_t rel_error=0;
+ Int_t counter=0;
+#endif
//Loop over the padrows:
- for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
+ for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
{
//Get the data on this padrow:
AliL3DigitData *digPt = tempPt->fDigitData;
if(i != (Int_t)tempPt->fRow)
{
- printf("AliL3HoughTransform::TransformCircle : Mismatching padrow numbering\n");
+ LOG(AliL3Log::kError,"AliL3HoughTransformerVhdl::TransformCircle","Data")
+ <<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
continue;
}
-
+
+ Float_t x = CalcX(row);
+ Float_t x2=x*x;
+ Float_t y=0,y2=0;
+ Float_t r2=0;
+ Float_t R2=0;
+ fLastPad=-1;
+
+
//Loop over the data on this padrow:
for(UInt_t j=0; j<tempPt->fNDigit; j++)
{
UShort_t charge = digPt[j].fCharge;
+
+ //check threshold
+ if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
+ continue;
+
UChar_t pad = digPt[j].fPad;
UShort_t time = digPt[j].fTime;
- if((Int_t)charge <= GetLowerThreshold())
- continue;
- Int_t sector,row;
- Float_t xyz[3];
-
- //Transform data to local cartesian coordinates:
- AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
- AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
-
- //Calculate the eta:
- Double_t eta = AliL3Transform::GetEta(xyz);
-
- //Get the corresponding index, which determines which histogram to fill:
- Int_t eta_index = GetEtaIndex(eta);
- if(eta_index < 0 || eta_index >= GetNEtaSegments())
+
+ if(fLastPad!=pad){ //only update if necessary
+ fLastIndex=fNEtas-1;
+ y = CalcY(pad,row);
+ y2 = y*y;
+ r2 = x2 + y2;
+ R2 = 1. / r2;
+ for(Int_t b=0; b<fNPhi0; b++)
+ fLUTKappa[b]=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
+
+ //Fill the histogram along the phirange
+ Float_t kappa=0;
+ Float_t kappaprime=0;
+ Float_t lastkappa=0;
+ Float_t lastkappaprime=0;
+ Int_t its=0;
+ Float_t R2=1/(x*x+y*y);
+ Float_t A=2*y*R2;
+ Float_t B=-2*x*R2;
+
+ //starting point
+ Float_t phi = fLUTphi0[fNPhi0/2];
+ kappa=A*cos(phi)+B*sin(phi);
+ kappaprime=B*cos(phi)-A*sin(phi);
+ its=fIts;
+ lastkappa=kappa;
+ lastkappaprime=kappaprime;
+ // hist->Fill(kappa,phi,charge);
+
+ for(Int_t b=fNPhi0/2+1; b<fNPhi0; b++){
+ Float_t exact_kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
+
+ phi=fLUTphi0[b];
+
+ if(its==0) { //initialize or re-adjust values
+ kappa=A*cos(phi)+B*sin(phi); //equals exact_kappa!
+ kappaprime=B*cos(phi)-A*sin(phi);
+ its=fIts;
+ } else {
+ //kappa=fCosEpsilon*lastkappa+fSinEpsilon*lastkappaprime;
+ //kappaprime=fCosEpsilon*lastkappaprime-fSinEpsilon*lastkappa;
+ kappa=lastkappa+fEpsilon*lastkappaprime;
+ kappaprime=lastkappaprime-fEpsilon*lastkappa;
+ its--;
+ }
+
+ lastkappa=kappa;
+ lastkappaprime=kappaprime;
+
+ // hist->Fill(kappa,phi,charge);
+ } // end positive running values
+
+ phi = fLUTphi0[fNPhi0/2];
+ kappa=A*cos(phi)+B*sin(phi);
+ kappaprime=B*cos(phi)-A*sin(phi);
+ its=fIts;
+ lastkappa=kappa;
+ lastkappaprime=kappaprime;
+ //hist->Fill(kappa,fLUTphi0[b],charge);
+
+ for(Int_t b=fNPhi0/2-1; b>=0; b--){
+ Float_t exact_kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
+
+ Float_t phi = fLUTphi0[b];
+
+ if(its==0) { //initialize or re-adjust values
+ kappa=A*cos(phi)+B*sin(phi); //equals exact_kappa!
+ kappaprime=B*cos(phi)-A*sin(phi);
+ its=fIts;
+ } else {
+ //kappa=fCosEpsilon*lastkappa-fSinEpsilon*lastkappaprime;
+ //kappaprime=fCosEpsilon*lastkappaprime+fSinEpsilon*lastkappa;
+ kappa=lastkappa-fEpsilon*lastkappaprime;
+ kappaprime=lastkappaprime+fEpsilon*lastkappa;
+ its--;
+ }
+
+ lastkappa=kappa;
+ lastkappaprime=kappaprime;
+
+ // hist->Fill(kappa,phi,charge);
+ }
+
+ fLastPad=pad;
+ }
+
+ Float_t z = CalcZ(time);
+
+ //find eta slice
+ Float_t rz2 = 1 + (x*x+y*y)/(z*z);
+ Int_t eta_index = FindIndex(rz2);
+
+ if(eta_index < 0 || eta_index >= fNEtas){
+ //LOG(AliL3Log::kWarning,"AliL3HoughTransformerVhdl::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<eta_index<<"."<<ENDLOG;
continue;
-
+ }
+
//Get the correct histogrampointer:
AliL3Histogram *hist = fParamSpace[eta_index];
- if(!hist)
- {
- printf("AliL3HoughTransformerVhdl::TransformCircle : Error getting histogram in index %d\n",eta_index);
- continue;
- }
+ if(!hist){
+ //LOG(AliL3Log::kWarning,"AliL3HoughTransformerVhdl::TransformCircle","Histograms")<<"Error getting histogram in index "<<eta_index<<"."<<ENDLOG;
+ continue;
+ }
- //Do the transformation:
- Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
- Float_t phi = AliL3Transform::GetPhi(xyz);
-
- //Fill the histogram along the phirange
- for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
- {
- Float_t phi0 = hist->GetBinCenterY(b);
- Float_t kappa = 2*sin(phi - phi0)/R;
- hist->Fill(kappa,phi0,charge);
+ for(Int_t b=0; b<fNPhi0; b++){
+
+#ifdef use_error
+ Float_t err=fabs((exact_kappa-kappa)/exact_kappa);
+
+ if(err>max_error) {
+ //cout << max_error << " - " << err << " " << kappa << " " << exact_kappa << " " << kappa/exact_kappa<< endl;
+ max_error=err;
}
+
+ rel_error+=err;
+ counter++;
+#endif
+ //hist->Fill(kappa,fLUTphi0[b],charge);
+
+ //cout << kappa << " " << fLUTphi0[b] << " " << charge << endl;
+ }
}
-
//Move the data pointer to the next padrow:
AliL3MemHandler::UpdateRowPointer(tempPt);
}
+
+#ifdef use_error
+ cout <<"Max Error: " << max_error << endl;
+ cout <<"Rel Error average: " << rel_error/counter << endl;
+#endif
+
}
-void AliL3HoughTransformerVhdl::TransformCircleC(Int_t row_range)
+void AliL3HoughTransformerVhdl::Print()
{
- //Circle transform, using combinations of every 2 points lying
- //on different padrows and within the same etaslice.
-
- AliL3DigitRowData *tempPt = GetDataPointer();
- if(!tempPt)
- LOG(AliL3Log::kError,"AliL3HoughTransformerVhdl::TransformCircleC","Data")
- <<"No input data "<<ENDLOG;
-
- Int_t counter=0;
- for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
- {
- counter += tempPt->fNDigit;
- AliL3MemHandler::UpdateRowPointer(tempPt);
- }
-
- struct Digit {
- Int_t row;
- Double_t r;
- Double_t phi;
- Int_t eta_index;
- Int_t charge;
- };
-
- Digit *digits = new Digit[counter];
- cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
-
- Int_t total_digits=counter;
- Int_t sector,row,tot_charge,pad,time,charge;
- Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
- Float_t xyz[3];
-
- counter=0;
- tempPt = GetDataPointer();
-
- for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
- {
- AliL3DigitData *digPt = tempPt->fDigitData;
- for(UInt_t di=0; di<tempPt->fNDigit; di++)
- {
- charge = digPt[di].fCharge;
- pad = digPt[di].fPad;
- time = digPt[di].fTime;
- AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
- AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
- eta = AliL3Transform::GetEta(xyz);
- digits[counter].row = i;
- digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
- digits[counter].phi = atan2(xyz[1],xyz[0]);
- digits[counter].eta_index = GetEtaIndex(eta);
- digits[counter].charge = charge;
- counter++;
- }
- AliL3MemHandler::UpdateRowPointer(tempPt);
- }
-
- for(Int_t i=0; i<total_digits; i++)
- {
- if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue;
- Int_t ind = digits[i].eta_index;
-
- for(Int_t j=i+1; j<total_digits; j++)
- {
- if(digits[i].row == digits[j].row) continue;
- if(digits[i].eta_index != digits[j].eta_index) continue;
- if(digits[i].row + row_range < digits[j].row) break;
-
- //Get the correct histogrampointer:
- AliL3Histogram *hist = fParamSpace[ind];
- if(!hist)
- {
- printf("AliL3HoughTransformerVhdl::TransformCircleC() : No histogram at index %d\n",ind);
- continue;
- }
-
- r1 = digits[i].r;
- phi1 = digits[i].phi;
- r2 = digits[j].r;
- phi2 = digits[j].phi;
- phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
- kappa = 2*sin(phi2-phi_0)/r2;
- tot_charge = digits[i].charge + digits[j].charge;
- hist->Fill(kappa,phi_0,tot_charge);
- }
- }
- delete [] digits;
+ AliL3HoughTransformerLUT::Print();
+
+ cout << "fEpsilon: " << fEpsilon << endl;
+ cout << "fIts: " << fIts << endl;
}
-void AliL3HoughTransformerVhdl::TransformLine()
+void AliL3HoughTransformerVhdl::PrintVhdl()
{
- //Do a line transform on the data.
-
-
- AliL3DigitRowData *tempPt = GetDataPointer();
- if(!tempPt)
- {
- LOG(AliL3Log::kError,"AliL3HoughTransformerVhdl::TransformLine","Data")
- <<"No input data "<<ENDLOG;
- return;
- }
-
- for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
- {
- AliL3DigitData *digPt = tempPt->fDigitData;
- if(i != (Int_t)tempPt->fRow)
- {
- printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
- continue;
- }
- for(UInt_t j=0; j<tempPt->fNDigit; j++)
- {
- UShort_t charge = digPt[j].fCharge;
- UChar_t pad = digPt[j].fPad;
- UShort_t time = digPt[j].fTime;
- if((Int_t)charge < GetLowerThreshold())
- continue;
- Int_t sector,row;
- Float_t xyz[3];
- AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
- AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
- Float_t eta = AliL3Transform::GetEta(xyz);
- Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
- if(eta_index < 0 || eta_index >= GetNEtaSegments())
- continue;
-
- //Get the correct histogram:
- AliL3Histogram *hist = fParamSpace[eta_index];
- if(!hist)
- {
- printf("AliL3HoughTransformerVhdl::TransformLine : Error getting histogram in index %d\n",eta_index);
- continue;
- }
- for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
- {
- Double_t theta = hist->GetBinCenterX(xbin);
- Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
- hist->Fill(theta,rho,charge);
- }
- }
- AliL3MemHandler::UpdateRowPointer(tempPt);
- }
+ cout << "fSlice := " << GetSlice() << ";" << endl;
+ cout << "fPatch := " << GetPatch() << ";" << endl;
+ // cout << "fSector := " << fSector << ";" << endl;
+ // cout << "fSectorRow := " << fSectorRow << ";" << endl;
+ // cout << "fMinRow := " << fMinRow << ";" << endl;
+ // cout << "fMaxRow := " << fMaxRow << ";" << endl;
+ // cout << "fNRows := " << fNRows << ";" << endl;
+ // cout << "fNEtas := " << fNEtas << ";" << endl;
+ // cout << "fNPhi0 := " << fNPhi0 << ";" << endl;
+ cout << "fZSign := " << fZSign << ";" << endl;
+ cout << "fZLengthPlusOff := " << fZLengthPlusOff << ";" << endl;
+ cout << "fPadPitch := " << fPadPitch << ";" << endl;
+ cout << "fTimeWidth := " << fTimeWidth << ";" << endl;
+ if(!fNRows) return;
+ cout << "fNLUTX :=" << fNRows << ";" << endl;
+ cout << "fLUTX := ( ";
+ for(Int_t i=0;i<fNRows-1;i++) cout << fLUTX[i] << ", ";
+ cout << fLUTX[fNRows-1] << " );\n" << endl;
+ cout << "fNLUTY := " << fNRows << ";" << endl;
+ cout << "fLUTY := ( ";
+ for(Int_t i=0;i<fNRows-1;i++) cout << fLUTY[i] << ", ";
+ cout << fLUTY[fNRows-1] << " );\n" << endl;
+
+ if(!fNEtas) return;
+ cout << "fNLUTEta := " << fNEtas << ";" << endl;
+ cout << "fLUTEta := (";
+ for(Int_t i=0;i<fNEtas-1;i++) cout << fLUTEta[i] << ", ";
+ cout << fLUTEta[fNEtas-1] << " );\n" << endl;
+
+ if(!fNPhi0) return;
+ cout << "fNxbin := " << fNxbin << endl;
+ cout << "fxmin := " << fXmin << endl;
+ cout << "fxmax := " << fXmax << endl;
+ cout << "fNybin := " << fNybin << endl;
+ cout << "fymin := " << fYmin << endl;
+ cout << "fymax := " << fYmax << endl;
+
+ cout << "\nfNLUTphi0 := " << fNPhi0 << ";" << endl;
+ cout << "fLUTphi0 := (";
+ for(Int_t i=0;i<fNPhi0-1;i++) cout << fLUTphi0[i] << ", ";
+ cout << fLUTphi0[fNPhi0-1] << " );\n" << endl;
+
+ cout << "fLUT2sinphi0 := (";
+ for(Int_t i=0;i<fNPhi0-1;i++) cout << fLUT2sinphi0[i] << ", ";
+ cout << fLUT2sinphi0[fNPhi0-1] << " );\n" << endl;
+
+ cout << "fLUT2cosphi0 := (";
+ for(Int_t i=0;i<fNPhi0;i++) cout << fLUT2cosphi0[i] << ", ";
+ cout << fLUT2cosphi0[fNPhi0-1] << " );\n" << endl;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff