HLT/comp/AliL3DataCompressor.cxx

   1 // @(#) $Id$
   2
   3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
   4 //*-- Copyright &copy ALICE HLT Group
   5 //_____________________________________________________________
   6 //
   7 //  AliL3DataCompression
   8 //
   9 // Interface class; binary <-> AliROOT handling of TPC data compression classes.
  10 //
  11
  12
  13 #include "AliL3StandardIncludes.h"
  14
  15 #include "AliL3Logging.h"
  16 #include "AliL3RootTypes.h"
  17 #include "AliL3Transform.h"
  18 #include "AliL3MemHandler.h"
  19 #include "AliL3SpacePointData.h"
  20 #include "AliL3CompressAC.h"
  21 #include "AliL3TrackArray.h"
  22 #include "AliL3ModelTrack.h"
  23 #include "AliL3Benchmark.h"
  24 #include "AliL3ClusterFitter.h"
  25
  26 #ifdef use_aliroot
  27 #include "AliL3FileHandler.h"
  28 #include <AliTPCcluster.h>
  29 #include <AliTPCParamSR.h>
  30 #include <AliTPCDigitsArray.h>
  31 #include <AliTPCClustersArray.h>
  32 #include <AliTPCClustersRow.h>
  33 #include <AliSimDigits.h>
  34 #include <AliTPC.h>
  35 #include <AliTPCv2.h>
  36 #include <AliRun.h>
  37 #endif
  38
  39 #ifdef use_root
  40 #include <TFile.h>
  41 #include <TDirectory.h>
  42 #include <TSystem.h>
  43 #include <TH2F.h>
  44 #endif
  45
  46 #include "AliL3DataCompressorHelper.h"
  47 #include "AliL3DataCompressor.h"
  48 #include <math.h>
  49
  50 #if __GNUC__ == 3
  51 using namespace std;
  52 #endif
  53
  54
  55 ClassImp(AliL3DataCompressor)
  56
  57 AliL3DataCompressor::AliL3DataCompressor()
  58 {
  59   // default constructor
  60   fBenchmark=0;
  61   fInputTracks=0;
  62   fKeepRemaining=kTRUE;
  63   fNoCompression=kFALSE;
  64   fEvent=0;
  65   fWriteClusterShape=kFALSE;
  66   fOutputFile=0;
  67   fCompRatioFile=0;
  68   fNusedClusters=0;
  69   fNunusedClusters=0;
  70   memset(fClusters,0,36*6*sizeof(AliL3SpacePointData*));
  71 }
  72
  73 AliL3DataCompressor::AliL3DataCompressor(Char_t *path,Bool_t keep,Bool_t writeshape)
  74 {
  75   // constructor
  76   strcpy(fPath,path);
  77   fBenchmark = new AliL3Benchmark();
  78   fInputTracks=0;
  79   fKeepRemaining=keep;
  80   fWriteClusterShape = writeshape;
  81   fEvent=0;
  82   fOutputFile=0;
  83   fNusedClusters=0;
  84   fNunusedClusters=0;
  85   fNoCompression=kFALSE;
  86   memset(fClusters,0,36*6*sizeof(AliL3SpacePointData*));
  87 #ifdef use_root
  88   Char_t name[1024];
  89   sprintf(name,"rm -f %s/comp/*",path);//Clean the directory
  90   gSystem->Exec(name);
  91 #endif
  92   OpenOutputFile();
  93 }
  94
  95 AliL3DataCompressor::~AliL3DataCompressor()
  96 {
  97   // destructor
  98   if(fInputTracks)
  99     delete fInputTracks;
 100   if(fBenchmark)
 101     delete fBenchmark;
 102   if(fClusters)
 103     {
 104       for(Int_t i=0; i<36; i++)
 105         for(Int_t j=0; j<6; j++)
 106           if(fClusters[i][j])
 107             delete fClusters[i][j];
 108     }
 109   CloseOutputFile();
 110 }
 111
 112 void AliL3DataCompressor::DoBench(Char_t *fname)
 113 {
 114   // does benchmarking
 115   fBenchmark->Analyze(fname);
 116 }
 117
 118 void AliL3DataCompressor::OpenOutputFile()
 119 {
 120   // opens the output file
 121 #ifndef use_aliroot
 122    LOG(AliL3Log::kError,"AliL3DataCompressor::OpenOutputFile","Version")
 123      <<"You have to compile with use_aliroot flag in order to use this function"<<ENDLOG;
 124 #else
 125   Char_t filename[1024];
 126
 127   sprintf(filename,"%s/comp/comprates.txt",fPath);
 128   fCompRatioFile = new ofstream(filename);
 129
 130   if(fOutputFile)
 131     if(fOutputFile->IsOpen())
 132       fOutputFile->Close();
 133
 134   sprintf(filename,"%s/alirunfile.root",fPath);
 135   TFile *f = TFile::Open(filename);
 136   AliTPCParam *param = (AliTPCParam*)f->Get(AliL3Transform::GetParamName());
 137   sprintf(filename,"%s/comp/AliTPCclusters.root",fPath);
 138   fOutputFile = TFile::Open(filename,"RECREATE");
 139   param->Write(param->GetTitle());
 140   f->Close();
 141 #endif
 142 }
 143
 144 void AliL3DataCompressor::CloseOutputFile()
 145 {
 146   // closes the output file
 147   if(fCompRatioFile)
 148     {
 149       fCompRatioFile->close();
 150       delete fCompRatioFile;
 151     }
 152
 153   if(!fOutputFile)
 154     return;
 155 #ifdef use_root
 156   if(!fOutputFile->IsOpen())
 157     return;
 158   fOutputFile->Close();
 159 #else
 160   fclose(fOutputFile);
 161 #endif
 162   fOutputFile=0;
 163 }
 164
 165 void AliL3DataCompressor::LoadData(Int_t event,Bool_t sp)
 166 {
 167   // Loads data
 168   fSinglePatch=sp;
 169   fEvent=event;
 170   AliL3MemHandler *clusterfile[36][6];
 171   Char_t fname[1024];
 172   for(Int_t s=0; s<=35; s++)
 173     {
 174       for(Int_t p=0; p<6; p++)
 175         {
 176           if(fClusters[s][p])
 177             delete fClusters[s][p];
 178           fClusters[s][p] = 0;
 179           clusterfile[s][p] = new AliL3MemHandler();
 180           if(fSinglePatch)
 181             sprintf(fname,"%s/cf/points_%d_%d_%d.raw",fPath,fEvent,s,-1);
 182           else
 183             sprintf(fname,"%s/cf/points_%d_%d_%d.raw",fPath,fEvent,s,p);
 184           clusterfile[s][p]->SetBinaryInput(fname);
 185
 186           fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
 187           clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
 188           clusterfile[s][p]->CloseBinaryInput();
 189
 190           if(fSinglePatch)
 191             break;
 192         }
 193     }
 194
 195   //cout<<endl<<"Reading from offlinecf"<<endl<<endl;
 196   sprintf(fname,"%s/cf/tracks_%d.raw",fPath,fEvent);
 197   AliL3MemHandler *tfile = new AliL3MemHandler();
 198   tfile->SetBinaryInput(fname);
 199
 200   if(fInputTracks)
 201     delete fInputTracks;
 202   fInputTracks = new AliL3TrackArray();
 203   tfile->Binary2TrackArray(fInputTracks);
 204   tfile->CloseBinaryInput();
 205   delete tfile;
 206 }
 207
 208 void AliL3DataCompressor::FillData(Int_t minHits,Bool_t expand)
 209 {
 210
 211   //Fill the track data into track and cluster structures, and write to file.
 212   //Preparation for compressing it.
 213
 214   cout<<"Filling data; "<<fInputTracks->GetNTracks()<<" tracks"<<endl;
 215   AliL3TrackArray *comptracks = new AliL3TrackArray("AliL3ModelTrack");
 216   fInputTracks->QSort();
 217   for(Int_t i=0; i<fInputTracks->GetNTracks(); i++)
 218     {
 219       AliL3Track *intrack = fInputTracks->GetCheckedTrack(i);
 220       if(!intrack) continue;
 221
 222       if(intrack->GetNHits()<minHits) break;
 223       if(intrack->GetPt()<0.1) continue;
 224
 225       intrack->CalculateHelix();
 226
 227       AliL3ModelTrack *outtrack = (AliL3ModelTrack*)comptracks->NextTrack();
 228       outtrack->SetNHits(intrack->GetNHits());
 229       outtrack->SetRowRange(intrack->GetFirstRow(),intrack->GetLastRow());
 230       outtrack->SetFirstPoint(intrack->GetFirstPointX(),intrack->GetFirstPointY(),intrack->GetFirstPointZ());
 231       outtrack->SetLastPoint(intrack->GetLastPointX(),intrack->GetLastPointY(),intrack->GetLastPointZ());
 232       outtrack->SetPt(intrack->GetPt());
 233       outtrack->SetPsi(intrack->GetPsi());
 234       outtrack->SetTgl(intrack->GetTgl());
 235       outtrack->SetCharge(intrack->GetCharge());
 236       outtrack->CalculateHelix();
 237       Int_t nhits = intrack->GetNHits();
 238       UInt_t *hitids = intrack->GetHitNumbers();
 239       Int_t origslice = (hitids[nhits-1]>>25)&0x7f;
 240       outtrack->Init(origslice,-1);
 241
 242       for(Int_t j=nhits-1; j>=0; j--)
 243         {
 244           UInt_t id=hitids[j];
 245           Int_t slice = (id>>25)&0x7f;
 246           Int_t patch = (id>>22)&0x7;
 247           UInt_t pos = id&0x3fffff;
 248
 249           //UInt_t size;
 250           AliL3SpacePointData *points = fClusters[slice][patch];//->GetDataPointer(size);
 251           Float_t xyz[3] = {points[pos].fX,points[pos].fY,points[pos].fZ};
 252           Int_t padrow = points[pos].fPadRow;
 253
 254           //Calculate the crossing point between track and padrow
 255           Float_t angle = 0; //Perpendicular to padrow in local coordinates
 256           AliL3Transform::Local2GlobalAngle(&angle,slice);
 257           if(!intrack->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow)))
 258             {
 259               cerr<<"AliL3DataCompressor::FillData : Error in crossing point calc on slice "<<slice<<" row "<<padrow<<endl;
 260               break;
 261               //outtrack->Print(kFALSE);
 262               //exit(5);
 263             }
 264
 265           Float_t xyzCross[3] = {intrack->GetPointX(),intrack->GetPointY(),intrack->GetPointZ()};
 266
 267           Int_t sector,row;
 268           AliL3Transform::Slice2Sector(slice,padrow,sector,row);
 269           AliL3Transform::Global2Raw(xyzCross,sector,row);
 270           AliL3Transform::Global2Raw(xyz,sector,row);
 271
 272           outtrack->SetPadHit(padrow,xyzCross[1]);
 273           outtrack->SetTimeHit(padrow,xyzCross[2]);
 274
 275           outtrack->SetCrossingAngleLUT(padrow,intrack->GetCrossingAngle(padrow,slice));
 276           outtrack->CalculateClusterWidths(padrow,kTRUE);
 277
 278           if(fWriteClusterShape)
 279             {
 280               Int_t patch = AliL3Transform::GetPatch(padrow);
 281               Float_t sigmaY2 = points[pos].fSigmaY2 / pow(AliL3Transform::GetPadPitchWidth(patch),2);
 282               Float_t sigmaZ2 = points[pos].fSigmaZ2 / pow(AliL3Transform::GetZWidth(),2);
 283               outtrack->SetCluster(padrow,xyz[1],xyz[2],points[pos].fCharge,sigmaY2,sigmaZ2,3);
 284             }
 285           else
 286             outtrack->SetCluster(padrow,xyz[1],xyz[2],points[pos].fCharge,0,0,3);
 287
 288           //IMPORTANT: Set the slice in which cluster is, you need it in AliL3ModelTrack::FillTrack!
 289           outtrack->GetClusterModel(padrow)->fSlice=slice;
 290           points[pos].fCharge = 0;//Mark this cluster as used.
 291           fNusedClusters++;
 292         }
 293       if(!expand)
 294         outtrack->SetNClusters(AliL3Transform::GetNRows(-1));
 295     }
 296   if(expand)
 297     ExpandTrackData(comptracks);
 298
 299   cout<<"Writing "<<comptracks->GetNTracks()<<" tracks to file"<<endl;
 300   AliL3Compress *comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
 301   comp->WriteFile(comptracks);
 302   delete comp;
 303   delete comptracks;
 304
 305 }
 306
 307 void AliL3DataCompressor::ExpandTrackData(AliL3TrackArray *tracks)
 308 {
 309   //Loop over tracks and try to assign unused clusters.
 310   //Only clusters which are closer than the max. residual are taken.
 311
 312   cout<<"Expanding "<<tracks->GetNTracks()<<" tracks"<<endl;
 313   for(Int_t i=0; i<tracks->GetNTracks(); i++)
 314     {
 315       AliL3ModelTrack *track = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
 316       if(!track) continue;
 317       if(track->GetNHits() == AliL3Transform::GetNRows()) continue;
 318
 319       Int_t nhits = track->GetNHits();
 320       //cout<<"Expanding track with "<<nhits<<" clusters"<<endl;
 321
 322       Int_t lastSlice=-1;
 323       for(Int_t padrow=AliL3Transform::GetNRows()-1; padrow>=0; padrow--)
 324         {
 325           if(track->IsPresent(padrow))
 326             {
 327               lastSlice = track->GetClusterModel(padrow)->fSlice;
 328               continue;
 329             }
 330
 331           if(lastSlice < 0) //the outer cluster is missing, so skip it - it will be written anyhow.
 332             continue;
 333
 334           //Check the slice of the next padrow:
 335           Int_t nextPadrow = padrow-1;
 336           Int_t nextSlice = -1;
 337           while(nextPadrow >=0)
 338             {
 339               if(track->IsPresent(nextPadrow))
 340                 {
 341                   nextSlice = track->GetClusterModel(nextPadrow)->fSlice;
 342                   break;
 343                 }
 344               nextPadrow--;
 345             }
 346           if(nextSlice>=0)
 347             if(nextSlice != lastSlice)//The track crosses a slice boundary here
 348               continue;
 349
 350           //UInt_t size;
 351           AliL3SpacePointData *points = fClusters[lastSlice][0];//->GetDataPointer(size);
 352
 353           Float_t angle = 0;
 354           AliL3Transform::Local2GlobalAngle(&angle,lastSlice);
 355           if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow)))
 356             continue;
 357           Float_t xyzCross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
 358           AliL3Transform::Global2LocHLT(xyzCross,lastSlice);
 359           Float_t mindist = 123456789;
 360           AliL3SpacePointData *closest=0;
 361           for(UInt_t j=0; j<fNcl[lastSlice][0]; j++)
 362             {
 363               if(points[j].fCharge == 0) continue;// || points[j].fPadRow != padrow) continue;
 364               if(points[j].fPadRow < padrow) continue;
 365               if(points[j].fPadRow > padrow) break;
 366               Float_t xyz[3] = {points[j].fX,points[j].fY,points[j].fZ};
 367               AliL3Transform::Global2LocHLT(xyz,lastSlice);
 368
 369               //Check for overflow:
 370               Int_t temp = (Int_t)rint((xyzCross[1]-xyz[1])/AliL3DataCompressorHelper::GetXYResidualStep(padrow));
 371               if( abs(temp) > 1<<(AliL3DataCompressorHelper::GetNPadBits()-1))
 372                 continue;
 373
 374               temp = (Int_t)rint((xyzCross[2]-xyz[2])/AliL3DataCompressorHelper::GetZResidualStep(padrow));
 375               if( abs(temp) > 1<<(AliL3DataCompressorHelper::GetNTimeBits()-1))
 376                 continue;
 377
 378               Float_t dist = sqrt( pow(xyzCross[1]-xyz[1],2) + pow(xyzCross[2]-xyz[2],2) );
 379               if(dist < mindist)
 380                 {
 381                   closest = &points[j];
 382                   mindist = dist;
 383                 }
 384             }
 385           if(closest) //there was a cluster assigned
 386             {
 387               Int_t sector,row;
 388               Float_t xyz[3] = {closest->fX,closest->fY,closest->fZ};
 389               AliL3Transform::Slice2Sector(lastSlice,padrow,sector,row);
 390               AliL3Transform::Local2Raw(xyzCross,sector,row);
 391               AliL3Transform::Global2Raw(xyz,sector,row);
 392
 393               track->SetPadHit(padrow,xyzCross[1]);
 394               track->SetTimeHit(padrow,xyzCross[2]);
 395
 396               if(fWriteClusterShape)
 397                 {
 398                   Float_t angle = track->GetCrossingAngle(padrow,lastSlice);
 399                   track->SetCrossingAngleLUT(padrow,angle);
 400                   track->CalculateClusterWidths(padrow,kTRUE);
 401                   Int_t patch = AliL3Transform::GetPatch(padrow);
 402                   Float_t sigmaY2 = closest->fSigmaY2 / pow(AliL3Transform::GetPadPitchWidth(patch),2);
 403                   Float_t sigmaZ2 = closest->fSigmaZ2 / pow(AliL3Transform::GetZWidth(),2);
 404                   track->SetCluster(padrow,xyz[1],xyz[2],closest->fCharge,sigmaY2,sigmaZ2,3);
 405                 }
 406               else
 407                 track->SetCluster(padrow,xyz[1],xyz[2],closest->fCharge,0,0,3);
 408               nhits++;
 409
 410               //IMPORTANT: Set the slice in which cluster is, you need it in AliL3ModelTrack::FillTrack!
 411               track->GetClusterModel(padrow)->fSlice=lastSlice;
 412               closest->fCharge = 0;//Mark this cluster as used.
 413             }
 414         }
 415       track->SetNClusters(AliL3Transform::GetNRows());
 416       //cout<<"Track was assigned "<<nhits<<" clusters"<<endl;
 417     }
 418
 419 }
 420
 421
 422
 423 void AliL3DataCompressor::DetermineMinBits()
 424 {
 425   //Make a pass through the modelled data (after FillData has been done) to determine
 426   //how many bits is needed to encode the residuals _without_ overflows.
 427
 428   AliL3Compress *comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
 429   comp->ReadFile('m');
 430   AliL3TrackArray *tracks = comp->GetTracks();
 431   if(tracks->GetNTracks()==0)
 432     {
 433       delete comp;
 434       return;
 435     }
 436
 437   Int_t maxtime=0,maxpad=0,maxsigma=0,maxcharge=0;
 438   Int_t dpad,dtime,charge,dsigmaY,dsigmaZ,npadbits,ntimebits,nchargebits,nshapebits=0;
 439   for(Int_t i=0; i<tracks->GetNTracks(); i++)
 440     {
 441       AliL3ModelTrack *track = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
 442       if(!track) continue;
 443       for(Int_t padrow=0; padrow<AliL3Transform::GetNRows(); padrow++)
 444         {
 445           if(!track->IsPresent(padrow)) continue;
 446           dpad = AliL3DataCompressorHelper::Abs(AliL3DataCompressorHelper::Nint(track->GetClusterModel(padrow)->fDPad));
 447           dtime = AliL3DataCompressorHelper::Abs(AliL3DataCompressorHelper::Nint(track->GetClusterModel(padrow)->fDTime));
 448           charge = AliL3DataCompressorHelper::Abs((Int_t)track->GetClusterModel(padrow)->fDCharge);
 449           dsigmaY = AliL3DataCompressorHelper::Abs(AliL3DataCompressorHelper::Nint(track->GetClusterModel(padrow)->fDSigmaY));
 450           dsigmaZ = AliL3DataCompressorHelper::Abs(AliL3DataCompressorHelper::Nint(track->GetClusterModel(padrow)->fDSigmaZ));
 451           if(dpad > maxpad)
 452             maxpad=dpad;
 453           if(dtime > maxtime)
 454             maxtime=dtime;
 455           if(charge > maxcharge)
 456             maxcharge=charge;
 457           if(dsigmaY > maxsigma)
 458             maxsigma=dsigmaY;
 459           if(dsigmaZ > maxsigma)
 460             maxsigma=dsigmaZ;
 461         }
 462     }
 463   cout<<"maxpad "<<maxpad<<" maxtime "<<maxtime<<" maxcharge "<<maxcharge<<endl;
 464   npadbits = (Int_t)ceil(log(Double_t(maxpad))/log(2.)) + 1; //need 1 extra bit to encode the sign
 465   ntimebits = (Int_t)ceil(log(Double_t(maxtime))/log(2.)) + 1;
 466   nchargebits = (Int_t)ceil(log(Double_t(maxcharge))/log(2.)); //Store as a absolute value
 467   if(fWriteClusterShape)
 468     nshapebits = (Int_t)ceil(log(Double_t(maxsigma))/log(2.)) + 1;
 469
 470   nchargebits = AliL3DataCompressorHelper::GetNChargeBits();
 471   cout<<"Updating bitnumbers; pad "<<npadbits<<" time "<<ntimebits<<" charge "<<nchargebits<<" shape "<<nshapebits<<endl;
 472   AliL3DataCompressorHelper::SetBitNumbers(npadbits,ntimebits,nchargebits,nshapebits);
 473 }
 474
 475 void AliL3DataCompressor::WriteRemaining(Bool_t select)
 476 {
 477   //Write remaining clusters (not assigned to any tracks) to file
 478
 479
 480   if(!fKeepRemaining)
 481     return;
 482
 483   if(select)
 484     SelectRemainingClusters();
 485
 486   if(!fSinglePatch)
 487     {
 488       cerr<<"AliL3Compressor::WriteRemaining : You have to modify this function when not running singlepatch"<<endl;
 489       return;
 490     }
 491   if(!fNoCompression)
 492     {
 493       cout<<"Compressing remaining clusters "<<endl;
 494       AliL3Compress *comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
 495       comp->CompressRemaining(fClusters,fNcl);
 496       delete comp;
 497       return;
 498     }
 499   else
 500     {
 501       cout<<"Writing remaining clusters"<<endl;
 502       Int_t nrows = AliL3Transform::GetNRows();
 503       Int_t *npoints = new Int_t[nrows];
 504       Char_t filename[1024];
 505       for(Int_t i=0; i<=35; i++)
 506         {
 507           for(Int_t patch=0; patch < 1; patch++)
 508             {
 509               sprintf(filename,"%s/comp/remains_%d_%d_%d.raw",fPath,fEvent,i,-1);
 510               FILE *outfile = fopen(filename,"w");
 511               if(!outfile)
 512                 {
 513                   cerr<<"AliL3DataCompressor::WriteRemaining : Cannot open file "<<filename<<endl;
 514                   exit(5);
 515                 }
 516
 517               AliL3SpacePointData *points = fClusters[i][patch];
 518
 519               memset(npoints,0,nrows*sizeof(Int_t));
 520
 521               for(UInt_t j=0; j<fNcl[i][patch]; j++)
 522                 {
 523                   if(points[j].fCharge == 0) continue; //has been used
 524                   npoints[points[j].fPadRow]++;
 525                 }
 526               Int_t size =0;
 527               Byte_t *data = 0;
 528               AliL3RemainingRow *tempPt=0;
 529
 530               Int_t lastRow = -2;
 531               Int_t localcounter=0;
 532
 533               for(UInt_t j=0; j<fNcl[i][patch]; j++)
 534                 {
 535                   if(points[j].fCharge == 0) continue; //has been used
 536
 537                   Int_t padrow = points[j].fPadRow;
 538                   if(padrow != lastRow)
 539                     {
 540                       if(lastRow != -2)
 541                         {
 542                           if(!tempPt)
 543                             {
 544                               cerr<<"AliL3DataCompressor::WriteRemaining : Zero row pointer "<<endl;
 545                               exit(5);
 546                             }
 547                           if(localcounter != tempPt->fNClusters)
 548                             {
 549                               cerr<<"AliL3DataCompressor::WriteRemaining : Mismatching clustercounter "<<localcounter<<" "
 550                                   <<(Int_t)tempPt->fNClusters<<endl;
 551                               exit(5);
 552                             }
 553                           //cout<<"Writing row "<<(int)tempPt->fPadRow<<" with "<<(int)tempPt->fNClusters<<" clusters"<<endl;
 554                           fwrite(tempPt,size,1,outfile);
 555                         }
 556                       if(data)
 557                         delete [] data;
 558                       size = sizeof(AliL3RemainingRow) + npoints[padrow]*sizeof(AliL3RemainingCluster);
 559                       data = new Byte_t[size];
 560                       tempPt = (AliL3RemainingRow*)data;
 561
 562                       localcounter=0;
 563                       tempPt->fPadRow = padrow;
 564                       tempPt->fNClusters = npoints[padrow];
 565                       lastRow = padrow;
 566                     }
 567                   if(localcounter >= npoints[padrow])
 568                     {
 569                       cerr<<"AliL3DataCompressor::WriteRemaining : Cluster counter out of range: "
 570                           <<localcounter<<" "<<npoints[padrow]<<endl;
 571                       exit(5);
 572                     }
 573
 574                   Float_t xyz[3] = {points[j].fX,points[j].fY,points[j].fZ};
 575                   Int_t sector,row;
 576                   AliL3Transform::Slice2Sector(i,padrow,sector,row);
 577                   AliL3Transform::Global2Raw(xyz,sector,row);
 578
 579                   Float_t padw = points[j].fSigmaY2 / pow(AliL3Transform::GetPadPitchWidth(AliL3Transform::GetPatch(padrow)),2);
 580                   Float_t timew = points[j].fSigmaZ2 / pow(AliL3Transform::GetZWidth(),2);
 581                   tempPt->fClusters[localcounter].fPad = xyz[1];
 582                   tempPt->fClusters[localcounter].fTime = xyz[2];
 583                   tempPt->fClusters[localcounter].fCharge = points[j].fCharge;
 584                   tempPt->fClusters[localcounter].fSigmaY2 = padw;
 585                   tempPt->fClusters[localcounter].fSigmaZ2 = timew;
 586                   localcounter++;
 587                   fNunusedClusters++;
 588                 }
 589
 590               //Write the last row:
 591               fwrite(tempPt,size,1,outfile);
 592               if(data)
 593                 delete [] data;
 594               fclose(outfile);
 595             }
 596         }
 597       delete [] npoints;
 598     }
 599 }
 600
 601 void AliL3DataCompressor::SelectRemainingClusters()
 602 {
 603   //Select which remaining clusters to write in addition to the compressed data.
 604   //In particular one can here make sure that "important" clusters are not missed:
 605   //The offline track finder perform seed finding in the outer padrows;
 606   //the first seeding is using pair of points on outermost padrow and
 607   //0.125*nrows more rows towards the vertex. The second seeding uses pair
 608   //of points on the outermost padrow-0.5*0.125*nrows and 0.125*nrows + 0.5*0.125*nrows
 609   //more rows towards the vertex. In order to evaluate the seeds, the track offline
 610   //track finder checks whether a certain amount of possible clusters (padrows) is
 611   //attached to the track, and then the kalman filtering starts.
 612   //To ensure a minimal loss off efficiency, all clusters in this region should be
 613   //intact.....
 614
 615   cout<<"Cleaning up clusters"<<endl;
 616   Int_t nrows = AliL3Transform::GetNRows();
 617   Int_t gap=(Int_t)(0.125*nrows), shift=(Int_t)(0.5*gap);
 618
 619   for(Int_t slice=0; slice<36; slice++)
 620     {
 621       AliL3SpacePointData *points = fClusters[slice][0];
 622       for(UInt_t i=0; i<fNcl[slice][0]; i++)
 623         {
 624           if(points[i].fCharge == 0) continue; //Already removed
 625           Int_t padrow = (Int_t)points[i].fPadRow;
 626
 627           //Check the widths (errors) of the cluster, and remove big bastards:
 628           Float_t padw = sqrt(points[i].fSigmaY2) / AliL3Transform::GetPadPitchWidth(AliL3Transform::GetPatch(padrow));
 629           Float_t timew = sqrt(points[i].fSigmaZ2) / AliL3Transform::GetZWidth();
 630           if(padw >= 2.55 || timew >= 2.55)//Because we use 1 byte to store
 631             {
 632               points[i].fCharge = 0;
 633               continue;
 634             }
 635
 636           Float_t xyz[3] = {points[i].fX,points[i].fY,points[i].fZ};
 637           Int_t sector,row;
 638           AliL3Transform::Slice2Sector(slice,padrow,sector,row);
 639           AliL3Transform::Global2Raw(xyz,sector,row);
 640
 641           if(padrow >= nrows-1-gap-shift) continue;//save all the clusters in this region
 642
 643           //if(padrow >= nrows-1-shift) continue;
 644
 645           //Save the clusters at the borders:
 646           //if(xyz[1] < 3 || xyz[1] >= AliL3Transform::GetNPads(padrow)-4)
 647           // continue;
 648
 649           //Save clusters on padrows used for offline seeding:
 650           if(padrow == nrows - 1 || padrow == nrows - 1 - gap ||                 //First seeding
 651              padrow == nrows - 1 - shift || padrow == nrows - 1 - gap - shift)   //Second seeding
 652             continue;
 653
 654           //Cluster did not meet any of the above criteria, so disregard it:
 655           points[i].fCharge = 0;
 656         }
 657     }
 658
 659 }
 660
 661 void AliL3DataCompressor::CompressAndExpand(Bool_t arithmeticCoding)
 662 {
 663   //Read tracks/clusters from file, compress data and uncompress it. Write compression rates to file.
 664   if(fNoCompression)
 665     return;
 666
 667   cout<<"Compressing and expanding data"<<endl;
 668   AliL3Compress *comp = 0;
 669   if(arithmeticCoding)
 670     comp = new AliL3CompressAC(-1,-1,fPath,fWriteClusterShape,fEvent);
 671   else
 672     comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
 673   comp->CompressFile();
 674   comp->ExpandFile();
 675   comp->PrintCompRatio(fCompRatioFile);
 676   delete comp;
 677
 678   ofstream &out = *fCompRatioFile;
 679   out<<AliL3DataCompressorHelper::GetNPadBits()<<' '<<AliL3DataCompressorHelper::GetNTimeBits()<<' '
 680      <<AliL3DataCompressorHelper::GetNChargeBits()<<' '<<AliL3DataCompressorHelper::GetNShapeBits()<<' '
 681      <<AliL3DataCompressorHelper::GetNPadBitsRemaining()<<' '<<AliL3DataCompressorHelper::GetNTimeBitsRemaining()<<' '
 682      <<AliL3DataCompressorHelper::GetNShapeBitsRemaining()<<endl;
 683   /*
 684   //Write the ratio between used and unused clusters to comp file:
 685   out<<fNusedClusters<<' '<<fNunusedClusters<<endl;
 686   */
 687 }
 688
 689
 690 void AliL3DataCompressor::RestoreData(Bool_t remainingOnly)
 691 {
 692   //Restore the uncompressed data together with the remaining clusters,
 693   //and write to a final cluster file which serves as an input to the
 694   //final offline tracker.
 695
 696 #ifndef use_aliroot
 697    LOG(AliL3Log::kError,"AliL3DataCompressor::RestoreData","Version")
 698      <<"You have to compile with use_aliroot flag in order to use this function"<<ENDLOG;
 699 #else
 700
 701   cout<<"Restoring data"<<endl;
 702
 703   const Int_t kmaxpoints=500000;
 704   TempCluster **clusters = new TempCluster*[36];
 705   Int_t *ncl = new Int_t[36];
 706   for(Int_t i=0; i<36; i++)
 707     {
 708       ncl[i]=0;
 709       clusters[i] = new TempCluster[kmaxpoints];
 710     }
 711
 712   if(!remainingOnly)
 713     ReadUncompressedData(clusters,ncl,kmaxpoints);
 714
 715   if(fKeepRemaining)
 716     ReadRemaining(clusters,ncl,kmaxpoints);
 717
 718   Char_t filename[1024];
 719   sprintf(filename,"%s/digitfile.root",fPath);
 720   TFile *rootfile = TFile::Open(filename);
 721   rootfile->cd();
 722   AliTPCParam *param = (AliTPCParam*)rootfile->Get(AliL3Transform::GetParamName());
 723
 724   AliTPCDigitsArray *darray = new AliTPCDigitsArray();
 725   darray->Setup(param);
 726   darray->SetClass("AliSimDigits");
 727   sprintf(filename,"TreeD_%s_%d",AliL3Transform::GetParamName(),fEvent);
 728   Bool_t ok = darray->ConnectTree(filename);
 729   if(!ok)
 730     {
 731       cerr<<"AliL3DataCompressor::RestoreData : Problems connecting tree"<<endl;
 732       return;
 733     }
 734
 735   fOutputFile->cd();
 736
 737   AliTPCClustersArray *carray = new AliTPCClustersArray();
 738   carray->Setup(param);
 739   carray->SetClusterType("AliTPCcluster");
 740   carray->MakeTree();
 741
 742   Int_t totcounter=0;
 743   for(Int_t slice=0; slice<=35; slice++)
 744     {
 745       TempCluster **clPt = new TempCluster*[kmaxpoints];
 746       cout<<"Sorting "<<ncl[slice]<<" clusters in slice "<<slice<<endl;
 747       for(Int_t i=0; i<ncl[slice]; i++)
 748         clPt[i] = &clusters[slice][i];
 749
 750       if(fNusedClusters)
 751         QSort(clPt,0,ncl[slice]);
 752
 753       //cout<<"padrow "<<clPt[i]->fPadrow<<" pad "<<clPt[i]->fPad<<" time "<<clPt[i]->fTime<<endl;
 754
 755       Int_t falseid=0;
 756       Int_t counter=0;
 757       for(Int_t padrow=AliL3Transform::GetFirstRow(-1); padrow<=AliL3Transform::GetLastRow(-1); padrow++)
 758         {
 759           Int_t sec,row;
 760           AliL3Transform::Slice2Sector(slice,padrow,sec,row);
 761           AliTPCClustersRow *clrow=carray->CreateRow(sec,row);
 762           AliSimDigits *digits = (AliSimDigits*)darray->LoadRow(sec,row);
 763           digits->ExpandBuffer();
 764           digits->ExpandTrackBuffer();
 765           Int_t patch = AliL3Transform::GetPatch(padrow);
 766           while(counter < ncl[slice] && clPt[counter]->fPadrow == padrow)
 767             {
 768               Float_t temp[3];
 769               AliL3Transform::Raw2Local(temp,sec,row,clPt[counter]->fPad,clPt[counter]->fTime);
 770
 771               AliTPCcluster *c = new AliTPCcluster();
 772               c->SetY(temp[1]);
 773               c->SetZ(temp[2]);
 774               c->SetQ(clPt[counter]->fCharge);
 775
 776               c->SetSigmaY2(clPt[counter]->fSigmaY2*pow(AliL3Transform::GetPadPitchWidth(patch),2));
 777               c->SetSigmaZ2(clPt[counter]->fSigmaZ2*pow(AliL3Transform::GetZWidth(),2));
 778               Int_t pad = AliL3DataCompressorHelper::Nint(clPt[counter]->fPad);
 779               Int_t time = AliL3DataCompressorHelper::Nint(clPt[counter]->fTime);
 780
 781               if(pad < 0)
 782                 pad=0;
 783               if(pad >= AliL3Transform::GetNPads(padrow))
 784                 pad = AliL3Transform::GetNPads(padrow)-1;
 785               if(time < 0 || time >= AliL3Transform::GetNTimeBins())
 786                 cerr<<"row "<<padrow<<" pad "<<pad<<" time "<<time<<endl;
 787
 788               for(Int_t lab=0; lab<3; lab++)
 789                 {
 790                   Int_t label = digits->GetTrackIDFast(time,pad,lab);
 791                   if(label > 1)
 792                     c->SetLabel(label-2,lab);
 793                   else if(label==0)
 794                     c->SetLabel(-2,lab);
 795                   else
 796                     c->SetLabel(-1,lab);
 797                   if(lab==0 && c->GetLabel(0) < 0)
 798                     {
 799                       falseid++;
 800                       //AliL3Transform::Local2Global(temp,slice);
 801                       //cout<<"slice "<<slice<<" padrow "<<padrow<<" y "<<temp[1]<<" z "<<temp[2]<<" label "<<c->GetLabel(0)<<endl;
 802                     }
 803                 }
 804               //cout<<"row "<<padrow<<" pad "<<clPt[counter]->fPad<<" time "<<clPt[counter]->fTime<<" sigmaY2 "<<c->GetSigmaY2()<<" sigmaZ2 "<<c->GetSigmaZ2()<<endl;
 805               clrow->InsertCluster(c);
 806               delete c;
 807               counter++;
 808               totcounter++;
 809             }
 810           carray->StoreRow(sec,row);
 811           carray->ClearRow(sec,row);
 812           darray->ClearRow(sec,row);
 813         }
 814       //cerr<<"Slice "<<slice<<" nclusters "<<counter<<" falseones "<<falseid<<endl;
 815       if(counter != ncl[slice])
 816         cerr<<"AliLDataCompressor::RestoreData : Mismatching cluster count :"<<counter<<" "<<ncl[slice]<<endl;
 817       delete [] clPt;
 818     }
 819
 820   cout<<"Writing "<<totcounter<<" clusters to rootfile "<<endl;
 821
 822   sprintf(filename,"TreeC_TPC_%d",fEvent);
 823   carray->GetTree()->SetName(filename);
 824   carray->GetTree()->Write();
 825   delete carray;
 826   delete darray;
 827   rootfile->Close();
 828
 829   for(Int_t i=0; i<36; i++)
 830     delete [] clusters[i];
 831   delete [] clusters;
 832   delete [] ncl;
 833 #endif
 834 }
 835
 836 void AliL3DataCompressor::ReadUncompressedData(TempCluster **clusters,Int_t *ncl,const Int_t kmaxpoints)
 837 {
 838   // Reads uncompressed data
 839   AliL3Compress *comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
 840   if(fNoCompression)
 841     {
 842       cout<<endl<<"Reading unmodified data, no compression has been done here!!!!"<<endl<<endl;
 843       comp->ReadFile('m');//Read the unmodified data (no compression has been done).
 844     }
 845   else
 846     {
 847       cout<<"Reading uncompressed tracks "<<endl;
 848       comp->ReadFile('u');
 849     }
 850
 851   AliL3TrackArray *tracks = comp->GetTracks();
 852
 853   //Float_t totcounter=0,pcounter=0,tcounter=0;
 854   Int_t charge;
 855   Float_t pad,time,sigmaY2,sigmaZ2;
 856   for(Int_t i=0; i<tracks->GetNTracks(); i++)
 857     {
 858       AliL3ModelTrack *track = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
 859       if(!track) continue;
 860       for(Int_t padrow=0; padrow < AliL3Transform::GetNRows(-1); padrow++)
 861         {
 862           if(!track->IsPresent(padrow)) continue;
 863           track->GetPad(padrow,pad);
 864           track->GetTime(padrow,time);
 865           track->GetClusterCharge(padrow,charge);
 866           track->GetSigmaY2(padrow,sigmaY2);
 867           track->GetSigmaZ2(padrow,sigmaZ2);
 868           Int_t slice = track->GetClusterModel(padrow)->fSlice;
 869           /*
 870             if(pad < -1 || pad > AliL3Transform::GetNPads(padrow) || time < -1 || time > AliL3Transform::GetNTimeBins())
 871             {
 872             cerr<<"AliL3DataCompressor::ReadUncompressData : Wrong pad "<<pad<<" or time "<<time<<" on row "<<padrow<<" track index "<<i<<endl;
 873             track->Print();
 874             exit(5);
 875             }
 876           */
 877           if(ncl[slice] >= kmaxpoints)
 878             {
 879               cerr<<"AliL3DataCompressor::ReadUncompressedData : Too many clusters"<<endl;
 880               exit(5);
 881             }
 882           clusters[slice][ncl[slice]].fPad = pad;
 883           clusters[slice][ncl[slice]].fTime = time;
 884           clusters[slice][ncl[slice]].fCharge = charge;
 885           clusters[slice][ncl[slice]].fSigmaY2 = sigmaY2;
 886           clusters[slice][ncl[slice]].fSigmaZ2 = sigmaZ2;
 887           clusters[slice][ncl[slice]].fPadrow = padrow;
 888           //cout<<"row "<<padrow<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<" sigmas "<<sigmaY2<<" "<<sigmaZ2<<endl;
 889           ncl[slice]++;
 890         }
 891     }
 892   delete comp;
 893 }
 894
 895 void AliL3DataCompressor::ReadRemaining(TempCluster **clusters,Int_t *ncl,const Int_t kmaxpoints)
 896 {
 897   // reads remaining clusters
 898   cout<<"Reading remaining clusters "<<endl;
 899   if(!fNoCompression)
 900     {
 901       AliL3Compress *comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
 902       comp->ExpandRemaining(clusters,ncl,kmaxpoints);
 903       delete comp;
 904       return;
 905     }
 906   else
 907     {
 908       AliL3MemHandler mem;
 909       Char_t filename[1024];
 910       for(Int_t slice=0; slice<=35; slice++)
 911         {
 912           for(Int_t p=0; p<1; p++)
 913             {
 914               sprintf(filename,"%s/comp/remains_%d_%d_%d.raw",fPath,fEvent,slice,-1);
 915
 916               mem.SetBinaryInput(filename);
 917               AliL3RemainingRow *tempPt = (AliL3RemainingRow*)mem.Allocate();
 918
 919               Int_t nrows=0;
 920               FILE *infile = mem.GetFilePointer();
 921               while(!feof(infile))
 922                 {
 923                   Byte_t *dPt = (Byte_t*)tempPt;
 924                   if(fread(tempPt,sizeof(AliL3RemainingRow),1,infile)!=1) break;
 925
 926                   dPt += sizeof(AliL3RemainingRow);
 927
 928                   Int_t size = sizeof(AliL3RemainingCluster)*tempPt->fNClusters;
 929
 930                   fread(dPt,size,1,infile);
 931                   dPt += size;
 932                   tempPt = (AliL3RemainingRow*)dPt;
 933                   nrows++;
 934                 }
 935
 936               mem.CloseBinaryInput();
 937               UInt_t dummy;
 938               tempPt = (AliL3RemainingRow*)mem.GetDataPointer(dummy);
 939
 940               for(Int_t i=0; i<nrows; i++)
 941                 {
 942                   AliL3RemainingCluster *points = tempPt->fClusters;
 943                   Int_t padrow = (Int_t)tempPt->fPadRow;
 944                   //Int_t sector,row;
 945                   //AliL3Transform::Slice2Sector(slice,padrow,sector,row);
 946                   //cout<<"Loading slice "<<slice<<" row "<<padrow<<" with "<<(Int_t)tempPt->fNClusters<<" clusters "<<endl;
 947                   for(Int_t j=0; j<tempPt->fNClusters; j++)
 948                     {
 949
 950                       //Float_t xyz[3] = {AliL3Transform::Row2X(padrow),points[j].fY,points[j].fZ};
 951                       //AliL3Transform::Local2Raw(xyz,sector,row);
 952
 953                       if(ncl[slice] >= kmaxpoints)
 954                         {
 955                           cerr<<"AliL3DataCompressor::ReadRemaining : Too many clusters"<<endl;
 956                           exit(5);
 957                         }
 958                       //cout<<"slice "<<slice<<" padrow "<<padrow<<" pad "<<xyz[1]<<" time "<<xyz[2]<<endl;
 959                       clusters[slice][ncl[slice]].fPad = points[j].fPad;
 960                       clusters[slice][ncl[slice]].fTime = points[j].fTime;
 961                       clusters[slice][ncl[slice]].fCharge = points[j].fCharge;
 962                       clusters[slice][ncl[slice]].fSigmaY2 = points[j].fSigmaY2;
 963                       clusters[slice][ncl[slice]].fSigmaZ2 = points[j].fSigmaZ2;
 964                       clusters[slice][ncl[slice]].fPadrow = padrow;
 965                       //cout<<"padrow "<<padrow<<" pad "<<clusters[slice][ncl[slice]].fPad<<" time "<<clusters[slice][ncl[slice]].fTime<<" charge "<<clusters[slice][ncl[slice]].fCharge<<" widths "<<clusters[slice][ncl[slice]].fSigmaY2<<" "<<clusters[slice][ncl[slice]].fSigmaZ2<<endl;
 966                       ncl[slice]++;
 967                     }
 968                   Byte_t *dPt = (Byte_t*)tempPt;
 969                   Int_t size = sizeof(AliL3RemainingRow) + tempPt->fNClusters*sizeof(AliL3RemainingCluster);
 970                   dPt += size;
 971                   tempPt = (AliL3RemainingRow*)dPt;
 972                 }
 973
 974               mem.Free();
 975             }
 976         }
 977     }
 978 }
 979
 980 void AliL3DataCompressor::QSort(TempCluster **a, Int_t first, Int_t last)
 981 {
 982   // Implementation of quick sort
 983   static TempCluster *tmp;
 984    static int i;           // "static" to save stack space
 985    int j;
 986
 987    while (last - first > 1) {
 988       i = first;
 989       j = last;
 990       for (;;) {
 991         while (++i < last && Compare(a[i], a[first]) < 0)
 992           ;
 993         while (--j > first && Compare(a[j], a[first]) > 0)
 994           ;
 995          if (i >= j)
 996             break;
 997
 998          tmp  = a[i];
 999          a[i] = a[j];
1000          a[j] = tmp;
1001       }
1002       if (j == first) {
1003          ++first;
1004          continue;
1005       }
1006       tmp = a[first];
1007       a[first] = a[j];
1008       a[j] = tmp;
1009       if (j - first < last - (j + 1)) {
1010          QSort(a, first, j);
1011          first = j + 1;   // QSort(j + 1, last);
1012       } else {
1013          QSort(a, j + 1, last);
1014          last = j;        // QSort(first, j);
1015       }
1016    }
1017 }
1018
1019 Int_t AliL3DataCompressor::Compare(TempCluster *a,TempCluster *b)
1020 {
1021   // compares two clusters
1022   if(a->fPadrow < b->fPadrow) return -1;
1023   if(a->fPadrow > b->fPadrow) return 1;
1024
1025   if(rint(a->fPad) == rint(b->fPad) && rint(a->fTime) == rint(b->fTime)) return 0;
1026
1027   if(rint(a->fPad) < rint(b->fPad)) return -1;
1028   if(rint(a->fPad) == rint(b->fPad) && rint(a->fTime) < rint(b->fTime)) return -1;
1029
1030   return 1;
1031 }
1032