diff --git a/.gitignore b/.gitignore index d1458d2..bff5b40 100644 --- a/.gitignore +++ b/.gitignore @@ -36,3 +36,6 @@ benchmark/* /profile/* /statprof/* /debug/* + +# Compiled packages +JetReconstructionCompiled diff --git a/compile/Project.toml b/compile/Project.toml new file mode 100644 index 0000000..1cd2875 --- /dev/null +++ b/compile/Project.toml @@ -0,0 +1,12 @@ +[deps] +ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" +JetReconstruction = "44e8cb2c-dfab-4825-9c70-d4808a591196" +PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d" + +[sources] +JetReconstruction = {path = ".."} + +[compat] +JetReconstruction = "0.4" +PackageCompiler = "2" +julia = "1.11" diff --git a/compile/README.md b/compile/README.md new file mode 100644 index 0000000..6ee62ea --- /dev/null +++ b/compile/README.md @@ -0,0 +1,87 @@ +# Compiling JetReconstruction.jl to a C-library + +Minimal C bindings for JetReconstruction.jl + +- [C-header](include/JetReconstruction.h) +- shared library compiled with [PackageCompiler.jl](https://github.com/JuliaLang/PackageCompiler.jl) + +## Building library + +To build the library, run the following command from the package root directory: + +```sh +julia --project=compile compile/build.jl +``` + +## Usage example + +### Example source file + +Example usage of C bindings in an application: + +```C +#include "JetReconstruction.h" +#include "julia_init.h" /*Should be automatically generated by PackageCompiler.jl and distributed together with the "JetReconstruction.h" header file*/ + +int main(int argc, char *argv[]) { + init_julia(argc, argv); /*initialization of julia runtime*/ + + /*Prepare array of pseudo jets*/ + size_t particles_len; + jetreconstruction_PseudoJet* particles; + /*Initialize with desired values*/ + + /*Call jet reconstruction*/ + jetreconstruction_JetAlgorithm algorithm = JETRECONSTRUCTION_JETALGORITHM_CA; + double R = 3.0; + jetreconstruction_RecoStrategy strategy = JETRECONSTRUCTION_RECOSTRATEGY_BEST; + + jetreconstruction_ClusterSequence cluster_seq; + int ret = jetreconstruction_jet_reconstruct(particles, len, algorithm, R, strategy, + &cluster_seq); + if (ret != JETRECONSTRUCTION_STATUSCODE_OK){ + /*An error occurred check the value or stderr for more information*/ + return 1; + } + + /*Use the cluster sequence in your application + then free memory allocations done by library*/ + jetreconstruction_ClusterSequence_free_members(&cluster_seq); + shutdown_julia(0); /*teardown of julia runtime*/ + return 0; +} + +``` + +### Example compilation + +To build an example application run the following command: + +```shell +cc -o jetreconstruction_test compile/test/jetreconstruction_test.c -IJetReconstructionCompiled/include -LJetReconstructionCompiled/lib -ljetreconstruction -ljulia +``` + +In case the compiled library resides in non-standard location, add its location to `LD_LIBRARY_PATH` when running example application: + +```shell +LD_LIBRARY_PATH=JetReconstructionCompiled/lib/:${LD_LIBRARY_PATH} ./jetreconstruction_test +``` + +### Compilation with CMake + +The JetReconstruction library comes with a CMake target `JetReconstruction::JetReconstruction`. Example usage in CMake file: + +```cmake +find_package(JetReconstruction REQUIRED) + +target_link_libraries(myTarget PUBLIC JetReconstruction::JetReconstruction) +``` + +## Limitations + +Currently it's not possible to create libraries for different platforms - no cross-compilation! + +The library is relocatable given the whole installation tree is moved, including libraries in the `lib/julia/` directory. + +It's advised to install the library in a separate directory to avoid possible conflicts. +The library must not be installed in the same directory as another Julia package compiled with `PackageCompiler.jl` as they would overwrite the package specific files in `share/julia`. diff --git a/compile/build.jl b/compile/build.jl new file mode 100644 index 0000000..717d753 --- /dev/null +++ b/compile/build.jl @@ -0,0 +1,65 @@ +using PackageCompiler +import ArgParse +import JetReconstruction + +function parse_args(args) + s = ArgParse.ArgParseSettings() + ArgParse.@add_arg_table s begin + "--source-dir" + help = "Directory containing the source files" + arg_type = String + default = splitdir(@__DIR__) |> first + + "--output-dir", "-o" + help = "Directory to save the compiled library" + arg_type = String + default = joinpath(splitdir(@__DIR__) |> first, "JetReconstructionCompiled") + end + + return ArgParse.parse_args(args, s) +end + +function configure_file(template_path::String, output_path::String, + replacements::Dict{String, String}) + template = read(template_path, String) + for (key, value) in replacements + template = replace(template, "@$(key)@" => value) + end + open(output_path, "w") do io + write(io, template) + end +end + +function (@main)(args) + parsed_args = parse_args(args) + source_dir = parsed_args["source-dir"] + output_dir = parsed_args["output-dir"] + + @info "Compiling package from $source_dir" + @info "Creating library in $output_dir" + PackageCompiler.create_library(source_dir, output_dir; + lib_name = "jetreconstruction", + header_files = [joinpath(@__DIR__, "include", + "JetReconstruction.h")], + precompile_execution_file = [joinpath(@__DIR__, + "precompile_execution.jl")], + incremental = false, + filter_stdlibs = true, + force = true) + + cmake_input = joinpath(@__DIR__, "cmake", "JetReconstruction") + cmake_output = joinpath(output_dir, "lib", "cmake", "JetReconstruction") + + version = pkgversion(JetReconstruction) + cmake_project_version = "$(version.major).$(version.minor).$(version.patch)" + + @info "Copying CMake file to $cmake_output" + mkpath(cmake_output) + cp_file(input_dir, basename, output_dir) = cp(joinpath(input_dir, basename), + joinpath(output_dir, basename)) + cp_file(cmake_input, "JetReconstructionConfig.cmake", cmake_output) + cp_file(cmake_input, "JetReconstructionTargets.cmake", cmake_output) + configure_file(joinpath(cmake_input, "JetReconstructionConfigVersion.cmake.in"), + joinpath(cmake_output, "JetReconstructionConfigVersion.cmake"), + Dict("PROJECT_VERSION" => cmake_project_version)) +end diff --git a/compile/cmake/JetReconstruction/JetReconstructionConfig.cmake b/compile/cmake/JetReconstruction/JetReconstructionConfig.cmake new file mode 100644 index 0000000..db499e6 --- /dev/null +++ b/compile/cmake/JetReconstruction/JetReconstructionConfig.cmake @@ -0,0 +1,39 @@ +# Config file for the JetReconstruction.jl package +# Manualy adjusted from standard cmake generated config file +get_filename_component(PACKAGE_PREFIX_DIR "${CMAKE_CURRENT_LIST_DIR}/../../../" ABSOLUTE) + +macro(set_and_check _var _file) + set(${_var} "${_file}") + if(NOT EXISTS "${_file}") + message(FATAL_ERROR "File or directory ${_file} referenced by variable ${_var} does not exist !") + endif() +endmacro() + +macro(check_required_components _NAME) + foreach(comp ${${_NAME}_FIND_COMPONENTS}) + if(NOT ${_NAME}_${comp}_FOUND) + if(${_NAME}_FIND_REQUIRED_${comp}) + set(${_NAME}_FOUND FALSE) + endif() + endif() + endforeach() +endmacro() + +#################################################################################### + +# - Create relocatable paths to headers. +# NOTE: Do not strictly need paths as all usage requirements are encoded in +# the imported targets created later. +set_and_check(JetReconstruction_INCLUDE_DIR "${PACKAGE_PREFIX_DIR}/include") + +# - Create path to installed read-only data files (e.g. yaml files) +set_and_check(JetReconstruction_DATA_DIR "${PACKAGE_PREFIX_DIR}/share/julia") + +# - Include the targets file to create the imported targets that a client can +# link to (libraries) or execute (programs) +include("${CMAKE_CURRENT_LIST_DIR}/JetReconstructionTargets.cmake") + +# print the default "Found:" message and check library location +include(FindPackageHandleStandardArgs) +get_property(TEST_JETRECONSTRUCTION_LIBRARY TARGET JetReconstruction::JetReconstruction PROPERTY LOCATION) +find_package_handle_standard_args(JetReconstruction DEFAULT_MSG CMAKE_CURRENT_LIST_FILE TEST_JETRECONSTRUCTION_LIBRARY) diff --git a/compile/cmake/JetReconstruction/JetReconstructionConfigVersion.cmake.in b/compile/cmake/JetReconstruction/JetReconstructionConfigVersion.cmake.in new file mode 100644 index 0000000..c696669 --- /dev/null +++ b/compile/cmake/JetReconstruction/JetReconstructionConfigVersion.cmake.in @@ -0,0 +1,65 @@ +# This is a basic version file for the Config-mode of find_package(). +# It is used by write_basic_package_version_file() as input file for configure_file() +# to create a version-file which can be installed along a config.cmake file. +# +# The created file sets PACKAGE_VERSION_EXACT if the current version string and +# the requested version string are exactly the same and it sets +# PACKAGE_VERSION_COMPATIBLE if the current version is >= requested version, +# but only if the requested major version is the same as the current one. +# The variable CVF_VERSION must be set before calling configure_file(). + + +set(PACKAGE_VERSION "@PROJECT_VERSION@") + +if(PACKAGE_VERSION VERSION_LESS PACKAGE_FIND_VERSION) + set(PACKAGE_VERSION_COMPATIBLE FALSE) +else() + + if(PACKAGE_VERSION MATCHES "^([0-9]+)\\.") + set(CVF_VERSION_MAJOR "${CMAKE_MATCH_1}") + if(NOT CVF_VERSION_MAJOR VERSION_EQUAL 0) + string(REGEX REPLACE "^0+" "" CVF_VERSION_MAJOR "${CVF_VERSION_MAJOR}") + endif() + else() + set(CVF_VERSION_MAJOR ${PACKAGE_VERSION}) + endif() + + if(PACKAGE_FIND_VERSION_RANGE) + # both endpoints of the range must have the expected major version + math (EXPR CVF_VERSION_MAJOR_NEXT "${CVF_VERSION_MAJOR} + 1") + if (NOT PACKAGE_FIND_VERSION_MIN_MAJOR STREQUAL CVF_VERSION_MAJOR + OR ((PACKAGE_FIND_VERSION_RANGE_MAX STREQUAL "INCLUDE" AND NOT PACKAGE_FIND_VERSION_MAX_MAJOR STREQUAL CVF_VERSION_MAJOR) + OR (PACKAGE_FIND_VERSION_RANGE_MAX STREQUAL "EXCLUDE" AND NOT PACKAGE_FIND_VERSION_MAX VERSION_LESS_EQUAL CVF_VERSION_MAJOR_NEXT))) + set(PACKAGE_VERSION_COMPATIBLE FALSE) + elseif(PACKAGE_FIND_VERSION_MIN_MAJOR STREQUAL CVF_VERSION_MAJOR + AND ((PACKAGE_FIND_VERSION_RANGE_MAX STREQUAL "INCLUDE" AND PACKAGE_VERSION VERSION_LESS_EQUAL PACKAGE_FIND_VERSION_MAX) + OR (PACKAGE_FIND_VERSION_RANGE_MAX STREQUAL "EXCLUDE" AND PACKAGE_VERSION VERSION_LESS PACKAGE_FIND_VERSION_MAX))) + set(PACKAGE_VERSION_COMPATIBLE TRUE) + else() + set(PACKAGE_VERSION_COMPATIBLE FALSE) + endif() + else() + if(PACKAGE_FIND_VERSION_MAJOR STREQUAL CVF_VERSION_MAJOR) + set(PACKAGE_VERSION_COMPATIBLE TRUE) + else() + set(PACKAGE_VERSION_COMPATIBLE FALSE) + endif() + + if(PACKAGE_FIND_VERSION STREQUAL PACKAGE_VERSION) + set(PACKAGE_VERSION_EXACT TRUE) + endif() + endif() +endif() + + +# if the installed or the using project don't have CMAKE_SIZEOF_VOID_P set, ignore it: +if("${CMAKE_SIZEOF_VOID_P}" STREQUAL "" OR "8" STREQUAL "") + return() +endif() + +# check that the installed version has the same 32/64bit-ness as the one which is currently searching: +if(NOT CMAKE_SIZEOF_VOID_P STREQUAL "8") + math(EXPR installedBits "8 * 8") + set(PACKAGE_VERSION "${PACKAGE_VERSION} (${installedBits}bit)") + set(PACKAGE_VERSION_UNSUITABLE TRUE) +endif() diff --git a/compile/cmake/JetReconstruction/JetReconstructionTargets.cmake b/compile/cmake/JetReconstruction/JetReconstructionTargets.cmake new file mode 100644 index 0000000..7d548f2 --- /dev/null +++ b/compile/cmake/JetReconstruction/JetReconstructionTargets.cmake @@ -0,0 +1,96 @@ +# Generated by CMake + +if("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}" LESS 2.8) + message(FATAL_ERROR "CMake >= 2.8.0 required") +endif() +if(CMAKE_VERSION VERSION_LESS "3.0.0") + message(FATAL_ERROR "CMake >= 3.0.0 required") +endif() +cmake_policy(PUSH) +cmake_policy(VERSION 3.0.0...3.28) + +# Commands may need to know the format version. +set(CMAKE_IMPORT_FILE_VERSION 1) + +# Protect against multiple inclusion, which would fail when already imported targets are added once more. +set(_cmake_targets_defined "") +set(_cmake_targets_not_defined "") +set(_cmake_expected_targets "") +foreach(_cmake_expected_target IN ITEMS JetReconstruction::JetReconstruction) + list(APPEND _cmake_expected_targets "${_cmake_expected_target}") + if(TARGET "${_cmake_expected_target}") + list(APPEND _cmake_targets_defined "${_cmake_expected_target}") + else() + list(APPEND _cmake_targets_not_defined "${_cmake_expected_target}") + endif() +endforeach() +unset(_cmake_expected_target) +if(_cmake_targets_defined STREQUAL _cmake_expected_targets) + unset(_cmake_targets_defined) + unset(_cmake_targets_not_defined) + unset(_cmake_expected_targets) + unset(CMAKE_IMPORT_FILE_VERSION) + cmake_policy(POP) + return() +endif() +if(NOT _cmake_targets_defined STREQUAL "") + string(REPLACE ";" ", " _cmake_targets_defined_text "${_cmake_targets_defined}") + string(REPLACE ";" ", " _cmake_targets_not_defined_text "${_cmake_targets_not_defined}") + message(FATAL_ERROR "Some (but not all) targets in this export set were already defined.\nTargets Defined: ${_cmake_targets_defined_text}\nTargets not yet defined: ${_cmake_targets_not_defined_text}\n") +endif() +unset(_cmake_targets_defined) +unset(_cmake_targets_not_defined) +unset(_cmake_expected_targets) + + +# Compute the installation prefix relative to this file. +get_filename_component(_IMPORT_PREFIX "${CMAKE_CURRENT_LIST_FILE}" PATH) +get_filename_component(_IMPORT_PREFIX "${_IMPORT_PREFIX}" PATH) +get_filename_component(_IMPORT_PREFIX "${_IMPORT_PREFIX}" PATH) +get_filename_component(_IMPORT_PREFIX "${_IMPORT_PREFIX}" PATH) +if(_IMPORT_PREFIX STREQUAL "/") + set(_IMPORT_PREFIX "") +endif() + +# Create imported target JetReconstruction::JetReconstruction +add_library(JetReconstruction::JetReconstruction SHARED IMPORTED) +set_target_properties(JetReconstruction::JetReconstruction PROPERTIES + INTERFACE_INCLUDE_DIRECTORIES "${_IMPORT_PREFIX}/include" + IMPORTED_LOCATION "${_IMPORT_PREFIX}/lib/libjetreconstruction.so" + IMPORTED_SONAME "libjetreconstruction.so" +) + +# Cleanup temporary variables. +set(_IMPORT_PREFIX) + +# Loop over all imported files and verify that they actually exist +foreach(_cmake_target IN LISTS _cmake_import_check_targets) + if(CMAKE_VERSION VERSION_LESS "3.28" + OR NOT DEFINED _cmake_import_check_xcframework_for_${_cmake_target} + OR NOT IS_DIRECTORY "${_cmake_import_check_xcframework_for_${_cmake_target}}") + foreach(_cmake_file IN LISTS "_cmake_import_check_files_for_${_cmake_target}") + if(NOT EXISTS "${_cmake_file}") + message(FATAL_ERROR "The imported target \"${_cmake_target}\" references the file + \"${_cmake_file}\" +but this file does not exist. Possible reasons include: +* The file was deleted, renamed, or moved to another location. +* An install or uninstall procedure did not complete successfully. +* The installation package was faulty and contained + \"${CMAKE_CURRENT_LIST_FILE}\" +but not all the files it references. +") + endif() + endforeach() + endif() + unset(_cmake_file) + unset("_cmake_import_check_files_for_${_cmake_target}") +endforeach() +unset(_cmake_target) +unset(_cmake_import_check_targets) + +# This file does not depend on other imported targets which have +# been exported from the same project but in a separate export set. + +# Commands beyond this point should not need to know the version. +set(CMAKE_IMPORT_FILE_VERSION) +cmake_policy(POP) diff --git a/compile/downstream/CMakeLists.txt b/compile/downstream/CMakeLists.txt new file mode 100644 index 0000000..32ea65e --- /dev/null +++ b/compile/downstream/CMakeLists.txt @@ -0,0 +1,15 @@ +cmake_minimum_required(VERSION 3.12) + +project( + JetReconstruction_downstream + LANGUAGES C + VERSION 0.1.0) + +find_package(JetReconstruction 0.4.3 REQUIRED) + +add_executable(jetreconstruction_test jetreconstruction_test.c) +target_link_libraries(jetreconstruction_test + PUBLIC JetReconstruction::JetReconstruction) + +enable_testing() +add_test(NAME JetReconstructionTest COMMAND jetreconstruction_test) diff --git a/compile/downstream/README.md b/compile/downstream/README.md new file mode 100644 index 0000000..9e29950 --- /dev/null +++ b/compile/downstream/README.md @@ -0,0 +1,3 @@ +# Downstream project example + +Example downstream project using C-bindings of JetReconstruction.jl diff --git a/compile/downstream/jetreconstruction_test.c b/compile/downstream/jetreconstruction_test.c new file mode 100644 index 0000000..60417a1 --- /dev/null +++ b/compile/downstream/jetreconstruction_test.c @@ -0,0 +1,87 @@ +#include "JetReconstruction.h" +#include "julia_init.h" +#include +#include +#include +#include + +void printPseudoJet(const jetreconstruction_PseudoJet *jet) { + assert(jet != NULL); + printf("PseudoJet(%f %f %f %f %ld %f %f %f %f)\n", jet->px, jet->py, jet->pz, + jet->E, jet->_cluster_hist_index, jet->_pt2, jet->_inv_pt2, jet->_rap, + jet->_phi); +} + +void printHistoryElement(const jetreconstruction_HistoryElement *history) { + assert(history != NULL); + printf("HistoryElement(%ld %ld %ld %ld %lf %lf)\n", history->parent1, + history->parent2, history->child, history->jetp_index, history->dij, + history->max_dij_so_far); +} + +void printClusterSequence(const jetreconstruction_ClusterSequence *sequence) { + printf("Cluster Sequence Information:\n" + "Algorithm: %d\n" + "Power: %f\n" + "R parameter: %f\n" + "Strategy: %d\n" + "Initial number of jets: %ld\n" + "Total event energy (Qtot): %f\n", + sequence->algorithm, sequence->power, sequence->R, sequence->strategy, + sequence->n_initial_jets, sequence->Qtot); + printf("Number of jets: %zu\n", sequence->jets_length); + if (sequence->jets != NULL) { + for (size_t i = 0; i < sequence->jets_length; i++) { + printPseudoJet(sequence->jets + i); + } + } + printf("History length: %zu\n", sequence->history_length); + if (sequence->history != NULL) { + for (size_t i = 0; i < sequence->history_length; i++) { + printHistoryElement(sequence->history + i); + } + } +} + +void printJetsResult(const jetreconstruction_JetsResult *results) { + assert(results != NULL); + for (size_t i = 0; i < results->length; ++i) { + printPseudoJet(results->data + i); + } +} + +int main(int argc, char *argv[]) { + clock_t start_time = clock(); + int sc = 0; + init_julia(argc, argv); + size_t len = 2; + jetreconstruction_PseudoJet particles[2]; + sc = jetreconstruction_PseudoJet_init(&particles[0], 0.0, 1.0, 2.0, 3.0); + assert(sc == JETRECONSTRUCTION_STATUSCODE_OK); + sc = jetreconstruction_PseudoJet_init(&particles[1], 1.0, 2.0, 3.0, 4.0); + assert(sc == JETRECONSTRUCTION_STATUSCODE_OK); + + jetreconstruction_JetAlgorithm algorithm = JETRECONSTRUCTION_JETALGORITHM_CA; + double R = 3.0; + jetreconstruction_RecoStrategy strategy = JETRECONSTRUCTION_RECOSTRATEGY_BEST; + + jetreconstruction_ClusterSequence cluster_seq; + sc = jetreconstruction_jet_reconstruct(particles, len, algorithm, R, strategy, + &cluster_seq); + assert(sc == JETRECONSTRUCTION_STATUSCODE_OK); + + printClusterSequence(&cluster_seq); + jetreconstruction_JetsResult result; + sc = jetreconstruction_exclusive_jets_njets(&cluster_seq, 2, &result); + assert(sc == JETRECONSTRUCTION_STATUSCODE_OK); + printJetsResult(&result); + + jetreconstruction_JetsResult_free_members(&result); + jetreconstruction_ClusterSequence_free_members(&cluster_seq); + shutdown_julia(0); + + clock_t end_time = clock(); + double time_spent = (double)(end_time - start_time) / CLOCKS_PER_SEC; + printf("Execution time: %f seconds\n", time_spent); + return 0; +} diff --git a/compile/include/JetReconstruction.h b/compile/include/JetReconstruction.h new file mode 100644 index 0000000..5c0c963 --- /dev/null +++ b/compile/include/JetReconstruction.h @@ -0,0 +1,344 @@ +/* + * SPDX-FileCopyrightText: 2022-2024 JetReconstruction.jl authors, CERN + * SPDX-License-Identifier: MIT + */ + +/** + * @file JetReconstruction.h + * @brief Header file for the C-bindings of JetReconstruction.jl. + * + * This file contains the definitions of data structures and functions used in + * JetReconstruction.jl. The library provides functionality for efficient jet + * reconstruction. + */ + +#ifndef JET_RECONSTRUCTION_JL_H_ +#define JET_RECONSTRUCTION_JL_H_ + +#include + +/*Special states in a history.*/ + +/**Special history state: Invalid child for this jet, meaning it did not + * recombine further */ +#define JETRECONSTRUCTION_INVALID -3 +/**Special history state: Original cluster so it has no parent */ +#define JETRECONSTRUCTION_NONEXISTENTPARENT -2 +/**Special history state: Cluster recombined with beam" */ +#define JETRECONSTRUCTION_BEAMJET -1 + +/** + * @enum jetreconstruction_StatusCode + * @brief Status codes for Jet Reconstruction operations. + * + * Common status codes for Jet Reconstruction operations. Most of the standard + * Julia + * https://docs.julialang.org/en/v1/manual/control-flow/#Exception-Handling + * exceptions are mapped to corresponding codes. + */ +typedef enum { + JETRECONSTRUCTION_STATUSCODE_OK = 0, /**< The operation succeeded. */ + JETRECONSTRUCTION_STATUSCODE_GENERICEXCEPTION = + 1, /**< An unspecified error, not covered by other status codes occurred. + */ + JETRECONSTRUCTION_STATUSCODE_ARGUMENTERROR = 2, + JETRECONSTRUCTION_STATUSCODE_BOUNDSERROR = 3, + JETRECONSTRUCTION_STATUSCODE_COMPOSITEEXCEPTION = 4, + JETRECONSTRUCTION_STATUSCODE_DIMENSIONMISMATCH = 5, + JETRECONSTRUCTION_STATUSCODE_DIVIDEERROR = 6, + JETRECONSTRUCTION_STATUSCODE_DOMAINERROR = 7, + JETRECONSTRUCTION_STATUSCODE_EOFERROR = 8, + JETRECONSTRUCTION_STATUSCODE_ERROREXCEPTION = 9, + JETRECONSTRUCTION_STATUSCODE_INEXACTERROR = 10, + JETRECONSTRUCTION_STATUSCODE_INITERROR = 11, + JETRECONSTRUCTION_STATUSCODE_INTERRUPTEXCEPTION = 12, + JETRECONSTRUCTION_STATUSCODE_INVALIDSTATEEXCEPTION = 13, + JETRECONSTRUCTION_STATUSCODE_KEYERROR = 14, + JETRECONSTRUCTION_STATUSCODE_LOADERROR = 15, + JETRECONSTRUCTION_STATUSCODE_OUTOFMEMORYERROR = 16, + JETRECONSTRUCTION_STATUSCODE_READONLYMEMORYERROR = 17, + /*JETRECONSTRUCTION_STATUSCODE_REMOTEEXCEPTION = 18, distributed only */ + JETRECONSTRUCTION_STATUSCODE_METHODERROR = 19, + JETRECONSTRUCTION_STATUSCODE_OVERFLOWERROR = 20, + /*JETRECONSTRUCTION_STATUSCODE_PARSEERROR = 21, meta only*/ + JETRECONSTRUCTION_STATUSCODE_SYSTEMERROR = 22, + JETRECONSTRUCTION_STATUSCODE_TYPEERROR = 23, + JETRECONSTRUCTION_STATUSCODE_UNDEFREFERROR = 24, + JETRECONSTRUCTION_STATUSCODE_UNDEFVARERROR = 25, + JETRECONSTRUCTION_STATUSCODE_STRINGINDEXERROR = 26 +} jetreconstruction_StatusCode; + +/** + * @enum jetreconstruction_JetAlgorithm + * @brief Enumeration representing different jet algorithms used in + * the JetReconstruction. + */ +typedef enum { + JETRECONSTRUCTION_JETALGORITHM_ANTIKT = 0, /**< The Anti-Kt algorithm. */ + JETRECONSTRUCTION_JETALGORITHM_CA = 1, /**< The Cambridge/Aachen algorithm. */ + JETRECONSTRUCTION_JETALGORITHM_KT = 2, /**< The Inclusive-Kt algorithm. */ + JETRECONSTRUCTION_JETALGORITHM_GENKT = + 3, /**< The Generalised Kt algorithm (with arbitrary power). */ + JETRECONSTRUCTION_JETALGORITHM_EEKT = + 4, /**< The Generalised e+e- kt algorithm. */ + JETRECONSTRUCTION_JETALGORITHM_DURHAM = + 5 /**< The e+e- kt algorithm, aka Durham. */ +} jetreconstruction_JetAlgorithm; + +/** + * @enum jetreconstruction_RecoStrategy + * @brief Enumeration representing the different strategies for jet + * reconstruction. + */ +typedef enum { + JETRECONSTRUCTION_RECOSTRATEGY_BEST = 0, /**< The best strategy. */ + JETRECONSTRUCTION_RECOSTRATEGY_N2PLAIN = 1, /**< The plain N2 strategy. */ + JETRECONSTRUCTION_RECOSTRATEGY_N2TILTED = 2 /**< The tiled N2 strategy. */ +} jetreconstruction_RecoStrategy; + +/** + * @struct jetreconstruction_PseudoJet + * @brief A struct representing a pseudojet, a four-momentum object used in + * jet reconstruction algorithms. Additional information for the link back into + * the history of the clustering is stored in the _cluster_hist_index field. + * There is caching of the more expensive calculations for rapidity and + * azimuthal angle. + */ +typedef struct { + double px; /**< The x-component of the momentum. */ + double py; /**< The y-component of the momentum. */ + double pz; /**< The z-component of the momentum. */ + double E; /**< The energy component of the momentum. */ + long _cluster_hist_index; /**< The index of the cluster history. */ + double _pt2; /**< The squared transverse momentum. */ + double _inv_pt2; /**< The inverse squared transverse momentum. */ + double _rap; /**< The rapidity. */ + double _phi; /**< The azimuthal angle. */ +} jetreconstruction_PseudoJet; + +/** + * @brief Initializes a jetreconstruction_PseudoJet object with given momentum. + * + * @param[in] ptr Pointer to the jetreconstruction_PseudoJet object to be + * initialized. + * @param[in] px The x-component of the momentum. + * @param[in] py The y-component of the momentum. + * @param[in] pz The z-component of the momentum. + * @param[in] E The energy component of the momentum. + * @return An integer status code indicating the success or failure. See common + * status codes jetreconstruction_StatusCode. + */ +jetreconstruction_StatusCode +jetreconstruction_PseudoJet_init(jetreconstruction_PseudoJet *ptr, double px, + double py, double pz, double E); + +/** + * @struct jetreconstruction_HistoryElement + * @brief A struct holding a record of jet mergers and finalisations. + */ +typedef struct { + long parent1; /**< Index in history where first parent of this jet was + created (@ref JETRECONSTRUCTION_NONEXISTENTPARENT if this jet + is an original particle) */ + long parent2; /**< Index in history where second parent of this jet was + created (@ref JETRECONSTRUCTION_NONEXISTENTPARENT if this jet + is an original particle); @ref JETRECONSTRUCTION_BEAMJET if + this history entry just labels the fact that the jet has + recombined with the beam */ + long child; /**< Index in history where the current jet is recombined with + another jet to form its child. It is Invalid + if this jet does not further recombine. */ + long jetp_index; /**< Index in the jets vector where we will find the + PseudoJet object corresponding to this jet (i.e. the + jet created at this entry of the history). NB: if this + element of the history corresponds to a beam + recombination, then @p jetp_index = @ref + JETRECONSTRUCTION_INVALID. + */ + double dij; /**< The distance corresponding to the recombination at this + stage of the clustering. */ + double max_dij_so_far; /**< The largest recombination distance seen so far + in the clustering history */ +} jetreconstruction_HistoryElement; + +/** + * @struct jetreconstruction_ClusterSequence + * @brief A struct holding the full history of a jet clustering sequence, + * including the final jets. + */ +typedef struct { + jetreconstruction_JetAlgorithm + algorithm; /**< The algorithm used for clustering. */ + double power; /**< The power value used for the clustering algorithm */ + double R; /**< The R parameter used for the clustering algorithm. */ + jetreconstruction_RecoStrategy + strategy; /**< The strategy used for clustering. */ + jetreconstruction_PseudoJet + *jets; /**< A pointer to an array of jetreconstruction_PseudoJet + containing the actual jets in the cluster sequence. */ + size_t jets_length; /**< The length of @ref jets array. */ + long n_initial_jets; /**< The initial number of particles used for exclusive + jets. */ + jetreconstruction_HistoryElement + *history; /**< A pointer to an array of jetreconstruction_HistoryElement + containing the branching history of the cluster sequence. + Each stage in the history indicates where to look in the jets + vector to get the physical PseudoJet. */ + size_t history_length; /**< The length of @ref history array. */ + double Qtot; /**< The total energy of the event. */ +} jetreconstruction_ClusterSequence; + +/** @private */ +void jetreconstruction_ClusterSequence_free_members_( + jetreconstruction_ClusterSequence *ptr); + +/** + * @brief Frees the members of a jetreconstruction_ClusterSequence structure. + * + * This function releases any resources held by the members of the given + * jetreconstruction_ClusterSequence structure. + * + * @param[in,out] ptr Pointer to the jetreconstruction_ClusterSequence + * structure. + */ +static inline void jetreconstruction_ClusterSequence_free_members( + jetreconstruction_ClusterSequence *ptr) { + jetreconstruction_ClusterSequence_free_members_(ptr); + ptr->jets = NULL; + ptr->jets_length = 0; + ptr->history = NULL; + ptr->history_length = 0; +} + +/** + * @brief Reconstructs jets from a collection of particles using a specified + * algorithm and strategy. + * + * This function reconstructs jets from a collection of particles using a + * specified algorithm and strategy. + * + * @note The memory allocated for members of @p result must be freed using + * the @ref jetreconstruction_ClusterSequence_free_members function after the + * @p result is no longer needed to prevent memory leaks. + * + * @param[in] particles Pointer to an array of pseudojet objects used for jet + * reconstruction. + * @param[in] particles_length The length of @p particles array. + * @param[in] algorithm The algorithm to use for jet reconstruction. + * @param[in] R The jet radius parameter. + * @param[in] strategy The jet reconstruction strategy to use. + * @param[out] result A pointer to which a cluster sequence containing the + * reconstructed jets and the merging history will be stored. + * @return An integer status code indicating the success or failure. See common + * status codes jetreconstruction_StatusCode. + */ +jetreconstruction_StatusCode jetreconstruction_jet_reconstruct( + const jetreconstruction_PseudoJet *particles, size_t particles_length, + jetreconstruction_JetAlgorithm algorithm, double R, + jetreconstruction_RecoStrategy strategy, + jetreconstruction_ClusterSequence *result); + +/** + * @struct jetreconstruction_JetsResult + * @brief A structure to hold the inclusive or exclusive jets. + */ +typedef struct { + jetreconstruction_PseudoJet + *data; /**< Pointer to an array of jetreconstruction_PseudoJet. */ + size_t length; /**< The length of the @ref data array. */ +} jetreconstruction_JetsResult; + +/** @private */ +void jetreconstruction_JetsResult_free_members_( + jetreconstruction_JetsResult *ptr); + +/** + * @brief Frees the members of a jetreconstruction_JetsResult structure. + * + * This function releases any resources held by the members of the given + * jetreconstruction_JetsResult structure. + * + * @param[in,out] ptr A pointer to the jetreconstruction_JetsResult structure. + */ +static inline void +jetreconstruction_JetsResult_free_members(jetreconstruction_JetsResult *ptr) { + jetreconstruction_JetsResult_free_members_(ptr); + ptr->data = NULL; + ptr->length = 0; +} + +/** + * @brief Return all exclusive jets of a jetreconstruction_ClusterSequence with + * a cut on the maximum distance parameter. + * + * This function computes the exclusive jets from a given + * jetreconstruction_ClusterSequence with a cut on the maximum distance + * parameter. + * + * @note The memory allocated for members of @p result must be freed using + * the @ref jetreconstruction_JetsResult_free_members function after the + * @p result is no longer needed to prevent memory leaks. + * + * @param[in] clustersequence A pointer to the jetreconstruction_ClusterSequence + * object containing the clustering history and jets. + * @param[in] dcut The distance parameter used to define the exclusive jets. + * @param[out] result A pointer to the jetreconstruction_JetsResult object where + * the resulting jets will be stored. + * @return An integer status code indicating the success or failure. See common + * status codes jetreconstruction_StatusCode. + */ +jetreconstruction_StatusCode jetreconstruction_exclusive_jets_dcut( + const jetreconstruction_ClusterSequence *clustersequence, double dcut, + jetreconstruction_JetsResult *result); + +/** + * @brief Return all exclusive jets of a jetreconstruction_ClusterSequence with + * a specific number of jets. + * + * This function computes a specific number of exclusive jets from a given + * jetreconstruction_ClusterSequence. + * + * @note The memory allocated for members of @p result must be freed using + * the @ref jetreconstruction_JetsResult_free_members function after the + * @p result is no longer needed to prevent memory leaks. + * + * @param[in] clustersequence A pointer to the jetreconstruction_ClusterSequence + * object containing the clustering history and jets. + * @param[in] njets The number of exclusive jets to be calculated. + * @param[out] result A pointer to the jetreconstruction_JetsResult object where + * the resulting jets will be stored. + * @return An integer status code indicating the success or failure. See common + * status codes jetreconstruction_StatusCode. + */ +jetreconstruction_StatusCode jetreconstruction_exclusive_jets_njets( + const jetreconstruction_ClusterSequence *clustersequence, size_t njets, + jetreconstruction_JetsResult *result); + +/** + * @brief Return all inclusive jets of a jetreconstruction_ClusterSequence with + * pt > @p ptmin. + * + * This function computes the inclusive jets from a given + * jetreconstruction_ClusterSequence. It iterates over the clustering history + * and checks the transverse momentum of each parent jet. If the transverse + * momentum is greater than or equal to @p ptmin, the jet is added to results + * structure. + * + * @note The memory allocated for members of @p result must be freed using + * the @ref jetreconstruction_JetsResult_free_members function after the + * @p result is no longer needed to prevent memory leaks. + * + * @param[in] clustersequence A pointer to the jetreconstruction_ClusterSequence + * object containing the clustering history and jets. + * @param[in] ptmin The minimum transverse momentum (pt) threshold for the + * inclusive jets. + * @param[out] result A pointer to the jetreconstruction_JetsResult object where + * the resulting jets will be stored. + * @return An integer status code indicating the success or failure. See common + * status codes jetreconstruction_StatusCode. + */ +jetreconstruction_StatusCode jetreconstruction_inclusive_jets( + const jetreconstruction_ClusterSequence *clustersequence, double ptmin, + jetreconstruction_JetsResult *result); + +#endif /* JET_RECONSTRUCTION_JL_H_ */ diff --git a/compile/precompile_execution.jl b/compile/precompile_execution.jl new file mode 100644 index 0000000..1a388ac --- /dev/null +++ b/compile/precompile_execution.jl @@ -0,0 +1,43 @@ +# Dummy code to call and cache compilation of all C-bindings +# The numeric values doesn't make sense and are used only to select correct dispatch + +using JetReconstruction +using JetReconstruction.C_JetReconstruction: jetreconstruction_PseudoJet_init, + C_ClusterSequence, + jetreconstruction_ClusterSequence_free_members_, + C_JetsResult, + jetreconstruction_JetsResult_free_members_, + jetreconstruction_jet_reconstruct, + jetreconstruction_inclusive_jets, + jetreconstruction_exclusive_jets_njets, + jetreconstruction_exclusive_jets_dcut + +pseudoJets_len = Csize_t(2) +pseudoJets_ptr = Ptr{PseudoJet}(Libc.malloc(pseudoJets_len * sizeof(PseudoJet))) + +jetreconstruction_PseudoJet_init(pseudoJets_ptr, 0.0, 1.0, 2.0, 3.0) +jetreconstruction_PseudoJet_init(pseudoJets_ptr + sizeof(PseudoJet), 1.0, 2.0, 3.0, 4.0) + +strategy = RecoStrategy.Best +algorithm = JetAlgorithm.CA +R = 2.0 + +clustersequence_ptr = Ptr{C_ClusterSequence{PseudoJet}}(Libc.malloc(sizeof(C_ClusterSequence{PseudoJet}))) +jetreconstruction_jet_reconstruct(pseudoJets_ptr, pseudoJets_len, algorithm, R, + RecoStrategy.Best, + clustersequence_ptr) + +results_ptr = Ptr{C_JetsResult{PseudoJet}}(Libc.malloc(sizeof(C_JetsResult{PseudoJet}))) +jetreconstruction_exclusive_jets_njets(clustersequence_ptr, Csize_t(2), results_ptr) +jetreconstruction_JetsResult_free_members_(results_ptr) + +jetreconstruction_exclusive_jets_dcut(clustersequence_ptr, 1.0, results_ptr) +jetreconstruction_JetsResult_free_members_(results_ptr) + +jetreconstruction_inclusive_jets(clustersequence_ptr, 0.0, results_ptr) +jetreconstruction_JetsResult_free_members_(results_ptr) + +jetreconstruction_ClusterSequence_free_members_(clustersequence_ptr) +Libc.free(pseudoJets_ptr) +Libc.free(clustersequence_ptr) +Libc.free(results_ptr) diff --git a/src/C_JetReconstruction/C_JetReconstruction.jl b/src/C_JetReconstruction/C_JetReconstruction.jl new file mode 100644 index 0000000..9aa06a1 --- /dev/null +++ b/src/C_JetReconstruction/C_JetReconstruction.jl @@ -0,0 +1,470 @@ +""" +Minimal C bindings for JetReconstruction.jl +""" +module C_JetReconstruction + +using ..JetReconstruction: JetAlgorithm, RecoStrategy, PseudoJet, ClusterSequence, + HistoryElement, jet_reconstruct, + exclusive_jets, inclusive_jets +using EnumX + +""" + @enumx StatusCode + +An enumeration of common status codes used in the JetReconstruction C-bindings. + +## Fields +- `OK` - The operation succeeded. +- `GenericException`- An unspecified error, not covered by other status codes occurred. +- The rest of the status codes have corresponding standard Julia exception. +""" +@enumx StatusCode OK=0 GenericException=1 ArgumentError=2 BoundsError=3 CompositeException=4 DimensionMismatch=5 DivideError=6 DomainError=7 EOFError=8 ErrorException=9 InexactError=10 InitError=11 InterruptException=12 InvalidStateException=13 KeyError=14 LoadError=15 OutOfMemoryError=16 ReadOnlyMemoryError=17 RemoteException=18 MethodError=19 OverflowError=20 ParseError=21 SystemError=22 TypeError=23 UndefRefError=24 UndefVarError=25 StringIndexError=26 + +function handle_exception(exception)::Cint + @error exception + return Cint(exception_to_enum(exception)) +end + +""" + exception_to_enum(::Any)::Cint + +Helper function matching Julia exception with C-style status code. +# Returns +- A numerical representation of the corresponding status code from the `StatusCode.T`. +""" +exception_to_enum(::Any) = Cint(StatusCode.GenericException) +exception_to_enum(::ArgumentError) = Cint(StatusCode.ArgumentError) +exception_to_enum(::BoundsError) = Cint(StatusCode.BoundsError) +exception_to_enum(::CompositeException) = Cint(StatusCode.CompositeException) +exception_to_enum(::DimensionMismatch) = Cint(StatusCode.DimensionMismatch) +exception_to_enum(::DivideError) = Cint(StatusCode.DivideError) +exception_to_enum(::DomainError) = Cint(StatusCode.DomainError) +exception_to_enum(::EOFError) = Cint(StatusCode.EOFError) +exception_to_enum(::ErrorException) = Cint(StatusCode.ErrorException) +exception_to_enum(::InexactError) = Cint(StatusCode.InexactError) +exception_to_enum(::InitError) = Cint(StatusCode.InitError) +exception_to_enum(::InterruptException) = Cint(StatusCode.InterruptException) +exception_to_enum(::InvalidStateException) = Cint(StatusCode.InvalidStateException) +exception_to_enum(::KeyError) = Cint(StatusCode.KeyError) +exception_to_enum(::LoadError) = Cint(StatusCode.LoadError) +exception_to_enum(::OutOfMemoryError) = Cint(StatusCode.OutOfMemoryError) +exception_to_enum(::ReadOnlyMemoryError) = Cint(StatusCode.ReadOnlyMemoryError) +# exception_to_enum(::Distributed.RemoteException) = Cint(StatusCode.RemoteException) +exception_to_enum(::MethodError) = Cint(StatusCode.MethodError) +exception_to_enum(::OverflowError) = Cint(StatusCode.OverflowError) +# exception_to_enum(::Meta.ParseError) = Cint(StatusCode.ParseError) +exception_to_enum(::SystemError) = Cint(StatusCode.SystemError) +exception_to_enum(::TypeError) = Cint(StatusCode.TypeError) +exception_to_enum(::UndefRefError) = Cint(StatusCode.UndefRefError) +exception_to_enum(::UndefVarError) = Cint(StatusCode.UndefVarError) +exception_to_enum(::StringIndexError) = Cint(StatusCode.StringIndexError) + +""" + unsafe_wrap_c_array(ptr::Ptr{T}, array_length::Csize_t) where {T} + +Wraps a C array into a Julia `Vector` for both bits and non-bits types. + +# Arguments +- `ptr::Ptr{T}`: A pointer to the C array. +- `array_length::Csize_t`: The length of the C array. + +# Returns +- A Julia `Vector{T}` containing the elements of the C array. + +# Safety +This function use 'unsafe' methods and has undefined behaviour +if pointer isn't valid or length isn't correct. +""" +function unsafe_wrap_c_array(ptr::Ptr{T}, array_length::Csize_t) where {T} + if isbitstype(T) + return unsafe_wrap(Vector{T}, ptr, array_length) + end + + vec = Vector{T}(undef, array_length) + for i in eachindex(vec) + @inbounds vec[i] = unsafe_load(ptr, i) + end + return vec +end + +""" + make_c_array(v::Vector{T}) where {T} + +Helper function for converting a Julia vector to a C-style array. +A C-style array is dynamically allocated and contents of input vector copied to it. + +# Arguments +- `v::Vector{T}`: A Julia vector of type `T`. + +# Returns +- `ptr::Ptr{T}`: A pointer to the allocated C-style array. +- `len::Int`: The length of the vector. + +# Throws +- `OutOfMemoryError`: If memory allocation fails. + +# Notes +- The caller is responsible for freeing the allocated memory using `Libc.free(ptr)`. +""" +function make_c_array(v::Vector{T}) where {T} + len = length(v) + ptr = Ptr{T}(Libc.malloc(len * sizeof(T))) + if ptr == C_NULL + throw(OutOfMemoryError("Libc.malloc failed to allocate memory")) + end + for i in 1:len + @inbounds unsafe_store!(ptr, v[i], i) + end + return ptr, len +end + +""" + jetreconstruction_PseudoJet_init(ptr::Ptr{PseudoJet}, px::Cdouble, py::Cdouble, pz::Cdouble, E::Cdouble) -> Cint + +C-binding for `PseudoJet` initialization. + +# Arguments +- `ptr::Ptr{PseudoJet}`: A pointer to the memory location where the `PseudoJet` object will be stored. +- `px::Cdouble`: The x-component of the momentum. +- `py::Cdouble`: The y-component of the momentum. +- `pz::Cdouble`: The z-component of the momentum. +- `E::Cdouble`: The energy of the jet. + +# Returns +- `Cint`: An integer status code indicating the success or failure. + +""" +Base.@ccallable function jetreconstruction_PseudoJet_init(ptr::Ptr{PseudoJet}, px::Cdouble, + py::Cdouble, pz::Cdouble, + E::Cdouble)::Cint + try + pseudojet = PseudoJet(px, py, pz, E) + unsafe_store!(ptr, pseudojet) + catch e + return handle_exception(e) + end + return Cint(StatusCode.OK) +end + +""" + struct C_ClusterSequence{T} + +A C-compatible struct corresponding to `ClusterSequence` +""" +struct C_ClusterSequence{T} + algorithm::JetAlgorithm.Algorithm + power::Cdouble + R::Cdouble + strategy::RecoStrategy.Strategy + jets::Ptr{T} + jets_length::Csize_t + n_initial_jets::Clong + history::Ptr{HistoryElement} + history_length::Csize_t + Qtot::Cdouble +end + +""" + free_members(ptr::Ptr{C_ClusterSequence{T}}) where {T} + +Internal function for freeing dynamically allocated `C_ClusterSequence` members memory. + +# Arguments +- `ptr::Ptr{C_ClusterSequence{T}}`: A pointer to a `C_ClusterSequence` structure. +""" +function free_members(ptr::Ptr{C_ClusterSequence{T}}) where {T} + if ptr != C_NULL + clusterseq = unsafe_load(ptr) + Libc.free(clusterseq.jets) + Libc.free(clusterseq.history) + # Struct is immutable so pointers can't be assigned NULL and lengths updated to zero (without making a copy) + end +end + +""" + jetreconstruction_ClusterSequence_free_members_(ptr::Ptr{C_ClusterSequence{PseudoJet}}) -> Cvoid + +C-binding for freeing the members of a `C_ClusterSequence` object pointed to by `ptr`. + +# Arguments +- `ptr::Ptr{C_ClusterSequence{PseudoJet}}`: A pointer to the `C_ClusterSequence` object whose members are to be freed. + +# Returns +- `Cvoid`: This function does not return a value. +""" +Base.@ccallable function jetreconstruction_ClusterSequence_free_members_(ptr::Ptr{C_ClusterSequence{PseudoJet}})::Cvoid + free_members(ptr) + return nothing +end + +""" + ClusterSequence(c::C_ClusterSequence{T}) where {T} + +Convert a `C_ClusterSequence` object to a `ClusterSequence` object. + + +# Arguments +- `c::C_ClusterSequence{T}`: The `C_ClusterSequence` object to be converted. + +# Returns +- `ClusterSequence{T}`: The converted `ClusterSequence` object. + +# Notes +- The array members of output are wrapping the array members of input. +- The input object must remain valid while the output object is used. +""" +function ClusterSequence{T}(c::C_ClusterSequence{T}) where {T} + jets_v = unsafe_wrap_c_array(c.jets, c.jets_length) + history_v = unsafe_wrap_c_array(c.history, c.history_length) + return ClusterSequence{T}(c.algorithm, c.power, c.R, c.strategy, jets_v, + c.n_initial_jets, + history_v, c.Qtot) +end + +""" + C_ClusterSequence(clustersequence::ClusterSequence{T}) where {T} + +Convert a `ClusterSequence` object to a `C_ClusterSequence` object. + +# Arguments +- `clustersequence::ClusterSequence{T}`: The `ClusterSequence` object to be converted. + +# Returns +- `C_ClusterSequence{T}`: The converted `C_ClusterSequence` object. + +# Notes +- The array members of input are deep-copied to output object, it's safe to use output even if input doesn't exist anymore. +- The memory allocated for output data-object members should be freed with the `free_members` function or equivalent. +""" +function C_ClusterSequence{T}(clustersequence::ClusterSequence{T}) where {T} + jets_ptr, jets_length = make_c_array(clustersequence.jets) + history_ptr, history_length = make_c_array(clustersequence.history) + return C_ClusterSequence{T}(clustersequence.algorithm, clustersequence.power, + clustersequence.R, clustersequence.strategy, jets_ptr, + jets_length, clustersequence.n_initial_jets, history_ptr, + history_length, clustersequence.Qtot) +end + +""" + c_jet_reconstruct(particles::Ptr{T}, + particles_length::Csize_t, + algorithm::JetAlgorithm.Algorithm, + R::Cdouble, + strategy::RecoStrategy.Strategy, + result::Ptr{C_ClusterSequence{U}})::Cint where {T, U} + +Internal helper functions for calling `jet_reconstruct` with C-compatible data-structers. + +# Arguments +- `particles::Ptr{T}`: Pointer to an array of pseudojet objects used for jet reconstruction. +- `particles_length::Csize_t`: The length of `particles`` array. +- `algorithm::JetAlgorithm.Algorithm`: The algorithm to use for jet reconstruction. +- `R::Cdouble`: The jet radius parameter.. +- `strategy::RecoStrategy.Strategy`: The jet reconstruction strategy to use. +- `result::Ptr{C_ClusterSequence{U}}`: A pointer to which a cluster sequence containing the reconstructed jets and the merging history will be stored. + +# Returns +- `Cint`: An integer status code indicating the success or failure. + +# Notes +- To avoid memory leaks the memory allocated for members of `result` should be freed with `free_members` function or equivalent. +""" +function c_jet_reconstruct(particles::Ptr{T}, + particles_length::Csize_t, + algorithm::JetAlgorithm.Algorithm, + R::Cdouble, + strategy::RecoStrategy.Strategy, + result::Ptr{C_ClusterSequence{U}}) where {T, U} + try + particles_v = unsafe_wrap_c_array(particles, particles_length) + clusterseq = jet_reconstruct(particles_v; p = nothing, algorithm = algorithm, R = R, + strategy = strategy) + c_clusterseq = C_ClusterSequence{U}(clusterseq) + unsafe_store!(result, c_clusterseq) + catch e + return handle_exception(e) + end + return Cint(StatusCode.OK) +end + +""" + jetreconstruction_jet_reconstruct(particles::Ptr{PseudoJet}, + particles_length::Csize_t, + algorithm::JetAlgorithm.Algorithm, + R::Cdouble, + strategy::RecoStrategy.Strategy, + result::Ptr{C_ClusterSequence{PseudoJet}})::Cint + +C-binding for `jet_reconstruct`. + +# Arguments +- `particles::Ptr{PseudoJet}`: Pointer to an array of pseudojet objects used for jet reconstruction. +- `particles_length::Csize_t`: The length of `particles`` array. +- `algorithm::JetAlgorithm.Algorithm`: The algorithm to use for jet reconstruction. +- `R::Cdouble`: The jet radius parameter.. +- `strategy::RecoStrategy.Strategy`: The jet reconstruction strategy to use. +- `result::Ptr{C_ClusterSequence{PseudoJet}}`: A pointer to which a cluster sequence containing the reconstructed jets and the merging history will be stored. + +# Returns +- `Cint`: An integer status code indicating the success or failure. + +# Notes +- To avoid memory leaks the memory allocated for members of `result` should be freed with `jetreconstruction_ClusterSequence_free_members_`. +""" +Base.@ccallable function jetreconstruction_jet_reconstruct(particles::Ptr{PseudoJet}, + particles_length::Csize_t, + algorithm::JetAlgorithm.Algorithm, + R::Cdouble, + strategy::RecoStrategy.Strategy, + result::Ptr{C_ClusterSequence{PseudoJet}})::Cint + return c_jet_reconstruct(particles, particles_length, algorithm, R, strategy, result) +end + +""" + struct C_JetsResult{T} + +A C-compatible wrapper for array-like data + +# Fields +- `data::Ptr{T}`: A pointer to the data of type `T`. +- `length::Csize_t`: The length of the data. +""" +struct C_JetsResult{T} + data::Ptr{T} + length::Csize_t +end + +""" + free_members(ptr::Ptr{C_JetsResult{T}}) where {T} + +Internal function for freeing dynamically allocated `C_JetsResult` members memory. + +# Arguments +- `ptr::Ptr{C_JetsResult{T}}`: A pointer to a `C_JetsResult` structure. +""" +function free_members(ptr::Ptr{C_JetsResult{T}}) where {T} + if ptr != C_NULL + result = unsafe_load(ptr) + Libc.free(result.data) + # Struct is immutable so pointer can't be assigned NULL and length updated to zero (without making a copy) + end +end + +""" + jetreconstruction_JetsResult_free_members_(ptr::Ptr{C_JetsResult{PseudoJet}})::Cvoid + +C-binding for freeing the members of a `C_JetsResult{PseudoJet}` object pointed to by `ptr`. + +# Arguments +- `ptr::Ptr{C_JetsResult{PseudoJet}}`: A pointer to the `C_JetsResult{PseudoJet}` structure whose members are to be freed. + +# Returns +- `Cvoid`: This function does not return any value. +""" +Base.@ccallable function jetreconstruction_JetsResult_free_members_(ptr::Ptr{C_JetsResult{PseudoJet}})::Cvoid + free_members(ptr) + return nothing +end + +""" + jets_selection(selector, clustersequence::Ptr{C_ClusterSequence{T}}, + result::Ptr{C_JetsResult{U}}; kwargs...)::Cint where {T, U} + +An internal helper function for calling calling functions selecting jets from a given cluster sequence. + +# Arguments +- `selector`: A function that takes a `ClusterSequence{T}` and returns a list of jets. +- `clustersequence::Ptr{C_ClusterSequence{T}}`: A pointer to a C-compatible cluster sequence. +- `result::Ptr{C_JetsResult{U}}`: A pointer to a C-compatible structure where the result will be stored. +- `kwargs...`: Additional keyword arguments to be passed to the selector function. + +# Returns +- `Cint`: An integer status code indicating the success or failure. +""" +function jets_selection(selector, clustersequence::Ptr{C_ClusterSequence{T}}, + result::Ptr{C_JetsResult{U}}; kwargs...)::Cint where {T, U} + try + c_clusterseq = unsafe_load(clustersequence) + clusterseq = ClusterSequence{T}(c_clusterseq) + jets_result = selector(clusterseq; T = U, kwargs...) + c_results = C_JetsResult{U}(make_c_array(jets_result)...) + unsafe_store!(result, c_results) + catch e + return handle_exception(e) + end + return Cint(StatusCode.OK) +end + +""" + jetreconstruction_exclusive_jets_dcut(clustersequence::Ptr{C_ClusterSequence{PseudoJet}}, + dcut::Cdouble, + result::Ptr{C_JetsResult{PseudoJet}}) -> Cint + +C-binding for `exclusive_jets` with a cut on the maximum distance parameter. + +# Arguments +- `clustersequence::Ptr{C_ClusterSequence{PseudoJet}}`: A pointer to the cluster sequence object containing the clustering history and jets. +- `dcut::Cdouble`: The distance parameter used to define the exclusive jets. +- `result::Ptr{C_JetsResult{PseudoJet}}`: A pointer to the results. + +# Returns +- `Cint`: An integer status code indicating the success or failure. + +# Notes +- To avoid memory leaks the memory allocated for members of `result` should be freed with `jetreconstruction_JetsResult_free_members_`. +""" +Base.@ccallable function jetreconstruction_exclusive_jets_dcut(clustersequence::Ptr{C_ClusterSequence{PseudoJet}}, + dcut::Cdouble, + result::Ptr{C_JetsResult{PseudoJet}})::Cint + return jets_selection(exclusive_jets, clustersequence, result; dcut = dcut) +end + +""" + jetreconstruction_exclusive_jets_njets(clustersequence::Ptr{C_ClusterSequence{PseudoJet}}, + njets::Csize_t, + result::Ptr{C_JetsResult{PseudoJet}}) -> Cint + +C-binding for `exclusive_jets` with a specific number of jets. + +# Arguments +- `clustersequence::Ptr{C_ClusterSequence{PseudoJet}}`: A pointer to the cluster sequence object containing the clustering history and jets. +- `njets::Csize_t`: The number of exclusive jets to be calculated. +- `result::Ptr{C_JetsResult{PseudoJet}}`: A pointer to the results. + +# Returns +- `Cint`: An integer status code indicating the success or failure. + +# Notes +- To avoid memory leaks the memory allocated for members of `result` should be freed with `jetreconstruction_JetsResult_free_members_`. +""" +Base.@ccallable function jetreconstruction_exclusive_jets_njets(clustersequence::Ptr{C_ClusterSequence{PseudoJet}}, + njets::Csize_t, + result::Ptr{C_JetsResult{PseudoJet}})::Cint + return jets_selection(exclusive_jets, clustersequence, result; njets = njets) +end + +""" + jetreconstruction_inclusive_jets(clustersequence::Ptr{C_ClusterSequence{PseudoJet}}, + ptmin::Cdouble, + result::Ptr{C_JetsResult{PseudoJet}}) -> Cint + +C-binding for `inclusive_jets`. + +# Arguments +- `clustersequence::Ptr{C_ClusterSequence{PseudoJet}}`: A pointer to the cluster sequence object containing the clustering history and jets. +- `ptmin::Cdouble`: The minimum transverse momentum (pt) threshold for the inclusive jets. +- `result::Ptr{C_JetsResult{PseudoJet}}`: A pointer to the results. + +# Returns +- `Cint`: An integer status code indicating the success or failure. + +# Notes +- To avoid memory leaks the memory allocated for members of `result` should be freed with `jetreconstruction_JetsResult_free_members_`. +""" +Base.@ccallable function jetreconstruction_inclusive_jets(clustersequence::Ptr{C_ClusterSequence{PseudoJet}}, + ptmin::Cdouble, + result::Ptr{C_JetsResult{PseudoJet}})::Cint + return jets_selection(inclusive_jets, clustersequence, result; ptmin = ptmin) +end + +end # module C_JetReconstruction diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl index 9929187..f7ff618 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -89,4 +89,7 @@ export jetsplot, animatereco include("JSONresults.jl") export FinalJet, FinalJets +# C bindings sub-module +include("C_JetReconstruction/C_JetReconstruction.jl") + end