HLT/hough/AliL3Hough.cxx

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