Skip to content

Commit de57d38

Browse files
committed
Added PseudoJet example in docs
1 parent 121d178 commit de57d38

File tree

1 file changed

+158
-0
lines changed

1 file changed

+158
-0
lines changed

docs/src/examples.md

Lines changed: 158 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -291,9 +291,167 @@ event = build_event_template()
291291
println("Event built: $(particles_size(event)) particles")
292292
```
293293

294+
## JetReconstruction Integration
295+
296+
HepMC3.jl can be used with [JetReconstruction.jl](https://github.com/JuliaHEP/JetReconstruction.jl) by converting HepMC3 particles to `PseudoJet` objects. The `PseudoJet` struct is a generic data structure that can be used for both proton-proton (pp) and electron-positron (e⁺e⁻) events. The difference lies in the jet reconstruction algorithm used:
297+
298+
- **For pp events**: Use algorithms like `JetAlgorithm.AntiKt`, `JetAlgorithm.Kt`, or `JetAlgorithm.CA` (Cambridge/Aachen)
299+
- **For e⁺e⁻ events**: Use algorithms like `JetAlgorithm.Durham` or `JetAlgorithm.EEKt` (generalized e⁺e⁻ kₜ)
300+
301+
This example demonstrates reading HepMC3 files and converting final state particles to the `PseudoJet` format:
302+
303+
```julia
304+
using HepMC3, JetReconstruction
305+
306+
"""
307+
read_final_state_particles(filename; max_events=-1)
308+
Read HepMC3 file and return JetReconstruction-compatible PseudoJet events.
309+
This function extracts final state particles (status == 1) from HepMC3 events
310+
and converts them to PseudoJet format for jet reconstruction.
311+
"""
312+
function read_final_state_particles(filename::String; max_events::Int=-1)
313+
# Read events using HepMC3 interface (handles compressed files automatically)
314+
events = read_hepmc_file_with_compression(filename; max_events=max_events)
315+
316+
# Convert to PseudoJets for JetReconstruction
317+
pseudojet_events = Vector{PseudoJet}[]
318+
319+
events_processed = 0
320+
for (event_idx, event_ptr) in enumerate(events)
321+
# Extract final state particles (status == 1)
322+
final_state = get_final_state_particles(event_ptr)
323+
324+
# Skip empty events
325+
if isempty(final_state)
326+
continue
327+
end
328+
329+
# Convert HepMC3 particles to PseudoJets
330+
input_particles = PseudoJet[]
331+
particle_index = 1
332+
333+
for particle in final_state
334+
props = get_particle_properties(particle)
335+
336+
# Create PseudoJet from four-momentum
337+
pseudojet = PseudoJet(
338+
props.momentum.px, # px
339+
props.momentum.py, # py
340+
props.momentum.pz, # pz
341+
props.momentum.e; # E
342+
cluster_hist_index = particle_index
343+
)
344+
push!(input_particles, pseudojet)
345+
particle_index += 1
346+
end
347+
348+
push!(pseudojet_events, input_particles)
349+
events_processed += 1
350+
351+
# Respect max_events limit
352+
if max_events > 0 && events_processed >= max_events
353+
break
354+
end
355+
end
356+
357+
@info "Processed $(length(pseudojet_events)) events from $filename"
358+
return pseudojet_events
359+
end
360+
361+
# Example: Read events and perform jet reconstruction
362+
filename = "events.hepmc3.zst" # Can be compressed or uncompressed
363+
364+
# Read final state particles
365+
pseudojet_events = read_final_state_particles(filename; max_events=100)
366+
367+
# Perform jet reconstruction on first event
368+
if !isempty(pseudojet_events)
369+
event_particles = pseudojet_events[1]
370+
371+
# For pp events: use anti-kT algorithm
372+
clusterseq_pp = jet_reconstruct(event_particles;
373+
algorithm = JetAlgorithm.AntiKt,
374+
R = 0.4)
375+
jets_pp = inclusive_jets(clusterseq_pp, ptmin = 5.0)
376+
377+
# For e⁺e⁻ events: use Durham algorithm (R parameter is ignored)
378+
clusterseq_ee = jet_reconstruct(event_particles;
379+
algorithm = JetAlgorithm.Durham)
380+
jets_ee = inclusive_jets(clusterseq_ee, ptmin = 5.0)
381+
382+
println("Event 1: $(length(event_particles)) particles")
383+
println("PP jets (anti-kT, R=0.4): $(length(jets_pp)) jets with pT > 5 GeV")
384+
println("e⁺e⁻ jets (Durham): $(length(jets_ee)) jets with pT > 5 GeV")
385+
end
386+
```
387+
388+
### Complete Pipeline Example
389+
390+
A complete workflow from HepMC3 file to jet analysis. Note that the same `PseudoJet` conversion works for both pp and e⁺e⁻ events; only the algorithm differs:
391+
392+
```julia
393+
using HepMC3, JetReconstruction
394+
395+
# For pp events
396+
function analyze_pp_jets_from_hepmc(filename::String; max_events::Int=10, R::Float64=0.4, ptmin::Float64=5.0)
397+
pseudojet_events = read_final_state_particles(filename; max_events=max_events)
398+
399+
println("Analyzing $(length(pseudojet_events)) pp events")
400+
println("=" ^ 60)
401+
402+
for (event_idx, event_particles) in enumerate(pseudojet_events)
403+
# Use anti-kT for pp events
404+
clusterseq = jet_reconstruct(event_particles;
405+
algorithm = JetAlgorithm.AntiKt,
406+
R = R)
407+
jets = inclusive_jets(clusterseq, ptmin = ptmin)
408+
409+
println("Event $event_idx: $(length(jets)) jets (anti-kT, R=$R)")
410+
end
411+
end
412+
413+
# For e⁺e⁻ events
414+
function analyze_ee_jets_from_hepmc(filename::String; max_events::Int=10, ptmin::Float64=5.0)
415+
pseudojet_events = read_final_state_particles(filename; max_events=max_events)
416+
417+
println("Analyzing $(length(pseudojet_events)) e⁺e⁻ events")
418+
println("=" ^ 60)
419+
420+
for (event_idx, event_particles) in enumerate(pseudojet_events)
421+
# Use Durham algorithm for e⁺e⁻ events (R parameter is ignored)
422+
clusterseq = jet_reconstruct(event_particles;
423+
algorithm = JetAlgorithm.Durham)
424+
jets = inclusive_jets(clusterseq, ptmin = ptmin)
425+
426+
println("Event $event_idx: $(length(jets)) jets (Durham)")
427+
end
428+
end
429+
430+
# Run analysis
431+
analyze_pp_jets_from_hepmc("pp_events.hepmc3.zst"; max_events=10, R=0.4, ptmin=5.0)
432+
analyze_ee_jets_from_hepmc("ee_events.hepmc3.zst"; max_events=10, ptmin=5.0)
433+
```
434+
435+
### Key Points
436+
437+
1. **Final State Extraction**: Use `get_final_state_particles()` to extract only stable particles (status == 1) from HepMC3 events.
438+
439+
2. **PseudoJet Conversion**: The `PseudoJet` struct is a generic data structure that works for both pp and e⁺e⁻ events. Convert HepMC3 particle four-momenta to `PseudoJet` objects using the same conversion function regardless of event type.
440+
441+
3. **Algorithm Selection**: The choice of jet reconstruction algorithm depends on the event type:
442+
- **pp events**: Use `JetAlgorithm.AntiKt`, `JetAlgorithm.Kt`, or `JetAlgorithm.CA` with an `R` parameter
443+
- **e⁺e⁻ events**: Use `JetAlgorithm.Durham` or `JetAlgorithm.EEKt` (the `R` parameter is ignored for Durham)
444+
445+
4. **Compressed Files**: The `read_hepmc_file_with_compression()` function automatically handles both compressed (`.zst`) and uncompressed HepMC3 files.
446+
447+
5. **Event Processing**: The conversion function processes events one at a time, allowing for efficient memory usage with large files.
448+
449+
This example demonstrates how to convert HepMC3 event data to the `PseudoJet` format expected by JetReconstruction.jl. The same conversion works for both pp and e⁺e⁻ events; only the jet reconstruction algorithm differs.
450+
294451
## See Also
295452

296453
- [Getting Started](@ref) for installation and basics
297454
- [Events](@ref) for event manipulation
298455
- [Particles](@ref) for particle operations
299456
- [File I/O](@ref) for reading and writing files
457+
- [JetReconstruction.jl](https://github.com/JuliaHEP/JetReconstruction.jl) for jet finding algorithms

0 commit comments

Comments
 (0)