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