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