Last changes from I. Belikov - the tracking performs track refitting
[u/mrichter/AliRoot.git] / HLT / comp / AliL3Compress.cxx
CommitLineData
5bf93292 1//$Id$
2
3// Author: Anders Vestbo <mailto:vestbo$fi.uib.no>
4//*-- Copyright &copy ASV
5
6#include <stdio.h>
7#include <stream.h>
8#include <stdlib.h>
029912b7 9#include <TH1.h>
10#include <TH2.h>
11#include <TRandom.h>
5bf93292 12
13#include "AliL3Compress.h"
14#include "AliL3TrackArray.h"
15#include "AliL3ModelTrack.h"
029912b7 16#include "AliL3Transform.h"
17#include "AliL3MemHandler.h"
5bf93292 18#include "bitio.h"
19
029912b7 20//_____________________________________________________________
21//
22// AliL3Compress
23//
24// Class for compressing and uncompressing data.
25
5bf93292 26ClassImp(AliL3Compress)
27
28AliL3Compress::AliL3Compress()
29{
95a00d93 30 fTracks=0;
029912b7 31 SetBitNumbers(7,7,10,4);
32 fSlice =0;
33 fPatch=0;
34 fDigits=0;
35 fDPt=0;
36}
37
38AliL3Compress::AliL3Compress(Int_t slice,Int_t patch,Int_t pad,Int_t time,Int_t charge,Int_t shape)
39{
40 fSlice=slice;
41 fPatch=patch;
42 SetBitNumbers(pad,time,charge,shape);
43 fTracks=0;
44 fDigits=0;
45 fDPt=0;
5bf93292 46}
47
48AliL3Compress::~AliL3Compress()
49{
95a00d93 50 if(fTracks)
51 delete fTracks;
029912b7 52 if(fDigits)
53 delete [] fDigits;
54 if(fDPt)
55 delete [] fDPt;
56}
57
58void AliL3Compress::SetBitNumbers(Int_t pad,Int_t time,Int_t charge,Int_t shape)
59{
60 fNumPadBits=pad;
61 fNumTimeBits=time;
62 fNumChargeBits=charge;
63 fNumShapeBits=shape;
5bf93292 64}
65
95a00d93 66void AliL3Compress::WriteFile(AliL3TrackArray *tracks,Char_t *filename)
5bf93292 67{
95a00d93 68 FILE *file = fopen(filename,"w");
5bf93292 69 Short_t ntracks = tracks->GetNTracks();
95a00d93 70 //cout<<"Writing "<<ntracks<<" tracks to file"<<endl;
71
72 Int_t count=0;
73 AliL3ClusterModel *clusters=0;
74 AliL3TrackModel *model=0;
5bf93292 75 for(Int_t i=0; i<ntracks; i++)
76 {
77 AliL3ModelTrack *track = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
78 if(!track) continue;
029912b7 79
80 track->FillModel();
95a00d93 81 model = track->GetModel();
82 if(model->fNClusters==0) continue;
83 clusters = track->GetClusters();
84 //cout<<"Writing "<<(int)model->fNClusters<<" clusters"<<endl;
5bf93292 85 if(fwrite(model,sizeof(AliL3TrackModel),1,file)!=1) break;
95a00d93 86 //cout<<"Writing "<<(int)model->fNClusters<<" clusters to file"<<endl;
87 if(fwrite(clusters,model->fNClusters*sizeof(AliL3ClusterModel),1,file)!=1) break;
029912b7 88 //track->Print();
95a00d93 89 count++;
029912b7 90
5bf93292 91 }
95a00d93 92 cout<<"Wrote "<<count<<" tracks "<<endl;
5bf93292 93 fclose(file);
94}
95
95a00d93 96void AliL3Compress::ReadFile(Char_t *filename)
5bf93292 97{
95a00d93 98 FILE *file = fopen(filename,"r");
029912b7 99 if(!file)
100 {
101 cerr<<"Cannot open file "<<filename<<endl;
102 return;
103 }
104
95a00d93 105 if(fTracks)
106 delete fTracks;
107 fTracks = new AliL3TrackArray("AliL3ModelTrack");
5bf93292 108
029912b7 109 cout<<"Reading file "<<filename<<endl;
5bf93292 110 while(!feof(file))
111 {
95a00d93 112 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->NextTrack();
029912b7 113 track->Init(fSlice,fPatch);
5bf93292 114 AliL3TrackModel *model = track->GetModel();
115 AliL3ClusterModel *clusters = track->GetClusters();
029912b7 116 //cout<<"Reading model "<<(int)model<<endl;
5bf93292 117 if(fread(model,sizeof(AliL3TrackModel),1,file)!=1) break;
029912b7 118 //cout<<"Reading clusters "<<(int)clusters<<endl;
5bf93292 119 if(fread(clusters,(model->fNClusters)*sizeof(AliL3ClusterModel),1,file)!=1) break;
029912b7 120 //cout<<"Filling track"<<endl;
121 track->FillTrack();
122 //track->Print();
5bf93292 123 }
029912b7 124
125 fTracks->RemoveLast();
95a00d93 126 cout<<"Read "<<fTracks->GetNTracks()<<" tracks from file"<<endl;
5bf93292 127 fclose(file);
128}
129
95a00d93 130void AliL3Compress::CompressFile(Char_t *infile,Char_t *outfile)
5bf93292 131{
132
95a00d93 133 BIT_FILE *output = OpenOutputBitFile(outfile);
134 FILE *input = fopen(infile,"r");
5bf93292 135
136 AliL3TrackModel track;
137 AliL3ClusterModel cluster;
138 Int_t temp;
029912b7 139 Short_t power;
140
141 Int_t timeo,pado,chargeo,shapeo;
142 timeo=pado=chargeo=shapeo=0;
5bf93292 143 while(!feof(input))
144 {
145 if(fread(&track,sizeof(AliL3TrackModel),1,input)!=1) break;
95a00d93 146
147 if(output->mask != 0x80) //Write the current byte to file.
148 {
149 if(putc(output->rack,output->file )!=output->rack)
150 cerr<<"AliL3Compress::ComressFile : Error writing to bitfile"<<endl;
151 output->mask=0x80;
029912b7 152 output->rack=0;
95a00d93 153 }
029912b7 154
155 //Write track parameters:
5bf93292 156 fwrite(&track,sizeof(AliL3TrackModel),1,output->file);
5bf93292 157 for(Int_t i=0; i<track.fNClusters; i++)
158 {
159 if(fread(&cluster,sizeof(AliL3ClusterModel),1,input)!=1) break;
029912b7 160
161 //Write empty flag:
162 temp = (Int_t)cluster.fPresent;
163 OutputBit(output,temp);
164 if(!temp) continue;
165
166 //Write time information:
5bf93292 167 temp = (Int_t)cluster.fDTime;
168 if(temp<0)
169 OutputBit(output,0);
170 else
171 OutputBit(output,1);
029912b7 172 power = 1<<fNumTimeBits;
173 if(abs(temp)>=power)
174 {
175 timeo++;
176 temp=power - 1;
177 }
178 OutputBits(output,abs(temp),fNumTimeBits);
179
180 //Write pad information:
5bf93292 181 temp = (Int_t)cluster.fDPad;
182 if(temp<0)
183 OutputBit(output,0);
184 else
185 OutputBit(output,1);
029912b7 186 power = 1<<fNumPadBits;
187 if(abs(temp)>=power)
188 {
189 pado++;
190 temp=power - 1;
191 }
192 OutputBits(output,abs(temp),fNumPadBits);
193
194 //Write charge information:
5bf93292 195 temp = (Int_t)cluster.fDCharge;
029912b7 196 if(temp<0)
197 OutputBit(output,0);
198 else
199 OutputBit(output,1);
200 power = 1<<fNumChargeBits;
201 if(abs(temp)>=power)
202 {
203 chargeo++;
204 temp=power - 1;
205 }
206 OutputBits(output,abs(temp),fNumChargeBits);
5bf93292 207
029912b7 208 //Write shape information:
209 temp = (Int_t)cluster.fDSigmaY2;
210 power = 1<<fNumShapeBits;
211 if(abs(temp) >= power)
212 {
213 shapeo++;
214 temp = power - 1;
215 }
216 OutputBits(output,abs(temp),fNumShapeBits);
217
218 temp = (Int_t)cluster.fDSigmaZ2;
219 if(abs(temp) >= power)
220 {
221 shapeo++;
222 temp=power - 1;
223 }
224 OutputBits(output,abs(temp),fNumShapeBits);
5bf93292 225 }
5bf93292 226 }
95a00d93 227
5bf93292 228 fclose(input);
229 CloseOutputBitFile(output);
029912b7 230
231 cout<<endl<<"There was following number of overflows: "<<endl
232 <<"Pad "<<pado<<endl
233 <<"Time "<<timeo<<endl
234 <<"Charge "<<chargeo<<endl
235 <<"Shape "<<shapeo<<endl;
5bf93292 236}
237
95a00d93 238void AliL3Compress::ExpandFile(Char_t *infile,Char_t *outfile)
5bf93292 239{
95a00d93 240 BIT_FILE *input = OpenInputBitFile(infile);
241 FILE *output = fopen(outfile,"w");
242
243 AliL3TrackModel trackmodel;
244 AliL3ClusterModel *clusters=0;
029912b7 245 Int_t count=0;
5bf93292 246
029912b7 247 clusters = new AliL3ClusterModel[(NumRows[fPatch])];
95a00d93 248 while(!feof(input->file))
5bf93292 249 {
029912b7 250 input->mask=0x80;//make sure we read a new byte from file.
95a00d93 251
029912b7 252 //Read and write track:
95a00d93 253 if(fread(&trackmodel,sizeof(AliL3TrackModel),1,input->file)!=1) break;
254 fwrite(&trackmodel,sizeof(AliL3TrackModel),1,output);
029912b7 255
256 for(Int_t i=0; i<NumRows[fPatch]; i++)
95a00d93 257 {
258 Int_t temp,sign;
029912b7 259
260 //Read empty flag:
95a00d93 261 temp = InputBit(input);
262 if(!temp)
263 {
264 clusters[i].fPresent=kFALSE;
265 continue;
266 }
267 clusters[i].fPresent=kTRUE;
029912b7 268
269 //Read time information:
95a00d93 270 sign=InputBit(input);
029912b7 271 temp = InputBits(input,fNumTimeBits);
95a00d93 272 if(!sign)
273 temp*=-1;
274 clusters[i].fDTime = temp;
029912b7 275
276 //Read pad information:
95a00d93 277 sign=InputBit(input);
029912b7 278 temp = InputBits(input,fNumPadBits);
95a00d93 279 if(!sign)
280 temp*=-1;
281 clusters[i].fDPad = temp;
029912b7 282
283 //Read charge information:
284 sign = InputBit(input);
285 temp=InputBits(input,fNumChargeBits);
286 if(!sign)
287 temp*=-1;
95a00d93 288 clusters[i].fDCharge = temp;
029912b7 289
290 //Read shape information:
291 temp = InputBits(input,fNumShapeBits);
292 clusters[i].fDSigmaY2 = temp;
293
294 temp = InputBits(input,fNumShapeBits);
295 clusters[i].fDSigmaZ2 = temp;
95a00d93 296 }
029912b7 297
298
299 count++;
95a00d93 300 fwrite(clusters,(trackmodel.fNClusters)*sizeof(AliL3ClusterModel),1,output);
029912b7 301
5bf93292 302 }
029912b7 303
304 delete [] clusters;
95a00d93 305 fclose(output);
5bf93292 306 CloseInputBitFile(input);
307}
029912b7 308
309void AliL3Compress::CreateDigitArray(Int_t maxnumber)
310{
311 fNUsed=0;
312 fNDigits = 0;
313 fMaxDigits=maxnumber;
314 if(fDigits) delete [] fDigits;
315 fDigits = new AliL3RandomDigitData[maxnumber];
316 if(fDPt) delete [] fDPt;
317 fDPt = new AliL3RandomDigitData*[maxnumber];
318}
319
320void AliL3Compress::RestoreData(Char_t *uncompfile)
321{
322
323 //Read the uncompressed file:
324 ReadFile(uncompfile);
325
326 CreateDigitArray(100000);
327
328 Float_t pad,time,sigmaY2,sigmaZ2;
329 Int_t charge;
330 for(Int_t j=NRows[fPatch][0]; j<=NRows[fPatch][1]; j++)
331 {
332 for(Int_t i=0; i<fTracks->GetNTracks(); i++)
333 {
334 AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->GetCheckedTrack(i);
335 if(!track) continue;
336 if(!track->GetPad(j,pad) ||
337 !track->GetTime(j,time) ||
338 !track->GetClusterCharge(j,charge) ||
339 !track->GetXYWidth(j,sigmaY2) ||
340 !track->GetZWidth(j,sigmaZ2))
341 continue;
342
343 CreateDigits(j,pad,time,charge,sigmaY2,sigmaZ2);
344 }
345 }
346
347 QSort(fDPt,0,fNDigits);
348}
349
350void AliL3Compress::PrintDigits()
351{
352 Int_t pad,time,charge,row;
353 for(Int_t i=0; i<fNDigits; i++)
354 {
355 row = fDPt[i]->fRow;
356 pad = fDPt[i]->fPad;
357 time = fDPt[i]->fTime;
358 charge = fDPt[i]->fCharge;
359 if(i>0 && row != fDPt[i-1]->fRow)
360 cout<<"---Padrow "<<row<<"---"<<endl;
361 cout<<"Pad "<<pad<<" time "<<time<<" charge "<<charge<<endl;
362 }
363}
364
365void AliL3Compress::WriteRestoredData(Char_t *remainfile,Char_t *restoredfile)
366{
367 //Get the remaining raw data array:
368 AliL3MemHandler *mem = new AliL3MemHandler();
369 mem->SetBinaryInput(remainfile);
370 UInt_t numdigits;
371 AliL3DigitRowData *origRow = mem->CompBinary2Memory(numdigits);
372 mem->CloseBinaryInput();
373
374 //Allocate memory for the merged data:
375 UInt_t size = mem->GetAllocatedSize() + fNDigits*sizeof(AliL3DigitData);
376 Byte_t *data = new Byte_t[size];
377 memset(data,0,size);
378 AliL3DigitRowData *tempRow = (AliL3DigitRowData*)data;
379
380 Int_t ndigits,action,charge;
381 UShort_t pad,time;
382 for(Int_t i=NRows[fPatch][0]; i<=NRows[fPatch][1]; i++)
383 {
384 tempRow->fRow = i;
385 ndigits=0;
386 AliL3DigitData *origDig = origRow->fDigitData;
387 AliL3DigitData *tempDig = tempRow->fDigitData;
388
389 while(1)
390 {
391 for(UInt_t j=0; j<origRow->fNDigit; j++)
392 {
393 pad = origDig[j].fPad;
394 time = origDig[j].fTime;
395 charge = origDig[j].fCharge;
396 while((action=ComparePoints(i,pad,time)) == 1)
397 {
398 tempDig[ndigits].fPad = fDPt[fNUsed]->fPad;
399 tempDig[ndigits].fTime = fDPt[fNUsed]->fTime;
400 tempDig[ndigits].fCharge = fDPt[fNUsed]->fCharge;
401 ndigits++;
402 fNUsed++;
403 }
404 if(action == 0)
405 {
406 tempDig[ndigits].fPad = pad;
407 tempDig[ndigits].fTime = time;
408 tempDig[ndigits].fCharge = charge;
409 ndigits++;
410 }
411 }
412 if(fNUsed >= fNDigits) break;
413 if(fDPt[fNUsed]->fRow != i) //we are on a new row
414 break;
415 tempDig[ndigits].fPad = fDPt[fNUsed]->fPad;
416 tempDig[ndigits].fTime = fDPt[fNUsed]->fTime;
417 tempDig[ndigits].fCharge = fDPt[fNUsed]->fCharge;
418 ndigits++;
419 fNUsed++;
420 }
421 tempRow->fNDigit = ndigits;
422 Int_t size = sizeof(AliL3DigitData)*tempRow->fNDigit + sizeof(AliL3DigitRowData);
423 Byte_t *byte_pt = (Byte_t*)tempRow;
424 byte_pt += size;
425 tempRow = (AliL3DigitRowData*)byte_pt;
426 mem->UpdateRowPointer(origRow);
427 }
428
429 mem->Free();
430 mem->SetBinaryOutput(restoredfile);
431 mem->Memory2CompBinary((UInt_t)NumRows[fPatch],(AliL3DigitRowData*)data);
432 mem->CloseBinaryOutput();
433
434 delete [] data;
435 delete mem;
436}
437
438void AliL3Compress::CreateDigits(Int_t row,Float_t pad,Float_t time,Int_t charge,Float_t sigmaY2,Float_t sigmaZ2)
439{
440 //Create raw data out of the cluster.
441
442 AliL3Transform *tr = new AliL3Transform();
443 TRandom *random = new TRandom();
444
445 Int_t entries=10000;
446 TH1F *hist1 = new TH1F("hist1","",tr->GetNPads(row),0,tr->GetNPads(row)-1);
447 TH1F *hist2 = new TH1F("hist2","",tr->GetNTimeBins(),0,tr->GetNTimeBins()-1);
448 TH2F *hist3 = new TH2F("hist3","",tr->GetNPads(row),0,tr->GetNPads(row)-1,tr->GetNTimeBins(),0,tr->GetNTimeBins()-1);
449
450 //Convert back the sigmas:
451 Float_t padw,timew;
452 if(fPatch < 3)
453 padw = tr->GetPadPitchWidthLow();
454 else
455 padw = tr->GetPadPitchWidthUp();
456 timew = tr->GetZWidth();
457
458 if(fPatch < 3)
459 sigmaY2 = sigmaY2/2.07;
460 sigmaY2 = sigmaY2/0.108;
461 sigmaY2 = sigmaY2/(padw*padw);
462 sigmaY2 = sigmaY2 - 1./12;
463
464 if(fPatch < 3)
465 sigmaZ2 = sigmaZ2/1.77;
466 sigmaZ2 = sigmaZ2/0.169;
467 sigmaZ2 = sigmaZ2/(timew*timew);
468 sigmaZ2 = sigmaZ2 - 1./12;
469
470 if(sigmaY2 <= 0 || sigmaZ2 <= 0)
471 {
472 cerr<<"AliL3Compress::CreateDigits() : Wrong sigmas : "<<sigmaY2<<" "<<sigmaZ2<<endl;
473 return;
474 }
475
476 //Create the distributions in pad and time:
477 for(Int_t i=0; i<entries; i++)
478 {
479 hist1->Fill(random->Gaus(pad,sqrt(sigmaY2)));
480 hist2->Fill(random->Gaus(time,sqrt(sigmaZ2)));
481 }
482
483 //Create the cluster:
484 Int_t bin1,bin2;
485 Double_t content1,content2,dpad,dtime;
486 for(Int_t i=0; i<hist1->GetEntries(); i++)
487 {
488 bin1 = hist1->GetBin(i);
489 content1 = hist1->GetBinContent(bin1);
490 if((Int_t)content1==0) continue;
491 content1 = charge*content1/entries;
492 dpad = hist1->GetBinCenter(bin1);
493 for(Int_t j=0; j<hist2->GetEntries(); j++)
494 {
495 bin2 = hist2->GetBin(j);
496 content2 = hist2->GetBinContent(bin2);
497 if((Int_t)content2==0) continue;
498 content2 = content1*content2/entries;
499 dtime = hist2->GetBinCenter(bin2);
500 hist3->Fill(dpad,dtime,content2);
501 }
502 }
503
504 //Fill it into the digit array:
505 for(Int_t i=0; i<hist3->GetNbinsX(); i++)
506 {
507 for(Int_t j=0; j<hist3->GetNbinsY(); j++)
508 {
509 bin1 = hist3->GetBin(i,j);
510 content1 = hist3->GetBinContent(bin1);
511 if((Int_t)content1 < 3) continue;
512 if(content1 >= 1024)
513 content1 = 1023;
514 if(fNDigits >= fMaxDigits)
515 {
516 cerr<<"AliL3Compress::CreateDigits() : Array index out of range : "<<fNDigits<<endl;
517 return;
518 }
519 fDigits[fNDigits].fCharge=(Int_t)content1;
520 fDigits[fNDigits].fRow = row;
521 fDigits[fNDigits].fPad = (Int_t)hist3->GetXaxis()->GetBinCenter(i);
522 fDigits[fNDigits].fTime = (Int_t)hist3->GetYaxis()->GetBinCenter(j);
523 fDPt[fNDigits] = &fDigits[fNDigits];
524 fNDigits++;
525 }
526 }
527
528 delete random;
529 delete hist1;
530 delete hist2;
531 delete hist3;
532 delete tr;
533}
534
535void AliL3Compress::QSort(AliL3RandomDigitData **a, Int_t first, Int_t last)
536{
537
538 // Sort array of AliL3RandomDigitData pointers using a quicksort algorithm.
539 // Uses CompareDigits() to compare objects.
540 // Thanks to Root!
541
542 static AliL3RandomDigitData *tmp;
543 static int i; // "static" to save stack space
544 int j;
545
546 while (last - first > 1) {
547 i = first;
548 j = last;
549 for (;;) {
550 while (++i < last && CompareDigits(a[i], a[first]) < 0)
551 ;
552 while (--j > first && CompareDigits(a[j], a[first]) > 0)
553 ;
554 if (i >= j)
555 break;
556
557 tmp = a[i];
558 a[i] = a[j];
559 a[j] = tmp;
560 }
561 if (j == first) {
562 ++first;
563 continue;
564 }
565 tmp = a[first];
566 a[first] = a[j];
567 a[j] = tmp;
568 if (j - first < last - (j + 1)) {
569 QSort(a, first, j);
570 first = j + 1; // QSort(j + 1, last);
571 } else {
572 QSort(a, j + 1, last);
573 last = j; // QSort(first, j);
574 }
575 }
576}