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 #if __GNUC__ == 3
  36 using namespace std;
  37 #endif
  38
  39 /** /class AliL3Hough
  40 //<pre>
  41 //_____________________________________________________________
  42 // AliL3Hough
  43 //
  44 // Interface class for the Hough transform
  45 //
  46 // Example how to use:
  47 //
  48 // AliL3Hough *hough = new AliL3Hough(path,kTRUE,NumberOfEtaSegments);
  49 // hough->ReadData(slice);
  50 // hough->Transform();
  51 // hough->FindTrackCandidates();
  52 //
  53 // AliL3TrackArray *tracks = hough->GetTracks(patch);
  54 //
  55 //</pre>
  56 */
  57
  58 ClassImp(AliL3Hough)
  59
  60 AliL3Hough::AliL3Hough()
  61 {
  62   //Constructor
  63
  64   fBinary        = kFALSE;
  65   fAddHistograms = kFALSE;
  66   fDoIterative   = kFALSE;
  67   fWriteDigits   = kFALSE;
  68   fUse8bits      = kFALSE;
  69
  70   fMemHandler       = 0;
  71   fHoughTransformer = 0;
  72   fEval             = 0;
  73   fPeakFinder       = 0;
  74   fTracks           = 0;
  75   fGlobalTracks     = 0;
  76   fMerger           = 0;
  77   fInterMerger      = 0;
  78   fGlobalMerger     = 0;
  79   fBenchmark        = 0;
  80
  81   fNEtaSegments     = 0;
  82   fNPatches         = 0;
  83   fVersion          = 0;
  84   fCurrentSlice     = 0;
  85   fEvent            = 0;
  86
  87   fKappaSpread = 6;
  88   fPeakRatio   = 0.5;
  89   fInputFile   = 0;
  90   fInputPtr    = 0;
  91
  92   SetTransformerParams();
  93   SetThreshold();
  94   SetNSaveIterations();
  95   SetPeakThreshold();
  96 #ifdef use_aliroot
  97   //just be sure that index is empty for new event
  98     AliL3FileHandler::CleanStaticIndex();
  99 #ifdef use_newio
 100     fRunLoader = 0;
 101 #endif
 102 #endif
 103 }
 104
 105 AliL3Hough::AliL3Hough(Char_t *path,Bool_t binary,Int_t n_eta_segments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr)
 106 {
 107   fBinary = binary;
 108   strcpy(fPath,path);
 109   fNEtaSegments  = n_eta_segments;
 110   fAddHistograms = kFALSE;
 111   fDoIterative   = kFALSE;
 112   fWriteDigits   = kFALSE;
 113   fUse8bits      = bit8;
 114   fVersion       = tv;
 115   fKappaSpread=6;
 116   fPeakRatio=0.5;
 117   if(!fBinary) {
 118     if(infile) {
 119       fInputFile = infile;
 120       fInputPtr = 0;
 121     }
 122     else {
 123       fInputFile = 0;
 124       fInputPtr = ptr;
 125     }
 126   }
 127   else {
 128     fInputFile = 0;
 129     fInputPtr = 0;
 130   }
 131 #ifdef use_aliroot
 132   //just be sure that index is empty for new event
 133     AliL3FileHandler::CleanStaticIndex();
 134 #ifdef use_newio
 135     fRunLoader = 0;
 136 #endif
 137 #endif
 138 }
 139
 140 AliL3Hough::~AliL3Hough()
 141 {
 142   //dtor
 143
 144   CleanUp();
 145   if(fMerger)
 146     delete fMerger;
 147   //cout << "Cleaned class merger " << endl;
 148   if(fInterMerger)
 149     delete fInterMerger;
 150   //cout << "Cleaned class inter " << endl;
 151   if(fPeakFinder)
 152     delete fPeakFinder;
 153   //cout << "Cleaned class peak " << endl;
 154   if(fGlobalMerger)
 155     delete fGlobalMerger;
 156   //cout << "Cleaned class global " << endl;
 157   if(fBenchmark)
 158     delete fBenchmark;
 159   //cout << "Cleaned class bench " << endl;
 160   if(fGlobalTracks)
 161     delete fGlobalTracks;
 162   //cout << "Cleaned class globaltracks " << endl;
 163 }
 164
 165 void AliL3Hough::CleanUp()
 166 {
 167   //Cleanup memory
 168
 169   for(Int_t i=0; i<fNPatches; i++)
 170     {
 171       if(fTracks[i]) delete fTracks[i];
 172       //cout << "Cleaned tracks " << i << endl;
 173       if(fEval[i]) delete fEval[i];
 174       //cout << "Cleaned eval " << i << endl;
 175       if(fHoughTransformer[i]) delete fHoughTransformer[i];
 176       //cout << "Cleaned traf " << i << endl;
 177       if(fMemHandler[i]) delete fMemHandler[i];
 178       //cout << "Cleaned mem " << i << endl;
 179     }
 180
 181   if(fTracks) delete [] fTracks;
 182   //cout << "Cleaned class tracks " << endl;
 183   if(fEval) delete [] fEval;
 184   //cout << "Cleaned class eval " << endl;
 185   if(fHoughTransformer) delete [] fHoughTransformer;
 186   //cout << "Cleaned cleass trafo " << endl;
 187   if(fMemHandler) delete [] fMemHandler;
 188   //cout << "Cleaned class mem " << endl;
 189 }
 190
 191 void AliL3Hough::Init(Char_t *path,Bool_t binary,Int_t n_eta_segments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr,Float_t zvertex)
 192 {
 193   fBinary = binary;
 194   strcpy(fPath,path);
 195   fNEtaSegments = n_eta_segments;
 196   fWriteDigits  = kFALSE;
 197   fUse8bits     = bit8;
 198   fVersion      = tv;
 199   if(!fBinary) {
 200     if(infile) {
 201       fInputFile = infile;
 202       fInputPtr = 0;
 203     }
 204     else {
 205       fInputFile = 0;
 206       fInputPtr = ptr;
 207     }
 208   }
 209   else {
 210     fInputFile = 0;
 211     fInputPtr = 0;
 212   }
 213   fZVertex = zvertex;
 214
 215   Init(); //do the rest
 216 }
 217
 218 void AliL3Hough::Init(Bool_t doit, Bool_t addhists)
 219 {
 220   fDoIterative   = doit;
 221   fAddHistograms = addhists;
 222
 223   fNPatches = AliL3Transform::GetNPatches();
 224
 225   fHoughTransformer = new AliL3HoughBaseTransformer*[fNPatches];
 226   fMemHandler = new AliL3MemHandler*[fNPatches];
 227
 228   fTracks = new AliL3TrackArray*[fNPatches];
 229   fEval = new AliL3HoughEval*[fNPatches];
 230
 231   fGlobalTracks = new AliL3TrackArray("AliL3HoughTrack");
 232
 233   for(Int_t i=0; i<fNPatches; i++)
 234     {
 235       switch (fVersion){ //choose Transformer
 236       case 1:
 237         fHoughTransformer[i] = new AliL3HoughTransformerLUT(0,i,fNEtaSegments);
 238         break;
 239       case 2:
 240         fHoughTransformer[i] = new AliL3HoughClusterTransformer(0,i,fNEtaSegments);
 241         break;
 242       case 3:
 243         fHoughTransformer[i] = new AliL3HoughTransformerVhdl(0,i,fNEtaSegments,fNSaveIterations);
 244         break;
 245       case 4:
 246         fHoughTransformer[i] = new AliL3HoughTransformerRow(0,i,fNEtaSegments,kFALSE,fZVertex);
 247         break;
 248       default:
 249         fHoughTransformer[i] = new AliL3HoughTransformer(0,i,fNEtaSegments,kFALSE,kFALSE);
 250       }
 251
 252       //      fHoughTransformer[i]->CreateHistograms(fNBinX[i],fLowPt[i],fNBinY[i],-fPhi[i],fPhi[i]);
 253       fHoughTransformer[i]->CreateHistograms(fNBinX[i],-fLowPt[i],fLowPt[i],fNBinY[i],-fPhi[i],fPhi[i]);
 254       //fHoughTransformer[i]->CreateHistograms(fLowPt[i],fUpperPt[i],fPtRes[i],fNBinY[i],fPhi[i]);
 255
 256       fHoughTransformer[i]->SetLowerThreshold(fThreshold[i]);
 257       fHoughTransformer[i]->SetUpperThreshold(100);
 258
 259       LOG(AliL3Log::kInformational,"AliL3Hough::Init","Version")
 260         <<"Initializing Hough transformer version "<<fVersion<<ENDLOG;
 261
 262       fEval[i] = new AliL3HoughEval();
 263       fTracks[i] = new AliL3TrackArray("AliL3HoughTrack");
 264       if(fUse8bits)
 265         fMemHandler[i] = new AliL3DataHandler();
 266       else
 267 #ifdef use_aliroot
 268         {
 269           if(!fInputFile) {
 270             if(!fInputPtr) {
 271               /* In case of reading digits file */
 272               fMemHandler[i] = new AliL3FileHandler(kTRUE); //use static index
 273               if(!fBinary) {
 274 #if use_newio
 275                 if(!fRunLoader) {
 276 #endif
 277                   Char_t filename[1024];
 278                   sprintf(filename,"%s/digitfile.root",fPath);
 279                   fMemHandler[i]->SetAliInput(filename);
 280 #if use_newio
 281                 }
 282                 else {
 283                   fMemHandler[i]->SetAliInput(fRunLoader);
 284                 }
 285 #endif
 286               }
 287             }
 288             else {
 289               /* In case of reading from DATE */
 290               fMemHandler[i] = new AliL3DDLDataFileHandler();
 291               fMemHandler[i]->SetReaderInput(fInputPtr,-1);
 292             }
 293           }
 294           else {
 295             /* In case of reading rawdata from ROOT file */
 296             fMemHandler[i] = new AliL3DDLDataFileHandler();
 297             fMemHandler[i]->SetReaderInput(fInputFile);
 298           }
 299         }
 300 #else
 301       fMemHandler[i] = new AliL3MemHandler();
 302 #endif
 303     }
 304
 305   fPeakFinder = new AliL3HoughMaxFinder("KappaPhi",50000);
 306   fMerger = new AliL3HoughMerger(fNPatches);
 307   fInterMerger = new AliL3HoughIntMerger();
 308   fGlobalMerger = 0;
 309   fBenchmark = new AliL3Benchmark();
 310 }
 311
 312 void AliL3Hough::SetTransformerParams(Float_t ptres,Float_t ptmin,Float_t ptmax,Int_t ny,Int_t patch)
 313 {
 314
 315   Int_t mrow;
 316   Float_t psi=0;
 317   if(patch==-1)
 318     mrow = 80;
 319   else
 320     mrow = AliL3Transform::GetLastRow(patch);
 321   if(ptmin)
 322     {
 323       Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(mrow),2) + pow(AliL3Transform::GetMaxY(mrow),2));
 324       Double_t kappa = -1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/ptmin;
 325       psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
 326       cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
 327     }
 328
 329   if(patch==-1)
 330     {
 331       Int_t i=0;
 332       while(i < 6)
 333         {
 334           fPtRes[i] = ptres;
 335           fLowPt[i] = ptmin;
 336           fUpperPt[i] = ptmax;
 337           fNBinY[i] = ny;
 338           fPhi[i] = psi;
 339           fNBinX[i]=0;
 340           i++;
 341         }
 342       return;
 343     }
 344
 345   fPtRes[patch] = ptres;
 346   fLowPt[patch] = ptmin;
 347   fUpperPt[patch] = ptmax;
 348   fNBinY[patch] = ny;
 349   fPhi[patch] = psi;
 350 }
 351 /*
 352 void AliL3Hough::SetTransformerParams(Int_t nx,Int_t ny,Float_t ptmin,Int_t patch)
 353 {
 354
 355   Int_t mrow=80;
 356   Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(mrow),2) + pow(AliL3Transform::GetMaxY(mrow),2));
 357   Double_t kappa = -1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/ptmin;
 358   Double_t psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
 359   cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
 360
 361   Int_t i=0;
 362   while(i < 6)
 363     {
 364       fLowPt[i] = ptmin;
 365       fNBinY[i] = ny;
 366       fNBinX[i] = nx;
 367       fPhi[i] = psi;
 368       i++;
 369     }
 370 }
 371 */
 372 void AliL3Hough::SetTransformerParams(Int_t nx,Int_t ny,Float_t ptmin,Int_t /*patch*/)
 373 {
 374
 375
 376   Int_t mrow=79;
 377   Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(mrow),2) + pow(AliL3Transform::GetMaxY(mrow),2));
 378   Double_t alpha1 = AliL3Transform::GetMaxY(mrow)/pow(lineradius,2);
 379   Double_t kappa = 1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/ptmin;
 380   Double_t psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
 381   //  cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
 382   AliL3HoughTrack track;
 383   track.SetTrackParameters(kappa,psi,1);
 384   Float_t hit[3];
 385   Int_t mrow2 = 158;
 386   track.GetCrossingPoint(mrow2,hit);
 387   Double_t lineradius2 = sqrt(pow(AliL3Transform::Row2X(mrow2),2) + pow(AliL3Transform::GetMaxY(mrow2),2));
 388   Double_t alpha2 = hit[1]/pow(lineradius2,2);
 389   //  cout<<"Calculated alphas range "<<alpha1<<" "<<alpha2<<" in patch "<<patch<<endl;
 390
 391   Int_t i=0;
 392   while(i < 6)
 393     {
 394       fLowPt[i] = 1.15*alpha1;
 395       fNBinY[i] = ny;
 396       fNBinX[i] = nx;
 397       fPhi[i] = 1.15*alpha2;
 398       i++;
 399     }
 400 }
 401
 402 void AliL3Hough::SetTransformerParams(Int_t nx,Int_t ny,Float_t lpt,Float_t phi)
 403 {
 404   Int_t i=0;
 405   while(i < 6)
 406     {
 407       fLowPt[i] = lpt;
 408       fNBinY[i] = ny;
 409       fNBinX[i] = nx;
 410       fPhi[i] = phi;
 411       i++;
 412     }
 413 }
 414
 415 void AliL3Hough::SetThreshold(Int_t t3,Int_t patch)
 416 {
 417   if(patch==-1)
 418     {
 419       Int_t i=0;
 420       while(i < 6)
 421         fThreshold[i++]=t3;
 422       return;
 423     }
 424   fThreshold[patch]=t3;
 425 }
 426
 427 void AliL3Hough::SetPeakThreshold(Int_t threshold,Int_t patch)
 428 {
 429   if(patch==-1)
 430     {
 431       Int_t i=0;
 432       while(i < 6)
 433         fPeakThreshold[i++]=threshold;
 434       return;
 435     }
 436   fPeakThreshold[patch]=threshold;
 437 }
 438
 439 void AliL3Hough::DoBench(Char_t *name)
 440 {
 441   fBenchmark->Analyze(name);
 442 }
 443
 444 void AliL3Hough::Process(Int_t minslice,Int_t maxslice)
 445 {
 446   //Process all slices [minslice,maxslice].
 447   fGlobalMerger = new AliL3HoughGlobalMerger(minslice,maxslice);
 448
 449   for(Int_t i=minslice; i<=maxslice; i++)
 450     {
 451       ReadData(i);
 452       Transform();
 453       if(fAddHistograms) {
 454         if(fVersion != 4)
 455           AddAllHistograms();
 456         else
 457           AddAllHistogramsRows();
 458       }
 459       FindTrackCandidates();
 460       //Evaluate();
 461       //fGlobalMerger->FillTracks(fTracks[0],i);
 462     }
 463 }
 464
 465 void AliL3Hough::ReadData(Int_t slice,Int_t eventnr)
 466 {
 467   //Read data from files, binary or root.
 468
 469 #ifdef use_aliroot
 470   if(fEvent!=eventnr) //just be sure that index is empty for new event
 471     AliL3FileHandler::CleanStaticIndex();
 472 #endif
 473   fCurrentSlice = slice;
 474
 475   for(Int_t i=0; i<fNPatches; i++)
 476     {
 477       fMemHandler[i]->Free();
 478       UInt_t ndigits=0;
 479       AliL3DigitRowData *digits =0;
 480       Char_t name[256];
 481       fMemHandler[i]->Init(slice,i);
 482       if(fBinary)//take input data from binary files
 483         {
 484           if(fUse8bits)
 485             sprintf(name,"%s/binaries/digits_c8_%d_%d_%d.raw",fPath,eventnr,slice,i);
 486           else
 487             sprintf(name,"%s/binaries/digits_%d_%d_%d.raw",fPath,eventnr,slice,i);
 488
 489           fMemHandler[i]->SetBinaryInput(name);
 490           digits = (AliL3DigitRowData *)fMemHandler[i]->CompBinary2Memory(ndigits);
 491           fMemHandler[i]->CloseBinaryInput();
 492         }
 493       else //read data from root file
 494         {
 495 #ifdef use_aliroot
 496           if(fEvent!=eventnr)
 497             fMemHandler[i]->FreeDigitsTree();//or else the new event is not loaded
 498           digits=(AliL3DigitRowData *)fMemHandler[i]->AliAltroDigits2Memory(ndigits,eventnr);
 499 #else
 500           cerr<<"You cannot read from rootfile now"<<endl;
 501 #endif
 502         }
 503
 504       //set input data and init transformer
 505       fHoughTransformer[i]->SetInputData(ndigits,digits);
 506       fHoughTransformer[i]->Init(slice,i,fNEtaSegments);
 507     }
 508
 509   fEvent=eventnr;
 510 }
 511
 512 void AliL3Hough::Transform(Int_t *row_range)
 513 {
 514   //Transform all data given to the transformer within the given slice
 515   //(after ReadData(slice))
 516
 517   Double_t initTime,cpuTime;
 518   initTime = GetCpuTime();
 519   for(Int_t i=0; i<fNPatches; i++)
 520     {
 521       // In case of Row transformer reset the arrays only once
 522       if((fVersion != 4) || (i == 0))
 523         fHoughTransformer[i]->Reset();//Reset the histograms
 524       fBenchmark->Start("Hough Transform");
 525       if(!row_range)
 526         fHoughTransformer[i]->TransformCircle();
 527       else
 528         fHoughTransformer[i]->TransformCircleC(row_range,1);
 529       fBenchmark->Stop("Hough Transform");
 530     }
 531   cpuTime = GetCpuTime() - initTime;
 532   LOG(AliL3Log::kInformational,"AliL3Hough::Transform()","Timing")
 533     <<"Transform done in average per patch of "<<cpuTime*1000/fNPatches<<" ms"<<ENDLOG;
 534 }
 535
 536 void AliL3Hough::MergePatches()
 537 {
 538   if(fAddHistograms) //Nothing to merge here
 539     return;
 540   fMerger->MergePatches(kTRUE);
 541 }
 542
 543 void AliL3Hough::MergeInternally()
 544 {
 545   if(fAddHistograms)
 546     fInterMerger->FillTracks(fTracks[0]);
 547   else
 548     fInterMerger->FillTracks(fMerger->GetOutTracks());
 549
 550   fInterMerger->MMerge();
 551 }
 552
 553 void AliL3Hough::ProcessSliceIter()
 554 {
 555   //Process current slice (after ReadData(slice)) iteratively.
 556
 557   if(!fAddHistograms)
 558     {
 559       for(Int_t i=0; i<fNPatches; i++)
 560         {
 561           ProcessPatchIter(i);
 562           fMerger->FillTracks(fTracks[i],i); //Copy tracks to merger
 563         }
 564     }
 565   else
 566     {
 567       for(Int_t i=0; i<10; i++)
 568         {
 569           Transform();
 570           AddAllHistograms();
 571           InitEvaluate();
 572           AliL3HoughBaseTransformer *tr = fHoughTransformer[0];
 573           for(Int_t j=0; j<fNEtaSegments; j++)
 574             {
 575               AliL3Histogram *hist = tr->GetHistogram(j);
 576               if(hist->GetNEntries()==0) continue;
 577               fPeakFinder->Reset();
 578               fPeakFinder->SetHistogram(hist);
 579               fPeakFinder->FindAbsMaxima();
 580               AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[0]->NextTrack();
 581               track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
 582               track->SetEtaIndex(j);
 583               track->SetEta(tr->GetEta(j,fCurrentSlice));
 584               for(Int_t k=0; k<fNPatches; k++)
 585                 {
 586                   fEval[i]->SetNumOfPadsToLook(2);
 587                   fEval[i]->SetNumOfRowsToMiss(2);
 588                   fEval[i]->RemoveFoundTracks();
 589                   /*
 590                   Int_t nrows=0;
 591                   if(!fEval[i]->LookInsideRoad(track,nrows))
 592                     {
 593                       fTracks[0]->Remove(fTracks[0]->GetNTracks()-1);
 594                       fTracks[0]->Compress();
 595                     }
 596                   */
 597                 }
 598             }
 599
 600         }
 601
 602     }
 603 }
 604
 605 void AliL3Hough::ProcessPatchIter(Int_t patch)
 606 {
 607   //Process patch in a iterative way.
 608   //transform + peakfinding + evaluation + transform +...
 609
 610   Int_t num_of_tries = 5;
 611   AliL3HoughBaseTransformer *tr = fHoughTransformer[patch];
 612   AliL3TrackArray *tracks = fTracks[patch];
 613   tracks->Reset();
 614   AliL3HoughEval *ev = fEval[patch];
 615   ev->InitTransformer(tr);
 616   //ev->RemoveFoundTracks();
 617   ev->SetNumOfRowsToMiss(3);
 618   ev->SetNumOfPadsToLook(2);
 619   AliL3Histogram *hist;
 620   for(Int_t t=0; t<num_of_tries; t++)
 621     {
 622       tr->Reset();
 623       tr->TransformCircle();
 624       for(Int_t i=0; i<fNEtaSegments; i++)
 625         {
 626           hist = tr->GetHistogram(i);
 627           if(hist->GetNEntries()==0) continue;
 628           fPeakFinder->Reset();
 629           fPeakFinder->SetHistogram(hist);
 630           fPeakFinder->FindAbsMaxima();
 631           //fPeakFinder->FindPeak1();
 632           AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->NextTrack();
 633           track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
 634           track->SetEtaIndex(i);
 635           track->SetEta(tr->GetEta(i,fCurrentSlice));
 636           /*
 637           Int_t nrows=0;
 638           if(!ev->LookInsideRoad(track,nrows))
 639             {
 640               tracks->Remove(tracks->GetNTracks()-1);
 641               tracks->Compress();
 642             }
 643           */
 644         }
 645     }
 646   fTracks[0]->QSort();
 647   LOG(AliL3Log::kInformational,"AliL3Hough::ProcessPatch","NTracks")
 648     <<AliL3Log::kDec<<"Found "<<tracks->GetNTracks()<<" tracks in patch "<<patch<<ENDLOG;
 649 }
 650
 651 void AliL3Hough::AddAllHistograms()
 652 {
 653   //Add the histograms within one etaslice.
 654   //Resulting histogram are in patch=0.
 655
 656   Double_t initTime,cpuTime;
 657   initTime = GetCpuTime();
 658   fBenchmark->Start("Add Histograms");
 659   for(Int_t i=0; i<fNEtaSegments; i++)
 660     {
 661       AliL3Histogram *hist0 = fHoughTransformer[0]->GetHistogram(i);
 662       for(Int_t j=1; j<fNPatches; j++)
 663         {
 664           AliL3Histogram *hist = fHoughTransformer[j]->GetHistogram(i);
 665           hist0->Add(hist);
 666         }
 667     }
 668   fBenchmark->Stop("Add Histograms");
 669   fAddHistograms = kTRUE;
 670   cpuTime = GetCpuTime() - initTime;
 671   LOG(AliL3Log::kInformational,"AliL3Hough::AddAllHistograms()","Timing")
 672     <<"Adding histograms in "<<cpuTime*1000<<" ms"<<ENDLOG;
 673 }
 674
 675 void AliL3Hough::AddAllHistogramsRows()
 676 {
 677   //Add the histograms within one etaslice.
 678   //Resulting histogram are in patch=0.
 679
 680   Double_t initTime,cpuTime;
 681   initTime = GetCpuTime();
 682   fBenchmark->Start("Add HistogramsRows");
 683
 684   UChar_t *tracknrows = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetTrackNRows();
 685   UChar_t *trackfirstrow = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetTrackFirstRow();
 686   UChar_t *tracklastrow = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetTrackLastRow();
 687
 688   for(Int_t i=0; i<fNEtaSegments; i++)
 689     {
 690       UChar_t *rowcount = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetRowCount(i);
 691       UChar_t *gapcount = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetGapCount(i);
 692       UChar_t *currentrowcount = ((AliL3HoughTransformerRow *)fHoughTransformer[0])->GetCurrentRowCount(i);
 693
 694       AliL3Histogram *hist = fHoughTransformer[0]->GetHistogram(i);
 695       Int_t xmin = hist->GetFirstXbin();
 696       Int_t xmax = hist->GetLastXbin();
 697       Int_t ymin = hist->GetFirstYbin();
 698       Int_t ymax = hist->GetLastYbin();
 699
 700       for(Int_t ybin=ymin; ybin<=ymax; ybin++)
 701         {
 702           for(Int_t xbin=xmin; xbin<=xmax; xbin++)
 703             {
 704               Int_t bin = hist->GetBin(xbin,ybin);
 705               if(tracklastrow[bin] > (currentrowcount[bin] + 1))
 706                  gapcount[bin]++;
 707               if(gapcount[bin] < MAX_N_GAPS)
 708                 if(rowcount[bin] >= MIN_TRACK_LENGTH)
 709                   if(((Int_t)rowcount[bin] + (Int_t)gapcount[bin])>=((Int_t)tracknrows[bin]-MAX_MISS_ROWS))
 710                     hist->AddBinContent(bin,(rowcount[bin]+trackfirstrow[bin]+159-tracklastrow[bin]));
 711             }
 712         }
 713     }
 714
 715   fBenchmark->Stop("Add HistogramsRows");
 716   fAddHistograms = kTRUE;
 717   cpuTime = GetCpuTime() - initTime;
 718   LOG(AliL3Log::kInformational,"AliL3Hough::AddAllHistogramsRows()","Timing")
 719     <<"Adding histograms in "<<cpuTime*1000<<" ms"<<ENDLOG;
 720 }
 721
 722 void AliL3Hough::AddTracks()
 723 {
 724   if(!fTracks[0])
 725     {
 726       cerr<<"AliL3Hough::AddTracks : No tracks"<<endl;
 727       return;
 728     }
 729   AliL3TrackArray *tracks = fTracks[0];
 730   for(Int_t i=0; i<tracks->GetNTracks(); i++)
 731     {
 732       AliL3Track *track = tracks->GetCheckedTrack(i);
 733       if(!track) continue;
 734       if(track->GetNHits()!=1) cerr<<"NHITS "<<track->GetNHits()<<endl;
 735       UInt_t *ids = track->GetHitNumbers();
 736       ids[0] = (fCurrentSlice&0x7f)<<25;
 737     }
 738
 739   fGlobalTracks->AddTracks(fTracks[0],0,fCurrentSlice);
 740 }
 741
 742 void AliL3Hough::FindTrackCandidatesRow()
 743 {
 744   if(fVersion != 4) {
 745     LOG(AliL3Log::kError,"AliL3Hough::FindTrackCandidatesRow()","")
 746       <<"Incompatible Peak Finder version!"<<ENDLOG;
 747     return;
 748   }
 749
 750   //Look for peaks in histograms, and find the track candidates
 751   Int_t n_patches;
 752   if(fAddHistograms)
 753     n_patches = 1; //Histograms have been added.
 754   else
 755     n_patches = fNPatches;
 756
 757   Double_t initTime,cpuTime;
 758   initTime = GetCpuTime();
 759   fBenchmark->Start("Find Maxima");
 760   for(Int_t i=0; i<n_patches; i++)
 761     {
 762       AliL3HoughBaseTransformer *tr = fHoughTransformer[i];
 763       fTracks[i]->Reset();
 764       fPeakFinder->Reset();
 765
 766       for(Int_t j=0; j<fNEtaSegments; j++)
 767         {
 768           AliL3Histogram *hist = tr->GetHistogram(j);
 769           if(hist->GetNEntries()==0) continue;
 770           fPeakFinder->SetHistogram(hist);
 771           fPeakFinder->SetEtaSlice(j);
 772 #ifdef do_mc
 773           LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","")
 774             <<"Starting "<<j<<" etaslice"<<ENDLOG;
 775 #endif
 776           fPeakFinder->SetThreshold(fPeakThreshold[i]);
 777           fPeakFinder->FindAdaptedRowPeaks(1,0,0);//Maxima finder for HoughTransformerRow
 778
 779           //fPeakFinder->FindMaxima(fPeakThreshold[i]); //Simple maxima finder
 780         }
 781
 782       for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
 783         {
 784           if(fPeakFinder->GetWeight(k) < 0) continue;
 785           AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[i]->NextTrack();
 786           Float_t psi = atan((fPeakFinder->GetXPeak(k)-fPeakFinder->GetYPeak(k))/(AliL3HoughTransformerRow::GetBeta1()-AliL3HoughTransformerRow::GetBeta2()));
 787           Float_t kappa = 2.0*(fPeakFinder->GetXPeak(k)*cos(psi)-AliL3HoughTransformerRow::GetBeta1()*sin(psi));
 788           //          track->SetTrackParameters(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetWeight(k));
 789           track->SetTrackParameters(kappa,psi,fPeakFinder->GetWeight(k));
 790           track->SetBinXY(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetXPeakSize(k),fPeakFinder->GetYPeakSize(k));
 791           Int_t etaindex = (fPeakFinder->GetStartEta(k)+fPeakFinder->GetEndEta(k))/2;
 792           track->SetEtaIndex(etaindex);
 793           Float_t starteta = tr->GetEta(fPeakFinder->GetStartEta(k),fCurrentSlice);
 794           Float_t endeta = tr->GetEta(fPeakFinder->GetEndEta(k),fCurrentSlice);
 795           track->SetEta((starteta+endeta)/2.0);
 796           track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
 797           track->SetSector(fCurrentSlice);
 798           track->SetSlice(fCurrentSlice);
 799 #ifdef do_mc
 800           Int_t label = tr->GetTrackID(etaindex,fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k));
 801           track->SetMCid(label);
 802           //      cout<<"Track found with label "<<label<<" at "<<fPeakFinder->GetXPeak(k)<<" "<<fPeakFinder->GetYPeak(k)<<" with weight "<<fPeakFinder->GetWeight(k)<<endl;
 803 #endif
 804         }
 805       LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","")
 806         <<"Found "<<fTracks[i]->GetNTracks()<<" tracks in patch "<<i<<ENDLOG;
 807       fTracks[i]->QSort();
 808     }
 809   fBenchmark->Stop("Find Maxima");
 810   cpuTime = GetCpuTime() - initTime;
 811   LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","Timing")
 812     <<"Maxima finding done in "<<cpuTime*1000<<" ms"<<ENDLOG;
 813 }
 814
 815 void AliL3Hough::FindTrackCandidates()
 816 {
 817   if(fVersion == 4) {
 818     LOG(AliL3Log::kError,"AliL3Hough::FindTrackCandidatesRow()","")
 819       <<"Incompatible Peak Finder version!"<<ENDLOG;
 820     return;
 821   }
 822
 823   Int_t n_patches;
 824   if(fAddHistograms)
 825     n_patches = 1; //Histograms have been added.
 826   else
 827     n_patches = fNPatches;
 828
 829   Double_t initTime,cpuTime;
 830   initTime = GetCpuTime();
 831   fBenchmark->Start("Find Maxima");
 832   for(Int_t i=0; i<n_patches; i++)
 833     {
 834       AliL3HoughBaseTransformer *tr = fHoughTransformer[i];
 835       fTracks[i]->Reset();
 836
 837       for(Int_t j=0; j<fNEtaSegments; j++)
 838         {
 839           AliL3Histogram *hist = tr->GetHistogram(j);
 840           if(hist->GetNEntries()==0) continue;
 841           fPeakFinder->Reset();
 842           fPeakFinder->SetHistogram(hist);
 843 #ifdef do_mc
 844           cout<<"Starting "<<j<<" etaslice"<<endl;
 845 #endif
 846           fPeakFinder->SetThreshold(fPeakThreshold[i]);
 847           fPeakFinder->FindAdaptedPeaks(fKappaSpread,fPeakRatio);
 848
 849           for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
 850             {
 851               AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[i]->NextTrack();
 852               track->SetTrackParameters(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetWeight(k));
 853               track->SetEtaIndex(j);
 854               track->SetEta(tr->GetEta(j,fCurrentSlice));
 855               track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
 856             }
 857         }
 858       cout<<"Found "<<fTracks[i]->GetNTracks()<<" tracks in patch "<<i<<endl;
 859       fTracks[i]->QSort();
 860     }
 861   fBenchmark->Stop("Find Maxima");
 862   cpuTime = GetCpuTime() - initTime;
 863   LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","Timing")
 864     <<"Maxima finding done in "<<cpuTime*1000<<" ms"<<ENDLOG;
 865 }
 866
 867 void AliL3Hough::InitEvaluate()
 868 {
 869   //Pass the transformer objects to the AliL3HoughEval objects:
 870   //This will provide the evaluation objects with all the necessary
 871   //data and parameters it needs.
 872
 873   for(Int_t i=0; i<fNPatches; i++)
 874     fEval[i]->InitTransformer(fHoughTransformer[i]);
 875 }
 876
 877 Int_t AliL3Hough::Evaluate(Int_t road_width,Int_t nrowstomiss)
 878 {
 879   //Evaluate the tracks, by looking along the road in the raw data.
 880   //If track does not cross all padrows - rows2miss, it is removed from the arrray.
 881   //If histograms were not added, the check is done locally in patch,
 882   //meaning that nrowstomiss is the number of padrows the road can miss with respect
 883   //to the number of rows in the patch.
 884   //If the histograms were added, the comparison is done globally in the _slice_,
 885   //meaing that nrowstomiss is the number of padrows the road can miss with
 886   //respect to the total number of padrows in the slice.
 887   //
 888   //Return value = number of tracks which were removed (only in case of fAddHistograms)
 889
 890   if(!fTracks[0])
 891     {
 892       LOG(AliL3Log::kError,"AliL3Hough::Evaluate","Track Array")
 893         <<"No tracks to work with..."<<ENDLOG;
 894       return 0;
 895     }
 896
 897   Int_t removed_tracks=0;
 898   AliL3TrackArray *tracks=0;
 899
 900   if(fAddHistograms)
 901     {
 902       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           track->SetNHits(0);
 908         }
 909     }
 910
 911   for(Int_t i=0; i<fNPatches; i++)
 912     EvaluatePatch(i,road_width,nrowstomiss);
 913
 914   //Here we check the tracks globally;
 915   //how many good rows (padrows with signal)
 916   //did it cross in the slice
 917   if(fAddHistograms)
 918     {
 919       for(Int_t j=0; j<tracks->GetNTracks(); j++)
 920         {
 921           AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
 922
 923           if(track->GetNHits() < AliL3Transform::GetNRows() - nrowstomiss)
 924             {
 925               tracks->Remove(j);
 926               removed_tracks++;
 927             }
 928         }
 929       tracks->Compress();
 930       tracks->QSort();
 931     }
 932
 933   return removed_tracks;
 934 }
 935
 936 void AliL3Hough::EvaluatePatch(Int_t i,Int_t road_width,Int_t nrowstomiss)
 937 {
 938   //Evaluate patch i.
 939
 940   fEval[i]->InitTransformer(fHoughTransformer[i]);
 941   fEval[i]->SetNumOfPadsToLook(road_width);
 942   fEval[i]->SetNumOfRowsToMiss(nrowstomiss);
 943   //fEval[i]->RemoveFoundTracks();
 944
 945   AliL3TrackArray *tracks=0;
 946
 947   if(!fAddHistograms)
 948     tracks = fTracks[i];
 949   else
 950     tracks = fTracks[0];
 951
 952   Int_t nrows=0;
 953   for(Int_t j=0; j<tracks->GetNTracks(); j++)
 954     {
 955       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
 956       if(!track)
 957         {
 958           LOG(AliL3Log::kWarning,"AliL3Hough::EvaluatePatch","Track array")
 959             <<"Track object missing!"<<ENDLOG;
 960           continue;
 961         }
 962       nrows=0;
 963       Int_t rowrange[2] = {AliL3Transform::GetFirstRow(i),AliL3Transform::GetLastRow(i)};
 964       Bool_t result = fEval[i]->LookInsideRoad(track,nrows,rowrange);
 965       if(fAddHistograms)
 966         {
 967           Int_t pre=track->GetNHits();
 968           track->SetNHits(pre+nrows);
 969         }
 970       else//the track crossed too few good padrows (padrows with signal) in the patch, so remove it
 971         {
 972           if(result == kFALSE)
 973             tracks->Remove(j);
 974         }
 975     }
 976
 977   tracks->Compress();
 978
 979 }
 980
 981 void AliL3Hough::MergeEtaSlices()
 982 {
 983   //Merge tracks found in neighbouring eta slices.
 984   //Removes the track with the lower weight.
 985
 986   fBenchmark->Start("Merge Eta-slices");
 987   AliL3TrackArray *tracks = fTracks[0];
 988   if(!tracks)
 989     {
 990       cerr<<"AliL3Hough::MergeEtaSlices : No tracks "<<endl;
 991       return;
 992     }
 993   for(Int_t j=0; j<tracks->GetNTracks(); j++)
 994     {
 995       AliL3HoughTrack *track1 = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
 996       if(!track1) continue;
 997       for(Int_t k=j+1; k<tracks->GetNTracks(); k++)
 998         {
 999           AliL3HoughTrack *track2 = (AliL3HoughTrack*)tracks->GetCheckedTrack(k);
1000           if(!track2) continue;
1001           if(abs(track1->GetEtaIndex() - track2->GetEtaIndex()) != 1) continue;
1002           if(fabs(track1->GetKappa()-track2->GetKappa()) < 0.006 &&
1003              fabs(track1->GetPsi()- track2->GetPsi()) < 0.1)
1004             {
1005               //cout<<"Merging track in slices "<<track1->GetEtaIndex()<<" "<<track2->GetEtaIndex()<<endl;
1006               if(track1->GetWeight() > track2->GetWeight())
1007                 tracks->Remove(k);
1008               else
1009                 tracks->Remove(j);
1010             }
1011         }
1012     }
1013   fBenchmark->Stop("Merge Eta-slices");
1014   tracks->Compress();
1015 }
1016
1017 void AliL3Hough::WriteTracks(Char_t *path)
1018 {
1019   //cout<<"AliL3Hough::WriteTracks : Sorting the tracsk"<<endl;
1020   //fGlobalTracks->QSort();
1021
1022   Char_t filename[1024];
1023   sprintf(filename,"%s/tracks_%d.raw",path,fEvent);
1024   AliL3MemHandler mem;
1025   mem.SetBinaryOutput(filename);
1026   mem.TrackArray2Binary(fGlobalTracks);
1027   mem.CloseBinaryOutput();
1028   fGlobalTracks->Reset();
1029 }
1030
1031 void AliL3Hough::WriteTracks(Int_t slice,Char_t *path)
1032 {
1033
1034   AliL3MemHandler mem;
1035   Char_t fname[100];
1036   if(fAddHistograms)
1037     {
1038       sprintf(fname,"%s/tracks_ho_%d_%d.raw",path,fEvent,slice);
1039       mem.SetBinaryOutput(fname);
1040       mem.TrackArray2Binary(fTracks[0]);
1041       mem.CloseBinaryOutput();
1042     }
1043   else
1044     {
1045       for(Int_t i=0; i<fNPatches; i++)
1046         {
1047           sprintf(fname,"%s/tracks_ho_%d_%d_%d.raw",path,fEvent,slice,i);
1048           mem.SetBinaryOutput(fname);
1049           mem.TrackArray2Binary(fTracks[i]);
1050           mem.CloseBinaryOutput();
1051         }
1052     }
1053 }
1054
1055 void AliL3Hough::WriteDigits(Char_t *outfile)
1056 {
1057 #ifdef use_aliroot
1058   //Write the current data to a new rootfile.
1059
1060   for(Int_t i=0; i<fNPatches; i++)
1061     {
1062       AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer[i]->GetDataPointer();
1063       fMemHandler[i]->AliDigits2RootFile(tempPt,outfile);
1064     }
1065 #else
1066   cerr<<"AliL3Hough::WriteDigits : You need to compile with AliROOT!"<<endl;
1067   return;
1068 #endif
1069 }
1070
1071 Double_t AliL3Hough::GetCpuTime()
1072 {
1073   //Return the Cputime in seconds.
1074  struct timeval tv;
1075  gettimeofday( &tv, NULL );
1076  return tv.tv_sec+(((Double_t)tv.tv_usec)/1000000.);
1077 }
1078