add submodule for c api, library compilation with PackageCompiler

.gitignore

@@ -36,3 +36,6 @@ benchmark/*
 /profile/*
 /statprof/*
 /debug/*
+
+# Compiled packages
+JetReconstructionCompiled

compile/Project.toml

@@ -0,0 +1,11 @@
+[deps]
+JetReconstruction = "44e8cb2c-dfab-4825-9c70-d4808a591196"
+PackageCompiler = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d"
+
+[sources]
+JetReconstruction = {path = ".."}
+
+[compat]
+JetReconstruction = "0.4"
+PackageCompiler = "2"
+julia = "1.11"

compile/README.md

@@ -0,0 +1,65 @@
+# Compiling JetReconstruction.jl to a C-library
+
+Minimal C bindings for JetReconstruction.jl
+
+- [C-header](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;
+  jetreconstruction_jet_reconstruct(particles, len, algorithm, R, strategy,
+                                    &cluster_seq);
+
+  /*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
+```
+

compile/build.jl

@@ -0,0 +1,16 @@
+using PackageCompiler
+
+output_directory = joinpath(splitdir(@__DIR__) |> first, "JetReconstructionCompiled")
+
+@info "Creating library in $output_directory"
+PackageCompiler.create_library(".", output_directory;
+                               lib_name = "jetreconstruction",
+                               header_files = [joinpath(@__DIR__, "include",
+                                                        "JetReconstruction.h")],
+                               #       precompile_execution_file = [joinpath(@__DIR__,
+                               #                                             "precompile_execution.jl")],
+                               #       precompile_statements_file = [jointpath(@__DIR__,
+                               #                                               "precompile_statements.jl")],
+                               incremental = false,
+                               filter_stdlibs = true,
+                               force = true)

compile/include/JetReconstruction.h

@@ -0,0 +1,145 @@
+/*
+ * SPDX-FileCopyrightText: 2022-2024 JetReconstruction.jl authors, CERN
+ * SPDX-License-Identifier: MIT
+ */
+
+#include <stddef.h>
+
+/*
+Enumeration representing different jet algorithms used in
+the JetReconstruction module.
+*/
+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;
+
+/*
+Scoped enumeration (using EnumX) 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;
+
+/* The `PseudoJet` struct represents 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;
+
+int jetreconstruction_PseudoJet_init(jetreconstruction_PseudoJet *ptr,
+                                     double px, double py, double pz, double E);
+/*
+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 (NonexistentParent if this jet is an original
+                     particle) */
+  long parent2;    /* Index in history where second parent of this jet was
+                     created (NonexistentParent if this jet is an original
+                     particle); 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 `jetp_index=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;
+
+/*
+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 (note that
+                               this value is always stored as a Float64 to be
+                    type stable) */
+  double R;      /* The R parameter used for the clustering algorithm. */
+  jetreconstruction_RecoStrategy
+      strategy; /* The strategy used for clustering. */
+  jetreconstruction_PseudoJet
+      *jets; /* The actual jets in the cluster sequence, which are of type `T
+  <: FourMomentum`. */
+  size_t jets_length;  /* Length of jets. */
+  long n_initial_jets; /* The initial number of particles used for exclusive
+                         jets. */
+  jetreconstruction_HistoryElement
+      *history;          /* 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; /* Length of history. */
+  double Qtot;           /* The total energy of the event. */
+} jetreconstruction_ClusterSequence;
+
+void jetreconstruction_ClusterSequence_free_members_(
+    jetreconstruction_ClusterSequence *ptr);
+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;
+}
+
+int jetreconstruction_jet_reconstruct(
+    const jetreconstruction_PseudoJet *particles, size_t particles_length,
+    jetreconstruction_JetAlgorithm algorithm, double R,
+    jetreconstruction_RecoStrategy strategy,
+    jetreconstruction_ClusterSequence *result);
+
+typedef struct {
+  jetreconstruction_PseudoJet *data;
+  size_t length;
+} jetreconstruction_JetsResult;
+
+void jetreconstruction_JetsResult_free_members_(
+    jetreconstruction_JetsResult *ptr);
+static inline void
+jetreconstruction_JetsResult_free_members(jetreconstruction_JetsResult *ptr) {
+  jetreconstruction_JetsResult_free_members_(ptr);
+  ptr->data = NULL;
+  ptr->length = 0;
+}
+
+int jetreconstruction_exclusive_jets_dcut(
+    const jetreconstruction_ClusterSequence *clustersequence, double dcut,
+    jetreconstruction_JetsResult *result);
+
+int jetreconstruction_exclusive_jets_njets(
+    const jetreconstruction_ClusterSequence *clustersequence, size_t njets,
+    jetreconstruction_JetsResult *result);
+
+int jetreconstruction_inclusive_jets(
+    const jetreconstruction_ClusterSequence *clustersequence, double ptmin,
+    jetreconstruction_JetsResult *result);

compile/precompile_execution.jl

@@ -0,0 +1 @@
+

compile/precompile_statements.jl

@@ -0,0 +1 @@
+

compile/test/jetreconstruction_test.c

@@ -0,0 +1,67 @@
+#include "JetReconstruction.h"
+#include "julia_init.h"
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+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);
+    }
+  }
+}
+
+int main(int argc, char *argv[]) {
+  init_julia(argc, argv);
+  size_t len = 2;
+  jetreconstruction_PseudoJet particles[2];
+  jetreconstruction_PseudoJet_init(&particles[0], 0.0, 1.0, 2.0, 3.0);
+  jetreconstruction_PseudoJet_init(&particles[1], 1.0, 2.0, 3.0, 4.0);
+
+  jetreconstruction_JetAlgorithm algorithm = JETRECONSTRUCTION_JETALGORITHM_CA;
+  double R = 3.0;
+  jetreconstruction_RecoStrategy strategy = JETRECONSTRUCTION_RECOSTRATEGY_BEST;
+
+  jetreconstruction_ClusterSequence cluster_seq;
+  jetreconstruction_jet_reconstruct(particles, len, algorithm, R, strategy,
+                                    &cluster_seq);
+
+  printClusterSequence(&cluster_seq);
+  jetreconstruction_JetsResult result;
+  jetreconstruction_exclusive_jets_njets(&cluster_seq, 2, &result);
+  jetreconstruction_JetsResult_free_members(&result);
+  jetreconstruction_ClusterSequence_free_members(&cluster_seq);
+  shutdown_julia(0);
+  return 0;
+}