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