New class AliESDEvent, backward compatibility with the old AliESD (Christian)
[u/mrichter/AliRoot.git] / HLT / hough / AliHLTHough.cxx
1 // @(#) $Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ALICE HLT Group
5
6 #include "AliHLTStandardIncludes.h"
7 #include <sys/time.h>
8
9 #include "AliHLTLogging.h"
10 #include "AliHLTHoughMerger.h"
11 #include "AliHLTHoughIntMerger.h"
12 #include "AliHLTHoughGlobalMerger.h"
13 #include "AliHLTHistogram.h"
14 #include "AliHLTHough.h"
15 #include "AliHLTHoughTransformer.h"
16 #ifdef HAVE_ALIHLTHOUGHCLUSTERTRANSFORMER
17 #include "AliHLTHoughClusterTransformer.h"
18 #endif //HAVE_ALIHLTHOUGHCLUSTERTRANSFORMER
19 #include "AliHLTHoughTransformerLUT.h"
20 #include "AliHLTHoughTransformerVhdl.h"
21 #include "AliHLTHoughTransformerRow.h"
22 #include "AliHLTHoughMaxFinder.h"
23 #include "AliHLTBenchmark.h"
24 #ifdef use_aliroot
25 #include "AliHLTFileHandler.h"
26 #else
27 #include "AliHLTMemHandler.h"
28 #endif
29 #include "AliHLTDataHandler.h"
30 #include "AliHLTDigitData.h"
31 #include "AliHLTHoughEval.h"
32 #include "AliHLTTransform.h"
33 #include "AliHLTTrackArray.h"
34 #include "AliHLTHoughTrack.h"
35 #include "AliHLTDDLDataFileHandler.h"
36 #include "AliHLTHoughKalmanTrack.h"
37
38 #include "TThread.h"
39
40 #if __GNUC__ >= 3
41 using namespace std;
42 #endif
43
44 /** /class AliHLTHough
45 //<pre>
46 //_____________________________________________________________
47 // AliHLTHough
48 //
49 // Interface class for the Hough transform
50 //
51 // Example how to use:
52 //
53 // AliHLTHough *hough = new AliHLTHough(path,kTRUE,NumberOfEtaSegments);
54 // hough->ReadData(slice);
55 // hough->Transform();
56 // hough->FindTrackCandidates();
57 // 
58 // AliHLTTrackArray *tracks = hough->GetTracks(patch);
59 //
60 //</pre>
61 */
62
63 ClassImp(AliHLTHough)
64
65 AliHLTHough::AliHLTHough()
66 {
67   //Constructor
68   
69   fBinary        = kFALSE;
70   fAddHistograms = kFALSE;
71   fDoIterative   = kFALSE; 
72   fWriteDigits   = kFALSE;
73   fUse8bits      = kFALSE;
74
75   fMemHandler       = 0;
76   fHoughTransformer = 0;
77   fEval             = 0;
78   fPeakFinder       = 0;
79   fTracks           = 0;
80   fGlobalTracks     = 0;
81   fMerger           = 0;
82   fInterMerger      = 0;
83   fGlobalMerger     = 0;
84   fBenchmark        = 0;
85   
86   fNEtaSegments     = 0;
87   fNPatches         = 0;
88   fLastPatch        =-1;
89   fVersion          = 0;
90   fCurrentSlice     = 0;
91   fEvent            = 0;
92   
93   fKappaSpread = 6;
94   fPeakRatio   = 0.5;
95   fInputFile   = 0;
96   fInputPtr    = 0;
97   fRawEvent    = 0;
98   
99   SetTransformerParams();
100   SetThreshold();
101   SetNSaveIterations();
102   SetPeakThreshold();
103 #ifdef use_aliroot
104   //just be sure that index is empty for new event
105     AliHLTFileHandler::CleanStaticIndex(); 
106 #ifdef use_newio
107     fRunLoader = 0;
108 #endif
109 #endif
110   fThread = 0;
111 }
112
113 AliHLTHough::AliHLTHough(Char_t *path,Bool_t binary,Int_t netasegments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr)
114 {
115   //Normal constructor
116   fBinary = binary;
117   strcpy(fPath,path);
118   fNEtaSegments  = netasegments;
119   fAddHistograms = kFALSE;
120   fDoIterative   = kFALSE; 
121   fWriteDigits   = kFALSE;
122   fUse8bits      = bit8;
123   fVersion       = tv;
124   fKappaSpread=6;
125   fPeakRatio=0.5;
126   if(!fBinary) {
127     if(infile) {
128       fInputFile = infile;
129       fInputPtr = 0;
130     }
131     else {
132       fInputFile = 0;
133       fInputPtr = ptr;
134     }
135   }
136   else {
137     fInputFile = 0;
138     fInputPtr = 0;
139   }
140 #ifdef use_aliroot
141   //just be sure that index is empty for new event
142     AliHLTFileHandler::CleanStaticIndex(); 
143 #ifdef use_newio
144     fRunLoader = 0;
145 #endif
146 #endif
147   fThread = 0;
148 }
149
150 AliHLTHough::~AliHLTHough()
151 {
152   //dtor
153
154   CleanUp();
155   if(fMerger)
156     delete fMerger;
157   //cout << "Cleaned class merger " << endl;
158   if(fInterMerger)
159     delete fInterMerger;
160   //cout << "Cleaned class inter " << endl;
161   if(fPeakFinder)
162     delete fPeakFinder;
163   //cout << "Cleaned class peak " << endl;
164   if(fGlobalMerger)
165     delete fGlobalMerger;
166   //cout << "Cleaned class global " << endl;
167   if(fBenchmark)
168     delete fBenchmark;
169   //cout << "Cleaned class bench " << endl;
170   if(fGlobalTracks)
171     delete fGlobalTracks;
172   //cout << "Cleaned class globaltracks " << endl;
173   if(fThread) {
174     //    fThread->Delete();
175     delete fThread;
176     fThread = 0;
177   }
178 }
179
180 void AliHLTHough::CleanUp()
181 {
182   //Cleanup memory
183   
184   for(Int_t i=0; i<fNPatches; i++)
185     {
186       if(fTracks[i]) delete fTracks[i];
187       //cout << "Cleaned tracks " << i << endl;
188       if(fEval[i]) delete fEval[i];
189       //cout << "Cleaned eval " << i << endl;
190       if(fHoughTransformer[i]) delete fHoughTransformer[i];
191       //cout << "Cleaned traf " << i << endl;
192       if(fMemHandler[i]) delete fMemHandler[i];
193       //cout << "Cleaned mem " << i << endl;
194     }
195   
196   if(fTracks) delete [] fTracks;
197   //cout << "Cleaned class tracks " << endl;
198   if(fEval) delete [] fEval;
199   //cout << "Cleaned class eval " << endl;
200   if(fHoughTransformer) delete [] fHoughTransformer;
201   //cout << "Cleaned cleass trafo " << endl;
202   if(fMemHandler) delete [] fMemHandler;
203   //cout << "Cleaned class mem " << endl;
204 }
205
206 void AliHLTHough::Init(Int_t netasegments,Int_t tv,AliRawEvent *rawevent,Float_t zvertex)
207 {
208   //Normal constructor
209   fNEtaSegments  = netasegments;
210   fVersion       = tv;
211   fRawEvent      = rawevent;
212   fZVertex       = zvertex;
213
214   Init();
215 }
216
217 void AliHLTHough::Init(Char_t *path,Bool_t binary,Int_t netasegments,Bool_t bit8,Int_t tv,Char_t *infile,Char_t *ptr,Float_t zvertex)
218 {
219   //Normal init of the AliHLTHough
220   fBinary = binary;
221   strcpy(fPath,path);
222   fNEtaSegments = netasegments;
223   fWriteDigits  = kFALSE;
224   fUse8bits     = bit8;
225   fVersion      = tv;
226   if(!fBinary) {
227     if(infile) {
228       fInputFile = infile;
229       fInputPtr = 0;
230     }
231     else {
232       fInputFile = 0;
233       fInputPtr = ptr;
234     }
235   }
236   else {
237     fInputFile = 0;
238     fInputPtr = 0;
239   }
240   fZVertex = zvertex;
241
242   Init(); //do the rest
243 }
244
245 void AliHLTHough::Init(Bool_t doit, Bool_t addhists)
246 {
247   // Init
248   fDoIterative   = doit; 
249   fAddHistograms = addhists;
250
251   fNPatches = AliHLTTransform::GetNPatches();
252   fHoughTransformer = new AliHLTHoughBaseTransformer*[fNPatches];
253   fMemHandler = new AliHLTMemHandler*[fNPatches];
254
255   fTracks = new AliHLTTrackArray*[fNPatches];
256   fEval = new AliHLTHoughEval*[fNPatches];
257   
258   fGlobalTracks = new AliHLTTrackArray("AliHLTHoughTrack");
259   
260   AliHLTHoughBaseTransformer *lasttransformer = 0;
261
262   for(Int_t i=0; i<fNPatches; i++)
263     {
264       switch (fVersion){ //choose Transformer
265       case 1: 
266         fHoughTransformer[i] = new AliHLTHoughTransformerLUT(0,i,fNEtaSegments);
267         break;
268       case 2:
269 #ifdef HAVE_ALIHLTHOUGHCLUSTERTRANSFORMER
270         fHoughTransformer[i] = new AliHLTHoughClusterTransformer(0,i,fNEtaSegments);
271 #endif //HAVE_ALIHLTHOUGHCLUSTERTRANSFORMER
272         break;
273       case 3:
274         fHoughTransformer[i] = new AliHLTHoughTransformerVhdl(0,i,fNEtaSegments,fNSaveIterations);
275         break;
276       case 4:
277         fHoughTransformer[i] = new AliHLTHoughTransformerRow(0,i,fNEtaSegments,kFALSE,fZVertex);
278         break;
279       default:
280         fHoughTransformer[i] = new AliHLTHoughTransformer(0,i,fNEtaSegments,kFALSE,kFALSE);
281       }
282
283       fHoughTransformer[i]->SetLastTransformer(lasttransformer);
284       lasttransformer = fHoughTransformer[i];
285       //      fHoughTransformer[i]->CreateHistograms(fNBinX[i],fLowPt[i],fNBinY[i],-fPhi[i],fPhi[i]);
286       fHoughTransformer[i]->CreateHistograms(fNBinX[i],-fLowPt[i],fLowPt[i],fNBinY[i],-fPhi[i],fPhi[i]);
287       //fHoughTransformer[i]->CreateHistograms(fLowPt[i],fUpperPt[i],fPtRes[i],fNBinY[i],fPhi[i]);
288
289       fHoughTransformer[i]->SetLowerThreshold(fThreshold[i]);
290       fHoughTransformer[i]->SetUpperThreshold(100);
291
292       LOG(AliHLTLog::kInformational,"AliHLTHough::Init","Version")
293         <<"Initializing Hough transformer version "<<fVersion<<ENDLOG;
294       
295       fEval[i] = new AliHLTHoughEval();
296       fTracks[i] = new AliHLTTrackArray("AliHLTHoughTrack");
297       if(fUse8bits)
298         fMemHandler[i] = new AliHLTDataHandler();
299       else
300 #ifdef use_aliroot
301         {
302           if(!fRawEvent) {
303             if(!fInputFile) {
304               if(!fInputPtr) {
305                 /* In case of reading digits file */
306                 fMemHandler[i] = new AliHLTFileHandler(kTRUE); //use static index
307                 if(!fBinary) {
308 #if use_newio
309                   if(!fRunLoader) {
310 #endif
311                     Char_t filename[1024];
312                     sprintf(filename,"%s/digitfile.root",fPath);
313                     fMemHandler[i]->SetAliInput(filename);
314 #if use_newio
315                   }
316                   else {
317                     fMemHandler[i]->SetAliInput(fRunLoader);
318                   }
319 #endif
320                 }
321               }
322               else {
323                 /* In case of reading from DATE */
324                 fMemHandler[i] = new AliHLTDDLDataFileHandler();
325                 fMemHandler[i]->SetReaderInput(fInputPtr,-1);
326               }
327             }
328             else {
329               /* In case of reading rawdata from ROOT file */
330               fMemHandler[i] = new AliHLTDDLDataFileHandler();
331               fMemHandler[i]->SetReaderInput(fInputFile);
332             }
333           }
334           else {
335             /* In case of reading rawdata using AliRawEvent */
336             fMemHandler[i] = new AliHLTDDLDataFileHandler();
337             fMemHandler[i]->SetReaderInput(fRawEvent);
338           }
339         }
340 #else
341       fMemHandler[i] = new AliHLTMemHandler();
342 #endif
343     }
344
345   fPeakFinder = new AliHLTHoughMaxFinder("KappaPhi",50000);
346   if(fVersion!=4) {
347     fMerger = new AliHLTHoughMerger(fNPatches);
348     fInterMerger = new AliHLTHoughIntMerger();
349   }
350   else {
351     fMerger = 0;
352     fInterMerger = 0;
353   }
354   fGlobalMerger = 0;
355   fBenchmark = new AliHLTBenchmark();
356 }
357
358 void AliHLTHough::SetTransformerParams(Float_t ptres,Float_t ptmin,Float_t ptmax,Int_t ny,Int_t patch)
359 {
360   // Setup the parameters for the Hough Transformer
361   // This includes the bin size and limits for
362   // the parameter space histograms
363
364   Int_t mrow;
365   Float_t psi=0;
366   if(patch==-1)
367     mrow = 80;
368   else
369     mrow = AliHLTTransform::GetLastRow(patch);
370   if(ptmin)
371     {
372       Double_t lineradius = sqrt(pow(AliHLTTransform::Row2X(mrow),2) + pow(AliHLTTransform::GetMaxY(mrow),2));
373       Double_t kappa = -1*AliHLTTransform::GetBField()*AliHLTTransform::GetBFact()/ptmin;
374       psi = AliHLTTransform::Deg2Rad(10) - asin(lineradius*kappa/2);
375       cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
376     }
377
378   if(patch==-1)
379     {
380       Int_t i=0;
381       while(i < 6)
382         {
383           fPtRes[i] = ptres;
384           fLowPt[i] = ptmin;
385           fUpperPt[i] = ptmax;
386           fNBinY[i] = ny;
387           fPhi[i] = psi;
388           fNBinX[i]=0;
389           i++;
390         }
391       return;
392     }
393
394   fPtRes[patch] = ptres;
395   fLowPt[patch] = ptmin;
396   fUpperPt[patch] = ptmax;
397   fNBinY[patch] = ny;
398   fPhi[patch] = psi;
399 }
400 /*
401 void AliHLTHough::SetTransformerParams(Int_t nx,Int_t ny,Float_t ptmin,Int_t patch)
402 {
403   // Setup the parameters for the Hough Transformer
404
405   Int_t mrow=80;
406   Double_t lineradius = sqrt(pow(AliHLTTransform::Row2X(mrow),2) + pow(AliHLTTransform::GetMaxY(mrow),2));
407   Double_t kappa = -1*AliHLTTransform::GetBField()*AliHLTTransform::GetBFact()/ptmin;
408   Double_t psi = AliHLTTransform::Deg2Rad(10) - asin(lineradius*kappa/2);
409   cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
410   
411   Int_t i=0;
412   while(i < 6)
413     {
414       fLowPt[i] = ptmin;
415       fNBinY[i] = ny;
416       fNBinX[i] = nx;
417       fPhi[i] = psi;
418       i++;
419     }
420 }
421 */
422 void AliHLTHough::SetTransformerParams(Int_t nx,Int_t ny,Float_t ptmin,Int_t /*patch*/)
423 {
424   // Setup the parameters for the Hough Transformer
425
426   Double_t lineradius = 1.0/(AliHLTHoughTransformerRow::GetBeta1()*sqrt(1.0+tan(AliHLTTransform::Pi()*10/180)*tan(AliHLTTransform::Pi()*10/180)));
427   Double_t alpha1 = AliHLTHoughTransformerRow::GetBeta1()*tan(AliHLTTransform::Pi()*10/180);
428   Double_t kappa = 1*AliHLTTransform::GetBField()*AliHLTTransform::GetBFact()/(ptmin*0.9);
429   Double_t psi = AliHLTTransform::Deg2Rad(10) - asin(lineradius*kappa/2);
430   //  cout<<"Calculated psi range "<<psi<<" in patch "<<patch<<endl;
431   Double_t alpha2 = alpha1 - (AliHLTHoughTransformerRow::GetBeta1()-AliHLTHoughTransformerRow::GetBeta2())*tan(psi);
432   //  cout<<"Calculated alphas range "<<alpha1<<" "<<alpha2<<" in patch "<<patch<<endl;
433
434   Int_t i=0;
435   while(i < 6)
436     {
437       fLowPt[i] = 1.1*alpha1;
438       fNBinY[i] = ny;
439       fNBinX[i] = nx;
440       fPhi[i] = alpha2;
441       i++;
442     }
443 }
444
445 void AliHLTHough::CalcTransformerParams(Float_t ptmin)
446 {
447   // Setup the parameters for the Row Hough Transformer
448   // Automatically adjusts the number of bins in X and Y in a way
449   // that the size of the hough bin is 2x (in X) and 2.5 (in Y) the
450   // size of the tpc pads
451
452   Double_t lineradius = 1.0/(AliHLTHoughTransformerRow::GetBeta1()*sqrt(1.0+tan(AliHLTTransform::Pi()*10/180)*tan(AliHLTTransform::Pi()*10/180)));
453   Double_t alpha1 = AliHLTHoughTransformerRow::GetBeta1()*tan(AliHLTTransform::Pi()*10/180);
454   Double_t kappa = 1*AliHLTTransform::GetBField()*AliHLTTransform::GetBFact()/(ptmin*0.9);
455   Double_t psi = AliHLTTransform::Deg2Rad(10) - asin(lineradius*kappa/2);
456   //  cout<<"Calculated psi range "<<psi<<endl;
457   Double_t alpha2 = alpha1 - (AliHLTHoughTransformerRow::GetBeta1()-AliHLTHoughTransformerRow::GetBeta2())*tan(psi);
458   alpha1 *= 1.1;
459   //  cout<<"Calculated alphas range "<<alpha1<<" "<<alpha2<<endl;
460
461   Double_t sizex = 2.0*AliHLTTransform::GetPadPitchWidthLow()*AliHLTHoughTransformerRow::GetBeta1()*AliHLTHoughTransformerRow::GetBeta1();
462   Double_t sizey = 2.5*AliHLTTransform::GetPadPitchWidthUp()*AliHLTHoughTransformerRow::GetBeta2()*AliHLTHoughTransformerRow::GetBeta2();
463
464   Int_t nx = 2*(Int_t)(alpha1/sizex)+1;
465   Int_t ny = 2*(Int_t)(alpha2/sizey)+1;
466   //  cout<<"Calculated number of bins "<<nx<<" "<<ny<<endl;
467
468   Int_t i=0;
469   while(i < 6)
470     {
471       fLowPt[i] = alpha1;
472       fNBinY[i] = ny;
473       fNBinX[i] = nx;
474       fPhi[i] = alpha2;
475       i++;
476     }
477 }
478
479 void AliHLTHough::SetTransformerParams(Int_t nx,Int_t ny,Float_t lpt,Float_t phi)
480 {
481   Int_t i=0;
482   while(i < 6)
483     {
484       fLowPt[i] = lpt;
485       fNBinY[i] = ny;
486       fNBinX[i] = nx;
487       fPhi[i] = phi;
488       i++;
489     }
490 }
491
492 void AliHLTHough::SetThreshold(Int_t t3,Int_t patch)
493 {
494   // Set digits threshold
495   if(patch==-1)
496     {
497       Int_t i=0;
498       while(i < 6)
499         fThreshold[i++]=t3;
500       return;
501     }
502   fThreshold[patch]=t3;
503 }
504
505 void AliHLTHough::SetPeakThreshold(Int_t threshold,Int_t patch)
506 {
507   // Set Peak Finder threshold
508   if(patch==-1)
509     {
510       Int_t i=0;
511       while(i < 6)
512         fPeakThreshold[i++]=threshold;
513       return;
514     }
515   fPeakThreshold[patch]=threshold;
516 }
517
518 void AliHLTHough::DoBench(Char_t *name)
519 {
520   fBenchmark->Analyze(name);
521 }
522
523 void AliHLTHough::Process(Int_t minslice,Int_t maxslice)
524 {
525   //Process all slices [minslice,maxslice].
526   fGlobalMerger = new AliHLTHoughGlobalMerger(minslice,maxslice);
527   
528   for(Int_t i=minslice; i<=maxslice; i++)
529     {
530       ReadData(i);
531       Transform();
532       if(fAddHistograms) {
533         if(fVersion != 4)
534           AddAllHistograms();
535         else
536           AddAllHistogramsRows();
537       }
538       FindTrackCandidates();
539       //Evaluate();
540       //fGlobalMerger->FillTracks(fTracks[0],i);
541     }
542 }
543
544 void AliHLTHough::ReadData(Int_t slice,Int_t eventnr)
545 {
546   //Read data from files, binary or root.
547   
548 #ifdef use_aliroot
549   if(fEvent!=eventnr) //just be sure that index is empty for new event
550     AliHLTFileHandler::CleanStaticIndex(); 
551 #endif
552   fCurrentSlice = slice;
553
554   for(Int_t i=0; i<fNPatches; i++)
555     {
556       fMemHandler[i]->Free();
557       UInt_t ndigits=0;
558       AliHLTDigitRowData *digits =0;
559       Char_t name[256];
560       fMemHandler[i]->Init(slice,i);
561       if(fBinary)//take input data from binary files
562         {
563           if(fUse8bits)
564             sprintf(name,"%s/binaries/digits_c8_%d_%d_%d.raw",fPath,eventnr,slice,i);
565           else
566             sprintf(name,"%s/binaries/digits_%d_%d_%d.raw",fPath,eventnr,slice,i);
567
568           fMemHandler[i]->SetBinaryInput(name);
569           digits = (AliHLTDigitRowData *)fMemHandler[i]->CompBinary2Memory(ndigits);
570           fMemHandler[i]->CloseBinaryInput();
571         }
572       else //read data from root file
573         {
574 #ifdef use_aliroot
575           if(fEvent!=eventnr)
576             fMemHandler[i]->FreeDigitsTree();//or else the new event is not loaded
577           digits=(AliHLTDigitRowData *)fMemHandler[i]->AliAltroDigits2Memory(ndigits,eventnr);
578 #else
579           cerr<<"You cannot read from rootfile now"<<endl;
580 #endif
581         }
582
583       //Set the pointer to the TPCRawStream in case of fast raw data reading
584       fHoughTransformer[i]->SetTPCRawStream(fMemHandler[i]->GetTPCRawStream());
585
586       //set input data and init transformer
587       fHoughTransformer[i]->SetInputData(ndigits,digits);
588       fHoughTransformer[i]->Init(slice,i,fNEtaSegments);
589     }
590
591   fEvent=eventnr;
592 }
593
594 void AliHLTHough::Transform(Int_t *rowrange)
595 {
596   //Transform all data given to the transformer within the given slice
597   //(after ReadData(slice))
598   
599   Double_t initTime,cpuTime;
600   initTime = GetCpuTime();
601   Int_t patchorder[6] = {5,2,0,1,3,4}; //The order in which patches are processed
602   //  Int_t patchorder[6] = {0,1,2,3,4,5}; //The order in which patches are processed
603   //  Int_t patchorder[6] = {5,4,3,2,1,0}; //The order in which patches are processed
604   //  Int_t patchorder[6] = {5,2,4,3,1,0}; //The order in which patches are processed
605   fLastPatch=-1;
606   for(Int_t i=0; i<fNPatches; i++)
607     {
608       // In case of Row transformer reset the arrays only once
609       if((fVersion != 4) || (i == 0)) {
610         fBenchmark->Start("Hough Reset");
611         fHoughTransformer[0]->Reset();//Reset the histograms
612         fBenchmark->Stop("Hough Reset");
613       }
614       fBenchmark->Start("Hough Transform");
615       PrepareForNextPatch(patchorder[i]);
616       if(!rowrange) {
617         char buf[256];
618         sprintf(buf,"Patch %d",patchorder[i]);
619         fBenchmark->Start(buf);
620         fHoughTransformer[patchorder[i]]->SetLastPatch(fLastPatch);
621         fHoughTransformer[patchorder[i]]->TransformCircle();
622         fBenchmark->Stop(buf);
623       }
624       else
625         fHoughTransformer[i]->TransformCircleC(rowrange,1);
626       fBenchmark->Stop("Hough Transform");
627       fLastPatch=patchorder[i];
628     }
629   cpuTime = GetCpuTime() - initTime;
630   LOG(AliHLTLog::kInformational,"AliHLTHough::Transform()","Timing")
631     <<"Transform done in average per patch of "<<cpuTime*1000/fNPatches<<" ms"<<ENDLOG;
632 }
633
634 void AliHLTHough::MergePatches()
635 {
636   // Merge patches if they are not summed
637   if(fAddHistograms) //Nothing to merge here
638     return;
639   fMerger->MergePatches(kTRUE);
640 }
641
642 void AliHLTHough::MergeInternally()
643 {
644   // Merge patches internally
645   if(fAddHistograms)
646     fInterMerger->FillTracks(fTracks[0]);
647   else
648     fInterMerger->FillTracks(fMerger->GetOutTracks());
649   
650   fInterMerger->MMerge();
651 }
652
653 void AliHLTHough::ProcessSliceIter()
654 {
655   //Process current slice (after ReadData(slice)) iteratively.
656   
657   if(!fAddHistograms)
658     {
659       for(Int_t i=0; i<fNPatches; i++)
660         {
661           ProcessPatchIter(i);
662           fMerger->FillTracks(fTracks[i],i); //Copy tracks to merger
663         }
664     }
665   else
666     {
667       for(Int_t i=0; i<10; i++)
668         {
669           Transform();
670           AddAllHistograms();
671           InitEvaluate();
672           AliHLTHoughBaseTransformer *tr = fHoughTransformer[0];
673           for(Int_t j=0; j<fNEtaSegments; j++)
674             {
675               AliHLTHistogram *hist = tr->GetHistogram(j);
676               if(hist->GetNEntries()==0) continue;
677               fPeakFinder->Reset();
678               fPeakFinder->SetHistogram(hist);
679               fPeakFinder->FindAbsMaxima();
680               AliHLTHoughTrack *track = (AliHLTHoughTrack*)fTracks[0]->NextTrack();
681               track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
682               track->SetEtaIndex(j);
683               track->SetEta(tr->GetEta(j,fCurrentSlice));
684               for(Int_t k=0; k<fNPatches; k++)
685                 {
686                   fEval[i]->SetNumOfPadsToLook(2);
687                   fEval[i]->SetNumOfRowsToMiss(2);
688                   fEval[i]->RemoveFoundTracks();
689                   /*
690                   Int_t nrows=0;
691                   if(!fEval[i]->LookInsideRoad(track,nrows))
692                     {
693                       fTracks[0]->Remove(fTracks[0]->GetNTracks()-1);
694                       fTracks[0]->Compress();
695                     }
696                   */
697                 }
698             }
699           
700         }
701       
702     }
703 }
704
705 void AliHLTHough::ProcessPatchIter(Int_t patch)
706 {
707   //Process patch in a iterative way. 
708   //transform + peakfinding + evaluation + transform +...
709
710   Int_t numoftries = 5;
711   AliHLTHoughBaseTransformer *tr = fHoughTransformer[patch];
712   AliHLTTrackArray *tracks = fTracks[patch];
713   tracks->Reset();
714   AliHLTHoughEval *ev = fEval[patch];
715   ev->InitTransformer(tr);
716   //ev->RemoveFoundTracks();
717   ev->SetNumOfRowsToMiss(3);
718   ev->SetNumOfPadsToLook(2);
719   AliHLTHistogram *hist;
720   for(Int_t t=0; t<numoftries; t++)
721     {
722       tr->Reset();
723       tr->TransformCircle();
724       for(Int_t i=0; i<fNEtaSegments; i++)
725         {
726           hist = tr->GetHistogram(i);
727           if(hist->GetNEntries()==0) continue;
728           fPeakFinder->Reset();
729           fPeakFinder->SetHistogram(hist);
730           fPeakFinder->FindAbsMaxima();
731           //fPeakFinder->FindPeak1();
732           AliHLTHoughTrack *track = (AliHLTHoughTrack*)tracks->NextTrack();
733           track->SetTrackParameters(fPeakFinder->GetXPeak(0),fPeakFinder->GetYPeak(0),fPeakFinder->GetWeight(0));
734           track->SetEtaIndex(i);
735           track->SetEta(tr->GetEta(i,fCurrentSlice));
736           /*
737           Int_t nrows=0;
738           if(!ev->LookInsideRoad(track,nrows))
739             {   
740               tracks->Remove(tracks->GetNTracks()-1);
741               tracks->Compress();
742             }
743           */
744         }
745     }
746   fTracks[0]->QSort();
747   LOG(AliHLTLog::kInformational,"AliHLTHough::ProcessPatch","NTracks")
748     <<AliHLTLog::kDec<<"Found "<<tracks->GetNTracks()<<" tracks in patch "<<patch<<ENDLOG;
749 }
750
751 void AliHLTHough::AddAllHistograms()
752 {
753   //Add the histograms within one etaslice.
754   //Resulting histogram are in patch=0.
755
756   Double_t initTime,cpuTime;
757   initTime = GetCpuTime();
758   fBenchmark->Start("Add Histograms");
759   for(Int_t i=0; i<fNEtaSegments; i++)
760     {
761       AliHLTHistogram *hist0 = fHoughTransformer[0]->GetHistogram(i);
762       for(Int_t j=1; j<fNPatches; j++)
763         {
764           AliHLTHistogram *hist = fHoughTransformer[j]->GetHistogram(i);
765           hist0->Add(hist);
766         }
767     }
768   fBenchmark->Stop("Add Histograms");
769   fAddHistograms = kTRUE;
770   cpuTime = GetCpuTime() - initTime;
771   LOG(AliHLTLog::kInformational,"AliHLTHough::AddAllHistograms()","Timing")
772     <<"Adding histograms in "<<cpuTime*1000<<" ms"<<ENDLOG;
773 }
774
775 void AliHLTHough::AddAllHistogramsRows()
776 {
777   //Add the histograms within one etaslice.
778   //Resulting histogram are in patch=0.
779
780   Double_t initTime,cpuTime;
781   initTime = GetCpuTime();
782   fBenchmark->Start("Add HistogramsRows");
783
784   UChar_t lastpatchlastrow = AliHLTTransform::GetLastRowOnDDL(fLastPatch)+1;
785
786   UChar_t *tracklastrow = ((AliHLTHoughTransformerRow *)fHoughTransformer[0])->GetTrackLastRow();
787
788   for(Int_t i=0; i<fNEtaSegments; i++)
789     {
790       UChar_t *gapcount = ((AliHLTHoughTransformerRow *)fHoughTransformer[0])->GetGapCount(i);
791       UChar_t *currentrowcount = ((AliHLTHoughTransformerRow *)fHoughTransformer[0])->GetCurrentRowCount(i);
792
793       AliHLTHistogram *hist = fHoughTransformer[0]->GetHistogram(i);
794       Int_t xmin = hist->GetFirstXbin();
795       Int_t xmax = hist->GetLastXbin();
796       Int_t ymin = hist->GetFirstYbin();
797       Int_t ymax = hist->GetLastYbin();
798       Int_t nxbins = hist->GetNbinsX()+2;
799
800       for(Int_t ybin=ymin; ybin<=ymax; ybin++)
801         {
802           for(Int_t xbin=xmin; xbin<=xmax; xbin++)
803             {
804               Int_t bin = xbin + ybin*nxbins; //Int_t bin = hist->GetBin(xbin,ybin);
805               if(gapcount[bin] < MAX_N_GAPS) {
806                 if(tracklastrow[bin] > lastpatchlastrow) {
807                   if(lastpatchlastrow > currentrowcount[bin])
808                     gapcount[bin] += (lastpatchlastrow-currentrowcount[bin]-1);
809                 }
810                 else {
811                   if(tracklastrow[bin] > currentrowcount[bin])
812                     gapcount[bin] += (tracklastrow[bin]-currentrowcount[bin]-1);
813                 }
814                 if(gapcount[bin] < MAX_N_GAPS)
815                   hist->AddBinContent(bin,(159-gapcount[bin]));
816               }
817             }
818         }
819     }
820
821   fBenchmark->Stop("Add HistogramsRows");
822   fAddHistograms = kTRUE;
823   cpuTime = GetCpuTime() - initTime;
824   LOG(AliHLTLog::kInformational,"AliHLTHough::AddAllHistogramsRows()","Timing")
825     <<"Adding histograms in "<<cpuTime*1000<<" ms"<<ENDLOG;
826 }
827
828 void AliHLTHough::PrepareForNextPatch(Int_t nextpatch)
829 {
830   // Prepare the parameter space for the processing of
831   // the next read patch. According to the already
832   // accumulated number of gaps in parameter space
833   // bins, the routine updates the dynamic
834   // pointers used in order to jump rapidly during the
835   // filling of the parameter space.
836
837   char buf[256];
838   sprintf(buf,"Prepare For Patch %d",nextpatch);
839   fBenchmark->Start(buf);
840
841   UChar_t lastpatchlastrow;
842   if(fLastPatch == -1)
843     lastpatchlastrow = 0;
844   else
845     lastpatchlastrow = AliHLTTransform::GetLastRowOnDDL(fLastPatch)+1;
846   UChar_t nextpatchfirstrow;
847   if(nextpatch==0)
848     nextpatchfirstrow = 0;
849   else
850     nextpatchfirstrow = AliHLTTransform::GetFirstRowOnDDL(nextpatch)-1;
851
852   UChar_t *trackfirstrow = ((AliHLTHoughTransformerRow *)fHoughTransformer[0])->GetTrackFirstRow();
853   UChar_t *tracklastrow = ((AliHLTHoughTransformerRow *)fHoughTransformer[0])->GetTrackLastRow();
854
855   for(Int_t i=0; i<fNEtaSegments; i++)
856     {
857       UChar_t *gapcount = ((AliHLTHoughTransformerRow *)fHoughTransformer[0])->GetGapCount(i);
858       UChar_t *currentrowcount = ((AliHLTHoughTransformerRow *)fHoughTransformer[0])->GetCurrentRowCount(i);
859       UChar_t *prevbin = ((AliHLTHoughTransformerRow *)fHoughTransformer[0])->GetPrevBin(i);
860       UChar_t *nextbin = ((AliHLTHoughTransformerRow *)fHoughTransformer[0])->GetNextBin(i);
861       UChar_t *nextrow = ((AliHLTHoughTransformerRow *)fHoughTransformer[0])->GetNextRow(i);
862
863       AliHLTHistogram *hist = fHoughTransformer[0]->GetHistogram(i);
864       Int_t xmin = hist->GetFirstXbin();
865       Int_t xmax = hist->GetLastXbin();
866       Int_t ymin = hist->GetFirstYbin();
867       Int_t ymax = hist->GetLastYbin();
868       Int_t nxbins = hist->GetNbinsX()+2;
869
870       if(fLastPatch != -1) {
871         UChar_t lastyvalue = 0;
872         Int_t endybin = ymin - 1;
873         for(Int_t ybin=nextrow[ymin]; ybin<=ymax; ybin = nextrow[++ybin])
874           {
875             UChar_t lastxvalue = 0;
876             UChar_t maxvalue = 0;
877             Int_t endxbin = xmin - 1;
878             for(Int_t xbin=xmin; xbin<=xmax; xbin++)
879               {
880                 Int_t bin = xbin + ybin*nxbins;
881                 UChar_t value = 0;
882                 if(gapcount[bin] < MAX_N_GAPS) {
883                   if(tracklastrow[bin] > lastpatchlastrow) {
884                     if(lastpatchlastrow > currentrowcount[bin])
885                       gapcount[bin] += (lastpatchlastrow-currentrowcount[bin]-1);
886                   }
887                   else {
888                     if(tracklastrow[bin] > currentrowcount[bin])
889                       gapcount[bin] += (tracklastrow[bin]-currentrowcount[bin]-1);
890                   }
891                   if(gapcount[bin] < MAX_N_GAPS) {
892                     value = 1;
893                     maxvalue = 1;
894                     if(trackfirstrow[bin] < nextpatchfirstrow)
895                       currentrowcount[bin] = nextpatchfirstrow;
896                     else
897                       currentrowcount[bin] = trackfirstrow[bin];
898                   }
899                 }
900                 if(value > 0)
901                   {
902                     nextbin[xbin + ybin*nxbins] = (UChar_t)xbin;
903                     prevbin[xbin + ybin*nxbins] = (UChar_t)xbin;
904                     if(value > lastxvalue)
905                       {
906                         UChar_t *tempnextbin = nextbin + endxbin + 1 + ybin*nxbins;
907                         memset(tempnextbin,(UChar_t)xbin,xbin-endxbin-1);
908                       }
909                     endxbin = xbin;
910                   }
911                 else
912                   {
913                     prevbin[xbin + ybin*nxbins] = (UChar_t)endxbin;
914                   }
915                 lastxvalue = value;
916               }
917             UChar_t *tempnextbin = nextbin + endxbin + 1 + ybin*nxbins;
918             memset(tempnextbin,(UChar_t)(xmax+1),xmax-endxbin);
919             if(maxvalue > 0)
920               {
921                 nextrow[ybin] = (UChar_t)ybin;
922                 if(maxvalue > lastyvalue)
923                   {
924                     UChar_t *tempnextrow = nextrow + endybin + 1;
925                     memset(tempnextrow,(UChar_t)ybin,ybin-endybin-1);
926                   }
927                 endybin = ybin;
928               }
929             lastyvalue = maxvalue;
930           }
931         UChar_t *tempnextrow = nextrow + endybin + 1;
932         memset(tempnextrow,(UChar_t)(ymax+1),ymax-endybin+1);
933       }
934       else {
935         UChar_t lastyvalue = 0;
936         Int_t endybin = ymin - 1;
937         for(Int_t ybin=ymin; ybin<=ymax; ybin++)
938           {
939             UChar_t maxvalue = 0;
940             for(Int_t xbin=xmin; xbin<=xmax; xbin++)
941               {
942                 Int_t bin = xbin + ybin*nxbins;
943                 if(gapcount[bin] < MAX_N_GAPS) {
944                   maxvalue = 1;
945                   if(trackfirstrow[bin] < nextpatchfirstrow)
946                     currentrowcount[bin] = nextpatchfirstrow;
947                   else
948                     currentrowcount[bin] = trackfirstrow[bin];
949                 }
950               }
951             if(maxvalue > 0)
952               {
953                 nextrow[ybin] = (UChar_t)ybin;
954                 if(maxvalue > lastyvalue)
955                   {
956                     UChar_t *tempnextrow = nextrow + endybin + 1;
957                     memset(tempnextrow,(UChar_t)ybin,ybin-endybin-1);
958                   }
959                 endybin = ybin;
960               }
961             lastyvalue = maxvalue;
962           }
963         UChar_t *tempnextrow = nextrow + endybin + 1;
964         memset(tempnextrow,(UChar_t)(ymax+1),ymax-endybin+1);
965       }
966     }
967
968   fBenchmark->Stop(buf);
969 }
970
971 void AliHLTHough::AddTracks()
972 {
973   // Add current slice slice tracks to the global list of found tracks
974   if(!fTracks[0])
975     {
976       cerr<<"AliHLTHough::AddTracks : No tracks"<<endl;
977       return;
978     }
979   AliHLTTrackArray *tracks = fTracks[0];
980   for(Int_t i=0; i<tracks->GetNTracks(); i++)
981     {
982       AliHLTTrack *track = tracks->GetCheckedTrack(i);
983       if(!track) continue;
984       if(track->GetNHits()!=1) cerr<<"NHITS "<<track->GetNHits()<<endl;
985       UInt_t *ids = track->GetHitNumbers();
986       ids[0] = (fCurrentSlice&0x7f)<<25;
987     }
988   
989   fGlobalTracks->AddTracks(fTracks[0],0,fCurrentSlice);
990 }
991
992 void AliHLTHough::FindTrackCandidatesRow()
993 {
994   // Find AliHLTHoughTransformerRow track candidates
995   if(fVersion != 4) {
996     LOG(AliHLTLog::kError,"AliHLTHough::FindTrackCandidatesRow()","")
997       <<"Incompatible Peak Finder version!"<<ENDLOG;
998     return;
999   }
1000
1001   //Look for peaks in histograms, and find the track candidates
1002   Int_t npatches;
1003   if(fAddHistograms)
1004     npatches = 1; //Histograms have been added.
1005   else
1006     npatches = fNPatches;
1007   
1008   Double_t initTime,cpuTime;
1009   initTime = GetCpuTime();
1010   fBenchmark->Start("Find Maxima");
1011   for(Int_t i=0; i<npatches; i++)
1012     {
1013       AliHLTHoughBaseTransformer *tr = fHoughTransformer[i];
1014       AliHLTHistogram *h = tr->GetHistogram(0);
1015       Float_t deltax = h->GetBinWidthX()*AliHLTHoughTransformerRow::GetDAlpha();
1016       Float_t deltay = h->GetBinWidthY()*AliHLTHoughTransformerRow::GetDAlpha();
1017       Float_t deltaeta = (tr->GetEtaMax()-tr->GetEtaMin())/tr->GetNEtaSegments()*AliHLTHoughTransformerRow::GetDEta();
1018       Float_t zvertex = tr->GetZVertex();
1019       fTracks[i]->Reset();
1020       fPeakFinder->Reset();
1021       
1022       for(Int_t j=0; j<fNEtaSegments; j++)
1023         {
1024           AliHLTHistogram *hist = tr->GetHistogram(j);
1025           if(hist->GetNEntries()==0) continue;
1026           fPeakFinder->SetHistogram(hist);
1027           fPeakFinder->SetEtaSlice(j);
1028           fPeakFinder->SetTrackLUTs(((AliHLTHoughTransformerRow *)tr)->GetTrackNRows(),((AliHLTHoughTransformerRow *)tr)->GetTrackFirstRow(),((AliHLTHoughTransformerRow *)tr)->GetTrackLastRow(),((AliHLTHoughTransformerRow *)tr)->GetNextRow(j));
1029 #ifdef do_mc
1030           LOG(AliHLTLog::kInformational,"AliHLTHough::FindTrackCandidates()","")
1031             <<"Starting "<<j<<" etaslice"<<ENDLOG;
1032 #endif
1033           fPeakFinder->SetThreshold(fPeakThreshold[i]);
1034           fPeakFinder->FindAdaptedRowPeaks(1,0,0);//Maxima finder for HoughTransformerRow
1035         }
1036   
1037       for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
1038         {
1039           //      if(fPeakFinder->GetWeight(k) < 0) continue;
1040           AliHLTHoughTrack *track = (AliHLTHoughTrack*)fTracks[i]->NextTrack();
1041           Double_t starteta = tr->GetEta(fPeakFinder->GetStartEta(k),fCurrentSlice);
1042           Double_t endeta = tr->GetEta(fPeakFinder->GetEndEta(k),fCurrentSlice);
1043           Double_t eta = (starteta+endeta)/2.0;
1044           track->SetTrackParametersRow(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),eta,fPeakFinder->GetWeight(k));
1045           track->SetPterr(deltax); track->SetPsierr(deltay); track->SetTglerr(deltaeta);
1046           track->SetBinXY(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetXPeakSize(k),fPeakFinder->GetYPeakSize(k));
1047           track->SetZ0(zvertex);
1048           Int_t etaindex = (fPeakFinder->GetStartEta(k)+fPeakFinder->GetEndEta(k))/2;
1049           track->SetEtaIndex(etaindex);
1050           Int_t rows[2];
1051           ((AliHLTHoughTransformerRow *)tr)->GetTrackLength(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),rows);
1052           track->SetRowRange(rows[0],rows[1]);
1053           track->SetSector(fCurrentSlice);
1054           track->SetSlice(fCurrentSlice);
1055 #ifdef do_mc
1056           Int_t label = tr->GetTrackID(etaindex,fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k));
1057           track->SetMCid(label);
1058 #endif
1059         }
1060       LOG(AliHLTLog::kInformational,"AliHLTHough::FindTrackCandidates()","")
1061         <<"Found "<<fTracks[i]->GetNTracks()<<" tracks in slice "<<fCurrentSlice<<ENDLOG;
1062       fTracks[i]->QSort();
1063     }
1064   fBenchmark->Stop("Find Maxima");
1065   cpuTime = GetCpuTime() - initTime;
1066   LOG(AliHLTLog::kInformational,"AliHLTHough::FindTrackCandidates()","Timing")
1067     <<"Maxima finding done in "<<cpuTime*1000<<" ms"<<ENDLOG;
1068 }
1069
1070 void AliHLTHough::FindTrackCandidates()
1071 {
1072   // Find AliHLTHoughTransformer track candidates
1073   if(fVersion == 4) {
1074     LOG(AliHLTLog::kError,"AliHLTHough::FindTrackCandidatesRow()","")
1075       <<"Incompatible Peak Finder version!"<<ENDLOG;
1076     return;
1077   }
1078
1079   Int_t npatches;
1080   if(fAddHistograms)
1081     npatches = 1; //Histograms have been added.
1082   else
1083     npatches = fNPatches;
1084   
1085   Double_t initTime,cpuTime;
1086   initTime = GetCpuTime();
1087   fBenchmark->Start("Find Maxima");
1088   for(Int_t i=0; i<npatches; i++)
1089     {
1090       AliHLTHoughBaseTransformer *tr = fHoughTransformer[i];
1091       fTracks[i]->Reset();
1092       
1093       for(Int_t j=0; j<fNEtaSegments; j++)
1094         {
1095           AliHLTHistogram *hist = tr->GetHistogram(j);
1096           if(hist->GetNEntries()==0) continue;
1097           fPeakFinder->Reset();
1098           fPeakFinder->SetHistogram(hist);
1099 #ifdef do_mc
1100           cout<<"Starting "<<j<<" etaslice"<<endl;
1101 #endif
1102           fPeakFinder->SetThreshold(fPeakThreshold[i]);
1103           fPeakFinder->FindAdaptedPeaks(fKappaSpread,fPeakRatio);
1104           
1105           for(Int_t k=0; k<fPeakFinder->GetEntries(); k++)
1106             {
1107               AliHLTHoughTrack *track = (AliHLTHoughTrack*)fTracks[i]->NextTrack();
1108               track->SetTrackParameters(fPeakFinder->GetXPeak(k),fPeakFinder->GetYPeak(k),fPeakFinder->GetWeight(k));
1109               track->SetEtaIndex(j);
1110               track->SetEta(tr->GetEta(j,fCurrentSlice));
1111               track->SetRowRange(AliHLTTransform::GetFirstRow(0),AliHLTTransform::GetLastRow(5));
1112             }
1113         }
1114       cout<<"Found "<<fTracks[i]->GetNTracks()<<" tracks in patch "<<i<<endl;
1115       fTracks[i]->QSort();
1116     }
1117   fBenchmark->Stop("Find Maxima");
1118   cpuTime = GetCpuTime() - initTime;
1119   LOG(AliHLTLog::kInformational,"AliHLTHough::FindTrackCandidates()","Timing")
1120     <<"Maxima finding done in "<<cpuTime*1000<<" ms"<<ENDLOG;
1121 }
1122
1123 void AliHLTHough::InitEvaluate()
1124 {
1125   //Pass the transformer objects to the AliHLTHoughEval objects:
1126   //This will provide the evaluation objects with all the necessary
1127   //data and parameters it needs.
1128   
1129   for(Int_t i=0; i<fNPatches; i++) 
1130     fEval[i]->InitTransformer(fHoughTransformer[i]);
1131 }
1132
1133 Int_t AliHLTHough::Evaluate(Int_t roadwidth,Int_t nrowstomiss)
1134 {
1135   //Evaluate the tracks, by looking along the road in the raw data.
1136   //If track does not cross all padrows - rows2miss, it is removed from the arrray.
1137   //If histograms were not added, the check is done locally in patch,
1138   //meaning that nrowstomiss is the number of padrows the road can miss with respect
1139   //to the number of rows in the patch.
1140   //If the histograms were added, the comparison is done globally in the _slice_, 
1141   //meaing that nrowstomiss is the number of padrows the road can miss with
1142   //respect to the total number of padrows in the slice.
1143   //
1144   //Return value = number of tracks which were removed (only in case of fAddHistograms)
1145   
1146   if(!fTracks[0])
1147     {
1148       LOG(AliHLTLog::kError,"AliHLTHough::Evaluate","Track Array")
1149         <<"No tracks to work with..."<<ENDLOG;
1150       return 0;
1151     }
1152   
1153   Int_t removedtracks=0;
1154   AliHLTTrackArray *tracks=0;
1155
1156   if(fAddHistograms)
1157     {
1158       tracks = fTracks[0];
1159       for(Int_t i=0; i<tracks->GetNTracks(); i++)
1160         {
1161           AliHLTTrack *track = tracks->GetCheckedTrack(i);
1162           if(!track) continue;
1163           track->SetNHits(0);
1164         }
1165     }
1166   
1167   for(Int_t i=0; i<fNPatches; i++)
1168     EvaluatePatch(i,roadwidth,nrowstomiss);
1169   
1170   //Here we check the tracks globally; 
1171   //how many good rows (padrows with signal) 
1172   //did it cross in the slice
1173   if(fAddHistograms) 
1174     {
1175       for(Int_t j=0; j<tracks->GetNTracks(); j++)
1176         {
1177           AliHLTHoughTrack *track = (AliHLTHoughTrack*)tracks->GetCheckedTrack(j);
1178           
1179           if(track->GetNHits() < AliHLTTransform::GetNRows() - nrowstomiss)
1180             {
1181               tracks->Remove(j);
1182               removedtracks++;
1183             }
1184         }
1185       tracks->Compress();
1186       tracks->QSort();
1187     }
1188     
1189   return removedtracks;
1190 }
1191
1192 void AliHLTHough::EvaluatePatch(Int_t i,Int_t roadwidth,Int_t nrowstomiss)
1193 {
1194   //Evaluate patch i.
1195   
1196   fEval[i]->InitTransformer(fHoughTransformer[i]);
1197   fEval[i]->SetNumOfPadsToLook(roadwidth);
1198   fEval[i]->SetNumOfRowsToMiss(nrowstomiss);
1199   //fEval[i]->RemoveFoundTracks();
1200   
1201   AliHLTTrackArray *tracks=0;
1202   
1203   if(!fAddHistograms)
1204     tracks = fTracks[i];
1205   else
1206     tracks = fTracks[0];
1207   
1208   Int_t nrows=0;
1209   for(Int_t j=0; j<tracks->GetNTracks(); j++)
1210     {
1211       AliHLTHoughTrack *track = (AliHLTHoughTrack*)tracks->GetCheckedTrack(j);
1212       if(!track)
1213         {
1214           LOG(AliHLTLog::kWarning,"AliHLTHough::EvaluatePatch","Track array")
1215             <<"Track object missing!"<<ENDLOG;
1216           continue;
1217         } 
1218       nrows=0;
1219       Int_t rowrange[2] = {AliHLTTransform::GetFirstRow(i),AliHLTTransform::GetLastRow(i)};
1220       Bool_t result = fEval[i]->LookInsideRoad(track,nrows,rowrange);
1221       if(fAddHistograms)
1222         {
1223           Int_t pre=track->GetNHits();
1224           track->SetNHits(pre+nrows);
1225         }
1226       else//the track crossed too few good padrows (padrows with signal) in the patch, so remove it
1227         {
1228           if(result == kFALSE)
1229             tracks->Remove(j);
1230         }
1231     }
1232   
1233   tracks->Compress();
1234
1235 }
1236
1237 void AliHLTHough::MergeEtaSlices()
1238 {
1239   //Merge tracks found in neighbouring eta slices.
1240   //Removes the track with the lower weight.
1241   
1242   fBenchmark->Start("Merge Eta-slices");
1243   AliHLTTrackArray *tracks = fTracks[0];
1244   if(!tracks)
1245     {
1246       cerr<<"AliHLTHough::MergeEtaSlices : No tracks "<<endl;
1247       return;
1248     }
1249   for(Int_t j=0; j<tracks->GetNTracks(); j++)
1250     {
1251       AliHLTHoughTrack *track1 = (AliHLTHoughTrack*)tracks->GetCheckedTrack(j);
1252       if(!track1) continue;
1253       for(Int_t k=j+1; k<tracks->GetNTracks(); k++)
1254         {
1255           AliHLTHoughTrack *track2 = (AliHLTHoughTrack*)tracks->GetCheckedTrack(k);
1256           if(!track2) continue;
1257           if(abs(track1->GetEtaIndex() - track2->GetEtaIndex()) != 1) continue;
1258           if(fabs(track1->GetKappa()-track2->GetKappa()) < 0.006 && 
1259              fabs(track1->GetPsi()- track2->GetPsi()) < 0.1)
1260             {
1261               //cout<<"Merging track in slices "<<track1->GetEtaIndex()<<" "<<track2->GetEtaIndex()<<endl;
1262               if(track1->GetWeight() > track2->GetWeight())
1263                 tracks->Remove(k);
1264               else
1265                 tracks->Remove(j);
1266             }
1267         }
1268     }
1269   fBenchmark->Stop("Merge Eta-slices");
1270   tracks->Compress();
1271 }
1272
1273 void AliHLTHough::WriteTracks(Char_t *path)
1274 {
1275   // Write found tracks into file
1276   //cout<<"AliHLTHough::WriteTracks : Sorting the tracsk"<<endl;
1277   //fGlobalTracks->QSort();
1278   
1279   Char_t filename[1024];
1280   sprintf(filename,"%s/tracks_%d.raw",path,fEvent);
1281   AliHLTMemHandler mem;
1282   mem.SetBinaryOutput(filename);
1283   mem.TrackArray2Binary(fGlobalTracks);
1284   mem.CloseBinaryOutput();
1285   fGlobalTracks->Reset();
1286 }
1287
1288 void AliHLTHough::WriteTracks(Int_t slice,Char_t *path)
1289 {
1290   // Write found tracks slice by slice into file
1291   
1292   AliHLTMemHandler mem;
1293   Char_t fname[100];
1294   if(fAddHistograms)
1295     {
1296       sprintf(fname,"%s/tracks_ho_%d_%d.raw",path,fEvent,slice);
1297       mem.SetBinaryOutput(fname);
1298       mem.TrackArray2Binary(fTracks[0]);
1299       mem.CloseBinaryOutput();
1300     }
1301   else 
1302     {
1303       for(Int_t i=0; i<fNPatches; i++)
1304         {
1305           sprintf(fname,"%s/tracks_ho_%d_%d_%d.raw",path,fEvent,slice,i);
1306           mem.SetBinaryOutput(fname);
1307           mem.TrackArray2Binary(fTracks[i]);
1308           mem.CloseBinaryOutput();
1309         }
1310     }
1311 }
1312
1313 #ifdef use_aliroot
1314 Int_t AliHLTHough::FillESD(AliESDEvent *esd)
1315 {
1316   // Fill the found hough transform tracks
1317   // into the ESD. The tracks are stored as
1318   // AliESDHLTtrack objects.
1319
1320   // No TPC PID so far,assuming pions
1321   Double_t prob[AliPID::kSPECIES];
1322   for(Int_t i=0;i<AliPID::kSPECIES;i++) {
1323     if(i==AliPID::kPion) prob[i]=1.0;
1324     else prob[i]=0.1;
1325   }
1326
1327   if(!fGlobalTracks) return 0;
1328   Int_t nglobaltracks = 0;
1329   for(Int_t i=0; i<fGlobalTracks->GetNTracks(); i++)
1330     {
1331       AliHLTHoughTrack *tpt = (AliHLTHoughTrack *)fGlobalTracks->GetCheckedTrack(i);
1332       if(!tpt) continue; 
1333       
1334       if(tpt->GetWeight()<0) continue;
1335       AliHLTHoughKalmanTrack *tpctrack = new AliHLTHoughKalmanTrack(*tpt);
1336       if(!tpctrack) continue;
1337       AliESDtrack *esdtrack2 = new AliESDtrack() ; 
1338       esdtrack2->UpdateTrackParams(tpctrack,AliESDtrack::kTPCin);
1339       esdtrack2->SetESDpid(prob);
1340       esd->AddTrack(esdtrack2);
1341       nglobaltracks++;
1342       delete esdtrack2;
1343       delete tpctrack;
1344     }
1345   return nglobaltracks;
1346 }
1347 #endif
1348
1349 void AliHLTHough::WriteDigits(Char_t *outfile)
1350 {
1351   //Write the current data to a new rootfile.
1352 #ifdef use_aliroot  
1353
1354   for(Int_t i=0; i<fNPatches; i++)
1355     {
1356       AliHLTDigitRowData *tempPt = (AliHLTDigitRowData*)fHoughTransformer[i]->GetDataPointer();
1357       fMemHandler[i]->AliDigits2RootFile(tempPt,outfile);
1358     }
1359 #else
1360   cerr<<"AliHLTHough::WriteDigits : You need to compile with AliROOT!"<<endl;
1361   return;
1362 #endif  
1363 }
1364
1365 Double_t AliHLTHough::GetCpuTime()
1366 {
1367   //Return the Cputime in seconds.
1368  struct timeval tv;
1369  gettimeofday( &tv, NULL );
1370  return tv.tv_sec+(((Double_t)tv.tv_usec)/1000000.);
1371 }
1372
1373 void *AliHLTHough::ProcessInThread(void *args)
1374 {
1375   // Called in case Hough transform tracking
1376   // is executed in a thread
1377
1378   AliHLTHough *instance = (AliHLTHough *)args;
1379   Int_t minslice = instance->GetMinSlice();
1380   Int_t maxslice = instance->GetMaxSlice();
1381   for(Int_t i=minslice; i<=maxslice; i++)
1382     {
1383       instance->ReadData(i,0);
1384       instance->Transform();
1385       instance->AddAllHistogramsRows();
1386       instance->FindTrackCandidatesRow();
1387       instance->AddTracks();
1388     }
1389   return (void *)0;
1390 }
1391
1392 void AliHLTHough::StartProcessInThread(Int_t minslice,Int_t maxslice)
1393 {
1394   // Starts the Hough transform tracking as a
1395   // separate thread. Takes as parameters the
1396   // range of TPC slices (sectors) to be reconstructed
1397
1398   if(!fThread) {
1399     char buf[255];
1400     sprintf(buf,"houghtrans_%d_%d",minslice,maxslice);
1401     SetMinMaxSlices(minslice,maxslice);
1402     //    fThread = new TThread(buf,(void (*) (void *))&ProcessInThread,(void *)this);
1403     fThread = new TThread(buf,&ProcessInThread,(void *)this);
1404     fThread->Run();
1405   }
1406   return;
1407 }
1408
1409 Int_t AliHLTHough::WaitForThreadFinish()
1410 {
1411   // Routine is used in case we run the
1412   // Hough transform tracking in several
1413   // threads and want to sync them before
1414   // writing the results to the ESD
1415
1416 #if ROOT_VERSION_CODE < 262403
1417   return TThread::Join(fThread->GetId());
1418 #else
1419   return fThread->Join(fThread->GetId());
1420 #endif
1421 }