HLT/hough/AliL3HoughTransformerRow.cxx

   1 // @(#) $Id$
   2
   3
   4 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   5 //*-- Copyright &copy ALICE HLT Group
   6
   7 #include "AliL3StandardIncludes.h"
   8
   9 #include "AliL3Logging.h"
  10 #include "AliL3MemHandler.h"
  11 #include "AliL3Transform.h"
  12 #include "AliL3DigitData.h"
  13 #include "AliL3HistogramAdaptive.h"
  14 #include "AliL3HoughTrack.h"
  15 #include "AliL3HoughTransformerRow.h"
  16
  17 #if GCCVERSION == 3
  18 using namespace std;
  19 #endif
  20
  21 //_____________________________________________________________
  22 // AliL3HoughTransformerRow
  23 //
  24 // TPC rows Hough transformation class
  25 //
  26
  27 ClassImp(AliL3HoughTransformerRow)
  28
  29 UChar_t **AliL3HoughTransformerRow::fgRowCount = 0;
  30 UChar_t **AliL3HoughTransformerRow::fgGapCount = 0;
  31 UChar_t **AliL3HoughTransformerRow::fgCurrentRowCount = 0;
  32 #ifdef do_mc
  33 AliL3TrackIndex **AliL3HoughTransformerRow::fgTrackID = 0;
  34 #endif
  35 UChar_t *AliL3HoughTransformerRow::fgTrackNRows = 0;
  36 UChar_t *AliL3HoughTransformerRow::fgTrackFirstRow = 0;
  37 UChar_t *AliL3HoughTransformerRow::fgTrackLastRow = 0;
  38
  39 Float_t AliL3HoughTransformerRow::fgBeta1 = AliL3Transform::Row2X(79)/pow(AliL3Transform::Row2X(79),2);
  40 Float_t AliL3HoughTransformerRow::fgBeta2 = (AliL3Transform::Row2X(158)+6.0)/pow((AliL3Transform::Row2X(158)+6.0),2);
  41
  42 AliL3HoughTransformerRow::AliL3HoughTransformerRow()
  43 {
  44   //Default constructor
  45   fParamSpace = 0;
  46   fDoMC = kFALSE;;
  47
  48   fLUTforwardZ=0;
  49   fLUTforwardZ2=0;
  50   fLUTbackwardZ=0;
  51   fLUTbackwardZ2=0;
  52 }
  53
  54 AliL3HoughTransformerRow::AliL3HoughTransformerRow(Int_t slice,Int_t patch,Int_t netasegments,Bool_t /*DoMC*/,Float_t zvertex) : AliL3HoughBaseTransformer(slice,patch,netasegments,zvertex)
  55 {
  56   //Normal constructor
  57   fParamSpace = 0;
  58   fDoMC = kFALSE;
  59   fDoMC=kFALSE;
  60 #ifdef do_mc
  61   fDoMC = kTRUE;
  62 #endif
  63
  64   fLUTforwardZ=0;
  65   fLUTforwardZ2=0;
  66   fLUTbackwardZ=0;
  67   fLUTbackwardZ2=0;
  68 }
  69
  70 AliL3HoughTransformerRow::~AliL3HoughTransformerRow()
  71 {
  72   //Destructor
  73   DeleteHistograms();
  74 #ifdef do_mc
  75   if(fgTrackID)
  76     {
  77       for(Int_t i=0; i<GetNEtaSegments(); i++)
  78         {
  79           if(!fgTrackID[i]) continue;
  80           delete fgTrackID[i];
  81         }
  82       delete [] fgTrackID;
  83       fgTrackID = 0;
  84     }
  85 #endif
  86
  87   if(fgRowCount)
  88     {
  89       for(Int_t i=0; i<GetNEtaSegments(); i++)
  90         {
  91           if(!fgRowCount[i]) continue;
  92           delete [] fgRowCount[i];
  93         }
  94       delete [] fgRowCount;
  95       fgRowCount = 0;
  96     }
  97   if(fgGapCount)
  98     {
  99       for(Int_t i=0; i<GetNEtaSegments(); i++)
 100         {
 101           if(!fgGapCount[i]) continue;
 102           delete [] fgGapCount[i];
 103         }
 104       delete [] fgGapCount;
 105       fgGapCount = 0;
 106     }
 107   if(fgCurrentRowCount)
 108     {
 109       for(Int_t i=0; i<GetNEtaSegments(); i++)
 110         {
 111           if(fgCurrentRowCount[i])
 112             delete [] fgCurrentRowCount[i];
 113         }
 114       delete [] fgCurrentRowCount;
 115       fgCurrentRowCount = 0;
 116     }
 117   if(fgTrackNRows)
 118     {
 119       delete [] fgTrackNRows;
 120       fgTrackNRows = 0;
 121     }
 122   if(fgTrackFirstRow)
 123     {
 124       delete [] fgTrackFirstRow;
 125       fgTrackFirstRow = 0;
 126     }
 127   if(fgTrackLastRow)
 128     {
 129       delete [] fgTrackLastRow;
 130       fgTrackLastRow = 0;
 131     }
 132 }
 133
 134 void AliL3HoughTransformerRow::DeleteHistograms()
 135 {
 136   // Clean up
 137   delete[] fLUTforwardZ;
 138   delete[] fLUTforwardZ2;
 139   delete[] fLUTbackwardZ;
 140   delete[] fLUTbackwardZ2;
 141
 142   fLUTforwardZ=0;
 143   fLUTforwardZ2=0;
 144   fLUTbackwardZ=0;
 145   fLUTbackwardZ2=0;
 146
 147   if(!fParamSpace)
 148     return;
 149   for(Int_t i=0; i<GetNEtaSegments(); i++)
 150     {
 151       if(!fParamSpace[i]) continue;
 152       delete fParamSpace[i];
 153     }
 154   delete [] fParamSpace;
 155 }
 156
 157 void AliL3HoughTransformerRow::CreateHistograms(Int_t nxbin,Float_t ptmin,
 158                                                 Int_t nybin,Float_t phimin,Float_t phimax)
 159 {
 160   //Create the histograms (parameter space).
 161   //These are 2D histograms, span by kappa (curvature of track)
 162   //and phi0 (emission angle with x-axis).
 163   //The arguments give the range and binning;
 164   //nxbin = #bins in kappa
 165   //nybin = #bins in phi0
 166   //ptmin = mimium Pt of track (corresponding to maximum kappa)
 167   //phi_min = mimimum phi0
 168   //phi_max = maximum phi0
 169
 170   Double_t x = AliL3Transform::GetBFact()*AliL3Transform::GetBField()/ptmin;
 171   //Double_t torad = AliL3Transform::Pi()/180;
 172   CreateHistograms(nxbin,-1.*x,x,nybin,phimin/**torad*/,phimax/**torad*/);
 173 }
 174
 175 void AliL3HoughTransformerRow::CreateHistograms(Int_t nxbin,Float_t xmin,Float_t xmax,
 176                                                 Int_t nybin,Float_t ymin,Float_t ymax)
 177 {
 178   //Create the histograms (parameter space)
 179   //nxbin = #bins in X
 180   //nybin = #bins in Y
 181   //xmin xmax ymin ymax = histogram limits in X and Y
 182   fParamSpace = new AliL3Histogram*[GetNEtaSegments()];
 183
 184   Char_t histname[256];
 185   for(Int_t i=0; i<GetNEtaSegments(); i++)
 186     {
 187       sprintf(histname,"paramspace_%d",i);
 188       fParamSpace[i] = new AliL3Histogram(histname,"",nxbin,xmin,xmax,nybin,ymin,ymax);
 189     }
 190
 191 #ifdef do_mc
 192   if(fDoMC)
 193     {
 194       AliL3Histogram *hist = fParamSpace[0];
 195       Int_t ncellsx = (hist->GetNbinsX()+3)/2;
 196       Int_t ncellsy = (hist->GetNbinsY()+3)/2;
 197       Int_t ncells = ncellsx*ncellsy;
 198       if(!fgTrackID)
 199         {
 200           LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 201             <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(AliL3TrackIndex)<<" bytes to fgTrackID"<<ENDLOG;
 202           fgTrackID = new AliL3TrackIndex*[GetNEtaSegments()];
 203           for(Int_t i=0; i<GetNEtaSegments(); i++)
 204             fgTrackID[i] = new AliL3TrackIndex[ncells];
 205         }
 206     }
 207 #endif
 208   AliL3Histogram *hist = fParamSpace[0];
 209   Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
 210   if(!fgRowCount)
 211     {
 212       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 213         <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fgRowCount"<<ENDLOG;
 214       fgRowCount = new UChar_t*[GetNEtaSegments()];
 215       for(Int_t i=0; i<GetNEtaSegments(); i++)
 216         fgRowCount[i] = new UChar_t[ncells];
 217     }
 218   if(!fgGapCount)
 219     {
 220       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 221         <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fgGapCount"<<ENDLOG;
 222       fgGapCount = new UChar_t*[GetNEtaSegments()];
 223       for(Int_t i=0; i<GetNEtaSegments(); i++)
 224         fgGapCount[i] = new UChar_t[ncells];
 225     }
 226   if(!fgCurrentRowCount)
 227     {
 228       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 229         <<"Transformer: Allocating "<<GetNEtaSegments()*ncells*sizeof(UChar_t)<<" bytes to fgCurrentRowCount"<<ENDLOG;
 230       fgCurrentRowCount = new UChar_t*[GetNEtaSegments()];
 231       for(Int_t i=0; i<GetNEtaSegments(); i++)
 232         fgCurrentRowCount[i] = new UChar_t[ncells];
 233     }
 234
 235   if(!fgTrackNRows)
 236     {
 237       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 238         <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fgTrackNRows"<<ENDLOG;
 239       fgTrackNRows = new UChar_t[ncells];
 240       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 241         <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fgTrackFirstRow"<<ENDLOG;
 242       fgTrackFirstRow = new UChar_t[ncells];
 243       LOG(AliL3Log::kInformational,"AliL3HoughTransformerRow::CreateHistograms()","")
 244         <<"Transformer: Allocating "<<ncells*sizeof(UChar_t)<<" bytes to fgTrackLastRow"<<ENDLOG;
 245       fgTrackLastRow = new UChar_t[ncells];
 246
 247       AliL3HoughTrack track;
 248       Int_t xmin = hist->GetFirstXbin();
 249       Int_t xmax = hist->GetLastXbin();
 250       Int_t ymin = hist->GetFirstYbin();
 251       Int_t ymax = hist->GetLastYbin();
 252       Int_t nxbins = hist->GetNbinsX()+2;
 253       for(Int_t ybin=ymin; ybin<=ymax; ybin++)
 254         {
 255           for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 256             {
 257               //cvetan: we get strange warning on gcc-2.95
 258               //warning: large integer implicitly truncated to unsigned type
 259               fgTrackNRows[xbin + ybin*nxbins] = 99999;
 260               for(Int_t deltay = -1; deltay <= 1; deltay += 2) {
 261                 for(Int_t deltax = -1; deltax <= 1; deltax += 2) {
 262
 263                   Float_t xtrack = hist->GetPreciseBinCenterX((Float_t)xbin+0.5*(Float_t)deltax);
 264                   Float_t ytrack = hist->GetPreciseBinCenterY((Float_t)ybin+0.5*(Float_t)deltay);
 265
 266                   Float_t psi = atan((xtrack-ytrack)/(fgBeta1-fgBeta2));
 267                   Float_t kappa = 2.0*(xtrack*cos(psi)-fgBeta1*sin(psi));
 268                   track.SetTrackParameters(kappa,psi,1);
 269                   Bool_t firstrow = kFALSE;
 270                   UInt_t maxfirstrow = 0;
 271                   UInt_t maxlastrow = 0;
 272                   UInt_t curfirstrow = 0;
 273                   UInt_t curlastrow = 0;
 274                   for(Int_t j=AliL3Transform::GetFirstRow(0); j<=AliL3Transform::GetLastRow(5); j++)
 275                     {
 276                       Float_t hit[3];
 277                       if(!track.GetCrossingPoint(j,hit)) continue;
 278                       AliL3Transform::LocHLT2Raw(hit,0,j);
 279                       if(hit[1]>=0 && hit[1]<AliL3Transform::GetNPads(j))
 280                         {
 281                           if(!firstrow) {
 282                             curfirstrow = j;
 283                             firstrow = kTRUE;
 284                           }
 285                           curlastrow = j;
 286                         }
 287                       else {
 288                         if(firstrow) {
 289                           firstrow = kFALSE;
 290                           if((curlastrow-curfirstrow) >= (maxlastrow-maxfirstrow)) {
 291                             maxfirstrow = curfirstrow;
 292                             maxlastrow = curlastrow;
 293                           }
 294                         }
 295                       }
 296                     }
 297                   if((curlastrow-curfirstrow) >= (maxlastrow-maxfirstrow)) {
 298                     maxfirstrow = curfirstrow;
 299                     maxlastrow = curlastrow;
 300                   }
 301                   if((maxlastrow-maxfirstrow) < fgTrackNRows[xbin + ybin*nxbins]) {
 302                     fgTrackNRows[xbin + ybin*nxbins] = maxlastrow-maxfirstrow;
 303                     fgTrackFirstRow[xbin + ybin*nxbins] = maxfirstrow;
 304                     fgTrackLastRow[xbin + ybin*nxbins] = maxlastrow;
 305                   }
 306                 }
 307               }
 308               //              cout<<" fgTrackNRows "<<xbin<<" "<<ybin<<" "<<(Int_t)fgTrackNRows[xbin + ybin*nxbins]<<" "<<endl;
 309             }
 310         }
 311     }
 312
 313   //create lookup table for z of the digits
 314   Int_t ntimebins = AliL3Transform::GetNTimeBins();
 315   fLUTforwardZ = new Float_t[ntimebins];
 316   fLUTforwardZ2 = new Float_t[ntimebins];
 317   fLUTbackwardZ = new Float_t[ntimebins];
 318   fLUTbackwardZ2 = new Float_t[ntimebins];
 319   for(Int_t i=0; i<ntimebins; i++){
 320     Float_t z;
 321     z=AliL3Transform::GetZFast(0,i,GetZVertex());
 322     fLUTforwardZ[i]=z;
 323     fLUTforwardZ2[i]=z*z;
 324     z=AliL3Transform::GetZFast(18,i,GetZVertex());
 325     fLUTbackwardZ[i]=z;
 326     fLUTbackwardZ2[i]=z*z;
 327   }
 328
 329 }
 330
 331 void AliL3HoughTransformerRow::Reset()
 332 {
 333   //Reset all the histograms. Should be done when processing new slice
 334
 335   if(!fParamSpace)
 336     {
 337       LOG(AliL3Log::kWarning,"AliL3HoughTransformer::Reset","Histograms")
 338         <<"No histograms to reset"<<ENDLOG;
 339       return;
 340     }
 341
 342   for(Int_t i=0; i<GetNEtaSegments(); i++)
 343     fParamSpace[i]->Reset();
 344
 345 #ifdef do_mc
 346   if(fDoMC)
 347     {
 348       AliL3Histogram *hist = fParamSpace[0];
 349       Int_t ncellsx = (hist->GetNbinsX()+3)/2;
 350       Int_t ncellsy = (hist->GetNbinsY()+3)/2;
 351       Int_t ncells = ncellsx*ncellsy;
 352       for(Int_t i=0; i<GetNEtaSegments(); i++)
 353         memset(fgTrackID[i],0,ncells*sizeof(AliL3TrackIndex));
 354     }
 355 #endif
 356   AliL3Histogram *hist = fParamSpace[0];
 357   Int_t ncells = (hist->GetNbinsX()+2)*(hist->GetNbinsY()+2);
 358   for(Int_t i=0; i<GetNEtaSegments(); i++)
 359     {
 360       memset(fgRowCount[i],0,ncells*sizeof(UChar_t));
 361       memset(fgGapCount[i],1,ncells*sizeof(UChar_t));
 362       //      memset(fgCurrentRowCount[i],160,ncells*sizeof(UChar_t));
 363       memcpy(fgCurrentRowCount[i],fgTrackFirstRow,ncells*sizeof(UChar_t));
 364     }
 365 }
 366
 367 Int_t AliL3HoughTransformerRow::GetEtaIndex(Double_t eta) const
 368 {
 369   //Return the histogram index of the corresponding eta.
 370
 371   Double_t etaslice = (GetEtaMax() - GetEtaMin())/GetNEtaSegments();
 372   Double_t index = (eta-GetEtaMin())/etaslice;
 373   return (Int_t)index;
 374 }
 375
 376 inline AliL3Histogram *AliL3HoughTransformerRow::GetHistogram(Int_t etaindex)
 377 {
 378   // Return a pointer to the histogram which contains etaindex eta slice
 379   if(!fParamSpace || etaindex >= GetNEtaSegments() || etaindex < 0)
 380     return 0;
 381   if(!fParamSpace[etaindex])
 382     return 0;
 383   return fParamSpace[etaindex];
 384 }
 385
 386 Double_t AliL3HoughTransformerRow::GetEta(Int_t etaindex,Int_t /*slice*/) const
 387 {
 388   // Return eta calculated in the middle of the eta slice
 389   Double_t etaslice = (GetEtaMax()-GetEtaMin())/GetNEtaSegments();
 390   Double_t eta=0;
 391   eta=(Double_t)((etaindex+0.5)*etaslice);
 392   return eta;
 393 }
 394
 395 struct AliL3EtaRow {
 396   UChar_t fStartPad; //First pad in the cluster
 397   UChar_t fEndPad; //Last pad in the cluster
 398   Bool_t fIsFound; //Is the cluster already found
 399   Float_t fStartY; //Y position of the first pad in the cluster
 400 #ifdef do_mc
 401   Int_t fMcLabels[MaxTrack]; //Array to store mc labels inside cluster
 402 #endif
 403 };
 404
 405 void AliL3HoughTransformerRow::TransformCircle()
 406 {
 407   //Do the Hough Transform
 408   Float_t beta1 = fgBeta1;
 409   Float_t beta2 = fgBeta2;
 410   Float_t beta1minusbeta2 = fgBeta1 - fgBeta2;
 411
 412   Int_t netasegments = GetNEtaSegments();
 413   Double_t etamin = GetEtaMin();
 414   Double_t etaslice = (GetEtaMax() - etamin)/netasegments;
 415
 416   Int_t lowerthreshold = GetLowerThreshold();
 417
 418   //Assumes that all the etaslice histos are the same!
 419   AliL3Histogram *h = fParamSpace[0];
 420   Float_t ymin = h->GetYmin();
 421   //Float_t y_max = h->GetYmax();
 422   Float_t histbin = h->GetBinWidthY();
 423   Int_t firstbin = h->GetFirstYbin();
 424   Int_t lastbin = h->GetLastYbin();
 425   Float_t xmin = h->GetXmin();
 426   Float_t xmax = h->GetXmax();
 427   Float_t xbin = (xmax-xmin)/h->GetNbinsX();
 428   Int_t firstbinx = h->GetFirstXbin()-1;
 429   Int_t lastbinx = h->GetLastXbin()+1;
 430   Int_t nbinx = h->GetNbinsX()+2;
 431
 432   UChar_t lastpad;
 433   Int_t lastetaindex;
 434   AliL3EtaRow *etaclust = new AliL3EtaRow[netasegments];
 435
 436   AliL3DigitRowData *tempPt = GetDataPointer();
 437   if(!tempPt)
 438     {
 439       LOG(AliL3Log::kError,"AliL3HoughTransformer::TransformCircle","Data")
 440         <<"No input data "<<ENDLOG;
 441       return;
 442     }
 443
 444   Int_t ipatch = GetPatch();
 445   Double_t padpitch = AliL3Transform::GetPadPitchWidth(ipatch);
 446   Int_t islice = GetSlice();
 447   Float_t *fLUTz;
 448   Float_t *fLUTz2;
 449   if(islice < 18) {
 450     fLUTz = fLUTforwardZ;
 451     fLUTz2 = fLUTforwardZ2;
 452   }
 453   else {
 454     fLUTz = fLUTbackwardZ;
 455     fLUTz2 = fLUTbackwardZ2;
 456   }
 457
 458   //Loop over the padrows:
 459   for(UChar_t i=AliL3Transform::GetFirstRow(ipatch); i<=AliL3Transform::GetLastRow(ipatch); i++)
 460     {
 461       Int_t npads = AliL3Transform::GetNPads((Int_t)i)-1;
 462
 463       lastpad = 255;
 464       //Flush eta clusters array
 465       memset(etaclust,0,netasegments*sizeof(AliL3EtaRow));
 466
 467       Float_t x = AliL3Transform::Row2X((Int_t)i);
 468       Float_t x2 = x*x;
 469       Float_t y=0,r2=0;
 470
 471       //Get the data on this padrow:
 472       AliL3DigitData *digPt = tempPt->fDigitData;
 473       if((Int_t)i != (Int_t)tempPt->fRow)
 474         {
 475           LOG(AliL3Log::kError,"AliL3HoughTransformerRow::TransformCircle","Data")
 476             <<"AliL3HoughTransform::TransformCircle : Mismatching padrow numbering "<<(Int_t)i<<" "<<(Int_t)tempPt->fRow<<ENDLOG;
 477           continue;
 478         }
 479       //      cout<<" Starting row "<<i<<endl;
 480       //Loop over the data on this padrow:
 481       for(UInt_t j=0; j<tempPt->fNDigit; j++)
 482         {
 483           UShort_t charge = digPt[j].fCharge;
 484           if((Int_t)charge <= lowerthreshold)
 485             continue;
 486           UChar_t pad = digPt[j].fPad;
 487           UShort_t time = digPt[j].fTime;
 488
 489           if(pad != lastpad)
 490             {
 491               y = (pad-0.5*npads)*padpitch;
 492               Float_t y2 = y*y;
 493               r2 = x2 + y2;
 494               lastetaindex = -1;
 495             }
 496
 497           //Transform data to local cartesian coordinates:
 498           Float_t z = fLUTz[(Int_t)time];
 499           Float_t z2 = fLUTz2[(Int_t)time];
 500           //Calculate the eta:
 501           Double_t r = sqrt(r2+z2);
 502           Double_t eta = 0.5 * log((r+z)/(r-z));
 503           //Get the corresponding index, which determines which histogram to fill:
 504           Int_t etaindex = (Int_t)((eta-etamin)/etaslice);
 505
 506 #ifndef do_mc
 507           if(etaindex == lastetaindex) continue;
 508 #endif
 509           //      cout<<" Digit at patch "<<ipatch<<" row "<<i<<" pad "<<(Int_t)pad<<" time "<<time<<" etaslice "<<etaindex<<endl;
 510
 511           if(etaindex < 0 || etaindex >= netasegments)
 512             continue;
 513
 514           if(!etaclust[etaindex].fIsFound)
 515             {
 516               etaclust[etaindex].fStartPad = pad;
 517               etaclust[etaindex].fEndPad = pad;
 518               etaclust[etaindex].fStartY = y - padpitch/2.0;
 519               etaclust[etaindex].fIsFound = 1;
 520 #ifdef do_mc
 521               if(fDoMC)
 522                 {
 523                   for(Int_t t=0; t<3; t++)
 524                     {
 525                       Int_t label = digPt[j].fTrackID[t];
 526                       if(label < 0) break;
 527                       UInt_t c;
 528                       for(c=0; c<(MaxTrack-1); c++)
 529                         if(etaclust[etaindex].fMcLabels[c] == label || etaclust[etaindex].fMcLabels[c] == 0)
 530                           break;
 531                       //                      if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
 532                       etaclust[etaindex].fMcLabels[c] = label;
 533                     }
 534                 }
 535 #endif
 536               continue;
 537             }
 538           else
 539             {
 540               if(pad <= (etaclust[etaindex].fEndPad + 1))
 541                 {
 542                   etaclust[etaindex].fEndPad = pad;
 543 #ifdef do_mc
 544                   if(fDoMC)
 545                     {
 546                       for(Int_t t=0; t<3; t++)
 547                         {
 548                           Int_t label = digPt[j].fTrackID[t];
 549                           if(label < 0) break;
 550                           UInt_t c;
 551                           for(c=0; c<(MaxTrack-1); c++)
 552                             if(etaclust[etaindex].fMcLabels[c] == label || etaclust[etaindex].fMcLabels[c] == 0)
 553                               break;
 554                           //                      if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
 555                           etaclust[etaindex].fMcLabels[c] = label;
 556                         }
 557                     }
 558 #endif
 559                 }
 560               else
 561                 {
 562                   Bool_t fillcluster = kTRUE;
 563                   if(fillcluster) {
 564
 565                   UChar_t *nrows = fgRowCount[etaindex];
 566                   UChar_t *ngaps = fgGapCount[etaindex];
 567                   UChar_t *currentrow = fgCurrentRowCount[etaindex];
 568                   UChar_t *lastrow = fgTrackLastRow;
 569
 570                   //Do the transformation:
 571                   Float_t starty = etaclust[etaindex].fStartY;
 572                   if(etaclust[etaindex].fStartPad == 0)
 573                     starty -= 0.0*padpitch;
 574                   Float_t r1 = x2 + starty*starty;
 575                   Float_t xoverr1 = x/r1;
 576                   Float_t startyoverr1 = starty/r1;
 577                   Float_t endy = (etaclust[etaindex].fEndPad-0.5*(npads-1))*padpitch;
 578                   if(etaclust[etaindex].fEndPad == npads)
 579                     endy += 0.0*padpitch;
 580                   Float_t r2 = x2 + endy*endy;
 581                   Float_t xoverr2 = x/r2;
 582                   Float_t endyoverr2 = endy/r2;
 583                   Float_t a1 = beta1minusbeta2/(xoverr1-beta2);
 584                   Float_t b1 = (xoverr1-beta1)/(xoverr1-beta2);
 585                   Float_t a2 = beta1minusbeta2/(xoverr2-beta2);
 586                   Float_t b2 = (xoverr2-beta1)/(xoverr2-beta2);
 587
 588                   Float_t kappa1 = (a1*startyoverr1+b1*ymin-xmin)/xbin;
 589                   Float_t deltaalpha1 = b1*histbin/xbin;
 590                   if(b1<0)
 591                     kappa1 += deltaalpha1;
 592                   Float_t kappa2 = (a2*endyoverr2+b2*ymin-xmin)/xbin;
 593                   Float_t deltaalpha2 = b2*histbin/xbin;
 594                   if(b2>=0)
 595                     kappa2 += deltaalpha2;
 596
 597                   //Fill the histogram along the phirange
 598                   for(Int_t b=firstbin; b<=lastbin; b++, kappa1 += deltaalpha1, kappa2 += deltaalpha2)
 599                     {
 600                       Int_t binx1 = 1 + (Int_t)kappa1;
 601                       if(binx1>lastbinx) continue;
 602                       if(binx1<firstbinx) binx1 = firstbinx;
 603                       Int_t binx2 = 1 + (Int_t)kappa2;
 604                       if(binx2<firstbinx) continue;
 605                       if(binx2>lastbinx) binx2 = lastbinx;
 606 #ifdef do_mc
 607                       if(binx2<binx1) {
 608                         LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::TransformCircle()","")
 609                           <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<x<<" "<<starty<<" "<<endy<<ENDLOG;
 610                       }
 611 #endif
 612                       Int_t tempbin = b*nbinx;
 613                       UChar_t *nrows2 = nrows + tempbin;
 614                       UChar_t *ngaps2 = ngaps + tempbin;
 615                       UChar_t *currentrow2 = currentrow + tempbin;
 616                       UChar_t *lastrow2 = lastrow + tempbin;
 617                       for(Int_t bin=binx1;bin<=binx2;bin++)
 618                         {
 619                           if(ngaps2[bin] < MAX_N_GAPS) {
 620                             if(i > (currentrow2[bin] + MAX_GAP_SIZE) && i < lastrow2[bin]) {
 621                               ngaps2[bin] = MAX_N_GAPS;
 622                               continue;
 623                             }
 624                             if(i > currentrow2[bin] && i < lastrow2[bin])
 625                               {
 626                                 nrows2[bin]++;
 627                                 if(i > (currentrow2[bin] + 1))
 628                                   ngaps2[bin]++;
 629                                 currentrow2[bin]=i;
 630                               }
 631 #ifdef do_mc
 632                             if(fDoMC)
 633                               {
 634                                 for(UInt_t t=0;t<(MaxTrack-1); t++)
 635                                   {
 636                                     Int_t label = etaclust[etaindex].fMcLabels[t];
 637                                     if(label == 0) break;
 638                                     UInt_t c;
 639                                     Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2)) + (Int_t)(bin/2);
 640                                     for(c=0; c<(MaxTrack-1); c++)
 641                                       if(fgTrackID[etaindex][tempbin2].fLabel[c] == label || fgTrackID[etaindex][tempbin2].fNHits[c] == 0)
 642                                         break;
 643                                     //                              if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
 644                                     fgTrackID[etaindex][tempbin2].fLabel[c] = label;
 645                                     if(fgTrackID[etaindex][tempbin2].fCurrentRow[c] != i) {
 646                                       fgTrackID[etaindex][tempbin2].fNHits[c]++;
 647                                       fgTrackID[etaindex][tempbin2].fCurrentRow[c] = i;
 648                                     }
 649                                   }
 650                               }
 651 #endif
 652                           }
 653                         }
 654                     }
 655                   //End of the transformation
 656
 657                   }
 658
 659                   etaclust[etaindex].fStartPad = pad;
 660                   etaclust[etaindex].fEndPad = pad;
 661                   etaclust[etaindex].fStartY = y - padpitch/2.0;
 662
 663 #ifdef do_mc
 664                   if(fDoMC)
 665                     {
 666                       memset(etaclust[etaindex].fMcLabels,0,MaxTrack);
 667                       for(Int_t t=0; t<3; t++)
 668                         {
 669                           Int_t label = digPt[j].fTrackID[t];
 670                           if(label < 0) break;
 671                           UInt_t c;
 672                           for(c=0; c<(MaxTrack-1); c++)
 673                             if(etaclust[etaindex].fMcLabels[c] == label || etaclust[etaindex].fMcLabels[c] == 0)
 674                               break;
 675                           //                      if(c == MaxTrack-1) cerr<<"AliL3HoughTransformer::TransformCircle : Cluster array reached maximum!! "<<c<<endl;
 676                           etaclust[etaindex].fMcLabels[c] = label;
 677                         }
 678                     }
 679 #endif
 680                 }
 681             }
 682           lastpad = pad;
 683           lastetaindex = etaindex;
 684         }
 685       //Write remaining clusters
 686       for(Int_t etaindex = 0;etaindex < netasegments;etaindex++)
 687         {
 688           //Check for empty row
 689           if((etaclust[etaindex].fStartPad == 0) && (etaclust[etaindex].fEndPad == 0)) continue;
 690
 691           UChar_t *nrows = fgRowCount[etaindex];
 692           UChar_t *ngaps = fgGapCount[etaindex];
 693           UChar_t *currentrow = fgCurrentRowCount[etaindex];
 694           UChar_t *lastrow = fgTrackLastRow;
 695
 696           //Do the transformation:
 697           Float_t starty = etaclust[etaindex].fStartY;
 698           if(etaclust[etaindex].fStartPad == 0)
 699             starty -= 0.0*padpitch;
 700           Float_t r1 = x2 + starty*starty;
 701           Float_t xoverr1 = x/r1;
 702           Float_t startyoverr1 = starty/r1;
 703           Float_t endy = (etaclust[etaindex].fEndPad-0.5*(npads-1))*padpitch;
 704           if(etaclust[etaindex].fEndPad == npads)
 705             endy += 0.0*padpitch;
 706           Float_t r2 = x2 + endy*endy;
 707           Float_t xoverr2 = x/r2;
 708           Float_t endyoverr2 = endy/r2;
 709           Float_t a1 = beta1minusbeta2/(xoverr1-beta2);
 710           Float_t b1 = (xoverr1-beta1)/(xoverr1-beta2);
 711           Float_t a2 = beta1minusbeta2/(xoverr2-beta2);
 712           Float_t b2 = (xoverr2-beta1)/(xoverr2-beta2);
 713
 714           Float_t kappa1 = (a1*startyoverr1+b1*ymin-xmin)/xbin;
 715           Float_t deltaalpha1 = b1*histbin/xbin;
 716           if(b1<0)
 717             kappa1 += deltaalpha1;
 718           Float_t kappa2 = (a2*endyoverr2+b2*ymin-xmin)/xbin;
 719           Float_t deltaalpha2 = b2*histbin/xbin;
 720           if(b2>=0)
 721             kappa2 += deltaalpha2;
 722
 723           //Fill the histogram along the phirange
 724           for(Int_t b=firstbin; b<=lastbin; b++, kappa1 += deltaalpha1, kappa2 += deltaalpha2)
 725             {
 726               Int_t binx1 = 1 + (Int_t)kappa1;
 727               if(binx1>lastbinx) continue;
 728               if(binx1<firstbinx) binx1 = firstbinx;
 729               Int_t binx2 = 1 + (Int_t)kappa2;
 730               if(binx2<firstbinx) continue;
 731               if(binx2>lastbinx) binx2 = lastbinx;
 732 #ifdef do_mc
 733               if(binx2<binx1) {
 734                 LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::TransformCircle()","")
 735                   <<"Wrong filling "<<binx1<<" "<<binx2<<" "<<i<<" "<<x<<" "<<starty<<" "<<endy<<ENDLOG;
 736               }
 737 #endif
 738               Int_t tempbin = b*nbinx;
 739               UChar_t *nrows2 = nrows + tempbin;
 740               UChar_t *ngaps2 = ngaps + tempbin;
 741               UChar_t *currentrow2 = currentrow + tempbin;
 742               UChar_t *lastrow2 = lastrow + tempbin;
 743               for(Int_t bin=binx1;bin<=binx2;bin++)
 744                 {
 745                   if(ngaps2[bin] < MAX_N_GAPS) {
 746                     if(i > (currentrow2[bin] + MAX_GAP_SIZE) && i < lastrow2[bin]) {
 747                       ngaps2[bin] = MAX_N_GAPS;
 748                       continue;
 749                     }
 750                     if(i > currentrow2[bin] && i < lastrow2[bin])
 751                       {
 752                         nrows2[bin]++;
 753                         if(i > (currentrow2[bin] + 1))
 754                           ngaps2[bin]++;
 755                         currentrow2[bin]=i;
 756                       }
 757 #ifdef do_mc
 758                     if(fDoMC)
 759                       {
 760                         for(UInt_t t=0;t<(MaxTrack-1); t++)
 761                           {
 762                             Int_t label = etaclust[etaindex].fMcLabels[t];
 763                             if(label == 0) break;
 764                             UInt_t c;
 765                             Int_t tempbin2 = ((Int_t)(b/2))*((Int_t)((nbinx+1)/2)) + (Int_t)(bin/2);
 766                             for(c=0; c<(MaxTrack-1); c++)
 767                               if(fgTrackID[etaindex][tempbin2].fLabel[c] == label || fgTrackID[etaindex][tempbin2].fNHits[c] == 0)
 768                                 break;
 769                             //                      if(c == MaxTrack-1) cout<<"AliL3HoughTransformer::TransformCircle : Array reached maximum!! "<<c<<endl;
 770                             fgTrackID[etaindex][tempbin2].fLabel[c] = label;
 771                             if(fgTrackID[etaindex][tempbin2].fCurrentRow[c] != i) {
 772                               fgTrackID[etaindex][tempbin2].fNHits[c]++;
 773                               fgTrackID[etaindex][tempbin2].fCurrentRow[c] = i;
 774                             }
 775                           }
 776                       }
 777 #endif
 778                   }
 779                 }
 780             }
 781           //End of the transformation
 782
 783         }
 784
 785       //Move the data pointer to the next padrow:
 786       AliL3MemHandler::UpdateRowPointer(tempPt);
 787     }
 788
 789   delete [] etaclust;
 790 }
 791
 792 Int_t AliL3HoughTransformerRow::GetTrackID(Int_t etaindex,Double_t kappa,Double_t psi)
 793 {
 794   // Returns the MC label for a given peak found in the Hough space
 795   if(!fDoMC)
 796     {
 797       LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID","Data")
 798         <<"Flag switched off"<<ENDLOG;
 799       return -1;
 800     }
 801
 802 #ifdef do_mc
 803   if(etaindex < 0 || etaindex > GetNEtaSegments())
 804     {
 805       LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID","Data")
 806         <<"Wrong etaindex "<<etaindex<<ENDLOG;
 807       return -1;
 808     }
 809   AliL3Histogram *hist = fParamSpace[etaindex];
 810   Int_t bin = hist->FindLabelBin(kappa,psi);
 811   if(bin==-1) {
 812     LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID()","")
 813       <<"Track candidate outside Hough space boundaries: "<<kappa<<" "<<psi<<ENDLOG;
 814     return -1;
 815   }
 816   Int_t label=-1;
 817   Int_t max=0;
 818   for(UInt_t i=0; i<(MaxTrack-1); i++)
 819     {
 820       Int_t nhits=fgTrackID[etaindex][bin].fNHits[i];
 821       if(nhits == 0) break;
 822       if(nhits > max)
 823         {
 824           max = nhits;
 825           label = fgTrackID[etaindex][bin].fLabel[i];
 826         }
 827     }
 828   Int_t label2=-1;
 829   Int_t max2=0;
 830   for(UInt_t i=0; i<(MaxTrack-1); i++)
 831     {
 832       Int_t nhits=fgTrackID[etaindex][bin].fNHits[i];
 833       if(nhits == 0) break;
 834       if(nhits > max2)
 835         {
 836           if(fgTrackID[etaindex][bin].fLabel[i]!=label) {
 837             max2 = nhits;
 838             label2 = fgTrackID[etaindex][bin].fLabel[i];
 839           }
 840         }
 841     }
 842   if(max2 !=0 ) {
 843     LOG(AliL3Log::kDebug,"AliL3HoughTransformerRow::GetTrackID()","")
 844       <<" TrackID"<<" label "<<label<<" max "<<max<<" label2 "<<label2<<" max2 "<<max2<<" "<<(Float_t)max2/(Float_t)max<<" "<<fgTrackID[etaindex][bin].fLabel[MaxTrack-1]<<" "<<(Int_t)fgTrackID[etaindex][bin].fNHits[MaxTrack-1]<<ENDLOG;
 845   }
 846   return label;
 847 #endif
 848   LOG(AliL3Log::kWarning,"AliL3HoughTransformerRow::GetTrackID()","")
 849     <<"Compile with do_mc flag!"<<ENDLOG;
 850   return -1;
 851 }
 852
 853 UChar_t *AliL3HoughTransformerRow::GetRowCount(Int_t etaindex)
 854 {
 855   return fgRowCount[etaindex];
 856 }
 857
 858 UChar_t *AliL3HoughTransformerRow::GetGapCount(Int_t etaindex)
 859 {
 860   return fgGapCount[etaindex];
 861 }
 862
 863 UChar_t *AliL3HoughTransformerRow::GetCurrentRowCount(Int_t etaindex)
 864 {
 865   return fgCurrentRowCount[etaindex];
 866 }
 867
 868 UChar_t *AliL3HoughTransformerRow::GetTrackNRows()
 869 {
 870   return fgTrackNRows;
 871 }
 872
 873
 874 UChar_t *AliL3HoughTransformerRow::GetTrackFirstRow()
 875 {
 876   return fgTrackFirstRow;
 877 }
 878
 879 UChar_t *AliL3HoughTransformerRow::GetTrackLastRow()
 880 {
 881   return fgTrackLastRow;
 882 }