Skip to content

Commit

Permalink
Merge pull request #242 from tudo-astroparticlephysics/UnitTest_overhaul
Browse files Browse the repository at this point in the history
UnitTest overhaul
  • Loading branch information
sudojan authored Jan 29, 2022
2 parents f38a4f2 + ac114c7 commit 8307c46
Show file tree
Hide file tree
Showing 22 changed files with 667 additions and 2,481 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/cpp.yml
Original file line number Diff line number Diff line change
Expand Up @@ -112,4 +112,4 @@ jobs:
path: PROPOSAL_BUILD
key: ${{ runner.os }}-cache-build-${{ matrix.compiler }}-${{ github.sha }}-key
- name: Run tests
run: ctest --test-dir PROPOSAL_BUILD -j2 --verbose -E "(Brems|Photo|Epair|Mupair)"
run: ctest --test-dir PROPOSAL_BUILD -j2 --verbose -E "(Brems.*Interpolant|Photo.*Interpolant|Epair.*Interpolant|Mupair.*Interpolant)"
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ crosssection::HardComponent::HardComponent(const ParticleDef& particle_def)
} else {
Logging::Get("proposal.parametrization")
->error(
"No HardComponent tables provided for the given particle %s",
"No HardComponent tables provided for the given particle {}",
particle_def.name.c_str());
}
}
Expand Down
10 changes: 7 additions & 3 deletions src/pyPROPOSAL/detail/secondaries.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -135,9 +135,13 @@ void init_secondaries(py::module& m)
secondaries::PhotoPairProduction> {}
.decl_rho_param(m_sub, "PhotoKochMotzForwardPeaked");

SecondariesBuilder<secondaries::WeakCooperSarkarMertsch,
secondaries::WeakInteraction> {}
.decl_param(m_sub, "WeakCooperSarkarMertsch");
py::class_<secondaries::WeakCooperSarkarMertsch, secondaries::WeakInteraction,
std::shared_ptr<secondaries::WeakCooperSarkarMertsch>>(m_sub,
"WeakCooperSarkarMertsch")
.def(py::init<ParticleDef, Medium>(), py::arg("particle_def"),
py::arg("medium"))
.def("calculate_relative_loss", &secondaries::WeakCooperSarkarMertsch::CalculateRelativeLoss,
py::arg("energy"), py::arg("rnd"), py::arg("component"));

py::class_<SecondariesCalculator, std::shared_ptr<SecondariesCalculator>>(
m_sub, "SecondariesCalculator")
Expand Down
183 changes: 22 additions & 161 deletions tests/Annihilation_TEST.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -31,98 +31,6 @@ ParticleDef getParticleDef(const std::string& name)
}
}

// TEST(Comparison, Comparison_equal_particle)
// {
// ParticleDef particle_def = EPlusDef();
// auto medium = std::make_shared<const Water>();
// double multiplier = 1.;

// Annihilation* Anni_A = new AnnihilationHeitler(particle_def, medium,
// multiplier); Annihilation* Anni_B = new AnnihilationHeitler(particle_def,
// medium, multiplier); EXPECT_TRUE(*Anni_A == *Anni_B);

// AnnihilationHeitler param_int(particle_def, medium, multiplier);
// EXPECT_TRUE(param_int == *Anni_A);

// AnnihilationIntegral* Int_A = new AnnihilationIntegral(param_int);
// CrossSectionIntegral* Int_B = new AnnihilationIntegral(param_int);
// EXPECT_TRUE(*Int_A == *Int_B);

// InterpolationDef InterpolDef;

// AnnihilationInterpolant* Interpol_A = new
// AnnihilationInterpolant(param_int, InterpolDef); CrossSectionInterpolant*
// Interpol_B = new AnnihilationInterpolant(param_int, InterpolDef);
// EXPECT_TRUE(*Interpol_A == *Interpol_B);

// delete Anni_A;
// delete Anni_B;
// delete Int_A;
// delete Int_B;
// delete Interpol_A;
// delete Interpol_B;
// }

// TEST(Comparison, Comparison_not_equal)
// {
// ParticleDef mu_def = MuMinusDef();
// ParticleDef e_minus_def = EMinusDef();
// ParticleDef e_plus_def = EPlusDef();
// auto medium_1 = std::make_shared<const Water>();
// auto medium_2 = std::make_shared<const Ice>();
// double multiplier_1 = 1.;
// double multiplier_2 = 2.;

// AnnihilationHeitler Anni_A(mu_def, medium_1, multiplier_1);
// AnnihilationHeitler Anni_B(e_minus_def, medium_1, multiplier_1);
// AnnihilationHeitler Anni_C(mu_def, medium_2, multiplier_1);
// AnnihilationHeitler Anni_D(e_plus_def, medium_1, multiplier_1);
// AnnihilationHeitler Anni_E(mu_def, medium_1, multiplier_2);
// EXPECT_TRUE(Anni_A != Anni_B);
// EXPECT_TRUE(Anni_A != Anni_C);
// EXPECT_TRUE(Anni_A != Anni_D);
// EXPECT_TRUE(Anni_A != Anni_E);

// AnnihilationIntegral Int_A(Anni_A);
// AnnihilationIntegral Int_B(Anni_B);
// EXPECT_TRUE(Int_A != Int_B);

// InterpolationDef InterpolDef;
// AnnihilationIntegral AnniIntegral_A(Anni_A);
// AnnihilationIntegral AnniIntegral_B(Anni_B);
// AnnihilationIntegral AnniIntegral_C(Anni_C);
// AnnihilationIntegral AnniIntegral_D(Anni_D);
// AnnihilationIntegral AnniIntegral_E(Anni_E);
// EXPECT_TRUE(AnniIntegral_A != AnniIntegral_B);
// EXPECT_TRUE(AnniIntegral_A != AnniIntegral_C);
// EXPECT_TRUE(AnniIntegral_A != AnniIntegral_D);
// EXPECT_TRUE(AnniIntegral_A != AnniIntegral_E);

// AnnihilationIntegral Integral_A(AnniIntegral_A);
// AnnihilationIntegral Integral_B(AnniIntegral_B);
// EXPECT_TRUE(Integral_A != Integral_B);
// }

// TEST(Assignment, Copyconstructor)
// {
// ParticleDef particle_def = EPlusDef();
// auto medium = std::make_shared<const Water>();
// double multiplier = 1.;

// AnnihilationHeitler Anni_A(particle_def, medium, multiplier);
// AnnihilationHeitler Anni_B = Anni_A;
// EXPECT_TRUE(Anni_A == Anni_B);

// AnnihilationIntegral Int_A(Anni_A);
// AnnihilationIntegral Int_B = Int_A;
// EXPECT_TRUE(Int_A == Int_B);

// InterpolationDef InterpolDef;
// AnnihilationInterpolant AnniInterpol_A(Anni_A, InterpolDef);
// AnnihilationInterpolant AnniInterpol_B = AnniInterpol_A;
// EXPECT_TRUE(AnniInterpol_A == AnniInterpol_B);
// }

TEST(Annihilation, Test_of_dNdx)
{
std::ifstream in;
Expand All @@ -138,27 +46,16 @@ TEST(Annihilation, Test_of_dNdx)

while (in >> particleName >> mediumName >> multiplier >> energy
>> parametrization >> dNdx_stored) {

ParticleDef particle_def = getParticleDef(particleName);
if (mediumName == "ice" || mediumName == "water")
mediumName += "PDG2001";
auto medium = CreateMedium(mediumName);

nlohmann::json config;
parametrization.erase(0, 12);
config["parametrization"] = parametrization;

auto cross = make_annihilation(particle_def, *medium, false, config);

if (energy <= particle_def.mass) {
#ifndef NDEBUG
EXPECT_DEATH(cross->CalculatedNdx(energy), "");
#endif
continue;
}
dNdx_new = cross->CalculatedNdx(energy) * medium->GetMassDensity();
dNdx_new *= medium->GetComponents()
.size(); // This has (probably) been a mistake in
// the old PROPOSAL version
dNdx_new = cross->CalculatedNdx(energy);
EXPECT_NEAR(dNdx_new, dNdx_stored, 1e-10 * dNdx_stored);
}
}
Expand All @@ -176,7 +73,6 @@ TEST(Annihilation, Test_Stochastic_Loss)
double rnd1;
double rnd2;
double stochastic_loss_stored;
double stochastic_loss_new;

std::cout.precision(16);
RandomGenerator::Get().SetSeed(0);
Expand All @@ -185,39 +81,21 @@ TEST(Annihilation, Test_Stochastic_Loss)
>> parametrization >> rnd1 >> rnd2 >> stochastic_loss_stored) {

ParticleDef particle_def = getParticleDef(particleName);
if (mediumName == "ice" || mediumName == "water")
mediumName += "PDG2001";
auto medium = CreateMedium(mediumName);

nlohmann::json config;
parametrization.erase(0, 12);
config["parametrization"] = parametrization;

auto cross = make_annihilation(particle_def, *medium, false, config);

if (energy <= particle_def.mass) {
#ifndef NDEBUG
EXPECT_DEATH(cross->CalculatedNdx(energy), "");
#endif
continue;
}
auto dNdx_full = cross->CalculatedNdx(energy);
auto components = medium->GetComponents();
double sum = 0;

for (auto comp : components) {
double dNdx_for_comp = cross->CalculatedNdx(energy, comp.GetHash());
sum += dNdx_for_comp;
if (sum > dNdx_full * rnd2) {
double rate_new = dNdx_for_comp * rnd1;
stochastic_loss_new = energy
* cross->CalculateStochasticLoss(
comp.GetHash(), energy, rate_new);
EXPECT_NEAR(stochastic_loss_new, stochastic_loss_stored,
1E-6 * stochastic_loss_stored);
break;
}
}
auto comp = components.at(int(rnd2 * components.size()));

double dNdx_for_comp = cross->CalculatedNdx(energy, comp.GetHash());

auto stochastic_loss = cross->CalculateStochasticLoss(comp.GetHash(), energy, rnd1 * dNdx_for_comp);
EXPECT_NEAR(stochastic_loss, stochastic_loss_stored, 1e-6 * stochastic_loss_stored);
EXPECT_EQ(stochastic_loss, 1.); // we always loose all energy
}
}

Expand All @@ -236,26 +114,21 @@ TEST(Annihilation, Test_of_dNdx_Interpolant)

while (in >> particleName >> mediumName >> multiplier >> energy
>> parametrization >> dNdx_stored) {

ParticleDef particle_def = getParticleDef(particleName);
if (mediumName == "ice" || mediumName == "water")
mediumName += "PDG2001";
auto medium = CreateMedium(mediumName);

nlohmann::json config;
parametrization.erase(0, 12);
config["parametrization"] = parametrization;

auto cross = make_annihilation(particle_def, *medium, true, config);

dNdx_new = cross->CalculatedNdx(energy) * medium->GetMassDensity();
dNdx_new *= medium->GetComponents()
.size(); // This has (probably) been a mistake in
// the old PROPOSAL version
EXPECT_NEAR(dNdx_new, dNdx_stored, 1e-4 * dNdx_stored); // 1e-10 -> 1e-4
dNdx_new = cross->CalculatedNdx(energy);
EXPECT_NEAR(dNdx_new, dNdx_stored, 1e-4 * dNdx_stored);
}
}

TEST(Annihilation, Test_of_e_interpol)
TEST(Annihilation, Test_of_e_Interpolant)
{
std::ifstream in;
getTestFile("Anni_e.txt", in);
Expand All @@ -268,40 +141,28 @@ TEST(Annihilation, Test_of_e_interpol)
double rnd1;
double rnd2;
double stochastic_loss_stored;
double stochastic_loss_new;

RandomGenerator::Get().SetSeed(0);

while (in >> particleName >> mediumName >> multiplier >> energy
>> parametrization >> rnd1 >> rnd2 >> stochastic_loss_stored) {

ParticleDef particle_def = getParticleDef(particleName);
if (mediumName == "ice" || mediumName == "water")
mediumName += "PDG2001";
auto medium = CreateMedium(mediumName);

nlohmann::json config;
parametrization.erase(0, 12);
config["parametrization"] = parametrization;

auto cross = make_annihilation(particle_def, *medium, true, config);
auto cross = make_annihilation(particle_def, *medium, false, config);

auto dNdx_full = cross->CalculatedNdx(energy);
auto components = medium->GetComponents();
double sum = 0;

for (auto comp : components) {
double dNdx_for_comp = cross->CalculatedNdx(energy, comp.GetHash());
sum += dNdx_for_comp;
if (sum > dNdx_full * rnd2) {
double rate_new = dNdx_for_comp * rnd1;
stochastic_loss_new = energy
* cross->CalculateStochasticLoss(
comp.GetHash(), energy, rate_new);
EXPECT_NEAR(stochastic_loss_new, stochastic_loss_stored,
1E-6 * stochastic_loss_stored);
break;
}
}
auto comp = components.at(int(rnd2 * components.size()));

double dNdx_for_comp = cross->CalculatedNdx(energy, comp.GetHash());

auto stochastic_loss = cross->CalculateStochasticLoss(comp.GetHash(), energy, rnd1 * dNdx_for_comp);
EXPECT_NEAR(stochastic_loss, stochastic_loss_stored, 1e-6 * stochastic_loss_stored);
EXPECT_EQ(stochastic_loss, 1.); // we always loose all energy
}
}

Expand Down
Loading

0 comments on commit 8307c46

Please sign in to comment.