diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc index be02d308c7b42deea043eae2413cc89ad843802f..7246ce77e9881131e721d89e3beb462fd6e2bcf9 100644 --- a/Documentation/Examples/vertical_EAS.cc +++ b/Documentation/Examples/vertical_EAS.cc @@ -74,6 +74,8 @@ void registerRandomStreams(const int seed) { int main(int argc, char** argv) { + logging::SetLevel(logging::level::info); + C8LOG_INFO("vertical_EAS"); if (argc < 4) { diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index f5f9040b05e343d83f019b1853f21ca53749b03c..8d2bd1743275ff6b334022d464001334ed7a92c5 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -126,9 +126,7 @@ namespace corsika::cascade { while (!fStack.IsEmpty()) { while (!fStack.IsEmpty()) { - C8LOG_TRACE(fmt::format("Stack: {}", fStack.as_string())); - for (const auto& tmp_p : fStack) - std::cout << " test " << tmp_p.GetPID() << std::endl; + C8LOG_TRACE(fmt::format("Stack: {}", fStack.as_string())); count_++; auto pNext = fStack.GetNextParticle(); C8LOG_DEBUG(fmt::format( @@ -189,7 +187,8 @@ namespace corsika::cascade { C8LOG_DEBUG( "total_lambda={} g/cm2, " ", next_interact={} g/cm2", - double((1. / total_inv_lambda) / 1_g * 1_cm * 1_cm), double(next_interact / 1_g * 1_cm * 1_cm)); + double((1. / total_inv_lambda) / 1_g * 1_cm * 1_cm), + double(next_interact / 1_g * 1_cm * 1_cm)); auto const* currentLogicalNode = vParticle.GetNode(); @@ -205,7 +204,7 @@ namespace corsika::cascade { // determine the maximum geometric step length from continuous processes LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step); - C8LOG_DEBUG("distance_max={} m", distance_max/1_m); + C8LOG_DEBUG("distance_max={} m", distance_max / 1_m); // determine combined total inverse decay time InverseTimeType const total_inv_lifetime = @@ -217,7 +216,7 @@ namespace corsika::cascade { C8LOG_DEBUG( "total_lifetime={} s" ", next_decay={} s", - (1/total_inv_lifetime)/1_s, next_decay/1_s); + (1 / total_inv_lifetime) / 1_s, next_decay / 1_s); // convert next_decay from time to length [m] LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() / @@ -227,7 +226,7 @@ namespace corsika::cascade { auto const min_distance = std::min({distance_interact, distance_decay, distance_max, geomMaxLength}); - C8LOG_DEBUG("transport particle by : {} m", min_distance/1_m); + C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m); // here the particle is actually moved along the trajectory to new position: // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step); diff --git a/Framework/Logging/Logging.h b/Framework/Logging/Logging.h index 774bffab01117476f62c66792328d00a9a1a2ce5..814f35de4de5c025c3bab8eaa7fdf51e135e7010 100644 --- a/Framework/Logging/Logging.h +++ b/Framework/Logging/Logging.h @@ -40,8 +40,9 @@ #define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_CRITICAL #endif +#include <spdlog/fmt/ostr.h> // will output whenerver a streaming operator is found +#include <spdlog/sinks/stdout_color_sinks.h> #include <spdlog/spdlog.h> -#include "spdlog/sinks/stdout_color_sinks.h" namespace corsika::logging { diff --git a/Framework/Particles/CMakeLists.txt b/Framework/Particles/CMakeLists.txt index 3c18f1572c6fd714cd6eb2aea979cedd9d817149..d512a409c9103dc054c37296c9f6b33ba624e86c 100644 --- a/Framework/Particles/CMakeLists.txt +++ b/Framework/Particles/CMakeLists.txt @@ -56,6 +56,7 @@ target_link_libraries ( CORSIKAparticles PUBLIC CORSIKAunits + spdlog::spdlog ) set_target_properties ( diff --git a/Framework/StackInterface/CombinedStack.h b/Framework/StackInterface/CombinedStack.h index 733f56b9a43b08a1045bb51eea44e0ffc251eccf..d8cecfd57bdf2fe43d549dffd1e4432818a6f961 100644 --- a/Framework/StackInterface/CombinedStack.h +++ b/Framework/StackInterface/CombinedStack.h @@ -8,10 +8,10 @@ #pragma once +#include <corsika/logging/Logging.h> #include <corsika/particles/ParticleProperties.h> #include <corsika/stack/Stack.h> #include <corsika/units/PhysicalUnits.h> -#include <corsika/logging/Logging.h> namespace corsika::stack { @@ -90,9 +90,8 @@ namespace corsika::stack { ///@} std::string as_string() const { - return fmt::format("[[{}][{}]]",PI_A::as_string(), PI_B::as_string()); + return fmt::format("[[{}][{}]]", PI_A::as_string(), PI_B::as_string()); } - }; /** @@ -115,7 +114,7 @@ namespace corsika::stack { unsigned int GetSize() const { return Stack1Impl::GetSize(); } unsigned int GetCapacity() const { return Stack1Impl::GetCapacity(); } - + /** * Function to copy particle at location i1 in stack to i2 */ diff --git a/Framework/StackInterface/SecondaryView.h b/Framework/StackInterface/SecondaryView.h index d69e6aabbddd241ab39136cdfa50032d077be10a..4d615a12005a11ff1a09364529c2392c8798fc0e 100644 --- a/Framework/StackInterface/SecondaryView.h +++ b/Framework/StackInterface/SecondaryView.h @@ -111,9 +111,6 @@ namespace corsika::stack { friend class ParticleBase<StackIterator>; - // template <template <typename> typename M1> - // friend class MSecondaryProducer; - private: /** * This is not accessible, since we don't want to allow creating a @@ -121,6 +118,7 @@ namespace corsika::stack { */ template <typename... Args> SecondaryView(Args... args) = delete; + SecondaryView() = delete; private: InnerStackTypeValue& inner_stack_; @@ -132,12 +130,26 @@ namespace corsika::stack { SecondaryView can only be constructed passing it a valid StackIterator to another Stack object **/ - SecondaryView(StackIteratorValue& vI) - : Stack<StackDataType&, ParticleInterface>(vI.GetStackData()) - , inner_stack_(vI.GetStack()) - , projectile_index_(vI.GetIndex()) { - C8LOG_TRACE("SecondaryView::SecondaryView"); - MSecondaryProducer<StackDataType, ParticleInterface>::new_view(vI); + SecondaryView(StackIteratorValue& particle) + : Stack<StackDataType&, ParticleInterface>(particle.GetStackData()) + , inner_stack_(particle.GetStack()) + , projectile_index_(particle.GetIndex()) { + C8LOG_TRACE("SecondaryView::SecondaryView(particle)"); + MSecondaryProducer<StackDataType, ParticleInterface>::new_view(particle); + } + /** + * Also allow to create a new View from a Projectile (StackIterator on View) + * + * Note, the view generated this way will be equivalent to the orignal view in + * terms of reference to the underlying data stack. It is not a "view to a view". + */ + SecondaryView(ViewType& view, StackIterator& projectile) + : Stack<StackDataType&, ParticleInterface>(view.GetStackData()) + , inner_stack_(view.inner_stack_) + , projectile_index_(view.GetIndexFromIterator(projectile.GetIndex())) { + C8LOG_TRACE("SecondaryView::SecondaryView(projectile)"); + StackIteratorValue particle(inner_stack_, projectile_index_); + MSecondaryProducer<StackDataType, ParticleInterface>::new_view(particle); } /** @@ -145,11 +157,19 @@ namespace corsika::stack { * SecondaryView is derived from. This projectile should not be * used to modify the Stack! */ - ConstStackIteratorValue parent() const { return ConstStackIteratorValue(inner_stack_, projectile_index_); } + /** + * This returns the projectile/parent in the original Stack, where this + * SecondaryView is derived from. This projectile should not be + * used to modify the Stack! + */ + StackIteratorValue asNewParent() const { + return StackIteratorValue(inner_stack_, projectile_index_); + } + /** * This return a projectile of this SecondaryView, which can be * used to modify the SecondaryView @@ -369,14 +389,15 @@ namespace corsika::stack { * performed. */ unsigned int GetIndexFromIterator(const unsigned int vI) const { - // this is too much: C8LOG_TRACE("SecondaryView::GetIndexFromIterator({})={}", vI, (vI?indices_[vI-1]:projectile_index_)); + // this is too much: C8LOG_TRACE("SecondaryView::GetIndexFromIterator({})={}", vI, + // (vI?indices_[vI-1]:projectile_index_)); if (vI == 0) return projectile_index_; return indices_[vI - 1]; } }; /** - * + * Class to handle the generation of new secondaries. */ template <class T1, template <class> class T2> class DefaultSecondaryProducer { diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h index 2158c0fd6ebdea20adcc02b09aed980c55e4c70e..6903aabd675bd166fc7871552397c315cfd50b11 100644 --- a/Framework/StackInterface/Stack.h +++ b/Framework/StackInterface/Stack.h @@ -8,14 +8,14 @@ #pragma once +#include <corsika/logging/Logging.h> #include <corsika/stack/StackIteratorInterface.h> #include <corsika/utl/MetaProgramming.h> -#include <corsika/logging/Logging.h> #include <stdexcept> +#include <string> #include <type_traits> #include <vector> -#include <string> /** All classes around management of particles on a stack. @@ -133,7 +133,7 @@ namespace corsika::stack { // template<typename>typename M2> template <class T2, template <class> class T3> class MSecondaryProducer> friend class SecondaryView; //<TStackData,MParticleInterface,M>; // access for - //SecondaryView + // SecondaryView friend class ParticleBase<StackIterator>; @@ -262,17 +262,19 @@ namespace corsika::stack { } std::string as_string() const { - std::string str(fmt::format("size {}, entries {}, deleted {} \n", getSize(), getEntries(), getDeleted())); + std::string str(fmt::format("size {}, entries {}, deleted {} \n", getSize(), + getEntries(), getDeleted())); // we make our own begin/end since we want ALL entries std::string new_line = " "; - for (unsigned int iPart = 0; iPart!=getSize(); ++iPart) { - ConstStackIterator itPart(*this, iPart); - str += fmt::format("{}{}{}", new_line, itPart.as_string(), (deleted_[itPart.GetIndex()]?" [deleted]":"")); - new_line = "\n "; + for (unsigned int iPart = 0; iPart != getSize(); ++iPart) { + ConstStackIterator itPart(*this, iPart); + str += fmt::format("{}{}{}", new_line, itPart.as_string(), + (deleted_[itPart.GetIndex()] ? " [deleted]" : "")); + new_line = "\n "; } return str; } - + /** * delete this particle */ diff --git a/Framework/StackInterface/StackIteratorInterface.h b/Framework/StackInterface/StackIteratorInterface.h index b0c86c4624a2baca42511cc26d2a2ca6b2176e6d..b790c4955f2fd62ed458a154fa7f9e594c97bac7 100644 --- a/Framework/StackInterface/StackIteratorInterface.h +++ b/Framework/StackInterface/StackIteratorInterface.h @@ -91,7 +91,7 @@ namespace corsika::stack { template <typename T, template <typename> typename T3> typename M2> // template <typename> typename M2> friend class SecondaryView; //<TStackData,TParticleInterface,M>; // access for - //SecondaryView + // SecondaryView template <typename T, template <typename> typename ParticleInterface> // friend class corsika::history::HistorySecondaryView; @@ -263,13 +263,13 @@ namespace corsika::stack { // TParticleInterface>; // // access for // SecondaryView - template <typename T1, //=TStackData, + template <typename T1, //=TStackData, template <typename> typename M1, //=TParticleInterface, // template <typename> typename M2> template <class T2, template <class> class T3> class MSecondaryProducer> friend class SecondaryView; //<TStackData,TParticleInterface,M>; // access for - //SecondaryView + // SecondaryView friend class StackIteratorInterface<TStackData, TParticleInterface, StackType>; diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt index 915aef2587a99c989ab635ec49712c9c46febdd1..856135c789ff02fa40fd77b3e18d93f68a8e3dfa 100644 --- a/Processes/Sibyll/CMakeLists.txt +++ b/Processes/Sibyll/CMakeLists.txt @@ -79,6 +79,7 @@ target_link_libraries ( CORSIKAenvironment CORSIKAsetup CORSIKArandom + CORSIKAlogging ) target_include_directories ( diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc index 1a814e0dd00f0684e228450541496e6c4dcdf85a..09580884c795ad462c488ab6f3c60734ebeb773b 100644 --- a/Processes/Sibyll/NuclearInteraction.cc +++ b/Processes/Sibyll/NuclearInteraction.cc @@ -19,10 +19,11 @@ #include <corsika/setup/SetupStack.h> #include <corsika/setup/SetupTrajectory.h> +#include <corsika/logging/Logging.h> + #include <set> +#include <sstream> -using std::cout; -using std::endl; using std::tuple; using std::vector; @@ -36,7 +37,8 @@ namespace corsika::process::sibyll { template <> NuclearInteraction<SetupEnvironment>::~NuclearInteraction() { - cout << "Nuclib::NuclearInteraction n=" << count_ << " Nnuc=" << nucCount_ << endl; + C8LOG_DEBUG( + fmt::format("Nuclib::NuclearInteraction n={} Nnuc={}", count_, nucCount_)); } template <> @@ -46,20 +48,22 @@ namespace corsika::process::sibyll { const int k = targetComponentsIndex_.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; + std::ostringstream table; + table << "Nuclear CrossSectionTable pCode=" << pCode << " :\n en/A "; + for (auto& j : pNuclei) table << std::setw(9) << j; + table << "\n"; // loop over energy bins - for (int i = 0; i < GetNEnergyBins(); ++i) { - cout << " " << i << " "; + for (unsigned int i = 0; i < GetNEnergyBins(); ++i) { + table << " " << i << " "; for (auto& n : pNuclei) { auto const j = GetNucleusA(n); - cout << " " << std::setprecision(5) << std::setw(8) - << cnucsignuc_.sigma[j - 1][k][i]; + table << " " << std::setprecision(5) << std::setw(8) + << cnucsignuc_.sigma[j - 1][k][i]; } - cout << endl; + table << "\n"; } + C8LOG_DEBUG(table.str()); } template <> @@ -82,23 +86,24 @@ namespace corsika::process::sibyll { return allElementsInUniverse; }); - cout << "NuclearInteraction: initializing nuclear cross sections..." << endl; + C8LOG_DEBUG("NuclearInteraction: initializing nuclear cross sections..."); // loop over target components, at most 4!! int k = -1; for (auto& ptarg : allElementsInUniverse) { ++k; - cout << "NuclearInteraction: init target component: " << ptarg << endl; + C8LOG_DEBUG(fmt::format("NuclearInteraction: init target component: {}", ptarg)); const int ib = GetNucleusA(ptarg); if (!hadronicInteraction_.IsValidTarget(ptarg)) { - cout << "NuclearInteraction::InitializeNuclearCrossSections: target nucleus? id=" - << ptarg << endl; + C8LOG_DEBUG(fmt::format( + "NuclearInteraction::InitializeNuclearCrossSections: target nucleus? id={}", + ptarg)); throw std::runtime_error( " target can not be handled by hadronic interaction model! "); } targetComponentsIndex_.insert(std::pair<Code, int>(ptarg, k)); // loop over energies, fNEnBins log. energy bins - for (int i = 0; i < GetNEnergyBins(); ++i) { + for (unsigned 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; @@ -109,7 +114,7 @@ namespace corsika::process::sibyll { 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 (unsigned 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, @@ -120,12 +125,11 @@ namespace corsika::process::sibyll { } } } - cout << "NuclearInteraction: cross sections for " << targetComponentsIndex_.size() - << " components initialized!" << endl; - for (auto& ptarg : allElementsInUniverse) { - cout << "cross section table: " << ptarg << endl; - PrintCrossSectionTable(ptarg); - } + C8LOG_DEBUG( + fmt::format("NuclearInteraction: cross sections for {} " + " components initialized!", + targetComponentsIndex_.size())); + for (auto& ptarg : allElementsInUniverse) { PrintCrossSectionTable(ptarg); } } template <> @@ -139,9 +143,9 @@ namespace corsika::process::sibyll { throw std::runtime_error("NuclearInteraction: energy outside tabulated range!"); const double e0 = elabnuc / 1_GeV; double sig; - cout << "ReadCrossSectionTable: " << ia << " " << ib << " " << e0 << endl; + C8LOG_DEBUG(fmt::format("ReadCrossSectionTable: {} {} {}", ia, ib, e0)); signuc2_(ia, ib, e0, sig); - cout << "ReadCrossSectionTable: sig=" << sig << endl; + C8LOG_DEBUG(fmt::format("ReadCrossSectionTable: sig={}", sig)); return sig * 1_mb; } @@ -156,19 +160,22 @@ namespace corsika::process::sibyll { throw std::runtime_error( "NuclearInteraction: GetCrossSection: particle not a nucleus!"); - auto const iBeamA = vP.GetNuclearA(); + const unsigned int iBeamA = vP.GetNuclearA(); HEPEnergyType LabEnergyPerNuc = vP.GetEnergy() / iBeamA; - cout << "NuclearInteraction: GetCrossSection: called with: beamNuclA= " << iBeamA - << " TargetId= " << TargetId << " LabEnergyPerNuc= " << LabEnergyPerNuc / 1_GeV - << endl; + C8LOG_DEBUG( + fmt::format("NuclearInteraction: GetCrossSection: called with: beamNuclA={} " + " TargetId={} LabEnergyPerNuc={}GeV ", + iBeamA, TargetId, LabEnergyPerNuc / 1_GeV)); // use nuclib to calc. nuclear cross sections // TODO: for now assumes air with hard coded composition // extend to arbitrary mixtures, requires smarter initialization // get nuclib projectile code: nucleon number if (iBeamA > GetMaxNucleusAProjectile() || iBeamA < 2) { - cout << "NuclearInteraction: beam nucleus outside allowed range for NUCLIB!" << endl - << "A=" << iBeamA << endl; + C8LOG_DEBUG( + "NuclearInteraction: beam nucleus outside allowed range for NUCLIB!" + "A=" + + std::to_string(iBeamA)); throw std::runtime_error( "NuclearInteraction: GetCrossSection: beam nucleus outside allowed range for " "NUCLIB!"); @@ -176,7 +183,7 @@ namespace corsika::process::sibyll { if (hadronicInteraction_.IsValidTarget(TargetId)) { auto const sigProd = ReadCrossSectionTable(iBeamA, TargetId, LabEnergyPerNuc); - cout << "cross section (mb): " << sigProd / 1_mb << endl; + C8LOG_DEBUG("cross section (mb): " + std::to_string(sigProd / 1_mb)); return std::make_tuple(sigProd, 0_mb); } else { throw std::runtime_error("target outside range."); @@ -219,7 +226,7 @@ namespace corsika::process::sibyll { // total momentum and energy HEPEnergyType Elab = vP.GetEnergy() + constants::nucleonMass; - int const nuclA = vP.GetNuclearA(); + unsigned int const nuclA = vP.GetNuclearA(); auto const ElabNuc = vP.GetEnergy() / nuclA; corsika::stack::MomentumVector pTotLab(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV}); @@ -230,13 +237,16 @@ namespace corsika::process::sibyll { const HEPEnergyType ECoM = sqrt( (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy 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 - << " beam pid:" << corsikaBeamId << endl - << " beam A: " << nuclA << endl - << " input energy per nucleon: " << ElabNuc / 1_GeV << endl - << " input energy CoM per nucleon: " << ECoMNN / 1_GeV << endl; + C8LOG_DEBUG( + fmt::format("NuclearInteraction: LambdaInt: \n" + " input energy: {}GeV\n" + " input energy CoM: {}GeV\n" + " beam pid: {}\n" + " beam A: {}\n" + " input energy per nucleon: {}GeV\n" + " input energy CoM per nucleon: {}GeV ", + Elab / 1_GeV, ECoM / 1_GeV, particles::GetName(corsikaBeamId), nuclA, + ElabNuc / 1_GeV, ECoMNN / 1_GeV)); // throw std::runtime_error("stop here"); // energy limits @@ -262,27 +272,30 @@ namespace corsika::process::sibyll { // loop over components in medium for (auto const targetId : mediumComposition.GetComponents()) { i++; - cout << "NuclearInteraction: get interaction length for target: " << targetId - << endl; + C8LOG_DEBUG("NuclearInteraction: get interaction length for target: " + + particles::GetName(targetId)); auto const [productionCrossSection, elaCrossSection] = GetCrossSection(vP, targetId); [[maybe_unused]] auto& dummy_elaCrossSection = elaCrossSection; - cout << "NuclearInteraction: " - << "IntLength: nuclib return (mb): " << productionCrossSection / 1_mb - << endl; + C8LOG_DEBUG( + "NuclearInteraction: " + "IntLength: nuclib return (mb): " + + std::to_string(productionCrossSection / 1_mb)); weightedProdCrossSection += w[i] * productionCrossSection; } - cout << "NuclearInteraction: " - << "IntLength: weighted CrossSection (mb): " << weightedProdCrossSection / 1_mb - << endl; + C8LOG_DEBUG( + "NuclearInteraction: " + "IntLength: weighted CrossSection (mb): " + + std::to_string(weightedProdCrossSection / 1_mb)); // calculate interaction length in medium GrammageType const int_length = mediumComposition.GetAverageMassNumber() * units::constants::u / weightedProdCrossSection; - cout << "NuclearInteraction: " - << "interaction length (g/cm2): " << int_length * (1_cm * 1_cm / (0.001_kg)) - << endl; + C8LOG_DEBUG( + "NuclearInteraction: " + "interaction length (g/cm2): " + + std::to_string(int_length * (1_cm * 1_cm / (0.001_kg)))); return int_length; } else { @@ -308,7 +321,8 @@ namespace corsika::process::sibyll { const auto ProjId = projectile.GetPID(); // TODO: calculate projectile mass in nuclearStackExtension // const auto ProjMass = projectile.GetMass(); - cout << "NuclearInteraction: DoInteraction: called with:" << ProjId << endl; + C8LOG_DEBUG("NuclearInteraction: DoInteraction: called with:" + + particles::GetName(ProjId)); // check if target-style nucleus (enum) if (ProjId != particles::Code::Nucleus) @@ -319,7 +333,8 @@ namespace corsika::process::sibyll { auto const ProjMass = projectile.GetNuclearZ() * particles::Proton::GetMass() + (projectile.GetNuclearA() - projectile.GetNuclearZ()) * particles::Neutron::GetMass(); - cout << "NuclearInteraction: projectile mass: " << ProjMass / 1_GeV << endl; + C8LOG_DEBUG("NuclearInteraction: projectile mass: " + + std::to_string(ProjMass / 1_GeV)); count_++; @@ -330,11 +345,12 @@ namespace corsika::process::sibyll { Point pOrig = projectile.GetPosition(); TimeType tOrig = projectile.GetTime(); - cout << "Interaction: position of interaction: " << pOrig.GetCoordinates() << endl; - cout << "Interaction: time: " << tOrig << endl; + C8LOG_DEBUG( + fmt::format("Interaction: position of interaction: {}", pOrig.GetCoordinates())); + C8LOG_DEBUG("Interaction: time: " + std::to_string(tOrig / 1_s)); // projectile nucleon number - const int kAProj = projectile.GetNuclearA(); + const unsigned int kAProj = projectile.GetNuclearA(); if (kAProj > GetMaxNucleusAProjectile()) throw std::runtime_error("Projectile nucleus too large for NUCLIB!"); @@ -344,18 +360,21 @@ namespace corsika::process::sibyll { auto const pProjectileLab = projectile.GetMomentum(); const FourVector PprojLab(eProjectileLab, pProjectileLab); - cout << "NuclearInteraction: eProj lab: " << eProjectileLab / 1_GeV << endl - << "NuclearInteraction: pProj lab: " << pProjectileLab.GetComponents() / 1_GeV - << endl; + C8LOG_DEBUG( + fmt::format("NuclearInteraction: eProj lab: {} " + "pProj lab: {} ", + eProjectileLab / 1_GeV, pProjectileLab.GetComponents() / 1_GeV)); + ; // define projectile nucleon HEPEnergyType const eProjectileNucLab = projectile.GetEnergy() / kAProj; auto const pProjectileNucLab = projectile.GetMomentum() / kAProj; const FourVector PprojNucLab(eProjectileNucLab, pProjectileNucLab); - cout << "NuclearInteraction: eProjNucleon lab: " << eProjectileNucLab / 1_GeV << endl - << "NuclearInteraction: pProjNucleon lab: " - << pProjectileNucLab.GetComponents() / 1_GeV << endl; + C8LOG_DEBUG(fmt::format( + "NuclearInteraction: eProjNucleon lab (GeV): {} " + "pProjNucleon lab (GeV): {}", + eProjectileNucLab / 1_GeV, pProjectileNucLab.GetComponents() / 1_GeV)); // define target // always a nucleon @@ -365,19 +384,21 @@ namespace corsika::process::sibyll { corsika::stack::MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV); const FourVector PtargNucLab(eTargetNucLab, pTargetNucLab); - cout << "NuclearInteraction: etarget lab: " << eTargetNucLab / 1_GeV << endl - << "NuclearInteraction: ptarget lab: " << pTargetNucLab.GetComponents() / 1_GeV - << endl; + C8LOG_DEBUG( + fmt::format("NuclearInteraction: etarget lab(GeV): {}" + "NuclearInteraction: ptarget lab(GeV): {} ", + eTargetNucLab / 1_GeV, pTargetNucLab.GetComponents() / 1_GeV)); // center-of-mass energy in nucleon-nucleon frame auto const PtotNN4 = PtargNucLab + PprojNucLab; HEPEnergyType EcmNN = PtotNN4.GetNorm(); - cout << "NuclearInteraction: nuc-nuc cm energy: " << EcmNN / 1_GeV << endl; + C8LOG_DEBUG("NuclearInteraction: nuc-nuc cm energy: " + + std::to_string(EcmNN / 1_GeV)); if (!hadronicInteraction_.IsValidCoMEnergy(EcmNN)) { - cout << "NuclearInteraction: nuc-nuc. CoM energy too low for hadronic " - "interaction model!" - << endl; + C8LOG_DEBUG( + "NuclearInteraction: nuc-nuc. CoM energy too low for hadronic " + "interaction model!"); throw std::runtime_error("NuclearInteraction: DoInteraction: energy too low!"); } @@ -389,14 +410,16 @@ namespace corsika::process::sibyll { // boost target auto const PtargNucCoM = boost.toCoM(PtargNucLab); - cout << "Interaction: ebeam CoM: " << PprojNucCoM.GetTimeLikeComponent() / 1_GeV - << endl - << "Interaction: pbeam CoM: " - << PprojNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl; - cout << "Interaction: etarget CoM: " << PtargNucCoM.GetTimeLikeComponent() / 1_GeV - << endl - << "Interaction: ptarget CoM: " - << PtargNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl; + C8LOG_DEBUG( + fmt::format("Interaction: ebeam CoM: {} " + ", pbeam CoM: {}", + PprojNucCoM.GetTimeLikeComponent() / 1_GeV, + PprojNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV)); + C8LOG_DEBUG( + fmt::format("Interaction: etarget CoM: {}" + ", ptarget CoM: {}", + PtargNucCoM.GetTimeLikeComponent() / 1_GeV, + PtargNucCoM.GetSpaceLikeComponents().GetComponents() / 1_GeV)); // sample target nucleon number // @@ -405,7 +428,7 @@ namespace corsika::process::sibyll { auto const* const currentNode = projectile.GetNode(); const auto& mediumComposition = currentNode->GetModelProperties().GetNuclearComposition(); - cout << "get nucleon-nucleus cross sections for target materials.." << endl; + C8LOG_DEBUG("get nucleon-nucleus cross sections for target materials.."); // get cross sections for target materials // using nucleon-target-nucleus cross section!!! /* @@ -417,8 +440,8 @@ namespace corsika::process::sibyll { for (size_t i = 0; i < compVec.size(); ++i) { auto const targetId = compVec[i]; - cout << "target component: " << targetId << endl; - cout << "beam id: " << beamId << endl; + C8LOG_DEBUG("target component: " + particles::GetName(targetId)); + C8LOG_DEBUG("beam id: " + particles::GetName(beamId)); const auto [sigProd, sigEla] = hadronicInteraction_.GetCrossSection(beamId, targetId, EcmNN); cross_section_of_components[i] = sigProd; @@ -427,7 +450,7 @@ namespace corsika::process::sibyll { const auto targetCode = mediumComposition.SampleTarget(cross_section_of_components, RNG_); - cout << "Interaction: target selected: " << targetCode << endl; + C8LOG_DEBUG("Interaction: target selected: " + particles::GetName(targetCode)); /* FOR NOW: allow nuclei with A<18 or protons only. when medium composition becomes more complex, approximations will have to be @@ -436,13 +459,13 @@ namespace corsika::process::sibyll { int kATarget = -1; if (IsNucleus(targetCode)) kATarget = GetNucleusA(targetCode); if (targetCode == particles::Proton::GetCode()) kATarget = 1; - cout << "NuclearInteraction: nuclib target code: " << kATarget << endl; + C8LOG_DEBUG("NuclearInteraction: nuclib target code: " + std::to_string(kATarget)); if (!hadronicInteraction_.IsValidTarget(targetCode)) throw std::runtime_error("target outside range. "); // end of target sampling // superposition - cout << "NuclearInteraction: sampling nuc. multiple interaction structure.. " << endl; + C8LOG_DEBUG("NuclearInteraction: sampling nuc. multiple interaction structure.. "); // get nucleon-nucleon cross section // (needed to determine number of nucleon-nucleon scatterings) const auto protonId = particles::Proton::GetCode(); @@ -454,24 +477,27 @@ namespace corsika::process::sibyll { // nuclear multiple scattering according to glauber (r.i.p.) int_nuc_(kATarget, kAProj, sigProd, sigEla); - cout << "number of nucleons in target : " << kATarget << endl - << "number of wounded nucleons in target : " << cnucms_.na << endl - << "number of nucleons in projectile : " << kAProj << endl - << "number of wounded nucleons in project. : " << cnucms_.nb << endl - << "number of inel. nuc.-nuc. interactions : " << cnucms_.ni << endl - << "number of elastic nucleons in target : " << cnucms_.nael << endl - << "number of elastic nucleons in project. : " << cnucms_.nbel << endl - << "impact parameter: " << cnucms_.b << endl; + C8LOG_DEBUG( + fmt::format("number of nucleons in target : {}\n" + "number of wounded nucleons in target : {}\n" + "number of nucleons in projectile : {}\n" + "number of wounded nucleons in project. : {}\n" + "number of inel. nuc.-nuc. interactions : {}\n" + "number of elastic nucleons in target : {}\n" + "number of elastic nucleons in project. : {}\n" + "impact parameter: {}", + kATarget, cnucms_.na, kAProj, cnucms_.nb, cnucms_.ni, cnucms_.nael, + cnucms_.nbel, cnucms_.b)); // calculate fragmentation - cout << "calculating nuclear fragments.." << endl; + C8LOG_DEBUG("calculating nuclear fragments.."); // number of interactions // include elastic const int nElasticNucleons = cnucms_.nbel; const int nInelNucleons = cnucms_.nb; const int nIntProj = nInelNucleons + nElasticNucleons; const double impactPar = cnucms_.b; // only needed to avoid passing common var. - int nFragments; + int nFragments = 0; // number of fragments is limited to 60 int AFragments[60]; // call fragmentation routine @@ -481,22 +507,22 @@ namespace corsika::process::sibyll { fragm_(kATarget, kAProj, nIntProj, impactPar, nFragments, AFragments); // this should not occur but well :) - if (nFragments > GetMaxNFragments()) + if (nFragments > (int)GetMaxNFragments()) throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!"); - cout << "number of fragments: " << nFragments << endl; + C8LOG_DEBUG("number of fragments: " + std::to_string(nFragments)); for (int j = 0; j < nFragments; ++j) - cout << "fragment: " << j << " A=" << AFragments[j] - << " px=" << fragments_.ppp[j][0] << " py=" << fragments_.ppp[j][1] - << " pz=" << fragments_.ppp[j][2] << endl; + C8LOG_DEBUG(fmt::format("fragment {}: A={} px={} py={} pz={}", j, AFragments[j], + fragments_.ppp[j][0], fragments_.ppp[j][1], + fragments_.ppp[j][2])); - cout << "adding nuclear fragments to particle stack.." << endl; + C8LOG_DEBUG("adding nuclear fragments to particle stack.."); // put nuclear fragments on corsika stack for (int j = 0; j < nFragments; ++j) { particles::Code specCode; - const auto nuclA = AFragments[j]; + const int nuclA = AFragments[j]; // get Z from stability line - const auto nuclZ = int(nuclA / 2.15 + 0.7); + const int nuclZ = int(nuclA / 2.15 + 0.7); // TODO: do we need to catch single nucleons?? if (nuclA == 1) @@ -510,20 +536,21 @@ namespace corsika::process::sibyll { particles::Proton::GetMass() * nuclZ + (nuclA - nuclZ) * particles::Neutron::GetMass(); // this neglects binding energy - cout << "NuclearInteraction: adding fragment: " << specCode << endl; - cout << "NuclearInteraction: A,Z: " << nuclA << "," << nuclZ << endl; - cout << "NuclearInteraction: mass: " << mass / 1_GeV << endl; + C8LOG_DEBUG("NuclearInteraction: adding fragment: " + particles::GetName(specCode)); + C8LOG_DEBUG("NuclearInteraction: A,Z: " + std::to_string(nuclA) + ", " + + std::to_string(nuclZ)); + C8LOG_DEBUG("NuclearInteraction: mass: " + std::to_string(mass / 1_GeV)); // CORSIKA 7 way // spectators inherit momentum from original projectile const double mass_ratio = mass / ProjMass; - cout << "NuclearInteraction: mass ratio " << mass_ratio << endl; + C8LOG_DEBUG("NuclearInteraction: mass ratio " + std::to_string(mass_ratio)); auto const Plab = PprojLab * mass_ratio; - cout << "NuclearInteraction: fragment momentum: " - << Plab.GetSpaceLikeComponents().GetComponents() / 1_GeV << endl; + C8LOG_DEBUG(fmt::format("NuclearInteraction: fragment momentum: {}", + Plab.GetSpaceLikeComponents().GetComponents() / 1_GeV)); if (nuclA == 1) // add nucleon @@ -545,7 +572,7 @@ namespace corsika::process::sibyll { // add elastic nucleons to corsika stack // TODO: the elastic interaction could be external like the inelastic interaction, // e.g. use existing ElasticModel - cout << "adding elastically scattered nucleons to particle stack.." << endl; + C8LOG_DEBUG("adding elastically scattered nucleons to particle stack.."); for (int j = 0; j < nElasticNucleons; ++j) { // TODO: sample proton or neutron auto const elaNucCode = particles::Code::Proton; @@ -564,23 +591,37 @@ namespace corsika::process::sibyll { } // add inelastic interactions - cout << "calculate inelastic nucleon-nucleon interactions.." << endl; + C8LOG_DEBUG("calculate inelastic nucleon-nucleon interactions.."); for (int j = 0; j < nInelNucleons; ++j) { // TODO: sample neutron or proton auto pCode = particles::Proton::GetCode(); // temporarily add to stack, will be removed after interaction in DoInteraction - cout << "inelastic interaction no. " << j << endl; - auto inelasticNucleon = projectile.AddSecondary( + C8LOG_DEBUG(fmt::format("inelastic interaction no. {}", j)); + setup::Stack nucleonStack; + // auto inelasticNucleon = projectile.AddSecondary( + auto inelasticNucleon = nucleonStack.AddParticle( tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{ pCode, PprojNucLab.GetTimeLikeComponent(), PprojNucLab.GetSpaceLikeComponents(), pOrig, tOrig}); - // create inelastic interaction - cout << "calling HadronicInteraction..." << endl; - //~ hadronicInteraction_.DoInteraction(inelasticNucleon); + inelasticNucleon.SetNode(projectile.GetNode()); + // create inelastic interaction for each nucleon + C8LOG_TRACE("calling HadronicInteraction..."); + // create new StackView for each of the nucleons + View nucleon_secondaries(inelasticNucleon); + // all inner hadronic event generator + hadronicInteraction_.DoInteraction(nucleon_secondaries); + // inelasticNucleon.Delete(); // this is just a temporary object + for (const auto& pSec : nucleon_secondaries) { + projectile.AddSecondary( + tuple<particles::Code, units::si::HEPEnergyType, + corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{ + pSec.GetPID(), pSec.GetEnergy(), pSec.GetMomentum(), pSec.GetPosition(), + pSec.GetTime()}); + } } - cout << "NuclearInteraction: DoInteraction: done" << endl; + C8LOG_DEBUG("NuclearInteraction: DoInteraction: done"); return process::EProcessReturn::eOk; } diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h index 3f5c4aeb6df7aedb7586d1a9247190e0e55e9698..4a561fad1298d78e8c45f03bd25711639366e1d6 100644 --- a/Processes/Sibyll/NuclearInteraction.h +++ b/Processes/Sibyll/NuclearInteraction.h @@ -41,9 +41,9 @@ namespace corsika::process::sibyll { corsika::units::si::HEPEnergyType GetMaxEnergyPerNucleonCoM() { return gMaxEnergyPerNucleonCoM_; } - int constexpr GetMaxNucleusAProjectile() { return gMaxNucleusAProjectile_; } - int constexpr GetMaxNFragments() { return gMaxNFragments_; } - int constexpr GetNEnergyBins() { return gNEnBins_; } + unsigned int constexpr GetMaxNucleusAProjectile() { return gMaxNucleusAProjectile_; } + unsigned int constexpr GetMaxNFragments() { return gMaxNFragments_; } + unsigned int constexpr GetNEnergyBins() { return gNEnBins_; } template <typename Particle> std::tuple<corsika::units::si::CrossSectionType, corsika::units::si::CrossSectionType> @@ -61,11 +61,11 @@ namespace corsika::process::sibyll { std::map<corsika::particles::Code, int> targetComponentsIndex_; corsika::random::RNG& RNG_ = corsika::random::RNGManager::GetInstance().GetRandomStream("sibyll"); - static constexpr int gNSample_ = + static constexpr unsigned 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 unsigned int gMaxNucleusAProjectile_ = 56; + static constexpr unsigned int gNEnBins_ = 6; + static constexpr unsigned 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_ = diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc index 45f3abaffa9c03a68a196647ef9b741d32f3628a..26814293241bfb877794687e7b6cf27790c947bb 100644 --- a/Processes/UrQMD/UrQMD.cc +++ b/Processes/UrQMD/UrQMD.cc @@ -8,6 +8,7 @@ #include <corsika/geometry/QuantityVector.h> #include <corsika/geometry/Vector.h> +#include <corsika/logging/Logging.h> #include <corsika/particles/ParticleProperties.h> #include <corsika/process/urqmd/UrQMD.h> #include <corsika/units/PhysicalUnits.h> @@ -39,6 +40,10 @@ CrossSectionType UrQMD::GetTabulatedCrossSection(particles::Code projectileCode, HEPEnergyType labEnergy) const { // translated to C++ from CORSIKA 7 subroutine cxtot_u + C8LOG_DEBUG("UrQMD::GetTabulatedCrossSection proj={}, targ={}, E={}GeV", + particles::GetName(projectileCode), particles::GetName(targetCode), + labEnergy / 1_GeV); + auto const kinEnergy = labEnergy - particles::GetMass(projectileCode); assert(kinEnergy >= HEPEnergyType::zero()); diff --git a/Setup/SetupStack.h b/Setup/SetupStack.h index b03e70abb5a954dd9f0c31a70586a71f6c508190..e72b191b38be385d2e3607806bdf9ba191ced9a8 100644 --- a/Setup/SetupStack.h +++ b/Setup/SetupStack.h @@ -86,7 +86,7 @@ namespace corsika::setup { typename corsika::setup::Stack::StackImpl, // CHECK with CLANG: corsika::setup::Stack::MPIType>; // corsika::setup::detail::StackWithGeometryInterface>; - corsika::setup::detail::StackWithHistoryInterface, StackViewProducer>; + corsika::setup::detail::StackWithHistoryInterface, StackViewProducer>; #elif defined(__GNUC__) || defined(__GNUG__) using TheStackView = corsika::stack::MakeView<corsika::setup::Stack, StackViewProducer>::type; diff --git a/Stack/DummyStack/DummyStack.h b/Stack/DummyStack/DummyStack.h index a6cfb02be6cba78e1544e0ab9c3cf1e542ee6d6d..26270b51cc0f8b9202ccdcba6dc268e6b29a3198 100644 --- a/Stack/DummyStack/DummyStack.h +++ b/Stack/DummyStack/DummyStack.h @@ -8,13 +8,13 @@ #pragma once -#include <corsika/particles/ParticleProperties.h> #include <corsika/logging/Logging.h> +#include <corsika/particles/ParticleProperties.h> #include <corsika/stack/Stack.h> #include <corsika/units/PhysicalUnits.h> -#include <tuple> #include <string> +#include <tuple> namespace corsika::stack { diff --git a/Stack/GeometryNodeStackExtension/GeometryNodeStackExtension.h b/Stack/GeometryNodeStackExtension/GeometryNodeStackExtension.h index 4b87872f2ffc5f4114a3140c60369b41422a5b22..d33997b8eecddca3d53be9fd4f1f7b1b95ad4835 100644 --- a/Stack/GeometryNodeStackExtension/GeometryNodeStackExtension.h +++ b/Stack/GeometryNodeStackExtension/GeometryNodeStackExtension.h @@ -10,7 +10,6 @@ #include <corsika/logging/Logging.h> #include <corsika/stack/Stack.h> -#include <corsika/logging/Logging.h> #include <tuple> #include <utility> diff --git a/Stack/History/HistoryObservationPlane.cc b/Stack/History/HistoryObservationPlane.cc index 16a3ae3d593b5e7de7d11eb079b806542b7f2df9..d9374eee7e00ee8177cd3027b86a9fac7bc6b9cd 100644 --- a/Stack/History/HistoryObservationPlane.cc +++ b/Stack/History/HistoryObservationPlane.cc @@ -8,7 +8,6 @@ #include <corsika/logging/Logging.h> #include <corsika/history/HistoryObservationPlane.hpp> -#include <corsika/logging/Logging.h> #include <boost/histogram/ostream.hpp> diff --git a/Stack/History/HistoryStackExtension.h b/Stack/History/HistoryStackExtension.h index 5a4881a6063092da95bdf95f72fd42c429161a8b..c9a0f289e0ab58bfb105678c44133208b710fa29 100644 --- a/Stack/History/HistoryStackExtension.h +++ b/Stack/History/HistoryStackExtension.h @@ -111,7 +111,8 @@ namespace corsika::history { } std::string as_string() const { - return fmt::format("i_parent={}, [evt: {}]", GetParentEventIndex(), (bool(GetEvent())?GetEvent()->as_string():"n/a")); + return fmt::format("i_parent={}, [evt: {}]", GetParentEventIndex(), + (bool(GetEvent()) ? GetEvent()->as_string() : "n/a")); } }; diff --git a/Stack/History/testHistoryView.cc b/Stack/History/testHistoryView.cc index 7e090e0b5921efaa6303ec25c0b8899796db0b62..d4c50e7a411b29b2efb93b3e125ff0565356b720 100644 --- a/Stack/History/testHistoryView.cc +++ b/Stack/History/testHistoryView.cc @@ -61,7 +61,8 @@ using TheTestStackView = corsika::stack::MakeView<TestStack, history::HistorySecondaryProducer>::type; #endif -using TestStackView = TheTestStackView; // history::HistorySecondaryView<TheTestStackView>; +using TestStackView = + TheTestStackView; // history::HistorySecondaryView<TheTestStackView>; TEST_CASE("HistoryStackExtension", "[stack]") {