IAP GITLAB

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

updated version of radio electron shower

parent 61c06a03
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
+ 164
− 191
View file @ d139637d
......@@ -72,10 +72,10 @@
using namespace corsika;
using namespace std;
void registerRandomStreams() {
RNGManager::getInstance().registerRandomStream("cascade");
RNGManager::getInstance().registerRandomStream("proposal");
RNGManager::getInstance().seedAll();
void registerRandomStreams(const int seed) {
RNGManager::getInstance().registerRandomStream("cascade");
RNGManager::getInstance().registerRandomStream("proposal");
RNGManager::getInstance().seedAll(seed);
}
template <typename TInterface>
......@@ -84,184 +84,157 @@ UniformRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
int main(int argc, char** argv) {
logging::set_level(logging::level::info);
if (argc != 2) {
std::cerr << "usage: em_shower <energy/GeV>" << std::endl;
return 1;
}
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
registerRandomStreams();
// setup environment, geometry
using EnvType = setup::Environment;
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
auto builder = make_layered_spherical_atmosphere_builder<
setup::EnvironmentInterface, MyExtraEnv>::create(center,
constants::EarthRadius::Mean, 1.000327,
Medium::AirDry1Atm,
MagneticFieldVector{rootCS, 50_uT,
0_T, 0_T});
builder.setNuclearComposition(
{{Code::Nitrogen, Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 2_km);
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
builder.addLinearLayer(1e9_cm, 112.8_km + constants::EarthRadius::Mean);
builder.assemble(env);
// setup particle stack, and add primary particle
setup::Stack 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;
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 = 1.4_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;
stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02
<< std::endl;
OutputManager output("radio_em_shower_outputs");
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env};
// 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_ = 100_m; radius_ <= 200_m; radius_ += 100_m) {
// for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
// auto phiRad_ = phi_ / 180. * M_PI;
// auto const point_ {Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_), showerCoreY_ + radius_ * sin(phiRad_), builder.getEarthRadius())};
// auto triggertime_ {(triggerpoint_ - point_).getNorm() / constants::c};
// const int rr_ = static_cast<int>(radius_ / 1_m);
// std::string name_ = "CoREAS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
// TimeDomainAntenna antenna_(name_, point_, triggertime_, duration_, sampleRate_);
// detectorCoREAS.addAntenna(antenna_);
// }
// }
//
// // setup ZHS antennas
// for (auto radius_ = 100_m; radius_ <= 200_m; radius_ += 100_m) {
// for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
// auto phiRad_ = phi_ / 180. * M_PI;
// auto const point_ {Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_), showerCoreY_ + radius_ * sin(phiRad_), builder.getEarthRadius())};
// auto triggertime_ {(triggerpoint_ - point_).getNorm() / constants::c};
// const int rr_ = static_cast<int>(radius_ / 1_m);
// std::string name_ = "ZHS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
// TimeDomainAntenna antenna_(name_, point_, triggertime_, duration_, sampleRate_);
// detectorZHS.addAntenna(antenna_);
// }
// }
// 4 dummy antennas for CoREAS
for (auto radius_ = 100_m; radius_ <= 200_m; radius_ += 100_m) {
for (auto phi_ = 0; phi_ <= 270; phi_ += 90) {
auto phiRad_ = phi_ / 180. * M_PI;
auto const point_{
Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_), showerCoreY_ + radius_ * sin(phiRad_),
builder.getEarthRadius())};
auto triggertime_{(triggerpoint_ - point_).getNorm() / constants::c};
const int rr_ = static_cast<int>(radius_ / 1_m);
std::string name_ = "CoREAS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
TimeDomainAntenna antenna_(name_, point_, triggertime_, duration_, sampleRate_);
detectorCoREAS.addAntenna(antenna_);
}
}
// 4 dummy antennas for ZHS
for (auto radius_ = 100_m; radius_ <= 200_m; radius_ += 100_m) {
for (auto phi_ = 0; phi_ <= 270; phi_ += 90) {
auto phiRad_ = phi_ / 180. * M_PI;
auto const point_{
Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_), showerCoreY_ + radius_ * sin(phiRad_),
builder.getEarthRadius())};
auto triggertime_{(triggerpoint_ - point_).getNorm() / constants::c};
const int rr_ = static_cast<int>(radius_ / 1_m);
std::string name_ = "ZHS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
TimeDomainAntenna antenna_(name_, point_, triggertime_, duration_, sampleRate_);
detectorZHS.addAntenna(antenna_);
}
}
// setup processes, decays and interactions
ParticleCut cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true);
corsika::proposal::Interaction emCascade(env);
corsika::proposal::ContinuousProcess emContinuous(env);
InteractionCounter emCascadeCounted(emCascade);
logging::set_level(logging::level::info);
if (argc != 4) {
std::cerr << "usage: radio_em_shower <energy/GeV> <concentric ring number> <seed>" << std::endl;
return 1;
}
int seed {std::stof(std::string(argv[3]))};
std::cout << "Seed: " << seed << std::endl;
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
registerRandomStreams(seed);
// setup environment, geometry
using EnvType = setup::Environment;
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
auto builder = make_layered_spherical_atmosphere_builder<
setup::EnvironmentInterface, MyExtraEnv>::create(center,
constants::EarthRadius::Mean, 1.000327,
Medium::AirDry1Atm,
MagneticFieldVector{rootCS, 50_uT,
0_T, 0_T});
builder.setNuclearComposition(
{{Code::Nitrogen, Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 2_km);
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
builder.addLinearLayer(1e9_cm, 112.8_km + constants::EarthRadius::Mean);
builder.assemble(env);
// setup particle stack, and add primary particle
setup::Stack 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;
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 = 1.4_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;
stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02
<< std::endl;
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_);
ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env};
// 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 phi_ = 0; phi_ <= 315; phi_ += 45) {
auto phiRad_ = phi_ / 180. * M_PI;
auto const point_ {Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_), showerCoreY_ + radius_ * sin(phiRad_), builder.getEarthRadius())};
auto triggertime_ {(triggerpoint_ - point_).getNorm() / constants::c};
std::string name_ = "CoREAS_R=" + std::to_string(rr_) + "_m--Phi=" + std::to_string(phi_) + "degrees";
TimeDomainAntenna antenna_(name_, point_, triggertime_, duration_, sampleRate_);
detectorCoREAS.addAntenna(antenna_);
}
// setup ZHS antennas
for (auto phi_ = 0; phi_ <= 315; phi_ += 45) {
auto phiRad_ = phi_ / 180. * M_PI;
auto const point_ {Point(rootCS, showerCoreX_ + radius_ * cos(phiRad_), showerCoreY_ + radius_ * sin(phiRad_), builder.getEarthRadius())};
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_, triggertime_, duration_, sampleRate_);
detectorZHS.addAntenna(antenna_);
}
// setup processes, decays and interactions
ParticleCut cut(5_MeV, 5_MeV, 100_GeV, 100_GeV, true);
corsika::proposal::Interaction emCascade(env);
corsika::proposal::ContinuousProcess emContinuous(env);
InteractionCounter emCascadeCounted(emCascade);
// TrackWriter trackWriter;
// output.add("tracks", trackWriter); // register TrackWriter
// long. profile; columns for photon, e+, e- still need to be added
LongitudinalProfile longprof{showerAxis};
// long. profile; columns for photon, e+, e- still need to be added
LongitudinalProfile longprof{showerAxis};
// initiate CoREAS
RadioProcess<decltype(detectorCoREAS), CoREAS<decltype(detectorCoREAS),
decltype(SimplePropagator(env))>, decltype(SimplePropagator(env))>
coreas(detectorCoREAS, env);
// 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);
// 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);
// 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);
// register ZHS with the output manager
output.add("ZHS", zhs);
// Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
......@@ -269,20 +242,20 @@ int main(int argc, char** argv) {
// "particles.dat");
// output.add("obsplane", observationLevel);
auto sequence = make_sequence(emCascadeCounted, emContinuous, cut, longprof, coreas, zhs);
// ,observationLevel, trackWriter);
// define air shower object, run simulation
setup::Tracking tracking;
Cascade EAS(env, tracking, sequence, output, stack);
auto sequence = make_sequence(emCascadeCounted, emContinuous, cut, longprof, coreas, zhs);
// ,observationLevel, trackWriter);
// 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.setNodes();
// EAS.forceInteraction();
// to fix the point of first interaction, uncomment the following two lines:
// EAS.setNodes();
// EAS.forceInteraction();
EAS.run();
EAS.run();
cut.showResults();
emContinuous.showResults();
cut.showResults();
emContinuous.showResults();
// observationLevel.showResults();
// const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
// cut.getEmEnergy() + emContinuous.getEnergyLost() +
......@@ -290,13 +263,13 @@ int main(int argc, char** argv) {
// cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
// << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
// observationLevel.reset();
cut.reset();
emContinuous.reset();
cut.reset();
emContinuous.reset();
auto const hists = emCascadeCounted.getHistogram();
save_hist(hists.labHist(), "inthist_lab_emShower.npz", true);
save_hist(hists.CMSHist(), "inthist_cms_emShower.npz", true);
longprof.save("longprof_emShower.txt");
auto const hists = emCascadeCounted.getHistogram();
save_hist(hists.labHist(), "inthist_lab_emShower.npz", true);
save_hist(hists.CMSHist(), "inthist_cms_emShower.npz", true);
longprof.save("longprof_emShower.txt");
output.endOfLibrary();
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