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