74 lines
2.2 KiB
C++
74 lines
2.2 KiB
C++
#include "yields.hh"
|
|
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <iomanip>
|
|
#include <iostream>
|
|
|
|
namespace {
|
|
|
|
constexpr double muonMassKeV = 105.658369*1000.;
|
|
constexpr double tolerance = 1.e-12;
|
|
|
|
struct YieldCase {
|
|
double energyKeV;
|
|
std::array<double,3> expected;
|
|
};
|
|
|
|
bool nearlyEqual(double actual, double expected)
|
|
{
|
|
return std::abs(actual-expected)<=tolerance;
|
|
}
|
|
|
|
} // namespace
|
|
|
|
int main()
|
|
{
|
|
// Reference values characterize the Gonin charge-state model at points
|
|
// spanning its documented 0.4--120 keV range. Components are mu+, Mu,
|
|
// and mu-, in the order returned by Yields::GetYields().
|
|
constexpr std::array<YieldCase,7> cases = {{
|
|
{0.4, {0.10683927921255008, 0.85160645015993075, 0.041554270627519131}},
|
|
{1.0, {0.19224098891496452, 0.77532099475070160, 0.032438016334333922}},
|
|
{5.0, {0.50213564424452095, 0.49076257493468323, 0.0071017808207957753}},
|
|
{10.0, {0.73317729701811696, 0.26571062487344432, 0.0011120781084387650}},
|
|
{30.0, {0.94261066802595883, 0.054679886047411919,0.0027094459266292200}},
|
|
{60.0, {0.95906983372572263, 0.029335370816802930,0.011594795457474402}},
|
|
{120.0, {0.97088138720749961, 0.0021919369643415010,0.026926675828158865}},
|
|
}};
|
|
|
|
Yields model;
|
|
bool passed = true;
|
|
|
|
for (const auto& testCase : cases) {
|
|
double actual[3] = {};
|
|
model.GetYields(testCase.energyKeV,muonMassKeV,actual);
|
|
|
|
double sum = 0.;
|
|
for (std::size_t i=0; i<testCase.expected.size(); ++i) {
|
|
sum += actual[i];
|
|
if (!std::isfinite(actual[i]) || actual[i]<0. || actual[i]>1.) {
|
|
std::cerr<<"non-physical yield at "<<testCase.energyKeV
|
|
<<" keV, component "<<i<<": "<<actual[i]<<'\n';
|
|
passed = false;
|
|
}
|
|
if (!nearlyEqual(actual[i],testCase.expected[i])) {
|
|
std::cerr<<std::setprecision(17)
|
|
<<"yield changed at "<<testCase.energyKeV
|
|
<<" keV, component "<<i<<": expected "
|
|
<<testCase.expected[i]<<", got "<<actual[i]<<'\n';
|
|
passed = false;
|
|
}
|
|
}
|
|
|
|
if (!nearlyEqual(sum,1.)) {
|
|
std::cerr<<std::setprecision(17)
|
|
<<"yields do not sum to one at "<<testCase.energyKeV
|
|
<<" keV: "<<sum<<'\n';
|
|
passed = false;
|
|
}
|
|
}
|
|
|
|
return passed ? 0 : 1;
|
|
}
|