Newer
Older
* (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.
*/
#include <corsika/framework/process/ProcessSequence.hpp>
#include <corsika/framework/process/SwitchProcessSequence.hpp>
#include <corsika/framework/core/Cascade.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/EnergyMomentumOperations.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/geometry/Sphere.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
#include <corsika/modules/writers/SubWriter.hpp>
#include <corsika/modules/writers/EnergyLossWriter.hpp>
#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>
The file Random.hpp implements 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>
#include <iostream>
#include <limits>
#include <typeinfo>
using namespace corsika;
using namespace std;
//
// The example main program for a particle cascade
//
int main() {
std::cout << "cascade_proton_example" << std::endl;
feenableexcept(FE_INVALID);
// initialize random number sequence(s)
RNGManager<>::getInstance().registerRandomStream("cascade");
auto& universe = *(env.getUniverse());
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
auto world = EnvType::createNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 150_km);
using MyHomogeneousModel = MediumPropertyModel<
UniformMagneticField<HomogeneousMedium<setup::EnvironmentInterface>>>;
Medium::AirDry1Atm, MagneticFieldVector(rootCS, 0_T, 0_T, 1_mT),
1_kg / (1_m * 1_m * 1_m), NuclearComposition({Code::Hydrogen}, {1.}));
// setup particle stack, and add primary particle
setup::Stack stack;
{
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);
cout << "input particle: " << beamCode << endl;
cout << "input angles: theta=" << theta << " phi=" << phi << endl;
cout << "input momentum: " << plab.getComponents() / 1_GeV << endl;
stack.addParticle(std::make_tuple(
beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)),
plab.normalized(), injectionPos, 0_ns));
}
// setup processes, decays and interactions
StackInspector<setup::Stack> stackInspect(1000, true, E0);
RNGManager<>::getInstance().registerRandomStream("sibyll");
RNGManager<>::getInstance().registerRandomStream("pythia");
corsika::pythia8::Interaction pythia;
corsika::pythia8::Decay decay;
ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env};
output.add("energyloss", dEdX);
BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss{dEdX};
ParticleCut<SubWriter<decltype(dEdX)>> cut(60_GeV, true, dEdX);
// RNGManager::getInstance().registerRandomStream("HadronicElasticModel");
// HadronicElasticModel::HadronicElasticInteraction
output.add("tracks", trackWriter); // register TrackWriter
// assemble all processes into an ordered process list
auto sequence = make_sequence(pythia, decay, eLoss, cut, trackWriter, stackInspect);
Cascade EAS(env, tracking, sequence, output, stack);