diff --git a/.clang-format b/.clang-format
index 0ac489959f851bf8930e40fc9eec547d23c4909b..afc401b287b98713c3647a8a4b1752b3c9c7c080 100644
--- a/.clang-format
+++ b/.clang-format
@@ -57,6 +57,7 @@ IncludeCategories:
Priority: 2
- Regex: '.*'
Priority: 3
+SortIncludes: false
IndentCaseLabels: true
IndentWidth: 2
IndentWrappedFunctionNames: false
@@ -76,7 +77,6 @@ PenaltyExcessCharacter: 1000000
PenaltyReturnTypeOnItsOwnLine: 200
PointerAlignment: Left
ReflowComments: true
-SortIncludes: true
SpaceAfterCStyleCast: false
SpaceBeforeAssignmentOperators: true
SpaceBeforeParens: ControlStatements
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 1bc50e665c00b23c53ef16812db92fba158991c2..464ce2402377ddacbf94a9147aa958ebbd9bdf3d 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -276,6 +276,8 @@ build_test-clang-8:
# normal pipeline for each commit
.example:
stage: example
+ before_script:
+ - apt-get -qq update && apt-get -qq install -y gdb
tags:
- corsika
script:
@@ -330,6 +332,8 @@ example-clang-8:
stage: build_test_example
tags:
- corsika
+ before_script:
+ - apt-get -qq update && apt-get -qq install -y gdb
script:
- cd build
- cmake --build . -- -j4
@@ -439,6 +443,8 @@ install-clang-8:
stage: optional
tags:
- corsika
+ before_script:
+ - apt-get -qq update && apt-get -qq install -y gdb
script:
- cd build
- cmake .. -DCMAKE_BUILD_TYPE=Release
diff --git a/.gitlab/merge_request_templates/Code Approval.md b/.gitlab/merge_request_templates/Code Approval.md
index fed15a38efcfdb1a0a4ded8050733598ec2856ca..ed076c02e7be42adb5da776b21546e985d2d60c5 100644
--- a/.gitlab/merge_request_templates/Code Approval.md
+++ b/.gitlab/merge_request_templates/Code Approval.md
@@ -1,15 +1,18 @@
-The code approval procedure is described in the [code approval procedure wiki](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/wikis/Code-Approval-Procedure)
+Issues: Closes #....
-- [ ] The MR is without WIP (work in progress) status
-- [ ] Make sure the most recent CI jobs (config, quality, build, tests) all run fine with no failures
+The code approval procedure is described in the wiki: [Code approval procedure wiki](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/wikis/Code-Approval-Procedure)
+
+- [ ] The MR is without `WIP/Draft` status
+- [ ] Make sure the most recent CI jobs (config, quality, build_test_example) all run fine with no failures
- if "check clang-format" failed: the code contributor has to run `./do-clang-format.py --apply` eventually with the `--all` option
- if "check copyright" failed the code contributor has to run`./do-copyright.py --add=20xy`
-- [ ] Also run all the extra optional jobs "coverage", "release-clang-8", "release-u-18.04" and make sure no problems occur
-- [ ] Check in the "coverage" job output that the coverage did not decreases (!). It should always stay, or increase. If it decreased --> ask contributor to add further unit tests, and check coverage report.
+- [ ] Make sure also the jobs with MR-label `ready for code review` succeed. This includes the optional jobs, in particular 'coverage', 'release-clang-8", "release-u-18.04" and make sure no problems occur. You may have to trigger a pipeline manually to check this.
+- [ ] Check in the "coverage" job output that the coverage did not decrease. It should always stay, or increase. If it decreased --> ask contributor to add further needed unit tests, and check coverage report.
- On the MR page, open the "Open in Web IDE" tool
- [ ] Check if the provided solution solves the Issue, discuss on gitlab
- [ ] Check that all changes are actually related to the issue
- [ ] There are no debug statements left, not even commented out
- [ ] Check all changes for coding rules and guidelines
-- When all above is done
- - [ ] **Add label "Code Review Finished"**
+- When all above is done:
+ - **Add MR label "Code Review Finished"**
+
diff --git a/CMakeModules/CorsikaUtilities.cmake b/CMakeModules/CorsikaUtilities.cmake
index d80269b3a24e549a7655b769099bd6579a79a112..96d29eeb614391ae0f4fd05a4892e7d02eeeef2d 100644
--- a/CMakeModules/CorsikaUtilities.cmake
+++ b/CMakeModules/CorsikaUtilities.cmake
@@ -174,13 +174,15 @@ endfunction (CORSIKA_ADD_TEST)
# specify the sources.
#
# Example: CORSIKA_ADD_EXAMPLE (testSomething
-# SOURCES source1.cc source2.cc someheader.h)
+# SOURCES source1.cc source2.cc someheader.h
+# RUN_OPTION "extra command line options"
+# RUN_IN_GDB)
#
# In all cases, you can further customize the target with
# target_link_libraries(testSomething ...) and so on.
#
function (CORSIKA_ADD_EXAMPLE)
- cmake_parse_arguments (PARSE_ARGV 1 C8_ADD_EXAMPLE "" "" "SOURCES;RUN_OPTIONS")
+ cmake_parse_arguments (PARSE_ARGV 1 C8_ADD_EXAMPLE "RUN_IN_GDB" "" "SOURCES;RUN_OPTIONS")
set (name ${ARGV0})
if (NOT C8_ADD_EXAMPLE_SOURCES)
@@ -204,12 +206,21 @@ function (CORSIKA_ADD_EXAMPLE)
add_custom_target (run_examples)
endif ()
add_dependencies (run_examples ${name})
+ if (C8_ADD_EXAMPLE_RUN_IN_GDB)
+ # convert cmake list into real string:
+ string (REPLACE ";" " " run_options_str "${run_options}")
+ # run the command in gdb and study backtrace a bit
+ set (CMD gdb -q --batch -ex "run ${run_options_str}" -ex bt -ex "info locals" -ex "up" -ex "info locals" -ex "up" -ex "info locals" -ex "up" -ex "info locals" -ex quit "${CMAKE_CURRENT_BINARY_DIR}/${name}")
+ else (C8_ADD_EXAMPLE_RUN_IN_GDB)
+ # just run the command as-is
+ set (CMD ${CMAKE_CURRENT_BINARY_DIR}/${name} ${run_options})
+ endif (C8_ADD_EXAMPLE_RUN_IN_GDB)
add_custom_command (TARGET run_examples
POST_BUILD
COMMAND ${CMAKE_COMMAND} -E echo ""
COMMAND ${CMAKE_COMMAND} -E echo "**************************************"
COMMAND ${CMAKE_COMMAND} -E echo "***** running example: ${name} " ${run_options} VERBATIM
- COMMAND ${CMAKE_CURRENT_BINARY_DIR}/${name} ${run_options} VERBATIM
+ COMMAND ${CMD} VERBATIM
WORKING_DIRECTORY ${PROJECT_BINARY_DIR}/example_outputs)
install (TARGETS ${name} DESTINATION share/examples)
endfunction (CORSIKA_ADD_EXAMPLE)
diff --git a/Documentation/Doxygen/Doxyfile.in b/Documentation/Doxygen/Doxyfile.in
index 77031e96dc0e9d9f1b17cb4fed4f0e2201deb284..b733d09258adfa7e7200033acda18823a8c8ca9d 100644
--- a/Documentation/Doxygen/Doxyfile.in
+++ b/Documentation/Doxygen/Doxyfile.in
@@ -7,11 +7,11 @@ GENERATE_XML = YES
OUTPUT_DIRECTORY = @CMAKE_CURRENT_BINARY_DIR@/
INPUT = @PROJECT_SOURCE_DIR@ @PROJECT_BINARY_DIR@/Framework
-EXCLUDE_PATTERNS = */ThirdParty/*/*
+EXCLUDE_PATTERNS = */ThirdParty/*/* */build*/corsika/*
FULL_PATH_NAMES = YES
STRIP_FROM_PATH = @PROJECT_SOURCE_DIR@
-FILE_PATTERNS = *.cc *.cpp *.cxx *.h *.dox *.inc *.md
+FILE_PATTERNS = *.cc *.cpp *.cxx *.h *.hpp *.dox *.inc *.md
EXTENSION_MAPPING = inc=C++
RECURSIVE = YES
diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt
index 17c711669f514a311dbaa9297139ca6d26959f8b..62c9d6bd3b9b05e28cf7985a4c9766a7ed4570c4 100644
--- a/Documentation/Examples/CMakeLists.txt
+++ b/Documentation/Examples/CMakeLists.txt
@@ -24,12 +24,12 @@ target_link_libraries (boundary_example
ProcessParticleCut
ProcessTrackingLine
ProcessPythia8
- ProcessProposal
CORSIKAprocesses
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
+ C8::ext::boost # boost::histogram
)
CORSIKA_ADD_EXAMPLE (cascade_example)
@@ -66,7 +66,6 @@ if (Pythia8_FOUND)
ProcessSibyll
ProcessPythia8
ProcessUrQMD
- ProcessSwitch
CORSIKAcascade
ProcessCONEXSourceCut
ProcessEnergyLoss
@@ -95,7 +94,6 @@ if (Pythia8_FOUND)
ProcessSibyll
ProcessPythia8
ProcessUrQMD
- ProcessSwitch
CORSIKAcascade
ProcessProposal
ProcessPythia8
@@ -123,6 +121,7 @@ CORSIKA_ADD_EXAMPLE (stopping_power stopping_power)
target_link_libraries (stopping_power
CORSIKAsetup
CORSIKAunits
+ CORSIKAlogging
ProcessEnergyLoss
CORSIKAparticles
CORSIKAgeometry
@@ -137,7 +136,7 @@ target_link_libraries (staticsequence_example
CORSIKAlogging)
-CORSIKA_ADD_EXAMPLE (em_shower RUN_OPTIONS 100.)
+CORSIKA_ADD_EXAMPLE (em_shower RUN_OPTIONS "100.")
target_link_libraries (em_shower
SuperStupidStack
CORSIKAunits
@@ -146,7 +145,6 @@ target_link_libraries (em_shower
ProcessSibyll
ProcessPythia8
ProcessUrQMD
- ProcessSwitch
CORSIKAcascade
ProcessCONEXSourceCut
ProcessEnergyLoss
diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index 40c60ebed365e24d8e4aed8595968581fa7252e7..9bf4e188ae18a1e668e86c9e8294d52a1c9642b1 100644
--- a/Documentation/Examples/boundary_example.cc
+++ b/Documentation/Examples/boundary_example.cc
@@ -34,6 +34,8 @@
#include <corsika/utl/CorsikaFenv.h>
+#include <corsika/logging/Logging.h>
+
#include <iostream>
#include <limits>
#include <typeinfo>
@@ -60,7 +62,9 @@ struct MyBoundaryCrossingProcess
EProcessReturn DoBoundaryCrossing(Particle& p,
typename Particle::BaseNodeType const& from,
typename Particle::BaseNodeType const& to) {
- std::cout << "boundary crossing! from: " << &from << "; to: " << &to << std::endl;
+
+ C8LOG_INFO("MyBoundaryCrossingProcess: crossing! from: {} to: {} ", fmt::ptr(&from),
+ fmt::ptr(&to));
auto const& name = particles::GetName(p.GetPID());
auto const start = p.GetPosition().GetCoordinates();
@@ -82,7 +86,9 @@ private:
//
int main() {
- std::cout << "boundary_example" << std::endl;
+ logging::SetLevel(logging::level::info);
+
+ C8LOG_INFO("boundary_example");
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
@@ -95,24 +101,23 @@ int main() {
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
- auto outerMedium = EnvType::CreateNode<Sphere>(
+ // create "world" as infinite sphere filled with protons
+ auto world = EnvType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
auto const props =
- outerMedium
- ->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
- 1_kg / (1_m * 1_m * 1_m),
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Proton},
- std::vector<float>{1.f}));
-
- auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5_km);
-
- innerMedium->SetModelProperties(props);
+ world->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
+ 1_kg / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(
+ std::vector<particles::Code>{particles::Code::Proton},
+ std::vector<float>{1.f}));
- outerMedium->AddChild(std::move(innerMedium));
+ // add a "target" sphere with 5km readius at 0,0,0
+ auto target = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5_km);
+ target->SetModelProperties(props);
- universe.AddChild(std::move(outerMedium));
+ world->AddChild(std::move(target));
+ universe.AddChild(std::move(world));
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
@@ -121,28 +126,28 @@ int main() {
process::sibyll::Interaction sibyll;
process::sibyll::Decay decay;
- process::particle_cut::ParticleCut cut(20_GeV, true, true);
+ process::particle_cut::ParticleCut cut(50_GeV, true, true);
- process::track_writer::TrackWriter trackWriter("tracks.dat");
+ process::track_writer::TrackWriter trackWriter("boundary_tracks.dat");
MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
// assemble all processes into an ordered process list
- auto sequence = sibyll << decay << cut << boundaryCrossing << trackWriter;
+ auto sequence = process::sequence(sibyll, decay, cut, boundaryCrossing, trackWriter);
// setup particle stack, and add primary particles
setup::Stack stack;
stack.Clear();
- const Code beamCode = Code::Proton;
- const HEPMassType mass = particles::GetMass(Code::Proton);
- const HEPEnergyType E0 = 50_TeV;
+ const Code beamCode = Code::MuPlus;
+ const HEPMassType mass = particles::GetMass(beamCode);
+ const HEPEnergyType E0 = 100_GeV;
std::uniform_real_distribution distTheta(0., 180.);
std::uniform_real_distribution distPhi(0., 360.);
std::mt19937 rng;
for (int i = 0; i < 100; ++i) {
- auto const theta = distTheta(rng);
- auto const phi = distPhi(rng);
+ double const theta = distTheta(rng);
+ double const phi = distPhi(rng);
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
@@ -155,9 +160,12 @@ int main() {
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
- cout << "input particle: " << beamCode << endl;
- cout << "input angles: theta=" << theta << " phi=" << phi << endl;
- cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
+ C8LOG_INFO(
+ "input particle: {} "
+ "input angles: theta={} phi={}"
+ "input momentum: {} GeV",
+ beamCode, theta, phi, plab.GetComponents() / 1_GeV);
+ // shoot particles from inside target out
Point pos(rootCS, 0_m, 0_m, 0_m);
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
@@ -170,10 +178,10 @@ int main() {
EAS.Run();
- cout << "Result: E0=" << E0 / 1_GeV << endl;
+ C8LOG_INFO("Result: E0={}GeV", E0 / 1_GeV);
cut.ShowResults();
- const HEPEnergyType Efinal =
- cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
- cout << "total energy (GeV): " << Efinal / 1_GeV << endl
- << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
+ [[maybe_unused]] const HEPEnergyType Efinal =
+ (cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy());
+ C8LOG_INFO("Total energy (GeV): {} relative difference (%): {}", Efinal / 1_GeV,
+ (Efinal / E0 - 1.) * 100);
}
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index fe8e05d87460d23d3517f4028a686198205dff30..64ed627e03c8c7687a2b75330dd50d46e1168627 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -35,6 +35,7 @@
#include <corsika/random/RNGManager.h>
#include <corsika/utl/CorsikaFenv.h>
+#include <corsika/logging/Logging.h>
#include <iostream>
#include <limits>
@@ -56,6 +57,8 @@ using namespace corsika::units::si;
//
int main() {
+ logging::SetLevel(logging::level::debug);
+
std::cout << "cascade_example" << std::endl;
const LengthType height_atmosphere = 112.8_km;
@@ -147,8 +150,8 @@ int main() {
process::energy_loss::EnergyLoss eLoss{showerAxis, cut.GetECut()};
// assemble all processes into an ordered process list
- auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut
- << trackWriter;
+ auto sequence =
+ process::sequence(stackInspect, sibyll, sibyllNuc, decay, eLoss, cut, trackWriter);
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
diff --git a/Documentation/Examples/cascade_proton_example.cc b/Documentation/Examples/cascade_proton_example.cc
index 6050f7fcfae52c7a5f77e76928abbf44113c276d..4ca053f377d67d2aa597b8e377575c4f960cd305 100644
--- a/Documentation/Examples/cascade_proton_example.cc
+++ b/Documentation/Examples/cascade_proton_example.cc
@@ -38,6 +38,8 @@
#include <corsika/utl/CorsikaFenv.h>
+#include <corsika/logging/Logging.h>
+
#include <iostream>
#include <limits>
#include <typeinfo>
@@ -59,6 +61,8 @@ using namespace corsika::units::si;
//
int main() {
+ logging::SetLevel(logging::level::info);
+
std::cout << "cascade_proton_example" << std::endl;
feenableexcept(FE_INVALID);
@@ -136,7 +140,7 @@ int main() {
// assemble all processes into an ordered process list
// auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;
- auto sequence = pythia << decay << cut << trackWriter << stackInspect;
+ auto sequence = process::sequence(pythia, decay, cut, trackWriter, stackInspect);
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
diff --git a/Documentation/Examples/em_shower.cc b/Documentation/Examples/em_shower.cc
index 1c5c07cd483242cceed18aa356f4167a0253abb5..c41d4574755db08489665f1bfde06fcded43d4c0 100644
--- a/Documentation/Examples/em_shower.cc
+++ b/Documentation/Examples/em_shower.cc
@@ -29,6 +29,8 @@
#include <corsika/utl/CorsikaFenv.h>
#include <corsika/process/interaction_counter/InteractionCounter.hpp>
+#include <corsika/logging/Logging.h>
+
#include <iomanip>
#include <iostream>
#include <limits>
@@ -54,6 +56,9 @@ void registerRandomStreams() {
}
int main(int argc, char** argv) {
+
+ logging::SetLevel(logging::level::info);
+
if (argc != 2) {
std::cerr << "usage: em_shower <energy/GeV>" << std::endl;
return 1;
@@ -144,8 +149,8 @@ int main(int argc, char** argv) {
process::observation_plane::ObservationPlane observationLevel(obsPlane,
"particles.dat");
- auto sequence = proposalCounted << em_continuous << longprof << cut << observationLevel
- << trackWriter;
+ auto sequence = process::sequence(proposalCounted, em_continuous, longprof, cut,
+ observationLevel, trackWriter);
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
@@ -169,11 +174,7 @@ int main(int argc, char** argv) {
em_continuous.Reset();
auto const hists = proposalCounted.GetHistogram();
- hists.saveLab("inthist_lab.txt");
- hists.saveCMS("inthist_cms.txt");
-
- longprof.save("longprof.txt");
-
- std::ofstream finish("finished");
- finish << "run completed without error" << std::endl;
+ hists.saveLab("inthist_lab_emShower.npz");
+ hists.saveCMS("inthist_cms_emShower.npz");
+ longprof.save("longprof_emShower.txt");
}
diff --git a/Documentation/Examples/staticsequence_example.cc b/Documentation/Examples/staticsequence_example.cc
index ca75fe4b22172488675fddca7601748acf8e3214..659a3d6c2d13e9f3840f80d162cb8dff3feadfac 100644
--- a/Documentation/Examples/staticsequence_example.cc
+++ b/Documentation/Examples/staticsequence_example.cc
@@ -23,43 +23,43 @@ using namespace std;
const int nData = 10;
-class Process1 : public BaseProcess<Process1> {
+class Process1 : public ContinuousProcess<Process1> {
public:
Process1() {}
- template <typename D, typename T, typename S>
- EProcessReturn DoContinuous(D& d, T&, S&) const {
+ template <typename D, typename T>
+ EProcessReturn DoContinuous(D& d, T&) const {
for (int i = 0; i < nData; ++i) d.p[i] += 1;
return EProcessReturn::eOk;
}
};
-class Process2 : public BaseProcess<Process2> {
+class Process2 : public ContinuousProcess<Process2> {
public:
Process2() {}
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+ template <typename D, typename T>
+ inline EProcessReturn DoContinuous(D& d, T&) const {
for (int i = 0; i < nData; ++i) d.p[i] -= 0.1 * i;
return EProcessReturn::eOk;
}
};
-class Process3 : public BaseProcess<Process3> {
+class Process3 : public ContinuousProcess<Process3> {
public:
Process3() {}
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D&, T&, S&) const {
+ template <typename D, typename T>
+ inline EProcessReturn DoContinuous(D&, T&) const {
return EProcessReturn::eOk;
}
};
-class Process4 : public BaseProcess<Process4> {
+class Process4 : public ContinuousProcess<Process4> {
public:
Process4(const double v)
: fV(v) {}
- template <typename D, typename T, typename S>
- inline EProcessReturn DoContinuous(D& d, T&, S&) const {
+ template <typename D, typename T>
+ inline EProcessReturn DoContinuous(D& d, T&) const {
for (int i = 0; i < nData; ++i) d.p[i] *= fV;
return EProcessReturn::eOk;
}
@@ -78,25 +78,31 @@ struct DummyTrajectory {};
void modular() {
- Process1 m1;
- Process2 m2;
- Process3 m3;
- Process4 m4(0.9);
+ // = 0
+ Process1 m1; // + 1.0
+ Process2 m2; // - (0.1*i)
+ Process3 m3; // * 1.0
+ Process4 m4(1.5); // * 1.5
- auto sequence = m1 << m2 << m3 << m4;
+ auto sequence = process::sequence(m1, m2, m3, m4);
DummyData particle;
DummyTrajectory track;
- const int n = 1000;
- for (int i = 0; i < n; ++i) { sequence.DoContinuous(particle, track); }
+ double check[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
- for (int i = 0; i < nData; ++i) {
- // cout << p.p[i] << endl;
- // assert(p.p[i] == n-i*100);
+ const int nEv = 10;
+ for (int iEv = 0; iEv < nEv; ++iEv) {
+ sequence.DoContinuous(particle, track);
+ for (int i = 0; i < nData; ++i) {
+ check[i] += 1. - 0.1 * i;
+ check[i] *= 1.5;
+ }
}
- cout << " done (nothing...) " << endl;
+ for (int i = 0; i < nData; ++i) { assert(particle.p[i] == check[i]); }
+
+ cout << " done (checking...) " << endl;
}
int main() {
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 042e64bc29f468afa32c45b654b9082adedc33dd..14eeea8b94560f4f7233c8d385606702fd5e6294 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -22,6 +22,7 @@
#include <corsika/geometry/Sphere.h>
#include <corsika/logging/Logging.h>
#include <corsika/process/ProcessSequence.h>
+#include <corsika/process/SwitchProcessSequence.h>
#include <corsika/process/StackProcess.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
@@ -34,7 +35,6 @@
#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/tracking_line/TrackingLine.h>
#include <corsika/process/urqmd/UrQMD.h>
#include <corsika/random/RNGManager.h>
@@ -60,6 +60,8 @@ using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
+using Particle = setup::Stack::StackIterator;
+
void registerRandomStreams(const int seed) {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
random::RNGManager::GetInstance().RegisterRandomStream("qgsjet");
@@ -149,23 +151,20 @@ int main(int argc, char** argv) {
std::cout << "point of injection: " << injectionPos.GetCoordinates() << std::endl;
if (A != 1) {
- stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType, unsigned short, unsigned short>{
- beamCode, E0, plab, injectionPos, 0_ns, A, Z});
+ stack.AddParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
} else {
stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::Proton, E0, plab, injectionPos, 0_ns});
+ std::make_tuple(particles::Code::Proton, E0, plab, injectionPos, 0_ns));
}
- std::cout << "shower axis length: " << (showerCore - injectionPos).norm() * 1.02
+ // we make the axis much longer than the inj-core distance since the
+ // profile will go beyond the core, depending on zenith angle
+ std::cout << "shower axis length: " << (showerCore - injectionPos).norm() * 1.5
<< std::endl;
environment::ShowerAxis const showerAxis{injectionPos,
- (showerCore - injectionPos) * 1.02, env};
+ (showerCore - injectionPos) * 1.5, env};
// setup processes, decays and interactions
@@ -200,7 +199,6 @@ int main(int argc, char** argv) {
decaySibyll.PrintDecayConfig();
- // PROPOSAL processs proposal{...};
process::particle_cut::ParticleCut cut{60_GeV, false, true};
process::proposal::Interaction proposal(env, cut.GetECut());
process::proposal::ContinuousProcess em_continuous(env, cut.GetECut());
@@ -218,21 +216,30 @@ int main(int argc, char** argv) {
process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
// assemble all processes into an ordered process list
-
- auto sibyllSequence = sibyllNucCounted << sibyllCounted;
- process::switch_process::SwitchProcess switchProcess(urqmdCounted, sibyllSequence,
- 55_GeV);
- auto decaySequence = decayPythia << decaySibyll;
-
- auto sequence = switchProcess << reset_particle_mass << decaySequence << proposalCounted
- << em_continuous << cut << longprof << observationLevel;
+ struct EnergySwitch {
+ HEPEnergyType cutE_;
+ EnergySwitch(HEPEnergyType cutE)
+ : cutE_(cutE) {}
+ process::SwitchResult operator()(const Particle& p) {
+ if (p.GetEnergy() < cutE_)
+ return process::SwitchResult::First;
+ else
+ return process::SwitchResult::Second;
+ }
+ };
+ auto hadronSequence =
+ process::select(urqmdCounted, process::sequence(sibyllNucCounted, sibyllCounted),
+ EnergySwitch(55_GeV));
+ auto decaySequence = process::sequence(decayPythia, decaySibyll);
+ auto sequence =
+ process::sequence(hadronSequence, reset_particle_mass, decaySequence,
+ proposalCounted, em_continuous, cut, observationLevel, longprof);
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
// to fix the point of first interaction, uncomment the following two lines:
- // EAS.SetNodes();
// EAS.forceInteraction();
EAS.Run();
@@ -252,11 +259,7 @@ int main(int argc, char** argv) {
auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
urqmdCounted.GetHistogram() + proposalCounted.GetHistogram();
- hists.saveLab("inthist_lab.npz");
- hists.saveCMS("inthist_cms.npz");
-
- longprof.save("longprof.txt");
-
- std::ofstream finish("finished");
- finish << "run completed without error" << std::endl;
+ hists.saveLab("inthist_lab_verticalEAS.npz");
+ hists.saveCMS("inthist_cms_verticalEAS.npz");
+ longprof.save("longprof_verticalEAS.txt");
}
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index 5536d2ea3e9aee1c2f3df00ed4b3b7db81d6e138..f0adcf74a8fcbd55abf8eaec6b75a42700630f6d 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -49,7 +49,7 @@ target_link_libraries (
CORSIKAgeometry
CORSIKAparticles
CORSIKAunits
- CORSIKArandom
+ CORSIKAlogging
)
target_include_directories (
diff --git a/Environment/ShowerAxis.cc b/Environment/ShowerAxis.cc
index fadc86ab352b5285021b4c6dfe4006e34ddab3a9..3cc5779e08a6c411a39a7d5aeec0831f5c1e4f2b 100644
--- a/Environment/ShowerAxis.cc
+++ b/Environment/ShowerAxis.cc
@@ -7,7 +7,9 @@
*/
#include <corsika/environment/ShowerAxis.h>
-#include <sstream>
+#include <corsika/logging/Logging.h>
+
+#include <string>
using namespace corsika::environment;
using namespace corsika::units::si;
@@ -20,19 +22,22 @@ GrammageType ShowerAxis::X(LengthType l) const {
decltype(X_.size()) const upper = lower + 1;
if (lower < 0) {
+ C8LOG_ERROR("cannot extrapolate to points behind point of injection l={} m", l / 1_m);
throw std::runtime_error("cannot extrapolate to points behind point of injection");
}
if (upper >= X_.size()) {
- std::stringstream errormsg;
- errormsg << "shower axis too short, cannot extrapolate (l / max_length_ = "
- << l / max_length_ << ")";
- throw std::runtime_error(errormsg.str().c_str());
+ const std::string err =
+ fmt::format("shower axis too short, cannot extrapolate (l / max_length_ = {} )",
+ l / max_length_);
+ C8LOG_ERROR(err);
+ throw std::runtime_error(err.c_str());
}
assert(0 <= lambda && lambda <= 1.);
- std::cout << l << ": " << lower << " " << lambda << " " << upper << std::endl;
+ C8LOG_TRACE("ShowerAxis::X l={} m, lower={}, lambda={}, upper={}", l / 1_m, lower,
+ lambda, upper);
// linear interpolation between X[lower] and X[upper]
return X_[upper] * lambda + X_[lower] * (1 - lambda);
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/Cascade/CMakeLists.txt b/Framework/Cascade/CMakeLists.txt
index 9bb6019340dc8578c4f9016091e0cc2c1de0162f..2a12099736e967a7bbcfc9db3395f46fc94590cd 100644
--- a/Framework/Cascade/CMakeLists.txt
+++ b/Framework/Cascade/CMakeLists.txt
@@ -51,6 +51,5 @@ target_link_libraries (
CORSIKAcascade
ProcessStackInspector
ProcessTrackingLine
- ProcessNullModel
CORSIKAtesting
)
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index b520acddde752d6e4b0db7d3cdf711f154b1689c..43a5275435480e45aaa3bb06ac824f2ba5636f40 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -74,11 +74,11 @@ namespace corsika::cascade {
private:
// Data members
- corsika::environment::Environment<MediumInterface> const& fEnvironment;
- TTracking& fTracking;
- TProcessList& fProcessSequence;
- TStack& fStack;
- corsika::random::RNG& fRNG =
+ corsika::environment::Environment<MediumInterface> const& environment_;
+ TTracking& tracking_;
+ TProcessList& process_sequence_;
+ TStack& stack_;
+ corsika::random::RNG& rng_ =
corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
unsigned int count_ = 0;
@@ -93,27 +93,15 @@ namespace corsika::cascade {
*/
Cascade(corsika::environment::Environment<MediumInterface> const& env, TTracking& tr,
TProcessList& pl, TStack& stack)
- : fEnvironment(env)
- , fTracking(tr)
- , fProcessSequence(pl)
- , fStack(stack)
+ : environment_(env)
+ , tracking_(tr)
+ , process_sequence_(pl)
+ , stack_(stack)
, count_(0) {
C8LOG_INFO(c8_ascii_);
-#ifdef WITH_HISTORY
- C8LOG_INFO(" - With full cascade HISTORY.");
-#endif
- }
-
- /**
- * set the nodes for all particles on the stack according to their numerical
- * position
- */
- void SetNodes() {
- std::for_each(fStack.begin(), fStack.end(), [&](auto& p) {
- auto const* numericalNode =
- fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
- p.SetNode(numericalNode);
- });
+ if constexpr (TStackView::has_event) {
+ C8LOG_INFO(" - With full cascade HISTORY.");
+ }
}
/**
@@ -121,20 +109,20 @@ namespace corsika::cascade {
* particles from the Stack until the Stack is empty.
*/
void Run() {
- SetNodes();
+ setNodes();
- while (!fStack.IsEmpty()) {
- while (!fStack.IsEmpty()) {
- C8LOG_TRACE(fmt::format("Stack: {}", fStack.as_string()));
+ while (!stack_.IsEmpty()) {
+ while (!stack_.IsEmpty()) {
+ C8LOG_TRACE("Stack: {}", stack_.as_string());
count_++;
- auto pNext = fStack.GetNextParticle();
- C8LOG_DEBUG(fmt::format(
+ auto pNext = stack_.GetNextParticle();
+ C8LOG_DEBUG(
"============== next particle : count={}, pid={}, "
", stack entries={}"
", stack deleted={}",
- count_, pNext.GetPID(), fStack.getEntries(), fStack.getDeleted()));
+ count_, pNext.GetPID(), stack_.getEntries(), stack_.getDeleted());
Step(pNext);
- fProcessSequence.DoStack(fStack);
+ process_sequence_.DoStack(stack_);
}
// do cascade equations, which can put new particles on Stack,
// thus, the double loop
@@ -149,11 +137,12 @@ namespace corsika::cascade {
*/
void forceInteraction() {
C8LOG_DEBUG("forced interaction!");
- auto vParticle = fStack.GetNextParticle();
+ setNodes();
+ auto vParticle = stack_.GetNextParticle();
TStackView secondaries(vParticle);
- interaction(vParticle, secondaries);
- fProcessSequence.DoSecondaries(secondaries);
- vParticle.Delete(); // todo: this should be reviewed, see below
+ interaction(secondaries);
+ process_sequence_.DoSecondaries(secondaries);
+ vParticle.Delete(); // primary particle has interacted and is gone
}
private:
@@ -172,16 +161,16 @@ namespace corsika::cascade {
using namespace corsika::units::si;
// determine geometric tracking
- auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(vParticle);
+ auto [step, geomMaxLength, nextVol] = tracking_.GetTrack(vParticle);
[[maybe_unused]] auto const& dummy_nextVol = nextVol;
// determine combined total interaction length (inverse)
InverseGrammageType const total_inv_lambda =
- fProcessSequence.GetTotalInverseInteractionLength(vParticle);
+ process_sequence_.GetInverseInteractionLength(vParticle);
// sample random exponential step length in grammage
corsika::random::ExponentialDistribution expDist(1 / total_inv_lambda);
- GrammageType const next_interact = expDist(fRNG);
+ GrammageType const next_interact = expDist(rng_);
C8LOG_DEBUG(
"total_lambda={} g/cm2, "
@@ -193,8 +182,8 @@ namespace corsika::cascade {
// assert that particle stays outside void Universe if it has no
// model properties set
- assert(currentLogicalNode != &*fEnvironment.GetUniverse() ||
- fEnvironment.GetUniverse()->HasModelProperties());
+ assert(currentLogicalNode != &*environment_.GetUniverse() ||
+ environment_.GetUniverse()->HasModelProperties());
// convert next_step from grammage to length
LengthType const distance_interact =
@@ -202,16 +191,16 @@ namespace corsika::cascade {
next_interact);
// determine the maximum geometric step length from continuous processes
- LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
+ LengthType const distance_max = process_sequence_.MaxStepLength(vParticle, step);
C8LOG_DEBUG("distance_max={} m", distance_max / 1_m);
// determine combined total inverse decay time
InverseTimeType const total_inv_lifetime =
- fProcessSequence.GetTotalInverseLifetime(vParticle);
+ process_sequence_.GetInverseLifetime(vParticle);
// sample random exponential decay time
corsika::random::ExponentialDistribution expDistDecay(1 / total_inv_lifetime);
- TimeType const next_decay = expDistDecay(fRNG);
+ TimeType const next_decay = expDistDecay(rng_);
C8LOG_DEBUG(
"total_lifetime={} s"
", next_decay={} s",
@@ -236,12 +225,11 @@ namespace corsika::cascade {
step.LimitEndTo(min_distance);
// apply all continuous processes on particle + track
- process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, step);
-
- if (status == process::EProcessReturn::eParticleAbsorbed) {
+ if (process_sequence_.DoContinuous(vParticle, step) ==
+ process::EProcessReturn::eParticleAbsorbed) {
C8LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
vParticle.GetPID(), vParticle.GetEnergy() / 1_GeV);
- vParticle.Delete();
+ if (!vParticle.isDeleted()) vParticle.Delete();
return;
}
@@ -271,10 +259,10 @@ namespace corsika::cascade {
[[maybe_unused]] auto projectile = secondaries.GetProjectile();
if (min_distance == distance_interact) {
- interaction(vParticle, secondaries);
+ interaction(secondaries);
} else {
assert(min_distance == distance_decay);
- decay(vParticle, secondaries);
+ decay(secondaries);
// make sure particle actually did decay if it should have done so
if (secondaries.getSize() == 1 &&
projectile.GetPID() == secondaries.GetNextParticle().GetPID())
@@ -283,7 +271,7 @@ namespace corsika::cascade {
particles::GetName(projectile.GetPID())));
}
- fProcessSequence.DoSecondaries(secondaries);
+ process_sequence_.DoSecondaries(secondaries);
vParticle.Delete();
} else { // step-length limitation within volume
@@ -293,10 +281,9 @@ namespace corsika::cascade {
[[maybe_unused]] auto const assertion = [&] {
auto const* numericalNodeAfterStep =
- fEnvironment.GetUniverse()->GetContainingNode(vParticle.GetPosition());
- C8LOG_TRACE(fmt::format(
- "Geometry check: numericalNodeAfterStep={} currentLogicalNode={}",
- fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode)));
+ environment_.GetUniverse()->GetContainingNode(vParticle.GetPosition());
+ C8LOG_TRACE("Geometry check: numericalNodeAfterStep={} currentLogicalNode={}",
+ fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode));
return numericalNodeAfterStep == currentLogicalNode;
};
@@ -306,45 +293,61 @@ namespace corsika::cascade {
/*
DoBoundary may delete the particle (or not)
- small caveat: any changes to vParticle, or even the production
+ caveat: any changes to vParticle, or even the production
of new secondaries is currently not passed to ParticleCut,
thus, particles outside the desired phase space may be produced.
+
+ todo: this must be fixed.
*/
- fProcessSequence.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
+ process_sequence_.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
}
}
- auto decay(Particle& particle, TStackView& view) {
+ process::EProcessReturn decay(TStackView& view) {
C8LOG_DEBUG("decay");
units::si::InverseTimeType const actual_decay_time =
- fProcessSequence.GetTotalInverseLifetime(particle);
+ process_sequence_.GetInverseLifetime(view.parent());
random::UniformRealDistribution<units::si::InverseTimeType> uniDist(
actual_decay_time);
- const auto sample_process = uniDist(fRNG);
- units::si::InverseTimeType inv_decay_count = units::si::InverseTimeType::zero();
- auto const returnCode =
- fProcessSequence.SelectDecay(particle, view, sample_process, inv_decay_count);
+ const auto sample_process = uniDist(rng_);
+ auto const returnCode = process_sequence_.SelectDecay(view, sample_process);
+ if (returnCode != process::EProcessReturn::eDecayed) {
+ C8LOG_WARN("Particle did not decay!");
+ }
SetEventType(view, history::EventType::Decay);
return returnCode;
}
- auto interaction(Particle& particle, TStackView& view) {
+ process::EProcessReturn interaction(TStackView& view) {
C8LOG_DEBUG("collide");
units::si::InverseGrammageType const current_inv_length =
- fProcessSequence.GetTotalInverseInteractionLength(particle);
+ process_sequence_.GetInverseInteractionLength(view.parent());
random::UniformRealDistribution<units::si::InverseGrammageType> uniDist(
current_inv_length);
- const auto sample_process = uniDist(fRNG);
- auto inv_lambda_count = units::si::InverseGrammageType::zero();
- auto const returnCode = fProcessSequence.SelectInteraction(
- particle, view, sample_process, inv_lambda_count);
+ const auto sample_process = uniDist(rng_);
+ auto const returnCode = process_sequence_.SelectInteraction(view, sample_process);
+ if (returnCode != process::EProcessReturn::eInteracted) {
+ C8LOG_WARN("Particle did not interace!");
+ }
SetEventType(view, history::EventType::Interaction);
return returnCode;
}
+ /**
+ * set the nodes for all particles on the stack according to their numerical
+ * position
+ */
+ void setNodes() {
+ std::for_each(stack_.begin(), stack_.end(), [&](auto& p) {
+ auto const* numericalNode =
+ environment_.GetUniverse()->GetContainingNode(p.GetPosition());
+ p.SetNode(numericalNode);
+ });
+ }
+
void SetEventType(TStackView& view, [[maybe_unused]] history::EventType eventType) {
if constexpr (TStackView::has_event) {
for (auto&& sec : view) { sec.GetEvent()->setEventType(eventType); }
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index 02f974f5245f487c29bd49d00f14216b36312f7d..7cf7b9224327fe3422a86b32e3515f1a44f4af48 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -11,7 +11,7 @@
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/null_model/NullModel.h>
+#include <corsika/process/NullModel.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
@@ -132,13 +132,13 @@ TEST_CASE("Cascade", "[Cascade]") {
tracking_line::TrackingLine tracking;
stack_inspector::StackInspector<TestCascadeStack> stackInspect(1, true, E0);
- null_model::NullModel nullModel;
+ process::NullModel nullModel;
const GrammageType X0 = 20_g / square(1_cm);
const HEPEnergyType Ecrit = 85_MeV;
ProcessSplit split(X0);
ProcessCut cut(Ecrit);
- auto sequence = nullModel << stackInspect << split << cut;
+ auto sequence = process::sequence(nullModel, stackInspect, split, cut);
TestCascadeStack stack;
stack.Clear();
stack.AddParticle(
@@ -161,7 +161,6 @@ TEST_CASE("Cascade", "[Cascade]") {
}
SECTION("forced interaction") {
- EAS.SetNodes();
EAS.forceInteraction();
CHECK(stack.getEntries() == 2);
CHECK(split.GetCalls() == 1);
diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h
index d217aa50ba2299136cd8358293be1eb35964e5dd..7abcc768b5444e15414f6880c6a56aa54b288087 100644
--- a/Framework/Particles/ParticleProperties.h
+++ b/Framework/Particles/ParticleProperties.h
@@ -54,8 +54,13 @@ namespace corsika::particles {
int constexpr GetNucleusA(Code const);
int constexpr GetNucleusZ(Code const);
+} // namespace corsika::particles
+
+// here we read the implementation of all those objects and function
#include <corsika/particles/GeneratedParticleProperties.inc>
+namespace corsika::particles {
+
/*!
* returns mass of particle in natural units
*/
diff --git a/Framework/Particles/pdxml_reader.py b/Framework/Particles/pdxml_reader.py
index 8c80bec9021cc2f798443c806f7bdd794d0752d5..16978d643142950a52993baab35d69987899a79e 100755
--- a/Framework/Particles/pdxml_reader.py
+++ b/Framework/Particles/pdxml_reader.py
@@ -447,7 +447,8 @@ def gen_classes(particle_db):
#
def inc_start():
string = ('// generated by pdxml_reader.py\n'
- '// MANUAL EDITS ON OWN RISK. THEY WILL BE OVERWRITTEN. \n')
+ '// MANUAL EDITS ON OWN RISK. THEY WILL BE OVERWRITTEN. \n'
+ 'namespace corsika::particles {\n')
return string
@@ -463,14 +464,14 @@ def detail_start():
#
#
def detail_end():
- string = "\n}//end namespace detail\n"
+ string = "\n} // end namespace detail\n"
return string
###############################################################
#
#
def inc_end():
- string = ""
+ string = "\n} // end namespace corsika::particles\n"
return string
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/Framework/ProcessSequence/BaseProcess.h b/Framework/ProcessSequence/BaseProcess.h
index 42934bbe9f6309ab437f16843ccddd696e1a3e46..c1a4e232e0e0689b67aaca62af99b0abcfb47c8f 100644
--- a/Framework/ProcessSequence/BaseProcess.h
+++ b/Framework/ProcessSequence/BaseProcess.h
@@ -21,6 +21,7 @@ namespace corsika::process {
ProcessSequence. Both, the ProcessSequence and all its elements
are of type BaseProcess<T>
+ \todo rename BaseProcess into just Process
*/
class _BaseProcess {};
diff --git a/Framework/ProcessSequence/CMakeLists.txt b/Framework/ProcessSequence/CMakeLists.txt
index 811e2c3864c4d617043586842c03bcc74db3bb6f..7a75548932875f31f19f4fdd0aab3fbf19cbf867 100644
--- a/Framework/ProcessSequence/CMakeLists.txt
+++ b/Framework/ProcessSequence/CMakeLists.txt
@@ -17,7 +17,10 @@ set (
StackProcess.h
DecayProcess.h
ProcessSequence.h
+ SwitchProcessSequence.h
ProcessReturn.h
+ ProcessTraits.h
+ NullModel.h
)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAprocesssequence ${CORSIKAprocesssequence_NAMESPACE} ${CORSIKAprocesssequence_HEADERS})
@@ -27,6 +30,7 @@ target_link_libraries (
CORSIKAprocesssequence
INTERFACE
CORSIKAsetup
+ CORSIKAlogging
)
target_include_directories (
@@ -53,8 +57,26 @@ target_link_libraries (
CORSIKA_ADD_TEST (testProcessSequence)
target_link_libraries (
testProcessSequence
- ProcessSwitch
CORSIKAgeometry
CORSIKAprocesssequence
CORSIKAtesting
)
+
+# --------------------
+# code unit testing
+CORSIKA_ADD_TEST (testNullModel)
+target_link_libraries (
+ testNullModel
+ CORSIKAsetup
+ CORSIKAgeometry
+ CORSIKAunits
+ CORSIKAtesting
+ )
+
+# # CORSIKA_ADD_TEST (testSwitchProcessSequence)
+# # target_link_libraries (
+# # testSwitchProcessSequence
+# # CORSIKAgeometry
+# # CORSIKAprocesssequence
+# # CORSIKAtesting
+# # )
diff --git a/Framework/ProcessSequence/ContinuousProcess.h b/Framework/ProcessSequence/ContinuousProcess.h
index 2c9e7bbc854b6000ec2e41f93c75a8aa509b771f..c4138a755264be14c3551ddbb310dceb2495231a 100644
--- a/Framework/ProcessSequence/ContinuousProcess.h
+++ b/Framework/ProcessSequence/ContinuousProcess.h
@@ -9,7 +9,6 @@
#pragma once
#include <corsika/process/BaseProcess.h>
-#include <corsika/process/ProcessReturn.h> // for convenience
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
diff --git a/Framework/ProcessSequence/DecayProcess.h b/Framework/ProcessSequence/DecayProcess.h
index eb0cf982f08a166948becba132ab6346baa60365..09ba807434de88dbbe1952ebf1c86de7e260d0ed 100644
--- a/Framework/ProcessSequence/DecayProcess.h
+++ b/Framework/ProcessSequence/DecayProcess.h
@@ -9,8 +9,6 @@
#pragma once
#include <corsika/process/BaseProcess.h>
-#include <corsika/process/ProcessReturn.h> // for convenience
-#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
@@ -35,12 +33,17 @@ namespace corsika::process {
EProcessReturn DoDecay(TParticle&);
template <typename TParticle>
- corsika::units::si::TimeType GetLifetime(TParticle& p);
+ corsika::units::si::TimeType GetLifetime(const TParticle&);
template <typename TParticle>
- corsika::units::si::InverseTimeType GetInverseLifetime(TParticle& vP) {
- return 1. / GetRef().GetLifetime(vP);
+ corsika::units::si::InverseTimeType GetInverseLifetime(const TParticle& particle) {
+ return 1. / GetRef().GetLifetime(particle);
}
+
+ /* template <typename TParticle>
+ corsika::units::si::InverseTimeType GetInverseInteractionLength(TParticle&& particle)
+ { auto p = std::move(particle); return 1. / GetRef().GetLifetime(p);
+ }*/
};
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/InteractionProcess.h b/Framework/ProcessSequence/InteractionProcess.h
index cfecccc8156d3b4519209ba5efe9c965f1f89b93..e70729db8d2ae7ba3149b1f42ade39d78aaf29b5 100644
--- a/Framework/ProcessSequence/InteractionProcess.h
+++ b/Framework/ProcessSequence/InteractionProcess.h
@@ -9,8 +9,6 @@
#pragma once
#include <corsika/process/BaseProcess.h>
-#include <corsika/process/ProcessReturn.h> // for convenience
-#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
@@ -35,11 +33,12 @@ namespace corsika::process {
EProcessReturn DoInteraction(TParticle&);
template <typename TParticle>
- corsika::units::si::GrammageType GetInteractionLength(TParticle& p);
+ corsika::units::si::GrammageType GetInteractionLength(const TParticle&);
template <typename TParticle>
- corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& p) {
- return 1. / GetRef().GetInteractionLength(p);
+ corsika::units::si::InverseGrammageType GetInverseInteractionLength(
+ const TParticle& particle) {
+ return 1. / GetRef().GetInteractionLength(particle);
}
};
diff --git a/Framework/ProcessSequence/NullModel.h b/Framework/ProcessSequence/NullModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..feb0debd90b9063cab60c48ae96e7c16ffbc5207
--- /dev/null
+++ b/Framework/ProcessSequence/NullModel.h
@@ -0,0 +1,22 @@
+/*
+ * (c) Copyright 2018 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
+
+#include <corsika/process/BaseProcess.h>
+
+namespace corsika::process {
+
+ class NullModel : public corsika::process::BaseProcess<NullModel> {
+
+ public:
+ NullModel() = default;
+ ~NullModel() = default;
+ };
+
+} // namespace corsika::process
diff --git a/Framework/ProcessSequence/ProcessReturn.h b/Framework/ProcessSequence/ProcessReturn.h
index ae1e45b951d7a2ff144c6d683bc7b88ff93eb831..c4e7d5890411e461b8d50c4b31d44e6803b6c8d2 100644
--- a/Framework/ProcessSequence/ProcessReturn.h
+++ b/Framework/ProcessSequence/ProcessReturn.h
@@ -39,8 +39,20 @@ namespace corsika::process {
return (static_cast<int>(a) & static_cast<int>(b)) != 0;
}
+ inline bool isOk(const EProcessReturn a) {
+ return static_cast<int>(a & EProcessReturn::eOk);
+ }
+
inline bool isAbsorbed(const EProcessReturn a) {
return static_cast<int>(a & EProcessReturn::eParticleAbsorbed);
}
+ inline bool isDecayed(const EProcessReturn a) {
+ return static_cast<int>(a & EProcessReturn::eDecayed);
+ }
+
+ inline bool isInteracted(const EProcessReturn a) {
+ return static_cast<int>(a & EProcessReturn::eInteracted);
+ }
+
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h
index ea125461368ae712407070ac56b6d993dff62ea4..adffdd9ebfe443cbcc590c472b18c2ed1f4ecdcd 100644
--- a/Framework/ProcessSequence/ProcessSequence.h
+++ b/Framework/ProcessSequence/ProcessSequence.h
@@ -9,6 +9,7 @@
#pragma once
#include <corsika/process/BaseProcess.h>
+#include <corsika/process/ProcessTraits.h>
#include <corsika/process/BoundaryCrossingProcess.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/process/DecayProcess.h>
@@ -16,7 +17,9 @@
#include <corsika/process/ProcessReturn.h>
#include <corsika/process/SecondariesProcess.h>
#include <corsika/process/StackProcess.h>
+#include <corsika/process/NullModel.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/logging/Logging.h>
#include <cmath>
#include <limits>
@@ -28,101 +31,80 @@ namespace corsika::process {
A compile time static list of processes. The compiler will
generate a new type based on template logic containing all the
- elements.
+ elements provided by the user.
+
+ TProcess1 and TProcess2 must both be derived from BaseProcess,
+ and are both references if possible (lvalue), otherwise (rvalue)
+ they are just classes. This allows us to handle both, rvalue as
+ well as lvalue Processes in the ProcessSequence.
\comment Using CRTP pattern,
https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
*/
+ template <typename TProcess1, typename TProcess2 = NullModel>
+ class ProcessSequence : public BaseProcess<ProcessSequence<TProcess1, TProcess2>> {
- // this is a marker to track which BaseProcess is also a ProcessSequence
- template <typename TClass>
- struct is_process_sequence : std::false_type {};
-
- template <typename TClass>
- bool constexpr is_process_sequence_v = is_process_sequence<TClass>::value;
-
- // we also need a marker to identify SwitchProcess
- namespace switch_process {
- template <typename TProcess1, typename TProcess2>
- class SwitchProcess; // fwd-decl.
- }
-
- // to detect SwitchProcesses inside the ProcessSequence
- template <typename TClass>
- struct is_switch_process : std::false_type {};
-
- template <typename TClass>
- bool constexpr is_switch_process_v = is_switch_process<TClass>::value;
-
- template <typename Process1, typename Process2>
- struct is_switch_process<switch_process::SwitchProcess<Process1, Process2>>
- : std::true_type {};
-
- /**
- T1 and T2 are both references if possible (lvalue), otherwise
- (rvalue) they are just classes. This allows us to handle both,
- rvalue as well as lvalue Processes in the ProcessSequence.
- */
- template <typename T1, typename T2>
- class ProcessSequence : public BaseProcess<ProcessSequence<T1, T2>> {
+ using TProcess1type = typename std::decay_t<TProcess1>;
+ using TProcess2type = typename std::decay_t<TProcess2>;
- using T1type = typename std::decay<T1>::type;
- using T2type = typename std::decay<T2>::type;
+ static bool constexpr t1ProcSeq = is_process_sequence_v<TProcess1type>;
+ static bool constexpr t2ProcSeq = is_process_sequence_v<TProcess2type>;
- static bool constexpr t1ProcSeq = is_process_sequence_v<T1type>;
- static bool constexpr t2ProcSeq = is_process_sequence_v<T2type>;
+ static bool constexpr t1SwitchProcSeq = is_switch_process_sequence_v<TProcess1type>;
+ static bool constexpr t2SwitchProcSeq = is_switch_process_sequence_v<TProcess2type>;
- static bool constexpr t1SwitchProc = is_switch_process_v<T1type>;
- static bool constexpr t2SwitchProc = is_switch_process_v<T2type>;
+ TProcess1 A_; // this is a reference, if possible
+ TProcess2 B_; // this is a reference, if possible
public:
- T1 A; // this is a reference, if possible
- T2 B; // this is a reference, if possible
+ ProcessSequence(TProcess1 in_A, TProcess2 in_B)
+ : A_(in_A)
+ , B_(in_B) {}
- ProcessSequence(T1 in_A, T2 in_B)
- : A(in_A)
- , B(in_B) {}
-
- template <typename Particle, typename VTNType>
- EProcessReturn DoBoundaryCrossing(Particle& p, VTNType const& from,
- VTNType const& to) {
+ template <typename TParticle, typename TVTNType>
+ EProcessReturn DoBoundaryCrossing(TParticle& particle, TVTNType const& from,
+ TVTNType const& to) {
EProcessReturn ret = EProcessReturn::eOk;
- if constexpr (std::is_base_of_v<BoundaryCrossingProcess<T1type>, T1type> ||
+ if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess1type>,
+ TProcess1type> ||
t1ProcSeq) {
- ret |= A.DoBoundaryCrossing(p, from, to);
+ ret |= A_.DoBoundaryCrossing(particle, from, to);
}
- if constexpr (std::is_base_of_v<BoundaryCrossingProcess<T2type>, T2type> ||
+ if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess2type>,
+ TProcess2type> ||
t2ProcSeq) {
- ret |= B.DoBoundaryCrossing(p, from, to);
+ ret |= B_.DoBoundaryCrossing(particle, from, to);
}
return ret;
}
template <typename TParticle, typename TTrack>
- EProcessReturn DoContinuous(TParticle& vP, TTrack& vT) {
+ inline EProcessReturn DoContinuous(TParticle& particle, TTrack& vT) {
EProcessReturn ret = EProcessReturn::eOk;
- if constexpr (std::is_base_of_v<ContinuousProcess<T1type>, T1type> || t1ProcSeq) {
- if (!isAbsorbed(ret)) { ret |= A.DoContinuous(vP, vT); }
+ if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>, TProcess1type> ||
+ t1ProcSeq) {
+ ret |= A_.DoContinuous(particle, vT);
}
- if constexpr (std::is_base_of_v<ContinuousProcess<T2type>, T2type> || t2ProcSeq) {
- if (!isAbsorbed(ret)) { ret |= B.DoContinuous(vP, vT); }
+ if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>, TProcess2type> ||
+ t2ProcSeq) {
+ if (!isAbsorbed(ret)) { ret |= B_.DoContinuous(particle, vT); }
}
return ret;
}
template <typename TSecondaries>
- EProcessReturn DoSecondaries(TSecondaries& vS) {
- EProcessReturn ret = EProcessReturn::eOk;
- if constexpr (std::is_base_of_v<SecondariesProcess<T1type>, T1type> || t1ProcSeq) {
- ret |= A.DoSecondaries(vS);
+ inline void DoSecondaries(TSecondaries& vS) {
+ if constexpr (std::is_base_of_v<SecondariesProcess<TProcess1type>, TProcess1type> ||
+ t1ProcSeq) {
+ A_.DoSecondaries(vS);
}
- if constexpr (std::is_base_of_v<SecondariesProcess<T2type>, T2type> || t2ProcSeq) {
- ret |= B.DoSecondaries(vS);
+ if constexpr (std::is_base_of_v<SecondariesProcess<TProcess2type>, TProcess2type> ||
+ t2ProcSeq) {
+ B_.DoSecondaries(vS);
}
- return ret;
}
/**
@@ -130,16 +112,18 @@ namespace corsika::process {
so they can be exectuted only each N steps. Often these are
"maintenacne processes" that do not need to run after each
single step of the simulations. In the CheckStep function it is
- tested if either A or B are StackProcess and if they are due
+ tested if either A_ or B_ are StackProcess and if they are due
for execution.
*/
- bool CheckStep() {
+ inline bool CheckStep() {
bool ret = false;
- if constexpr (std::is_base_of_v<StackProcess<T1type>, T1type> || t1ProcSeq) {
- ret |= A.CheckStep();
+ if constexpr (std::is_base_of_v<StackProcess<TProcess1type>, TProcess1type> ||
+ (t1ProcSeq && !t1SwitchProcSeq)) {
+ ret |= A_.CheckStep();
}
- if constexpr (std::is_base_of_v<StackProcess<T2type>, T2type> || t2ProcSeq) {
- ret |= B.CheckStep();
+ if constexpr (std::is_base_of_v<StackProcess<TProcess2type>, TProcess2type> ||
+ (t2ProcSeq && !t2SwitchProcSeq)) {
+ ret |= B_.CheckStep();
}
return ret;
}
@@ -148,189 +132,257 @@ namespace corsika::process {
Execute the StackProcess-es in the ProcessSequence
*/
template <typename TStack>
- EProcessReturn DoStack(TStack& vS) {
- EProcessReturn ret = EProcessReturn::eOk;
- if constexpr (std::is_base_of_v<StackProcess<T1type>, T1type> || t1ProcSeq) {
- if (A.CheckStep()) { ret |= A.DoStack(vS); }
+ inline void DoStack(TStack& stack) {
+ if constexpr (std::is_base_of_v<StackProcess<TProcess1type>, TProcess1type> ||
+ (t1ProcSeq && !t1SwitchProcSeq)) {
+ if (A_.CheckStep()) { A_.DoStack(stack); }
}
- if constexpr (std::is_base_of_v<StackProcess<T2type>, T2type> || t2ProcSeq) {
- if (B.CheckStep()) { ret |= B.DoStack(vS); }
+ if constexpr (std::is_base_of_v<StackProcess<TProcess2type>, TProcess2type> ||
+ (t2ProcSeq && !t2SwitchProcSeq)) {
+ if (B_.CheckStep()) { B_.DoStack(stack); }
}
- return ret;
}
template <typename TParticle, typename TTrack>
- corsika::units::si::LengthType MaxStepLength(TParticle& vP, TTrack& vTrack) {
+ inline corsika::units::si::LengthType MaxStepLength(TParticle& particle,
+ TTrack& vTrack) {
corsika::units::si::LengthType
max_length = // if no other process in the sequence implements it
std::numeric_limits<double>::infinity() * corsika::units::si::meter;
- if constexpr (std::is_base_of_v<ContinuousProcess<T1type>, T1type> || t1ProcSeq) {
- corsika::units::si::LengthType const len = A.MaxStepLength(vP, vTrack);
+ if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>, TProcess1type> ||
+ t1ProcSeq) {
+ corsika::units::si::LengthType const len = A_.MaxStepLength(particle, vTrack);
max_length = std::min(max_length, len);
}
- if constexpr (std::is_base_of_v<ContinuousProcess<T2type>, T2type> || t2ProcSeq) {
- corsika::units::si::LengthType const len = B.MaxStepLength(vP, vTrack);
+ if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>, TProcess2type> ||
+ t2ProcSeq) {
+ corsika::units::si::LengthType const len = B_.MaxStepLength(particle, vTrack);
max_length = std::min(max_length, len);
}
return max_length;
}
template <typename TParticle>
- corsika::units::si::GrammageType GetTotalInteractionLength(TParticle& vP) {
- return 1. / GetInverseInteractionLength(vP);
- }
-
- template <typename TParticle>
- corsika::units::si::InverseGrammageType GetTotalInverseInteractionLength(
- TParticle& vP) {
- return GetInverseInteractionLength(vP);
+ inline corsika::units::si::GrammageType GetInteractionLength(TParticle&& particle) {
+ return 1. / GetInverseInteractionLength(particle);
}
template <typename TParticle>
- corsika::units::si::InverseGrammageType GetInverseInteractionLength(TParticle& vP) {
+ inline corsika::units::si::InverseGrammageType GetInverseInteractionLength(
+ TParticle&& particle) {
using namespace corsika::units::si;
- InverseGrammageType tot = 0 * meter * meter / gram;
+ InverseGrammageType tot = 0 * meter * meter / gram; // default value
- if constexpr (std::is_base_of_v<InteractionProcess<T1type>, T1type> || t1ProcSeq ||
- t1SwitchProc) {
- tot += A.GetInverseInteractionLength(vP);
+ if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>, TProcess1type> ||
+ t1ProcSeq) {
+ tot += A_.GetInverseInteractionLength(particle);
}
- if constexpr (std::is_base_of_v<InteractionProcess<T2type>, T2type> || t2ProcSeq ||
- t2SwitchProc) {
- tot += B.GetInverseInteractionLength(vP);
+ if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>, TProcess2type> ||
+ t2ProcSeq) {
+ tot += B_.GetInverseInteractionLength(particle);
}
return tot;
}
- template <typename TParticle, typename TSecondaries>
- EProcessReturn SelectInteraction(
- TParticle& vP, TSecondaries& vS,
- [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_select,
- corsika::units::si::InverseGrammageType& lambda_inv_count) {
+ template <typename TSecondaryView>
+ inline EProcessReturn SelectInteraction(
+ TSecondaryView& view,
+ [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_select,
+ [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_sum =
+ corsika::units::si::InverseGrammageType::zero()) {
+
+ // TODO: add check for lambda_inv_select>lambda_inv_tot
- if constexpr (t1ProcSeq || t1SwitchProc) {
+ if constexpr (t1ProcSeq) {
// if A is a process sequence --> check inside
const EProcessReturn ret =
- A.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
- // if A did succeed, stop routine
+ A_.SelectInteraction(view, lambda_inv_select, lambda_inv_sum);
+ // if A_ did succeed, stop routine. Not checking other static branch B_.
if (ret != EProcessReturn::eOk) { return ret; }
- } else if constexpr (std::is_base_of_v<InteractionProcess<T1type>, T1type>) {
+ } else if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>,
+ TProcess1type>) {
// if this is not a ContinuousProcess --> evaluate probability
- lambda_inv_count += A.GetInverseInteractionLength(vP);
+ const auto particle = view.parent();
+ lambda_inv_sum += A_.GetInverseInteractionLength(particle);
// check if we should execute THIS process and then EXIT
- if (lambda_select < lambda_inv_count) {
- A.DoInteraction(vS);
+ if (lambda_inv_select < lambda_inv_sum) {
+ A_.DoInteraction(view);
return EProcessReturn::eInteracted;
}
- } // end branch A
+ } // end branch A_
- if constexpr (t2ProcSeq || t2SwitchProc) {
- // if A is a process sequence --> check inside
- const EProcessReturn ret =
- B.SelectInteraction(vP, vS, lambda_select, lambda_inv_count);
- // if A did succeed, stop routine
- if (ret != EProcessReturn::eOk) { return ret; }
- } else if constexpr (std::is_base_of_v<InteractionProcess<T2type>, T2type>) {
+ if constexpr (t2ProcSeq) {
+ // if B_ is a process sequence --> check inside
+ return B_.SelectInteraction(view, lambda_inv_select, lambda_inv_sum);
+ } else if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>,
+ TProcess2type>) {
// if this is not a ContinuousProcess --> evaluate probability
- lambda_inv_count += B.GetInverseInteractionLength(vP);
+ lambda_inv_sum += B_.GetInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
- if (lambda_select < lambda_inv_count) {
- B.DoInteraction(vS);
+ if (lambda_inv_select < lambda_inv_sum) {
+ B_.DoInteraction(view);
return EProcessReturn::eInteracted;
}
- } // end branch A
+ } // end branch B_
return EProcessReturn::eOk;
}
template <typename TParticle>
- corsika::units::si::TimeType GetTotalLifetime(TParticle& p) {
- return 1. / GetInverseLifetime(p);
- }
-
- template <typename TParticle>
- corsika::units::si::InverseTimeType GetTotalInverseLifetime(TParticle& p) {
- return GetInverseLifetime(p);
+ inline corsika::units::si::TimeType GetLifetime(TParticle& particle) {
+ return 1. / GetInverseLifetime(particle);
}
template <typename TParticle>
- corsika::units::si::InverseTimeType GetInverseLifetime(TParticle& p) {
+ inline corsika::units::si::InverseTimeType GetInverseLifetime(TParticle&& particle) {
using namespace corsika::units::si;
- corsika::units::si::InverseTimeType tot = 0 / second;
+ corsika::units::si::InverseTimeType tot = 0 / second; // default value
- if constexpr (std::is_base_of_v<DecayProcess<T1type>, T1type> || t1ProcSeq) {
- tot += A.GetInverseLifetime(p);
+ if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>, TProcess1type> ||
+ t1ProcSeq) {
+ tot += A_.GetInverseLifetime(particle);
}
- if constexpr (std::is_base_of_v<DecayProcess<T2type>, T2type> || t2ProcSeq) {
- tot += B.GetInverseLifetime(p);
+ if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>, TProcess2type> ||
+ t2ProcSeq) {
+ tot += B_.GetInverseLifetime(particle);
}
return tot;
}
// select decay process
- template <typename TParticle, typename TSecondaries>
- EProcessReturn SelectDecay(
- TParticle& vP, TSecondaries& vS,
- [[maybe_unused]] corsika::units::si::InverseTimeType decay_select,
- corsika::units::si::InverseTimeType& decay_inv_count) {
+ template <typename TSecondaryView>
+ inline EProcessReturn SelectDecay(
+ TSecondaryView& view,
+ [[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_select,
+ [[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_sum =
+ corsika::units::si::InverseTimeType::zero()) {
+
+ // TODO: add check for decay_inv_select>decay_inv_tot
+
if constexpr (t1ProcSeq) {
- // if A is a process sequence --> check inside
- const EProcessReturn ret = A.SelectDecay(vP, vS, decay_select, decay_inv_count);
- // if A did succeed, stop routine
+ // if A_ is a process sequence --> check inside
+ const EProcessReturn ret = A_.SelectDecay(view, decay_inv_select, decay_inv_sum);
+ // if A_ did succeed, stop routine here (not checking other static branch B_)
if (ret != EProcessReturn::eOk) { return ret; }
- } else if constexpr (std::is_base_of_v<DecayProcess<T1type>, T1type>) {
+ } else if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>,
+ TProcess1type>) {
// if this is not a ContinuousProcess --> evaluate probability
- decay_inv_count += A.GetInverseLifetime(vP);
+ decay_inv_sum += A_.GetInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
- if (decay_select < decay_inv_count) { // more pedagogical: rndm_select <
- // decay_inv_count / decay_inv_tot
- A.DoDecay(vS);
+ if (decay_inv_select < decay_inv_sum) { // more pedagogical: rndm_select <
+ // decay_inv_sum / decay_inv_tot
+ A_.DoDecay(view);
return EProcessReturn::eDecayed;
}
- } // end branch A
+ } // end branch A_
if constexpr (t2ProcSeq) {
- // if A is a process sequence --> check inside
- const EProcessReturn ret = B.SelectDecay(vP, vS, decay_select, decay_inv_count);
- // if A did succeed, stop routine
- if (ret != EProcessReturn::eOk) { return ret; }
- } else if constexpr (std::is_base_of_v<DecayProcess<T2type>, T2type>) {
+ // if B_ is a process sequence --> check inside
+ return B_.SelectDecay(view, decay_inv_select, decay_inv_sum);
+ } else if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>,
+ TProcess2type>) {
// if this is not a ContinuousProcess --> evaluate probability
- decay_inv_count += B.GetInverseLifetime(vP);
+ decay_inv_sum += B_.GetInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
- if (decay_select < decay_inv_count) {
- B.DoDecay(vS);
+ if (decay_inv_select < decay_inv_sum) {
+ B_.DoDecay(view);
return EProcessReturn::eDecayed;
}
- } // end branch B
+ } // end branch B_
return EProcessReturn::eOk;
}
};
- // the << operator assembles many BaseProcess, ContinuousProcess, and
- // Interaction/DecayProcess objects into a ProcessSequence, all combinatorics
- // must be allowed, this is why we define a macro to define all
- // combinations here:
-
- // enable the << operator to construct ProcessSequence from two
- // Processes, only if poth Processes derive from BaseProcesses
+ /**
+ * \function sequence
+ *
+ * to construct ProcessSequences in a flexible and dynamic way the
+ * `sequence` factory functions are provided
+ *
+ * Any objects of type
+ * - BaseProcess,
+ * - ContinuousProcess, and
+ * - Interaction/DecayProcess,
+ * - StackProcess,
+ * - SecondariesProcess
+ * can be assembled into a ProcessSequence, all
+ * combinatorics are allowed.
+
+ * The sequence function checks that all its arguments are all of
+ * types derived from BaseProcess. Also the ProcessSequence itself
+ * is derived from type BaseProcess
+ **/
+
+ template <typename... TProcesses, typename TProcess1>
+ inline typename std::enable_if_t<
+ std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess1>>,
+ typename std::decay_t<TProcess1>>,
+ ProcessSequence<TProcess1, decltype(sequence(std::declval<TProcesses>()...))>>
+ sequence(TProcess1&& vA, TProcesses&&... vBs) {
+ return ProcessSequence<TProcess1, decltype(sequence(std::declval<TProcesses>()...))>(
+ vA, sequence(std::forward<TProcesses>(vBs)...));
+ }
template <typename TProcess1, typename TProcess2>
- inline typename std::enable_if<
- std::is_base_of<BaseProcess<typename std::decay<TProcess1>::type>,
- typename std::decay<TProcess1>::type>::value &&
- std::is_base_of<BaseProcess<typename std::decay<TProcess2>::type>,
- typename std::decay<TProcess2>::type>::value,
- ProcessSequence<TProcess1, TProcess2>>::type
- operator<<(TProcess1&& vA, TProcess2&& vB) {
+ inline typename std::enable_if_t<
+ std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess1>>,
+ typename std::decay_t<TProcess1>> &&
+ std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess2>>,
+ typename std::decay_t<TProcess2>>,
+ ProcessSequence<TProcess1, TProcess2>>
+ sequence(TProcess1&& vA, TProcess2&& vB) {
return ProcessSequence<TProcess1, TProcess2>(vA, vB);
}
- /// marker to identify objectas ProcessSequence
+ /**
+ * also allow a single Process in ProcessSequence, accompany by
+ * `NullModel`
+ **/
+ template <typename TProcess>
+ inline typename std::enable_if_t<
+ std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess>>,
+ typename std::decay_t<TProcess>>,
+ ProcessSequence<TProcess, NullModel>>
+ sequence(TProcess&& vA) {
+ return ProcessSequence<TProcess, NullModel>(vA, NullModel());
+ }
+
+ /**
+ * traits marker to identify objectas ProcessSequence
+ **/
+ template <typename TProcess1, typename TProcess2>
+ struct is_process_sequence<ProcessSequence<TProcess1, TProcess2>> : std::true_type {
+ // only switch on for BaseProcesses
+ template <typename std::enable_if_t<
+ std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess1>>,
+ typename std::decay_t<TProcess1>> &&
+ std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess2>>,
+ typename std::decay_t<TProcess2>>,
+ int>>
+ is_process_sequence() {}
+ };
+
+ /**
+ * traits marker to identify objects containing any StackProcesses
+ **/
+ namespace detail {
+ // need helper alias to achieve this:
+ template <typename TProcess1, typename TProcess2,
+ typename = typename std::enable_if_t<
+ contains_stack_process_v<TProcess1> ||
+ std::is_base_of_v<StackProcess<typename std::decay_t<TProcess1>>,
+ typename std::decay_t<TProcess1>> ||
+ contains_stack_process_v<TProcess2> ||
+ std::is_base_of_v<StackProcess<typename std::decay_t<TProcess2>>,
+ typename std::decay_t<TProcess2>>,
+ int>>
+ using enable_if_stack = ProcessSequence<TProcess1, TProcess2>;
+ } // namespace detail
+
template <typename TProcess1, typename TProcess2>
- struct is_process_sequence<corsika::process::ProcessSequence<TProcess1, TProcess2>>
+ struct contains_stack_process<detail::enable_if_stack<TProcess1, TProcess2>>
: std::true_type {};
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/ProcessTraits.h b/Framework/ProcessSequence/ProcessTraits.h
new file mode 100644
index 0000000000000000000000000000000000000000..e0cb00ebfde8827bbb098dc34d13df1c94a72eb2
--- /dev/null
+++ b/Framework/ProcessSequence/ProcessTraits.h
@@ -0,0 +1,41 @@
+/*
+ * (c) Copyright 2018 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::process {
+
+ /**
+ * A traits marker to track which BaseProcess is also a ProcessSequence
+ **/
+ template <typename TClass>
+ struct is_process_sequence : std::false_type {};
+
+ template <typename TClass>
+ bool constexpr is_process_sequence_v = is_process_sequence<TClass>::value;
+
+ /**
+ * A traits marker to identiy a BaseProcess that is also SwitchProcessesSequence
+ **/
+
+ template <typename TClass>
+ struct is_switch_process_sequence : std::false_type {};
+
+ template <typename TClass>
+ bool constexpr is_switch_process_sequence_v = is_switch_process_sequence<TClass>::value;
+
+ /**
+ * A traits marker to identify ProcessSequence that contain a StackProcess
+ **/
+ template <typename TClass>
+ struct contains_stack_process : std::false_type {};
+
+ template <typename TClass>
+ bool constexpr contains_stack_process_v = contains_stack_process<TClass>::value;
+
+} // namespace corsika::process
diff --git a/Framework/ProcessSequence/SecondariesProcess.h b/Framework/ProcessSequence/SecondariesProcess.h
index 355a465c78bbdf4008b63f18e091a58321b4a949..a9bc3832a2e8a82949a6d03c3b9c1c5d44a49828 100644
--- a/Framework/ProcessSequence/SecondariesProcess.h
+++ b/Framework/ProcessSequence/SecondariesProcess.h
@@ -9,8 +9,6 @@
#pragma once
#include <corsika/process/BaseProcess.h>
-#include <corsika/process/ProcessReturn.h> // for convenience
-#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
@@ -30,7 +28,7 @@ namespace corsika::process {
/// here starts the interface-definition part
// -> enforce TDerived to implement DoSecondaries...
template <typename TSecondaries>
- inline EProcessReturn DoSecondaries(TSecondaries&);
+ inline void DoSecondaries(TSecondaries&);
};
} // namespace corsika::process
diff --git a/Framework/ProcessSequence/StackProcess.h b/Framework/ProcessSequence/StackProcess.h
index 64a5dff97b6d2df4df5837b3bcc0e2b3052938e6..ce5c0998df246ed5ceb20bed618b40814a79077f 100644
--- a/Framework/ProcessSequence/StackProcess.h
+++ b/Framework/ProcessSequence/StackProcess.h
@@ -9,8 +9,6 @@
#pragma once
#include <corsika/process/BaseProcess.h>
-#include <corsika/process/ProcessReturn.h> // for convenience
-#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
@@ -36,7 +34,7 @@ namespace corsika::process {
/// here starts the interface-definition part
// -> enforce TDerived to implement DoStack...
template <typename TStack>
- inline EProcessReturn DoStack(TStack&);
+ inline void DoStack(TStack&);
int GetStep() const { return fIStep; }
bool CheckStep() { return !((++fIStep) % fNStep); }
diff --git a/Framework/ProcessSequence/SwitchProcessSequence.h b/Framework/ProcessSequence/SwitchProcessSequence.h
new file mode 100644
index 0000000000000000000000000000000000000000..f1b1dfcf5eee747f2007d207788ad2f686285766
--- /dev/null
+++ b/Framework/ProcessSequence/SwitchProcessSequence.h
@@ -0,0 +1,393 @@
+/*
+ * (c) Copyright 2018 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
+
+#include <corsika/process/BaseProcess.h>
+#include <corsika/process/ProcessTraits.h>
+#include <corsika/process/BoundaryCrossingProcess.h>
+#include <corsika/process/ContinuousProcess.h>
+#include <corsika/process/DecayProcess.h>
+#include <corsika/process/InteractionProcess.h>
+#include <corsika/process/ProcessReturn.h>
+#include <corsika/process/SecondariesProcess.h>
+#include <corsika/process/StackProcess.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <cmath>
+#include <limits>
+#include <type_traits>
+
+namespace corsika::process {
+
+ // enum for the process switch selection: identify if First or
+ // Second process branch should be used.
+ enum class SwitchResult { First, Second };
+
+ /**
+ \class SwitchProcessSequence
+
+ A compile time static list of processes that uses an internal
+ TSelect class to switch between different versions of processes
+ (or process sequence).
+
+ TProcess1 and TProcess2 must be derived from BaseProcess and are
+ both references if possible (lvalue), otherwise (rvalue) they are
+ just classes. This allows us to handle both, rvalue as well as
+ lvalue Processes in the SwitchProcessSequence.
+
+ TSelect has to implement a `operator()(const Particle&)` and has to
+ return either SwitchResult::First or SwitchResult::Second. Note:
+ TSelect may absolutely also use random numbers to sample between
+ its results. This can be used to achieve arbitrarily smooth
+ transition or mixtures of processes.
+
+ Warning: do not put StackProcess into a SwitchProcessSequence
+ since this makes no sense. The StackProcess acts on an entire
+ particle stack and not on indiviidual particles.
+
+ \comment See also class ProcessSequence
+ **/
+
+ template <typename TProcess1, typename TProcess2, typename TSelect>
+ class SwitchProcessSequence
+ : public BaseProcess<SwitchProcessSequence<TProcess1, TProcess2, TSelect>> {
+
+ using TProcess1type = typename std::decay_t<TProcess1>;
+ using TProcess2type = typename std::decay_t<TProcess2>;
+
+ static bool constexpr t1ProcSeq = is_process_sequence_v<TProcess1type>;
+ static bool constexpr t2ProcSeq = is_process_sequence_v<TProcess2type>;
+
+ // make sure only BaseProcess types TProcess1/2 are passed
+ static_assert(std::is_base_of_v<BaseProcess<TProcess1type>, TProcess1type>,
+ "can only use process derived from BaseProcess in "
+ "SwitchProcessSequence, for Process 1");
+ static_assert(std::is_base_of_v<BaseProcess<TProcess2type>, TProcess2type>,
+ "can only use process derived from BaseProcess in "
+ "SwitchProcessSequence, for Process 2");
+
+ // make sure TSelect is a function
+ static_assert(!std::is_function_v<TSelect>, "TSelect must be a function type");
+
+ // make sure none of TProcess1/2 is a StackProcess
+ static_assert(!std::is_base_of_v<StackProcess<TProcess1type>, TProcess1type>,
+ "cannot use StackProcess in SwitchProcessSequence, for Process 1");
+ static_assert(!std::is_base_of_v<StackProcess<TProcess2type>, TProcess2type>,
+ "cannot use StackProcess in SwitchProcessSequence, for Process 2");
+
+ // if TProcess1/2 are already ProcessSequences, make sure they do not contain
+ // StackProcess
+ // if constexpr (t1ProcSeq) {
+ static_assert(!contains_stack_process_v<TProcess1type>,
+ "cannot use StackProcess in SwitchProcessSequence, remove from "
+ "ProcessSequence 1");
+ //}
+
+ // if constexpr (t2ProcSeq)
+ static_assert(!contains_stack_process_v<TProcess2type>,
+ "cannot use StackProcess in SwitchProcessSequence, remove from "
+ "ProcessSequence 2");
+
+ TSelect select_; // this is a reference, if possible
+
+ TProcess1 A_; // this is a reference, if possible
+ TProcess2 B_; // this is a reference, if possible
+
+ public:
+ SwitchProcessSequence(TProcess1 in_A, TProcess2 in_B, TSelect sel)
+ : select_(sel)
+ , A_(in_A)
+ , B_(in_B) {}
+
+ template <typename TParticle, typename TVTNType>
+ EProcessReturn DoBoundaryCrossing(TParticle& particle, TVTNType const& from,
+ TVTNType const& to) {
+ switch (select_(particle)) {
+ case SwitchResult::First: {
+ if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess1type>,
+ TProcess1type> ||
+ t1ProcSeq) {
+ return A_.DoBoundaryCrossing(particle, from, to);
+ }
+ break;
+ }
+ case SwitchResult::Second: {
+ if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess2type>,
+ TProcess2type> ||
+ t2ProcSeq) {
+ return B_.DoBoundaryCrossing(particle, from, to);
+ }
+ break;
+ }
+ }
+ return EProcessReturn::eOk;
+ }
+
+ template <typename TParticle, typename TTrack>
+ inline EProcessReturn DoContinuous(TParticle& particle, TTrack& vT) {
+ switch (select_(particle)) {
+ case SwitchResult::First: {
+ if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>,
+ TProcess1type> ||
+ t1ProcSeq) {
+ return A_.DoContinuous(particle, vT);
+ }
+ break;
+ }
+ case SwitchResult::Second: {
+ if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>,
+ TProcess2type> ||
+ t2ProcSeq) {
+ return B_.DoContinuous(particle, vT);
+ }
+ break;
+ }
+ }
+ return EProcessReturn::eOk;
+ }
+
+ template <typename TSecondaries>
+ inline void DoSecondaries(TSecondaries& vS) {
+ const auto& particle = vS.parent();
+ switch (select_(particle)) {
+ case SwitchResult::First: {
+ if constexpr (std::is_base_of_v<SecondariesProcess<TProcess1type>,
+ TProcess1type> ||
+ t1ProcSeq) {
+ A_.DoSecondaries(vS);
+ }
+ break;
+ }
+ case SwitchResult::Second: {
+ if constexpr (std::is_base_of_v<SecondariesProcess<TProcess2type>,
+ TProcess2type> ||
+ t2ProcSeq) {
+ B_.DoSecondaries(vS);
+ }
+ break;
+ }
+ }
+ }
+
+ template <typename TParticle, typename TTrack>
+ inline corsika::units::si::LengthType MaxStepLength(TParticle& particle,
+ TTrack& vTrack) {
+ switch (select_(particle)) {
+ case SwitchResult::First: {
+ if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>,
+ TProcess1type> ||
+ t1ProcSeq) {
+ return A_.MaxStepLength(particle, vTrack);
+ }
+ break;
+ }
+ case SwitchResult::Second: {
+ if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>,
+ TProcess2type> ||
+ t2ProcSeq) {
+ return B_.MaxStepLength(particle, vTrack);
+ }
+ break;
+ }
+ }
+
+ // if no other process in the sequence implements it
+ return std::numeric_limits<double>::infinity() * corsika::units::si::meter;
+ }
+
+ template <typename TParticle>
+ inline corsika::units::si::GrammageType GetInteractionLength(TParticle&& particle) {
+ return 1. / GetInverseInteractionLength(particle);
+ }
+
+ template <typename TParticle>
+ inline corsika::units::si::InverseGrammageType GetInverseInteractionLength(
+ TParticle&& particle) {
+ using namespace corsika::units::si;
+
+ switch (select_(particle)) {
+ case SwitchResult::First: {
+ if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>,
+ TProcess1type> ||
+ t1ProcSeq) {
+ return A_.GetInverseInteractionLength(particle);
+ }
+ break;
+ }
+ case SwitchResult::Second: {
+ if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>,
+ TProcess2type> ||
+ t2ProcSeq) {
+ return B_.GetInverseInteractionLength(particle);
+ }
+ break;
+ }
+ }
+ return 0 * meter * meter / gram; // default value
+ }
+
+ template <typename TSecondaryView>
+ inline EProcessReturn SelectInteraction(
+ TSecondaryView& view,
+ [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_select,
+ [[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_sum =
+ corsika::units::si::InverseGrammageType::zero()) {
+ switch (select_(view.parent())) {
+ case SwitchResult::First: {
+ if constexpr (t1ProcSeq) {
+ // if A_ is a process sequence --> check inside
+ const EProcessReturn ret =
+ A_.SelectInteraction(view, lambda_inv_select, lambda_inv_sum);
+ // if A_ did succeed, stop routine. Not checking other static branch B_.
+ if (ret != EProcessReturn::eOk) { return ret; }
+ } else if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>,
+ TProcess1type>) {
+ // if this is not a ContinuousProcess --> evaluate probability
+ lambda_inv_sum += A_.GetInverseInteractionLength(view.parent());
+ // check if we should execute THIS process and then EXIT
+ if (lambda_inv_select < lambda_inv_sum) {
+ A_.DoInteraction(view);
+ return EProcessReturn::eInteracted;
+ }
+ } // end branch A_
+ break;
+ }
+
+ case SwitchResult::Second: {
+
+ if constexpr (t2ProcSeq) {
+ // if B_ is a process sequence --> check inside
+ return B_.SelectInteraction(view, lambda_inv_select, lambda_inv_sum);
+ } else if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>,
+ TProcess2type>) {
+ // if this is not a ContinuousProcess --> evaluate probability
+ lambda_inv_sum += B_.GetInverseInteractionLength(view.parent());
+ // check if we should execute THIS process and then EXIT
+ if (lambda_inv_select < lambda_inv_sum) {
+ B_.DoInteraction(view);
+ return EProcessReturn::eInteracted;
+ }
+ } // end branch B_
+ break;
+ }
+ }
+ return EProcessReturn::eOk;
+ }
+
+ template <typename TParticle>
+ inline corsika::units::si::TimeType GetLifetime(TParticle&& particle) {
+ return 1. / GetInverseLifetime(particle);
+ }
+
+ template <typename TParticle>
+ inline corsika::units::si::InverseTimeType GetInverseLifetime(TParticle&& particle) {
+ using namespace corsika::units::si;
+
+ switch (select_(particle)) {
+ case SwitchResult::First: {
+ if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>, TProcess1type> ||
+ t1ProcSeq) {
+ return A_.GetInverseLifetime(particle);
+ }
+ break;
+ }
+
+ case SwitchResult::Second: {
+ if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>, TProcess2type> ||
+ t2ProcSeq) {
+ return B_.GetInverseLifetime(particle);
+ }
+ break;
+ }
+ }
+ return 0 / second; // default value
+ }
+
+ // select decay process
+ template <typename TSecondaryView>
+ inline EProcessReturn SelectDecay(
+ TSecondaryView& view,
+ [[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_select,
+ [[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_sum =
+ corsika::units::si::InverseTimeType::zero()) {
+ switch (select_(view.parent())) {
+ case SwitchResult::First: {
+ if constexpr (t1ProcSeq) {
+ // if A_ is a process sequence --> check inside
+ const EProcessReturn ret =
+ A_.SelectDecay(view, decay_inv_select, decay_inv_sum);
+ // if A_ did succeed, stop routine here (not checking other static branch B_)
+ if (ret != EProcessReturn::eOk) { return ret; }
+ } else if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>,
+ TProcess1type>) {
+ // if this is not a ContinuousProcess --> evaluate probability
+ decay_inv_sum += A_.GetInverseLifetime(view.parent());
+ // check if we should execute THIS process and then EXIT
+ if (decay_inv_select < decay_inv_sum) {
+ // more pedagogical: rndm_select < decay_inv_sum / decay_inv_tot
+ A_.DoDecay(view);
+ return EProcessReturn::eDecayed;
+ }
+ } // end branch A_
+ break;
+ }
+
+ case SwitchResult::Second: {
+
+ if constexpr (t2ProcSeq) {
+ // if B_ is a process sequence --> check inside
+ return B_.SelectDecay(view, decay_inv_select, decay_inv_sum);
+ } else if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>,
+ TProcess2type>) {
+ // if this is not a ContinuousProcess --> evaluate probability
+ decay_inv_sum += B_.GetInverseLifetime(view.parent());
+ // check if we should execute THIS process and then EXIT
+ if (decay_inv_select < decay_inv_sum) {
+ B_.DoDecay(view);
+ return EProcessReturn::eDecayed;
+ }
+ } // end branch B_
+ break;
+ }
+ }
+ return EProcessReturn::eOk;
+ }
+ };
+
+ // the method `select(proc1,proc1,selector)` assembles many
+ // BaseProcesses, and ProcessSequences into a SwitchProcessSequence,
+ // all combinatorics must be allowed, this is why we define a macro
+ // to define all combinations here:
+
+ // Both, Processes1 and Processes2, must derive from BaseProcesses
+
+ template <typename TProcess1, typename TProcess2, typename TSelect>
+ inline typename std::enable_if_t<
+ std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess1>>,
+ typename std::decay_t<TProcess1>> &&
+ std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess2>>,
+ typename std::decay_t<TProcess2>>,
+ SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
+ select(TProcess1&& vA, TProcess2&& vB, TSelect selector) {
+ return SwitchProcessSequence<TProcess1, TProcess2, TSelect>(vA, vB, selector);
+ }
+
+ /// traits marker to identify objectas ProcessSequence
+ template <typename TProcess1, typename TProcess2, typename TSelect>
+ struct is_process_sequence<
+ corsika::process::SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
+ : std::true_type {};
+
+ /// traits marker to identify objectas SwitchProcessSequence
+ template <typename TProcess1, typename TProcess2, typename TSelect>
+ struct is_switch_process_sequence<
+ corsika::process::SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
+ : std::true_type {};
+
+} // namespace corsika::process
diff --git a/Framework/ProcessSequence/testNullModel.cc b/Framework/ProcessSequence/testNullModel.cc
new file mode 100644
index 0000000000000000000000000000000000000000..33bfda3d01cda939f6205271cb2d00d3e6d8cb7e
--- /dev/null
+++ b/Framework/ProcessSequence/testNullModel.cc
@@ -0,0 +1,29 @@
+/*
+ * (c) Copyright 2018 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.
+ */
+
+#include <catch2/catch.hpp>
+
+#include <corsika/process/NullModel.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+using namespace corsika::units::si;
+using namespace corsika::process;
+using namespace corsika;
+
+TEST_CASE("NullModel", "[processes]") {
+
+ SECTION("interface") { [[maybe_unused]] NullModel model; }
+}
diff --git a/Framework/ProcessSequence/testProcessSequence.cc b/Framework/ProcessSequence/testProcessSequence.cc
index b51876886acbd22a80dd20aca4abe5d50c23e216..875cbc5821fe5e67da0ece3977f13446f2a0659a 100644
--- a/Framework/ProcessSequence/testProcessSequence.cc
+++ b/Framework/ProcessSequence/testProcessSequence.cc
@@ -11,9 +11,11 @@
#include <array>
#include <iomanip>
#include <iostream>
+#include <typeinfo>
#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/switch_process/SwitchProcess.h>
+#include <corsika/process/SwitchProcessSequence.h>
+#include <corsika/units/PhysicalUnits.h>
using namespace corsika;
using namespace corsika::units::si;
@@ -22,7 +24,12 @@ using namespace std;
static const int nData = 10;
-int globalCount = 0;
+int globalCount = 0; // simple counter
+
+int checkDecay = 0; // use this as a bit field
+int checkInteract = 0; // use this as a bit field
+int checkSec = 0; // use this as a bit field
+int checkCont = 0; // use this as a bit field
class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> {
int fV = 0;
@@ -38,7 +45,8 @@ public:
template <typename D, typename T>
inline EProcessReturn DoContinuous(D& d, T&) const {
cout << "ContinuousProcess1::DoContinuous" << endl;
- for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
+ checkCont |= 1;
+ for (int i = 0; i < nData; ++i) d.data_[i] += 0.933;
return EProcessReturn::eOk;
}
};
@@ -56,7 +64,27 @@ public:
template <typename D, typename T>
inline EProcessReturn DoContinuous(D& d, T&) const {
cout << "ContinuousProcess2::DoContinuous" << endl;
- for (int i = 0; i < nData; ++i) d.p[i] += 0.933;
+ checkCont |= 2;
+ for (int i = 0; i < nData; ++i) d.data_[i] += 0.111;
+ return EProcessReturn::eOk;
+ }
+};
+
+class ContinuousProcess3 : public ContinuousProcess<ContinuousProcess3> {
+ int fV = 0;
+
+public:
+ ContinuousProcess3(const int v)
+ : fV(v) {
+ cout << "globalCount: " << globalCount << ", fV: " << fV << std::endl;
+ globalCount++;
+ }
+
+ template <typename D, typename T>
+ inline EProcessReturn DoContinuous(D& d, T&) const {
+ cout << "ContinuousProcess3::DoContinuous" << endl;
+ checkCont |= 4;
+ for (int i = 0; i < nData; ++i) d.data_[i] += 0.333;
return EProcessReturn::eOk;
}
};
@@ -69,12 +97,18 @@ public:
globalCount++;
}
- template <typename D, typename S>
- inline EProcessReturn DoInteraction(D& d, S&) const {
- for (int i = 0; i < nData; ++i) d.p[i] += 1 + i;
+ template <typename TView>
+ inline EProcessReturn DoInteraction(TView& v) const {
+ checkInteract |= 1;
+ for (int i = 0; i < nData; ++i) v.parent().data_[i] += 1 + i;
return EProcessReturn::eOk;
}
+ template <typename TParticle>
+ corsika::units::si::GrammageType GetInteractionLength(TParticle&) const {
+ return 10_g / square(1_cm);
+ }
+
private:
int fV;
};
@@ -89,15 +123,17 @@ public:
globalCount++;
}
- template <typename Particle>
- inline EProcessReturn DoInteraction(Particle&) const {
+ template <typename TView>
+ inline EProcessReturn DoInteraction(TView& v) const {
+ checkInteract |= 2;
+ for (int i = 0; i < nData; ++i) v.parent().data_[i] /= 1.1;
cout << "Process2::DoInteraction" << endl;
return EProcessReturn::eOk;
}
template <typename Particle>
GrammageType GetInteractionLength(Particle&) const {
cout << "Process2::GetInteractionLength" << endl;
- return 3_g / (1_cm * 1_cm);
+ return 20_g / (1_cm * 1_cm);
}
};
@@ -111,15 +147,17 @@ public:
globalCount++;
}
- template <typename Particle>
- inline EProcessReturn DoInteraction(Particle&) const {
+ template <typename TView>
+ inline EProcessReturn DoInteraction(TView& v) const {
+ checkInteract |= 4;
+ for (int i = 0; i < nData; ++i) v.parent().data_[i] *= 1.01;
cout << "Process3::DoInteraction" << endl;
return EProcessReturn::eOk;
}
template <typename Particle>
GrammageType GetInteractionLength(Particle&) const {
cout << "Process3::GetInteractionLength" << endl;
- return 1_g / (1_cm * 1_cm);
+ return 30_g / (1_cm * 1_cm);
}
};
@@ -135,11 +173,14 @@ public:
template <typename D, typename T>
inline EProcessReturn DoContinuous(D& d, T&) const {
- for (int i = 0; i < nData; ++i) { d.p[i] /= 1.2; }
+ std::cout << "Base::DoContinuous" << std::endl;
+ checkCont |= 8;
+ for (int i = 0; i < nData; ++i) { d.data_[i] /= 1.2; }
return EProcessReturn::eOk;
}
- template <typename Particle>
- EProcessReturn DoInteraction(Particle&) const {
+ template <typename TView>
+ EProcessReturn DoInteraction(TView&) const {
+ checkInteract |= 8;
return EProcessReturn::eOk;
}
};
@@ -156,8 +197,28 @@ public:
TimeType GetLifetime(Particle&) const {
return 1_s;
}
+ template <typename TView>
+ EProcessReturn DoDecay(TView&) const {
+ checkDecay |= 1;
+ return EProcessReturn::eOk;
+ }
+};
+
+class Decay2 : public DecayProcess<Decay2> {
+
+public:
+ Decay2(const int) {
+ cout << "Decay2()" << endl;
+ globalCount++;
+ }
+
template <typename Particle>
- EProcessReturn DoDecay(Particle&) const {
+ TimeType GetLifetime(Particle&) const {
+ return 2_s;
+ }
+ template <typename TView>
+ EProcessReturn DoDecay(TView&) const {
+ checkDecay |= 2;
return EProcessReturn::eOk;
}
};
@@ -178,10 +239,17 @@ public:
struct DummyStack {};
struct DummyData {
- double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+ double data_[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
};
struct DummyTrajectory {};
+struct DummyView {
+ DummyView(DummyData& p)
+ : p_(p) {}
+ DummyData& p_;
+ DummyData& parent() { return p_; }
+};
+
TEST_CASE("Process Sequence", "[Process Sequence]") {
SECTION("Check construction") {
@@ -195,7 +263,17 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
Process4 m4(3);
CHECK(globalCount == 4);
- [[maybe_unused]] auto sequence = m1 << m2 << m3 << m4;
+ auto sequence1 = process::sequence(m1, m2, m3, m4);
+ CHECK(is_process_sequence_v<decltype(sequence1)> == true);
+ CHECK(is_process_sequence_v<decltype(m2)> == false);
+ CHECK(is_switch_process_sequence_v<decltype(sequence1)> == false);
+ CHECK(is_switch_process_sequence_v<decltype(m2)> == false);
+
+ auto sequence2 = process::sequence(m1, m2, m3);
+ CHECK(is_process_sequence_v<decltype(sequence2)> == true);
+
+ auto sequence3 = process::sequence(m4);
+ CHECK(is_process_sequence_v<decltype(sequence3)> == true);
}
SECTION("interaction length") {
@@ -206,11 +284,13 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
DummyData particle;
- auto sequence2 = cp1 << m2 << m3;
- GrammageType const tot = sequence2.GetTotalInteractionLength(particle);
- InverseGrammageType const tot_inv =
- sequence2.GetTotalInverseInteractionLength(particle);
+ auto sequence2 = sequence(cp1, m2, m3);
+ GrammageType const tot = sequence2.GetInteractionLength(particle);
+ InverseGrammageType const tot_inv = sequence2.GetInverseInteractionLength(particle);
cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
+
+ CHECK(tot / 1_g * square(1_cm) == 12);
+ CHECK(tot_inv * 1_g / square(1_cm) == 1. / 12);
globalCount = 0;
}
@@ -223,22 +303,24 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
DummyData particle;
- auto sequence2 = cp1 << m2 << m3 << d3;
- TimeType const tot = sequence2.GetTotalLifetime(particle);
- InverseTimeType const tot_inv = sequence2.GetTotalInverseLifetime(particle);
+ auto sequence2 = sequence(cp1, m2, m3, d3);
+ TimeType const tot = sequence2.GetLifetime(particle);
+ InverseTimeType const tot_inv = sequence2.GetInverseLifetime(particle);
cout << "lambda_tot=" << tot << "; lambda_tot_inv=" << tot_inv << endl;
+ CHECK(tot / 1_s == 1);
+ CHECK(tot_inv * 1_s == 1.);
globalCount = 0;
}
SECTION("sectionTwo") {
globalCount = 0;
- ContinuousProcess1 cp1(0);
- ContinuousProcess2 cp2(1);
- Process2 m2(2);
- Process3 m3(3);
+ ContinuousProcess1 cp1(0); // += 0.933
+ ContinuousProcess2 cp2(1); // += 0.111
+ Process2 m2(2); // /= 1.1
+ Process3 m3(3); // *= 1.01
- auto sequence2 = cp1 << m2 << m3 << cp2;
+ auto sequence2 = sequence(cp1, m2, m3, cp2);
DummyData particle;
DummyTrajectory track;
@@ -247,15 +329,18 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
globalCount = 0;
cout << "-->docont" << endl;
- sequence2.DoContinuous(particle, track);
- cout << "-->dodisc" << endl;
- cout << "-->done" << endl;
+ // validation data
+ double test_data[nData] = {0};
const int nLoop = 5;
cout << "Running loop with n=" << nLoop << endl;
- for (int i = 0; i < nLoop; ++i) { sequence2.DoContinuous(particle, track); }
+ for (int iLoop = 0; iLoop < nLoop; ++iLoop) {
+ for (int i = 0; i < nData; ++i) { test_data[i] += 0.933 + 0.111; }
+ sequence2.DoContinuous(particle, track);
+ }
for (int i = 0; i < nData; i++) {
- cout << "data[" << i << "]=" << particle.p[i] << endl;
+ cout << "data_[" << i << "]=" << particle.data_[i] << endl;
+ CHECK(particle.data_[i] == Approx(test_data[i]).margin(1e-9));
}
cout << "done" << endl;
}
@@ -266,7 +351,7 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
Stack1 s1(1);
Stack1 s2(2);
- auto sequence = s1 << s2;
+ auto sequence = process::sequence(s1, s2);
DummyStack stack;
@@ -275,17 +360,120 @@ TEST_CASE("Process Sequence", "[Process Sequence]") {
CHECK(s1.GetCount() == 20);
CHECK(s2.GetCount() == 10);
+
+ ContinuousProcess2 cp2(1); // += 0.111
+ Process2 m2(2); // /= 1.1
+ auto sequence2 = process::sequence(cp2, m2);
+ auto sequence3 = process::sequence(cp2, m2, s1);
+
+ CHECK(is_process_sequence_v<decltype(sequence2)> == true);
+ CHECK(is_process_sequence_v<decltype(sequence3)> == true);
+ CHECK(contains_stack_process_v<decltype(sequence2)> == false);
+ CHECK(contains_stack_process_v<decltype(sequence3)> == true);
}
}
-/*
- Note: there is a fine-grained dedicated test-suite for SwitchProcess
- in Processes/SwitchProcess/testSwtichProcess
- */
-TEST_CASE("SwitchProcess") {
- globalCount = 0;
- Process1 p1(0);
- Process2 p2(1);
- switch_process::SwitchProcess s(p1, p2, 10_GeV);
- CHECK(is_switch_process_v<decltype(s)>);
+TEST_CASE("Switch Process Sequence", "[Process Sequence]") {
+
+ SECTION("Check construction") {
+
+ struct TestSelect {
+ corsika::process::SwitchResult operator()(const DummyData& p) const {
+ std::cout << "TestSelect data=" << p.data_[0] << std::endl;
+ if (p.data_[0] > 0) return corsika::process::SwitchResult::First;
+ return corsika::process::SwitchResult::Second;
+ }
+ };
+ TestSelect select;
+
+ auto sequence1 = process::sequence(Process1(0), ContinuousProcess2(0), Decay1(0));
+ auto sequence2 = process::sequence(ContinuousProcess3(0), Process2(0), Decay2(0));
+
+ auto sequence =
+ process::sequence(ContinuousProcess1(0), Process3(0),
+ SwitchProcessSequence(sequence1, sequence2, select));
+
+ auto sequence_alt = process::sequence(
+ ContinuousProcess1(0), Process3(0),
+ process::select(process::sequence(Process1(0), ContinuousProcess2(0), Decay1(0)),
+ process::sequence(ContinuousProcess3(0), Process2(0), Decay2(0)),
+ select));
+
+ // check that same process sequence can be build in different ways
+ CHECK(typeid(sequence) == typeid(sequence_alt));
+ CHECK(is_process_sequence_v<decltype(sequence)> == true);
+ CHECK(is_process_sequence_v<decltype(
+ SwitchProcessSequence(sequence1, sequence2, select))> == true);
+
+ DummyData particle;
+ DummyTrajectory track;
+ DummyView view(particle);
+
+ checkDecay = 0;
+ checkInteract = 0;
+ checkSec = 0;
+ checkCont = 0;
+ particle.data_[0] = 100; // data positive
+ sequence.DoContinuous(particle, track);
+ CHECK(checkInteract == 0);
+ CHECK(checkDecay == 0);
+ CHECK(checkCont == 0b011);
+ CHECK(checkSec == 0);
+
+ checkDecay = 0;
+ checkInteract = 0;
+ checkSec = 0;
+ checkCont = 0;
+ particle.data_[0] = -100; // data negative
+ sequence_alt.DoContinuous(particle, track);
+ CHECK(checkInteract == 0);
+ CHECK(checkDecay == 0);
+ CHECK(checkCont == 0b101);
+ CHECK(checkSec == 0);
+
+ // 1/(30g/cm2) is Process3
+ corsika::units::si::InverseGrammageType lambda_select = .9 / 30. * square(1_cm) / 1_g;
+ corsika::units::si::InverseTimeType time_select = 0.1 / second;
+
+ checkDecay = 0;
+ checkInteract = 0;
+ checkSec = 0;
+ checkCont = 0;
+ particle.data_[0] = 100; // data positive
+ sequence.SelectInteraction(view, lambda_select);
+ sequence.SelectDecay(view, time_select);
+ CHECK(checkInteract == 0b100); // this is Process3
+ CHECK(checkDecay == 0b001); // this is Decay1
+ CHECK(checkCont == 0);
+ CHECK(checkSec == 0);
+ lambda_select = 1.01 / 30. * square(1_cm) / 1_g;
+ checkInteract = 0;
+ sequence.SelectInteraction(view, lambda_select);
+ CHECK(checkInteract == 0b001); // this is Process1
+
+ checkDecay = 0;
+ checkInteract = 0;
+ checkSec = 0;
+ checkCont = 0;
+ particle.data_[0] = -100; // data negative
+ sequence.SelectInteraction(view, lambda_select);
+ sequence.SelectDecay(view, time_select);
+ CHECK(checkInteract == 0b010); // this is Process2
+ CHECK(checkDecay == 0b010); // this is Decay2
+ CHECK(checkCont == 0);
+ CHECK(checkSec == 0);
+
+ checkDecay = 0;
+ checkInteract = 0;
+ checkSec = 0;
+ checkCont = 0;
+ particle.data_[0] = -100; // data negative
+ sequence.DoSecondaries(view);
+ Stack1 stack(0);
+ sequence.DoStack(stack);
+ CHECK(checkInteract == 0);
+ CHECK(checkDecay == 0);
+ CHECK(checkCont == 0);
+ CHECK(checkSec == 0);
+ }
}
diff --git a/Framework/ProcessSequence/testSwitchProcessSequence.cc b/Framework/ProcessSequence/testSwitchProcessSequence.cc
new file mode 100644
index 0000000000000000000000000000000000000000..a778f2efb5dd587b3fc2341dbb7cb5e205c671ab
--- /dev/null
+++ b/Framework/ProcessSequence/testSwitchProcessSequence.cc
@@ -0,0 +1,247 @@
+/*
+ * (c) Copyright 2018 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.
+ */
+
+#include <corsika/process/switch_process/SwitchProcess.h>
+#include <corsika/stack/SecondaryView.h>
+#include <corsika/stack/Stack.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <catch2/catch.hpp>
+
+#include <algorithm>
+#include <random>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units::si;
+
+class TestStackData {
+
+public:
+ // these functions are needed for the Stack interface
+ void Clear() { fData.clear(); }
+ unsigned int GetSize() const { return fData.size(); }
+ unsigned int GetCapacity() const { return fData.size(); }
+ void Copy(int i1, int i2) { fData[i2] = fData[i1]; }
+ void Swap(int i1, int i2) { std::swap(fData[i1], fData[i2]); }
+
+ // custom data access function
+ void SetData(unsigned int i, HEPEnergyType v) { fData[i] = v; }
+ HEPEnergyType GetData(const unsigned int i) const { return fData[i]; }
+
+ // these functions are also needed by the Stack interface
+ void IncrementSize() { fData.resize(fData.size() + 1); }
+ void DecrementSize() {
+ if (fData.size() > 0) { fData.pop_back(); }
+ }
+
+ // custom private data section
+private:
+ std::vector<HEPEnergyType> fData;
+};
+
+/**
+ * From static_cast of a StackIterator over entries in the
+ * TestStackData class you get and object of type
+ * TestParticleInterface defined here
+ *
+ * It provides Get/Set methods to read and write data to the "Data"
+ * storage of TestStackData obtained via
+ * "StackIteratorInterface::GetStackData()", given the index of the
+ * iterator "StackIteratorInterface::GetIndex()"
+ *
+ */
+template <typename StackIteratorInterface>
+class TestParticleInterface
+ : public corsika::stack::ParticleBase<StackIteratorInterface> {
+
+public:
+ using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
+ using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
+
+ /*
+ The SetParticleData methods are called for creating new entries
+ on the stack. You can specifiy various parametric versions to
+ perform this task:
+ */
+
+ // default version for particle-creation from input data
+ void SetParticleData(const std::tuple<HEPEnergyType> v) { SetEnergy(std::get<0>(v)); }
+ void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/,
+ std::tuple<HEPEnergyType> v) {
+ SetEnergy(std::get<0>(v));
+ }
+
+ // here are the fundamental methods for access to TestStackData data
+ void SetEnergy(HEPEnergyType v) { GetStackData().SetData(GetIndex(), v); }
+ HEPEnergyType GetEnergy() const { return GetStackData().GetData(GetIndex()); }
+};
+
+using SimpleStack = corsika::stack::Stack<TestStackData, TestParticleInterface>;
+
+// see issue 161
+#if defined(__clang__)
+using StackTestView = corsika::stack::SecondaryView<TestStackData, TestParticleInterface>;
+#elif defined(__GNUC__) || defined(__GNUG__)
+using StackTestView = corsika::stack::MakeView<SimpleStack>::type;
+#endif
+
+auto constexpr kgMSq = 1_kg / (1_m * 1_m);
+
+template <int N>
+struct DummyProcess : InteractionProcess<DummyProcess<N>> {
+
+ template <typename TParticle>
+ corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const {
+ return N * kgMSq;
+ }
+
+ template <typename TSecondaries>
+ corsika::process::EProcessReturn DoInteraction(TSecondaries& vSec) {
+ // to figure out which process was selected in the end, we produce N
+ // secondaries for DummyProcess<N>
+
+ for (int i = 0; i < N; ++i) {
+ vSec.AddSecondary(std::tuple<HEPEnergyType>{vSec.GetEnergy() / N});
+ }
+
+ return EProcessReturn::eOk;
+ }
+};
+
+using DummyLowEnergyProcess = DummyProcess<1>;
+using DummyHighEnergyProcess = DummyProcess<2>;
+using DummyAdditionalProcess = DummyProcess<3>;
+
+TEST_CASE("SwitchProcess from InteractionProcess") {
+ DummyLowEnergyProcess low;
+ DummyHighEnergyProcess high;
+ DummyAdditionalProcess proc;
+
+ switch_process::SwitchProcess switchProcess(low, high, 1_TeV);
+ auto seq = switchProcess << proc;
+
+ SimpleStack stack;
+
+ SECTION("low energy") {
+ stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
+ auto p = stack.GetNextParticle();
+
+ // low energy process returns 1 kg/m²
+ SECTION("interaction length") {
+ CHECK(switchProcess.GetInteractionLength(p) / kgMSq == Approx(1));
+ CHECK(seq.GetInteractionLength(p) / kgMSq == Approx(3. / 4));
+ }
+ }
+
+ SECTION("high energy") {
+ stack.AddParticle(std::tuple<HEPEnergyType>{4_TeV});
+ auto p = stack.GetNextParticle();
+
+ // high energy process returns 2 kg/m²
+ SECTION("interaction length") {
+ CHECK(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2));
+ CHECK(seq.GetInteractionLength(p) / kgMSq == Approx(6. / 5));
+ }
+
+ // high energy process creates 2 secondaries
+ SECTION("SelectInteraction") {
+ typename SimpleStack::ParticleType theParticle =
+ stack.GetNextParticle(); // as in corsika::Cascade
+ StackTestView view(theParticle);
+ auto projectile = view.GetProjectile();
+
+ InverseGrammageType invLambda = 0 / kgMSq;
+ switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
+
+ CHECK(view.getSize() == 2);
+ }
+ }
+}
+
+TEST_CASE("SwitchProcess from ProcessSequence") {
+ DummyProcess<1> innerA;
+ DummyProcess<2> innerB;
+ DummyProcess<3> outer;
+ DummyProcess<4> additional;
+
+ auto seq = innerA << innerB;
+
+ switch_process::SwitchProcess switchProcess(seq, outer, 1_TeV);
+ auto completeSeq = switchProcess << additional;
+
+ SimpleStack stack;
+
+ SECTION("low energy") {
+ stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
+ auto p = stack.GetNextParticle();
+
+ SECTION("interaction length") {
+ CHECK(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2. / 3));
+ CHECK(completeSeq.GetInteractionLength(p) / kgMSq == Approx(4. / 7));
+ }
+
+ SECTION("SelectInteraction") {
+ std::vector<int> numberOfSecondaries;
+
+ for (int i = 0; i < 1000; ++i) {
+ typename SimpleStack::ParticleType theParticle =
+ stack.GetNextParticle(); // as in corsika::Cascade
+ StackTestView view(theParticle);
+ auto projectile = view.GetProjectile();
+
+ double r = i / 1000.;
+ InverseGrammageType invLambda = r * 7. / 4 / kgMSq;
+
+ InverseGrammageType accumulator = 0 / kgMSq;
+ completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
+
+ numberOfSecondaries.push_back(view.getSize());
+ }
+
+ auto const mean =
+ std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
+ numberOfSecondaries.size();
+ CHECK(mean == Approx(12. / 7.).margin(0.01));
+ }
+ }
+
+ SECTION("high energy") {
+ stack.AddParticle(std::tuple<HEPEnergyType>{3.0_TeV});
+ auto p = stack.GetNextParticle();
+
+ SECTION("interaction length") {
+ CHECK(switchProcess.GetInteractionLength(p) / kgMSq == Approx(3));
+ CHECK(completeSeq.GetInteractionLength(p) / kgMSq == Approx(12. / 7.));
+ }
+
+ SECTION("SelectInteraction") {
+ std::vector<int> numberOfSecondaries;
+
+ for (int i = 0; i < 1000; ++i) {
+ typename SimpleStack::ParticleType theParticle =
+ stack.GetNextParticle(); // as in corsika::Cascade
+ StackTestView view(theParticle);
+ auto projectile = view.GetProjectile();
+
+ double r = i / 1000.;
+ InverseGrammageType invLambda = r * 7. / 12. / kgMSq;
+
+ InverseGrammageType accumulator = 0 / kgMSq;
+ completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
+
+ numberOfSecondaries.push_back(view.getSize());
+ }
+
+ auto const mean =
+ std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
+ numberOfSecondaries.size();
+ CHECK(mean == Approx(24. / 7.).margin(0.01));
+ }
+ }
+}
diff --git a/Processes/SwitchProcess/testSwitchProcess.cc b/Framework/ProcessSequence/testSwitchProcessSequence.h
similarity index 94%
rename from Processes/SwitchProcess/testSwitchProcess.cc
rename to Framework/ProcessSequence/testSwitchProcessSequence.h
index 9ee148e1855edc2b45c81594425e8786e2434d09..709dc8510684c9e7c4dc8ebbb91fa20611ac17ad 100644
--- a/Processes/SwitchProcess/testSwitchProcess.cc
+++ b/Framework/ProcessSequence/testSwitchProcessSequence.h
@@ -137,19 +137,6 @@ TEST_CASE("SwitchProcess from InteractionProcess") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(1));
REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(3. / 4));
}
-
- // low energy process creates 1 secondary
- //~ SECTION("SelectInteraction") {
- //~ typename SimpleStack::ParticleType theParticle =
- //~ stack.GetNextParticle(); // as in corsika::Cascade
- //~ StackTestView view(theParticle);
- //~ auto projectile = view.GetProjectile();
-
- //~ InverseGrammageType invLambda = 0 / kgMSq;
- //~ switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
-
- //~ REQUIRE(view.GetSize() == 1);
- //~ }
}
SECTION("high energy") {
diff --git a/Framework/Random/CMakeLists.txt b/Framework/Random/CMakeLists.txt
index f9c72763ca232ed24db5ff5dea1a407eb24b6a7c..f909f77f1f32c88f77a52482839af3cf43e2949a 100644
--- a/Framework/Random/CMakeLists.txt
+++ b/Framework/Random/CMakeLists.txt
@@ -20,9 +20,8 @@ CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKArandom ${CORSIKArandom_NAMESPACE} ${CO
target_link_libraries (
CORSIKArandom
- INTERFACE
CORSIKAutilities
- CORSIKAunits
+ CORSIKAlogging
)
target_include_directories (
@@ -50,4 +49,5 @@ target_link_libraries (
testRandom
CORSIKArandom
CORSIKAtesting
+ CORSIKAunits
)
diff --git a/Framework/Random/RNGManager.cc b/Framework/Random/RNGManager.cc
index 238b38636e4ef5ec6ab9a3e219ba044d7a503845..8bbcf15dd391a744d69e017f7570e5163a0308a2 100644
--- a/Framework/Random/RNGManager.cc
+++ b/Framework/Random/RNGManager.cc
@@ -7,6 +7,8 @@
*/
#include <corsika/random/RNGManager.h>
+#include <corsika/logging/Logging.h>
+
#include <sstream>
void corsika::random::RNGManager::RegisterRandomStream(std::string const& pStreamName) {
@@ -42,14 +44,21 @@ std::stringstream corsika::random::RNGManager::dumpState() const {
}
void corsika::random::RNGManager::SeedAll(uint64_t vSeed) {
- for (auto& entry : rngs) { entry.second.seed(vSeed++); }
+ for (auto& entry : rngs) {
+ auto seed = vSeed++;
+ C8LOG_TRACE("Random seed stream {} seed {}", entry.first, seed);
+ entry.second.seed(seed);
+ }
}
void corsika::random::RNGManager::SeedAll() {
std::random_device rd;
-
+ std::seed_seq sseq{rd(), rd(), rd(), rd(), rd(), rd()};
for (auto& entry : rngs) {
- std::seed_seq sseq{rd(), rd(), rd(), rd(), rd(), rd()};
- entry.second.seed(sseq);
+ std::vector<std::uint32_t> seeds(1);
+ sseq.generate(seeds.begin(), seeds.end());
+ std::uint32_t seed = seeds[0];
+ C8LOG_TRACE("Random seed stream {} seed {}", entry.first, seed);
+ entry.second.seed(seed);
}
}
diff --git a/Framework/Random/UniformRealDistribution.h b/Framework/Random/UniformRealDistribution.h
index 7e092076da1ea7d179df5313a9d767f9e38faf9f..0c63e85428a1e2556a81995b386ba46493546866 100644
--- a/Framework/Random/UniformRealDistribution.h
+++ b/Framework/Random/UniformRealDistribution.h
@@ -8,7 +8,6 @@
#pragma once
-#include <corsika/units/PhysicalUnits.h>
#include <random>
namespace corsika::random {
diff --git a/Framework/StackInterface/SecondaryView.h b/Framework/StackInterface/SecondaryView.h
index dcdcab8127e6217bfc4f75d68787ae92b80a9178..709a8c4feb658659f9369f938490b7101e91b1da 100644
--- a/Framework/StackInterface/SecondaryView.h
+++ b/Framework/StackInterface/SecondaryView.h
@@ -168,8 +168,9 @@ namespace corsika::stack {
* SecondaryView is derived from. This projectile should not be
* used to modify the Stack!
*/
- ConstStackIteratorValue parent() const {
- return ConstStackIteratorValue(inner_stack_, projectile_index_);
+ StackIteratorValue parent()
+ const { // todo: check if this can't be ConstStackIteratorValue
+ return StackIteratorValue(inner_stack_, projectile_index_);
}
/**
diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h
index e166fccba8a8acdb957222b08b42bdcf63a87193..85ec6418fe0898573b3db16d9d52e867404d3b76 100644
--- a/Framework/StackInterface/Stack.h
+++ b/Framework/StackInterface/Stack.h
@@ -301,9 +301,11 @@ namespace corsika::stack {
/**
* check if this particle was already deleted
*/
- bool isDeleted(const StackIterator& p) { return isDeleted(p.GetIndex()); }
+ bool isDeleted(const StackIterator& p) const { return isDeleted(p.GetIndex()); }
bool isDeleted(const ConstStackIterator& p) const { return isDeleted(p.GetIndex()); }
- bool isDeleted(const ParticleInterfaceType& p) { return isDeleted(p.GetIterator()); }
+ bool isDeleted(const ParticleInterfaceType& p) const {
+ return isDeleted(p.GetIterator());
+ }
/**
* Function to ultimatively remove the last entry from the stack,
diff --git a/Framework/StackInterface/StackIteratorInterface.h b/Framework/StackInterface/StackIteratorInterface.h
index 2229c17bc395a2f12b83aae0ae5c89162c95fdd3..6c4901084843c568016450205435884d3e6bb972 100644
--- a/Framework/StackInterface/StackIteratorInterface.h
+++ b/Framework/StackInterface/StackIteratorInterface.h
@@ -64,7 +64,7 @@ namespace corsika::stack {
For two examples see stack_example.cc, or the
corsika::processes::sibyll::SibStack class
- */
+ **/
template <typename TStackData, template <typename> typename TParticleInterface,
typename StackType = Stack<TStackData, TParticleInterface>>
@@ -105,6 +105,10 @@ namespace corsika::stack {
StackIteratorInterface() = delete;
public:
+ StackIteratorInterface(StackIteratorInterface&& rhs)
+ : index_(std::move(rhs.index_))
+ , data_(std::move(rhs.data_)) {}
+
StackIteratorInterface(StackIteratorInterface const& vR)
: index_(vR.index_)
, data_(vR.data_) {}
@@ -118,7 +122,7 @@ namespace corsika::stack {
/** iterator must always point to data, with an index:
@param data reference to the stack [rw]
@param index index on stack
- */
+ **/
StackIteratorInterface(StackType& data, const unsigned int index)
: index_(index)
, data_(&data) {}
@@ -129,7 +133,7 @@ namespace corsika::stack {
@param args variadic list of data to initialize stack entry, this must be
consistent with the definition of the user-provided
ParticleInterfaceType::SetParticleData(...) function
- */
+ **/
template <typename... Args>
StackIteratorInterface(StackType& data, const unsigned int index, const Args... args)
: index_(index)
@@ -146,7 +150,7 @@ namespace corsika::stack {
@param args variadic list of data to initialize stack entry, this must be
consistent with the definition of the user-provided
ParticleInterfaceType::SetParticleData(...) function
- */
+ **/
template <typename... Args>
StackIteratorInterface(StackType& data, const unsigned int index,
StackIteratorInterface& parent, const Args... args)
@@ -160,7 +164,7 @@ namespace corsika::stack {
public:
/** @name Iterator interface
@{
- */
+ **/
StackIteratorInterface& operator++() {
do {
++index_;
@@ -195,14 +199,15 @@ namespace corsika::stack {
/**
* Convert iterator to value type, where value type is the user-provided particle
* readout class
- */
+ **/
ParticleInterfaceType& operator*() {
return static_cast<ParticleInterfaceType&>(*this);
}
+
/**
* Convert iterator to const value type, where value type is the user-provided
* particle readout class
- */
+ **/
const ParticleInterfaceType& operator*() const {
return static_cast<const ParticleInterfaceType&>(*this);
}
@@ -212,7 +217,7 @@ namespace corsika::stack {
/**
* @name Stack data access
* @{
- */
+ **/
/// Get current particle index
inline unsigned int GetIndex() const { return index_; }
/// Get current particle Stack object
@@ -234,7 +239,7 @@ namespace corsika::stack {
@class ConstStackIteratorInterface
This is the iterator class for const-access to stack data
- */
+ **/
template <typename TStackData, template <typename> typename TParticleInterface,
typename StackType = Stack<TStackData, TParticleInterface>>
@@ -275,6 +280,10 @@ namespace corsika::stack {
ConstStackIteratorInterface() = delete;
public:
+ ConstStackIteratorInterface(ConstStackIteratorInterface&& rhs)
+ : index_(std::move(rhs.index_))
+ , data_(std::move(rhs.data_)) {}
+
ConstStackIteratorInterface(const StackType& data, const unsigned int index)
: index_(index)
, data_(&data) {}
@@ -290,13 +299,13 @@ namespace corsika::stack {
\endverbatim
See documentation of StackIteratorInterface for more details.
- */
+ **/
bool isDeleted() const { return GetStack().isDeleted(*this); }
public:
/** @name Iterator interface
- */
+ **/
///@{
ConstStackIteratorInterface& operator++() {
do {
@@ -339,7 +348,7 @@ namespace corsika::stack {
protected:
/** @name Stack data access
Only the const versions for read-only access
- */
+ **/
///@{
inline unsigned int GetIndex() const { return index_; }
inline const StackType& GetStack() const { return *data_; }
diff --git a/Framework/Utilities/CMakeLists.txt b/Framework/Utilities/CMakeLists.txt
index 5618ab4a8cf35be1e36630eac8d3f8d0ab627e95..95170d4a35b01b6f44eaf6aff11f47ebf0881e00 100644
--- a/Framework/Utilities/CMakeLists.txt
+++ b/Framework/Utilities/CMakeLists.txt
@@ -37,6 +37,7 @@ set (
UTILITIES_DEPENDS
CORSIKAgeometry
CORSIKAunits
+ CORSIKAlogging
C8::ext::boost # so far only for MetaProgramming
C8::ext::eigen3 # for COMboost
)
diff --git a/Framework/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc
index c6688199ab5d543cfd1670b93e7416cbd8a961f7..5bb4be9da2a2784dc70efa0a7b6b5b1233d3b3f4 100644
--- a/Framework/Utilities/COMBoost.cc
+++ b/Framework/Utilities/COMBoost.cc
@@ -12,6 +12,7 @@
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
#include <corsika/utl/sgn.h>
+#include <corsika/logging/Logging.h>
#include <cmath>
@@ -38,9 +39,8 @@ COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pproj
setBoost(coshEta, sinhEta);
- std::cout << "COMBoost (1-beta)=" << 1 - sinhEta / coshEta << " gamma=" << coshEta
- << std::endl;
- std::cout << " det = " << boost_.determinant() - 1 << std::endl;
+ C8LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta, coshEta,
+ boost_.determinant() - 1);
}
COMBoost::COMBoost(geometry::Vector<units::si::hepmomentum_d> const& momentum,
diff --git a/Framework/Utilities/COMBoost.h b/Framework/Utilities/COMBoost.h
index 0c78d49acb6ef06dd41c1e64c7f5048172c1e34a..73b7c93001321dff22902477dc701aee043e9040 100644
--- a/Framework/Utilities/COMBoost.h
+++ b/Framework/Utilities/COMBoost.h
@@ -11,6 +11,7 @@
#include <corsika/geometry/CoordinateSystem.h>
#include <corsika/geometry/FourVector.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/logging/Logging.h>
#include <Eigen/Dense>
@@ -69,10 +70,10 @@ namespace corsika::utl {
Eigen::Vector2d com;
com << (Ecm * (1 / 1_GeV)), (pCM.eVector(2) * (1 / 1_GeV).magnitude());
- std::cout << "COMBoost::fromCoM Ecm=" << Ecm / 1_GeV << " GeV, "
- << " pcm = " << pCM / 1_GeV << " (norm = " << pCM.norm() / 1_GeV
- << " GeV), invariant mass = " << p.GetNorm() / 1_GeV << " GeV"
- << std::endl;
+ C8LOG_TRACE(
+ "COMBoost::fromCoM Ecm={} GeV"
+ " pcm={} GeV (norm = {} GeV), invariant mass={} GeV",
+ Ecm / 1_GeV, pCM / 1_GeV, pCM.norm() / 1_GeV, p.GetNorm() / 1_GeV);
auto const boostedZ = inverseBoost_ * com;
auto const E_lab = boostedZ(0) * 1_GeV;
@@ -84,10 +85,11 @@ namespace corsika::utl {
FourVector f(E_lab, pLab);
- std::cout << "COMBoost::fromCoM --> Elab=" << E_lab / 1_GeV << "GeV, "
- << " plab = " << pLab.GetComponents() << " (norm =" << pLab.norm() / 1_GeV
- << " GeV), invariant mass = " << f.GetNorm() / 1_GeV << " GeV"
- << std::endl;
+ C8LOG_TRACE("COMBoost::fromCoM --> Elab={} GeV",
+ " plab={} GeV (norm={} GeV) "
+ " GeV), invariant mass = {}",
+ E_lab / 1_GeV, f.GetNorm() / 1_GeV, pLab.GetComponents(),
+ pLab.norm() / 1_GeV);
return f;
}
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/CMakeLists.txt b/Processes/CMakeLists.txt
index 5deaf153c5944aa5fe0800db28809ff5131f512b..4e9f59df367346d51a9072e6a7bfe95ec87fb4c1 100644
--- a/Processes/CMakeLists.txt
+++ b/Processes/CMakeLists.txt
@@ -2,7 +2,6 @@
add_subdirectory (AnalyticProcessors)
add_subdirectory (ExampleProcessors)
-add_subdirectory (NullModel)
# tracking
add_subdirectory (TrackingLine)
# hadron interaction models
@@ -31,13 +30,11 @@ add_subdirectory (ParticleCut)
add_subdirectory (OnShellCheck)
# meta-processes
add_subdirectory (InteractionCounter)
-add_subdirectory (SwitchProcess)
##########################################
# add_custom_target(CORSIKAprocesses)
add_library (CORSIKAprocesses INTERFACE)
-add_dependencies (CORSIKAprocesses ProcessNullModel)
add_dependencies (CORSIKAprocesses ProcessSibyll)
add_dependencies (CORSIKAprocesses ProcessProposal)
if (Pythia8_FOUND)
diff --git a/Processes/EnergyLoss/CMakeLists.txt b/Processes/EnergyLoss/CMakeLists.txt
index 62fca1ca017b3dbea3a4a331381596b87931ac01..3df3978d5fdb4a23b3548e5c68b8a55f6c11f4d8 100644
--- a/Processes/EnergyLoss/CMakeLists.txt
+++ b/Processes/EnergyLoss/CMakeLists.txt
@@ -32,6 +32,7 @@ target_link_libraries (
CORSIKAgeometry
CORSIKAenvironment
CORSIKAsetup
+ CORSIKAlogging
)
target_include_directories (
diff --git a/Processes/EnergyLoss/EnergyLoss.cc b/Processes/EnergyLoss/EnergyLoss.cc
index 69dcd7877e49b760980b15056c657a0468bf546f..336294d8e97ee9be5fdc3c82b46b30365849e826 100644
--- a/Processes/EnergyLoss/EnergyLoss.cc
+++ b/Processes/EnergyLoss/EnergyLoss.cc
@@ -13,6 +13,8 @@
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/logging/Logging.h>
+
#include <corsika/geometry/Line.h>
#include <cmath>
@@ -85,19 +87,18 @@ HEPEnergyType EnergyLoss::BetheBloch(SetupParticle const& p, GrammageType const
double const beta2 = (gamma2 - 1) / gamma2; // 1-1/gamma2 (1-1/gamma)*(1+1/gamma);
// (gamma_2-1)/gamma_2 = (1-1/gamma2);
double constexpr c2 = 1; // HEP convention here c=c2=1
- cout << "BetheBloch beta2=" << beta2 << " gamma2=" << gamma2 << endl;
+ C8LOG_DEBUG("BetheBloch beta2={}, gamma2={}", beta2, gamma2);
[[maybe_unused]] double const eta2 = beta2 / (1 - beta2);
HEPMassType const Wmax =
2 * me * c2 * beta2 * gamma2 / (1 + 2 * gamma * me / m + me2 / m2);
// approx, but <<1% HEPMassType const Wmax = 2*me*c2*beta2*gamma2; for HEAVY
// PARTICLES Wmax ~ 2me v2 for non-relativistic particles
- cout << "BetheBloch Wmax=" << Wmax << endl;
+ C8LOG_DEBUG("BetheBloch Wmax={}", Wmax);
// Sternheimer parameterization, density corrections towards high energies
// NOTE/TODO: when Cbar is 0 it needs to be approximated from parameterization ->
// MISSING
- cout << "BetheBloch p.GetMomentum().GetNorm()/m=" << p.GetMomentum().GetNorm() / m
- << endl;
+ C8LOG_DEBUG("BetheBloch p.GetMomentum().GetNorm()/m{}=", p.GetMomentum().GetNorm() / m);
double const x = log10(p.GetMomentum().GetNorm() / m);
double delta = 0;
if (x >= x1) {
@@ -107,7 +108,7 @@ HEPEnergyType EnergyLoss::BetheBloch(SetupParticle const& p, GrammageType const
} else if (x < x0) { // and IF conductor (otherwise, this is 0)
delta = delta0 * pow(100, 2 * (x - x0));
}
- cout << "BetheBloch delta=" << delta << endl;
+ C8LOG_DEBUG("BetheBloch delta={}", delta);
// with further low energies correction, accurary ~1% down to beta~0.05 (1MeV for p)
@@ -141,9 +142,6 @@ HEPEnergyType EnergyLoss::BetheBloch(SetupParticle const& p, GrammageType const
double const y2 = Z * Z * alpha * alpha / beta2;
double const bloch = -y2 * (1.202 - y2 * (1.042 - 0.855 * y2 + 0.343 * y2 * y2));
- // cout << "BetheBloch Erel=" << Erel << " barkas=" << barkas << " bloch=" << bloch <<
- // endl;
-
double const aux = 2 * me * c2 * beta2 * gamma2 * Wmax / (Ieff * Ieff);
return -K * Z2 * ZoverA / beta2 *
(0.5 * log(aux) - beta2 - Cadj / Z - delta / 2 + barkas + bloch) * dX;
@@ -168,24 +166,18 @@ process::EProcessReturn EnergyLoss::DoContinuous(SetupParticle& p, SetupTrack co
GrammageType const dX =
p.GetNode()->GetModelProperties().IntegratedGrammage(t, t.GetLength());
- cout << "EnergyLoss " << p.GetPID() << ", z=" << p.GetChargeNumber()
- << ", dX=" << dX / 1_g * square(1_cm) << "g/cm2" << endl;
+ C8LOG_DEBUG("EnergyLoss pid={}, z={}, dX={} g/cm2", p.GetPID(), p.GetChargeNumber(),
+ 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;
- cout << "EnergyLoss dE=" << dE / 1_MeV << "MeV, "
- << " E=" << E / 1_GeV << "GeV, Ekin=" << Ekin / 1_GeV << ", Enew=" << Enew / 1_GeV
- << "GeV" << endl;
- auto status = process::EProcessReturn::eOk;
- if (E < emCut_) {
- Enew = emCut_;
- status = process::EProcessReturn::eParticleAbsorbed;
- }
+ C8LOG_DEBUG("EnergyLoss dE={} MeV, E={} GeV, Ekin={} GeV, Enew={} GeV", dE / 1_MeV,
+ E / 1_GeV, Ekin / 1_GeV, Enew / 1_GeV);
p.SetEnergy(Enew);
MomentumUpdate(p, Enew);
FillProfile(t, dE);
- return status;
+ return process::EProcessReturn::eOk;
}
LengthType EnergyLoss::MaxStepLength(SetupParticle const& vParticle,
@@ -228,34 +220,35 @@ 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;
- std::cout << "energy deposit of " << -dE << " between " << grammageStart << " and "
- << grammageEnd << std::endl;
+ 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;
- std::cout << "filling bin " << bin << " with weight " << weight << ": " << increment
- << std::endl;
- }
+ 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); }
+ }
- std::cout << "total energy added to histogram: " << energyCount << std::endl;
+ C8LOG_DEBUG("total energy added to histogram: {} ", energyCount);
}
void EnergyLoss::PrintProfile() const {
diff --git a/Processes/LongitudinalProfile/LongitudinalProfile.cc b/Processes/LongitudinalProfile/LongitudinalProfile.cc
index 0e174ee7338dfe2f5fb47bbbf6facfe811fe2ced..f9fd05c1baf2e450fc22e791a60e788a51ceda57 100644
--- a/Processes/LongitudinalProfile/LongitudinalProfile.cc
+++ b/Processes/LongitudinalProfile/LongitudinalProfile.cc
@@ -11,6 +11,8 @@
#include <corsika/particles/ParticleProperties.h>
+#include <corsika/logging/Logging.h>
+
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
@@ -39,6 +41,10 @@ corsika::process::EProcessReturn LongitudinalProfile::DoContinuous(Particle cons
GrammageType const grammageStart = shower_axis_.projectedX(vTrack.GetPosition(0));
GrammageType const grammageEnd = shower_axis_.projectedX(vTrack.GetPosition(1));
+ C8LOG_INFO(
+ "pos1={} m, pos2={}, X={} g/cm2", vTrack.GetPosition(0).GetCoordinates() / 1_m,
+ vTrack.GetPosition(1).GetCoordinates() / 1_m, grammageStart / 1_g * square(1_cm));
+
const int binStart = std::ceil(grammageStart / dX_);
const int binEnd = std::floor(grammageEnd / dX_);
diff --git a/Processes/NullModel/NullModel.cc b/Processes/NullModel/NullModel.cc
deleted file mode 100644
index 1a3afc793fd60c8cd5085c03cb444d3d1a399221..0000000000000000000000000000000000000000
--- a/Processes/NullModel/NullModel.cc
+++ /dev/null
@@ -1,31 +0,0 @@
-/*
- * (c) Copyright 2018 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.
- */
-
-#include <corsika/process/null_model/NullModel.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-using namespace corsika;
-namespace corsika::process::null_model {
-
- NullModel::NullModel(units::si::LengthType maxStepLength)
- : fMaxStepLength(maxStepLength) {}
-
- template <>
- process::EProcessReturn NullModel::DoContinuous(setup::Stack::ParticleType&,
- setup::Trajectory&) const {
- return process::EProcessReturn::eOk;
- }
-
- template <>
- units::si::LengthType NullModel::MaxStepLength(setup::Stack::ParticleType&,
- setup::Trajectory&) const {
- return fMaxStepLength;
- }
-
-} // namespace corsika::process::null_model
diff --git a/Processes/NullModel/NullModel.h b/Processes/NullModel/NullModel.h
deleted file mode 100644
index ba80599f3e67cc9104b28b0ef35921b101841bb3..0000000000000000000000000000000000000000
--- a/Processes/NullModel/NullModel.h
+++ /dev/null
@@ -1,30 +0,0 @@
-/*
- * (c) Copyright 2018 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
-
-#include <corsika/process/BaseProcess.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::process::null_model {
-
- class NullModel : public corsika::process::BaseProcess<NullModel> {
- corsika::units::si::LengthType const fMaxStepLength;
-
- public:
- NullModel(corsika::units::si::LengthType maxStepLength =
- corsika::units::si::meter * std::numeric_limits<double>::infinity());
-
- template <typename Particle, typename Track>
- process::EProcessReturn DoContinuous(Particle&, Track&) const;
-
- template <typename Particle, typename Track>
- corsika::units::si::LengthType MaxStepLength(Particle&, Track&) const;
- };
-
-} // namespace corsika::process::null_model
diff --git a/Processes/NullModel/testNullModel.cc b/Processes/NullModel/testNullModel.cc
deleted file mode 100644
index 653393e11e38e380424e53d618befe7fab4b2ace..0000000000000000000000000000000000000000
--- a/Processes/NullModel/testNullModel.cc
+++ /dev/null
@@ -1,53 +0,0 @@
-/*
- * (c) Copyright 2018 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.
- */
-
-#include <catch2/catch.hpp>
-
-#include <corsika/process/null_model/NullModel.h>
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-using namespace corsika::units::si;
-using namespace corsika::process::null_model;
-using namespace corsika;
-
-TEST_CASE("NullModel", "[processes]") {
-
- auto const& dummyCS =
- geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- geometry::Point const origin(dummyCS, {0_m, 0_m, 0_m});
- geometry::Vector<units::si::SpeedType::dimension_type> v(dummyCS, 0_m / second,
- 0_m / second, 1_m / second);
- geometry::Line line(origin, v);
- geometry::Trajectory<geometry::Line> track(line, 10_s);
-
- setup::Stack stack;
- setup::Stack::ParticleType particle = stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::Electron, 100_GeV,
- corsika::stack::MomentumVector(dummyCS, {0_GeV, 0_GeV, -1_GeV}),
- geometry::Point(dummyCS, {0_m, 0_m, 10_km}), 0_ns});
- SECTION("interface") {
-
- NullModel model(10_m);
-
- [[maybe_unused]] const process::EProcessReturn ret =
- model.DoContinuous(particle, track);
- LengthType const length = model.MaxStepLength(particle, track);
-
- CHECK((length / 10_m) == Approx(1));
- }
-}
diff --git a/Processes/ObservationPlane/ObservationPlane.cc b/Processes/ObservationPlane/ObservationPlane.cc
index cb7ee95bd6241c942189417da34107a6e19efbd0..ab3a8023215528084dd242b656597f74e9989d23 100644
--- a/Processes/ObservationPlane/ObservationPlane.cc
+++ b/Processes/ObservationPlane/ObservationPlane.cc
@@ -7,9 +7,9 @@
*/
#include <corsika/process/observation_plane/ObservationPlane.h>
+#include <corsika/logging/Logging.h>
#include <fstream>
-#include <iostream>
using namespace corsika::process::observation_plane;
using namespace corsika::units::si;
@@ -25,7 +25,7 @@ ObservationPlane::ObservationPlane(geometry::Plane const& obsPlane,
}
corsika::process::EProcessReturn ObservationPlane::DoContinuous(
- setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
+ setup::Stack::ParticleType& particle, setup::Trajectory const& trajectory) {
TimeType const timeOfIntersection =
(plane_.GetCenter() - trajectory.GetR0()).dot(plane_.GetNormal()) /
trajectory.GetV0().dot(plane_.GetNormal());
@@ -45,6 +45,7 @@ corsika::process::EProcessReturn ObservationPlane::DoContinuous(
if (deleteOnHit_) {
count_ground_++;
energy_ground_ += energy;
+ particle.Delete();
return process::EProcessReturn::eParticleAbsorbed;
} else {
return process::EProcessReturn::eOk;
@@ -62,15 +63,19 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
}
auto const pointOfIntersection = trajectory.GetPosition(timeOfIntersection);
- return (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
+ auto dist = (trajectory.GetR0() - pointOfIntersection).norm() * 1.0001;
+ C8LOG_TRACE("ObservationPlane::MaxStepLength l={} m", dist / 1_m);
+ return dist;
}
void ObservationPlane::ShowResults() const {
- std::cout << " ******************************" << std::endl
- << " ObservationPlane: " << std::endl;
- std::cout << " energy in ground (GeV) : " << energy_ground_ / 1_GeV << std::endl
- << " no. of particles in ground : " << count_ground_ << std::endl
- << " ******************************" << std::endl;
+ C8LOG_INFO(
+ " ******************************\n"
+ " ObservationPlane: \n"
+ " energy in ground (GeV) : {}\n"
+ " no. of particles in ground : {}\n"
+ " ******************************",
+ energy_ground_ / 1_GeV, count_ground_);
}
void ObservationPlane::Reset() {
diff --git a/Processes/ObservationPlane/ObservationPlane.h b/Processes/ObservationPlane/ObservationPlane.h
index 7223defb1748b5c67a04969316725fd09c65b8d9..86c5ba0f2100043f5a7c4b9c2ba250b1e3a57a20 100644
--- a/Processes/ObservationPlane/ObservationPlane.h
+++ b/Processes/ObservationPlane/ObservationPlane.h
@@ -29,7 +29,7 @@ namespace corsika::process::observation_plane {
ObservationPlane(geometry::Plane const&, std::string const&, bool = true);
corsika::process::EProcessReturn DoContinuous(
- corsika::setup::Stack::ParticleType const& vParticle,
+ corsika::setup::Stack::ParticleType& vParticle,
corsika::setup::Trajectory const& vTrajectory);
corsika::units::si::LengthType MaxStepLength(
diff --git a/Processes/ParticleCut/ParticleCut.cc b/Processes/ParticleCut/ParticleCut.cc
index 72a97ee065d654799c7ae72cb4742a0d8d0a1452..728c7d25bb149ea20fb492628c9326f489426e62 100644
--- a/Processes/ParticleCut/ParticleCut.cc
+++ b/Processes/ParticleCut/ParticleCut.cc
@@ -27,9 +27,9 @@ namespace corsika::process {
if (vP.GetPID() == particles::Code::Nucleus) {
// calculate energy per nucleon
auto const ElabNuc = energyLab / vP.GetNuclearA();
- return (ElabNuc <= fECut);
+ return (ElabNuc <= energy_cut_);
} else {
- return (energyLab <= fECut);
+ return (energyLab <= energy_cut_);
}
}
@@ -87,13 +87,12 @@ namespace corsika::process {
return false; // this particle will not be removed/cut
}
- EProcessReturn ParticleCut::DoSecondaries(corsika::setup::StackView& vS) {
+ void ParticleCut::DoSecondaries(corsika::setup::StackView& vS) {
auto particle = vS.begin();
while (particle != vS.end()) {
if (checkCutParticle(particle)) { particle.Delete(); }
++particle; // next entry in SecondaryView
}
- return EProcessReturn::eOk;
}
process::EProcessReturn ParticleCut::DoContinuous(
@@ -102,13 +101,15 @@ namespace corsika::process {
C8LOG_TRACE("ParticleCut::DoContinuous");
if (checkCutParticle(particle)) {
C8LOG_TRACE("removing during continuous");
+ particle.Delete();
+ // signal to upstream code that this particle was deleted
return process::EProcessReturn::eParticleAbsorbed;
}
return process::EProcessReturn::eOk;
}
ParticleCut::ParticleCut(const units::si::HEPEnergyType eCut, bool em, bool inv)
- : fECut(eCut)
+ : energy_cut_(eCut)
, bCutEm(em)
, bCutInv(inv) {
fEmEnergy = 0_GeV;
diff --git a/Processes/ParticleCut/ParticleCut.h b/Processes/ParticleCut/ParticleCut.h
index 4705051d4261028c49caca1af39066563fd1bff6..b0e8977530a8fd7a800908bd6eb48519a411d9f0 100644
--- a/Processes/ParticleCut/ParticleCut.h
+++ b/Processes/ParticleCut/ParticleCut.h
@@ -12,6 +12,7 @@
#include <corsika/process/ContinuousProcess.h>
#include <corsika/process/SecondariesProcess.h>
#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
@@ -19,7 +20,7 @@ namespace corsika::process {
class ParticleCut : public process::SecondariesProcess<ParticleCut>,
public corsika::process::ContinuousProcess<ParticleCut> {
- units::si::HEPEnergyType const fECut;
+ units::si::HEPEnergyType const energy_cut_;
bool bCutEm;
bool bCutInv;
@@ -32,7 +33,7 @@ namespace corsika::process {
public:
ParticleCut(const units::si::HEPEnergyType eCut, bool em, bool inv);
- EProcessReturn DoSecondaries(corsika::setup::StackView&);
+ void DoSecondaries(corsika::setup::StackView&);
EProcessReturn DoContinuous(corsika::setup::Stack::ParticleType& vParticle,
corsika::setup::Trajectory const& vTrajectory);
@@ -42,7 +43,7 @@ namespace corsika::process {
return units::si::meter * std::numeric_limits<double>::infinity();
}
- units::si::HEPEnergyType GetECut() const { return fECut; }
+ units::si::HEPEnergyType GetECut() const { return energy_cut_; }
units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
units::si::HEPEnergyType GetCutEnergy() const { return fEnergy; }
units::si::HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
diff --git a/Processes/ParticleCut/testParticleCut.cc b/Processes/ParticleCut/testParticleCut.cc
index 2092636fa372628080769d22c52679289d9d045b..921eecc3ea2e59095f276bb03d7faec4ab43174a 100644
--- a/Processes/ParticleCut/testParticleCut.cc
+++ b/Processes/ParticleCut/testParticleCut.cc
@@ -37,15 +37,16 @@ TEST_CASE("ParticleCut", "[processes]") {
const HEPEnergyType Eabove = 1_TeV;
const HEPEnergyType Ebelow = 10_GeV;
// list of arbitrary particles
- std::vector<particles::Code> particleList = {
+ const std::vector<particles::Code> particleList = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short,
particles::Code::Electron, particles::Code::MuPlus, particles::Code::NuE,
- particles::Code::Neutron};
+ particles::Code::Neutron, particles::Code::NuMu};
SECTION("cut on particle type: inv") {
ParticleCut cut(20_GeV, false, true);
+ CHECK(cut.GetECut() == 20_GeV);
// add primary particle to stack
auto particle = stack.AddParticle(
@@ -62,15 +63,15 @@ TEST_CASE("ParticleCut", "[processes]") {
// add secondaries, all with energies above the threshold
// only cut is by species
for (auto proType : particleList)
- projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType>{
+ projectile.AddSecondary(std::make_tuple(
proType, Eabove, corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
cut.DoSecondaries(view);
CHECK(view.getEntries() == 9);
+ CHECK(cut.GetNumberInvParticles() == 2);
+ CHECK(cut.GetInvEnergy() / 1_GeV == 2000);
}
SECTION("cut on particle type: em") {
@@ -79,11 +80,9 @@ TEST_CASE("ParticleCut", "[processes]") {
// add primary particle to stack
auto particle = stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::Proton, Eabove,
- corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+ std::make_tuple(particles::Code::Proton, Eabove,
+ corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
// view on secondary particles
corsika::setup::StackView view(particle);
// ref. to primary particle through the secondary view.
@@ -91,16 +90,16 @@ TEST_CASE("ParticleCut", "[processes]") {
auto projectile = view.GetProjectile();
// add secondaries, all with energies above the threshold
// only cut is by species
- for (auto proType : particleList)
- projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType>{
+ for (auto proType : particleList) {
+ projectile.AddSecondary(std::make_tuple(
proType, Eabove, corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
-
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
+ }
cut.DoSecondaries(view);
- CHECK(view.getEntries() == 9);
+ CHECK(view.getEntries() == 10);
+ CHECK(cut.GetNumberEmParticles() == 1);
+ CHECK(cut.GetEmEnergy() / 1_GeV == 1000);
}
SECTION("cut low energy") {
@@ -108,11 +107,9 @@ TEST_CASE("ParticleCut", "[processes]") {
// add primary particle to stack
auto particle = stack.AddParticle(
- std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- particles::Code::Proton, Eabove,
- corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+ std::make_tuple(particles::Code::Proton, Eabove,
+ corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
// view on secondary particles
setup::StackView view{particle};
// ref. to primary particle through the secondary view.
@@ -121,15 +118,52 @@ TEST_CASE("ParticleCut", "[processes]") {
// add secondaries, all with energies below the threshold
// only cut is by species
for (auto proType : particleList)
- projectile.AddSecondary(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType>{
+ projectile.AddSecondary(std::make_tuple(
proType, Ebelow, corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
- geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns));
+ unsigned short A = 18;
+ unsigned short Z = 8;
+ projectile.AddSecondary(
+ std::make_tuple(particles::Code::Nucleus, Eabove * A,
+ corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns, A, Z));
+ projectile.AddSecondary(
+ std::make_tuple(particles::Code::Nucleus, Ebelow * A,
+ corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns, A, Z));
+
+ cut.DoSecondaries(view);
+
+ CHECK(view.getEntries() == 1);
+ CHECK(view.getSize() == 13);
+ }
+
+ SECTION("cut on time") {
+ ParticleCut cut(20_GeV, false, false);
+ const TimeType too_late = 1_s;
+ // add primary particle to stack
+ auto particle = stack.AddParticle(
+ std::make_tuple(particles::Code::Proton, Eabove,
+ corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), 1_ns));
+ // view on secondary particles
+ corsika::setup::StackView view(particle);
+ // ref. to primary particle through the secondary view.
+ // only this way the secondary view is populated
+ auto projectile = view.GetProjectile();
+ // add secondaries, all with energies above the threshold
+ // only cut is by species
+ for (auto proType : particleList) {
+ projectile.AddSecondary(std::make_tuple(
+ proType, Eabove, corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 0_GeV}),
+ geometry::Point(rootCS, 0_m, 0_m, 0_m), too_late));
+ }
cut.DoSecondaries(view);
- REQUIRE(view.getEntries() == 0);
- REQUIRE(view.getSize() == 10);
+ CHECK(view.getEntries() == 0);
+ CHECK(cut.GetCutEnergy() / 1_GeV == 11000);
+ cut.Reset();
+ CHECK(cut.GetCutEnergy() == 0_GeV);
}
}
diff --git a/Processes/Proposal/CMakeLists.txt b/Processes/Proposal/CMakeLists.txt
index cc0be5438fc75bb48f32b0501957d41cf1ce38a2..08a37a68b441536a1fba53314165b00991d64d19 100644
--- a/Processes/Proposal/CMakeLists.txt
+++ b/Processes/Proposal/CMakeLists.txt
@@ -22,6 +22,7 @@ TARGET_LINK_LIBRARIES (
CORSIKAunits
CORSIKAthirdparty
CORSIKAgeometry
+ CORSIKAlogging
CORSIKAenvironment
ProcessParticleCut
PROPOSAL::PROPOSAL
diff --git a/Processes/Proposal/ContinuousProcess.cc b/Processes/Proposal/ContinuousProcess.cc
index 59ec7e87f56d92cc4433abf090d838ccd24fc6c7..371f28262bd3b7a5e5f67fc2d2415c221cd43ffc 100644
--- a/Processes/Proposal/ContinuousProcess.cc
+++ b/Processes/Proposal/ContinuousProcess.cc
@@ -15,6 +15,7 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
+#include <corsika/logging/Logging.h>
#include <iostream>
@@ -106,19 +107,6 @@ namespace corsika::process::proposal {
// if the particle has a charge take multiple scattering into account
if (vP.GetChargeNumber() != 0) Scatter(vP, dE, dX);
-
- // Update the energy and absorbe the particle if it's below the energy
- // threshold, because it will no longer propagated.
- if (final_energy <= emCut_) {
- vP.SetEnergy(emCut_);
- vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm());
- return process::EProcessReturn::eParticleAbsorbed;
- }
- if (final_energy <= emCut_) {
- vP.SetEnergy(emCut_);
- vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm());
- return process::EProcessReturn::eParticleAbsorbed;
- }
vP.SetEnergy(final_energy);
vP.SetMomentum(vP.GetMomentum() * vP.GetEnergy() / vP.GetMomentum().GetNorm());
return process::EProcessReturn::eOk;
@@ -132,7 +120,7 @@ namespace corsika::process::proposal {
if (!CanInteract(vP.GetPID()))
return units::si::meter * std::numeric_limits<double>::infinity();
- // Limit the step size of a conitnous loss. The maximal continuous loss seems to be a
+ // Limit the step size of a conitnuous loss. The maximal continuous loss seems to be a
// hyper parameter which must be adjusted.
//
// in any case: never go below 0.99*emCut_ This needs to be
@@ -152,7 +140,10 @@ namespace corsika::process::proposal {
1_g / square(1_cm);
// return it in distance aequivalent
- return vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage);
+ auto dist = vP.GetNode()->GetModelProperties().ArclengthFromGrammage(vT, grammage);
+ C8LOG_TRACE("PROPOSAL::MaxStepLength X={} g/cm2, l={} m ",
+ grammage / 1_g * square(1_cm), dist / 1_m);
+ return dist;
}
void ContinuousProcess::ShowResults() const {
diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc
index 4caf02487afa2d7a89d600a33a2c46eefb00c294..92920dfbacdcaa809df665db59ae57bc70d8b46c 100644
--- a/Processes/Pythia/Decay.cc
+++ b/Processes/Pythia/Decay.cc
@@ -28,7 +28,8 @@ using Track = Trajectory;
namespace corsika::process::pythia {
- Decay::Decay() {
+ Decay::Decay(const bool print_listing)
+ : print_listing_(print_listing) {
// set random number generator in pythia
Pythia8::RndmEngine* rndm = new corsika::process::pythia::Random();
@@ -231,8 +232,10 @@ namespace corsika::process::pythia {
else
cout << "Pythia::Decay: particles after decay: " << event.size() << endl;
- // list final state
- event.list();
+ if (print_listing_) {
+ // list final state
+ event.list();
+ }
// loop over final state
for (int i = 0; i < event.size(); ++i)
diff --git a/Processes/Pythia/Decay.h b/Processes/Pythia/Decay.h
index fb4708df15e3c7948a0d352d2ebcd326af68aa53..51a43b7e3824ec399eefa4bf96893144cb2db8f4 100644
--- a/Processes/Pythia/Decay.h
+++ b/Processes/Pythia/Decay.h
@@ -11,6 +11,8 @@
#include <Pythia8/Pythia.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/DecayProcess.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/geometry/FourVector.h>
namespace corsika::process {
@@ -21,9 +23,10 @@ namespace corsika::process {
class Decay : public corsika::process::DecayProcess<Decay> {
int fCount = 0;
bool handleAllDecays_ = true;
+ bool print_listing_ = false;
public:
- Decay();
+ Decay(const bool print_listing = false);
Decay(std::set<particles::Code>);
~Decay();
diff --git a/Processes/Pythia/Interaction.cc b/Processes/Pythia/Interaction.cc
index fe641a62629a5bfc2b75f4ce9d0ada6c984322ce..7a0dbfa64878c77d3fbed0b6fd58798219dc0b38 100644
--- a/Processes/Pythia/Interaction.cc
+++ b/Processes/Pythia/Interaction.cc
@@ -29,9 +29,9 @@ namespace corsika::process::pythia {
typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
- Interaction::~Interaction() {}
+ Interaction::Interaction(const bool print_listing)
+ : print_listing_(print_listing) {
- Interaction::Interaction() {
cout << "Pythia::Interaction n=" << fCount << endl;
using random::RNGManager;
@@ -164,7 +164,7 @@ namespace corsika::process::pythia {
}
template <>
- units::si::GrammageType Interaction::GetInteractionLength(Particle& p) {
+ units::si::GrammageType Interaction::GetInteractionLength(const Particle& p) {
using namespace units;
using namespace units::si;
@@ -365,8 +365,11 @@ namespace corsika::process::pythia {
// link to pythia stack
Pythia8::Event& event = fPythia.event;
- // print final state
- event.list();
+
+ if (print_listing_) {
+ // list final state
+ event.list();
+ }
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV;
diff --git a/Processes/Pythia/Interaction.h b/Processes/Pythia/Interaction.h
index 4b06aeb6d58f6bf034b49b054376e540fa279628..ef2bdfbcc5566951434c73df86ec621e561f48b9 100644
--- a/Processes/Pythia/Interaction.h
+++ b/Processes/Pythia/Interaction.h
@@ -22,10 +22,11 @@ namespace corsika::process::pythia {
int fCount = 0;
bool fInitialized = false;
+ bool print_listing_ = false;
public:
- Interaction();
- ~Interaction();
+ Interaction(const bool print_listing = false);
+ ~Interaction() = default;
void SetParticleListStable(std::vector<particles::Code> const&);
void SetUnstable(const corsika::particles::Code);
@@ -47,7 +48,7 @@ namespace corsika::process::pythia {
const corsika::units::si::HEPEnergyType CoMenergy);
template <typename TParticle>
- corsika::units::si::GrammageType GetInteractionLength(TParticle&);
+ corsika::units::si::GrammageType GetInteractionLength(const TParticle&);
/**
In this function PYTHIA is called to produce one event. The
diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
index 62d0b45819be56c8fcbec7b67360a17c96650b73..8297739b793e33d7b4fe9b65bfd66def992f372e 100644
--- a/Processes/Sibyll/Decay.cc
+++ b/Processes/Sibyll/Decay.cc
@@ -23,11 +23,12 @@ using namespace corsika::setup;
using SetupView = corsika::setup::StackView;
using SetupProjectile = corsika::setup::StackView::ParticleType;
-using SetupParticle = corsika::setup::Stack::ParticleType;
+using Particle = corsika::setup::Stack::ParticleType;
namespace corsika::process::sibyll {
- Decay::Decay() {
+ Decay::Decay(const bool sibyll_printout_on)
+ : sibyll_listing_(sibyll_printout_on) {
// switch off decays to avoid internal decay chains
SetAllStable();
// handle all decays by default
@@ -40,7 +41,7 @@ namespace corsika::process::sibyll {
SetAllStable();
}
- Decay::~Decay() { C8LOG_DEBUG(fmt::format("Sibyll::Decay n={}", fCount)); }
+ Decay::~Decay() { C8LOG_DEBUG("Sibyll::Decay n={}", count_); }
bool Decay::CanHandleDecay(const particles::Code vParticleCode) {
using namespace corsika::particles;
@@ -60,7 +61,7 @@ namespace corsika::process::sibyll {
void Decay::SetHandleDecay(const particles::Code vParticleCode) {
handleAllDecays_ = false;
- C8LOG_DEBUG(fmt::format("Sibyll::Decay: set to handle decay of {}", vParticleCode));
+ C8LOG_DEBUG("Sibyll::Decay: set to handle decay of {}", vParticleCode);
if (Decay::CanHandleDecay(vParticleCode))
handledDecays_.insert(vParticleCode);
else
@@ -102,13 +103,13 @@ namespace corsika::process::sibyll {
}
void Decay::SetUnstable(const particles::Code vCode) {
- C8LOG_DEBUG(fmt::format("Sibyll::Decay: setting {} unstable. ", vCode));
+ C8LOG_DEBUG("Sibyll::Decay: setting {} unstable. ", vCode);
const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
}
void Decay::SetStable(const particles::Code vCode) {
- C8LOG_DEBUG(fmt::format("Sibyll::Decay: setting {} stable. ", vCode));
+ C8LOG_DEBUG("Sibyll::Decay: setting {} stable. ", vCode);
const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
}
@@ -121,28 +122,26 @@ 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);
- C8LOG_DEBUG(
- fmt::format("Decay: Sibyll decay configuration: {} is {}", vCode,
- (s_csydec_.idb[absSibCode - 1] <= 0) ? "stable" : "unstable"));
+ 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");
}
void Decay::PrintDecayConfig() {
- C8LOG_DEBUG(fmt::format("Sibyll::Decay: decay configuration:"));
+ C8LOG_DEBUG("Sibyll::Decay: decay configuration:");
if (handleAllDecays_) {
- C8LOG_DEBUG(fmt::format(
- " all particles known to Sibyll are handled by Sibyll::Decay!"));
+ C8LOG_DEBUG(" all particles known to Sibyll are handled by Sibyll::Decay!");
} else {
- for (auto& pCode : handledDecays_) {
- C8LOG_DEBUG(fmt::format(" Decay of {} is handled by Sibyll!", pCode));
+ for ([[maybe_unused]] auto& pCode : handledDecays_) {
+ C8LOG_DEBUG(" Decay of {} is handled by Sibyll!", pCode);
}
}
}
template <>
- units::si::TimeType Decay::GetLifetime(SetupParticle const& vP) {
+ units::si::TimeType Decay::GetLifetime(Particle const& vP) {
using namespace units::si;
const particles::Code pid = vP.GetPID();
@@ -155,19 +154,19 @@ 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(fmt::format("Sibyll::Decay: code: {} ", vP.GetPID()));
- C8LOG_DEBUG(fmt::format("Sibyll::Decay: MinStep: t0: {} ", t0));
- C8LOG_DEBUG(fmt::format("Sibyll::Decay: MinStep: energy: {} GeV ", E / 1_GeV));
- C8LOG_DEBUG(fmt::format("Sibyll::Decay: momentum: {} GeV ",
- vP.GetMomentum().GetComponents() / 1_GeV));
- C8LOG_DEBUG(fmt::format("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());
- C8LOG_DEBUG(fmt::format("Sibyll::Decay: sib mass: {}", get_sibyll_mass2(sib_id)));
- C8LOG_DEBUG(fmt::format("Sibyll::Decay: MinStep: gamma: {}", gamma));
- C8LOG_DEBUG(fmt::format("Sibyll::Decay: MinStep: tau {} s: ", lifetime / 1_s));
+ C8LOG_DEBUG("Sibyll::Decay: code: {} ", vP.GetPID());
+ C8LOG_DEBUG("Sibyll::Decay: MinStep: t0: {} ", t0);
+ C8LOG_DEBUG("Sibyll::Decay: MinStep: energy: {} GeV ", E / 1_GeV);
+ C8LOG_DEBUG("Sibyll::Decay: momentum: {} GeV ",
+ 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);
+ [[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);
return lifetime;
} else
@@ -186,7 +185,7 @@ namespace corsika::process::sibyll {
if (!IsDecayHandled(pCode))
throw std::runtime_error("STOP! Sibyll not configured to execute this decay!");
- fCount++;
+ count_++;
SibStack ss;
ss.Clear();
@@ -206,12 +205,14 @@ namespace corsika::process::sibyll {
PrintDecayConfig(pCode);
// call sibyll decay
- C8LOG_DEBUG(fmt::format("Decay: calling Sibyll decay routine.."));
+ C8LOG_DEBUG("Decay: calling Sibyll decay routine..");
decsib_();
- // print output
- int print_unit = 6;
- sib_list_(print_unit);
+ if (sibyll_listing_) {
+ // print output
+ int print_unit = 6;
+ sib_list_(print_unit);
+ }
// reset to stable
SetStable(pCode);
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 42867ed2b560487901489c0df1bb153380948147..e5862344e86c13d6dbc11ee78246987fe3fcf198 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -20,11 +20,12 @@ namespace corsika::process {
namespace sibyll {
class Decay : public corsika::process::DecayProcess<Decay> {
- int fCount = 0;
+ int count_ = 0;
bool handleAllDecays_ = true;
+ bool sibyll_listing_ = false;
public:
- Decay();
+ Decay(const bool sibyll_listing = false);
Decay(std::set<particles::Code>);
~Decay();
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index b028ca8769e3761c8adf6453b6162478271a756e..aa222b7f65ced91a80b5e6402a0ae5a45926c6dc 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -25,7 +25,7 @@ using std::tuple;
using namespace corsika;
using namespace corsika::setup;
-using SetupParticle = setup::Stack::StackIterator;
+using Particle = setup::Stack::StackIterator;
using SetupView = setup::StackView;
using Track = Trajectory;
@@ -33,7 +33,8 @@ namespace corsika::process::sibyll {
bool Interaction::initialized_ = false;
- Interaction::Interaction() {
+ Interaction::Interaction(const bool sibyll_printout_on)
+ : sibyll_listing_(sibyll_printout_on) {
using random::RNGManager;
// initialize Sibyll
@@ -81,8 +82,7 @@ namespace corsika::process::sibyll {
}
template <>
- units::si::GrammageType Interaction::GetInteractionLength(
- SetupParticle const& vP) const {
+ units::si::GrammageType Interaction::GetInteractionLength(Particle const& vP) const {
using namespace units;
using namespace units::si;
@@ -172,187 +172,186 @@ namespace corsika::process::sibyll {
auto const projectile = view.GetProjectile();
const auto corsikaBeamId = projectile.GetPID();
- C8LOG_DEBUG(
- fmt::format("ProcessSibyll: "
- "DoInteraction: {} interaction? ",
- corsikaBeamId, process::sibyll::CanInteract(corsikaBeamId)));
if (particles::IsNucleus(corsikaBeamId)) {
// nuclei handled by different process, this should not happen
throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!");
}
- if (process::sibyll::CanInteract(corsikaBeamId)) {
- // position and time of interaction, not used in Sibyll
- Point const pOrig = projectile.GetPosition();
- TimeType const tOrig = projectile.GetTime();
+ // position and time of interaction, not used in Sibyll
+ Point const pOrig = projectile.GetPosition();
+ TimeType const tOrig = projectile.GetTime();
- // define projectile
- HEPEnergyType const eProjectileLab = projectile.GetEnergy();
- auto const pProjectileLab = projectile.GetMomentum();
- const CoordinateSystem& originalCS = pProjectileLab.GetCoordinateSystem();
+ // define projectile
+ HEPEnergyType const eProjectileLab = projectile.GetEnergy();
+ auto const pProjectileLab = projectile.GetMomentum();
+ const CoordinateSystem& originalCS = pProjectileLab.GetCoordinateSystem();
- // define target
- // for Sibyll is always a single nucleon
- // FOR NOW: target is always at rest
- const auto eTargetLab = 0_GeV + constants::nucleonMass;
- const auto pTargetLab = MomentumVector(originalCS, 0_GeV, 0_GeV, 0_GeV);
- const FourVector PtargLab(eTargetLab, pTargetLab);
+ C8LOG_DEBUG(
+ "ProcessSibyll: "
+ "DoInteraction: pid {} interaction ",
+ corsikaBeamId);
- C8LOG_DEBUG(
- fmt::format("Interaction: ebeam lab: {} GeV"
- "Interaction: pbeam lab: {} GeV",
- eProjectileLab / 1_GeV, pProjectileLab.GetComponents()));
- C8LOG_DEBUG(
- fmt::format("Interaction: etarget lab: {} GeV "
- "Interaction: ptarget lab: {} GeV",
- eTargetLab / 1_GeV, pTargetLab.GetComponents() / 1_GeV));
+ // define target
+ // for Sibyll is always a single nucleon
+ // FOR NOW: target is always at rest
+ const auto eTargetLab = 0_GeV + constants::nucleonMass;
+ const auto pTargetLab = MomentumVector(originalCS, 0_GeV, 0_GeV, 0_GeV);
+ const FourVector PtargLab(eTargetLab, pTargetLab);
- const FourVector PprojLab(eProjectileLab, pProjectileLab);
+ C8LOG_DEBUG(
+ "Interaction: ebeam lab: {} GeV"
+ "Interaction: pbeam lab: {} GeV",
+ eProjectileLab / 1_GeV, pProjectileLab.GetComponents());
+ C8LOG_DEBUG(
+ "Interaction: etarget lab: {} GeV "
+ "Interaction: ptarget lab: {} GeV",
+ eTargetLab / 1_GeV, pTargetLab.GetComponents() / 1_GeV);
- // define target kinematics in lab frame
- // define boost to and from CoM frame
- // CoM frame definition in Sibyll projectile: +z
- COMBoost const boost(PprojLab, constants::nucleonMass);
- auto const& csPrime = boost.GetRotatedCS();
+ const FourVector PprojLab(eProjectileLab, pProjectileLab);
- // just for show:
- // boost projecticle
- auto const PprojCoM = boost.toCoM(PprojLab);
+ // define target kinematics in lab frame
+ // define boost to and from CoM frame
+ // CoM frame definition in Sibyll projectile: +z
+ COMBoost const boost(PprojLab, constants::nucleonMass);
+ auto const& csPrime = boost.GetRotatedCS();
- // boost target
- auto const PtargCoM = boost.toCoM(PtargLab);
+ // just for show:
+ // boost projecticle
+ [[maybe_unused]] auto const PprojCoM = boost.toCoM(PprojLab);
- C8LOG_DEBUG(
- fmt::format("Interaction: ebeam CoM: {} GeV "
- "Interaction: pbeam CoM: {} GeV ",
- PprojCoM.GetTimeLikeComponent() / 1_GeV,
- PprojCoM.GetSpaceLikeComponents().GetComponents(csPrime) / 1_GeV));
- C8LOG_DEBUG(
- fmt::format("Interaction: etarget CoM: {} GeV "
- "Interaction: ptarget CoM: {} GeV ",
- PtargCoM.GetTimeLikeComponent() / 1_GeV,
- PtargCoM.GetSpaceLikeComponents().GetComponents(csPrime) / 1_GeV));
-
- C8LOG_DEBUG(fmt::format("Interaction: position of interaction: {} ",
- pOrig.GetCoordinates()));
- C8LOG_DEBUG(fmt::format("Interaction: time: {} ", tOrig));
-
- HEPEnergyType Etot = eProjectileLab + eTargetLab;
- MomentumVector Ptot = projectile.GetMomentum();
- // invariant mass, i.e. cm. energy
- HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
-
- // sample target mass number
- auto const* currentNode = projectile.GetNode();
- auto const& mediumComposition =
- currentNode->GetModelProperties().GetNuclearComposition();
- // get cross sections for target materials
- /*
- Here we read the cross section from the interaction model again,
- should be passed from GetInteractionLength if possible
- */
- //#warning reading interaction cross section again, should not be necessary
- auto const& compVec = mediumComposition.GetComponents();
- std::vector<CrossSectionType> cross_section_of_components(compVec.size());
-
- for (size_t i = 0; i < compVec.size(); ++i) {
- auto const targetId = compVec[i];
- const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm);
- [[maybe_unused]] const auto& dummy_sigEla = sigEla;
- cross_section_of_components[i] = sigProd;
- }
+ // boost target
+ [[maybe_unused]] auto const PtargCoM = boost.toCoM(PtargLab);
- const auto targetCode =
- mediumComposition.SampleTarget(cross_section_of_components, RNG_);
- C8LOG_DEBUG(fmt::format("Interaction: target selected: {} ", targetCode));
- /*
- FOR NOW: allow nuclei with A<18 or protons only.
- when medium composition becomes more complex, approximations will have to be
- allowed air in atmosphere also contains some Argon.
- */
- int targetSibCode = -1;
- if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode);
- if (targetCode == particles::Proton::GetCode()) targetSibCode = 1;
- C8LOG_DEBUG(fmt::format("Interaction: sibyll code: {}", targetSibCode));
- if (targetSibCode > maxTargetMassNumber_ || targetSibCode < 1)
- throw std::runtime_error(
- "Sibyll target outside range. Only nuclei with A<18 or protons are "
- "allowed.");
+ C8LOG_DEBUG(
+ "Interaction: ebeam CoM: {} GeV "
+ "Interaction: pbeam CoM: {} GeV ",
+ PprojCoM.GetTimeLikeComponent() / 1_GeV,
+ PprojCoM.GetSpaceLikeComponents().GetComponents(csPrime) / 1_GeV);
+ C8LOG_DEBUG(
+ "Interaction: etarget CoM: {} GeV "
+ "Interaction: ptarget CoM: {} GeV ",
+ PtargCoM.GetTimeLikeComponent() / 1_GeV,
+ PtargCoM.GetSpaceLikeComponents().GetComponents(csPrime) / 1_GeV);
+
+ C8LOG_DEBUG("Interaction: position of interaction: {} ", pOrig.GetCoordinates());
+ C8LOG_DEBUG("Interaction: time: {} ", tOrig);
+
+ HEPEnergyType Etot = eProjectileLab + eTargetLab;
+ MomentumVector Ptot = projectile.GetMomentum();
+ // invariant mass, i.e. cm. energy
+ HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
+
+ // sample target mass number
+ auto const* currentNode = projectile.GetNode();
+ auto const& mediumComposition =
+ currentNode->GetModelProperties().GetNuclearComposition();
+ // get cross sections for target materials
+ /*
+ Here we read the cross section from the interaction model again,
+ should be passed from GetInteractionLength if possible
+ */
+ //#warning reading interaction cross section again, should not be necessary
+ auto const& compVec = mediumComposition.GetComponents();
+ std::vector<CrossSectionType> cross_section_of_components(compVec.size());
+
+ for (size_t i = 0; i < compVec.size(); ++i) {
+ auto const targetId = compVec[i];
+ const auto [sigProd, sigEla] = GetCrossSection(corsikaBeamId, targetId, Ecm);
+ [[maybe_unused]] const auto& dummy_sigEla = sigEla;
+ cross_section_of_components[i] = sigProd;
+ }
+
+ const auto targetCode =
+ mediumComposition.SampleTarget(cross_section_of_components, RNG_);
+ C8LOG_DEBUG("Interaction: target selected: {} ", targetCode);
+ /*
+ FOR NOW: allow nuclei with A<18 or protons only.
+ when medium composition becomes more complex, approximations will have to be
+ allowed air in atmosphere also contains some Argon.
+ */
+ int targetSibCode = -1;
+ if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode);
+ if (targetCode == particles::Proton::GetCode()) targetSibCode = 1;
+ C8LOG_DEBUG("Interaction: sibyll code: {}", targetSibCode);
+ if (targetSibCode > maxTargetMassNumber_ || targetSibCode < 1)
+ throw std::runtime_error(
+ "Sibyll target outside range. Only nuclei with A<18 or protons are "
+ "allowed.");
- // beam id for sibyll
- const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId);
+ // beam id for sibyll
+ const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId);
+ C8LOG_DEBUG(
+ "Interaction: "
+ " DoInteraction: E(GeV): {} "
+ " Ecm(GeV): {} ",
+ eProjectileLab / 1_GeV, Ecm / 1_GeV);
+ if (Ecm > GetMaxEnergyCoM())
+ throw std::runtime_error("Interaction::DoInteraction: CoM energy too high!");
+ // FR: removed eProjectileLab < 8.5_GeV ||
+ if (Ecm < GetMinEnergyCoM()) {
C8LOG_DEBUG(
- fmt::format("Interaction: "
- " DoInteraction: E(GeV): {} "
- " Ecm(GeV): {} ",
- eProjectileLab / 1_GeV, Ecm / 1_GeV));
- if (Ecm > GetMaxEnergyCoM())
- throw std::runtime_error("Interaction::DoInteraction: CoM energy too high!");
- // FR: removed eProjectileLab < 8.5_GeV ||
- if (Ecm < GetMinEnergyCoM()) {
- C8LOG_DEBUG(
- fmt::format("Interaction: "
- " DoInteraction: should have dropped particle.. "
- "THIS IS AN ERROR"));
- throw std::runtime_error("energy too low for SIBYLL");
- } else {
- count_++;
- // Sibyll does not know about units..
- const double sqs = Ecm / 1_GeV;
- // running sibyll, filling stack
- sibyll_(kBeam, targetSibCode, sqs);
+ "Interaction: "
+ " DoInteraction: should have dropped particle.. "
+ "THIS IS AN ERROR");
+ throw std::runtime_error("energy too low for SIBYLL");
+ } else {
+ count_++;
+ // Sibyll does not know about units..
+ const double sqs = Ecm / 1_GeV;
+ // running sibyll, filling stack
+ sibyll_(kBeam, targetSibCode, sqs);
+ if (sibyll_listing_) {
// print final state
int print_unit = 6;
sib_list_(print_unit);
nucCount_ += get_nwounded() - 1;
+ }
- // add particles from sibyll to stack
- // link to sibyll stack
- SibStack ss;
-
- MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
- HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
- for (auto& psib : ss) {
-
- // abort on particles that have decayed in Sibyll. Should not happen!
- if (psib.HasDecayed())
- throw std::runtime_error("found particle that decayed in SIBYLL!");
-
- // transform 4-momentum to lab. frame
- // note that the momentum needs to be rotated back
- auto const tmp = psib.GetMomentum().GetComponents();
- auto const pCoM = Vector<hepmomentum_d>(csPrime, tmp);
- HEPEnergyType const eCoM = psib.GetEnergy();
- auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
- auto const p3lab = Plab.GetSpaceLikeComponents();
- assert(p3lab.GetCoordinateSystem() == originalCS); // just to be sure!
-
- // add to corsika stack
- auto pnew = view.AddSecondary(
- make_tuple(process::sibyll::ConvertFromSibyll(psib.GetPID()),
- Plab.GetTimeLikeComponent(), p3lab, pOrig, tOrig));
-
- Plab_final += pnew.GetMomentum();
- Elab_final += pnew.GetEnergy();
- Ecm_final += psib.GetEnergy();
- }
- C8LOG_DEBUG(fmt::format(
- "conservation (all GeV):"
- "Ecm_initial(per nucleon)={}, Ecm_final(per nucleon)={}, "
- "Elab_initial={}, Elab_final={}, "
- "diff (%)={}, "
- "E in nucleons={}, "
- "Plab_initial={}, "
- "Plab_final={} ",
- Ecm / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Etot / 1_GeV,
- Elab_final / 1_GeV, (Elab_final / Etot / get_nwounded() - 1) * 100,
- constants::nucleonMass * get_nwounded() / 1_GeV,
- (pProjectileLab / 1_GeV).GetComponents(),
- (Plab_final / 1_GeV).GetComponents()));
+ // add particles from sibyll to stack
+ // link to sibyll stack
+ SibStack ss;
+
+ MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
+ for (auto& psib : ss) {
+
+ // abort on particles that have decayed in Sibyll. Should not happen!
+ if (psib.HasDecayed())
+ throw std::runtime_error("found particle that decayed in SIBYLL!");
+
+ // transform 4-momentum to lab. frame
+ // note that the momentum needs to be rotated back
+ auto const tmp = psib.GetMomentum().GetComponents();
+ auto const pCoM = Vector<hepmomentum_d>(csPrime, tmp);
+ HEPEnergyType const eCoM = psib.GetEnergy();
+ auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
+ auto const p3lab = Plab.GetSpaceLikeComponents();
+ assert(p3lab.GetCoordinateSystem() == originalCS); // just to be sure!
+
+ // add to corsika stack
+ auto pnew = view.AddSecondary(
+ make_tuple(process::sibyll::ConvertFromSibyll(psib.GetPID()),
+ Plab.GetTimeLikeComponent(), p3lab, pOrig, tOrig));
+
+ Plab_final += pnew.GetMomentum();
+ Elab_final += pnew.GetEnergy();
+ Ecm_final += psib.GetEnergy();
}
+ C8LOG_DEBUG(
+ "conservation (all GeV):"
+ "Ecm_initial(per nucleon)={}, Ecm_final(per nucleon)={}, "
+ "Elab_initial={}, Elab_final={}, "
+ "diff (%)={}, "
+ "E in nucleons={}, "
+ "Plab_initial={}, "
+ "Plab_final={} ",
+ Ecm / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV, Etot / 1_GeV,
+ Elab_final / 1_GeV, (Elab_final / Etot / get_nwounded() - 1) * 100,
+ constants::nucleonMass * get_nwounded() / 1_GeV,
+ (pProjectileLab / 1_GeV).GetComponents(), (Plab_final / 1_GeV).GetComponents());
}
return process::EProcessReturn::eOk;
}
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index cccc971bd36517f955f23da88d6b7a588fd4341d..52a43309d06efe3f0fbc622b84dc29a82c9d7328 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -21,9 +21,10 @@ namespace corsika::process::sibyll {
int count_ = 0;
int nucCount_ = 0;
static bool initialized_; ///! flag to assure init is done only once
+ bool sibyll_listing_;
public:
- Interaction();
+ Interaction(const bool sibyll_printout_on = false);
~Interaction();
void SetAllStable();
@@ -45,7 +46,7 @@ namespace corsika::process::sibyll {
const corsika::units::si::HEPEnergyType) const;
template <typename TParticle>
- corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const;
+ corsika::units::si::GrammageType GetInteractionLength(const TParticle&) const;
/**
In this function SIBYLL is called to produce one event. The
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/Processes/StackInspector/StackInspector.cc b/Processes/StackInspector/StackInspector.cc
index da4a91baafaf6253137956d57297cf198ba238c4..6bd1bceb3c912abbaf971ba18039504a6422c0b0 100644
--- a/Processes/StackInspector/StackInspector.cc
+++ b/Processes/StackInspector/StackInspector.cc
@@ -40,7 +40,7 @@ template <typename TStack>
StackInspector<TStack>::~StackInspector() {}
template <typename TStack>
-process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) {
+void StackInspector<TStack>::DoStack(const TStack& vS) {
[[maybe_unused]] int i = 0;
HEPEnergyType Etot = 0_GeV;
@@ -63,7 +63,7 @@ process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) {
const std::chrono::duration<double> elapsed_seconds = now - StartTime_;
std::time_t const now_time = std::chrono::system_clock::to_time_t(now);
auto const dE = E0_ - Etot;
- if (dE < dE_threshold_) return process::EProcessReturn::eOk;
+ if (dE < dE_threshold_) return;
double const progress = dE / E0_;
double const eta_seconds = elapsed_seconds.count() / progress;
@@ -77,7 +77,7 @@ process::EProcessReturn StackInspector<TStack>::DoStack(const TStack& vS) {
<< ", nStep=" << GetStep() << ", stackEntries=" << vS.getEntries()
<< ", Estack=" << Etot / 1_GeV << " GeV"
<< ", ETA=" << std::put_time(std::localtime(&eta_time), "%T") << endl;
- return process::EProcessReturn::eOk;
+ return;
}
#include <corsika/cascade/testCascade.h>
diff --git a/Processes/StackInspector/StackInspector.h b/Processes/StackInspector/StackInspector.h
index 280ff15e6b97117d9936996a1953490063d688d2..6d99a4164844c4fae10fb0e29b15afc15a67d8aa 100644
--- a/Processes/StackInspector/StackInspector.h
+++ b/Processes/StackInspector/StackInspector.h
@@ -30,7 +30,7 @@ namespace corsika::process {
const corsika::units::si::HEPEnergyType vE0);
~StackInspector();
- EProcessReturn DoStack(const TStack&);
+ void DoStack(const TStack&);
/**
* To set a new E0, for example when a new shower event is started
diff --git a/Processes/StackInspector/testStackInspector.cc b/Processes/StackInspector/testStackInspector.cc
index 89edc2a44e8f56f2d9022f2d0ecb688a7e41881f..d76d9fd09b0bd57b37f81daa3f57106eea151910 100644
--- a/Processes/StackInspector/testStackInspector.cc
+++ b/Processes/StackInspector/testStackInspector.cc
@@ -46,7 +46,6 @@ TEST_CASE("StackInspector", "[processes]") {
SECTION("interface") {
StackInspector<TestCascadeStack> model(1, true, E0);
-
- [[maybe_unused]] const process::EProcessReturn ret = model.DoStack(stack);
+ model.DoStack(stack);
}
}
diff --git a/Processes/TrackingLine/CMakeLists.txt b/Processes/TrackingLine/CMakeLists.txt
index 8411487ff9df4088ea26b7754f1fcd0c5fd0867f..dd626f8778c5a4139fe042fc651bb5e20eb315fe 100644
--- a/Processes/TrackingLine/CMakeLists.txt
+++ b/Processes/TrackingLine/CMakeLists.txt
@@ -21,7 +21,6 @@ set_target_properties (
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
-# PUBLIC_HEADER "${MODEL_HEADERS}"
)
# target dependencies on other libraries (also the header onlys)
@@ -33,6 +32,7 @@ target_link_libraries (
CORSIKAunits
CORSIKAenvironment
CORSIKAgeometry
+ CORSIKAlogging
)
target_include_directories (
diff --git a/Processes/TrackingLine/TrackingLine.cc b/Processes/TrackingLine/TrackingLine.cc
index ee76db62c7db07a0988a462367808530f6c636df..23efd27760e873d6cfd9ba564008bd78a95012af 100644
--- a/Processes/TrackingLine/TrackingLine.cc
+++ b/Processes/TrackingLine/TrackingLine.cc
@@ -12,6 +12,7 @@
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
+#include <corsika/logging/Logging.h>
#include <limits>
#include <stdexcept>
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 48bdc3a45da5d1fafcbc3ae144b7916ebd601797..5e20d920c7e0f875618b06943c244ec0e23d329f 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -14,6 +14,8 @@
#include <corsika/geometry/Trajectory.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/logging/Logging.h>
+
#include <optional>
#include <type_traits>
#include <utility>
@@ -38,7 +40,7 @@ namespace corsika::process {
class TrackingLine {
public:
- TrackingLine(){};
+ TrackingLine() = default;
template <typename Particle> // was Stack previously, and argument was
// Stack::StackIterator
@@ -49,27 +51,22 @@ namespace corsika::process {
p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
auto const currentPosition = p.GetPosition();
- std::cout << "TrackingLine pid: " << p.GetPID()
- << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
- std::cout << "TrackingLine pos: "
- << currentPosition.GetCoordinates()
- // << " [" << p.GetNode()->GetModelProperties().GetName() << "]"
- << std::endl;
- std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
- std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
- << " GeV " << std::endl;
- std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
+ C8LOG_DEBUG(
+ "TrackingLine pid: {}"
+ " , E = {} GeV",
+ p.GetPID(), p.GetEnergy() / 1_GeV);
+ C8LOG_DEBUG("TrackingLine pos: {}", currentPosition.GetCoordinates());
+ C8LOG_DEBUG("TrackingLine E: {} GeV", p.GetEnergy() / 1_GeV);
+ C8LOG_DEBUG("TrackingLine p: {} GeV", p.GetMomentum().GetComponents() / 1_GeV);
+ C8LOG_DEBUG("TrackingLine v: {} ", velocity.GetComponents());
// to do: include effect of magnetic field
geometry::Line line(currentPosition, velocity);
auto const* currentLogicalVolumeNode = p.GetNode();
- //~ auto const* currentNumericalVolumeNode =
- //~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
- auto const numericallyInside =
- currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
- std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false");
+ C8LOG_DEBUG("numericallyInside = {} ",
+ currentLogicalVolumeNode->GetVolume().Contains(currentPosition));
auto const& children = currentLogicalVolumeNode->GetChildNodes();
auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
@@ -85,14 +82,11 @@ namespace corsika::process {
if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
auto const [t1, t2] = *opt;
- std::cout << "intersection times: " << t1 / 1_s << "; "
- << t2 / 1_s
- // << " " << vtn.GetModelProperties().GetName()
- << std::endl;
+ C8LOG_DEBUG("intersection times: {} s; {} s", t1 / 1_s, t2 / 1_s);
if (t1.magnitude() > 0)
intersections.emplace_back(t1, &vtn);
else if (t2.magnitude() > 0)
- std::cout << "inside other volume" << std::endl;
+ C8LOG_DEBUG("inside other volume");
}
};
@@ -123,10 +117,7 @@ namespace corsika::process {
min = minIter->first;
}
- std::cout << " t-intersect: "
- << min
- // << " " << minIter->second->GetModelProperties().GetName()
- << std::endl;
+ C8LOG_DEBUG(" t-intersect: {} ", min);
return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
velocity.norm() * min, minIter->second);
diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc
index 26814293241bfb877794687e7b6cf27790c947bb..9d78d20ae7f209cdc62d0c368ee336a70dabd49a 100644
--- a/Processes/UrQMD/UrQMD.cc
+++ b/Processes/UrQMD/UrQMD.cc
@@ -19,7 +19,6 @@
#include <cmath>
#include <fstream>
#include <functional>
-#include <iostream>
#include <random>
#include <sstream>
@@ -40,10 +39,6 @@ CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code projectileCode,
HEPEnergyType labEnergy) const {
// translated to C++ from CORSIKA 7 subroutine cxtot_u
- C8LOG_DEBUG("UrQMD::GetTabulatedCrossSection proj={}, targ={}, E={}GeV",
- particles::GetName(projectileCode), particles::GetName(targetCode),
- labEnergy / 1_GeV);
-
auto const kinEnergy = labEnergy - particles::GetMass(projectileCode);
assert(kinEnergy >= HEPEnergyType::zero());
@@ -91,10 +86,10 @@ CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code projectileCode,
case particles::Code::K0Bar:
projectileIndex = 8;
break;
- default:
- std::cout << "WARNING: UrQMD cross-section not tabulated for " << projectileCode
- << std::endl;
+ default: {
+ C8LOG_WARN("WARNING: UrQMD cross-section not tabulated for {}", projectileCode);
return CrossSectionType::zero();
+ }
}
int targetIndex;
@@ -120,6 +115,10 @@ CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code projectileCode,
xs_interp_support_table_[projectileIndex][targetIndex][je + i - 1 - 1] * w[i];
}
+ C8LOG_DEBUG("UrQMD::GetTabulatedCrossSection proj={}, targ={}, E={}GeV, sigma={}",
+ particles::GetName(projectileCode), particles::GetName(targetCode),
+ labEnergy / 1_GeV, result);
+
return result;
}
@@ -256,6 +255,9 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupView& view) {
auto const& projectilePosition = projectile.GetPosition();
auto const projectileTime = projectile.GetTime();
+ C8LOG_DEBUG("UrQMD::DoInteraction pid={} E={} GeV", projectileCode,
+ projectileEnergyLab / 1_GeV);
+
// sample target particle
auto const& mediumComposition =
projectile.GetNode()->GetModelProperties().GetNuclearComposition();
@@ -348,14 +350,14 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupView& view) {
auto const energy = sqrt(momentum.squaredNorm() + square(particles::GetMass(code)));
momentum.rebase(originalCS); // transform back into standard lab frame
- std::cout << i << " " << code << " " << momentum.GetComponents() << std::endl;
+ C8LOG_DEBUG(" Secondary {} code {} p={} GeV", i, code,
+ momentum.GetComponents() / 1_GeV);
view.AddSecondary(
- std::tuple<particles::Code, HEPEnergyType, stack::MomentumVector, geometry::Point,
- TimeType>{code, energy, momentum, projectilePosition, projectileTime});
+ std::make_tuple(code, energy, momentum, projectilePosition, projectileTime));
}
- std::cout << "UrQMD generated " << sys_.npart << " secondaries!" << std::endl;
+ C8LOG_DEBUG("UrQMD generated {} secondaries!", sys_.npart);
return process::EProcessReturn::eOk;
}
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()));
+ }
}