IAP GITLAB

Skip to content
Snippets Groups Projects
cascade_proton_example.cpp 5.34 KiB
Newer Older
ralfulrich's avatar
ralfulrich committed
/*
 * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
ralfulrich's avatar
ralfulrich committed
 *
 * 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.
 */

ralfulrich's avatar
ralfulrich committed
#include <corsika/framework/process/ProcessSequence.hpp>
ralfulrich's avatar
ralfulrich committed
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/core/Cascade.hpp>
ralfulrich's avatar
ralfulrich committed
#include <corsika/framework/core/PhysicalUnits.hpp>
ralfulrich's avatar
ralfulrich committed
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
ralfulrich's avatar
ralfulrich committed
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
ralfulrich's avatar
ralfulrich committed

#include <corsika/output/OutputManager.hpp>
ralfulrich's avatar
ralfulrich committed
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
ralfulrich's avatar
ralfulrich committed

#include <corsika/media/Environment.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>

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

#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/StackInspector.hpp>
#include <corsika/modules/Sibyll.hpp>
#include <corsika/modules/ParticleCut.hpp>
#include <corsika/modules/TrackWriter.hpp>
#include <corsika/modules/HadronicElasticModel.hpp>
#include <corsika/modules/Pythia8.hpp>
ralfulrich's avatar
ralfulrich committed
/*
  NOTE, WARNING, ATTENTION
  The file Random.hpp implements the hooks of external modules to the C8 random
ralfulrich's avatar
ralfulrich committed
  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>
ralfulrich's avatar
ralfulrich committed

#include <iostream>
#include <limits>
#include <typeinfo>

using namespace corsika;
using namespace std;

//
// The example main program for a particle cascade
//
int main() {
ralfulrich's avatar
ralfulrich committed
  logging::set_level(logging::level::info);

  std::cout << "cascade_proton_example" << std::endl;

ralfulrich's avatar
ralfulrich committed
  feenableexcept(FE_INVALID);
  // initialize random number sequence(s)
  RNGManager<>::getInstance().registerRandomStream("cascade");
ralfulrich's avatar
ralfulrich committed

ralfulrich's avatar
ralfulrich committed
  OutputManager output("cascade_proton_outputs");

ralfulrich's avatar
ralfulrich committed
  // setup environment, geometry
  using EnvType = setup::Environment;
ralfulrich's avatar
ralfulrich committed
  EnvType env;
ralfulrich's avatar
ralfulrich committed
  auto& universe = *(env.getUniverse());
  CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
  auto world = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
ralfulrich's avatar
ralfulrich committed

ralfulrich's avatar
ralfulrich committed
  using MyHomogeneousModel = MediumPropertyModel<
      UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>;
ralfulrich's avatar
ralfulrich committed

  world->setModelProperties<MyHomogeneousModel>(
      Medium::AirDry1Atm, MagneticFieldVector(rootCS, 0_T, 0_T, 1_mT),
ralfulrich's avatar
ralfulrich committed
      1_kg / (1_m * 1_m * 1_m), NuclearComposition({Code::Hydrogen}, {1.}));
ralfulrich's avatar
ralfulrich committed

  universe.addChild(std::move(world));
ralfulrich's avatar
ralfulrich committed

  // setup particle stack, and add primary particle
  setup::Stack stack;
ralfulrich's avatar
ralfulrich committed
  stack.clear();
ralfulrich's avatar
ralfulrich committed
  const Code beamCode = Code::Proton;
ralfulrich's avatar
ralfulrich committed
  const HEPMassType mass = Proton::mass;
  const HEPEnergyType E0 = 200_GeV;
ralfulrich's avatar
ralfulrich committed
  double theta = 0.;
  double phi = 0.;

  Point injectionPos(rootCS, 0_m, 0_m, 0_m);
ralfulrich's avatar
ralfulrich committed
  {
    auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
      return sqrt(Elab * Elab - m * 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);
ralfulrich's avatar
ralfulrich committed
    auto plab = MomentumVector(rootCS, {px, py, pz});
ralfulrich's avatar
ralfulrich committed
    cout << "input particle: " << beamCode << endl;
    cout << "input angles: theta=" << theta << " phi=" << phi << endl;
ralfulrich's avatar
ralfulrich committed
    cout << "input momentum: " << plab.getComponents() / 1_GeV << endl;
ralfulrich's avatar
ralfulrich committed
    stack.addParticle(std::make_tuple(
        beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
        plab.normalized(), injectionPos, 0_ns));
ralfulrich's avatar
ralfulrich committed
  }

  // setup processes, decays and interactions
  setup::Tracking tracking;
ralfulrich's avatar
ralfulrich committed
  StackInspector<setup::Stack> stackInspect(1000, true, E0);
  RNGManager<>::getInstance().registerRandomStream("sibyll");
  RNGManager<>::getInstance().registerRandomStream("pythia");
ralfulrich's avatar
ralfulrich committed
  corsika::pythia8::Interaction pythia;
  corsika::pythia8::Decay decay;
ralfulrich's avatar
ralfulrich committed

  ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
ralfulrich's avatar
ralfulrich committed
  EnergyLossWriter dEdX{showerAxis};
ralfulrich's avatar
ralfulrich committed
  output.add("energyloss", dEdX);

  BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss{dEdX};
ralfulrich's avatar
ralfulrich committed
  ParticleCut<SubWriter<decltype(dEdX)>> cut(60_GeV, true, dEdX);
ralfulrich's avatar
ralfulrich committed
  cut.printThresholds();
ralfulrich's avatar
ralfulrich committed

  // RNGManager::getInstance().registerRandomStream("HadronicElasticModel");
  // HadronicElasticModel::HadronicElasticInteraction
ralfulrich's avatar
ralfulrich committed
  // hadronicElastic(env);

  TrackWriter trackWriter;
ralfulrich's avatar
ralfulrich committed
  output.add("tracks", trackWriter); // register TrackWriter

ralfulrich's avatar
ralfulrich committed
  // assemble all processes into an ordered process list
  auto sequence = make_sequence(pythia, decay, eLoss, cut, trackWriter, stackInspect);
ralfulrich's avatar
ralfulrich committed

  // define air shower object, run simulation
  Cascade EAS(env, tracking, sequence, output, stack);
  output.startOfShower();
ralfulrich's avatar
ralfulrich committed
  EAS.run();
  output.endOfShower();
ralfulrich's avatar
ralfulrich committed
  CORSIKA_LOG_INFO("Done");
ralfulrich's avatar
ralfulrich committed

  output.endOfLibrary();
ralfulrich's avatar
ralfulrich committed
}