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