corsika.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.
 */

#define TRACE

/* 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/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 <iomanip>
#include <iostream>
#include <limits>
#include <string>

/*
  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>
#include <corsika/modules/qgsjetII/Random.hpp>

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);
}

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) {

  logging::set_level(logging::level::info);

  CORSIKA_LOG_INFO("corsika");

  if (argc < 5) {
    std::cerr << "usage: corsika <A> <Z> <energy/GeV> <Nevt> [seed] \n"
                 "       if A=0, Z is interpreted as PDG code \n"
                 "       if no seed is given, a random seed is chosen \n"
              << std::endl;
    return 1;
  }
  feenableexcept(FE_INVALID);

  int seed = 0;
  int number_showers = std::stoi(std::string(argv[4]));

  if (argc > 5) { seed = std::stoi(std::string(argv[5])); }

  // 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,
                                                       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, 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);

  // pre-setup particle stack
  unsigned short const A = std::stoi(std::string(argv[1]));
  Code beamCode;
  HEPEnergyType mass;
  unsigned short Z = 0;
  if (A > 0) {
    beamCode = Code::Nucleus;
    Z = std::stoi(std::string(argv[2]));
    mass = get_nucleus_mass(A, Z);
  } else {
    int pdg = std::stoi(std::string(argv[2]));
    beamCode = convert_from_PDG(PDGCode(pdg));
    mass = get_mass(beamCode);
  }
  HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
  double theta = 20.;
  double phi = 180.;
  auto const thetaRad = theta / 180. * M_PI;
  auto const phiRad = phi / 180. * M_PI;

  auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
    return sqrt((Elab - m) * (Elab + 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(thetaRad, phiRad, 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
       << ", norm = " << plab.getNorm() << endl;

  auto const observationHeight = 0_km + builder.getEarthRadius();
  auto const injectionHeight = 111.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) * cos(phiRad),
                                    -sin(thetaRad) * sin(phiRad), cos(thetaRad)}} *
                       t;

  std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;

  // 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.2
            << std::endl;

  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env};

  // create the output manager that we then register outputs with
  OutputManager output("corsika_outputs");

  // setup processes, decays and interactions

  // corsika::qgsjetII::Interaction qgsjet;
  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{50_GeV, 50_GeV, 50_GeV, 50_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);

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

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

  auto sequence =
      make_sequence(stackInspect, hadronSequence, decaySequence, emCascadeCounted,
                    emContinuous, cut, trackWriter, observationLevel, longprof);

  // define air shower object, run simulation
  setup::Tracking tracking;
  setup::Stack stack;
  Cascade EAS(env, tracking, sequence, output, stack);

  output.startOfLibrary();

  for (int i_shower = 1; i_shower < number_showers + 1; i_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";

    CORSIKA_LOG_INFO("Shower {} / {} ", i_shower, number_showers);

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

    if (A > 1) {
      stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, A, Z));

    } else {
      if (A == 1) {
        if (Z == 1) {
          stack.addParticle(std::make_tuple(Code::Proton, plab, injectionPos, 0_ns));
        } else if (Z == 0) {
          stack.addParticle(std::make_tuple(Code::Neutron, plab, injectionPos, 0_ns));
        } else {
          std::cerr << "illegal parameters" << std::endl;
          return EXIT_FAILURE;
        }
      } else {
        stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
      }
    }

    // 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();

    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);
    output.endOfShower();
  }
  output.endOfLibrary();
}