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 302 additions and 676 deletions
examples/hybrid_MC.cpp → examples/cascade_examples/mc_conex.cpp
View file @ 0c1d3948
......@@ -6,50 +6,43 @@
* 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/framework/process/InteractionCounter.hpp>
/* clang-format on */
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/Step.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
// #include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/modules/writers/PrimaryWriter.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/UrQMD.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/CONEX.hpp>
#include <corsika/setup/SetupStack.hpp>
......@@ -63,6 +56,12 @@
using namespace corsika;
using namespace std;
//
// An example of running an EAS where the hadronic cascade is
// handled by sibyll+URQMD and the EM cascade is treated with
// CONEX + Bethe Bloch (as opposed to PROPOSAL).
//
/**
* Random number stream initialization
*
......@@ -70,17 +69,19 @@ using namespace std;
*/
void registerRandomStreams(uint64_t seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("qgsjet");
RNGManager<>::getInstance().registerRandomStream("sibyll");
RNGManager<>::getInstance().registerRandomStream("conex");
RNGManager<>::getInstance().registerRandomStream("epos");
RNGManager<>::getInstance().registerRandomStream("proposal");
RNGManager<>::getInstance().registerRandomStream("pythia");
RNGManager<>::getInstance().registerRandomStream("qgsjet");
RNGManager<>::getInstance().registerRandomStream("sibyll");
RNGManager<>::getInstance().registerRandomStream("urqmd");
RNGManager<>::getInstance().registerRandomStream("proposal");
RNGManager<>::getInstance().registerRandomStream("conex");
if (seed == 0) {
std::random_device rd;
seed = rd();
CORSIKA_LOG_INFO("new random seed (auto) {}", seed);
CORSIKA_LOG_INFO("random seed (auto) {} ", seed);
} else {
CORSIKA_LOG_INFO("random seed {} ", seed);
}
RNGManager<>::getInstance().setSeed(seed);
}
......@@ -141,8 +142,12 @@ private:
/**
* Selection of environment interface implementation:
*/
using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
using EnvType = Environment<EnvironmentInterface>;
template <typename T>
using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
using StackType = setup::Stack<EnvType>;
using TrackingType = setup::Tracking;
int main(int argc, char** argv) {
......@@ -157,7 +162,6 @@ int main(int argc, char** argv) {
" if no seed is given, a random seed is chosen");
return 1;
}
feenableexcept(FE_INVALID);
uint64_t seed = 0;
if (argc > 4) seed = std::stol(std::string(argv[4]));
......@@ -165,20 +169,15 @@ int main(int argc, char** argv) {
registerRandomStreams(seed);
// setup environment, geometry
using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
using EnvType = Environment<EnvironmentInterface>;
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
// build a Linsley US Standard atmosphere into `env`
MagneticFieldVector bField{rootCS, 50_uT, 0_T, 0_T};
create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm,
MagneticFieldVector{rootCS, 0_T, 50_uT, 0_T});
env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm, bField);
// setup particle stack, and add primary particle
setup::Stack<EnvType> stack;
stack.clear();
unsigned short const A = std::stoi(std::string(argv[1]));
unsigned short const Z = std::stoi(std::string(argv[2]));
Code const beamCode = get_nucleus_code(A, Z);
......@@ -194,12 +193,6 @@ int main(int argc, char** argv) {
auto const [px, py, pz] = momentumComponents(thetaRad, P0);
auto plab = MomentumVector(rootCS, {px, py, pz});
CORSIKA_LOG_INFO(
"input particle: {}, "
"input angles: theta={}, "
"input momentum: {} GeV, "
", norm={}",
beamCode, theta, plab.getComponents() / 1_GeV, plab.getNorm());
auto const observationHeight = 0_km + constants::EarthRadius::Mean;
auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
......@@ -211,28 +204,25 @@ int main(int argc, char** argv) {
showerCore +
Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
CORSIKA_LOG_INFO("point of injection: {} ", injectionPos.getCoordinates());
stack.addParticle(std::make_tuple(
Code::Proton, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
plab.normalized(), injectionPos, 0_ns));
CORSIKA_LOG_INFO("shower axis length: {} m",
(showerCore - injectionPos).getNorm() * 1.02);
OutputManager output("hybrid_MC_outputs");
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000};
auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
CORSIKA_LOG_INFO("Primary particle: {}", beamCode);
CORSIKA_LOG_INFO("Zenith angle: {} (rad)", theta);
CORSIKA_LOG_INFO("Momentum: {} (GeV)", plab.getComponents() / 1_GeV);
CORSIKA_LOG_INFO("Propagation dir: {}", plab.getNorm());
CORSIKA_LOG_INFO("Injection point: {}", injectionPos.getCoordinates());
CORSIKA_LOG_INFO("shower axis length: {} ",
(showerCore - injectionPos).getNorm() * 1.02);
// setup processes, decays and interactions
// SETUP WRITERS
corsika::sibyll::Interaction sibyll{env};
InteractionCounter sibyllCounted{sibyll};
corsika::pythia8::Decay decayPythia;
OutputManager output("hybrid_MC_outputs");
// register energy losses as output
EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
output.add("energyloss", dEdX);
// create a track writer and register it with the output manager
......@@ -242,19 +232,29 @@ int main(int argc, char** argv) {
ParticleCut<SubWriter<decltype(dEdX)>> cut(3_GeV, false, dEdX);
BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss(dEdX);
LongitudinalWriter profile{showerAxis, 10_g / square(1_cm)};
LongitudinalWriter profile{showerAxis, nAxisBins, dX};
output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane<TrackingType> observationLevel(obsPlane,
DirectionVector(rootCS, {1., 0., 0.}));
output.add("particles", observationLevel);
PrimaryWriter<TrackingType, ParticleWriterParquet> primaryWriter(observationLevel);
output.add("primary", primaryWriter);
// SETUP PROCESSES, DECAYS, INTERACTIONS
corsika::sibyll::Interaction sibyll{env};
InteractionCounter sibyllCounted{sibyll};
corsika::pythia8::Decay decayPythia;
CONEXhybrid // SubWriter<decltype(dEdX>, SubWriter<decltype(profile)>>
conex_model(center, showerAxis, t, injectionHeight, E0, get_PDG(Code::Proton), dEdX,
profile);
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane<setup::Tracking> observationLevel(
obsPlane, DirectionVector(rootCS, {1., 0., 0.}), "particles.dat");
output.add("obsplane", observationLevel);
corsika::urqmd::UrQMD urqmd_model;
InteractionCounter urqmdCounted{urqmd_model};
......@@ -265,7 +265,7 @@ int main(int argc, char** argv) {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
bool operator()(const setup::Stack<EnvType>::particle_type& p) const {
bool operator()(const StackType::particle_type& p) const {
return (p.getEnergy() < cutE_);
}
};
......@@ -273,17 +273,27 @@ int main(int argc, char** argv) {
auto sequence = make_sequence(hadronSequence, decayPythia, eLoss, cut, conex_model,
longprof, observationLevel, trackCheck);
StackType stack;
stack.clear();
// define air shower object, run simulation
setup::Tracking tracking;
TrackingType tracking;
Cascade EAS(env, tracking, sequence, output, stack);
output.startOfLibrary();
auto const primaryProperties = std::make_tuple(
Code::Proton, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
plab.normalized(), injectionPos, 0_ns);
stack.addParticle(primaryProperties);
primaryWriter.recordPrimary(primaryProperties);
// to fix the point of first interaction, uncomment the following two lines:
// EAS.SetNodes();
// EAS.forceInteraction();
output.startOfLibrary();
EAS.run();
output.endOfLibrary();
const HEPEnergyType Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
CORSIKA_LOG_INFO(
......@@ -296,5 +306,7 @@ int main(int argc, char** argv) {
save_hist(hists.labHist(), "inthist_lab_hybrid.npz", true);
save_hist(hists.CMSHist(), "inthist_cms_hybrid.npz", true);
output.endOfLibrary();
CORSIKA_LOG_INFO("done");
}
examples/radio_em_shower.cpp → examples/cascade_examples/radio_em_shower.cpp
View file @ 0c1d3948
......@@ -21,10 +21,11 @@
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/modules/writers/PrimaryWriter.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
......@@ -64,6 +65,13 @@
using namespace corsika;
using namespace std;
//
// An example of running an casacde which generates radio input for a
// cascade that has hadronic interactions turned off. Currently, this
// will produce output for only one antenna, but there are lines below,
// which are commented out, to enable a full star-shape pattern of antennas
//
void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("proposal");
......@@ -72,14 +80,21 @@ void registerRandomStreams(int seed) {
if (seed == 0) {
std::random_device rd;
seed = rd();
cout << "new random seed (auto) " << seed << endl;
CORSIKA_LOG_INFO("random seed (auto) {} ", seed);
} else {
CORSIKA_LOG_INFO("random seed {} ", seed);
}
RNGManager<>::getInstance().setSeed(seed);
}
using EnvironmentInterface =
IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
using EnvType = Environment<EnvironmentInterface>;
template <typename TInterface>
using MyExtraEnv =
UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
using StackType = setup::Stack<EnvType>;
using TrackingType = setup::Tracking;
int main(int argc, char** argv) {
......@@ -93,22 +108,21 @@ int main(int argc, char** argv) {
}
int seed{static_cast<int>(std::stof(std::string(argv[2])))};
std::cout << "Seed: " << seed << std::endl;
CORSIKA_LOG_INFO("Seed {}", seed);
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
registerRandomStreams(seed);
// setup environment, geometry
using EnvironmentInterface =
IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
using EnvType = Environment<EnvironmentInterface>;
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
double const refractive_index = 1.000327;
MagneticFieldVector const bField{rootCS, 50_uT, 0_T, 0_T};
create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
env, AtmosphereId::LinsleyUSStd, center, 1.000327, Medium::AirDry1Atm,
MagneticFieldVector{rootCS, 50_uT, 0_T, 0_T});
env, AtmosphereId::LinsleyUSStd, center, refractive_index, Medium::AirDry1Atm,
bField);
std::unordered_map<Code, HEPEnergyType> energy_resolution = {
{Code::Electron, 5_MeV},
......@@ -118,26 +132,19 @@ int main(int argc, char** argv) {
for (auto [pcode, energy] : energy_resolution)
set_energy_production_threshold(pcode, energy);
// setup particle stack, and add primary particle
setup::Stack<EnvType> stack;
stack.clear();
const Code beamCode = Code::Electron;
Code const beamCode = Code::Electron;
auto const mass = get_mass(beamCode);
const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
double theta = 0.;
HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[1]));
double const theta = 0.;
auto const thetaRad = theta / 180. * M_PI;
HEPMomentumType P0 = calculate_momentum(E0, mass);
HEPMomentumType const P0 = calculate_momentum(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 = 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.getNorm() << endl;
auto const plab = MomentumVector(rootCS, {px, py, pz});
auto const observationHeight = 1.4_km + constants::EarthRadius::Mean;
auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
......@@ -148,104 +155,103 @@ int main(int argc, char** argv) {
Point const injectionPos =
showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
stack.addParticle(std::make_tuple(
beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
plab.normalized(), injectionPos, 0_ns));
CORSIKA_LOG_INFO("shower axis length: {} ",
(showerCore - injectionPos).getNorm() * 1.02);
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000};
auto const dX = 10_g / square(1_cm); // Binning of the writers along the shower axis
uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
// setup the radio antennas
TimeType const groundHitTime{(showerCore - injectionPos).getNorm() / constants::c};
std::string outname_ = "radio_em_shower_outputs"; // + std::to_string(rr_);
OutputManager output(outname_);
// Radio antennas and relevant information
// the antenna time variables
const TimeType duration_{1e-6_s};
const InverseTimeType sampleRate_{1e+9_Hz};
TimeType const duration{1e-6_s};
InverseTimeType const sampleRate{1e+9_Hz};
// the detector (aka antenna collection) for CoREAS and ZHS
AntennaCollection<TimeDomainAntenna> detectorCoREAS;
AntennaCollection<TimeDomainAntenna> detectorZHS;
auto const showerCoreX_{showerCore.getCoordinates().getX()};
auto const showerCoreY_{showerCore.getCoordinates().getY()};
auto const injectionPosX_{injectionPos.getCoordinates().getX()};
auto const injectionPosY_{injectionPos.getCoordinates().getY()};
auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
auto const showerCoreX{showerCore.getCoordinates().getX()};
auto const showerCoreY{showerCore.getCoordinates().getY()};
auto const injectionPosX{injectionPos.getCoordinates().getX()};
auto const injectionPosY{injectionPos.getCoordinates().getY()};
auto const injectionPosZ{injectionPos.getCoordinates().getZ()};
auto const triggerpoint{Point(rootCS, injectionPosX, injectionPosY, injectionPosZ)};
CORSIKA_LOG_INFO("Primary particle: {}", beamCode);
CORSIKA_LOG_INFO("Zenith angle: {} (rad)", theta);
CORSIKA_LOG_INFO("Momentum: {} (GeV)", plab.getComponents() / 1_GeV);
CORSIKA_LOG_INFO("Propagation dir: {}", plab.getNorm());
CORSIKA_LOG_INFO("Injection point: {}", injectionPos.getCoordinates());
CORSIKA_LOG_INFO("shower axis length: {} ",
(showerCore - injectionPos).getNorm() * 1.02);
CORSIKA_LOG_INFO("Trigger Point is: {}", triggerpoint);
// // setup CoREAS antennas - use the for loop for star shape pattern
// for (auto radius_1 = 25_m; radius_1 <= 500_m; radius_1 += 25_m) {
// for (auto phi_1 = 0; phi_1 <= 315; phi_1 += 45) {
auto radius_1 = 200_m;
auto phi_1 = 45;
auto phiRad_1 = phi_1 / 180. * M_PI;
auto rr_1 = static_cast<int>(radius_1 / 1_m);
auto const point_1{Point(rootCS, showerCoreX_ + radius_1 * cos(phiRad_1),
showerCoreY_ + radius_1 * sin(phiRad_1),
constants::EarthRadius::Mean)};
std::cout << "Antenna point: " << point_1 << std::endl;
auto triggertime_1{(triggerpoint_ - point_1).getNorm() / constants::c};
std::string name_1 =
"CoREAS_R=" + std::to_string(rr_1) + "_m--Phi=" + std::to_string(phi_1) + "degrees";
TimeDomainAntenna antenna_1(name_1, point_1, rootCS, triggertime_1, duration_,
sampleRate_, triggertime_1);
detectorCoREAS.addAntenna(antenna_1);
// for (auto radius_coreas = 25_m; radius_coreas <= 500_m; radius_coreas += 25_m) {
// for (auto phi_coreas = 0; phi_coreas <= 315; phi_coreas += 45) {
auto radius_coreas = 200_m;
auto phi_coreas = 45;
auto phiRad_coreas = phi_coreas / 180. * M_PI;
auto rr_coreas = static_cast<int>(radius_coreas / 1_m);
auto const point_coreas{Point(rootCS, showerCoreX + radius_coreas * cos(phiRad_coreas),
showerCoreY + radius_coreas * sin(phiRad_coreas),
constants::EarthRadius::Mean)};
std::cout << "Antenna point: " << point_coreas << std::endl;
CORSIKA_LOG_INFO("Antenna point {}", injectionPos.getCoordinates());
auto triggertime_coreas{(triggerpoint - point_coreas).getNorm() / constants::c};
std::string name_coreas = "CoREAS_R=" + std::to_string(rr_coreas) +
"_m--Phi=" + std::to_string(phi_coreas) + "degrees";
TimeDomainAntenna antenna_coreas(name_coreas, point_coreas, rootCS, triggertime_coreas,
duration, sampleRate, triggertime_coreas);
detectorCoREAS.addAntenna(antenna_coreas);
// }
// }
// // setup ZHS antennas - use the for loop for star shape pattern
// for (auto radius_ = 25_m; radius_ <= 500_m; radius_ += 25_m) {
// for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
auto radius_ = 200_m;
auto phi_ = 45;
auto phiRad_ = phi_ / 180. * M_PI;
auto rr_ = static_cast<int>(radius_ / 1_m);
auto const point_{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_),
showerCoreY_ + radius_ * sin(phiRad_),
constants::EarthRadius::Mean)};
auto triggertime_{(triggerpoint_ - point_).getNorm() / constants::c};
std::string name_ =
"ZHS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
TimeDomainAntenna antenna_(name_, point_, rootCS, triggertime_, duration_, sampleRate_,
triggertime_);
detectorZHS.addAntenna(antenna_);
// for (auto radius_zhs = 25_m; radius_zhs <= 500_m; radius_zhs += 25_m) {
// for (auto phi_zhs = 0; phi_zhs <= 315; phi_zhs += 45) {
auto radius_zhs = 200_m;
auto phi_zhs = 45;
auto phiRad_zhs = phi_zhs / 180. * M_PI;
auto rr_zhs = static_cast<int>(radius_zhs / 1_m);
auto const point_zhs{Point(rootCS, showerCoreX + radius_zhs * cos(phiRad_zhs),
showerCoreY + radius_zhs * sin(phiRad_zhs),
constants::EarthRadius::Mean)};
auto triggertime_zhs{(triggerpoint - point_zhs).getNorm() / constants::c};
std::string name_zhs = "ZHS_R=" + std::to_string(rr_zhs) +
"_m--Phi=" + std::to_string(phi_zhs) + "degrees";
TimeDomainAntenna antenna_zhs(name_zhs, point_zhs, rootCS, triggertime_zhs, duration,
sampleRate, triggertime_zhs);
detectorZHS.addAntenna(antenna_zhs);
// }
// }
// setup processes, decays and interactions
EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
// register energy losses as output
output.add("dEdX", dEdX);
ParticleCut<SubWriter<decltype(dEdX)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true, dEdX);
corsika::sophia::InteractionModel sophia;
EnergyLossWriter energyloss{showerAxis, dX, nAxisBins};
ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true,
energyloss);
corsika::sibyll::Interaction sibyll{env};
corsika::sophia::InteractionModel sophia;
HEPEnergyType heThresholdNN = 80_GeV;
corsika::proposal::Interaction emCascade(
env, sophia, sibyll.getHadronInteractionModel(), heThresholdNN);
corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
// BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
corsika::proposal::ContinuousProcess<SubWriter<decltype(energyloss)>> emContinuous(
env, energyloss);
// NOT possible right now, due to interface differenc in PROPOSAL
// InteractionCounter emCascadeCounted(emCascade);
OutputManager output("radio_em_shower_outputs");
output.add("energyloss", energyloss);
TrackWriter tracks;
output.add("tracks", tracks);
// long. profile
LongitudinalWriter profile{showerAxis, 10_g / square(1_cm)};
LongitudinalWriter profile{showerAxis, nAxisBins, dX};
output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
......@@ -269,15 +275,27 @@ int main(int argc, char** argv) {
output.add("ZHS", zhs);
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
ObservationPlane<TrackingType, ParticleWriterParquet> observationLevel{
obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
output.add("particles", observationLevel);
// auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, coreas, zhs);
PrimaryWriter<TrackingType, ParticleWriterParquet> primaryWriter(observationLevel);
output.add("primary", primaryWriter);
auto sequence = make_sequence(emCascade, emContinuous, longprof, coreas, zhs, tracks,
observationLevel, cut);
// define air shower object, run simulation
setup::Tracking tracking;
TrackingType tracking;
auto const primaryProperties = std::make_tuple(
beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
plab.normalized(), injectionPos, 0_ns);
// setup particle stack, and add primary particle
StackType stack;
stack.clear();
stack.addParticle(primaryProperties);
primaryWriter.recordPrimary(primaryProperties);
output.startOfLibrary();
Cascade EAS(env, tracking, sequence, output, stack);
......@@ -287,7 +305,8 @@ int main(int argc, char** argv) {
EAS.run();
HEPEnergyType const Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
HEPEnergyType const Efinal =
energyloss.getEnergyLost() + observationLevel.getEnergyGround();
CORSIKA_LOG_INFO(
"total energy budget (GeV): {}, "
......@@ -295,4 +314,6 @@ int main(int argc, char** argv) {
Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
output.endOfLibrary();
return EXIT_SUCCESS;
}
examples/water.cpp → examples/cascade_examples/water.cpp
View file @ 0c1d3948
......@@ -24,15 +24,16 @@
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/modules/writers/PrimaryWriter.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/Sophia.hpp>
#include <corsika/modules/UrQMD.hpp>
#include <corsika/modules/FLUKA.hpp>
#include <corsika/modules/tracking/TrackingStraight.hpp>
#include <corsika/output/OutputManager.hpp>
......@@ -41,6 +42,7 @@
#include <corsika/stack/VectorStack.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
#include <CLI/App.hpp>
#include <CLI/Config.hpp>
......@@ -50,9 +52,9 @@ using namespace corsika;
using IMediumType = IMediumPropertyModel<IMediumModel>;
using EnvType = Environment<IMediumType>;
using StackActive = setup::Stack<EnvType>;
using StackView = StackActive::stack_view_type;
using Particle = StackActive::particle_type;
using StackType = setup::Stack<EnvType>;
using TrackingType = tracking_line::Tracking;
using Particle = StackType::particle_type;
void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
......@@ -61,7 +63,7 @@ void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("sophia");
RNGManager<>::getInstance().registerRandomStream("epos");
RNGManager<>::getInstance().registerRandomStream("pythia");
RNGManager<>::getInstance().registerRandomStream("urqmd");
RNGManager<>::getInstance().registerRandomStream("fluka");
RNGManager<>::getInstance().registerRandomStream("proposal");
if (seed == 0) {
std::random_device rd;
......@@ -75,12 +77,12 @@ void registerRandomStreams(int seed) {
int main(int argc, char** argv) {
// * process input
Code primaryType;
HEPEnergyType e0, eCut;
Code beamCode;
HEPEnergyType E0, eCut;
int A, Z, n_event;
int randomSeed;
std::string output_dir;
CLI::App app{"Neutrino event generator"};
CLI::App app{"Cascade in water"};
// we start by definining a sub-group for the primary ID
auto opt_Z = app.add_option("-Z", Z, "Atomic number for primary")
->check(CLI::Range(0, 26))
......@@ -110,8 +112,25 @@ int main(int argc, char** argv) {
app.add_option("-s", randomSeed, "Seed for random number")
->check(CLI::NonNegativeNumber)
->default_val(0);
// parse the command line options into the variables
CLI11_PARSE(app, argc, argv);
std::string_view const loglevel = app["-v"]->as<std::string_view>();
if (loglevel == "warn") {
logging::set_level(logging::level::warn);
} else if (loglevel == "info") {
logging::set_level(logging::level::info);
} else if (loglevel == "debug") {
logging::set_level(logging::level::debug);
} else if (loglevel == "trace") {
#ifndef _C8_DEBUG_
CORSIKA_LOG_ERROR("trace log level requires a Debug build.");
return 1;
#endif
logging::set_level(logging::level::trace);
}
// check that we got either PDG or A/Z
// this can be done with option_groups but the ordering
// gets all messed up
......@@ -121,37 +140,24 @@ int main(int argc, char** argv) {
return 1;
}
}
// initialize random number sequence(s)
registerRandomStreams(randomSeed);
// check if we want to use a PDG code instead
if (app.count("--pdg") > 0) {
primaryType = convert_from_PDG(PDGCode(app["--pdg"]->as<int>()));
beamCode = convert_from_PDG(PDGCode(app["--pdg"]->as<int>()));
} else {
// check manually for proton and neutrons
if ((A == 1) && (Z == 1))
primaryType = Code::Proton;
beamCode = Code::Proton;
else if ((A == 1) && (Z == 0))
primaryType = Code::Neutron;
beamCode = Code::Neutron;
else
primaryType = get_nucleus_code(A, Z);
beamCode = get_nucleus_code(A, Z);
}
e0 = app["-E"]->as<double>() * 1_GeV;
eCut = app["--eCut"]->as<double>() * 1_GeV;
std::string_view const loglevel = app["-v"]->as<std::string_view>();
if (loglevel == "warn") {
logging::set_level(logging::level::warn);
} else if (loglevel == "info") {
logging::set_level(logging::level::info);
} else if (loglevel == "debug") {
logging::set_level(logging::level::debug);
} else if (loglevel == "trace") {
#ifndef _C8_DEBUG_
CORSIKA_LOG_ERROR("trace log level requires a Debug build.");
return 1;
#endif
logging::set_level(logging::level::trace);
}
registerRandomStreams(randomSeed);
// * environment and universe
EnvType env;
......@@ -163,27 +169,26 @@ int main(int argc, char** argv) {
Point const center{rootCS, 0_m, 0_m, 0_m};
auto sphere = std::make_unique<Sphere>(center, 100_m);
auto node = std::make_unique<VolumeTreeNode<IMediumType>>(std::move(sphere));
// Hydrogen is not supported by UrQMD yet. See #456
auto comp = NuclearComposition({{Code::Oxygen}, {1.0}});
NuclearComposition const nuclearComposition{{Code::Hydrogen, Code::Oxygen},
{2.0 / 3.0, 1.0 / 3.0}};
// density of sea water
auto density = 1.02_g / (1_cm * 1_cm * 1_cm);
auto water_medium =
std::make_shared<MediumPropertyModel<HomogeneousMedium<IMediumType>>>(
Medium::WaterLiquid, density, comp);
Medium::WaterLiquid, density, nuclearComposition);
node->setModelProperties(water_medium);
universe->addChild(std::move(node));
}
// * detector geometry
// * make downward-going shower axis and a observation plane in x-y-plane
auto injectorLength = 50_m;
Point const injectorPos = Point(rootCS, {0_m, 0_m, injectorLength});
auto const& injectCS = make_translation(rootCS, injectorPos.getCoordinates());
Point const injectionPos = Point(rootCS, {0_m, 0_m, injectorLength});
auto const& injectCS = make_translation(rootCS, injectionPos.getCoordinates());
DirectionVector upVec(rootCS, {0., 0., 1.});
DirectionVector leftVec(rootCS, {1., 0., 0.});
DirectionVector downVec(rootCS, {0., 0., -1.});
// * observation plane
std::vector<ObservationPlane<tracking_line::Tracking>> obsPlanes;
std::vector<ObservationPlane<TrackingType, ParticleWriterParquet>> obsPlanes;
const int nPlane = 5;
for (int i = 0; i < nPlane - 1; i++) {
Point planeCenter{injectCS, {0_m, 0_m, -(i + 1) * 3_m}};
......@@ -193,12 +198,14 @@ int main(int argc, char** argv) {
Plane{Point{injectCS, {0_m, 0_m, -50_m}}, upVec}, leftVec, true);
// * longitutional profile
ShowerAxis const showerAxis{injectorPos, 1.2 * injectorLength * downVec, env};
LongitudinalWriter longiWriter{showerAxis, 5500, 1_g / square(1_cm)};
ShowerAxis const showerAxis{injectionPos, 1.2 * injectorLength * downVec, env};
auto const dX = 1_g / square(1_cm); // Binning of the writers along the shower axis
uint const nAxisBins = showerAxis.getMaximumX() / dX + 1; // Get maximum number of bins
LongitudinalWriter longiWriter{showerAxis, nAxisBins, dX};
LongitudinalProfile<SubWriter<decltype(longiWriter)>> longprof{longiWriter};
// * energy loss profile
EnergyLossWriter dEdX{showerAxis, 1_g / square(1_cm), 5500};
EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
// * physical process list
// particle production threshold
......@@ -207,19 +214,19 @@ int main(int argc, char** argv) {
ParticleCut<SubWriter<decltype(dEdX)>> cut(emCut, emCut, hadCut, hadCut, true, dEdX);
// hadronic interactions
HEPEnergyType heHadronModelThreshold = 63.1_GeV;
HEPEnergyType heHadronModelThreshold = std::pow(10, 1.9) * 1_GeV;
corsika::sibyll::Interaction sibyll(env);
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted(urqmd);
corsika::fluka::Interaction leIntModel{env};
InteractionCounter leIntCounted{leIntModel};
struct EnergySwitch {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
bool operator()(const Particle& p) const { return (p.getKineticEnergy() < cutE_); }
};
// auto lowModel = make_sequence(urqmd);
auto hadronSequence =
make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, sibyll);
make_select(EnergySwitch(heHadronModelThreshold), leIntCounted, sibyll);
// decay process
corsika::pythia8::Decay decayPythia;
......@@ -243,19 +250,41 @@ int main(int argc, char** argv) {
// hard coded
auto obsPlaneSequence =
make_sequence(obsPlanes[0], obsPlanes[1], obsPlanes[2], obsPlanes[3], obsPlanes[4]);
output.add("longi_profile", longiWriter);
output.add("energy_loss", dEdX);
PrimaryWriter<TrackingType, ParticleWriterParquet> primaryWriter(obsPlanes.back());
output.add("primary", primaryWriter);
output.add("profile", longiWriter);
output.add("energyloss", dEdX);
// * the final process sequence
auto sequence = make_sequence(physics_sequence, longprof, obsPlaneSequence, cut);
// * tracking and stack
tracking_line::Tracking tracking;
StackActive stack;
TrackingType tracking;
StackType stack;
// * cascade manager
Cascade EAS(env, tracking, sequence, output, stack);
E0 = app["-E"]->as<double>() * 1_GeV;
HEPEnergyType mass = get_mass(beamCode);
// convert Elab to Plab
HEPMomentumType P0 = calculate_momentum(E0, mass);
auto plab = MomentumVector(rootCS, P0 * downVec.getNorm());
// print our primary parameters all in one place
if (app["--pdg"]->count() > 0) {
CORSIKA_LOG_INFO("Primary PDG ID: {}", app["--pdg"]->as<int>());
} else {
CORSIKA_LOG_INFO("Primary Z/A: {}/{}", Z, A);
}
CORSIKA_LOG_INFO("Primary Energy: {}", E0);
CORSIKA_LOG_INFO("Primary Momentum: {}", P0);
CORSIKA_LOG_INFO("Primary Direction: {}", plab.getNorm());
CORSIKA_LOG_INFO("Point of Injection: {}", injectionPos.getCoordinates());
CORSIKA_LOG_INFO("Shower Axis Length: {}", injectorLength);
// * main loop
output.startOfLibrary();
for (int i_shower = 0; i_shower < n_event; i_shower++) {
......@@ -263,8 +292,11 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO("Event: {} / {}", i_shower, n_event);
// * inject primary
auto primary = stack.addParticle(std::make_tuple(
primaryType, e0 - get_mass(primaryType), downVec, injectorPos, 0_ns));
auto const primaryProperties = std::make_tuple(
beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
plab.normalized(), injectionPos, 0_ns);
auto primary = stack.addParticle(primaryProperties);
stack.addParticle(primaryProperties);
EAS.run();
......@@ -273,8 +305,8 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO(
"total energy budget (TeV): {:.2f} (dEdX={:.2f} ground={:.2f}), "
"relative difference (%): {:.3f}",
e0 / 1_TeV, dEdX.getEnergyLost() / 1_TeV, obsPlaneFinal.getEnergyGround() / 1_TeV,
(Efinal / e0 - 1.) * 100.);
E0 / 1_TeV, dEdX.getEnergyLost() / 1_TeV, obsPlaneFinal.getEnergyGround() / 1_TeV,
(Efinal / E0 - 1.) * 100.);
}
output.endOfLibrary();
......
examples/cascade_proton_example.cpp deleted 100644 → 0
View file @ db6fe545
/*
* (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/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/HadronicElasticModel.hpp>
#include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <iostream>
#include <limits>
#include <typeinfo>
using namespace corsika;
using namespace std;
//
// The example main program for a particle cascade
//
int main() {
logging::set_level(logging::level::warn);
std::cout << "cascade_proton_example" << std::endl;
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
RNGManager<>::getInstance().registerRandomStream("cascade");
OutputManager output("cascade_proton_outputs");
// setup environment, geometry
using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
using EnvType = Environment<EnvironmentInterface>;
EnvType env;
auto& universe = *(env.getUniverse());
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
auto world = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
using MyHomogeneousModel =
MediumPropertyModel<UniformMagneticField<HomogeneousMedium<EnvironmentInterface>>>;
world->setModelProperties<MyHomogeneousModel>(
Medium::AirDry1Atm, MagneticFieldVector(rootCS, 0_T, 0_T, 1_mT),
1_kg / (1_m * 1_m * 1_m), NuclearComposition({Code::Proton}, {1.}));
universe.addChild(std::move(world));
// setup particle stack, and add primary particle
setup::Stack<EnvType> stack;
stack.clear();
const Code beamCode = Code::Proton;
const HEPMassType mass = Proton::mass;
const HEPEnergyType E0 = 200_GeV;
double theta = 0.;
double phi = 0.;
Point injectionPos(rootCS, 0_m, 0_m, 0_m);
{
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 = 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::make_tuple(
beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
plab.normalized(), injectionPos, 0_ns));
}
// setup processes, decays and interactions
setup::Tracking tracking;
StackInspector<setup::Stack<EnvType>> stackInspect(1000, true, E0);
RNGManager<>::getInstance().registerRandomStream("pythia");
RNGManager<>::getInstance().registerRandomStream("sibyll");
corsika::sibyll::Interaction sibyll{env};
// corsika::pythia8::Interaction pythia;
corsika::pythia8::Decay decay;
ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
EnergyLossWriter dEdX{showerAxis};
output.add("energyloss", dEdX);
BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss{dEdX};
ParticleCut<SubWriter<decltype(dEdX)>> cut(60_GeV, true, dEdX);
cut.printThresholds();
// RNGManager::getInstance().registerRandomStream("HadronicElasticModel");
// HadronicElasticModel::HadronicElasticInteraction
// hadronicElastic(env);
TrackWriter trackWriter;
output.add("tracks", trackWriter); // register TrackWriter
// assemble all processes into an ordered process list
auto sequence = make_sequence(sibyll, decay, eLoss, trackWriter, stackInspect, cut);
// define air shower object, run simulation
Cascade EAS(env, tracking, sequence, output, stack);
output.startOfLibrary();
EAS.run();
output.endOfLibrary();
CORSIKA_LOG_INFO("Done");
}
examples/boundary_example.cpp → examples/framework_examples/boundary_crossing.cpp
View file @ 0c1d3948
......@@ -37,6 +37,13 @@
using namespace corsika;
using namespace std;
//
// This example shows how to make a custom process which deletes particles that
// cross a particular boundary (a sphere in this case)
// For a plane boundary, this can be implemented by adding an ObservationPlane
// or Observation Volume object instead (the standard "absorbing" geometry objects)
//
template <bool deleteParticle>
struct MyBoundaryCrossingProcess
: public BoundaryCrossingProcess<MyBoundaryCrossingProcess<deleteParticle>> {
......@@ -65,9 +72,6 @@ private:
std::ofstream file_;
};
//
// The example main program for a particle cascade
//
int main() {
logging::set_level(logging::level::warn);
......
examples/particle_list_example.cpp → examples/framework_examples/known_particles.cpp
View file @ 0c1d3948
......@@ -18,11 +18,14 @@ using namespace corsika;
using namespace std;
//
// The example main program for a particle list
// This example prints out all of the particles that are
// known to CORSIKA 8 and the corresponding IDs in the
// hadronic model codes.
//
int main() {
logging::set_level(logging::level::warn);
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] %v");
logging::info(
......@@ -30,7 +33,7 @@ int main() {
"------------------------------------------\n"
" particles in CORSIKA\n"
"------------------------------------------\n");
int const width = 20 + 10 + 10 + 10 + 15 + 15 + 17;
int const width = 20 + 10 + 10 + 10 + 16 + 16 + 17;
logging::info(
"Name | "
"PDG-id | "
......@@ -46,7 +49,7 @@ int main() {
? to_string(corsika::sibyll::getSibyllMass(p) / 1_GeV)
: "");
auto const qgs_id = corsika::qgsjetII::convertToQgsjetII(p);
logging::info("{:20} | {:10} | {:10} | {:10} | {:>15.5} | {:>15.5} |", p,
logging::info("{:20} | {:10} | {:10} | {:10} | {:>16.5} | {:>16.5} |", p,
static_cast<int>(get_PDG(p)),
(sib_id != corsika::sibyll::SibyllCode::Unknown
? to_string(static_cast<int>(sib_id))
......
examples/staticsequence_example.cpp → examples/framework_examples/static_sequence.cpp
View file @ 0c1d3948
......@@ -131,4 +131,4 @@ int main() {
modular();
return 0;
}
\ No newline at end of file
}
examples/stopping_power.cpp → examples/physics_examples/stopping_power.cpp
View file @ 0c1d3948
......@@ -9,7 +9,6 @@
#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMediumModel.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
......@@ -30,7 +29,7 @@ using namespace std;
//
int main() {
logging::set_level(logging::level::warn);
logging::set_level(logging::level::info);
CORSIKA_LOG_INFO("stopping_power");
......@@ -52,7 +51,9 @@ int main() {
setup::Stack<EnvType> stack;
std::ofstream file("dEdX.dat");
std::string fileName = "dEdX.dat";
CORSIKA_LOG_INFO("Writing to file {}", fileName);
std::ofstream file(fileName);
file << "# beta*gamma, dE/dX / MeV/(g/cm²)" << std::endl;
for (HEPEnergyType E0 = 200_MeV; E0 < 1_PeV; E0 *= 1.05) {
......@@ -79,5 +80,5 @@ int main() {
HEPEnergyType dE = eLoss.getTotalEnergyLoss(p, 1_g / square(1_cm));
file << P0 / mass << "\t" << -dE / 1_MeV << std::endl;
}
CORSIKA_LOG_INFO("finished writing dEdX.dat");
CORSIKA_LOG_INFO("finished writing {}", fileName);
}
examples/clover_leaf.cpp → examples/physics_examples/synchrotron_clover_leaf.cpp
View file @ 0c1d3948
......@@ -24,7 +24,6 @@
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/UniformRefractiveIndex.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
......@@ -35,14 +34,10 @@
#include <corsika/modules/radio/antennas/Antenna.hpp>
#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
#include <corsika/modules/radio/detectors/AntennaCollection.hpp>
#include <corsika/modules/radio/propagators/NumericalIntegratingPropagator.hpp>
#include <corsika/modules/radio/propagators/DummyTestPropagator.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/writers/PrimaryWriter.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/TimeCut.hpp>
// #include <corsika/modules/TrackWriter.hpp>
#include <iomanip>
#include <iostream>
......@@ -57,6 +52,7 @@ using namespace std;
// A simple shower to get the electric field trace of an electron (& a positron)
// in order to reveal the "clover leaf" pattern of energy fluence to the ground.
//
int main() {
logging::set_level(logging::level::info);
......@@ -70,7 +66,7 @@ int main() {
auto seed = rd();
RNGManager<>::getInstance().setSeed(seed);
OutputManager output("1 electron - 1 positron");
OutputManager output("clover_leaf_outputs");
// create a suitable environment
using IModelInterface =
......@@ -105,7 +101,6 @@ int main() {
auto const injectionPosY_{injectionPos.getCoordinates().getY()};
auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
// the antenna characteristics
const TimeType duration_{2e-6_s}; // 0.8e-4_s
......@@ -113,7 +108,6 @@ int main() {
// the detectors
AntennaCollection<TimeDomainAntenna> detectorCoREAS;
// AntennaCollection<TimeDomainAntenna> detectorZHS;
std::string name_center = "CoREAS_R=0_m--Phi=0degrees";
auto triggertime_center{((triggerpoint_ - center).getNorm() / constants::c) - 500_ns};
......@@ -127,7 +121,7 @@ int main() {
auto rr_1 = static_cast<int>(radius_1 / 1_m);
auto const point_1{Point(rootCS, centerX + radius_1 * cos(phiRad_1),
centerY + radius_1 * sin(phiRad_1), centerZ)};
std::cout << "Antenna point: " << point_1 << std::endl;
CORSIKA_LOG_INFO("Antenna point: {}", point_1);
auto triggertime_1{((triggerpoint_ - point_1).getNorm() / constants::c) - 500_ns};
std::string name_1 = "CoREAS_R=" + std::to_string(rr_1) +
"_m--Phi=" + std::to_string(phi_1) + "degrees";
......@@ -186,16 +180,8 @@ int main() {
coreas(detectorCoREAS, SP);
output.add("CoREAS", coreas);
// RadioProcess<decltype(detectorZHS), ZHS<decltype(detectorZHS),
// decltype(SP)>, decltype(SP)>
// zhs(detectorZHS, SP);
// output.add("ZHS", zhs);
TimeCut cut(period / 4);
// TrackWriter trackWriter;
// output.add("tracks", trackWriter); // register TrackWriter
// assemble all processes into an ordered process list
auto sequence = make_sequence(coreas, cut);
......
examples/vertical_EAS.cpp deleted 100644 → 0
View file @ db6fe545
/*
* (c) Copyright 2022 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.
*/
#define TRACE
/* 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/framework/process/InteractionCounter.hpp>
/* clang-format on */
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>
#include <corsika/framework/geometry/Plane.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/process/InteractionCounter.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/utility/SaveBoostHistogram.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/output/OutputManager.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/GeomagneticModel.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMagneticFieldModel.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/Epos.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/Sophia.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/UrQMD.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
#include <iomanip>
#include <iostream>
#include <limits>
#include <string>
using namespace corsika;
using namespace std;
using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
using EnvType = Environment<EnvironmentInterface>;
using Particle = setup::Stack<EnvType>::particle_type;
void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("sibyll");
// RNGManager<>::getInstance().registerRandomStream("sophia");
RNGManager<>::getInstance().registerRandomStream("pythia");
RNGManager<>::getInstance().registerRandomStream("urqmd");
RNGManager<>::getInstance().registerRandomStream("proposal");
if (seed == 0) {
std::random_device rd;
seed = rd();
cout << "new random seed (auto) " << seed << endl;
}
RNGManager<>::getInstance().setSeed(seed);
}
template <typename T>
using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
// argv : 1.number of nucleons, 2.number of protons,
// 3.total energy in GeV, 4.number of showers,
// 5.seed (0 by default to generate random values for all)
int main(int argc, char** argv) {
logging::set_level(logging::level::warn);
CORSIKA_LOG_INFO("vertical_EAS");
if (argc < 4) {
std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> [seed] \n"
" if A=0, Z is interpreted as PDG code \n"
" if no seed is given, a random seed is chosen \n"
<< 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
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
// build a Linsley US Standard atmosphere into `env`
create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm,
MagneticFieldVector{rootCS, 20.4_uT, 0_T, -43.23_uT});
// Uncomment if you want to use PROPOSAL
// std::unordered_map<Code, HEPEnergyType> energy_resolution = {
// {Code::Electron, 2_MeV},
// {Code::Positron, 2_MeV},
// {Code::Photon, 2_MeV},
// };
// for (auto [pcode, energy] : energy_resolution)
// set_energy_production_threshold(pcode, energy);
// pre-setup particle stack
unsigned short const A = std::stoi(std::string(argv[1]));
Code beamCode;
unsigned short Z = 0;
if (A > 0) {
Z = std::stoi(std::string(argv[2]));
if (A == 1 && Z == 0)
beamCode = Code::Neutron;
else if (A == 1 && Z == 1)
beamCode = Code::Proton;
else
beamCode = get_nucleus_code(A, Z);
} else {
int pdg = std::stoi(std::string(argv[2]));
beamCode = convert_from_PDG(PDGCode(pdg));
}
HEPEnergyType mass = get_mass(beamCode);
HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
double theta = 0.;
double phi = 180.;
auto const thetaRad = theta / 180. * constants::pi;
auto const phiRad = phi / 180. * constants::pi;
HEPMomentumType P0 = calculate_momentum(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(thetaRad, phiRad, P0);
auto plab = 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
<< ", norm = " << plab.getNorm() << endl;
auto const observationHeight = 0.0_km + constants::EarthRadius::Mean;
auto const injectionHeight = 111.75_km + constants::EarthRadius::Mean;
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 + DirectionVector{rootCS,
{-sin(thetaRad) * cos(phiRad),
-sin(thetaRad) * sin(phiRad), cos(thetaRad)}} *
t;
std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
// 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).getNorm() * 1.5
<< std::endl;
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env, false,
1000};
// create the output manager that we then register outputs with
OutputManager output("vertical_EAS_outputs");
// EnergyLossWriterParquet dEdX_output{showerAxis, 10_g / square(1_cm), 200};
// register energy losses as output
// output.add("dEdX", dEdX_output);
// register profile output
// construct the overall energy loss writer and register it
EnergyLossWriter dEdX{showerAxis};
output.add("energyloss", dEdX);
// construct the continuous energy loss model
BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
// construct a particle cut - cuts are set to values close to reality, put
// higher values for faster runs
ParticleCut<SubWriter<decltype(dEdX)>> cut{2_MeV, 2_MeV, 2_GeV, 300_MeV, true, dEdX};
// setup longitudinal profile
LongitudinalWriter longProf{showerAxis};
output.add("profile", longProf);
LongitudinalProfile<SubWriter<decltype(longProf)>> profile{longProf};
// create a track writer and register it with the output manager
TrackWriter<TrackWriterParquet> trackWriter;
output.add("tracks", trackWriter);
// setup processes, decays and interactions
corsika::sibyll::Interaction sibyll{env};
InteractionCounter sibyllCounted{sibyll};
HEPEnergyType heThresholdNN = 60_GeV;
// PROPOSAL is disabled for this example
// corsika::sophia::InteractionModel sophia;
// corsika::proposal::Interaction emCascade(env, sophia,
// sibyll.getHadronInteractionModel(), heThresholdNN);
// corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>>
// emContinuous(env, dEdX);
corsika::pythia8::Decay decayPythia;
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted{urqmd};
StackInspector<setup::Stack<EnvType>> stackInspect(50000, false, E0);
// assemble all processes into an ordered process list
struct EnergySwitch {
HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {}
bool operator()(const Particle& p) const { return (p.getEnergy() < cutE_); }
};
auto hadronSequence = make_select(EnergySwitch(55_GeV), urqmdCounted, sibyllCounted);
// directory for outputs
string const labHist_file = "inthist_lab_verticalEAS.npz";
string const cMSHist_file = "inthist_cms_verticalEAS.npz";
// setup particle stack, and add primary particle
setup::Stack<EnvType> stack;
stack.clear();
stack.addParticle(std::make_tuple(
beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
plab.normalized(), injectionPos, 0_ns));
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{
obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
output.add("particles", observationLevel);
auto sequence = make_sequence(stackInspect, hadronSequence, decayPythia, emContinuous,
profile, observationLevel, cut);
// define air shower object, run simulation
setup::Tracking tracking;
output.startOfLibrary();
Cascade EAS(env, tracking, sequence, output, stack);
EAS.run();
auto const hists = sibyllCounted.getHistogram() + urqmdCounted.getHistogram();
save_hist(hists.labHist(), labHist_file, true);
save_hist(hists.CMSHist(), cMSHist_file, true);
output.endOfLibrary();
}

data privacy - impressum