IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 4814f83b authored by Nikos Karastathis's avatar Nikos Karastathis :ocean:
Browse files

small updates to em_shower.cpp and vertical_EAS.cpp

parent b07c2c87
No related branches found
No related tags found
1 merge request!437Resolve "Examples need some polishing"
Pipeline #8546 canceled
/* /*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * 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 * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
...@@ -55,9 +55,9 @@ ...@@ -55,9 +55,9 @@
NOTE, WARNING, ATTENTION NOTE, WARNING, ATTENTION
The .../Random.hpppp implement the hooks of external modules to the C8 random The .../Random.hpppp implement the hooks of external modules to the C8 random
number generator. It has to occur excatly ONCE per linked number generator. It has to occur exactly ONCE per linked
executable. If you include the header below multiple times and executable. If you include the header below multiple times and
link this togehter, it will fail. link this together, it will fail.
*/ */
#include <corsika/modules/Random.hpp> #include <corsika/modules/Random.hpp>
...@@ -84,12 +84,14 @@ int main(int argc, char** argv) { ...@@ -84,12 +84,14 @@ int main(int argc, char** argv) {
logging::set_level(logging::level::info); logging::set_level(logging::level::info);
if (argc != 3) { if (argc != 3) {
std::cerr << "usage: em_shower <energy/GeV> <theta>" << std::endl; std::cerr << "usage: em_shower <energy/GeV> [seed] - put 0 for random seed" << std::endl;
return 1; return 1;
} }
feenableexcept(FE_INVALID); feenableexcept(FE_INVALID);
int seed = 0;
if (argc > 2) { seed = std::stoi(std::string(argv[2])); }
// initialize random number sequence(s) // initialize random number sequence(s)
int seed = 2723141261;
registerRandomStreams(seed); registerRandomStreams(seed);
// setup environment, geometry // setup environment, geometry
...@@ -119,7 +121,7 @@ int main(int argc, char** argv) { ...@@ -119,7 +121,7 @@ int main(int argc, char** argv) {
const Code beamCode = Code::Electron; 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])); const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
double theta = 1. * std::stof(std::string(argv[2])); double theta = 0.;
auto const thetaRad = theta / 180. * M_PI; auto const thetaRad = theta / 180. * M_PI;
HEPMomentumType P0 = calculate_momentum(E0, mass); HEPMomentumType P0 = calculate_momentum(E0, mass);
...@@ -155,8 +157,7 @@ int main(int argc, char** argv) { ...@@ -155,8 +157,7 @@ int main(int argc, char** argv) {
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env, ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000}; false, 1000};
std::string outname_ {"em_shower_outputs"}; OutputManager output("em_shower_outputs");
OutputManager output(outname_ + std::to_string(std::stof(std::string(argv[1]))) + "theta-" + std::to_string(std::stof(std::string(argv[2]))));
EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200}; EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
// register energy losses as output // register energy losses as output
...@@ -174,8 +175,8 @@ int main(int argc, char** argv) { ...@@ -174,8 +175,8 @@ int main(int argc, char** argv) {
// 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);
// TrackWriter tracks; TrackWriter tracks;
// output.add("tracks", tracks); output.add("tracks", tracks);
// long. profile // long. profile
LongitudinalWriter profile{showerAxis, 10_g / square(1_cm)}; LongitudinalWriter profile{showerAxis, 10_g / square(1_cm)};
......
/* /*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
* *
* This software is distributed under the terms of the GNU General Public * 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 * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
...@@ -68,9 +68,9 @@ ...@@ -68,9 +68,9 @@
NOTE, WARNING, ATTENTION NOTE, WARNING, ATTENTION
The .../Random.hpppp implement the hooks of external modules to the C8 random The .../Random.hpppp implement the hooks of external modules to the C8 random
number generator. It has to occur excatly ONCE per linked number generator. It has to occur exactly ONCE per linked
executable. If you include the header below multiple times and executable. If you include the header below multiple times and
link this togehter, it will fail. link this together, it will fail.
*/ */
#include <corsika/modules/Random.hpp> #include <corsika/modules/Random.hpp>
...@@ -105,7 +105,7 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>; ...@@ -105,7 +105,7 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) { int main(int argc, char** argv) {
logging::set_level(logging::level::warn); logging::set_level(logging::level::info);
CORSIKA_LOG_INFO("vertical_EAS"); CORSIKA_LOG_INFO("vertical_EAS");
...@@ -135,13 +135,14 @@ int main(int argc, char** argv) { ...@@ -135,13 +135,14 @@ int main(int argc, char** argv) {
env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm, env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm,
MagneticFieldVector{rootCS, 20.4_uT, 0_T, 43.23_uT}); MagneticFieldVector{rootCS, 20.4_uT, 0_T, 43.23_uT});
std::unordered_map<Code, HEPEnergyType> energy_resolution = { // Uncomment if you want to use PROPOSAL
{Code::Electron, 2_MeV}, // std::unordered_map<Code, HEPEnergyType> energy_resolution = {
{Code::Positron, 2_MeV}, // {Code::Electron, 2_MeV},
{Code::Photon, 2_MeV}, // {Code::Positron, 2_MeV},
}; // {Code::Photon, 2_MeV},
for (auto [pcode, energy] : energy_resolution) // };
set_energy_production_threshold(pcode, energy); // for (auto [pcode, energy] : energy_resolution)
// set_energy_production_threshold(pcode, energy);
// pre-setup particle stack // pre-setup particle stack
unsigned short const A = std::stoi(std::string(argv[1])); unsigned short const A = std::stoi(std::string(argv[1]));
...@@ -216,7 +217,8 @@ int main(int argc, char** argv) { ...@@ -216,7 +217,8 @@ int main(int argc, char** argv) {
// construct the continuous energy loss model // construct the continuous energy loss model
BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX}; BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
// construct a particle cut // 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}; ParticleCut<SubWriter<decltype(dEdX)>> cut{2_MeV, 2_MeV, 2_GeV, 300_MeV, true, dEdX};
// setup longitudinal profile // setup longitudinal profile
...@@ -226,8 +228,8 @@ int main(int argc, char** argv) { ...@@ -226,8 +228,8 @@ int main(int argc, char** argv) {
LongitudinalProfile<SubWriter<decltype(longProf)>> profile{longProf}; LongitudinalProfile<SubWriter<decltype(longProf)>> profile{longProf};
// create a track writer and register it with the output manager // create a track writer and register it with the output manager
// TrackWriter<TrackWriterParquet> trackWriter; TrackWriter<TrackWriterParquet> trackWriter;
// output.add("tracks", trackWriter); output.add("tracks", trackWriter);
// setup processes, decays and interactions // setup processes, decays and interactions
...@@ -235,6 +237,7 @@ int main(int argc, char** argv) { ...@@ -235,6 +237,7 @@ int main(int argc, char** argv) {
InteractionCounter sibyllCounted{sibyll}; InteractionCounter sibyllCounted{sibyll};
HEPEnergyType heThresholdNN = 60_GeV; HEPEnergyType heThresholdNN = 60_GeV;
// PROPOSAL is disabled for this example
// corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(), heThresholdNN); // corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(), heThresholdNN);
// corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX); // corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment