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