From 14a82546becef13bc8d006b2ee86c4ec6bca6bd8 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Tue, 28 Feb 2023 16:06:54 +0100
Subject: [PATCH] call event generation
---
.../detail/modules/fluka/InteractionModel.inl | 78 ++++++++++++++-----
corsika/modules/fluka/InteractionModel.hpp | 4 +
modules/fluka/FLUKA.hpp | 24 +++++-
tests/modules/testFluka.cpp | 35 +++++++--
4 files changed, 116 insertions(+), 25 deletions(-)
diff --git a/corsika/detail/modules/fluka/InteractionModel.inl b/corsika/detail/modules/fluka/InteractionModel.inl
index 530cafbe0..2798d9fb5 100644
--- a/corsika/detail/modules/fluka/InteractionModel.inl
+++ b/corsika/detail/modules/fluka/InteractionModel.inl
@@ -16,6 +16,7 @@
#include <utility>
#include <boost/iterator/zip_iterator.hpp>
+#include <Eigen/Dense>
#include <corsika/media/Environment.hpp>
#include <corsika/media/NuclearComposition.hpp>
@@ -29,9 +30,11 @@
namespace corsika::fluka {
template <typename TEnvironment>
inline InteractionModel::InteractionModel(TEnvironment const& env)
- : materials_{genFlukaMaterials(env)} {
+ : materials_{genFlukaMaterials(env)}
+ , cumsgx_{std::make_unique<double[]>(materials_.size() * 3)} {
for (auto const& [code, matno] : materials_) {
- std::cout << get_name(code, full_name{}) << " " << matno << std::endl;
+ CORSIKA_LOGGER_DEBUG(logger_, "FLUKA material initialization: {} -> {}",
+ get_name(code, full_name{}), matno);
}
}
@@ -89,11 +92,8 @@ namespace corsika::fluka {
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';
+ CORSIKA_LOGGER_DEBUG(logger_, fmt::format("Elab = {} GeV", Elab * invGeV));
+ CORSIKA_LOGGER_DEBUG(logger_, fmt::format("EkinLab = {} GeV", EkinLab * invGeV));
CrossSectionType const xs = ::fluka::sgmxyz_(&flukaCodeProj, &flukaMaterial, &EkinLab,
&labMomentum, &iflxyz_) *
@@ -106,7 +106,59 @@ namespace corsika::fluka {
Code const projectileId,
Code const targetId,
FourMomentum const& projectileP4,
- FourMomentum const& targetP4) {}
+ FourMomentum const& targetP4) {
+
+ auto const flukaCodeProj =
+ static_cast<FLUKACodeIntType>(convertToFluka(projectileId));
+ auto const flukaMaterial = getMaterialIndex(targetId);
+
+ HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
+ if (!isValid(projectileId, flukaMaterial, sqrtS)) {
+ std::string const errmsg = fmt::format(
+ "Event generation with invalid configuration requested: proj: {}, target: {}",
+ get_name(projectileId, full_name{}), get_name(targetId, full_name{}));
+ CORSIKA_LOGGER_CRITICAL(logger_, errmsg);
+ throw std::runtime_error{errmsg.c_str()};
+ }
+
+ COMBoost const targetRestBoost{targetP4.getSpaceLikeComponents(), get_mass(targetId)};
+ FourMomentum const projectileLab4mom = targetRestBoost.toCoM(projectileP4);
+ HEPEnergyType const Elab = projectileLab4mom.getTimeLikeComponent();
+ auto constexpr invGeV = 1 / 1_GeV;
+ double const EkinLab = (Elab - get_mass(projectileId)) * invGeV;
+
+ auto const plab = projectileLab4mom.getSpaceLikeComponents();
+ auto const& cs = plab.getCoordinateSystem();
+ auto const labMomentum = plab.getNorm();
+ double const labMomentumGeV = labMomentum * invGeV;
+
+ auto const direction = (plab / labMomentum).getComponents().getEigenVector();
+
+ ::fluka::evtxyz_(&flukaCodeProj, &flukaMaterial, &EkinLab, &labMomentumGeV,
+ &direction[0], &direction[1], &direction[2], &iflxyz_, cumsgx_.get(),
+ cumsgx_.get() + materials_.size(),
+ cumsgx_.get() + materials_.size() * 2);
+
+ //~ extern struct {
+ //~ int nevhep; // event number
+ //~ int nhep; // number of entries
+ //~ hepmc_array<int> isthep; // status code
+ //~ hepmc_array<int> idhep; // PDG particle id
+ //~ hepmc_array<int[2]> jmohep; // position of first, second mother
+ //~ hepmc_array<int[2]> jdahep; // position of first, last daughter
+ //~ hepmc_array<double[5]> phep; // 4-momemtum, mass (GeV)
+ //~ hepmc_array<double[4]> vhep; // vertex, production time in mm
+ //~ } hepevt_;
+
+ for (int i = 0; i < ::fluka::hepevt_.nhep; ++i) {
+ int const pdg = ::fluka::hepevt_.idhep[i];
+ int const status = ::fluka::hepevt_.isthep[i];
+ auto const mom = QuantityVector<hepenergy_d>{
+ Eigen::Map<Eigen::Vector3d>(&::fluka::hepevt_.phep[i][0]) * invGeV.magnitude()};
+
+ std::cout << pdg << '\t' << status << '\t' << mom << std::endl;
+ }
+ }
template <typename TEnvironment>
inline std::vector<std::pair<Code, int>> InteractionModel::genFlukaMaterials(
@@ -146,16 +198,6 @@ namespace corsika::fluka {
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
- * Iflxyz = 11 -> inelastic + elastic
- * Iflxyz =100 -> only emd
- * Iflxyz =101 -> inelastic + emd
- * Iflxyz =110 -> elastic + emd
- * Iflxyz =111 -> inelastic + elastic + emd
- */
-
std::fill(&nelmfl[0], &nelmfl[nElements], 1);
std::fill(&wfelml[0], &wfelml[nElements], 1.);
std::transform(allElementsInUniverse.cbegin(), allElementsInUniverse.cend(),
diff --git a/corsika/modules/fluka/InteractionModel.hpp b/corsika/modules/fluka/InteractionModel.hpp
index d2ecf12f9..f662040bc 100644
--- a/corsika/modules/fluka/InteractionModel.hpp
+++ b/corsika/modules/fluka/InteractionModel.hpp
@@ -10,11 +10,13 @@
#include <vector>
#include <utility>
+#include <memory>
#include <corsika/media/Environment.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/framework/geometry/FourVector.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
namespace corsika::fluka {
@@ -39,6 +41,8 @@ namespace corsika::fluka {
private:
std::vector<std::pair<Code, int>> const
materials_; //!< map target Code to FLUKA material no.
+ std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_FLUKA_Interaction");
+ std::unique_ptr<double[]> cumsgx_; //!< dump for evtxyz cumsg*
template <typename TEnvironment>
static std::vector<std::pair<Code, int>> genFlukaMaterials(TEnvironment const&);
diff --git a/modules/fluka/FLUKA.hpp b/modules/fluka/FLUKA.hpp
index b7f9c3e73..e24f6d22e 100644
--- a/modules/fluka/FLUKA.hpp
+++ b/modules/fluka/FLUKA.hpp
@@ -149,8 +149,9 @@ namespace fluka {
* *
* *
*----------------------------------------------------------------------*/
- //~ void evtxyz_ (int& , MMAT , EKIN , PPROJ , TXX, TYY, TZZ,
- //~ & IFLXYZ, CUMSGI, CUMSGE, CUMSGM )
+ void evtxyz_(int const* KPROJ, int const* MMAT, double const* EKIN, double const* PPROJ,
+ double const* TXX, double const* TYY, double const* TZZ, int const* IFLXYZ,
+ double CUMSGI[], double CUMSGE[], double CUMSGM[]);
/*----------------------------------------------------------------------*
* *
@@ -171,5 +172,24 @@ namespace fluka {
*----------------------------------------------------------------------*/
double sgmxyz_(int const* KPROJ, int const* MMAT, double const* EKIN,
double const* PPROJ, int const* IFLXYZ);
+
+ /*----------------------------------------------------------------------*
+ * *
+ * Copyright (C) 2023-2023 by Alfredo Ferrari & Paola Sala *
+ * All Rights Reserved. *
+ * *
+ * FiLL HEP common: *
+ * *
+ * Authors: Alfredo Ferrari & Paola Sala *
+ * *
+ * *
+ * Created on 06 February 2023 by Alfredo Ferrari & Paola Sala *
+ * Private Private *
+ * *
+ * Last change on 07-Feb-23 by Alfredo Ferrari *
+ * Private *
+ * *
+ *----------------------------------------------------------------------*/
+ void fllhep_();
}
} // namespace fluka
diff --git a/tests/modules/testFluka.cpp b/tests/modules/testFluka.cpp
index 818ce0963..0770bf06c 100644
--- a/tests/modules/testFluka.cpp
+++ b/tests/modules/testFluka.cpp
@@ -24,6 +24,9 @@
#include <fstream>
#include <iomanip>
+#include <SetupTestStack.hpp>
+#include <SetupTestEnvironment.hpp>
+
using namespace corsika;
TEST_CASE("FLUKACodeConversion") {
@@ -35,18 +38,20 @@ TEST_CASE("FLUKACodeConversion") {
}
TEST_CASE("FLUKA") {
- using DummyEnvironmentInterface = IMediumModel;
+ using DummyEnvironmentInterface =
+ IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>;
using DummyEnvironment = Environment<DummyEnvironmentInterface>;
- using MyHomogeneousModel = HomogeneousMedium<DummyEnvironmentInterface>;
+ using MyHomogeneousModel = MediumPropertyModel<
+ UniformMagneticField<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(
+ Medium::AirDry1Atm, Vector(cs, 0_T, 0_T, 0_T), 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}));
+ std::vector<double>{.25, .25, .25, .25}});
corsika::fluka::InteractionModel flukaModel{env};
@@ -90,4 +95,24 @@ TEST_CASE("FLUKA") {
target4mom) > 0_mb);
}
}
+
+ {
+ auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
+ auto const& cs = *csPtr;
+ auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+ Code::Hydrogen, 1_GeV, (DummyEnvironment::BaseNodeType* const)nodePtr, *csPtr);
+ { [[maybe_unused]] auto const& dummy_StackPtr = stackPtr; }
+
+ auto const projectileCode = Code::PiPlus;
+ auto const targetCode = Code::Nitrogen;
+ auto const p = 20_GeV;
+ 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}};
+
+ flukaModel.doInteraction(*secViewPtr, projectileCode, targetCode, projectile4mom,
+ target4mom);
+ }
}
--
GitLab