HLT/hough/AliL3HoughTest.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 #include "AliL3HoughTest.h"
   8 #include "AliL3ModelTrack.h"
   9 #include "AliL3Transform.h"
  10 #include "AliL3Histogram.h"
  11 #include "AliL3TrackArray.h"
  12 #include "AliL3HoughTrack.h"
  13
  14 #include <TRandom.h>
  15 #include <TMath.h>
  16 #include <TH2.h>
  17 #include <TH3.h>
  18 #include <TAxis.h>
  19
  20 #if __GNUC__ == 3
  21 using namespace std;
  22 #endif
  23
  24 //_____________________________________________________________
  25 // AliL3HoughTest
  26
  27
  28 ClassImp(AliL3HoughTest)
  29
  30 AliL3HoughTest::AliL3HoughTest()
  31 {
  32   //ctor
  33   fData=0;
  34 }
  35
  36
  37 AliL3HoughTest::~AliL3HoughTest()
  38 {
  39   //dtor
  40   if(fData)
  41     delete [] fData;
  42 }
  43
  44 Bool_t AliL3HoughTest::GenerateTrackData(Double_t pt,Double_t psi,Double_t tgl,Int_t sign,Int_t patch,Int_t minhits)
  45 {
  46   //Generate digits according to given track parameters?
  47   fCurrentPatch=patch;
  48   if(fData)
  49     delete fData;
  50   fData = new AliL3SimData[AliL3Transform::GetNRows(patch)];
  51   memset(fData,0,AliL3Transform::GetNRows(patch)*sizeof(AliL3SimData));
  52
  53   AliL3ModelTrack *track = new AliL3ModelTrack();
  54   track->Init(0,patch);
  55   track->SetPt(pt);
  56   track->SetPsi(psi);
  57   track->SetTgl(tgl);
  58   track->SetCharge(sign);
  59   track->SetFirstPoint(0,0,0);
  60   track->CalculateHelix();
  61
  62   Int_t temp[200];
  63   //  Int_t temp2[AliL3Transform::GetNTimeBins()];
  64   Int_t * temp2 = new Int_t[AliL3Transform::GetNTimeBins()];
  65   Int_t entries=100;
  66   Int_t clustercharge=100;
  67   Int_t hitcounter=0;
  68   for(Int_t i=AliL3Transform::GetFirstRow(patch); i<=AliL3Transform::GetLastRow(patch); i++)
  69     {
  70
  71       Float_t xyz[3];
  72       if(!track->GetCrossingPoint(i,xyz))
  73         continue;
  74
  75       Int_t rowindex = i - AliL3Transform::GetFirstRow(patch);
  76       Int_t sector,row;
  77       AliL3Transform::Slice2Sector(0,i,sector,row);
  78       AliL3Transform::Local2Raw(xyz,sector,row);
  79
  80       if(xyz[1] < 0 || xyz[1] >= AliL3Transform::GetNPads(i) || xyz[2] < 0 || xyz[2] >= AliL3Transform::GetNTimeBins())
  81         continue;
  82       hitcounter++;
  83       track->SetPadHit(i,xyz[1]);
  84       track->SetTimeHit(i,xyz[2]);
  85
  86       Double_t beta = track->GetCrossingAngle(i);
  87       track->SetCrossingAngleLUT(i,beta);
  88       track->CalculateClusterWidths(i);
  89
  90       memset(temp,0,200*sizeof(Int_t));
  91       memset(temp2,0,AliL3Transform::GetNTimeBins()*sizeof(Int_t));
  92       Double_t xysigma = sqrt(track->GetParSigmaY2(i));
  93       Double_t zsigma = sqrt(track->GetParSigmaZ2(i));
  94
  95       Int_t minpad=200,j;
  96       Int_t mintime = 1000;
  97       for(j=0; j<entries; j++)
  98         {
  99           Int_t pad = TMath::Nint(gRandom->Gaus(xyz[1],xysigma));
 100           Int_t time = TMath::Nint(gRandom->Gaus(xyz[2],zsigma));
 101           if(pad < 0 || pad >= AliL3Transform::GetNPads(i) || time < 0 || time >= AliL3Transform::GetNTimeBins())
 102             continue;
 103           temp[pad]++;
 104           temp2[time]++;
 105           if(pad < minpad)
 106             minpad=pad;
 107           if(time < mintime)
 108             mintime=time;
 109         }
 110       Int_t npads=0;
 111       for(j=0; j<200; j++)
 112         {
 113           if(temp[j]==0) continue;
 114
 115           Int_t index = j - minpad;
 116
 117           if(index < 0 || index >= 10)
 118             {
 119               cerr<<"AliL3HoughTest::GenerateTrackData : Wrong index "<<index<<endl;
 120               exit(5);
 121             }
 122           npads++;
 123           Int_t seqcharge = clustercharge*temp[j]/entries;
 124           Int_t ntimes=0;
 125           for(Int_t k=0; k<AliL3Transform::GetNTimeBins(); k++)
 126             {
 127               if(temp2[k]==0) continue;
 128               Int_t tindex = k - mintime;
 129               if(tindex < 0 || tindex >= 10)
 130                 {
 131                   cerr<<"AliL3HoughTest::GenerateTrackData : Wrong timeindex "<<tindex<<" "<<k<<" "<<mintime<<endl;
 132                   exit(5);
 133                 }
 134               Int_t charge = seqcharge*temp2[k]/entries;
 135               if(charge < 3 )
 136                 continue;
 137               if(charge > 1023)
 138                 charge=1023;
 139               //cout<<"row "<<i<<" pad "<<j<<" time "<<k<<" charge "<<charge<<endl;
 140               ntimes++;
 141               fData[rowindex].fPads[index][tindex]=charge;
 142             }
 143         }
 144       fData[rowindex].fMinpad=minpad;
 145       fData[rowindex].fMintime=mintime;
 146       fData[rowindex].fNPads=npads;
 147     }
 148   delete track;
 149   delete [] temp2;
 150   if(hitcounter < minhits)
 151     return kFALSE;
 152   return kTRUE;
 153 }
 154
 155 void AliL3HoughTest::Transform2Circle(AliL3Histogram *hist)
 156 {
 157   //Hough trasnform
 158   if(!fData)
 159     {
 160       cerr<<"AliL3HoughTest::Transform : No data"<<endl;
 161       return;
 162     }
 163   Float_t r,phi,phi0,kappa,xyz[3];
 164   Int_t pad,time,charge,sector,row;
 165   for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
 166     {
 167       Int_t rowindex = i - AliL3Transform::GetFirstRow(fCurrentPatch);
 168       AliL3Transform::Slice2Sector(0,i,sector,row);
 169       for(Int_t j=0; j<fData[rowindex].fNPads; j++)
 170         {
 171           pad = j + fData[rowindex].fMinpad;
 172           for(Int_t k=0; k<10; k++)
 173             {
 174               time = k + fData[rowindex].fMintime;
 175               charge = fData[rowindex].fPads[j][k];
 176               if(charge == 0) continue;
 177               AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
 178
 179               r = sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 180
 181               phi = AliL3Transform::GetPhi(xyz);
 182
 183               for(Int_t k=hist->GetFirstYbin(); k<=hist->GetLastYbin(); k++)
 184                 {
 185                   phi0 = hist->GetBinCenterY(k);
 186                   kappa = 2*sin(phi-phi0)/r;
 187                   hist->Fill(kappa,phi0,charge);
 188                 }
 189             }
 190         }
 191     }
 192 }
 193
 194 void AliL3HoughTest::Transform2CircleC(AliL3Histogram *hist)
 195 {
 196   //Hough transform
 197   if(!fData)
 198     {
 199       cerr<<"AliL3HoughTest::TransformC : No data"<<endl;
 200       return;
 201     }
 202   Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
 203   Float_t r1,r2,phi1,phi2,phi0,kappa,hit[3],hit2[3];
 204   for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
 205     {
 206       Int_t rowindex1 = i - AliL3Transform::GetFirstRow(fCurrentPatch);
 207       for(Int_t d1=0; d1<fData[rowindex1].fNPads; d1++)
 208         {
 209           pad1 = d1 + fData[rowindex1].fMinpad;
 210           for(Int_t j=0; j<10; j++)
 211             {
 212               time1 = j + fData[rowindex1].fMintime;
 213               charge1 = fData[rowindex1].fPads[d1][j];
 214               if(charge1==0) continue;
 215               AliL3Transform::Slice2Sector(0,i,sector,row);
 216               AliL3Transform::Raw2Local(hit,sector,row,pad1,time1);
 217               r1 = sqrt(hit[0]*hit[0]+hit[1]*hit[1]);
 218               phi1 = atan2(hit[1],hit[0]);
 219
 220               for(Int_t j=i+1; j<=AliL3Transform::GetLastRow(fCurrentPatch); j++)
 221                 {
 222                   Int_t rowindex2 = j - AliL3Transform::GetFirstRow(fCurrentPatch);
 223                   for(Int_t d2=0; d2<fData[rowindex2].fNPads; d2++)
 224                     {
 225                       pad2 = d2 + fData[rowindex2].fMinpad;
 226                       for(Int_t k=0; k<10; k++)
 227                         {
 228                           time2 = k + fData[rowindex2].fMintime;
 229                           charge2 = fData[rowindex2].fPads[d2][k];
 230                           if(charge2==0) continue;
 231                           AliL3Transform::Slice2Sector(0,j,sector,row);
 232                           AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
 233                           r2 = sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1]);
 234                           phi2 = atan2(hit2[1],hit2[0]);
 235                           phi0 = atan( (r2*sin(phi1) - r1*sin(phi2)) / (r2*cos(phi1) - r1*cos(phi2)) );
 236
 237                           kappa = 2*sin(phi1-phi0) / r1;
 238                           hist->Fill(kappa,phi0,charge1+charge2);
 239                         }
 240                     }
 241                 }
 242               return;
 243             }
 244         }
 245     }
 246 }
 247
 248 void AliL3HoughTest::Transform2CircleF(AliL3Histogram *hist)
 249 {
 250   //Fix one point in the middle of the tpc
 251
 252   if(!fData)
 253     {
 254       cerr<<"AliL3HoughTest::TransformF : No data"<<endl;
 255       return;
 256     }
 257   Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
 258   Float_t r1,r2,phi1,phi2,phi0,kappa,hit[3],hit2[3];
 259
 260   Int_t rowindex1 = 80 - AliL3Transform::GetFirstRow(fCurrentPatch);
 261   for(Int_t d1=1; d1<fData[rowindex1].fNPads; d1++)
 262     {
 263       pad1 = d1 + fData[rowindex1].fMinpad;
 264
 265       AliL3Transform::Slice2Sector(0,80,sector,row);
 266       AliL3Transform::Raw2Local(hit,sector,row,pad1,0);
 267       r1 = sqrt(hit[0]*hit[0]+hit[1]*hit[1]);
 268       phi1 = atan2(hit[1],hit[0]);
 269
 270       for(Int_t j=0; j<10; j++)
 271         {
 272           time1 = j + fData[rowindex1].fMintime;
 273           charge1 = fData[rowindex1].fPads[d1][j];
 274           if(charge1==0) continue;
 275
 276           for(Int_t j=0; j<=AliL3Transform::GetLastRow(fCurrentPatch); j++)
 277             {
 278               if(j==80) continue;
 279               Int_t rowindex2 = j - AliL3Transform::GetFirstRow(fCurrentPatch);
 280               for(Int_t d2=0; d2<fData[rowindex2].fNPads; d2++)
 281                 {
 282                   pad2 = d2 + fData[rowindex2].fMinpad;
 283                   for(Int_t k=0; k<10; k++)
 284                     {
 285                       time2 = k + fData[rowindex2].fMintime;
 286                       charge2 = fData[rowindex2].fPads[d2][k];
 287                       if(charge2==0) continue;
 288                       AliL3Transform::Slice2Sector(0,j,sector,row);
 289                       AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
 290                       r2 = sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1]);
 291                       phi2 = atan2(hit2[1],hit2[0]);
 292                       phi0 = atan( (r2*sin(phi1) - r1*sin(phi2)) / (r2*cos(phi1) - r1*cos(phi2)) );
 293
 294                       kappa = 2*sin(phi1-phi0) / r1;
 295                       hist->Fill(kappa,phi0,charge1+charge2);
 296                       //cout<<"Filling "<<kappa<<" psi "<<phi0<<" charge "<<charge1<<endl;
 297                     }
 298                 }
 299             }
 300         }
 301       return;
 302     }
 303 }
 304
 305 void AliL3HoughTest::Transform2Line(AliL3Histogram *hist,Int_t *rowrange)
 306 {
 307   //Hough transform
 308   if(!fData)
 309     {
 310       cerr<<"AliL3HoughTest::Transform2Line : No data"<<endl;
 311       return;
 312     }
 313
 314   Int_t pad,time,charge,sector,row;
 315   Float_t hit[3],theta,rho;
 316   for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
 317     {
 318       if(i<rowrange[0])
 319         continue;
 320       if(i>rowrange[1])
 321         break;
 322       Int_t rowindex = i - AliL3Transform::GetFirstRow(fCurrentPatch);
 323       for(Int_t d=0; d<fData[rowindex].fNPads; d++)
 324         {
 325           pad = d + fData[rowindex].fMinpad;
 326           for(Int_t j=0; j<10; j++)
 327             {
 328               time = j + fData[rowindex].fMintime;
 329               charge = fData[rowindex].fPads[d][j];
 330               if(charge==0) continue;
 331               AliL3Transform::Slice2Sector(0,i,sector,row);
 332               AliL3Transform::Raw2Local(hit,sector,row,pad,time);
 333
 334               hit[0] = hit[0] - AliL3Transform::Row2X(rowrange[0]);
 335
 336               for(Int_t xbin=hist->GetFirstXbin(); xbin<hist->GetLastXbin(); xbin++)
 337                 {
 338                   theta = hist->GetBinCenterX(xbin);
 339                   rho = hit[0]*cos(theta) + hit[1]*sin(theta);
 340                   hist->Fill(theta,rho,charge);
 341                 }
 342             }
 343         }
 344     }
 345 }
 346
 347 void AliL3HoughTest::Transform2LineC(AliL3Histogram *hist,Int_t *rowrange)
 348 {
 349   //Hough transform?
 350   if(!fData)
 351     {
 352       cerr<<"AliL3HoughTest::Transform2Line : No data"<<endl;
 353       return;
 354     }
 355
 356   Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
 357   Float_t theta,rho,hit[3],hit2[3];
 358   for(Int_t i=rowrange[0]; i<=rowrange[1]; i++)
 359     {
 360       Int_t rowindex1 = i - AliL3Transform::GetFirstRow(fCurrentPatch);
 361       for(Int_t d1=0; d1<fData[rowindex1].fNPads; d1++)
 362         {
 363           pad1 = d1 + fData[rowindex1].fMinpad;
 364           for(Int_t j=0; j<10; j++)
 365             {
 366               time1 = j + fData[rowindex1].fMintime;
 367               charge1 = fData[rowindex1].fPads[d1][j];
 368               if(charge1==0) continue;
 369               AliL3Transform::Slice2Sector(0,i,sector,row);
 370               AliL3Transform::Raw2Local(hit,sector,row,pad1,time1);
 371
 372               hit[0] = hit[0] - AliL3Transform::Row2X(rowrange[0]);
 373
 374               for(Int_t i2=i+1; i2<=rowrange[1]; i2++)
 375                 {
 376                   Int_t rowindex2 = i2 - AliL3Transform::GetFirstRow(fCurrentPatch);
 377                   for(Int_t d2=0; d2<fData[rowindex2].fNPads; d2++)
 378                     {
 379                       pad2 = d2 + fData[rowindex2].fMinpad;
 380                       for(Int_t k=0; k<10; k++)
 381                         {
 382                           time2 = k + fData[rowindex2].fMintime;
 383                           charge2 = fData[rowindex2].fPads[d2][k];
 384                           if(charge2==0) continue;
 385                           AliL3Transform::Slice2Sector(0,i2,sector,row);
 386                           AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
 387
 388                           hit2[0] = hit2[0] - AliL3Transform::Row2X(rowrange[0]);
 389
 390                           theta = atan2(hit2[0]-hit[0],hit[1]-hit2[1]);
 391                           rho = hit[0]*cos(theta)+hit[1]*sin(theta);
 392                           hist->Fill(theta,rho,charge1+charge2);
 393                         }
 394                     }
 395                 }
 396             }
 397         }
 398     }
 399 }
 400
 401 void AliL3HoughTest::FillImage(TH2 *hist,Int_t displayrow)
 402 {
 403   //Draw digits data
 404   if(!fData)
 405     {
 406       cerr<<"AliL3HoughTest::FillImage : No data to fill"<<endl;
 407       return;
 408     }
 409
 410   for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
 411     {
 412       Int_t rowindex = i - AliL3Transform::GetFirstRow(fCurrentPatch);
 413       if(displayrow >=0)
 414         if(i != displayrow) continue;
 415
 416       //cout<<"row "<<i<<" fNPads "<<fData[rowindex].fNPads<<endl;
 417       for(Int_t j=0; j<fData[rowindex].fNPads; j++)
 418         {
 419           Int_t pad = j + fData[rowindex].fMinpad;
 420           for(Int_t k=0; k<10; k++)
 421             {
 422               Int_t time = k + fData[rowindex].fMintime;
 423               Int_t charge = fData[rowindex].fPads[j][k];
 424               if(charge==0) continue;
 425               //cout<<i<<" "<<pad<<" "<<time<<" "<<charge<<endl;
 426               Float_t xyz[3];
 427               Int_t sector,row;
 428               AliL3Transform::Slice2Sector(0,i,sector,row);
 429               AliL3Transform::Raw2Local(xyz,sector,row,pad,time);
 430               if(displayrow>=0)
 431                 hist->Fill(pad,time,charge);
 432               else
 433                 hist->Fill(xyz[0],xyz[1],charge);
 434             }
 435         }
 436       if(displayrow>=0)
 437         break;
 438     }
 439 }
 440
 441 void AliL3HoughTest::Transform2Line3D(TH3 *hist,Int_t *rowrange,Float_t *phirange)
 442 {
 443   //HT in 3D?
 444   if(!fData)
 445     {
 446       cerr<<"AliL3HoughTest::Transform2Line : No data"<<endl;
 447       return;
 448     }
 449
 450   Int_t pad,time,charge,sector,row;
 451   Float_t hit[3],theta,rho,r,delta;
 452   for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
 453     {
 454
 455       if(i<rowrange[0])
 456         continue;
 457       if(i>rowrange[1])
 458         break;
 459
 460       Int_t rowindex = i - AliL3Transform::GetFirstRow(fCurrentPatch);
 461       for(Int_t d=0; d<fData[rowindex].fNPads; d++)
 462         {
 463           pad = d + fData[rowindex].fMinpad;
 464           for(Int_t j=0; j<10; j++)
 465             {
 466               time = j + fData[rowindex].fMintime;
 467               charge = fData[rowindex].fPads[d][j];
 468               if(charge==0) continue;
 469               AliL3Transform::Slice2Sector(0,i,sector,row);
 470               AliL3Transform::Raw2Local(hit,sector,row,pad,time);
 471
 472               Float_t phi = AliL3Transform::GetPhi(hit);
 473               if(phi < phirange[0] || phi > phirange[1])
 474                 continue;
 475
 476               hit[0] = hit[0] - AliL3Transform::Row2X(rowrange[0]);
 477               Float_t x = hit[0] + AliL3Transform::Row2X(rowrange[0]);
 478               r = sqrt(x*x + hit[1]*hit[1]);
 479               delta = atan(hit[2]/r);
 480
 481               for(Int_t xbin=hist->GetXaxis()->GetFirst(); xbin<=hist->GetXaxis()->GetLast(); xbin++)
 482                 {
 483                   theta = hist->GetXaxis()->GetBinCenter(xbin);
 484                   rho = hit[0]*cos(theta) + hit[1]*sin(theta);
 485                   hist->Fill(theta,rho,delta,charge);
 486                 }
 487             }
 488         }
 489     }
 490 }
 491
 492 void AliL3HoughTest::Transform2LineC3D(TH3 *hist,Int_t *rowrange)
 493 {
 494   //HT in 3D?
 495   if(!fData)
 496     {
 497       cerr<<"AliL3HoughTest::Transform2Line : No data"<<endl;
 498       return;
 499     }
 500
 501   Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
 502   Float_t theta,rho,hit[3],hit2[3],r1,r2,delta,delta1,delta2;
 503   for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
 504     {
 505       if(i<rowrange[0])
 506         continue;
 507       if(i>rowrange[1])
 508         break;
 509       Int_t rowindex1 = i - AliL3Transform::GetFirstRow(fCurrentPatch);
 510       for(Int_t d1=0; d1<fData[rowindex1].fNPads; d1++)
 511         {
 512           pad1 = d1 + fData[rowindex1].fMinpad;
 513           for(Int_t j=0; j<10; j++)
 514             {
 515               time1 = j + fData[rowindex1].fMintime;
 516               charge1 = fData[rowindex1].fPads[d1][j];
 517               if(charge1==0) continue;
 518               AliL3Transform::Slice2Sector(0,i,sector,row);
 519               AliL3Transform::Raw2Local(hit,sector,row,pad1,time1);
 520               r1 = sqrt(hit[0]*hit[0]+hit[1]*hit[1]);
 521               delta1 = atan(hit[2]/r1);
 522               hit[0] = hit[0] - AliL3Transform::Row2X(rowrange[0]);
 523
 524               for(Int_t i2=i+1; i2<=rowrange[1]; i2++)
 525                 {
 526                   Int_t rowindex2 = i2 - AliL3Transform::GetFirstRow(fCurrentPatch);
 527                   for(Int_t d2=0; d2<fData[rowindex2].fNPads; d2++)
 528                     {
 529                       pad2 = d2 + fData[rowindex2].fMinpad;
 530                       for(Int_t k=0; k<10; k++)
 531                         {
 532                           time2 = k + fData[rowindex2].fMintime;
 533                           charge2 = fData[rowindex2].fPads[d2][k];
 534                           if(charge2==0) continue;
 535                           AliL3Transform::Slice2Sector(0,i2,sector,row);
 536                           AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
 537                           r2 = sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1]);
 538                           delta2 = atan(hit2[2]/r2);
 539                           delta = (charge1*delta1 + charge2*delta2)/(charge1+charge2);
 540                           hit2[0] = hit2[0] - AliL3Transform::Row2X(rowrange[0]);
 541
 542                           theta = atan2(hit2[0]-hit[0],hit[1]-hit2[1]);
 543                           rho = hit[0]*cos(theta)+hit[1]*sin(theta);
 544                           hist->Fill(theta,rho,delta,charge1+charge2);
 545                         }
 546                     }
 547                 }
 548             }
 549         }
 550     }
 551 }
 552
 553 void AliL3HoughTest::TransformLines2Circle(TH3 *hist,AliL3TrackArray *tracks)
 554 {
 555   //Another HT?
 556   for(Int_t i=0; i<tracks->GetNTracks(); i++)
 557     {
 558       AliL3HoughTrack *tr = (AliL3HoughTrack*)tracks->GetCheckedTrack(i);
 559       if(!tr) continue;
 560       Int_t middlerow = (Int_t)(tr->GetFirstRow()+(tr->GetLastRow()-tr->GetFirstRow())/2);
 561       Float_t hit[3];
 562       tr->GetLineCrossingPoint(middlerow,hit);
 563       hit[0] += AliL3Transform::Row2X(tr->GetFirstRow());
 564       Float_t r = sqrt(hit[0]*hit[0] + hit[1]*hit[1]);
 565       hit[2] = r*tr->GetTgl();
 566       Float_t phi = atan2(hit[1],hit[0]);
 567       Float_t theta = tr->GetPsiLine() - AliL3Transform::Pi()/2;
 568       Float_t psi = 2*phi - theta;
 569       Float_t kappa = 2/r*sin(phi-psi);
 570       Float_t delta = atan(hit[2]/r);
 571       hist->Fill(kappa,psi,delta,tr->GetWeight());
 572
 573     }
 574 }
 575
 576 void AliL3HoughTest::Transform2Center(AliL3Histogram *hist)
 577 {
 578   //Choose parameter space to be center of curvature.
 579
 580   if(!fData)
 581     {
 582       cerr<<"AliL3HoughTest::TransformC : No data"<<endl;
 583       return;
 584     }
 585   Int_t pad1,pad2,time1,time2,charge1,charge2,sector,row;
 586   Float_t r1,r2,phi1,phi2,hit[3],hit2[3];
 587   //Float_t phi_0,kappa;
 588   for(Int_t i=AliL3Transform::GetFirstRow(fCurrentPatch); i<=AliL3Transform::GetLastRow(fCurrentPatch); i++)
 589     {
 590       Int_t rowindex1 = i - AliL3Transform::GetFirstRow(fCurrentPatch);
 591       for(Int_t d1=0; d1<fData[rowindex1].fNPads; d1++)
 592         {
 593           pad1 = d1 + fData[rowindex1].fMinpad;
 594           for(Int_t j=0; j<10; j++)
 595             {
 596               time1 = j + fData[rowindex1].fMintime;
 597               charge1 = fData[rowindex1].fPads[d1][j];
 598               if(charge1==0) continue;
 599               AliL3Transform::Slice2Sector(0,i,sector,row);
 600               AliL3Transform::Raw2Local(hit,sector,row,pad1,time1);
 601               r1 = sqrt(hit[0]*hit[0]+hit[1]*hit[1])/2;
 602               phi1 = atan2(hit[1],hit[0]);
 603
 604               for(Int_t j=i+1; j<=AliL3Transform::GetLastRow(fCurrentPatch); j++)
 605                 {
 606                   Int_t rowindex2 = j - AliL3Transform::GetFirstRow(fCurrentPatch);
 607                   for(Int_t d2=0; d2<fData[rowindex2].fNPads; d2++)
 608                     {
 609                       pad2 = d2 + fData[rowindex2].fMinpad;
 610                       for(Int_t k=0; k<10; k++)
 611                         {
 612                           time2 = k + fData[rowindex2].fMintime;
 613                           charge2 = fData[rowindex2].fPads[d2][k];
 614                           if(charge2==0) continue;
 615                           AliL3Transform::Slice2Sector(0,j,sector,row);
 616                           AliL3Transform::Raw2Local(hit2,sector,row,pad2,time2);
 617                           r2 = sqrt(hit2[0]*hit2[0]+hit2[1]*hit2[1])/2;
 618                           phi2 = atan2(hit2[1],hit2[0]);
 619                           Float_t yc = (r2-(r1/cos(phi1))*cos(phi2))/(sin(phi2)-tan(phi1)*cos(phi2));
 620                           Float_t xc = r1/cos(phi1) - (r2*tan(phi1)-r1*sin(phi1)*cos(phi2))/(sin(phi2)-tan(phi1)*cos(phi2));
 621                           hist->Fill(xc,yc,charge1+charge2);
 622                         }
 623                     }
 624                 }
 625             }
 626         }
 627     }
 628 }
 629
 630
 631 void AliL3HoughTest::FindAbsMaxima(TH3 *hist,Int_t zsearch,Float_t &maxx,Float_t &maxy,Float_t &maxz,Int_t &maxvalue) const
 632 {
 633   //Find peaks in the Hough space
 634   TH1 *h1 = hist->Project3D("z");
 635
 636   TAxis *z = hist->GetZaxis();
 637   Float_t zpeak[50];
 638   Int_t n=0,i,zbin[50];
 639   for(i=z->GetFirst()+1; i<z->GetLast()-1; i++)
 640     {
 641       int bin1 = h1->GetBin(i-1);
 642       int bin2 = h1->GetBin(i);
 643       int bin3 = h1->GetBin(i+1);
 644       if(h1->GetBinContent(bin2) > h1->GetBinContent(bin1) && h1->GetBinContent(bin2) > h1->GetBinContent(bin3))
 645         {
 646           zbin[n]=bin2;
 647           zpeak[n++]=h1->GetBinCenter(bin2);
 648         }
 649     }
 650
 651   Int_t zrange[2] = {z->GetFirst(),z->GetLast()};
 652   z->SetRange(zbin[0]-zsearch,zbin[0]+zsearch);
 653
 654   TH1 *h2 = hist->Project3D("yx");
 655   z->SetRange(zrange[0],zrange[1]);
 656
 657   TAxis *x = h2->GetXaxis();
 658   TAxis *y = h2->GetYaxis();
 659   maxvalue=0;
 660   for(i=0; i<n; i++)
 661     {
 662       for(Int_t xbin=x->GetFirst(); xbin<=x->GetLast(); xbin++)
 663         {
 664           Float_t xvalue = x->GetBinCenter(xbin);
 665           for(Int_t ybin=y->GetFirst(); ybin<=y->GetLast(); ybin++)
 666             {
 667               Float_t yvalue = y->GetBinCenter(ybin);
 668               Int_t value = (Int_t)h2->GetBinContent(xbin,ybin);
 669               if(value > maxvalue)
 670                 {
 671                   maxvalue=value;
 672                   maxx = xvalue;
 673                   maxy = yvalue;
 674                   maxz = zpeak[i];
 675                 }
 676             }
 677         }
 678     }
 679
 680   delete h1;
 681   delete h2;
 682 }
 683