HLT/misc/AliTransBit.cxx

   1 ///////////////////////////////////////////////////////////////////////////////
   2 //                                                                           //
   3 //  Time Projection Chamber ADC bit compresion lookup table                  //
   4 //                                                                           //
   5 //
   6 //  Origin: Marian Ivanov , GSI Darmstadt                                   //
   7 //                                                                          //
   8
   9 //
  10 /*
  11   Conversion equation:
  12   For AliTransBit_v1
  13   dy/dx= Int_t(1+x/fX0)
  14
  15   For AliTransBit_v2
  16   y  =(2**bit0) ln(1+x/fX0) / ln(1+(2**bit1-1)/fX0)
  17
  18   where x0 is calculated in function GetOptimumX0()
  19
  20 */
  21
  22 // Example Session                                                          //
  23 /*
  24
  25 Int_t b0=10;  // original number of bits
  26   Int_t b1=8;   // compressed
  27
  28   AliTransBit_v1 trans;
  29   Int_t x0=TMath::Nint(TMath::Exp(b0*TMath::Log(2)));
  30   Int_t x1=TMath::Nint(TMath::Exp(b1*TMath::Log(2)));
  31   trans.SetBits(b0,b1);
  32   trans.FindOptimumX0();
  33   trans.Update();
  34   cout<<trans.Get0to1(x0-2)<<"\n";
  35   cout<<trans.Get1to0(x1-2)<<"\n";
  36
  37   // to produce table of
  38   for( Int_t i=0;i<x1;i++) cout<<i<<"\t"<<trans.Get1to0(i)<<"\n"; > table1to0.txt
  39   for( Int_t i=0;i<x0;i++) cout<<i<<"\t"<<trans.Get0to1(i)<<"\n"; > table0to1.txt
  40
  41   for( Int_t i=0;i<x1-1;i++) cout<<i<<"\t"<<trans.Get1to0(i+1)-trans.Get1to0(i)<<"\n"; > tabled.txt
  42
  43
  44
  45 */
  46 //Begin_Html                                                                //
  47 /*                                                                          //
  48 <img src="gif/AliTPCTransBit.gif">
  49 */
  50 //End_Html
  51 //                                                                           //
  52 //                                                                          //
  53 ///////////////////////////////////////////////////////////////////////////////
  54
  55
  56 #include "AliTransBit.h"
  57
  58 #include <iostream.h>
  59 #include <math.h>
  60
  61 ClassImp(AliTransBit)
  62 ClassImp(AliTransBit_v1)
  63 ClassImp(AliTransBit_v2)
  64
  65
  66 AliTransBit::AliTransBit()
  67 {
  68   //
  69   //
  70   fTable0 = 0;
  71   fTable1 = 0;
  72   fBit0   = 10;
  73   fBit1   = 8;
  74   fX0     = 0;
  75 }
  76
  77 AliTransBit::~AliTransBit()
  78 {
  79   //
  80   //default destructor
  81   if (fTable0!=0) delete [] fTable0;
  82   if (fTable1!=0) delete [] fTable1;
  83 }
  84
  85
  86
  87
  88
  89 Double_t AliTransBit_v1::FindOptimumX0()
  90 {
  91   //
  92   //find x0 for which derivation at xder1 is equal 1
  93   //
  94   Int_t x0=(Int_t)rint(exp(fBit0*log(2)));
  95   Int_t x1=(Int_t)rint(exp(fBit1*log(2)));
  96
  97   fX0 = ((x1-2)*(x1-2)/2.)/(x0-x1-1);  //starting fX0
  98   Int_t digit=0;
  99   Int_t j;
 100   Double_t weight=1;
 101   for(j=0;(j<50)&&(digit!=(x0-1));j++){
 102     Int_t olddigit=digit;
 103     digit=0;
 104     for (Int_t  i=0; i<x1-1;i++) digit+= Int_t(1+Double_t(i)/fX0);
 105     fX0*= (1.+weight*Double_t(digit)/Double_t(x0-1))/(1.+weight);
 106     if ( ((olddigit-(x0-1)) * (digit-(x0-1))) <0 )
 107       weight*=0.25;
 108     // cout<<j<<"\t"<<fX0<<"\t"<<digit<<"\n";
 109   }
 110   //  cout<<"End of iteration "<<j<<"\n";
 111   //cout<<"digit..."<<digit<<"\n";
 112   return fX0;
 113 }
 114
 115
 116 void AliTransBit_v1::Update()
 117 {
 118   //
 119   //construct lookup tables for loosy compresion from
 120   if (fX0<1) fX0 = FindOptimumX0();
 121   Int_t x0=(Int_t)rint(exp(fBit0*log(2)));
 122   Int_t x1=(Int_t)rint(exp(fBit1*log(2)));
 123
 124   //fTable0 - conversion from bit0 coding to bit1 coding
 125   if (fTable0!=0) delete fTable0;
 126   fTable0= new Int_t[x0];
 127   //fTable1 - conversion table from bit1 to bit0 coding
 128   if (fTable1!=0) delete [] fTable1;
 129   fTable1= new Int_t[x1];
 130   //
 131   Int_t digit=0;
 132   for (Int_t  i=0; i<x1-1;i++){
 133     Int_t ddig=Int_t(1+Double_t(i)/fX0);
 134     for (Int_t j=0;j<ddig;j++) if ((digit+j)<x0) fTable0[digit+j] = i;
 135     fTable1[i]= digit+ddig/2;
 136     digit+= ddig;
 137   }
 138   for (Int_t i=digit;i<x0-1;i++) fTable0[i]=x1-2;
 139   fTable1[x1-1]=x0-1; //overflow
 140   fTable0[x0-1]=x1-1; //overflow
 141   return;
 142 }
 143
 144
 145 Double_t AliTransBit_v2::FindOptimumX0()
 146 {
 147   //
 148   //find x0 for which derivation at xder1 is equal 1
 149   //
 150   const Float_t xder1=1;
 151   const Float_t dx=0.1;
 152
 153   Float_t x0=exp(fBit0*log(2));
 154   Float_t x1=exp(fBit1*log(2));
 155   Float_t deriv = 0;
 156   Float_t x;
 157   for (x=x1;( (x>1)&&(deriv<1)) ;x-=dx)
 158     {
 159       deriv = (x1-1)/( log(1.+x0/x) *x *(1+xder1/x));
 160     }
 161   x+=dx/2.;
 162   fX0 = x;
 163   return fX0;
 164 }
 165
 166
 167 void AliTransBit_v2::Update()
 168 {
 169   //
 170   //construct lookup tables for loosy compresion from
 171   if (fX0<1) fX0 = FindOptimumX0();
 172   //Float_t x0=(Int_t)rint(exp(fBit0*log(2)));
 173   //Float_t x1=(Int_t)rint(exp(fBit1*log(2)));
 174   Int_t x0=(Int_t)rint(exp(fBit0*log(2)));
 175   Int_t x1=(Int_t)rint(exp(fBit1*log(2)));
 176
 177   //fTable0 - conversion from bit0 coding to bit1 coding
 178   if (fTable0!=0) delete fTable0;
 179   fTable0= new Int_t[x0];
 180
 181   //fTable1 - conversion table from bit1 to bit0 coding
 182   if (fTable1!=0) delete [] fTable1;
 183   fTable1= new Int_t[x1];
 184
 185   //
 186   Int_t i;
 187
 188   for ( i=0; i<x0;i++)
 189       fTable0[i] =(Int_t)rint((x1-0.501)*log(1.+Float_t(i)/fX0)/
 190                                  log(1.+(x0-1)/fX0));
 191
 192   Int_t old0=-1;
 193   Int_t old1=-1;
 194   Int_t new1;
 195   for ( i=0; i<x0;i++){
 196       new1 = fTable0[i];
 197       if (new1!=old1){
 198         if (old1>=0) fTable1[old1] =(old0+i)/2;
 199         old0 = i;
 200         old1 = new1;
 201       }
 202   }
 203   fTable1[old1]=(Int_t)rint((Float_t)(old0+x0)/2);
 204
 205   return;
 206 }