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