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 29 additions and 1809 deletions
# Guidelines for code development, structure, formating etc.
The CORSIKA Project very much welcomes contributions. Here we outlined
The CORSIKA Project very much welcomes contributions. Here we outline
how you can find the right place to contribute, and how to do that.
Connect to https://gitlab.ikp.kit.edu and corsika-devel@lists.kit.edu (self-register at https://www.lists.kit.edu/sympa/subscribe/corsika-devel) to get in touch with the project.
The CORSIKA Project decides on the [GUIDELINES](CONTRIBUTING.md) and can decide to
Connect to https://gitlab.iap.kit.edu and corsika-devel@lists.kit.edu (self-register at https://www.lists.kit.edu/sympa/subscribe/corsika-devel) to get in touch with the project.
The CORSIKA Project decides on the [contributing guidelines](CONTRIBUTING.md) and can decide to
change/improve them.
# How to contribute
- We organize all development via `Issues` that may be feature requests,
- We organize all development via gitlab `Issues` that may be feature requests,
ideas, discussions, or bugs fix requests.
- New issues can be created, or existing issues
picked up or contributed to.
......@@ -19,7 +19,8 @@ change/improve them.
created directly via the gitlab web interface.
- Proposed code to close one issue (located in a specific git
branch) is reviewed, discussed, and eventually merged
into the master branch to close the issue.
into the master branch via a merge-request (MR) to close the issue.
- all merge request will undergo a code review, and must be approved before merge, in order to ensure high code qualtiy: [Code Approval Procedure](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/-/wikis/Code-Approval-Procedure)
## Code formatting
......@@ -31,9 +32,9 @@ formatting. We provide a script `do-clang-format.sh`, which can be
very useful. But we urge everybody to integrate `clang-format` already
on the level of your source code editor. See [the official
page](https://clang.llvm.org/docs/ClangFormat.html) for information
about `clang-format` and its editor integration. At least: run
`do-clang-format.sh` prior to any `git add` command. Code with
improper formatting will not be accepted for merging.
about `clang-format` and its editor integrations. At least: run
`do-clang-format.sh` prior to any `git add/commit` command. Code with
improper formatting will not be accepted for merging. It will trigger automatic warnings by the continuous integration (CI) system.
The definition of source code format is written down in the file
[.clang-format](.clang-format) and can be changed, if the CORSIKA
......@@ -42,109 +43,35 @@ e.g. [link1](https://clangformat.com/) or
[link2](https://zed0.co.uk/clang-format-configurator/).
## Naming conventions
While `clang-format` does the structural formatting, we still need to agree on naming conventions:
- Classes and structs start with capital letters
- Class member variables start with "f"
- Any static variable has a "g" prefix. A static member variable starts with "fg"
- Class member functions start with capital letters
- Any class getter begins with "Get", and setter with "Set". Logical getters start with "Is" or "Has".
- enums should be "enum class" and start with a capital "E"
- Function parameter names start with "v"
- We use namespaces to avoid clashes and to structure code
- *Everything* is part of the corsika namespace
- All classes and objects are encapsulated into suited sub-namespaces,
thus corsika::geometry, corsika::processes, corsika::units, etc.
- Namespace names do not use capital letters.
- Every header file is copied during build and install into
"include/corsika/[namespace]" which also means, each header file
can only provide definitions for a _single_ namespace. It is one
main purpose of namespaces to structure the location of header
files.
- Each header file uses an include protection that includes at
least the namespace name, and header file name, thus, `#ifndef
__include_geometry_Point_h__` or `#ifndef __geometry_Point_h__`,
or similar are acceptable.
- Header files should always be included with `<..>`, thus,
`#include <corsika/geometry/Point.h>` since the build system
will always provide the correct include directives (and files
anyway cannot be found in file-system paths, which generally do
not follow the namespace naming conventions outlined
here).
- Header files are named after the main class (or object) they
define. This also means each header file name starts with a
capital letter.
## Coding rules
- Code may not introduce any compiler errors, or warnings (latter: on current CI runners)
- All unit tests must succeed at all times. If tests fail, code is not merged.
- We use C++17 concepts wherever useful and helpful
- On any major error or malfunction we throw an exception. This is needed and required for complex physics and shower debugging.
- We never catch exceptions for error handling, there might be very few special exceptions from this. We need to discuss such cases.
- Everything that should not change should be `const`
- Everything that does not need to be visible to the outside of a class/struct should be `private` or `protected`
- We prefer the use of references, wherever useful
- There cannot be any pointer in the interface of any class or object
exposed to outside users, there might be pointers for very special cases
inside of classes.
- When you contribute new code, or extend existing code, at the same time provide unit-tests for all functionality.
- When you contribute new physics algorithms, in addition you also need to provide a validation module
- Code must be documented with `doxygen` commands
## CMAKE formatting
- command are lower cases, e.g. ```set (...)```
- variables upper case, e.g. ```set (VAR1 Text)```
Since cmake itself lacks structure almost entirely:
- put a space between command and start of parenthesis, e.g. ```command (...)```
- add two spaces for logical indent
```
if (condition)
do something
endif (condition)
```
- break long lines to start with new keyword in new line (indented)
```
install (
FILES ${CORSIKAstackinterface_HEADERS}
DESTINATION include/${CORSIKAstackinterface_NAMESPACE}
)
```
- add plenty of comments to all cmake code
- use expressive variables and functions
## Coding rules, conventions and guidelines
Please read the [Coding wiki page](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/-/wikis/Coding-Conventions-and-Guidelines).
## Release versioning scheme
Releases of CORSIKA are thought to be the baseline for larger scale
testing, and full production. The releases are numbered as x.y.z,
starting with x=8 form the gitlab c++ version. X will only be
incremented for major design or physics changes. The y index is
updated for new official releases that normally contain improved or
enhanced physics performance, and may also contain normal interface
changes to accomodate improvements. The z index can be updated more
Releases of CORSIKA 8 are thought to be the baseline for larger scale
validation, and full production. The releases are numbered as CORSIKA 8 vx.y.z,
starting with x=1. The x index is updated for new releases that normally contain improved or
enhanced physics performance, and also interface
changes to accomodate improvements. The y and z indices can be updated more
frequently for bug fixes or new features. Changes in z will not
contain (major) interface changes, thus, production code will remain
fully compatible within changes of z. Special releases of CORSIKA can
also have a tag name from git, e.g. as in the "milestone1" release.
contain interface changes, thus, production code will remain
fully compatible within changes of z. Special releases of CORSIKA might
also have a release name.
# How to become scientific author of the CORSIKA Project
The CORSIKA Project decides on who becomes scientific author. The
following conditions are clearly sufficient, but not all of them are
following conditions are sufficient, but not all of them are
required all the time:
- responsibility for a particular functionality or software/management part
- have read and follow these [GUIDELINES](CONTRIBUTING.md)
- have read and follow our [contributing guidelines](CONTRIBUTING.md)
- active in the CORSIKA Project, that means responsive to
discussions and problems in corsika-devel@list.kit.edu or on https//gitlab.ikp.kit.edu,
of relevant *Issues*,
or in (phone) meetings
- agreement to the [COLLABORATION_AGREEMENT](COLLABORATION_AGREEMENT.md) is strictly required
- the members of the CORSIKA Project panel agree
discussions and problems in corsika-devel@list.kit.edu or on https//gitlab.iap.kit.edu,
of relevant *Issues*, in (phone) meetings or our Mattermost workspace
- agreement to the [using and collaborating agreement](USING_COLLABORATING.md)
- the members of the CORSIKA Project must agree
add_subdirectory (Doxygen)
add_subdirectory (Examples)
find_package (Doxygen OPTIONAL_COMPONENTS dot mscgen dia)
if (DOXYGEN_FOUND)
if (NOT DOXYGEN_DOT_EXECUTABLE)
message (FATAL_ERROR "Found doxygen but not 'dot' command, please install graphviz or set DOXYGEN_DOT_EXECUTABLE")
endif()
set (DOXYGEN_IN ${CMAKE_CURRENT_SOURCE_DIR}/Doxyfile.in)
set (DOXYGEN_OUT ${CMAKE_CURRENT_BINARY_DIR}/Doxyfile)
set (DOXYGEN_GENERATE_HTML YES)
set (DOXYGEN_GENERATE_MAN YES)
configure_file (${DOXYGEN_IN} ${DOXYGEN_OUT} @ONLY)
message ("Start doxygen with \"make doxygen\"")
# note the option ALL which allows to build the docs together with the application
add_custom_target (doxygen # ALL
COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_OUT}
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
COMMENT "Generating API documentation with Doxygen"
VERBATIM)
add_custom_command(TARGET doxygen POST_BUILD
COMMAND cd ${CMAKE_CURRENT_BINARY_DIR}/latex; pdflatex refman.tex
)
install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/html DESTINATION share/doc OPTIONAL)
install(FILES ${CMAKE_CURRENT_BINARY_DIR}/latex/refman.pdf DESTINATION share/doc OPTIONAL)
else (DOXYGEN_FOUND)
message ("Doxygen need to be installed to generate the doxygen documentation")
endif (DOXYGEN_FOUND)
PROJECT_NAME = CORSIKA 8
PROJECT_NUMBER = 0.0.0
GENERATE_HTML = YES
GENERATE_LATEX = YES
GENERATE_XML = YES
OUTPUT_DIRECTORY = @CMAKE_CURRENT_BINARY_DIR@/
INPUT = @PROJECT_SOURCE_DIR@ @PROJECT_BINARY_DIR@/Framework
EXCLUDE_PATTERNS = */ThirdParty/*/*
FULL_PATH_NAMES = YES
STRIP_FROM_PATH = @PROJECT_SOURCE_DIR@
FILE_PATTERNS = *.cc *.cpp *.cxx *.h *.dox *.inc *.md
EXTENSION_MAPPING = inc=C++
RECURSIVE = YES
SOURCE_BROWSER = YES
# INLINE_SOURCES
GENERATE_TREEVIEW = YES
CLASS_DIAGRAMS = NO
HAVE_DOT = YES
CLASS_GRAPH = YES
COLLABORATION_GRAPH = NO
UML_LOOK = NO
TEMPLATE_RELATIONS = YES
INCLUDE_GRAPH = YES
GRAPHICAL_HIERARCHY = YES
DIRECTORY_GRAPH = YES
GENERATE_LEGEND = YES
USE_MATHJAX = YES
SEARCHENGINE = YES
CORSIKA_ADD_TEST (helix_example)
target_link_libraries (helix_example CORSIKAgeometry CORSIKAunits)
install (TARGETS helix_example DESTINATION share/examples)
CORSIKA_ADD_TEST (geometry_example)
target_link_libraries (geometry_example CORSIKAgeometry CORSIKAunits)
install (TARGETS geometry_example DESTINATION share/examples)
CORSIKA_ADD_TEST (logger_example)
target_link_libraries (logger_example CORSIKAunits CORSIKAlogging)
install (TARGETS logger_example DESTINATION share/examples)
CORSIKA_ADD_TEST (stack_example)
target_link_libraries (stack_example SuperStupidStack CORSIKAunits
CORSIKAlogging)
# address sanitizer is making this example too slow, so we only do "undefined"
CORSIKA_ADD_TEST (cascade_example SANITIZE "undefined")
target_link_libraries (cascade_example
SuperStupidStack
CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
CORSIKAcascade
ProcessEnergyLoss
ProcessStackInspector
ProcessParticleCut
ProcessTrackWriter
ProcessTrackingLine
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
install (TARGETS cascade_example DESTINATION share/examples)
CORSIKA_ADD_TEST (boundary_example)
target_link_libraries (boundary_example
SuperStupidStack
CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
CORSIKAcascade
ProcessTrackWriter
ProcessParticleCut
ProcessTrackingLine
CORSIKAprocesses
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
install (TARGETS boundary_example DESTINATION share/examples)
if (Pythia8_FOUND)
CORSIKA_ADD_TEST (cascade_proton_example)
target_link_libraries (cascade_proton_example
SuperStupidStack
CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
ProcessPythia
CORSIKAcascade
ProcessEnergyLoss
ProcessTrackWriter
ProcessStackInspector
ProcessTrackingLine
ProcessParticleCut
ProcessHadronicElasticModel
ProcessStackInspector
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
install (TARGETS cascade_proton_example DESTINATION share/examples)
endif()
CORSIKA_ADD_TEST(vertical_EAS)
target_link_libraries (vertical_EAS
SuperStupidStack
CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
ProcessPythia
ProcessUrQMD
ProcessSwitch
CORSIKAcascade
ProcessEnergyLoss
ProcessObservationPlane
ProcessTrackWriter
ProcessTrackingLine
ProcessParticleCut
ProcessStackInspector
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
install (TARGETS vertical_EAS DESTINATION share/examples)
CORSIKA_ADD_TEST(stopping_power)
target_link_libraries (stopping_power
SuperStupidStack
CORSIKAunits
ProcessEnergyLoss
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
)
install (TARGETS vertical_EAS DESTINATION share/examples)
CORSIKA_ADD_TEST (staticsequence_example)
target_link_libraries (staticsequence_example
CORSIKAprocesssequence
CORSIKAunits
CORSIKAgeometry
CORSIKAlogging)
install (TARGETS staticsequence_example DESTINATION share/examples)
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
#include <corsika/utl/CorsikaFenv.h>
#include <iostream>
#include <limits>
#include <typeinfo>
using namespace corsika;
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;
using namespace std;
using namespace corsika::units::si;
template <bool deleteParticle>
struct MyBoundaryCrossingProcess
: public BoundaryCrossingProcess<MyBoundaryCrossingProcess<deleteParticle>> {
MyBoundaryCrossingProcess(std::string const& filename) { fFile.open(filename); }
template <typename Particle>
EProcessReturn DoBoundaryCrossing(Particle& p,
typename Particle::BaseNodeType const& from,
typename Particle::BaseNodeType const& to) {
std::cout << "boundary crossing! from: " << &from << "; to: " << &to << std::endl;
auto const& name = particles::GetName(p.GetPID());
auto const start = p.GetPosition().GetCoordinates();
fFile << name << " " << start[0] / 1_m << ' ' << start[1] / 1_m << ' '
<< start[2] / 1_m << '\n';
if constexpr (deleteParticle) { p.Delete(); }
return EProcessReturn::eOk;
}
void Init() {}
private:
std::ofstream fFile;
};
//
// The example main program for a particle cascade
//
int main() {
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
using EnvType = Environment<setup::IEnvironmentModel>;
EnvType env;
auto& universe = *(env.GetUniverse());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto outerMedium = EnvType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
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::Proton},
std::vector<float>{1.f}));
auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5_km);
innerMedium->SetModelProperties(props);
outerMedium->AddChild(std::move(innerMedium));
universe.AddChild(std::move(outerMedium));
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
process::sibyll::Interaction sibyll;
process::sibyll::Decay decay;
process::particle_cut::ParticleCut cut(20_GeV);
process::track_writer::TrackWriter trackWriter("tracks.dat");
MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
// assemble all processes into an ordered process list
auto sequence = sibyll << decay << cut << boundaryCrossing << trackWriter;
// setup particle stack, and add primary particles
setup::Stack stack;
stack.Clear();
const Code beamCode = Code::Proton;
const HEPMassType mass = particles::GetMass(Code::Proton);
const HEPEnergyType E0 = 50_TeV;
std::uniform_real_distribution distTheta(0., 180.);
std::uniform_real_distribution distPhi(0., 360.);
std::mt19937 rng;
for (int i = 0; i < 100; ++i) {
auto const theta = distTheta(rng);
auto const phi = distPhi(rng);
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
beamCode, E0, plab, pos, 0_ns});
}
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Init();
EAS.Run();
cout << "Result: E0=" << E0 / 1_GeV << endl;
cut.ShowResults();
const HEPEnergyType Efinal =
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
cout << "total energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
#include <corsika/utl/CorsikaFenv.h>
#include <iostream>
#include <limits>
using namespace corsika;
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;
using namespace std;
using namespace corsika::units::si;
//
// The example main program for a particle cascade
//
int main() {
const LengthType height_atmosphere = 112.8_km;
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
using EnvType = environment::Environment<setup::IEnvironmentModel>;
EnvType env;
auto& universe = *(env.GetUniverse());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto outerMedium = EnvType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
// fraction of oxygen
const float fox = 0.20946;
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 innerMedium = EnvType::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));
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
const Code beamCode = Code::Nucleus;
const int nuclA = 4;
const int nuclZ = int(nuclA / 2.15 + 0.7);
const HEPMassType mass = GetNucleusMass(nuclA, nuclZ);
const HEPEnergyType E0 = nuclA * 1_TeV;
double theta = 0.;
double phi = 0.;
Point const injectionPos(
rootCS, 0_m, 0_m,
height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
beamCode, E0, plab, injectionPos, 0_ns, nuclA, nuclZ});
}
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
process::sibyll::Decay decay;
// cascade with only HE model ==> HE cut
process::particle_cut::ParticleCut cut(80_GeV);
process::track_writer::TrackWriter trackWriter("tracks.dat");
process::energy_loss::EnergyLoss eLoss(
injectionPos, geometry::Vector<dimensionless_d>(rootCS, {0, 0, -1}));
// assemble all processes into an ordered process list
auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut
<< trackWriter;
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Init();
EAS.Run();
eLoss.PrintProfile(); // print longitudinal profile
cut.ShowResults();
const HEPEnergyType Efinal =
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
<< "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/pythia/Interaction.h>
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
#include <corsika/utl/CorsikaFenv.h>
#include <boost/type_index.hpp>
using boost::typeindex::type_id_with_cvr;
#include <iostream>
#include <limits>
#include <typeinfo>
using namespace corsika;
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;
using namespace std;
using namespace corsika::units::si;
//
// The example main program for a particle cascade
//
int main() {
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
using EnvType = Environment<setup::IEnvironmentModel>;
EnvType env;
auto& universe = *(env.GetUniverse());
auto theMedium =
EnvType::CreateNode<Sphere>(Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = HomogeneousMedium<IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
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();
const Code beamCode = Code::Proton;
const HEPMassType mass = particles::Proton::GetMass();
const HEPEnergyType E0 = 100_GeV;
double theta = 0.;
double phi = 0.;
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt(Elab * Elab - m * m);
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
beamCode, E0, plab, pos, 0_ns});
}
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
const std::vector<particles::Code> trackedHadrons = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
// process::sibyll::Interaction sibyll(env);
process::pythia::Interaction pythia;
// process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
// process::sibyll::Decay decay(trackedHadrons);
process::pythia::Decay decay(trackedHadrons);
process::particle_cut::ParticleCut cut(20_GeV);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
process::track_writer::TrackWriter trackWriter("tracks.dat");
// assemble all processes into an ordered process list
// auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;
auto sequence = pythia << decay << cut << trackWriter << stackInspect;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Init();
EAS.Run();
cout << "Result: E0=" << E0 / 1_GeV << endl;
cut.ShowResults();
const HEPEnergyType Efinal =
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
cout << "total energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <cstdlib>
#include <iostream>
#include <typeinfo>
using namespace corsika;
using namespace corsika::geometry;
using namespace corsika::units::si;
int main() {
// define the root coordinate system
geometry::CoordinateSystem& root =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// another CS defined by a translation relative to the root CS
CoordinateSystem cs2 = root.translate({0_m, 0_m, 1_m});
// rotations are possible, too; parameters are axis vector and angle
CoordinateSystem cs3 =
root.rotate(QuantityVector<length_d>{1_m, 0_m, 0_m}, 90 * degree_angle);
// now let's define some geometrical objects:
Point const p1(root, {0_m, 0_m, 0_m}); // the origin of the root CS
Point const p2(cs2, {0_m, 0_m, 0_m}); // the origin of cs2
Vector<length_d> const diff =
p2 -
p1; // the distance between the points, basically the translation vector given above
auto const norm = diff.squaredNorm(); // squared length with the right dimension
// print the components of the vector as given in the different CS
std::cout << "p2-p1 components in root: " << diff.GetComponents(root) << std::endl;
std::cout << "p2-p1 components in cs2: " << diff.GetComponents(cs2)
<< std::endl; // by definition invariant under translations
std::cout << "p2-p1 components in cs3: " << diff.GetComponents(cs3)
<< std::endl; // but not under rotations
std::cout << "p2-p1 norm^2: " << norm << std::endl;
assert(norm == 1 * meter * meter);
Sphere s(p1, 10_m); // define a sphere around a point with a radius
std::cout << "p1 inside s: " << s.Contains(p2) << std::endl;
assert(s.Contains(p2) == 1);
Sphere s2(p1, 3_um); // another sphere
std::cout << "p1 inside s2: " << s2.Contains(p2) << std::endl;
assert(s2.Contains(p2) == 0);
// let's try parallel projections:
auto const v1 = Vector<length_d>(root, {1_m, 1_m, 0_m});
auto const v2 = Vector<length_d>(root, {1_m, 0_m, 0_m});
auto const v3 = v1.parallelProjectionOnto(v2);
// cross product
auto const cross =
v1.cross(v2).normalized(); // normalized() returns dimensionless, normalized vectors
// if a CS is not given as parameter for getComponents(), the components
// in the "home" CS are returned
std::cout << "v1: " << v1.GetComponents() << std::endl;
std::cout << "v2: " << v2.GetComponents() << std::endl;
std::cout << "parallel projection of v1 onto v2: " << v3.GetComponents() << std::endl;
std::cout << "normalized cross product of v1 x v2" << cross.GetComponents()
<< std::endl;
return EXIT_SUCCESS;
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/geometry/Helix.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <array>
#include <cstdlib>
#include <iostream>
using namespace corsika;
using namespace corsika::geometry;
using namespace corsika::units::si;
int main() {
geometry::CoordinateSystem& root =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
Point const r0(root, {0_m, 0_m, 0_m});
auto const omegaC = 2 * M_PI * 1_Hz;
Vector<speed_d> vPar(root, {0_m / second, 0_m / second, 10_cm / second});
Vector<speed_d> vPerp(root, {1_m / second, 0_m / second, 0_m / second});
Helix h(r0, omegaC, vPar, vPerp);
auto constexpr t0 = 0_s;
auto constexpr t1 = 1_s;
auto constexpr dt = 1_ms;
auto constexpr n = long((t1 - t0) / dt) + 1;
auto arr = std::make_unique<std::array<std::array<double, 4>, n>>();
auto& positions = *arr;
for (auto [t, i] = std::tuple{t0, 0}; t < t1; t += dt, ++i) {
auto const r = h.GetPosition(t).GetCoordinates();
positions[i][0] = t / 1_s;
positions[i][1] = r[0] / 1_m;
positions[i][2] = r[1] / 1_m;
positions[i][3] = r[2] / 1_m;
}
std::cout << positions[n - 2][0] << " " << positions[n - 2][1] << " "
<< positions[n - 2][2] << " " << positions[n - 2][3] << std::endl;
return EXIT_SUCCESS;
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/logging/Logger.h>
#include <boost/format.hpp>
#include <fstream>
#include <iostream>
#include <string>
using namespace std;
using namespace corsika::logging;
int main() {
{
cout << "writing to \"another.log\"" << endl;
ofstream logfile("another.log");
sink::SinkStream unbuffered_sink(logfile);
sink::BufferedSinkStream sink(logfile, sink::StdBuffer(10000));
Logger<MessageOn, sink::BufferedSinkStream> info("\033[32m", "info", sink);
Logger<MessageOn, sink::BufferedSinkStream> err("\033[31m", "error", sink);
// logger<ostream,messageconst,StdBuffer> info(std::cout, StdBuffer(10000));
/*
Logging& logs = Logging::GetInstance();
logs.AddLogger<>("info", info);
auto& log_1 = logs.GetLogger("info"); // no so useful, since type of log_1 is
std::any
*/
for (int i = 0; i < 10000; ++i) {
LOG(info, "irgendwas", " ", string("and more"), " ",
boost::format("error: %i message: %s. done."), i, "stupido");
LOG(err, "Fehler");
}
}
{
sink::NoSink off;
Logger<MessageOff> info("", "", off);
for (int i = 0; i < 10000; ++i) {
LOG(info, "irgendwas", string("and more"),
boost::format("error: %i message: %s. done."), i, "stupido", "a-number:", 8.99,
"ENDE");
}
}
return 0;
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/particles/ParticleProperties.h>
#include <corsika/stack/super_stupid/SuperStupidStack.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <cassert>
#include <iomanip>
#include <iostream>
using namespace corsika;
using namespace corsika::units::si;
using namespace corsika::stack;
using namespace corsika::geometry;
using namespace std;
void fill(corsika::stack::super_stupid::SuperStupidStack& s) {
const geometry::CoordinateSystem& rootCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
for (int i = 0; i < 11; ++i) {
s.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Electron, 1.5_GeV * i,
corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 1_GeV}),
geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
}
}
void read(corsika::stack::super_stupid::SuperStupidStack& s) {
assert(s.GetSize() == 11); // stack has 11 particles
HEPEnergyType total_energy;
int i = 0;
for (auto& p : s) {
total_energy += p.GetEnergy();
// particles are electrons with 1.5 GeV energy times i
assert(p.GetPID() == particles::Code::Electron);
assert(p.GetEnergy() == 1.5_GeV * (i++));
}
}
int main() {
corsika::stack::super_stupid::SuperStupidStack s;
fill(s);
read(s);
return 0;
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <array>
#include <iomanip>
#include <iostream>
#include <corsika/process/ProcessSequence.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
using namespace corsika;
using namespace corsika::units::si;
using namespace corsika::process;
using namespace std;
const int nData = 10;
class Process1 : public BaseProcess<Process1> {
public:
Process1() {}
template <typename D, typename T, typename S>
EProcessReturn DoContinuous(D& d, T&, S&) const {
for (int i = 0; i < nData; ++i) d.p[i] += 1;
return EProcessReturn::eOk;
}
};
class Process2 : public BaseProcess<Process2> {
public:
Process2() {}
template <typename D, typename T, typename S>
inline EProcessReturn DoContinuous(D& d, T&, S&) const {
for (int i = 0; i < nData; ++i) d.p[i] -= 0.1 * i;
return EProcessReturn::eOk;
}
};
class Process3 : public BaseProcess<Process3> {
public:
Process3() {}
template <typename D, typename T, typename S>
inline EProcessReturn DoContinuous(D&, T&, S&) const {
return EProcessReturn::eOk;
}
};
class Process4 : public BaseProcess<Process4> {
public:
Process4(const double v)
: fV(v) {}
template <typename D, typename T, typename S>
inline EProcessReturn DoContinuous(D& d, T&, S&) const {
for (int i = 0; i < nData; ++i) d.p[i] *= fV;
return EProcessReturn::eOk;
}
private:
double fV;
};
struct DummyData {
double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
};
struct DummyStack {
void clear() {}
};
struct DummyTrajectory {};
void modular() {
Process1 m1;
Process2 m2;
Process3 m3;
Process4 m4(0.9);
auto sequence = m1 << m2 << m3 << m4;
DummyData particle;
DummyTrajectory track;
const int n = 1000;
for (int i = 0; i < n; ++i) { sequence.DoContinuous(particle, track); }
for (int i = 0; i < nData; ++i) {
// cout << p.p[i] << endl;
// assert(p.p[i] == n-i*100);
}
cout << " done (nothing...) " << endl;
}
int main() {
modular();
return 0;
}
/*
* (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/environment/Environment.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/CorsikaFenv.h>
#include <fstream>
#include <iostream>
#include <limits>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::particles;
using namespace corsika::geometry;
using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
//
// This example demonstrates the energy loss of muons as function of beta*gamma (=p/m)
//
int main() {
feenableexcept(FE_INVALID);
// setup environment, geometry
using EnvType = Environment<setup::IEnvironmentModel>;
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const injectionPos(
rootCS, 0_m, 0_m,
112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
Vector<dimensionless_d> showerAxis(rootCS, {0, 0, -1});
process::energy_loss::EnergyLoss eLoss(injectionPos, showerAxis);
setup::Stack stack;
std::ofstream file("dEdX.dat");
file << "# beta*gamma, dE/dX / eV/(g/cm²)" << std::endl;
for (HEPEnergyType E0 = 300_MeV; E0 < 1_PeV; E0 *= 1.05) {
stack.Clear();
const Code beamCode = Code::MuPlus;
const HEPMassType mass = GetMass(beamCode);
double theta = 0.;
double phi = 0.;
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
beamCode, E0, plab, injectionPos, 0_ns});
auto const p = stack.GetNextParticle();
HEPEnergyType dE = eLoss.TotalEnergyLoss(p, 1_g / square(1_cm));
file << P0 / mass << "\t" << -dE / 1_eV << std::endl;
}
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/StackProcess.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/observation_plane/ObservationPlane.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/switch_process/SwitchProcess.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/FlatExponential.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Plane.h>
#include <corsika/geometry/Sphere.h>
//~ #include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/urqmd/UrQMD.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
#include <corsika/utl/CorsikaFenv.h>
#include <iomanip>
#include <iostream>
#include <limits>
#include <typeinfo>
using namespace corsika;
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;
using namespace std;
using namespace corsika::units::si;
void registerRandomStreams() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
random::RNGManager::GetInstance().RegisterRandomStream("UrQMD");
random::RNGManager::GetInstance().SeedAll();
}
int main() {
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
registerRandomStreams();
// setup environment, geometry
using EnvType = Environment<setup::IEnvironmentModel>;
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto& universe = *(env.GetUniverse());
auto theMedium = EnvType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
auto constexpr temperature = 295_K; // AIRES default temperature for isothermal model
auto constexpr lambda =
-constants::R * temperature / (constants::g_sub_n * 28.966_g / mole);
// 1036 g/cm² taken from AIRES code
auto constexpr rho0 = -1036_g / square(1_cm) / lambda;
using FlatExp = environment::FlatExponential<environment::IMediumModel>;
theMedium->SetModelProperties<FlatExp>(
Point{rootCS, 0_m, 0_m, 0_m}, Vector<dimensionless_d>{rootCS, {0., 0., 1.}}, rho0,
lambda,
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Nitrogen,
particles::Code::Oxygen},
std::vector<float>{
0.7847f,
1.f - 0.7847f})); // values taken from AIRES manual, Ar removed for now
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
const Code beamCode = Code::Proton;
auto const mass = particles::GetMass(beamCode);
const HEPEnergyType E0 = 0.1_PeV;
double theta = 0.;
double phi = 0.;
Point const injectionPos(
rootCS, 0_m, 0_m, 112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
// {
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
beamCode, E0, plab, injectionPos, 0_ns});
// }
Line const line(injectionPos, plab.normalized() * 1_m * 1_Hz);
auto const velocity = line.GetV0().norm();
auto const observationHeight = 1.425_km;
setup::Trajectory const showerAxis(line, (112.8_km - observationHeight) / velocity);
auto const grammage = theMedium->GetModelProperties().IntegratedGrammage(
showerAxis, (112.8_km - observationHeight));
std::cout << "Grammage to ground: " << grammage / (1_g / square(1_cm)) << " g/cm²"
<< std::endl;
universe.AddChild(std::move(theMedium));
// setup processes, decays and interactions
const std::vector<particles::Code> trackedHadrons = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
//~ process::sibyll::Decay decay(trackedHadrons);
process::pythia::Decay decay(trackedHadrons);
process::particle_cut::ParticleCut cut(5_GeV);
process::track_writer::TrackWriter trackWriter("tracks.dat");
process::energy_loss::EnergyLoss eLoss(showerAxis);
Plane const obsPlane(Point(rootCS, 0_m, 0_m, observationHeight),
Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
process::observation_plane::ObservationPlane observationLevel(obsPlane,
"particles.dat");
// assemble all processes into an ordered process list
process::UrQMD::UrQMD urqmd;
auto sibyllSequence = sibyll << sibyllNuc;
process::switch_process::SwitchProcess switchProcess(urqmd, sibyllSequence, 55_GeV);
auto sequence = switchProcess << decay << eLoss << cut << observationLevel
<< trackWriter;
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Init();
EAS.Run();
eLoss.PrintProfile(); // print longitudinal profile
cut.ShowResults();
const HEPEnergyType Efinal =
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
<< "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
std::ofstream finish("finished");
finish << "run completed without error" << std::endl;
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_Environment_BaseExponential_h_
#define _include_Environment_BaseExponential_h_
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
#include <limits>
namespace corsika::environment {
/**
* This class provides the grammage/length conversion functionality for
* (locally) flat exponential atmospheres.
*/
template <class TDerived>
class BaseExponential {
protected:
units::si::MassDensityType const fRho0;
units::si::LengthType const fLambda;
units::si::InverseLengthType const fInvLambda;
geometry::Point const fP0;
auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
// clang-format off
/**
* For a (normalized) axis \f$ \vec{a} \f$, the grammage along a non-orthogonal line with (normalized)
* direction \f$ \vec{u} \f$ is given by
* \f[
* X = \frac{\varrho_0 \lambda}{\vec{u} \cdot \vec{a}} \left( \exp\left( \vec{u} \cdot \vec{a} \frac{l}{\lambda} \right) - 1 \right)
* \f], where \f$ \varrho_0 \f$ is the density at the starting point.
*
* If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
* \f[
* X = \varrho_0 l;
* \f]
*/
// clang-format on
units::si::GrammageType IntegratedGrammage(
geometry::Trajectory<geometry::Line> const& vLine, units::si::LengthType vL,
geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
if (vL == units::si::LengthType::zero()) { return units::si::GrammageType::zero(); }
auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
if (uDotA == 0) {
return vL * rhoStart;
} else {
return rhoStart * (fLambda / uDotA) * (exp(uDotA * vL * fInvLambda) - 1);
}
}
// clang-format off
/**
* For a (normalized) axis \f$ \vec{a} \f$, the length of a non-orthogonal line with (normalized)
* direction \f$ \vec{u} \f$ corresponding to grammage \f$ X \f$ is given by
* \f[
* l = \begin{cases}
* \frac{\lambda}{\vec{u} \cdot \vec{a}} \log\left(Y \right), & \text{if} Y := 0 > 1 +
* \vec{u} \cdot \vec{a} \frac{X}{\rho_0 \lambda}
* \infty & \text{else,}
* \end{cases}
* \f] where \f$ \varrho_0 \f$ is the density at the starting point.
*
* If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
* \f[
* l = \frac{X}{\varrho_0}
* \f]
*/
// clang-format on
units::si::LengthType ArclengthFromGrammage(
geometry::Trajectory<geometry::Line> const& vLine,
units::si::GrammageType vGrammage,
geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
if (uDotA == 0) {
return vGrammage / rhoStart;
} else {
auto const logArg = vGrammage * fInvLambda * uDotA / rhoStart + 1;
if (logArg > 0) {
return fLambda / uDotA * log(logArg);
} else {
return std::numeric_limits<typename decltype(
vGrammage)::value_type>::infinity() *
units::si::meter;
}
}
}
public:
BaseExponential(geometry::Point const& vP0, units::si::MassDensityType vRho,
units::si::LengthType vLambda)
: fRho0(vRho)
, fLambda(vLambda)
, fInvLambda(1 / vLambda)
, fP0(vP0) {}
};
} // namespace corsika::environment
#endif
set (
ENVIRONMENT_HEADERS
VolumeTreeNode.h
IMediumModel.h
NuclearComposition.h
HomogeneousMedium.h
InhomogeneousMedium.h
HomogeneousMedium.h
LinearApproximationIntegrator.h
DensityFunction.h
Environment.h
NameModel.h
BaseExponential.h
FlatExponential.h
SlidingPlanarExponential.h
)
set (
ENVIRONMENT_NAMESPACE
corsika/environment
)
add_library (CORSIKAenvironment INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAenvironment ${ENVIRONMENT_NAMESPACE} ${ENVIRONMENT_HEADERS})
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
CORSIKAenvironment
INTERFACE
CORSIKAgeometry
CORSIKAparticles
CORSIKAunits
CORSIKArandom
)
target_include_directories (
CORSIKAenvironment
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include>
)
install (
TARGETS CORSIKAenvironment
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
PUBLIC_HEADER DESTINATION include/${ENVIRONMENT_NAMESPACE}
)
# --------------------
# code unit testing
CORSIKA_ADD_TEST(testEnvironment)
target_link_libraries (
testEnvironment
CORSIKAsetup
CORSIKAenvironment
CORSIKAtesting
)
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_environment_DensityFunction_h_
#define _include_environment_DensityFunction_h_
#include <corsika/environment/LinearApproximationIntegrator.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
namespace corsika::environment {
template <class TDerivableRho,
template <typename> class TIntegrator = LinearApproximationIntegrator>
class DensityFunction
: public TIntegrator<DensityFunction<TDerivableRho, TIntegrator>> {
friend class TIntegrator<DensityFunction<TDerivableRho, TIntegrator>>;
TDerivableRho fRho; //!< functor for density
public:
DensityFunction(TDerivableRho rho)
: fRho(rho) {}
corsika::units::si::MassDensityType EvaluateAt(
corsika::geometry::Point const& p) const {
return fRho(p);
}
};
} // namespace corsika::environment
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_environment_Environment_h
#define _include_environment_Environment_h
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Sphere.h>
#include <limits>
namespace corsika::environment {
struct Universe : public corsika::geometry::Sphere {
Universe(corsika::geometry::CoordinateSystem const& pCS)
: corsika::geometry::Sphere(
corsika::geometry::Point{pCS, 0 * corsika::units::si::meter,
0 * corsika::units::si::meter,
0 * corsika::units::si::meter},
corsika::units::si::meter * std::numeric_limits<double>::infinity()) {}
bool Contains(corsika::geometry::Point const&) const override { return true; }
};
template <typename IEnvironmentModel>
class Environment {
public:
using BaseNodeType = VolumeTreeNode<IEnvironmentModel>;
Environment()
: fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem()}
, fUniverse(std::make_unique<BaseNodeType>(
std::make_unique<Universe>(fCoordinateSystem))) {}
// using IEnvironmentModel = corsika::setup::IEnvironmentModel;
auto& GetUniverse() { return fUniverse; }
auto const& GetUniverse() const { return fUniverse; }
auto const& GetCoordinateSystem() const { return fCoordinateSystem; }
// factory method for creation of VolumeTreeNodes
template <class TVolumeType, typename... TVolumeArgs>
static auto CreateNode(TVolumeArgs&&... args) {
static_assert(std::is_base_of_v<corsika::geometry::Volume, TVolumeType>,
"unusable type provided, needs to be derived from "
"\"corsika::geometry::Volume\"");
return std::make_unique<BaseNodeType>(
std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
}
private:
corsika::geometry::CoordinateSystem const& fCoordinateSystem;
typename BaseNodeType::VTNUPtr fUniverse;
};
// using SetupBaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
// using SetupEnvironment = Environment<corsika::setup::IEnvironmentModel>;
} // namespace corsika::environment
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_Environment_FlatExponential_h_
#define _include_Environment_FlatExponential_h_
#include <corsika/environment/BaseExponential.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::environment {
//clang-format off
/**
* flat exponential density distribution with
* \f[
* \varrho(r) = \varrho_0 \exp\left( \frac{1}{\lambda} (r - p) \cdot
* \vec{a} \right).
* \f]
* \f$ \vec{a} \f$ denotes the axis and should be normalized to avoid degeneracy
* with the scale parameter \f$ \lambda \f$.
*/
//clang-format on
template <class T>
class FlatExponential : public BaseExponential<FlatExponential<T>>, public T {
geometry::Vector<units::si::dimensionless_d> const fAxis;
NuclearComposition const fNuclComp;
using Base = BaseExponential<FlatExponential<T>>;
public:
FlatExponential(geometry::Point const& vP0,
geometry::Vector<units::si::dimensionless_d> const& vAxis,
units::si::MassDensityType vRho, units::si::LengthType vLambda,
NuclearComposition vNuclComp)
: Base(vP0, vRho, vLambda)
, fAxis(vAxis)
, fNuclComp(vNuclComp) {}
units::si::MassDensityType GetMassDensity(geometry::Point const& vP) const override {
return Base::fRho0 * exp(Base::fInvLambda * (vP - Base::fP0).dot(fAxis));
}
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
units::si::GrammageType IntegratedGrammage(
geometry::Trajectory<geometry::Line> const& vLine,
units::si::LengthType vTo) const override {
return Base::IntegratedGrammage(vLine, vTo, fAxis);
}
units::si::LengthType ArclengthFromGrammage(
geometry::Trajectory<geometry::Line> const& vLine,
units::si::GrammageType vGrammage) const override {
return Base::ArclengthFromGrammage(vLine, vGrammage, fAxis);
}
};
} // namespace corsika::environment
#endif