IAP GITLAB

Skip to content
Snippets Groups Projects
Commit cd8e728b authored by Felix Riehn's avatar Felix Riehn
Browse files

Merge branch '132-generalize-nuclear-cross-sections' into 'master' 

Resolve "generalize nuclear cross sections"

Closes #132

See merge request AirShowerPhysics/corsika!84
parents 92dbf62a 60abf6dc
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 276 additions and 87 deletions
Environment/NuclearComposition.h
+ 3
− 3
View file @ cd8e728b
...@@ -35,15 +35,15 @@ namespace corsika::environment { ...@@ -35,15 +35,15 @@ namespace corsika::environment {
public: public:
using value_type = double; using value_type = double;
using iterator_category = std::input_iterator_tag; using iterator_category = std::input_iterator_tag;
using pointer = double*; using pointer = value_type*;
using reference = double&; using reference = value_type&;
using difference_type = ptrdiff_t; using difference_type = ptrdiff_t;
WeightProviderIterator(AConstIterator a, BConstIterator b) WeightProviderIterator(AConstIterator a, BConstIterator b)
: fAIter(a) : fAIter(a)
, fBIter(b) {} , fBIter(b) {}
double operator*() const { return ((*fAIter) * (*fBIter)).magnitude(); } value_type operator*() const { return ((*fAIter) * (*fBIter)).magnitude(); }
WeightProviderIterator& operator++() { // prefix ++ WeightProviderIterator& operator++() { // prefix ++
++fAIter; ++fAIter;
......
Environment/VolumeTreeNode.h
+ 20
− 0
View file @ cd8e728b
...@@ -66,6 +66,24 @@ namespace corsika::environment { ...@@ -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) { void AddChild(VTNUPtr pChild) {
pChild->fParentNode = this; pChild->fParentNode = this;
fChildNodes.push_back(std::move(pChild)); fChildNodes.push_back(std::move(pChild));
...@@ -88,6 +106,8 @@ namespace corsika::environment { ...@@ -88,6 +106,8 @@ namespace corsika::environment {
auto const& GetModelProperties() const { return *fModelProperties; } auto const& GetModelProperties() const { return *fModelProperties; }
IMPSharedPtr GetModelPropertiesPtr() const { return fModelProperties; }
template <typename TModelProperties, typename... Args> template <typename TModelProperties, typename... Args>
auto SetModelProperties(Args&&... args) { auto SetModelProperties(Args&&... args) {
static_assert(std::is_base_of_v<IModelProperties, TModelProperties>, static_assert(std::is_base_of_v<IModelProperties, TModelProperties>,
......
Processes/Sibyll/Interaction.cc
+ 20
− 24
View file @ cd8e728b
...@@ -53,7 +53,7 @@ namespace corsika::process::sibyll { ...@@ -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, Interaction::GetCrossSection(const particles::Code BeamId,
const particles::Code TargetId, const particles::Code TargetId,
const units::si::HEPEnergyType CoMenergy) { const units::si::HEPEnergyType CoMenergy) {
...@@ -61,23 +61,24 @@ namespace corsika::process::sibyll { ...@@ -61,23 +61,24 @@ namespace corsika::process::sibyll {
double sigProd, sigEla, dummy, dum1, dum3, dum4; double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3]; double dumdif[3];
const int iBeam = process::sibyll::GetSibyllXSCode(BeamId); 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; const double dEcm = CoMenergy / 1_GeV;
int iTarget = -1;
if (particles::IsNucleus(TargetId)) { if (particles::IsNucleus(TargetId)) {
iTarget = particles::GetNucleusA(TargetId); const int iTarget = particles::GetNucleusA(TargetId);
if (iTarget > 18 || iTarget == 0) if (iTarget > fMaxTargetMassNumber || iTarget == 0)
throw std::runtime_error( throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 are allowed."); "Sibyll target outside range. Only nuclei with A<18 are allowed.");
sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla); sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
} else if (TargetId == particles::Proton::GetCode()) { } else if (TargetId == particles::Proton::GetCode()) {
sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4); sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
iTarget = 1;
} else { } else {
// no interaction in sibyll possible, return infinite cross section? or throw? // no interaction in sibyll possible, return infinite cross section? or throw?
sigProd = std::numeric_limits<double>::infinity(); sigProd = std::numeric_limits<double>::infinity();
sigEla = 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 <> template <>
...@@ -97,14 +98,11 @@ namespace corsika::process::sibyll { ...@@ -97,14 +98,11 @@ namespace corsika::process::sibyll {
// read from cross section code table // read from cross section code table
const bool kInteraction = process::sibyll::CanInteract(corsikaBeamId); 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 // FOR NOW: assume target is at rest
MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV}); MomentumVector pTarget(rootCS, {0_GeV, 0_GeV, 0_GeV});
// total momentum and energy // 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}); MomentumVector pTotLab(rootCS, {0_GeV, 0_GeV, 0_GeV});
pTotLab += p.GetMomentum(); pTotLab += p.GetMomentum();
pTotLab += pTarget; pTotLab += pTarget;
...@@ -119,7 +117,8 @@ namespace corsika::process::sibyll { ...@@ -119,7 +117,8 @@ namespace corsika::process::sibyll {
<< " beam pid:" << p.GetPID() << endl; << " beam pid:" << p.GetPID() << endl;
// TODO: move limits into variables // 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 // get target from environment
/* /*
...@@ -134,7 +133,6 @@ namespace corsika::process::sibyll { ...@@ -134,7 +133,6 @@ namespace corsika::process::sibyll {
// determine average interaction length // determine average interaction length
// weighted sum // weighted sum
int i = -1; int i = -1;
double avgTargetMassNumber = 0.;
si::CrossSectionType weightedProdCrossSection = 0_mbarn; si::CrossSectionType weightedProdCrossSection = 0_mbarn;
// get weights of components from environment/medium // get weights of components from environment/medium
const auto w = mediumComposition.GetFractions(); const auto w = mediumComposition.GetFractions();
...@@ -143,7 +141,7 @@ namespace corsika::process::sibyll { ...@@ -143,7 +141,7 @@ namespace corsika::process::sibyll {
i++; i++;
cout << "Interaction: get interaction length for target: " << targetId << endl; cout << "Interaction: get interaction length for target: " << targetId << endl;
auto const [productionCrossSection, elaCrossSection, numberOfNucleons] = auto const [productionCrossSection, elaCrossSection ] =
GetCrossSection(corsikaBeamId, targetId, ECoM); GetCrossSection(corsikaBeamId, targetId, ECoM);
[[maybe_unused]] auto elaCrossSectionCopy = [[maybe_unused]] auto elaCrossSectionCopy =
elaCrossSection; // ONLY TO AVOID COMPILER WARNING elaCrossSection; // ONLY TO AVOID COMPILER WARNING
...@@ -152,7 +150,6 @@ namespace corsika::process::sibyll { ...@@ -152,7 +150,6 @@ namespace corsika::process::sibyll {
<< " IntLength: sibyll return (mb): " << productionCrossSection / 1_mbarn << " IntLength: sibyll return (mb): " << productionCrossSection / 1_mbarn
<< endl; << endl;
weightedProdCrossSection += w[i] * productionCrossSection; weightedProdCrossSection += w[i] * productionCrossSection;
avgTargetMassNumber += w[i] * numberOfNucleons;
} }
cout << "Interaction: " cout << "Interaction: "
<< "IntLength: weighted CrossSection (mb): " << "IntLength: weighted CrossSection (mb): "
...@@ -161,7 +158,7 @@ namespace corsika::process::sibyll { ...@@ -161,7 +158,7 @@ namespace corsika::process::sibyll {
// calculate interaction length in medium // calculate interaction length in medium
//#warning check interaction length units //#warning check interaction length units
GrammageType const int_length = GrammageType const int_length =
avgTargetMassNumber * units::constants::u / weightedProdCrossSection; mediumComposition.GetAverageMassNumber() * units::constants::u / weightedProdCrossSection;
cout << "Interaction: " cout << "Interaction: "
<< "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm << "interaction length (g/cm2): " << int_length / (0.001_kg) * 1_cm * 1_cm
<< endl; << endl;
...@@ -205,10 +202,8 @@ namespace corsika::process::sibyll { ...@@ -205,10 +202,8 @@ namespace corsika::process::sibyll {
// define target // define target
// for Sibyll is always a single nucleon // 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 // 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 auto pTargetLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargLab(eTargetLab, pTargetLab); const FourVector PtargLab(eTargetLab, pTargetLab);
...@@ -227,7 +222,7 @@ namespace corsika::process::sibyll { ...@@ -227,7 +222,7 @@ namespace corsika::process::sibyll {
// define target kinematics in lab frame // define target kinematics in lab frame
// define boost to and from CoM frame // define boost to and from CoM frame
// CoM frame definition in Sibyll projectile: +z // CoM frame definition in Sibyll projectile: +z
COMBoost const boost(PprojLab, nucleon_mass); COMBoost const boost(PprojLab, constants::nucleonMass);
// just for show: // just for show:
// boost projecticle // boost projecticle
...@@ -268,11 +263,9 @@ namespace corsika::process::sibyll { ...@@ -268,11 +263,9 @@ namespace corsika::process::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) { for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i]; auto const targetId = compVec[i];
const auto [sigProd, sigEla, nNuc] = const auto [sigProd, sigEla ] =
GetCrossSection(corsikaBeamId, targetId, Ecm); GetCrossSection(corsikaBeamId, targetId, Ecm);
cross_section_of_components[i] = sigProd; cross_section_of_components[i] = sigProd;
[[maybe_unused]] int ideleteme =
nNuc; // to avoid not used warning in array binding
[[maybe_unused]] auto sigElaCopy = [[maybe_unused]] auto sigElaCopy =
sigEla; // to avoid not used warning in array binding sigEla; // to avoid not used warning in array binding
} }
...@@ -289,7 +282,7 @@ namespace corsika::process::sibyll { ...@@ -289,7 +282,7 @@ namespace corsika::process::sibyll {
if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode); if (IsNucleus(targetCode)) targetSibCode = GetNucleusA(targetCode);
if (targetCode == particles::Proton::GetCode()) targetSibCode = 1; if (targetCode == particles::Proton::GetCode()) targetSibCode = 1;
cout << "Interaction: sibyll code: " << targetSibCode << endl; cout << "Interaction: sibyll code: " << targetSibCode << endl;
if (targetSibCode > 18 || targetSibCode < 1) if (targetSibCode > fMaxTargetMassNumber || targetSibCode < 1)
throw std::runtime_error( throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 or protons are " "Sibyll target outside range. Only nuclei with A<18 or protons are "
"allowed."); "allowed.");
...@@ -300,7 +293,10 @@ namespace corsika::process::sibyll { ...@@ -300,7 +293,10 @@ namespace corsika::process::sibyll {
cout << "Interaction: " cout << "Interaction: "
<< " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV << " DoInteraction: E(GeV):" << eProjectileLab / 1_GeV
<< " Ecm(GeV): " << Ecm / 1_GeV << endl; << " 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: " cout << "Interaction: "
<< " DoInteraction: should have dropped particle.. " << " DoInteraction: should have dropped particle.. "
<< "THIS IS AN ERROR" << endl; << "THIS IS AN ERROR" << endl;
......
Processes/Sibyll/Interaction.h
+ 18
− 6
View file @ cd8e728b
...@@ -37,13 +37,20 @@ namespace corsika::process::sibyll { ...@@ -37,13 +37,20 @@ namespace corsika::process::sibyll {
void Init(); void Init();
bool WasInitialized() { return fInitialized; } bool WasInitialized() { return fInitialized; }
bool ValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) { bool IsValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) {
using namespace corsika::units::si; return (fMinEnergyCoM <= ecm) && (ecm <= fMaxEnergyCoM);
return (10_GeV < ecm) && (ecm < 1_PeV);
} }
int GetMaxTargetMassNumber() { return fMaxTargetMassNumber; }
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType, corsika::units::si::HEPEnergyType GetMinEnergyCoM() { return fMinEnergyCoM; }
int> 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, GetCrossSection(const corsika::particles::Code BeamId,
const corsika::particles::Code TargetId, const corsika::particles::Code TargetId,
const corsika::units::si::HEPEnergyType CoMenergy); const corsika::units::si::HEPEnergyType CoMenergy);
...@@ -63,6 +70,11 @@ namespace corsika::process::sibyll { ...@@ -63,6 +70,11 @@ namespace corsika::process::sibyll {
corsika::environment::Environment const& fEnvironment; corsika::environment::Environment const& fEnvironment;
corsika::random::RNG& fRNG = corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm"); 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 } // namespace corsika::process::sibyll
......
Processes/Sibyll/NuclearInteraction.cc
+ 142
− 53
View file @ cd8e728b
...@@ -22,6 +22,8 @@ ...@@ -22,6 +22,8 @@
#include <corsika/setup/SetupStack.h> #include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h> #include <corsika/setup/SetupTrajectory.h>
#include <set>
using std::cout; using std::cout;
using std::endl; using std::endl;
using std::tuple; using std::tuple;
...@@ -41,6 +43,7 @@ namespace corsika::process::sibyll { ...@@ -41,6 +43,7 @@ namespace corsika::process::sibyll {
NuclearInteraction::~NuclearInteraction() { NuclearInteraction::~NuclearInteraction() {
cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount << endl; cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount << endl;
fTargetComponentsIndex.clear();
} }
void NuclearInteraction::Init() { void NuclearInteraction::Init() {
...@@ -51,25 +54,136 @@ namespace corsika::process::sibyll { ...@@ -51,25 +54,136 @@ namespace corsika::process::sibyll {
// TODO: safe to run multiple initializations? // TODO: safe to run multiple initializations?
if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init(); 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 // initialize nuclib
// TODO: make sure this does not overlap with sibyll // TODO: make sure this does not overlap with sibyll
nuc_nuc_ini_(); 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 <> template <>
tuple<units::si::CrossSectionType, units::si::CrossSectionType> tuple<units::si::CrossSectionType, units::si::CrossSectionType>
NuclearInteraction::GetCrossSection(Particle& p, const particles::Code TargetId) { NuclearInteraction::GetCrossSection(Particle& p, const particles::Code TargetId) {
using namespace units::si; using namespace units::si;
double sigProd; if ( p.GetPID() != particles::Code::Nucleus)
auto const pCode = p.GetPID();
if (pCode != particles::Code::Nucleus)
throw std::runtime_error( throw std::runtime_error(
"NuclearInteraction: GetCrossSection: particle not a nucleus!"); "NuclearInteraction: GetCrossSection: particle not a nucleus!");
auto const iBeam = p.GetNuclearA(); auto const iBeamA = p.GetNuclearA();
HEPEnergyType LabEnergyPerNuc = p.GetEnergy() / iBeam; HEPEnergyType LabEnergyPerNuc = p.GetEnergy() / iBeamA;
cout << "NuclearInteraction: GetCrossSection: called with: beamNuclA= " << iBeam cout << "NuclearInteraction: GetCrossSection: called with: beamNuclA= " << iBeamA
<< " TargetId= " << TargetId << " LabEnergyPerNuc= " << LabEnergyPerNuc / 1_GeV << " TargetId= " << TargetId << " LabEnergyPerNuc= " << LabEnergyPerNuc / 1_GeV
<< endl; << endl;
...@@ -77,38 +191,20 @@ namespace corsika::process::sibyll { ...@@ -77,38 +191,20 @@ namespace corsika::process::sibyll {
// TODO: for now assumes air with hard coded composition // TODO: for now assumes air with hard coded composition
// extend to arbitrary mixtures, requires smarter initialization // extend to arbitrary mixtures, requires smarter initialization
// get nuclib projectile code: nucleon number // 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 cout << "NuclearInteraction: beam nucleus outside allowed range for NUCLIB!" << endl
<< "A=" << iBeam << endl; << "A=" << iBeamA << endl;
throw std::runtime_error( throw std::runtime_error(
"NuclearInteraction: GetCrossSection: beam nucleus outside allowed range for " "NuclearInteraction: GetCrossSection: beam nucleus outside allowed range for "
"NUCLIB!"); "NUCLIB!");
} }
const double dElabNuc = LabEnergyPerNuc / 1_GeV; if (fHadronicInteraction.IsValidTarget(TargetId)) {
// TODO: these limitations are still sibyll specific. auto const sigProd = ReadCrossSectionTable( iBeamA, TargetId, LabEnergyPerNuc );
// available target nuclei depends on the hadronic interaction model and the cout << "cross section (mb): " << sigProd / 1_mbarn << endl;
// initialization return std::make_tuple(sigProd, 0_mbarn);
if (dElabNuc < 10.) } else {
throw std::runtime_error("NuclearInteraction: GetCrossSection: energy too low!"); throw std::runtime_error( "target outside range.");
// 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);
} }
return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mbarn, return std::make_tuple(std::numeric_limits<double>::infinity() * 1_mbarn,
std::numeric_limits<double>::infinity() * 1_mbarn); std::numeric_limits<double>::infinity() * 1_mbarn);
...@@ -137,14 +233,12 @@ namespace corsika::process::sibyll { ...@@ -137,14 +233,12 @@ namespace corsika::process::sibyll {
"projectiles should use NuclearStackExtension!"); "projectiles should use NuclearStackExtension!");
// read from cross section code table // 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 // FOR NOW: assume target is at rest
corsika::stack::MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV}); corsika::stack::MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
// total momentum and energy // total momentum and energy
HEPEnergyType Elab = p.GetEnergy() + nucleon_mass; HEPEnergyType Elab = p.GetEnergy() + constants::nucleonMass;
int const nuclA = p.GetNuclearA(); int const nuclA = p.GetNuclearA();
auto const ElabNuc = p.GetEnergy() / nuclA; auto const ElabNuc = p.GetEnergy() / nuclA;
...@@ -155,7 +249,7 @@ namespace corsika::process::sibyll { ...@@ -155,7 +249,7 @@ namespace corsika::process::sibyll {
// calculate cm. energy // calculate cm. energy
const HEPEnergyType ECoM = sqrt( const HEPEnergyType ECoM = sqrt(
(Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy (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" cout << "NuclearInteraction: LambdaInt: \n"
<< " input energy: " << Elab / 1_GeV << endl << " input energy: " << Elab / 1_GeV << endl
<< " input energy CoM: " << ECoM / 1_GeV << endl << " input energy CoM: " << ECoM / 1_GeV << endl
...@@ -167,7 +261,7 @@ namespace corsika::process::sibyll { ...@@ -167,7 +261,7 @@ namespace corsika::process::sibyll {
// energy limits // energy limits
// TODO: values depend on hadronic interaction model !! this is sibyll specific // 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 // get target from environment
/* /*
...@@ -264,7 +358,8 @@ namespace corsika::process::sibyll { ...@@ -264,7 +358,8 @@ namespace corsika::process::sibyll {
// projectile nucleon number // projectile nucleon number
const int kAProj = p.GetNuclearA(); // GetNucleusA(ProjId); 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 // kinematics
// define projectile nucleus // define projectile nucleus
...@@ -287,10 +382,8 @@ namespace corsika::process::sibyll { ...@@ -287,10 +382,8 @@ namespace corsika::process::sibyll {
// define target // define target
// always a nucleon // always a nucleon
auto constexpr nucleon_mass =
0.5 * (particles::Proton::GetMass() + particles::Neutron::GetMass());
// target is always at rest // target is always at rest
const auto eTargetNucLab = 0_GeV + nucleon_mass; const auto eTargetNucLab = 0_GeV + constants::nucleonMass;
const auto pTargetNucLab = const auto pTargetNucLab =
corsika::stack::MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV); corsika::stack::MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargNucLab(eTargetNucLab, pTargetNucLab); const FourVector PtargNucLab(eTargetNucLab, pTargetNucLab);
...@@ -304,7 +397,7 @@ namespace corsika::process::sibyll { ...@@ -304,7 +397,7 @@ namespace corsika::process::sibyll {
HEPEnergyType EcmNN = PtotNN4.GetNorm(); HEPEnergyType EcmNN = PtotNN4.GetNorm();
cout << "NuclearInteraction: nuc-nuc cm energy: " << EcmNN / 1_GeV << endl; 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 " cout << "NuclearInteraction: nuc-nuc. CoM energy too low for hadronic "
"interaction model!" "interaction model!"
<< endl; << endl;
...@@ -312,7 +405,7 @@ namespace corsika::process::sibyll { ...@@ -312,7 +405,7 @@ namespace corsika::process::sibyll {
} }
// define boost to NUCLEON-NUCLEON frame // define boost to NUCLEON-NUCLEON frame
COMBoost const boost(PprojNucLab, nucleon_mass); COMBoost const boost(PprojNucLab, constants::nucleonMass);
// boost projecticle // boost projecticle
auto const PprojNucCoM = boost.toCoM(PprojNucLab); auto const PprojNucCoM = boost.toCoM(PprojNucLab);
...@@ -349,11 +442,10 @@ namespace corsika::process::sibyll { ...@@ -349,11 +442,10 @@ namespace corsika::process::sibyll {
auto const targetId = compVec[i]; auto const targetId = compVec[i];
cout << "target component: " << targetId << endl; cout << "target component: " << targetId << endl;
cout << "beam id: " << beamId << endl; cout << "beam id: " << beamId << endl;
const auto [sigProd, sigEla, nNuc] = const auto [sigProd, sigEla] =
fHadronicInteraction.GetCrossSection(beamId, targetId, EcmNN); fHadronicInteraction.GetCrossSection(beamId, targetId, EcmNN);
cross_section_of_components[i] = sigProd; cross_section_of_components[i] = sigProd;
[[maybe_unused]] auto sigElaCopy = sigEla; // ONLY TO AVOID COMPILER WARNINGS [[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); const auto targetCode = mediumComposition.SampleTarget(cross_section_of_components, fRNG);
...@@ -367,10 +459,8 @@ namespace corsika::process::sibyll { ...@@ -367,10 +459,8 @@ namespace corsika::process::sibyll {
if (IsNucleus(targetCode)) kATarget = GetNucleusA(targetCode); if (IsNucleus(targetCode)) kATarget = GetNucleusA(targetCode);
if (targetCode == particles::Proton::GetCode()) kATarget = 1; if (targetCode == particles::Proton::GetCode()) kATarget = 1;
cout << "NuclearInteraction: nuclib target code: " << kATarget << endl; cout << "NuclearInteraction: nuclib target code: " << kATarget << endl;
if (kATarget > 18 || kATarget < 1) if(!fHadronicInteraction.IsValidTarget(targetCode))
throw std::runtime_error( throw std::runtime_error("target outside range. ");
"Sibyll target outside range. Only nuclei with A<18 or protons are "
"allowed.");
// end of target sampling // end of target sampling
// superposition // superposition
...@@ -378,9 +468,8 @@ namespace corsika::process::sibyll { ...@@ -378,9 +468,8 @@ namespace corsika::process::sibyll {
// get nucleon-nucleon cross section // get nucleon-nucleon cross section
// (needed to determine number of nucleon-nucleon scatterings) // (needed to determine number of nucleon-nucleon scatterings)
const auto protonId = particles::Proton::GetCode(); const auto protonId = particles::Proton::GetCode();
const auto [prodCrossSection, elaCrossSection, dum] = const auto [prodCrossSection, elaCrossSection] =
fHadronicInteraction.GetCrossSection(protonId, protonId, EcmNN); fHadronicInteraction.GetCrossSection(protonId, protonId, EcmNN);
[[maybe_unused]] auto dumCopy = dum; // ONLY TO AVOID COMPILER WARNING
const double sigProd = prodCrossSection / 1_mbarn; const double sigProd = prodCrossSection / 1_mbarn;
const double sigEla = elaCrossSection / 1_mbarn; const double sigEla = elaCrossSection / 1_mbarn;
// sample number of interactions (only input variables, output in common cnucms) // sample number of interactions (only input variables, output in common cnucms)
...@@ -414,7 +503,7 @@ namespace corsika::process::sibyll { ...@@ -414,7 +503,7 @@ namespace corsika::process::sibyll {
fragm_(kATarget, kAProj, nIntProj, impactPar, nFragments, AFragments); fragm_(kATarget, kAProj, nIntProj, impactPar, nFragments, AFragments);
// this should not occur but well :) // this should not occur but well :)
if (nFragments > 60) if (nFragments > GetMaxNFragments())
throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!"); throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
cout << "number of fragments: " << nFragments << endl; cout << "number of fragments: " << nFragments << endl;
......
Processes/Sibyll/NuclearInteraction.h
+ 19
− 1
View file @ cd8e728b
...@@ -40,7 +40,14 @@ namespace corsika::process::sibyll { ...@@ -40,7 +40,14 @@ namespace corsika::process::sibyll {
corsika::process::sibyll::Interaction& hadint); corsika::process::sibyll::Interaction& hadint);
~NuclearInteraction(); ~NuclearInteraction();
void Init(); 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> template <typename Particle>
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType> std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
GetCrossSection(Particle& p, const corsika::particles::Code TargetId); GetCrossSection(Particle& p, const corsika::particles::Code TargetId);
...@@ -54,8 +61,19 @@ namespace corsika::process::sibyll { ...@@ -54,8 +61,19 @@ namespace corsika::process::sibyll {
private: private:
corsika::environment::Environment const& fEnvironment; corsika::environment::Environment const& fEnvironment;
corsika::process::sibyll::Interaction& fHadronicInteraction; corsika::process::sibyll::Interaction& fHadronicInteraction;
std::map<corsika::particles::Code,int> fTargetComponentsIndex;
corsika::random::RNG& fRNG = corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm"); 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 } // namespace corsika::process::sibyll
......
Processes/Sibyll/nuclib.h
+ 11
− 0
View file @ cd8e728b
...@@ -34,6 +34,13 @@ extern struct { ...@@ -34,6 +34,13 @@ extern struct {
*/ */
extern struct { double ppp[60][3]; } fragments_; 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 // NUCLIB
// subroutine to initiate nuclib // subroutine to initiate nuclib
...@@ -46,5 +53,9 @@ void int_nuc_(const int&, const int&, const double&, const double&); ...@@ -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 fragm_(const int&, const int&, const int&, const double&, int&, int*);
void signuc_(const int&, const double&, double&); 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 #endif
Processes/Sibyll/signuc.f
+ 43
− 0
View file @ cd8e728b
...@@ -274,3 +274,46 @@ C ADD INELASTIC AND QUASI-ELASTIC CROSS-SECTIONS ...@@ -274,3 +274,46 @@ C ADD INELASTIC AND QUASI-ELASTIC CROSS-SECTIONS
RETURN RETURN
END 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
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