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
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 2391 deletions
Documentation/Examples/CMakeLists.txt deleted 100644 → 0
View file @ dbf6e995
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 (stack_example)
target_link_libraries (stack_example CORSIKAsetup CORSIKAunits)
# address sanitizer is making this example too slow, so we only do "undefined"
CORSIKA_ADD_EXAMPLE (boundary_example)
target_link_libraries (boundary_example
CORSIKAsetup
CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
CORSIKAcascade
ProcessTrackWriter
ProcessParticleCut
ProcessTrackingLine
ProcessPythia8
CORSIKAprocesses
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
C8::ext::boost # boost::histogram
)
CORSIKA_ADD_EXAMPLE (cascade_example)
target_link_libraries (cascade_example
CORSIKAsetup
CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
ProcessProposal
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
CORSIKAsetup
CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
ProcessPythia8
ProcessUrQMD
CORSIKAcascade
ProcessEnergyLoss
ProcessTrackWriter
ProcessStackInspector
ProcessTrackingLine
ProcessParticleCut
ProcessOnShellCheck
ProcessHadronicElasticModel
ProcessStackInspector
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
CORSIKA_ADD_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000.)
target_link_libraries (vertical_EAS
CORSIKAsetup
CORSIKAunits
CORSIKAlogging
CORSIKArandom
CORSIKAhistory
ProcessSibyll
ProcessPythia8
ProcessUrQMD
CORSIKAcascade
ProcessProposal
ProcessPythia8
ProcessObservationPlane
ProcessInteractionCounter
ProcessTrackWriter
ProcessEnergyLoss
ProcessTrackingLine
ProcessParticleCut
ProcessOnShellCheck
ProcessStackInspector
ProcessLongitudinalProfile
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
CORSIKAhistory # for HistoryObservationPlane
)
CORSIKA_ADD_EXAMPLE (hybrid_MC RUN_OPTIONS 4 2 10000.)
target_link_libraries (hybrid_MC
CORSIKAsetup
CORSIKAunits
CORSIKAlogging
CORSIKArandom
CORSIKAhistory
ProcessCONEXSourceCut
ProcessInteractionCounter
ProcessSibyll
ProcessPythia8
ProcessUrQMD
CORSIKAcascade
ProcessPythia8
ProcessObservationPlane
ProcessInteractionCounter
ProcessTrackWriter
ProcessEnergyLoss
ProcessTrackingLine
ProcessParticleCut
ProcessOnShellCheck
ProcessStackInspector
ProcessLongitudinalProfile
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
CORSIKAhistory # for HistoryObservationPlane
)
endif()
CORSIKA_ADD_EXAMPLE (stopping_power stopping_power)
target_link_libraries (stopping_power
CORSIKAsetup
CORSIKAunits
CORSIKAlogging
ProcessEnergyLoss
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
)
CORSIKA_ADD_EXAMPLE (staticsequence_example)
target_link_libraries (staticsequence_example
CORSIKAprocesssequence
CORSIKAunits
CORSIKAgeometry
CORSIKAlogging)
CORSIKA_ADD_EXAMPLE (em_shower RUN_OPTIONS "100.")
target_link_libraries (em_shower
SuperStupidStack
CORSIKAunits
CORSIKAlogging
CORSIKArandom
ProcessSibyll
ProcessPythia8
ProcessUrQMD
CORSIKAcascade
ProcessEnergyLoss
ProcessObservationPlane
ProcessInteractionCounter
ProcessTrackWriter
ProcessProposal
ProcessTrackingLine
ProcessParticleCut
ProcessOnShellCheck
ProcessStackInspector
ProcessLongitudinalProfile
CORSIKAprocesses
CORSIKAcascade
CORSIKAparticles
CORSIKAgeometry
CORSIKAenvironment
CORSIKAprocesssequence
)
Documentation/Examples/cascade_example.cc deleted 100644 → 0
View file @ dbf6e995
/*
* (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/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 <corsika/logging/Logging.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::geometry;
using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
//
// The example main program for a particle cascade
//
int main() {
logging::SetLevel(logging::level::info);
std::cout << "cascade_example" << std::endl;
const LengthType height_atmosphere = 112.8_km;
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
setup::Environment env;
auto& universe = *(env.GetUniverse());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto world = setup::Environment::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel =
environment::MediumPropertyModel<environment::UniformMagneticField<
environment::HomogeneousMedium<setup::EnvironmentInterface>>>;
// fraction of oxygen
const float fox = 0.20946;
auto const props = world->SetModelProperties<MyHomogeneousModel>(
environment::Medium::AirDry1Atm, Vector(rootCS, 0_T, 0_T, 0_T),
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 =
setup::Environment::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
innerMedium->SetModelProperties(props);
world->AddChild(std::move(innerMedium));
universe.AddChild(std::move(world));
// 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
setup::Tracking tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
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, true, true);
process::track_writer::TrackWriter trackWriter("tracks.dat");
process::energy_loss::EnergyLoss eLoss{showerAxis, cut.GetECut()};
// assemble all processes into an ordered process list
auto sequence =
process::sequence(stackInspect, sibyll, sibyllNuc, decay, eLoss, cut, trackWriter);
// define air shower object, run simulation
cascade::Cascade EAS(env, tracking, sequence, stack);
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;
cut.Reset();
}
Documentation/Examples/cascade_proton_example.cc deleted 100644 → 0
View file @ dbf6e995
/*
* (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/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 <corsika/logging/Logging.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::geometry;
using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
//
// The example main program for a particle cascade
//
int main() {
logging::SetLevel(logging::level::trace);
std::cout << "cascade_proton_example" << std::endl;
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
using EnvType = setup::Environment;
EnvType env;
auto& universe = *(env.GetUniverse());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
auto theMedium = EnvType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel =
environment::MediumPropertyModel<environment::UniformMagneticField<
environment::HomogeneousMedium<setup::EnvironmentInterface>>>;
theMedium->SetModelProperties<MyHomogeneousModel>(
environment::Medium::AirDry1Atm, geometry::Vector(rootCS, 0_T, 0_T, 1_T),
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));
// 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
setup::Tracking tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
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, true, true);
// 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 = process::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.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/em_shower.cc deleted 100644 → 0
View file @ dbf6e995
/*
* (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/cascade/Cascade.h>
#include <corsika/environment/Environment.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/longitudinal_profile/LongitudinalProfile.h>
#include <corsika/process/observation_plane/ObservationPlane.h>
#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/proposal/ContinuousProcess.h>
#include <corsika/process/proposal/Interaction.h>
#include <corsika/process/track_writer/TrackWriter.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 <corsika/process/interaction_counter/InteractionCounter.hpp>
#include <corsika/logging/Logging.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::geometry;
using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
void registerRandomStreams() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
random::RNGManager::GetInstance().RegisterRandomStream("proposal");
random::RNGManager::GetInstance().SeedAll();
}
template <typename T>
using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
int main(int argc, char** argv) {
logging::SetLevel(logging::level::info);
if (argc != 2) {
std::cerr << "usage: em_shower <energy/GeV>" << std::endl;
return 1;
}
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
registerRandomStreams();
// setup environment, geometry
using EnvType = setup::Environment;
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
auto builder = environment::make_layered_spherical_atmosphere_builder<
setup::EnvironmentInterface,
MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
environment::Medium::AirDry1Atm,
geometry::Vector{rootCS, 0_T, 0_T, 1_T});
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::Electron;
auto const mass = particles::GetMass(beamCode);
const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
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 = 1.4_km + builder.getEarthRadius();
auto const injectionHeight = 112.75_km + builder.getEarthRadius();
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
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;
stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
beamCode, 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
// PROPOSAL processs proposal{...};
process::particle_cut::ParticleCut cut(10_GeV, false, true);
process::proposal::Interaction proposal(env, cut.GetECut());
process::proposal::ContinuousProcess em_continuous(env, cut.GetECut());
process::interaction_counter::InteractionCounter proposalCounted(proposal);
process::track_writer::TrackWriter trackWriter("tracks.dat");
// long. profile; columns for gamma, e+, e- still need to be added
process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
process::observation_plane::ObservationPlane observationLevel(
obsPlane, Vector<dimensionless_d>(rootCS, {1., 0., 0.}), "particles.dat");
auto sequence = process::sequence(proposalCounted, em_continuous, longprof, cut,
observationLevel, trackWriter);
// define air shower object, run simulation
setup::Tracking tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
// to fix the point of first interaction, uncomment the following two lines:
// EAS.SetNodes();
// EAS.forceInteraction();
EAS.Run();
cut.ShowResults();
em_continuous.showResults();
observationLevel.ShowResults();
const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
cut.GetEmEnergy() + em_continuous.energyLost() +
observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.Reset();
cut.Reset();
em_continuous.reset();
auto const hists = proposalCounted.GetHistogram();
hists.saveLab("inthist_lab_emShower.npz");
hists.saveCMS("inthist_cms_emShower.npz");
longprof.save("longprof_emShower.txt");
}
Documentation/Examples/hybrid_MC.cc deleted 100644 → 0
View file @ dbf6e995
/*
* (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.
*/
/* clang-format off */
// InteractionCounter used boost/histogram, which
// fails if boost/type_traits have been included before. Thus, we have
// to include it first...
#include <corsika/process/interaction_counter/InteractionCounter.hpp>
/* clang-format on */
#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/logging/Logging.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/SwitchProcessSequence.h>
#include <corsika/process/StackProcess.h>
#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
#include <corsika/process/energy_loss/EnergyLoss.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/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>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::geometry;
using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
void registerRandomStreams(const int seed) {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
random::RNGManager::GetInstance().RegisterRandomStream("qgsjet");
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
random::RNGManager::GetInstance().RegisterRandomStream("urqmd");
random::RNGManager::GetInstance().RegisterRandomStream("proposal");
if (seed == 0)
random::RNGManager::GetInstance().SeedAll();
else
random::RNGManager::GetInstance().SeedAll(seed);
}
template <typename T>
using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
int main(int argc, char** argv) {
logging::SetLevel(logging::level::info);
C8LOG_INFO("hybrid_MC");
if (argc < 4) {
std::cerr << "usage: hybrid_MC <A> <Z> <energy/GeV> [seed]" << std::endl;
std::cerr << " if no seed is given, a random seed is chosen" << std::endl;
return 1;
}
feenableexcept(FE_INVALID);
int seed = 0;
if (argc > 4) seed = std::stoi(std::string(argv[4]));
// initialize random number sequence(s)
registerRandomStreams(seed);
// setup environment, geometry
using EnvType = setup::Environment;
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
auto builder = environment::make_layered_spherical_atmosphere_builder<
setup::EnvironmentInterface,
MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
environment::Medium::AirDry1Atm,
geometry::Vector{rootCS, 0_T, 0_T, 1_T});
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(-units::static_pow<2>(sin(thetaRad) * observationHeight) +
units::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, false, true};
process::energy_loss::EnergyLoss eLoss{showerAxis, cut.GetECut()};
corsika::process::conex_source_cut::CONEXSourceCut conex(
center, showerAxis, t, injectionHeight, E0,
particles::GetPDG(particles::Code::Proton));
process::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
process::observation_plane::ObservationPlane observationLevel(
obsPlane, Vector<dimensionless_d>(rootCS, {1., 0., 0.}), "particles.dat");
process::UrQMD::UrQMD urqmd;
process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
// assemble all processes into an ordered process list
struct EnergySwitch {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
process::SwitchResult operator()(const setup::Stack::ParticleType& p) {
if (p.GetEnergy() < cutE_)
return process::SwitchResult::First;
else
return process::SwitchResult::Second;
}
};
auto hadronSequence =
process::select(urqmdCounted, process::sequence(sibyllNucCounted, sibyllCounted),
EnergySwitch(55_GeV));
auto decaySequence = process::sequence(decayPythia, decaySibyll);
auto sequence = process::sequence(hadronSequence, reset_particle_mass, decaySequence,
eLoss, cut, conex, longprof, observationLevel);
// define air shower object, run simulation
setup::Tracking tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
// to fix the point of first interaction, uncomment the following two lines:
// EAS.SetNodes();
// EAS.forceInteraction();
EAS.Run();
cut.ShowResults();
eLoss.showResults();
observationLevel.ShowResults();
const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
cut.GetEmEnergy() + eLoss.energyLost() +
observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.Reset();
cut.Reset();
eLoss.reset();
conex.SolveCE();
auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
urqmdCounted.GetHistogram();
hists.saveLab("inthist_lab.txt");
hists.saveCMS("inthist_cms.txt");
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;
}
Documentation/Examples/particle_list_example.cc deleted 100644 → 0
View file @ dbf6e995
/*
* (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() {
std::cout << "particle_list_example" << std::endl;
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 @ dbf6e995
/*
* (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.getEntries() == 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() {
std::cout << "stack_example" << std::endl;
corsika::stack::super_stupid::SuperStupidStack s;
fill(s);
read(s);
return 0;
}
Documentation/Examples/staticsequence_example.cc deleted 100644 → 0
View file @ dbf6e995
/*
* (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 ContinuousProcess<Process1> {
public:
Process1() {}
template <typename D, typename T>
EProcessReturn DoContinuous(D& d, T&) const {
for (int i = 0; i < nData; ++i) d.p[i] += 1;
return EProcessReturn::eOk;
}
};
class Process2 : public ContinuousProcess<Process2> {
public:
Process2() {}
template <typename D, typename T>
inline EProcessReturn DoContinuous(D& d, T&) const {
for (int i = 0; i < nData; ++i) d.p[i] -= 0.1 * i;
return EProcessReturn::eOk;
}
};
class Process3 : public ContinuousProcess<Process3> {
public:
Process3() {}
template <typename D, typename T>
inline EProcessReturn DoContinuous(D&, T&) const {
return EProcessReturn::eOk;
}
};
class Process4 : public ContinuousProcess<Process4> {
public:
Process4(const double v)
: fV(v) {}
template <typename D, typename T>
inline EProcessReturn DoContinuous(D& d, T&) 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() {
// = 0
Process1 m1; // + 1.0
Process2 m2; // - (0.1*i)
Process3 m3; // * 1.0
Process4 m4(1.5); // * 1.5
auto sequence = process::sequence(m1, m2, m3, m4);
DummyData particle;
DummyTrajectory track;
double check[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
const int nEv = 10;
for (int iEv = 0; iEv < nEv; ++iEv) {
sequence.DoContinuous(particle, track);
for (int i = 0; i < nData; ++i) {
check[i] += 1. - 0.1 * i;
check[i] *= 1.5;
}
}
for (int i = 0; i < nData; ++i) { assert(particle.p[i] == check[i]); }
cout << " done (checking...) " << endl;
}
int main() {
std::cout << "staticsequence_example" << std::endl;
modular();
return 0;
}
Documentation/Examples/stopping_power.cc deleted 100644 → 0
View file @ dbf6e995
/*
* (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() {
std::cout << "stopping_power" << std::endl;
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, 300_MeV};
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 @ dbf6e995
/*
* (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.
*/
/* clang-format off */
// InteractionCounter used boost/histogram, which
// fails if boost/type_traits have been included before. Thus, we have
// to include it first...
#include <corsika/process/interaction_counter/InteractionCounter.hpp>
/* clang-format on */
#include <corsika/cascade/Cascade.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/FlatExponential.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/IMagneticFieldModel.h>
#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/ShowerAxis.h>
#include <corsika/environment/SlidingPlanarExponential.h>
#include <corsika/environment/UniformMagneticField.h>
#include <corsika/geometry/Plane.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/logging/Logging.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/SwitchProcessSequence.h>
#include <corsika/process/StackProcess.h>
#include <corsika/process/energy_loss/EnergyLoss.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/track_writer/TrackWriter.h>
#include <corsika/process/proposal/ContinuousProcess.h>
#include <corsika/process/proposal/Interaction.h>
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.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>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::geometry;
using namespace corsika::environment;
using namespace std;
using namespace corsika::units::si;
using Particle = setup::Stack::StackIterator;
void registerRandomStreams(const int seed) {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
random::RNGManager::GetInstance().RegisterRandomStream("qgsjet");
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
random::RNGManager::GetInstance().RegisterRandomStream("urqmd");
random::RNGManager::GetInstance().RegisterRandomStream("proposal");
if (seed == 0)
random::RNGManager::GetInstance().SeedAll();
else
random::RNGManager::GetInstance().SeedAll(seed);
}
template <typename T>
using MyExtraEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
int main(int argc, char** argv) {
logging::SetLevel(logging::level::trace);
C8LOG_INFO("vertical_EAS");
if (argc < 4) {
std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> [seed]" << std::endl;
std::cerr << " if no seed is given, a random seed is chosen" << std::endl;
return 1;
}
feenableexcept(FE_INVALID);
int seed = 0;
if (argc > 4) seed = std::stoi(std::string(argv[4]));
// initialize random number sequence(s)
registerRandomStreams(seed);
// setup environment, geometry
using EnvType = setup::Environment;
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
auto builder = environment::make_layered_spherical_atmosphere_builder<
setup::EnvironmentInterface,
MyExtraEnv>::create(center, units::constants::EarthRadius::Mean,
environment::Medium::AirDry1Atm,
geometry::Vector{rootCS, 0_T, 5000_mT, 0_T});
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(-units::static_pow<2>(sin(thetaRad) * observationHeight) +
units::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::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
} else {
if (Z == 1) {
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});
} else if (Z == 0) {
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType>{particles::Code::Neutron, E0,
plab, injectionPos, 0_ns});
} else {
std::cerr << "illegal parameters" << std::endl;
return EXIT_FAILURE;
}
}
// we make the axis much longer than the inj-core distance since the
// profile will go beyond the core, depending on zenith angle
std::cout << "shower axis length: " << (showerCore - injectionPos).norm() * 1.5
<< std::endl;
environment::ShowerAxis const showerAxis{injectionPos,
(showerCore - injectionPos) * 1.5, env};
// setup processes, decays and interactions
process::particle_cut::ParticleCut cut{60_GeV, false, true};
process::proposal::Interaction proposal(env, cut.GetECut());
process::proposal::ContinuousProcess em_continuous(env, cut.GetECut());
process::interaction_counter::InteractionCounter proposalCounted(proposal);
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::on_shell_check::OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
process::track_writer::TrackWriter trackWriter("tracks.dat");
process::longitudinal_profile::LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, Vector<dimensionless_d>(rootCS, {0., 0., 1.}));
process::observation_plane::ObservationPlane observationLevel(
obsPlane, Vector<dimensionless_d>(rootCS, {1., 0., 0.}), "particles.dat");
process::UrQMD::UrQMD urqmd;
process::interaction_counter::InteractionCounter urqmdCounted{urqmd};
// assemble all processes into an ordered process list
struct EnergySwitch {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
process::SwitchResult operator()(const Particle& p) {
if (p.GetEnergy() < cutE_)
return process::SwitchResult::First;
else
return process::SwitchResult::Second;
}
};
auto hadronSequence =
process::select(urqmdCounted, process::sequence(sibyllNucCounted, sibyllCounted),
EnergySwitch(55_GeV));
auto decaySequence = process::sequence(decayPythia, decaySibyll);
stack_inspector::StackInspector<setup::Stack> stackInspect(1000, false, E0);
auto sequence =
process::sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
proposalCounted, em_continuous, cut, trackWriter, observationLevel, longprof);
// define air shower object, run simulation
setup::Tracking tracking;
cascade::Cascade EAS(env, tracking, sequence, stack);
// to fix the point of first interaction, uncomment the following two lines:
// EAS.forceInteraction();
EAS.Run();
cut.ShowResults();
em_continuous.showResults();
observationLevel.ShowResults();
const HEPEnergyType Efinal = cut.GetCutEnergy() + cut.GetInvEnergy() +
cut.GetEmEnergy() + em_continuous.energyLost() +
observationLevel.GetEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.Reset();
cut.Reset();
em_continuous.reset();
auto const hists = sibyllCounted.GetHistogram() + sibyllNucCounted.GetHistogram() +
urqmdCounted.GetHistogram() + proposalCounted.GetHistogram();
hists.saveLab("inthist_lab_verticalEAS.npz");
hists.saveCMS("inthist_cms_verticalEAS.npz");
longprof.save("longprof_verticalEAS.txt");
}
Environment/CMakeLists.txt deleted 100644 → 0
View file @ dbf6e995
add_custom_command (
OUTPUT ${PROJECT_BINARY_DIR}/Environment/GeneratedMediaProperties.inc
COMMAND ${PROJECT_SOURCE_DIR}/Environment/readProperties.py
${PROJECT_SOURCE_DIR}/Environment/properties8.dat
DEPENDS readProperties.py
properties8.dat
WORKING_DIRECTORY
${PROJECT_BINARY_DIR}/Environment
COMMENT "Read NIST properties8 data file and produce C++ source code GeneratedMediaProperties.inc"
VERBATIM
)
set (
ENVIRONMENT_SOURCES
ShowerAxis.cc
)
set (
ENVIRONMENT_HEADERS
VolumeTreeNode.h
IEmpty.hpp
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
NoMagneticField.h
IRefractiveIndexModel.h
UniformRefractiveIndex.h
IMediumPropertyModel.h
MediumProperties.h
MediumPropertyModel.h
${PROJECT_BINARY_DIR}/Environment/GeneratedMediaProperties.inc # this one is auto-generated
)
set (
ENVIRONMENT_NAMESPACE
corsika/environment
)
add_library (CORSIKAenvironment STATIC ${ENVIRONMENT_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAenvironment ${ENVIRONMENT_NAMESPACE} ${ENVIRONMENT_HEADERS})
set_target_properties (
CORSIKAenvironment
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
PUBLIC_HEADER "${ENVIRONMENT_HEADERS}"
)
# ....................................................
# since GeneratedMediaProperties.inc is an automatically produced file in the build directory,
# create a symbolic link into the source tree, so that it can be found and edited more easily
# this is not needed for the build to succeed! .......
add_custom_command (
OUTPUT ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedParticleProperties.inc
COMMAND ${CMAKE_COMMAND} -E create_symlink ${PROJECT_BINARY_DIR}/include/corsika/environment/GeneratedMediaProperties.inc ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedMediaProperties.inc
COMMENT "Generate link in source-dir: ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedMediaProperties.inc"
)
add_custom_target (SourceDirLinkMedia DEPENDS ${PROJECT_BINARY_DIR}/Environment/GeneratedMediaProperties.inc)
add_dependencies (CORSIKAenvironment SourceDirLinkMedia)
# .....................................................
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
CORSIKAenvironment
CORSIKAgeometry
CORSIKAparticles
CORSIKAunits
CORSIKAlogging
)
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 @ dbf6e995
/*
* (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/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))) {}
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;
};
} // namespace corsika::environment
Environment/FlatExponential.h deleted 100644 → 0
View file @ dbf6e995
/*
* (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/particles/ParticleProperties.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/setup/SetupTrajectory.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(setup::Trajectory const& vLine,
units::si::LengthType vTo) const override {
return Base::IntegratedGrammage(vLine, vTo, fAxis);
}
units::si::LengthType ArclengthFromGrammage(
setup::Trajectory 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 @ dbf6e995
/*
* (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 <corsika/setup/SetupTrajectory.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::setup::Trajectory const&,
corsika::units::si::LengthType pTo) const override {
using namespace corsika::units::si;
return pTo * fDensity;
}
corsika::units::si::LengthType ArclengthFromGrammage(
corsika::setup::Trajectory const&,
corsika::units::si::GrammageType pGrammage) const override {
return pGrammage / fDensity;
}
};
} // namespace corsika::environment
Environment/IMediumModel.h deleted 100644 → 0
View file @ dbf6e995
/*
* (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/units/PhysicalUnits.h>
#include <corsika/setup/SetupTrajectory.h>
namespace corsika::environment {
class IMediumModel {
public:
virtual ~IMediumModel() = default; // LCOV_EXCL_LINE
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::setup::Trajectory const&,
corsika::units::si::LengthType) const = 0;
virtual corsika::units::si::LengthType ArclengthFromGrammage(
corsika::setup::Trajectory const&,
corsika::units::si::GrammageType) const = 0;
virtual NuclearComposition const& GetNuclearComposition() const = 0;
};
} // namespace corsika::environment
Environment/InhomogeneousMedium.h deleted 100644 → 0
View file @ dbf6e995
/*
* (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::setup::Trajectory const& pLine,
corsika::units::si::LengthType pTo) const override {
return fDensityFunction.IntegrateGrammage(pLine, pTo);
}
corsika::units::si::LengthType ArclengthFromGrammage(
corsika::setup::Trajectory const& pLine,
corsika::units::si::GrammageType pGrammage) const override {
return fDensityFunction.ArclengthFromGrammage(pLine, pGrammage);
}
};
} // namespace corsika::environment
Environment/LayeredSphericalAtmosphereBuilder.h deleted 100644 → 0
View file @ dbf6e995
/*
* (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/FlatExponential.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/SlidingPlanarExponential.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h>
#include <functional>
#include <memory>
#include <stack>
#include <tuple>
#include <type_traits>
namespace corsika::environment {
// forward-decl
template <typename TMediumInterface, template <typename> typename MExtraEnvirnoment>
struct make_layered_spherical_atmosphere_builder;
/**
* Helper class to setup concentric spheres of layered atmosphere
* with spcified density profiles (exponential, linear, ...).
*
* This can be used most importantly to replicate CORSIKA7
* atmospheres.
*
* Each layer by definition has a density profile and a (constant)
* nuclear composition model.
*
*/
namespace detail {
struct NoExtraModelInner {};
template <typename M>
struct NoExtraModel {};
template <template <typename> typename M>
struct has_extra_models : std::true_type {};
template <>
struct has_extra_models<NoExtraModel> : std::false_type {};
} // namespace detail
template <typename TMediumInterface = environment::IMediumModel,
template <typename> typename TMediumModelExtra = detail::NoExtraModel,
typename... TModelArgs>
class LayeredSphericalAtmosphereBuilder {
std::unique_ptr<NuclearComposition> composition_;
geometry::Point center_;
units::si::LengthType previousRadius_{units::si::LengthType::zero()};
units::si::LengthType earthRadius_;
std::tuple<TModelArgs...> const additionalModelArgs_;
std::stack<typename VolumeTreeNode<TMediumInterface>::VTNUPtr>
layers_; // innermost layer first
void checkRadius(units::si::LengthType r) const {
if (r <= previousRadius_) {
throw std::runtime_error("radius must be greater than previous");
}
}
LayeredSphericalAtmosphereBuilder() = delete;
LayeredSphericalAtmosphereBuilder(const LayeredSphericalAtmosphereBuilder&) = delete;
LayeredSphericalAtmosphereBuilder(const LayeredSphericalAtmosphereBuilder&&) = delete;
LayeredSphericalAtmosphereBuilder& operator=(
const LayeredSphericalAtmosphereBuilder&) = delete;
// friend, to allow construction
template <typename, template <typename> typename>
friend struct make_layered_spherical_atmosphere_builder;
protected:
LayeredSphericalAtmosphereBuilder(TModelArgs... args, corsika::geometry::Point center,
units::si::LengthType earthRadius)
: center_(center)
, earthRadius_(earthRadius)
, additionalModelArgs_{args...} {}
public:
void setNuclearComposition(NuclearComposition composition) {
composition_ = std::make_unique<NuclearComposition>(composition);
}
void addExponentialLayer(units::si::GrammageType b, units::si::LengthType c,
units::si::LengthType upperBoundary) {
using namespace units::si;
auto const radius = earthRadius_ + upperBoundary;
checkRadius(radius);
previousRadius_ = radius;
auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
std::make_unique<geometry::Sphere>(center_, radius));
auto const rho0 = b / c;
std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
// helper lambda in which the last 5 arguments to make_shared<...> are bound
auto lastBound = [&](auto... argPack) {
return std::make_shared<
TMediumModelExtra<environment::SlidingPlanarExponential<TMediumInterface>>>(
argPack..., center_, rho0, -c, *composition_, earthRadius_);
};
// now unpack the additional arguments
auto model = std::apply(lastBound, additionalModelArgs_);
node->SetModelProperties(std::move(model));
} else {
node->template SetModelProperties<SlidingPlanarExponential<TMediumInterface>>(
center_, rho0, -c, *composition_, earthRadius_);
}
layers_.push(std::move(node));
}
void addLinearLayer(units::si::LengthType c, units::si::LengthType upperBoundary) {
using namespace units::si;
auto const radius = earthRadius_ + upperBoundary;
checkRadius(radius);
previousRadius_ = radius;
auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
std::make_unique<geometry::Sphere>(center_, radius));
units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
auto const rho0 = b / c;
std::cout << "rho0 = " << rho0;
if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
// helper lambda in which the last 2 arguments to make_shared<...> are bound
auto lastBound = [&](auto... argPack) {
return std::make_shared<
TMediumModelExtra<environment::HomogeneousMedium<TMediumInterface>>>(
argPack..., rho0, *composition_);
};
// now unpack the additional arguments
auto model = std::apply(lastBound, additionalModelArgs_);
node->SetModelProperties(std::move(model));
} else {
node->template SetModelProperties<
environment::HomogeneousMedium<TMediumInterface>>(rho0, *composition_);
}
layers_.push(std::move(node));
}
int size() const { return layers_.size(); }
void assemble(Environment<TMediumInterface>& 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<TMediumInterface> assemble() {
Environment<TMediumInterface> env;
assemble(env);
return env;
}
/**
* Get the current Earth radius.
*/
units::si::LengthType getEarthRadius() const { return earthRadius_; }
}; // end class LayeredSphericalAtmosphereBuilder
/**
* \class make_layered_spherical_atmosphere_builder
*
* Helper class to create LayeredSphericalAtmosphereBuilder, the
* extra environment models have to be passed as template-template
* argument to make_layered_spherical_atmosphere_builder, the member
* function `create` does then take an unspecified number of extra
* parameters to internalize those models for all layers later
* produced.
**/
template <typename TMediumInterface = environment::IMediumModel,
template <typename> typename MExtraEnvirnoment = detail::NoExtraModel>
struct make_layered_spherical_atmosphere_builder {
template <typename... TArgs>
static auto create(geometry::Point const& center, units::si::LengthType earthRadius,
TArgs... args) {
return environment::LayeredSphericalAtmosphereBuilder<TMediumInterface,
MExtraEnvirnoment, TArgs...>{
std::forward<TArgs>(args)..., center, earthRadius};
}
};
} // namespace corsika::environment
Environment/LinearApproximationIntegrator.h deleted 100644 → 0
View file @ dbf6e995
/*
* (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/setup/SetupTrajectory.h>
namespace corsika::environment {
template <class TDerived>
class LinearApproximationIntegrator {
auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
public:
auto IntegrateGrammage(corsika::setup::Trajectory 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.GetDirection(0));
return (c0 + 0.5 * c1 * length) * length;
}
auto ArclengthFromGrammage(corsika::setup::Trajectory 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.GetDirection(0));
return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
}
auto MaximumLength(corsika::setup::Trajectory const& line,
[[maybe_unused]] double relError) const {
using namespace corsika::units::si;
[[maybe_unused]] auto const c1 = GetImplementation().fRho.SecondDerivative(
line.GetPosition(0), line.GetDirection(0));
// todo: provide a real, working implementation
return 1_m * std::numeric_limits<double>::infinity();
}
};
} // namespace corsika::environment
Environment/MediumProperties.h deleted 100644 → 0
View file @ dbf6e995
This diff is collapsed.
Environment/MediumTypes.h deleted 100644 → 0
View file @ dbf6e995
/*
* (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.
*/
#pragma once
namespace corsika::environment {
/**
* Medium types are useful most importantly for effective models
* like energy losses. a particular medium (mixture of components)
* may have specif properties not reflected by its mixture of
* components.
*/
enum class EMediumType { eUnknown, eAir, eWater, eIce, eRock };
} // namespace corsika::environment

data privacy - impressum