[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,129 @@
+# Jet Substructure
+
+## Structures
+
+### `MassDropTagger`
+
+The `MassDropTagger` structure is used for tagging jets based on mass drop conditions, which helps in identifying subjets within a jet that undergo a significant drop in mass.
+
+**Fields**:
+
+- `mu::Float64`: Maximum allowed mass ratio for a jet to pass the tagging condition.
+- `y::Float64`: Minimum kT distance threshold for parent jet separation.
+
+---
+
+### `SoftDropTagger`
+
+The `SoftDropTagger` instance is used to apply soft-drop grooming to jets, removing soft, wide-angle radiation. This approach is commonly used in jet grooming to reduce contamination from soft particles.
+
+**Fields**:
+
+- `zcut::Float64`: Minimum allowed energy fraction for subjets.
+- `b::Float64`: Angular exponent controlling soft radiation suppression.
+- `cluster_rad::Float64`: New radius used to recluster components of the jet. Defaults to `1.0` if no value is specified.
+
+---
+
+### `JetFilter`
+
+The `JetFilter` structure is used to filter jets based on a specific radius and the number of hardest subjets. This technique reduces contamination from peripheral soft particles.
+
+**Fields**:
+
+- `filter_radius::Float64`: Radius parameter used to recluster subjets.
+- `num_hardest_jets::Int`: Number of hardest subjets retained in the filtered result.
+
+---
+
+### `JetTrim`
+
+`JetTrim` instance is used to trim jets by removing soft, large-angle components from the jet. This is useful in cleaning up jets to remove softer particles at wide angles.
+
+**Fields**:
+
+- `trim_radius::Float64`: Radius used for reclustering in trimming.
+- `trim_fraction::Float64`: Minimum momentum fraction for retained subjets.
+- `recluster_method::JetAlgorithm.Algorithm`: Method identifier for reclustering.
+
+---
+
+## Functions
+
+### `mass_drop`
+
+```julia
+mass_drop(jet::PseudoJet, clusterseq::ClusterSequence, tag::MassDropTagger) -> PseudoJet
+```
+
+The `mass_drop` function identifies subjets in a jet that pass the mass drop tagging condition. It iterates through the clustering history of the jet, stopping at the first jet that satisfies the mass and distance thresholds.
+
+**Arguments** :
+
+* `jet`: `PseudoJet` instance representing the jet to be tagged.
+* `clusterseq`: `ClusterSequence` with jet clustering history.
+* `tag`: `MassDropTagger` instance providing mass drop parameters.
+
+**Returns** :
+
+`PseudoJet`: The jet (or subjet) that satisfies the mass drop condition, or a zero-momentum `PseudoJet` if no tagging occurs.
+
+---
+
+### `soft_drop`
+
+```julia
+soft_drop(jet::PseudoJet, clusterseq::ClusterSequence, tag::SoftDropTagger) -> PseudoJet
+```
+
+The `soft_drop` 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.
+
+**Arguments** :
+
+* `jet`: `PseudoJet` instance to groom.
+* `clusterseq`: `ClusterSequence` containing jet history.
+* `tag`: `SoftDropTagger` instance with soft-drop parameters.
+
+**Returns** :
+
+`PseudoJet`: Groomed jet or zero-momentum `PseudoJet` if grooming fails.
+
+---
+
+### `jet_filtering`
+
+```julia-repl
+jet_filtering(jet::PseudoJet, clusterseq::ClusterSequence, filter::JetFilter) -> PseudoJet
+```
+
+The `jet_filtering` function filters a jet to retain only the hardest subjets based on a specified radius and number. This helps in refining the jet structure by reducing soft particle contamination.
+
+**Arguments** :
+
+* `jet`: `PseudoJet` instance representing the jet to filter.
+* `clusterseq`: `ClusterSequence` containing jet history.
+* `filter`: `JetFilter` instance specifying radius and number of subjets.
+
+**Returns** :
+
+`PseudoJet`: Filtered jet composed of the hardest subjets.
+
+---
+
+### `jet_trimming`
+
+```julia
+jet_trimming(jet::PseudoJet, clusterseq::ClusterSequence, trim::JetTrim) -> PseudoJet
+```
+
+The `jet_trimming` 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.
+
+**Arguments** :
+
+* `jet`: `PseudoJet` instance representing the jet to trim.
+* `clusterseq`: `ClusterSequence` containing jet history.
+* `trim`: `JetTrim` instance specifying trimming parameters.
+
+**Returns** :
+
+`PseudoJet`: Trimmed jet composed of retained subjets.

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

@@ -69,6 +69,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")