HLT/hough/AliL3HoughTransformerLUT.cxx

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