|
| 1 | +// Exercises the core phase-space generator (PhaseSpace::Generate) that GooFit |
| 2 | +// uses to produce normalization/toy samples. Runs on whichever Thrust backend |
| 3 | +// the build selects. |
| 4 | +#include <catch2/catch.hpp> |
| 5 | + |
| 6 | +#include <cmath> |
| 7 | +#include <vector> |
| 8 | + |
| 9 | +#include <mcbooster/GContainers.h> |
| 10 | +#include <mcbooster/Generate.h> |
| 11 | +#include <mcbooster/Vector4R.h> |
| 12 | + |
| 13 | +using mcbooster::GReal_t; |
| 14 | +using mcbooster::PhaseSpace; |
| 15 | +using mcbooster::Vector4R; |
| 16 | + |
| 17 | +namespace { |
| 18 | +// D0 -> K- pi+ pi0, a realistic 3-body decay. |
| 19 | +const GReal_t kMotherMass = 1.86484; |
| 20 | +const std::vector<GReal_t> kMasses{0.493677, 0.139570, 0.134977}; |
| 21 | +const mcbooster::GLong_t kNEvents = 5000; |
| 22 | +} // namespace |
| 23 | + |
| 24 | +TEST_CASE("PhaseSpace reports its configuration", "[generate]") { |
| 25 | + PhaseSpace phsp(kMotherMass, kMasses, kNEvents); |
| 26 | + CHECK(phsp.GetNDaughters() == 3); |
| 27 | + CHECK(phsp.GetNEvents() == kNEvents); |
| 28 | +} |
| 29 | + |
| 30 | +TEST_CASE("Generated daughters have the requested masses", "[generate]") { |
| 31 | + PhaseSpace phsp(kMotherMass, kMasses, kNEvents); |
| 32 | + phsp.Generate(Vector4R(kMotherMass, 0.0, 0.0, 0.0)); |
| 33 | + |
| 34 | + for(int d = 0; d < 3; d++) { |
| 35 | + mcbooster::Particles_h daughters = phsp.GetDaughters(d); |
| 36 | + for(int i = 0; i < 50; i++) { |
| 37 | + INFO("daughter " << d << ", event " << i); |
| 38 | + CHECK(daughters[i].mass() == Approx(kMasses[d]).epsilon(1e-5)); |
| 39 | + } |
| 40 | + } |
| 41 | +} |
| 42 | + |
| 43 | +TEST_CASE("Generation conserves energy and momentum", "[generate]") { |
| 44 | + PhaseSpace phsp(kMotherMass, kMasses, kNEvents); |
| 45 | + phsp.Generate(Vector4R(kMotherMass, 0.0, 0.0, 0.0)); |
| 46 | + |
| 47 | + mcbooster::Particles_h d0 = phsp.GetDaughters(0); |
| 48 | + mcbooster::Particles_h d1 = phsp.GetDaughters(1); |
| 49 | + mcbooster::Particles_h d2 = phsp.GetDaughters(2); |
| 50 | + |
| 51 | + for(int i = 0; i < 50; i++) { |
| 52 | + Vector4R total = d0[i] + d1[i] + d2[i]; |
| 53 | + INFO("event " << i); |
| 54 | + CHECK(total.get(0) == Approx(kMotherMass).epsilon(1e-5)); // energy -> mother mass |
| 55 | + CHECK(total.get(1) == Approx(0.0).margin(1e-5)); // px |
| 56 | + CHECK(total.get(2) == Approx(0.0).margin(1e-5)); // py |
| 57 | + CHECK(total.get(3) == Approx(0.0).margin(1e-5)); // pz |
| 58 | + CHECK(total.mass() == Approx(kMotherMass).epsilon(1e-5)); |
| 59 | + } |
| 60 | +} |
| 61 | + |
| 62 | +TEST_CASE("Weights are finite and positive", "[generate]") { |
| 63 | + PhaseSpace phsp(kMotherMass, kMasses, kNEvents); |
| 64 | + phsp.Generate(Vector4R(kMotherMass, 0.0, 0.0, 0.0)); |
| 65 | + |
| 66 | + mcbooster::RealVector_h weights = phsp.GetWeights(); |
| 67 | + REQUIRE(weights.size() == static_cast<size_t>(kNEvents)); |
| 68 | + |
| 69 | + GReal_t sum = 0.0; |
| 70 | + for(size_t i = 0; i < weights.size(); i++) { |
| 71 | + INFO("event " << i); |
| 72 | + REQUIRE(std::isfinite(weights[i])); |
| 73 | + REQUIRE(weights[i] > 0.0); |
| 74 | + sum += weights[i]; |
| 75 | + } |
| 76 | + CHECK(sum > 0.0); |
| 77 | +} |
| 78 | + |
| 79 | +TEST_CASE("Unweighting keeps a subset of events", "[generate]") { |
| 80 | + PhaseSpace phsp(kMotherMass, kMasses, kNEvents); |
| 81 | + phsp.Generate(Vector4R(kMotherMass, 0.0, 0.0, 0.0)); |
| 82 | + |
| 83 | + mcbooster::GULong_t kept = phsp.Unweight(); |
| 84 | + CHECK(kept > 0); |
| 85 | + CHECK(kept <= static_cast<mcbooster::GULong_t>(kNEvents)); |
| 86 | +} |
0 commit comments