1 #include <ecl/digitization/shaperdsp.h>
6 using namespace Belle2::ECL;
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};
12 o << s.t <<
" " << s.c0 <<
" " << s.c0 <<
" " << s.c1 <<
" " << s.s1 <<
" " << s.e0 <<
" " << s.e1 <<
" " << s.es <<
" " << s.ed;
19 o << s0 <<
" " << s.et0 <<
" " << s.et1;
23 void ShaperDSP_t::sv123shift_t::init(
double _t,
const ShaperDSP_t& _p)
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);
38 double c0r = r.c0 * c0 - r.s0 * s0;
39 double s0r = r.c0 * s0 + r.s0 * c0;
43 double c1r = r.c1 * c1 - r.s1 * s1;
44 double s1r = r.c1 * s1 + r.s1 * c1;
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;
72 void ShaperDSP_t::shaperdspshift_t::init(
double _t,
const ShaperDSP_t& _p)
74 sv123shift_t::init(_t, _p);
75 et0 = exp(-t * _p.
_dt0);
76 et1 = exp(-t * _p.
_dt1);
83 sv123shift_t::operator +=(r0);
107 dd_t
operator +(
const dd_t& c0,
const dd_t& c1)
109 return dd_t(c0.first + c1.first, c0.second + c1.second);
114 return dd_t(a * c.first, a * c.second);
117 void ShaperDSP_t::Sv123_init(
double t01,
double tb1,
double t02,
double tb2,
double td1,
double ts1)
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;
124 dr = (ts1 - td1) / td1;
125 if (std::abs(dr) >= 0.0000001) {
137 dr = ((t01 - t02) * (t01 - t02) + (tb1 - tb2) * (tb1 - tb2)) / (t01 * t01 + tb1 * tb1);
140 if (dr < 0.0000000001) {
142 dks0 = dks1 * 1.00001;
144 dks0 = dks1 * 0.99999;
160 dzna = (dksm2 + dwm * dwm) * (dksm2 + dwp * dwp);
162 das0 = dw1 * (dksm2 + dwp * dwm);
163 dac0 = -2 * dksm * dw0 * dw1;
164 das1 = dw0 * (dksm2 - dwp * dwm);
169 dzn0 = dcs0 * dcs0 + dsn0 * dsn0;
171 dsn0s = (dsn0 * das0 - dcs0 * dac0) / dzn0;
172 dcs0s = (dcs0 * das0 + dsn0 * dac0) / dzn0;
176 dzn0 = dcs0 * dcs0 + dsn0 * dsn0;
178 dsn1s = (dsn0 * das1 - dcs0 * dac1) / dzn0;
179 dcs1s = (dcs0 * das1 + dsn0 * dac1) / dzn0;
183 dzn0 = dcs0 * dcs0 + dsn0 * dsn0;
185 dsn0d = (dsn0 * das0 - dcs0 * dac0) / dzn0;
186 dcs0d = (dcs0 * das0 + dsn0 * dac0) / dzn0;
190 dzn0 = dcs0 * dcs0 + dsn0 * dsn0;
192 dsn1d = (dsn0 * das1 - dcs0 * dac1) / dzn0;
193 dcs1d = (dcs0 * das1 + dsn0 * dac1) / dzn0;
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);
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;
218 return (f0 + f1) + (fd - fs);
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;
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;
233 double f = (f0 + f1) + (fd - fs);
234 double fp = (f0p + f1p) - (fdp - fsp);
240 if (t0.validshift(_tp)) {
241 double ft0p = Sv123(t0 + _tp), ft0 = 0, ft0m = 0;
244 if (t0.validshift(_tm))
245 ft0m = Sv123(t0 + _tm);
247 return _w0 * ft0 + _w1 * (ft0p + ft0m);
254 if (t0.validshift(_tp)) {
255 dd_t ft0p = ddSv123(t0 + _tp), ft0 = dd_t(0, 0), ft0m = dd_t(0, 0);
258 if (t0.validshift(_tm))
259 ft0m = ddSv123(t0 + _tm);
261 return _w0 * ft0 + _w1 * (ft0p + ft0m);
269 double f = Sv123_filtered(t00);
271 double z = t0.t * _dt1;
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);
285 dd_t f = ddSv123_filtered(t00);
287 double z = t0.t * _dt1;
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);
295 double up = -_dt0 * u + ((_dt1 * t0.et0) * (t0.et1 * z2)) * ((z2 * z) * (1 / 120.));
296 f.second -= _ccc * up;
301 void ShaperDSP_t::init(
const double* s,
double unitscale)
309 Sv123_init(t01, tb1, t02, tb2, td1, ts1);
325 _tp.init(_filterdt, *
this);
326 _tm.init(-_filterdt, *
this);
329 void ShaperDSP_t::init(
const double* s)
334 void ShaperDSP_t::init(
const std::vector<double>& s,
double unitscale)
337 init(s.data(), unitscale);
339 init(_defs, unitscale);
342 double ShaperDSP_t::operator()(
double t)
const
345 return ShaperDSP(t0);
348 double ShaperDSP_t::operator()(
double* x,
double*)
350 return (*
this)(x[0]);
353 void ShaperDSP_t::settimestride(
double dt)
355 _tstride.init(dt , *
this);
358 void ShaperDSP_t::setseedoffset(
double dt)
360 _toffset.init(dt , *
this);
363 void ShaperDSP_t::settimeseed(
double t0)
365 _tzero.init(t0 - _toff, *
this);
368 void ShaperDSP_t::nextseed()
373 void ShaperDSP_t::fillarray(
int n,
double* s)
const
377 for (
int i = 0; i < n; i++, t0 += _tstride) {
378 s[i] = ShaperDSP(t0);
382 void ShaperDSP_t::fillarray(
int n, dd_t* s)
const
386 for (
int i = 0; i < n; i++, t0 += _tstride) {
387 s[i] = ddShaperDSP(t0);
391 void ShaperDSP_t::fillarray(
double t,
int n,
double* s)
const
395 for (
int i = 0; i < n; i++, t0 += _tstride) {
396 s[i] = ShaperDSP(t0);
400 void ShaperDSP_t::fillarray(
double t,
int n, dd_t* s)
const
404 for (
int i = 0; i < n; i++, t0 += _tstride) {
405 s[i] = ddShaperDSP(t0);
409 void ShaperDSP_t::fillvector(std::vector<double>& s)
const
411 fillarray(s.size(), s.data());
414 void ShaperDSP_t::fillvector(std::vector<dd_t>& s)
const
416 fillarray(s.size(), s.data());
419 void ShaperDSP_t::fillvector(
double t, std::vector<double>& s)
const
421 fillarray(t, s.size(), s.data());
424 void ShaperDSP_t::fillvector(
double t, std::vector<dd_t>& s)
const
426 fillarray(t, s.size(), s.data());