Added FCC example (#5)

* First FCC analysis example

* Added a  Jupyter notebook

* updated documentation

Project.toml

@@ -6,21 +6,22 @@ version = "0.1.0"
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 Corpuscles = "e78a0372-c628-4773-9c8d-eb17d149bf93"
+FHist = "68837c9b-b678-4cd5-9925-8a54edc8f695"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 UnROOT = "3cd96dde-e98d-4713-81e9-a4a1b0235ce9"
 YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
 
 [compat]
-julia = "1"
 Accessors = "0.1"
 Corpuscles = "2"
 Graphs = "1"
 StaticArrays = "1.9"
 StructArrays = "0.6"
-YAML = "0.4"
 UnROOT = "=0.10.22"
+YAML = "0.4"
+julia = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

docs/src/api.md

@@ -5,13 +5,13 @@ Documentation for `EDM4hep.jl` public interface.
 ## Index - Types
 ```@index
 Pages   = ["api.md"]
-Modules = [EDM4hep, EDM4hep.RootIO]
+Modules = [EDM4hep, EDM4hep.RootIO, EDM4hep.Histograms]
 Order   = [:type]
 ```
 ## Index - Functions
 ```@index
 Pages   = ["api.md"]
-Modules = [EDM4hep, EDM4hep.RootIO]
+Modules = [EDM4hep, EDM4hep.RootIO, EDM4hep.Histograms]
 Order   = [:function]
 ```
 
@@ -24,12 +24,12 @@ Order = [:module]
 This is the list of all types defined for EDM4hep using the PODIO yaml file.
 
 ```@autodocs
-Modules = [EDM4hep, EDM4hep.RootIO]
+Modules = [EDM4hep, EDM4hep.RootIO, EDM4hep.Histograms]
 Order = [:type]
 ```
 ## Functions
 ```@autodocs
-Modules = [EDM4hep, EDM4hep.RootIO]
+Modules = [EDM4hep, EDM4hep.RootIO, EDM4hep.Histograms]
 Order = [:function]
 ```
 

docs/src/index.md

@@ -37,6 +37,7 @@ There are a number of issues and problems still to be resolved. We keep track of
 - Generate doc string with member information - *DONE*
 - Generate accessors for one-to-many relations, vector members - *DONE* 
 - Support latest version (RC2) of RNTuple format (waiting for a file being generated)
+- Support for multi-threading (i.e. be able to add `@threads` in the event loop)
 
 
 ## Tests
@@ -77,6 +78,9 @@ for p in mcps
     end
 end
 ```
+### examples/FCC/analysis_mH-recoil.jl
+This is basically the example `higgs/mH-recoil/mumu` from [FCCAnalyses](https://github.com/HEP-FCC/FCCAnalyses). It shows in a realistic manner hoe to develop analysis functions using the Data Types from EDM4hep to created high-level analysis functions. 
+The Jupyter notebook [analysis_mH-recoil.ipynb](./assets/analysis_mH-recoil.pdf) shows the same example in form of a notebook. 
 
 ## EDM4hep Data Model
 This is the diagram for the EDM4hep datamodel including relationships. 

examples/FCC/Project.toml

@@ -0,0 +1,9 @@
+[deps]
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+EDM4hep = "eb32b910-dde9-4347-8fce-cd6be3498f0c"
+FHist = "68837c9b-b678-4cd5-9925-8a54edc8f695"
+IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
+LorentzVectorHEP = "f612022c-142a-473f-8cfd-a09cf3793c6c"
+Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"

examples/FCC/analysis_functions.jl

@@ -0,0 +1,41 @@
+using EDM4hep
+using LorentzVectorHEP
+using Combinatorics
+using IterTools
+
+export resonanceBuilder, recoilBuilder
+
+"""
+    resonanceBuilder(rmass::AbstractFloat, legs::AbstractVector{ReconstructedParticle})
+
+Returns a container with the best resonance of 2 by 2 combinatorics of the `legs` container
+sorted by closest to the input `rmass` in absolute value.
+"""
+function resonanceBuilder(rmass::AbstractFloat, legs::AbstractVector{ReconstructedParticle})
+    result = ReconstructedParticle[]
+    length(legs) < 2 && return result
+    for (a,b) in combinations(legs, 2)
+        lv = LorentzVector(a.energy, a.momentum...) + LorentzVector(b.energy, b.momentum...)
+        rcharge = a.charge + b.charge
+        push!(result, ReconstructedParticle(mass=mass(lv), momentum=(lv.x, lv.y, lv.z), charge=rcharge))
+    end
+    sort!(result, lt =  (a,b) -> abs(rmass-a.mass) < abs(rmass-b.mass))
+    return result[1:1]  # take the best one
+end
+
+"""
+    recoilBuilder(comenergy::AbstractFloat, legs::AbstractVector{ReconstructedParticle})
+
+    build the recoil from an arbitrary list of input `ReconstructedParticle`s and the center of mass energy.
+"""
+function recoilBuilder(comenergy::AbstractFloat, in::AbstractVector{ReconstructedParticle})
+    result = ReconstructedParticle[]
+    isempty(in) && return result
+    recoil_lv = LorentzVector(comenergy, 0, 0, 0)
+    for p in in
+        recoil_lv -= LorentzVector(p.mass, p.momentum...)
+    end
+    push!(result, ReconstructedParticle(mass=mass(recoil_lv), momentum=(recoil_lv.x, recoil_lv.y, recoil_lv.z)))
+    return result
+end
+

examples/FCC/analysis_mH-recoil.jl

@@ -0,0 +1,65 @@
+using Revise
+using EDM4hep
+using EDM4hep.RootIO
+using EDM4hep.SystemOfUnits
+using EDM4hep.Histograms
+using Parameters
+
+using Plots
+
+include("analysis_functions.jl")
+
+f = "root://eospublic.cern.ch//eos/experiment/fcc/ee/generation/DelphesEvents/winter2023/IDEA/p8_ee_ZZ_ecm240/events_000189367.root"
+
+reader = RootIO.Reader(f);
+events = RootIO.get(reader, "events");
+
+@with_kw struct Histograms
+    mz          = H1D("m_{Z} [GeV]",125,0,250, unit=:GeV)
+    mz_zoom     = H1D("m_{Z} [GeV]",40,80,100, unit=:GeV)
+    lr_m        = H1D("Z leptonic recoil [GeV]", 100, 0, 200, unit=:GeV)
+    lr_m_zoom   = H1D("Z leptonic recoil [GeV]", 200, 80, 160, unit=:GeV)
+    lr_m_zoom1  = H1D("Z leptonic recoil [GeV]", 100, 120, 140, unit=:GeV)
+    lr_m_zoom2  = H1D("Z leptonic recoil [GeV]", 200, 120, 140, unit=:GeV)
+    lr_m_zoom3  = H1D("Z leptonic recoil [GeV]", 400, 120, 140, unit=:GeV)
+    lr_m_zoom4  = H1D("Z leptonic recoil [GeV]", 800, 120, 140, unit=:GeV)
+    lr_m_zoom5  = H1D("Z leptonic recoil [GeV]", 2000, 120, 140, unit=:GeV)
+    lr_m_zoom6  = H1D("Z leptonic recoil [GeV]", 100, 130.3, 132.5, unit=:GeV)
+end
+function do_plot(histos::Histograms)
+    img = plot(layout=(3,4), show=true, size=(1400,1000))
+    for (i,fn) in enumerate(fieldnames(Histograms))
+        h = getfield(histos, fn)
+        plot!(subplot=i, h.hist, title=h.title, show=true, cgrad=:plasma)
+    end
+    return img
+end
+
+myhists = Histograms()
+
+for evt in events
+    recps = RootIO.get(reader, evt, "ReconstructedParticles");
+    muons = RootIO.get(reader, evt, "Muon#0"; btype=ObjectID{ReconstructedParticle})
+    sel_muons = filter(x -> pₜ(x) > 10GeV, muons)
+    zed_leptonic = resonanceBuilder(91GeV, sel_muons)
+    zed_leptonic_recoil = recoilBuilder(240GeV, zed_leptonic)
+    if length(zed_leptonic) == 1    #  Filter to have exactly one Z candidate
+        Zcand_m        = zed_leptonic[1].mass
+        Zcand_recoil_m = zed_leptonic_recoil[1].mass
+        Zcand_q        = zed_leptonic[1].charge
+        if 80GeV <= Zcand_m <= 100GeV
+            #---Fill histograms now--------------------------------------
+            push!(myhists.mz, Zcand_m)
+            push!(myhists.mz_zoom, Zcand_m)
+            push!(myhists.lr_m, Zcand_recoil_m)
+            push!(myhists.lr_m_zoom1, Zcand_recoil_m)
+            push!(myhists.lr_m_zoom2, Zcand_recoil_m)
+            push!(myhists.lr_m_zoom3, Zcand_recoil_m)
+            push!(myhists.lr_m_zoom4, Zcand_recoil_m)
+            push!(myhists.lr_m_zoom5, Zcand_recoil_m)
+            push!(myhists.lr_m_zoom6, Zcand_recoil_m)
+        end
+    end
+end
+
+do_plot(myhists)

src/Datatypes.jl

@@ -35,3 +35,10 @@ function Base.getproperty(obj::MCParticle, sym::Symbol)
         return getfield(obj, sym)
     end
 end
+
+#--------------------------------------------------------------------------------------------------
+#----Utility functions for ReconstructedParticle---------------------------------------------------
+#--------------------------------------------------------------------------------------------------
+export pₜ
+pₜ( o::ReconstructedParticle) = √(o.momentum.x^2 + o.momentum.y^2)
+

src/EDM4hep.jl

@@ -12,5 +12,7 @@ module EDM4hep
     include("Datatypes.jl")
     include("EDStore.jl")
     include("RootIO.jl")
-
+    include("SystemOfUnits.jl")
+    include("Histograms.jl")
+
 end # module EDM4hep

src/Histograms.jl

@@ -0,0 +1,65 @@
+module Histograms
+    using FHist
+    using EDM4hep
+    using EDM4hep.SystemOfUnits
+
+    export H1D, H2D, H3D
+
+    const Hist1DType = typeof(Hist1D(bins=range(0,1,10)))
+    const Hist2DType = typeof(Hist2D(bins=(range(0,1,10),range(0,1,10))))
+    const Hist3DType = typeof(Hist3D(bins=(range(0,1,10),range(0,1,10),range(0,1,10))))
+
+    _getvalue(unit::Symbol) = getfield(EDM4hep.SystemOfUnits, unit)
+    _getvalue(units::Tuple{Symbol,Symbol}) = (getfield(EDM4hep.SystemOfUnits, units[1]), getfield(EDM4hep.SystemOfUnits, units[2]))
+    _getvalue(units::Tuple{Symbol,Symbol,Symbol}) = (getfield(EDM4hep.SystemOfUnits, units[1]), getfield(EDM4hep.SystemOfUnits, units[2]), getfield(EDM4hep.SystemOfUnits, units[3]))
+
+    """
+    H1D(title::String, nbins::Int, min::Float, max::Float, unit::Symbol)
+        Create a 1-dimentional histgram carrying the title and units
+    """
+    struct H1D
+        title::String
+        hist::Hist1DType
+        usym::Symbol
+        uval::Float64
+        H1D(title, nbins, min, max; unit=:nounit) = new(title, Hist1D(Float64; bins=range(min,max,nbins+1), overflow=false), unit, _getvalue(unit))
+    end
+
+    Base.push!(h::H1D, v, w=1) = atomic_push!(h.hist, v/h.uval, w)
+    Base.merge!(h1::H1D, h2::H1D) = merge!(h1.hist, h2.hist)
+    Base.empty!(h1::H1D) = empty!(h1.hist)
+
+    """
+    H2D(title::String, xbins::Int, xmin::Float, xmax::Float, ybins::Int, ymin::Float, ymax::Float, unit::Tuple{Symbol,Symbol})
+        Create a 2-dimentional histgram carrying the title and units
+    """
+    struct H2D
+        title::String
+        hist::Hist2DType
+        unit::Tuple{Symbol,Symbol}
+        uval::Tuple{Float64,Float64}
+        H2D(title, xbins, xmin, xmax, ybins, ymin, ymax; units=(:nounit, :nounit)) = new(title, Hist2D(Float64;bins=(range(xmin,xmax,xbins+1), range(ymin,ymax, ybins+1)), overflow=true), units, _getvalue(units))
+    end
+
+    Base.push!(h::H2D, u, v, w=1) = atomic_push!(h.hist, u/h.uval[1], v/h.uval[2], w)
+    Base.merge!(h1::H2D, h2::H2D) = merge!(h1.hist, h2.hist) 
+    Base.empty!(h1::H2D) = empty!(h1.hist)
+
+    """
+    H3D(title::String, xbins::Int, xmin::Float, xmax::Float, ybins::Int, ymin::Float, ymax::Float, zbins::Int, zmin::Float, zmax::Float, unit::Tuple{Symbol,Symbol,Symbol})
+        Create a 2-dimentional histgram carrying the title and units
+    """
+    struct H3D
+        title::String
+        hist::Hist3DType
+        unit::Tuple{Symbol,Symbol,Symbol}
+        uval::Tuple{Float64,Float64,Float64}
+        H3D(title, xbins, xmin, xmax, ybins, ymin, ymax, zbins, zmin, zmax; units=(:nounit, :nounit, :nounit)) = 
+            new(title, Hist3D(Float64;bins=(range(xmin,xmax,xbins+1), range(ymin,ymax, ybins+1), range(zmin,zmax, zbins+1)), overflow=true), units, _getvalue(units))
+    end
+
+    Base.push!(h::H3D, x, y, z, w=1) = atomic_push!(h.hist, x/h.uval[1], y/h.uval[2], z/h.uval[3], w)
+    Base.merge!(h1::H3D, h2::H3D) = merge!(h1.hist, h2.hist) 
+    Base.empty!(h1::H3D) = empty!(h1.hist)
+
+end

src/RootIO.jl

@@ -360,10 +360,10 @@ module RootIO
     """
     function get(reader::Reader, evt::UnROOT.LazyEvent, bname::String; btype::Type=Any, register=true)
         btype = btype === Any ? reader.btypes[bname] : btype     # Allow the user to force the actual type
-        if btype == ObjectID
-            register=false                    # Do not register a collection of ObjectIDs
-            sa = _get(reader, evt, bname, ObjectID{EDM4hep.POD}, register)
-            return convert.(eltype(eltype(sa)), sa)
+        if btype <: ObjectID
+            ED = eltype(btype)
+            sa = _get(reader, evt, bname, btype, false)
+            return StructArray(ED[convert(ED, oid) for oid in sa])     # explicitly convert to the pointed objects
         else
             sa = _get(reader, evt, bname, btype, register)
             return sa