IAP GITLAB

Skip to content
Snippets Groups Projects

Compare revisions

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

Source

Select target project
No results found

Target

Select target project
  • AirShowerPhysics/corsika
  • rulrich/corsika
  • AAAlvesJr/corsika
  • Andre/corsika
  • arrabito/corsika
  • Nikos/corsika
  • olheiser73/corsika
  • AirShowerPhysics/papers/corsika
  • pranav/corsika
9 results
Show changes
Showing
with 0 additions and 1547 deletions
add_subdirectory (Doxygen)
add_subdirectory (Examples)
add_executable (helix_example helix_example.cc)
target_compile_options(helix_example PRIVATE -g) # do not skip asserts
target_link_libraries(helix_example CORSIKAgeometry CORSIKAunits)
install (TARGETS helix_example DESTINATION share/examples)
CORSIKA_ADD_TEST (helix_example)
add_executable (geometry_example geometry_example.cc)
target_compile_options(geometry_example PRIVATE -g) # do not skip asserts
target_link_libraries (geometry_example CORSIKAgeometry CORSIKAunits)
install (TARGETS geometry_example DESTINATION share/examples)
CORSIKA_ADD_TEST (geometry_example)
add_executable (logger_example logger_example.cc)
target_compile_options(logger_example PRIVATE -g) # do not skip asserts
target_link_libraries (logger_example CORSIKAunits CORSIKAlogging)
install (TARGETS logger_example DESTINATION share/examples)
CORSIKA_ADD_TEST (logger_example)
add_executable (stack_example stack_example.cc)
target_compile_options(stack_example PRIVATE -g) # do not skip asserts
target_link_libraries (stack_example SuperStupidStack CORSIKAunits CORSIKAlogging)
CORSIKA_ADD_TEST (stack_example)
add_executable (cascade_example cascade_example.cc)
target_compile_options(cascade_example PRIVATE -g) # do not skip asserts
target_link_libraries (cascade_example SuperStupidStack CORSIKAunits CORSIKAlogging
CORSIKArandom
ProcessSibyll
CORSIKAcascade
ProcessStackInspector
ProcessTrackWriter
ProcessHadronicElasticModel
CORSIKAprocesses
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
install (TARGETS cascade_example DESTINATION share/examples)
CORSIKA_ADD_TEST (cascade_example)
add_executable (staticsequence_example staticsequence_example.cc)
target_compile_options(staticsequence_example PRIVATE -g) # do not skip asserts
target_link_libraries (staticsequence_example
CORSIKAprocesssequence
CORSIKAunits
CORSIKAgeometry
CORSIKAlogging)
install (TARGETS staticsequence_example DESTINATION share/examples)
CORSIKA_ADD_TEST (staticsequence_example)
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
#include <corsika/utl/CorsikaFenv.h>
#include <boost/type_index.hpp>
using boost::typeindex::type_id_with_cvr;
#include <iostream>
#include <limits>
#include <typeinfo>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
class ProcessCut : public corsika::process::ContinuousProcess<ProcessCut> {
HEPEnergyType fECut;
HEPEnergyType fEnergy = 0_GeV;
HEPEnergyType fEmEnergy = 0_GeV;
int fEmCount = 0;
HEPEnergyType fInvEnergy = 0_GeV;
int fInvCount = 0;
public:
ProcessCut(const HEPEnergyType cut)
: fECut(cut) {}
template <typename Particle>
bool isBelowEnergyCut(Particle& p) const {
// nuclei
if (p.GetPID() == corsika::particles::Code::Nucleus) {
auto const ElabNuc = p.GetEnergy() / p.GetNuclearA();
auto const EcmNN = sqrt(2. * ElabNuc * 0.93827_GeV);
if (ElabNuc < fECut || EcmNN < 10_GeV)
return true;
else
return false;
} else {
// TODO: center-of-mass energy hard coded
const HEPEnergyType Ecm = sqrt(2. * p.GetEnergy() * 0.93827_GeV);
if (p.GetEnergy() < fECut || Ecm < 10_GeV)
return true;
else
return false;
}
}
bool isEmParticle(Code pCode) const {
bool is_em = false;
// FOR NOW: switch
switch (pCode) {
case Code::Electron:
is_em = true;
break;
case Code::Positron:
is_em = true;
break;
case Code::Gamma:
is_em = true;
break;
default:
break;
}
return is_em;
}
void defineEmParticles() const {
// create bool array identifying em particles
}
bool isInvisible(Code pCode) const {
bool is_inv = false;
// FOR NOW: switch
switch (pCode) {
case Code::NuE:
is_inv = true;
break;
case Code::NuEBar:
is_inv = true;
break;
case Code::NuMu:
is_inv = true;
break;
case Code::NuMuBar:
is_inv = true;
break;
case Code::MuPlus:
is_inv = true;
break;
case Code::MuMinus:
is_inv = true;
break;
case Code::Neutron:
is_inv = true;
break;
case Code::AntiNeutron:
is_inv = true;
break;
default:
break;
}
return is_inv;
}
template <typename Particle>
LengthType MaxStepLength(Particle& p, setup::Trajectory&) const {
cout << "ProcessCut: MinStep: pid: " << p.GetPID() << endl;
cout << "ProcessCut: MinStep: energy (GeV): " << p.GetEnergy() / 1_GeV << endl;
const Code pid = p.GetPID();
if (isEmParticle(pid) || isInvisible(pid) || isBelowEnergyCut(p)) {
cout << "ProcessCut: MinStep: next cut: " << 0. << endl;
return 0_m;
} else {
LengthType next_step = 1_m * std::numeric_limits<double>::infinity();
cout << "ProcessCut: MinStep: next cut: " << next_step << endl;
return next_step;
}
}
template <typename Particle, typename Stack>
EProcessReturn DoContinuous(Particle& p, setup::Trajectory&, Stack&) {
const Code pid = p.GetPID();
HEPEnergyType energy = p.GetEnergy();
cout << "ProcessCut: DoContinuous: " << pid << " E= " << energy
<< ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV" << endl;
EProcessReturn ret = EProcessReturn::eOk;
if (isEmParticle(pid)) {
cout << "removing em. particle..." << endl;
fEmEnergy += energy;
fEmCount += 1;
// p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} else if (isInvisible(pid)) {
cout << "removing inv. particle..." << endl;
fInvEnergy += energy;
fInvCount += 1;
// p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
} else if (isBelowEnergyCut(p)) {
cout << "removing low en. particle..." << endl;
fEnergy += energy;
// p.Delete();
ret = EProcessReturn::eParticleAbsorbed;
}
return ret;
}
void Init() {
fEmEnergy = 0. * 1_GeV;
fEmCount = 0;
fInvEnergy = 0. * 1_GeV;
fInvCount = 0;
fEnergy = 0. * 1_GeV;
// defineEmParticles();
}
void ShowResults() {
cout << " ******************************" << endl
<< " ParticleCut: " << endl
<< " energy in em. component (GeV): " << fEmEnergy / 1_GeV << endl
<< " no. of em. particles injected: " << fEmCount << endl
<< " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << endl
<< " no. of inv. particles injected: " << fInvCount << endl
<< " energy below particle cut (GeV): " << fEnergy / 1_GeV << endl
<< " ******************************" << endl;
}
HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
HEPEnergyType GetCutEnergy() const { return fEnergy; }
HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
};
//
// The example main program for a particle cascade
//
int main() {
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
corsika::random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
corsika::environment::Environment env;
auto& universe = *(env.GetUniverse());
auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
// fraction of oxygen
const float fox = 0.20946;
using MyHomogeneousModel =
corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
corsika::environment::NuclearComposition(
std::vector<corsika::particles::Code>{corsika::particles::Code::Nitrogen,
corsika::particles::Code::Oxygen},
std::vector<float>{(float)1. - fox, fox}));
universe.AddChild(std::move(theMedium));
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
// setup processes, decays and interactions
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
corsika::process::sibyll::Interaction sibyll(env);
corsika::process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
corsika::process::sibyll::Decay decay;
ProcessCut cut(20_GeV);
// corsika::random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// corsika::process::HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
corsika::process::TrackWriter::TrackWriter trackWriter("tracks.dat");
// assemble all processes into an ordered process list
// auto sequence = p0 << sibyll << decay << hadronicElastic << cut << trackWriter;
auto sequence = p0 << sibyll << sibyllNuc << decay << cut << trackWriter;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
const Code beamCode = Code::Nucleus;
const int nuclA = 56;
const int nuclZ = int(nuclA / 2.15 + 0.7);
const HEPMassType mass = corsika::particles::Proton::GetMass() * nuclZ +
(nuclA - nuclZ) * corsika::particles::Neutron::GetMass();
const HEPEnergyType E0 =
nuclA *
100_GeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
double theta = 0.;
double phi = 0.;
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt(Elab * Elab - m * m);
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = stack::super_stupid::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
stack.AddParticle(beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ);
}
// define air shower object, run simulation
corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Init();
EAS.Run();
cout << "Result: E0=" << E0 / 1_GeV << endl;
cut.ShowResults();
const HEPEnergyType Efinal =
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
cout << "total energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/logging/Logger.h>
#include <boost/format.hpp>
#include <fstream>
#include <iostream>
#include <string>
using namespace std;
using namespace corsika::logging;
int main() {
{
cout << "writing to \"another.log\"" << endl;
ofstream logfile("another.log");
sink::SinkStream unbuffered_sink(logfile);
sink::BufferedSinkStream sink(logfile, sink::StdBuffer(10000));
Logger<MessageOn, sink::BufferedSinkStream> info("\033[32m", "info", sink);
Logger<MessageOn, sink::BufferedSinkStream> err("\033[31m", "error", sink);
// logger<ostream,messageconst,StdBuffer> info(std::cout, StdBuffer(10000));
/*
Logging& logs = Logging::GetInstance();
logs.AddLogger<>("info", info);
auto& log_1 = logs.GetLogger("info"); // no so useful, since type of log_1 is
std::any
*/
for (int i = 0; i < 10000; ++i) {
LOG(info, "irgendwas", " ", string("and more"), " ",
boost::format("error: %i message: %s. done."), i, "stupido");
LOG(err, "Fehler");
}
}
{
sink::NoSink off;
Logger<MessageOff> info("", "", off);
for (int i = 0; i < 10000; ++i) {
LOG(info, "irgendwas", string("and more"),
boost::format("error: %i message: %s. done."), i, "stupido", "a-number:", 8.99,
"ENDE");
}
}
return 0;
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/particles/ParticleProperties.h>
#include <corsika/stack/super_stupid/SuperStupidStack.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <cassert>
#include <iomanip>
#include <iostream>
using namespace corsika;
using namespace corsika::units::si;
using namespace corsika::stack;
using namespace corsika::geometry;
using namespace std;
void fill(corsika::stack::super_stupid::SuperStupidStack& s) {
const geometry::CoordinateSystem& rootCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
for (int i = 0; i < 11; ++i) {
s.AddParticle(corsika::particles::Code::Electron, 1.5_GeV * i,
stack::super_stupid::MomentumVector(rootCS, {0_GeV, 0_GeV, 1_GeV}),
geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns);
}
}
void read(corsika::stack::super_stupid::SuperStupidStack& s) {
assert(s.GetSize() == 11); // stack has 11 particles
HEPEnergyType total_energy;
int i = 0;
for (auto& p : s) {
total_energy += p.GetEnergy();
// particles are electrons with 1.5 GeV energy times i
assert(p.GetPID() == corsika::particles::Code::Electron);
assert(p.GetEnergy() == 1.5_GeV * (i++));
}
}
int main() {
corsika::stack::super_stupid::SuperStupidStack s;
fill(s);
read(s);
return 0;
}
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <array>
#include <iomanip>
#include <iostream>
#include <corsika/process/ProcessSequence.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
using namespace corsika;
using namespace corsika::units::si;
using namespace corsika::process;
using namespace std;
const int nData = 10;
class Process1 : public BaseProcess<Process1> {
public:
Process1() {}
template <typename D, typename T, typename S>
EProcessReturn DoContinuous(D& d, T&, S&) const {
for (int i = 0; i < nData; ++i) d.p[i] += 1;
return EProcessReturn::eOk;
}
};
class Process2 : public BaseProcess<Process2> {
public:
Process2() {}
template <typename D, typename T, typename S>
inline EProcessReturn DoContinuous(D& d, T&, S&) const {
for (int i = 0; i < nData; ++i) d.p[i] -= 0.1 * i;
return EProcessReturn::eOk;
}
};
class Process3 : public BaseProcess<Process3> {
public:
Process3() {}
template <typename D, typename T, typename S>
inline EProcessReturn DoContinuous(D&, T&, S&) const {
return EProcessReturn::eOk;
}
};
class Process4 : public BaseProcess<Process4> {
public:
Process4(const double v)
: fV(v) {}
template <typename D, typename T, typename S>
inline EProcessReturn DoContinuous(D& d, T&, S&) const {
for (int i = 0; i < nData; ++i) d.p[i] *= fV;
return EProcessReturn::eOk;
}
private:
double fV;
};
struct DummyData {
double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
};
struct DummyStack {
void clear() {}
};
struct DummyTrajectory {};
void modular() {
Process1 m1;
Process2 m2;
Process3 m3;
Process4 m4(0.9);
auto sequence = m1 << m2 << m3 << m4;
DummyData p;
DummyStack s;
DummyTrajectory t;
const int n = 1000;
for (int i = 0; i < n; ++i) { sequence.DoContinuous(p, t, s); }
for (int i = 0; i < nData; ++i) {
// cout << p.p[i] << endl;
// assert(p.p[i] == n-i*100);
}
cout << " done (nothing...) " << endl;
}
int main() {
modular();
return 0;
}
set (
ENVIRONMENT_HEADERS
VolumeTreeNode.h
IMediumModel.h
NuclearComposition.h
HomogeneousMedium.h
Environment.h
)
set (
ENVIRONMENT_NAMESPACE
corsika/environment
)
add_library (CORSIKAenvironment INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAenvironment ${ENVIRONMENT_NAMESPACE} ${ENVIRONMENT_HEADERS})
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
CORSIKAenvironment
INTERFACE
CORSIKAgeometry
CORSIKAparticles
CORSIKAunits
CORSIKArandom
)
target_include_directories (
CORSIKAenvironment
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include>
)
install (
TARGETS CORSIKAenvironment
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
PUBLIC_HEADER DESTINATION include/${ENVIRONMENT_NAMESPACE}
)
# --------------------
# code unit testing
add_executable (testEnvironment testEnvironment.cc)
target_link_libraries (
testEnvironment
CORSIKAsetup
CORSIKAenvironment
CORSIKAthirdparty # for catch2
)
CORSIKA_ADD_TEST(testGeometry)
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_Environment_h
#define _include_Environment_h
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/setup/SetupEnvironment.h>
#include <limits>
namespace corsika::environment {
using BaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
struct Universe : public corsika::geometry::Sphere {
Universe(corsika::geometry::CoordinateSystem const& pCS)
: corsika::geometry::Sphere(
corsika::geometry::Point{pCS, 0 * corsika::units::si::meter,
0 * corsika::units::si::meter,
0 * corsika::units::si::meter},
corsika::units::si::meter * std::numeric_limits<double>::infinity()) {}
bool Contains(corsika::geometry::Point const&) const override { return true; }
};
// template <typename IEnvironmentModel>
class Environment {
public:
Environment()
: fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem()}
, fUniverse(std::make_unique<BaseNodeType>(
std::make_unique<Universe>(fCoordinateSystem))) {}
using IEnvironmentModel = corsika::setup::IEnvironmentModel;
auto& GetUniverse() { return fUniverse; }
auto const& GetUniverse() const { return fUniverse; }
auto const& GetCoordinateSystem() const { return fCoordinateSystem; }
// factory method for creation of VolumeTreeNodes
template <class TVolumeType, typename... TVolumeArgs>
static auto CreateNode(TVolumeArgs&&... args) {
static_assert(std::is_base_of_v<corsika::geometry::Volume, TVolumeType>,
"unusable type provided, needs to be derived from "
"\"corsika::geometry::Volume\"");
return std::make_unique<BaseNodeType>(
std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
}
private:
corsika::geometry::CoordinateSystem const& fCoordinateSystem;
BaseNodeType::VTNUPtr fUniverse;
};
} // namespace corsika::environment
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
/*
* (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_HomogeneousMedium_h_
#define _include_HomogeneousMedium_h_
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
#include <cassert>
/**
* a homogeneous medium
*/
namespace corsika::environment {
template <class T>
class HomogeneousMedium : public T {
corsika::units::si::MassDensityType const fDensity;
NuclearComposition const fNuclComp;
public:
HomogeneousMedium(corsika::units::si::MassDensityType pDensity,
NuclearComposition pNuclComp)
: fDensity(pDensity)
, fNuclComp(pNuclComp){};
corsika::units::si::MassDensityType GetMassDensity(
corsika::geometry::Point const&) const override {
return fDensity;
}
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::LengthType pTo) const override {
using namespace corsika::units::si;
return pTo * fDensity;
}
corsika::units::si::LengthType ArclengthFromGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::GrammageType pGrammage) const override {
return pGrammage / fDensity;
}
};
} // namespace corsika::environment
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_IMediumModel_h
#define _include_IMediumModel_h
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <random>
namespace corsika::environment {
class IMediumModel {
public:
virtual ~IMediumModel() = default;
virtual corsika::units::si::MassDensityType GetMassDensity(
corsika::geometry::Point const&) const = 0;
// todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
// approach for now, only lines are supported
virtual corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::LengthType) const = 0;
virtual corsika::units::si::LengthType ArclengthFromGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const&,
corsika::units::si::GrammageType) const = 0;
virtual NuclearComposition const& GetNuclearComposition() const = 0;
template <class TRNG>
corsika::particles::Code SampleTarget(
std::vector<corsika::units::si::CrossSectionType> const& sigma,
TRNG& randomStream) const {
using namespace corsika::units::si;
auto const& nuclComp = GetNuclearComposition();
auto const& fractions = nuclComp.GetFractions();
assert(sigma.size() == fractions.size());
std::vector<float> weights(fractions.size());
for (size_t i = 0; i < fractions.size(); ++i) {
weights[i] = fractions[i] * sigma[i].magnitude();
}
std::discrete_distribution channelDist(weights.cbegin(), weights.cend());
const int iChannel = channelDist(randomStream);
return nuclComp.GetComponents()[iChannel];
}
};
} // namespace corsika::environment
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_NuclearComposition_h
#define _include_NuclearComposition_h
#include <corsika/particles/ParticleProperties.h>
#include <cassert>
#include <numeric>
#include <stdexcept>
#include <vector>
namespace corsika::environment {
class NuclearComposition {
std::vector<float> const fNumberFractions; //!< relative fractions of number density
std::vector<corsika::particles::Code> const
fComponents; //!< particle codes of consitutents
double const fAvgMassNumber;
public:
NuclearComposition(std::vector<corsika::particles::Code> pComponents,
std::vector<float> pFractions)
: fNumberFractions(pFractions)
, fComponents(pComponents)
, fAvgMassNumber(std::inner_product(
pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0.,
[](double x, double y) { return x + y; },
[](auto const& compID, auto const& fraction) {
return corsika::particles::GetNucleusA(compID) * fraction;
})) {
assert(pComponents.size() == pFractions.size());
auto const sumFractions =
std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
throw std::runtime_error("element fractions do not add up to 1");
}
}
auto size() const { return fNumberFractions.size(); }
auto const& GetFractions() const { return fNumberFractions; }
auto const& GetComponents() const { return fComponents; }
auto const GetAverageMassNumber() const { return fAvgMassNumber; }
};
} // namespace corsika::environment
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#ifndef _include_VolumeTreeNode_H
#define _include_VolumeTreeNode_H
#include <corsika/geometry/Volume.h>
#include <memory>
#include <vector>
namespace corsika::environment {
class Empty {}; //<! intended for usage as default template argument
template <typename IModelProperties = Empty>
class VolumeTreeNode {
public:
using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>;
using IMPSharedPtr = std::shared_ptr<IModelProperties>;
using VolUPtr = std::unique_ptr<corsika::geometry::Volume>;
VolumeTreeNode(VolUPtr pVolume = nullptr)
: fGeoVolume(std::move(pVolume)) {}
bool Contains(corsika::geometry::Point const& p) const {
return fGeoVolume->Contains(p);
}
VolumeTreeNode<IModelProperties> const* Excludes(
corsika::geometry::Point const& p) const {
auto exclContainsIter =
std::find_if(fExcludedNodes.cbegin(), fExcludedNodes.cend(),
[&](auto const& s) { return bool(s->Contains(p)); });
return exclContainsIter != fExcludedNodes.cend() ? *exclContainsIter : nullptr;
}
/** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given
* \class Point \arg p, or nullptr iff \arg p is not contained in this volume.
*/
VolumeTreeNode<IModelProperties> const* GetContainingNode(
corsika::geometry::Point const& p) const {
if (!Contains(p)) { return nullptr; }
if (auto const childContainsIter =
std::find_if(fChildNodes.cbegin(), fChildNodes.cend(),
[&](auto const& s) { return bool(s->Contains(p)); });
childContainsIter == fChildNodes.cend()) // not contained in any of the children
{
if (auto const exclContainsIter = Excludes(p)) // contained in any excluded nodes
{
return exclContainsIter->GetContainingNode(p);
} else {
return this;
}
} else {
return (*childContainsIter)->GetContainingNode(p);
}
}
void AddChild(VTNUPtr pChild) {
pChild->fParentNode = this;
fChildNodes.push_back(std::move(pChild));
// It is a bad idea to return an iterator to the inserted element
// because it might get invalidated when the vector needs to grow
// later and the caller won't notice.
}
void ExcludeOverlapWith(VTNUPtr const& pNode) {
fExcludedNodes.push_back(pNode.get());
}
auto* GetParent() const { return fParentNode; };
auto const& GetChildNodes() const { return fChildNodes; }
auto const& GetExcludedNodes() const { return fExcludedNodes; }
auto const& GetVolume() const { return *fGeoVolume; }
auto const& GetModelProperties() const { return *fModelProperties; }
template <typename TModelProperties, typename... Args>
auto SetModelProperties(Args&&... args) {
static_assert(std::is_base_of_v<IModelProperties, TModelProperties>,
"unusable type provided");
fModelProperties = std::make_shared<TModelProperties>(std::forward<Args>(args)...);
return fModelProperties;
}
void SetModelProperties(IMPSharedPtr ptr) { fModelProperties = ptr; }
template <class MediumType, typename... Args>
static auto CreateMedium(Args&&... args) {
static_assert(std::is_base_of_v<IMediumModel, MediumType>,
"unusable type provided, needs to be derived from \"IMediumModel\"");
return std::make_shared<MediumType>(std::forward<Args>(args)...);
}
private:
std::vector<VTNUPtr> fChildNodes;
std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
VolumeTreeNode<IModelProperties> const* fParentNode = nullptr;
VolUPtr fGeoVolume;
IMPSharedPtr fModelProperties;
};
} // namespace corsika::environment
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/particles/ParticleProperties.h>
#include <catch2/catch.hpp>
using namespace corsika::geometry;
using namespace corsika::environment;
using namespace corsika::units::si;
TEST_CASE("HomogeneousMedium") {
NuclearComposition const protonComposition(
std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
std::vector<float>{1.f});
HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
}
add_subdirectory (Utilities)
add_subdirectory (Units)
add_subdirectory (Geometry)
add_subdirectory (Particles)
add_subdirectory (Logging)
add_subdirectory (StackInterface)
add_subdirectory (ProcessSequence)
add_subdirectory (Cascade)
add_subdirectory (Random)
# namespace of library -> location of header files
set (
CORSIKAcascade_NAMESPACE
corsika/cascade
)
# header files of this library
set (
CORSIKAcascade_HEADERS
Cascade.h
)
add_library (CORSIKAcascade INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAcascade ${CORSIKAcascade_NAMESPACE} ${CORSIKAcascade_HEADERS})
#target_link_libraries (
# CORSIKAcascade
# CORSIKAparticles
# CORSIKAunits
# CORSIKAthirdparty # for catch2
# )
# include directive for upstream code
target_include_directories (
CORSIKAcascade
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/>
)
# install library
install (
FILES ${CORSIKAcascade_HEADERS}
DESTINATION include/${CORSIKAcascade_NAMESPACE}
)
# ----------------
# code unit testing
add_executable (
testCascade
testCascade.cc
)
target_link_libraries (
testCascade
# CORSIKAutls
CORSIKArandom
ProcessSibyll
CORSIKAcascade
ProcessStackInspector
CORSIKAstackinterface
CORSIKAprocesses
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
CORSIKAunits
CORSIKAthirdparty # for catch2
)
CORSIKA_ADD_TEST(testCascade)
/**
Here are have to explain the corsika::cascade::Cascade class and its
functionality.
*/
\ No newline at end of file
/*
* (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/environment/Environment.h>
#include <corsika/process/ProcessReturn.h>
#include <corsika/random/RNGManager.h>
#include <corsika/random/UniformRealDistribution.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <cmath>
#include <iostream>
/**
* The cascade namespace assembles all objects needed to simulate full particles cascades.
*/
namespace corsika::cascade {
/**
* \class Cascade
*
* The Cascade class is constructed from template arguments making
* it very versatile. Via the template arguments physics models are
* plugged into the cascade simulation.
*
* <b>Tracking</b> must be a class according to the
* TrackingInterface providing the functions: <code>void
* Init();</code> and <code>auto GetTrack(Particle const& p)</auto>,
* where the latter has a return type of <code>
* geometry::Trajectory<corsika::geometry::Line or Helix> </code>
*
* <b>ProcessList</b> must be a ProcessSequence.
* TimeOfIntersection(corsika::geometry::Line const& line,
*
* <b>Stack</b> is the storage object for particle data, i.e. with
* Particle class type <code>Stack::ParticleType</code>
*
*
*/
template <typename Tracking, typename ProcessList, typename Stack>
class Cascade {
using Particle = typename Stack::ParticleType;
// we only want fully configured objects
Cascade() = delete;
public:
/**
* Cascade class cannot be default constructed, but needs a valid
* list of physics processes for configuration at construct time.
*/
Cascade(corsika::environment::Environment const& env, Tracking& tr, ProcessList& pl,
Stack& stack)
: fEnvironment(env)
, fTracking(tr)
, fProcessSequence(pl)
, fStack(stack) {}
/**
* The Init function is called before the actual cascade simulations.
* All components of the Cascade simulation must be configured here.
*/
void Init() {
fTracking.Init();
fProcessSequence.Init();
fStack.Init();
}
/**
* The Run function is the main simulation loop, which processes
* particles from the Stack until the Stack is empty.
*/
void Run() {
while (!fStack.IsEmpty()) {
while (!fStack.IsEmpty()) {
auto pNext = fStack.GetNextParticle();
Step(pNext);
}
// do cascade equations, which can put new particles on Stack,
// thus, the double loop
// DoCascadeEquations();
}
}
private:
/**
* The Step function is executed for each particle from the
* stack. It will calcualte geometric transport of the particles,
* and apply continuous and stochastic processes to it, which may
* lead to energy losses, scattering, absorption, decays and the
* production of secondary particles.
*
* New particles produced in one step are subject to further
* processing, e.g. thinning, etc.
*/
void Step(Particle& particle) {
using namespace corsika::units::si;
// determine geometric tracking
corsika::setup::Trajectory step = fTracking.GetTrack(particle);
// determine combined total interaction length (inverse)
InverseGrammageType const total_inv_lambda =
fProcessSequence.GetTotalInverseInteractionLength(particle, step);
// sample random exponential step length in grammage
std::exponential_distribution expDist(total_inv_lambda * (1_g / (1_m * 1_m)));
GrammageType const next_interact = (1_g / (1_m * 1_m)) * expDist(fRNG);
std::cout << "total_inv_lambda=" << total_inv_lambda
<< ", next_interact=" << next_interact << std::endl;
// convert next_step from grammage to length
auto const* currentNode =
fEnvironment.GetUniverse()->GetContainingNode(particle.GetPosition());
if (currentNode == &*fEnvironment.GetUniverse()) {
throw std::runtime_error("particle entered void universe");
}
LengthType const distance_interact =
currentNode->GetModelProperties().ArclengthFromGrammage(step, next_interact);
// determine the maximum geometric step length
LengthType const distance_max = fProcessSequence.MaxStepLength(particle, step);
std::cout << "distance_max=" << distance_max << std::endl;
// determine combined total inverse decay time
InverseTimeType const total_inv_lifetime =
fProcessSequence.GetTotalInverseLifetime(particle);
// sample random exponential decay time
std::exponential_distribution expDistDecay(total_inv_lifetime * 1_s);
TimeType const next_decay = 1_s * expDistDecay(fRNG);
std::cout << "total_inv_lifetime=" << total_inv_lifetime
<< ", next_decay=" << next_decay << std::endl;
// convert next_decay from time to length [m]
LengthType const distance_decay = next_decay * particle.GetMomentum().norm() /
particle.GetEnergy() *
corsika::units::constants::c;
// take minimum of geometry, interaction, decay for next step
auto const min_distance =
std::min({distance_interact, distance_decay, distance_max});
std::cout << " move particle by : " << min_distance << std::endl;
// here the particle is actually moved along the trajectory to new position:
// std::visit(corsika::setup::ParticleUpdate<Particle>{particle}, step);
particle.SetPosition(step.PositionFromArclength(min_distance));
// .... also update time, momentum, direction, ...
step.LimitEndTo(min_distance);
// apply all continuous processes on particle + track
corsika::process::EProcessReturn status =
fProcessSequence.DoContinuous(particle, step, fStack);
if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
std::cout << "Cascade: delete absorbed particle " << particle.GetPID() << " "
<< particle.GetEnergy() / 1_GeV << "GeV" << std::endl;
particle.Delete();
return;
}
std::cout << "sth. happening before geometric limit ? "
<< ((min_distance < distance_max) ? "yes" : "no") << std::endl;
if (min_distance < distance_max) { // interaction to happen within geometric limit
// check whether decay or interaction limits this step
if (min_distance == distance_interact) {
std::cout << "collide" << std::endl;
InverseGrammageType const actual_inv_length =
fProcessSequence.GetTotalInverseInteractionLength(particle, step);
corsika::random::UniformRealDistribution<InverseGrammageType> uniDist(
actual_inv_length);
const auto sample_process = uniDist(fRNG);
InverseGrammageType inv_lambda_count = 0. * meter * meter / gram;
fProcessSequence.SelectInteraction(particle, fStack, sample_process,
inv_lambda_count);
} else {
std::cout << "decay" << std::endl;
InverseTimeType const actual_decay_time =
fProcessSequence.GetTotalInverseLifetime(particle);
corsika::random::UniformRealDistribution<InverseTimeType> uniDist(
actual_decay_time);
const auto sample_process = uniDist(fRNG);
InverseTimeType inv_decay_count = 0 / second;
fProcessSequence.SelectDecay(particle, fStack, sample_process, inv_decay_count);
}
}
}
private:
corsika::environment::Environment const& fEnvironment;
Tracking& fTracking;
ProcessList& fProcessSequence;
Stack& fStack;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
};
} // namespace corsika::cascade
#endif
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
namespace cascade;
void Cascade::Step(auto& sequence, Particle& particle) {
double nextStep = sequence.MinStepLength(particle);
Trajectory trajectory = sequence.Transport(particle, nextStep);
sequence.DoContinuous(particle, trajectory);
sequence.DoDiscrete(particle);
}
/*
* (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 <limits>
#include <corsika/environment/Environment.h>
#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/stack/super_stupid/SuperStupidStack.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/setup/SetupStack.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;
corsika::environment::Environment MakeDummyEnv() {
corsika::environment::Environment env; // dummy environment
auto& universe = *(env.GetUniverse());
auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel =
corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_g / (1_m * 1_m * 1_m),
corsika::environment::NuclearComposition(
std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
std::vector<float>{1.}));
universe.AddChild(std::move(theMedium));
return env;
}
class ProcessSplit : public corsika::process::ContinuousProcess<ProcessSplit> {
int fCount = 0;
int fCalls = 0;
HEPEnergyType fEcrit;
public:
ProcessSplit(HEPEnergyType e)
: fEcrit(e) {}
template <typename Particle, typename T>
LengthType MaxStepLength(Particle&, T&) const {
return 1_m;
}
template <typename Particle, typename T, typename Stack>
EProcessReturn DoContinuous(Particle& p, T&, Stack&) {
fCalls++;
HEPEnergyType E = p.GetEnergy();
if (E < fEcrit) {
p.Delete();
fCount++;
} else {
p.SetEnergy(E / 2);
p.AddSecondary(p.GetPID(), E / 2, p.GetMomentum(), p.GetPosition(), p.GetTime());
}
return EProcessReturn::eOk;
}
void Init() {
fCount = 0;
fCalls = 0;
}
int GetCount() const { return fCount; }
int GetCalls() const { return fCalls; }
private:
};
TEST_CASE("Cascade", "[Cascade]") {
corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
rmng.RegisterRandomStream("cascade");
auto env = MakeDummyEnv();
tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
const HEPEnergyType Ecrit = 85_MeV;
ProcessSplit p1(Ecrit);
auto sequence = p0 << p1;
setup::Stack stack;
corsika::cascade::Cascade EAS(env, tracking, sequence, stack);
CoordinateSystem const& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
stack.Clear();
HEPEnergyType E0 = 100_GeV;
stack.AddParticle(
particles::Code::Electron, E0,
corsika::stack::super_stupid::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
Point(rootCS, {0_m, 0_m, 10_km}), 0_ns);
EAS.Init();
EAS.Run();
CHECK(p1.GetCount() == 2048);
CHECK(p1.GetCalls() == 4095);
/*
SECTION("sectionTwo") {
for (int i = 0; i < 0; ++i) {
stack.Clear();
auto particle = stack.NewParticle();
HEPEnergyType E0 = 100_GeV * pow(10, i);
particle.SetEnergy(E0);
EAS.Init();
EAS.Run();
// cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
}
}
*/
}