HLT/hough/AliL3HoughTransformerLUT.cxx

   1 // @(#) $Id$
   2
   3 // Author: Constantin Loizides <mailto:loizides@ikf.uni-frankfurt.de>
   4 //*-- Copyright &copy ALICE HLT Group
   5 /** /class AliL3HoughTransformerLUT
   6 //<pre>
   7 //_____________________________________________________________
   8 // AliL3HoughTransformerLUT
   9 //
  10 // Hough transformation class using Look-UP-Tables
  11 //
  12 //</pre>
  13 */
  14
  15 #include "AliL3StandardIncludes.h"
  16
  17 #include "AliL3Logging.h"
  18 #include "AliL3HoughTransformerLUT.h"
  19 #include "AliL3Transform.h"
  20 #include "AliL3MemHandler.h"
  21 #include "AliL3DigitData.h"
  22 #include "AliL3Histogram.h"
  23 #include "AliL3HistogramAdaptive.h"
  24
  25 #if __GNUC__ >= 3
  26 using namespace std;
  27 #endif
  28
  29 ClassImp(AliL3HoughTransformerLUT)
  30
  31 AliL3HoughTransformerLUT::AliL3HoughTransformerLUT() : AliL3HoughBaseTransformer()
  32 {
  33   // default contructor
  34   fParamSpace=0;
  35   fMinRow=0;
  36   fMaxRow=0;
  37
  38   fNRows=0;
  39   fNEtas=0;
  40   fNPhi0=0;
  41   fSector=0;
  42   fSlice=0;
  43   fSectorRow=0;
  44   fZSign=0;
  45   fZLengthPlusOff=0;
  46   fTimeWidth=0;
  47   fPadPitch=0;
  48   fEtaSlice=0;
  49
  50   fLUTX=0;
  51   fLUTY=0;
  52   fLUTEta=0;
  53   fLUTEtaReal=0;
  54   fLUTphi0=0;
  55   fLUT2sinphi0=0;
  56   fLUT2cosphi0=0;
  57   fLUTKappa=0;
  58
  59   fLastPad=0;
  60   fLastIndex=0;
  61   fAccCharge=0;
  62   fX=fY=0.;
  63 }
  64
  65 AliL3HoughTransformerLUT::AliL3HoughTransformerLUT(Int_t slice,Int_t patch,Int_t nEtaSegments) : AliL3HoughBaseTransformer(slice,patch,nEtaSegments)
  66 {
  67   // constructor
  68   fParamSpace=0;
  69   fMinRow=0;
  70   fMaxRow=0;
  71
  72   fNRows=0;
  73   fNEtas=0;
  74   fNPhi0=0;
  75   fSector=0;
  76   fSectorRow=0;
  77   fZSign=0;
  78   fZLengthPlusOff=0;
  79   fTimeWidth=0;
  80   fPadPitch=0;
  81   fEtaSlice=0;
  82
  83   fLUTX=0;
  84   fLUTY=0;
  85   fLUTEta=0;
  86   fLUTEtaReal=0;
  87   fLUTphi0=0;
  88   fLUT2sinphi0=0;
  89   fLUT2cosphi0=0;
  90   fLUTKappa=0;
  91
  92   fLastPad=0;
  93   fLastIndex=0;
  94   fAccCharge=0;
  95   fX=fY=0.;
  96
  97   Init(slice,patch,nEtaSegments);
  98 }
  99
 100 AliL3HoughTransformerLUT::~AliL3HoughTransformerLUT()
 101 {
 102   // destructor
 103   DeleteHistograms();
 104
 105   if(fNRows){
 106     delete[] fLUTX;
 107     delete[] fLUTY;
 108     fNRows=0;
 109   }
 110
 111   if(fNEtas){
 112     delete[] fLUTEta;
 113     delete[] fLUTEtaReal;
 114     fNEtas=0;
 115   }
 116 }
 117
 118 void AliL3HoughTransformerLUT::DeleteHistograms()
 119 {
 120   // deletes all histograms
 121   if(fNPhi0){
 122     delete[] fLUT2sinphi0;
 123     delete[] fLUT2cosphi0;
 124     delete[] fLUTphi0;
 125     delete[] fLUTKappa;
 126     fNPhi0=0;
 127     fLUTphi0=0;
 128     fLUT2sinphi0=0;
 129     fLUT2cosphi0=0;
 130     fLUTKappa=0;
 131   }
 132
 133   if(!fParamSpace){
 134     for(Int_t i=0; i<fNEtas; i++)
 135       {
 136         if(!fParamSpace[i]) continue;
 137         delete fParamSpace[i];
 138       }
 139     delete [] fParamSpace;
 140   }
 141 }
 142
 143 void AliL3HoughTransformerLUT::Init(Int_t slice,Int_t patch,Int_t nEtaSegments,Int_t /*nSeqs*/)
 144 {
 145   // Initialization
 146   AliL3HoughBaseTransformer::Init(slice,patch,nEtaSegments);
 147
 148   //delete old LUT tables
 149   if(fNRows){
 150     fNRows=0;
 151     delete[] fLUTX;
 152     delete[] fLUTY;
 153   }
 154   if(fNEtas){
 155     delete[] fLUTEta;
 156     delete[] fLUTEtaReal;
 157     fNEtas=0;
 158   }
 159
 160   //member values
 161   fMinRow=AliL3Transform::GetFirstRow(patch);;
 162   fMaxRow=AliL3Transform::GetLastRow(patch);;
 163   fNRows=AliL3Transform::GetNRows(patch);;
 164   fNEtas=nEtaSegments;
 165   if(fNEtas!=GetNEtaSegments()) {
 166     cerr << "AliL3HoughTransformerLUT::Init -> Error: Number of Etas must be equal!" << endl;
 167     exit(1);
 168   }
 169
 170   AliL3Transform::Slice2Sector(slice,fMinRow,fSector,fSectorRow);
 171   fZSign = slice < 18 ? 1:-1; //see AliL3Transformer
 172   fSlice = slice;
 173   fZLengthPlusOff=AliL3Transform::GetZLength()+AliL3Transform::GetZOffset();
 174   fTimeWidth=AliL3Transform::GetZWidth();
 175
 176   fPadPitch=0;
 177   if(fSector<AliL3Transform::GetNSectorLow())
 178     fPadPitch=AliL3Transform::GetPadPitchWidthLow();
 179   else
 180     fPadPitch=AliL3Transform::GetPadPitchWidthUp();
 181
 182   Float_t etaMax=GetEtaMax();
 183   Float_t etaMin=GetEtaMin();
 184   fEtaSlice=(etaMax-etaMin)/nEtaSegments;
 185
 186   //lookup tables for X and Y
 187   fLUTX=new Float_t[fNRows];
 188   fLUTY=new Float_t[fNRows];
 189   for(Int_t rr=0;rr<fNRows;rr++){
 190     fLUTX[rr]=Float_t(AliL3Transform::Row2X(rr+fMinRow));
 191     fLUTY[rr]=Float_t(0.5*(AliL3Transform::GetNPads(rr+fMinRow)-1)*fPadPitch);
 192     //cout << rr << ": " << (Float_t)fLUTX[rr] << " " << (Float_t)fLUTY[rr] << endl;
 193   }
 194
 195   //lookup tables for rz2s <=> etas
 196   //will only be used ifdef FULLLUT
 197   fLUTEta=new Float_t[fNEtas];
 198   fLUTEtaReal=new Float_t[fNEtas];
 199   for(Int_t rr=0;rr<fNEtas;rr++){
 200     Float_t eta=etaMin+(rr+1)*fEtaSlice;
 201     fLUTEta[rr]=CalcRoverZ2(eta);
 202     fLUTEtaReal[rr]=eta-0.5*fEtaSlice;
 203     //cout << rr << ": " << eta << " " << fLUTEtaReal[rr] << " " << GetEta(rr,fSlice) << " - " << fLUTEta[rr] << endl;
 204   }
 205 }
 206
 207 void AliL3HoughTransformerLUT::CreateHistograms(Float_t ptmin,Float_t ptmax,Float_t ptres,
 208                                              Int_t nybin,Float_t psi)
 209 {
 210   //Create histograms.
 211   //_Only_ to be used in case of the adaptive histograms!
 212   //phimax is given in radians!!
 213
 214   if(ptmin > ptmax)
 215     {
 216       cerr<<"AliL3HoughTransformerLUT::CreateHistograms: Error in ptrange "<<ptmin<<" "<<ptmax<<endl;
 217       return;
 218     }
 219   if(psi < 0)
 220     {
 221       cerr<<"AliL3HoughTransformerLUT::CreateHistograms: Wrong psi-angle "<<psi<<endl;
 222       return;
 223     }
 224
 225   fParamSpace = new AliL3Histogram*[fNEtas];
 226   Char_t histname[256];
 227   for(Int_t i=0; i<GetNEtaSegments(); i++)
 228     {
 229       sprintf(histname,"paramspace_%d",i);
 230       fParamSpace[i] = new AliL3HistogramAdaptive(histname,ptmin,ptmax,ptres,nybin,-psi,psi);
 231     }
 232
 233   //create lookup table for sin and cos
 234   fNPhi0=nybin;
 235   fLUTphi0=new Float_t[fNPhi0];
 236   fLUT2sinphi0=new Float_t[fNPhi0];
 237   fLUT2cosphi0=new Float_t[fNPhi0];
 238   fLUTKappa=new Float_t[fNPhi0];
 239   AliL3Histogram *hist=fParamSpace[0];
 240   Int_t i=0;
 241   for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
 242     {
 243       Float_t phi0 = hist->GetBinCenterY(b);
 244       fLUTphi0[i]=phi0;
 245       fLUT2sinphi0[i]=2.*sin(phi0);
 246       fLUT2cosphi0[i]=2.*cos(phi0);
 247       fLUTKappa[i]=0.;
 248       i++;
 249       //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
 250     }
 251 }
 252
 253 void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Float_t ptMin,Int_t nybin,Float_t phimin,Float_t phimax)
 254 {
 255   //Create the histograms (parameter space).
 256   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
 257   //The arguments give the range and binning;
 258   //nxbin = #bins in kappa
 259   //nybin = #bins in phi0
 260   //ptMin = mimium Pt of track (corresponding to maximum kappa)
 261   //phi_min = mimimum phi0 (degrees)
 262   //phi_max = maximum phi0 (degrees)
 263
 264   Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/ptMin;
 265   //Double_t torad = AliL3Transform::Pi()/180;
 266   CreateHistograms(nxbin,-1.*x,x,nybin,phimin/**torad*/,phimax/**torad*/);
 267 }
 268
 269 void AliL3HoughTransformerLUT::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,Int_t nybin,Float_t ymin,Float_t ymax)
 270 {
 271   fParamSpace = new AliL3Histogram*[fNEtas];
 272
 273   Char_t histname[256];
 274   for(Int_t i=0; i<fNEtas; i++)
 275     {
 276       sprintf(histname,"paramspace_%d",i);
 277       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
 278     }
 279
 280   //create lookup table for sin and cos
 281   fNPhi0=nybin;
 282   fLUTphi0=new Float_t[fNPhi0];
 283   fLUT2sinphi0=new Float_t[fNPhi0];
 284   fLUT2cosphi0=new Float_t[fNPhi0];
 285   fLUTKappa=new Float_t[fNPhi0];
 286
 287   AliL3Histogram *hist=fParamSpace[0];
 288   Int_t i=0;
 289   for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
 290     {
 291       Float_t phi0 = hist->GetBinCenterY(b);
 292       fLUTphi0[i]=phi0;
 293       fLUT2sinphi0[i]=2.*sin(phi0);
 294       fLUT2cosphi0[i]=2.*cos(phi0);
 295       fLUTKappa[i]=0.;
 296       //cout << i << ": " << fLUTphi0[i] << " " << fLUT2sinphi0[i] << " " << fLUT2cosphi0[i] << endl;
 297       i++;
 298     }
 299 }
 300
 301 void AliL3HoughTransformerLUT::Reset()
 302 {
 303   //Reset all the histograms
 304
 305   if(!fParamSpace)
 306     {
 307       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
 308         <<"No histograms to reset"<<ENDLOG;
 309       return;
 310     }
 311
 312   for(Int_t i=0; i<fNEtas; i++)
 313     fParamSpace[i]->Reset();
 314 }
 315
 316 Int_t AliL3HoughTransformerLUT::GetEtaIndex(Double_t eta) const
 317 {
 318   //Return the histogram index of the corresponding eta.
 319
 320 #ifdef FULLLUT
 321   /* try to imitate a circuit -> should
 322      go into the VHDL implementation of transformer */
 323   Float_t rz2=CalcRoverZ2(eta);
 324   return FindIndex(rz2);
 325 #else /* optimize for speed on the computer */
 326   Double_t index = (eta-GetEtaMin())/fEtaSlice;
 327   return (Int_t)index;
 328 #endif
 329 }
 330
 331 AliL3Histogram *AliL3HoughTransformerLUT::GetHistogram(Int_t etaIndex)
 332 {
 333   // gets hitogram
 334   if(!fParamSpace || etaIndex >= fNEtas || etaIndex < 0)
 335     return 0;
 336   if(!fParamSpace[etaIndex])
 337     return 0;
 338   return fParamSpace[etaIndex];
 339 }
 340
 341 Double_t AliL3HoughTransformerLUT::GetEta(Int_t etaIndex,Int_t slice) const
 342 {
 343   // gets eta
 344   if(etaIndex >= fNEtas || etaIndex < 0){
 345     LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::GetEta","Index") << "Index out of range."<<ENDLOG;
 346     return 0.;
 347   }
 348   if(slice != fSlice){
 349     LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::GetEta","Index") << "Given slice does not match internal slice."<<ENDLOG;
 350     return 0.;
 351   }
 352
 353   return(fLUTEtaReal[etaIndex]);
 354 }
 355
 356 void AliL3HoughTransformerLUT::TransformCircle()
 357 {
 358   //Transform the input data with a circle HT.
 359   //The function loops over all the data, and transforms each pixel with the equation:
 360   //
 361   //kappa = 2/R^2 * (y*cos(phi0) - x*sin(phi0) )
 362   //
 363   //where R^2 = x^2 + y^2
 364   //
 365   //After a day of testing it is proven that this h
 366
 367   AliL3DigitRowData *tempPt = GetDataPointer();
 368   if(!tempPt)
 369     {
 370       LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
 371         <<"No input data "<<ENDLOG;
 372       return;
 373     }
 374
 375   Int_t lowch=GetLowerThreshold();
 376   Int_t highch=GetUpperThreshold();
 377
 378   //Loop over the padrows:
 379   for(Int_t i=fMinRow, row=0; i<=fMaxRow; i++, row++)
 380     {
 381       //Get the data on this padrow:
 382       AliL3DigitData *digPt = tempPt->fDigitData;
 383       if(i != (Int_t)tempPt->fRow)
 384         {
 385           LOG(AliL3Log::kError,"AliL3HoughTransformerLUT::TransformCircle","Data")
 386             <<"AliL3HoughTransformerLUT::TransformCircle : Mismatching padrow numbering "<<i<<" != "<<(Int_t)tempPt->fRow<<ENDLOG;
 387           continue;
 388         }
 389
 390
 391       //start a new row
 392       fLastPad=-1;
 393       //store x for this row
 394       fX = CalcX(row);
 395       //accumulate charge per etaslice
 396       fAccCharge=0;
 397       //last histogram
 398       fLastIndex=-1;
 399
 400       //Loop over the data on this padrow:
 401       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 402         {
 403           UShort_t charge = digPt[j].fCharge;
 404
 405           //check threshold
 406           if((Int_t)charge <=  lowch)
 407             continue;
 408           if ((Int_t)charge > highch)
 409             charge=highch;
 410
 411           UChar_t pad = digPt[j].fPad;
 412           UShort_t time = digPt[j].fTime;
 413
 414           if(fLastPad!=pad){ //only update if necessary
 415
 416             //if there is accumalated charge,
 417             //put it in the old histogram
 418             //using the old X and Y values
 419             if(fAccCharge>0) AddCurveToHistogram(-1); //fLastIndex will be set to -1
 420
 421 #ifdef FULLLUT
 422             fLastIndex=fNEtas-1;
 423 #endif
 424
 425             //calculate Y for the new pad
 426             fY = CalcY(pad,row);
 427             //remember new pad
 428             fLastPad=pad;
 429           }
 430
 431 #ifdef FULLLUT
 432           //find eta slice
 433           Float_t z = CalcZ(time);
 434
 435           Float_t z2=z*z;
 436           Float_t rz2 = 1 + r2/z2;
 437           Int_t etaIndex = FindIndex(rz2,fLastIndex);
 438 #else
 439           Float_t z = CalcZ(time);
 440           Double_t r = sqrt(fX*fX+fY*fY+z*z);
 441           Double_t eta = 0.5 * log((r+z)/(r-z));
 442           Int_t etaIndex = GetEtaIndex(eta);
 443 #endif
 444           if(etaIndex < 0 || etaIndex >= fNEtas){
 445             LOG(AliL3Log::kWarning,"AliL3HoughTransformerLUT::TransformCircle","Histograms")<<"No histograms corresponding to eta index value of "<<etaIndex<<"."<<ENDLOG;
 446             continue;
 447           }
 448
 449 #ifndef FULLLUT
 450           if(fLastIndex==-1) fLastIndex=etaIndex;
 451 #endif
 452
 453           if(fLastIndex!=etaIndex){ //enter old values first
 454             AddCurveToHistogram(etaIndex);
 455           }
 456
 457           fAccCharge+=charge;
 458           //fAccCharge+=1;
 459         }
 460
 461       //in case there is charge, store it before moving to the next row
 462       if(fAccCharge>0) AddCurveToHistogram(-1);
 463
 464       //Move the data pointer to the next padrow:
 465       AliL3MemHandler::UpdateRowPointer(tempPt);
 466     }
 467 }
 468
 469 void AliL3HoughTransformerLUT::Print()
 470 {
 471   // debug printout
 472   cout << "fSlice: " << GetSlice() << endl;
 473   cout << "fPatch: " << GetPatch() << endl;
 474   cout << "fSector: " << fSector << endl;
 475   cout << "fSectorRow: " << fSectorRow << endl;
 476   cout << "fMinRow: " << fMinRow << endl;
 477   cout << "fMaxRow: " << fMaxRow << endl;
 478   cout << "fNRows: " << fNRows << endl;
 479   cout << "fNEtas: " << fNEtas << endl;
 480   cout << "fNPhi0: " << fNPhi0 << endl;
 481   cout << "fZSign: " << fZSign << endl;
 482   cout << "fZLengthPlusOff: " << fZLengthPlusOff << endl;
 483   cout << "fPadPitch: " << fPadPitch << endl;
 484   cout << "fTimeWidth: " << fTimeWidth << endl;
 485
 486   if(!fNRows) return;
 487   cout << "fLUTX " << fNRows << endl;
 488   for(Int_t i=0;i<fNRows;i++) cout << "fLUTX[" << i << "]=" << (Float_t)fLUTX[i] << endl;
 489   cout << "fLUTY " << fNRows << endl;
 490   for(Int_t i=0;i<fNRows;i++) cout << "fLUTY[" << i << "]=" << fLUTY[i] << endl;
 491   if(!fNEtas) return;
 492   cout << "fLUTEta " << fNEtas << endl;
 493   for(Int_t i=0;i<fNEtas;i++) cout << "fLUTEta[" << i << "]=" << fLUTEta[i] << endl;
 494   cout << "fLUTEtaReal " << fNEtas << endl;
 495   for(Int_t i=0;i<fNEtas;i++) cout << "fLUTEtaReal[" << i << "]=" << fLUTEtaReal[i] << endl;
 496   if(!fNPhi0) return;
 497   cout << "fLUTphi0 " << fNPhi0 << endl;
 498   for(Int_t i=0;i<fNPhi0;i++) cout << "fLUTPhi0[" << i << "]=" << fLUTphi0[i] << endl;
 499   cout << "fLUT2sinphi0 " << fNPhi0 << endl;
 500   for(Int_t i=0;i<fNPhi0;i++) cout << "fLUT2sinphi0[" << i << "]=" << fLUT2sinphi0[i] << endl;
 501   cout << "fLUT2cosphi0 " << fNPhi0 << endl;
 502   for(Int_t i=0;i<fNPhi0;i++) cout << "fLUT2cosphi0[" << i << "]=" << fLUT2cosphi0[i] << endl;
 503 }