[WIP] Implementation of substructure modules

docs/make.jl

@@ -17,6 +17,7 @@ makedocs(sitename = "JetReconstruction.jl",
              "Visualisation" => "visualisation.md",
              "Particle Inputs" => "particles.md",
              "Reconstruction Strategies" => "strategy.md",
+             "Substructure" => "substructure.md",
              "Reference Docs" => Any["Public API" => "lib/public.md",
                                      "Internal API" => "lib/internal.md"],
              "Extras" => Any["Serialisation" => "extras/serialisation.md"]

docs/src/substructure.md

@@ -0,0 +1,88 @@
+# Jet Substructure
+Jet substructure techniques provide powerful tools for analyzing and refining the properties of jets. Below are some of the key jet substructure functions, that are available.
+
+### Mass Drop Tagging
+
+```julia
+mass_drop(jet::PseudoJet, clusterseq::ClusterSequence, tag::MassDropTagger) -> PseudoJet
+```
+
+The [mass_drop](@ref) function identifies subjets in a jet that pass the mass drop tagging condition. To use the `mass_drop` function:
+
+-  Configure the [MassDropTagger](@ref) with parameters tailored to the expected mass drop and distance thresholds for analysis.
+``` julia
+tagger = MassDropTagger(mu=0.67, y=0.09)
+```
+- Apply the `mass_drop` function to the jet and its clustering sequence.
+``` julia 
+tagged_jet = mass_drop(jet, clusterseq, tagger)
+```
+- If the jet is tagged successfully, the function returns the identified subjet. Else it returns `PseudoJet(0.0, 0.0, 0.0, 0.0)`.
+---
+
+### Soft Drop Tagging
+
+```julia
+soft_drop(jet::PseudoJet, clusterseq::ClusterSequence, tag::SoftDropTagger) -> PseudoJet
+```
+
+The [soft_drop](@ref) function applies soft-drop grooming to remove soft, wide-angle radiation from jets. It reclusters the jet with a specified radius and clustering method, iteratively checking the soft-drop condition on subjets. To use the `soft_drop` function:
+
+- Construct the [SoftDropTagger](@ref) by defining the parameters for the grooming process, such as the energy fraction (`zcut`) and angular exponent (`b`)
+``` julia
+tagger = SoftDropTagger(zcut=0.1, b=2.0)
+```
+By default, the reclustering radius is set to `1.0` which can be modified by defing the `tagger` object as
+``` julia
+tagger = SoftDropTagger(zcut=0.1, b=2.0, cluster_rad=0.4)
+```
+- Apply the `soft_drop` function to the jet and its clustering sequence.
+``` julia 
+tagged_jet = soft_drop(jet, clusterseq, tagger)
+```
+- If the jet is tagged successfully, the function returns the identified subjet. Else it returns `nothing`.
+
+---
+
+### Jet Filtering
+
+```julia
+jet_filtering(jet::PseudoJet, clusterseq::ClusterSequence, filter::JetFilter) -> PseudoJet
+```
+
+The [jet_filtering](@ref) function filters a jet to retain only the hardest subjets based on a specified radius and number. To use the function:
+
+- Set the radius for subjet clustering and the number of hardest subjets to retain in a [JetFilter](@ref) method.
+```julia
+filter = JetFilter(filter_radius=0.3, num_hardest_jets=3)
+```
+- Apply the `jet_filtering` function to refine the jet.
+```julia
+filtered_jet = jet_filtering(jet, clusterseq, filter)
+```
+- The function returns the filtered jet that retains only the most significant subjets, reducing noise.
+---
+
+### Jet Trimming
+
+```julia
+jet_trimming(jet::PseudoJet, clusterseq::ClusterSequence, trim::JetTrim) -> PseudoJet
+```
+
+The [jet_trimming](@ref) function trims a jet by removing subjets with transverse momentum below a specified fraction of the main jet's momentum. This method cleans up jets by removing soft particles. To use this function:
+
+- Configure the trimming radius and momentum fraction threshold in [JetTrim](@ref)
+```julia
+trim = JetTrim(trim_radius=0.3, trim_fraction=0.3, recluster_method=JetAlgorithm.CA)
+```
+It is to be noted that the `jet_trimming` function reclusters the constituents of the jet using either `C/A` or `kT
+` algorithm, which needs to be specified.
+- Apply the `jet_trimming` function to clean the jet.
+```julia
+trimmed_jet = jet_trimming(jet, clusterseq, trim)
+```
+- The function returns the trimmed jet. 
+
+---
+
+For a more detailed guide to use these substructure modules, one can refer to the provided examples in the `examples/Substructure` directory.

examples/Substructure/JetGrooming.jl

@@ -0,0 +1,39 @@
+#! /usr/bin/env julia
+using JetReconstruction
+
+input_file = joinpath(dirname(pathof(JetReconstruction)),
+                      "..", "test", "data", "events.pp13TeV.hepmc3.gz")
+events = read_final_state_particles(input_file)
+
+# Event to pick
+event_no = 1
+
+cluster_seq = jet_reconstruct(events[event_no], p = 0, R = 1.0)
+jets = inclusive_jets(cluster_seq; ptmin = 5.0, T = PseudoJet)
+
+r = 0.3     # recluster radius
+n = 3       # number of hard jets to consider
+
+filter = JetFilter(r, n)
+
+@info "Jet Filtering: recluster radius = $r, hard subjets to consider = $n"
+for jet in jets
+    filtered = jet_filtering(jet, cluster_seq, filter)
+
+    println("Original jet: pt = $(JetReconstruction.pt(jet)), rap = $(JetReconstruction.rapidity(jet)), phi = $(JetReconstruction.phi(jet)), E = $(jet.E)")
+    println("Filtered jet: pt = $(JetReconstruction.pt(filtered)), rap = $(JetReconstruction.rapidity(filtered)), phi = $(JetReconstruction.phi(filtered)), E = $(filtered.E)\n")
+end
+
+r = 0.3                     # recluster radius
+f = 0.3                     # trim fraction
+m = JetAlgorithm.CA         # recluster method
+
+trim = JetTrim(r, f, m)
+
+@info "Jet Trimming: recluster radius = $r, trim fraction = $f, recluster method = $m"
+for jet in jets
+    trimmed = jet_trimming(jet, cluster_seq, trim)
+
+    println("Original jet: pt = $(JetReconstruction.pt(jet)), rap = $(JetReconstruction.rapidity(jet)), phi = $(JetReconstruction.phi(jet)), E = $(jet.E)")
+    println("Trimmed jet: pt = $(JetReconstruction.pt(trimmed)), rap = $(JetReconstruction.rapidity(trimmed)), phi = $(JetReconstruction.phi(trimmed)), E = $(trimmed.E)\n")
+end

examples/Substructure/JetTagging.jl

@@ -0,0 +1,36 @@
+#! /usr/bin/env julia
+using JetReconstruction
+
+input_file = joinpath(dirname(pathof(JetReconstruction)),
+                      "..", "test", "data", "events.pp13TeV.hepmc3.gz")
+events = read_final_state_particles(input_file)
+
+# Event to pick
+event_no = 1
+
+cluster_seq = jet_reconstruct(events[event_no], p = 0, R = 1.0)
+jets = inclusive_jets(cluster_seq; ptmin = 5.0, T = PseudoJet)
+
+μ = 0.67        # jet mass ratio
+y = 0.09        # symmetry cut
+
+MDtagger = MassDropTagger(μ, y)
+
+@info "Mass Drop Tagging: μ = $μ, y = $y"
+for jet in jets
+    tagged = mass_drop(jet, cluster_seq, MDtagger)
+    println("Original jet: pt = $(JetReconstruction.pt(jet)), rap = $(JetReconstruction.rapidity(jet)), phi = $(JetReconstruction.phi(jet)), E = $(jet.E)")
+    println("Tagged jet: pt = $(JetReconstruction.pt(tagged)), rap = $(JetReconstruction.rapidity(tagged)), phi = $(JetReconstruction.phi(tagged)), E = $(tagged.E)\n")
+end
+
+z = 0.1         # soft drop threshold
+b = 2.0         # angular exponent
+
+SDtagger = SoftDropTagger(z, b)
+
+@info "Soft Drop Tagging: recluster radius = $(SDtagger.cluster_rad), zcut = $z, b = $b"
+for jet in jets
+    tagged = soft_drop(jet, cluster_seq, SDtagger)
+    println("Original jet: pt = $(JetReconstruction.pt(jet)), rap = $(JetReconstruction.rapidity(jet)), phi = $(JetReconstruction.phi(jet)), E = $(jet.E)")
+    println("Tagged jet: pt = $(JetReconstruction.pt(tagged)), rap = $(JetReconstruction.rapidity(tagged)), phi = $(JetReconstruction.phi(tagged)), E = $(tagged.E)\n")
+end

examples/Substructure/Project.toml

@@ -0,0 +1,2 @@
+[deps]
+JetReconstruction = "44e8cb2c-dfab-4825-9c70-d4808a591196"

examples/Substructure/README.md

@@ -0,0 +1,14 @@
+# Examples for Substructure Modules
+
+The `JetGrooming.jl` file shows the usage of `jet_filtering` and `jet_trimming` functions while the `JetTagging.jl` file demonstrates how to use `mass_drop` and `soft_drop` functions. 
+
+To use these examples run
+
+```julia
+julia --project JetTagging.jl
+...
+julia --project JetGrooming.jl
+...
+```
+
+The parameters of tagging and grooming and the input files can be easily changed in the scripts.

src/JetReconstruction.jl

@@ -70,6 +70,11 @@ export ee_genkt_algorithm
 include("GenericAlgo.jl")
 export jet_reconstruct
 
+## Substructure modules
+include("Substructure.jl")
+export MassDropTagger, SoftDropTagger, JetFilter, JetTrim, mass_drop, soft_drop,
+       jet_filtering, jet_trimming
+
 # Simple HepMC3 reader
 include("HepMC3.jl")