HLT/hough/AliL3HoughTransformerVhdl.cxx

   1 //$Id$
   2
   3 // Author: Constantin Loizides <mailto:loizides@fi.uib.no>
   4 //*-- Copyright&Copy CL
   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
  15 #if GCCVERSION == 3
  16 using namespace std;
  17 #endif
  18
  19 /** \class AliL3HoughTransformerVhdl
  20 // <pre>
  21 //_____________________________________________________________
  22 // AliL3HoughTransformerVhdl
  23 //
  24 // Hough transformation class for VHDL comparism.
  25 //
  26 //</pre>
  27 */
  28
  29 ClassImp(AliL3HoughTransformerVhdl)
  30
  31
  32 AliL3HoughTransformerVhdl::AliL3HoughTransformerVhdl()
  33   : AliL3HoughTransformerLUT()
  34
  35 {
  36   fEpsilon=0;
  37   fSinEpsilon=0;
  38   fCosEpsilon=1;
  39   fIts=0;
  40 }
  41
  42 AliL3HoughTransformerVhdl::AliL3HoughTransformerVhdl(Int_t slice,Int_t patch,Int_t n_eta_segments,Int_t n_its)
  43   : AliL3HoughTransformerLUT(slice,patch,n_eta_segments)
  44 {
  45   fEpsilon=0;
  46   fSinEpsilon=0;
  47   fCosEpsilon=1;
  48   fIts=n_its;
  49 }
  50
  51 AliL3HoughTransformerVhdl::~AliL3HoughTransformerVhdl()
  52 {
  53
  54 }
  55
  56 void AliL3HoughTransformerVhdl::CreateHistograms(Int_t nxbin,Double_t pt_min,Int_t nybin,Double_t phimin,Double_t phimax)
  57 {
  58   AliL3HoughTransformerLUT::CreateHistograms(nxbin,pt_min,nybin,phimin,phimax);
  59 }
  60
  61 void AliL3HoughTransformerVhdl::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
  62                                                  Int_t nybin,Double_t ymin,Double_t ymax)
  63 {
  64   AliL3HoughTransformerLUT::CreateHistograms(nxbin,xmin,xmax,nybin,ymin,ymax);
  65
  66   fEpsilon=(ymax-ymin)/nybin;
  67   fSinEpsilon=sin(fEpsilon);
  68   fCosEpsilon=cos(fEpsilon);
  69
  70   //  cout << fEpsilon << " - " << (xmax-xmin)/nxbin << endl;
  71 }
  72
  73 void AliL3HoughTransformerVhdl::Init(Int_t slice,Int_t patch,Int_t n_eta_segments,Int_t n_its)
  74 {
  75   AliL3HoughTransformerLUT::Init(slice,patch,n_eta_segments);
  76 }
  77
  78 void AliL3HoughTransformerVhdl::TransformCircle()
  79 {
  80   //Transform the input data with a circle HT.
  81   //The function loops over all the data, and transforms each pixel with the equation:
  82   //
  83   //kappa      = lastkappa +- epsilon * lastkappaprime
  84   //
  85   //kappaprime = lastkappaprime -+ epsilon * lastkappa
  86   //
  87   //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
  88   //which histogram in which the pixel should be transformed, the eta-value is calcluated
  89   //and the proper histogram index is found by GetEtaIndex(eta).
  90
  91   AliL3DigitRowData *tempPt = GetDataPointer();
  92   if(!tempPt)
  93     {
  94       LOG(AliL3Log::kError,"AliL3HoughTransformerVhdl::TransformCircle","Data")
  95         <<"No input data "<<ENDLOG;
  96       return;
  97     }
  98
  99 //#define use_error
 100 #ifdef use_error
 101       Float_t max_error=0;
 102       Float_t rel_error=0;
 103       Int_t counter=0;
 104 #endif
 105
 106   //Loop over the padrows:
 107   for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
 108     {
 109       //Get the data on this padrow:
 110       AliL3DigitData *digPt = tempPt->fDigitData;
 111       if(i != (Int_t)tempPt->fRow)
 112         {
 113           LOG(AliL3Log::kError,"AliL3HoughTransformerVhdl::TransformCircle","Data")
 114             <<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
 115           continue;
 116         }
 117
 118       //Loop over the data on this padrow:
 119       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 120         {
 121           UShort_t charge = digPt[j].fCharge;
 122
 123           //check threshold
 124           if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
 125             continue;
 126
 127           UChar_t pad = digPt[j].fPad;
 128           UShort_t time = digPt[j].fTime;
 129
 130           Float_t x = CalcX(row);
 131           Float_t y = CalcY(pad,row);
 132           Float_t z = CalcZ(time);
 133
 134           //find eta slice
 135           Float_t rz2 = 1 + (x*x+y*y)/(z*z);
 136           Int_t eta_index = FindIndex(rz2);
 137
 138           if(eta_index < 0 || eta_index >= fNEtas){
 139             //LOG(AliL3Log::kWarning,"AliL3HoughTransformerVhdl::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<eta_index<<"."<<ENDLOG;
 140             continue;
 141           }
 142
 143           //Get the correct histogrampointer:
 144           AliL3Histogram *hist = fParamSpace[eta_index];
 145           if(!hist){
 146             //LOG(AliL3Log::kWarning,"AliL3HoughTransformerVhdl::TransformCircle","Histograms")<<"Error getting histogram in index "<<eta_index<<"."<<ENDLOG;
 147             continue;
 148           }
 149
 150           //Fill the histogram along the phirange
 151           Float_t kappa=0;
 152           Float_t kappaprime=0;
 153           Float_t lastkappa=0;
 154           Float_t lastkappaprime=0;
 155           Int_t its=0;
 156           Float_t R2=1/(x*x+y*y);
 157           Float_t A=2*y*R2;
 158           Float_t B=-2*x*R2;
 159
 160 #if 1
 161           //starting point
 162           Float_t phi = fLUTphi0[fNPhi0/2];
 163           kappa=A*cos(phi)+B*sin(phi);
 164           kappaprime=B*cos(phi)-A*sin(phi);
 165           its=fIts;
 166           lastkappa=kappa;
 167           lastkappaprime=kappaprime;
 168           hist->Fill(kappa,phi,charge);
 169
 170           for(Int_t b=fNPhi0/2+1; b<fNPhi0; b++){
 171             Float_t exact_kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
 172
 173             phi=fLUTphi0[b];
 174
 175             if(its==0) { //initialize or re-adjust values
 176               kappa=A*cos(phi)+B*sin(phi); //equals exact_kappa!
 177               kappaprime=B*cos(phi)-A*sin(phi);
 178               its=fIts;
 179             } else {
 180               //kappa=fCosEpsilon*lastkappa+fSinEpsilon*lastkappaprime;
 181               //kappaprime=fCosEpsilon*lastkappaprime-fSinEpsilon*lastkappa;
 182               kappa=lastkappa+fEpsilon*lastkappaprime;
 183               kappaprime=lastkappaprime-fEpsilon*lastkappa;
 184               its--;
 185             }
 186
 187             lastkappa=kappa;
 188             lastkappaprime=kappaprime;
 189
 190             hist->Fill(kappa,phi,charge);
 191           } // end positive running values
 192
 193           phi = fLUTphi0[fNPhi0/2];
 194           kappa=A*cos(phi)+B*sin(phi);
 195           kappaprime=B*cos(phi)-A*sin(phi);
 196           its=fIts;
 197           lastkappa=kappa;
 198           lastkappaprime=kappaprime;
 199           //hist->Fill(kappa,fLUTphi0[b],charge);
 200
 201           for(Int_t b=fNPhi0/2-1; b>=0; b--){
 202             Float_t exact_kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
 203
 204             Float_t phi = fLUTphi0[b];
 205
 206             if(its==0) { //initialize or re-adjust values
 207               kappa=A*cos(phi)+B*sin(phi); //equals exact_kappa!
 208               kappaprime=B*cos(phi)-A*sin(phi);
 209               its=fIts;
 210             } else {
 211               //kappa=fCosEpsilon*lastkappa-fSinEpsilon*lastkappaprime;
 212               //kappaprime=fCosEpsilon*lastkappaprime+fSinEpsilon*lastkappa;
 213               kappa=lastkappa-fEpsilon*lastkappaprime;
 214               kappaprime=lastkappaprime+fEpsilon*lastkappa;
 215               its--;
 216             }
 217
 218             lastkappa=kappa;
 219             lastkappaprime=kappaprime;
 220
 221             hist->Fill(kappa,phi,charge);
 222           }
 223
 224 #else
 225           for(Int_t b=0; b<fNPhi0; b++){
 226             Float_t exact_kappa=R2*(y*fLUT2cosphi0[b]-x*fLUT2sinphi0[b]);
 227
 228             Float_t phi = fLUTphi0[b];
 229
 230             if(its==0) { //initialize or re-adjust values
 231               kappa=A*cos(phi)+B*sin(phi); //equals exact_kappa!
 232               kappaprime=B*cos(phi)-A*sin(phi);
 233               its=fIts;
 234             } else {
 235               //kappa=fCosEpsilon*lastkappa+fSinEpsilon*lastkappaprime;
 236               //kappaprime=fCosEpsilon*lastkappaprime-fSinEpsilon*lastkappa;
 237               kappa=lastkappa+fEpsilon*lastkappaprime;
 238               kappaprime=lastkappaprime-fEpsilon*lastkappa;
 239               its--;
 240             }
 241
 242 #ifdef use_error
 243             Float_t err=fabs((exact_kappa-kappa)/exact_kappa);
 244
 245             if(err>max_error) {
 246               //cout << max_error << " - " << err << " " << kappa << " " << exact_kappa << " " << kappa/exact_kappa<< endl;
 247               max_error=err;
 248             }
 249
 250             rel_error+=err;
 251             counter++;
 252 #endif
 253
 254             lastkappa=kappa;
 255             lastkappaprime=kappaprime;
 256
 257             hist->Fill(kappa,fLUTphi0[b],charge);
 258
 259             //cout << kappa << " " << fLUTphi0[b] << " " << charge << endl;
 260           }
 261 #endif /*test*/
 262
 263         }
 264       //Move the data pointer to the next padrow:
 265       AliL3MemHandler::UpdateRowPointer(tempPt);
 266     }
 267
 268 #ifdef use_error
 269   cout <<"Max Error: " << max_error << endl;
 270   cout <<"Rel Error average: " << rel_error/counter << endl;
 271 #endif
 272 }
 273
 274 void AliL3HoughTransformerVhdl::Print()
 275 {
 276   AliL3HoughTransformerLUT::Print();
 277
 278   cout << "fEpsilon: " << fEpsilon << endl;
 279   cout << "fIts: " << fIts << endl;
 280 }