From def128ebebbc23607b85c03795b4e840e5c7b714 Mon Sep 17 00:00:00 2001
From: Nikos Karastathis <n.karastathis@kit.edu>
Date: Mon, 22 Mar 2021 16:46:11 +0100
Subject: [PATCH] first try of synchnotron radiation example (need to tweak
StraightPropagator.hpp first)
---
examples/radio_shower2.cpp | 330 +++++++++++--------------------------
1 file changed, 96 insertions(+), 234 deletions(-)
diff --git a/examples/radio_shower2.cpp b/examples/radio_shower2.cpp
index e7620cc88..e7c3e0398 100644
--- a/examples/radio_shower2.cpp
+++ b/examples/radio_shower2.cpp
@@ -6,48 +6,25 @@
* the license.
*/
-/* clang-format off */
-// InteractionCounter used boost/histogram, which
-// fails if boost/type_traits have been included before. Thus, we have
-// to include it first...
-#include <corsika/framework/process/InteractionCounter.hpp>
-/* clang-format on */
-#include <corsika/framework/geometry/Plane.hpp>
-#include <corsika/framework/geometry/Sphere.hpp>
-#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/process/ProcessSequence.hpp>
-#include <corsika/framework/process/SwitchProcessSequence.hpp>
-#include <corsika/framework/process/InteractionCounter.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+
#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <corsika/framework/core/Cascade.hpp>
-#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/core/Logging.hpp>
#include <corsika/media/Environment.hpp>
-#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
-#include <corsika/media/IMagneticFieldModel.hpp>
-#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
-#include <corsika/media/UniformRefractiveIndex.hpp>
-#include <corsika/media/ShowerAxis.hpp>
-#include <corsika/media/SlidingPlanarExponential.hpp>
-#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/setup/SetupEnvironment.hpp>
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
#include <corsika/modules/radio/RadioProcess.hpp>
#include <corsika/modules/radio/CoREAS.hpp>
@@ -58,14 +35,13 @@
#include <corsika/modules/radio/propagators/SignalPath.hpp>
#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
-
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-#include <iomanip>
-#include <iostream>
-#include <limits>
-#include <string>
+#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>
/*
NOTE, WARNING, ATTENTION
@@ -78,241 +54,127 @@
#include <corsika/modules/sibyll/Random.hpp>
#include <corsika/modules/urqmd/Random.hpp>
+#include <iostream>
+#include <limits>
+#include <typeinfo>
+
using namespace corsika;
using namespace std;
-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");
- RNGManager::getInstance().registerRandomStream("proposal");
-
- if (seed == 0)
- RNGManager::getInstance().seedAll();
- else
- RNGManager::getInstance().seedAll(seed);
-}
-
+//
+// The example main program for a particle cascade
+//
+int main() {
-int main(int argc, char** argv) {
+ logging::set_level(logging::level::info);
+ corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
- corsika_logger->set_pattern("[%n:%^%-8l%$] %s:%#: %v");
- logging::set_level(logging::level::trace);
+ std::cout << "cascade_proton_example" << std::endl;
- CORSIKA_LOG_INFO("vertical_EAS");
-
- if (argc < 4) {
- std::cerr << "usage: vertical_EAS <A> <Z> <energy/GeV> [seed]" << std::endl;
- std::cerr << " if no seed is given, a random seed is chosen" << std::endl;
- return 1;
- }
feenableexcept(FE_INVALID);
-
- int seed = 0;
- if (argc > 4) seed = std::stoi(std::string(argv[4]));
// initialize random number sequence(s)
- registerRandomStreams(seed);
+ RNGManager::getInstance().registerRandomStream("cascade");
- // setup environment
- using IModelInterface = IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
- using AtmModel = UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<HomogeneousMedium
- <IModelInterface>>>>;
- using EnvType = Environment<AtmModel>;
+ // setup environment, geometry
+ using EnvType = setup::Environment;
EnvType env;
+ auto& universe = *(env.getUniverse());
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
- Point const center{rootCS, 0_m, 0_m, 0_m};
+
// the antenna location
- const auto point1{Point(env.getCoordinateSystem(), 50_m, 50_m, 50_m)};
- const auto point2{Point(env.getCoordinateSystem(), 25_m, 25_m, 25_m)};
+ const auto point1{Point(rootCS, 50_m, 50_m, 0_m)};
+ const auto point2{Point(rootCS, 50_m, -50_m, 0_m)};
+ const auto point3{Point(rootCS, -50_m, 50_m, 0_m)};
+ const auto point4{Point(rootCS, -50_m, -50_m, 0_m)};
+
// the antennas
- TimeDomainAntenna ant1("antenna1", point1, 0_s, 100_s, 1/1e-8_s);
- TimeDomainAntenna ant2("antenna2", point2, 0_s, 100_s, 1/1e-8_s);
+ TimeDomainAntenna ant1("antenna1", point1, 0_s, 1_s, 1/1e+6_s);
+ TimeDomainAntenna ant2("antenna2", point2, 0_s, 1_s, 1/1e+6_s);
+ TimeDomainAntenna ant3("antenna3", point3, 0_s, 1_s, 1/1e+6_s);
+ TimeDomainAntenna ant4("antenna4", point4, 0_s, 1_s, 1/1e+6_s);
+
// the detector
AntennaCollection<TimeDomainAntenna> detector;
detector.addAntenna(ant1);
detector.addAntenna(ant2);
- // a refractive index
- const double ri_{1.000327};
- // the constant density
- const auto density{19.2_g / cube(1_cm)};
- // the composition we use for the homogeneous medium
- NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
- std::vector<float>{1.f});
- // create magnetic field vector
- Vector B0(rootCS, 0_T,50_uT, 0_T);
- // create the medium
- auto Medium = EnvType::createNode<Sphere>(
- center, 1_km * std::numeric_limits<double>::infinity());
- // set the properties
- auto const props = Medium->setModelProperties<AtmModel>(ri_,
- Medium::AirDry1Atm, B0, density, protonComposition);
- env.getUniverse()->addChild(std::move(Medium));
+ detector.addAntenna(ant3);
+ detector.addAntenna(ant4);
+
+ auto world = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
+
+ using MyHomogeneousModel = MediumPropertyModel<
+ UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>;
+
+ world->setModelProperties<MyHomogeneousModel>(
+ Medium::AirDry1Atm, MagneticFieldVector(rootCS, 0_T, 0_T, 1_T),
+ 1_kg / (1_m * 1_m * 1_m),
+ NuclearComposition(std::vector<Code>{Code::Hydrogen},
+ std::vector<float>{(float)1.}));
+
+ universe.addChild(std::move(world));
// setup particle stack, and add primary particle
setup::Stack stack;
stack.clear();
- const Code beamCode = Code::Nucleus;
- unsigned short const A = std::stoi(std::string(argv[1]));
- unsigned short Z = std::stoi(std::string(argv[2]));
- auto const mass = get_nucleus_mass(A, Z);
- const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[3]));
+ const Code beamCode = Code::Electron;
+ const HEPMassType mass = Electron::mass;
+ const HEPEnergyType E0 = 1000_GeV;
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 = 112.75_km + builder.getEarthRadius();
- auto const t = -observationHeight * cos(thetaRad) +
- sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
- static_pow<2>(injectionHeight));
- Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
- Point const injectionPos =
- showerCore + 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, E0, plab, injectionPos, 0_ns, A, Z));
-
- } else {
- if (Z == 1) {
- stack.addParticle(std::make_tuple(Code::Proton, E0, plab, injectionPos, 0_ns));
- } else if (Z == 0) {
- stack.addParticle(std::make_tuple(Code::Neutron, E0, plab, injectionPos, 0_ns));
- } else {
- std::cerr << "illegal parameters" << std::endl;
- return EXIT_FAILURE;
- }
+ 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, E0, 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;
-
- ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.5, env};
-
// setup processes, decays and interactions
+ setup::Tracking tracking;
+ StackInspector<setup::Stack> stackInspect(1, true, E0);
- 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, 60_GeV, 60_GeV, 60_GeV, true};
- corsika::proposal::Interaction emCascade(env);
- corsika::proposal::ContinuousProcess emContinuous(env);
- InteractionCounter emCascadeCounted(emCascade);
+
+// ParticleCut cut(60_GeV, true, true);
// put radio process here
RadioProcess<decltype(detector), CoREAS<decltype(detector),
decltype(StraightPropagator(env))>, decltype(StraightPropagator(env))>
coreas(detector, env);
- OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
- TrackWriter trackWriter("tracks.dat");
-
- LongitudinalProfile longprof{showerAxis};
- Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
- ObservationPlane observationLevel(obsPlane, DirectionVector(rootCS, {1., 0., 0.}),
- "particles.dat");
- corsika::urqmd::UrQMD urqmd;
- InteractionCounter urqmdCounted{urqmd};
- StackInspector<setup::Stack> stackInspect(1000, false, E0);
+ TrackWriter trackWriter("tracks.dat");
+ ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
+// BetheBlochPDG eLoss{showerAxis};
// assemble all processes into an ordered process list
- struct EnergySwitch {
- HEPEnergyType cutE_;
- EnergySwitch(HEPEnergyType cutE)
- : cutE_(cutE) {}
- SwitchResult operator()(const Particle& p) {
- if (p.getEnergy() < cutE_)
- return SwitchResult::First;
- else
- return SwitchResult::Second;
- }
- };
- auto hadronSequence = make_select(
- urqmdCounted, make_sequence(sibyllNucCounted, sibyllCounted), EnergySwitch(55_GeV));
- auto decaySequence = make_sequence(decayPythia, decaySibyll);
- auto sequence =
- make_sequence(stackInspect, hadronSequence, reset_particle_mass, decaySequence,
- emContinuous, cut, coreas, trackWriter, observationLevel, longprof);
+ // auto sequence = make_sequence(sibyll, sibyllNuc, decay, eLoss, cut, trackWriter,
+ // stackInspect); auto sequence = make_sequence(sibyll, decay, eLoss, cut, trackWriter,
+ // stackInspect);
+ auto sequence = make_sequence(coreas, trackWriter, stackInspect);
// define air shower object, run simulation
- setup::Tracking tracking;
Cascade EAS(env, tracking, sequence, stack);
-
- // to fix the point of first interaction, uncomment the following two lines:
- // EAS.forceInteraction();
-
EAS.run();
+ cout << "Result: E0=" << E0 / 1_GeV << endl;
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();
+ const HEPEnergyType Efinal =
+ cut.getCutEnergy() + cut.getInvEnergy() + cut.getEmEnergy();
+ cout << "total energy (GeV): " << Efinal / 1_GeV << endl
+ << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
// get radio output
coreas.writeOutput();
-
- auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
- urqmdCounted.getHistogram();
-
- save_hist(hists.labHist(), "inthist_lab_verticalEAS.npz", true);
- save_hist(hists.CMSHist(), "inthist_cms_verticalEAS.npz", true);
- longprof.save("longprof_verticalEAS.txt");
-
-}
\ No newline at end of file
+}
--
GitLab