HLT/hough/AliL3HoughTransformer.cxx

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