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