IAP GITLAB

Commit f16cd74d authored by Ralf Ulrich's avatar Ralf Ulrich

style

parent d191c4c4
Pipeline #3133 canceled with stages
/*
* (c) Copyright 2018 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 <corsika/stack/dummy/DummyStack.h>
using namespace corsika;
using namespace corsika::stack;
#include <catch2/catch.hpp>
#include <tuple>
TEST_CASE("DummyStack", "[stack]") {
using TestStack = dummy::DummyStack;
dummy::NoData noData;
SECTION("write node") {
TestStack s;
s.AddParticle(std::tuple<dummy::NoData>{noData});
CHECK(s.getEntries() == 1);
}
SECTION("stack fill and cleanup") {
TestStack s;
// add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
for (int i = 0; i < 99; ++i) { s.AddParticle(std::tuple<dummy::NoData>{noData}); }
CHECK(s.getEntries() == 99);
for (int i = 0; i < 99; ++i) s.GetNextParticle().Delete();
CHECK(s.getEntries() == 0);
}
}
set (GeometryNodeStackExtension_HEADERS GeometryNodeStackExtension.h)
set (GeometryNodeStackExtension_NAMESPACE corsika/stack/node)
add_library (GeometryNodeStackExtension INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (GeometryNodeStackExtension ${GeometryNodeStackExtension_NAMESPACE} ${GeometryNodeStackExtension_HEADERS})
target_link_libraries (
GeometryNodeStackExtension
INTERFACE
CORSIKAstackinterface
CORSIKAunits
CORSIKAparticles
CORSIKAgeometry
)
target_include_directories (
GeometryNodeStackExtension
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include>
)
install (
FILES
${GeometryNodeStackExtension_HEADERS}
DESTINATION
include/${GeometryNodeStackExtension_NAMESPACE}
)
# ----------------
# code unit testing
CORSIKA_ADD_TEST(testGeometryNodeStackExtension)
target_link_libraries (
testGeometryNodeStackExtension
GeometryNodeStackExtension
CORSIKAtesting
)
/*
* (c) Copyright 2018 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/logging/Logging.h>
#include <corsika/stack/Stack.h>
#include <tuple>
#include <utility>
#include <vector>
namespace corsika::stack::node {
/**
* @class GeometryDataInterface
*
* 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
*/
template <typename T, typename TEnvType>
class GeometryDataInterface : public T {
protected:
using T::GetStack;
using T::GetStackData;
public:
using T::GetIndex;
using BaseNodeType = typename TEnvType::BaseNodeType;
public:
// default version for particle-creation from input data
void SetParticleData(const std::tuple<BaseNodeType const*> v) {
SetNode(std::get<0>(v));
}
void SetParticleData(GeometryDataInterface& parent,
const std::tuple<BaseNodeType const*>) {
SetNode(parent.GetNode()); // copy Node from parent particle!
}
void SetParticleData() { SetNode(nullptr); }
void SetParticleData(GeometryDataInterface& parent) {
SetNode(parent.GetNode()); // copy Node from parent particle!
}
std::string as_string() const { return fmt::format("node={}", fmt::ptr(GetNode())); }
void SetNode(BaseNodeType const* v) { GetStackData().SetNode(GetIndex(), v); }
BaseNodeType const* GetNode() const { return GetStackData().GetNode(GetIndex()); }
};
// definition of stack-data object to store geometry information
template <typename TEnvType>
/**
* @class GeometryData
*
* definition of stack-data object to store geometry information
*/
class GeometryData {
public:
using BaseNodeType = typename TEnvType::BaseNodeType;
// these functions are needed for the Stack interface
void Clear() { fNode.clear(); }
unsigned int GetSize() const { return fNode.size(); }
unsigned int GetCapacity() const { return fNode.size(); }
void Copy(const int i1, const int i2) { fNode[i2] = fNode[i1]; }
void Swap(const int i1, const int i2) { std::swap(fNode[i1], fNode[i2]); }
// custom data access function
void SetNode(const int i, BaseNodeType const* v) { fNode[i] = v; }
BaseNodeType const* GetNode(const int i) const { return fNode[i]; }
// these functions are also needed by the Stack interface
void IncrementSize() { fNode.push_back(nullptr); }
void DecrementSize() {
if (fNode.size() > 0) { fNode.pop_back(); }
}
// custom private data section
private:
std::vector<const BaseNodeType*> fNode;
};
template <typename T, typename TEnv>
struct MakeGeometryDataInterface {
typedef GeometryDataInterface<T, TEnv> type;
};
} // namespace corsika::stack::node
/*
* (c) Copyright 2018 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 <corsika/stack/CombinedStack.h>
#include <corsika/stack/dummy/DummyStack.h>
#include <corsika/stack/node/GeometryNodeStackExtension.h>
using namespace corsika;
using namespace corsika::stack;
#include <catch2/catch.hpp>
#include <iostream>
using namespace std;
// this is our dummy environment, it only knows its trivial BaseNodeType
class DummyEnv {
public:
typedef int BaseNodeType;
};
// the GeometryNode stack needs to know the type of geometry-nodes from the DummyEnv:
template <typename TStackIter>
using DummyGeometryDataInterface =
typename corsika::stack::node::MakeGeometryDataInterface<TStackIter, DummyEnv>::type;
// combine dummy stack with geometry information for tracking
template <typename TStackIter>
using StackWithGeometryInterface =
corsika::stack::CombinedParticleInterface<dummy::DummyStack::MPIType,
DummyGeometryDataInterface, TStackIter>;
using TestStack =
corsika::stack::CombinedStack<typename stack::dummy::DummyStack::StackImpl,
stack::node::GeometryData<DummyEnv>,
StackWithGeometryInterface>;
TEST_CASE("GeometryNodeStackExtension", "[stack]") {
dummy::NoData noData;
SECTION("write node") {
const int data = 5;
TestStack s;
s.AddParticle(std::make_tuple(noData), std::tuple<const int*>{&data});
CHECK(s.getEntries() == 1);
}
SECTION("write/read node") {
const int data = 15;
TestStack s;
auto p = s.AddParticle(std::make_tuple(noData));
p.SetNode(&data);
CHECK(s.getEntries() == 1);
const auto pout = s.GetNextParticle();
CHECK(*(pout.GetNode()) == 15);
}
SECTION("stack fill and cleanup") {
const int data = 16;
TestStack s;
// add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
for (int i = 0; i < 99; ++i) {
auto p = s.AddParticle(std::tuple<dummy::NoData>{noData});
p.SetNode(&data);
}
CHECK(s.getEntries() == 99);
double v = 0;
for (int i = 0; i < 99; ++i) {
auto p = s.GetNextParticle();
v += *(p.GetNode());
p.Delete();
}
CHECK(v == 99 * data);
CHECK(s.getEntries() == 0);
}
}
set (
HISTORY_HEADERS
EventType.hpp
Event.hpp
HistorySecondaryProducer.hpp
HistoryObservationPlane.hpp
HistoryStackExtension.hpp
SecondaryParticle.hpp
)
set (
HISTORY_NAMESPACE
corsika/stack/history
)
if (WITH_HISTORY)
set (
HISTORY_SOURCES
HistoryObservationPlane.cpp
)
add_library (CORSIKAhistory STATIC ${HISTORY_SOURCES})
set (IS_INTERFACE "")
else (WITH_HISTORY)
add_library (CORSIKAhistory INTERFACE)
set (IS_INTERFACE "INTERFACE")
endif (WITH_HISTORY)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAhistory ${HISTORY_NAMESPACE} ${HISTORY_HEADERS})
target_link_libraries (
CORSIKAhistory
${IS_INTERFACE}
CORSIKAparticles
CORSIKAgeometry
CORSIKAprocesssequence # for HistoryObservationPlane
CORSIKAsetup # for HistoryObservationPlane
CORSIKAlogging
SuperStupidStack
NuclearStackExtension # for testHistoryView.cc
C8::ext::boost # for HistoryObservationPlane
)
target_include_directories (
CORSIKAhistory
INTERFACE
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include/include>
)
install (
FILES ${HISTORY_HEADERS}
DESTINATION include/${HISTORY_NAMESPACE}
)
# ----------------
# code unit testing
CORSIKA_ADD_TEST(testHistoryStack)
target_link_libraries (
testHistoryStack
CORSIKAhistory
CORSIKAtesting
DummyStack
)
CORSIKA_ADD_TEST(testHistoryView)
target_link_libraries (
testHistoryView
CORSIKAhistory
CORSIKAtesting
)
/*
* (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.
*/
#pragma once
#include <corsika/logging/Logging.h>
#include <corsika/particles/ParticleProperties.h>
#include <corsika/stack/history/EventType.hpp>
#include <corsika/stack/history/SecondaryParticle.hpp>
#include <memory>
#include <optional>
#include <vector>
namespace corsika::history {
class Event;
using EventPtr = std::shared_ptr<history::Event>;
class Event {
size_t projectile_index_ = 0; //!< index of projectile on stack
std::vector<SecondaryParticle> secondaries_;
EventPtr parent_event_;
EventType type_ = EventType::Uninitialized;
std::optional<corsika::particles::Code> targetCode_;
public:
Event() = default;
void setParentEvent(EventPtr const& evt) { parent_event_ = evt; }
bool hasParentEvent() const { return bool(parent_event_); }
EventPtr& parentEvent() { return parent_event_; }
EventPtr const& parentEvent() const { return parent_event_; }
void setProjectileIndex(size_t i) { projectile_index_ = i; }
size_t projectileIndex() const { return projectile_index_; }
size_t addSecondary(units::si::HEPEnergyType energy,
geometry::Vector<units::si::hepmomentum_d> const& momentum,
particles::Code pid) {
secondaries_.emplace_back(energy, momentum, pid);
return secondaries_.size() - 1;
}
std::vector<SecondaryParticle> const& secondaries() const { return secondaries_; }
void setTargetCode(const particles::Code t) { targetCode_ = t; }
std::string as_string() const {
return fmt::format("hasParent={}, projIndex={}, Nsec={}", hasParentEvent(),
projectile_index_, secondaries_.size());
}
EventType eventType() const { return type_; }
void setEventType(EventType t) { type_ = t; }
};
} // namespace corsika::history
/*
* (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.
*/
#pragma once
namespace corsika::history {
enum class EventType { Uninitialized, Interaction, Decay };
}
/*
* (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 <corsika/logging/Logging.h>
#include <corsika/stack/history/HistoryObservationPlane.hpp>
#include <boost/histogram/ostream.hpp>
#include <iomanip>
#include <iostream>
using namespace corsika::units::si;
using namespace corsika::history;
using namespace corsika;
HistoryObservationPlane::HistoryObservationPlane(setup::Stack const& stack,
geometry::Plane const& obsPlane,
bool deleteOnHit)
: stack_{stack}
, plane_{obsPlane}
, deleteOnHit_{deleteOnHit} {}
corsika::process::EProcessReturn HistoryObservationPlane::DoContinuous(
setup::Stack::ParticleType const& particle, setup::Trajectory const& trajectory) {
TimeType const timeOfIntersection =
(plane_.GetCenter() - trajectory.GetLine().GetR0()).dot(plane_.GetNormal()) /
trajectory.GetLine().GetV0().dot(plane_.GetNormal());
if (timeOfIntersection < TimeType::zero()) { return process::EProcessReturn::eOk; }
if (plane_.IsAbove(trajectory.GetLine().GetR0()) ==
plane_.IsAbove(trajectory.GetPosition(1))) {
return process::EProcessReturn::eOk;
}
C8LOG_DEBUG(fmt::format("HistoryObservationPlane: Particle detected: pid={}",
particle.GetPID()));
auto const pid = particle.GetPID();
if (particles::IsMuon(pid)) { fillHistoryHistogram(particle); }
if (deleteOnHit_) {
return process::EProcessReturn::eParticleAbsorbed;
} else {
return process::EProcessReturn::eOk;
}
}
LengthType HistoryObservationPlane::MaxStepLength(setup::Stack::ParticleType const&,
setup::Trajectory const& trajectory) {
TimeType const timeOfIntersection =
(plane_.GetCenter() - trajectory.GetLine().GetR0()).dot(plane_.GetNormal()) /
trajectory.GetLine().GetV0().dot(plane_.GetNormal());
if (timeOfIntersection < TimeType::zero()) {
return std::numeric_limits<double>::infinity() * 1_m;
}
auto const pointOfIntersection = trajectory.GetLine().GetPosition(timeOfIntersection);
return (trajectory.GetLine().GetR0() - pointOfIntersection).norm() * 1.0001;
}
void HistoryObservationPlane::fillHistoryHistogram(
setup::Stack::ParticleType const& muon) {
double const muon_energy = muon.GetEnergy() / 1_GeV;
int genctr{0};
Event const* event = muon.GetEvent().get();
while (event) {
auto const projectile = stack_.cfirst() + event->projectileIndex();
if (event->eventType() == EventType::Interaction) {
genctr++;
double const projEnergy = projectile.GetEnergy() / 1_GeV;
int const pdg = static_cast<int>(particles::GetPDG(projectile.GetPID()));
histogram_(muon_energy, projEnergy, pdg);
}
event = event->parentEvent().get(); // projectile.GetEvent().get();
}
}
/*
* (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.
*/
#pragma once
#include <corsika/geometry/Plane.h>
#include <corsika/process/ContinuousProcess.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
#include <boost/histogram.hpp>
#include <functional>
namespace corsika::history {
namespace detail {
inline auto hist_factory() {
namespace bh = boost::histogram;
namespace bha = bh::axis;
auto h = bh::make_histogram(
bha::regular<float, bha::transform::log>{11 * 5, 1e0, 1e11, "muon energy"},
bha::regular<float, bha::transform::log>{11 * 5, 1e0, 1e11,
"projectile energy"},
bha::category<int, bh::use_default, bha::option::growth_t>{
{211, -211, 2212, -2212}, "projectile PDG"});
return h;
}
} // namespace detail
class HistoryObservationPlane
: public corsika::process::ContinuousProcess<HistoryObservationPlane> {
public:
HistoryObservationPlane(setup::Stack const&, geometry::Plane const&, bool = true);
corsika::units::si::LengthType MaxStepLength(
corsika::setup::Stack::ParticleType const&,
corsika::setup::Trajectory const& vTrajectory);
corsika::process::EProcessReturn DoContinuous(
corsika::setup::Stack::ParticleType const& vParticle,
corsika::setup::Trajectory const& vTrajectory);
auto const& histogram() const { return histogram_; }
private:
void fillHistoryHistogram(setup::Stack::ParticleType const&);
setup::Stack const& stack_;
geometry::Plane const plane_;
bool const deleteOnHit_;
decltype(detail::hist_factory()) histogram_ = detail::hist_factory();
};
} // namespace corsika::history
/*
* (c) Copyright 2018 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/stack/SecondaryView.h>
#include <corsika/stack/history/Event.hpp>
#include <corsika/logging/Logging.h>
#include <boost/type_index.hpp>
#include <memory>
#include <type_traits>
#include <utility>
namespace corsika::history {
//! mix-in class for SecondaryView that fills secondaries into an \class Event
template <class T1, template <class> class T2>
class HistorySecondaryProducer {
public:
EventPtr event_;
static bool constexpr has_event{true};
public:
template <typename Particle>
HistorySecondaryProducer(Particle const& p)