IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 2d21e488 authored by ralfulrich's avatar ralfulrich Committed by Ralf Ulrich
Browse files

cleanup examples

parent b5843a0a
No related branches found
No related tags found
1 merge request!352Better examples
Hide whitespace changes
Inline Side-by-side
Showing
with 94 additions and 212 deletions
corsika/detail/modules/proposal/ContinuousProcess.inl
+ 4
− 3
View file @ 2d21e488
/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
......@@ -132,9 +133,9 @@ namespace corsika::proposal {
std::max(energy * 0.9, // either 10% relative loss max., or
get_kinetic_energy_threshold(
code) // energy thresholds globally defined for individual particles
* 0.9999 // need to go slightly below global e-cut to assure removal
// in ParticleCut. This does not matter since at cut-time
// the entire energy is removed.
* 0.9999 // need to go slightly below global e-cut to assure removal in
// ParticleCut. This does not matter since at cut-time the
// entire energy is removed.
);
// solving the track integral for giving energy lim
......
corsika/detail/modules/proposal/Interaction.inl
+ 2
− 1
View file @ 2d21e488
......@@ -98,7 +98,8 @@ namespace corsika::proposal {
auto vec = QuantityVector(vecProposal.GetX() * E, vecProposal.GetY() * E,
vecProposal.GetZ() * E);
auto p = MomentumVector(labCS, vec);
auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
auto sec_code =
convert_from_PDG(static_cast<PDGCode>(s.type));
view.addSecondary(
std::make_tuple(sec_code, p, projectile.getPosition(), projectile.getTime()));
}
......
examples/CMakeLists.txt
+ 5
− 1
View file @ 2d21e488
......@@ -35,7 +35,7 @@ CORSIKA_REGISTER_EXAMPLE (cascade_proton_example)
add_executable (vertical_EAS vertical_EAS.cpp)
target_link_libraries (vertical_EAS CORSIKA8::CORSIKA8)
CORSIKA_REGISTER_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000. 1)
CORSIKA_REGISTER_EXAMPLE (vertical_EAS RUN_OPTIONS 4 2 10000.)
add_executable (stopping_power stopping_power.cpp)
target_link_libraries (stopping_power CORSIKA8::CORSIKA8)
......@@ -56,3 +56,7 @@ CORSIKA_REGISTER_EXAMPLE (em_shower RUN_OPTIONS "100.")
add_executable (hybrid_MC hybrid_MC.cpp)
target_link_libraries (hybrid_MC CORSIKA8::CORSIKA8)
CORSIKA_REGISTER_EXAMPLE (hybrid_MC RUN_OPTIONS 4 2 10000.)
add_executable (corsika corsika.cpp)
target_link_libraries (corsika CORSIKA8::CORSIKA8)
CORSIKA_REGISTER_EXAMPLE (corsika RUN_OPTIONS 4 2 10000. 1)
examples/boundary_example.cpp
+ 1
− 10
View file @ 2d21e488
......@@ -26,7 +26,6 @@
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/ParticleCut.hpp>
......@@ -38,9 +37,6 @@
executable. If you include the header below multiple times and
link this togehter, it will fail.
*/
#include <corsika/modules/sibyll/Random.hpp>
#include <corsika/modules/urqmd/Random.hpp>
#include <iostream>
#include <limits>
#include <typeinfo>
......@@ -83,7 +79,6 @@ private:
int main() {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CORSIKA_LOG_INFO("boundary_example");
......@@ -120,10 +115,6 @@ int main() {
// setup processes, decays and interactions
setup::Tracking tracking;
RNGManager::getInstance().registerRandomStream("sibyll");
corsika::sibyll::Interaction sibyll;
corsika::sibyll::Decay decay;
ParticleCut cut(50_GeV, true, true);
TrackWriter trackWriter;
......@@ -132,7 +123,7 @@ int main() {
MyBoundaryCrossingProcess<true> boundaryCrossing("crossings.dat");
// assemble all processes into an ordered process list
auto sequence = make_sequence(sibyll, decay, cut, boundaryCrossing, trackWriter);
auto sequence = make_sequence(cut, boundaryCrossing, trackWriter);
// setup particle stack, and add primary particles
setup::Stack stack;
......
examples/cascade_proton_example.cpp
+ 0
− 1
View file @ 2d21e488
......@@ -59,7 +59,6 @@ using namespace std;
int main() {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
std::cout << "cascade_proton_example" << std::endl;
......
examples/em_shower.cpp
+ 3
− 2
View file @ 2d21e488
......@@ -139,7 +139,8 @@ int main(int argc, char** argv) {
<< std::endl;
OutputManager output("em_shower_outputs");
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env};
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000.};
// setup processes, decays and interactions
......@@ -189,4 +190,4 @@ int main(int argc, char** argv) {
longprof.save("longprof_emShower.txt");
output.endOfLibrary();
}
}
\ No newline at end of file
examples/hybrid_MC.cpp
+ 2
− 1
View file @ 2d21e488
......@@ -173,7 +173,8 @@ int main(int argc, char** argv) {
<< std::endl;
OutputManager output("hybrid_MC_outputs");
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env};
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env, true,
1000};
// setup processes, decays and interactions
......
examples/stack_example.cpp
+ 2
− 2
View file @ 2d21e488
......@@ -44,11 +44,11 @@ void read(VectorStack& s) {
int main() {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
std::cout << "stack_example" << std::endl;
CORSIKA_LOG_INFO("stack_example");
VectorStack s;
fill(s);
read(s);
CORSIKA_LOG_INFO("done");
return 0;
}
examples/stopping_power.cpp
+ 7
− 9
View file @ 2d21e488
......@@ -30,7 +30,6 @@ using namespace std;
int main() {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CORSIKA_LOG_INFO("stopping_power");
......@@ -48,15 +47,16 @@ int main() {
rootCS, 0_m, 0_m,
112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
ShowerAxis showerAxis{injectionPos, Vector<length_d>{rootCS, 0_m, 0_m, 1_m}, env};
ShowerAxis showerAxis{injectionPos, Vector<length_d>{rootCS, 0_m, 0_m, 1_m}, env, false,
100};
BetheBlochPDG eLoss{showerAxis};
setup::Stack stack;
std::ofstream file("dEdX.dat");
file << "# beta*gamma, dE/dX / eV/(g/cm²)" << std::endl;
file << "# beta*gamma, dE/dX / MeV/(g/cm²)" << std::endl;
for (HEPEnergyType E0 = 300_MeV; E0 < 1_PeV; E0 *= 1.05) {
for (HEPEnergyType E0 = 200_MeV; E0 < 1_PeV; E0 *= 1.05) {
stack.clear();
const Code beamCode = Code::MuPlus;
const HEPMassType mass = get_mass(beamCode);
......@@ -72,16 +72,14 @@ int main() {
-ptot * cos(theta));
};
auto const [px, py, pz] =
momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0);
momentumComponents(theta / 180. * constants::pi, phi / 180. * constants::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, plab, injectionPos, 0_ns));
auto const p = stack.getNextParticle();
HEPEnergyType dE = eLoss.getTotalEnergyLoss(p, 1_g / square(1_cm));
file << P0 / mass << "\t" << -dE / 1_eV << std::endl;
file << P0 / mass << "\t" << -dE / 1_MeV << std::endl;
}
CORSIKA_LOG_INFO("finished writing dEdX.dat");
}
examples/vertical_EAS.cpp
+ 68
− 182
View file @ 2d21e488
......@@ -43,15 +43,12 @@
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/OnShellCheck.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/Pythia8.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/UrQMD.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/QGSJetII.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
......@@ -71,7 +68,6 @@
*/
#include <corsika/modules/sibyll/Random.hpp>
#include <corsika/modules/urqmd/Random.hpp>
#include <corsika/modules/qgsjetII/Random.hpp>
using namespace corsika;
using namespace std;
......@@ -80,7 +76,6 @@ using Particle = setup::Stack::particle_type;
void registerRandomStreams(const int seed) {
RNGManager::getInstance().registerRandomStream("cascade");
RNGManager::getInstance().registerRandomStream("qgsjet");
RNGManager::getInstance().registerRandomStream("sibyll");
RNGManager::getInstance().registerRandomStream("pythia");
RNGManager::getInstance().registerRandomStream("urqmd");
......@@ -105,8 +100,8 @@ int main(int argc, char** argv) {
CORSIKA_LOG_INFO("vertical_EAS");
if (argc < 5) {
std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> <Nevt> [seed] \n"
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;
......@@ -115,9 +110,8 @@ int main(int argc, char** argv) {
feenableexcept(FE_INVALID);
int seed = 0;
int number_showers = std::stoi(std::string(argv[4]));
if (argc > 5) { seed = std::stoi(std::string(argv[5])); }
if (argc > 4) { seed = std::stoi(std::string(argv[4])); }
// initialize random number sequence(s)
registerRandomStreams(seed);
......@@ -162,8 +156,8 @@ int main(int argc, char** argv) {
HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
double theta = 0.;
double phi = 180.;
auto const thetaRad = theta / 180. * M_PI;
auto const phiRad = phi / 180. * M_PI;
auto const thetaRad = theta / 180. * constants::pi;
auto const phiRad = phi / 180. * constants::pi;
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
......@@ -200,14 +194,14 @@ int main(int argc, char** argv) {
std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.5
<< std::endl;
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
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");
// setup processes, decays and interactions
// corsika::qgsjetII::Interaction qgsjet;
corsika::sibyll::Interaction sibyll;
InteractionCounter sibyllCounted(sibyll);
......@@ -240,15 +234,6 @@ int main(int argc, char** argv) {
decaySibyll.printDecayConfig();
ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
corsika::proposal::Interaction emCascade(env);
corsika::proposal::ContinuousProcess emContinuous(env);
InteractionCounter emCascadeCounted(emCascade);
OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
corsika::urqmd::UrQMD urqmd;
InteractionCounter urqmdCounted{urqmd};
......@@ -265,170 +250,71 @@ int main(int argc, char** argv) {
make_sequence(sibyllNucCounted, sibyllCounted));
auto decaySequence = make_sequence(decayPythia, decaySibyll);
for (int i_shower = 1; i_shower < number_showers + 1; i_shower++) {
// directory for outputs
string const labHist_file = "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz";
string const cMSHist_file = "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz";
string const longprof_file = "longprof_verticalEAS_" + to_string(i_shower) + ".txt";
std::cout << std::endl;
std::cout << "Shower " << i_shower << "/" << number_showers << std::endl;
// setup particle stack, and add primary particle
setup::Stack stack;
stack.clear();
unsigned short const A = std::stoi(std::string(argv[1]));
Code beamCode;
HEPEnergyType mass;
unsigned short Z = 0;
if (A > 0) {
beamCode = Code::Nucleus;
Z = std::stoi(std::string(argv[2]));
mass = get_nucleus_mass(A, Z);
} else {
int pdg = std::stoi(std::string(argv[2]));
beamCode = convert_from_PDG(PDGCode(pdg));
mass = get_mass(beamCode);
}
HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
double theta = 0.;
auto const thetaRad = theta / 180. * M_PI;
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
return sqrt((Elab - m) * (Elab + m));
};
HEPMomentumType P0 = elab2plab(E0, mass);
auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
};
auto const [px, py, pz] = momentumComponents(thetaRad, P0);
auto plab = 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 = 0_km + builder.getEarthRadius();
auto const injectionHeight = 111.75_km + builder.getEarthRadius();
auto const t = (injectionHeight - observationHeight) / cos(thetaRad);
Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
Point const injectionPos =
showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
if (A > 1) {
stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, A, Z));
// directory for outputs
string const labHist_file = "inthist_lab_verticalEAS.npz";
string const cMSHist_file = "inthist_cms_verticalEAS.npz";
string const longprof_file = "longprof_verticalEAS.txt";
} else {
if (A == 1) {
if (Z == 1) {
stack.addParticle(std::make_tuple(Code::Proton, plab, injectionPos, 0_ns));
} else if (Z == 0) {
stack.addParticle(std::make_tuple(Code::Neutron, plab, injectionPos, 0_ns));
} else {
std::cerr << "illegal parameters" << std::endl;
return EXIT_FAILURE;
}
// setup particle stack, and add primary particle
setup::Stack stack;
stack.clear();
if (A > 1) {
stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, A, Z));
} else {
if (A == 1) {
if (Z == 1) {
stack.addParticle(std::make_tuple(Code::Proton, plab, injectionPos, 0_ns));
} else if (Z == 0) {
stack.addParticle(std::make_tuple(Code::Neutron, plab, injectionPos, 0_ns));
} else {
stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
std::cerr << "illegal parameters" << std::endl;
return EXIT_FAILURE;
}
} else {
stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
}
}
// 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;
BetheBlochPDG emContinuous(showerAxis);
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env,
false};
// setup processes, decays and interactions
// corsika::qgsjetII::Interaction qgsjet;
corsika::sibyll::Interaction sibyll;
InteractionCounter sibyllCounted(sibyll);
corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
InteractionCounter sibyllNucCounted(sibyllNuc);
corsika::pythia8::Decay decayPythia;
// use sibyll decay routine for decays of particles unknown to pythia
corsika::sibyll::Decay decaySibyll{{
Code::N1440Plus,
Code::N1440MinusBar,
Code::N1440_0,
Code::N1440_0Bar,
Code::N1710Plus,
Code::N1710MinusBar,
Code::N1710_0,
Code::N1710_0Bar,
Code::Pi1300Plus,
Code::Pi1300Minus,
Code::Pi1300_0,
Code::KStar0_1430_0,
Code::KStar0_1430_0Bar,
Code::KStar0_1430_Plus,
Code::KStar0_1430_MinusBar,
}};
decaySibyll.printDecayConfig();
ParticleCut cut{60_GeV, true, true};
// corsika::proposal::Interaction emCascade(env);
// corsika::proposal::ContinuousProcess emContinuous(env);
// InteractionCounter emCascadeCounted(emCascade);
BetheBlochPDG emContinuous(showerAxis);
OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
TrackWriter trackWriter;
output.add("tracks", trackWriter); // register TrackWriter
LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
// register the observation plane with the output
output.add("particles", observationLevel);
auto sequence = make_sequence( // emCascadeCounted,
stackInspect, hadronSequence, reset_particle_mass, decaySequence,
// emContinuous,
BetheBlochPDG(showerAxis), cut, trackWriter, observationLevel, longprof);
// define air shower object, run simulation
setup::Tracking tracking;
Cascade EAS(env, tracking, sequence, output, stack);
// to fix the point of first interaction, uncomment the following two lines:
// EAS.forceInteraction();
EAS.run();
cut.showResults();
// emContinuous.showResults();
observationLevel.showResults();
const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
cut.getEmEnergy() + // emContinuous.getEnergyLost() +
observationLevel.getEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.reset();
cut.reset();
// emContinuous.reset();
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
save_hist(hists.labHist(), labHist_file, true);
save_hist(hists.CMSHist(), cMSHist_file, true);
longprof.save(longprof_file);
output.endOfLibrary();
}
TrackWriter trackWriter;
output.add("tracks", trackWriter); // register TrackWriter
LongitudinalProfile longprof{showerAxis};
Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
// register the observation plane with the output
output.add("particles", observationLevel);
auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence, emContinuous,
cut, trackWriter, observationLevel, longprof);
// define air shower object, run simulation
setup::Tracking tracking;
Cascade EAS(env, tracking, sequence, output, stack);
EAS.run();
cut.showResults();
// emContinuous.showResults();
observationLevel.showResults();
const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
cut.getEmEnergy() + // emContinuous.getEnergyLost() +
observationLevel.getEnergyGround();
cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
observationLevel.reset();
cut.reset();
// emContinuous.reset();
auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
urqmdCounted.getHistogram();
save_hist(hists.labHist(), labHist_file, true);
save_hist(hists.CMSHist(), cMSHist_file, true);
longprof.save(longprof_file);
output.endOfLibrary();
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum