3 // Author: Anders Vestbo <mailto:vestbo@fi.uib.no>
4 //*-- Copyright © ALICE HLT Group
5 //_____________________________________________________________
7 // AliL3DataCompression
9 // Interface class; binary <-> AliROOT handling of TPC data compression classes.
13 #include "AliL3StandardIncludes.h"
15 #include "AliL3Logging.h"
16 #include "AliL3RootTypes.h"
17 #include "AliL3Transform.h"
18 #include "AliL3MemHandler.h"
19 #include "AliL3SpacePointData.h"
20 #include "AliL3CompressAC.h"
21 #include "AliL3TrackArray.h"
22 #include "AliL3ModelTrack.h"
23 #include "AliL3Benchmark.h"
24 #include "AliL3ClusterFitter.h"
27 #include "AliL3FileHandler.h"
28 #include <AliTPCcluster.h>
29 #include <AliTPCParamSR.h>
30 #include <AliTPCDigitsArray.h>
31 #include <AliTPCClustersArray.h>
32 #include <AliTPCClustersRow.h>
33 #include <AliSimDigits.h>
41 #include <TDirectory.h>
46 #include "AliL3DataCompressorHelper.h"
47 #include "AliL3DataCompressor.h"
55 ClassImp(AliL3DataCompressor)
57 AliL3DataCompressor::AliL3DataCompressor()
59 // default constructor
63 fNoCompression=kFALSE;
65 fWriteClusterShape=kFALSE;
70 memset(fClusters,0,36*6*sizeof(AliL3SpacePointData*));
73 AliL3DataCompressor::AliL3DataCompressor(Char_t *path,Bool_t keep,Bool_t writeshape)
77 fBenchmark = new AliL3Benchmark();
80 fWriteClusterShape = writeshape;
85 fNoCompression=kFALSE;
86 memset(fClusters,0,36*6*sizeof(AliL3SpacePointData*));
89 sprintf(name,"rm -f %s/comp/*",path);//Clean the directory
95 AliL3DataCompressor::~AliL3DataCompressor()
104 for(Int_t i=0; i<36; i++)
105 for(Int_t j=0; j<6; j++)
107 delete fClusters[i][j];
112 void AliL3DataCompressor::DoBench(Char_t *fname)
115 fBenchmark->Analyze(fname);
118 void AliL3DataCompressor::OpenOutputFile()
120 // opens the output file
122 LOG(AliL3Log::kError,"AliL3DataCompressor::OpenOutputFile","Version")
123 <<"You have to compile with use_aliroot flag in order to use this function"<<ENDLOG;
125 Char_t filename[1024];
127 sprintf(filename,"%s/comp/comprates.txt",fPath);
128 fCompRatioFile = new ofstream(filename);
131 if(fOutputFile->IsOpen())
132 fOutputFile->Close();
134 sprintf(filename,"%s/alirunfile.root",fPath);
135 TFile *f = TFile::Open(filename);
136 AliTPCParam *param = (AliTPCParam*)f->Get(AliL3Transform::GetParamName());
137 sprintf(filename,"%s/comp/AliTPCclusters.root",fPath);
138 fOutputFile = TFile::Open(filename,"RECREATE");
139 param->Write(param->GetTitle());
144 void AliL3DataCompressor::CloseOutputFile()
146 // closes the output file
149 fCompRatioFile->close();
150 delete fCompRatioFile;
156 if(!fOutputFile->IsOpen())
158 fOutputFile->Close();
165 void AliL3DataCompressor::LoadData(Int_t event,Bool_t sp)
170 AliL3MemHandler *clusterfile[36][6];
172 for(Int_t s=0; s<=35; s++)
174 for(Int_t p=0; p<6; p++)
177 delete fClusters[s][p];
179 clusterfile[s][p] = new AliL3MemHandler();
181 sprintf(fname,"%s/cf/points_%d_%d_%d.raw",fPath,fEvent,s,-1);
183 sprintf(fname,"%s/cf/points_%d_%d_%d.raw",fPath,fEvent,s,p);
184 clusterfile[s][p]->SetBinaryInput(fname);
186 fClusters[s][p] = (AliL3SpacePointData*)clusterfile[s][p]->Allocate();
187 clusterfile[s][p]->Binary2Memory(fNcl[s][p],fClusters[s][p]);
188 clusterfile[s][p]->CloseBinaryInput();
195 //cout<<endl<<"Reading from offlinecf"<<endl<<endl;
196 sprintf(fname,"%s/cf/tracks_%d.raw",fPath,fEvent);
197 AliL3MemHandler *tfile = new AliL3MemHandler();
198 tfile->SetBinaryInput(fname);
202 fInputTracks = new AliL3TrackArray();
203 tfile->Binary2TrackArray(fInputTracks);
204 tfile->CloseBinaryInput();
208 void AliL3DataCompressor::FillData(Int_t minHits,Bool_t expand)
211 //Fill the track data into track and cluster structures, and write to file.
212 //Preparation for compressing it.
214 cout<<"Filling data; "<<fInputTracks->GetNTracks()<<" tracks"<<endl;
215 AliL3TrackArray *comptracks = new AliL3TrackArray("AliL3ModelTrack");
216 fInputTracks->QSort();
217 for(Int_t i=0; i<fInputTracks->GetNTracks(); i++)
219 AliL3Track *intrack = fInputTracks->GetCheckedTrack(i);
220 if(!intrack) continue;
222 if(intrack->GetNHits()<minHits) break;
223 if(intrack->GetPt()<0.1) continue;
225 intrack->CalculateHelix();
227 AliL3ModelTrack *outtrack = (AliL3ModelTrack*)comptracks->NextTrack();
228 outtrack->SetNHits(intrack->GetNHits());
229 outtrack->SetRowRange(intrack->GetFirstRow(),intrack->GetLastRow());
230 outtrack->SetFirstPoint(intrack->GetFirstPointX(),intrack->GetFirstPointY(),intrack->GetFirstPointZ());
231 outtrack->SetLastPoint(intrack->GetLastPointX(),intrack->GetLastPointY(),intrack->GetLastPointZ());
232 outtrack->SetPt(intrack->GetPt());
233 outtrack->SetPsi(intrack->GetPsi());
234 outtrack->SetTgl(intrack->GetTgl());
235 outtrack->SetCharge(intrack->GetCharge());
236 outtrack->CalculateHelix();
237 Int_t nhits = intrack->GetNHits();
238 UInt_t *hitids = intrack->GetHitNumbers();
239 Int_t origslice = (hitids[nhits-1]>>25)&0x7f;
240 outtrack->Init(origslice,-1);
242 for(Int_t j=nhits-1; j>=0; j--)
245 Int_t slice = (id>>25)&0x7f;
246 Int_t patch = (id>>22)&0x7;
247 UInt_t pos = id&0x3fffff;
250 AliL3SpacePointData *points = fClusters[slice][patch];//->GetDataPointer(size);
251 Float_t xyz[3] = {points[pos].fX,points[pos].fY,points[pos].fZ};
252 Int_t padrow = points[pos].fPadRow;
254 //Calculate the crossing point between track and padrow
255 Float_t angle = 0; //Perpendicular to padrow in local coordinates
256 AliL3Transform::Local2GlobalAngle(&angle,slice);
257 if(!intrack->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow)))
259 cerr<<"AliL3DataCompressor::FillData : Error in crossing point calc on slice "<<slice<<" row "<<padrow<<endl;
261 //outtrack->Print(kFALSE);
265 Float_t xyzCross[3] = {intrack->GetPointX(),intrack->GetPointY(),intrack->GetPointZ()};
268 AliL3Transform::Slice2Sector(slice,padrow,sector,row);
269 AliL3Transform::Global2Raw(xyzCross,sector,row);
270 AliL3Transform::Global2Raw(xyz,sector,row);
272 outtrack->SetPadHit(padrow,xyzCross[1]);
273 outtrack->SetTimeHit(padrow,xyzCross[2]);
275 outtrack->SetCrossingAngleLUT(padrow,intrack->GetCrossingAngle(padrow,slice));
276 outtrack->CalculateClusterWidths(padrow,kTRUE);
278 if(fWriteClusterShape)
280 Int_t patch = AliL3Transform::GetPatch(padrow);
281 Float_t sigmaY2 = points[pos].fSigmaY2 / pow(AliL3Transform::GetPadPitchWidth(patch),2);
282 Float_t sigmaZ2 = points[pos].fSigmaZ2 / pow(AliL3Transform::GetZWidth(),2);
283 outtrack->SetCluster(padrow,xyz[1],xyz[2],points[pos].fCharge,sigmaY2,sigmaZ2,3);
286 outtrack->SetCluster(padrow,xyz[1],xyz[2],points[pos].fCharge,0,0,3);
288 //IMPORTANT: Set the slice in which cluster is, you need it in AliL3ModelTrack::FillTrack!
289 outtrack->GetClusterModel(padrow)->fSlice=slice;
290 points[pos].fCharge = 0;//Mark this cluster as used.
294 outtrack->SetNClusters(AliL3Transform::GetNRows(-1));
297 ExpandTrackData(comptracks);
299 cout<<"Writing "<<comptracks->GetNTracks()<<" tracks to file"<<endl;
300 AliL3Compress *comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
301 comp->WriteFile(comptracks);
307 void AliL3DataCompressor::ExpandTrackData(AliL3TrackArray *tracks)
309 //Loop over tracks and try to assign unused clusters.
310 //Only clusters which are closer than the max. residual are taken.
312 cout<<"Expanding "<<tracks->GetNTracks()<<" tracks"<<endl;
313 for(Int_t i=0; i<tracks->GetNTracks(); i++)
315 AliL3ModelTrack *track = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
317 if(track->GetNHits() == AliL3Transform::GetNRows()) continue;
319 Int_t nhits = track->GetNHits();
320 //cout<<"Expanding track with "<<nhits<<" clusters"<<endl;
323 for(Int_t padrow=AliL3Transform::GetNRows()-1; padrow>=0; padrow--)
325 if(track->IsPresent(padrow))
327 lastSlice = track->GetClusterModel(padrow)->fSlice;
331 if(lastSlice < 0) //the outer cluster is missing, so skip it - it will be written anyhow.
334 //Check the slice of the next padrow:
335 Int_t nextPadrow = padrow-1;
336 Int_t nextSlice = -1;
337 while(nextPadrow >=0)
339 if(track->IsPresent(nextPadrow))
341 nextSlice = track->GetClusterModel(nextPadrow)->fSlice;
347 if(nextSlice != lastSlice)//The track crosses a slice boundary here
351 AliL3SpacePointData *points = fClusters[lastSlice][0];//->GetDataPointer(size);
354 AliL3Transform::Local2GlobalAngle(&angle,lastSlice);
355 if(!track->CalculateReferencePoint(angle,AliL3Transform::Row2X(padrow)))
357 Float_t xyzCross[3] = {track->GetPointX(),track->GetPointY(),track->GetPointZ()};
358 AliL3Transform::Global2LocHLT(xyzCross,lastSlice);
359 Float_t mindist = 123456789;
360 AliL3SpacePointData *closest=0;
361 for(UInt_t j=0; j<fNcl[lastSlice][0]; j++)
363 if(points[j].fCharge == 0) continue;// || points[j].fPadRow != padrow) continue;
364 if(points[j].fPadRow < padrow) continue;
365 if(points[j].fPadRow > padrow) break;
366 Float_t xyz[3] = {points[j].fX,points[j].fY,points[j].fZ};
367 AliL3Transform::Global2LocHLT(xyz,lastSlice);
369 //Check for overflow:
370 Int_t temp = (Int_t)rint((xyzCross[1]-xyz[1])/AliL3DataCompressorHelper::GetXYResidualStep(padrow));
371 if( abs(temp) > 1<<(AliL3DataCompressorHelper::GetNPadBits()-1))
374 temp = (Int_t)rint((xyzCross[2]-xyz[2])/AliL3DataCompressorHelper::GetZResidualStep(padrow));
375 if( abs(temp) > 1<<(AliL3DataCompressorHelper::GetNTimeBits()-1))
378 Float_t dist = sqrt( pow(xyzCross[1]-xyz[1],2) + pow(xyzCross[2]-xyz[2],2) );
381 closest = &points[j];
385 if(closest) //there was a cluster assigned
388 Float_t xyz[3] = {closest->fX,closest->fY,closest->fZ};
389 AliL3Transform::Slice2Sector(lastSlice,padrow,sector,row);
390 AliL3Transform::Local2Raw(xyzCross,sector,row);
391 AliL3Transform::Global2Raw(xyz,sector,row);
393 track->SetPadHit(padrow,xyzCross[1]);
394 track->SetTimeHit(padrow,xyzCross[2]);
396 if(fWriteClusterShape)
398 Float_t angle = track->GetCrossingAngle(padrow,lastSlice);
399 track->SetCrossingAngleLUT(padrow,angle);
400 track->CalculateClusterWidths(padrow,kTRUE);
401 Int_t patch = AliL3Transform::GetPatch(padrow);
402 Float_t sigmaY2 = closest->fSigmaY2 / pow(AliL3Transform::GetPadPitchWidth(patch),2);
403 Float_t sigmaZ2 = closest->fSigmaZ2 / pow(AliL3Transform::GetZWidth(),2);
404 track->SetCluster(padrow,xyz[1],xyz[2],closest->fCharge,sigmaY2,sigmaZ2,3);
407 track->SetCluster(padrow,xyz[1],xyz[2],closest->fCharge,0,0,3);
410 //IMPORTANT: Set the slice in which cluster is, you need it in AliL3ModelTrack::FillTrack!
411 track->GetClusterModel(padrow)->fSlice=lastSlice;
412 closest->fCharge = 0;//Mark this cluster as used.
415 track->SetNClusters(AliL3Transform::GetNRows());
416 //cout<<"Track was assigned "<<nhits<<" clusters"<<endl;
423 void AliL3DataCompressor::DetermineMinBits()
425 //Make a pass through the modelled data (after FillData has been done) to determine
426 //how many bits is needed to encode the residuals _without_ overflows.
428 AliL3Compress *comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
430 AliL3TrackArray *tracks = comp->GetTracks();
431 if(tracks->GetNTracks()==0)
437 Int_t maxtime=0,maxpad=0,maxsigma=0,maxcharge=0;
438 Int_t dpad,dtime,charge,dsigmaY,dsigmaZ,npadbits,ntimebits,nchargebits,nshapebits=0;
439 for(Int_t i=0; i<tracks->GetNTracks(); i++)
441 AliL3ModelTrack *track = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
443 for(Int_t padrow=0; padrow<AliL3Transform::GetNRows(); padrow++)
445 if(!track->IsPresent(padrow)) continue;
446 dpad = AliL3DataCompressorHelper::Abs(AliL3DataCompressorHelper::Nint(track->GetClusterModel(padrow)->fDPad));
447 dtime = AliL3DataCompressorHelper::Abs(AliL3DataCompressorHelper::Nint(track->GetClusterModel(padrow)->fDTime));
448 charge = AliL3DataCompressorHelper::Abs((Int_t)track->GetClusterModel(padrow)->fDCharge);
449 dsigmaY = AliL3DataCompressorHelper::Abs(AliL3DataCompressorHelper::Nint(track->GetClusterModel(padrow)->fDSigmaY));
450 dsigmaZ = AliL3DataCompressorHelper::Abs(AliL3DataCompressorHelper::Nint(track->GetClusterModel(padrow)->fDSigmaZ));
455 if(charge > maxcharge)
457 if(dsigmaY > maxsigma)
459 if(dsigmaZ > maxsigma)
463 cout<<"maxpad "<<maxpad<<" maxtime "<<maxtime<<" maxcharge "<<maxcharge<<endl;
464 npadbits = (Int_t)ceil(log(Double_t(maxpad))/log(2.)) + 1; //need 1 extra bit to encode the sign
465 ntimebits = (Int_t)ceil(log(Double_t(maxtime))/log(2.)) + 1;
466 nchargebits = (Int_t)ceil(log(Double_t(maxcharge))/log(2.)); //Store as a absolute value
467 if(fWriteClusterShape)
468 nshapebits = (Int_t)ceil(log(Double_t(maxsigma))/log(2.)) + 1;
470 nchargebits = AliL3DataCompressorHelper::GetNChargeBits();
471 cout<<"Updating bitnumbers; pad "<<npadbits<<" time "<<ntimebits<<" charge "<<nchargebits<<" shape "<<nshapebits<<endl;
472 AliL3DataCompressorHelper::SetBitNumbers(npadbits,ntimebits,nchargebits,nshapebits);
475 void AliL3DataCompressor::WriteRemaining(Bool_t select)
477 //Write remaining clusters (not assigned to any tracks) to file
484 SelectRemainingClusters();
488 cerr<<"AliL3Compressor::WriteRemaining : You have to modify this function when not running singlepatch"<<endl;
493 cout<<"Compressing remaining clusters "<<endl;
494 AliL3Compress *comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
495 comp->CompressRemaining(fClusters,fNcl);
501 cout<<"Writing remaining clusters"<<endl;
502 Int_t nrows = AliL3Transform::GetNRows();
503 Int_t *npoints = new Int_t[nrows];
504 Char_t filename[1024];
505 for(Int_t i=0; i<=35; i++)
507 for(Int_t patch=0; patch < 1; patch++)
509 sprintf(filename,"%s/comp/remains_%d_%d_%d.raw",fPath,fEvent,i,-1);
510 FILE *outfile = fopen(filename,"w");
513 cerr<<"AliL3DataCompressor::WriteRemaining : Cannot open file "<<filename<<endl;
517 AliL3SpacePointData *points = fClusters[i][patch];
519 memset(npoints,0,nrows*sizeof(Int_t));
521 for(UInt_t j=0; j<fNcl[i][patch]; j++)
523 if(points[j].fCharge == 0) continue; //has been used
524 npoints[points[j].fPadRow]++;
528 AliL3RemainingRow *tempPt=0;
531 Int_t localcounter=0;
533 for(UInt_t j=0; j<fNcl[i][patch]; j++)
535 if(points[j].fCharge == 0) continue; //has been used
537 Int_t padrow = points[j].fPadRow;
538 if(padrow != lastRow)
544 cerr<<"AliL3DataCompressor::WriteRemaining : Zero row pointer "<<endl;
547 if(localcounter != tempPt->fNClusters)
549 cerr<<"AliL3DataCompressor::WriteRemaining : Mismatching clustercounter "<<localcounter<<" "
550 <<(Int_t)tempPt->fNClusters<<endl;
553 //cout<<"Writing row "<<(int)tempPt->fPadRow<<" with "<<(int)tempPt->fNClusters<<" clusters"<<endl;
554 fwrite(tempPt,size,1,outfile);
558 size = sizeof(AliL3RemainingRow) + npoints[padrow]*sizeof(AliL3RemainingCluster);
559 data = new Byte_t[size];
560 tempPt = (AliL3RemainingRow*)data;
563 tempPt->fPadRow = padrow;
564 tempPt->fNClusters = npoints[padrow];
567 if(localcounter >= npoints[padrow])
569 cerr<<"AliL3DataCompressor::WriteRemaining : Cluster counter out of range: "
570 <<localcounter<<" "<<npoints[padrow]<<endl;
574 Float_t xyz[3] = {points[j].fX,points[j].fY,points[j].fZ};
576 AliL3Transform::Slice2Sector(i,padrow,sector,row);
577 AliL3Transform::Global2Raw(xyz,sector,row);
579 Float_t padw = points[j].fSigmaY2 / pow(AliL3Transform::GetPadPitchWidth(AliL3Transform::GetPatch(padrow)),2);
580 Float_t timew = points[j].fSigmaZ2 / pow(AliL3Transform::GetZWidth(),2);
581 tempPt->fClusters[localcounter].fPad = xyz[1];
582 tempPt->fClusters[localcounter].fTime = xyz[2];
583 tempPt->fClusters[localcounter].fCharge = points[j].fCharge;
584 tempPt->fClusters[localcounter].fSigmaY2 = padw;
585 tempPt->fClusters[localcounter].fSigmaZ2 = timew;
590 //Write the last row:
591 fwrite(tempPt,size,1,outfile);
601 void AliL3DataCompressor::SelectRemainingClusters()
603 //Select which remaining clusters to write in addition to the compressed data.
604 //In particular one can here make sure that "important" clusters are not missed:
605 //The offline track finder perform seed finding in the outer padrows;
606 //the first seeding is using pair of points on outermost padrow and
607 //0.125*nrows more rows towards the vertex. The second seeding uses pair
608 //of points on the outermost padrow-0.5*0.125*nrows and 0.125*nrows + 0.5*0.125*nrows
609 //more rows towards the vertex. In order to evaluate the seeds, the track offline
610 //track finder checks whether a certain amount of possible clusters (padrows) is
611 //attached to the track, and then the kalman filtering starts.
612 //To ensure a minimal loss off efficiency, all clusters in this region should be
615 cout<<"Cleaning up clusters"<<endl;
616 Int_t nrows = AliL3Transform::GetNRows();
617 Int_t gap=(Int_t)(0.125*nrows), shift=(Int_t)(0.5*gap);
619 for(Int_t slice=0; slice<36; slice++)
621 AliL3SpacePointData *points = fClusters[slice][0];
622 for(UInt_t i=0; i<fNcl[slice][0]; i++)
624 if(points[i].fCharge == 0) continue; //Already removed
625 Int_t padrow = (Int_t)points[i].fPadRow;
627 //Check the widths (errors) of the cluster, and remove big bastards:
628 Float_t padw = sqrt(points[i].fSigmaY2) / AliL3Transform::GetPadPitchWidth(AliL3Transform::GetPatch(padrow));
629 Float_t timew = sqrt(points[i].fSigmaZ2) / AliL3Transform::GetZWidth();
630 if(padw >= 2.55 || timew >= 2.55)//Because we use 1 byte to store
632 points[i].fCharge = 0;
636 Float_t xyz[3] = {points[i].fX,points[i].fY,points[i].fZ};
638 AliL3Transform::Slice2Sector(slice,padrow,sector,row);
639 AliL3Transform::Global2Raw(xyz,sector,row);
641 if(padrow >= nrows-1-gap-shift) continue;//save all the clusters in this region
643 //if(padrow >= nrows-1-shift) continue;
645 //Save the clusters at the borders:
646 //if(xyz[1] < 3 || xyz[1] >= AliL3Transform::GetNPads(padrow)-4)
649 //Save clusters on padrows used for offline seeding:
650 if(padrow == nrows - 1 || padrow == nrows - 1 - gap || //First seeding
651 padrow == nrows - 1 - shift || padrow == nrows - 1 - gap - shift) //Second seeding
654 //Cluster did not meet any of the above criteria, so disregard it:
655 points[i].fCharge = 0;
661 void AliL3DataCompressor::CompressAndExpand(Bool_t arithmeticCoding)
663 //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 = 0;
670 comp = new AliL3CompressAC(-1,-1,fPath,fWriteClusterShape,fEvent);
672 comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
673 comp->CompressFile();
675 comp->PrintCompRatio(fCompRatioFile);
678 ofstream &out = *fCompRatioFile;
679 out<<AliL3DataCompressorHelper::GetNPadBits()<<' '<<AliL3DataCompressorHelper::GetNTimeBits()<<' '
680 <<AliL3DataCompressorHelper::GetNChargeBits()<<' '<<AliL3DataCompressorHelper::GetNShapeBits()<<' '
681 <<AliL3DataCompressorHelper::GetNPadBitsRemaining()<<' '<<AliL3DataCompressorHelper::GetNTimeBitsRemaining()<<' '
682 <<AliL3DataCompressorHelper::GetNShapeBitsRemaining()<<endl;
684 //Write the ratio between used and unused clusters to comp file:
685 out<<fNusedClusters<<' '<<fNunusedClusters<<endl;
690 void AliL3DataCompressor::RestoreData(Bool_t remainingOnly)
692 //Restore the uncompressed data together with the remaining clusters,
693 //and write to a final cluster file which serves as an input to the
694 //final offline tracker.
697 LOG(AliL3Log::kError,"AliL3DataCompressor::RestoreData","Version")
698 <<"You have to compile with use_aliroot flag in order to use this function"<<ENDLOG;
701 cout<<"Restoring data"<<endl;
703 const Int_t kmaxpoints=500000;
704 TempCluster **clusters = new TempCluster*[36];
705 Int_t *ncl = new Int_t[36];
706 for(Int_t i=0; i<36; i++)
709 clusters[i] = new TempCluster[kmaxpoints];
713 ReadUncompressedData(clusters,ncl,kmaxpoints);
716 ReadRemaining(clusters,ncl,kmaxpoints);
718 Char_t filename[1024];
719 sprintf(filename,"%s/digitfile.root",fPath);
720 TFile *rootfile = TFile::Open(filename);
722 AliTPCParam *param = (AliTPCParam*)rootfile->Get(AliL3Transform::GetParamName());
724 AliTPCDigitsArray *darray = new AliTPCDigitsArray();
725 darray->Setup(param);
726 darray->SetClass("AliSimDigits");
727 sprintf(filename,"TreeD_%s_%d",AliL3Transform::GetParamName(),fEvent);
728 Bool_t ok = darray->ConnectTree(filename);
731 cerr<<"AliL3DataCompressor::RestoreData : Problems connecting tree"<<endl;
737 AliTPCClustersArray *carray = new AliTPCClustersArray();
738 carray->Setup(param);
739 carray->SetClusterType("AliTPCcluster");
743 for(Int_t slice=0; slice<=35; slice++)
745 TempCluster **clPt = new TempCluster*[kmaxpoints];
746 cout<<"Sorting "<<ncl[slice]<<" clusters in slice "<<slice<<endl;
747 for(Int_t i=0; i<ncl[slice]; i++)
748 clPt[i] = &clusters[slice][i];
751 QSort(clPt,0,ncl[slice]);
753 //cout<<"padrow "<<clPt[i]->fPadrow<<" pad "<<clPt[i]->fPad<<" time "<<clPt[i]->fTime<<endl;
757 for(Int_t padrow=AliL3Transform::GetFirstRow(-1); padrow<=AliL3Transform::GetLastRow(-1); padrow++)
760 AliL3Transform::Slice2Sector(slice,padrow,sec,row);
761 AliTPCClustersRow *clrow=carray->CreateRow(sec,row);
762 AliSimDigits *digits = (AliSimDigits*)darray->LoadRow(sec,row);
763 digits->ExpandBuffer();
764 digits->ExpandTrackBuffer();
765 Int_t patch = AliL3Transform::GetPatch(padrow);
766 while(counter < ncl[slice] && clPt[counter]->fPadrow == padrow)
769 AliL3Transform::Raw2Local(temp,sec,row,clPt[counter]->fPad,clPt[counter]->fTime);
771 AliTPCcluster *c = new AliTPCcluster();
774 c->SetQ(clPt[counter]->fCharge);
776 c->SetSigmaY2(clPt[counter]->fSigmaY2*pow(AliL3Transform::GetPadPitchWidth(patch),2));
777 c->SetSigmaZ2(clPt[counter]->fSigmaZ2*pow(AliL3Transform::GetZWidth(),2));
778 Int_t pad = AliL3DataCompressorHelper::Nint(clPt[counter]->fPad);
779 Int_t time = AliL3DataCompressorHelper::Nint(clPt[counter]->fTime);
783 if(pad >= AliL3Transform::GetNPads(padrow))
784 pad = AliL3Transform::GetNPads(padrow)-1;
785 if(time < 0 || time >= AliL3Transform::GetNTimeBins())
786 cerr<<"row "<<padrow<<" pad "<<pad<<" time "<<time<<endl;
788 for(Int_t lab=0; lab<3; lab++)
790 Int_t label = digits->GetTrackIDFast(time,pad,lab);
792 c->SetLabel(label-2,lab);
797 if(lab==0 && c->GetLabel(0) < 0)
800 //AliL3Transform::Local2Global(temp,slice);
801 //cout<<"slice "<<slice<<" padrow "<<padrow<<" y "<<temp[1]<<" z "<<temp[2]<<" label "<<c->GetLabel(0)<<endl;
804 //cout<<"row "<<padrow<<" pad "<<clPt[counter]->fPad<<" time "<<clPt[counter]->fTime<<" sigmaY2 "<<c->GetSigmaY2()<<" sigmaZ2 "<<c->GetSigmaZ2()<<endl;
805 clrow->InsertCluster(c);
810 carray->StoreRow(sec,row);
811 carray->ClearRow(sec,row);
812 darray->ClearRow(sec,row);
814 //cerr<<"Slice "<<slice<<" nclusters "<<counter<<" falseones "<<falseid<<endl;
815 if(counter != ncl[slice])
816 cerr<<"AliLDataCompressor::RestoreData : Mismatching cluster count :"<<counter<<" "<<ncl[slice]<<endl;
820 cout<<"Writing "<<totcounter<<" clusters to rootfile "<<endl;
822 sprintf(filename,"TreeC_TPC_%d",fEvent);
823 carray->GetTree()->SetName(filename);
824 carray->GetTree()->Write();
829 for(Int_t i=0; i<36; i++)
830 delete [] clusters[i];
836 void AliL3DataCompressor::ReadUncompressedData(TempCluster **clusters,Int_t *ncl,const Int_t kmaxpoints)
838 // Reads uncompressed data
839 AliL3Compress *comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
842 cout<<endl<<"Reading unmodified data, no compression has been done here!!!!"<<endl<<endl;
843 comp->ReadFile('m');//Read the unmodified data (no compression has been done).
847 cout<<"Reading uncompressed tracks "<<endl;
851 AliL3TrackArray *tracks = comp->GetTracks();
853 //Float_t totcounter=0,pcounter=0,tcounter=0;
855 Float_t pad,time,sigmaY2,sigmaZ2;
856 for(Int_t i=0; i<tracks->GetNTracks(); i++)
858 AliL3ModelTrack *track = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
860 for(Int_t padrow=0; padrow < AliL3Transform::GetNRows(-1); padrow++)
862 if(!track->IsPresent(padrow)) continue;
863 track->GetPad(padrow,pad);
864 track->GetTime(padrow,time);
865 track->GetClusterCharge(padrow,charge);
866 track->GetSigmaY2(padrow,sigmaY2);
867 track->GetSigmaZ2(padrow,sigmaZ2);
868 Int_t slice = track->GetClusterModel(padrow)->fSlice;
870 if(pad < -1 || pad > AliL3Transform::GetNPads(padrow) || time < -1 || time > AliL3Transform::GetNTimeBins())
872 cerr<<"AliL3DataCompressor::ReadUncompressData : Wrong pad "<<pad<<" or time "<<time<<" on row "<<padrow<<" track index "<<i<<endl;
877 if(ncl[slice] >= kmaxpoints)
879 cerr<<"AliL3DataCompressor::ReadUncompressedData : Too many clusters"<<endl;
882 clusters[slice][ncl[slice]].fPad = pad;
883 clusters[slice][ncl[slice]].fTime = time;
884 clusters[slice][ncl[slice]].fCharge = charge;
885 clusters[slice][ncl[slice]].fSigmaY2 = sigmaY2;
886 clusters[slice][ncl[slice]].fSigmaZ2 = sigmaZ2;
887 clusters[slice][ncl[slice]].fPadrow = padrow;
888 //cout<<"row "<<padrow<<" pad "<<pad<<" time "<<time<<" charge "<<charge<<" sigmas "<<sigmaY2<<" "<<sigmaZ2<<endl;
895 void AliL3DataCompressor::ReadRemaining(TempCluster **clusters,Int_t *ncl,const Int_t kmaxpoints)
897 // reads remaining clusters
898 cout<<"Reading remaining clusters "<<endl;
901 AliL3Compress *comp = new AliL3Compress(-1,-1,fPath,fWriteClusterShape,fEvent);
902 comp->ExpandRemaining(clusters,ncl,kmaxpoints);
909 Char_t filename[1024];
910 for(Int_t slice=0; slice<=35; slice++)
912 for(Int_t p=0; p<1; p++)
914 sprintf(filename,"%s/comp/remains_%d_%d_%d.raw",fPath,fEvent,slice,-1);
916 mem.SetBinaryInput(filename);
917 AliL3RemainingRow *tempPt = (AliL3RemainingRow*)mem.Allocate();
920 FILE *infile = mem.GetFilePointer();
923 Byte_t *dPt = (Byte_t*)tempPt;
924 if(fread(tempPt,sizeof(AliL3RemainingRow),1,infile)!=1) break;
926 dPt += sizeof(AliL3RemainingRow);
928 Int_t size = sizeof(AliL3RemainingCluster)*tempPt->fNClusters;
930 fread(dPt,size,1,infile);
932 tempPt = (AliL3RemainingRow*)dPt;
936 mem.CloseBinaryInput();
938 tempPt = (AliL3RemainingRow*)mem.GetDataPointer(dummy);
940 for(Int_t i=0; i<nrows; i++)
942 AliL3RemainingCluster *points = tempPt->fClusters;
943 Int_t padrow = (Int_t)tempPt->fPadRow;
945 //AliL3Transform::Slice2Sector(slice,padrow,sector,row);
946 //cout<<"Loading slice "<<slice<<" row "<<padrow<<" with "<<(Int_t)tempPt->fNClusters<<" clusters "<<endl;
947 for(Int_t j=0; j<tempPt->fNClusters; j++)
950 //Float_t xyz[3] = {AliL3Transform::Row2X(padrow),points[j].fY,points[j].fZ};
951 //AliL3Transform::Local2Raw(xyz,sector,row);
953 if(ncl[slice] >= kmaxpoints)
955 cerr<<"AliL3DataCompressor::ReadRemaining : Too many clusters"<<endl;
958 //cout<<"slice "<<slice<<" padrow "<<padrow<<" pad "<<xyz[1]<<" time "<<xyz[2]<<endl;
959 clusters[slice][ncl[slice]].fPad = points[j].fPad;
960 clusters[slice][ncl[slice]].fTime = points[j].fTime;
961 clusters[slice][ncl[slice]].fCharge = points[j].fCharge;
962 clusters[slice][ncl[slice]].fSigmaY2 = points[j].fSigmaY2;
963 clusters[slice][ncl[slice]].fSigmaZ2 = points[j].fSigmaZ2;
964 clusters[slice][ncl[slice]].fPadrow = padrow;
965 //cout<<"padrow "<<padrow<<" pad "<<clusters[slice][ncl[slice]].fPad<<" time "<<clusters[slice][ncl[slice]].fTime<<" charge "<<clusters[slice][ncl[slice]].fCharge<<" widths "<<clusters[slice][ncl[slice]].fSigmaY2<<" "<<clusters[slice][ncl[slice]].fSigmaZ2<<endl;
968 Byte_t *dPt = (Byte_t*)tempPt;
969 Int_t size = sizeof(AliL3RemainingRow) + tempPt->fNClusters*sizeof(AliL3RemainingCluster);
971 tempPt = (AliL3RemainingRow*)dPt;
980 void AliL3DataCompressor::QSort(TempCluster **a, Int_t first, Int_t last)
982 // Implementation of quick sort
983 static TempCluster *tmp;
984 static int i; // "static" to save stack space
987 while (last - first > 1) {
991 while (++i < last && Compare(a[i], a[first]) < 0)
993 while (--j > first && Compare(a[j], a[first]) > 0)
1009 if (j - first < last - (j + 1)) {
1011 first = j + 1; // QSort(j + 1, last);
1013 QSort(a, j + 1, last);
1014 last = j; // QSort(first, j);
1019 Int_t AliL3DataCompressor::Compare(TempCluster *a,TempCluster *b)
1021 // compares two clusters
1022 if(a->fPadrow < b->fPadrow) return -1;
1023 if(a->fPadrow > b->fPadrow) return 1;
1025 if(rint(a->fPad) == rint(b->fPad) && rint(a->fTime) == rint(b->fTime)) return 0;
1027 if(rint(a->fPad) < rint(b->fPad)) return -1;
1028 if(rint(a->fPad) == rint(b->fPad) && rint(a->fTime) < rint(b->fTime)) return -1;