Update to the current version in the Bergen CVS. Most important
[u/mrichter/AliRoot.git] / 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