IAP GITLAB

Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • AirShowerPhysics/corsika
  • rulrich/corsika
  • AAAlvesJr/corsika
  • Andre/corsika
  • arrabito/corsika
  • Nikos/corsika
  • olheiser73/corsika
  • AirShowerPhysics/papers/corsika
  • pranav/corsika
9 results
Show changes
Showing
with 6 additions and 12678 deletions
/*
* (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 <chrono>
#include <utility>
namespace corsika::analytics {
/// Wraps and measures the runtime of a single function type object
/**
*
* @tparam TFunc funtion pointer that should be wrapped
* @tparam TClock type of the clock that should be used for measurements
* @tparam TDuration type of std::duration to measure the elapsed time
*/
template <typename TFunc, typename TClock = std::chrono::high_resolution_clock,
typename TDuration = std::chrono::microseconds>
class FunctionTimer {
private:
typename TClock::time_point start_;
TDuration timeDiff_;
TFunc function_;
public:
/// Constructs the wrapper with the given functionpointer
FunctionTimer(TFunc f)
: function_(f) {}
template <typename... TArgs>
auto operator()(TArgs&&... args) -> std::invoke_result_t<TFunc, TArgs...> {
start_ = TClock::now();
auto tmp = function_(std::forward<TArgs>(args)...);
timeDiff_ = std::chrono::duration_cast<TDuration>(TClock::now() - start_);
return tmp;
}
inline TDuration getTime() const { return timeDiff_; }
};
} // namespace corsika::analytics
/*
* (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.
*/
/**
* @File Logging.h
*
* CORSIKA8 logging utilities.
*
* See testLogging.cc for a complete set of examples for
* how the logging functions should be used.
*/
#pragma once
// Configure some behaviour of sdlog.
// This must be done before spdlog is included.
// use the coarse system clock. This is *much* faster
// but introduces a real time error of O(10 ms).
#define SPDLOG_CLOCK_COARSE
// do not create a default logger (we provide our own "corsika" logger)
#define SPDLOG_DISABLE_DEFAULT_LOGGER
// use __PRETTY_FUNCTION__ instead of __FUNCTION__ where
// printing function names in trace statements. This is much
// nicer than __FUNCTION__ under GCC/clang.
#define SPDLOG_FUNCTION __PRETTY_FUNCTION__
// if this is a Debug build, include debug messages in objects
#ifdef DEBUG
// trace is the highest level of logging (ALL messages will be printed)
#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_TRACE
#else // otherwise, remove everything but "critical" messages
#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_CRITICAL
#endif
#include <spdlog/fmt/ostr.h> // will output whenerver a streaming operator is found
#include <spdlog/sinks/stdout_color_sinks.h>
#include <spdlog/spdlog.h>
namespace corsika::logging {
// bring spdlog into the corsika::logging namespace
using namespace spdlog;
/*
* The default pattern for CORSIKA8 loggers.
*/
const std::string default_pattern{"[%n:%^%-8l%$] %v"};
/**
* Create a new C8-style logger.
*
* Use this if you are explicitly (and can guarantee) that you
* are creating a logger for the first time. It is recommended
* that for regular usage, the `GetLogger` function is used instead
* as that will also create the logger if it has not yet been created.
*
* Calling `CreateLogger` twice to create the same logger will
* result in an spdlog duplicate exception.
*
* @param name The unique name of the logger.
* @param defaultlog If True, set this as the default logger.
* @returns The constructed and formatted logger.
*/
inline auto CreateLogger(std::string const& name, bool const defaultlog = false) {
// create the logger
// this is currently a colorized multi-threading safe logger
auto logger = spdlog::stdout_color_mt(name);
// set the default C8 format
logger->set_pattern(default_pattern);
// if defaultlog is True, we set this as the default spdlog logger.
if (defaultlog) { spdlog::set_default_logger(logger); }
// and return the logger
return logger;
}
/**
* Get a smart pointer to an existing logger.
*
* This should be the default method for code to obtain a
* logger. If the logger *does not* exist, it is *created* and
* returned to the caller.
*
* This should be preferred over `CreateLogger`.
*
* @param name The name of the logger to get.
* @param defaultlog If True, make this the default logger.
* @returns The constructed and formatted logger.
*/
inline auto GetLogger(std::string const& name, bool const defaultlog = false) {
// attempt to get the logger from the registry
auto logger = spdlog::get(name);
// weg found the logger, so just return it
if (logger) {
return logger;
} else { // logger was not found so create it
return CreateLogger(name, defaultlog);
}
}
/**
* The default "corsika" logger.
*/
inline auto corsika = GetLogger("corsika", true);
/**
* Set the default log level for all *newly* created loggers.
*
* @param name The minimum log level required to print.
*
*/
inline auto SetDefaultLevel(level::level_enum const minlevel) -> void {
spdlog::set_level(minlevel);
}
/**
* Set the log level for the "corsika" logger.
*
* @param name The minimum log level required to print.
*
*/
inline auto SetLevel(level::level_enum const minlevel) -> void {
corsika->set_level(minlevel);
}
/**
* Set the log level for a specific logger.
*
* @param logger The logger to set the level of.
* @param name The minimum log level required to print.
*
*/
template <typename TLogger>
inline auto SetLevel(TLogger& logger, level::level_enum const minlevel) -> void {
logger->set_level(minlevel);
}
/**
* Add the source (filename, line no) info to the logger.
*
* @param logger The logger to set the level of.
*
*/
template <typename TLogger>
inline auto AddSourceInfo(TLogger& logger) -> void {
logger->set_pattern("[%n:%^%-8l%$(%s:%!:%#)] %v");
}
/**
* Reset the logging pattern to the default.
*
* @param logger The logger to set the level of.
*
*/
template <typename TLogger>
inline auto ResetPattern(TLogger& logger) -> void {
logger->set_pattern(default_pattern);
}
// define our macro-style loggers
#define C8LOG_TRACE SPDLOG_TRACE
#define C8LOG_DEBUG SPDLOG_DEBUG
#define C8LOG_INFO SPDLOG_INFO
#define C8LOG_WARN SPDLOG_WARN
#define C8LOG_ERROR SPDLOG_ERROR
#define C8LOG_CRITICAL SPDLOG_CRITICAL
// and the specific logger versions
#define C8LOG_LOGGER_TRACE SPDLOG_LOGGER_TRACE
#define C8LOG_LOGGER_DEBUG SPDLOG_LOGGER_DEBUG
#define C8LOG_LOGGER_INFO SPDLOG_LOGGER_INFO
#define C8LOG_LOGGER_WARN SPDLOG_LOGGER_WARN
#define C8LOG_LOGGER_ERROR SPDLOG_LOGGER_ERROR
#define C8LOG_LOGGER_CRITICAL SPDLOG_LOGGER_CRITICAL
} // namespace corsika::logging
This diff is collapsed.
#!/usr/bin/env python3
import sys, math, itertools, re, csv, pprint
import xml.etree.ElementTree as ET
from collections import OrderedDict
import pickle
import io
GeVfm = 0.19732696312541853
c_speed_of_light = 29.9792458e10 # mm / s
# for nuclear masses
mneutron = 0.9395654133 # GeV
mproton = 0.9382720813 # GeV
##############################################################
#
# reading xml input data, return line by line particle data
#
def parsePythia(filename):
tree = ET.parse(filename)
root = tree.getroot()
for particle in root.iter("particle"):
name = particle.attrib["name"]
antiName = "Unknown"
if ("antiName" in particle.attrib):
antiName = particle.attrib["antiName"]
pdg_id = int(particle.attrib["id"])
mass = float(particle.attrib["m0"]) # GeV
electric_charge = int(particle.attrib["chargeType"]) # in units of e/3
ctau = 0.
if pdg_id in (11, 12, 14, 16, 22, 2212): # these are the stable particles !
ctau = float('Inf')
elif 'tau0' in particle.attrib:
ctau = float(particle.attrib['tau0']) # mm / c
elif 'mWidth' in particle.attrib:
ctau = GeVfm / float(particle.attrib['mWidth']) * 1e-15 * 1000.0 # mm / s
elif pdg_id in (0, 423, 433, 4312, 4322, 5112, 5222): # those are certainly not stable....
ctau = 0.
else:
print ("missing lifetime: " + str(pdg_id) + " " + str(name))
sys.exit(1)
yield (pdg_id, name, mass, electric_charge, antiName, ctau/c_speed_of_light)
# TODO: read decay channels from child elements
if "antiName" in particle.attrib:
yield (-pdg_id, antiName, mass, -electric_charge, name, ctau/c_speed_of_light)
##############################################################
#
# reading xml input data, return line by line particle data
#
def parseNuclei(filename):
tree = ET.parse(filename)
root = tree.getroot()
for particle in root.iter("particle"):
name = particle.attrib["name"]
antiName = "Unknown"
if ("antiName" in particle.attrib):
antiName = particle.attrib["antiName"]
pdg_id = int(particle.attrib["id"])
A = int(particle.attrib["A"])
Z = int(particle.attrib["Z"])
# mass in GeV
if ("mass" in particle.attrib):
mass = particle.attrib["mass"]
else:
mass = (A-Z)*mneutron + Z*mproton
electric_charge = Z*3 # in units of e/3
ctau = float('Inf')
yield (pdg_id, name, mass, electric_charge, antiName, ctau/c_speed_of_light, A, Z)
##############################################################
#
# returns dict with particle codes and class names
#
def class_names(filename):
tree = ET.parse(filename)
root = tree.getroot()
map = {}
for particle in root.iter("particle"):
name = particle.attrib["classname"]
pdg_id = int(particle.attrib["pdgID"])
map[pdg_id] = name
return map
##############################################################
#
# Automatically produce a string qualifying as C++ class name
#
# This function produces names of type "DeltaPlusPlus"
#
def c_identifier_camel(name):
orig = name
name = name[0].upper() + name[1:].lower() # all lower case
for c in "() /": # replace funny characters
name = name.replace(c, "_")
name = name.replace("bar", "Bar")
name = name.replace("*", "Star")
name = name.replace("'", "Prime")
name = name.replace("+", "Plus")
name = name.replace("-", "Minus")
# move "Bar" to end of name
ibar = name.find('Bar')
if ibar > 0 and ibar < len(name)-3:
name = name[:ibar] + name[ibar+3:] + 'Bar'
# cleanup "_"s
while True:
tmp = name.replace("__", "_")
if tmp == name:
break
else:
name = tmp
name.strip("_")
# remove all "_", if this does not by accident concatenate two numbers
istart = 0
while True:
i = name.find('_', istart)
if i < 1 or i > len(name)-1:
break
istart = i
if name[i-1].isdigit() and name[i+1].isdigit():
# there is a number on both sides
break
name = name[:i] + name[i+1:]
# and last, for example: make NuE out of Nue
if name[i-1].islower() and name[i].islower():
if i < len(name)-1:
name = name[:i] + name[i].upper() + name[i+1:]
else:
name = name[:i] + name[i].upper()
# check if name is valid C++ identifier
pattern = re.compile(r'^[a-zA-Z_][a-zA-Z_0-9]*$')
if pattern.match(name):
return name
else:
raise Exception("could not generate C identifier for '{:s}': result '{:s}'".format(orig, name))
##########################################################
#
# returns dict containing all data from pythia-xml input
#
def read_pythia_db(filename, particle_db, classnames):
counter = itertools.count(len(particle_db))
for (pdg, name, mass, electric_charge, antiName, lifetime) in parsePythia(filename):
c_id = "Unknown"
if pdg in classnames:
c_id = classnames[pdg]
else:
c_id = c_identifier_camel(name) # the camel case names
hadron =abs(pdg) > 100
if c_id in particle_db.keys():
raise RuntimeError("particle '{:s}' already known (new PDG id {:d}, stored PDG id: {:d})".format(c_id, pdg, particle_db[c_id]['pdg']))
particle_db[c_id] = {
"name" : name,
"antiName" : antiName,
"pdg" : pdg,
"mass" : mass, # in GeV
"electric_charge" : electric_charge, # in e/3
"lifetime" : lifetime,
"ngc_code" : next(counter),
"isNucleus" : False,
"isHadron" : hadron,
}
return particle_db
##########################################################
#
# returns dict containing all data from pythia-xml input
#
def read_nuclei_db(filename, particle_db, classnames):
counter = itertools.count(len(particle_db))
for (pdg, name, mass, electric_charge, antiName, lifetime, A, Z) in parseNuclei(filename):
c_id = "Unknown"
if pdg in classnames:
c_id = classnames[pdg]
else:
c_id = c_identifier_camel(name)
particle_db[c_id] = {
"name" : name,
"antiName" : antiName,
"pdg" : pdg,
"mass" : mass, # in GeV
"electric_charge" : electric_charge, # in e/3
"lifetime" : lifetime,
"ngc_code" : next(counter),
"A" : A,
"Z" : Z,
"isNucleus" : True,
"isHadron" : True,
}
return particle_db
###############################################################
#
# build conversion table PDG -> ngc
#
def gen_conversion_PDG_ngc(particle_db):
# todo: find a optimum value, think about cache miss with array vs lookup time with map
P_MAX = 500 # the maximum PDG code that is still filled into the table
conversionDict = dict()
conversionTable = [None] * (2*P_MAX + 1)
for cId, p in particle_db.items():
pdg = p['pdg']
if abs(pdg) < P_MAX:
if conversionTable[pdg + P_MAX]:
raise Exception("table entry already occupied")
else:
conversionTable[pdg + P_MAX] = cId
else:
if pdg in conversionDict.keys():
raise Exception(f"map entry {pdg} already occupied")
else:
conversionDict[pdg] = cId
output = io.StringIO()
def oprint(*args, **kwargs):
print(*args, **kwargs, file=output)
oprint(f"static std::array<Code, {len(conversionTable)}> constexpr conversionArray {{")
for ngc in conversionTable:
oprint(" Code::{0},".format(ngc if ngc else "Unknown"))
oprint("};")
oprint()
oprint("static std::map<PDGCode, Code> const conversionMap {")
for ngc in conversionDict.values():
oprint(f" {{PDGCode::{ngc}, Code::{ngc}}},")
oprint("};")
oprint()
return output.getvalue()
###############################################################
#
# return string with enum of all internal particle codes
#
def gen_internal_enum(particle_db):
string = ("enum class Code : CodeIntType {\n"
" FirstParticle = 1, // if you want to loop over particles, you want to start with \"1\" \n") # identifier for eventual loops...
for k in filter(lambda k: "ngc_code" in particle_db[k], particle_db):
last_ngc_id = particle_db[k]['ngc_code']
string += " {key:s} = {code:d},\n".format(key = k, code = last_ngc_id)
string += (" LastParticle = {:d},\n" # identifier for eventual loops...
"}};").format(last_ngc_id + 1)
if last_ngc_id > 0x7fff: # does not fit into int16_t
raise Exception("Integer overflow in internal particle code definition prevented!")
return string
###############################################################
#
# return string with enum of all PDG particle codes
#
def gen_pdg_enum(particle_db):
string = "enum class PDGCode : PDGCodeType {\n"
for cId in particle_db:
pdgCode = particle_db[cId]['pdg']
string += " {key:s} = {code:d},\n".format(key = cId, code = pdgCode)
string += " };\n"
return string
###############################################################
#
# return string with all data arrays
#
def gen_properties(particle_db):
# number of particles, size of tables
string = "static constexpr std::size_t size = {size:d};\n".format(size = len(particle_db))
string += "\n"
# all particle initializer_list
string += "constexpr std::initializer_list<Code> all_particles = {"
for k in particle_db:
string += " Code::{name:s},\n".format(name = k)
string += "};\n"
string += "\n"
# particle masses table
string += "static constexpr std::array<corsika::units::si::HEPMassType const, size> masses = {\n"
for p in particle_db.values():
string += " {mass:e} * 1e9 * corsika::units::si::electronvolt, // {name:s}\n".format(mass = p['mass'], name = p['name'])
string += "};\n\n"
# PDG code table
string += "static constexpr std::array<PDGCode, size> pdg_codes = {\n"
for p in particle_db.keys():
string += f" PDGCode::{p},\n"
string += "};\n"
# name string table
string += "static const std::array<std::string const, size> names = {\n"
for p in particle_db.values():
string += " \"{name:s}\",\n".format(name = p['name'])
string += "};\n"
# electric charges table
string += "static constexpr std::array<int16_t, size> electric_charges = {\n"
for p in particle_db.values():
string += " {charge:d},\n".format(charge = p['electric_charge'])
string += "};\n"
# anti-particle table
# string += "static constexpr std::array<size, size> anti_particle = {\n"
# for p in particle_db.values():
# string += " {anti:d},\n".format(charge = p['anti_particle'])
# string += "};\n"
# lifetime
#string += "static constexpr std::array<corsika::units::si::TimeType const, size> lifetime = {\n"
string += "static constexpr std::array<double const, size> lifetime = {\n"
for p in particle_db.values():
if p['lifetime'] == float("Inf") :
string += " std::numeric_limits<double>::infinity(), \n" # * corsika::units::si::second, \n"
else :
string += " {tau:e}, \n".format(tau = p['lifetime'])
#string += " {tau:e} * corsika::units::si::second, \n".format(tau = p['lifetime'])
string += "};\n"
# is Hadron flag
string += "static constexpr std::array<bool, size> isHadron = {\n"
for p in particle_db.values():
value = 'false'
if p['isHadron']:
value = 'true'
string += " {val},\n".format(val = value)
string += "};\n"
### nuclear data ###
# nucleus mass number A
string += "static constexpr std::array<int16_t, size> nucleusA = {\n"
for p in particle_db.values():
A = p.get('A', 0)
string += " {val},\n".format(val = A)
string += "};\n"
# nucleus charge number Z
string += "static constexpr std::array<int16_t, size> nucleusZ = {\n"
for p in particle_db.values():
Z = p.get('Z', 0)
string += " {val},\n".format(val = Z)
string += "};\n"
return string
###############################################################
#
# return string with a list of classes for all particles
#
def gen_classes(particle_db):
string = "// list of C++ classes to access particle properties"
for cname in particle_db:
antiP = 'Unknown'
for cname_anti in particle_db:
if (particle_db[cname_anti]['name'] == particle_db[cname]['antiName']):
antiP = cname_anti
break
string += "\n"
string += "/** @class " + cname + "\n\n"
string += " * Particle properties are taken from the PYTHIA8 ParticleData.xml file:<br>\n"
string += " * - pdg=" + str(particle_db[cname]['pdg']) +"\n"
string += " * - mass=" + str(particle_db[cname]['mass']) + " GeV \n"
string += " * - charge= " + str(particle_db[cname]['electric_charge']/3) + " \n"
string += " * - name=" + str(cname) + "\n"
string += " * - anti=" + str(antiP) + "\n"
if (particle_db[cname]['isNucleus']):
string += " * - nuclear A=" + str(particle_db[cname]['A']) + "\n"
string += " * - nuclear Z=" + str(particle_db[cname]['Z']) + "\n"
string += "*/\n\n"
string += "class " + cname + " {\n"
string += " public:\n"
string += " static constexpr Code GetCode() { return Type; }\n"
string += " static constexpr corsika::units::si::HEPMassType GetMass() { return corsika::particles::GetMass(Type); }\n"
string += " static constexpr corsika::units::si::ElectricChargeType GetCharge() { return corsika::particles::GetCharge(Type); }\n"
string += " static constexpr int16_t GetChargeNumber() { return corsika::particles::GetChargeNumber(Type); }\n"
string += " static std::string const& GetName() { return corsika::particles::GetName(Type); }\n"
string += " static constexpr Code GetAntiParticle() { return AntiType; }\n"
string += " static constexpr bool IsNucleus() { return corsika::particles::IsNucleus(Type); }\n"
string += " static constexpr int16_t GetNucleusA() { return corsika::particles::GetNucleusA(Type); }\n"
string += " static constexpr int16_t GetNucleusZ() { return corsika::particles::GetNucleusZ(Type); }\n"
string += " static constexpr Code Type = Code::" + cname + ";\n"
string += " static constexpr Code AntiType = Code::" + antiP + ";\n"
string += " private:\n"
string += " static constexpr CodeIntType TypeIndex = static_cast<CodeIntType const>(Type);\n"
string += "};\n"
return string
###############################################################
#
#
def inc_start():
string = ('// generated by pdxml_reader.py\n'
'// MANUAL EDITS ON OWN RISK. THEY WILL BE OVERWRITTEN. \n'
'namespace corsika::particles {\n')
return string
###############################################################
#
#
def detail_start():
string = ('namespace detail {\n\n')
return string
###############################################################
#
#
def detail_end():
string = "\n} // end namespace detail\n"
return string
###############################################################
#
#
def inc_end():
string = "\n} // end namespace corsika::particles\n"
return string
###################################################################
#
# Serialize particle_db into file
#
def serialize_particle_db(particle_db, file):
pickle.dump(particle_db, file)
###################################################################
#
# Main function
#
if __name__ == "__main__":
if len(sys.argv) != 4:
print("usage: {:s} <Pythia8.xml> <Nuclei.xml> <ClassNames.xml>".format(sys.argv[0]), file=sys.stderr)
sys.exit(1)
print("\n pdxml_reader.py: automatically produce particle properties from input files\n")
names = class_names(sys.argv[3])
particle_db = OrderedDict()
read_pythia_db(sys.argv[1], particle_db, names)
read_nuclei_db(sys.argv[2], particle_db, names)
with open("GeneratedParticleProperties.inc", "w") as f:
print(inc_start(), file=f)
print(gen_internal_enum(particle_db), file=f)
print(gen_pdg_enum(particle_db), file=f)
print(detail_start(), file=f)
print(gen_properties(particle_db), file=f)
print(gen_conversion_PDG_ngc(particle_db), file=f)
print(detail_end(), file=f)
print(gen_classes(particle_db), file=f)
print(inc_end(), file=f)
with open("particle_db.pkl", "wb") as f:
serialize_particle_db(particle_db, f)
/*
* (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/process/BaseProcess.h>
namespace corsika::process {
class NullModel : public corsika::process::BaseProcess<NullModel> {
public:
NullModel() = default;
~NullModel() = default;
};
} // namespace corsika::process
/*
* (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
namespace corsika::process {
/**
* A traits marker to track which BaseProcess is also a ProcessSequence
**/
template <typename TClass>
struct is_process_sequence : std::false_type {};
template <typename TClass>
bool constexpr is_process_sequence_v = is_process_sequence<TClass>::value;
/**
* A traits marker to identiy a BaseProcess that is also SwitchProcessesSequence
**/
template <typename TClass>
struct is_switch_process_sequence : std::false_type {};
template <typename TClass>
bool constexpr is_switch_process_sequence_v = is_switch_process_sequence<TClass>::value;
/**
* A traits marker to identify ProcessSequence that contain a StackProcess
**/
template <typename TClass>
struct contains_stack_process : std::false_type {};
template <typename TClass>
bool constexpr contains_stack_process_v = contains_stack_process<TClass>::value;
} // namespace corsika::process
/*
* (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/process/BaseProcess.h>
#include <corsika/process/ProcessTraits.h>
#include <corsika/process/BoundaryCrossingProcess.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/process/DecayProcess.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/ProcessReturn.h>
#include <corsika/process/SecondariesProcess.h>
#include <corsika/process/StackProcess.h>
#include <corsika/units/PhysicalUnits.h>
#include <cmath>
#include <limits>
#include <type_traits>
namespace corsika::process {
// enum for the process switch selection: identify if First or
// Second process branch should be used.
enum class SwitchResult { First, Second };
/**
\class SwitchProcessSequence
A compile time static list of processes that uses an internal
TSelect class to switch between different versions of processes
(or process sequence).
TProcess1 and TProcess2 must be derived from BaseProcess and are
both references if possible (lvalue), otherwise (rvalue) they are
just classes. This allows us to handle both, rvalue as well as
lvalue Processes in the SwitchProcessSequence.
TSelect has to implement a `operator()(const Particle&)` and has to
return either SwitchResult::First or SwitchResult::Second. Note:
TSelect may absolutely also use random numbers to sample between
its results. This can be used to achieve arbitrarily smooth
transition or mixtures of processes.
Warning: do not put StackProcess into a SwitchProcessSequence
since this makes no sense. The StackProcess acts on an entire
particle stack and not on indiviidual particles.
\comment See also class ProcessSequence
**/
template <typename TProcess1, typename TProcess2, typename TSelect>
class SwitchProcessSequence
: public BaseProcess<SwitchProcessSequence<TProcess1, TProcess2, TSelect>> {
using TProcess1type = typename std::decay_t<TProcess1>;
using TProcess2type = typename std::decay_t<TProcess2>;
static bool constexpr t1ProcSeq = is_process_sequence_v<TProcess1type>;
static bool constexpr t2ProcSeq = is_process_sequence_v<TProcess2type>;
// make sure only BaseProcess types TProcess1/2 are passed
static_assert(std::is_base_of_v<BaseProcess<TProcess1type>, TProcess1type>,
"can only use process derived from BaseProcess in "
"SwitchProcessSequence, for Process 1");
static_assert(std::is_base_of_v<BaseProcess<TProcess2type>, TProcess2type>,
"can only use process derived from BaseProcess in "
"SwitchProcessSequence, for Process 2");
// make sure none of TProcess1/2 is a StackProcess
static_assert(!std::is_base_of_v<StackProcess<TProcess1type>, TProcess1type>,
"cannot use StackProcess in SwitchProcessSequence, for Process 1");
static_assert(!std::is_base_of_v<StackProcess<TProcess2type>, TProcess2type>,
"cannot use StackProcess in SwitchProcessSequence, for Process 2");
// if TProcess1/2 are already ProcessSequences, make sure they do not contain
// any StackProcess
static_assert(!contains_stack_process_v<TProcess1type>,
"cannot use StackProcess in SwitchProcessSequence, remove from "
"ProcessSequence 1");
static_assert(!contains_stack_process_v<TProcess2type>,
"cannot use StackProcess in SwitchProcessSequence, remove from "
"ProcessSequence 2");
TSelect select_; // this is a reference, if possible
TProcess1 A_; // this is a reference, if possible
TProcess2 B_; // this is a reference, if possible
public:
SwitchProcessSequence(TProcess1 in_A, TProcess2 in_B, TSelect sel)
: select_(sel)
, A_(in_A)
, B_(in_B) {}
template <typename TParticle, typename TVTNType>
EProcessReturn DoBoundaryCrossing(TParticle& particle, TVTNType const& from,
TVTNType const& to) {
switch (select_(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess1type>,
TProcess1type> ||
t1ProcSeq) {
return A_.DoBoundaryCrossing(particle, from, to);
}
break;
}
case SwitchResult::Second: {
if constexpr (std::is_base_of_v<BoundaryCrossingProcess<TProcess2type>,
TProcess2type> ||
t2ProcSeq) {
return B_.DoBoundaryCrossing(particle, from, to);
}
break;
}
}
return EProcessReturn::eOk;
}
template <typename TParticle, typename TTrack>
inline EProcessReturn DoContinuous(TParticle& particle, TTrack& vT) {
switch (select_(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>,
TProcess1type> ||
t1ProcSeq) {
return A_.DoContinuous(particle, vT);
}
break;
}
case SwitchResult::Second: {
if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>,
TProcess2type> ||
t2ProcSeq) {
return B_.DoContinuous(particle, vT);
}
break;
}
}
return EProcessReturn::eOk;
}
template <typename TSecondaries>
inline void DoSecondaries(TSecondaries& vS) {
const auto& particle = vS.parent();
switch (select_(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<SecondariesProcess<TProcess1type>,
TProcess1type> ||
t1ProcSeq) {
A_.DoSecondaries(vS);
}
break;
}
case SwitchResult::Second: {
if constexpr (std::is_base_of_v<SecondariesProcess<TProcess2type>,
TProcess2type> ||
t2ProcSeq) {
B_.DoSecondaries(vS);
}
break;
}
}
}
template <typename TParticle, typename TTrack>
inline corsika::units::si::LengthType MaxStepLength(TParticle& particle,
TTrack& vTrack) {
switch (select_(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<ContinuousProcess<TProcess1type>,
TProcess1type> ||
t1ProcSeq) {
return A_.MaxStepLength(particle, vTrack);
}
break;
}
case SwitchResult::Second: {
if constexpr (std::is_base_of_v<ContinuousProcess<TProcess2type>,
TProcess2type> ||
t2ProcSeq) {
return B_.MaxStepLength(particle, vTrack);
}
break;
}
}
// if no other process in the sequence implements it
return std::numeric_limits<double>::infinity() * corsika::units::si::meter;
}
template <typename TParticle>
inline corsika::units::si::GrammageType GetInteractionLength(TParticle&& particle) {
return 1. / GetInverseInteractionLength(particle);
}
template <typename TParticle>
inline corsika::units::si::InverseGrammageType GetInverseInteractionLength(
TParticle&& particle) {
using namespace corsika::units::si;
switch (select_(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>,
TProcess1type> ||
t1ProcSeq) {
return A_.GetInverseInteractionLength(particle);
}
break;
}
case SwitchResult::Second: {
if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>,
TProcess2type> ||
t2ProcSeq) {
return B_.GetInverseInteractionLength(particle);
}
break;
}
}
return 0 * meter * meter / gram; // default value
}
template <typename TSecondaryView>
inline EProcessReturn SelectInteraction(
TSecondaryView& view,
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_select,
[[maybe_unused]] corsika::units::si::InverseGrammageType lambda_inv_sum =
corsika::units::si::InverseGrammageType::zero()) {
switch (select_(view.parent())) {
case SwitchResult::First: {
if constexpr (t1ProcSeq) {
// if A_ is a process sequence --> check inside
const EProcessReturn ret =
A_.SelectInteraction(view, lambda_inv_select, lambda_inv_sum);
// if A_ did succeed, stop routine. Not checking other static branch B_.
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of_v<InteractionProcess<TProcess1type>,
TProcess1type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += A_.GetInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
if (lambda_inv_select < lambda_inv_sum) {
A_.DoInteraction(view);
return EProcessReturn::eInteracted;
}
} // end branch A_
break;
}
case SwitchResult::Second: {
if constexpr (t2ProcSeq) {
// if B_ is a process sequence --> check inside
return B_.SelectInteraction(view, lambda_inv_select, lambda_inv_sum);
} else if constexpr (std::is_base_of_v<InteractionProcess<TProcess2type>,
TProcess2type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += B_.GetInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
if (lambda_inv_select < lambda_inv_sum) {
B_.DoInteraction(view);
return EProcessReturn::eInteracted;
}
} // end branch B_
break;
}
}
return EProcessReturn::eOk;
}
template <typename TParticle>
inline corsika::units::si::TimeType GetLifetime(TParticle&& particle) {
return 1. / GetInverseLifetime(particle);
}
template <typename TParticle>
inline corsika::units::si::InverseTimeType GetInverseLifetime(TParticle&& particle) {
using namespace corsika::units::si;
switch (select_(particle)) {
case SwitchResult::First: {
if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>, TProcess1type> ||
t1ProcSeq) {
return A_.GetInverseLifetime(particle);
}
break;
}
case SwitchResult::Second: {
if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>, TProcess2type> ||
t2ProcSeq) {
return B_.GetInverseLifetime(particle);
}
break;
}
}
return 0 / second; // default value
}
// select decay process
template <typename TSecondaryView>
inline EProcessReturn SelectDecay(
TSecondaryView& view,
[[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_select,
[[maybe_unused]] corsika::units::si::InverseTimeType decay_inv_sum =
corsika::units::si::InverseTimeType::zero()) {
switch (select_(view.parent())) {
case SwitchResult::First: {
if constexpr (t1ProcSeq) {
// if A_ is a process sequence --> check inside
const EProcessReturn ret =
A_.SelectDecay(view, decay_inv_select, decay_inv_sum);
// if A_ did succeed, stop routine here (not checking other static branch B_)
if (ret != EProcessReturn::eOk) { return ret; }
} else if constexpr (std::is_base_of_v<DecayProcess<TProcess1type>,
TProcess1type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += A_.GetInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
if (decay_inv_select < decay_inv_sum) {
// more pedagogical: rndm_select < decay_inv_sum / decay_inv_tot
A_.DoDecay(view);
return EProcessReturn::eDecayed;
}
} // end branch A_
break;
}
case SwitchResult::Second: {
if constexpr (t2ProcSeq) {
// if B_ is a process sequence --> check inside
return B_.SelectDecay(view, decay_inv_select, decay_inv_sum);
} else if constexpr (std::is_base_of_v<DecayProcess<TProcess2type>,
TProcess2type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += B_.GetInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
if (decay_inv_select < decay_inv_sum) {
B_.DoDecay(view);
return EProcessReturn::eDecayed;
}
} // end branch B_
break;
}
}
return EProcessReturn::eOk;
}
};
// the method `select(proc1,proc1,selector)` assembles many
// BaseProcesses, and ProcessSequences into a SwitchProcessSequence,
// all combinatorics must be allowed, this is why we define a macro
// to define all combinations here:
// Both, Processes1 and Processes2, must derive from BaseProcesses
template <typename TProcess1, typename TProcess2, typename TSelect>
inline typename std::enable_if_t<
std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess1>>,
typename std::decay_t<TProcess1>> &&
std::is_base_of_v<BaseProcess<typename std::decay_t<TProcess2>>,
typename std::decay_t<TProcess2>>,
SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
select(TProcess1&& vA, TProcess2&& vB, TSelect selector) {
return SwitchProcessSequence<TProcess1, TProcess2, TSelect>(vA, vB, selector);
}
/// traits marker to identify objectas ProcessSequence
template <typename TProcess1, typename TProcess2, typename TSelect>
struct is_process_sequence<
corsika::process::SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
: std::true_type {};
/// traits marker to identify objectas SwitchProcessSequence
template <typename TProcess1, typename TProcess2, typename TSelect>
struct is_switch_process_sequence<
corsika::process::SwitchProcessSequence<TProcess1, TProcess2, TSelect>>
: std::true_type {};
} // namespace corsika::process
/*
* (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.
*/
#include <catch2/catch.hpp>
#include <corsika/process/NullModel.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
using namespace corsika::units::si;
using namespace corsika::process;
using namespace corsika;
TEST_CASE("NullModel", "[processes]") {
SECTION("interface") { [[maybe_unused]] NullModel model; }
}
/*
* (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.
*/
#include <corsika/process/switch_process/SwitchProcess.h>
#include <corsika/stack/SecondaryView.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
#include <catch2/catch.hpp>
#include <algorithm>
#include <random>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units::si;
class TestStackData {
public:
// these functions are needed for the Stack interface
void Clear() { fData.clear(); }
unsigned int GetSize() const { return fData.size(); }
unsigned int GetCapacity() const { return fData.size(); }
void Copy(int i1, int i2) { fData[i2] = fData[i1]; }
void Swap(int i1, int i2) { std::swap(fData[i1], fData[i2]); }
// custom data access function
void SetData(unsigned int i, HEPEnergyType v) { fData[i] = v; }
HEPEnergyType GetData(const unsigned int i) const { return fData[i]; }
// these functions are also needed by the Stack interface
void IncrementSize() { fData.resize(fData.size() + 1); }
void DecrementSize() {
if (fData.size() > 0) { fData.pop_back(); }
}
// custom private data section
private:
std::vector<HEPEnergyType> fData;
};
/**
* From static_cast of a StackIterator over entries in the
* TestStackData class you get and object of type
* TestParticleInterface defined here
*
* It provides Get/Set methods to read and write data to the "Data"
* storage of TestStackData obtained via
* "StackIteratorInterface::GetStackData()", given the index of the
* iterator "StackIteratorInterface::GetIndex()"
*
*/
template <typename StackIteratorInterface>
class TestParticleInterface
: public corsika::stack::ParticleBase<StackIteratorInterface> {
public:
using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
/*
The SetParticleData methods are called for creating new entries
on the stack. You can specifiy various parametric versions to
perform this task:
*/
// default version for particle-creation from input data
void SetParticleData(const std::tuple<HEPEnergyType> v) { SetEnergy(std::get<0>(v)); }
void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/,
std::tuple<HEPEnergyType> v) {
SetEnergy(std::get<0>(v));
}
// here are the fundamental methods for access to TestStackData data
void SetEnergy(HEPEnergyType v) { GetStackData().SetData(GetIndex(), v); }
HEPEnergyType GetEnergy() const { return GetStackData().GetData(GetIndex()); }
};
using SimpleStack = corsika::stack::Stack<TestStackData, TestParticleInterface>;
// see issue 161
#if defined(__clang__)
using StackTestView = corsika::stack::SecondaryView<TestStackData, TestParticleInterface>;
#elif defined(__GNUC__) || defined(__GNUG__)
using StackTestView = corsika::stack::MakeView<SimpleStack>::type;
#endif
auto constexpr kgMSq = 1_kg / (1_m * 1_m);
template <int N>
struct DummyProcess : InteractionProcess<DummyProcess<N>> {
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const {
return N * kgMSq;
}
template <typename TSecondaries>
corsika::process::EProcessReturn DoInteraction(TSecondaries& vSec) {
// to figure out which process was selected in the end, we produce N
// secondaries for DummyProcess<N>
for (int i = 0; i < N; ++i) {
vSec.AddSecondary(std::tuple<HEPEnergyType>{vSec.GetEnergy() / N});
}
return EProcessReturn::eOk;
}
};
using DummyLowEnergyProcess = DummyProcess<1>;
using DummyHighEnergyProcess = DummyProcess<2>;
using DummyAdditionalProcess = DummyProcess<3>;
TEST_CASE("SwitchProcess from InteractionProcess") {
DummyLowEnergyProcess low;
DummyHighEnergyProcess high;
DummyAdditionalProcess proc;
switch_process::SwitchProcess switchProcess(low, high, 1_TeV);
auto seq = switchProcess << proc;
SimpleStack stack;
SECTION("low energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
auto p = stack.GetNextParticle();
// low energy process returns 1 kg/m²
SECTION("interaction length") {
CHECK(switchProcess.GetInteractionLength(p) / kgMSq == Approx(1));
CHECK(seq.GetInteractionLength(p) / kgMSq == Approx(3. / 4));
}
}
SECTION("high energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{4_TeV});
auto p = stack.GetNextParticle();
// high energy process returns 2 kg/m²
SECTION("interaction length") {
CHECK(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2));
CHECK(seq.GetInteractionLength(p) / kgMSq == Approx(6. / 5));
}
// high energy process creates 2 secondaries
SECTION("SelectInteraction") {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
InverseGrammageType invLambda = 0 / kgMSq;
switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
CHECK(view.getSize() == 2);
}
}
}
TEST_CASE("SwitchProcess from ProcessSequence") {
DummyProcess<1> innerA;
DummyProcess<2> innerB;
DummyProcess<3> outer;
DummyProcess<4> additional;
auto seq = innerA << innerB;
switch_process::SwitchProcess switchProcess(seq, outer, 1_TeV);
auto completeSeq = switchProcess << additional;
SimpleStack stack;
SECTION("low energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
auto p = stack.GetNextParticle();
SECTION("interaction length") {
CHECK(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2. / 3));
CHECK(completeSeq.GetInteractionLength(p) / kgMSq == Approx(4. / 7));
}
SECTION("SelectInteraction") {
std::vector<int> numberOfSecondaries;
for (int i = 0; i < 1000; ++i) {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
double r = i / 1000.;
InverseGrammageType invLambda = r * 7. / 4 / kgMSq;
InverseGrammageType accumulator = 0 / kgMSq;
completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
numberOfSecondaries.push_back(view.getSize());
}
auto const mean =
std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
numberOfSecondaries.size();
CHECK(mean == Approx(12. / 7.).margin(0.01));
}
}
SECTION("high energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{3.0_TeV});
auto p = stack.GetNextParticle();
SECTION("interaction length") {
CHECK(switchProcess.GetInteractionLength(p) / kgMSq == Approx(3));
CHECK(completeSeq.GetInteractionLength(p) / kgMSq == Approx(12. / 7.));
}
SECTION("SelectInteraction") {
std::vector<int> numberOfSecondaries;
for (int i = 0; i < 1000; ++i) {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
double r = i / 1000.;
InverseGrammageType invLambda = r * 7. / 12. / kgMSq;
InverseGrammageType accumulator = 0 / kgMSq;
completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
numberOfSecondaries.push_back(view.getSize());
}
auto const mean =
std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
numberOfSecondaries.size();
CHECK(mean == Approx(24. / 7.).margin(0.01));
}
}
}
/*
* (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.
*/
#include <corsika/process/switch_process/SwitchProcess.h>
#include <corsika/stack/SecondaryView.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
#include <catch2/catch.hpp>
#include <algorithm>
#include <random>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units::si;
class TestStackData {
public:
// these functions are needed for the Stack interface
void Clear() { fData.clear(); }
unsigned int GetSize() const { return fData.size(); }
unsigned int GetCapacity() const { return fData.size(); }
void Copy(int i1, int i2) { fData[i2] = fData[i1]; }
void Swap(int i1, int i2) { std::swap(fData[i1], fData[i2]); }
// custom data access function
void SetData(unsigned int i, HEPEnergyType v) { fData[i] = v; }
HEPEnergyType GetData(const unsigned int i) const { return fData[i]; }
// these functions are also needed by the Stack interface
void IncrementSize() { fData.resize(fData.size() + 1); }
void DecrementSize() {
if (fData.size() > 0) { fData.pop_back(); }
}
// custom private data section
private:
std::vector<HEPEnergyType> fData;
};
/**
* From static_cast of a StackIterator over entries in the
* TestStackData class you get and object of type
* TestParticleInterface defined here
*
* It provides Get/Set methods to read and write data to the "Data"
* storage of TestStackData obtained via
* "StackIteratorInterface::GetStackData()", given the index of the
* iterator "StackIteratorInterface::GetIndex()"
*
*/
template <typename StackIteratorInterface>
class TestParticleInterface
: public corsika::stack::ParticleBase<StackIteratorInterface> {
public:
using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
/*
The SetParticleData methods are called for creating new entries
on the stack. You can specifiy various parametric versions to
perform this task:
*/
// default version for particle-creation from input data
void SetParticleData(const std::tuple<HEPEnergyType> v) { SetEnergy(std::get<0>(v)); }
void SetParticleData(TestParticleInterface<StackIteratorInterface>& /*parent*/,
std::tuple<HEPEnergyType> v) {
SetEnergy(std::get<0>(v));
}
// here are the fundamental methods for access to TestStackData data
void SetEnergy(HEPEnergyType v) { GetStackData().SetData(GetIndex(), v); }
HEPEnergyType GetEnergy() const { return GetStackData().GetData(GetIndex()); }
};
using SimpleStack = corsika::stack::Stack<TestStackData, TestParticleInterface>;
// see issue 161
#if defined(__clang__)
using StackTestView = corsika::stack::SecondaryView<TestStackData, TestParticleInterface>;
#elif defined(__GNUC__) || defined(__GNUG__)
using StackTestView = corsika::stack::MakeView<SimpleStack>::type;
#endif
auto constexpr kgMSq = 1_kg / (1_m * 1_m);
template <int N>
struct DummyProcess : InteractionProcess<DummyProcess<N>> {
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const {
return N * kgMSq;
}
template <typename TSecondaries>
corsika::process::EProcessReturn DoInteraction(TSecondaries& vSec) {
// to figure out which process was selected in the end, we produce N
// secondaries for DummyProcess<N>
for (int i = 0; i < N; ++i) {
vSec.AddSecondary(std::tuple<HEPEnergyType>{vSec.GetEnergy() / N});
}
return EProcessReturn::eOk;
}
};
using DummyLowEnergyProcess = DummyProcess<1>;
using DummyHighEnergyProcess = DummyProcess<2>;
using DummyAdditionalProcess = DummyProcess<3>;
TEST_CASE("SwitchProcess from InteractionProcess") {
DummyLowEnergyProcess low;
DummyHighEnergyProcess high;
DummyAdditionalProcess proc;
switch_process::SwitchProcess switchProcess(low, high, 1_TeV);
auto seq = switchProcess << proc;
SimpleStack stack;
SECTION("low energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
auto p = stack.GetNextParticle();
// low energy process returns 1 kg/m²
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(1));
REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(3. / 4));
}
}
SECTION("high energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{4_TeV});
auto p = stack.GetNextParticle();
// high energy process returns 2 kg/m²
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2));
REQUIRE(seq.GetTotalInteractionLength(p) / kgMSq == Approx(6. / 5));
}
// high energy process creates 2 secondaries
SECTION("SelectInteraction") {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
InverseGrammageType invLambda = 0 / kgMSq;
switchProcess.SelectInteraction(p, projectile, 0.01 / kgMSq, invLambda);
REQUIRE(view.getSize() == 2);
}
}
}
TEST_CASE("SwitchProcess from ProcessSequence") {
DummyProcess<1> innerA;
DummyProcess<2> innerB;
DummyProcess<3> outer;
DummyProcess<4> additional;
auto seq = innerA << innerB;
switch_process::SwitchProcess switchProcess(seq, outer, 1_TeV);
auto completeSeq = switchProcess << additional;
SimpleStack stack;
SECTION("low energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{0.5_TeV});
auto p = stack.GetNextParticle();
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(2. / 3));
REQUIRE(completeSeq.GetTotalInteractionLength(p) / kgMSq == Approx(4. / 7));
}
SECTION("SelectInteraction") {
std::vector<int> numberOfSecondaries;
for (int i = 0; i < 1000; ++i) {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
double r = i / 1000.;
InverseGrammageType invLambda = r * 7. / 4 / kgMSq;
InverseGrammageType accumulator = 0 / kgMSq;
completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
numberOfSecondaries.push_back(view.getSize());
}
auto const mean =
std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
numberOfSecondaries.size();
REQUIRE(mean == Approx(12. / 7.).margin(0.01));
}
}
SECTION("high energy") {
stack.AddParticle(std::tuple<HEPEnergyType>{3.0_TeV});
auto p = stack.GetNextParticle();
SECTION("interaction length") {
REQUIRE(switchProcess.GetInteractionLength(p) / kgMSq == Approx(3));
REQUIRE(completeSeq.GetTotalInteractionLength(p) / kgMSq == Approx(12. / 7.));
}
SECTION("SelectInteraction") {
std::vector<int> numberOfSecondaries;
for (int i = 0; i < 1000; ++i) {
typename SimpleStack::ParticleType theParticle =
stack.GetNextParticle(); // as in corsika::Cascade
StackTestView view(theParticle);
auto projectile = view.GetProjectile();
double r = i / 1000.;
InverseGrammageType invLambda = r * 7. / 12. / kgMSq;
InverseGrammageType accumulator = 0 / kgMSq;
completeSeq.SelectInteraction(p, projectile, invLambda, accumulator);
numberOfSecondaries.push_back(view.getSize());
}
auto const mean =
std::accumulate(numberOfSecondaries.cbegin(), numberOfSecondaries.cend(), 0.) /
numberOfSecondaries.size();
REQUIRE(mean == Approx(24. / 7.).margin(0.01));
}
}
}
/*
* (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/utl/CorsikaData.h>
#include <cstdlib>
#include <stdexcept>
#include <string>
std::string corsika::utl::CorsikaData(std::string const& key) {
if (auto const* p = std::getenv("CORSIKA_DATA"); p != nullptr) {
auto const path = std::string(p) + "/" + key;
return path;
} else {
throw std::runtime_error("CORSIKA_DATA not set");
}
}
/*
* (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.
*/
#ifndef CORSIKA_CORSIKADATA_H
#define CORSIKA_CORSIKADATA_H
#include <string>
namespace corsika::utl {
/**
* returns the full path of the file \p filename within the CORSIKA_DATA directory
*/
std::string CorsikaData(std::string const& filename);
} // namespace corsika::utl
#endif // CORSIKA_CORSIKADATA_H
/*
* (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 <cmath>
#include "quartic.h"
//---------------------------------------------------------------------------
// solve cubic equation x^3 + a*x^2 + b*x + c
// x - array of size 3
// In case 3 real roots: => x[0], x[1], x[2], return 3
// 2 real roots: x[0], x[1], return 2
// 1 real root : x[0], x[1] ± i*x[2], return 1
unsigned int solveP3(double* x, double a, double b, double c) {
double a2 = a * a;
double q = (a2 - 3 * b) / 9;
double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54;
double r2 = r * r;
double q3 = q * q * q;
double A, B;
if (r2 < q3) {
double t = r / sqrt(q3);
if (t < -1) t = -1;
if (t > 1) t = 1;
t = acos(t);
a /= 3;
q = -2 * sqrt(q);
x[0] = q * cos(t / 3) - a;
x[1] = q * cos((t + M_2PI) / 3) - a;
x[2] = q * cos((t - M_2PI) / 3) - a;
return 3;
} else {
A = -pow(fabs(r) + sqrt(r2 - q3), 1. / 3);
if (r < 0) A = -A;
B = (0 == A ? 0 : q / A);
a /= 3;
x[0] = (A + B) - a;
x[1] = -0.5 * (A + B) - a;
x[2] = 0.5 * sqrt(3.) * (A - B);
if (fabs(x[2]) < eps) {
x[2] = x[1];
return 2;
}
return 1;
}
}
//---------------------------------------------------------------------------
// solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
// Attention - this function returns dynamically allocated array. It has to be released
// afterwards.
DComplex* solve_quartic(double a, double b, double c, double d) {
double a3 = -b;
double b3 = a * c - 4. * d;
double c3 = -a * a * d - c * c + 4. * b * d;
// cubic resolvent
// y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0
double x3[3];
unsigned int iZeroes = solveP3(x3, a3, b3, c3);
double q1, q2, p1, p2, D, sqD, y;
y = x3[0];
// The essence - choosing Y with maximal absolute value.
if (iZeroes != 1) {
if (fabs(x3[1]) > fabs(y)) y = x3[1];
if (fabs(x3[2]) > fabs(y)) y = x3[2];
}
// h1+h2 = y && h1*h2 = d <=> h^2 -y*h + d = 0 (h === q)
D = y * y - 4 * d;
if (fabs(D) < eps) // in other words - D==0
{
q1 = q2 = y * 0.5;
// g1+g2 = a && g1+g2 = b-y <=> g^2 - a*g + b-y = 0 (p === g)
D = a * a - 4 * (b - y);
if (fabs(D) < eps) // in other words - D==0
p1 = p2 = a * 0.5;
else {
sqD = sqrt(D);
p1 = (a + sqD) * 0.5;
p2 = (a - sqD) * 0.5;
}
} else {
sqD = sqrt(D);
q1 = (y + sqD) * 0.5;
q2 = (y - sqD) * 0.5;
// g1+g2 = a && g1*h2 + g2*h1 = c ( && g === p ) Krammer
p1 = (a * q1 - c) / (q1 - q2);
p2 = (c - a * q2) / (q1 - q2);
}
DComplex* retval = new DComplex[4];
// solving quadratic eq. - x^2 + p1*x + q1 = 0
D = p1 * p1 - 4 * q1;
if (D < 0.0) {
retval[0].real(-p1 * 0.5);
retval[0].imag(sqrt(-D) * 0.5);
retval[1] = std::conj(retval[0]);
} else {
sqD = sqrt(D);
retval[0].real((-p1 + sqD) * 0.5);
retval[1].real((-p1 - sqD) * 0.5);
}
// solving quadratic eq. - x^2 + p2*x + q2 = 0
D = p2 * p2 - 4 * q2;
if (D < 0.0) {
retval[2].real(-p2 * 0.5);
retval[2].imag(sqrt(-D) * 0.5);
retval[3] = std::conj(retval[2]);
} else {
sqD = sqrt(D);
retval[2].real((-p2 + sqD) * 0.5);
retval[3].real((-p2 - sqD) * 0.5);
}
return retval;
}
/***************************************************************************
* Copyright (C) 2016 by Саша Миленковић *
* sasa.milenkovic.xyz@gmail.com *
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* ( http://www.gnu.org/licenses/gpl-3.0.en.html ) *
* *
* You should have received a copy of the GNU General Public License *
* along with this program; if not, write to the *
* Free Software Foundation, Inc., *
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
***************************************************************************/
#ifndef QUARTIC_H_INCLUDED
#define QUARTIC_H_INCLUDED
#include <complex>
const double PI = 3.141592653589793238463L;
const double M_2PI = 2 * PI;
const double eps = 1e-12;
typedef std::complex<double> DComplex;
//---------------------------------------------------------------------------
// useful for testing
inline DComplex polinom_2(DComplex x, double a, double b) {
// Horner's scheme for x*x + a*x + b
return x * (x + a) + b;
}
//---------------------------------------------------------------------------
// useful for testing
inline DComplex polinom_3(DComplex x, double a, double b, double c) {
// Horner's scheme for x*x*x + a*x*x + b*x + c;
return x * (x * (x + a) + b) + c;
}
//---------------------------------------------------------------------------
// useful for testing
inline DComplex polinom_4(DComplex x, double a, double b, double c, double d) {
// Horner's scheme for x*x*x*x + a*x*x*x + b*x*x + c*x + d;
return x * (x * (x * (x + a) + b) + c) + d;
}
//---------------------------------------------------------------------------
// x - array of size 3
// In case 3 real roots: => x[0], x[1], x[2], return 3
// 2 real roots: x[0], x[1], return 2
// 1 real root : x[0], x[1] ± i*x[2], return 1
unsigned int solveP3(double* x, double a, double b, double c);
//---------------------------------------------------------------------------
// solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
// Attention - this function returns dynamically allocated array. It has to be released
// afterwards.
DComplex* solve_quartic(double a, double b, double c, double d);
#endif // QUARTIC_H_INCLUDED
This diff is collapsed.
Maintainers of the CORSIKA8 project are collaborators actively taking
care of code contributions, quality control, development, and related
discussions.
General and infrastructure:
- Ralf Ulrich <ralf.ulrich@kit.edu>, KIT
- Maximilian Reininghaus <maximilian.reininghaus@kit.edu>, KIT
- Hans Dembinski <hdembins@mpi-hd.mpg.de>, Dortmund
- Antonio Augusto Alves Junior <antonio.junior@kit.edu>, KIT
High performance, GPU:
- Dominik Baack <dominik.baack@tu-dortmund.de>, Dortmund
- Antonio Augusto Alves Junior <antonio.junior@kit.edu>, KIT
- Luisa Arrabito <arrabito@in2p3.fr>, Montpellier
Electromagnetic models:
- Jean-Marco Alameddine <jean-marco.alameddine@udo.edu>, Dortmund
- Jan Soedingrekso <jan.soedingrekso@tu-dortmund.de>, Dortmund
- Maximilian Sackel <maximilian.sackel@udo.edu>, Dortmund
Hadron models:
- Felix Riehn <friehn@lip.pt>, Santiago/Lisbon
- Anatoli Fedynitch <anatoli.fedynitch@icecube.wisc.edu> ICRR Tokyo
Radio:
- Remy Prechelt <prechelt@hawaii.edu>
- Tim Huege <tim.huege@kit.edu>, KIT
Testing, containers:
- Lukas Nellen <lukas@nucleares.unam.mx>
set (
MODEL_SOURCES
)
set (
MODEL_HEADERS
ExecTime.h
ExecTimeImpl.h
ImplBoundary.h
ImplContinuous.h
ImplDecay.h
ImplInteraction.h
ImplSecondaries.h
)
set (
MODEL_NAMESPACE
corsika/process/analytic_processors
)
add_library (AnalyticProcessors INTERFACE) #STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (AnalyticProcessors ${MODEL_NAMESPACE} ${MODEL_HEADERS})
target_link_libraries(AnalyticProcessors
INTERFACE
CORSIKAlogging)
target_include_directories (
AnalyticProcessors
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
TARGETS AnalyticProcessors
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
)
# --------------------
# code unit testing
CORSIKA_ADD_TEST (testExecTime
SOURCES
testExecTime.cc
${MODEL_HEADERS}
)
target_link_libraries (
testExecTime
AnalyticProcessors
ExampleProcessors
CORSIKAtesting
C8::ext::eigen3
CORSIKAthirdparty # for catch2
)
/*
* (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 <chrono>
#include <type_traits>
#include <corsika/logging/Logging.h>
#include <corsika/process/BoundaryCrossingProcess.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/process/DecayProcess.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/SecondariesProcess.h>
#include <corsika/process/StackProcess.h>
#include <corsika/process/analytic_processors/ImplBoundary.h>
#include <corsika/process/analytic_processors/ImplContinuous.h>
#include <corsika/process/analytic_processors/ImplDecay.h>
#include <corsika/process/analytic_processors/ImplInteraction.h>
#include <corsika/process/analytic_processors/ImplSecondaries.h>
namespace corsika::process {
namespace analytic_processors {
/// Time measurement of individual processes
/** This class allowes to log the runtime spend in all default calls to the process.
* No distinction is made between individual function calls of the process, the
* runtime is accumulated and or avaraged without differentiation.
*
* The class is currently only implemented for BoundaryProcess, ContinuousProcess,
* DecayProcess, InteractionProcess and SecondariesProcess and captures only the
* according functions of the base class given as template parameter. Trying to access
* BoundaryProcess functions with a DecayProcess as template parameter will currently
* give long errormessages.
*
* Inherits all functionality of the class that should be measured, this includes
* functions like getters and setters
*/
template <typename T>
class ExecTime
: public Boundary<T, std::is_base_of<corsika::process::BoundaryCrossingProcess<
typename T::TProcessType>,
T>::value>,
public Continuous<T, std::is_base_of<corsika::process::ContinuousProcess<
typename T::TProcessType>,
T>::value>,
public Decay<
T, std::is_base_of<corsika::process::DecayProcess<typename T::TProcessType>,
T>::value>,
public Interaction<T, std::is_base_of<corsika::process::InteractionProcess<
typename T::TProcessType>,
T>::value>,
public Secondaries<T, std::is_base_of<corsika::process::SecondariesProcess<
typename T::TProcessType>,
T>::value> {
static_assert(std::is_base_of<corsika::process::_BaseProcess, T>::value,
"error message");
public:
~ExecTime() {
C8LOG_INFO("Accumulated time spend in process {} is {} µs", typeid(T).name(),
this->sumTime());
}
};
} // namespace analytic_processors
} // namespace corsika::process
\ No newline at end of file
/*
* (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 <chrono>
#include <type_traits>
#include <corsika/logging/Logging.h>
#include <corsika/process/BoundaryCrossingProcess.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/process/DecayProcess.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/SecondariesProcess.h>
#include <corsika/process/StackProcess.h>
#include <corsika/process/analytic_processors/ImplBoundary.h>
#include <corsika/process/analytic_processors/ImplContinuous.h>
#include <corsika/process/analytic_processors/ImplDecay.h>
#include <corsika/process/analytic_processors/ImplInteraction.h>
#include <corsika/process/analytic_processors/ImplSecondaries.h>
namespace corsika::process {
namespace analytic_processors {
namespace detail {
/// Process type independent functionality of the Process runtime measurement class
/// ExecTime
/** Inherits all functionality of the class that should be sampled, this includes
* special functions for getters and setters
*
*
*
*/
template <typename T>
class ExecTimeImpl : public T {
private:
std::chrono::high_resolution_clock::time_point startTime_;
std::chrono::duration<double, std::micro> cumulatedTime_;
volatile double mean_;
volatile double mean2_;
volatile double min_;
volatile double max_;
volatile long long n_;
protected:
public:
using _T = T;
ExecTimeImpl() {
min_ = std::numeric_limits<double>::max();
cumulatedTime_ = std::chrono::duration<double, std::micro>(0);
max_ = 0;
mean_ = 0;
mean2_ = 0;
n_ = 0;
}
/// Starts the internal
inline void start() { startTime_ = std::chrono::high_resolution_clock::now(); }
/// Stops the internal timer and updates measurements
inline void stop() {
auto end = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::micro> timeDiv =
std::chrono::duration_cast<std::chrono::duration<double, std::micro> >(
end - startTime_);
this->update(timeDiv);
}
/// Updates the floating mean and variance as well as the global min and max of
/// the sampled runtimes
void update(std::chrono::duration<double, std::micro> timeDif) {
cumulatedTime_ += timeDif;
n_ = n_ + 1;
if (max_ < timeDif.count()) max_ = timeDif.count();
if (timeDif.count() < min_) min_ = timeDif.count();
double delta = timeDif.count() - mean_;
mean_ += delta / static_cast<double>(n_);
double delta2 = timeDif.count() - mean_;
mean2_ += delta * delta2;
}
inline double mean() const { return mean_; }
inline double min() const { return min_; }
inline double max() const { return max_; }
inline double var() const { return mean2_ / n_; }
inline double sumTime() const { return cumulatedTime_.count(); }
};
} // namespace detail
} // namespace analytic_processors
} // namespace corsika::process
\ No newline at end of file
/*
* (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/process/analytic_processors/ExecTimeImpl.h>
#include <corsika/analytics/ClassTimer.h>
namespace corsika::process {
namespace analytic_processors {
namespace detail {
template <typename T>
class ExecTimeImpl;
}
/// Base for Boundary Implementation
template <class T, bool TCheck>
class Boundary;
/// Specialisation if class is not BoundaryProcess
template <class T>
class Boundary<T, false> {};
/// Specialisation if class is a BoundaryProcess
template <class T>
class Boundary<T, true> : public detail::ExecTimeImpl<T> {
private:
public:
template <typename Particle, typename VTNType>
EProcessReturn DoBoundaryCrossing(Particle& p, VTNType const& from,
VTNType const& to) {
// Use of the ClassTimer function -> see ClassTimer for documentation
auto tc = corsika::analytics::ClassTimer<
EProcessReturn (detail::ExecTimeImpl<T>::_T::*)(Particle&, VTNType const&,
VTNType const&),
&detail::ExecTimeImpl<T>::_T::template DoBoundaryCrossing<Particle, VTNType>>(
*this);
EProcessReturn r = tc.call(p, from, to);
this->update(
std::chrono::duration_cast<std::chrono::duration<double, std::micro>>(
tc.getTime()));
return r;
}
};
} // namespace analytic_processors
} // namespace corsika::process
\ No newline at end of file