diff --git a/Environment/NuclearComposition.h b/Environment/NuclearComposition.h
index 5f89bd2cadb50e48f01931abafbc4ee94305cf8b..dee18a4463c340698bbff372cc504d4744748375 100644
--- a/Environment/NuclearComposition.h
+++ b/Environment/NuclearComposition.h
@@ -35,15 +35,15 @@ namespace corsika::environment {
public:
using value_type = double;
using iterator_category = std::input_iterator_tag;
- using pointer = double*;
- using reference = double&;
+ using pointer = value_type*;
+ using reference = value_type&;
using difference_type = ptrdiff_t;
WeightProviderIterator(AConstIterator a, BConstIterator b)
: fAIter(a)
, fBIter(b) {}
- double operator*() const { return ((*fAIter) * (*fBIter)).magnitude(); }
+ value_type operator*() const { return ((*fAIter) * (*fBIter)).magnitude(); }
WeightProviderIterator& operator++() { // prefix ++
++fAIter;
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
index 3331c41bf5e060bb904f5a78f78485af0440d04b..9cd3be1a9e534c5f6b74bb4e584dac1c78e0d9bf 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -66,6 +66,24 @@ namespace corsika::environment {
}
}
+ /**
+ * Traverses the VolumeTree pre- or post-order and calls the functor \p func for each
+ * node. \p func takes a reference to VolumeTreeNode as argument. The return value \p
+ * func is ignored.
+ */
+ template <typename TCallable, bool preorder = true>
+ void walk(TCallable func) {
+ if constexpr (preorder) {
+ func(*this);
+ std::for_each(fChildNodes.begin(), fChildNodes.end(),
+ [&](auto& v) { v->walk(func); });
+ } else {
+ std::for_each(fChildNodes.begin(), fChildNodes.end(),
+ [&](auto& v) { v->walk(func); });
+ func(*this);
+ }
+ }
+
void AddChild(VTNUPtr pChild) {
pChild->fParentNode = this;
fChildNodes.push_back(std::move(pChild));
@@ -88,6 +106,8 @@ namespace corsika::environment {
auto const& GetModelProperties() const { return *fModelProperties; }
+ IMPSharedPtr GetModelPropertiesPtr() const { return fModelProperties; }
+
template <typename TModelProperties, typename... Args>
auto SetModelProperties(Args&&... args) {
static_assert(std::is_base_of_v<IModelProperties, TModelProperties>,
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 5ccc1b743ce35f1e72a1ebd1e5cccba95a8d25e9..0e65c074bba270dcd4db98866e079893d80bf622 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -53,7 +53,7 @@ namespace corsika::process::sibyll {
}
}
- tuple<units::si::CrossSectionType, units::si::CrossSectionType, int>
+ tuple<units::si::CrossSectionType, units::si::CrossSectionType>
Interaction::GetCrossSection(const particles::Code BeamId,
const particles::Code TargetId,
const units::si::HEPEnergyType CoMenergy) {
@@ -61,23 +61,24 @@ namespace corsika::process::sibyll {
double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3];
const int iBeam = process::sibyll::GetSibyllXSCode(BeamId);
+ if( !IsValidCoMEnergy(CoMenergy) ){
+ throw std::runtime_error("Interaction: GetCrossSection: CoM energy outside range for Sibyll!");
+ }
const double dEcm = CoMenergy / 1_GeV;
- int iTarget = -1;
if (particles::IsNucleus(TargetId)) {
- iTarget = particles::GetNucleusA(TargetId);
- if (iTarget > 18 || iTarget == 0)
+ const int iTarget = particles::GetNucleusA(TargetId);
+ if (iTarget > fMaxTargetMassNumber || iTarget == 0)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 are allowed.");
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
} else if (TargetId == particles::Proton::GetCode()) {
sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
- iTarget = 1;
} else {
// no interaction in sibyll possible, return infinite cross section? or throw?
sigProd = std::numeric_limits<double>::infinity();
sigEla = std::numeric_limits<double>::infinity();
}
- return std::make_tuple(sigProd * 1_mbarn, sigEla * 1_mbarn, iTarget);
+ return std::make_tuple(sigProd * 1_mbarn, sigEla * 1_mbarn);
}
template <>
@@ -97,14 +98,11 @@ namespace corsika::process::sibyll {
// read from cross section code table
const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId);
- const HEPMassType nucleon_mass =
- 0.5 * (particles::Proton::GetMass() + particles::Neutron::GetMass());
-
// FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy
- HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
+ HEPEnergyType Elab = p.GetEnergy() + constants::nucleonMass;
MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
pTotLab += p.GetMomentum();
pTotLab += pTarget;
@@ -119,7 +117,8 @@ namespace corsika::process::sibyll {
<< " beam pid:" << p.GetPID() << endl;
// TODO: move limits into variables
- if (kInteraction && Elab >= 8.5_GeV && ECoM >= 10_GeV) {
+ // FR: removed && Elab >= 8.5_GeV
+ if (kInteraction && IsValidCoMEnergy(ECoM)) {
// get target from environment
/*
@@ -134,7 +133,6 @@ namespace corsika::process::sibyll {
// determine average interaction length
// weighted sum
int i = -1;
- double avgTargetMassNumber = 0.;
si::CrossSectionType weightedProdCrossSection = 0_mbarn;
// get weights of components from environment/medium
const auto w = mediumComposition.GetFractions();
@@ -143,7 +141,7 @@ namespace corsika::process::sibyll {
i++;
cout << "Interaction: get interaction length for target: " << targetId << endl;
- auto const [productionCrossSection, elaCrossSection, numberOfNucleons] =
+ auto const [productionCrossSection, elaCrossSection ] =
GetCrossSection(corsikaBeamId, targetId, ECoM);
[[maybe_unused]] auto elaCrossSectionCopy =
elaCrossSection; // ONLY TO AVOID COMPILER WARNING
@@ -152,7 +150,6 @@ namespace corsika::process::sibyll {
<< " IntLength: sibyll return (mb): " << productionCrossSection / 1_mbarn
<< endl;
weightedProdCrossSection += w[i] * productionCrossSection;
- avgTargetMassNumber += w[i] * numberOfNucleons;
}
cout << "Interaction: "
<< "IntLength: weighted CrossSection (mb): "
@@ -161,7 +158,7 @@ namespace corsika::process::sibyll {
// calculate interaction length in medium
//#warning check interaction length units
GrammageType const int_length =
- avgTargetMassNumber * units::constants::u / weightedProdCrossSection;
+ mediumComposition.GetAverageMassNumber() * units::constants::u / weightedProdCrossSection;
cout << "Interaction: "
<< "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
<< endl;
@@ -205,10 +202,8 @@ namespace corsika::process::sibyll {
// define target
// for Sibyll is always a single nucleon
- auto constexpr nucleon_mass =
- 0.5 * (particles::Proton::GetMass() + particles::Neutron::GetMass());
// FOR NOW: target is always at rest
- const auto eTargetLab = 0_GeV + nucleon_mass;
+ const auto eTargetLab = 0_GeV + constants::nucleonMass;
const auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargLab(eTargetLab, pTargetLab);
@@ -227,7 +222,7 @@ namespace corsika::process::sibyll {
// define target kinematics in lab frame
// define boost to and from CoM frame
// CoM frame definition in Sibyll projectile: +z
- COMBoost const boost(PprojLab, nucleon_mass);
+ COMBoost const boost(PprojLab, constants::nucleonMass);
// just for show:
// boost projecticle
@@ -268,11 +263,9 @@ namespace corsika::process::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
- const auto [sigProd, sigEla, nNuc] =
+ const auto [sigProd, sigEla ] =
GetCrossSection(corsikaBeamId, targetId, Ecm);
cross_section_of_components[i] = sigProd;
- [[maybe_unused]] int ideleteme =
- nNuc; // to avoid not used warning in array binding
[[maybe_unused]] auto sigElaCopy =
sigEla; // to avoid not used warning in array binding
}
@@ -289,7 +282,7 @@ namespace corsika::process::sibyll {
if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode);
if (targetCode == particles::Proton::GetCode()) targetSibCode = 1;
cout << "Interaction: sibyll code: " << targetSibCode << endl;
- if (targetSibCode > 18 || targetSibCode < 1)
+ if (targetSibCode > fMaxTargetMassNumber || targetSibCode < 1)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 or protons are "
"allowed.");
@@ -300,7 +293,10 @@ namespace corsika::process::sibyll {
cout << "Interaction: "
<< " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
<< " Ecm(GeV): " << Ecm / 1_GeV << endl;
- if (eProjectileLab < 8.5_GeV || Ecm < 10_GeV) {
+ if( Ecm > GetMaxEnergyCoM() )
+ throw std::runtime_error("Interaction::DoInteraction: CoM energy too high!");
+ // FR: removed eProjectileLab < 8.5_GeV ||
+ if ( Ecm < GetMinEnergyCoM() ) {
cout << "Interaction: "
<< " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << endl;
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 40463f6657c2b6e877be5a018b7cfcf67ae6f884..509f50a61861d30e0797e86d5190abcf469cb5b8 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -37,13 +37,20 @@ namespace corsika::process::sibyll {
void Init();
bool WasInitialized() { return fInitialized; }
- bool ValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
- using namespace corsika::units::si;
- return (10_GeV < ecm) && (ecm < 1_PeV);
+ bool IsValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
+ return (fMinEnergyCoM <= ecm) && (ecm <= fMaxEnergyCoM);
}
-
- std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType,
- int>
+ int GetMaxTargetMassNumber() { return fMaxTargetMassNumber; }
+ corsika::units::si::HEPEnergyType GetMinEnergyCoM() { return fMinEnergyCoM; }
+ corsika::units::si::HEPEnergyType GetMaxEnergyCoM() { return fMaxEnergyCoM; }
+ bool IsValidTarget(corsika::particles::Code TargetId)
+ {
+ return ( corsika::particles::GetNucleusA(TargetId) < fMaxTargetMassNumber )
+ && corsika::particles::IsNucleus( TargetId );
+ }
+
+
+ std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(const corsika::particles::Code BeamId,
const corsika::particles::Code TargetId,
const corsika::units::si::HEPEnergyType CoMenergy);
@@ -63,6 +70,11 @@ namespace corsika::process::sibyll {
corsika::environment::Environment const& fEnvironment;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+ const corsika::units::si::HEPEnergyType fMinEnergyCoM =
+ 10. * 1e9 * corsika::units::si::electronvolt;
+ const corsika::units::si::HEPEnergyType fMaxEnergyCoM =
+ 1.e6 * 1e9 * corsika::units::si::electronvolt;
+ const int fMaxTargetMassNumber = 18;
};
} // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index 6a37363235f66c455206e16c0134a9b623cf4ff4..0e65121d478942af665824f1348777a3764e21da 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -22,6 +22,8 @@
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
+#include <set>
+
using std::cout;
using std::endl;
using std::tuple;
@@ -41,6 +43,7 @@ namespace corsika::process::sibyll {
NuclearInteraction::~NuclearInteraction() {
cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount << endl;
+ fTargetComponentsIndex.clear();
}
void NuclearInteraction::Init() {
@@ -51,25 +54,136 @@ namespace corsika::process::sibyll {
// TODO: safe to run multiple initializations?
if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init();
+ // check compatibility of energy ranges, someone could try to use low-energy model..
+ if(!fHadronicInteraction.IsValidCoMEnergy(GetMinEnergyPerNucleonCoM())||
+ !fHadronicInteraction.IsValidCoMEnergy(GetMaxEnergyPerNucleonCoM()))
+ throw std::runtime_error("NuclearInteraction: hadronic interaction model incompatible!");
+
// 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;
+
+ auto& universe = *(fEnvironment.GetUniverse());
+
+ auto const allElementsInUniverse = std::invoke([&]() {
+ std::set<particles::Code> allElementsInUniverse;
+ auto collectElements = [&](auto& vtn) {
+ if (auto const mp = vtn.GetModelPropertiesPtr();
+ mp != nullptr) { // do not query Universe it self, it has no ModelProperties
+ auto const& comp = mp->GetNuclearComposition().GetComponents();
+ for (auto const c : comp)
+ allElementsInUniverse.insert(c);
+ // std::for_each(comp.cbegin(), comp.cend(),
+ // [&](particles::Code c) { allElementsInUniverse.insert(c); });
+ }
+ };
+ universe.walk(collectElements);
+ return allElementsInUniverse;
+ });
+
+
+ cout << "NuclearInteraction: initializing nuclear cross sections..." << endl;
+
+ // loop over target components, at most 4!!
+ int k =-1;
+ for(auto &ptarg: allElementsInUniverse){
+ ++k;
+ cout << "NuclearInteraction: init target component: " << ptarg << endl;
+ const int ib = GetNucleusA( ptarg );
+ if(!fHadronicInteraction.IsValidTarget( ptarg )){
+ cout << "NuclearInteraction::InitializeNuclearCrossSections: target nucleus? id="
+ << ptarg << endl;
+ throw std::runtime_error(" target can not be handled by hadronic interaction model! ");
+ }
+ fTargetComponentsIndex.insert( std::pair<Code,int>(ptarg, k) );
+ // loop over energies, fNEnBins log. energy bins
+ for(int i=0; i<GetNEnergyBins(); ++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 ] =
+ fHadronicInteraction.GetCrossSection(protonId, protonId, Ecm);
+ const double dsig = siginel / 1_mbarn;
+ const double dsigela = sigela / 1_mbarn;
+ // loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
+ for(int j=1; j<fMaxNucleusAProjectile; ++j){
+ const int jj = j+1;
+ double sig_out, dsig_out, sigqe_out, dsigqe_out;
+ 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 for " << fTargetComponentsIndex.size() << " components initialized!" << endl;
+ for(auto &ptarg: allElementsInUniverse){
+ cout << "cross section table: " << ptarg << endl;
+ PrintCrossSectionTable( ptarg );
+ }
+
}
- // TODO: remove number of nucleons, avg target mass is available in environment
+ void NuclearInteraction::PrintCrossSectionTable( corsika::particles::Code pCode)
+ {
+ using namespace corsika::particles;
+ const int k = fTargetComponentsIndex.at( pCode );
+ Code pNuclei [] = {Code::Helium, Code::Lithium7, Code::Oxygen,
+ Code::Neon, Code::Argon, Code::Iron};
+ cout << "en/A ";
+ for(auto &j: pNuclei)
+ cout << std::setw(9) << j;
+ cout << endl;
+
+ // loop over energy bins
+ for(int i=0; i<GetNEnergyBins(); ++i){
+ cout << " " << i << " ";
+ for(auto &n: pNuclei){
+ auto const j= GetNucleusA( n );
+ cout << " " << std::setprecision(5) << std::setw(8) << cnucsignuc_.sigma[j-1][k][i];
+ }
+ cout << endl;
+ }
+ }
+
+ units::si::CrossSectionType NuclearInteraction::ReadCrossSectionTable(const int ia, particles::Code pTarget, units::si::HEPEnergyType elabnuc)
+ {
+ using namespace corsika::particles;
+ using namespace units::si;
+ const int ib = fTargetComponentsIndex.at( pTarget )+1; // table index in fortran
+ auto const ECoMNuc = sqrt( 2. * corsika::units::constants::nucleonMass * elabnuc );
+ if( ECoMNuc < GetMinEnergyPerNucleonCoM() || ECoMNuc > GetMaxEnergyPerNucleonCoM() )
+ throw std::runtime_error("NuclearInteraction: energy outside tabulated range!");
+ const double e0 = elabnuc / 1_GeV;
+ double sig;
+ cout << "ReadCrossSectionTable: " << ia << " " << ib << " " << e0 << endl;
+ signuc2_(ia,ib,e0,sig);
+ cout << "ReadCrossSectionTable: sig=" << sig << endl;
+ return sig * 1_mbarn;
+ }
+
+ // TODO: remove elastic cross section?
template <>
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
NuclearInteraction::GetCrossSection(Particle& p, const particles::Code TargetId) {
using namespace units::si;
- double sigProd;
- auto const pCode = p.GetPID();
- if (pCode != particles::Code::Nucleus)
+ if ( p.GetPID() != particles::Code::Nucleus)
throw std::runtime_error(
"NuclearInteraction: GetCrossSection: particle not a nucleus!");
- auto const iBeam = p.GetNuclearA();
- HEPEnergyType LabEnergyPerNuc = p.GetEnergy() / iBeam;
- cout << "NuclearInteraction: GetCrossSection: called with: beamNuclA= " << iBeam
+ auto const iBeamA = p.GetNuclearA();
+ HEPEnergyType LabEnergyPerNuc = p.GetEnergy() / iBeamA;
+ cout << "NuclearInteraction: GetCrossSection: called with: beamNuclA= " << iBeamA
<< " TargetId= " << TargetId << " LabEnergyPerNuc= " << LabEnergyPerNuc / 1_GeV
<< endl;
@@ -77,38 +191,20 @@ namespace corsika::process::sibyll {
// TODO: for now assumes air with hard coded composition
// extend to arbitrary mixtures, requires smarter initialization
// get nuclib projectile code: nucleon number
- if (iBeam > 56 || iBeam < 2) {
+ if (iBeamA > GetMaxNucleusAProjectile() || iBeamA < 2) {
cout << "NuclearInteraction: beam nucleus outside allowed range for NUCLIB!" << endl
- << "A=" << iBeam << endl;
+ << "A=" << iBeamA << endl;
throw std::runtime_error(
"NuclearInteraction: GetCrossSection: beam nucleus outside allowed range for "
"NUCLIB!");
}
- const double dElabNuc = LabEnergyPerNuc / 1_GeV;
- // TODO: these limitations are still sibyll specific.
- // available target nuclei depends on the hadronic interaction model and the
- // initialization
- if (dElabNuc < 10.)
- throw std::runtime_error("NuclearInteraction: GetCrossSection: energy too low!");
-
- // TODO: these limitations are still sibyll specific.
- // available target nuclei depends on the hadronic interaction model and the
- // initialization
- if (particles::IsNucleus(TargetId)) {
- const int iTarget = particles::GetNucleusA(TargetId);
- if (iTarget > 18 || iTarget == 0)
- throw std::runtime_error(
- "Sibyll target outside range. Only nuclei with A<18 are allowed.");
- cout << "NuclearInteraction: calling signuc.." << endl;
- cout << "WARNING: using hard coded cross section for Nucleus-Air with "
- "SIBYLL! (fix me!)"
- << endl;
- // TODO: target id is not used because cross section is still hard coded and fixed
- // to air.
- signuc_(iBeam, dElabNuc, sigProd);
- cout << "cross section: " << sigProd << endl;
- return std::make_tuple(sigProd * 1_mbarn, 0_mbarn);
+ if (fHadronicInteraction.IsValidTarget(TargetId)) {
+ auto const sigProd = ReadCrossSectionTable( iBeamA, TargetId, LabEnergyPerNuc );
+ cout << "cross section (mb): " << sigProd / 1_mbarn << endl;
+ return std::make_tuple(sigProd, 0_mbarn);
+ } else {
+ throw std::runtime_error( "target outside range.");
}
return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mbarn,
std::numeric_limits<double>::infinity() * 1_mbarn);
@@ -137,14 +233,12 @@ namespace corsika::process::sibyll {
"projectiles should use NuclearStackExtension!");
// read from cross section code table
- const HEPMassType nucleon_mass =
- 0.5 * (particles::Proton::GetMass() + particles::Neutron::GetMass());
// FOR NOW: assume target is at rest
corsika::stack::MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
// total momentum and energy
- HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
+ HEPEnergyType Elab = p.GetEnergy() + constants::nucleonMass;
int const nuclA = p.GetNuclearA();
auto const ElabNuc = p.GetEnergy() / nuclA;
@@ -155,7 +249,7 @@ namespace corsika::process::sibyll {
// calculate cm. energy
const HEPEnergyType ECoM = sqrt(
(Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
- auto const ECoMNN = sqrt(2. * ElabNuc * nucleon_mass);
+ auto const ECoMNN = sqrt(2. * ElabNuc * constants::nucleonMass);
cout << "NuclearInteraction: LambdaInt: \n"
<< " input energy: " << Elab / 1_GeV << endl
<< " input energy CoM: " << ECoM / 1_GeV << endl
@@ -167,7 +261,7 @@ namespace corsika::process::sibyll {
// energy limits
// TODO: values depend on hadronic interaction model !! this is sibyll specific
- if (ElabNuc >= 8.5_GeV && ECoMNN >= 10_GeV) {
+ if (ElabNuc >= 8.5_GeV && ECoMNN >= fMinEnergyPerNucleonCoM && ECoMNN < fMaxEnergyPerNucleonCoM) {
// get target from environment
/*
@@ -264,7 +358,8 @@ namespace corsika::process::sibyll {
// projectile nucleon number
const int kAProj = p.GetNuclearA(); // GetNucleusA(ProjId);
- if (kAProj > 56) throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
+ if (kAProj > GetMaxNucleusAProjectile() )
+ throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
// kinematics
// define projectile nucleus
@@ -287,10 +382,8 @@ namespace corsika::process::sibyll {
// define target
// always a nucleon
- auto constexpr nucleon_mass =
- 0.5 * (particles::Proton::GetMass() + particles::Neutron::GetMass());
// target is always at rest
- const auto eTargetNucLab = 0_GeV + nucleon_mass;
+ const auto eTargetNucLab = 0_GeV + constants::nucleonMass;
const auto pTargetNucLab =
corsika::stack::MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargNucLab(eTargetNucLab, pTargetNucLab);
@@ -304,7 +397,7 @@ namespace corsika::process::sibyll {
HEPEnergyType EcmNN = PtotNN4.GetNorm();
cout << "NuclearInteraction: nuc-nuc cm energy: " << EcmNN / 1_GeV << endl;
- if (!fHadronicInteraction.ValidCoMEnergy(EcmNN)) {
+ if (!fHadronicInteraction.IsValidCoMEnergy(EcmNN)) {
cout << "NuclearInteraction: nuc-nuc. CoM energy too low for hadronic "
"interaction model!"
<< endl;
@@ -312,7 +405,7 @@ namespace corsika::process::sibyll {
}
// define boost to NUCLEON-NUCLEON frame
- COMBoost const boost(PprojNucLab, nucleon_mass);
+ COMBoost const boost(PprojNucLab, constants::nucleonMass);
// boost projecticle
auto const PprojNucCoM = boost.toCoM(PprojNucLab);
@@ -349,11 +442,10 @@ namespace corsika::process::sibyll {
auto const targetId = compVec[i];
cout << "target component: " << targetId << endl;
cout << "beam id: " << beamId << endl;
- const auto [sigProd, sigEla, nNuc] =
+ const auto [sigProd, sigEla] =
fHadronicInteraction.GetCrossSection(beamId, targetId, EcmNN);
cross_section_of_components[i] = sigProd;
[[maybe_unused]] auto sigElaCopy = sigEla; // ONLY TO AVOID COMPILER WARNINGS
- [[maybe_unused]] auto sigNucCopy = nNuc; // ONLY TO AVOID COMPILER WARNINGS
}
const auto targetCode = mediumComposition.SampleTarget(cross_section_of_components, fRNG);
@@ -367,10 +459,8 @@ namespace corsika::process::sibyll {
if (IsNucleus(targetCode)) kATarget = GetNucleusA(targetCode);
if (targetCode == particles::Proton::GetCode()) kATarget = 1;
cout << "NuclearInteraction: nuclib target code: " << kATarget << endl;
- if (kATarget > 18 || kATarget < 1)
- throw std::runtime_error(
- "Sibyll target outside range. Only nuclei with A<18 or protons are "
- "allowed.");
+ if(!fHadronicInteraction.IsValidTarget(targetCode))
+ throw std::runtime_error("target outside range. ");
// end of target sampling
// superposition
@@ -378,9 +468,8 @@ namespace corsika::process::sibyll {
// get nucleon-nucleon cross section
// (needed to determine number of nucleon-nucleon scatterings)
const auto protonId = particles::Proton::GetCode();
- const auto [prodCrossSection, elaCrossSection, dum] =
+ const auto [prodCrossSection, elaCrossSection] =
fHadronicInteraction.GetCrossSection(protonId, protonId, EcmNN);
- [[maybe_unused]] auto dumCopy = dum; // ONLY TO AVOID COMPILER WARNING
const double sigProd = prodCrossSection / 1_mbarn;
const double sigEla = elaCrossSection / 1_mbarn;
// sample number of interactions (only input variables, output in common cnucms)
@@ -414,7 +503,7 @@ namespace corsika::process::sibyll {
fragm_(kATarget, kAProj, nIntProj, impactPar, nFragments, AFragments);
// this should not occur but well :)
- if (nFragments > 60)
+ if (nFragments > GetMaxNFragments())
throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
cout << "number of fragments: " << nFragments << endl;
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
index cf1f6593330fddc72e867dd7d40387747d6ee9be..ad1d0bc7e06e064d859f0817d0128d6b1d73170e 100644
--- a/Processes/Sibyll/NuclearInteraction.h
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -40,7 +40,14 @@ namespace corsika::process::sibyll {
corsika::process::sibyll::Interaction& hadint);
~NuclearInteraction();
void Init();
-
+ void InitializeNuclearCrossSections();
+ void PrintCrossSectionTable(corsika::particles::Code);
+ corsika::units::si::CrossSectionType ReadCrossSectionTable(const int, corsika::particles::Code, corsika::units::si::HEPEnergyType);
+ corsika::units::si::HEPEnergyType GetMinEnergyPerNucleonCoM() { return fMinEnergyPerNucleonCoM; }
+ corsika::units::si::HEPEnergyType GetMaxEnergyPerNucleonCoM() { return fMaxEnergyPerNucleonCoM; }
+ int GetMaxNucleusAProjectile() {return fMaxNucleusAProjectile; }
+ int GetMaxNFragments() { return fMaxNFragments; }
+ int GetNEnergyBins() { return fNEnBins; }
template <typename Particle>
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(Particle& p, const corsika::particles::Code TargetId);
@@ -54,8 +61,19 @@ namespace corsika::process::sibyll {
private:
corsika::environment::Environment const& fEnvironment;
corsika::process::sibyll::Interaction& fHadronicInteraction;
+ std::map<corsika::particles::Code,int> fTargetComponentsIndex;
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+ const int fNSample = 500; // number of samples in MC estimation of cross section
+ const int fMaxNucleusAProjectile = 56;
+ const int fNEnBins = 6;
+ const int fMaxNFragments = 60;
+ // energy limits defined by table used for cross section in signuc.f
+ // 10**1 GeV to 10**6 GeV
+ const corsika::units::si::HEPEnergyType fMinEnergyPerNucleonCoM =
+ 10. * 1e9 * corsika::units::si::electronvolt;
+ const corsika::units::si::HEPEnergyType fMaxEnergyPerNucleonCoM =
+ 1.e6 * 1e9 * corsika::units::si::electronvolt;
};
} // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/nuclib.h b/Processes/Sibyll/nuclib.h
index e4b89b52bed8bab3f1898183b2398c5ef3c2d3b3..ddd1856cd2eea972e257ba2958f6705ae15bcf96 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 b21bce8cbb4893c90b1bfca4f5c1a0ca617fcb4f..b88f2ae0bd37f76bcb693692ac0fea0911b053bf 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
+