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 }