diff --git a/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/analysis.f b/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/analysis.f
new file mode 100644
index 00000000..511311d6
--- /dev/null
+++ b/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/analysis.f
@@ -0,0 +1,83 @@
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_begin(nwgt,weights_info)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      integer nwgt
+      character*(*) weights_info(*)
+
+      call HwU_inithist(nwgt,weights_info)
+      call HwU_book(1,'mll', 6, 0d0, 6d0)
+      return
+      end
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_end(dummy)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      double precision dummy
+      call HwU_write_file
+      return
+      end
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      include 'nexternal.inc'
+      include 'cuts.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      integer ibody
+      double precision p(0:4,nexternal)
+      double precision wgts(*)
+      double precision ppl(0:3), pplb(0:3), ppv(0:3), xmll
+      double precision obs
+      integer i
+
+      double precision p_reco(0:4,nexternal)
+      integer iPDG_reco(nexternal),grid
+
+      call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
+     $                       p, iPDG, p_reco, iPDG_reco)
+
+      do i = nincoming+1, nexternal
+        if (iPDG_reco(i).eq.13) ppl(0:3)=p_reco(0:3,i)
+        if (iPDG_reco(i).eq.-13) pplb(0:3)=p_reco(0:3,i)
+      enddo
+      do i=0,3
+        ppv(i)=ppl(i)+pplb(i)
+      enddo
+
+      xmll=sqrt(ppv(0)**2-ppv(1)**2-ppv(2)**2-ppv(3)**2)
+      obs=-1d0
+
+      if (xmll.lt.1000d0) then
+        obs = -2d0
+      else if (xmll.lt.1500d0) then
+        obs=0.5d0
+      else if (xmll.lt.2000d0) then
+        obs=1.5d0
+      else if (xmll.lt.3000d0) then
+        obs=2.5d0
+      else if (xmll.lt.4000d0) then
+        obs=3.5d0
+      else if (xmll.lt.5000d0) then
+        obs=4.5d0
+      else if (xmll.lt.7000d0) then
+        obs=5.5d0
+      else
+        obs=-3d0
+      endif
+
+      if (obs.gt.0d0) then
+        call HwU_fill(1,obs,wgts)
+      else
+        write (*,*) "error: event outside bins", obs
+        stop 1
+      endif
+
+ 999  return
+      end
diff --git a/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/change_scale_to_mll.patch b/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/change_scale_to_mll.patch
new file mode 100644
index 00000000..55b0ba10
--- /dev/null
+++ b/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/change_scale_to_mll.patch
@@ -0,0 +1,67 @@
+--- CMSDY2D11/SubProcesses/setscales.f  2020-05-21 17:23:55.126143088 +0200
++++ CMSDY2D11/SubProcesses/setscales.f.new  2020-05-21 17:27:26.262700419 +0200
+@@ -527,6 +527,17 @@
+       integer i,j
+       character*80 temp_scale_id
+       common/ctemp_scale_id/temp_scale_id
++      integer iPDG_reco(nexternal)
++      double precision ppl(0:3), pplb(0:3), ppv(0:3), xmll
++      double precision p_reco(0:4,nexternal), p_in(0:4,nexternal)
++c     les houches accord stuff to identify particles
++c
++      integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
++     &        icolup(2,nexternal,maxflow),niprocs
++      common /c_leshouche_inc/idup,mothup,icolup,niprocs
++c Masses of external particles
++      double precision pmass(nexternal)
++      common/to_mass/pmass
+ c
+       tmp=0
+       if(ickkw.eq.-1)then
+@@ -568,10 +579,42 @@
+ cc                 dynamical_scale_choice = 10                                   cc
+ cc      in the run_card (run_card.dat)                                           cc
+ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+-         write(*,*) "User-defined scale not set"
+-         stop 1
+-         temp_scale_id='User-defined dynamical scale' ! use a meaningful string
+-         tmp = 0d0
++         temp_scale_id='Mll' ! use a meaningful string
++         tmp = -1d0
++         do i=1,nexternal
++           p_in(0:3,i) = pp(0:3,i)
++           p_in(4,i) = pmass(i)
++         enddo
++         call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
++     $                          p_in, idup(1,1), p_reco, iPDG_reco)
++
++         do j = nincoming+1, nexternal
++           if (iPDG_reco(j).eq.13) ppl(0:3)=p_reco(0:3,j)
++           if (iPDG_reco(j).eq.-13) pplb(0:3)=p_reco(0:3,j)
++         enddo
++         do i=0,3
++           ppv(i)=ppl(i)+pplb(i)
++         enddo
++
++         xmll=sqrt(ppv(0)**2-ppv(1)**2-ppv(2)**2-ppv(3)**2)
++
++         if (xmll.lt.1000d0) then
++           write (*,*) "error: event outside bins", xmll
++         else if (xmll.lt.1500d0) then
++           tmp=0.5*(1000d0+1500d0)
++         else if (xmll.lt.2000d0) then
++           tmp=0.5*(1500d0+2000d0)
++         else if (xmll.lt.3000d0) then
++           tmp=0.5*(2000d0+3000d0)
++         else if (xmll.lt.4000d0) then
++           tmp=0.5*(3000d0+4000d0)
++         else if (xmll.lt.5000d0) then
++           tmp=0.5*(4000d0+5000d0)
++         else if (xmll.lt.7000d0) then
++           tmp=0.5*(5000d0+7000d0)
++         else
++           write (*,*) "error: event outside bins", xmll
++         endif
+ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+ cc      USER-DEFINED SCALE: END OF USER CODE                                     cc
+ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
diff --git a/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/launch.txt b/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/launch.txt
new file mode 100644
index 00000000..a1f0ce6a
--- /dev/null
+++ b/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/launch.txt
@@ -0,0 +1,28 @@
+launch @OUTPUT@
+fixed_order = ON
+set maxjetflavor 5
+set gf @GF@
+set mh @MH@
+set mt @MT@
+set mw @MW@
+set mz @MZ@
+set wh @WH@
+set wt @WT@
+set ww @WW@
+set wz @WZ@
+set ebeam1 6500
+set ebeam2 6500
+set pdlabel lhapdf
+set lhaid @LHAPDF_ID@
+set dynamical_scale_choice 10
+set reweight_scale True
+set ptl = 10.0
+#user_defined_cut set ptl1min = 22.0
+set etal = 2.4
+set mll_sf = 1000.0
+#user_defined_cut set mmllmax = 7000.0
+set rphreco = 0.1
+set req_acc_FO 0.0001
+set pineappl True
+done
+quit
diff --git a/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/metadata.txt b/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/metadata.txt
new file mode 100644
index 00000000..5e3a34c4
--- /dev/null
+++ b/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/metadata.txt
@@ -0,0 +1,6 @@
+x1_label=Mll
+x1_label_tex=$M_{\ell\bar{\ell}}$
+x1_unit=GeV
+y_label=dsig/dMll
+y_label_tex=$\frac{\mathrm{d}\sigma}{\mathrm{d}M_{\ell\bar{\ell}}}$
+y_unit=pb/GeV
diff --git a/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/output.txt b/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/output.txt
new file mode 100644
index 00000000..2e44624d
--- /dev/null
+++ b/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/output.txt
@@ -0,0 +1,7 @@
+set complex_mass_scheme True
+import model loop_qcd_qed_sm_Gmu
+define p = p b b~
+define j = p
+generate p p > mu+ mu- [QCD QED]
+output @OUTPUT@
+quit
diff --git a/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/postrun.sh b/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/postrun.sh
new file mode 100755
index 00000000..5db22126
--- /dev/null
+++ b/nnpdf31_proc/PDF4LHC_DY_13_TEV_21_PHENO/postrun.sh
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+pineappl remap "${GRID}" "${GRID}".tmp '1000,1500,2000,3000,4000,5000,7000'
+mv "${GRID}".tmp "${GRID}"
diff --git a/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/analysis.f b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/analysis.f
new file mode 100644
index 00000000..f9988def
--- /dev/null
+++ b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/analysis.f
@@ -0,0 +1,83 @@
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_begin(nwgt,weights_info)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      integer nwgt
+      character*(*) weights_info(*)
+
+      call HwU_inithist(nwgt,weights_info)
+      call HwU_book(1,'mll', 7, 0d0, 7d0)
+      return
+      end
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_end(dummy)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      double precision dummy
+      call HwU_write_file
+      return
+      end
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+      subroutine analysis_fill(p,istatus,ipdg,wgts,ibody)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      implicit none
+      include 'nexternal.inc'
+      include 'cuts.inc'
+      integer istatus(nexternal)
+      integer iPDG(nexternal)
+      integer ibody
+      double precision p(0:4,nexternal)
+      double precision wgts(*)
+      double precision ppl(0:3), pplb(0:3), ppv(0:3), xmll
+      double precision obs
+      integer i
+
+      double precision p_reco(0:4,nexternal)
+      integer iPDG_reco(nexternal),grid
+
+      call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
+     $                       p, iPDG, p_reco, iPDG_reco)
+
+      do i = nincoming+1, nexternal
+        if (iPDG_reco(i).eq.14) ppl(0:3)=p_reco(0:3,i)
+        if (iPDG_reco(i).eq.-13) pplb(0:3)=p_reco(0:3,i)
+      enddo
+      do i=0,3
+        ppv(i)=ppl(i)+pplb(i)
+      enddo
+
+      xmll=sqrt(ppv(0)**2-ppv(1)**2-ppv(2)**2-ppv(3)**2)
+      obs=-1d0
+
+      if (xmll.lt.1000d0) then
+        obs=0.5d0
+      else if (xmll.lt.1500d0) then
+        obs=1.5d0
+      else if (xmll.lt.2000d0) then
+        obs=2.5d0
+      else if (xmll.lt.3000d0) then
+        obs=3.5d0
+      else if (xmll.lt.4000d0) then
+        obs=4.5d0
+      else if (xmll.lt.5000d0) then
+        obs=5.5d0
+      else if (xmll.lt.7000d0) then
+        obs=6.5d0
+      else
+        obs=-3d0
+      endif
+
+      if (obs.gt.0d0) then
+        call HwU_fill(1,obs,wgts)
+      else
+        !write (*,*) "error: event outside bins", obs
+        !stop 1
+      endif
+
+ 999  return
+      end
diff --git a/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/bias_large_invariant_mass.patch b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/bias_large_invariant_mass.patch
new file mode 100644
index 00000000..91543da4
--- /dev/null
+++ b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/bias_large_invariant_mass.patch
@@ -0,0 +1,49 @@
+--- PDF4LHC_WP_13_TEV_21_PHENO/SubProcesses/cuts.f.old	2022-01-18 14:44:46.095088533 +0100
++++ PDF4LHC_WP_13_TEV_21_PHENO/SubProcesses/cuts.f	2022-01-18 14:49:52.421021878 +0100
+@@ -975,25 +975,33 @@
+ c      'bias' = event_norm      
+ c
+       implicit none
++      include 'cuts.inc'
+       include 'nexternal.inc'
+       double precision bias_wgt,p(0:3,nexternal),H_T
+       integer ipdg(nexternal),i
++      double precision ppl(0:3), pplb(0:3), ppv(0:3)
++      double precision p_reco(0:4,nexternal)
++      integer iPDG_reco(nexternal),grid,xmlnu
++
+ 
+       bias_wgt=1d0
+ 
+-c How to enhance the tails is very process dependent. For example for
+-c top quark production one could use:
+-c      do i=1,nexternal
+-c         if (ipdg(i).eq.6) then
+-c            bias_wgt=sqrt(p(1,i)**2+p(2,i)**2)**3
+-c         endif
+-c      enddo
+-c Or to use H_T^2 one does     
+-c      H_T=0d0
+-c      do i=3,nexternal
+-c         H_T=H_T+sqrt(max(0d0,(p(0,i)+p(3,i))*(p(0,i)-p(3,i))))
+-c      enddo
+-c      bias_wgt=H_T**2
++      call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
++     $                       p, iPDG, p_reco, iPDG_reco)
++
++      do i = nincoming+1, nexternal
++        if (iPDG_reco(i).eq.14) ppl(0:3)=p_reco(0:3,i)
++        if (iPDG_reco(i).eq.-13) pplb(0:3)=p_reco(0:3,i)
++      enddo
++      do i=0,3
++        ppv(i)=ppl(i)+pplb(i)
++      enddo
++
++      xmlnu=sqrt(ppv(0)**2-ppv(1)**2-ppv(2)**2-ppv(3)**2)
++      if (xmlnu < 1000.0d0) then
++        bias_wgt=1d-10
++      endif
++      bias_wgt=xmlnu**3
+       return
+       end
+ 
diff --git a/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/change_scale_to_mll.patch b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/change_scale_to_mll.patch
new file mode 100644
index 00000000..e0a9ca84
--- /dev/null
+++ b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/change_scale_to_mll.patch
@@ -0,0 +1,67 @@
+--- CMSDY2D11/SubProcesses/setscales.f  2020-05-21 17:23:55.126143088 +0200
++++ CMSDY2D11/SubProcesses/setscales.f.new  2020-05-21 17:27:26.262700419 +0200
+@@ -527,6 +527,17 @@
+       integer i,j
+       character*80 temp_scale_id
+       common/ctemp_scale_id/temp_scale_id
++      integer iPDG_reco(nexternal)
++      double precision ppl(0:3), pplb(0:3), ppv(0:3), xmll
++      double precision p_reco(0:4,nexternal), p_in(0:4,nexternal)
++c     les houches accord stuff to identify particles
++c
++      integer idup(nexternal,maxproc),mothup(2,nexternal,maxproc),
++     &        icolup(2,nexternal,maxflow),niprocs
++      common /c_leshouche_inc/idup,mothup,icolup,niprocs
++c Masses of external particles
++      double precision pmass(nexternal)
++      common/to_mass/pmass
+ c
+       tmp=0
+       if(ickkw.eq.-1)then
+@@ -568,10 +579,42 @@
+ cc                 dynamical_scale_choice = 10                                   cc
+ cc      in the run_card (run_card.dat)                                           cc
+ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+-         write(*,*) "User-defined scale not set"
+-         stop 1
+-         temp_scale_id='User-defined dynamical scale' ! use a meaningful string
+-         tmp = 0d0
++         temp_scale_id='Mll' ! use a meaningful string
++         tmp = -1d0
++         do i=1,nexternal
++           p_in(0:3,i) = pp(0:3,i)
++           p_in(4,i) = pmass(i)
++         enddo
++         call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
++     $                          p_in, idup(1,1), p_reco, iPDG_reco)
++
++         do j = nincoming+1, nexternal
++           if (iPDG_reco(j).eq.14) ppl(0:3)=p_reco(0:3,j)
++           if (iPDG_reco(j).eq.-13) pplb(0:3)=p_reco(0:3,j)
++         enddo
++         do i=0,3
++           ppv(i)=ppl(i)+pplb(i)
++         enddo
++
++         xmll=sqrt(ppv(0)**2-ppv(1)**2-ppv(2)**2-ppv(3)**2)
++
++         if (xmll.lt.0000d0) then
++           tmp=0.5*(0d0+1000d0)
++         else if (xmll.lt.1500d0) then
++           tmp=0.5*(1000d0+1500d0)
++         else if (xmll.lt.2000d0) then
++           tmp=0.5*(1500d0+2000d0)
++         else if (xmll.lt.3000d0) then
++           tmp=0.5*(2000d0+3000d0)
++         else if (xmll.lt.4000d0) then
++           tmp=0.5*(3000d0+4000d0)
++         else if (xmll.lt.5000d0) then
++           tmp=0.5*(4000d0+5000d0)
++         else if (xmll.lt.7000d0) then
++           tmp=0.5*(5000d0+7000d0)
++         else
++           !write (*,*) "error: event outside bins", xmll
++         endif
+ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+ cc      USER-DEFINED SCALE: END OF USER CODE                                     cc
+ ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
diff --git a/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/launch.txt b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/launch.txt
new file mode 100644
index 00000000..3b1b520e
--- /dev/null
+++ b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/launch.txt
@@ -0,0 +1,35 @@
+launch @OUTPUT@
+fixed_order = ON
+set maxjetflavor 5
+set gf @GF@
+set mh @MH@
+set mt @MT@
+set mw @MW@
+set mz @MZ@
+set wh @WH@
+set wt @WT@
+set ww @WW@
+set wz @WZ@
+set ebeam1 6500
+set ebeam2 6500
+set pdlabel lhapdf
+set lhaid @LHAPDF_ID@
+set dynamical_scale_choice 10
+set reweight_scale True
+set event_norm = bias
+set ptl = 22.0
+set etal = 2.4
+#user_defined_cut set mmlnumax = 7000.0
+set rphreco = 0.1
+# with a bias function we degrade the uncertainty of the integrated cross
+# section, so don't optimize this quantity; instead, set the integration
+# parameters explicitly
+set req_acc_FO -1
+set npoints_fo_grid 500000 # step 0 points
+set niters_fo_grid 100 # step 0 iterations
+set npoints_fo 100000000 # step 1 points
+set niters_fo 50 # step 1 iterations
+#
+set pineappl True
+done
+quit
diff --git a/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/metadata.txt b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/metadata.txt
new file mode 100644
index 00000000..d03cc13d
--- /dev/null
+++ b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/metadata.txt
@@ -0,0 +1,6 @@
+x1_label=Mlnu
+x1_label_tex=$M_{\bar{\ell}\nu_\ell}$
+x1_unit=GeV
+y_label=dsig/dMlnu
+y_label_tex=$\frac{\mathrm{d}\sigma}{\mathrm{d}M_{\bar{\ell}\nu_\ell}}$
+y_unit=pb/GeV
diff --git a/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/output.txt b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/output.txt
new file mode 100644
index 00000000..e33260d5
--- /dev/null
+++ b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/output.txt
@@ -0,0 +1,7 @@
+set complex_mass_scheme True
+import model loop_qcd_qed_sm_Gmu
+define p = p b b~
+define j = p
+generate p p > mu+ vm [QCD QED]
+output @OUTPUT@
+quit
diff --git a/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/postrun.sh b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/postrun.sh
new file mode 100755
index 00000000..6d1fec43
--- /dev/null
+++ b/nnpdf31_proc/PDF4LHC_WP_13_TEV_21_PHENO/postrun.sh
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+pineappl remap "${GRID}" "${GRID}".tmp '0,1000,1500,2000,3000,4000,5000,7000'
+mv "${GRID}".tmp "${GRID}"
diff --git a/runcardsrunner/external/mg5/cuts_code/mmlnumax.f b/runcardsrunner/external/mg5/cuts_code/mmlnumax.f
new file mode 100644
index 00000000..ca5dd82a
--- /dev/null
+++ b/runcardsrunner/external/mg5/cuts_code/mmlnumax.f
@@ -0,0 +1,13 @@
+c     cut for mmlnumax (SFOS lepton pairs)
+      do i=1,nexternal-1
+        if (is_a_lm(i) .or. is_a_lp(i)) then
+          do j=1,nexternal
+          if (abs(ipdg_reco(i) + ipdg_reco(j)) .eq. 1) then
+              if (invm2_04(p_reco(0,i),p_reco(0,j),1d0) .gt. {}**2) then
+                passcuts_user=.false.
+                return
+              endif
+            endif
+          enddo
+        endif
+      enddo
diff --git a/runcardsrunner/external/mg5/cuts_code/mmlnumin.f b/runcardsrunner/external/mg5/cuts_code/mmlnumin.f
new file mode 100644
index 00000000..9dd7c380
--- /dev/null
+++ b/runcardsrunner/external/mg5/cuts_code/mmlnumin.f
@@ -0,0 +1,13 @@
+c     cut for mmlnumax (SFOS lepton pairs)
+      do i=1,nexternal-1
+        if (is_a_lm(i) .or. is_a_lp(i)) then
+          do j=1,nexternal
+          if (abs(ipdg_reco(i) + ipdg_reco(j)) .eq. 1) then
+              if (invm2_04(p_reco(0,i),p_reco(0,j),1d0) .lt. {}**2) then
+                passcuts_user=.false.
+                return
+              endif
+            endif
+          enddo
+        endif
+      enddo