IAP GITLAB

Commit 04291651 authored by Ralf Ulrich's avatar Ralf Ulrich

major refactory update

parent 3c70ce61
......@@ -68,17 +68,17 @@ check-clang-format:
allow_failure: true
### CodeQuality tool ####
include:
- template: Code-Quality.gitlab-ci.yml
code_quality:
stage: quality
rules:
- if: '$CODE_QUALITY_DISABLED'
when: never
- if: '$CI_PIPELINE_SOURCE == "merge_request_event"' # Run code quality job in merge request pipelines
- if: '$CI_COMMIT_BRANCH == $CI_DEFAULT_BRANCH' # Run code quality job in pipelines on the master branch (but not in other branch pipelines)
- if: '$CI_COMMIT_TAG' # Run code quality job in pipelines for tags
#include:
# - template: Code-Quality.gitlab-ci.yml
#
#code_quality:
# stage: quality
# rules:
# - if: '$CODE_QUALITY_DISABLED'
# when: never
# - if: '$CI_PIPELINE_SOURCE == "merge_request_event"' # Run code quality job in merge request pipelines
# - if: '$CI_COMMIT_BRANCH == $CI_DEFAULT_BRANCH' # Run code quality job in pipelines on the master branch (but not in other branch pipelines)
# - if: '$CI_COMMIT_TAG' # Run code quality job in pipelines for tags
####### CONFIG ##############
......
cmake_minimum_required (VERSION 3.9)
#
#+++++++++++++++++++++++++++++
# prevent in-source builds and give warning message
#
......@@ -9,7 +9,7 @@ if ("${CMAKE_BINARY_DIR}" STREQUAL "${CMAKE_SOURCE_DIR}")
NOTE: cmake will now create CMakeCache.txt and CMakeFiles/*.
You must delete them, or cmake will refuse to work.")
endif ()
#
#+++++++++++++++++++++++++++++
# project name
#
......@@ -19,20 +19,20 @@ project (
DESCRIPTION "CORSIKA C++ project (alpha status)"
LANGUAGES CXX
)
#
#+++++++++++++++++++++++++++++
# as long as there still are modules using it:
#
enable_language (Fortran)
set (CMAKE_Fortran_FLAGS "-std=legacy -Wfunction-elimination")
#
#+++++++++++++++++++++++++++++
# warn user if system is not UNIX
#
if (NOT UNIX)
message (FATAL_ERROR "| CORSIKA8 > This is an unsupported system.")
endif ()
#
#+++++++++++++++++++++++++++++
# cmake path dir, and cmake config
#
......@@ -42,12 +42,12 @@ include (CorsikaUtilities) # extra cmake function
set (CMAKE_VERBOSE_MAKEFILE OFF) # this can be done with `make VERBOSE=1`
# ignore many irrelevant Up-to-date messages during install
set (CMAKE_INSTALL_MESSAGE LAZY)
#
#+++++++++++++++++++++++++++++
# Setup hardware and infrastructure dependent defines
#
include(CorsikaDefines)
#
#+++++++++++++++++++++++++++++
# check if compiler is C++17 compliant
#
......@@ -59,7 +59,7 @@ endif ()
# set CXX compile flags and options and warning settings
set (CMAKE_CXX_STANDARD 17)
set (CMAKE_CXX_EXTENSIONS OFF)
#
#+++++++++++++++++++++++++++++
# Compiler and linker flags, settings
#
......@@ -70,12 +70,12 @@ set (CMAKE_CXX_EXTENSIONS OFF)
set (CMAKE_CXX_FLAGS "-Wall -pedantic -Wextra -Wno-ignored-qualifiers")
set (CMAKE_Fortran_FLAGS "-std=legacy -Wfunction-elimination")
set (DEFAULT_BUILD_TYPE "Release")
#
# clang produces a lot of unecessary warnings without this:
add_compile_options (
$<$<OR:$<CXX_COMPILER_ID:Clang>,$<CXX_COMPILER_ID:AppleClang>>:-Wno-nonportable-include-path>
)
#
#+++++++++++++++++++++++++++++
# Build types settings
#
......@@ -85,7 +85,7 @@ set (CMAKE_EXE_LINKER_FLAGS_COVERAGE "--coverage")
set (CMAKE_SHARED_LINKER_FLAGS_COVERAGE "--coverage")
# set a flag to inform code that we are in debug mode
set (CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -DDEBUG")
#
#+++++++++++++++++++++++++++++
# Build type selection
#
......@@ -111,7 +111,7 @@ else (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
endif ()
message (STATUS "Build type is: ${CMAKE_BUILD_TYPE}")
endif (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
#
#+++++++++++++++++++++++++++++
# Setup external dependencies.
#
......@@ -137,7 +137,7 @@ conan_cmake_run (CONANFILE conanfile.txt
# add cnpy temporarily. As long as SaveBoostHistogram does not save to parquet directly
#
add_subdirectory (externals/cnpy)
#
#+++++++++++++++++++++++++++++
# Coverage
#
......@@ -172,7 +172,7 @@ if (CMAKE_BUILD_TYPE STREQUAL Coverage)
)
add_custom_target (coverage DEPENDS coverage-report)
endif ()
#
#+++++++++++++++++++++++++++++
# CTest config
#
......@@ -188,25 +188,25 @@ if (NOT DEFINED ENV{CI})
# add call to ./do-copyright.py to run as unit-test-case
add_test (NAME copyright_notices COMMAND ./do-copyright.py WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
endif (NOT DEFINED ENV{CI})
#
#+++++++++++++++++++++++++++++
# Options
#
# HISTORY option selection
option (WITH_HISTORY "Flag to switch on/off HISTORY" ON)
#
#+++++++++++++++++++++++++++++
# Externals
#
set (Python_ADDITIONAL_VERSIONS 3)
find_package (PythonInterp 3 REQUIRED)
#
#+++++++++++++++++++++++++++++
# get phys_units
#
find_package (PhysUnits REQUIRED)
#
#
#+++++++++++++++++++++++++++++
# CORSIKA8
#
......@@ -229,19 +229,19 @@ target_link_libraries (
cnpy # for SaveBoostHistogram
stdc++fs # experimental::filesystem
)
# those are needed, since some headers (namely GeneratedParticleProperties.inc) are produced by python script from ParticleData.xml
# "src" is needed, since some headers (namely GeneratedParticleProperties.inc ec.) are produced by python script from e.g. ParticleData.xml
add_subdirectory (src)
add_subdirectory (documentation)
#
#
#+++++++++++++++++++++++++++++
# =~~~~~~~~~~~~~~~~~~~~~~~~~=
# = Add of subdirectories =
# =~~~~~~~~~~~~~~~~~~~~~~~~~=
#
#+++++++++++++++++++++++++++++
# modules
#
set (CORSIKA_DATA_WITH_TEST ON) # we want to run the corsika-data unit test
add_subdirectory (modules/data)
add_subdirectory (modules/pythia)
add_subdirectory (modules/sibyll)
......@@ -249,16 +249,23 @@ add_subdirectory (modules/qgsjetII)
add_subdirectory (modules/urqmd)
add_subdirectory (modules/conex)
#
# we really need a better proposal CMakeList.txt files:
set (ADD_PYTHON OFF)
file(READ modules/CMakeLists_PROPOSAL.txt overwrite_CMakeLists_PROPOSAL_txt)
file(WRITE modules/proposal/CMakeLists.txt ${overwrite_CMakeLists_PROPOSAL_txt})
add_subdirectory (modules/proposal)
target_link_libraries (CORSIKA8 INTERFACE PROPOSAL)
#+++++++++++++++++++++++++++++++
# unit testing
#
add_subdirectory (tests)
#
#+++++++++++++++++++++++++++++++
# examples
#
add_subdirectory (examples)
#
#+++++++++++++++++++++++++++++++
#
# final summary output
......
/*
* (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/framework/process/InteractionHistogram.hpp>
namespace corsika {
template <class TCountedProcess>
InteractionCounter<TCountedProcess>::InteractionCounter(TCountedProcess& process)
: process_(process) {}
template <class TCountedProcess>
template <typename TSecondaryView>
inline void InteractionCounter<TCountedProcess>::doInteraction(TSecondaryView& view) {
auto const projectile = view.getProjectile();
auto const massNumber = projectile.getNode()
->getModelProperties()
.getNuclearComposition()
.getAverageMassNumber();
auto const massTarget = massNumber * constants::nucleonMass;
if (auto const projectile_id = projectile.getPID(); projectile_id == Code::Nucleus) {
auto const A = projectile.getNuclearA();
auto const Z = projectile.getNuclearZ();
histogram_.fill(projectile_id, projectile.getEnergy(), massTarget, A, Z);
} else {
histogram_.fill(projectile_id, projectile.getEnergy(), massTarget);
}
process_.doInteraction(view);
}
template <class TCountedProcess>
template <typename TParticle>
inline GrammageType InteractionCounter<TCountedProcess>::getInteractionLength(
TParticle const& particle) const {
return process_.getInteractionLength(particle);
}
template <class TCountedProcess>
inline InteractionHistogram const& InteractionCounter<TCountedProcess>::getHistogram()
const {
return histogram_;
}
} // namespace corsika
/*
* (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 <boost/histogram.hpp>
namespace corsika {
namespace detail {
inline auto hist_factory(unsigned int const bin_number, double const e_low,
double const e_high) {
namespace bh = boost::histogram;
namespace bha = bh::axis;
auto h = bh::make_histogram(
bha::category<int, bh::use_default, bha::option::growth_t>{{2212, 2112},
"projectile PDG"},
bha::regular<float, bha::transform::log>{bin_number, (float)e_low,
(float)e_high, "energy/eV"});
return h;
}
} // namespace detail
} // namespace corsika
/*
* (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 <fstream>
#include <string>
#include <corsika/framework/process/InteractionHistogram.hpp>
#include <corsika/detail/framework/process/InteractionHistogram.hpp> // for detail namespace
namespace corsika {
InteractionHistogram::InteractionHistogram()
: inthist_cms_{detail::hist_factory(num_bins_cms, lower_edge_cms, upper_edge_cms)}
, inthist_lab_{detail::hist_factory(num_bins_lab, lower_edge_lab, upper_edge_lab)} {
}
void InteractionHistogram::fill(Code projectile_id, HEPEnergyType lab_energy,
HEPEnergyType mass_target, int A, int Z) {
auto constexpr inv_eV = 1 / 1_eV;
if (projectile_id == Code::Nucleus) {
auto const sqrtS = sqrt(A * A * (constants::nucleonMass * constants::nucleonMass) +
mass_target * mass_target + 2 * lab_energy * mass_target);
int32_t const pdg = 1'000'000'000l + Z * 10'000l + A * 10l;
inthist_lab_(pdg, lab_energy * inv_eV);
inthist_cms_(pdg, sqrtS * inv_eV);
} else {
auto const projectile_mass = get_mass(projectile_id);
auto const sqrtS = sqrt(projectile_mass * projectile_mass +
mass_target * mass_target + 2 * lab_energy * mass_target);
inthist_cms_(static_cast<int>(get_PDG(projectile_id)), sqrtS * inv_eV);
inthist_lab_(static_cast<int>(get_PDG(projectile_id)), lab_energy * inv_eV);
}
}
void InteractionHistogram::saveLab(std::string const& filename, SaveMode mode) const {
corsika::save_hist(inthist_lab_, filename, mode);
}
void InteractionHistogram::saveCMS(std::string const& filename, SaveMode mode) const {
corsika::save_hist(inthist_cms_, filename, mode);
}
InteractionHistogram& InteractionHistogram::operator+=(
InteractionHistogram const& other) {
inthist_lab_ += other.inthist_lab_;
inthist_cms_ += other.inthist_cms_;
return *this;
}
InteractionHistogram InteractionHistogram::operator+(InteractionHistogram other) const {
other.inthist_lab_ += inthist_lab_;
other.inthist_cms_ += inthist_cms_;
return other;
}
} // namespace corsika
\ No newline at end of file
......@@ -102,14 +102,15 @@ namespace corsika {
template <typename TProcess1, typename TProcess2, typename TSelect>
template <typename TParticle, typename TTrack>
inline LengthType SwitchProcessSequence<TProcess1, TProcess2, TSelect>::maxStepLength(
TParticle& particle, TTrack& vTrack) {
inline LengthType SwitchProcessSequence<TProcess1, TProcess2,
TSelect>::getMaxStepLength(TParticle& particle,
TTrack& vTrack) {
switch (select_(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<ContinuousProcess<process1_type>,
process1_type> ||
t1ProcSeq) {
return A_.maxStepLength(particle, vTrack);
return A_.getMaxStepLength(particle, vTrack);
}
break;
}
......@@ -117,7 +118,7 @@ namespace corsika {
if constexpr (std::is_base_of_v<ContinuousProcess<process2_type>,
process2_type> ||
t2ProcSeq) {
return B_.maxStepLength(particle, vTrack);
return B_.getMaxStepLength(particle, vTrack);
}
break;
}
......
......@@ -12,7 +12,7 @@
extern "C" {
#warning No enabling/disabling of floating point exceptions - platform needs better implementation
inline int feenableexcept(int excepts) { return -1; }
inline int feenableexcept(int /*excepts*/) { return -1; }
inline int fedisableexcept(int excepts) { return -1; }
inline int fedisableexcept(int /*excepts*/) { return -1; }
}
/*
* (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/framework/core/ParticleProperties.hpp>
#include <corsika/framework/logging/Logging.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
#include <cmath>
#include <iomanip>
#include <limits>
namespace corsika {
LongitudinalProfile::LongitudinalProfile(ShowerAxis const& shower_axis, GrammageType dX)
: dX_(dX)
, shower_axis_{shower_axis}
, profiles_{static_cast<unsigned int>(shower_axis.getMaximumX() / dX_) + 1} {}
template <typename TParticle, typename TTrack>
ProcessReturn LongitudinalProfile::doContinuous(TParticle const& vP,
TTrack const& vTrack) {
auto const pid = vP.getPID();
GrammageType const grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0));
GrammageType const grammageEnd = shower_axis_.getProjectedX(vTrack.getPosition(1));
CORSIKA_LOG_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_);
for (int b = binStart; b <= binEnd; ++b) {
if (pid == Code::Gamma) {
profiles_.at(b)[ProfileIndex::Gamma]++;
} else if (pid == Code::Positron) {
profiles_.at(b)[ProfileIndex::Positron]++;
} else if (pid == Code::Electron) {
profiles_.at(b)[ProfileIndex::Electron]++;
} else if (pid == Code::MuPlus) {
profiles_.at(b)[ProfileIndex::MuPlus]++;
} else if (pid == Code::MuMinus) {
profiles_.at(b)[ProfileIndex::MuMinus]++;
} else if (is_hadron(pid)) {
profiles_.at(b)[ProfileIndex::Hadron]++;
}
}
return ProcessReturn::Ok;
}
void LongitudinalProfile::save(std::string const& filename, const int width,
const int precision) {
std::ofstream f{filename};
f << "# X / g·cm¯², gamma, e+, e-, mu+, mu-, all hadrons" << std::endl;
for (size_t b = 0; b < profiles_.size(); ++b) {
f << std::setprecision(5) << std::setw(11) << b * (dX_ / (1_g / 1_cm / 1_cm));
for (auto const& N : profiles_.at(b)) {
f << std::setw(width) << std::setprecision(precision) << std::scientific << N;
}
f << std::endl;
}
}
} // namespace corsika
\ No newline at end of file
......@@ -10,10 +10,10 @@
#include <fstream>
namespace corsika::observation_plane {
namespace corsika {
ObservationPlane::ObservationPlane(corsika::Plane const& obsPlane,
std::string const& filename, bool deleteOnHit)
ObservationPlane::ObservationPlane(Plane const& obsPlane, std::string const& filename,
bool deleteOnHit)
: plane_(obsPlane)
, outputStream_(filename)
, deleteOnHit_(deleteOnHit)
......@@ -37,8 +37,8 @@ namespace corsika::observation_plane {
}
const auto energy = particle.getEnergy();
outputStream_ << static_cast<int>(corsika::get_PDG(particle.getPID())) << ' '
<< energy / 1_eV << ' '
outputStream_ << static_cast<int>(get_PDG(particle.getPID())) << ' ' << energy / 1_eV
<< ' '
<< (trajectory.getPosition(1) - plane_.getCenter()).getNorm() / 1_m
<< std::endl;
......@@ -52,7 +52,7 @@ namespace corsika::observation_plane {
}
}
corsika::LengthType ObservationPlane::getMaxStepLength(
LengthType ObservationPlane::getMaxStepLength(
corsika::setup::Stack::particle_type const&,
corsika::setup::Trajectory const& trajectory) {
......@@ -85,4 +85,4 @@ namespace corsika::observation_plane {
count_ground_ = 0;
}
} // namespace corsika::observation_plane
} // namespace corsika
/*
* (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/framework/geometry/FourVector.hpp>
#include <corsika/modules/OnShellCheck.hpp>
namespace corsika {
OnShellCheck::OnShellCheck(const double vMassTolerance, const double vEnergyTolerance,
const bool vError)
: mass_tolerance_(vMassTolerance)
, energy_tolerance_(vEnergyTolerance)
, throw_error_(vError) {
std::cout << "OnShellCheck: mass tolerance is set to " << mass_tolerance_ * 100 << "%"
<< std::endl
<< " energy tolerance is set to " << energy_tolerance_ * 100
<< "%" << std::endl;
}
OnShellCheck::~OnShellCheck() {
std::cout << "OnShellCheck: summary" << std::endl
<< " particles shifted: " << int(count_) << std::endl;
if (count_)
std::cout << " average energy shift (%): " << average_shift_ / count_ * 100.
<< std::endl
<< " max. energy shift (%): " << max_shift_ * 100. << std::endl;
}
template <typename TView>
void OnShellCheck::doSecondaries(TView& vS) {
for (auto& p : vS) {
auto const pid = p.getPID();
if (!is_hadron(pid) || is_nucleus(pid)) continue;
auto const e_original = p.getEnergy();
auto const p_original = p.getMomentum();
auto const Plab = FourVector(e_original, p_original);
auto const m_kinetic = Plab.getNorm();
auto const m_corsika = get_mass(pid);
auto const m_err_abs = abs(m_kinetic - m_corsika);
if (m_err_abs >= mass_tolerance_ * m_corsika) {
const HEPEnergyType e_shifted =
sqrt(p_original.getSquaredNorm() + m_corsika * m_corsika);
auto const e_shift_relative = (e_shifted / e_original - 1);
count_ = count_ + 1;
average_shift_ += abs(e_shift_relative);
if (abs(e_shift_relative) > max_shift_) max_shift_ = abs(e_shift_relative);
std::cout << "OnShellCheck: shift particle mass for " << pid << std::endl
<< std::setw(40) << std::setfill(' ')
<< "corsika mass (GeV): " << m_corsika / 1_GeV << std::endl
<< std::setw(40) << std::setfill(' ')
<< "kinetic mass (GeV): " << m_kinetic / 1_GeV << std::endl
<< std::setw(40) << std::setfill(' ')
<< "m_kin-m_cor (GeV): " << m_err_abs / 1_GeV << std::endl
<< std::setw(40) << std::setfill(' ')
<< "mass tolerance (GeV): " << (m_corsika * mass_tolerance_) / 1_GeV
<< std::endl;
/*
For now we warn if the necessary shift is larger than 1%.
we could promote this to an error.
*/
if (abs(e_shift_relative) > energy_tolerance_) {
std::cout << "OnShellCheck: warning! shifted particle energy by "
<< e_shift_relative * 100 << " %" << std::endl;
if (throw_error_)
throw std::runtime_error(
"OnShellCheck: error! shifted energy by large amount!");
}
// reset energy
p.setEnergy(e_shifted);
} else
std::cout << "OnShellCheck: particle mass for " << pid << " OK" << std::endl;
}
}
} // namespace corsika
......@@ -10,7 +10,7 @@
#include <corsika/modules/ParticleCut.hpp>
namespace corsika::particle_cut {
namespace corsika {
ParticleCut::ParticleCut(const HEPEnergyType eCut, bool em, bool inv)
: energy_cut_(eCut)
......@@ -130,4 +130,4 @@ namespace corsika::particle_cut {
energy_ = 0_GeV;
}
} // namespace corsika::particle_cut
} // namespace corsika
......@@ -21,7 +21,7 @@
#include <iostream>
#include <limits>
namespace corsika::stack_inspector {
namespace corsika {
template <typename TStack>