IAP GITLAB

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

Merge branch 'sibyll_small_cleanup_2' into 'master' 

Sibyll small cleanup

See merge request AirShowerPhysics/corsika!173
parents 4a991297 1fb97319
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Processes/Sibyll/CMakeLists.txt
+ 2
− 2
View file @ 91cb1961
......@@ -57,8 +57,8 @@ add_custom_command (
COMMAND ${CMAKE_COMMAND} -E create_symlink ${PROJECT_BINARY_DIR}/include/corsika/process/sibyll/Generated.inc ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc
COMMENT "Generate link in source-dir: ${CMAKE_CURRENT_SOURCE_DIR}/Generated.inc"
)
add_custom_target (SourceDirLink2 DEPENDS ${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc)
add_dependencies (ProcessSibyll SourceDirLink2)
add_custom_target (SourceDirLinkSib DEPENDS ${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc)
add_dependencies (ProcessSibyll SourceDirLinkSib)
# .....................................................
......
Processes/Sibyll/Interaction.cc
+ 10
− 10
View file @ 91cb1961
......@@ -37,7 +37,7 @@ namespace corsika::process::sibyll {
Interaction::Interaction() {}
Interaction::~Interaction() {
cout << "Sibyll::Interaction n=" << fCount << " Nnuc=" << fNucCount << endl;
cout << "Sibyll::Interaction n=" << count_ << " Nnuc=" << nucCount_ << endl;
}
void Interaction::Init() {
......@@ -45,9 +45,9 @@ namespace corsika::process::sibyll {
using random::RNGManager;
// initialize Sibyll
if (!fInitialized) {
if (!initialized_) {
sibyll_ini_();
fInitialized = true;
initialized_ = true;
}
}
......@@ -94,7 +94,7 @@ namespace corsika::process::sibyll {
const double dEcm = CoMenergy / 1_GeV;
if (particles::IsNucleus(TargetId)) {
const int iTarget = particles::GetNucleusA(TargetId);
if (iTarget > fMaxTargetMassNumber || iTarget == 0)
if (iTarget > maxTargetMassNumber_ || 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);
......@@ -281,7 +281,7 @@ namespace corsika::process::sibyll {
}
const auto targetCode =
mediumComposition.SampleTarget(cross_section_of_components, fRNG);
mediumComposition.SampleTarget(cross_section_of_components, RNG_);
cout << "Interaction: target selected: " << targetCode << endl;
/*
FOR NOW: allow nuclei with A<18 or protons only.
......@@ -292,7 +292,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 > fMaxTargetMassNumber || targetSibCode < 1)
if (targetSibCode > maxTargetMassNumber_ || targetSibCode < 1)
throw std::runtime_error(
"Sibyll target outside range. Only nuclei with A<18 or protons are "
"allowed.");
......@@ -312,17 +312,17 @@ namespace corsika::process::sibyll {
<< "THIS IS AN ERROR" << endl;
throw std::runtime_error("energy too low for SIBYLL");
} else {
fCount++;
count_++;
// Sibyll does not know about units..
const double sqs = Ecm / 1_GeV;
// running sibyll, filling stack
sibyll_(kBeam, targetSibCode, sqs);
if (fInternalDecays) {
if (internalDecays_) {
// particles that decay internally will never appear on the corsika stack
// switch on all decays except for the particles we want to take part in the
// tracking
SetAllUnstable();
SetStable(fTrackedParticles);
SetStable(trackedParticles_);
decsib_();
// reset
SetAllStable();
......@@ -330,7 +330,7 @@ namespace corsika::process::sibyll {
// print final state
int print_unit = 6;
sib_list_(print_unit);
fNucCount += get_nwounded() - 1;
nucCount_ += get_nwounded() - 1;
// add particles from sibyll to stack
// link to sibyll stack
......
Processes/Sibyll/Interaction.h
+ 15
− 15
View file @ 91cb1961
......@@ -21,9 +21,9 @@ namespace corsika::process::sibyll {
class Interaction : public corsika::process::InteractionProcess<Interaction> {
int fCount = 0;
int fNucCount = 0;
bool fInitialized = false;
int count_ = 0;
int nucCount_ = 0;
bool initialized_ = false;
public:
Interaction();
......@@ -39,15 +39,15 @@ namespace corsika::process::sibyll {
void SetAllUnstable();
void SetAllStable();
bool WasInitialized() { return fInitialized; }
bool WasInitialized() { return initialized_; }
bool IsValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) const {
return (fMinEnergyCoM <= ecm) && (ecm <= fMaxEnergyCoM);
return (minEnergyCoM_ <= ecm) && (ecm <= maxEnergyCoM_);
}
int GetMaxTargetMassNumber() const { return fMaxTargetMassNumber; }
corsika::units::si::HEPEnergyType GetMinEnergyCoM() const { return fMinEnergyCoM; }
corsika::units::si::HEPEnergyType GetMaxEnergyCoM() const { return fMaxEnergyCoM; }
int GetMaxTargetMassNumber() const { return maxTargetMassNumber_; }
corsika::units::si::HEPEnergyType GetMinEnergyCoM() const { return minEnergyCoM_; }
corsika::units::si::HEPEnergyType GetMaxEnergyCoM() const { return maxEnergyCoM_; }
bool IsValidTarget(corsika::particles::Code TargetId) const {
return (corsika::particles::GetNucleusA(TargetId) < fMaxTargetMassNumber) &&
return (corsika::particles::GetNucleusA(TargetId) < maxTargetMassNumber_) &&
corsika::particles::IsNucleus(TargetId);
}
......@@ -67,10 +67,10 @@ namespace corsika::process::sibyll {
corsika::process::EProcessReturn DoInteraction(TProjectile&);
private:
corsika::random::RNG& fRNG =
corsika::random::RNG& RNG_ =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
// FOR NOW keep trackedParticles private, could be configurable
std::vector<particles::Code> const fTrackedParticles = {
std::vector<particles::Code> const trackedParticles_ = {
particles::Code::PiPlus, particles::Code::PiMinus,
particles::Code::Pi0, particles::Code::KMinus,
particles::Code::KPlus, particles::Code::K0Long,
......@@ -82,12 +82,12 @@ namespace corsika::process::sibyll {
particles::Code::DMinus, particles::Code::D0,
particles::Code::MuMinus, particles::Code::MuPlus,
particles::Code::D0Bar};
const bool fInternalDecays = true;
const corsika::units::si::HEPEnergyType fMinEnergyCoM =
const bool internalDecays_ = true;
const corsika::units::si::HEPEnergyType minEnergyCoM_ =
10. * 1e9 * corsika::units::si::electronvolt;
const corsika::units::si::HEPEnergyType fMaxEnergyCoM =
const corsika::units::si::HEPEnergyType maxEnergyCoM_ =
1.e6 * 1e9 * corsika::units::si::electronvolt;
const int fMaxTargetMassNumber = 18;
const int maxTargetMassNumber_ = 18;
};
} // namespace corsika::process::sibyll
......
Processes/Sibyll/NuclearInteraction.cc
+ 26
− 26
View file @ 91cb1961
......@@ -39,19 +39,19 @@ namespace corsika::process::sibyll {
template <>
NuclearInteraction<SetupEnvironment>::NuclearInteraction(
process::sibyll::Interaction& hadint, SetupEnvironment const& env)
: fEnvironment(env)
, fHadronicInteraction(hadint) {}
: environment_(env)
, hadronicInteraction_(hadint) {}
template <>
NuclearInteraction<SetupEnvironment>::~NuclearInteraction() {
cout << "Nuclib::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount << endl;
cout << "Nuclib::NuclearInteraction n=" << count_ << " Nnuc=" << nucCount_ << endl;
}
template <>
void NuclearInteraction<SetupEnvironment>::PrintCrossSectionTable(
corsika::particles::Code pCode) {
using namespace corsika::particles;
const int k = fTargetComponentsIndex.at(pCode);
const int k = targetComponentsIndex_.at(pCode);
Code pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
Code::Neon, Code::Argon, Code::Iron};
cout << "en/A ";
......@@ -75,7 +75,7 @@ namespace corsika::process::sibyll {
using namespace corsika::particles;
using namespace units::si;
auto& universe = *(fEnvironment.GetUniverse());
auto& universe = *(environment_.GetUniverse());
auto const allElementsInUniverse = std::invoke([&]() {
std::set<particles::Code> allElementsInUniverse;
......@@ -98,13 +98,13 @@ namespace corsika::process::sibyll {
++k;
cout << "NuclearInteraction: init target component: " << ptarg << endl;
const int ib = GetNucleusA(ptarg);
if (!fHadronicInteraction.IsValidTarget(ptarg)) {
if (!hadronicInteraction_.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));
targetComponentsIndex_.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
......@@ -113,14 +113,14 @@ namespace corsika::process::sibyll {
// get p-p cross sections
auto const protonId = Code::Proton;
auto const [siginel, sigela] =
fHadronicInteraction.GetCrossSection(protonId, protonId, Ecm);
hadronicInteraction_.GetCrossSection(protonId, protonId, Ecm);
const double dsig = siginel / 1_mb;
const double dsigela = sigela / 1_mb;
// loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
for (int j = 1; j < gMaxNucleusAProjectile; ++j) {
for (int j = 1; j < gMaxNucleusAProjectile_; ++j) {
const int jj = j + 1;
double sig_out, dsig_out, sigqe_out, dsigqe_out;
sigma_mc_(jj, ib, dsig, dsigela, gNSample, sig_out, dsig_out, sigqe_out,
sigma_mc_(jj, ib, dsig, dsigela, gNSample_, sig_out, dsig_out, sigqe_out,
dsigqe_out);
// write to table
cnucsignuc_.sigma[j][k][i] = sig_out;
......@@ -128,7 +128,7 @@ namespace corsika::process::sibyll {
}
}
}
cout << "NuclearInteraction: cross sections for " << fTargetComponentsIndex.size()
cout << "NuclearInteraction: cross sections for " << targetComponentsIndex_.size()
<< " components initialized!" << endl;
for (auto& ptarg : allElementsInUniverse) {
cout << "cross section table: " << ptarg << endl;
......@@ -140,11 +140,11 @@ namespace corsika::process::sibyll {
void NuclearInteraction<SetupEnvironment>::Init() {
// initialize hadronic interaction module
// TODO: safe to run multiple initializations?
if (!fHadronicInteraction.WasInitialized()) fHadronicInteraction.Init();
if (!hadronicInteraction_.WasInitialized()) hadronicInteraction_.Init();
// check compatibility of energy ranges, someone could try to use low-energy model..
if (!fHadronicInteraction.IsValidCoMEnergy(GetMinEnergyPerNucleonCoM()) ||
!fHadronicInteraction.IsValidCoMEnergy(GetMaxEnergyPerNucleonCoM()))
if (!hadronicInteraction_.IsValidCoMEnergy(GetMinEnergyPerNucleonCoM()) ||
!hadronicInteraction_.IsValidCoMEnergy(GetMaxEnergyPerNucleonCoM()))
throw std::runtime_error(
"NuclearInteraction: hadronic interaction model incompatible!");
......@@ -161,7 +161,7 @@ namespace corsika::process::sibyll {
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
const int ib = targetComponentsIndex_.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!");
......@@ -202,7 +202,7 @@ namespace corsika::process::sibyll {
"NUCLIB!");
}
if (fHadronicInteraction.IsValidTarget(TargetId)) {
if (hadronicInteraction_.IsValidTarget(TargetId)) {
auto const sigProd = ReadCrossSectionTable(iBeamA, TargetId, LabEnergyPerNuc);
cout << "cross section (mb): " << sigProd / 1_mb << endl;
return std::make_tuple(sigProd, 0_mb);
......@@ -269,8 +269,8 @@ namespace corsika::process::sibyll {
// energy limits
// TODO: values depend on hadronic interaction model !! this is sibyll specific
if (ElabNuc >= 8.5_GeV && ECoMNN >= gMinEnergyPerNucleonCoM &&
ECoMNN < gMaxEnergyPerNucleonCoM) {
if (ElabNuc >= 8.5_GeV && ECoMNN >= gMinEnergyPerNucleonCoM_ &&
ECoMNN < gMaxEnergyPerNucleonCoM_) {
// get target from environment
/*
......@@ -347,7 +347,7 @@ namespace corsika::process::sibyll {
(vP.GetNuclearA() - vP.GetNuclearZ()) * particles::Neutron::GetMass();
cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << endl;
fCount++;
count_++;
const CoordinateSystem& rootCS =
RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
......@@ -360,7 +360,7 @@ namespace corsika::process::sibyll {
cout << "Interaction: time: " << tOrig << endl;
// projectile nucleon number
const int kAProj = vP.GetNuclearA(); // GetNucleusA(ProjId);
const int kAProj = vP.GetNuclearA();
if (kAProj > GetMaxNucleusAProjectile())
throw std::runtime_error("Projectile nucleus too large for NUCLIB!");
......@@ -400,7 +400,7 @@ namespace corsika::process::sibyll {
HEPEnergyType EcmNN = PtotNN4.GetNorm();
cout << "NuclearInteraction: nuc-nuc cm energy: " << EcmNN / 1_GeV << endl;
if (!fHadronicInteraction.IsValidCoMEnergy(EcmNN)) {
if (!hadronicInteraction_.IsValidCoMEnergy(EcmNN)) {
cout << "NuclearInteraction: nuc-nuc. CoM energy too low for hadronic "
"interaction model!"
<< endl;
......@@ -446,13 +446,13 @@ namespace corsika::process::sibyll {
cout << "target component: " << targetId << endl;
cout << "beam id: " << beamId << endl;
const auto [sigProd, sigEla] =
fHadronicInteraction.GetCrossSection(beamId, targetId, EcmNN);
hadronicInteraction_.GetCrossSection(beamId, targetId, EcmNN);
cross_section_of_components[i] = sigProd;
[[maybe_unused]] auto sigElaCopy = sigEla; // ONLY TO AVOID COMPILER WARNINGS
}
const auto targetCode =
mediumComposition.SampleTarget(cross_section_of_components, fRNG);
mediumComposition.SampleTarget(cross_section_of_components, RNG_);
cout << "Interaction: target selected: " << targetCode << endl;
/*
FOR NOW: allow nuclei with A<18 or protons only.
......@@ -463,7 +463,7 @@ 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 (!fHadronicInteraction.IsValidTarget(targetCode))
if (!hadronicInteraction_.IsValidTarget(targetCode))
throw std::runtime_error("target outside range. ");
// end of target sampling
......@@ -473,7 +473,7 @@ namespace corsika::process::sibyll {
// (needed to determine number of nucleon-nucleon scatterings)
const auto protonId = particles::Proton::GetCode();
const auto [prodCrossSection, elaCrossSection] =
fHadronicInteraction.GetCrossSection(protonId, protonId, EcmNN);
hadronicInteraction_.GetCrossSection(protonId, protonId, EcmNN);
const double sigProd = prodCrossSection / 1_mb;
const double sigEla = elaCrossSection / 1_mb;
// sample number of interactions (only input variables, output in common cnucms)
......@@ -602,7 +602,7 @@ namespace corsika::process::sibyll {
PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig});
// create inelastic interaction
cout << "calling HadronicInteraction..." << endl;
fHadronicInteraction.DoInteraction(inelasticNucleon);
hadronicInteraction_.DoInteraction(inelasticNucleon);
}
cout << "NuclearInteraction: DoInteraction: done" << endl;
......
Processes/Sibyll/NuclearInteraction.h
+ 17
− 17
View file @ 91cb1961
......@@ -27,8 +27,8 @@ namespace corsika::process::sibyll {
class NuclearInteraction
: public corsika::process::InteractionProcess<NuclearInteraction<TEnvironment>> {
int fCount = 0;
int fNucCount = 0;
int count_ = 0;
int nucCount_ = 0;
public:
NuclearInteraction(corsika::process::sibyll::Interaction&, TEnvironment const&);
......@@ -39,14 +39,14 @@ namespace corsika::process::sibyll {
corsika::units::si::CrossSectionType ReadCrossSectionTable(
const int, corsika::particles::Code, corsika::units::si::HEPEnergyType);
corsika::units::si::HEPEnergyType GetMinEnergyPerNucleonCoM() {
return gMinEnergyPerNucleonCoM;
return gMinEnergyPerNucleonCoM_;
}
corsika::units::si::HEPEnergyType GetMaxEnergyPerNucleonCoM() {
return gMaxEnergyPerNucleonCoM;
return gMaxEnergyPerNucleonCoM_;
}
int constexpr GetMaxNucleusAProjectile() { return gMaxNucleusAProjectile; }
int constexpr GetMaxNFragments() { return gMaxNFragments; }
int constexpr GetNEnergyBins() { return gNEnBins; }
int constexpr GetMaxNucleusAProjectile() { return gMaxNucleusAProjectile_; }
int constexpr GetMaxNFragments() { return gMaxNFragments_; }
int constexpr GetNEnergyBins() { return gNEnBins_; }
template <typename Particle>
std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType>
......@@ -59,21 +59,21 @@ namespace corsika::process::sibyll {
corsika::process::EProcessReturn DoInteraction(Projectile&);
private:
TEnvironment const& fEnvironment;
corsika::process::sibyll::Interaction& fHadronicInteraction;
std::map<corsika::particles::Code, int> fTargetComponentsIndex;
corsika::random::RNG& fRNG =
TEnvironment const& environment_;
corsika::process::sibyll::Interaction& hadronicInteraction_;
std::map<corsika::particles::Code, int> targetComponentsIndex_;
corsika::random::RNG& RNG_ =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
static constexpr int gNSample =
static constexpr int gNSample_ =
500; // number of samples in MC estimation of cross section
static constexpr int gMaxNucleusAProjectile = 56;
static constexpr int gNEnBins = 6;
static constexpr int gMaxNFragments = 60;
static constexpr int gMaxNucleusAProjectile_ = 56;
static constexpr int gNEnBins_ = 6;
static constexpr int gMaxNFragments_ = 60;
// energy limits defined by table used for cross section in signuc.f
// 10**1 GeV to 10**6 GeV
static constexpr corsika::units::si::HEPEnergyType gMinEnergyPerNucleonCoM =
static constexpr corsika::units::si::HEPEnergyType gMinEnergyPerNucleonCoM_ =
10. * 1e9 * corsika::units::si::electronvolt;
static constexpr corsika::units::si::HEPEnergyType gMaxEnergyPerNucleonCoM =
static constexpr corsika::units::si::HEPEnergyType gMaxEnergyPerNucleonCoM_ =
1.e6 * 1e9 * corsika::units::si::electronvolt;
};
......
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