HLT/hough/AliL3HoughTransformerVhdl.cxx

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