IAP GITLAB

Commit ab32aa88 authored by Ralf Ulrich's avatar Ralf Ulrich

cascade and some exmples and modules

parent c2a13032
......@@ -114,11 +114,10 @@ endif (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
#
#+++++++++++++++++++++++++++++
# Setup external dependencies.
#
# the version of the cmake-conan project we use for build.
set (CMAKE_CONAN_VERSION "v0.13")
#
# download the cmake-conan integration
if (NOT EXISTS "${CMAKE_BINARY_DIR}/conan.cmake")
message (STATUS "Downloading conan.cmake from https://github.com/conan-io/cmake-conan")
......@@ -126,15 +125,18 @@ if (NOT EXISTS "${CMAKE_BINARY_DIR}/conan.cmake")
"${CMAKE_BINARY_DIR}/conan.cmake"
TLS_VERIFY ON)
endif ()
#
# add the cmake-conan integration
include (${CMAKE_BINARY_DIR}/conan.cmake)
#
# download and build all dependencies
conan_cmake_run (CONANFILE conanfile.txt
BASIC_SETUP CMAKE_TARGETS
BUILD missing)
#
# add cnpy temporarily. As long as SaveBoostHistogram does not save to parquet directly
#
add_subdirectory (externals/cnpy)
#
#+++++++++++++++++++++++++++++
# Coverage
......@@ -216,12 +218,16 @@ target_include_directories (
$<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>
$<INSTALL_INTERFACE:${CMAKE_INSTALL_PREFIX}>
)
# since CORSIKA8 is a header only library we must specify all link dependencies here:
target_link_libraries (
CORSIKA8
INTERFACE
PhysUnits
CONAN_PKG::eigen
CONAN_PKG::spdlog
CONAN_PKG::boost
cnpy # for SaveBoostHistogram
stdc++fs # experimental::filesystem
)
# those are needed, since some headers (namely GeneratedParticleProperties.inc) are produced by python script from ParticleData.xml
add_subdirectory (src)
......
This diff is collapsed.
......@@ -90,10 +90,12 @@ namespace corsika {
}
}
inline HEPMassType nucleus_mass(const int A, const int Z) {
auto const absA = std::abs(A);
auto const absZ = std::abs(Z);
return get_mass(Code::Proton) * absZ + (absA - absZ) * get_mass(Code::Neutron);
inline HEPMassType get_nucleus_mass(unsigned int const A, unsigned int const Z) {
return get_mass(Code::Proton) * Z + (A - Z) * get_mass(Code::Neutron);
}
std::initializer_list<Code> constexpr get_all_particles() {
return particle::detail::all_particles;
}
} // namespace corsika
......@@ -20,7 +20,7 @@ namespace corsika {
}
template <typename TTimeType, typename TSpaceVecType>
TSpaceVecType& FourVector<TTimeType, TSpaceVecType>::spaceLikeComponents() {
TSpaceVecType& FourVector<TTimeType, TSpaceVecType>::getSpaceLikeComponents() {
return spaceLike_;
}
......
......@@ -20,8 +20,8 @@ namespace corsika {
return normal_.dot(vP - center_) > LengthType::zero();
}
inline LengthType Plane::getDistanceTo(corsika::Point const& vP) const {
return (normal_ * (vP - center_).dot(normal_)).norm();
inline LengthType Plane::getDistanceTo(Point const& vP) const {
return (normal_ * (vP - center_).dot(normal_)).getNorm();
}
inline Point const& Plane::getCenter() const { return center_; }
......
......@@ -174,7 +174,7 @@ namespace corsika {
template <typename TDimension>
auto Vector<TDimension>::operator/(double const p) const {
return Vector<TDimension>(BaseVector<TDimension>::getCoordinateSystem(),
BaseVector<TDimension>::quantityVector() / p);
BaseVector<TDimension>::getQuantityVector() / p);
}
template <typename TDimension>
......
......@@ -25,9 +25,10 @@
namespace corsika {
template <typename TProcess1, typename TProcess2>
template <typename Particle, typename VTNType>
template <typename TParticle>
ProcessReturn ProcessSequence<TProcess1, TProcess2>::doBoundaryCrossing(
Particle& particle, VTNType const& from, VTNType const& to) {
TParticle& particle, typename TParticle::node_type const& from,
typename TParticle::node_type const& to) {
ProcessReturn ret = ProcessReturn::Ok;
if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process1_type>,
......@@ -103,19 +104,19 @@ namespace corsika {
template <typename TProcess1, typename TProcess2>
template <typename TParticle, typename TTrack>
LengthType ProcessSequence<TProcess1, TProcess2>::maxStepLength(TParticle& particle,
TTrack& vTrack) {
LengthType ProcessSequence<TProcess1, TProcess2>::getMaxStepLength(TParticle& particle,
TTrack& vTrack) {
LengthType max_length = // if no other process in the sequence implements it
std::numeric_limits<double>::infinity() * meter;
if constexpr (std::is_base_of_v<ContinuousProcess<process1_type>, process1_type> ||
t1ProcSeq) {
LengthType const len = A_.maxStepLength(particle, vTrack);
LengthType const len = A_.getMaxStepLength(particle, vTrack);
max_length = std::min(max_length, len);
}
if constexpr (std::is_base_of_v<ContinuousProcess<process2_type>, process2_type> ||
t2ProcSeq) {
LengthType const len = B_.maxStepLength(particle, vTrack);
LengthType const len = B_.getMaxStepLength(particle, vTrack);
max_length = std::min(max_length, len);
}
return max_length;
......
......@@ -54,8 +54,8 @@ namespace corsika {
///@}
template <template <typename> class TParticleInterfaceA,
template <typename> class TParticleInterfaceB, typename TStackIterator>
inline std::string CombinedParticleInterface<TParticleInterfaceA, TParticleInterfaceB, TStackIterator>::as_string() const {
return fmt::format("[[{}][{}]]", pi_a_type::as_string(), pi_b_type::as_string());
inline std::string CombinedParticleInterface<TParticleInterfaceA, TParticleInterfaceB, TStackIterator>::asString() const {
return fmt::format("[[{}][{}]]", pi_a_type::asString(), pi_b_type::asString());
}
......
......@@ -187,14 +187,14 @@ namespace corsika {
template <typename TStackDataType, template <typename> typename TParticleInterface,
template <typename T1, template <class> class T2> class MSecondaryProducer>
std::string SecondaryView<TStackDataType, TParticleInterface,
MSecondaryProducer>::as_string() const {
MSecondaryProducer>::asString() const {
std::string str(fmt::format("size {}\n", getSize()));
// we make our own begin/end since we want ALL entries
std::string new_line = " ";
for (unsigned int iPart = 0; iPart != getSize(); ++iPart) {
const_stack_view_iterator itPart(*this, iPart);
str += fmt::format(
"{}{}{}", new_line, itPart.as_string(),
"{}{}{}", new_line, itPart.asString(),
(inner_stack_.deleted_[getIndexFromIterator(itPart.getIndex())] ? " [deleted]"
: ""));
new_line = "\n ";
......
This diff is collapsed.
......@@ -19,14 +19,14 @@
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/logging/Logging.hpp>
// using namespace corsika::units::si;
namespace corsika {
COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pprojectile,
HEPMassType const massTarget)
inline COMBoost::COMBoost(FourVector<HEPEnergyType, MomentumVector> const& Pprojectile,
HEPMassType const massTarget)
: originalCS_{Pprojectile.getSpaceLikeComponents().getCoordinateSystem()}
, rotatedCS_{make_rotationToZ(originalCS_, Pprojectile.getSpaceLikeComponents())} {
, rotatedCS_{
make_rotationToZ(Pprojectile.getSpaceLikeComponents().getCoordinateSystem(),
Pprojectile.getSpaceLikeComponents())} {
auto const pProjectile = Pprojectile.getSpaceLikeComponents();
auto const pProjNormSquared = pProjectile.getSquaredNorm();
auto const pProjNorm = sqrt(pProjNormSquared);
......@@ -42,13 +42,13 @@ namespace corsika {
setBoost(coshEta, sinhEta);
CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta, coshEta,
boost_.determinant() - 1);
CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
coshEta, boost_.determinant() - 1);
}
COMBoost::COMBoost(Vector<hepmomentum_d> const& momentum, HEPEnergyType mass)
inline COMBoost::COMBoost(MomentumVector const& momentum, HEPEnergyType mass)
: originalCS_{momentum.getCoordinateSystem()}
, rotatedCS_{make_rotationToZ(originalCS_, momentum)} {
, rotatedCS_{make_rotationToZ(momentum.getCoordinateSystem(), momentum)} {
auto const squaredNorm = momentum.getSquaredNorm();
auto const norm = sqrt(squaredNorm);
auto const sinhEta = -norm / mass;
......@@ -57,7 +57,7 @@ namespace corsika {
}
template <typename FourVector>
FourVector COMBoost::toCoM(FourVector const& p) const {
inline FourVector COMBoost::toCoM(FourVector const& p) const {
auto pComponents = p.getSpaceLikeComponents().getComponents(rotatedCS_);
Eigen::Vector3d eVecRotated = pComponents.getEigenVector();
Eigen::Vector2d lab;
......@@ -70,11 +70,11 @@ namespace corsika {
eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
return FourVector(E_CoM, Vector<hepmomentum_d>(rotatedCS_, eVecRotated));
return FourVector(E_CoM, MomentumVector(rotatedCS_, eVecRotated));
}
template <typename FourVector>
FourVector COMBoost::fromCoM(FourVector const& p) const {
inline FourVector COMBoost::fromCoM(FourVector const& p) const {
auto pCM = p.getSpaceLikeComponents().getComponents(rotatedCS_);
auto const Ecm = p.getTimeLikeComponent();
......@@ -105,7 +105,7 @@ namespace corsika {
return f;
}
void COMBoost::setBoost(double coshEta, double sinhEta) {
inline void COMBoost::setBoost(double coshEta, double sinhEta) {
boost_ << coshEta, sinhEta, sinhEta, coshEta;
inverseBoost_ << coshEta, -sinhEta, -sinhEta, coshEta;
}
......
......@@ -8,7 +8,7 @@
#pragma once
#include <cnpy/cnpy.hpp>
#include <cnpy.hpp>
#include <boost/histogram.hpp>
......@@ -23,7 +23,7 @@ namespace corsika {
template <class Axes, class Storage>
inline void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
std::string const& filename, SaveMode mode = SaveMode::append) {
std::string const& filename, SaveMode mode) {
unsigned const rank = h.rank();
// append vs. overwrite
......
......@@ -12,68 +12,112 @@
#include <corsika/framework/logging/Logging.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/SlidingPlanarExponential.hpp>
namespace corsika {
void LayeredSphericalAtmosphereBuilder::checkRadius(LengthType r) const {
template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
typename... TModelArgs>
void LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra,
TModelArgs...>::checkRadius(LengthType r) const {
if (r <= previousRadius_) {
throw std::runtime_error("radius must be greater than previous");
}
}
void LayeredSphericalAtmosphereBuilder::setNuclearComposition(
NuclearComposition composition) {
template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
typename... TModelArgs>
void LayeredSphericalAtmosphereBuilder<
TMediumInterface, TMediumModelExtra,
TModelArgs...>::setNuclearComposition(NuclearComposition const& composition) {
composition_ = std::make_unique<NuclearComposition>(composition);
}
void LayeredSphericalAtmosphereBuilder::addExponentialLayer(GrammageType b,
LengthType c,
LengthType upperBoundary) {
auto const radius = seaLevel_ + upperBoundary;
template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
typename... TModelArgs>
void LayeredSphericalAtmosphereBuilder<
TMediumInterface, TMediumModelExtra,
TModelArgs...>::addExponentialLayer(GrammageType b, LengthType c,
LengthType upperBoundary) {
auto const radius = earthRadius_ + upperBoundary;
checkRadius(radius);
previousRadius_ = radius;
auto node = std::make_unique<VolumeTreeNode<IMediumModel>>(
auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
std::make_unique<Sphere>(center_, radius));
auto const rho0 = b / c;
CORSIKA_LOG_INFO("rho0 = {}, c = {}", rho0, c);
node->setModelProperties<SlidingPlanarExponential<IMediumModel>>(
center_, rho0, -c, *composition_, seaLevel_);
if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
// helper lambda in which the last 5 arguments to make_shared<...> are bound
auto lastBound = [&](auto... argPack) {
return std::make_shared<
TMediumModelExtra<SlidingPlanarExponential<TMediumInterface>>>(
argPack..., center_, rho0, -c, *composition_, earthRadius_);
};
// now unpack the additional arguments
auto model = std::apply(lastBound, additionalModelArgs_);
node->setModelProperties(std::move(model));
} else {
node->template setModelProperties<SlidingPlanarExponential<TMediumInterface>>(
center_, rho0, -c, *composition_, earthRadius_);
}
layers_.push(std::move(node));
}
void LayeredSphericalAtmosphereBuilder::addLinearLayer(LengthType c,
LengthType upperBoundary) {
auto const radius = seaLevel_ + upperBoundary;
template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
typename... TModelArgs>
void LayeredSphericalAtmosphereBuilder<
TMediumInterface, TMediumModelExtra,
TModelArgs...>::addLinearLayer(LengthType c, LengthType upperBoundary) {
auto const radius = earthRadius_ + upperBoundary;
checkRadius(radius);
previousRadius_ = radius;
GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
std::make_unique<Sphere>(center_, radius));
units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
auto const rho0 = b / c;
CORSIKA_LOG_INFO("rho0 = {}", rho0);
if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
// helper lambda in which the last 2 arguments to make_shared<...> are bound
auto lastBound = [&](auto... argPack) {
return std::make_shared<TMediumModelExtra<HomogeneousMedium<TMediumInterface>>>(
argPack..., rho0, *composition_);
};
auto node = std::make_unique<VolumeTreeNode<IMediumModel>>(
std::make_unique<Sphere>(center_, radius));
node->setModelProperties<HomogeneousMedium<IMediumModel>>(rho0, *composition_);
// now unpack the additional arguments
auto model = std::apply(lastBound, additionalModelArgs_);
node->setModelProperties(std::move(model));
} else {
node->template setModelProperties<HomogeneousMedium<TMediumInterface>>(
rho0, *composition_);
}
layers_.push(std::move(node));
}
Environment<IMediumModel> LayeredSphericalAtmosphereBuilder::assemble() {
Environment<IMediumModel> env;
template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
typename... TModelArgs>
Environment<TMediumInterface> LayeredSphericalAtmosphereBuilder<
TMediumInterface, TMediumModelExtra, TModelArgs...>::assemble() {
Environment<TMediumInterface> env;
assemble(env);
return env;
}
void LayeredSphericalAtmosphereBuilder::assemble(Environment<IMediumModel>& env) {
template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
typename... TModelArgs>
void LayeredSphericalAtmosphereBuilder<
TMediumInterface, TMediumModelExtra,
TModelArgs...>::assemble(Environment<TMediumInterface>& env) {
auto& universe = env.getUniverse();
auto* outmost = universe.get();
......@@ -86,4 +130,14 @@ namespace corsika {
}
}
template <typename TMediumInterface, template <typename> typename MExtraEnvirnoment>
struct make_layered_spherical_atmosphere_builder {
template <typename... TArgs>
static auto create(Point const& center, LengthType earthRadius, TArgs... args) {
return LayeredSphericalAtmosphereBuilder<TMediumInterface, MExtraEnvirnoment,
TArgs...>{std::forward<TArgs>(args)...,
center, earthRadius};
}
};
} // namespace corsika
......@@ -46,7 +46,7 @@ namespace corsika {
}
template <typename TFunction>
inline double NuclearComposition::getWeightedSum(TFunction const& func) const {
inline auto NuclearComposition::getWeightedSum(TFunction const& func) const {
using ResultQuantity = decltype(func(*components_.cbegin()));
auto const prod = [&](auto const compID, auto const fraction) {
......
......@@ -14,20 +14,30 @@
namespace corsika::particle_cut {
ParticleCut::ParticleCut(const HEPEnergyType eCut, bool em, bool inv)
: energy_cut_(eCut)
, doCutEm_(em)
, doCutInv_(inv)
, energy_(0_GeV)
, em_energy_(0_GeV)
, em_count_(0)
, inv_energy_(0_GeV)
, inv_count_(0) {}
template <typename TParticle>
bool ParticleCut::ParticleIsBelowEnergyCut(TParticle const& vP) const {
auto const energyLab = vP.GetEnergy();
bool ParticleCut::isBelowEnergyCut(TParticle const& vP) const {
auto const energyLab = vP.getEnergy();
// nuclei
if (vP.GetPID() == corsika::Code::Nucleus) {
if (vP.getPID() == Code::Nucleus) {
// calculate energy per nucleon
auto const ElabNuc = energyLab / vP.GetNuclearA();
return (ElabNuc < fECut);
auto const ElabNuc = energyLab / vP.getNuclearA();
return (ElabNuc < energy_cut_);
} else {
return (energyLab < fECut);
return (energyLab < energy_cut_);
}
}
bool ParticleCut::ParticleIsEmParticle(Code vCode) const {
bool ParticleCut::isEmParticle(Code vCode) const {
// FOR NOW: switch
switch (vCode) {
case Code::Gamma:
......@@ -39,7 +49,7 @@ namespace corsika::particle_cut {
}
}
bool ParticleCut::ParticleIsInvisible(Code vCode) const {
bool ParticleCut::isInvisible(Code vCode) const {
switch (vCode) {
case Code::NuE:
case Code::NuEBar:
......@@ -52,57 +62,74 @@ namespace corsika::particle_cut {
}
}
template <typename TParticle>
bool ParticleCut::checkCutParticle(const TParticle& particle) {
const Code pid = particle.getPID();
HEPEnergyType energy = particle.getEnergy();
CORSIKA_LOG_DEBUG(fmt::format("ParticleCut: checking {}, E= {} GeV, EcutTot={} GeV",
pid, energy / 1_GeV,
(em_energy_ + inv_energy_ + energy_) / 1_GeV));
if (doCutEm_ && isEmParticle(pid)) {
CORSIKA_LOG_DEBUG("removing em. particle...");
em_energy_ += energy;
em_count_ += 1;
return true;
} else if (doCutInv_ && isInvisible(pid)) {
CORSIKA_LOG_DEBUG("removing inv. particle...");
inv_energy_ += energy;
inv_count_ += 1;
return true;
} else if (isBelowEnergyCut(particle)) {
CORSIKA_LOG_DEBUG("removing low en. particle...");
energy_ += energy;
return true;
} else if (particle.getTime() > 10_ms) {
CORSIKA_LOG_DEBUG("removing OLD particle...");
energy_ += energy;
return true;
}
return false; // this particle will not be removed/cut
}
void ParticleCut::doSecondaries(corsika::setup::StackView& vS) {
auto particle = vS.begin();
while (particle != vS.end()) {
if (checkCutParticle(particle)) { particle.erase(); }
++particle; // next entry in SecondaryView
}
}
auto p = vS.begin();
while (p != vS.end()) {
const Code pid = p.GetPID();
HEPEnergyType energy = p.GetEnergy();
std::cout << "ProcessCut: DoSecondaries: " << pid << " E= " << energy
<< ", EcutTot=" << (fEmEnergy + fInvEnergy + fEnergy) / 1_GeV << " GeV"
<< std::endl;
if (ParticleIsEmParticle(pid)) {
std::cout << "removing em. particle..." << std::endl;
fEmEnergy += energy;
fEmCount += 1;
p.Delete();
} else if (ParticleIsInvisible(pid)) {
std::cout << "removing inv. particle..." << std::endl;
fInvEnergy += energy;
fInvCount += 1;
p.Delete();
} else if (ParticleIsBelowEnergyCut(p)) {
std::cout << "removing low en. particle..." << std::endl;
fEnergy += energy;
p.Delete();
} else if (p.GetTime() > 10_ms) {
std::cout << "removing OLD particle..." << std::endl;
fEnergy += energy;
p.Delete();
} else {
++p; // next entry in SecondaryView
}
ProcessReturn ParticleCut::doContinuous(corsika::setup::Stack::particle_type& particle,
corsika::setup::Trajectory const&) {
CORSIKA_LOG_TRACE("ParticleCut::DoContinuous");
if (checkCutParticle(particle)) {
CORSIKA_LOG_TRACE("removing during continuous");
particle.erase();
// signal to upstream code that this particle was deleted
return ProcessReturn::ParticleAbsorbed;
}
return ProcessReturn::Ok;
}
void ParticleCut::Init() {
fEmEnergy = 0_GeV;
fEmCount = 0;
fInvEnergy = 0_GeV;
fInvCount = 0;
fEnergy = 0_GeV;
// defineEmParticles();
void ParticleCut::showResults() {
CORSIKA_LOG_INFO(
" ******************************\n"
" energy in em. component (GeV): {} \n "
" no. of em. particles injected: {} \n "
" energy in inv. component (GeV): {} \n "
" no. of inv. particles injected: {} \n "
" energy below particle cut (GeV): {} \n"
" ******************************",
em_energy_ / 1_GeV, em_count_, inv_energy_ / 1_GeV, inv_count_, energy_ / 1_GeV);
}
void ParticleCut::ShowResults() {
std::cout << " ******************************" << std::endl
<< " ParticleCut: " << std::endl
<< " energy in em. component (GeV): " << fEmEnergy / 1_GeV << std::endl
<< " no. of em. particles injected: " << fEmCount << std::endl
<< " energy in inv. component (GeV): " << fInvEnergy / 1_GeV << std::endl
<< " no. of inv. particles injected: " << fInvCount << std::endl
<< " energy below particle cut (GeV): " << fEnergy / 1_GeV << std::endl
<< " ******************************" << std::endl;
void ParticleCut::reset() {
em_energy_ = 0_GeV;
em_count_ = 0;
inv_energy_ = 0_GeV;
inv_count_ = 0;
energy_ = 0_GeV;
}
} // namespace corsika::particle_cut
......@@ -43,19 +43,17 @@ namespace corsika::stack_inspector {
HEPEnergyType Etot = 0_GeV;
for (const auto& iterP : vS) {
HEPEnergyType E = iterP.GetEnergy();
HEPEnergyType E = iterP.getEnergy();
Etot += E;
if (ReportStack_) {
corsika::CoordinateSystem& rootCS =
corsika::RootCoordinateSystem::getInstance()
.GetRootCoordinateSystem(); // for printout
auto pos = iterP.GetPosition().GetCoordinates(rootCS);
CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem(); // for printout
auto pos = iterP.getPosition().getCoordinates(rootCS);
std::cout << "StackInspector: i=" << std::setw(5) << std::fixed << (i++)
<< ", id=" << std::setw(30) << iterP.GetPID() << " E=" << std::setw(15)
<< ", id=" << std::setw(30) << iterP.getPID() << " E=" << std::setw(15)
<< std::scientific << (E / 1_GeV) << " GeV, "
<< " pos=" << pos << " node = " << iterP.GetNode();
if (iterP.GetPID() == Code::Nucleus)
std::cout << " nuc_ref=" << iterP.GetNucleusRef();
<< " pos=" << pos << " node = " << iterP.getNode();
if (iterP.getPID() == Code::Nucleus)
std::cout << " nuc_ref=" << iterP.getNucleusRef();
std::cout << std::endl;
}
}
......@@ -75,15 +73,9 @@ namespace corsika::stack_inspector {
<< " time=" << std::put_time(std::localtime(&now_time), "%T")
<< ", running=" << elapsed_seconds.count() << " seconds"
<< " (" << std::setw(3) << int(progress * 100) << "%)"
<< ", nStep=" << getStep() << ", stackSize=" << vS.GetSize()
<< ", nStep=" << getStep() << ", stackSize=" << vS.getSize()
<< ", Estack=" << Etot / 1_GeV << " GeV"
<< ", ETA=" << std::put_time(std::localtime(&eta_time), "%T") << std::endl;
}
template <typename TStack>
void StackInspector<TStack>::Init() {
<