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