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 @@ ...@@ -6,50 +6,43 @@
* the license. * the license.
*/ */
/* clang-format off */ #include <corsika/framework/core/Cascade.hpp>
// InteractionCounter used boost/histogram, which #include <corsika/framework/core/EnergyMomentumOperations.hpp>
// 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/Logging.hpp> #include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp> #include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/Step.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/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/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.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/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp> #include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/MediumPropertyModel.hpp> #include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/ShowerAxis.hpp> #include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp> #include <corsika/media/UniformMagneticField.hpp>
#include <corsika/modules/BetheBlochPDG.hpp> #include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp> #include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp> #include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/ParticleCut.hpp> #include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/Pythia8.hpp> #include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/Sibyll.hpp> #include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/UrQMD.hpp> #include <corsika/modules/UrQMD.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/CONEX.hpp> #include <corsika/modules/CONEX.hpp>
#include <corsika/setup/SetupStack.hpp> #include <corsika/setup/SetupStack.hpp>
...@@ -63,6 +56,12 @@ ...@@ -63,6 +56,12 @@
using namespace corsika; using namespace corsika;
using namespace std; 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 * Random number stream initialization
* *
...@@ -70,17 +69,19 @@ using namespace std; ...@@ -70,17 +69,19 @@ using namespace std;
*/ */
void registerRandomStreams(uint64_t seed) { void registerRandomStreams(uint64_t seed) {
RNGManager<>::getInstance().registerRandomStream("cascade"); RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("qgsjet"); RNGManager<>::getInstance().registerRandomStream("conex");
RNGManager<>::getInstance().registerRandomStream("sibyll");
RNGManager<>::getInstance().registerRandomStream("epos"); RNGManager<>::getInstance().registerRandomStream("epos");
RNGManager<>::getInstance().registerRandomStream("proposal");
RNGManager<>::getInstance().registerRandomStream("pythia"); RNGManager<>::getInstance().registerRandomStream("pythia");
RNGManager<>::getInstance().registerRandomStream("qgsjet");
RNGManager<>::getInstance().registerRandomStream("sibyll");
RNGManager<>::getInstance().registerRandomStream("urqmd"); RNGManager<>::getInstance().registerRandomStream("urqmd");
RNGManager<>::getInstance().registerRandomStream("proposal");
RNGManager<>::getInstance().registerRandomStream("conex");
if (seed == 0) { if (seed == 0) {
std::random_device rd; std::random_device rd;
seed = 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); RNGManager<>::getInstance().setSeed(seed);
} }
...@@ -141,8 +142,12 @@ private: ...@@ -141,8 +142,12 @@ private:
/** /**
* Selection of environment interface implementation: * Selection of environment interface implementation:
*/ */
using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
using EnvType = Environment<EnvironmentInterface>;
template <typename T> template <typename T>
using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>; using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
using StackType = setup::Stack<EnvType>;
using TrackingType = setup::Tracking;
int main(int argc, char** argv) { int main(int argc, char** argv) {
...@@ -157,7 +162,6 @@ 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"); " if no seed is given, a random seed is chosen");
return 1; return 1;
} }
feenableexcept(FE_INVALID);
uint64_t seed = 0; uint64_t seed = 0;
if (argc > 4) seed = std::stol(std::string(argv[4])); if (argc > 4) seed = std::stol(std::string(argv[4]));
...@@ -165,20 +169,15 @@ int main(int argc, char** argv) { ...@@ -165,20 +169,15 @@ int main(int argc, char** argv) {
registerRandomStreams(seed); registerRandomStreams(seed);
// setup environment, geometry // setup environment, geometry
using EnvironmentInterface = IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
using EnvType = Environment<EnvironmentInterface>;
EnvType env; EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem(); CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m}; Point const center{rootCS, 0_m, 0_m, 0_m};
// build a Linsley US Standard atmosphere into `env` // build a Linsley US Standard atmosphere into `env`
MagneticFieldVector bField{rootCS, 50_uT, 0_T, 0_T};
create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>( create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm, env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm, bField);
MagneticFieldVector{rootCS, 0_T, 50_uT, 0_T});
// 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 A = std::stoi(std::string(argv[1]));
unsigned short const Z = std::stoi(std::string(argv[2])); unsigned short const Z = std::stoi(std::string(argv[2]));
Code const beamCode = get_nucleus_code(A, Z); Code const beamCode = get_nucleus_code(A, Z);
...@@ -194,12 +193,6 @@ int main(int argc, char** argv) { ...@@ -194,12 +193,6 @@ int main(int argc, char** argv) {
auto const [px, py, pz] = momentumComponents(thetaRad, P0); auto const [px, py, pz] = momentumComponents(thetaRad, P0);
auto plab = MomentumVector(rootCS, {px, py, pz}); 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 observationHeight = 0_km + constants::EarthRadius::Mean;
auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean; auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
...@@ -211,28 +204,25 @@ int main(int argc, char** argv) { ...@@ -211,28 +204,25 @@ int main(int argc, char** argv) {
showerCore + showerCore +
Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t; 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, ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000}; 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}; OutputManager output("hybrid_MC_outputs");
InteractionCounter sibyllCounted{sibyll};
corsika::pythia8::Decay decayPythia;
// register energy losses as output // register energy losses as output
EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200}; EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
output.add("energyloss", dEdX); output.add("energyloss", dEdX);
// create a track writer and register it with the output manager // create a track writer and register it with the output manager
...@@ -242,19 +232,29 @@ int main(int argc, char** argv) { ...@@ -242,19 +232,29 @@ int main(int argc, char** argv) {
ParticleCut<SubWriter<decltype(dEdX)>> cut(3_GeV, false, dEdX); ParticleCut<SubWriter<decltype(dEdX)>> cut(3_GeV, false, dEdX);
BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss(dEdX); BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss(dEdX);
LongitudinalWriter profile{showerAxis, 10_g / square(1_cm)}; LongitudinalWriter profile{showerAxis, nAxisBins, dX};
output.add("profile", profile); output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{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)>> CONEXhybrid // SubWriter<decltype(dEdX>, SubWriter<decltype(profile)>>
conex_model(center, showerAxis, t, injectionHeight, E0, get_PDG(Code::Proton), dEdX, conex_model(center, showerAxis, t, injectionHeight, E0, get_PDG(Code::Proton), dEdX,
profile); 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; corsika::urqmd::UrQMD urqmd_model;
InteractionCounter urqmdCounted{urqmd_model}; InteractionCounter urqmdCounted{urqmd_model};
...@@ -265,7 +265,7 @@ int main(int argc, char** argv) { ...@@ -265,7 +265,7 @@ int main(int argc, char** argv) {
HEPEnergyType cutE_; HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE) EnergySwitch(HEPEnergyType cutE)
: cutE_(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_); return (p.getEnergy() < cutE_);
} }
}; };
...@@ -273,17 +273,27 @@ int main(int argc, char** argv) { ...@@ -273,17 +273,27 @@ int main(int argc, char** argv) {
auto sequence = make_sequence(hadronSequence, decayPythia, eLoss, cut, conex_model, auto sequence = make_sequence(hadronSequence, decayPythia, eLoss, cut, conex_model,
longprof, observationLevel, trackCheck); longprof, observationLevel, trackCheck);
StackType stack;
stack.clear();
// define air shower object, run simulation // define air shower object, run simulation
setup::Tracking tracking; TrackingType tracking;
Cascade EAS(env, tracking, sequence, output, stack); 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: // to fix the point of first interaction, uncomment the following two lines:
// EAS.SetNodes(); // EAS.SetNodes();
// EAS.forceInteraction(); // EAS.forceInteraction();
output.startOfLibrary();
EAS.run(); EAS.run();
output.endOfLibrary();
const HEPEnergyType Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround(); const HEPEnergyType Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
CORSIKA_LOG_INFO( CORSIKA_LOG_INFO(
...@@ -296,5 +306,7 @@ int main(int argc, char** argv) { ...@@ -296,5 +306,7 @@ int main(int argc, char** argv) {
save_hist(hists.labHist(), "inthist_lab_hybrid.npz", true); save_hist(hists.labHist(), "inthist_lab_hybrid.npz", true);
save_hist(hists.CMSHist(), "inthist_cms_hybrid.npz", true); save_hist(hists.CMSHist(), "inthist_cms_hybrid.npz", true);
output.endOfLibrary();
CORSIKA_LOG_INFO("done"); CORSIKA_LOG_INFO("done");
} }
examples/radio_em_shower.cpp → examples/cascade_examples/radio_em_shower.cpp
View file @ 0c1d3948
...@@ -21,10 +21,11 @@ ...@@ -21,10 +21,11 @@
#include <corsika/framework/utility/CorsikaFenv.hpp> #include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/framework/utility/SaveBoostHistogram.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/EnergyLossWriter.hpp>
#include <corsika/modules/writers/LongitudinalWriter.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/Environment.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp> #include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
...@@ -64,6 +65,13 @@ ...@@ -64,6 +65,13 @@
using namespace corsika; using namespace corsika;
using namespace std; 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) { void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade"); RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("proposal"); RNGManager<>::getInstance().registerRandomStream("proposal");
...@@ -72,14 +80,21 @@ void registerRandomStreams(int seed) { ...@@ -72,14 +80,21 @@ void registerRandomStreams(int seed) {
if (seed == 0) { if (seed == 0) {
std::random_device rd; std::random_device rd;
seed = 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); RNGManager<>::getInstance().setSeed(seed);
} }
using EnvironmentInterface =
IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
using EnvType = Environment<EnvironmentInterface>;
template <typename TInterface> template <typename TInterface>
using MyExtraEnv = using MyExtraEnv =
UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>; UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
using StackType = setup::Stack<EnvType>;
using TrackingType = setup::Tracking;
int main(int argc, char** argv) { int main(int argc, char** argv) {
...@@ -93,22 +108,21 @@ 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])))}; 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); feenableexcept(FE_INVALID);
// initialize random number sequence(s) // initialize random number sequence(s)
registerRandomStreams(seed); registerRandomStreams(seed);
// setup environment, geometry // setup environment, geometry
using EnvironmentInterface =
IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
using EnvType = Environment<EnvironmentInterface>;
EnvType env; EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem(); CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m}; 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>( create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
env, AtmosphereId::LinsleyUSStd, center, 1.000327, Medium::AirDry1Atm, env, AtmosphereId::LinsleyUSStd, center, refractive_index, Medium::AirDry1Atm,
MagneticFieldVector{rootCS, 50_uT, 0_T, 0_T}); bField);
std::unordered_map<Code, HEPEnergyType> energy_resolution = { std::unordered_map<Code, HEPEnergyType> energy_resolution = {
{Code::Electron, 5_MeV}, {Code::Electron, 5_MeV},
...@@ -118,26 +132,19 @@ int main(int argc, char** argv) { ...@@ -118,26 +132,19 @@ int main(int argc, char** argv) {
for (auto [pcode, energy] : energy_resolution) for (auto [pcode, energy] : energy_resolution)
set_energy_production_threshold(pcode, energy); set_energy_production_threshold(pcode, energy);
// setup particle stack, and add primary particle Code const beamCode = Code::Electron;
setup::Stack<EnvType> stack;
stack.clear();
const Code beamCode = Code::Electron;
auto const mass = get_mass(beamCode); auto const mass = get_mass(beamCode);
const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1])); HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[1]));
double theta = 0.; double const theta = 0.;
auto const thetaRad = theta / 180. * M_PI; 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) { auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad)); return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
}; };
auto const [px, py, pz] = momentumComponents(thetaRad, P0); auto const [px, py, pz] = momentumComponents(thetaRad, P0);
auto plab = MomentumVector(rootCS, {px, py, pz}); auto const 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 observationHeight = 1.4_km + constants::EarthRadius::Mean; auto const observationHeight = 1.4_km + constants::EarthRadius::Mean;
auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean; auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
...@@ -148,104 +155,103 @@ int main(int argc, char** argv) { ...@@ -148,104 +155,103 @@ int main(int argc, char** argv) {
Point const injectionPos = Point const injectionPos =
showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t; 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, ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000}; 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}; 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 // Radio antennas and relevant information
// the antenna time variables // the antenna time variables
const TimeType duration_{1e-6_s}; TimeType const duration{1e-6_s};
const InverseTimeType sampleRate_{1e+9_Hz}; InverseTimeType const sampleRate{1e+9_Hz};
// the detector (aka antenna collection) for CoREAS and ZHS // the detector (aka antenna collection) for CoREAS and ZHS
AntennaCollection<TimeDomainAntenna> detectorCoREAS; AntennaCollection<TimeDomainAntenna> detectorCoREAS;
AntennaCollection<TimeDomainAntenna> detectorZHS; AntennaCollection<TimeDomainAntenna> detectorZHS;
auto const showerCoreX_{showerCore.getCoordinates().getX()}; auto const showerCoreX{showerCore.getCoordinates().getX()};
auto const showerCoreY_{showerCore.getCoordinates().getY()}; auto const showerCoreY{showerCore.getCoordinates().getY()};
auto const injectionPosX_{injectionPos.getCoordinates().getX()}; auto const injectionPosX{injectionPos.getCoordinates().getX()};
auto const injectionPosY_{injectionPos.getCoordinates().getY()}; auto const injectionPosY{injectionPos.getCoordinates().getY()};
auto const injectionPosZ_{injectionPos.getCoordinates().getZ()}; auto const injectionPosZ{injectionPos.getCoordinates().getZ()};
auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)}; auto const triggerpoint{Point(rootCS, injectionPosX, injectionPosY, injectionPosZ)};
std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
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 // // 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 radius_coreas = 25_m; radius_coreas <= 500_m; radius_coreas += 25_m) {
// for (auto phi_1 = 0; phi_1 <= 315; phi_1 += 45) { // for (auto phi_coreas = 0; phi_coreas <= 315; phi_coreas += 45) {
auto radius_1 = 200_m; auto radius_coreas = 200_m;
auto phi_1 = 45; auto phi_coreas = 45;
auto phiRad_1 = phi_1 / 180. * M_PI; auto phiRad_coreas = phi_coreas / 180. * M_PI;
auto rr_1 = static_cast<int>(radius_1 / 1_m); auto rr_coreas = static_cast<int>(radius_coreas / 1_m);
auto const point_1{Point(rootCS, showerCoreX_ + radius_1 * cos(phiRad_1), auto const point_coreas{Point(rootCS, showerCoreX + radius_coreas * cos(phiRad_coreas),
showerCoreY_ + radius_1 * sin(phiRad_1), showerCoreY + radius_coreas * sin(phiRad_coreas),
constants::EarthRadius::Mean)}; constants::EarthRadius::Mean)};
std::cout << "Antenna point: " << point_1 << std::endl; std::cout << "Antenna point: " << point_coreas << std::endl;
auto triggertime_1{(triggerpoint_ - point_1).getNorm() / constants::c}; CORSIKA_LOG_INFO("Antenna point {}", injectionPos.getCoordinates());
std::string name_1 = auto triggertime_coreas{(triggerpoint - point_coreas).getNorm() / constants::c};
"CoREAS_R=" + std::to_string(rr_1) + "_m--Phi=" + std::to_string(phi_1) + "degrees"; std::string name_coreas = "CoREAS_R=" + std::to_string(rr_coreas) +
TimeDomainAntenna antenna_1(name_1, point_1, rootCS, triggertime_1, duration_, "_m--Phi=" + std::to_string(phi_coreas) + "degrees";
sampleRate_, triggertime_1); TimeDomainAntenna antenna_coreas(name_coreas, point_coreas, rootCS, triggertime_coreas,
detectorCoREAS.addAntenna(antenna_1); duration, sampleRate, triggertime_coreas);
detectorCoREAS.addAntenna(antenna_coreas);
// } // }
// } // }
// // setup ZHS antennas - use the for loop for star shape pattern // // setup ZHS antennas - use the for loop for star shape pattern
// for (auto radius_ = 25_m; radius_ <= 500_m; radius_ += 25_m) { // for (auto radius_zhs = 25_m; radius_zhs <= 500_m; radius_zhs += 25_m) {
// for (auto phi_ = 0; phi_ <= 315; phi_ += 45) { // for (auto phi_zhs = 0; phi_zhs <= 315; phi_zhs += 45) {
auto radius_ = 200_m; auto radius_zhs = 200_m;
auto phi_ = 45; auto phi_zhs = 45;
auto phiRad_ = phi_ / 180. * M_PI; auto phiRad_zhs = phi_zhs / 180. * M_PI;
auto rr_ = static_cast<int>(radius_ / 1_m); auto rr_zhs = static_cast<int>(radius_zhs / 1_m);
auto const point_{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_), auto const point_zhs{Point(rootCS, showerCoreX + radius_zhs * cos(phiRad_zhs),
showerCoreY_ + radius_ * sin(phiRad_), showerCoreY + radius_zhs * sin(phiRad_zhs),
constants::EarthRadius::Mean)}; constants::EarthRadius::Mean)};
auto triggertime_{(triggerpoint_ - point_).getNorm() / constants::c}; auto triggertime_zhs{(triggerpoint - point_zhs).getNorm() / constants::c};
std::string name_ = std::string name_zhs = "ZHS_R=" + std::to_string(rr_zhs) +
"ZHS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees"; "_m--Phi=" + std::to_string(phi_zhs) + "degrees";
TimeDomainAntenna antenna_(name_, point_, rootCS, triggertime_, duration_, sampleRate_, TimeDomainAntenna antenna_zhs(name_zhs, point_zhs, rootCS, triggertime_zhs, duration,
triggertime_); sampleRate, triggertime_zhs);
detectorZHS.addAntenna(antenna_); detectorZHS.addAntenna(antenna_zhs);
// } // }
// } // }
// setup processes, decays and interactions // setup processes, decays and interactions
EnergyLossWriter energyloss{showerAxis, dX, nAxisBins};
EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200}; ParticleCut<SubWriter<decltype(energyloss)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true,
// register energy losses as output energyloss);
output.add("dEdX", dEdX);
ParticleCut<SubWriter<decltype(dEdX)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true, dEdX);
corsika::sophia::InteractionModel sophia;
corsika::sibyll::Interaction sibyll{env}; corsika::sibyll::Interaction sibyll{env};
corsika::sophia::InteractionModel sophia;
HEPEnergyType heThresholdNN = 80_GeV; HEPEnergyType heThresholdNN = 80_GeV;
corsika::proposal::Interaction emCascade( corsika::proposal::Interaction emCascade(
env, sophia, sibyll.getHadronInteractionModel(), heThresholdNN); env, sophia, sibyll.getHadronInteractionModel(), heThresholdNN);
corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX); corsika::proposal::ContinuousProcess<SubWriter<decltype(energyloss)>> emContinuous(
// BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX}; env, energyloss);
// NOT possible right now, due to interface differenc in PROPOSAL // NOT possible right now, due to interface differenc in PROPOSAL
// InteractionCounter emCascadeCounted(emCascade); // InteractionCounter emCascadeCounted(emCascade);
OutputManager output("radio_em_shower_outputs");
output.add("energyloss", energyloss);
TrackWriter tracks; TrackWriter tracks;
output.add("tracks", tracks); output.add("tracks", tracks);
// long. profile LongitudinalWriter profile{showerAxis, nAxisBins, dX};
LongitudinalWriter profile{showerAxis, 10_g / square(1_cm)};
output.add("profile", profile); output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile}; LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
...@@ -269,15 +275,27 @@ int main(int argc, char** argv) { ...@@ -269,15 +275,27 @@ int main(int argc, char** argv) {
output.add("ZHS", zhs); output.add("ZHS", zhs);
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.})); Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{ ObservationPlane<TrackingType, ParticleWriterParquet> observationLevel{
obsPlane, DirectionVector(rootCS, {1., 0., 0.})}; obsPlane, DirectionVector(rootCS, {1., 0., 0.})};
output.add("particles", observationLevel); 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, auto sequence = make_sequence(emCascade, emContinuous, longprof, coreas, zhs, tracks,
observationLevel, cut); observationLevel, cut);
// define air shower object, run simulation // 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(); output.startOfLibrary();
Cascade EAS(env, tracking, sequence, output, stack); Cascade EAS(env, tracking, sequence, output, stack);
...@@ -287,7 +305,8 @@ int main(int argc, char** argv) { ...@@ -287,7 +305,8 @@ int main(int argc, char** argv) {
EAS.run(); EAS.run();
HEPEnergyType const Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround(); HEPEnergyType const Efinal =
energyloss.getEnergyLost() + observationLevel.getEnergyGround();
CORSIKA_LOG_INFO( CORSIKA_LOG_INFO(
"total energy budget (GeV): {}, " "total energy budget (GeV): {}, "
...@@ -295,4 +314,6 @@ int main(int argc, char** argv) { ...@@ -295,4 +314,6 @@ int main(int argc, char** argv) {
Efinal / 1_GeV, (Efinal / E0 - 1) * 100); Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
output.endOfLibrary(); output.endOfLibrary();
return EXIT_SUCCESS;
} }
examples/water.cpp → examples/cascade_examples/water.cpp
View file @ 0c1d3948
...@@ -24,15 +24,16 @@ ...@@ -24,15 +24,16 @@
#include <corsika/media/ShowerAxis.hpp> #include <corsika/media/ShowerAxis.hpp>
#include <corsika/modules/ObservationPlane.hpp> #include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/LongitudinalProfile.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/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/PROPOSAL.hpp>
#include <corsika/modules/ParticleCut.hpp> #include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/Pythia8.hpp> #include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/Sibyll.hpp> #include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/Sophia.hpp> #include <corsika/modules/Sophia.hpp>
#include <corsika/modules/UrQMD.hpp> #include <corsika/modules/FLUKA.hpp>
#include <corsika/modules/tracking/TrackingStraight.hpp> #include <corsika/modules/tracking/TrackingStraight.hpp>
#include <corsika/output/OutputManager.hpp> #include <corsika/output/OutputManager.hpp>
...@@ -41,6 +42,7 @@ ...@@ -41,6 +42,7 @@
#include <corsika/stack/VectorStack.hpp> #include <corsika/stack/VectorStack.hpp>
#include <corsika/setup/SetupStack.hpp> #include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
#include <CLI/App.hpp> #include <CLI/App.hpp>
#include <CLI/Config.hpp> #include <CLI/Config.hpp>
...@@ -50,9 +52,9 @@ using namespace corsika; ...@@ -50,9 +52,9 @@ using namespace corsika;
using IMediumType = IMediumPropertyModel<IMediumModel>; using IMediumType = IMediumPropertyModel<IMediumModel>;
using EnvType = Environment<IMediumType>; using EnvType = Environment<IMediumType>;
using StackActive = setup::Stack<EnvType>; using StackType = setup::Stack<EnvType>;
using StackView = StackActive::stack_view_type; using TrackingType = tracking_line::Tracking;
using Particle = StackActive::particle_type; using Particle = StackType::particle_type;
void registerRandomStreams(int seed) { void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade"); RNGManager<>::getInstance().registerRandomStream("cascade");
...@@ -61,7 +63,7 @@ void registerRandomStreams(int seed) { ...@@ -61,7 +63,7 @@ void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("sophia"); RNGManager<>::getInstance().registerRandomStream("sophia");
RNGManager<>::getInstance().registerRandomStream("epos"); RNGManager<>::getInstance().registerRandomStream("epos");
RNGManager<>::getInstance().registerRandomStream("pythia"); RNGManager<>::getInstance().registerRandomStream("pythia");
RNGManager<>::getInstance().registerRandomStream("urqmd"); RNGManager<>::getInstance().registerRandomStream("fluka");
RNGManager<>::getInstance().registerRandomStream("proposal"); RNGManager<>::getInstance().registerRandomStream("proposal");
if (seed == 0) { if (seed == 0) {
std::random_device rd; std::random_device rd;
...@@ -75,12 +77,12 @@ void registerRandomStreams(int seed) { ...@@ -75,12 +77,12 @@ void registerRandomStreams(int seed) {
int main(int argc, char** argv) { int main(int argc, char** argv) {
// * process input // * process input
Code primaryType; Code beamCode;
HEPEnergyType e0, eCut; HEPEnergyType E0, eCut;
int A, Z, n_event; int A, Z, n_event;
int randomSeed; int randomSeed;
std::string output_dir; 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 // we start by definining a sub-group for the primary ID
auto opt_Z = app.add_option("-Z", Z, "Atomic number for primary") auto opt_Z = app.add_option("-Z", Z, "Atomic number for primary")
->check(CLI::Range(0, 26)) ->check(CLI::Range(0, 26))
...@@ -110,8 +112,25 @@ int main(int argc, char** argv) { ...@@ -110,8 +112,25 @@ int main(int argc, char** argv) {
app.add_option("-s", randomSeed, "Seed for random number") app.add_option("-s", randomSeed, "Seed for random number")
->check(CLI::NonNegativeNumber) ->check(CLI::NonNegativeNumber)
->default_val(0); ->default_val(0);
// parse the command line options into the variables
CLI11_PARSE(app, argc, argv); 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 // check that we got either PDG or A/Z
// this can be done with option_groups but the ordering // this can be done with option_groups but the ordering
// gets all messed up // gets all messed up
...@@ -121,37 +140,24 @@ int main(int argc, char** argv) { ...@@ -121,37 +140,24 @@ int main(int argc, char** argv) {
return 1; return 1;
} }
} }
// initialize random number sequence(s)
registerRandomStreams(randomSeed);
// check if we want to use a PDG code instead // check if we want to use a PDG code instead
if (app.count("--pdg") > 0) { if (app.count("--pdg") > 0) {
primaryType = convert_from_PDG(PDGCode(app["--pdg"]->as<int>())); beamCode = convert_from_PDG(PDGCode(app["--pdg"]->as<int>()));
} else { } else {
// check manually for proton and neutrons // check manually for proton and neutrons
if ((A == 1) && (Z == 1)) if ((A == 1) && (Z == 1))
primaryType = Code::Proton; beamCode = Code::Proton;
else if ((A == 1) && (Z == 0)) else if ((A == 1) && (Z == 0))
primaryType = Code::Neutron; beamCode = Code::Neutron;
else 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; 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 // * environment and universe
EnvType env; EnvType env;
...@@ -163,27 +169,26 @@ int main(int argc, char** argv) { ...@@ -163,27 +169,26 @@ int main(int argc, char** argv) {
Point const center{rootCS, 0_m, 0_m, 0_m}; Point const center{rootCS, 0_m, 0_m, 0_m};
auto sphere = std::make_unique<Sphere>(center, 100_m); auto sphere = std::make_unique<Sphere>(center, 100_m);
auto node = std::make_unique<VolumeTreeNode<IMediumType>>(std::move(sphere)); auto node = std::make_unique<VolumeTreeNode<IMediumType>>(std::move(sphere));
// Hydrogen is not supported by UrQMD yet. See #456 NuclearComposition const nuclearComposition{{Code::Hydrogen, Code::Oxygen},
auto comp = NuclearComposition({{Code::Oxygen}, {1.0}}); {2.0 / 3.0, 1.0 / 3.0}};
// density of sea water // density of sea water
auto density = 1.02_g / (1_cm * 1_cm * 1_cm); auto density = 1.02_g / (1_cm * 1_cm * 1_cm);
auto water_medium = auto water_medium =
std::make_shared<MediumPropertyModel<HomogeneousMedium<IMediumType>>>( std::make_shared<MediumPropertyModel<HomogeneousMedium<IMediumType>>>(
Medium::WaterLiquid, density, comp); Medium::WaterLiquid, density, nuclearComposition);
node->setModelProperties(water_medium); node->setModelProperties(water_medium);
universe->addChild(std::move(node)); universe->addChild(std::move(node));
} }
// * detector geometry // * make downward-going shower axis and a observation plane in x-y-plane
auto injectorLength = 50_m; auto injectorLength = 50_m;
Point const injectorPos = Point(rootCS, {0_m, 0_m, injectorLength}); Point const injectionPos = Point(rootCS, {0_m, 0_m, injectorLength});
auto const& injectCS = make_translation(rootCS, injectorPos.getCoordinates()); auto const& injectCS = make_translation(rootCS, injectionPos.getCoordinates());
DirectionVector upVec(rootCS, {0., 0., 1.}); DirectionVector upVec(rootCS, {0., 0., 1.});
DirectionVector leftVec(rootCS, {1., 0., 0.}); DirectionVector leftVec(rootCS, {1., 0., 0.});
DirectionVector downVec(rootCS, {0., 0., -1.}); DirectionVector downVec(rootCS, {0., 0., -1.});
// * observation plane std::vector<ObservationPlane<TrackingType, ParticleWriterParquet>> obsPlanes;
std::vector<ObservationPlane<tracking_line::Tracking>> obsPlanes;
const int nPlane = 5; const int nPlane = 5;
for (int i = 0; i < nPlane - 1; i++) { for (int i = 0; i < nPlane - 1; i++) {
Point planeCenter{injectCS, {0_m, 0_m, -(i + 1) * 3_m}}; Point planeCenter{injectCS, {0_m, 0_m, -(i + 1) * 3_m}};
...@@ -193,12 +198,14 @@ int main(int argc, char** argv) { ...@@ -193,12 +198,14 @@ int main(int argc, char** argv) {
Plane{Point{injectCS, {0_m, 0_m, -50_m}}, upVec}, leftVec, true); Plane{Point{injectCS, {0_m, 0_m, -50_m}}, upVec}, leftVec, true);
// * longitutional profile // * longitutional profile
ShowerAxis const showerAxis{injectorPos, 1.2 * injectorLength * downVec, env}; ShowerAxis const showerAxis{injectionPos, 1.2 * injectorLength * downVec, env};
LongitudinalWriter longiWriter{showerAxis, 5500, 1_g / square(1_cm)}; 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}; LongitudinalProfile<SubWriter<decltype(longiWriter)>> longprof{longiWriter};
// * energy loss profile // * energy loss profile
EnergyLossWriter dEdX{showerAxis, 1_g / square(1_cm), 5500}; EnergyLossWriter dEdX{showerAxis, dX, nAxisBins};
// * physical process list // * physical process list
// particle production threshold // particle production threshold
...@@ -207,19 +214,19 @@ int main(int argc, char** argv) { ...@@ -207,19 +214,19 @@ int main(int argc, char** argv) {
ParticleCut<SubWriter<decltype(dEdX)>> cut(emCut, emCut, hadCut, hadCut, true, dEdX); ParticleCut<SubWriter<decltype(dEdX)>> cut(emCut, emCut, hadCut, hadCut, true, dEdX);
// hadronic interactions // hadronic interactions
HEPEnergyType heHadronModelThreshold = 63.1_GeV; HEPEnergyType heHadronModelThreshold = std::pow(10, 1.9) * 1_GeV;
corsika::sibyll::Interaction sibyll(env); corsika::sibyll::Interaction sibyll(env);
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted(urqmd); corsika::fluka::Interaction leIntModel{env};
InteractionCounter leIntCounted{leIntModel};
struct EnergySwitch { struct EnergySwitch {
HEPEnergyType cutE_; HEPEnergyType cutE_;
EnergySwitch(HEPEnergyType cutE) EnergySwitch(HEPEnergyType cutE)
: cutE_(cutE) {} : cutE_(cutE) {}
bool operator()(const Particle& p) const { return (p.getKineticEnergy() < cutE_); } bool operator()(const Particle& p) const { return (p.getKineticEnergy() < cutE_); }
}; };
// auto lowModel = make_sequence(urqmd);
auto hadronSequence = auto hadronSequence =
make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, sibyll); make_select(EnergySwitch(heHadronModelThreshold), leIntCounted, sibyll);
// decay process // decay process
corsika::pythia8::Decay decayPythia; corsika::pythia8::Decay decayPythia;
...@@ -243,19 +250,41 @@ int main(int argc, char** argv) { ...@@ -243,19 +250,41 @@ int main(int argc, char** argv) {
// hard coded // hard coded
auto obsPlaneSequence = auto obsPlaneSequence =
make_sequence(obsPlanes[0], obsPlanes[1], obsPlanes[2], obsPlanes[3], obsPlanes[4]); 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 // * the final process sequence
auto sequence = make_sequence(physics_sequence, longprof, obsPlaneSequence, cut); auto sequence = make_sequence(physics_sequence, longprof, obsPlaneSequence, cut);
// * tracking and stack // * tracking and stack
tracking_line::Tracking tracking; TrackingType tracking;
StackActive stack; StackType stack;
// * cascade manager // * cascade manager
Cascade EAS(env, tracking, sequence, output, stack); 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 // * main loop
output.startOfLibrary(); output.startOfLibrary();
for (int i_shower = 0; i_shower < n_event; i_shower++) { for (int i_shower = 0; i_shower < n_event; i_shower++) {
...@@ -263,8 +292,11 @@ int main(int argc, char** argv) { ...@@ -263,8 +292,11 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO("Event: {} / {}", i_shower, n_event); CORSIKA_LOG_INFO("Event: {} / {}", i_shower, n_event);
// * inject primary // * inject primary
auto primary = stack.addParticle(std::make_tuple( auto const primaryProperties = std::make_tuple(
primaryType, e0 - get_mass(primaryType), downVec, injectorPos, 0_ns)); beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
plab.normalized(), injectionPos, 0_ns);
auto primary = stack.addParticle(primaryProperties);
stack.addParticle(primaryProperties);
EAS.run(); EAS.run();
...@@ -273,8 +305,8 @@ int main(int argc, char** argv) { ...@@ -273,8 +305,8 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO( CORSIKA_LOG_INFO(
"total energy budget (TeV): {:.2f} (dEdX={:.2f} ground={:.2f}), " "total energy budget (TeV): {:.2f} (dEdX={:.2f} ground={:.2f}), "
"relative difference (%): {:.3f}", "relative difference (%): {:.3f}",
e0 / 1_TeV, dEdX.getEnergyLost() / 1_TeV, obsPlaneFinal.getEnergyGround() / 1_TeV, E0 / 1_TeV, dEdX.getEnergyLost() / 1_TeV, obsPlaneFinal.getEnergyGround() / 1_TeV,
(Efinal / e0 - 1.) * 100.); (Efinal / E0 - 1.) * 100.);
} }
output.endOfLibrary(); 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 @@ ...@@ -37,6 +37,13 @@
using namespace corsika; using namespace corsika;
using namespace std; 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> template <bool deleteParticle>
struct MyBoundaryCrossingProcess struct MyBoundaryCrossingProcess
: public BoundaryCrossingProcess<MyBoundaryCrossingProcess<deleteParticle>> { : public BoundaryCrossingProcess<MyBoundaryCrossingProcess<deleteParticle>> {
...@@ -65,9 +72,6 @@ private: ...@@ -65,9 +72,6 @@ private:
std::ofstream file_; std::ofstream file_;
}; };
//
// The example main program for a particle cascade
//
int main() { int main() {
logging::set_level(logging::level::warn); 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; ...@@ -18,11 +18,14 @@ using namespace corsika;
using namespace std; 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() { int main() {
logging::set_level(logging::level::warn); logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] %v"); corsika_logger->set_pattern("[%n:%^%-8l%$] %v");
logging::info( logging::info(
...@@ -30,7 +33,7 @@ int main() { ...@@ -30,7 +33,7 @@ int main() {
"------------------------------------------\n" "------------------------------------------\n"
" particles in CORSIKA\n" " particles in CORSIKA\n"
"------------------------------------------\n"); "------------------------------------------\n");
int const width = 20 + 10 + 10 + 10 + 15 + 15 + 17; int const width = 20 + 10 + 10 + 10 + 16 + 16 + 17;
logging::info( logging::info(
"Name | " "Name | "
"PDG-id | " "PDG-id | "
...@@ -46,7 +49,7 @@ int main() { ...@@ -46,7 +49,7 @@ int main() {
? to_string(corsika::sibyll::getSibyllMass(p) / 1_GeV) ? to_string(corsika::sibyll::getSibyllMass(p) / 1_GeV)
: ""); : "");
auto const qgs_id = corsika::qgsjetII::convertToQgsjetII(p); 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)), static_cast<int>(get_PDG(p)),
(sib_id != corsika::sibyll::SibyllCode::Unknown (sib_id != corsika::sibyll::SibyllCode::Unknown
? to_string(static_cast<int>(sib_id)) ? 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() { ...@@ -131,4 +131,4 @@ int main() {
modular(); modular();
return 0; 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 @@ ...@@ -9,7 +9,6 @@
#include <corsika/media/Environment.hpp> #include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp> #include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMediumModel.hpp> #include <corsika/media/IMediumModel.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/framework/geometry/Sphere.hpp> #include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/modules/BetheBlochPDG.hpp> #include <corsika/modules/BetheBlochPDG.hpp>
...@@ -30,7 +29,7 @@ using namespace std; ...@@ -30,7 +29,7 @@ using namespace std;
// //
int main() { int main() {
logging::set_level(logging::level::warn); logging::set_level(logging::level::info);
CORSIKA_LOG_INFO("stopping_power"); CORSIKA_LOG_INFO("stopping_power");
...@@ -52,7 +51,9 @@ int main() { ...@@ -52,7 +51,9 @@ int main() {
setup::Stack<EnvType> stack; 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; file << "# beta*gamma, dE/dX / MeV/(g/cm²)" << std::endl;
for (HEPEnergyType E0 = 200_MeV; E0 < 1_PeV; E0 *= 1.05) { for (HEPEnergyType E0 = 200_MeV; E0 < 1_PeV; E0 *= 1.05) {
...@@ -79,5 +80,5 @@ int main() { ...@@ -79,5 +80,5 @@ int main() {
HEPEnergyType dE = eLoss.getTotalEnergyLoss(p, 1_g / square(1_cm)); HEPEnergyType dE = eLoss.getTotalEnergyLoss(p, 1_g / square(1_cm));
file << P0 / mass << "\t" << -dE / 1_MeV << std::endl; 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 @@ ...@@ -24,7 +24,6 @@
#include <corsika/media/MediumPropertyModel.hpp> #include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp> #include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/UniformRefractiveIndex.hpp> #include <corsika/media/UniformRefractiveIndex.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/setup/SetupStack.hpp> #include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp> #include <corsika/setup/SetupTrajectory.hpp>
...@@ -35,14 +34,10 @@ ...@@ -35,14 +34,10 @@
#include <corsika/modules/radio/antennas/Antenna.hpp> #include <corsika/modules/radio/antennas/Antenna.hpp>
#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp> #include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
#include <corsika/modules/radio/detectors/AntennaCollection.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/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/TimeCut.hpp>
// #include <corsika/modules/TrackWriter.hpp>
#include <iomanip> #include <iomanip>
#include <iostream> #include <iostream>
...@@ -57,6 +52,7 @@ using namespace std; ...@@ -57,6 +52,7 @@ using namespace std;
// A simple shower to get the electric field trace of an electron (& a positron) // 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. // in order to reveal the "clover leaf" pattern of energy fluence to the ground.
// //
int main() { int main() {
logging::set_level(logging::level::info); logging::set_level(logging::level::info);
...@@ -70,7 +66,7 @@ int main() { ...@@ -70,7 +66,7 @@ int main() {
auto seed = rd(); auto seed = rd();
RNGManager<>::getInstance().setSeed(seed); RNGManager<>::getInstance().setSeed(seed);
OutputManager output("1 electron - 1 positron"); OutputManager output("clover_leaf_outputs");
// create a suitable environment // create a suitable environment
using IModelInterface = using IModelInterface =
...@@ -105,7 +101,6 @@ int main() { ...@@ -105,7 +101,6 @@ int main() {
auto const injectionPosY_{injectionPos.getCoordinates().getY()}; auto const injectionPosY_{injectionPos.getCoordinates().getY()};
auto const injectionPosZ_{injectionPos.getCoordinates().getZ()}; auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)}; auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
// the antenna characteristics // the antenna characteristics
const TimeType duration_{2e-6_s}; // 0.8e-4_s const TimeType duration_{2e-6_s}; // 0.8e-4_s
...@@ -113,7 +108,6 @@ int main() { ...@@ -113,7 +108,6 @@ int main() {
// the detectors // the detectors
AntennaCollection<TimeDomainAntenna> detectorCoREAS; AntennaCollection<TimeDomainAntenna> detectorCoREAS;
// AntennaCollection<TimeDomainAntenna> detectorZHS;
std::string name_center = "CoREAS_R=0_m--Phi=0degrees"; std::string name_center = "CoREAS_R=0_m--Phi=0degrees";
auto triggertime_center{((triggerpoint_ - center).getNorm() / constants::c) - 500_ns}; auto triggertime_center{((triggerpoint_ - center).getNorm() / constants::c) - 500_ns};
...@@ -127,7 +121,7 @@ int main() { ...@@ -127,7 +121,7 @@ int main() {
auto rr_1 = static_cast<int>(radius_1 / 1_m); auto rr_1 = static_cast<int>(radius_1 / 1_m);
auto const point_1{Point(rootCS, centerX + radius_1 * cos(phiRad_1), auto const point_1{Point(rootCS, centerX + radius_1 * cos(phiRad_1),
centerY + radius_1 * sin(phiRad_1), centerZ)}; 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}; auto triggertime_1{((triggerpoint_ - point_1).getNorm() / constants::c) - 500_ns};
std::string name_1 = "CoREAS_R=" + std::to_string(rr_1) + std::string name_1 = "CoREAS_R=" + std::to_string(rr_1) +
"_m--Phi=" + std::to_string(phi_1) + "degrees"; "_m--Phi=" + std::to_string(phi_1) + "degrees";
...@@ -186,16 +180,8 @@ int main() { ...@@ -186,16 +180,8 @@ int main() {
coreas(detectorCoREAS, SP); coreas(detectorCoREAS, SP);
output.add("CoREAS", coreas); 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); TimeCut cut(period / 4);
// TrackWriter trackWriter;
// output.add("tracks", trackWriter); // register TrackWriter
// assemble all processes into an ordered process list // assemble all processes into an ordered process list
auto sequence = make_sequence(coreas, cut); 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