HLT/hough/AliL3HoughTransformer.cxx

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