Checking in for the weekend. Compressing/uncompressing works. Restoring data - buildi...
[u/mrichter/AliRoot.git] / HLT / comp / AliL3Compress.cxx
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>
9 #include <TH1.h>
10 #include <TH2.h>
11 #include <TRandom.h>
12
13 #include "AliL3Compress.h"
14 #include "AliL3TrackArray.h"
15 #include "AliL3ModelTrack.h"
16 #include "AliL3Transform.h"
17 #include "AliL3MemHandler.h"
18 #include "bitio.h"
19
20 //_____________________________________________________________
21 //
22 //  AliL3Compress
23 //
24 // Class for compressing and uncompressing data.
25
26 ClassImp(AliL3Compress)
27
28 AliL3Compress::AliL3Compress()
29 {
30   fTracks=0;
31   SetBitNumbers(7,7,10,4);
32   fSlice =0;
33   fPatch=0;
34   fDigits=0;
35   fDPt=0;
36 }
37
38 AliL3Compress::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;
46 }
47
48 AliL3Compress::~AliL3Compress()
49 {
50   if(fTracks)
51     delete fTracks;
52   if(fDigits)
53     delete [] fDigits;
54   if(fDPt)
55     delete [] fDPt;
56 }
57
58 void 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;
64 }
65
66 void AliL3Compress::WriteFile(AliL3TrackArray *tracks,Char_t *filename)
67 {
68   FILE *file = fopen(filename,"w");
69   Short_t ntracks = tracks->GetNTracks();
70   //cout<<"Writing "<<ntracks<<" tracks to file"<<endl;
71     
72   Int_t count=0;
73   AliL3ClusterModel *clusters=0;
74   AliL3TrackModel *model=0;
75   for(Int_t i=0; i<ntracks; i++)
76     {
77       AliL3ModelTrack *track = (AliL3ModelTrack*)tracks->GetCheckedTrack(i);
78       if(!track) continue;
79             
80       track->FillModel();
81       model = track->GetModel();
82       if(model->fNClusters==0) continue;
83       clusters = track->GetClusters();
84       //cout<<"Writing "<<(int)model->fNClusters<<" clusters"<<endl;
85       if(fwrite(model,sizeof(AliL3TrackModel),1,file)!=1) break;
86       //cout<<"Writing "<<(int)model->fNClusters<<" clusters to file"<<endl;
87       if(fwrite(clusters,model->fNClusters*sizeof(AliL3ClusterModel),1,file)!=1) break;
88       //track->Print();
89       count++;
90       
91     }
92   cout<<"Wrote "<<count<<" tracks "<<endl;
93   fclose(file);
94 }
95
96 void AliL3Compress::ReadFile(Char_t *filename)
97 {
98   FILE *file = fopen(filename,"r");
99   if(!file)
100     {
101       cerr<<"Cannot open file "<<filename<<endl;
102       return;
103     }
104
105   if(fTracks)
106     delete fTracks;
107   fTracks = new AliL3TrackArray("AliL3ModelTrack");
108   
109   cout<<"Reading file "<<filename<<endl;
110   while(!feof(file))
111     {
112       AliL3ModelTrack *track = (AliL3ModelTrack*)fTracks->NextTrack();
113       track->Init(fSlice,fPatch);
114       AliL3TrackModel *model = track->GetModel();
115       AliL3ClusterModel *clusters = track->GetClusters();
116       //cout<<"Reading model "<<(int)model<<endl;
117       if(fread(model,sizeof(AliL3TrackModel),1,file)!=1) break;
118       //cout<<"Reading clusters "<<(int)clusters<<endl;
119       if(fread(clusters,(model->fNClusters)*sizeof(AliL3ClusterModel),1,file)!=1) break;
120       //cout<<"Filling track"<<endl;
121       track->FillTrack();
122       //track->Print();
123     }
124
125   fTracks->RemoveLast();
126   cout<<"Read "<<fTracks->GetNTracks()<<" tracks from file"<<endl;
127   fclose(file);
128 }
129
130 void AliL3Compress::CompressFile(Char_t *infile,Char_t *outfile)
131 {
132   
133   BIT_FILE *output = OpenOutputBitFile(outfile);
134   FILE *input = fopen(infile,"r");
135
136   AliL3TrackModel track;
137   AliL3ClusterModel cluster;
138   Int_t temp;
139   Short_t power;
140   
141   Int_t timeo,pado,chargeo,shapeo;
142   timeo=pado=chargeo=shapeo=0;
143   while(!feof(input))
144     {
145       if(fread(&track,sizeof(AliL3TrackModel),1,input)!=1) break;
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;
152           output->rack=0;
153         }
154       
155       //Write track parameters:
156       fwrite(&track,sizeof(AliL3TrackModel),1,output->file);
157       for(Int_t i=0; i<track.fNClusters; i++)
158         {
159           if(fread(&cluster,sizeof(AliL3ClusterModel),1,input)!=1) break;
160           
161           //Write empty flag:
162           temp = (Int_t)cluster.fPresent;
163           OutputBit(output,temp);
164           if(!temp) continue;
165           
166           //Write time information:
167           temp = (Int_t)cluster.fDTime;
168           if(temp<0)
169             OutputBit(output,0);
170           else
171             OutputBit(output,1);
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:
181           temp = (Int_t)cluster.fDPad;
182           if(temp<0)
183             OutputBit(output,0);
184           else
185             OutputBit(output,1);
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:
195           temp = (Int_t)cluster.fDCharge;
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);
207           
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);
225         }
226     }
227   
228   fclose(input);
229   CloseOutputBitFile(output);
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;
236 }
237
238 void AliL3Compress::ExpandFile(Char_t *infile,Char_t *outfile)
239 {
240   BIT_FILE *input = OpenInputBitFile(infile);
241   FILE *output = fopen(outfile,"w");
242
243   AliL3TrackModel trackmodel;
244   AliL3ClusterModel *clusters=0;
245   Int_t count=0;
246   
247   clusters = new AliL3ClusterModel[(NumRows[fPatch])];
248   while(!feof(input->file))
249     {
250       input->mask=0x80;//make sure we read a new byte from file.
251       
252       //Read and write track:
253       if(fread(&trackmodel,sizeof(AliL3TrackModel),1,input->file)!=1) break;
254       fwrite(&trackmodel,sizeof(AliL3TrackModel),1,output);
255       
256       for(Int_t i=0; i<NumRows[fPatch]; i++)
257         {
258           Int_t temp,sign;
259           
260           //Read empty flag:
261           temp = InputBit(input);
262           if(!temp) 
263             {
264               clusters[i].fPresent=kFALSE;
265               continue;
266             }
267           clusters[i].fPresent=kTRUE;
268           
269           //Read time information:
270           sign=InputBit(input);
271           temp = InputBits(input,fNumTimeBits);
272           if(!sign)
273             temp*=-1;
274           clusters[i].fDTime = temp;
275           
276           //Read pad information:
277           sign=InputBit(input);
278           temp = InputBits(input,fNumPadBits);
279           if(!sign)
280             temp*=-1;
281           clusters[i].fDPad = temp;
282           
283           //Read charge information:
284           sign = InputBit(input);
285           temp=InputBits(input,fNumChargeBits);
286           if(!sign)
287             temp*=-1;
288           clusters[i].fDCharge = temp;
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;
296         }
297       
298
299       count++;
300       fwrite(clusters,(trackmodel.fNClusters)*sizeof(AliL3ClusterModel),1,output);
301       
302     }
303   
304   delete [] clusters;
305   fclose(output);
306   CloseInputBitFile(input);
307 }
308
309 void 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
320 void 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
350 void 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
365 void 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
438 void 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
535 void 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 }