From 63528c40e3f1341e564e63b3c41e7512c71fc46b Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Tue, 11 May 2021 19:08:30 +0100
Subject: [PATCH] corsika/detail/modules/epos/Interaction.inl
---
corsika/modules/epos/Interaction.hpp | 29 ++++--
corsika/modules/epos/Random.hpp | 8 +-
modules/epos/epos.cpp | 14 ++-
modules/epos/epos.hpp | 135 ++++++++++++++++++++++++++-
tests/modules/testEpos.cpp | 37 +++++++-
5 files changed, 207 insertions(+), 16 deletions(-)
diff --git a/corsika/modules/epos/Interaction.hpp b/corsika/modules/epos/Interaction.hpp
index 7c5e179dd..9505e6145 100644
--- a/corsika/modules/epos/Interaction.hpp
+++ b/corsika/modules/epos/Interaction.hpp
@@ -17,29 +17,46 @@
namespace corsika::epos {
class Interaction : public InteractionProcess<Interaction> {
-
+ std::string data_path_;
+ unsigned int count_ = 0;
+
public:
- Interaction();
+ Interaction(const std::string& dataPath = "");
~Interaction();
- //! returns production and elastic cross section for hadrons in sibyll. Inputs are:
+ //! returns production and elastic cross section for hadrons in epos. Inputs are:
//! CorsikaId of beam particle, CorsikaId of target particle and center-of-mass
//! energy. Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
- CrossSectionType getCrossSection(Code const, Code const, HEPEnergyType const) const;
+ std::tuple<CrossSectionType, CrossSectionType> getCrossSection(
+ Code const, Code const, HEPEnergyType const) const;
template <typename TParticle>
GrammageType getInteractionLength(TParticle const&) const;
/**
- In this function SIBYLL is called to produce one event. The
+ In this function EPOSLHC is called to produce one event. The
event is copied (and boosted) into the shower lab frame.
*/
template <typename TSecondaries>
void doInteraction(TSecondaries&);
+ bool isValidCoMEnergy(HEPEnergyType const ecm) const {
+ return (minEnergyCoM_ <= ecm) && (ecm <= maxEnergyCoM_);
+ }
+ //! eposlhc only accepts nuclei with 4<=A<=18 as targets, or protons aka Hydrogen or
+ //! neutrons (p,n == nucleon)
+ bool isValidTarget(Code const TargetId) const {
+ return false;
+ }
+
+ void initialize_eposlhc_c7();
+
private:
default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("epos");
- unsigned int count_;
+ HEPEnergyType const minEnergyCoM_ = -10. * 1e9 * electronvolt;
+ HEPEnergyType const maxEnergyCoM_ = -1.e6 * 1e9 * electronvolt;
+ int const maxTargetMassNumber_ = -1;
+ int const minNuclearTargetA_ = -10;
};
} // namespace corsika::epos
diff --git a/corsika/modules/epos/Random.hpp b/corsika/modules/epos/Random.hpp
index 1e7e7bc7c..c3fec35ab 100644
--- a/corsika/modules/epos/Random.hpp
+++ b/corsika/modules/epos/Random.hpp
@@ -14,17 +14,17 @@
/**
* \file epos/Random.hpp
*
- * This file is an integral part of the sibyll interface. It must be
- * linked to the executable linked to sibyll exactly once
+ * This file is an integral part of the epos interface. It must be
+ * linked to the executable linked to epos exactly once
*
*/
namespace epos {
- double rndm_interface() {
+ float rndm_interface() {
static corsika::default_prng_type& rng =
corsika::RNGManager::getInstance().getRandomStream("epos");
- std::uniform_real_distribution<double> dist;
+ std::uniform_real_distribution<float> dist;
return dist(rng);
}
diff --git a/modules/epos/epos.cpp b/modules/epos/epos.cpp
index 757fba685..1836776c7 100644
--- a/modules/epos/epos.cpp
+++ b/modules/epos/epos.cpp
@@ -1,4 +1,5 @@
#include <epos.hpp>
+#include <iostream>
namespace epos {
@@ -8,7 +9,18 @@ namespace epos {
// this is needed as linker object, but it is not needed to do anything
void ranfcv_(double&) {}
- double rangen_() { return ::epos::rndm_interface(); }
+ float rangen_() { return ::epos::rndm_interface(); }
+ datadir::datadir(const std::string& dir) {
+ if (dir.length() > 500) {
+ std::cerr << "Epos error, will cut datadir \"" << dir
+ << "\" to 500 characters: " << std::endl;
+ }
+ int i = 0;
+ for (i = 0; i < std::min(500, int(dir.length())); ++i) data[i] = dir[i];
+ data[i + 0] = ' ';
+ data[i + 1] = '\0';
+ length = dir.length();
+ }
}
diff --git a/modules/epos/epos.hpp b/modules/epos/epos.hpp
index 47c2939d2..d7799369a 100644
--- a/modules/epos/epos.hpp
+++ b/modules/epos/epos.hpp
@@ -9,6 +9,7 @@
#pragma once
#include <array>
+#include <string>
/**
* \file epos.hpp
@@ -24,13 +25,15 @@ namespace epos {
*
* CORSIKA8, for example, has to provide an implementation of this.
**/
- extern double rndm_interface();
+ extern float rndm_interface();
extern "C" {
void aaset_(int&);
void atitle_();
- double LHCparameters_();
+ void ainit_();
+ void aepos_(int&);
+ double lhcparameters_();
void hdecin_(bool&);
void hnbspd_(int&);
void hnbpajini_();
@@ -59,8 +62,7 @@ namespace epos {
void ranfini_(double&, int&, int&);
void ranfcv_(double&);
// void ranfgt(int&);
- double rangen_();
-
+ float rangen_();
// common blocks as
// defined in epos.inc
@@ -95,6 +97,12 @@ namespace epos {
int model;
} appli_;
+ extern struct {
+ int iapplxs;
+ int modelxs;
+ } xsappli_;
+
+
extern struct {
int nevent;
int nfull;
@@ -102,6 +110,21 @@ namespace epos {
int ninicon;
} events_;
+ extern struct {
+ int neventxs;
+ int iframexs;
+ } xsevent_;
+
+ // common/metr1/iospec,iocova,iopair,iozero,ioflac,iomom
+ extern struct {
+ int iospec;
+ int iocova;
+ int iopair;
+ int iozero;
+ int ioflac;
+ int iomom;
+ } metr1_;
+
extern struct {
int ifrade;
int iframe;
@@ -272,5 +295,109 @@ namespace epos {
// bimevt, kolevt, koievt, pmxevt, egyevt, npjevt , ntgevt, npnevt, nppevt, ntnevt,
// ntpevt, jpnevt, jppevt, jtnevt, jtpevt , xbjevt, qsqevt, nglevt, zppevt, zptevt,
// minfra, maxfra, kohevt
+
+ // common/cseed/seedi,seedj,seedj2,seedc,iseqini,iseqsim
+ extern struct {
+ double seedi;
+ double seedj;
+ double seedj2;
+ double seedc;
+ int iseqini;
+ int iseqsim;
+ } cseed_;
+
+ // integer istore,istmax,irescl,ntrymx,nclean,iopdg,ioidch
+ // common/othe1/istore,istmax,gaumx,irescl,ntrymx,nclean,iopdg,ioidch
+ extern struct {
+ int istore;
+ int istmax;
+ int gaumx;
+ int irescl;
+ int ntrymx;
+ int nclean;
+ int iopdg;
+ int ioidch;
+ } othe1_;
+
+ // character*500 fnch,fnhi,fndt,fnii,fnid,fnie,fnrj,fnmt
+ // * ,fngrv,fncp,fnnx,fncs,fndr,fnhpf
+ // common/fname/ fnch, fnhi, fndt, fnii, fnid, fnie, fnrj, fnmt
+ // * ,fngrv,fncp,fnnx,fncs,fndr,fnhpf
+ extern struct {
+ char fnch[500];
+ char fnhi[500];
+ char fndt[500];
+ char fnii[500];
+ char fnid[500];
+ char fnie[500];
+ char fnrj[500];
+ char fnmt[500];
+ char fngrv[500];
+ char fncp[500];
+ char fnnx[500];
+ char fncs[500];
+ char fndr[500];
+ char fnhpf[500];
+
+ } fname_;
+
+ // integer nfnch,nfnhi,nfndt,nfnii,nfnid,nfnie,nfnrj,nfnmt
+ // *,nfngrv,nfncp,nfnnx,nfncs,nfndr,nfnhpf
+ // common/nfname/nfnch,nfnhi,nfndt,nfnii,nfnid,nfnie,nfnrj,nfnmt
+ // *,nfngrv,nfncp,nfnnx,nfncs,nfndr,nfnhpf
+ extern struct {
+ int nfnch;
+ int nfnhi;
+ int nfndt;
+ int nfnii;
+ int nfnid;
+ int nfnie;
+ int nfnrj;
+ int nfnmt;
+ int nfngrv;
+ int nfncp;
+ int nfnnx;
+ int nfncs;
+ int nfndr;
+ int nfnhpf;
+ } nfname_;
+
+ // integer iprmpt,ish,ishsub,irandm,irewch,iecho,modsho,idensi
+ // common/prnt1/iprmpt,ish,ishsub,irandm,irewch,iecho,modsho,idensi
+ extern struct {
+ int iprmpt;
+ int ish;
+ int ishsub;
+ int irandm;
+ int irewch;
+ int iecho;
+ int modsho;
+ int idensi;
+ } prnt1_;
+ unsigned int constexpr mmry = 1;
+ unsigned int constexpr mxptl = 200000 / mmry;
+ // real pptl,tivptl,xorptl
+ // integer nptl,iorptl,idptl,istptl,ifrptl,jorptl,ibptl,ityptl
+ // common/cptl/nptl,pptl(5,mxptl),iorptl(mxptl),idptl(mxptl)
+ extern struct {
+ int nptl;
+ float pptl[mxptl][5];
+ int iorptl[mxptl];
+ int idptl[mxptl];
+ } cptl_;
+
+ /**
+ Small helper class to provide a data-directory name in the format eposlhc expects
+ */
+ class datadir {
+ private:
+ datadir operator=(const std::string& dir);
+ datadir operator=(const datadir&);
+
+ public:
+ datadir(const std::string& dir);
+ char data[500];
+ int length;
+ };
}
} // namespace epos
diff --git a/tests/modules/testEpos.cpp b/tests/modules/testEpos.cpp
index f0b01ce36..e00ca7dbd 100644
--- a/tests/modules/testEpos.cpp
+++ b/tests/modules/testEpos.cpp
@@ -81,7 +81,38 @@ TEST_CASE("EposInterface", "[processes]") {
auto const& cs = *csPtr;
[[maybe_unused]] auto const& env_dummy = env;
- // RNGManager::getInstance().registerRandomStream("epos");
+ RNGManager::getInstance().registerRandomStream("epos");
+
+ SECTION("InteractionInterface - random number"){
+ auto const rndm = ::epos::rangen_();
+ CHECK(rndm>0);
+ CHECK(rndm<1);
+ }
+
+ SECTION("InteractionInterface - valid targets") {
+
+ Interaction model;
+ // eposlhc accepts protons or nuclei with 4<=A<=18 as targets
+ CHECK_FALSE(model.isValidTarget(Code::Electron));
+ CHECK(model.isValidTarget(Code::Hydrogen));
+ CHECK_FALSE(model.isValidTarget(Code::Deuterium));
+ CHECK(model.isValidTarget(Code::Helium));
+ CHECK_FALSE(model.isValidTarget(Code::Helium3));
+ CHECK_FALSE(model.isValidTarget(Code::Iron));
+ CHECK(model.isValidTarget(Code::Oxygen));
+
+ // hydrogen target == proton target == neutron target
+ auto const [xs_prod_pp, xs_ela_pp] =
+ model.getCrossSection(Code::Proton, Code::Proton, 100_GeV);
+ auto const [xs_prod_pn, xs_ela_pn] =
+ model.getCrossSection(Code::Proton, Code::Neutron, 100_GeV);
+ auto const [xs_prod_pHydrogen, xs_ela_pHydrogen] =
+ model.getCrossSection(Code::Proton, Code::Hydrogen, 100_GeV);
+ CHECK(xs_prod_pp == xs_prod_pHydrogen);
+ CHECK(xs_prod_pp == xs_prod_pn);
+ CHECK(xs_ela_pp == xs_ela_pHydrogen);
+ CHECK(xs_ela_pn == xs_ela_pHydrogen);
+ }
SECTION("InteractionInterface - low energy") {
@@ -96,6 +127,10 @@ TEST_CASE("EposInterface", "[processes]") {
Interaction model;
model.doInteraction(view);
+
+ [[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
+ CHECK(length / 1_g * 1_cm * 1_cm == Approx(88.7).margin(0.1));
+
}
}
--
GitLab