diff --git a/CMakeLists.txt b/CMakeLists.txt
index c06fdca0aba3e9c600b3a14919ccaff83f21c148..d29d6566625b54239eda5b76578b97498cc97c5d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -83,7 +83,8 @@ find_package( PhysUnits REQUIRED )
#+++++++++++++++++++++++++++++
# get Pythia
#
-find_package( Pythia8 REQUIRED )
+#find_package( Pythia8 REQUIRED )
+add_subdirectory(dependencies/pythia)
if(Pythia_FOUND)
include_directories(${Pythia_INCLUDE_DIR})
endif(Pythia_FOUND)
diff --git a/corsika/detail/framework/random/RNGManager.inl b/corsika/detail/framework/random/RNGManager.inl
index 8bc88135b1d0122e1668f85c34780d1104b91ece..749107cd9206adb3e3e3fc78b9e666f59f7b7f27 100644
--- a/corsika/detail/framework/random/RNGManager.inl
+++ b/corsika/detail/framework/random/RNGManager.inl
@@ -10,50 +10,42 @@
#pragma once
-
namespace corsika {
- void RNGManager::RegisterRandomStream(std::string const& pStreamName)
- {
- RNG rng;
+ inline void RNGManager::RegisterRandomStream(std::string const& pStreamName) {
+ RNG rng;
- if (auto const& it = seeds.find(pStreamName); it != seeds.end()) {
- rng.seed(it->second);
- }
+ if (auto const& it = seeds.find(pStreamName); it != seeds.end()) {
+ rng.seed(it->second);
+ }
- rngs[pStreamName] = std::move(rng);
- }
+ rngs[pStreamName] = std::move(rng);
+ }
- RNG& RNGManager::GetRandomStream
- (
- std::string const& pStreamName) {
- return rngs.at(pStreamName);
- }
+ inline RNG& RNGManager::GetRandomStream(std::string const& pStreamName) {
+ return rngs.at(pStreamName);
+ }
- std::stringstream RNGManager::dumpState() const
- {
- std::stringstream buffer;
- for (auto const& [streamName, rng] : rngs) {
- buffer << '"' << streamName << "\" = \"" << rng << '"' << std::endl;
- }
+ inline std::stringstream RNGManager::dumpState() const {
+ std::stringstream buffer;
+ for (auto const& [streamName, rng] : rngs) {
+ buffer << '"' << streamName << "\" = \"" << rng << '"' << std::endl;
+ }
- return buffer;
- }
+ return buffer;
+ }
- void RNGManager::SeedAll(uint64_t vSeed)
- {
- for (auto& entry : rngs) { entry.second.seed(vSeed++); }
- }
+ inline void RNGManager::SeedAll(uint64_t vSeed) {
+ for (auto& entry : rngs) { entry.second.seed(vSeed++); }
+ }
- void RNGManager::SeedAll()
- {
- std::random_device rd;
+ inline void RNGManager::SeedAll() {
+ std::random_device rd;
- for (auto& entry : rngs) {
- std::seed_seq sseq{rd(), rd(), rd(), rd(), rd(), rd()};
- entry.second.seed(sseq);
- }
- }
+ for (auto& entry : rngs) {
+ std::seed_seq sseq{rd(), rd(), rd(), rd(), rd(), rd()};
+ entry.second.seed(sseq);
+ }
+ }
} // namespace corsika
-
diff --git a/corsika/detail/framework/utility/COMBoost.inl b/corsika/detail/framework/utility/COMBoost.inl
index da63c03fe2e524fad74bb89c72723d3a873557d3..f0e9f9e4bd525738f0e42a30e070642246a0a3ca 100644
--- a/corsika/detail/framework/utility/COMBoost.inl
+++ b/corsika/detail/framework/utility/COMBoost.inl
@@ -10,117 +10,114 @@
#pragma once
+#include <Eigen/Dense>
+#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/CoordinateSystem.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/utility/sgn.hpp>
-#include <Eigen/Dense>
-#include <corsika/framework/core/PhysicalUnits.hpp>
namespace corsika {
- auto const& COMBoost::GetRotationMatrix() const { return fRotation; }
-
- //! transforms a 4-momentum from lab frame to the center-of-mass frame
- template <typename FourVector>
- FourVector COMBoost::toCoM(const FourVector& p) const {
- using namespace corsika::units::si;
- auto pComponents = p.GetSpaceLikeComponents().GetComponents(fCS);
- Eigen::Vector3d eVecRotated = fRotation * pComponents.eVector;
- Eigen::Vector2d lab;
-
- lab << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
- (eVecRotated(2) * (1 / 1_GeV).magnitude());
-
- auto const boostedZ = fBoost * lab;
- auto const E_CoM = boostedZ(0) * 1_GeV;
-
- eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
-
- return FourVector(E_CoM,
- corsika::Vector<corsika::units::si::hepmomentum_d>(fCS, eVecRotated));
+ auto const& COMBoost::GetRotationMatrix() const { return fRotation; }
+
+ //! transforms a 4-momentum from lab frame to the center-of-mass frame
+ template <typename FourVector>
+ inline FourVector COMBoost::toCoM(const FourVector& p) const {
+ using namespace corsika::units::si;
+ auto pComponents = p.GetSpaceLikeComponents().GetComponents(fCS);
+ Eigen::Vector3d eVecRotated = fRotation * pComponents.eVector;
+ Eigen::Vector2d lab;
+
+ lab << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
+ (eVecRotated(2) * (1 / 1_GeV).magnitude());
+
+ auto const boostedZ = fBoost * lab;
+ auto const E_CoM = boostedZ(0) * 1_GeV;
+
+ eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
+
+ return FourVector(
+ E_CoM, corsika::Vector<corsika::units::si::hepmomentum_d>(fCS, eVecRotated));
+ }
+
+ //! transforms a 4-momentum from the center-of-mass frame back to lab frame
+ template <typename FourVector>
+ inline FourVector COMBoost::fromCoM(const FourVector& p) const {
+ using namespace corsika::units::si;
+ Eigen::Vector2d com;
+ com << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
+ (p.GetSpaceLikeComponents().GetComponents().eVector(2) * (1 / 1_GeV).magnitude());
+
+ auto const plab = p.GetSpaceLikeComponents().GetComponents();
+ std::cout << "COMBoost::fromCoM Ecm=" << p.GetTimeLikeComponent() / 1_GeV << " GeV, "
+ << " pcm = " << plab / 1_GeV << " (norm = " << plab.norm() / 1_GeV
+ << " GeV), invariant mass = " << p.GetNorm() / 1_GeV << " GeV" << std::endl;
+
+ auto const boostedZ = fInverseBoost * com;
+ auto const E_lab = boostedZ(0) * 1_GeV;
+
+ auto pLab = p.GetSpaceLikeComponents().GetComponents();
+ pLab.eVector(2) = boostedZ(1) * (1_GeV).magnitude();
+ pLab.eVector = fRotation.transpose() * pLab.eVector;
+
+ FourVector f(E_lab, corsika::Vector(fCS, pLab));
+
+ std::cout << "COMBoost::fromCoM --> Elab=" << E_lab / 1_GeV << "GeV, "
+ << " pcm = " << pLab / 1_GeV << " (norm =" << pLab.norm() / 1_GeV
+ << " GeV), invariant mass = " << f.GetNorm() / 1_GeV << " GeV" << std::endl;
+
+ return f;
+ }
+
+ inline COMBoost::COMBoost(
+ FourVector<corsika::units::si::HEPEnergyType,
+ Vector<corsika::units::si::hepmomentum_d>> const& Pprojectile,
+ const corsika::units::si::HEPMassType massTarget)
+ : fCS(Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()) {
+ auto const pProjectile = Pprojectile.GetSpaceLikeComponents();
+ auto const pProjNorm = pProjectile.norm();
+ auto const a = (pProjectile / pProjNorm).GetComponents().eVector;
+ auto const a1 = a(0), a2 = a(1);
+
+ auto const s = sgn(a(2));
+ auto const c = 1 / (1 + s * a(2));
+
+ Eigen::Matrix3d A, B;
+
+ if (s > 0) {
+ A << 1, 0, -a1, // comment to prevent clang-format
+ 0, 1, -a2, // .
+ a1, a2, 1; // .
+ B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
+ -a1 * a2 * c, -a2 * a2 * c, 0, // .
+ 0, 0, -(a1 * a1 + a2 * a2) * c; // .
+
+ } else {
+ A << 1, 0, a1, // comment to prevent clang-format
+ 0, -1, -a2, // .
+ a1, a2, -1; // .
+ B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
+ +a1 * a2 * c, +a2 * a2 * c, 0, // .
+ 0, 0, (a1 * a1 + a2 * a2) * c; // .
}
- //! transforms a 4-momentum from the center-of-mass frame back to lab frame
- template <typename FourVector>
- FourVector COMBoost::fromCoM(const FourVector& p) const {
- using namespace corsika::units::si;
- Eigen::Vector2d com;
- com << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
- (p.GetSpaceLikeComponents().GetComponents().eVector(2) *
- (1 / 1_GeV).magnitude());
-
- auto const plab = p.GetSpaceLikeComponents().GetComponents();
- std::cout << "COMBoost::fromCoM Ecm=" << p.GetTimeLikeComponent() / 1_GeV
- << " GeV, "
- << " pcm = " << plab / 1_GeV << " (norm = " << plab.norm() / 1_GeV
- << " GeV), invariant mass = " << p.GetNorm() / 1_GeV << " GeV"
- << std::endl;
-
- auto const boostedZ = fInverseBoost * com;
- auto const E_lab = boostedZ(0) * 1_GeV;
-
- auto pLab = p.GetSpaceLikeComponents().GetComponents();
- pLab.eVector(2) = boostedZ(1) * (1_GeV).magnitude();
- pLab.eVector = fRotation.transpose() * pLab.eVector;
-
- FourVector f(E_lab, corsika::Vector(fCS, pLab));
-
- std::cout << "COMBoost::fromCoM --> Elab=" << E_lab / 1_GeV << "GeV, "
- << " pcm = " << pLab / 1_GeV << " (norm =" << pLab.norm() / 1_GeV
- << " GeV), invariant mass = " << f.GetNorm() / 1_GeV << " GeV"
- << std::endl;
-
- return f;
- }
-
- COMBoost::COMBoost(FourVector<corsika::units::si::HEPEnergyType, Vector<corsika::units::si::hepmomentum_d>> const& Pprojectile,
- const corsika::units::si::HEPMassType massTarget)
- : fCS(Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()) {
- auto const pProjectile = Pprojectile.GetSpaceLikeComponents();
- auto const pProjNorm = pProjectile.norm();
- auto const a = (pProjectile / pProjNorm).GetComponents().eVector;
- auto const a1 = a(0), a2 = a(1);
-
- auto const s = sgn(a(2));
- auto const c = 1 / (1 + s * a(2));
+ fRotation = A + B;
- Eigen::Matrix3d A, B;
+ // calculate boost
+ double const beta = pProjNorm / (Pprojectile.GetTimeLikeComponent() + massTarget);
- if (s > 0) {
- A << 1, 0, -a1, // comment to prevent clang-format
- 0, 1, -a2, // .
- a1, a2, 1; // .
- B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
- -a1 * a2 * c, -a2 * a2 * c, 0, // .
- 0, 0, -(a1 * a1 + a2 * a2) * c; // .
+ /* Accurracy matters here, beta = 1 - epsilon for ultra-relativistic boosts */
+ double const coshEta = 1 / std::sqrt((1 + beta) * (1 - beta));
+ //~ double const coshEta = 1 / std::sqrt((1-beta*beta));
+ double const sinhEta = -beta * coshEta;
- } else {
- A << 1, 0, a1, // comment to prevent clang-format
- 0, -1, -a2, // .
- a1, a2, -1; // .
- B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
- +a1 * a2 * c, +a2 * a2 * c, 0, // .
- 0, 0, (a1 * a1 + a2 * a2) * c; // .
- }
+ std::cout << "COMBoost (1-beta)=" << 1 - beta << " gamma=" << coshEta << std::endl;
+ std::cout << " det = " << fRotation.determinant() - 1 << std::endl;
- fRotation = A + B;
+ fBoost << coshEta, sinhEta, sinhEta, coshEta;
- // calculate boost
- double const beta = pProjNorm / (Pprojectile.GetTimeLikeComponent() + massTarget);
+ fInverseBoost << coshEta, -sinhEta, -sinhEta, coshEta;
+ }
- /* Accurracy matters here, beta = 1 - epsilon for ultra-relativistic boosts */
- double const coshEta = 1 / std::sqrt((1 + beta) * (1 - beta));
- //~ double const coshEta = 1 / std::sqrt((1-beta*beta));
- double const sinhEta = -beta * coshEta;
-
- std::cout << "COMBoost (1-beta)=" << 1 - beta << " gamma=" << coshEta << std::endl;
- std::cout << " det = " << fRotation.determinant() - 1 << std::endl;
-
- fBoost << coshEta, sinhEta, sinhEta, coshEta;
-
- fInverseBoost << coshEta, -sinhEta, -sinhEta, coshEta;
- }
} // namespace corsika
diff --git a/corsika/detail/modules/pythia8/Decay.inl b/corsika/detail/modules/pythia8/Decay.inl
new file mode 100644
index 0000000000000000000000000000000000000000..b407bee40862ac4dffecdbeaaf8e05392305e223
--- /dev/null
+++ b/corsika/detail/modules/pythia8/Decay.inl
@@ -0,0 +1,140 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/modules/pythia8/Pythia8.hpp>
+#include <corsika/modules/pythia8/Decay.hpp>
+#include <corsika/modules/pythia8/Random.hpp>
+
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+namespace corsika::pythia8 {
+
+ Decay::Decay(std::vector<corsika::Code> pParticles)
+ : fTrackedParticles(pParticles) {}
+
+ Decay::~Decay() { std::cout << "Pythia::Decay n=" << fCount << std::endl; }
+
+ void Decay::Init() {
+
+ Decay::SetParticleListStable(fTrackedParticles);
+
+ // set random number generator in pythia
+ Pythia8::RndmEngine* rndm = new corsika::pythia8::Random();
+ fPythia.setRndmEnginePtr(rndm);
+
+ fPythia.readString("Next:numberShowInfo = 0");
+ fPythia.readString("Next:numberShowProcess = 0");
+ fPythia.readString("Next:numberShowEvent = 0");
+
+ fPythia.readString("Print:quiet = on");
+
+ fPythia.readString("ProcessLevel:all = off");
+ fPythia.readString("ProcessLevel:resonanceDecays = off");
+
+ fPythia.particleData.readString("59:m0 = 101.00");
+
+ fPythia.init();
+ }
+
+ void Decay::SetParticleListStable(const std::vector<corsika::Code> particleList) {
+ for (auto p : particleList) Decay::SetStable(p);
+ }
+
+ void Decay::SetUnstable(const corsika::Code pCode) {
+ std::cout << "Pythia::Decay: setting " << pCode << " unstable.." << std::endl;
+ fPythia.particleData.mayDecay(static_cast<int>(corsika::GetPDG(pCode)), true);
+ }
+
+ void Decay::SetStable(const corsika::Code pCode) {
+ std::cout << "Pythia::Decay: setting " << pCode << " stable.." << std::endl;
+ fPythia.particleData.mayDecay(static_cast<int>(corsika::GetPDG(pCode)), false);
+ }
+
+ template <typename TParticle>
+ units::si::TimeType Decay::GetLifetime(TParticle const& p) {
+ using namespace units::si;
+
+ HEPEnergyType E = p.GetEnergy();
+ HEPMassType m = p.GetMass();
+
+ const double gamma = E / m;
+
+ const TimeType t0 = corsika::GetLifetime(p.GetPID());
+ auto const lifetime = gamma * t0;
+
+ return lifetime;
+ }
+
+ template <typename TProjectile>
+ void Decay::DoDecay(TProjectile& vP) {
+ using corsika::Point;
+ using namespace units::si;
+
+ auto const decayPoint = vP.GetPosition();
+ auto const t0 = vP.GetTime();
+
+ // coordinate system, get global frame of reference
+ corsika::CoordinateSystem& rootCS =
+ corsika::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+ fCount++;
+
+ // pythia stack
+ Pythia8::Event& event = fPythia.event;
+ event.reset();
+
+ // set particle unstable
+ Decay::SetUnstable(vP.GetPID());
+
+ // input particle PDG
+ auto const pdgCode = static_cast<int>(corsika::GetPDG(vP.GetPID()));
+
+ auto const pcomp = vP.GetMomentum().GetComponents();
+ double px = pcomp[0] / 1_GeV;
+ double py = pcomp[1] / 1_GeV;
+ double pz = pcomp[2] / 1_GeV;
+ double en = vP.GetEnergy() / 1_GeV;
+ double m = corsika::GetMass(vP.GetPID()) / 1_GeV;
+
+ // add particle to pythia stack
+ event.append(pdgCode, 1, 0, 0, px, py, pz, en, m);
+
+ if (!fPythia.next())
+ std::cout << "Pythia::Decay: decay failed!" << std::endl;
+ else
+ std::cout << "Pythia::Decay: particles after decay: " << event.size() << std::endl;
+
+ // list final state
+ event.list();
+
+ // loop over final state
+ for (int i = 0; i < event.size(); ++i)
+ if (event[i].isFinal()) {
+ auto const pyId =
+ corsika::ConvertFromPDG(static_cast<corsika::PDGCode>(event[i].id()));
+ HEPEnergyType pyEn = event[i].e() * 1_GeV;
+ MomentumVector pyP(rootCS, {event[i].px() * 1_GeV, event[i].py() * 1_GeV,
+ event[i].pz() * 1_GeV});
+
+ std::cout << "particle: id=" << pyId << " momentum=" << pyP.GetComponents() / 1_GeV
+ << " energy=" << pyEn << std::endl;
+
+ vP.AddSecondary(
+ std::tuple<corsika::Code, units::si::HEPEnergyType,
+ corsika::MomentumVector, corsika::Point, units::si::TimeType>{
+ pyId, pyEn, pyP, decayPoint, t0});
+ }
+
+ // set particle stable
+ Decay::SetStable(vP.GetPID());
+ }
+
+} // namespace corsika::process::pythia
diff --git a/dependencies/Pythia/Interaction.cc b/corsika/detail/modules/pythia8/Interaction.inl
similarity index 59%
rename from dependencies/Pythia/Interaction.cc
rename to corsika/detail/modules/pythia8/Interaction.inl
index dbe514cf272059e48ac0ddcdb8ce6a303d39dc51..33b649912723afe43e35b9f5698d93fc836ef27a 100644
--- a/dependencies/Pythia/Interaction.cc
+++ b/corsika/detail/modules/pythia8/Interaction.inl
@@ -1,40 +1,39 @@
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
+ * See file AUTHORS for a list of contributors.
+ *
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
-#include <dependencies/Pythia/Interaction.h>
+#pragma once
+
+#include <corsika/modules/pythia8/Interaction.hpp>
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/NuclearComposition.h>
#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
#include <corsika/setup/SetupStack.hpp>
#include <tuple>
-using std::cout;
-using std::endl;
-using std::tuple;
-
-using namespace corsika;
-using namespace corsika::setup;
-using SecondaryView = corsika::setup::StackView;
+using Projectile = corsika::setup::StackView::ParticleType;
using Particle = corsika::setup::Stack::ParticleType;
-namespace corsika::pythia {
+namespace corsika::pythia8 {
typedef corsika::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
- Interaction::~Interaction() {}
+ Interaction::~Interaction() {
+ std::cout << "Pythia::Interaction n=" << fCount << std::endl;
+ }
- Interaction::Interaction() {
- cout << "Pythia::Interaction n=" << fCount << endl;
+ void Interaction::Init() {
- using random::RNGManager;
+ using corsika::RNGManager;
// initialize Pythia
if (!fInitialized) {
@@ -50,12 +49,16 @@ namespace corsika::pythia {
if (fInternalDecays) {
// define which particles are passed to corsika, i.e. which particles make it into
// history even very shortlived particles like charm or pi0 are of interest here
- const std::vector<Code> HadronsWeWantTrackedByCorsika = {
- Code::PiPlus, Code::PiMinus, Code::Pi0, Code::KMinus,
- Code::KPlus, Code::K0Long, Code::K0Short, Code::SigmaPlus,
- Code::SigmaMinus, Code::Lambda0, Code::Xi0, Code::XiMinus,
- Code::OmegaMinus, Code::DPlus, Code::DMinus, Code::D0,
- Code::D0Bar};
+ const std::vector<corsika::Code> HadronsWeWantTrackedByCorsika = {
+ corsika::Code::PiPlus, corsika::Code::PiMinus,
+ corsika::Code::Pi0, corsika::Code::KMinus,
+ corsika::Code::KPlus, corsika::Code::K0Long,
+ corsika::Code::K0Short, corsika::Code::SigmaPlus,
+ corsika::Code::SigmaMinus, corsika::Code::Lambda0,
+ corsika::Code::Xi0, corsika::Code::XiMinus,
+ corsika::Code::OmegaMinus, corsika::Code::DPlus,
+ corsika::Code::DMinus, corsika::Code::D0,
+ corsika::Code::D0Bar};
Interaction::SetParticleListStable(HadronsWeWantTrackedByCorsika);
}
@@ -67,36 +70,39 @@ namespace corsika::pythia {
}
}
- void Interaction::SetParticleListStable(std::vector<Code> const& particleList) {
+ void Interaction::SetParticleListStable(
+ std::vector<corsika::Code> const& particleList) {
for (auto p : particleList) Interaction::SetStable(p);
}
- void Interaction::SetUnstable(const Code pCode) {
- cout << "Pythia::Interaction: setting " << pCode << " unstable.." << endl;
- fPythia.particleData.mayDecay(static_cast<int>(GetPDG(pCode)), true);
+ void Interaction::SetUnstable(const corsika::Code pCode) {
+ std::cout << "Pythia::Interaction: setting " << pCode << " unstable.." << std::endl;
+ fPythia.particleData.mayDecay(static_cast<int>(corsika::GetPDG(pCode)), true);
}
- void Interaction::SetStable(const Code pCode) {
- cout << "Pythia::Interaction: setting " << pCode << " stable.." << endl;
- fPythia.particleData.mayDecay(static_cast<int>(GetPDG(pCode)), false);
+ void Interaction::SetStable(const corsika::Code pCode) {
+ std::cout << "Pythia::Interaction: setting " << pCode << " stable.." << std::endl;
+ fPythia.particleData.mayDecay(static_cast<int>(corsika::GetPDG(pCode)), false);
}
void Interaction::ConfigureLabFrameCollision(
- const Code BeamId, const Code TargetId, const units::si::HEPEnergyType BeamEnergy) {
+ const corsika::Code BeamId, const corsika::Code TargetId,
+ const units::si::HEPEnergyType BeamEnergy) {
using namespace units::si;
// Pythia configuration of the current event
// very clumsy. I am sure this can be done better..
// set beam
// beam id for pythia
- auto const pdgBeam = static_cast<int>(GetPDG(BeamId));
+ auto const pdgBeam = static_cast<int>(corsika::GetPDG(BeamId));
std::stringstream stBeam;
stBeam << "Beams:idA = " << pdgBeam;
fPythia.readString(stBeam.str());
// set target
- auto pdgTarget = static_cast<int>(GetPDG(TargetId));
+ auto pdgTarget = static_cast<int>(corsika::GetPDG(TargetId));
// replace hydrogen with proton, otherwise pythia goes into heavy ion mode!
- if (TargetId == Code::Hydrogen) pdgTarget = static_cast<int>(GetPDG(Code::Proton));
+ if (TargetId == corsika::Code::Hydrogen)
+ pdgTarget = static_cast<int>(corsika::GetPDG(corsika::Code::Proton));
std::stringstream stTarget;
stTarget << "Beams:idB = " << pdgTarget;
fPythia.readString(stTarget.str());
@@ -119,17 +125,17 @@ namespace corsika::pythia {
pCode == corsika::Code::PiMinus || pCode == corsika::Code::PiPlus;
}
- tuple<units::si::CrossSectionType, units::si::CrossSectionType>
- Interaction::GetCrossSection(const Code BeamId, const Code TargetId,
+ std::tuple<units::si::CrossSectionType, units::si::CrossSectionType>
+ Interaction::GetCrossSection(const corsika::Code BeamId, const corsika::Code TargetId,
const units::si::HEPEnergyType CoMenergy) {
using namespace units::si;
// interaction possible in pythia?
- if (TargetId == Code::Proton || TargetId == Code::Hydrogen) {
+ if (TargetId == corsika::Code::Proton || TargetId == corsika::Code::Hydrogen) {
if (CanInteract(BeamId) && ValidCoMEnergy(CoMenergy)) {
// input particle PDG
- auto const pdgCodeBeam = static_cast<int>(GetPDG(BeamId));
- auto const pdgCodeTarget = static_cast<int>(GetPDG(TargetId));
+ auto const pdgCodeBeam = static_cast<int>(corsika::GetPDG(BeamId));
+ auto const pdgCodeTarget = static_cast<int>(corsika::GetPDG(TargetId));
const double ecm = CoMenergy / 1_GeV;
// calculate cross section
@@ -157,24 +163,23 @@ namespace corsika::pythia {
using namespace units;
using namespace units::si;
- using namespace geometry;
// coordinate system, get global frame of reference
CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
- const Code corsikaBeamId = p.GetPID();
+ const corsika::Code corsikaBeamId = p.GetPID();
// beam particles for pythia : 1, 2, 3 for p, pi, k
// read from cross section code table
const bool kInteraction = CanInteract(corsikaBeamId);
// FOR NOW: assume target is at rest
- process::pythia::MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
+ corsika::MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy
HEPEnergyType Elab = p.GetEnergy() + constants::nucleonMass;
- process::pythia::MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
+ corsika::MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
pTotLab += p.GetMomentum();
pTotLab += pTarget;
auto const pTotLabNorm = pTotLab.norm();
@@ -182,10 +187,10 @@ namespace corsika::pythia {
const HEPEnergyType ECoM = sqrt(
(Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
- cout << "Interaction: LambdaInt: \n"
- << " input energy: " << p.GetEnergy() / 1_GeV << endl
- << " beam can interact:" << kInteraction << endl
- << " beam pid:" << p.GetPID() << endl;
+ std::cout << "Interaction: LambdaInt: \n"
+ << " input energy: " << p.GetEnergy() / 1_GeV << std::endl
+ << " beam can interact:" << kInteraction << std::endl
+ << " beam pid:" << p.GetPID() << std::endl;
// TODO: move limits into variables
if (kInteraction && Elab >= 8.5_GeV && ValidCoMEnergy(ECoM)) {
@@ -206,17 +211,17 @@ namespace corsika::pythia {
return std::get<0>(this->GetCrossSection(corsikaBeamId, vTargetID, ECoM));
});
- cout << "Interaction: IntLength: weighted CrossSection (mb): "
- << weightedProdCrossSection / 1_mb << endl
- << "Interaction: IntLength: average mass number: "
- << mediumComposition.GetAverageMassNumber() << endl;
+ std::cout << "Interaction: IntLength: weighted CrossSection (mb): "
+ << weightedProdCrossSection / 1_mb << std::endl
+ << "Interaction: IntLength: average mass number: "
+ << mediumComposition.GetAverageMassNumber() << std::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;
+ std::cout << "Interaction: "
+ << "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
+ << std::endl;
return int_length;
}
@@ -230,33 +235,29 @@ namespace corsika::pythia {
*/
template <>
- process::EProcessReturn Interaction::DoInteraction(SecondaryView& view) {
+ corsika::EProcessReturn Interaction::DoInteraction(Projectile& vP) {
using namespace units;
- using namespace utl;
using namespace units::si;
- using namespace geometry;
-
- auto const projectile = view.GetProjectile();
- const auto corsikaBeamId = projectile.GetPID();
- cout << "Pythia::Interaction: "
- << "DoInteraction: " << corsikaBeamId << " interaction? "
- << process::pythia::Interaction::CanInteract(corsikaBeamId) << endl;
+ const auto corsikaBeamId = vP.GetPID();
+ std::cout << "Pythia::Interaction: "
+ << "DoInteraction: " << corsikaBeamId << " interaction? "
+ << corsika::pythia8::Interaction::CanInteract(corsikaBeamId) << std::endl;
- if (IsNucleus(corsikaBeamId)) {
+ if (corsika::IsNucleus(corsikaBeamId)) {
// nuclei handled by different process, this should not happen
throw std::runtime_error("Nuclear projectile are not handled by PYTHIA!");
}
- if (process::pythia::Interaction::CanInteract(corsikaBeamId)) {
+ if (corsika::pythia8::Interaction::CanInteract(corsikaBeamId)) {
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// position and time of interaction, not used in Sibyll
- Point pOrig = projectile.GetPosition();
- TimeType tOrig = projectile.GetTime();
+ Point pOrig = vP.GetPosition();
+ TimeType tOrig = vP.GetTime();
// define target
// FOR NOW: target is always at rest
@@ -265,14 +266,15 @@ namespace corsika::pythia {
const FourVector PtargLab(eTargetLab, pTargetLab);
// define projectile
- HEPEnergyType const eProjectileLab = projectile.GetEnergy();
- auto const pProjectileLab = projectile.GetMomentum();
+ HEPEnergyType const eProjectileLab = vP.GetEnergy();
+ auto const pProjectileLab = vP.GetMomentum();
- cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << endl
- << "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
- << endl;
- cout << "Interaction: etarget lab: " << eTargetLab / 1_GeV << endl
- << "Interaction: ptarget lab: " << pTargetLab.GetComponents() / 1_GeV << endl;
+ std::cout << "Interaction: ebeam lab: " << eProjectileLab / 1_GeV << std::endl
+ << "Interaction: pbeam lab: " << pProjectileLab.GetComponents() / 1_GeV
+ << std::endl;
+ std::cout << "Interaction: etarget lab: " << eTargetLab / 1_GeV << std::endl
+ << "Interaction: ptarget lab: " << pTargetLab.GetComponents() / 1_GeV
+ << std::endl;
const FourVector PprojLab(eProjectileLab, pProjectileLab);
@@ -288,25 +290,26 @@ namespace corsika::pythia {
// boost target
auto const PtargCoM = boost.toCoM(PtargLab);
- cout << "Interaction: ebeam CoM: " << PprojCoM.GetTimeLikeComponent() / 1_GeV
- << endl
- << "Interaction: pbeam CoM: "
- << PprojCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
- cout << "Interaction: etarget CoM: " << PtargCoM.GetTimeLikeComponent() / 1_GeV
- << endl
- << "Interaction: ptarget CoM: "
- << PtargCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl;
+ std::cout << "Interaction: ebeam CoM: " << PprojCoM.GetTimeLikeComponent() / 1_GeV
+ << std::endl
+ << "Interaction: pbeam CoM: "
+ << PprojCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << std::endl;
+ std::cout << "Interaction: etarget CoM: " << PtargCoM.GetTimeLikeComponent() / 1_GeV
+ << std::endl
+ << "Interaction: ptarget CoM: "
+ << PtargCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << std::endl;
- cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl;
- cout << "Interaction: time: " << tOrig << endl;
+ std::cout << "Interaction: position of interaction: " << pOrig.GetCoordinates()
+ << std::endl;
+ std::cout << "Interaction: time: " << tOrig << std::endl;
HEPEnergyType Etot = eProjectileLab + eTargetLab;
- MomentumVector Ptot = projectile.GetMomentum();
+ MomentumVector Ptot = vP.GetMomentum();
// invariant mass, i.e. cm. energy
HEPEnergyType Ecm = sqrt(Etot * Etot - Ptot.squaredNorm());
// sample target mass number
- const auto* currentNode = projectile.GetNode();
+ const auto* currentNode = vP.GetNode();
const auto& mediumComposition =
currentNode->GetModelProperties().GetNuclearComposition();
// get cross sections for target materials
@@ -327,20 +330,21 @@ namespace corsika::pythia {
const auto corsikaTargetId =
mediumComposition.SampleTarget(cross_section_of_components, fRNG);
- cout << "Interaction: target selected: " << corsikaTargetId << endl;
+ std::cout << "Interaction: target selected: " << corsikaTargetId << std::endl;
- if (corsikaTargetId != Code::Hydrogen && corsikaTargetId != Code::Neutron &&
- corsikaTargetId != Code::Proton)
+ if (corsikaTargetId != corsika::Code::Hydrogen &&
+ corsikaTargetId != corsika::Code::Neutron &&
+ corsikaTargetId != corsika::Code::Proton)
throw std::runtime_error("DoInteraction: wrong target for PYTHIA");
- cout << "Interaction: "
- << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
- << " Ecm(GeV): " << Ecm / 1_GeV << endl;
+ std::cout << "Interaction: "
+ << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
+ << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
if (eProjectileLab < 8.5_GeV || !ValidCoMEnergy(Ecm)) {
- cout << "Interaction: "
- << " DoInteraction: should have dropped particle.. "
- << "THIS IS AN ERROR" << endl;
+ std::cout << "Interaction: "
+ << " DoInteraction: should have dropped particle.. "
+ << "THIS IS AN ERROR" << std::endl;
throw std::runtime_error("energy too low for PYTHIA");
} else {
@@ -363,27 +367,28 @@ namespace corsika::pythia {
// skip particles that have decayed in pythia
if (!p8p.isFinal()) continue;
- auto const pyId = ConvertFromPDG(static_cast<PDGCode>(p8p.id()));
+ auto const pyId =
+ corsika::ConvertFromPDG(static_cast<corsika::PDGCode>(p8p.id()));
const MomentumVector pyPlab(
rootCS, {p8p.px() * 1_GeV, p8p.py() * 1_GeV, p8p.pz() * 1_GeV});
HEPEnergyType const pyEn = p8p.e() * 1_GeV;
// add to corsika stack
- auto pnew = view.AddSecondary(
- tuple<Code, units::si::HEPEnergyType, stack::MomentumVector,
- geometry::Point, units::si::TimeType>{pyId, pyEn, pyPlab, pOrig,
- tOrig});
+ auto pnew = vP.AddSecondary(
+ std::tuple<corsika::Code, units::si::HEPEnergyType, corsika::MomentumVector,
+ corsika::Point, units::si::TimeType>{pyId, pyEn, pyPlab, pOrig,
+ tOrig});
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
}
- cout << "conservation (all GeV): "
- << "Elab_final=" << Elab_final / 1_GeV
- << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
+ std::cout << "conservation (all GeV): "
+ << "Elab_final=" << Elab_final / 1_GeV
+ << ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << std::endl;
}
}
- return process::EProcessReturn::eOk;
+ return EProcessReturn::eOk;
}
-} // namespace corsika::pythia
+} // namespace corsika::pythia8
diff --git a/dependencies/Pythia/Random.cc b/corsika/detail/modules/pythia8/Random.inl
similarity index 63%
rename from dependencies/Pythia/Random.cc
rename to corsika/detail/modules/pythia8/Random.inl
index 099a1fd2a78b0437df4059f06ff2b17c2d9f94e2..cfd19bb19cb768a79c303760e34733f5e23f0f15 100644
--- a/dependencies/Pythia/Random.cc
+++ b/corsika/detail/modules/pythia8/Random.inl
@@ -1,15 +1,19 @@
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
+ * See file AUTHORS for a list of contributors.
+ *
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
-#include <dependencies/Pythia/Random.h>
+#pragma once
+
+#include <corsika/modules/pythia8/Random.hpp>
-namespace corsika::pythia {
+namespace corsika::pythia8 {
double Random::flat() { return fDist(fRNG); }
-} // namespace corsika::pythia
+} // namespace corsika::pythia8
diff --git a/corsika/framework/utility/COMBoost.hpp b/corsika/framework/utility/COMBoost.hpp
index 04232fe6eb88d8cfbe15c1685d800f80302f99ae..17e15352b8f82ff6ef3bdbcb02682465fdae3b17 100644
--- a/corsika/framework/utility/COMBoost.hpp
+++ b/corsika/framework/utility/COMBoost.hpp
@@ -37,15 +37,15 @@ namespace corsika {
corsika::Vector<corsika::units::si::hepmomentum_d>>& Pprojectile,
const corsika::units::si::HEPEnergyType massTarget);
- auto const& GetRotationMatrix() const;
+ inline auto const& GetRotationMatrix() const;
//! transforms a 4-momentum from lab frame to the center-of-mass frame
template <typename FourVector>
- FourVector toCoM(const FourVector& p) const ;
+ inline FourVector toCoM(const FourVector& p) const ;
//! transforms a 4-momentum from the center-of-mass frame back to lab frame
template <typename FourVector>
- FourVector fromCoM(const FourVector& p) const ;
+ inline FourVector fromCoM(const FourVector& p) const ;
};
} // namespace corsika
diff --git a/corsika/modules/pythia8/Decay.hpp b/corsika/modules/pythia8/Decay.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d4b35c869be5865323ef784b88f99dfd582701fe
--- /dev/null
+++ b/corsika/modules/pythia8/Decay.hpp
@@ -0,0 +1,48 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/sequence/DecayProcess.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+
+#include <corsika/modules/pythia8/Pythia8.hpp>
+
+
+namespace corsika::pythia8 {
+
+ typedef corsika::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
+
+ class Decay : public corsika::DecayProcess<Decay> {
+ const std::vector<corsika::Code> fTrackedParticles;
+ int fCount = 0;
+
+ public:
+ Decay(std::vector<corsika::Code>);
+ ~Decay();
+ void Init();
+
+ void SetParticleListStable(const std::vector<corsika::Code>);
+ void SetUnstable(const corsika::Code);
+ void SetStable(const corsika::Code);
+
+ template <typename TParticle>
+ corsika::units::si::TimeType GetLifetime(TParticle const&);
+
+ template <typename TProjectile>
+ void DoDecay(TProjectile&);
+
+ private:
+ Pythia8::Pythia fPythia;
+ };
+
+} // namespace corsika::pythia8
+
+#include <corsika/detail/modules/pythia8/Decay.inl>
diff --git a/dependencies/Pythia/Interaction.h b/corsika/modules/pythia8/Interaction.hpp
similarity index 68%
rename from dependencies/Pythia/Interaction.h
rename to corsika/modules/pythia8/Interaction.hpp
index c5006ea6804d61f9c65ec27cfbf4d824e01712f9..ddc488d4a6d0e021daa0f3d0286d0f137bacf0d9 100644
--- a/dependencies/Pythia/Interaction.h
+++ b/corsika/modules/pythia8/Interaction.hpp
@@ -1,6 +1,8 @@
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
+ * See file AUTHORS for a list of contributors.
+ *
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
@@ -8,15 +10,15 @@
#pragma once
-#include <Pythia8/Pythia.h>
-
+#include <corsika/modules/pythia8/Pythia8.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/sequence/InteractionProcess.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
#include <tuple>
-namespace corsika::pythia {
+namespace corsika::pythia8 {
class Interaction : public corsika::InteractionProcess<Interaction> {
@@ -24,10 +26,12 @@ namespace corsika::pythia {
bool fInitialized = false;
public:
- Interaction();
+ Interaction() {}
~Interaction();
- void SetParticleListStable(std::vector<Code> const&);
+ void Init();
+
+ void SetParticleListStable(std::vector<corsika::Code> const&);
void SetUnstable(const corsika::Code);
void SetStable(const corsika::Code);
@@ -38,10 +42,12 @@ namespace corsika::pythia {
}
bool CanInteract(const corsika::Code);
- void ConfigureLabFrameCollision(const corsika::Code, const corsika::Code,
+ void ConfigureLabFrameCollision(const corsika::Code,
+ const corsika::Code,
const corsika::units::si::HEPEnergyType);
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
- GetCrossSection(const corsika::Code BeamId, const corsika::Code TargetId,
+ GetCrossSection(const corsika::Code BeamId,
+ const corsika::Code TargetId,
const corsika::units::si::HEPEnergyType CoMenergy);
template <typename TParticle>
@@ -52,15 +58,17 @@ namespace corsika::pythia {
event is copied (and boosted) into the shower lab frame.
*/
- template <typename TSecondaryView>
- corsika::EProcessReturn DoInteraction(TSecondaryView&);
+ template <typename TProjectile>
+ corsika::EProcessReturn DoInteraction(TProjectile&);
private:
- corsika::random::RNG& fRNG =
- corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
+ corsika::RNG& fRNG =
+ corsika::RNGManager::GetInstance().GetRandomStream("pythia");
Pythia8::Pythia fPythia;
Pythia8::SigmaTotal fSigma;
const bool fInternalDecays = true;
};
-} // namespace corsika::pythia
+} // namespace corsika::pythia8
+
+#include <corsika/detail/modules/pythia8/Interaction.inl>
diff --git a/corsika/modules/pythia8/Pythia8.hpp b/corsika/modules/pythia8/Pythia8.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4ea2cef62885cd9cac3a2869cdc753c7d4c1738e
--- /dev/null
+++ b/corsika/modules/pythia8/Pythia8.hpp
@@ -0,0 +1,3 @@
+#pragma once
+
+#include <Pythia8/Pythia.h>
diff --git a/dependencies/Pythia/Random.h b/corsika/modules/pythia8/Random.hpp
similarity index 64%
rename from dependencies/Pythia/Random.h
rename to corsika/modules/pythia8/Random.hpp
index 59bee7a356b5afbeb425f843bcd1dafff98e25ff..b7294bcb4756cedcee6505a64b6b62206feaed44 100644
--- a/dependencies/Pythia/Random.h
+++ b/corsika/modules/pythia8/Random.hpp
@@ -1,6 +1,8 @@
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
+ * See file AUTHORS for a list of contributors.
+ *
* This software is distributed under the terms of the GNU General Public
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
@@ -9,20 +11,21 @@
#pragma once
#include <Pythia8/Pythia.h>
-#include <corsika/framework/random/RNGManager.hpp>
-namespace corsika::process {
+#include <corsika/framework/random/RNGManager.hpp>
- namespace pythia {
+namespace corsika::pythia8 {
class Random : public Pythia8::RndmEngine {
double flat();
private:
std::uniform_real_distribution<double> fDist;
- corsika::random::RNG& fRNG =
- corsika::random::RNGManager::GetInstance().GetRandomStream("pythia");
+ corsika::RNG& fRNG =
+ corsika::RNGManager::GetInstance().GetRandomStream("pythia");
};
- } // namespace pythia
-} // namespace corsika::process
+} // namespace corsika::pythia8
+
+
+#include <corsika/detail/modules/pythia8/Random.inl>
diff --git a/dependencies/Pythia/CMakeLists.txt b/dependencies/Pythia/CMakeLists.txt
deleted file mode 100644
index 516afe3af6a4c81b3302b97082a35c0efcc4fa4e..0000000000000000000000000000000000000000
--- a/dependencies/Pythia/CMakeLists.txt
+++ /dev/null
@@ -1,73 +0,0 @@
-set (
- MODEL_SOURCES
- Decay.cc
- Random.cc
- Interaction.cc
- )
-
-set (
- MODEL_HEADERS
- Decay.h
- Random.h
- Interaction.h
- )
-
-set (
- MODEL_NAMESPACE
- dependencies/Pythia
- )
-
-add_library (ProcessPythia8 STATIC ${MODEL_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessPythia8 ${MODEL_NAMESPACE} ${MODEL_HEADERS})
-
-
-set_target_properties (
- ProcessPythia8
- PROPERTIES
- VERSION ${PROJECT_VERSION}
- SOVERSION 1
- )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
- ProcessPythia8
- CORSIKAprocesssequence
- CORSIKAparticles
- CORSIKAutilities
- CORSIKAunits
- CORSIKAthirdparty
- CORSIKAgeometry
- CORSIKAenvironment
- CORSIKAsetup
- CORSIKArandom
- C8::ext::pythia8
- )
-
-target_include_directories (
- ProcessPythia8
- INTERFACE
- $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
- $<INSTALL_INTERFACE:include/include>
- )
-
-install (
- TARGETS ProcessPythia8
- LIBRARY DESTINATION lib
- ARCHIVE DESTINATION lib
- )
-
-
-# --------------------
-# code unit testing
-
-CORSIKA_ADD_TEST(testPythia8
- SOURCES
- testPythia8.cc
- ${MODEL_HEADERS}
-)
-target_link_libraries (
- testPythia8
- ProcessPythia8
- CORSIKAtesting
- )
-
diff --git a/dependencies/Pythia/Decay.cc b/dependencies/Pythia/Decay.cc
deleted file mode 100644
index 690593913f9969814dcc0eb8d894019210235744..0000000000000000000000000000000000000000
--- a/dependencies/Pythia/Decay.cc
+++ /dev/null
@@ -1,262 +0,0 @@
-/*
- * (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 <Pythia8/Pythia.h>
-#include <dependencies/Pythia/Decay.h>
-#include <dependencies/Pythia/Random.h>
-
-#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/framework/utility/COMBoost.hpp>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-using std::cout;
-using std::endl;
-using std::tuple;
-using std::vector;
-
-using namespace corsika;
-using namespace corsika::setup;
-using View = corsika::StackView;
-using Particle = corsika::Stack::ParticleType;
-using Track = Trajectory;
-
-namespace corsika::pythia {
-
- Decay::Decay() {
-
- // set random number generator in pythia
- Pythia8::RndmEngine* rndm = new corsika::pythia::Random();
- fPythia.setRndmEnginePtr(rndm);
-
- /*
- issue xyz: definition of particles and decay channels use the same mechanism in
- corsika and pythia we should force pythia to use the file in corsika.
- */
- // bool ParticleData::reInit(string startFile, bool xmlFormat = true)
- // read in particle data from Corsika 8
- // fPythia.particleData.reInit("/home/felix/ngcorsika/corsika-build/include/corsika/framework/coreParticleData.xml");
- // fPythia.particleData.checkTable();
-
- fPythia.readString("Next:numberShowInfo = 0");
- fPythia.readString("Next:numberShowProcess = 0");
- fPythia.readString("Next:numberShowEvent = 0");
-
- fPythia.readString("Print:quiet = on");
- fPythia.readString("Check:particleData = 0");
-
- /*
- switching off event check in pythia is needed to allow decays that are off-shell
- according to the mass definition in pythia.
- the consistency of particle masses between event generators is an unsolved issues
- */
- cout << "Pythia::Init: switching off event checking in pythia.." << endl;
- fPythia.readString("Check:event = 1");
-
- fPythia.readString("ProcessLevel:all = off");
- fPythia.readString("ProcessLevel:resonanceDecays = off");
-
- // making sure
- SetStable(Code::Pi0);
-
- // fPythia.particleData.readString("59:m0 = 101.00");
-
- if (!fPythia.init())
- throw std::runtime_error("Pythia::Decay: Initialization failed!");
- }
-
- Decay::Decay(std::set<Code> vHandled)
- : handleAllDecays_(false)
- , handledDecays_(vHandled) {}
-
- Decay::~Decay() { cout << "Pythia::Decay n=" << fCount << endl; }
-
- bool Decay::CanHandleDecay(const Code vParticleCode) {
- using namespace corsika::particles;
- // if known to pythia and not proton, electron 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 (CanDecay(vParticleCode)) // non-zero for particles known to sibyll
- return true;
- else
- return false;
- }
-
- void Decay::SetHandleDecay(const Code vParticleCode) {
- handleAllDecays_ = false;
- cout << "Pythia::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 pythia!");
- }
-
- void Decay::SetHandleDecay(const vector<Code> vParticleList) {
- handleAllDecays_ = false;
- for (auto p : vParticleList) SetHandleDecay(p);
- }
-
- bool Decay::IsDecayHandled(const corsika::Code vParticleCode) {
- if (handleAllDecays_ && CanHandleDecay(vParticleCode))
- return true;
- else
- return handledDecays_.find(vParticleCode) != Decay::handledDecays_.end();
- }
-
- bool Decay::IsStable(const Code vCode) {
- return fPythia.particleData.canDecay(static_cast<int>(GetPDG(vCode)));
- }
-
- void Decay::PrintDecayConfig(const Code vCode) {
- cout << "Decay: Pythia decay configuration:" << endl;
- cout << vCode << " is ";
- if (IsStable(vCode))
- cout << "stable" << endl;
- else
- cout << "unstable" << endl;
- }
-
- void Decay::PrintDecayConfig() {
- cout << "Pythia::Decay: decay configuration:" << endl;
- if (handleAllDecays_)
- cout << " all particles known to Pythia are handled by Pythia::Decay!" << endl;
- else
- for (auto& pCode : handledDecays_)
- cout << "Decay of " << pCode << " is handled by Pythia!" << endl;
- }
-
- void Decay::SetStable(const vector<Code> particleList) {
- for (auto p : particleList) Decay::SetStable(p);
- }
-
- bool Decay::CanDecay(const Code pCode) {
- std::cout << "Pythia::Decay: checking if particle: " << pCode
- << " can decay in PYTHIA? ";
- const bool ans = fPythia.particleData.canDecay(static_cast<int>(GetPDG(pCode)));
- std::cout << ans << std::endl;
- return ans;
- }
-
- void Decay::SetUnstable(const Code pCode) {
- cout << "Pythia::Decay: setting " << pCode << " unstable.." << endl;
- fPythia.particleData.mayDecay(static_cast<int>(GetPDG(pCode)), true);
- }
-
- void Decay::SetStable(const Code pCode) {
- cout << "Pythia::Decay: setting " << pCode << " stable.." << endl;
- fPythia.particleData.mayDecay(static_cast<int>(GetPDG(pCode)), false);
- }
-
- template <>
- units::si::TimeType Decay::GetLifetime(Particle const& vP) {
- using namespace units::si;
-
- const auto pid = vP.GetPID();
- if (CanDecay(pid)) {
- HEPEnergyType E = vP.GetEnergy();
- HEPMassType m = vP.GetMass();
-
- const double gamma = E / m;
-
- const TimeType t0 = GetLifetime(pid);
- auto const lifetime = gamma * t0;
- cout << "Pythia::Decay: code: " << vP.GetPID() << endl;
- cout << "Pythia::Decay: MinStep: t0: " << t0 << endl;
- cout << "Pythia::Decay: MinStep: energy: " << E / 1_GeV << " GeV" << endl;
- cout << "Pythia::Decay: momentum: " << vP.GetMomentum().GetComponents() / 1_GeV
- << " GeV" << endl;
- cout << "Pythia::Decay: MinStep: gamma: " << gamma << endl;
- cout << "Pythia::Decay: MinStep: tau: " << lifetime << endl;
-
- return lifetime;
- } else
- return std::numeric_limits<double>::infinity() * 1_s;
- }
-
- template <>
- void Decay::DoDecay(View& view) {
- using geometry::Point;
- using namespace units;
- using namespace units::si;
-
- auto const projectile = view.GetProjectile();
-
- auto const& decayPoint = projectile.GetPosition();
- auto const t0 = projectile.GetTime();
-
- auto const& labMomentum = projectile.GetMomentum();
- geometry::CoordinateSystem const& labCS = labMomentum.GetCoordinateSystem();
-
- // define target kinematics in lab frame
- // define boost to and from CoM frame
- // CoM frame definition in Pythia projectile: +z
- utl::COMBoost const boost(labMomentum, projectile.GetMass());
- auto const& rotatedCS = boost.GetRotatedCS();
-
- fCount++;
-
- // pythia stack
- Pythia8::Event& event = fPythia.event;
- event.reset();
-
- auto const particleId = projectile.GetPID();
-
- // set particle unstable
- Decay::SetUnstable(particleId);
-
- // input particle PDG
- auto const pdgCode = static_cast<int>(GetPDG(particleId));
-
- double constexpr px = 0;
- double constexpr py = 0;
- double constexpr pz = 0;
- double const en = projectile.GetMass() / 1_GeV;
- double const m = en;
-
- // add particle to pythia stack
- event.append(pdgCode, 1, 0, 0, px, py, pz, en, m);
-
- if (!fPythia.next())
- throw std::runtime_error("Pythia::Decay: decay failed!");
- else
- cout << "Pythia::Decay: particles after decay: " << event.size() << endl;
-
- // list final state
- event.list();
-
- // loop over final state
- for (int i = 0; i < event.size(); ++i)
- if (event[i].isFinal()) {
- auto const pyId = ConvertFromPDG(static_cast<PDGCode>(event[i].id()));
- HEPEnergyType const Erest = event[i].e() * 1_GeV;
- MomentumVector const pRest(
- rotatedCS,
- {event[i].px() * 1_GeV, event[i].py() * 1_GeV, event[i].pz() * 1_GeV});
- geometry::FourVector const fourMomRest{Erest, pRest};
- auto const fourMomLab = boost.fromCoM(fourMomRest);
-
- cout << "particle: id=" << pyId << " momentum="
- << fourMomLab.GetSpaceLikeComponents().GetComponents(labCS) / 1_GeV
- << " energy=" << fourMomLab.GetTimeLikeComponent() << endl;
-
- view.AddSecondary(
- tuple<Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
- geometry::Point, units::si::TimeType>{
- pyId, fourMomLab.GetTimeLikeComponent(),
- fourMomLab.GetSpaceLikeComponents(), decayPoint, t0});
- }
-
- // set particle stable
- Decay::SetStable(particleId);
- }
-
-} // namespace corsika::pythia
diff --git a/dependencies/Pythia/Decay.h b/dependencies/Pythia/Decay.h
deleted file mode 100644
index 48ce89f8e47c2ae5717fcfd60514d0959a8087ad..0000000000000000000000000000000000000000
--- a/dependencies/Pythia/Decay.h
+++ /dev/null
@@ -1,77 +0,0 @@
-/*
- * (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 <Pythia8/Pythia.h>
-#include <corsika/framework/coreParticleProperties.h>
-#include <corsika/process/DecayProcess.h>
-
-namespace corsika {
-
- namespace pythia {
-
- typedef corsika::Vector<corsika::units::si::hepmomentum_d> MomentumVector;
-
- class Decay : public corsika::DecayProcess<Decay> {
- int fCount = 0;
- bool handleAllDecays_ = true;
-
- public:
- Decay();
- Decay(std::set<Code>);
- ~Decay();
-
- // is Pythia::Decay set to handle the decay of this particle?
- bool IsDecayHandled(const corsika::Code);
-
- // is decay possible in principle?
- bool CanHandleDecay(const corsika::Code);
-
- // set Pythia::Decay to handle the decay of this particle!
- void SetHandleDecay(const corsika::Code);
- // set Pythia::Decay to handle the decay of this list of particles!
- void SetHandleDecay(const std::vector<Code>);
- // set Pythia::Decay to handle all particle decays
- void SetHandleAllDecays();
-
- // print internal configuration for this particle
- void PrintDecayConfig(const corsika::Code);
- // print configuration of decays in corsika
- void PrintDecayConfig();
-
- bool CanDecay(const corsika::Code);
-
- /**
- In this function PYTHIA is asked for the lifetime of the input particle.
- Unknown particles should return an infinite lifetime so that another decay process
- can act on the particle later on.
- */
-
- template <typename TParticle>
- corsika::units::si::TimeType GetLifetime(TParticle const&);
-
- /**
- In this function PYTHIA is called to execute the decay of the input particle.
- */
-
- template <typename TSecondaryView>
- void DoDecay(TSecondaryView&);
-
- private:
- void SetUnstable(const corsika::Code);
- void SetStable(const corsika::Code);
- void SetStable(const std::vector<Code>);
- bool IsStable(const corsika::Code);
-
- Pythia8::Pythia fPythia;
- std::set<Code> handledDecays_;
- };
-
- } // namespace pythia
-} // namespace corsika
diff --git a/dependencies/Pythia/testPythia8.cc b/dependencies/Pythia/testPythia8.cc
deleted file mode 100644
index f495f56ecf0472ba35e7145976b5c2826fec02c4..0000000000000000000000000000000000000000
--- a/dependencies/Pythia/testPythia8.cc
+++ /dev/null
@@ -1,185 +0,0 @@
-/*
- * (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 <Pythia8/Pythia.h>
-#include <dependencies/Pythia/Decay.h>
-#include <dependencies/Pythia/Interaction.h>
-
-#include <corsika/framework/random/RNGManager.hpp>
-
-#include <corsika/framework/coreParticleProperties.h>
-
-#include <corsika/framework/geometry/Point.h>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-
-#include <corsika/utl/CorsikaFenv.h>
-#include <catch2/catch.hpp>
-
-TEST_CASE("Pythia", "[processes]") {
-
- SECTION("linking pythia") {
- using namespace Pythia8;
- using std::cout;
- using std::endl;
-
- // Generator. Process selection. LHC initialization. Histogram.
- Pythia pythia;
-
- pythia.readString("Next:numberShowInfo = 0");
- pythia.readString("Next:numberShowProcess = 0");
- pythia.readString("Next:numberShowEvent = 0");
-
- pythia.readString("ProcessLevel:all = off");
-
- pythia.init();
-
- Event& event = pythia.event;
- event.reset();
-
- pythia.particleData.mayDecay(321, true);
- double pz = 100.;
- double m = 0.49368;
- event.append(321, 1, 0, 0, 0., 0., 100., sqrt(pz * pz + m * m), m);
-
- if (!pythia.next())
- cout << "decay failed!" << endl;
- else
- cout << "particles after decay: " << event.size() << endl;
- event.list();
-
- // loop over final state
- for (int i = 0; i < pythia.event.size(); ++i)
- if (pythia.event[i].isFinal()) {
- cout << "particle: id=" << pythia.event[i].id() << endl;
- }
- }
-
- SECTION("pythia interface") {
- using namespace corsika;
-
- random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-
- process::pythia::Decay model;
- }
-}
-
-#include <corsika/framework/geometry/Point.h>
-#include <corsika/framework/geometry/RootCoordinateSystem.h>
-#include <corsika/framework/geometry/Vector.h>
-
-#include <corsika/framework/core/PhysicalUnits.hpp>
-
-#include <corsika/framework/coreParticleProperties.h>
-#include <corsika/setup/SetupStack.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.h>
-
-using namespace corsika;
-using namespace corsika::units::si;
-
-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("pythia process") {
-
- // setup environment, geometry
- environment::Environment<environment::IMediumModel> env;
-
- geometry::CoordinateSystem const& cs = env.GetCoordinateSystem();
-
- auto theMedium =
- environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
- geometry::Point{cs, 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<Code>{Code::Hydrogen},
- std::vector<float>{1.}));
-
- auto const* nodePtr = theMedium.get(); // save the medium for later use before moving it
-
- SECTION("pythia decay") {
- feenableexcept(FE_INVALID);
- setup::Stack stack;
- const HEPEnergyType E0 = 10_GeV;
- HEPMomentumType P0 = sqrt(E0 * E0 - PiPlus::GetMass() * PiPlus::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<Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
- geometry::Point, units::si::TimeType>{Code::PiPlus, E0, plab, pos,
- 0_ns});
-
- random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-
- setup::StackView view(particle);
-
- process::pythia::Decay model;
-
- [[maybe_unused]] const TimeType time = model.GetLifetime(particle);
- model.DoDecay(view);
- CHECK(stack.getEntries() == 3);
- auto const pSum = sumMomentum(view, cs);
- CHECK((pSum - plab).norm() / 1_GeV == Approx(0).margin(1e-4));
- CHECK((pSum.norm() - plab.norm()) / 1_GeV == Approx(0).margin(1e-4));
- }
-
- SECTION("pythia decay config") {
- process::pythia::Decay model({Code::PiPlus, Code::PiMinus});
- REQUIRE(model.IsDecayHandled(Code::PiPlus));
- REQUIRE(model.IsDecayHandled(Code::PiMinus));
- REQUIRE_FALSE(model.IsDecayHandled(Code::KPlus));
-
- const std::vector<Code> particleTestList = {Code::PiPlus, Code::PiMinus, Code::KPlus,
- Code::Lambda0Bar, Code::D0Bar};
-
- // setup decays
- model.SetHandleDecay(particleTestList);
- for (auto& pCode : particleTestList) REQUIRE(model.IsDecayHandled(pCode));
-
- // individually
- model.SetHandleDecay(Code::KMinus);
-
- // possible decays
- REQUIRE_FALSE(model.CanHandleDecay(Code::Proton));
- REQUIRE_FALSE(model.CanHandleDecay(Code::Electron));
- REQUIRE(model.CanHandleDecay(Code::PiPlus));
- REQUIRE(model.CanHandleDecay(Code::MuPlus));
- }
-
- SECTION("pythia interaction") {
-
- setup::Stack stack;
- const HEPEnergyType E0 = 100_GeV;
- HEPMomentumType P0 = sqrt(E0 * E0 - PiPlus::GetMass() * PiPlus::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<Code, units::si::HEPEnergyType, corsika::stack::MomentumVector,
- geometry::Point, units::si::TimeType>{Code::PiPlus, E0, plab, pos,
- 0_ns});
- particle.SetNode(nodePtr);
- setup::StackView view(particle);
-
- process::pythia::Interaction model;
-
- [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(view);
- [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
- }
-}
diff --git a/dependencies/pythia/CMakeLists.txt b/dependencies/pythia/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..701303282c3de35ae9515cdb3fb26dd4b333516f
--- /dev/null
+++ b/dependencies/pythia/CMakeLists.txt
@@ -0,0 +1,75 @@
+include(ExternalProject)
+# eventually add AUTO here, too:
+set (ThirdPartyChoiceValues "C8;SYSTEM" CACHE STRING
+ "List of possible values for the ThirdParty package choice")
+mark_as_advanced (ThirdPartyChoiceValues)
+
+##############################################################################
+# check for Pythia8: either use C8 or system-level installation
+
+message ("***** Configuring Pythia8 version")
+
+set (USE_PYTHIA8_C8 "C8" CACHE STRING
+ "Selection of pythia8 package. Can be \'C8\' or \'SYSTEM\'. Default: \'C8\'.")
+set_property (CACHE USE_PYTHIA8_C8 PROPERTY STRINGS ${ThirdPartyChoiceValues})
+if (NOT (${USE_PYTHIA8_C8} IN_LIST ThirdPartyChoiceValues))
+ message (SEND_ERROR "Illegal USE_PYTHIA8_C8=\"${USE_PYTHIA8_C8}\" can only be one of: ${ThirdPartyChoiceValues}")
+endif (NOT (${USE_PYTHIA8_C8} IN_LIST ThirdPartyChoiceValues))
+message (STATUS "USE_PYTHIA8_C8='${USE_PYTHIA8_C8}'")
+
+add_library (C8::ext::pythia8 STATIC IMPORTED GLOBAL)
+if ("x_${USE_PYTHIA8_C8}" STREQUAL "x_SYSTEM")
+
+ find_package (Pythia8 REQUIRED)
+ message (STATUS "Using system-level Pythia8 version ${Pythia8_VERSION} at ${Pythia8_INCLUDE_DIR}")
+ set_target_properties (
+ C8::ext::pythia8 PROPERTIES
+ IMPORTED_LOCATION ${Pythia8_LIBRARY}
+ IMPORTED_LINK_INTERFACE_LIBRARIES dl
+ INTERFACE_INCLUDE_DIRECTORIES ${Pythia8_INCLUDE_DIR}
+ )
+ set (Pythia8_FOUND 1 PARENT_SCOPE)
+
+else ()
+
+ set (_C8_Pythia8_VERSION "8235")
+ message (STATUS "Building ThirdParty/pythia8 using pythia${_C8_Pythia8_VERSION}-stripped.tar.bz2")
+ message (STATUS "This will take a bit.....")
+ ExternalProject_Add (pythia8
+ URL ${CMAKE_CURRENT_SOURCE_DIR}/pythia${_C8_Pythia8_VERSION}-stripped.tar.bz2
+ URL_MD5 83132880c0594b808bd7fd37fb642606
+ SOURCE_DIR ${CMAKE_CURRENT_BINARY_DIR}/pythia8/source
+ INSTALL_DIR ${CMAKE_CURRENT_BINARY_DIR}/pythia8/install
+ CONFIGURE_COMMAND ./configure --cxx-common=-Wno-deprecated-copy --prefix=${CMAKE_CURRENT_BINARY_DIR}/pythia8/install
+ BUILD_IN_SOURCE ON
+ EXCLUDE_FROM_ALL TRUE
+ )
+ set (Pythia8_FOUND 1 PARENT_SCOPE)
+ ExternalProject_Get_Property (pythia8 INSTALL_DIR)
+ set (Pythia8_VERSION ${_C8_Pythia8_VERSION} CACHE STRING "Version of Pythia8")
+ set (Pythia8_PREFIX ${INSTALL_DIR})
+ set (Pythia8_INCLUDE_DIR ${Pythia8_PREFIX}/include)
+ set (Pythia8_LIBRARY ${Pythia8_PREFIX}/lib/libpythia8.a)
+ add_dependencies (C8::ext::pythia8 pythia8)
+
+ # create include directory at config time
+ file (MAKE_DIRECTORY ${Pythia8_INCLUDE_DIR})
+
+ set (Pythia8_INSTALL_DIR ${CMAKE_INSTALL_PREFIX}/share/externals/pythia8)
+ install (DIRECTORY ${INSTALL_DIR}/ DESTINATION ${Pythia8_INSTALL_DIR})
+
+ set_target_properties (
+ C8::ext::pythia8 PROPERTIES
+ IMPORTED_LOCATION ${Pythia8_LIBRARY}
+ IMPORTED_LINK_INTERFACE_LIBRARIES dl
+ INTERFACE_INCLUDE_DIRECTORIES
+ $<BUILD_INTERFACE:${Pythia8_INCLUDE_DIR}>
+ )
+
+ message (STATUS "Pythia8_INCLUDE_DIR ${Pythia8_INCLUDE_DIR}")
+endif ()
+
+
+##### add pythia8 to CORSIKA8 build
+add_dependencies (CORSIKA8 C8::ext::pythia8)
+target_link_libraries (CORSIKA8 INTERFACE C8::ext::pythia8)
diff --git a/dependencies/pythia/pythia8235-stripped.tar.bz2 b/dependencies/pythia/pythia8235-stripped.tar.bz2
new file mode 100644
index 0000000000000000000000000000000000000000..da280f738809e54eace85622485d13492448250e
Binary files /dev/null and b/dependencies/pythia/pythia8235-stripped.tar.bz2 differ
diff --git a/dependencies/sibyll/CMakeLists.txt b/dependencies/sibyll/CMakeLists.txt
index 32c967165b1e2cd6fa924348cbf9f699f40ac0de..86c9f9791a679644212deb69cb58cb9db46613ac 100644
--- a/dependencies/sibyll/CMakeLists.txt
+++ b/dependencies/sibyll/CMakeLists.txt
@@ -44,4 +44,6 @@ target_link_libraries (
gfortran
)
+# add sibyll to corsika8 build
add_dependencies (CORSIKA8 Sibyll_static)
+target_link_libraries (CORSIKA8 INTERFACE Sibyll_static)
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 7aa9031a498a25bff9be6169e66969dbc0b8ac95..8572b734653e4f9cf16ec6ffa64c5de21e6f3700 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -5,9 +5,8 @@ set (test_modules_sources
testNullModel.cpp
testObservationPlane.cpp
testParticleCut.cpp
- #testPythia.cpp
+ testPythia8.cpp
#testQGSJetII.cpp
- #testSibyll.cpp
#testStackInspector.cpp
#testSwitchProcess.cpp
#testTrackingLine.cpp
@@ -16,7 +15,7 @@ set (test_modules_sources
add_executable (testModules ${test_modules_sources})
-target_link_libraries (testModules CORSIKA8 Sibyll_static Catch2)
+target_link_libraries (testModules CORSIKA8 Catch2)
target_compile_options (testModules PRIVATE -g) # do not skip asserts
diff --git a/tests/modules/testPythia.cpp b/tests/modules/testPythia.cpp
deleted file mode 100644
index f229b907f92968d7c26ff7bed70b4476ce56282d..0000000000000000000000000000000000000000
--- a/tests/modules/testPythia.cpp
+++ /dev/null
@@ -1,165 +0,0 @@
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <Pythia8/Pythia.h>
-#include <corsika/process/pythia/Decay.h>
-#include <corsika/process/pythia/Interaction.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>
-
-TEST_CASE("Pythia", "[processes]") {
-
- SECTION("linking pythia") {
- using namespace Pythia8;
- using std::cout;
- using std::endl;
-
- // Generator. Process selection. LHC initialization. Histogram.
- Pythia pythia;
-
- pythia.readString("Next:numberShowInfo = 0");
- pythia.readString("Next:numberShowProcess = 0");
- pythia.readString("Next:numberShowEvent = 0");
-
- pythia.readString("ProcessLevel:all = off");
-
- pythia.init();
-
- Event& event = pythia.event;
- event.reset();
-
- pythia.particleData.mayDecay(321, true);
- double pz = 100.;
- double m = 0.49368;
- event.append(321, 1, 0, 0, 0., 0., 100., sqrt(pz * pz + m * m), m);
-
- if (!pythia.next())
- cout << "decay failed!" << endl;
- else
- cout << "particles after decay: " << event.size() << endl;
- event.list();
-
- // loop over final state
- for (int i = 0; i < pythia.event.size(); ++i)
- if (pythia.event[i].isFinal()) {
- cout << "particle: id=" << pythia.event[i].id() << endl;
- }
- }
-
- SECTION("pythia interface") {
- using namespace corsika;
-
- const std::vector<particles::Code> particleList = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
-
- random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-
- process::pythia::Decay model(particleList);
-
- model.Init();
- }
-}
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-
-#include <corsika/units/PhysicalUnits.h>
-
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/setup/SetupStack.h>
-#include <corsika/setup/SetupTrajectory.h>
-
-#include <corsika/environment/Environment.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.h>
-
-using namespace corsika;
-using namespace corsika::units::si;
-
-TEST_CASE("pythia process") {
-
- // setup environment, geometry
- environment::Environment<environment::IMediumModel> env;
-
- geometry::CoordinateSystem const& cs = env.GetCoordinateSystem();
-
- auto theMedium =
- environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
- geometry::Point{cs, 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::Hydrogen},
- std::vector<float>{1.}));
-
- auto const* nodePtr = theMedium.get(); // save the medium for later use before moving it
-
- SECTION("pythia decay") {
-
- setup::Stack stack;
- const HEPEnergyType E0 = 10_GeV;
- HEPMomentumType P0 =
- sqrt(E0 * E0 - particles::PiPlus::GetMass() * particles::PiPlus::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::PiPlus, E0, plab, pos, 0_ns});
-
- const std::vector<particles::Code> particleList = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
-
- random::RNGManager::GetInstance().RegisterRandomStream("pythia");
-
- corsika::stack::SecondaryView view(particle);
- auto projectile = view.GetProjectile();
-
- process::pythia::Decay model(particleList);
- model.Init();
- model.DoDecay(projectile);
- [[maybe_unused]] const TimeType time = model.GetLifetime(particle);
- }
-
- SECTION("pythia interaction") {
-
- setup::Stack stack;
- const HEPEnergyType E0 = 100_GeV;
- HEPMomentumType P0 =
- sqrt(E0 * E0 - particles::PiPlus::GetMass() * particles::PiPlus::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::PiPlus, E0, plab, pos, 0_ns});
- particle.SetNode(nodePtr);
- corsika::stack::SecondaryView view(particle);
- auto projectile = view.GetProjectile();
-
- process::pythia::Interaction model;
- model.Init();
- [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
- [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
- }
-}
diff --git a/tests/modules/testPythia8.cpp b/tests/modules/testPythia8.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..07cc2d930cb0378e81c352e8c1d22b5471f63ffd
--- /dev/null
+++ b/tests/modules/testPythia8.cpp
@@ -0,0 +1,163 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/modules/pythia8/Decay.hpp>
+#include <corsika/modules/pythia8/Interaction.hpp>
+
+#include <corsika/framework/random/RNGManager.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/framework/geometry/Point.hpp>
+
+#include <catch2/catch.hpp>
+
+TEST_CASE("Pythia", "[processes]") {
+
+ SECTION("linking pythia") {
+ using namespace Pythia8;
+ using std::cout;
+ using std::endl;
+
+ // Generator. Process selection. LHC initialization. Histogram.
+ Pythia pythia;
+
+ pythia.readString("Next:numberShowInfo = 0");
+ pythia.readString("Next:numberShowProcess = 0");
+ pythia.readString("Next:numberShowEvent = 0");
+
+ pythia.readString("ProcessLevel:all = off");
+
+ pythia.init();
+
+ Event& event = pythia.event;
+ event.reset();
+
+ pythia.particleData.mayDecay(321, true);
+ double pz = 100.;
+ double m = 0.49368;
+ event.append(321, 1, 0, 0, 0., 0., 100., sqrt(pz * pz + m * m), m);
+
+ if (!pythia.next())
+ cout << "decay failed!" << endl;
+ else
+ cout << "particles after decay: " << event.size() << endl;
+ event.list();
+
+ // loop over final state
+ for (int i = 0; i < pythia.event.size(); ++i)
+ if (pythia.event[i].isFinal()) {
+ cout << "particle: id=" << pythia.event[i].id() << endl;
+ }
+ }
+
+ SECTION("pythia interface") {
+ using namespace corsika;
+
+ const std::vector<corsika::Code> particleList = {
+ corsika::Code::PiPlus, corsika::Code::PiMinus, corsika::Code::KPlus,
+ corsika::Code::KMinus, corsika::Code::K0Long, corsika::Code::K0Short};
+
+ corsika::RNGManager::GetInstance().RegisterRandomStream("pythia");
+
+ corsika::pythia8::Decay model(particleList);
+
+ model.Init();
+ }
+}
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+
+using namespace corsika;
+using namespace corsika::units::si;
+
+TEST_CASE("pythia process") {
+
+ // setup environment, geometry
+ corsika::Environment<corsika::IMediumModel> env;
+
+ corsika::CoordinateSystem const& cs = env.GetCoordinateSystem();
+
+ auto theMedium =
+ corsika::Environment<corsika::IMediumModel>::CreateNode<corsika::Sphere>(
+ corsika::Point{cs, 0_m, 0_m, 0_m},
+ 1_km * std::numeric_limits<double>::infinity());
+
+ using MyHomogeneousModel = corsika::HomogeneousMedium<corsika::IMediumModel>;
+ theMedium->SetModelProperties<MyHomogeneousModel>(
+ 1_kg / (1_m * 1_m * 1_m),
+ corsika::NuclearComposition(std::vector<corsika::Code>{corsika::Code::Hydrogen},
+ std::vector<float>{1.}));
+
+ auto const* nodePtr = theMedium.get(); // save the medium for later use before moving it
+
+ SECTION("pythia decay") {
+
+ setup::Stack stack;
+ const HEPEnergyType E0 = 10_GeV;
+ HEPMomentumType P0 =
+ sqrt(E0 * E0 - corsika::PiPlus::GetMass() * corsika::PiPlus::GetMass());
+ auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+ corsika::Point pos(cs, 0_m, 0_m, 0_m);
+ auto particle = stack.AddParticle(
+ std::tuple<corsika::Code, units::si::HEPEnergyType, corsika::MomentumVector,
+ corsika::Point, units::si::TimeType>{corsika::Code::PiPlus, E0, plab,
+ pos, 0_ns});
+
+ const std::vector<corsika::Code> particleList = {
+ corsika::Code::PiPlus, corsika::Code::PiMinus, corsika::Code::KPlus,
+ corsika::Code::KMinus, corsika::Code::K0Long, corsika::Code::K0Short};
+
+ corsika::RNGManager::GetInstance().RegisterRandomStream("pythia");
+
+ corsika::SecondaryView view(particle);
+ auto projectile = view.GetProjectile();
+
+ corsika::pythia8::Decay model(particleList);
+ model.Init();
+ model.DoDecay(projectile);
+ [[maybe_unused]] const TimeType time = model.GetLifetime(particle);
+ }
+
+ SECTION("pythia interaction") {
+
+ setup::Stack stack;
+ const HEPEnergyType E0 = 100_GeV;
+ HEPMomentumType P0 =
+ sqrt(E0 * E0 - corsika::PiPlus::GetMass() * corsika::PiPlus::GetMass());
+ auto plab = corsika::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+ corsika::Point pos(cs, 0_m, 0_m, 0_m);
+ auto particle = stack.AddParticle(
+ std::tuple<corsika::Code, units::si::HEPEnergyType, corsika::MomentumVector,
+ corsika::Point, units::si::TimeType>{corsika::Code::PiPlus, E0, plab,
+ pos, 0_ns});
+ particle.SetNode(nodePtr);
+ corsika::SecondaryView view(particle);
+ auto projectile = view.GetProjectile();
+
+ corsika::pythia8::Interaction model;
+ model.Init();
+ [[maybe_unused]] const corsika::EProcessReturn ret = model.DoInteraction(projectile);
+ [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
+ }
+}