IAP GITLAB

Skip to content
Snippets Groups Projects
Commit cd66e584 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

Merge branch '41-implement-simple-homogenous-atmosphere-model' of...

Merge branch '41-implement-simple-homogenous-atmosphere-model' of gitlab.ikp.kit.edu:AirShowerPhysics/corsika into 41-implement-simple-homogenous-atmosphere-model
parents 03d6cedf efaf7afc
No related branches found
No related tags found
No related merge requests found
Showing
with 816 additions and 63 deletions
**/*~
**/*.bak
**/*log
build/
Framework/Particles/GeneratedParticleProperties.inc
flymd.html
flymd.md
Maximilian Reininghaus, maximilian.reininghaus@kit.edu
Ralf Ulrich, ralf.ulrich@kit.edu
Everything
\ No newline at end of file
* milestone1 release: So 7. Okt 15:51:07 CEST 2018
......@@ -7,6 +7,8 @@ project (
LANGUAGES CXX
)
enable_language (Fortran)
# ignore many irrelevant Up-to-date messages during install
set (CMAKE_INSTALL_MESSAGE LAZY)
......@@ -15,13 +17,28 @@ set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/CMakeModules)
include (CorsikaUtilities) # a few cmake function
# enable warnings and disallow non-standard language
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -pedantic -Wextra")
add_definitions(-Wall -pedantic -Wextra -Wno-ignored-qualifiers)
# --std=c++17
set (CMAKE_CXX_STANDARD 17)
enable_testing ()
set (CTEST_OUTPUT_ON_FAILURE 1)
# Set a default build type if none was specified
set(default_build_type "Release")
if(EXISTS "${CMAKE_SOURCE_DIR}/.git")
set(default_build_type "Debug")
endif()
if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
message(STATUS "Setting build type to '${default_bluild_type}' as no other was specified.")
set(CMAKE_BUILD_TYPE "${default_build_type}" CACHE
STRING "Choose the type of build." FORCE)
# Set the possible values of build type for cmake-gui
set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS
"Debug" "Release" "MinSizeRel" "RelWithDebInfo")
endif()
# unit testing coverage, does not work yet
#include (CodeCoverage)
##set(COVERAGE_LCOV_EXCLUDES 'Documentation/*')
......@@ -32,21 +49,22 @@ set (CTEST_OUTPUT_ON_FAILURE 1)
# #-j ${PROCESSOR_COUNT}
# # DEPENDENCIES corsika
# )
#add_custom_target (corsika_pre_build)
#add_custom_command (TARGET corsika_pre_build PRE_BUILD COMMAND "${PROJECT_SOURCE_DIR}/pre_compile.py")
# dependencies
find_package (Boost 1.40 COMPONENTS program_options REQUIRED)
find_package (Boost 1.60 REQUIRED)
find_package (Eigen3 REQUIRED)
#find_package (HDF5) # not yet needed
# order of subdirectories
# order of subdirectories
add_subdirectory (ThirdParty)
#add_subdirectory (Utilities)
add_subdirectory (Framework)
add_subdirectory (Environment)
add_subdirectory (Stack)
add_subdirectory (Setup)
add_subdirectory (Processes)
add_subdirectory (Documentation)
add_subdirectory (Main)
......@@ -55,3 +55,26 @@ function (CORSIKA_COPY_HEADERS_TO_NAMESPACE for_library in_namespace)
endfunction (CORSIKA_COPY_HEADERS_TO_NAMESPACE)
#
# use: CORSIKA_ADD_FILES_ABSOLUTE varname
#
# add list of filenames with absolute paths (pointing to CMAKE_SOURCE_DIR) to ${varname} in PARAENT_SCOPE
#
macro (CORSIKA_ADD_FILES_ABSOLUTE varname)
file (RELATIVE_PATH _relPath "${PROJECT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}")
foreach (_src ${ARGN})
if (_relPath)
list (APPEND "${varname}" "${CMAKE_SOURCE_DIR}/${_relPath}/${_src}")
else()
list (APPEND "${varname}" "${CMAKE_SOURCE_DIR}/${_src}")
endif()
endforeach()
if (_relPath)
# propagate SRCS to parent directory
set ("${varname}" "${${varname}}" PARENT_SCOPE)
endif()
endmacro(CORSIKA_ADD_FILES_ABSOLUTE)
# Collaboration agreement
The CORSIKA project very much welcomes all collaboration and
contributions. The aim of the CORSIKA project is to create a
scientific software framework as a fundamental tool for research. The
collaboration agreement and the licensing model are based on the
guidelines layed out by HSF
[[1]](https://hepsoftwarefoundation.org/activities/licensing.html) or
CERN
[[2]](https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&ved=2ahUKEwiLqKG00dXdAhUOZFAKHdIwAh4QFjAAegQIARAC&url=https%3A%2F%2Findico.cern.ch%2Fcategory%2F4251%2Fattachments%2F101%2F505%2FOSL-2012-01-Open_Source_Licences_at_CERN-Short_version.pdf&usg=AOvVaw1n4S0PQCSeE6wbdfdhKDqF),
[[3]](http://legal.web.cern.ch/licensing/software), and follow the
examples of other big scientific software projects.
The CORSIKA project consists of the contributions from the scientific
community and individuals in a best effort to deliver the best
possible computing performance and physics output.
The MCnet guidelines developed by [www.montecarlonet.org](www.montecarlonet.org)
are copied in [MCNET_GUIDELINES](MCNET_GUIDELINES) -- they provide a very good
additional scope that contributors should read and follow.
All possible
liability and licensing question are only handled by the adopted
software license.
## The software license of the CORSIKA project
The license adopted for the CORSIKA project is the explicit copyleft
license GPLv3, as copied in full in the file
[LICENSE](LICENSE). Each source file of the CORSIKA project contains a
short statement of the copyright and this license. Each binary or
source code release of CORSIKA must contain the file LICENSE. The
code, documentation and content in the folder [ThirdParty](ThirdParty)
is not integral part of the CORSIKA project and can be based on or
include other licenses, which must be compatible with GPLv3. Check the
content of this folder for details. It depends on the configuration of
the build system to what extend this code is used to build CORSIKA.
## Who is the "copyright holder"
For legal reasons and the ability to maintain the CORSIKA project
effectively over a very long lifespan of several decades, all
contributors are required to transfer their copyright to the CORSIKA
Project. Of course you will be duly credited and your name will appear
on the contributors page and in the [AUTHORS](AUTHORS) file shipped
with every binary and source distribution. The copyright transfer is
necessary to be able to effectively defend the project in case of
litigation. The copyright holder may change, if decided by the CORSIKA
Project. The current copyright holder is the CORSIKA Project
corsika-project@lists.kit.edu, with the current chair person Ralf Ulrich (KIT) ralf.ulrich@kit.edu.
## Definition of a "contributor"
Contributor is a person of whom at least one merge request was
accepted for the master branch of the CORSIKA project at
[https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika).
All contributors will be co-listed and credited as (software) authors
of the CORSIKA project, as well as listed indefinitely in the
[AUTHORS](AUTHORS) file. Contributors should add their name and
contact data to the [AUTHORS](AUTHORS) file, as part of one of their
merge requests. This file is always distributed together with all
source and binary releases. If you contribute to any non-master
branch, you can add your name to the [AUTHORS](AUTHORS) file of this
particular branch, but all official releases are normally performed
via the master branch.
If you want to contribute, you need to read
[GUIDELINES](GUIDELINES.md) and comply with these rules, or help to
improve them.
## Definition and working mode of the CORSIKA Project panel
The CORSIKA Project panel makes all decisions for the CORSIKA
Project. It can also change the
[COLLABORATION\_AGREEMENT](COLLABORATION\_AGREEMENT.md), the
[GUIDELINES](GUIDELINES.md) or any other structure or document relevant for the CORSIKA Project.
The CORSIKA Project *panel* consists (October 2018) of
* Ralph Engel (KIT)
* Dieter Heck (KIT)
* Tanguy Pierog (KIT)
* Maximilian Reininghaus (KIT)
* Markus Roth (KIT)
* Ralf Ulrich (KIT)
* Michael Unger (KIT)
* Darko Veberic (KIT)
and can be contacted via corsika-project@lists.kit.edu. The chair
person of the CORSIKA Project is Ralf Ulrich (KIT). Members of the
CORSIKA Project *panel* are *Maintainer* of the CORSIKA Project in
gitlab at
[https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika),
and have special responsibilities for this reason.
## Changing to a different license, for parts, or the complete project
The CORSIKA Project panel can change the license for parts or the entire project.
## Planning and performing releases
The CORSIKA Project panel decides on releases of the software, and about the content of it.
## Changes to the Collaboration Agreement
The CORSIKA Project panel decides on changes to the Collaboration
agreement.
# Contributing
The CORSIKA project is committed to fostering and preserving a
diverse, welcoming community; all participants are expected to
respect that.
- [Getting Started](#getting-started)
- [Very first steps](#very-first-steps)
- [Opening issues](#opening-issues)
- [Participating in discussions](#participating-in-discussions)
- [Improving and reviewing docs](#improving-and-reviewing-docs)
- [Reviewing and testing changes](#reviewing-and-testing-changes)
- [Proposing and making changes](#proposing-and-making-changes)
- [Finding something to work on](#finding-something-to-work-on)
- [Before you start](#before-you-start-work)
- [Before you open your PR](#before-you-open-your-pr)
- [Review process](#review-process)
- [Getting more involved](#getting-more-involved)
## Getting started
Not sure where to start? If you haven't already, take a look at the
[docs](http://xi-editor.github.io/xi-editor/docs.html) to get a better
sense of the project. Read through some issues and some open PRs, to
get a sense for the habits of existing contributors. Drop by the #xi
channel on [irc.mozilla.org](https://mozilla.logbot.info/xi) to follow
ongoing discussions or ask questions. Clone the repos you're
interested in, and make sure you can build and run the tests. If you
can't, open an issue, and someone will try to help. Once you're up and
running, there are a a number of ways to participate:
### Opening issues
If you have a question or a feature request or think you've found a bug,
please open an issue. When opening an issue, include any details that
might be relevant: for a bug this might be the steps required to
reproduce; for a feature request it might be a detailed explanation of
the behaviour you are imagining, an outline of how it would be used,
and/or examples of how this feature is used in other editors.
#### Before you open an issue
Before opening an issue, **try to identify where the issue belongs**.
Is it a problem with the frontend or with core? The frontend is
responsible for drawing windows and UI, and handling events; the core
is responsible for most everything else. Issues with the frontend
should be opened in that frontend's repository, and issues with
core should be opened in the
[xi-editor](https://github.com/xi-editor/xi-editor/issues) repo.
Finally, before opening an issue, **use github's search bar** to make
sure there isn't an existing (open or closed) issue for your particular
problem.
### Participating in discussions
An _explicit_ goal of xi-editor is to be an educational resource.
Everyone is encouraged to participate in discussion issues (issues with
the 'discussion' or 'planning' labels), and we expect people
participating in discussions to be respectful of the fact that we all
have different backgrounds and levels of experience. Similarly, if
something is confusing, feel free to ask for clarification! If you're
confused, other people probably are as well.
### Improving and reviewing docs
If the docs are unclear or incomplete, please open an issue or a PR to
improve them. This is an especially valuable form of feedback from new
contributors, who are seeing things for the first time, and will be best
positioned to identify areas that need improvement.
### Reviewing and testing changes
One of the best ways to get more familiar with the project is by reading
other people's pull requests. If there's something in a commit that you
don't understand, this is a great time to ask for clarification. Testing
changes is also very helpful, especially for bug fixes or feature
additions. Check out a change and try it out; does it work? Can you find
edge cases? Manual testing is very valuable. For more information on
reviews, see [code review process](#review-process).
## Proposing and making changes
### Finding something to work on
If you're looking for something to work on, a good first step is to browse
the [issues](https://github.com/xi-editor/xi-editor/issues). Specifically,
issues that are labeled
[help wanted](https://github.com/xi-editor/xi-editor/issues?q=is%3Aissue+is%3Aopen+label%3A%22help+wanted%22) and/or
[easy](https://github.com/xi-editor/xi-editor/issues?q=is%3Aissue+is%3Aopen+label%3Aeasy)
are good places to start. If you can't find anything there, feel free to ask
on IRC, or play around with the editor and try to identify something that
_you_ think is missing.
### Before you start work
Before starting to work on an issue, consider the following:
- _Is it a bugfix or small change?_ If you notice a small bug somewhere,
and you believe you have a fix, feel free to open a pull request directly.
- _Is it a feature?_ If you have an idea for a new editor feature that is
along the lines of something that already exists (for instance, adding a
new command to reverse the letters in a selected region) _consider_
opening a short issue beforehand, describing the feature you have in mind.
Other contributors might be able to identify possible issues or
refinements. This isn't _necessary_, but it might end up saving you work,
and it means you will get to close an issue when your PR gets merged,
which feels good.
- _Is it a major feature, affecting for instance the behaviour or appearance
of a frontend, or the API or architecture of core?_ Before working on a
large change, please open a discussion/proposal issue. This should describe
the problem you're trying to solve, and the approach you're considering;
think of this as a 'lite' version of Rust's
[RFC](https://github.com/rust-lang/rfcs) process.
### Before you open your PR
Before pressing the 'Create pull request' button,
- _Run the tests_. It's easy to accidentally break something with even a small
change, so always run the tests locally before submitting (or updating) a PR.
You can run all checks locally with the `xi-editor/rust/run_all_checks`. script.
- _Add a message for your reviewers_. When submitting a PR, take advantage
of the opportunity to include a message. Your goal here should be to help
your reviewers. Are there any parts of your change that you're uncertain
about? Are there any non-obvious explanations for some of your decisions?
If your change changes some behaviour, how might it be tested?
- ***Be your own first reviewer***. On the page where you enter your message,
you have a final opportunity to see your PR _as it will be seen by your
reviewers_. This is a great opportunity to give it one last review, yourself.
Imagine that it is someone else's work, that you're reviewing: what comments
would you have? If you spot a typo or a problem, you can push an update in
place, without losing your PR message or other state.
- _Add yourself to the AUTHORS file_. If this is your first substantive pull
request in this repo, feel free to add yourself to the AUTHORS file.
### Review process
Every non-trivial pull request goes through review. Everyone is welcome to
participate in review; review is an excellent time to ask questions about
code or design decisions that you don't understand.
All pull requests must be approved by an appropriate reviewer before they
are merged. For bug fixes and smaller changes, this can be anyone who has
commit rights to the repo in question; larger changes (changes which add a
feature, or significantly change behaviour or API) should also be approved by
a maintainer.
Before being merged, a change must pass
[CI](https://en.wikipedia.org/wiki/Continuous_integration).
#### Responsibilites of the approving reviewer
If you approve a change, it is expected that you:
- understand what the change is trying to do, and how it is doing it
- have manually built and tested the change, to verify it works as intended
- believe the change generally matches the idioms, formatting rules,
and overall coding style of the relevant repo
- are ready and able to help resolve any problems that may be introduced by
merging the change.
If a PR is made by a contributor who has write access to the repo in question,
they are responsible for merging/rebasing the PR after it has been approved;
otherwise it will be merged by the reviewer.
If a patch adds or modifies behaviour that is observable in the client,
the reviewer should build the patch and verify that it works as expected.
### After submitting your change
You've done all this, and submitted your patch. What now?
_Read other PRs_. If you're waiting for a review, it's likely that other
pull requests are waiting for review as well. This can be a good time
to go and take a look at what other work is happening in the project;
and if another PR has review comments, it might provide a clue to the
type of feedback you might expect.
_Patience_. As a general goal, we try to respond to all pull requests
within a few days, and to do preliminary review within a week, but we
don't always succeed. If you've opened a PR and haven't heard from
anyone, feel free to comment on it, or stop by the IRC channel, to ask
if anyone has had a chance to take a look. It's very possible that it's
been lost in the shuffle.
## Getting more involved
If you are participating in the xi-editor project, you may receive
additional privileges:
_Organization membership_: If you are regularly making contributions
to a xi project, in any of the forms outlined above, we will be happy to
add you to the xi-editor organization, which will give you the ability
to do things like add labels to issues and view active projects.
_Contributor_: If you are regularly making substantive contributions
to a specific xi project, we will be happy to add you as a contributor
to the repo in question. Contributors are encouraged to review and
approve changes, respond to issues, and generally help to maintain
the project in question.
_Maintainer_: If you are making substantive contributions to multiple
repos over an extended period, you are regularly reviewing the work of
other contributors, and you are actively participating in planning and
discussion, you may, (at the discretion of @raphlinus) be invited to
take on the role of _maintainer_. Maintainers are responsible for
coordinating the general direction of the project, resolving
architectural questions, and doing the day to day work of moving the
project forward.
find_package (Doxygen REQUIRED dot OPTIONAL_COMPONENTS mscgen dia)
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)
......@@ -28,4 +32,3 @@ else (DOXYGEN_FOUND)
message ("Doxygen need to be installed to generate the doxygen documentation")
endif (DOXYGEN_FOUND)
......@@ -14,10 +14,24 @@ add_executable (stack_example stack_example.cc)
target_link_libraries (stack_example SuperStupidStack CORSIKAunits CORSIKAlogging)
install (TARGETS stack_example DESTINATION share/examples)
add_executable (cascade_example cascade_example.cc)
target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogging
CORSIKArandom
CORSIKAsibyll
CORSIKAcascade
ProcessStackInspector
CORSIKAprocesses
CORSIKAparticles
CORSIKAgeometry
CORSIKAprocesssequence
)
install (TARGETS cascade_example DESTINATION share/examples)
add_executable (staticsequence_example staticsequence_example.cc)
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/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/random/RNGManager.h>
#include <corsika/cascade/SibStack.h>
#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/units/PhysicalUnits.h>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
#include <iostream>
#include <typeinfo>
using namespace std;
static int fCount = 0;
class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
public:
ProcessSplit() {}
template <typename Particle>
double MinStepLength(Particle& p, setup::Trajectory&) const {
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
int kBeam = 1;
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
// target nuclei: A < 18
// FOR NOW: assume target is oxygen
int kTarget = 1;
double beamEnergy = p.GetEnergy() / 1_GeV;
std::cout << "ProcessSplit: "
<< "MinStep: en: " << beamEnergy << " pid:" << kBeam << std::endl;
double prodCrossSection, dummy, dum1, dum2, dum3, dum4;
double dumdif[3];
if (kTarget == 1)
sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif, dum3, dum4);
else
sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy);
std::cout << "ProcessSplit: "
<< "MinStep: sibyll return: " << prodCrossSection << std::endl;
CrossSectionType sig = prodCrossSection * 1_mbarn;
std::cout << "ProcessSplit: "
<< "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
std::cout << "ProcessSplit: "
<< "nucleon mass " << nucleon_mass << std::endl;
// calculate interaction length in medium
double int_length = kTarget * (nucleon_mass / 1_g) / (sig / 1_cmeter / 1_cmeter);
// pick random step lenth
std::cout << "ProcessSplit: "
<< "interaction length (g/cm2): " << int_length << std::endl;
// add exponential sampling
int a = 0;
const double next_step = -int_length * log(s_rndm_(a));
/*
what are the units of the output? slant depth or 3space length?
*/
std::cout << "ProcessSplit: "
<< "next step (g/cm2): " << next_step << std::endl;
return next_step;
}
template <typename Particle, typename Stack>
EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
// corsika::utls::ignore(p);
return EProcessReturn::eOk;
}
template <typename Particle, typename Stack>
void DoDiscrete(Particle& p, Stack& s) const {
cout << "DoDiscrete: " << p.GetPID() << " interaction? "
<< process::sibyll::CanInteract(p.GetPID()) << endl;
if (process::sibyll::CanInteract(p.GetPID())) {
// get energy of particle from stack
/*
stack is in GeV in lab. frame
convert to GeV in cm. frame
(assuming proton at rest as target AND
assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
*/
EnergyType E = p.GetEnergy();
EnergyType Ecm = sqrt(2. * E * 0.93827_GeV);
int kBeam = process::sibyll::ConvertToSibyllRaw(p.GetPID());
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
// FOR NOW: set target to proton
int kTarget = 1; // p.GetPID();
std::cout << "ProcessSplit: "
<< " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV
<< std::endl;
if (E < 8.5_GeV || Ecm < 10_GeV) {
std::cout << "ProcessSplit: "
<< " DoDiscrete: dropping particle.." << std::endl;
p.Delete();
fCount++;
} else {
// Sibyll does not know about units..
double sqs = Ecm / 1_GeV;
// running sibyll, filling stack
sibyll_(kBeam, kTarget, sqs);
// running decays
// decsib_();
// print final state
int print_unit = 6;
sib_list_(print_unit);
// delete current particle
p.Delete();
// add particles from sibyll to stack
// link to sibyll stack
SibStack ss;
/*
get transformation between Stack-frame and SibStack-frame
for EAS Stack-frame is lab. frame, could be different for CRMC-mode
the transformation should be derived from the input momenta
in general transformation is rotation + boost
*/
const EnergyType proton_mass = 0.93827_GeV;
const double gamma = (E + proton_mass) / (Ecm);
const double gambet = sqrt(E * E - proton_mass * proton_mass) / Ecm;
// SibStack does not know about momentum yet so we need counter to access momentum
// array in Sibyll
int i = -1;
for (auto& p : ss) {
++i;
// transform to lab. frame, primitve
const double en_lab = gambet * s_plist_.p[2][i] + gamma * p.GetEnergy();
// add to corsika stack
auto pnew = s.NewParticle();
pnew.SetEnergy(en_lab * 1_GeV);
pnew.SetPID(process::sibyll::ConvertFromSibyll(p.GetPID()));
}
}
} else
p.Delete();
}
void Init() {
fCount = 0;
// define reference frame? --> defines boosts between corsika stack and model stack.
// initialize random numbers for sibyll
// FOR NOW USE SIBYLL INTERNAL !!!
// rnd_ini_();
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
;
const std::string str_name = "s_rndm";
rmng.RegisterRandomStream(str_name);
// // corsika::random::RNG srng;
// auto srng = rmng.GetRandomStream("s_rndm");
// test random number generator
std::cout << "ProcessSplit: "
<< " test sequence of random numbers." << std::endl;
int a = 0;
for (int i = 0; i < 8; ++i) std::cout << i << " " << s_rndm_(a) << std::endl;
// initialize Sibyll
sibyll_ini_();
// set particles stable / unstable
// use stack to loop over particles
setup::Stack ds;
ds.NewParticle().SetPID(Code::Proton);
ds.NewParticle().SetPID(Code::Neutron);
ds.NewParticle().SetPID(Code::PiPlus);
ds.NewParticle().SetPID(Code::PiMinus);
ds.NewParticle().SetPID(Code::KPlus);
ds.NewParticle().SetPID(Code::KMinus);
ds.NewParticle().SetPID(Code::K0Long);
ds.NewParticle().SetPID(Code::K0Short);
for (auto& p : ds) {
int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
// set particle stable by setting table value negative
cout << "ProcessSplit: Init: setting " << p.GetPID() << "(" << s_id << ")"
<< " stable in Sibyll .." << endl;
s_csydec_.idb[s_id] = -s_csydec_.idb[s_id - 1];
p.Delete();
}
}
int GetCount() { return fCount; }
private:
};
double s_rndm_(int&) {
static corsika::random::RNG& rmng =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
;
return rmng() / (double)rmng.max();
}
int main() {
tracking_line::TrackingLine<setup::Stack> tracking;
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
const auto sequence = p0 + p1;
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
stack.Clear();
auto particle = stack.NewParticle();
EnergyType E0 = 100_GeV;
particle.SetEnergy(E0);
particle.SetPID(Code::Proton);
EAS.Init();
EAS.Run();
cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
}
#include <corsika/geometry/CoordinateSystem.h>
/**
* (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>
......@@ -8,12 +19,14 @@
#include <iostream>
#include <typeinfo>
using namespace corsika;
using namespace corsika::geometry;
using namespace corsika::units::si;
int main() {
// define the root coordinate system
CoordinateSystem root;
geometry::CoordinateSystem& root =
geometry::RootCoordinateSystem::GetInstance().GetRootCS();
// another CS defined by a translation relative to the root CS
CoordinateSystem cs2 = root.translate({0_m, 0_m, 1_m});
......
#include <corsika/geometry/CoordinateSystem.h>
/**
* (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() {
CoordinateSystem root;
geometry::CoordinateSystem& root =
geometry::RootCoordinateSystem::GetInstance().GetRootCS();
Point const r0(root, {0_m, 0_m, 0_m});
auto const omegaC = 2 * M_PI * 1_Hz;
......
/**
* (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>
......
/**
* (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 <iomanip>
......
/**
* (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/setup/SetupTrajectory.h> // TODO: try to break this dependence later
using corsika::setup::Trajectory;
#include <corsika/units/PhysicalUnits.h> // dito
using namespace corsika::units::si;
using namespace std;
using namespace corsika::process;
......@@ -11,7 +27,7 @@ class Process1 : public BaseProcess<Process1> {
public:
Process1() {}
template <typename D, typename T, typename S>
EProcessReturn DoContinuous(D& d, T& t, S& s) const {
EProcessReturn DoContinuous(D& d, T&, S&) const {
for (int i = 0; i < 10; ++i) d.p[i] += 1;
return EProcessReturn::eOk;
}
......@@ -22,7 +38,7 @@ public:
Process2() {}
template <typename D, typename T, typename S>
inline EProcessReturn DoContinuous(D& d, T& t, S& s) const {
inline EProcessReturn DoContinuous(D&, T&, S&) const {
// for (int i=0; i<10; ++i) d.p[i] *= 2;
return EProcessReturn::eOk;
}
......@@ -34,7 +50,7 @@ public:
Process3() {}
template <typename D, typename T, typename S>
inline EProcessReturn DoContinuous(D& d, T& t, S& s) const {
inline EProcessReturn DoContinuous(D& /*d*/, T& /*t*/, S& /*s*/) const {
// for (int i=0; i<10; ++i) d.p[i] += fV;
return EProcessReturn::eOk;
}
......@@ -48,7 +64,7 @@ public:
// Process4(const int v) : fV(v) {}
Process4() {}
template <typename D, typename T, typename S>
inline EProcessReturn DoContinuous(D& d, T& t, S& s) const {
inline EProcessReturn DoContinuous(D& /*d*/, T& /*t*/, S& /*s*/) const {
// for (int i=0; i<10; ++i) d.p[i] /= fV;
return EProcessReturn::eOk;
}
......@@ -61,7 +77,6 @@ struct DummyData {
double p[10];
};
struct DummyStack {};
struct DummyTrajectory {};
void modular() {
......@@ -73,8 +88,8 @@ void modular() {
const auto sequence = m1 + m2 + m3 + m4;
DummyData p;
DummyTrajectory t;
DummyStack s;
Trajectory t;
const int n = 100000000;
for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
......
/**
* (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.
*/
add_subdirectory (Utilities)
add_subdirectory (Units)
add_subdirectory (Geometry)
add_subdirectory (Particles)
......
......@@ -9,16 +9,33 @@ set (
set (
CORSIKAcascade_HEADERS
Cascade.h
sibyll2.3c.h
SibStack.h
)
#set (
# CORSIKAcascade_SOURCES
# TrackingStep.cc
# Cascade.cc
# )
set (
CORSIKAsibyll_NAMESPACE
corsika/cascade
)
set (
CORSIKAsibyll_SOURCES
sibyll2.3c.f
gasdev.f
)
add_library (CORSIKAsibyll STATIC ${CORSIKAsibyll_SOURCES})
#add_library (CORSIKAcascade STATIC ${CORSIKAcascade_SOURCES})
add_library (CORSIKAcascade INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAcascade ${CORSIKAcascade_NAMESPACE} ${CORSIKAcascade_HEADERS})
#target_link_libraries (
......@@ -51,11 +68,16 @@ add_executable (
target_link_libraries (
testCascade
# CORSIKAutls
CORSIKArandom
CORSIKAsibyll
CORSIKAcascade
ProcessStackInspector
CORSIKAstackinterface
CORSIKAprocesses
CORSIKAparticles
CORSIKAgeometry
CORSIKAprocesssequence
SuperStupidStack
CORSIKAunits
CORSIKAthirdparty # for catch2
)
......
/**
* (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>
using namespace corsika::cascade;
template <typename ProcessList, typename Particle, typename Trajectory, typename Stack>
Cascade<ProcessList, Particle, Trajectory, Stack>::Cascade() {
// kkk;
// kk;
}
template <typename ProcessList, typename Particle, typename Trajectory, typename Stack>
template <typename ProcessList, typename Stack> //, typename Trajectory>
void Cascade::Init() {
fStack.Init();
fProcesseList.Init();
}
template <typename ProcessList, typename Particle, typename Trajectory, typename Stack>
template <typename ProcessList, typename Stack> //, typename Trajectory>
void Cascade::Run() {
if (!fStack.IsEmpty()) {
if (!fStack.IsEmpty()) {
......@@ -27,10 +32,12 @@ void Cascade::Run() {
}
}
template <typename Sequence, typename Trajectory>
template <typename ProcessList, typename Stack> //, typename Trajectory>
void Cascade::Step(Particle& particle) {
double nextStep = fProcesseList.MinStepLength(particle);
Trajectory trajectory = fProcesseList.Transport(particle, nextStep);
corsika::geometry::LineTrajectory trajectory =
fProcesseList.Transport(particle, nextStep);
sequence.DoContinuous(particle, trajectory);
// whats going on here? Everywhere else DoDiscrete is passed a Stack reference as well
sequence.DoDiscrete(particle);
}
#ifndef _include_Cascade_h_
#define _include_Cascade_h_
#include <corsika/geometry/LineTrajectory.h> // to be removed
#include <corsika/geometry/Point.h> // to be removed
/**
* (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_corsika_cascade_Cascade_h_
#define _include_corsika_cascade_Cascade_h_
#include <corsika/process/ProcessReturn.h>
#include <corsika/units/PhysicalUnits.h>
#include <type_traits>
#include <corsika/setup/SetupTrajectory.h>
using namespace corsika::units::si;
namespace corsika::cascade {
template <typename Trajectory, typename ProcessList, typename Stack>
template <typename Tracking, typename ProcessList, typename Stack>
class Cascade {
typedef typename Stack::ParticleType Particle;
Cascade() = delete;
public:
Cascade(ProcessList& pl, Stack& stack)
: fProcesseList(pl)
, fStack(stack) {}
Cascade(Tracking& tr, ProcessList& pl, Stack& stack)
: fTracking(tr)
, fProcesseList(pl)
, fStack(stack) {
// static_assert(std::is_member_function_pointer<decltype(&ProcessList::DoDiscrete)>::value,
//"ProcessList has not function DoDiscrete.");
// static_assert(std::is_member_function_pointer<decltype(&ProcessList::DoContinuous)>::value,
// "ProcessList has not function DoContinuous.");
}
void Init() {
fStack.Init();
fTracking.Init();
fProcesseList.Init();
fStack.Init();
}
void Run() {
while (!fStack.IsEmpty()) {
while (!fStack.IsEmpty()) {
// Particle& p = *fStack.GetNextParticle();
EnergyType Emin;
typename Stack::StackIterator pMin(fStack, 0);
bool first = true;
for (typename Stack::StackIterator ip = fStack.begin(); ip != fStack.end();
++ip) {
if (first || ip.GetEnergy() < Emin) {
first = false;
pMin = ip;
Emin = pMin.GetEnergy();
}
}
Step(pMin);
Particle& pNext = *fStack.GetNextParticle();
Step(pNext);
}
// do cascade equations, which can put new particles on Stack,
// thus, the double loop
......@@ -50,24 +60,25 @@ namespace corsika::cascade {
}
void Step(Particle& particle) {
double nextStep = fProcesseList.MinStepLength(particle);
corsika::geometry::CoordinateSystem root;
Trajectory trajectory(
corsika::geometry::Point(root, {0_m, 0_m, 0_m}),
corsika::geometry::Vector<corsika::units::si::SpeedType::dimension_type>(
root, 0 * 1_m / second, 0 * 1_m / second, 1 * 1_m / second));
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
fProcesseList.MinStepLength(particle, step);
/// here the particle is actually moved along the trajectory to new position:
std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
corsika::process::EProcessReturn status =
fProcesseList.DoContinuous(particle, trajectory, fStack);
fProcesseList.DoContinuous(particle, step, fStack);
if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
fStack.Delete(particle);
fStack.Delete(particle); // TODO: check if this is really needed
} else {
fProcesseList.DoDiscrete(particle, fStack);
}
}
private:
Stack& fStack;
Tracking& fTracking;
ProcessList& fProcesseList;
Stack& fStack;
};
} // namespace corsika::cascade
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment