Replica number problem in PDFset for PDF uncertainty calculation

Hi,

My first suggestions would be to use the offline method to compute such uncertainties instead of the online method:
In the run_card, you can set "store_rwgt_info" on "True" and do not ask for such PDF reweighing.
Then when your job is finished, you can do the following:
./bin/aMCatNLO
Then type
systematics run_01 --pdf=XXXXX
For more help on the valid options associate to this command, you can type "help systematics"

The second options is to keep to do it with the online method, but then you need to hack the code to increase the number of supported PDF.
You can just edit the file SubProcesses/reweight0.inc and change the value of
c Cross sections after variations of scales and PDFs
c
      integer maxscales,maxPDFs
      parameter (maxscales=9)
      parameter (maxPDFs=200)
To the number that you need.

Cheers,

Olivier

> On 6 Jun 2017, at 07:17, Sung Hak Lim <email address hidden> wrote:
>
> New question #641514 on MadGraph5_aMC@NLO:
> https://answers.launchpad.net/mg5amcnlo/+question/641514
>
> Dear MG5 team.
>
> May I know current MG5 version is compatible to NNPDFset with larger replica, such as NNPDF30_nnlo_as_0118_1000?
>
> I just got following errors.
>
> tail: 읽기를 위해 '/data/starlight/mgoutput/ttbar_nlo/SubProcesses/P0_gg_ttx/GF3/events.lhe.rwgt'을(를) 열 수 없음: 그런 파일이나 디렉터리가 없습니다 << it says that events.lhe.rwgt is not found.
> Error detected in "launch auto "
> write debug file /data/starlight/mgoutput/ttbar_nlo/run_20_tag_1_debug.log
> If you need help with this issue please contact us on https://answers.launchpad.net/mg5amcnlo
> aMCatNLOError : An error occurred during reweight. Check the'reweight_xsec_events.output' files inside the 'SubProcesses/P*/G*/ directories for details
>
> I also checked that this error disappears when we set PDFset member number in NNPDF30_nnlo_as_0118_1000.info below 200. However, it reappears for larger PDFset member number 201.
>
> Thanks for further help and have a nice day!
>
> Best wishes,
> Sung Hak Lim
>
> ==========================
>
> MG5 version: 2.5.4
> Process: p p > t t~ [QCD]
> Used pdfset lhaid: 230000 295000
>
>
> reweight_xsec_events.output:
>
> ==== LHAPDF6 USING DEFAULT-TYPE LHAGLUE INTERFACE ====
> LHAPDF 6.1.6 loading /home/starlight/local/share/LHAPDF/NNPDF23_nlo_as_0119/NNPDF23_nlo_as_0119_0000.dat
> NNPDF23_nlo_as_0119 PDF set, member #0, version 6; LHAPDF ID = 230000
> LHAPDF 6.1.6 loading /home/starlight/local/share/LHAPDF/NNPDF30_nnlo_as_0118_1000/NNPDF30_nnlo_as_0118_1000_0000.dat
> NNPDF30_nnlo_as_0118_1000 PDF set, member #0, version 1; LHAPDF ID = 295000
> Thanks for using LHAPDF 6.1.6. Please make sure to cite the paper:
> Eur.Phys.J. C75 (2015) 3, 132 (http://arxiv.org/abs/1412.7420)
> A PDF is used, so alpha_s(MZ) is going to be modified
> Old value of alpha_s from param_card: 0.11799999999999999
> New value of alpha_s from PDF lhapdf : 0.11900242160818149
> using LHAPDF
> Enter event file name
> Enter 1 to save all cross sections on tape
> 0 otherwise
> Including scale uncertainties for dynamical_scale_choice -1
> - renormalisation scale variation: 1.0000000000000000 2.0000000000000000 0.50000000000000000
> - factorisation scale variation: 1.0000000000000000 2.0000000000000000 0.50000000000000000
> Including central PDF with uncertainties for NNPDF23_nlo_as_0119
> Including central PDF with uncertainties for NNPDF30_nnlo_as_0118_1000 \
> Too many PDFs: increase maxPDFs in reweight0.inc to 1
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hi Olivier,

Let me try.

Thanks for answers and let me know you the result.

Best wishes,
Sung Hak

Hi Olivier,

I just tried the first method, and it's working well.

However, when I use the same PDF for the online and offline method for a sanity check, the pdf variations that I've got are different.

Here are the outputs I've got:

=================
Online
=================
   --------------------------------------------------------------
      Summary:
      Process p p > t t~ [QCD]
      Run at p-p collider (6500.0 + 6500.0 GeV)
      Number of events generated: 10000
      Total cross section: 6.935e+02 +- 3.9e+00 pb
   --------------------------------------------------------------
      Scale variation (computed from LHE events):
          Dynamical_scale_choice -1 (envelope of 9 values):
              6.884e+02 pb +9.0% -10.2%
      PDF variation (computed from LHE events):
          NNPDF23_nlo_as_0119 (101 members; using replicas method):
              6.884e+02 pb +1.5% -1.5%
   --------------------------------------------------------------

=======================
Offline
=======================
INFO: # events generated with PDF: NNPDF23_nlo_as_0119 (230000)
INFO: #Will Compute 145 weights per event.
INFO: #***************************************************************************
#
# original cross-section: 687.654961899
# scale variation: +-45.5% -55.1%
# central scheme variation: + 0% -61.8%
# PDF variation: +0.753% -0.753%
#
# dynamical scheme # 1 : 309.923 +6.75% -9.27% # \sum ET
# dynamical scheme # 2 : 308.619 +11.4% -11.3% # \sum\sqrt{m^2+pt^2}
# dynamical scheme # 3 : 343.827 +9.02% -10.2% # 0.5 \sum\sqrt{m^2+pt^2}
# dynamical scheme # 4 : 262.73 +6.46% -8.6% # \sqrt{\hat s}
#***************************************************************************

May I know what's happening behind this difference?

Also in offline method, dynamical schemes seem to giving too different cross section compared to original cross section....

Best wishes,
Sung Hak

Hi,

Your offline seems indeed completly wrong for a NLO accurate sample.
your scale variation in particular is clearly wrong. Did you correctly set the flag at generation time to store the details need for the computation? and/or is this run a LO run?

In my case, I have perfect agreement between those two methods at NLO. (for the exact same process/ PDF set):
INFO: # events generated with PDF: NNPDF23_nlo_as_0119 (230000)
INFO: #Will Compute 145 weights per event.
INFO: #***************************************************************************
#
# original cross-section: 683.560892025
# scale variation: +9.6% -10.6%
# central scheme variation: + 0% -10.6%
# PDF variation: +1.46% -1.46%
#
# dynamical scheme # 1 : 620.883 +8.05% -9.83% # \sum ET
# dynamical scheme # 2 : 611.377 +11.8% -11.5% # \sum\sqrt{m^2+pt^2}
# dynamical scheme # 3 : 683.561 +9.6% -10.6% # 0.5 \sum\sqrt{m^2+pt^2}
# dynamical scheme # 4 : 615.592 +13.6% -12.4% # \sqrt{\hat s}
#***************************************************************************

Cheers,

Olivier

Hi,

I also agree that something goes wrong in my offline case.

It's NLO run and let me reconfigure my MG5 system again.

Best wishes,
Sung Hak

Hi,

I tested this again with independent desktop and installing all pacakges using install command in MadGraph itself.

I also tested various things, changing gcc, check other process (for example, p p > z z [QCD]) , reinstall lhapdf, testing in other collegue's desktop configuration

Still, when I run systematics run in offline mode, I got all weights 1/2 than my expectation.

Here I attached my run_card when doing event generation with store_rwgt_info = True.

#***********************************************************************
# MadGraph5_aMC@NLO *
# *
# run_card.dat aMC@NLO *
# *
# This file is used to set the parameters of the run. *
# *
# Some notation/conventions: *
# *
# Lines starting with a hash (#) are info or comments *
# *
# mind the format: value = variable ! comment *
# *
# Some of the values of variables can be list. These can either be *
# comma or space separated. *
#***********************************************************************
#
#*******************
# Running parameters
#*******************
#
#***********************************************************************
# Tag name for the run (one word) *
#***********************************************************************
  tag_1 = run_tag ! name of the run
#***********************************************************************
# Number of LHE events (and their normalization) and the required *
# (relative) accuracy on the Xsec. *
# These values are ignored for fixed order runs *
#***********************************************************************
  10000 = nevents ! Number of unweighted events requested
  -1.0 = req_acc ! Required accuracy (-1=auto determined from nevents)
  -1 = nevt_job ! Max number of events per job in event generation.
                 ! (-1= no split).
#***********************************************************************
# Normalize the weights of LHE events such that they sum or average to *
# the total cross section *
#***********************************************************************
  average = event_norm ! average or sum
#***********************************************************************
# Number of points per itegration channel (ignored for aMC@NLO runs) *
#***********************************************************************
  0.01 = req_acc_fo ! Required accuracy (-1=ignored, and use the
                     ! number of points and iter. below)
# These numbers are ignored except if req_acc_FO is equal to -1
  5000 = npoints_fo_grid ! number of points to setup grids
  4 = niters_fo_grid ! number of iter. to setup grids
  10000 = npoints_fo ! number of points to compute Xsec
  6 = niters_fo ! number of iter. to compute Xsec
#***********************************************************************
# Random number seed *
#***********************************************************************
  0 = iseed ! rnd seed (0=assigned automatically=default))
#***********************************************************************
# Collider type and energy *
#***********************************************************************
  1 = lpp1 ! beam 1 type (0 = no PDF)
  1 = lpp2 ! beam 2 type (0 = no PDF)
  6500.0 = ebeam1 ! beam 1 energy in GeV
  6500.0 = ebeam2 ! beam 2 energy in GeV
#***********************************************************************
# PDF choice: this automatically fixes also alpha_s(MZ) and its evol. *
#***********************************************************************
  lhapdf = pdlabel ! PDF set
  230000 = lhaid ! If pdlabel=lhapdf, this is the lhapdf number. Only
              ! numbers for central PDF sets are allowed. Can be a list;
              ! PDF sets beyond the first are included via reweighting.
#***********************************************************************
# Include the NLO Monte Carlo subtr. terms for the following parton *
# shower (HERWIG6 | HERWIGPP | PYTHIA6Q | PYTHIA6PT | PYTHIA8) *
# WARNING: PYTHIA6PT works only for processes without FSR!!!! *
#***********************************************************************
  PYTHIA8 = parton_shower
  1.0 = shower_scale_factor ! multiply default shower starting
                                  ! scale by this factor
#***********************************************************************
# Renormalization and factorization scales *
# (Default functional form for the non-fixed scales is the sum of *
# the transverse masses divided by two of all final state particles *
# and partons. This can be changed in SubProcesses/set_scales.f or via *
# dynamical_scale_choice option) *
#***********************************************************************
  False = fixed_ren_scale ! if .true. use fixed ren scale
  False = fixed_fac_scale ! if .true. use fixed fac scale
  91.118 = mur_ref_fixed ! fixed ren reference scale
  91.118 = muf_ref_fixed ! fixed fact reference scale
  -1 = dynamical_scale_choice ! Choose one (or more) of the predefined
           ! dynamical choices. Can be a list; scale choices beyond the
           ! first are included via reweighting
  1.0 = mur_over_ref ! ratio of current muR over reference muR
  1.0 = muf_over_ref ! ratio of current muF over reference muF
#***********************************************************************
# Reweight variables for scale dependence and PDF uncertainty *
#***********************************************************************
  1.0, 2.0, 0.5 = rw_rscale ! muR factors to be included by reweighting
  1.0, 2.0, 0.5 = rw_fscale ! muF factors to be included by reweighting
  True = reweight_scale ! Reweight to get scale variation using the
            ! rw_rscale and rw_fscale factors. Should be a list of
            ! booleans of equal length to dynamical_scale_choice to
            ! specify for which choice to include scale dependence.
  False = reweight_pdf ! Reweight to get PDF uncertainty. Should be a
            ! list booleans of equal length to lhaid to specify for
            ! which PDF set to include the uncertainties.
#***********************************************************************
# Store reweight information in the LHE file for off-line model- *
# parameter reweighting at NLO+PS accuracy *
#***********************************************************************
  True = store_rwgt_info ! Store info for reweighting in LHE file
#***********************************************************************
# ickkw parameter: *
# 0: No merging *
# 3: FxFx Merging - WARNING! Applies merging only at the hard-event *
# level. After showering an MLM-type merging should be applied as *
# well. See http://amcatnlo.cern.ch/FxFx_merging.htm for details. *
# 4: UNLOPS merging (with pythia8 only). No interface from within *
# MG5_aMC available, but available in Pythia8. *
# -1: NNLL+NLO jet-veto computation. See arxiv:1412.8408 [hep-ph]. *
#***********************************************************************
  0 = ickkw
#***********************************************************************
#
#***********************************************************************
# BW cutoff (M+/-bwcutoff*Gamma). Determines which resonances are *
# written in the LHE event file *
#***********************************************************************
  15.0 = bwcutoff
#***********************************************************************
# Cuts on the jets. Jet clustering is performed by FastJet. *
# - When matching to a parton shower, these generation cuts should be *
# considerably softer than the analysis cuts. *
# - More specific cuts can be specified in SubProcesses/cuts.f *
#***********************************************************************
  1.0 = jetalgo ! FastJet jet algorithm (1=kT, 0=C/A, -1=anti-kT)
  0.7 = jetradius ! The radius parameter for the jet algorithm
  10.0 = ptj ! Min jet transverse momentum
  -1.0 = etaj ! Max jet abs(pseudo-rap) (a value .lt.0 means no cut)
#***********************************************************************
# Cuts on the charged leptons (e+, e-, mu+, mu-, tau+ and tau-) *
# More specific cuts can be specified in SubProcesses/cuts.f *
#***********************************************************************
  0.0 = ptl ! Min lepton transverse momentum
  -1.0 = etal ! Max lepton abs(pseudo-rap) (a value .lt.0 means no cut)
  0.0 = drll ! Min distance between opposite sign lepton pairs
  0.0 = drll_sf ! Min distance between opp. sign same-flavor lepton pairs
  0.0 = mll ! Min inv. mass of all opposite sign lepton pairs
  30.0 = mll_sf ! Min inv. mass of all opp. sign same-flavor lepton pairs
#***********************************************************************
# Photon-isolation cuts, according to hep-ph/9801442. When ptgmin=0, *
# all the other parameters are ignored. *
# More specific cuts can be specified in SubProcesses/cuts.f *
#***********************************************************************
  20.0 = ptgmin ! Min photon transverse momentum
  -1.0 = etagamma ! Max photon abs(pseudo-rap)
  0.4 = r0gamma ! Radius of isolation code
  1.0 = xn ! n parameter of eq.(3.4) in hep-ph/9801442
  1.0 = epsgamma ! epsilon_gamma parameter of eq.(3.4) in hep-ph/9801442
  True = isoem ! isolate photons from EM energy (photons and leptons)
#***********************************************************************
# For aMCfast+APPLGRID use in PDF fitting (http://amcfast.hepforge.org)*
#***********************************************************************
  0 = iappl ! aMCfast switch (0=OFF, 1=prepare grids, 2=fill grids)
#***********************************************************************

I'm still trying to keep solving this now.

Best wishes,
Sung Hak

Hi,

I run the following script (./bin/mg5_aMC PATH_TO_FILE)
With the file containing the following lines (which should correspond to the run_card that you send me):
shower=OFF
set pdlabel lhapdf
set lhaid 230000
set store_rwgt_info T
set parton_shower PYTHIA8
set mll_sf 30
set ptgmin 20
set isoem T
launch -i
systematics run_01
exit
launch
set store_wgt_info F
set reweight_pdf True

You will see the results below both for 2.5.5 and 2.5.4 version, but in both case, I do not reproduce any problem.
Since the one which is wrong in your case is the offline method, changing gcc is irrelevant. Indeed that code is pure
Python. (With the python interface to lhapdf which is c++). So my guess is that your problem is related to the lhapdf.
In my case, lhapdf version is 6.1.6 and was installed via "install lhapdf" of madgrah.
What is your version of lhapdf?

Cheers,

Olivier

And here is the result that I get for 2.5.5

1) first run without pdf systematics
   --------------------------------------------------------------
      Summary:
      Process p p > z z [QCD]
      Run at p-p collider (6500.0 + 6500.0 GeV)
      Number of events generated: 10000
      Total cross section: 1.375e+01 +- 6.4e-02 pb
   --------------------------------------------------------------
      Scale variation (computed from LHE events):
          Dynamical_scale_choice -1 (envelope of 9 values):
              1.374e+01 pb +3.3% -3.5%
   --------------------------------------------------------------

2) evaluation of the systematics
INFO: #***************************************************************************
#
# original cross-section: 13.7318660283
# scale variation: +3.27% -3.5%
# central scheme variation: +1.3% -0.991%
# PDF variation: +1.18% -1.18%
#
# dynamical scheme # 1 : 13.5958 +2.66% -2.93% # \sum ET
# dynamical scheme # 2 : 13.7051 +2.87% -2.91% # \sum\sqrt{m^2+pt^2}
# dynamical scheme # 3 : 13.7319 +3.27% -3.5% # 0.5 \sum\sqrt{m^2+pt^2}
# dynamical scheme # 4 : 13.9103 +3.29% -3.2% # \sqrt{\hat s}
#***************************************************************************

3) second run with pdf systematics
   --------------------------------------------------------------
      Summary:
      Process p p > z z [QCD]
      Run at p-p collider (6500.0 + 6500.0 GeV)
      Number of events generated: 10000
      Total cross section: 1.368e+01 +- 5.2e-02 pb
   --------------------------------------------------------------
      Scale variation (computed from LHE events):
          Dynamical_scale_choice -1 (envelope of 9 values):
              1.383e+01 pb +3.7% -3.9%
      PDF variation (computed from LHE events):
          NNPDF23_nlo_as_0119 (101 members; using replicas method):
              1.384e+01 pb +1.2% -1.2%
   --------------------------------------------------------------

And here is the result that I get for 2.5.4:

      Summary:
      Process p p > z z [QCD]
      Run at p-p collider (6500.0 + 6500.0 GeV)
      Number of events generated: 10000
      Total cross section: 1.375e+01 +- 6.4e-02 pb
   --------------------------------------------------------------
      Scale variation (computed from LHE events):
          Dynamical_scale_choice -1 (envelope of 9 values):
              1.374e+01 pb +3.3% -3.5%
   --------------------------------------------------------------

#***************************************************************************
#
# original cross-section: 13.7318660283
# scale variation: +3.27% -3.5%
# central scheme variation: +0.00857% -1.84%
# PDF variation: +1.18% -1.18%
#
# dynamical scheme # 1 : 13.5511 +2.56% -2.86% # \sum ET
# dynamical scheme # 2 : 13.7061 +2.88% -2.92% # \sum\sqrt{m^2+pt^2}
# dynamical scheme # 3 : 13.733 +3.28% -3.51% # 0.5 \sum\sqrt{m^2+pt^2}
# dynamical scheme # 4 : 13.4793 +2.15% -2.42% # \sqrt{\hat s}
#***************************************************************************

On 12 Jun 2017, at 10:36, Sung Hak Lim <<email address hidden><mailto:<email address hidden>>> wrote:

store_rwgt_info

Hi Olivier,

I also just noticed that ./bin/internal/systematics.py runs within python.

I also agree that the error could comes in is lhapdf part in python.

I'm using LHAPDF 6.1.6 now.

I tested all cases with MG5 2.5.5 and MG5 2.5.4, installing lhapdf by MG5 and by myself. but all give same errors.

Now I'm looking for chance of my python misconfiguration. Let's see..

Best wishes,
Sung Hak