Added support for NEWIO, merged cern-hlt tree, updated to latest track candidate...
[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
40#include "AliL3DataCompressor.h"
41
42#if GCCVERSION == 3
43using namespace std;
44#endif
45
5b3f37f6 46//_____________________________________________________________
47//
48// AliL3DataCompression
49//
50// Interface class; binary <-> AliROOT handling of TPC data compression classes.
51//
52
53
54ClassImp(AliL3DataCompressor)
55
3e87ef69 56Int_t AliL3DataCompressor::fNumTimeBits = 12;
57Int_t AliL3DataCompressor::fNumPadBits = 12;
58Int_t AliL3DataCompressor::fNumChargeBits = 14;
59Int_t AliL3DataCompressor::fNumShapeBits = 14;
0a86fbb7 60Float_t AliL3DataCompressor::fXYResidualStep1 = 0.03;
61Float_t AliL3DataCompressor::fXYResidualStep2 = 0.03;
62Float_t AliL3DataCompressor::fXYResidualStep3 = 0.03;
63Float_t AliL3DataCompressor::fZResidualStep1 = 0.05;
64Float_t AliL3DataCompressor::fZResidualStep2 = 0.05;
65Float_t AliL3DataCompressor::fZResidualStep3 = 0.05;
3e87ef69 66Float_t AliL3DataCompressor::fXYWidthStep = 0.005;
67Float_t AliL3DataCompressor::fZWidthStep = 0.005;
68Int_t AliL3DataCompressor::fClusterCharge = 100;
69
5b3f37f6 70AliL3DataCompressor::AliL3DataCompressor()
71{
bd53cfb7 72 fBenchmark=0;
3e87ef69 73 fInputTracks=0;
74 fKeepRemaining=kTRUE;
75 fEvent=0;
76 fWriteClusterShape=kFALSE;
77 fOutputFile=0;
78 fCompRatioFile=0;
b2a02bce 79 fNusedClusters=0;
80 fNunusedClusters=0;
3e87ef69 81 memset(fClusters,0,36*6*sizeof(AliL3SpacePointData*));
5b3f37f6 82}
83
3e87ef69 84AliL3DataCompressor::AliL3DataCompressor(Char_t *path,Bool_t keep,Bool_t writeshape)
5b3f37f6 85{
5b3f37f6 86 strcpy(fPath,path);
bd53cfb7 87 fBenchmark = new AliL3Benchmark();
3e87ef69 88 fInputTracks=0;
89 fKeepRemaining=keep;
90 fWriteClusterShape = writeshape;
91 fEvent=0;
92 fOutputFile=0;
b2a02bce 93 fNusedClusters=0;
94 fNunusedClusters=0;
3e87ef69 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();
5b3f37f6 102}
103
104AliL3DataCompressor::~AliL3DataCompressor()
105{
3e87ef69 106 if(fInputTracks)
107 delete fInputTracks;
bd53cfb7 108 if(fBenchmark)
109 delete fBenchmark;
3e87ef69 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();
bd53cfb7 118}
119
120void AliL3DataCompressor::DoBench(Char_t *fname)
121{
122 fBenchmark->Analyze(fname);
5b3f37f6 123}
124
3e87ef69 125void 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
0a86fbb7 133void 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
141void 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
149const Float_t AliL3DataCompressor::GetXYResidualStep(Int_t row)
3e87ef69 150{
0a86fbb7 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
164const 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 }
3e87ef69 177}
178
179void AliL3DataCompressor::OpenOutputFile()
5b3f37f6 180{
3e87ef69 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];
5b3f37f6 186
3e87ef69 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
204void AliL3DataCompressor::CloseOutputFile()
205{
206 if(fCompRatioFile)
5b3f37f6 207 {
3e87ef69 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}
5b3f37f6 223
3e87ef69 224void 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);
5b3f37f6 243
3e87ef69 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();
5b3f37f6 247
3e87ef69 248 if(fSinglePatch)
249 break;
5b3f37f6 250 }
251 }
252
3e87ef69 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;
5b3f37f6 263}
264
3e87ef69 265void AliL3DataCompressor::FillData(Int_t min_hits,Bool_t expand)
5b3f37f6 266{
5b3f37f6 267
3e87ef69 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++)
5b3f37f6 275 {
3e87ef69 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--)
5b3f37f6 298 {
3e87ef69 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;
5b3f37f6 308
3e87ef69 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;
0a86fbb7 315 break;
316 //outtrack->Print(kFALSE);
317 //exit(5);
3e87ef69 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);
5b3f37f6 342
3e87ef69 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.
b2a02bce 346 fNusedClusters++;
5b3f37f6 347 }
3e87ef69 348 if(!expand)
349 outtrack->SetNClusters(AliL3Transform::GetNRows(-1));
5b3f37f6 350 }
3e87ef69 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
5b3f37f6 361}
362
3e87ef69 363void AliL3DataCompressor::ExpandTrackData(AliL3TrackArray *tracks)
5b3f37f6 364{
3e87ef69 365 //Loop over tracks and try to assign unused clusters.
366 //Only clusters which are closer than the max. residual are taken.
5b3f37f6 367
3e87ef69 368 cout<<"Expanding "<<tracks->GetNTracks()<<" tracks"<<endl;
369 for(Int_t i=0; i<tracks->GetNTracks(); i++)
5b3f37f6 370 {
3e87ef69 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--)
5b3f37f6 380 {
3e87ef69 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:
0a86fbb7 426 Int_t temp = (Int_t)rint((xyz_cross[1]-xyz[1])/GetXYResidualStep(padrow));
3e87ef69 427 if( abs(temp) > 1<<(GetNPadBits()-1))
428 continue;
429
0a86fbb7 430 temp = (Int_t)rint((xyz_cross[2]-xyz[2])/GetZResidualStep(padrow));
3e87ef69 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 }
5b3f37f6 470 }
3e87ef69 471 track->SetNClusters(AliL3Transform::GetNRows());
472 //cout<<"Track was assigned "<<nhits<<" clusters"<<endl;
5b3f37f6 473 }
474
475}
476
3e87ef69 477void AliL3DataCompressor::WriteRemaining(Bool_t select)
5b3f37f6 478{
3e87ef69 479 //Write remaining clusters (not assigned to any tracks) to file
5b3f37f6 480
5b3f37f6 481
3e87ef69 482 if(!fKeepRemaining)
483 return;
5b3f37f6 484
3e87ef69 485 if(select)
486 SelectRemainingClusters();
487
488 Char_t filename[1024];
489
490 if(!fSinglePatch)
5b3f37f6 491 {
3e87ef69 492 cerr<<"AliL3Compressor::WriteRemaining : You have to modify this function when not running singlepatch"<<endl;
5b3f37f6 493 return;
3e87ef69 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++)
5b3f37f6 500 {
3e87ef69 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++;
b2a02bce 577 fNunusedClusters++;
3e87ef69 578 }
b2a02bce 579
3e87ef69 580 //Write the last row:
3e87ef69 581 fwrite(tempPt,size,1,outfile);
582 if(data)
583 delete [] data;
584 fclose(outfile);
585 }
5b3f37f6 586 }
3e87ef69 587 delete [] npoints;
588}
589
590void AliL3DataCompressor::SelectRemainingClusters()
591{
592 //Select which remaining clusters to write in addition to the compressed data.
0a86fbb7 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
3e87ef69 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);
5b3f37f6 607
3e87ef69 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
0a86fbb7 622 if(padrow >= nrows-1-gap-shift) continue;//save all the clusters in this region
623
3e87ef69 624 //if(padrow >= nrows-1-shift) continue;
625
626 //Save the clusters at the borders:
0a86fbb7 627 //if(xyz[1] < 3 || xyz[1] >= AliL3Transform::GetNPads(padrow)-4)
628 // continue;
3e87ef69 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 }
5b3f37f6 639
5b3f37f6 640}
641
3e87ef69 642void 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;
b2a02bce 651
652 //Write the ratio between used and unused clusters to comp file:
653 ofstream &out = *fCompRatioFile;
654 out<<fNusedClusters<<' '<<fNunusedClusters<<endl;
3e87ef69 655}
656
657
0a86fbb7 658void AliL3DataCompressor::RestoreData(Bool_t remaining_only)
5b3f37f6 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
3e87ef69 662 //final offline tracker.
5b3f37f6 663
3e87ef69 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
0a86fbb7 680 if(!remaining_only)
681 ReadUncompressedData(clusters,ncl,maxpoints);
5b3f37f6 682
3e87ef69 683 if(fKeepRemaining)
684 ReadRemaining(clusters,ncl,maxpoints);
bd53cfb7 685
3e87ef69 686 Char_t filename[1024];
5b3f37f6 687 sprintf(filename,"%s/digitfile.root",fPath);
688 TFile *rootfile = TFile::Open(filename);
689 rootfile->cd();
690 AliTPCParam *param = (AliTPCParam*)rootfile->Get(AliL3Transform::GetParamName());
3e87ef69 691
5b3f37f6 692 AliTPCDigitsArray *darray = new AliTPCDigitsArray();
693 darray->Setup(param);
694 darray->SetClass("AliSimDigits");
3e87ef69 695 sprintf(filename,"TreeD_%s_%d",AliL3Transform::GetParamName(),fEvent);
5b3f37f6 696 Bool_t ok = darray->ConnectTree(filename);
697 if(!ok)
698 {
699 cerr<<"AliL3DataCompressor::RestoreData : Problems connecting tree"<<endl;
700 return;
701 }
3e87ef69 702
703 fOutputFile->cd();
704
705 AliTPCClustersArray *carray = new AliTPCClustersArray();
706 carray->Setup(param);
707 carray->SetClusterType("AliTPCcluster");
708 carray->MakeTree();
5b3f37f6 709
3e87ef69 710 Int_t totcounter=0;
711 for(Int_t slice=0; slice<=35; slice++)
5b3f37f6 712 {
3e87ef69 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++)
5b3f37f6 725 {
3e87ef69 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)
5b3f37f6 734 {
3e87ef69 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);
5b3f37f6 747
3e87ef69 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;
bd53cfb7 754
3e87ef69 755 for(Int_t lab=0; lab<3; lab++)
bd53cfb7 756 {
3e87ef69 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)
bd53cfb7 765 {
3e87ef69 766 falseid++;
767 //AliL3Transform::Local2Global(temp,slice);
768 //cout<<"slice "<<slice<<" padrow "<<padrow<<" y "<<temp[1]<<" z "<<temp[2]<<" label "<<c->GetLabel(0)<<endl;
bd53cfb7 769 }
bd53cfb7 770 }
3e87ef69 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
803void 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
856void 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;
bd53cfb7 878
3e87ef69 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++)
5b3f37f6 902 {
5b3f37f6 903
3e87ef69 904 Float_t xyz[3] = {AliL3Transform::Row2X(padrow),points[j].fY,points[j].fZ};
5b3f37f6 905
3e87ef69 906 AliL3Transform::Local2Raw(xyz,sector,row);
bd53cfb7 907
3e87ef69 908 if(ncl[slice] >= maxpoints)
5b3f37f6 909 {
3e87ef69 910 cerr<<"AliL3DataCompressor::ReadRemaining : Too many clusters"<<endl;
911 exit(5);
5b3f37f6 912 }
3e87ef69 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]++;
bd53cfb7 921 }
3e87ef69 922 Byte_t *dPt = (Byte_t*)tempPt;
923 Int_t size = sizeof(AliL3RemainingRow) + tempPt->fNClusters*sizeof(AliL3RemainingCluster);
924 dPt += size;
925 tempPt = (AliL3RemainingRow*)dPt;
5b3f37f6 926 }
3e87ef69 927
928 mem.Free();
5b3f37f6 929 }
930 }
3e87ef69 931}
932
933void 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
971Int_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;
5b3f37f6 997
3e87ef69 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;
5b3f37f6 1000
3e87ef69 1001 return 1;
5b3f37f6 1002}
1003