IAP GITLAB

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

update radio_em_shower.cpp after rebasing

parent a6b865b5
No related branches found
No related tags found
1 merge request!329Radio interface
Hide whitespace changes
Inline Side-by-side
examples/radio_em_shower.cpp
+ 241
− 236
View file @ 7c3d873d
/*
* (c) Copyright 2020 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.
*/
* (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.
*/
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
......@@ -39,6 +39,7 @@
#include <corsika/modules/ObservationPlane.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/PROPOSAL.hpp>
#include <corsika/modules/radio/RadioProcess.hpp>
......@@ -62,248 +63,252 @@
#include <typeinfo>
/*
NOTE, WARNING, ATTENTION
NOTE, WARNING, ATTENTION
The .../Random.hpppp implement the hooks of external modules to the C8 random
number generator. It has to occur excatly ONCE per linked
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>
The .../Random.hpppp implement the hooks of external modules to the C8 random
number generator. It has to occur excatly ONCE per linked
executable. If you include the header below multiple times and
link this togehter, it will fail.
*/
#include <corsika/modules/Random.hpp>
using namespace corsika;
using namespace std;
void registerRandomStreams(int seed) {
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("proposal");
if (seed == 0) {
std::random_device rd;
seed = rd();
cout << "new random seed (auto) " << seed << endl;
}
RNGManager<>::getInstance().setSeed(seed);
RNGManager<>::getInstance().registerRandomStream("cascade");
RNGManager<>::getInstance().registerRandomStream("proposal");
RNGManager<>::getInstance().registerRandomStream("sibyll");
if (seed == 0) {
std::random_device rd;
seed = rd();
cout << "new random seed (auto) " << seed << endl;
}
RNGManager<>::getInstance().setSeed(seed);
}
template <typename TInterface>
using MyExtraEnv =
UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
int main(int argc, char** argv) {
logging::set_level(logging::level::info);
if (argc != 3) {
std::cerr
<< "usage: radio_em_shower <energy/GeV> <seed> - put seed=0 to use random seed"
<< std::endl;
return 1;
}
int seed{static_cast<int>(std::stof(std::string(argv[2])))};
std::cout << "Seed: " << seed << std::endl;
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
registerRandomStreams(seed);
// setup environment, geometry
using EnvironmentInterface =
IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
using EnvType = Environment<EnvironmentInterface>;
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
env, AtmosphereId::LinsleyUSStd, center, 1.000327, Medium::AirDry1Atm,
MagneticFieldVector{rootCS, 50_uT, 0_T, 0_T});
std::unordered_map<Code, HEPEnergyType> energy_resolution = {
{Code::Electron, 10_MeV},
{Code::Positron, 10_MeV},
{Code::Photon, 10_MeV},
};
for (auto [pcode, energy] : energy_resolution)
set_energy_production_threshold(pcode, energy);
// setup particle stack, and add primary particle
setup::Stack<EnvType> stack;
stack.clear();
const Code beamCode = Code::Electron;
auto const mass = get_mass(beamCode);
const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
double theta = 0.;
auto const thetaRad = theta / 180. * M_PI;
HEPMomentumType P0 = calculate_momentum(E0, mass);
auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
};
auto const [px, py, pz] = momentumComponents(thetaRad, P0);
auto plab = MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << endl;
cout << "input momentum: " << plab.getComponents() / 1_GeV
<< ", norm = " << plab.getNorm() << endl;
auto const observationHeight = 1.4_km + constants::EarthRadius::Mean;
auto const injectionHeight = 112.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), 0, cos(thetaRad)}} * t;
std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
stack.addParticle(std::make_tuple(
beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
plab.normalized(), injectionPos, 0_ns));
CORSIKA_LOG_INFO("shower axis length: {} ",
(showerCore - injectionPos).getNorm() * 1.02);
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000};
TimeType const groundHitTime{(showerCore - injectionPos).getNorm() / constants::c};
// int ring_number {std::stof(std::string(argv[2]))};
// std::cout << "Ring number : " << ring_number << std::endl;
// auto const radius_ {ring_number * 25_m};
// std::cout << "Radius = " << radius_ << std::endl;
// const int rr_ = static_cast<int>(radius_ / 1_m);
std::string outname_ = "radio_em_shower_outputs"; // + std::to_string(rr_);
OutputManager output(outname_);
// Radio objects
// the antenna time variables
const TimeType duration_{1e-6_s};
const InverseTimeType sampleRate_{1e+9_Hz};
// the detector (aka antenna collection) for CoREAS and ZHS
AntennaCollection<TimeDomainAntenna> detectorCoREAS;
AntennaCollection<TimeDomainAntenna> detectorZHS;
auto const showerCoreX_{showerCore.getCoordinates().getX()};
auto const showerCoreY_{showerCore.getCoordinates().getY()};
auto const injectionPosX_{injectionPos.getCoordinates().getX()};
auto const injectionPosY_{injectionPos.getCoordinates().getY()};
auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
// // setup CoREAS antennas
for (auto radius_1 = 25_m; radius_1 <= 500_m; radius_1 += 25_m) {
for (auto phi_1 = 0; phi_1 <= 315; phi_1 += 45) {
// auto radius_1 = 200_m;
// auto phi_1 = 45;
auto phiRad_1 = phi_1 / 180. * M_PI;
auto rr_1 = static_cast<int>(radius_1 / 1_m);
auto const point_1{Point(rootCS, showerCoreX_ + radius_1 * cos(phiRad_1),
showerCoreY_ + radius_1 * sin(phiRad_1),
constants::EarthRadius::Mean)};
std::cout << "Antenna point: " << point_1 << std::endl;
auto triggertime_1{(triggerpoint_ - point_1).getNorm() / constants::c};
std::string name_1 = "CoREAS_R=" + std::to_string(rr_1) +
"_m--Phi=" + std::to_string(phi_1) + "degrees";
TimeDomainAntenna antenna_1(name_1, point_1, rootCS, triggertime_1, duration_, sampleRate_,
triggertime_1);
detectorCoREAS.addAntenna(antenna_1);
}
}
// primary particle times -> t ground
// setup ZHS antennas
for (auto radius_ = 25_m; radius_ <= 500_m; radius_ += 25_m) {
for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
// auto radius_ = 200_m;
// auto phi_ = 45;
auto phiRad_ = phi_ / 180. * M_PI;
auto rr_ = static_cast<int>(radius_ / 1_m);
auto const point_{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_),
showerCoreY_ + radius_ * sin(phiRad_),
logging::set_level(logging::level::info);
if (argc != 3) {
std::cerr
<< "usage: radio_em_shower <energy/GeV> <seed> - put seed=0 to use random seed"
<< std::endl;
return 1;
}
int seed{static_cast<int>(std::stof(std::string(argv[2])))};
std::cout << "Seed: " << seed << std::endl;
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
registerRandomStreams(seed);
// setup environment, geometry
using EnvironmentInterface =
IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
using EnvType = Environment<EnvironmentInterface>;
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
env, AtmosphereId::LinsleyUSStd, center, 1.000327, Medium::AirDry1Atm,
MagneticFieldVector{rootCS, 50_uT, 0_T, 0_T});
std::unordered_map<Code, HEPEnergyType> energy_resolution = {
{Code::Electron, 5_MeV},
{Code::Positron, 5_MeV},
{Code::Photon, 5_MeV},
};
for (auto [pcode, energy] : energy_resolution)
set_energy_production_threshold(pcode, energy);
// setup particle stack, and add primary particle
setup::Stack<EnvType> stack;
stack.clear();
const Code beamCode = Code::Electron;
auto const mass = get_mass(beamCode);
const HEPEnergyType E0 = 1_GeV * std::stof(std::string(argv[1]));
double theta = 0.;
auto const thetaRad = theta / 180. * M_PI;
HEPMomentumType P0 = calculate_momentum(E0, mass);
auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
};
auto const [px, py, pz] = momentumComponents(thetaRad, P0);
auto plab = MomentumVector(rootCS, {px, py, pz});
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << endl;
cout << "input momentum: " << plab.getComponents() / 1_GeV
<< ", norm = " << plab.getNorm() << endl;
auto const observationHeight = 1.4_km + constants::EarthRadius::Mean;
auto const injectionHeight = 112.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), 0, cos(thetaRad)}} * t;
std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
stack.addParticle(std::make_tuple(
beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
plab.normalized(), injectionPos, 0_ns));
CORSIKA_LOG_INFO("shower axis length: {} ",
(showerCore - injectionPos).getNorm() * 1.02);
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env,
false, 1000};
TimeType const groundHitTime{(showerCore - injectionPos).getNorm() / constants::c};
// int ring_number {std::stof(std::string(argv[2]))};
// std::cout << "Ring number : " << ring_number << std::endl;
// auto const radius_ {ring_number * 25_m};
// std::cout << "Radius = " << radius_ << std::endl;
// const int rr_ = static_cast<int>(radius_ / 1_m);
std::string outname_ = "radio_em_shower_outputs"; // + std::to_string(rr_);
OutputManager output(outname_);
// Radio objects
// the antenna time variables
const TimeType duration_{1e-6_s};
const InverseTimeType sampleRate_{1e+9_Hz};
// the detector (aka antenna collection) for CoREAS and ZHS
AntennaCollection<TimeDomainAntenna> detectorCoREAS;
AntennaCollection<TimeDomainAntenna> detectorZHS;
auto const showerCoreX_{showerCore.getCoordinates().getX()};
auto const showerCoreY_{showerCore.getCoordinates().getY()};
auto const injectionPosX_{injectionPos.getCoordinates().getX()};
auto const injectionPosY_{injectionPos.getCoordinates().getY()};
auto const injectionPosZ_{injectionPos.getCoordinates().getZ()};
auto const triggerpoint_{Point(rootCS, injectionPosX_, injectionPosY_, injectionPosZ_)};
std::cout << "Trigger Point is: " << triggerpoint_ << std::endl;
// // setup CoREAS antennas
for (auto radius_1 = 25_m; radius_1 <= 500_m; radius_1 += 25_m) {
for (auto phi_1 = 0; phi_1 <= 315; phi_1 += 45) {
// auto radius_1 = 200_m;
// auto phi_1 = 45;
auto phiRad_1 = phi_1 / 180. * M_PI;
auto rr_1 = static_cast<int>(radius_1 / 1_m);
auto const point_1{Point(rootCS, showerCoreX_ + radius_1 * cos(phiRad_1),
showerCoreY_ + radius_1 * sin(phiRad_1),
constants::EarthRadius::Mean)};
auto triggertime_{(triggerpoint_ - point_).getNorm() / constants::c};
std::string name_ =
"ZHS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
TimeDomainAntenna antenna_(name_, point_, rootCS, triggertime_, duration_, sampleRate_,
triggertime_);
detectorZHS.addAntenna(antenna_);
}
}
// ----------------------- Radio objects
// --------------------------------------------------------------------
// setup processes, decays and interactions
EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
// register energy losses as output
output.add("dEdX", dEdX);
ParticleCut<SubWriter<decltype(dEdX)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true, dEdX);
corsika::proposal::Interaction emCascade(env);
corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
// BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
// NOT possible right now, due to interface differenc in PROPOSAL
// InteractionCounter emCascadeCounted(emCascade);
TrackWriter tracks;
output.add("tracks", tracks);
// long. profile
LongitudinalWriter profile{showerAxis, 10_g / square(1_cm), 200};
output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
// initiate CoREAS
RadioProcess<decltype(detectorCoREAS),
CoREAS<decltype(detectorCoREAS), decltype(SimplePropagator(env))>,
decltype(SimplePropagator(env))>
coreas(detectorCoREAS, env);
// register CoREAS with the output manager
output.add("CoREAS", coreas);
// initiate ZHS
RadioProcess<decltype(detectorZHS),
ZHS<decltype(detectorZHS), decltype(SimplePropagator(env))>,
decltype(SimplePropagator(env))>
zhs(detectorZHS, env);
// // register ZHS with the output manager
output.add("ZHS", zhs);
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(emCascade, emContinuous, longprof, cut, coreas, zhs);
auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, coreas, zhs,
observationLevel, tracks);
// define air shower object, run simulation
setup::Tracking tracking;
output.startOfLibrary();
Cascade EAS(env, tracking, sequence, output, stack);
// to fix the point of first interaction, uncomment the following two lines:
// EAS.forceInteraction();
EAS.run();
HEPEnergyType const Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
CORSIKA_LOG_INFO(
"total energy budget (GeV): {}, "
"relative difference (%): {}",
Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
output.endOfLibrary();
std::cout << "Antenna point: " << point_1 << std::endl;
auto triggertime_1{(triggerpoint_ - point_1).getNorm() / constants::c};
std::string name_1 = "CoREAS_R=" + std::to_string(rr_1) +
"_m--Phi=" + std::to_string(phi_1) + "degrees";
TimeDomainAntenna antenna_1(name_1, point_1, rootCS, triggertime_1, duration_, sampleRate_,
triggertime_1);
detectorCoREAS.addAntenna(antenna_1);
}
}
// primary particle times -> t ground
// setup ZHS antennas
for (auto radius_ = 25_m; radius_ <= 500_m; radius_ += 25_m) {
for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
// auto radius_ = 200_m;
// auto phi_ = 45;
auto phiRad_ = phi_ / 180. * M_PI;
auto rr_ = static_cast<int>(radius_ / 1_m);
auto const point_{Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_),
showerCoreY_ + radius_ * sin(phiRad_),
constants::EarthRadius::Mean)};
auto triggertime_{(triggerpoint_ - point_).getNorm() / constants::c};
std::string name_ =
"ZHS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
TimeDomainAntenna antenna_(name_, point_, rootCS, triggertime_, duration_, sampleRate_,
triggertime_);
detectorZHS.addAntenna(antenna_);
}
}
// ----------------------- Radio objects
// --------------------------------------------------------------------
// setup processes, decays and interactions
EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200};
// register energy losses as output
output.add("dEdX", dEdX);
ParticleCut<SubWriter<decltype(dEdX)>> cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true, dEdX);
corsika::sibyll::Interaction sibyll{env};
HEPEnergyType heThresholdNN = 80_GeV;
corsika::proposal::Interaction emCascade(env, sibyll.getHadronInteractionModel(),
heThresholdNN);
corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
// BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
// NOT possible right now, due to interface differenc in PROPOSAL
// InteractionCounter emCascadeCounted(emCascade);
TrackWriter tracks;
output.add("tracks", tracks);
// long. profile
LongitudinalWriter profile{showerAxis, 10_g / square(1_cm), 200};
output.add("profile", profile);
LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile};
// initiate CoREAS
RadioProcess<decltype(detectorCoREAS),
CoREAS<decltype(detectorCoREAS), decltype(SimplePropagator(env))>,
decltype(SimplePropagator(env))>
coreas(detectorCoREAS, env);
// register CoREAS with the output manager
output.add("CoREAS", coreas);
// initiate ZHS
RadioProcess<decltype(detectorZHS),
ZHS<decltype(detectorZHS), decltype(SimplePropagator(env))>,
decltype(SimplePropagator(env))>
zhs(detectorZHS, env);
// // register ZHS with the output manager
output.add("ZHS", zhs);
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(emCascade, emContinuous, longprof, cut, coreas, zhs);
auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, coreas, zhs,
observationLevel, tracks);
// define air shower object, run simulation
setup::Tracking tracking;
output.startOfLibrary();
Cascade EAS(env, tracking, sequence, output, stack);
// to fix the point of first interaction, uncomment the following two lines:
// EAS.forceInteraction();
EAS.run();
HEPEnergyType const Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround();
CORSIKA_LOG_INFO(
"total energy budget (GeV): {}, "
"relative difference (%): {}",
Efinal / 1_GeV, (Efinal / E0 - 1) * 100);
output.endOfLibrary();
}
\ No newline at end of file
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