HLT/hough/AliL3HoughTransformer.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6 #include "AliL3MemHandler.h"
   7 #include "AliL3Logging.h"
   8 #include "AliL3HoughTransformer.h"
   9 #include "AliL3Transform.h"
  10 #include "AliL3DigitData.h"
  11 #include "AliL3Histogram.h"
  12
  13 //_____________________________________________________________
  14 // AliL3HoughTransformer
  15 //
  16 // Hough transformation class
  17 //
  18
  19 ClassImp(AliL3HoughTransformer)
  20
  21 AliL3HoughTransformer::AliL3HoughTransformer()
  22 {
  23   //Default constructor
  24   fParamSpace = 0;
  25 #ifdef do_mc
  26   fTrackID = 0;
  27 #endif
  28 }
  29
  30 AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
  31 {
  32   //Normal constructor
  33   fParamSpace = 0;
  34 #ifdef do_mc
  35   fTrackID = 0;
  36 #endif
  37 }
  38
  39 AliL3HoughTransformer::~AliL3HoughTransformer()
  40 {
  41   DeleteHistograms();
  42 #ifdef do_mc
  43   if(fTrackID)
  44     {
  45       for(Int_t i=0; i<GetNEtaSegments(); i++)
  46         {
  47           if(!fTrackID[i]) continue;
  48           delete fTrackID[i];
  49         }
  50       delete [] fTrackID;
  51     }
  52 #endif
  53 }
  54
  55 //void AliL3HoughTransformer::Init(Int_t slice=0,Int_t patch=0,Int_t n_eta_segments=100){}
  56
  57 void AliL3HoughTransformer::DeleteHistograms()
  58 {
  59   if(!fParamSpace)
  60     return;
  61   for(Int_t i=0; i<GetNEtaSegments(); i++)
  62     {
  63       if(!fParamSpace[i]) continue;
  64       delete fParamSpace[i];
  65     }
  66   delete [] fParamSpace;
  67 }
  68
  69 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min,
  70                                              Int_t nybin,Double_t phimin,Double_t phimax)
  71 {
  72   //Create the histograms (parameter space).
  73   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
  74   //The arguments give the range and binning;
  75   //nxbin = #bins in kappa
  76   //nybin = #bins in phi0
  77   //pt_min = mimium Pt of track (corresponding to maximum kappa)
  78   //phi_min = mimimum phi0 (degrees)
  79   //phi_max = maximum phi0 (degrees)
  80
  81   Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
  82   Double_t torad = AliL3Transform::Pi()/180;
  83   CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
  84 }
  85
  86 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
  87                                              Int_t nybin,Double_t ymin,Double_t ymax)
  88 {
  89
  90   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
  91
  92   Char_t histname[256];
  93   for(Int_t i=0; i<GetNEtaSegments(); i++)
  94     {
  95       sprintf(histname,"paramspace_%d",i);
  96       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
  97     }
  98
  99 #ifdef do_mc
 100   Int_t ncells = (nxbin+2)*(nybin+2);
 101   cout<<"Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
 102   fTrackID = new TrackIndex*[GetNEtaSegments()];
 103   for(Int_t i=0; i<GetNEtaSegments(); i++)
 104     fTrackID[i] = new TrackIndex[ncells];
 105 #endif
 106 }
 107
 108 void AliL3HoughTransformer::Reset()
 109 {
 110   //Reset all the histograms
 111
 112   if(!fParamSpace)
 113     {
 114       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
 115         <<"No histograms to reset"<<ENDLOG;
 116       return;
 117     }
 118
 119   for(Int_t i=0; i<GetNEtaSegments(); i++)
 120     fParamSpace[i]->Reset();
 121 #ifdef do_mc
 122   AliL3Histogram *hist = fParamSpace[0];
 123   Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
 124   for(Int_t i=0; i<GetNEtaSegments(); i++)
 125     memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
 126 #endif
 127 }
 128
 129 Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta)
 130 {
 131   //Return the histogram index of the corresponding eta.
 132
 133   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
 134   Double_t index = (eta-GetEtaMin())/etaslice;
 135   return (Int_t)index;
 136 }
 137
 138 inline AliL3Histogram *AliL3HoughTransformer::GetHistogram(Int_t eta_index)
 139 {
 140   if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
 141     return 0;
 142   if(!fParamSpace[eta_index])
 143     return 0;
 144   return fParamSpace[eta_index];
 145 }
 146
 147 Double_t AliL3HoughTransformer::GetEta(Int_t eta_index,Int_t slice)
 148 {
 149   Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
 150   Double_t eta=(Double_t)((eta_index+0.5)*eta_slice);
 151   if(slice>17) eta*=-1;
 152   return eta;
 153 }
 154
 155 void AliL3HoughTransformer::TransformCircle()
 156 {
 157   //Transform the input data with a circle HT.
 158   //The function loops over all the data, and transforms each pixel with the equations:
 159   //
 160   //kappa = 2/R*sin(phi - phi0)
 161   //
 162   //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
 163   //
 164   //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
 165   //which histogram in which the pixel should be transformed, the eta-value is calcluated
 166   //and the proper histogram index is found by GetEtaIndex(eta).
 167
 168
 169   AliL3DigitRowData *tempPt = GetDataPointer();
 170   if(!tempPt)
 171     {
 172       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
 173         <<"No input data "<<ENDLOG;
 174       return;
 175     }
 176
 177   //Loop over the padrows:
 178   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 179     {
 180       //Get the data on this padrow:
 181       AliL3DigitData *digPt = tempPt->fDigitData;
 182       if(i != (Int_t)tempPt->fRow)
 183         {
 184           cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl;
 185           continue;
 186         }
 187
 188       //Loop over the data on this padrow:
 189       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 190         {
 191           UShort_t charge = digPt[j].fCharge;
 192           UChar_t pad = digPt[j].fPad;
 193           UShort_t time = digPt[j].fTime;
 194           if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
 195             continue;
 196           Int_t sector,row;
 197           Float_t xyz[3];
 198
 199           //Transform data to local cartesian coordinates:
 200           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 201           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 202
 203           //Calculate the eta:
 204           Double_t eta = AliL3Transform::GetEta(xyz);
 205
 206           //Get the corresponding index, which determines which histogram to fill:
 207           Int_t eta_index = GetEtaIndex(eta);
 208           if(eta_index < 0 || eta_index >= GetNEtaSegments())
 209             continue;
 210
 211           //Get the correct histogrampointer:
 212           AliL3Histogram *hist = fParamSpace[eta_index];
 213           if(!hist)
 214             {
 215               printf("AliL3HoughTransformer::TransformCircle : Error getting histogram in index %d\n",eta_index);
 216               continue;
 217             }
 218
 219           //Do the transformation:
 220           Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 221           Float_t phi = AliL3Transform::GetPhi(xyz);
 222
 223           //Fill the histogram along the phirange
 224           for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
 225             {
 226               Float_t phi0 = hist->GetBinCenterY(b);
 227               Float_t kappa = 2*sin(phi - phi0)/R;
 228               hist->Fill(kappa,phi0,charge);
 229 #ifdef do_mc
 230               Int_t bin = hist->FindBin(kappa,phi0);
 231               for(Int_t t=0; t<3; t++)
 232                 {
 233                   Int_t label = digPt[j].fTrackID[t];
 234                   if(label < 0) break;
 235                   UInt_t c;
 236                   for(c=0; c<MaxTrack; c++)
 237                     if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0)
 238                       break;
 239                   if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
 240                   fTrackID[eta_index][bin].fLabel[c] = label;
 241                   fTrackID[eta_index][bin].fNHits[c]++;
 242                 }
 243 #endif
 244             }
 245         }
 246
 247       //Move the data pointer to the next padrow:
 248       AliL3MemHandler::UpdateRowPointer(tempPt);
 249     }
 250 }
 251
 252 void AliL3HoughTransformer::TransformCircleC(Int_t row_range)
 253 {
 254   //Circle transform, using combinations of every 2 points lying
 255   //on different padrows and within the same etaslice.
 256
 257   AliL3DigitRowData *tempPt = GetDataPointer();
 258   if(!tempPt)
 259     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
 260       <<"No input data "<<ENDLOG;
 261
 262   Int_t counter=0;
 263   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 264     {
 265       counter += tempPt->fNDigit;
 266       AliL3MemHandler::UpdateRowPointer(tempPt);
 267     }
 268
 269   struct Digit {
 270     Int_t row;
 271     Double_t r;
 272     Double_t phi;
 273     Int_t eta_index;
 274     Int_t charge;
 275     Int_t trackID[3];
 276   };
 277
 278   Digit *digits = new Digit[counter];
 279   cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
 280
 281   Int_t total_digits=counter;
 282   Int_t sector,row,tot_charge,pad,time,charge;
 283   Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
 284   Float_t xyz[3];
 285
 286   counter=0;
 287   tempPt = GetDataPointer();
 288
 289   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 290     {
 291       AliL3DigitData *digPt = tempPt->fDigitData;
 292       for(UInt_t di=0; di<tempPt->fNDigit; di++)
 293         {
 294           charge = digPt[di].fCharge;
 295           pad = digPt[di].fPad;
 296           time = digPt[di].fTime;
 297           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 298           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 299           eta = AliL3Transform::GetEta(xyz);
 300           digits[counter].row = i;
 301           digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 302           digits[counter].phi = atan2(xyz[1],xyz[0]);
 303           digits[counter].eta_index = GetEtaIndex(eta);
 304           digits[counter].charge = charge;
 305           digits[counter].trackID[0] = digPt[di].fTrackID[0];
 306           digits[counter].trackID[1] = digPt[di].fTrackID[1];
 307           digits[counter].trackID[2] = digPt[di].fTrackID[2];
 308           counter++;
 309         }
 310       AliL3MemHandler::UpdateRowPointer(tempPt);
 311     }
 312
 313   for(Int_t i=0; i<total_digits; i++)
 314     {
 315       if(digits[i].eta_index < 0 || digits[i].eta_index >= GetNEtaSegments()) continue;
 316       Int_t ind = digits[i].eta_index;
 317
 318       for(Int_t j=i+1; j<total_digits; j++)
 319         {
 320           if(digits[i].row == digits[j].row) continue;
 321           if(digits[i].eta_index != digits[j].eta_index) continue;
 322           if(digits[i].row + row_range < digits[j].row) break;
 323
 324           //Get the correct histogrampointer:
 325           AliL3Histogram *hist = fParamSpace[ind];
 326           if(!hist)
 327             {
 328               printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",ind);
 329               continue;
 330             }
 331
 332           r1 = digits[i].r;
 333           phi1 = digits[i].phi;
 334           r2 = digits[j].r;
 335           phi2 = digits[j].phi;
 336           phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
 337           kappa = 2*sin(phi2-phi_0)/r2;
 338           tot_charge = digits[i].charge + digits[j].charge;
 339           hist->Fill(kappa,phi_0,tot_charge);
 340 #ifdef do_mc
 341           Int_t bin = hist->FindBin(kappa,phi_0);
 342           for(Int_t l=0; l<3; l++)
 343             {
 344               for(Int_t m=0; m<3; m++)
 345                 {
 346                   if(digits[i].trackID[l] == digits[j].trackID[m])
 347                     {
 348                       Int_t label = digits[i].trackID[l];
 349                       if(label < 0) continue;
 350                       UInt_t c;
 351                       for(c=0; c<MaxTrack; c++)
 352                         if(fTrackID[ind][bin].fLabel[c] == label || fTrackID[ind][bin].fNHits[c] == 0)
 353                           break;
 354                       if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircleC : Array reached maximum!! "<<c<<endl;
 355                       fTrackID[ind][bin].fLabel[c] = label;
 356                       fTrackID[ind][bin].fNHits[c]++;
 357                     }
 358                 }
 359             }
 360 #endif
 361         }
 362     }
 363   delete [] digits;
 364 }
 365
 366 void AliL3HoughTransformer::TransformLine()
 367 {
 368   //Do a line transform on the data.
 369
 370
 371   AliL3DigitRowData *tempPt = GetDataPointer();
 372   if(!tempPt)
 373     {
 374       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data")
 375         <<"No input data "<<ENDLOG;
 376       return;
 377     }
 378
 379   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 380     {
 381       AliL3DigitData *digPt = tempPt->fDigitData;
 382       if(i != (Int_t)tempPt->fRow)
 383         {
 384           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
 385           continue;
 386         }
 387       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 388         {
 389           UShort_t charge = digPt[j].fCharge;
 390           UChar_t pad = digPt[j].fPad;
 391           UShort_t time = digPt[j].fTime;
 392           if(charge < GetLowerThreshold())
 393             continue;
 394           Int_t sector,row;
 395           Float_t xyz[3];
 396           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 397           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 398           Float_t eta = AliL3Transform::GetEta(xyz);
 399           Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
 400           if(eta_index < 0 || eta_index >= GetNEtaSegments())
 401             continue;
 402
 403           //Get the correct histogram:
 404           AliL3Histogram *hist = fParamSpace[eta_index];
 405           if(!hist)
 406             {
 407               printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
 408               continue;
 409             }
 410           for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
 411             {
 412               Double_t theta = hist->GetBinCenterX(xbin);
 413               Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
 414               hist->Fill(theta,rho,charge);
 415             }
 416         }
 417       AliL3MemHandler::UpdateRowPointer(tempPt);
 418     }
 419
 420 }
 421
 422 Int_t AliL3HoughTransformer::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
 423 {
 424
 425 #ifdef do_mc
 426   if(eta_index < 0 || eta_index > GetNEtaSegments())
 427     {
 428       cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl;
 429       return -1;
 430     }
 431   AliL3Histogram *hist = fParamSpace[eta_index];
 432   Int_t bin = hist->FindBin(kappa,psi);
 433   Int_t label=-1;
 434   Int_t max=0;
 435   for(UInt_t i=0; i<MaxTrack; i++)
 436     {
 437       Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
 438       if(nhits == 0) break;
 439       if(nhits > max)
 440         {
 441           max = nhits;
 442           label = fTrackID[eta_index][bin].fLabel[i];
 443         }
 444     }
 445   return label;
 446 #endif
 447   cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
 448   return -1;
 449 }
 450