From c1dd792630d5de1b69d15f3e7d8bc1905fa0cb93 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Wed, 7 Oct 2020 22:10:02 +0200
Subject: [PATCH] better interface for CONEX, better output
---
.../Examples/cascade_proton_example.cc | 7 ++--
Documentation/Examples/hybrid_MC.cc | 23 +++++------
Documentation/Examples/vertical_EAS.cc | 4 +-
.../LayeredSphericalAtmosphereBuilder.h | 4 +-
Environment/MediumTypes.h | 6 +--
Environment/testEnvironment.cc | 4 +-
Processes/CONEXSourceCut/CONEXSourceCut.cc | 41 ++++++++++---------
Processes/CONEXSourceCut/CONEXSourceCut.h | 5 +--
Processes/EnergyLoss/EnergyLoss.cc | 12 ++++++
Processes/EnergyLoss/EnergyLoss.h | 5 +++
.../testInteractionCounter.cc | 7 ++--
Processes/Proposal/ContinuousProcess.cc | 4 +-
Processes/Proposal/ContinuousProcess.h | 6 +--
Processes/Pythia/testPythia8.cc | 15 ++++---
Processes/QGSJetII/testQGSJetII.cc | 8 ++--
Processes/Sibyll/NuclearInteraction.cc | 9 ++--
Processes/Sibyll/testSibyll.cc | 22 ++++++----
Processes/UrQMD/testUrQMD.cc | 26 +++++++-----
Setup/SetupStack.h | 24 +++++------
19 files changed, 129 insertions(+), 103 deletions(-)
diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc
index 51433db11..44e2bfa01 100644
--- a/Documentation/Examples/cascade_proton_example.cc
+++ b/Documentation/Examples/cascade_proton_example.cc
@@ -74,14 +74,13 @@ int main() {
auto& universe = *(env.GetUniverse());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
- auto theMedium =
- EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
+ auto theMedium = EnvType::CreateNode<Sphere>(
+ Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel =
environment::UniformMediumType<environment::UniformMagneticField<
environment::HomogeneousMedium<setup::EnvironmentInterface>>>;
-
+
theMedium->SetModelProperties<MyHomogeneousModel>(
environment::EMediumType::eAir, geometry::Vector(rootCS, 0_T, 0_T, 1_T),
1_kg / (1_m * 1_m * 1_m),
diff --git a/Documentation/Examples/hybrid_MC.cc b/Documentation/Examples/hybrid_MC.cc
index 407107f63..bad8cee91 100644
--- a/Documentation/Examples/hybrid_MC.cc
+++ b/Documentation/Examples/hybrid_MC.cc
@@ -23,18 +23,17 @@
#include <corsika/logging/Logging.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/StackProcess.h>
+#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
#include <corsika/process/observation_plane/ObservationPlane.h>
#include <corsika/process/on_shell_check/OnShellCheck.h>
-#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/switch_process/SwitchProcess.h>
-#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/process/urqmd/UrQMD.h>
#include <corsika/random/RNGManager.h>
@@ -76,7 +75,6 @@ void registerRandomStreams(const int seed) {
template <typename T>
using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
-
int main(int argc, char** argv) {
logging::SetLevel(logging::level::info);
@@ -100,7 +98,8 @@ int main(int argc, char** argv) {
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
- environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface> builder{center};
+ environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface> builder{
+ center};
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
@@ -216,8 +215,7 @@ int main(int argc, char** argv) {
process::energy_loss::EnergyLoss eLoss{showerAxis, cut.GetECut()};
corsika::process::conex_source_cut::CONEXSourceCut conex(
- center, showerAxis, t, injectionHeight, E0,
- particles::GetPDG(particles::Code::Proton));
+ center, showerAxis, t, injectionHeight, E0, particles::Code::Proton);
process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
@@ -237,8 +235,8 @@ int main(int argc, char** argv) {
55_GeV);
auto decaySequence = decayPythia << decaySibyll;
- auto sequence = switchProcess << reset_particle_mass << decaySequence
- << eLoss << cut << conex << longprof << observationLevel;
+ auto sequence = switchProcess << reset_particle_mass << decaySequence << eLoss << cut
+ << conex << longprof << observationLevel;
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
@@ -251,22 +249,21 @@ int main(int argc, char** argv) {
EAS.Run();
cut.ShowResults();
- //eLoss.ShowResults();
+ eLoss.showResults();
observationLevel.ShowResults();
const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
- cut.GetEmEnergy() + //eLoss.GetEnergyLost() +
+ cut.GetEmEnergy() + eLoss.energyLost() +
observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.Reset();
cut.Reset();
- //eLoss.Reset();
+ eLoss.reset();
- //conex.addParticle(0, E0, 0_eV, emPosition, momentum.normalized(), 0_s);
conex.SolveCE();
auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
- urqmdCounted.GetHistogram();
+ urqmdCounted.GetHistogram();
hists.saveLab("inthist_lab.txt");
hists.saveCMS("inthist_cms.txt");
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 6c602705b..2d2a1b7cb 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -78,7 +78,6 @@ void registerRandomStreams(const int seed) {
template <typename T>
using MEnv = environment::UniformMediumType<environment::UniformMagneticField<T>>;
-
int main(int argc, char** argv) {
logging::SetLevel(logging::level::info);
@@ -102,7 +101,8 @@ int main(int argc, char** argv) {
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
- environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface> builder{center};
+ environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface> builder{
+ center};
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
diff --git a/Environment/LayeredSphericalAtmosphereBuilder.h b/Environment/LayeredSphericalAtmosphereBuilder.h
index cbb705bbf..01f4855af 100644
--- a/Environment/LayeredSphericalAtmosphereBuilder.h
+++ b/Environment/LayeredSphericalAtmosphereBuilder.h
@@ -7,12 +7,12 @@
*/
#include <corsika/environment/Environment.h>
+#include <corsika/environment/FlatExponential.h>
+#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/SlidingPlanarExponential.h>
#include <corsika/environment/VolumeTreeNode.h>
-#include <corsika/environment/FlatExponential.h>
-#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalConstants.h>
diff --git a/Environment/MediumTypes.h b/Environment/MediumTypes.h
index b2f7a0aaa..bc44560fd 100644
--- a/Environment/MediumTypes.h
+++ b/Environment/MediumTypes.h
@@ -14,9 +14,9 @@ namespace corsika::environment {
* Medium types are useful most importantly for effective models
* like energy losses. a particular medium (mixture of components)
* may have specif properties not reflected by its mixture of
- * components.
+ * components.
*/
-
+
enum class EMediumType { eUnknown, eAir, eWater, eIce, eRock };
-}
+} // namespace corsika::environment
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index 99aec33f0..5480aec94 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -327,12 +327,12 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
auto const R = builder.getEarthRadius();
// check magnetic field at several locations
- const Point pTest(gCS, -10_m, 4_m, R+35_m);
+ const Point pTest(gCS, -10_m, 4_m, R + 35_m);
CHECK(B0.GetComponents(gCS) == univ->GetContainingNode(pTest)
->GetModelProperties()
.GetMagneticField(pTest)
.GetComponents(gCS));
- const Point pTest2(gCS, 10_m, -4_m, R+15_km);
+ const Point pTest2(gCS, 10_m, -4_m, R + 15_km);
CHECK(B1.GetComponents(gCS) == univ->GetContainingNode(pTest2)
->GetModelProperties()
.GetMagneticField(pTest2)
diff --git a/Processes/CONEXSourceCut/CONEXSourceCut.cc b/Processes/CONEXSourceCut/CONEXSourceCut.cc
index 94059e886..f8d7b0911 100644
--- a/Processes/CONEXSourceCut/CONEXSourceCut.cc
+++ b/Processes/CONEXSourceCut/CONEXSourceCut.cc
@@ -6,10 +6,12 @@
* the license.
*/
+#include <corsika/logging/Logging.h>
#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
#include <corsika/process/conex_source_cut/CONEX_f.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalConstants.h>
+
#include <algorithm>
#include <fstream>
#include <iomanip>
@@ -23,32 +25,28 @@ using namespace corsika::setup;
corsika::process::EProcessReturn CONEXSourceCut::DoSecondaries(
corsika::setup::StackView& vS) {
auto p = vS.begin();
-
while (p != vS.end()) {
Code const pid = p.GetPID();
-
- auto const it = std::find_if(egs_em_codes_.cbegin(), egs_em_codes_.cend(),
- [=](auto const& p) { return pid == p.first; });
- if (it != egs_em_codes_.cend()) {
- // EM particle
-
- auto const egs_pid = it->second;
-
- addParticle(egs_pid, p.GetEnergy(), p.GetMass(), p.GetPosition(),
- p.GetMomentum().normalized(), p.GetTime());
-
+ if (addParticle(pid, p.GetEnergy(), p.GetMass(), p.GetPosition(),
+ p.GetMomentum().normalized(), p.GetTime())) {
p.Delete();
}
++p;
}
-
return corsika::process::EProcessReturn::eOk;
}
-void CONEXSourceCut::addParticle(int egs_pid, HEPEnergyType energy, HEPEnergyType mass,
- geometry::Point const& position,
+bool CONEXSourceCut::addParticle(particles::Code pid, HEPEnergyType energy,
+ HEPEnergyType mass, geometry::Point const& position,
geometry::Vector<dimensionless_d> const& direction,
TimeType t) {
+
+ auto const it = std::find_if(egs_em_codes_.cbegin(), egs_em_codes_.cend(),
+ [=](auto const& p) { return pid == p.first; });
+ if (it == egs_em_codes_.cend()) { return false; }
+
+ // EM particle
+ auto const egs_pid = it->second;
std::cout << "position conexObs: " << position.GetCoordinates(conexObservationCS_)
<< std::endl;
@@ -122,6 +120,8 @@ void CONEXSourceCut::addParticle(int egs_pid, HEPEnergyType energy, HEPEnergyTyp
int n = 1, i = 1;
conex::cegs4_(n, i);
+
+ return true;
}
void CONEXSourceCut::SolveCE() {
@@ -160,9 +160,12 @@ void CONEXSourceCut::SolveCE() {
conex::get_shower_hadron_(icuth, nX, Hadrons[0]);
std::ofstream file{"conex_output.txt"};
+ file << fmt::format("#{:>8} {:>13} {:>13} {:>13} {:>13} {:>13} {:>13} {:>13}\n", "X",
+ "N", "dEdX", "Mu", "dMu", "Gamma", "El", "Had");
for (int i = 0; i < nX; ++i) {
- file << X[i] << " " << N[i] << " " << dEdX[i] << " " << Mu[i] << " " << dMu[i] << " "
- << Gamma[i] << " " << Electrons[i] << " " << Hadrons[i] << std::endl;
+ file << fmt::format(
+ " {:>8.2f} {:>13.3} {:>13.3} {:>13.3} {:>13.3} {:>13.3} {:>13.3} {:>13.3}\n",
+ X[i], N[i], dEdX[i], Mu[i], dMu[i], Gamma[i], Electrons[i], Hadrons[i]);
}
std::ofstream fitout{"conex_fit.txt"};
@@ -186,7 +189,7 @@ CONEXSourceCut::CONEXSourceCut(geometry::Point center,
units::si::LengthType groundDist,
units::si::LengthType injectionHeight,
units::si::HEPEnergyType primaryEnergy,
- particles::PDGCode primaryID)
+ particles::Code primaryPID)
: center_{center}
, showerAxis_{showerAxis}
, groundDist_{groundDist}
@@ -278,7 +281,7 @@ CONEXSourceCut::CONEXSourceCut(geometry::Point center,
std::cout << "theta (deg) = " << theta << "; phi (deg) = " << phi << std::endl;
- int ipart = static_cast<int>(primaryID);
+ int ipart = static_cast<int>(particles::GetPDG(primaryPID));
auto rng = corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
double dimpact = 0.; // valid only if shower core is fixed on the observation plane; for
diff --git a/Processes/CONEXSourceCut/CONEXSourceCut.h b/Processes/CONEXSourceCut/CONEXSourceCut.h
index 121d8c163..5923cdbf0 100644
--- a/Processes/CONEXSourceCut/CONEXSourceCut.h
+++ b/Processes/CONEXSourceCut/CONEXSourceCut.h
@@ -31,13 +31,12 @@ namespace corsika::process {
CONEXSourceCut(geometry::Point center, environment::ShowerAxis const& showerAxis,
units::si::LengthType groundDist,
units::si::LengthType injectionHeight,
- units::si::HEPEnergyType primaryEnergy,
- particles::PDGCode primaryID);
+ units::si::HEPEnergyType primaryEnergy, particles::Code pid);
corsika::process::EProcessReturn DoSecondaries(corsika::setup::StackView&);
void SolveCE();
- void addParticle(int egs_pid, units::si::HEPEnergyType energy,
+ bool addParticle(particles::Code pid, units::si::HEPEnergyType energy,
units::si::HEPEnergyType mass, geometry::Point const& position,
geometry::Vector<units::si::dimensionless_d> const& direction,
units::si::TimeType t);
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index 336294d8e..66c8a76e2 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -265,3 +265,15 @@ void EnergyLoss::PrintProfile() const {
HEPEnergyType EnergyLoss::GetTotal() const {
return std::accumulate(profile_.cbegin(), profile_.cend(), HEPEnergyType::zero());
}
+
+void EnergyLoss::showResults() const {
+ using namespace corsika::units::si; // required for operator::_MeV
+ std::cout << " ******************************" << std::endl
+ << " PROCESS::ContinuousProcess: " << std::endl;
+ std::cout << " energy lost dE (GeV) : " << energy_lost_ / 1_GeV << std::endl;
+}
+
+void EnergyLoss::reset() {
+ using namespace corsika::units::si; // required for operator::_MeV
+ energy_lost_ = 0_GeV;
+}
diff --git a/Processes/EnergyLoss/EnergyLoss.h b/Processes/EnergyLoss/EnergyLoss.h
index 89d20bf05..f60e16b35 100644
--- a/Processes/EnergyLoss/EnergyLoss.h
+++ b/Processes/EnergyLoss/EnergyLoss.h
@@ -45,6 +45,10 @@ namespace corsika::process::energy_loss {
static units::si::HEPEnergyType TotalEnergyLoss(setup::Stack::ParticleType const&,
const units::si::GrammageType);
+ void showResults() const;
+ void reset();
+ corsika::units::si::HEPEnergyType energyLost() const { return energy_lost_; }
+
private:
void FillProfile(setup::Trajectory const&, units::si::HEPEnergyType);
@@ -57,6 +61,7 @@ namespace corsika::process::energy_loss {
environment::ShowerAxis const& shower_axis_;
corsika::units::si::HEPEnergyType emCut_;
std::vector<units::si::HEPEnergyType> profile_; // longitudinal profile
+ units::si::HEPEnergyType energy_lost_ = 0 * units::si::electronvolt;
};
units::si::GrammageType const dX_threshold_ = std::invoke([]() {
diff --git a/Processes/InteractionCounter/testInteractionCounter.cc b/Processes/InteractionCounter/testInteractionCounter.cc
index f8ded748c..762c615cd 100644
--- a/Processes/InteractionCounter/testInteractionCounter.cc
+++ b/Processes/InteractionCounter/testInteractionCounter.cc
@@ -27,7 +27,6 @@ using namespace corsika::process::interaction_counter;
using namespace corsika::units;
using namespace corsika::units::si;
-
struct DummyProcess {
template <typename TParticle>
GrammageType GetInteractionLength([[maybe_unused]] TParticle const& particle) {
@@ -56,7 +55,8 @@ TEST_CASE("InteractionCounter") {
SECTION("DoInteraction nucleus") {
unsigned short constexpr A = 14, Z = 7;
- auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::Nucleus, A, Z, 105_TeV, nodePtr, *csPtr);
+ auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::Nucleus, A,
+ Z, 105_TeV, nodePtr, *csPtr);
REQUIRE(stackPtr->getEntries() == 1);
REQUIRE(secViewPtr->getEntries() == 0);
@@ -76,7 +76,8 @@ TEST_CASE("InteractionCounter") {
SECTION("DoInteraction Lambda") {
auto constexpr code = particles::Code::Lambda0;
auto constexpr codeInt = static_cast<particles::CodeIntType>(code);
- auto [stackPtr, secViewPtr] = setup::testing::setupStack(code, 0,0, 105_TeV, nodePtr, *csPtr);
+ auto [stackPtr, secViewPtr] =
+ setup::testing::setupStack(code, 0, 0, 105_TeV, nodePtr, *csPtr);
REQUIRE(stackPtr->getEntries() == 1);
REQUIRE(secViewPtr->getEntries() == 0);
diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
index 5d758be3b..fb4d528fa 100644
--- a/Processes/Proposal/ContinuousProcess.cc
+++ b/Processes/Proposal/ContinuousProcess.cc
@@ -146,14 +146,14 @@ namespace corsika::process::proposal {
return dist;
}
- void ContinuousProcess::ShowResults() const {
+ void ContinuousProcess::showResults() const {
using namespace corsika::units::si; // required for operator::_MeV
std::cout << " ******************************" << std::endl
<< " PROCESS::ContinuousProcess: " << std::endl;
std::cout << " energy lost dE (GeV) : " << energy_lost_ / 1_GeV << std::endl;
}
- void ContinuousProcess::Reset() {
+ void ContinuousProcess::reset() {
using namespace corsika::units::si; // required for operator::_MeV
energy_lost_ = 0_GeV;
}
diff --git a/Processes/Proposal/ContinuousProcess.h b/Processes/Proposal/ContinuousProcess.h
index 7d9a88d28..6d6641e52 100644
--- a/Processes/Proposal/ContinuousProcess.h
+++ b/Processes/Proposal/ContinuousProcess.h
@@ -72,8 +72,8 @@ namespace corsika::process::proposal {
template <typename Particle, typename Track>
corsika::units::si::LengthType MaxStepLength(Particle const&, Track const&);
- void ShowResults() const;
- void Reset();
- corsika::units::si::HEPEnergyType GetEnergyLost() const { return energy_lost_; }
+ void showResults() const;
+ void reset();
+ corsika::units::si::HEPEnergyType energyLost() const { return energy_lost_; }
};
} // namespace corsika::process::proposal
diff --git a/Processes/Pythia/testPythia8.cc b/Processes/Pythia/testPythia8.cc
index f64de2b0c..063d318bb 100644
--- a/Processes/Pythia/testPythia8.cc
+++ b/Processes/Pythia/testPythia8.cc
@@ -77,15 +77,15 @@ TEST_CASE("Pythia", "[processes]") {
#include <corsika/units/PhysicalUnits.h>
#include <corsika/particles/ParticleProperties.h>
+#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
-#include <corsika/setup/SetupEnvironment.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
-#include <corsika/environment/UniformMediumType.h>
#include <corsika/environment/UniformMagneticField.h>
+#include <corsika/environment/UniformMediumType.h>
using namespace corsika;
using namespace corsika::units::si;
@@ -105,12 +105,14 @@ TEST_CASE("pythia process") {
auto const& cs = *csPtr;
[[maybe_unused]] auto const& env_dummy = env;
[[maybe_unused]] auto const& node_dummy = nodePtr;
-
+
SECTION("pythia decay") {
feenableexcept(FE_INVALID);
const HEPEnergyType P0 = 10_GeV;
- auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::PiPlus, 0, 0, P0, nodePtr, *csPtr);
- const auto plab = corsika::stack::MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
+ auto [stackPtr, secViewPtr] =
+ setup::testing::setupStack(particles::Code::PiPlus, 0, 0, P0, nodePtr, *csPtr);
+ const auto plab = corsika::stack::MomentumVector(
+ cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
auto& stack = *stackPtr;
auto& view = *secViewPtr;
auto particle = stackPtr->first();
@@ -154,7 +156,8 @@ TEST_CASE("pythia process") {
SECTION("pythia interaction") {
feenableexcept(FE_INVALID);
- auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::PiPlus, 0, 0, 100_GeV, nodePtr, *csPtr);
+ auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::PiPlus, 0,
+ 0, 100_GeV, nodePtr, *csPtr);
auto& view = *secViewPtr;
auto particle = stackPtr->first();
diff --git a/Processes/QGSJetII/testQGSJetII.cc b/Processes/QGSJetII/testQGSJetII.cc
index c281d9908..0dd569b50 100644
--- a/Processes/QGSJetII/testQGSJetII.cc
+++ b/Processes/QGSJetII/testQGSJetII.cc
@@ -110,15 +110,15 @@ TEST_CASE("QgsjetII", "[processes]") {
#include <corsika/units/PhysicalUnits.h>
#include <corsika/particles/ParticleProperties.h>
-#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupEnvironment.h>
+#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
-#include <corsika/environment/UniformMediumType.h>
#include <corsika/environment/UniformMagneticField.h>
+#include <corsika/environment/UniformMediumType.h>
using namespace corsika::units::si;
using namespace corsika::units;
@@ -134,8 +134,8 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
SECTION("InteractionInterface") {
- auto [stackPtr, secViewPtr] =
- setup::testing::setupStack(particles::Code::Proton, 0,0, 110_GeV, nodePtr, *csPtr);
+ auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::Proton, 0,
+ 0, 110_GeV, nodePtr, *csPtr);
setup::StackView& view = *(secViewPtr.get());
auto particle = stackPtr->first();
auto projectile = secViewPtr->GetProjectile();
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index 51d33a614..4c22491a4 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -133,7 +133,8 @@ namespace corsika::process::sibyll {
}
template <>
- units::si::CrossSectionType NuclearInteraction<setup::Environment>::ReadCrossSectionTable(
+ units::si::CrossSectionType
+ NuclearInteraction<setup::Environment>::ReadCrossSectionTable(
const int ia, particles::Code pTarget, units::si::HEPEnergyType elabnuc) {
using namespace corsika::particles;
using namespace units::si;
@@ -153,8 +154,8 @@ namespace corsika::process::sibyll {
template <>
template <>
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
- NuclearInteraction<setup::Environment>::GetCrossSection(Particle const& vP,
- const particles::Code TargetId) {
+ NuclearInteraction<setup::Environment>::GetCrossSection(
+ Particle const& vP, const particles::Code TargetId) {
using namespace units::si;
if (vP.GetPID() != particles::Code::Nucleus)
throw std::runtime_error(
@@ -615,7 +616,7 @@ namespace corsika::process::sibyll {
template <>
NuclearInteraction<setup::Environment>::NuclearInteraction(
- process::sibyll::Interaction& hadint, setup::Environment const& env)
+ process::sibyll::Interaction& hadint, setup::Environment const& env)
: environment_(env)
, hadronicInteraction_(hadint) {
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 69cde670f..950977e70 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -106,12 +106,14 @@ TEST_CASE("SibyllInterface", "[processes]") {
SECTION("InteractionInterface - low energy") {
const HEPEnergyType P0 = 60_GeV;
- auto [stack, viewPtr] = setup::testing::setupStack(particles::Code::Proton, 0,0, P0, nodePtr, cs);
- const auto plab = corsika::stack::MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
+ auto [stack, viewPtr] =
+ setup::testing::setupStack(particles::Code::Proton, 0, 0, P0, nodePtr, cs);
+ const auto plab = corsika::stack::MomentumVector(
+ cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
setup::StackView& view = *viewPtr;
-
+
auto particle = stack->first();
-
+
Interaction model;
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view);
@@ -180,13 +182,14 @@ TEST_CASE("SibyllInterface", "[processes]") {
CHECK((pSum - plab).norm() / 1_GeV == Approx(0).margin(plab.norm() * 0.05 / 1_GeV));
CHECK(pSum.norm() / P0 == Approx(1).margin(0.05));
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
- CHECK(length/1_g*1_cm*1_cm == Approx(88.7).margin(0.1));
+ CHECK(length / 1_g * 1_cm * 1_cm == Approx(88.7).margin(0.1));
CHECK(view.getSize() == 20);
}
SECTION("NuclearInteractionInterface") {
- auto [stack, viewPtr] = setup::testing::setupStack(particles::Code::Nucleus, 4, 2, 500_GeV, nodePtr, cs);
+ auto [stack, viewPtr] =
+ setup::testing::setupStack(particles::Code::Nucleus, 4, 2, 500_GeV, nodePtr, cs);
setup::StackView& view = *viewPtr;
auto particle = stack->first();
@@ -195,15 +198,16 @@ TEST_CASE("SibyllInterface", "[processes]") {
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view);
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
- CHECK(length/1_g*1_cm*1_cm == Approx(44.2).margin(.1));
+ CHECK(length / 1_g * 1_cm * 1_cm == Approx(44.2).margin(.1));
CHECK(view.getSize() == 11);
}
SECTION("DecayInterface") {
- auto [stackPtr, viewPtr] = setup::testing::setupStack(particles::Code::Lambda0, 0,0, 10_GeV, nodePtr, cs);
+ auto [stackPtr, viewPtr] =
+ setup::testing::setupStack(particles::Code::Lambda0, 0, 0, 10_GeV, nodePtr, cs);
setup::StackView& view = *viewPtr;
- auto& stack = *stackPtr;
+ auto& stack = *stackPtr;
auto particle = stack.first();
Decay model;
diff --git a/Processes/UrQMD/testUrQMD.cc b/Processes/UrQMD/testUrQMD.cc
index 85616887e..17a010c2f 100644
--- a/Processes/UrQMD/testUrQMD.cc
+++ b/Processes/UrQMD/testUrQMD.cc
@@ -53,7 +53,6 @@ auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS)
return sum;
}
-
TEST_CASE("UrQMD") {
SECTION("conversion") {
REQUIRE_THROWS(process::UrQMD::ConvertFromUrQMD(106, 0));
@@ -67,7 +66,8 @@ TEST_CASE("UrQMD") {
UrQMD urqmd;
SECTION("interaction length") {
- auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Nitrogen);
+ auto [env, csPtr, nodePtr] =
+ setup::testing::setupEnvironment(particles::Code::Nitrogen);
auto const& cs = *csPtr;
[[maybe_unused]] auto const& env_dummy = env;
[[maybe_unused]] auto const& node_dummy = nodePtr;
@@ -78,7 +78,7 @@ TEST_CASE("UrQMD") {
particles::Code::K0, particles::Code::K0Bar, particles::Code::K0Long};
for (auto code : validProjectileCodes) {
- auto [stack, view] = setup::testing::setupStack(code, 0,0, 100_GeV, nodePtr, cs);
+ auto [stack, view] = setup::testing::setupStack(code, 0, 0, 100_GeV, nodePtr, cs);
REQUIRE(stack->getEntries() == 1);
REQUIRE(view->getEntries() == 0);
@@ -89,12 +89,14 @@ TEST_CASE("UrQMD") {
}
SECTION("nucleus projectile") {
- auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Oxygen);
+ auto [env, csPtr, nodePtr] =
+ setup::testing::setupEnvironment(particles::Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env; // against warnings
[[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
unsigned short constexpr A = 14, Z = 7;
- auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::Nucleus, A, Z, 400_GeV, nodePtr, *csPtr);
+ auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::Nucleus, A,
+ Z, 400_GeV, nodePtr, *csPtr);
REQUIRE(stackPtr->getEntries() == 1);
REQUIRE(secViewPtr->getEntries() == 0);
@@ -113,12 +115,13 @@ TEST_CASE("UrQMD") {
}
SECTION("\"special\" projectile") {
- auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Oxygen);
+ auto [env, csPtr, nodePtr] =
+ setup::testing::setupEnvironment(particles::Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env; // against warnings
[[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
- auto [stackPtr, secViewPtr] =
- setup::testing::setupStack(particles::Code::PiPlus, 0,0, 400_GeV, nodePtr, *csPtr);
+ auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::PiPlus, 0,
+ 0, 400_GeV, nodePtr, *csPtr);
REQUIRE(stackPtr->getEntries() == 1);
REQUIRE(secViewPtr->getEntries() == 0);
@@ -139,12 +142,13 @@ TEST_CASE("UrQMD") {
}
SECTION("K0Long projectile") {
- auto [env, csPtr, nodePtr] = setup::testing::setupEnvironment(particles::Code::Oxygen);
+ auto [env, csPtr, nodePtr] =
+ setup::testing::setupEnvironment(particles::Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env; // against warnings
[[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
- auto [stackPtr, secViewPtr] =
- setup::testing::setupStack(particles::Code::K0Long,0,0, 400_GeV, nodePtr, *csPtr);
+ auto [stackPtr, secViewPtr] = setup::testing::setupStack(particles::Code::K0Long, 0,
+ 0, 400_GeV, nodePtr, *csPtr);
REQUIRE(stackPtr->getEntries() == 1);
REQUIRE(secViewPtr->getEntries() == 0);
diff --git a/Setup/SetupStack.h b/Setup/SetupStack.h
index 9b4cb8712..f06aaf760 100644
--- a/Setup/SetupStack.h
+++ b/Setup/SetupStack.h
@@ -131,17 +131,17 @@ namespace corsika::setup {
} // namespace corsika::setup
- /**
- * standard stack setup for unit tests. This can be moved to "test"
- * directory, when available.
- */
+/**
+ * standard stack setup for unit tests. This can be moved to "test"
+ * directory, when available.
+ */
namespace corsika::setup::testing {
inline auto setupStack(particles::Code vProjectileType, int vA, int vZ,
- units::si::HEPEnergyType vMomentum,
- const setup::Environment::BaseNodeType* vNodePtr,
- geometry::CoordinateSystem const& cs) {
+ units::si::HEPEnergyType vMomentum,
+ const setup::Environment::BaseNodeType* vNodePtr,
+ geometry::CoordinateSystem const& cs) {
using namespace corsika;
using namespace corsika::units::si;
@@ -157,18 +157,16 @@ namespace corsika::setup::testing {
auto particle = stack->AddParticle(
std::make_tuple(particles::Code::Nucleus, E0, pLab, origin, 0_ns, vA, vZ));
particle.SetNode(vNodePtr);
- return std::make_tuple(
- std::move(stack),
- std::make_unique<setup::StackView>(particle));
+ return std::make_tuple(std::move(stack),
+ std::make_unique<setup::StackView>(particle));
} else { // not a nucleus
HEPEnergyType const E0 = sqrt(
units::static_pow<2>(particles::GetMass(vProjectileType)) + pLab.squaredNorm());
auto particle =
stack->AddParticle(std::make_tuple(vProjectileType, E0, pLab, origin, 0_ns));
particle.SetNode(vNodePtr);
- return std::make_tuple(
- std::move(stack),
- std::make_unique<setup::StackView>(particle));
+ return std::make_tuple(std::move(stack),
+ std::make_unique<setup::StackView>(particle));
}
}
--
GitLab