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, DistributionsDefine 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) ≈ 1trueDefine 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.0025676779122516535For splines
nw = QCDNUM.zmfillw()
splint_params = QCDNUM.SPLINTParams();
quark_coeffs = QuarkCoefficients();
ZMFILLW: start weight calculations 4 38 0 0
ZMFILLW: calculations completedDefine 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 ~/.julia/packages/QCDNUM/gx4Mo/src/splint/initialisation.jl:36Check 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) # charge = 1
input_xsecp = @cfunction(my_funcp, Float64, (Ref{Int32}, Ref{Int32}, Ref{UInt8}))
my_funcm = get_input_xsec_func(-1) # charge = -1
input_xsecm = @cfunction(my_funcm, Float64, (Ref{Int32}, Ref{Int32}, Ref{UInt8}))Ptr{Nothing} @0x00007f86cd1f7cd0plot
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, 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(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, xbins_M_begin[i], xbins_M_end[i],
q2bins_M_begin[i], q2bins_M_end[i], 318.0, 4)
IntXsec_eM[i] = QCDNUM.dsp_ints2(iaF_eM, xbins_M_begin[i], xbins_M_end[i],
q2bins_M_begin[i], q2bins_M_end[i], 318.0, 4)
end1 for e-p and 0 for e+p
ePp = 0;
eMp = 1;
TM_eP = get_TM_elements(ePp);
TM_eM = get_TM_elements(eMp);
K_eP = get_K_elements(ePp);
K_eM = get_K_elements(eMp);
nbins_out = size(TM_eP)[2];
counts_pred_eP = zeros(nbins_out);
counts_pred_eM = zeros(nbins_out);
for j in 1:nbins_out
for i in 1:nbins
counts_pred_eP[j] += TM_eP[i, j] * (1.0 / K_eP[i]) * IntXsec_eP[i]
counts_pred_eM[j] += TM_eM[i, j] * (1.0 / K_eM[i]) * IntXsec_eM[i]
end
end
counts_pred_eM153-element Vector{Float64}:
782.0747989157136
543.341827637819
407.6562006131687
343.27554129757675
221.58522672609953
194.89598293573863
140.65397196560687
278.93744351583194
853.6716286836335
696.8839416798988
⋮
13.826820303419813
7.897518762639923
4.452160806595315
2.3743347283583263
1.413392410229477
0.9056333285465725
0.09650371692146183
52.806123532812435
2.065405204640709This page was generated using Literate.jl.