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