2 #include <unordered_map>
15 StD(
const std::vector<double>& v,
const double mean)
19 for (
const auto& x : v)
29 return std::numeric_limits<double>::quiet_NaN();
38 return std::numeric_limits<double>::quiet_NaN();
50 const double deltaXMax,
51 const double deltaXRef,
56 const double lgERef = 19;
57 const double lambda = lambda0 + deltaX * (lambda1 + deltaX * lambda2);
60 std::pow(10., (lgE - lgERef) * gamma) *
61 std::pow((deltaX - deltaX0) / (deltaXRef - deltaX0), (deltaXMax - deltaX0) / lambda) *
62 std::exp((deltaXRef - deltaX) / lambda);
73 const double rr = r / rG;
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;
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;
100 const double x = r / r0 * cos(psi) * tan(theta);
101 return 1 + x * (1 + 0.5 * x * l);
112 const double cp = cos(psi);
116 theta * (M1*M1 * 3./2 * cp +
117 theta * M3 * 1./6 * (15 * cp*cp - 2)));
123 Gaussian(
const double x,
const double mu,
const double sigma)
137 const std::vector<double>& nvec)
140 for (
unsigned int i = 0, n = nvec.size(); i < n; ++i)
141 val += nvec[i] *
Gaussian(x, delta * i, sigma);
148 const std::vector<double>& sigmavec,
149 const std::vector<double>& muvec,
150 const std::vector<double>& nvec)
153 for (
unsigned int i = 0, n = nvec.size(); i < n; ++i)
154 val += nvec[i] *
Gaussian(x, muvec[i], sigmavec[i]);
161 const double DX,
const double Xvert,
const double hs,
const double hs2)
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;
178 const double DX,
const double Xvert,
const double hs,
const double hs2,
const double height)
190 const double xMax =
AnalyticXmax(r, psi, theta, DX, Xvert, hs, hs2, height);
192 const double hProj = r * cos(psi) * sin(theta);
195 const double delta = (hXmax - (hProj+height))/cos(theta);
196 const double val =
std::sqrt(delta*delta+r*r) - delta;
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);
215 const double s2 = sigma * sigma;
216 const double r = (x - mu) / sigma;
218 - 0.5*std::log(2*M_PI * s2)
220 - std::log(0.5*(1 - std::erfc(mu / (
std::sqrt(2)*sigma))));
239 return 0.5 + 0.5 * std::erf((std::log(x) - mu) / (
std::sqrt(2) * sigma));
246 return log(mode / ns_bin) + sigma*sigma;
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)))));
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)));
279 return std::log(timeQuantile/ns_bin) -
utl::kSqrt2 * sigma *
ErfInv(2*quantile-1);
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.;
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;
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);
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));
348 GammaPDF(
const double x,
const double mu,
const double theta)
350 return pow(x, mu - 1) * std::exp(-(x / theta)) / (tgamma(mu) *
pow(theta, mu));
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)
358 const double projHeight = r*cos(psi)*sin(theta);
359 return (Xvert / cos(theta)) * std::exp(-(projHeight + (height - 1400)) / (hs - hs2 * (projHeight+height))) - Xmax;
365 const double DX,
const double Xvert,
const double hs,
const double hs2,
const double height)
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;
375 const double hs,
const double hs2)
377 const double logdepth = log(xVertGround/(depth * cos(theta)));
378 const double val = hs / (1/logdepth + hs2) + 1400.;
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);
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);
411 const double corr_sth = 0 * sst;
412 const double corr_lgE = (-0.0462 + (lgE-19.) * 0.08688);
413 return -(corr_sth+corr_lgE);
417 const double corr_sth = 0 * sst;
419 const double corr_lgE = (-0.0462 + (lgE-19.) * 0.08688);
420 return -(corr_sth + corr_lgE);
427 GHGauss(
const double r,
const double psi,
const double theta,
const double Xmax,
const double lambda)
438 return realEvent ? (797.16 + (lgE-19.) * 47.1) : (797.16 + (lgE-19.) * 47.1);
449 return realEvent ? 0.26 + 0.0179 + (lgE - 19.) * -0.0176 : 0.0179 + (lgE - 19.) * -0.0176;
455 XmaxFromRmu(
const double Rmu,
const double lgE,
const bool realEvent)
459 const double Xmax0 =
avgXmax0(lgE, realEvent);
460 const double lnRmu0 =
avglnRmu0(lgE, realEvent);
461 const double lambda = 21.28 + (lgE - 19.) * -0.645;
462 const double beta = 0.065 + (lgE - 19.) * -0.011;
464 const double deltalnRmu = (log(Rmu) - lnRmu0) > 0.24 ? 0.24 : log(Rmu) - lnRmu0;
466 const double Xmax = Xmax0 - lambda / beta * (deltalnRmu);
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;
483 XmaxDistribution(
const double Rmu,
const double Xmax,
const double lgE,
const bool realEvent)
487 const double mXmax =
XmaxFromRmu(Rmu, lgE, realEvent);
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)));
505 double lnA(
const double Rmu,
const double Xmax,
const double lgE,
const bool realEvent)
507 const double lnRmu = Rmu<1 ? Rmu-1 : log(Rmu);
509 const double Xmax0 =
avgXmax0(lgE, realEvent);
510 const double lnRmu0 =
avglnRmu0(lgE, realEvent);
512 const double lambda = 21.28 + (lgE - 19.) * -0.645;
513 const double beta = 0.065 + (lgE - 19.) * -0.011;
520 const double avgSigmaXmax = 70 + 40/(log(56)*lambda) *(Xmax-Xmax0);
521 const double avgSigmaLnRmu = 0.138 - 0.05/(log(56)*beta) * (lnRmu-lnRmu0);
523 const double sigmaXmax = avgSigmaXmax < 20 ? 20 : avgSigmaXmax;
524 const double sigmaLnRmu = avgSigmaLnRmu < 0.08 ? 0.08 : avgSigmaLnRmu;
526 const double phi0 =
utl::Sqr(beta*sigmaXmax/sigmaLnRmu)/lambda;
528 const double lnA = (phi0 * (lnRmu - lnRmu0) - beta * (Xmax-Xmax0))/(beta*(lambda+phi0));
double TwoFormAzimuthCorrectionFactor(const double theta, const double psi, const double M1, const double M3)
constexpr double kSqrt2Pi
double StandardDeviation(const std::vector< double > &v, const double mean)
double BiasCorrectionXmax(const double sst, const double lgE, const bool realEvent)
double BiasCorrectionRmu(const double sst, const double lgE, const bool realEvent)
double SpikyFunction(const double x, const double sigma, const double delta, const std::vector< double > &nvec)
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)
double ErfcInv(const double x)
double XmaxFromRmu(const double Rmu, const double lgE, const bool realEvent)
double LogarithmOfTruncatedGaussianPDF(const double x, const double mu, const double sigma)
double avglnRmu0(const double lgE, const bool realEvent)
double LogNormalMuFromModeAndSigma(const double mode, const double sigma, const double ns_bin)
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)
double ActivationFunctionBackground(const double f, const double fThrsh, const double cThrsh, const double bg)
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)
double XmaxDistribution(const double Rmu, const double Xmax, const double lgE, const bool realEvent)
double LogNormalMuFromExpAndStd(const double exp_val, const double std, const double ns_bin)
double LogNormalSigmaFromExpAndStd(const double exp_val, const double std, const double ns_bin)
double Gaussian(const double x, const double mu, const double sigma)
a t3: algo, second, usecond and a vector of <t3stat>
double lnA(const double Rmu, const double Xmax, const double lgE, const bool realEvent)
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)
double SigmaXmaxFromRmu(const double Rmu, const double lgE, const bool realEvent)
double GHGauss(const double r, const double psi, const double theta, const double Xmax, const double lambda)
double LogNormalTimeQuantileFromMuAndSigma(const double mu, const double sigma, const double quantile, const double ns_bin)
double avgXmax0(const double lgE, const bool realEvent)
double LogNormalMuFromTimeQuantileAndSigma(const double timeQuantile, const double quantile, const double sigma, const double ns_bin)
double AnalyticHeightFromDepth(const double depth, const double theta, const double xVertGround, const double hs, const double hs2)
double Mean(const std::vector< double > &v)
double LogNormalPDF(const double x, const double mu, const double sigma)
double LogNormalCDF(const double x, const double mu, const double sigma)
double GammaPDF(const double x, const double mu, const double theta)
double StD(const std::vector< double > &v, const double mean)
double LogarithmOfRayleighPDF(const double x, const double sigma)
double AzimuthCorrectionFactor(const double theta, const double psi, const double r, const double r0, const double l)
double NKG(const double r, const double n, const double rG, const double s)
double SignalStartTime(const double r, const double psi, const double theta, const double DX, const double Xvert, const double hs, const double hs2)