Examples

Simple QCD evolution and tabulation

 1#!/usr/bin/env python3
 2"""Example using HOPPET in Python to perform PDF evolution and tabulation."""
 3import hoppet as hp
 4import numpy as np
 5
 6# This is the PDF subroutine. Notice that the signature is different
 7# from the one in the C++ code to allow for the python interface, that
 8# handles pointers in a more complicated way. In the actual example we
 9# use the internal BenchmarkPDFunpol which is identical
10
11#def hera_lhc(x, Q):
12#    N_g = 1.7
13#    N_ls = 0.387975
14#    N_uv = 5.107200
15#    N_dv = 3.064320
16#    N_db = N_ls / 2
17#
18#    uv = N_uv * pow(x, 0.8) * pow((1 - x), 3)
19#    dv = N_dv * pow(x, 0.8) * pow((1 - x), 4)
20#    dbar = N_db * pow(x, -0.1) * pow(1 - x, 6)
21#    ubar = dbar * (1 - x)
22#
23#    pdf = np.zeros(13) # Initialise the array with zeros
24#
25#    pdf[ 0+6] = N_g * pow(x, -0.1) * pow(1 - x, 5)
26#    pdf[-3+6] = 0.2 * (dbar + ubar)
27#    pdf[ 3+6] = pdf[-3 + 6]
28#    pdf[ 2+6] = uv + ubar
29#    pdf[-2+6] = ubar
30#    pdf[ 1+6] = dv + dbar
31#    pdf[-1+6] = dbar
32#    
33## Overkill
34#    pdf[ 4+6] = 0.0
35#    pdf[ 5+6] = 0.0
36#    pdf[ 6+6] = 0.0
37#    pdf[-4+6] = 0.0
38#    pdf[-5+6] = 0.0
39#    pdf[-6+6] = 0.0
40#
41#    return pdf
42
43def main():
44    dy = 0.1    
45    nloop = 3
46    # Start hoppet
47    hp.Start(dy, nloop)
48    
49    asQ0 = 0.35
50    Q0 = np.sqrt(2.0)
51    # Do the evolution (replace hp.BenchmarkPDFunpol with a function
52    # that acts like hera_lhc above to use your own initial condition)
53    hp.Evolve(asQ0, Q0, nloop, 1.0, hp.BenchmarkPDFunpol, Q0)
54
55    # Evaluate the PDFs at some x values and print them
56    xvals = [1e-5,1e-4,1e-3,1e-2,0.1,0.3,0.5,0.7,0.9]
57    Q = 100.0
58
59    print("")
60    print("           Evaluating PDFs at Q =",Q, " GeV")
61    print("    x      u-ubar      d-dbar    2(ubr+dbr)    c+cbar       gluon")
62    for ix in range(9):
63        pdf_array = hp.Eval(xvals[ix], Q)
64        print("{:7.1E} {:11.4E} {:11.4E} {:11.4E} {:11.4E} {:11.4E}".format(
65            xvals[ix],
66            pdf_array[6 + 2] - pdf_array[6 - 2], 
67            pdf_array[6 + 1] - pdf_array[6 - 1], 
68            2 * (pdf_array[6 - 1] + pdf_array[6 - 2]),
69            pdf_array[6 - 4] + pdf_array[6 + 4],
70            pdf_array[6 + 0]
71        ))
72
73    #hp.WriteLHAPDFGrid("test_python",0)
74
75    hp.DeleteAll()
76
77if __name__ == "__main__":
78    main()

QED evolution

  1#!/usr/bin/env python3
  2"""Example using HOPPET in Python to perform PDF evolution and
  3tabulation, including QED effects."""
  4import hoppet as hp
  5import numpy as np
  6
  7def hera_lhc(x, Q):
  8    """This routine provides the initial conditions for the PDF evolution,
  9    setting up the individual flavours at the specified x and Q values,
 10    and returning them in a numpy array. 
 11    
 12    In this instance, the Q value is ignored, but it is required by the
 13    underlying HOPPET interface.
 14
 15    Notice that the signature is different from that in the C++ code to
 16    allow for the python interface, which handles pointers in a more
 17    complicated way. 
 18    """
 19    N_g = 1.7
 20    N_ls = 0.387975
 21    N_uv = 5.107200
 22    N_dv = 3.064320
 23    N_db = N_ls / 2
 24
 25    gluon = N_g * pow(x, -0.1) * pow(1 - x, 5)
 26    uv = N_uv * pow(x, 0.8) * pow((1 - x), 3)
 27    dv = N_dv * pow(x, 0.8) * pow((1 - x), 4)
 28    dbar = N_db * pow(x, -0.1) * pow(1 - x, 6)
 29    ubar = dbar * (1 - x)
 30
 31    pdf = np.zeros(18) # Initialise the array with zeros
 32
 33    pdf[ 0+6] = 0.99 * gluon
 34    pdf[-3+6] = 0.2 * (dbar + ubar)
 35    pdf[ 3+6] = pdf[-3 + 6]
 36    pdf[ 2+6] = uv + 0.9 * ubar
 37    pdf[-2+6] = 0.9 * ubar
 38    pdf[ 1+6] = dv + 0.9 * dbar
 39    pdf[-1+6] = 0.9 * dbar
 40
 41    # qed part
 42    pdf[  8+6] = gluon / 100.0 + ubar / 10.0  # photon 
 43    pdf[  9+6] = dbar  / 10.0  # electron + positron
 44    pdf[ 10+6] = ubar  / 10.0  # mu+ + mu-
 45    pdf[ 11+6] = dbar  / 10.0  # tau+ + tau-
 46    
 47    # Make sure other flavours are zero
 48    pdf[ 4+6] = 0.0
 49    pdf[ 5+6] = 0.0
 50    pdf[ 6+6] = 0.0
 51    pdf[-4+6] = 0.0
 52    pdf[-5+6] = 0.0
 53    pdf[-6+6] = 0.0
 54
 55    return pdf
 56
 57def main():
 58    dy = 0.1    
 59    nloop = 3
 60    
 61    # set QED
 62    use_qcd_qed  = True
 63    use_Plq_nnlo = False
 64    hp.SetQED(True, use_qcd_qed, use_Plq_nnlo)
 65    
 66    # Start hoppet
 67    hp.StartExtended(12.0, dy, 1.0, 28000.0, dy/4.0, nloop, -6, hp.factscheme_MSbar)
 68
 69    # Set heavy flavour scheme
 70    mc = 1.414213563
 71    mb = 4.5
 72    mt = 175.0
 73    hp.SetVFN(mc, mb, mt)
 74
 75    
 76    asQ0 = 0.35
 77    Q0 = np.sqrt(2.0)
 78    # Do the evolution. 
 79    hp.Evolve(asQ0, Q0, nloop, 1.0, hera_lhc, Q0)
 80
 81    # Evaluate the PDFs at some x values and print them
 82    xvals = [1e-5,1e-4,1e-3,1e-2,0.1,0.3,0.5,0.7,0.9]
 83    Q = 100.0
 84    
 85    print("")
 86    print("           Evaluating PDFs at Q =",Q, " GeV")
 87    print("    x      u-ubar      d-dbar    2(ubr+dbr)    c+cbar       gluon      photon     e+ + e-    mu+ + mu-   tau+ + tau-")
 88    for ix in range(9):
 89        pdf_array = hp.Eval(xvals[ix], Q)
 90        print("{:7.1E} {:11.4E} {:11.4E} {:11.4E} {:11.4E} {:11.4E} {:11.4E} {:11.4E} {:11.4E} {:11.4E}".format(
 91            xvals[ix],
 92            pdf_array[6 + 2] - pdf_array[6 - 2], 
 93            pdf_array[6 + 1] - pdf_array[6 - 1], 
 94            2 * (pdf_array[6 - 1] + pdf_array[6 - 2]),
 95            pdf_array[6 - 4] + pdf_array[6 + 4],
 96            pdf_array[6 + 0], pdf_array[6 + 8], pdf_array[6 + 9], pdf_array[6 + 10], pdf_array[6 + 11]
 97        ))
 98
 99    #hp.WriteLHAPDFGrid("test_python",0)
100
101    hp.DeleteAll()
102
103if __name__ == "__main__":
104    main()

Structure function tabulation

  1#!/usr/bin/env python3
  2"""Example using HOPPET in Python to compute structure functions."""
  3import hoppet as hp
  4import numpy as np
  5
  6def main():
  7    Qmax = 13000.0
  8    nloop_coefs = 4
  9    xmur = 1.0
 10    xmuf = 1.0
 11    sc_choice = hp.scale_choice_Q # Uses Q as the central scale choice
 12    Qmin = 1.0
 13
 14    # Set heavy flavour scheme
 15    mc = 1.414213563  # sqrt(2.0_dp) + epsilon
 16    mb = 4.5
 17    mt = 175.0
 18    hp.SetPoleMassVFN(mc, mb, mt)
 19
 20    # Streamlined initialization
 21    # including  parameters for x-grid
 22    order = -6 # interpolation order, not perturbative order in alphas!
 23    ymax  = 16.0
 24    dy    = 0.05  # dble_val_opt("-dy",0.1_dp)
 25    dlnlnQ = dy/4.0
 26    nloop = 3 
 27    minQval = min(xmuf*Qmin, Qmin)
 28    maxQval = max(xmuf*Qmax, Qmax)
 29
 30    # initialise the grid and dglap holder, using the streamlined
 31    # interface for simplicity
 32    hp.StartExtended(ymax,dy,minQval,maxQval,dlnlnQ,nloop,order,hp.factscheme_MSbar)
 33
 34    # Setup all constants and parameters needed by the structure functions
 35    nflav= 5
 36    hp.StartStrFctExtended(nloop_coefs,nflav,sc_choice,1.0,True,80.377,91.1876)
 37    
 38    asQ0 = 0.35
 39    Q0 = np.sqrt(2.0)
 40    # Do the evolution. 
 41    hp.Evolve(asQ0, Q0, nloop, 1.0, hp.BenchmarkPDFunpol, Q0)
 42
 43    hp.InitStrFct(nloop_coefs, True, xmur, xmuf)
 44
 45    # write out the structure functions
 46    ymax = np.log(1e5) 
 47    Q = 100.0
 48    write_f1_py(Q, ymax, 10, Q, Q, True)
 49    write_f2_py(Q, ymax, 10, Q, Q, True)
 50    write_f3_py(Q, ymax, 10, Q, Q, True)
 51
 52    hp.DeleteAll()
 53
 54def write_f1_py(Qtest, ymax, ny, muF=None, muR=None, use_sep_orders=False):
 55    print(f"# Q = {Qtest:10.4f}")
 56    if use_sep_orders:
 57        print("#     x        F1Wp(LO)    F1Wm(LO)    F1Wp(NLO)   F1Wm(NLO)  F1Wp(NNLO)  F1Wm(NNLO)"
 58              "  F1Wp(N3LO)  F1Wm(N3LO)    F1Z(LO)    F1Z(NLO)    F1Z(NNLO)"
 59              "   F1Z(N3LO)    F1γ(LO)    F1γ(NLO)    F1γ(NNLO)   F1γ(N3LO)   F1γZ(LO)"
 60              "   F1γZ(NLO)   F1γZ(NNLO)  F1γZ(N3LO)")
 61    else:
 62        print("# x F1Wp F1Wm F1Z F1γ F1γZ")
 63
 64    for iy in range(ny, 0, -1):
 65        ytest = iy * ymax / ny
 66        xval = np.exp(-ytest)
 67        if use_sep_orders:
 68            res_lo = hp.StrFctLO(xval, Qtest, muR, muF)
 69            print(f"{xval:12.4E} {res_lo[hp.iF1Wp]:11.4E} {res_lo[hp.iF1Wm]:11.4E}", end=' ')
 70            F1Z_LO = res_lo[hp.iF1Z]
 71            F1EM_LO = res_lo[hp.iF1EM]
 72            F1gZ_LO = res_lo[hp.iF1gZ]
 73
 74            res_nlo = hp.StrFctNLO(xval, Qtest, muR, muF)
 75            print(f"{res_nlo[hp.iF1Wp]:11.4E} {res_nlo[hp.iF1Wm]:11.4E}", end=' ')
 76            F1Z_NLO = res_nlo[hp.iF1Z]
 77            F1EM_NLO = res_nlo[hp.iF1EM]
 78            F1gZ_NLO = res_nlo[hp.iF1gZ]
 79
 80            res_nnlo = hp.StrFctNNLO(xval, Qtest, muR, muF)
 81            print(f"{res_nnlo[hp.iF1Wp]:11.4E} {res_nnlo[hp.iF1Wm]:11.4E}", end=' ')
 82            F1Z_NNLO = res_nnlo[hp.iF1Z]
 83            F1EM_NNLO = res_nnlo[hp.iF1EM]
 84            F1gZ_NNLO = res_nnlo[hp.iF1gZ]
 85
 86            res_n3lo = hp.StrFctN3LO(xval, Qtest, muR, muF)
 87            print(f"{res_n3lo[hp.iF1Wp]:11.4E} {res_n3lo[hp.iF1Wm]:11.4E}", end=' ')
 88            F1Z_N3LO = res_n3lo[hp.iF1Z]
 89            F1EM_N3LO = res_n3lo[hp.iF1EM]
 90            F1gZ_N3LO = res_n3lo[hp.iF1gZ]
 91
 92            print(f"{F1Z_LO:11.4E} {F1Z_NLO:11.4E} {F1Z_NNLO:11.4E} {F1Z_N3LO:11.4E} "
 93                  f"{F1EM_LO:11.4E} {F1EM_NLO:11.4E} {F1EM_NNLO:11.4E} {F1EM_N3LO:11.4E} "
 94                  f"{F1gZ_LO:11.4E} {F1gZ_NLO:11.4E} {F1gZ_NNLO:11.4E} {F1gZ_N3LO:11.4E}")
 95        else:
 96            res = hp.StrFct(xval, Qtest, muR, muF)
 97            print(f"{xval:12.4E} {res[hp.iF1Wp]:11.4E} {res[hp.iF1Wm]:11.4E} {res[hp.iF1Z]:11.4E} {res[hp.iF1EM]:11.4E} {res[hp.iF1gZ]:11.4E}")
 98    print()
 99    print()
100
101def write_f2_py(Qtest, ymax, ny, muF=None, muR=None, use_sep_orders=False):
102    print(f"# Q =           {Qtest:10.4f}")
103    if use_sep_orders:
104        print("#     x        F2Wp(LO)    F2Wm(LO)    F2Wp(NLO)   F2Wm(NLO)  F2Wp(NNLO)  F2Wm(NNLO)"
105              "  F2Wp(N3LO)  F2Wm(N3LO)    F2Z(LO)    F2Z(NLO)    F2Z(NNLO)"
106              "   F2Z(N3LO)    F2γ(LO)    F2γ(NLO)    F2γ(NNLO)   F2γ(N3LO)   F2γZ(LO)"
107              "   F2γZ(NLO)   F2γZ(NNLO)  F2γZ(N3LO)")
108    else:
109        print("# x F2Wp F2Wm F2Z F2γ F2γZ")
110
111    for iy in range(ny, 0, -1):
112        ytest = iy * ymax / ny
113        xval = np.exp(-ytest)
114        if use_sep_orders:
115            res_lo = hp.StrFctLO(xval, Qtest, muR, muF)
116            print(f"{xval:12.4E} {res_lo[hp.iF2Wp]:11.4E} {res_lo[hp.iF2Wm]:11.4E}", end=' ')
117            F2Z_LO = res_lo[hp.iF2Z]
118            F2EM_LO = res_lo[hp.iF2EM]
119            F2gZ_LO = res_lo[hp.iF2gZ]
120
121            res_nlo = hp.StrFctNLO(xval, Qtest, muR, muF)
122            print(f"{res_nlo[hp.iF2Wp]:11.4E} {res_nlo[hp.iF2Wm]:11.4E}", end=' ')
123            F2Z_NLO = res_nlo[hp.iF2Z]
124            F2EM_NLO = res_nlo[hp.iF2EM]
125            F2gZ_NLO = res_nlo[hp.iF2gZ]
126
127            res_nnlo = hp.StrFctNNLO(xval, Qtest, muR, muF)
128            print(f"{res_nnlo[hp.iF2Wp]:11.4E} {res_nnlo[hp.iF2Wm]:11.4E}", end=' ')
129            F2Z_NNLO = res_nnlo[hp.iF2Z]
130            F2EM_NNLO = res_nnlo[hp.iF2EM]
131            F2gZ_NNLO = res_nnlo[hp.iF2gZ]
132
133            res_n3lo = hp.StrFctN3LO(xval, Qtest, muR, muF)
134            print(f"{res_n3lo[hp.iF2Wp]:11.4E} {res_n3lo[hp.iF2Wm]:11.4E}", end=' ')
135            F2Z_N3LO = res_n3lo[hp.iF2Z]
136            F2EM_N3LO = res_n3lo[hp.iF2EM]
137            F2gZ_N3LO = res_n3lo[hp.iF2gZ]
138
139            print(f"{F2Z_LO:11.4E} {F2Z_NLO:11.4E} {F2Z_NNLO:11.4E} {F2Z_N3LO:11.4E} "
140                  f"{F2EM_LO:11.4E} {F2EM_NLO:11.4E} {F2EM_NNLO:11.4E} {F2EM_N3LO:11.4E} "
141                  f"{F2gZ_LO:11.4E} {F2gZ_NLO:11.4E} {F2gZ_NNLO:11.4E} {F2gZ_N3LO:11.4E}")
142        else:
143            res = hp.StrFct(xval, Qtest, muR, muF)
144            print(f"{xval:12.4E} {res[hp.iF2Wp]:11.4E} {res[hp.iF2Wm]:11.4E} {res[hp.iF2Z]:11.4E} {res[hp.iF2EM]:11.4E} {res[hp.iF2gZ]:11.4E}")
145    print()
146    print()
147
148def write_f3_py(Qtest, ymax, ny, muF=None, muR=None, use_sep_orders=False):
149    print(f"# Q =                     {Qtest:10.4f}")
150    if use_sep_orders:
151        print("#     x        F3Wp(LO)    F3Wm(LO)    F3Wp(NLO)   F3Wm(NLO)  F3Wp(NNLO)  F3Wm(NNLO)"
152              "  F3Wp(N3LO)  F3Wm(N3LO)    F3Z(LO)    F3Z(NLO)    F3Z(NNLO)"
153              "   F3Z(N3LO)    F3γZ(LO)"
154              "   F3γZ(NLO)  F3γZ(NNLO)  F3γZ(N3LO)")
155    else:
156        print("# x F3Wp F3Wm F3Z F3γZ")
157
158    for iy in range(ny, 0, -1):
159        ytest = iy * ymax / ny
160        xval = np.exp(-ytest)
161        if use_sep_orders:
162            res_lo = hp.StrFctLO(xval, Qtest, muR, muF)
163            print(f"{xval:12.4E} {res_lo[hp.iF3Wp]:11.4E} {res_lo[hp.iF3Wm]:11.4E}", end=' ')
164            F3Z_LO = res_lo[hp.iF3Z]
165            F3gZ_LO = res_lo[hp.iF3gZ]
166
167            res_nlo = hp.StrFctNLO(xval, Qtest, muR, muF)
168            print(f"{res_nlo[hp.iF3Wp]:11.4E} {res_nlo[hp.iF3Wm]:11.4E}", end=' ')
169            F3Z_NLO = res_nlo[hp.iF3Z]
170            F3gZ_NLO = res_nlo[hp.iF3gZ]
171
172            res_nnlo = hp.StrFctNNLO(xval, Qtest, muR, muF)
173            print(f"{res_nnlo[hp.iF3Wp]:11.4E} {res_nnlo[hp.iF3Wm]:11.4E}", end=' ')
174            F3Z_NNLO = res_nnlo[hp.iF3Z]
175            F3gZ_NNLO = res_nnlo[hp.iF3gZ]
176    
177            res_n3lo = hp.StrFctN3LO(xval, Qtest, muR, muF)
178            print(f"{res_n3lo[hp.iF3Wp]:11.4E} {res_n3lo[hp.iF3Wm]:11.4E}", end=' ')
179            F3Z_N3LO = res_n3lo[hp.iF3Z]
180            F3gZ_N3LO = res_n3lo[hp.iF3gZ]
181
182            print(f"{F3Z_LO:11.4E} {F3Z_NLO:11.4E} {F3Z_NNLO:11.4E} {F3Z_N3LO:11.4E} "
183                  f"{F3gZ_LO:11.4E} {F3gZ_NLO:11.4E} {F3gZ_NNLO:11.4E} {F3gZ_N3LO:11.4E}")
184        else:
185            res = hp.StrFct(xval, Qtest, muR, muF)
186            print(f"{xval:12.4E} {res[hp.iF3Wp]:11.4E} {res[hp.iF3Wm]:11.4E} {res[hp.iF3Z]:11.4E} {res[hp.iF3gZ]:11.4E}")
187    print()
188    print()
189
190if __name__ == "__main__":
191    main()

Reading LHAPDF into a tabulation and fast evaluation

  1#! /usr/bin/env python3
  2"""
  3This is a small example that uses the Hoppet and LHAPDF python
  4interfaces. It loads an LHAPDF set and assigns it to Hoppet.
  5"""
  6
  7import hoppet as hp
  8from hoppet import lhapdf
  9import lhapdf
 10import numpy as np
 11import argparse
 12import time
 13
 14def main():
 15    # Get commandline
 16    parser = argparse.ArgumentParser(description="Check of an LHAPDF grid against HOPPET evolution.")
 17    parser.add_argument('-pdf', type=str, default="PDF4LHC21_40", help='LHAPDF set name (default PDF4LHC21_40)')
 18    parser.add_argument('-yorder', type=int, default=2, help='order for interpolation in y=ln1/x (default=2)')
 19    parser.add_argument('-lnlnQorder', type=int, default=2, help='order for interpolation in lnlnQ (default=2)')
 20
 21    args = parser.parse_args()
 22
 23    hlinfo = hp.lhapdf.load(args.pdf)
 24    p_lhapdf = lhapdf.mkPDF(args.pdf, 0)
 25
 26    print("Checking that hlinfo is valid:", hlinfo is not None)
 27
 28    # Overwrite the yorder and lnlnQorder interpolation orders
 29    hp.SetYLnlnQInterpOrders(args.yorder, args.lnlnQorder)
 30
 31    # Evaluate the PDFs at some x values and print them
 32    xvals = [1e-5,1e-4,1e-3,1e-2,0.1,0.3,0.5,0.7,0.9]
 33    Q = 100.0
 34
 35    print("")
 36    print("           Evaluating PDFs at Q =",Q, " GeV")
 37    print("-------------- HOPPET ------------------")
 38    if not hlinfo.QED:
 39        print("    x      u-ubar      d-dbar    2(ubr+dbr)    c+cbar       gluon")
 40        for ix in range(9):
 41            pdf_array = hp.Eval(xvals[ix], Q)
 42            print("{:7.1E} {:11.4E} {:11.4E} {:11.4E} {:11.4E} {:11.4E}".format(
 43                xvals[ix],
 44                pdf_array[hp.iflv_u] - pdf_array[hp.iflv_ubar], 
 45                pdf_array[hp.iflv_d] - pdf_array[hp.iflv_dbar], 
 46                2 * (pdf_array[hp.iflv_dbar] + pdf_array[hp.iflv_ubar]),
 47                pdf_array[hp.iflv_cbar] + pdf_array[hp.iflv_c],
 48                pdf_array[hp.iflv_g]
 49            ))
 50        print("-------------- LHAPDF ------------------")
 51        for ix in range(9):
 52            pdf_dict = p_lhapdf.xfxQ(xvals[ix], Q)
 53            #print(type(pdf_dict)
 54            print("{:7.1E} {:11.4E} {:11.4E} {:11.4E} {:11.4E} {:11.4E}".format(
 55                xvals[ix],
 56                pdf_dict[2] - pdf_dict[-2], 
 57                pdf_dict[1] - pdf_dict[-1], 
 58                2 * (pdf_dict[-1] + pdf_dict[-2]),
 59                pdf_dict[-4] + pdf_dict[4],
 60                pdf_dict[21]
 61            ))
 62    else:
 63        print("    x      u-ubar      d-dbar    2(ubr+dbr)    c+cbar       gluon      photon")
 64        for ix in range(9):
 65            pdf_array = hp.Eval(xvals[ix], Q)
 66            print("{:7.1E} {:11.4E} {:11.4E} {:11.4E} {:11.4E} {:11.4E} {:11.4E}".format(
 67                xvals[ix],
 68                pdf_array[hp.iflv_u] - pdf_array[hp.iflv_ubar], 
 69                pdf_array[hp.iflv_d] - pdf_array[hp.iflv_dbar], 
 70                2 * (pdf_array[hp.iflv_dbar] + pdf_array[hp.iflv_ubar]),
 71                pdf_array[hp.iflv_cbar] + pdf_array[hp.iflv_c],
 72                pdf_array[hp.iflv_g], pdf_array[hp.iflv_photon]
 73            ))
 74        print("-------------- LHAPDF ------------------")
 75        for ix in range(9):
 76            pdf_dict = p_lhapdf.xfxQ(xvals[ix], Q)
 77            #print(type(pdf_dict)
 78            print("{:7.1E} {:11.4E} {:11.4E} {:11.4E} {:11.4E} {:11.4E} {:11.4E}".format(
 79                xvals[ix],
 80                pdf_dict[2] - pdf_dict[-2], 
 81                pdf_dict[1] - pdf_dict[-1], 
 82                2 * (pdf_dict[-1] + pdf_dict[-2]),
 83                pdf_dict[-4] + pdf_dict[4],
 84                pdf_dict[21], pdf_dict[22]
 85            ))
 86    
 87
 88
 89    #hp.WriteLHAPDFGrid("test_python",0)
 90
 91    # Define grids for timing test
 92    npoints = 1000
 93    xvals_timing = np.logspace(np.log10(1e-5), np.log10(0.9), npoints)
 94    Qvals_timing = np.logspace(np.log10(1.0), np.log10(1000.0), npoints)
 95
 96    # Timing HOPPET
 97    start_hoppet = time.perf_counter()
 98    for Q in Qvals_timing:
 99        for x in xvals_timing:
100            pdf_array = hp.Eval(x, Q)
101    end_hoppet = time.perf_counter()
102    print(f"HOPPET evaluation time {(end_hoppet - start_hoppet)/npoints/npoints*1e9:.2f} ns")
103
104    # Timing LHAPDF
105    # Load the PDF from LHAPDF
106    start_lhapdf = time.perf_counter()
107    for Q in Qvals_timing:
108        for x in xvals_timing:
109            pdf_dict = p_lhapdf.xfxQ(x, Q)
110            # The following call is faster but works only with a patch to LHAPDF
111            # pdf_list = p_lhapdf.xfxQv(x, Q)
112    end_lhapdf = time.perf_counter()
113    print(f"LHAPDF evaluation time {(end_lhapdf - start_lhapdf)/npoints/npoints*1e9:.2f} ns")
114
115    # AlphaS timings
116    start_hoppet_as = time.perf_counter()
117    for Q in Qvals_timing:
118        as_hoppet = hp.AlphaS(Q)
119    end_hoppet_as = time.perf_counter()
120    print(f"HOPPET alphaS time {(end_hoppet_as - start_hoppet_as)/npoints*1e9:.2f} ns") 
121    start_lhapdf_as = time.perf_counter()
122    for Q in Qvals_timing:
123        as_lhapdf = p_lhapdf.alphasQ(Q)
124    end_lhapdf_as = time.perf_counter()
125    print(f"LHAPDF alphaS time {(end_lhapdf_as - start_lhapdf_as)/npoints*1e9:.2f} ns")
126    
127    hp.DeleteAll()
128
129
130if __name__ == "__main__":
131    main()