diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index b632dda7765f723ffd79fc7dd1e14807fe4f57b5..d654067f1e6fc4c683c4f8aedb0a32877df2b860 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -22,7 +22,7 @@ variables: # # multi-step pipeline for each commit # -# Mote: "Draft/WIP:" merge request, non-Draft/non-WIP merge requests and "Ready for Code Review" MR are all +# Note: "Draft/WIP:" merge request, non-Draft/non-WIP merge requests and "Ready for Code Review" MR are all # handled explicitly # stages: diff --git a/CMakeLists.txt b/CMakeLists.txt index 1cb48d5e73caba1ed6df27fa8a191facd5186287..1a3992a0ae584b2da249e1ff191e0e3f82a4e8fd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -138,19 +138,35 @@ if (CMAKE_BUILD_TYPE STREQUAL Coverage) # compile coverage under -O0 to avoid any optimization, function elimation etc. add_compile_options ("-O0") + # search for local lcov + find_program (c8_lcov_bin lcov) + if (NOT c8_lcov_bin) + set (c8_lcov_bin "${PROJECT_SOURCE_DIR}/externals/lcov/bin/lcov") + message ("use C8 version of lcov ${c8_lcov_bin}") + endif () + + # search for local genhtml + find_program (c8_genhtml_bin genhtml) + if (NOT c8_genhtml_bin) + set (c8_genhtml_bin "${PROJECT_SOURCE_DIR}/externals/lcov/bin/genhtml") + message ("use C8 version of genhtml ${c8_genhtml_bin}") + endif () + set (GCOV gcov CACHE STRING "gcov executable" FORCE) - set (LCOV_BIN_DIR "${PROJECT_SOURCE_DIR}/externals/lcov/bin") + + # collect coverage data add_custom_command ( OUTPUT raw-coverage.info COMMAND ${CMAKE_COMMAND} -E echo "Note: you need to run ctest at least once to generate the coverage data" - COMMAND ${LCOV_BIN_DIR}/lcov --gcov-tool=${GCOV} --directory . --capture --output-file raw-coverage.info + COMMAND ${c8_lcov_bin} --gcov-tool=${GCOV} --rc lcov_branch_coverage=1 + --directory . --capture --output-file raw-coverage.info ) # remove uninteresting entries add_custom_command ( OUTPUT coverage.info - COMMAND ${LCOV_BIN_DIR}/lcov -q --remove raw-coverage.info "*/usr/*" "/usr/*" --output-file coverage2.info - COMMAND ${LCOV_BIN_DIR}/lcov --remove coverage2.info + COMMAND ${c8_lcov_bin} -q --remove raw-coverage.info "*/usr/*" "/usr/*" --output-file coverage2.info + COMMAND ${c8_lcov_bin} --remove coverage2.info "*/externals/*" "*/tests/*" "*/sibyll2.3d.cpp" "*/.conan/*" "*/include/Pythia8/*" "*/install/*" "${CMAKE_SOURCE_DIR}/modules/*" "${CMAKE_BINARY_DIR}/modules/*" @@ -161,7 +177,7 @@ if (CMAKE_BUILD_TYPE STREQUAL Coverage) # generate html report add_custom_command ( OUTPUT coverage-report - COMMAND ${LCOV_BIN_DIR}/genhtml --demangle-cpp coverage.info -o coverage-report + COMMAND ${c8_genhtml_bin} --branch-coverage --show-details --title "CORSIKA 8 test coverage" --legend --demangle-cpp coverage.info -o coverage-report DEPENDS coverage.info ) add_custom_target (coverage DEPENDS coverage-report) diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl index 82b6887a3edcf18ff4f2044adc0c1d33efb9f9fe..78f32c2fd34438b0836b02627fea377a0102b2c6 100644 --- a/corsika/detail/framework/core/Cascade.inl +++ b/corsika/detail/framework/core/Cascade.inl @@ -33,8 +33,32 @@ namespace corsika { + template <typename TTracking, typename TProcessList, typename TOutput, typename TStack> + inline Cascade<TTracking, TProcessList, TOutput, TStack>::Cascade( + Environment<medium_interface_type> const& env, TTracking& tr, TProcessList& pl, + TOutput& out, TStack& stack) + : environment_(env) + , tracking_(tr) + , sequence_(pl) + , output_(out) + , stack_(stack) + , forceInteraction_(false) { + CORSIKA_LOG_INFO(c8_ascii_); + CORSIKA_LOG_INFO("This is CORSIKA {}.{}.{}.{}", CORSIKA_RELEASE_NUMBER, + CORSIKA_MAJOR_NUMBER, CORSIKA_MINOR_NUMBER, CORSIKA_PATCH_NUMBER); + CORSIKA_LOG_INFO("Tracking algorithm: {} (version {})", TTracking::getName(), + TTracking::getVersion()); + if constexpr (stack_view_type::has_event) { + CORSIKA_LOG_INFO("Stack - with full cascade HISTORY."); + } + } + template <typename TTracking, typename TProcessList, typename TOutput, typename TStack> inline void Cascade<TTracking, TProcessList, TOutput, TStack>::run() { + + // trigger the start of the outputs for this shower + output_.startOfShower(); + setNodes(); // put each particle on stack in correct environment volume while (!stack_.isEmpty()) { @@ -60,39 +84,14 @@ namespace corsika { // thus, the double loop sequence_.doCascadeEquations(stack_); } + + // indicate end of shower + output_.endOfShower(); } template <typename TTracking, typename TProcessList, typename TOutput, typename TStack> inline void Cascade<TTracking, TProcessList, TOutput, TStack>::forceInteraction() { - CORSIKA_LOG_TRACE("forced interaction!"); - setNodes(); - auto particle = stack_.getNextParticle(); - stack_view_type secondaries(particle); - - auto const* currentLogicalNode = particle.getNode(); - // assert that particle stays outside void Universe if it has no - // model properties set - assert((currentLogicalNode != &*environment_.getUniverse() || - environment_.getUniverse()->hasModelProperties()) && - "FATAL: The environment model has no valid properties set!"); - NuclearComposition const& composition = - currentLogicalNode->getModelProperties().getNuclearComposition(); - - // determine projectile - HEPEnergyType const Elab = particle.getEnergy(); - FourMomentum const projectileP4{Elab, particle.getMomentum()}; - // determine cross section in material - CrossSectionType const sigma = - composition.getWeightedSum([=](Code const targetId) -> CrossSectionType { - FourMomentum const targetP4( - get_mass(targetId), - MomentumVector(particle.getMomentum().getCoordinateSystem(), - {0_GeV, 0_GeV, 0_GeV})); - return sequence_.getCrossSection(particle, targetId, targetP4); - }); - interaction(secondaries, projectileP4, composition, sigma); - sequence_.doSecondaries(secondaries); - particle.erase(); // primary particle is done + forceInteraction_ = true; } template <typename TTracking, typename TProcessList, typename TOutput, typename TStack> @@ -126,6 +125,16 @@ namespace corsika { return sequence_.getCrossSection(particle, targetId, targetP4); }); + if (forceInteraction_) { + CORSIKA_LOG_TRACE("forced interaction!"); + forceInteraction_ = false; // just one (first) interaction + stack_view_type secondaries(particle); + interaction(secondaries, projectileP4, composition, total_cx); + sequence_.doSecondaries(secondaries); + particle.erase(); // primary particle is done + return; + } + // calculate interaction length in medium GrammageType const total_lambda = (composition.getAverageMassNumber() * constants::u) / total_cx; @@ -302,7 +311,7 @@ namespace corsika { sequence_.doSecondaries(secondaries); particle.erase(); - } // namespace corsika + } template <typename TTracking, typename TProcessList, typename TOutput, typename TStack> inline ProcessReturn Cascade<TTracking, TProcessList, TOutput, TStack>::decay( diff --git a/corsika/detail/framework/core/ParticleProperties.inl b/corsika/detail/framework/core/ParticleProperties.inl index 73c2e6f89cc1df311bc312a386b87945e5a1d6e5..fc9fbf636241f4c83d364cb9a6da80c458a3ecc8 100644 --- a/corsika/detail/framework/core/ParticleProperties.inl +++ b/corsika/detail/framework/core/ParticleProperties.inl @@ -14,7 +14,7 @@ namespace corsika { - inline HEPEnergyType constexpr calculate_kinetic_energy_threshold(Code const code) { + inline HEPEnergyType constexpr get_kinetic_energy_threshold(Code const code) { if (is_nucleus(code)) return particle::detail::threshold_nuclei; return particle::detail::thresholds[static_cast<CodeIntType>(code)]; } @@ -32,6 +32,8 @@ namespace corsika { return particle::detail::masses[static_cast<CodeIntType>(code)]; } + inline bool constexpr is_charged(Code const c) { return get_charge_number(c) != 0; } + inline bool constexpr is_nucleus(Code const code) { return code >= Code::Nucleus; } inline Code constexpr get_nucleus_code(size_t const A, @@ -60,12 +62,18 @@ namespace corsika { A * 10); // 10LZZZAAAI } + inline PDGCode constexpr get_PDG(unsigned int const A, unsigned int const Z) { + return PDGCode(1000000000 + Z * 10000 + A * 10); + } + inline int16_t constexpr get_charge_number(Code const code) { if (is_nucleus(code)) return get_nucleus_Z(code); return particle::detail::electric_charges[static_cast<CodeIntType>(code)]; } inline ElectricChargeType constexpr get_charge(Code const code) { + if (code == Code::Nucleus) + throw std::runtime_error("charge of particle::Nucleus undefined"); return get_charge_number(code) * constants::e; } diff --git a/corsika/detail/modules/LongitudinalProfile.inl b/corsika/detail/modules/LongitudinalProfile.inl index 5c9b54c60f7bd913ed78d36c1965ea4df56161ab..542a82c073022e89d2b289f6b2f2697323aac5fb 100644 --- a/corsika/detail/modules/LongitudinalProfile.inl +++ b/corsika/detail/modules/LongitudinalProfile.inl @@ -14,66 +14,31 @@ #include <cmath> #include <iomanip> #include <limits> +#include <utility> namespace corsika { - inline LongitudinalProfile::LongitudinalProfile(ShowerAxis const& shower_axis, - GrammageType dX) - : dX_(dX) - , shower_axis_{shower_axis} - , profiles_{static_cast<unsigned int>(shower_axis.getMaximumX() / dX_) + 1} {} + template <typename TOutput> + template <typename... TArgs> + inline LongitudinalProfile<TOutput>::LongitudinalProfile(TArgs&&... args) + : TOutput(std::forward<TArgs>(args)...) {} + template <typename TOutput> template <typename TParticle, typename TTrack> - inline ProcessReturn LongitudinalProfile::doContinuous(TParticle const& vP, - TTrack const& vTrack, - bool const) { - auto const pid = vP.getPID(); - - GrammageType const grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0)); - GrammageType const grammageEnd = shower_axis_.getProjectedX(vTrack.getPosition(1)); - - CORSIKA_LOG_DEBUG("longprof: pos1={} m, pos2={}, X1={} g/cm2, X2={} g/cm2", - vTrack.getPosition(0).getCoordinates() / 1_m, - vTrack.getPosition(1).getCoordinates() / 1_m, - grammageStart / 1_g * square(1_cm), - grammageEnd / 1_g * square(1_cm)); - - // Note: particle may go also "upward", thus, grammageEnd<grammageStart - int const binStart = std::ceil(grammageStart / dX_); - int const binEnd = std::floor(grammageEnd / dX_); - - for (int b = binStart; b <= binEnd; ++b) { - if (pid == Code::Photon) { - profiles_.at(b)[ProfileIndex::Photon]++; - } else if (pid == Code::Positron) { - profiles_.at(b)[ProfileIndex::Positron]++; - } else if (pid == Code::Electron) { - profiles_.at(b)[ProfileIndex::Electron]++; - } else if (pid == Code::MuPlus) { - profiles_.at(b)[ProfileIndex::MuPlus]++; - } else if (pid == Code::MuMinus) { - profiles_.at(b)[ProfileIndex::MuMinus]++; - } else if (is_hadron(pid)) { - profiles_.at(b)[ProfileIndex::Hadron]++; - } else if (is_neutrino(pid)) { - profiles_.at(b)[ProfileIndex::Invisible]++; - } - } + inline ProcessReturn LongitudinalProfile<TOutput>::doContinuous( + TParticle const& particle, TTrack const& track, bool const) { + auto const pid = particle.getPID(); + this->write(track, pid, 1.0); // weight hardcoded so far return ProcessReturn::Ok; } - inline void LongitudinalProfile::save(std::string const& filename, const int width, - const int precision) { - CORSIKA_LOG_DEBUG("Write longprof to {}", filename); - std::ofstream f{filename}; - f << "# X / g·cm¯², photon, e+, e-, mu+, mu-, all hadrons, neutrinos" << std::endl; - for (size_t b = 0; b < profiles_.size(); ++b) { - f << std::setprecision(5) << std::setw(11) << b * (dX_ / (1_g / 1_cm / 1_cm)); - for (auto const& N : profiles_.at(b)) { - f << std::setw(width) << std::setprecision(precision) << std::scientific << N; - } - f << std::endl; - } + template <typename TOutput> + inline YAML::Node LongitudinalProfile<TOutput>::getConfig() const { + YAML::Node node; + node["type"] = "LongitudinalProfile"; + + return node; } + } // namespace corsika diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl index 7100db181d85ba603be2bda29eba9758a0bff8aa..8f4a658b504d896f674943f4729cc2a30336329f 100644 --- a/corsika/detail/modules/ObservationPlane.inl +++ b/corsika/detail/modules/ObservationPlane.inl @@ -9,15 +9,16 @@ namespace corsika { template <typename TTracking, typename TOutput> + template <typename... TArgs> ObservationPlane<TTracking, TOutput>::ObservationPlane(Plane const& obsPlane, DirectionVector const& x_axis, - bool deleteOnHit) - : plane_(obsPlane) - , deleteOnHit_(deleteOnHit) - , energy_ground_(0_GeV) - , count_ground_(0) + bool const deleteOnHit, + TArgs&&... args) + : TOutput(std::forward<TArgs>(args)...) + , plane_(obsPlane) , xAxis_(x_axis.normalized()) - , yAxis_(obsPlane.getNormal().cross(xAxis_)) {} + , yAxis_(obsPlane.getNormal().cross(xAxis_)) + , deleteOnHit_(deleteOnHit) {} template <typename TTracking, typename TOutput> template <typename TParticle, typename TTrajectory> @@ -50,19 +51,18 @@ namespace corsika { Vector const displacement = pointOfIntersection - plane_.getCenter(); // add our particles to the output file stream + double const weight = 1.; // particle.getWeight() this->write(particle.getPID(), energy, displacement.dot(xAxis_), - displacement.dot(yAxis_), particle.getTime()); + displacement.dot(yAxis_), 0_m, weight); CORSIKA_LOG_TRACE("Particle detected absorbed={}", deleteOnHit_); if (deleteOnHit_) { - count_ground_++; - energy_ground_ += energy; return ProcessReturn::ParticleAbsorbed; } else { return ProcessReturn::Ok; } - } + } // namespace corsika template <typename TTracking, typename TOutput> template <typename TParticle, typename TTrajectory> @@ -90,17 +90,6 @@ namespace corsika { return trajectory.getLength(fractionOfIntersection); } - template <typename TTracking, typename TOutput> - inline void ObservationPlane<TTracking, TOutput>::showResults() const { - CORSIKA_LOG_INFO( - " ******************************\n" - " ObservationPlane: \n" - " energy an ground (GeV) : {}\n" - " no. of particles at ground : {}\n" - " ******************************", - energy_ground_ / 1_GeV, count_ground_); - } - template <typename TTracking, typename TOutput> inline YAML::Node ObservationPlane<TTracking, TOutput>::getConfig() const { using namespace units::si; @@ -110,7 +99,7 @@ namespace corsika { // basic info node["type"] = "ObservationPlane"; - node["units"] = "m"; // add default units for values + node["units"]["length"] = "m"; // add default units for values // the center of the plane auto const center{plane_.getCenter()}; @@ -144,10 +133,4 @@ namespace corsika { return node; } - template <typename TTracking, typename TOutput> - inline void ObservationPlane<TTracking, TOutput>::reset() { - energy_ground_ = 0_GeV; - count_ground_ = 0; - } - } // namespace corsika diff --git a/corsika/detail/modules/ParticleCut.inl b/corsika/detail/modules/ParticleCut.inl index 0358bbb5a158ccf76765fac902280c3832618e7a..ac72089c166c3a183424b9bb0e6931e565e0634a 100644 --- a/corsika/detail/modules/ParticleCut.inl +++ b/corsika/detail/modules/ParticleCut.inl @@ -12,19 +12,19 @@ namespace corsika { - inline ParticleCut::ParticleCut(HEPEnergyType const eEleCut, - HEPEnergyType const ePhoCut, - HEPEnergyType const eHadCut, HEPEnergyType const eMuCut, - bool const inv) - : doCutEm_(false) - , doCutInv_(inv) - , energy_cut_(0_GeV) - , energy_timecut_(0_GeV) - , energy_emcut_(0_GeV) - , energy_invcut_(0_GeV) - , em_count_(0) - , inv_count_(0) - , energy_count_() { + template <typename TOutput> + template <typename... TArgs> + inline ParticleCut<TOutput>::ParticleCut(HEPEnergyType const eEleCut, + HEPEnergyType const ePhoCut, + HEPEnergyType const eHadCut, + HEPEnergyType const eMuCut, bool const inv, + TArgs&&... outputArgs) + : TOutput(std::forward<TArgs>(outputArgs)...) + , cut_electrons_(eEleCut) + , cut_photons_(ePhoCut) + , cut_muons_(eMuCut) + , cut_hadrons_(eHadCut) + , doCutInv_(inv) { for (auto p : get_all_particles()) { if (is_hadron(p)) // nuclei are also hadrons set_kinetic_energy_threshold(p, eHadCut); @@ -43,137 +43,96 @@ namespace corsika { eEleCut / 1_GeV, ePhoCut / 1_GeV, eHadCut / 1_GeV, eMuCut / 1_GeV); } - inline ParticleCut::ParticleCut(HEPEnergyType const eHadCut, HEPEnergyType const eMuCut, - bool const inv) - : doCutEm_(true) - , doCutInv_(inv) - , energy_cut_(0_GeV) - , energy_timecut_(0_GeV) - , energy_emcut_(0_GeV) - , energy_invcut_(0_GeV) - , em_count_(0) - , inv_count_(0) - , energy_count_(0) { - - for (auto p : get_all_particles()) { - if (is_hadron(p)) - set_kinetic_energy_threshold(p, eHadCut); - else if (is_muon(p)) - set_kinetic_energy_threshold(p, eMuCut); - } - set_kinetic_energy_threshold(Code::Nucleus, eHadCut); - CORSIKA_LOG_DEBUG( - "setting thresholds: hadrons = {} GeV, " - "muons = {} GeV", - eHadCut / 1_GeV, eMuCut / 1_GeV); - } - - inline ParticleCut::ParticleCut(HEPEnergyType const eCut, bool const em, bool const inv) - : doCutEm_(em) - , doCutInv_(inv) - , energy_cut_(0_GeV) - , energy_timecut_(0_GeV) - , energy_emcut_(0_GeV) - , energy_invcut_(0_GeV) - , em_count_(0) - , inv_count_(0) - , energy_count_(0) { - for (auto p : get_all_particles()) set_kinetic_energy_threshold(p, eCut); + template <typename TOutput> + template <typename... TArgs> + inline ParticleCut<TOutput>::ParticleCut(HEPEnergyType const eCut, bool const inv, + TArgs&&... outputArgs) + : TOutput(std::forward<TArgs>(outputArgs)...) + , doCutInv_(inv) { + for (auto p : get_all_particles()) { set_kinetic_energy_threshold(p, eCut); } set_kinetic_energy_threshold(Code::Nucleus, eCut); - CORSIKA_LOG_DEBUG("setting kinetic energy threshold for all particles to {} GeV", - eCut / 1_GeV); + CORSIKA_LOG_DEBUG("setting kinetic energy threshold {} GeV", eCut / 1_GeV); } - inline ParticleCut::ParticleCut( - std::unordered_map<Code const, HEPEnergyType const> const& eCuts, bool const em, - bool const inv) - : doCutEm_(em) - , doCutInv_(inv) - , energy_cut_(0_GeV) - , energy_timecut_(0_GeV) - , energy_emcut_(0_GeV) - , energy_invcut_(0_GeV) - , em_count_(0) - , inv_count_(0) - , energy_count_(0) { + template <typename TOutput> + template <typename... TArgs> + inline ParticleCut<TOutput>::ParticleCut( + std::unordered_map<Code const, HEPEnergyType const> const& eCuts, bool const inv, + TArgs&&... args) + : TOutput(std::forward<TArgs>(args)...) + , doCutInv_(inv) { set_kinetic_energy_thresholds(eCuts); CORSIKA_LOG_DEBUG("setting threshold particles individually"); } + template <typename TOutput> template <typename TParticle> - inline bool ParticleCut::isBelowEnergyCut(TParticle const& vP) const { + inline bool ParticleCut<TOutput>::isBelowEnergyCut(TParticle const& vP) const { auto const energyLab = vP.getKineticEnergy(); auto const pid = vP.getPID(); // nuclei if (is_nucleus(pid)) { // calculate energy per nucleon auto const ElabNuc = energyLab / get_nucleus_A(pid); - return (ElabNuc < calculate_kinetic_energy_threshold(pid)); + return (ElabNuc < get_kinetic_energy_threshold(pid)); } else { - return (energyLab < calculate_kinetic_energy_threshold(pid)); + return (energyLab < get_kinetic_energy_threshold(pid)); } } - inline bool ParticleCut::isInvisible(Code const& vCode) const { - return is_neutrino(vCode); - } - + template <typename TOutput> template <typename TParticle> - inline bool ParticleCut::checkCutParticle(TParticle const& particle) { + inline bool ParticleCut<TOutput>::checkCutParticle(TParticle const& particle) { Code const pid = particle.getPID(); HEPEnergyType const kine_energy = particle.getKineticEnergy(); HEPEnergyType const energy = particle.getEnergy(); CORSIKA_LOG_DEBUG( - "ParticleCut: checking {} ({}), E_kin= {} GeV, EcutTot={} GeV, E={} GeV, m={} " + "ParticleCut: checking {} ({}), E_kin= {} GeV, E={} GeV, m={} " "GeV", - pid, particle.getPDG(), kine_energy / 1_GeV, - (energy_emcut_ + energy_invcut_ + energy_cut_ + energy_timecut_) / 1_GeV, - energy / 1_GeV, particle.getMass() / 1_GeV); + pid, particle.getPDG(), kine_energy / 1_GeV, energy / 1_GeV, + particle.getMass() / 1_GeV); CORSIKA_LOG_DEBUG("p={}", particle.asString()); - if (doCutEm_ && is_em(pid)) { - CORSIKA_LOG_DEBUG("removing em. particle..."); - energy_emcut_ += kine_energy; - em_count_ += 1; - energy_event_ += kine_energy; - return true; - } else if (doCutInv_ && is_neutrino(pid)) { + if (doCutInv_ && is_neutrino(pid)) { CORSIKA_LOG_DEBUG("removing inv. particle..."); - energy_invcut_ += kine_energy; - inv_count_ += 1; - energy_event_ += kine_energy; return true; } else if (isBelowEnergyCut(particle)) { CORSIKA_LOG_DEBUG("removing low en. particle..."); - energy_cut_ += kine_energy; - energy_count_ += 1; - energy_event_ += kine_energy; return true; } else if (particle.getTime() > 10_ms) { CORSIKA_LOG_DEBUG("removing OLD particle..."); - energy_timecut_ += kine_energy; - energy_event_ += kine_energy; return true; + } else { + for (auto const& cut : cuts_) { + if (pid == cut.first && kine_energy < cut.second) { return true; } + } } return false; // this particle will not be removed/cut } + template <typename TOutput> template <typename TStackView> - inline void ParticleCut::doSecondaries(TStackView& vS) { - energy_event_ = 0_GeV; // per event counting for printout + inline void ParticleCut<TOutput>::doSecondaries(TStackView& vS) { + HEPEnergyType energy_event = 0_GeV; // per event counting for printout auto particle = vS.begin(); while (particle != vS.end()) { - if (checkCutParticle(particle)) { particle.erase(); } + if (checkCutParticle(particle)) { + HEPEnergyType Ekin = particle.getKineticEnergy(); + this->write(particle.getPosition(), particle.getPID(), Ekin); + particle.erase(); + } ++particle; // next entry in SecondaryView } - CORSIKA_LOG_DEBUG("Event cut: {} GeV", energy_event_ / 1_GeV); + CORSIKA_LOG_DEBUG("Event cut: {} GeV", energy_event / 1_GeV); } + template <typename TOutput> template <typename TParticle, typename TTrajectory> - inline ProcessReturn ParticleCut::doContinuous(TParticle& particle, TTrajectory const&, - bool const) { - CORSIKA_LOG_TRACE("ParticleCut::DoContinuous"); + inline ProcessReturn ParticleCut<TOutput>::doContinuous(TParticle& particle, + TTrajectory const&, + bool const) { if (checkCutParticle(particle)) { + this->write(particle.getPosition(), particle.getPID(), particle.getKineticEnergy()); CORSIKA_LOG_TRACE("removing during continuous"); // signal to upstream code that this particle was deleted return ProcessReturn::ParticleAbsorbed; @@ -181,34 +140,38 @@ namespace corsika { return ProcessReturn::Ok; } - inline void ParticleCut::printThresholds() { - for (auto p : get_all_particles()) { - auto const Eth = calculate_kinetic_energy_threshold(p); - CORSIKA_LOG_INFO("kinetic energy threshold for particle {} is {} GeV", p, - Eth / 1_GeV); - } - } + template <typename TOutput> + inline void ParticleCut<TOutput>::printThresholds() const { + + CORSIKA_LOG_DEBUG("kinetic energy threshold for electrons is {} GeV", + cut_electrons_ / 1_GeV); + CORSIKA_LOG_DEBUG("kinetic energy threshold for photons is {} GeV", + cut_photons_ / 1_GeV); + CORSIKA_LOG_DEBUG("kinetic energy threshold for muons is {} GeV", cut_muons_ / 1_GeV); + CORSIKA_LOG_DEBUG("kinetic energy threshold for hadros is {} GeV", + cut_hadrons_ / 1_GeV); - inline void ParticleCut::showResults() { - CORSIKA_LOG_INFO( - "\n ******************************\n " - " kinetic energy removed by cut of electromagnetic (GeV): {} (number: {})\n " - " kinetic energy removed by cut of invisible (GeV): {} (number: {})\n " - " kinetic energy removed by kinetic energy cut (GeV): {} (number: {}) \n " - " kinetic energy removed by time cut (GeV): {} \n" - " ******************************", - energy_emcut_ / 1_GeV, em_count_, energy_invcut_ / 1_GeV, inv_count_, - energy_cut_ / 1_GeV, energy_count_, energy_timecut_ / 1_GeV); + for (auto const& cut : cuts_) { + CORSIKA_LOG_DEBUG("kinetic energy threshold for particle {} is {} GeV", cut.first, + cut.second / 1_GeV); + } } - inline void ParticleCut::reset() { - energy_emcut_ = 0_GeV; - em_count_ = 0; - energy_invcut_ = 0_GeV; - inv_count_ = 0; - energy_cut_ = 0_GeV; - energy_count_ = 0; - energy_timecut_ = 0_GeV; + template <typename TOutput> + inline YAML::Node ParticleCut<TOutput>::getConfig() const { + + YAML::Node node; + node["type"] = "ParticleCut"; + node["units"]["energy"] = "GeV"; + node["cut_electrons"] = cut_electrons_ / 1_GeV; + node["cut_photons"] = cut_photons_ / 1_GeV; + node["cut_muons"] = cut_muons_ / 1_GeV; + node["cut_hadrons"] = cut_hadrons_ / 1_GeV; + node["cut_invisibles"] = doCutInv_; + for (auto const& cut : cuts_) { + node[fmt::format("cut_{}", cut.first)] = cut.second / 1_GeV; + } + return node; } } // namespace corsika diff --git a/corsika/detail/modules/StackInspector.inl b/corsika/detail/modules/StackInspector.inl index 06705ac8f6be0892750e8c3942e4710d273b6dc8..f801abf624705b5bb5891bd34304c3ed0570936b 100644 --- a/corsika/detail/modules/StackInspector.inl +++ b/corsika/detail/modules/StackInspector.inl @@ -54,20 +54,24 @@ namespace corsika { } } - auto const now = std::chrono::system_clock::now(); - std::chrono::duration<double> const elapsed_seconds = now - StartTime_; - std::time_t const now_time = std::chrono::system_clock::to_time_t(now); + std::chrono::system_clock::time_point const now = std::chrono::system_clock::now(); + std::chrono::duration<double> const elapsed_seconds = now - StartTime_; // seconds auto const dE = E0_ - Etot; if (dE < dE_threshold_) return; double const progress = dE / E0_; + // for printout + std::time_t const now_time = std::chrono::system_clock::to_time_t(now); std::time_t const start_time = std::chrono::system_clock::to_time_t(StartTime_); if (progress > 0) { double const eta_seconds = elapsed_seconds.count() / progress; - std::time_t const eta_time = std::chrono::system_clock::to_time_t( - StartTime_ + std::chrono::seconds((int)eta_seconds)); + std::chrono::system_clock::time_point const eta = + StartTime_ + std::chrono::seconds((int)eta_seconds); + + // for printout + std::time_t const eta_time = std::chrono::system_clock::to_time_t(eta); int const yday0 = std::localtime(&start_time)->tm_yday; int const yday1 = std::localtime(&eta_time)->tm_yday; diff --git a/corsika/detail/modules/TrackWriter.inl b/corsika/detail/modules/TrackWriter.inl index a1102ce2cf7ee18de3902baff0c141617dc099f7..39bb395acd18eb4d4b46f80f17833a89a47ebf3b 100644 --- a/corsika/detail/modules/TrackWriter.inl +++ b/corsika/detail/modules/TrackWriter.inl @@ -15,41 +15,41 @@ namespace corsika { - template <typename TOutputWriter> - inline TrackWriter<TOutputWriter>::TrackWriter() {} + template <typename TOutput> + inline TrackWriter<TOutput>::TrackWriter() {} - template <typename TOutputWriter> + template <typename TOutput> template <typename TParticle, typename TTrack> - inline ProcessReturn TrackWriter<TOutputWriter>::doContinuous(TParticle const& vP, - TTrack const& vT, - bool const) { + inline ProcessReturn TrackWriter<TOutput>::doContinuous(TParticle const& vP, + TTrack const& vT, bool const) { auto const start = vT.getPosition(0).getCoordinates(); auto const end = vT.getPosition(1).getCoordinates(); // write the track to the file - this->write(vP.getPID(), vP.getEnergy(), start, vP.getTime() - vT.getDuration(), end, - vP.getTime()); + TOutput::write(vP.getPID(), vP.getEnergy(), vP.getWeight(), start, + vP.getTime() - vT.getDuration(), end, vP.getTime()); return ProcessReturn::Ok; } - template <typename TOutputWriter> + template <typename TOutput> template <typename TParticle, typename TTrack> - inline LengthType TrackWriter<TOutputWriter>::getMaxStepLength(TParticle const&, - TTrack const&) { + inline LengthType TrackWriter<TOutput>::getMaxStepLength(TParticle const&, + TTrack const&) { return meter * std::numeric_limits<double>::infinity(); } - template <typename TOutputWriter> - YAML::Node TrackWriter<TOutputWriter>::getConfig() const { + template <typename TOutput> + YAML::Node TrackWriter<TOutput>::getConfig() const { using namespace units::si; YAML::Node node; // add default units for values node["type"] = "TrackWriter"; - node["units"] = "GeV | m | s"; + node["units"]["energy"] = "GeV"; + node["units"]["length"] = "m"; return node; } diff --git a/corsika/detail/modules/conex/CONEXhybrid.inl b/corsika/detail/modules/conex/CONEXhybrid.inl index 00e79005e236ef415bbdc5174d0615dec90f5049..36f5f93cee38427c58970af0455ca076266f8b3f 100644 --- a/corsika/detail/modules/conex/CONEXhybrid.inl +++ b/corsika/detail/modules/conex/CONEXhybrid.inl @@ -22,10 +22,14 @@ namespace corsika { - inline CONEXhybrid::CONEXhybrid(Point center, ShowerAxis const& showerAxis, - LengthType groundDist, LengthType injectionHeight, - HEPEnergyType primaryEnergy, PDGCode primaryPDG) - : center_{center} + template <typename TOutputE, typename TOutputN> + inline CONEXhybrid<TOutputE, TOutputN>::CONEXhybrid( + Point const& center, ShowerAxis const& showerAxis, LengthType groundDist, + LengthType injectionHeight, HEPEnergyType primaryEnergy, PDGCode primaryPDG, + TOutputE& args1, TOutputN& args2) + : SubWriter<TOutputE>(args1) + , SubWriter<TOutputN>(args2) + , center_{center} , showerAxis_{showerAxis} , groundDist_{groundDist} , injectionHeight_{injectionHeight} @@ -57,8 +61,7 @@ namespace corsika { return b.normalized(); })} - , y_sf_{showerAxis_.getDirection().cross(x_sf_)} - , energy_em_(0_GeV) { + , y_sf_{showerAxis_.getDirection().cross(x_sf_)} { CORSIKA_LOG_DEBUG("x_sf (conexObservationCS): {}", x_sf_.getComponents(conexObservationCS_)); @@ -99,7 +102,8 @@ namespace corsika { configPath.c_str(), configPath.size()); } - inline void CONEXhybrid::initCascadeEquations() { + template <typename TOutputE, typename TOutputN> + inline void CONEXhybrid<TOutputE, TOutputN>::initCascadeEquations() { // set phi, theta Vector<length_d> ez{conexObservationCS_, {0._m, 0._m, -1_m}}; @@ -132,8 +136,9 @@ namespace corsika { ::conex::conexrun_(ipart, eprima, theta, phi, xminp, dimpact, ioseed.data()); } + template <typename TOutputE, typename TOutputN> template <typename TStackView> - inline void CONEXhybrid::doSecondaries(TStackView& vS) { + inline void CONEXhybrid<TOutputE, TOutputN>::doSecondaries(TStackView& vS) { auto p = vS.begin(); while (p != vS.end()) { Code const pid = p.getPID(); @@ -145,9 +150,10 @@ namespace corsika { } } - inline bool CONEXhybrid::addParticle(Code pid, HEPEnergyType energy, HEPEnergyType mass, - Point const& position, - DirectionVector const& direction, TimeType t) { + template <typename TOutputE, typename TOutputN> + inline bool CONEXhybrid<TOutputE, TOutputN>::addParticle( + Code pid, HEPEnergyType energy, HEPEnergyType mass, Point const& position, + DirectionVector const& direction, TimeType t, double weight) { auto const it = std::find_if(egs_em_codes_.cbegin(), egs_em_codes_.cend(), [=](auto const& p) { return pid == p.first; }); @@ -182,14 +188,11 @@ namespace corsika { double const v = direction.dot(x_sf_).magnitude(); double const w = direction.dot(showerAxis_.getDirection()).magnitude(); - double const weight = 1; // NEEDS TO BE CHANGED WHEN WE HAVE WEIGHTS! - // generation, TO BE CHANGED WHEN WE HAVE THAT INFORMATION AVAILABLE int const latchin = 1; double const E = energy / 1_GeV; double const m = mass / 1_GeV; - energy_em_ += energy; CORSIKA_LOG_DEBUG("CONEXhybrid: removing {} {:5e} GeV", egs_pid, energy); @@ -232,8 +235,9 @@ namespace corsika { return true; } + template <typename TOutputE, typename TOutputN> template <typename TStack> - inline void CONEXhybrid::doCascadeEquations(TStack&) { + inline void CONEXhybrid<TOutputE, TOutputN>::doCascadeEquations(TStack&) { ::conex::conexcascade_(); @@ -268,33 +272,25 @@ namespace corsika { ::conex::get_shower_electron_(icute, nX, Electrons[0]); ::conex::get_shower_hadron_(icuth, nX, Hadrons[0]); - std::ofstream file{"conex_output.txt"}; - file << fmt::format("#{:>10} {:>13} {:>13} {:>13} {:>13} {:>13} {:>13} {:>13}\n", "X", - "N", "dEdX", "Mu", "dMu", "Photon", "El", "Had"); + // make sure CONEX binning is same to C8: + GrammageType dX = (X[1] - X[0]) * 1_g / square(1_cm); + for (int i = 0; i < nX; ++i) { - file << fmt::format(" {:>10.2f} {:.5e} {:.5e} {:.5e} {:.5e} {:.5e} {:.5e} {:.5e}\n", - X[i], N[i], dEdX[i], Mu[i], dMu[i], Photon[i], Electrons[i], - Hadrons[i]); + GrammageType curX = X[i] * 1_g / square(1_cm); + SubWriter<TOutputE>::write(curX, curX + dX, + Code::Unknown, // this is sum of all dEdX + dEdX[i] * 1_GeV / 1_g * square(1_cm) * dX); + SubWriter<TOutputN>::write(curX, curX + dX, Code::Photon, Photon[i]); + SubWriter<TOutputN>::write(curX, curX + dX, Code::Proton, Hadrons[i]); + SubWriter<TOutputN>::write(curX, curX + dX, Code::Electron, Electrons[i]); + SubWriter<TOutputN>::write(curX, curX + dX, Code::MuMinus, Mu[i]); } - - std::ofstream fitout{"conex_fit.txt"}; - fitout << fitpars[1 - 1] << " # log10(eprima/eV)" << std::endl; - fitout << fitpars[2 - 1] << " # theta" << std::endl; - fitout << fitpars[3 - 1] << " # X1 (first interaction)" << std::endl; - fitout << fitpars[4 - 1] << " # Nmax" << std::endl; - fitout << fitpars[5 - 1] << " # X0" << std::endl; - fitout << fitpars[6 - 1] << " # P1" << std::endl; - fitout << fitpars[7 - 1] << " # P2" << std::endl; - fitout << fitpars[8 - 1] << " # P3" << std::endl; - fitout << fitpars[9 - 1] << " # chi^2 / sqrt(Nmax)" << std::endl; - fitout << fitpars[10 - 1] << " # Xmax" << std::endl; - fitout << fitpars[11 - 1] << " # phi" << std::endl; - fitout << fitpars[12 - 1] << " # inelasticity 1st int." << std::endl; - fitout << fitpars[13 - 1] << " # ???" << std::endl; } - inline HEPEnergyType CONEXhybrid::getEnergyEM() const { return energy_em_; } + template <typename TOutputE, typename TOutputN> + inline YAML::Node CONEXhybrid<TOutputE, TOutputN>::getConfig() const { - inline void CONEXhybrid::reset() { energy_em_ = 0_GeV; } + return YAML::Node(); + } } // namespace corsika diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl index 066997364c8d7843cdf912e893294b3eada2a825..d163f761a8af34c2f5ef7d6156db8af005209c09 100644 --- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl +++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl @@ -19,19 +19,15 @@ namespace corsika { - auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { - return sqrt((Elab - m) * (Elab + m)); - }; - - inline BetheBlochPDG::BetheBlochPDG(ShowerAxis const& shower_axis) - : dX_(10_g / square(1_cm)) // profile binning - , dX_threshold_(0.0001_g / square(1_cm)) - , shower_axis_(shower_axis) - , profile_(int(shower_axis.getMaximumX() / dX_) + 1) {} + template <typename TOutput> + template <typename... TArgs> + inline BetheBlochPDG<TOutput>::BetheBlochPDG(TArgs&&... args) + : TOutput(std::forward<TArgs>(args)...) {} + template <typename TOutput> template <typename TParticle> - inline HEPEnergyType BetheBlochPDG::getBetheBloch(TParticle const& p, - GrammageType const dX) { + inline HEPEnergyType BetheBlochPDG<TOutput>::getBetheBloch(TParticle const& p, + GrammageType const dX) { // all these are material constants and have to come through Environment // right now: values for nitrogen_D @@ -124,24 +120,28 @@ namespace corsika { } // radiation losses according to PDG 2018, ch. 33 ref. [5] + template <typename TOutput> template <typename TParticle> - inline HEPEnergyType BetheBlochPDG::getRadiationLosses(TParticle const& vP, - GrammageType const vDX) { + inline HEPEnergyType BetheBlochPDG<TOutput>::getRadiationLosses( + TParticle const& vP, GrammageType const vDX) { // simple-minded hard-coded value for b(E) inspired by data from // http://pdg.lbl.gov/2018/AtomicNuclearProperties/ for N and O. auto constexpr b = 3.0 * 1e-6 * square(1_cm) / 1_g; return -vP.getEnergy() * b * vDX; } + template <typename TOutput> template <typename TParticle> - inline HEPEnergyType BetheBlochPDG::getTotalEnergyLoss(TParticle const& vP, - GrammageType const vDX) { + inline HEPEnergyType BetheBlochPDG<TOutput>::getTotalEnergyLoss( + TParticle const& vP, GrammageType const vDX) { return getBetheBloch(vP, vDX) + getRadiationLosses(vP, vDX); } + template <typename TOutput> template <typename TParticle, typename TTrajectory> - inline ProcessReturn BetheBlochPDG::doContinuous(TParticle& p, TTrajectory const& t, - bool const) { + inline ProcessReturn BetheBlochPDG<TOutput>::doContinuous(TParticle& particle, + TTrajectory const& track, + bool const) { // if this step was limiting the CORSIKA stepping, the particle is lost /* see Issue https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/389 @@ -152,25 +152,28 @@ namespace corsika { } */ - if (p.getChargeNumber() == 0) return ProcessReturn::Ok; - - GrammageType const dX = p.getNode()->getModelProperties().getIntegratedGrammage(t); - CORSIKA_LOG_TRACE("EnergyLoss pid={}, z={}, dX={} g/cm2", p.getPID(), - p.getChargeNumber(), dX / 1_g * square(1_cm)); - HEPEnergyType dE = getTotalEnergyLoss(p, dX); - auto E = p.getEnergy(); - [[maybe_unused]] const auto Ekin = E - p.getMass(); - auto Enew = E + dE; - CORSIKA_LOG_TRACE("EnergyLoss dE={} MeV, E={} GeV, Ekin={} GeV, Enew={} GeV", - dE / 1_MeV, E / 1_GeV, Ekin / 1_GeV, Enew / 1_GeV); - p.setEnergy(Enew); // kinetic energy on stack - fillProfile(t, dE); + if (particle.getChargeNumber() == 0) return ProcessReturn::Ok; + + GrammageType const dX = + particle.getNode()->getModelProperties().getIntegratedGrammage(track); + CORSIKA_LOG_TRACE("EnergyLoss pid={}, z={}, dX={} g/cm2", particle.getPID(), + particle.getChargeNumber(), dX / 1_g * square(1_cm)); + HEPEnergyType const dE = getTotalEnergyLoss(particle, dX); + [[maybe_unused]] const auto Ekin = particle.getKineticEnergy(); + auto EkinNew = Ekin + dE; + CORSIKA_LOG_TRACE("EnergyLoss dE={} MeV, Ekin={} GeV, EkinNew={} GeV", dE / 1_MeV, + Ekin / 1_GeV, EkinNew / 1_GeV); + particle.setKineticEnergy(EkinNew); + + // also send to output + TOutput::write(track, particle.getPID(), -dE); return ProcessReturn::Ok; } + template <typename TOutput> template <typename TParticle, typename TTrajectory> - inline LengthType BetheBlochPDG::getMaxStepLength(TParticle const& vParticle, - TTrajectory const& vTrack) const { + inline LengthType BetheBlochPDG<TOutput>::getMaxStepLength( + TParticle const& vParticle, TTrajectory const& vTrack) const { if (vParticle.getChargeNumber() == 0) { return meter * std::numeric_limits<double>::infinity(); } @@ -180,7 +183,7 @@ namespace corsika { auto const energy = vParticle.getKineticEnergy(); auto const energy_lim = std::max(energy * 0.9, // either 10% relative loss max., or - calculate_kinetic_energy_threshold( + get_kinetic_energy_threshold( vParticle.getPID()) // energy thresholds globally defined // for individual particles * 0.99999 // need to go slightly below global e-cut to assure @@ -193,80 +196,12 @@ namespace corsika { vTrack, maxGrammage); } - template <typename TTrajectory> - inline void BetheBlochPDG::fillProfile(TTrajectory const& vTrack, - const HEPEnergyType dE) { - - GrammageType grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0)); - GrammageType grammageEnd = shower_axis_.getProjectedX(vTrack.getPosition(1)); - - if (grammageStart > grammageEnd) { // particle going upstream - std::swap(grammageStart, grammageEnd); - } - - GrammageType const deltaX = grammageEnd - grammageStart; - - if (deltaX < dX_threshold_) return; - - // only register the range that is covered by the profile - int const maxBin = int(profile_.size() - 1); - int binStart = grammageStart / dX_; - if (binStart < 0) binStart = 0; - if (binStart > maxBin) binStart = maxBin; - int binEnd = grammageEnd / dX_; - if (binEnd < 0) binEnd = 0; - if (binEnd > maxBin) binEnd = maxBin; - - CORSIKA_LOG_TRACE("energy deposit of -dE={} between {} and {}", -dE, grammageStart, - grammageEnd); - - auto energyCount = HEPEnergyType::zero(); - - auto const factor = -dE / deltaX; - auto fill = [&](int const bin, GrammageType const weight) { - auto const increment = factor * weight; - profile_[bin] += increment; - energyCount += increment; - - CORSIKA_LOG_TRACE("filling bin {} with weight {} : {} ", bin, weight, increment); - }; - - // fill longitudinal profile - if (binStart == binEnd) { - fill(binStart, deltaX); - } else { - fill(binStart, ((1 + binStart) * dX_ - grammageStart)); - fill(binEnd, (grammageEnd - binEnd * dX_)); - for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); } - } - - CORSIKA_LOG_TRACE("total energy added to histogram: {} ", energyCount); - } - - inline void BetheBlochPDG::printProfile() const { - std::ofstream file("EnergyLossProfile.dat"); - file << "# EnergyLoss profile" << std::endl - << "# lower X bin edge [g/cm2] dE/dX [GeV/g/cm2]\n"; - double const deltaX = dX_ / 1_g * square(1_cm); - for (size_t i = 0; i < profile_.size(); ++i) { - file << std::scientific << std::setw(15) << i * deltaX << std::setw(15) - << profile_.at(i) / (deltaX * 1_GeV) << '\n'; - } - file.close(); - } + template <typename TOutput> + inline YAML::Node BetheBlochPDG<TOutput>::getConfig() const { - inline HEPEnergyType BetheBlochPDG::getTotal() const { - return std::accumulate(profile_.cbegin(), profile_.cend(), HEPEnergyType::zero()); + YAML::Node node; + node["type"] = "BetheBlochPDG"; + return node; } - inline void BetheBlochPDG::showResults() const { - CORSIKA_LOG_INFO( - " ******************************\n" - " PROCESS::ContinuousProcess: \n" - " energy lost dE (GeV) : {}\n ", - energy_lost_ / 1_GeV); - } - - inline void BetheBlochPDG::reset() { energy_lost_ = 0_GeV; } - } // namespace corsika diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl index 89f571b69d620a3e893e6b14ebe6b48913f02bd5..77f77e7795a9268e58da251e1b851aad2d0e3d2e 100644 --- a/corsika/detail/modules/proposal/ContinuousProcess.inl +++ b/corsika/detail/modules/proposal/ContinuousProcess.inl @@ -18,8 +18,9 @@ namespace corsika::proposal { - inline void ContinuousProcess::buildCalculator(Code code, - NuclearComposition const& comp) { + template <typename TOutput> + inline void ContinuousProcess<TOutput>::buildCalculator( + Code code, NuclearComposition const& comp) { // search crosssection builder for given particle auto p_cross = cross.find(code); if (p_cross == cross.end()) @@ -29,7 +30,7 @@ namespace corsika::proposal { // take some minutes if you have to build the tables and cannot read the // from disk auto const emCut = - calculate_kinetic_energy_threshold(code) + + get_kinetic_energy_threshold(code) + get_mass(code); //! energy thresholds globally defined for individual particles auto c = p_cross->second(media.at(comp.getHash()), emCut); @@ -48,24 +49,29 @@ namespace corsika::proposal { calc[std::make_pair(comp.getHash(), code)] = std::move(calculator); } - template <typename TEnvironment> - inline ContinuousProcess::ContinuousProcess(TEnvironment const& _env) - : ProposalProcessBase(_env) {} + template <typename TOutput> + template <typename TEnvironment, typename... TOutputArgs> + inline ContinuousProcess<TOutput>::ContinuousProcess(TEnvironment const& _env, + TOutputArgs&&... args) + : TOutput(args...) + , ProposalProcessBase(_env) {} + template <typename TOutput> template <typename TParticle> - inline void ContinuousProcess::scatter(TParticle& vP, HEPEnergyType const& loss, - GrammageType const& grammage) { + inline void ContinuousProcess<TOutput>::scatter(TParticle& particle, + HEPEnergyType const& loss, + GrammageType const& grammage) { // get or build corresponding calculators - auto c = getCalculator(vP, calc); + auto c = getCalculator(particle, calc); // Cast corsika vector to proposal vector - auto vP_dir = vP.getDirection(); - auto d = vP_dir.getComponents(); + auto particle_dir = particle.getDirection(); + auto d = particle_dir.getComponents(); auto direction = PROPOSAL::Cartesian3D(d.getX().magnitude(), d.getY().magnitude(), d.getZ().magnitude()); - auto E_f = vP.getEnergy() - loss; + auto E_f = particle.getEnergy() - loss; // draw random numbers required for scattering process std::uniform_real_distribution<double> distr(0., 1.); @@ -74,27 +80,28 @@ namespace corsika::proposal { // calculate deflection based on particle energy, loss auto deflection = (c->second).scatter->CalculateMultipleScattering( - grammage / 1_g * square(1_cm), vP.getEnergy() / 1_MeV, E_f / 1_MeV, rnd); + grammage / 1_g * square(1_cm), particle.getEnergy() / 1_MeV, E_f / 1_MeV, rnd); [[maybe_unused]] auto [unused1, final_direction] = PROPOSAL::multiple_scattering::ScatterInitialDirection(direction, deflection); // update particle direction after continuous loss caused by multiple // scattering - vP.setDirection( - {vP_dir.getCoordinateSystem(), + particle.setDirection( + {particle_dir.getCoordinateSystem(), {final_direction.GetX(), final_direction.GetY(), final_direction.GetZ()}}); } + template <typename TOutput> template <typename TParticle, typename TTrajectory> - inline ProcessReturn ContinuousProcess::doContinuous(TParticle& vP, - TTrajectory const& vT, - bool const) { + inline ProcessReturn ContinuousProcess<TOutput>::doContinuous(TParticle& vP, + TTrajectory const& track, + bool const) { if (!canInteract(vP.getPID())) return ProcessReturn::Ok; - if (vT.getLength() == 0_m) return ProcessReturn::Ok; + if (track.getLength() == 0_m) return ProcessReturn::Ok; // calculate passed grammage - auto dX = vP.getNode()->getModelProperties().getIntegratedGrammage(vT); + auto dX = vP.getNode()->getModelProperties().getIntegratedGrammage(track); // get or build corresponding track integral calculator and solve the // integral @@ -103,17 +110,21 @@ namespace corsika::proposal { vP.getEnergy() / 1_MeV, dX / 1_g * 1_cm * 1_cm) * 1_MeV; auto dE = vP.getEnergy() - final_energy; - energy_lost_ += dE; // if the particle has a charge take multiple scattering into account if (vP.getChargeNumber() != 0) scatter(vP, dE, dX); vP.setEnergy(final_energy); // on the stack, this is just kinetic energy, E-m + + // also send to output + TOutput::write(track, vP.getPID(), dE); + return ProcessReturn::Ok; } + template <typename TOutput> template <typename TParticle, typename TTrajectory> - inline LengthType ContinuousProcess::getMaxStepLength(TParticle const& vP, - TTrajectory const& vT) { + inline LengthType ContinuousProcess<TOutput>::getMaxStepLength( + TParticle const& vP, TTrajectory const& track) { auto const code = vP.getPID(); if (!canInteract(code)) return meter * std::numeric_limits<double>::infinity(); @@ -123,7 +134,7 @@ namespace corsika::proposal { auto const energy = vP.getEnergy(); auto const energy_lim = std::max(energy * 0.9, // either 10% relative loss max., or - calculate_kinetic_energy_threshold( + get_kinetic_energy_threshold( code) // energy thresholds globally defined for individual particles * 0.9999 // need to go slightly below global e-cut to assure removal // in ParticleCut. This does not matter since at cut-time @@ -137,20 +148,16 @@ namespace corsika::proposal { square(1_cm); // return it in distance aequivalent - auto dist = vP.getNode()->getModelProperties().getArclengthFromGrammage(vT, grammage); + auto dist = + vP.getNode()->getModelProperties().getArclengthFromGrammage(track, grammage); CORSIKA_LOG_TRACE("PROPOSAL::getMaxStepLength X={} g/cm2, l={} m ", grammage / 1_g * square(1_cm), dist / 1_m); return dist; } - inline void ContinuousProcess::showResults() const { - CORSIKA_LOG_DEBUG( - " ******************************\n" - " PROCESS::ContinuousProcess: \n" - " energy lost dE (GeV) : {}", - energy_lost_ / 1_GeV); + template <typename TOutput> + inline YAML::Node ContinuousProcess<TOutput>::getConfig() const { + return YAML::Node(); } - inline void ContinuousProcess::reset() { energy_lost_ = 0_GeV; } - } // namespace corsika::proposal diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl index af6d603f6874113c4430c157085557956ca19581..d8916cbe976e77b8865d7d0edaacc6bfc0135e0b 100644 --- a/corsika/detail/modules/proposal/Interaction.inl +++ b/corsika/detail/modules/proposal/Interaction.inl @@ -33,7 +33,7 @@ namespace corsika::proposal { // take some minutes if you have to build the tables and cannot read the // from disk auto const emCut = - calculate_kinetic_energy_threshold(code) + + get_kinetic_energy_threshold(code) + get_mass(code); //! energy thresholds globally defined for individual particles auto c = p_cross->second(media.at(comp.getHash()), emCut); diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl index 9e070f1c62df8860f9bbb4d89df26f9c87728e09..1772df106dcb04d8744eac83959e833f5fa27f9a 100644 --- a/corsika/detail/modules/urqmd/UrQMD.inl +++ b/corsika/detail/modules/urqmd/UrQMD.inl @@ -13,6 +13,7 @@ #include <corsika/framework/core/ParticleProperties.hpp> #include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/framework/core/EnergyMomentumOperations.hpp> #include <corsika/framework/geometry/QuantityVector.hpp> #include <corsika/framework/geometry/Vector.hpp> #include <corsika/framework/utility/COMBoost.hpp> @@ -280,9 +281,11 @@ namespace corsika::urqmd { CORSIKA_LOG_DEBUG(" {} {} {} ", i, code, momentum.getComponents()); HEPEnergyType const mass = get_mass(code); - HEPEnergyType Ekin = sqrt(momentum.getSquaredNorm() + mass * mass) - mass; + HEPEnergyType const Ekin = calculate_kinetic_energy(momentum.getNorm(), mass); if (Ekin <= 0_GeV) { - CORSIKA_LOG_WARN("Negative kinetic energy {} {}. Skipping.", code, Ekin); + if (Ekin < 0_GeV) { + CORSIKA_LOG_WARN("Negative kinetic energy {} {}. Skipping.", code, Ekin); + } view.addSecondary( std::make_tuple(code, 0_eV, DirectionVector{originalCS, {0, 0, 0}})); } else { diff --git a/corsika/detail/modules/writers/EnergyLossWriter.inl b/corsika/detail/modules/writers/EnergyLossWriter.inl new file mode 100644 index 0000000000000000000000000000000000000000..ddf79f372c4971721d218aa7902f3d32151feeff --- /dev/null +++ b/corsika/detail/modules/writers/EnergyLossWriter.inl @@ -0,0 +1,233 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/framework/utility/FindXmax.hpp> + +#include <corsika/media/ShowerAxis.hpp> + +#include <exception> + +namespace corsika { + + template <typename TOutput> + inline EnergyLossWriter<TOutput>::EnergyLossWriter(ShowerAxis const& axis, + GrammageType dX, + unsigned int const nBins, + GrammageType dX_threshold) + : TOutput(dEdX_output::ProfileIndexNames) + , showerAxis_(axis) + , dX_(dX) + , nBins_(nBins) + , dX_threshold_(dX_threshold) {} + + template <typename TOutput> + inline void EnergyLossWriter<TOutput>::startOfLibrary( + boost::filesystem::path const& directory) { + TOutput::startOfLibrary(directory); + } + + template <typename TOutput> + inline void EnergyLossWriter<TOutput>::startOfShower(unsigned int const showerId) { + TOutput::startOfShower(showerId); + // reset profile + profile_.clear(); + profile_.resize(nBins_); + } + + template <typename TOutput> + inline void EnergyLossWriter<TOutput>::endOfLibrary() { + TOutput::endOfLibrary(); + } + + template <typename TOutput> + inline void EnergyLossWriter<TOutput>::endOfShower(unsigned int const showerId) { + + int iRow{0}; + for (dEdX_output::Profile const& row : profile_) { + // here: write to underlying writer (e.g. parquet) + TOutput::write(showerId, iRow * dX_, row); + iRow++; + } + + TOutput::endOfShower(showerId); + } + + template <typename TOutput> + template <typename TTrack> + inline void EnergyLossWriter<TOutput>::write(TTrack const& track, Code const PID, + HEPEnergyType const dE) { + + GrammageType grammageStart = showerAxis_.getProjectedX(track.getPosition(0)); + GrammageType grammageEnd = showerAxis_.getProjectedX(track.getPosition(1)); + + if (grammageStart > grammageEnd) { // particle going upstream + std::swap(grammageStart, grammageEnd); + } + + GrammageType const deltaX = grammageEnd - grammageStart; + + CORSIKA_LOGGER_TRACE( + TOutput::getLogger(), + "dE={} GeV, grammageStart={} g/cm2, End={}g /cm2, deltaX={} g/cm2", dE / 1_GeV, + grammageStart / 1_g * square(1_cm), grammageEnd / 1_g * square(1_cm), + deltaX / 1_g * square(1_cm)); + + if (deltaX < dX_threshold_) { + CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "Point-like dE"); + this->write(track.getPosition(0), PID, dE); + return; + } + + // only register the range that is covered by the profile + int const maxBin = int(profile_.size() - 1); + int binStart = grammageStart / dX_; + if (binStart < 0) binStart = 0; + if (binStart > maxBin) binStart = maxBin; + int binEnd = grammageEnd / dX_; + if (binEnd < 0) binEnd = 0; + if (binEnd > maxBin) binEnd = maxBin; + + CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "maxBin={}, binStart={}, binEnd={}", + maxBin, binStart, binEnd); + + auto energyCount = HEPEnergyType::zero(); + + auto const factor = dE / deltaX; // [ energy / grammage ] + auto fill = [&](int const bin, GrammageType const weight) { + auto const increment = factor * weight; + CORSIKA_LOGGER_TRACE(TOutput::getLogger(), + "filling bin={} with weight {} : dE={} GeV ", bin, weight, + increment / 1_GeV); + profile_[bin][static_cast<int>(dEdX_output::ProfileIndex::Total)] += increment; + energyCount += increment; + }; + + // fill longitudinal profile + if (binStart == binEnd) { + fill(binStart, deltaX); + } else { + fill(binStart, ((1 + binStart) * dX_ - grammageStart)); + fill(binEnd, (grammageEnd - binEnd * dX_)); + for (int bin = binStart + 1; bin < binEnd; ++bin) { fill(bin, dX_); } + } + + CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "total energy added to histogram: {} GeV ", + energyCount / 1_GeV); + } + + template <typename TOutput> + inline void EnergyLossWriter<TOutput>::write(Point const& point, Code const, + HEPEnergyType const dE) { + GrammageType grammage = showerAxis_.getProjectedX(point); + int const maxBin = int(profile_.size() - 1); + int bin = grammage / dX_; + if (bin < 0) bin = 0; + if (bin > maxBin) bin = maxBin; + + CORSIKA_LOGGER_TRACE(TOutput::getLogger(), "add local energy loss bin={} dE={} GeV ", + bin, dE / 1_GeV); + + profile_[bin][static_cast<int>(dEdX_output::ProfileIndex::Total)] += dE; + } + + template <typename TOutput> + inline void EnergyLossWriter<TOutput>::write(GrammageType const Xstart, + GrammageType const Xend, Code const, + HEPEnergyType const dE) { + double const bstart = Xstart / dX_; + double const bend = Xend / dX_; + + if (abs(bstart - floor(bstart + 0.5)) > 1e-2 || + abs(bend - floor(bend + 0.5)) > 1e-2 || abs(bend - bstart - 1) > 1e-2) { + CORSIKA_LOGGER_ERROR( + TOutput::getLogger(), + "CascadeEquation (CONEX) and Corsika8 dX grammage binning are not the same! " + "Xstart={} Xend={} dX={} g/cm2", + Xstart / 1_g * square(1_cm), Xend / 1_g * square(1_cm), + dX_ / 1_g * square(1_cm)); + throw std::runtime_error( + "CONEX and Corsika8 dX grammage binning are not the same!"); + } + + size_t const bin = size_t((bend + bstart) / 2); + CORSIKA_LOGGER_TRACE(TOutput::getLogger(), + "add binned energy loss {} {} bin={} dE={} GeV ", bstart, bend, + bin, dE / 1_GeV); + if (bin >= profile_.size()) { + CORSIKA_LOGGER_WARN(TOutput::getLogger(), + "Grammage bin {} outside of profile {}. skipping.", bin, + profile_.size()); + return; + } + profile_[bin][static_cast<int>(dEdX_output::ProfileIndex::Total)] += dE; + } + + template <typename TOutput> + inline HEPEnergyType EnergyLossWriter<TOutput>::getEnergyLost() const { + HEPEnergyType tot = HEPEnergyType::zero(); + for (dEdX_output::Profile const& row : profile_) + tot += row.at(static_cast<int>(dEdX_output::ProfileIndex::Total)); + return tot; + } + + template <typename TOutput> + inline YAML::Node EnergyLossWriter<TOutput>::getConfig() const { + + YAML::Node node; + + node["type"] = "EnergyLoss"; + node["units"]["energy"] = "GeV"; + node["units"]["grammage"] = "g/cm^2"; + node["bin-size"] = dX_ / (1_g / square(1_cm)); + node["nbins"] = nBins_; + node["grammage_threshold"] = dX_threshold_ / (1_g / square(1_cm)); + + return node; + } + + template <typename TOutput> + inline YAML::Node EnergyLossWriter<TOutput>::getSummary() const { + + // determined Xmax and dEdXmax from quadratic interpolation + double maximum = 0; + size_t iMaximum = 0; + for (size_t i = 0; i < profile_.size() - 3; ++i) { + double value = + (profile_[i + 0].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) + + profile_[i + 1].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) + + profile_[i + 2].at(static_cast<int>(dEdX_output::ProfileIndex::Total))) / + 1_GeV; + if (value > maximum) { + maximum = value; + iMaximum = i; + } + } + + double const dX = dX_ / 1_g * square(1_cm); + + auto [Xmax, dEdXmax] = FindXmax::interpolateProfile( + dX * (0.5 + iMaximum), dX * (1.5 + iMaximum), dX * (2.5 + iMaximum), + profile_[iMaximum + 0].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) / + 1_GeV, + profile_[iMaximum + 1].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) / + 1_GeV, + profile_[iMaximum + 2].at(static_cast<int>(dEdX_output::ProfileIndex::Total)) / + 1_GeV); + + YAML::Node summary; + summary["sum_dEdX"] = getEnergyLost() / 1_GeV; + summary["Xmax"] = Xmax; + summary["dEdXmax"] = dEdXmax; + return summary; + } + +} // namespace corsika diff --git a/corsika/detail/modules/writers/EnergyLossWriterParquet.inl b/corsika/detail/modules/writers/EnergyLossWriterParquet.inl new file mode 100644 index 0000000000000000000000000000000000000000..55b301f509bc094a878e529e8ebb58360dc6568a --- /dev/null +++ b/corsika/detail/modules/writers/EnergyLossWriterParquet.inl @@ -0,0 +1,73 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/framework/utility/FindXmax.hpp> + +#include <corsika/media/ShowerAxis.hpp> + +#include <exception> + +namespace corsika { + + template <size_t NColumns> + inline EnergyLossWriterParquet<NColumns>::EnergyLossWriterParquet( + std::array<const char*, NColumns> const& colNames) + : columns_(colNames) {} + + template <size_t NColumns> + inline void EnergyLossWriterParquet<NColumns>::startOfLibrary( + boost::filesystem::path const& directory) { + + // setup the streamer + output_.initStreamer((directory / "dEdX.parquet").string()); + + // enable compression with the default level + // output_.enableCompression(); + + // build the schema + output_.addField("X", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + for (auto const& col : columns_) { + output_.addField(col, parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + } + + // and build the streamer + output_.buildStreamer(); + } + + template <size_t NColumns> + inline void EnergyLossWriterParquet<NColumns>::write( + unsigned int const showerId, GrammageType const grammage, + std::array<HEPEnergyType, NColumns> const& data) { + + // and write the data into the column + *(output_.getWriter()) << showerId + << static_cast<float>(grammage / 1_g * square(1_cm)); + for (HEPEnergyType const dedx : data) { + *(output_.getWriter()) << static_cast<float>(dedx / 1_GeV); + } + *(output_.getWriter()) << parquet::EndRow; + } + + template <size_t NColumns> + inline void EnergyLossWriterParquet<NColumns>::startOfShower(unsigned int const) {} + + template <size_t NColumns> + inline void EnergyLossWriterParquet<NColumns>::endOfShower(unsigned int const) {} + + template <size_t NColumns> + inline void EnergyLossWriterParquet<NColumns>::endOfLibrary() { + output_.closeStreamer(); + } + +} // namespace corsika diff --git a/corsika/detail/modules/writers/LongitudinalProfileWriterParquet.inl b/corsika/detail/modules/writers/LongitudinalProfileWriterParquet.inl new file mode 100644 index 0000000000000000000000000000000000000000..2e5801fee1882cc0ee44c522ae63fc1b205a5f4b --- /dev/null +++ b/corsika/detail/modules/writers/LongitudinalProfileWriterParquet.inl @@ -0,0 +1,75 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/framework/utility/FindXmax.hpp> + +#include <corsika/media/ShowerAxis.hpp> + +#include <string> +#include <exception> + +namespace corsika { + + template <size_t NColumns> + inline LongitudinalProfileWriterParquet<NColumns>::LongitudinalProfileWriterParquet( + std::array<const char*, NColumns> const& columns) + : columns_(columns) {} + + template <size_t NColumns> + inline void LongitudinalProfileWriterParquet<NColumns>::startOfLibrary( + boost::filesystem::path const& directory) { + // setup the streamer + output_.initStreamer((directory / "profile.parquet").string()); + + // enable compression with the default level + // output_.enableCompression(); + + // build the schema + output_.addField("X", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + for (auto const& col : columns_) { + output_.addField(col, parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + } + + // and build the streamer + output_.buildStreamer(); + } + + template <size_t NColumns> + inline void LongitudinalProfileWriterParquet<NColumns>::startOfShower( + unsigned int const) {} + + template <size_t NColumns> + inline void LongitudinalProfileWriterParquet<NColumns>::endOfShower( + unsigned int const) {} + + template <size_t NColumns> + inline void LongitudinalProfileWriterParquet<NColumns>::endOfLibrary() { + output_.closeStreamer(); + } + + template <size_t NColumns> + inline void LongitudinalProfileWriterParquet<NColumns>::write( + unsigned int const showerId, GrammageType const grammage, + std::array<double, NColumns> const& data) { + + // and write the data into the column + *(output_.getWriter()) << showerId + << static_cast<float>(grammage / 1_g * square(1_cm)); + for (double const weight : data) { + *(output_.getWriter()) << static_cast<float>(weight); + } + *(output_.getWriter()) << parquet::EndRow; + } + +} // namespace corsika diff --git a/corsika/detail/modules/writers/LongitudinalWriter.inl b/corsika/detail/modules/writers/LongitudinalWriter.inl new file mode 100644 index 0000000000000000000000000000000000000000..1095fcdb021d6d78b2aacb51bcda519ee820dfd6 --- /dev/null +++ b/corsika/detail/modules/writers/LongitudinalWriter.inl @@ -0,0 +1,168 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/framework/utility/FindXmax.hpp> + +#include <corsika/media/ShowerAxis.hpp> + +#include <exception> + +namespace corsika { + + template <typename TOutput> + inline LongitudinalWriter<TOutput>::LongitudinalWriter(ShowerAxis const& axis, + GrammageType dX, + size_t const nBins) + : TOutput(number_profile::ProfileIndexNames) + , showerAxis_(axis) + , dX_(dX) + , nBins_(nBins) {} + + template <typename TOutput> + inline void LongitudinalWriter<TOutput>::startOfLibrary( + boost::filesystem::path const& directory) { + TOutput::startOfLibrary(directory); + } + + template <typename TOutput> + inline void LongitudinalWriter<TOutput>::startOfShower(unsigned int const showerId) { + TOutput::startOfShower(showerId); + // reset profile + profile_.clear(); + profile_.resize(nBins_); + } + + template <typename TOutput> + inline void LongitudinalWriter<TOutput>::endOfLibrary() { + TOutput::endOfLibrary(); + } + + template <typename TOutput> + inline void LongitudinalWriter<TOutput>::endOfShower(unsigned int const showerId) { + + int iRow{0}; + for (number_profile::ProfileData const& row : profile_) { + // here: write to underlying writer (e.g. parquet) + TOutput::write(showerId, iRow * dX_, row); + iRow++; + } + TOutput::endOfShower(showerId); + } + + template <typename TOutput> + template <typename TTrack> + inline void LongitudinalWriter<TOutput>::write(TTrack const& track, Code const pid, + double const weight) { + GrammageType const grammageStart = showerAxis_.getProjectedX(track.getPosition(0)); + GrammageType const grammageEnd = showerAxis_.getProjectedX(track.getPosition(1)); + + // Note: particle may go also "upward", thus, grammageEnd<grammageStart + int const binStart = std::ceil(grammageStart / dX_); + int const binEnd = std::floor(grammageEnd / dX_); + + CORSIKA_LOGGER_TRACE(TOutput::getLogger(), + "grammageStart={} End={} binStart={}, end={}", + grammageStart / 1_g * square(1_cm), + grammageEnd / 1_g * square(1_cm), binStart, binEnd); + + for (int bin = binStart; bin <= binEnd; ++bin) { + if (pid == Code::Photon) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Photon)] += + weight; + } else if (pid == Code::Positron) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Positron)] += + weight; + } else if (pid == Code::Electron) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Electron)] += + weight; + } else if (pid == Code::MuPlus) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuPlus)] += + weight; + } else if (pid == Code::MuMinus) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuMinus)] += + weight; + } else if (is_hadron(pid)) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Hadron)] += + weight; + } + if (is_charged(pid)) { + profile_[bin][static_cast<int>(number_profile::ProfileIndex::Charged)] += weight; + } + } + } + + template <typename TOutput> + inline void LongitudinalWriter<TOutput>::write(GrammageType const Xstart, + GrammageType const Xend, Code const pid, + double const weight) { + double const bstart = Xstart / dX_; + double const bend = Xend / dX_; + + if (abs(bstart - floor(bstart + 0.5)) > 1e-2 || + abs(bend - floor(bend + 0.5)) > 1e-2 || abs(bend - bstart - 1) > 1e-2) { + CORSIKA_LOGGER_ERROR(TOutput::getLogger(), + "CONEX and Corsika8 dX grammage binning are not the same! " + "Xstart={} Xend={} dX={}", + Xstart / 1_g * square(1_cm), Xend / 1_g * square(1_cm), + dX_ / 1_g * square(1_cm)); + throw std::runtime_error( + "CONEX and Corsika8 dX grammage binning are not the same!"); + } + + size_t const bin = size_t((bend + bstart) / 2); + if (bin >= profile_.size()) { + CORSIKA_LOGGER_WARN(TOutput::getLogger(), + "Grammage bin {} outside of profile {}. skipping.", bin, + profile_.size()); + return; + } + + if (pid == Code::Photon) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Photon)] += weight; + } else if (pid == Code::Positron) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Positron)] += + weight; + } else if (pid == Code::Electron) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Electron)] += + weight; + } else if (pid == Code::MuPlus) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuPlus)] += weight; + } else if (pid == Code::MuMinus) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::MuMinus)] += weight; + } else if (is_hadron(pid)) { + profile_.at(bin)[static_cast<int>(number_profile::ProfileIndex::Hadron)] += weight; + } + if (is_charged(pid)) { + profile_[bin][static_cast<int>(number_profile::ProfileIndex::Charged)] += weight; + } + } + + template <typename TOutput> + inline YAML::Node LongitudinalWriter<TOutput>::getConfig() const { + + YAML::Node node; + + node["type"] = "LongitudinalProfile"; + node["units"]["grammage"] = "g/cm^2"; + node["bin-size"] = dX_ / (1_g / square(1_cm)); + node["nbins"] = nBins_; + + return node; + } + + template <typename TOutput> + inline YAML::Node LongitudinalWriter<TOutput>::getSummary() const { + YAML::Node summary; + return summary; + } + +} // namespace corsika diff --git a/corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl b/corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl deleted file mode 100644 index 5eef28a41c8d3e487501157fd79a64d7e37c934e..0000000000000000000000000000000000000000 --- a/corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl +++ /dev/null @@ -1,57 +0,0 @@ -/* - * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu - * - * This software is distributed under the terms of the GNU General Public - * Licence version 3 (GPL Version 3). See file LICENSE for a full version of - * the license. - */ - -#pragma once - -namespace corsika { - - inline ObservationPlaneWriterParquet::ObservationPlaneWriterParquet() - : output_() {} - - inline void ObservationPlaneWriterParquet::startOfLibrary( - boost::filesystem::path const& directory) { - - // setup the streamer - output_.initStreamer((directory / "particles.parquet").string()); - - // enable compression with the default level - output_.enableCompression(); - - // build the schema - output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32, - parquet::ConvertedType::INT_32); - output_.addField("energy", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, - parquet::ConvertedType::NONE); - output_.addField("x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, - parquet::ConvertedType::NONE); - output_.addField("y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, - parquet::ConvertedType::NONE); - output_.addField("t", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, - parquet::ConvertedType::NONE); - - // and build the streamer - output_.buildStreamer(); - } - - inline void ObservationPlaneWriterParquet::endOfShower() { ++shower_; } - - inline void ObservationPlaneWriterParquet::endOfLibrary() { output_.closeStreamer(); } - - inline void ObservationPlaneWriterParquet::write(Code const& pid, - HEPEnergyType const& energy, - LengthType const& x, - LengthType const& y, - TimeType const& t) { - // write the next row - we must write `shower_` first. - *(output_.getWriter()) << shower_ << static_cast<int>(get_PDG(pid)) - << static_cast<float>(energy / 1_GeV) - << static_cast<float>(x / 1_m) << static_cast<float>(y / 1_m) - << static_cast<float>(t / 1_ns) << parquet::EndRow; - } - -} // namespace corsika diff --git a/corsika/detail/modules/writers/ParticleWriterParquet.inl b/corsika/detail/modules/writers/ParticleWriterParquet.inl new file mode 100644 index 0000000000000000000000000000000000000000..4c6a49adadc103d38e2dbf5ecd278fcf262b189f --- /dev/null +++ b/corsika/detail/modules/writers/ParticleWriterParquet.inl @@ -0,0 +1,99 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/framework/core/Logging.hpp> + +namespace corsika { + + inline ParticleWriterParquet::ParticleWriterParquet() + : output_() + , showerId_(0) + , totalEnergy_(0_eV) {} + + inline void ParticleWriterParquet::startOfLibrary( + boost::filesystem::path const& directory) { + + // setup the streamer + output_.initStreamer((directory / "particles.parquet").string()); + + // enable compression with the default level + // output_.enableCompression(); + + // build the schema + output_.addField("pdg", parquet::Repetition::REQUIRED, parquet::Type::INT32, + parquet::ConvertedType::INT_32); + output_.addField("energy", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + output_.addField("x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + output_.addField("y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + output_.addField("z", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + output_.addField("weight", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); + + // and build the streamer + output_.buildStreamer(); + + showerId_ = 0; + } + + inline void ParticleWriterParquet::startOfShower(unsigned int const showerId) { + showerId_ = showerId; + totalEnergy_ = 0_eV; + countHadrons_ = 0; + countOthers_ = 0; + countEM_ = 0; + countMuons_ = 0; + } + + inline void ParticleWriterParquet::endOfShower(unsigned int const) {} + + inline void ParticleWriterParquet::endOfLibrary() { output_.closeStreamer(); } + + inline void ParticleWriterParquet::write(Code const& pid, HEPEnergyType const& energy, + LengthType const& x, LengthType const& y, + LengthType const& z, double const weight) { + + // write the next row - we must write `shower_` first. + *(output_.getWriter()) << showerId_ << static_cast<int>(get_PDG(pid)) + << static_cast<float>(energy / 1_GeV) + << static_cast<float>(x / 1_m) << static_cast<float>(y / 1_m) + << static_cast<float>(z / 1_m) << static_cast<float>(weight) + << parquet::EndRow; + + totalEnergy_ += energy; + + if (is_hadron(pid)) { + ++countHadrons_; + } else if (is_muon(pid)) { + ++countMuons_; + } else if (is_em(pid)) { + ++countEM_; + } else { + ++countOthers_; + } + } + + /** + * Return collected library-level summary for output. + */ + YAML::Node ParticleWriterParquet::getSummary() const { + YAML::Node summary; + summary["Eground"] = totalEnergy_ / 1_GeV; + summary["hadrons"] = countHadrons_; + summary["muons"] = countMuons_; + summary["em"] = countEM_; + summary["others"] = countOthers_; + return summary; + } + +} // namespace corsika diff --git a/corsika/detail/modules/writers/TrackWriterParquet.inl b/corsika/detail/modules/writers/TrackWriterParquet.inl index 90e0ba2c720e0ec079e853c9c9fac62f4a2ddba7..88dfa8e03c7160802638d1b8afc9ad9bf1887ece 100644 --- a/corsika/detail/modules/writers/TrackWriterParquet.inl +++ b/corsika/detail/modules/writers/TrackWriterParquet.inl @@ -11,7 +11,8 @@ namespace corsika { inline TrackWriterParquet::TrackWriterParquet() - : output_() {} + : output_() + , showerId_(0) {} inline void TrackWriterParquet::startOfLibrary( boost::filesystem::path const& directory) { @@ -24,6 +25,8 @@ namespace corsika { parquet::ConvertedType::INT_32); output_.addField("energy", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, parquet::ConvertedType::NONE); + output_.addField("weight", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, + parquet::ConvertedType::NONE); output_.addField("start_x", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, parquet::ConvertedType::NONE); output_.addField("start_y", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, @@ -43,24 +46,31 @@ namespace corsika { // and build the streamer output_.buildStreamer(); + showerId_ = 0; + } + + inline void TrackWriterParquet::startOfShower(unsigned int const showerId) { + showerId_ = showerId; } - inline void TrackWriterParquet::endOfShower() { ++shower_; } + inline void TrackWriterParquet::endOfShower(unsigned int const) {} inline void TrackWriterParquet::endOfLibrary() { output_.closeStreamer(); } - inline void TrackWriterParquet::write(Code const& pid, HEPEnergyType const& energy, + inline void TrackWriterParquet::write(Code const pid, HEPEnergyType const energy, + double const weight, QuantityVector<length_d> const& start, - TimeType const& t_start, + TimeType const t_start, QuantityVector<length_d> const& end, - TimeType const& t_end) { + TimeType const t_end) { // write the next row - we must write `shower_` first. // clang-format off *(output_.getWriter()) - << shower_ + << showerId_ << static_cast<int>(get_PDG(pid)) << static_cast<float>(energy / 1_GeV) + << static_cast<float>(weight) << static_cast<float>(start[0] / 1_m) << static_cast<float>(start[1] / 1_m) << static_cast<float>(start[2] / 1_m) diff --git a/corsika/detail/output/BaseOutput.inl b/corsika/detail/output/BaseOutput.inl deleted file mode 100644 index 24ec861ac9c2b3083ea834261f7fa0a57fd8559e..0000000000000000000000000000000000000000 --- a/corsika/detail/output/BaseOutput.inl +++ /dev/null @@ -1,15 +0,0 @@ -/* - * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu - * - * This software is distributed under the terms of the GNU General Public - * Licence version 3 (GPL Version 3). See file LICENSE for a full version of - * the license. - */ -#pragma once - -namespace corsika { - - inline BaseOutput::BaseOutput() - : shower_(0) {} - -} // namespace corsika diff --git a/corsika/detail/output/OutputManager.inl b/corsika/detail/output/OutputManager.inl index 1fa7e0267c66af797cdf851d45c16f2146527776..af813534ea6761f502ab2f9ac33a5a8fc7ef3433 100644 --- a/corsika/detail/output/OutputManager.inl +++ b/corsika/detail/output/OutputManager.inl @@ -22,23 +22,78 @@ namespace corsika { - inline void OutputManager::writeNode(YAML::Node const& node, - boost::filesystem::path const& path) const { + inline OutputManager::OutputManager( + std::string const& name, + boost::filesystem::path const& dir = boost::filesystem::current_path()) + : root_(dir / name) + , name_(name) + , count_(0) { - // construct a YAML emitter for this config file - YAML::Emitter out; + // check if this directory already exists + if (boost::filesystem::exists(root_)) { + CORSIKA_LOGGER_ERROR(logger_, + "Output directory '{}' already exists! Do not overwrite!.", + root_.string()); + throw std::runtime_error("Output directory already exists."); + } - // and write the node to the output - out << node; + // construct the directory for this library + boost::filesystem::create_directories(root_); + + CORSIKA_LOGGER_INFO(logger_, fmt::format("Output library: \"{}\"", root_.string())); + writeYAML(getConfig(), root_ / ("config.yaml")); + } - // open the output file - this is <output name>.yaml - boost::filesystem::ofstream file(path); + template <typename TOutput> + inline void OutputManager::add(std::string const& name, TOutput& output) { + + if (state_ != OutputState::NoInit) { + // if "add" is called after the ouptput has started, this is an ERROR. + CORSIKA_LOGGER_ERROR( + logger_, "Cannot add more outputs to OutputManager after output was started."); + throw std::runtime_error( + "Cannot add more outputs to OutputManager after output was started."); + } - // dump the YAML to the file - file << out.c_str() << std::endl; + // check if that name is already in the map + if (outputs_.count(name) > 0) { + CORSIKA_LOGGER_ERROR( + logger_, "'{}' is already registered. All outputs must have unique names.", + name); + throw std::runtime_error("Output already exists. Do not overwrite!"); + } + + // if we get here, the name is not already in the map + // so we create the output and register it into the map + outputs_.insert(std::make_pair(name, std::ref(output))); + + // create the directory for this process. + boost::filesystem::create_directory(root_ / name); } - inline void OutputManager::writeTopLevelConfig() const { + inline OutputManager::~OutputManager() { + + if (state_ == OutputState::ShowerInProgress) { + // if this the destructor is called before the shower has been explicitly + // ended, print a warning and end the shower before continuing. + CORSIKA_LOGGER_WARN(logger_, + "OutputManager was destroyed before endOfShower() called." + " The last shower in this libray may be incomplete."); + endOfShower(); + } + + // write the top level summary file (summary.yaml) + writeSummary(); + + // if we are being destructed but EndOfLibrary() has not been called, + // make sure that we gracefully close all the outputs. This is a supported + // method of operation so we don't issue a warning here + if (state_ == OutputState::LibraryReady) { endOfLibrary(); } + } + + inline int OutputManager::getEventId() const { return count_; } + + inline YAML::Node OutputManager::getConfig() const { YAML::Node config; @@ -47,16 +102,15 @@ namespace corsika { config["creator"] = "CORSIKA8"; // a tag to identify C8 libraries config["version"] = "8.0.0-prealpha"; // the current version - // write the node to a file - writeNode(config, root_ / ("config.yaml")); + return config; } - inline void OutputManager::writeTopLevelSummary() const { + inline YAML::Node OutputManager::getSummary() const { - YAML::Node config; + YAML::Node summary; // the total number of showers contained in the library - config["showers"] = count_; + summary["showers"] = count_; // this next section handles writing some time and duration information @@ -81,87 +135,17 @@ namespace corsika { auto runtime{end_time - start_time}; // add the time and duration info - config["start time"] = timeToString(start_time); - config["end time"] = timeToString(end_time); - config["runtime"] = fmt::format("{:%H:%M:%S}", runtime); - - // write the node to a file - writeNode(config, root_ / ("summary.yaml")); - } - - inline void OutputManager::initOutput(std::string const& name) const { - // construct the path to this directory - auto const path{root_ / name}; - - // create the directory for this process. - boost::filesystem::create_directory(path); - - // get the config for this output - auto config = outputs_.at(name).get().getConfig(); - - // and assign the name for this output - config["name"] = name; - - // write the config for this output to the file - writeNode(config, path / "config.yaml"); - } - - inline OutputManager::OutputManager( - std::string const& name, - boost::filesystem::path const& dir = boost::filesystem::current_path()) - : name_(name) - , root_(dir / name) { - - // check if this directory already exists - if (boost::filesystem::exists(root_)) { - logger->error("Output directory '{}' already exists! Do not overwrite!.", - root_.string()); - throw std::runtime_error("Output directory already exists."); - } + summary["start time"] = timeToString(start_time); + summary["end time"] = timeToString(end_time); + summary["runtime"] = fmt::format("{:%H:%M:%S}", runtime); - // construct the directory for this library - boost::filesystem::create_directories(root_); - - // write the top level config file - writeTopLevelConfig(); + return summary; } - inline OutputManager::~OutputManager() { + inline void OutputManager::writeSummary() const { - if (state_ == OutputState::ShowerInProgress) { - // if this the destructor is called before the shower has been explicitly - // ended, print a warning and end the shower before continuing. - logger->warn( - "OutputManager was destroyed before endOfShower() called." - " The last shower in this libray may be incomplete."); - endOfShower(); - } - - // write the top level summary file (summary.yaml) - writeTopLevelSummary(); - - // if we are being destructed but EndOfLibrary() has not been called, - // make sure that we gracefully close all the outputs. This is a supported - // method of operation so we don't issue a warning here - if (state_ == OutputState::LibraryReady) { endOfLibrary(); } - } - - template <typename TOutput> - inline void OutputManager::add(std::string const& name, TOutput& output) { - - // check if that name is already in the map - if (outputs_.count(name) > 0) { - logger->error("'{}' is already registered. All outputs must have unique names.", - name); - throw std::runtime_error("Output already exists. Do not overwrite!"); - } - - // if we get here, the name is not already in the map - // so we create the output and register it into the map - outputs_.insert(std::make_pair(name, std::ref(output))); - - // and initialize this output - initOutput(name); + // write the node to a file + writeYAML(getSummary(), root_ / ("summary.yaml")); } inline void OutputManager::startOfLibrary() { @@ -176,15 +160,22 @@ namespace corsika { // we now forward this signal to all of our outputs for (auto& [name, output] : outputs_) { - // construct the path to this output subdirectory - auto const path{root_ / name}; - // and start the library - output.get().startOfLibrary(path); + output.get().startOfLibrary(root_ / name); + + // get the config from this output + auto config = output.get().getConfig(); + + // add the name keyword + config["name"] = name; + + // write the output configuration to config.yaml in the output directory + writeYAML(config, root_ / name / ("config.yaml")); } // we have now started running state_ = OutputState::LibraryReady; + count_ = 0; // event counter } inline void OutputManager::startOfShower() { @@ -194,15 +185,7 @@ namespace corsika { if (state_ == OutputState::NoInit) { startOfLibrary(); } // now start the event for all the outputs - for (auto& [name, output] : outputs_) { - { - [[maybe_unused]] auto const& dummy_name = name; - } - output.get().startOfShower(); - } - - // increment our shower count - ++count_; + for (auto& [name, output] : outputs_) { output.get().startOfShower(count_); } // and transition to the in progress state state_ = OutputState::ShowerInProgress; @@ -210,15 +193,13 @@ namespace corsika { inline void OutputManager::endOfShower() { - for (auto& [name, output] : outputs_) { - { - [[maybe_unused]] auto const& dummy_name = name; - } - output.get().endOfShower(); - } + for (auto& [name, output] : outputs_) { output.get().endOfShower(count_); } // switch back to the initialized state state_ = OutputState::LibraryReady; + + // increment our shower count + ++count_; } inline void OutputManager::endOfLibrary() { @@ -230,12 +211,11 @@ namespace corsika { // write the summary for each output and forward the endOfLibrary call() for (auto& [name, output] : outputs_) { - - // we get the summary for each output as a YAML node - auto summary{outputs_.at(name).get().getSummary()}; - - // write the summary for this output to the file - writeNode(summary, root_ / name / "summary.yaml"); + // save eventual YAML summary + YAML::Node const summary = output.get().getSummary(); + if (!summary.IsNull()) { + writeYAML(output.get().getSummary(), root_ / name / ("summary.yaml")); + } // and forward the end of library call output.get().endOfLibrary(); @@ -243,6 +223,6 @@ namespace corsika { // and the library has finished state_ = OutputState::LibraryFinished; - } + } // namespace corsika } // namespace corsika diff --git a/corsika/detail/output/ParquetStreamer.inl b/corsika/detail/output/ParquetStreamer.inl index d0bb2f9e01a42fef34c00fd2d41c189e60f252e2..ddef90eedb8afb992d13f4dadcd994c0e1501c09 100644 --- a/corsika/detail/output/ParquetStreamer.inl +++ b/corsika/detail/output/ParquetStreamer.inl @@ -23,7 +23,7 @@ namespace corsika { // add run and event tags to the file addField("shower", parquet::Repetition::REQUIRED, parquet::Type::INT32, - parquet::ConvertedType::INT_32); + parquet::ConvertedType::UINT_32); } template <typename... TArgs> @@ -31,9 +31,9 @@ namespace corsika { fields_.push_back(parquet::schema::PrimitiveNode::Make(args...)); } - inline void ParquetStreamer::enableCompression(int const /*level*/) { - // builder_.compression(parquet::Compression::ZSTD); - // builder_.compression_level(level); + inline void ParquetStreamer::enableCompression(int const level) { + builder_.compression(parquet::Compression::LZ4); + builder_.compression_level(level); } inline void ParquetStreamer::buildStreamer() { diff --git a/corsika/detail/output/YAMLStreamer.inl b/corsika/detail/output/YAMLStreamer.inl new file mode 100644 index 0000000000000000000000000000000000000000..bdca88ac26a168c26367edebad941d06bdf7297f --- /dev/null +++ b/corsika/detail/output/YAMLStreamer.inl @@ -0,0 +1,29 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +namespace corsika { + + inline void YAMLStreamer::writeYAML(YAML::Node const& node, + boost::filesystem::path const& path) const { + + // construct a YAML emitter for this config file + YAML::Emitter out; + + // and write the node to the output + out << node; + + // open the output file - this is <output name>.yaml + boost::filesystem::ofstream file(path); + + // dump the YAML to the file + file << out.c_str() << std::endl; + } + +} // namespace corsika diff --git a/corsika/detail/setup/SetupStack.hpp b/corsika/detail/setup/SetupStack.hpp index 7a1327b2e54f6dfc81824e001d2a290a4f612e8b..2b40ccbdb0eafcc2c48ead0dc437c1314d6e2286 100644 --- a/corsika/detail/setup/SetupStack.hpp +++ b/corsika/detail/setup/SetupStack.hpp @@ -22,7 +22,8 @@ namespace corsika { namespace setup::detail { // ------------------------------------------ - // add geometry node tracking data to stack: + // add geometry node data to stack. This is fundamentally needed + // for robust tracking through multiple volumes. // the GeometryNode stack needs to know the type of geometry-nodes from the // environment: @@ -42,16 +43,36 @@ namespace corsika { DefaultSecondaryProducer>; // ------------------------------------------ - // Add [optional] history data to stack, too: + // add weight data to stack. This is fundamentally needed + // for thinning. + + // the "pure" weight stack (interface) + template <typename TStackIter> + using SetupWeightDataInterface = + typename weights::MakeWeightDataInterface<TStackIter>::type; + + // combine geometry-node-vector data stack with weight information for tracking + template <typename TStackIter> + using StackWithWeightInterface = + CombinedParticleInterface<StackWithGeometry::pi_type, SetupWeightDataInterface, + TStackIter>; + + // the combined stack data: particle + geometry + weight + using StackWithWeight = + CombinedStack<typename StackWithGeometry::stack_data_type, weights::WeightData, + StackWithWeightInterface, DefaultSecondaryProducer>; + + // ------------------------------------------ + // Add [OPTIONAL] history data to stack, too. + // This keeps the entire lineage of particles in memory. - // combine dummy stack with geometry information for tracking template <typename TStackIter> using StackWithHistoryInterface = - CombinedParticleInterface<StackWithGeometry::pi_type, + CombinedParticleInterface<StackWithWeight::pi_type, history::HistoryEventDataInterface, TStackIter>; using StackWithHistory = - CombinedStack<typename StackWithGeometry::stack_data_type, + CombinedStack<typename StackWithWeight::stack_data_type, history::HistoryEventData, StackWithHistoryInterface, history::HistorySecondaryProducer>; diff --git a/corsika/framework/core/Cascade.hpp b/corsika/framework/core/Cascade.hpp index e9b6fbaf2599fc3b41bcf1ad0976ef2e426e7c85..398fdcb33541bcd5322e1dcfe415d4999439582e 100644 --- a/corsika/framework/core/Cascade.hpp +++ b/corsika/framework/core/Cascade.hpp @@ -73,21 +73,7 @@ namespace corsika { ~Cascade() = default; Cascade& operator=(Cascade const&) = default; Cascade(Environment<medium_interface_type> const& env, TTracking& tr, - TProcessList& pl, TOutput& out, TStack& stack) - : environment_(env) - , tracking_(tr) - , sequence_(pl) - , output_(out) - , stack_(stack) { - CORSIKA_LOG_INFO(c8_ascii_); - CORSIKA_LOG_INFO("This is CORSIKA {}.{}.{}.{}", CORSIKA_RELEASE_NUMBER, - CORSIKA_MAJOR_NUMBER, CORSIKA_MINOR_NUMBER, CORSIKA_PATCH_NUMBER); - CORSIKA_LOG_INFO("Tracking algorithm: {} (version {})", TTracking::getName(), - TTracking::getVersion()); - if constexpr (stack_view_type::has_event) { - CORSIKA_LOG_INFO("Stack - with full cascade HISTORY."); - } - } + TProcessList& pl, TOutput& out, TStack& stack); //! @} /** @@ -135,6 +121,7 @@ namespace corsika { TOutput& output_; TStack& stack_; default_prng_type& rng_ = RNGManager<>::getInstance().getRandomStream("cascade"); + bool forceInteraction_; unsigned int count_ = 0; // but this here temporarily. Should go into dedicated file later: diff --git a/corsika/framework/core/ParticleProperties.hpp b/corsika/framework/core/ParticleProperties.hpp index 1b344954194c056005874eff5a28d0f15a1468ab..dd03aa106c309b6abe3e7f58b69922ebbb8ae875 100644 --- a/corsika/framework/core/ParticleProperties.hpp +++ b/corsika/framework/core/ParticleProperties.hpp @@ -92,7 +92,7 @@ namespace corsika { int16_t constexpr get_charge_number(Code const); //!< electric charge in units of e ElectricChargeType constexpr get_charge(Code const); //!< electric charge HEPMassType constexpr get_mass(Code const); //!< mass - HEPEnergyType constexpr calculate_kinetic_energy_threshold( + HEPEnergyType constexpr get_kinetic_energy_threshold( Code const); //!< get kinetic energy threshold below which the particle is //!< discarded, by default set to zero void constexpr set_kinetic_energy_threshold( @@ -109,7 +109,7 @@ namespace corsika { //! Particle code according to PDG, "Monte Carlo Particle Numbering Scheme" PDGCode constexpr get_PDG(Code const); - + PDGCode constexpr get_PDG(unsigned int const A, unsigned int const Z); std::string_view constexpr get_name(Code const); //!< name of the particle as string TimeType constexpr get_lifetime(Code const); //!< lifetime @@ -117,6 +117,7 @@ namespace corsika { bool constexpr is_em(Code const); //!< true if particle is electron, positron or photon bool constexpr is_muon(Code const); //!< true if particle is mu+ or mu- bool constexpr is_neutrino(Code const); //!< true if particle is (anti-) neutrino + bool constexpr is_charged(Code const); //!< true if particle is charged /** * @brief Creates the Code for a nucleus of type 10LZZZAAAI. diff --git a/corsika/framework/utility/FindXmax.hpp b/corsika/framework/utility/FindXmax.hpp new file mode 100644 index 0000000000000000000000000000000000000000..5b3e1e6b946408b1c81d5c99785fda6d265fbea7 --- /dev/null +++ b/corsika/framework/utility/FindXmax.hpp @@ -0,0 +1,120 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/framework/core/Logging.hpp> + +#include <tuple> + +namespace corsika { + + /** + * Interpolates profiles around maximum. + * + * The maximum of profiles can be robustly estimated from the three maximal points by + * quadratic interpolation. This code is copied from CONEX and is awaiting full adaption + * into CORSIKA8. This is a temporary solution (just copied). + */ + + class FindXmax { + + static bool invert3by3(double A[3][3]) { + const double kSmall = 1e-80; + + double determinant = A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1]) - + A[0][1] * (A[1][0] * A[2][2] - A[1][2] * A[2][0]) + + A[0][2] * (A[1][0] * A[2][1] - A[1][1] * A[2][0]); + + double absDet = fabs(determinant); + + if (absDet < kSmall) { + CORSIKA_LOG_WARN("invert3by3: Error-matrix singular (absDet={})", absDet); + return false; + } + + double B[3][3]; + + B[0][0] = A[1][1] * A[2][2] - A[1][2] * A[2][1]; + B[1][0] = A[1][2] * A[2][0] - A[2][2] * A[1][0]; + B[2][0] = A[1][0] * A[2][1] - A[1][1] * A[2][0]; + + B[0][1] = A[0][2] * A[2][1] - A[2][2] * A[0][1]; + B[1][1] = A[0][0] * A[2][2] - A[2][0] * A[0][2]; + B[2][1] = A[0][1] * A[2][0] - A[0][0] * A[2][1]; + + B[0][2] = A[0][1] * A[1][2] - A[1][1] * A[0][2]; + B[1][2] = A[0][2] * A[1][0] - A[1][2] * A[0][0]; + B[2][2] = A[0][0] * A[1][1] - A[0][1] * A[1][0]; + + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) A[i][j] = B[i][j] / determinant; + + return true; + } + + /**************************************************** + * + * solves linear system + * + * / y[0] \ / A[0][0] A[0][1] A[0][2] \ / x[0] \ + * | y[1] | = | A[1][0] A[1][1] A[1][2] | * | x[1] | + * \ y[2] / \ A[2][0] A[2][1] A[2][2] / \ x[2] / + * + * + * Input: y[3] and A[3][3] + * Output: returns true when succeded (i.e. A is not singular) + * A is overwritten with its inverse + * + * M.Unger 12/1/05 + * + ****************************************************/ + static bool solve3by3(double y[3], double A[3][3], double x[3]) { + if (invert3by3(A)) { + + for (int i = 0; i < 3; i++) + x[i] = A[i][0] * y[0] + A[i][1] * y[1] + A[i][2] * y[2]; + + return true; + + } else + return false; + } + + public: + static std::tuple<double, double> interpolateProfile(double x1, double x2, double x3, + double y1, double y2, + double y3) { + + // quadratic "fit" around maximum to get dEdXmax and Xmax + double x[3] = {x1, x2, x3}; + double y[3] = {y1, y2, y3}; + + double A[3][3]; + A[0][0] = x[0] * x[0]; + A[0][1] = x[0]; + A[0][2] = 1.; + A[1][0] = x[1] * x[1]; + A[1][1] = x[1]; + A[1][2] = 1.; + A[2][0] = x[2] * x[2]; + A[2][1] = x[2]; + A[2][2] = 1.; + + double a[3]; + + solve3by3(y, A, a); + + if (a[0] < 0.) { + double const Xmax = -a[1] / (2. * a[0]); + return std::make_tuple(Xmax, a[0] * Xmax * Xmax + a[1] * Xmax + a[2]); + } + return std::make_tuple(0, 0); + } + }; +} // namespace corsika \ No newline at end of file diff --git a/corsika/modules/LongitudinalProfile.hpp b/corsika/modules/LongitudinalProfile.hpp index b5d199fb87949e5456e47b097eea65248d8d5be2..af9c89a768cf726fc9de604bd6c821fc5734011f 100644 --- a/corsika/modules/LongitudinalProfile.hpp +++ b/corsika/modules/LongitudinalProfile.hpp @@ -12,6 +12,8 @@ #include <corsika/framework/process/ContinuousProcess.hpp> #include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/modules/writers/LongitudinalProfileWriterParquet.hpp> + #include <array> #include <fstream> #include <limits> @@ -33,11 +35,13 @@ namespace corsika { * boundaries. */ - class LongitudinalProfile : public ContinuousProcess<LongitudinalProfile> { + template <typename TOutput> + class LongitudinalProfile : public ContinuousProcess<LongitudinalProfile<TOutput>>, + public TOutput { public: - LongitudinalProfile(ShowerAxis const&, - GrammageType dX = 10_g / square(1_cm)); // profile binning); + template <typename... TArgs> + LongitudinalProfile(TArgs&&... args); template <typename TParticle, typename TTrack> ProcessReturn doContinuous( @@ -49,22 +53,7 @@ namespace corsika { return meter * std::numeric_limits<double>::infinity(); } - void save(std::string const&, int const width = 14, int const precision = 6); - - private: - GrammageType const dX_; - ShowerAxis const& shower_axis_; - using ProfileEntry = std::array<uint32_t, 7>; - enum ProfileIndex { - Photon = 0, - Positron = 1, - Electron = 2, - MuPlus = 3, - MuMinus = 4, - Hadron = 5, - Invisible = 6, - }; - std::vector<ProfileEntry> profiles_; // longitudinal profile + YAML::Node getConfig() const; }; } // namespace corsika diff --git a/corsika/modules/ObservationPlane.hpp b/corsika/modules/ObservationPlane.hpp index 6d71d8f7611e98b7b38044782538a06c3cb69a64..997e86782eb84ebf6dbbf7d5a632eb38580ed60b 100644 --- a/corsika/modules/ObservationPlane.hpp +++ b/corsika/modules/ObservationPlane.hpp @@ -11,32 +11,49 @@ #include <corsika/framework/core/PhysicalUnits.hpp> #include <corsika/framework/geometry/Plane.hpp> #include <corsika/framework/process/ContinuousProcess.hpp> -#include <corsika/modules/writers/ObservationPlaneWriterParquet.hpp> +#include <corsika/modules/writers/ParticleWriterParquet.hpp> +#include <corsika/modules/writers/WriterOff.hpp> namespace corsika { /** - @ingroup Modules - @{ - - The ObservationPlane writes PDG codes, energies, and distances of particles to the - central point of the plane into its output file. The particles are considered - "absorbed" afterwards. - - **Note/Limitation:** as discussed in - https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/397 - you cannot put two ObservationPlanes exactly on top of each - other. Even if one of them is "permeable". You have to put a - small gap in between the two plane in such a scenario, or develop - another more specialized output class. + * @ingroup Modules + * @{ + * + * The ObservationPlane writes PDG codes, energies, and distances of particles to the + * central point of the plane into its output file. The particles are considered + * "absorbed" afterwards. + * + * The default output format is parquet. + * + * **Note/Limitation:** as discussed in + * https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/397 + * you cannot put two ObservationPlanes exactly on top of each + * other. Even if one of them is "permeable". You have to put a + * small gap in between the two plane in such a scenario, or develop + * another more specialized output class. */ - template <typename TTracking, typename TOutputWriter = ObservationPlaneWriterParquet> - class ObservationPlane - : public ContinuousProcess<ObservationPlane<TTracking, TOutputWriter>>, - public TOutputWriter { + template <typename TTracking, typename TOutput = ParticleWriterParquet> + class ObservationPlane : public ContinuousProcess<ObservationPlane<TTracking, TOutput>>, + public TOutput { + + using TOutput::write; public: - ObservationPlane(Plane const&, DirectionVector const&, bool = true); + /** + * Construct a new Observation Plane object. + * + * @tparam TArgs + * @param plane The plane. + * @param x_dir The x-direction/axis. + * @param absorbing Flag to make the plane absorbing. + * @param outputArgs + */ + template <typename... TArgs> + ObservationPlane(Plane const& plane, DirectionVector const& x_dir, + bool const absorbing = true, TArgs&&... outputArgs); + + ~ObservationPlane() {} template <typename TParticle, typename TTrajectory> ProcessReturn doContinuous(TParticle& vParticle, TTrajectory& vTrajectory, @@ -45,18 +62,13 @@ namespace corsika { template <typename TParticle, typename TTrajectory> LengthType getMaxStepLength(TParticle const&, TTrajectory const& vTrajectory); - void showResults() const; - void reset(); - HEPEnergyType getEnergyGround() const { return energy_ground_; } YAML::Node getConfig() const; private: Plane const plane_; - bool const deleteOnHit_; - HEPEnergyType energy_ground_; - unsigned int count_ground_; DirectionVector const xAxis_; DirectionVector const yAxis_; + bool const deleteOnHit_; }; //! @} } // namespace corsika diff --git a/corsika/modules/ParticleCut.hpp b/corsika/modules/ParticleCut.hpp index 159186e1df41de1ed27bd1c8a8d6b3d34ea1a0a1..4cc08634ad7f3218f31c44659173d69876921203 100644 --- a/corsika/modules/ParticleCut.hpp +++ b/corsika/modules/ParticleCut.hpp @@ -15,82 +15,99 @@ #include <corsika/framework/process/SecondariesProcess.hpp> #include <corsika/framework/process/ContinuousProcess.hpp> +#include <corsika/modules/writers/WriterOff.hpp> + namespace corsika { /** - simple ParticleCut process. Goes through the secondaries of an interaction and - removes particles according to their kinetic energy. Particles with a time delay of - more than 10ms are removed as well. Invisible particles (neutrinos) can be removed if - selected. The threshold value is set to 0 by default but in principle can be configured - for each particle. Special constructors for cuts by the following groups are - implemented: (electrons,positrons), photons, hadrons and muons. - **/ - class ParticleCut : public SecondariesProcess<ParticleCut>, - public ContinuousProcess<ParticleCut> { + * ParticleCut process to kill particles. + * + * Goes through the secondaries of an interaction and + * removes particles according to their kinetic energy. Particles with a time delay of + * more than 10ms are removed as well. Invisible particles (neutrinos) can be removed if + * selected. The threshold value is set to 0 by default but in principle can be + * configured for each particle. Special constructors for cuts by the following groups + * are implemented: (electrons,positrons), photons, hadrons and muons. + */ + template <typename TOutput = WriterOff> + class ParticleCut : public SecondariesProcess<ParticleCut<TOutput>>, + public ContinuousProcess<ParticleCut<TOutput>>, + public TOutput { public: /** * particle cut with kinetic energy thresholds for electrons, photons, - * hadrons (including nuclei with energy per nucleon) and muons - * invisible particles (neutrinos) can be cut or not - **/ + * hadrons (including nuclei with energy per nucleon) and muons + * invisible particles (neutrinos) can be cut or not. + * + * @param outputArgs - optional arguments of TOutput writer + */ + template <typename... TArgs> ParticleCut(HEPEnergyType const eEleCut, HEPEnergyType const ePhoCut, - HEPEnergyType const eHadCut, HEPEnergyType const eMuCut, bool const inv); - - //! simple cut. hadrons and muons are cut by threshold. EM particles are all - //! discarded. - ParticleCut(HEPEnergyType const eHadCut, HEPEnergyType const euCut, bool const inv); + HEPEnergyType const eHadCut, HEPEnergyType const eMuCut, bool const inv, + TArgs&&... args); - //! simplest cut. all particles have same threshold. EM particles can be set to be - //! discarded altogether. - ParticleCut(HEPEnergyType const eCut, bool const em, bool const inv); + /** + * particle cut with kinetic energy thresholds for all particles. + * + * @param outputArgs - optional arguments of TOutput writer + */ + template <typename... TArgs> + ParticleCut(HEPEnergyType const eCut, bool const inv, TArgs&&... OutputArgs); - //! threshold for specific particles redefined. EM and invisible particles can be set - //! to be discarded altogether. + /** + * Threshold for specific particles redefined. EM and invisible particles can be set + * to be discarded altogether. + * + * @param outputArgs - optional arguments of TOutput writer + */ + template <typename... TArgs> ParticleCut(std::unordered_map<Code const, HEPEnergyType const> const& eCuts, - bool const em, bool const inv); + bool const inv, TArgs&&... outputArgs); + /** + * Cut particles which are secondaries from discrete processes. + * + * @tparam TStackView + */ template <typename TStackView> void doSecondaries(TStackView&); + /** + * Cut particles during contunuous processes (energy losses etc). + * + * @tparam TParticle + * @tparam TTrajectory + * @param vParticle + * @param vTrajectory + * @param limitFlag + * @return ProcessReturn + */ template <typename TParticle, typename TTrajectory> ProcessReturn doContinuous( TParticle& vParticle, TTrajectory const& vTrajectory, const bool limitFlag = false); // this is not used for ParticleCut + /** + * Limit on continuous step length imposed by ParticleCut: none. + * + * @tparam TParticle + * @tparam TTrajectory + * @return LengthType + */ template <typename TParticle, typename TTrajectory> LengthType getMaxStepLength(TParticle const&, TTrajectory const&) { return meter * std::numeric_limits<double>::infinity(); } - void printThresholds(); - void showResults(); // LCOV_EXCL_LINE - void reset(); + void printThresholds() const; - HEPEnergyType getElectronKineticECut() const { - return calculate_kinetic_energy_threshold(Code::Electron); - } - HEPEnergyType getPhotonKineticECut() const { - return calculate_kinetic_energy_threshold(Code::Photon); - } - HEPEnergyType getMuonKineticECut() const { - return calculate_kinetic_energy_threshold(Code::MuPlus); - } - HEPEnergyType getHadronKineticECut() const { - return calculate_kinetic_energy_threshold(Code::Proton); - } - //! returns total energy of particles that were removed by cut for invisible particles - HEPEnergyType getInvEnergy() const { return energy_invcut_; } - //! returns total energy of particles that were removed by cut in time - HEPEnergyType getTimeCutEnergy() const { return energy_timecut_; } - //! returns total energy of particles that were removed by cut in kinetic energy - HEPEnergyType getCutEnergy() const { return energy_cut_; } - //! returns total energy of particles that were removed by cut for electromagnetic - //! particles - HEPEnergyType getEmEnergy() const { return energy_emcut_; } - //! returns number of electromagnetic particles - unsigned int getNumberEmParticles() const { return em_count_; } - //! returns number of invisible particles - unsigned int getNumberInvParticles() const { return inv_count_; } + HEPEnergyType getElectronKineticECut() const { return cut_electrons_; } + HEPEnergyType getPhotonKineticECut() const { return cut_photons_; } + HEPEnergyType getMuonKineticECut() const { return cut_muons_; } + HEPEnergyType getHadronKineticECut() const { return cut_hadrons_; } + + //! get configuration of this node, for output + YAML::Node getConfig() const override; private: template <typename TParticle> @@ -99,22 +116,14 @@ namespace corsika { template <typename TParticle> bool isBelowEnergyCut(TParticle const&) const; - //! defines which particles are invisible, by default only neutrinos - bool isInvisible(Code const&) const; - private: - bool doCutEm_; + HEPEnergyType cut_electrons_; + HEPEnergyType cut_photons_; + HEPEnergyType cut_muons_; + HEPEnergyType cut_hadrons_; bool doCutInv_; - HEPEnergyType energy_cut_ = 0 * electronvolt; - HEPEnergyType energy_timecut_ = 0 * electronvolt; - HEPEnergyType energy_emcut_ = 0 * electronvolt; - HEPEnergyType energy_invcut_ = 0 * electronvolt; - unsigned int em_count_ = 0; - unsigned int inv_count_ = 0; - unsigned int energy_count_ = 0; - - HEPEnergyType energy_event_; // per event sum - }; + std::unordered_map<Code const, HEPEnergyType const> cuts_; + }; // namespace corsika } // namespace corsika diff --git a/corsika/modules/StackInspector.hpp b/corsika/modules/StackInspector.hpp index 6723ec3d9bedb0ff3097463133cc5d1a5e8dbe7c..6e723a6ced41be1e34a970991d3e55fc40c6af88 100644 --- a/corsika/modules/StackInspector.hpp +++ b/corsika/modules/StackInspector.hpp @@ -48,7 +48,7 @@ namespace corsika { bool ReportStack_; HEPEnergyType E0_; const HEPEnergyType dE_threshold_ = 1_eV; - decltype(std::chrono::system_clock::now()) StartTime_; + std::chrono::system_clock::time_point StartTime_; }; } // namespace corsika diff --git a/corsika/modules/TrackWriter.hpp b/corsika/modules/TrackWriter.hpp index 79e1dfc91bbb52d08300699fc3e6ee1e098feb8a..544fbd4346af91891c424a9e205443bf30539481 100644 --- a/corsika/modules/TrackWriter.hpp +++ b/corsika/modules/TrackWriter.hpp @@ -10,12 +10,25 @@ #include <corsika/framework/process/ContinuousProcess.hpp> #include <corsika/modules/writers/TrackWriterParquet.hpp> +#include <corsika/modules/writers/WriterOff.hpp> namespace corsika { - template <typename TOutputWriter = TrackWriterParquet> - class TrackWriter : public ContinuousProcess<TrackWriter<TOutputWriter>>, - public TOutputWriter { + /** + * @ingroup Modules + * @{ + * + * To write 3D track data to disk. + * + * Since the only sole purpose of this module is to generate track + * output on disk, the default output mode is not "WriterOff" but + * directly TrackWriterParquet. It can of course be changed. + * + * @tparam TOutput with default TrackWriterParquet + */ + + template <typename TOutput = TrackWriterParquet> + class TrackWriter : public ContinuousProcess<TrackWriter<TOutput>>, public TOutput { public: TrackWriter(); @@ -27,8 +40,12 @@ namespace corsika { LengthType getMaxStepLength(TParticle const&, TTrack const&); YAML::Node getConfig() const; + + private: }; + //! @} + } // namespace corsika #include <corsika/detail/modules/TrackWriter.inl> diff --git a/corsika/modules/conex/CONEXhybrid.hpp b/corsika/modules/conex/CONEXhybrid.hpp index 0d3878e474212cab10c6fefaf81ad8b23067d145..f20ec852889524003d901924f837a9dad5b28ceb 100644 --- a/corsika/modules/conex/CONEXhybrid.hpp +++ b/corsika/modules/conex/CONEXhybrid.hpp @@ -17,26 +17,63 @@ #include <corsika/framework/geometry/Vector.hpp> #include <corsika/media/ShowerAxis.hpp> +#include <corsika/modules/writers/SubWriter.hpp> +#include <corsika/modules/writers/EnergyLossWriter.hpp> +#include <corsika/modules/writers/LongitudinalWriter.hpp> + #include <corsika/modules/conex/CONEX_f.hpp> +/** + * @file CONEXhybrid.hpp + */ + namespace corsika { namespace conex { LengthType constexpr earthRadius{6371315 * meter}; } // namespace conex - class CONEXhybrid : public CascadeEquationsProcess<CONEXhybrid>, - public SecondariesProcess<CONEXhybrid> { + /** + * Access to the CONEX model. + * + * The fortran version of CONEX is interfaced. This is a SecondariesProcess since it + * skims the Stack for input particles (specific energies, or types) in order to fill + * internal CONEX source tables. CONEX is the called after the Cascade is finished to do + * a CascadeEquations step. This typically produces output data, and in addition can + * even generate new secondaries on the main Stack. + * + * Note that the output is processed by `SubWriter<TOutputE/N>`. The SubWriters have to + * be initialized with valid objects of type TOutputE/N at construction time.If no + * output is wished, `WriterOff` can be used for this purpose. + * + * @tparam TOutputE -- Output writer for dEdX data. + * @tparam TOutputN -- Output writer for particle number profile data. + */ + template <typename TOutputE, typename TOutputN> + class CONEXhybrid : public CascadeEquationsProcess<CONEXhybrid<TOutputE, TOutputN>>, + public SecondariesProcess<CONEXhybrid<TOutputE, TOutputN>>, + public SubWriter<TOutputE>, + public SubWriter<TOutputN> { public: - CONEXhybrid(Point center, ShowerAxis const& showerAxis, LengthType groundDist, - LengthType injectionHeight, HEPEnergyType primaryEnergy, PDGCode pdg); - /** - * Main entry point to pass new particle data towards CONEX. If a - * particles is selected for CONEX, it is removed from the CORSIKA - * 8 stack. + * Constructor with output sink specified. + * + * The output will be generated by the output sink. + * + * @param center center of earth. + * @param showerAxis shower axis to convert geometry to grammage. + * @param groundDist distance to ground. + * @param injectionHeight height of injection. + * @param primaryEnergy energy of primary particle + * @param pdg type of primary particle. + * @param outputE object to initialized SubWriter<TOutputE> + * @param outputN object to initialized SubWriter<TOutputN> */ + CONEXhybrid(Point const& center, ShowerAxis const& showerAxis, LengthType groundDist, + LengthType injectionHeight, HEPEnergyType primaryEnergy, PDGCode pdg, + TOutputE& outputE, TOutputN& outputN); + template <typename TStackView> void doSecondaries(TStackView&); @@ -47,7 +84,7 @@ namespace corsika { void initCascadeEquations(); /** - * Cascade equations are solved basoned on the data in the tables + * Cascade equations are solved basoned on the data in the tables. */ template <typename TStack> void doCascadeEquations(TStack& stack); @@ -58,12 +95,11 @@ namespace corsika { */ bool addParticle(Code pid, HEPEnergyType energy, HEPEnergyType mass, Point const& position, Vector<dimensionless_d> const& direction, - TimeType t); + TimeType t, double weight = 1); CoordinateSystemPtr const& getObserverCS() const { return conexObservationCS_; } - HEPEnergyType getEnergyEM() const; - void reset(); + YAML::Node getConfig() const final override; private: // data members @@ -81,7 +117,6 @@ namespace corsika { CoordinateSystemPtr const conexObservationCS_; //!< CONEX observation frame DirectionVector const x_sf_, y_sf_; //!< unit vectors of CONEX shower frame, z_sf is shower axis direction - HEPEnergyType energy_em_; }; } // namespace corsika diff --git a/corsika/modules/energy_loss/BetheBlochPDG.hpp b/corsika/modules/energy_loss/BetheBlochPDG.hpp index 5f606b82db38b532e52c24c2399911b194f658d9..4c2f56de71804458f204c9680a73c0950da5483a 100644 --- a/corsika/modules/energy_loss/BetheBlochPDG.hpp +++ b/corsika/modules/energy_loss/BetheBlochPDG.hpp @@ -12,7 +12,8 @@ #include <corsika/framework/geometry/Point.hpp> #include <corsika/framework/geometry/Vector.hpp> #include <corsika/framework/process/ContinuousProcess.hpp> -#include <corsika/media/ShowerAxis.hpp> + +#include <corsika/modules/writers/WriterOff.hpp> #include <map> @@ -33,20 +34,22 @@ namespace corsika { * */ - class BetheBlochPDG : public ContinuousProcess<BetheBlochPDG> { + template <typename TOutput = WriterOff> + class BetheBlochPDG : public ContinuousProcess<BetheBlochPDG<TOutput>>, public TOutput { using MeVgcm2 = decltype(1e6 * electronvolt / gram * square(1e-2 * meter)); public: - BetheBlochPDG(ShowerAxis const& showerAxis); + template <typename... TOutputArgs> + BetheBlochPDG(TOutputArgs&&... args); - /** clang-format-off + /** * Interface function of ContinuousProcess. * - * \param particle The particle to process in its current state - * \param track The trajectory in space of this particle, on which doContinuous should - * act - * \param limitFlag flag to identify, if BetheBlochPDG::getMaxStepLength is the + * @param particle The particle to process in its current state + * @param track The trajectory in space of this particle, on which doContinuous + *should act + * @param limitFlag flag to identify, if BetheBlochPDG::getMaxStepLength is the * globally limiting factor (or not) clang-format-on **/ template <typename TParticle, typename TTrajectory> @@ -54,10 +57,8 @@ namespace corsika { bool const limitFlag); template <typename TParticle, typename TTrajectory> - LengthType getMaxStepLength(TParticle const&, - TTrajectory const&) - const; //! limited by the energy threshold! By default the limit is the particle - //! rest mass, i.e. kinetic energy is zero + LengthType getMaxStepLength(TParticle const&, TTrajectory const&) const; + template <typename TParticle> static HEPEnergyType getBetheBloch(TParticle const&, const GrammageType); @@ -67,24 +68,7 @@ namespace corsika { template <typename TParticle> static HEPEnergyType getTotalEnergyLoss(TParticle const&, const GrammageType); - void showResults() const; - void reset(); - HEPEnergyType getEnergyLost() const { return energy_lost_; } - void printProfile() const; - HEPEnergyType getTotal() const; - - private: - template <typename TParticle> - void updateMomentum(TParticle&, HEPEnergyType Enew); - - template <typename TTrajectory> - void fillProfile(TTrajectory const&, HEPEnergyType); - - GrammageType const dX_ = 10_g / square(1_cm); // profile binning - GrammageType const dX_threshold_ = 0.0001_g / square(1_cm); - ShowerAxis const& shower_axis_; - HEPEnergyType energy_lost_ = HEPEnergyType::zero(); - std::vector<HEPEnergyType> profile_; // longitudinal profile + YAML::Node getConfig() const override; }; } // namespace corsika diff --git a/corsika/modules/proposal/ContinuousProcess.hpp b/corsika/modules/proposal/ContinuousProcess.hpp index bf421186c1e68c2522060f601fbc5116bfc2bc7b..ec25320177cb3adb8c560e6c8af54eb23606f070 100644 --- a/corsika/modules/proposal/ContinuousProcess.hpp +++ b/corsika/modules/proposal/ContinuousProcess.hpp @@ -29,9 +29,11 @@ namespace corsika::proposal { //! use of interpolation tables which are runtime intensive calculation, but can be //! reused by setting the \param PROPOSAL::InterpolationDef::path_to_tables variable. //! + template <typename TOutput = WriterOff> class ContinuousProcess - : public corsika::ContinuousProcess<proposal::ContinuousProcess>, - ProposalProcessBase { + : public corsika::ContinuousProcess<proposal::ContinuousProcess<TOutput>>, + public ProposalProcessBase, + public TOutput { struct Calculator { std::unique_ptr<PROPOSAL::Displacement> disp; @@ -41,8 +43,6 @@ namespace corsika::proposal { std::unordered_map<calc_key_t, Calculator, hash> calc; //!< Stores the displacement and scattering calculators. - HEPEnergyType energy_lost_ = 0 * electronvolt; - //! //! Build the displacement and scattering calculators and add it to calc. //! @@ -53,8 +53,8 @@ namespace corsika::proposal { //! Produces the continuous loss calculator for leptons based on nuclear //! compositions and stochastic description limited by the particle cut. //! - template <typename TEnvironment> - ContinuousProcess(TEnvironment const&); + template <typename TEnvironment, typename... TOutputArgs> + ContinuousProcess(TEnvironment const&, TOutputArgs&&...); //! //! Multiple Scattering of the lepton. Stochastic deflection is not yet taken into @@ -81,9 +81,10 @@ namespace corsika::proposal { template <typename TParticle, typename TTrack> LengthType getMaxStepLength(TParticle const&, TTrack const&); - void showResults() const; - void reset(); - HEPEnergyType getEnergyLost() const { return energy_lost_; } + /** + * Provide the config as YAML object to be stored on disk as output. + */ + YAML::Node getConfig() const; }; } // namespace corsika::proposal diff --git a/corsika/modules/writers/EnergyLossWriter.hpp b/corsika/modules/writers/EnergyLossWriter.hpp new file mode 100644 index 0000000000000000000000000000000000000000..56ff1da673cbe29996b1a15c96d947d879d39f94 --- /dev/null +++ b/corsika/modules/writers/EnergyLossWriter.hpp @@ -0,0 +1,173 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/output/BaseOutput.hpp> +#include <corsika/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/media/ShowerAxis.hpp> +#include <corsika/modules/writers/WriterOff.hpp> +#include <corsika/modules/writers/EnergyLossWriterParquet.hpp> + +#include <vector> +#include <array> + +/** + * @file EnergyLossWriter.hpp + */ + +namespace corsika { + + // clang-format-off + /** + * The energy loss writer can be used to pool several energy loss processes into one + * output file/stream. + * + * Typically many processes/modules can lead to energy losses in the shower. The + * EnergyLossWriter can be used in combination with the SubWriter class to collect all + * of them into a single output stream: + * + * \code {.cpp} + * # showerAxis must be a ShowerAxis object + * # the X binning can be specified. + * EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200}; + * # add to OutputManager: + * output.add("energyloss", dEdX); + * # add SubWriters, e.g. Bethe-Bloch: + * BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss{dEdX}; + * ... + * \endcode + * + * The EnergyLossWriter processes data on single-particle-level. The final output + * writer, e.g. EnergyLossWriterParquet, processes data on profile-level (bins in X). + * The default output option is parquet format. + * + * @tparam TOutput + */ + // clang-format-on + + /** + * Local helper namespace to store number and names of dEdX profile columns. + */ + namespace dEdX_output { + + /** + * Definition of longitudinal profile columns. + */ + enum class ProfileIndex { Total, Entries }; + + /** + * Number of columns (static). + */ + size_t constexpr NColumns = static_cast<int>(ProfileIndex::Entries); + + /** + * Names of columns in output. + */ + static std::array<char const*, NColumns> constexpr ProfileIndexNames{{"total"}}; + + /** + * Data type to store column data. + */ + typedef std::array<HEPEnergyType, NColumns> Profile; + + } // namespace dEdX_output + + /** + * The EnergyLossWriter can be used to pool the dEdX energy loss of several + * processes/modules into one output file/stream. + * + * Typically several processes/modules can lead to energy losses along the shower axis + * in the shower. The EnergyLossWriter can be used in combination with the SubWriter + * class to collect all of them into a single output stream: + * + * \code {.cpp} + * # showerAxis must be a ShowerAxis object + * # the X binning can be specified. + * EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200}; + * # add to OutputManager: + * output.add("energyloss", dEdX); + * # add SubWriters, e.g. BetheBlochPDG, CONEX: + * BetheBlochPDG<SubWriter<decltype(dEdX)>> long{dEdX}; + * CONEXhybrid<SubWriter<decltype(dEdX)>> conex{..., dEdX}; + * ... + * \endcode + * + * The default output option is parquet format. + * + * @tparam TOutput + */ + + template <typename TOutput = EnergyLossWriterParquet<dEdX_output::NColumns>> + class EnergyLossWriter : public TOutput { + + public: + /** + * Construct a new writer. + */ + EnergyLossWriter(ShowerAxis const& axis, + GrammageType dX = 10_g / square(1_cm), // profile binning + unsigned int const nBins = 200, // number of bins + GrammageType dX_threshold = 0.0001_g / + square(1_cm)); // ignore too short tracks + + void startOfLibrary(boost::filesystem::path const& directory) final override; + + void startOfShower(unsigned int const showerId) final override; + + void endOfShower(unsigned int const showerId) final override; + + void endOfLibrary() final override; + + /** + * Add continuous energy loss. + */ + template <typename TTrack> + void write(TTrack const& track, Code const PID, HEPEnergyType const dE); + + /** + * Add localized energy loss. + */ + void write(Point const& point, Code const PID, HEPEnergyType const dE); + + /** + * Add binned energy loss. + */ + void write(GrammageType const Xstart, GrammageType const Xend, Code const PID, + HEPEnergyType const dE); + + /** + * Get total observed energy loss. + * + * @return HEPEnergyType The total energy. + */ + HEPEnergyType getEnergyLost() const; + + /** + * Return a summary. + */ + YAML::Node getSummary() const; + + /** + * Return the configuration of this output. + */ + YAML::Node getConfig() const; + + private: + ShowerAxis const& showerAxis_; ///< conversion between geometry and grammage + GrammageType dX_; ///< binning of profile. + size_t nBins_; ///< number of profile bins. + GrammageType dX_threshold_; ///< too short tracks are discarded. + std::vector<dEdX_output::Profile> profile_; // longitudinal profile + + }; // namespace corsika + +} // namespace corsika + +#include <corsika/detail/modules/writers/EnergyLossWriter.inl> diff --git a/corsika/modules/writers/EnergyLossWriterParquet.hpp b/corsika/modules/writers/EnergyLossWriterParquet.hpp new file mode 100644 index 0000000000000000000000000000000000000000..ef58f56c6cdedd784dbefb62a0ccd21d487f48b7 --- /dev/null +++ b/corsika/modules/writers/EnergyLossWriterParquet.hpp @@ -0,0 +1,80 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/output/BaseOutput.hpp> +#include <corsika/output/ParquetStreamer.hpp> +#include <corsika/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/media/ShowerAxis.hpp> + +#include <array> + +namespace corsika { + + /** + * The actual writer to save dEdX data to disk. + * + * The purpose of this class is not to collect single-particle-level energy loss data. + * But to write entire binned profiles to disk at the end of a shower. + * To fill the shower data, you have to use EnergyLossWriter in combination with + * SubWriter. The EnergyLossWriterParquet is the default output mode of the + * EnergyLossWriter. + * + * @tparam NColumn -- the number of columns written to output. column names and data + * must be provided consistently. + */ + + template <size_t NColumns> + class EnergyLossWriterParquet : public BaseOutput { + + public: + /** + * Construct a new writer. + */ + EnergyLossWriterParquet(std::array<const char*, NColumns> const& colNames); + + /** + * Called at the start of each library. + */ + void startOfLibrary(boost::filesystem::path const& directory) override; + + /** + * Called at the start of each shower. + */ + void startOfShower(unsigned int const showerId) override; + + /** + * Called at the end of each shower. + */ + void endOfShower(unsigned int const showerId) override; + + /** + * Called at the end of each library. + * + * This must also increment the run number since we override + * the default behaviour of BaseOutput. + */ + void endOfLibrary() override; + + /** + * Write energy lost to the file. + */ + void write(unsigned int const showerId, GrammageType const grammage, + std::array<HEPEnergyType, NColumns> const& data); + + private: + std::array<const char*, NColumns> columns_; ///< column names + ParquetStreamer output_; ///< The primary output file. + + }; // namespace corsika + +} // namespace corsika + +#include <corsika/detail/modules/writers/EnergyLossWriterParquet.inl> diff --git a/corsika/modules/writers/LongitudinalProfileWriterParquet.hpp b/corsika/modules/writers/LongitudinalProfileWriterParquet.hpp new file mode 100644 index 0000000000000000000000000000000000000000..bba7000ff30f73260740c66ab35a3539f97259b6 --- /dev/null +++ b/corsika/modules/writers/LongitudinalProfileWriterParquet.hpp @@ -0,0 +1,76 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/output/BaseOutput.hpp> +#include <corsika/output/ParquetStreamer.hpp> +#include <corsika/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> + +#include <array> + +namespace corsika { + + /** + * The actual writer to save longitudinal profile data to disk. + * + * The purpose of this class is not to collect single-particle-level profile data, + * but to write entire binned profiles to disk at the end of a shower. + * To fill the shower data, you have to use ProfileWriter in combination with + * SubWriter. The LongitudinalProfileWriterParquet is the default output mode of the + * LongitudinalWriter. + */ + + template <size_t NColumns> + class LongitudinalProfileWriterParquet : public BaseOutput { + + public: + /** + * Construct a new writer. + */ + LongitudinalProfileWriterParquet(std::array<const char*, NColumns> const& colNames); + + /** + * Called at the start of each library. + */ + void startOfLibrary(boost::filesystem::path const& directory) override; + + /** + * Called at the start of each shower. + */ + void startOfShower(unsigned int const showerId) override; + + /** + * Called at the end of each shower. + */ + void endOfShower(unsigned int const showerId) override; + + /** + * Called at the end of each library. + * + * This must also increment the run number since we override + * the default behaviour of BaseOutput. + */ + void endOfLibrary() override; + + /** + * Add profile to disk. + */ + void write(unsigned int const showerId, GrammageType const grammage, + std::array<double, NColumns> const& data); + + private: + std::array<const char*, NColumns> columns_; //!< column names + ParquetStreamer output_; //!< The parquet streamer for this process. + + }; // namespace corsika + +} // namespace corsika + +#include <corsika/detail/modules/writers/LongitudinalProfileWriterParquet.inl> diff --git a/corsika/modules/writers/LongitudinalWriter.hpp b/corsika/modules/writers/LongitudinalWriter.hpp new file mode 100644 index 0000000000000000000000000000000000000000..ac00dc3f0b25a98d1e138f7ee75bd8adf3033869 --- /dev/null +++ b/corsika/modules/writers/LongitudinalWriter.hpp @@ -0,0 +1,145 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/output/BaseOutput.hpp> +#include <corsika/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/media/ShowerAxis.hpp> +#include <corsika/modules/writers/WriterOff.hpp> +#include <corsika/modules/writers/LongitudinalProfileWriterParquet.hpp> + +#include <vector> +#include <array> + +/** + * @file LongitudinalWriter.hpp + */ + +namespace corsika { + + /** + * Local helper namespace to store number and names of particle number profile columns. + */ + namespace number_profile { + + /** + * Definition of longitudinal profile columns. + */ + enum class ProfileIndex { + Charged, + Hadron, + Photon, + Electron, + Positron, + MuPlus, + MuMinus, + Entries + }; + + /** + * Number of columns (static). + */ + size_t constexpr NColumns = static_cast<size_t>(ProfileIndex::Entries); + + /** + * Names of columns in output. + */ + static std::array<char const*, NColumns> constexpr ProfileIndexNames{ + {"charged", "hadron", "photon", "electron", "positron", "muplus", "muminus"}}; + + /** + * Data type to store column data. + */ + typedef std::array<double, NColumns> ProfileData; + } // namespace number_profile + + // clang-format-off + /** + * The LongitudinalWriter can be used to pool the particle counts of several + * longitudinal profile processes into one output file/stream. + * + * Typically several processes/modules can lead to particle counts along the shower axis + * in the shower. The LongitudianalWriter can be used in combination with the SubWriter + * class to collect all of them into a single output stream: + * + * \code {.cpp} + * # showerAxis must be a ShowerAxis object + * # the X binning can be specified. + * LongitudinalWriter profile{showerAxis, 10_g / square(1_cm), 200}; + * # add to OutputManager: + * output.add("profile", profile); + * # add SubWriters, e.g. LongitudinalProfile, CONEX: + * LongitudinalProfile<SubWriter<decltype(profile)>> long{profile}; + * CONEXhybrid<SubWriter<decltype(profile)>> conex{..., profile}; + * ... + * \endcode + * + * The default output option is parquet format. + * + * @tparam TOutput + */ + // clang-format-on + + template <typename TOutput = LongitudinalProfileWriterParquet<number_profile::NColumns>> + class LongitudinalWriter : public TOutput { + + public: + /** + * Construct a new writer. + */ + LongitudinalWriter(ShowerAxis const& axis, + GrammageType dX = 10_g / square(1_cm), // profile binning + size_t const nBins = 200); // number of bins + + void startOfLibrary(boost::filesystem::path const& directory) final override; + + void startOfShower(unsigned int const showerId) final override; + + void endOfShower(unsigned int const showerId) final override; + + void endOfLibrary() final override; + + /** + * Add continuous profile. + */ + template <typename TTrack> + void write(TTrack const& track, Code const pid, double const weight); + + /** + * Add binned profile. + */ + void write(GrammageType const Xstart, GrammageType const Xend, Code const pid, + double const weight); + + /** + * Return a summary. + */ + YAML::Node getSummary() const; + + /** + * Return the configuration of this output. + */ + YAML::Node getConfig() const; + + number_profile::ProfileData const& getProfile( + number_profile::ProfileIndex index) const { + return profile_.at(static_cast<int>(index)); + } + + private: + ShowerAxis const& showerAxis_; ///< conversion between geometry and grammage + GrammageType dX_; ///< binning of profile. + size_t nBins_; ///< number of profile bins. + std::vector<number_profile::ProfileData> profile_; // longitudinal profile + }; + +} // namespace corsika + +#include <corsika/detail/modules/writers/LongitudinalWriter.inl> diff --git a/corsika/modules/writers/ObservationPlaneWriterParquet.hpp b/corsika/modules/writers/ParticleWriterParquet.hpp similarity index 50% rename from corsika/modules/writers/ObservationPlaneWriterParquet.hpp rename to corsika/modules/writers/ParticleWriterParquet.hpp index 5963aa1002377749c3c2a9876d2a5bee8887b074..12451f4eddbcb55e4e14ca5a1cc21a9e124b9844 100644 --- a/corsika/modules/writers/ObservationPlaneWriterParquet.hpp +++ b/corsika/modules/writers/ParticleWriterParquet.hpp @@ -15,27 +15,28 @@ namespace corsika { - class ObservationPlaneWriterParquet : public BaseOutput { - - ParquetStreamer output_; ///< The primary output file. + class ParticleWriterParquet : public BaseOutput { public: /** * Construct an ObservationPlane. - * - * @param name The name of this output. */ - ObservationPlaneWriterParquet(); + ParticleWriterParquet(); /** * Called at the start of each library. */ void startOfLibrary(boost::filesystem::path const& directory) final override; + /** + * Called at the beginning of each shower. + */ + void startOfShower(unsigned int const showerId) final override; + /** * Called at the end of each shower. */ - void endOfShower() final override; + void endOfShower(unsigned int const showerId) final override; /** * Called at the end of each library. @@ -45,16 +46,36 @@ namespace corsika { */ void endOfLibrary() final override; - protected: /** - * Write a particle to the file. + * Write a PDG/corsika::Code particle to the file. */ void write(Code const& pid, units::si::HEPEnergyType const& energy, units::si::LengthType const& x, units::si::LengthType const& y, - units::si::TimeType const& t); + units::si::LengthType const& z, const double weight); + + /** + * Return collected library-level summary for output. + */ + YAML::Node getSummary() const final override; + + /** + * If plane is absorbing particles: return the total energy absorbed. + */ + HEPEnergyType getEnergyGround() const { return totalEnergy_; } + + private: + ParquetStreamer output_; ///< The primary output file. + unsigned int showerId_; ///< current shower Id + + double countHadrons_ = 0; ///< count hadrons hitting plane + double countMuons_ = 0; ///< count muons hitting plane + double countEM_ = 0; ///< count EM particles hitting plane. + double countOthers_ = 0; ///< count othe types of particles hitting plane + + HEPEnergyType totalEnergy_; ///< energy absorbed in ground. - }; // class ObservationPlaneWriterParquet + }; // class ParticleWriterParquet } // namespace corsika -#include <corsika/detail/modules/writers/ObservationPlaneWriterParquet.inl> +#include <corsika/detail/modules/writers/ParticleWriterParquet.inl> diff --git a/corsika/modules/writers/SubWriter.hpp b/corsika/modules/writers/SubWriter.hpp new file mode 100644 index 0000000000000000000000000000000000000000..d75017715114c97f9b845b4caeb63b8fbeae94cd --- /dev/null +++ b/corsika/modules/writers/SubWriter.hpp @@ -0,0 +1,38 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/output/Configurable.hpp> + +#include <utility> + +namespace corsika { + + template <typename TOutput> + class SubWriter : public Configurable { + + TOutput& output_; + + public: + SubWriter(TOutput& output) + : output_(output) {} + + virtual ~SubWriter() {} + + /* + * Write data to the underlying writer. + */ + template <typename... TArgs> + void write(TArgs&&... args) { + output_.write(std::forward<TArgs>(args)...); + } + + }; // class SubWriter + +} // namespace corsika diff --git a/corsika/modules/writers/TrackWriterParquet.hpp b/corsika/modules/writers/TrackWriterParquet.hpp index 8473cbf7f4a17c649ddae2f5b20ffaf7dc400616..2d7ac7b2824804d3634dce11395b938ebe7cd566 100644 --- a/corsika/modules/writers/TrackWriterParquet.hpp +++ b/corsika/modules/writers/TrackWriterParquet.hpp @@ -31,10 +31,15 @@ namespace corsika { */ void startOfLibrary(boost::filesystem::path const& directory) final override; + /** + * Called at the start of each shower. + */ + void startOfShower(unsigned int const showerId) final override; + /** * Called at the end of each shower. */ - void endOfShower() final override; + void endOfShower(unsigned int const showerId) final override; /** * Called at the end of each library. @@ -44,16 +49,16 @@ namespace corsika { */ void endOfLibrary() final override; - protected: /** * Write a track to the file. */ - void write(Code const& pid, units::si::HEPEnergyType const& energy, - QuantityVector<length_d> const& start, TimeType const& t_start, - QuantityVector<length_d> const& end, TimeType const& t_end); + void write(Code const pid, HEPEnergyType const energy, double const weight, + QuantityVector<length_d> const& start, TimeType const t_start, + QuantityVector<length_d> const& end, TimeType const t_end); private: ParquetStreamer output_; ///< The primary output file. + unsigned int showerId_; ///< event Id counter }; // class TrackWriterParquet diff --git a/corsika/modules/writers/WriterOff.hpp b/corsika/modules/writers/WriterOff.hpp new file mode 100644 index 0000000000000000000000000000000000000000..3d752cc9a80de7a85b0705c3d6dac9dac9d8722d --- /dev/null +++ b/corsika/modules/writers/WriterOff.hpp @@ -0,0 +1,50 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <corsika/output/BaseOutput.hpp> + +namespace corsika { + + /** + * Generic class to switch off any output. + * + * The 'write' method is catch-all and does nothing. + */ + + class WriterOff : public BaseOutput { + + public: + /** + * The purpose of this catch-all constructor is to discard all parameters. + * @tparam TArgs + */ + template <typename... TArgs> + WriterOff(TArgs&&...) {} + + virtual ~WriterOff() {} + + void startOfLibrary(boost::filesystem::path const&) override {} + + void endOfShower(unsigned int const) override {} + + void endOfLibrary() override {} + + /** + * The purpose of this catch-all method is to discard all data. + * @tparam TArgs + */ + template <typename... TArgs> + void write(TArgs&&...) {} + + virtual YAML::Node getConfig() const { return YAML::Node(); } + + }; // class WriterOff + +} // namespace corsika diff --git a/corsika/output/BaseOutput.hpp b/corsika/output/BaseOutput.hpp index 0024c5b7a918e4705fbbadd3bdcdd59dfcab0449..7d9165dc107321ee42212d4d4cb24008f74e41a7 100644 --- a/corsika/output/BaseOutput.hpp +++ b/corsika/output/BaseOutput.hpp @@ -7,8 +7,9 @@ */ #pragma once +#include <corsika/framework/core/Logging.hpp> +#include <corsika/output/Configurable.hpp> #include <boost/filesystem.hpp> - #include <yaml-cpp/yaml.h> namespace corsika { @@ -17,10 +18,11 @@ namespace corsika { * This is the base class for all outputs so that they * can be stored in homogeneous containers. */ - class BaseOutput { + class BaseOutput : public Configurable { protected: - BaseOutput(); + BaseOutput() = default; + virtual ~BaseOutput() = default; public: /** @@ -30,33 +32,37 @@ namespace corsika { /** * Called at the start of each event/shower. + * + * @param showerId Shower counter. */ - virtual void startOfShower() {} + virtual void startOfShower(unsigned int const /*showerId*/) {} /** * Called at the end of each event/shower. + * + * @param showerId Shower counter. */ - virtual void endOfShower() = 0; + virtual void endOfShower(unsigned int const showerId) = 0; // LCOV_EXCL_LINE /** * Called at the end of each run. */ - virtual void endOfLibrary() = 0; + virtual void endOfLibrary() = 0; // LCOV_EXCL_LINE /** - * Get the configuration of this output. + * Flag to indicate readiness. */ - virtual YAML::Node getConfig() const = 0; + bool isInit() const { return is_init_; } /** - * Get any summary information for the entire library. + * The output logger. */ - virtual YAML::Node getSummary() { return YAML::Node(); } + static auto getLogger() { return logger_; } /** - * Flag to indicate readiness. + * Provide YAML Summary for this BaseOutput. */ - bool isInit() const { return is_init_; } + virtual YAML::Node getSummary() const { return YAML::Node(); } protected: /** @@ -64,12 +70,9 @@ namespace corsika { */ void setInit(bool const v) { is_init_ = v; } - int shower_{0}; ///< The current event number. - private: - bool is_init_{false}; ///< flag to indicate readiness + bool is_init_{false}; ///< flag to indicate readiness + inline static auto logger_{get_logger("output")}; ///< A custom logger. }; } // namespace corsika - -#include <corsika/detail/output/BaseOutput.inl> diff --git a/corsika/output/Configurable.hpp b/corsika/output/Configurable.hpp new file mode 100644 index 0000000000000000000000000000000000000000..7259d4464e700e65fe217be1b8e5c9b373191f6a --- /dev/null +++ b/corsika/output/Configurable.hpp @@ -0,0 +1,34 @@ +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ +#pragma once + +#include <corsika/framework/core/Logging.hpp> +#include <corsika/output/Configurable.hpp> +#include <boost/filesystem.hpp> +#include <yaml-cpp/yaml.h> + +namespace corsika { + + /** + * This is the base class for all classes that have + * YAML representations of their configurations. + */ + class Configurable { + + protected: + Configurable() = default; + virtual ~Configurable() = default; + + public: + /** + * Provide YAML configuration for this BaseOutput. + */ + virtual YAML::Node getConfig() const = 0; // LCOV_EXCL_LINE + }; + +} // namespace corsika diff --git a/corsika/output/DummyOutputManager.hpp b/corsika/output/DummyOutputManager.hpp index fc5cf4647d5daa227fceceda5961bb051a840058..e7bc1fa76645fad68362e9f1e1b0335d8018f17d 100644 --- a/corsika/output/DummyOutputManager.hpp +++ b/corsika/output/DummyOutputManager.hpp @@ -12,7 +12,7 @@ namespace corsika { /*! * An output manager that does nothing. */ - class DummyOutputManager final { + class DummyOutputManager { public: /** diff --git a/corsika/output/NoOutput.hpp b/corsika/output/NoOutput.hpp index 788149b7047581f57398822fed2c1570189eedb4..cda72d0f1647b9ed79fa211c72548f70599f5810 100644 --- a/corsika/output/NoOutput.hpp +++ b/corsika/output/NoOutput.hpp @@ -30,12 +30,12 @@ namespace corsika { /** * Called at the start of each event/shower. */ - void startOfShower() final override {} + void startOfShower(unsigned int const) final override {} /** * Called at the end of each event/shower. */ - void endOfShower() final override {} + void endOfShower(unsigned int const) final override {} /** * Called at the end of each run. @@ -50,7 +50,7 @@ namespace corsika { /** * Get any summary information for the entire library. */ - YAML::Node getSummary() final override { return YAML::Node(); }; + YAML::Node getSummary() const final override { return YAML::Node(); }; protected: void write(Code const&, units::si::HEPEnergyType const&, units::si::LengthType const&, @@ -58,5 +58,3 @@ namespace corsika { }; } // namespace corsika - -#include <corsika/detail/output/BaseOutput.inl> diff --git a/corsika/output/OutputManager.hpp b/corsika/output/OutputManager.hpp index c3e570749a0c3e04702a65cadb95b13a0900461b..8cc9a63600355507d6a9c2188a19b41a1a50c70a 100644 --- a/corsika/output/OutputManager.hpp +++ b/corsika/output/OutputManager.hpp @@ -10,15 +10,16 @@ #include <chrono> #include <string> #include <boost/filesystem.hpp> -#include <corsika/output/BaseOutput.hpp> #include <corsika/framework/core/Logging.hpp> +#include <corsika/output/BaseOutput.hpp> +#include <corsika/output/YAMLStreamer.hpp> namespace corsika { /*! * Manages CORSIKA 8 output streams. */ - class OutputManager final { + class OutputManager : public YAMLStreamer { /** * Indicates the current state of this manager. @@ -30,65 +31,51 @@ namespace corsika { LibraryFinished, }; - OutputState state_{OutputState::NoInit}; ///< The current state of this manager. - std::string const name_; ///< The name of this simulation file. - boost::filesystem::path const root_; ///< The top-level directory for the output. - int count_{0}; ///< The current ID of this shower. - std::chrono::time_point<std::chrono::system_clock> const start_time{ - std::chrono::system_clock::now()}; ///< The time the manager is created. - inline static auto logger{get_logger("output")}; ///< A custom logger. - /** - * The outputs that have been registered with us. - */ - std::map<std::string, std::reference_wrapper<BaseOutput>> outputs_; - + public: /** - * Write a YAML-node to a file. + * Construct an OutputManager instance with a name in a given directory. + * + * @param name The name of this output collection. + * @param dir The directory where the output directory will be stored. */ - void writeNode(YAML::Node const& node, boost::filesystem::path const& path) const; + OutputManager(std::string const& name, boost::filesystem::path const& dir); /** - * Write the top-level config of this simulation. + * Handle graceful closure of the outputs upon destruction. */ - void writeTopLevelConfig() const; + ~OutputManager(); - /** - * Initialize the "registered" output with a given name. - */ - void initOutput(std::string const& name) const; + template <typename TOutput> + void add(std::string const& name, TOutput& output); /** - * Write the top-level summary of this library. + * Produces the summary YAML. + * + * @return YAML::Node */ - void writeTopLevelSummary() const; + YAML::Node getSummary() const; - public: /** - * Construct an OutputManager instance with a name in a given directory. + * Produces the config YAML. * - * @param name The name of this output collection. - * @param dir The directory where the output directory will be stored. + * @return YAML::Node */ - OutputManager(std::string const& name, boost::filesystem::path const& dir); + YAML::Node getConfig() const; + private: /** - * Handle graceful closure of the outputs upon destruction. + * Write the top-level config of this simulation. */ - ~OutputManager(); + void writeConfig() const; /** - * Register an existing output to this manager. - * - * @param name The unique name of this output. - * @param output The output module. + * Write the top-level summary of this library. */ - template <typename TOutput> - void add(std::string const& name, TOutput& output); + void writeSummary() const; + public: /** * Called at the start of each library. - * - * This iteratively calls startOfLibrary on each registered output. */ void startOfLibrary(); @@ -110,6 +97,24 @@ namespace corsika { */ void endOfLibrary(); + /** + * Return current event number. + */ + int getEventId() const; + + private: + boost::filesystem::path root_; ///< The unique output directory. + OutputState state_{OutputState::NoInit}; ///< The current state of this manager. + std::string const name_; ///< The name of this simulation file. + int count_{0}; ///< The current ID of this shower. + std::chrono::time_point<std::chrono::system_clock> const start_time{ + std::chrono::system_clock::now()}; ///< The time the manager is created. + inline static auto logger_{get_logger("output")}; ///< A custom logger. + /** + * The outputs that have been registered here. + */ + std::map<std::string, std::reference_wrapper<BaseOutput>> outputs_; + }; // class OutputManager } // namespace corsika diff --git a/corsika/output/YAMLStreamer.hpp b/corsika/output/YAMLStreamer.hpp new file mode 100644 index 0000000000000000000000000000000000000000..e0b333ccc22382622190a70a046bcfe5a28362c2 --- /dev/null +++ b/corsika/output/YAMLStreamer.hpp @@ -0,0 +1,30 @@ + +/* + * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#pragma once + +#include <boost/filesystem.hpp> +#include <yaml-cpp/yaml.h> +#include <string> + +namespace corsika { + + /** + * This class automates the construction of simple tabular + * YAML files using the YAML::StreamWriter. + */ + class YAMLStreamer { + public: + inline void writeYAML(YAML::Node const& node, + boost::filesystem::path const& path) const; + + }; // class YAMLStreamer +} // namespace corsika + +#include <corsika/detail/output/YAMLStreamer.inl> diff --git a/corsika/setup/SetupStack.hpp b/corsika/setup/SetupStack.hpp index 305089ec5a7577d8816fda6d8004cf0fcfcd52f5..9028c2ab175127d07f3c430bbbba5edb5c82722d 100644 --- a/corsika/setup/SetupStack.hpp +++ b/corsika/setup/SetupStack.hpp @@ -29,9 +29,9 @@ namespace corsika::setup { #else // WITH_HISTORY /* - * the version without history + * the version without history (and geometry data and weights) */ - using Stack = detail::StackWithGeometry; + using Stack = detail::StackWithWeight; #endif diff --git a/corsika/stack/WeightStackExtension.hpp b/corsika/stack/WeightStackExtension.hpp index b2e6081f2638c0fe116e5ad87d18af0502e71605..bbe636db1537ea357c34cd2ae9d75fc548f98c67 100644 --- a/corsika/stack/WeightStackExtension.hpp +++ b/corsika/stack/WeightStackExtension.hpp @@ -23,7 +23,7 @@ namespace corsika::weights { * Corresponding defintion of a stack-readout object, the iteractor * dereference operator will deliver access to these function * defintion of a stack-readout object, the iteractor dereference - * operator will deliver access to these function + * operator will deliver access to these function. */ /** @@ -53,7 +53,7 @@ namespace corsika::weights { /** * @class WeightData * - * definition of stack-data object to store geometry information + * definition of stack-data object to store geometry information. */ class WeightData { diff --git a/documentation/tracking.rst b/documentation/tracking.rst new file mode 100644 index 0000000000000000000000000000000000000000..12227b3a8a0695b85ddbef7b32425c1b9376e551 --- /dev/null +++ b/documentation/tracking.rst @@ -0,0 +1,17 @@ +Tracking +======== + +Tracking is the process of moving a kinetic object from one point A to another point B. +Physical forces and laws should be applied explicitly or numerically in this process. + +There are two main aspects of tracking: + +1. The actual proposal of a new point B, thus, a movement. +2. The description of the movement within a 3D geometry with boundaries and different volumina. + +Besides this, there furthermore is also the determination of the effects of the +movement (e.g. magnetic fields, refractivity), or the accumulation of local +quantities (e.g. grammage). +For these purposes typically numerical integration routines are needed. + + diff --git a/examples/boundary_example.cpp b/examples/boundary_example.cpp index a11da70c6559c8d593a75dcaad2fe006b0c67927..bffda918f64b3fc7ef79193f6470f9d104f9e711 100644 --- a/examples/boundary_example.cpp +++ b/examples/boundary_example.cpp @@ -118,7 +118,7 @@ int main() { // setup processes, decays and interactions setup::Tracking tracking; - ParticleCut cut(50_GeV, true, true); + ParticleCut cut(50_GeV, true); TrackWriter trackWriter; output.add("tracks", trackWriter); // register TrackWriter @@ -143,10 +143,7 @@ int main() { double const theta = distTheta(rng); double const phi = distPhi(rng); - auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { - return sqrt((Elab - m) * (Elab + m)); - }; - HEPMomentumType P0 = elab2plab(E0, mass); + HEPMomentumType const P0 = calculate_momentum(E0, mass); auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) { return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi), -ptot * cos(theta)); @@ -169,16 +166,9 @@ int main() { // define air shower object, run simulation Cascade EAS(env, tracking, sequence, output, stack); - output.startOfShower(); + output.startOfLibrary(); EAS.run(); - output.endOfShower(); - - CORSIKA_LOG_INFO("Result: E0={}GeV", E0 / 1_GeV); - cut.showResults(); - [[maybe_unused]] const HEPEnergyType Efinal = - (cut.getCutEnergy() + cut.getInvEnergy() + cut.getEmEnergy()); - CORSIKA_LOG_INFO("Total energy (GeV): {} relative difference (%): {}", Efinal / 1_GeV, - (Efinal / E0 - 1.) * 100); - output.endOfLibrary(); + + CORSIKA_LOG_INFO("Done"); } diff --git a/examples/cascade_example.cpp b/examples/cascade_example.cpp index 2f0d254472a477890d382ffe34ab56d6c6dabad8..b43428ae50e482f0760b4bb133bc0cd4fbfa1c14 100644 --- a/examples/cascade_example.cpp +++ b/examples/cascade_example.cpp @@ -17,6 +17,8 @@ #include <corsika/framework/utility/CorsikaFenv.hpp> #include <corsika/output/OutputManager.hpp> +#include <corsika/modules/writers/SubWriter.hpp> +#include <corsika/modules/writers/EnergyLossWriter.hpp> #include <corsika/media/Environment.hpp> #include <corsika/media/HomogeneousMedium.hpp> @@ -58,9 +60,9 @@ int main() { logging::set_level(logging::level::info); - std::cout << "cascade_example" << std::endl; + CORSIKA_LOG_INFO("cascade_example"); - const LengthType height_atmosphere = 112.8_km; + LengthType const height_atmosphere = 112.8_km; feenableexcept(FE_INVALID); // initialize random number sequence(s) @@ -109,26 +111,24 @@ int main() { OutputManager output("cascade_outputs"); - ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env}; + theta *= M_PI / 180.; + phi *= M_PI / 180.; + DirectionVector const direction( + rootCS, {sin(theta) * cos(phi), sin(theta) * sin(phi), -cos(theta)}); + + ShowerAxis const showerAxis{injectionPos, direction * 100_km, env}; + EnergyLossWriter dEdX{showerAxis}; + output.add("energyloss", dEdX); { - auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { - return sqrt((Elab - m) * (Elab + m)); - }; - HEPMomentumType P0 = elab2plab(E0, mass); - auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) { - return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi), - -ptot * cos(theta)); - }; - auto const [px, py, pz] = - momentumComponents(theta / 180. * M_PI, phi / 180. * M_PI, P0); - auto plab = MomentumVector(rootCS, {px, py, pz}); - cout << "input particle: " << beamCode << endl; - cout << "input angles: theta=" << theta << " phi=" << phi << endl; - cout << "input momentum: " << plab.getComponents() / 1_GeV << endl; + HEPMomentumType const P0 = calculate_momentum(E0, mass); + auto plab = direction * P0; + CORSIKA_LOG_INFO("input particle: {}", beamCode); + CORSIKA_LOG_INFO("input angles: theta={} phi={}", theta, phi); + CORSIKA_LOG_INFO("input momentum: {}", plab.getComponents() / 1_GeV); stack.addParticle(std::make_tuple( - beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)), - plab.normalized(), injectionPos, 0_ns)); + beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)), direction, + injectionPos, 0_ns)); } // setup processes, decays and interactions @@ -140,14 +140,14 @@ int main() { corsika::sibyll::Interaction sibyll; corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env); corsika::sibyll::Decay decay; + // cascade with only HE model ==> HE cut - ParticleCut cut(80_GeV, true, true); + ParticleCut<SubWriter<decltype(dEdX)>> cut(80_GeV, true, dEdX); + BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss{dEdX}; TrackWriter trackWriter; output.add("tracks", trackWriter); // register TrackWriter - BetheBlochPDG eLoss{showerAxis}; - // assemble all processes into an ordered process list auto sequence = make_sequence(stackInspect, make_sequence(sibyllNuc, sibyll), decay, eLoss, cut, trackWriter); @@ -155,20 +155,17 @@ int main() { // define air shower object, run simulation Cascade EAS(env, tracking, sequence, output, stack); - output.startOfShower(); + output.startOfLibrary(); EAS.run(); - output.endOfShower(); - - eLoss.printProfile(); // print longitudinal profile - - cut.showResults(); - const HEPEnergyType Efinal = - cut.getCutEnergy() + cut.getInvEnergy() + cut.getEmEnergy(); - cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl - << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl; - cout << "total dEdX energy (GeV): " << eLoss.getTotal() / 1_GeV << endl - << "relative difference (%): " << eLoss.getTotal() / E0 * 100 << endl; - cut.reset(); - output.endOfLibrary(); + + const HEPEnergyType Efinal = dEdX.getEnergyLost(); + CORSIKA_LOG_INFO( + "\n" + "total cut energy (GeV) : {}\n" + "relative difference (%): {}\n" + "total dEdX energy (GeV): {}\n" + "relative difference (%): {}\n", + Efinal / 1_GeV, (Efinal / E0 - 1) * 100, dEdX.getEnergyLost() / 1_GeV, + dEdX.getEnergyLost() / E0 * 100); } diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp index 9934211b74786c9c9b0f7cdf3c70d5cf15c1005c..7f41162b2feaf3efa785711697d17243a0e6ca89 100644 --- a/examples/cascade_proton_example.cpp +++ b/examples/cascade_proton_example.cpp @@ -17,6 +17,8 @@ #include <corsika/framework/utility/CorsikaFenv.hpp> #include <corsika/output/OutputManager.hpp> +#include <corsika/modules/writers/SubWriter.hpp> +#include <corsika/modules/writers/EnergyLossWriter.hpp> #include <corsika/media/Environment.hpp> #include <corsika/media/HomogeneousMedium.hpp> @@ -124,7 +126,13 @@ int main() { RNGManager<>::getInstance().registerRandomStream("pythia"); corsika::pythia8::Interaction pythia; corsika::pythia8::Decay decay; - ParticleCut cut(60_GeV, true, true); + + ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env}; + EnergyLossWriter dEdX{showerAxis}; + output.add("energyloss", dEdX); + + BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss{dEdX}; + ParticleCut<SubWriter<decltype(dEdX)>> cut(60_GeV, true, dEdX); cut.printThresholds(); // RNGManager::getInstance().registerRandomStream("HadronicElasticModel"); @@ -134,24 +142,14 @@ int main() { TrackWriter trackWriter; output.add("tracks", trackWriter); // register TrackWriter - ShowerAxis const showerAxis{injectionPos, Vector{rootCS, 0_m, 0_m, -100_km}, env}; - BetheBlochPDG eLoss{showerAxis}; - // assemble all processes into an ordered process list auto sequence = make_sequence(pythia, decay, eLoss, cut, trackWriter, stackInspect); // define air shower object, run simulation Cascade EAS(env, tracking, sequence, output, stack); - output.startOfShower(); + output.startOfLibrary(); EAS.run(); - output.endOfShower(); - - cout << "Result: E0=" << E0 / 1_GeV << endl; - cut.showResults(); - const HEPEnergyType Efinal = - cut.getCutEnergy() + cut.getInvEnergy() + cut.getEmEnergy(); - cout << "total energy (GeV): " << Efinal / 1_GeV << endl - << "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl; - output.endOfLibrary(); + + CORSIKA_LOG_INFO("Done"); } diff --git a/examples/corsika.cpp b/examples/corsika.cpp index 6cf3051ad33a9ef51aa7bd8c91b8ed6a7e3dda81..ffda4726879bea8baeb3b1b319dc49d842c851b9 100644 --- a/examples/corsika.cpp +++ b/examples/corsika.cpp @@ -27,7 +27,9 @@ #include <corsika/framework/random/RNGManager.hpp> #include <corsika/output/OutputManager.hpp> -#include <corsika/output/NoOutput.hpp> +#include <corsika/modules/writers/SubWriter.hpp> +#include <corsika/modules/writers/EnergyLossWriter.hpp> +#include <corsika/modules/writers/LongitudinalWriter.hpp> #include <corsika/media/Environment.hpp> #include <corsika/media/FlatExponential.hpp> @@ -62,7 +64,6 @@ #include <CLI/Config.hpp> #include <iomanip> -#include <iostream> #include <limits> #include <string> @@ -155,10 +156,6 @@ int main(int argc, char** argv) { app.add_flag("--force-interaction", force_interaction, "Force the location of the first interaction.") ->group("Misc."); - app.add_option("-v,--verbosity", "Verbosity level") - ->default_str("info") - ->check(CLI::IsMember({"warn", "info", "debug", "trace"})) - ->group("Misc."); app.add_option("-v,--verbosity", "Verbosity level: warn, info, debug, trace.") ->default_val("info") ->check(CLI::IsMember({"warn", "info", "debug", "trace"})) @@ -189,8 +186,7 @@ int main(int argc, char** argv) { // gets all messed up if (app.count("--pdg") == 0) { if ((app.count("-A") == 0) || (app.count("-Z") == 0)) { - std::cerr << "If --pdg is not provided, then both -A and -Z are required." - << std::endl; + CORSIKA_LOG_ERROR("If --pdg is not provided, then both -A and -Z are required."); return 1; } } @@ -248,7 +244,7 @@ int main(int argc, char** argv) { auto const phiRad = app["--azimuth"]->as<double>() / 180. * M_PI; // convert Elab to Plab - HEPMomentumType P0 = sqrt((E0 - mass) * (E0 + mass)); + HEPMomentumType P0 = calculate_momentum(E0, mass); // convert the momentum to the zenith and azimuth angle of the primary auto const [px, py, pz] = @@ -278,10 +274,13 @@ int main(int argc, char** argv) { // create the output manager that we then register outputs with OutputManager output(app["--filename"]->as<std::string>()); - /* === START: SETUP PROCESS LIST === */ - // corsika::epos::Interaction heModel; - // corsika::qgsjetII::Interaction heModel; - // InteractionCounter heModelCounted(heModel); + // register energy losses as output + EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200}; + output.add("energyloss", dEdX); + + // create a track writer and register it with the output manager + TrackWriter tracks; + output.add("tracks", tracks); corsika::sibyll::Interaction sibyll; InteractionCounter sibyllCounted(sibyll); @@ -314,23 +313,23 @@ int main(int argc, char** argv) { // decaySibyll.printDecayConfig(); - HEPEnergyType const emcut = 1_GeV; - HEPEnergyType const hadcut = 1_GeV; - ParticleCut cut(emcut, emcut, hadcut, hadcut, true); + HEPEnergyType const emcut = 50_GeV; + HEPEnergyType const hadcut = 50_GeV; + ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, true, dEdX); corsika::proposal::Interaction emCascade(env); // NOT available for PROPOSAL due to interface trouble: - // InteractionCounter emCascadeCounted(emCascade); - // corsika::proposal::ContinuousProcess emContinuous(env); - BetheBlochPDG emContinuous(showerAxis); - - // cut.printThresholds(); + // InteractionCounter emCascadeCounted(emCascade); + // corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env); + BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX}; - LongitudinalProfile longprof(showerAxis); + LongitudinalWriter profile{showerAxis, 10_g / square(1_cm), 200}; + output.add("profile", profile); + LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile}; corsika::urqmd::UrQMD urqmd; InteractionCounter urqmdCounted(urqmd); - StackInspector<setup::Stack> stackInspect(50000, false, E0); + StackInspector<setup::Stack> stackInspect(10000, false, E0); // assemble all processes into an ordered process list struct EnergySwitch { @@ -342,15 +341,13 @@ int main(int argc, char** argv) { auto hadronSequence = make_select(EnergySwitch(63.1_GeV), urqmdCounted, heModelCounted); auto decaySequence = make_sequence(decayPythia, decaySibyll); - // track writer - TrackWriter trackWriter; - output.add("tracks", trackWriter); // register TrackWriter + TrackWriter trackWriter{tracks}; // observation plane Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.})); - ObservationPlane<setup::Tracking> observationLevel( - obsPlane, DirectionVector(rootCS, {1., 0., 0.})); - // register the observation plane with the output + ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{ + obsPlane, DirectionVector(rootCS, {1., 0., 0.})}; + // register ground particle output output.add("particles", observationLevel); // assemble the final process sequence @@ -383,14 +380,10 @@ int main(int argc, char** argv) { CORSIKA_LOG_INFO("Shower {} / {} ", i_shower, nevent); - // trigger the start of the outputs for this shower - output.startOfShower(); - // directory for outputs string const outdir(app["--filename"]->as<std::string>()); string const labHist_file = outdir + "/inthist_lab_" + to_string(i_shower) + ".npz"; string const cMSHist_file = outdir + "/inthist_cms_" + to_string(i_shower) + ".npz"; - string const longprof_file = outdir + "/longprof_" + to_string(i_shower) + ".txt"; // setup particle stack, and add primary particle stack.clear(); @@ -409,17 +402,14 @@ int main(int argc, char** argv) { // run the shower EAS.run(); - cut.showResults(); - // emContinuous.showResults(); - observationLevel.showResults(); - const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() + - cut.getEmEnergy() + // emContinuous.getEnergyLost() + - observationLevel.getEnergyGround(); - cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl - << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl; - observationLevel.reset(); - cut.reset(); - // emContinuous.reset(); + HEPEnergyType const Efinal = + dEdX.getEnergyLost() + observationLevel.getEnergyGround(); + + CORSIKA_LOG_INFO( + "total energy budget (GeV): {} (dEdX={} ground={}), " + "relative difference (%): {}", + Efinal / 1_GeV, dEdX.getEnergyLost() / 1_GeV, + observationLevel.getEnergyGround() / 1_GeV, (Efinal / E0 - 1) * 100); // auto const hists = heModelCounted.getHistogram() + urqmdCounted.getHistogram(); auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() + @@ -427,10 +417,6 @@ int main(int argc, char** argv) { save_hist(hists.labHist(), labHist_file, true); save_hist(hists.CMSHist(), cMSHist_file, true); - longprof.save(longprof_file); - - // trigger the output manager to save this shower to disk - output.endOfShower(); } // and finalize the output on disk diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp index fd6191d84a5f24495b711b9ad794745368b25967..b4a2c39d414290318a9ceed010a96cf25ddfbf02 100644 --- a/examples/em_shower.cpp +++ b/examples/em_shower.cpp @@ -22,6 +22,9 @@ #include <corsika/framework/utility/SaveBoostHistogram.hpp> #include <corsika/output/OutputManager.hpp> +#include <corsika/modules/writers/SubWriter.hpp> +#include <corsika/modules/writers/EnergyLossWriter.hpp> +#include <corsika/modules/writers/LongitudinalWriter.hpp> #include <corsika/media/Environment.hpp> #include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp> @@ -107,10 +110,7 @@ int main(int argc, char** argv) { double theta = 0.; auto const thetaRad = theta / 180. * M_PI; - auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { - return sqrt((Elab - m) * (Elab + m)); - }; - HEPMomentumType P0 = elab2plab(E0, mass); + HEPMomentumType P0 = calculate_momentum(E0, mass); auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) { return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad)); }; @@ -137,60 +137,61 @@ int main(int argc, char** argv) { beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)), plab.normalized(), injectionPos, 0_ns)); - std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02 - << std::endl; + CORSIKA_LOG_INFO("shower axis length: {} ", + (showerCore - injectionPos).getNorm() * 1.02); - OutputManager output("em_shower_outputs"); ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env, false, 1000}; + OutputManager output("em_shower_outputs"); + + EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200}; + // register energy losses as output + output.add("dEdX", dEdX); + // setup processes, decays and interactions - ParticleCut cut(60_GeV, 60_GeV, 100_PeV, 100_PeV, true); + ParticleCut<SubWriter<decltype(dEdX)>> cut(60_GeV, 60_GeV, 100_PeV, 100_PeV, true, + dEdX); corsika::proposal::Interaction emCascade(env); - corsika::proposal::ContinuousProcess emContinuous(env); + corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX); + // BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX}; // NOT possible right now, due to interface differenc in PROPOSAL // InteractionCounter emCascadeCounted(emCascade); - TrackWriter trackWriter; - output.add("tracks", trackWriter); // register TrackWriter + TrackWriter tracks; + output.add("tracks", tracks); - // long. profile; columns for photon, e+, e- still need to be added - LongitudinalProfile longprof(showerAxis); + // long. profile + LongitudinalWriter profile{showerAxis, 10_g / square(1_cm), 200}; + output.add("profile", profile); + LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile}; Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.})); - ObservationPlane<setup::Tracking> observationLevel( - obsPlane, DirectionVector(rootCS, {1., 0., 0.}), "particles.dat"); - output.add("obsplane", observationLevel); + ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{ + obsPlane, DirectionVector(rootCS, {1., 0., 0.})}; + output.add("particles", observationLevel); - auto sequence = make_sequence(emCascade, emContinuous, longprof, cut, observationLevel, - trackWriter); + auto sequence = + make_sequence(emCascade, emContinuous, longprof, cut, observationLevel, tracks); // define air shower object, run simulation setup::Tracking tracking; + + output.startOfLibrary(); Cascade EAS(env, tracking, sequence, output, stack); // to fix the point of first interaction, uncomment the following two lines: - // EAS.setNodes(); // EAS.forceInteraction(); - output.startOfShower(); EAS.run(); - output.endOfShower(); - - cut.showResults(); - emContinuous.showResults(); - observationLevel.showResults(); - const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() + - cut.getEmEnergy() + emContinuous.getEnergyLost() + - observationLevel.getEnergyGround(); - cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl - << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl; - observationLevel.reset(); - cut.reset(); - emContinuous.reset(); - - longprof.save("longprof_emShower.txt"); + + HEPEnergyType const Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround(); + + CORSIKA_LOG_INFO( + "total energy budget (GeV): {}, " + "relative difference (%): {}", + Efinal / 1_GeV, (Efinal / E0 - 1) * 100); output.endOfLibrary(); } diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp index 2de965757d4e0c3c0e0df83af5ab2cf47aa9055b..4e695063468ad88ca2159c0642f57eb54117eced 100644 --- a/examples/hybrid_MC.cpp +++ b/examples/hybrid_MC.cpp @@ -27,6 +27,9 @@ #include <corsika/framework/random/RNGManager.hpp> #include <corsika/output/OutputManager.hpp> +#include <corsika/modules/writers/SubWriter.hpp> +#include <corsika/modules/writers/EnergyLossWriter.hpp> +#include <corsika/modules/writers/LongitudinalWriter.hpp> #include <corsika/media/Environment.hpp> #include <corsika/media/FlatExponential.hpp> @@ -40,7 +43,7 @@ #include <corsika/modules/BetheBlochPDG.hpp> #include <corsika/modules/LongitudinalProfile.hpp> #include <corsika/modules/ObservationPlane.hpp> -#include <corsika/modules/OnShellCheck.hpp> +#include <corsika/modules/TrackWriter.hpp> #include <corsika/modules/ParticleCut.hpp> #include <corsika/modules/Pythia8.hpp> #include <corsika/modules/Sibyll.hpp> @@ -69,6 +72,11 @@ using namespace corsika; using namespace std; +/** + * Random number stream initialization + * + * @param seed + */ void registerRandomStreams(uint64_t seed) { RNGManager<>::getInstance().registerRandomStream("cascade"); RNGManager<>::getInstance().registerRandomStream("qgsjet"); @@ -81,17 +89,34 @@ void registerRandomStreams(uint64_t seed) { if (seed == 0) { std::random_device rd; seed = rd(); - cout << "new random seed (auto) " << seed << endl; + CORSIKA_LOG_INFO("new random seed (auto) {}", seed); } RNGManager<>::getInstance().setSeed(seed); } +/** + * New (for demonstration) ContinuousProcess which will check if a particles has traversed + * below the observation level. + */ class TrackCheck : public ContinuousProcess<TrackCheck> { public: + /** + * Construct a new Track Check object. + * + * @param plane -- the actual observation level + */ TrackCheck(Plane const& plane) : plane_(plane) {} + /** + * The doContinous method to check a particular particle. + * + * @tparam TParticle + * @tparam TTrack + * @param particle + * @return ProcessReturn + */ template <typename TParticle, typename TTrack> ProcessReturn doContinuous(TParticle const& particle, TTrack const&, bool const) { auto const delta = particle.getPosition() - plane_.getCenter(); @@ -105,6 +130,13 @@ public: return ProcessReturn::Ok; } + /** + * No limit on tracking step length imposed here, of course. + * + * @tparam TParticle + * @tparam TTrack + * @return LengthType + */ template <typename TParticle, typename TTrack> LengthType getMaxStepLength(TParticle const&, TTrack const&) const { return std::numeric_limits<double>::infinity() * 1_m; @@ -114,6 +146,9 @@ private: Plane plane_; }; +/** + * Selection of environment interface implementation: + */ template <typename T> using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>; @@ -124,8 +159,10 @@ int main(int argc, char** argv) { CORSIKA_LOG_INFO("hybrid_MC"); if (argc < 4) { - std::cerr << "usage: hybrid_MC <A> <Z> <energy/GeV> [seed]" << std::endl; - std::cerr << " if no seed is given, a random seed is chosen" << std::endl; + CORSIKA_LOG_ERROR( + "\n" + "usage: hybrid_MC <A> <Z> <energy/GeV> [seed] \n" + " if no seed is given, a random seed is chosen"); return 1; } feenableexcept(FE_INVALID); @@ -157,20 +194,19 @@ int main(int argc, char** argv) { double theta = 0.; auto const thetaRad = theta / 180. * M_PI; - auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { - return sqrt((Elab - m) * (Elab + m)); - }; - HEPMomentumType P0 = elab2plab(E0, mass); + HEPMomentumType P0 = calculate_momentum(E0, mass); auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) { return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad)); }; auto const [px, py, pz] = momentumComponents(thetaRad, P0); auto plab = MomentumVector(rootCS, {px, py, pz}); - cout << "input particle: " << beamCode << endl; - cout << "input angles: theta=" << theta << endl; - cout << "input momentum: " << plab.getComponents() / 1_GeV - << ", norm = " << plab.getNorm() << endl; + CORSIKA_LOG_INFO( + "input particle: {}, " + "input angles: theta={}, " + "input momentum: {} GeV, " + ", norm={}", + beamCode, theta, plab.getComponents() / 1_GeV, plab.getNorm()); auto const observationHeight = 0_km + constants::EarthRadius::Mean; auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean; @@ -182,18 +218,18 @@ int main(int argc, char** argv) { showerCore + Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t; - std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl; + CORSIKA_LOG_INFO("point of injection: {} ", injectionPos.getCoordinates()); stack.addParticle(std::make_tuple( Code::Proton, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)), plab.normalized(), injectionPos, 0_ns)); - std::cout << "shower axis length: " << (showerCore - injectionPos).getNorm() * 1.02 - << std::endl; + CORSIKA_LOG_INFO("shower axis length: {} m", + (showerCore - injectionPos).getNorm() * 1.02); OutputManager output("hybrid_MC_outputs"); - ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env, true, - 1000}; + ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env, + false, 1000}; // setup processes, decays and interactions @@ -228,15 +264,24 @@ int main(int argc, char** argv) { decaySibyll.printDecayConfig(); - ParticleCut cut(3_GeV, false, true); - BetheBlochPDG eLoss(showerAxis); + // register energy losses as output + EnergyLossWriter dEdX{showerAxis, 10_g / square(1_cm), 200}; + output.add("energyloss", dEdX); + + // create a track writer and register it with the output manager + TrackWriter<TrackWriterParquet> tracks; + output.add("tracks", tracks); - CONEXhybrid conex_model(center, showerAxis, t, injectionHeight, E0, - get_PDG(Code::Proton)); + ParticleCut<SubWriter<decltype(dEdX)>> cut(3_GeV, false, dEdX); + BetheBlochPDG<SubWriter<decltype(dEdX)>> eLoss(dEdX); - OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false); + LongitudinalWriter profile{showerAxis, 10_g / square(1_cm), 200}; + output.add("profile", profile); + LongitudinalProfile<SubWriter<decltype(profile)>> longprof{profile}; - LongitudinalProfile longprof(showerAxis); + CONEXhybrid // SubWriter<decltype(dEdX>, SubWriter<decltype(profile)>> + conex_model(center, showerAxis, t, injectionHeight, E0, get_PDG(Code::Proton), dEdX, + profile); Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.})); ObservationPlane<setup::Tracking> observationLevel( @@ -246,6 +291,8 @@ int main(int argc, char** argv) { corsika::urqmd::UrQMD urqmd_model; InteractionCounter urqmdCounted{urqmd_model}; + TrackCheck trackCheck(obsPlane); + // assemble all processes into an ordered process list struct EnergySwitch { HEPEnergyType cutE_; @@ -258,8 +305,8 @@ int main(int argc, char** argv) { auto hadronSequence = make_select(EnergySwitch(55_GeV), urqmdCounted, make_sequence(sibyllNucCounted, sibyllCounted)); auto decaySequence = make_sequence(decayPythia, decaySibyll); - auto sequence = make_sequence(hadronSequence, reset_particle_mass, decaySequence, eLoss, - cut, conex_model, longprof, observationLevel); + auto sequence = make_sequence(hadronSequence, decaySequence, eLoss, cut, conex_model, + longprof, observationLevel, trackCheck); // define air shower object, run simulation setup::Tracking tracking; @@ -269,21 +316,15 @@ int main(int argc, char** argv) { // EAS.SetNodes(); // EAS.forceInteraction(); - output.startOfShower(); + output.startOfLibrary(); EAS.run(); - output.endOfShower(); - - cut.showResults(); - eLoss.showResults(); - observationLevel.showResults(); - const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() + - cut.getEmEnergy() + eLoss.getEnergyLost() + - observationLevel.getEnergyGround(); - cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl - << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl; - observationLevel.reset(); - cut.reset(); - eLoss.reset(); + output.endOfLibrary(); + + const HEPEnergyType Efinal = dEdX.getEnergyLost() + observationLevel.getEnergyGround(); + CORSIKA_LOG_INFO( + "total cut energy (GeV): {}, " + "relative difference (%): {}", + Efinal / 1_GeV, (Efinal / E0 - 1) * 100); auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() + urqmdCounted.getHistogram(); @@ -291,10 +332,5 @@ int main(int argc, char** argv) { save_hist(hists.labHist(), "inthist_lab_hybrid.npz", true); save_hist(hists.CMSHist(), "inthist_cms_hybrid.npz", true); - longprof.save("longprof.txt"); - - std::ofstream finish("finished"); - finish << "run completed without error" << std::endl; - - output.endOfLibrary(); + CORSIKA_LOG_INFO("done"); } diff --git a/examples/mars.cpp b/examples/mars.cpp index 12a4039bb89c983b50408a52369d4b0d5ef4af52..84cf468233b48f08cc67f6d481d627dcbda1101b 100644 --- a/examples/mars.cpp +++ b/examples/mars.cpp @@ -29,7 +29,9 @@ #include <corsika/framework/random/RNGManager.hpp> #include <corsika/output/OutputManager.hpp> -#include <corsika/output/NoOutput.hpp> +#include <corsika/modules/writers/SubWriter.hpp> +#include <corsika/modules/writers/EnergyLossWriter.hpp> +#include <corsika/modules/writers/LongitudinalWriter.hpp> #include <corsika/media/Environment.hpp> #include <corsika/media/FlatExponential.hpp> @@ -45,7 +47,6 @@ #include <corsika/modules/BetheBlochPDG.hpp> #include <corsika/modules/LongitudinalProfile.hpp> #include <corsika/modules/ObservationPlane.hpp> -#include <corsika/modules/OnShellCheck.hpp> #include <corsika/modules/StackInspector.hpp> #include <corsika/modules/TrackWriter.hpp> #include <corsika/modules/ParticleCut.hpp> @@ -208,8 +209,7 @@ int main(int argc, char** argv) { // gets all messed up if (app.count("--pdg") == 0) { if ((app.count("-A") == 0) || (app.count("-Z") == 0)) { - std::cerr << "If --pdg is not provided, then both -A and -Z are required." - << std::endl; + CORSIKA_LOG_ERROR("If --pdg is not provided, then both -A and -Z are required."); return 1; } } @@ -305,12 +305,20 @@ int main(int argc, char** argv) { // we make the axis much longer than the inj-core distance since the // profile will go beyond the core, depending on zenith angle - ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env}; /* === END: CONSTRUCT GEOMETRY === */ // create the output manager that we then register outputs with OutputManager output(app["--filename"]->as<std::string>()); + ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env}; + + EnergyLossWriter dEdX{showerAxis}; + output.add("energyloss", dEdX); + + HEPEnergyType const emcut = 1_GeV; + HEPEnergyType const hadcut = 1_GeV; + ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, true, dEdX); + /* === START: SETUP PROCESS LIST === */ corsika::sibyll::Interaction sibyll; InteractionCounter sibyllCounted(sibyll); @@ -343,16 +351,16 @@ int main(int argc, char** argv) { // decaySibyll.printDecayConfig(); - HEPEnergyType const emcut = 1_GeV; - HEPEnergyType const hadcut = 1_GeV; - ParticleCut cut(emcut, emcut, hadcut, hadcut, true); corsika::proposal::Interaction emCascade(env); // NOT possible right now, due to interface difference for PROPOSAL: // InteractionCounter emCascadeCounted(emCascade); - // corsika::proposal::ContinuousProcess emContinuous(env); - BetheBlochPDG emContinuous(showerAxis); + // corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> + // emContinuous(env,dEdX); + BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX}; - LongitudinalProfile longprof{showerAxis, 1_g / square(1_cm)}; + LongitudinalWriter longprof{showerAxis}; + output.add("profile", longprof); + LongitudinalProfile<SubWriter<decltype(longprof)>> profile{longprof}; corsika::urqmd::UrQMD urqmd; InteractionCounter urqmdCounted{urqmd}; @@ -382,7 +390,7 @@ int main(int argc, char** argv) { // assemble the final process sequence auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence, emCascade, emContinuous, - cut, trackWriter, observationLevel, longprof); + cut, trackWriter, observationLevel, profile); /* === END: SETUP PROCESS LIST === */ // create the cascade object using the default stack and tracking implementation @@ -415,7 +423,6 @@ int main(int argc, char** argv) { string const outdir(app["--filename"]->as<std::string>()); string const labHist_file = outdir + "/inthist_lab_" + to_string(i_shower) + ".npz"; string const cMSHist_file = outdir + "/inthist_cms_" + to_string(i_shower) + ".npz"; - string const longprof_file = outdir + "/longprof_" + to_string(i_shower) + ".txt"; // setup particle stack, and add primary particle stack.clear(); @@ -428,24 +435,19 @@ int main(int argc, char** argv) { // run the shower EAS.run(); - cut.showResults(); - // emContinuous.showResults(); - observationLevel.showResults(); - const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() + - cut.getEmEnergy() + // emContinuous.getEnergyLost() + - observationLevel.getEnergyGround(); - cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl - << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl; - observationLevel.reset(); - cut.reset(); - // emContinuous.reset(); + HEPEnergyType const Efinal = + dEdX.getEnergyLost() + observationLevel.getEnergyGround(); + + CORSIKA_LOG_INFO( + "total energy budget (GeV): {}, " + "relative difference (%): {}", + Efinal / 1_GeV, (Efinal / E0 - 1) * 100); auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() + urqmdCounted.getHistogram(); save_hist(hists.labHist(), labHist_file, true); save_hist(hists.CMSHist(), cMSHist_file, true); - longprof.save(longprof_file); // trigger the output manager to save this shower to disk output.endOfShower(); diff --git a/examples/stopping_power.cpp b/examples/stopping_power.cpp index 8ca094f6e5c9f26b436236e21f71c12a7f5c2b7a..8d0641960940c34e00d110780735b4d06ecdfb9a 100644 --- a/examples/stopping_power.cpp +++ b/examples/stopping_power.cpp @@ -48,9 +48,7 @@ int main() { rootCS, 0_m, 0_m, 112.8_km); // this is the CORSIKA 7 start of atmosphere/universe - ShowerAxis showerAxis{injectionPos, Vector<length_d>{rootCS, 0_m, 0_m, 1_m}, env, false, - 100}; - BetheBlochPDG eLoss{showerAxis}; + BetheBlochPDG eLoss; setup::Stack stack; @@ -64,10 +62,7 @@ int main() { double theta = 0.; double phi = 0.; - auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { - return sqrt((Elab - m) * (Elab + m)); - }; - HEPMomentumType P0 = elab2plab(E0, mass); + HEPMomentumType P0 = calculate_momentum(E0, mass); auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) { return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi), -ptot * cos(theta)); diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp index f38b1b2904d5eb777e0449571e892af58d82ccea..ba69ffe2814e599a3f7bf6eb582a864beffc62c5 100644 --- a/examples/vertical_EAS.cpp +++ b/examples/vertical_EAS.cpp @@ -29,7 +29,9 @@ #include <corsika/framework/random/RNGManager.hpp> #include <corsika/output/OutputManager.hpp> -#include <corsika/output/NoOutput.hpp> +#include <corsika/modules/writers/SubWriter.hpp> +#include <corsika/modules/writers/EnergyLossWriter.hpp> +#include <corsika/modules/writers/LongitudinalWriter.hpp> #include <corsika/media/Environment.hpp> #include <corsika/media/FlatExponential.hpp> @@ -152,10 +154,7 @@ int main(int argc, char** argv) { auto const thetaRad = theta / 180. * constants::pi; auto const phiRad = phi / 180. * constants::pi; - auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) { - return sqrt((Elab - m) * (Elab + m)); - }; - HEPMomentumType P0 = elab2plab(E0, mass); + HEPMomentumType P0 = calculate_momentum(E0, mass); auto momentumComponents = [](double theta, double phi, HEPMomentumType ptot) { return std::make_tuple(ptot * sin(theta) * cos(phi), ptot * sin(theta) * sin(phi), -ptot * cos(theta)); @@ -193,6 +192,31 @@ int main(int argc, char** argv) { // create the output manager that we then register outputs with OutputManager output("vertical_EAS_outputs"); + // EnergyLossWriterParquet dEdX_output{showerAxis, 10_g / square(1_cm), 200}; + // register energy losses as output + // output.add("dEdX", dEdX_output); + // register profile output + + // construct the overall energy loss writer and register it + EnergyLossWriter dEdX{showerAxis}; + output.add("energyloss", dEdX); + + // construct the continuous energy loss model + BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX}; + + // construct a particle cut + ParticleCut<SubWriter<decltype(dEdX)>> cut{E0, E0, 60_GeV, 60_GeV, true, dEdX}; + + // setup longitudinal profile + LongitudinalWriter longProf{showerAxis}; + output.add("profile", longProf); + + LongitudinalProfile<SubWriter<decltype(longProf)>> profile{longProf}; + + // create a track writer and register it with the output manager + TrackWriter<TrackWriterParquet> trackWriter; + output.add("tracks", trackWriter); + // setup processes, decays and interactions corsika::sibyll::Interaction sibyll; @@ -226,8 +250,6 @@ int main(int argc, char** argv) { decaySibyll.printDecayConfig(); - ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true}; - corsika::urqmd::UrQMD urqmd; InteractionCounter urqmdCounted{urqmd}; StackInspector<setup::Stack> stackInspect(50000, false, E0); @@ -246,7 +268,6 @@ int main(int argc, char** argv) { // directory for outputs string const labHist_file = "inthist_lab_verticalEAS.npz"; string const cMSHist_file = "inthist_cms_verticalEAS.npz"; - string const longprof_file = "longprof_verticalEAS.txt"; // setup particle stack, and add primary particle setup::Stack stack; @@ -256,21 +277,16 @@ int main(int argc, char** argv) { beamCode, calculate_kinetic_energy(plab.getNorm(), get_mass(beamCode)), plab.normalized(), injectionPos, 0_ns)); - BetheBlochPDG emContinuous(showerAxis); - - TrackWriter trackWriter; - output.add("tracks", trackWriter); // register TrackWriter - - LongitudinalProfile longprof{showerAxis}; - Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.})); - ObservationPlane<setup::Tracking, NoOutput> observationLevel( - obsPlane, DirectionVector(rootCS, {1., 0., 0.})); - // register the observation plane with the output + ObservationPlane<setup::Tracking, ParticleWriterParquet> observationLevel{ + obsPlane, DirectionVector(rootCS, {1., 0., 0.})}; output.add("particles", observationLevel); + // auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence, + // emContinuous, + // cut, trackWriter, observationLevel, longprof); auto sequence = make_sequence(stackInspect, hadronSequence, decaySequence, emContinuous, - cut, trackWriter, observationLevel, longprof); + cut, trackWriter, observationLevel, profile); // define air shower object, run simulation setup::Tracking tracking; @@ -279,24 +295,11 @@ int main(int argc, char** argv) { EAS.run(); output.endOfShower(); - cut.showResults(); - // emContinuous.showResults(); - observationLevel.showResults(); - const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() + - cut.getEmEnergy() + // emContinuous.getEnergyLost() + - observationLevel.getEnergyGround(); - cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl - << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl; - observationLevel.reset(); - cut.reset(); - // emContinuous.reset(); - auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() + urqmdCounted.getHistogram(); save_hist(hists.labHist(), labHist_file, true); save_hist(hists.CMSHist(), cMSHist_file, true); - longprof.save(longprof_file); output.endOfLibrary(); } diff --git a/python/corsika/io/outputs/__init__.py b/python/corsika/io/outputs/__init__.py index fa38acdde3c2dc8dd3cc351c16d27159cecd1efb..d39d9a8b4bd68465ab1d799f10c15bb20032cbc9 100644 --- a/python/corsika/io/outputs/__init__.py +++ b/python/corsika/io/outputs/__init__.py @@ -9,6 +9,16 @@ from .observation_plane import ObservationPlane from .track_writer import TrackWriter +from .longitudinal_profile import LongitudinalProfile +from .bethe_bloch import BetheBlochPDG +from .particle_cut import ParticleCut from .output import Output -__all__ = ["Output", "ObservationPlane", "TrackWriter"] +__all__ = [ + "Output", + "ObservationPlane", + "TrackWriter", + "LongitudinalProfile", + "BetheBlochPDG", + "ParticleCut", +] diff --git a/python/corsika/io/outputs/bethe_bloch.py b/python/corsika/io/outputs/bethe_bloch.py new file mode 100644 index 0000000000000000000000000000000000000000..b12d29036de5553d139c9cee0f7c8120f4fd4f7b --- /dev/null +++ b/python/corsika/io/outputs/bethe_bloch.py @@ -0,0 +1,87 @@ +""" + Read data written by BetheBlochPDG. + + (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" +import logging +import os.path as op +from typing import Any + +import pyarrow.parquet as pq + +from .output import Output + + +class BetheBlochPDG(Output): + """ + Read particle data from an BetheBlochPDG. + """ + + def __init__(self, path: str): + """ + Load the particle data into a parquet table. + + Parameters + ---------- + path: str + The path to the directory containing this output. + """ + super().__init__(path) + + # try and load our data + try: + self.__data = pq.read_table(op.join(path, "energyloss.parquet")) + except Exception as e: + logging.getLogger("corsika").warn( + f"An error occured loading a BetheBlochPDG: {e}" + ) + + def is_good(self) -> bool: + """ + Returns true if this output has been read successfully + and has the correct files/state/etc. + + Returns + ------- + bool: + True if this is a good output. + """ + return self.__data is not None + + def astype(self, dtype: str = "pandas", **kwargs: Any) -> Any: + """ + Load the particle data from this bethe bloch instance. + + All additional keyword arguments are passed to `parquet.read_table` + + Parameters + ---------- + dtype: str + The data format to return the data in (i.e. numpy, pandas, etc.) + + Returns + ------- + Any: + The return type of this method is determined by `dtype`. + """ + if dtype == "arrow": + return self.__data + elif dtype == "pandas": + return self.__data.to_pandas() + else: + raise ValueError( + ( + f"Unknown format '{dtype}' for BetheBlochPDG. " + "We currently only support ['arrow', 'pandas']." + ) + ) + + def __repr__(self) -> str: + """ + Return a string representation of this class. + """ + return f"BetheBlochPDG('{self.config['name']}')" diff --git a/python/corsika/io/outputs/longitudinal_profile.py b/python/corsika/io/outputs/longitudinal_profile.py new file mode 100644 index 0000000000000000000000000000000000000000..00000567d827a12e94c198a0ef1102b1c010a74e --- /dev/null +++ b/python/corsika/io/outputs/longitudinal_profile.py @@ -0,0 +1,87 @@ +""" + Read data written by LongitudinalProfile. + + (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" +import logging +import os.path as op +from typing import Any + +import pyarrow.parquet as pq + +from .output import Output + + +class LongitudinalProfile(Output): + """ + Read particle data from an LongitudinalProfile. + """ + + def __init__(self, path: str): + """ + Load the particle data into a parquet table. + + Parameters + ---------- + path: str + The path to the directory containing this output. + """ + super().__init__(path) + + # try and load our data + try: + self.__data = pq.read_table(op.join(path, "profile.parquet")) + except Exception as e: + logging.getLogger("corsika").warn( + f"An error occured loading a LongitudinalProfile: {e}" + ) + + def is_good(self) -> bool: + """ + Returns true if this output has been read successfully + and has the correct files/state/etc. + + Returns + ------- + bool: + True if this is a good output. + """ + return self.__data is not None + + def astype(self, dtype: str = "pandas", **kwargs: Any) -> Any: + """ + Load the particle data from this longitudinal profile. + + All additional keyword arguments are passed to `parquet.read_table` + + Parameters + ---------- + dtype: str + The data format to return the data in (i.e. numpy, pandas, etc.) + + Returns + ------- + Any: + The return type of this method is determined by `dtype`. + """ + if dtype == "arrow": + return self.__data + elif dtype == "pandas": + return self.__data.to_pandas() + else: + raise ValueError( + ( + f"Unknown format '{dtype}' for LongitudinalProfile. " + "We currently only support ['arrow', 'pandas']." + ) + ) + + def __repr__(self) -> str: + """ + Return a string representation of this class. + """ + return f"LongitudinalProfile('{self.config['name']}')" diff --git a/python/corsika/io/outputs/output.py b/python/corsika/io/outputs/output.py index 83a4c7928b786b9bad4301bec39b7e03aaf567f8..fed4368aac5d1cb8be7c1cbe6a955f1f799df4df 100644 --- a/python/corsika/io/outputs/output.py +++ b/python/corsika/io/outputs/output.py @@ -145,6 +145,9 @@ class Output(ABC): """ Load the top-level summary from a given library path. + If there is not a summary file for this output, then an + empty dictionary is returned. + Parameters ---------- path: str @@ -154,12 +157,11 @@ class Output(ABC): ------- dict: The summary as a python dictionary. - - Raises - ------ - FileNotFoundError - If the summary file cannot be found - """ - with open(op.join(path, "summary.yaml"), "r") as f: - return yaml.load(f, Loader=yaml.Loader) + + # if the summary file doesn't exist, we just an empty dict + if not op.exists(op.join(path, "summary.yaml")): + return {} + else: + with open(op.join(path, "summary.yaml"), "r") as f: + return yaml.load(f, Loader=yaml.Loader) diff --git a/python/corsika/io/outputs/particle_cut.py b/python/corsika/io/outputs/particle_cut.py new file mode 100644 index 0000000000000000000000000000000000000000..a55b5fe718cf788fdba051dfa891e77b2d75ede7 --- /dev/null +++ b/python/corsika/io/outputs/particle_cut.py @@ -0,0 +1,87 @@ +""" + Read data written by ParticleCut. + + (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + + This software is distributed under the terms of the GNU General Public + Licence version 3 (GPL Version 3). See file LICENSE for a full version of + the license. +""" +import logging +import os.path as op +from typing import Any + +import pyarrow.parquet as pq + +from .output import Output + + +class ParticleCut(Output): + """ + Read particle data from an ParticleCut. + """ + + def __init__(self, path: str): + """ + Load the particle data into a parquet table. + + Parameters + ---------- + path: str + The path to the directory containing this output. + """ + super().__init__(path) + + # try and load our data + try: + self.__data = pq.read_table(op.join(path, "energyloss.parquet")) + except Exception as e: + logging.getLogger("corsika").warn( + f"An error occured loading a ParticleCut: {e}" + ) + + def is_good(self) -> bool: + """ + Returns true if this output has been read successfully + and has the correct files/state/etc. + + Returns + ------- + bool: + True if this is a good output. + """ + return self.__data is not None + + def astype(self, dtype: str = "pandas", **kwargs: Any) -> Any: + """ + Load the particle data from this particle cut instance. + + All additional keyword arguments are passed to `parquet.read_table` + + Parameters + ---------- + dtype: str + The data format to return the data in (i.e. numpy, pandas, etc.) + + Returns + ------- + Any: + The return type of this method is determined by `dtype`. + """ + if dtype == "arrow": + return self.__data + elif dtype == "pandas": + return self.__data.to_pandas() + else: + raise ValueError( + ( + f"Unknown format '{dtype}' for ParticleCut. " + "We currently only support ['arrow', 'pandas']." + ) + ) + + def __repr__(self) -> str: + """ + Return a string representation of this class. + """ + return f"ParticleCut('{self.config['name']}')" diff --git a/tests/framework/CMakeLists.txt b/tests/framework/CMakeLists.txt index 7d67ad5aafd1f64314fa67ff4d9ed95b34ba442e..caaec9e4cc04763a30ceec719107e57f9b91a578 100644 --- a/tests/framework/CMakeLists.txt +++ b/tests/framework/CMakeLists.txt @@ -22,6 +22,7 @@ set (test_framework_sources testInteractionCounter.cpp testInteractionLengthModifier.cpp testSolver.cpp + testEnergyMomentum.cpp ) CORSIKA_ADD_TEST (testFramework SOURCES ${test_framework_sources}) diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp index 624dc041575d1304c4e7bdcb012810e7ca3cc478..d6bc33f5929ffd6c71eed9636d9166135ad18037 100644 --- a/tests/framework/testCascade.cpp +++ b/tests/framework/testCascade.cpp @@ -72,6 +72,7 @@ class DummyTracking { public: template <typename TParticle> auto getTrack(TParticle const& particle) { + calls_++; VelocityVector const initialVelocity = particle.getMomentum() / particle.getEnergy() * constants::c; Line const theLine = Line(particle.getPosition(), initialVelocity); @@ -81,8 +82,10 @@ public: // trajectory: just go ahead forever particle.getNode()); // next volume node } + int getCalls() const { return calls_; } static std::string getName() { return "DummyTracking"; } static std::string getVersion() { return "1.0.0"; } + int calls_ = 0; }; class ProcessSplit : public InteractionProcess<ProcessSplit> { @@ -174,16 +177,22 @@ TEST_CASE("Cascade", "[Cascade]") { SECTION("full cascade") { EAS.run(); + CHECK(tracking.getCalls() == 2047); CHECK(cut.getCount() == 2048); CHECK(cut.getCalls() == 2047); // final particle is still on stack and not yet deleted CHECK(split.getCalls() == 2047); } SECTION("forced interaction") { - EAS.setNodes(); + CHECK(tracking.getCalls() == 0); EAS.forceInteraction(); - CHECK(stack.getEntries() == 2); - CHECK(stack.getSize() == 3); - CHECK(split.getCalls() == 1); + CHECK(tracking.getCalls() == 0); + CHECK(stack.getEntries() == 1); + CHECK(stack.getSize() == 1); + EAS.run(); + CHECK(tracking.getCalls() == 2046); // one LESS than without forceInteraction + CHECK(stack.getEntries() == 0); + CHECK(stack.getSize() == 13); + CHECK(split.getCalls() == 2047); } } diff --git a/tests/framework/testEnergyMomentum.cpp b/tests/framework/testEnergyMomentum.cpp new file mode 100644 index 0000000000000000000000000000000000000000..efef1a9f860566f361a664e5bfd21114b625e730 --- /dev/null +++ b/tests/framework/testEnergyMomentum.cpp @@ -0,0 +1,36 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <catch2/catch.hpp> + +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/framework/core/EnergyMomentumOperations.hpp> + +using namespace corsika; + +TEST_CASE("EnergyMomentumOperations") { + + logging::set_level(logging::level::info); + + HEPMomentumType const p0 = 30_GeV; + HEPMassType const m0 = 40_GeV; + HEPEnergyType const e0 = 50_GeV; + + CHECK(calculate_mass_sqr(e0, p0) / (m0 * m0) == Approx(1)); + CHECK(calculate_mass(e0, p0) / m0 == Approx(1)); + + CHECK(calculate_momentum_sqr(e0, m0) / (p0 * p0) == Approx(1)); + CHECK(calculate_momentum(e0, m0) / p0 == Approx(1)); + + CHECK(calculate_total_energy_sqr(p0, m0) / (e0 * e0) == Approx(1)); + CHECK(calculate_total_energy(p0, m0) / e0 == Approx(1)); + + CHECK(calculate_kinetic_energy(p0, m0) / (e0 - m0) == Approx(1)); + CHECK(calculate_lab_energy((100_GeV * 100_GeV), 1_GeV, 100_MeV) / 49999.99_GeV == + Approx(1).margin(0.1)); +} diff --git a/tests/framework/testLogging.cpp b/tests/framework/testLogging.cpp index bf92cf8df00e8f12d64c0d5779f59b7976c4e09d..5b4fc2233fbb5a480d25e8ab12145b6f5887ab20 100644 --- a/tests/framework/testLogging.cpp +++ b/tests/framework/testLogging.cpp @@ -107,7 +107,7 @@ TEST_CASE("Logging", "[Logging]") { // these print with the "loggerE" logger CORSIKA_LOGGER_INFO(logger, "(8) test macro style logging"); - CORSIKA_LOGGER_WARN(logger, "(8) test macro style logging"); + CORSIKA_LOGGER_WARN(logger, "(8) test macro {} logging", "style"); // reset the logging pattern logging::reset_pattern(logger); diff --git a/tests/framework/testParticles.cpp b/tests/framework/testParticles.cpp index 3ec511b627161dec1a09389f8fdad962e3bb13b6..eaa22922212be3badea612ecb262e13249025eb1 100644 --- a/tests/framework/testParticles.cpp +++ b/tests/framework/testParticles.cpp @@ -41,6 +41,9 @@ TEST_CASE("ParticleProperties", "[Particles]") { CHECK(Positron::charge / constants::e == Approx(+1)); CHECK(get_charge(Positron::anti_code) / constants::e == Approx(-1)); CHECK(Photon::charge / constants::e == 0.); + CHECK_FALSE(is_charged(Code::Photon)); + CHECK(is_charged(Code::Iron)); + CHECK(is_charged(Code::PiPlus)); } SECTION("Names") { @@ -87,11 +90,11 @@ TEST_CASE("ParticleProperties", "[Particles]") { SECTION("Energy threshold") { //! by default energy thresholds are set to zero - CHECK(calculate_kinetic_energy_threshold(Electron::code) == 0_GeV); + CHECK(get_kinetic_energy_threshold(Electron::code) == 0_GeV); set_kinetic_energy_threshold(Electron::code, 10_GeV); - CHECK_FALSE(calculate_kinetic_energy_threshold(Code::Electron) == 1_GeV); - CHECK(calculate_kinetic_energy_threshold(Code::Electron) == 10_GeV); + CHECK_FALSE(get_kinetic_energy_threshold(Code::Electron) == 1_GeV); + CHECK(get_kinetic_energy_threshold(Code::Electron) == 10_GeV); } SECTION("Particle groups: electromagnetic") { diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt index a4e1c1b63c1b86df66e347d3ab28ca294985e97e..00f28c41cf60e57530d63d248b7208c21c3859a6 100644 --- a/tests/modules/CMakeLists.txt +++ b/tests/modules/CMakeLists.txt @@ -8,7 +8,6 @@ set (test_modules_sources testPythia8.cpp testUrQMD.cpp testCONEX.cpp - ## testOnShellCheck.cpp testParticleCut.cpp testSibyll.cpp testEpos.cpp diff --git a/tests/modules/conex_fit_REF.txt b/tests/modules/conex_fit_REF.txt deleted file mode 100644 index b9d071745e5daa1b371971d0e698a90477d9db1c..0000000000000000000000000000000000000000 --- a/tests/modules/conex_fit_REF.txt +++ /dev/null @@ -1,13 +0,0 @@ -15 # log10(eprima/eV) -60 # theta -1e+30 # X1 (first interaction) -1156.19 # Nmax -471.873 # X0 -48.4425 # P1 --0.00311659 # P2 -2.89514e-06 # P3 -0.000124814 # chi^2 / sqrt(Nmax) -833.734 # Xmax --0 # phi --1 # inelasticity 1st int. -100000 # ??? diff --git a/tests/modules/conex_output_REF.txt b/tests/modules/conex_output_REF.txt deleted file mode 100644 index 191f7fe5d90deccaf36dee714ef9176b46b38962..0000000000000000000000000000000000000000 --- a/tests/modules/conex_output_REF.txt +++ /dev/null @@ -1,208 +0,0 @@ -# X N dEdX Mu dMu Photon El Had - 0.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 10.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 20.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 30.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 40.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 50.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 60.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 70.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 80.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 90.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 100.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 110.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 120.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 130.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 140.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 150.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 160.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 170.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 180.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 190.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 200.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 210.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 220.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 230.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 240.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 250.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 260.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 270.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 280.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 290.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 300.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 310.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 320.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 330.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 340.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 350.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 360.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 370.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 380.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 390.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 400.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 410.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 420.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 430.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 440.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 450.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 460.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 470.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 480.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 490.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 500.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 510.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 520.00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 530.00 0.00000e+00 5.85428e-04 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00 - 540.00 2.00000e+00 7.23527e-03 0.00000e+00 0.00000e+00 1.00000e+00 1.00000e+00 0.00000e+00 - 550.00 3.81570e+00 1.48925e-02 4.38481e-07 4.38546e-07 1.00348e+01 2.09494e+00 2.25571e-03 - 560.00 7.86774e+00 2.91240e-02 2.46741e-03 2.56643e-03 2.24061e+01 4.37781e+00 7.06821e-03 - 570.00 1.46622e+01 5.15114e-02 7.23749e-03 5.02553e-03 4.18262e+01 8.22011e+00 1.32594e-02 - 580.00 2.49442e+01 8.40423e-02 1.43304e-02 7.56723e-03 7.16470e+01 1.40716e+01 2.01365e-02 - 590.00 3.95547e+01 1.28990e-01 2.35811e-02 1.00024e-02 1.15230e+02 2.24381e+01 2.73808e-02 - 600.00 5.92836e+01 1.88109e-01 3.47690e-02 1.22697e-02 1.75963e+02 3.38002e+01 3.48011e-02 - 610.00 8.47925e+01 2.62836e-01 4.76464e-02 1.43355e-02 2.57126e+02 4.85679e+01 4.22585e-02 - 620.00 1.16562e+02 3.53494e-01 6.19492e-02 1.61758e-02 3.61702e+02 6.70496e+01 4.96374e-02 - 630.00 1.54849e+02 4.58394e-01 7.74028e-02 1.77721e-02 4.92175e+02 8.94260e+01 5.68348e-02 - 640.00 1.99661e+02 5.79214e-01 9.37289e-02 1.91123e-02 6.50346e+02 1.15733e+02 6.37563e-02 - 650.00 2.50742e+02 7.15031e-01 1.10652e-01 2.01905e-02 8.37183e+02 1.45850e+02 7.03164e-02 - 660.00 3.07580e+02 8.64314e-01 1.27904e-01 2.10073e-02 1.05271e+03 1.79506e+02 7.64387e-02 - 670.00 3.69423e+02 1.02499e+00 1.45234e-01 2.15693e-02 1.29596e+03 2.16284e+02 8.20578e-02 - 680.00 4.35321e+02 1.19454e+00 1.62408e-01 2.18886e-02 1.56498e+03 2.55644e+02 8.71198e-02 - 690.00 5.04157e+02 1.37009e+00 1.79217e-01 2.19816e-02 1.85688e+03 2.96943e+02 9.15828e-02 - 700.00 5.74706e+02 1.54853e+00 1.95474e-01 2.18680e-02 2.16794e+03 3.39463e+02 9.54178e-02 - 710.00 6.45681e+02 1.72668e+00 2.11024e-01 2.15699e-02 2.49374e+03 3.82446e+02 9.86077e-02 - 720.00 7.15785e+02 1.90135e+00 2.25734e-01 2.11107e-02 2.82933e+03 4.25116e+02 1.01147e-01 - 730.00 7.83757e+02 2.06950e+00 2.39502e-01 2.05143e-02 3.16942e+03 4.66716e+02 1.03042e-01 - 740.00 8.48418e+02 2.22833e+00 2.52250e-01 1.98043e-02 3.50857e+03 5.06525e+02 1.04307e-01 - 750.00 9.08702e+02 2.37527e+00 2.63924e-01 1.90035e-02 3.84134e+03 5.43889e+02 1.04967e-01 - 760.00 9.63684e+02 2.50823e+00 2.74491e-01 1.81332e-02 4.16251e+03 5.78228e+02 1.05051e-01 - 770.00 1.01260e+03 2.62546e+00 2.83941e-01 1.72132e-02 4.46718e+03 6.09058e+02 1.04598e-01 - 780.00 1.05486e+03 2.72563e+00 2.92276e-01 1.62613e-02 4.75092e+03 6.35992e+02 1.03648e-01 - 790.00 1.09003e+03 2.80785e+00 2.99518e-01 1.52932e-02 5.00989e+03 6.58744e+02 1.02245e-01 - 800.00 1.11787e+03 2.87163e+00 3.05698e-01 1.43228e-02 5.24086e+03 6.77133e+02 1.00437e-01 - 810.00 1.13829e+03 2.91686e+00 3.10857e-01 1.33618e-02 5.44130e+03 6.91071e+02 9.82693e-02 - 820.00 1.15134e+03 2.94378e+00 3.15045e-01 1.24199e-02 5.60935e+03 7.00563e+02 9.57902e-02 - 830.00 1.15722e+03 2.95295e+00 3.18317e-01 1.15052e-02 5.74386e+03 7.05695e+02 9.30460e-02 - 840.00 1.15623e+03 2.94522e+00 3.20732e-01 1.06240e-02 5.84432e+03 7.06625e+02 9.00815e-02 - 850.00 1.14877e+03 2.92161e+00 3.22352e-01 9.78106e-03 5.91085e+03 7.03570e+02 8.69393e-02 - 860.00 1.13531e+03 2.88337e+00 3.23239e-01 8.97989e-03 5.94413e+03 6.96795e+02 8.36597e-02 - 870.00 1.11638e+03 2.83181e+00 3.23457e-01 8.22282e-03 5.94531e+03 6.86605e+02 8.02802e-02 - 880.00 1.09255e+03 2.76838e+00 3.23066e-01 7.51114e-03 5.91600e+03 6.73331e+02 7.68350e-02 - 890.00 1.06441e+03 2.69455e+00 3.22126e-01 6.84533e-03 5.85811e+03 6.57319e+02 7.33554e-02 - 900.00 1.03257e+03 2.61178e+00 3.20696e-01 6.22346e-03 5.77386e+03 6.38924e+02 6.98695e-02 - 910.00 9.97596e+02 2.52155e+00 3.18830e-01 5.64994e-03 5.66565e+03 6.18502e+02 6.64019e-02 - 920.00 9.60081e+02 2.42526e+00 3.16581e-01 5.11810e-03 5.53603e+03 5.96399e+02 6.29744e-02 - 930.00 9.20573e+02 2.32427e+00 3.13997e-01 4.62818e-03 5.38760e+03 5.72952e+02 5.96059e-02 - 940.00 8.79591e+02 2.21986e+00 3.11123e-01 4.17830e-03 5.22300e+03 5.48479e+02 5.63122e-02 - 950.00 8.37617e+02 2.11321e+00 3.08002e-01 3.76640e-03 5.04482e+03 5.23279e+02 5.31067e-02 - 960.00 7.95094e+02 2.00540e+00 3.04671e-01 3.39029e-03 4.85560e+03 4.97627e+02 5.00003e-02 - 970.00 7.52424e+02 1.89740e+00 3.01165e-01 3.04773e-03 4.65775e+03 4.71774e+02 4.70016e-02 - 980.00 7.09961e+02 1.79009e+00 2.97517e-01 2.73647e-03 4.45357e+03 4.45946e+02 4.41172e-02 - 990.00 6.68019e+02 1.68421e+00 2.93755e-01 2.45427e-03 4.24517e+03 4.20342e+02 4.13518e-02 - 1000.00 6.26870e+02 1.58043e+00 2.89906e-01 2.19892e-03 4.03453e+03 3.95135e+02 3.87087e-02 - 1010.00 5.86742e+02 1.47931e+00 2.85992e-01 1.96832e-03 3.82342e+03 3.70476e+02 3.61895e-02 - 1020.00 5.47828e+02 1.38130e+00 2.82033e-01 1.76043e-03 3.61345e+03 3.46490e+02 3.37946e-02 - 1030.00 5.10281e+02 1.28677e+00 2.78049e-01 1.57333e-03 3.40603e+03 3.23282e+02 3.15236e-02 - 1040.00 4.74226e+02 1.19603e+00 2.74054e-01 1.40519e-03 3.20238e+03 3.00934e+02 2.93748e-02 - 1050.00 4.39753e+02 1.10930e+00 2.70063e-01 1.25431e-03 3.00358e+03 2.79511e+02 2.73459e-02 - 1060.00 4.06926e+02 1.02672e+00 2.66088e-01 1.11909e-03 2.81051e+03 2.59060e+02 2.54341e-02 - 1070.00 3.75785e+02 9.48386e-01 2.62139e-01 9.98042e-04 2.62389e+03 2.39612e+02 2.36358e-02 - 1080.00 3.46349e+02 8.74343e-01 2.58225e-01 8.89812e-04 2.44433e+03 2.21186e+02 2.19472e-02 - 1090.00 3.18617e+02 8.04583e-01 2.54353e-01 7.93136e-04 2.27225e+03 2.03787e+02 2.03642e-02 - 1100.00 2.92572e+02 7.39061e-01 2.50531e-01 7.06860e-04 2.10799e+03 1.87411e+02 1.88823e-02 - 1110.00 2.68184e+02 6.77698e-01 2.46763e-01 6.29929e-04 1.95177e+03 1.72044e+02 1.74971e-02 - 1120.00 2.45412e+02 6.20389e-01 2.43053e-01 5.61380e-04 1.80368e+03 1.57665e+02 1.62040e-02 - 1130.00 2.24206e+02 5.67007e-01 2.39405e-01 5.00341e-04 1.66377e+03 1.44248e+02 1.49982e-02 - 1140.00 2.04507e+02 5.17406e-01 2.35821e-01 4.46019e-04 1.53198e+03 1.31761e+02 1.38752e-02 - 1150.00 1.86253e+02 4.71427e-01 2.32304e-01 3.97698e-04 1.40820e+03 1.20166e+02 1.28304e-02 - 1160.00 1.69377e+02 4.28904e-01 2.28855e-01 3.54733e-04 1.29227e+03 1.09427e+02 1.18594e-02 - 1170.00 1.53809e+02 3.89663e-01 2.25475e-01 3.16542e-04 1.18395e+03 9.95019e+01 1.09577e-02 - 1180.00 1.39478e+02 3.53525e-01 2.22164e-01 2.82602e-04 1.08302e+03 9.03487e+01 1.01211e-02 - 1190.00 1.26312e+02 3.20312e-01 2.18922e-01 2.52446e-04 9.89185e+02 8.19248e+01 9.34560e-03 - 1200.00 1.14240e+02 2.89845e-01 2.15750e-01 2.25652e-04 9.02145e+02 7.41873e+01 8.62722e-03 - 1210.00 1.03191e+02 2.61948e-01 2.12646e-01 2.01848e-04 8.21584e+02 6.70937e+01 7.96222e-03 - 1220.00 9.30964e+01 2.36449e-01 2.09611e-01 1.80696e-04 7.47176e+02 6.06020e+01 7.34704e-03 - 1230.00 8.38898e+01 2.13181e-01 2.06643e-01 1.61899e-04 6.78587e+02 5.46717e+01 6.77829e-03 - 1240.00 7.55069e+01 1.91985e-01 2.03742e-01 1.45190e-04 6.15484e+02 4.92631e+01 6.25275e-03 - 1250.00 6.78859e+01 1.72705e-01 2.00906e-01 1.30335e-04 5.57535e+02 4.43383e+01 5.76739e-03 - 1260.00 6.09683e+01 1.55195e-01 1.98130e-01 1.17262e-04 5.04414e+02 3.98611e+01 5.31903e-03 - 1270.00 5.46984e+01 1.39317e-01 1.95414e-01 1.05646e-04 4.55802e+02 3.57970e+01 4.90530e-03 - 1280.00 4.90237e+01 1.24937e-01 1.92757e-01 9.50727e-05 4.11389e+02 3.21130e+01 4.52379e-03 - 1290.00 4.38947e+01 1.11933e-01 1.90160e-01 8.56932e-05 3.70878e+02 2.87785e+01 4.17206e-03 - 1300.00 3.92651e+01 1.00189e-01 1.87622e-01 7.73455e-05 3.33983e+02 2.57642e+01 3.84789e-03 - 1310.00 3.50918e+01 8.95957e-02 1.85139e-01 6.99040e-05 3.00431e+02 2.30432e+01 3.54919e-03 - 1320.00 3.13345e+01 8.00532e-02 1.82712e-01 6.32628e-05 2.69963e+02 2.05899e+01 3.27403e-03 - 1330.00 2.79558e+01 7.14672e-02 1.80339e-01 5.73299e-05 2.42335e+02 1.83808e+01 3.02062e-03 - 1340.00 2.49210e+01 6.37509e-02 1.78018e-01 5.20245e-05 2.17316e+02 1.63939e+01 2.78728e-03 - 1350.00 2.21984e+01 5.68239e-02 1.75748e-01 4.72756e-05 1.94688e+02 1.46089e+01 2.57247e-03 - 1360.00 1.97584e+01 5.06124e-02 1.73528e-01 4.30204e-05 1.74250e+02 1.30072e+01 2.37474e-03 - 1370.00 1.75741e+01 4.50483e-02 1.71356e-01 3.92036e-05 1.55812e+02 1.15715e+01 2.19276e-03 - 1380.00 1.56207e+01 4.00694e-02 1.69230e-01 3.57762e-05 1.39199e+02 1.02860e+01 2.02529e-03 - 1390.00 1.38756e+01 3.56186e-02 1.67151e-01 3.26950e-05 1.24247e+02 9.13605e+00 1.87120e-03 - 1400.00 1.23181e+01 3.16437e-02 1.65116e-01 2.99219e-05 1.10806e+02 8.10852e+00 1.72942e-03 - 1410.00 1.09293e+01 2.80972e-02 1.63124e-01 2.74231e-05 9.87358e+01 7.19123e+00 1.59898e-03 - 1420.00 9.69209e+00 2.49358e-02 1.61173e-01 2.51687e-05 8.79087e+01 6.37315e+00 1.47897e-03 - 1430.00 8.59095e+00 2.21202e-02 1.59264e-01 2.31324e-05 7.82066e+01 5.64421e+00 1.36856e-03 - 1440.00 7.61175e+00 1.96149e-02 1.57394e-01 2.12909e-05 6.95215e+01 4.99530e+00 1.26699e-03 - 1450.00 6.74172e+00 1.73874e-02 1.55562e-01 1.96233e-05 6.17545e+01 4.41814e+00 1.17354e-03 - 1460.00 5.96934e+00 1.54087e-02 1.53768e-01 1.81115e-05 5.48152e+01 3.90523e+00 1.08757e-03 - 1470.00 5.28418e+00 1.36522e-02 1.52010e-01 1.67391e-05 4.86213e+01 3.44980e+00 1.00846e-03 - 1480.00 4.67687e+00 1.20943e-02 1.50287e-01 1.54917e-05 4.30977e+01 3.04574e+00 9.35681e-04 - 1490.00 4.13896e+00 1.07135e-02 1.48599e-01 1.43566e-05 3.81763e+01 2.68754e+00 8.68707e-04 - 1500.00 3.66287e+00 9.49054e-03 1.46944e-01 1.33223e-05 3.37953e+01 2.37024e+00 8.07074e-04 - 1510.00 3.24178e+00 8.40815e-03 1.45321e-01 1.23789e-05 2.98986e+01 2.08938e+00 7.50350e-04 - 1520.00 2.86960e+00 7.45080e-03 1.43730e-01 1.15172e-05 2.64355e+01 1.84096e+00 6.98138e-04 - 1530.00 2.54085e+00 6.60460e-03 1.42170e-01 1.07293e-05 2.33603e+01 1.62139e+00 6.50074e-04 - 1540.00 2.25066e+00 5.85710e-03 1.40639e-01 1.00080e-05 2.06317e+01 1.42745e+00 6.05822e-04 - 1550.00 1.99465e+00 5.19718e-03 1.39137e-01 9.34691e-06 1.82125e+01 1.25627e+00 5.65075e-04 - 1560.00 1.76892e+00 4.61490e-03 1.37664e-01 8.74037e-06 1.60693e+01 1.10527e+00 5.27549e-04 - 1570.00 1.57000e+00 4.10140e-03 1.36218e-01 8.18326e-06 1.41718e+01 9.72171e-01 4.92983e-04 - 1580.00 1.39481e+00 3.64879e-03 1.34798e-01 7.67099e-06 1.24932e+01 8.54913e-01 4.61139e-04 - 1590.00 1.24058e+00 3.25004e-03 1.33405e-01 7.19948e-06 1.10092e+01 7.51676e-01 4.31797e-04 - 1600.00 1.10487e+00 2.89890e-03 1.32037e-01 6.76504e-06 9.69811e+00 6.60838e-01 4.04754e-04 - 1610.00 9.85502e-01 2.58981e-03 1.30694e-01 6.36435e-06 8.54061e+00 5.80955e-01 3.79826e-04 - 1620.00 8.80554e-01 2.31783e-03 1.29374e-01 5.99444e-06 7.51933e+00 5.10743e-01 3.56842e-04 - 1630.00 7.88315e-01 2.07858e-03 1.28079e-01 5.65261e-06 6.61881e+00 4.49065e-01 3.35645e-04 - 1640.00 7.07268e-01 1.86819e-03 1.26806e-01 5.33645e-06 5.82525e+00 3.94912e-01 3.16092e-04 - 1650.00 6.36076e-01 1.68321e-03 1.25555e-01 5.04376e-06 5.12635e+00 3.47387e-01 2.98051e-04 - 1660.00 5.73552e-01 1.52061e-03 1.24326e-01 4.77255e-06 4.51117e+00 3.05700e-01 2.81401e-04 - 1670.00 5.18650e-01 1.37770e-03 1.23118e-01 4.51430e-06 3.96997e+00 2.69148e-01 2.66030e-04 - 1680.00 4.70446e-01 1.25210e-03 1.21931e-01 4.28779e-06 3.49412e+00 2.37114e-01 2.51835e-04 - 1690.00 4.28124e-01 1.14171e-03 1.20764e-01 4.07094e-06 3.07594e+00 2.09051e-01 2.38723e-04 - 1700.00 3.90965e-01 1.04469e-03 1.19616e-01 3.86932e-06 2.70861e+00 1.84474e-01 2.26607e-04 - 1710.00 3.58335e-01 9.59401e-04 1.18487e-01 3.68171e-06 2.38611e+00 1.62958e-01 2.15408e-04 - 1720.00 3.29677e-01 8.84409e-04 1.17377e-01 3.50700e-06 2.10309e+00 1.44128e-01 2.05053e-04 - 1730.00 3.04500e-01 8.18446e-04 1.16285e-01 3.34418e-06 1.85482e+00 1.27653e-01 1.95475e-04 - 1740.00 2.82373e-01 7.60401e-04 1.15211e-01 3.19232e-06 1.63712e+00 1.13242e-01 1.86612e-04 - 1750.00 2.62916e-01 7.09294e-04 1.14155e-01 3.05057e-06 1.44629e+00 1.00639e-01 1.78408e-04 - 1760.00 2.45797e-01 6.64265e-04 1.13115e-01 2.91817e-06 1.27908e+00 8.96182e-02 1.70811e-04 - 1770.00 2.30724e-01 6.24560e-04 1.12092e-01 2.79438e-06 1.13261e+00 7.99830e-02 1.63773e-04 - 1780.00 2.17439e-01 5.89514e-04 1.11085e-01 2.67858e-06 1.00434e+00 7.15593e-02 1.57251e-04 - 1790.00 2.05719e-01 5.58545e-04 1.10094e-01 2.57016e-06 8.92035e-01 6.41948e-02 1.51202e-04 - 1800.00 1.95366e-01 5.31144e-04 1.09118e-01 2.46858e-06 7.93727e-01 5.77558e-02 1.45591e-04 - 1810.00 1.86208e-01 5.06864e-04 1.08157e-01 2.37332e-06 7.07685e-01 5.21255e-02 1.40384e-04 - 1820.00 1.78094e-01 4.85312e-04 1.07211e-01 2.28394e-06 6.32386e-01 4.72011e-02 1.35548e-04 - 1830.00 1.70892e-01 4.66147e-04 1.06280e-01 2.19999e-06 5.66493e-01 4.28930e-02 1.31056e-04 - 1840.00 1.64488e-01 4.49069e-04 1.05363e-01 2.12111e-06 5.08829e-01 3.91226e-02 1.26880e-04 - 1850.00 1.58779e-01 4.33815e-04 1.04459e-01 2.04691e-06 4.58364e-01 3.58211e-02 1.22996e-04 - 1860.00 1.53678e-01 4.20156e-04 1.03569e-01 1.97707e-06 4.14192e-01 3.29287e-02 1.19381e-04 - 1870.00 1.49108e-01 4.07891e-04 1.02692e-01 1.91129e-06 3.75521e-01 3.03928e-02 1.16015e-04 - 1880.00 1.45002e-01 3.96847e-04 1.01828e-01 1.84927e-06 3.41655e-01 2.81676e-02 1.12879e-04 - 1890.00 1.41301e-01 3.86870e-04 1.00976e-01 1.79077e-06 3.11984e-01 2.62131e-02 1.09955e-04 - 1900.00 1.37954e-01 3.77827e-04 1.00137e-01 1.73553e-06 2.85978e-01 2.44944e-02 1.07226e-04 - 1910.00 1.34917e-01 3.69600e-04 9.93102e-02 1.68334e-06 2.63168e-01 2.29812e-02 1.04679e-04 - 1920.00 1.32151e-01 3.62090e-04 9.84951e-02 1.63399e-06 2.43147e-01 2.16468e-02 1.02300e-04 - 1930.00 1.29621e-01 3.55208e-04 9.76915e-02 1.58728e-06 2.25560e-01 2.04683e-02 1.00075e-04 - 1940.00 1.27300e-01 3.48876e-04 9.68993e-02 1.54304e-06 2.10094e-01 1.94254e-02 9.79931e-05 - 1950.00 1.25160e-01 3.43027e-04 9.61182e-02 1.50111e-06 1.96479e-01 1.85007e-02 9.60438e-05 - 1960.00 1.23180e-01 3.37603e-04 9.53480e-02 1.46134e-06 1.84475e-01 1.76789e-02 9.42171e-05 - 1970.00 1.21340e-01 3.32552e-04 9.45885e-02 1.42358e-06 1.73878e-01 1.69467e-02 9.25038e-05 - 1980.00 1.19623e-01 3.27831e-04 9.38394e-02 1.38771e-06 1.64506e-01 1.62925e-02 9.08956e-05 - 1990.00 1.18015e-01 3.23400e-04 9.31006e-02 1.35360e-06 1.56201e-01 1.57064e-02 8.93848e-05 - 2000.00 1.16502e-01 3.19225e-04 9.23718e-02 1.32114e-06 1.48827e-01 1.51796e-02 8.79642e-05 - 2010.00 1.15074e-01 3.15278e-04 9.16528e-02 1.29022e-06 1.42265e-01 1.47044e-02 8.66273e-05 - 2020.00 1.13721e-01 3.11533e-04 9.09435e-02 1.26076e-06 1.36409e-01 1.42744e-02 8.53680e-05 - 2030.00 1.12434e-01 3.07967e-04 9.02437e-02 1.23266e-06 1.31169e-01 1.38838e-02 8.41806e-05 - 2040.00 1.11206e-01 3.04561e-04 8.95531e-02 1.20583e-06 1.26468e-01 1.35276e-02 8.30601e-05 - 2050.00 1.10031e-01 3.01297e-04 8.88716e-02 1.18021e-06 1.22235e-01 1.32015e-02 8.20016e-05 - 2060.00 1.08903e-01 0.00000e+00 8.81990e-02 1.15571e-06 1.18412e-01 1.29018e-02 8.10008e-05 diff --git a/tests/modules/testCONEX.cpp b/tests/modules/testCONEX.cpp index d04f3ac808c051defdd4de6c4ca6bafa9d59adad..9596a8acab4bbed0f42a059d855488900e7e8b97 100644 --- a/tests/modules/testCONEX.cpp +++ b/tests/modules/testCONEX.cpp @@ -21,6 +21,7 @@ #include <corsika/modules/CONEX.hpp> #include <corsika/modules/Sibyll.hpp> +#include <corsika/modules/writers/WriterOff.hpp> #include <corsika/framework/random/RNGManager.hpp> @@ -48,7 +49,7 @@ using MExtraEnvirnoment = MediumPropertyModel<UniformMagneticField<T>>; struct DummyStack {}; -TEST_CASE("CONEXSourceCut") { +TEST_CASE("CONEX") { logging::set_level(logging::level::info); @@ -93,8 +94,7 @@ TEST_CASE("CONEXSourceCut") { static_pow<2>(injectionHeight)); Point const showerCore{rootCS, 0_m, 0_m, observationHeight}; Point const injectionPos = - showerCore + - Vector<dimensionless_d>{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t; + showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t; ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.02, env}; @@ -102,7 +102,16 @@ TEST_CASE("CONEXSourceCut") { corsika::sibyll::Interaction sibyll; [[maybe_unused]] corsika::sibyll::NuclearInteractionModel sibyllNuc(sibyll, env); - CONEXhybrid conex(center, showerAxis, t, injectionHeight, E0, get_PDG(Code::Proton)); + EnergyLossWriter<WriterOff> w1(showerAxis); + LongitudinalWriter<WriterOff> w2(showerAxis); + CONEXhybrid<decltype(w1), decltype(w2)> conex(center, showerAxis, t, injectionHeight, + E0, get_PDG(Code::Proton), w1, w2); + // initialize writers + w1.startOfLibrary("test"); + w1.startOfShower(0); + w2.startOfLibrary("test"); + w2.startOfShower(0); + // init conex conex.initCascadeEquations(); HEPEnergyType const Eem{1_PeV}; @@ -120,45 +129,17 @@ TEST_CASE("CONEXSourceCut") { emPosition.getCoordinates(conex.getObserverCS()), emPosition.getCoordinates(rootCS)); - conex.addParticle(Code::Proton, Eem, 0_eV, emPosition, momentum.normalized(), 0_s); + conex.addParticle(Code::Proton, Eem, Proton::mass, emPosition, momentum.normalized(), + 0_s); // supperimpose a photon auto const momentumPhoton = showerAxis.getDirection() * 1_TeV; conex.addParticle(Code::Photon, 1_TeV, 0_eV, emPosition, momentumPhoton.normalized(), 0_s); DummyStack stack; conex.doCascadeEquations(stack); -} - -#include <algorithm> -#include <iterator> -#include <string> -#include <fstream> - -TEST_CASE("ConexOutput", "[output validation]") { - - logging::set_level(logging::level::info); - - auto file = GENERATE(as<std::string>{}, "conex_fit", "conex_output"); - - SECTION(std::string("check saved data, ") + file + ".txt") { - - // compare to reference data - std::ifstream file1(file + ".txt"); - std::ifstream file1ref(refDataDir + "/" + file + "_REF.txt"); - - std::istreambuf_iterator<char> begin1(file1); - std::istreambuf_iterator<char> begin1ref(file1ref); - std::istreambuf_iterator<char> end; + CHECK(w1.getEnergyLost() / 1_TeV == Approx(1.0).epsilon(0.1)); - while (begin1 != end && begin1ref != end) { - CHECK(*begin1 == *begin1ref); - ++begin1; - ++begin1ref; - } - CHECK(begin1 == end); - CHECK(begin1ref == end); - file1.close(); - file1ref.close(); - } + auto const cfg = conex.getConfig(); + CHECK(cfg.size() == 0); } diff --git a/tests/modules/testObservationPlane.cpp b/tests/modules/testObservationPlane.cpp index 6411d9738b5cef34e36290091d5adb57621b6d5a..5432c0946fcb9c56dd509ee80b710dfdfd08c93e 100644 --- a/tests/modules/testObservationPlane.cpp +++ b/tests/modules/testObservationPlane.cpp @@ -17,8 +17,6 @@ #include <corsika/framework/core/ParticleProperties.hpp> #include <corsika/framework/core/PhysicalUnits.hpp> -#include <corsika/output/NoOutput.hpp> - #include <SetupTestEnvironment.hpp> #include <SetupTestStack.hpp> #include <SetupTestTrajectory.hpp> @@ -55,8 +53,8 @@ TEST_CASE("ObservationPlane", "interface") { SECTION("horizontal plane") { Plane const obsPlane(Point(cs, {10_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.})); - ObservationPlane<setup::Tracking, NoOutput> obs(obsPlane, - DirectionVector(cs, {0., 1., 0.})); + ObservationPlane<setup::Tracking, WriterOff> obs(obsPlane, + DirectionVector(cs, {0., 1., 0.})); LengthType const length = obs.getMaxStepLength(particle, no_used_track); ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true); @@ -84,7 +82,7 @@ TEST_CASE("ObservationPlane", "interface") { SECTION("transparent plane") { Plane const obsPlane(Point(cs, {1_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.})); - ObservationPlane<setup::Tracking, NoOutput> obs( + ObservationPlane<setup::Tracking, WriterOff> obs( obsPlane, DirectionVector(cs, {0., 0., 1.}), false); LengthType const length = obs.getMaxStepLength(particle, no_used_track); @@ -105,8 +103,8 @@ TEST_CASE("ObservationPlane", "interface") { Plane const obsPlane(Point(cs, {10_m, 5_m, 5_m}), DirectionVector(cs, {1, 0.1, -0.05})); - ObservationPlane<setup::Tracking, NoOutput> obs(obsPlane, - DirectionVector(cs, {0., 1., 0.})); + ObservationPlane<setup::Tracking, WriterOff> obs(obsPlane, + DirectionVector(cs, {0., 1., 0.})); LengthType const length = obs.getMaxStepLength(particle, no_used_track); ProcessReturn const ret = obs.doContinuous(particle, no_used_track, true); @@ -117,7 +115,7 @@ TEST_CASE("ObservationPlane", "interface") { SECTION("output") { Plane const obsPlane(Point(cs, {1_m, 0_m, 0_m}), DirectionVector(cs, {1., 0., 0.})); - ObservationPlane<setup::Tracking, NoOutput> obs( + ObservationPlane<setup::Tracking, WriterOff> obs( obsPlane, DirectionVector(cs, {0., 0., 1.}), false); auto const cfg = obs.getConfig(); CHECK(cfg["type"]); diff --git a/tests/modules/testParticleCut.cpp b/tests/modules/testParticleCut.cpp index ef2b6c6365056a51a2ea974e548353a2683f5baa..e12ac460144b571c70951498179b464c40917819 100644 --- a/tests/modules/testParticleCut.cpp +++ b/tests/modules/testParticleCut.cpp @@ -52,7 +52,7 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") { SECTION("cut on particle type: inv") { // particle cut with 20GeV threshold for all, also cut invisible - ParticleCut cut(20_GeV, false, true); + ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true); CHECK(cut.getHadronKineticECut() == 20_GeV); // add primary particle to stack @@ -75,13 +75,11 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") { cut.doSecondaries(view); CHECK(view.getEntries() == 9); - CHECK(cut.getNumberInvParticles() == 2); - CHECK(cut.getInvEnergy() / 1_GeV == Approx(2000.)); } SECTION("cut on particle type: em") { - ParticleCut cut(20_GeV, true, false); + ParticleCut cut(1_EeV, 1_EeV, 1_GeV, 1_GeV, false); // add primary particle to stack auto particle = stack.addParticle(std::make_tuple( @@ -100,12 +98,10 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") { cut.doSecondaries(view); CHECK(view.getEntries() == 10); - CHECK(cut.getNumberEmParticles() == 1); - CHECK(cut.getEmEnergy() / 1_GeV == 1000.); } SECTION("cut low energy") { - ParticleCut cut(20_GeV, true, true); + ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true); // add primary particle to stack auto particle = stack.addParticle(std::make_tuple( @@ -167,16 +163,16 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") { } SECTION("cut low energy: reset thresholds of arbitrary set of particles") { - ParticleCut cut({{Code::Electron, 5_MeV}, {Code::Positron, 50_MeV}}, false, true); - CHECK(calculate_kinetic_energy_threshold(Code::Electron) != - calculate_kinetic_energy_threshold(Code::Positron)); - CHECK_FALSE(calculate_kinetic_energy_threshold(Code::Electron) == Electron::mass); + ParticleCut cut({{Code::Electron, 5_MeV}, {Code::Positron, 50_MeV}}, false); + CHECK(get_kinetic_energy_threshold(Code::Electron) != + get_kinetic_energy_threshold(Code::Positron)); + CHECK_FALSE(get_kinetic_energy_threshold(Code::Electron) == Electron::mass); // test default values still correct - CHECK(calculate_kinetic_energy_threshold(Code::Proton) == 5_GeV); + CHECK(get_kinetic_energy_threshold(Code::Proton) == 5_GeV); } SECTION("cut on time") { - ParticleCut cut(20_GeV, false, false); + ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, false); const TimeType too_late = 1_s; // add primary particle to stack @@ -196,20 +192,15 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") { cut.doSecondaries(view); CHECK(view.getEntries() == 0); - CHECK(cut.getTimeCutEnergy() == 11 * Eabove); - CHECK(cut.getCutEnergy() == 0_GeV); - cut.reset(); - CHECK(cut.getTimeCutEnergy() == 0_GeV); } setup::Trajectory const track = setup::testing::make_track<setup::Trajectory>( Line{point0, VelocityVector{rootCS, {0_m / second, 0_m / second, -constants::c}}}, 12_m / constants::c); - SECTION("cut on DoContinous, just invisibles") { + SECTION("cut on doContinous, just invisibles") { - ParticleCut cut(20_GeV, false, true); - CHECK(cut.getHadronKineticECut() == 20_GeV); + ParticleCut cut(20_GeV, 20_GeV, 20_GeV, 20_GeV, true); // add particles, all with energies above the threshold // only cut is by species @@ -224,7 +215,5 @@ TEST_CASE("ParticleCut", "process,continuous,secondary") { } CHECK(stack.getEntries() == 9); - CHECK(cut.getNumberInvParticles() == 2); - CHECK(cut.getInvEnergy() / 1_GeV == 2000); } } diff --git a/tests/output/CMakeLists.txt b/tests/output/CMakeLists.txt index 98a43b4c218bf8d9e21bb7b9e56ac1275cbc9b0e..f30b7b595682989f1fbf1c47bb043633a839e18e 100644 --- a/tests/output/CMakeLists.txt +++ b/tests/output/CMakeLists.txt @@ -5,6 +5,9 @@ set (test_output_sources testParquetStreamer.cpp testWriterObservationPlane.cpp testWriterTracks.cpp + testWriterEnergyLoss.cpp + testWriterLongitudinal.cpp + testWriterOff.cpp ) CORSIKA_ADD_TEST (testOutput SOURCES ${test_output_sources}) diff --git a/tests/output/testOutputManager.cpp b/tests/output/testOutputManager.cpp index c9194abfa210e0347e97c50a0dc4add954b03e25..ee886f526bda289f3ddbd352682cab341852a696 100644 --- a/tests/output/testOutputManager.cpp +++ b/tests/output/testOutputManager.cpp @@ -20,43 +20,38 @@ using namespace corsika; struct DummyNoOutput : public NoOutput { void check() { NoOutput::startOfLibrary("./"); - NoOutput::startOfShower(); - NoOutput::endOfShower(); + NoOutput::startOfShower(0); + NoOutput::endOfShower(0); NoOutput::endOfLibrary(); - NoOutput::getConfig(); - NoOutput::getSummary(); } void checkWrite() { NoOutput::write(Code::Unknown, 1_eV, 1_m, 1_m, 1_ns); } }; struct DummyOutput : public BaseOutput { - mutable bool isConfig_ = false; - mutable bool isSummary_ = false; bool startLibrary_ = false; bool startShower_ = false; bool endLibrary_ = false; bool endShower_ = false; - void startOfLibrary(boost::filesystem::path const&) { startLibrary_ = true; } + void startOfLibrary(boost::filesystem::path const&) override { startLibrary_ = true; } - void startOfShower() { - BaseOutput::startOfShower(); + YAML::Node getConfig() const final override { return YAML::Node(); } + + void startOfShower(unsigned int const shower = 0) override { + BaseOutput::startOfShower(shower); + setInit(true); startShower_ = true; } - void endOfShower() { endShower_ = true; } - - void endOfLibrary() { endLibrary_ = true; } + void endOfShower(unsigned int const) override { endShower_ = true; } - YAML::Node getConfig() const { - isConfig_ = true; - return YAML::Node(); - } + void endOfLibrary() override { endLibrary_ = true; } - YAML::Node getSummary() { - isSummary_ = true; - return BaseOutput::getSummary(); + YAML::Node getSummary() const final override { + YAML::Node summary; + summary["test"] = "test"; + return summary; } }; @@ -83,14 +78,13 @@ TEST_CASE("OutputManager") { "test", test)); // should emit warning which cannot be catched, but no action or failure - CHECK(test.isConfig_); - test.isConfig_ = false; - output.startOfLibrary(); CHECK(test.startLibrary_); test.startLibrary_ = false; + CHECK_FALSE(test.isInit()); output.startOfShower(); + CHECK(test.isInit()); CHECK(test.startShower_); test.startShower_ = false; @@ -100,11 +94,7 @@ TEST_CASE("OutputManager") { output.endOfLibrary(); CHECK(test.endLibrary_); - CHECK(test.isSummary_); - test.isSummary_ = false; test.endLibrary_ = false; - - CHECK(boost::filesystem::exists("./out_test/check/test/summary.yaml")); } SECTION("auto-write") { @@ -122,6 +112,11 @@ TEST_CASE("OutputManager") { DummyOutput test; output->add("test", test); output->startOfLibrary(); + + // cannot add more after library started + DummyOutput test2; + CHECK_THROWS(output->add("test2", test2)); + output->startOfShower(); // check support for closing automatically @@ -144,33 +139,21 @@ TEST_CASE("OutputManager") { OutputManager output("check", "./out_test"); CHECK_THROWS(new OutputManager("check", "./out_test")); - // CHECK_THROWS(output.startOfShower()); - // CHECK_THROWS(output.endOfShower()); CHECK_THROWS(output.endOfLibrary()); output.startOfLibrary(); CHECK_THROWS(output.startOfLibrary()); - // CHECK_THROWS(output.endOfShower()); - // CHECK_THROWS(output.endOfLibrary()); output.startOfShower(); CHECK_THROWS(output.startOfLibrary()); - // CHECK_THROWS(output.startOfShower()); - // CHECK_THROWS(output.endOfLibrary()); output.endOfShower(); CHECK_THROWS(output.startOfLibrary()); - // CHECK_THROWS(output.startOfShower()); - // CHECK_THROWS(output.endOfShower()); output.endOfLibrary(); - - // CHECK_THROWS(output.endOfShower()); - // CHECK_THROWS(output.startOfShower()); - // CHECK_THROWS(output.endOfLibrary()); } SECTION("NoOutput") { @@ -180,5 +163,7 @@ TEST_CASE("OutputManager") { nothing.check(); nothing.checkWrite(); + nothing.getConfig(); + nothing.getSummary(); } -} \ No newline at end of file +} diff --git a/tests/output/testParquetStreamer.cpp b/tests/output/testParquetStreamer.cpp index 387e5cade02d9e9cf1f1c1532cc1594bffb65d4a..13ade23599eafbb8403dad9d55a595c4aea220d4 100644 --- a/tests/output/testParquetStreamer.cpp +++ b/tests/output/testParquetStreamer.cpp @@ -36,17 +36,17 @@ TEST_CASE("ParquetStreamer") { test.addField("testfloat", parquet::Repetition::REQUIRED, parquet::Type::FLOAT, parquet::ConvertedType::NONE); - test.enableCompression(1); + // test.enableCompression(1); needs to be enabled via conan test.buildStreamer(); CHECK(test.isInit()); + unsigned int testId = 2; int testint = 1; double testfloat = 2.0; std::shared_ptr<parquet::StreamWriter> writer = test.getWriter(); - (*writer) << static_cast<int>(testint) << static_cast<int>(testint) - << static_cast<float>(testfloat) << parquet::EndRow; + (*writer) << testId << testint << static_cast<float>(testfloat) << parquet::EndRow; test.closeStreamer(); CHECK_THROWS(test.getWriter()); diff --git a/tests/output/testWriterEnergyLoss.cpp b/tests/output/testWriterEnergyLoss.cpp new file mode 100644 index 0000000000000000000000000000000000000000..43a09f441fdef0c6bbb8e1c0277c96e8722ed271 --- /dev/null +++ b/tests/output/testWriterEnergyLoss.cpp @@ -0,0 +1,132 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <catch2/catch.hpp> + +#include <boost/filesystem.hpp> + +#include <corsika/modules/writers/EnergyLossWriter.hpp> + +#include <corsika/media/HomogeneousMedium.hpp> +#include <corsika/media/ShowerAxis.hpp> + +#include <corsika/framework/geometry/StraightTrajectory.hpp> +#include <corsika/framework/geometry/Point.hpp> +#include <corsika/framework/geometry/Line.hpp> +#include <corsika/framework/geometry/CoordinateSystem.hpp> + +#include <corsika/framework/core/Logging.hpp> + +using namespace corsika; + +const auto density = 1_kg / (1_m * 1_m * 1_m); + +auto setupEnvironment(Code vTargetCode) { + // setup environment, geometry + auto env = std::make_unique<Environment<IMediumModel>>(); + auto& universe = *(env->getUniverse()); + const CoordinateSystemPtr& cs = env->getCoordinateSystem(); + + auto theMedium = Environment<IMediumModel>::createNode<Sphere>( + Point{cs, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity()); + + using MyHomogeneousModel = HomogeneousMedium<IMediumModel>; + theMedium->setModelProperties<MyHomogeneousModel>( + density, NuclearComposition({vTargetCode}, {1.})); + + auto const* nodePtr = theMedium.get(); + universe.addChild(std::move(theMedium)); + + return std::make_tuple(std::move(env), &cs, nodePtr); +} + +class TestEnergyLoss : public corsika::EnergyLossWriter<> { +public: + TestEnergyLoss(corsika::ShowerAxis const& axis) + : EnergyLossWriter(axis) {} +}; + +TEST_CASE("EnergyLossWriter") { + + // best to run unit-tests with TRACE output + logging::set_level(logging::level::trace); + + auto [env, csPtr, nodePtr] = setupEnvironment(Code::Nitrogen); + auto const& cs = *csPtr; + [[maybe_unused]] auto const& env_dummy = env; + [[maybe_unused]] auto const& node_dummy = nodePtr; + + auto const observationHeight = 0_km; + auto const injectionHeight = 10_km; + auto const t = -observationHeight + injectionHeight; + Point const showerCore{cs, 0_m, 0_m, observationHeight}; + Point const injectionPos = showerCore + DirectionVector{cs, {0, 0, 1}} * t; + + ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos), *env, + false, // -> throw exceptions + 1000}; // -> number of bins + + // preparation + if (boost::filesystem::exists("./output_dir_eloss")) { + boost::filesystem::remove_all("./output_dir_eloss"); + } + boost::filesystem::create_directory("./output_dir_eloss"); + + TestEnergyLoss test(showerAxis); + test.startOfLibrary("./output_dir_eloss"); + test.startOfShower(0); + + CHECK(test.getEnergyLost() / 1_GeV == Approx(0)); + test.write(100_g / square(1_cm), 110_g / square(1_cm), Code::Photon, 100_GeV); + CHECK(test.getEnergyLost() / 1_GeV == Approx(100)); + test.write(Point(cs, {0_m, 0_m, observationHeight + 10_km}), Code::Proton, 100_GeV); + CHECK(test.getEnergyLost() / 1_GeV == Approx(200)); + + // generate straight simple track + CoordinateSystemPtr rootCS = get_root_CoordinateSystem(); + Point r0(rootCS, {0_m, 0_m, 6.555_km}); + SpeedType const V0 = constants::c; + VelocityVector v0(rootCS, {0_m / second, 0_m / second, -V0}); + Line const line(r0, v0); + auto const time = 1000_ns; + StraightTrajectory track(line, time); + StraightTrajectory trackShort(line, time / 5e3); // short + StraightTrajectory trackPointLike(line, time / 1e7); // ultra short + StraightTrajectory trackInverse({track.getPosition(1), -v0}, time); + // test write + test.write(track, Code::Proton, 100_GeV); + test.write(trackInverse, Code::Proton, 100_GeV); // equivalent + test.write(trackShort, Code::Proton, 100_GeV); // this is in a single bin + test.write(trackPointLike, Code::Proton, + 100_GeV); // this is just a located point-like dE + + // incompatible binning + CHECK_THROWS(test.write(100_g / square(1_cm), // extra line break by purpose + 120_g / square(1_cm), Code::Photon, 100_GeV)); + test.write(100000_g / square(1_cm), 100010_g / square(1_cm), Code::Photon, + 100_GeV); // this doesn't throw, but it skips + + test.endOfShower(0); + test.endOfLibrary(); + + CHECK(boost::filesystem::exists("./output_dir_eloss/dEdX.parquet")); + + auto const config = test.getConfig(); + CHECK(config["type"].as<std::string>() == "EnergyLoss"); + CHECK(config["units"]["energy"].as<std::string>() == "GeV"); + CHECK(config["units"]["grammage"].as<std::string>() == "g/cm^2"); + CHECK(config["bin-size"].as<double>() == 10.); + CHECK(config["nbins"].as<int>() == 200); + CHECK(config["grammage_threshold"].as<double>() == Approx(0.0001)); + + auto const summary = test.getSummary(); + CHECK(summary["sum_dEdX"].as<double>() == 600); + // makes not yet sense: + // CHECK(summary["Xmax"].as<double>() == 200); + // CHECK(summary["dEdXmax"].as<double>() == 200); +} diff --git a/tests/output/testWriterLongitudinal.cpp b/tests/output/testWriterLongitudinal.cpp new file mode 100644 index 0000000000000000000000000000000000000000..beb00767635a7709773a3d9c6381b0fb424d863f --- /dev/null +++ b/tests/output/testWriterLongitudinal.cpp @@ -0,0 +1,126 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <catch2/catch.hpp> + +#include <boost/filesystem.hpp> + +#include <corsika/modules/writers/LongitudinalWriter.hpp> + +#include <corsika/media/HomogeneousMedium.hpp> +#include <corsika/media/ShowerAxis.hpp> + +#include <corsika/framework/geometry/StraightTrajectory.hpp> +#include <corsika/framework/geometry/Point.hpp> +#include <corsika/framework/geometry/Line.hpp> +#include <corsika/framework/geometry/CoordinateSystem.hpp> + +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/framework/core/Logging.hpp> + +#include <string> + +using namespace corsika; + +const auto density = 1_kg / (1_m * 1_m * 1_m); + +auto setupEnvironment2(Code vTargetCode) { + // setup environment, geometry + auto env = std::make_unique<Environment<IMediumModel>>(); + auto& universe = *(env->getUniverse()); + const CoordinateSystemPtr& cs = env->getCoordinateSystem(); + + auto theMedium = Environment<IMediumModel>::createNode<Sphere>( + Point{cs, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity()); + + using MyHomogeneousModel = HomogeneousMedium<IMediumModel>; + theMedium->setModelProperties<MyHomogeneousModel>( + density, NuclearComposition({vTargetCode}, {1.})); + + auto const* nodePtr = theMedium.get(); + universe.addChild(std::move(theMedium)); + + return std::make_tuple(std::move(env), &cs, nodePtr); +} + +class TestLongitudinal : public corsika::LongitudinalWriter<> { +public: + TestLongitudinal(corsika::ShowerAxis const& axis) + : LongitudinalWriter(axis) {} +}; + +TEST_CASE("LongitudinalWriter") { + + logging::set_level(logging::level::info); + + auto [env, csPtr, nodePtr] = setupEnvironment2(Code::Nitrogen); + auto const& cs = *csPtr; + [[maybe_unused]] auto const& env_dummy = env; + [[maybe_unused]] auto const& node_dummy = nodePtr; + + auto const observationHeight = 0_km; + auto const injectionHeight = 10_km; + auto const t = -observationHeight + injectionHeight; + Point const showerCore{cs, 0_m, 0_m, observationHeight}; + Point const injectionPos = showerCore + DirectionVector{cs, {0, 0, 1}} * t; + + ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos), *env, + false, // -> throw exceptions + 1000}; // -> number of bins + + // preparation + if (boost::filesystem::exists("./output_dir_long")) { + boost::filesystem::remove_all("./output_dir_long"); + } + boost::filesystem::create_directory("./output_dir_long"); + + TestLongitudinal test(showerAxis); + test.startOfLibrary("./output_dir_long"); + test.startOfShower(0); + + // generate straight simple track + CoordinateSystemPtr rootCS = get_root_CoordinateSystem(); + Point r0(rootCS, {0_km, 0_m, 8_km}); + SpeedType const V0 = constants::c; + VelocityVector v0(rootCS, {0_m / second, 0_m / second, -V0}); + Line const line(r0, v0); + auto const time = 1000_ns; + StraightTrajectory track(line, time); + // test write + test.write(track, Code::Proton, 1.0); + test.write(track, Code::Photon, 1.0); + test.write(track, Code::Electron, 1.0); + test.write(track, Code::Positron, 1.0); + test.write(track, Code::MuPlus, 1.0); + test.write(track, Code::MuMinus, 1.0); + + test.write(10_g / square(1_cm), 20_g / square(1_cm), Code::PiPlus, 1.0); + test.write(10_g / square(1_cm), 20_g / square(1_cm), Code::Electron, 1.0); + test.write(10_g / square(1_cm), 20_g / square(1_cm), Code::Positron, 1.0); + test.write(10_g / square(1_cm), 20_g / square(1_cm), Code::Photon, 1.0); + test.write(10_g / square(1_cm), 20_g / square(1_cm), Code::MuPlus, 1.0); + test.write(10_g / square(1_cm), 20_g / square(1_cm), Code::MuMinus, 1.0); + + // wrong binning + CHECK_THROWS(test.write(10_g / square(1_cm), 10.1_g / square(1_cm), Code::PiPlus, 1.0)); + test.write(100000_g / square(1_cm), 100010_g / square(1_cm), Code::PiPlus, + 1.0); // this doesn't throw, it just skips + + test.endOfShower(0); + test.endOfLibrary(); + + CHECK(boost::filesystem::exists("./output_dir_long/profile.parquet")); + + auto const config = test.getConfig(); + CHECK(config["type"].as<std::string>() == "LongitudinalProfile"); + CHECK(config["units"]["grammage"].as<std::string>() == "g/cm^2"); + CHECK(config["bin-size"].as<double>() == 10.); + CHECK(config["nbins"].as<int>() == 200); + + auto const summary = test.getSummary(); // nothing to check yet +} diff --git a/tests/output/testWriterObservationPlane.cpp b/tests/output/testWriterObservationPlane.cpp index 7289ccc81e142ad7fbce0bd4fb63cb5fb24c0b2d..c58861f52cd3f80bddf9ebd51ea4554aa39a1491 100644 --- a/tests/output/testWriterObservationPlane.cpp +++ b/tests/output/testWriterObservationPlane.cpp @@ -10,19 +10,23 @@ #include <boost/filesystem.hpp> -#include <corsika/modules/writers/ObservationPlaneWriterParquet.hpp> +#include <corsika/modules/writers/ParticleWriterParquet.hpp> #include <corsika/framework/core/Logging.hpp> #include <corsika/framework/geometry/QuantityVector.hpp> using namespace corsika; -struct TestWriterPlane : public ObservationPlaneWriterParquet { +struct TestWriterPlane : public ParticleWriterParquet { YAML::Node getConfig() const { return YAML::Node(); } void checkWrite() { - ObservationPlaneWriterParquet::write(Code::Unknown, 1_eV, 2_m, 3_m, 4_ns); + ParticleWriterParquet::write(Code::Unknown, 1_GeV, 2_m, 3_m, 0_m, 1.0); + ParticleWriterParquet::write(Code::Proton, 1_GeV, 2_m, 3_m, 0_m, 1.0); + ParticleWriterParquet::write(Code::MuPlus, 1_GeV, 2_m, 3_m, 0_m, 1.0); + ParticleWriterParquet::write(Code::MuMinus, 1_GeV, 2_m, 3_m, 0_m, 1.0); + ParticleWriterParquet::write(Code::Photon, 1_GeV, 2_m, 3_m, 0_m, 1.0); } }; @@ -40,11 +44,22 @@ TEST_CASE("ObservationPlaneWriterParquet") { TestWriterPlane test; test.startOfLibrary("./output_dir"); - test.startOfShower(); + test.startOfShower(0); + + // write a few particles test.checkWrite(); - test.endOfShower(); + + test.endOfShower(0); test.endOfLibrary(); CHECK(boost::filesystem::exists("./output_dir/particles.parquet")); + + auto const summary = test.getSummary(); + + CHECK(summary["Eground"].as<double>() == Approx(5)); + CHECK(summary["hadrons"].as<int>() == Approx(1)); + CHECK(summary["muons"].as<int>() == Approx(2)); + CHECK(summary["em"].as<int>() == Approx(1)); + CHECK(summary["others"].as<int>() == Approx(1)); } } diff --git a/tests/output/testWriterOff.cpp b/tests/output/testWriterOff.cpp new file mode 100644 index 0000000000000000000000000000000000000000..1f572da3f42d19a69db53c72b72cb7758e8f6690 --- /dev/null +++ b/tests/output/testWriterOff.cpp @@ -0,0 +1,44 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * This software is distributed under the terms of the GNU General Public + * Licence version 3 (GPL Version 3). See file LICENSE for a full version of + * the license. + */ + +#include <catch2/catch.hpp> + +#include <boost/filesystem.hpp> + +#include <corsika/modules/writers/WriterOff.hpp> + +#include <corsika/framework/core/Logging.hpp> + +using namespace corsika; + +/* +class TestEnergyLoss : public corsika::EnergyLossWriter<> { +public: + TestEnergyLoss(corsika::ShowerAxis const& axis) + : EnergyLossWriter(axis) {} + + YAML::Node getConfig() const { return YAML::Node(); } +}; +*/ + +TEST_CASE("WriterOff") { + + logging::set_level(logging::level::info); + + WriterOff test("irrelevant", 3); + WriterOff test2(); + + test.startOfLibrary("./output_dir_eloss"); + test.startOfShower(0); + test.endOfShower(0); + test.endOfLibrary(); + + auto const config = test.getConfig(); + + auto const summary = test.getSummary(); +} diff --git a/tests/output/testWriterTracks.cpp b/tests/output/testWriterTracks.cpp index 43fbf88704e423027a157c099f472e44b9ff63f9..f1d109fbe2d7c6ddeed81e2fa35552a674b6618a 100644 --- a/tests/output/testWriterTracks.cpp +++ b/tests/output/testWriterTracks.cpp @@ -21,8 +21,8 @@ struct TestWriterTrack : public TrackWriterParquet { YAML::Node getConfig() const { return YAML::Node(); } void checkWrite() { - TrackWriterParquet::write(Code::Unknown, 1_eV, {2_m, 3_m, 4_m}, 5_s, {6_m, 7_m, 8_m}, - 9_s); + TrackWriterParquet::write(Code::Unknown, 1_eV, 1.0, {2_m, 3_m, 4_m}, 1_ns, + {5_m, 6_m, 7_m}, 2_ns); } }; @@ -40,9 +40,9 @@ TEST_CASE("TrackWriterParquet") { TestWriterTrack test; test.startOfLibrary("./output_dir_tracks"); - test.startOfShower(); + test.startOfShower(0); test.checkWrite(); - test.endOfShower(); + test.endOfShower(0); test.endOfLibrary(); CHECK(boost::filesystem::exists("./output_dir_tracks/tracks.parquet"));