L3 becomes HLT
[u/mrichter/AliRoot.git] / HLT / comp / AliHLTDataCompressor.cxx
1 // @(#) $Id$
2
3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright &copy ALICE HLT Group
5 //_____________________________________________________________
6 //
7 //  AliHLTDataCompression
8 //
9 // Interface class; binary <-> AliROOT handling of TPC data compression classes.
10 //
11
12
13 #include "AliHLTStandardIncludes.h"
14
15 #include "AliHLTLogging.h"
16 #include "AliHLTRootTypes.h"
17 #include "AliHLTTransform.h"
18 #include "AliHLTMemHandler.h"
19 #include "AliHLTSpacePointData.h"
20 #include "AliHLTCompressAC.h"
21 #include "AliHLTTrackArray.h"
22 #include "AliHLTModelTrack.h"
23 #include "AliHLTBenchmark.h"
24 #include "AliHLTClusterFitter.h"
25
26 #ifdef use_aliroot
27 #include "AliHLTFileHandler.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 "AliHLTDataCompressorHelper.h"
47 #include "AliHLTDataCompressor.h"
48 #include <math.h>
49
50 #if __GNUC__ == 3
51 using namespace std;
52 #endif
53
54
55 ClassImp(AliHLTDataCompressor)
56
57 AliHLTDataCompressor::AliHLTDataCompressor()
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(AliHLTSpacePointData*));
71 }
72
73 AliHLTDataCompressor::AliHLTDataCompressor(Char_t *path,Bool_t keep,Bool_t writeshape)
74 {
75   // constructor
76   strcpy(fPath,path);
77   fBenchmark = new AliHLTBenchmark();
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(AliHLTSpacePointData*));
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 AliHLTDataCompressor::~AliHLTDataCompressor()
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 AliHLTDataCompressor::DoBench(Char_t *fname)
113 {
114   // does benchmarking
115   fBenchmark->Analyze(fname);
116 }
117
118 void AliHLTDataCompressor::OpenOutputFile()
119 {
120   // opens the output file
121 #ifndef use_aliroot
122    LOG(AliHLTLog::kError,"AliHLTDataCompressor::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(AliHLTTransform::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 AliHLTDataCompressor::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 AliHLTDataCompressor::LoadData(Int_t event,Bool_t sp)
166 {
167   // Loads data
168   fSinglePatch=sp;
169   fEvent=event;
170   AliHLTMemHandler *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 AliHLTMemHandler();
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] = (AliHLTSpacePointData*)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   AliHLTMemHandler *tfile = new AliHLTMemHandler();
198   tfile->SetBinaryInput(fname);
199   
200   if(fInputTracks)
201     delete fInputTracks;
202   fInputTracks = new AliHLTTrackArray();
203   tfile->Binary2TrackArray(fInputTracks);
204   tfile->CloseBinaryInput();
205   delete tfile;
206 }
207
208 void AliHLTDataCompressor::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   AliHLTTrackArray *comptracks = new AliHLTTrackArray("AliHLTModelTrack");
216   fInputTracks->QSort();
217   for(Int_t i=0; i<fInputTracks->GetNTracks(); i++)
218     {
219       AliHLTTrack *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       AliHLTModelTrack *outtrack = (AliHLTModelTrack*)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           AliHLTSpacePointData *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           AliHLTTransform::Local2GlobalAngle(&angle,slice);
257           if(!intrack->CalculateReferencePoint(angle,AliHLTTransform::Row2X(padrow)))
258             {
259               cerr<<"AliHLTDataCompressor::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           AliHLTTransform::Slice2Sector(slice,padrow,sector,row);
269           AliHLTTransform::Global2Raw(xyzCross,sector,row);
270           AliHLTTransform::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 = AliHLTTransform::GetPatch(padrow);
281               Float_t sigmaY2 = points[pos].fSigmaY2 / pow(AliHLTTransform::GetPadPitchWidth(patch),2);
282               Float_t sigmaZ2 = points[pos].fSigmaZ2 / pow(AliHLTTransform::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 AliHLTModelTrack::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(AliHLTTransform::GetNRows(-1));
295     }
296   if(expand)
297     ExpandTrackData(comptracks);
298   
299   cout<<"Writing "<<comptracks->GetNTracks()<<" tracks to file"<<endl;
300   AliHLTCompress *comp = new AliHLTCompress(-1,-1,fPath,fWriteClusterShape,fEvent);
301   comp->WriteFile(comptracks);
302   delete comp;
303   delete comptracks;
304   
305 }
306
307 void AliHLTDataCompressor::ExpandTrackData(AliHLTTrackArray *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       AliHLTModelTrack *track = (AliHLTModelTrack*)tracks->GetCheckedTrack(i);
316       if(!track) continue;
317       if(track->GetNHits() == AliHLTTransform::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=AliHLTTransform::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           AliHLTSpacePointData *points = fClusters[lastSlice][0];//->GetDataPointer(size);
352           
353           Float_t angle = 0;
354           AliHLTTransform::Local2GlobalAngle(&angle,lastSlice);
355           if(!track->CalculateReferencePoint(angle,AliHLTTransform::Row2X(padrow)))
356             continue;
357           Float_t xyzCross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
358           AliHLTTransform::Global2LocHLT(xyzCross,lastSlice);
359           Float_t mindist = 123456789;
360           AliHLTSpacePointData *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               AliHLTTransform::Global2LocHLT(xyz,lastSlice);
368               
369               //Check for overflow:
370               Int_t temp = (Int_t)rint((xyzCross[1]-xyz[1])/AliHLTDataCompressorHelper::GetXYResidualStep(padrow));
371               if( abs(temp) > 1<<(AliHLTDataCompressorHelper::GetNPadBits()-1))
372                 continue;
373               
374               temp = (Int_t)rint((xyzCross[2]-xyz[2])/AliHLTDataCompressorHelper::GetZResidualStep(padrow));
375               if( abs(temp) > 1<<(AliHLTDataCompressorHelper::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               AliHLTTransform::Slice2Sector(lastSlice,padrow,sector,row);
390               AliHLTTransform::Local2Raw(xyzCross,sector,row);
391               AliHLTTransform::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 = AliHLTTransform::GetPatch(padrow);
402                   Float_t sigmaY2 = closest->fSigmaY2 / pow(AliHLTTransform::GetPadPitchWidth(patch),2);
403                   Float_t sigmaZ2 = closest->fSigmaZ2 / pow(AliHLTTransform::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 AliHLTModelTrack::FillTrack!
411               track->GetClusterModel(padrow)->fSlice=lastSlice;
412               closest->fCharge = 0;//Mark this cluster as used.
413             }
414         }
415       track->SetNClusters(AliHLTTransform::GetNRows());
416       //cout<<"Track was assigned "<<nhits<<" clusters"<<endl;
417     }
418   
419 }
420
421
422
423 void AliHLTDataCompressor::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   AliHLTCompress *comp = new AliHLTCompress(-1,-1,fPath,fWriteClusterShape,fEvent);
429   comp->ReadFile('m');
430   AliHLTTrackArray *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       AliHLTModelTrack *track = (AliHLTModelTrack*)tracks->GetCheckedTrack(i);
442       if(!track) continue;
443       for(Int_t padrow=0; padrow<AliHLTTransform::GetNRows(); padrow++)
444         {
445           if(!track->IsPresent(padrow)) continue;
446           dpad = AliHLTDataCompressorHelper::Abs(AliHLTDataCompressorHelper::Nint(track->GetClusterModel(padrow)->fDPad));
447           dtime = AliHLTDataCompressorHelper::Abs(AliHLTDataCompressorHelper::Nint(track->GetClusterModel(padrow)->fDTime));
448           charge = AliHLTDataCompressorHelper::Abs((Int_t)track->GetClusterModel(padrow)->fDCharge);
449           dsigmaY = AliHLTDataCompressorHelper::Abs(AliHLTDataCompressorHelper::Nint(track->GetClusterModel(padrow)->fDSigmaY));
450           dsigmaZ = AliHLTDataCompressorHelper::Abs(AliHLTDataCompressorHelper::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 = AliHLTDataCompressorHelper::GetNChargeBits();
471   cout<<"Updating bitnumbers; pad "<<npadbits<<" time "<<ntimebits<<" charge "<<nchargebits<<" shape "<<nshapebits<<endl;
472   AliHLTDataCompressorHelper::SetBitNumbers(npadbits,ntimebits,nchargebits,nshapebits);
473 }
474
475 void AliHLTDataCompressor::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<<"AliHLTCompressor::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       AliHLTCompress *comp = new AliHLTCompress(-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 = AliHLTTransform::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<<"AliHLTDataCompressor::WriteRemaining : Cannot open file "<<filename<<endl;
514                   exit(5);
515                 }
516
517               AliHLTSpacePointData *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               AliHLTRemainingRow *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<<"AliHLTDataCompressor::WriteRemaining : Zero row pointer "<<endl;
545                               exit(5);
546                             }
547                           if(localcounter != tempPt->fNClusters)
548                             {
549                               cerr<<"AliHLTDataCompressor::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(AliHLTRemainingRow) + npoints[padrow]*sizeof(AliHLTRemainingCluster);
559                       data = new Byte_t[size];
560                       tempPt = (AliHLTRemainingRow*)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<<"AliHLTDataCompressor::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                   AliHLTTransform::Slice2Sector(i,padrow,sector,row);
577                   AliHLTTransform::Global2Raw(xyz,sector,row);
578                   
579                   Float_t padw = points[j].fSigmaY2 / pow(AliHLTTransform::GetPadPitchWidth(AliHLTTransform::GetPatch(padrow)),2);
580                   Float_t timew = points[j].fSigmaZ2 / pow(AliHLTTransform::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 AliHLTDataCompressor::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 = AliHLTTransform::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       AliHLTSpacePointData *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) / AliHLTTransform::GetPadPitchWidth(AliHLTTransform::GetPatch(padrow));
629           Float_t timew = sqrt(points[i].fSigmaZ2) / AliHLTTransform::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           AliHLTTransform::Slice2Sector(slice,padrow,sector,row);
639           AliHLTTransform::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] >= AliHLTTransform::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 AliHLTDataCompressor::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   AliHLTCompress *comp = 0;
669   if(arithmeticCoding)
670     comp = new AliHLTCompressAC(-1,-1,fPath,fWriteClusterShape,fEvent);
671   else
672     comp = new AliHLTCompress(-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<<AliHLTDataCompressorHelper::GetNPadBits()<<' '<<AliHLTDataCompressorHelper::GetNTimeBits()<<' '
680      <<AliHLTDataCompressorHelper::GetNChargeBits()<<' '<<AliHLTDataCompressorHelper::GetNShapeBits()<<' '
681      <<AliHLTDataCompressorHelper::GetNPadBitsRemaining()<<' '<<AliHLTDataCompressorHelper::GetNTimeBitsRemaining()<<' '
682      <<AliHLTDataCompressorHelper::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 AliHLTDataCompressor::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(AliHLTLog::kError,"AliHLTDataCompressor::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(AliHLTTransform::GetParamName());
723
724   AliTPCDigitsArray *darray = new AliTPCDigitsArray();
725   darray->Setup(param);
726   darray->SetClass("AliSimDigits");
727   sprintf(filename,"TreeD_%s_%d",AliHLTTransform::GetParamName(),fEvent);
728   Bool_t ok = darray->ConnectTree(filename);
729   if(!ok)
730     {
731       cerr<<"AliHLTDataCompressor::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=AliHLTTransform::GetFirstRow(-1); padrow<=AliHLTTransform::GetLastRow(-1); padrow++)
758         {
759           Int_t sec,row;
760           AliHLTTransform::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 = AliHLTTransform::GetPatch(padrow);
766           while(counter < ncl[slice] && clPt[counter]->fPadrow == padrow)
767             {
768               Float_t temp[3];
769               AliHLTTransform::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(AliHLTTransform::GetPadPitchWidth(patch),2));
777               c->SetSigmaZ2(clPt[counter]->fSigmaZ2*pow(AliHLTTransform::GetZWidth(),2));
778               Int_t pad = AliHLTDataCompressorHelper::Nint(clPt[counter]->fPad);
779               Int_t time = AliHLTDataCompressorHelper::Nint(clPt[counter]->fTime);
780               
781               if(pad < 0)
782                 pad=0;
783               if(pad >= AliHLTTransform::GetNPads(padrow))
784                 pad = AliHLTTransform::GetNPads(padrow)-1;
785               if(time < 0 || time >= AliHLTTransform::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                       //AliHLTTransform::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 AliHLTDataCompressor::ReadUncompressedData(TempCluster **clusters,Int_t *ncl,const Int_t kmaxpoints)
837 {
838   // Reads uncompressed data  
839   AliHLTCompress *comp = new AliHLTCompress(-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   AliHLTTrackArray *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       AliHLTModelTrack *track = (AliHLTModelTrack*)tracks->GetCheckedTrack(i);
859       if(!track) continue;
860       for(Int_t padrow=0; padrow < AliHLTTransform::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 > AliHLTTransform::GetNPads(padrow) || time < -1 || time > AliHLTTransform::GetNTimeBins())
871             {
872             cerr<<"AliHLTDataCompressor::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<<"AliHLTDataCompressor::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 AliHLTDataCompressor::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       AliHLTCompress *comp = new AliHLTCompress(-1,-1,fPath,fWriteClusterShape,fEvent);
902       comp->ExpandRemaining(clusters,ncl,kmaxpoints);
903       delete comp;
904       return;
905     }
906   else
907     {
908       AliHLTMemHandler 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               AliHLTRemainingRow *tempPt = (AliHLTRemainingRow*)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(AliHLTRemainingRow),1,infile)!=1) break;
925               
926                   dPt += sizeof(AliHLTRemainingRow);
927               
928                   Int_t size = sizeof(AliHLTRemainingCluster)*tempPt->fNClusters;
929               
930                   fread(dPt,size,1,infile);
931                   dPt += size;
932                   tempPt = (AliHLTRemainingRow*)dPt;
933                   nrows++;
934                 }
935           
936               mem.CloseBinaryInput();
937               UInt_t dummy;
938               tempPt = (AliHLTRemainingRow*)mem.GetDataPointer(dummy);
939           
940               for(Int_t i=0; i<nrows; i++)
941                 {
942                   AliHLTRemainingCluster *points = tempPt->fClusters;
943                   Int_t padrow = (Int_t)tempPt->fPadRow;
944                   //Int_t sector,row;
945                   //AliHLTTransform::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] = {AliHLTTransform::Row2X(padrow),points[j].fY,points[j].fZ};
951                       //AliHLTTransform::Local2Raw(xyz,sector,row);
952                   
953                       if(ncl[slice] >= kmaxpoints)
954                         {
955                           cerr<<"AliHLTDataCompressor::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(AliHLTRemainingRow) + tempPt->fNClusters*sizeof(AliHLTRemainingCluster);
970                   dPt += size;
971                   tempPt = (AliHLTRemainingRow*)dPt;
972                 }
973           
974               mem.Free();
975             }
976         }
977     }
978 }
979
980 void AliHLTDataCompressor::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 AliHLTDataCompressor::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