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