diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000000000000000000000000000000000000..762ac959fa307f188f4cad00297a4ec9b7aeadd3 --- /dev/null +++ b/.clang-format @@ -0,0 +1,89 @@ +--- +Language: Cpp +AccessModifierOffset: -2 +AlignAfterOpenBracket: Align +AlignConsecutiveAssignments: false +AlignConsecutiveDeclarations: false +AlignEscapedNewlinesLeft: true +AlignOperands: true +AlignTrailingComments: true +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortBlocksOnASingleLine: true +AllowShortCaseLabelsOnASingleLine: false +AllowShortFunctionsOnASingleLine: All +AllowShortIfStatementsOnASingleLine: true +AllowShortLoopsOnASingleLine: true +AlwaysBreakAfterDefinitionReturnType: None +AlwaysBreakAfterReturnType: None +AlwaysBreakBeforeMultilineStrings: true +AlwaysBreakTemplateDeclarations: true +BinPackArguments: true +BinPackParameters: true +BraceWrapping: + AfterClass: false + AfterControlStatement: false + AfterEnum: false + AfterFunction: false + AfterNamespace: true + AfterObjCDeclaration: false + AfterStruct: false + AfterUnion: false + BeforeCatch: false + BeforeElse: false + IndentBraces: false +BreakBeforeBinaryOperators: None +BreakBeforeBraces: Attach +BreakBeforeTernaryOperators: true +BreakConstructorInitializersBeforeComma: true +# BreakInheritanceListBeforeComma: true +ColumnLimit: 90 +CommentPragmas: '^ IWYU pragma:' +ConstructorInitializerAllOnOneLineOrOnePerLine: false +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: true +DerivePointerAlignment: false +DisableFormat: false +ExperimentalAutoDetectBinPacking: false +ForEachMacros: [ foreach, Q_FOREACH, BOOST_FOREACH ] +IncludeCategories: + - Regex: '^<.*\.h>' + Priority: 1 + - Regex: '^<.*' + Priority: 2 + - Regex: '.*' + Priority: 3 +IndentCaseLabels: true +IndentWidth: 2 +IndentWrappedFunctionNames: false +KeepEmptyLinesAtTheStartOfBlocks: true +MacroBlockBegin: '' +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: All +ObjCBlockIndentWidth: 2 +ObjCSpaceAfterProperty: false +ObjCSpaceBeforeProtocolList: false +PenaltyBreakBeforeFirstCallParameter: 1 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakString: 1000 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 200 +PointerAlignment: Left +ReflowComments: true +SortIncludes: true +SpaceAfterCStyleCast: false +SpaceBeforeAssignmentOperators: true +SpaceBeforeParens: ControlStatements +SpaceInEmptyParentheses: false +SpacesBeforeTrailingComments: 1 +SpacesInAngles: false +SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false +SpacesInParentheses: false +SpacesInSquareBrackets: false +Standard: Auto +TabWidth: 8 +UseTab: Never +... diff --git a/.gitignore b/.gitignore index ccf6c21528c26c60f368479b82c9f2f82aaa466f..9656bc848ff52f63938067173e1372348527c92a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ **/*~ build/ +Framework/Particles/GeneratedParticleProperties.inc diff --git a/CMakeLists.txt b/CMakeLists.txt index 7dbe1c5a4f4f1def70927102343d158f662727e3..5a38c4ecb422e5ea823f864deaf0df7acceb5ac9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,29 +1,35 @@ cmake_minimum_required (VERSION 3.4.3) -project (corsika VERSION 8.0.0 DESCRIPTION "CORSIKA C++ project" LANGUAGES CXX) +project ( + corsika + VERSION 8.0.0 + DESCRIPTION "CORSIKA C++ project" + LANGUAGES CXX + ) # ignore many irrelevant Up-to-date messages during install set (CMAKE_INSTALL_MESSAGE LAZY) # directory for local cmake modules set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/CMakeModules) +include (CorsikaUtilities) # a few cmake function # --std=c++17 set (CMAKE_CXX_STANDARD 17) enable_testing () set (CTEST_OUTPUT_ON_FAILURE 1) -# testing coverage -include(CodeCoverage) -#set(COVERAGE_LCOV_EXCLUDES 'Documentation/*') -#setup_target_for_coverage(${PROJECT_NAME}_coverage ${PROJECT_TEST_NAME} coverage) -SETUP_TARGET_FOR_COVERAGE_GCOVR_HTML( - NAME corsika_coverage - EXECUTABLE ctest - #-j ${PROCESSOR_COUNT} -# DEPENDENCIES corsika - ) - +# unit testing coverage, does not work yet +#include (CodeCoverage) +##set(COVERAGE_LCOV_EXCLUDES 'Documentation/*') +##setup_target_for_coverage(${PROJECT_NAME}_coverage ${PROJECT_TEST_NAME} coverage) +#SETUP_TARGET_FOR_COVERAGE_GCOVR_HTML ( +# NAME corsika_coverage +# EXECUTABLE ctest +# #-j ${PROCESSOR_COUNT} +# # DEPENDENCIES corsika +# ) + #add_custom_target (corsika_pre_build) #add_custom_command (TARGET corsika_pre_build PRE_BUILD COMMAND "${PROJECT_SOURCE_DIR}/pre_compile.py") @@ -34,8 +40,9 @@ find_package (Eigen3 REQUIRED) # order of subdirectories add_subdirectory (ThirdParty) -add_subdirectory (Documentation) add_subdirectory (Framework) -#add_subdirectory (Processes) +add_subdirectory (Stack) +add_subdirectory (Processes) +add_subdirectory (Documentation) add_subdirectory (Main) diff --git a/CMakeModules/CorsikaUtilities.cmake b/CMakeModules/CorsikaUtilities.cmake new file mode 100644 index 0000000000000000000000000000000000000000..f98a7585468a68c1845192c4ff779debb5396a76 --- /dev/null +++ b/CMakeModules/CorsikaUtilities.cmake @@ -0,0 +1,57 @@ + +# +# takes a list of input files and prepends a path +# + +function (CORSIKA_PREPEND_PATH return prefix) + set (listVar "") + foreach (f ${ARGN}) + list (APPEND listVar "${prefix}/${f}") + endforeach (f) + set (${return} "${listVar}" PARENT_SCOPE) +endfunction (CORSIKA_PREPEND_PATH) + + +# +# use: CORSIKA_COPY_HEADERS_TO_NAMESPACE theLib theNamesapce header1.h header2.h ... +# +# creates a dependence of theLib on the presence of all listed header files in the +# build-directory include/theNamespace directory +# +# if needed, create symbolic links from the source files to this build-directory location +# +# any path information from input filenames is stripped, IF path was specified it is used for the link destination, if NOT the link is relative to the CMAKE_CURRENT_SOURCE_DIR +# +function (CORSIKA_COPY_HEADERS_TO_NAMESPACE for_library in_namespace) + set (HEADERS_BUILD "") + foreach (HEADER ${ARGN}) + # find eventual path, and handle it specificly + get_filename_component (barename ${HEADER} NAME) + get_filename_component (baredir ${HEADER} DIRECTORY) + # remove path, prepend build-directory destination + list (APPEND HEADERS_BUILD "${PROJECT_BINARY_DIR}/include/${in_namespace}/${barename}") + # define source location, use path if specified, otherwise CMAKE_CURRENT_SOURCE_DIR + set (FROM_DIR ${CMAKE_CURRENT_SOURCE_DIR}) + if (NOT "${baredir}" STREQUAL "") + set (FROM_DIR ${baredir}) + endif () + # define command to perform the symbolic linking + add_custom_command ( + OUTPUT ${PROJECT_BINARY_DIR}/include/${in_namespace}/${barename} + COMMAND ${CMAKE_COMMAND} -E make_directory ${PROJECT_BINARY_DIR}/include/${in_namespace} + COMMAND ${CMAKE_COMMAND} -E create_symlink ${FROM_DIR}/${barename} ${PROJECT_BINARY_DIR}/include/${in_namespace}/${barename} + COMMENT "Generate link: ${PROJECT_BINARY_DIR}/include/${in_namespace}/${barename}" + ) + endforeach (HEADER) + + # main target for this build step, depends on header files in build area + add_custom_target ( + ${for_library}_BUILD + DEPENDS ${HEADERS_BUILD} + ) + + # connect the _BUILD target to the main (external) target + add_dependencies (${for_library} ${for_library}_BUILD) + +endfunction (CORSIKA_COPY_HEADERS_TO_NAMESPACE) + diff --git a/Documentation/Doxygen/Doxyfile.in b/Documentation/Doxygen/Doxyfile.in index 42bcd4f961972d8dcba4dd253261608e2da90bfc..47bcf61eb06572fab45ea308f3bf9091bbae512f 100644 --- a/Documentation/Doxygen/Doxyfile.in +++ b/Documentation/Doxygen/Doxyfile.in @@ -2,22 +2,24 @@ PROJECT_NAME = CORSIKA PROJECT_NUMBER = 8.0.0 OUTPUT_DIRECTORY = @CMAKE_CURRENT_BINARY_DIR@/ -INPUT = @CMAKE_CURRENT_SOURCE_DIR@/../.. -EXCLUDE_PATTERN = "*/ThirdParty/*/*" +INPUT = @PROJECT_SOURCE_DIR@ @PROJECT_BINARY_DIR@/Framework +EXCLUDE_PATTERNS = */ThirdParty/*/* GENERATE_HTML = YES GENERATE_LATEX = YES -FILE_PATTERNS = *.cc *.cpp *.cxx *.h *.dox +FILE_PATTERNS = *.cc *.cpp *.cxx *.h *.dox *.inc +EXTENSION_MAPPING = inc=C++ RECURSIVE = YES SOURCE_BROWSER = YES # INLINE_SOURCES -CLASS_DIAGRAMS = YES +CLASS_DIAGRAMS = NO HAVE_DOT = YES CLASS_GRAPH = YES -UML_LOOK = YES +COLLABORATION_GRAPH = NO +UML_LOOK = NO TEMPLATE_RELATIONS = YES INCLUDE_GRAPH = YES GRAPHICAL_HIERARCHY = YES diff --git a/Documentation/Examples/CMakeLists.txt b/Documentation/Examples/CMakeLists.txt index cfd96473717afafb31fd4cc52fe266fe48fef024..d4efb2fc768ac5226f9160fcf2454ff37f8db92f 100644 --- a/Documentation/Examples/CMakeLists.txt +++ b/Documentation/Examples/CMakeLists.txt @@ -11,10 +11,10 @@ target_link_libraries (logger_example CORSIKAunits CORSIKAlogging) install (TARGETS logger_example DESTINATION share/examples) add_executable (stack_example stack_example.cc) -target_link_libraries (stack_example CORSIKAstack CORSIKAunits CORSIKAlogging) +target_link_libraries (stack_example SuperStupidStack CORSIKAunits CORSIKAlogging) install (TARGETS stack_example DESTINATION share/examples) add_executable (staticsequence_example staticsequence_example.cc) -target_link_libraries (staticsequence_example CORSIKAstack CORSIKAprocesssequence CORSIKAunits CORSIKAlogging) +target_link_libraries (staticsequence_example SuperStupidStack CORSIKAprocesssequence CORSIKAunits CORSIKAlogging) install (TARGETS staticsequence_example DESTINATION share/examples) diff --git a/Documentation/Examples/geometry_example.cc b/Documentation/Examples/geometry_example.cc index b2c2404922b4f66f4d53749bac90b1094815d1fb..9c459db161a116047b87e1ffd18e00c753fef374 100644 --- a/Documentation/Examples/geometry_example.cc +++ b/Documentation/Examples/geometry_example.cc @@ -1,63 +1,66 @@ -#include <Geometry/Vector.h> -#include <Geometry/Sphere.h> -#include <Geometry/Point.h> -#include <Geometry/CoordinateSystem.h> -#include <Units/PhysicalUnits.h> +#include <corsika/geometry/CoordinateSystem.h> +#include <corsika/geometry/Point.h> +#include <corsika/geometry/Sphere.h> +#include <corsika/geometry/Vector.h> +#include <corsika/units/PhysicalUnits.h> +#include <cstdlib> #include <iostream> #include <typeinfo> -#include <cstdlib> -using namespace phys::units; -using namespace phys::units::io; // support stream << unit -using namespace phys::units::literals; // support unit literals like 5_m; - -int main() -{ - // define the root coordinate system - CoordinateSystem root; - - // another CS defined by a translation relative to the root CS - CoordinateSystem cs2 = root.translate({0_m, 0_m, 1_m}); - - // rotations are possible, too; parameters are axis vector and angle - CoordinateSystem cs3 = root.rotate({1_m, 0_m, 0_m}, 90 * degree_angle); - - // now let's define some geometrical objects: - Point const p1(root, {0_m, 0_m, 0_m}); // the origin of the root CS - Point const p2(cs2, {0_m, 0_m, 0_m}); // the origin of cs2 - - Vector<length_d> const diff = p2 - p1; // the distance between the points, basically the translation vector given above - auto const norm = diff.squaredNorm(); // squared length with the right dimension - - // print the components of the vector as given in the different CS - std::cout << "p2-p1 components in root: " << diff.GetComponents(root) << std::endl; - std::cout << "p2-p1 components in cs2: " << diff.GetComponents(cs2) << std::endl; // by definition invariant under translations - std::cout << "p2-p1 components in cs3: " << diff.GetComponents(cs3) << std::endl; // but not under rotations - std::cout << "p2-p1 norm^2: " << norm << std::endl; - - Sphere s(p1, 10_m); // define a sphere around a point with a radius - std::cout << "p1 inside s: " << s.isInside(p2) << std::endl; - - Sphere s2(p1, 3_um); // another sphere - std::cout << "p1 inside s2: " << s2.isInside(p2) << std::endl; - - - // let's try parallel projections: - auto const v1 = Vector<length_d>(root, {1_m, 1_m, 0_m}); - auto const v2 = Vector<length_d>(root, {1_m, 0_m, 0_m}); - - auto const v3 = v1.parallelProjectionOnto(v2); - - // cross product - auto const cross = v1.cross(v2).normalized(); // normalized() returns dimensionless, normalized vectors - - // if a CS is not given as parameter for getComponents(), the components - // in the "home" CS are returned - std::cout << "v1: " << v1.GetComponents() << std::endl; - std::cout << "v2: " <<v2.GetComponents() << std::endl; - std::cout << "parallel projection of v1 onto v2: " << v3.GetComponents() << std::endl; - std::cout << "normalized cross product of v1 x v2" << cross.GetComponents() << std::endl; - - return EXIT_SUCCESS; +using namespace corsika::geometry; +using namespace corsika::units; + +int main() { + // define the root coordinate system + CoordinateSystem root; + + // another CS defined by a translation relative to the root CS + CoordinateSystem cs2 = root.translate({0_m, 0_m, 1_m}); + + // rotations are possible, too; parameters are axis vector and angle + CoordinateSystem cs3 = root.rotate({1_m, 0_m, 0_m}, 90 * degree_angle); + + // now let's define some geometrical objects: + Point const p1(root, {0_m, 0_m, 0_m}); // the origin of the root CS + Point const p2(cs2, {0_m, 0_m, 0_m}); // the origin of cs2 + + Vector<length_d> const diff = + p2 - + p1; // the distance between the points, basically the translation vector given above + auto const norm = diff.squaredNorm(); // squared length with the right dimension + + // print the components of the vector as given in the different CS + std::cout << "p2-p1 components in root: " << diff.GetComponents(root) << std::endl; + std::cout << "p2-p1 components in cs2: " << diff.GetComponents(cs2) + << std::endl; // by definition invariant under translations + std::cout << "p2-p1 components in cs3: " << diff.GetComponents(cs3) + << std::endl; // but not under rotations + std::cout << "p2-p1 norm^2: " << norm << std::endl; + + Sphere s(p1, 10_m); // define a sphere around a point with a radius + std::cout << "p1 inside s: " << s.isInside(p2) << std::endl; + + Sphere s2(p1, 3_um); // another sphere + std::cout << "p1 inside s2: " << s2.isInside(p2) << std::endl; + + // let's try parallel projections: + auto const v1 = Vector<length_d>(root, {1_m, 1_m, 0_m}); + auto const v2 = Vector<length_d>(root, {1_m, 0_m, 0_m}); + + auto const v3 = v1.parallelProjectionOnto(v2); + + // cross product + auto const cross = + v1.cross(v2).normalized(); // normalized() returns dimensionless, normalized vectors + + // if a CS is not given as parameter for getComponents(), the components + // in the "home" CS are returned + std::cout << "v1: " << v1.GetComponents() << std::endl; + std::cout << "v2: " << v2.GetComponents() << std::endl; + std::cout << "parallel projection of v1 onto v2: " << v3.GetComponents() << std::endl; + std::cout << "normalized cross product of v1 x v2" << cross.GetComponents() + << std::endl; + + return EXIT_SUCCESS; } diff --git a/Documentation/Examples/helix_example.cc b/Documentation/Examples/helix_example.cc index ae7250b16ffefbcf22dcb20f360ceb14debdad54..5caf35a43cf7d7a225161c02c4ecacf6a35475f5 100644 --- a/Documentation/Examples/helix_example.cc +++ b/Documentation/Examples/helix_example.cc @@ -1,46 +1,45 @@ -#include <Units/PhysicalUnits.h> -#include <Geometry/Vector.h> -#include <Geometry/CoordinateSystem.h> -#include <Geometry/Point.h> -#include <Geometry/Helix.h> +#include <corsika/geometry/CoordinateSystem.h> +#include <corsika/geometry/Helix.h> +#include <corsika/geometry/Point.h> +#include <corsika/geometry/Vector.h> +#include <corsika/units/PhysicalUnits.h> +#include <array> #include <cstdlib> #include <iostream> -#include <array> -using namespace phys::units; -using namespace phys::units::io; // support stream << unit -using namespace phys::units::literals; // support unit literals like 5_m; - -int main() -{ - CoordinateSystem root; - - Point const r0(root, {0_m, 0_m, 0_m}); - auto const omegaC = 2 * M_PI * 1_Hz; - Vector<speed_d> vPar(root, {0_m / second, 0_m / second, 10_cm / second}); - Vector<speed_d> vPerp(root, {1_m / second, 0_m / second, 0_m / second}); - - Helix h(r0, omegaC, vPar, vPerp); - - auto constexpr t0 = 0_s; - auto constexpr t1 = 1_s; - auto constexpr dt = 1_us; - auto constexpr n = long((t1 - t0) / dt) + 1; - - auto arr = std::make_unique<std::array<std::array<double, 4>, n>>(); - auto& positions = *arr; - - for (auto [t, i] = std::tuple{t0, 0}; t < t1; t += dt, ++i) - { - auto const r = h.GetPosition(t).GetCoordinates(); - - positions[i][0] = t / 1_s; - positions[i][1] = r[0] / 1_m; - positions[i][2] = r[1] / 1_m; - positions[i][3] = r[2] / 1_m; - } - - std::cout << positions[n-2][0] << " " << positions[n-2][1] << " " << positions[n-2][2] << " " << positions[n-2][3] << std::endl; - - return EXIT_SUCCESS; +using namespace corsika::geometry; +using namespace corsika::units; + +int main() { + + CoordinateSystem root; + + Point const r0(root, {0_m, 0_m, 0_m}); + auto const omegaC = 2 * M_PI * 1_Hz; + Vector<speed_d> vPar(root, {0_m / second, 0_m / second, 10_cm / second}); + Vector<speed_d> vPerp(root, {1_m / second, 0_m / second, 0_m / second}); + + Helix h(r0, omegaC, vPar, vPerp); + + auto constexpr t0 = 0_s; + auto constexpr t1 = 1_s; + auto constexpr dt = 1_us; + auto constexpr n = long((t1 - t0) / dt) + 1; + + auto arr = std::make_unique<std::array<std::array<double, 4>, n>>(); + auto& positions = *arr; + + for (auto [t, i] = std::tuple{t0, 0}; t < t1; t += dt, ++i) { + auto const r = h.GetPosition(t).GetCoordinates(); + + positions[i][0] = t / 1_s; + positions[i][1] = r[0] / 1_m; + positions[i][2] = r[1] / 1_m; + positions[i][3] = r[2] / 1_m; + } + + std::cout << positions[n - 2][0] << " " << positions[n - 2][1] << " " + << positions[n - 2][2] << " " << positions[n - 2][3] << std::endl; + + return EXIT_SUCCESS; } diff --git a/Documentation/Examples/logger_example.cc b/Documentation/Examples/logger_example.cc index d8212e6612dc2f45eb73990571d5d4c4125562ce..94fcb588a2b4ebe225dd7750a17d5beb66abda38 100644 --- a/Documentation/Examples/logger_example.cc +++ b/Documentation/Examples/logger_example.cc @@ -1,45 +1,46 @@ -#include <Logging/Logger.h> - -#include <string> -#include <iostream> -#include <fstream> - +#include <corsika/logging/Logger.h> #include <boost/format.hpp> +#include <fstream> +#include <iostream> +#include <string> using namespace std; +using namespace corsika::logging; -int -main() -{ +int main() { { cout << "writing to \"another.log\"" << endl; ofstream logfile("another.log"); - logger::sink::SinkStream unbuffered_sink(logfile); - logger::sink::BufferedSinkStream sink(logfile, logger::sink::StdBuffer(10000)); - logger::Logger<logger::MessageOn, logger::sink::BufferedSinkStream> info("\033[32m", "info", sink); - logger::Logger<logger::MessageOn, logger::sink::BufferedSinkStream> err("\033[31m", "error", sink); - //logger<ostream,messageconst,StdBuffer> info(std::cout, StdBuffer(10000)); - + sink::SinkStream unbuffered_sink(logfile); + sink::BufferedSinkStream sink(logfile, sink::StdBuffer(10000)); + Logger<MessageOn, sink::BufferedSinkStream> info("\033[32m", "info", sink); + Logger<MessageOn, sink::BufferedSinkStream> err("\033[31m", "error", sink); + // logger<ostream,messageconst,StdBuffer> info(std::cout, StdBuffer(10000)); + /* Logging& logs = Logging::GetInstance(); logs.AddLogger<>("info", info); - auto& log_1 = logs.GetLogger("info"); // no so useful, since type of log_1 is std::any + auto& log_1 = logs.GetLogger("info"); // no so useful, since type of log_1 is + std::any */ - - for (int i=0; i<100000; ++i) { - LOG(info, "irgendwas"," ", string("and more")," ", boost::format("error: %i message: %s. done."), i, "stupido"); + + for (int i = 0; i < 100000; ++i) { + LOG(info, "irgendwas", " ", string("and more"), " ", + boost::format("error: %i message: %s. done."), i, "stupido"); LOG(err, "Fehler"); } } - + { - logger::sink::NoSink off; - logger::Logger<logger::MessageOff> info("", "", off); - - for (int i=0; i<100000; ++i) { - LOG(info, "irgendwas", string("and more"), boost::format("error: %i message: %s. done."), i, "stupido", "a-number:", 8.99, "ENDE" ); + sink::NoSink off; + Logger<MessageOff> info("", "", off); + + for (int i = 0; i < 100000; ++i) { + LOG(info, "irgendwas", string("and more"), + boost::format("error: %i message: %s. done."), i, "stupido", "a-number:", 8.99, + "ENDE"); } } - + return 0; } diff --git a/Documentation/Examples/stack_example.cc b/Documentation/Examples/stack_example.cc index d27467f56999a113442befdb644ee3e8893cda63..d54f204341bb878bdc663690c515be9c201b4bc3 100644 --- a/Documentation/Examples/stack_example.cc +++ b/Documentation/Examples/stack_example.cc @@ -1,34 +1,36 @@ -#include <ParticleStack/StackOne.h> - -#include <iostream> +#include <corsika/particles/ParticleProperties.h> +#include <corsika/stack/super_stupid/SuperStupidStack.h> #include <iomanip> +#include <iostream> using namespace std; +// using namespace corsika::literals; +// using namespace corsika::io; + +using namespace corsika::units; +using namespace corsika::stack; -void fill(stack::StackOne& s) -{ - for (int i=0; i<11; ++i) { +void fill(corsika::stack::super_stupid::SuperStupidStack& s) { + for (int i = 0; i < 11; ++i) { auto p = s.NewParticle(); - p.SetId(i); - p.SetEnergy(1.5*i); + p.SetId(corsika::particles::Code::Electron); + p.SetEnergy(1.5_GeV * i); } } -void -read(stack::StackOne& s) -{ +void read(corsika::stack::super_stupid::SuperStupidStack& s) { cout << "found Stack with " << s.GetSize() << " particles. " << endl; - double Etot = 0; + EnergyType Etot; for (auto p : s) { Etot += p.GetEnergy(); + cout << "particle: " << p.GetId() << " with " << p.GetEnergy() / 1_GeV << " GeV" + << endl; } - cout << "Etot=" << Etot << endl; + cout << "Etot=" << Etot << " = " << Etot / 1_GeV << " GeV" << endl; } -int -main() -{ - stack::StackOne s; +int main() { + corsika::stack::super_stupid::SuperStupidStack s; fill(s); read(s); return 0; diff --git a/Documentation/Examples/staticsequence_example.cc b/Documentation/Examples/staticsequence_example.cc index c8bfe54d61d216baddb9a5ec9d2199cacae369f7..a28150c23d3057a7a10bc94fcb3aeb17d812dd70 100644 --- a/Documentation/Examples/staticsequence_example.cc +++ b/Documentation/Examples/staticsequence_example.cc @@ -1,99 +1,83 @@ -#include <iostream> -#include <iomanip> #include <array> +#include <iomanip> +#include <iostream> -#include <ProcessSequence/ProcessSequence.h> +#include <corsika/process/ProcessSequence.h> using namespace std; +using namespace corsika::process; - -class Process1 : public processes::Base <Process1> -{ +class Process1 : public corsika::process::BaseProcess<Process1> { public: Process1() {} - template<typename D> void DoContinuous(D& d) const { - for (int i=0; i<10; ++i) d.p[i] += 1; + template <typename D> + void DoContinuous(D& d) const { + for (int i = 0; i < 10; ++i) d.p[i] += 1; } }; -class Process2 : public processes::Base <Process2> -{ +class Process2 : public corsika::process::BaseProcess<Process2> { public: Process2() {} - - template<typename D> inline void DoContinuous(D& d) const { - //for (int i=0; i<10; ++i) d.p[i] *= 2; + + template <typename D> + inline void DoContinuous(D& d) const { + // for (int i=0; i<10; ++i) d.p[i] *= 2; } }; -class Process3 : public processes::Base <Process3> -{ +class Process3 : public BaseProcess<Process3> { public: - //Process3(const int v) :fV(v) {} + // Process3(const int v) :fV(v) {} Process3() {} - - template<typename D> inline void DoContinuous(D& d) const { - //for (int i=0; i<10; ++i) d.p[i] += fV; + + template <typename D> + inline void DoContinuous(D& d) const { + // for (int i=0; i<10; ++i) d.p[i] += fV; } - + private: - //int fV; + // int fV; }; -class Process4 : public processes::Base <Process4> -{ +class Process4 : public BaseProcess<Process4> { public: - //Process4(const int v) : fV(v) {} - Process4() {} - template<typename D> inline void DoContinuous(D& d) const { - //for (int i=0; i<10; ++i) d.p[i] /= fV; + // Process4(const int v) : fV(v) {} + Process4() {} + template <typename D> + inline void DoContinuous(D& d) const { + // for (int i=0; i<10; ++i) d.p[i] /= fV; } - + private: - //int fV; + // int fV; }; - class Data { public: std::array<double, 10> p{{0.}}; }; - - -void -modular() -{ +void modular() { Data d0; - + Process1 m1; Process2 m2; Process3 m3; Process4 m4; - - const auto sequence = m1 + m2 + m3 + m4; - + + const auto sequence = m1 + m2 + m3 + m4; + const int n = 100000000; - for (int i=0; i<n; ++i) { - sequence.DoContinuous(d0); - } - + for (int i = 0; i < n; ++i) { sequence.DoContinuous(d0); } + double s = 0; - for (int i=0; i<10; ++i) { - s += d0.p[i]; - } - + for (int i = 0; i < 10; ++i) { s += d0.p[i]; } + cout << scientific << " v=" << s << " n=" << n << endl; } - - - -int -main() -{ +int main() { modular(); return 0; } - - diff --git a/Framework/CMakeLists.txt b/Framework/CMakeLists.txt index 9de69dc03263d359c3c9fd99aa26eb64cc6f122a..2fdb7c85aded8c665c56e249786444dcd274b186 100644 --- a/Framework/CMakeLists.txt +++ b/Framework/CMakeLists.txt @@ -4,5 +4,4 @@ add_subdirectory (Geometry) add_subdirectory (Particles) add_subdirectory (Logging) add_subdirectory (StackInterface) -add_subdirectory (ParticleStack) add_subdirectory (ProcessSequence) diff --git a/Framework/Cascade/Cascade.cc b/Framework/Cascade/Cascade.cc index ffd4607039c933165ade08c624b332f5fddbf9b2..f4496f536634ba9b82a0a740199ad38b2f2054e4 100644 --- a/Framework/Cascade/Cascade.cc +++ b/Framework/Cascade/Cascade.cc @@ -2,26 +2,20 @@ namespace cascade; -template<typename Sequence, typename Trajectory> -void -Cascade::Cascade() -{ +template <typename Sequence, typename Trajectory> +void Cascade::Cascade() { + kkk; + kk; } - -template<typename Sequence, typename Trajectory> -void -Cascade::Init() -{ +template <typename Sequence, typename Trajectory> +void Cascade::Init() { fStack.Init(); fProcesseList.Init(); } - -template<typename Sequence, typename Trajectory> -void -Cascade::Run() -{ +template <typename Sequence, typename Trajectory> +void Cascade::Run() { if (!fStack.IsEmpty()) { if (!fStack.IsEmpty()) { Particle& p = fStack.GetNextParticle(); @@ -33,14 +27,10 @@ Cascade::Run() } } - -template<typename Sequence, typename Trajectory> -void -Cascade::Step(Particle& particle) -{ +template <typename Sequence, typename Trajectory> +void Cascade::Step(Particle& particle) { double nextStep = fProcesseList.MinStepLength(particle); Trajectory trajectory = fProcesseList.Transport(particle, nextStep); sequence.DoContinuous(particle, trajectory); sequence.DoDiscrete(particle); } - diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 8513a95f809632f495548666ff25d24edbfe24ae..a8af260b512c91a62d12a776fb74b3ecf5520917 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -3,9 +3,9 @@ namespace cascade { - template<typename Processes, typename Trajectory, typename Stack> + template <typename Processes, typename Trajectory, typename Stack> class Cascade { - + public: Cascade(); @@ -16,9 +16,8 @@ namespace cascade { private: Stack fStack; Processes fProcesseList; - }; - -} + +} // namespace cascade #endif diff --git a/Framework/Cascade/Step.cc b/Framework/Cascade/Step.cc index 5a04f06a60b25f0a35ec27e93ba917f478fac104..60349c748f8438857a5dc4ffda338063a180634f 100644 --- a/Framework/Cascade/Step.cc +++ b/Framework/Cascade/Step.cc @@ -2,13 +2,9 @@ namespace cascade; - -void -Cascade::Step(auto& sequence, Particle& particle) -{ +void Cascade::Step(auto& sequence, Particle& particle) { double nextStep = sequence.MinStepLength(particle); Trajectory trajectory = sequence.Transport(particle, nextStep); sequence.DoContinuous(particle, trajectory); sequence.DoDiscrete(particle); } - diff --git a/Framework/Geometry/BaseVector.h b/Framework/Geometry/BaseVector.h index 02edb6fa8199e0d70a0eea1b326ca6f07b8bd5b3..f2079e91cc953b34499c2146da68c4cf88dfed8d 100644 --- a/Framework/Geometry/BaseVector.h +++ b/Framework/Geometry/BaseVector.h @@ -1,25 +1,27 @@ #ifndef _include_BASEVECTOR_H_ #define _include_BASEVECTOR_H_ -#include <Geometry/QuantityVector.h> -#include <Geometry/CoordinateSystem.h> +#include <corsika/geometry/CoordinateSystem.h> +#include <corsika/geometry/QuantityVector.h> -/*! - * Common base class for Vector and Point. Currently it does basically nothing. - */ +namespace corsika::geometry { -template <typename dim> -class BaseVector -{ -protected: + /*! + * Common base class for Vector and Point. Currently it does basically nothing. + */ + + template <typename dim> + class BaseVector { + protected: QuantityVector<dim> qVector; CoordinateSystem const* cs; - -public: - BaseVector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector) : - qVector(pQVector), cs(&pCS) - { - } -}; + + public: + BaseVector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector) + : qVector(pQVector) + , cs(&pCS) {} + }; + +} // namespace corsika::geometry #endif diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt index 16dce73058acacd30e22672ff8eb59acd2b3ab88..c55295c897dcf0253293b3342d5042aa77cc8e26 100644 --- a/Framework/Geometry/CMakeLists.txt +++ b/Framework/Geometry/CMakeLists.txt @@ -1,31 +1,68 @@ -set (GEOMETRY_SOURCES CoordinateSystem.cc) -set (GEOMETRY_HEADERS Vector.h Point.h Sphere.h CoordinateSystem.h Helix.h) +set ( + GEOMETRY_SOURCES + CoordinateSystem.cc + ) + +set ( + GEOMETRY_HEADERS + Vector.h + Point.h + Sphere.h + CoordinateSystem.h + Helix.h + BaseVector.h + QuantityVector.h + ) + +set ( + GEOMETRY_NAMESPACE + corsika/geometry + ) add_library (CORSIKAgeometry STATIC ${GEOMETRY_SOURCES}) +CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAgeometry ${GEOMETRY_NAMESPACE} ${GEOMETRY_HEADERS}) -set_target_properties (CORSIKAgeometry PROPERTIES VERSION ${PROJECT_VERSION}) -set_target_properties (CORSIKAgeometry PROPERTIES SOVERSION 1) +set_target_properties ( + CORSIKAgeometry + PROPERTIES + VERSION ${PROJECT_VERSION} + SOVERSION 1 + PUBLIC_HEADER "${GEOMETRY_HEADERS}" + ) -set_target_properties (CORSIKAgeometry PROPERTIES PUBLIC_HEADER "${GEOMETRY_HEADERS}") - -# target dependencies on other libraries (also header only) -target_link_libraries (CORSIKAgeometry CORSIKAunits) +# target dependencies on other libraries (also the header onlys) +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> - ) +target_include_directories ( + CORSIKAgeometry + PRIVATE ${EIGEN3_INCLUDE_DIR} + INTERFACE ${EIGEN3_INCLUDE_DIR} + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> + $<INSTALL_INTERFACE:include/include> + ) -install (TARGETS CORSIKAgeometry - LIBRARY DESTINATION lib - ARCHIVE DESTINATION lib - PUBLIC_HEADER DESTINATION include/Geometry) +install ( + TARGETS CORSIKAgeometry + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib + PUBLIC_HEADER DESTINATION include/${GEOMETRY_NAMESPACE} + ) -# code testing +# -------------------- +# code unit 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) + +target_link_libraries ( + testGeometry + CORSIKAgeometry + CORSIKAunits + CORSIKAthirdparty # for catch2 + ) + +add_test (NAME testGeometry COMMAND testGeometry) diff --git a/Framework/Geometry/CoordinateSystem.cc b/Framework/Geometry/CoordinateSystem.cc index 9c9f97e39f4088a5b2adcc2c1ac1ac24e69f33d1..db0dea292aa5757f636c374da7d2c8155e34287f 100644 --- a/Framework/Geometry/CoordinateSystem.cc +++ b/Framework/Geometry/CoordinateSystem.cc @@ -1,54 +1,45 @@ -#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; +#include <corsika/geometry/CoordinateSystem.h> + +using namespace corsika::geometry; + +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; } diff --git a/Framework/Geometry/CoordinateSystem.h b/Framework/Geometry/CoordinateSystem.h index 238054ca7db7698f302efcd09cc90abc14a688b5..b8a7afc25be0f3ab800fd7d4f4a69ae9646578e1 100644 --- a/Framework/Geometry/CoordinateSystem.h +++ b/Framework/Geometry/CoordinateSystem.h @@ -1,72 +1,64 @@ #ifndef _include_COORDINATESYSTEM_H_ #define _include_COORDINATESYSTEM_H_ -#include <Geometry/QuantityVector.h> -#include <Units/PhysicalUnits.h> - +#include <corsika/geometry/QuantityVector.h> +#include <corsika/units/PhysicalUnits.h> #include <Eigen/Dense> typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform; typedef Eigen::Translation<double, 3> EigenTranslation; -class CoordinateSystem -{ +namespace corsika::geometry { + + using corsika::units::length_d; + + 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); + + 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 rotate(QuantityVector<phys::units::length_d> axis, double angle) const - { - EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())}; - - return CoordinateSystem(*this, rotation); + + auto translate(QuantityVector<length_d> vector) const { + EigenTransform const translation{EigenTranslation(vector.eVector)}; + + return CoordinateSystem(*this, translation); } - - 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 rotate(QuantityVector<phys::units::length_d> axis, double angle) const { + EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())}; + + return CoordinateSystem(*this, rotation); } - - auto const* GetReference() const - { - return reference; + + 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& GetTransform() const - { - return transf; - } -}; + + auto const* GetReference() const { return reference; } + + auto const& GetTransform() const { return transf; } + }; + +} // namespace corsika::geometry #endif diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h index 10480c00fca7bec80153c76c52410feb626e3922..c392356d46e64d63d156ca76f848845e916cde36 100644 --- a/Framework/Geometry/Helix.h +++ b/Framework/Geometry/Helix.h @@ -1,41 +1,48 @@ #ifndef _include_HELIX_H_ #define _include_HELIX_H_ -#include <Geometry/Vector.h> -#include <Geometry/Point.h> -#include <Units/PhysicalUnits.h> - +#include <corsika/geometry/Point.h> +#include <corsika/geometry/Vector.h> +#include <corsika/units/PhysicalUnits.h> #include <cmath> -class Helix // TODO: inherit from to-be-implemented "Trajectory" -{ - using SpeedVec = Vector<Speed::dimension_type>; - +namespace corsika::geometry { + + using corsika::units::frequency_d; + using corsika::units::FrequencyType; + using corsika::units::quantity; + using corsika::units::SpeedType; + using corsika::units::TimeType; + + class Helix // TODO: inherit from to-be-implemented "Trajectory" + { + using SpeedVec = Vector<SpeedType::dimension_type>; + Point const r0; - Frequency const omegaC; + FrequencyType const omegaC; SpeedVec const vPar; SpeedVec vPerp, uPerp; - - Length const radius; - -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(vPerp.cross(vPar.normalized())), - radius(pvPar.norm() / abs(pOmegaC)) - { - - } - - auto GetPosition(Time t) const - { - return r0 + vPar * t + (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC; - } - - auto GetRadius() const - { - return radius; + + LengthType const radius; + + public: + Helix(Point const& pR0, quantity<frequency_d> pOmegaC, SpeedVec const& pvPar, + SpeedVec const& pvPerp) + : r0(pR0) + , omegaC(pOmegaC) + , vPar(pvPar) + , vPerp(pvPerp) + , uPerp(vPerp.cross(vPar.normalized())) + , radius(pvPar.norm() / abs(pOmegaC)) {} + + auto GetPosition(TimeType t) const { + return r0 + vPar * t + + (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC; } -}; + + auto GetRadius() const { return radius; } + }; + +} // namespace corsika::geometry #endif diff --git a/Framework/Geometry/LineTrajectory.h b/Framework/Geometry/LineTrajectory.h index f9efee6e67d523a5ecaadaadfb40f76fa72acb5f..4d56a598f58b8a5073b583654b1377c09d4d92f6 100644 --- a/Framework/Geometry/LineTrajectory.h +++ b/Framework/Geometry/LineTrajectory.h @@ -1,26 +1,26 @@ #ifndef _include_LINETRAJECTORY_H #define _include_LINETRAJECTORY_H -#include <Framework/Geometry/Point.h> -#include <Framework/Geometry/Vector.h> #include <Units/PhysicalUnits.h> +#include <corsika/Point.h> +#include <corsika/Vector.h> -class LineTrajectory // TODO: inherit from Trajectory -{ +namesapce corsika { + + class LineTrajectory // TODO: inherit from Trajectory + { using SpeedVec = Vector<Speed::dimension_type>; - + Point const r0; SpeedVec const v0; - - LineTrajectory(Point const& pR0, SpeedVec const& pV0) : - r0(r0), v0(pV0) - { - } - - auto GetPosition(Time t) const - { - return r0 + v0 * t; - } -}; + + LineTrajectory(Point const& pR0, SpeedVec const& pV0) + : r0(r0) + , v0(pV0) {} + + auto GetPosition(Time t) const { return r0 + v0 * t; } + }; + +} // end namesapce #endif diff --git a/Framework/Geometry/Point.h b/Framework/Geometry/Point.h index f5efef09d77a24ee0587fb42f5be1d24150958f9..108fe87f4aadd5339d93160a0b5f0efaab80199a 100644 --- a/Framework/Geometry/Point.h +++ b/Framework/Geometry/Point.h @@ -1,68 +1,65 @@ #ifndef _include_POINT_H_ #define _include_POINT_H_ -#include <Geometry/BaseVector.h> -#include <Geometry/QuantityVector.h> -#include <Geometry/Vector.h> -#include <Units/PhysicalUnits.h> +#include <corsika/geometry/BaseVector.h> +#include <corsika/geometry/QuantityVector.h> +#include <corsika/geometry/Vector.h> +#include <corsika/units/PhysicalUnits.h> -/*! - * A Point represents a point in position space. It is defined by its - * coordinates with respect to some CoordinateSystem. - */ -class Point : public BaseVector<phys::units::length_d> -{ -public: - Point(CoordinateSystem const& pCS, QuantityVector<phys::units::length_d> pQVector) : - BaseVector<phys::units::length_d>(pCS, pQVector) - { - } +namespace corsika::geometry { - Point(CoordinateSystem const& cs, Length x, Length y, Length z) : - BaseVector<phys::units::length_d>(cs, {x, y, z}) - { - } - - auto GetCoordinates() const - { + using corsika::units::length_d; + using corsika::units::LengthType; + + /*! + * A Point represents a point in position space. It is defined by its + * coordinates with respect to some CoordinateSystem. + */ + class Point : public BaseVector<phys::units::length_d> { + public: + Point(CoordinateSystem const& pCS, QuantityVector<phys::units::length_d> pQVector) + : BaseVector<phys::units::length_d>(pCS, pQVector) {} + + Point(CoordinateSystem const& cs, LengthType x, LengthType y, LengthType 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); + } } - - 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); - } - } - + /*! * transforms the Point into another CoordinateSystem by changing its * coordinates interally */ - void rebase(CoordinateSystem const& pCS) - { - BaseVector<phys::units::length_d>::qVector = GetCoordinates(pCS); - BaseVector<phys::units::length_d>::cs = &pCS; + 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)); + + 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)); } - + /*! * returns the distance Vector between two points */ - 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)); + 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)); } -}; + }; + +} // namespace corsika::geometry #endif diff --git a/Framework/Geometry/QuantityVector.h b/Framework/Geometry/QuantityVector.h index 4f9209f8e2ee76bb6549f892712e7cf2cea563fc..c5af419ee4df97ec731a7fb0c95123fc2c843de5 100644 --- a/Framework/Geometry/QuantityVector.h +++ b/Framework/Geometry/QuantityVector.h @@ -1,139 +1,118 @@ #ifndef _include_QUANTITYVECTOR_H_ #define _include_QUANTITYVECTOR_H_ -#include <Units/PhysicalUnits.h> +#include <corsika/units/PhysicalUnits.h> #include <Eigen/Dense> + #include <iostream> #include <utility> -/*! - * A QuantityVector is a three-component container based on Eigen::Vector3d - * with a phys::units::dimension. Arithmethic operators are defined that - * propagate the dimensions by dimensional analysis. - */ +namespace corsika { -template <typename dim> -class QuantityVector -{ -protected: - using Quantity = phys::units::quantity<dim, double>; //< the phys::units::quantity corresponding to the dimension - -public: + /*! + * A QuantityVector is a three-component container based on Eigen::Vector3d + * with a phys::units::dimension. Arithmethic operators are defined that + * propagate the dimensions by dimensional analysis. + */ + + template <typename dim> + class QuantityVector { + protected: + using Quantity = phys::units::quantity<dim, double>; //< the phys::units::quantity + // corresponding to the dimension + + public: Eigen::Vector3d eVector; //!< the actual container where the raw numbers are stored - + typedef dim dimension; //!< should be a phys::units::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]); + 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 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 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); } - - 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()); - } + 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) - { - eVector /= p; - return *this; + auto operator/(phys::units::quantity<ScalarDim, double> const p) const { + return (*this) * (1 / p); } - - auto& operator*=(double const p) - { - eVector *= p; - return *this; - } - - auto& operator+=(QuantityVector<dim> const& pQVec) - { - eVector += pQVec.eVector; - return *this; + + 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) { + eVector /= p; + return *this; } - - auto& operator-=(QuantityVector<dim> const& pQVec) - { - eVector -= pQVec.eVector; - return *this; + + auto& operator*=(double const p) { + eVector *= p; + return *this; } - - auto& operator-() const - { - return QuantityVector<dim>(-eVector); + + auto& operator+=(QuantityVector<dim> const& pQVec) { + eVector += pQVec.eVector; + return *this; } - - auto normalized() const - { - return (*this) * (1 / norm()); + + 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()); } + }; + +} // end namespace corsika 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; +auto& operator<<(std::ostream& os, corsika::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 diff --git a/Framework/Geometry/Sphere.h b/Framework/Geometry/Sphere.h index cc87e1088161c0ea6800bad116c154e23deff015..552dfde6f1341a7c0bbfb3903e27c8fbab5e1e05 100644 --- a/Framework/Geometry/Sphere.h +++ b/Framework/Geometry/Sphere.h @@ -1,26 +1,25 @@ #ifndef _include_SPHERE_H_ #define _include_SPHERE_H_ -#include <Geometry/Point.h> -#include <Units/PhysicalUnits.h> +#include <corsika/geometry/Point.h> +#include <corsika/units/PhysicalUnits.h> -class Sphere -{ +namespace corsika::geometry { + + class Sphere { 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(); + LengthType const radius; + + public: + Sphere(Point const& pCenter, LengthType const pRadius) + : center(pCenter) + , radius(pRadius) {} + + auto isInside(Point const& p) const { + return radius * radius > (center - p).squaredNorm(); } + }; -}; +} // namespace corsika::geometry #endif diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h index 3c48f9df8d52df91ab95be553d113a619059a766..fcebc10f73a12ad37494a59d32614c88e9697a31 100644 --- a/Framework/Geometry/Vector.h +++ b/Framework/Geometry/Vector.h @@ -1,89 +1,75 @@ #ifndef _include_VECTOR_H_ #define _include_VECTOR_H_ -#include <Geometry/BaseVector.h> -#include <Geometry/QuantityVector.h> -#include <Units/PhysicalUnits.h> +#include <corsika/geometry/BaseVector.h> +#include <corsika/geometry/QuantityVector.h> + +#include <corsika/units/PhysicalUnits.h> /*! * A Vector represents a 3-vector in Euclidean space. It is defined by components * given in a specific CoordinateSystem. It has a physical dimension ("unit") * as part of its type, so you cannot mix up e.g. electric with magnetic fields * (but you could calculate their cross-product to get an energy flux vector). - * + * * When transforming coordinate systems, a Vector is subject to the rotational * part only and invariant under translations. */ -template <typename dim> -class Vector : public BaseVector<dim> -{ +namespace corsika::geometry { + + 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)) - { - } - + 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)) {} + /*! * returns a QuantityVector with the components given in the "home" * CoordinateSystem of the Vector */ - auto GetComponents() const - { - return BaseVector<dim>::qVector; - } - + auto GetComponents() const { return BaseVector<dim>::qVector; } + /*! * returns a QuantityVector with the components given in an arbitrary * CoordinateSystem */ - 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); - } + 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); + } } - + /*! * transforms the Vector into another CoordinateSystem by changing * its components internally */ - void rebase(CoordinateSystem const& pCS) - { - BaseVector<dim>::qVector = GetComponents(pCS); - BaseVector<dim>::cs = &pCS; + void rebase(CoordinateSystem const& pCS) { + BaseVector<dim>::qVector = GetComponents(pCS); + BaseVector<dim>::cs = &pCS; } - + /*! * returns the norm/length of the Vector. Before using this method, * think about whether squaredNorm() might be cheaper for your computation. */ - auto norm() const - { - return BaseVector<dim>::qVector.norm(); - } - + auto norm() const { return BaseVector<dim>::qVector.norm(); } + /*! * returns the squared norm of the Vector. Before using this method, * think about whether norm() might be cheaper for your computation. */ - auto squaredNorm() const - { - return BaseVector<dim>::qVector.squaredNorm(); - } - + auto squaredNorm() const { return BaseVector<dim>::qVector.squaredNorm(); } + /*! * returns a Vector \f$ \vec{v}_{\parallel} \f$ which is the parallel projection * of this vector \f$ \vec{v}_1 \f$ along another Vector \f$ \vec{v}_2 \f$ given by @@ -92,112 +78,100 @@ public: * \f] */ 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()))); + 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 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(*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-(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; + + 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); - } + 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/(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/(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& operator-=(Vector<dim> const& pVec) { + BaseVector<dim>::qVector -= pVec.GetComponents(*BaseVector<dim>::cs); + return *this; } - - auto normalized() const - { - return (*this) * (1 / norm()); + + 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); - } - } - -}; + 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); + } + } + }; + +} // namespace corsika::geometry #endif diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc index 95baf8d2a9405978ceb4c77658326cd668be10f5..af824df1c37a31182c8a31dd2d2fbe7fb64afd6f 100644 --- a/Framework/Geometry/testGeometry.cc +++ b/Framework/Geometry/testGeometry.cc @@ -1,107 +1,116 @@ -#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one cpp file +#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 <Geometry/CoordinateSystem.h> #include <Geometry/Point.h> #include <Geometry/Sphere.h> +#include <corsika/units/PhysicalUnits.h> #include <cmath> using namespace phys::units; using namespace phys::units::literals; -TEST_CASE( "transformations between CoordinateSystems") -{ - CoordinateSystem rootCS; - - REQUIRE ( CoordinateSystem::GetTransformation(rootCS, rootCS).isApprox(EigenTransform::Identity()) ); - - QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m}; - Point p1(rootCS, coordinates); - - QuantityVector<magnetic_flux_density_d> components{1.*tesla, 0.*tesla, 0.*tesla}; - Vector<magnetic_flux_density_d> v1(rootCS, components); - - REQUIRE ( (p1.GetCoordinates() - coordinates).norm().magnitude() == Approx(0) ); - REQUIRE ( (p1.GetCoordinates(rootCS) - coordinates).norm().magnitude() == Approx(0) ); - - SECTION ( "unconnected CoordinateSystems" ) { - CoordinateSystem rootCS2; - REQUIRE_THROWS(CoordinateSystem::GetTransformation(rootCS, rootCS2)); - } - - SECTION( "translations" ) { - QuantityVector<length_d> const translationVector{0_m, 4_m, 0_m}; - - CoordinateSystem translatedCS = rootCS.translate(translationVector); - - REQUIRE (translatedCS.GetReference() == &rootCS); - - REQUIRE ( (p1.GetCoordinates(translatedCS) + translationVector).norm().magnitude() == Approx(0) ); - - // Vectors are not subject to translations - REQUIRE ((v1.GetComponents(rootCS) - v1.GetComponents(translatedCS)).norm().magnitude() == Approx( 0 )); - - Point p2(translatedCS, {0_m, 0_m, 0_m}); - REQUIRE ( ((p2 - p1).GetComponents() - translationVector).norm().magnitude() == Approx( 0 )); - } - - SECTION( "multiple translations" ) { - QuantityVector<length_d> const tv1{0_m, 5_m, 0_m}; - CoordinateSystem cs2 = rootCS.translate(tv1); - - QuantityVector<length_d> const tv2{3_m, 0_m, 0_m}; - CoordinateSystem cs3 = rootCS.translate(tv2); - - QuantityVector<length_d> const tv3{0_m, 0_m, 2_m}; - CoordinateSystem cs4 = cs3.translate(tv3); - - REQUIRE( cs4.GetReference()->GetReference() == &rootCS ); - - REQUIRE( CoordinateSystem::GetTransformation(cs3, cs2).isApprox(rootCS.translate({3_m, -5_m, 0_m}).GetTransform()) ); - REQUIRE( CoordinateSystem::GetTransformation(cs2, cs3).isApprox(rootCS.translate({-3_m, +5_m, 0_m}).GetTransform()) ); - } - - SECTION( "rotations" ) { - QuantityVector<length_d> const axis{0_m, 0_m, 1_km}; - double const angle = 90. / 180. * M_PI; - - CoordinateSystem rotatedCS = rootCS.rotate(axis, angle); - REQUIRE (rotatedCS.GetReference() == &rootCS); - - REQUIRE( v1.GetComponents(rotatedCS)[1].magnitude() == Approx((-1. * tesla).magnitude()) ); - - // vector norm invariant under rotation - REQUIRE( v1.GetComponents(rotatedCS).norm().magnitude() == Approx( v1.GetComponents(rootCS).norm().magnitude() ) ); - } - - SECTION("multiple rotations") { - QuantityVector<length_d> const zAxis{0_m, 0_m, 1_km}; - QuantityVector<length_d> const yAxis{0_m, 7_nm, 0_m}; - QuantityVector<length_d> const xAxis{2_m, 0_nm, 0_m}; - - double const angle = 90. / 180. * M_PI; - - CoordinateSystem rotated1 = rootCS.rotate(zAxis, angle); - CoordinateSystem rotated2 = rotated1.rotate(yAxis, angle); - CoordinateSystem rotated3 = rotated2.rotate(zAxis, -angle); - - CoordinateSystem combined = rootCS.rotate(xAxis, -angle); - - auto comp1 = v1.GetComponents(rootCS); - auto comp3 = v1.GetComponents(combined); - REQUIRE ((comp1 - comp3).norm().magnitude() == Approx(0)); - } +TEST_CASE("transformations between CoordinateSystems") { + CoordinateSystem rootCS; + + REQUIRE(CoordinateSystem::GetTransformation(rootCS, rootCS) + .isApprox(EigenTransform::Identity())); + + QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m}; + Point p1(rootCS, coordinates); + + QuantityVector<magnetic_flux_density_d> components{1. * tesla, 0. * tesla, 0. * tesla}; + Vector<magnetic_flux_density_d> v1(rootCS, components); + + REQUIRE((p1.GetCoordinates() - coordinates).norm().magnitude() == Approx(0)); + REQUIRE((p1.GetCoordinates(rootCS) - coordinates).norm().magnitude() == Approx(0)); + + SECTION("unconnected CoordinateSystems") { + CoordinateSystem rootCS2; + REQUIRE_THROWS(CoordinateSystem::GetTransformation(rootCS, rootCS2)); + } + + SECTION("translations") { + QuantityVector<length_d> const translationVector{0_m, 4_m, 0_m}; + + CoordinateSystem translatedCS = rootCS.translate(translationVector); + + REQUIRE(translatedCS.GetReference() == &rootCS); + + REQUIRE((p1.GetCoordinates(translatedCS) + translationVector).norm().magnitude() == + Approx(0)); + + // Vectors are not subject to translations + REQUIRE( + (v1.GetComponents(rootCS) - v1.GetComponents(translatedCS)).norm().magnitude() == + Approx(0)); + + Point p2(translatedCS, {0_m, 0_m, 0_m}); + REQUIRE(((p2 - p1).GetComponents() - translationVector).norm().magnitude() == + Approx(0)); + } + + SECTION("multiple translations") { + QuantityVector<length_d> const tv1{0_m, 5_m, 0_m}; + CoordinateSystem cs2 = rootCS.translate(tv1); + + QuantityVector<length_d> const tv2{3_m, 0_m, 0_m}; + CoordinateSystem cs3 = rootCS.translate(tv2); + + QuantityVector<length_d> const tv3{0_m, 0_m, 2_m}; + CoordinateSystem cs4 = cs3.translate(tv3); + + REQUIRE(cs4.GetReference()->GetReference() == &rootCS); + + REQUIRE(CoordinateSystem::GetTransformation(cs3, cs2).isApprox( + rootCS.translate({3_m, -5_m, 0_m}).GetTransform())); + REQUIRE(CoordinateSystem::GetTransformation(cs2, cs3).isApprox( + rootCS.translate({-3_m, +5_m, 0_m}).GetTransform())); + } + + SECTION("rotations") { + QuantityVector<length_d> const axis{0_m, 0_m, 1_km}; + double const angle = 90. / 180. * M_PI; + + CoordinateSystem rotatedCS = rootCS.rotate(axis, angle); + REQUIRE(rotatedCS.GetReference() == &rootCS); + + REQUIRE(v1.GetComponents(rotatedCS)[1].magnitude() == + Approx((-1. * tesla).magnitude())); + + // vector norm invariant under rotation + REQUIRE(v1.GetComponents(rotatedCS).norm().magnitude() == + Approx(v1.GetComponents(rootCS).norm().magnitude())); + } + + SECTION("multiple rotations") { + QuantityVector<length_d> const zAxis{0_m, 0_m, 1_km}; + QuantityVector<length_d> const yAxis{0_m, 7_nm, 0_m}; + QuantityVector<length_d> const xAxis{2_m, 0_nm, 0_m}; + + double const angle = 90. / 180. * M_PI; + + CoordinateSystem rotated1 = rootCS.rotate(zAxis, angle); + CoordinateSystem rotated2 = rotated1.rotate(yAxis, angle); + CoordinateSystem rotated3 = rotated2.rotate(zAxis, -angle); + + CoordinateSystem combined = rootCS.rotate(xAxis, -angle); + + auto comp1 = v1.GetComponents(rootCS); + auto comp3 = v1.GetComponents(combined); + REQUIRE((comp1 - comp3).norm().magnitude() == Approx(0)); + } } TEST_CASE("Sphere") { - CoordinateSystem rootCS; - Point center(rootCS, {0_m, 3_m, 4_m}); - Sphere sphere(center, 5_m); - - SECTION ("isInside") { - REQUIRE_FALSE(sphere.isInside(Point(rootCS, {100_m, 0_m, 0_m}))); - - REQUIRE(sphere.isInside(Point(rootCS, {2_m, 3_m, 4_m}))); - } + CoordinateSystem rootCS; + Point center(rootCS, {0_m, 3_m, 4_m}); + Sphere sphere(center, 5_m); + + SECTION("isInside") { + REQUIRE_FALSE(sphere.isInside(Point(rootCS, {100_m, 0_m, 0_m}))); + + REQUIRE(sphere.isInside(Point(rootCS, {2_m, 3_m, 4_m}))); + } } diff --git a/Framework/Logging/BufferedSink.h b/Framework/Logging/BufferedSink.h index 9934744dc1d069e603ec4f72f4ea25652eb97e4c..7ad36927fec13087d8875af96279013a82003a4e 100644 --- a/Framework/Logging/BufferedSink.h +++ b/Framework/Logging/BufferedSink.h @@ -1,10 +1,9 @@ #ifndef _include_BufferedSink_h_ #define _include_BufferedSink_h_ -namespace logger { +namespace corsika::logging { namespace sink { - /** Output buffer template. NoBuffer does nothingk. @@ -19,47 +18,53 @@ namespace logger { */ /** Output buffer template. StdBuffer records fSize characters in - local memeory before passing it on to further output stages. + 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; } + 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> + 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(); + 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); } - if (!fBuffer.Test(msg)) - fOutput << msg; - else - fBuffer.Add(msg); - } - void Close() { fOutput << fBuffer.GetString(); } + void Close() { fOutput << fBuffer.GetString(); } + private: - TStream& fOutput; - TBuffer fBuffer; + TStream& fOutput; + TBuffer fBuffer; }; - + typedef BufferedSink<std::ostream, StdBuffer> BufferedSinkStream; - }// end namespace -} // end namespace - + } // namespace sink +} // namespace corsika::logging + #endif diff --git a/Framework/Logging/CMakeLists.txt b/Framework/Logging/CMakeLists.txt index 581df617b38eb288e0f198d60d2cfa18c1e82a7d..d7af35a8fb5f5fb7baa743262ed9491392cc7816 100644 --- a/Framework/Logging/CMakeLists.txt +++ b/Framework/Logging/CMakeLists.txt @@ -1,18 +1,55 @@ - add_library (CORSIKAlogging INTERFACE) -target_include_directories (CORSIKAlogging INTERFACE ${PROJECT_SOURCE_DIR}/Framework) +# namespace of library -> location of header files +set ( + CORSIKAlogging_NAMESPACE + corsika/logging + ) + +# header files of this library +set ( + CORSIKAlogging_HEADERS + Logger.h + Sink.h + MessageOn.h + MessageOff.h + NoSink.h + Sink.h + BufferedSink.h + ) + +CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAlogging ${CORSIKAlogging_NAMESPACE} ${CORSIKAlogging_HEADERS}) + +# include directive for upstream code +target_include_directories ( + CORSIKAlogging + INTERFACE + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> + $<INSTALL_INTERFACE:include/> + ) + +# install library +install ( + FILES ${CORSIKAlogging_HEADERS} + DESTINATION include/${CORSIKAlogging_NAMESPACE} + ) -target_include_directories (CORSIKAlogging INTERFACE $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Framework> - $<INSTALL_INTERFACE:include/Framework> - ) +# ---------------- +# code unit testing +add_executable ( + testLogging + testLogging.cc + ) -install (FILES Logger.h Sink.h MessageOn.h MessageOff.h NoSink.h Sink.h BufferedSink.h - DESTINATION include/Logging) +target_link_libraries ( + testLogging + CORSIKAlogging + CORSIKAthirdparty # for catch2 + ) -# code testing -add_executable (testLogging testLogging.cc) -target_link_libraries (testLogging CORSIKAlogging CORSIKAthirdparty) # for catch2 -add_test (NAME testLogging COMMAND testLogging -o report.xml -r junit) +add_test ( + NAME testLogging + COMMAND testLogging + ) diff --git a/Framework/Logging/Logger.h b/Framework/Logging/Logger.h index b98654d3f224711bf1d3315564f316c63544759e..fca23d70caa28a899af91f025852df0089c10b81 100644 --- a/Framework/Logging/Logger.h +++ b/Framework/Logging/Logger.h @@ -1,66 +1,83 @@ +/** + @File Logger.h + + Everything around logfile generation and text output. + */ + #ifndef _include_logger_h_ #define _include_logger_h_ #include <iosfwd> -#include <string> #include <sstream> +#include <string> #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> - +#include <corsika/logging/BufferedSink.h> +#include <corsika/logging/MessageOff.h> +#include <corsika/logging/MessageOn.h> +#include <corsika/logging/NoSink.h> +#include <corsika/logging/Sink.h> using namespace std; using namespace boost; -/** - Everything around logfile generation and text output. -*/ +namespace corsika::logging { -namespace logger { - /** + @class 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> + 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 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> + 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 +} // namespace corsika::logging + +/** + * @def LOG(...) + * + * This is the main interface to the logging facilities. If Logger + * object are defined (e.g. log1) use as + * @example LOG(log1, "var1=", variable1int, "var2=", variabl2double) + * for arbitrary long sequence + * of arguments. This may also include boost::format objects the + * output is concatenated, if log1 is switched off at compile time, + * the whole LOG command is optimized away by the compiler. + */ -#define LOG(__LOGGER,...) \ - __LOGGER.Log(__LOGGER.GetName(), __FILE__,":", __LINE__, " (", __func__, ") -> ", ##__VA_ARGS__); +#define LOG(__LOGGER, ...) \ + __LOGGER.Log(__LOGGER.GetName(), __FILE__, ":", __LINE__, " (", __func__, ") -> ", \ + ##__VA_ARGS__); #endif - diff --git a/Framework/Logging/MessageOff.h b/Framework/Logging/MessageOff.h index ef62e2577eedae957767888156bbe04474e71c0b..d60070002cbcb25e76eb641b19597bc6c2fac752 100644 --- a/Framework/Logging/MessageOff.h +++ b/Framework/Logging/MessageOff.h @@ -1,7 +1,7 @@ #ifndef _include_MessageOff_h_ #define _include_MessageOff_h_ -namespace logger { +namespace corsika::logging { /** Helper class to ignore all arguments to MessagesOn::Message and @@ -9,11 +9,12 @@ namespace logger { */ class MessageOff { protected: - template<typename First, typename ... Strings> std::string Message(const First& arg, const Strings&... rest) { + template <typename First, typename... Strings> + std::string Message(const First& arg, const Strings&... rest) { return ""; } }; -} // end namespace +} // namespace corsika::logging #endif diff --git a/Framework/Logging/MessageOn.h b/Framework/Logging/MessageOn.h index 43b847448a5b92d3b3c48d5768db874663a9d8a1..a25928674a73ccbe1837190043a2688be4118ea1 100644 --- a/Framework/Logging/MessageOn.h +++ b/Framework/Logging/MessageOn.h @@ -1,8 +1,7 @@ #ifndef _include_MessageOn_h_ #define _include_MessageOn_h_ - -namespace logger { +namespace corsika::logging { /** Helper class to convert all input arguments of MessageOn::Message @@ -11,44 +10,51 @@ namespace logger { class MessageOn { protected: std::string Message() { return "\n"; } - - template<typename First, typename ... Strings> std::string Message(const First& arg, const Strings&... rest) { + + 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) { + + 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) { + + 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) { + + 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) { + + 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) { + 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) { + + 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 +} // namespace corsika::logging #endif diff --git a/Framework/Logging/NoSink.h b/Framework/Logging/NoSink.h index 13e37f3b9a511182f792d81922e7dfba4c62d8ef..d0ba53649f2f4654e9b19fd75c1e84445eed84f6 100644 --- a/Framework/Logging/NoSink.h +++ b/Framework/Logging/NoSink.h @@ -1,7 +1,7 @@ #ifndef _include_NoSink_h_ #define _include_NoSink_h_ -namespace logger { +namespace corsika::logging { namespace sink { @@ -10,7 +10,7 @@ namespace logger { inline void Close() {} }; - }// end namespace -} // end namespace - + } // namespace sink +} // namespace corsika::logging + #endif diff --git a/Framework/Logging/Sink.h b/Framework/Logging/Sink.h index 6ee615176036f03c74a1deb50493bfdcd6cddd1b..f51b70294ebdf4ac63362448e5b851d7a52b62b3 100644 --- a/Framework/Logging/Sink.h +++ b/Framework/Logging/Sink.h @@ -1,37 +1,37 @@ #ifndef _include_Sink_h_ #define _include_Sink_h_ -namespace logger { +namespace corsika::logging { /** a sink for the logger must implement the two functions - operator<<(const std::string&) + operator<<(const std::string&) and Close() - + See example: NoSink */ - + namespace sink { - + /** - Definition of Sink for log output. - */ - template<typename TStream> + 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() {} + Sink(TStream& out) + : fOutput(out) {} + void operator<<(const std::string& msg) { fOutput << msg; } + void Close() {} + private: - TStream& fOutput; + TStream& fOutput; }; - + typedef Sink<std::ostream> SinkStream; - }// end namespace -} // end namespace - + } // namespace sink +} // namespace corsika::logging + #endif diff --git a/Framework/Logging/testLogging.cc b/Framework/Logging/testLogging.cc index 7d2ccb56977a70d1b953e70df4f8a7ca42e5fdd7..d28c46de441f66801367ef39f1f103ac3997b057 100644 --- a/Framework/Logging/testLogging.cc +++ b/Framework/Logging/testLogging.cc @@ -1,19 +1,13 @@ -#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one cpp file +#include <corsika/logging/Logger.h> + +#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one + // cpp file #include <catch2/catch.hpp> -#include <Logging/Logger.h> +TEST_CASE("Logging", "[Logging]") { + SECTION("sectionOne") {} -TEST_CASE( "Logging", "[Logging]" ) -{ - SECTION( "sectionOne" ) - { - } - - SECTION( "sectionTwo" ) - { - } + SECTION("sectionTwo") {} - SECTION( "sectionThree" ) - { - } + SECTION("sectionThree") {} } diff --git a/Framework/ParticleProperties/ParticleProperties.hpp b/Framework/ParticleProperties/ParticleProperties.hpp deleted file mode 100644 index f627966915e96af83dd282bc75cf16e65ca1b280..0000000000000000000000000000000000000000 --- a/Framework/ParticleProperties/ParticleProperties.hpp +++ /dev/null @@ -1,45 +0,0 @@ -#ifndef PARTICLEPROPERTIES_HPP -#define PARTICLEPROPERTIES_HPP -#include <array> -#include <cstdint> -#include <iostream> - -namespace ParticleProperties { - -typedef int16_t PDGCode; - -#include "generated_particle_properties.inc" - -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 getElectricChargeQN(InternalParticleCode const p) -{ - return electric_charge[static_cast<uint8_t const>(p)]; -} - -auto constexpr getElectricCharge(InternalParticleCode const p) -{ - return getElectricChargeQN(p) * (phys::units::e / 3.); -} - -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; -} - -} -#endif diff --git a/Framework/ParticleStack/CMakeLists.txt b/Framework/ParticleStack/CMakeLists.txt deleted file mode 100644 index 93b5280ec328b02eebe3ad5f10a1c73ed37e6c7c..0000000000000000000000000000000000000000 --- a/Framework/ParticleStack/CMakeLists.txt +++ /dev/null @@ -1,18 +0,0 @@ - -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) diff --git a/Framework/ParticleStack/StackIterator.h b/Framework/ParticleStack/StackIterator.h deleted file mode 100644 index 616ee53b89be6f786278e39ba38236c2c6ac9324..0000000000000000000000000000000000000000 --- a/Framework/ParticleStack/StackIterator.h +++ /dev/null @@ -1,70 +0,0 @@ -#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 diff --git a/Framework/ParticleStack/StackOne.h b/Framework/ParticleStack/StackOne.h deleted file mode 100644 index 0d8ab00680f334bd01bf7f072131ec933b65b33a..0000000000000000000000000000000000000000 --- a/Framework/ParticleStack/StackOne.h +++ /dev/null @@ -1,80 +0,0 @@ -#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 diff --git a/Framework/Particles/CMakeLists.txt b/Framework/Particles/CMakeLists.txt index b408a7690153aca39fd0f9453fa940d483068ef1..b7e2b58afca6ee46f33df5952bae4407e5243b88 100644 --- a/Framework/Particles/CMakeLists.txt +++ b/Framework/Particles/CMakeLists.txt @@ -1,41 +1,81 @@ -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) +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 + ParticleProperties.h + ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc # this one is auto-generated + ) -#set (PARTICLES_SOURCES ) -set (PARTICLE_HEADERS Particles.h ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc) +set ( + PARTICLE_NAMESPACE + corsika/particles + ) add_library (CORSIKAparticles INTERFACE) -#set_target_properties (CORSIKAparticles PROPERTIES VERSION ${PROJECT_VERSION}) -#set_target_properties (CORSIKAparticles PROPERTIES SOVERSION 1) +CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAparticles ${PARTICLE_NAMESPACE} ${PARTICLE_HEADERS}) -#set_target_properties (CORSIKAparticles PROPERTIES PUBLIC_HEADER "${PARTICLES_HEADERS}") - -# target dependencies on other libraries (also header only) -#target_link_libraries (CORSIKAparticles CORSIKAunits) +# .................................................... +# since GeneratedParticleProperties.inc is an automatically produced file in the build directory, +# create a symbolic link into the source tree, so that it can be found and edited more easily +# this is not needed for the build to succeed! ....... +add_custom_command ( + OUTPUT ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedParticleProperties.inc + COMMAND ${CMAKE_COMMAND} -E create_symlink ${PROJECT_BINARY_DIR}/include/corsika/GeneratedParticleProperties.inc ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedParticleProperties.inc + COMMENT "Generate link in source-dir: ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedParticleProperties.inc" + ) +add_custom_target (SourceDirLink DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedParticleProperties.inc) +add_dependencies (CORSIKAparticles SourceDirLink) +# ..................................................... -#target_include_directories (CORSIKAparticles PRIVATE ${EIGEN3_INCLUDE_DIR}) -target_include_directories (CORSIKAparticles INTERFACE ${EIGEN3_INCLUDE_DIR}) -target_include_directories (CORSIKAparticles INTERFACE $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Framework> - $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/Framework> - $<INSTALL_INTERFACE:include/Framework> - ) -#install (TARGETS CORSIKAparticles -# LIBRARY DESTINATION lib -# ARCHIVE DESTINATION lib -# PUBLIC_HEADER DESTINATION include/Particles) +target_include_directories ( + CORSIKAparticles + INTERFACE + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> + $<INSTALL_INTERFACE:include> + ) -install (FILES ${PARTICLE_HEADERS} DESTINATION include/Particles) +install ( + FILES + ${PARTICLE_HEADERS} + DESTINATION + include/${PARTICLE_NAMESPACE} + ) -# code 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) +# -------------------- +# 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 + ) diff --git a/Framework/Particles/ParticleClassNames.xml b/Framework/Particles/ParticleClassNames.xml index 8af3e990e20de3d17f64be65c3e5816220a78a6c..637f6b69597e0c8f6b8eceecfe6158f2e96dd460 100644 --- a/Framework/Particles/ParticleClassNames.xml +++ b/Framework/Particles/ParticleClassNames.xml @@ -9,5 +9,14 @@ <particle pdgID="-11" classname="Positron" /> <particle pdgID="22" classname="Gamma" /> + <particle id="130" classname="K0Long"> </particle> + <particle id="310" classname="K0Short"> </particle> + + <particle id="2112" classname="Neutron"> </particle> + <particle id="-2112" classname="AntiNeutron"> </particle> + + <particle id="2212" classname="Proton"> </particle> + <particle id="-2212" classname="AntiProton"> </particle> + </chapter> diff --git a/Framework/Particles/ParticleData.xml b/Framework/Particles/ParticleData.xml index e32cb637984b7407c3b409473ee1de17ad924a13..16e7878abe42147c0d2ec3a285465c975ecec93c 100644 --- a/Framework/Particles/ParticleData.xml +++ b/Framework/Particles/ParticleData.xml @@ -4716,14 +4716,15 @@ m0="9.89860" mWidth="0.01000" mMin="9.85500" mMax="9.89500"> <channel onMode="1" bRatio="1.0000000" meMode="91" products="21 21"/> </particle> - +--> + <particle id="13122" name="Lambda(1405)0" antiName="Lambda(1405)bar0" spinType="2" chargeType="0" colType="0" m0="1.40510" mWidth="0.05000" mMin="1.32000" mMax="1.80000"> <channel onMode="1" bRatio="0.3333000" products="3222 -211"/> <channel onMode="1" bRatio="0.3333000" products="3112 211"/> <channel onMode="1" bRatio="0.3334000" products="3212 111"/> </particle> - +<!-- <particle id="14122" name="Lambda_c(2593)+" antiName="Lambda_c(2593)-" spinType="2" chargeType="3" colType="0" m0="2.59225" mWidth="0.00260" mMin="2.57000" mMax="2.63000"> <channel onMode="1" bRatio="0.2400000" products="4222 -211"/> diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h new file mode 100644 index 0000000000000000000000000000000000000000..7df92b661d891cf252cbef52205e1818e71a834d --- /dev/null +++ b/Framework/Particles/ParticleProperties.h @@ -0,0 +1,60 @@ +/** + @file Particles.h + + Interface to particle properties + */ + +#ifndef _include_Particle_h_ +#define _include_Particle_h_ + +#include <array> +#include <cstdint> +#include <iostream> + +#include <corsika/units/PhysicalUnits.h> +#include <corsika/particles/GeneratedParticleProperties.inc> + +/** + * @namespace particle + * + * The properties of all elementary particles is stored here. The data + * is taken from the Pythia ParticleData.xml file. + * + */ + +namespace corsika::particles { + + /** + * @function GetMass + * + * return mass of particle + */ + auto constexpr GetMass(Code const p) { return masses[static_cast<uint8_t const>(p)]; } + + auto constexpr GetPDG(Code const p) { return pdg_codes[static_cast<uint8_t const>(p)]; } + + auto constexpr GetElectricChargeNumber(Code const p) { + return electric_charge[static_cast<uint8_t const>(p)] / 3; + } + + auto constexpr GetElectricCharge(Code const p) { + return GetElectricChargeNumber(p) * (corsika::units::constants::e); + } + + auto const GetName(Code const p) { return names[static_cast<uint8_t const>(p)]; } + + namespace io { + + std::ostream& operator<<(std::ostream& stream, Code const p) { + stream << GetName(p); + return stream; + } + + } // namespace io + +} // namespace corsika::particles + +// to inject the operator<< into the root namespace +using namespace corsika::particles::io; + +#endif diff --git a/Framework/Particles/Particles.h b/Framework/Particles/Particles.h deleted file mode 100644 index af521d7191791cb4f9f8df05d5270d0b2d233b2d..0000000000000000000000000000000000000000 --- a/Framework/Particles/Particles.h +++ /dev/null @@ -1,45 +0,0 @@ -#ifndef _include_Particle_h_ -#define _include_Particle_h_ - -#include <array> -#include <cstdint> -#include <iostream> - -#include <Particles/GeneratedParticleProperties.inc> - -namespace ParticleProperties { - - 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 GetElectricChargeQN(InternalParticleCode const p) - { - return electric_charge[static_cast<uint8_t const>(p)]; - } - - auto constexpr GetElectricCharge(InternalParticleCode const p) - { - return GetElectricChargeQN(p) * (phys::units::e / 3.); - } - - 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; - } - -} - -#endif diff --git a/Framework/Particles/pdxml_reader.py b/Framework/Particles/pdxml_reader.py index b1c0d3ee031733a90e694e391dad0c76bbc86e33..46958ea055785689515d4cb6aa6a127dbe125a7d 100755 --- a/Framework/Particles/pdxml_reader.py +++ b/Framework/Particles/pdxml_reader.py @@ -62,16 +62,16 @@ def class_names(filename): # # Automatically produce a string qualifying as C++ class name # - +# This function produces names of type "DELTA_PLUS_PLUS" +# def c_identifier(name): orig = name name = name.upper() - for c in "() ": + 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") @@ -91,8 +91,66 @@ def c_identifier(name): raise Exception("could not generate C identifier for '{:s}'".format(orig)) +############################################################## +# +# Automatically produce a string qualifying as C++ class name +# +# This function produces names of type "DeltaPlusPlus" +# +def c_identifier_cammel(name): + orig = name + name = name[0].upper() + name[1:].lower() # all lower case + for c in "() /": # replace funny characters + name = name.replace(c, "_") + + name = name.replace("bar", "Bar") + name = name.replace("*", "Star") + name = name.replace("'", "Prime") + name = name.replace("+", "Plus") + name = name.replace("-", "Minus") + + # move "Bar" to end of name + ibar = name.find('Bar') + if (ibar>0 and ibar<len(name)-3) : + name = name[:ibar] + name[ibar+3:] + str('Bar') + + # cleanup "_"s + while True: + tmp = name.replace("__", "_") + if tmp == name: + break + else: + name = tmp + name.strip("_") + + # remove all "_", if this does not by accident concatenate two number + istart = 0 + while True: + i = name.find('_', istart) + if (i<1 or i>len(name)-1): + break + istart = i + if (name[i-1].isdigit() and name[i+1].isdigit()): + # there is a number on both sides + break + name = name[:i] + name[i+1:] + # and last, for example: make NuE out of Nue + if (name[i-1].islower() and name[i].islower()) : + if (i<len(name)-1) : + name = name[:i] + name[i].upper() + name[i+1:] + else : + name = name[:i] + name[i].upper() + + # check if name is valid C++ identifier + pattern = re.compile(r'^[a-zA-Z_][a-zA-Z_0-9]*$') + if pattern.match(name): + return name + else: + raise Exception("could not generate C identifier for '{:s}' name='{:s}'".format(orig, name)) + + -######################################################### +# ######################################################## # # returns dict containing all data from pythia-xml input # @@ -105,7 +163,7 @@ def build_pythia_db(filename, classnames): if (pdg in classnames): c_id = classnames[pdg] else: - c_id = c_identifier(name) + c_id = c_identifier_cammel(name) # the cammel case names #~ print(name, c_id, sep='\t', file=sys.stderr) #~ enums += "{:s} = {:d}, ".format(c_id, corsika_id) @@ -128,10 +186,15 @@ def build_pythia_db(filename, classnames): # def gen_internal_enum(pythia_db): - string = "enum class InternalParticleCode : uint8_t {\n" - + string = "enum class Code : uint8_t {\n" + + string += " FirstParticle = 1, // if you want to loop over particles, you want to start with \"1\" \n" # identifier for eventual loops... + last_ngc_id = 0 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']) + last_ngc_id = pythia_db[k]['ngc_code'] + string += " {key:s} = {code:d},\n".format(key = k, code = last_ngc_id) + + string += " LastParticle = " + str(last_ngc_id+1) + ",\n" # identifier for eventual loops... string += "};" return string @@ -203,17 +266,27 @@ def gen_classes(pythia_db): break string += "\n"; - string += ("struct " + cname + " {\n" - " static constexpr InternalParticleCode GetType() { return Type; }\n" - " static auto constexpr GetMass() { return masses[TypeIndex]; }\n" - " static auto constexpr GetCharge() { return phys::units::e*electric_charge[TypeIndex]/3; }\n" - " static auto const GetName() { return names[TypeIndex]; }\n" - " static auto constexpr GetAntiParticle() { return AntiType; }\n" - " static auto constexpr Type = InternalParticleCode::") + cname + ";\n" + \ - " static auto constexpr AntiType = InternalParticleCode::" + antiP + ";\n" + \ - (" private:\n" - " static constexpr uint8_t TypeIndex = static_cast<uint8_t const>(Type);\n" - "};\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 Code GetCode() { return Type; }\n" + string += " static quantity<energy_d> GetMass() { return masses[TypeIndex]; }\n" + string += " static quantity<electric_charge_d> GetCharge() { return corsika::units::constants::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 Code GetAntiParticle() { return AntiType; }\n" + string += " static const Code Type = Code::" + cname + ";\n" + string += " static const Code AntiType = Code::" + antiP + ";\n" + string += " private:\n" + string += " static const uint8_t TypeIndex = static_cast<uint8_t const>(Type);\n" + string += "};\n" return string @@ -223,16 +296,22 @@ def gen_classes(pythia_db): # def inc_start(): - string = ("" - "#ifndef _include_GeneratedParticleDataTable_h_\n" - "#define _include_GeneratedParticleDataTable_h_\n\n" - "#include <array>\n" - "#include <cstdint>\n" - "#include <iostream>\n\n" - "using namespace phys::units;\n" - "using namespace phys::units::literals;\n\n" - "namespace ParticleProperties {\n\n" - "typedef int16_t PDGCode;\n\n") + string = "" + string += "#ifndef _include_GeneratedParticleDataTable_h_\n" + string += "#define _include_GeneratedParticleDataTable_h_\n\n" + string += "#include <corsika/units/PhysicalUnits.h>\n" + string += "#include <corsika/units/PhysicalConstants.h>\n" + string += "#include <array>\n" + string += "#include <cstdint>\n" +# string += "#include <iostream>\n\n" + string += "namespace corsika { \n\n" +# string += "using namespace literals; \n" + string += "namespace particles { \n\n" + string += "using corsika::units::energy_d;\n" + string += "using corsika::units::electric_charge_d;\n" + string += "using corsika::units::quantity;\n" + string += "using corsika::units::operator\"\"_GeV;\n" + string += "typedef int16_t PDGCode;\n\n" return string @@ -243,6 +322,7 @@ def inc_start(): def inc_end(): string = "" string += "\n}\n\n" + string += "\n}\n\n" string += "#endif\n" return string @@ -254,16 +334,18 @@ def inc_end(): if __name__ == "__main__": - if len(sys.argv) != 3: - print("usage: pdxml_reader.py Pythia8.xml ClassNames.xml", file=sys.stderr) + 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") + 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) diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc index 233706325cdb4b19c08423a58c07f3c3026e47fc..476264171bb8be56f715f0eee3ee9fe74f77116b 100644 --- a/Framework/Particles/testParticles.cc +++ b/Framework/Particles/testParticles.cc @@ -1,30 +1,24 @@ -#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 <Particles/Particles.h> +#include <corsika/particles/ParticleProperties.h> +#include <corsika/units/PhysicalUnits.h> -using namespace phys::units; -using namespace phys::units::literals; +#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one + // cpp file +#include <catch2/catch.hpp> -using namespace ParticleProperties; +using namespace corsika::units; +using namespace corsika::particles; -TEST_CASE( "Particles", "[Particles]" ) -{ - SECTION( "Types" ) - { - REQUIRE( Electron::GetType()==InternalParticleCode::Electron ); - } +TEST_CASE("Particles", "[Particles]") { - 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-" ); - } + SECTION("Types") { REQUIRE(Electron::GetCode() == Code::Electron); } + SECTION("Data") { + REQUIRE(Electron::GetMass() / 0.511_MeV == Approx(1)); + REQUIRE(Electron::GetMass() / GetMass(Code::Electron) == Approx(1)); + REQUIRE(Electron::GetCharge() / constants::e == Approx(-1)); + REQUIRE(Positron::GetCharge() / constants::e == Approx(+1)); + REQUIRE(GetElectricCharge(Positron::GetAntiParticle()) / constants::e == Approx(-1)); + REQUIRE(Electron::GetName() == "e-"); + } } diff --git a/Framework/ProcessSequence/CMakeLists.txt b/Framework/ProcessSequence/CMakeLists.txt index 966af83432bb3b662c9d83fae9be3111547b3382..dd84758450610555e4fc2fc0d62f701b5306df24 100644 --- a/Framework/ProcessSequence/CMakeLists.txt +++ b/Framework/ProcessSequence/CMakeLists.txt @@ -1,10 +1,51 @@ add_library (CORSIKAprocesssequence INTERFACE) -target_include_directories (CORSIKAprocesssequence INTERFACE $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Framework> - $<INSTALL_INTERFACE:include/Framework> - ) +# namespace of library -> location of header files +set ( + CORSIKAprocesssequence_NAMESPACE + corsika/process + ) -install (FILES ProcessSequence.h - DESTINATION include/ProcessSequence) +# header files of this library +set ( + CORSIKAprocesssequence_HEADERS + ProcessSequence.h + ) +CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAprocesssequence ${CORSIKAprocesssequence_NAMESPACE} ${CORSIKAprocesssequence_HEADERS}) + +# include directive for upstream code +target_include_directories ( + CORSIKAprocesssequence + INTERFACE + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> + $<INSTALL_INTERFACE:include/> + ) + +# install library +install ( + FILES ${CORSIKAprocesssequence_HEADERS} + DESTINATION include/${CORSIKAprocesssequence_NAMESPACE} + ) + +# ---------------- +# code unit testing +#add_executable ( +# testLogging +# testLogging.cc +# ) + +#target_link_libraries ( +# testLogging +# CORSIKAprocesssequence +# CORSIKAthirdparty # for catch2 +# ) + +#add_test ( +# NAME testLogging +# COMMAND testLogging +# ) + + + diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h index 967f6f347aacfd6789465417fa995ba7687f4402..621b7bb04cf38b97f34e6bc32d1f61ae1cf9710b 100644 --- a/Framework/ProcessSequence/ProcessSequence.h +++ b/Framework/ProcessSequence/ProcessSequence.h @@ -3,26 +3,21 @@ #include <iostream> #include <typeinfo> -using namespace std; -namespace processes { +namespace corsika::process { /** - /class Base - + /class BaseProcess + The structural base type of a process object in a ProcessSequence. Both, the ProcessSequence and all its elements - are of type Base<T> + are of type BaseProcess<T> */ - + template <typename derived> - struct Base - { - const derived& GetRef() const - { - return static_cast<const derived&>(*this); - } + struct BaseProcess { + const derived& GetRef() const { return static_cast<const derived&>(*this); } }; /** @@ -30,55 +25,57 @@ namespace processes { A compile time static list of processes. The compiler will generate a new type based on template logic containing all the - elements. + elements. - \comment Using CRTP pattern, https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern + \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> > - { + class ProcessSequence : public BaseProcess<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); } - + : 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() ); + inline const ProcessSequence<T1, T2> operator+(const BaseProcess<T1>& A, + const BaseProcess<T2>& B) { + return ProcessSequence<T1, T2>(A.GetRef(), B.GetRef()); } - - + /* template <typename T1> struct depth_lhs { static const int num = 0; }; - + // terminating condition @@ -89,29 +86,26 @@ namespace processes { 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) { @@ -122,7 +116,7 @@ namespace processes { } }; */ - + /* template<typename T1, typename T2> const Process& @@ -148,7 +142,6 @@ namespace processes { } */ -} // end namespace - -#endif +} // namespace corsika::process +#endif diff --git a/Framework/StackInterface/CMakeLists.txt b/Framework/StackInterface/CMakeLists.txt index f37313a26f26c2595eb43682195a85304f00822e..45a3ebb4de10a1cb739932805c6f4db15153fc8c 100644 --- a/Framework/StackInterface/CMakeLists.txt +++ b/Framework/StackInterface/CMakeLists.txt @@ -1,10 +1,31 @@ +set ( + CORSIKAstackinterface_HEADERS + Stack.h + StackIterator.h + ) -add_library (CORSIKAstackinterface INTERFACE) +set ( + CORSIKAstackinterface_NAMESPACE + corsika/stack + ) -target_include_directories (CORSIKAstackinterface INTERFACE $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Framework> - $<INSTALL_INTERFACE:include/Framework> - ) +add_library ( + CORSIKAstackinterface + INTERFACE + ) -install (FILES Stack.h StackIterator.h - DESTINATION include/StackInterface) +CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAstackinterface ${CORSIKAstackinterface_NAMESPACE} ${CORSIKAstackinterface_HEADERS}) +target_include_directories ( + CORSIKAstackinterface + INTERFACE + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> + $<INSTALL_INTERFACE:include> + ) + +install ( + FILES + ${CORSIKAstackinterface_HEADERS} + DESTINATION + include/${CORSIKAstackinterface_NAMESPACE} + ) diff --git a/Framework/StackInterface/Stack.h b/Framework/StackInterface/Stack.h index 9b9065d7e0b0933710cf3baf28361b6f0109a546..42187c1aeec22b455edc2194b311e05e01dd9d4e 100644 --- a/Framework/StackInterface/Stack.h +++ b/Framework/StackInterface/Stack.h @@ -1,54 +1,56 @@ #ifndef _include_Stack_h__ #define _include_Stack_h__ -#include <StackInterface/StackIterator.h> // to help application programmres +#include <corsika/stack/StackIterator.h> // include here, to help application programmres /** All classes around management of particles on a stack. */ -namespace stack { +namespace corsika::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> + + 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: + using DataImpl::IncrementSize; + + 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); } - + 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); } + 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); } + iterator NewParticle() { + IncrementSize(); + return iterator(*this, GetSize() - 1); + } /// delete last particle on stack by decrementing stack size void DeleteLast() { DecrementSize(); } }; -} // end namespace - +} // namespace corsika::stack + #endif - diff --git a/Framework/StackInterface/StackIterator.h b/Framework/StackInterface/StackIterator.h index b2e9a8b6c523b8330731e729e5a7321a130fe3c5..4a951f00af33ddb865c6c7b7723a248bc44331b5 100644 --- a/Framework/StackInterface/StackIterator.h +++ b/Framework/StackInterface/StackIterator.h @@ -1,17 +1,18 @@ #ifndef _include_StackIterator_h__ #define _include_StackIterator_h__ -#include <iostream> #include <iomanip> +#include <iostream> -namespace stack { +namespace corsika::stack { // forward decl. - template<class Stack, class Particle> class StackIteratorInfo; + 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 @@ -19,78 +20,97 @@ namespace stack { This allows to write code like \verbatim - for (auto& p : theStack) { p.SetEnergy(newEnergy); } + 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 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 - { + + template <typename Stack, typename Particle> + class StackIterator : public Particle { friend Stack; friend Particle; - friend StackIteratorInfo<Stack,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; } + 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); } + // 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 ... + + This is the class where custom ... Internal helper class for StackIterator. Document better... */ - - template<typename _Stack, typename Particle> + + template <typename _Stack, typename Particle> class StackIteratorInfo { - - friend Particle; + + 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(); } + 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 - +} // namespace corsika::stack + #endif diff --git a/Framework/Units/CMakeLists.txt b/Framework/Units/CMakeLists.txt index 7881afc4f9264391915ca73b5a71dd21b06f9054..7948859e264449cb29bda5c3b1f141a86dbd2d65 100644 --- a/Framework/Units/CMakeLists.txt +++ b/Framework/Units/CMakeLists.txt @@ -1,14 +1,23 @@ add_library (CORSIKAunits INTERFACE) +set (CORSIKAunits_NAMESPACE corsika/units) +set ( + CORSIKAunits_HEADERS + PhysicalUnits.h + PhysicalConstants.h + ) + +CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAunits ${CORSIKAunits_NAMESPACE} ${CORSIKAunits_HEADERS}) + target_include_directories (CORSIKAunits INTERFACE - $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/Framework> + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> $<BUILD_INTERFACE:${PROJECT_SOURCE_DIR}/ThirdParty> $<INSTALL_INTERFACE:include> ) -install (FILES PhysicalUnits.h DESTINATION include/Units) +install (FILES ${CORSIKAunits_HEADERS} DESTINATION include/${CORSIKAunits_NAMESPACE}) # code testing add_executable (testUnits testUnits.cc) diff --git a/Framework/Units/PhysicalConstants.h b/Framework/Units/PhysicalConstants.h new file mode 100644 index 0000000000000000000000000000000000000000..a666f15ef3bfc61d7abfb589a9763597db4f7d62 --- /dev/null +++ b/Framework/Units/PhysicalConstants.h @@ -0,0 +1,60 @@ +/** + * \file PhysicalConstants + * + * \brief Several physical constants. + * \author Michael S. Kenniston, Martin Moene + * \date 7 September 2013 + * \since 0.4 + * + * Copyright 2013 Universiteit Leiden. All rights reserved. + * + * Copyright (c) 2001 by Michael S. Kenniston. For the most + * recent version check www.xnet.com/~msk/quantity. Permission is granted + * to use this code without restriction so long as this copyright + * notice appears in all source files. + * + * This code is provided as-is, with no warrantee of correctness. + * + * Distributed under the Boost Software License, Version 1.0. (See accompanying + * file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + * + * + * + */ + +#ifndef INCLUDE_PHYSICAL_CONSTANTS_H +#define INCLUDE_PHYSICAL_CONSTANTS_H + +#include <phys/units/quantity.hpp> + +namespace corsika::units::constants { + + using namespace phys::units; + + // acceleration of free-fall, standard + constexpr phys::units::quantity<phys::units::acceleration_d> g_sub_n{ + phys::units::Rep(9.80665L) * phys::units::meter / + phys::units::square(phys::units::second)}; + + // Avogadro constant + constexpr quantity<dimensions<0, 0, 0, 0, 0, -1> > N_sub_A{Rep(6.02214199e+23L) / mole}; + // electronvolt + constexpr quantity<energy_d> eV{Rep(1.60217733e-19L) * joule}; + + // elementary charge + constexpr quantity<electric_charge_d> e{Rep(1.602176462e-19L) * coulomb}; + + // Planck constant + constexpr quantity<dimensions<2, 1, -1> > h{Rep(6.62606876e-34L) * joule * second}; + + // speed of light in a vacuum + constexpr quantity<speed_d> c{Rep(299792458L) * meter / second}; + + // unified atomic mass unit + constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram}; + + // etc. + +} // namespace corsika::units::constants + +#endif // PHYS_UNITS_PHYSICAL_CONSTANTS_HPP_INCLUDED diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h index 1aa55dec1a80c2f2098e117e3f22c028c4ae187c..9aeefb15124b182add95860a197b4bebc213d98b 100644 --- a/Framework/Units/PhysicalUnits.h +++ b/Framework/Units/PhysicalUnits.h @@ -1,33 +1,43 @@ #ifndef _include_PhysicalUnits_h_ #define _include_PhysicalUnits_h_ -#include <phys/units/quantity.hpp> +#include <corsika/units/PhysicalConstants.h> + #include <phys/units/io.hpp> -#include <phys/units/physical_constants.hpp> +#include <phys/units/quantity.hpp> /** - @file PhysicalUnits - - Define _XeV literals, alowing 10_GeV in the code. -*/ - -/*using namespace phys::units::io; -using namespace phys::units::literals;*/ + * @file PhysicalUnits + * + * Define _XeV literals, alowing 10_GeV in the code. + */ namespace phys { namespace units { namespace literals { - QUANTITY_DEFINE_SCALING_LITERALS(eV, energy_d, magnitude(eV) ) + QUANTITY_DEFINE_SCALING_LITERALS(eV, energy_d, + magnitude(corsika::units::constants::eV)) } - } -} + } // namespace units +} // namespace phys -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>; -using Energy = phys::units::quantity<phys::units::energy_d, double>; +namespace corsika::units { -#endif + using namespace phys::units; + using namespace phys::units::literals; + // namespace literals = phys::units::literals; + using LengthType = phys::units::quantity<phys::units::length_d, double>; + using TimeType = phys::units::quantity<phys::units::time_interval_d, double>; + using SpeedType = phys::units::quantity<phys::units::speed_d, double>; + using FrequencyType = phys::units::quantity<phys::units::frequency_d, double>; + using ElectricChargeType = + phys::units::quantity<phys::units::electric_charge_d, double>; + using EnergyType = phys::units::quantity<phys::units::energy_d, double>; + +} // end namespace corsika::units + +// we want to call the operator<< without namespace... I think +using namespace phys::units::io; + +#endif diff --git a/Framework/Units/testUnits.cc b/Framework/Units/testUnits.cc index 5cf6e93c87bfb9304dd48e9ac4282699553b1ad2..080aabe8c9a08b38bcb89f9055b4a9951705254c 100644 --- a/Framework/Units/testUnits.cc +++ b/Framework/Units/testUnits.cc @@ -1,73 +1,72 @@ -#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one cpp file +#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 <corsika/units/PhysicalUnits.h> #include <array> -using namespace phys::units; -using namespace phys::units::literals; +using namespace corsika::units; -TEST_CASE( "PhysicalUnits", "[Units]" ) { +TEST_CASE("PhysicalUnits", "[Units]") { SECTION("Consistency") { - REQUIRE( 1_m/1_m == Approx(1) ); - //REQUIRE_FALSE( 1_m/1_s == 1 ); // static assert + 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 ); + REQUIRE(E1 == 10_GeV); - quantity<phys::units::length_d> l1 = 10_nm; - - quantity<phys::units::length_d> arr0[5]; + LengthType l1 = 10_nm; + + LengthType 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 }; + LengthType arr1[2] = {{1_mm}, {2_cm}}; + + std::array<EnergyType, 4> arr2; // empty array + + std::array<EnergyType, 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) ); + 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) ); + 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 EnergyType 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 ); + REQUIRE(E3 == 180_GeV); } - } diff --git a/Main/shower.cc b/Main/shower.cc index a0f78b5a8bcea889e582c07b6e11f1600010b2c7..7bc09673a29ed82e49f1939a52b3f5c86f31adad 100644 --- a/Main/shower.cc +++ b/Main/shower.cc @@ -1,10 +1,3 @@ #include <utl/Stack.h> -int -main(int argc, char** argv) -{ - - - - return 0; -} +int main(int argc, char** argv) { return 0; } diff --git a/Processes/CMakeLists.txt b/Processes/CMakeLists.txt index f1e2ce31a5222add07914f66e1f8b9d5b48be415..9b81c6d9857687cb808691565a1b5c9b6b413cfd 100644 --- a/Processes/CMakeLists.txt +++ b/Processes/CMakeLists.txt @@ -1 +1,2 @@ add_subdirectory (NullModel) +add_subdirectory (Sibyll) diff --git a/Processes/NullModel/CMakeLists.txt b/Processes/NullModel/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..ff8b026d993b5ccf9f98c2352a571c9e4dd3e018 --- /dev/null +++ b/Processes/NullModel/CMakeLists.txt @@ -0,0 +1,61 @@ + +set ( + MODEL_SOURCES + NullModel.cc + ) + +set ( + MODEL_HEADERS + NullModel.h + ) + +set ( + MODEL_NAMESPACE + corsika/process/null_model + ) + +add_library (ProcessNullModel STATIC ${MODEL_SOURCES}) +CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessNullModel ${MODEL_NAMESPACE} ${MODEL_HEADERS}) + +set_target_properties ( + ProcessNullModel + PROPERTIES + VERSION ${PROJECT_VERSION} + SOVERSION 1 +# PUBLIC_HEADER "${MODEL_HEADERS}" + ) + +# target dependencies on other libraries (also the header onlys) +target_link_libraries ( + ProcessNullModel + CORSIKAunits + ) + +target_include_directories ( + ProcessNullModel + INTERFACE + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> + $<INSTALL_INTERFACE:include/include> + ) + +install ( + TARGETS ProcessNullModel + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib +# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE} + ) + + +# -------------------- +# code unit testing +add_executable (testNullModel testNullModel.cc) + +target_link_libraries ( + testNullModel + CORSIKAgeometry + CORSIKAunits + CORSIKAthirdparty # for catch2 + ) + +add_test (NAME testNullModel COMMAND testNullModel) + diff --git a/Processes/NullModel/NullModel.cc b/Processes/NullModel/NullModel.cc index 96972982496c08a681e56b1a4cd68d604d266801..4da71dc16066f65af4623db7a45338ac9c8fde0c 100644 --- a/Processes/NullModel/NullModel.cc +++ b/Processes/NullModel/NullModel.cc @@ -1,2 +1,13 @@ -#include <Processes/NullModel/NullModel.h> +#include <corsika/process/null_model/NullModel.h> +using namespace corsika::process::null_model; + +NullModel::NullModel() {} + +NullModel::~NullModel() {} + +void NullModel::init() {} + +void NullModel::run() {} + +double NullModel::GetStepLength() { return 0; } diff --git a/Processes/NullModel/NullModel.h b/Processes/NullModel/NullModel.h index a78e14d14eeeb8b9f7b1739e7d29860a4f506b05..98782765696e7d0a87d19611bba0d9fc1e9f4c7f 100644 --- a/Processes/NullModel/NullModel.h +++ b/Processes/NullModel/NullModel.h @@ -1,20 +1,23 @@ #ifndef _Physics_NullModel_NullModel_h_ #define _Physics_NullModel_NullModel_h_ -namespace processes { +namespace corsika::process { - class NullModel { + namespace null_model { - public: - NullModel(); - ~NullModel(); - - void init(); - void run(); - double GetStepLength(); - }; - -} + class NullModel { -#endif + public: + NullModel(); + ~NullModel(); + + void init(); + void run(); + double GetStepLength(); + }; + + } // namespace null_model +} // namespace corsika::process + +#endif diff --git a/Processes/NullModel/testNullModel.cc b/Processes/NullModel/testNullModel.cc new file mode 100644 index 0000000000000000000000000000000000000000..311c732cec4b3cbdb8bcb4d4886b37df0e940696 --- /dev/null +++ b/Processes/NullModel/testNullModel.cc @@ -0,0 +1,12 @@ +#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one + // cpp file +#include <catch2/catch.hpp> + +#include <corsika/units/PhysicalUnits.h> + +TEST_CASE("NullModel", "[processes]") { + + SECTION("bla") {} + + SECTION("blubb") {} +} diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..041fdb8c34e8947d9ecc10b19961b7a2a653519b --- /dev/null +++ b/Processes/Sibyll/CMakeLists.txt @@ -0,0 +1,61 @@ + +set ( + MODEL_SOURCES + ParticleConversion.cc + ) + +set ( + MODEL_HEADERS + ParticleConversion.h + ) + +set ( + MODEL_NAMESPACE + corsika/process/sibyll + ) + +add_library (ProcessSibyll STATIC ${MODEL_SOURCES}) +CORSIKA_COPY_HEADERS_TO_NAMESPACE (ProcessSibyll ${MODEL_NAMESPACE} ${MODEL_HEADERS}) + +set_target_properties ( + ProcessSibyll + PROPERTIES + VERSION ${PROJECT_VERSION} + SOVERSION 1 +# PUBLIC_HEADER "${MODEL_HEADERS}" + ) + +# target dependencies on other libraries (also the header onlys) +target_link_libraries ( + ProcessSibyll + CORSIKAunits + ) + +target_include_directories ( + ProcessSibyll + INTERFACE + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> + $<INSTALL_INTERFACE:include/include> + ) + +install ( + TARGETS ProcessSibyll + LIBRARY DESTINATION lib + ARCHIVE DESTINATION lib +# PUBLIC_HEADER DESTINATION include/${MODEL_NAMESPACE} + ) + + +# -------------------- +# code unit testing +add_executable (testSibyll testSibyll.cc) + +target_link_libraries ( + testSibyll + CORSIKAgeometry + CORSIKAunits + CORSIKAthirdparty # for catch2 + ) + +add_test (NAME testSibyll COMMAND testSibyll) + diff --git a/Processes/Sibyll/ParticleConversion.cc b/Processes/Sibyll/ParticleConversion.cc new file mode 100644 index 0000000000000000000000000000000000000000..6aca2fcf0a07a1ee89eeb859fbb6b40b37e96ec3 --- /dev/null +++ b/Processes/Sibyll/ParticleConversion.cc @@ -0,0 +1,9 @@ +#include <corsika/particles/ParticleProperties.h> +#include <corsika/process/sibyll/ParticleConversion.h> + +using namespace corsika::process::sibyll; + +// const std::map<sibyll::PID, ParticleProperties::InternalParticleCode> +// process::sibyll::Sibyll2Corsika = { +// {PID::E_MINUS, InternalParticleCode::Electron}, +//}; diff --git a/Processes/Sibyll/ParticleConversion.h b/Processes/Sibyll/ParticleConversion.h new file mode 100644 index 0000000000000000000000000000000000000000..f39c1297964d33846bedcb6160203f11f882c95d --- /dev/null +++ b/Processes/Sibyll/ParticleConversion.h @@ -0,0 +1,199 @@ +#ifndef _include_processes_sibyll_particles_h_ +#define _include_processes_sibyll_particles_h_ + +#include <corsika/particles/ParticleProperties.h> + +#include <map> + +namespace corsika::process { + + namespace sibyll { + + enum class PID : int { + E_MINUS = 3, + E_PLUS = 2, + NU_E = 15, + NU_E_BAR = 16, + MU_MINUS = 5, + MU_PLUS = 4, + NU_MU = 17, + NU_MU_BAR = 18, + TAU_MINUS = 91, + TAU_PLUS = 90, + NU_TAU = 92, + NU_TAU_BAR = 93, + GAMMA = 1, + PI_0 = 6, + RHO_0 = 27, + K_L_0 = 11, + PI_PLUS = 7, + PI_MINUS = 8, + RHO_PLUS = 25, + RHO_MINUS = 26, + ETA = 23, + OMEGA = 32, + K_S_0 = 12, + K_STAR_0 = 30, + K_STAR_BAR_0 = 31, + K_PLUS = 9, + K_MINUS = 10, + K_STAR_PLUS = 28, + K_STAR_MINUS = 29, + D_PLUS = 59, + D_MINUS = 60, + D_STAR_PLUS = 78, + D_STAR_MINUS = 79, + D_0 = 71, + D_BAR_0 = 72, + D_STAR_0 = 80, + D_STAR_BAR_0 = 81, + D_S_PLUS = 74, + D_S_MINUS = 75, + D_STAR_S_PLUS = 76, + D_STAR_S_MINUS = 77, + ETA_C = 73, + N_0 = 14, + N_BAR_0 = -14, + DELTA_0 = 42, + DELTA_BAR_0 = -42, + P_PLUS = 13, + P_BAR_MINUS = -13, + DELTA_PLUS = 41, + DELTA_BAR_MINUS = -41, + DELTA_PLUS_PLUS = 40, + DELTA_BAR_MINUS_MINUS = -40, + SIGMA_MINUS = 36, + SIGMA_BAR_PLUS = -36, + LAMBDA_0 = 39, + LAMBDA_BAR_0 = -39, + SIGMA_0 = 35, + SIGMA_BAR_0 = -35, + SIGMA_PLUS = 34, + SIGMA_BAR_MINUS = -34, + XI_MINUS = 38, + XI_BAR_PLUS = -38, + XI_0 = 37, + XI_BAR_0 = -37, + OMEGA_MINUS = 49, + OMEGA_BAR_PLUS = -49, + SIGMA_C_0 = 86, + SIGMA_C_BAR_0 = -86, + SIGMA_STAR_C_0 = 96, + SIGMA_STAR_C_BAR_0 = -96, + LAMBDA_C_PLUS = 89, + LAMBDA_C_BAR_MINUS = -89, + XI_C_0 = 88, + XI_C_BAR_0 = -88, + SIGMA_C_PLUS = 85, + SIGMA_C_BAR_MINUS = -85, + SIGMA_STAR_C_PLUS = 95, + SIGMA_STAR_C_BAR_MINUS = -95, + SIGMA_C_PLUS_PLUS = 84, + SIGMA_C_BAR_MINUS_MINUS = -84, + SIGMA_STAR_C_PLUS_PLUS = 94, + SIGMA_STAR_C_BAR_MINUS_MINUS = -94, + XI_C_PLUS = 87, + XI_C_BAR_MINUS = -87, + OMEGA_C_0 = 99, + OMEGA_C_BAR_0 = -99, + J_PSI = 83, + VOID = 0, + }; + + static const std::map<sibyll::PID, corsika::particles::Code> Sibyll2Corsika = { + {PID::E_MINUS, corsika::particles::Code::Electron}, + {PID::E_PLUS, corsika::particles::Code::Positron}, + {PID::NU_E, corsika::particles::Code::NuE}, + {PID::NU_E_BAR, corsika::particles::Code::NuEBar}, + {PID::MU_MINUS, corsika::particles::Code::MuMinus}, + {PID::MU_PLUS, corsika::particles::Code::MuPlus}, + {PID::NU_MU, corsika::particles::Code::NuMu}, + {PID::NU_MU_BAR, corsika::particles::Code::NuMuBar}, + {PID::TAU_MINUS, corsika::particles::Code::TauMinus}, + /* + TAU_PLUS = 90, + NU_TAU = 92, + NU_TAU_BAR = 93, + GAMMA = 1, + PI_0 = 6, + RHO_0 = 27, + K_L_0 = 11, + PI_PLUS = 7, + PI_MINUS = 8, + RHO_PLUS = 25, + RHO_MINUS = 26, + ETA = 23, + OMEGA = 32, + K_S_0 = 12, + K_STAR_0 = 30, + K_STAR_BAR_0 = 31, + K_PLUS = 9, + K_MINUS = 10, + K_STAR_PLUS = 28, + K_STAR_MINUS = 29, + D_PLUS = 59, + D_MINUS = 60, + D_STAR_PLUS = 78, + D_STAR_MINUS = 79, + D_0 = 71, + D_BAR_0 = 72, + D_STAR_0 = 80, + D_STAR_BAR_0 = 81, + D_S_PLUS = 74, + D_S_MINUS = 75, + D_STAR_S_PLUS = 76, + D_STAR_S_MINUS = 77, + ETA_C = 73, + N_0 = 14, + N_BAR_0 = -14, + DELTA_0 = 42, + DELTA_BAR_0 = -42, + P_PLUS = 13, + P_BAR_MINUS = -13, + DELTA_PLUS = 41, + DELTA_BAR_MINUS = -41, + DELTA_PLUS_PLUS = 40, + DELTA_BAR_MINUS_MINUS = -40, + SIGMA_MINUS = 36, + SIGMA_BAR_PLUS = -36, + LAMBDA_0 = 39, + LAMBDA_BAR_0 = -39, + SIGMA_0 = 35, + SIGMA_BAR_0 = -35, + SIGMA_PLUS = 34, + SIGMA_BAR_MINUS = -34, + XI_MINUS = 38, + XI_BAR_PLUS = -38, + XI_0 = 37, + XI_BAR_0 = -37, + OMEGA_MINUS = 49, + OMEGA_BAR_PLUS = -49, + SIGMA_C_0 = 86, + SIGMA_C_BAR_0 = -86, + SIGMA_STAR_C_0 = 96, + SIGMA_STAR_C_BAR_0 = -96, + LAMBDA_C_PLUS = 89, + LAMBDA_C_BAR_MINUS = -89, + XI_C_0 = 88, + XI_C_BAR_0 = -88, + SIGMA_C_PLUS = 85, + SIGMA_C_BAR_MINUS = -85, + SIGMA_STAR_C_PLUS = 95, + SIGMA_STAR_C_BAR_MINUS = -95, + SIGMA_C_PLUS_PLUS = 84, + SIGMA_C_BAR_MINUS_MINUS = -84, + SIGMA_STAR_C_PLUS_PLUS = 94, + SIGMA_STAR_C_BAR_MINUS_MINUS = -94, + XI_C_PLUS = 87, + XI_C_BAR_MINUS = -87, + OMEGA_C_0 = 99, + OMEGA_C_BAR_0 = -99, + J_PSI = 83, + VOID = 0,*/ + }; + + } // namespace sibyll + +} // namespace corsika::process + +#endif diff --git a/Processes/Sibyll/Particles.cc b/Processes/Sibyll/Particles.cc deleted file mode 100644 index fc449b0223badebad1249a48bf60af3b3e8b5fea..0000000000000000000000000000000000000000 --- a/Processes/Sibyll/Particles.cc +++ /dev/null @@ -1,9 +0,0 @@ -#include <Sibyll/Particle.h> - -using namespace processes; - -const std::map<sibyll::PID, ParticleProperties::InternalParticleCode> -processes::sibyll::Sibyll2Corsika = { - {sibyll::PID::E_MINUS, InternalParticleCode::Electron}, -}; - diff --git a/Processes/Sibyll/Particles.h b/Processes/Sibyll/Particles.h deleted file mode 100644 index aa6ddfd5073a7ab742704b87a4c158ee4ced5892..0000000000000000000000000000000000000000 --- a/Processes/Sibyll/Particles.h +++ /dev/null @@ -1,199 +0,0 @@ -#ifndef _include_processes_sibyll_particles_h_ -#define _include_processes_sibyll_particles_h_ - -#include <Particles/Particles.h> - -#include <map> - -namespace processes { - - namespace sibyll { - - enum class PID : int { - E_MINUS = 3, - E_PLUS = 2, - NU_E=15, - NU_E_BAR=16, - MU_MINUS=5, - MU_PLUS=4, - NU_MU=17, - NU_MU_BAR=18, - TAU_MINUS=91, - TAU_PLUS=90, - NU_TAU=92, - NU_TAU_BAR=93, - GAMMA=1, - /*etc etc etc - PI_0 6 - RHO_0 27 - K_L_0 11 - PI_PLUS 7 - PI_MINUS 8 - RHO_PLUS 25 - RHO_MINUS 26 - ETA 23 - OMEGA 32 - K_S_0 12 - K_STAR_0 30 - K_STAR_BAR_0 31 - K_PLUS 9 - K_MINUS 10 - K_STAR_PLUS 28 - K_STAR_MINUS 29 - D_PLUS 59 - D_MINUS 60 - D_STAR_PLUS 78 - D_STAR_MINUS 79 - D_0 71 - D_BAR_0 72 - D_STAR_0 80 - D_STAR_BAR_0 81 - D_S_PLUS 74 - D_S_MINUS 75 - D_STAR_S_PLUS 76 - D_STAR_S_MINUS 77 - ETA_C 73 - N_0 14 - N_BAR_0 -14 - DELTA_0 42 - DELTA_BAR_0 -42 - P_PLUS 13 - P_BAR_MINUS -13 - DELTA_PLUS 41 - DELTA_BAR_MINUS -41 - DELTA_PLUS_PLUS 40 - DELTA_BAR_MINUS_MINUS -40 - SIGMA_MINUS 36 - SIGMA_BAR_PLUS -36 - LAMBDA_0 39 - LAMBDA_BAR_0 -39 - SIGMA_0 35 - SIGMA_BAR_0 -35 - SIGMA_PLUS 34 - SIGMA_BAR_MINUS -34 - XI_MINUS 38 - XI_BAR_PLUS -38 - XI_0 37 - XI_BAR_0 -37 - OMEGA_MINUS 49 - OMEGA_BAR_PLUS -49 - SIGMA_C_0 86 - SIGMA_C_BAR_0 -86 - SIGMA_STAR_C_0 96 - SIGMA_STAR_C_BAR_0 -96 - LAMBDA_C_PLUS 89 - LAMBDA_C_BAR_MINUS -89 - XI_C_0 88 - XI_C_BAR_0 -88 - SIGMA_C_PLUS 85 - SIGMA_C_BAR_MINUS -85 - SIGMA_STAR_C_PLUS 95 - SIGMA_STAR_C_BAR_MINUS -95 - SIGMA_C_PLUS_PLUS 84 - SIGMA_C_BAR_MINUS_MINUS -84 - SIGMA_STAR_C_PLUS_PLUS 94 - SIGMA_STAR_C_BAR_MINUS_MINUS -94 - XI_C_PLUS 87 - XI_C_BAR_MINUS -87 - OMEGA_C_0 99 - OMEGA_C_BAR_0 -99 - J_PSI 83 - VOID 0 - */ - }; - - static const std::map<sibyll::PID, ParticleProperties::InternalParticleCode> Sibyll2Corsika; - - -/* -E_MINUS 3 -E_PLUS 2 -NU_E 15 -NU_E_BAR 16 -MU_MINUS 5 -MU_PLUS 4 -NU_MU 17 -NU_MU_BAR 18 -TAU_MINUS 91 -TAU_PLUS 90 -NU_TAU 92 -NU_TAU_BAR 93 -GAMMA 1 -PI_0 6 -RHO_0 27 -K_L_0 11 -PI_PLUS 7 -PI_MINUS 8 -RHO_PLUS 25 -RHO_MINUS 26 -ETA 23 -OMEGA 32 -K_S_0 12 -K_STAR_0 30 -K_STAR_BAR_0 31 -K_PLUS 9 -K_MINUS 10 -K_STAR_PLUS 28 -K_STAR_MINUS 29 -D_PLUS 59 -D_MINUS 60 -D_STAR_PLUS 78 -D_STAR_MINUS 79 -D_0 71 -D_BAR_0 72 -D_STAR_0 80 -D_STAR_BAR_0 81 -D_S_PLUS 74 -D_S_MINUS 75 -D_STAR_S_PLUS 76 -D_STAR_S_MINUS 77 -ETA_C 73 -N_0 14 -N_BAR_0 -14 -DELTA_0 42 -DELTA_BAR_0 -42 -P_PLUS 13 -P_BAR_MINUS -13 -DELTA_PLUS 41 -DELTA_BAR_MINUS -41 -DELTA_PLUS_PLUS 40 -DELTA_BAR_MINUS_MINUS -40 -SIGMA_MINUS 36 -SIGMA_BAR_PLUS -36 -LAMBDA_0 39 -LAMBDA_BAR_0 -39 -SIGMA_0 35 -SIGMA_BAR_0 -35 -SIGMA_PLUS 34 -SIGMA_BAR_MINUS -34 -XI_MINUS 38 -XI_BAR_PLUS -38 -XI_0 37 -XI_BAR_0 -37 -OMEGA_MINUS 49 -OMEGA_BAR_PLUS -49 -SIGMA_C_0 86 -SIGMA_C_BAR_0 -86 -SIGMA_STAR_C_0 96 -SIGMA_STAR_C_BAR_0 -96 -LAMBDA_C_PLUS 89 -LAMBDA_C_BAR_MINUS -89 -XI_C_0 88 -XI_C_BAR_0 -88 -SIGMA_C_PLUS 85 -SIGMA_C_BAR_MINUS -85 -SIGMA_STAR_C_PLUS 95 -SIGMA_STAR_C_BAR_MINUS -95 -SIGMA_C_PLUS_PLUS 84 -SIGMA_C_BAR_MINUS_MINUS -84 -SIGMA_STAR_C_PLUS_PLUS 94 -SIGMA_STAR_C_BAR_MINUS_MINUS -94 -XI_C_PLUS 87 -XI_C_BAR_MINUS -87 -OMEGA_C_0 99 -OMEGA_C_BAR_0 -99 -J_PSI 83 -VOID 0 -*/ - -#endif diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc new file mode 100644 index 0000000000000000000000000000000000000000..4c5bed4df1b42c351587569ed5a0b192f19f1f2d --- /dev/null +++ b/Processes/Sibyll/testSibyll.cc @@ -0,0 +1,27 @@ + +#include <corsika/particles/ParticleProperties.h> +#include <corsika/process/sibyll/ParticleConversion.h> +#include <corsika/units/PhysicalUnits.h> + +#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; + +TEST_CASE("Sibyll", "[processes]") { + + SECTION("ParticleConversion") { + REQUIRE(corsika::particles::Electron::GetCode() == + process::sibyll::Sibyll2Corsika.at(process::sibyll::PID::E_MINUS)); + } + + SECTION("Data") { + REQUIRE(corsika::particles::GetName(process::sibyll::Sibyll2Corsika.at( + process::sibyll::PID::E_PLUS)) == "e+"); + } + + SECTION("bla") {} + + SECTION("blubb") {} +} diff --git a/Stack/CMakeLists.txt b/Stack/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..4971f73f8bcad97257ea4fb90fc28f25679677d0 --- /dev/null +++ b/Stack/CMakeLists.txt @@ -0,0 +1,2 @@ + +add_subdirectory (SuperStupidStack) diff --git a/Stack/SuperStupidStack/CMakeLists.txt b/Stack/SuperStupidStack/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..6b00c993c4c712a4a2107ebd1539d21431da1752 --- /dev/null +++ b/Stack/SuperStupidStack/CMakeLists.txt @@ -0,0 +1,29 @@ + +set (SuperStupidStack_HEADERS SuperStupidStack.h) +set (SuperStupidStack_NAMESPACE corsika/stack/super_stupid) + +add_library (SuperStupidStack INTERFACE) + +CORSIKA_COPY_HEADERS_TO_NAMESPACE (SuperStupidStack ${SuperStupidStack_NAMESPACE} ${SuperStupidStack_HEADERS}) + +target_link_libraries ( + SuperStupidStack + INTERFACE + CORSIKAstackinterface + CORSIKAunits + CORSIKAparticles + ) + +target_include_directories ( + SuperStupidStack + INTERFACE + $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include> + $<INSTALL_INTERFACE:include> + ) + +install ( + FILES + ${SuperStupidStack_HEADERS} + DESTINATION + include/${SuperStupidStack_NAMESPACE} + ) diff --git a/Stack/SuperStupidStack/SuperStupidStack.h b/Stack/SuperStupidStack/SuperStupidStack.h new file mode 100644 index 0000000000000000000000000000000000000000..1f02e9087b82f18ef8674c149aad73786d0d3cda --- /dev/null +++ b/Stack/SuperStupidStack/SuperStupidStack.h @@ -0,0 +1,95 @@ +#ifndef _include_superstupidstack_h_ +#define _include_superstupidstack_h_ + +#include <string> +#include <vector> + +#include <corsika/particles/ParticleProperties.h> +#include <corsika/stack/Stack.h> +#include <corsika/units/PhysicalUnits.h> + +namespace corsika::stack { + + namespace super_stupid { + + using corsika::particles::Code; + using corsika::units::EnergyType; + using corsika::units::operator""_GeV; // literals; + + /** + * Example of a particle object on the stack. + */ + + template <typename _Stack> + class ParticleRead : public StackIteratorInfo<_Stack, ParticleRead<_Stack> > { + + using StackIteratorInfo<_Stack, ParticleRead>::GetIndex; + using StackIteratorInfo<_Stack, ParticleRead>::GetStack; + + public: + void SetId(const Code id) { GetStack().SetId(GetIndex(), id); } + void SetEnergy(const EnergyType& e) { GetStack().SetEnergy(GetIndex(), e); } + + Code GetId() const { return GetStack().GetId(GetIndex()); } + const EnergyType& GetEnergy() const { return GetStack().GetEnergy(GetIndex()); } + }; + + /** + * + * Memory implementation of the most simple (stupid) particle stack object. + */ + + class SuperStupidStackImpl { + + public: + void Clear() { + fDataE.clear(); + fDataId.clear(); + } + + int GetSize() const { return fDataId.size(); } + int GetCapacity() const { return fDataId.size(); } + + void SetId(const int i, const Code id) { fDataId[i] = id; } + void SetEnergy(const int i, const EnergyType& e) { fDataE[i] = e; } + + const Code GetId(const int i) const { return fDataId[i]; } + const EnergyType& GetEnergy(const int i) const { return fDataE[i]; } + + /** + * Function to copy particle at location i2 in stack to i1 + */ + void Copy(const int i1, const int i2) { + fDataE[i2] = fDataE[i1]; + fDataId[i2] = fDataId[i1]; + } + + protected: + void IncrementSize() { + fDataE.push_back(0_GeV); + fDataId.push_back(Code::unknown); + } + void DecrementSize() { + if (fDataE.size() > 0) { + fDataE.pop_back(); + fDataId.pop_back(); + } + } + + private: + /// the actual memory to store particle data + + std::vector<Code> fDataId; + std::vector<EnergyType> fDataE; + + }; // end class SuperStupidStackImpl + + typedef StackIterator<SuperStupidStackImpl, ParticleRead<SuperStupidStackImpl> > + Particle; + typedef Stack<SuperStupidStackImpl, Particle> SuperStupidStack; + + } // namespace super_stupid + +} // namespace corsika::stack + +#endif diff --git a/ThirdParty/ThirdParty.dox b/ThirdParty/ThirdParty.dox index 956b180f5c8936762d0083b75c9bc599314a7769..ec0e3de3826ccab8b9685df1b54c357e86957cfc 100644 --- a/ThirdParty/ThirdParty.dox +++ b/ThirdParty/ThirdParty.dox @@ -25,4 +25,8 @@ The catch2 unit testing library is from: https://github.com/catchorg/Catch2 Licence: BSL-1.0 (https://github.com/martinmoene/PhysUnits-CT-Cpp11/blob/master/LICENSE_1_0.txt) References: https://github.com/catchorg/Catch2 +@section eigen3 + +eigen3 .... + */ diff --git a/ThirdParty/eigen-eigen-b3f3d4950030.tar b/ThirdParty/eigen-eigen-b3f3d4950030.tar new file mode 100644 index 0000000000000000000000000000000000000000..4c050d1cc7e51abd98fc0c03b81b350c87b2b5f3 Binary files /dev/null and b/ThirdParty/eigen-eigen-b3f3d4950030.tar differ diff --git a/corsika.dox b/corsika.dox index cf0e6d393182618ba93c65a99eec5dd0781debe5..d48248f7fdb75dec91c53fff786f5d65d4b98abf 100644 --- a/corsika.dox +++ b/corsika.dox @@ -1,17 +1,30 @@ /** @mainpage CORSIKA air shower simulations framework -Documentation and reference guide for the CORSIKA8 software framework -for air shower simulations. CORSIKA is the most comprehensive -framework for air shower simulations. CORSIKA is characterized by -simulating particle cascades in astrophysical environments. The impact -of stochastic and continuous processes on the cascade development is -simulated. To boost computational efficiency different techniques are -provided, like thinning or cascade equations. +Documentation and reference guide for the CORSIKA8 (CORSIKA version 8) +software framework for air shower simulations. CORSIKA8 is developed +at <a +href="https://gitlab.ikp.kit.edu/AirShowerPhysics">https://gitlab.ikp.kit.edu</a>. If +you got the code from somewhere else, consider to switch to the +official development repository. If you want to report bugs, or want +to suggest features or future development, please submit an "issue" on +this gitlab server. + +CORSIKA is the most comprehensive framework for simulating particle +cascades in astrophysical environments, for example extensive air +showers. The impact of stochastic and continuous processes on the +cascade development is simulated. To boost computational efficiency +different techniques are provided, like thinning or cascade equations. The software makes extensive use of static design patterns and compiler optimization. Thus, the most fundamental configuration decision of the user must be performed at compile time. At runtime only spcific parameters can still be changed. +When you contribute to CORSIKA, follow the guidlines outlined here: + +<a +href="https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/wikis/home">Development +wiki text</a> + */ diff --git a/do-clang-format.sh b/do-clang-format.sh new file mode 100755 index 0000000000000000000000000000000000000000..e84335c520be53e1f011cf0466d8006b91085f83 --- /dev/null +++ b/do-clang-format.sh @@ -0,0 +1,3 @@ +clang-format -i -style=file `find . -iregex '^.*\.\(cc\|h\)$' -not -path './ThirdParty/*'` + +