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