IAP GITLAB

Commit 50b12aae authored by Ralf Ulrich's avatar Ralf Ulrich

Merge branch 'add_medium_type' into 'master'

Add medium type, as new property (air, water, rock...)

See merge request AirShowerPhysics/corsika!265
parents 03aa719a d96f1a7a
......@@ -7,3 +7,4 @@ flymd.html
flymd.md
Testing
tags
Environment/GeneratedMediaProperties.inc
......@@ -77,7 +77,7 @@ endif (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
# enable warnings and disallow non-standard language
# configure the various build types here, too
# FYI: optimizer flags: -O2 would not trade speed for size, neither O2/3 use fast-math
# debug: O0, relwithdebinfo: 02, release: O3, minsizerel: Os (all defaults)
# debug: O0, relwithdebinfo: 02, release: O3, minsizerel: Os (all defaults), coverage -> O0
set (CMAKE_CXX_FLAGS "-Wall -pedantic -Wextra -Wno-ignored-qualifiers")
set (CMAKE_Fortran_FLAGS "-std=legacy -Wfunction-elimination")
......@@ -104,7 +104,9 @@ if (CMAKE_BUILD_TYPE STREQUAL Coverage)
find_package (Perl REQUIRED)
# compile coverage under -O2 to remove unused functions
add_compile_options ("-O2")
# add_compile_options ("-O2")
# compile coverage under -O0 to avoid any optimization, function elimation etc.
add_compile_options ("-O0")
set (GCOV gcov CACHE STRING "gcov executable" FORCE)
set (LCOV_BIN_DIR "${PROJECT_SOURCE_DIR}/ThirdParty/lcov/bin")
......
......@@ -67,7 +67,6 @@ if (Pythia8_FOUND)
ProcessPythia8
ProcessUrQMD
CORSIKAcascade
ProcessCONEXSourceCut
ProcessEnergyLoss
ProcessTrackWriter
ProcessStackInspector
......@@ -106,7 +105,38 @@ if (Pythia8_FOUND)
ProcessOnShellCheck
ProcessStackInspector
ProcessLongitudinalProfile
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
CORSIKAhistory # for HistoryObservationPlane
)
CORSIKA_ADD_EXAMPLE (hybrid_MC RUN_OPTIONS 4 2 10000.)
target_link_libraries (hybrid_MC
CORSIKAsetup
CORSIKAunits
CORSIKAlogging
CORSIKArandom
CORSIKAhistory
ProcessCONEXSourceCut
ProcessInteractionCounter
ProcessSibyll
ProcessPythia8
ProcessUrQMD
CORSIKAcascade
ProcessPythia8
ProcessObservationPlane
ProcessInteractionCounter
ProcessTrackWriter
ProcessEnergyLoss
ProcessTrackingLine
ProcessParticleCut
ProcessOnShellCheck
ProcessStackInspector
ProcessLongitudinalProfile
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
......@@ -146,7 +176,6 @@ target_link_libraries (em_shower
ProcessPythia8
ProcessUrQMD
CORSIKAcascade
ProcessCONEXSourceCut
ProcessEnergyLoss
ProcessObservationPlane
ProcessInteractionCounter
......
......@@ -45,7 +45,6 @@ using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;
......@@ -95,7 +94,7 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
using EnvType = Environment<setup::IEnvironmentModel>;
using EnvType = setup::Environment;
EnvType env;
auto& universe = *(env.GetUniverse());
......@@ -105,12 +104,16 @@ int main() {
auto world = EnvType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
auto const 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}));
using MyHomogeneousModel =
environment::MediumPropertyModel<environment::UniformMagneticField<
environment::HomogeneousMedium<setup::EnvironmentInterface>>>;
auto const props = world->SetModelProperties<MyHomogeneousModel>(
environment::Medium::AirDry1Atm, Vector(rootCS, 0_T, 0_T, 0_T),
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Proton},
std::vector<float>{1.f}));
// 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);
......
......@@ -45,7 +45,6 @@ using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;
......@@ -57,7 +56,7 @@ using namespace corsika::units::si;
//
int main() {
logging::SetLevel(logging::level::debug);
logging::SetLevel(logging::level::info);
std::cout << "cascade_example" << std::endl;
......@@ -68,33 +67,34 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
using EnvType = environment::Environment<setup::IEnvironmentModel>;
EnvType env;
setup::Environment env;
auto& universe = *(env.GetUniverse());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto outerMedium = EnvType::CreateNode<Sphere>(
auto world = setup::Environment::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel =
environment::MediumPropertyModel<environment::UniformMagneticField<
environment::HomogeneousMedium<setup::EnvironmentInterface>>>;
// fraction of oxygen
const float fox = 0.20946;
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 const props = world->SetModelProperties<MyHomogeneousModel>(
environment::Medium::AirDry1Atm, Vector(rootCS, 0_T, 0_T, 0_T),
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);
auto innerMedium =
setup::Environment::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));
world->AddChild(std::move(innerMedium));
universe.AddChild(std::move(world));
// setup particle stack, and add primary particle
setup::Stack stack;
......
......@@ -49,7 +49,6 @@ using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;
......@@ -70,24 +69,26 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
using EnvType = Environment<setup::IEnvironmentModel>;
using EnvType = setup::Environment;
EnvType env;
auto& universe = *(env.GetUniverse());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto theMedium = EnvType::CreateNode<Sphere>(
Point{rootCS, 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::MediumPropertyModel<environment::UniformMagneticField<
environment::HomogeneousMedium<setup::EnvironmentInterface>>>;
using MyHomogeneousModel = HomogeneousMedium<IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
environment::Medium::AirDry1Atm, geometry::Vector(rootCS, 0_T, 0_T, 1_T),
1_kg / (1_m * 1_m * 1_m),
NuclearComposition(std::vector<particles::Code>{particles::Code::Hydrogen},
std::vector<float>{(float)1.}));
universe.AddChild(std::move(theMedium));
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
......
......@@ -42,7 +42,6 @@ using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;
......@@ -55,6 +54,9 @@ void registerRandomStreams() {
random::RNGManager::GetInstance().SeedAll();
}
template <typename T>
using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
int main(int argc, char** argv) {
logging::SetLevel(logging::level::info);
......@@ -68,11 +70,15 @@ int main(int argc, char** argv) {
registerRandomStreams();
// setup environment, geometry
using EnvType = Environment<setup::IEnvironmentModel>;
using EnvType = setup::Environment;
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
environment::LayeredSphericalAtmosphereBuilder builder{center};
auto builder = environment::make_layered_spherical_atmosphere_builder<
setup::EnvironmentInterface,
MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
environment::Medium::AirDry1Atm,
geometry::Vector{rootCS, 0_T, 0_T, 1_T});
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
......@@ -82,7 +88,6 @@ int main(int argc, char** argv) {
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
builder.addLinearLayer(1e9_cm, 112.8_km);
builder.assemble(env);
// setup particle stack, and add primary particle
......@@ -162,16 +167,16 @@ int main(int argc, char** argv) {
EAS.Run();
cut.ShowResults();
em_continuous.ShowResults();
em_continuous.showResults();
observationLevel.ShowResults();
const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
cut.GetEmEnergy() + em_continuous.GetEnergyLost() +
cut.GetEmEnergy() + em_continuous.energyLost() +
observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.Reset();
cut.Reset();
em_continuous.Reset();
em_continuous.reset();
auto const hists = proposalCounted.GetHistogram();
hists.saveLab("inthist_lab_emShower.npz");
......
This diff is collapsed.
......@@ -53,7 +53,6 @@ using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;
......@@ -76,6 +75,9 @@ void registerRandomStreams(const int seed) {
random::RNGManager::GetInstance().SeedAll(seed);
}
template <typename T>
using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
int main(int argc, char** argv) {
logging::SetLevel(logging::level::info);
......@@ -95,11 +97,15 @@ int main(int argc, char** argv) {
registerRandomStreams(seed);
// setup environment, geometry
using EnvType = Environment<setup::IEnvironmentModel>;
using EnvType = setup::Environment;
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
environment::LayeredSphericalAtmosphereBuilder builder{center};
auto builder = environment::make_layered_spherical_atmosphere_builder<
setup::EnvironmentInterface,
MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
environment::Medium::AirDry1Atm,
geometry::Vector{rootCS, 0_T, 0_T, 1_T});
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
......@@ -109,7 +115,6 @@ int main(int argc, char** argv) {
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
builder.addLinearLayer(1e9_cm, 112.8_km);
builder.assemble(env);
// setup particle stack, and add primary particle
......@@ -245,16 +250,16 @@ int main(int argc, char** argv) {
EAS.Run();
cut.ShowResults();
em_continuous.ShowResults();
em_continuous.showResults();
observationLevel.ShowResults();
const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
cut.GetEmEnergy() + em_continuous.GetEnergyLost() +
cut.GetEmEnergy() + em_continuous.energyLost() +
observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.Reset();
cut.Reset();
em_continuous.Reset();
em_continuous.reset();
auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
urqmdCounted.GetHistogram() + proposalCounted.GetHistogram();
......
add_custom_command (
OUTPUT ${PROJECT_BINARY_DIR}/Environment/GeneratedMediaProperties.inc
COMMAND ${PROJECT_SOURCE_DIR}/Environment/readProperties.py
${PROJECT_SOURCE_DIR}/Environment/properties8.dat
DEPENDS readProperties.py
properties8.dat
WORKING_DIRECTORY
${PROJECT_BINARY_DIR}/Environment
COMMENT "Read NIST properties8 data file and produce C++ source code GeneratedMediaProperties.inc"
VERBATIM
)
set (
ENVIRONMENT_SOURCES
LayeredSphericalAtmosphereBuilder.cc
ShowerAxis.cc
)
......@@ -25,6 +37,10 @@ set (
UniformMagneticField.h
IRefractiveIndexModel.h
UniformRefractiveIndex.h
IMediumPropertyModel.h
MediumProperties.h
MediumPropertyModel.h
${PROJECT_BINARY_DIR}/Environment/GeneratedMediaProperties.inc # this one is auto-generated
)
set (
......@@ -43,6 +59,19 @@ set_target_properties (
PUBLIC_HEADER "${ENVIRONMENT_HEADERS}"
)
# ....................................................
# since GeneratedMediaProperties.inc is an automatically produced file in the build directory,
# create a symbolic link into the source tree, so that it can be found and edited more easily
# this is not needed for the build to succeed! .......
add_custom_command (
OUTPUT ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedParticleProperties.inc
COMMAND ${CMAKE_COMMAND} -E create_symlink ${PROJECT_BINARY_DIR}/include/corsika/environment/GeneratedMediaProperties.inc ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedMediaProperties.inc
COMMENT "Generate link in source-dir: ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedMediaProperties.inc"
)
add_custom_target (SourceDirLinkMedia DEPENDS ${PROJECT_BINARY_DIR}/Environment/GeneratedMediaProperties.inc)
add_dependencies (CORSIKAenvironment SourceDirLinkMedia)
# .....................................................
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
CORSIKAenvironment
......
......@@ -38,8 +38,6 @@ namespace corsika::environment {
, fUniverse(std::make_unique<BaseNodeType>(
std::make_unique<Universe>(fCoordinateSystem))) {}
// using IEnvironmentModel = corsika::setup::IEnvironmentModel;
auto& GetUniverse() { return fUniverse; }
auto const& GetUniverse() const { return fUniverse; }
......@@ -61,7 +59,4 @@ namespace corsika::environment {
typename BaseNodeType::VTNUPtr fUniverse;
};
// using SetupBaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
// using SetupEnvironment = Environment<corsika::setup::IEnvironmentModel>;
} // namespace corsika::environment
/*
* (c) Copyright 2020 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/environment/MediumProperties.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::environment {
/**
* An interface for type of media, needed e.g. to determine energy losses
*
* This is the base interface for media types.
*
*/
template <typename Model>
class IMediumPropertyModel : public Model {
public:
/**
* Evaluate the medium type at a given location.
*
* @param point The location to evaluate at.
* @returns The media type
*/
virtual Medium medium(corsika::geometry::Point const&) const = 0;
/**
* A virtual default destructor.
*/
virtual ~IMediumPropertyModel() = default;
}; // END: class IMediumTypeModel
} // namespace corsika::environment
/*
* (c) Copyright 2020 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/environment/FlatExponential.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/SlidingPlanarExponential.h>
using namespace corsika::environment;
void LayeredSphericalAtmosphereBuilder::checkRadius(units::si::LengthType r) const {
if (r <= previousRadius_) {
throw std::runtime_error("radius must be greater than previous");
}
}
void LayeredSphericalAtmosphereBuilder::setNuclearComposition(
NuclearComposition composition) {
composition_ = std::make_unique<NuclearComposition>(composition);
}
void LayeredSphericalAtmosphereBuilder::addExponentialLayer(
units::si::GrammageType b, units::si::LengthType c,
units::si::LengthType upperBoundary) {
auto const radius = earthRadius_ + upperBoundary;
checkRadius(radius);
previousRadius_ = radius;
auto node = std::make_unique<VolumeTreeNode<IMediumModel>>(
std::make_unique<geometry::Sphere>(center_, radius));
auto const rho0 = b / c;
std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
node->SetModelProperties<SlidingPlanarExponential<IMediumModel>>(
center_, rho0, -c, *composition_, earthRadius_);
layers_.push(std::move(node));
}
void LayeredSphericalAtmosphereBuilder::addLinearLayer(
units::si::LengthType c, units::si::LengthType upperBoundary) {
using namespace units::si;
auto const radius = earthRadius_ + upperBoundary;
checkRadius(radius);
previousRadius_ = radius;
units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
auto const rho0 = b / c;
std::cout << "rho0 = " << rho0;
auto node = std::make_unique<VolumeTreeNode<environment::IMediumModel>>(
std::make_unique<geometry::Sphere>(center_, radius));
node->SetModelProperties<HomogeneousMedium<IMediumModel>>(rho0, *composition_);
layers_.push(std::move(node));
}
Environment<IMediumModel> LayeredSphericalAtmosphereBuilder::assemble() {
Environment<IMediumModel> env;
assemble(env);
return env;
}
void LayeredSphericalAtmosphereBuilder::assemble(Environment<IMediumModel>& env) {
auto& universe = env.GetUniverse();
auto* outmost = universe.get();
while (!layers_.empty()) {
auto l = std::move(layers_.top());
auto* tmp = l.get();
outmost->AddChild(std::move(l));
layers_.pop();
outmost = tmp;
}
}
......@@ -7,52 +7,213 @@
*/
#include <corsika/environment/Environment.h>
#include <corsika/environment/FlatExponential.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/SlidingPlanarExponential.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h>
#include <functional>
#include <memory>
#include <stack>
#include <tuple>
#include <type_traits>
namespace corsika::environment {
// forward-decl
template <typename TMediumInterface, template <typename> typename MExtraEnvirnoment>
struct make_layered_spherical_atmosphere_builder;
/**
* Helper class to setup concentric spheres of layered atmosphere
* with spcified density profiles (exponential, linear, ...).
*
* This can be used most importantly to replicate CORSIKA7
* atmospheres.
*
* Each layer by definition has a density profile and a (constant)
* nuclear composition model.
*
*/
namespace detail {
struct NoExtraModelInner {};
template <typename M>
struct NoExtraModel {};
template <template <typename> typename M>
struct has_extra_models : std::true_type {};
template <>
struct has_extra_models<NoExtraModel> : std::false_type {};
} // namespace detail
template <typename TMediumInterface = environment::IMediumModel,
template <typename> typename TMediumModelExtra = detail::NoExtraModel,
typename... TModelArgs>
class LayeredSphericalAtmosphereBuilder {
std::unique_ptr<NuclearComposition> composition_;
geometry::Point center_;
units::si::LengthType previousRadius_{units::si::LengthType::zero()};
units::si::LengthType earthRadius_;