From ae81c13aff136adabf91bf63c1ac8cb2875198ce Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Wed, 30 Jan 2019 12:07:44 +0100
Subject: [PATCH] added nuclearstackextension
---
Stack/NuclearStackExtension/CMakeLists.txt | 49 ++++
.../NuclearStackExtension.h | 235 ++++++++++++++++++
.../testNuclearStackExtension.cc | 192 ++++++++++++++
3 files changed, 476 insertions(+)
create mode 100644 Stack/NuclearStackExtension/CMakeLists.txt
create mode 100644 Stack/NuclearStackExtension/NuclearStackExtension.h
create mode 100644 Stack/NuclearStackExtension/testNuclearStackExtension.cc
diff --git a/Stack/NuclearStackExtension/CMakeLists.txt b/Stack/NuclearStackExtension/CMakeLists.txt
new file mode 100644
index 00000000..f2f34cfc
--- /dev/null
+++ b/Stack/NuclearStackExtension/CMakeLists.txt
@@ -0,0 +1,49 @@
+
+set (NuclearStackExtension_HEADERS NuclearStackExtension.h)
+set (NuclearStackExtension_NAMESPACE corsika/stack/nuclear_extension)
+
+add_library (NuclearStackExtension INTERFACE)
+
+CORSIKA_COPY_HEADERS_TO_NAMESPACE (NuclearStackExtension ${NuclearStackExtension_NAMESPACE} ${NuclearStackExtension_HEADERS})
+
+target_link_libraries (
+ NuclearStackExtension
+ INTERFACE
+ CORSIKAstackinterface
+ CORSIKAunits
+ CORSIKAparticles
+ CORSIKAgeometry
+ )
+
+target_include_directories (
+ NuclearStackExtension
+ INTERFACE
+ $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+ $<INSTALL_INTERFACE:include>
+ )
+
+install (
+ FILES
+ ${NuclearStackExtension_HEADERS}
+ DESTINATION
+ include/${NuclearStackExtension_NAMESPACE}
+ )
+
+# ----------------
+# code unit testing
+add_executable (
+ testNuclearStackExtension
+ testNuclearStackExtension.cc
+ )
+
+target_link_libraries (
+ testNuclearStackExtension
+ SuperStupidStack
+ NuclearStackExtension
+ # CORSIKAutls
+ ProcessStackInspector
+ CORSIKAgeometry
+ CORSIKAunits
+ CORSIKAthirdparty # for catch2
+ )
+CORSIKA_ADD_TEST(testNuclearStackExtension)
diff --git a/Stack/NuclearStackExtension/NuclearStackExtension.h b/Stack/NuclearStackExtension/NuclearStackExtension.h
new file mode 100644
index 00000000..1d5955fd
--- /dev/null
+++ b/Stack/NuclearStackExtension/NuclearStackExtension.h
@@ -0,0 +1,235 @@
+
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_stack_nuclearstackextension_h_
+#define _include_stack_nuclearstackextension_h_
+
+#include <corsika/stack/super_stupid/SuperStupidStack.h>
+
+#include <vector>
+
+namespace corsika::stack {
+
+ namespace nuclear_extension {
+
+ using corsika::stack::super_stupid::MomentumVector;
+
+ template <typename StackIteratorInterface>
+ class ParticleInterface
+ : public corsika::stack::super_stupid::ParticleInterface<StackIteratorInterface> {
+
+ protected:
+ // using corsika::stack::ParticleBase<StackIteratorInterface>::GetStack;
+ //using corsika::stack::super_stupid::ParticleInterface<StackIteratorInterface>::GetStackData;
+ using corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData;
+ using corsika::stack::ParticleBase<StackIteratorInterface>::GetIndex;
+
+ public:
+ void SetParticleData(const corsika::particles::Code vDataPID,
+ const corsika::units::si::HEPEnergyType vDataE,
+ const corsika::stack::super_stupid::MomentumVector& vMomentum,
+ const corsika::geometry::Point& vPosition,
+ const corsika::units::si::TimeType vTime,
+ const int vA = 0,
+ const int vZ = 0) {
+ if (vDataPID == corsika::particles::Code::Nucleus) {
+ if (vA==0 || vZ==0) {
+ std::ostringstream err;
+ err << "NuclearStackExtension: no A and Z specified for new Nucleus!";
+ throw std::runtime_error(err.str());
+ }
+ SetNuclearRef(corsika::stack::ParticleBase<StackIteratorInterface>::GetStackData().GetNucleusNextRef()); // store this nucleus data ref
+ SetNuclearA(vA);
+ SetNuclearZ(vZ);
+ } else {
+ SetNuclearRef(-1); // this is not a nucleus
+ }
+ corsika::stack::super_stupid::ParticleInterface<StackIteratorInterface>::SetParticleData(vDataPID,
+ vDataE,
+ vMomentum,
+ vPosition,
+ vTime);
+ }
+
+ void SetParticleData(ParticleInterface<StackIteratorInterface>& parent,
+ const corsika::particles::Code vDataPID,
+ const corsika::units::si::HEPEnergyType vDataE,
+ const corsika::stack::super_stupid::MomentumVector& vMomentum,
+ const corsika::geometry::Point& vPosition,
+ const corsika::units::si::TimeType vTime,
+ const int vA = 0,
+ const int vZ = 0) {
+ SetParticleData(vDataPID,
+ vDataE,
+ vMomentum,
+ vPosition,
+ vTime,
+ vA,
+ vZ);
+ }
+
+ /**
+ * @name individual setters
+ * @{
+ */
+ void SetNuclearA(const int vA) { GetStackData().SetNuclearA(GetIndex(), vA); }
+ void SetNuclearZ(const int vZ) { GetStackData().SetNuclearZ(GetIndex(), vZ); }
+ /// @}
+
+ /**
+ * @name individual getters
+ * @{
+ */
+ int GetNuclearA() const { return GetStackData().GetNuclearA(GetIndex()); }
+ int GetNuclearZ() const { return GetStackData().GetNuclearZ(GetIndex()); }
+ /// @}
+
+ protected:
+ void SetNuclearRef(const int vR) { GetStackData().SetNuclearRef(GetIndex(), vR); }
+ };
+
+
+ /**
+ * Memory implementation of the most simple (stupid) particle stack object.
+ */
+
+ class NuclearStackExtensionImpl
+ : public corsika::stack::super_stupid::SuperStupidStackImpl {
+
+ public:
+ void Init() { corsika::stack::super_stupid::SuperStupidStackImpl::Init(); }
+
+ void Clear() {
+ corsika::stack::super_stupid::SuperStupidStackImpl::Clear();
+ fNuclearRef.clear();
+ fNuclearA.clear();
+ fNuclearZ.clear();
+ }
+
+ int GetSize() const { return fNuclearRef.size(); }
+ int GetCapacity() const { return fNuclearRef.size(); }
+
+ void SetNuclearA(const int i, const int vA) { fNuclearA[GetNucleusRef(i)] = vA; }
+ void SetNuclearZ(const int i, const int vZ) { fNuclearZ[GetNucleusRef(i)] = vZ; }
+ void SetNuclearRef(const int i, const int v) { fNuclearRef[i] = v; }
+
+ int GetNuclearA(const int i) const { return fNuclearA[GetNucleusRef(i)]; }
+ int GetNuclearZ(const int i) const { return fNuclearZ[GetNucleusRef(i)]; }
+ // this function will create new storage for Nuclear Properties, and return the reference to it
+ int GetNucleusNextRef() { fNuclearA.push_back(0); fNuclearZ.push_back(0); return fNuclearA.size()-1; }
+
+ int GetNucleusRef(const int i) const {
+ if (fNuclearRef[i]>=0)
+ return fNuclearRef[i];
+ std::ostringstream err;
+ err << "NuclearStackExtension: no nucleus at ref=" << i;
+ throw std::runtime_error(err.str());
+ }
+
+ /**
+ * Function to copy particle at location i1 in stack to i2
+ */
+ void Copy(const int i1, const int i2) {
+ corsika::stack::super_stupid::SuperStupidStackImpl::Copy(i1, i2);
+ const int ref1 = fNuclearRef[i1];
+ const int ref2 = fNuclearRef[i2];
+ if (ref2<0) {
+ if (ref1>=0) {
+ // i1 is nucleus, i2 is not
+ fNuclearRef[i2] = GetNucleusNextRef();
+ fNuclearA[fNuclearRef[i2]] = fNuclearA[ref1];
+ fNuclearZ[fNuclearRef[i2]] = fNuclearZ[ref1];
+ } else {
+ // neither i1 nor i2 are nuclei
+ }
+ } else {
+ if (ref1>=0) {
+ // both are nuclei, i2 is overwritten with nucleus i1
+ fNuclearA[ref2] = fNuclearA[ref1];
+ fNuclearZ[ref2] = fNuclearZ[ref1];
+ } else {
+ // i2 is overwritten with non-nucleus i1
+ fNuclearA.erase(fNuclearA.begin() + ref2);
+ fNuclearZ.erase(fNuclearZ.begin() + ref2);
+ const int n = fNuclearRef.size();
+ for (int i=0; i<n; ++i) {
+ if (fNuclearRef[i]>=ref2) {
+ fNuclearRef[i] -= 1;
+ }
+ }
+ }
+ }
+ }
+
+ /**
+ * Function to copy particle at location i2 in stack to i1
+ */
+ void Swap(const int i1, const int i2) {
+ corsika::stack::super_stupid::SuperStupidStackImpl::Swap(i1, i2);
+ const int ref1 = fNuclearRef[i1];
+ const int ref2 = fNuclearRef[i2];
+ if (ref2<0) {
+ if (ref1>=0) {
+ // i1 is nucleus, i2 is not
+ std::swap(fNuclearRef[i2], fNuclearRef[i1]);
+ } else {
+ // neither i1 nor i2 are nuclei
+ }
+ } else {
+ if (ref1>=0) {
+ // both are nuclei, i2 is overwritten with nucleus i1
+ std::swap(fNuclearRef[i2], fNuclearRef[i1]);
+ } else {
+ // i2 is overwritten with non-nucleus i1
+ std::swap(fNuclearRef[i2], fNuclearRef[i1]);
+ }
+ }
+ }
+
+ protected:
+ void IncrementSize() {
+ corsika::stack::super_stupid::SuperStupidStackImpl::IncrementSize();
+ fNuclearRef.push_back(-1);
+ }
+
+ void DecrementSize() {
+ corsika::stack::super_stupid::SuperStupidStackImpl::DecrementSize();
+ if (fNuclearRef.size() > 0) {
+ const int ref = fNuclearRef.back();
+ fNuclearRef.pop_back();
+ if (ref>=0) {
+ fNuclearA.erase(fNuclearA.begin() + ref);
+ fNuclearZ.erase(fNuclearZ.begin() + ref);
+ const int n = fNuclearRef.size();
+ for (int i=0; i<n; ++i) {
+ if (fNuclearRef[i] >= ref) {
+ fNuclearRef[i] -= 1;
+ }
+ }
+ }
+ }
+ }
+
+ private:
+ /// the actual memory to store particle data
+
+ std::vector<int> fNuclearRef;
+ std::vector<int> fNuclearA;
+ std::vector<int> fNuclearZ;
+
+ }; // end class SuperStupidStackImpl
+
+ typedef Stack<NuclearStackExtensionImpl, ParticleInterface> NuclearStackExtension;
+
+ } // namespace nuclear_extension
+} // namespace corsika::stack
+
+#endif
diff --git a/Stack/NuclearStackExtension/testNuclearStackExtension.cc b/Stack/NuclearStackExtension/testNuclearStackExtension.cc
new file mode 100644
index 00000000..f8fac4bd
--- /dev/null
+++ b/Stack/NuclearStackExtension/testNuclearStackExtension.cc
@@ -0,0 +1,192 @@
+
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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/geometry/RootCoordinateSystem.h>
+#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
+#include <corsika/units/PhysicalUnits.h>
+
+using namespace corsika::geometry;
+using namespace corsika::units::si;
+
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+ // cpp file
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::stack::nuclear_extension;
+
+#include <iostream>
+using namespace std;
+
+TEST_CASE("NuclearStackExtension", "[stack]") {
+
+ geometry::CoordinateSystem& dummyCS =
+ geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+ SECTION("write non nucleus") {
+ NuclearStackExtension s;
+ s.AddParticle(particles::Code::Electron, 1.5_GeV,
+ 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") {
+ NuclearStackExtension s;
+ s.AddParticle(particles::Code::Nucleus, 1.5_GeV,
+ MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 10, 10);
+ REQUIRE(s.GetSize() == 1);
+ }
+
+ SECTION("write invalid nucleus") {
+ NuclearStackExtension s;
+ REQUIRE_THROWS(s.AddParticle(particles::Code::Nucleus, 1.5_GeV,
+ MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, 0, 0));
+ }
+
+ SECTION("read non nucleus") {
+ NuclearStackExtension s;
+ s.AddParticle(particles::Code::Electron, 1.5_GeV,
+ 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") {
+ NuclearStackExtension s;
+ s.AddParticle(particles::Code::Nucleus, 1.5_GeV,
+ 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.GetNuclearA() == 10);
+ REQUIRE(pout.GetNuclearZ() == 9);
+ }
+
+ SECTION("read invalid nucleus") {
+ NuclearStackExtension s;
+ s.AddParticle(particles::Code::Electron, 1.5_GeV,
+ 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") {
+
+ NuclearStackExtension 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(particles::Code::Nucleus, 1.5_GeV,
+ MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, i, i/2);
+ } else {
+ s.AddParticle(particles::Code::Electron, 1.5_GeV,
+ MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s);
+ }
+ }
+
+ REQUIRE(s.GetSize() == 99);
+
+ for (int i = 0; i < 99; ++i) s.GetNextParticle().Delete();
+
+ REQUIRE(s.GetSize() == 0);
+ }
+
+
+ SECTION("stack operations") {
+
+ NuclearStackExtension 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(particles::Code::Nucleus, i*15_GeV,
+ MomentumVector(dummyCS, {1_GeV, 1_GeV, 1_GeV}),
+ Point(dummyCS, {1 * meter, 1 * meter, 1 * meter}), 100_s, i, i/2);
+ } else {
+ s.AddParticle(particles::Code::Electron, i*1.5_GeV,
+ 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);
+ 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);
+ }
+
+ // 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.Copy(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);
+
+ }
+
+}
--
GitLab