IAP GITLAB

Commit 3dfcaef7 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus 🖖

begin making Environment templated

parent b47afb90
......@@ -16,8 +16,10 @@
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/NameModel.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
......@@ -214,16 +216,25 @@ public:
HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
};
struct MyBoundaryCrossingProcess
: public BoundaryCrossingProcess<MyBoundaryCrossingProcess> {
template <typename Particle>
EProcessReturn DoBoundaryCrossing(Particle&, environment::BaseNodeType const& from,
environment::BaseNodeType const& to) {
std::cout << "boundary crossing! from:" << &from << "; to: " << &to << std::endl;
return EProcessReturn::eOk;
}
void Init() {}
struct MyBoundaryCrossingProcess : public BoundaryCrossingProcess<MyBoundaryCrossingProcess> {
//~ environment::BaseNodeType const& fA, fB;
MyBoundaryCrossingProcess() {}
//~ MyBoundaryCrossingProcess(environment::BaseNodeType const& a, environment::BaseNodeType const& b) : fA(a), fB(b) {}
template <typename Particle>
EProcessReturn DoBoundaryCrossing(Particle& p, environment::BaseNodeType const& from, environment::BaseNodeType const& to) {
std::cout << "boundary crossing! from:" << &from << "; to: " << &to << std::endl;
//~ if ((&fA == &from && &fB == &to) || (&fA == &to && &fB == &from)) {
p.Delete();
//~ }
return EProcessReturn::eOk;
}
void Init() {}
};
//
......@@ -244,24 +255,25 @@ int main() {
// fraction of oxygen
const float fox = 0.20946;
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
outerMedium->SetModelProperties<MyHomogeneousModel>(
using MyHomogeneousModel = environment::HomogeneousMedium<setup::IEnvironmentModel>;
outerMedium->SetModelProperties<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 = environment::Environment::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 10_m);
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
2000_m);
innerMedium->SetModelProperties<MyHomogeneousModel>(
innerMedium->SetModelProperties<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}));
outerMedium->AddChild(std::move(innerMedium));
universe.AddChild(std::move(outerMedium));
......@@ -279,26 +291,26 @@ int main() {
ProcessCut cut(20_GeV);
random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
process::HadronicElasticModel::HadronicElasticInteraction hadronicElastic(env);
process::HadronicElasticModel::HadronicElasticInteraction
hadronicElastic(env);
process::TrackWriter::TrackWriter trackWriter("tracks.dat");
// assemble all processes into an ordered process list
// auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
auto sequence = p0 << sibyll << sibyllNuc << decay << cut << MyBoundaryCrossingProcess()
<< trackWriter;
auto sequence = p0 /*<< sibyll << sibyllNuc << decay << cut*/ << hadronicElastic << MyBoundaryCrossingProcess() << trackWriter;
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
const Code beamCode = Code::Nucleus;
const Code beamCode = Code::Proton;
const int nuclA = 56;
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 * 100_GeV;
double theta = 27.234;
double phi = 0.;
const HEPEnergyType E0 = nuclA * 100_TeV;
double theta = 0;
double phi = 0;
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
......@@ -316,10 +328,12 @@ 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);
for (int l = 0; l < 100; ++l) {
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ});
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
)
set (
......
......@@ -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
*
* See file AUTHORS for a list of contributors.
*
* 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.
*/
#ifndef _include_NameModel_h
#define _include_NameModel_h
#include <string>
#include <utility>
namespace corsika::environment {
template <typename T>
struct NameModel : public T {
template <typename... Args>
NameModel(std::string const& name, Args&&... args) : T(std::forward<Args>(args)...), fName(name) {}
std::string const& GetName() const {
return fName;
}
private:
std::string fName;
};
} // namespace corsika::environment
#endif
......@@ -21,9 +21,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>;
......@@ -31,6 +32,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);
}
......@@ -45,7 +47,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 {
......@@ -89,12 +91,12 @@ namespace corsika::environment {
auto const& GetModelProperties() const { return *fModelProperties; }
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)
......
......@@ -22,6 +22,7 @@
#include <cmath>
#include <iostream>
#include <type_traits>
/**
* The cascade namespace assembles all objects needed to simulate full particles cascades.
......@@ -53,6 +54,8 @@ 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;
......@@ -62,7 +65,7 @@ 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,
Cascade(corsika::environment::Environment<MediumInterface> const& env, Tracking& tr, ProcessList& pl,
Stack& stack)
: fEnvironment(env)
, fTracking(tr)
......@@ -239,7 +242,7 @@ namespace corsika::cascade {
}
private:
corsika::environment::Environment const& fEnvironment;
corsika::environment::Environment<MediumInterface> const& fEnvironment;
Tracking& fTracking;
ProcessList& fProcessSequence;
Stack& fStack;
......
......@@ -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,12 @@ using namespace corsika::geometry;
#include <iostream>
using namespace std;
environment::Environment MakeDummyEnv() {
environment::Environment env; // dummy environment
auto MakeDummyEnv() {
TestEnvironmentType env; // dummy environment
auto& universe = *(env.GetUniverse());
auto theMedium = environment::Environment::CreateNode<Sphere>(
auto theMedium = TestEnvironmentType::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
......@@ -113,14 +113,14 @@ TEST_CASE("Cascade", "[Cascade]") {
rmng.RegisterRandomStream("cascade");
auto env = MakeDummyEnv();
tracking_line::TrackingLine<setup::Stack, setup::Trajectory> tracking(env);
tracking_line::TrackingLine tracking;
stack_inspector::StackInspector<setup::Stack> p0(true);
stack_inspector::StackInspector<TestCascadeStack> p0(true);
const HEPEnergyType Ecrit = 85_MeV;
ProcessSplit p1(Ecrit);
auto sequence = p0 << p1;
setup::Stack stack;
TestCascadeStack stack;
cascade::Cascade EAS(env, tracking, sequence, stack);
CoordinateSystem const& rootCS =
......
#ifndef _Framework_Cascade_testCascade_h
#define _Framework_Cascade_testCascade_h
#include <corsika/environment/Environment.h>
#include <corsika/setup/SetupStack.h>
using TestEnvironmentType = corsika::environment::Environment<corsika::environment::IMediumModel>;
template <typename T>
using SetupGeometryDataInterface = GeometryDataInterface<T, TestEnvironmentType>;
// combine particle data stack with geometry information for tracking
template <typename StackIter>
using StackWithGeometryInterface =
corsika::stack::CombinedParticleInterface<corsika::setup::detail::ParticleDataStack::PIType,
SetupGeometryDataInterface, StackIter>;
using TestCascadeStack = corsika::stack::CombinedStack<typename corsika::setup::detail::ParticleDataStack::StackImpl, GeometryData<TestEnvironmentType>, StackWithGeometryInterface>;
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
......@@ -13,6 +12,7 @@
#define _include_TRAJECTORY_H
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Line.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::geometry {
......
......@@ -22,11 +22,11 @@ namespace corsika::process {
/**
* This method is called when a particle crosses the boundary between the nodes
* \arg from and \arg to.
* \p from and \p to.
*/
template <typename Particle>
EProcessReturn DoBoundaryCrossing(Particle&, environment::BaseNodeType const& from,
environment::BaseNodeType const& to);
template <typename Particle, typename VTNType>
EProcessReturn DoBoundaryCrossing(Particle&, VTNType const& from,
VTNType const& to);
};
template <class T>
......
......@@ -71,9 +71,9 @@ namespace corsika::process {
// example for a trait-based call:
// void Hello() const { detail::CallHello<T1,T2>::Call(A, B); }
template <typename Particle>
EProcessReturn DoBoundaryCrossing(Particle& p, environment::BaseNodeType const& from,
environment::BaseNodeType const& to) {
template <typename Particle, typename VTNType>
EProcessReturn DoBoundaryCrossing(Particle& p, VTNType const& from,
VTNType const& to) {
EProcessReturn ret = EProcessReturn::eOk;
if constexpr (std::is_base_of<BoundaryCrossingProcess<T1type>, T1type>::value ||
......
......@@ -90,7 +90,7 @@ namespace corsika::stack {
// protected:
/**
@name Access to underlying stack dfata, these are service
@name Access to underlying stack fData, these are service
function for user classes. User code can only rely on GetIndex
and GetStackData to retrieve data
@{
......
......@@ -35,16 +35,14 @@ namespace corsika::process::HadronicElasticModel {
environment::Environment const& env, units::si::CrossSectionType x,
units::si::CrossSectionType y)
: fX(x)
, fY(y)
, fEnvironment(env) {}
, fY(y) {}
template <>
units::si::GrammageType HadronicElasticInteraction::GetInteractionLength(
Particle const& p, Track&) {
using namespace units::si;
if (p.GetPID() == particles::Code::Proton) {
auto const* currentNode =
fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
auto const* currentNode = p.GetNode();
auto const& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
......
......@@ -50,14 +50,12 @@ namespace corsika::process::HadronicElasticModel {
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream(
"HadronicElasticModel");
corsika::environment::Environment const& fEnvironment;
inveV2 B(eV2 s) const;
corsika::units::si::CrossSectionType CrossSection(SquaredHEPEnergyType s) const;
public:
HadronicElasticInteraction(corsika::environment::Environment const&,
// x & y values taken from DL for pp collisions
HadronicElasticInteraction(// x & y values taken from DL for pp collisions
units::si::CrossSectionType x = 0.0217 * units::si::barn,
units::si::CrossSectionType y = 0.05608 * units::si::barn);
void Init();
......
......@@ -34,8 +34,7 @@ using Track = Trajectory;
namespace corsika::process::sibyll {
Interaction::Interaction(environment::Environment const& env)
: fEnvironment(env) {}
Interaction::Interaction() {}
Interaction::~Interaction() {
cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount << endl;
......@@ -127,8 +126,7 @@ namespace corsika::process::sibyll {
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
const auto currentNode =
fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
auto const* currentNode = p.GetNode();
const auto mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// determine average interaction length
......@@ -254,8 +252,8 @@ namespace corsika::process::sibyll {
HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
// sample target mass number
const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
const auto& mediumComposition =
auto const* currentNode = p.GetNode();
auto const& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// get cross sections for target materials
/*
......@@ -268,13 +266,9 @@ namespace corsika::process::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
const auto [sigProd, sigEla, nNuc] =
[[maybe_unsused]] const auto [sigProd, sigEla, nNuc] =
GetCrossSection(corsikaBeamId, targetId, Ecm);
cross_section_of_components[i] = sigProd;
[[maybe_unused]] int ideleteme =
nNuc; // to avoid not used warning in array binding
[[maybe_unused]] auto sigElaCopy =
sigEla; // to avoid not used warning in array binding
}
const auto targetCode = mediumComposition.SampleTarget(
......
......@@ -18,10 +18,6 @@
#include <corsika/units/PhysicalUnits.h>
#include <tuple>
namespace corsika::environment {
class Environment;
}
namespace corsika::process::sibyll {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
......@@ -31,7 +27,7 @@ namespace corsika::process::sibyll {
bool fInitialized = false;
public:
Interaction(corsika::environment::Environment const& env);
Interaction();
~Interaction();
void Init();
......@@ -60,7 +56,6 @@ namespace corsika::process::sibyll {
corsika::process::EProcessReturn DoInteraction(Particle&, Stack&);
private:
corsika::environment::Environment const& fEnvironment;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
};
......
......@@ -34,10 +34,8 @@ using Track = Trajectory;
namespace corsika::process::sibyll {
NuclearInteraction::NuclearInteraction(environment::Environment const& env,
process::sibyll::Interaction& hadint)
: fEnvironment(env)
, fHadronicInteraction(hadint) {}
NuclearInteraction::NuclearInteraction(process::sibyll::Interaction& hadint)
: fHadronicInteraction(hadint) {}
NuclearInteraction::~NuclearInteraction() {
cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount << endl;
......@@ -175,9 +173,8 @@ namespace corsika::process::sibyll {
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
const auto currentNode =
fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
const auto mediumComposition =
auto const* const currentNode = p.GetNode();
auto const& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// determine average interaction length
// weighted sum
......@@ -332,7 +329,7 @@ namespace corsika::process::sibyll {
//
// proton stand-in for nucleon
const auto beamId = particles::Proton::GetCode();
const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
auto const* const currentNode = p.GetNode();
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
cout << "get nucleon-nucleus cross sections for target materials.." << endl;
......
......@@ -16,10 +16,6 @@
#include <corsika/process/InteractionProcess.h>
#include <corsika/random/RNGManager.h>
namespace corsika::environment {
class Environment;
}
namespace corsika::process::sibyll {
class Interaction; // fwd-decl
......@@ -36,8 +32,7 @@ namespace corsika::process::sibyll {
int fNucCount = 0;
public:
NuclearInteraction(corsika::environment::Environment const& env,
corsika::process::sibyll::Interaction& hadint);
NuclearInteraction(corsika::process::sibyll::Interaction& hadint);
~NuclearInteraction();
void Init();
......@@ -52,7 +47,6 @@ namespace corsika::process::sibyll {
corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s);
private:
corsika::environment::Environment const& fEnvironment;
corsika::process::sibyll::Interaction& fHadronicInteraction;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
......
......@@ -85,10 +85,10 @@ using namespace corsika::units;
TEST_CASE("SibyllInterface", "[processes]") {
// setup environment, geometry
environment::Environment env;
environment::Environment<environment::IMediumModel> env;
auto& universe = *(env.GetUniverse());
auto theMedium = environment::Environment::CreateNode<geometry::Sphere>(
auto theMedium = environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
......@@ -123,7 +123,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, pos, 0_ns});
Interaction model(env);
Interaction model;