IAP GITLAB

Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • AirShowerPhysics/corsika
  • rulrich/corsika
  • AAAlvesJr/corsika
  • Andre/corsika
  • arrabito/corsika
  • Nikos/corsika
  • olheiser73/corsika
  • AirShowerPhysics/papers/corsika
  • pranav/corsika
9 results
Show changes
Showing
with 0 additions and 1046 deletions
set (GEOMETRY_SOURCES CoordinateSystem.cc)
set (GEOMETRY_HEADERS Vector.h Point.h Sphere.h CoordinateSystem.h Helix.h)
add_library (CORSIKAgeometry STATIC ${GEOMETRY_SOURCES})
set_target_properties (CORSIKAgeometry PROPERTIES VERSION ${PROJECT_VERSION})
set_target_properties (CORSIKAgeometry PROPERTIES SOVERSION 1)
set_target_properties (CORSIKAgeometry PROPERTIES PUBLIC_HEADER "${GEOMETRY_HEADERS}")
# target dependencies on other libraries (also header only)
target_link_libraries (CORSIKAgeometry CORSIKAunits)
target_include_directories (CORSIKAgeometry PRIVATE ${EIGEN3_INCLUDE_DIR})
target_include_directories (CORSIKAgeometry INTERFACE ${EIGEN3_INCLUDE_DIR})
target_include_directories (CORSIKAgeometry INTERFACE $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Framework>
$<INSTALL_INTERFACE:include/Framework>
)
install (TARGETS CORSIKAgeometry
LIBRARY DESTINATION lib
ARCHIVE DESTINATION lib
PUBLIC_HEADER DESTINATION include/Geometry)
# code testing
add_executable (testGeometry testGeometry.cc)
target_link_libraries (testGeometry CORSIKAgeometry CORSIKAunits CORSIKAthirdparty) # for catch2
add_test (NAME testGeometry COMMAND testGeometry -o report.xml -r junit)
#include <Geometry/CoordinateSystem.h>
EigenTransform CoordinateSystem::getTransformation(CoordinateSystem const& c1, CoordinateSystem const& c2)
{
CoordinateSystem const* a{&c1};
CoordinateSystem const* b{&c2};
CoordinateSystem const* commonBase{nullptr};
while (a != b && b != nullptr)
{
a = &c1;
while (a != b && a != nullptr)
{
a = a->getReference();
}
if (a == b)
break;
b = b->getReference();
}
if (a == b && a != nullptr)
{
commonBase = a;
}
else
{
throw std::string("no connection between coordinate systems found!");
}
EigenTransform t = EigenTransform::Identity();
auto* p = &c1;
while (p != commonBase)
{
t = p->getTransform() * t;
p = p->getReference();
}
p = &c2;
while (p != commonBase)
{
t = p->getTransform().inverse(Eigen::TransformTraits::Isometry) * t;
p = p->getReference();
}
return t;
}
#ifndef _include_COORDINATESYSTEM_H_
#define _include_COORDINATESYSTEM_H_
#include <Geometry/QuantityVector.h>
#include <Units/PhysicalUnits.h>
#include <Eigen/Dense>
typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
typedef Eigen::Translation<double, 3> EigenTranslation;
class CoordinateSystem
{
CoordinateSystem const* reference = nullptr;
EigenTransform transf;
CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf) :
reference(&reference),
transf(transf)
{
}
public:
static EigenTransform getTransformation(CoordinateSystem const& c1, CoordinateSystem const& c2);
CoordinateSystem() : // for creating the root CS
transf(EigenTransform::Identity())
{
}
auto& operator=(const CoordinateSystem& pCS)
{
reference = pCS.reference;
transf = pCS.transf;
return *this;
}
auto translate(QuantityVector<phys::units::length_d> vector) const
{
EigenTransform const translation{EigenTranslation(vector.eVector)};
return CoordinateSystem(*this, translation);
}
auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const
{
EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
return CoordinateSystem(*this, rotation);
}
auto translateAndRotate(QuantityVector<phys::units::length_d> translation, QuantityVector<phys::units::length_d> axis, double angle)
{
EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) * EigenTranslation(translation.eVector)};
return CoordinateSystem(*this, transf);
}
auto const* getReference() const
{
return reference;
}
auto const& getTransform() const
{
return transf;
}
};
#endif
#ifndef _include_HELIX_H_
#define _include_HELIX_H_
#include <Geometry/Vector.h>
#include <Geometry/Point.h>
#include <Units/PhysicalUnits.h>
class Helix // TODO: inherit from to-be-implemented "Trajectory"
{
using SpeedVec = Vector<phys::units::speed_d>;
Point const r0;
phys::units::quantity<phys::units::frequency_d> const omegaC;
SpeedVec const vPar;
SpeedVec vPerp, uPerp;
public:
Helix(Point const pR0, phys::units::quantity<phys::units::frequency_d> pOmegaC,
SpeedVec const pvPar, SpeedVec const pvPerp) :
r0(pR0), omegaC(pOmegaC), vPar(pvPar), vPerp(pvPerp), uPerp(vPar.normalized().cross(vPerp))
{
}
Point getPosition(phys::units::quantity<phys::units::time_interval_d> t) const
{
return r0 + vPar * t + (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
}
};
#endif
#ifndef _include_POINT_H_
#define _include_POINT_H_
#include <Geometry/BaseVector.h>
#include <Geometry/QuantityVector.h>
#include <Geometry/Vector.h>
#include <Units/PhysicalUnits.h>
class Point : public BaseVector<phys::units::length_d>
{
using Length = phys::units::quantity<phys::units::length_d, double>;
public:
Point(CoordinateSystem const& pCS, QuantityVector<phys::units::length_d> pQVector) :
BaseVector<phys::units::length_d>(pCS, pQVector)
{
}
Point(CoordinateSystem const& cs, Length x, Length y, Length z) :
BaseVector<phys::units::length_d>(cs, {x, y, z})
{
}
auto getCoordinates() const
{
return BaseVector<phys::units::length_d>::qVector;
}
auto getCoordinates(CoordinateSystem const& pCS) const
{
if (&pCS == BaseVector<phys::units::length_d>::cs)
{
return BaseVector<phys::units::length_d>::qVector;
}
else
{
return QuantityVector<phys::units::length_d>(CoordinateSystem::getTransformation(*BaseVector<phys::units::length_d>::cs, pCS) * BaseVector<phys::units::length_d>::qVector.eVector);
}
}
void rebase(CoordinateSystem const& pCS)
{
BaseVector<phys::units::length_d>::qVector = getCoordinates(pCS);
BaseVector<phys::units::length_d>::cs = &pCS;
}
Point operator+(Vector<phys::units::length_d> const& pVec) const
{
return Point(*BaseVector<phys::units::length_d>::cs, getCoordinates() + pVec.getComponents(*BaseVector<phys::units::length_d>::cs));
}
Vector<phys::units::length_d> operator-(Point const& pB) const
{
auto& cs = *BaseVector<phys::units::length_d>::cs;
return Vector<phys::units::length_d>(cs, getCoordinates() - pB.getCoordinates(cs));
}
};
#endif
#ifndef _include_QUANTITYVECTOR_H_
#define _include_QUANTITYVECTOR_H_
#include <Units/PhysicalUnits.h>
#include <Eigen/Dense>
#include <iostream>
#include <utility>
template <typename dim>
class QuantityVector
{
protected:
using Quantity = phys::units::quantity<dim, double>;
//using QuantitySquared = decltype(std::declval<Quantity>() * std::declval<Quantity>());
public:
Eigen::Vector3d eVector;
typedef dim dimension;
QuantityVector(Quantity a, Quantity b, Quantity c) :
eVector{a.magnitude(), b.magnitude(), c.magnitude()}
{
}
QuantityVector(Eigen::Vector3d pBareVector) :
eVector(pBareVector)
{
}
auto operator[](size_t index) const
{
return Quantity(phys::units::detail::magnitude_tag, eVector[index]);
}
auto norm() const
{
return Quantity(phys::units::detail::magnitude_tag, eVector.norm());
}
auto squaredNorm() const
{
using QuantitySquared = decltype(std::declval<Quantity>() * std::declval<Quantity>());
return QuantitySquared(phys::units::detail::magnitude_tag, eVector.squaredNorm());
}
auto operator+(QuantityVector<dim> const& pQVec) const
{
return QuantityVector<dim>(eVector + pQVec.eVector);
}
auto operator-(QuantityVector<dim> const& pQVec) const
{
return QuantityVector<dim>(eVector - pQVec.eVector);
}
template <typename ScalarDim>
auto operator*(phys::units::quantity<ScalarDim, double> const p) const
{
using ResQuantity = phys::units::detail::Product<ScalarDim, dim, double, double>;
if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not a "Quantity" anymore
{
return QuantityVector<phys::units::dimensionless_d>(eVector * p.magnitude());
}
else
{
return QuantityVector<typename ResQuantity::dimension_type>(eVector * p.magnitude());
}
}
template <typename ScalarDim>
auto operator/(phys::units::quantity<ScalarDim, double> const p) const
{
return (*this) * (1 / p);
}
auto operator*(double const p) const
{
return QuantityVector<dim>(eVector * p);
}
auto operator/(double const p) const
{
return QuantityVector<dim>(eVector / p);
}
auto& operator/=(double const p) const
{
eVector /= p;
return *this;
}
auto& operator*=(double const p)
{
eVector *= p;
return *this;
}
auto& operator+=(QuantityVector<dim> const& pQVec)
{
eVector += pQVec.eVector;
return *this;
}
auto& operator-=(QuantityVector<dim> const& pQVec)
{
eVector -= pQVec.eVector;
return *this;
}
auto& operator-() const
{
return QuantityVector<dim>(-eVector);
}
auto normalized() const
{
return (*this) * (1 / norm());
}
};
template <typename dim>
auto& operator<<(std::ostream& os, QuantityVector<dim> qv)
{
using Quantity = phys::units::quantity<dim, double>;
os << '(' << qv.eVector(0) << ' ' << qv.eVector(1) << ' ' << qv.eVector(2)
<< ") " << phys::units::to_unit_symbol<dim, double>(Quantity(phys::units::detail::magnitude_tag, 1));
return os;
}
#endif
#ifndef _include_SPHERE_H_
#define _include_SPHERE_H_
#include <Geometry/Point.h>
#include <Units/PhysicalUnits.h>
class Sphere
{
using Length = phys::units::quantity<phys::units::length_d, double>;
Point center;
Length const radius;
public:
Sphere(Point const& pCenter, Length const pRadius) :
center(pCenter), radius(pRadius)
{
}
auto isInside(Point const& p) const
{
return radius*radius > (center - p).squaredNorm();
}
};
#endif
#ifndef _include_VECTOR_H_
#define _include_VECTOR_H_
#include <Geometry/BaseVector.h>
#include <Geometry/QuantityVector.h>
#include <Units/PhysicalUnits.h>
template <typename dim>
class Vector : public BaseVector<dim>
{
using Quantity = phys::units::quantity<dim, double>;
public:
Vector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector) :
BaseVector<dim>(pCS, pQVector)
{
}
Vector(CoordinateSystem const& cs, Quantity x, Quantity y, Quantity z) :
BaseVector<dim>(cs, QuantityVector<dim>(x, y, z))
{
}
auto getComponents() const
{
return BaseVector<dim>::qVector;
}
auto getComponents(CoordinateSystem const& pCS) const
{
if (&pCS == BaseVector<dim>::cs)
{
return BaseVector<dim>::qVector;
}
else
{
return QuantityVector<dim>(CoordinateSystem::getTransformation(*BaseVector<dim>::cs, pCS).linear() * BaseVector<dim>::qVector.eVector);
}
}
void rebase(CoordinateSystem const& pCS)
{
BaseVector<dim>::qVector = getComponents(pCS);
BaseVector<dim>::cs = &pCS;
}
auto norm() const
{
return BaseVector<dim>::qVector.norm();
}
auto squaredNorm() const
{
return BaseVector<dim>::qVector.squaredNorm();
}
template <typename dim2>
auto parallelProjectionOnto(Vector<dim2> const& pVec, CoordinateSystem const& pCS) const
{
auto const ourCompVec = getComponents(pCS);
auto const otherCompVec = pVec.getComponents(pCS);
auto const& a = ourCompVec.eVector;
auto const& b = otherCompVec.eVector;
return Vector<dim>(pCS, QuantityVector<dim>(b * ((a.dot(b)) / b.squaredNorm())));
}
template <typename dim2>
auto parallelProjectionOnto(Vector<dim2> const& pVec) const
{
return parallelProjectionOnto<dim2>(pVec, *BaseVector<dim>::cs);
}
auto operator+(Vector<dim> const& pVec) const
{
auto const components = getComponents(*BaseVector<dim>::cs) + pVec.getComponents(*BaseVector<dim>::cs);
return Vector<dim>(*BaseVector<dim>::cs, components);
}
auto operator-(Vector<dim> const& pVec) const
{
auto const components = getComponents() - pVec.getComponents(*BaseVector<dim>::cs);
return Vector<dim>(*BaseVector<dim>::cs, components);
}
auto& operator*=(double const p)
{
BaseVector<dim>::qVector *= p;
return *this;
}
template <typename ScalarDim>
auto operator*(phys::units::quantity<ScalarDim, double> const p) const
{
using ProdQuantity = phys::units::detail::Product<dim, ScalarDim, double, double>;
if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless, not a "Quantity" anymore
{
return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
}
else
{
return Vector<typename ProdQuantity::dimension_type>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
}
}
template <typename ScalarDim>
auto operator/(phys::units::quantity<ScalarDim, double> const p) const
{
return (*this) * (1 / p);
}
auto operator*(double const p) const
{
return Vector<dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
}
auto operator/(double const p) const
{
return Vector<dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector / p);
}
auto& operator+=(Vector<dim> const& pVec)
{
BaseVector<dim>::qVector += pVec.getComponents(*BaseVector<dim>::cs);
return *this;
}
auto& operator-=(Vector<dim> const& pVec)
{
BaseVector<dim>::qVector -= pVec.getComponents(*BaseVector<dim>::cs);
return *this;
}
auto& operator-() const
{
return Vector<dim>(*BaseVector<dim>::cs, - BaseVector<dim>::qVector);
}
auto normalized() const
{
return (*this) * (1 / norm());
}
template <typename dim2>
auto cross(Vector<dim2> pV) const
{
auto const c1 = getComponents().eVector;
auto const c2 = pV.getComponents(*BaseVector<dim>::cs).eVector;
auto const bareResult = c1.cross(c2);
using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless, not a "Quantity" anymore
{
return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs, bareResult);
}
else
{
return Vector<typename ProdQuantity::dimension_type>(*BaseVector<dim>::cs, bareResult);
}
}
};
#endif
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one cpp file
#include <catch2/catch.hpp>
#include <Units/PhysicalUnits.h>
using namespace phys::units;
TEST_CASE( "PhysicalUnits", "[Units]" )
{
SECTION( "sectionOne" )
{
REQUIRE( 1_m/1_m == 1 );
}
SECTION( "sectionTwo" )
{
REQUIRE_FALSE( 1_m/1_m == 2 );
}
SECTION( "sectionThree" )
{
REQUIRE( 1_s/1_s == 2 );
}
}
#ifndef _include_BufferedSink_h_
#define _include_BufferedSink_h_
namespace logger {
namespace sink {
/**
Output buffer template. NoBuffer does nothingk.
*/
/*
struct NoBuffer {
inline bool Test(const std::string&) const { return false; }
inline std::string GetString() const { return std::string(""); }
inline void Clear() {}
inline void Add(const std::string&) {}
};
*/
/**
Output buffer template. StdBuffer records fSize characters in
local memeory before passing it on to further output stages.
*/
struct StdBuffer {
StdBuffer(const int size) : fSize(size) {}
inline bool Test(const std::string& s) { return int(fBuffer.tellp())+s.length() < fSize; }
inline std::string GetString() const { return fBuffer.str(); }
inline void Clear() { fBuffer.str(""); }
inline void Add(const std::string& s) { fBuffer << s; }
private:
int fSize;
std::ostringstream fBuffer;
};
/**
Definition of Sink for log output.
*/
template<typename TStream, typename TBuffer=StdBuffer>
class BufferedSink {
public:
BufferedSink(TStream& out, TBuffer buffer = {} ) : fOutput(out), fBuffer(std::move(buffer)) {}
void operator<<(const std::string& msg) {
if (!fBuffer.Test(msg)) {
fOutput << fBuffer.GetString();
fBuffer.Clear();
}
if (!fBuffer.Test(msg))
fOutput << msg;
else
fBuffer.Add(msg);
}
void Close() { fOutput << fBuffer.GetString(); }
private:
TStream& fOutput;
TBuffer fBuffer;
};
typedef BufferedSink<std::ostream, StdBuffer> BufferedSinkStream;
}// end namespace
} // end namespace
#endif
add_library (CORSIKAlogging INTERFACE)
target_include_directories (CORSIKAlogging INTERFACE ${PROJECT_SOURCE_DIR}/Framework)
target_include_directories (CORSIKAlogging INTERFACE $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Framework>
$<INSTALL_INTERFACE:include/Framework>
)
install (FILES Logger.h Sink.h MessageOn.h MessageOff.h NoSink.h Sink.h BufferedSink.h
DESTINATION include/Logging)
#ifndef _include_logger_h_
#define _include_logger_h_
#include <iosfwd>
#include <string>
#include <sstream>
#include <typeinfo>
#include <boost/format.hpp>
#include <Logging/MessageOn.h>
#include <Logging/MessageOff.h>
#include <Logging/Sink.h>
#include <Logging/NoSink.h>
#include <Logging/BufferedSink.h>
using namespace std;
using namespace boost;
/**
Everything around logfile generation and text output.
*/
namespace logger {
/**
Defines one stream to accept messages, and to wrote those into
TSink. The helper class MessageOn will convert input at
compile-time into message strings. The helper class MessageOff,
will just do nothing and will be optimized away at compile time.
*/
template<typename MSG=MessageOn, typename TSink=sink::NoSink>
class Logger : private MSG {
using MSG::Message;
public:
// Logger() : fName("") {}
Logger(const std::string color, const std::string name, TSink& sink) : fSink(sink), fName(color+"["+name+"]\033[39m ") {}
~Logger() { fSink.Close(); }
// Logger(const Logger&) = delete;
/**
Function to add string-concatenation of all inputs to output
sink.
*/
template<typename ... Strings>
void Log(const Strings&... inputs) {
fSink << MSG::Message(inputs...);
}
const std::string& GetName() const { return fName; }
private:
TSink& fSink;
std::string fName;
};
} // end namesapce
#define LOG(__LOGGER,...) \
__LOGGER.Log(__LOGGER.GetName(), __FILE__,":", __LINE__, " (", __func__, ") -> ", ##__VA_ARGS__);
#endif
CXXFLAGS+=-I. --std=c++14 -O3
all: test test_off
#test: test.o
# $(CXX) $(CXXFLAGS) $(LDFLAGS) $^ -o $@
clean:
rm -rf *.o test test_off *.dat *.log
#ifndef _include_MessageOff_h_
#define _include_MessageOff_h_
namespace logger {
/**
Helper class to ignore all arguments to MessagesOn::Message and
always return empty string "".
*/
class MessageOff {
protected:
template<typename First, typename ... Strings> std::string Message(const First& arg, const Strings&... rest) {
return "";
}
};
} // end namespace
#endif
#ifndef _include_MessageOn_h_
#define _include_MessageOn_h_
namespace logger {
/**
Helper class to convert all input arguments of MessageOn::Message
into string-concatenated version and return this as string.
*/
class MessageOn {
protected:
std::string Message() { return "\n"; }
template<typename First, typename ... Strings> std::string Message(const First& arg, const Strings&... rest) {
std::ostringstream ss;
ss << arg << Message(rest...);
return ss.str();
}
template<typename ... Strings> std::string Message(const int& arg, const Strings&... rest) {
return std::to_string(arg) + Message(rest...);
}
template<typename ... Strings> std::string Message(const double& arg, const Strings&... rest) {
return std::to_string(arg) + Message(rest...);
}
template<typename ... Strings> std::string Message(char const * arg, const Strings&... rest) {
return std::string(arg) + Message(rest...);
}
template<typename ... Strings> std::string Message(const std::string& arg, const Strings&... rest) {
return arg + Message(rest...);
}
// ----------------------
// boost format
template<typename ... Strings> std::string Message(const boost::format& fmt, const Strings&... rest) {
boost::format FMT(fmt);
return bformat(FMT, rest...);
}
template<typename Arg, typename ... Strings> std::string bformat(boost::format& fmt, const Arg& arg, const Strings&... rest) {
fmt % arg;
return bformat(fmt, rest...);
}
std::string bformat(boost::format& fmt) { return fmt.str() + "\n"; }
};
}// end namesapce
#endif
#ifndef _include_NoSink_h_
#define _include_NoSink_h_
namespace logger {
namespace sink {
struct NoSink {
inline void operator<<(const std::string&) {}
inline void Close() {}
};
}// end namespace
} // end namespace
#endif
#ifndef _include_Sink_h_
#define _include_Sink_h_
namespace logger {
/**
a sink for the logger must implement the two functions
operator<<(const std::string&)
and
Close()
See example: NoSink
*/
namespace sink {
/**
Definition of Sink for log output.
*/
template<typename TStream>
class Sink {
public:
Sink(TStream& out) : fOutput(out) {}
void operator<<(const std::string& msg) {
fOutput << msg;
}
void Close() {}
private:
TStream& fOutput;
};
typedef Sink<std::ostream> SinkStream;
}// end namespace
} // end namespace
#endif
set (STACK_HEADERS StackOne.h)
add_library (CORSIKAstack INTERFACE)
#set_target_properties (CORSIKAstack PROPERTIES VERSION ${PROJECT_VERSION})
#set_target_properties (CORSIKAstack PROPERTIES SOVERSION 1)
#set_target_properties (CORSIKAstack PROPERTIES PUBLIC_HEADER "${STACK_HEADERS}")
#target_link_libraries (CORSIKAstackinterface CORSIKAunits)
#target_include_directories (CORSIKAstack PRIVATE ${EIGEN3_INCLUDE_DIR})
target_include_directories (CORSIKAstack INTERFACE $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Framework>
$<INSTALL_INTERFACE:include/Framework>
)
install (FILES StackOne.h DESTINATION include/Stack)
#ifndef _include_StackIterator_h__
#define _include_StackIterator_h__
#include <iostream>
#include <iomanip>
namespace stack {
template<class Stack, class Particle> class StackIteratorInfo;
/**
The main interface to iterator over objects on a stack.
*/
template<typename Stack, typename Particle>
class StackIterator : public Particle
{
friend Stack;
friend Particle;
friend StackIteratorInfo<Stack,Particle>;
private:
int fIndex;
//#warning stacks should not be copied because of this:
Stack* fData;
public:
StackIterator() : fData(0), fIndex(0) { }
StackIterator(Stack& data, const int index) : fData(&data), fIndex(index) { }
StackIterator(const StackIterator& mit) : fData(mit.fData), fIndex(mit.fIndex) { }
StackIterator& operator++() { ++fIndex; return *this; }
StackIterator operator++(int) { StackIterator tmp(*this); ++fIndex; return tmp; }
bool operator==(const StackIterator& rhs) { return fIndex == rhs.fIndex; }
bool operator!=(const StackIterator& rhs) { return fIndex != rhs.fIndex; }
StackIterator& operator*() { return *this; }
const StackIterator& operator*() const { return *this; }
protected:
int GetIndex() const { return fIndex; }
Stack& GetStack() { return *fData; }
const Stack& GetStack() const { return *fData; }
inline StackIterator<Stack,Particle>& base_ref() { return static_cast<StackIterator<Stack, Particle>&>(*this); }
inline const StackIterator<Stack,Particle>& base_ref() const { return static_cast<const StackIterator<Stack, Particle>&>(*this); }
};
/**
Internal helper class for StackIterator
*/
template<class _Stack, class Particle>
class StackIteratorInfo {
friend Particle;
private:
StackIteratorInfo() {}
protected:
inline _Stack& Stack() { return static_cast<StackIterator<_Stack, Particle>*>(this)->GetStack(); }
inline int Index() const { return static_cast<const StackIterator<_Stack, Particle>*>(this)->GetIndex(); }
inline const _Stack& Stack() const { return static_cast<const StackIterator<_Stack, Particle>*>(this)->GetStack(); }
};
} // end namespace stack
#endif
#ifndef _include_stackone_h_
#define _include_stackone_h_
#include <vector>
#include <string>
#include <StackInterface/Stack.h>
namespace stack {
/**
Example of a particle object on the stack.
*/
template<typename _Stack>
class ParticleReadOne : public StackIteratorInfo<_Stack, ParticleReadOne<_Stack> >
{
using StackIteratorInfo<_Stack, ParticleReadOne>::GetIndex;
using StackIteratorInfo<_Stack, ParticleReadOne>::GetStack;
public:
void SetId(const int id) { GetStack().SetId(GetIndex(), id); }
void SetEnergy(const double e) { GetStack().SetEnergy(GetIndex(), e); }
int GetId() const { GetStack().GetId(GetIndex()); }
double GetEnergy() const { GetStack().GetEnergy(GetIndex()); }
double GetPDG() const { return 0; } // ConvertToPDG(GetId()); }
void SetPDG(double v) { GetStack().SetId(0, 0); } //fIndex, ConvertFromPDG(v)); }
};
/**
Memory implementation of the most simple particle stack object.
*/
class StackOneImpl
{
private:
/// the actual memory to store particle data
std::vector<int> fId;
std::vector<double> fData;
public:
void Clear() { fData.clear(); }
int GetSize() const { return fData.size(); }
int GetCapacity() const { return fData.size(); }
void SetId(const int i, const int id) { fId[i] = id; }
void SetEnergy(const int i, const double e) { fData[i] = e; }
const int GetId(const int i) const { return fId[i]; }
const double GetEnergy(const int i) const { return fData[i]; }
/**
Function to copy particle at location i2 in stack to i1
*/
void Copy(const int i1, const int i2) {
fData[i2] = fData[i1];
fId[i2] = fId[i1];
}
protected:
void IncrementSize() { fData.push_back(0.); fId.push_back(0.); }
void DecrementSize() { if (fData.size()>0) { fData.pop_back(); fId.pop_back(); } }
};
typedef StackIterator<StackOneImpl, ParticleReadOne<StackOneImpl> > ParticleOne;
typedef Stack<StackOneImpl, ParticleOne> StackOne;
} // end namespace
#endif