Added Kappa information.
[u/mrichter/AliRoot.git] / HLT / hough / AliL3HoughTransformerVhdl.cxx
index 801e3ac1ea9655bce83439514168de06c983b275..02bdcb7a11f7a16817402b7c12a38d54349bdfd7 100644 (file)
 //$Id$
 
 // Author: Constantin Loizides <mailto:loizides@fi.uib.no>
-//*-- Copyright &copy 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
@@ -126,208 +103,253 @@ void AliL3HoughTransformerVhdl::TransformCircle()
        <<"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(charge <= GetThreshold())
-           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;
+         }
+
+         for(Int_t b=0; b<fNPhi0; b++){ 
 
-         //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);
+#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(charge < GetThreshold())
-           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 << "-- Kappa -- " << endl;
+  cout << "fNxbin := " << fNxbin << endl;
+  cout << "fxmin := " << fXmin << endl;
+  cout << "fxmax := " << fXmax << endl;
+  Float_t kdiff=(fXmax-fXmin)/fNxbin;
+  Float_t k=fXmin;
+  cout << "fKappa := ";
+  for(Int_t i=0;i<fNxbin;i++){
+    k+=kdiff;
+    cout << k << ", ";
+  }
+  cout << k+kdiff << " );\n" << endl;
+
+  cout << "-- Phi --" << 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 << "\nfNLUT2sinphi0 := " << fNPhi0 << ";" << endl;
+  cout << "fLUT2sinphi0 := (";
+  for(Int_t i=0;i<fNPhi0-1;i++) cout << fLUT2sinphi0[i] << ", ";
+  cout << fLUT2sinphi0[fNPhi0-1] << " );\n" << endl; 
+
+  cout << "\nfNLUT2cosphi0 := " << fNPhi0 << ";" << endl;
+  cout << "fLUT2cosphi0 := (";
+  for(Int_t i=0;i<fNPhi0;i++) cout << fLUT2cosphi0[i] << ", ";
+  cout << fLUT2cosphi0[fNPhi0-1] << " );\n" << endl; 
+
+  //cout << "\nfEpsilon := " << fEpsilon << endl;
+  //cout << "\nfSinEpsilon := " << fSinEpsilon << endl;
+  //cout << "\nfCosEpsilon := " << fCosEpsilon << endl;
 }