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
Hide whitespace changes
Inline Side-by-side
Showing
with 178 additions and 1689 deletions
Processes/StackInspector/testStackInspector.cc deleted 100644 → 0
View file @ 6f2ab230
/*
* (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.
*/
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/cascade/testCascade.h>
using namespace corsika::units::si;
using namespace corsika::process::stack_inspector;
using namespace corsika;
using namespace corsika::geometry;
TEST_CASE("StackInspector", "[processes]") {
auto const& rootCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
geometry::Point const origin(rootCS, {0_m, 0_m, 0_m});
geometry::Vector<units::si::SpeedType::dimension_type> v(rootCS, 0_m / second,
0_m / second, 1_m / second);
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
TestCascadeStack stack;
stack.Clear();
HEPEnergyType E0 = 100_GeV;
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Electron, E0,
corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
Point(rootCS, {0_m, 0_m, 10_km}), 0_ns});
SECTION("interface") {
StackInspector<TestCascadeStack> model(1, true);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret = model.DoStack(stack);
}
}
Processes/TrackWriter/CMakeLists.txt deleted 100644 → 0
View file @ 6f2ab230
set (
MODEL_SOURCES
TrackWriter.cc
)
set (
MODEL_HEADERS
TrackWriter.h
)
set (
MODEL_NAMESPACE
corsika/process/track_writer
)
add_library (ProcessTrackWriter STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackWriter ${MODEL_NAMESPACE} ${MODEL_HEADERS})
set_target_properties (
ProcessTrackWriter
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
# PUBLIC_HEADER "${MODEL_HEADERS}"
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessTrackWriter
CORSIKAunits
CORSIKAparticles
CORSIKAgeometry
CORSIKAsetup
)
target_include_directories (
ProcessTrackWriter
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
TARGETS ProcessTrackWriter
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
)
# --------------------
# code unit testing
#add_executable (testNullModel testNullModel.cc)
#target_link_libraries (
# testNullModel ProcessNullModel
# CORSIKAsetup
# CORSIKAgeometry
# CORSIKAunits
# CORSIKAthirdparty # for catch2
# )
# CORSIKA_ADD_TEST(testNullModel)
Processes/TrackWriter/TrackWriter.cc deleted 100644 → 0
View file @ 6f2ab230
/*
* (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/process/track_writer/TrackWriter.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <limits>
using namespace corsika::setup;
using Particle = Stack::ParticleType;
using Track = Trajectory;
namespace corsika::process::TrackWriter {
void TrackWriter::Init() {
using namespace std::string_literals;
fFile.open(fFilename);
fFile << "# PID, E / eV, start coordinates / m, displacement vector to end / m "s
<< '\n';
}
template <>
process::EProcessReturn TrackWriter::DoContinuous(Particle& vP, Track& vT) {
using namespace units::si;
auto const start = vT.GetPosition(0).GetCoordinates();
auto const delta = vT.GetPosition(1).GetCoordinates() - start;
auto const& name = particles::GetName(vP.GetPID());
fFile << name << " " << vP.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
<< start[1] / 1_m << ' ' << start[2] / 1_m << " " << delta[0] / 1_m << ' '
<< delta[1] / 1_m << ' ' << delta[2] / 1_m << '\n';
return process::EProcessReturn::eOk;
}
template <>
units::si::LengthType TrackWriter::MaxStepLength(Particle&, Track&) {
return units::si::meter * std::numeric_limits<double>::infinity();
}
} // namespace corsika::process::TrackWriter
Processes/TrackWriter/TrackWriter.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (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 _Processes_TrackWriter_h_
#define _Processes_TrackWriter_h_
#include <corsika/process/ContinuousProcess.h>
#include <corsika/units/PhysicalUnits.h>
#include <fstream>
#include <string>
namespace corsika::process::TrackWriter {
class TrackWriter : public corsika::process::ContinuousProcess<TrackWriter> {
public:
TrackWriter(std::string const& filename)
: fFilename(filename) {}
void Init();
template <typename Particle, typename Track>
corsika::process::EProcessReturn DoContinuous(Particle&, Track&);
template <typename Particle, typename Track>
corsika::units::si::LengthType MaxStepLength(Particle&, Track&);
private:
std::string const fFilename;
std::ofstream fFile;
};
} // namespace corsika::process::TrackWriter
#endif
Processes/TrackingLine/CMakeLists.txt deleted 100644 → 0
View file @ 6f2ab230
set (
MODEL_HEADERS
TrackingLine.h
)
set (
MODEL_SOURCES
TrackingLine.cc
)
set (
MODEL_NAMESPACE
corsika/process/tracking_line
)
add_library (ProcessTrackingLine STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackingLine ${MODEL_NAMESPACE} ${MODEL_HEADERS})
set_target_properties (
ProcessTrackingLine
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
# PUBLIC_HEADER "${MODEL_HEADERS}"
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessTrackingLine
CORSIKAsetup
CORSIKAutilities
CORSIKAenvironment
CORSIKAunits
CORSIKAenvironment
CORSIKAgeometry
)
target_include_directories (
ProcessTrackingLine
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
# #-- -- -- -- -- -- -- -- -- --
# #code unit testing
add_executable (testTrackingLine testTrackingLine.cc)
target_link_libraries (
testTrackingLine
ProcessTrackingLine
CORSIKAthirdparty # for catch2
)
CORSIKA_ADD_TEST(testTrackingLine)
Processes/TrackingLine/TrackingLine.cc deleted 100644 → 0
View file @ 6f2ab230
/*
* (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/environment/Environment.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/QuantityVector.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <algorithm>
#include <iostream>
#include <stdexcept>
#include <utility>
using namespace corsika;
namespace corsika::process::tracking_line {
std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
TimeOfIntersection(corsika::geometry::Line const& line,
geometry::Sphere const& sphere) {
auto const delta = line.GetR0() - sphere.GetCenter();
auto const v = line.GetV0();
auto const vSqNorm = v.squaredNorm();
auto const R = sphere.GetRadius();
auto const vDotDelta = v.dot(delta);
auto const discriminant =
vDotDelta * vDotDelta - vSqNorm * (delta.squaredNorm() - R * R);
//~ std::cout << "discriminant: " << discriminant << std::endl;
//~ std::cout << "alpha: " << alpha << std::endl;
//~ std::cout << "beta: " << beta << std::endl;
if (discriminant.magnitude() > 0) {
auto const sqDisc = sqrt(discriminant);
auto const invDenom = 1 / vSqNorm;
return std::make_pair((-vDotDelta - sqDisc) * invDenom,
(-vDotDelta + sqDisc) * invDenom);
} else {
return {};
}
}
} // namespace corsika::process::tracking_line
Processes/TrackingLine/TrackingLine.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (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_processes_TrackingLine_h_
#define _include_corsika_processes_TrackingLine_h_
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <optional>
#include <type_traits>
#include <utility>
namespace corsika::environment {
template <typename IEnvironmentModel>
class Environment;
template <typename IEnvironmentModel>
class VolumeTreeNode;
} // namespace corsika::environment
namespace corsika::process {
namespace tracking_line {
std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
TimeOfIntersection(corsika::geometry::Line const& line,
geometry::Sphere const& sphere);
class TrackingLine {
public:
TrackingLine(){};
template <typename Particle> // was Stack previously, and argument was
// Stack::StackIterator
auto GetTrack(Particle const& p) {
using namespace corsika::units::si;
using namespace corsika::geometry;
geometry::Vector<SpeedType::dimension_type> const velocity =
p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
auto const currentPosition = p.GetPosition();
std::cout << "TrackingLine pid: " << p.GetPID()
<< " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
std::cout << "TrackingLine pos: "
<< currentPosition.GetCoordinates()
// << " [" << p.GetNode()->GetModelProperties().GetName() << "]"
<< std::endl;
std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
<< " GeV " << std::endl;
std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
// to do: include effect of magnetic field
geometry::Line line(currentPosition, velocity);
auto const* currentLogicalVolumeNode = p.GetNode();
//~ auto const* currentNumericalVolumeNode =
//~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
auto const numericallyInside =
currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false");
auto const& children = currentLogicalVolumeNode->GetChildNodes();
auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
// for entering from outside
auto addIfIntersects = [&](auto const& vtn) {
auto const& volume = vtn.GetVolume();
auto const& sphere = dynamic_cast<geometry::Sphere const&>(
volume); // for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
auto const [t1, t2] = *opt;
std::cout << "intersection times: " << t1 / 1_s << "; "
<< t2 / 1_s
// << " " << vtn.GetModelProperties().GetName()
<< std::endl;
if (t1.magnitude() > 0)
intersections.emplace_back(t1, &vtn);
else if (t2.magnitude() > 0)
std::cout << "inside other volume" << std::endl;
}
};
for (auto const& child : children) { addIfIntersects(*child); }
for (auto const* ex : excluded) { addIfIntersects(*ex); }
{
auto const& sphere = dynamic_cast<geometry::Sphere const&>(
currentLogicalVolumeNode->GetVolume());
// for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
[[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line, sphere);
[[maybe_unused]] auto dummy_t1 = t1;
intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent());
}
auto const minIter = std::min_element(
intersections.cbegin(), intersections.cend(),
[](auto const& a, auto const& b) { return a.first < b.first; });
TimeType min;
if (minIter == intersections.cend()) {
min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
// to handle the numerics properly
throw std::runtime_error("no intersection with anything!");
} else {
min = minIter->first;
}
std::cout << " t-intersect: "
<< min
// << " " << minIter->second->GetModelProperties().GetName()
<< std::endl;
return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
velocity.norm() * min, minIter->second);
}
};
} // namespace tracking_line
} // namespace corsika::process
#endif
Processes/TrackingLine/testTrackingLine.cc deleted 100644 → 0
View file @ 6f2ab230
/*
* (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/process/tracking_line/TrackingLine.h>
#include <testTrackingLineStack.h> // test-build, and include file is obtained from CMAKE_CURRENT_SOURCE_DIR
#include <corsika/environment/Environment.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
#include <corsika/setup/SetupTrajectory.h>
using corsika::setup::Trajectory;
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::geometry;
#include <iostream>
using namespace std;
using namespace corsika::units::si;
TEST_CASE("TrackingLine") {
environment::Environment<environment::Empty> env; // dummy environment
auto const& cs = env.GetCoordinateSystem();
tracking_line::TrackingLine tracking;
SECTION("intersection with sphere") {
Point const origin(cs, {0_m, 0_m, -5_m});
Point const center(cs, {0_m, 0_m, 10_m});
Sphere const sphere(center, 1_m);
Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
0_m / second, 1_m / second);
Line line(origin, v);
setup::Trajectory traj(line, 12345_s);
auto const opt =
tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
REQUIRE(opt.has_value());
auto [t1, t2] = opt.value();
REQUIRE(t1 / 14_s == Approx(1));
REQUIRE(t2 / 16_s == Approx(1));
auto const optNoIntersection =
tracking_line::TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
REQUIRE_FALSE(optNoIntersection.has_value());
}
SECTION("maximally possible propagation") {
auto& universe = *(env.GetUniverse());
auto const radius = 20_m;
auto theMedium = environment::Environment<environment::Empty>::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
auto const* theMediumPtr = theMedium.get();
universe.AddChild(std::move(theMedium));
TestTrackingLineStack stack;
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::MuPlus,
1_GeV,
{cs, {0_GeV, 0_GeV, 1_GeV}},
{cs, {0_m, 0_m, 0_km}},
0_ns});
auto p = stack.GetNextParticle();
p.SetNode(theMediumPtr);
Point const origin(cs, {0_m, 0_m, 0_m});
Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
0_m / second, 1_m / second);
Line line(origin, v);
auto const [traj, geomMaxLength, nextVol] = tracking.GetTrack(p);
[[maybe_unused]] auto dummy_geomMaxLength = geomMaxLength;
[[maybe_unused]] auto dummy_nextVol = nextVol;
REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius))
.GetComponents(cs)
.norm()
.magnitude() == Approx(0).margin(1e-4));
}
}
Processes/TrackingLine/testTrackingLineStack.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (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_process_trackinling_teststack_h_
#define _include_process_trackinling_teststack_h_
#include <corsika/environment/Environment.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/units/PhysicalUnits.h>
using TestEnvironmentType =
corsika::environment::Environment<corsika::environment::Empty>;
template <typename T>
using SetupGeometryDataInterface = GeometryDataInterface<T, TestEnvironmentType>;
// combine particle data stack with geometry information for tracking
template <typename StackIter>
using StackWithGeometryInterface = corsika::stack::CombinedParticleInterface<
corsika::setup::detail::ParticleDataStack::PIType, SetupGeometryDataInterface,
StackIter>;
using TestTrackingLineStack = corsika::stack::CombinedStack<
typename corsika::setup::detail::ParticleDataStack::StackImpl,
GeometryData<TestEnvironmentType>, StackWithGeometryInterface>;
#endif
Processes/UrQMD/CMakeLists.txt deleted 100644 → 0
View file @ 6f2ab230
set (
MODEL_SOURCES
UrQMD.cc
urqmdInterface.F
addpart.f
angdis.f
anndec.f
blockres.f
boxprg.f
cascinit.f
coload.f
dectim.f
delpart.f
detbal.f
dwidth.f
error.f
getmass.f
getspin.f
init.f
iso.f
ityp2pdg.f
jdecay2.f
make22.f
numrec.f
output.f
paulibl.f
proppot.f
saveinfo.f
scatter.f
siglookup.f
string.f
tabinit.f
urqmd.f
whichres.f
)
set (
MODEL_HEADERS
UrQMD.h
)
set (
MODEL_NAMESPACE
corsika/process/urqmd
)
add_library (ProcessUrQMD STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessUrQMD ${MODEL_NAMESPACE} ${MODEL_HEADERS})
set_target_properties (
ProcessUrQMD
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
# PUBLIC_HEADER "${MODEL_HEADERS}"
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessUrQMD
CORSIKAprocesssequence
CORSIKAparticles
CORSIKAunits
CORSIKAgeometry
CORSIKArandom
CORSIKAsetup
)
target_include_directories (
ProcessUrQMD
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
TARGETS ProcessUrQMD
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
)
# --------------------
# code unit testing
add_executable (testUrQMD
testUrQMD.cc
${MODEL_HEADERS}
)
target_link_libraries (
testUrQMD
ProcessUrQMD
CORSIKAsetup
CORSIKArandom
CORSIKAgeometry
CORSIKAunits
CORSIKAthirdparty # for catch2
)
CORSIKA_ADD_TEST(testUrQMD)
Processes/UrQMD/UrQMD.cc deleted 100644 → 0
View file @ 6f2ab230
#include <corsika/geometry/QuantityVector.h>
#include <corsika/geometry/Vector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/urqmd/UrQMD.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <algorithm>
#include <functional>
#include <iostream>
#include <random>
using namespace corsika::process::UrQMD;
using namespace corsika::units::si;
UrQMD::UrQMD() { iniurqmd_(); }
using SetupStack = corsika::setup::Stack;
using SetupParticle = corsika::setup::Stack::StackIterator;
using SetupProjectile = corsika::setup::StackView::StackIterator;
using SetupTrack = corsika::setup::Trajectory;
template <typename TParticle> // need template here, as this is called both with
// SetupParticle as well as SetupProjectile
CrossSectionType UrQMD::GetCrossSection(TParticle const& vProjectile,
corsika::particles::Code vTargetCode) const {
using namespace units::si;
// TODO: energy cuts, return 0 for non-hadrons
auto const projectileCode = vProjectile.GetPID();
auto const projectileEnergyLab = vProjectile.GetEnergy();
// the following is a translation of ptsigtot() into C++
if (projectileCode != particles::Code::Nucleus &&
!IsNucleus(vTargetCode)) { // both particles are "special"
auto const mProj = particles::GetMass(projectileCode);
auto const mTar = particles::GetMass(vTargetCode);
double sqrtS =
sqrt(units::si::detail::static_pow<2>(mProj) +
units::si::detail::static_pow<2>(mTar) + 2 * projectileEnergyLab * mTar) *
(1 / 1_GeV);
// we must set some UrQMD globals first...
auto const [ityp, iso3] = ConvertToUrQMD(projectileCode);
inputs_.spityp[0] = ityp;
inputs_.spiso3[0] = iso3;
auto const [itypTar, iso3Tar] = ConvertToUrQMD(vTargetCode);
inputs_.spityp[1] = itypTar;
inputs_.spiso3[1] = iso3Tar;
int one = 1;
int two = 2;
return sigtot_(one, two, sqrtS) * 1_mb;
} else {
int const Ap =
(projectileCode == particles::Code::Nucleus) ? vProjectile.GetNuclearA() : 1;
int const At = IsNucleus(vTargetCode) ? particles::GetNucleusA(vTargetCode) : 1;
double const maxImpact = nucrad_(Ap) + nucrad_(At) + 2 * options_.CTParam[30 - 1];
return 10_mb * M_PI * units::si::detail::static_pow<2>(maxImpact);
// is a constant cross-section really reasonable?
}
}
bool UrQMD::CanInteract(particles::Code vCode) const {
// According to the manual, UrQMD can use all mesons, baryons and nucleons
// which are modeled also as input particles. I think it is safer to accept
// only the usual long-lived species as input.
// TODO: Charmed mesons should be added to the list, too
static particles::Code const validProjectileCodes[] = {
particles::Code::Nucleus, particles::Code::Proton, particles::Code::AntiProton,
particles::Code::Neutron, particles::Code::AntiNeutron, particles::Code::PiPlus,
particles::Code::PiMinus, particles::Code::KPlus, particles::Code::KMinus,
particles::Code::K0, particles::Code::K0Bar};
return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
vCode) != std::cend(validProjectileCodes);
}
GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle, SetupTrack&) const {
if (!CanInteract(vParticle.GetPID())) {
// we could do the canInteract check in GetCrossSection, too but if
// we do it here we have the advantage of avoiding the loop
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
auto const& mediumComposition =
vParticle.GetNode()->GetModelProperties().GetNuclearComposition();
using namespace std::placeholders;
CrossSectionType const weightedProdCrossSection = mediumComposition.WeightedSum(
std::bind(&UrQMD::GetCrossSection<decltype(vParticle)>, this, vParticle, _1));
return mediumComposition.GetAverageMassNumber() * units::constants::u /
weightedProdCrossSection;
}
corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjectile) {
using namespace units::si;
auto const projectileCode = vProjectile.GetPID();
auto const projectileEnergyLab = vProjectile.GetEnergy();
auto const& projectileMomentumLab = vProjectile.GetMomentum();
auto const& projectilePosition = vProjectile.GetPosition();
auto const projectileTime = vProjectile.GetTime();
// sample target particle
auto const& mediumComposition =
vProjectile.GetNode()->GetModelProperties().GetNuclearComposition();
auto const componentCrossSections = std::invoke([&]() {
auto const& components = mediumComposition.GetComponents();
std::vector<CrossSectionType> crossSections;
crossSections.reserve(components.size());
for (auto const c : components) {
crossSections.push_back(GetCrossSection(vProjectile, c));
}
return crossSections;
});
auto const targetCode = mediumComposition.SampleTarget(componentCrossSections, fRNG);
auto const targetA = particles::GetNucleusA(targetCode);
auto const targetZ = particles::GetNucleusZ(targetCode);
inputs_.nevents = 1;
sys_.eos = 0; // could be configurable in principle
inputs_.outsteps = 1;
sys_.nsteps = 1;
// initialization regarding projectile
if (particles::Code::Nucleus == projectileCode) {
// is this everything?
inputs_.prspflg = 0;
sys_.Ap = vProjectile.GetNuclearA();
sys_.Zp = vProjectile.GetNuclearZ();
rsys_.ebeam = (projectileEnergyLab - vProjectile.GetMass()) * (1 / 1_GeV) /
vProjectile.GetNuclearA();
rsys_.bdist = nucrad_(targetA) + nucrad_(sys_.Ap) + 2 * options_.CTParam[30 - 1];
int const id = 1;
cascinit_(sys_.Zp, sys_.Ap, id);
} else {
inputs_.prspflg = 1;
sys_.Ap = 1; // even for non-baryons this has to be set, see vanilla UrQMD.f
rsys_.bdist = nucrad_(targetA) + nucrad_(1) + 2 * options_.CTParam[30 - 1];
rsys_.ebeam = (projectileEnergyLab - vProjectile.GetMass()) * (1 / 1_GeV);
auto const [ityp, iso3] = ConvertToUrQMD(projectileCode);
// todo: conversion of K_long/short into strong eigenstates;
inputs_.spityp[0] = ityp;
inputs_.spiso3[0] = iso3;
}
// initilazation regarding target
if (particles::IsNucleus(targetCode)) {
sys_.Zt = targetZ;
sys_.At = targetA;
inputs_.trspflg = 0; // nucleus as target
int const id = 2;
cascinit_(sys_.Zt, sys_.At, id);
} else {
inputs_.trspflg = 1; // special particle as target
auto const [ityp, iso3] = ConvertToUrQMD(targetCode);
inputs_.spityp[1] = ityp;
inputs_.spiso3[1] = iso3;
}
int iflb = 0; // flag for retrying interaction in case of empty event, 0 means retry
urqmd_(iflb);
// now retrieve secondaries from UrQMD
auto const& originalCS = projectileMomentumLab.GetCoordinateSystem();
geometry::CoordinateSystem const zAxisFrame =
originalCS.RotateToZ(projectileMomentumLab);
for (int i = 0; i < sys_.npart; ++i) {
auto const code = ConvertFromUrQMD(isys_.ityp[i], isys_.iso3[i]);
// "coor_.p0[i] * 1_GeV" is likely off-shell as UrQMD doesn't preserve masses well
auto momentum = geometry::Vector(
zAxisFrame,
geometry::QuantityVector<dimensionless_d>{coor_.px[i], coor_.py[i], coor_.pz[i]} *
1_GeV);
auto const energy = sqrt(momentum.squaredNorm() + square(particles::GetMass(code)));
momentum.rebase(originalCS); // transform back into standard lab frame
std::cout << i << " " << code << " " << momentum.GetComponents() << std::endl;
vProjectile.AddSecondary(
std::tuple<particles::Code, HEPEnergyType, stack::MomentumVector, geometry::Point,
TimeType>{code, energy, momentum, projectilePosition, projectileTime});
}
std::cout << "UrQMD generated " << sys_.npart << " secondaries!" << std::endl;
return process::EProcessReturn::eOk;
}
/**
* the random number generator function of UrQMD
*/
double corsika::process::UrQMD::ranf_(int&) {
static corsika::random::RNG& rng =
corsika::random::RNGManager::GetInstance().GetRandomStream("UrQMD");
static std::uniform_real_distribution<double> dist;
return dist(rng);
}
corsika::particles::Code corsika::process::UrQMD::ConvertFromUrQMD(int vItyp, int vIso3) {
int const pdgInt =
pdgid_(vItyp, vIso3); // use the conversion function provided by UrQMD
if (pdgInt == 0) { // pdgid_ returns 0 on error
throw std::runtime_error("UrQMD pdgid() returned 0");
}
auto const pdg = static_cast<particles::PDGCode>(pdgInt);
return particles::ConvertFromPDG(pdg);
}
std::pair<int, int> corsika::process::UrQMD::ConvertToUrQMD(
corsika::particles::Code code) {
static const std::map<int, std::pair<int, int>> mapPDGToUrQMD{
// data mostly from github.com/afedynitch/ParticleDataTool
{22, {100, 0}}, // gamma
{111, {101, 0}}, // pi0
{211, {101, 2}}, // pi+
{-211, {101, -2}}, // pi-
{321, {106, 1}}, // K+
{-321, {-106, -1}}, // K-
{311, {106, -1}}, // K0
{-311, {-106, 1}}, // K0bar
{2212, {1, 1}}, // p
{2112, {1, -1}}, // n
{221, {102, 0}}, // eta
{213, {104, 2}}, // rho+
{-213, {104, -2}}, // rho-
{113, {104, 0}}, // rho0
{323, {108, 2}}, // K*+
{-323, {108, -2}}, // K*-
{313, {108, 0}}, // K*0
{-313, {-108, 0}}, // K*0-bar
{223, {103, 0}}, // omega
{333, {109, 0}}, // phi
{3222, {40, 2}}, // Sigma+
{3212, {40, 0}}, // Sigma0
{3112, {40, -2}}, // Sigma-
{3322, {49, 0}}, // Xi0
{3312, {49, -1}}, // Xi-
{3122, {27, 0}}, // Lambda0
{2224, {17, 4}}, // Delta++
{2214, {17, 2}}, // Delta+
{2114, {17, 0}}, // Delta0
{1114, {17, -2}}, // Delta-
{3224, {41, 2}}, // Sigma*+
{3214, {41, 0}}, // Sigma*0
{3114, {41, -2}}, // Sigma*-
{3324, {50, 0}}, // Xi*0
{3314, {50, -1}}, // Xi*-
{3334, {55, 0}}, // Omega-
{411, {133, 2}}, // D+
{-411, {133, -2}}, // D-
{421, {133, 0}}, // D0
{-421, {-133, 0}}, // D0-bar
{441, {107, 0}}, // etaC
{431, {138, 1}}, // Ds+
{-431, {138, -1}}, // Ds-
{433, {139, 1}}, // Ds*+
{-433, {139, -1}}, // Ds*-
{413, {134, 1}}, // D*+
{-413, {134, -1}}, // D*-
{10421, {134, 0}}, // D*0
{-10421, {-134, 0}}, // D*0-bar
{443, {135, 0}}, // jpsi
};
return mapPDGToUrQMD.at(static_cast<int>(GetPDG(code)));
}
Processes/UrQMD/UrQMD.h deleted 100644 → 0
View file @ 6f2ab230
#ifndef _Processes_UrQMD_UrQMD_h
#define _Processes_UrQMD_UrQMD_h
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/random/RNGManager.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <array>
#include <utility>
namespace corsika::process::UrQMD {
class UrQMD : public corsika::process::InteractionProcess<UrQMD> {
public:
UrQMD();
corsika::units::si::GrammageType GetInteractionLength(
corsika::setup::Stack::StackIterator&, corsika::setup::Trajectory&) const;
template <typename TParticle>
corsika::units::si::CrossSectionType GetCrossSection(TParticle const&,
corsika::particles::Code) const;
corsika::process::EProcessReturn DoInteraction(
corsika::setup::StackView::StackIterator&);
bool CanInteract(particles::Code) const;
private:
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("UrQMD");
};
namespace constants {
// from coms.f
int constexpr nmax = 500;
// from options.f
int constexpr numcto = 400;
int constexpr numctp = 400;
// from inputs.f
int constexpr aamax = 300;
} // namespace constants
template <typename T>
using nmaxArray = std::array<T, constants::nmax>;
using nmaxIntArray = nmaxArray<int>;
using nmaxDoubleArray = nmaxArray<double>;
extern "C" {
void iniurqmd_();
double ranf_(int&);
void cascinit_(int const&, int const&, int const&);
double nucrad_(int const&);
void urqmd_(int&);
int pdgid_(int&, int&);
double sigtot_(int&, int&, double&);
// defined in coms.f
extern struct {
int npart, nbar, nmes, ctag, nsteps, uid_cnt, ranseed, event;
int Ap; // projectile mass number (in case of nucleus)
int At; // target mass number (in case of nucleus)
int Zp; // projectile charge number (in case of nucleus)
int Zt; // target charge number (in case of nucleus)
int eos, dectag, NHardRes, NSoftRes, NDecRes, NElColl, NBlColl;
} sys_;
extern struct {
double time, acttime, bdist, bimp, bmin;
double ebeam; // lab-frame energy of projectile
double ecm;
} rsys_;
// defined in coms.f
extern struct {
nmaxIntArray spin, ncoll, charge, ityp, lstcoll, iso3, origin, strid, uid;
} isys_;
// defined in coor.f
extern struct {
nmaxDoubleArray r0, rx, ry, rz, p0, px, py, pz, fmass, rww, dectime;
} coor_;
// defined in inputs.f
extern struct {
int nevents;
std::array<int, 2> spityp; // particle codes of: [0]: projectile, [1]: target
int prspflg; // projectile special flag
int trspflg; // target special flag, set to 1 unless target is nucleus > H
std::array<int, 2> spiso3; // particle codes of: [0]: projectile, [1]: target
int outsteps, bflag, srtflag, efuncflag, nsrt, npb, firstev;
} inputs_;
// defined in inputs.f
extern struct {
double srtmin, srtmax, pbeam, betann, betatar, betapro, pbmin, pbmax;
} input2_;
// defined in options.f
extern struct {
std::array<double, constants::numcto> CTOption;
std::array<double, constants::numctp> CTParam;
} options_;
extern struct {
int fixedseed, bf13, bf14, bf15, bf16, bf17, bf18, bf19, bf20;
} loptions_;
// defined in urqmdInterface.F
extern struct { std::array<double, 3> xs, bim; } cxs_u2_;
}
/**
* convert CORSIKA code to UrQMD code tuple
*
* In the current implementation a detour via the PDG code is made.
*/
std::pair<int, int> ConvertToUrQMD(particles::Code);
particles::Code ConvertFromUrQMD(int vItyp, int vIso3);
} // namespace corsika::process::UrQMD
#endif
Processes/UrQMD/testUrQMD.cc deleted 100644 → 0
View file @ 6f2ab230
/*
* (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/process/urqmd/UrQMD.h>
#include <corsika/random/RNGManager.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/CorsikaFenv.h>
#include <corsika/particles/ParticleProperties.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 <tuple>
#include <utility>
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::process::UrQMD;
using namespace corsika::units::si;
template <typename TStackView>
auto sumCharge(TStackView const& view) {
int totalCharge = 0;
for (auto const& p : view) { totalCharge += particles::GetChargeNumber(p.GetPID()); }
return totalCharge;
}
template <typename TStackView>
auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS) {
geometry::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
for (auto const& p : view) { sum += p.GetMomentum(); }
return sum;
}
auto setupEnvironment(particles::Code vTargetCode) {
// setup environment, geometry
auto env = std::make_unique<environment::Environment<environment::IMediumModel>>();
auto& universe = *(env->GetUniverse());
const geometry::CoordinateSystem& cs = env->GetCoordinateSystem();
auto theMedium =
environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
geometry::Point{cs, 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(std::vector<particles::Code>{vTargetCode},
std::vector<float>{1.}));
auto const* nodePtr = theMedium.get();
universe.AddChild(std::move(theMedium));
return std::make_tuple(std::move(env), &cs, nodePtr);
}
template <typename TNodeType>
auto setupStack(int vA, int vZ, HEPEnergyType vMomentum, TNodeType* vNodePtr,
geometry::CoordinateSystem const& cs) {
auto stack = std::make_unique<setup::Stack>();
auto constexpr mN = corsika::units::constants::nucleonMass;
geometry::Point const origin(cs, {0_m, 0_m, 0_m});
corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
HEPEnergyType const E0 =
sqrt(units::si::detail::static_pow<2>(mN * vA) + pLab.squaredNorm());
auto particle =
stack->AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, E0, pLab, origin, 0_ns, vA, vZ});
particle.SetNode(vNodePtr);
return std::make_tuple(
std::move(stack),
std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
}
template <typename TNodeType>
auto setupStack(particles::Code vProjectileType, HEPEnergyType vMomentum,
TNodeType* vNodePtr, geometry::CoordinateSystem const& cs) {
auto stack = std::make_unique<setup::Stack>();
geometry::Point const origin(cs, {0_m, 0_m, 0_m});
corsika::stack::MomentumVector const pLab(cs, {vMomentum, 0_GeV, 0_GeV});
HEPEnergyType const E0 =
sqrt(units::si::detail::static_pow<2>(particles::GetMass(vProjectileType)) +
pLab.squaredNorm());
auto particle = stack->AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
vProjectileType, E0, pLab, origin, 0_ns});
particle.SetNode(vNodePtr);
return std::make_tuple(
std::move(stack),
std::make_unique<decltype(corsika::stack::SecondaryView(particle))>(particle));
}
TEST_CASE("UrQMD") {
SECTION("conversion") {
REQUIRE_THROWS(process::UrQMD::ConvertFromUrQMD(106, 0));
REQUIRE(process::UrQMD::ConvertFromUrQMD(101, 0) == particles::Code::Pi0);
REQUIRE(process::UrQMD::ConvertToUrQMD(particles::Code::PiPlus) ==
std::make_pair<int, int>(101, 2));
}
feenableexcept(FE_INVALID);
corsika::random::RNGManager::GetInstance().RegisterRandomStream("UrQMD");
UrQMD urqmd;
SECTION("cross sections") {
auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Unknown);
auto const& cs = *csPtr;
particles::Code validProjectileCodes[] = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::Proton,
particles::Code::Neutron, particles::Code::KPlus, particles::Code::KMinus,
particles::Code::K0, particles::Code::K0Bar};
for (auto code : validProjectileCodes) {
auto [stack, view] = setupStack(code, 100_GeV, nodePtr, cs);
REQUIRE(stack->GetSize() == 1);
// simple check whether the cross-section is non-vanishing
REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Proton) /
1_mb >
0);
REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Nitrogen) /
1_mb >
0);
REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Oxygen) /
1_mb >
0);
REQUIRE(urqmd.GetCrossSection(view->GetProjectile(), particles::Code::Argon) /
1_mb >
0);
}
}
SECTION("nucleon projectile") {
auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
unsigned short constexpr A = 14, Z = 7;
auto [stackPtr, secViewPtr] = setupStack(A, Z, 400_GeV, nodePtr, *csPtr);
// must be assigned to variable, cannot be used as rvalue?!
auto projectile = secViewPtr->GetProjectile();
auto const projectileMomentum = projectile.GetMomentum();
[[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
REQUIRE(sumCharge(*secViewPtr) ==
Z + particles::GetChargeNumber(particles::Code::Oxygen));
auto const secMomSum =
sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());
REQUIRE((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
Approx(0).margin(1e-2));
}
SECTION("\"special\" projectile") {
auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
auto [stackPtr, secViewPtr] =
setupStack(particles::Code::PiPlus, 400_GeV, nodePtr, *csPtr);
// must be assigned to variable, cannot be used as rvalue?!
auto projectile = secViewPtr->GetProjectile();
auto const projectileMomentum = projectile.GetMomentum();
[[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
REQUIRE(sumCharge(*secViewPtr) ==
particles::GetChargeNumber(particles::Code::PiPlus) +
particles::GetChargeNumber(particles::Code::Oxygen));
auto const secMomSum =
sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());
REQUIRE((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
Approx(0).margin(1e-2));
}
}
README.md
View file @ 00b7c3f3
# CORSIKA 8 Framework for Particle Cascades in Astroparticle Physics
# CORSIKA 8 Framework for Particle Cascades in Astroparticle Physics
The purpose of CORSIKA is to simulate any particle cascades in
astroparticle physics or astrophysical context. A lot of emphasis is
The purpose of CORSIKA 8 is to simulate any particle cascades in
astroparticle physics or astrophysical context. A lot of emphasis has been
put on modularity, flexibility, completeness, validation and
correctness. To boost computational efficiency different techniques
correctness. To boost computational efficiency, different techniques
are provided, like thinning or cascade equations. The aim is that
CORSIKA remains the most comprehensive framework for simulating
CORSIKA 8 remains the most comprehensive framework for simulating
particle cascades with stochastic and continuous processes.
The software makes extensive use of static design patterns and
compiler optimization. Thus, the most fundamental configuration
decision of the user must be performed at compile time. At run time
only specific model parameters can still be changed.
decisions of the user must be performed at compile time. At run time,
model parameters can still be changed.
CORSIKA 8 is released under the GPLv3 license. This does not exclude
that specific CORSIKA 8 versions can be released for specific purposes
under different licensing. See [license
file](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/LICENSE)
CORSIKA 8 is by default released under the BSD 3-Clause License. See [license
file](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/blob/master/LICENSE)
which is part of every release and the source code.
If you use, or want to refer to, CORSIKA 8 please cite ["Towards a Next
Generation of CORSIKA: A Framework for the Simulation of Particle
Cascades in Astroparticle Physics", Comput.Softw.Big Sci. 3 (2019)
2](https://doi.org/10.1007/s41781-018-0013-0). We kindly ask (and
expect) any relevant improvement or addition to be offered or
Cascades in Astroparticle Physics", Comput. Softw. Big Sci. 3 (2019)
2](https://doi.org/10.1007/s41781-018-0013-0) as well as
["Simulating radio emission from particle cascades with CORSIKA 8", Astropart. Phys. 166 (2025)
103072](https://doi.org/10.1016/j.astropartphys.2024.103072).
We kindly ask (and require) any relevant improvement or addition to be offered or
contributed to the main CORSIKA 8 repository for the benefit of the
whole community.
When you plan to contribute to CORSIKA 8 check the guidelines outlined here:
CORSIKA 8 makes use of various third-party code, in particular interaction
models. Please check the [using and collaborating
agreement](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/blob/master/USING_COLLABORATING.md)
for further information on this topic.
If you plan to contribute to CORSIKA 8, please check the guidelines outlined here:
[coding
guidelines](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/CONTRIBUTING.md). Code
that fails the review by the CORSIKA author group must be improved
guidelines](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/blob/master/CONTRIBUTING.md). Code
that fails the review by the CORSIKA 8 author group must be improved
before it can be merged in the official code base. After your code has
been accepted and merged you become a contributor of the CORSIKA 8
project and you should include yourself in the
[AUTHORS](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/AUTHORS)
file.
IMPORTANT: Before you contribute, you need to read and agree to the
[collaboration
agreement](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/COLLABORATION_AGREEMENT.md). The
agreement can be discussed, and eventually improved.
been accepted and merged, you become a contributor of the CORSIKA 8
project (code author).
We also want to point you to the [MCnet
guidelines](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/MCNET_GUIDELINES),
which are very useful also for us.
IMPORTANT: Before you contribute, you need to read and agree to the conditions set out in the
[using and collaborating
agreement](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/blob/master/USING_COLLABORATING.md).
The agreement can be discussed, and eventually improved if necessary.
## Get in contact
* Connect to https://gitlab.ikp.kit.edu; register yourself and join the "Air Shower Physics" group
* Join our chat threads using Mattermost via this [invite link](https://mattermost.hzdr.de/signup_user_complete/?id=xtdd8jyt6trbiezt71gaz3z4ge&md=link&sbr=su). Click the `GitLab` button, then `Sign in with Helmholtz ID`. You will be able to make an account by either finding your institution, or using your e.g. ORCID, GitHub, or Google account.
* Connect to https://gitlab.iap.kit.edu, register yourself and join the "Air Shower Physics" group. Write to us on Mattermost (in the User Questions channel), or directly contact one of the [steering comittee members](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/-/wikis/Steering-Committee) in case there are problems with that.
* Connect to corsika-devel@lists.kit.edu (self-register at
https://www.lists.kit.edu/sympa/subscribe/corsika-devel) to get in
touch with the project
touch with the project.
## Installation
CORSIKA 8 is tested regularly at least on gcc7.3.0 and clang-6.0.0.
Additional software prerequisites: eigen3, boost, cmake, g++, git.
On a bare Ubuntu 18.04, just add:
CORSIKA 8 is tested regularly at least on `gcc11.0.0` and `clang-14.0.0`.
### Prerequisites
You will also need:
- Python 3 (supported versions are Python >= 3.6), with pip
- cmake > 3.4
- git
- g++, gfortran, binutils, make
- optional: FLUKA (see below)
On a bare Ubuntu machine, just add:
``` shell
sudo apt-get install python3 python3-pip cmake g++ gfortran git doxygen graphviz
```
sudo apt-get install libeigen3-dev libboost-dev cmake g++ git
### Creating a virtual environment and Conan
It is recommended that you install CORSIKA 8 and its dependencies within a python3 virtual environment.
To do so, you can run the following.
``` shell
# Create the environment using your native python3 binary
python3 -m venv /path/to/new/virtual/environment/corsika-8
# Load the environment (should be run each time you open a new terminal)
source /path/to/new/virtual/environment/corsika-8/bin/activate
```
Follow these steps to download and install CORSIKA 8 milestone2
You will need to load the environment each time that you open a new terminal.
CORSIKA 8 uses the [conan](https://conan.io/) package manager to
manage our dependencies. Currently, version 2.50.0 or higher is required.
**Note**: if you are NOT using a virtual environment, you may want to use the `pip install --user` flag.
``` shell
pip install conan particle==0.25.1 numpy
```
git clone git@gitlab.ikp.kit.edu:AirShowerPhysics/corsika.git
cd corsika
git checkout milestone2
mkdir ../corsika-build
cd ../corsika-build
cmake ../corsika -DCMAKE_INSTALL_PREFIX=../corsika-install
make -j8
### Enabling FLUKA support
For legal reasons we do not distribute/bundle FLUKA together with CORSIKA 8.
As FLUKA is the standard low-energy hadronic interaction model for CORSIKA 8, you have to download
it separately from (http://www.fluka.org/), which requires registering there as FLUKA user.
The following should be done *before* compiling CORSIKA 8:
1. Note your system's version of gfortran (`gfortran --version`) and glibc (`ldd --version`)
2. Download the FLUKA __binary__, ensuring that it matches the versions you found above
3. Download the FLUKA __data file__ (will be named something similar to __fluka20xy.z-data.tar.gz__).
4. Un-tar the files that you downloaded using `tar -xf <filename>`
5. Set environmental variables `export FLUFOR=gfortran` and `export FLUPRO=<path to where you unzipped the files>`. Note that the `FLUPRO` directory should contain __libflukahp.a__. Both of these variables will have to be set every time you open a new terminal.
6. Go to the `FLUPRO` directory and run `make`. This will compile an exe, __flukahp__, in your current directory.
7. Follow the normal steps to compile CORSIKA 8 (see below).
When you later install CORSIKA 8, you should see a message during the __cmake__ step indicating the FLUKA was correctly found.
``` shell
libflukahp.a found in directory <some location here> via FLUPRO environment variable
FLUKA support is enabled.
```
### Compiling CORSIKA 8
Once Conan is installed and FLUKA provided, follow these steps to download and install CORSIKA 8:
``` shell
cd ./top/directory/for/corsika/installation
git clone --recursive git@gitlab.iap.kit.edu:AirShowerPhysics/corsika.git
# Or for https: git clone --recursive https://gitlab.iap.kit.edu/AirShowerPhysics/corsika.git
mkdir corsika-build
cd corsika-build
../corsika/conan-install.sh --source-directory ../corsika --release-with-debug
# conan-install.sh takes required options from command line to install dependencies for 'Debug', 'Release' and 'RelWithDebInfo' builds.
../corsika/corsika-cmake.sh -c "-DCMAKE_BUILD_TYPE="RelWithDebInfo" -DWITH_FLUKA=ON -DCMAKE_INSTALL_PREFIX=../corsika-install"
make -j4 #The number should match the number of available cores on your machine
make install
make test
```
and if you want to see how the first simple hadron cascade develops, see `Documentation/Examples/cascade_example.cc` for a starting point.
Run the cascade_example with:
## Alternate installation using docker containers
There are docker containers prepared that bring all the environment and packages you need to run CORSIKA. See [docker hub](https://hub.docker.com/repository/docker/corsika/devel) for a complete overview.
### Prerequisites
You only need docker, e.g. on Ubuntu: `sudo apt-get install docker` and of course root access.
### Compiling
Follow these steps to download and install CORSIKA 8, master development version
```shell
cd ./top/directory/for/corsika/installation
git clone --recursive git@gitlab.iap.kit.edu:AirShowerPhysics/corsika.git
sudo docker run -v $PWD:/corsika -it corsika/devel:clang-8 /bin/bash
mkdir corsika-build
cd corsika-build
../corsika/conan-install.sh --source-directory ../corsika --release-with-debug
# conan-install.sh takes required options from command line to install dependencies for 'Debug', 'Release' and 'RelWithDebInfo' builds.
../corsika/corsika-cmake.sh -c "-DCMAKE_BUILD_TYPE="RelWithDebInfo" -DWITH_FLUKA=ON -DCMAKE_INSTALL_PREFIX=../corsika-install"
make -j4 #The number should match the number of available cores on your machine
make install
```
cd ../corsika-install
share/examples/cascade_example
## Running Unit Tests
To run the unit tests, do the following.
```shell
cd ./corsika-build
ctest -j4 #The number should match the number of available cores on your machine
```
## Running applications and examples
### Standard applications
Applications for standard use-cases are located in the `applications` directory.
These are example scripts that can be used directly or slightly modified for your use case.
See [applications/README.md] for more.
The applications are compiled automatically after running `make` and will appear your `corsika-build/bin` directory.
After running `make install` the binaries will also be copied into your `corsika-install/bin` directory as well.
For example, from inside your `corsika-install/bin` directory, run
```shell
c8_air_shower --pdg 2212 -E 1e5 -f my_shower
```
This will run a vertical 100 TeV proton shower and will create and put the output into `./my_shower`.
Visualize output (needs gnuplot installed):
### Building the examples
Unlike the applications, the examples must be compiled as a second step.
From your top corsika directory, (the one that includes `corsika-build` and `corsika-install`) run
```shell
export CONAN_DEPENDENCIES=$PWD/corsika-install/lib/cmake/dependencies
cmake -DCMAKE_TOOLCHAIN_FILE=${CONAN_DEPENDENCIES}/conan_toolchain.cmake -DCMAKE_PREFIX_PATH=${CONAN_DEPENDENCIES} -DCMAKE_POLICY_DEFAULT_CMP0091=NEW -DCMAKE_BUILD_TYPE=RelWithDebInfo -Dcorsika_DIR=$PWD/corsika-build -DWITH_FLUKA=ON -S $PWD/corsika/examples -B $PWD/corsika-build-examples
cd corsika-build-examples
make -j4 #The number should match the number of available cores on your machine
```
bash share/tools/plot_tracks.sh tracks.dat
firefox tracks.dat.gif
You can run the examples by using the binaries in `corsika-build-examples/bin/`.
For example:
```shell
corsika-build-examples/bin/known_particles
```
This will print out all of the particles that are known by CORSIKA.
### Generating doxygen documentation
To generate the documentation, you need doxygen and graphviz. On Ubuntu 18.04, do:
```
To generate the documentation, you need doxygen and graphviz. If you work with
the docker corsika/devel containers this is already included.
Otherwise, e.g. on Ubuntu machines, do:
```shell
sudo apt-get install doxygen graphviz
```
Switch to the corsika build directory and do
```
make doxygen
Switch to the `corsika-build` directory and do
```shell
make docs
make install
```
browse with firefox:
open with firefox:
```shell
firefox ../corsika-install/share/corsika/doc/html/index.html
```
firefox ../corsika-install/share/doc/html/index.html
```
SCIENTIFIC_AUTHORS deleted 100644 → 0
View file @ 6f2ab230
See https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/wikis/ICRC2019
Setup/CMakeLists.txt deleted 100644 → 0
View file @ 6f2ab230
set (
SETUP_HEADERS
SetupStack.h
SetupLogger.h
SetupEnvironment.h
SetupTrajectory.h
)
set (
SETUP_NAMESPACE
corsika/setup
)
add_library (CORSIKAsetup INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAsetup ${SETUP_NAMESPACE} ${SETUP_HEADERS})
target_link_libraries (
CORSIKAsetup
INTERFACE
CORSIKAgeometry
SuperStupidStack
NuclearStackExtension
)
target_include_directories (
CORSIKAsetup
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
FILES ${SETUP_HEADERS}
DESTINATION include/${SETUP_NAMESPACE}
)
Setup/SetupEnvironment.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (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_setup_environment_h_
#define _include_corsika_setup_environment_h_
#include <corsika/environment/Environment.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NameModel.h>
namespace corsika::setup {
using IEnvironmentModel = environment::IMediumModel;
using SetupEnvironment = environment::Environment<IEnvironmentModel>;
} // namespace corsika::setup
#endif
Setup/SetupLogger.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (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_setup_logger_h_
#define _include_corsika_setup_logger_h_
namespace corsika {}
#endif
Setup/SetupStack.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (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 _corsika_setup_setupstack_h_
#define _corsika_setup_setupstack_h_
// the basic particle data stack:
#include <corsika/stack/super_stupid/SuperStupidStack.h>
// extension with nuclear data for Code::Nucleus
#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
// extension with geometry information for tracking
#include <corsika/environment/Environment.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/stack/CombinedStack.h>
#include <tuple>
#include <utility>
#include <vector>
// definition of stack-data object to store geometry information
template <typename TEnvType>
/**
* @class GeometryData
*
* definition of stack-data object to store geometry information
*/
class GeometryData {
public:
using BaseNodeType = typename TEnvType::BaseNodeType;
// these functions are needed for the Stack interface
void Init() {}
void Clear() { fNode.clear(); }
unsigned int GetSize() const { return fNode.size(); }
unsigned int GetCapacity() const { return fNode.size(); }
void Copy(const int i1, const int i2) { fNode[i2] = fNode[i1]; }
void Swap(const int i1, const int i2) { std::swap(fNode[i1], fNode[i2]); }
// custom data access function
void SetNode(const int i, BaseNodeType const* v) { fNode[i] = v; }
auto const* GetNode(const int i) const { return fNode[i]; }
// these functions are also needed by the Stack interface
void IncrementSize() { fNode.push_back(nullptr); }
void DecrementSize() {
if (fNode.size() > 0) { fNode.pop_back(); }
}
// custom private data section
private:
std::vector<const BaseNodeType*> fNode;
};
/**
* @class GeometryDataInterface
*
* corresponding defintion of a stack-readout object, the iteractor
* dereference operator will deliver access to these function
// defintion of a stack-readout object, the iteractor dereference
// operator will deliver access to these function
*/
template <typename T, typename TEnvType>
class GeometryDataInterface : public T {
public:
using T::GetIndex;
using T::GetStackData;
using T::SetParticleData;
using BaseNodeType = typename TEnvType::BaseNodeType;
// default version for particle-creation from input data
void SetParticleData(const std::tuple<BaseNodeType const*> v) {
SetNode(std::get<0>(v));
}
void SetParticleData(GeometryDataInterface& parent,
const std::tuple<BaseNodeType const*>) {
SetNode(parent.GetNode()); // copy Node from parent particle!
}
void SetParticleData() { SetNode(nullptr); }
void SetParticleData(GeometryDataInterface& parent) {
SetNode(parent.GetNode()); // copy Node from parent particle!
}
void SetNode(BaseNodeType const* v) { GetStackData().SetNode(GetIndex(), v); }
auto const* GetNode() const { return GetStackData().GetNode(GetIndex()); }
};
namespace corsika::setup {
namespace detail {
//
// this is an auxiliary help typedef, which I don't know how to put
// into NuclearStackExtension.h where it belongs...
template <typename StackIter>
using ExtendedParticleInterfaceType =
corsika::stack::nuclear_extension::NuclearParticleInterface<
corsika::stack::super_stupid::SuperStupidStack::PIType, StackIter>;
//
// the particle data stack with extra nuclear information:
using ParticleDataStack = corsika::stack::nuclear_extension::NuclearStackExtension<
corsika::stack::super_stupid::SuperStupidStack, ExtendedParticleInterfaceType>;
template <typename T>
using SetupGeometryDataInterface = GeometryDataInterface<T, setup::SetupEnvironment>;
// combine particle data stack with geometry information for tracking
template <typename StackIter>
using StackWithGeometryInterface =
corsika::stack::CombinedParticleInterface<ParticleDataStack::PIType,
SetupGeometryDataInterface, StackIter>;
using StackWithGeometry =
corsika::stack::CombinedStack<typename ParticleDataStack::StackImpl,
GeometryData<setup::SetupEnvironment>,
StackWithGeometryInterface>;
} // namespace detail
// this is the REAL stack we use:
using Stack = detail::StackWithGeometry;
/*
See Issue 161
unfortunately clang does not support this in the same way (yet) as
gcc, so we have to distinguish here. If clang cataches up, we
could remove the clang branch here and also in
corsika::Cascade. The gcc code is much more generic and
universal. If we could do the gcc version, we won't had to define
StackView globally, we could do it with MakeView whereever it is
actually needed. Keep an eye on this!
*/
#if defined(__clang__)
using StackView =
corsika::stack::SecondaryView<typename corsika::setup::Stack::StackImpl,
corsika::setup::detail::StackWithGeometryInterface>;
#elif defined(__GNUC__) || defined(__GNUG__)
using StackView = corsika::stack::MakeView<corsika::setup::Stack>::type;
#endif
} // namespace corsika::setup
#endif
Setup/SetupTrajectory.h deleted 100644 → 0
View file @ 6f2ab230
/*
* (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 _corsika_setup_setuptrajectory_h_
#define _corsika_setup_setuptrajectory_h_
#include <corsika/geometry/Helix.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/units/PhysicalUnits.h>
// #include <variant>
namespace corsika::setup {
/// definition of Trajectory base class, to be used in tracking and cascades
typedef corsika::geometry::Trajectory<corsika::geometry::Line> Trajectory;
/*
typedef std::variant<std::monostate, corsika::geometry::Trajectory<Line>,
corsika::geometry::Trajectory<Helix>>
Trajectory;
/// helper visitor to modify Particle by moving along Trajectory
template <typename Particle>
class ParticleUpdate {
Particle& fP;
public:
ParticleUpdate(Particle& p)
: fP(p) {}
void operator()(std::monostate const&) {}
template <typename T>
void operator()(T const& trajectory) {
fP.SetPosition(trajectory.GetPosition(1));
}
};
/// helper visitor to modify Particle by moving along Trajectory
class GetDuration {
public:
corsika::units::si::TimeType operator()(std::monostate const&) {
return 0 * corsika::units::si::second;
}
template <typename T>
corsika::units::si::TimeType operator()(T const& trajectory) {
return trajectory.GetDuration();
}
};
*/
} // namespace corsika::setup
#endif

data privacy - impressum