Merge pull request #363 from tudo-astroparticlephysics/323-no-interac…

…tiontypedecay-type-added-to-secondaries-if-particle-decays-at-rest-mass

Add order of `Type` enum in Propagator

src/PROPOSAL/PROPOSAL/Propagator.h

@@ -72,8 +72,8 @@ class Propagator {
 
     ParticleDef p_def;
     enum Type : int {
-        MinimalE = 0,
-        Decay = 1,
+        Decay = 0,
+        MinimalE = 1,
         Stochastic = 2,
     };
     enum AdvancementType : int {

tests/Propagator_TEST.cxx

@@ -4,6 +4,7 @@
 #include "PROPOSAL/propagation_utility/TimeBuilder.h"
 #include "PROPOSAL/propagation_utility/InteractionBuilder.h"
 #include "PROPOSAL/propagation_utility/ContRandBuilder.h"
+#include "PROPOSAL/propagation_utility/DecayBuilder.h"
 #include "PROPOSAL/density_distr/density_homogeneous.h"
 #include "PROPOSAL/scattering/ScatteringFactory.h"
 #include "PROPOSAL/geometry/Sphere.h"
@@ -57,6 +58,46 @@ TEST(Propagator, min_energy)
     }
 }
 
+TEST(Propagator, RestMassDecay)
+{
+    // make sure that all muons decay at the end of their propagation, even if they have reached their rest mass
+    // (see issue https://github.com/tudo-astroparticlephysics/PROPOSAL/issues/323)
+    auto p_def = MuMinusDef();
+    auto medium = Ice();
+    auto cuts = std::make_shared<EnergyCutSettings>(INF, 0.05, true);
+    auto cross = GetStdCrossSections(p_def, medium, cuts, true);
+
+    auto collection = PropagationUtility::Collection();
+    collection.interaction_calc = make_interaction(cross, true);
+    collection.displacement_calc = make_displacement(cross, true);
+    collection.time_calc = make_time(cross, p_def, true);
+    collection.cont_rand = make_contrand(cross, true);
+    collection.scattering = make_scattering(MultipleScatteringType::Highland, {}, p_def, medium);
+    collection.decay_calc = make_decay(cross, p_def, true);
+
+    auto prop_utility = PropagationUtility(collection);
+
+    auto density_distr = std::make_shared<Density_homogeneous>(medium);
+    auto world = std::make_shared<Sphere>(Cartesian3D(0, 0, 0), 1e20);
+
+    auto sector = std::make_tuple(world, prop_utility, density_distr);
+    std::vector<Sector> sec_vec = {sector};
+
+    auto prop = Propagator(p_def, sec_vec);
+
+    auto init_state = ParticleState();
+    init_state.energy = 150;
+    init_state.position = Cartesian3D(0, 0, 0);
+    init_state.direction = Cartesian3D(0, 0, 1);
+
+    for (size_t i=0; i<1000; i++) {
+        auto sec = prop.Propagate(init_state);
+        EXPECT_FALSE(sec.GetDecayProducts().empty()); // check that decay products are produced
+    }
+
+}
+
+
 int main(int argc, char** argv)
 {
     ::testing::InitGoogleTest(&argc, argv);