Belle II Software  release-05-02-19
shaperdsp.cc
1 #include <ecl/digitization/shaperdsp.h>
2 #include <cmath>
3 #include <iostream>
4 
5 using namespace std;
6 using namespace Belle2::ECL;
7 
8 const double ShaperDSP_t::_defs[] = {0.5, 0.6483, 0.4017, 0.3741, 0.8494, 0.00144547, 4.7071, 0.8156, 0.5556, 0.2752}; // default parameters
9 
10 ostream& operator<<(ostream& o, const ShaperDSP_t::sv123shift_t& s)
11 {
12  o << s.t << " " << s.c0 << " " << s.c0 << " " << s.c1 << " " << s.s1 << " " << s.e0 << " " << s.e1 << " " << s.es << " " << s.ed;
13  return o;
14 }
15 
16 ostream& operator<<(ostream& o, const ShaperDSP_t::shaperdspshift_t& s)
17 {
18  const ShaperDSP_t::sv123shift_t& s0 = static_cast<const ShaperDSP_t::sv123shift_t&>(s);
19  o << s0 << " " << s.et0 << " " << s.et1;
20  return o;
21 }
22 
23 void ShaperDSP_t::sv123shift_t::init(double _t, const ShaperDSP_t& _p)
24 {
25  t = _t;
26  sincos(t * _p._dw0, &s0, &c0);
27  sincos(t * _p._dw1, &s1, &c1);
28  e0 = exp(-t * _p._dks0);
29  e1 = exp(-t * _p._dks1);
30  es = exp(-t * _p._ds);
31  ed = exp(-t * _p._dd);
32 }
33 
34 ShaperDSP_t::sv123shift_t& ShaperDSP_t::sv123shift_t::operator +=(const ShaperDSP_t::sv123shift_t& r)
35 {
36  t += r.t;
37 
38  double c0r = r.c0 * c0 - r.s0 * s0;
39  double s0r = r.c0 * s0 + r.s0 * c0;
40  c0 = c0r;
41  s0 = s0r;
42 
43  double c1r = r.c1 * c1 - r.s1 * s1;
44  double s1r = r.c1 * s1 + r.s1 * c1;
45  c1 = c1r;
46  s1 = s1r;
47 
48  e0 *= r.e0;
49  e1 *= r.e1;
50  es *= r.es;
51  ed *= r.ed;
52  return *this;
53 }
54 
56 {
57  sv123shift_t a;
58  a.t = t + r.t;
59 
60  a.c0 = r.c0 * c0 - r.s0 * s0;
61  a.s0 = r.c0 * s0 + r.s0 * c0;
62  a.c1 = r.c1 * c1 - r.s1 * s1;
63  a.s1 = r.c1 * s1 + r.s1 * c1;
64 
65  a.e0 = e0 * r.e0;
66  a.e1 = e1 * r.e1;
67  a.es = es * r.es;
68  a.ed = ed * r.ed;
69  return a;
70 }
71 
72 void ShaperDSP_t::shaperdspshift_t::init(double _t, const ShaperDSP_t& _p)
73 {
74  sv123shift_t::init(_t, _p);
75  et0 = exp(-t * _p._dt0);
76  et1 = exp(-t * _p._dt1);
77 }
78 
79 ShaperDSP_t::shaperdspshift_t& ShaperDSP_t::shaperdspshift_t::operator +=(const ShaperDSP_t::shaperdspshift_t& r)
80 {
81  const sv123shift_t& r0 = static_cast<const ShaperDSP_t::sv123shift_t&>(r);
82 
83  sv123shift_t::operator +=(r0);
84  et0 *= r.et0;
85  et1 *= r.et1;
86  return *this;
87 }
88 
90 {
91  const sv123shift_t& r0 = static_cast<const ShaperDSP_t::sv123shift_t&>(r);
93  sv123shift_t& a0 = static_cast<ShaperDSP_t::sv123shift_t&>(a);
94 
95  a0 = sv123shift_t::operator+(r0);
96  a.et0 = et0 * r.et0;
97  a.et1 = et1 * r.et1;
98  return a;
99 }
100 
102 {
103  ShaperDSP_t::shaperdspshift_t a = c0.operator + (c1);
104  return a;
105 }
106 
107 dd_t operator +(const dd_t& c0, const dd_t& c1)
108 {
109  return dd_t(c0.first + c1.first, c0.second + c1.second);
110 }
111 
112 dd_t operator *(double a, const dd_t& c)
113 {
114  return dd_t(a * c.first, a * c.second);
115 }
116 
117 void ShaperDSP_t::Sv123_init(double t01, double tb1, double t02, double tb2, double td1, double ts1)
118 {
119  double dks0, dks1, dksm,
120  dw0, dw1, dwp, dwm, das1, dac1, das0, dac0, dzna, dksm2, ds, dd,
121  dcs0, dsn0, dzn0, td, ts, dr,
122  dcs0s, dsn0s, dcs0d, dsn0d, dcs1s, dsn1s, dcs1d, dsn1d;
123 
124  dr = (ts1 - td1) / td1;
125  if (std::abs(dr) >= 0.0000001) {
126  td = td1;
127  ts = ts1;
128  } else {
129  td = td1;
130  if (ts1 > td1) {
131  ts = td1 * 1.00001;
132  } else {
133  ts = td1 * 0.99999;
134  }
135  }
136 
137  dr = ((t01 - t02) * (t01 - t02) + (tb1 - tb2) * (tb1 - tb2)) / (t01 * t01 + tb1 * tb1);
138  dks0 = 1 / t01;
139  dks1 = 1 / t02;
140  if (dr < 0.0000000001) {
141  if (dks0 > dks1) {
142  dks0 = dks1 * 1.00001;
143  } else {
144  dks0 = dks1 * 0.99999;
145  }
146  }
147 
148  dksm = dks1 - dks0;
149 
150  ds = 1 / ts;
151  dd = 1 / td;
152 
153  dw0 = 1 / tb1;
154  dw1 = 1 / tb2;
155  dwp = dw0 + dw1;
156  dwm = dw1 - dw0;
157 
158  dksm2 = dksm * dksm;
159 
160  dzna = (dksm2 + dwm * dwm) * (dksm2 + dwp * dwp);
161 
162  das0 = dw1 * (dksm2 + dwp * dwm);
163  dac0 = -2 * dksm * dw0 * dw1;
164  das1 = dw0 * (dksm2 - dwp * dwm);
165  dac1 = -dac0;
166 
167  dsn0 = (ds - dks0);
168  dcs0 = -dw0;
169  dzn0 = dcs0 * dcs0 + dsn0 * dsn0;
170 
171  dsn0s = (dsn0 * das0 - dcs0 * dac0) / dzn0;
172  dcs0s = (dcs0 * das0 + dsn0 * dac0) / dzn0;
173 
174  dsn0 = (ds - dks1);
175  dcs0 = -dw1;
176  dzn0 = dcs0 * dcs0 + dsn0 * dsn0;
177 
178  dsn1s = (dsn0 * das1 - dcs0 * dac1) / dzn0;
179  dcs1s = (dcs0 * das1 + dsn0 * dac1) / dzn0;
180 
181  dsn0 = (dd - dks0);
182  dcs0 = -dw0;
183  dzn0 = dcs0 * dcs0 + dsn0 * dsn0;
184 
185  dsn0d = (dsn0 * das0 - dcs0 * dac0) / dzn0;
186  dcs0d = (dcs0 * das0 + dsn0 * dac0) / dzn0;
187 
188  dsn0 = (dd - dks1);
189  dcs0 = -dw1;
190  dzn0 = dcs0 * dcs0 + dsn0 * dsn0;
191 
192  dsn1d = (dsn0 * das1 - dcs0 * dac1) / dzn0;
193  dcs1d = (dcs0 * das1 + dsn0 * dac1) / dzn0;
194 
195  _dw0 = dw0;
196  _dw1 = dw1;
197  _dks0 = dks0;
198  _dks1 = dks1;
199  _ds = ds;
200  _dd = dd;
201 
202  double sc = 1 / (dzna * (ts - td));
203  _cs0 = sc * (dsn0s - dsn0d);
204  _cc0 = sc * (dcs0s - dcs0d);
205  _cs1 = sc * (dsn1s - dsn1d);
206  _cc1 = sc * (dcs1s - dcs1d);
207  _ces = sc * (dcs0s + dcs1s);
208  _ced = sc * (dcs0d + dcs1d);
209 }
210 
211 double ShaperDSP_t::Sv123(const sv123shift_t& c) const
212 {
213  double f0 = (_cs0 * c.s0 + _cc0 * c.c0) * c.e0;
214  double f1 = (_cs1 * c.s1 + _cc1 * c.c1) * c.e1;
215  double fs = _ces * c.es;
216  double fd = _ced * c.ed;
217 
218  return (f0 + f1) + (fd - fs);
219 }
220 
221 dd_t ShaperDSP_t::ddSv123(const sv123shift_t& c) const
222 {
223  double f0 = (_cs0 * c.s0 + _cc0 * c.c0) * c.e0;
224  double f1 = (_cs1 * c.s1 + _cc1 * c.c1) * c.e1;
225  double fs = _ces * c.es;
226  double fd = _ced * c.ed;
227 
228  double f0p = (_dw0 * c.e0) * (_cs0 * c.c0 - _cc0 * c.s0) - (_dks0 * f0);
229  double f1p = (_dw1 * c.e1) * (_cs1 * c.c1 - _cc1 * c.s1) - (_dks1 * f1);
230  double fsp = _ds * fs;
231  double fdp = _dd * fd;
232 
233  double f = (f0 + f1) + (fd - fs);
234  double fp = (f0p + f1p) - (fdp - fsp);
235  return dd_t(f, fp);
236 }
237 
238 double ShaperDSP_t::Sv123_filtered(const sv123shift_t& t0) const
239 {
240  if (t0.validshift(_tp)) {
241  double ft0p = Sv123(t0 + _tp), ft0 = 0, ft0m = 0;
242  if (t0.t > 0) {
243  ft0 = Sv123(t0);
244  if (t0.validshift(_tm))
245  ft0m = Sv123(t0 + _tm);
246  }
247  return _w0 * ft0 + _w1 * (ft0p + ft0m);
248  }
249  return 0;
250 }
251 
252 dd_t ShaperDSP_t::ddSv123_filtered(const sv123shift_t& t0) const
253 {
254  if (t0.validshift(_tp)) {
255  dd_t ft0p = ddSv123(t0 + _tp), ft0 = dd_t(0, 0), ft0m = dd_t(0, 0);
256  if (t0.t > 0) {
257  ft0 = ddSv123(t0);
258  if (t0.validshift(_tm))
259  ft0m = ddSv123(t0 + _tm);
260  }
261  return _w0 * ft0 + _w1 * (ft0p + ft0m);
262  }
263  return dd_t(0, 0);
264 }
265 
266 double ShaperDSP_t::ShaperDSP(const shaperdspshift_t& t0) const
267 {
268  const sv123shift_t& t00 = static_cast<const ShaperDSP_t::sv123shift_t&>(t0);
269  double f = Sv123_filtered(t00);
270  if (t0.t > 0) {
271  double z = t0.t * _dt1;
272  double z2 = z * z;
273  double odd = 1 + z2 * (1 / 6. + z2 * (1 / 120.));
274  double evn = 1 + z2 * (1 / 2. + z2 * (1 / 24.));
275  double texp = evn + z * odd;
276  double df = _ccc * t0.et0 * (1 - t0.et1 * texp);
277  f -= df;
278  }
279  return f;
280 }
281 
282 dd_t ShaperDSP_t::ddShaperDSP(const shaperdspshift_t& t0) const
283 {
284  const sv123shift_t& t00 = static_cast<const ShaperDSP_t::sv123shift_t&>(t0);
285  dd_t f = ddSv123_filtered(t00);
286  if (t0.t > 0) {
287  double z = t0.t * _dt1;
288  double z2 = z * z;
289  double odd = 1 + z2 * (1 / 6. + z2 * (1 / 120.));
290  double evn = 1 + z2 * (1 / 2. + z2 * (1 / 24.));
291  double texp = evn + z * odd;
292  double u = t0.et0 * (1 - t0.et1 * texp);
293  f.first -= _ccc * u;
294 
295  double up = -_dt0 * u + ((_dt1 * t0.et0) * (t0.et1 * z2)) * ((z2 * z) * (1 / 120.));
296  f.second -= _ccc * up;
297  }
298  return f;
299 }
300 
301 void ShaperDSP_t::init(const double* s, double unitscale)
302 {
303  double t01 = s[2];
304  double tb1 = s[3];
305  double t02 = s[7];
306  double tb2 = s[8];
307  double td1 = s[1];
308  double ts1 = s[4];
309  Sv123_init(t01, tb1, t02, tb2, td1, ts1);
310  _dt0 = 1 / s[6];
311  _dt1 = 1 / s[2];
312 
313  _toff = s[0];
314  _w0 = 1.0 - s[9];
315  _w1 = 0.5 * s[9];
316  _ccc = s[5];
317 
318  _cs0 *= unitscale;
319  _cc0 *= unitscale;
320  _cs1 *= unitscale;
321  _cc1 *= unitscale;
322  _ces *= unitscale;
323  _ced *= unitscale;
324 
325  _tp.init(_filterdt, *this);
326  _tm.init(-_filterdt, *this);
327 }
328 
329 void ShaperDSP_t::init(const double* s)
330 {
331  init(s, 27.7221);
332 }
333 
334 void ShaperDSP_t::init(const std::vector<double>& s, double unitscale)
335 {
336  if (s.size() == 10)
337  init(s.data(), unitscale);
338  else
339  init(_defs, unitscale);
340 }
341 
342 double ShaperDSP_t::operator()(double t) const
343 {
344  shaperdspshift_t t0 = shaperdspshift_t(t - _toff, *this);
345  return ShaperDSP(t0);
346 }
347 
348 double ShaperDSP_t::operator()(double* x, double*)
349 {
350  return (*this)(x[0]);
351 }
352 
353 void ShaperDSP_t::settimestride(double dt)
354 {
355  _tstride.init(dt , *this);
356 }
357 
358 void ShaperDSP_t::setseedoffset(double dt)
359 {
360  _toffset.init(dt , *this);
361 }
362 
363 void ShaperDSP_t::settimeseed(double t0)
364 {
365  _tzero.init(t0 - _toff, *this);
366 }
367 
368 void ShaperDSP_t::nextseed()
369 {
370  _tzero += _toffset;
371 }
372 
373 void ShaperDSP_t::fillarray(int n, double* s) const
374 {
375  shaperdspshift_t t0 = _tzero;
376 
377  for (int i = 0; i < n; i++, t0 += _tstride) {
378  s[i] = ShaperDSP(t0);
379  }
380 }
381 
382 void ShaperDSP_t::fillarray(int n, dd_t* s) const
383 {
384  shaperdspshift_t t0 = _tzero;
385 
386  for (int i = 0; i < n; i++, t0 += _tstride) {
387  s[i] = ddShaperDSP(t0);
388  }
389 }
390 
391 void ShaperDSP_t::fillarray(double t, int n, double* s) const
392 {
393  shaperdspshift_t t0 = shaperdspshift_t(t - _toff, *this);
394 
395  for (int i = 0; i < n; i++, t0 += _tstride) {
396  s[i] = ShaperDSP(t0);
397  }
398 }
399 
400 void ShaperDSP_t::fillarray(double t, int n, dd_t* s) const
401 {
402  shaperdspshift_t t0 = shaperdspshift_t(t - _toff, *this);
403 
404  for (int i = 0; i < n; i++, t0 += _tstride) {
405  s[i] = ddShaperDSP(t0);
406  }
407 }
408 
409 void ShaperDSP_t::fillvector(std::vector<double>& s) const
410 {
411  fillarray(s.size(), s.data());
412 }
413 
414 void ShaperDSP_t::fillvector(std::vector<dd_t>& s) const
415 {
416  fillarray(s.size(), s.data());
417 }
418 
419 void ShaperDSP_t::fillvector(double t, std::vector<double>& s) const
420 {
421  fillarray(t, s.size(), s.data());
422 }
423 
424 void ShaperDSP_t::fillvector(double t, std::vector<dd_t>& s) const
425 {
426  fillarray(t, s.size(), s.data());
427 }
Belle2::operator<<
std::ostream & operator<<(std::ostream &output, const IntervalOfValidity &iov)
Definition: IntervalOfValidity.cc:196
Belle2::ECL::ShaperDSP_t
Class include function that calculate electronic response from energy deposit
Definition: shaperdsp.h:30
Belle2::ECL::ShaperDSP_t::_dw1
double _dw1
circular frequency of the second Bessel stage
Definition: shaperdsp.h:96
Belle2::ECL::ShaperDSP_t::_dks0
double _dks0
decrement of the first Bessel stage
Definition: shaperdsp.h:98
Belle2::ECL::ShaperDSP_t::_dd
double _dd
inverse time of the differential stage
Definition: shaperdsp.h:104
Belle2::ECL::ShaperDSP_t::shaperdspshift_t
struct for a shift of the shaper dsp
Definition: shaperdsp.h:63
Belle2::operator*
B2Vector3< DataType > operator*(DataType a, const B2Vector3< DataType > &p)
non-memberfunction Scaling of 3-vectors with a real number
Definition: B2Vector3.h:528
Belle2::ECL::ShaperDSP_t::sv123shift_t
struct to encapsulate the electronic response from energy deposit
Definition: shaperdsp.h:34
Belle2::ECL::ShaperDSP_t::_dks1
double _dks1
decrement of the second Bessel stage
Definition: shaperdsp.h:100
Belle2::ECL::ShaperDSP_t::_dt0
double _dt0
coefficient for first exponent factor
Definition: shaperdsp.h:106
Belle2::ECL::ShaperDSP_t::_dt1
double _dt1
coefficient for second exponent factor
Definition: shaperdsp.h:108
Belle2::operator+
B2Vector3< DataType > operator+(const TVector3 &a, const B2Vector3< DataType > &b)
non-memberfunction for adding a TVector3 to a B2Vector3
Definition: B2Vector3.h:535
Belle2::ECL::ShaperDSP_t::_ds
double _ds
inverse scintillation decay time
Definition: shaperdsp.h:102
Belle2::ECL::ShaperDSP_t::_dw0
double _dw0
circular frequency of the first Bessel stage
Definition: shaperdsp.h:94