HLT/hough/AliL3Hough.cxx

   1 //$Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ASV
   5
   6
   7 #include <string.h>
   8 #include <sys/time.h>
   9
  10 #include "AliL3HoughMerger.h"
  11 #include "AliL3HoughIntMerger.h"
  12 #include "AliL3HoughGlobalMerger.h"
  13 #include "AliL3Logging.h"
  14 #include "AliL3Histogram.h"
  15 #include "AliL3Hough.h"
  16 #include "AliL3HoughTransformer.h"
  17 #include "AliL3HoughTransformerVhdl.h"
  18 #include "AliL3HoughMaxFinder.h"
  19 #ifdef use_aliroot
  20 #include "AliL3FileHandler.h"
  21 #else
  22 #include "AliL3MemHandler.h"
  23 #endif
  24 #include "AliL3DataHandler.h"
  25 #include "AliL3DigitData.h"
  26 #include "AliL3HoughEval.h"
  27 #include "AliL3Transform.h"
  28 #include "AliL3TrackArray.h"
  29 #include "AliL3HoughTrack.h"
  30
  31
  32 /** /class AliL3Hough
  33 //<pre>
  34 //_____________________________________________________________
  35 // AliL3Hough
  36 //
  37 // Interface class for the Hough transform
  38 //
  39 // Example how to use:
  40 //
  41 // AliL3Hough *hough = new AliL3Hough(path,kTRUE,NumberOfEtaSegments);
  42 // hough->ReadData(slice);
  43 // hough->Transform();
  44 // hough->FindTrackCandidates();
  45 //
  46 // AliL3TrackArray *tracks = hough->GetTracks(patch);
  47 //</pre>
  48 */
  49
  50 ClassImp(AliL3Hough)
  51
  52 AliL3Hough::AliL3Hough()
  53 {
  54   //Constructor
  55
  56   fBinary        = kFALSE;
  57   fAddHistograms = kFALSE;
  58   fDoIterative   = kFALSE;
  59   fWriteDigits   = kFALSE;
  60   fUse8bits      = kFALSE;
  61
  62   fMemHandler       = 0;
  63   fHoughTransformer = 0;
  64   fEval             = 0;
  65   fPeakFinder       = 0;
  66   fTracks           = 0;
  67   fMerger           = 0;
  68   fInterMerger      = 0;
  69   fGlobalMerger     = 0;
  70
  71   fNEtaSegments     = 0;
  72   fNPatches         = 0;
  73   fVersion          = 0;
  74   fCurrentSlice     = 0;
  75
  76   SetTransformerParams();
  77   SetThreshold();
  78 }
  79
  80 AliL3Hough::AliL3Hough(Char_t *path,Bool_t binary,Int_t n_eta_segments=100,Bool_t bit8=kFALSE,Int_t tv=0)
  81 {
  82   //Default ctor.
  83
  84   fBinary = binary;
  85   strcpy(fPath,path);
  86   fNEtaSegments  = n_eta_segments;
  87   fAddHistograms = kFALSE;
  88   fDoIterative   = kFALSE;
  89   fWriteDigits   = kFALSE;
  90   fUse8bits      = bit8;
  91   fVersion       = tv;
  92 }
  93
  94 AliL3Hough::~AliL3Hough()
  95 {
  96   //dtor
  97
  98   CleanUp();
  99   if(fMerger)
 100     delete fMerger;
 101   if(fInterMerger)
 102     delete fInterMerger;
 103   if(fPeakFinder)
 104     delete fPeakFinder;
 105   if(fGlobalMerger)
 106     delete fGlobalMerger;
 107 }
 108
 109 void AliL3Hough::CleanUp()
 110 {
 111   //Cleanup memory
 112
 113   for(Int_t i=0; i<fNPatches; i++)
 114     {
 115       if(fTracks[i]) delete fTracks[i];
 116       if(fEval[i]) delete fEval[i];
 117       if(fHoughTransformer[i]) delete fHoughTransformer[i];
 118       if(fMemHandler[i]) delete fMemHandler[i];
 119     }
 120
 121   /*
 122   if(fTracks) delete [] fTracks;
 123   if(fEval) delete [] fEval;
 124   if(fHoughTransformer) delete [] fHoughTransformer;
 125   if(fMemHandler) delete [] fMemHandler;
 126   */
 127 }
 128
 129 void AliL3Hough::Init(Char_t *path,Bool_t binary,Int_t n_eta_segments=100,Bool_t bit8=kFALSE,Int_t tv=0)
 130 {
 131   fBinary = binary;
 132   strcpy(fPath,path);
 133   fNEtaSegments = n_eta_segments;
 134   fWriteDigits  = kFALSE;
 135   fUse8bits     = bit8;
 136   fVersion      = tv;
 137
 138   Init(); //do the rest
 139 }
 140
 141 void AliL3Hough::Init(Bool_t doit=kFALSE, Bool_t addhists=kFALSE)
 142 {
 143   fDoIterative   = doit;
 144   fAddHistograms = addhists;
 145
 146   AliL3Transform::Init(fPath);
 147   fNPatches = AliL3Transform::GetNPatches();
 148
 149   fHoughTransformer = new AliL3HoughBaseTransformer*[fNPatches];
 150   fMemHandler = new AliL3MemHandler*[fNPatches];
 151   fTracks = new AliL3TrackArray*[fNPatches];
 152   fEval = new AliL3HoughEval*[fNPatches];
 153
 154   for(Int_t i=0; i<fNPatches; i++)
 155     {
 156       //VESTBO: Changes here, slice information is in principle not needed
 157       //by the HoughTransformers. You always gave here a one, I changed
 158       //it to zero and have Init function with is called in ReadData (see
 159       //a few lines below.). I did not have time to think about it longer
 160       //as there was that strange bug, but we should either consider not
 161       //to give slice at all or to give it correcty (thats what I am doing
 162       //a the moment.)
 163
 164       switch (fVersion){ //choose Transformer
 165       case 1:
 166         fHoughTransformer[i] = new AliL3HoughTransformerVhdl(0,i,fNEtaSegments);
 167         break;
 168       default:
 169         fHoughTransformer[i] = new AliL3HoughTransformer(0,i,fNEtaSegments);
 170       }
 171
 172       fHoughTransformer[i]->CreateHistograms(fNBinX,fLowPt,fNBinY,-fPhi,fPhi);
 173       fHoughTransformer[i]->SetLowerThreshold(fThreshold);
 174
 175       LOG(AliL3Log::kInformational,"AliL3Hough::Init","Version")
 176         <<"Initializing Hough transformer version "<<fVersion<<ENDLOG;
 177
 178       fEval[i] = new AliL3HoughEval();
 179       fTracks[i] = new AliL3TrackArray("AliL3HoughTrack");
 180       if(fUse8bits)
 181         fMemHandler[i] = new AliL3DataHandler();
 182       else
 183 #ifdef use_aliroot
 184         {
 185           fMemHandler[i] = new AliL3FileHandler();
 186           if(!fBinary)
 187             {
 188               Char_t filename[100];
 189               sprintf(filename,"%s/digitfile",fPath);
 190               fMemHandler[i]->SetAliInput(filename);
 191             }
 192         }
 193 #else
 194       fMemHandler[i] = new AliL3MemHandler();
 195 #endif
 196     }
 197
 198   fPeakFinder = new AliL3HoughMaxFinder("KappaPhi",1000);
 199   fMerger = new AliL3HoughMerger(fNPatches);
 200   fInterMerger = new AliL3HoughIntMerger();
 201   fGlobalMerger = 0;
 202 }
 203
 204 void AliL3Hough::Process(Int_t minslice,Int_t maxslice)
 205 {
 206   //Process all slices [minslice,maxslice].
 207   fGlobalMerger = new AliL3HoughGlobalMerger(minslice,maxslice);
 208
 209   for(Int_t i=minslice; i<=maxslice; i++)
 210     {
 211       ReadData(i);
 212       Transform();
 213       if(fAddHistograms)
 214         AddAllHistograms();
 215       FindTrackCandidates();
 216       Evaluate();
 217       fGlobalMerger->FillTracks(fTracks[0],i);
 218     }
 219 }
 220
 221 void AliL3Hough::ReadData(Int_t slice,Int_t eventnr=0)
 222 {
 223   //Read data from files, binary or root.
 224
 225   fCurrentSlice = slice;
 226   for(Int_t i=0; i<fNPatches; i++)
 227     {
 228       fMemHandler[i]->Free();
 229       UInt_t ndigits=0;
 230       AliL3DigitRowData *digits =0;
 231       Char_t name[256];
 232       fMemHandler[i]->Init(slice,i);
 233       if(fBinary)//take input data from binary files
 234         {
 235           if(fUse8bits)
 236             sprintf(name,"%sdigits_c8_%d_%d.raw",fPath,slice,i);
 237           else
 238             sprintf(name,"%sdigits_%d_%d.raw",fPath,slice,i);
 239           fMemHandler[i]->SetBinaryInput(name);
 240           digits = (AliL3DigitRowData *)fMemHandler[i]->CompBinary2Memory(ndigits);
 241           fMemHandler[i]->CloseBinaryInput();
 242         }
 243       else //read data from root file
 244         {
 245 #ifdef use_aliroot
 246           digits=(AliL3DigitRowData *)fMemHandler[i]->AliDigits2Memory(ndigits,eventnr);
 247           fMemHandler[i]->FreeDigitsTree();
 248 #else
 249           cerr<<"You cannot read from rootfile now"<<endl;
 250 #endif
 251         }
 252       fHoughTransformer[i]->SetInputData(ndigits,digits);
 253
 254       //VESTBO: due to bug commented out!!!
 255       //fHoughTransformer[i]->Init(slice,i,fNEtaSegments);
 256       //fHoughTransformer[i]->Print();
 257     }
 258 }
 259
 260 void AliL3Hough::Transform(Int_t row_range=-1)
 261 {
 262   //Transform all data given to the transformer within the given slice
 263   //(after ReadData(slice))
 264
 265   Double_t initTime,cpuTime;
 266   initTime = GetCpuTime();
 267   for(Int_t i=0; i<fNPatches; i++)
 268     {
 269       fHoughTransformer[i]->Reset();//Reset the histograms
 270       if(row_range < 0)
 271         fHoughTransformer[i]->TransformCircle();
 272       else
 273         fHoughTransformer[i]->TransformCircleC(row_range);
 274     }
 275   cpuTime = GetCpuTime() - initTime;
 276   LOG(AliL3Log::kInformational,"AliL3Hough::Transform()","Timing")
 277     <<"Transform done in average per patch of "<<cpuTime*1000/fNPatches<<" ms"<<ENDLOG;
 278 }
 279
 280 void AliL3Hough::MergePatches()
 281 {
 282   if(fAddHistograms) //Nothing to merge here
 283     return;
 284   fMerger->MergePatches(kTRUE);
 285 }
 286
 287 void AliL3Hough::MergeInternally()
 288 {
 289   if(fAddHistograms)
 290     fInterMerger->FillTracks(fTracks[0]);
 291   else
 292     fInterMerger->FillTracks(fMerger->GetOutTracks());
 293
 294   fInterMerger->MMerge();
 295 }
 296
 297 void AliL3Hough::ProcessSliceIter()
 298 {
 299   //Process current slice (after ReadData(slice)) iteratively.
 300
 301   for(Int_t i=0; i<fNPatches; i++)
 302     {
 303       ProcessPatchIter(i);
 304       fMerger->FillTracks(fTracks[i],i); //Copy tracks to merger
 305     }
 306 }
 307
 308 void AliL3Hough::ProcessPatchIter(Int_t patch)
 309 {
 310   //Process patch in a iterative way.
 311   //transform + peakfinding + evaluation + transform +...
 312
 313   Int_t num_of_tries = 10;
 314   AliL3HoughBaseTransformer *tr = fHoughTransformer[patch];
 315   AliL3TrackArray *tracks = fTracks[patch];
 316   tracks->Reset();
 317   AliL3HoughEval *ev = fEval[patch];
 318   ev->InitTransformer(tr);
 319   ev->RemoveFoundTracks();
 320   ev->SetNumOfRowsToMiss(2);
 321   ev->SetNumOfPadsToLook(2);
 322   AliL3Histogram *hist;
 323   for(Int_t t=0; t<num_of_tries; t++)
 324     {
 325       tr->Reset();
 326       tr->TransformCircle();
 327       for(Int_t i=0; i<fNEtaSegments; i++)
 328         {
 329           hist = tr->GetHistogram(i);
 330           if(hist->GetNEntries()==0) continue;
 331           fPeakFinder->SetHistogram(hist);
 332           //Int_t n=1;
 333           Float_t x,y;
 334           //fPeakFinder->FindAbsMaxima(*x,*y);
 335           fPeakFinder->FindPeak(3,0.95,5,x,y);
 336           AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->NextTrack();
 337           track->SetTrackParameters(x,y,1);
 338           if(!ev->LookInsideRoad(track,i))
 339             {
 340               tracks->Remove(tracks->GetNTracks()-1);
 341               tracks->Compress();
 342             }
 343         }
 344     }
 345   LOG(AliL3Log::kInformational,"AliL3Hough::ProcessPatch","NTracks")
 346     <<AliL3Log::kDec<<"Found "<<tracks->GetNTracks()<<" tracks in patch "<<patch<<ENDLOG;
 347 }
 348
 349 void AliL3Hough::AddAllHistograms()
 350 {
 351   //Add the histograms within one etaslice.
 352   //Resulting histogram are in patch=0.
 353
 354   Double_t initTime,cpuTime;
 355   initTime = GetCpuTime();
 356   for(Int_t i=0; i<fNEtaSegments; i++)
 357     {
 358       AliL3Histogram *hist0 = fHoughTransformer[0]->GetHistogram(i);
 359       for(Int_t j=1; j<fNPatches; j++)
 360         {
 361           AliL3Histogram *hist = fHoughTransformer[j]->GetHistogram(i);
 362           hist0->Add(hist);
 363         }
 364     }
 365   fAddHistograms = kTRUE;
 366   cpuTime = GetCpuTime() - initTime;
 367   LOG(AliL3Log::kInformational,"AliL3Hough::AddAllHistograms()","Timing")
 368     <<"Adding histograms in "<<cpuTime*1000<<" ms"<<ENDLOG;
 369 }
 370
 371 void AliL3Hough::FindTrackCandidates()
 372 {
 373   //Look for peaks in histograms, and find the track candidates
 374
 375   Int_t n_patches;
 376   if(fAddHistograms)
 377     n_patches = 1; //Histograms have been added.
 378   else
 379     n_patches = fNPatches;
 380
 381   Double_t initTime,cpuTime;
 382   initTime = GetCpuTime();
 383
 384   for(Int_t i=0; i<n_patches; i++)
 385     {
 386       AliL3HoughBaseTransformer *tr = fHoughTransformer[i];
 387       fTracks[i]->Reset();
 388
 389       for(Int_t j=0; j<fNEtaSegments; j++)
 390         {
 391           AliL3Histogram *hist = tr->GetHistogram(j);
 392           if(hist->GetNEntries()==0) continue;
 393           fPeakFinder->Reset();
 394           fPeakFinder->SetHistogram(hist);
 395           //fPeakFinder->FindPeak1(3,1);
 396           fPeakFinder->FindMaxima(0,0); //Simple maxima finder
 397           //fPeakFinder->FindAbsMaxima();
 398
 399           for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
 400             {
 401               if(fPeakFinder->GetWeight(k) == 0) continue;
 402               AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[i]->NextTrack();
 403               track->SetTrackParameters(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetWeight(k));
 404               track->SetEtaIndex(j);
 405               track->SetEta(tr->GetEta(j,fCurrentSlice));
 406               track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
 407             }
 408         }
 409       fTracks[i]->QSort();
 410     }
 411   cpuTime = GetCpuTime() - initTime;
 412   LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","Timing")
 413     <<"Maxima finding done in "<<cpuTime*1000<<" ms"<<ENDLOG;
 414 }
 415
 416 void AliL3Hough::InitEvaluate()
 417 {
 418   //Pass the transformer objects to the AliL3HoughEval objects:
 419   //This will provide the evaluation objects with all the necessary
 420   //data and parameters it needs.
 421
 422   for(Int_t i=0; i<fNPatches; i++)
 423     fEval[i]->InitTransformer(fHoughTransformer[i]);
 424 }
 425
 426 Int_t AliL3Hough::Evaluate(Int_t road_width,Int_t nrowstomiss)
 427 {
 428   //Evaluate the tracks, by looking along the road in the raw data.
 429   //If track does not cross all padrows - rows2miss, it is removed from the arrray.
 430   //If histograms were not added, the check is done locally in patch,
 431   //meaning that nrowstomiss is the number of padrows the road can miss with respect
 432   //to the number of rows in the patch.
 433   //If the histograms were added, the comparison is done globally in the _slice_,
 434   //meaing that nrowstomiss is the number of padrows the road can miss with
 435   //respect to the total number of padrows in the slice.
 436   //
 437   //Return value = number of tracks which were removed (only in case of fAddHistograms)
 438
 439   if(!fTracks[0])
 440     {
 441       LOG(AliL3Log::kError,"AliL3Hough::Evaluate","Track Array")
 442         <<"No tracks to work with..."<<ENDLOG;
 443       return 0;
 444     }
 445
 446   Int_t removed_tracks=0;
 447   AliL3TrackArray *tracks=0;
 448
 449   if(fAddHistograms)
 450     {
 451       tracks = fTracks[0];
 452       for(Int_t i=0; i<tracks->GetNTracks(); i++)
 453         {
 454           AliL3Track *track = tracks->GetCheckedTrack(i);
 455           if(!track) continue;
 456           track->SetNHits(0);
 457         }
 458     }
 459
 460   for(Int_t i=0; i<fNPatches; i++)
 461     EvaluatePatch(i,road_width,nrowstomiss);
 462
 463   //Here we check the tracks globally;
 464   //how many good rows (padrows with signal)
 465   //did it cross in the slice
 466   if(fAddHistograms)
 467
 468     {
 469       for(Int_t j=0; j<tracks->GetNTracks(); j++)
 470         {
 471           AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
 472
 473           if(track->GetNHits() < AliL3Transform::GetNRows() - nrowstomiss)
 474             {
 475               tracks->Remove(j);
 476               removed_tracks++;
 477             }
 478         }
 479       tracks->Compress();
 480       tracks->QSort();
 481     }
 482
 483   return removed_tracks;
 484 }
 485
 486 void AliL3Hough::EvaluatePatch(Int_t i,Int_t road_width,Int_t nrowstomiss)
 487 {
 488   //Evaluate patch i.
 489
 490   fEval[i]->InitTransformer(fHoughTransformer[i]);
 491   fEval[i]->SetNumOfPadsToLook(road_width);
 492   fEval[i]->SetNumOfRowsToMiss(nrowstomiss);
 493   //fEval[i]->RemoveFoundTracks();
 494
 495   AliL3TrackArray *tracks=0;
 496
 497   if(!fAddHistograms)
 498     tracks = fTracks[i];
 499   else
 500     tracks = fTracks[0];
 501
 502   Int_t nrows=0;
 503   for(Int_t j=0; j<tracks->GetNTracks(); j++)
 504     {
 505       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
 506       if(!track)
 507         {
 508           LOG(AliL3Log::kWarning,"AliL3Hough::EvaluatePatch","Track array")
 509             <<"Track object missing!"<<ENDLOG;
 510           continue;
 511         }
 512       nrows=0;
 513       Bool_t result = fEval[i]->LookInsideRoad(track,nrows);
 514       if(fAddHistograms)
 515         {
 516           Int_t pre=track->GetNHits();
 517           track->SetNHits(pre+nrows);
 518         }
 519       //else//the track crossed too few good padrows (padrows with signal) in the patch, so remove it
 520       //{
 521       if(result == kFALSE)
 522         tracks->Remove(j);
 523       //}
 524     }
 525
 526   tracks->Compress();
 527   /*
 528     if(!fAddHistograms)
 529     {
 530     tracks->Compress();
 531     tracks->QSort();
 532     fMerger->FillTracks(tracks,i); //Copy tracks to the track merger
 533     }
 534   */
 535
 536 }
 537
 538 void AliL3Hough::EvaluateWithEta()
 539 {
 540   if(!fTracks[0])
 541     {
 542       printf("AliL3Hough::EvaluateWithEta: NO TRACKS\n");
 543       return;
 544     }
 545   printf("Number of tracks before evaluation %d\n",fTracks[0]->GetNTracks());
 546
 547   for(Int_t i=0; i<fNPatches; i++)
 548     {
 549       fEval[i]->InitTransformer(fHoughTransformer[i]);
 550       fEval[i]->FindEta(fTracks[0]);
 551     }
 552   fMerger->FillTracks(fTracks[0],0);
 553 }
 554
 555 void AliL3Hough::WriteTracks(Int_t slice,Char_t *path)
 556 {
 557   //Write the tracks in slice
 558
 559   AliL3MemHandler *mem = new AliL3MemHandler();
 560   Char_t fname[100];
 561   if(fAddHistograms)
 562     {
 563       sprintf(fname,"%s/tracks_ho_%d.raw",path,slice);
 564       mem->SetBinaryOutput(fname);
 565       mem->TrackArray2Binary(fTracks[0]);
 566       mem->CloseBinaryOutput();
 567     }
 568   else
 569     {
 570       for(Int_t i=0; i<fNPatches; i++)
 571         {
 572           sprintf(fname,"%s/tracks_ho_%d_%d.raw",path,slice,i);
 573           mem->SetBinaryOutput(fname);
 574           mem->TrackArray2Binary(fTracks[i]);
 575           mem->CloseBinaryOutput();
 576         }
 577     }
 578   delete mem;
 579
 580 }
 581
 582 void AliL3Hough::WriteDigits(Char_t *outfile)
 583 {
 584 #ifdef use_aliroot
 585   //Write the current data to a new rootfile.
 586
 587   for(Int_t i=0; i<fNPatches; i++)
 588     {
 589       AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer[i]->GetDataPointer();
 590       fMemHandler[i]->AliDigits2RootFile(tempPt,outfile);
 591     }
 592 #else
 593   cerr<<"AliL3Hough::WriteDigits : You need to compile with AliROOT!"<<endl;
 594   return;
 595 #endif
 596 }
 597
 598 Double_t AliL3Hough::GetCpuTime()
 599 {
 600   //Return the Cputime in seconds.
 601  struct timeval tv;
 602  gettimeofday( &tv, NULL );
 603  return tv.tv_sec+(((Double_t)tv.tv_usec)/1000000.);
 604  //return (Double_t)(clock()) / CLOCKS_PER_SEC;
 605 }
 606