diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 1e240145b65201a9b1b6448e75d7ea1cdc1b4cf9..64ed627e03c8c7687a2b75330dd50d46e1168627 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/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;
diff --git a/Framework/Analytics/testClassTimer.cc b/Framework/Analytics/testClassTimer.cc
index 290eabc0e3e8362c126ca2b6a13ee3b47e1b19b8..95f82b4b2caa93394634e5c39ee542dacda94f2f 100644
--- a/Framework/Analytics/testClassTimer.cc
+++ b/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);
}
}
diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc
index e01f0b1fb8f93ab0db33c4cf4b270eab19720fbb..8b3f7184c6b31d5b447dbdb173847cdfc0d080ae 100644
--- a/Framework/Particles/testParticles.cc
+++ b/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));
}
}
diff --git a/Processes/AnalyticProcessors/testExecTime.cc b/Processes/AnalyticProcessors/testExecTime.cc
index dbda248c889279e58a0feb4be6650e02eb02bb7a..e2fc710e7a201096184edac399befe01a8d04356 100644
--- a/Processes/AnalyticProcessors/testExecTime.cc
+++ b/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));
}
}
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index 42efc08a51303fb3bc347d527053cd9a74e0247b..336294d8e97ee9be5fdc3c82b46b30365849e826 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/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);
}
diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
index a85541adca9b3ba6a21382ec552129907cd4f9ad..8297739b793e33d7b4fe9b65bfd66def992f372e 100644
--- a/Processes/Sibyll/Decay.cc
+++ b/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);
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index f020922f077f2e0b70093746b66099f7dddc0a0e..aa222b7f65ced91a80b5e6402a0ae5a45926c6dc 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/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 "
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index c53beeae7a6934dd090af9ea5302deacf342b681..de77d5d33a32b1cb4461abb88552ee14bea415a0 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/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: {} "
diff --git a/Stack/NuclearStackExtension/NuclearStackExtension.h b/Stack/NuclearStackExtension/NuclearStackExtension.h
index 4409a215eed0c331e922b311235f37e10ded24d8..c2ac33e2e0f25079dbd535ed52360bbb0c1e6244 100644
--- a/Stack/NuclearStackExtension/NuclearStackExtension.h
+++ b/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); }
diff --git a/Stack/NuclearStackExtension/testNuclearStackExtension.cc b/Stack/NuclearStackExtension/testNuclearStackExtension.cc
index 6896953a6079ba0c66b4ba630cf96a2c78a3bd48..72017c91a05c70d9c062fd72c20f0ff31cc4e246 100644
--- a/Stack/NuclearStackExtension/testNuclearStackExtension.cc
+++ b/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()));
+ }
}