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