From 8b729687ef16caeadc2b26f501120f6a8f7ea7e7 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Thu, 9 Feb 2023 13:51:41 +0100
Subject: [PATCH] lots of things corrected, cross-sections working now
---
.../detail/modules/fluka/InteractionModel.inl | 176 ++++++++++-------
corsika/modules/FLUKA.hpp | 31 +++
corsika/modules/fluka/InteractionModel.hpp | 53 ++---
.../modules}/fluka/ParticleConversion.hpp | 0
modules/fluka/FLUKA.hpp | 185 ++++++++++++------
src/modules/fluka/fluka_codes.dat | 153 +++++----------
tests/modules/testFluka.cpp | 85 ++++++--
7 files changed, 399 insertions(+), 284 deletions(-)
create mode 100644 corsika/modules/FLUKA.hpp
rename {modules => corsika/modules}/fluka/ParticleConversion.hpp (100%)
diff --git a/corsika/detail/modules/fluka/InteractionModel.inl b/corsika/detail/modules/fluka/InteractionModel.inl
index 197adb54b..530cafbe0 100644
--- a/corsika/detail/modules/fluka/InteractionModel.inl
+++ b/corsika/detail/modules/fluka/InteractionModel.inl
@@ -5,10 +5,11 @@
* Licence version 3 (GPL Version 3). See file LICENSE for a full version of
* the license.
*/
-
+
#pragma once
#include <algorithm>
+#include <cstdlib>
#include <vector>
#include <iterator>
#include <set>
@@ -23,66 +24,94 @@
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <FLUKA.hpp>
-#include <ParticleConversion.hpp>
+#include <corsika/modules/fluka/ParticleConversion.hpp>
namespace corsika::fluka {
template <typename TEnvironment>
- inline InteractionModel::InteractionModel(TEnvironment const& env) :materials_{genFlukaMaterials(env)} {
- for (auto const& [code, matno]: materials_) {
- std::cout << code << " " << matno << std::endl;
- }
+ inline InteractionModel::InteractionModel(TEnvironment const& env)
+ : materials_{genFlukaMaterials(env)} {
+ for (auto const& [code, matno] : materials_) {
+ std::cout << get_name(code, full_name{}) << " " << matno << std::endl;
+ }
+ }
+
+ inline bool InteractionModel::isValid(Code projectileID, int material,
+ HEPEnergyType /*sqrtS*/) const {
+ if (!fluka::canInteract(projectileID)) {
+ // invalid projectile
+ //~ std::cout << "invalid projecile (cannot interact)";
+ return false;
}
-
- inline bool InteractionModel::isValid(Code projectileID, Code targetID, HEPEnergyType /*sqrtS*/) const {
- if (fluka::canInteract(projectileID)) {
- // invalid projectile
- return false;
- }
-
- if (getMaterialIndex(targetID) == -1) {
- // unknown material
- return false;
- }
-
- // TODO: check validity range of sqrtS
- return true;
+
+ if (material < 0) {
+ // invalid/uninitialized target
+ return false;
+ }
+
+ // TODO: check validity range of sqrtS
+ return true;
+ }
+
+ inline bool InteractionModel::isValid(Code projectileID, Code targetID,
+ HEPEnergyType sqrtS) const {
+ return isValid(projectileID, getMaterialIndex(targetID), sqrtS);
}
-
+
inline int InteractionModel::getMaterialIndex(Code targetID) const {
- if (auto it = std::find(materials_.cbegin(), materials_.cend(), [=](std::pair<Code, int> const& p) {return p.first == targetID;});
+ auto compare = [=](std::pair<Code, int> const& p) { return p.first == targetID; };
+ if (auto it = std::find_if(materials_.cbegin(), materials_.cend(), compare);
it == materials_.cend()) {
- return -1;
+ return -1;
} else {
- return it->second;
+ return it->second;
}
}
-
+
inline CrossSectionType InteractionModel::getCrossSection(
Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
FourMomentum const& targetP4) const {
+ auto const flukaCodeProj =
+ static_cast<FLUKACodeIntType>(convertToFluka(projectileId));
+ auto const flukaMaterial = getMaterialIndex(targetId);
+
HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
- if (!isValid(projectileId, targetId, sqrtS)) { return CrossSectionType::zero(); }
-
- COMBoost const targetRestBoost{projectileP4.getSpaceLikeComponents(), get_mass(targetId)};
+ if (!isValid(projectileId, flukaMaterial, sqrtS)) { return CrossSectionType::zero(); }
+
+ COMBoost const targetRestBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
FourMomentum const projectileLab4mom = targetRestBoost.toCoM(projectileP4);
- //~ projec
-
- int const iflxyz = 1;
-
- // SIGREA = SGMXYZ ( KPROJ, MMAT, EKIN, PPROJ, iflxyz );
+ HEPEnergyType const Elab = projectileLab4mom.getTimeLikeComponent();
+ auto constexpr invGeV = 1 / 1_GeV;
+ double const EkinLab = (Elab - get_mass(projectileId)) * invGeV;
+ double const labMomentum =
+ projectileLab4mom.getSpaceLikeComponents().getNorm() * invGeV;
+
+ auto const plab = projectileLab4mom.getSpaceLikeComponents();
+
+ std::cout << targetRestBoost.boost_ << '\n';
+ std::cout << targetRestBoost.rotatedCS_->getTransform().matrix() << '\n';
+ std::cout << "Elab / GeV = " << Elab * invGeV << '\n';
+ std::cout << "plab = " << plab << '\n';
+ std::cout << "EkinLab / GeV = " << EkinLab << '\n';
+ std::cout << "labMomentum / GeV = "
+ << projectileLab4mom.getSpaceLikeComponents() * invGeV << '\n';
+
+ CrossSectionType const xs = ::fluka::sgmxyz_(&flukaCodeProj, &flukaMaterial, &EkinLab,
+ &labMomentum, &iflxyz_) *
+ 1_mb;
+ return xs;
}
-
- template <typename TSecondaryView>
- inline void InteractionModel::doInteraction(
- TSecondaryView& view, Code const projectileId, Code const targetId,
- FourMomentum const& projectileP4, FourMomentum const& targetP4) {
-
- }
-
-
+
+ template <typename TSecondaryView>
+ inline void InteractionModel::doInteraction(TSecondaryView& view,
+ Code const projectileId,
+ Code const targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) {}
+
template <typename TEnvironment>
- inline std::vector<std::pair<Code, int>> InteractionModel::genFlukaMaterials(TEnvironment const& env) {
- auto const& universe = *(env.getUniverse());
+ inline std::vector<std::pair<Code, int>> InteractionModel::genFlukaMaterials(
+ TEnvironment const& env) {
+ auto const& universe = *(env.getUniverse());
// generate complete list of all nuclei types in universe
auto const allElementsInUniverse = std::invoke([&]() {
@@ -97,26 +126,26 @@ namespace corsika::fluka {
universe.walk(collectElements);
return allElementsInUniverse;
});
-
- /*
- * We define one material per element/isotope we have in C8. Cross-section averaging and
- * target isotope selection happen in C8.
+
+ /*
+ * We define one material per element/isotope we have in C8. Cross-section averaging
+ * and target isotope selection happen in C8.
*/
-
+
int const nElements = allElementsInUniverse.size();
auto nelmfl = std::make_unique<int[]>(nElements);
std::vector<int> izelfl;
izelfl.reserve(nElements);
auto wfelml = std::make_unique<double[]>(nElements);
- auto const mxelfl = nElements;
+ auto const mxelfl = nElements;
double const pptmax = 1e11; // GeV
- double const ef2dp3 = 0; // GeV, 0 means default is used
- double const df2dp3 = -1; // default
+ double const ef2dp3 = 0; // GeV, 0 means default is used
+ double const df2dp3 = -1; // default
bool const lprint = true;
auto mtflka = std::make_unique<int[]>(mxelfl);
- char crvrck[8+1] = "76466879"; // magic number that FLUKA uses to see if it's the right version
-
-
+ char crvrck[8 + 1] =
+ "76466879"; // magic number that FLUKA uses to see if it's the right version
+
/*
* Iflxyz = 1 -> only inelastic
* Iflxyz = 10 -> only elastic
@@ -126,31 +155,34 @@ namespace corsika::fluka {
* Iflxyz =110 -> elastic + emd
* Iflxyz =111 -> inelastic + elastic + emd
*/
- int const iflxyz_ = 1;
-
+
std::fill(&nelmfl[0], &nelmfl[nElements], 1);
std::fill(&wfelml[0], &wfelml[nElements], 1.);
std::transform(allElementsInUniverse.cbegin(), allElementsInUniverse.cend(),
- std::back_inserter(izelfl), get_nucleus_Z);
-
+ std::back_inserter(izelfl), get_nucleus_Z);
+
+ // check if FLUPRO is set
+ if (std::getenv("FLUPRO") == nullptr) {
+ throw std::runtime_error{"FLUPRO not set; required to initialize FLUKA"};
+ }
+
// call FLUKA
::fluka::stpxyz_(&nElements, nelmfl.get(), izelfl.data(), wfelml.get(), &mxelfl,
- &pptmax, &ef2dp3, &df2dp3, &iflxyz_, &lprint, mtflka.get(), crvrck);
-
+ &pptmax, &ef2dp3, &df2dp3, &iflxyz_, &lprint, mtflka.get(), crvrck);
+
// now create & fill vector of (C8 Code, FLUKA mat. no.) pairs
std::vector<std::pair<Code, int>> mapping;
mapping.reserve(nElements);
-
- auto it = boost::make_zip_iterator(boost::make_tuple(
- allElementsInUniverse.begin(), &mtflka[0]));
- auto end = boost::make_zip_iterator(boost::make_tuple(
- allElementsInUniverse.end(), &mtflka[nElements]));
+
+ auto it = boost::make_zip_iterator(
+ boost::make_tuple(allElementsInUniverse.begin(), &mtflka[0]));
+ auto const end = boost::make_zip_iterator(
+ boost::make_tuple(allElementsInUniverse.end(), &mtflka[nElements]));
for (; it != end; ++it) {
- boost::tuple<Code const&, int&> tup = *it;
- mapping.emplace_back(tup.get<0>(), tup.get<1>());
+ boost::tuple<Code const&, int&> tup = *it;
+ mapping.emplace_back(tup.get<0>(), tup.get<1>());
}
-
- return mapping;
- }
-}
+ return mapping;
+ }
+} // namespace corsika::fluka
diff --git a/corsika/modules/FLUKA.hpp b/corsika/modules/FLUKA.hpp
new file mode 100644
index 000000000..02555ea5e
--- /dev/null
+++ b/corsika/modules/FLUKA.hpp
@@ -0,0 +1,31 @@
+/*
+ * (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/modules/fluka/ParticleConversion.hpp>
+#include <corsika/modules/fluka/InteractionModel.hpp>
+#include <corsika/framework/process/InteractionProcess.hpp>
+
+/**
+ * @file Epos.hpp
+ *
+ * Includes all the parts of the EPOS model. Defines the InteractionProcess<TModel>
+ * classes needed for the ProcessSequence.
+ */
+
+namespace corsika::fluka {
+ /**
+ * fluka::Interaction is the process for ProcessSequence.
+ *
+ * The fluka::InteractionModel is wrapped as an InteractionProcess here in order
+ * to provide all the functions for ProcessSequence.
+ */
+ class Interaction : public fluka::InteractionModel,
+ public corsika::InteractionProcess<Interaction> {};
+} // namespace corsika::fluka
diff --git a/corsika/modules/fluka/InteractionModel.hpp b/corsika/modules/fluka/InteractionModel.hpp
index 836ec2105..d2ecf12f9 100644
--- a/corsika/modules/fluka/InteractionModel.hpp
+++ b/corsika/modules/fluka/InteractionModel.hpp
@@ -8,6 +8,9 @@
#pragma once
+#include <vector>
+#include <utility>
+
#include <corsika/media/Environment.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
@@ -16,31 +19,33 @@
namespace corsika::fluka {
class InteractionModel {
- public:
- template <typename TEnvironment>
- InteractionModel(TEnvironment const&);
-
- CrossSectionType getCrossSection(
- Code projectileId, Code targetId, FourMomentum const& projectileP4,
- FourMomentum const& targetP4) const;
-
- bool isValid(Code projectileID, Code targetID, HEPEnergyType sqrtS) const;
-
- int getMaterialIndex(Code targetID) const;
-
- template <typename TSecondaryView>
- void doInteraction(
- TSecondaryView& view, Code const projectileId, Code const targetId,
- FourMomentum const& projectileP4, FourMomentum const& targetP4);
-
- //~ static int const iflxyz_ = 1; //!< select interaction types (see fluka.h); hardcoded to inel. for now
- private:
- std::vector<std::pair<Code, int>> const materials_; //!< map target Code to FLUKA material no.
-
- template <typename TEnvironment>
- static std::vector<std::pair<Code, int>> genFlukaMaterials(TEnvironment const&);
+ public:
+ template <typename TEnvironment>
+ InteractionModel(TEnvironment const&);
+
+ CrossSectionType getCrossSection(Code projectileId, Code targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const;
+
+ bool isValid(Code projectileID, Code targetID, HEPEnergyType sqrtS) const;
+ bool isValid(Code projectileID, int material, HEPEnergyType sqrtS) const;
+
+ int getMaterialIndex(Code targetID) const;
+
+ template <typename TSecondaryView>
+ void doInteraction(TSecondaryView& view, Code const projectileId, Code const targetId,
+ FourMomentum const& projectileP4, FourMomentum const& targetP4);
+
+ private:
+ std::vector<std::pair<Code, int>> const
+ materials_; //!< map target Code to FLUKA material no.
+
+ template <typename TEnvironment>
+ static std::vector<std::pair<Code, int>> genFlukaMaterials(TEnvironment const&);
// TODO: random number stream
};
-}
+
+ inline static int const iflxyz_ = 1;
+} // namespace corsika::fluka
#include <corsika/detail/modules/fluka/InteractionModel.inl>
diff --git a/modules/fluka/ParticleConversion.hpp b/corsika/modules/fluka/ParticleConversion.hpp
similarity index 100%
rename from modules/fluka/ParticleConversion.hpp
rename to corsika/modules/fluka/ParticleConversion.hpp
diff --git a/modules/fluka/FLUKA.hpp b/modules/fluka/FLUKA.hpp
index 40457737a..76c542b03 100644
--- a/modules/fluka/FLUKA.hpp
+++ b/modules/fluka/FLUKA.hpp
@@ -9,68 +9,123 @@
#pragma once
namespace fluka {
-extern "C" {
-/*----------------------------------------------------------------------*
-* *
-* Copyright (C) 2022-2022 by Alfredo Ferrari & Paola Sala *
-* All Rights Reserved. *
-* *
-* *
-* SeTuP for X/Y/Z interaction types: *
-* *
-* Authors: Alfredo Ferrari & Paola Sala *
-* *
-* *
-* Created on 24 October 2022 by Alfredo Ferrari & Paola Sala *
-* Private Private *
-* *
-* Last change on 01-Nov-22 by Alfredo Ferrari *
-* Private *
-* *
-* Input variables: *
-* *
-* Nmatfl = Number of requested Fluka materials *
-* Nelmfl(m) = Number of elements in the m_th requested Fluka *
-* material, Nelmfl(m)=1 means simple element, no *
-* compound *
-* Izelfl(n) = Cumulative array of the Z of the elements con- *
-* stituting the requested Fluka materials *
-* (the Z for the elements of the m_th materials *
-* start at n = Sum_i=1^m-1 [Nelmfl(i)] + 1) *
-* Wfelml(n) = Cumulative array of the weight fractions of the*
-* elements constituting the requested Fluka *
-* materials (the weight fractions for the elem- *
-* ents of the m_th materials start at *
-* n = Sum_i=1^m-1 [Nelmfl(i)] + 1) *
-* Mxelfl = Dimension of the Iz/Wfelml arrays, it must be *
-* Mxelfl >= Sum_i=1^nmatfl [Nelmfl(i)] *
-* Pptmax = Maximum momentum (GeV/c) to be used in initia- *
-* lization (optional) *
-* Ef2dp3 = Transition energy (GeV) from Fluka (Peanut) to *
-* Dpmjet3 for hA interactions (optional) *
-* Df2dp3 = Smearing (+/-Df2dp3) (GeV) of the transition *
-* energy from Fluka (Peanut) to Dpmjet3 for hA *
-* interactions (optional) *
-* Lprint = Material printout flag *
-* *
-* Output variables: *
-* *
-* Mtflka(m) = Fluka material number corresponding to the m_th*
-* requested material *
-* *
-*----------------------------------------------------------------------*
-*/
-void stpxyz_(int const* NMATFL,
- int const NELMFL[],
- int const IZELFL[],
- double const WFELFL[],
- int const* MXELFL,
- double const* PPTMAX,
- double const* EF2DP3,
- double const* DF2DP3,
- int const* IFLXYZ,
- bool const* LPRINT,
- int* MTFLKA,
- char CRVRCK[8]);
-}
-}
+ extern "C" {
+ /* Iflxyz must be consistent in all calls, to Stpxyz, Sgmxyz, Evtxyz
+ * Iflxyz = 1 -> only inelastic
+ * Iflxyz = 10 -> only elastic
+ * Iflxyz = 11 -> inelastic + elastic
+ * Iflxyz =100 -> only emd
+ * Iflxyz =101 -> inelastic + emd
+ * Iflxyz =110 -> elastic + emd
+ * Iflxyz =111 -> inelastic + elastic + emd
+ */
+
+ /*----------------------------------------------------------------------*
+ * *
+ * Copyright (C) 2022-2022 by Alfredo Ferrari & Paola Sala *
+ * All Rights Reserved. *
+ * *
+ * *
+ * SeTuP for X/Y/Z interaction types: *
+ * *
+ * Authors: Alfredo Ferrari & Paola Sala *
+ * *
+ * *
+ * Created on 24 October 2022 by Alfredo Ferrari & Paola Sala *
+ * Private Private *
+ * *
+ * Last change on 01-Nov-22 by Alfredo Ferrari *
+ * Private *
+ * *
+ * Input variables: *
+ * *
+ * Nmatfl = Number of requested Fluka materials *
+ * Nelmfl(m) = Number of elements in the m_th requested Fluka *
+ * material, Nelmfl(m)=1 means simple element, no *
+ * compound *
+ * Izelfl(n) = Cumulative array of the Z of the elements con- *
+ * stituting the requested Fluka materials *
+ * (the Z for the elements of the m_th materials *
+ * start at n = Sum_i=1^m-1 [Nelmfl(i)] + 1) *
+ * Wfelml(n) = Cumulative array of the weight fractions of the*
+ * elements constituting the requested Fluka *
+ * materials (the weight fractions for the elem- *
+ * ents of the m_th materials start at *
+ * n = Sum_i=1^m-1 [Nelmfl(i)] + 1) *
+ * Mxelfl = Dimension of the Iz/Wfelml arrays, it must be *
+ * Mxelfl >= Sum_i=1^nmatfl [Nelmfl(i)] *
+ * Pptmax = Maximum momentum (GeV/c) to be used in initia- *
+ * lization (optional) *
+ * Ef2dp3 = Transition energy (GeV) from Fluka (Peanut) to *
+ * Dpmjet3 for hA interactions (optional) *
+ * Df2dp3 = Smearing (+/-Df2dp3) (GeV) of the transition *
+ * energy from Fluka (Peanut) to Dpmjet3 for hA *
+ * interactions (optional) *
+ * Lprint = Material printout flag *
+ * *
+ * Output variables: *
+ * *
+ * Mtflka(m) = Fluka material number corresponding to the m_th*
+ * requested material *
+ * *
+ *----------------------------------------------------------------------*/
+ void stpxyz_(int const* NMATFL, int const NELMFL[], int const IZELFL[],
+ double const WFELFL[], int const* MXELFL, double const* PPTMAX,
+ double const* EF2DP3, double const* DF2DP3, int const* IFLXYZ,
+ bool const* LPRINT, int* MTFLKA, char CRVRCK[8]);
+
+ /*----------------------------------------------------------------------*
+ * *
+ * Copyright (C) 2022-2022 by Alfredo Ferrari & Paola Sala *
+ * All Rights Reserved. *
+ * *
+ * *
+ * EVenT for X/Y/Z interaction types: *
+ * *
+ * Authors: Alfredo Ferrari & Paola Sala *
+ * *
+ * *
+ * Created on 24 October 2022 by Alfredo Ferrari & Paola Sala *
+ * Private Private *
+ * *
+ * Last change on 28-Oct-22 by Alfredo Ferrari *
+ * Private *
+ * *
+ * Output variables: *
+ * *
+ * Cumsgi(i) = cumulative macroscopic inelastic cross section *
+ * (cm^2/g x rho) for the i_th element of the mmat*
+ * material *
+ * Cumsge(i) = cumulative macroscopic elastic cross section *
+ * (cm^2/g x rho) for the i_th element of the mmat*
+ * material *
+ * Cumsgm(i) = cumulative macroscopic emd cross section *
+ * (cm^2/g x rho) for the i_th element of the mmat*
+ * material *
+ * *
+ * *
+ *----------------------------------------------------------------------*/
+ //~ void evtxyz_ (int& , MMAT , EKIN , PPROJ , TXX, TYY, TZZ,
+ //~ & IFLXYZ, CUMSGI, CUMSGE, CUMSGM )
+
+ /*----------------------------------------------------------------------*
+ * *
+ * Copyright (C) 2022-2022 by Alfredo Ferrari & Paola Sala *
+ * All Rights Reserved. *
+ * *
+ * SiGMa (mb) for X/Y/Z interaction types: *
+ * *
+ * Authors: Alfredo Ferrari & Paola Sala *
+ * *
+ * *
+ * Created on 24 October 2022 by Alfredo Ferrari & Paola Sala *
+ * Private Private *
+ * *
+ * Last change on 27-Oct-22 by Alfredo Ferrari *
+ * Private *
+ * *
+ *----------------------------------------------------------------------*/
+ double sgmxyz_(int const* KPROJ, int const* MMAT, double const* EKIN,
+ double const* PPROJ, int const* IFLXYZ);
+ }
+} // namespace fluka
diff --git a/src/modules/fluka/fluka_codes.dat b/src/modules/fluka/fluka_codes.dat
index 71bad2741..b7f42677d 100644
--- a/src/modules/fluka/fluka_codes.dat
+++ b/src/modules/fluka/fluka_codes.dat
@@ -1,123 +1,60 @@
+# FLUKA particle code conversion table
+# C8 identifier, FLUKA code, has cross-section/can interact
+Proton 1 yes
+AntiProton 2 yes
+Electron 3
+Positron 4
+NuE 5
+NuEBar 6
Photon 7
-Electron 4
-Positron 3
-
+Neutron 8 yes
+AntiNeutron 9 yes
MuPlus 10
MuMinus 11
-Pi0 23
+K0Long 12 yes
PiPlus 13 yes
PiMinus 14 yes
-K0Long 12 yes
KPlus 15 yes
KMinus 16 yes
-Neutron 8 yes
-Proton 1 yes
-AntiProton 2 yes
-K0Short 19 yes
-# Eta 0
Lambda0 17 yes
-SigmaPlus 21
-Sigma0 22
-SigmaMinus 20
-Xi0 34
-XiMinus 36
-OmegaMinus 38
-AntiNeutron 9 yes
Lambda0Bar 18 yes
-SigmaMinusBar 31
-Sigma0Bar 32
-SigmaPlusBar 33
-Xi0Bar 34
-XiPlusBar 37
-OmegaPlusBar 39
-
-# Omega 0
-# Rho0 0
-# RhoPlus 0
-# RhoMinus 0
-# DeltaPlusPlus 0
-# DeltaPlus 0
-# Delta0 0
-# DeltaMinus 0
-# DeltaMinusMinusBar 0
-# DeltaMinusBar 0
-# Delta0Bar 0
-# DeltaPlusBar 0
-# KStar0 0
-# KStarPlus 0
-# KStarMinusBar 0
-# KStar0Bar 0
-NuE 5
-NuEBar 6
+K0Short 19 yes
+SigmaMinus 20 yes
+SigmaPlus 21 yes
+Sigma0 22 yes
+Pi0 23 yes
+K0 24
+K0Bar 25
NuMu 27
NuMuBar 28
-
-D0 16
-DPlus 25
-DMinus 24
-D0Bar 15
-# DSPlus 0
-# DSMinusBar 0
-# EtaC 0
-# DStar0 0
-# DStarPlus 0
-# DStarMinusBar 0
-# DStar0Bar 0
-# DStarSPlus 0
-# DStarsMinusBar 0
-# 0
-# JPsi 0
+SigmaMinusBar 31 yes
+Sigma0Bar 32 yes
+SigmaPlusBar 33 yes
+Xi0 34 yes
+Xi0Bar 35 yes
+XiMinus 36 yes
+XiPlusBar 37 yes
+OmegaMinus 38 yes
+OmegaPlusBar 39 yes
TauPlus 41
TauMinus 42
NuTau 43
NuTauBar 44
-# 0
-# 0
-LambdaCPlus 17
-XiCPlus 34
-XiC0 36
-SigmaCPlus 21
-SigmaCPlus 22
-SigmaC0 20
-XiPrimeCPlus 34
-XiPrimeC0 36
-OmegaC0 38
-
-LambdaCMinusBar 18
-XiCMinusBar 35
-XiC0Bar 37
-SigmaCMinusMinusBar 31
-SigmaCMinusBar 32
-SigmaC0Bar 33
-XiPrimeCMinusBar 35
-XiPrimeC0Bar 37
-OmegaC0Bar 39
-
-# SigmaStarCPlusPlus 0
-# SigmaStarCPlus 0
-# SigmaStarC0 0
-
-# SigmaStarCBar-- 0
-# SigmaStarCMinusBar 0
-# SigmaStarC0Bar 0
-
-# B0 0
-# BPlus 0
-# BMinusBar 0
-# B0Bar 0
-# BS0 0
-# BS0Bar 0
-# BCPlus 0
-# BCMinusBar 0
-# LambdaB0 0
-# SigmaBMinus 0
-# SigmaBPlus 0
-# XiB0 0
-# XiBMinus 0
-# OmegaBMinus 0
-# LambdaB0Bar 0
-# SigmaBPlusBar 0
-# SigmaBMinusBar 0
-# XiB0Bar 0
-# XiBPlusBar 0
-# OmegaBPlusBar 0
+DPlus 45
+DMinus 46
+D0 47
+D0Bar 48
+DsPlus 49
+DsMinus 50
+LambdaCPlus 51
+XiCPlus 52
+XiC0 53
+XiPrimeCPlus 54
+XiPrimeC0 55
+OmegaC0 56
+LambdaCMinusBar 57
+XiCMinusBar 58
+XiC0Bar 59
+XiPrimeCMinusBar 60
+XiPrimeC0Bar 61
+OmegaC0Bar 62
diff --git a/tests/modules/testFluka.cpp b/tests/modules/testFluka.cpp
index 96f265bf6..c154cbbde 100644
--- a/tests/modules/testFluka.cpp
+++ b/tests/modules/testFluka.cpp
@@ -6,10 +6,13 @@
* the license.
*/
-#include <corsika/modules/fluka/InteractionModel.hpp>
+#include <corsika/modules/FLUKA.hpp>
//~ #include <SetupTestEnvironment.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+
#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/framework/geometry/CoordinateSystem.hpp>
#include <corsika/media/Environment.hpp>
@@ -18,21 +21,73 @@
#include <catch2/catch.hpp>
+#include <fstream>
+#include <iomanip>
+
using namespace corsika;
+TEST_CASE("FLUKACodeConversion") {
+ REQUIRE(corsika::fluka::convertToFluka(Code::PiPlus) ==
+ corsika::fluka::FLUKACode::PiPlus);
+ REQUIRE(corsika::fluka::convertToFlukaRaw(Code::PiPlus) == 13);
+ REQUIRE(corsika::fluka::convertToFlukaRaw(Code::Proton) == 1);
+ REQUIRE(corsika::fluka::convertToFlukaRaw(Code::Lambda0) == 17);
+}
+
TEST_CASE("FLUKA") {
- using DummyEnvironmentInterface = IMediumModel;
- using DummyEnvironment = Environment<DummyEnvironmentInterface>;
- using MyHomogeneousModel = HomogeneousMedium<DummyEnvironmentInterface>;
-
- DummyEnvironment env;
- auto& universe = *env.getUniverse();
- CoordinateSystemPtr const& cs = env.getCoordinateSystem();
- universe.setModelProperties<MyHomogeneousModel>(
- 1_kg / (1_m * 1_m * 1_m),
- NuclearComposition(std::vector<Code>{Code::Hydrogen, Code::Oxygen, Code::Nitrogen, Code::Argon},
- std::vector<double>{1./4., 1./4., 1./4., 1./4.}));
-
- corsika::fluka::InteractionModel flukaModel{env};
-
+ using DummyEnvironmentInterface = IMediumModel;
+ using DummyEnvironment = Environment<DummyEnvironmentInterface>;
+ using MyHomogeneousModel = HomogeneousMedium<DummyEnvironmentInterface>;
+
+ DummyEnvironment env;
+ auto& universe = *env.getUniverse();
+ CoordinateSystemPtr const& cs = env.getCoordinateSystem();
+ universe.setModelProperties<MyHomogeneousModel>(
+ 1_kg / (1_m * 1_m * 1_m),
+ NuclearComposition(
+ std::vector<Code>{Code::Hydrogen, Code::Oxygen, Code::Nitrogen, Code::Argon},
+ std::vector<double>{.25, .25, .25, .25}));
+
+ corsika::fluka::InteractionModel flukaModel{env};
+
+ // test getMaterialIndex
+ REQUIRE(flukaModel.getMaterialIndex(Code::Hydrogen) > 0);
+ REQUIRE(flukaModel.getMaterialIndex(Code::Oxygen) > 0);
+ REQUIRE(flukaModel.getMaterialIndex(Code::Nitrogen) > 0);
+ REQUIRE(flukaModel.getMaterialIndex(Code::Argon) > 0);
+ REQUIRE(flukaModel.getMaterialIndex(Code::Uranium) < 0);
+
+ { // test getCrossSection for allowed projectile/target combinations
+ auto combinationsOK = std::vector{
+ std::tuple{Code::PiMinus, Code::Hydrogen},
+ std::tuple{Code::PiMinus, Code::Nitrogen},
+ std::tuple{Code::PiMinus, Code::Oxygen},
+ std::tuple{Code::KMinus, Code::Oxygen},
+ std::tuple{Code::K0Long, Code::Oxygen},
+ std::tuple{Code::K0Short, Code::Oxygen},
+ std::tuple{Code::Lambda0, Code::Oxygen},
+ std::tuple{Code::SigmaPlus, Code::Oxygen},
+ std::tuple{Code::Proton, Code::Oxygen},
+ std::tuple{Code::AntiProton, Code::Oxygen},
+ std::tuple{Code::KMinus, Code::Hydrogen},
+ std::tuple{Code::K0Long, Code::Hydrogen},
+ std::tuple{Code::K0Short, Code::Hydrogen},
+ std::tuple{Code::Lambda0, Code::Hydrogen},
+ std::tuple{Code::SigmaPlus, Code::Hydrogen},
+ std::tuple{Code::Proton, Code::Hydrogen},
+ std::tuple{Code::AntiProton, Code::Hydrogen},
+ };
+
+ HEPEnergyType const p = 100_GeV;
+ for (auto const& [projectileCode, targetCode] : combinations) {
+ auto const projectile4mom =
+ FourVector{calculate_total_energy(p, get_mass(projectileCode)),
+ MomentumVector{cs, 0_eV, 0_eV, p}};
+ auto const target4mom =
+ FourVector{get_mass(targetCode), MomentumVector{cs, 0_eV, 0_eV, 0_eV}};
+
+ CHECK(flukaModel.getCrossSection(projectileCode, targetCode, projectile4mom,
+ target4mom) > 0_mb);
+ }
+ }
}
--
GitLab