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