diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 615c2e1bdd71b2bfc12da6cad825f4157540de78..18c03158ab46f9909d06860f04cb2a173a88ae00 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -216,9 +216,9 @@ public:
HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
};
+template <bool deleteParticle>
struct MyBoundaryCrossingProcess
- : public BoundaryCrossingProcess<MyBoundaryCrossingProcess> {
- //~ environment::BaseNodeType const& fA, fB;
+ : public BoundaryCrossingProcess<MyBoundaryCrossingProcess<deleteParticle>> {
MyBoundaryCrossingProcess(std::string const& filename) { fFile.open(filename); }
@@ -226,12 +226,7 @@ struct MyBoundaryCrossingProcess
EProcessReturn DoBoundaryCrossing(Particle& p,
typename Particle::BaseNodeType const& from,
typename Particle::BaseNodeType const& to) {
- std::cout << "boundary crossing! from: " << from.GetModelProperties().GetName()
- << "; to: " << to.GetModelProperties().GetName() << std::endl;
-
- //~ if ((&fA == &from && &fB == &to) || (&fA == &to && &fB == &from)) {
- p.Delete();
- //~ }
+ std::cout << "boundary crossing! from: " << &from << "; to: " << &to << std::endl;
auto const& name = particles::GetName(p.GetPID());
auto const start = p.GetPosition().GetCoordinates();
@@ -239,25 +234,15 @@ struct MyBoundaryCrossingProcess
fFile << name << " " << start[0] / 1_m << ' ' << start[1] / 1_m << ' '
<< start[2] / 1_m << '\n';
+ if constexpr (deleteParticle) { p.Delete(); }
+
return EProcessReturn::eOk;
}
- std::ofstream fFile;
-
void Init() {}
-};
-template <class T>
-class TheNameModel : public T {
- std::string const fName;
-
-public:
- template <typename... Args>
- TheNameModel(std::string const& name, Args&&... args)
- : T(std::forward<Args>(args)...)
- , fName(name) {}
-
- std::string const& GetName() const override { return fName; }
+private:
+ std::ofstream fFile;
};
//
@@ -273,41 +258,29 @@ int main() {
EnvType env;
auto& universe = *(env.GetUniverse());
- auto outerMedium =
- EnvType::CreateNode<Sphere>(Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
- 1_km * std::numeric_limits<double>::infinity());
+ const CoordinateSystem& rootCS = env.GetCoordinateSystem();
+
+ auto outerMedium = EnvType::CreateNode<Sphere>(
+ Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
// fraction of oxygen
const float fox = 0.20946;
- outerMedium->SetModelProperties<
- TheNameModel<environment::HomogeneousMedium<setup::IEnvironmentModel>>>(
- "outer", 1_kg / (1_m * 1_m * 1_m),
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Nitrogen,
- particles::Code::Oxygen},
- std::vector<float>{(float)1. - fox, fox}));
-
- auto innerMedium = EnvType::CreateNode<Sphere>(
- Point{env.GetCoordinateSystem(), 0_m, 0_m, -500_m}, 5000_m);
-
- innerMedium->SetModelProperties<
- TheNameModel<environment::HomogeneousMedium<setup::IEnvironmentModel>>>(
- "inner", 1_kg / (1_m * 1_m * 1_m),
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Nitrogen,
- particles::Code::Oxygen},
- std::vector<float>{(float)1. - fox, fox}));
+ 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::Nitrogen,
+ particles::Code::Oxygen},
+ std::vector<float>{1.f - fox, fox}));
+
+ auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
+
+ innerMedium->SetModelProperties(props);
outerMedium->AddChild(std::move(innerMedium));
universe.AddChild(std::move(outerMedium));
- universe.SetModelProperties<
- TheNameModel<environment::HomogeneousMedium<setup::IEnvironmentModel>>>(
- "Universe", 0_kg / (1_m * 1_m * 1_m),
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Proton}, std::vector<float>{1.}));
-
- const CoordinateSystem& rootCS = env.GetCoordinateSystem();
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
@@ -319,17 +292,11 @@ int main() {
process::sibyll::Decay decay;
ProcessCut cut(20_GeV);
- random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
- process::HadronicElasticModel::HadronicElasticInteraction hadronicElastic;
-
process::TrackWriter::TrackWriter trackWriter("tracks.dat");
- MyBoundaryCrossingProcess boundaryCrossing("crossings.dat");
+ MyBoundaryCrossingProcess<false> boundaryCrossing("crossings.dat");
// assemble all processes into an ordered process list
- // auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
- auto sequence =
- /*p0 <<*/ sibyll << sibyllNuc << decay << cut /* << hadronicElastic */
- << boundaryCrossing << trackWriter;
+ auto sequence = sibyll << sibyllNuc << decay << cut << boundaryCrossing << trackWriter;
// setup particle stack, and add primary particle
setup::Stack stack;
@@ -339,30 +306,29 @@ int main() {
const int nuclZ = int(nuclA / 2.15 + 0.7);
const HEPMassType mass = particles::Proton::GetMass() * nuclZ +
(nuclA - nuclZ) * particles::Neutron::GetMass();
- const HEPEnergyType E0 = nuclA * 10_TeV;
- double phi = 0;
- for (double theta = 0; theta < 180; theta += 20) {
- auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
- return sqrt(Elab * Elab - m * m);
- };
- HEPMomentumType P0 = elab2plab(E0, mass);
- auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
- return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
- -ptot * cos(theta));
- };
- 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;
- Point pos(rootCS, 0_m, 0_m, 0_m);
- for (int l = 0; l < 1; ++l) {
- stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
- corsika::stack::MomentumVector, geometry::Point,
- units::si::TimeType>{beamCode, E0, plab, pos, 0_ns});
- }
- }
+ const HEPEnergyType E0 = nuclA * 100_TeV;
+ double const phi = 0;
+ double const theta = 0;
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt(Elab * Elab - m * m);
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
+ -ptot * cos(theta));
+ };
+ 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;
+ Point pos(rootCS, 0_m, 0_m, 0_m);
+
+ stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ beamCode, E0, plab, pos, 0_ns});
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index a8185a0f10639e2615b64a06a8eea4f3eb6982ab..2d5088ed343530c945f9f16fbd124440088063e3 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -20,6 +20,7 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
+#include <cassert>
#include <cmath>
#include <iostream>
#include <type_traits>
@@ -91,9 +92,6 @@ namespace corsika::cascade {
auto const* numericalNode =
fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
p.SetNode(numericalNode);
-
- std::cout << "initial node " << p.GetNode()->GetModelProperties().GetName()
- << std::endl;
});
}
@@ -229,19 +227,15 @@ namespace corsika::cascade {
} else { // step-length limitation within volume
std::cout << "step-length limitation" << std::endl;
}
- auto const* numericalNodeAfterStep =
- fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
- std::cout << "nodes: " << currentLogicalNode->GetModelProperties().GetName()
- << " " << numericalNodeAfterStep->GetModelProperties().GetName()
- << std::endl;
+ auto const assertion = [&] {
+ auto const* numericalNodeAfterStep =
+ fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
+ return numericalNodeAfterStep == currentLogicalNode;
+ };
- if (numericalNodeAfterStep != currentLogicalNode) {
- std::cout << "position " << particle.GetPosition().GetCoordinates()
- << std::endl;
- throw std::runtime_error("numerical and logical nodes don't match");
- }
- } else { // boundary crossing
+ assert(assertion()); // numerical and logical nodes don't match
+ } else { // boundary crossing
std::cout << "boundary crossing! next node = " << nextVol << std::endl;
particle.SetNode(nextVol);
fProcessSequence.DoBoundaryCrossing(particle, *currentLogicalNode, *nextVol);
@@ -255,7 +249,7 @@ namespace corsika::cascade {
Stack& fStack;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
- };
+ }; // namespace corsika::cascade
} // namespace corsika::cascade
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 318bb5c2369129ff44e780d3b6fb5c4c2446f7f5..8f9bf45120562a7e31a12683e26e8f3a56931450 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -99,6 +99,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
+ auto const* nodePtr = theMedium.get();
universe.AddChild(std::move(theMedium));
const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
@@ -123,6 +124,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, pos, 0_ns});
+ particle.SetNode(nodePtr);
Interaction model;
@@ -147,6 +149,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2});
+ particle.SetNode(nodePtr);
Interaction hmodel;
NuclearInteraction model(hmodel);
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 4db45da2f7f5730e2123fa0bddd2a64b52bfe98c..e847a3fd7ebef23583c632cc5b71af474c883428 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -43,7 +43,7 @@ namespace corsika::process {
template <typename Particle> // was Stack previously, and argument was
// Stack::StackIterator
- auto GetTrack(Particle const& p) {
+ auto GetTrack(Particle const& p) {
using namespace corsika::units::si;
using namespace corsika::geometry;
geometry::Vector<SpeedType::dimension_type> const velocity =
@@ -52,8 +52,10 @@ namespace corsika::process {
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() << "]\n";
+ 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;
@@ -84,8 +86,10 @@ 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;
+ std::cout << "intersection times: " << t1 / 1_s << "; "
+ << t2 / 1_s
+ // << " " << vtn.GetModelProperties().GetName()
+ << std::endl;
if (t1.magnitude() > 0)
intersections.emplace_back(t1, &vtn);
else if (t2.magnitude() > 0)
@@ -119,10 +123,10 @@ namespace corsika::process {
min = minIter->first;
}
- std::cout << " t-intersect: " << min << " "
- << minIter->second->GetModelProperties().GetName() << std::endl;
- //~ std::cout << "point of intersection: " <<
- //line.GetPosition(min).GetCoordinates() << std::endl;
+ std::cout << " t-intersect: "
+ << min
+ // << " " << minIter->second->GetModelProperties().GetName()
+ << std::endl;
return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
velocity.norm() * min, minIter->second);
diff --git a/Setup/SetupEnvironment.h b/Setup/SetupEnvironment.h
index f55cda3507ec882e10301e3bce95e3e5d1719914..83397da951f2d7594fea3e2c5a3dfa1d2b88da19 100644
--- a/Setup/SetupEnvironment.h
+++ b/Setup/SetupEnvironment.h
@@ -12,13 +12,13 @@
#ifndef _include_corsika_setup_environment_h_
#define _include_corsika_setup_environment_h_
+#include <corsika/environment/Environment.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NameModel.h>
-#include <corsika/environment/Environment.h>
namespace corsika::setup {
- using IEnvironmentModel = environment::NameModel<environment::IMediumModel>;
+ using IEnvironmentModel = environment::IMediumModel;
using SetupEnvironment = environment::Environment<IEnvironmentModel>;
-}
+} // namespace corsika::setup
#endif