Added support for NEWIO, merged cern-hlt tree, updated to latest track candidate...
[u/mrichter/AliRoot.git] / HLT / hough / AliL3Hough.cxx
1 // @(#) $Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ALICE HLT Group
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 #include "AliL3Benchmark.h"
21 #ifdef use_aliroot
22 #include "AliL3FileHandler.h"
23 #else
24 #include "AliL3MemHandler.h"
25 #endif
26 #include "AliL3DataHandler.h"
27 #include "AliL3DigitData.h"
28 #include "AliL3HoughEval.h"
29 #include "AliL3Transform.h"
30 #include "AliL3TrackArray.h"
31 #include "AliL3HoughTrack.h"
32 #include "AliL3DDLDataFileHandler.h"
33
34 #if GCCVERSION == 3
35 using namespace std;
36 #endif
37
38 /** /class AliL3Hough
39 //<pre>
40 //_____________________________________________________________
41 // AliL3Hough
42 //
43 // Interface class for the Hough transform
44 //
45 // Example how to use:
46 //
47 // AliL3Hough *hough = new AliL3Hough(path,kTRUE,NumberOfEtaSegments);
48 // hough->ReadData(slice);
49 // hough->Transform();
50 // hough->FindTrackCandidates();
51 // 
52 // AliL3TrackArray *tracks = hough->GetTracks(patch);
53 //</pre>
54 */
55
56 ClassImp(AliL3Hough)
57
58 AliL3Hough::AliL3Hough()
59 {
60   //Constructor
61   
62   fBinary        = kFALSE;
63   fAddHistograms = kFALSE;
64   fDoIterative   = kFALSE; 
65   fWriteDigits   = kFALSE;
66   fUse8bits      = kFALSE;
67
68   fMemHandler       = 0;
69   fHoughTransformer = 0;
70   fEval             = 0;
71   fPeakFinder       = 0;
72   fTracks           = 0;
73   fGlobalTracks     = 0;
74   fMerger           = 0;
75   fInterMerger      = 0;
76   fGlobalMerger     = 0;
77   fBenchmark        = 0;
78   
79   fNEtaSegments     = 0;
80   fNPatches         = 0;
81   fVersion          = 0;
82   fCurrentSlice     = 0;
83   fEvent            = 0;
84
85   SetTransformerParams();
86   SetThreshold();
87   SetNSaveIterations();
88   SetPeakThreshold();
89 }
90
91 AliL3Hough::AliL3Hough(Char_t *path,Bool_t binary,Int_t n_eta_segments,Bool_t bit8,Int_t tv,Char_t *infile)
92 {
93   //Default ctor.
94
95   fBinary = binary;
96   strcpy(fPath,path);
97   fNEtaSegments  = n_eta_segments;
98   fAddHistograms = kFALSE;
99   fDoIterative   = kFALSE; 
100   fWriteDigits   = kFALSE;
101   fUse8bits      = bit8;
102   fVersion       = tv;
103   if(!fBinary)
104     fInputFile = infile;
105   else
106     fInputFile = 0;
107 }
108
109 AliL3Hough::~AliL3Hough()
110 {
111   //dtor
112
113   CleanUp();
114   if(fMerger)
115     delete fMerger;
116   if(fInterMerger)
117     delete fInterMerger;
118   if(fPeakFinder)
119     delete fPeakFinder;
120   if(fGlobalMerger)
121     delete fGlobalMerger;
122   if(fBenchmark)
123     delete fBenchmark;
124   if(fGlobalTracks)
125     delete fGlobalTracks;
126 }
127
128 void AliL3Hough::CleanUp()
129 {
130   //Cleanup memory
131   
132   for(Int_t i=0; i<fNPatches; i++)
133     {
134       if(fTracks[i]) delete fTracks[i];
135       if(fEval[i]) delete fEval[i];
136       if(fHoughTransformer[i]) delete fHoughTransformer[i];
137       if(fMemHandler[i]) delete fMemHandler[i];
138     }
139   
140   if(fTracks) delete [] fTracks;
141   if(fEval) delete [] fEval;
142   if(fHoughTransformer) delete [] fHoughTransformer;
143   if(fMemHandler) delete [] fMemHandler;
144   
145 }
146
147 void AliL3Hough::Init(Char_t *path,Bool_t binary,Int_t n_eta_segments,Bool_t bit8,Int_t tv,Char_t *infile)
148 {
149   fBinary = binary;
150   strcpy(fPath,path);
151   fNEtaSegments = n_eta_segments;
152   fWriteDigits  = kFALSE;
153   fUse8bits     = bit8;
154   fVersion      = tv;
155   if(!fBinary)
156     fInputFile = infile;
157   else
158     fInputFile = 0;
159
160   Init(); //do the rest
161 }
162
163 void AliL3Hough::Init(Bool_t doit, Bool_t addhists)
164 {
165   fDoIterative   = doit; 
166   fAddHistograms = addhists;
167
168   fNPatches = AliL3Transform::GetNPatches();
169   
170   fHoughTransformer = new AliL3HoughBaseTransformer*[fNPatches];
171   fMemHandler = new AliL3MemHandler*[fNPatches];
172
173   fTracks = new AliL3TrackArray*[fNPatches];
174   fEval = new AliL3HoughEval*[fNPatches];
175   
176   fGlobalTracks = new AliL3TrackArray("AliL3HoughTrack");
177   
178   for(Int_t i=0; i<fNPatches; i++)
179     {
180       switch (fVersion){ //choose Transformer
181       case 1: 
182         fHoughTransformer[i] = new AliL3HoughTransformerLUT(0,i,fNEtaSegments);
183         break;
184       case 2:
185         fHoughTransformer[i] = new AliL3HoughClusterTransformer(0,i,fNEtaSegments);
186         break;
187       case 3:
188         fHoughTransformer[i] = new AliL3HoughTransformerVhdl(0,i,fNEtaSegments,fNSaveIterations);
189         break;
190       default:
191         fHoughTransformer[i] = new AliL3HoughTransformer(0,i,fNEtaSegments,kFALSE,kFALSE);
192       }
193
194       fHoughTransformer[i]->CreateHistograms(fNBinX[i],fLowPt[i],fNBinY[i],-fPhi[i],fPhi[i]);
195       //fHoughTransformer[i]->CreateHistograms(fLowPt[i],fUpperPt[i],fPtRes[i],fNBinY[i],fPhi[i]);
196       fHoughTransformer[i]->SetLowerThreshold(fThreshold[i]);
197       fHoughTransformer[i]->SetUpperThreshold(100);
198
199       LOG(AliL3Log::kInformational,"AliL3Hough::Init","Version")
200         <<"Initializing Hough transformer version "<<fVersion<<ENDLOG;
201       
202       fEval[i] = new AliL3HoughEval();
203       fTracks[i] = new AliL3TrackArray("AliL3HoughTrack");
204       if(fUse8bits)
205         fMemHandler[i] = new AliL3DataHandler();
206       else
207 #ifdef use_aliroot
208         {
209           if(!fInputFile) {
210             /* In case of reading digits file */
211             fMemHandler[i] = new AliL3FileHandler();
212             if(!fBinary) {
213               Char_t filename[1024];
214               sprintf(filename,"%s/digitfile.root",fPath);
215               fMemHandler[i]->SetAliInput(filename);
216             }
217           }
218           else {
219             /* In case of reading rawdata from ROOT file */
220             fMemHandler[i] = new AliL3DDLDataFileHandler();
221             fMemHandler[i]->SetReaderInput(fInputFile);
222           }
223         }
224 #else
225       fMemHandler[i] = new AliL3MemHandler();
226 #endif
227     }
228
229   fPeakFinder = new AliL3HoughMaxFinder("KappaPhi",1000);
230   fMerger = new AliL3HoughMerger(fNPatches);
231   fInterMerger = new AliL3HoughIntMerger();
232   fGlobalMerger = 0;
233   fBenchmark = new AliL3Benchmark();
234 }
235
236 void AliL3Hough::SetTransformerParams(Float_t ptres,Float_t ptmin,Float_t ptmax,Int_t ny,Int_t patch)
237 {
238
239   Int_t mrow;
240   Float_t psi=0;
241   if(patch==-1)
242     mrow = 80;
243   else
244     mrow = AliL3Transform::GetLastRow(patch);
245   if(ptmin)
246     {
247       Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(mrow),2) + pow(AliL3Transform::GetMaxY(mrow),2));
248       Double_t kappa = -1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/ptmin;
249       psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
250       cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
251     }
252
253   if(patch==-1)
254     {
255       Int_t i=0;
256       while(i < 6)
257         {
258           fPtRes[i] = ptres;
259           fLowPt[i] = ptmin;
260           fUpperPt[i] = ptmax;
261           fNBinY[i] = ny;
262           fPhi[i] = psi;
263           fNBinX[i]=0;
264           i++;
265         }
266       return;
267     }
268
269   fPtRes[patch] = ptres;
270   fLowPt[patch] = ptmin;
271   fUpperPt[patch] = ptmax;
272   fNBinY[patch] = ny;
273   fPhi[patch] = psi;
274 }
275
276 void AliL3Hough::SetTransformerParams(Int_t nx,Int_t ny,Float_t ptmin,Int_t patch)
277 {
278
279   Int_t mrow=80;
280   Double_t lineradius = sqrt(pow(AliL3Transform::Row2X(mrow),2) + pow(AliL3Transform::GetMaxY(mrow),2));
281   Double_t kappa = -1*AliL3Transform::GetBField()*AliL3Transform::GetBFact()/ptmin;
282   Double_t psi = AliL3Transform::Deg2Rad(10) - asin(lineradius*kappa/2);
283   cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
284   
285   Int_t i=0;
286   while(i < 6)
287     {
288       fLowPt[i] = ptmin;
289       fNBinY[i] = ny;
290       fNBinX[i] = nx;
291       fPhi[i] = psi;
292       i++;
293     }
294 }
295
296 void AliL3Hough::SetTransformerParams(Int_t nx,Int_t ny,Float_t lpt,Float_t phi)
297 {
298   Int_t i=0;
299   while(i < 6)
300     {
301       fLowPt[i] = lpt;
302       fNBinY[i] = ny;
303       fNBinX[i] = nx;
304       fPhi[i] = phi;
305       i++;
306     }
307 }
308
309 void AliL3Hough::SetThreshold(Int_t t3,Int_t patch)
310 {
311   if(patch==-1)
312     {
313       Int_t i=0;
314       while(i < 6)
315         fThreshold[i++]=t3;
316       return;
317     }
318   fThreshold[patch]=t3;
319 }
320
321 void AliL3Hough::SetPeakThreshold(Int_t threshold,Int_t patch)
322 {
323   if(patch==-1)
324     {
325       Int_t i=0;
326       while(i < 6)
327         fPeakThreshold[i++]=threshold;
328       return;
329     }
330   fPeakThreshold[patch]=threshold;
331 }
332
333 void AliL3Hough::DoBench(Char_t *name)
334 {
335   fBenchmark->Analyze(name);
336 }
337
338 void AliL3Hough::Process(Int_t minslice,Int_t maxslice)
339 {
340   //Process all slices [minslice,maxslice].
341   fGlobalMerger = new AliL3HoughGlobalMerger(minslice,maxslice);
342   
343   for(Int_t i=minslice; i<=maxslice; i++)
344     {
345       ReadData(i);
346       Transform();
347       if(fAddHistograms)
348         AddAllHistograms();
349       FindTrackCandidates();
350       //Evaluate();
351       //fGlobalMerger->FillTracks(fTracks[0],i);
352     }
353 }
354
355 void AliL3Hough::ReadData(Int_t slice,Int_t eventnr)
356 {
357   //Read data from files, binary or root.
358   
359   fEvent=eventnr;
360   fCurrentSlice = slice;
361   for(Int_t i=0; i<fNPatches; i++)
362     {
363       fMemHandler[i]->Free();
364       UInt_t ndigits=0;
365       AliL3DigitRowData *digits =0;
366       Char_t name[256];
367       fMemHandler[i]->Init(slice,i);
368       if(fBinary)//take input data from binary files
369         {
370           if(fUse8bits)
371             sprintf(name,"%s/binaries/digits_c8_%d_%d_%d.raw",fPath,eventnr,slice,i);
372           else
373             sprintf(name,"%s/binaries/digits_%d_%d_%d.raw",fPath,eventnr,slice,i);
374
375           fMemHandler[i]->SetBinaryInput(name);
376           digits = (AliL3DigitRowData *)fMemHandler[i]->CompBinary2Memory(ndigits);
377           fMemHandler[i]->CloseBinaryInput();
378         }
379       else //read data from root file
380         {
381 #ifdef use_aliroot
382           digits=(AliL3DigitRowData *)fMemHandler[i]->AliAltroDigits2Memory(ndigits,eventnr);
383 #else
384           cerr<<"You cannot read from rootfile now"<<endl;
385 #endif
386         }
387
388       //set input data and init transformer
389       fHoughTransformer[i]->SetInputData(ndigits,digits);
390       fHoughTransformer[i]->Init(slice,i,fNEtaSegments);
391     }
392 }
393
394 void AliL3Hough::Transform(Int_t *row_range)
395 {
396   //Transform all data given to the transformer within the given slice
397   //(after ReadData(slice))
398   
399   Double_t initTime,cpuTime;
400   initTime = GetCpuTime();
401   for(Int_t i=0; i<fNPatches; i++)
402     {
403       fHoughTransformer[i]->Reset();//Reset the histograms
404       fBenchmark->Start("Hough Transform");
405       if(!row_range)
406         fHoughTransformer[i]->TransformCircle();
407       else
408         fHoughTransformer[i]->TransformCircleC(row_range,1);
409       fBenchmark->Stop("Hough Transform");
410     }
411   cpuTime = GetCpuTime() - initTime;
412   LOG(AliL3Log::kInformational,"AliL3Hough::Transform()","Timing")
413     <<"Transform done in average per patch of "<<cpuTime*1000/fNPatches<<" ms"<<ENDLOG;
414 }
415
416 void AliL3Hough::MergePatches()
417 {
418   if(fAddHistograms) //Nothing to merge here
419     return;
420   fMerger->MergePatches(kTRUE);
421 }
422
423 void AliL3Hough::MergeInternally()
424 {
425   if(fAddHistograms)
426     fInterMerger->FillTracks(fTracks[0]);
427   else
428     fInterMerger->FillTracks(fMerger->GetOutTracks());
429   
430   fInterMerger->MMerge();
431 }
432
433 void AliL3Hough::ProcessSliceIter()
434 {
435   //Process current slice (after ReadData(slice)) iteratively.
436   
437   if(!fAddHistograms)
438     {
439       for(Int_t i=0; i<fNPatches; i++)
440         {
441           ProcessPatchIter(i);
442           fMerger->FillTracks(fTracks[i],i); //Copy tracks to merger
443         }
444     }
445   else
446     {
447       for(Int_t i=0; i<10; i++)
448         {
449           Transform();
450           AddAllHistograms();
451           InitEvaluate();
452           AliL3HoughBaseTransformer *tr = fHoughTransformer[0];
453           for(Int_t j=0; j<fNEtaSegments; j++)
454             {
455               AliL3Histogram *hist = tr->GetHistogram(j);
456               if(hist->GetNEntries()==0) continue;
457               fPeakFinder->Reset();
458               fPeakFinder->SetHistogram(hist);
459               fPeakFinder->FindAbsMaxima();
460               AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[0]->NextTrack();
461               track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
462               track->SetEtaIndex(j);
463               track->SetEta(tr->GetEta(j,fCurrentSlice));
464               for(Int_t k=0; k<fNPatches; k++)
465                 {
466                   fEval[i]->SetNumOfPadsToLook(2);
467                   fEval[i]->SetNumOfRowsToMiss(2);
468                   fEval[i]->RemoveFoundTracks();
469                   /*
470                   Int_t nrows=0;
471                   if(!fEval[i]->LookInsideRoad(track,nrows))
472                     {
473                       fTracks[0]->Remove(fTracks[0]->GetNTracks()-1);
474                       fTracks[0]->Compress();
475                     }
476                   */
477                 }
478             }
479           
480         }
481       
482     }
483 }
484
485 void AliL3Hough::ProcessPatchIter(Int_t patch)
486 {
487   //Process patch in a iterative way. 
488   //transform + peakfinding + evaluation + transform +...
489
490   Int_t num_of_tries = 5;
491   AliL3HoughBaseTransformer *tr = fHoughTransformer[patch];
492   AliL3TrackArray *tracks = fTracks[patch];
493   tracks->Reset();
494   AliL3HoughEval *ev = fEval[patch];
495   ev->InitTransformer(tr);
496   //ev->RemoveFoundTracks();
497   ev->SetNumOfRowsToMiss(3);
498   ev->SetNumOfPadsToLook(2);
499   AliL3Histogram *hist;
500   for(Int_t t=0; t<num_of_tries; t++)
501     {
502       tr->Reset();
503       tr->TransformCircle();
504       for(Int_t i=0; i<fNEtaSegments; i++)
505         {
506           hist = tr->GetHistogram(i);
507           if(hist->GetNEntries()==0) continue;
508           fPeakFinder->Reset();
509           fPeakFinder->SetHistogram(hist);
510           fPeakFinder->FindAbsMaxima();
511           //fPeakFinder->FindPeak1();
512           AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->NextTrack();
513           track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
514           track->SetEtaIndex(i);
515           track->SetEta(tr->GetEta(i,fCurrentSlice));
516           /*
517           Int_t nrows=0;
518           if(!ev->LookInsideRoad(track,nrows))
519             {   
520               tracks->Remove(tracks->GetNTracks()-1);
521               tracks->Compress();
522             }
523           */
524         }
525     }
526   fTracks[0]->QSort();
527   LOG(AliL3Log::kInformational,"AliL3Hough::ProcessPatch","NTracks")
528     <<AliL3Log::kDec<<"Found "<<tracks->GetNTracks()<<" tracks in patch "<<patch<<ENDLOG;
529 }
530
531 void AliL3Hough::AddAllHistograms()
532 {
533   //Add the histograms within one etaslice.
534   //Resulting histogram are in patch=0.
535
536   Double_t initTime,cpuTime;
537   initTime = GetCpuTime();
538   fBenchmark->Start("Add Histograms");
539   for(Int_t i=0; i<fNEtaSegments; i++)
540     {
541       AliL3Histogram *hist0 = fHoughTransformer[0]->GetHistogram(i);
542       for(Int_t j=1; j<fNPatches; j++)
543         {
544           AliL3Histogram *hist = fHoughTransformer[j]->GetHistogram(i);
545           hist0->Add(hist);
546         }
547     }
548   fBenchmark->Stop("Add Histograms");
549   fAddHistograms = kTRUE;
550   cpuTime = GetCpuTime() - initTime;
551   LOG(AliL3Log::kInformational,"AliL3Hough::AddAllHistograms()","Timing")
552     <<"Adding histograms in "<<cpuTime*1000<<" ms"<<ENDLOG;
553 }
554
555 void AliL3Hough::AddTracks()
556 {
557   if(!fTracks[0])
558     {
559       cerr<<"AliL3Hough::AddTracks : No tracks"<<endl;
560       return;
561     }
562   AliL3TrackArray *tracks = fTracks[0];
563   for(Int_t i=0; i<tracks->GetNTracks(); i++)
564     {
565       AliL3Track *track = tracks->GetCheckedTrack(i);
566       if(!track) continue;
567       if(track->GetNHits()!=1) cerr<<"NHITS "<<track->GetNHits()<<endl;
568       UInt_t *ids = track->GetHitNumbers();
569       ids[0] = (fCurrentSlice&0x7f)<<25;
570     }
571   
572   fGlobalTracks->AddTracks(fTracks[0],0,fCurrentSlice);
573 }
574
575 void AliL3Hough::FindTrackCandidates()
576 {
577   //Look for peaks in histograms, and find the track candidates
578   
579   Int_t n_patches;
580   if(fAddHistograms)
581     n_patches = 1; //Histograms have been added.
582   else
583     n_patches = fNPatches;
584   
585   Double_t initTime,cpuTime;
586   initTime = GetCpuTime();
587   fBenchmark->Start("Find Maxima");
588   for(Int_t i=0; i<n_patches; i++)
589     {
590       AliL3HoughBaseTransformer *tr = fHoughTransformer[i];
591       fTracks[i]->Reset();
592       
593       for(Int_t j=0; j<fNEtaSegments; j++)
594         {
595           AliL3Histogram *hist = tr->GetHistogram(j);
596           if(hist->GetNEntries()==0) continue;
597           fPeakFinder->Reset();
598           fPeakFinder->SetHistogram(hist);
599           
600           //fPeakFinder->FindAbsMaxima();
601           fPeakFinder->SetThreshold(fPeakThreshold[i]);
602           fPeakFinder->FindAdaptedPeaks(6);
603           //fPeakFinder->FindMaxima(fPeakThreshold[i]); //Simple maxima finder
604           
605           for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
606             {
607               AliL3HoughTrack *track = (AliL3HoughTrack*)fTracks[i]->NextTrack();
608               track->SetTrackParameters(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetWeight(k));
609               track->SetEtaIndex(j);
610               track->SetEta(tr->GetEta(j,fCurrentSlice));
611               track->SetRowRange(AliL3Transform::GetFirstRow(0),AliL3Transform::GetLastRow(5));
612             }
613         }
614       cout<<"Found "<<fTracks[i]->GetNTracks()<<" tracks in patch "<<i<<endl;
615       fTracks[i]->QSort();
616     }
617   fBenchmark->Stop("Find Maxima");
618   cpuTime = GetCpuTime() - initTime;
619   LOG(AliL3Log::kInformational,"AliL3Hough::FindTrackCandidates()","Timing")
620     <<"Maxima finding done in "<<cpuTime*1000<<" ms"<<ENDLOG;
621 }
622
623 void AliL3Hough::InitEvaluate()
624 {
625   //Pass the transformer objects to the AliL3HoughEval objects:
626   //This will provide the evaluation objects with all the necessary
627   //data and parameters it needs.
628   
629   for(Int_t i=0; i<fNPatches; i++) 
630     fEval[i]->InitTransformer(fHoughTransformer[i]);
631 }
632
633 Int_t AliL3Hough::Evaluate(Int_t road_width,Int_t nrowstomiss)
634 {
635   //Evaluate the tracks, by looking along the road in the raw data.
636   //If track does not cross all padrows - rows2miss, it is removed from the arrray.
637   //If histograms were not added, the check is done locally in patch,
638   //meaning that nrowstomiss is the number of padrows the road can miss with respect
639   //to the number of rows in the patch.
640   //If the histograms were added, the comparison is done globally in the _slice_, 
641   //meaing that nrowstomiss is the number of padrows the road can miss with
642   //respect to the total number of padrows in the slice.
643   //
644   //Return value = number of tracks which were removed (only in case of fAddHistograms)
645   
646   if(!fTracks[0])
647     {
648       LOG(AliL3Log::kError,"AliL3Hough::Evaluate","Track Array")
649         <<"No tracks to work with..."<<ENDLOG;
650       return 0;
651     }
652   
653   Int_t removed_tracks=0;
654   AliL3TrackArray *tracks=0;
655
656   if(fAddHistograms)
657     {
658       tracks = fTracks[0];
659       for(Int_t i=0; i<tracks->GetNTracks(); i++)
660         {
661           AliL3Track *track = tracks->GetCheckedTrack(i);
662           if(!track) continue;
663           track->SetNHits(0);
664         }
665     }
666   
667   for(Int_t i=0; i<fNPatches; i++)
668     EvaluatePatch(i,road_width,nrowstomiss);
669   
670   //Here we check the tracks globally; 
671   //how many good rows (padrows with signal) 
672   //did it cross in the slice
673   if(fAddHistograms) 
674     {
675       for(Int_t j=0; j<tracks->GetNTracks(); j++)
676         {
677           AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
678           
679           if(track->GetNHits() < AliL3Transform::GetNRows() - nrowstomiss)
680             {
681               tracks->Remove(j);
682               removed_tracks++;
683             }
684         }
685       tracks->Compress();
686       tracks->QSort();
687     }
688     
689   return removed_tracks;
690 }
691
692 void AliL3Hough::EvaluatePatch(Int_t i,Int_t road_width,Int_t nrowstomiss)
693 {
694   //Evaluate patch i.
695   
696   fEval[i]->InitTransformer(fHoughTransformer[i]);
697   fEval[i]->SetNumOfPadsToLook(road_width);
698   fEval[i]->SetNumOfRowsToMiss(nrowstomiss);
699   //fEval[i]->RemoveFoundTracks();
700   
701   AliL3TrackArray *tracks=0;
702   
703   if(!fAddHistograms)
704     tracks = fTracks[i];
705   else
706     tracks = fTracks[0];
707   
708   Int_t nrows=0;
709   for(Int_t j=0; j<tracks->GetNTracks(); j++)
710     {
711       AliL3HoughTrack *track = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
712       if(!track)
713         {
714           LOG(AliL3Log::kWarning,"AliL3Hough::EvaluatePatch","Track array")
715             <<"Track object missing!"<<ENDLOG;
716           continue;
717         } 
718       nrows=0;
719       Int_t rowrange[2] = {AliL3Transform::GetFirstRow(i),AliL3Transform::GetLastRow(i)};
720       Bool_t result = fEval[i]->LookInsideRoad(track,nrows,rowrange);
721       if(fAddHistograms)
722         {
723           Int_t pre=track->GetNHits();
724           track->SetNHits(pre+nrows);
725         }
726       else//the track crossed too few good padrows (padrows with signal) in the patch, so remove it
727         {
728           if(result == kFALSE)
729             tracks->Remove(j);
730         }
731     }
732   
733   tracks->Compress();
734
735 }
736
737 void AliL3Hough::MergeEtaSlices()
738 {
739   //Merge tracks found in neighbouring eta slices.
740   //Removes the track with the lower weight.
741   
742   fBenchmark->Start("Merge Eta-slices");
743   AliL3TrackArray *tracks = fTracks[0];
744   if(!tracks)
745     {
746       cerr<<"AliL3Hough::MergeEtaSlices : No tracks "<<endl;
747       return;
748     }
749   for(Int_t j=0; j<tracks->GetNTracks(); j++)
750     {
751       AliL3HoughTrack *track1 = (AliL3HoughTrack*)tracks->GetCheckedTrack(j);
752       if(!track1) continue;
753       for(Int_t k=j+1; k<tracks->GetNTracks(); k++)
754         {
755           AliL3HoughTrack *track2 = (AliL3HoughTrack*)tracks->GetCheckedTrack(k);
756           if(!track2) continue;
757           if(abs(track1->GetEtaIndex() - track2->GetEtaIndex()) != 1) continue;
758           if(fabs(track1->GetKappa()-track2->GetKappa()) < 0.006 && 
759              fabs(track1->GetPsi()- track2->GetPsi()) < 0.1)
760             {
761               //cout<<"Merging track in slices "<<track1->GetEtaIndex()<<" "<<track2->GetEtaIndex()<<endl;
762               if(track1->GetWeight() > track2->GetWeight())
763                 tracks->Remove(k);
764               else
765                 tracks->Remove(j);
766             }
767         }
768     }
769   fBenchmark->Stop("Merge Eta-slices");
770   tracks->Compress();
771 }
772
773 void AliL3Hough::WriteTracks(Char_t *path)
774 {
775   //cout<<"AliL3Hough::WriteTracks : Sorting the tracsk"<<endl;
776   //fGlobalTracks->QSort();
777   
778   Char_t filename[1024];
779   sprintf(filename,"%s/tracks_%d.raw",path,fEvent);
780   AliL3MemHandler mem;
781   mem.SetBinaryOutput(filename);
782   mem.TrackArray2Binary(fGlobalTracks);
783   mem.CloseBinaryOutput();
784   fGlobalTracks->Reset();
785 }
786
787 void AliL3Hough::WriteTracks(Int_t slice,Char_t *path)
788 {
789   
790   AliL3MemHandler mem;
791   Char_t fname[100];
792   if(fAddHistograms)
793     {
794       sprintf(fname,"%s/tracks_ho_%d_%d.raw",path,fEvent,slice);
795       mem.SetBinaryOutput(fname);
796       mem.TrackArray2Binary(fTracks[0]);
797       mem.CloseBinaryOutput();
798     }
799   else 
800     {
801       for(Int_t i=0; i<fNPatches; i++)
802         {
803           sprintf(fname,"%s/tracks_ho_%d_%d_%d.raw",path,fEvent,slice,i);
804           mem.SetBinaryOutput(fname);
805           mem.TrackArray2Binary(fTracks[i]);
806           mem.CloseBinaryOutput();
807         }
808     }
809 }
810
811 void AliL3Hough::WriteDigits(Char_t *outfile)
812 {
813 #ifdef use_aliroot  
814   //Write the current data to a new rootfile.
815
816   for(Int_t i=0; i<fNPatches; i++)
817     {
818       AliL3DigitRowData *tempPt = (AliL3DigitRowData*)fHoughTransformer[i]->GetDataPointer();
819       fMemHandler[i]->AliDigits2RootFile(tempPt,outfile);
820     }
821 #else
822   cerr<<"AliL3Hough::WriteDigits : You need to compile with AliROOT!"<<endl;
823   return;
824 #endif  
825 }
826
827 Double_t AliL3Hough::GetCpuTime()
828 {
829   //Return the Cputime in seconds.
830  struct timeval tv;
831  gettimeofday( &tv, NULL );
832  return tv.tv_sec+(((Double_t)tv.tv_usec)/1000000.);
833 }
834