cascade_experiment.cc

/*
 * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
 *
 * See file AUTHORS for a list of contributors.
 *
 * 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/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
#include <corsika/process/energy_loss/EnergyLoss.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/process/tracking_line/TrackingLine.h>

#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>

#include <corsika/environment/Environment.h>
#include <corsika/environment/BaseExponential.h>
#include <corsika/environment/FlatExponential.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>

#include <corsika/geometry/Sphere.h>

#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/sibyll/Interaction.h>

#include <corsika/process/particle_cut/ParticleCut.h>
#include <corsika/process/track_writer/TrackWriter.h>

#include <corsika/units/PhysicalUnits.h>

#include <corsika/random/RNGManager.h>

#include <corsika/utl/CorsikaFenv.h>

#include <iostream>
#include <limits>

using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;

using namespace std;
using namespace corsika::units::si;

//
// The example main program for a particle cascade
//
int main() {

  const LengthType height_atmosphere = 112.8_km;

  feenableexcept(FE_INVALID);
  // initialize random number sequence(s)
  random::RNGManager::GetInstance().RegisterRandomStream("cascade");

  // setup environment, geometry
  using EnvType = environment::Environment<corsika::environment::IMediumModel>;
  EnvType env;
  auto& universe = *(env.GetUniverse());

  const CoordinateSystem& rootCS = env.GetCoordinateSystem();

  // FlatExponential Modell
  auto theMedium = EnvType::CreateNode<Sphere>(
          Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());

  auto constexpr temperature = 295_K; // AIRES default temperature for isothermal model
  auto constexpr lambda =
          -constants::R * temperature / (constants::g_sub_n * 28.966_g / mole);
  // 1036 g/cm² taken from AIRES code
  auto constexpr rho0 = -1036_g / square(1_cm) / lambda;
  using FlatExp = environment::FlatExponential<environment::IMediumModel>;
  theMedium->SetModelProperties<FlatExp>(
    Point{rootCS, 0_m, 0_m, 0_m}, Vector<dimensionless_d>{rootCS, {0., 0., 1.}}, rho0,
          lambda,
          environment::NuclearComposition(
                  std::vector<particles::Code>{particles::Code::Nitrogen,
                                               particles::Code::Oxygen},
                  std::vector<float>{
                          0.7847f,
                          1.f - 0.7847f}
          )
  ); // values taken from AIRES manual, Ar removed for now

  universe.AddChild(std::move(theMedium));

  // setup particle stack
  setup::Stack stack;
  stack.Clear();

  // primary particle
  const Code beamCode = Code::Nucleus;
  const int nuclA = 4;
  const int nuclZ = int(nuclA / 2.15 + 0.7);
  const HEPMassType mass = GetNucleusMass(nuclA, nuclZ);

  HEPEnergyType E0;
  double theta;

  // particle angle
  for (unsigned int j = 0; j <= 70; j += 10) {
    theta = static_cast<double>(j);

    // particle energy
    for (HEPEnergyType i = 10_GeV; i <= 100_TeV; i += 100_GeV) {
        E0 = nuclA * i;

        // run number
        for (unsigned int run = 0; run < 1000; run++) {
          std::cout << "### work with theta=" << theta
                    << " E0=" << E0
                    << " run=" << run
                    << std::endl;

          // initialize random number sequence(s) new
          std::stringstream seed;
          seed << "cascade" << i << j << run;
          random::RNGManager::GetInstance().RegisterRandomStream(seed.str());

          double phi = 0.;

          {
            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(theta / 180. * M_PI, phi / 180. * M_PI, P0);
            auto plab = corsika::stack::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 << endl;
            Point pos(rootCS, 0_m, 0_m,
                      height_atmosphere); // this is the CORSIKA 7 start of atmosphere/universe
            stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
                    corsika::stack::MomentumVector, geometry::Point,
                    units::si::TimeType, unsigned short, unsigned short>{
                    beamCode, E0, plab, pos, 0_ns, nuclA, nuclZ});
          }

          // setup processes, decays and interactions
          tracking_line::TrackingLine tracking;
          stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);

          random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
          random::RNGManager::GetInstance().RegisterRandomStream("pythia");

          process::sibyll::Interaction sibyll;
          process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
          process::sibyll::Decay decay;

          // cascade with only HE model ==> HE cut
          process::particle_cut::ParticleCut cut(1_GeV);

          process::track_writer::TrackWriter trackWriter("tracks.dat");
          process::energy_loss::EnergyLoss eLoss;

          // assemble all processes into an ordered process list
          auto sequence = stackInspect << sibyll
                                       << sibyllNuc
                                       << decay
                                       << eLoss
                                       << cut
                                       << trackWriter;

          // define air shower object, run simulation
          cascade::Cascade EAS(env, tracking, sequence, stack);
          EAS.Init();
          EAS.Run();

          eLoss.PrintProfile(); // print longitudinal profile

          cut.ShowResults();
          const HEPEnergyType Efinal =
                  cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
          cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
               << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
          cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
               << "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
        }
      }
  }
}