HLT/hough/AliL3HoughTransformerLUT.cxx

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