HLT/hough/AliL3HoughTransformer.cxx

   1 // @(#) $Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ALICE HLT Group
   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   fEtaOverlap=kFALSE;
  33 #ifdef do_mc
  34   fTrackID = 0;
  35 #endif
  36 }
  37
  38 AliL3HoughTransformer::AliL3HoughTransformer(Int_t slice,Int_t patch,Int_t n_eta_segments,Bool_t DoEtaOverlap,Bool_t DoMC) : AliL3HoughBaseTransformer(slice,patch,n_eta_segments)
  39 {
  40   //Normal constructor
  41   fParamSpace = 0;
  42   fDoMC = kFALSE;
  43   fEtaOverlap = DoEtaOverlap;
  44   fDoMC=kFALSE;
  45   if(DoMC==kTRUE)
  46     {
  47       if(patch==0)
  48         fDoMC = kTRUE;
  49       else
  50         fDoMC = kFALSE;
  51     }
  52 #ifdef do_mc
  53   fTrackID = 0;
  54 #endif
  55 }
  56
  57 AliL3HoughTransformer::~AliL3HoughTransformer()
  58 {
  59   DeleteHistograms();
  60 #ifdef do_mc
  61   if(fTrackID)
  62     {
  63       for(Int_t i=0; i<GetNEtaSegments(); i++)
  64         {
  65           if(!fTrackID[i]) continue;
  66           delete fTrackID[i];
  67         }
  68       delete [] fTrackID;
  69     }
  70 #endif
  71 }
  72
  73 void AliL3HoughTransformer::DeleteHistograms()
  74 {
  75   if(!fParamSpace)
  76     return;
  77   for(Int_t i=0; i<GetNEtaSegments(); i++)
  78     {
  79       if(!fParamSpace[i]) continue;
  80       delete fParamSpace[i];
  81     }
  82   delete [] fParamSpace;
  83 }
  84
  85 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t pt_min,
  86                                              Int_t nybin,Double_t phimin,Double_t phimax)
  87 {
  88   //Create the histograms (parameter space).
  89   //These are 2D histograms, span by kappa (curvature of track) and phi0 (emission angle with x-axis).
  90   //The arguments give the range and binning;
  91   //nxbin = #bins in kappa
  92   //nybin = #bins in phi0
  93   //pt_min = mimium Pt of track (corresponding to maximum kappa)
  94   //phi_min = mimimum phi0 (degrees)
  95   //phi_max = maximum phi0 (degrees)
  96
  97   Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/pt_min;
  98   Double_t torad = AliL3Transform::Pi()/180;
  99   CreateHistograms(nxbin,-1.*x,x,nybin,phimin*torad,phimax*torad);
 100 }
 101
 102 void AliL3HoughTransformer::CreateHistograms(Int_t nxbin,Double_t xmin,Double_t xmax,
 103                                              Int_t nybin,Double_t ymin,Double_t ymax)
 104 {
 105
 106   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
 107
 108   Char_t histname[256];
 109   for(Int_t i=0; i<GetNEtaSegments(); i++)
 110     {
 111       sprintf(histname,"paramspace_%d",i);
 112       //fParamSpace[i] = new AliL3HistogramAdaptive(histname,0.1,1.5,0.05,nybin,ymin,ymax);
 113       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
 114     }
 115
 116 #ifdef do_mc
 117   if(fDoMC)
 118     {
 119       AliL3Histogram *hist = fParamSpace[0];
 120       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
 121       cout<<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(TrackIndex)<<" bytes to fTrackID"<<endl;
 122       fTrackID = new TrackIndex*[GetNEtaSegments()];
 123       for(Int_t i=0; i<GetNEtaSegments(); i++)
 124         fTrackID[i] = new TrackIndex[ncells];
 125     }
 126 #endif
 127 }
 128
 129 void AliL3HoughTransformer::Reset()
 130 {
 131   //Reset all the histograms
 132
 133   if(!fParamSpace)
 134     {
 135       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
 136         <<"No histograms to reset"<<ENDLOG;
 137       return;
 138     }
 139
 140   for(Int_t i=0; i<GetNEtaSegments(); i++)
 141     fParamSpace[i]->Reset();
 142 #ifdef do_mc
 143   if(fDoMC)
 144     {
 145       AliL3Histogram *hist = fParamSpace[0];
 146       Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
 147       for(Int_t i=0; i<GetNEtaSegments(); i++)
 148         memset(fTrackID[i],0,ncells*sizeof(TrackIndex));
 149     }
 150 #endif
 151 }
 152
 153 Int_t AliL3HoughTransformer::GetEtaIndex(Double_t eta)
 154 {
 155   //Return the histogram index of the corresponding eta.
 156
 157   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
 158   Double_t index = (eta-GetEtaMin())/etaslice;
 159   return (Int_t)index;
 160 }
 161
 162 void AliL3HoughTransformer::GetEtaIndexes(Double_t eta,Int_t *indexes)
 163 {
 164   //Return histogram indexes in case of overlapping etaslices.
 165
 166   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
 167   Int_t index = (Int_t)((eta-GetEtaMin())/etaslice);
 168   if(index%2 == 0)
 169     {
 170       indexes[0] = index;
 171       indexes[1] = index - 1;
 172     }
 173   else
 174     {
 175       indexes[0] = index - 1;
 176       indexes[1] = index;
 177     }
 178 }
 179
 180 inline AliL3Histogram *AliL3HoughTransformer::GetHistogram(Int_t eta_index)
 181 {
 182   if(!fParamSpace || eta_index >= GetNEtaSegments() || eta_index < 0)
 183     return 0;
 184   if(!fParamSpace[eta_index])
 185     return 0;
 186   return fParamSpace[eta_index];
 187 }
 188
 189 Double_t AliL3HoughTransformer::GetEta(Int_t eta_index,Int_t slice)
 190 {
 191   Double_t eta_slice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
 192   Double_t eta=0;
 193   if(fEtaOverlap)
 194     {
 195       Int_t index = eta_index + 1;
 196       eta=(Double_t)((index)*eta_slice);
 197     }
 198   else
 199     eta=(Double_t)((eta_index+0.5)*eta_slice);
 200   return eta;
 201 }
 202
 203
 204 void AliL3HoughTransformer::TransformCircle()
 205 {
 206   //Transform the input data with a circle HT.
 207   //The function loops over all the data, and transforms each pixel with the equations:
 208   //
 209   //kappa = 2/R*sin(phi - phi0)
 210   //
 211   //where R = sqrt(x*x +y*y), and phi = arctan(y/x)
 212   //
 213   //Each pixel then transforms into a curve in the (kappa,phi0)-space. In order to find
 214   //which histogram in which the pixel should be transformed, the eta-value is calculated
 215   //and the proper histogram index is found by GetEtaIndex(eta).
 216
 217
 218   AliL3DigitRowData *tempPt = GetDataPointer();
 219   if(!tempPt)
 220     {
 221       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
 222         <<"No input data "<<ENDLOG;
 223       return;
 224     }
 225
 226   //Loop over the padrows:
 227   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 228     {
 229       //Get the data on this padrow:
 230       AliL3DigitData *digPt = tempPt->fDigitData;
 231       if(i != (Int_t)tempPt->fRow)
 232         {
 233           cerr<<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<i<<" "<<(Int_t)tempPt->fRow<<endl;
 234           continue;
 235         }
 236
 237       //Loop over the data on this padrow:
 238       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 239         {
 240           UShort_t charge = digPt[j].fCharge;
 241           UChar_t pad = digPt[j].fPad;
 242           UShort_t time = digPt[j].fTime;
 243           if((Int_t)charge <= GetLowerThreshold() || (Int_t)charge > GetUpperThreshold())
 244             continue;
 245
 246           Int_t sector,row;
 247           Float_t xyz[3];
 248
 249           //Transform data to local cartesian coordinates:
 250           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 251           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 252
 253           //Calculate the eta:
 254           Double_t eta = AliL3Transform::GetEta(xyz);
 255
 256           //Get the corresponding index, which determines which histogram to fill:
 257           Int_t eta_index = GetEtaIndex(eta);
 258
 259           if(eta_index < 0 || eta_index >= GetNEtaSegments())
 260             continue;
 261
 262           //Get the correct histogrampointer:
 263           AliL3Histogram *hist = fParamSpace[eta_index];
 264           if(!hist)
 265             {
 266               cerr<<"AliL3HoughTransformer::TransformCircle : Error getting histogram in index "<<eta_index<<endl;
 267               continue;
 268             }
 269
 270           //Do the transformation:
 271           Float_t R = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 272           Float_t phi = AliL3Transform::GetPhi(xyz);
 273
 274           //Fill the histogram along the phirange
 275           for(Int_t b=hist->GetFirstYbin(); b<=hist->GetLastYbin(); b++)
 276             {
 277               Float_t phi0 = hist->GetBinCenterY(b);
 278               Float_t kappa = 2*sin(phi - phi0)/R;
 279               hist->Fill(kappa,phi0,charge);
 280               //hist->Fill(kappa,phi0,1);
 281 #ifdef do_mc
 282               if(fDoMC)
 283                 {
 284                   Int_t bin = hist->FindBin(kappa,phi0);
 285                   for(Int_t t=0; t<3; t++)
 286                     {
 287                       Int_t label = digPt[j].fTrackID[t];
 288                       if(label < 0) break;
 289                       UInt_t c;
 290                       for(c=0; c<MaxTrack; c++)
 291                         if(fTrackID[eta_index][bin].fLabel[c] == label || fTrackID[eta_index][bin].fNHits[c] == 0)
 292                           break;
 293                       if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
 294                       fTrackID[eta_index][bin].fLabel[c] = label;
 295                       fTrackID[eta_index][bin].fNHits[c]++;
 296                     }
 297                 }
 298 #endif
 299             }
 300         }
 301
 302       //Move the data pointer to the next padrow:
 303       AliL3MemHandler::UpdateRowPointer(tempPt);
 304     }
 305 }
 306
 307
 308 struct Digit {
 309   Int_t row;
 310   Double_t r;
 311   Double_t phi;
 312   Int_t charge;
 313   Digit *next;
 314 };
 315 struct EtaContainer {
 316   Digit *first;
 317   Digit *last;
 318 };
 319 void AliL3HoughTransformer::TransformCircleC(Int_t *row_range,Int_t every)
 320 {
 321   //Circle transform, using combinations of every 2 points lying
 322   //on different padrows and within the same etaslice.
 323
 324   AliL3DigitRowData *tempPt = GetDataPointer();
 325   if(!tempPt)
 326     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
 327       <<"No input data "<<ENDLOG;
 328
 329   Int_t minrow = AliL3Transform::GetFirstRow(GetPatch());
 330   Int_t maxrow = AliL3Transform::GetLastRow(GetPatch());
 331   if(row_range)
 332     {
 333       minrow = row_range[0];
 334       maxrow = row_range[1];
 335       if(minrow < AliL3Transform::GetFirstRow(GetPatch()) || minrow >= AliL3Transform::GetLastRow(GetPatch()))
 336         minrow = AliL3Transform::GetFirstRow(GetPatch());
 337       if(maxrow < AliL3Transform::GetFirstRow(GetPatch()) || maxrow >= AliL3Transform::GetLastRow(GetPatch()))
 338         maxrow = AliL3Transform::GetLastRow(GetPatch());
 339       if(minrow > maxrow || maxrow==minrow)
 340         {
 341           cerr<<"AliL3HoughTransformer::TransformCircleC : Bad row range "<<minrow<<" "<<maxrow<<endl;
 342           return;
 343         }
 344     }
 345   else
 346     {
 347       minrow = AliL3Transform::GetFirstRow(GetPatch());
 348       maxrow = AliL3Transform::GetLastRow(GetPatch());
 349     }
 350
 351   Int_t counter=0;
 352   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 353     {
 354       counter += tempPt->fNDigit;
 355       AliL3MemHandler::UpdateRowPointer(tempPt);
 356     }
 357
 358   Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1);
 359   EtaContainer *etaPt = new EtaContainer[bound];
 360   memset(etaPt,0,bound*sizeof(EtaContainer));
 361
 362   Digit *digits = new Digit[counter];
 363   cout<<"Allocating "<<counter*sizeof(Digit)<<" bytes to digitsarray"<<endl;
 364   memset(digits,0,counter*sizeof(Digit));
 365
 366   Int_t sector,row,tot_charge,pad,time,charge;
 367   Double_t r1,r2,phi1,phi2,eta,kappa,phi_0;
 368   Float_t xyz[3];
 369
 370   counter=0;
 371   tempPt = GetDataPointer();
 372
 373   cout<<"Calculating digits in patch "<<GetPatch()<<endl;
 374   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 375     {
 376       AliL3DigitData *digPt = tempPt->fDigitData;
 377       for(UInt_t di=0; di<tempPt->fNDigit; di++)
 378         {
 379           charge = digPt[di].fCharge;
 380           pad = digPt[di].fPad;
 381           time = digPt[di].fTime;
 382           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 383           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 384           eta = AliL3Transform::GetEta(xyz);
 385
 386           digits[counter].row = i;
 387           digits[counter].r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 388           digits[counter].phi = atan2(xyz[1],xyz[0]);
 389           digits[counter].charge = charge;
 390
 391           if(!fEtaOverlap)
 392             {
 393               Int_t eta_index = GetEtaIndex(eta);
 394
 395               Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + eta_index;
 396
 397               if(index > 0 && index < bound)
 398                 {
 399                   if(etaPt[index].first == 0)
 400                     etaPt[index].first = &digits[counter];
 401                   else
 402                     (etaPt[index].last)->next = &digits[counter];
 403                   etaPt[index].last = &digits[counter];
 404                 }
 405             }
 406           else
 407             {
 408               Int_t eta_index[2];
 409               GetEtaIndexes(eta,eta_index);
 410               Int_t index[2];
 411               index[0] = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + eta_index[0];
 412               index[1] = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + eta_index[1];
 413               if(index[0] == index[1])
 414                 {
 415                   cerr<<"Same etaindexes "<<index[0]<<" "<<index[1]<<endl;
 416                   exit(5);
 417                 }
 418
 419               Int_t ind = index[0];
 420               if(ind > 0 && ind < bound)
 421                 {
 422                   if(etaPt[ind].first == 0)
 423                     etaPt[ind].first = &digits[counter];
 424                   else
 425                     (etaPt[ind].last)->next = &digits[counter];
 426                   etaPt[ind].last = &digits[counter];
 427                 }
 428
 429               ind = index[1];
 430               if(ind > 0 && ind < bound)
 431                 {
 432                   if(etaPt[ind].first == 0)
 433                     etaPt[ind].first = &digits[counter];
 434                   else
 435                     (etaPt[ind].last)->next = &digits[counter];
 436                   etaPt[ind].last = &digits[counter];
 437                 }
 438             }
 439
 440           counter++;
 441         }
 442       AliL3MemHandler::UpdateRowPointer(tempPt);
 443     }
 444
 445   cout<<"Doing the combinatorics"<<endl;
 446
 447   Digit *dPt1,*dPt2;
 448
 449   for(Int_t e=0; e<GetNEtaSegments(); e++)
 450     {
 451       for(Int_t i=minrow; i<=maxrow; i+=every)
 452         {
 453           Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e;
 454
 455           for(dPt1 = (Digit*)etaPt[index1].first; dPt1 != 0; dPt1 = (Digit*)dPt1->next)
 456             {
 457               for(Int_t j=i+every; j<=maxrow; j+=every)
 458                 {
 459                   Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e;
 460
 461                   for(dPt2 = (Digit*)etaPt[index2].first; dPt2 != 0; dPt2 = (Digit*)dPt2->next)
 462                     {
 463                       if(dPt1->row == dPt2->row)
 464                         {
 465                           cerr<<"same row; indexes "<<index1<<" "<<index2<<endl;
 466                           exit(5);
 467                         }
 468
 469                       //Get the correct histogrampointer:
 470                       AliL3Histogram *hist = fParamSpace[e];
 471                       if(!hist)
 472                         {
 473                           printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",i);
 474                           continue;
 475                         }
 476
 477                       //Do the transform:
 478                       r1 = dPt1->r;
 479                       phi1 = dPt1->phi;
 480                       r2 = dPt2->r;
 481                       phi2 = dPt2->phi;
 482                       phi_0 = atan( (r2*sin(phi1)-r1*sin(phi2))/(r2*cos(phi1)-r1*cos(phi2)) );
 483                       kappa = 2*sin(phi2-phi_0)/r2;
 484                       tot_charge = dPt1->charge + dPt2->charge;
 485                       hist->Fill(kappa,phi_0,tot_charge);
 486
 487                     }
 488                 }
 489             }
 490         }
 491     }
 492
 493   cout<<"done"<<endl;
 494   delete [] etaPt;
 495   delete [] digits;
 496
 497 }
 498
 499 void AliL3HoughTransformer::TransformLine(Int_t *row_range,Float_t *phirange)
 500 {
 501   //Do a line transform on the data.
 502
 503
 504   AliL3DigitRowData *tempPt = GetDataPointer();
 505   if(!tempPt)
 506     {
 507       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformLine","Data")
 508         <<"No input data "<<ENDLOG;
 509       return;
 510     }
 511
 512   Int_t minrow = AliL3Transform::GetFirstRow(GetPatch());
 513   Int_t maxrow = AliL3Transform::GetLastRow(GetPatch());
 514   if(row_range)
 515     {
 516       minrow = row_range[0];
 517       maxrow = row_range[1];
 518       if(minrow < AliL3Transform::GetFirstRow(GetPatch()) || minrow >= AliL3Transform::GetLastRow(GetPatch()))
 519         minrow = AliL3Transform::GetFirstRow(GetPatch());
 520       if(maxrow < AliL3Transform::GetFirstRow(GetPatch()) || maxrow >= AliL3Transform::GetLastRow(GetPatch()))
 521         maxrow = AliL3Transform::GetLastRow(GetPatch());
 522       if(minrow > maxrow || maxrow==minrow)
 523         {
 524           cerr<<"AliL3HoughTransformer::TransformCircleC : Bad row range "<<minrow<<" "<<maxrow<<endl;
 525           return;
 526         }
 527     }
 528
 529   for(Int_t i=minrow; i<=maxrow; i++)
 530     {
 531       AliL3DigitData *digPt = tempPt->fDigitData;
 532       if(i != (Int_t)tempPt->fRow)
 533         {
 534           printf("AliL3HoughTransform::TransformLine : Mismatching padrow numbering\n");
 535           continue;
 536         }
 537       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 538         {
 539           UShort_t charge = digPt[j].fCharge;
 540           UChar_t pad = digPt[j].fPad;
 541           UShort_t time = digPt[j].fTime;
 542           if(charge < GetLowerThreshold())
 543             continue;
 544           Int_t sector,row;
 545           Float_t xyz[3];
 546           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 547           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 548
 549           if(phirange)
 550             {
 551               Float_t phi = AliL3Transform::GetPhi(xyz);
 552               if(phi < phirange[0] || phi > phirange[1])
 553                 continue;
 554             }
 555           Float_t eta = AliL3Transform::GetEta(xyz);
 556           Int_t eta_index = GetEtaIndex(eta);//(Int_t)(eta/etaslice);
 557           if(eta_index < 0 || eta_index >= GetNEtaSegments())
 558             continue;
 559
 560           xyz[0] = xyz[0] - AliL3Transform::Row2X(minrow);
 561
 562           //Get the correct histogram:
 563           AliL3Histogram *hist = fParamSpace[eta_index];
 564           if(!hist)
 565             {
 566               printf("AliL3HoughTransformer::TransformLine : Error getting histogram in index %d\n",eta_index);
 567               continue;
 568             }
 569           for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
 570             {
 571               Double_t theta = hist->GetBinCenterX(xbin);
 572               Double_t rho = xyz[0]*cos(theta) + xyz[1]*sin(theta);
 573               hist->Fill(theta,rho,charge);
 574             }
 575         }
 576       AliL3MemHandler::UpdateRowPointer(tempPt);
 577     }
 578
 579 }
 580
 581 struct LDigit {
 582   Int_t row;
 583   Int_t charge;
 584   Float_t y;
 585   LDigit *next;
 586 };
 587 struct LEtaContainer {
 588   LDigit *first;
 589   LDigit *last;
 590 };
 591 void AliL3HoughTransformer::TransformLineC(Int_t *rowrange,Float_t *phirange)
 592 {
 593   AliL3DigitRowData *tempPt = GetDataPointer();
 594   if(!tempPt)
 595     LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircleC","Data")
 596       <<"No input data "<<ENDLOG;
 597
 598
 599   Int_t counter=0;
 600   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 601     {
 602       counter += tempPt->fNDigit;
 603       AliL3MemHandler::UpdateRowPointer(tempPt);
 604     }
 605
 606   Int_t bound = (GetNEtaSegments()+1)*(AliL3Transform::GetNRows(GetPatch())+1);
 607   LEtaContainer *etaPt = new LEtaContainer[bound];
 608   memset(etaPt,0,bound*sizeof(LEtaContainer));
 609
 610   LDigit *digits = new LDigit[counter];
 611   cout<<"Allocating "<<counter*sizeof(LDigit)<<" bytes to digitsarray"<<endl;
 612   memset(digits,0,counter*sizeof(LDigit));
 613
 614   Int_t sector,row;
 615   Float_t xyz[3];
 616
 617   counter=0;
 618   tempPt = GetDataPointer();
 619
 620   cout<<"Calculating digits in patch "<<GetPatch()<<endl;
 621   for(Int_t i=AliL3Transform::GetFirstRow(GetPatch()); i<=AliL3Transform::GetLastRow(GetPatch()); i++)
 622     {
 623       AliL3DigitData *digPt = tempPt->fDigitData;
 624       for(UInt_t di=0; di<tempPt->fNDigit; di++)
 625         {
 626           Int_t charge = digPt[di].fCharge;
 627           Int_t pad = digPt[di].fPad;
 628           Int_t time = digPt[di].fTime;
 629           AliL3Transform::Slice2Sector(GetSlice(),i,sector,row);
 630           AliL3Transform::Raw2Local(xyz,sector,row,(Int_t)pad,(Int_t)time);
 631           Double_t eta = AliL3Transform::GetEta(xyz);
 632
 633           Float_t phi = atan2(xyz[1],xyz[0]);
 634           if(phi < phirange[0] || phi > phirange[1]) continue;
 635
 636           digits[counter].row = i;
 637           digits[counter].y = xyz[1];
 638           digits[counter].charge = charge;
 639
 640           Int_t eta_index = GetEtaIndex(eta);
 641           Int_t index = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + eta_index;
 642
 643           if(index > 0 && index < bound)
 644             {
 645               if(etaPt[index].first == 0)
 646                 etaPt[index].first = &digits[counter];
 647               else
 648                 (etaPt[index].last)->next = &digits[counter];
 649               etaPt[index].last = &digits[counter];
 650             }
 651           counter++;
 652         }
 653       AliL3MemHandler::UpdateRowPointer(tempPt);
 654     }
 655
 656   cout<<"Doing the combinatorics"<<endl;
 657
 658   LDigit *dPt1,*dPt2;
 659
 660   for(Int_t e=0; e<GetNEtaSegments(); e++)
 661     {
 662       for(Int_t i=rowrange[0]; i<=rowrange[1]; i++)
 663         {
 664           Int_t index1 = (GetNEtaSegments()+1)*(i-AliL3Transform::GetFirstRow(GetPatch())) + e;
 665
 666           for(dPt1 = (LDigit*)etaPt[index1].first; dPt1 != 0; dPt1 = (LDigit*)dPt1->next)
 667             {
 668               for(Int_t j=i+1; j<=rowrange[1]; j++)
 669                 {
 670                   Int_t index2 = (GetNEtaSegments()+1)*(j-AliL3Transform::GetFirstRow(GetPatch())) + e;
 671
 672                   for(dPt2 = (LDigit*)etaPt[index2].first; dPt2 != 0; dPt2 = (LDigit*)dPt2->next)
 673                     {
 674                       if(dPt1->row == dPt2->row)
 675                         {
 676                           cerr<<"same row; indexes "<<index1<<" "<<index2<<endl;
 677                           exit(5);
 678                         }
 679
 680                       //Get the correct histogrampointer:
 681                       AliL3Histogram *hist = fParamSpace[e];
 682                       if(!hist)
 683                         {
 684                           printf("AliL3HoughTransformer::TransformCircleC() : No histogram at index %d\n",i);
 685                           continue;
 686                         }
 687
 688                       //Do the transform:
 689                       float x1 = AliL3Transform::Row2X(dPt1->row) - AliL3Transform::Row2X(rowrange[0]);
 690                       float x2 = AliL3Transform::Row2X(dPt2->row) - AliL3Transform::Row2X(rowrange[0]);
 691                       float y1 = dPt1->y;
 692                       float y2 = dPt2->y;
 693                       float theta = atan2(x2-x1,y1-y2);
 694                       float rho = x1*cos(theta)+y1*sin(theta);
 695                       hist->Fill(theta,rho,1);//dPt1->charge+dPt2->charge);
 696                     }
 697                 }
 698             }
 699         }
 700     }
 701
 702   cout<<"done"<<endl;
 703   delete [] etaPt;
 704   delete [] digits;
 705 }
 706
 707 Int_t AliL3HoughTransformer::GetTrackID(Int_t eta_index,Double_t kappa,Double_t psi)
 708 {
 709   if(!fDoMC)
 710     {
 711       cerr<<"AliL3HoughTransformer::GetTrackID : Flag switched off"<<endl;
 712       return -1;
 713     }
 714
 715 #ifdef do_mc
 716   if(eta_index < 0 || eta_index > GetNEtaSegments())
 717     {
 718       cerr<<"AliL3HoughTransformer::GetTrackID : Wrong etaindex "<<eta_index<<endl;
 719       return -1;
 720     }
 721   AliL3Histogram *hist = fParamSpace[eta_index];
 722   Int_t bin = hist->FindBin(kappa,psi);
 723   Int_t label=-1;
 724   Int_t max=0;
 725   for(UInt_t i=0; i<MaxTrack; i++)
 726     {
 727       Int_t nhits=fTrackID[eta_index][bin].fNHits[i];
 728       if(nhits == 0) break;
 729       if(nhits > max)
 730         {
 731           max = nhits;
 732           label = fTrackID[eta_index][bin].fLabel[i];
 733         }
 734     }
 735   //nhits = max;
 736   return label;
 737 #endif
 738   cout<<"AliL3HoughTransformer::GetTrackID : Compile with do_mc flag!"<<endl;
 739   return -1;
 740 }
 741