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