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