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