Forward model

Here, we go through an example of simulating the full forward model, from the prior definition to the expected number of events in different bins of the detector response.

using QCDNUM, PartonDensity
using Plots, Printf, NaNMath, Parameters, Random, Distributions

zeus_include_path = string(chop(pathof(PartonDensity), tail=20), "data/ZEUS_I1787035/ZEUS_I1787035.jl")

const MD_DOCS = include(zeus_include_path)

PartonDensity.MetaDataZEUS{Main.var"#f_cross_section_to_counts#1"{Array{Float64, 3}, Array{Float64, 3}, Matrix{Float64}, Matrix{Float64}, Float64, Float64, Vector{Float64}, Vector{Float64}}}("ZEUS_I1787035", 141.44, 185.018, 0.018, 0.018, 318.0, [400.0, 400.0, 400.0, 400.0, 400.0, 400.0, 400.0, 400.0, 400.0, 400.0  …  8000.0, 8000.0, 8000.0, 8000.0, 8000.0, 11000.0, 11000.0, 11000.0, 20000.0, 20000.0], [650.0, 650.0, 650.0, 650.0, 650.0, 650.0, 650.0, 650.0, 650.0, 650.0  …  11000.0, 11000.0, 11000.0, 11000.0, 11000.0, 20000.0, 20000.0, 20000.0, 30000.0, 30000.0], [0.0035, 0.007, 0.01, 0.012, 0.015, 0.017, 0.019, 0.022, 0.024, 0.026  …  0.45, 0.51, 0.57, 0.64, 0.78, 0.13, 0.256, 0.6, 0.25, 0.6], [0.007, 0.01, 0.012, 0.015, 0.017, 0.019, 0.022, 0.024, 0.026, 0.028  …  0.51, 0.57, 0.64, 0.78, 0.99999, 0.256, 0.6, 0.99999, 0.6, 0.99999], [743, 580, 441, 416, 283, 248, 227, 504, 789, 681  …  27, 19, 12, 8, 5, 4, 1, 1, 120, 8], [567, 418, 321, 274, 207, 176, 152, 371, 599, 497  …  15, 14, 9, 8, 7, 2, 0, 0, 41, 3], Main.var"#f_cross_section_to_counts#1"{Array{Float64, 3}, Array{Float64, 3}, Matrix{Float64}, Matrix{Float64}, Float64, Float64, Vector{Float64}, Vector{Float64}}([0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], [0.0 0.0 … 0.0 0.0; 0.003688148812 0.0018500689892 … 0.0 0.0; … ; 0.0 0.0 … 1.0453516999999999 0.0005606785471999999; 0.0 0.0 … 0.0046826205620000005 0.065218845], [0.0 0.0 … 0.0 0.0; 0.00234097344 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.94503136 0.000449623616; 0.0 0.0 … 0.08914963199999999 0.0282215232], 185.018, 141.44, [1.07559, 1.00494, 0.997145, 0.994087, 0.991775, 0.990224, 0.98391, 0.991145, 0.993283, 0.976963  …  0.897826, 0.889232, 0.890147, 0.892718, 0.88395, 0.921764, 0.882477, 0.891006, 0.911363, 0.877797], [1.07535, 1.00332, 1.0004, 0.990515, 0.982716, 0.991055, 0.988341, 0.981917, 0.980651, 0.98678  …  0.857518, 0.837175, 0.875777, 0.897924, 0.935441, 0.933215, 0.899119, 0.910332, 1.1251, 1.12007]))

Define input PDFs

We can use DirichletPDFParams or ValencePDFParams, as long as we do so according to the PDF parametrisation and priors docs.

random_seed = 42

weights = [3.0, 1.0, 5.0, 5.0, 1.0, 1.0, 1.0, 0.5, 0.5]
θ = rand(MersenneTwister(random_seed), Dirichlet(weights))
pdf_params = DirichletPDFParams(K_u=4.0, K_d=6.0, λ_g1=0.7, λ_g2=-0.4,
    K_g=6.0, λ_q=-0.5, K_q=5.0, θ=θ);

Plot the input PDFs

plot_input_pdfs(pdf_params)

Sanity check that sum = 1

int_xtotx(pdf_params) ≈ 1

true

Define QCDNUM grids, weights and settings

grid = QCDNUM.GridParams(x_min=[1.0e-3], x_weights=[1], nx=100,
    qq_bounds=[1.0e2, 3.0e4], qq_weights=[1.0, 1.0], nq=50, spline_interp=3)

qcdnum_params = QCDNUM.EvolutionParams(order=2, α_S=0.118, q0=100.0, grid_params=grid,
    n_fixed_flav=5, iqc=1, iqb=1, iqt=1, weight_type=1);

Initialise

QCDNUM.init()

  +---------------------------------------------------------------------+
  |                                                                     |
  |    If you use QCDNUM, please refer to:                              |
  |                                                                     |
  |    M. Botje, Comput. Phys. Commun. 182(2011)490, arXiV:1005.1481    |
  |                                                                     |
  +---------------------------------------------------------------------+

Evolve the PDFs using QCDNUM

Define input PDF function

• See QCDNUM docs under evolfg

There are functions available to generate the necessary input PDF function in the correct format for QCDNUM.jl (see get_input_pdf_func()), along with the mapping between this input function and quark species (see input_pdf_map).

Get function and PDF input map to fully describe the QCDNUM.InputPDF

my_func = get_input_pdf_func(pdf_params)
input_pdf = QCDNUM.InputPDF(func=my_func, map=input_pdf_map)

QCDNUM.InputPDF
  func: _input_pdfs (function of type PartonDensity.var"#_input_pdfs#22"{DirichletPDFParams{Float64, Vector{Float64}}})
  cfunc: Base.CFunction
  map: Array{Float64}((156,)) [0.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

Evolve the PDF over the specified grid

ϵ = QCDNUM.evolve(input_pdf, qcdnum_params)

0.0023831637076243695

For splines

nw = QCDNUM.zmfillw()
splint_params = QCDNUM.SPLINTParams();
quark_coeffs = QuarkCoefficients();

 ZMFILLW: start weight calculations   4  38   0   0
 ZMFILLW: calculations completed

Define initial spline

QCDNUM.ssp_spinit(splint_params.nuser)

ia = QCDNUM.isp_s2make(splint_params.nsteps_x, splint_params.nsteps_q);
xnd = QCDNUM.ssp_unodes(ia, splint_params.nnodes_x, 0);
qnd = QCDNUM.ssp_vnodes(ia, splint_params.nnodes_q, 0);

┌ Warning: SPLINT is already initialised, skipping call to QCDNUM.ssp_spinit()
└ @ QCDNUM /opt/myjulia/packages/QCDNUM/gx4Mo/src/splint/initialisation.jl:36

Check nodes and erase

QCDNUM.ssp_nprint(ia);
QCDNUM.ssp_erase(ia);

  N   IX     XNODE    NQMA      IQ     QNODE    NXMI
  1    1  0.10000E-02    0       1  0.10000E+03    0
  2    5  0.13183E-02    0       2  0.11234E+03    0
  3   10  0.18621E-02    0      12  0.35982E+03    0
  4   15  0.26303E-02    0      22  0.11525E+04    0
  5   20  0.37154E-02    0      32  0.36911E+04    0
  6   25  0.52481E-02    0      42  0.11822E+05    0
  7   30  0.74131E-02    0      50  0.30000E+05    0
  8   35  0.10471E-01    0       -       -         -
  9   40  0.14791E-01    0       -       -         -
 10   45  0.20893E-01    0       -       -         -
 11   50  0.29512E-01    0       -       -         -
 12   55  0.41687E-01    0       -       -         -
 13   60  0.58884E-01    0       -       -         -
 14   65  0.83176E-01    0       -       -         -
 15   70  0.11749E+00    0       -       -         -
 16   75  0.16596E+00    0       -       -         -
 17   80  0.23442E+00    0       -       -         -
 18   85  0.33113E+00    0       -       -         -
 19   90  0.46774E+00    0       -       -         -
 20   95  0.66069E+00    0       -       -         -
 21  100  0.93325E+00    0       -       -         -
 22  101  0.10000E+01    0       -       -         -

Set nodes and fill spline with structure function

iaFLup = QCDNUM.isp_s2user(xnd, splint_params.nnodes_x, qnd, splint_params.nnodes_q);
QCDNUM.ssp_s2f123(iaFLup, qcdnum_params.output_pdf_loc, quark_coeffs.proup, 1, 0.0);

iaF2up = QCDNUM.isp_s2user(xnd, splint_params.nnodes_x, qnd, splint_params.nnodes_q);
QCDNUM.ssp_s2f123(iaF2up, qcdnum_params.output_pdf_loc, quark_coeffs.proup, 2, 0.0);

iaF3up = QCDNUM.isp_s2user(xnd, splint_params.nnodes_x, qnd, splint_params.nnodes_q);
QCDNUM.ssp_s2f123(iaF3up, qcdnum_params.output_pdf_loc, quark_coeffs.valup, 3, 0.0);

iaFLdn = QCDNUM.isp_s2user(xnd, splint_params.nnodes_x, qnd, splint_params.nnodes_q);
QCDNUM.ssp_s2f123(iaFLdn, qcdnum_params.output_pdf_loc, quark_coeffs.prodn, 1, 0.0);

iaF2dn = QCDNUM.isp_s2user(xnd, splint_params.nnodes_x, qnd, splint_params.nnodes_q);
QCDNUM.ssp_s2f123(iaF2dn, qcdnum_params.output_pdf_loc, quark_coeffs.prodn, 2, 0.0);

iaF3dn = QCDNUM.isp_s2user(xnd, splint_params.nnodes_x, qnd, splint_params.nnodes_q);
QCDNUM.ssp_s2f123(iaF3dn, qcdnum_params.output_pdf_loc, quark_coeffs.valdn, 3, 0.0);

store spline addresses

QCDNUM.ssp_uwrite(splint_params.spline_addresses.F2up, Float64(iaF2up));
QCDNUM.ssp_uwrite(splint_params.spline_addresses.F2dn, Float64(iaF2dn));
QCDNUM.ssp_uwrite(splint_params.spline_addresses.F3up, Float64(iaF3up));
QCDNUM.ssp_uwrite(splint_params.spline_addresses.F3dn, Float64(iaF3dn));
QCDNUM.ssp_uwrite(splint_params.spline_addresses.FLup, Float64(iaFLup));
QCDNUM.ssp_uwrite(splint_params.spline_addresses.FLdn, Float64(iaFLdn));

my_funcp = get_input_xsec_func(1, MD_DOCS) # charge = 1
input_xsecp = @cfunction(my_funcp, Float64, (Ref{Int32}, Ref{Int32}, Ref{UInt8}))

my_funcm = get_input_xsec_func(-1, MD_DOCS) # charge = -1
input_xsecm = @cfunction(my_funcm, Float64, (Ref{Int32}, Ref{Int32}, Ref{UInt8}))

Ptr{Nothing} @0x00007fad241c4e80

plot

g = qcdnum_params.grid_params
xsec_on_grid = zeros(g.nx, g.nq);

for ix = 1:g.nx
    for iq = 1:g.nq
        xsec_on_grid[ix, iq] = _fun_xsec_i(1, MD_DOCS, ix, iq) # charge = 1
    end
end

qcdnum_x_grid = QCDNUM.gxcopy(g.nx)
qcdnum_qq_grid = QCDNUM.gqcopy(g.nq)
p1 = heatmap(qcdnum_x_grid, qcdnum_qq_grid, NaNMath.log10.(xsec_on_grid[:, :]'))
plot(p1, xlabel="x", ylabel="q2",
    xaxis=:log, yaxis=:log)

plot(qcdnum_x_grid, NaNMath.log10.(xsec_on_grid[:, 4]),
    label="Q2=141 (input scale)", lw=3)
plot!(qcdnum_x_grid, NaNMath.log10.(xsec_on_grid[:, 22]), label="Q2=1152", lw=3)
plot!(qcdnum_x_grid, NaNMath.log10.(xsec_on_grid[:, 35]), label="Q2=5233", lw=3)
plot!(qcdnum_x_grid, NaNMath.log10.(xsec_on_grid[:, 41]), label="Q2=10523", lw=3)
plot!(xaxis=:log, legend=:bottomleft, xlabel="x",
    ylabel="log10(cross section spline input)", ylims=(-7, 5))

iaF_eP = QCDNUM.isp_s2make(1, 2);
QCDNUM.ssp_uwrite(splint_params.spline_addresses.F_eP, Float64(iaF_eP));
QCDNUM.ssp_s2fill(iaF_eP, input_xsecp, splint_params.rscut);

iaF_eM = QCDNUM.isp_s2make(1, 2);
QCDNUM.ssp_uwrite(splint_params.spline_addresses.F_eM, Float64(iaF_eM));
QCDNUM.ssp_s2fill(iaF_eM, input_xsecm, splint_params.rscut);

plot spline

spline = zeros(g.nx, g.nq);

for ix = 1:g.nx
    for iq = 1:g.nq
        spline[ix, iq] = QCDNUM.dsp_funs2(iaF_eP, qcdnum_x_grid[ix],
            qcdnum_qq_grid[iq], 1)
    end
end

p1 = heatmap(qcdnum_x_grid, qcdnum_qq_grid, NaNMath.log10.(spline[:, :]'))
plot(p1, xlabel="x", ylabel="q2",
    xaxis=:log, yaxis=:log)

Integrate over the cross section spline and find expected events numbers

Here, we neglect any possible contribution from systematic errors

nbins = size(MD_DOCS.m_xbins_M_begin)[1]
IntXsec_eP = zeros(nbins);
IntXsec_eM = zeros(nbins);
for i in 1:nbins
    IntXsec_eP[i] = QCDNUM.dsp_ints2(iaF_eP, MD_DOCS.m_xbins_M_begin[i], MD_DOCS.m_xbins_M_end[i], MD_DOCS.m_q2bins_M_begin[i], MD_DOCS.m_q2bins_M_end[i], MD_DOCS.sqrtS, 4)
    IntXsec_eM[i] = QCDNUM.dsp_ints2(iaF_eM, MD_DOCS.m_xbins_M_begin[i], MD_DOCS.m_xbins_M_end[i], MD_DOCS.m_q2bins_M_begin[i], MD_DOCS.m_q2bins_M_end[i], MD_DOCS.sqrtS, 4)
end

counts_pred_eP, counts_pred_eM = MD_DOCS.f_cross_section_to_counts(MD_DOCS.Ld_ePp, MD_DOCS.Ld_eMp, IntXsec_eP, IntXsec_eM)

counts_pred_eM

153-element Vector{Float64}:
 651.6430308810627
 462.3474434974716
 353.8250890402993
 304.3408226976592
 200.58796784717183
 180.09268590111463
 132.68825877655811
 276.3002068237758
 704.4763532747882
 587.3200097391523
   ⋮
  14.133048867826203
   8.186102087178547
   4.6763201873520055
   2.522324578546638
   1.5184036052034076
   0.9856294574716068
   0.1073669154427297
  53.9465198773069
   2.211269559843675

This page was generated using Literate.jl.