improved some tests, improved EnergyLoss::FillProfile

Documentation/Examples/cascade_example.cc

@@ -57,7 +57,7 @@ using namespace corsika::units::si;
 //
 int main() {
 
-  logging::SetLevel(logging::level::info);
+  logging::SetLevel(logging::level::debug);
 
   std::cout << "cascade_example" << std::endl;
 

Framework/Analytics/testClassTimer.cc

@@ -110,7 +110,7 @@ TEST_CASE("Analytics", "[Timer]") {
 
     tc.call();
 
-    REQUIRE(tc.getTime().count() == Approx(100000).margin(1000));
+    CHECK(tc.getTime().count() == Approx(100000).margin(10000));
   }
 
   SECTION("Measure runtime of a function with arguments") {
@@ -120,7 +120,7 @@ TEST_CASE("Analytics", "[Timer]") {
 
     tc.call(1);
 
-    REQUIRE(tc.getTime().count() == Approx(100000).margin(1000));
+    CHECK(tc.getTime().count() == Approx(100000).margin(10000));
   }
 
   SECTION("Measure runtime of a const function without arguments") {
@@ -131,17 +131,17 @@ TEST_CASE("Analytics", "[Timer]") {
 
     tc.call();
 
-    REQUIRE(tc.getTime().count() == Approx(100000).margin(1000));
+    CHECK(tc.getTime().count() == Approx(100000).margin(10000));
   }
 
   SECTION("Measure runtime of function inside class") {
 
     auto test = foo();
-    REQUIRE(test.inside() == 123);
+    CHECK(test.inside() == 123);
   }
 
   SECTION("Measure runtime of function inside class") {
     auto test = fooT3<fooT1>();
-    REQUIRE(test.inside_t(1, 'a', 'b') == 123);
+    CHECK(test.inside_t(1, 'a', 'b') == 123);
   }
 }

Framework/Particles/testParticles.cc

@@ -18,123 +18,126 @@ using namespace corsika::particles;
 TEST_CASE("ParticleProperties", "[Particles]") {
 
   SECTION("Types") {
-    REQUIRE(Electron::GetCode() == Code::Electron);
-    REQUIRE(Positron::GetCode() == Code::Positron);
-    REQUIRE(Proton::GetCode() == Code::Proton);
-    REQUIRE(Neutron::GetCode() == Code::Neutron);
-    REQUIRE(Gamma::GetCode() == Code::Gamma);
-    REQUIRE(PiPlus::GetCode() == Code::PiPlus);
+    CHECK(Electron::GetCode() == Code::Electron);
+    CHECK(Positron::GetCode() == Code::Positron);
+    CHECK(Proton::GetCode() == Code::Proton);
+    CHECK(Neutron::GetCode() == Code::Neutron);
+    CHECK(Gamma::GetCode() == Code::Gamma);
+    CHECK(PiPlus::GetCode() == Code::PiPlus);
   }
 
   SECTION("Masses") {
-    REQUIRE(Electron::GetMass() / (511_keV) == Approx(1));
-    REQUIRE(Electron::GetMass() / GetMass(Code::Electron) == Approx(1));
+    CHECK(Electron::GetMass() / (511_keV) == Approx(1));
+    CHECK(Electron::GetMass() / GetMass(Code::Electron) == Approx(1));
 
-    REQUIRE((Proton::GetMass() + Neutron::GetMass()) /
-                corsika::units::constants::nucleonMass ==
-            Approx(2));
+    CHECK((Proton::GetMass() + Neutron::GetMass()) /
+              corsika::units::constants::nucleonMass ==
+          Approx(2));
   }
 
   SECTION("Charges") {
-    REQUIRE(Electron::GetCharge() / constants::e == Approx(-1));
-    REQUIRE(Positron::GetCharge() / constants::e == Approx(+1));
-    REQUIRE(GetCharge(Positron::GetAntiParticle()) / constants::e == Approx(-1));
+    CHECK(Electron::GetCharge() / constants::e == Approx(-1));
+    CHECK(Positron::GetCharge() / constants::e == Approx(+1));
+    CHECK(GetCharge(Positron::GetAntiParticle()) / constants::e == Approx(-1));
   }
 
   SECTION("Names") {
-    REQUIRE(Electron::GetName() == "e-");
-    REQUIRE(PiMinus::GetName() == "pi-");
-    REQUIRE(Nucleus::GetName() == "nucleus");
-    REQUIRE(Iron::GetName() == "iron");
+    CHECK(Electron::GetName() == "e-");
+    CHECK(PiMinus::GetName() == "pi-");
+    CHECK(Nucleus::GetName() == "nucleus");
+    CHECK(Iron::GetName() == "iron");
   }
 
   SECTION("PDG") {
-    REQUIRE(GetPDG(Code::PiPlus) == PDGCode::PiPlus);
-    REQUIRE(GetPDG(Code::DPlus) == PDGCode::DPlus);
-    REQUIRE(GetPDG(Code::NuMu) == PDGCode::NuMu);
-    REQUIRE(GetPDG(Code::NuE) == PDGCode::NuE);
-    REQUIRE(GetPDG(Code::MuMinus) == PDGCode::MuMinus);
-
-    REQUIRE(static_cast<int>(GetPDG(Code::PiPlus)) == 211);
-    REQUIRE(static_cast<int>(GetPDG(Code::DPlus)) == 411);
-    REQUIRE(static_cast<int>(GetPDG(Code::NuMu)) == 14);
-    REQUIRE(static_cast<int>(GetPDG(Code::NuEBar)) == -12);
-    REQUIRE(static_cast<int>(GetPDG(Code::MuMinus)) == 13);
+    CHECK(GetPDG(Code::PiPlus) == PDGCode::PiPlus);
+    CHECK(GetPDG(Code::DPlus) == PDGCode::DPlus);
+    CHECK(GetPDG(Code::NuMu) == PDGCode::NuMu);
+    CHECK(GetPDG(Code::NuE) == PDGCode::NuE);
+    CHECK(GetPDG(Code::MuMinus) == PDGCode::MuMinus);
+
+    CHECK(static_cast<int>(GetPDG(Code::PiPlus)) == 211);
+    CHECK(static_cast<int>(GetPDG(Code::DPlus)) == 411);
+    CHECK(static_cast<int>(GetPDG(Code::NuMu)) == 14);
+    CHECK(static_cast<int>(GetPDG(Code::NuEBar)) == -12);
+    CHECK(static_cast<int>(GetPDG(Code::MuMinus)) == 13);
   }
 
   SECTION("Conversion PDG -> internal") {
-    REQUIRE(ConvertFromPDG(PDGCode::KStarMinus) == Code::KStarMinus);
-    REQUIRE(ConvertFromPDG(PDGCode::MuPlus) == Code::MuPlus);
-    REQUIRE(ConvertFromPDG(PDGCode::SigmaStarCMinusBar) == Code::SigmaStarCMinusBar);
+    CHECK(ConvertFromPDG(PDGCode::KStarMinus) == Code::KStarMinus);
+    CHECK(ConvertFromPDG(PDGCode::MuPlus) == Code::MuPlus);
+    CHECK(ConvertFromPDG(PDGCode::SigmaStarCMinusBar) == Code::SigmaStarCMinusBar);
   }
 
   SECTION("Lifetimes") {
-    REQUIRE(GetLifetime(Code::Electron) ==
-            std::numeric_limits<double>::infinity() * corsika::units::si::second);
-    REQUIRE(GetLifetime(Code::DPlus) < GetLifetime(Code::Gamma));
-    REQUIRE(GetLifetime(Code::RhoPlus) / corsika::units::si::second ==
-            (Approx(4.414566727909413e-24).epsilon(1e-3)));
-    REQUIRE(GetLifetime(Code::SigmaMinusBar) / corsika::units::si::second ==
-            (Approx(8.018880848563575e-11).epsilon(1e-5)));
-    REQUIRE(GetLifetime(Code::MuPlus) / corsika::units::si::second ==
-            (Approx(2.1970332555864364e-06).epsilon(1e-5)));
+    CHECK(GetLifetime(Code::Electron) ==
+          std::numeric_limits<double>::infinity() * corsika::units::si::second);
+    CHECK(GetLifetime(Code::DPlus) < GetLifetime(Code::Gamma));
+    CHECK(GetLifetime(Code::RhoPlus) / corsika::units::si::second ==
+          (Approx(4.414566727909413e-24).epsilon(1e-3)));
+    CHECK(GetLifetime(Code::SigmaMinusBar) / corsika::units::si::second ==
+          (Approx(8.018880848563575e-11).epsilon(1e-5)));
+    CHECK(GetLifetime(Code::MuPlus) / corsika::units::si::second ==
+          (Approx(2.1970332555864364e-06).epsilon(1e-5)));
   }
 
   SECTION("Particle groups: electromagnetic") {
-    REQUIRE(IsEM(Code::Gamma));
-    REQUIRE(IsEM(Code::Electron));
-    REQUIRE_FALSE(IsEM(Code::MuPlus));
-    REQUIRE_FALSE(IsEM(Code::NuE));
-    REQUIRE_FALSE(IsEM(Code::Proton));
-    REQUIRE_FALSE(IsEM(Code::PiPlus));
-    REQUIRE_FALSE(IsEM(Code::Oxygen));
+    CHECK(IsEM(Code::Gamma));
+    CHECK(IsEM(Code::Electron));
+    CHECK_FALSE(IsEM(Code::MuPlus));
+    CHECK_FALSE(IsEM(Code::NuE));
+    CHECK_FALSE(IsEM(Code::Proton));
+    CHECK_FALSE(IsEM(Code::PiPlus));
+    CHECK_FALSE(IsEM(Code::Oxygen));
   }
 
   SECTION("Particle groups: hadrons") {
-    REQUIRE_FALSE(IsHadron(Code::Gamma));
-    REQUIRE_FALSE(IsHadron(Code::Electron));
-    REQUIRE_FALSE(IsHadron(Code::MuPlus));
-    REQUIRE_FALSE(IsHadron(Code::NuE));
-    REQUIRE(IsHadron(Code::Proton));
-    REQUIRE(IsHadron(Code::PiPlus));
-    REQUIRE(IsHadron(Code::Oxygen));
-    REQUIRE(IsHadron(Code::Nucleus));
+    CHECK_FALSE(IsHadron(Code::Gamma));
+    CHECK_FALSE(IsHadron(Code::Electron));
+    CHECK_FALSE(IsHadron(Code::MuPlus));
+    CHECK_FALSE(IsHadron(Code::NuE));
+    CHECK(IsHadron(Code::Proton));
+    CHECK(IsHadron(Code::PiPlus));
+    CHECK(IsHadron(Code::Oxygen));
+    CHECK(IsHadron(Code::Nucleus));
   }
 
   SECTION("Particle groups: muons") {
-    REQUIRE_FALSE(IsMuon(Code::Gamma));
-    REQUIRE_FALSE(IsMuon(Code::Electron));
-    REQUIRE(IsMuon(Code::MuPlus));
-    REQUIRE_FALSE(IsMuon(Code::NuE));
-    REQUIRE_FALSE(IsMuon(Code::Proton));
-    REQUIRE_FALSE(IsMuon(Code::PiPlus));
-    REQUIRE_FALSE(IsMuon(Code::Oxygen));
+    CHECK_FALSE(IsMuon(Code::Gamma));
+    CHECK_FALSE(IsMuon(Code::Electron));
+    CHECK(IsMuon(Code::MuPlus));
+    CHECK_FALSE(IsMuon(Code::NuE));
+    CHECK_FALSE(IsMuon(Code::Proton));
+    CHECK_FALSE(IsMuon(Code::PiPlus));
+    CHECK_FALSE(IsMuon(Code::Oxygen));
   }
 
   SECTION("Particle groups: neutrinos") {
-    REQUIRE_FALSE(IsNeutrino(Code::Gamma));
-    REQUIRE_FALSE(IsNeutrino(Code::Electron));
-    REQUIRE_FALSE(IsNeutrino(Code::MuPlus));
-    REQUIRE(IsNeutrino(Code::NuE));
-    REQUIRE_FALSE(IsNeutrino(Code::Proton));
-    REQUIRE_FALSE(IsNeutrino(Code::PiPlus));
-    REQUIRE_FALSE(IsNeutrino(Code::Oxygen));
+    CHECK_FALSE(IsNeutrino(Code::Gamma));
+    CHECK_FALSE(IsNeutrino(Code::Electron));
+    CHECK_FALSE(IsNeutrino(Code::MuPlus));
+    CHECK(IsNeutrino(Code::NuE));
+    CHECK_FALSE(IsNeutrino(Code::Proton));
+    CHECK_FALSE(IsNeutrino(Code::PiPlus));
+    CHECK_FALSE(IsNeutrino(Code::Oxygen));
   }
 
   SECTION("Nuclei") {
-    REQUIRE_FALSE(IsNucleus(Code::Gamma));
-    REQUIRE(IsNucleus(Code::Argon));
-    REQUIRE_FALSE(IsNucleus(Code::Proton));
-    REQUIRE(IsNucleus(Code::Hydrogen));
-    REQUIRE(Argon::IsNucleus());
-    REQUIRE_FALSE(EtaC::IsNucleus());
-
-    REQUIRE(GetNucleusA(Code::Hydrogen) == 1);
-    REQUIRE(GetNucleusA(Code::Tritium) == 3);
-    REQUIRE(Hydrogen::GetNucleusZ() == 1);
-    REQUIRE(Tritium::GetNucleusA() == 3);
-
-    REQUIRE_THROWS(GetNucleusA(Code::Nucleus));
-    REQUIRE_THROWS(GetNucleusZ(Code::Nucleus));
+    CHECK_FALSE(IsNucleus(Code::Gamma));
+    CHECK(IsNucleus(Code::Argon));
+    CHECK_FALSE(IsNucleus(Code::Proton));
+    CHECK(IsNucleus(Code::Hydrogen));
+    CHECK(Argon::IsNucleus());
+    CHECK_FALSE(EtaC::IsNucleus());
+
+    CHECK(GetNucleusA(Code::Hydrogen) == 1);
+    CHECK(GetNucleusA(Code::Tritium) == 3);
+    CHECK(Hydrogen::GetNucleusZ() == 1);
+    CHECK(Tritium::GetNucleusA() == 3);
+
+    // Nucleus is a generic object, it has no specific properties
+    CHECK_THROWS(GetNucleusA(Code::Nucleus));
+    CHECK_THROWS(GetNucleusZ(Code::Nucleus));
+    CHECK_THROWS(GetMass(Code::Nucleus));
+    CHECK_THROWS(GetCharge(Code::Nucleus));
   }
 }

Processes/AnalyticProcessors/testExecTime.cc

@@ -33,86 +33,86 @@ TEST_CASE("Timing process", "[proccesses][analytic_processors ExecTime]") {
   SECTION("BoundaryCrossing") {
     ExecTime<DummyBoundaryCrossingProcess<10>> execTime;
     auto start = std::chrono::steady_clock::now();
-    REQUIRE(execTime.DoBoundaryCrossing(tmp, 0, 0) == EProcessReturn::eOk);
+    CHECK(execTime.DoBoundaryCrossing(tmp, 0, 0) == EProcessReturn::eOk);
     auto end = std::chrono::steady_clock::now();
-    REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
-            Approx(10).margin(5));
+    CHECK(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
+          Approx(10).margin(5));
 
     for (int i = 0; i < 100; i++) execTime.DoBoundaryCrossing(tmp, 0, 0);
 
-    REQUIRE(execTime.mean() == Approx(10 * 1000).margin(2 * 1000));
+    CHECK(execTime.mean() == Approx(10 * 1000).margin(2 * 1000));
 
-    REQUIRE(execTime.sumTime() == Approx(10 * 100 * 1000).margin((10 * 100) * 1000));
+    CHECK(execTime.sumTime() == Approx(10 * 100 * 1000).margin((10 * 100) * 1000));
 
-    REQUIRE(fabs(execTime.var()) < 20000);
+    CHECK(fabs(execTime.var()) < 100000);
   }
 
   SECTION("Continuous") {
     ExecTime<DummyContinuousProcess<50>> execTime;
     auto start = std::chrono::steady_clock::now();
-    REQUIRE(execTime.DoContinuous(tmp, tmp) == EProcessReturn::eOk);
+    CHECK(execTime.DoContinuous(tmp, tmp) == EProcessReturn::eOk);
     auto end = std::chrono::steady_clock::now();
-    REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
-            Approx(50).margin(5));
+    CHECK(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
+          Approx(50).margin(5));
 
     for (int i = 0; i < 100; i++) execTime.DoContinuous(tmp, tmp);
 
-    REQUIRE(execTime.mean() == Approx(50 * 1000).margin(2 * 1000));
+    CHECK(execTime.mean() == Approx(50 * 1000).margin(2 * 1000));
 
-    REQUIRE(execTime.sumTime() == Approx(50 * 100 * 1000).margin((10 * 100) * 1000));
+    CHECK(execTime.sumTime() == Approx(50 * 100 * 1000).margin((10 * 100) * 1000));
 
-    REQUIRE(fabs(execTime.var()) < 20000);
+    CHECK(fabs(execTime.var()) < 100000);
   }
 
   SECTION("Decay") {
     ExecTime<DummyDecayProcess<10>> execTime;
     auto start = std::chrono::steady_clock::now();
-    REQUIRE(execTime.DoDecay(tmp) == EProcessReturn::eOk);
+    CHECK(execTime.DoDecay(tmp) == EProcessReturn::eOk);
     auto end = std::chrono::steady_clock::now();
-    REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
-            Approx(10).margin(5));
+    CHECK(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
+          Approx(10).margin(5));
 
     for (int i = 0; i < 100; i++) execTime.DoDecay(tmp);
 
-    REQUIRE(execTime.mean() == Approx(10 * 1000).margin(2 * 1000));
+    CHECK(execTime.mean() == Approx(10 * 1000).margin(2 * 1000));
 
-    REQUIRE(execTime.sumTime() == Approx(10 * 100 * 100).margin((10 * 100) * 1000));
+    CHECK(execTime.sumTime() == Approx(10 * 100 * 100).margin((10 * 100) * 1000));
 
-    REQUIRE(fabs(execTime.var()) < 20000);
+    CHECK(fabs(execTime.var()) < 100000);
   }
 
   SECTION("Interaction") {
     ExecTime<DummyInteractionProcess<10>> execTime;
     auto start = std::chrono::steady_clock::now();
-    REQUIRE(execTime.DoInteraction(tmp) == EProcessReturn::eOk);
+    CHECK(execTime.DoInteraction(tmp) == EProcessReturn::eOk);
     auto end = std::chrono::steady_clock::now();
-    REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
-            Approx(10).margin(5));
+    CHECK(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
+          Approx(10).margin(5));
 
     for (int i = 0; i < 100; i++) execTime.DoInteraction(tmp);
 
-    REQUIRE(execTime.mean() == Approx(10 * 1000).margin(2 * 1000));
+    CHECK(execTime.mean() == Approx(10 * 1000).margin(2 * 1000));
 
-    REQUIRE(execTime.sumTime() == Approx(10 * 100 * 1000).margin((10 * 100) * 1000));
+    CHECK(execTime.sumTime() == Approx(10 * 100 * 1000).margin((10 * 100) * 1000));
 
-    REQUIRE(fabs(execTime.var()) < 20000);
+    CHECK(fabs(execTime.var()) < 100000);
   }
 
   SECTION("Secondaries") {
     ExecTime<DummySecondariesProcess<10>> execTime;
     auto start = std::chrono::steady_clock::now();
-    REQUIRE(execTime.DoSecondaries(tmp) == EProcessReturn::eOk);
+    CHECK(execTime.DoSecondaries(tmp) == EProcessReturn::eOk);
     auto end = std::chrono::steady_clock::now();
-    REQUIRE(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
-            Approx(10).margin(5));
+    CHECK(std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() ==
+          Approx(10).margin(5));
 
     for (int i = 0; i < 100; i++) execTime.DoSecondaries(tmp);
 
-    REQUIRE(execTime.mean() == Approx(10 * 1000).margin(2 * 1000));
+    CHECK(execTime.mean() == Approx(10 * 1000).margin(2 * 1000));
 
-    REQUIRE(execTime.sumTime() == Approx(10 * 100 * 1000).margin((10 * 100) * 1000));
+    CHECK(execTime.sumTime() == Approx(10 * 100 * 1000).margin((10 * 100) * 1000));
 
-    REQUIRE(fabs(execTime.var()) < 20000);
+    CHECK(fabs(execTime.var()) < 100000);
   }
 
   SECTION("TestMeanAlgo") {
@@ -128,7 +128,7 @@ TEST_CASE("Timing process", "[proccesses][analytic_processors ExecTime]") {
 
     std::vector<double> elems;
 
-    for (int i = 0; i < 1000000; i++) {
+    for (int i = 0; i < 1000; i++) {
       double timeDiv = distribution(generator);
 
       elems.push_back(timeDiv);
@@ -148,7 +148,7 @@ TEST_CASE("Timing process", "[proccesses][analytic_processors ExecTime]") {
       fMean2 += delta * delta2;
     }
 
-    REQUIRE(fN == 1000000);
+    CHECK(fN == 1000);
 
     double mean = 0;
     std::for_each(elems.begin(), elems.end(), [&](double i) { mean += i; });
@@ -159,10 +159,10 @@ TEST_CASE("Timing process", "[proccesses][analytic_processors ExecTime]") {
                   [&](double i) { var += (mean - i) * (mean - i); });
     var = var / fN;
 
-    REQUIRE(mean == Approx(10000.0).margin(10));
-    REQUIRE(var == Approx(200.0 * 200).margin(200));
+    CHECK(mean == Approx(10000.0).margin(10));
+    CHECK(var == Approx(200.0 * 200).margin(2000));
 
-    REQUIRE(fMean2 / fN == Approx(200 * 200).margin(200)); // Varianz
-    REQUIRE(fMean == Approx(10000).margin(10));
+    CHECK(fMean2 / fN == Approx(200 * 200).margin(2000)); // Varianz
+    CHECK(fMean == Approx(10000).margin(10));
   }
 }

Processes/EnergyLoss/EnergyLoss.cc

@@ -170,7 +170,7 @@ process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p, SetupTrack co
               dX / 1_g * square(1_cm));
   HEPEnergyType dE = TotalEnergyLoss(p, dX);
   auto E = p.GetEnergy();
-  const auto Ekin = E - p.GetMass();
+  [[maybe_unused]] const auto Ekin = E - p.GetMass();
   auto Enew = E + dE;
   C8LOG_DEBUG("EnergyLoss  dE={} MeV, E={} GeV, Ekin={} GeV, Enew={} GeV", dE / 1_MeV,
               E / 1_GeV, Ekin / 1_GeV, Enew / 1_GeV);
@@ -220,31 +220,33 @@ void EnergyLoss::FillProfile(SetupTrack const& vTrack, const HEPEnergyType dE) {
   GrammageType const grammageEnd = shower_axis_.projectedX(vTrack.GetPosition(1));
   const auto deltaX = grammageEnd - grammageStart;
 
-  const int binStart = grammageStart / dX_;
-  const int binEnd = grammageEnd / dX_;
+  int binStart = grammageStart / dX_;
+  if (binStart < 0) return;
+  int binEnd = grammageEnd / dX_;
+  if (binEnd > int(profile_.size() - 1)) return;
+  if (deltaX < dX_threshold_) return;
 
   C8LOG_DEBUG("energy deposit of -dE={} between {} and {}", -dE, grammageStart,
               grammageEnd);
 
   auto energyCount = HEPEnergyType::zero();
 
-  auto fill = [&](int bin, GrammageType weight) {
-    if (deltaX > dX_threshold_) {
-      auto const increment = -dE * weight / deltaX;
-      profile_[bin] += increment;
-      energyCount += increment;
+  auto fill = [&](const int bin, const double weight) {
+    auto const increment = -dE * weight;
+    profile_[bin] += increment;
+    energyCount += increment;
 
-      C8LOG_DEBUG("filling bin {} with weight {} : {} ", bin, weight, increment);
-    }
+    C8LOG_DEBUG("filling bin {} with weight {} : {} ", bin, weight, increment);
   };
 
   // fill longitudinal profile
-  fill(binStart, (1 + binStart) * dX_ - grammageStart);
-  fill(binEnd, grammageEnd - binEnd * dX_);
-
-  if (binStart == binEnd) { fill(binStart, -dX_); }
-
-  for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); }
+  if (binStart == binEnd) {
+    fill(binStart, 1);
+  } else {
+    fill(binStart, ((1 + binStart) * dX_ - grammageStart) / deltaX);
+    fill(binEnd, (grammageEnd - binEnd * dX_) / deltaX);
+    for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, 1); }
+  }
 
   C8LOG_DEBUG("total energy added to histogram: {} ", energyCount);
 }

Processes/Sibyll/Decay.cc

@@ -122,9 +122,9 @@ namespace corsika::process::sibyll {
     for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = abs(s_csydec_.idb[i]);
   }
 
-  void Decay::PrintDecayConfig(const particles::Code vCode) {
-    const int sibCode = process::sibyll::ConvertToSibyllRaw(vCode);
-    const int absSibCode = abs(sibCode);
+  void Decay::PrintDecayConfig([[maybe_unused]] const particles::Code vCode) {
+    [[maybe_unused]] const int sibCode = process::sibyll::ConvertToSibyllRaw(vCode);
+    [[maybe_unused]] const int absSibCode = abs(sibCode);
     C8LOG_DEBUG("Decay: Sibyll decay configuration: {} is {}", vCode,
                 (s_csydec_.idb[absSibCode - 1] <= 0) ? "stable" : "unstable");
   }
@@ -134,7 +134,7 @@ namespace corsika::process::sibyll {
     if (handleAllDecays_) {
       C8LOG_DEBUG("     all particles known to Sibyll are handled by Sibyll::Decay!");
     } else {
-      for (auto& pCode : handledDecays_) {
+      for ([[maybe_unused]] auto& pCode : handledDecays_) {
         C8LOG_DEBUG("      Decay of {}  is handled by Sibyll!", pCode);
       }
     }
@@ -154,7 +154,7 @@ namespace corsika::process::sibyll {
       const TimeType t0 = particles::GetLifetime(vP.GetPID());
       auto const lifetime = gamma * t0;
 
-      const auto mkin =
+      [[maybe_unused]] const auto mkin =
           (E * E - vP.GetMomentum().squaredNorm()); // delta_mass(vP.GetMomentum(), E, m);
       C8LOG_DEBUG("Sibyll::Decay: code: {} ", vP.GetPID());
       C8LOG_DEBUG("Sibyll::Decay: MinStep: t0: {} ", t0);
@@ -163,7 +163,7 @@ namespace corsika::process::sibyll {
                   vP.GetMomentum().GetComponents() / 1_GeV);
       C8LOG_DEBUG("Sibyll::Decay: momentum: shell mass-kin. inv. mass {} {}",
                   mkin / 1_GeV / 1_GeV, m / 1_GeV * m / 1_GeV);
-      auto sib_id = process::sibyll::ConvertToSibyllRaw(vP.GetPID());
+      [[maybe_unused]] auto sib_id = process::sibyll::ConvertToSibyllRaw(vP.GetPID());
       C8LOG_DEBUG("Sibyll::Decay: sib mass: {}", get_sibyll_mass2(sib_id));
       C8LOG_DEBUG("Sibyll::Decay: MinStep: gamma:  {}", gamma);
       C8LOG_DEBUG("Sibyll::Decay: MinStep: tau {} s: ", lifetime / 1_s);

Processes/Sibyll/Interaction.cc

@@ -218,10 +218,10 @@ namespace corsika::process::sibyll {
 
     // just for show:
     // boost projecticle
-    auto const PprojCoM = boost.toCoM(PprojLab);
+    [[maybe_unused]] auto const PprojCoM = boost.toCoM(PprojLab);
 
     // boost target
-    auto const PtargCoM = boost.toCoM(PtargLab);
+    [[maybe_unused]] auto const PtargCoM = boost.toCoM(PtargLab);
 
     C8LOG_DEBUG(
         "Interaction: ebeam CoM: {} GeV "

Processes/Sibyll/NuclearInteraction.cc

@@ -235,7 +235,7 @@ namespace corsika::process::sibyll {
     pTotLab += pTarget;
     auto const pTotLabNorm = pTotLab.norm();
     // calculate cm. energy
-    const HEPEnergyType ECoM = sqrt(
+    [[maybe_unused]] const HEPEnergyType ECoM = sqrt(
         (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
     auto const ECoMNN = sqrt(2. * ElabNuc * constants::nucleonMass);
     C8LOG_DEBUG(
@@ -406,10 +406,10 @@ namespace corsika::process::sibyll {
     // define boost to NUCLEON-NUCLEON frame
     COMBoost const boost(PprojNucLab, constants::nucleonMass);
     // boost projecticle
-    auto const PprojNucCoM = boost.toCoM(PprojNucLab);
+    [[maybe_unused]] auto const PprojNucCoM = boost.toCoM(PprojNucLab);
 
     // boost target
-    auto const PtargNucCoM = boost.toCoM(PtargNucLab);
+    [[maybe_unused]] auto const PtargNucCoM = boost.toCoM(PtargNucLab);
 
     C8LOG_DEBUG(
         fmt::format("Interaction: ebeam CoM: {} "

Stack/NuclearStackExtension/NuclearStackExtension.h

@@ -190,7 +190,9 @@ namespace corsika::stack {
         return InnerParticleInterface<StackIteratorInterface>::GetChargeNumber();
       }
 
-      int GetNucleusRef() const { return GetStackData().GetNucleusRef(GetIndex()); }
+      int GetNucleusRef() const {
+        return GetStackData().GetNucleusRef(GetIndex());
+      } // LCOV_EXCL_LINE
 
     protected:
       void SetNucleusRef(const int vR) { GetStackData().SetNucleusRef(GetIndex(), vR); }

Stack/NuclearStackExtension/testNuclearStackExtension.cc

@@ -119,6 +119,7 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
 
     ParticleDataStack s;
     // add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
+    // i=9, 19, 29, etc. are nuclei
     for (int i = 0; i < 99; ++i) {
       if ((i + 1) % 10 == 0) {
         s.AddParticle(std::make_tuple(
@@ -167,6 +168,21 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
       CHECK(p93.GetTime() == 100_s);
     }
 
+    // copy
+    {
+      s.Copy(s.begin() + 89, s.begin() + 79); // nuclei to nuclei
+      const auto& p89 = s.cbegin() + 89;
+      const auto& p79 = s.cbegin() + 79;
+
+      CHECK(p89.GetPID() == particles::Code::Nucleus);
+      CHECK(p89.GetEnergy() == 89 * 15_GeV);
+      CHECK(p89.GetTime() == 100_s);
+
+      CHECK(p79.GetPID() == particles::Code::Nucleus);
+      CHECK(p79.GetEnergy() == 89 * 15_GeV);
+      CHECK(p79.GetTime() == 100_s);
+    }
+
     // swap
     {
       s.Swap(s.begin() + 11, s.begin() + 10);
@@ -206,4 +222,37 @@ TEST_CASE("NuclearStackExtension", "[stack]") {
     for (int i = 0; i < 99; ++i) s.last().Delete();
     CHECK(s.getEntries() == 0);
   }
+
+  SECTION("not allowed") {
+    NuclearStackExtension<corsika::stack::super_stupid::SuperStupidStack,
+                          ExtendedParticleInterfaceType>
+        s;
+
+    // not valid:
+    CHECK_THROWS(s.AddParticle(std::make_tuple(
+        particles::Code::Oxygen, 1.5_GeV,
+        corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+        Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 16, 8)));
+
+    // valid
+    auto particle = s.AddParticle(
+        std::make_tuple(particles::Code::Nucleus, 1.5_GeV,
+                        corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+                        Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 9));
+
+    // not valid
+    CHECK_THROWS(particle.AddSecondary(std::make_tuple(
+        particles::Code::Oxygen, 1.5_GeV,
+        corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+        Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 16, 8)));
+
+    // add a another nucleus, so there are two now
+    s.AddParticle(
+        std::make_tuple(particles::Code::Nucleus, 1.5_GeV,
+                        corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+                        Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 9));
+
+    // not valid, since end() is not a valid entry
+    CHECK_THROWS(s.Swap(s.begin(), s.end()));
+  }
 }