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