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