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 #include "TMath.h"
  58 #include "iostream.h"
  59
  60 ClassImp(AliTransBit)
  61 ClassImp(AliTransBit_v1)
  62 ClassImp(AliTransBit_v2)
  63
  64
  65 AliTransBit::AliTransBit()
  66 {
  67   //
  68   //
  69   fTable0 = 0;
  70   fTable1 = 0;
  71   fBit0   = 10;
  72   fBit1   = 8;
  73   fX0     = 0;
  74 }
  75
  76 AliTransBit::~AliTransBit()
  77 {
  78   //
  79   //default destructor
  80   if (fTable0!=0) delete [] fTable0;
  81   if (fTable1!=0) delete [] fTable1;
  82 }
  83
  84
  85
  86
  87
  88 Double_t AliTransBit_v1::FindOptimumX0()
  89 {
  90   //
  91   //find x0 for which derivation at xder1 is equal 1
  92   //
  93   Int_t x0=TMath::Nint(TMath::Exp(fBit0*TMath::Log(2)));
  94   Int_t x1=TMath::Nint(TMath::Exp(fBit1*TMath::Log(2)));
  95
  96   fX0 = ((x1-2)*(x1-2)/2.)/(x0-x1-1);  //starting fX0
  97   Int_t digit=0;
  98   Int_t j;
  99   Double_t weight=1;
 100   for(j=0;(j<50)&&(digit!=(x0-1));j++){
 101     Int_t olddigit=digit;
 102     digit=0;
 103     for (Int_t  i=0; i<x1-1;i++) digit+= Int_t(1+Double_t(i)/fX0);
 104     fX0*= (1.+weight*Double_t(digit)/Double_t(x0-1))/(1.+weight);
 105     if ( ((olddigit-(x0-1)) * (digit-(x0-1))) <0 )
 106       weight*=0.25;
 107     // cout<<j<<"\t"<<fX0<<"\t"<<digit<<"\n";
 108   }
 109   //  cout<<"End of iteration "<<j<<"\n";
 110   //cout<<"digit..."<<digit<<"\n";
 111   return fX0;
 112 }
 113
 114
 115 void AliTransBit_v1::Update()
 116 {
 117   //
 118   //construct lookup tables for loosy compresion from
 119   if (fX0<1) fX0 = FindOptimumX0();
 120   Int_t x0=TMath::Nint(TMath::Exp(fBit0*TMath::Log(2)));
 121   Int_t x1=TMath::Nint(TMath::Exp(fBit1*TMath::Log(2)));
 122
 123   //fTable0 - conversion from bit0 coding to bit1 coding
 124   if (fTable0!=0) delete fTable0;
 125   fTable0= new Int_t[x0];
 126   //fTable1 - conversion table from bit1 to bit0 coding
 127   if (fTable1!=0) delete [] fTable1;
 128   fTable1= new Int_t[x1];
 129   //
 130   Int_t digit=0;
 131   for (Int_t  i=0; i<x1-1;i++){
 132     Int_t ddig=Int_t(1+Double_t(i)/fX0);
 133     for (Int_t j=0;j<ddig;j++) if ((digit+j)<x0) fTable0[digit+j] = i;
 134     fTable1[i]= digit+ddig/2;
 135     digit+= ddig;
 136   }
 137   for (Int_t i=digit;i<x0-1;i++) fTable0[i]=x1-2;
 138   fTable1[x1-1]=x0-1; //overflow
 139   fTable0[x0-1]=x1-1; //overflow
 140   return;
 141 }
 142
 143
 144 Double_t AliTransBit_v2::FindOptimumX0()
 145 {
 146   //
 147   //find x0 for which derivation at xder1 is equal 1
 148   //
 149   const Float_t xder1=1;
 150   const Float_t dx=0.1;
 151
 152   Float_t x0=TMath::Exp(fBit0*TMath::Log(2));
 153   Float_t x1=TMath::Exp(fBit1*TMath::Log(2));
 154   Float_t deriv = 0;
 155   Float_t x;
 156   for (x=x1;( (x>1)&&(deriv<1)) ;x-=dx)
 157     {
 158       deriv = (x1-1)/( TMath::Log(1.+x0/x) *x *(1+xder1/x));
 159     }
 160   x+=dx/2.;
 161   fX0 = x;
 162   return fX0;
 163 }
 164
 165
 166 void AliTransBit_v2::Update()
 167 {
 168   //
 169   //construct lookup tables for loosy compresion from
 170   if (fX0<1) fX0 = FindOptimumX0();
 171   //Float_t x0=TMath::Nint(TMath::Exp(fBit0*TMath::Log(2)));
 172   //Float_t x1=TMath::Nint(TMath::Exp(fBit1*TMath::Log(2)));
 173   Int_t x0=TMath::Nint(TMath::Exp(fBit0*TMath::Log(2)));
 174   Int_t x1=TMath::Nint(TMath::Exp(fBit1*TMath::Log(2)));
 175
 176   //fTable0 - conversion from bit0 coding to bit1 coding
 177   if (fTable0!=0) delete fTable0;
 178   fTable0= new Int_t[x0];
 179
 180   //fTable1 - conversion table from bit1 to bit0 coding
 181   if (fTable1!=0) delete [] fTable1;
 182   fTable1= new Int_t[x1];
 183
 184   //
 185   Int_t i;
 186
 187   for ( i=0; i<x0;i++)
 188       fTable0[i] =TMath::Nint((x1-0.501)*TMath::Log(1.+Float_t(i)/fX0)/
 189                                  TMath::Log(1.+(x0-1)/fX0));
 190
 191   Int_t old0=-1;
 192   Int_t old1=-1;
 193   Int_t new1;
 194   for ( i=0; i<x0;i++){
 195       new1 = fTable0[i];
 196       if (new1!=old1){
 197         if (old1>=0) fTable1[old1] =(old0+i)/2;
 198         old0 = i;
 199         old1 = new1;
 200       }
 201   }
 202   fTable1[old1]=TMath::Nint((Float_t)(old0+x0)/2);
 203
 204   return;
 205 }