HLT/TPCLib/Ref/AliHLTTPCModeller.cxx

   1 // @(#) $Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ALICE HLT Group
   5
   6 #include "AliHLTTPCStandardIncludes.h"
   7
   8 #include "AliHLTTPCLogging.h"
   9 #include "AliHLTTPCModeller.h"
  10 #include "AliHLTTPCMemHandler.h"
  11 #include "AliHLTTPCTrackArray.h"
  12 #include "AliHLTTPCModelTrack.h"
  13 #include "AliHLTTPCDigitData.h"
  14 #include "AliHLTTPCTransform.h"
  15 #include "AliHLTTPCSpacePointData.h"
  16
  17 #ifdef use_aliroot
  18 #include "AliHLTTPCFileHandler.h"
  19 #endif
  20
  21 #if GCCVERSION == 3
  22 using namespace std;
  23 #endif
  24
  25 //_____________________________________________________________
  26 // AliHLTTPCModeller
  27 //
  28 // Class for modeling TPC data.
  29 //
  30 // This performs the cluster finding, based on track parameters.
  31 // Basically it propagates the tracks to all padrows, and looks
  32 // for a corresponding cluster. For the moment only cog is calculated,
  33 // and no deconvolution is done.
  34
  35 ClassImp(AliHLTTPCModeller)
  36
  37 AliHLTTPCModeller::AliHLTTPCModeller()
  38 {
  39   fMemHandler=0;
  40   fTracks=0;
  41   fRow=0;
  42   fTrackThreshold=0;
  43   SetOverlap();
  44   SetTrackThreshold();
  45   SetSearchRange();
  46   SetMaxClusterRange(0,0);
  47   fDebug=kFALSE;
  48 }
  49
  50
  51 AliHLTTPCModeller::~AliHLTTPCModeller()
  52 {
  53   if(fMemHandler)
  54     delete fMemHandler;
  55   if(fTracks)
  56     delete fTracks;
  57   if(fRow)
  58     delete [] fRow;
  59 }
  60
  61 void AliHLTTPCModeller::Init(Int_t slice,Int_t patch,Char_t *trackdata,Char_t *path,Bool_t houghtracks,Bool_t binary)
  62 {
  63   fSlice = slice;
  64   fPatch = patch;
  65   fHoughTracks=houghtracks;
  66
  67   sprintf(fPath,"%s",path);
  68
  69   fTracks = new AliHLTTPCTrackArray("AliHLTTPCModelTrack");
  70
  71   Char_t fname[100];
  72   AliHLTTPCMemHandler *file = new AliHLTTPCMemHandler();
  73   if(!houghtracks)
  74     sprintf(fname,"%s/tracks_tr_%d_0.raw",trackdata,fSlice); //output tracks from the tracker (no merging)
  75   else
  76     sprintf(fname,"%s/tracks_ho_%d.raw",trackdata,fSlice);
  77   //sprintf(fname,"%s/tracks_ho_%d_%d.raw",trackdata,fSlice,fPatch);
  78   if(!file->SetBinaryInput(fname))
  79     {
  80       cerr<<"AliHLTTPCModeller::Init : Error opening trackfile: "<<fname<<endl;
  81       return;
  82     }
  83   file->Binary2TrackArray(fTracks);
  84   file->CloseBinaryInput();
  85   delete file;
  86
  87   if(!houghtracks)
  88     fTracks->QSort();
  89
  90   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
  91     {
  92       AliHLTTPCModelTrack *track = (AliHLTTPCModelTrack*)fTracks->GetCheckedTrack(i);
  93       if(!track) continue;
  94       track->Init(fSlice,fPatch);
  95
  96       //Only if the tracks has been merged across sector boundaries:
  97       //if(!houghtracks)
  98       //track->Rotate(fSlice,kTRUE); //!!!!!!!!!!!!!!!!!!!
  99
 100       track->CalculateHelix();
 101     }
 102
 103   Int_t ntimes = AliHLTTPCTransform::GetNTimeBins()+1;
 104   Int_t npads = AliHLTTPCTransform::GetNPads(AliHLTTPCTransform::GetLastRow(fPatch))+1;//Max num of pads.
 105   Int_t bounds = ntimes*npads;
 106   fRow = new Digit[bounds];
 107
 108
 109   UInt_t ndigits=0;
 110   AliHLTTPCDigitRowData *digits=0;
 111 #ifdef use_aliroot
 112   fMemHandler = new AliHLTTPCFileHandler();
 113   fMemHandler->Init(slice,patch);
 114   if(binary == kFALSE)
 115     {
 116       sprintf(fname,"%s/digitfile.root",fPath);
 117       fMemHandler->SetAliInput(fname);
 118       digits = fMemHandler->AliAltroDigits2Memory(ndigits);
 119     }
 120   else
 121     {
 122       sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
 123       if(!fMemHandler->SetBinaryInput(fname))
 124         {
 125           cerr<<"AliHLTTPCModeller::Init : Error opening file "<<fname<<endl;
 126           return;
 127         }
 128       digits=(AliHLTTPCDigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
 129     }
 130 #else
 131   fMemHandler = new AliHLTTPCMemHandler();
 132   fMemHandler->Init(slice,patch);
 133   if(binary == kFALSE)
 134     {
 135       cerr<<"AliHLTTPCModeller::Init : Compile with AliROOT if you want rootfile as input"<<endl;
 136       return;
 137     }
 138   else
 139     {
 140       sprintf(fname,"%sdigits_%d_%d.raw",fPath,fSlice,fPatch);
 141       if(!fMemHandler->SetBinaryInput(fname))
 142         {
 143           cerr<<"AliHLTTPCModeller::Init : Error opening file "<<fname<<endl;
 144           return;
 145         }
 146     }
 147   digits=(AliHLTTPCDigitRowData*)fMemHandler->CompBinary2Memory(ndigits);
 148 #endif
 149
 150   SetInputData(digits);
 151 }
 152
 153 void AliHLTTPCModeller::FindClusters()
 154 {
 155   if(fDebug)
 156       {
 157 #if 0
 158     cout<<"AliHLTTPCModeller::FindClusters : Processing slice "<<fSlice<<" patch "<<fPatch<<endl;
 159 #endif
 160       }
 161   if(!fTracks)
 162     {
 163       cerr<<"AliHLTTPCModeller::Process : No tracks"<<endl;
 164       return;
 165     }
 166   if(!fRowData)
 167     {
 168       cerr<<"AliHLTTPCModeller::Process : No data "<<endl;
 169       return;
 170     }
 171
 172   AliHLTTPCDigitRowData *rowPt = fRowData;
 173   AliHLTTPCDigitData *digPt=0;
 174
 175   Int_t pad,time;
 176   Short_t charge;
 177   Cluster cluster;
 178   ClusterRegion region[200];
 179
 180   for(Int_t i=AliHLTTPCTransform::GetFirstRow(fPatch); i<=AliHLTTPCTransform::GetLastRow(fPatch); i++)
 181     {
 182       if(i != (Int_t)rowPt->fRow)
 183         {
 184           cerr<<"AliHLTTPCModeller::FindClusters : Mismatching rownumbering "<<i<<" "<<rowPt->fRow<<endl;
 185           return;
 186         }
 187       fCurrentPadRow = i;
 188       memset((void*)fRow,0,(AliHLTTPCTransform::GetNTimeBins()+1)*(AliHLTTPCTransform::GetNPads(i)+1)*sizeof(Digit));
 189       digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
 190       //cout<<"Loading row "<<i<<" with "<<(Int_t)rowPt->fNDigit<<" digits"<<endl;
 191       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 192         {
 193           pad = digPt[j].fPad;
 194           time = digPt[j].fTime;
 195           charge = digPt[j].fCharge;
 196           fRow[(AliHLTTPCTransform::GetNTimeBins()+1)*pad + time].fCharge = charge;
 197           fRow[(AliHLTTPCTransform::GetNTimeBins()+1)*pad + time].fUsed = kFALSE;
 198           //cout<<"Row "<<i<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
 199         }
 200
 201       for(Int_t k=0; k<fTracks->GetNTracks(); k++)
 202         {
 203           AliHLTTPCModelTrack *track = (AliHLTTPCModelTrack*)fTracks->GetCheckedTrack(k);
 204           if(!track) continue;
 205
 206           if(track->GetPadHit(i)<0 || track->GetTimeHit(i)<0 || track->GetNOverlaps(i)>0)//track->GetOverlap(i)>=0)
 207             {
 208               //cout<<"Track "<<k<<" is empty on row "<<i<<" "<<track->GetPadHit(i)<<" "<<track->GetTimeHit(i)<<endl;
 209               track->SetCluster(i,0,0,0,0,0,0); //The track has left the patch, or it is overlapping
 210               continue;
 211             }
 212
 213           Int_t minpad,mintime,maxpad,maxtime;
 214           minpad = mintime = 999;
 215           maxpad = maxtime = 0;
 216
 217           memset(&cluster,0,sizeof(Cluster));
 218           LocateCluster(track,region,minpad,maxpad);//,mintime,maxtime);
 219           if(maxpad - minpad + 1 > fMaxPads ||  // maxtime - mintime + 1 > fMaxTimebins ||
 220              maxpad - minpad < 1)               //  || maxtime - mintime < 1)
 221             {
 222               //cout<<"Cluster not found on row "<<i<<" maxpad "<<maxpad<<" minpad "<<minpad<<" maxtime "<<maxtime<<" mintime "<<mintime
 223               //  <<" padhit "<<track->GetPadHit(i)<<" timehit "<<track->GetTimeHit(i)<<endl;
 224
 225               track->SetCluster(i,0,0,0,0,0,0);
 226               continue;
 227             }
 228
 229           Int_t npads=0;
 230           for(pad=minpad; pad<=maxpad; pad++)
 231             {
 232               Int_t ntimes=0;
 233               for(time=region[pad].mintime; time<=region[pad].maxtime; time++)
 234                 {
 235                   charge = fRow[(AliHLTTPCTransform::GetNTimeBins()+1)*pad+time].fCharge;
 236                   if(!charge) continue;
 237                   if(fRow[(AliHLTTPCTransform::GetNTimeBins()+1)*pad+time].fUsed == kTRUE)
 238                     continue;
 239                   ntimes++;
 240
 241                   //Update the cluster parameters with this timebin
 242                   cluster.fTime += time*charge;
 243                   cluster.fPad += pad*charge;
 244                   cluster.fCharge += charge;
 245                   cluster.fSigmaY2 += pad*pad*charge;
 246                   cluster.fSigmaZ2 += time*time*charge;
 247                   fRow[(AliHLTTPCTransform::GetNTimeBins()+1)*pad+time].fUsed = kTRUE;
 248                 }
 249               if(ntimes)
 250                 npads++;
 251             }
 252           FillCluster(track,&cluster,i,npads);
 253         }
 254       FillZeros(rowPt);
 255       fMemHandler->UpdateRowPointer(rowPt);
 256     }
 257   //cout<<"done processing"<<endl;
 258
 259
 260   //Debug:
 261   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 262     {
 263       AliHLTTPCModelTrack *track = (AliHLTTPCModelTrack*)fTracks->GetCheckedTrack(i);
 264       if(!track) continue;
 265       if(track->GetNClusters() != AliHLTTPCTransform::GetNRows(fPatch))
 266         cerr<<endl<<"Mismatching hitcounts; nclusters: "<<track->GetNClusters()<<" nrows "<<AliHLTTPCTransform::GetNRows(fPatch)<<endl<<endl;
 267     }
 268
 269 }
 270
 271
 272 void AliHLTTPCModeller::LocateCluster(AliHLTTPCModelTrack *track,ClusterRegion *region,Int_t &padmin,Int_t &padmax)
 273 {
 274   //Set the cluster range
 275   //This method searches for _all_ nonzeros timebins which are neigbours.
 276   //This makes it rather impractical when dealing with high occupancy,
 277   //because then you might have very large "cluster" areas from low
 278   //pt electrons/noise.
 279
 280   Int_t row=fCurrentPadRow,charge,prtmin=0,prtmax=999;
 281   Int_t hitpad = (Int_t)rint(track->GetPadHit(row));
 282   Int_t hittime = (Int_t)rint(track->GetTimeHit(row));
 283   Int_t tmin = hittime;
 284   Int_t tmax = tmin;
 285
 286   Int_t clustercharge=0;
 287   Int_t pad=hitpad;
 288   Bool_t pm = kTRUE;
 289   Int_t npads=0,middlemax=tmax,middlemin=tmin;
 290   while(1)
 291     {
 292       Bool_t padpr=kFALSE;
 293       Int_t time = hittime;
 294       Bool_t tm = kTRUE;
 295       if(pad < 0)
 296         {
 297           padmin = 0;
 298           pad = hitpad+1;
 299           pm = kFALSE;
 300           prtmin = middlemin;
 301           prtmax = middlemax;
 302           continue;
 303         }
 304       else if(pad >= AliHLTTPCTransform::GetNPads(row))
 305         {
 306           padmax = AliHLTTPCTransform::GetNPads(row)-1;
 307           break;
 308         }
 309
 310       tmin = 999;
 311       tmax = 0;
 312       //if(row==0)
 313       //cout<<"Starting to look in pad "<<pad<<" time "<<time<<endl;
 314       while(1)
 315         {
 316           if(time < 0)
 317             {
 318               time = hittime+1;
 319               tm = kFALSE;
 320             }
 321           else if(time >= AliHLTTPCTransform::GetNTimeBins())
 322             {
 323               //timemax = AliHLTTPCTransform::GetNTimeBins()-1;
 324               break;
 325             }
 326           charge = fRow[(AliHLTTPCTransform::GetNTimeBins()+1)*pad+time].fCharge;
 327           //if(row==0)
 328           //cout<<"charge "<<charge<<" at pad "<<pad<<" time "<<time<<endl;
 329           if(charge>0)
 330             {
 331               clustercharge+=charge;
 332               padpr = kTRUE;
 333               if(time < tmin)
 334                 tmin = time;
 335               if(time > tmax)
 336                 tmax = time;
 337               if(tm)
 338                 time--;
 339               else
 340                 time++;
 341             }
 342           else
 343             {
 344               if(tm)
 345                 {
 346                   //if(abs(time - hittime) < fTimeSearch && padpr == kFALSE)//Keep looking
 347                   if(time > prtmin && npads!=0)
 348                     time--;
 349                   else
 350                     {
 351                       time = hittime+1;
 352                       tm=kFALSE;
 353                     }
 354                 }
 355               //else if(abs(time-hittime) < fTimeSearch && padpr == kFALSE)//Keep looking
 356               else if(time < prtmax && npads != 0)
 357                 time++;
 358               else
 359                 break;
 360             }
 361         }
 362       if(npads==0)
 363         {
 364           middlemax = tmax;
 365           middlemin = tmin;
 366         }
 367       //      if(row==0)
 368       //cout<<"tmax "<<tmax<<" tmin "<<tmin<<" prtmin "<<prtmin<<" ptrmax "<<prtmax<<endl;
 369
 370       if(padpr && tmax >= prtmin && tmin <= prtmax)//Sequence is overlapping with the previous
 371         {
 372           //if(row==0)
 373           //cout<<"Incrementing pad "<<endl;
 374           npads++;
 375
 376           region[pad].mintime=tmin;
 377           region[pad].maxtime=tmax;
 378
 379           /*
 380           if(tmin < timemin)
 381             timemin=tmin;
 382           if(tmax > timemax)
 383             timemax=tmax;
 384           */
 385           if(pad < padmin)
 386             padmin = pad;
 387           if(pad > padmax)
 388             padmax = pad;
 389           if(pm)
 390             pad--;
 391           else
 392             pad++;
 393
 394           prtmin = tmin;
 395           prtmax = tmax;
 396         }
 397       else
 398         {
 399           if(pm)
 400             {
 401               if(abs(pad-hitpad)<fPadSearch && clustercharge == 0)
 402                 pad--;
 403               else
 404                 {
 405                   //if(row==0)
 406                   //cout<<"Setting new pad "<<hitpad+1<<endl;
 407                   pad = hitpad+1;
 408                   pm = kFALSE;
 409                   prtmin = middlemin;
 410                   prtmax = middlemax;
 411                   continue;
 412                 }
 413             }
 414           else
 415             {
 416               if(abs(pad-hitpad)<fPadSearch && clustercharge==0)
 417                 pad++;
 418               else
 419                 break;
 420             }
 421         }
 422     }
 423
 424 }
 425
 426
 427 void AliHLTTPCModeller::FillCluster(AliHLTTPCModelTrack *track,Cluster *cluster,Int_t row,Int_t npads)
 428 {
 429   if(cluster->fCharge==0)
 430     {
 431       track->SetCluster(row,0,0,0,0,0,0);
 432       return;
 433     }
 434   Float_t fcharge = (Float_t)cluster->fCharge;
 435   Float_t fpad = ((Float_t)cluster->fPad/fcharge);
 436   Float_t ftime = ((Float_t)cluster->fTime/fcharge);
 437   Float_t sigmaY2,sigmaZ2;
 438   CalcClusterWidth(cluster,sigmaY2,sigmaZ2);
 439   track->SetCluster(row,fpad,ftime,fcharge,sigmaY2,sigmaZ2,npads);
 440 #ifdef do_mc
 441   Int_t trackID[3];
 442   GetTrackID((Int_t)rint(fpad),(Int_t)rint(ftime),trackID);
 443   track->SetClusterLabel(row,trackID);
 444 #endif
 445 }
 446
 447
 448
 449 void AliHLTTPCModeller::FillZeros(AliHLTTPCDigitRowData *rowPt,Bool_t reversesign)
 450 {
 451   //Fill zero where data has been used.
 452
 453   AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
 454   for(UInt_t j=0; j<rowPt->fNDigit; j++)
 455     {
 456       Int_t pad = digPt[j].fPad;
 457       Int_t time = digPt[j].fTime;
 458       if(fRow[(AliHLTTPCTransform::GetNTimeBins()+1)*pad+time].fUsed==kTRUE)
 459         {
 460           if(reversesign)
 461             {
 462               if(digPt[j].fCharge < 1024)
 463                 digPt[j].fCharge += 1024;
 464             }
 465           else
 466             digPt[j].fCharge = 0;
 467         }
 468     }
 469 }
 470
 471 void AliHLTTPCModeller::WriteRemaining()
 472 {
 473   //Write remaining (nonzero) digits to file.
 474
 475   AliHLTTPCDigitRowData *rowPt;
 476   rowPt = (AliHLTTPCDigitRowData*)fRowData;
 477   Int_t digitcount=0;
 478   Int_t ndigits[(AliHLTTPCTransform::GetNRows(fPatch))];
 479   for(Int_t i=AliHLTTPCTransform::GetFirstRow(fPatch); i<=AliHLTTPCTransform::GetLastRow(fPatch); i++)
 480     {
 481       AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
 482       ndigits[(i-AliHLTTPCTransform::GetFirstRow(fPatch))]=0;
 483       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 484         {
 485           if(digPt[j].fCharge==0) continue;
 486           digitcount++;
 487           ndigits[(i-AliHLTTPCTransform::GetFirstRow(fPatch))]++;
 488         }
 489       //cout<<"Difference "<<(int)ndigits[(i-AliHLTTPCTransform::GetFirstRow(fPatch))]<<" "<<(int)rowPt->fNDigit<<endl;
 490       fMemHandler->UpdateRowPointer(rowPt);
 491     }
 492
 493   Int_t size = digitcount*sizeof(AliHLTTPCDigitData) + AliHLTTPCTransform::GetNRows(fPatch)*sizeof(AliHLTTPCDigitRowData);
 494   Byte_t *data = new Byte_t[size];
 495   memset(data,0,size);
 496   AliHLTTPCDigitRowData *tempPt = (AliHLTTPCDigitRowData*)data;
 497   rowPt = (AliHLTTPCDigitRowData*)fRowData;
 498
 499   for(Int_t i=AliHLTTPCTransform::GetFirstRow(fPatch); i<=AliHLTTPCTransform::GetLastRow(fPatch); i++)
 500     {
 501       Int_t localcount=0;
 502       tempPt->fRow = i;
 503       tempPt->fNDigit = ndigits[(i-AliHLTTPCTransform::GetFirstRow(fPatch))];
 504       AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
 505       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 506         {
 507           if(digPt[j].fCharge==0) continue;
 508           if(localcount >= ndigits[(i-AliHLTTPCTransform::GetFirstRow(fPatch))])
 509             {
 510               cerr<<"AliHLTTPCModeller::WriteRemaining : Digitarray out of range!!"<<endl;
 511               return;
 512             }
 513           tempPt->fDigitData[localcount].fCharge = digPt[j].fCharge;
 514           tempPt->fDigitData[localcount].fPad = digPt[j].fPad;
 515           tempPt->fDigitData[localcount].fTime = digPt[j].fTime;
 516
 517           localcount++;
 518         }
 519       if(ndigits[(i-AliHLTTPCTransform::GetFirstRow(fPatch))] != localcount)
 520         {
 521           cerr<<"AliHLTTPCModeller::WriteRemaining : Mismatch in digitcount"<<endl;
 522           return;
 523         }
 524       fMemHandler->UpdateRowPointer(rowPt);
 525       Byte_t *tmp = (Byte_t*)tempPt;
 526       Int_t size = sizeof(AliHLTTPCDigitRowData) + ndigits[(i-AliHLTTPCTransform::GetFirstRow(fPatch))]*sizeof(AliHLTTPCDigitData);
 527       tmp += size;
 528       tempPt = (AliHLTTPCDigitRowData*)tmp;
 529     }
 530
 531   Char_t fname[100];
 532   AliHLTTPCMemHandler *mem = new AliHLTTPCMemHandler();
 533   sprintf(fname,"%s/comp/remains_%d_%d.raw",fPath,fSlice,fPatch);
 534   mem->Init(fSlice,fPatch);
 535   mem->SetBinaryOutput(fname);
 536   mem->Memory2CompBinary((UInt_t)AliHLTTPCTransform::GetNRows(fPatch),(AliHLTTPCDigitRowData*)data);
 537   mem->CloseBinaryOutput();
 538   delete mem;
 539   delete [] data;
 540 }
 541
 542 void AliHLTTPCModeller::RemoveBadTracks()
 543 {
 544   //Remove tracsk which should not be included in the compression scheme.
 545
 546   for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 547     {
 548       AliHLTTPCModelTrack *track = (AliHLTTPCModelTrack*)fTracks->GetCheckedTrack(i);
 549       if(!track) continue;
 550
 551       if(track->GetPt() < 0.08)
 552         {
 553           fTracks->Remove(i);
 554           continue;
 555         }
 556
 557       if(!fHoughTracks)
 558         if(track->GetNHits() < fTrackThreshold)
 559           fTracks->Remove(i);
 560     }
 561   fTracks->Compress();
 562
 563 }
 564
 565 void AliHLTTPCModeller::CalculateCrossingPoints()
 566 {
 567   if(fDebug)
 568       {
 569 #if 0
 570     cout<<"Calculating crossing points on "<<fTracks->GetNTracks()<<" tracks"<<endl;
 571 #endif
 572       }
 573   if(!fTracks)
 574     {
 575       cerr<<"AliHLTTPCModeller::CalculateCrossingPoints(): No tracks"<<endl;
 576       return;
 577     }
 578   Float_t hit[3];
 579
 580   Int_t sector,row;
 581   for(Int_t i=AliHLTTPCTransform::GetLastRow(fPatch); i>=AliHLTTPCTransform::GetFirstRow(fPatch); i--)
 582     {
 583       for(Int_t j=0; j<fTracks->GetNTracks(); j++)
 584         {
 585           AliHLTTPCModelTrack *track = (AliHLTTPCModelTrack*)fTracks->GetCheckedTrack(j);
 586           if(!track) continue;
 587
 588           if(!track->GetCrossingPoint(i,hit))
 589             {
 590               //cerr<<"AliHLTTPCModeller::CalculateCrossingPoints : Track "<<j<<" does not intersect row "<<i<<" :"<<endl<<
 591               //        " pt "<<track->GetPt()<<
 592               //        " tgl "<<track->GetTgl()<<" psi "<<track->GetPsi()<<" charge "<<track->GetCharge()<<endl;
 593               //fTracks->Remove(j);
 594               track->SetPadHit(i,-1);
 595               track->SetTimeHit(i,-1);
 596               continue;
 597             }
 598           //cout<<"X "<<hit[0]<<" Y "<<hit[1]<<" Z "<<hit[2]<<" tgl "<<track->GetTgl()<<endl;
 599
 600           AliHLTTPCTransform::Slice2Sector(fSlice,i,sector,row);
 601           AliHLTTPCTransform::Local2Raw(hit,sector,row);
 602           //cout<<"Pad "<<hit[1]<<" time "<<hit[2]<<" in sector "<<sector<<" row "<<row<<endl;
 603           if(hit[1]<0 || hit[1]>AliHLTTPCTransform::GetNPads(i) ||
 604              hit[2]<0 || hit[2]>AliHLTTPCTransform::GetNTimeBins())
 605             {//Track is leaving the patch, so flag the track hits (<0)
 606               track->SetPadHit(i,-1);
 607               track->SetTimeHit(i,-1);
 608               continue;
 609             }
 610
 611           track->SetPadHit(i,hit[1]);
 612           track->SetTimeHit(i,hit[2]);
 613           track->CalculateClusterWidths(i);
 614
 615           Double_t beta = track->GetCrossingAngle(i);
 616           track->SetCrossingAngleLUT(i,beta);
 617
 618           //if(hit[1]<0 || hit[2]>445)
 619           //if(hit[2]<0 || hit[2]>445)
 620           //cout<<"pad "<<hit[1]<<" time "<<hit[2]<<" pt "<<track->GetPt()<<" psi "<<track->GetPsi()<<" tgl "<<track->GetTgl()<<" firstpoint "<<track->GetFirstPointX()<<" "<<track->GetFirstPointY()<<" "<<track->GetFirstPointZ()<<endl;
 621           //cout<<"Crossing pad "<<hit[1]<<" time "<<hit[2]<<endl;
 622         }
 623     }
 624   fTracks->Compress();
 625   if(fDebug)
 626       {
 627 #if 0
 628     cout<<"And there are "<<fTracks->GetNTracks()<<" tracks remaining"<<endl;
 629 #endif
 630       }
 631 }
 632
 633 void AliHLTTPCModeller::CheckForOverlaps(Float_t dangle,Int_t *rowrange)
 634 {
 635   //Flag the tracks that overlap
 636
 637   if(fDebug)
 638       {
 639 #if 0
 640     cout<<"Checking for overlaps on "<<fTracks->GetNTracks()<<endl;
 641 #endif
 642       }
 643   Int_t counter=0;
 644
 645   for(Int_t k=AliHLTTPCTransform::GetFirstRow(fPatch); k<=AliHLTTPCTransform::GetLastRow(fPatch); k++)
 646     {
 647       if(rowrange)
 648         {
 649           if(k < rowrange[0]) continue;
 650           if(k > rowrange[1]) break;
 651         }
 652       for(Int_t i=0; i<fTracks->GetNTracks(); i++)
 653         {
 654           AliHLTTPCModelTrack *track1 = (AliHLTTPCModelTrack*)fTracks->GetCheckedTrack(i);
 655           if(!track1) continue;
 656           if(track1->GetPadHit(k)<0 || track1->GetTimeHit(k)<0) continue;
 657
 658           for(Int_t j=i+1; j<fTracks->GetNTracks(); j++)
 659             {
 660               AliHLTTPCModelTrack *track2 = (AliHLTTPCModelTrack*)fTracks->GetCheckedTrack(j);
 661               if(!track2) continue;
 662               if(track2->GetPadHit(k)<0 || track2->GetTimeHit(k)<0) continue;
 663
 664               if(abs((Int_t)rint(track1->GetPadHit(k))-(Int_t)rint(track2->GetPadHit(k))) <= fPadOverlap &&
 665                  abs((Int_t)rint(track1->GetTimeHit(k))-(Int_t)rint(track2->GetTimeHit(k))) <= fTimeOverlap)
 666                 {
 667                   if(dangle>0 && fabs(track1->GetCrossingAngleLUT(k) - track2->GetCrossingAngleLUT(k)) < dangle)
 668                     fTracks->Remove(j);
 669
 670                   //cout<<"row "<<k<<" "<<i<<" "<<j<<" "<<track1->GetPadHit(k)<<" "<<track2->GetPadHit(k)<<" "<<fabs(track1->GetCrossingAngleLUT(k) - track2->GetCrossingAngleLUT(k))<<endl;
 671
 672                   else
 673                     track1->SetOverlap(k,j);
 674                   counter++;
 675                 }
 676             }
 677         }
 678     }
 679   fTracks->Compress();
 680   if(fDebug)
 681       {
 682 #if 0
 683     cout<<"and there are "<<fTracks->GetNTracks()<<" track left"<<endl;
 684 #endif
 685       }
 686   //cout<<"found "<<counter<<" done"<<endl;
 687 }
 688
 689
 690 void AliHLTTPCModeller::CalcClusterWidth(Cluster *cl,Float_t &sigmaY2,Float_t &sigmaZ2)
 691 {
 692
 693   Float_t padw,timew;
 694
 695   padw = AliHLTTPCTransform::GetPadPitchWidth(fPatch);
 696
 697   Float_t charge = (Float_t)cl->fCharge;
 698   Float_t pad = (Float_t)cl->fPad/charge;
 699   Float_t time = (Float_t)cl->fTime/charge;
 700   Float_t s2 = (Float_t)cl->fSigmaY2/charge - pad*pad;
 701
 702   //Save the sigmas in pad and time:
 703
 704   sigmaY2 = (s2);// + 1./12);//*padw*padw;
 705
 706   /*Constants added by offline
 707     if(s2 != 0)
 708     {
 709     sigmaY2 = sigmaY2*0.108;
 710     if(fPatch<3)
 711     sigmaY2 = sigmaY2*2.07;
 712     }
 713   */
 714
 715   s2 = (Float_t)cl->fSigmaZ2/charge - time*time;
 716   timew = AliHLTTPCTransform::GetZWidth();
 717   sigmaZ2 = (s2);// +1./12);//*timew*timew;
 718
 719
 720
 721   /*
 722     Constants added by offline
 723     if(s2 != 0)
 724     {
 725     sigmaZ2 = sigmaZ2*0.169;
 726     if(fPatch < 3)
 727     sigmaZ2 = sigmaZ2*1.77;
 728     }
 729   */
 730 }
 731
 732 void AliHLTTPCModeller::GetTrackID(Int_t pad,Int_t time,Int_t *trackID)
 733 {
 734 #ifdef do_mc
 735   AliHLTTPCDigitRowData *rowPt = (AliHLTTPCDigitRowData*)fRowData;
 736
 737   trackID[0]=trackID[1]=trackID[2]=-2;
 738
 739   for(Int_t i=AliHLTTPCTransform::GetFirstRow(fPatch); i<=AliHLTTPCTransform::GetLastRow(fPatch); i++)
 740     {
 741       if(rowPt->fRow < (UInt_t)fCurrentPadRow)
 742         {
 743           AliHLTTPCMemHandler::UpdateRowPointer(rowPt);
 744           continue;
 745         }
 746       AliHLTTPCDigitData *digPt = (AliHLTTPCDigitData*)rowPt->fDigitData;
 747       for(UInt_t j=0; j<rowPt->fNDigit; j++)
 748         {
 749           Int_t cpad = digPt[j].fPad;
 750           Int_t ctime = digPt[j].fTime;
 751           if(cpad != pad) continue;
 752           if(ctime != time) continue;
 753           //if(cpad != pad && ctime != ctime) continue;
 754           //cout<<"Reading row "<<fCurrentRow<<" pad "<<cpad<<" time "<<ctime<<" trackID "<<digPt[j].fTrackID[0]<<endl;
 755           trackID[0] = digPt[j].fTrackID[0];
 756           trackID[1] = digPt[j].fTrackID[1];
 757           trackID[2] = digPt[j].fTrackID[2];
 758           break;
 759           //cout<<"Reading trackID "<<trackID[0]<<endl;
 760         }
 761       break;
 762     }
 763 #endif
 764   return;
 765 }
 766