0e8b6510df2afad188026eb7da4d964fb827b1cc
[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
9 #include "AliL3HoughMerger.h"
10 #include "AliL3HoughIntMerger.h"
11 #include "AliL3HoughGlobalMerger.h"
12 #include "AliL3Logging.h"
13 #include "AliL3Histogram.h"
14 #include "AliL3Hough.h"
15 #include "AliL3HoughTransformer.h"
16 #include "AliL3HoughMaxFinder.h"
17 #ifdef use_aliroot
18 #include "AliL3FileHandler.h"
19 #else
20 #include "AliL3MemHandler.h"
21 #endif
22 #include "AliL3DigitData.h"
23 #include "AliL3HoughEval.h"
24 #include "AliL3Transform.h"
25 #include "AliL3Defs.h"
26 #include "AliL3TrackArray.h"
27 #include "AliL3HoughTrack.h"
28
29
30 //_____________________________________________________________
31 // AliL3Hough
32 //
33 // Interface class for the Hough transform
34 //
35 // Example how to use:
36 //
37 // AliL3Hough *hough = new AliL3Hough(path,kTRUE,NumberOfEtaSegments);
38 // hough->ReadData(slice);
39 // hough->Transform();
40 // hough->FindTrackCandidates();
41 // 
42 // AliL3TrackArray *tracks = hough->GetTracks(patch);
43
44 ClassImp(AliL3Hough)
45
46 AliL3Hough::AliL3Hough()
47 {
48   //Constructor
49   
50   fBinary = kFALSE;
51   fNEtaSegments = 0;
52   fAddHistograms = kFALSE;
53   fDoIterative = kFALSE; 
54   fWriteDigits=kFALSE;
55   fNPatches=0;
56   fMemHandler = 0;
57   fHoughTransformer = 0;
58   fEval = 0;
59   fPeakFinder = 0;
60   fTracks = 0;
61   fMerger = 0;
62   fInterMerger = 0;
63   fGlobalMerger = 0;
64   fTransform = 0;
65 }
66
67
68 AliL3Hough::AliL3Hough(Char_t *path,Bool_t binary,Int_t n_eta_segments)
69 {
70   //Default ctor.
71
72   fBinary = binary;
73   strcpy(fPath,path);
74   fNEtaSegments = n_eta_segments;
75   fAddHistograms = kFALSE;
76   fDoIterative = kFALSE; 
77   fWriteDigits = kFALSE;
78   Init();
79 }
80
81
82 AliL3Hough::~AliL3Hough()
83 {
84   //dtor
85
86   CleanUp();
87   if(fMerger)
88     delete fMerger;
89   if(fInterMerger)
90     delete fInterMerger;
91   if(fPeakFinder)
92     delete fPeakFinder;
93   if(fGlobalMerger)
94     delete fGlobalMerger;
95   if(fTransform)
96     delete fTransform;
97 }
98
99 void AliL3Hough::CleanUp()
100 {
101   //Cleanup memory
102   
103   for(Int_t i=0; i<fNPatches; i++)
104     {
105       if(fTracks[i]) delete fTracks[i];
106       if(fEval[i]) delete fEval[i];
107       if(fHoughTransformer[i]) delete fHoughTransformer[i];
108       if(fMemHandler[i]) delete fMemHandler[i];
109     }
110   
111   /*    
112         if(fTracks) delete [] fTracks;
113         if(fEval) delete [] fEval;
114         if(fHoughTransformer) delete [] fHoughTransformer;
115         if(fMemHandler) delete [] fMemHandler;
116   */
117 }
118
119 void AliL3Hough::Init()
120 {
121   fPeakThreshold = 0;
122   fNPatches = NPatches;
123   fTransform = new AliL3Transform(fPath);
124   fHoughTransformer = new AliL3HoughBaseTransformer*[fNPatches];
125 #ifdef use_aliroot
126   fMemHandler = new AliL3FileHandler*[fNPatches];
127 #else
128   fMemHandler = new AliL3MemHandler*[fNPatches];
129 #endif
130   fTracks = new AliL3TrackArray*[fNPatches];
131   fEval = new AliL3HoughEval*[fNPatches];
132   for(Int_t i=0; i<fNPatches; i++)
133     {
134       fHoughTransformer[i] = new AliL3HoughTransformer(1,i,fNEtaSegments);
135       fHoughTransformer[i]->SetTransformer(fTransform);
136       //fHoughTransformer[i]->CreateHistograms(64,-0.003,0.003,64,-0.26,0.26);
137       fHoughTransformer[i]->CreateHistograms(64,0.1,64,-30,30);
138       fHoughTransformer[i]->SetThreshold(3);
139       fEval[i] = new AliL3HoughEval();
140       fTracks[i] = new AliL3TrackArray("AliL3HoughTrack");
141 #ifdef use_aliroot
142       fMemHandler[i] = new AliL3FileHandler();
143       if(!fBinary)
144         fMemHandler[i]->SetAliInput(fPath);
145 #else
146       fMemHandler[i] = new AliL3MemHandler();
147 #endif
148       
149     }
150   fPeakFinder = new AliL3HoughMaxFinder("KappaPhi");
151   fMerger = new AliL3HoughMerger(fNPatches);
152   fInterMerger = new AliL3HoughIntMerger();
153   fGlobalMerger = 0;
154 }
155
156 void AliL3Hough::Process(Int_t minslice,Int_t maxslice)
157 {
158   //Process all slices [minslice,maxslice].
159   fGlobalMerger = new AliL3HoughGlobalMerger(minslice,maxslice);
160   
161   for(Int_t i=minslice; i<=maxslice; i++)
162     {
163       ReadData(i);
164       Transform();
165       if(fAddHistograms)
166         AddAllHistograms();
167       FindTrackCandidates();
168       Evaluate();
169       fGlobalMerger->FillTracks(fTracks[0],i);
170     }
171   
172   
173 }
174
175 void AliL3Hough::ReadData(Int_t slice)
176 {
177   //Read data from files, binary or root.
178
179   for(Int_t i=0; i<fNPatches; i++)
180     {
181       fMemHandler[i]->Free();
182       UInt_t ndigits=0;
183       AliL3DigitRowData *digits =0;
184       Char_t name[256];
185       if(fBinary)//take input data from binary files
186         {
187           sprintf(name,"%sdigits_%d_%d.raw",fPath,slice,i);
188           fMemHandler[i]->SetBinaryInput(name);
189           digits = (AliL3DigitRowData *)fMemHandler[i]->CompBinary2Memory(ndigits);
190           fMemHandler[i]->CloseBinaryInput();
191         }
192       else //read data from root file
193         {
194 #ifdef use_aliroot
195           fMemHandler[i]->Init(slice,i,NRows[i]);
196           digits=(AliL3DigitRowData *)fMemHandler[i]->AliDigits2Memory(ndigits); 
197 #else
198           cerr<<"You cannot read from rootfile now"<<endl;
199 #endif
200         }
201       fHoughTransformer[i]->SetInputData(ndigits,digits);
202     }
203 }
204
205 void AliL3Hough::Transform(Int_t row_range)
206 {
207   //Transform all data given to the transformer within the given slice
208   //(after ReadData(slice))
209
210   for(Int_t i=0; i<fNPatches; i++)
211     {
212       fHoughTransformer[i]->Reset();//Reset the histograms
213       if(row_range < 0)
214         fHoughTransformer[i]->TransformCircle();
215       else
216         fHoughTransformer[i]->TransformCircleC(row_range);
217     }
218 }
219
220 void AliL3Hough::MergePatches()
221 {
222   if(fAddHistograms) //Nothing to merge here
223     return;
224   fMerger->SetTransformer(fTransform);
225   fMerger->MergePatches(kTRUE);
226 }
227
228 void AliL3Hough::MergeInternally()
229 {
230   if(fAddHistograms)
231     fInterMerger->FillTracks(fTracks[0]);
232   else
233     fInterMerger->FillTracks(fMerger->GetOutTracks());
234   
235   fInterMerger->MMerge();
236 }
237
238 void AliL3Hough::ProcessSliceIter()
239 {
240   //Process current slice (after ReadData(slice)) iteratively.
241   
242   for(Int_t i=0; i<fNPatches; i++)
243     {
244       ProcessPatchIter(i);
245       fMerger->FillTracks(fTracks[i],i); //Copy tracks to merger
246     }
247   
248 }
249
250 void AliL3Hough::ProcessPatchIter(Int_t patch)
251 {
252   //Process patch in a iterative way. 
253   //transform + peakfinding + evaluation + transform +...
254
255   Int_t num_of_tries = 10;
256   AliL3HoughBaseTransformer *tr = fHoughTransformer[patch];
257   AliL3TrackArray *tracks = fTracks[patch];
258   tracks->Reset();
259   AliL3HoughEval *ev = fEval[patch];
260   ev->InitTransformer(tr);
261   ev->RemoveFoundTracks();
262   ev->SetNumOfRowsToMiss(2);
263   ev->SetNumOfPadsToLook(2);
264   AliL3Histogram *hist;
265   for(Int_t t=0; t<num_of_tries; t++)
266     {
267       tr->Reset();
268       tr->TransformCircle();
269       for(Int_t i=0; i<fNEtaSegments; i++)
270         {
271           hist = tr->GetHistogram(i);
272           if(hist->GetNEntries()==0) continue;
273           fPeakFinder->SetHistogram(hist);
274           //Int_t n=1;
275           Float_t x,y;
276           //fPeakFinder->FindAbsMaxima(*x,*y);
277           fPeakFinder->FindPeak(3,0.95,5,x,y);
278           AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->NextTrack();
279           track->SetTrackParameters(x,y,1);
280           if(!ev->LookInsideRoad(track,i))
281             {   
282               tracks->Remove(tracks->GetNTracks()-1);
283               tracks->Compress();
284             }
285         }
286     }
287   LOG(AliL3Log::kInformational,"AliL3Hough::ProcessPatch","NTracks")
288     <<AliL3Log::kDec<<"Found "<<tracks->GetNTracks()<<" tracks in patch "<<patch<<ENDLOG;
289 }
290
291
292 void AliL3Hough::AddAllHistograms()
293 {
294   //Add the histograms within one etaslice.
295   //Resulting histogram are in patch=0.
296
297   for(Int_t i=0; i<fNEtaSegments; i++)
298     {
299       AliL3Histogram *hist0 = fHoughTransformer[0]->GetHistogram(i);
300       for(Int_t j=1; j<fNPatches; j++)
301         {
302           AliL3Histogram *hist = fHoughTransformer[j]->GetHistogram(i);
303           hist0->Add(hist);
304         }
305     }
306   fAddHistograms = kTRUE;
307 }
308
309 void AliL3Hough::FindTrackCandidates()
310 {
311   //Look for peaks in histograms, and find the track candidates
312   
313   Int_t n_patches;
314   if(fAddHistograms)
315     n_patches = 1; //Histograms has been added.
316   else
317     n_patches = fNPatches;
318
319   
320   for(Int_t i=0; i<n_patches; i++)
321     {
322       AliL3HoughBaseTransformer *tr = fHoughTransformer[i];
323       Double_t eta_slice = (tr->GetEtaMax()-tr->GetEtaMin()/tr->GetNEtaSegments());
324       fTracks[i]->Reset();
325       for(Int_t j=0; j<fNEtaSegments; j++)
326         {
327           AliL3Histogram *hist = tr->GetHistogram(j);
328           if(hist->GetNEntries()==0) continue;
329           fPeakFinder->SetHistogram(hist);
330           fPeakFinder->SetThreshold(fPeakThreshold);
331           Int_t n=20;
332           Float_t x[n];
333           Float_t y[n];
334           Int_t weight[n];
335           //fPeakFinder->FindPeak1(x,y,weight,n,2,1);
336           fPeakFinder->FindMaxima(x,y,weight,n);
337           for(Int_t k=0; k<n; k++)
338             {
339               if(weight[k] == 0) continue;
340               
341               AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[i]->NextTrack();
342               track->SetTrackParameters(x[k],y[k],weight[k]);
343               track->SetEtaIndex(j);
344               track->SetEta((Double_t)(j*eta_slice));
345               track->SetRowRange(NRows[0][0],NRows[5][1]);
346             }
347         }
348       fTracks[i]->QSort();
349     }
350 }
351
352 void AliL3Hough::Evaluate(Int_t road_width)
353 {
354   //Evaluate the tracks, by looking along the road in the raw data.
355   
356
357   if(!fTracks[0])
358     {
359       LOG(AliL3Log::kError,"AliL3Hough::Evaluate","Track Array")
360         <<"No tracks to work with..."<<ENDLOG;
361       return;
362     }
363   
364   printf("Number of tracks before evaluation %d\n",fTracks[0]->GetNTracks());
365   AliL3TrackArray *tracks;
366   for(Int_t i=0; i<fNPatches; i++)
367     {
368       fEval[i]->InitTransformer(fHoughTransformer[i]);
369       continue;
370       fEval[i]->SetNumOfRowsToMiss(2);
371       fEval[i]->SetNumOfPadsToLook(road_width);
372       if(fAddHistograms)
373         tracks = fTracks[0];
374       else
375         tracks = fTracks[i];
376       for(Int_t j=0; j<tracks->GetNTracks(); j++)
377         {
378           AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
379           if(!track)
380             {
381               LOG(AliL3Log::kWarning,"AliL3Hough::Evaluate","Track array")
382                 <<"Track object missing!"<<ENDLOG;
383               continue;
384             }
385           
386           if(!fEval[i]->LookInsideRoad(track,track->GetEtaIndex()))
387             tracks->Remove(j);
388           if(fAddHistograms)
389             track->SetRowRange(NRows[0][0],NRows[5][1]);//All rows included
390         }
391       tracks->Compress();
392       tracks->QSort(); //Sort the tracks according to weight
393       
394       if(!fAddHistograms)
395         fMerger->FillTracks(tracks,i); //Copy tracks to the track merger
396     }
397   
398 }
399
400 void AliL3Hough::EvaluateWithEta()
401 {
402   if(!fTracks[0])
403     {
404       printf("AliL3Hough::EvaluateWithEta: NO TRACKS\n");
405       return;
406     }
407   printf("Number of tracks before evaluation %d\n",fTracks[0]->GetNTracks());
408  
409   for(Int_t i=0; i<fNPatches; i++)
410     {
411       fEval[i]->InitTransformer(fHoughTransformer[i]);
412       fEval[i]->FindEta(fTracks[0]);
413     }
414   fMerger->FillTracks(fTracks[0],0);
415 }
416
417 void AliL3Hough::WriteTracks(Char_t *path)
418 {
419   AliL3MemHandler *mem = new AliL3MemHandler();
420   Char_t fname[100];
421   if(fAddHistograms)
422     {
423       sprintf(fname,"%s/tracks.raw",path);
424       mem->SetBinaryOutput(fname);
425       mem->TrackArray2Binary(fTracks[0]);
426       mem->CloseBinaryOutput();
427     }
428   else 
429     {
430       for(Int_t i=0; i<fNPatches; i++)
431         {
432           sprintf(fname,"%s/tracks_%d.raw",path,i);
433           mem->SetBinaryOutput(fname);
434           mem->TrackArray2Binary(fTracks[i]);
435           mem->CloseBinaryOutput();
436         }
437     }
438   delete mem;
439   
440 }
441 #ifdef use_aliroot
442 void AliL3Hough::WriteDigits(Char_t *outfile)
443 {
444   //Write the current data to a new rootfile.
445
446   for(Int_t i=0; i<fNPatches; i++)
447     {
448       AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer[i]->GetDataPointer();
449       fMemHandler[i]->AliDigits2RootFile(tempPt,outfile);
450     }
451   
452 }
453 #endif