first attempt of radio shower example

corsika/modules/radio/RadioProcess.hpp

@@ -83,7 +83,9 @@ namespace corsika {
       // important for controlling the runtime of radio (by ignoring particles
       // that aren't going to contribute i.e. heavy hadrons)
       //if (valid(particle, track)) {
-      if (particle.is_em) { return this->implementation().simulate(particle, track); }
+      if (particle == Code::Electron || particle == Code::Positron) {
+        return this->implementation().simulate(particle, track);
+      }
       //}
     }
 

examples/radio_shower.cpp

@@ -49,6 +49,15 @@
 #include <corsika/modules/UrQMD.hpp>
 #include <corsika/modules/PROPOSAL.hpp>
 
+#include <corsika/modules/radio/RadioProcess.hpp>
+#include <corsika/modules/radio/CoREAS.hpp>
+#include <corsika/modules/radio/antennas/TimeDomainAntenna.hpp>
+#include <corsika/modules/radio/detectors/RadioDetector.hpp>
+#include <corsika/modules/radio/propagators/StraightPropagator.hpp>
+#include <corsika/modules/radio/propagators/SignalPath.hpp>
+#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
+
+
 #include <corsika/setup/SetupStack.hpp>
 #include <corsika/setup/SetupTrajectory.hpp>
 
@@ -110,15 +119,198 @@ int main(int argc, char** argv) {
   // initialize random number sequence(s)
   registerRandomStreams(seed);
 
-  // setup environment (idea 2)
+  // setup environment
   using EnvType = setup::Environment;
   EnvType env;
   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)};
+  // 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);
+  // the detector
+  std::vector<TimeDomainAntenna> detector;
+  detector.push_back(ant1);
+  detector.push_back(ant2);
   auto builder = make_layered_spherical_atmosphere_builder<
       setup::EnvironmentInterface, MyExtraEnv>::create(center,
                                                        constants::EarthRadius::Mean, 1.,
                                                        Medium::AirDry1Atm,
                                                        MagneticFieldVector{rootCS, 0_T,
                                                                            50_uT, 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, 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);
+  builder.assemble(env);
+
+  // 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]));
+  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;
+    }
+  }
+
+  // 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
+
+  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);
+  // put radio here
+  CoREAS<TimeDomainAntenna, StraightPropagator(env)> coreas;
+
+  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);
+
+  // 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);
+
+  // 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();
+
+  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();
+
+  // get radio pulse
+  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");
+
 }