Functions.cc
Go to the documentation of this file.
1 // LOCAL
2 #include <unordered_map>
3 #include <src/Functions.h>
4 
5 // OFFLINE
6 #include <utl/Math.h>
7 #include <math.h>
8 #include <cmath>
9 
10 
11 namespace un2 {
12 
13  inline
14  double
15  StD(const std::vector<double>& v, const double mean)
16  {
17  // assumes v.size() >= 2
18  double var = 0;
19  for (const auto& x : v)
20  var += utl::Sqr(x - mean);
21  return std::sqrt(var / (v.size() - 1));
22  }
23 
24 
25  double
26  StandardDeviation(const std::vector<double>& v, const double mean)
27  {
28  if (v.size() < 2)
29  return std::numeric_limits<double>::quiet_NaN();
30  return StD(v, mean);
31  }
32 
33 
34  double
35  StandardDeviation(const std::vector<double>& v)
36  {
37  if (v.size() < 2)
38  return std::numeric_limits<double>::quiet_NaN();
39  return StD(v, Mean(v));
40  }
41 
42 
43  // universality signal model paper 2016, ave engel roth schulz, eq. (12)
44  double
45  ModifiedGH(/* x */ const double deltaX,
46  /*p[0]*/ const double sRef,
47  /*p[1]*/ const double lgE,
48  /*p[2]*/ const double gamma,
49  /*p[3]*/ const double deltaX0,
50  /*p[4]*/ const double deltaXMax,
51  /*p[5]*/ const double deltaXRef,
52  /*p[6]*/ const double lambda0,
53  /*p[7]*/ const double lambda1,
54  /*p[8]*/ const double lambda2)
55  {
56  const double lgERef = 19;
57  const double lambda = lambda0 + deltaX * (lambda1 + deltaX * lambda2);
58  const double s0 =
59  sRef *
60  std::pow(10., (lgE - lgERef) * gamma) *
61  std::pow((deltaX - deltaX0) / (deltaXRef - deltaX0), (deltaXMax - deltaX0) / lambda) *
62  std::exp((deltaXRef - deltaX) / lambda);
63  return s0;
64  }
65 
66 
67  double
68  NKG(/* x */ const double r,
69  /*p[0]*/ const double n,
70  /*p[1]*/ const double rG,
71  /*p[2]*/ const double s)
72  {
73  const double rr = r / rG;
74  const double rho =
75  n / (utl::kTwoPi*rG*rG) * tgamma(4.5 - s) / (tgamma(s) * tgamma(4.5 - 2*s)) *
76  pow(rr, s - 2) * pow(1 + rr, s - 4.5);
77  return (rho > 0) ? rho : 0;
78  }
79 
80 
81  double
83  const double fThrsh,
84  const double cThrsh,
85  const double bg)
86  {
87  const double x = log(fThrsh / f) / log(cThrsh);
88  const double w = atan(x) / M_PI + 0.5;
89  return f * (1 - w) + bg * w;
90  }
91 
92 
93  double
94  AzimuthCorrectionFactor(const double theta,
95  const double psi,
96  const double r,
97  const double r0,
98  const double l)
99  {
100  const double x = r / r0 * cos(psi) * tan(theta);
101  return 1 + x * (1 + 0.5 * x * l);
102  }
103 
104 
105  // this is NOT empiric
106  double
107  TwoFormAzimuthCorrectionFactor(const double theta,
108  const double psi,
109  const double M1,
110  const double M3)
111  {
112  const double cp = cos(psi);
113  const double val =
114  1 +
115  theta * cp * (M1 +
116  theta * (M1*M1 * 3./2 * cp +
117  theta * M3 * 1./6 * (15 * cp*cp - 2)));
118  return val;
119  }
120 
121 
122  double
123  Gaussian(const double x, const double mu, const double sigma)
124  {
125  return std::exp(-0.5 * utl::Sqr((x - mu) / sigma)) / (utl::kSqrt2Pi * sigma);
126  }
127 
128 
129  // the spikyFunction is supposed to cover the shape of the SSD trace as a sum of gauss peaks
130  // of a constant width (sigma), separated in mu by a constant step size delta OR of individual
131  // step size, given by a vector of individual positions on the t-axis. the height of
132  // the individual gaussians is scaled by nvec[i], where nvec.size() dictates the number of peaks
133  double
134  SpikyFunction(const double x,
135  const double sigma,
136  const double delta,
137  const std::vector<double>& nvec)
138  {
139  double val = 0;
140  for (unsigned int i = 0, n = nvec.size(); i < n; ++i)
141  val += nvec[i] * Gaussian(x, delta * i, sigma);
142  return val;
143  }
144 
145 
146  double
147  SpikyFunction(const double x,
148  const std::vector<double>& sigmavec,
149  const std::vector<double>& muvec,
150  const std::vector<double>& nvec)
151  {
152  double val = 0;
153  for (unsigned int i = 0, n = nvec.size(); i < n; ++i)
154  val += nvec[i] * Gaussian(x, muvec[i], sigmavec[i]);
155  return val;
156  }
157 
158 
159  double
160  SignalStartTime(const double r, const double psi, const double theta,
161  const double DX, const double Xvert, const double hs, const double hs2)
162  {
163  const double ct = std::cos(theta);
164  const double tt = std::tan(theta);
165  const double t1 = hs - hs2 * cos(psi) * sin(theta);
166  const double t2 = std::log(Xvert / (ct * AnalyticXmax(r, psi, theta, DX)));
167  const double rsp = r * std::sin(psi);
168  const double rcp = r * std::cos(psi);
169  const double t3 = rcp - t1 * tt * t2;
170  const double val =
171  -std::sqrt((t1*t1 * t2*t2 / (ct*ct))) + std::sqrt(t1*t1 * t2*t2 + rsp*rsp + t3*t3);
172  return val;
173  }
174 
175 
176  double
177  SignalStartTimePlaneFront(const double r, const double psi, const double theta,
178  const double DX, const double Xvert, const double hs, const double hs2, const double height)
179  {
180  // TODO: refactor
181  /* old and shitty
182  const double ehs = hs - hs2 * cos(psi) * sin(theta);
183  const double rs = r * sin(psi);
184  const double lo = std::log(Xvert/(cos(theta)*AnalyticXmax(r, psi, theta, DX)));
185  const double t = r*cos(psi) - ehs * tan(theta) * lo;
186  const double rf = std::sqrt(ehs*ehs * lo*lo + rs*rs + t*t);
187 
188  return rf - std::sqrt(rf*rf - r*r);
189  */
190  const double xMax = AnalyticXmax(r, psi, theta, DX, Xvert, hs, hs2, height);
191 
192  const double hProj = r * cos(psi) * sin(theta);
193  const double hXmax = AnalyticHeightFromDepth(xMax, theta, Xvert, hs, hs2);
194 
195  const double delta = (hXmax - (hProj+height))/cos(theta);
196  const double val = std::sqrt(delta*delta+r*r) - delta;
197  return val;
198  }
199 
200 
201  double
202  LogarithmOfRayleighPDF(const double x, const double sigma)
203  {
204  if (x < 0)
205  return std::numeric_limits<double>::min();
206  const double s2 = sigma * sigma;
207  const double val = std::log(x/s2) - x*x / (2*s2);
208  return val;
209  }
210 
211 
212  double
213  LogarithmOfTruncatedGaussianPDF(const double x, const double mu, const double sigma)
214  {
215  const double s2 = sigma * sigma;
216  const double r = (x - mu) / sigma;
217  const double val =
218  - 0.5*std::log(2*M_PI * s2)
219  - 0.5*r*r
220  - std::log(0.5*(1 - std::erfc(mu / (std::sqrt(2)*sigma))));
221  return val;
222  }
223 
224 
225  double
226  LogNormalPDF(const double x, const double mu, const double sigma)
227  {
228  if (x <= 0)
229  return 0;
230  return std::exp(-0.5 * utl::Sqr((log(x) - mu) / sigma)) / (utl::kSqrt2Pi * sigma * x);
231  }
232 
233 
234  double
235  LogNormalCDF(const double x, const double mu, const double sigma)
236  {
237  if (x <= 0)
238  return 0;
239  return 0.5 + 0.5 * std::erf((std::log(x) - mu) / (std::sqrt(2) * sigma));
240  }
241 
242 
243  double
244  LogNormalMuFromModeAndSigma(const double mode, const double sigma, const double ns_bin)
245  {
246  return log(mode / ns_bin) + sigma*sigma;
247  }
248 
249 
250  double
251  LogNormalMuFromExpAndStd(const double exp_val, const double std, const double ns_bin)
252  {
253  const double le = log(exp_val / ns_bin);
254  return 0.5 * (2 * le - log(1 + std::exp(-2*(le - log(std / ns_bin)))));
255  }
256 
257 
258  double
259  LogNormalSigmaFromExpAndStd(const double exp_val, const double std, const double ns_bin)
260  {
261  const double sn = std / ns_bin;
262  const double en = exp_val / ns_bin;
263  return std::sqrt(log((sn*sn + en*en) / (en*en)));
264  }
265 
266  double
267  LogNormalTimeQuantileFromMuAndSigma(const double mu, const double sigma, const double quantile, const double ns_bin)
268  {
269  return ns_bin * std::exp(mu + utl::kSqrt2 * sigma * ErfInv(2*quantile-1));
270  }
271 
272 
273  double
274  LogNormalMuFromTimeQuantileAndSigma(const double timeQuantile, const double quantile, const double sigma, const double ns_bin)
275  {
276 
277  // old and wrong version
278  //return std::log(timeQuantile/ns_bin) + utl::kSqrt2 * sigma * ErfcInv(2*quantile);
279  return std::log(timeQuantile/ns_bin) - utl::kSqrt2 * sigma * ErfInv(2*quantile-1);
280 
281  }
282 
283  // stolen from http://libit.sourceforge.net/math_8c-source.html
284  double ErfInv(double x)
285  {
286  const double erfinv_a3 = -0.140543331;
287  const double erfinv_a2 = 0.914624893;
288  const double erfinv_a1 = -1.645349621;
289  const double erfinv_a0 = 0.886226899;
290  const double erfinv_b4 = 0.012229801;
291  const double erfinv_b3 = -0.329097515;
292  const double erfinv_b2 = 1.442710462;
293  const double erfinv_b1 = -2.118377725;
294  const double erfinv_b0 = 1;
295  const double erfinv_c3 = 1.641345311;
296  const double erfinv_c2 = 3.429567803;
297  const double erfinv_c1 = -1.62490649;
298  const double erfinv_c0 = -1.970840454;
299  const double erfinv_d2 = 1.637067800;
300  const double erfinv_d1 = 3.543889200;
301  const double erfinv_d0 = 1.;
302  double x2, r, y;
303  int sign_x;
304 
305  if (x < -1 || x > 1)
306  return NAN;
307 
308  if (x == 0)
309  return 0;
310 
311  if (x > 0)
312  sign_x = 1;
313  else {
314  sign_x = -1;
315  x = -x;
316  }
317 
318  if (x <= 0.7) {
319 
320  x2 = x * x;
321  r =
322  x * (((erfinv_a3 * x2 + erfinv_a2) * x2 + erfinv_a1) * x2 + erfinv_a0);
323  r /= (((erfinv_b4 * x2 + erfinv_b3) * x2 + erfinv_b2) * x2 +
324  erfinv_b1) * x2 + erfinv_b0;
325  }
326  else {
327  y = sqrt (-log ((1 - x) / 2));
328  r = (((erfinv_c3 * y + erfinv_c2) * y + erfinv_c1) * y + erfinv_c0);
329  r /= ((erfinv_d2 * y + erfinv_d1) * y + erfinv_d0);
330  }
331 
332  r = r * sign_x;
333  x = x * sign_x;
334 
335  r -= (erf (r) - x) / (2 / sqrt (M_PI) * std::exp(-r * r));
336  r -= (erf (r) - x) / (2 / sqrt (M_PI) * std::exp(-r * r));
337 
338  return r;
339  }
340 
341  double ErfcInv(const double x)
342  {
343  return ErfInv(1-x);
344  }
345 
346 
347  double
348  GammaPDF(const double x, const double mu, const double theta)
349  {
350  return pow(x, mu - 1) * std::exp(-(x / theta)) / (tgamma(mu) * pow(theta, mu));
351  }
352 
353 
354  double
355  AnalyticDX(const double r, const double psi, const double theta,
356  const double Xmax, const double Xvert, const double hs, const double hs2, const double height)
357  {
358  const double projHeight = r*cos(psi)*sin(theta);
359  return (Xvert / cos(theta)) * std::exp(-(projHeight + (height - 1400)) / (hs - hs2 * (projHeight+height))) - Xmax;
360  }
361 
362 
363  double
364  AnalyticXmax(const double r, const double psi, const double theta,
365  const double DX, const double Xvert, const double hs, const double hs2, const double height)
366  {
367  const double projHeight = r*cos(psi)*sin(theta);
368  const double val = (Xvert / cos(theta)) * std::exp(-(projHeight + (height - 1400)) / (hs - hs2 * (projHeight+height))) - DX;
369  return (val <= 0.1) ? 0.1 : val;
370  }
371 
372 
373  double
374  AnalyticHeightFromDepth(const double depth, const double theta, const double xVertGround,
375  const double hs, const double hs2)
376  {
377  const double logdepth = log(xVertGround/(depth * cos(theta)));
378  const double val = hs / (1/logdepth + hs2) + 1400.;
379  return val;
380 
381  }
382 
383  // dependend on sin²θ and log10(recE0/eV)
384  double BiasCorrectionXmax(const double sst, const double lgE, const bool realEvent)
385  {
386 
387  if (realEvent) {
388 
389  //const double corr_sth = -15.67427143194855 + sst * -105.37047967258333 + sst*sst * 982.876974428145 + sst*sst*sst * -1321.9250208741992;
390  const double corr_lgE = (10 - 16.7 *(lgE-19.));
391  const double corr_sth = (0 + sst * 0 + sst*sst * 0);
392  return -(corr_sth + corr_lgE);
393 
394  } else {
395 
396  const double corr_sth = -(21.500192617693102 + sst * 112.80972683392125 + sst*sst * -261.03423400887186);
397  const double corr_lgE = (10 - 16.7 *(lgE-19.));
398  return -(corr_sth + corr_lgE);
399 
400  }
401 
402  }
403 
404  double BiasCorrectionRmu(const double sst, const double lgE, const bool realEvent)
405  {
406 
407  if (realEvent) {
408 
409  //const double corr_sth = 0.03050172683372075 + sst * 1.220036864771516 + sst*sst * -7.288933408076791 + sst*sst*sst * 8.92870121242658;
410  //const double corr_lgE = -1.452508874274341 + lgE * 0.07534666688375187;
411  const double corr_sth = 0 * sst;
412  const double corr_lgE = (-0.0462 + (lgE-19.) * 0.08688);
413  return -(corr_sth+corr_lgE);
414 
415  } else {
416 
417  const double corr_sth = 0 * sst;
418  //const double corr_lgE = -(-2.3137564320935122 + lgE * 0.11851604224552405);
419  const double corr_lgE = (-0.0462 + (lgE-19.) * 0.08688);
420  return -(corr_sth + corr_lgE);
421 
422  }
423 
424  }
425 
426  double
427  GHGauss(const double r, const double psi, const double theta, const double Xmax, const double lambda)
428  {
429  return std::exp((-utl::Sqr(AnalyticDX(r, psi, theta, Xmax)) +
430  utl::Sqr(AnalyticDX(r, M_PI/2, theta, Xmax))) / (2 * Xmax * lambda));
431  }
432 
433  double
434  avgXmax0(const double lgE, const bool realEvent)
435  {
436 
437  // return realEvent ? 54.66 * lgE - 275 : 54.66 * lgE - 247.82;
438  return realEvent ? (797.16 + (lgE-19.) * 47.1) : (797.16 + (lgE-19.) * 47.1);
439 
440  }
441 
442  double
443  avglnRmu0(const double lgE, const bool realEvent)
444  {
445 
446  // based on EPOS-LHC
447  // 0.26 is what the first measurements using hybrids were implying
448  // should be updated soon
449  return realEvent ? 0.26 + 0.0179 + (lgE - 19.) * -0.0176 : 0.0179 + (lgE - 19.) * -0.0176;
450 
451  }
452 
453 
454  double
455  XmaxFromRmu(const double Rmu, const double lgE, const bool realEvent)
456  {
457 
458  /* returns the median value for Xmax as a function of Rmu*/
459  const double Xmax0 = avgXmax0(lgE, realEvent);
460  const double lnRmu0 = avglnRmu0(lgE, realEvent);
461  const double lambda = 21.28 + (lgE - 19.) * -0.645; // these values must be extracted from simulations
462  const double beta = 0.065 + (lgE - 19.) * -0.011; // using the universality parametrization at hand!!!
463 
464  const double deltalnRmu = (log(Rmu) - lnRmu0) > 0.24 ? 0.24 : log(Rmu) - lnRmu0;
465 
466  const double Xmax = Xmax0 - lambda / beta * (deltalnRmu);
467  return Xmax;
468 
469  }
470 
471 
472  double
473  SigmaXmaxFromRmu(const double Rmu, const double lgE, const bool realEvent)
474  {
475  const double lnRmu0 = avglnRmu0(lgE, realEvent);
476  const double deltalnRmu = (log(Rmu) - lnRmu0) > 0.24 ? 0.24 : log(Rmu) - lnRmu0;
477  const double sigmaXmax = 90. - 70. * (deltalnRmu)/0.24;
478  return sigmaXmax;
479 
480  }
481 
482  double
483  XmaxDistribution(const double Rmu, const double Xmax, const double lgE, const bool realEvent)
484  {
485  // approx. Xmax distribution as a function of Rmu for the constrained fit
486  // you can try either one of the following three
487  const double mXmax = XmaxFromRmu(Rmu, lgE, realEvent);
488  const double sigmaXmax = SigmaXmaxFromRmu(Rmu, lgE, realEvent);
489 
490  // gumbel
491  const double beta = sigmaXmax / M_PI * pow(6, 0.5);
492  const double mu = mXmax + log(log(2)) * beta;
493  const double z = (Xmax - mu)/beta;
494  return 1/beta*exp(-(z+exp(-z)));
495 
496  // triangle
497  //const double tri = 700 - abs(mXmax - Xmax);
498  //return tri
499 
500  // gaussian
501  //return 1./utl::kSqrt2Pi/sigmaXmax * std::exp(-0.5*utl::Sqr((mXmax - Xmax) / sigmaXmax));
502 
503  }
504 
505  double lnA(const double Rmu, const double Xmax, const double lgE, const bool realEvent)
506  {
507  const double lnRmu = Rmu<1 ? Rmu-1 : log(Rmu); // avoid crazy outliers using ln approximation
508 
509  const double Xmax0 = avgXmax0(lgE, realEvent);
510  const double lnRmu0 = avglnRmu0(lgE, realEvent);
511 
512  const double lambda = 21.28 + (lgE - 19.) * -0.645; // these values must be extracted from simulations
513  const double beta = 0.065 + (lgE - 19.) * -0.011; // using the universality parametrization at hand!!!
514 
515  // calculate the average expected dispersion for Xmax and lnRmu from the approx. values of
516  // sigma(Xmax) and sigma(lnRmu) for iron and proton, linearly interpolate in between
517  // a lower limit is used so that the distribution does not become artificially thin
518  // no upper sigmaXmax limit is needed, in this case the Xmax fit is then just less constrained
519  // which is ok for protons, which do flucutate severly
520  const double avgSigmaXmax = 70 + 40/(log(56)*lambda) *(Xmax-Xmax0);
521  const double avgSigmaLnRmu = 0.138 - 0.05/(log(56)*beta) * (lnRmu-lnRmu0);
522 
523  const double sigmaXmax = avgSigmaXmax < 20 ? 20 : avgSigmaXmax;
524  const double sigmaLnRmu = avgSigmaLnRmu < 0.08 ? 0.08 : avgSigmaLnRmu;
525 
526  const double phi0 = utl::Sqr(beta*sigmaXmax/sigmaLnRmu)/lambda;
527 
528  const double lnA = (phi0 * (lnRmu - lnRmu0) - beta * (Xmax-Xmax0))/(beta*(lambda+phi0));
529  return lnA;
530 
531  }
532 
533 }
double TwoFormAzimuthCorrectionFactor(const double theta, const double psi, const double M1, const double M3)
Definition: Functions.cc:107
constexpr double kSqrt2Pi
Definition: MathConstants.h:32
double StandardDeviation(const std::vector< double > &v, const double mean)
Definition: Functions.cc:26
double BiasCorrectionXmax(const double sst, const double lgE, const bool realEvent)
Definition: Functions.cc:384
double BiasCorrectionRmu(const double sst, const double lgE, const bool realEvent)
Definition: Functions.cc:404
double SpikyFunction(const double x, const double sigma, const double delta, const std::vector< double > &nvec)
Definition: Functions.cc:134
constexpr T Sqr(const T &x)
double AnalyticXmax(const double r, const double psi, const double theta, const double DX, const double Xvert, const double hs, const double hs2, const double height)
Definition: Functions.cc:364
double ErfcInv(const double x)
Definition: Functions.cc:341
double XmaxFromRmu(const double Rmu, const double lgE, const bool realEvent)
Definition: Functions.cc:455
double LogarithmOfTruncatedGaussianPDF(const double x, const double mu, const double sigma)
Definition: Functions.cc:213
double avglnRmu0(const double lgE, const bool realEvent)
Definition: Functions.cc:443
double LogNormalMuFromModeAndSigma(const double mode, const double sigma, const double ns_bin)
Definition: Functions.cc:244
double SignalStartTimePlaneFront(const double r, const double psi, const double theta, const double DX, const double Xvert, const double hs, const double hs2, const double height)
Definition: Functions.cc:177
double ActivationFunctionBackground(const double f, const double fThrsh, const double cThrsh, const double bg)
Definition: Functions.cc:82
double pow(const double x, const unsigned int i)
double ModifiedGH(const double deltaX, const double sRef, const double lgE, const double gamma, const double deltaX0, const double deltaXMax, const double deltaXRef, const double lambda0, const double lambda1, const double lambda2)
Definition: Functions.cc:45
double ErfInv(double x)
Definition: Functions.cc:284
constexpr double s
Definition: AugerUnits.h:163
double XmaxDistribution(const double Rmu, const double Xmax, const double lgE, const bool realEvent)
Definition: Functions.cc:483
double LogNormalMuFromExpAndStd(const double exp_val, const double std, const double ns_bin)
Definition: Functions.cc:251
double LogNormalSigmaFromExpAndStd(const double exp_val, const double std, const double ns_bin)
Definition: Functions.cc:259
constexpr double kTwoPi
Definition: MathConstants.h:27
double Gaussian(const double x, const double mu, const double sigma)
Definition: Functions.cc:123
a t3: algo, second, usecond and a vector of &lt;t3stat&gt;
Definition: XbT2.h:29
a second level trigger
Definition: XbT2.h:8
constexpr double kSqrt2
Definition: MathConstants.h:30
double lnA(const double Rmu, const double Xmax, const double lgE, const bool realEvent)
Definition: Functions.cc:505
double AnalyticDX(const double r, const double psi, const double theta, const double Xmax, const double Xvert, const double hs, const double hs2, const double height)
Definition: Functions.cc:355
double SigmaXmaxFromRmu(const double Rmu, const double lgE, const bool realEvent)
Definition: Functions.cc:473
double GHGauss(const double r, const double psi, const double theta, const double Xmax, const double lambda)
Definition: Functions.cc:427
double LogNormalTimeQuantileFromMuAndSigma(const double mu, const double sigma, const double quantile, const double ns_bin)
Definition: Functions.cc:267
double avgXmax0(const double lgE, const bool realEvent)
Definition: Functions.cc:434
double LogNormalMuFromTimeQuantileAndSigma(const double timeQuantile, const double quantile, const double sigma, const double ns_bin)
Definition: Functions.cc:274
double AnalyticHeightFromDepth(const double depth, const double theta, const double xVertGround, const double hs, const double hs2)
Definition: Functions.cc:374
double Mean(const std::vector< double > &v)
Definition: Functions.h:31
double LogNormalPDF(const double x, const double mu, const double sigma)
Definition: Functions.cc:226
double LogNormalCDF(const double x, const double mu, const double sigma)
Definition: Functions.cc:235
double GammaPDF(const double x, const double mu, const double theta)
Definition: Functions.cc:348
double StD(const std::vector< double > &v, const double mean)
Definition: Functions.cc:15
double LogarithmOfRayleighPDF(const double x, const double sigma)
Definition: Functions.cc:202
double AzimuthCorrectionFactor(const double theta, const double psi, const double r, const double r0, const double l)
Definition: Functions.cc:94
double NKG(const double r, const double n, const double rG, const double s)
Definition: Functions.cc:68
double SignalStartTime(const double r, const double psi, const double theta, const double DX, const double Xvert, const double hs, const double hs2)
Definition: Functions.cc:160

, generated on Tue Sep 26 2023.