Belle II Software  release-05-01-25
ECLCompress.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2017 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Alexei Sibidanov *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <ecl/digitization/ECLCompress.h>
12 #include <ecl/digitization/EclConfiguration.h>
13 #include <math.h>
14 
15 using namespace Belle2;
16 using namespace ECL;
17 
18 unsigned int Belle2::ECL::ilog2(unsigned int v) // find the log base 2 of 32-bit v
19 {
20  static const unsigned char MultiplyDeBruijnBitPosition[32] = {
21  0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
22  8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31
23  };
24 
25  v |= v >> 1; // first round down to one less than a power of 2
26  v |= v >> 2;
27  v |= v >> 4;
28  v |= v >> 8;
29  v |= v >> 16;
30 
31  return MultiplyDeBruijnBitPosition[(v * 0x07C4ACDDU) >> 27];
32 }
33 
34 void ECLBaseCompress::compress(BitStream& out, const int* adc)
35 {
36  int amin = adc[0], amax = adc[0];
37  for (unsigned int i = 0; i < EclConfiguration::m_nsmp; i++) {
38  amax = std::max(amax, adc[i]);
39  amin = std::min(amin, adc[i]);
40  }
41 
42  amin &= 0x3ffff;
43  out.putNBits(amin, 18);
44  unsigned int w = ilog2(amax - amin) + 1;
45  w &= 0x1f;
46  out.putNBits(w, 5);
47 
48  for (unsigned int i = 0; i < EclConfiguration::m_nsmp; i++) {
49  unsigned int d = adc[i] - amin;
50  out.putNBits(d, w);
51  }
52 }
53 
55 {
56  unsigned int amin = out.getNBits(18), w = out.getNBits(5);
57  for (unsigned int i = 0; i < EclConfiguration::m_nsmp; i++)
58  adc[i] = amin + out.getNBits(w);
59 }
60 
61 void ECLDeltaCompress::compress(BitStream& out, const int* adc)
62 {
63  int a0 = adc[0], dmin = 0, dmax = 0;
64  for (int k = 1; k < EclConfiguration::m_nsmp; k++) {
65  int d = adc[k] - a0;
66  dmin = std::min(dmin, d);
67  dmax = std::max(dmax, d);
68  a0 = adc[k];
69  }
70 
71  unsigned int w = ilog2(std::max(2 * dmax, -2 * dmin)) + 1;
72  w &= 0xf;
73  out.putNBits(w, 4);
74 
75  a0 = adc[0] & 0x3ffff;
76  out.putNBits(a0, 18);
77 
78  unsigned int base = 1 << (w - 1);
79  for (unsigned int i = 1; i < EclConfiguration::m_nsmp; i++) {
80  unsigned int d = int(adc[i]) - a0 + base;
81  out.putNBits(d, w);
82  a0 = adc[i];
83  }
84 }
85 
87 {
88  unsigned int w = out.getNBits(4);
89  adc[0] = out.getNBits(18);
90  unsigned int base = 1 << (w - 1);
91  for (unsigned int i = 1; i < EclConfiguration::m_nsmp; i++) {
92  adc[i] = adc[i - 1] + out.getNBits(w) - base;
93  }
94 }
95 
96 namespace {
102  width_t widths_scale025[] = {
103  {7, 9, 18, 32},// 7.72447
104  {5, 7, 18, 32},// 5.45839
105  {5, 7, 18, 32},// 5.41033
106  {5, 7, 18, 32},// 5.32805
107  {5, 7, 18, 32},// 5.24319
108  {4, 6, 18, 32},// 5.01064
109  {4, 6, 18, 32},// 4.75667
110  {4, 6, 18, 32},// 4.53016
111  {4, 6, 18, 32},// 4.33823
112  {3, 5, 18, 32},// 4.08838
113  {3, 5, 18, 32},// 3.7205
114  {3, 4, 18, 32},// 3.41795
115  {2, 4, 18, 32},// 3.17582
116  {2, 4, 18, 32},// 2.76922
117  {2, 3, 18, 32},// 2.42437
118  {2, 3, 18, 32},// 2.20985
119  {2, 3, 18, 32},// 2.09761
120  {2, 3, 18, 32},// 2.04576
121  {2, 3, 18, 32},// 2.02288
122  {2, 4, 18, 32},// 2.01152
123  {2, 4, 18, 32},// 2.00643
124  {2, 4, 18, 32},// 2.00386
125  {2, 4, 18, 32},// 2.00244
126  {2, 4, 18, 32},// 2.00171
127  {2, 4, 18, 32},// 2.0012
128  {2, 4, 18, 32},// 2.00081
129  {2, 4, 18, 32},// 2.00055
130  {2, 4, 18, 32},// 2.00036
131  {2, 4, 18, 32},// 2.00023
132  {2, 3, 18, 32},// 2.00012
133  {2, 3, 18, 32} // 2.00004
134  };
135 
141  width_t widths_phs2_scale10[] = {
142  {5, 7, 9, 32},// 5.82104
143  {4, 6, 8, 32},// 4.76806
144  {4, 6, 8, 32},// 4.70815
145  {4, 6, 8, 32},// 4.61517
146  {3, 5, 7, 32},// 4.42656
147  {3, 5, 7, 32},// 4.22157
148  {3, 5, 7, 32},// 4.01412
149  {3, 5, 7, 32},// 3.80959
150  {2, 4, 6, 32},// 3.60224
151  {2, 4, 6, 32},// 3.31705
152  {2, 4, 6, 32},// 3.03457
153  {2, 3, 5, 32},// 2.71501
154  {2, 3, 5, 32},// 2.45094
155  {2, 3, 5, 32},// 2.25788
156  {2, 3, 5, 32},// 2.13303
157  {2, 3, 5, 32},// 2.06428
158  {2, 3, 5, 32},// 2.02847
159  {2, 3, 5, 32},// 2.01253
160  {1, 2, 4, 32},// 1.86085
161  {1, 2, 4, 32},// 1.68465
162  {1, 2, 4, 32},// 1.53003
163  {1, 2, 4, 32},// 1.38031
164  {1, 2, 4, 32},// 1.27103
165  {1, 2, 4, 32},// 1.18264
166  {1, 2, 4, 32},// 1.11546
167  {1, 2, 4, 32},// 1.07223
168  {1, 2, 4, 32},// 1.04641
169  {1, 2, 4, 32},// 1.03003
170  {1, 2, 4, 32},// 1.01772
171  {1, 2, 3, 32},// 1.01304
172  {1, 2, 4, 32},// 1.0107
173  };
174 }
175 
182 void stream_int(BitStream& OUT, int x, const width_t& w)
183 {
184  int ax = abs(x), m0 = 1 << (w.w0 - 1), m1 = 1 << (w.w1 - 1), m2 = 1 << (w.w2 - 1);
185  if (ax < m0) { // integer fits into w0 bits
186  OUT.putNBits(x, w.w0);
187  } else if (ax < m1) {// integer fits into w1 bits
188  OUT.putNBits((x << w.w0) | m0, w.w1 + w.w0); // first stream prefix showing that we are switching to the next bit width format
189  } else if (ax < m2) {// integer fits into w2 bits
190  OUT.putNBits((m1 << w.w0) | m0, w.w1 + w.w0);
191  OUT.putNBits(x, w.w2);
192  } else {// integer fits into w3 bits
193  OUT.putNBits((m1 << w.w0) | m0, w.w1 + w.w0);
194  OUT.putNBits(m2, w.w2);
195  OUT.putNBits(x, w.w3);
196  }
197 }
198 
204 int fetch_int(BitStream& IN, const width_t& w)
205 {
206  int m0 = 1 << (w.w0 - 1), m1 = 1 << (w.w1 - 1), m2 = 1 << (w.w2 - 1);
207  int t = IN.getNBits(w.w0);
208  if (t == m0) {
209  t = IN.getNBits(w.w1);
210  if (t == m1) {
211  t = IN.getNBits(w.w2);
212  if (t == m2) {
213  t = IN.getNBits(w.w3);
214  } else {
215  t = (t << (32 - w.w2)) >> (32 - w.w2);
216  }
217  } else {
218  t = (t << (32 - w.w1)) >> (32 - w.w1);
219  }
220  } else {
221  t = (t << (32 - w.w0)) >> (32 - w.w0);
222  }
223  return t;
224 }
225 
226 extern "C" {
231  void e10_31(const double* I, double* O);
232 
237  void e01_31(const double* I, double* O);
238 }
239 
240 ECLDCTCompress::ECLDCTCompress(double scale, double c0, width_t* w): m_scale(scale), m_c0(c0), m_widths(w) {}
241 
242 void ECLDCTCompress::compress(BitStream& OUT, const int* a)
243 {
244  const int N = EclConfiguration::m_nsmp;
245  double buf[N], out[N];
246  for (int k = 0; k < N; k++) buf[k] = a[k];
247  e10_31(buf, out);
248  for (int k = 0; k < N; k++) out[k] *= 1.0 / (2 * N);
249 
250  int km = N;
251  for (; km > 16; --km) if (lrint(out[km - 1]*m_scale) != 0) break;
252  OUT.putNBits(N - km, 4);
253  out[0] -= m_c0;
254  for (int k = 0; k < km; k++) {
255  int t = lrint(out[k] * m_scale);
256  stream_int(OUT, t, m_widths[k]);
257  }
258 }
259 
261 {
262  const double iscale = 1 / m_scale;
263  const int N = EclConfiguration::m_nsmp;
264  int nz = in.getNBits(4);
265  for (int i = 0; i < N - nz; i++) adc[i] = fetch_int(in, m_widths[i]);
266  for (int i = N - nz; i < N; i++) adc[i] = 0;
267 
268  double c[N], out[N];
269  for (int k = 0; k < N; k++) c[k] = adc[k] * iscale;
270  c[0] += m_c0;
271 
272  e01_31(c, out);
273 
274  for (int k = 0; k < N; k++) adc[k] = lrint(out[k]);
275 }
276 
277 ECLCompress* Belle2::ECL::selectAlgo(int compAlgo)
278 {
279  ECLCompress* comp = NULL;
280  if (compAlgo == 1) {
281  comp = new ECLBaseCompress;
282  } else if (compAlgo == 2) {
283  comp = new ECLDeltaCompress;
284  } else if (compAlgo == 3) {
285  comp = new ECLDCTCompress(1, 3012, widths_phs2_scale10);
286  } else if (compAlgo == 4) {
287  comp = new ECLDCTCompress(0.25, 3144, widths_scale025);
288  }
289  return comp;
290 }
Belle2::ECL::ECLDCTCompress::m_scale
const double m_scale
Scale factor for quantization.
Definition: ECLCompress.h:168
Belle2::ECL::ECLDCTCompress::m_c0
const double m_c0
Average waveform amplitude.
Definition: ECLCompress.h:169
Belle2::ECL::ECLBaseCompress
ECL waveform compression/decompression to/from the BitStream storage with the BASE algorithm.
Definition: ECLCompress.h:139
e10_31
void e10_31(const R *I, R *O)
DCT-II or "the" DCT transformation of 31-point signal This function contains 320 FP additions,...
Definition: DiscreteCosineTransform_31points.cc:117
Belle2::ECL::EclConfiguration::m_nsmp
static constexpr int m_nsmp
number of ADC measurements for signal fitting
Definition: EclConfiguration.h:51
e01_31
void e01_31(const R *I, R *O)
DCT-III or "the inverse" DCT transformation of 31-point signal This function contains 320 FP addition...
Definition: DiscreteCosineTransform_31points.cc:589
Belle2::ECL::BitStream
Bit stream struct.
Definition: ECLCompress.h:39
Belle2::ECL::ECLDeltaCompress::compress
void compress(BitStream &out, const int *adc) override
Compress the ECL waveform.
Definition: ECLCompress.cc:61
Belle2::ECL::width_t
Bit widths for the prefix coding to encode integers which are mainly concentrated around zero and pro...
Definition: ECLCompress.h:108
Belle2::ECL::ECLDeltaCompress
ECL waveform compression/decompression to/from the BitStream storage with the DELTA algorithm.
Definition: ECLCompress.h:148
Belle2::ECL::ECLDCTCompress::m_widths
const width_t * m_widths
Bit widths for the DCT coefficients for prefix encoding.
Definition: ECLCompress.h:170
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECL::ECLCompress
Abstract class (interface) for ECL waveform compression/decompression to/from the BitStream storage.
Definition: ECLCompress.h:118
Belle2::ECL::BitStream::getNBits
unsigned int getNBits(unsigned int n)
Fetch n bits.
Definition: ECLCompress.h:74
Belle2::ECL::ECLDCTCompress
ECL waveform compression/decompression to/from the BitStream storage based on the Discrete Cosine Tra...
Definition: ECLCompress.h:157
Belle2::ECL::ECLDCTCompress::ECLDCTCompress
ECLDCTCompress(double scale, double c0, width_t *w)
Constructor for DCT based compression algorithm.
Definition: ECLCompress.cc:240
Belle2::ECL::ECLBaseCompress::uncompress
void uncompress(BitStream &out, int *adc) override
Decompress the ECL waveform.
Definition: ECLCompress.cc:54
Belle2::ECL::BitStream::putNBits
void putNBits(unsigned int value, unsigned int n)
Push n least significant bits of "value" to the stream.
Definition: ECLCompress.h:61
Belle2::ECL::ECLDCTCompress::uncompress
void uncompress(BitStream &out, int *adc) override
Decompress the ECL waveform.
Definition: ECLCompress.cc:260
Belle2::ECL::ECLDeltaCompress::uncompress
void uncompress(BitStream &out, int *adc) override
Decompress the ECL waveform.
Definition: ECLCompress.cc:86
Belle2::ECL::ECLBaseCompress::compress
void compress(BitStream &out, const int *adc) override
Compress the ECL waveform.
Definition: ECLCompress.cc:34
Belle2::ECL::ECLDCTCompress::compress
void compress(BitStream &out, const int *adc) override
Compress the ECL waveform.
Definition: ECLCompress.cc:242