From a9e4841230c9d344277bf43802e03bed4fe31c36 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Wed, 13 Feb 2019 19:05:56 +0000
Subject: [PATCH] started nucleus initialization
---
Processes/Sibyll/NuclearInteraction.cc | 84 ++++++++++++++++++++++++++
Processes/Sibyll/NuclearInteraction.h | 5 +-
Processes/Sibyll/nuclib.h | 11 ++++
Processes/Sibyll/signuc.f | 43 +++++++++++++
4 files changed, 142 insertions(+), 1 deletion(-)
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index 6a3736323..64013f054 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -54,8 +54,92 @@ namespace corsika::process::sibyll {
// initialize nuclib
// TODO: make sure this does not overlap with sibyll
nuc_nuc_ini_();
+
+ // initialize cross sections
+ InitializeNuclearCrossSections();
+ }
+
+ void NuclearInteraction::InitializeNuclearCrossSections()
+ {
+ using namespace corsika::particles;
+ using namespace units::si;
+ // TODO: get composition of target volumes from environment
+ // now: hard coded list for air
+ constexpr Code target_nuclei[] = { Code::Oxygen, Code::Nitrogen };
+
+ int fNSample = 500; // number of samples in MC estimation of cross section
+
+ cout << "NuclearInteraction: initializing nuclear cross sections..." << endl;
+
+ // loop over target components, at most 4!!
+ int k =-1;
+ for(auto &ptarg: target_nuclei){
+ ++k;
+ cout << "NuclearInteraction: init target component: " << ptarg << endl;
+ const int ib = GetNucleusA( ptarg );
+ // fill map,list or what ever
+ // TODO
+ // loop over energies, 6 log. energy bins
+ for(int i=0; i<6; ++i){
+ // hard coded energy grid, has to be aligned to definition in signuc2!!, no comment..
+ const units::si::HEPEnergyType Ecm = pow(10., 1. + 1.*i ) * 1_GeV;
+ // get p-p cross sections
+ auto const protonId = Code::Proton;
+ auto const [siginel, sigela, dum ] =
+ fHadronicInteraction.GetCrossSection(protonId, protonId, Ecm);
+ const double dsig = siginel / 1_mbarn;
+ const double dsigela = sigela / 1_mbarn;
+ // loop over projectiles
+ for(int j=0; j<56; ++j){
+ const int jj = j+1;
+ double sig_out, dsig_out, sigqe_out, dsigqe_out;
+ // cout << "energy="<< Ecm/1_GeV << " sig_pp=" << dsig
+ // << " (A,B)=" << jj << "," << ib
+ // << " sig="<<sig_out << " err(%)=" << dsig_out/sig_out*100.
+ // << endl;
+ sigma_mc_(jj,ib,dsig,dsigela,fNSample, sig_out, dsig_out, sigqe_out, dsigqe_out);
+ // write to table
+ cnucsignuc_.sigma[j][k][i] = sig_out;
+ cnucsignuc_.sigqe[j][k][i] = sigqe_out;
+ }
+ }
+ }
+ cout << "NuclearInteraction: cross sections initialized!" << endl;
+ PrintCrossSectionTable(0);
+ PrintCrossSectionTable(1);
+ PrintCrossSectionTable(2);
+ throw std::runtime_error("stop here");
}
+ void NuclearInteraction::PrintCrossSectionTable(int k)
+ {
+ int nuclA [] = {2,4,8,20,35,56};
+ cout << "en/A ";
+ for(int j=0; j<6; ++j)
+ cout << std::setw(8) << nuclA[j];
+ cout << endl;
+
+ for(int i=0; i<6; ++i){
+ cout << " " << i << " ";
+ for(int j=0; j<6; ++j){
+ cout << " " << std::setprecision(5) << std::setw(7) << cnucsignuc_.sigma[nuclA[j]-1][k][i];
+ }
+ cout << endl;
+ }
+ }
+
+ units::si::CrossSectionType NuclearInteraction::ReadCrossSectionTable(particles::Code pBeam, particles::Code pTarget, units::si::HEPEnergyType elabnuc)
+ {
+ using namespace corsika::particles;
+ using namespace units::si;
+ const int ia = GetNucleusA(pBeam);
+ const int ib = GetNucleusA(pTarget);
+ const double e0 = elabnuc/ 1_GeV;
+ double sig;
+ signuc2_(ia,ib,e0,sig);
+ return sig * 1_mbarn;
+ }
+
// TODO: remove number of nucleons, avg target mass is available in environment
template <>
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index cf1f65933..0653a0526 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -40,7 +40,10 @@ namespace corsika::process::sibyll {
corsika::process::sibyll::Interaction& hadint);
~NuclearInteraction();
void Init();
-
+ void InitializeNuclearCrossSections();
+ void PrintCrossSectionTable(int);
+ corsika::units::si::CrossSectionType ReadCrossSectionTable(corsika::particles::Code, corsika::particles::Code, corsika::units::si::HEPEnergyType);
+
template <typename Particle>
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(Particle& p, const corsika::particles::Code TargetId);
diff --git a/Processes/Sibyll/nuclib.h b/Processes/Sibyll/nuclib.h
index e4b89b52b..ddd1856cd 100644
--- a/Processes/Sibyll/nuclib.h
+++ b/Processes/Sibyll/nuclib.h
@@ -34,6 +34,13 @@ extern struct {
*/
extern struct { double ppp[60][3]; } fragments_;
+ // COMMON /cnucsignuc/SIGMA(6,4,56), SIGQE(6,4,56)
+ extern struct
+ {
+ double sigma[56][4][6];
+ double sigqe[56][4][6];
+ } cnucsignuc_;
+
// NUCLIB
// subroutine to initiate nuclib
@@ -46,5 +53,9 @@ void int_nuc_(const int&, const int&, const double&, const double&);
void fragm_(const int&, const int&, const int&, const double&, int&, int*);
void signuc_(const int&, const double&, double&);
+
+void signuc2_(const int&, const int&, const double&, double&);
+
+void sigma_mc_(const int&, const int&, const double&, const double&, const int&, double&, double&, double&, double&);
}
#endif
diff --git a/Processes/Sibyll/signuc.f b/Processes/Sibyll/signuc.f
index b21bce8cb..b88f2ae0b 100644
--- a/Processes/Sibyll/signuc.f
+++ b/Processes/Sibyll/signuc.f
@@ -274,3 +274,46 @@ C ADD INELASTIC AND QUASI-ELASTIC CROSS-SECTIONS
RETURN
END
+
+
+ SUBROUTINE SIGNUC2(IA,IB,E0,SSIGNUC)
+
+C-----------------------------------------------------------------------
+C SIG(MA) NUC(LEUS) 2
+C INPUT: IA : projectile mass number
+C IB : position in cross section table for target nucleus
+C specified when tables are filled
+C E0 : Energy per nucleon in GeV in the lab.
+C OUTPUT: inelastic cross section
+C-----------------------------------------------------------------------
+ IMPLICIT NONE
+ DOUBLE PRECISION SIGMA,SIGQE
+ COMMON /cnucsignuc/SIGMA(6,4,56), SIGQE(6,4,56)
+ DIMENSION AA(6)
+ DOUBLE PRECISION AA,DA,AMIN,ABEAM,S1,S2,ASQS
+ DOUBLE PRECISION E0,SSIGNUC
+ INTEGER IA,IB,J,JE,NE
+ DOUBLE PRECISION QUAD_INT
+ EXTERNAL QUAD_INT
+ SAVE
+ DATA NE /6/, AMIN /1.D0/, DA /1.D0/
+ DATA AA /1.D0,2.D0,3.D0,4.D0,5.D0,6.D0/
+
+C ENERGY E0 IN GEV
+ ASQS = 0.5D0*LOG10(1.876D0*E0)
+ JE = MIN( INT( (ASQS-AMIN)/DA )+1, NE-2 )
+ ABEAM = DBLE(IA)
+C INELASTIC CROSS-SECTION
+ S1 = QUAD_INT( ASQS, AA(JE),AA(JE+1),AA(JE+2),
+ + SIGMA(JE,IB,IA),SIGMA(JE+1,IB,IA),
+ + SIGMA(JE+2,IB,IA) )
+C QUASI ELASTIC CROSS-SECTION
+ S2 = QUAD_INT( ASQS, AA(JE),AA(JE+1),AA(JE+2),
+ + SIGQE(JE,IB,IA),SIGQE(JE+1,IB,IA),
+ + SIGQE(JE+2,IB,IA) )
+C ADD INELASTIC AND QUASI-ELASTIC CROSS-SECTIONS
+ SSIGNUC = S1 + S2
+
+ RETURN
+ END
+
--
GitLab