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 6 additions and 1830 deletions
#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one cpp file
#include <catch2/catch.hpp>
#include <fwk/Logger.h>
TEST_CASE( "Logging", "[Logging]" )
{
SECTION( "sectionOne" )
{
}
SECTION( "sectionTwo" )
{
}
SECTION( "sectionThree" )
{
}
}
set (PUBLIC_HEADER_FILES 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_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include>
)
#add_custom_command (
# CORSIKAstack
# PRE_BUILD
# COMMAND ${CMAKE_COMMAND} -E copy_if_different ${PUBLIC_HEADER_FILES} ${PROJECT_BINARY_DIR}/include/corsika})
# )
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>
//#include <corsika/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 { return GetStack().GetId(GetIndex()); }
double GetEnergy() const { return 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
add_custom_command (
OUTPUT ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc
COMMAND ${PROJECT_SOURCE_DIR}/Framework/Particles/pdxml_reader.py ${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleData.xml ${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleClassNames.xml
DEPENDS
pdxml_reader.py
ParticleData.xml
ParticleClassNames.xml
WORKING_DIRECTORY
${PROJECT_BINARY_DIR}/Framework/Particles/
COMMENT "Read PYTHIA8 particle data and produce C++ source code GeneratedParticleProperties.inc"
VERBATIM
)
# all public header files of library, includes automatic generated file(s)
set (
PARTICLE_HEADERS
Particles.h
${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc # this one is auto-generated
)
# we need 2nd list with just the header files from the source code directory
set (
PARTICLE_HEADERS_SOURCE
Particles.h
)
set (
PARTICLE_NAMESPACE
fwk
)
add_library (CORSIKAparticles INTERFACE)
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAparticles ${PARTICLE_NAMESPACE} ${PARTICLE_HEADERS})
target_include_directories (
CORSIKAparticles
INTERFACE
# ${EIGEN3_INCLUDE_DIR}
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include>
)
install (
FILES
${PARTICLE_HEADERS}
DESTINATION
include/Particles
)
# --------------------
# code unit testing
add_executable (
testParticles
testParticles.cc
${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc
)
target_link_libraries (
testParticles
CORSIKAparticles
CORSIKAunits
CORSIKAthirdparty # for catch2
)
add_test (
NAME
testParticles
COMMAND
testParticles
)
<chapter name="Particle Names">
<!-- For selected particles it is possible to specify the C++ class
names for a specific unique PDG code -->
<particle id="0" classname="unknown"> </particle> <!-- VOID in Pythia8 -->
<particle id="11" classname="Electron"> </particle>
<particle id="-11" classname="Positron"> </particle>
<particle id="22" classname="Gamma"> </particle>
</chapter>
/**
@file Particles.h
Interface to particle properties
*/
#ifndef _include_Particle_h_
#define _include_Particle_h_
#include <array>
#include <cstdint>
#include <iostream>
#include <fwk/GeneratedParticleProperties.inc>
namespace fwk {
/**
* @namespace particle
*
* The properties of all elementary particles is stored here. The data
* is taken from the Pythia ParticleData.xml file.
*
*/
namespace particle {
/**
* @function GetMass
*
* return mass of particle
*/
auto constexpr GetMass(InternalParticleCode const p)
{
return masses[static_cast<uint8_t const>(p)];
}
auto constexpr GetPDG(InternalParticleCode const p)
{
return pdg_codes[static_cast<uint8_t const>(p)];
}
auto constexpr GetElectricChargeNumber(InternalParticleCode const p)
{
return electric_charge[static_cast<uint8_t const>(p)] / 3;
}
auto constexpr GetElectricCharge(InternalParticleCode const p)
{
return GetElectricChargeNumber(p) * (phys::units::e);
}
auto const GetName(InternalParticleCode const p)
{
return names[static_cast<uint8_t const>(p)];
}
std::ostream& operator<< (std::ostream& stream, InternalParticleCode const p)
{
stream << GetName(p);
return stream;
}
} // end namespace
} // end namespace
#endif
#!/usr/bin/env python3
import sys, math, itertools, re, csv, pprint
import xml.etree.ElementTree as ET
from collections import OrderedDict
##############################################################
#
# reading xml input data, return line by line particle data
#
def parse(filename):
tree = ET.parse(filename)
root = tree.getroot()
for particle in root.iter("particle"):
name = particle.attrib["name"]
antiName = "Unknown"
# print (str(particle.attrib))
if ("antiName" in particle.attrib):
antiName = particle.attrib["antiName"]
# print ("found anti: " + name + " " + antiName + "\n" )
pdg_id = int(particle.attrib["id"])
mass = float(particle.attrib["m0"]) # GeV
electric_charge = int(particle.attrib["chargeType"]) # in units of e/3
decay_width = float(particle.attrib.get("mWidth", 0)) # GeV
lifetime = float(particle.attrib.get("tau0", math.inf)) # mm / c
yield (pdg_id, name, mass, electric_charge, antiName)
# TODO: read decay channels from child elements
if "antiName" in particle.attrib:
yield (-pdg_id, antiName, mass, -electric_charge, name)
##############################################################
#
# returns dict with particle codes and class names
#
def class_names(filename):
tree = ET.parse(filename)
root = tree.getroot()
map = {}
for particle in root.iter("particle"):
name = particle.attrib["classname"]
pdg_id = int(particle.attrib["id"])
map[pdg_id] = name
return map
##############################################################
#
# Automatically produce a string qualifying as C++ class name
#
def c_identifier(name):
orig = name
name = name.upper()
for c in "() ":
name = name.replace(c, "_")
name = name.replace("BAR", "_BAR")
name = name.replace("0", "_0")
name = name.replace("/", "_")
name = name.replace("*", "_STAR")
name = name.replace("'", "_PRIME")
name = name.replace("+", "_PLUS")
name = name.replace("-", "_MINUS")
while True:
tmp = name.replace("__", "_")
if tmp == name:
break
else:
name = tmp
pattern = re.compile(r'^[A-Z_][A-Z_0-9]*$')
if pattern.match(name):
return name.strip("_")
else:
raise Exception("could not generate C identifier for '{:s}'".format(orig))
# ########################################################
#
# returns dict containing all data from pythia-xml input
#
def build_pythia_db(filename, classnames):
particle_db = OrderedDict()
for (pdg, name, mass, electric_charge, antiName) in parse(filename):
c_id = "unknown"
if (pdg in classnames):
c_id = classnames[pdg]
else:
c_id = c_identifier(name)
#~ print(name, c_id, sep='\t', file=sys.stderr)
#~ enums += "{:s} = {:d}, ".format(c_id, corsika_id)
particle_db[c_id] = {
"name" : name,
"antiName" : antiName,
"pdg" : pdg,
"mass" : mass, # in GeV
"electric_charge" : electric_charge # in e/3
}
return particle_db
###############################################################
#
# return string with enum of all internal particle codes
#
def gen_internal_enum(pythia_db):
string = "enum class InternalParticleCode : uint8_t {\n"
for k in filter(lambda k: "ngc_code" in pythia_db[k], pythia_db):
string += " {key:s} = {code:d},\n".format(key = k, code = pythia_db[k]['ngc_code'])
string += "};"
return string
###############################################################
#
# return string with all data arrays
#
def gen_properties(pythia_db):
# number of particles, size of tables
string = "static constexpr std::size_t size = {size:d};\n".format(size = len(pythia_db))
string += "\n"
# particle masses table
string += "static constexpr std::array<quantity<energy_d> const, size> masses = {\n"
for p in pythia_db.values():
string += " {mass:f}_GeV, // {name:s}\n".format(mass = p['mass'], name = p['name'])
string += "};\n\n"
# PDG code table
string += "static constexpr std::array<PDGCode const, size> pdg_codes = {\n"
for p in pythia_db.values():
string += " {pdg:d}, // {name:s}\n".format(pdg = p['pdg'], name = p['name'])
string += "};\n"
# name string table
string += "static const std::array<std::string const, size> names = {\n"
for p in pythia_db.values():
string += " \"{name:s}\",\n".format(name = p['name'])
string += "};\n"
# electric charges table
string += "static constexpr std::array<int16_t, size> electric_charge = {\n"
for p in pythia_db.values():
string += " {charge:d},\n".format(charge = p['electric_charge'])
string += "};\n"
# anti-particle table
# string += "static constexpr std::array<size, size> anti_particle = {\n"
# for p in pythia_db.values():
# string += " {anti:d},\n".format(charge = p['anti_particle'])
# string += "};\n"
return string
###############################################################
#
# return string with a list of classes for all particles
#
def gen_classes(pythia_db):
string = "// list of C++ classes to access particle properties"
for cname in pythia_db:
antiP = 'unknown'
for cname_anti in pythia_db:
if (pythia_db[cname_anti]['name'] == pythia_db[cname]['antiName']):
antiP = cname_anti
break
string += "\n";
string += "/** @class " + cname + "\n\n"
string += " * Particle properties are taken from the PYTHIA8 ParticleData.xml file:<br>\n"
string += " * - pdg=" + str(pythia_db[cname]['pdg']) +"\n"
string += " * - mass=" + str(pythia_db[cname]['mass']) + " GeV \n"
string += " * - charge= " + str(pythia_db[cname]['electric_charge']/3) + " \n"
string += " * - name=" + str(cname) + "\n"
string += " * - anti=" + str(antiP) + "\n"
string += "*/\n\n"
string += "class " + cname + "{\n"
string += " public:\n"
string += " static InternalParticleCode GetType() { return Type; }\n"
string += " static quantity<energy_d> GetMass() { return masses[TypeIndex]; }\n"
string += " static quantity<electric_charge_d> GetCharge() { return phys::units::e*electric_charge[TypeIndex]/3; }\n"
string += " static int GetChargeNumber() { return electric_charge[TypeIndex]/3; }\n"
string += " static std::string GetName() { return names[TypeIndex]; }\n"
string += " static InternalParticleCode GetAntiParticle() { return AntiType; }\n"
string += " static const InternalParticleCode Type = InternalParticleCode::" + cname + ";\n"
string += " static const InternalParticleCode AntiType = InternalParticleCode::" + antiP + ";\n"
string += " private:\n"
string += " static const uint8_t TypeIndex = static_cast<uint8_t const>(Type);\n"
string += "};\n"
return string
###############################################################
#
#
def inc_start():
string = ""
string += "#ifndef _include_GeneratedParticleDataTable_h_\n"
string += "#define _include_GeneratedParticleDataTable_h_\n\n"
string += "#include <array>\n"
string += "#include <cstdint>\n"
string += "#include <iostream>\n\n"
string += "using namespace phys::units;\n"
string += "using namespace phys::units::literals;\n\n"
string += "namespace fwk { \n\n"
string += "namespace particle { \n\n"
string += "typedef int16_t PDGCode;\n\n"
return string
###############################################################
#
#
def inc_end():
string = ""
string += "\n}\n\n"
string += "\n}\n\n"
string += "#endif\n"
return string
###################################################################
#
# Main function
#
if __name__ == "__main__":
if (len(sys.argv)!=3) :
print ("pdxml_reader.py Pythia8.xml ClassNames.xml")
sys.exit(0)
names = class_names(sys.argv[2])
pythia_db = build_pythia_db(sys.argv[1], names)
print ("\n pdxml_reader.py: Automatically produce particle-properties from PYTHIA8 xml file\n")
counter = itertools.count(0)
not_modeled = []
for p in pythia_db:
pythia_db[p]['ngc_code'] = next(counter)
with open("GeneratedParticleProperties.inc", "w") as f:
print(inc_start(), file=f)
print(gen_internal_enum(pythia_db), file=f)
print(gen_properties(pythia_db), file=f)
print(gen_classes(pythia_db), file=f)
print(inc_end(), file=f)
#~ print(pdg_id_table, mass_table, name_table, enums, sep='\n\n')
#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>
#include <fwk/Particles.h>
using namespace phys::units;
using namespace phys::units::literals;
using namespace fwk::particle;
TEST_CASE( "Particles", "[Particles]" )
{
SECTION( "Types" )
{
REQUIRE( Electron::GetType()==InternalParticleCode::Electron );
}
SECTION( "Data" )
{
REQUIRE( Electron::GetMass()/0.511_MeV==Approx(1) );
REQUIRE( Electron::GetMass()/GetMass(InternalParticleCode::Electron)==Approx(1) );
REQUIRE( Electron::GetCharge()/phys::units::e==Approx(-1) );
REQUIRE( Positron::GetCharge()/phys::units::e==Approx(+1) );
REQUIRE( GetElectricCharge(Positron::GetAntiParticle())/phys::units::e==Approx(-1) );
REQUIRE( Electron::GetName() == "e-" );
}
}
add_library (CORSIKAprocesssequence INTERFACE)
target_include_directories (CORSIKAprocesssequence INTERFACE $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Framework>
$<INSTALL_INTERFACE:include/Framework>
)
install (FILES ProcessSequence.h
DESTINATION include/ProcessSequence)
#ifndef _include_ProcessSequence_h_
#define _include_ProcessSequence_h_
#include <iostream>
#include <typeinfo>
using namespace std;
namespace processes {
/**
/class Base
The structural base type of a process object in a
ProcessSequence. Both, the ProcessSequence and all its elements
are of type Base<T>
*/
template <typename derived>
struct Base
{
const derived& GetRef() const
{
return static_cast<const derived&>(*this);
}
};
/**
\class ProcessSequence
A compile time static list of processes. The compiler will
generate a new type based on template logic containing all the
elements.
\comment Using CRTP pattern, https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
*/
template <typename T1, typename T2>
class ProcessSequence : public Base <ProcessSequence<T1,T2> >
{
public:
const T1& A;
const T2& B;
ProcessSequence(const T1& in_A, const T2& in_B)
: A(in_A)
, B(in_B)
{ }
template<typename D>
inline void DoContinuous(D& d) const { A.DoContinuous(d); B.DoContinuous(d); } // add trajectory
template<typename D>
inline double MinStepLength(D& d) const { return min(A.MinStepLength(d), B.MinStepLength(d)); }
//template<typename D>
//inline Trajectory Transport(D& d, double& length) const { A.Transport(d, length); B.Transport(d, length); }
template<typename D>
inline void DoDiscrete(D& d) const { A.DoDiscrete(d); B.DoDiscrete(d); }
};
template <typename T1, typename T2>
inline
const ProcessSequence<T1,T2>
operator+ (const Base<T1>& A, const Base<T2>& B)
{
return ProcessSequence<T1,T2>( A.GetRef(), B.GetRef() );
}
/*
template <typename T1>
struct depth_lhs
{
static const int num = 0;
};
// terminating condition
template <typename T1, typename T2>
struct depth_lhs< Sequence<T1,T2> >
{
// try to expand the left node (T1) which might be a Sequence type
static const int num = 1 + depth_lhs<T1>::num;
};
*/
/*
template <typename T1>
struct mat_ptrs
{
static const int num = 0;
inline static void
get_ptrs(const Process** ptrs, const T1& X)
{
ptrs[0] = reinterpret_cast<const Process*>(&X);
}
};
template <typename T1, typename T2>
struct mat_ptrs< Sequence<T1,T2> >
{
static const int num = 1 + mat_ptrs<T1>::num;
inline static void
get_ptrs(const Process** in_ptrs, const Sequence<T1,T2>& X)
{
// traverse the left node
mat_ptrs<T1>::get_ptrs(in_ptrs, X.A);
// get address of the matrix on the right node
in_ptrs[num] = reinterpret_cast<const Process*>(&X.B);
}
};
*/
/*
template<typename T1, typename T2>
const Process&
Process::operator=(const Sequence<T1,T2>& X)
{
int N = 1 + depth_lhs< Sequence<T1,T2> >::num;
const Process* ptrs[N];
mat_ptrs< Sequence<T1,T2> >::get_ptrs(ptrs, X);
int r = ptrs[0]->rows;
int c = ptrs[0]->cols;
// ... check that all matrices have the same size ...
set_size(r, c);
for(int j=0; j<r*c; ++j)
{
double sum = ptrs[0]->data[j];
for(int i=1; i<N; ++i)
{
sum += ptrs[i]->data[j];
}
data[j] = sum;
}
return *this;
}
*/
} // end namespace
#endif
add_library (CORSIKAstack INTERFACE)
target_include_directories (CORSIKAstack INTERFACE $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Framework>
$<INSTALL_INTERFACE:include/Framework>
)
install (FILES Stack.h StackIterator.h
DESTINATION include/Stack)
set (
PUBLIC_HEADER_FILES
Stack.h
StackIterator.h
)
add_library (
CORSIKAstackinterface
INTERFACE
)
target_include_directories (
CORSIKAstackinterface
INTERFACE
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/include>
$<INSTALL_INTERFACE:include>
)
#add_custom_command (
# TARGET
# CORSIKAstackinterface
# PRE_BUILD
# COMMAND ${CMAKE_COMMAND} -E copy_if_different ${PUBLIC_HEADER_FILES} ${PROJECT_BINARY_DIR}/include/corsika
# )
install (
FILES
${PUBLIC_HEADER_FILES}
DESTINATION
include/corsika
)
#ifndef _include_Stack_h__
#define _include_Stack_h__
#include <StackInterface/StackIterator.h> // to help application programmres
//#include <corsika/StackIterator.h> // to help application programmres
/**
All classes around management of particles on a stack.
*/
namespace stack {
/**
Interface definition of a Stack object. The Stack implements the
std-type begin/end function to allow integration in normal for
loops etc.
*/
template<typename DataImpl, typename Particle>
class Stack : public DataImpl {
public:
using DataImpl::GetCapacity;
using DataImpl::GetSize;
using DataImpl::Clear;
using DataImpl::Copy;
using DataImpl::IncrementSize;
using DataImpl::DecrementSize;
public:
typedef Particle iterator;
typedef const Particle const_iterator;
/// these are functions required by std containers and std loops
iterator begin() { return iterator(*this, 0); }
iterator end() { return iterator(*this, GetSize()); }
iterator last() { return iterator(*this, GetSize()-1); }
/// these are functions required by std containers and std loops
const_iterator cbegin() const { return const_iterator(*this, 0); }
const_iterator cend() const { return const_iterator(*this, GetSize()); }
const_iterator clast() const { return const_iterator(*this, GetSize()-1); }
/// increase stack size, create new particle at end of stack
iterator NewParticle() { IncrementSize(); return iterator(*this, GetSize()-1); }
/// delete last particle on stack by decrementing stack size
void DeleteLast() { DecrementSize(); }
};
} // end namespace
#endif
#ifndef _include_StackIterator_h__
#define _include_StackIterator_h__
#include <iostream>
#include <iomanip>
namespace stack {
// forward decl.
template<class Stack, class Particle> class StackIteratorInfo;
/**
@class StackIterator
The StackIterator is the main interface to iterator over
particles on a stack. At the same time StackIterator is a
Particle object by itself, thus there is no difference between
type and ref_type for convenience of the physicist.
This allows to write code like
\verbatim
for (auto& p : theStack) { p.SetEnergy(newEnergy); }
\endverbatim
The template argument Stack determines the type of Stack object
the data is stored in. A pointer to the Stack object is part of
the StackIterator. In addition to Stack the iterator only knows
the index fIndex in the Stack data.
The template argument Particles acts as a policy to provide
readout function of Particle data from the stack. The Particle
class must know how to retrieve information from the Stack data
for a particle entry at any index fIndex.
*/
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; }
// this is probably not needed rigth now:
//inline StackIterator<Stack,Particle>& BaseRef() { return static_cast<StackIterator<Stack, Particle>&>(*this); }
//inline const StackIterator<Stack,Particle>& BaseRef() const { return static_cast<const StackIterator<Stack, Particle>&>(*this); }
};
/**
@class StackIteratorInfo
This is the class where custom ...
Internal helper class for StackIterator. Document better...
*/
template<typename _Stack, typename Particle>
class StackIteratorInfo {
friend Particle;
private:
StackIteratorInfo() {}
protected:
inline _Stack& GetStack() { return static_cast<StackIterator<_Stack, Particle>*>(this)->GetStack(); }
inline int GetIndex() const { return static_cast<const StackIterator<_Stack, Particle>*>(this)->GetIndex(); }
inline const _Stack& GetStack() const { return static_cast<const StackIterator<_Stack, Particle>*>(this)->GetStack(); }
};
} // end namespace stack
#endif
add_library (CORSIKAunits INTERFACE)
target_include_directories (CORSIKAunits
INTERFACE
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Framework>
$<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/ThirdParty>
$<INSTALL_INTERFACE:include>
)
install (FILES PhysicalUnits.h DESTINATION include/Units)
# code testing
add_executable (testUnits testUnits.cc)
target_link_libraries (testUnits CORSIKAunits CORSIKAthirdparty) # for catch2
add_test(NAME testUnits COMMAND testUnits)
#ifndef _include_PhysicalUnits_h_
#define _include_PhysicalUnits_h_
#include <phys/units/quantity.hpp>
#include <phys/units/io.hpp>
#include <phys/units/physical_constants.hpp>
/**
@file PhysicalUnits
Define _XeV literals, alowing 10_GeV in the code.
*/
/*using namespace phys::units::io;
using namespace phys::units::literals;*/
namespace phys {
namespace units {
namespace literals {
QUANTITY_DEFINE_SCALING_LITERALS(eV, energy_d, magnitude(eV) )
}
}
}
using Length = phys::units::quantity<phys::units::length_d, double>;
using Time = phys::units::quantity<phys::units::time_interval_d, double>;
using Speed = phys::units::quantity<phys::units::speed_d, double>;
using Frequency = phys::units::quantity<phys::units::frequency_d, double>;
using ElectricCharge = phys::units::quantity<phys::units::electric_charge_d, double>;
#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>
#include <array>
using namespace phys::units;
using namespace phys::units::literals;
TEST_CASE( "PhysicalUnits", "[Units]" ) {
SECTION("Consistency") {
REQUIRE( 1_m/1_m == Approx(1) );
//REQUIRE_FALSE( 1_m/1_s == 1 ); // static assert
}
SECTION("Constructors") {
auto E1 = 10_GeV;
REQUIRE( E1==10_GeV );
quantity<phys::units::length_d> l1 = 10_nm;
quantity<phys::units::length_d> arr0[5];
arr0[0] = 5_m;
quantity<phys::units::length_d> arr1[2] = { {1_mm}, {2_cm} };
std::array<quantity<phys::units::energy_d>, 4> arr2; // empty array
std::array<quantity<phys::units::energy_d>, 4> arr3 = { 1_GeV, 1_eV, 5_MeV };
}
SECTION("Powers in literal units") {
REQUIRE( 1_s/1_ds == Approx(1e1) );
REQUIRE( 1_m/1_cm == Approx(1e2) );
REQUIRE( 1_m/1_mm == Approx(1e3) );
REQUIRE( 1_V/1_uV == Approx(1e6) );
REQUIRE( 1_s/1_ns == Approx(1e9) );
REQUIRE( 1_eV/1_peV == Approx(1e12) );
REQUIRE( 1_A/1_fA == Approx(1e15) );
REQUIRE( 1_mol/1_amol == Approx(1e18) );
REQUIRE( 1_K/1_zK == Approx(1e21) );
REQUIRE( 1_K/1_yK == Approx(1e24) );
REQUIRE( 1_A/1_hA == Approx(1e-2) );
REQUIRE( 1_m/1_km == Approx(1e-3) );
REQUIRE( 1_m/1_Mm == Approx(1e-6) );
REQUIRE( 1_V/1_GV == Approx(1e-9) );
REQUIRE( 1_s/1_Ts == Approx(1e-12) );
REQUIRE( 1_eV/1_PeV == Approx(1e-15) );
REQUIRE( 1_A/1_EA == Approx(1e-18) );
REQUIRE( 1_K/1_ZK == Approx(1e-21) );
REQUIRE( 1_mol/1_Ymol == Approx(1e-24) );
}
SECTION("Powers and units") {
REQUIRE( 1 * ampere / 1_A == Approx(1e0) );
REQUIRE( mega * bar / bar == Approx(1e6) );
}
SECTION("Formulas") {
const quantity<phys::units::energy_d> E2 = 20_GeV * 2;
REQUIRE( E2==40_GeV );
const double lgE = log10(E2/1_GeV);
REQUIRE( lgE==Approx(log10(40.)) );
const auto E3 = E2 + 100_GeV + pow(10, lgE) * 1_GeV;
REQUIRE( E3==180_GeV );
}
}
This diff is collapsed.