Weird results from virtual reweighter

Asked by Andrin Kessler

Hi,

I am using the virt_reweighter python script to reweigh some Z boson production events. I largely get the expected results but every run there are a few outliers.

The process p p > Z > mu+ mu- [QCD], for instance, mostly gives the correct numbers but every once in a while there is a weight that is about 10^2 or 10^3 times too large (also of negative sign sometimes). This results in an incorrect ratio of the reweighted quantity over the Born weight (should analytically be ~2.49). The events where this happens don't seem special to my, admittedly untrained, eyes.
I am using a dynamical scale of mu = Q and massless quarks (ickkw = -1). To do this I have brute forced this choice in the setscales.f script in a somewhat inelegant way but judging from the event files this part of the code seems to work for every event.

On a side note, the process p p > Z [QCD] also contains a few zero value weights after the reweighting script is run which I don't understand but am less concerned about at the moment.

I include the slightly altered function in setscales.f, the run card I use, and a few examples of events where the reweighting worked and a few where it failed.

Thank you for your help!
Andrin

Within setscales.f:

      function scale_global_reference(pp)
c This is a function of the kinematic configuration pp, which returns
c a scale to be used as a reference for renormalization scale
      implicit none
      include 'genps.inc'
      include 'nexternal.inc'
      include 'run.inc'
      include 'cuts.inc'
      double precision ppv(0:3)
      double precision scale_global_reference,pp(0:3,nexternal)
      double precision tmp,pt,et,dot,xm2,sumdot,xmt2,ptmp(0:3)
      external pt,et,dot,sumdot
      integer i,j
      character*80 temp_scale_id
      common/ctemp_scale_id/temp_scale_id
c
      tmp=0
      if(ickkw.eq.-1) dynamical_scale_choice=10
      if(.false.)then
c Special for analytic resummation in veto'ed cross sections:
         tmp=ptj
         temp_scale_id='NLO+NNLL veto scale: ptj_max'
      elseif(dynamical_scale_choice.eq.1) then
c Total transverse energy of the event.
          tmp=0d0
          do i=3,nexternal
             tmp=tmp+et(pp(0,i))
          enddo
          temp_scale_id='sum_i eT(i), i=final state'
      elseif(dynamical_scale_choice.eq.2) then
c sum of the transverse mass divide
c m^2+pt^2=p(0)^2-p(3)^2=(p(0)+p(3))*(p(0)-p(3))
          tmp=0d0
          do i=3,nexternal
            tmp=tmp+dsqrt(max(0d0,(pp(0,i)+pp(3,i))*(pp(0,i)-pp(3,i))))
          enddo
          temp_scale_id='sum_i mT(i), i=final state'
      elseif(dynamical_scale_choice.eq.3.or.dynamical_scale_choice.eq.-1) then
c sum of the transverse mass divide by 2
c m^2+pt^2=p(0)^2-p(3)^2=(p(0)+p(3))*(p(0)-p(3))
          tmp=0d0
          do i=3,nexternal
            tmp=tmp+dsqrt(max(0d0,(pp(0,i)+pp(3,i))*(pp(0,i)-pp(3,i))))
          enddo
          tmp=tmp/2d0
          temp_scale_id='H_T/2 := sum_i mT(i)/2, i=final state'
      elseif(dynamical_scale_choice.eq.0) then
c fixed scale
          tmp=muR_ref_fixed
          temp_scale_id='fixed scale'
      elseif(dynamical_scale_choice.eq.10) then
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc USER-DEFINED SCALE: ENTER YOUR CODE HERE cc
cc to use this code you must set cc
cc dynamical_scale_choice = 10 cc
cc in the run_card (run_card.dat) cc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        temp_scale_id='User-defined dynamical scale: mu = Q' ! use a meaningful string
        tmp = 0d0
        do i = 0,3
          ppv(i) = 0.0
          do j = 3, nexternal-1
            ppv(i) = ppv(i) + pp(i,j)
          enddo
        enddo

        tmp = dsqrt(dot(ppv,ppv))
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cc USER-DEFINED SCALE: END OF USER CODE cc
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      else
        write(*,*)'Unknown option in scale_global_reference',dynamical_scale_choice
        stop
      endif

      scale_global_reference=tmp

Run card:

<LesHouchesEvents version="3.0">
  <header>
  <MG5ProcCard>
#************************************************************
#* MadGraph5_aMC@NLO *
#* *
#* * * *
#* * * * * *
#* * * * * 5 * * * * *
#* * * * * *
#* * * *
#* *
#* *
#* VERSION 2.7.0 2020-01-20 *
#* *
#* The MadGraph5_aMC@NLO Development Team - Find us at *
#* https://server06.fynu.ucl.ac.be/projects/madgraph *
#* *
#************************************************************
#* *
#* Command File for MadGraph5_aMC@NLO *
#* *
#* run as ./bin/mg5_aMC filename *
#* *
#************************************************************
set default_unset_couplings 99
set group_subprocesses Auto
set ignore_six_quark_processes False
set loop_optimized_output True
set loop_color_flows False
set gauge unitary
set complex_mass_scheme False
set max_npoint_for_channel 0
import model sm
define p = g u c d s u~ c~ d~ s~
define j = g u c d s u~ c~ d~ s~
define l+ = e+ mu+
define l- = e- mu-
define vl = ve vm vt
define vl~ = ve~ vm~ vt~
import model loop_sm-no_b_mass
define p = 21 2 4 1 3 -2 -4 -1 -3 5 -5 # pass to 5 flavors
define j = p
generate p p > Z > mu+ mu- [QCD]
output ppzmupmum2
  </MG5ProcCard>
  <slha>
######################################################################
## PARAM_CARD AUTOMATICALY GENERATED BY MG5 FOLLOWING UFO MODEL ####
######################################################################
## ##
## Width set on Auto will be computed following the information ##
## present in the decay.py files of the model. ##
## See arXiv:1402.1178 for more details. ##
## ##
######################################################################

###################################
## INFORMATION FOR LOOP
###################################
Block loop
    1 9.118800e+01 # MU_R

###################################
## INFORMATION FOR MASS
###################################
Block mass
    6 1.730000e+02 # MT
   15 1.777000e+00 # MTA
   23 9.118800e+01 # MZ
   25 1.250000e+02 # MH
## Dependent parameters, given by model restrictions.
## Those values should be edited following the
## analytical expression. MG5 ignores those values
## but they are important for interfacing the output of MG5
## to external program such as Pythia.
  1 0.000000e+00 # d : 0.0
  2 0.000000e+00 # u : 0.0
  3 0.000000e+00 # s : 0.0
  4 0.000000e+00 # c : 0.0
  5 0.000000e+00 # b : 0.0
  11 0.000000e+00 # e- : 0.0
  12 0.000000e+00 # ve : 0.0
  13 0.000000e+00 # mu- : 0.0
  14 0.000000e+00 # vm : 0.0
  16 0.000000e+00 # vt : 0.0
  21 0.000000e+00 # g : 0.0
  22 0.000000e+00 # a : 0.0
  24 8.041900e+01 # w+ : cmath.sqrt(MZ__exp__2/2. + cmath.sqrt(MZ__exp__4/4. - (aEW*cmath.pi*MZ__exp__2)/(Gf*sqrt__2)))

###################################
## INFORMATION FOR SMINPUTS
###################################
Block sminputs
    1 1.325070e+02 # aEWM1
    2 1.166390e-05 # Gf
    3 1.180000e-01 # aS

###################################
## INFORMATION FOR YUKAWA
###################################
Block yukawa
    6 1.730000e+02 # ymt
   15 1.777000e+00 # ymtau

###################################
## INFORMATION FOR DECAY
###################################
DECAY 6 1.491500e+00 # WT
DECAY 23 2.441404e+00 # WZ
DECAY 24 2.047600e+00 # WW
DECAY 25 6.382339e-03 # WH
## Dependent parameters, given by model restrictions.
## Those values should be edited following the
## analytical expression. MG5 ignores those values
## but they are important for interfacing the output of MG5
## to external program such as Pythia.
DECAY 1 0.000000e+00 # d : 0.0
DECAY 2 0.000000e+00 # u : 0.0
DECAY 3 0.000000e+00 # s : 0.0
DECAY 4 0.000000e+00 # c : 0.0
DECAY 5 0.000000e+00 # b : 0.0
DECAY 11 0.000000e+00 # e- : 0.0
DECAY 12 0.000000e+00 # ve : 0.0
DECAY 13 0.000000e+00 # mu- : 0.0
DECAY 14 0.000000e+00 # vm : 0.0
DECAY 15 0.000000e+00 # ta- : 0.0
DECAY 16 0.000000e+00 # vt : 0.0
DECAY 21 0.000000e+00 # g : 0.0
DECAY 22 0.000000e+00 # a : 0.0
#===========================================================
# QUANTUM NUMBERS OF NEW STATE(S) (NON SM PDG CODE)
#===========================================================

Block QNUMBERS 82 # gh
        1 0 # 3 times electric charge
        2 1 # number of spin states (2S+1)
        3 8 # colour rep (1: singlet, 3: triplet, 8: octet)
        4 1 # Particle/Antiparticle distinction (0=own anti)
  </slha>
  <MGRunCard>
#***********************************************************************
# 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. *
# *
# To display additional parameter, you can use the command: *
# update to_full *
#***********************************************************************
#
#*******************
# 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
 -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 ! valid settings: average, sum, bias
#***********************************************************************
# 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 *
#***********************************************************************
         33 = iseed
#***********************************************************************
# 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
 25200 = 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!!!! *
#***********************************************************************
  HERWIG6 = 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
 10 = 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 *
#***********************************************************************
 False = 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]. *
#***********************************************************************
 -1 = ickkw
#***********************************************************************
#
#***********************************************************************
# BW cutoff (M+/-bwcutoff*Gamma). Determines which resonances are *
# written in the LHE event file *
#***********************************************************************
 1000000.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)
#***********************************************************************
# Cuts associated to MASSIVE particles identified by their PDG codes. *
# All cuts are applied to both particles and anti-particles, so use *
# POSITIVE PDG CODES only. Example of the syntax is {6 : 100} or *
# {6:100, 25:200} for multiple particles *
#***********************************************************************
  {} = pt_min_pdg ! Min pT for a massive particle
  {} = pt_max_pdg ! Max pT for a massive particle
  {} = mxx_min_pdg ! inv. mass for any pair of (anti)particles
#***********************************************************************
# For aMCfast+APPLGRID use in PDF fitting (http://amcfast.hepforge.org)*
#***********************************************************************
 0 = iappl ! aMCfast switch (0=OFF, 1=prepare grids, 2=fill grids)
#***********************************************************************
  </MGRunCard>
  <scalesfunctionalform>
    muR User-defined dynamical scale
    muF1 User-defined dynamical scale
    muF2 User-defined dynamical scale
    QES User-defined dynamical scale
  </scalesfunctionalform>
  <MonteCarloMasses>
       1 0.320000E+00
       2 0.320000E+00
       3 0.500000E+00
       4 0.155000E+01
       5 0.495000E+01
      11 0.510999E-03
      13 0.105658E+00
      15 0.177682E+01
      21 0.750000E+00
  </MonteCarloMasses>
  <initrwgt>
    <weightgroup name='scale_variation -1' combine='envelope'>
      <weight id='1001'> dyn= -1 muR=0.10000E+01 muF=0.10000E+01 </weight>
      <weight id='1002'> dyn= -1 muR=0.20000E+01 muF=0.10000E+01 </weight>
      <weight id='1003'> dyn= -1 muR=0.50000E+00 muF=0.10000E+01 </weight>
      <weight id='1004'> dyn= -1 muR=0.10000E+01 muF=0.20000E+01 </weight>
      <weight id='1005'> dyn= -1 muR=0.20000E+01 muF=0.20000E+01 </weight>
      <weight id='1006'> dyn= -1 muR=0.50000E+00 muF=0.20000E+01 </weight>
      <weight id='1007'> dyn= -1 muR=0.10000E+01 muF=0.50000E+00 </weight>
      <weight id='1008'> dyn= -1 muR=0.20000E+01 muF=0.50000E+00 </weight>
      <weight id='1009'> dyn= -1 muR=0.50000E+00 muF=0.50000E+00 </weight>
    </weightgroup>
  </initrwgt>
  </header>
  <init>
   2212 2212 0.65000000E+04 0.65000000E+04 -1 -1 25200 25200 -4 1
 0.16394779E+04 0.58588773E+01 0.16394779E+04 0
  </init>

Two working events:

<event>
 5 0 +4.0869030e+03 8.92889780e+01 7.54677110e-03 1.18377930e-01
        1 -1 0 0 501 0 +0.0000000000e+00 +0.0000000000e+00 +2.2288570000e+02 2.2288570000e+02 0.0000000000e+00 0.0000e+00 9.0000e+00
       -1 -1 0 0 0 501 +0.0000000000e+00 +0.0000000000e+00 -8.9423879000e+00 8.9423879000e+00 0.0000000000e+00 0.0000e+00 9.0000e+00
       23 2 1 2 0 0 +0.0000000000e+00 +0.0000000000e+00 +2.1394331000e+02 2.3182809000e+02 8.9288978000e+01 0.0000e+00 0.0000e+00
      -13 1 3 3 0 0 -3.2064456000e+01 -1.5905795000e+01 +1.7625204000e+02 1.7984967000e+02 0.0000000000e+00 0.0000e+00 9.0000e+00
       13 1 3 3 0 0 +3.2064456000e+01 +1.5905795000e+01 +3.7691277000e+01 5.1978418000e+01 0.0000000000e+00 0.0000e+00 9.0000e+00
#aMCatNLO 1 0 0 0 0 0.00000000E+00 0.00000000E+00 9 5 0 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00
<rwgt>
<wgt id='1001'> 0.16394E+04 </wgt>
<wgt id='1002'> 0.16394E+04 </wgt>
<wgt id='1003'> 0.16394E+04 </wgt>
<wgt id='1004'> 0.18213E+04 </wgt>
<wgt id='1005'> 0.18213E+04 </wgt>
<wgt id='1006'> 0.18213E+04 </wgt>
<wgt id='1007'> 0.14413E+04 </wgt>
<wgt id='1008'> 0.14413E+04 </wgt>
<wgt id='1009'> 0.14413E+04 </wgt>
</rwgt>
</event>
<event>
 5 0 +4.0869040e+03 9.03552890e+01 7.54677110e-03 1.18164620e-01
        1 -1 0 0 501 0 +0.0000000000e+00 +0.0000000000e+00 +1.2012059000e+02 1.2012059000e+02 0.0000000000e+00 0.0000e+00 9.0000e+00
       -1 -1 0 0 0 501 +0.0000000000e+00 +0.0000000000e+00 -1.6991421000e+01 1.6991421000e+01 0.0000000000e+00 0.0000e+00 9.0000e+00
       23 2 1 2 0 0 +0.0000000000e+00 +0.0000000000e+00 +1.0312917000e+02 1.3711201000e+02 9.0355289000e+01 0.0000e+00 0.0000e+00
      -13 1 3 3 0 0 +7.2441661000e+00 -3.8988041000e+01 +8.4409249000e+01 9.3260209000e+01 0.0000000000e+00 0.0000e+00 9.0000e+00
       13 1 3 3 0 0 -7.2441661000e+00 +3.8988041000e+01 +1.8719923000e+01 4.3851805000e+01 0.0000000000e+00 0.0000e+00 9.0000e+00
#aMCatNLO 1 0 0 0 0 0.00000000E+00 0.00000000E+00 9 5 0 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00
<rwgt>
<wgt id='1001'> 0.16394E+04 </wgt>
<wgt id='1002'> 0.16394E+04 </wgt>
<wgt id='1003'> 0.16394E+04 </wgt>
<wgt id='1004'> 0.18258E+04 </wgt>
<wgt id='1005'> 0.18258E+04 </wgt>
<wgt id='1006'> 0.18258E+04 </wgt>
<wgt id='1007'> 0.14364E+04 </wgt>
<wgt id='1008'> 0.14364E+04 </wgt>
<wgt id='1009'> 0.14364E+04 </wgt>
</rwgt>
</event>

Two failed events with weights +5.0756756e+06 and -6.1385739e+06:

<event>
 5 0 +5.0756756e+06 8.69829150e+01 7.54677110e-03 1.18849150e-01
       -3 -1 0 0 0 501 +0.0000000000e+00 +0.0000000000e+00 +3.3183067000e+03 3.3183067000e+03 0.0000000000e+00 0.0000e+00 9.0000e+00
        3 -1 0 0 501 0 +0.0000000000e+00 +0.0000000000e+00 -5.7002171000e-01 5.7002171000e-01 0.0000000000e+00 0.0000e+00 9.0000e+00
       23 2 1 2 0 0 +0.0000000000e+00 +0.0000000000e+00 +3.3177367000e+03 3.3188768000e+03 8.6982915000e+01 0.0000e+00 0.0000e+00
      -13 1 3 3 0 0 -6.1189663000e+00 +4.2667179000e+01 +1.4377739000e+03 1.4384199000e+03 0.0000000000e+00 0.0000e+00 9.0000e+00
       13 1 3 3 0 0 +6.1189663000e+00 -4.2667179000e+01 +1.8799628000e+03 1.8804569000e+03 0.0000000000e+00 0.0000e+00 9.0000e+00
#aMCatNLO 1 0 0 0 0 0.00000000E+00 0.00000000E+00 9 5 0 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00
<rwgt>
<wgt id='1001'> 0.16395E+04 </wgt>
<wgt id='1002'> 0.16395E+04 </wgt>
<wgt id='1003'> 0.16395E+04 </wgt>
<wgt id='1004'> 0.17570E+04 </wgt>
<wgt id='1005'> 0.17570E+04 </wgt>
<wgt id='1006'> 0.17570E+04 </wgt>
<wgt id='1007'> 0.14848E+04 </wgt>
<wgt id='1008'> 0.14848E+04 </wgt>
<wgt id='1009'> 0.14848E+04 </wgt>
</rwgt>
</event>
<event>
 5 0 -6.1385739e+06 9.53447640e+01 7.54677110e-03 1.17207690e-01
        2 -1 0 0 501 0 +0.0000000000e+00 +0.0000000000e+00 +3.2519728000e+03 3.2519728000e+03 0.0000000000e+00 0.0000e+00 9.0000e+00
       -2 -1 0 0 0 501 +0.0000000000e+00 +0.0000000000e+00 -6.9885455000e-01 6.9885455000e-01 0.0000000000e+00 0.0000e+00 9.0000e+00
       23 2 1 2 0 0 +0.0000000000e+00 +0.0000000000e+00 +3.2512740000e+03 3.2526717000e+03 9.5344764000e+01 0.0000e+00 0.0000e+00
      -13 1 3 3 0 0 +1.5191224000e+01 +1.3061496000e+01 +1.4988705000e+02 1.5122005000e+02 0.0000000000e+00 0.0000e+00 9.0000e+00
       13 1 3 3 0 0 -1.5191224000e+01 -1.3061496000e+01 +3.1013869000e+03 3.1014516000e+03 0.0000000000e+00 0.0000e+00 9.0000e+00
#aMCatNLO 1 0 0 0 0 0.00000000E+00 0.00000000E+00 9 5 0 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00
<rwgt>
<wgt id='1001'> 0.16395E+04 </wgt>
<wgt id='1002'> 0.16395E+04 </wgt>
<wgt id='1003'> 0.16395E+04 </wgt>
<wgt id='1004'> 0.17024E+04 </wgt>
<wgt id='1005'> 0.17024E+04 </wgt>
<wgt id='1006'> 0.17024E+04 </wgt>
<wgt id='1007'> 0.15562E+04 </wgt>
<wgt id='1008'> 0.15562E+04 </wgt>
<wgt id='1009'> 0.15562E+04 </wgt>
</rwgt>
</event>

Question information

Language:
English Edit question
Status:
Solved
For:
MadGraph5_aMC@NLO Edit question
Assignee:
Rikkert Frederix Edit question
Solved by:
Andrin Kessler
Solved:
Last query:
Last reply:
Revision history for this message
Rikkert Frederix (frederix) said :
#1

Dear Andrin,

I'm confused. What are you trying to do? As explained in arXiv:1412.8408, the ickkw=-1 is used for the jet-veto implementation according to scheme B of that paper. Why are you using it to generate events and then to reweight them with the virtual corrections?

Best regards,
Rikkert

Revision history for this message
Rikkert Frederix (frederix) said :
#2

Sorry, forget my previous message. I think I understood. I'll have a look at get back to you.

Revision history for this message
Rikkert Frederix (frederix) said :
#3

Dear Andrin,

Okay, I've reproduced your events and tried to reweight them. All goes smooth without problems. I tried both

p p > mu+ mu- [QCD]

and

p p > z > mu+ mu- [QCD].

(for the latter, I hacked the channelToShellName() function of the virt_reweighter.py so that it finds the correct folders/directories. Did you do the same?).

I always get something like this, i.e., where the 'max' is very close to the 'avg':

INFO:root:Virtual weight statistics,
  avg : 4.09E+03
  cum. var. : 1.82E+01
  max : 4.13E+03
INFO:root:Virtual cross-section obtained
  4.091777E+03 +/- 4.09E+01

Is there anything else for which you did not use the defaults?

Best,
Rikkert

Revision history for this message
Andrin Kessler (nirdna) said :
#4

Dear Rikkert,

Yes, I also altered the channelToShellName() function in the case of the intermediate particle.

Unfortunately, I will only have access to my machine tomorrow again but I will check the numbers (at least the average seems to be the same from memory) and report back when I have them.

The only other non-default options I could think of are the loop_sm-no_b_mass model I used to generate the events, and I also arbitrarily increased the bwcutoff parameter so that the eventfiles dont exclude any intermediate Z bosons regardless of kinematics since that introduced problems for some other code. I assumed this is unrelated.

I used the options

order = LO
fixed_order = OFF
shower = OFF
madspin = OFF
reweight = OFF

when generating the events.

Thanks for your efforts and best wishes,
Andrin

Revision history for this message
Andrin Kessler (nirdna) said :
#5

Dear Rikkert,

I ran the reweighter a few times and get these results:

INFO:root:Virtual weight statistics,
  avg : -7.72E+03
  cum. var. : 4.88E+06
  max : 2.34E+08
INFO:root:Virtual cross-section obtained
  -7.723259E+03 +/- 4.88E+04

So not at all the same. I do at least get the same result every time for the same events. This was done with the above linked 10'000 events run.

I have done a 100 event test run before and there i get :

INFO:root:Virtual weight statistics,
  avg : 1.47E+06
  cum. var. : 1.43E+07
  max : 1.45E+08
INFO:root:Virtual cross-section obtained
  1.466396E+06 +/- 1.45E+06

I note again that there are only a handful of events where this happens but the large results obviously dominate the average.

Best,
Andrin

Revision history for this message
Rikkert Frederix (frederix) said :
#6

Hi Andrin,

Humm.. this is strange. Could you try changing the reduction library that's used by MadLoop. I.e., changing

#MLReductionLib
!6|7|1

to

#MLReductionLib
6

(or any of the other 3 numbers there; note also that the exclamation mark needs to be removed) in the MadLoopParams.dat file? I think you probably need to change this file in both the ./Cards/ directory and in the ./SubProcesses/ directory. Do all reduction libraries give problems, or only specific ones?

Best,
Rikkert

Revision history for this message
Andrin Kessler (nirdna) said :
#7

Dear Rikkert,

Using 6 i get the same result as before.

7 gives the following error (I probably lack this library):

INFO:root: ... Done.
Traceback (most recent call last):
  File "./virt_reweighter.py", line 317, in <module>
    mu_r=event.scale)
  File "/home/andrin/Dev/pdfQFT/MG5_aMC_v2_6_6/bin/ppzmupmum2/../../madgraph/various/process_checks.py", line 1848, in get_me_value
    (ret_code, last_non_empty, error)
madgraph.MadGraph5Error: The MadLoop stability checker crashed with return code = 0, and last output:

stdout: mdl_WH = 6.3823389999999999E-003

stderr: STOP No available loop reduction lib is provided. Make sure MLReductionLib is correct.

Where the second to last line stdout: mdl_WH is either the value above, zero, or reads "internal Params".

Using 1 seems to work fine! I get:

INFO:root:Virtual weight statistics,
  avg : 4.09E+03
  cum. var. : 2.55E+00
  max : 4.32E+03
INFO:root:Virtual cross-section obtained
  4.086930E+03 +/- 4.09E+01

Best,
Andrin

Revision history for this message
Rikkert Frederix (frederix) said :
#8

Hi Andrin,

Okay, so it seems some issue with Ninja. Do you have a system-wide installation for ninja, or one local to MG5_aMC? If it's the former, I guess there are some incompatibilities with your version and the MG5_aMC approved version. If it's the latter, maybe try re-installing it, so that you get for sure the latest version. Should be possible with an 'install ninja' from the command prompt. Note that ninja uses internally OneLOop. Hence, it might also be an incompatibility with OneLOop instead of Ninja...

Anyway, good that it works with CutTools as reduction library. It's not the most efficient, but for this process that should be just fine.

best,
Rikkert

Revision history for this message
Andrin Kessler (nirdna) said :
#9

Dear Rikkert,

For my purposes this will do just fine but I will try a full reinstall soon anyway and keep an eye out for errors of this kind.

Thank you for your help!

Best wishes,
Andrin