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