IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 14a82546 authored by Maximilian Reininghaus's avatar Maximilian Reininghaus :vulcan:
Browse files

call event generation

parent 4bcd7062
No related branches found
No related tags found
1 merge request!468Resolve "Add FLUKA"
Hide whitespace changes
Inline Side-by-side
corsika/detail/modules/fluka/InteractionModel.inl
+ 60
− 18
View file @ 14a82546
......@@ -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(),
......
corsika/modules/fluka/InteractionModel.hpp
+ 4
− 0
View file @ 14a82546
......@@ -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&);
......
modules/fluka/FLUKA.hpp
+ 22
− 2
View file @ 14a82546
......@@ -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
tests/modules/testFluka.cpp
+ 30
− 5
View file @ 14a82546
......@@ -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);
}
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum