mars.cpp

/*
 * (c) Copyright 2018 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.
 */

/* 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/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/utility/CorsikaFenv.hpp>
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/geometry/PhysicalGeometry.hpp>

#include <corsika/output/OutputManager.hpp>
#include <corsika/output/NoOutput.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/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.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/modules/QGSJetII.hpp>

#include <corsika/setup/SetupStack.hpp>
#include <corsika/setup/SetupTrajectory.hpp>

#include <CLI/App.hpp>
#include <CLI/Formatter.hpp>
#include <CLI/Config.hpp>

#include <iomanip>
#include <iostream>
#include <limits>
#include <string>

/*
  NOTE, WARNING, ATTENTION

  The .../Random.hpp 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;

using Particle = setup::Stack::particle_type;

typedef decltype(1 * pascal) PressureType;
typedef decltype(1 * degree_celsius) TemperatureType;

class MarsAtmModel {
public:
  MarsAtmModel() = delete;
  MarsAtmModel(PressureType a, InverseLengthType b, TemperatureType c,
               decltype(1 * degree_celsius / 1_m) d)
      : a_(a)
      , b_(b)
      , c_(c)
      , d_(d) {}

  MassDensityType operator()(LengthType height) const {
    PressureType const pressure = a_ * exp(-b_ * height);
    TemperatureType const temperature = -c_ - d_ * height + 273.1_K; // in K
    constexpr decltype(square(1_m) / (square(1_s) * 1_K)) constant =
        1000 * 0.1921 * square(1_m) / (square(1_s) * 1_K);
    return pressure / (constant * temperature);
  }

private:
  PressureType a_;
  InverseLengthType b_;
  TemperatureType c_;
  decltype(1_K / 1_m) d_;
};

void registerRandomStreams(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) {
    std::random_device rd;
    seed = rd();
    cout << "new random seed (auto) " << seed << endl;
  }
  RNGManager<>::getInstance().setSeed(seed);
}

template <typename T>
using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;

// argv : 1.number of nucleons, 2.number of protons,
//        3.total energy in GeV, 4.number of showers,
//        5.seed (0 by default to generate random values for all)

int main(int argc, char** argv) {

  // the main command line description
  CLI::App app{"Simulate standard (downgoing) showers with CORSIKA 8."};

  // some options that we want to fill in
  int A, Z, nevent = 0;

  // the following section adds the options to the parser

  // we start by definining a sub-group for the primary ID
  auto opt_Z = app.add_option("-Z", Z, "Atomic number for primary")
                   ->check(CLI::Range(0, 26))
                   ->group("Primary");
  auto opt_A = app.add_option("-A", A, "Atomic mass number for primary")
                   ->needs(opt_Z)
                   ->check(CLI::Range(1, 58))
                   ->group("Primary");
  app.add_option("-p,--pdg", "PDG code for primary.")
      ->excludes(opt_A)
      ->excludes(opt_Z)
      ->group("Primary");
  // the remainding options
  app.add_option("-E,--energy", "Primary energy in GeV")
      ->required()
      ->check(CLI::PositiveNumber)
      ->group("Primary");
  app.add_option("-z,--zenith", "Primary zenith angle (deg)")
      ->required()
      ->default_val(0.)
      ->check(CLI::Range(0, 90))
      ->group("Primary");
  app.add_option("-a,--azimuth", "Primary azimuth angle (deg)")
      ->default_val(0.)
      ->check(CLI::Range(0, 360))
      ->group("Primary");
  app.add_option("-N,--nevent", nevent, "The number of events/showers to run.")
      ->required()
      ->check(CLI::PositiveNumber)
      ->group("Library/Output");
  app.add_option("-f,--filename", "Filename for output library.")
      ->required()
      ->default_val("corsika_library")
      ->check(CLI::NonexistentPath)
      ->group("Library/Output");
  app.add_option("-s,--seed", "The random number seed.")
      ->default_val(12351739)
      ->check(CLI::NonNegativeNumber)
      ->group("Misc.");
  app.add_flag("--force-interaction", "Force the location of the first interaction.")
      ->group("Misc.");
  app.add_option("-v,--verbosity", "Verbosity level: warn, info, debug, trace.")
      ->default_val("info")
      ->check(CLI::IsMember({"warn", "info", "debug", "trace"}))
      ->group("Misc.");

  // parse the command line options into the variables
  CLI11_PARSE(app, argc, argv);

  if (app.count("--verbosity")) {
    string const loglevel = app["verbosity"]->as<string>();
    if (loglevel == "warn") {
      logging::set_level(logging::level::warn);
    } else if (loglevel == "info") {
      logging::set_level(logging::level::info);
    } else if (loglevel == "debug") {
      logging::set_level(logging::level::debug);
    } else if (loglevel == "trace") {
#ifndef DEBUG
      CORSIKA_LOG_ERROR("trace log level requires a Debug build.");
      return 1;
#endif
      logging::set_level(logging::level::trace);
    }
  }

  // check that we got either PDG or A/Z
  // this can be done with option_groups but the ordering
  // gets all messed up
  if (app.count("--pdg") == 0) {
    if ((app.count("-A") == 0) || (app.count("-Z") == 0)) {
      std::cerr << "If --pdg is not provided, then both -A and -Z are required."
                << std::endl;
      return 1;
    }
  }

  // initialize random number sequence(s)
  registerRandomStreams(app["--seed"]->as<int>());

  /* === START: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
  using EnvType = setup::Environment;
  EnvType env;
  CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
  Point const center{rootCS, 0_m, 0_m, 0_m};
  LengthType const radiusMars = 3389.5_km;
  auto builder =
      make_layered_spherical_atmosphere_builder<setup::EnvironmentInterface, MyExtraEnv>::
          create(center,
                 radiusMars,                                   // Mars
                 Medium::AirDry1Atm,                           // Mars, close enough
                 MagneticFieldVector{rootCS, 0_T, 0_uT, 0_T}); // Mars

  builder.setNuclearComposition(                             // Mars
      {{Code::Nitrogen, Code::Oxygen}, {1. / 3., 2. / 3.}}); // simplified
  //{{Code::Carbon, Code::Oxygen, // 95.97 CO2
  //          Code::Nitrogen},            // 1.89 N2 + 1.93 Argon + 0.146 O2
  //       {0.9597 / 3, 0.9597 * 2 / 3,
  //      1 - 0.9597}}); // values taken from AIRES manual, Ar removed for now

  MarsAtmModel layer1(0.699e3 * pascal, 0.00009 / 1_m, 31.0 * degree_celsius,
                      0.000998 * 1 * degree_celsius / 1_m);
  MarsAtmModel layer2(0.699e3 * pascal, 0.00009 / 1_m, 23.4 * degree_celsius,
                      0.00222 * 1 * degree_celsius / 1_m);

  builder.addTabularLayer(layer1, 100, 100_m, 7_km);
  builder.addTabularLayer(layer2, 300, 500_m, 100_km);
  builder.addLinearLayer(1e9_cm, 112.8_km);
  builder.assemble(env);
  /* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */

  ofstream atmout("mars.dat");
  for (LengthType h = 0_m; h < 110_km; h += 100_m) {
    Point const ptest{rootCS, 0_m, 0_m, builder.getPlanetRadius() + h};
    auto rho =
        env.getUniverse()->getContainingNode(ptest)->getModelProperties().getMassDensity(
            ptest);
    atmout << h / 1_m << " " << rho / 1_kg * cube(1_m) << "\n";
  }
  atmout.close();

  /* === START: CONSTRUCT PRIMARY PARTICLE === */

  // parse the primary ID as a PDG or A/Z code
  Code beamCode;
  HEPEnergyType mass;

  // check if we want to use a PDG code instead
  if (app.count("--pdg") > 0) {
    beamCode = convert_from_PDG(PDGCode(app["--pdg"]->as<int>()));
    mass = get_mass(beamCode);
  } else {
    // check manually for proton and neutrons
    if ((A == 0) && (Z == 1)) beamCode = Code::Proton;
    if ((A == 1) && (Z == 1)) beamCode = Code::Neutron;
    mass = get_nucleus_mass(A, Z);
  }

  // particle energy
  HEPEnergyType const E0 = 1_GeV * app["--energy"]->as<float>();

  // direction of the shower in (theta, phi) space
  auto const thetaRad = app["--zenith"]->as<float>() / 180. * M_PI;
  auto const phiRad = app["--azimuth"]->as<float>() / 180. * M_PI;

  // convert Elab to Plab
  HEPMomentumType P0 = sqrt((E0 - mass) * (E0 + mass));

  // convert the momentum to the zenith and azimuth angle of the primary
  auto const [px, py, pz] =
      std::make_tuple(P0 * sin(thetaRad) * cos(phiRad), P0 * sin(thetaRad) * sin(phiRad),
                      -P0 * cos(thetaRad));
  auto plab = MomentumVector(rootCS, {px, py, pz});
  /* === END: CONSTRUCT PRIMARY PARTICLE === */

  /* === START: CONSTRUCT GEOMETRY === */
  auto const observationHeight = 0_km + builder.getPlanetRadius();
  auto const injectionHeight = 111.75_km + builder.getPlanetRadius();
  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) * cos(phiRad),
                                    -sin(thetaRad) * sin(phiRad), cos(thetaRad)}} *
                       t;

  // we make the axis much longer than the inj-core distance since the
  // profile will go beyond the core, depending on zenith angle
  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env};
  /* === END: CONSTRUCT GEOMETRY === */

  // create the output manager that we then register outputs with
  OutputManager output(app["--filename"]->as<std::string>());

  /* === START: SETUP PROCESS LIST === */
  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{1_GeV, 1_GeV, 1_GeV, 1_GeV, false};
  corsika::proposal::Interaction emCascade(env);
  corsika::proposal::ContinuousProcess emContinuous(env);
  InteractionCounter emCascadeCounted(emCascade);

  LongitudinalProfile longprof{showerAxis};

  corsika::urqmd::UrQMD urqmd;
  InteractionCounter urqmdCounted{urqmd};
  StackInspector<setup::Stack> stackInspect(5000, false, E0);

  // assemble all processes into an ordered process list
  struct EnergySwitch {
    HEPEnergyType cutE_;
    EnergySwitch(HEPEnergyType cutE)
        : cutE_(cutE) {}
    bool operator()(const Particle& p) { return (p.getKineticEnergy() < cutE_); }
  };
  auto hadronSequence = make_select(EnergySwitch(80_GeV), urqmdCounted,
                                    make_sequence(sibyllNucCounted, sibyllCounted));
  auto decaySequence = make_sequence(decayPythia, decaySibyll);

  // track writer
  TrackWriter trackWriter;
  output.add("tracks", trackWriter); // register TrackWriter

  // observation plane
  Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
  ObservationPlane<setup::Tracking, NoOutput> observationLevel(
      obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
  // register the observation plane with the output
  output.add("particles", observationLevel);

  // assemble the final process sequence
  auto sequence =
      make_sequence(stackInspect, hadronSequence, decaySequence, emCascadeCounted,
                    emContinuous, cut, trackWriter, observationLevel, longprof);
  /* === END: SETUP PROCESS LIST === */

  // create the cascade object using the default stack and tracking implementation
  setup::Tracking tracking;
  setup::Stack stack;
  Cascade EAS(env, tracking, sequence, output, stack);

  // print our primary parameters all in one place
  if (app["--pdg"]->count() > 0) {
    CORSIKA_LOG_INFO("Primary PDG ID: {}", app["--pdg"]->as<int>());
  } else {
    CORSIKA_LOG_INFO("Primary Z/A: {}/{}", Z, A);
  }
  CORSIKA_LOG_INFO("Primary Energy: {}", E0);
  CORSIKA_LOG_INFO("Primary Momentum: {}", P0);
  CORSIKA_LOG_INFO("Point of Injection: {}", injectionPos.getCoordinates());
  CORSIKA_LOG_INFO("Shower Axis Length: {}", (showerCore - injectionPos).getNorm() * 1.2);

  // trigger the output manager to open the library for writing
  output.startOfLibrary();

  // loop over each shower
  for (int i_shower = 1; i_shower < nevent + 1; i_shower++) {
    CORSIKA_LOG_INFO("Shower {} / {} ", i_shower, nevent);

    // trigger the start of the outputs for this shower
    output.startOfShower();

    // 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";

    // setup particle stack, and add primary particle
    stack.clear();

    // add the desired particle to the stack
    if (A > 1) {
      stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, A, Z));
    } else {
      stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
    }

    // if we want to fix the first location of the shower
    if (app["--force-interaction"]) EAS.forceInteraction();

    // run the shower
    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);

    // trigger the output manager to save this shower to disk
    output.endOfShower();
  }

  // and finalize the output on disk
  output.endOfLibrary();
}