IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 02cc30ee authored by ralfulrich's avatar ralfulrich
Browse files

various small ctest improvemnts, started to add testHistory

parent 3da64964
No related branches found
No related tags found
1 merge request!254History
#include <memory>
#include <utility>
#include <variant>
using ParticleCode = std::variant<corsika::particles::ParticleCode, std::pair<int, int>>;
class SecondaryParticle {
corsika::units::si::HEPEnergyType const energy_;
corsika::geometry::Vector<corsika::units::si::hepmomentum_d> const momentum_;
ParticleCode const pid_;
// what else...?
// - polarization
public:
SecondaryParticle(energy, momentum, pid)
: energy_{energy}
, momentum_{momentum}
, pid_{pid} {}
};
class Event; // forward declaration
using EventPtr = std::shared_ptr<Event>;
struct ProjectileParticle {
std::pair<EventPtr, size_t> origin;
corsika::geometry::Point position;
coriska::units::si::TimeType time;
corsika::geometry::Vector<corsika::units::si::hepmomentum_d> momentum;
corsika::units::si::HEPEnergyType energy;
// weight
};
struct Event {
corsika::particles::ParticleCode targetCode; // cannot be const...
ProjectileParticle projectile;
std::vector<SecondaryParticle> secondaries;
};
template <typename TStackView>
class EventBuilder {
TStackView& stackView_;
EventPtr event_;
public:
EventBuilder(TStackView& view)
: stackView_{stackView}
, event_{std::make_shared<Event>()} {}
template <typename... Args>
void AddSecondary(Args&&... args) {
auto const s = stackView.AddSecondary(std::forward<Args>(args), event);
event->secondaries.emplace_back(s.GetEnergy(), s.GetMomentum(), s.GetParticleID());
}
void SetProjectile() {}
void SetTarget(corsika::particles::ParticleCode targetCode) {
event_->targetCode = targetCode;
}
};
......@@ -287,5 +287,5 @@ TEST_CASE("SwitchProcess") {
Process1 p1(0);
Process2 p2(1);
switch_process::SwitchProcess s(p1, p2, 10_GeV);
REQUIRE(is_switch_process_v<decltype(s)>);
CHECK(is_switch_process_v<decltype(s)>);
}
......@@ -2,4 +2,4 @@ add_subdirectory (DummyStack)
add_subdirectory (SuperStupidStack)
add_subdirectory (NuclearStackExtension)
add_subdirectory (GeometryNodeStackExtension)
#add_subdirectory (History)
add_subdirectory (History)
......@@ -36,6 +36,7 @@ namespace corsika::stack::node {
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));
......@@ -55,6 +56,8 @@ namespace corsika::stack::node {
BaseNodeType const* GetNode() const { return GetStackData().GetNode(GetIndex()); }
};
// definition of stack-data object to store geometry information
template <typename TEnvType>
......
......@@ -28,7 +28,6 @@ public:
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<
......
......@@ -35,9 +35,6 @@ install (
CORSIKA_ADD_TEST(testHistory)
target_link_libraries (
testHistory
SuperStupidStack
CORSIKAhistory
CORSIKAparticles
CORSIKAunits
CORSIKAtesting
)
......@@ -8,29 +8,22 @@
#pragma once
// the basic particle data stack:
#include <corsika/stack/super_stupid/SuperStupidStack.h>
#include <corsika/stack/Stack.h>
// extension with nuclear data for Code::Nucleus
#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
// extension with nuclear data AND volume node ref
#include <corsika/setup/GeometryNodeStackExtension.h>
#include <memory>
#include <tuple>
#include <utility>
class Event {};
#include <vector>
#include <memory>
namespace corsika::history {
/**
* @class HistoryData
*
* definition of stack-data object to store history information
* this is vector with shared_ptr<Event>
* definition of stack-data object to store history information this
* is vector with shared_ptr<TEvent>, where TEvent is a free
* template parameter for customization.
*/
template <typename TEvent>
class HistoryData {
public:
......@@ -42,8 +35,8 @@ namespace corsika::history {
void Swap(const int i1, const int i2) { std::swap(fEvent[i1], fEvent[i2]); }
// custom data access function
void SetEvent(const int i, std::shared_ptr<Event> v) { fEvent[i] = v; }
std::shared_ptr<Event> GetEvent(const int i) const { return fEvent[i]; }
void SetEvent(const int i, std::shared_ptr<TEvent> v) { fEvent[i] = v; }
std::shared_ptr<TEvent> GetEvent(const int i) const { return fEvent[i]; }
// these functions are also needed by the Stack interface
void IncrementSize() { fEvent.push_back(nullptr); }
......@@ -53,7 +46,7 @@ namespace corsika::history {
// custom private data section
private:
std::vector<std::shared_ptr<Event>> fEvent;
std::vector<std::shared_ptr<TEvent>> fEvent;
};
/**
......@@ -64,62 +57,36 @@ namespace corsika::history {
// defintion of a stack-readout object, the iteractor dereference
// operator will deliver access to these function
*/
template <typename T>
template <typename T, typename TEvent>
class HistoryDataInterface : public T {
protected:
using T::GetStack;
using T::GetStackData;
public:
using T::GetIndex;
using T::GetStackData;
using T::SetParticleData;
public:
using T::GetIndex;
public:
// default version for particle-creation from input data
void SetParticleData(const std::tuple<Event const*> v) { SetEvent(std::get<0>(v)); }
void SetParticleData(const std::tuple<TEvent const*> v) { SetEvent(std::get<0>(v)); }
void SetParticleData(HistoryDataInterface& parent,
const std::tuple<std::shared_ptr<Event>>) {
const std::tuple<std::shared_ptr<TEvent>>) {
SetEvent(parent.GetEvent()); // copy Event from parent particle!
}
void SetParticleData() { SetEvent(nullptr); }
void SetParticleData(HistoryDataInterface& parent) {
SetEvent(parent.GetEvent()); // copy Event from parent particle!
}
void SetEvent(std::shared_ptr<Event> v) { GetStackData().SetEvent(GetIndex(), v); }
std::shared_ptr<Event> GetEvent() const {
void SetEvent(std::shared_ptr<TEvent> v) { GetStackData().SetEvent(GetIndex(), v); }
std::shared_ptr<TEvent> GetEvent() const {
return GetStackData().GetEvent(GetIndex());
}
};
namespace detail {
//
// this is an auxiliary help typedef, which I don't know how to put
// into NuclearStackExtension.h where it belongs...
template <typename StackIter>
using ExtendedParticleInterfaceType =
corsika::stack::nuclear_extension::NuclearParticleInterface<
corsika::stack::super_stupid::SuperStupidStack::PIType, StackIter>;
//
// the particle data stack with extra nuclear information:
using ParticleDataStack = corsika::stack::nuclear_extension::NuclearStackExtension<
corsika::stack::super_stupid::SuperStupidStack, ExtendedParticleInterfaceType>;
template <typename T>
using SetupHistoryDataInterface = HistoryDataInterface<T>;
// combine particle data stack with history information for tracking
template <typename StackIter>
using StackWithHistoryInterface =
corsika::stack::CombinedParticleInterface<ParticleDataStack::PIType,
SetupHistoryDataInterface, StackIter>;
using StackWithHistory =
corsika::stack::CombinedStack<typename ParticleDataStack::StackImpl, HistoryData,
StackWithHistoryInterface>;
} // namespace detail
template <typename InnerStack, template <typename> typename _PI>
using NuclearStackExtension =
Stack<NuclearStackExtensionImpl<typename InnerStack::StackImpl>, _PI>;
template <typename T, typename TEvent>
struct MakeHistoryDataInterface {
typedef HistoryDataInterface<T, TEvent> type;
};
} // namespace corsika::history
......@@ -6,236 +6,51 @@
* the license.
*/
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/history/HistoryStackExtension.h>
#include <corsika/stack/super_stupid/SuperStupidStack.h>
#include <corsika/stack/dummy/DummyStack.h>
#include <corsika/stack/CombinedStack.h>
using namespace corsika;
using namespace corsika::history;
using namespace corsika::stack::nuclear_extension;
using namespace corsika::geometry;
using namespace corsika::units::si;
using namespace corsika::stack;
#include <catch2/catch.hpp>
// this is an auxiliary help typedef, which I don't know how to put
// into NuclearStackExtension.h where it belongs...
template <typename StackIter>
using ExtendedParticleInterfaceType =
corsika::stack::nuclear_extension::NuclearParticleInterface<
corsika::stack::super_stupid::SuperStupidStack::template PIType, StackIter>;
using ExtStack = NuclearStackExtension<corsika::stack::super_stupid::SuperStupidStack,
ExtendedParticleInterfaceType>;
#include <iostream>
using namespace std;
TEST_CASE("HistoryStackExtension", "[stack]") {
// this is our dummy environment, it only knows its trivial BaseNodeType
class DummyEvent {
public:
int id;
};
geometry::CoordinateSystem& dummyCS =
geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
// the GeometryNode stack needs to know the type of geometry-nodes from the DummyEnv:
template <typename TStackIter>
using DummyHistoryDataInterface = typename history::MakeHistoryDataInterface<TStackIter, DummyEvent>::type;
SECTION("write non nucleus") {
HistoryStackExtension<corsika::stack::super_stupid::SuperStupidStack,
ExtendedParticleInterfaceType>
s;
s.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Electron, 1.5_GeV,
corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
REQUIRE(s.GetSize() == 1);
}
SECTION("write nucleus") {
HistoryStackExtension<corsika::stack::super_stupid::SuperStupidStack,
ExtendedParticleInterfaceType>
s;
s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, 1.5_GeV,
corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 10});
REQUIRE(s.GetSize() == 1);
}
// combine dummy stack with geometry information for tracking
template <typename TStackIter>
using StackWithHistoryInterface = corsika::stack::CombinedParticleInterface<
dummy::DummyStack::PIType,
DummyHistoryDataInterface, TStackIter>;
SECTION("write invalid nucleus") {
ExtStack s;
REQUIRE_THROWS(
s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, 1.5_GeV,
corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 0, 0}));
}
using TestStack = corsika::stack::CombinedStack<
typename stack::dummy::DummyStack::StackImpl,
history::HistoryData<DummyEvent>,
StackWithHistoryInterface>;
SECTION("read non nucleus") {
ExtStack s;
s.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Electron, 1.5_GeV,
corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
const auto pout = s.GetNextParticle();
REQUIRE(pout.GetPID() == particles::Code::Electron);
REQUIRE(pout.GetEnergy() == 1.5_GeV);
REQUIRE(pout.GetTime() == 100_s);
}
SECTION("read nucleus") {
ExtStack s;
s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, 1.5_GeV,
corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 9});
const auto pout = s.GetNextParticle();
REQUIRE(pout.GetPID() == particles::Code::Nucleus);
REQUIRE(pout.GetEnergy() == 1.5_GeV);
REQUIRE(pout.GetTime() == 100_s);
REQUIRE(pout.GetHistoryA() == 10);
REQUIRE(pout.GetNuclearZ() == 9);
}
SECTION("read invalid nucleus") {
ExtStack s;
s.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
particles::Code::Electron, 1.5_GeV,
corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
const auto pout = s.GetNextParticle();
REQUIRE_THROWS(pout.GetNuclearA());
REQUIRE_THROWS(pout.GetNuclearZ());
}
SECTION("stack fill and cleanup") {
TEST_CASE("HistoryStackExtension", "[stack]") {
ExtStack s;
// add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
for (int i = 0; i < 99; ++i) {
if ((i + 1) % 10 == 0) {
s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, 1.5_GeV,
corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, i, i / 2});
} else {
s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType>{
particles::Code::Electron, 1.5_GeV,
corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
}
}
const dummy::NoData noData;
SECTION("write event") {
REQUIRE(s.GetSize() == 99);
for (int i = 0; i < 99; ++i) s.GetNextParticle().Delete();
REQUIRE(s.GetSize() == 0);
TestStack s;
s.AddParticle(std::tuple<dummy::NoData>{noData});
REQUIRE(s.GetSize() == 1);
}
SECTION("stack operations") {
ExtStack s;
// add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
for (int i = 0; i < 99; ++i) {
if ((i + 1) % 10 == 0) {
s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType, unsigned short, unsigned short>{
particles::Code::Nucleus, i * 15_GeV,
corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, i, i / 2});
} else {
s.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point,
units::si::TimeType>{
particles::Code::Electron, i * 1.5_GeV,
corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
}
}
// copy
{
s.Copy(s.begin() + 9, s.begin() + 10); // nuclei to non-nuclei
const auto& p9 = s.cbegin() + 9;
const auto& p10 = s.cbegin() + 10;
REQUIRE(p9.GetPID() == particles::Code::Nucleus);
REQUIRE(p9.GetEnergy() == 9 * 15_GeV);
REQUIRE(p9.GetTime() == 100_s);
REQUIRE(p9.GetNuclearA() == 9);
REQUIRE(p9.GetNuclearZ() == 9 / 2);
REQUIRE(p10.GetPID() == particles::Code::Nucleus);
REQUIRE(p10.GetEnergy() == 9 * 15_GeV);
REQUIRE(p10.GetTime() == 100_s);
REQUIRE(p10.GetNuclearA() == 9);
REQUIRE(p10.GetNuclearZ() == 9 / 2);
}
// copy
{
s.Copy(s.begin() + 93, s.begin() + 9); // non-nuclei to nuclei
const auto& p93 = s.cbegin() + 93;
const auto& p9 = s.cbegin() + 9;
REQUIRE(p9.GetPID() == particles::Code::Electron);
REQUIRE(p9.GetEnergy() == 93 * 1.5_GeV);
REQUIRE(p9.GetTime() == 100_s);
REQUIRE(p93.GetPID() == particles::Code::Electron);
REQUIRE(p93.GetEnergy() == 93 * 1.5_GeV);
REQUIRE(p93.GetTime() == 100_s);
}
// swap
{
s.Swap(s.begin() + 11, s.begin() + 10);
const auto& p11 = s.cbegin() + 11; // now: nucleus
const auto& p10 = s.cbegin() + 10; // now: electron
REQUIRE(p11.GetPID() == particles::Code::Nucleus);
REQUIRE(p11.GetEnergy() == 9 * 15_GeV);
REQUIRE(p11.GetTime() == 100_s);
REQUIRE(p11.GetNuclearA() == 9);
REQUIRE(p11.GetNuclearZ() == 9 / 2);
REQUIRE(p10.GetPID() == particles::Code::Electron);
REQUIRE(p10.GetEnergy() == 11 * 1.5_GeV);
REQUIRE(p10.GetTime() == 100_s);
}
// swap two nuclei
{
s.Swap(s.begin() + 29, s.begin() + 59);
const auto& p29 = s.cbegin() + 29;
const auto& p59 = s.cbegin() + 59;
REQUIRE(p29.GetPID() == particles::Code::Nucleus);
REQUIRE(p29.GetEnergy() == 59 * 15_GeV);
REQUIRE(p29.GetTime() == 100_s);
REQUIRE(p29.GetNuclearA() == 59);
REQUIRE(p29.GetNuclearZ() == 59 / 2);
REQUIRE(p59.GetPID() == particles::Code::Nucleus);
REQUIRE(p59.GetEnergy() == 29 * 15_GeV);
REQUIRE(p59.GetTime() == 100_s);
REQUIRE(p59.GetNuclearA() == 29);
REQUIRE(p59.GetNuclearZ() == 29 / 2);
}
for (int i = 0; i < 99; ++i) s.DeleteLast();
REQUIRE(s.GetSize() == 0);
}
//REQUIRE(pout.GetPID() == particles::Code::Electron);
}
......@@ -37,14 +37,14 @@ TEST_CASE("SuperStupidStack", "[stack]") {
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
// read
REQUIRE(s.GetSize() == 1);
CHECK(s.GetSize() == 1);
auto pout = s.GetNextParticle();
REQUIRE(pout.GetPID() == particles::Code::Electron);
REQUIRE(pout.GetEnergy() == 1.5_GeV);
// REQUIRE(pout.GetMomentum() == stack::MomentumVector(dummyCS, {1_GeV,
// 1_GeV, 1_GeV})); REQUIRE(pout.GetPosition() == Point(dummyCS, {1 * meter, 1 *
CHECK(pout.GetPID() == particles::Code::Electron);
CHECK(pout.GetEnergy() == 1.5_GeV);
// CHECK(pout.GetMomentum() == stack::MomentumVector(dummyCS, {1_GeV,
// 1_GeV, 1_GeV})); CHECK(pout.GetPosition() == Point(dummyCS, {1 * meter, 1 *
// meter, 1 * meter}));
REQUIRE(pout.GetTime() == 100_s);
CHECK(pout.GetTime() == 100_s);
}
SECTION("write+delete") {
......@@ -59,10 +59,10 @@ TEST_CASE("SuperStupidStack", "[stack]") {
corsika::stack::MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s});
REQUIRE(s.GetSize() == 99);
CHECK(s.GetSize() == 99);
for (int i = 0; i < 99; ++i) s.GetNextParticle().Delete();
REQUIRE(s.GetSize() == 0);
CHECK(s.GetSize() == 0);
}
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment