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