IAP GITLAB

Commit 86808ad9 authored by Ralf Ulrich's avatar Ralf Ulrich

Resolve "Handling of boundary crossings in geometry tree"

parent cd9f46dc
......@@ -34,7 +34,6 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
ProcessStackInspector
ProcessTrackWriter
ProcessTrackingLine
ProcessHadronicElasticModel
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
......@@ -45,6 +44,24 @@ target_link_libraries (cascade_example SuperStupidStack CORSIKAunits
install (TARGETS cascade_example DESTINATION share/examples)
CORSIKA_ADD_TEST (cascade_example)
add_executable (boundary_example boundary_example.cc)
target_compile_options(boundary_example PRIVATE -g) # do not skip asserts
target_link_libraries (boundary_example SuperStupidStack CORSIKAunits CORSIKAlogging
CORSIKArandom
ProcessSibyll
CORSIKAcascade
ProcessTrackWriter
ProcessTrackingLine
CORSIKAprocesses
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
install (TARGETS boundary_example DESTINATION share/examples)
CORSIKA_ADD_TEST (boundary_example)
add_executable (cascade_proton_example cascade_proton_example.cc)
target_compile_options(cascade_proton_example PRIVATE -g) # do not skip asserts
target_link_libraries (cascade_proton_example SuperStupidStack CORSIKAunits
......@@ -63,6 +80,7 @@ target_link_libraries (cascade_proton_example SuperStupidStack CORSIKAunits
CORSIKAenvironment
CORSIKAprocesssequence
)
install (TARGETS cascade_proton_example DESTINATION share/examples)
CORSIKA_ADD_TEST (cascade_proton_example)
......
This diff is collapsed.
......@@ -12,10 +12,10 @@
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
......@@ -39,9 +39,6 @@
#include <corsika/utl/CorsikaFenv.h>
#include <boost/type_index.hpp>
using boost::typeindex::type_id_with_cvr;
#include <iostream>
#include <limits>
#include <typeinfo>
......@@ -229,28 +226,36 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
environment::Environment env;
using EnvType = environment::Environment<setup::IEnvironmentModel>;
EnvType env;
auto& universe = *(env.GetUniverse());
auto theMedium = environment::Environment::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;
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
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}));
universe.AddChild(std::move(theMedium));
auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
innerMedium->SetModelProperties(props);
outerMedium->AddChild(std::move(innerMedium));
universe.AddChild(std::move(outerMedium));
// setup processes, decays and interactions
tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
tracking_line::TrackingLine tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(true);
const std::vector<particles::Code> trackedHadrons = {
......@@ -259,29 +264,16 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
process::sibyll::Interaction sibyll(env);
process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
process::sibyll::Decay decay(trackedHadrons);
// random::RNGManager::GetInstance().RegisterRandomStream("pythia");
// process::pythia::Decay decay(trackedHadrons);
ProcessCut cut(20_GeV);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
process::TrackWriter::TrackWriter trackWriter("tracks.dat");
process::EnergyLoss::EnergyLoss eLoss;
// assemble all processes into an ordered process list
// auto sequence = stackInspect << sibyll << decay << hadronicElastic << cut <<
// trackWriter; auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss <<
// cut << trackWriter;
// auto sequence = sibyll << sibyllNuc << decay << eLoss << cut;
auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
auto sequence = sibyll << sibyllNuc << decay << cut << trackWriter;
// setup particle stack, and add primary particle
setup::Stack stack;
......
......@@ -229,12 +229,13 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
environment::Environment env;
using EnvType = environment::Environment<setup::IEnvironmentModel>;
EnvType env;
auto& universe = *(env.GetUniverse());
auto theMedium = environment::Environment::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
auto theMedium =
EnvType::CreateNode<Sphere>(Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
......@@ -248,18 +249,17 @@ int main() {
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
// setup processes, decays and interactions
tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
tracking_line::TrackingLine tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(true);
const std::vector<particles::Code> trackedHadrons = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
// process::sibyll::Interaction sibyll(env);
process::pythia::Interaction pythia(env);
process::pythia::Interaction pythia;
// process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
// process::sibyll::Decay decay(trackedHadrons);
process::pythia::Decay decay(trackedHadrons);
......@@ -274,8 +274,8 @@ int main() {
// assemble all processes into an ordered process list
// auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
// auto sequence = p0 << sibyll << sibyllNuc << decay << cut << trackWriter;
auto sequence = p0 << pythia << decay << cut << trackWriter;
auto sequence = stackInspect << pythia << decay << cut << trackWriter;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
......@@ -285,7 +285,7 @@ int main() {
stack.Clear();
const Code beamCode = Code::Proton;
const HEPMassType mass = particles::Proton::GetMass();
const HEPEnergyType E0 = 100_GeV;
const HEPEnergyType E0 = 100_GeV;
double theta = 0.;
double phi = 0.;
......@@ -305,10 +305,10 @@ int main() {
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});
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
......
......@@ -9,6 +9,7 @@ set (
LinearApproximationIntegrator.h
DensityFunction.h
Environment.h
NameModel.h
FlatExponential.h
)
......
......@@ -9,20 +9,17 @@
* the license.
*/
#ifndef _include_Environment_h
#define _include_Environment_h
#ifndef _include_environment_Environment_h
#define _include_environment_Environment_h
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/setup/SetupEnvironment.h>
#include <limits>
namespace corsika::environment {
using BaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
struct Universe : public corsika::geometry::Sphere {
Universe(corsika::geometry::CoordinateSystem const& pCS)
: corsika::geometry::Sphere(
......@@ -34,16 +31,18 @@ namespace corsika::environment {
bool Contains(corsika::geometry::Point const&) const override { return true; }
};
// template <typename IEnvironmentModel>
template <typename IEnvironmentModel>
class Environment {
public:
using BaseNodeType = VolumeTreeNode<IEnvironmentModel>;
Environment()
: fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem()}
, fUniverse(std::make_unique<BaseNodeType>(
std::make_unique<Universe>(fCoordinateSystem))) {}
using IEnvironmentModel = corsika::setup::IEnvironmentModel;
// using IEnvironmentModel = corsika::setup::IEnvironmentModel;
auto& GetUniverse() { return fUniverse; }
auto const& GetUniverse() const { return fUniverse; }
......@@ -63,9 +62,12 @@ namespace corsika::environment {
private:
corsika::geometry::CoordinateSystem const& fCoordinateSystem;
BaseNodeType::VTNUPtr fUniverse;
typename BaseNodeType::VTNUPtr fUniverse;
};
// using SetupBaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
// using SetupEnvironment = Environment<corsika::setup::IEnvironmentModel>;
} // namespace corsika::environment
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
......@@ -9,29 +8,20 @@
* the license.
*/
#ifndef _include_process_trackinling_teststack_h_
#define _include_process_trackinling_teststack_h_
typedef corsika::units::si::hepmomentum_d MOMENTUM;
#ifndef _include_NameModel_h
#define _include_NameModel_h
struct DummyParticle {
HEPEnergyType fEnergy;
Vector<MOMENTUM> fMomentum;
Point fPosition;
#include <string>
#include <utility>
DummyParticle(HEPEnergyType pEnergy, Vector<MOMENTUM> pMomentum, Point pPosition)
: fEnergy(pEnergy)
, fMomentum(pMomentum)
, fPosition(pPosition) {}
namespace corsika::environment {
auto GetEnergy() const { return fEnergy; }
auto GetMomentum() const { return fMomentum; }
auto GetPosition() const { return fPosition; }
auto GetPID() const { return corsika::particles::Code::Unknown; }
};
template <typename T>
struct NameModel : public T {
virtual std::string const& GetName() const = 0;
virtual ~NameModel() = default;
};
struct DummyStack {
using ParticleType = DummyParticle;
};
} // namespace corsika::environment
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
......@@ -12,6 +11,7 @@
#ifndef _include_VolumeTreeNode_H
#define _include_VolumeTreeNode_H
#include <corsika/environment/IMediumModel.h>
#include <corsika/geometry/Volume.h>
#include <memory>
#include <vector>
......@@ -20,9 +20,10 @@ namespace corsika::environment {
class Empty {}; //<! intended for usage as default template argument
template <typename IModelProperties = Empty>
template <typename TModelProperties = Empty>
class VolumeTreeNode {
public:
using IModelProperties = TModelProperties;
using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>;
using IMPSharedPtr = std::shared_ptr<IModelProperties>;
using VolUPtr = std::unique_ptr<corsika::geometry::Volume>;
......@@ -30,6 +31,7 @@ namespace corsika::environment {
VolumeTreeNode(VolUPtr pVolume = nullptr)
: fGeoVolume(std::move(pVolume)) {}
//! convenience function equivalent to Volume::Contains
bool Contains(corsika::geometry::Point const& p) const {
return fGeoVolume->Contains(p);
}
......@@ -44,7 +46,7 @@ namespace corsika::environment {
}
/** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given
* \class Point \arg p, or nullptr iff \arg p is not contained in this volume.
* \class Point \p p, or nullptr iff \p p is not contained in this volume.
*/
VolumeTreeNode<IModelProperties> const* GetContainingNode(
corsika::geometry::Point const& p) const {
......@@ -106,14 +108,14 @@ namespace corsika::environment {
auto const& GetModelProperties() const { return *fModelProperties; }
IMPSharedPtr GetModelPropertiesPtr() const { return fModelProperties; }
auto HasModelProperties() const { return fModelProperties != nullptr; }
template <typename TModelProperties, typename... Args>
template <typename ModelProperties, typename... Args>
auto SetModelProperties(Args&&... args) {
static_assert(std::is_base_of_v<IModelProperties, TModelProperties>,
static_assert(std::is_base_of_v<IModelProperties, ModelProperties>,
"unusable type provided");
fModelProperties = std::make_shared<TModelProperties>(std::forward<Args>(args)...);
fModelProperties = std::make_shared<ModelProperties>(std::forward<Args>(args)...);
return fModelProperties;
}
......
......@@ -9,6 +9,7 @@ set (
set (
CORSIKAcascade_HEADERS
Cascade.h
testCascade.h
)
add_library (CORSIKAcascade INTERFACE)
......
......@@ -20,9 +20,11 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <cassert>
#include <cmath>
#include <iostream>
#include <limits>
#include <type_traits>
/**
* The cascade namespace assembles all objects needed to simulate full particles cascades.
......@@ -38,14 +40,12 @@ namespace corsika::cascade {
* plugged into the cascade simulation.
*
* <b>Tracking</b> must be a class according to the
* TrackingInterface providing the functions: <code>void
* Init();</code> and <code>auto GetTrack(Particle const& p)</auto>,
* where the latter has a return type of <code>
* geometry::Trajectory<corsika::geometry::Line or Helix> </code>
*
* <b>ProcessList</b> must be a ProcessSequence.
* TimeOfIntersection(corsika::geometry::Line const& line,
* TrackingInterface providing the functions:
* <code>auto GetTrack(Particle const& p)</auto>,
* with the return type <code>geometry::Trajectory<corsika::geometry::Line>
* </code>
*
* <b>ProcessList</b> must be a ProcessSequence. *
* <b>Stack</b> is the storage object for particle data, i.e. with
* Particle class type <code>Stack::ParticleType</code>
*
......@@ -56,6 +56,9 @@ namespace corsika::cascade {
template <typename Tracking, typename ProcessList, typename Stack>
class Cascade {
using Particle = typename Stack::ParticleType;
using VolumeTreeNode =
std::remove_pointer_t<decltype(((Particle*)nullptr)->GetNode())>;
using MediumInterface = typename VolumeTreeNode::IModelProperties;
// we only want fully configured objects
Cascade() = delete;
......@@ -65,8 +68,8 @@ namespace corsika::cascade {
* Cascade class cannot be default constructed, but needs a valid
* list of physics processes for configuration at construct time.
*/
Cascade(corsika::environment::Environment const& env, Tracking& tr, ProcessList& pl,
Stack& stack)
Cascade(corsika::environment::Environment<MediumInterface> const& env, Tracking& tr,
ProcessList& pl, Stack& stack)
: fEnvironment(env)
, fTracking(tr)
, fProcessSequence(pl)
......@@ -77,16 +80,29 @@ namespace corsika::cascade {
* All components of the Cascade simulation must be configured here.
*/
void Init() {
fTracking.Init();
fProcessSequence.Init();
fStack.Init();
}
/**
* 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);
});
}
/**
* The Run function is the main simulation loop, which processes
* particles from the Stack until the Stack is empty.
*/
void Run() {
SetNodes();
while (!fStack.IsEmpty()) {
while (!fStack.IsEmpty()) {
auto pNext = fStack.GetNextParticle();
......@@ -113,7 +129,7 @@ namespace corsika::cascade {
using namespace corsika::units::si;
// determine geometric tracking
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
auto [step, geomMaxLength, nextVol] = fTracking.GetTrack(particle);
// determine combined total interaction length (inverse)
InverseGrammageType const total_inv_lambda =
......@@ -126,19 +142,17 @@ namespace corsika::cascade {
std::cout << "total_inv_lambda=" << total_inv_lambda
<< ", next_interact=" << next_interact << std::endl;
// convert next_step from grammage to length
auto const* currentNode =
fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
auto const* currentLogicalNode = particle.GetNode();
if (currentNode == &*fEnvironment.GetUniverse()) {
throw std::runtime_error("particle entered void universe");
}
// assert that particle stays outside void Universe if it has no
// model properties set
assert(currentLogicalNode != &*fEnvironment.GetUniverse() ||
fEnvironment.GetUniverse()->HasModelProperties());
// convert next_step from grammage to length
LengthType const distance_interact =
(next_interact == std::numeric_limits<double>::infinity() * 1_g / (1_m * 1_m))
? std::numeric_limits<double>::infinity() * 1_m
: currentNode->GetModelProperties().ArclengthFromGrammage(step,
next_interact);
currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
next_interact);
// determine the maximum geometric step length
LengthType const distance_max = fProcessSequence.MaxStepLength(particle, step);
......@@ -161,7 +175,7 @@ namespace corsika::cascade {
// take minimum of geometry, interaction, decay for next step
auto const min_distance =
std::min({distance_interact, distance_decay, distance_max});
std::min({distance_interact, distance_decay, distance_max, geomMaxLength});
std::cout << " move particle by : " << min_distance << std::endl;
......@@ -172,13 +186,6 @@ namespace corsika::cascade {
step.LimitEndTo(min_distance);
// particle.GetNode(); // previous VolumeNode
particle.SetNode(
currentNode); // NOTE @Max : here we need to distinguish: IF particle step is
// limited by tracking (via fTracking.GetTrack()), THEN we need
// to check/update VolumeNodes. In all other cases it is
// guaranteed that we are still in the same volume
// apply all continuous processes on particle + track
corsika::process::EProcessReturn status =
fProcessSequence.DoContinuous(particle, step, fStack);
......@@ -191,11 +198,9 @@ namespace corsika::cascade {
}
std::cout << "sth. happening before geometric limit ? "
<< ((min_distance < distance_max) ? "yes" : "no") << std::endl;
if (min_distance < distance_max) { // interaction to happen within geometric limit
// check whether decay or interaction limits this step
<< ((min_distance < geomMaxLength) ? "yes" : "no") << std::endl;
if (min_distance < geomMaxLength) { // interaction to happen within geometric limit
if (min_distance == distance_interact) {
std::cout << "collide" << std::endl;
......@@ -208,7 +213,7 @@ namespace corsika::cascade {
InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
fProcessSequence.SelectInteraction(particle, step, fStack, sample_process,
inv_lambda_count);
} else {
} else if (min_distance == distance_decay) {
std::cout << "decay" << std::endl;
InverseTimeType const actual_decay_time =
fProcessSequence.GetTotalInverseLifetime(particle);
......@@ -218,18 +223,33 @@ namespace corsika::cascade {
const auto sample_process = uniDist(fRNG);
InverseTimeType inv_decay_count = 0 / second;
fProcessSequence.SelectDecay(particle, fStack, sample_process, inv_decay_count);
} else { // step-length limitation within volume
std::cout << "step-length limitation" << std::endl;
}
auto const assertion = [&] {
auto const* numericalNodeAfterStep =
fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
return numericalNodeAfterStep == currentLogicalNode;
};
assert(assertion()); // numerical and logical nodes don't match
} else { // boundary crossing, step is limited by volume boundary
std::cout << "boundary crossing! next node = " << nextVol << std::endl;
particle.SetNode(nextVol);
// DoBoundary may delete the particle (or not)
fProcessSequence.DoBoundaryCrossing(particle, *currentLogicalNode, *nextVol);
}
}
private:
corsika::environment::Environment const& fEnvironment;
corsika::environment::Environment<MediumInterface> const& fEnvironment;
Tracking& fTracking;
ProcessList& fProcessSequence;
Stack& fStack;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
};
}; // namespace corsika::cascade
} // namespace corsika::cascade
......
......@@ -11,7 +11,7 @@
#include <limits>
#include <corsika/environment/Environment.h>
#include <corsika/cascade/testCascade.h>
#include <corsika/cascade/Cascade.h>
......@@ -25,7 +25,6 @@
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
......@@ -46,11 +45,11 @@ using namespace corsika::geometry;
#include <iostream>
using namespace std;
environment::Environment MakeDummyEnv() {
environment::Environment env; // dummy environment
auto MakeDummyEnv() {
T