From 92d0d8f0078bf5ea8f2ddce7a6655410824a5b2a Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sat, 24 Apr 2021 21:25:38 +0200
Subject: [PATCH] also re-enable all tests and examples
---
corsika/detail/modules/TrackWriter.inl | 21 ++---
corsika/detail/output/ParquetStreamer.inl | 6 +-
corsika/output/NoOutput.hpp | 59 ++++++++++++
examples/vertical_EAS.cpp | 104 ++++++++++++----------
tests/modules/testObservationPlane.cpp | 10 ++-
5 files changed, 134 insertions(+), 66 deletions(-)
create mode 100644 corsika/output/NoOutput.hpp
diff --git a/corsika/detail/modules/TrackWriter.inl b/corsika/detail/modules/TrackWriter.inl
index f77e763ac..8af3c31dc 100644
--- a/corsika/detail/modules/TrackWriter.inl
+++ b/corsika/detail/modules/TrackWriter.inl
@@ -15,13 +15,14 @@
namespace corsika {
- template <typename TOutput>
- inline TrackWriter<TOutput>::TrackWriter() {}
+ template <typename TOutputWriter>
+ inline TrackWriter<TOutputWriter>::TrackWriter() {}
- template <typename TOutput>
+ template <typename TOutputWriter>
template <typename TParticle, typename TTrack>
- inline ProcessReturn TrackWriter<TOutput>::doContinuous(TParticle const& vP,
- TTrack const& vT, bool const) {
+ inline ProcessReturn TrackWriter<TOutputWriter>::doContinuous(TParticle const& vP,
+ TTrack const& vT,
+ bool const) {
auto const start = vT.getPosition(0).getCoordinates();
auto const end = vT.getPosition(1).getCoordinates();
@@ -32,15 +33,15 @@ namespace corsika {
return ProcessReturn::Ok;
}
- template <typename TOutput>
+ template <typename TOutputWriter>
template <typename TParticle, typename TTrack>
- inline LengthType TrackWriter<TOutput>::getMaxStepLength(TParticle const&,
- TTrack const&) {
+ inline LengthType TrackWriter<TOutputWriter>::getMaxStepLength(TParticle const&,
+ TTrack const&) {
return meter * std::numeric_limits<double>::infinity();
}
- template <typename TOutput>
- YAML::Node TrackWriter<TOutput>::getConfig() const {
+ template <typename TOutputWriter>
+ YAML::Node TrackWriter<TOutputWriter>::getConfig() const {
using namespace units::si;
YAML::Node node;
diff --git a/corsika/detail/output/ParquetStreamer.inl b/corsika/detail/output/ParquetStreamer.inl
index 8f29e3a4e..0a128f992 100644
--- a/corsika/detail/output/ParquetStreamer.inl
+++ b/corsika/detail/output/ParquetStreamer.inl
@@ -30,9 +30,9 @@ namespace corsika {
fields_.push_back(parquet::schema::PrimitiveNode::Make(args...));
}
- void ParquetStreamer::enableCompression(int const level) {
- builder_.compression(parquet::Compression::ZSTD);
- builder_.compression_level(level);
+ void ParquetStreamer::enableCompression(int const /*level*/) {
+ // builder_.compression(parquet::Compression::ZSTD);
+ // builder_.compression_level(level);
}
void ParquetStreamer::buildStreamer() {
diff --git a/corsika/output/NoOutput.hpp b/corsika/output/NoOutput.hpp
new file mode 100644
index 000000000..67710e2d0
--- /dev/null
+++ b/corsika/output/NoOutput.hpp
@@ -0,0 +1,59 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+namespace corsika {
+
+ /**
+ * This is the base class for all outputs so that they
+ * can be stored in homogeneous containers.
+ */
+ class NoOutput : public BaseOutput {
+
+ protected:
+ NoOutput() {}
+
+ public:
+ /**
+ * Called at the start of each run.
+ */
+ void startOfLibrary(boost::filesystem::path const&) final override {}
+
+ /**
+ * Called at the start of each event/shower.
+ */
+ void startOfShower() final override {}
+
+ /**
+ * Called at the end of each event/shower.
+ */
+ void endOfShower() final override {}
+
+ /**
+ * Called at the end of each run.
+ */
+ void endOfLibrary() final override {}
+
+ /**
+ * Get the configuration of this output.
+ */
+ YAML::Node getConfig() const { return YAML::Node(); };
+
+ /**
+ * Get any summary information for the entire library.
+ */
+ YAML::Node getSummary() final override { return YAML::Node(); };
+
+ protected:
+ void write(Code const&, units::si::HEPEnergyType const&, units::si::LengthType const&,
+ units::si::LengthType const&) {}
+ };
+
+} // namespace corsika
+
+#include <corsika/detail/output/BaseOutput.inl>
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 3ada2b875..fa08acc5e 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -136,6 +136,7 @@ int main(int argc, char** argv) {
{{Code::Nitrogen, Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+ builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 2_km);
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
@@ -143,22 +144,6 @@ int main(int argc, char** argv) {
builder.addLinearLayer(1e9_cm, 112.8_km + constants::EarthRadius::Mean);
builder.assemble(env);
- CORSIKA_LOG_DEBUG(
- "environment setup: universe={}, layer1={}, layer2={}, layer3={}, layer4={}, "
- "layer5={}",
- fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 130_km, 0_m, 0_m}))),
- fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 110_km, 0_m, 0_m}))),
- fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 50_km, 0_m, 0_m}))),
- fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 20_km, 0_m, 0_m}))),
- fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 5_km, 0_m, 0_m}))),
- fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 2_km, 0_m, 0_m}))));
-
// pre-setup particle stack
unsigned short const A = std::stoi(std::string(argv[1]));
Code beamCode;
@@ -216,6 +201,9 @@ int main(int argc, char** argv) {
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
+ // create the output manager that we then register outputs with
+ OutputManager output("vertical_EAS_outputs");
+
// setup processes, decays and interactions
// corsika::qgsjetII::Interaction qgsjet;
@@ -282,8 +270,6 @@ int main(int argc, char** argv) {
string const labHist_file = "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz";
string const cMSHist_file = "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz";
string const longprof_file = "longprof_verticalEAS_" + to_string(i_shower) + ".txt";
- string const tracks_file = "tracks_" + to_string(i_shower) + ".dat";
- string const particles_file = "particles_" + to_string(i_shower) + ".dat";
std::cout << std::endl;
std::cout << "Shower " << i_shower << "/" << number_showers << std::endl;
@@ -291,9 +277,50 @@ int main(int argc, char** argv) {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.clear();
+ unsigned short const A = std::stoi(std::string(argv[1]));
+ Code beamCode;
+ HEPEnergyType mass;
+ unsigned short Z = 0;
+ if (A > 0) {
+ beamCode = Code::Nucleus;
+ Z = std::stoi(std::string(argv[2]));
+ mass = get_nucleus_mass(A, Z);
+ } else {
+ int pdg = std::stoi(std::string(argv[2]));
+ beamCode = convert_from_PDG(PDGCode(pdg));
+ mass = get_mass(beamCode);
+ }
+ HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
+ double theta = 0.;
+ auto const thetaRad = theta / 180. * M_PI;
+
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
+ };
+
+ auto const [px, py, pz] = momentumComponents(thetaRad, P0);
+ auto plab = MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << endl;
+ cout << "input momentum: " << plab.getComponents() / 1_GeV
+ << ", norm = " << plab.getNorm() << endl;
+
+ auto const observationHeight = 0_km + builder.getEarthRadius();
+ auto const injectionHeight = 111.75_km + builder.getEarthRadius();
+ auto const t = (injectionHeight - observationHeight) / cos(thetaRad);
+ Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+ Point const injectionPos =
+ showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
+
+ std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
if (A > 1) {
stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
+
} else {
if (A == 1) {
if (Z == 1) {
@@ -314,10 +341,8 @@ int main(int argc, char** argv) {
std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
<< std::endl;
- ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
-
- // create the output manager that we then register outputs with
- OutputManager output("vertical_EAS_outputs");
+ ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env,
+ false};
// setup processes, decays and interactions
@@ -353,10 +378,11 @@ int main(int argc, char** argv) {
decaySibyll.printDecayConfig();
- ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
- corsika::proposal::Interaction emCascade(env);
- corsika::proposal::ContinuousProcess emContinuous(env);
- InteractionCounter emCascadeCounted(emCascade);
+ ParticleCut cut{60_GeV, true, true};
+ // corsika::proposal::Interaction emCascade(env);
+ // corsika::proposal::ContinuousProcess emContinuous(env);
+ // InteractionCounter emCascadeCounted(emCascade);
+ BetheBlochPDG emContinuous(showerAxis);
OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
TrackWriter trackWriter;
@@ -369,26 +395,6 @@ int main(int argc, char** argv) {
// register the observation plane with the output
output.add("obsplane", observationLevel);
- corsika::urqmd::UrQMD urqmd;
- InteractionCounter urqmdCounted{urqmd};
- StackInspector<setup::Stack> stackInspect(50000, false, E0);
-
- // assemble all processes into an ordered process list
- struct EnergySwitch {
- HEPEnergyType cutE_;
- EnergySwitch(HEPEnergyType cutE)
- : cutE_(cutE) {}
- SwitchResult operator()(const Particle& p) {
- if (p.getEnergy() < cutE_)
- return SwitchResult::First;
- else
- return SwitchResult::Second;
- }
- };
- auto hadronSequence =
- make_select(urqmdCounted, make_sequence(sibyllNucCounted, sibyllCounted),
- EnergySwitch(55_GeV));
- auto decaySequence = make_sequence(decayPythia, decaySibyll);
auto sequence =
make_sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
emContinuous, cut, trackWriter, observationLevel, longprof);
@@ -403,16 +409,16 @@ int main(int argc, char** argv) {
EAS.run();
cut.showResults();
- emContinuous.showResults();
+ // emContinuous.showResults();
observationLevel.showResults();
const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
- cut.getEmEnergy() + emContinuous.getEnergyLost() +
+ cut.getEmEnergy() + // emContinuous.getEnergyLost() +
observationLevel.getEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.reset();
cut.reset();
- emContinuous.reset();
+ // emContinuous.reset();
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
diff --git a/tests/modules/testObservationPlane.cpp b/tests/modules/testObservationPlane.cpp
index 1c0b65fc8..03de5d47d 100644
--- a/tests/modules/testObservationPlane.cpp
+++ b/tests/modules/testObservationPlane.cpp
@@ -17,6 +17,8 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/output/NoOutput.hpp>
+
#include <SetupTestEnvironment.hpp>
#include <SetupTestStack.hpp>
#include <SetupTestTrajectory.hpp>
@@ -54,7 +56,7 @@ TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") {
SECTION("horizontal plane") {
Plane const obsPlane(Point(cs, {10_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.}));
- ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 1., 0.}));
+ ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 1., 0.}));
LengthType const length = obs.getMaxStepLength(particle, no_used_track);
ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true);
@@ -65,10 +67,10 @@ TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") {
SECTION("transparent plane") {
Plane const obsPlane(Point(cs, {1_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.}));
- ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 0., 1.}));
+ ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 0., 1.}));
LengthType const length = obs.getMaxStepLength(particle, no_used_track);
- ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true);
+ ProcessReturn const ret = obs.doContinuous(particle, no_used_track, false);
CHECK(length / 1_m == Approx(1).margin(1e-4));
CHECK(ret == ProcessReturn::Ok);
@@ -83,7 +85,7 @@ TEST_CASE("ObservationPlane", "[proccesses][observation_plane]") {
Plane const obsPlane(Point(cs, {10_m, 5_m, 5_m}),
DirectionVector(cs, {1, 0.1, -0.05}));
- ObservationPlane obs(obsPlane, DirectionVector(cs, {0., 1., 0.}));
+ ObservationPlane<NoOutput> obs(obsPlane, DirectionVector(cs, {0., 1., 0.}));
LengthType const length = obs.getMaxStepLength(particle, no_used_track);
ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true);
--
GitLab