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 0 additions and 3546 deletions
Processes/QGSJetII/CMakeLists.txt deleted 100644 → 0
View file @ cc68b005
if (NOT DEFINED ENV{CORSIKA_DATA})
message (WARNING "WARNING: QGSJetII will not work without corsika-data. Set CORSIKA_DATA path!")
endif ()
add_custom_command (
OUTPUT ${PROJECT_BINARY_DIR}/Processes/QGSJetII/Generated.inc
COMMAND ${PROJECT_SOURCE_DIR}/Processes/QGSJetII/code_generator.py
${PROJECT_BINARY_DIR}/Framework/Particles/particle_db.pkl
${PROJECT_SOURCE_DIR}/Processes/QGSJetII/qgsjet-II-04-codes.dat
DEPENDS code_generator.py
qgsjet-II-04-codes.dat
CORSIKAparticles
WORKING_DIRECTORY
${PROJECT_BINARY_DIR}/Processes/QGSJetII/
COMMENT "Generate conversion tables for particle codes QGSJetII <-> CORSIKA"
VERBATIM
)
set (
MODEL_SOURCES
ParticleConversion.cc
Interaction.cc
qgsjet-II-04.f
qgsjet-II-04.cc
)
set (
MODEL_HEADERS
ParticleConversion.h
qgsjet-II-04.h
QGSJetIIStack.h
QGSJetIIFragmentsStack.h
Interaction.h
${PROJECT_BINARY_DIR}/Processes/QGSJetII/Generated.inc
)
set (
MODEL_NAMESPACE
corsika/process/qgsjetII
)
add_library (ProcessQGSJetII STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessQGSJetII ${MODEL_NAMESPACE} ${MODEL_HEADERS})
# ....................................................
# since Generated.inc is an automatically produced file in the build directory,
# create a symbolic link into the source tree, so that it can be found and edited more easily
# this is not needed for the build to succeed! .......
add_custom_command (
OUTPUT ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc
COMMAND ${CMAKE_COMMAND} -E create_symlink ${PROJECT_BINARY_DIR}/include/corsika/process/qgsjetII/Generated.inc ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc
COMMENT "Generate link in source-dir: ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc"
)
add_custom_target (SourceDirLinkQgsII DEPENDS ${PROJECT_BINARY_DIR}/Processes/QGSJetII/Generated.inc)
add_dependencies (ProcessQGSJetII SourceDirLinkQgsII)
# .....................................................
set_target_properties (
ProcessQGSJetII
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessQGSJetII
CORSIKAprocesssequence
CORSIKAparticles
CORSIKAutilities
CORSIKAunits
CORSIKAthirdparty
CORSIKAgeometry
CORSIKAenvironment
CORSIKAsetup
CORSIKArandom
CorsikaData
)
target_include_directories (
ProcessQGSJetII
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
TARGETS ProcessQGSJetII
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
)
# --------------------
# code unit testing
CORSIKA_ADD_TEST(testQGSJetII
SOURCES
testQGSJetII.cc
${MODEL_HEADERS}
)
target_link_libraries (
testQGSJetII
ProcessQGSJetII
CORSIKAtesting
)
# also provide QGSJetII large data tables for testing in build tree (links)
# -> changed that to use environment variable "${CORSIKA_DATA}/QGSJetII"
# add_custom_command (
# OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/qgsdat-II-04 ${CMAKE_CURRENT_BINARY_DIR}/sectnu-II-04
# COMMAND ${CMAKE_COMMAND} -E create_symlink ${CMAKE_CURRENT_SOURCE_DIR}/qgsdat-II-04 ${CMAKE_CURRENT_BINARY_DIR}/qgsdat-II-04
# COMMAND ${CMAKE_COMMAND} -E create_symlink ${CMAKE_CURRENT_SOURCE_DIR}/sectnu-II-04 ${CMAKE_CURRENT_BINARY_DIR}/sectnu-II-04
# COMMENT "Generate link in source-dir: qgsjet-II-04"
# )
# add_custom_target (QGSJetDataLink DEPENDS
# ${CMAKE_CURRENT_BINARY_DIR}/qgsdat-II-04
# ${CMAKE_CURRENT_BINARY_DIR}/sectnu-II-04
# )
# add_dependencies (testQGSJetII QGSJetDataLink)
Processes/QGSJetII/Interaction.cc deleted 100644 → 0
View file @ cc68b005
/*
* (c) Copyright 2020 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/process/qgsjetII/Interaction.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/FourVector.h>
#include <corsika/geometry/QuantityVector.h>
#include <corsika/process/qgsjetII/QGSJetIIFragmentsStack.h>
#include <corsika/process/qgsjetII/QGSJetIIStack.h>
#include <corsika/process/qgsjetII/qgsjet-II-04.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/utl/COMBoost.h>
#include <random>
#include <sstream>
#include <string>
#include <tuple>
using std::cout;
using std::endl;
using std::ostringstream;
using std::string;
using std::tuple;
using namespace corsika;
using namespace corsika::setup;
using SetupParticle = setup::Stack::StackIterator;
using SetupProjectile = setup::StackView::StackIterator;
using Track = Trajectory;
namespace corsika::process::qgsjetII {
Interaction::Interaction(const string& dataPath)
: data_path_(dataPath) {
if (dataPath == "") {
if (std::getenv("CORSIKA_DATA")) {
data_path_ = string(std::getenv("CORSIKA_DATA")) + "/QGSJetII/";
cout << "Searching for QGSJetII data tables in " << data_path_ << endl;
}
}
// initialize QgsjetII
if (!initialized_) {
qgset_();
datadir DIR(data_path_);
qgaini_(DIR.data);
initialized_ = true;
}
}
Interaction::~Interaction() { cout << "QgsjetII::Interaction n=" << count_ << endl; }
units::si::CrossSectionType Interaction::GetCrossSection(
const particles::Code beamId, const particles::Code targetId,
const units::si::HEPEnergyType Elab, const unsigned int Abeam,
const unsigned int targetA) const {
using namespace units::si;
double sigProd = std::numeric_limits<double>::infinity();
if (process::qgsjetII::CanInteract(beamId)) {
const int xsCode = process::qgsjetII::GetQgsjetIIXSCodeRaw(beamId);
int iTarget = 1;
if (particles::IsNucleus(targetId)) {
iTarget = targetA;
if (iTarget > maxMassNumber_ || iTarget <= 0) {
std::ostringstream txt;
txt << "QgsjetII target outside range. iTarget=" << iTarget;
throw std::runtime_error(txt.str().c_str());
}
}
int iProjectile = 1;
if (particles::IsNucleus(beamId)) {
iProjectile = Abeam;
if (iProjectile > maxMassNumber_ || iProjectile <= 0)
throw std::runtime_error("QgsjetII target outside range. ");
}
cout << "QgsjetII::GetCrossSection Elab=" << Elab << " xs-code=" << xsCode
<< " iProjectile=" << iProjectile << " iTarget=" << iTarget << endl;
sigProd = qgsect_(Elab / 1_GeV, xsCode, iProjectile, iTarget);
cout << "QgsjetII::GetCrossSection sigProd=" << sigProd << endl;
}
return sigProd * 1_mb;
}
template <>
units::si::GrammageType Interaction::GetInteractionLength(
SetupParticle const& vP) const {
using namespace units;
using namespace units::si;
using namespace geometry;
// coordinate system, get global frame of reference
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
const particles::Code corsikaBeamId = vP.GetPID();
// beam particles for qgsjetII : 1, 2, 3 for p, pi, k
// read from cross section code table
const bool kInteraction = process::qgsjetII::CanInteract(corsikaBeamId);
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy
HEPEnergyType Elab = vP.GetEnergy();
cout << "Interaction: LambdaInt: \n"
<< " input energy: " << vP.GetEnergy() / 1_GeV << endl
<< " beam can interact:" << kInteraction << endl
<< " beam pid:" << vP.GetPID() << endl;
if (kInteraction) {
int Abeam = 0;
if (particles::IsNucleus(vP.GetPID())) Abeam = vP.GetNuclearA();
// 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..
*/
auto const* currentNode = vP.GetNode();
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
si::CrossSectionType weightedProdCrossSection = mediumComposition.WeightedSum(
[=](particles::Code targetID) -> si::CrossSectionType {
int targetA = 0;
if (corsika::particles::IsNucleus(targetID))
targetA = particles::GetNucleusA(targetID);
return GetCrossSection(corsikaBeamId, targetID, Elab, Abeam, targetA);
});
cout << "Interaction: "
<< "IntLength: weighted CrossSection (mb): " << weightedProdCrossSection / 1_mb
<< endl;
// calculate interaction length in medium
GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
units::constants::u / weightedProdCrossSection;
cout << "Interaction: "
<< "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
<< endl;
return int_length;
}
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
/**
In this function QGSJETII is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
*/
template <>
process::EProcessReturn Interaction::DoInteraction(SetupProjectile& vP) {
using namespace units;
using namespace utl;
using namespace units::si;
using namespace geometry;
const auto corsikaBeamId = vP.GetPID();
cout << "ProcessQgsjetII: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
<< process::qgsjetII::CanInteract(corsikaBeamId) << endl;
if (process::qgsjetII::CanInteract(corsikaBeamId)) {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in QgsjetII
Point pOrig = vP.GetPosition();
TimeType tOrig = vP.GetTime();
// define target
// for QgsjetII is always a single nucleon
// FOR NOW: target is always at rest
const auto targetEnergyLab = 0_GeV + constants::nucleonMass;
const auto targetMomentumLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargLab(targetEnergyLab, targetMomentumLab);
// define projectile
HEPEnergyType const projectileEnergyLab = vP.GetEnergy();
auto const projectileMomentumLab = vP.GetMomentum();
int beamA = 1;
if (particles::IsNucleus(corsikaBeamId)) beamA = vP.GetNuclearA();
const HEPEnergyType projectileEnergyLabPerNucleon = projectileEnergyLab / beamA;
cout << "Interaction: ebeam lab: " << projectileEnergyLab / 1_GeV << endl
<< "Interaction: pbeam lab: " << projectileMomentumLab.GetComponents() / 1_GeV
<< endl;
cout << "Interaction: etarget lab: " << targetEnergyLab / 1_GeV << endl
<< "Interaction: ptarget lab: " << targetMomentumLab.GetComponents() / 1_GeV
<< endl;
cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "Interaction: time: " << tOrig << endl;
// sample target mass number
auto const* currentNode = vP.GetNode();
auto const& 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
*/
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];
int targetA = 0;
if (corsika::particles::IsNucleus(targetId))
targetA = particles::GetNucleusA(targetId);
const auto sigProd =
GetCrossSection(corsikaBeamId, targetId, projectileEnergyLab, beamA, targetA);
cross_section_of_components[i] = sigProd;
}
const auto targetCode =
mediumComposition.SampleTarget(cross_section_of_components, rng_);
cout << "Interaction: target selected: " << targetCode << endl;
int targetMassNumber = 1; // proton
if (particles::IsNucleus(targetCode)) { // nucleus
targetMassNumber = particles::GetNucleusA(targetCode);
if (targetMassNumber > maxMassNumber_)
throw std::runtime_error("QgsjetII target mass outside range.");
} else {
if (targetCode != particles::Proton::GetCode())
throw std::runtime_error("QgsjetII Taget not possible.");
}
cout << "Interaction: target qgsjetII code/A: " << targetMassNumber << endl;
int projectileMassNumber = 1; // "1" means "hadron"
QgsjetIIHadronType qgsjet_hadron_type =
process::qgsjetII::GetQgsjetIIHadronType(corsikaBeamId);
if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
projectileMassNumber = vP.GetNuclearA();
if (projectileMassNumber > maxMassNumber_)
throw std::runtime_error("QgsjetII projectile mass outside range.");
std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
std::uniform_int_distribution select(0, 1);
qgsjet_hadron_type = nucleons[select(rng_)];
} else {
// from conex: replace pi0 or rho0 with pi+/pi- in alternating sequence
if (qgsjet_hadron_type == QgsjetIIHadronType::NeutralLightMesonType) {
qgsjet_hadron_type = alternate_;
alternate_ = (alternate_ == QgsjetIIHadronType::PiPlusType
? QgsjetIIHadronType::PiMinusType
: QgsjetIIHadronType::PiPlusType);
}
}
cout << "Interaction: projectile qgsjetII code/A: " << projectileMassNumber << " "
<< corsikaBeamId << endl;
int qgsjet_hadron_type_int = static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
cout << "Interaction: "
<< " DoInteraction: E(GeV):" << projectileEnergyLab / 1_GeV << endl;
count_++;
qgini_(projectileEnergyLab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber,
targetMassNumber);
qgconf_();
// bookkeeping
MomentumVector Plab_final(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV;
// to read the secondaries
// define rotation to and from CoM frame
// CoM frame definition in QgsjetII projectile: +z
auto const& originalCS = projectileMomentumLab.GetCoordinateSystem();
geometry::CoordinateSystem const zAxisFrame =
originalCS.RotateToZ(projectileMomentumLab);
// nuclear projectile fragments
QGSJetIIFragmentsStack qfs;
for (auto& fragm : qfs) {
particles::Code idFragm = particles::Code::Nucleus;
const int A = fragm.GetFragmentSize();
int Z = 0;
switch (A) {
case 1: { // proton/neutron
std::uniform_real_distribution<double> select;
if (select(rng_) > 0.5) {
idFragm = particles::Code::Proton;
Z = 1;
} else {
idFragm = particles::Code::Neutron;
Z = 0;
}
const HEPMassType nucleonMass = particles::GetMass(idFragm);
auto momentum = geometry::Vector(
zAxisFrame, corsika::geometry::QuantityVector<hepmomentum_d>{
0.0_GeV, 0.0_GeV,
sqrt((projectileEnergyLabPerNucleon + nucleonMass) *
(projectileEnergyLabPerNucleon - nucleonMass))});
auto const energy = sqrt(momentum.squaredNorm() + square(nucleonMass));
momentum.rebase(originalCS); // transform back into standard lab frame
std::cout << "secondary fragment> id=" << idFragm
<< " p=" << momentum.GetComponents() << std::endl;
auto pnew = vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{idFragm, energy, momentum,
pOrig, tOrig});
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
} break;
case 2: // deuterium
Z = 1;
break;
case 3: // tritium
Z = 1;
break;
case 4: // helium
Z = 2;
break;
default: // nucleus
{
Z = int(A / 2.15 + 0.7);
}
}
if (idFragm == particles::Code::Nucleus) { // thus: not p or n
const HEPMassType nucleusMass =
particles::Proton::GetMass() * Z + particles::Neutron::GetMass() * (A - Z);
auto momentum = geometry::Vector(
zAxisFrame, geometry::QuantityVector<hepmomentum_d>{
0.0_GeV, 0.0_GeV,
sqrt((projectileEnergyLabPerNucleon * A + nucleusMass) *
(projectileEnergyLabPerNucleon * A - nucleusMass))});
auto const energy = sqrt(momentum.squaredNorm() + square(nucleusMass));
momentum.rebase(originalCS); // transform back into standard lab frame
std::cout << "secondary fragment> id=" << idFragm
<< " p=" << momentum.GetComponents() << " A=" << A << " Z=" << Z
<< std::endl;
auto pnew = vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType, unsigned short, unsigned short>{
idFragm, energy, momentum, pOrig, tOrig, A, Z});
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
}
}
// secondaries
QGSJetIIStack qs;
for (auto& psec : qs) {
auto momentum = psec.GetMomentum(zAxisFrame);
auto const energy = psec.GetEnergy();
momentum.rebase(originalCS); // transform back into standard lab frame
std::cout << "secondary fragment> id="
<< process::qgsjetII::ConvertFromQgsjetII(psec.GetPID())
<< " p=" << momentum.GetComponents() << std::endl;
auto pnew = vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{
process::qgsjetII::ConvertFromQgsjetII(psec.GetPID()), energy, momentum,
pOrig, tOrig});
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
}
cout << "conservation (all GeV): Ecm_final= n/a" /* << Ecm_final / 1_GeV*/ << endl
<< "Elab_final=" << Elab_final / 1_GeV
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents()
<< ", N_wounded,targ="
<< QGSJetIIFragmentsStackData::GetWoundedNucleonsTarget()
<< ", N_wounded,proj="
<< QGSJetIIFragmentsStackData::GetWoundedNucleonsProjectile()
<< ", N_fragm,proj=" << qfs.GetSize() << endl;
}
return process::EProcessReturn::eOk;
}
} // namespace corsika::process::qgsjetII
Processes/QGSJetII/Interaction.h deleted 100644 → 0
View file @ cc68b005
/*
* (c) Copyright 2020 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.
*/
#pragma once
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/process/qgsjetII/ParticleConversion.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
#include <string>
namespace corsika::process::qgsjetII {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
std::string data_path_;
int count_ = 0;
bool initialized_ = false;
QgsjetIIHadronType alternate_ =
QgsjetIIHadronType::PiPlusType; // for pi0, rho0 projectiles
public:
Interaction(const std::string& dataPath = "");
~Interaction();
bool WasInitialized() { return initialized_; }
int GetMaxTargetMassNumber() const { return maxMassNumber_; }
bool IsValidTarget(corsika::particles::Code TargetId) const {
return (corsika::particles::GetNucleusA(TargetId) < maxMassNumber_) &&
corsika::particles::IsNucleus(TargetId);
}
corsika::units::si::CrossSectionType GetCrossSection(
const corsika::particles::Code, const corsika::particles::Code,
const corsika::units::si::HEPEnergyType, const unsigned int Abeam = 0,
const unsigned int Atarget = 0) const;
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const;
/**
In this function QGSJETII is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
*/
template <typename TProjectile>
corsika::process::EProcessReturn DoInteraction(TProjectile&);
private:
corsika::random::RNG& rng_ =
corsika::random::RNGManager::GetInstance().GetRandomStream("qgsjet");
static constexpr int maxMassNumber_ = 208;
};
} // namespace corsika::process::qgsjetII
Processes/QGSJetII/ParticleConversion.cc deleted 100644 → 0
View file @ cc68b005
/*
* (c) Copyright 2020 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/particles/ParticleProperties.h>
#include <corsika/process/qgsjetII/ParticleConversion.h>
using namespace corsika::process::qgsjetII;
Processes/QGSJetII/QGSJetIIFragmentsStack.h deleted 100644 → 0
View file @ cc68b005
/*
* (c) Copyright 2020 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.
*/
#pragma once
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/process/qgsjetII/ParticleConversion.h>
#include <corsika/process/qgsjetII/qgsjet-II-04.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process::qgsjetII {
class QGSJetIIFragmentsStackData {
public:
void Dump() const {}
void Clear() {
qgarr13_.nsf = 0;
qgarr55_.nwt = 0;
}
unsigned int GetSize() const { return qgarr13_.nsf; }
unsigned int GetCapacity() const { return iapmax; }
static unsigned int GetWoundedNucleonsTarget() { return qgarr55_.nwt; }
static unsigned int GetWoundedNucleonsProjectile() { return qgarr55_.nwp; }
int GetFragmentSize(const unsigned int i) const { return qgarr13_.iaf[i]; }
void SetFragmentSize(const unsigned int i, const int v) { qgarr13_.iaf[i] = v; }
void Copy(const unsigned int i1, const unsigned int i2) {
qgarr13_.iaf[i2] = qgarr13_.iaf[i1];
}
void Swap(const unsigned int i1, const unsigned int i2) {
std::swap(qgarr13_.iaf[i1], qgarr13_.iaf[i2]);
}
void IncrementSize() { qgarr13_.nsf++; }
void DecrementSize() {
if (qgarr13_.nsf > 0) { qgarr13_.nsf--; }
}
};
template <typename StackIteratorInterface>
class FragmentsInterface : public corsika::stack::ParticleBase<StackIteratorInterface> {
using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
public:
void SetParticleData(const int vSize) { SetFragmentSize(vSize); }
void SetParticleData(FragmentsInterface<StackIteratorInterface>& /*parent*/,
const int vSize) {
SetFragmentSize(vSize);
}
void SetFragmentSize(const int v) { GetStackData().SetFragmentSize(GetIndex(), v); }
double GetFragmentSize() const { return GetStackData().GetFragmentSize(GetIndex()); }
};
typedef corsika::stack::Stack<QGSJetIIFragmentsStackData, FragmentsInterface>
QGSJetIIFragmentsStack;
} // end namespace corsika::process::qgsjetII
Processes/QGSJetII/QGSJetIIStack.h deleted 100644 → 0
View file @ cc68b005
/*
* (c) Copyright 2020 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.
*/
#pragma once
#include <corsika/geometry/CoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/process/qgsjetII/ParticleConversion.h>
#include <corsika/process/qgsjetII/qgsjet-II-04.h>
#include <corsika/stack/Stack.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::process::qgsjetII {
typedef corsika::geometry::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
class QGSJetIIStackData {
public:
void Dump() const {}
void Clear() {
qgarr12_.nsp = 0;
qgarr13_.nsf = 0;
qgarr55_.nwt = 0;
}
unsigned int GetSize() const { return qgarr12_.nsp; }
unsigned int GetCapacity() const { return nptmax; }
void SetId(const unsigned int i, const int v) { qgarr14_.ich[i] = v; }
void SetEnergy(const unsigned int i, const corsika::units::si::HEPEnergyType v) {
using namespace corsika::units::si;
qgarr14_.esp[i][0] = v / 1_GeV;
}
void SetMomentum(const unsigned int i, const MomentumVector& v) {
using namespace corsika::units::si;
auto tmp = v.GetComponents();
qgarr14_.esp[i][2] = tmp[0] / 1_GeV;
qgarr14_.esp[i][3] = tmp[1] / 1_GeV;
qgarr14_.esp[i][1] = tmp[2] / 1_GeV;
}
int GetId(const unsigned int i) const { return qgarr14_.ich[i]; }
corsika::units::si::HEPEnergyType GetEnergy(const int i) const {
using namespace corsika::units::si;
return qgarr14_.esp[i][0] * 1_GeV;
}
MomentumVector GetMomentum(const unsigned int i,
const corsika::geometry::CoordinateSystem& CS) const {
using namespace corsika::units::si;
geometry::QuantityVector<hepmomentum_d> components = {qgarr14_.esp[i][2] * 1_GeV,
qgarr14_.esp[i][3] * 1_GeV,
qgarr14_.esp[i][1] * 1_GeV};
return MomentumVector(CS, components);
}
void Copy(const unsigned int i1, const unsigned int i2) {
qgarr14_.ich[i2] = qgarr14_.ich[i1];
for (unsigned int i = 0; i < 4; ++i) qgarr14_.esp[i2][i] = qgarr14_.esp[i1][i];
}
void Swap(const unsigned int i1, const unsigned int i2) {
std::swap(qgarr14_.ich[i1], qgarr14_.ich[i2]);
for (unsigned int i = 0; i < 4; ++i)
std::swap(qgarr14_.esp[i1][i], qgarr14_.esp[i2][i]);
}
void IncrementSize() { qgarr12_.nsp++; }
void DecrementSize() {
if (qgarr12_.nsp > 0) { qgarr12_.nsp--; }
}
};
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, const corsika::units::si::HEPEnergyType vE,
const MomentumVector& vP,
const corsika::units::si::HEPMassType) {
SetPID(vID);
SetEnergy(vE);
SetMomentum(vP);
}
void SetParticleData(ParticleInterface<StackIteratorInterface>& /*parent*/,
const int vID, const corsika::units::si::HEPEnergyType vE,
const MomentumVector& vP,
const corsika::units::si::HEPMassType) {
SetPID(vID);
SetEnergy(vE);
SetMomentum(vP);
}
void SetEnergy(const corsika::units::si::HEPEnergyType v) {
GetStackData().SetEnergy(GetIndex(), v);
}
corsika::units::si::HEPEnergyType GetEnergy() const {
return GetStackData().GetEnergy(GetIndex());
}
void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
corsika::process::qgsjetII::QgsjetIICode GetPID() const {
return static_cast<corsika::process::qgsjetII::QgsjetIICode>(
GetStackData().GetId(GetIndex()));
}
MomentumVector GetMomentum(const corsika::geometry::CoordinateSystem& CS) const {
return GetStackData().GetMomentum(GetIndex(), CS);
}
void SetMomentum(const MomentumVector& v) {
GetStackData().SetMomentum(GetIndex(), v);
}
};
typedef corsika::stack::Stack<QGSJetIIStackData, ParticleInterface> QGSJetIIStack;
} // end namespace corsika::process::qgsjetII
Processes/QGSJetII/testQGSJetII.cc deleted 100644 → 0
View file @ cc68b005
/*
* (c) Copyright 2020 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/process/qgsjetII/Interaction.h>
#include <corsika/process/qgsjetII/ParticleConversion.h>
#include <corsika/random/RNGManager.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::process::qgsjetII;
using namespace corsika::units::si;
template <typename TStackView>
auto sumCharge(TStackView const& view) {
int totalCharge = 0;
for (auto const& p : view) { totalCharge += particles::GetChargeNumber(p.GetPID()); }
return totalCharge;
}
template <typename TStackView>
auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS) {
geometry::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
for (auto const& p : view) { sum += p.GetMomentum(); }
return sum;
}
TEST_CASE("QgsjetII", "[processes]") {
SECTION("QgsjetII -> Corsika") {
CHECK(particles::PiPlus::GetCode() == process::qgsjetII::ConvertFromQgsjetII(
process::qgsjetII::QgsjetIICode::PiPlus));
}
SECTION("Corsika -> QgsjetII") {
CHECK(process::qgsjetII::ConvertToQgsjetII(particles::PiMinus::GetCode()) ==
process::qgsjetII::QgsjetIICode::PiMinus);
CHECK(process::qgsjetII::ConvertToQgsjetIIRaw(particles::Proton::GetCode()) == 2);
}
SECTION("canInteractInQgsjetII") {
CHECK(process::qgsjetII::CanInteract(particles::Proton::GetCode()));
CHECK(process::qgsjetII::CanInteract(particles::Code::KPlus));
CHECK(process::qgsjetII::CanInteract(particles::Nucleus::GetCode()));
// CHECK(process::qgsjetII::CanInteract(particles::Helium::GetCode()));
CHECK_FALSE(process::qgsjetII::CanInteract(particles::EtaC::GetCode()));
CHECK_FALSE(process::qgsjetII::CanInteract(particles::SigmaC0::GetCode()));
}
SECTION("cross-section type") {
CHECK(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Neutron) ==
process::qgsjetII::QgsjetIIXSClass::Baryons);
CHECK(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::K0Long) ==
process::qgsjetII::QgsjetIIXSClass::Kaons);
CHECK(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::Proton) ==
process::qgsjetII::QgsjetIIXSClass::Baryons);
CHECK(process::qgsjetII::GetQgsjetIIXSCode(particles::Code::PiMinus) ==
process::qgsjetII::QgsjetIIXSClass::LightMesons);
}
}
#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>
#include <corsika/process/qgsjetII/qgsjet-II-04.h>
using namespace corsika::units::si;
using namespace corsika::units;
TEST_CASE("QgsjetIIInterface", "[processes]") {
// setup environment, geometry
environment::Environment<environment::IMediumModel> env;
auto& universe = *(env.GetUniverse());
auto theMedium =
environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
auto const* nodePtr = theMedium.get();
universe.AddChild(std::move(theMedium));
const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
corsika::random::RNGManager::GetInstance().RegisterRandomStream("qgsjet");
SECTION("InteractionInterface") {
setup::Stack stack;
const HEPEnergyType E0 = 100_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, pos, 0_ns});
particle.SetNode(nodePtr);
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
auto const projectileMomentum = projectile.GetMomentum();
Interaction model;
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
CHECK(length / (1_g / square(1_cm)) == Approx(93.47).margin(0.1));
/***********************************
It as turned out already two times (#291 and #307) that the detailed output of
QGSJetII event generation depends on the gfortran version used. This is not reliable
and cannot be tested in a unit test here. One related problem was already found (#291)
and is realted to undefined behaviour in the evaluation of functions in logical
expressions. It is not clear if #307 is the same issue.
CHECK(view.GetSize() == 14);
CHECK(sumCharge(view) == 2);
************************************/
auto const secMomSum = sumMomentum(view, projectileMomentum.GetCoordinateSystem());
CHECK((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
Approx(0).margin(1e-2));
}
}
Processes/Sibyll/CMakeLists.txt deleted 100644 → 0
View file @ cc68b005
add_custom_command (
OUTPUT ${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
COMMAND ${PROJECT_SOURCE_DIR}/Processes/Sibyll/code_generator.py
${PROJECT_BINARY_DIR}/Framework/Particles/particle_db.pkl
${PROJECT_SOURCE_DIR}/Processes/Sibyll/sibyll_codes.dat
DEPENDS code_generator.py
sibyll_codes.dat
CORSIKAparticles
WORKING_DIRECTORY
${PROJECT_BINARY_DIR}/Processes/Sibyll/
COMMENT "Generate conversion tables for particle codes SIBYLL <-> CORSIKA"
VERBATIM
)
set (
MODEL_SOURCES
ParticleConversion.cc
Interaction.cc
Decay.cc
NuclearInteraction.cc
sibyll2.3d.f
nuclib.f
signuc.f
sibyll2.3d.cc
gasdev.f
)
set (
MODEL_HEADERS
ParticleConversion.h
sibyll2.3d.h
nuclib.h
SibStack.h
Decay.h
Interaction.h
NuclearInteraction.h
${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
)
set (
MODEL_NAMESPACE
corsika/process/sibyll
)
add_library (ProcessSibyll STATIC ${MODEL_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessSibyll ${MODEL_NAMESPACE} ${MODEL_HEADERS})
# ....................................................
# since Generated.inc is an automatically produced file in the build directory,
# create a symbolic link into the source tree, so that it can be found and edited more easily
# this is not needed for the build to succeed! .......
add_custom_command (
OUTPUT ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc
COMMAND ${CMAKE_COMMAND} -E create_symlink ${PROJECT_BINARY_DIR}/include/corsika/process/sibyll/Generated.inc ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc
COMMENT "Generate link in source-dir: ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc"
)
add_custom_target (SourceDirLinkSib DEPENDS ${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc)
add_dependencies (ProcessSibyll SourceDirLinkSib)
# .....................................................
set_target_properties (
ProcessSibyll
PROPERTIES
VERSION ${PROJECT_VERSION}
SOVERSION 1
)
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessSibyll
CORSIKAprocesssequence
CORSIKAparticles
CORSIKAutilities
CORSIKAunits
CORSIKAthirdparty
CORSIKAgeometry
CORSIKAenvironment
CORSIKAsetup
CORSIKArandom
)
target_include_directories (
ProcessSibyll
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
TARGETS ProcessSibyll
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
)
# --------------------
# code unit testing
CORSIKA_ADD_TEST(testSibyll
SOURCES
testSibyll.cc
${MODEL_HEADERS}
)
target_link_libraries (
testSibyll
ProcessSibyll
CORSIKAtesting
)
Processes/Sibyll/Decay.cc deleted 100644 → 0
View file @ cc68b005
/*
* (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/process/sibyll/Decay.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
using std::cout;
using std::endl;
using std::tuple;
using std::vector;
using namespace corsika;
using namespace corsika::setup;
using SetupView = corsika::setup::StackView;
using SetupProjectile = corsika::setup::StackView::ParticleType;
using SetupParticle = corsika::setup::Stack::ParticleType;
namespace corsika::process::sibyll {
Decay::Decay() {
// switch off decays to avoid internal decay chains
SetAllStable();
// handle all decays by default
handleAllDecays_ = true;
}
Decay::Decay(std::set<particles::Code> vHandled)
: handleAllDecays_(false)
, handledDecays_(vHandled) {
SetAllStable();
}
Decay::~Decay() { cout << "Sibyll::Decay n=" << fCount << endl; }
bool Decay::CanHandleDecay(const particles::Code vParticleCode) {
using namespace corsika::particles;
// if known to sibyll and not proton or neutrino it can decay
if (vParticleCode == Code::Proton || vParticleCode == Code::AntiProton ||
vParticleCode == Code::NuE || vParticleCode == Code::NuMu ||
vParticleCode == Code::NuTau || vParticleCode == Code::NuEBar ||
vParticleCode == Code::NuMuBar || vParticleCode == Code::NuTauBar ||
vParticleCode == Code::Electron || vParticleCode == Code::Positron)
return false;
else if (process::sibyll::ConvertToSibyllRaw(
vParticleCode)) // non-zero for particles known to sibyll
return true;
else
return false;
}
void Decay::SetHandleDecay(const particles::Code vParticleCode) {
handleAllDecays_ = false;
cout << "Sibyll::Decay: set to handle decay of " << vParticleCode << endl;
if (Decay::CanHandleDecay(vParticleCode))
handledDecays_.insert(vParticleCode);
else
throw std::runtime_error("this decay can not be handled by sibyll!");
}
void Decay::SetHandleDecay(const vector<particles::Code> vParticleList) {
handleAllDecays_ = false;
for (auto p : vParticleList) Decay::SetHandleDecay(p);
}
bool Decay::IsDecayHandled(const corsika::particles::Code vParticleCode) {
if (handleAllDecays_ && Decay::CanHandleDecay(vParticleCode))
return true;
else
return Decay::handledDecays_.find(vParticleCode) != Decay::handledDecays_.end()
? true
: false;
}
void Decay::SetStable(const vector<particles::Code> vParticleList) {
for (auto p : vParticleList) Decay::SetStable(p);
}
void Decay::SetUnstable(const vector<particles::Code> vParticleList) {
for (auto p : vParticleList) Decay::SetUnstable(p);
}
bool Decay::IsStable(const particles::Code vCode) {
return abs(process::sibyll::ConvertToSibyllRaw(vCode)) <= 0 ? true : false;
}
bool Decay::IsUnstable(const particles::Code vCode) {
return abs(process::sibyll::ConvertToSibyllRaw(vCode)) > 0 ? true : false;
}
void Decay::SetDecay(const particles::Code vCode, const bool vMakeUnstable) {
vMakeUnstable ? SetUnstable(vCode) : SetStable(vCode);
}
void Decay::SetUnstable(const particles::Code vCode) {
cout << "Sibyll::Decay: setting " << vCode << " unstable.." << endl;
const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
}
void Decay::SetStable(const particles::Code vCode) {
cout << "Sibyll::Decay: setting " << vCode << " stable.." << endl;
const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
}
void Decay::SetAllStable() {
for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = -1 * abs(s_csydec_.idb[i]);
}
void Decay::SetAllUnstable() {
for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = abs(s_csydec_.idb[i]);
}
void Decay::PrintDecayConfig(const particles::Code vCode) {
cout << "Decay: Sibyll decay configuration:" << endl;
const int sibCode = process::sibyll::ConvertToSibyllRaw(vCode);
const int absSibCode = abs(sibCode);
cout << vCode << " is ";
if (s_csydec_.idb[absSibCode - 1] <= 0)
cout << "stable" << endl;
else
cout << "unstable" << endl;
}
void Decay::PrintDecayConfig() {
cout << "Sibyll::Decay: decay configuration:" << endl;
if (handleAllDecays_)
cout << " all particles known to Sibyll are handled by Sibyll::Decay!" << endl;
else
for (auto& pCode : handledDecays_)
cout << "Decay of " << pCode << " is handled by Sibyll!" << endl;
}
template <>
units::si::TimeType Decay::GetLifetime(SetupParticle const& vP) {
using namespace units::si;
const particles::Code pid = vP.GetPID();
if (Decay::IsDecayHandled(pid)) {
HEPEnergyType E = vP.GetEnergy();
HEPMassType m = vP.GetMass();
const double gamma = E / m;
const TimeType t0 = particles::GetLifetime(vP.GetPID());
auto const lifetime = gamma * t0;
const auto mkin =
(E * E - vP.GetMomentum().squaredNorm()); // delta_mass(vP.GetMomentum(), E, m);
cout << "Sibyll::Decay: code: " << vP.GetPID() << endl;
cout << "Sibyll::Decay: MinStep: t0: " << t0 << endl;
cout << "Sibyll::Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
cout << "Sibyll::Decay: momentum: " << vP.GetMomentum().GetComponents() / 1_GeV
<< " GeV" << endl;
cout << "Sibyll::Decay: momentum: shell mass-kin. inv. mass "
<< mkin / 1_GeV / 1_GeV << " " << m / 1_GeV * m / 1_GeV << endl;
auto sib_id = process::sibyll::ConvertToSibyllRaw(vP.GetPID());
cout << "Sibyll::Decay: sib mass: " << get_sibyll_mass2(sib_id) << endl;
cout << "Sibyll::Decay: MinStep: gamma: " << gamma << endl;
cout << "Sibyll::Decay: MinStep: tau (s): " << lifetime / 1_s << endl;
return lifetime;
} else
return std::numeric_limits<double>::infinity() * 1_s;
}
template <>
void Decay::DoDecay(SetupProjectile& vP) {
using geometry::Point;
using namespace units::si;
const particles::Code pCode = vP.GetPID();
// check if sibyll is configured to handle this decay!
if (!IsDecayHandled(pCode))
throw std::runtime_error("STOP! Sibyll not configured to execute this decay!");
fCount++;
SibStack ss;
ss.Clear();
// copy particle to sibyll stack
ss.AddParticle(process::sibyll::ConvertToSibyllRaw(pCode), vP.GetEnergy(),
vP.GetMomentum(),
// setting particle mass with Corsika values, may be inconsistent
// with sibyll internal values
particles::GetMass(pCode));
// remember position
Point const decayPoint = vP.GetPosition();
TimeType const t0 = vP.GetTime();
// remember if particles is unstable
// auto const priorIsUnstable = IsUnstable(pCode);
// switch on decay for this particle
SetUnstable(pCode);
PrintDecayConfig(pCode);
// call sibyll decay
cout << "Decay: calling Sibyll decay routine.." << endl;
decsib_();
// print output
int print_unit = 6;
sib_list_(print_unit);
// reset to stable
SetStable(pCode);
// copy particles from sibyll stack to corsika
for (auto& psib : ss) {
// FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
if (psib.HasDecayed()) continue;
// add to corsika stack
vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
geometry::Point, units::si::TimeType>{
process::sibyll::ConvertFromSibyll(psib.GetPID()), psib.GetEnergy(),
psib.GetMomentum(), decayPoint, t0});
}
// empty sibyll stack
ss.Clear();
}
} // namespace corsika::process::sibyll
Processes/Sibyll/Decay.h deleted 100644 → 0
View file @ cc68b005
/*
* (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.
*/
#pragma once
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/DecayProcess.h>
#include <corsika/process/SecondariesProcess.h>
#include <set>
#include <vector>
namespace corsika::process {
namespace sibyll {
class Decay : public corsika::process::DecayProcess<Decay> {
int fCount = 0;
bool handleAllDecays_ = true;
public:
Decay();
Decay(std::set<particles::Code>);
~Decay();
void PrintDecayConfig(const corsika::particles::Code);
void PrintDecayConfig();
void SetHadronsUnstable();
// is Sibyll::Decay set to handle the decay of this particle?
bool IsDecayHandled(const corsika::particles::Code);
// is decay possible in principle?
bool CanHandleDecay(const corsika::particles::Code);
// set Sibyll::Decay to handle the decay of this particle!
void SetHandleDecay(const corsika::particles::Code);
// set Sibyll::Decay to handle the decay of this list of particles!
void SetHandleDecay(const std::vector<particles::Code>);
// set Sibyll::Decay to handle all particle decays
void SetHandleAllDecay();
template <typename TParticle>
corsika::units::si::TimeType GetLifetime(TParticle const&);
/**
In this function SIBYLL is called to produce to decay the input particle.
*/
template <typename TProjectile>
void DoDecay(TProjectile&);
template <typename TParticleView>
EProcessReturn DoSecondaries(TParticleView&);
private:
// internal routines to set particles stable and unstable in the COMMON blocks in
// sibyll
void SetStable(const std::vector<particles::Code>);
void SetUnstable(const std::vector<particles::Code>);
void SetStable(const corsika::particles::Code);
void SetUnstable(const corsika::particles::Code);
// internally set all particles to decay/not to decay
void SetAllUnstable();
void SetAllStable();
// will this particle be stable in sibyll ?
bool IsStable(const corsika::particles::Code);
// will this particle decay in sibyll ?
bool IsUnstable(const corsika::particles::Code);
// set particle with input code to decay or not
void SetDecay(const particles::Code, const bool);
std::set<particles::Code> handledDecays_;
};
} // namespace sibyll
} // namespace corsika::process
Processes/Sibyll/Interaction.cc deleted 100644 → 0
View file @ cc68b005
/*
* (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/process/sibyll/Interaction.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/FourVector.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/SibStack.h>
#include <corsika/process/sibyll/sibyll2.3d.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/utl/COMBoost.h>
#include <tuple>
using std::cout;
using std::endl;
using std::tuple;
using namespace corsika;
using namespace corsika::setup;
using SetupParticle = setup::Stack::StackIterator;
using SetupProjectile = setup::StackView::StackIterator;
using Track = Trajectory;
namespace corsika::process::sibyll {
bool Interaction::initialized_ = false;
Interaction::Interaction() {
using random::RNGManager;
// initialize Sibyll
if (!initialized_) {
sibyll_ini_();
initialized_ = true;
}
}
Interaction::~Interaction() {
cout << "Sibyll::Interaction n=" << count_ << " Nnuc=" << nucCount_ << endl;
}
void Interaction::SetAllStable() {
for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = -1 * abs(s_csydec_.idb[i]);
}
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
Interaction::GetCrossSection(const particles::Code BeamId,
const particles::Code TargetId,
const units::si::HEPEnergyType CoMenergy) const {
using namespace units::si;
double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3];
const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
if (!IsValidCoMEnergy(CoMenergy)) {
throw std::runtime_error(
"Interaction: GetCrossSection: CoM energy outside range for Sibyll!");
}
const double dEcm = CoMenergy / 1_GeV;
if (particles::IsNucleus(TargetId)) {
const int iTarget = particles::GetNucleusA(TargetId);
if (iTarget > maxTargetMassNumber_ || 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 == particles::Proton::GetCode()) {
sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
} 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_mb, sigEla * 1_mb);
}
template <>
units::si::GrammageType Interaction::GetInteractionLength(
SetupParticle const& vP) const {
using namespace units;
using namespace units::si;
using namespace geometry;
// coordinate system, get global frame of reference
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
const particles::Code corsikaBeamId = vP.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);
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy
HEPEnergyType Elab = vP.GetEnergy() + constants::nucleonMass;
MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
pTotLab += vP.GetMomentum();
pTotLab += pTarget;
auto const pTotLabNorm = pTotLab.norm();
// calculate cm. energy
const HEPEnergyType ECoM = sqrt(
(Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
cout << "Interaction: LambdaInt: \n"
<< " input energy: " << vP.GetEnergy() / 1_GeV << endl
<< " beam can interact:" << kInteraction << endl
<< " beam pid:" << vP.GetPID() << endl;
// TODO: move limits into variables
// FR: removed && Elab >= 8.5_GeV
if (kInteraction && IsValidCoMEnergy(ECoM)) {
// 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..
*/
auto const* currentNode = vP.GetNode();
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
si::CrossSectionType weightedProdCrossSection = mediumComposition.WeightedSum(
[=](particles::Code targetID) -> si::CrossSectionType {
return std::get<0>(this->GetCrossSection(corsikaBeamId, targetID, ECoM));
});
cout << "Interaction: "
<< "IntLength: weighted CrossSection (mb): " << weightedProdCrossSection / 1_mb
<< endl;
// calculate interaction length in medium
GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
units::constants::u / weightedProdCrossSection;
cout << "Interaction: "
<< "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
<< 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 <>
process::EProcessReturn Interaction::DoInteraction(SetupProjectile& vP) {
using namespace utl;
using namespace units;
using namespace units::si;
using namespace geometry;
const auto corsikaBeamId = vP.GetPID();
cout << "ProcessSibyll: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
<< process::sibyll::CanInteract(corsikaBeamId) << endl;
if (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)) {
// position and time of interaction, not used in Sibyll
Point const pOrig = vP.GetPosition();
TimeType const tOrig = vP.GetTime();
// define projectile
HEPEnergyType const eProjectileLab = vP.GetEnergy();
auto const pProjectileLab = vP.GetMomentum();
const CoordinateSystem& originalCS = pProjectileLab.GetCoordinateSystem();
// define target
// for Sibyll is always a single nucleon
// FOR NOW: target is always at rest
const auto eTargetLab = 0_GeV + constants::nucleonMass;
const auto pTargetLab = MomentumVector(originalCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargLab(eTargetLab, pTargetLab);
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, constants::nucleonMass);
auto const& csPrime = boost.GetRotatedCS();
// 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(csPrime) / 1_GeV << endl;
cout << "Interaction: etarget CoM: " << PtargCoM.GetTimeLikeComponent() / 1_GeV
<< endl
<< "Interaction: ptarget CoM: "
<< PtargCoM.GetSpaceLikeComponents().GetComponents(csPrime) / 1_GeV << endl;
cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "Interaction: time: " << tOrig << endl;
HEPEnergyType Etot = eProjectileLab + eTargetLab;
MomentumVector Ptot = vP.GetMomentum();
// invariant mass, i.e. cm. energy
HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
// sample target mass number
auto const* currentNode = vP.GetNode();
auto const& 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<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] = GetCrossSection(corsikaBeamId, targetId, Ecm);
[[maybe_unused]] const auto& dummy_sigEla = sigEla;
cross_section_of_components[i] = sigProd;
}
const auto targetCode =
mediumComposition.SampleTarget(cross_section_of_components, RNG_);
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 == particles::Proton::GetCode()) targetSibCode = 1;
cout << "Interaction: sibyll code: " << targetSibCode << endl;
if (targetSibCode > maxTargetMassNumber_ || 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);
cout << "Interaction: "
<< " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
<< " Ecm(GeV): " << Ecm / 1_GeV << endl;
if (Ecm > GetMaxEnergyCoM())
throw std::runtime_error("Interaction::DoInteraction: CoM energy too high!");
// FR: removed eProjectileLab < 8.5_GeV ||
if (Ecm < GetMinEnergyCoM()) {
cout << "Interaction: "
<< " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << endl;
throw std::runtime_error("energy too low for SIBYLL");
} else {
count_++;
// Sibyll does not know about units..
const double sqs = Ecm / 1_GeV;
// running sibyll, filling stack
sibyll_(kBeam, targetSibCode, sqs);
// print final state
int print_unit = 6;
sib_list_(print_unit);
nucCount_ += get_nwounded() - 1;
// add particles from sibyll to stack
// link to sibyll stack
SibStack ss;
MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
// abort on particles that have decayed in Sibyll. Should not happen!
if (psib.HasDecayed())
throw std::runtime_error("found particle that decayed in SIBYLL!");
// transform 4-momentum to lab. frame
// note that the momentum needs to be rotated back
auto const tmp = psib.GetMomentum().GetComponents();
auto const pCoM = Vector<hepmomentum_d>(csPrime, tmp);
HEPEnergyType const eCoM = psib.GetEnergy();
auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
auto const p3lab = Plab.GetSpaceLikeComponents();
assert(p3lab.GetCoordinateSystem() == originalCS); // just to be sure!
// add to corsika stack
auto pnew = vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{
process::sibyll::ConvertFromSibyll(psib.GetPID()),
Plab.GetTimeLikeComponent(), p3lab, pOrig, tOrig});
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
Ecm_final += psib.GetEnergy();
}
cout << "conservation (all GeV):" << endl
<< "Ecm_initial(per nucleon)=" << Ecm / 1_GeV << " Ecm_final(per nucleon)="
<< Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV << endl
<< "Elab_initial=" << Etot / 1_GeV << " Elab_final=" << Elab_final / 1_GeV
<< " diff (%)=" << (Elab_final / Etot / get_nwounded() - 1) * 100
<< " E in nucleons=" << constants::nucleonMass * get_nwounded() / 1_GeV
<< endl
<< "Plab_initial=" << (pProjectileLab / 1_GeV).GetComponents()
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
}
}
return process::EProcessReturn::eOk;
}
} // namespace corsika::process::sibyll
Processes/Sibyll/Interaction.h deleted 100644 → 0
View file @ cc68b005
/*
* (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.
*/
#pragma once
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
#include <tuple>
namespace corsika::process::sibyll {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
int count_ = 0;
int nucCount_ = 0;
static bool initialized_; ///! flag to assure init is done only once
public:
Interaction();
~Interaction();
void SetAllStable();
static bool WasInitialized() { return initialized_; }
bool IsValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) const {
return (minEnergyCoM_ <= ecm) && (ecm <= maxEnergyCoM_);
}
int GetMaxTargetMassNumber() const { return maxTargetMassNumber_; }
corsika::units::si::HEPEnergyType GetMinEnergyCoM() const { return minEnergyCoM_; }
corsika::units::si::HEPEnergyType GetMaxEnergyCoM() const { return maxEnergyCoM_; }
bool IsValidTarget(corsika::particles::Code TargetId) const {
return (corsika::particles::GetNucleusA(TargetId) < maxTargetMassNumber_) &&
corsika::particles::IsNucleus(TargetId);
}
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(const corsika::particles::Code, const corsika::particles::Code,
const corsika::units::si::HEPEnergyType) const;
template <typename TParticle>
corsika::units::si::GrammageType GetInteractionLength(TParticle const&) const;
/**
In this function SIBYLL is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
*/
template <typename TProjectile>
corsika::process::EProcessReturn DoInteraction(TProjectile&);
private:
corsika::random::RNG& RNG_ =
corsika::random::RNGManager::GetInstance().GetRandomStream("sibyll");
const corsika::units::si::HEPEnergyType minEnergyCoM_ =
10. * 1e9 * corsika::units::si::electronvolt;
const corsika::units::si::HEPEnergyType maxEnergyCoM_ =
1.e6 * 1e9 * corsika::units::si::electronvolt;
const int maxTargetMassNumber_ = 18;
};
} // namespace corsika::process::sibyll
Processes/Sibyll/NuclearInteraction.cc deleted 100644 → 0
View file @ cc68b005
/*
* (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/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/FourVector.h>
#include <corsika/process/sibyll/nuclib.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/COMBoost.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <set>
using std::cout;
using std::endl;
using std::tuple;
using std::vector;
using namespace corsika;
using namespace corsika::setup;
using Particle = Stack::ParticleType; // StackIterator; // ParticleType;
using Projectile = StackView::StackIterator; // StackView::ParticleType;
using Track = Trajectory;
namespace corsika::process::sibyll {
template <>
NuclearInteraction<SetupEnvironment>::~NuclearInteraction() {
cout << "Nuclib::NuclearInteraction n=" << count_ << " Nnuc=" << nucCount_ << endl;
}
template <>
void NuclearInteraction<SetupEnvironment>::PrintCrossSectionTable(
corsika::particles::Code pCode) {
using namespace corsika::particles;
const int k = targetComponentsIndex_.at(pCode);
Code pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
Code::Neon, Code::Argon, Code::Iron};
cout << "en/A ";
for (auto& j : pNuclei) cout << std::setw(9) << j;
cout << endl;
// loop over energy bins
for (int i = 0; i < GetNEnergyBins(); ++i) {
cout << " " << i << " ";
for (auto& n : pNuclei) {
auto const j = GetNucleusA(n);
cout << " " << std::setprecision(5) << std::setw(8)
<< cnucsignuc_.sigma[j - 1][k][i];
}
cout << endl;
}
}
template <>
void NuclearInteraction<SetupEnvironment>::InitializeNuclearCrossSections() {
using namespace corsika::particles;
using namespace units::si;
auto& universe = *(environment_.GetUniverse());
auto const allElementsInUniverse = std::invoke([&]() {
std::set<particles::Code> allElementsInUniverse;
auto collectElements = [&](auto& vtn) {
if (vtn.HasModelProperties()) {
auto const& comp =
vtn.GetModelProperties().GetNuclearComposition().GetComponents();
for (auto const c : comp) allElementsInUniverse.insert(c);
}
};
universe.walk(collectElements);
return allElementsInUniverse;
});
cout << "NuclearInteraction: initializing nuclear cross sections..." << endl;
// loop over target components, at most 4!!
int k = -1;
for (auto& ptarg : allElementsInUniverse) {
++k;
cout << "NuclearInteraction: init target component: " << ptarg << endl;
const int ib = GetNucleusA(ptarg);
if (!hadronicInteraction_.IsValidTarget(ptarg)) {
cout << "NuclearInteraction::InitializeNuclearCrossSections: target nucleus? id="
<< ptarg << endl;
throw std::runtime_error(
" target can not be handled by hadronic interaction model! ");
}
targetComponentsIndex_.insert(std::pair<Code, int>(ptarg, k));
// loop over energies, fNEnBins log. energy bins
for (int i = 0; i < GetNEnergyBins(); ++i) {
// hard coded energy grid, has to be aligned to definition in signuc2!!, no
// comment..
const units::si::HEPEnergyType Ecm = pow(10., 1. + 1. * i) * 1_GeV;
// get p-p cross sections
auto const protonId = Code::Proton;
auto const [siginel, sigela] =
hadronicInteraction_.GetCrossSection(protonId, protonId, Ecm);
const double dsig = siginel / 1_mb;
const double dsigela = sigela / 1_mb;
// loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
for (int j = 1; j < gMaxNucleusAProjectile_; ++j) {
const int jj = j + 1;
double sig_out, dsig_out, sigqe_out, dsigqe_out;
sigma_mc_(jj, ib, dsig, dsigela, gNSample_, sig_out, dsig_out, sigqe_out,
dsigqe_out);
// write to table
cnucsignuc_.sigma[j][k][i] = sig_out;
cnucsignuc_.sigqe[j][k][i] = sigqe_out;
}
}
}
cout << "NuclearInteraction: cross sections for " << targetComponentsIndex_.size()
<< " components initialized!" << endl;
for (auto& ptarg : allElementsInUniverse) {
cout << "cross section table: " << ptarg << endl;
PrintCrossSectionTable(ptarg);
}
}
template <>
units::si::CrossSectionType NuclearInteraction<SetupEnvironment>::ReadCrossSectionTable(
const int ia, particles::Code pTarget, units::si::HEPEnergyType elabnuc) {
using namespace corsika::particles;
using namespace units::si;
const int ib = targetComponentsIndex_.at(pTarget) + 1; // table index in fortran
auto const ECoMNuc = sqrt(2. * corsika::units::constants::nucleonMass * elabnuc);
if (ECoMNuc < GetMinEnergyPerNucleonCoM() || ECoMNuc > GetMaxEnergyPerNucleonCoM())
throw std::runtime_error("NuclearInteraction: energy outside tabulated range!");
const double e0 = elabnuc / 1_GeV;
double sig;
cout << "ReadCrossSectionTable: " << ia << " " << ib << " " << e0 << endl;
signuc2_(ia, ib, e0, sig);
cout << "ReadCrossSectionTable: sig=" << sig << endl;
return sig * 1_mb;
}
// TODO: remove elastic cross section?
template <>
template <>
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
NuclearInteraction<SetupEnvironment>::GetCrossSection(Particle const& vP,
const particles::Code TargetId) {
using namespace units::si;
if (vP.GetPID() != particles::Code::Nucleus)
throw std::runtime_error(
"NuclearInteraction: GetCrossSection: particle not a nucleus!");
auto const iBeamA = vP.GetNuclearA();
HEPEnergyType LabEnergyPerNuc = vP.GetEnergy() / iBeamA;
cout << "NuclearInteraction: GetCrossSection: called with: beamNuclA= " << iBeamA
<< " TargetId= " << TargetId << " LabEnergyPerNuc= " << LabEnergyPerNuc / 1_GeV
<< 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 (iBeamA > GetMaxNucleusAProjectile() || iBeamA < 2) {
cout << "NuclearInteraction: beam nucleus outside allowed range for NUCLIB!" << endl
<< "A=" << iBeamA << endl;
throw std::runtime_error(
"NuclearInteraction: GetCrossSection: beam nucleus outside allowed range for "
"NUCLIB!");
}
if (hadronicInteraction_.IsValidTarget(TargetId)) {
auto const sigProd = ReadCrossSectionTable(iBeamA, TargetId, LabEnergyPerNuc);
cout << "cross section (mb): " << sigProd / 1_mb << endl;
return std::make_tuple(sigProd, 0_mb);
} else {
throw std::runtime_error("target outside range.");
}
return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mb,
std::numeric_limits<double>::infinity() * 1_mb);
}
template <>
template <>
units::si::GrammageType NuclearInteraction<SetupEnvironment>::GetInteractionLength(
Particle const& vP) {
using namespace units;
using namespace units::si;
using namespace geometry;
// coordinate system, get global frame of reference
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
const particles::Code corsikaBeamId = vP.GetPID();
if (corsikaBeamId != particles::Code::Nucleus) {
// check if target-style nucleus (enum), these are not allowed as projectile
if (particles::IsNucleus(corsikaBeamId))
throw std::runtime_error(
"NuclearInteraction: GetInteractionLength: Wrong nucleus type. Nuclear "
"projectiles should use NuclearStackExtension!");
else {
// no nuclear interaction
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
}
// read from cross section code table
// FOR NOW: assume target is at rest
corsika::stack::MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
// total momentum and energy
HEPEnergyType Elab = vP.GetEnergy() + constants::nucleonMass;
int const nuclA = vP.GetNuclearA();
auto const ElabNuc = vP.GetEnergy() / nuclA;
corsika::stack::MomentumVector pTotLab(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
pTotLab += vP.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 * constants::nucleonMass);
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 >= gMinEnergyPerNucleonCoM_ &&
ECoMNN < gMaxEnergyPerNucleonCoM_) {
// 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..
*/
auto const* const currentNode = vP.GetNode();
auto const& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// determine average interaction length
// weighted sum
int i = -1;
si::CrossSectionType weightedProdCrossSection = 0_mb;
// 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(vP, targetId);
[[maybe_unused]] auto& dummy_elaCrossSection = elaCrossSection;
cout << "NuclearInteraction: "
<< "IntLength: nuclib return (mb): " << productionCrossSection / 1_mb
<< endl;
weightedProdCrossSection += w[i] * productionCrossSection;
}
cout << "NuclearInteraction: "
<< "IntLength: weighted CrossSection (mb): " << weightedProdCrossSection / 1_mb
<< endl;
// calculate interaction length in medium
GrammageType const int_length = mediumComposition.GetAverageMassNumber() *
units::constants::u / weightedProdCrossSection;
cout << "NuclearInteraction: "
<< "interaction length (g/cm2): " << int_length * (1_cm * 1_cm / (0.001_kg))
<< endl;
return int_length;
} else {
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
}
template <>
template <>
process::EProcessReturn NuclearInteraction<SetupEnvironment>::DoInteraction(
Projectile& vP) {
// this routine superimposes different nucleon-nucleon interactions
// in a nucleus-nucleus interaction, based the SIBYLL routine SIBNUC
using namespace units;
using namespace utl;
using namespace units::si;
using namespace geometry;
const auto ProjId = vP.GetPID();
// TODO: calculate projectile mass in nuclearStackExtension
// const auto ProjMass = vP.GetMass();
cout << "NuclearInteraction: DoInteraction: called with:" << ProjId << endl;
// check if target-style nucleus (enum)
if (ProjId != particles::Code::Nucleus)
throw std::runtime_error(
"NuclearInteraction: DoInteraction: Wrong nucleus type. Nuclear projectiles "
"should use NuclearStackExtension!");
auto const ProjMass =
vP.GetNuclearZ() * particles::Proton::GetMass() +
(vP.GetNuclearA() - vP.GetNuclearZ()) * particles::Neutron::GetMass();
cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << endl;
count_++;
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in NUCLIB
Point pOrig = vP.GetPosition();
TimeType tOrig = vP.GetTime();
cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
cout << "Interaction: time: " << tOrig << endl;
// projectile nucleon number
const int kAProj = vP.GetNuclearA();
if (kAProj > GetMaxNucleusAProjectile())
throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
// kinematics
// define projectile nucleus
HEPEnergyType const eProjectileLab = vP.GetEnergy();
auto const pProjectileLab = vP.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 = vP.GetEnergy() / kAProj;
auto const pProjectileNucLab = vP.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
// target is always at rest
const auto eTargetNucLab = 0_GeV + constants::nucleonMass;
const auto pTargetNucLab =
corsika::stack::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 (!hadronicInteraction_.IsValidCoMEnergy(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, constants::nucleonMass);
// 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 = particles::Proton::GetCode();
auto const* const currentNode = vP.GetNode();
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();
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] =
hadronicInteraction_.GetCrossSection(beamId, targetId, EcmNN);
cross_section_of_components[i] = sigProd;
[[maybe_unused]] auto sigElaCopy = sigEla; // ONLY TO AVOID COMPILER WARNINGS
}
const auto targetCode =
mediumComposition.SampleTarget(cross_section_of_components, RNG_);
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 == particles::Proton::GetCode()) kATarget = 1;
cout << "NuclearInteraction: nuclib target code: " << kATarget << endl;
if (!hadronicInteraction_.IsValidTarget(targetCode))
throw std::runtime_error("target outside range. ");
// 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 = particles::Proton::GetCode();
const auto [prodCrossSection, elaCrossSection] =
hadronicInteraction_.GetCrossSection(protonId, protonId, EcmNN);
const double sigProd = prodCrossSection / 1_mb;
const double sigEla = elaCrossSection / 1_mb;
// 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 > GetMaxNFragments())
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) {
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 = particles::Code::Proton;
else
specCode = particles::Code::Nucleus;
// TODO: mass of nuclei?
const HEPMassType mass =
particles::Proton::GetMass() * nuclZ +
(nuclA - nuclZ) * 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
vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{
specCode, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(),
pOrig, tOrig});
else
// add nucleus
vP.AddSecondary(tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
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 = particles::Code::Proton;
// CORSIKA 7 way
// elastic nucleons inherit momentum from original projectile
// neglecting momentum transfer in interaction
const double mass_ratio = particles::GetMass(elaNucCode) / ProjMass;
auto const Plab = PprojLab * mass_ratio;
vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
geometry::Point, units::si::TimeType>{
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 = particles::Proton::GetCode();
// temporarily add to stack, will be removed after interaction in DoInteraction
cout << "inelastic interaction no. " << j << endl;
auto inelasticNucleon = vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
geometry::Point, units::si::TimeType>{
pCode, PprojNucLab.GetTimeLikeComponent(),
PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig});
// create inelastic interaction
cout << "calling HadronicInteraction..." << endl;
hadronicInteraction_.DoInteraction(inelasticNucleon);
}
cout << "NuclearInteraction: DoInteraction: done" << endl;
return process::EProcessReturn::eOk;
}
template <>
NuclearInteraction<SetupEnvironment>::NuclearInteraction(
process::sibyll::Interaction& hadint, SetupEnvironment const& env)
: environment_(env)
, hadronicInteraction_(hadint) {
// initialize hadronic interaction module
// TODO: safe to run multiple initializations?
// check compatibility of energy ranges, someone could try to use low-energy model..
if (!hadronicInteraction_.IsValidCoMEnergy(GetMinEnergyPerNucleonCoM()) ||
!hadronicInteraction_.IsValidCoMEnergy(GetMaxEnergyPerNucleonCoM()))
throw std::runtime_error(
"NuclearInteraction: hadronic interaction model incompatible!");
// initialize nuclib
// TODO: make sure this does not overlap with sibyll
nuc_nuc_ini_();
// initialize cross sections
InitializeNuclearCrossSections();
}
} // namespace corsika::process::sibyll
Processes/Sibyll/NuclearInteraction.h deleted 100644 → 0
View file @ cc68b005
/*
* (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.
*/
#pragma once
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/InteractionProcess.h>
#include <corsika/random/RNGManager.h>
namespace corsika::process::sibyll {
class Interaction; // fwd-decl
/**
*
*
**/
template <class TEnvironment>
class NuclearInteraction
: public corsika::process::InteractionProcess<NuclearInteraction<TEnvironment>> {
int count_ = 0;
int nucCount_ = 0;
public:
NuclearInteraction(corsika::process::sibyll::Interaction&, TEnvironment const&);
~NuclearInteraction();
void InitializeNuclearCrossSections();
void PrintCrossSectionTable(corsika::particles::Code);
corsika::units::si::CrossSectionType ReadCrossSectionTable(
const int, corsika::particles::Code, corsika::units::si::HEPEnergyType);
corsika::units::si::HEPEnergyType GetMinEnergyPerNucleonCoM() {
return gMinEnergyPerNucleonCoM_;
}
corsika::units::si::HEPEnergyType GetMaxEnergyPerNucleonCoM() {
return gMaxEnergyPerNucleonCoM_;
}
int constexpr GetMaxNucleusAProjectile() { return gMaxNucleusAProjectile_; }
int constexpr GetMaxNFragments() { return gMaxNFragments_; }
int constexpr GetNEnergyBins() { return gNEnBins_; }
template <typename Particle>
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(Particle const& p, const corsika::particles::Code TargetId);
template <typename Particle>
corsika::units::si::GrammageType GetInteractionLength(Particle const&);
template <typename Projectile>
corsika::process::EProcessReturn DoInteraction(Projectile&);
private:
TEnvironment const& environment_;
corsika::process::sibyll::Interaction& hadronicInteraction_;
std::map<corsika::particles::Code, int> targetComponentsIndex_;
corsika::random::RNG& RNG_ =
corsika::random::RNGManager::GetInstance().GetRandomStream("sibyll");
static constexpr int gNSample_ =
500; // number of samples in MC estimation of cross section
static constexpr int gMaxNucleusAProjectile_ = 56;
static constexpr int gNEnBins_ = 6;
static constexpr int gMaxNFragments_ = 60;
// energy limits defined by table used for cross section in signuc.f
// 10**1 GeV to 10**6 GeV
static constexpr corsika::units::si::HEPEnergyType gMinEnergyPerNucleonCoM_ =
10. * 1e9 * corsika::units::si::electronvolt;
static constexpr corsika::units::si::HEPEnergyType gMaxEnergyPerNucleonCoM_ =
1.e6 * 1e9 * corsika::units::si::electronvolt;
};
} // namespace corsika::process::sibyll
Processes/Sibyll/ParticleConversion.cc deleted 100644 → 0
View file @ cc68b005
/*
* (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/particles/ParticleProperties.h>
#include <corsika/process/sibyll/ParticleConversion.h>
using namespace corsika::process::sibyll;
corsika::units::si::HEPMassType corsika::process::sibyll::GetSibyllMass(
corsika::particles::Code const pCode) {
using namespace corsika::units;
using namespace corsika::units::si;
if (pCode == corsika::particles::Code::Nucleus)
throw std::runtime_error("Cannot GetMass() of particle::Nucleus -> unspecified");
auto sCode = ConvertToSibyllRaw(pCode);
if (sCode == 0)
throw std::runtime_error("GetSibyllMass: unknown particle!");
else
return sqrt(get_sibyll_mass2(sCode)) * 1_GeV;
}
Processes/Sibyll/SibStack.h deleted 100644 → 0
View file @ cc68b005
/*
* (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.
*/
#pragma once
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/process/sibyll/sibyll2.3d.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 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
Processes/Sibyll/rndm_dbl.f deleted 100644 → 0
View file @ cc68b005
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.3d.cc deleted 100644 → 0
View file @ cc68b005
/*
* (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/process/sibyll/sibyll2.3d.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[abs(id) - 1]; }
double s_rndm_(int&) {
static corsika::random::RNG& rng =
corsika::random::RNGManager::GetInstance().GetRandomStream("sibyll");
std::uniform_real_distribution<double> dist;
return dist(rng);
}
Processes/Sibyll/sibyll_codes.dat deleted 100644 → 0
View file @ cc68b005
# 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"
# The unknown particle is to handle all particles that are not known to SIBYLL.
# It is important that sibyll-id does not overlap with any existing sibyll particle!
# Be careful
Unknown 0 0 CannotInteract
# Here is the list of particles known to sibyll
Electron 3 0 CannotInteract
Positron 2 0 CannotInteract
NuE 15 0 CannotInteract
NuEBar 16 0 CannotInteract
MuMinus 5 0 CannotInteract
MuPlus 4 0 CannotInteract
NuMu 17 0 CannotInteract
NuMuBar 18 0 CannotInteract
TauMinus 91 0 CannotInteract
TauPlus 90 0 CannotInteract
NuTau 92 0 CannotInteract
NuTauBar 93 0 CannotInteract
Gamma 1 0 CannotInteract
Pi0 6 1 Pion
# rho0 could interact but sibyll has no cross section/interaction length. was used for gamma had int
Rho0 27 0 CannotInteract
K0Long 11 1 Kaon
K0 21 0 Kaon
K0Bar 22 0 Kaon
PiPlus 7 1 Pion
PiMinus 8 1 Pion
RhoPlus 25 0 CannotInteract
RhoMinus 26 0 CannotInteract
Eta 23 0 CannotInteract
EtaPrime 24 0 CannotInteract
Pi1300Plus 62 0 CannotInteract
Pi1300Minus 63 0 CannotInteract
Pi1300_0 61 0 CannotInteract
Omega 32 0 CannotInteract
K0Short 12 1 Kaon
KStar0 30 0 CannotInteract
KStar0Bar 31 0 CannotInteract
KPlus 9 1 Kaon
KMinus 10 1 Kaon
KStarPlus 28 0 CannotInteract
KStarMinus 29 0 CannotInteract
KStar0_1430_0 66 0 CannotInteract
KStar0_1430_0Bar 67 0 CannotInteract
KStar0_1430_Plus 64 0 CannotInteract
KStar0_1430_MinusBar 65 0 CannotInteract
DPlus 59 1 Kaon
DMinus 60 1 Kaon
DStarPlus 78 0 CannotInteract
DStarMinus 79 0 CannotInteract
D0 71 1 Kaon
D0Bar 72 1 Kaon
DStar0 80 0 CannotInteract
DStar0Bar 81 0 CannotInteract
DsPlus 74 1 Kaon
DsMinus 75 1 Kaon
DStarSPlus 76 0 CannotInteract
DStarSMinus 77 0 CannotInteract
EtaC 73 0 CannotInteract
Neutron 14 1 Baryon
AntiNeutron -14 1 Baryon
Delta0 42 0 CannotInteract
Delta0Bar -42 0 CannotInteract
DeltaMinus 43 0 CannotInteract
DeltaPlusBar -43 0 CannotInteract
Proton 13 1 Baryon
AntiProton -13 1 Baryon
N1440Plus 51 0 CannotInteract
N1440MinusBar -51 0 CannotInteract
N1440_0 52 0 CannotInteract
N1440_0Bar -52 0 CannotInteract
N1710Plus 53 0 CannotInteract
N1710MinusBar -53 0 CannotInteract
N1710_0 54 0 CannotInteract
N1710_0Bar -54 0 CannotInteract
DeltaPlus 41 0 CannotInteract
DeltaMinusBar -41 0 CannotInteract
DeltaPlusPlus 40 0 CannotInteract
DeltaMinusMinusBar -40 0 CannotInteract
SigmaMinus 36 1 Baryon
SigmaPlusBar -36 1 Baryon
SigmaStarMinus 46 0 CannotInteract
SigmaStarPlusBar -46 0 CannotInteract
SigmaStarPlus 44 0 CannotInteract
SigmaStarMinusBar -44 0 CannotInteract
SigmaStar0 45 0 CannotInteract
SigmaStar0Bar -45 0 CannotInteract
Lambda0 39 1 Baryon
Lambda0Bar -39 1 Baryon
Sigma0 35 1 Baryon
Sigma0Bar -35 1 Baryon
SigmaPlus 34 1 Baryon
SigmaMinusBar -34 1 Baryon
XiMinus 38 1 Baryon
XiPlusBar -38 1 Baryon
Xi0 37 1 Baryon
Xi0Bar -37 1 Baryon
XiStarMinus 48 0 CannotInteract
XiStarPlusBar -48 0 CannotInteract
XiStar0 47 0 CannotInteract
XiStar0Bar -47 0 CannotInteract
OmegaMinus 49 0 CannotInteract
OmegaPlusBar -49 0 CannotInteract
SigmaC0 86 0 CannotInteract
SigmaC0Bar -86 0 CannotInteract
SigmaStarC0 96 0 CannotInteract
SigmaStarC0Bar -96 0 CannotInteract
LambdaCPlus 89 1 Baryon
LambdaCMinusBar -89 1 Baryon
XiC0 88 1 Baryon
XiC0Bar -88 1 Baryon
SigmaCPlus 85 0 CannotInteract
SigmaCMinusBar -85 0 CannotInteract
SigmaStarCPlus 95 0 CannotInteract
SigmaStarCMinusBar -95 0 CannotInteract
SigmaCPlusPlus 84 0 CannotInteract
SigmaCMinusMinusBar -84 0 CannotInteract
SigmaStarCPlusPlus 94 0 CannotInteract
SigmaStarCMinusMinusBar -94 0 CannotInteract
XiCPlus 87 1 Baryon
XiCMinusBar -87 1 Baryon
XiStarCPlus 97 0 CannotInteract
XiStarCMinusBar -97 0 CannotInteract
XiStarC0 98 0 CannotInteract
XiStarC0Bar -98 0 CannotInteract
OmegaC0 99 0 CannotInteract
OmegaC0Bar -99 0 CannotInteract
Jpsi 83 0 CannotInteract
Phi 33 0 CannotInteract
Processes/Sibyll/testSibyll.cc deleted 100644 → 0
View file @ cc68b005
/*
* (c) Copyright 2019 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/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>
#include <catch2/catch.hpp>
using namespace corsika;
using namespace corsika::process::sibyll;
using namespace corsika::units;
using namespace corsika::units::si;
TEST_CASE("Sibyll", "[processes]") {
SECTION("Sibyll -> Corsika") {
CHECK(particles::Electron::GetCode() ==
process::sibyll::ConvertFromSibyll(process::sibyll::SibyllCode::Electron));
}
SECTION("Corsika -> Sibyll") {
CHECK(process::sibyll::ConvertToSibyll(particles::Electron::GetCode()) ==
process::sibyll::SibyllCode::Electron);
CHECK(process::sibyll::ConvertToSibyllRaw(particles::Proton::GetCode()) == 13);
CHECK(process::sibyll::ConvertToSibyll(particles::XiStarC0::GetCode()) ==
process::sibyll::SibyllCode::XiStarC0);
}
SECTION("canInteractInSibyll") {
CHECK(process::sibyll::CanInteract(particles::Proton::GetCode()));
CHECK(process::sibyll::CanInteract(particles::Code::XiCPlus));
CHECK_FALSE(process::sibyll::CanInteract(particles::Electron::GetCode()));
CHECK_FALSE(process::sibyll::CanInteract(particles::SigmaC0::GetCode()));
CHECK_FALSE(process::sibyll::CanInteract(particles::Nucleus::GetCode()));
CHECK_FALSE(process::sibyll::CanInteract(particles::Helium::GetCode()));
}
SECTION("cross-section type") {
CHECK(process::sibyll::GetSibyllXSCode(particles::Code::Electron) == 0);
CHECK(process::sibyll::GetSibyllXSCode(particles::Code::K0Long) == 3);
CHECK(process::sibyll::GetSibyllXSCode(particles::Code::SigmaPlus) == 1);
CHECK(process::sibyll::GetSibyllXSCode(particles::Code::PiMinus) == 2);
}
SECTION("sibyll mass") {
CHECK_FALSE(process::sibyll::GetSibyllMass(particles::Code::Electron) == 0_GeV);
}
}
#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>
#include <corsika/process/sibyll/sibyll2.3d.h>
using namespace corsika::units::si;
using namespace corsika::units;
template <typename TStackView>
auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS) {
geometry::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
for (auto const& p : view) { sum += p.GetMomentum(); }
return sum;
}
TEST_CASE("SibyllInterface", "[processes]") {
// setup environment, geometry
environment::Environment<environment::IMediumModel> env;
auto& universe = *(env.GetUniverse());
auto theMedium =
environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
theMedium->SetModelProperties<MyHomogeneousModel>(
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
auto const* nodePtr = theMedium.get();
universe.AddChild(std::move(theMedium));
const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
SECTION("InteractionInterface - low energy") {
setup::Stack stack;
const HEPEnergyType E0 = 60_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = corsika::stack::MomentumVector(cs, {P0, 0_eV, 0_eV});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, pos, 0_ns});
particle.SetNode(nodePtr);
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
Interaction model;
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
auto const pSum = sumMomentum(view, cs);
/*
Interactions between hadrons (h) and nuclei (A) in Sibyll are treated in the
hadron-nucleon center-of-mass frame (hnCoM). The incoming hadron (h) and
nucleon (N) are assumed massless, such that the energy and momentum in the hnCoM are
: E_i_cm = 0.5 * SQS and P_i_cm = +- 0.5 * SQS where i is either the projectile
hadron or the target nucleon and SQS is the hadron-nucleon center-of-mass energy.
The true energies and momenta, accounting for the hadron masses, are: E_i = ( S +
m_i**2 - m_j**2 ) / (2 * SQS) and Pcm = +-
sqrt( (S-(m_j+m_i)**2) * (s-(m_j-m_i)**2) ) / (2*SQS) where m_i is the projectiles
mass and m_j is the target particles mass. In terms of lab. frame variables Pcm =
m_j * Plab_i / SQS, where Plab_i is the momentum of the projectile (i) in the lab.
and m_j is the mass of the target, i.e. the particle at rest (usually a nucleon).
Any hadron-nucleus event can contain several nucleon interactions. In case of Nw
(number of wounded nucleons) nucleons interacting in the hadron-nucleus interaction,
the total energy and momentum in the hadron(i)-nucleon(N) center-of-mass frame are:
momentum: p_projectile + p_nucleon_1 + p_nucleon_2 + .... p_nucleon_Nw = -(Nw-1) *
Pcm with center-of-mass momentum Pcm = p_projectile = - p_nucleon_i. For the energy:
E_projectile + E_nucleon_1 + ... E_nucleon_Nw = E_projectile + Nw * E_nucleon.
Using the above definitions of center-of-mass energies and momenta this leads to the
total energy: E_tot = SQS/2 * (1+Nw) + (m_N**2-m_i**2)/(2*SQS) * (Nw-1) and P_tot
= -m_N * Plab_i / SQS * (Nw-1).
A Lorentztransformation of these quantities to the lab. frame recovers Plab_i for
the total momentum, so momentum is exactly conserved, and Elab_i + Nw * m_N for the
total energy. Not surprisingly the total energy differs from the total energy before
the collision by the mass of the additional nucleons (Nw-1)*m_N. In relative terms
the additional energy is entirely negligible and as it is not kinetic energy there
is zero influence on the shower development.
Due to the ommission of the hadron masses in Sibyll, the total energy and momentum
in the center-of-mass system after the collision are just: E_tot = SQS/2 * (1+Nw)
and P_tot = SQS/2 * (1-Nw). After the Lorentztransformation the total momentum in
the lab. thus differs from the initial value by (1-Nw)/2 * ( m_N + m_i**2 / (2 *
Plab_i) ) and momentum is NOT conserved. Note however that the second term quickly
vanishes as the lab. momentum of the projectile increases. The first term is fixed
as it depends only on the number of additional nucleons, in relative terms it is
always small at high energies.
For this reason the numerical precision in these tests is limited to 5% to still
pass at low energies and no absolute check is implemented, e.g.
CHECK(pSum.GetComponents(cs).GetX() / P0 == Approx(1).margin(0.05));
CHECK((pSum - plab).norm()/1_GeV == Approx(0).margin(plab.norm() * 0.05/1_GeV));
/FR'2020
See also:
Issue 272 / MR 204
https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/merge_requests/204
*/
CHECK(pSum.GetComponents(cs).GetX() / P0 == Approx(1).margin(0.05));
CHECK(pSum.GetComponents(cs).GetY() / 1_GeV == Approx(0).margin(1e-4));
CHECK(pSum.GetComponents(cs).GetZ() / 1_GeV == Approx(0).margin(1e-4));
CHECK((pSum - plab).norm() / 1_GeV == Approx(0).margin(plab.norm() * 0.05 / 1_GeV));
CHECK(pSum.norm() / P0 == Approx(1).margin(0.05));
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
SECTION("InteractionInterface - high energy") {
setup::Stack stack;
const HEPEnergyType E0 = 60_EeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = corsika::stack::MomentumVector(cs, {P0, 0_eV, 0_eV});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Proton, E0, plab, pos, 0_ns});
particle.SetNode(nodePtr);
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
Interaction model;
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
auto const pSum = sumMomentum(view, cs);
CHECK(pSum.GetComponents(cs).GetX() / P0 == Approx(1).margin(0.001));
CHECK(pSum.GetComponents(cs).GetY() / 1_GeV == Approx(0).margin(1e-4));
CHECK(pSum.GetComponents(cs).GetZ() / 1_GeV == Approx(0).margin(1e-4));
CHECK((pSum - plab).norm() / 1_GeV == Approx(0).margin(plab.norm() * 0.001 / 1_GeV));
CHECK(pSum.norm() / P0 == Approx(1).margin(0.05));
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
SECTION("NuclearInteractionInterface") {
setup::Stack stack;
const HEPEnergyType E0 = 400_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle =
stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, E0, plab, pos, 0_ns, 4, 2});
particle.SetNode(nodePtr);
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
Interaction hmodel;
NuclearInteraction model(hmodel, env);
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
SECTION("DecayInterface") {
setup::Stack stack;
const HEPEnergyType E0 = 10_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::Proton::GetMass() * particles::Proton::GetMass());
auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Lambda0, E0, plab, pos, 0_ns});
corsika::stack::SecondaryView view(particle);
auto projectile = view.GetProjectile();
Decay model;
model.PrintDecayConfig();
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
/*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(projectile);
// run checks
// lambda decays into proton and pi- or neutron and pi+
CHECK(stack.GetSize() == 3);
}
SECTION("DecayConfiguration") {
Decay model({particles::Code::PiPlus, particles::Code::PiMinus});
CHECK(model.IsDecayHandled(particles::Code::PiPlus));
CHECK(model.IsDecayHandled(particles::Code::PiMinus));
CHECK_FALSE(model.IsDecayHandled(particles::Code::KPlus));
const std::vector<particles::Code> particleTestList = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
particles::Code::Lambda0Bar, particles::Code::D0Bar};
// setup decays
model.SetHandleDecay(particleTestList);
for (auto& pCode : particleTestList) CHECK(model.IsDecayHandled(pCode));
// individually
model.SetHandleDecay(particles::Code::KMinus);
// possible decays
CHECK_FALSE(model.CanHandleDecay(particles::Code::Proton));
CHECK_FALSE(model.CanHandleDecay(particles::Code::Electron));
CHECK(model.CanHandleDecay(particles::Code::PiPlus));
CHECK(model.CanHandleDecay(particles::Code::MuPlus));
}
}

data privacy - impressum