Fewer events generated with custom cuts

I am attempting to generate events for the process p p > z w+ j j, z > l+ l-, w+ > l+ vl applying a custom cut on the pT of the Z boson. However, after implementing this cut in the dummy_fct.f file, I found that for moderately high values of the cut, MadGraph struggles to produce a sufficient number of events.

To verify that my custom cut is correctly implemented, I generated events for the simpler process p p > z j, z > l+ l- and applied the same cut. I compared these to events generated using the ptllmin and ptllmax cuts in the run_card (which should be equivalent). The results are consistent: with ptllmin = 3000. and ptllmax = 3100., I obtain a cross-section of 3.757e-08 ± 1.9e-10 , while my custom ptzmin and ptzmax cuts yield 3.734e-08 ± 5.7e-09. Additionally, I confirmed the proper implementation of my custom cuts using a short Python script to analyze the LHE files.

However, when I repeat this check for the p p > z j j, z > l+ l-, I encounter a discrepancy: with the ptll cut, events are generated quickly, and MadGraph has no trouble producing 10000 events. In contrast, using my custom cut significantly slows down the event generation, and MadGraph only manages to generate around 1000 events.

Given that the custom pT cut on Z should theoretically be equivalent to the ptll cut in the p p -> Z( Z -> l+ l-) j j process, why does MadGraph struggle to generate sufficient events with the custom cut?

Here I attach what my dummy function looks like:

logical FUNCTION dummy_cuts(P)

      IMPLICIT NONE
c
c Constants
c
      include 'genps.inc'
      include 'nexternal.inc'
      include 'cuts.inc'
      include 'run.inc'
C
C ARGUMENTS
C
      REAL*8 P(0:3,nexternal)
C
C PARAMETERS
C
      real*8 PI
      parameter( PI = 3.14159265358979323846d0 )
c
c particle identification
      LOGICAL IS_A_J(NEXTERNAL),IS_A_L(NEXTERNAL)
      LOGICAL IS_A_B(NEXTERNAL),IS_A_A(NEXTERNAL),IS_A_ONIUM(NEXTERNAL)
      LOGICAL IS_A_NU(NEXTERNAL),IS_HEAVY(NEXTERNAL)
      logical do_cuts(nexternal)
      COMMON /TO_SPECISA/IS_A_J,IS_A_A,IS_A_L,IS_A_B,IS_A_NU,IS_HEAVY,IS_A_ONIUM, do_cuts

c identifying the PDG id of the particles
c first, idup is created and then is filled with les houches info
      include 'maxamps.inc'
      integer idup(nexternal,maxproc,maxsproc)
      integer mothup(2,nexternal)
      integer icolup(2,nexternal,maxflow,maxsproc)
      include 'leshouche.inc'

c kinematic functions
      REAL*8 SumDot,pt,PtDot
c
      real*8 z_mass, best_mass_diff, tmp_mass, best_pt
      integer best_i, best_j, i, j

      dummy_cuts=.true.
c
      z_mass = 91.188d0
      best_mass_diff = 1.0d10**2
      best_i = -1
      best_j = -1
c
c Loop through all pairs of leptons of external particles
      do i = 3, nexternal-1
            if (.not. IS_A_L(i)) cycle
            do j = i+1, nexternal
            if (.not. IS_A_L(j)) cycle
c
c Check if the leptons have the same flavor
            if (idup(i,1,1) .ne. (-idup(j,1,1))) cycle
c
c Calculate the invariant mass of the lepton pair
            tmp_mass = dsqrt(SumDot(P(:,i),P(:,j),1.d0))
c
c Check if this pair is the best candidate for the Z boson
            if (abs(tmp_mass - z_mass) .LT. best_mass_diff) then
                  best_mass_diff = abs(tmp_mass - z_mass)
                  best_i = i
                  best_j = j
            endif
            enddo
      enddo
c
c Check if a valid pair was found
      if (best_i == -1 .or. best_j == -1) then
            return
      endif

c
c Calculate the transverse momentum (pT) of the sum of the momenta of the best pair
      best_pt = PtDot(P(:,best_i), P(:,best_j))

      if (ptzmax .le. ptzmin) then
            write(*,*) "Error: Inconsistent ptzmin and ptzmax values. They must be positive and ptzmax must be larger than ptzmin."
            return
      endif

     if ((ptzmin .ge. 0.0d0 .and. best_pt .lt. ptzmin*dabs(ptzmin)) .or.
     & (ptzmax .gt. 0.0d0 .and. best_pt .gt. ptzmax*dabs(ptzmax))) then
            dummy_cuts = .false.
            return
      endif

      return
      end

And here is the the full tag file of the Zjj process with my custom cuts:

<LesHouchesEvents version="3.0">
<header>
<!--
#*********************************************************************
# *
# MadGraph5_aMC@NLO *
# *
# Going Beyond *
# *
# http://madgraph.hep.uiuc.edu *
# http://madgraph.phys.ucl.ac.be *
# http://amcatnlo.cern.ch *
# *
# The MadGraph5_aMC@NLO team *
# *
#....................................................................*
# *
# This file contains all the information necessary to reproduce *
# the events generated: *
# *
# 1. software version *
# 2. proc_card : code generation info including model *
# 3. param_card : model primary parameters in the LH format *
# 4. run_card : running parameters (collider and cuts) *
# 5. pythia_card : present only if pythia has been run *
# 6. pgs_card : present only if pgs has been run *
# 7. delphes_cards : present only if delphes has been run *
# *
# *
#*********************************************************************
-->
<MGVersion>
3.5.3
</MGVersion>
<MG5ProcCard>
<![CDATA[
#************************************************************
#* MadGraph5_aMC@NLO *
#* *
#* * * *
#* * * * * *
#* * * * * 5 * * * * *
#* * * * * *
#* * * *
#* *
#* *
#* VERSION 3.5.3 2023-12-23 *
#* *
#* 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 group_subprocesses Auto
set ignore_six_quark_processes False
set low_mem_multicore_nlo_generation False
set complex_mass_scheme False
set include_lepton_initiated_processes False
set gauge unitary
set loop_optimized_output True
set loop_color_flows False
set max_npoint_for_channel 0
set default_unset_couplings 99
set max_t_for_channel 99
set zerowidth_tchannel True
set nlo_mixed_expansion True
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~
generate p p > z j j, z > l+ l-
output PP_Zjj
]]>
</MG5ProcCard>
<MGProcCard>
#*********************************************************************
# MadGraph/MadEvent *
# http://madgraph.hep.uiuc.edu *
# *
# proc_card.dat *
#*********************************************************************
# *
# This Files is generated by MADGRAPH 5 *
# *
# WARNING: This Files is generated for MADEVENT (compatibility issue)*
# This files is NOT a valid MG4 proc_card.dat *
# Running this in MG4 will NEVER reproduce the result of MG5*
# *
#*********************************************************************
#*********************************************************************
# Process(es) requested : mg2 input *
#*********************************************************************
# Begin PROCESS # This is TAG. Do not modify this line
p p > z j j , z > l+ l- #Process
# Be carefull the coupling are here in MG5 convention

end_coup # End the couplings input

done # this tells MG there are no more procs
# End PROCESS # This is TAG. Do not modify this line
#*********************************************************************
# Model information *
#*********************************************************************
# Begin MODEL # This is TAG. Do not modify this line
sm
# End MODEL # This is TAG. Do not modify this line
#*********************************************************************
# Start multiparticle definitions *
#*********************************************************************
# Begin MULTIPARTICLES # This is TAG. Do not modify this line

# End MULTIPARTICLES # This is TAG. Do not modify this line
</MGProcCard>
<MGRunCard>
<![CDATA[
#*********************************************************************
# MadGraph5_aMC@NLO *
# *
# run_card.dat MadEvent *
# *
# This file is used to set the parameters of the run. *
# *
# Some notation/conventions: *
# *
# Lines starting with a '# ' are info or comments *
# *
# mind the format: value = variable ! comment *
# *
# To display more options, you can type the command: *
# update to_full *
#*********************************************************************
#
#*********************************************************************
# Tag name for the run (one word) *
#*********************************************************************
  tag_1 = run_tag ! name of the run
#*********************************************************************
# Number of events and rnd seed *
# Warning: Do not generate more than 1M events in a single run *
#*********************************************************************
  10000 = nevents ! Number of unweighted events requested
 129 = iseed ! rnd seed (0=assigned automatically=default))
#*********************************************************************
# Collider type and energy *
# lpp: 0=No PDF, 1=proton, -1=antiproton, *
# 2=elastic photon of proton/ion beam *
# +/-3=PDF of electron/positron beam *
# +/-4=PDF of muon/antimuon beam *
#*********************************************************************
     1 = lpp1 ! beam 1 type
     1 = lpp2 ! beam 2 type
     7000.0 = ebeam1 ! beam 1 total energy in GeV
     7000.0 = ebeam2 ! beam 2 total energy in GeV
# To see polarised beam options: type "update beam_pol"

#*********************************************************************
# PDF CHOICE: this automatically fixes alpha_s and its evol. *
# pdlabel: lhapdf=LHAPDF (installation needed) [1412.7420] *
# iww=Improved Weizsaecker-Williams Approx.[hep-ph/9310350] *
# eva=Effective W/Z/A Approx. [2111.02442] *
# edff=EDFF in gamma-UPC [eq.(11) in 2207.03012] *
# chff=ChFF in gamma-UPC [eq.(13) in 2207.03012] *
# none=No PDF, same as lhapdf with lppx=0 *
#*********************************************************************
     lhapdf = pdlabel ! PDF set
     315000 = lhaid ! if pdlabel=lhapdf, this is the lhapdf number
# To see heavy ion options: type "update ion_pdf"
#*********************************************************************
# Renormalization and factorization scales *
#*********************************************************************
 False = fixed_ren_scale ! if .true. use fixed ren scale
     False = fixed_fac_scale ! if .true. use fixed fac scale
 91.188 = scale ! fixed ren scale
 91.188 = dsqrt_q2fact1 ! fixed fact scale for pdf1
 91.188 = dsqrt_q2fact2 ! fixed fact scale for pdf2
 2 = dynamical_scale_choice ! Choose one of the preselected dynamical choices
 1.0 = scalefact ! scale factor for event-by-event scales

#*********************************************************************
# Type and output format
#*********************************************************************
  False = gridpack !True = setting up the grid pack
  -1.0 = time_of_flight ! threshold (in mm) below which the invariant livetime is not written (-1 means not written)
  average = event_norm ! average/sum. Normalization of the weight in the LHEF
# To see MLM/CKKW merging options: type "update MLM" or "update CKKW"

#*********************************************************************
#
#*********************************************************************
# Phase-Space Optimization strategy (basic options)
#*********************************************************************
   0 = nhel ! using helicities importance sampling or not.
                             ! 0: sum over helicity, 1: importance sampling
   2 = sde_strategy ! default integration strategy (hep-ph/2021.00773)
                             ! 1 is old strategy (using amp square)
        ! 2 is new strategy (using only the denominator)
# To see advanced option for Phase-Space optimization: type "update psoptim"
#*********************************************************************
# Customization (custom cuts/scale/bias/...) *
# list of files containing fortran function that overwrite default *
#*********************************************************************
  = custom_fcts ! List of files containing user hook function
#*******************************
# Parton level cuts definition *
#*******************************
  0.0 = dsqrt_shat ! minimal shat for full process
#
#
#*********************************************************************
# BW cutoff (M+/-bwcutoff*Gamma) ! Define on/off-shell for "$" and decay
#*********************************************************************
  15.0 = bwcutoff ! (M+/-bwcutoff*Gamma)
 #*********************************************************************
 # Apply pt/E/eta/dr/mij/kt_durham cuts on decay products or not
 # (note that etmiss/ptll/ptheavy/ht/sorted cuts always apply)
 #*********************************************************************
   False = cut_decays ! Cut decay products
#*********************************************************************
# Standard Cuts *
#*********************************************************************
# Minimum and maximum pt's (for max, -1 means no cut) *
#*********************************************************************
 20.0 = ptj ! minimum pt for the jets
 10.0 = ptl ! minimum pt for the charged leptons
 -1.0 = ptjmax ! maximum pt for the jets
 -1.0 = ptlmax ! maximum pt for the charged leptons
 {} = pt_min_pdg ! pt cut for other particles (use pdg code). Applied on particle and anti-particle
 {} = pt_max_pdg ! pt cut for other particles (syntax e.g. {6: 100, 25: 50})
#
# For display option for energy cut in the partonic center of mass frame type 'update ecut'
#
#*********************************************************************
# Maximum and minimum absolute rapidity (for max, -1 means no cut) *
#*********************************************************************
 5.0 = etaj ! max rap for the jets
 2.5 = etal ! max rap for the charged leptons
 0.0 = etalmin ! main rap for the charged leptons
 {} = eta_min_pdg ! rap cut for other particles (use pdg code). Applied on particle and anti-particle
 {} = eta_max_pdg ! rap cut for other particles (syntax e.g. {6: 2.5, 23: 5})
#*********************************************************************
# Minimum and maximum DeltaR distance *
#*********************************************************************
 0.4 = drjj ! min distance between jets
 0.4 = drll ! min distance between leptons
 0.4 = drjl ! min distance between jet and lepton
 -1.0 = drjjmax ! max distance between jets
 -1.0 = drllmax ! max distance between leptons
 -1.0 = drjlmax ! max distance between jet and lepton
#*********************************************************************
# Minimum and maximum invariant mass for pairs *
#*********************************************************************
 0.0 = mmjj ! min invariant mass of a jet pair
 0.0 = mmll ! min invariant mass of l+l- (same flavour) lepton pair
 -1.0 = mmjjmax ! max invariant mass of a jet pair
 -1.0 = mmllmax ! max invariant mass of l+l- (same flavour) lepton pair
 {} = mxx_min_pdg ! min invariant mass of a pair of particles X/X~ (e.g. {6:250})
 {'default': False} = mxx_only_part_antipart ! if True the invariant mass is applied only
                       ! to pairs of particle/antiparticle and not to pairs of the same pdg codes.
 #*********************************************************************
 # Minimum and maximum invariant mass for all letpons *
 #*********************************************************************
 0.0 = mmnl ! min invariant mass for all letpons (l+- and vl)
 -1.0 = mmnlmax ! max invariant mass for all letpons (l+- and vl)
 #*********************************************************************
 # Minimum and maximum pt for 4-momenta sum of leptons / neutrino *
 # for pair of lepton includes only same flavor, opposite charge
 #*********************************************************************
 0.0 = ptllmin ! Minimum pt for 4-momenta sum of leptons(l and vl)
 -1.0 = ptllmax ! Maximum pt for 4-momenta sum of leptons(l and vl)
#*********************************************************************
# Inclusive cuts *
#*********************************************************************
 0.0 = xptj ! minimum pt for at least one jet
 0.0 = xptl ! minimum pt for at least one charged lepton
 #*********************************************************************
 # Control the pt's of the jets sorted by pt *
 #*********************************************************************
 0.0 = ptj1min ! minimum pt for the leading jet in pt
 0.0 = ptj2min ! minimum pt for the second jet in pt
 -1.0 = ptj1max ! maximum pt for the leading jet in pt
 -1.0 = ptj2max ! maximum pt for the second jet in pt
 0 = cutuse ! reject event if fails any (0) / all (1) jet pt cuts
 #*********************************************************************
 # Control the pt's of leptons sorted by pt *
 #*********************************************************************
 0.0 = ptl1min ! minimum pt for the leading lepton in pt
 0.0 = ptl2min ! minimum pt for the second lepton in pt
 -1.0 = ptl1max ! maximum pt for the leading lepton in pt
 -1.0 = ptl2max ! maximum pt for the second lepton in pt
 #*********************************************************************
 # Control the Ht(k)=Sum of k leading jets *
 #*********************************************************************
 0.0 = htjmin ! minimum jet HT=Sum(jet pt)
 -1.0 = htjmax ! maximum jet HT=Sum(jet pt)
 0.0 = ihtmin !inclusive Ht for all partons (including b)
 -1.0 = ihtmax !inclusive Ht for all partons (including b)
 #*********************************************************************
 # WBF cuts *
 #*********************************************************************
 0.0 = xetamin ! minimum rapidity for two jets in the WBF case
 0.0 = deltaeta ! minimum rapidity for two jets in the WBF case
#*********************************************************************
# maximal pdg code for quark to be considered as a light jet *
# (otherwise b cuts are applied) *
#*********************************************************************
 4 = maxjetflavor ! Maximum jet pdg code
#*********************************************************************
#
#*********************************************************************
# Store info for systematics studies *
# WARNING: Do not use for interference type of computation *
#*********************************************************************
   True = use_syst ! Enable systematics studies
#
systematics = systematics_program ! none, systematics [python], SysCalc [depreceted, C++]
['--mur=0.5,1,2', '--muf=0.5,1,2', '--pdf=errorset'] = systematics_arguments ! see: https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Systematics#Systematicspythonmodule

1 = hard_survey
3000 = ptzmin
3100 = ptzmax

]]>
</MGRunCard>
<slha>
######################################################################
## PARAM_CARD AUTOMATICALY GENERATED BY MG5 ####
######################################################################
###################################
## INFORMATION FOR MASS
###################################
BLOCK MASS #
      5 4.700000e+00 # mb
      6 1.730000e+02 # mt
      15 1.777000e+00 # mta
      23 9.118800e+01 # mz
      25 1.250000e+02 # mh
      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
      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.180021e-01 # as (note that parameter not used if you use a pdf set)
###################################
## INFORMATION FOR YUKAWA
###################################
BLOCK YUKAWA #
      5 4.700000e+00 # ymb
      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
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
</slha>
<MGGenerationInfo>
# Number of Events : 10000
# Integrated weight (pb) : 2.0902233753e-08
</MGGenerationInfo>
</header>
</LesHouchesEvents>

Any help would be very appreciated!

- Sergio R