HLT/hough/AliL3Hough.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 <sys/time.h>
   8
   9 #include "AliL3Logging.h"
  10 #include "AliL3HoughMerger.h"
  11 #include "AliL3HoughIntMerger.h"
  12 #include "AliL3HoughGlobalMerger.h"
  13 #include "AliL3Histogram.h"
  14 #include "AliL3Hough.h"
  15 #include "AliL3HoughTransformer.h"
  16 #include "AliL3HoughClusterTransformer.h"
  17 #include "AliL3HoughTransformerLUT.h"
  18 #include "AliL3HoughTransformerVhdl.h"
  19 #include "AliL3HoughTransformerRow.h"
  20 #include "AliL3HoughMaxFinder.h"
  21 #include "AliL3Benchmark.h"
  22 #ifdef use_aliroot
  23 #include "AliL3FileHandler.h"
  24 #else
  25 #include "AliL3MemHandler.h"
  26 #endif
  27 #include "AliL3DataHandler.h"
  28 #include "AliL3DigitData.h"
  29 #include "AliL3HoughEval.h"
  30 #include "AliL3Transform.h"
  31 #include "AliL3TrackArray.h"
  32 #include "AliL3HoughTrack.h"
  33 #include "AliL3DDLDataFileHandler.h"
  34
  35 #include "TThread.h"
  36
  37 #if __GNUC__ == 3
  38 using namespace std;
  39 #endif
  40
  41 /** /class AliL3Hough
  42 //<pre>
  43 //_____________________________________________________________
  44 // AliL3Hough
  45 //
  46 // Interface class for the Hough transform
  47 //
  48 // Example how to use:
  49 //
  50 // AliL3Hough *hough = new AliL3Hough(path,kTRUE,NumberOfEtaSegments);
  51 // hough->ReadData(slice);
  52 // hough->Transform();
  53 // hough->FindTrackCandidates();
  54 //
  55 // AliL3TrackArray *tracks = hough->GetTracks(patch);
  56 //
  57 //</pre>
  58 */
  59
  60 ClassImp(AliL3Hough)
  61
  62 AliL3Hough::AliL3Hough()
  63 {
  64   //Constructor
  65
  66   fBinary        = kFALSE;
  67   fAddHistograms = kFALSE;
  68   fDoIterative   = kFALSE;
  69   fWriteDigits   = kFALSE;
  70   fUse8bits      = kFALSE;
  71
  72   fMemHandler       = 0;
  73   fHoughTransformer = 0;
  74   fEval             = 0;
  75   fPeakFinder       = 0;
  76   fTracks           = 0;
  77   fGlobalTracks     = 0;
  78   fMerger           = 0;
  79   fInterMerger      = 0;
  80   fGlobalMerger     = 0;
  81   fBenchmark        = 0;
  82
  83   fNEtaSegments     = 0;
  84   fNPatches         = 0;
  85   fLastPatch        =-1;
  86   fVersion          = 0;
  87   fCurrentSlice     = 0;
  88   fEvent            = 0;
  89
  90   fKappaSpread = 6;
  91   fPeakRatio   = 0.5;
  92   fInputFile   = 0;
  93   fInputPtr    = 0;
  94   fRawEvent    = 0;
  95
  96   SetTransformerParams();
  97   SetThreshold();
  98   SetNSaveIterations();
  99   SetPeakThreshold();
 100 #ifdef use_aliroot
 101   //just be sure that index is empty for new event
 102     AliL3FileHandler::CleanStaticIndex();
 103 #ifdef use_newio
 104     fRunLoader = 0;
 105 #endif
 106 #endif
 107   fThread = 0;
 108 }
 109
 110 AliL3Hough::AliL3Hough(Char_t *path,Bool_t binary,Int_t netasegments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr)
 111 {
 112   //Normal constructor
 113   fBinary = binary;
 114   strcpy(fPath,path);
 115   fNEtaSegments  = netasegments;
 116   fAddHistograms = kFALSE;
 117   fDoIterative   = kFALSE;
 118   fWriteDigits   = kFALSE;
 119   fUse8bits      = bit8;
 120   fVersion       = tv;
 121   fKappaSpread=6;
 122   fPeakRatio=0.5;
 123   if(!fBinary) {
 124     if(infile) {
 125       fInputFile = infile;
 126       fInputPtr = 0;
 127     }
 128     else {
 129       fInputFile = 0;
 130       fInputPtr = ptr;
 131     }
 132   }
 133   else {
 134     fInputFile = 0;
 135     fInputPtr = 0;
 136   }
 137 #ifdef use_aliroot
 138   //just be sure that index is empty for new event
 139     AliL3FileHandler::CleanStaticIndex();
 140 #ifdef use_newio
 141     fRunLoader = 0;
 142 #endif
 143 #endif
 144   fThread = 0;
 145 }
 146
 147 AliL3Hough::~AliL3Hough()
 148 {
 149   //dtor
 150
 151   CleanUp();
 152   if(fMerger)
 153     delete fMerger;
 154   //cout << "Cleaned class merger " << endl;
 155   if(fInterMerger)
 156     delete fInterMerger;
 157   //cout << "Cleaned class inter " << endl;
 158   if(fPeakFinder)
 159     delete fPeakFinder;
 160   //cout << "Cleaned class peak " << endl;
 161   if(fGlobalMerger)
 162     delete fGlobalMerger;
 163   //cout << "Cleaned class global " << endl;
 164   if(fBenchmark)
 165     delete fBenchmark;
 166   //cout << "Cleaned class bench " << endl;
 167   if(fGlobalTracks)
 168     delete fGlobalTracks;
 169   //cout << "Cleaned class globaltracks " << endl;
 170   if(fThread) {
 171     //    fThread->Delete();
 172     delete fThread;
 173     fThread = 0;
 174   }
 175 }
 176
 177 void AliL3Hough::CleanUp()
 178 {
 179   //Cleanup memory
 180
 181   for(Int_t i=0; i<fNPatches; i++)
 182     {
 183       if(fTracks[i]) delete fTracks[i];
 184       //cout << "Cleaned tracks " << i << endl;
 185       if(fEval[i]) delete fEval[i];
 186       //cout << "Cleaned eval " << i << endl;
 187       if(fHoughTransformer[i]) delete fHoughTransformer[i];
 188       //cout << "Cleaned traf " << i << endl;
 189       if(fMemHandler[i]) delete fMemHandler[i];
 190       //cout << "Cleaned mem " << i << endl;
 191     }
 192
 193   if(fTracks) delete [] fTracks;
 194   //cout << "Cleaned class tracks " << endl;
 195   if(fEval) delete [] fEval;
 196   //cout << "Cleaned class eval " << endl;
 197   if(fHoughTransformer) delete [] fHoughTransformer;
 198   //cout << "Cleaned cleass trafo " << endl;
 199   if(fMemHandler) delete [] fMemHandler;
 200   //cout << "Cleaned class mem " << endl;
 201 }
 202
 203 void AliL3Hough::Init(Int_t netasegments,Int_t tv,AliRawEvent *rawevent,Float_t zvertex)
 204 {
 205   //Normal constructor
 206   fNEtaSegments  = netasegments;
 207   fVersion       = tv;
 208   fRawEvent      = rawevent;
 209   fZVertex       = zvertex;
 210
 211   Init();
 212 }
 213
 214 void AliL3Hough::Init(Char_t *path,Bool_t binary,Int_t netasegments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr,Float_t zvertex)
 215 {
 216   //Normal init of the AliL3Hough
 217   fBinary = binary;
 218   strcpy(fPath,path);
 219   fNEtaSegments = netasegments;
 220   fWriteDigits  = kFALSE;
 221   fUse8bits     = bit8;
 222   fVersion      = tv;
 223   if(!fBinary) {
 224     if(infile) {
 225       fInputFile = infile;
 226       fInputPtr = 0;
 227     }
 228     else {
 229       fInputFile = 0;
 230       fInputPtr = ptr;
 231     }
 232   }
 233   else {
 234     fInputFile = 0;
 235     fInputPtr = 0;
 236   }
 237   fZVertex = zvertex;
 238
 239   Init(); //do the rest
 240 }
 241
 242 void AliL3Hough::Init(Bool_t doit, Bool_t addhists)
 243 {
 244   // Init
 245   fDoIterative   = doit;
 246   fAddHistograms = addhists;
 247
 248   fNPatches = AliL3Transform::GetNPatches();
 249   fHoughTransformer = new AliL3HoughBaseTransformer*[fNPatches];
 250   fMemHandler = new AliL3MemHandler*[fNPatches];
 251
 252   fTracks = new AliL3TrackArray*[fNPatches];
 253   fEval = new AliL3HoughEval*[fNPatches];
 254
 255   fGlobalTracks = new AliL3TrackArray("AliL3HoughTrack");
 256
 257   AliL3HoughBaseTransformer *lasttransformer = 0;
 258
 259   for(Int_t i=0; i<fNPatches; i++)
 260     {
 261       switch (fVersion){ //choose Transformer
 262       case 1:
 263         fHoughTransformer[i] = new AliL3HoughTransformerLUT(0,i,fNEtaSegments);
 264         break;
 265       case 2:
 266         fHoughTransformer[i] = new AliL3HoughClusterTransformer(0,i,fNEtaSegments);
 267         break;
 268       case 3:
 269         fHoughTransformer[i] = new AliL3HoughTransformerVhdl(0,i,fNEtaSegments,fNSaveIterations);
 270         break;
 271       case 4:
 272         fHoughTransformer[i] = new AliL3HoughTransformerRow(0,i,fNEtaSegments,kFALSE,fZVertex);
 273         break;
 274       default:
 275         fHoughTransformer[i] = new AliL3HoughTransformer(0,i,fNEtaSegments,kFALSE,kFALSE);
 276       }
 277
 278       fHoughTransformer[i]->SetLastTransformer(lasttransformer);
 279       lasttransformer = fHoughTransformer[i];
 280       //      fHoughTransformer[i]->CreateHistograms(fNBinX[i],fLowPt[i],fNBinY[i],-fPhi[i],fPhi[i]);
 281       fHoughTransformer[i]->CreateHistograms(fNBinX[i],-fLowPt[i],fLowPt[i],fNBinY[i],-fPhi[i],fPhi[i]);
 282       //fHoughTransformer[i]->CreateHistograms(fLowPt[i],fUpperPt[i],fPtRes[i],fNBinY[i],fPhi[i]);
 283
 284       fHoughTransformer[i]->SetLowerThreshold(fThreshold[i]);
 285       fHoughTransformer[i]->SetUpperThreshold(100);
 286
 287       LOG(AliL3Log::kInformational,"AliL3Hough::Init","Version")
 288         <<"Initializing Hough transformer version "<<fVersion<<ENDLOG;
 289
 290       fEval[i] = new AliL3HoughEval();
 291       fTracks[i] = new AliL3TrackArray("AliL3HoughTrack");
 292       if(fUse8bits)
 293         fMemHandler[i] = new AliL3DataHandler();
 294       else
 295 #ifdef use_aliroot
 296         {
 297           if(!fRawEvent) {
 298             if(!fInputFile) {
 299               if(!fInputPtr) {
 300                 /* In case of reading digits file */
 301                 fMemHandler[i] = new AliL3FileHandler(kTRUE); //use static index
 302                 if(!fBinary) {
 303 #if use_newio
 304                   if(!fRunLoader) {
 305 #endif
 306                     Char_t filename[1024];
 307                     sprintf(filename,"%s/digitfile.root",fPath);
 308                     fMemHandler[i]->SetAliInput(filename);
 309 #if use_newio
 310                   }
 311                   else {
 312                     fMemHandler[i]->SetAliInput(fRunLoader);
 313                   }
 314 #endif
 315                 }
 316               }
 317               else {
 318                 /* In case of reading from DATE */
 319                 fMemHandler[i] = new AliL3DDLDataFileHandler();
 320                 fMemHandler[i]->SetReaderInput(fInputPtr,-1);
 321               }
 322             }
 323             else {
 324               /* In case of reading rawdata from ROOT file */
 325               fMemHandler[i] = new AliL3DDLDataFileHandler();
 326               fMemHandler[i]->SetReaderInput(fInputFile);
 327             }
 328           }
 329           else {
 330             /* In case of reading rawdata using AliRawEvent */
 331             fMemHandler[i] = new AliL3DDLDataFileHandler();
 332             fMemHandler[i]->SetReaderInput(fRawEvent);
 333           }
 334         }
 335 #else
 336       fMemHandler[i] = new AliL3MemHandler();
 337 #endif
 338     }
 339
 340   fPeakFinder = new AliL3HoughMaxFinder("KappaPhi",50000);
 341   fMerger = new AliL3HoughMerger(fNPatches);
 342   fInterMerger = new AliL3HoughIntMerger();
 343   fGlobalMerger = 0;
 344   fBenchmark = new AliL3Benchmark();
 345 }
 346
 347 void AliL3Hough::SetTransformerParams(Float_t ptres,Float_t ptmin,Float_t ptmax,Int_t ny,Int_t patch)
 348 {
 349   // Setup the parameters for the Hough Transformer
 350   Int_t mrow;
 351   Float_t psi=0;
 352   if(patch==-1)
 353     mrow = 80;
 354   else
 355     mrow = AliL3Transform::GetLastRow(patch);
 356   if(ptmin)
 357     {
 358       Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(mrow),2) + pow(AliL3Transform::GetMaxY(mrow),2));
 359       Double_t kappa = -1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/ptmin;
 360       psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
 361       cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
 362     }
 363
 364   if(patch==-1)
 365     {
 366       Int_t i=0;
 367       while(i < 6)
 368         {
 369           fPtRes[i] = ptres;
 370           fLowPt[i] = ptmin;
 371           fUpperPt[i] = ptmax;
 372           fNBinY[i] = ny;
 373           fPhi[i] = psi;
 374           fNBinX[i]=0;
 375           i++;
 376         }
 377       return;
 378     }
 379
 380   fPtRes[patch] = ptres;
 381   fLowPt[patch] = ptmin;
 382   fUpperPt[patch] = ptmax;
 383   fNBinY[patch] = ny;
 384   fPhi[patch] = psi;
 385 }
 386 /*
 387 void AliL3Hough::SetTransformerParams(Int_t nx,Int_t ny,Float_t ptmin,Int_t patch)
 388 {
 389   // Setup the parameters for the Hough Transformer
 390
 391   Int_t mrow=80;
 392   Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(mrow),2) + pow(AliL3Transform::GetMaxY(mrow),2));
 393   Double_t kappa = -1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/ptmin;
 394   Double_t psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
 395   cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
 396
 397   Int_t i=0;
 398   while(i < 6)
 399     {
 400       fLowPt[i] = ptmin;
 401       fNBinY[i] = ny;
 402       fNBinX[i] = nx;
 403       fPhi[i] = psi;
 404       i++;
 405     }
 406 }
 407 */
 408 void AliL3Hough::SetTransformerParams(Int_t nx,Int_t ny,Float_t ptmin,Int_t /*patch*/)
 409 {
 410   // Setup the parameters for the Hough Transformer
 411
 412
 413   Int_t mrow=79;
 414   Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(mrow),2) + pow(AliL3Transform::GetMaxY(mrow),2));
 415   Double_t alpha1 = AliL3Transform::GetMaxY(mrow)/pow(lineradius,2);
 416   Double_t kappa = 1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/ptmin;
 417   Double_t psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
 418   //  cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
 419   AliL3HoughTrack track;
 420   track.SetTrackParameters(kappa,psi,1);
 421   Float_t hit[3];
 422   Int_t mrow2 = 158;
 423   track.GetCrossingPoint(mrow2,hit);
 424   Double_t lineradius2 = sqrt(pow(AliL3Transform::Row2X(mrow2),2) + pow(AliL3Transform::GetMaxY(mrow2),2));
 425   Double_t alpha2 = hit[1]/pow(lineradius2,2);
 426   //  cout<<"Calculated alphas range "<<alpha1<<" "<<alpha2<<" in patch "<<patch<<endl;
 427
 428   Int_t i=0;
 429   while(i < 6)
 430     {
 431       fLowPt[i] = 1.15*alpha1;
 432       fNBinY[i] = ny;
 433       fNBinX[i] = nx;
 434       fPhi[i] = 1.15*alpha2;
 435       i++;
 436     }
 437 }
 438
 439 void AliL3Hough::SetTransformerParams(Int_t nx,Int_t ny,Float_t lpt,Float_t phi)
 440 {
 441   Int_t i=0;
 442   while(i < 6)
 443     {
 444       fLowPt[i] = lpt;
 445       fNBinY[i] = ny;
 446       fNBinX[i] = nx;
 447       fPhi[i] = phi;
 448       i++;
 449     }
 450 }
 451
 452 void AliL3Hough::SetThreshold(Int_t t3,Int_t patch)
 453 {
 454   // Set digits threshold
 455   if(patch==-1)
 456     {
 457       Int_t i=0;
 458       while(i < 6)
 459         fThreshold[i++]=t3;
 460       return;
 461     }
 462   fThreshold[patch]=t3;
 463 }
 464
 465 void AliL3Hough::SetPeakThreshold(Int_t threshold,Int_t patch)
 466 {
 467   // Set Peak Finder threshold
 468   if(patch==-1)
 469     {
 470       Int_t i=0;
 471       while(i < 6)
 472         fPeakThreshold[i++]=threshold;
 473       return;
 474     }
 475   fPeakThreshold[patch]=threshold;
 476 }
 477
 478 void AliL3Hough::DoBench(Char_t *name)
 479 {
 480   fBenchmark->Analyze(name);
 481 }
 482
 483 void AliL3Hough::Process(Int_t minslice,Int_t maxslice)
 484 {
 485   //Process all slices [minslice,maxslice].
 486   fGlobalMerger = new AliL3HoughGlobalMerger(minslice,maxslice);
 487
 488   for(Int_t i=minslice; i<=maxslice; i++)
 489     {
 490       ReadData(i);
 491       Transform();
 492       if(fAddHistograms) {
 493         if(fVersion != 4)
 494           AddAllHistograms();
 495         else
 496           AddAllHistogramsRows();
 497       }
 498       FindTrackCandidates();
 499       //Evaluate();
 500       //fGlobalMerger->FillTracks(fTracks[0],i);
 501     }
 502 }
 503
 504 void AliL3Hough::ReadData(Int_t slice,Int_t eventnr)
 505 {
 506   //Read data from files, binary or root.
 507
 508 #ifdef use_aliroot
 509   if(fEvent!=eventnr) //just be sure that index is empty for new event
 510     AliL3FileHandler::CleanStaticIndex();
 511 #endif
 512   fCurrentSlice = slice;
 513
 514   for(Int_t i=0; i<fNPatches; i++)
 515     {
 516       fMemHandler[i]->Free();
 517       UInt_t ndigits=0;
 518       AliL3DigitRowData *digits =0;
 519       Char_t name[256];
 520       fMemHandler[i]->Init(slice,i);
 521       if(fBinary)//take input data from binary files
 522         {
 523           if(fUse8bits)
 524             sprintf(name,"%s/binaries/digits_c8_%d_%d_%d.raw",fPath,eventnr,slice,i);
 525           else
 526             sprintf(name,"%s/binaries/digits_%d_%d_%d.raw",fPath,eventnr,slice,i);
 527
 528           fMemHandler[i]->SetBinaryInput(name);
 529           digits = (AliL3DigitRowData *)fMemHandler[i]->CompBinary2Memory(ndigits);
 530           fMemHandler[i]->CloseBinaryInput();
 531         }
 532       else //read data from root file
 533         {
 534 #ifdef use_aliroot
 535           if(fEvent!=eventnr)
 536             fMemHandler[i]->FreeDigitsTree();//or else the new event is not loaded
 537           digits=(AliL3DigitRowData *)fMemHandler[i]->AliAltroDigits2Memory(ndigits,eventnr);
 538 #else
 539           cerr<<"You cannot read from rootfile now"<<endl;
 540 #endif
 541         }
 542
 543       //Set the pointer to the TPCRawStream in case of fast raw data reading
 544       fHoughTransformer[i]->SetTPCRawStream(fMemHandler[i]->GetTPCRawStream());
 545
 546       //set input data and init transformer
 547       fHoughTransformer[i]->SetInputData(ndigits,digits);
 548       fHoughTransformer[i]->Init(slice,i,fNEtaSegments);
 549     }
 550
 551   fEvent=eventnr;
 552 }
 553
 554 void AliL3Hough::Transform(Int_t *rowrange)
 555 {
 556   //Transform all data given to the transformer within the given slice
 557   //(after ReadData(slice))
 558
 559   Double_t initTime,cpuTime;
 560   initTime = GetCpuTime();
 561   Int_t patchorder[6] = {5,2,0,1,3,4}; //The order in which patches are processed
 562   //  Int_t patchorder[6] = {0,1,2,3,4,5}; //The order in which patches are processed
 563   //  Int_t patchorder[6] = {5,4,3,2,1,0}; //The order in which patches are processed
 564   //  Int_t patchorder[6] = {5,2,4,3,1,0}; //The order in which patches are processed
 565   fLastPatch=-1;
 566   for(Int_t i=0; i<fNPatches; i++)
 567     {
 568       // In case of Row transformer reset the arrays only once
 569       if((fVersion != 4) || (i == 0)) {
 570         fBenchmark->Start("Hough Reset");
 571         fHoughTransformer[0]->Reset();//Reset the histograms
 572         fBenchmark->Stop("Hough Reset");
 573       }
 574       fBenchmark->Start("Hough Transform");
 575       PrepareForNextPatch(patchorder[i]);
 576       if(!rowrange) {
 577         char buf[256];
 578         sprintf(buf,"Patch %d",patchorder[i]);
 579         fBenchmark->Start(buf);
 580         fHoughTransformer[patchorder[i]]->SetLastPatch(fLastPatch);
 581         fHoughTransformer[patchorder[i]]->TransformCircle();
 582         fBenchmark->Stop(buf);
 583       }
 584       else
 585         fHoughTransformer[i]->TransformCircleC(rowrange,1);
 586       fBenchmark->Stop("Hough Transform");
 587       fLastPatch=patchorder[i];
 588     }
 589   cpuTime = GetCpuTime() - initTime;
 590   LOG(AliL3Log::kInformational,"AliL3Hough::Transform()","Timing")
 591     <<"Transform done in average per patch of "<<cpuTime*1000/fNPatches<<" ms"<<ENDLOG;
 592 }
 593
 594 void AliL3Hough::MergePatches()
 595 {
 596   // Merge patches if they are not summed
 597   if(fAddHistograms) //Nothing to merge here
 598     return;
 599   fMerger->MergePatches(kTRUE);
 600 }
 601
 602 void AliL3Hough::MergeInternally()
 603 {
 604   // Merge patches internally
 605   if(fAddHistograms)
 606     fInterMerger->FillTracks(fTracks[0]);
 607   else
 608     fInterMerger->FillTracks(fMerger->GetOutTracks());
 609
 610   fInterMerger->MMerge();
 611 }
 612
 613 void AliL3Hough::ProcessSliceIter()
 614 {
 615   //Process current slice (after ReadData(slice)) iteratively.
 616
 617   if(!fAddHistograms)
 618     {
 619       for(Int_t i=0; i<fNPatches; i++)
 620         {
 621           ProcessPatchIter(i);
 622           fMerger->FillTracks(fTracks[i],i); //Copy tracks to merger
 623         }
 624     }
 625   else
 626     {
 627       for(Int_t i=0; i<10; i++)
 628         {
 629           Transform();
 630           AddAllHistograms();
 631           InitEvaluate();
 632           AliL3HoughBaseTransformer *tr = fHoughTransformer[0];
 633           for(Int_t j=0; j<fNEtaSegments; j++)
 634             {
 635               AliL3Histogram *hist = tr->GetHistogram(j);
 636               if(hist->GetNEntries()==0) continue;
 637               fPeakFinder->Reset();
 638               fPeakFinder->SetHistogram(hist);
 639               fPeakFinder->FindAbsMaxima();
 640               AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[0]->NextTrack();
 641               track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
 642               track->SetEtaIndex(j);
 643               track->SetEta(tr->GetEta(j,fCurrentSlice));
 644               for(Int_t k=0; k<fNPatches; k++)
 645                 {
 646                   fEval[i]->SetNumOfPadsToLook(2);
 647                   fEval[i]->SetNumOfRowsToMiss(2);
 648                   fEval[i]->RemoveFoundTracks();
 649                   /*
 650                   Int_t nrows=0;
 651                   if(!fEval[i]->LookInsideRoad(track,nrows))
 652                     {
 653                       fTracks[0]->Remove(fTracks[0]->GetNTracks()-1);
 654                       fTracks[0]->Compress();
 655                     }
 656                   */
 657                 }
 658             }
 659
 660         }
 661
 662     }
 663 }
 664
 665 void AliL3Hough::ProcessPatchIter(Int_t patch)
 666 {
 667   //Process patch in a iterative way.
 668   //transform + peakfinding + evaluation + transform +...
 669
 670   Int_t numoftries = 5;
 671   AliL3HoughBaseTransformer *tr = fHoughTransformer[patch];
 672   AliL3TrackArray *tracks = fTracks[patch];
 673   tracks->Reset();
 674   AliL3HoughEval *ev = fEval[patch];
 675   ev->InitTransformer(tr);
 676   //ev->RemoveFoundTracks();
 677   ev->SetNumOfRowsToMiss(3);
 678   ev->SetNumOfPadsToLook(2);
 679   AliL3Histogram *hist;
 680   for(Int_t t=0; t<numoftries; t++)
 681     {
 682       tr->Reset();
 683       tr->TransformCircle();
 684       for(Int_t i=0; i<fNEtaSegments; i++)
 685         {
 686           hist = tr->GetHistogram(i);
 687           if(hist->GetNEntries()==0) continue;
 688           fPeakFinder->Reset();
 689           fPeakFinder->SetHistogram(hist);
 690           fPeakFinder->FindAbsMaxima();
 691           //fPeakFinder->FindPeak1();
 692           AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->NextTrack();
 693           track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
 694           track->SetEtaIndex(i);
 695           track->SetEta(tr->GetEta(i,fCurrentSlice));
 696           /*
 697           Int_t nrows=0;
 698           if(!ev->LookInsideRoad(track,nrows))
 699             {
 700               tracks->Remove(tracks->GetNTracks()-1);
 701               tracks->Compress();
 702             }
 703           */
 704         }
 705     }
 706   fTracks[0]->QSort();
 707   LOG(AliL3Log::kInformational,"AliL3Hough::ProcessPatch","NTracks")
 708     <<AliL3Log::kDec<<"Found "<<tracks->GetNTracks()<<" tracks in patch "<<patch<<ENDLOG;
 709 }
 710
 711 void AliL3Hough::AddAllHistograms()
 712 {
 713   //Add the histograms within one etaslice.
 714   //Resulting histogram are in patch=0.
 715
 716   Double_t initTime,cpuTime;
 717   initTime = GetCpuTime();
 718   fBenchmark->Start("Add Histograms");
 719   for(Int_t i=0; i<fNEtaSegments; i++)
 720     {
 721       AliL3Histogram *hist0 = fHoughTransformer[0]->GetHistogram(i);
 722       for(Int_t j=1; j<fNPatches; j++)
 723         {
 724           AliL3Histogram *hist = fHoughTransformer[j]->GetHistogram(i);
 725           hist0->Add(hist);
 726         }
 727     }
 728   fBenchmark->Stop("Add Histograms");
 729   fAddHistograms = kTRUE;
 730   cpuTime = GetCpuTime() - initTime;
 731   LOG(AliL3Log::kInformational,"AliL3Hough::AddAllHistograms()","Timing")
 732     <<"Adding histograms in "<<cpuTime*1000<<" ms"<<ENDLOG;
 733 }
 734
 735 void AliL3Hough::AddAllHistogramsRows()
 736 {
 737   //Add the histograms within one etaslice.
 738   //Resulting histogram are in patch=0.
 739
 740   Double_t initTime,cpuTime;
 741   initTime = GetCpuTime();
 742   fBenchmark->Start("Add HistogramsRows");
 743
 744   UChar_t lastpatchlastrow = AliL3Transform::GetLastRowOnDDL(fLastPatch)+1;
 745
 746   UChar_t *tracklastrow = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetTrackLastRow();
 747
 748   for(Int_t i=0; i<fNEtaSegments; i++)
 749     {
 750       UChar_t *gapcount = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetGapCount(i);
 751       UChar_t *currentrowcount = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetCurrentRowCount(i);
 752
 753       AliL3Histogram *hist = fHoughTransformer[0]->GetHistogram(i);
 754       Int_t xmin = hist->GetFirstXbin();
 755       Int_t xmax = hist->GetLastXbin();
 756       Int_t ymin = hist->GetFirstYbin();
 757       Int_t ymax = hist->GetLastYbin();
 758       Int_t nxbins = hist->GetNbinsX()+2;
 759
 760       for(Int_t ybin=ymin; ybin<=ymax; ybin++)
 761         {
 762           for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 763             {
 764               Int_t bin = xbin + ybin*nxbins; //Int_t bin = hist->GetBin(xbin,ybin);
 765               if(gapcount[bin] < MAX_N_GAPS) {
 766                 if(tracklastrow[bin] > lastpatchlastrow) {
 767                   if(lastpatchlastrow > currentrowcount[bin])
 768                     gapcount[bin] += (lastpatchlastrow-currentrowcount[bin]-1);
 769                 }
 770                 else {
 771                   if(tracklastrow[bin] > currentrowcount[bin])
 772                     gapcount[bin] += (tracklastrow[bin]-currentrowcount[bin]-1);
 773                 }
 774                 if(gapcount[bin] < MAX_N_GAPS)
 775                   hist->AddBinContent(bin,(159-gapcount[bin]));
 776               }
 777             }
 778         }
 779     }
 780
 781   fBenchmark->Stop("Add HistogramsRows");
 782   fAddHistograms = kTRUE;
 783   cpuTime = GetCpuTime() - initTime;
 784   LOG(AliL3Log::kInformational,"AliL3Hough::AddAllHistogramsRows()","Timing")
 785     <<"Adding histograms in "<<cpuTime*1000<<" ms"<<ENDLOG;
 786 }
 787
 788 void AliL3Hough::PrepareForNextPatch(Int_t nextpatch)
 789 {
 790   char buf[256];
 791   sprintf(buf,"Prepare For Patch %d",nextpatch);
 792   fBenchmark->Start(buf);
 793
 794   UChar_t lastpatchlastrow;
 795   if(fLastPatch == -1)
 796     lastpatchlastrow = 0;
 797   else
 798     lastpatchlastrow = AliL3Transform::GetLastRowOnDDL(fLastPatch)+1;
 799   UChar_t nextpatchfirstrow;
 800   if(nextpatch==0)
 801     nextpatchfirstrow = 0;
 802   else
 803     nextpatchfirstrow = AliL3Transform::GetFirstRowOnDDL(nextpatch)-1;
 804
 805   UChar_t *trackfirstrow = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetTrackFirstRow();
 806   UChar_t *tracklastrow = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetTrackLastRow();
 807
 808   for(Int_t i=0; i<fNEtaSegments; i++)
 809     {
 810       UChar_t *gapcount = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetGapCount(i);
 811       UChar_t *currentrowcount = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetCurrentRowCount(i);
 812       UChar_t *prevbin = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetPrevBin(i);
 813       UChar_t *nextbin = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetNextBin(i);
 814       UChar_t *nextrow = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetNextRow(i);
 815
 816       AliL3Histogram *hist = fHoughTransformer[0]->GetHistogram(i);
 817       Int_t xmin = hist->GetFirstXbin();
 818       Int_t xmax = hist->GetLastXbin();
 819       Int_t ymin = hist->GetFirstYbin();
 820       Int_t ymax = hist->GetLastYbin();
 821       Int_t nxbins = hist->GetNbinsX()+2;
 822
 823       UChar_t lastyvalue = 0;
 824       Int_t endybin = ymin - 1;
 825       for(Int_t ybin=ymin; ybin<=ymax; ybin++)
 826         {
 827           UChar_t lastxvalue = 0;
 828           UChar_t maxvalue = 0;
 829           Int_t endxbin = xmin - 1;
 830           for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 831             {
 832               Int_t bin = xbin + ybin*nxbins;
 833               UChar_t value = 0;
 834               if(gapcount[bin] < MAX_N_GAPS) {
 835                 value = 1;
 836                 maxvalue = 1;
 837                 if(tracklastrow[bin] > lastpatchlastrow) {
 838                   if(lastpatchlastrow > currentrowcount[bin])
 839                     gapcount[bin] += (lastpatchlastrow-currentrowcount[bin]-1);
 840                 }
 841                 else {
 842                   if(tracklastrow[bin] > currentrowcount[bin])
 843                     gapcount[bin] += (tracklastrow[bin]-currentrowcount[bin]-1);
 844                 }
 845                 if(trackfirstrow[bin] < nextpatchfirstrow)
 846                   currentrowcount[bin] = nextpatchfirstrow;
 847                 else
 848                   currentrowcount[bin] = trackfirstrow[bin];
 849               }
 850               if(fLastPatch != -1) {
 851                 if(value > 0)
 852                   {
 853                     nextbin[xbin + ybin*nxbins] = (UChar_t)xbin;
 854                     prevbin[xbin + ybin*nxbins] = (UChar_t)xbin;
 855                     if(value > lastxvalue)
 856                       {
 857                         UChar_t *tempnextbin = nextbin + endxbin + 1 + ybin*nxbins;
 858                         memset(tempnextbin,(UChar_t)xbin,xbin-endxbin-1);
 859                       }
 860                     endxbin = xbin;
 861                   }
 862                 else
 863                   {
 864                     prevbin[xbin + ybin*nxbins] = (UChar_t)endxbin;
 865                   }
 866                 lastxvalue = value;
 867               }
 868             }
 869           if(fLastPatch != -1) {
 870             UChar_t *tempnextbin = nextbin + endxbin + 1 + ybin*nxbins;
 871             memset(tempnextbin,(UChar_t)(xmax+1),xmax-endxbin);
 872             if(maxvalue > 0)
 873               {
 874                 nextrow[ybin] = (UChar_t)ybin;
 875                 if(maxvalue > lastyvalue)
 876                   {
 877                     UChar_t *tempnextrow = nextrow + endybin + 1;
 878                     memset(tempnextrow,(UChar_t)ybin,ybin-endybin-1);
 879                   }
 880                 endybin = ybin;
 881               }
 882             lastyvalue = maxvalue;
 883           }
 884         }
 885       if(fLastPatch != -1) {
 886         UChar_t *tempnextrow = nextrow + endybin + 1;
 887         memset(tempnextrow,(UChar_t)(ymax+1),ymax-endybin);
 888       }
 889     }
 890
 891   fBenchmark->Stop(buf);
 892 }
 893
 894 void AliL3Hough::AddTracks()
 895 {
 896   // Add current slice slice tracks to the global list of found tracks
 897   if(!fTracks[0])
 898     {
 899       cerr<<"AliL3Hough::AddTracks : No tracks"<<endl;
 900       return;
 901     }
 902   AliL3TrackArray *tracks = fTracks[0];
 903   for(Int_t i=0; i<tracks->GetNTracks(); i++)
 904     {
 905       AliL3Track *track = tracks->GetCheckedTrack(i);
 906       if(!track) continue;
 907       if(track->GetNHits()!=1) cerr<<"NHITS "<<track->GetNHits()<<endl;
 908       UInt_t *ids = track->GetHitNumbers();
 909       ids[0] = (fCurrentSlice&0x7f)<<25;
 910     }
 911
 912   fGlobalTracks->AddTracks(fTracks[0],0,fCurrentSlice);
 913 }
 914
 915 void AliL3Hough::FindTrackCandidatesRow()
 916 {
 917   // Find AliL3HoughTransformerRow track candidates
 918   if(fVersion != 4) {
 919     LOG(AliL3Log::kError,"AliL3Hough::FindTrackCandidatesRow()","")
 920       <<"Incompatible Peak Finder version!"<<ENDLOG;
 921     return;
 922   }
 923
 924   //Look for peaks in histograms, and find the track candidates
 925   Int_t npatches;
 926   if(fAddHistograms)
 927     npatches = 1; //Histograms have been added.
 928   else
 929     npatches = fNPatches;
 930
 931   Double_t initTime,cpuTime;
 932   initTime = GetCpuTime();
 933   fBenchmark->Start("Find Maxima");
 934   for(Int_t i=0; i<npatches; i++)
 935     {
 936       AliL3HoughBaseTransformer *tr = fHoughTransformer[i];
 937       fTracks[i]->Reset();
 938       fPeakFinder->Reset();
 939
 940       for(Int_t j=0; j<fNEtaSegments; j++)
 941         {
 942           AliL3Histogram *hist = tr->GetHistogram(j);
 943           if(hist->GetNEntries()==0) continue;
 944           fPeakFinder->SetHistogram(hist);
 945           fPeakFinder->SetEtaSlice(j);
 946           fPeakFinder->SetTrackLUTs(((AliL3HoughTransformerRow *)tr)->fTrackNRows,((AliL3HoughTransformerRow *)tr)->fTrackFirstRow,((AliL3HoughTransformerRow *)tr)->fTrackLastRow);
 947 #ifdef do_mc
 948           LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","")
 949             <<"Starting "<<j<<" etaslice"<<ENDLOG;
 950 #endif
 951           fPeakFinder->SetThreshold(fPeakThreshold[i]);
 952           fPeakFinder->FindAdaptedRowPeaks(1,0,0);//Maxima finder for HoughTransformerRow
 953
 954           //fPeakFinder->FindMaxima(fPeakThreshold[i]); //Simple maxima finder
 955         }
 956
 957       for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
 958         {
 959           if(fPeakFinder->GetWeight(k) < 0) continue;
 960           AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[i]->NextTrack();
 961           Float_t psi = atan((fPeakFinder->GetXPeak(k)-fPeakFinder->GetYPeak(k))/(AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2()));
 962           Float_t kappa = 2.0*(fPeakFinder->GetXPeak(k)*cos(psi)-AliL3HoughTransformerRow::GetBeta1()*sin(psi));
 963           //          track->SetTrackParameters(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetWeight(k));
 964           track->SetTrackParameters(kappa,psi,fPeakFinder->GetWeight(k));
 965           track->SetBinXY(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetXPeakSize(k),fPeakFinder->GetYPeakSize(k));
 966           Int_t etaindex = (fPeakFinder->GetStartEta(k)+fPeakFinder->GetEndEta(k))/2;
 967           track->SetEtaIndex(etaindex);
 968           Float_t starteta = tr->GetEta(fPeakFinder->GetStartEta(k),fCurrentSlice);
 969           Float_t endeta = tr->GetEta(fPeakFinder->GetEndEta(k),fCurrentSlice);
 970           track->SetEta((starteta+endeta)/2.0);
 971           track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
 972           track->SetSector(fCurrentSlice);
 973           track->SetSlice(fCurrentSlice);
 974 #ifdef do_mc
 975           Int_t label = tr->GetTrackID(etaindex,fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k));
 976           track->SetMCid(label);
 977           //      cout<<"Track found with label "<<label<<" at "<<fPeakFinder->GetXPeak(k)<<" "<<fPeakFinder->GetYPeak(k)<<" with weight "<<fPeakFinder->GetWeight(k)<<endl;
 978 #endif
 979         }
 980       LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","")
 981         <<"Found "<<fTracks[i]->GetNTracks()<<" tracks in slice "<<fCurrentSlice<<ENDLOG;
 982       fTracks[i]->QSort();
 983     }
 984   fBenchmark->Stop("Find Maxima");
 985   cpuTime = GetCpuTime() - initTime;
 986   LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","Timing")
 987     <<"Maxima finding done in "<<cpuTime*1000<<" ms"<<ENDLOG;
 988 }
 989
 990 void AliL3Hough::FindTrackCandidates()
 991 {
 992   // Find AliL3HoughTransformer track candidates
 993   if(fVersion == 4) {
 994     LOG(AliL3Log::kError,"AliL3Hough::FindTrackCandidatesRow()","")
 995       <<"Incompatible Peak Finder version!"<<ENDLOG;
 996     return;
 997   }
 998
 999   Int_t npatches;
1000   if(fAddHistograms)
1001     npatches = 1; //Histograms have been added.
1002   else
1003     npatches = fNPatches;
1004
1005   Double_t initTime,cpuTime;
1006   initTime = GetCpuTime();
1007   fBenchmark->Start("Find Maxima");
1008   for(Int_t i=0; i<npatches; i++)
1009     {
1010       AliL3HoughBaseTransformer *tr = fHoughTransformer[i];
1011       fTracks[i]->Reset();
1012
1013       for(Int_t j=0; j<fNEtaSegments; j++)
1014         {
1015           AliL3Histogram *hist = tr->GetHistogram(j);
1016           if(hist->GetNEntries()==0) continue;
1017           fPeakFinder->Reset();
1018           fPeakFinder->SetHistogram(hist);
1019 #ifdef do_mc
1020           cout<<"Starting "<<j<<" etaslice"<<endl;
1021 #endif
1022           fPeakFinder->SetThreshold(fPeakThreshold[i]);
1023           fPeakFinder->FindAdaptedPeaks(fKappaSpread,fPeakRatio);
1024
1025           for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
1026             {
1027               AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[i]->NextTrack();
1028               track->SetTrackParameters(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetWeight(k));
1029               track->SetEtaIndex(j);
1030               track->SetEta(tr->GetEta(j,fCurrentSlice));
1031               track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
1032             }
1033         }
1034       cout<<"Found "<<fTracks[i]->GetNTracks()<<" tracks in patch "<<i<<endl;
1035       fTracks[i]->QSort();
1036     }
1037   fBenchmark->Stop("Find Maxima");
1038   cpuTime = GetCpuTime() - initTime;
1039   LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","Timing")
1040     <<"Maxima finding done in "<<cpuTime*1000<<" ms"<<ENDLOG;
1041 }
1042
1043 void AliL3Hough::InitEvaluate()
1044 {
1045   //Pass the transformer objects to the AliL3HoughEval objects:
1046   //This will provide the evaluation objects with all the necessary
1047   //data and parameters it needs.
1048
1049   for(Int_t i=0; i<fNPatches; i++)
1050     fEval[i]->InitTransformer(fHoughTransformer[i]);
1051 }
1052
1053 Int_t AliL3Hough::Evaluate(Int_t roadwidth,Int_t nrowstomiss)
1054 {
1055   //Evaluate the tracks, by looking along the road in the raw data.
1056   //If track does not cross all padrows - rows2miss, it is removed from the arrray.
1057   //If histograms were not added, the check is done locally in patch,
1058   //meaning that nrowstomiss is the number of padrows the road can miss with respect
1059   //to the number of rows in the patch.
1060   //If the histograms were added, the comparison is done globally in the _slice_,
1061   //meaing that nrowstomiss is the number of padrows the road can miss with
1062   //respect to the total number of padrows in the slice.
1063   //
1064   //Return value = number of tracks which were removed (only in case of fAddHistograms)
1065
1066   if(!fTracks[0])
1067     {
1068       LOG(AliL3Log::kError,"AliL3Hough::Evaluate","Track Array")
1069         <<"No tracks to work with..."<<ENDLOG;
1070       return 0;
1071     }
1072
1073   Int_t removedtracks=0;
1074   AliL3TrackArray *tracks=0;
1075
1076   if(fAddHistograms)
1077     {
1078       tracks = fTracks[0];
1079       for(Int_t i=0; i<tracks->GetNTracks(); i++)
1080         {
1081           AliL3Track *track = tracks->GetCheckedTrack(i);
1082           if(!track) continue;
1083           track->SetNHits(0);
1084         }
1085     }
1086
1087   for(Int_t i=0; i<fNPatches; i++)
1088     EvaluatePatch(i,roadwidth,nrowstomiss);
1089
1090   //Here we check the tracks globally;
1091   //how many good rows (padrows with signal)
1092   //did it cross in the slice
1093   if(fAddHistograms)
1094     {
1095       for(Int_t j=0; j<tracks->GetNTracks(); j++)
1096         {
1097           AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
1098
1099           if(track->GetNHits() < AliL3Transform::GetNRows() - nrowstomiss)
1100             {
1101               tracks->Remove(j);
1102               removedtracks++;
1103             }
1104         }
1105       tracks->Compress();
1106       tracks->QSort();
1107     }
1108
1109   return removedtracks;
1110 }
1111
1112 void AliL3Hough::EvaluatePatch(Int_t i,Int_t roadwidth,Int_t nrowstomiss)
1113 {
1114   //Evaluate patch i.
1115
1116   fEval[i]->InitTransformer(fHoughTransformer[i]);
1117   fEval[i]->SetNumOfPadsToLook(roadwidth);
1118   fEval[i]->SetNumOfRowsToMiss(nrowstomiss);
1119   //fEval[i]->RemoveFoundTracks();
1120
1121   AliL3TrackArray *tracks=0;
1122
1123   if(!fAddHistograms)
1124     tracks = fTracks[i];
1125   else
1126     tracks = fTracks[0];
1127
1128   Int_t nrows=0;
1129   for(Int_t j=0; j<tracks->GetNTracks(); j++)
1130     {
1131       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
1132       if(!track)
1133         {
1134           LOG(AliL3Log::kWarning,"AliL3Hough::EvaluatePatch","Track array")
1135             <<"Track object missing!"<<ENDLOG;
1136           continue;
1137         }
1138       nrows=0;
1139       Int_t rowrange[2] = {AliL3Transform::GetFirstRow(i),AliL3Transform::GetLastRow(i)};
1140       Bool_t result = fEval[i]->LookInsideRoad(track,nrows,rowrange);
1141       if(fAddHistograms)
1142         {
1143           Int_t pre=track->GetNHits();
1144           track->SetNHits(pre+nrows);
1145         }
1146       else//the track crossed too few good padrows (padrows with signal) in the patch, so remove it
1147         {
1148           if(result == kFALSE)
1149             tracks->Remove(j);
1150         }
1151     }
1152
1153   tracks->Compress();
1154
1155 }
1156
1157 void AliL3Hough::MergeEtaSlices()
1158 {
1159   //Merge tracks found in neighbouring eta slices.
1160   //Removes the track with the lower weight.
1161
1162   fBenchmark->Start("Merge Eta-slices");
1163   AliL3TrackArray *tracks = fTracks[0];
1164   if(!tracks)
1165     {
1166       cerr<<"AliL3Hough::MergeEtaSlices : No tracks "<<endl;
1167       return;
1168     }
1169   for(Int_t j=0; j<tracks->GetNTracks(); j++)
1170     {
1171       AliL3HoughTrack *track1 = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
1172       if(!track1) continue;
1173       for(Int_t k=j+1; k<tracks->GetNTracks(); k++)
1174         {
1175           AliL3HoughTrack *track2 = (AliL3HoughTrack*)tracks->GetCheckedTrack(k);
1176           if(!track2) continue;
1177           if(abs(track1->GetEtaIndex() - track2->GetEtaIndex()) != 1) continue;
1178           if(fabs(track1->GetKappa()-track2->GetKappa()) < 0.006 &&
1179              fabs(track1->GetPsi()- track2->GetPsi()) < 0.1)
1180             {
1181               //cout<<"Merging track in slices "<<track1->GetEtaIndex()<<" "<<track2->GetEtaIndex()<<endl;
1182               if(track1->GetWeight() > track2->GetWeight())
1183                 tracks->Remove(k);
1184               else
1185                 tracks->Remove(j);
1186             }
1187         }
1188     }
1189   fBenchmark->Stop("Merge Eta-slices");
1190   tracks->Compress();
1191 }
1192
1193 void AliL3Hough::WriteTracks(Char_t *path)
1194 {
1195   // Write found tracks into file
1196   //cout<<"AliL3Hough::WriteTracks : Sorting the tracsk"<<endl;
1197   //fGlobalTracks->QSort();
1198
1199   Char_t filename[1024];
1200   sprintf(filename,"%s/tracks_%d.raw",path,fEvent);
1201   AliL3MemHandler mem;
1202   mem.SetBinaryOutput(filename);
1203   mem.TrackArray2Binary(fGlobalTracks);
1204   mem.CloseBinaryOutput();
1205   fGlobalTracks->Reset();
1206 }
1207
1208 void AliL3Hough::WriteTracks(Int_t slice,Char_t *path)
1209 {
1210   // Write found tracks slice by slice into file
1211
1212   AliL3MemHandler mem;
1213   Char_t fname[100];
1214   if(fAddHistograms)
1215     {
1216       sprintf(fname,"%s/tracks_ho_%d_%d.raw",path,fEvent,slice);
1217       mem.SetBinaryOutput(fname);
1218       mem.TrackArray2Binary(fTracks[0]);
1219       mem.CloseBinaryOutput();
1220     }
1221   else
1222     {
1223       for(Int_t i=0; i<fNPatches; i++)
1224         {
1225           sprintf(fname,"%s/tracks_ho_%d_%d_%d.raw",path,fEvent,slice,i);
1226           mem.SetBinaryOutput(fname);
1227           mem.TrackArray2Binary(fTracks[i]);
1228           mem.CloseBinaryOutput();
1229         }
1230     }
1231 }
1232
1233 #ifdef use_aliroot
1234 Int_t AliL3Hough::FillESD(AliESD *esd)
1235 {
1236   if(!fGlobalTracks) return 0;
1237   Int_t nglobaltracks = 0;
1238   for(Int_t i=0; i<fGlobalTracks->GetNTracks(); i++)
1239     {
1240       AliL3HoughTrack *tpt = (AliL3HoughTrack *)fGlobalTracks->GetCheckedTrack(i);
1241       if(!tpt) continue;
1242
1243       AliESDHLTtrack *esdtrack = new AliESDHLTtrack();
1244
1245       esdtrack->SetRowRange(tpt->GetFirstRow(),tpt->GetLastRow());
1246       esdtrack->SetNHits(tpt->GetNHits());
1247       esdtrack->SetFirstPoint(tpt->GetFirstPointX(),tpt->GetFirstPointY(),tpt->GetFirstPointZ());
1248       esdtrack->SetLastPoint(tpt->GetLastPointX(),tpt->GetLastPointY(),tpt->GetLastPointZ());
1249       esdtrack->SetPt(tpt->GetPt());
1250       esdtrack->SetPsi(tpt->GetPsi());
1251       esdtrack->SetTgl(tpt->GetTgl());
1252       esdtrack->SetCharge(tpt->GetCharge());
1253       esdtrack->SetMCid(tpt->GetMCid());
1254       esdtrack->SetWeight(tpt->GetWeight());
1255       esdtrack->SetSector(tpt->GetSector());
1256       esdtrack->SetBinXY(tpt->GetBinX(),tpt->GetBinY(),tpt->GetSizeX(),tpt->GetSizeY());
1257       esdtrack->SetPID(tpt->GetPID());
1258       esdtrack->ComesFromMainVertex(tpt->ComesFromMainVertex());
1259
1260       esd->AddHLTHoughTrack(esdtrack);
1261       nglobaltracks++;
1262       delete esdtrack;
1263     }
1264   return nglobaltracks;
1265 }
1266 #endif
1267
1268 void AliL3Hough::WriteDigits(Char_t *outfile)
1269 {
1270   //Write the current data to a new rootfile.
1271 #ifdef use_aliroot
1272
1273   for(Int_t i=0; i<fNPatches; i++)
1274     {
1275       AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer[i]->GetDataPointer();
1276       fMemHandler[i]->AliDigits2RootFile(tempPt,outfile);
1277     }
1278 #else
1279   cerr<<"AliL3Hough::WriteDigits : You need to compile with AliROOT!"<<endl;
1280   return;
1281 #endif
1282 }
1283
1284 Double_t AliL3Hough::GetCpuTime()
1285 {
1286   //Return the Cputime in seconds.
1287  struct timeval tv;
1288  gettimeofday( &tv, NULL );
1289  return tv.tv_sec+(((Double_t)tv.tv_usec)/1000000.);
1290 }
1291
1292 void *AliL3Hough::ProcessInThread(void *args)
1293 {
1294   AliL3Hough *instance = (AliL3Hough *)args;
1295   Int_t minslice = instance->GetMinSlice();
1296   Int_t maxslice = instance->GetMaxSlice();
1297   for(Int_t i=minslice; i<=maxslice; i++)
1298     {
1299       instance->ReadData(i,0);
1300       instance->Transform();
1301       instance->AddAllHistogramsRows();
1302       instance->FindTrackCandidatesRow();
1303       instance->AddTracks();
1304     }
1305   return (void *)0;
1306 }
1307
1308 void AliL3Hough::StartProcessInThread(Int_t minslice,Int_t maxslice)
1309 {
1310   if(!fThread) {
1311     char buf[255];
1312     sprintf(buf,"houghtrans_%d_%d",minslice,maxslice);
1313     SetMinMaxSlices(minslice,maxslice);
1314     //    fThread = new TThread(buf,(void (*) (void *))&ProcessInThread,(void *)this);
1315     fThread = new TThread(buf,&ProcessInThread,(void *)this);
1316     fThread->Run();
1317   }
1318   return;
1319 }
1320
1321 Int_t AliL3Hough::WaitForThreadFinish()
1322 {
1323   return TThread::Join(fThread->GetId());
1324 }