Merge branch...

Merge branch '463-urqmd-generate-zero-kinematic-energy-secondries-and-induce-wrong-energy-cut-counting' into 'master'

Book kinetic energy in ParticleCut

Closes #463

See merge request AirShowerPhysics/corsika!403

corsika/detail/modules/ParticleCut.inl

@@ -133,26 +133,26 @@ namespace corsika {
     CORSIKA_LOG_DEBUG("p={}", particle.asString());
     if (doCutEm_ && is_em(pid)) {
       CORSIKA_LOG_DEBUG("removing em. particle...");
-      energy_emcut_ += energy;
+      energy_emcut_ += kine_energy;
       em_count_ += 1;
-      energy_event_ += energy;
+      energy_event_ += kine_energy;
       return true;
     } else if (doCutInv_ && is_neutrino(pid)) {
       CORSIKA_LOG_DEBUG("removing inv. particle...");
-      energy_invcut_ += energy;
+      energy_invcut_ += kine_energy;
       inv_count_ += 1;
-      energy_event_ += energy;
+      energy_event_ += kine_energy;
       return true;
     } else if (isBelowEnergyCut(particle)) {
       CORSIKA_LOG_DEBUG("removing low en. particle...");
-      energy_cut_ += energy;
+      energy_cut_ += kine_energy;
       energy_count_ += 1;
-      energy_event_ += energy;
+      energy_event_ += kine_energy;
       return true;
     } else if (particle.getTime() > 10_ms) {
       CORSIKA_LOG_DEBUG("removing OLD particle...");
-      energy_timecut_ += energy;
-      energy_event_ += energy;
+      energy_timecut_ += kine_energy;
+      energy_event_ += kine_energy;
       return true;
     }
     return false; // this particle will not be removed/cut
@@ -192,10 +192,10 @@ namespace corsika {
   inline void ParticleCut::showResults() {
     CORSIKA_LOG_INFO(
         "\n ******************************\n "
-        " energy removed by cut of electromagnetic (GeV): {} (number: {})\n "
-        " energy removed by cut of invisible (GeV): {} (number: {})\n "
-        " energy removed by kinetic energy cut (GeV): {} (number: {}) \n "
-        " energy removed by time cut (GeV): {} \n"
+        " kinetic energy removed by cut of electromagnetic (GeV): {} (number: {})\n "
+        " kinetic energy removed by cut of invisible (GeV): {} (number: {})\n "
+        " kinetic energy removed by kinetic energy cut (GeV): {} (number: {}) \n "
+        " kinetic energy removed by time cut (GeV): {} \n"
         " ******************************",
         energy_emcut_ / 1_GeV, em_count_, energy_invcut_ / 1_GeV, inv_count_,
         energy_cut_ / 1_GeV, energy_count_, energy_timecut_ / 1_GeV);

tests/modules/testParticleCut.cpp

@@ -76,7 +76,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
 
     CHECK(view.getEntries() == 9);
     CHECK(cut.getNumberInvParticles() == 2);
-    CHECK(cut.getInvEnergy() / 1_GeV == 2000);
+    CHECK(cut.getInvEnergy() / 1_GeV == Approx(2000.));
   }
 
   SECTION("cut on particle type: em") {
@@ -101,16 +101,15 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
 
     CHECK(view.getEntries() == 10);
     CHECK(cut.getNumberEmParticles() == 1);
-    CHECK(cut.getEmEnergy() / 1_GeV == Approx(1000. + 0.000511));
+    CHECK(cut.getEmEnergy() / 1_GeV == 1000.);
   }
 
   SECTION("cut low energy") {
     ParticleCut cut(20_GeV, true, true);
 
     // add primary particle to stack
-    auto particle = stack.addParticle(std::make_tuple(Code::Proton, Eabove - Proton::mass,
-                                                      DirectionVector(rootCS, {1, 0, 0}),
-                                                      point0, 0_ns));
+    auto particle = stack.addParticle(std::make_tuple(
+        Code::Proton, Eabove, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
     // view on secondary particles
     setup::StackView view(particle);
     // ref. to primary particle through the secondary view.
@@ -151,12 +150,10 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
     // add secondaries
     projectile.addSecondary(std::make_tuple(
         Code::Photon, 3_MeV, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
-    projectile.addSecondary(std::make_tuple(Code::Electron, 3_MeV - Electron::mass,
-                                            DirectionVector(rootCS, {1, 0, 0}), point0,
-                                            0_ns));
-    projectile.addSecondary(std::make_tuple(Code::PiPlus, 4_GeV - PiPlus::mass,
-                                            DirectionVector(rootCS, {1, 0, 0}), point0,
-                                            0_ns));
+    projectile.addSecondary(std::make_tuple(
+        Code::Electron, 3_MeV, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
+    projectile.addSecondary(std::make_tuple(
+        Code::PiPlus, 4_GeV, DirectionVector(rootCS, {1, 0, 0}), point0, 0_ns));
 
     unsigned short A = 18;
     unsigned short Z = 8;
@@ -197,9 +194,8 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") {
     // add secondaries, all with energies above the threshold
     // only cut is by time
     for (auto proType : particleList) {
-      projectile.addSecondary(std::make_tuple(proType, Eabove - get_mass(proType),
-                                              DirectionVector(rootCS, {1, 0, 0}), point0,
-                                              too_late));
+      projectile.addSecondary(std::make_tuple(
+          proType, Eabove, DirectionVector(rootCS, {1, 0, 0}), point0, too_late));
     }
     cut.doSecondaries(view);