HLT/TPCLib/tracking/AliHLTTPCHough.cxx

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