IAP GITLAB

Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • AirShowerPhysics/corsika
  • rulrich/corsika
  • AAAlvesJr/corsika
  • Andre/corsika
  • arrabito/corsika
  • Nikos/corsika
  • olheiser73/corsika
  • AirShowerPhysics/papers/corsika
  • pranav/corsika
9 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 178 additions and 2603 deletions
Processes/Sibyll/Interaction.h deleted 100644 → 0
View file @ 036d63e8
/*
* (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.
*/
#ifndef _corsika_process_sibyll_interaction_h_
#define _corsika_process_sibyll_interaction_h_
#include <corsika/process/InteractionProcess.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/FourVector.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
namespace corsika::process::sibyll {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
int fCount = 0;
int fNucCount = 0;
bool fInitialized = false;
public:
Interaction(corsika::environment::Environment const& env)
: fEnvironment(env) {}
~Interaction() {
std::cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount
<< std::endl;
}
void Init() {
using corsika::random::RNGManager;
using std::cout;
using std::endl;
// initialize Sibyll
if (!fInitialized) {
sibyll_ini_();
fInitialized = true;
}
}
bool WasInitialized() { return fInitialized; }
bool ValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
using namespace corsika::units::si;
return (10_GeV < ecm) && (ecm < 1_PeV);
}
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType,
int>
GetCrossSection(const corsika::particles::Code BeamId,
const corsika::particles::Code TargetId,
const corsika::units::si::HEPEnergyType CoMenergy) {
using namespace corsika::units::si;
double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3];
const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
const double dEcm = CoMenergy / 1_GeV;
int iTarget = -1;
if (corsika::particles::IsNucleus(TargetId)) {
iTarget = corsika::particles::GetNucleusA(TargetId);
if (iTarget > 18 || iTarget == 0)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 are allowed.");
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
} else if (TargetId == corsika::particles::Proton::GetCode()) {
sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
iTarget = 1;
} else {
// no interaction in sibyll possible, return infinite cross section? or throw?
sigProd = std::numeric_limits<double>::infinity();
sigEla = std::numeric_limits<double>::infinity();
}
return std::make_tuple(sigProd * 1_mbarn, sigEla * 1_mbarn, iTarget);
}
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&) {
using namespace corsika::units;
using namespace corsika::units::si;
using namespace corsika::geometry;
using std::cout;
using std::endl;
// coordinate system, get global frame of reference
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
const particles::Code corsikaBeamId = p.GetPID();
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
corsika::particles::Neutron::GetMass());
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy
HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
pTotLab += p.GetMomentum();
pTotLab += pTarget;
auto const pTotLabNorm = pTotLab.norm();
// calculate cm. energy
const HEPEnergyType ECoM = sqrt(
(Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
std::cout << "Interaction: LambdaInt: \n"
<< " input energy: " << p.GetEnergy() / 1_GeV << endl
<< " beam can interact:" << kInteraction << endl
<< " beam pid:" << p.GetPID() << endl;
// TODO: move limits into variables
if (kInteraction && Elab >= 8.5_GeV && ECoM >= 10_GeV) {
// get target from environment
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
const auto currentNode =
fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
const auto mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// determine average interaction length
// weighted sum
int i = -1;
double avgTargetMassNumber = 0.;
si::CrossSectionType weightedProdCrossSection = 0_mbarn;
// get weights of components from environment/medium
const auto w = mediumComposition.GetFractions();
// loop over components in medium
for (auto const targetId : mediumComposition.GetComponents()) {
i++;
cout << "Interaction: get interaction length for target: " << targetId << endl;
auto const [productionCrossSection, elaCrossSection, numberOfNucleons] =
GetCrossSection(corsikaBeamId, targetId, ECoM);
[[maybe_unused]] auto elaCrossSectionCopy =
elaCrossSection; // ONLY TO AVOID COMPILER WARNING
std::cout << "Interaction: "
<< " IntLength: sibyll return (mb): "
<< productionCrossSection / 1_mbarn << std::endl;
weightedProdCrossSection += w[i] * productionCrossSection;
avgTargetMassNumber += w[i] * numberOfNucleons;
}
cout << "Interaction: "
<< "IntLength: weighted CrossSection (mb): "
<< weightedProdCrossSection / 1_mbarn << endl;
// calculate interaction length in medium
//#warning check interaction length units
GrammageType const int_length =
avgTargetMassNumber * corsika::units::constants::u / weightedProdCrossSection;
std::cout << "Interaction: "
<< "interaction length (g/cm2): "
<< int_length / (0.001_kg) * 1_cm * 1_cm << std::endl;
return int_length;
}
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
/**
In this function SIBYLL is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
*/
template <typename Particle, typename Stack>
corsika::process::EProcessReturn DoInteraction(Particle& p, Stack&) {
using namespace corsika::units;
using namespace corsika::utl;
using namespace corsika::units::si;
using namespace corsika::geometry;
using std::cout;
using std::endl;
const auto corsikaBeamId = p.GetPID();
cout << "ProcessSibyll: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
<< process::sibyll::CanInteract(corsikaBeamId) << endl;
if (corsika::particles::IsNucleus(corsikaBeamId)) {
// nuclei handled by different process, this should not happen
throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!");
}
if (process::sibyll::CanInteract(corsikaBeamId)) {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in Sibyll
Point pOrig = p.GetPosition();
TimeType tOrig = p.GetTime();
// define target
// for Sibyll is always a single nucleon
auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
corsika::particles::Neutron::GetMass());
// FOR NOW: target is always at rest
const auto eTargetLab = 0_GeV + nucleon_mass;
const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargLab(eTargetLab, pTargetLab);
// define projectile
HEPEnergyType const eProjectileLab = p.GetEnergy();
auto const pProjectileLab = p.GetMomentum();
cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
<< "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
<< endl;
cout << "Interaction: etarget lab: " << eTargetLab / 1_GeV << endl
<< "Interaction: ptarget lab: " << pTargetLab.GetComponents() / 1_GeV
<< endl;
const FourVector PprojLab(eProjectileLab, pProjectileLab);
// define target kinematics in lab frame
// define boost to and from CoM frame
// CoM frame definition in Sibyll projectile: +z
COMBoost const boost(PprojLab, nucleon_mass);
// just for show:
// boost projecticle
auto const PprojCoM = boost.toCoM(PprojLab);
// boost target
auto const PtargCoM = boost.toCoM(PtargLab);
cout << "Interaction: ebeam CoM: " << PprojCoM.GetTimeLikeComponent() / 1_GeV
<< endl
<< "Interaction: pbeam CoM: "
<< PprojCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
cout << "Interaction: etarget CoM: " << PtargCoM.GetTimeLikeComponent() / 1_GeV
<< endl
<< "Interaction: ptarget CoM: "
<< PtargCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
cout << "Interaction: position of interaction: " << pOrig.GetCoordinates()
<< endl;
cout << "Interaction: time: " << tOrig << endl;
HEPEnergyType Etot = eProjectileLab + eTargetLab;
MomentumVector Ptot = p.GetMomentum();
// invariant mass, i.e. cm. energy
HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
// sample target mass number
const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// get cross sections for target materials
/*
Here we read the cross section from the interaction model again,
should be passed from GetInteractionLength if possible
*/
//#warning reading interaction cross section again, should not be necessary
auto const& compVec = mediumComposition.GetComponents();
std::vector<si::CrossSectionType> cross_section_of_components(compVec.size());
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
const auto [sigProd, sigEla, nNuc] =
GetCrossSection(corsikaBeamId, targetId, Ecm);
cross_section_of_components[i] = sigProd;
[[maybe_unused]] int ideleteme =
nNuc; // to avoid not used warning in array binding
[[maybe_unused]] auto sigElaCopy =
sigEla; // to avoid not used warning in array binding
}
const auto targetCode = currentNode->GetModelProperties().SampleTarget(
cross_section_of_components, fRNG);
cout << "Interaction: target selected: " << targetCode << endl;
/*
FOR NOW: allow nuclei with A<18 or protons only.
when medium composition becomes more complex, approximations will have to be
allowed air in atmosphere also contains some Argon.
*/
int targetSibCode = -1;
if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode);
if (targetCode == corsika::particles::Proton::GetCode()) targetSibCode = 1;
cout << "Interaction: sibyll code: " << targetSibCode << endl;
if (targetSibCode > 18 || targetSibCode < 1)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 or protons are "
"allowed.");
// beam id for sibyll
const int kBeam = process::sibyll::ConvertToSibyllRaw(corsikaBeamId);
std::cout << "Interaction: "
<< " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
<< " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
if (eProjectileLab < 8.5_GeV || Ecm < 10_GeV) {
std::cout << "Interaction: "
<< " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << std::endl;
throw std::runtime_error("energy too low for SIBYLL");
// p.Delete(); delete later... different process
} else {
fCount++;
// Sibyll does not know about units..
const double sqs = Ecm / 1_GeV;
// running sibyll, filling stack
sibyll_(kBeam, targetSibCode, sqs);
// running decays
// setTrackedParticlesStable();
decsib_();
// print final state
int print_unit = 6;
sib_list_(print_unit);
fNucCount += get_nwounded() - 1;
// delete current particle
p.Delete();
// add particles from sibyll to stack
// link to sibyll stack
SibStack ss;
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
// skip particles that have decayed in Sibyll
if (psib.HasDecayed()) continue;
// transform energy to lab. frame
auto const pCoM = psib.GetMomentum();
HEPEnergyType const eCoM = psib.GetEnergy();
auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
// add to corsika stack
auto pnew = p.AddSecondary(process::sibyll::ConvertFromSibyll(psib.GetPID()),
Plab.GetTimeLikeComponent(),
Plab.GetSpaceLikeComponents(), pOrig, tOrig);
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
Ecm_final += psib.GetEnergy();
}
std::cout << "conservation (all GeV): Ecm_final=" << Ecm_final / 1_GeV
<< std::endl
<< "Elab_final=" << Elab_final / 1_GeV
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents()
<< std::endl;
}
}
return process::EProcessReturn::eOk;
}
private:
corsika::environment::Environment const& fEnvironment;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
};
} // namespace corsika::process::sibyll
#endif
Processes/Sibyll/NuclearInteraction.h deleted 100644 → 0
View file @ 036d63e8
/*
* (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.
*/
#ifndef _corsika_process_sibyll_nuclearinteraction_h_
#define _corsika_process_sibyll_nuclearinteraction_h_
#include <corsika/process/InteractionProcess.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/sibyll/nuclib.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
namespace corsika::process::sibyll {
class NuclearInteraction
: public corsika::process::InteractionProcess<NuclearInteraction> {
int fCount = 0;
int fNucCount = 0;
public:
NuclearInteraction(corsika::environment::Environment const& env,
corsika::process::sibyll::Interaction& hadint)
: fEnvironment(env)
, fHadronicInteraction(hadint) {}
~NuclearInteraction() {
std::cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount
<< std::endl;
}
void Init() {
using corsika::random::RNGManager;
using std::cout;
using std::endl;
// initialize hadronic interaction module
// TODO: safe to run multiple initializations?
if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init();
// initialize nuclib
// TODO: make sure this does not overlap with sibyll
nuc_nuc_ini_();
}
// TODO: remove number of nucleons, avg target mass is available in environment
template <typename Particle>
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(Particle const& p, const corsika::particles::Code TargetId) {
using namespace corsika::units::si;
double sigProd;
auto const pCode = p.GetPID();
if (pCode != corsika::particles::Code::Nucleus)
throw std::runtime_error(
"NuclearInteraction: GetCrossSection: particle not a nucleus!");
auto const iBeam = p.GetNuclearA();
HEPEnergyType LabEnergyPerNuc = p.GetEnergy() / iBeam;
std::cout << "NuclearInteraction: GetCrossSection: called with: beamNuclA= "
<< iBeam << " TargetId= " << TargetId
<< " LabEnergyPerNuc= " << LabEnergyPerNuc / 1_GeV << std::endl;
// use nuclib to calc. nuclear cross sections
// TODO: for now assumes air with hard coded composition
// extend to arbitrary mixtures, requires smarter initialization
// get nuclib projectile code: nucleon number
if (iBeam > 56 || iBeam < 2) {
std::cout << "NuclearInteraction: beam nucleus outside allowed range for NUCLIB!"
<< std::endl
<< "A=" << iBeam << std::endl;
throw std::runtime_error(
"NuclearInteraction: GetCrossSection: beam nucleus outside allowed range for "
"NUCLIB!");
}
const double dElabNuc = LabEnergyPerNuc / 1_GeV;
// TODO: these limitations are still sibyll specific.
// available target nuclei depends on the hadronic interaction model and the
// initialization
if (dElabNuc < 10.)
throw std::runtime_error("NuclearInteraction: GetCrossSection: energy too low!");
// TODO: these limitations are still sibyll specific.
// available target nuclei depends on the hadronic interaction model and the
// initialization
if (corsika::particles::IsNucleus(TargetId)) {
const int iTarget = corsika::particles::GetNucleusA(TargetId);
if (iTarget > 18 || iTarget == 0)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 are allowed.");
std::cout << "NuclearInteraction: calling signuc.." << std::endl;
std::cout << "WARNING: using hard coded cross section for Nucleus-Air with "
"SIBYLL! (fix me!)"
<< std::endl;
// TODO: target id is not used because cross section is still hard coded and fixed
// to air.
signuc_(iBeam, dElabNuc, sigProd);
std::cout << "cross section: " << sigProd << std::endl;
return std::make_tuple(sigProd * 1_mbarn, 0_mbarn);
}
return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mbarn,
std::numeric_limits<double>::infinity() * 1_mbarn);
}
template <typename Particle, typename Track>
corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&) {
using namespace corsika::units;
using namespace corsika::units::si;
using namespace corsika::geometry;
using std::cout;
using std::endl;
// coordinate system, get global frame of reference
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
const particles::Code corsikaBeamId = p.GetPID();
if (!corsika::particles::IsNucleus(corsikaBeamId)) {
// no nuclear interaction
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
// check if target-style nucleus (enum)
if (corsikaBeamId != corsika::particles::Code::Nucleus)
throw std::runtime_error(
"NuclearInteraction: GetInteractionLength: Wrong nucleus type. Nuclear "
"projectiles should use NuclearStackExtension!");
// read from cross section code table
const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
corsika::particles::Neutron::GetMass());
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
// total momentum and energy
HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
int const nuclA = p.GetNuclearA();
auto const ElabNuc = p.GetEnergy() / nuclA;
MomentumVector pTotLab(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
pTotLab += p.GetMomentum();
pTotLab += pTarget;
auto const pTotLabNorm = pTotLab.norm();
// calculate cm. energy
const HEPEnergyType ECoM = sqrt(
(Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
auto const ECoMNN = sqrt(2. * ElabNuc * nucleon_mass);
std::cout << "NuclearInteraction: LambdaInt: \n"
<< " input energy: " << Elab / 1_GeV << endl
<< " input energy CoM: " << ECoM / 1_GeV << endl
<< " beam pid:" << corsikaBeamId << endl
<< " beam A: " << nuclA << endl
<< " input energy per nucleon: " << ElabNuc / 1_GeV << endl
<< " input energy CoM per nucleon: " << ECoMNN / 1_GeV << endl;
// throw std::runtime_error("stop here");
// energy limits
// TODO: values depend on hadronic interaction model !! this is sibyll specific
if (ElabNuc >= 8.5_GeV && ECoMNN >= 10_GeV) {
// get target from environment
/*
the target should be defined by the Environment,
ideally as full particle object so that the four momenta
and the boosts can be defined..
*/
const auto currentNode =
fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
const auto mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// determine average interaction length
// weighted sum
int i = -1;
si::CrossSectionType weightedProdCrossSection = 0_mbarn;
// get weights of components from environment/medium
const auto w = mediumComposition.GetFractions();
// loop over components in medium
for (auto const targetId : mediumComposition.GetComponents()) {
i++;
cout << "NuclearInteraction: get interaction length for target: " << targetId
<< endl;
auto const [productionCrossSection, elaCrossSection] =
GetCrossSection(p, targetId);
[[maybe_unused]] auto elaCrossSectionCopy = elaCrossSection;
std::cout << "NuclearInteraction: "
<< "IntLength: nuclib return (mb): "
<< productionCrossSection / 1_mbarn << std::endl;
weightedProdCrossSection += w[i] * productionCrossSection;
}
cout << "NuclearInteraction: "
<< "IntLength: weighted CrossSection (mb): "
<< weightedProdCrossSection / 1_mbarn << endl;
// calculate interaction length in medium
GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
corsika::units::constants::u /
weightedProdCrossSection;
std::cout << "NuclearInteraction: "
<< "interaction length (g/cm2): "
<< int_length / (0.001_kg) * 1_cm * 1_cm << std::endl;
return int_length;
} else {
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
}
template <typename Particle, typename Stack>
corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) {
// this routine superimposes different nucleon-nucleon interactions
// in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
using namespace corsika::units;
using namespace corsika::utl;
using namespace corsika::units::si;
using namespace corsika::geometry;
using std::cout;
using std::endl;
const auto ProjId = p.GetPID();
// TODO: calculate projectile mass in nuclearStackExtension
// const auto ProjMass = p.GetMass();
cout << "NuclearInteraction: DoInteraction: called with:" << ProjId << endl;
if (!IsNucleus(ProjId)) {
cout << "WARNING: non nuclear projectile in NUCLIB!" << endl;
// this should not happen
// throw std::runtime_error("Non nuclear projectile in NUCLIB!");
return process::EProcessReturn::eOk;
}
// check if target-style nucleus (enum)
if (ProjId != corsika::particles::Code::Nucleus)
throw std::runtime_error(
"NuclearInteraction: DoInteraction: Wrong nucleus type. Nuclear projectiles "
"should use NuclearStackExtension!");
auto const ProjMass =
p.GetNuclearZ() * corsika::particles::Proton::GetMass() +
(p.GetNuclearA() - p.GetNuclearZ()) * corsika::particles::Neutron::GetMass();
cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << endl;
fCount++;
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in NUCLIB
Point pOrig = p.GetPosition();
TimeType tOrig = p.GetTime();
cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "Interaction: time: " << tOrig << endl;
// projectile nucleon number
const int kAProj = p.GetNuclearA(); // GetNucleusA(ProjId);
if (kAProj > 56)
throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
// kinematics
// define projectile nucleus
HEPEnergyType const eProjectileLab = p.GetEnergy();
auto const pProjectileLab = p.GetMomentum();
const FourVector PprojLab(eProjectileLab, pProjectileLab);
cout << "NuclearInteraction: eProj lab: " << eProjectileLab / 1_GeV << endl
<< "NuclearInteraction: pProj lab: " << pProjectileLab.GetComponents() / 1_GeV
<< endl;
// define projectile nucleon
HEPEnergyType const eProjectileNucLab = p.GetEnergy() / kAProj;
auto const pProjectileNucLab = p.GetMomentum() / kAProj;
const FourVector PprojNucLab(eProjectileNucLab, pProjectileNucLab);
cout << "NuclearInteraction: eProjNucleon lab: " << eProjectileNucLab / 1_GeV
<< endl
<< "NuclearInteraction: pProjNucleon lab: "
<< pProjectileNucLab.GetComponents() / 1_GeV << endl;
// define target
// always a nucleon
auto constexpr nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
corsika::particles::Neutron::GetMass());
// target is always at rest
const auto eTargetNucLab = 0_GeV + nucleon_mass;
const auto pTargetNucLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargNucLab(eTargetNucLab, pTargetNucLab);
cout << "NuclearInteraction: etarget lab: " << eTargetNucLab / 1_GeV << endl
<< "NuclearInteraction: ptarget lab: " << pTargetNucLab.GetComponents() / 1_GeV
<< endl;
// center-of-mass energy in nucleon-nucleon frame
auto const PtotNN4 = PtargNucLab + PprojNucLab;
HEPEnergyType EcmNN = PtotNN4.GetNorm();
cout << "NuclearInteraction: nuc-nuc cm energy: " << EcmNN / 1_GeV << endl;
if (!fHadronicInteraction.ValidCoMEnergy(EcmNN)) {
cout << "NuclearInteraction: nuc-nuc. CoM energy too low for hadronic "
"interaction model!"
<< endl;
throw std::runtime_error("NuclearInteraction: DoInteraction: energy too low!");
}
// define boost to NUCLEON-NUCLEON frame
COMBoost const boost(PprojNucLab, nucleon_mass);
// boost projecticle
auto const PprojNucCoM = boost.toCoM(PprojNucLab);
// boost target
auto const PtargNucCoM = boost.toCoM(PtargNucLab);
cout << "Interaction: ebeam CoM: " << PprojNucCoM.GetTimeLikeComponent() / 1_GeV
<< endl
<< "Interaction: pbeam CoM: "
<< PprojNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
cout << "Interaction: etarget CoM: " << PtargNucCoM.GetTimeLikeComponent() / 1_GeV
<< endl
<< "Interaction: ptarget CoM: "
<< PtargNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
// sample target nucleon number
//
// proton stand-in for nucleon
const auto beamId = corsika::particles::Proton::GetCode();
const auto currentNode = fEnvironment.GetUniverse()->GetContainingNode(pOrig);
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
cout << "get nucleon-nucleus cross sections for target materials.." << endl;
// get cross sections for target materials
// using nucleon-target-nucleus cross section!!!
/*
Here we read the cross section from the interaction model again,
should be passed from GetInteractionLength if possible
*/
auto const& compVec = mediumComposition.GetComponents();
std::vector<si::CrossSectionType> cross_section_of_components(compVec.size());
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
cout << "target component: " << targetId << endl;
cout << "beam id: " << beamId << endl;
const auto [sigProd, sigEla, nNuc] =
fHadronicInteraction.GetCrossSection(beamId, targetId, EcmNN);
cross_section_of_components[i] = sigProd;
[[maybe_unused]] auto sigElaCopy = sigEla; // ONLY TO AVOID COMPILER WARNINGS
[[maybe_unused]] auto sigNucCopy = nNuc; // ONLY TO AVOID COMPILER WARNINGS
}
const auto targetCode = currentNode->GetModelProperties().SampleTarget(
cross_section_of_components, fRNG);
cout << "Interaction: target selected: " << targetCode << endl;
/*
FOR NOW: allow nuclei with A<18 or protons only.
when medium composition becomes more complex, approximations will have to be
allowed air in atmosphere also contains some Argon.
*/
int kATarget = -1;
if (IsNucleus(targetCode)) kATarget = GetNucleusA(targetCode);
if (targetCode == corsika::particles::Proton::GetCode()) kATarget = 1;
cout << "NuclearInteraction: nuclib target code: " << kATarget << endl;
if (kATarget > 18 || kATarget < 1)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 or protons are "
"allowed.");
// end of target sampling
// superposition
cout << "NuclearInteraction: sampling nuc. multiple interaction structure.. "
<< endl;
// get nucleon-nucleon cross section
// (needed to determine number of nucleon-nucleon scatterings)
const auto protonId = corsika::particles::Proton::GetCode();
const auto [prodCrossSection, elaCrossSection, dum] =
fHadronicInteraction.GetCrossSection(protonId, protonId, EcmNN);
[[maybe_unused]] auto dumCopy = dum; // ONLY TO AVOID COMPILER WARNING
const double sigProd = prodCrossSection / 1_mbarn;
const double sigEla = elaCrossSection / 1_mbarn;
// sample number of interactions (only input variables, output in common cnucms)
// nuclear multiple scattering according to glauber (r.i.p.)
int_nuc_(kATarget, kAProj, sigProd, sigEla);
cout << "number of nucleons in target : " << kATarget << endl
<< "number of wounded nucleons in target : " << cnucms_.na << endl
<< "number of nucleons in projectile : " << kAProj << endl
<< "number of wounded nucleons in project. : " << cnucms_.nb << endl
<< "number of inel. nuc.-nuc. interactions : " << cnucms_.ni << endl
<< "number of elastic nucleons in target : " << cnucms_.nael << endl
<< "number of elastic nucleons in project. : " << cnucms_.nbel << endl
<< "impact parameter: " << cnucms_.b << endl;
// calculate fragmentation
cout << "calculating nuclear fragments.." << endl;
// number of interactions
// include elastic
const int nElasticNucleons = cnucms_.nbel;
const int nInelNucleons = cnucms_.nb;
const int nIntProj = nInelNucleons + nElasticNucleons;
const double impactPar = cnucms_.b; // only needed to avoid passing common var.
int nFragments;
// number of fragments is limited to 60
int AFragments[60];
// call fragmentation routine
// input: target A, projectile A, number of int. nucleons in projectile, impact
// parameter (fm) output: nFragments, AFragments in addition the momenta ar stored
// in pf in common fragments, neglected
fragm_(kATarget, kAProj, nIntProj, impactPar, nFragments, AFragments);
// this should not occur but well :)
if (nFragments > 60)
throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
cout << "number of fragments: " << nFragments << endl;
for (int j = 0; j < nFragments; ++j)
cout << "fragment: " << j << " A=" << AFragments[j]
<< " px=" << fragments_.ppp[j][0] << " py=" << fragments_.ppp[j][1]
<< " pz=" << fragments_.ppp[j][2] << endl;
cout << "adding nuclear fragments to particle stack.." << endl;
// put nuclear fragments on corsika stack
for (int j = 0; j < nFragments; ++j) {
corsika::particles::Code specCode;
const auto nuclA = AFragments[j];
// get Z from stability line
const auto nuclZ = int(nuclA / 2.15 + 0.7);
// TODO: do we need to catch single nucleons??
if (nuclA == 1)
// TODO: sample neutron or proton
specCode = corsika::particles::Code::Proton;
else
specCode = corsika::particles::Code::Nucleus;
// TODO: mass of nuclei?
const HEPMassType mass =
corsika::particles::Proton::GetMass() * nuclZ +
(nuclA - nuclZ) *
corsika::particles::Neutron::GetMass(); // this neglects binding energy
cout << "NuclearInteraction: adding fragment: " << specCode << endl;
cout << "NuclearInteraction: A,Z: " << nuclA << "," << nuclZ << endl;
cout << "NuclearInteraction: mass: " << mass / 1_GeV << endl;
// CORSIKA 7 way
// spectators inherit momentum from original projectile
const double mass_ratio = mass / ProjMass;
cout << "NuclearInteraction: mass ratio " << mass_ratio << endl;
auto const Plab = PprojLab * mass_ratio;
cout << "NuclearInteraction: fragment momentum: "
<< Plab.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
if (nuclA == 1)
// add nucleon
p.AddSecondary(specCode, Plab.GetTimeLikeComponent(),
Plab.GetSpaceLikeComponents(), pOrig, tOrig);
else
// add nucleus
p.AddSecondary(specCode, Plab.GetTimeLikeComponent(),
Plab.GetSpaceLikeComponents(), pOrig, tOrig, nuclA, nuclZ);
}
// add elastic nucleons to corsika stack
// TODO: the elastic interaction could be external like the inelastic interaction,
// e.g. use existing ElasticModel
cout << "adding elastically scattered nucleons to particle stack.." << endl;
for (int j = 0; j < nElasticNucleons; ++j) {
// TODO: sample proton or neutron
auto const elaNucCode = corsika::particles::Code::Proton;
// CORSIKA 7 way
// elastic nucleons inherit momentum from original projectile
// neglecting momentum transfer in interaction
const double mass_ratio = corsika::particles::GetMass(elaNucCode) / ProjMass;
auto const Plab = PprojLab * mass_ratio;
p.AddSecondary(elaNucCode, Plab.GetTimeLikeComponent(),
Plab.GetSpaceLikeComponents(), pOrig, tOrig);
}
// add inelastic interactions
cout << "calculate inelastic nucleon-nucleon interactions.." << endl;
for (int j = 0; j < nInelNucleons; ++j) {
// TODO: sample neutron or proton
auto pCode = corsika::particles::Proton::GetCode();
// temporarily add to stack, will be removed after interaction in DoInteraction
cout << "inelastic interaction no. " << j << endl;
auto inelasticNucleon =
p.AddSecondary(pCode, PprojNucLab.GetTimeLikeComponent(),
PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig);
// create inelastic interaction
cout << "calling HadronicInteraction..." << endl;
fHadronicInteraction.DoInteraction(inelasticNucleon, s);
}
// delete parent particle
p.Delete();
cout << "NuclearInteraction: DoInteraction: done" << endl;
return process::EProcessReturn::eOk;
}
private:
corsika::environment::Environment const& fEnvironment;
corsika::process::sibyll::Interaction& fHadronicInteraction;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
};
} // namespace corsika::process::sibyll
#endif
Processes/Sibyll/ParticleConversion.cc deleted 100644 → 0
View file @ 036d63e8
/*
* (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/particles/ParticleProperties.h>
#include <corsika/process/sibyll/ParticleConversion.h>
using namespace corsika::process::sibyll;
// const std::map<sibyll::PID, ParticleProperties::InternalParticleCode>
// process::sibyll::Sibyll2Corsika = {
// {PID::E_MINUS, InternalParticleCode::Electron},
//};
Processes/Sibyll/ParticleConversion.h deleted 100644 → 0
View file @ 036d63e8
/*
* (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.
*/
#ifndef _include_processes_sibyll_particles_h_
#define _include_processes_sibyll_particles_h_
#include <corsika/particles/ParticleProperties.h>
#include <bitset2/bitset2.hpp>
#include <map>
namespace corsika::process::sibyll {
enum class SibyllCode : int8_t;
using SibyllCodeIntType = std::underlying_type<SibyllCode>::type;
#include <corsika/process/sibyll/Generated.inc>
bool constexpr KnownBySibyll(corsika::particles::Code pCode) {
return isKnown[static_cast<corsika::particles::CodeIntType>(pCode)];
}
bool constexpr CanInteract(corsika::particles::Code pCode) {
return canInteract[static_cast<corsika::particles::CodeIntType>(pCode)];
}
SibyllCode constexpr ConvertToSibyll(corsika::particles::Code pCode) {
return static_cast<SibyllCode>(
corsika2sibyll[static_cast<corsika::particles::CodeIntType>(pCode)]);
}
corsika::particles::Code constexpr ConvertFromSibyll(SibyllCode pCode) {
return sibyll2corsika[static_cast<SibyllCodeIntType>(pCode) - minSibyll];
}
int constexpr ConvertToSibyllRaw(corsika::particles::Code pCode) {
return static_cast<int>(ConvertToSibyll(pCode));
}
int constexpr GetSibyllXSCode(corsika::particles::Code pCode) {
return corsika2sibyllXStype[static_cast<corsika::particles::CodeIntType>(pCode)];
}
} // namespace corsika::process::sibyll
#endif
Processes/Sibyll/SibStack.h deleted 100644 → 0
View file @ 036d63e8
/*
* (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.
*/
#ifndef _include_sibstack_h_
#define _include_sibstack_h_
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/sibyll2.3c.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process::sibyll {
typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
class SibStackData {
public:
void Init();
void Dump() const {}
void Clear() { s_plist_.np = 0; }
unsigned int GetSize() const { return s_plist_.np; }
unsigned int GetCapacity() const { return 8000; }
void SetId(const unsigned int i, const int v) { s_plist_.llist[i] = v; }
void SetEnergy(const unsigned int i, const corsika::units::si::HEPEnergyType v) {
using namespace corsika::units::si;
s_plist_.p[3][i] = v / 1_GeV;
}
void SetMass(const unsigned int i, const corsika::units::si::HEPMassType v) {
using namespace corsika::units::si;
s_plist_.p[4][i] = v / 1_GeV;
}
void SetMomentum(const unsigned int i, const MomentumVector& v) {
using namespace corsika::units::si;
auto tmp = v.GetComponents();
for (int idx = 0; idx < 3; ++idx) s_plist_.p[idx][i] = tmp[idx] / 1_GeV;
}
int GetId(const unsigned int i) const { return s_plist_.llist[i]; }
corsika::units::si::HEPEnergyType GetEnergy(const int i) const {
using namespace corsika::units::si;
return s_plist_.p[3][i] * 1_GeV;
}
corsika::units::si::HEPEnergyType GetMass(const unsigned int i) const {
using namespace corsika::units::si;
return s_plist_.p[4][i] * 1_GeV;
}
MomentumVector GetMomentum(const unsigned int i) const {
using corsika::geometry::CoordinateSystem;
using corsika::geometry::QuantityVector;
using corsika::geometry::RootCoordinateSystem;
using namespace corsika::units::si;
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
QuantityVector<hepmomentum_d> components = {
s_plist_.p[0][i] * 1_GeV, s_plist_.p[1][i] * 1_GeV, s_plist_.p[2][i] * 1_GeV};
return MomentumVector(rootCS, components);
}
void Copy(const unsigned int i1, const unsigned int i2) {
s_plist_.llist[i2] = s_plist_.llist[i1];
for (unsigned int i = 0; i < 5; ++i) s_plist_.p[i][i2] = s_plist_.p[i][i1];
}
void Swap(const unsigned int i1, const unsigned int i2) {
std::swap(s_plist_.llist[i1], s_plist_.llist[i2]);
for (unsigned int i = 0; i < 5; ++i)
std::swap(s_plist_.p[i][i1], s_plist_.p[i][i2]);
}
void IncrementSize() { s_plist_.np++; }
void DecrementSize() {
if (s_plist_.np > 0) { s_plist_.np--; }
}
};
template <typename StackIteratorInterface>
class ParticleInterface : public corsika::stack::ParticleBase<StackIteratorInterface> {
using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
public:
void SetParticleData(const int vID, // corsika::process::sibyll::SibyllCode vID,
const corsika::units::si::HEPEnergyType vE,
const MomentumVector& vP,
const corsika::units::si::HEPMassType vM) {
SetPID(vID);
SetEnergy(vE);
SetMomentum(vP);
SetMass(vM);
}
void SetParticleData(ParticleInterface<StackIteratorInterface>& /*parent*/,
const int vID, // corsika::process::sibyll::SibyllCode vID,
const corsika::units::si::HEPEnergyType vE,
const MomentumVector& vP,
const corsika::units::si::HEPMassType vM) {
SetPID(vID);
SetEnergy(vE);
SetMomentum(vP);
SetMass(vM);
}
void SetEnergy(const corsika::units::si::HEPEnergyType v) {
GetStackData().SetEnergy(GetIndex(), v);
}
corsika::units::si::HEPEnergyType GetEnergy() const {
return GetStackData().GetEnergy(GetIndex());
}
bool HasDecayed() const { return abs(GetStackData().GetId(GetIndex())) > 100; }
void SetMass(const corsika::units::si::HEPMassType v) {
GetStackData().SetMass(GetIndex(), v);
}
corsika::units::si::HEPEnergyType GetMass() const {
return GetStackData().GetMass(GetIndex());
}
void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
corsika::process::sibyll::SibyllCode GetPID() const {
return static_cast<corsika::process::sibyll::SibyllCode>(
GetStackData().GetId(GetIndex()));
}
MomentumVector GetMomentum() const { return GetStackData().GetMomentum(GetIndex()); }
void SetMomentum(const MomentumVector& v) {
GetStackData().SetMomentum(GetIndex(), v);
}
};
typedef corsika::stack::Stack<SibStackData, ParticleInterface> SibStack;
} // end namespace corsika::process::sibyll
#endif
Processes/Sibyll/rndm_dbl.f deleted 100644 → 0
View file @ 036d63e8
C***********************************************************************
C
C interface to PHOJET double precision random number generator
C for SIBYLL \FR'14
C
C***********************************************************************
DOUBLE PRECISION FUNCTION S_RNDM(IDUMMY)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DUMMY = dble(IDUMMY)
S_RNDM= PHO_RNDM(DUMMY)
END
C***********************************************************************
C
C initialization routine for double precision random number generator
C calls PHO_RNDIN \FR'14
C
C***********************************************************************
SUBROUTINE RND_INI
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /RNDMGAS/ ISET
ISET = 0
CALL PHO_RNDIN(12,34,56,78)
END
DOUBLE PRECISION FUNCTION GASDEV(Idum)
C***********************************************************************
C Gaussian deviation
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IMPLICIT INTEGER(I-N)
COMMON /RNDMGAS/ ISET
SAVE
DATA ISET/0/
gasdev=idum
IF (ISET.EQ.0) THEN
1 V1=2.D0*S_RNDM(0)-1.D0
V2=2.D0*S_RNDM(1)-1.D0
R=V1**2+V2**2
IF(R.GE.1.D0)GO TO 1
FAC=SQRT(-2.D0*LOG(R)/R)
GSET=V1*FAC
GASDEV=V2*FAC
ISET=1
ELSE
GASDEV=GSET
ISET=0
ENDIF
RETURN
END
C***********************************************************************
DOUBLE PRECISION FUNCTION PHO_RNDM(DUMMY)
C***********************************************************************
C
C random number generator
C
C initialization by call to PHO_RNDIN needed!
C
C the algorithm is taken from
C G.Marsaglia, A.Zaman: 'Toward a unversal random number generator'
C Florida State Univ. preprint FSU-SCRI-87-70
C
C implementation by K. Hahn (Dec. 88), changed to include possibility
C of saving / reading generator registers to / from file (R.E. 10/98)
C
C generator should not depend on the hardware (if a real has
C at least 24 significant bits in internal representation),
C the period is about 2**144,
C
C internal registers:
C U(97),C,CD,CM,I,J - seed values as initialized in PHO_RNDIN
C
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
COMMON /PORAND/ U(97),C,CD,CM,I,J
100 CONTINUE
RNDMI = DUMMY
RNDMI = U(I)-U(J)
IF ( RNDMI.LT.0.D0 ) RNDMI = RNDMI+1.D0
U(I) = RNDMI
I = I-1
IF ( I.EQ.0 ) I = 97
J = J-1
IF ( J.EQ.0 ) J = 97
C = C-CD
IF ( C.LT.0.D0 ) C = C+CM
RNDMI = RNDMI-C
IF ( RNDMI.LT.0.D0 ) RNDMI = RNDMI+1.D0
IF((ABS(RNDMI).LT.0.D0).OR.(ABS(RNDMI-1.D0).LT.1.D-10)) GOTO 100
PHO_RNDM = RNDMI
END
CDECK ID>, PHO_RNDIN
SUBROUTINE PHO_RNDIN(NA1,NA2,NA3,NB1)
C***********************************************************************
C
C initialization of PHO_RNDM, has to be called before using PHO_RNDM
C
C input:
C NA1,NA2,NA3,NB1 - values for initializing the generator
C NA? must be in 1..178 and not all 1;
C 12,34,56 are the standard values
C NB1 must be in 1..168;
C 78 is the standard value
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
COMMON /PORAND/ U(97),C,CD,CM,I,J
MA1 = NA1
MA2 = NA2
MA3 = NA3
MB1 = NB1
I = 97
J = 33
DO 20 II2 = 1,97
S = 0.D0
T = 0.5D0
DO 10 II1 = 1,24
MAT = MOD(MOD(MA1*MA2,179)*MA3,179)
MA1 = MA2
MA2 = MA3
MA3 = MAT
MB1 = MOD(53*MB1+1,169)
IF ( MOD(MB1*MAT,64).GE.32 ) S = S+T
T = 0.5D0*T
10 CONTINUE
U(II2) = S
20 CONTINUE
C = 362436.D0/16777216.D0
CD = 7654321.D0/16777216.D0
CM = 16777213.D0/16777216.D0
END
CDECK ID>, PHO_RNDSI
SUBROUTINE PHO_RNDSI(UIN,CIN,CDIN,CMIN,IIN,JIN)
C***********************************************************************
C
C updates internal random number generator registers using
C registers given as arguments
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
DIMENSION UIN(97)
COMMON /PORAND/ U(97),C,CD,CM,I,J
DO 10 KKK = 1,97
U(KKK) = UIN(KKK)
10 CONTINUE
C = CIN
CD = CDIN
CM = CMIN
I = IIN
J = JIN
END
CDECK ID>, PHO_RNDSO
SUBROUTINE PHO_RNDSO(UOUT,COUT,CDOUT,CMOUT,IOUT,JOUT)
C***********************************************************************
C
C copies internal registers from randon number generator
C to arguments
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
DIMENSION UOUT(97)
COMMON /PORAND/ U(97),C,CD,CM,I,J
DO 10 KKK = 1,97
UOUT(KKK) = U(KKK)
10 CONTINUE
COUT = C
CDOUT = CD
CMOUT = CM
IOUT = I
JOUT = J
END
CDECK ID>, PHO_RNDTE
SUBROUTINE PHO_RNDTE(IO)
C***********************************************************************
C
C test of random number generator PHO_RNDM
C
C input:
C IO defines output
C 0 output only if an error is detected
C 1 output independend on an error
C
C uses PHO_RNDSI and PHO_RNDSO to bring the random number generator
C to same status as it had before the test run
C
C***********************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
SAVE
C input/output channels
INTEGER LI,LO
COMMON /POINOU/ LI,LO
DIMENSION UU(97)
DIMENSION U(6),X(6),D(6)
DATA U / 6533892.D0 , 14220222.D0 , 7275067.D0 ,
& 6172232.D0 , 8354498.D0 , 10633180.D0 /
CALL PHO_RNDSO(UU,CC,CCD,CCM,II,JJ)
CALL PHO_RNDIN(12,34,56,78)
DO 10 II1 = 1,20000
XX = PHO_RNDM(SD)
10 CONTINUE
SD = 0.D0
DO 20 II2 = 1,6
X(II2) = 4096.D0*(4096.D0*PHO_RNDM(XX))
D(II2) = X(II2)-U(II2)
SD = SD+ABS(D(II2))
20 CONTINUE
CALL PHO_RNDSI(UU,CC,CCD,CCM,II,JJ)
IF ((IO.EQ.1).OR.(ABS(SD).GT.0.D-10)) THEN
WRITE(LO,50) (U(I),X(I),D(I),I=1,6)
ENDIF
50 FORMAT(/,' PHO_RNDTE: test of the random number generator:',/,
& ' expected value calculated value difference',/,
& 6(F17.1,F20.1,F15.3,/),
& ' generator has the same status as before calling PHO_RNDTE',/)
END
CDECK ID>, PHO_RNDST
SUBROUTINE PHO_RNDST(MODE,FILENA)
C***********************************************************************
C
C read / write random number generator status from / to file
C
C input: MODE 1 read registers from file
C 2 dump registers to file
C
C FILENA file name
C
C***********************************************************************
IMPLICIT NONE
SAVE
INTEGER MODE
CHARACTER*(*) FILENA
C input/output channels
INTEGER LI,LO
COMMON /POINOU/ LI,LO
DOUBLE PRECISION UU,CC,CCD,CCM
DIMENSION UU(97)
INTEGER I,II,JJ
CHARACTER*80 CH_DUMMY
IF(MODE.EQ.1) THEN
WRITE(LO,'(/,1X,2A,A,/)') 'PHO_RNDST: ',
& 'reading random number registers from file ',FILENA
OPEN(12,FILE=FILENA,ERR=1010,STATUS='OLD')
READ(12,*,ERR=1010) CH_DUMMY
DO I=1,97
READ(12,*,ERR=1010) UU(I)
ENDDO
READ(12,*,ERR=1010) CC
READ(12,*,ERR=1010) CCD
READ(12,*,ERR=1010) CCM
READ(12,*,ERR=1010) II,JJ
CLOSE(12)
CALL PHO_RNDSI(UU,CC,CCD,CCM,II,JJ)
ELSE IF(MODE.EQ.2) THEN
WRITE(LO,'(/,1X,2A,A,/)') 'PHO_RNDST: ',
& 'dumping random number registers to file ',FILENA
OPEN(12,FILE=FILENA,ERR=1010,STATUS='UNKNOWN')
CALL PHO_RNDSO(UU,CC,CCD,CCM,II,JJ)
WRITE(12,'(1X,A)',ERR=1020) 'random number status registers:'
DO I=1,97
WRITE(12,'(1X,1P,E28.20)',ERR=1020) UU(I)
ENDDO
WRITE(12,'(1X,1P,E28.20)',ERR=1020) CC
WRITE(12,'(1X,1P,E28.20)',ERR=1020) CCD
WRITE(12,'(1X,1P,E28.20)',ERR=1020) CCM
WRITE(12,'(1X,2I4)',ERR=1020) II,JJ
CLOSE(12)
ELSE
WRITE(LO,'(/,1X,2A,I6,/)') 'PHO_RNDST: ',
& 'called with invalid mode, nothing done (mode)',MODE
ENDIF
RETURN
1010 CONTINUE
WRITE(LO,'(1X,2A,A,/)') 'PHO_RNDST: ',
& 'cannot open or read file ',FILENA
RETURN
1020 CONTINUE
WRITE(LO,'(1X,A,A,/)') 'PHO_RNDST: ',
& 'cannot open or write file ',FILENA
RETURN
END
C----------------------------------------
C standard generator
C----------------------------------------
REAL FUNCTION S_RNDM_std(IDUMMY)
C...Generator from the LUND montecarlo
C...Purpose: to generate random numbers uniformly distributed between
C...0 and 1, excluding the endpoints.
COMMON/LUDATR/MRLU(6),RRLU(100)
SAVE /LUDATR/
EQUIVALENCE (MRLU1,MRLU(1)),(MRLU2,MRLU(2)),(MRLU3,MRLU(3)),
&(MRLU4,MRLU(4)),(MRLU5,MRLU(5)),(MRLU6,MRLU(6)),
&(RRLU98,RRLU(98)),(RRLU99,RRLU(99)),(RRLU00,RRLU(100))
C... Initialize generation from given seed.
S_RNDM_std = real(idummy)
IF(MRLU2.EQ.0) THEN
IF (MRLU1 .EQ. 0) MRLU1 = 19780503 ! initial seed
IJ=MOD(MRLU1/30082,31329)
KL=MOD(MRLU1,30082)
I=MOD(IJ/177,177)+2
J=MOD(IJ,177)+2
K=MOD(KL/169,178)+1
L=MOD(KL,169)
DO 110 II=1,97
S=0.
T=0.5
DO 100 JJ=1,24
M=MOD(MOD(I*J,179)*K,179)
I=J
J=K
K=M
L=MOD(53*L+1,169)
IF(MOD(L*M,64).GE.32) S=S+T
T=0.5*T
100 CONTINUE
RRLU(II)=S
110 CONTINUE
TWOM24=1.
DO 120 I24=1,24
TWOM24=0.5*TWOM24
120 CONTINUE
RRLU98=362436.*TWOM24
RRLU99=7654321.*TWOM24
RRLU00=16777213.*TWOM24
MRLU2=1
MRLU3=0
MRLU4=97
MRLU5=33
ENDIF
C...Generate next random number.
130 RUNI=RRLU(MRLU4)-RRLU(MRLU5)
IF(RUNI.LT.0.) RUNI=RUNI+1.
RRLU(MRLU4)=RUNI
MRLU4=MRLU4-1
IF(MRLU4.EQ.0) MRLU4=97
MRLU5=MRLU5-1
IF(MRLU5.EQ.0) MRLU5=97
RRLU98=RRLU98-RRLU99
IF(RRLU98.LT.0.) RRLU98=RRLU98+RRLU00
RUNI=RUNI-RRLU98
IF(RUNI.LT.0.) RUNI=RUNI+1.
IF(RUNI.LE.0.OR.RUNI.GE.1.) GOTO 130
C...Update counters. Random number to output.
MRLU3=MRLU3+1
IF(MRLU3.EQ.1000000000) THEN
MRLU2=MRLU2+1
MRLU3=0
ENDIF
S_RNDM_std=RUNI
END
Processes/Sibyll/sibyll2.3c.cc deleted 100644 → 0
View file @ 036d63e8
/*
* (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/process/sibyll/sibyll2.3c.h>
#include <corsika/random/RNGManager.h>
#include <random>
int get_nwounded() { return s_chist_.nwd; }
double get_sibyll_mass2(int& id) { return s_mass1_.am2[id]; }
double s_rndm_(int&) {
static corsika::random::RNG& rng =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
std::uniform_real_distribution<double> dist;
return dist(rng);
}
Processes/Sibyll/sibyll_codes.dat deleted 100644 → 0
View file @ 036d63e8
# input file for particle conversion to/from SIBYLL
# the format of this file is: "corsika-identifier" "sibyll-id" "can-interact-in-sibyll" "cross-section-type"
Electron 3 0 0
Positron 2 0 0
NuE 15 0 0
NuEBar 16 0 0
MuMinus 5 0 0
MuPlus 4 0 0
NuMu 17 0 0
NuMuBar 18 0 0
TauMinus 91 0 0
TauPlus 90 0 0
NuTau 92 0 0
NuTauBar 93 0 0
Gamma 1 0 0
Pi0 6 1 2
# rho0 could interact but sibyll has no cross section/interaction length. was used for gamma had int
Rho0 27 0 0
K0Long 11 1 3
PiPlus 7 1 2
PiMinus 8 1 2
RhoPlus 25 0 0
RhoMinus 26 0 0
Eta 23 0 0
Omega 32 0 0
K0Short 12 1 3
KStar0 30 0 0
KStar0Bar 31 0 0
KPlus 9 1 3
KMinus 10 1 3
KStarPlus 28 0 0
KStarMinus 29 0 0
DPlus 59 1 0
DMinus 60 1 0
DStarPlus 78 0 0
DStarMinus 79 0 0
D0 71 1 0
D0Bar 72 1 0
DStar0 80 0 0
DStar0Bar 81 0 0
DsPlus 74 1 0
DsMinus 75 1 0
DStarSPlus 76 0 0
DStarSMinus 77 0 0
EtaC 73 0 0
Neutron 14 1 1
AntiNeutron -14 1 1
Delta0 42 0 0
Delta0Bar -42 0 0
Proton 13 1 1
AntiProton -13 1 1
DeltaPlus 41 0 0
DeltaMinusBar -41 0 0
DeltaPlusPlus 40 0 0
DeltaMinusMinusBar -40 0 0
SigmaMinus 36 1 1
SigmaPlusBar -36 1 1
Lambda0 39 1 1
Lambda0Bar -39 1 1
Sigma0 35 1 1
Sigma0Bar -35 1 1
SigmaPlus 34 1 1
SigmaMinusBar -34 1 1
XiMinus 38 1 1
XiPlusBar -38 1 1
Xi0 37 1 1
Xi0Bar -37 1 1
OmegaMinus 49 0 0
OmegaPlusBar -49 0 0
SigmaC0 86 0 0
SigmaC0Bar -86 0 0
SigmaStarC0 96 0 0
SigmaStarC0Bar -96 0 0
LambdaCPlus 89 1 1
LambdaCMinusBar -89 1 1
XiC0 88 1 1
XiC0Bar -88 1 1
SigmaCPlus 85 0 0
SigmaCMinusBar -85 0 0
SigmaStarCPlus 95 0 0
SigmaStarCMinusBar -95 0 0
SigmaCPlusPlus 84 0 0
SigmaCMinusMinusBar -84 0 0
SigmaStarCPlusPlus 94 0 0
SigmaStarCMinusMinusBar -94 0 0
XiCPlus 87 1 1
XiCMinusBar -87 1 1
OmegaC0 99 0 0
OmegaC0Bar -99 0 0
Jpsi 83 0 0
#Unknown 0 0 0
Processes/Sibyll/testSibyll.cc deleted 100644 → 0
View file @ 036d63e8
/*
* (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/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/random/RNGManager.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::process::sibyll;
TEST_CASE("Sibyll", "[processes]") {
SECTION("Sibyll -> Corsika") {
REQUIRE(corsika::particles::Electron::GetCode() ==
process::sibyll::ConvertFromSibyll(process::sibyll::SibyllCode::Electron));
}
SECTION("Corsika -> Sibyll") {
REQUIRE(process::sibyll::ConvertToSibyll(corsika::particles::Electron::GetCode()) ==
process::sibyll::SibyllCode::Electron);
REQUIRE(process::sibyll::ConvertToSibyllRaw(corsika::particles::Proton::GetCode()) ==
13);
}
SECTION("KnownBySibyll") {
REQUIRE(process::sibyll::KnownBySibyll(corsika::particles::Electron::GetCode()));
REQUIRE_FALSE(
process::sibyll::KnownBySibyll(corsika::particles::XiPrimeC0::GetCode()));
}
SECTION("canInteractInSibyll") {
REQUIRE(process::sibyll::CanInteract(corsika::particles::Proton::GetCode()));
REQUIRE(process::sibyll::CanInteract(corsika::particles::Code::XiCPlus));
REQUIRE_FALSE(process::sibyll::CanInteract(corsika::particles::Electron::GetCode()));
REQUIRE_FALSE(process::sibyll::CanInteract(corsika::particles::SigmaC0::GetCode()));
}
SECTION("cross-section type") {
REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::Electron) == 0);
REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::K0Long) == 3);
REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::SigmaPlus) == 1);
REQUIRE(process::sibyll::GetSibyllXSCode(corsika::particles::Code::PiMinus) == 2);
}
}
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
using namespace corsika::units::si;
using namespace corsika::units;
TEST_CASE("SibyllInterface", "[processes]") {
// setup environment, geometry
corsika::environment::Environment env;
auto& universe = *(env.GetUniverse());
auto theMedium = corsika::environment::Environment::CreateNode<geometry::Sphere>(
geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel =
corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
corsika::environment::NuclearComposition(
std::vector<corsika::particles::Code>{corsika::particles::Code::Oxygen},
std::vector<float>{1.}));
universe.AddChild(std::move(theMedium));
const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
geometry::Point const origin(cs, {0_m, 0_m, 0_m});
geometry::Vector<corsika::units::si::SpeedType::dimension_type> v(
cs, 0_m / second, 0_m / second, 1_m / second);
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
corsika::random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
SECTION("InteractionInterface") {
setup::Stack stack;
const HEPEnergyType E0 = 100_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(particles::Code::Proton, E0, plab, pos, 0_ns);
Interaction model(env);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
model.DoInteraction(particle, stack);
[[maybe_unused]] const GrammageType length =
model.GetInteractionLength(particle, track);
}
SECTION("NuclearInteractionInterface") {
setup::Stack stack;
const HEPEnergyType E0 = 400_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2);
Interaction hmodel(env);
NuclearInteraction model(env, hmodel);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
model.DoInteraction(particle, stack);
[[maybe_unused]] const GrammageType length =
model.GetInteractionLength(particle, track);
}
SECTION("DecayInterface") {
setup::Stack stack;
const HEPEnergyType E0 = 10_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = stack::super_stupid::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(particles::Code::Proton, E0, plab, pos, 0_ns);
Decay model;
model.Init();
/*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(particle,
stack);
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
}
}
Processes/StackInspector/CMakeLists.txt deleted 100644 → 0
View file @ 036d63e8
set (
MODEL_SOURCES
StackInspector.cc
)
set (
MODEL_HEADERS
StackInspector.h
)
set (
MODEL_NAMESPACE
corsika/process/stack_inspector
)
add_library (ProcessStackInspector STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessStackInspector ${MODEL_NAMESPACE} ${MODEL_HEADERS})
set_target_properties (
ProcessStackInspector
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
# PUBLIC_HEADER "${MODEL_HEADERS}"
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessStackInspector
CORSIKAunits
CORSIKAgeometry
CORSIKAsetup
CORSIKAlogging
)
target_include_directories (
ProcessStackInspector
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
TARGETS ProcessStackInspector
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
)
# --------------------
# code unit testing
add_executable (testStackInspector testStackInspector.cc)
target_link_libraries (
testStackInspector
ProcessStackInspector
CORSIKAsetup
CORSIKAgeometry
CORSIKAunits
CORSIKAthirdparty # for catch2
)
CORSIKA_ADD_TEST(testStackInspector)
Processes/StackInspector/StackInspector.cc deleted 100644 → 0
View file @ 036d63e8
/*
* (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/geometry/RootCoordinateSystem.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/logging/Logger.h>
#include <corsika/setup/SetupTrajectory.h>
#include <iostream>
#include <limits>
using namespace std;
using namespace corsika;
using namespace corsika::particles;
using namespace corsika::units::si;
using namespace corsika::process::stack_inspector;
template <typename Stack>
StackInspector<Stack>::StackInspector(const bool aReport)
: fReport(aReport)
, fCountStep(0) {}
template <typename Stack>
StackInspector<Stack>::~StackInspector() {}
template <typename Stack>
process::EProcessReturn StackInspector<Stack>::DoContinuous(Particle&, setup::Trajectory&,
Stack& s) {
if (!fReport) return process::EProcessReturn::eOk;
[[maybe_unused]] int i = 0;
HEPEnergyType Etot = 0_GeV;
for (auto& iterP : s) {
HEPEnergyType E = iterP.GetEnergy();
Etot += E;
geometry::CoordinateSystem& rootCS = geometry::RootCoordinateSystem::GetInstance()
.GetRootCoordinateSystem(); // for printout
auto pos = iterP.GetPosition().GetCoordinates(rootCS);
cout << "StackInspector: i=" << setw(5) << fixed << (i++) << ", id=" << setw(30)
<< iterP.GetPID() << " E=" << setw(15) << scientific << (E / 1_GeV) << " GeV, "
<< " pos=" << pos;
// if (iterP.GetPID()==Code::Nucleus)
// cout << " nuc_ref=" << iterP.GetNucleusRef();
cout << endl;
}
fCountStep++;
cout << "StackInspector: nStep=" << fCountStep << " stackSize=" << s.GetSize()
<< " Estack=" << Etot / 1_GeV << " GeV" << endl;
return process::EProcessReturn::eOk;
}
template <typename Stack>
corsika::units::si::LengthType StackInspector<Stack>::MaxStepLength(Particle&,
setup::Trajectory&) {
return std::numeric_limits<double>::infinity() * meter;
}
template <typename Stack>
void StackInspector<Stack>::Init() {
fCountStep = 0;
}
#include <corsika/setup/SetupStack.h>
template class process::stack_inspector::StackInspector<setup::Stack>;
Processes/StackInspector/StackInspector.h deleted 100644 → 0
View file @ 036d63e8
/*
* (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.
*/
#ifndef _Physics_StackInspector_StackInspector_h_
#define _Physics_StackInspector_StackInspector_h_
#include <corsika/process/ContinuousProcess.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process {
namespace stack_inspector {
template <typename Stack>
class StackInspector
: public corsika::process::ContinuousProcess<StackInspector<Stack>> {
typedef typename Stack::ParticleType Particle;
public:
StackInspector(const bool aReport);
~StackInspector();
void Init();
EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack& s);
//~ template <typename Particle>
corsika::units::si::LengthType MaxStepLength(Particle&,
corsika::setup::Trajectory&);
private:
bool fReport;
int fCountStep = 0;
};
} // namespace stack_inspector
} // namespace corsika::process
#endif
Processes/StackInspector/testStackInspector.cc deleted 100644 → 0
View file @ 036d63e8
/*
* (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.
*/
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
#include <corsika/process/stack_inspector/StackInspector.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
using namespace corsika::units::si;
using namespace corsika::process::stack_inspector;
using namespace corsika;
TEST_CASE("StackInspector", "[processes]") {
auto const& rootCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
geometry::Point const origin(rootCS, {0_m, 0_m, 0_m});
geometry::Vector<corsika::units::si::SpeedType::dimension_type> v(
rootCS, 0_m / second, 0_m / second, 1_m / second);
geometry::Line line(origin, v);
geometry::Trajectory<geometry::Line> track(line, 10_s);
setup::Stack stack;
auto particle = stack.AddParticle(
particles::Code::Electron, 10_GeV,
corsika::stack::super_stupid::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
geometry::Point(rootCS, {0_m, 0_m, 10_km}), 0_ns);
SECTION("interface") {
StackInspector<setup::Stack> model(true);
model.Init();
[[maybe_unused]] const process::EProcessReturn ret =
model.DoContinuous(particle, track, stack);
[[maybe_unused]] const LengthType length = model.MaxStepLength(particle, track);
}
}
Processes/TrackWriter/CMakeLists.txt deleted 100644 → 0
View file @ 036d63e8
set (
MODEL_SOURCES
TrackWriter.cc
)
set (
MODEL_HEADERS
TrackWriter.h
)
set (
MODEL_NAMESPACE
corsika/process/track_writer
)
add_library (ProcessTrackWriter STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackWriter ${MODEL_NAMESPACE} ${MODEL_HEADERS})
set_target_properties (
ProcessTrackWriter
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
# PUBLIC_HEADER "${MODEL_HEADERS}"
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessTrackWriter
CORSIKAunits
CORSIKAparticles
CORSIKAgeometry
CORSIKAsetup
)
target_include_directories (
ProcessTrackWriter
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
TARGETS ProcessTrackWriter
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE}
)
# --------------------
# code unit testing
#add_executable (testNullModel testNullModel.cc)
#target_link_libraries (
# testNullModel ProcessNullModel
# CORSIKAsetup
# CORSIKAgeometry
# CORSIKAunits
# CORSIKAthirdparty # for catch2
# )
# CORSIKA_ADD_TEST(testNullModel)
Processes/TrackWriter/TrackWriter.cc deleted 100644 → 0
View file @ 036d63e8
/*
* (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/process/track_writer/TrackWriter.h>
#include <string>
void corsika::process::TrackWriter::TrackWriter::Init() {
using namespace std::string_literals;
fFile.open(fFilename);
fFile << "# PID, E / eV, start coordinates / m, displacement vector to end / m "s
<< '\n';
}
Processes/TrackWriter/TrackWriter.h deleted 100644 → 0
View file @ 036d63e8
/*
* (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.
*/
#ifndef _Processes_TrackWriter_h_
#define _Processes_TrackWriter_h_
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <fstream>
#include <limits>
#include <string>
namespace corsika::process::TrackWriter {
class TrackWriter : public corsika::process::ContinuousProcess<TrackWriter> {
public:
TrackWriter(std::string const& filename)
: fFilename(filename) {}
void Init();
template <typename Particle, typename Track, typename Stack>
corsika::process::EProcessReturn DoContinuous(Particle& p, Track& t, Stack&) {
using namespace corsika::units::si;
auto const start = t.GetPosition(0).GetCoordinates();
auto const delta = t.GetPosition(1).GetCoordinates() - start;
auto const& name = corsika::particles::GetName(p.GetPID());
fFile << name << " " << p.GetEnergy() / 1_eV << ' ' << start[0] / 1_m << ' '
<< start[1] / 1_m << ' ' << start[2] / 1_m << " " << delta[0] / 1_m << ' '
<< delta[1] / 1_m << ' ' << delta[2] / 1_m << '\n';
return corsika::process::EProcessReturn::eOk;
}
template <typename Particle, typename Track>
corsika::units::si::LengthType MaxStepLength(Particle&, Track&) {
return corsika::units::si::meter * std::numeric_limits<double>::infinity();
}
private:
std::string const fFilename;
std::ofstream fFile;
};
} // namespace corsika::process::TrackWriter
#endif
Processes/TrackingLine/CMakeLists.txt deleted 100644 → 0
View file @ 036d63e8
set (
MODEL_HEADERS
TrackingLine.h
)
set (
MODEL_NAMESPACE
corsika/process/tracking_line
)
add_library (ProcessTrackingLine INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessTrackingLine ${MODEL_NAMESPACE} ${MODEL_HEADERS})
target_include_directories (
ProcessTrackingLine
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
# #-- -- -- -- -- -- -- -- -- --
# #code unit testing
add_executable (testTrackingLine testTrackingLine.cc)
target_link_libraries (
testTrackingLine
ProcessTrackingLine
CORSIKAsetup
CORSIKAutilities
CORSIKAenvironment
CORSIKAunits
CORSIKAenvironment
CORSIKAgeometry
CORSIKAthirdparty # for catch2
)
CORSIKA_ADD_TEST(testTrackingLine)
Processes/TrackingLine/TrackingLine.h deleted 100644 → 0
View file @ 036d63e8
/*
* (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.
*/
#ifndef _include_corsika_processes_TrackingLine_h_
#define _include_corsika_processes_TrackingLine_h_
#include <corsika/geometry/Point.h>
#include <corsika/geometry/QuantityVector.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/environment/Environment.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <algorithm>
#include <iostream>
#include <optional>
#include <stdexcept>
#include <utility>
namespace corsika::process {
namespace tracking_line {
template <typename Stack>
class TrackingLine { //
using Particle = typename Stack::ParticleType;
corsika::environment::Environment const& fEnvironment;
public:
std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
TimeOfIntersection(corsika::geometry::Line const& line,
geometry::Sphere const& sphere) {
using namespace corsika::units::si;
auto const& cs = fEnvironment.GetCoordinateSystem();
geometry::Point const origin(cs, 0_m, 0_m, 0_m);
auto const r0 = (line.GetR0() - origin);
auto const v0 = line.GetV0();
auto const c0 = (sphere.GetCenter() - origin);
auto const alpha = r0.dot(v0) - 2 * v0.dot(c0);
auto const beta = c0.squaredNorm() + r0.squaredNorm() + 2 * c0.dot(r0) -
sphere.GetRadius() * sphere.GetRadius();
auto const discriminant = alpha * alpha - 4 * beta * v0.squaredNorm();
//~ std::cout << "discriminant: " << discriminant << std::endl;
//~ std::cout << "alpha: " << alpha << std::endl;
//~ std::cout << "beta: " << beta << std::endl;
if (discriminant.magnitude() > 0) {
(-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm());
return std::make_pair((-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm()),
(-alpha + sqrt(discriminant)) / (2 * v0.squaredNorm()));
} else {
return {};
}
}
TrackingLine(corsika::environment::Environment const& pEnv)
: fEnvironment(pEnv) {}
void Init() {}
auto GetTrack(Particle const& p) {
using std::cout;
using std::endl;
using namespace corsika::units::si;
using namespace corsika::geometry;
geometry::Vector<SpeedType::dimension_type> const velocity =
p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
auto const currentPosition = p.GetPosition();
std::cout << "TrackingLine pid: " << p.GetPID()
<< " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
std::cout << "TrackingLine pos: " << currentPosition.GetCoordinates()
<< std::endl;
std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
<< " GeV " << std::endl;
std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
// to do: include effect of magnetic field
geometry::Line line(currentPosition, velocity);
auto const* currentVolumeNode =
fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
auto const& children = currentVolumeNode->GetChildNodes();
auto const& excluded = currentVolumeNode->GetExcludedNodes();
std::vector<TimeType> intersectionTimes;
auto addIfIntersects = [&](auto& vtn) {
auto const& volume = vtn.GetVolume();
auto const& sphere = dynamic_cast<geometry::Sphere const&>(
volume); // for the moment we are a bit bold here and assume
// everything is a sphere, crashes with exception if not
if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
auto const [t1, t2] = *opt;
if (t1.magnitude() >= 0)
intersectionTimes.push_back(t1);
else if (t2.magnitude() >= 0)
intersectionTimes.push_back(t2);
}
};
for (auto const& child : children) { addIfIntersects(*child); }
for (auto const* child : excluded) { addIfIntersects(*child); }
addIfIntersects(*currentVolumeNode);
auto const minIter =
std::min_element(intersectionTimes.cbegin(), intersectionTimes.cend());
TimeType min;
if (minIter == intersectionTimes.cend()) {
min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
// to handle the numerics properly
//~ throw std::runtime_error("no intersection with anything!");
} else {
min = *minIter;
}
std::cout << " t-intersect: " << min << std::endl;
return geometry::Trajectory<corsika::geometry::Line>(line, min);
}
};
} // namespace tracking_line
} // namespace corsika::process
#endif
Processes/TrackingLine/testTrackingLine.cc deleted 100644 → 0
View file @ 036d63e8
/*
* (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/environment/Environment.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/tracking_line/TrackingLine.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Sphere.h>
#include <corsika/geometry/Vector.h>
#include <corsika/setup/SetupTrajectory.h>
using corsika::setup::Trajectory;
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
// cpp file
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::geometry;
#include <iostream>
using namespace std;
using namespace corsika::units::si;
typedef corsika::units::si::hepmomentum_d MOMENTUM;
struct DummyParticle {
HEPEnergyType fEnergy;
Vector<MOMENTUM> fMomentum;
Point fPosition;
DummyParticle(HEPEnergyType pEnergy, Vector<MOMENTUM> pMomentum, Point pPosition)
: fEnergy(pEnergy)
, fMomentum(pMomentum)
, fPosition(pPosition) {}
auto GetEnergy() const { return fEnergy; }
auto GetMomentum() const { return fMomentum; }
auto GetPosition() const { return fPosition; }
auto GetPID() const { return corsika::particles::Code::Unknown; }
};
struct DummyStack {
using ParticleType = DummyParticle;
};
TEST_CASE("TrackingLine") {
corsika::environment::Environment env; // dummy environment
auto const& cs = env.GetCoordinateSystem();
tracking_line::TrackingLine<DummyStack> tracking(env);
SECTION("intersection with sphere") {
Point const origin(cs, {0_m, 0_m, 0_m});
Point const center(cs, {0_m, 0_m, 10_m});
Sphere const sphere(center, 1_m);
Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
0_m / second, 1_m / second);
Line line(origin, v);
geometry::Trajectory<Line> traj(line, 12345_s);
auto const opt =
tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
REQUIRE(opt.has_value());
auto [t1, t2] = opt.value();
REQUIRE(t1 / 9_s == Approx(1));
REQUIRE(t2 / 11_s == Approx(1));
auto const optNoIntersection =
tracking.TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
REQUIRE_FALSE(optNoIntersection.has_value());
}
SECTION("maximally possible propagation") {
auto& universe = *(env.GetUniverse());
//~ std::cout << env.GetUniverse().get() << std::endl;
DummyParticle p(1_GeV, Vector<MOMENTUM>(cs, 0_GeV, 0_GeV, 1_GeV),
Point(cs, 0_m, 0_m, 0_m));
auto const radius = 20_m;
auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
universe.AddChild(std::move(theMedium));
Point const origin(cs, {0_m, 0_m, 0_m});
Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m / second,
0_m / second, 1_m / second);
Line line(origin, v);
auto const traj = tracking.GetTrack(p);
REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius))
.GetComponents(cs)
.norm()
.magnitude() == Approx(0).margin(1e-4));
}
}
README.md
View file @ 00b7c3f3
# CORSIKA 8 Framework for Particle Cascades in Astroparticle Physics
# CORSIKA 8 Framework for Particle Cascades in Astroparticle Physics
The purpose of CORSIKA is to simulate any particle cascades in
astroparticle physics or astrophysical context. A lot of emphasis is
The purpose of CORSIKA 8 is to simulate any particle cascades in
astroparticle physics or astrophysical context. A lot of emphasis has been
put on modularity, flexibility, completeness, validation and
correctness. To boost computational efficiency different techniques
correctness. To boost computational efficiency, different techniques
are provided, like thinning or cascade equations. The aim is that
CORSIKA remains the most comprehensive framework for simulating
CORSIKA 8 remains the most comprehensive framework for simulating
particle cascades with stochastic and continuous processes.
The software makes extensive use of static design patterns and
compiler optimization. Thus, the most fundamental configuration
decision of the user must be performed at compile time. At run time
only specific model parameters can still be changed.
decisions of the user must be performed at compile time. At run time,
model parameters can still be changed.
CORSIKA 8 is released under the GPLv3 license. This does not exclude
that specific CORSIKA 8 versions can be released for specific purposes
under different licensing. See [license
file](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/LICENSE)
CORSIKA 8 is by default released under the BSD 3-Clause License. See [license
file](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/blob/master/LICENSE)
which is part of every release and the source code.
If you use, or want to refer to, CORSIKA 8 please cite ["Towards a Next
Generation of CORSIKA: A Framework for the Simulation of Particle
Cascades in Astroparticle Physics", Comput.Softw.Big Sci. 3 (2019)
2](https://doi.org/10.1007/s41781-018-0013-0). We kindly ask (and
expect) any relevant improvement or addition to be offered or
Cascades in Astroparticle Physics", Comput. Softw. Big Sci. 3 (2019)
2](https://doi.org/10.1007/s41781-018-0013-0) as well as
["Simulating radio emission from particle cascades with CORSIKA 8", Astropart. Phys. 166 (2025)
103072](https://doi.org/10.1016/j.astropartphys.2024.103072).
We kindly ask (and require) any relevant improvement or addition to be offered or
contributed to the main CORSIKA 8 repository for the benefit of the
whole community.
When you plan to contribute to CORSIKA 8 check the guidelines outlined here:
CORSIKA 8 makes use of various third-party code, in particular interaction
models. Please check the [using and collaborating
agreement](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/blob/master/USING_COLLABORATING.md)
for further information on this topic.
If you plan to contribute to CORSIKA 8, please check the guidelines outlined here:
[coding
guidelines](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/CONTRIBUTING.md). Code
that fails the review by the CORSIKA author group must be improved
guidelines](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/blob/master/CONTRIBUTING.md). Code
that fails the review by the CORSIKA 8 author group must be improved
before it can be merged in the official code base. After your code has
been accepted and merged you become a contributor of the CORSIKA 8
project and you should include yourself in the
[AUTHORS](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/AUTHORS)
file.
IMPORTANT: Before you contribute, you need to read and agree to the
[collaboration
agreement](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/COLLABORATION_AGREEMENT.md). The
agreement can be discussed, and eventually improved.
been accepted and merged, you become a contributor of the CORSIKA 8
project (code author).
We also want to point you to the [MCnet
guidelines](https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/blob/master/MCNET_GUIDELINES),
which are very useful also for us.
IMPORTANT: Before you contribute, you need to read and agree to the conditions set out in the
[using and collaborating
agreement](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/blob/master/USING_COLLABORATING.md).
The agreement can be discussed, and eventually improved if necessary.
## Get in contact
* Connect to https://gitlab.ikp.kit.edu; register yourself and join the "Air Shower Physics" group
* Join our chat threads using Mattermost via this [invite link](https://mattermost.hzdr.de/signup_user_complete/?id=xtdd8jyt6trbiezt71gaz3z4ge&md=link&sbr=su). Click the `GitLab` button, then `Sign in with Helmholtz ID`. You will be able to make an account by either finding your institution, or using your e.g. ORCID, GitHub, or Google account.
* Connect to https://gitlab.iap.kit.edu, register yourself and join the "Air Shower Physics" group. Write to us on Mattermost (in the User Questions channel), or directly contact one of the [steering comittee members](https://gitlab.iap.kit.edu/AirShowerPhysics/corsika/-/wikis/Steering-Committee) in case there are problems with that.
* Connect to corsika-devel@lists.kit.edu (self-register at
https://www.lists.kit.edu/sympa/subscribe/corsika-devel) to get in
touch with the project
touch with the project.
## Installation
CORSIKA 8 is tested regularly at least on gcc7.3.0 and clang-6.0.0.
Additional software prerequisites: eigen3, boost, cmake, g++, git.
On a bare Ubuntu 18.04, just add:
CORSIKA 8 is tested regularly at least on `gcc11.0.0` and `clang-14.0.0`.
### Prerequisites
You will also need:
- Python 3 (supported versions are Python >= 3.6), with pip
- cmake > 3.4
- git
- g++, gfortran, binutils, make
- optional: FLUKA (see below)
On a bare Ubuntu machine, just add:
``` shell
sudo apt-get install python3 python3-pip cmake g++ gfortran git doxygen graphviz
```
sudo apt-get install libeigen3-dev libboost-dev cmake g++ git
### Creating a virtual environment and Conan
It is recommended that you install CORSIKA 8 and its dependencies within a python3 virtual environment.
To do so, you can run the following.
``` shell
# Create the environment using your native python3 binary
python3 -m venv /path/to/new/virtual/environment/corsika-8
# Load the environment (should be run each time you open a new terminal)
source /path/to/new/virtual/environment/corsika-8/bin/activate
```
Follow these steps to download and install CORSIKA 8 milestone2
You will need to load the environment each time that you open a new terminal.
CORSIKA 8 uses the [conan](https://conan.io/) package manager to
manage our dependencies. Currently, version 2.50.0 or higher is required.
**Note**: if you are NOT using a virtual environment, you may want to use the `pip install --user` flag.
``` shell
pip install conan particle==0.25.1 numpy
```
git clone git@gitlab.ikp.kit.edu:AirShowerPhysics/corsika.git
cd corsika
git checkout milestone2
mkdir ../corsika-build
cd ../corsika-build
cmake ../corsika -DCMAKE_INSTALL_PREFIX=../corsika-install
make -j8
### Enabling FLUKA support
For legal reasons we do not distribute/bundle FLUKA together with CORSIKA 8.
As FLUKA is the standard low-energy hadronic interaction model for CORSIKA 8, you have to download
it separately from (http://www.fluka.org/), which requires registering there as FLUKA user.
The following should be done *before* compiling CORSIKA 8:
1. Note your system's version of gfortran (`gfortran --version`) and glibc (`ldd --version`)
2. Download the FLUKA __binary__, ensuring that it matches the versions you found above
3. Download the FLUKA __data file__ (will be named something similar to __fluka20xy.z-data.tar.gz__).
4. Un-tar the files that you downloaded using `tar -xf <filename>`
5. Set environmental variables `export FLUFOR=gfortran` and `export FLUPRO=<path to where you unzipped the files>`. Note that the `FLUPRO` directory should contain __libflukahp.a__. Both of these variables will have to be set every time you open a new terminal.
6. Go to the `FLUPRO` directory and run `make`. This will compile an exe, __flukahp__, in your current directory.
7. Follow the normal steps to compile CORSIKA 8 (see below).
When you later install CORSIKA 8, you should see a message during the __cmake__ step indicating the FLUKA was correctly found.
``` shell
libflukahp.a found in directory <some location here> via FLUPRO environment variable
FLUKA support is enabled.
```
### Compiling CORSIKA 8
Once Conan is installed and FLUKA provided, follow these steps to download and install CORSIKA 8:
``` shell
cd ./top/directory/for/corsika/installation
git clone --recursive git@gitlab.iap.kit.edu:AirShowerPhysics/corsika.git
# Or for https: git clone --recursive https://gitlab.iap.kit.edu/AirShowerPhysics/corsika.git
mkdir corsika-build
cd corsika-build
../corsika/conan-install.sh --source-directory ../corsika --release-with-debug
# conan-install.sh takes required options from command line to install dependencies for 'Debug', 'Release' and 'RelWithDebInfo' builds.
../corsika/corsika-cmake.sh -c "-DCMAKE_BUILD_TYPE="RelWithDebInfo" -DWITH_FLUKA=ON -DCMAKE_INSTALL_PREFIX=../corsika-install"
make -j4 #The number should match the number of available cores on your machine
make install
make test
```
and if you want to see how the first simple hadron cascade develops, see `Documentation/Examples/cascade_example.cc` for a starting point.
Run the cascade_example with:
## Alternate installation using docker containers
There are docker containers prepared that bring all the environment and packages you need to run CORSIKA. See [docker hub](https://hub.docker.com/repository/docker/corsika/devel) for a complete overview.
### Prerequisites
You only need docker, e.g. on Ubuntu: `sudo apt-get install docker` and of course root access.
### Compiling
Follow these steps to download and install CORSIKA 8, master development version
```shell
cd ./top/directory/for/corsika/installation
git clone --recursive git@gitlab.iap.kit.edu:AirShowerPhysics/corsika.git
sudo docker run -v $PWD:/corsika -it corsika/devel:clang-8 /bin/bash
mkdir corsika-build
cd corsika-build
../corsika/conan-install.sh --source-directory ../corsika --release-with-debug
# conan-install.sh takes required options from command line to install dependencies for 'Debug', 'Release' and 'RelWithDebInfo' builds.
../corsika/corsika-cmake.sh -c "-DCMAKE_BUILD_TYPE="RelWithDebInfo" -DWITH_FLUKA=ON -DCMAKE_INSTALL_PREFIX=../corsika-install"
make -j4 #The number should match the number of available cores on your machine
make install
```
cd ../corsika-install
share/examples/cascade_example
## Running Unit Tests
To run the unit tests, do the following.
```shell
cd ./corsika-build
ctest -j4 #The number should match the number of available cores on your machine
```
## Running applications and examples
### Standard applications
Applications for standard use-cases are located in the `applications` directory.
These are example scripts that can be used directly or slightly modified for your use case.
See [applications/README.md] for more.
The applications are compiled automatically after running `make` and will appear your `corsika-build/bin` directory.
After running `make install` the binaries will also be copied into your `corsika-install/bin` directory as well.
For example, from inside your `corsika-install/bin` directory, run
```shell
c8_air_shower --pdg 2212 -E 1e5 -f my_shower
```
This will run a vertical 100 TeV proton shower and will create and put the output into `./my_shower`.
Visualize output (needs gnuplot installed):
### Building the examples
Unlike the applications, the examples must be compiled as a second step.
From your top corsika directory, (the one that includes `corsika-build` and `corsika-install`) run
```shell
export CONAN_DEPENDENCIES=$PWD/corsika-install/lib/cmake/dependencies
cmake -DCMAKE_TOOLCHAIN_FILE=${CONAN_DEPENDENCIES}/conan_toolchain.cmake -DCMAKE_PREFIX_PATH=${CONAN_DEPENDENCIES} -DCMAKE_POLICY_DEFAULT_CMP0091=NEW -DCMAKE_BUILD_TYPE=RelWithDebInfo -Dcorsika_DIR=$PWD/corsika-build -DWITH_FLUKA=ON -S $PWD/corsika/examples -B $PWD/corsika-build-examples
cd corsika-build-examples
make -j4 #The number should match the number of available cores on your machine
```
bash share/tools/plot_tracks.sh tracks.dat
firefox tracks.dat.gif
You can run the examples by using the binaries in `corsika-build-examples/bin/`.
For example:
```shell
corsika-build-examples/bin/known_particles
```
This will print out all of the particles that are known by CORSIKA.
### Generating doxygen documentation
To generate the documentation, you need doxygen and graphviz. On Ubuntu 18.04, do:
```
To generate the documentation, you need doxygen and graphviz. If you work with
the docker corsika/devel containers this is already included.
Otherwise, e.g. on Ubuntu machines, do:
```shell
sudo apt-get install doxygen graphviz
```
Switch to the corsika build directory and do
```
make doxygen
Switch to the `corsika-build` directory and do
```shell
make docs
make install
```
browse with firefox:
open with firefox:
```shell
firefox ../corsika-install/share/corsika/doc/html/index.html
```
firefox ../corsika-install/share/doc/html/index.html
```

data privacy - impressum