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 0 additions and 1950 deletions
Documentation/Examples/CMakeLists.txt deleted 100644 → 0
View file @ 0e5bce9d
CORSIKA_ADD_EXAMPLE (helix_example)
target_link_libraries (helix_example CORSIKAgeometry CORSIKAunits)
CORSIKA_ADD_EXAMPLE (particle_list_example)
target_link_libraries (particle_list_example CORSIKAparticles CORSIKAunits CORSIKAprocesses ProcessSibyll ProcessQGSJetII)
CORSIKA_ADD_EXAMPLE (geometry_example)
target_link_libraries (geometry_example CORSIKAgeometry CORSIKAunits)
CORSIKA_ADD_EXAMPLE (logger_example)
target_link_libraries (logger_example CORSIKAunits CORSIKAlogging)
CORSIKA_ADD_EXAMPLE (stack_example)
target_link_libraries (stack_example SuperStupidStack CORSIKAunits
CORSIKAlogging)
# address sanitizer is making this example too slow, so we only do "undefined"
CORSIKA_ADD_EXAMPLE (boundary_example)
target_link_libraries (boundary_example
SuperStupidStack
CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
CORSIKAcascade
ProcessTrackWriter
ProcessParticleCut
ProcessTrackingLine
CORSIKAprocesses
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
CORSIKA_ADD_EXAMPLE (cascade_example)
target_link_libraries (cascade_example
SuperStupidStack
CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
CORSIKAcascade
ProcessEnergyLoss
ProcessTrackWriter
ProcessStackInspector
ProcessTrackingLine
ProcessParticleCut
ProcessHadronicElasticModel
ProcessStackInspector
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
if (Pythia8_FOUND)
CORSIKA_ADD_EXAMPLE (cascade_proton_example)
target_link_libraries (cascade_proton_example
SuperStupidStack
CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
ProcessPythia8
CORSIKAcascade
ProcessEnergyLoss
ProcessTrackWriter
ProcessStackInspector
ProcessTrackingLine
ProcessParticleCut
ProcessHadronicElasticModel
ProcessStackInspector
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
CORSIKA_ADD_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000.)
target_link_libraries (vertical_EAS
SuperStupidStack
CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
ProcessPythia8
ProcessUrQMD
ProcessSwitch
CORSIKAcascade
ProcessCONEXSourceCut
ProcessEnergyLoss
ProcessObservationPlane
ProcessInteractionCounter
ProcessTrackWriter
ProcessTrackingLine
ProcessParticleCut
ProcessOnShellCheck
ProcessStackInspector
ProcessLongitudinalProfile
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
endif()
CORSIKA_ADD_EXAMPLE (stopping_power stopping_power)
target_link_libraries (stopping_power
SuperStupidStack
CORSIKAunits
ProcessEnergyLoss
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
)
CORSIKA_ADD_EXAMPLE (staticsequence_example)
target_link_libraries (staticsequence_example
CORSIKAprocesssequence
CORSIKAunits
CORSIKAgeometry
CORSIKAlogging)
Documentation/Examples/boundary_example.cc deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
#include <corsika/utl/CorsikaFenv.h>
#include <iostream>
#include <limits>
#include <typeinfo>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
template <bool deleteParticle>
struct MyBoundaryCrossingProcess
: public BoundaryCrossingProcess<MyBoundaryCrossingProcess<deleteParticle>> {
MyBoundaryCrossingProcess(std::string const& filename) { fFile.open(filename); }
template <typename Particle>
EProcessReturn DoBoundaryCrossing(Particle& p,
typename Particle::BaseNodeType const& from,
typename Particle::BaseNodeType const& to) {
std::cout << "boundary crossing! from: " << &from << "; to: " << &to << std::endl;
auto const& name = particles::GetName(p.GetPID());
auto const start = p.GetPosition().GetCoordinates();
fFile << name << " " << start[0] / 1_m << ' ' << start[1] / 1_m << ' '
<< start[2] / 1_m << '\n';
if constexpr (deleteParticle) { p.Delete(); }
return EProcessReturn::eOk;
}
void Init() {}
private:
std::ofstream fFile;
};
//
// The example main program for a particle cascade
//
int main() {
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
using EnvType = Environment<setup::IEnvironmentModel>;
EnvType env;
auto& universe = *(env.GetUniverse());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto outerMedium = EnvType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
auto const props =
outerMedium
->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Proton},
std::vector<float>{1.f}));
auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5_km);
innerMedium->SetModelProperties(props);
outerMedium->AddChild(std::move(innerMedium));
universe.AddChild(std::move(outerMedium));
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
process::sibyll::Interaction sibyll;
process::sibyll::Decay decay;
process::particle_cut::ParticleCut cut(20_GeV);
process::track_writer::TrackWriter trackWriter("tracks.dat");
MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
// assemble all processes into an ordered process list
auto sequence = sibyll << decay << cut << boundaryCrossing << trackWriter;
// setup particle stack, and add primary particles
setup::Stack stack;
stack.Clear();
const Code beamCode = Code::Proton;
const HEPMassType mass = particles::GetMass(Code::Proton);
const HEPEnergyType E0 = 50_TeV;
std::uniform_real_distribution distTheta(0., 180.);
std::uniform_real_distribution distPhi(0., 360.);
std::mt19937 rng;
for (int i = 0; i < 100; ++i) {
auto const theta = distTheta(rng);
auto const phi = distPhi(rng);
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
beamCode, E0, plab, pos, 0_ns});
}
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Init();
EAS.Run();
cout << "Result: E0=" << E0 / 1_GeV << endl;
cut.ShowResults();
const HEPEnergyType Efinal =
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
cout << "total energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
}
Documentation/Examples/cascade_example.cc deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/ShowerAxis.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
#include <corsika/utl/CorsikaFenv.h>
#include <iostream>
#include <limits>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
//
// The example main program for a particle cascade
//
int main() {
const LengthType height_atmosphere = 112.8_km;
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
using EnvType = environment::Environment<setup::IEnvironmentModel>;
EnvType env;
auto& universe = *(env.GetUniverse());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto outerMedium = EnvType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
// fraction of oxygen
const float fox = 0.20946;
auto const props =
outerMedium
->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Nitrogen,
particles::Code::Oxygen},
std::vector<float>{1.f - fox, fox}));
auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
innerMedium->SetModelProperties(props);
outerMedium->AddChild(std::move(innerMedium));
universe.AddChild(std::move(outerMedium));
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
const Code beamCode = Code::Nucleus;
const int nuclA = 4;
const int nuclZ = int(nuclA / 2.15 + 0.7);
const HEPMassType mass = GetNucleusMass(nuclA, nuclZ);
const HEPEnergyType E0 = nuclA * 1_TeV;
double theta = 0.;
double phi = 0.;
Point const injectionPos(
rootCS, 0_m, 0_m,
height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe
ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -5000_km}, env};
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
beamCode, E0, plab, injectionPos, 0_ns, nuclA, nuclZ});
}
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
process::sibyll::Decay decay;
// cascade with only HE model ==> HE cut
process::particle_cut::ParticleCut cut(80_GeV);
process::track_writer::TrackWriter trackWriter("tracks.dat");
process::energy_loss::EnergyLoss eLoss{showerAxis};
// assemble all processes into an ordered process list
auto sequence = stackInspect << sibyll << sibyllNuc << decay << eLoss << cut
<< trackWriter;
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Init();
EAS.Run();
eLoss.PrintProfile(); // print longitudinal profile
cut.ShowResults();
const HEPEnergyType Efinal =
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
<< "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
}
Documentation/Examples/cascade_proton_example.cc deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/hadronic_elastic_model/HadronicElasticModel.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/pythia/Interaction.h>
#include <corsika/process/track_writer/TrackWriter.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/random/RNGManager.h>
#include <corsika/utl/CorsikaFenv.h>
#include <iostream>
#include <limits>
#include <typeinfo>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
//
// The example main program for a particle cascade
//
int main() {
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
using EnvType = Environment<setup::IEnvironmentModel>;
EnvType env;
auto& universe = *(env.GetUniverse());
auto theMedium =
EnvType::CreateNode<Sphere>(Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = HomogeneousMedium<IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
NuclearComposition(std::vector<particles::Code>{particles::Code::Hydrogen},
std::vector<float>{(float)1.}));
universe.AddChild(std::move(theMedium));
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
const Code beamCode = Code::Proton;
const HEPMassType mass = particles::Proton::GetMass();
const HEPEnergyType E0 = 100_GeV;
double theta = 0.;
double phi = 0.;
{
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt(Elab * Elab - m * m);
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
Point pos(rootCS, 0_m, 0_m, 0_m);
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
beamCode, E0, plab, pos, 0_ns});
}
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
// process::sibyll::Interaction sibyll(env);
process::pythia::Interaction pythia;
// process::sibyll::NuclearInteraction sibyllNuc(env, sibyll);
// process::sibyll::Decay decay;
process::pythia::Decay decay;
process::particle_cut::ParticleCut cut(20_GeV);
// random::RNGManager::GetInstance().RegisterRandomStream("HadronicElasticModel");
// process::HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
process::track_writer::TrackWriter trackWriter("tracks.dat");
// assemble all processes into an ordered process list
// auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;
auto sequence = pythia << decay << cut << trackWriter << stackInspect;
// cout << "decltype(sequence)=" << type_id_with_cvr<decltype(sequence)>().pretty_name()
// << "\n";
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Init();
EAS.Run();
cout << "Result: E0=" << E0 / 1_GeV << endl;
cut.ShowResults();
const HEPEnergyType Efinal =
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
cout << "total energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
}
Documentation/Examples/logger_example.cc deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/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;
}
Documentation/Examples/particle_list_example.cc deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/qgsjetII/ParticleConversion.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/setup/SetupEnvironment.h>
#include <corsika/units/PhysicalUnits.h>
#include <iomanip>
#include <iostream>
#include <string>
using namespace corsika::units;
using namespace corsika::units::si;
using namespace corsika::particles;
using namespace std;
//
// The example main program for a particle list
//
int main() {
cout << "------------------------------------------"
<< "particles in CORSIKA"
<< "------------------------------------------" << endl;
cout << std::setw(20) << "Name"
<< " | " << std::setw(10) << "PDG-id"
<< " | " << std::setw(10) << "SIBYLL-id"
<< " | " << std::setw(10) << "QGSJETII-id"
<< " | " << std::setw(18) << "PDG-mass (GeV)"
<< " | " << std::setw(18) << "SIBYLL-mass (GeV)"
<< " | " << endl;
cout << std::setw(104) << std::setfill('-') << "-" << endl;
for (auto p : getAllParticles()) {
if (!IsNucleus(p)) {
corsika::process::sibyll::SibyllCode sib_id =
corsika::process::sibyll::ConvertToSibyll(p);
auto const sib_mass =
(sib_id != corsika::process::sibyll::SibyllCode::Unknown
? to_string(corsika::process::sibyll::GetSibyllMass(p) / 1_GeV)
: "--");
auto const qgs_id = corsika::process::qgsjetII::ConvertToQgsjetII(p);
cout << std::setw(20) << std::setfill(' ') << p << " | " << std::setw(10)
<< static_cast<int>(GetPDG(p)) << " | " << std::setw(10)
<< (sib_id != corsika::process::sibyll::SibyllCode::Unknown
? to_string(static_cast<int>(sib_id))
: "--")
<< " | " << std::setw(10)
<< (qgs_id != corsika::process::qgsjetII::QgsjetIICode::Unknown
? to_string(static_cast<int>(qgs_id))
: "--")
<< " | " << std::setw(18) << std::setprecision(5) << GetMass(p) / 1_GeV
<< " | " << std::setw(18) << std::setprecision(5) << sib_mass << " | " << endl;
}
}
cout << std::setw(104) << std::setfill('-') << "-" << endl;
}
Documentation/Examples/stack_example.cc deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/particles/ParticleProperties.h>
#include <corsika/stack/super_stupid/SuperStupidStack.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <cassert>
#include <iomanip>
#include <iostream>
using namespace corsika;
using namespace corsika::units::si;
using namespace corsika::stack;
using namespace corsika::geometry;
using namespace std;
void fill(corsika::stack::super_stupid::SuperStupidStack& s) {
const geometry::CoordinateSystem& rootCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
for (int i = 0; i < 11; ++i) {
s.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Electron, 1.5_GeV * i,
corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, 1_GeV}),
geometry::Point(rootCS, 0_m, 0_m, 0_m), 0_ns});
}
}
void read(corsika::stack::super_stupid::SuperStupidStack& s) {
assert(s.GetSize() == 11); // stack has 11 particles
HEPEnergyType total_energy;
int i = 0;
for (auto& p : s) {
total_energy += p.GetEnergy();
// particles are electrons with 1.5 GeV energy times i
assert(p.GetPID() == particles::Code::Electron);
assert(p.GetEnergy() == 1.5_GeV * (i++));
}
}
int main() {
corsika::stack::super_stupid::SuperStupidStack s;
fill(s);
read(s);
return 0;
}
Documentation/Examples/staticsequence_example.cc deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <array>
#include <iomanip>
#include <iostream>
#include <corsika/process/ProcessSequence.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
using namespace corsika;
using namespace corsika::units::si;
using namespace corsika::process;
using namespace std;
const int nData = 10;
class Process1 : public BaseProcess<Process1> {
public:
Process1() {}
template <typename D, typename T, typename S>
EProcessReturn DoContinuous(D& d, T&, S&) const {
for (int i = 0; i < nData; ++i) d.p[i] += 1;
return EProcessReturn::eOk;
}
};
class Process2 : public BaseProcess<Process2> {
public:
Process2() {}
template <typename D, typename T, typename S>
inline EProcessReturn DoContinuous(D& d, T&, S&) const {
for (int i = 0; i < nData; ++i) d.p[i] -= 0.1 * i;
return EProcessReturn::eOk;
}
};
class Process3 : public BaseProcess<Process3> {
public:
Process3() {}
template <typename D, typename T, typename S>
inline EProcessReturn DoContinuous(D&, T&, S&) const {
return EProcessReturn::eOk;
}
};
class Process4 : public BaseProcess<Process4> {
public:
Process4(const double v)
: fV(v) {}
template <typename D, typename T, typename S>
inline EProcessReturn DoContinuous(D& d, T&, S&) const {
for (int i = 0; i < nData; ++i) d.p[i] *= fV;
return EProcessReturn::eOk;
}
private:
double fV;
};
struct DummyData {
double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
};
struct DummyStack {
void clear() {}
};
struct DummyTrajectory {};
void modular() {
Process1 m1;
Process2 m2;
Process3 m3;
Process4 m4(0.9);
auto sequence = m1 << m2 << m3 << m4;
DummyData particle;
DummyTrajectory track;
const int n = 1000;
for (int i = 0; i < n; ++i) { sequence.DoContinuous(particle, track); }
for (int i = 0; i < nData; ++i) {
// cout << p.p[i] << endl;
// assert(p.p[i] == n-i*100);
}
cout << " done (nothing...) " << endl;
}
int main() {
modular();
return 0;
}
Documentation/Examples/stopping_power.cc deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/ShowerAxis.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/CorsikaFenv.h>
#include <fstream>
#include <iostream>
#include <limits>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::particles;
using namespace corsika::geometry;
using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
//
// This example demonstrates the energy loss of muons as function of beta*gamma (=p/m)
//
int main() {
feenableexcept(FE_INVALID);
// setup environment, geometry
using EnvType = Environment<IMediumModel>;
EnvType env;
env.GetUniverse()->SetModelProperties<HomogeneousMedium<IMediumModel>>(
1_g / cube(1_cm), NuclearComposition{{particles::Code::Unknown}, {1.f}});
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const injectionPos(
rootCS, 0_m, 0_m,
112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
environment::ShowerAxis showerAxis{injectionPos,
Vector<length_d>{rootCS, 0_m, 0_m, 1_m}, env};
process::energy_loss::EnergyLoss eLoss{showerAxis};
setup::Stack stack;
std::ofstream file("dEdX.dat");
file << "# beta*gamma, dE/dX / eV/(g/cm²)" << std::endl;
for (HEPEnergyType E0 = 300_MeV; E0 < 1_PeV; E0 *= 1.05) {
stack.Clear();
const Code beamCode = Code::MuPlus;
const HEPMassType mass = GetMass(beamCode);
double theta = 0.;
double phi = 0.;
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi),
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << endl;
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
beamCode, E0, plab, injectionPos, 0_ns});
auto const p = stack.GetNextParticle();
HEPEnergyType dE = eLoss.TotalEnergyLoss(p, 1_g / square(1_cm));
file << P0 / mass << "\t" << -dE / 1_eV << std::endl;
}
}
Documentation/Examples/vertical_EAS.cc deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/cascade/Cascade.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/FlatExponential.h>
#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/ShowerAxis.h>
#include <corsika/geometry/Plane.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/StackProcess.h>
#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/interaction_counter/InteractionCounter.h>
#include <corsika/process/longitudinal_profile/LongitudinalProfile.h>
#include <corsika/process/observation_plane/ObservationPlane.h>
#include <corsika/process/on_shell_check/OnShellCheck.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/switch_process/SwitchProcess.h>
#include <corsika/process/tracking_line/TrackingLine.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 <corsika/utl/CorsikaFenv.h>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
#include <typeinfo>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
void registerRandomStreams() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
random::RNGManager::GetInstance().RegisterRandomStream("qgran");
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
random::RNGManager::GetInstance().RegisterRandomStream("UrQMD");
random::RNGManager::GetInstance().SeedAll();
}
int main(int argc, char** argv) {
if (argc != 4) {
std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV>" << std::endl;
return 1;
}
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
registerRandomStreams();
// setup environment, geometry
using EnvType = Environment<setup::IEnvironmentModel>;
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
environment::LayeredSphericalAtmosphereBuilder builder{center};
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
builder.addLinearLayer(1e9_cm, 112.8_km);
builder.assemble(env);
// setup particle stack, and add primary particle
setup::Stack stack;
stack.Clear();
const Code beamCode = Code::Nucleus;
unsigned short const A = std::stoi(std::string(argv[1]));
unsigned short Z = std::stoi(std::string(argv[2]));
auto const mass = particles::GetNucleusMass(A, Z);
const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[3]));
double theta = 0.;
auto const thetaRad = theta / 180. * M_PI;
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
};
auto const [px, py, pz] = momentumComponents(thetaRad, P0);
auto plab = corsika::stack::MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << endl;
cout << "input momentum: " << plab.GetComponents() / 1_GeV << ", norm = " << plab.norm()
<< endl;
auto const observationHeight = 0_km + builder.getEarthRadius();
auto const injectionHeight = 112.75_km + builder.getEarthRadius();
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-si::detail::static_pow<2>(sin(thetaRad) * observationHeight) +
si::detail::static_pow<2>(injectionHeight));
Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
Point const injectionPos =
showerCore +
Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
std::cout << "point of injection: " << injectionPos.GetCoordinates() << std::endl;
if (A != 1) {
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
beamCode, E0, plab, injectionPos, 0_ns, A, Z});
} else {
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, injectionPos, 0_ns});
}
std::cout << "shower axis length: " << (showerCore - injectionPos).norm() * 1.02
<< std::endl;
environment::ShowerAxis const showerAxis{injectionPos,
(showerCore - injectionPos) * 1.02, env};
// setup processes, decays and interactions
process::sibyll::Interaction sibyll;
process::interaction_counter::InteractionCounter sibyllCounted(sibyll);
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
process::interaction_counter::InteractionCounter sibyllNucCounted(sibyllNuc);
process::pythia::Decay decayPythia;
// use sibyll decay routine for decays of particles unknown to pythia
process::sibyll::Decay decaySibyll{{
Code::N1440Plus,
Code::N1440MinusBar,
Code::N1440_0,
Code::N1440_0Bar,
Code::N1710Plus,
Code::N1710MinusBar,
Code::N1710_0,
Code::N1710_0Bar,
Code::Pi1300Plus,
Code::Pi1300Minus,
Code::Pi1300_0,
Code::KStar0_1430_0,
Code::KStar0_1430_0Bar,
Code::KStar0_1430_Plus,
Code::KStar0_1430_MinusBar,
}};
decaySibyll.PrintDecayConfig();
process::particle_cut::ParticleCut cut{60_GeV};
process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
process::energy_loss::EnergyLoss eLoss(showerAxis);
process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
process::observation_plane::ObservationPlane observationLevel(obsPlane,
"particles.dat");
process::UrQMD::UrQMD urqmd;
process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
sibyllNuc.Init();
sibyll.Init();
process::conex_source_cut::CONEXSourceCut conexSource(
center, showerAxis, t, injectionHeight, E0,
particles::GetPDG(particles::Code::Proton));
// assemble all processes into an ordered process list
auto sibyllSequence = sibyllNucCounted << sibyllCounted;
process::switch_process::SwitchProcess switchProcess(urqmdCounted, sibyllSequence,
55_GeV);
auto decaySequence = decayPythia << decaySibyll;
auto sequence = switchProcess << reset_particle_mass << decaySequence << conexSource << longprof
<< eLoss << cut << observationLevel;
// define air shower object, run simulation
tracking_line::TrackingLine tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
EAS.Init();
// to fix the point of first interaction, uncomment the following two lines:
// EAS.SetNodes();
// EAS.forceInteraction();
EAS.Run();
eLoss.PrintProfile(); // print longitudinal profile
conexSource.SolveCE();
cut.ShowResults();
const HEPEnergyType Efinal =
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
<< "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
urqmdCounted.GetHistogram();
hists.saveLab("inthist_lab.txt");
hists.saveCMS("inthist_cms.txt");
longprof.save("longprof.txt");
std::ofstream finish("finished");
finish << "run completed without error" << std::endl;
}
Environment/CMakeLists.txt deleted 100644 → 0
View file @ 0e5bce9d
set (
ENVIRONMENT_SOURCES
LayeredSphericalAtmosphereBuilder.cc
ShowerAxis.cc
)
set (
ENVIRONMENT_HEADERS
VolumeTreeNode.h
IMediumModel.h
NuclearComposition.h
HomogeneousMedium.h
InhomogeneousMedium.h
HomogeneousMedium.h
LinearApproximationIntegrator.h
DensityFunction.h
Environment.h
NameModel.h
BaseExponential.h
FlatExponential.h
SlidingPlanarExponential.h
LayeredSphericalAtmosphereBuilder.h
ShowerAxis.h
IMagneticFieldModel.h
UniformMagneticField.h
IRefractiveIndexModel.h
UniformRefractiveIndex.h
)
set (
ENVIRONMENT_NAMESPACE
corsika/environment
)
add_library (CORSIKAenvironment STATIC ${ENVIRONMENT_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAenvironment ${ENVIRONMENT_NAMESPACE} ${ENVIRONMENT_HEADERS})
set_target_properties (
CORSIKAgeometry
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
PUBLIC_HEADER "${ENVIRONMENT_HEADERS}"
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
CORSIKAenvironment
CORSIKAgeometry
CORSIKAparticles
CORSIKAunits
CORSIKArandom
)
target_include_directories (
CORSIKAenvironment
PUBLIC
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include>
)
install (
TARGETS CORSIKAenvironment
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
PUBLIC_HEADER DESTINATION include/${ENVIRONMENT_NAMESPACE}
)
# --------------------
# code unit testing
CORSIKA_ADD_TEST(testEnvironment)
target_link_libraries (
testEnvironment
CORSIKAsetup
CORSIKAenvironment
CORSIKAtesting
)
CORSIKA_ADD_TEST(testShowerAxis)
target_link_libraries (
testShowerAxis
CORSIKAsetup
CORSIKAenvironment
CORSIKAtesting
)
Environment/Environment.h deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Sphere.h>
#include <limits>
namespace corsika::environment {
struct Universe : public corsika::geometry::Sphere {
Universe(corsika::geometry::CoordinateSystem const& pCS)
: corsika::geometry::Sphere(
corsika::geometry::Point{pCS, 0 * corsika::units::si::meter,
0 * corsika::units::si::meter,
0 * corsika::units::si::meter},
corsika::units::si::meter * std::numeric_limits<double>::infinity()) {}
bool Contains(corsika::geometry::Point const&) const override { return true; }
};
template <typename IEnvironmentModel>
class Environment {
public:
using BaseNodeType = VolumeTreeNode<IEnvironmentModel>;
Environment()
: fCoordinateSystem{corsika::geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem()}
, fUniverse(std::make_unique<BaseNodeType>(
std::make_unique<Universe>(fCoordinateSystem))) {}
// using IEnvironmentModel = corsika::setup::IEnvironmentModel;
auto& GetUniverse() { return fUniverse; }
auto const& GetUniverse() const { return fUniverse; }
auto const& GetCoordinateSystem() const { return fCoordinateSystem; }
// factory method for creation of VolumeTreeNodes
template <class TVolumeType, typename... TVolumeArgs>
static auto CreateNode(TVolumeArgs&&... args) {
static_assert(std::is_base_of_v<corsika::geometry::Volume, TVolumeType>,
"unusable type provided, needs to be derived from "
"\"corsika::geometry::Volume\"");
return std::make_unique<BaseNodeType>(
std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
}
private:
corsika::geometry::CoordinateSystem const& fCoordinateSystem;
typename BaseNodeType::VTNUPtr fUniverse;
};
// using SetupBaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
// using SetupEnvironment = Environment<corsika::setup::IEnvironmentModel>;
} // namespace corsika::environment
Environment/FlatExponential.h deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/environment/BaseExponential.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::environment {
//clang-format off
/**
* flat exponential density distribution with
* \f[
* \varrho(r) = \varrho_0 \exp\left( \frac{1}{\lambda} (r - p) \cdot
* \vec{a} \right).
* \f]
* \f$ \vec{a} \f$ denotes the axis and should be normalized to avoid degeneracy
* with the scale parameter \f$ \lambda \f$.
*/
//clang-format on
template <class T>
class FlatExponential : public BaseExponential<FlatExponential<T>>, public T {
geometry::Vector<units::si::dimensionless_d> const fAxis;
NuclearComposition const fNuclComp;
using Base = BaseExponential<FlatExponential<T>>;
public:
FlatExponential(geometry::Point const& vP0,
geometry::Vector<units::si::dimensionless_d> const& vAxis,
units::si::MassDensityType vRho, units::si::LengthType vLambda,
NuclearComposition vNuclComp)
: Base(vP0, vRho, vLambda)
, fAxis(vAxis)
, fNuclComp(vNuclComp) {}
units::si::MassDensityType GetMassDensity(geometry::Point const& vP) const override {
return Base::fRho0 * exp(Base::fInvLambda * (vP - Base::fP0).dot(fAxis));
}
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
units::si::GrammageType IntegratedGrammage(
geometry::Trajectory<geometry::Line> const& vLine,
units::si::LengthType vTo) const override {
return Base::IntegratedGrammage(vLine, vTo, fAxis);
}
units::si::LengthType ArclengthFromGrammage(
geometry::Trajectory<geometry::Line> const& vLine,
units::si::GrammageType vGrammage) const override {
return Base::ArclengthFromGrammage(vLine, vGrammage, fAxis);
}
};
} // namespace corsika::environment
Environment/HomogeneousMedium.h deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/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
Environment/IMediumModel.h deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/units/PhysicalUnits.h>
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;
};
} // namespace corsika::environment
Environment/InhomogeneousMedium.h deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/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>
/**
* A general inhomogeneous medium. The mass density distribution TDensityFunction must be
* a \f$C^2\f$-function.
*/
namespace corsika::environment {
template <class T, class TDensityFunction>
class InhomogeneousMedium : public T {
NuclearComposition const fNuclComp;
TDensityFunction const fDensityFunction;
public:
template <typename... Args>
InhomogeneousMedium(NuclearComposition pNuclComp, Args&&... rhoArgs)
: fNuclComp(pNuclComp)
, fDensityFunction(rhoArgs...) {}
corsika::units::si::MassDensityType GetMassDensity(
corsika::geometry::Point const& p) const override {
return fDensityFunction.EvaluateAt(p);
}
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& pLine,
corsika::units::si::LengthType pTo) const override {
return fDensityFunction.IntegrateGrammage(pLine, pTo);
}
corsika::units::si::LengthType ArclengthFromGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& pLine,
corsika::units::si::GrammageType pGrammage) const override {
return fDensityFunction.ArclengthFromGrammage(pLine, pGrammage);
}
};
} // namespace corsika::environment
Environment/LayeredSphericalAtmosphereBuilder.cc deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/environment/FlatExponential.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/SlidingPlanarExponential.h>
using namespace corsika::environment;
void LayeredSphericalAtmosphereBuilder::checkRadius(units::si::LengthType r) const {
if (r <= previousRadius_) {
throw std::runtime_error("radius must be greater than previous");
}
}
void LayeredSphericalAtmosphereBuilder::setNuclearComposition(
NuclearComposition composition) {
composition_ = std::make_unique<NuclearComposition>(composition);
}
void LayeredSphericalAtmosphereBuilder::addExponentialLayer(
units::si::GrammageType b, units::si::LengthType c,
units::si::LengthType upperBoundary) {
auto const radius = earthRadius_ + upperBoundary;
checkRadius(radius);
previousRadius_ = radius;
auto node = std::make_unique<VolumeTreeNode<IMediumModel>>(
std::make_unique<geometry::Sphere>(center_, radius));
auto const rho0 = b / c;
std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
node->SetModelProperties<SlidingPlanarExponential<IMediumModel>>(
center_, rho0, -c, *composition_, earthRadius_);
layers_.push(std::move(node));
}
void LayeredSphericalAtmosphereBuilder::addLinearLayer(
units::si::LengthType c, units::si::LengthType upperBoundary) {
using namespace units::si;
auto const radius = earthRadius_ + upperBoundary;
checkRadius(radius);
previousRadius_ = radius;
units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
auto const rho0 = b / c;
std::cout << "rho0 = " << rho0;
auto node = std::make_unique<VolumeTreeNode<environment::IMediumModel>>(
std::make_unique<geometry::Sphere>(center_, radius));
node->SetModelProperties<HomogeneousMedium<IMediumModel>>(rho0, *composition_);
layers_.push(std::move(node));
}
Environment<IMediumModel> LayeredSphericalAtmosphereBuilder::assemble() {
Environment<IMediumModel> env;
assemble(env);
return env;
}
void LayeredSphericalAtmosphereBuilder::assemble(Environment<IMediumModel>& env) {
auto& universe = env.GetUniverse();
auto* outmost = universe.get();
while (!layers_.empty()) {
auto l = std::move(layers_.top());
auto* tmp = l.get();
outmost->AddChild(std::move(l));
layers_.pop();
outmost = tmp;
}
}
Environment/LayeredSphericalAtmosphereBuilder.h deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#include <corsika/environment/Environment.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h>
#include <memory>
#include <stack>
namespace corsika::environment {
class LayeredSphericalAtmosphereBuilder {
std::unique_ptr<NuclearComposition> composition_;
geometry::Point center_;
units::si::LengthType previousRadius_{units::si::LengthType::zero()};
units::si::LengthType earthRadius_;
std::stack<VolumeTreeNode<environment::IMediumModel>::VTNUPtr>
layers_; // innermost layer first
void checkRadius(units::si::LengthType) const;
public:
LayeredSphericalAtmosphereBuilder(
corsika::geometry::Point center,
units::si::LengthType earthRadius = units::constants::EarthRadius::Mean)
: center_(center)
, earthRadius_(earthRadius) {}
void setNuclearComposition(NuclearComposition);
void addExponentialLayer(units::si::GrammageType, units::si::LengthType,
units::si::LengthType);
auto size() const { return layers_.size(); }
void addLinearLayer(units::si::LengthType, units::si::LengthType);
void assemble(Environment<IMediumModel>&);
Environment<IMediumModel> assemble();
/**
* Get the current Earth radius.
*/
units::si::LengthType getEarthRadius() const { return earthRadius_; };
};
} // namespace corsika::environment
Environment/LinearApproximationIntegrator.h deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <limits>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
namespace corsika::environment {
template <class TDerived>
class LinearApproximationIntegrator {
auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
public:
auto IntegrateGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
corsika::units::si::LengthType length) const {
auto const c0 = GetImplementation().EvaluateAt(line.GetPosition(0));
auto const c1 = GetImplementation().fRho.FirstDerivative(
line.GetPosition(0), line.NormalizedDirection());
return (c0 + 0.5 * c1 * length) * length;
}
auto ArclengthFromGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
corsika::units::si::GrammageType grammage) const {
auto const c0 = GetImplementation().fRho(line.GetPosition(0));
auto const c1 = GetImplementation().fRho.FirstDerivative(
line.GetPosition(0), line.NormalizedDirection());
return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
}
auto MaximumLength(corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
[[maybe_unused]] double relError) const {
using namespace corsika::units::si;
[[maybe_unused]] auto const c1 = GetImplementation().fRho.SecondDerivative(
line.GetPosition(0), line.NormalizedDirection());
// todo: provide a real, working implementation
return 1_m * std::numeric_limits<double>::infinity();
}
};
} // namespace corsika::environment
Environment/NuclearComposition.h deleted 100644 → 0
View file @ 0e5bce9d
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
*
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
#pragma once
#include <corsika/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
#include <cassert>
#include <functional>
#include <numeric>
#include <random>
#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;
template <class AConstIterator, class BConstIterator>
class WeightProviderIterator {
AConstIterator fAIter;
BConstIterator fBIter;
public:
using value_type = double;
using iterator_category = std::input_iterator_tag;
using pointer = value_type*;
using reference = value_type&;
using difference_type = ptrdiff_t;
WeightProviderIterator(AConstIterator a, BConstIterator b)
: fAIter(a)
, fBIter(b) {}
value_type operator*() const { return ((*fAIter) * (*fBIter)).magnitude(); }
WeightProviderIterator& operator++() { // prefix ++
++fAIter;
++fBIter;
return *this;
}
auto operator==(WeightProviderIterator other) { return fAIter == other.fAIter; }
auto operator!=(WeightProviderIterator other) { return !(*this == other); }
};
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.,
std::plus<double>(), [](auto const compID, auto const fraction) -> double {
if (particles::IsNucleus(compID)) {
return particles::GetNucleusA(compID) * fraction;
} else {
return particles::GetMass(compID) /
units::si::ConvertSIToHEP(units::constants::u) * 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");
}
}
template <typename TFunction>
auto WeightedSum(TFunction func) const {
using ResultQuantity = decltype(func(*fComponents.cbegin()));
auto const prod = [&](auto const compID, auto const fraction) {
return func(compID) * fraction;
};
if constexpr (phys::units::is_quantity_v<ResultQuantity>) {
return std::inner_product(
fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
ResultQuantity::zero(), // .zero() is defined for quantity types only
std::plus<ResultQuantity>(), prod);
} else {
return std::inner_product(
fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
ResultQuantity(0), // in other cases we have to use a bare 0
std::plus<ResultQuantity>(), prod);
}
}
auto size() const { return fNumberFractions.size(); }
auto const& GetFractions() const { return fNumberFractions; }
auto const& GetComponents() const { return fComponents; }
auto const GetAverageMassNumber() const { return fAvgMassNumber; }
template <class TRNG>
corsika::particles::Code SampleTarget(
std::vector<corsika::units::si::CrossSectionType> const& sigma,
TRNG& randomStream) const {
using namespace corsika::units::si;
assert(sigma.size() == fNumberFractions.size());
std::discrete_distribution channelDist(
WeightProviderIterator<decltype(fNumberFractions.begin()),
decltype(sigma.begin())>(fNumberFractions.begin(),
sigma.begin()),
WeightProviderIterator<decltype(fNumberFractions.begin()),
decltype(sigma.end())>(fNumberFractions.end(),
sigma.end()));
auto const iChannel = channelDist(randomStream);
return fComponents[iChannel];
}
};
} // namespace corsika::environment

data privacy - impressum