IAP GITLAB

Commit 9ef3d775 authored by Ralf Ulrich's avatar Ralf Ulrich

logging, test, examples

parent 4b28369d
Pipeline #3172 passed with stages
in 20 minutes and 53 seconds
......@@ -337,8 +337,11 @@ namespace corsika::pythia8 {
// link to pythia stack
Pythia8::Event& event = pythia_.event;
// print final state
event.list();
if (print_listing_) {
// print final state
event.list();
}
MomentumVector Plab_final(labCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
HEPEnergyType Elab_final = 0_GeV;
......
......@@ -153,10 +153,12 @@ namespace corsika::qgsjetII {
event is copied (and boosted) into the shower lab frame.
*/
template <typename TParticle>
void Interaction::doInteraction(TParticle& vP) {
template <typename TView>
void Interaction::doInteraction(TView& view) {
auto const projectile = view.getProjectile();
const auto corsikaBeamId = vP.getPID();
const auto corsikaBeamId = projectile.getPID();
std::cout << "ProcessQgsjetII: "
<< "DoInteraction: " << corsikaBeamId << " interaction? "
<< corsika::qgsjetII::canInteract(corsikaBeamId) << std::endl;
......@@ -166,22 +168,24 @@ namespace corsika::qgsjetII {
CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
// position and time of interaction, not used in QgsjetII
Point pOrig = vP.getPosition();
TimeType tOrig = vP.getTime();
Point pOrig = projectile.getPosition();
TimeType tOrig = projectile.getTime();
// define target
// for QgsjetII is always a single nucleon
// FOR NOW: target is always at rest
const auto targetEnergyLab = 0_GeV + constants::nucleonMass;
const auto targetMomentumLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
const FourVector PtargLab(targetEnergyLab, targetMomentumLab);
auto const targetEnergyLab = 0_GeV + constants::nucleonMass;
auto const targetMomentumLab = MomentumVector(rootCS, 0_GeV, 0_GeV, 0_GeV);
FourVector const PtargLab(targetEnergyLab, targetMomentumLab);
// define projectile
HEPEnergyType const projectileEnergyLab = vP.getEnergy();
auto const projectileMomentumLab = vP.getMomentum();
HEPEnergyType const projectileEnergyLab = projectile.getEnergy();
auto const projectileMomentumLab = projectile.getMomentum();
int beamA = 0;
if (is_nucleus(corsikaBeamId)) beamA = vP.getNuclearA();
if (is_nucleus(corsikaBeamId)) beamA = projectile.getNuclearA();
HEPEnergyType const projectileEnergyLabPerNucleon = projectileEnergyLab / beamA;
std::cout << "Interaction: ebeam lab: " << projectileEnergyLab / 1_GeV << std::endl
<< "Interaction: pbeam lab: "
......@@ -195,7 +199,7 @@ namespace corsika::qgsjetII {
std::cout << "Interaction: time: " << tOrig << std::endl;
// sample target mass number
auto const* currentNode = vP.getNode();
auto const* currentNode = projectile.getNode();
auto const& mediumComposition =
currentNode->getModelProperties().getNuclearComposition();
// get cross sections for target materials
......@@ -219,19 +223,37 @@ namespace corsika::qgsjetII {
mediumComposition.sampleTarget(cross_section_of_components, rng_);
std::cout << "Interaction: target selected: " << targetCode << std::endl;
int targetQgsCode = -1;
if (is_nucleus(targetCode)) targetQgsCode = get_nucleus_A(targetCode);
if (targetCode == Code::Proton) targetQgsCode = 1;
std::cout << "Interaction: target qgsjetII code/A: " << targetQgsCode << std::endl;
if (targetQgsCode > maxMassNumber_ || targetQgsCode < 1)
throw std::runtime_error("QgsjetII target outside range.");
int projQgsCode = 1;
if (is_nucleus(corsikaBeamId)) projQgsCode = vP.getNuclearA();
std::cout << "Interaction: projectile qgsjetII code/A: " << projQgsCode << " "
<< corsikaBeamId << std::endl;
if (projQgsCode > maxMassNumber_ || projQgsCode < 1)
throw std::runtime_error("QgsjetII target outside range.");
int targetMassNumber = 1; // proton
if (is_nucleus(targetCode)) { // nucleus
targetMassNumber = get_nucleus_A(targetCode);
if (targetMassNumber > maxMassNumber_)
throw std::runtime_error("QgsjetII target mass outside range.");
} else {
if (targetCode != Proton::code)
throw std::runtime_error("QgsjetII Taget not possible.");
}
std::cout << "Interaction: target qgsjetII code/A: " << targetMassNumber << std::endl;
int projectileMassNumber = 1; // "1" means "hadron"
QgsjetIIHadronType qgsjet_hadron_type =
qgsjetII::getQgsjetIIHadronType(corsikaBeamId);
if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
projectileMassNumber = projectile.getNuclearA();
if (projectileMassNumber > maxMassNumber_)
throw std::runtime_error("QgsjetII projectile mass outside range.");
std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
std::uniform_int_distribution select(0, 1);
qgsjet_hadron_type = nucleons[select(rng_)];
} else {
// from conex: replace pi0 or rho0 with pi+/pi- in alternating sequence
if (qgsjet_hadron_type == QgsjetIIHadronType::NeutralLightMesonType) {
qgsjet_hadron_type = alternate_;
alternate_ = (alternate_ == QgsjetIIHadronType::PiPlusType
? QgsjetIIHadronType::PiMinusType
: QgsjetIIHadronType::PiPlusType);
}
}
// beam id for qgsjetII
int kBeam = 2; // default: proton Shouldn't we randomize neutron/proton for nuclei?
......@@ -254,9 +276,9 @@ namespace corsika::qgsjetII {
std::cout << "Interaction: "
<< " DoInteraction: E(GeV):" << projectileEnergyLab / 1_GeV << std::endl;
count_++;
qgini_(projectileEnergyLab / 1_GeV, kBeam, projQgsCode, targetQgsCode);
// this is from CRMC, is this REALLY needed ???
qgini_(projectileEnergyLab / 1_GeV, kBeam, projQgsCode, targetQgsCode);
int qgsjet_hadron_type_int = static_cast<QgsjetIICodeIntType>(qgsjet_hadron_type);
qgini_(projectileEnergyLab / 1_GeV, qgsjet_hadron_type_int, projectileMassNumber,
targetMassNumber);
qgconf_();
// bookkeeping
......@@ -278,21 +300,30 @@ namespace corsika::qgsjetII {
int Z = 0;
switch (A) {
case 1: { // proton/neutron
idFragm = Code::Proton;
auto momentum =
Vector(zAxisFrame, QuantityVector<hepmomentum_d>{
0.0_GeV, 0.0_GeV,
sqrt((projectileEnergyLab + Proton::mass) *
(projectileEnergyLab - Proton::mass))});
auto const energy =
sqrt(momentum.getSquaredNorm() + square(get_mass(idFragm)));
std::uniform_real_distribution<double> select;
if (select(rng_) > 0.5) {
idFragm = Code::Proton;
Z = 1;
} else {
idFragm = Code::Neutron;
Z = 0;
}
const HEPMassType nucleonMass = get_mass(idFragm);
auto momentum = Vector(
zAxisFrame, corsika::QuantityVector<hepmomentum_d>{
0.0_GeV, 0.0_GeV,
sqrt((projectileEnergyLabPerNucleon + nucleonMass) *
(projectileEnergyLabPerNucleon - nucleonMass))});
auto const energy = sqrt(momentum.getSquaredNorm() + square(nucleonMass));
momentum.rebase(originalCS); // transform back into standard lab frame
std::cout << "secondary fragment> id=" << idFragm
<< " p=" << momentum.getComponents() << std::endl;
auto pnew =
vP.addSecondary(std::make_tuple(idFragm, energy, momentum, pOrig, tOrig));
auto pnew = view.addSecondary(
std::make_tuple(idFragm, energy, momentum,
pOrig, tOrig));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
} break;
......@@ -311,21 +342,23 @@ namespace corsika::qgsjetII {
}
}
if (idFragm == Code::Nucleus) {
if (idFragm == Code::Nucleus) { // thus: not p or n
const HEPMassType nucleusMass =
Proton::mass * Z + Neutron::mass * (A - Z);
auto momentum = Vector(
zAxisFrame, QuantityVector<hepmomentum_d>{
0.0_GeV, 0.0_GeV,
sqrt((projectileEnergyLab + constants::nucleonMass * A) *
(projectileEnergyLab - constants::nucleonMass * A))});
sqrt((projectileEnergyLabPerNucleon * A + nucleusMass) *
(projectileEnergyLabPerNucleon * A - nucleusMass))});
auto const energy =
sqrt(momentum.getSquaredNorm() + square(constants::nucleonMass * A));
auto const energy = sqrt(momentum.getSquaredNorm() + square(nucleusMass));
momentum.rebase(originalCS); // transform back into standard lab frame
std::cout << "secondary fragment> id=" << idFragm
<< " p=" << momentum.getComponents() << " A=" << A << " Z=" << Z
<< std::endl;
auto pnew = vP.addSecondary(
std::make_tuple(idFragm, energy, momentum, pOrig, tOrig, A, Z));
auto pnew = view.addSecondary(
std::make_tuple(
idFragm, energy, momentum, pOrig, tOrig, A, Z));
Plab_final += pnew.getMomentum();
Elab_final += pnew.getEnergy();
}
......@@ -342,7 +375,7 @@ namespace corsika::qgsjetII {
std::cout << "secondary fragment> id="
<< corsika::qgsjetII::convertFromQgsjetII(psec.getPID())
<< " p=" << momentum.getComponents() << std::endl;
auto pnew = vP.addSecondary(
auto pnew = view.addSecondary(
std::make_tuple(corsika::qgsjetII::convertFromQgsjetII(psec.getPID()), energy,
momentum, pOrig, tOrig));
Plab_final += pnew.getMomentum();
......
......@@ -46,6 +46,11 @@
namespace corsika {
/*
* The default pattern for CORSIKA8 loggers.
*/
const std::string default_pattern{"[%n:%^%-8l%$] %v"};
/**
* Create a new C8-style logger.
*
......@@ -83,7 +88,8 @@ namespace corsika {
/**
* The default "corsika" logger.
*/
inline std::shared_ptr<spdlog::logger> corsika_logger = get_logger("corsika", true);
static inline std::shared_ptr<spdlog::logger> corsika_logger = get_logger("corsika", true);
// corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// many of these free functions are special to the logging
// infrastructure so we hide them in the corsika::logging namespace.
......@@ -140,4 +146,4 @@ namespace corsika {
} // namespace corsika
#include <corsika/detail/framework/logging/Logging.inl>
#include <corsika/detail/framework/core/Logging.inl>
......@@ -27,7 +27,7 @@ namespace corsika {
* RootCoordinateSystem
*/
static CoordinateSystemPtr& get_root_CoordinateSystem() {
static inline CoordinateSystemPtr& get_root_CoordinateSystem() {
static CoordinateSystemPtr rootCS(new CoordinateSystem); // THIS IS IT
return rootCS;
}
......
......@@ -8,4 +8,12 @@
#pragma once
namespace corsika {}
namespace corsika::setup {
/**
* The default "corsika" logger.
*/
static std::shared_ptr<spdlog::logger> corsika_logger = get_logger("corsika", true);
// corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
}
......@@ -79,7 +79,8 @@ private:
//
int main() {
// logging::SetLevel(logging::level::info);
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CORSIKA_LOG_INFO("boundary_example");
......
......@@ -35,7 +35,7 @@
/*
NOTE, WARNING, ATTENTION
The .../Random.hpppp implement the hooks of external modules to the C8 random
The .../Random.hpp implement the hooks of external modules to the C8 random
number generator. It has to occur excatly ONCE per linked
executable. If you include the header below multiple times and
link this togehter, it will fail.
......@@ -54,7 +54,9 @@ using namespace std;
//
int main() {
logging::set_level(logging::level::info);
//logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::trace);
std::cout << "cascade_example" << std::endl;
......
......@@ -58,6 +58,7 @@ using namespace std;
int main() {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
std::cout << "cascade_proton_example" << std::endl;
......
......@@ -65,6 +65,7 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
if (argc != 2) {
std::cerr << "usage: em_shower <energy/GeV>" << std::endl;
......
......@@ -20,6 +20,9 @@ using namespace corsika;
int main() {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CORSIKA_LOG_INFO("geometry_example");
// define the root coordinate system
......
......@@ -20,6 +20,9 @@ using namespace corsika;
int main() {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CORSIKA_LOG_INFO("helix_example");
CoordinateSystemPtr const& root = get_root_CoordinateSystem();
......
......@@ -85,6 +85,7 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CORSIKA_LOG_INFO("hybrid_MC");
......
......@@ -34,7 +34,8 @@ using namespace std;
//
int main() {
logging::set_level(spdlog::level::info);
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::info(
"------------------------------------------\n"
......
......@@ -43,6 +43,9 @@ void read(VectorStack& s) {
int main() {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
std::cout << "stack_example" << std::endl;
VectorStack s;
fill(s);
......
......@@ -104,6 +104,9 @@ void modular() {
int main() {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
std::cout << "staticsequence_example" << std::endl;
modular();
......
......@@ -29,7 +29,10 @@ using namespace std;
//
int main() {
std::cout << "stopping_power" << std::endl;
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CORSIKA_LOG_INFO("stopping_power");
feenableexcept(FE_INVALID);
......
......@@ -90,6 +90,7 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) {
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::info);
CORSIKA_LOG_INFO("vertical_EAS");
......
......@@ -38,6 +38,10 @@ auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
};
TEST_CASE("rotation") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1_GeV;
HEPMassType const targetMass = 1.0e300_eV;
......@@ -166,6 +170,10 @@ TEST_CASE("rotation") {
}
TEST_CASE("boosts") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// define target kinematics in lab frame
HEPMassType const targetMass = 1_GeV;
MomentumVector pTargetLab{rootCS, {0_eV, 0_eV, 0_eV}};
......@@ -319,6 +327,10 @@ TEST_CASE("boosts") {
}
TEST_CASE("rest frame") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
HEPMassType const projectileMass = 1_GeV;
HEPMomentumType const P0 = 1_TeV;
MomentumVector pProjectileLab{rootCS, {0_GeV, P0, 0_GeV}};
......
......@@ -15,6 +15,7 @@
#include <corsika/modules/StackInspector.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
......@@ -136,6 +137,7 @@ public:
TEST_CASE("Cascade", "[Cascade]") {
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::trace);
HEPEnergyType E0 = 100_GeV;
......
......@@ -7,6 +7,7 @@
*/
#include <corsika/framework/analytics/ClassTimer.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <catch2/catch.hpp>
......@@ -103,6 +104,10 @@ public:
};
TEST_CASE("ClassTimer", "[Timer]") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Measure runtime of a function without arguments") {
auto test = foo();
......
......@@ -11,6 +11,7 @@
#include <corsika/framework/stack/CombinedStack.hpp>
#include <corsika/framework/stack/SecondaryView.hpp>
#include <corsika/framework/stack/Stack.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <testTestStack.hpp> // for testing: simple stack. This is a
// test-build, and inluce file is obtained from CMAKE_CURRENT_SOURCE_DIR
......@@ -87,6 +88,9 @@ using StackTest = CombinedStack<TestStackData, TestStackData2, CombinedTestInter
TEST_CASE("Combined Stack", "[stack]") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// helper function for sum over stack data
auto sum = [](const StackTest& stack) {
double v = 0;
......@@ -283,6 +287,9 @@ using StackTest2 = CombinedStack<typename StackTest::stack_implementation_type,
TEST_CASE("Combined Stack - multi", "[stack]") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("create secondaries") {
StackTest2 s;
......@@ -378,6 +385,9 @@ using Particle2 = typename StackTest2::particle_type;
TEST_CASE("Combined Stack - secondary view") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("create secondaries via secondaryview") {
StackTest2 stack;
......
......@@ -21,6 +21,9 @@ static void handle_fpe(int /*signo*/) { gRESULT = 0; }
TEST_CASE("CorsikaFenv", "[fenv]") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Enable all exceptions") { feenableexcept(FE_ALL_EXCEPT); }
signal(SIGFPE, handle_fpe);
......
......@@ -20,6 +20,9 @@ using namespace corsika;
TEST_CASE("four vectors") {
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::info);
// this is just needed as a baseline
CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
......
......@@ -7,6 +7,7 @@
*/
#include <corsika/framework/analytics/FunctionTimer.hpp>
#include <corsika/framework/core/Logging.hpp>
#include <catch2/catch.hpp>
......@@ -30,6 +31,10 @@ public:
};
TEST_CASE("FunctionTimer", "[Timer]") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Measure runtime of a free function") {
auto test = corsika::FunctionTimer(testFunc);
......
......@@ -27,6 +27,10 @@ using namespace corsika::testing;
double constexpr absMargin = 1.0e-8;
TEST_CASE("transformations between CoordinateSystems") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
......
......@@ -21,6 +21,9 @@ double constexpr absMargin = 1.0e-8;
TEST_CASE("Helix class") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
const CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
Point r0(rootCS, {0_m, 0_m, 0_m});
......
......@@ -38,7 +38,8 @@ struct DummyProcess {
TEST_CASE("InteractionCounter", "[process]") {
logging::set_level(logging::level::debug);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::info);
DummyProcess d;
InteractionCounter countedProcess(d);
......@@ -100,6 +101,9 @@ TEST_CASE("InteractionCounter", "[process]") {
TEST_CASE("InteractionCounterOutput", "[output validation]") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
auto file = GENERATE(as<std::string>{}, "testInteractionCounter_file1",
"testInteractionCounter_file2");
......
......@@ -13,12 +13,15 @@
using namespace corsika;
TEST_CASE("Logging", "[Logging]") {
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("top level functions using default corsika logger") {
logging::info("This is an info message!");
logging::warn("This is a warning message!");
logging::debug("This is a debug message!");
logging::error("This is an error message!");
logging::critical("This is a critical error message!");
logging::info("(1) This is an info message!");
logging::warn("(1) This is a warning message!");
logging::debug("(1) This is a debug message!");
logging::error("(1) This is an error message!");
logging::critical("(1) This is a critical error message!");
}
SECTION("create a specific logger") {
......@@ -30,21 +33,21 @@ TEST_CASE("Logging", "[Logging]") {
logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// and make sure we can log with this created object
logger->info("This is an info message!");
logger->warn("This is a warning message!");
logger->debug("This is a debug message!");
logger->error("This is an error message!");
logger->critical("This is a critical error message!");
logger->info("(2) This is an info message!");
logger->warn("(2) This is a warning message!");
logger->debug("(2) This is a debug message!");
logger->error("(2) This is an error message!");
logger->critical("(2) This is a critical error message!");
// get a reference to the logger using Get
auto other = get_logger("loggerA");