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