diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt index 9fe9c0093cfd4f2151a54c7a3ebcf99ec4230ba8..6b304129abbac66f4a62c67f7f015c1d0c3f7c35 100644 --- a/Processes/Sibyll/CMakeLists.txt +++ b/Processes/Sibyll/CMakeLists.txt @@ -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) # ..................................................... diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc index def7ca1cc996c2629613d8bef2f76efe234ee6f2..553792e8f945ed6fb414c55f0d898650012da9cf 100644 --- a/Processes/Sibyll/Interaction.cc +++ b/Processes/Sibyll/Interaction.cc @@ -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 diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h index a1b385474703ae7ef877ee375bc4315c6b8625bb..62302d09190c1044c65ccf52289de0915645dc9c 100644 --- a/Processes/Sibyll/Interaction.h +++ b/Processes/Sibyll/Interaction.h @@ -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 diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc index 384330c636a45a3dbeed8fc11ba2193e9dd1fd31..a79fd80c2a5bc751bf247b4a01cc725986f9501f 100644 --- a/Processes/Sibyll/NuclearInteraction.cc +++ b/Processes/Sibyll/NuclearInteraction.cc @@ -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; diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h index 4d12647f21b00415f8992fa796513fe9ee1eed5f..199901cb8cd3b1061b6b2c2a30ea1f37e77d0e03 100644 --- a/Processes/Sibyll/NuclearInteraction.h +++ b/Processes/Sibyll/NuclearInteraction.h @@ -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; };