From 16c98b1fe6524982a7a7e905bed6e0cd8c67d3a2 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Tue, 11 May 2021 19:09:49 +0100
Subject: [PATCH] finished corsika style init and event init
---
corsika/detail/modules/epos/Interaction.inl | 236 +++++++++++++-------
1 file changed, 155 insertions(+), 81 deletions(-)
diff --git a/corsika/detail/modules/epos/Interaction.inl b/corsika/detail/modules/epos/Interaction.inl
index f6c1790c5..4dedf78cf 100644
--- a/corsika/detail/modules/epos/Interaction.inl
+++ b/corsika/detail/modules/epos/Interaction.inl
@@ -21,6 +21,7 @@
#include <epos.hpp>
+#include <string>
#include <tuple>
using namespace corsika;
@@ -28,81 +29,121 @@ using SetupParticle = setup::Stack::stack_iterator_type;
namespace corsika::epos {
- inline Interaction::Interaction() {
-
- // if(ilowegy.ne.1.or.MCleModel.eq.4)xsegymin=dble(0.5*egymin**2)
- // if(MCModel.eq.4)xsegymax=min(xsegymax,dble(0.5*egymax**2))
- // nrnody=nrnodyxs
- // do i=1,nrnody
- // nody(i)= nodyxs(i)
- // enddo
- // #if !__CXCORSIKA__ && !__CORSIKA8__
-
- // inicnt=inicnt+1
- // c isetcs=2 ! epos cross-section from tabulated calculation (h-A and AA)
- ::epos::hadr6_.isetcs=3; // epos cross-section from tabulated simulations
- // (h-A and A-A)
- ::epos::nucl6_.infragm=2; // --> model how projectiles fragment!
- ::epos::hadr6_.isigma=0; // !do not print out the
- // cross section on screen ionudi=1
-
- // nfnii=nxsfnii ! epos file name
- // fnii=xsfnii
- // nfnid=nxsfnid
- // fnid=xsfnid
- // nfnie=nxsfnie
- // fnie=xsfnie
- // nfnrj=nxsfnrj
- // fnrj=xsfnrj
- // nfncs=nxsfncs
- // fncs=xsfncs
- // nfnch=nxsfnch
- // fnch=xsfnch
-
- // c air
- /* for (int i=1; i<=3; ++i) {
- ::epos::nxsair_.airanxs[i]=aira[i];
- ::epos::nxsair_.airznxs[i]=airz[i];
- ::epos::nxsair_.airwnxs[i]=airw[i];
- }
- ::epos::nxsair_.airavanxs=airava;
- ::epos::nxsair_.airavznxs=airavz;
- */
-
- ::epos::appli_.iappl = iapplxs;
- ::epos::events_.nevent = neventxs;
- ::epos::othe2_.iframe = iframexs;
-
- // if(fnch(1:nfnch).ne.'none')
- // & open(ifcx,file=fnch(1:nfnch),status='unknown')
-
- // call iclass(idproj,iclpro)
- // call iclass(idtarg,icltar)
- // if(inicnt.eq.1)then
- // call ranfgt(seedp) !not to change the seed ...
- ::epos::hdecin_(false);
- ::epos::hnbspd_(iospec);
- // ktnbod=0
- ::epos::hnbpajini_();
- // if(iclegy2.gt.1)then
- // egyfac=(egymax*1.0001/egylow)**(1./float(iclegy2-1))
- // else
- // egyfac=1.
- // endif
- // endif
- // maproj=mamx !to set difnuc up to the maximum mass
- // call conini
- // call psaini
- // call ranfst(seedp) ! ... after this initialization
+ inline Interaction::Interaction(const std::string& dataPath)
+ : data_path_(dataPath) {
+ if (dataPath == "") {
+ if (std::getenv("CORSIKA_DATA")) {
+ data_path_ = std::string(std::getenv("CORSIKA_DATA")) + "/EPOS/";
+ CORSIKA_LOG_DEBUG("Searching for EPOSLHC data tables in {}", data_path_);
+ }
+ }
+
+ // initialize Eposlhc
+ static bool initialized = false;
+ if (!initialized) {
+ initialize_eposlhc_c7();
+ initialized = true;
+ }
+ }
+
+
+ inline void Interaction::initialize_eposlhc_c7() {
+
+ // corsika7 ini
+ int iarg = 0;
+ ::epos::aaset_(iarg);
+ //::epos::atitle_();
+
+ //::epos::prnt1_.ish = 3; // debug level in epos
+ //::epos::prnt1_.ifch = 6; // output unit
+ //::epos::prnt1_.iecho = 1;
+
+ // dummy set seeds for random number generator in epos. need to fool epos checks...
+ // we will use external generator
+ ::epos::cseed_.seedi = 1;
+ ::epos::cseed_.seedj = 1;
+ ::epos::cseed_.seedc = 1;
+
+ ::epos::enrgy_.egymin = 6.;
+ ::epos::enrgy_.egymax = 2.e6;
+
+ ::epos::lhcparameters_();
+
+ ::epos::hadr6_.isigma = 0;// do not show cross section
+ ::epos::hadr6_.isetcs = 3; /* !option to obtain pomeron parameters
+ ! 0.....determine parameters but do not use Kfit
+ ! 1.....determine parameters and use Kfit
+ ! else..get from table
+ ! should be sufficiently detailed
+ ! say iclegy1=1,iclegy2=99
+ ! table is always done, more or less detailed!!!
+ !and option to use cross section tables
+ ! 2....tabulation
+ ! 3....simulation
+ */
+ ::epos::cjinti_.ionudi =
+ 1; // !include quasi elastic events but strict calculation of xs
+ ::epos::cjinti_.iorsce = 0; // !color exchange turned on(1) or off(0)
+ ::epos::cjinti_.iorsdf = 3; // !droplet formation turned on(>0) or off(0)
+ ::epos::cjinti_.iorshh = 0; // !other hadron-hadron int. turned on(1) or off(0)
+
+ ::epos::othe1_.istore = 0; // do not produce epos output file
+ ::epos::nucl6_.infragm = 0; // keep free nucleons in fragmentation
+
+ // set paths to tables in corsika data
+ ::epos::datadir BASE(data_path_);
+ strcpy(::epos::fname_.fnnx, BASE.data);
+ ::epos::nfname_.nfnnx = BASE.length;
+
+ ::epos::datadir TL(data_path_ + "epos.initl");
+ strcpy(::epos::fname_.fnii, TL.data);
+ ::epos::nfname_.nfnii = TL.length;
+
+ ::epos::datadir EV(data_path_ + "epos.iniev");
+ strcpy(::epos::fname_.fnie, EV.data);
+ ::epos::nfname_.nfnie = EV.length;
+
+ ::epos::datadir RJ(data_path_ + "epos.inirj"); // lhcparameters adds ".lhc"
+ strcpy(::epos::fname_.fnrj, RJ.data);
+ ::epos::nfname_.nfnrj = RJ.length;
+
+ ::epos::datadir CS(data_path_ + "epos.inics"); // lhcparameters adds ".lhc"
+ strcpy(::epos::fname_.fncs, CS.data);
+ ::epos::nfname_.nfncs = CS.length;
+
+ //::epos::fname_.fnid="/home/felix/ngcorsika/corsika/modules/data/EPOS/epos.inidi";
+ // EPOPAR input ../epos/epos.param !initialization input file for epos
+ // EPOPAR fname inics ../epos/epos.inics !initialization input file for epos
+ // EPOPAR fname iniev ../epos/epos.iniev !initialization input file for epos
+ // EPOPAR fname initl ../epos/epos.initl !initialization input file for epos
+ // EPOPAR fname inirj ../epos/epos.inirj !initialization input file for epos
+ // EPOPAR fname inihy ../epos/epos.ini1b !initialization input file for epos
+
+ // dummy event
+ ::epos::hadr25_.idprojin = 1120;
+ ::epos::hadr25_.idtargin = 1120;
+ //#if __CONEX__ && __EPOS__ && __HIGHMEM__
+ // maproj = 250
+ //#else
+ ::epos::nucl1_.maproj = 56;
+ // #endif
+ ::epos::nucl1_.laproj = 28;
+ ::epos::nucl1_.matarg = 14;
+ ::epos::nucl1_.latarg = 1;
+ ::epos::hadr1_.pnll = 200.;
+ ::epos::lept1_.engy = -1.;
+
+ ::epos::ainit_();
+
}
inline Interaction::~Interaction() {
CORSIKA_LOG_DEBUG("Epos::Interaction n={} ", count_);
}
- inline corsika::CrossSectionType Interaction::getCrossSection(
- const corsika::Code BeamId, const corsika::Code TargetId,
- const corsika::HEPEnergyType CoMenergy) const {}
+ inline std::tuple<corsika::CrossSectionType, corsika::CrossSectionType>
+ Interaction::getCrossSection(const corsika::Code BeamId, const corsika::Code TargetId,
+ const corsika::HEPEnergyType CoMenergy) const {}
template <>
inline corsika::GrammageType Interaction::getInteractionLength(
@@ -111,23 +152,18 @@ namespace corsika::epos {
return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
}
- /**
- In this function SIBYLL is called to produce one event. The
- event is copied (and boosted) into the shower lab frame.
- */
-
template <typename TSecondaryView>
inline void Interaction::doInteraction(TSecondaryView& view) {
auto const projectile = view.getProjectile();
- const auto corsikaBeamId = projectile.getPID();
+ auto const corsikaBeamId = projectile.getPID();
if (corsika::is_nucleus(corsikaBeamId)) {
// nuclei handled by different process, this should not happen
- throw std::runtime_error("Nuclear projectile are not handled by SIBYLL!");
+ throw std::runtime_error("Nuclear projectile are not handled by EPOSLHC!");
}
- // position and time of interaction, not used in Sibyll
+ // position and time of interaction, not used in Epos
Point const pOrig = projectile.getPosition();
TimeType const tOrig = projectile.getTime();
@@ -138,7 +174,7 @@ namespace corsika::epos {
CORSIKA_LOG_DEBUG(
"ProcessEPOS: "
- "DoInteraction: pid {} interaction ",
+ "DoInteraction: {} interaction ",
corsikaBeamId);
// define target
@@ -207,7 +243,8 @@ namespace corsika::epos {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
- auto const sigProd = getCrossSection(corsikaBeamId, targetId, Ecm);
+ [[maybe_unused]] auto const [sigProd, sigEla] =
+ getCrossSection(corsikaBeamId, targetId, Ecm);
cross_section_of_components[i] = sigProd;
}
@@ -215,8 +252,45 @@ namespace corsika::epos {
mediumComposition.sampleTarget(cross_section_of_components, RNG_);
CORSIKA_LOG_DEBUG("Interaction: target selected: {} ", targetCode);
- //
-
+ // from corsika7 interface
+ // NEXLNK-part
+
+ //projectile
+ // if(is_nucleus(corsikaBeamId)){
+ // ::epos::hadr25_.idprojin =1120;
+ // ::epos::nucl1_.laproj = projectile.get_nucleus_Z(); // Z
+ // ::epos::nucl1_.maproj = projectile.get_nucleus_A();; // A
+ // } else {
+ ::epos::hadr25_.idprojin = 1120 ; // id "NEXUS code"
+ ::epos::nucl1_.laproj = -1; // Z (-1 for hadron)
+ ::epos::nucl1_.maproj = 1; // A
+ //}
+ // target
+ //if(is_nucleus(targetCode)){
+ ::epos::hadr25_.idtargin = 1120; // id "NEXUS code"
+ ::epos::nucl1_.latarg = 1; // Z (-1 with id 1220 for neutron)
+ ::epos::nucl1_.matarg = 1; // A
+
+ // hadron-nucleon momentum
+ ::epos::hadr1_.pnll = float(200);
+
+ // C SET ENGY NEGATIVE TO FORCE CALCULATION IN LAB FRAME
+ ::epos::lept1_.engy = -1.;
+ ::epos::enrgy_.ecms = -1.;
+ ::epos::enrgy_.elab = -1.;
+ ::epos::enrgy_.ekin = -1.;
+
+ //C INTIALIZE ENERGY AND PARTICLE DEPENDENT PORTION OF EPOS/NEXUS
+ //C AT THE FIRST CALL: READ ALSO DATA SETS
+ ::epos::ainit_();
+
+ // create event
+ int iarg = 1;
+ ::epos::aepos_(iarg);
+
+ // NSTORE-part
+ std::cout << "npart: " << ::epos::cptl_.nptl << std::endl;
+
/*
maproj=maprojxs
laproj=laprojxs
--
GitLab