Are 1D luminosity plots correct?

In the context of the PDF4LHC benchmarking exercise, I've realised that there is something that I don't understand with the luminosity plots generated with validphys. If I plot the 1D gg luminosity for sqrts=14 TeV I obtain different results w.r.t. those obtained with APFELweb. I attach the plots obtained with APFELweb and with validphys. lumi_APFEL.pdf lumi_validphys.pdf I think that the APFELweb result is correct because if I write a small python script by scratch I indeed obtain that the ordering of the ratios around mH is NNPDF31>MSHT20~CT18NNLO (as in the APFELweb plot) - while it's the other way around in the validphys plot. May this be a bug? Note that this is important, because all the luminosity plots in the NNPDF4.0 paper may be affected.

assigned to @enocera

assigned to @enocera

assigned to @enocera

assigned to @enocera

added bug label

changed title from Are luminosity plots correct? to Are 1D luminosity plots correct?

@scarrazza Would it be possible to make the NNPDF4.0 sets available in APFELweb? If I try to upload the NNPDF40_nnlo_as_01180 grid I receive an error message about disk space issues. Thanks.

Created by: Zaharid

Can you share the simple python script?

I think that my question is why the APFELweb and validphys output differ. That being said, here's the python script (I was asked to compute the lumi at mH, so it returns only one value). Note that the luminosity definition used differs from the standard one by a PDF-independent kinematic factor, which however cancels out when luminosity ratios are taken into account.

import lhapdf as lh
import scipy.integrate as integrate

pdfset = ["NNPDF40_nnlo_as_01180",
          "NNPDF31_nnlo_as_0118",
          "MSHT20nnlo_as118",
          "CT18NNLO"]

mH    = 125.5 #GeV
sqrts = 14000 # GeV
tau   = mH * mH / ( sqrts * sqrts )

x_min = tau
x_max = 1.0

lh.setVerbosity(0)

for iset in pdfset:
    lumi = 0.
    p = lh.mkPDF(iset,0)
    f = lambda x1: 1/x1 * p.xfxQ(21, x1, mH) * p.xfxQ(21, tau/x1, mH) 
    lumi = integrate.quad(f, x_min, x_max, epsrel=0.05, limit=10)[0]
    
    print(iset, lumi)

    if iset=="NNPDF40_nnlo_as_01180" or iset=="NNPDF31_nnlo_as_0118":
        lumi = 0.
        for irep in range(1,100,1):
            p = lh.mkPDF(iset,irep)
            lumi = lumi + integrate.quad(f, x_min, x_max, epsrel=0.05, limit=10)[0]
        print(iset,lumi/100)

And for completeness here is the vp runcard that I used to generate the luminosity plot

meta:
  title: Luminosity plot
  author: Emanuele R. Nocera
  keywords: NNPDF4.0

pdfs:
    - {id: NNPDF31_nnlo_as_0118, label: "NNPDF3.1 (NNLO)"}
    - {id: CT18NNLO, label: "CT18 (NNLO)"}
    - {id: MSHT20nnlo_as118, label: "MSHT20 (NNLO)"}

sqrts: 14000 # GeV
  
first:
  lumi_channel: "gg"
  ymin: 0.8
  ymax: 1.2

PDFscalespecs:
  - xscale: log

template_text: |
  {@with PDFscalespecs@}
  {@xscaletitle@} scale
  =====================
  {@first plot_lumi1d(normalize_to=1)@}
  {@endwith@}

actions_:
  - report(main=True)

Created by: Zaharid

If I change the script to have the same normalization as vp, I find that the results are the same as vp, i.e. the cteq lumi is higher

import lhapdf as lh
import scipy.integrate as integrate

pdfset = ["NNPDF40_nnlo_as_01180",
          "NNPDF31_nnlo_as_0118",
          #"MSHT20nnlo_as118",
          "CT18NNLO"]

mH    = 125.5 #GeV
sqrts = 14000 # GeV
tau   = mH * mH / ( sqrts * sqrts )

x_min = tau
x_max = 1.0

lh.setVerbosity(0)

for iset in pdfset:
    lumi = 0.
    p = lh.mkPDF(iset,0)
    f = lambda x1: 1/x1/(tau/x1)/sqrts/sqrts * p.xfxQ(21, x1, mH) * p.xfxQ(21, tau/x1, mH) 
    lumi = integrate.quad(f, x_min, x_max, epsrel=0.05, limit=10)[0]
    
    print(iset, lumi*1e4)

    if iset=="NNPDF40_nnlo_as_01180" or iset=="NNPDF31_nnlo_as_0118":
        lumi = 0.
        for irep in range(1,100,1):
            p = lh.mkPDF(iset,irep)
            lumi = lumi + integrate.quad(f, x_min, x_max, epsrel=0.05, limit=10)[0]
        print(iset,lumi/100*1e4)

Is this

f = lambda x1: 1/x1/(tau/x1)/sqrts/sqrts * p.xfxQ(21, x1, mH) * p.xfxQ(21, tau/x1, mH) 

correct? If I look at Eq.(1) of https://arxiv.org/pdf/1607.01831.pdf I would expect

f = lambda x1: 1/x1 * 1/sqrts/sqrts * 1/x1 * p.xfxQ(21, x1, mH) * 1/(tau/x1) * p.xfxQ(21, tau/x1, mH) 

Isn't there a factor 1/x1 missing in your (and validphys) definition of the luminosity?

@Zaharid I'm sorry to insist, but I believe that this https://github.com/NNPDF/nnpdf/blob/622fb9c6b60276d2e98986212f49797fb6d59df3/validphys2/src/validphys/gridvalues.py#L170 should read

return res/x1/x2/s/x1

If I do this, I indeed obtain the attached plot, which is similar to the APFELweb one (and the correct one, I believe) modulo improving the smoothness of the integration. gglumi_corrected.pdf

Created by: Zaharid

Appears to be wrong indeed...

Just to be sure: are you saying that you agree with me? If so, let me open a PR.

Created by: Zaharid

The lumi 1d values in vp look wrong (they are not symmetric under x1 -> x2, i ->j and also not what that paper says). But I get the impression that that 2d version of evaluate_luminosity should be symmetric.

Created by: Zaharid

Incidentally I am pretty sure the current version was "fixed" somehow from a previous state, but can't find the discussion.

The paper I'm referring to (https://arxiv.org/pdf/1607.01831.pdf) is the one we cite in the NNPDF3.1 paper. In the NNPDF4.0 paper we don't cite any reference for the luminosity expressions.

BTW, I think that Eq.(2) in https://arxiv.org/pdf/1607.01831.pdf (the 2D luminosity) makes sense. I also believe that Eq.(3) (the 2D luminosity integrated) makes sense as well. And Eq.(3) is equivalent to Eq.(1) (after working out some simple algebra).

And I'd say that Eq.(1) in https://arxiv.org/pdf/1607.01831.pdf is the same as Eq.(9.1) in the NNPDF4.0 paper.

mentioned in merge request !1443 (merged)

closed