diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 9be7959..de5d879 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -18,7 +18,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.8'
+          - '1.9'
           - '1'
           - 'nightly'
         os:
diff --git a/LICENSE b/LICENSE
index ab061ea..64a4127 100644
--- a/LICENSE
+++ b/LICENSE
@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2022 Atell Krasnopolski
+Copyright (c) 2022-2024 The Authors, CERN
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal
diff --git a/Project.toml b/Project.toml
index 677ece7..9394b00 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,26 +1,31 @@
 name = "JetReconstruction"
 uuid = "44e8cb2c-dfab-4825-9c70-d4808a591196"
 authors = ["Atell Krasnopolski <delta_atell@protonmail.com>", "Graeme A Stewart <graeme.andrew.stewart@cern.ch", "Philippe Gras <philippe.gras@cern.ch>"]
-version = "0.2.0"
+version = "0.3.0"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
-ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
-Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
-FlameGraphs = "08572546-2f56-4bcf-ba4e-bab62c3a3f89"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
 LorentzVectorHEP = "f612022c-142a-473f-8cfd-a09cf3793c6c"
 LorentzVectors = "3f54b04b-17fc-5cd4-9758-90c048d965e3"
-Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
-ProfileSVG = "132c30aa-f267-4189-9183-c8a63c7e05e6"
-StatProfilerHTML = "a8a75453-ed82-57c9-9e16-4cd1196ecbf5"
+
+[weakdeps]
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+
+[extensions]
+JetVisualisation = "Makie"
 
 [compat]
+Accessors = "0.1.36"
+EnumX = "1.0.4"
+JSON = "0.21.4"
+Makie = "0.21"
 LorentzVectorHEP = "0.1.6"
-julia = "1.8"
+LorentzVectors = "0.4.3"
+julia = "1.9"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/README.md b/README.md
index 8e18c8b..40640d4 100644
--- a/README.md
+++ b/README.md
@@ -10,7 +10,7 @@ Algorithms used are based on the C++ FastJet package (<https://fastjet.fr>,
 [hep-ph/0512210](https://arxiv.org/abs/hep-ph/0512210),
 [arXiv:1111.6097](https://arxiv.org/abs/1111.6097)), reimplemented natively in Julia.
 
-The algorithms include anti-${k}_\text{T}$, Cambridge/Aachen and inclusive $k_\text{T}$.
+The algorithms include anti-$`{k}_\text{T}`$, Cambridge/Aachen and inclusive $`k_\text{T}`$.
 
 ### Interface
 
@@ -71,9 +71,11 @@ plain_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +)
 
 Note that there is no `strategy` option in these interfaces.
 
-### Example
+### Examples
 
-See the `examples/jetreco.jl` script for a full example of how to call jet reconstruction.
+In the examples directory there are a number of example scripts.
+
+See the `jetreco.jl` script for an example of how to call jet reconstruction.
 
 ```sh
 julia --project=. examples/jetreco.jl --maxevents=100 --nsamples=1 --strategy=N2Plain test/data/events.hepmc3
@@ -82,21 +84,40 @@ julia --project=. examples/jetreco.jl --maxevents=100 --nsamples=1 --strategy=N2
 ...
 ```
 
-The example also shows how to use `JetReconstruction.HepMC3` to read HepMC3 ASCII files (via the `read_final_state_particles()` wrapper).
+There are options to explicitly set the algorithm (use `--help` to see these).
+
+The example also shows how to use `JetReconstruction.HepMC3` to read HepMC3
+ASCII files (via the `read_final_state_particles()` wrapper).
+
+Further examples, which show visualisation, timing measurements, profiling, etc.
+are given - see the `README.md` file in the examples directory.
+
+Note that due to additional dependencies the `Project.toml` file for the
+examples is different from the package itself.
 
 ### Plotting
 
-![illustration](img/illustration.jpeg)
+![illustration](img/jetvis.png)
+
+To visualise the clustered jets as a 3d bar plot (see illustration above) we now
+use `Makie.jl`. See the `jetsplot` function in `ext/JetVisualisation.jl` and its
+documentation for more. There are two worked examples in the `examples`
+directory.
 
-To visualise the clustered jets as a 3d bar plot (see illustration above) we now use `Makie.jl`. See the `jetsplot` function and its documentation for more. 
+The plotting code is a package extension and will load if the one of the `Makie`
+modules is loaded in the environment.
 
 ### Serialisation
 
-The package also provides methods such as `loadjets`, `loadjets!`, and `savejets` that one can use to save and load objects on/from disk easily in a very flexible format. See documentation for more.
+The package also provides methods such as `loadjets`, `loadjets!`, and
+`savejets` that one can use to save and load objects on/from disk easily in a
+very flexible format. See documentation for more.
 
 ## Reference
 
-Although it has been developed further since the CHEP2023 conference, the CHEP conference proceedings, [arXiv:2309.17309](https://arxiv.org/abs/2309.17309), should be cited if you use this package:
+Although it has been developed further since the CHEP2023 conference, the CHEP
+conference proceedings, [arXiv:2309.17309](https://arxiv.org/abs/2309.17309),
+should be cited if you use this package:
 
 ```bibtex
 @misc{stewart2023polyglot,
diff --git a/examples/.gitignore b/examples/.gitignore
new file mode 100644
index 0000000..2a92ec2
--- /dev/null
+++ b/examples/.gitignore
@@ -0,0 +1,3 @@
+# Ignore visualisation outputs
+*.pdf
+*.png
diff --git a/examples/Project.toml b/examples/Project.toml
new file mode 100644
index 0000000..fd546ef
--- /dev/null
+++ b/examples/Project.toml
@@ -0,0 +1,12 @@
+[deps]
+ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
+FlameGraphs = "08572546-2f56-4bcf-ba4e-bab62c3a3f89"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+JetReconstruction = "44e8cb2c-dfab-4825-9c70-d4808a591196"
+LorentzVectorHEP = "f612022c-142a-473f-8cfd-a09cf3793c6c"
+Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
+ProfileSVG = "132c30aa-f267-4189-9183-c8a63c7e05e6"
+StatProfilerHTML = "a8a75453-ed82-57c9-9e16-4cd1196ecbf5"
+WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008"
diff --git a/examples/README.md b/examples/README.md
new file mode 100644
index 0000000..3f4f3e8
--- /dev/null
+++ b/examples/README.md
@@ -0,0 +1,31 @@
+# JetReconstruction.jl Examples
+
+This directory has a number of example files that show how to used the
+`JetReconstruction.jl` package.
+
+Because of extra dependencies in these scripts, one must use the `Project.toml`
+file in this directory.
+
+## `jetreco.jl`
+
+This is a basic jet reconstruction example that shows how to call the package to
+perform a jet reconstruction, with different algorithms and (optionally)
+strategy, producing exclusive and inclusive jet selections.
+
+## `instrumented-jetreco.jl`
+
+This is a more sophisticated example that allows performance measurements to be
+made of the reconstruction, as well as profiling (flamegraphs and memory
+profiling).
+
+## `visualise-jets.jl`
+
+This script will produce a PNG/PDF showing the results of a jet reconstruction.
+This is a 3D plot where all the initial energy deposits are visualised, with
+colours that indicate in which final cluster the deposit ended up in.
+
+## `visualise-jets.ipynb`
+
+This notebook will produce a visualisation of jet reconstruction in the browser.
+This is a 3D plot where all the initial energy deposits are visualised, with
+colours that indicate in which final cluster the deposit ended up in.
diff --git a/examples/instrumented-jetreco.jl b/examples/instrumented-jetreco.jl
new file mode 100644
index 0000000..11dfe88
--- /dev/null
+++ b/examples/instrumented-jetreco.jl
@@ -0,0 +1,284 @@
+#! /usr/bin/env julia
+"""
+Example of running jet reconstruction on a HepMC3 file, which also
+produces timing measurements.
+
+Options are available to profile the code (flamegraphs and allocations).
+"""
+
+using FlameGraphs: FlameGraphs
+using ProfileSVG: ProfileSVG
+
+using ArgParse
+using Profile
+using Colors
+using StatProfilerHTML
+using Logging
+using JSON
+
+using LorentzVectorHEP
+using JetReconstruction
+
+# Parsing for algorithm and strategy enums
+include(joinpath(@__DIR__, "parse-options.jl"))
+
+function profile_code(profile, jet_reconstruction, events, niters; R = 0.4, p = -1, strategy = RecoStrategy.N2Tiled)
+    Profile.init(n = 5 * 10^6, delay = 0.00001)
+    profile_events(events) = begin
+        for evt in events
+            jet_reconstruction(evt, R = R, p = p, strategy = strategy)
+        end
+    end
+    profile_events(events[1:1])
+    @profile for i ∈ 1:niters
+        profile_events(events)
+    end
+    profile_path = joinpath("profile", profile, "profsvg.svg")
+    mkpath(dirname(profile_path))
+    statprofilehtml(path = dirname(profile_path))
+    fcolor = FlameGraphs.FlameColors(
+        reverse(colormap("Blues", 15))[1:5],
+        colorant"slategray4",
+        colorant"gray95",
+        reverse(colormap("Reds", 15))[1:5],
+        reverse(sequential_palette(39, 10; s = 38, b = 2))[1:5],#brownish pallette
+    )
+    ProfileSVG.save(
+        fcolor,
+        profile_path,
+        combine = true,
+        timeunit = :ms,
+        font = "Arial, Helvetica, sans-serif",
+    )
+    println(
+        "Flame graph from ProfileSVG.jl at file://",
+        abspath(profile_path),
+        "\n",
+        """
+        \tRed tint:          Runtime dispatch
+        \tBrown/yellow tint: Garbage collection
+        \tBlue tint:         OK
+        """,
+    )
+end
+"""
+Top level call funtion for demonstrating the use of jet reconstruction
+
+This uses the "jet_reconstruct" wrapper, so algorithm switching
+happens inside the JetReconstruction package itself.
+
+Some other ustilities are also supported here, such as profiling and
+serialising the reconstructed jet outputs.
+"""
+function jet_process(
+    events::Vector{Vector{PseudoJet}};
+    distance::Real = 0.4,
+    algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt,
+    ptmin::Real = 5.0,
+    dcut = nothing,
+    njets = nothing,
+    strategy::RecoStrategy.Strategy,
+    nsamples::Integer = 1,
+    gcoff::Bool = false,
+    profile = nothing,
+    alloc::Bool = false,
+    dump::Union{String, Nothing} = nothing,
+)
+    @info "Will process $(size(events)[1]) events"
+
+    # Map algorithm to power
+    power = JetReconstruction.algorithm2power[algorithm]
+
+    # If we are dumping the results, setup the JSON structure
+    if !isnothing(dump)
+        jet_collection = FinalJets[]
+    end
+
+    # Warmup code if we are doing a multi-sample timing run
+    if nsamples > 1 || !isnothing(profile)
+        @info "Doing initial warm-up run"
+        for event in events
+            _ = inclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin)
+        end
+    end
+
+    if !isnothing(profile)
+        profile_code(profile, jet_reconstruct, events, nsamples, R = distance, p = power, strategy = strategy)
+        return nothing
+    end
+
+    if alloc
+        println("Memory allocation statistics:")
+        @timev for event in events
+            _ = inclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin)
+        end
+        return nothing
+    end
+
+    # Now setup timers and run the loop
+    GC.gc()
+    cummulative_time = 0.0
+    cummulative_time2 = 0.0
+    lowest_time = typemax(Float64)
+    for irun ∈ 1:nsamples
+        gcoff && GC.enable(false)
+        t_start = time_ns()
+        for (ievt, event) in enumerate(events)
+            if !isnothing(njets)
+                finaljets = exclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), njets = njets)
+            elseif !isnothing(dcut)
+                finaljets = exclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), dcut = dcut)
+            else
+                finaljets = inclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin)
+            end
+            # Only print the jet content once
+            if irun == 1
+                @info begin
+                    jet_output = "Event $(ievt)\n"
+                    sort!(finaljets, by = x -> pt(x), rev = true)
+                    for (ijet, jet) in enumerate(finaljets)
+                        jet_output *= " $(ijet) - $(jet)\n"
+                    end
+                    "$(jet_output)"
+                end
+                if !isnothing(dump)
+                    push!(jet_collection, FinalJets(ievt, finaljets))
+                end
+            end
+        end
+        t_stop = time_ns()
+        gcoff && GC.enable(true)
+        dt_μs = convert(Float64, t_stop - t_start) * 1.e-3
+        if nsamples > 1
+            @info "$(irun)/$(nsamples) $(dt_μs)"
+        end
+        cummulative_time += dt_μs
+        cummulative_time2 += dt_μs^2
+        lowest_time = dt_μs < lowest_time ? dt_μs : lowest_time
+    end
+
+    mean = cummulative_time / nsamples
+    cummulative_time2 /= nsamples
+    if nsamples > 1
+        sigma = sqrt(nsamples / (nsamples - 1) * (cummulative_time2 - mean^2))
+    else
+        sigma = 0.0
+    end
+    mean /= length(events)
+    sigma /= length(events)
+    lowest_time /= length(events)
+    # Why also record the lowest time? 
+    # 
+    # The argument is that on a "busy" machine, the run time of an application is
+    # always TrueRunTime+Overheads, where Overheads is a nuisance parameter that
+    # adds jitter, depending on the other things the machine is doing. Therefore
+    # the minimum value is (a) more stable and (b) reflects better the intrinsic
+    # code performance.
+    println("Processed $(length(events)) events $(nsamples) times")
+    println("Average time per event $(mean) ± $(sigma) μs")
+    println("Lowest time per event $lowest_time μs")
+
+    if !isnothing(dump)
+        open(dump, "w") do io
+            JSON.print(io, jet_collection, 2)
+        end
+    end
+end
+
+parse_command_line(args) = begin
+    s = ArgParseSettings(autofix_names = true)
+    @add_arg_table! s begin
+        "--maxevents", "-n"
+        help = "Maximum number of events to read. -1 to read all events from the  file."
+        arg_type = Int
+        default = -1
+
+        "--skip", "-s"
+        help = "Number of events to skip at beginning of the file."
+        arg_type = Int
+        default = 0
+
+        "--ptmin"
+        help = "Minimum p_t for final inclusive jets (energy unit is the same as the input clusters, usually GeV)"
+        arg_type = Float64
+        default = 5.0
+
+        "--exclusive-dcut"
+        help = "Return all exclusive jets where further merging would have d>d_cut"
+        arg_type = Float64
+
+        "--exclusive-njets"
+        help = "Return all exclusive jets once clusterisation has produced n jets"
+        arg_type = Int
+
+        "--distance", "-R"
+        help = "Distance parameter for jet merging"
+        arg_type = Float64
+        default = 0.4
+
+        "--algorithm", "-A"
+        help = """Algorithm to use for jet reconstruction: $(join(JetReconstruction.AllJetRecoAlgorithms, ", "))"""
+        arg_type = JetAlgorithm.Algorithm
+        default = JetAlgorithm.AntiKt
+
+        "--strategy", "-S"
+        help = """Strategy for the algorithm, valid values: $(join(JetReconstruction.AllJetRecoStrategies, ", "))"""
+        arg_type = RecoStrategy.Strategy
+        default = RecoStrategy.Best
+
+        "--nsamples", "-m"
+        help = "Number of measurement points to acquire."
+        arg_type = Int
+        default = 1
+
+        "--gcoff"
+        help = "Turn off Julia garbage collector during each time measurement."
+        action = :store_true
+
+        "--profile"
+        help = """Profile code and generate a flame graph; the results 
+        will be stored in 'profile/PROFILE' subdirectory, where PROFILE is the value
+        given to this option)"""
+
+        "--alloc"
+        help = "Provide memory allocation statistics."
+        action = :store_true
+
+        "--dump"
+        help = "Write list of reconstructed jets to a JSON formatted file"
+
+        "--info"
+        help = "Print info level log messages"
+        action = :store_true
+
+        "--debug"
+        help = "Print debug level log messages"
+        action = :store_true
+
+        "file"
+        help = "HepMC3 event file in HepMC3 to read."
+        required = true
+    end
+    return parse_args(args, s; as_symbols = true)
+end
+
+main() = begin
+    args = parse_command_line(ARGS)
+    if args[:debug]
+        logger = ConsoleLogger(stdout, Logging.Debug)
+    elseif args[:info]
+        logger = ConsoleLogger(stdout, Logging.Info)
+    else
+        logger = ConsoleLogger(stdout, Logging.Warn)
+    end
+    global_logger(logger)
+    events::Vector{Vector{PseudoJet}} =
+        read_final_state_particles(args[:file], maxevents = args[:maxevents], skipevents = args[:skip])
+    jet_process(events, distance = args[:distance], algorithm = args[:algorithm], strategy = args[:strategy],
+        ptmin = args[:ptmin], dcut = args[:exclusive_dcut], njets = args[:exclusive_njets],
+        nsamples = args[:nsamples], gcoff = args[:gcoff], profile = args[:profile],
+        alloc = args[:alloc], dump = args[:dump])
+    nothing
+end
+
+main()
diff --git a/examples/jetreco.jl b/examples/jetreco.jl
index 65d65e9..811345c 100644
--- a/examples/jetreco.jl
+++ b/examples/jetreco.jl
@@ -1,12 +1,9 @@
 #! /usr/bin/env julia
 """
-Wrapper to run jet reco code feeding in the standard set of HepMC events that
-are used for benchmarking jet reco.
+Simple example of a jet reconstruction code that reads in a text HepMC3 file
+and performs standard jet reconstruction on the final state particles.
 """
 
-using FlameGraphs: FlameGraphs
-using ProfileSVG: ProfileSVG
-
 using ArgParse
 using Profile
 using Colors
@@ -17,274 +14,137 @@ using JSON
 using LorentzVectorHEP
 using JetReconstruction
 
-function profile_code(profile, jet_reconstruction, events, niters; R=0.4, p=-1, strategy=RecoStrategy.N2Tiled)
-	Profile.init(n = 5*10^6, delay = 0.00001)
-	profile_events(events) = begin
-		for evt in events
-			jet_reconstruction(evt, R=R, p=p, strategy=strategy)
-		end
-	end
-	profile_events(events[1:1])
-	@profile for i ∈ 1:niters
-		profile_events(events)
-	end
-	profile_path = joinpath("profile", profile, "profsvg.svg")
-	mkpath(dirname(profile_path))
-	statprofilehtml(path=dirname(profile_path))
-	fcolor = FlameGraphs.FlameColors(
-		reverse(colormap("Blues", 15))[1:5],
-		colorant"slategray4",
-		colorant"gray95",
-		reverse(colormap("Reds", 15))[1:5],
-		reverse(sequential_palette(39, 10; s = 38, b = 2))[1:5],#brownish pallette
-	)
-	ProfileSVG.save(
-		fcolor,
-		profile_path,
-		combine = true,
-		timeunit = :ms,
-		font = "Arial, Helvetica, sans-serif",
-	)
-	println(
-		"Flame graph from ProfileSVG.jl at file://",
-		abspath(profile_path),
-		"\n",
-		"""
-\tRed tint:          Runtime dispatch
-\tBrown/yellow tint: Garbage collection
-\tBlue tint:         OK
-""",
-	)
-end
+# Parsing for algorithm and strategy enums
+include(joinpath(@__DIR__, "parse-options.jl"))
+
 """
 Top level call funtion for demonstrating the use of jet reconstruction
 
 This uses the "jet_reconstruct" wrapper, so algorithm switching
 happens inside the JetReconstruction package itself.
 
-Some other ustilities are also supported here, such as profiling and
-serialising the reconstructed jet outputs.
+Final jets can be serialised if the "dump" option is given
 """
 function jet_process(
-	events::Vector{Vector{PseudoJet}};
-	distance::Real = 0.4,
-	power::Integer = -1,
-	ptmin::Real = 5.0,
-	dcut = nothing,
-	njets = nothing,
-	strategy::RecoStrategy.Strategy,
-	nsamples::Integer = 1,
-	gcoff::Bool = false,
-	profile = nothing,
-    alloc::Bool = false,
-	dump::Union{String, Nothing} = nothing,
+    events::Vector{Vector{PseudoJet}};
+    distance::Real = 0.4,
+    algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt,
+    ptmin::Real = 5.0,
+    dcut = nothing,
+    njets = nothing,
+    strategy::RecoStrategy.Strategy,
+    dump::Union{String, Nothing} = nothing,
 )
-	@info "Will process $(size(events)[1]) events"
-
-	# If we are dumping the results, setup the JSON structure
-	if !isnothing(dump)
-		jet_collection = FinalJets[]
-	end
+    @info "Will process $(size(events)[1]) events"
 
-	# Warmup code if we are doing a multi-sample timing run
-	if nsamples > 1 || !isnothing(profile)
-		@info "Doing initial warm-up run"
-		for event in events
-			_ = inclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin)
-		end
-	end
+    # If we are dumping the results, setup the JSON structure
+    if !isnothing(dump)
+        jet_collection = FinalJets[]
+    end
 
-	if !isnothing(profile)
-		profile_code(profile, jet_reconstruct, events, nsamples, R = distance, p = power, strategy = strategy)
-		return nothing
-	end
+    # A friendly label for the algorithm and final jet selection
+    if !isnothing(njets)
+        @info "Running exclusive jets with n_jets = $(njets)"
+    elseif !isnothing(dcut)
+        @info "Running exclusive jets with dcut = $(dcut)"
+    else
+        @info "Running inclusive jets with ptmin = $(ptmin)"
+    end
 
-    if alloc
-        println("Memory allocation statistics:")
-        @timev for event in events
-			_ = inclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin)
+    # Map algorithm to power
+    power = JetReconstruction.algorithm2power[algorithm]
+
+    # Now run over each event
+    for (ievt, event) in enumerate(events)
+        if !isnothing(njets)
+            finaljets = exclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), njets = njets)
+        elseif !isnothing(dcut)
+            finaljets = exclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), dcut = dcut)
+        else
+            finaljets = inclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin)
+        end
+        @info begin
+            jet_output = "Event $(ievt)\n"
+            sort!(finaljets, by = x -> pt(x), rev = true)
+            for (ijet, jet) in enumerate(finaljets)
+                jet_output *= " $(ijet) - $(jet)\n"
+            end
+            "$(jet_output)"
+        end
+        if !isnothing(dump)
+            push!(jet_collection, FinalJets(ievt, finaljets))
         end
-        return nothing
     end
 
-	# Now setup timers and run the loop
-	GC.gc()
-	cummulative_time = 0.0
-	cummulative_time2 = 0.0
-	lowest_time = typemax(Float64)
-	for irun ∈ 1:nsamples
-		gcoff && GC.enable(false)
-		t_start = time_ns()
-		for (ievt, event) in enumerate(events)
-			if !isnothing(njets)
-				finaljets = exclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), njets=njets)
-			elseif !isnothing(dcut)
-				finaljets = exclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), dcut=dcut)
-			else
-				finaljets = inclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin)
-			end
-			# Only print the jet content once
-			if irun == 1
-				@info begin
-					jet_output = "Event $(ievt)\n"
-					sort!(finaljets, by = x -> pt(x), rev=true)
-					for (ijet, jet) in enumerate(finaljets)
-						jet_output *= " $(ijet) - $(jet)\n"
-					end
-					"$(jet_output)"
-				end
-				if !isnothing(dump)
-					push!(jet_collection, FinalJets(ievt, finaljets))
-				end
-			end
-		end
-		t_stop = time_ns()
-		gcoff && GC.enable(true)
-		dt_μs = convert(Float64, t_stop - t_start) * 1.e-3
-        if nsamples > 1
-			@info "$(irun)/$(nsamples) $(dt_μs)"
-		end
-		cummulative_time += dt_μs
-		cummulative_time2 += dt_μs^2
-		lowest_time = dt_μs < lowest_time ? dt_μs : lowest_time
-	end
-
-	mean = cummulative_time / nsamples
-	cummulative_time2 /= nsamples
-	if nsamples > 1
-		sigma = sqrt(nsamples / (nsamples - 1) * (cummulative_time2 - mean^2))
-	else
-		sigma = 0.0
-	end
-	mean /= length(events)
-	sigma /= length(events)
-	lowest_time /= length(events)
-	# Why also record the lowest time? 
-	# 
-	# The argument is that on a "busy" machine, the run time of an application is
-	# always TrueRunTime+Overheads, where Overheads is a nuisance parameter that
-	# adds jitter, depending on the other things the machine is doing. Therefore
-	# the minimum value is (a) more stable and (b) reflects better the intrinsic
-	# code performance.
-	println("Processed $(length(events)) events $(nsamples) times")
-	println("Average time per event $(mean) ± $(sigma) μs")
-	println("Lowest time per event $lowest_time μs")
-
-	if !isnothing(dump)
-		open(dump, "w") do io
-			JSON.print(io, jet_collection, 2)
-		end
-	end
+    if !isnothing(dump)
+        open(dump, "w") do io
+            JSON.print(io, jet_collection, 2)
+        end
+    end
 end
 
 parse_command_line(args) = begin
-	s = ArgParseSettings(autofix_names = true)
-	@add_arg_table! s begin
-		"--maxevents", "-n"
-		help = "Maximum number of events to read. -1 to read all events from the  file."
-		arg_type = Int
-		default = -1
-
-		"--skip", "-s"
-		help = "Number of events to skip at beginning of the file."
-		arg_type = Int
-		default = 0
-
-		"--ptmin"
-		help = "Minimum p_t for final inclusive jets (energy unit is the same as the input clusters, usually GeV)"
-		arg_type = Float64
-		default = 5.0
-
-		"--exclusive-dcut"
-		help = "Return all exclusive jets where further merging would have d>d_cut"
-		arg_type = Float64
-
-		"--exclusive-njets"
-		help = "Return all exclusive jets once clusterisation has produced n jets"
-		arg_type = Int
-
-		"--distance", "-R"
-		help = "Distance parameter for jet merging"
-		arg_type = Float64
-		default = 0.4
-
-		"--power"
-		help = "Distance measure momentum power (-1 - antikt; 0 - Cambridge/Aachen; 1 - inclusive k_t)"
-		arg_type = Int
-		default = -1
-
-		"--strategy"
-		help = """Strategy for the algorithm, valid values: $(join(JetReconstruction.AllJetRecoStrategies, ", "))"""
-		arg_type = RecoStrategy.Strategy
-		default = RecoStrategy.Best
-
-		"--nsamples", "-m"
-		help = "Number of measurement points to acquire."
-		arg_type = Int
-		default = 1
-
-		"--gcoff"
-		help = "Turn off Julia garbage collector during each time measurement."
-		action = :store_true
-
-		"--profile"
-		help = """Profile code and generate a flame graph; the results 
-		will be stored in 'profile/PROFILE' subdirectory, where PROFILE is the value
-		given to this option)"""
-
-        "--alloc"
-        help = "Provide memory allocation statistics."
-        action = :store_true
-
-		"--dump"
-		help = "Write list of reconstructed jets to a JSON formatted file"
-
-		"--info"
-		help = "Print info level log messages"
-		action = :store_true
-
-		"--debug"
-		help = "Print debug level log messages"
-		action = :store_true
-
-		"file"
-		help = "HepMC3 event file in HepMC3 to read."
-		required = true
-	end
-	return parse_args(args, s; as_symbols = true)
+    s = ArgParseSettings(autofix_names = true)
+    @add_arg_table! s begin
+        "--maxevents", "-n"
+        help = "Maximum number of events to read. -1 to read all events from the  file."
+        arg_type = Int
+        default = -1
+
+        "--skip", "-s"
+        help = "Number of events to skip at beginning of the file."
+        arg_type = Int
+        default = 0
+
+        "--ptmin"
+        help = "Minimum p_t for final inclusive jets (energy unit is the same as the input clusters, usually GeV)"
+        arg_type = Float64
+        default = 5.0
+
+        "--exclusive-dcut"
+        help = "Return all exclusive jets where further merging would have d>d_cut"
+        arg_type = Float64
+
+        "--exclusive-njets"
+        help = "Return all exclusive jets once clusterisation has produced n jets"
+        arg_type = Int
+
+        "--distance", "-R"
+        help = "Distance parameter for jet merging"
+        arg_type = Float64
+        default = 0.4
+
+        "--algorithm", "-A"
+        help = """Algorithm to use for jet reconstruction: $(join(JetReconstruction.AllJetRecoAlgorithms, ", "))"""
+        arg_type = JetAlgorithm.Algorithm
+        default = JetAlgorithm.AntiKt
+
+        "--strategy", "-S"
+        help = """Strategy for the algorithm, valid values: $(join(JetReconstruction.AllJetRecoStrategies, ", "))"""
+        arg_type = RecoStrategy.Strategy
+        default = RecoStrategy.Best
+
+        "--dump"
+        help = "Write list of reconstructed jets to a JSON formatted file"
+
+        "file"
+        help = "HepMC3 event file in HepMC3 to read."
+        required = true
+    end
+    return parse_args(args, s; as_symbols = true)
 end
 
 
-function ArgParse.parse_item(::Type{RecoStrategy.Strategy}, x::AbstractString)
-	s = tryparse(RecoStrategy.Strategy, x)
-	if s === nothing
-		throw(ErrorException("Invalid value for strategy: $(x)"))
-	end
-	s
-end
-
 main() = begin
-	args = parse_command_line(ARGS)
-	if args[:debug]
-		logger = ConsoleLogger(stdout, Logging.Debug)
-	elseif args[:info]
-		logger = ConsoleLogger(stdout, Logging.Info)
-	else
-		logger = ConsoleLogger(stdout, Logging.Warn)
-	end
-	global_logger(logger)
-	events::Vector{Vector{PseudoJet}} =
-		read_final_state_particles(args[:file], maxevents = args[:maxevents], skipevents = args[:skip])
-	jet_process(events, distance = args[:distance], power = args[:power], strategy = args[:strategy],
-		ptmin = args[:ptmin], dcut = args[:exclusive_dcut], njets = args[:exclusive_njets],
-		nsamples = args[:nsamples], gcoff = args[:gcoff], profile = args[:profile],
-		alloc = args[:alloc], dump = args[:dump])
-	nothing
+    args = parse_command_line(ARGS)
+    logger = ConsoleLogger(stdout, Logging.Info)
+    global_logger(logger)
+    events::Vector{Vector{PseudoJet}} =
+        read_final_state_particles(args[:file], maxevents = args[:maxevents], skipevents = args[:skip])
+    jet_process(events, distance = args[:distance], algorithm = args[:algorithm], strategy = args[:strategy],
+        ptmin = args[:ptmin], dcut = args[:exclusive_dcut], njets = args[:exclusive_njets],
+        dump = args[:dump])
+    nothing
 end
 
-# Running through the debugger in VS Code actually has
-# ARGS[0] = "/some/path/.vscode/extensions/julialang.language-julia-1.47.2/scripts/debugger/run_debugger.jl",
-# so then the program does nothing at all if it only tests for abspath(PROGRAM_FILE) == @__FILE__
-# In addition, deep debugging with Infiltrator needs to start in an interactive session,
-# so execute main() unconditionally
 main()
diff --git a/examples/parse-options.jl b/examples/parse-options.jl
new file mode 100644
index 0000000..c94bdfc
--- /dev/null
+++ b/examples/parse-options.jl
@@ -0,0 +1,20 @@
+"""
+Add `parse_item` code for interpeting `JetAlgorithm.Algorithm` and `RecoStrategy.Strategy` types
+from the command line.
+"""
+
+function ArgParse.parse_item(::Type{RecoStrategy.Strategy}, x::AbstractString)
+	s = tryparse(RecoStrategy.Strategy, x)
+	if s === nothing
+		throw(ErrorException("Invalid value for strategy: $(x)"))
+	end
+	s
+end
+
+function ArgParse.parse_item(::Type{JetAlgorithm.Algorithm}, x::AbstractString)
+	s = tryparse(JetAlgorithm.Algorithm, x)
+	if s === nothing
+		throw(ErrorException("Invalid value for algorithm: $(x)"))
+	end
+	s
+end
\ No newline at end of file
diff --git a/examples/visualise-jets.ipynb b/examples/visualise-jets.ipynb
new file mode 100644
index 0000000..ae7e1bb
--- /dev/null
+++ b/examples/visualise-jets.ipynb
@@ -0,0 +1,153 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Jet Reconstruction Visualisation\n",
+    "\n",
+    "Use Makie to visualise the results of jet reconstruction"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using JetReconstruction\n",
+    "using WGLMakie"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "-1"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "input_file = joinpath(dirname(pathof(JetReconstruction)), \"..\", \"test\", \"data\", \"events.hepmc3\")\n",
+    "event_number = 5\n",
+    "power = JetReconstruction.algorithm2power[JetAlgorithm.AntiKt]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "┌ Info: Total Events: 1\n",
+      "└ @ JetReconstruction /Users/graemes/code/JetReconstruction.jl/src/Utils.jl:25\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "1-element Vector{Vector{PseudoJet}}:\n",
+       " [Pseudojet(px: 0.10735113569266963 py: 0.2004614564282816 pz: 4.422690640225503 E: 4.430731468492096; pt: 0.22739626612557606 eta: 3.5020136230471057 phi: 1.0791380602175356, m: 0.13957000000001873), Pseudojet(px: 0.0011180093991561963 py: 0.012301879085520844 pz: 22.359575220942624 E: 22.365027980969874; pt: 0.01235257762782233 eta: 4.506078564441087 phi: 1.4801641127937162, m: 0.49367999999995554), Pseudojet(px: -0.33118180618811155 py: -0.24231713636363833 pz: -3.25942622170417 E: 3.2881207494491758; pt: 0.4103644518297046 eta: -2.715069658291437 phi: 3.773261626597981, m: 0.1395699999999775), Pseudojet(px: 0.2845031139591477 py: -0.2609196772690977 pz: 0.5824183002026259 E: 0.7125390944710178; pt: 0.3860325113752232 eta: 1.148884934008853 phi: 5.540999045397645, m: 0.13957000000000067), Pseudojet(px: 0.277302069395142 py: 0.03634561140067778 pz: 39.80187826849783 E: 39.80592231168471; pt: 0.27967381207205877 eta: 4.943811191960002 phi: 0.13032576256213585, m: 0.49367999999988865), Pseudojet(px: 0.5585153132771401 py: -0.3956483216555057 pz: -3.612052304986627 E: 3.6789779271494787; pt: 0.6844537600115003 eta: -2.345409125258054 phi: 5.666848658975102, m: 0.13956999999996506), Pseudojet(px: 0.18720823164796177 py: -0.1318330234187746 pz: -0.01476649422315373 E: 0.2685604258116496; pt: 0.22896914215782896 eta: -0.055039390302584666 phi: 5.669642585890234, m: 0.13957000000000006), Pseudojet(px: -0.284009089213228 py: 0.19114033359969312 pz: 0.5990537728620824 E: 0.7039467292094088; pt: 0.34233870637765934 eta: 1.259742297825951 phi: 2.549212870154946, m: 0.13957000000000028), Pseudojet(px: -0.1198632609880622 py: -1.01567218528733 pz: 12.315390709550222 E: 12.358571329368008; pt: 1.0227204844438351 eta: 3.1740560094308856 phi: 4.59491858779022, m: 0.1395699999999233), Pseudojet(px: 0.14193817088911215 py: -0.10210725343936139 pz: -491.6296473818766 E: 491.63057629580845; pt: 0.17484946542748359 eta: -6.936184139290346 phi: 5.659571757044644, m: 0.9395699999616959)  …  Pseudojet(px: -0.22734940597450876 py: -0.19391176346017772 pz: 1.2342343780389635 E: 1.269891343516611; pt: 0.29881352781491993 eta: 2.125875213203082 phi: 3.8477817997568358, m: -4.012260442327338e-8), Pseudojet(px: -0.523471846359122 py: -0.48320925068178905 pz: 2.179296432537861 E: 2.2927814755766254; pt: 0.7124001360717755 eta: 1.836968693576574 phi: 3.887016739844139, m: -6.664001874625056e-8), Pseudojet(px: -0.2312767181096859 py: -0.21992677138349426 pz: 0.4294329015812787 E: 0.5350414208930659; pt: 0.31914997275693896 eta: 1.1059220830216552 phi: 3.9018412914459106, m: -1.235539015514925e-8), Pseudojet(px: -0.11069390032673052 py: -0.08757514264029723 pz: 0.3836884534292784 E: 0.4088268269976469; pt: 0.14114724644148235 eta: 1.7254081403827335 phi: 3.8109120397691947, m: -1.3301947021005397e-8), Pseudojet(px: -0.10329088004039041 py: -0.06428883857518583 pz: 0.8355555651306337 E: 0.8558241337948111; pt: 0.12166371959160473 eta: 2.212114260143481 phi: 3.6983242717426164, m: 0.13957000000000566), Pseudojet(px: -0.47893890239446585 py: -0.7722751051190008 pz: 7.740016906120406 E: 7.794429600820265; pt: 0.9087306037618517 eta: 2.8271088536617697 phi: 4.157273223909725, m: 0.13957000000024747), Pseudojet(px: 0.26463547776785407 py: 0.1378250615977995 pz: -0.33116508213450935 E: 0.44575553313785565; pt: 0.29837507217905723 eta: -0.9569868503353292 phi: 0.480157443201187, m: -3.725290298461914e-9), Pseudojet(px: 0.12023014344988193 py: -0.02895697531429787 pz: -0.14561313580629956 E: 0.19104182561071778; pt: 0.12366807920127183 eta: -1.0014573565077023 phi: 6.046840349964005, m: 4.562530187486071e-9), Pseudojet(px: 0.23193210651528592 py: -0.20537892641762634 pz: -0.18706512331287167 E: 0.3618927545684421; pt: 0.3097951023645705 eta: -0.5721105945608865 phi: 5.55843189731101, m: 9.125060374972142e-9), Pseudojet(px: 0.13205569774138806 py: -0.1001219907287217 pz: -0.23680434557793556 E: 0.2890318640185567; pt: 0.16572000583347526 eta: -1.1546901232363953 phi: 5.634469512021886, m: 5.890201144234053e-9)]"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "# We read only one event\n",
+    "events::Vector{Vector{PseudoJet}} = read_final_state_particles(input_file, maxevents = event_number, skipevents = event_number)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "ClusterSequence(JetReconstruction.JetAlgorithm.AntiKt, JetReconstruction.RecoStrategy.N2Tiled, PseudoJet[Pseudojet(px: 0.10735113569266963 py: 0.2004614564282816 pz: 4.422690640225503 E: 4.430731468492096; pt: 0.22739626612557606 eta: 3.5020136230471057 phi: 1.0791380602175356, m: 0.13957000000001873), Pseudojet(px: 0.0011180093991561963 py: 0.012301879085520844 pz: 22.359575220942624 E: 22.365027980969874; pt: 0.01235257762782233 eta: 4.506078564441087 phi: 1.4801641127937162, m: 0.49367999999995554), Pseudojet(px: -0.33118180618811155 py: -0.24231713636363833 pz: -3.25942622170417 E: 3.2881207494491758; pt: 0.4103644518297046 eta: -2.715069658291437 phi: 3.773261626597981, m: 0.1395699999999775), Pseudojet(px: 0.2845031139591477 py: -0.2609196772690977 pz: 0.5824183002026259 E: 0.7125390944710178; pt: 0.3860325113752232 eta: 1.148884934008853 phi: 5.540999045397645, m: 0.13957000000000067), Pseudojet(px: 0.277302069395142 py: 0.03634561140067778 pz: 39.80187826849783 E: 39.80592231168471; pt: 0.27967381207205877 eta: 4.943811191960002 phi: 0.13032576256213585, m: 0.49367999999988865), Pseudojet(px: 0.5585153132771401 py: -0.3956483216555057 pz: -3.612052304986627 E: 3.6789779271494787; pt: 0.6844537600115003 eta: -2.345409125258054 phi: 5.666848658975102, m: 0.13956999999996506), Pseudojet(px: 0.18720823164796177 py: -0.1318330234187746 pz: -0.01476649422315373 E: 0.2685604258116496; pt: 0.22896914215782896 eta: -0.055039390302584666 phi: 5.669642585890234, m: 0.13957000000000006), Pseudojet(px: -0.284009089213228 py: 0.19114033359969312 pz: 0.5990537728620824 E: 0.7039467292094088; pt: 0.34233870637765934 eta: 1.259742297825951 phi: 2.549212870154946, m: 0.13957000000000028), Pseudojet(px: -0.1198632609880622 py: -1.01567218528733 pz: 12.315390709550222 E: 12.358571329368008; pt: 1.0227204844438351 eta: 3.1740560094308856 phi: 4.59491858779022, m: 0.1395699999999233), Pseudojet(px: 0.14193817088911215 py: -0.10210725343936139 pz: -491.6296473818766 E: 491.63057629580845; pt: 0.17484946542748359 eta: -6.936184139290346 phi: 5.659571757044644, m: 0.9395699999616959)  …  Pseudojet(px: -0.39835711999954054 py: 0.2670425498120344 pz: 0.7823069850508971 E: 0.9340296824136827; pt: 0.4795832758389738 eta: 1.2129464793726437 phi: 2.551037707350037, m: 0.17437634675018016), Pseudojet(px: -0.951550846981664 py: -0.9977730565509817 pz: 2.4230718853428095 E: 3.73223181209067; pt: 1.3787675970846611 eta: 0.7739642228927174 phi: 3.9506982929567993, m: 2.481386074530168), Pseudojet(px: 0.575045650586158 py: -0.5793693320717683 pz: 119.64892585365241 E: 119.65220155609016; pt: 0.8163003878495616 eta: 5.5994728271350915 phi: 5.494041817072272, m: 0.3428308084295518), Pseudojet(px: 0.7430438203838197 py: -0.628041726281588 pz: 141.8604114139369 E: 141.8648142562351; pt: 0.9729082839411632 eta: 5.536755612159961 phi: 5.581467943888408, m: 0.5501335253218786), Pseudojet(px: -0.7419306814989711 py: -1.2366670329604335 pz: -624.3010720899828 E: 624.3038227364076; pt: 1.4421534199108947 eta: -6.512850727112166 phi: 4.17201082701939, m: 1.163900383803333), Pseudojet(px: 0.8187724552223374 py: -1.067476764067555 pz: 231.86864600370112 E: 231.87426160333024; pt: 1.3453233720020452 eta: 5.66076861307267 phi: 5.366692093104101, m: 0.8912348314063477), Pseudojet(px: -0.744694588399963 py: -0.009592334481453038 pz: 208.62889484542603 E: 208.63058856570524; pt: 0.7447563647750817 eta: 6.207268055736434 phi: 3.154472838529462, m: 0.38994718343853185), Pseudojet(px: -0.1192134634148563 py: -0.3533286513298661 pz: 28.166510818816608 E: 28.169333465655836; pt: 0.37289809024711745 eta: 4.950705618562694 phi: 4.386981945468061, m: 0.14129121379207538), Pseudojet(px: 0.23303548004166347 py: -0.13392810235789748 pz: 912.9974388192057 E: 912.9976433113742; pt: 0.26877922456811276 eta: 8.002430623603196 phi: 5.761568088926931, m: 0.5487799398343952), Pseudojet(px: -0.08919052127372945 py: -0.11399069359282013 pz: -13.401707836816184 E: 13.4025408354216; pt: 0.1447370972171675 eta: -5.189519498733623 phi: 4.04844938780284, m: 0.0371358916999187)], 627, JetReconstruction.HistoryElement[JetReconstruction.HistoryElement(-2, -2, 920, 1, 0.0, 0.0), JetReconstruction.HistoryElement(-2, -2, 875, 2, 0.0, 0.0), JetReconstruction.HistoryElement(-2, -2, 990, 3, 0.0, 0.0), JetReconstruction.HistoryElement(-2, -2, 665, 4, 0.0, 0.0), JetReconstruction.HistoryElement(-2, -2, 1247, 5, 0.0, 0.0), JetReconstruction.HistoryElement(-2, -2, 1059, 6, 0.0, 0.0), JetReconstruction.HistoryElement(-2, -2, 1045, 7, 0.0, 0.0), JetReconstruction.HistoryElement(-2, -2, 1224, 8, 0.0, 0.0), JetReconstruction.HistoryElement(-2, -2, 948, 9, 0.0, 0.0), JetReconstruction.HistoryElement(-2, -2, 1216, 10, 0.0, 0.0)  …  JetReconstruction.HistoryElement(248, -1, -3, -3, 9.280205444134888, 9.280205444134888), JetReconstruction.HistoryElement(270, -1, -3, -3, 10.433687082239864, 10.433687082239864), JetReconstruction.HistoryElement(5, -1, -3, -3, 12.784872340278474, 12.784872340278474), JetReconstruction.HistoryElement(1243, -1, -3, -3, 13.842311134655677, 13.842311134655677), JetReconstruction.HistoryElement(117, -1, -3, -3, 23.637588808292218, 23.637588808292218), JetReconstruction.HistoryElement(67, 408, 1251, 1203, 30.63667250108433, 30.63667250108433), JetReconstruction.HistoryElement(1250, -1, -3, -3, 47.73536891403872, 47.73536891403872), JetReconstruction.HistoryElement(11, -1, -3, -3, 65.27369292060092, 65.27369292060092), JetReconstruction.HistoryElement(494, -1, -3, -3, 396.63695275112633, 396.63695275112633), JetReconstruction.HistoryElement(429, -1, -3, -3, 29349.05857739363, 29349.05857739363)], 13000.000000397144)"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "cs = jet_reconstruct(events[1], R = 1.0, p = power)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div class=\"bonito-fragment\" id=\"560b6b52-3454-4de9-a7d0-e3eb34f27d29\" data-jscall-id=\"root\">\n",
+       "  <div>\n",
+       "    <script src=\"http://localhost:9384/assets/57a24233568ab04a755c3e4d5bc6580eb4504a70-Bonito.bundled.js\" type=\"module\"></script>\n",
+       "    <style></style>\n",
+       "    <div></div>\n",
+       "  </div>\n",
+       "  <div>\n",
+       "    <script type=\"module\">    Bonito.lock_loading(() => {\n",
+       "        return Bonito.fetch_binary('http://localhost:9384/assets/aa745754d65135ffedf6b982959baeee68633bca-4487588625566090765.bin').then(msgs=> Bonito.init_session('560b6b52-3454-4de9-a7d0-e3eb34f27d29', msgs, 'root'));\n",
+       "    })\n",
+       "</script>\n",
+       "    <script type=\"module\">    import('http://localhost:9384/assets/fe1a4bb6992a553d47533fa13ef2db822ae1d57f-Websocket.bundled.js').then(WS => {\n",
+       "        WS.setup_connection({proxy_url: 'http://localhost:9384', session_id: '560b6b52-3454-4de9-a7d0-e3eb34f27d29', compression_enabled: false})\n",
+       "    })\n",
+       "</script>\n",
+       "    <div style=\"width: 100%; height: 100%\" data-jscall-id=\"1\">\n",
+       "      <canvas style=\"display: block\" data-jscall-id=\"2\" tabindex=\"0\"></canvas>\n",
+       "    </div>\n",
+       "  </div>\n",
+       "</div>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "jetreco = jetsplot(events[1], cs; Module = WGLMakie)\n",
+    "display(jetreco);"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 1.10.3",
+   "language": "julia",
+   "name": "julia-1.10"
+  },
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.10.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/examples/visualise-jets.jl b/examples/visualise-jets.jl
new file mode 100644
index 0000000..77a7b90
--- /dev/null
+++ b/examples/visualise-jets.jl
@@ -0,0 +1,76 @@
+#! /usr/bin/env julia
+"""Use the visualisation tools to plot the reconstructed jets"""
+
+using ArgParse
+using Logging
+using CairoMakie
+
+using JetReconstruction
+
+# Parsing for algorithm and strategy enums
+include(joinpath(@__DIR__, "parse-options.jl"))
+
+function main()
+
+    s = ArgParseSettings(autofix_names = true)
+    @add_arg_table s begin
+        "--event", "-e"
+        help = "Event in file to visualise"
+        arg_type = Int
+        default = 1
+
+        "--ptmin"
+        help = "Minimum p_t for final inclusive jets (energy unit is the same as the input clusters, usually GeV)"
+        arg_type = Float64
+        default = 5.0
+
+        "--exclusive-dcut"
+        help = "Return all exclusive jets where further merging would have d>d_cut"
+        arg_type = Float64
+
+        "--exclusive-njets"
+        help = "Return all exclusive jets once clusterisation has produced n jets"
+        arg_type = Int
+
+        "--distance", "-R"
+        help = "Distance parameter for jet merging"
+        arg_type = Float64
+        default = 0.4
+
+        "--algorithm", "-A"
+        help = """Algorithm to use for jet reconstruction: $(join(JetReconstruction.AllJetRecoAlgorithms, ", "))"""
+        arg_type = JetAlgorithm.Algorithm
+        default = JetAlgorithm.AntiKt
+
+        "--strategy", "-S"
+        help = """Strategy for the algorithm, valid values: $(join(JetReconstruction.AllJetRecoStrategies, ", "))"""
+        arg_type = RecoStrategy.Strategy
+        default = RecoStrategy.Best
+
+        "file"
+        help = "HepMC3 event file in HepMC3 to read"
+        required = true
+
+        "output"
+        help = "File for output image"
+        default = "jetvis.png"
+    end
+    args = parse_args(ARGS, s; as_symbols = true)
+
+
+    logger = ConsoleLogger(stdout, Logging.Info)
+    global_logger(logger)
+
+    events::Vector{Vector{PseudoJet}} =
+        read_final_state_particles(args[:file], maxevents = args[:event], skipevents = args[:event])
+
+    power = JetReconstruction.algorithm2power[args[:algorithm]]
+    cs = jet_reconstruct(events[1], R = args[:distance], p = power, strategy = args[:strategy])
+
+    plt = jetsplot(events[1], cs; Module = CairoMakie)
+    save(args[:output], plt)
+
+    @info "Saved jet visualisation to $(args[:output])"
+end
+
+main()
diff --git a/src/JetVis.jl b/ext/JetVisualisation.jl
similarity index 62%
rename from src/JetVis.jl
rename to ext/JetVisualisation.jl
index bcaa20a..f5eb57b 100644
--- a/src/JetVis.jl
+++ b/ext/JetVisualisation.jl
@@ -1,17 +1,19 @@
 ## Jet visualisation
-# not a submodule
+
+module JetVisualisation
+
+using JetReconstruction
+using Makie
 
 function get_all_ancestors(idx, cs::ClusterSequence)
-	if cs.history[idx].parent1 == NonexistentParent
-		return [cs.history[idx].jetp_index]
-	#elseif cs.history[idx].parent2 == BeamJet
-	#	return 
-	else
-		branch1 = get_all_ancestors(cs.history[idx].parent1, cs)
-		cs.history[idx].parent2 == BeamJet && return branch1
-		branch2 = get_all_ancestors(cs.history[idx].parent2, cs)
-		return [branch1; branch2]
-	end
+    if cs.history[idx].parent1 == JetReconstruction.NonexistentParent
+        return [cs.history[idx].jetp_index]
+    else
+        branch1 = get_all_ancestors(cs.history[idx].parent1, cs)
+        cs.history[idx].parent2 == JetReconstruction.BeamJet && return branch1
+        branch2 = get_all_ancestors(cs.history[idx].parent2, cs)
+        return [branch1; branch2]
+    end
 end
 
 """
@@ -42,14 +44,14 @@ using WGLMakie
 jetsplot(my_objects, cs, Module=Main) #default
 ```
 """
-function jetsplot(objects, cs::ClusterSequence; barsize_phi=0.1, barsize_eta=0.1, colormap=:glasbey_hv_n256, Module=Main)
-	idx_arrays = Vector{Int}[]
+function JetReconstruction.jetsplot(objects, cs::ClusterSequence; barsize_phi = 0.1, barsize_eta = 0.1, colormap = :glasbey_hv_n256, Module = CairoMakie)
+    idx_arrays = Vector{Int}[]
     for elt in cs.history
-        elt.parent2 == BeamJet || continue
-		push!(idx_arrays, get_all_ancestors(elt.parent1, cs))
+        elt.parent2 == JetReconstruction.BeamJet || continue
+        push!(idx_arrays, get_all_ancestors(elt.parent1, cs))
     end
 
-	jetsplot(objects, idx_arrays; barsize_phi, barsize_eta, colormap, Module)
+    jetsplot(objects, idx_arrays; barsize_phi, barsize_eta, colormap, Module)
 end
 
 """
@@ -82,26 +84,28 @@ using WGLMakie
 jetsplot(my_objects, my_colour_arrays, Module=Main) #default
 ```
 """
-function jetsplot(objects, idx_arrays; barsize_phi=0.1, barsize_eta=0.1, colormap=:glasbey_hv_n256, Module=Main)
-	cs = fill(0, length(objects)) # colours
-	for i in 1:length(idx_arrays), j in idx_arrays[i]
-		cs[j] = i
-	end
-
-	pts = sqrt.(pt2.(objects))
-
-	Module.meshscatter(
-		Module.Point3f.(phi.(objects), rapidity.(objects), 0pts);
-	  	color = cs,
-		markersize = Module.Vec3f.(barsize_phi, barsize_eta, pts),
-		colormap = colormap,
-		marker = Module.Rect3f(Module.Vec3f(0), Module.Vec3f(1)),
-	 	figure = (resolution=(700,600),),
-		axis = (
-			type = Module.Axis3, perspectiveness = 0.5, azimuth = 2.6, elevation=0.5,
-                        xlabel = "ϕ", ylabel = "η", zlabel = "kt",
-		        limits = (nothing, nothing, nothing, nothing, 0, findmax(pts)[1]+10)
-		),
-	    shading=false
-	)
+function JetReconstruction.jetsplot(objects, idx_arrays; barsize_phi = 0.1, barsize_eta = 0.1, colormap = :glasbey_hv_n256, Module = Main)
+    cs = fill(0, length(objects)) # colours
+    for i in 1:length(idx_arrays), j in idx_arrays[i]
+        cs[j] = i
+    end
+
+    pts = sqrt.(JetReconstruction.pt2.(objects))
+
+    Module.meshscatter(
+        Module.Point3f.(JetReconstruction.phi.(objects), JetReconstruction.rapidity.(objects), 0pts);
+        color = cs,
+        markersize = Module.Vec3f.(barsize_phi, barsize_eta, pts),
+        colormap = colormap,
+        marker = Module.Rect3f(Module.Vec3f(0), Module.Vec3f(1)),
+        figure = (size = (700, 600),),
+        axis = (
+            type = Module.Axis3, perspectiveness = 0.5, azimuth = 2.6, elevation = 0.5,
+            xlabel = "ϕ", ylabel = "η", zlabel = "kt",
+            limits = (nothing, nothing, nothing, nothing, 0, findmax(pts)[1] + 10),
+        ),
+        shading = NoShading,
+    )
+end
+
 end
diff --git a/img/jetvis.png b/img/jetvis.png
new file mode 100644
index 0000000..9429a97
Binary files /dev/null and b/img/jetvis.png differ
diff --git a/src/AlgorithmStrategyEnums.jl b/src/AlgorithmStrategyEnums.jl
index 1bd9365..ed24065 100644
--- a/src/AlgorithmStrategyEnums.jl
+++ b/src/AlgorithmStrategyEnums.jl
@@ -6,13 +6,16 @@ using EnumX
 const AllJetRecoStrategies = [String(Symbol(x)) for x in instances(RecoStrategy.Strategy)]
 
 # Algorithm emun
-@enumx T = Algorithm JetAlgorithm AntiKt Cambridge Kt EEKt Durham
+@enumx T = Algorithm JetAlgorithm AntiKt CA Kt EEKt Durham
 const AllJetRecoAlgorithms = [String(Symbol(x)) for x in instances(JetAlgorithm.Algorithm)]
 
-# Map from power values to algorithms
+# Map from algorithms to power values used 
 power2algorithm = Dict(-1 => JetAlgorithm.AntiKt,
-    0 => JetAlgorithm.Cambridge,
+    0 => JetAlgorithm.CA,
     1 => JetAlgorithm.Kt)
+algorithm2power = Dict(JetAlgorithm.AntiKt => -1,
+    JetAlgorithm.CA => 0,
+    JetAlgorithm.Kt => 1)
 
 # Map from string to an enum value (used for CLI parsing)
 Base.tryparse(E::Type{<:Enum}, str::String) =
diff --git a/src/ClusterSequence.jl b/src/ClusterSequence.jl
index 7037881..adf98c4 100644
--- a/src/ClusterSequence.jl
+++ b/src/ClusterSequence.jl
@@ -170,7 +170,7 @@ function exclusive_jets(clusterseq::ClusterSequence; dcut = nothing, njets = not
     end
 
     # Check that an algorithm was used that makes sense for exclusive jets
-    if !(clusterseq.algorithm ∈ (JetAlgorithm.Cambridge, JetAlgorithm.Kt, JetAlgorithm.EEKt, JetAlgorithm.Durham))
+    if !(clusterseq.algorithm ∈ (JetAlgorithm.CA, JetAlgorithm.Kt, JetAlgorithm.EEKt, JetAlgorithm.Durham))
         throw(ArgumentError("Algorithm used is not suitable for exclusive jets ($(clusterseq.algorithm))"))
     end
 
@@ -189,10 +189,10 @@ function exclusive_jets(clusterseq::ClusterSequence; dcut = nothing, njets = not
 
     excl_jets = LorentzVectorCyl[]
     for j in stop_point:length(clusterseq.history)
-        @info "Search $j ($(clusterseq.history[j].parent1) + $(clusterseq.history[j].parent2))"
+        @debug "Search $j ($(clusterseq.history[j].parent1) + $(clusterseq.history[j].parent2))"
         for parent in (clusterseq.history[j].parent1, clusterseq.history[j].parent2)
             if (parent < stop_point && parent > 0)
-                @info "Added exclusive jet index $(clusterseq.history[parent].jetp_index)"
+                @debug "Added exclusive jet index $(clusterseq.history[parent].jetp_index)"
                 jet = clusterseq.jets[clusterseq.history[parent].jetp_index]
                 push!(excl_jets, LorentzVectorCyl(pt(jet), rapidity(jet), phi(jet), mass(jet)))
             end
@@ -206,7 +206,7 @@ end
 """Return all number of exclusive jets of a ClusterSequence that are above a certain dcut value"""
 function n_exclusive_jets(clusterseq::ClusterSequence; dcut::AbstractFloat)
     # Check that an algorithm was used that makes sense for exclusive jets
-    if !(clusterseq.algorithm ∈ (JetAlgorithm.Cambridge, JetAlgorithm.Kt, JetAlgorithm.EEKt, JetAlgorithm.Durham))
+    if !(clusterseq.algorithm ∈ (JetAlgorithm.CA, JetAlgorithm.Kt, JetAlgorithm.EEKt, JetAlgorithm.Durham))
         throw(ArgumentError("Algorithm used is not suitable for exclusive jets ($(clusterseq.algorithm))"))
     end
 
@@ -214,7 +214,7 @@ function n_exclusive_jets(clusterseq::ClusterSequence; dcut::AbstractFloat)
     # first time max_dij_so_far > dcut)
     i_dcut = length(clusterseq.history)
     for i_history ∈ length(clusterseq.history):-1:1
-        @info "Examining $i_history, max_dij=$(clusterseq.history[i_history].max_dij_so_far)"
+        @debug "Examining $i_history, max_dij=$(clusterseq.history[i_history].max_dij_so_far)"
         if clusterseq.history[i_history].max_dij_so_far <= dcut
             i_dcut = i_history
             break
@@ -224,15 +224,4 @@ function n_exclusive_jets(clusterseq::ClusterSequence; dcut::AbstractFloat)
     # The number of jets is then given by this formula
     length(clusterseq.history) - i_dcut
 
-    # int i = _history.size() - 1; // last jet
-    # while (i >= 0) {
-    #   if (_history[i].max_dij_so_far <= dcut) {break;}
-    #   i--;
-    # }
-    # int stop_point = i + 1;
-    # // relation between stop_point, njets assumes one extra jet disappears
-    # // at each clustering.
-    # int njets = 2*_initial_n - stop_point;
-    # return njets;
-
 end
\ No newline at end of file
diff --git a/src/JSONresults.jl b/src/JSONresults.jl
index 8126ab6..e52788c 100644
--- a/src/JSONresults.jl
+++ b/src/JSONresults.jl
@@ -10,4 +10,3 @@ struct FinalJets
     jetid::Int64
     jets::Vector{FinalJet}
 end
-
diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl
index fd36ab7..bb1af4e 100644
--- a/src/JetReconstruction.jl
+++ b/src/JetReconstruction.jl
@@ -61,12 +61,12 @@ export savejets, loadjets!, loadjets
 include("Utils.jl")
 export read_final_state_particles, read_final_state_particles_lv, pseudojets2vectors, final_jets
 
-# jet visualisation
-include("JetVis.jl")
+# Jet visualisation is an extension, see ext/JetVisualisation.jl
+function jetsplot() end
 export jetsplot
 
 # JSON results
 include("JSONresults.jl")
-export FinalJet, FinalJets, JSON3
+export FinalJet, FinalJets
 
 end
diff --git a/src/Pseudojet.jl b/src/Pseudojet.jl
index f37b14f..1d3a739 100644
--- a/src/Pseudojet.jl
+++ b/src/Pseudojet.jl
@@ -24,7 +24,7 @@ const _invalid_rap = -1.e200
 # \class PseudoJet
 # Class to contain pseudojets, including minimal information of use to
 # jet-clustering routines.
-mutable struct PseudoJet<:FourMomentum
+mutable struct PseudoJet <: FourMomentum
     # construct a pseudojet from explicit components
     px::Float64
     py::Float64
@@ -37,20 +37,20 @@ mutable struct PseudoJet<:FourMomentum
     _phi::Float64
 end
 
-    
+
 PseudoJet(px::Float64, py::Float64, pz::Float64, E::Float64,
-          _cluster_hist_index::Int,
-          pt2::Float64) = PseudoJet(px,
-                                    py, pz, E, _cluster_hist_index,
-                                    pt2, 1. / pt2, _invalid_rap, _invalid_phi)
+    _cluster_hist_index::Int,
+    pt2::Float64) = PseudoJet(px,
+    py, pz, E, _cluster_hist_index,
+    pt2, 1.0 / pt2, _invalid_rap, _invalid_phi)
 
 PseudoJet(px::Float64, py::Float64,
-          pz::Float64, E::Float64) = PseudoJet(px, py, pz, E, 0, px^2 + py^2)
+    pz::Float64, E::Float64) = PseudoJet(px, py, pz, E, 0, px^2 + py^2)
 
 import Base.show
 show(io::IO, jet::PseudoJet) = begin
     print(io, "Pseudojet(px: ", jet.px, " py: ", jet.py, " pz: ", jet.pz, " E: ", jet.E, "; ",
-          "pt: ", sqrt(jet._pt2), " eta: ", rapidity(jet), " phi: ", phi(jet), ", m: ", m(jet), ")")
+        "pt: ", sqrt(jet._pt2), " eta: ", rapidity(jet), " phi: ", phi(jet), ", m: ", m(jet), ")")
 end
 
 
@@ -60,7 +60,7 @@ set_momentum!(j::PseudoJet, px, py, pz, E) = begin
     j.pz = pz
     j.E = E
     j._pt2 = px^2 + py^2
-    j._inv_pt2 = 1.0/j._pt2
+    j._inv_pt2 = 1.0 / j._pt2
     j._rap = _invalid_eta
     j._phi = _invalid_phi
 end
@@ -69,10 +69,10 @@ _ensure_valid_rap_phi(p::PseudoJet) = p._phi == _invalid_phi && _set_rap_phi!(p)
 
 _set_rap_phi!(p::PseudoJet) = begin
 
-    p._phi = p._pt2 == 0.0 ? 0.0 : atan(p.py,p.px)
+    p._phi = p._pt2 == 0.0 ? 0.0 : atan(p.py, p.px)
     if p._phi < 0.0
-         p._phi += 2π
-     elseif p._phi >= 2π
+        p._phi += 2π
+    elseif p._phi >= 2π
         p._phi -= 2π  # can happen if phi=-|eps<1e-15|?
     end
 
@@ -82,15 +82,17 @@ _set_rap_phi!(p::PseudoJet) = begin
         #    different rapidities (so as to lift the degeneracy between
         #                         them) [this can be relevant at parton-level]
         MaxRapHere = _MaxRap + abs(p.pz)
-        p._rap = p.pz >= 0.0 ?  MaxRapHere : -MaxRapHere
+        p._rap = p.pz >= 0.0 ? MaxRapHere : -MaxRapHere
     else
         # get the rapidity in a way that's modestly insensitive to roundoff
         # error when things pz,E are large (actually the best we can do without
         # explicit knowledge of mass)
         effective_m2 = max(0.0, m2(p)) # force non tachyonic mass
         E_plus_pz    = p.E + abs(p.pz) # the safer of p+, p-
-        p._rap = 0.5*log((p._pt2 + effective_m2)/(E_plus_pz*E_plus_pz))
-        if p.pz > 0 p._rap = - p._rap; end
+        p._rap       = 0.5 * log((p._pt2 + effective_m2) / (E_plus_pz * E_plus_pz))
+        if p.pz > 0
+            p._rap = -p._rap
+        end
     end
     nothing
 end
@@ -113,7 +115,7 @@ pt2(p::PseudoJet) = p._pt2
 pt(p::PseudoJet) = sqrt(p._pt2)
 
 "Returns the squared invariant mass"
-m2(p::PseudoJet) = (p.E + p.pz)*(p.E-p.pz) - p._pt2
+m2(p::PseudoJet) = (p.E + p.pz) * (p.E - p.pz) - p._pt2
 
 "Returns the magnitude of the momentum, |p|"
 mag(p::PseudoJet) = sqrt(muladd(p.px, p.px, p.py^2) + p.pz^2)
@@ -140,7 +142,7 @@ const η = eta
 # (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP, SciKitHEP Particle and ROOT)
 m(p::PseudoJet) = begin
     x = m2(p)
-    x < 0. ? -sqrt(-x) : sqrt(x)
+    x < 0.0 ? -sqrt(-x) : sqrt(x)
 end
 
 # Ensure we have accessors for jet parameters
@@ -155,5 +157,5 @@ import Base.+;
 
 +(j1::PseudoJet, j2::PseudoJet) = begin
     PseudoJet(j1.px + j2.px, j1.py + j2.py,
-                                            j1.pz + j2.pz, j1.E + j2.E)
+        j1.pz + j2.pz, j1.E + j2.E)
 end