Merge lp:~maddevelopers/mg5amcnlo/3.3.0 into lp:mg5amcnlo

[Olivier Mattelaer]

This includes:
 - Bug fixing from 2.9.6
 - EWA branch
 - New_clustering branch
 - tagged photon
 - lepton_from_proton branch
 - Move pythia8 support to pythya8.3

Known Bug:
 - jet veto cross-section (since 3.x.x)
 - issue with heft at NLO (since 3.x.x as well)

[lp:~maddevelopers/mg5amcnlo/3.3.0]

979. By olivier-mattelaer on 2021-11-08

merge with 2.7.0_gpu version 368->371: fixing randomness in ordering

980. By olivier-mattelaer on 2021-11-08

update IOTest (now should be reproducible version to version)

981. By olivier-mattelaer on 2021-11-10

force set to be ordered in debug mode (result more and more reproducible)

[Olivier Mattelaer (olivier-mattelaer)]

I have added a bunch of "python hack --applied only to debug mode--" to force set to be ordered and avoid the huge list of false positive in the IOTest.

If they are no objection to this, I will release it Thursday or Friday.

Olivier

[lp:~maddevelopers/mg5amcnlo/3.3.0]

982. By olivier-mattelaer on 2021-11-12

fixing some py2 ordering issue, fixing side effect of ordered set, fixing tests

[Olivier Mattelaer (olivier-mattelaer)]

let's release this.

review: Approve

=== modified file '.bzrignore'
--- .bzrignore	2013-03-17 12:01:13 +0000
+++ .bzrignore	2021-11-12 09:03:34 +0000
@@ -1,3 +1,4 @@
 models/RS/.DS_Store
 help
 ./tests/input_files/my_smallsm
+
 
=== modified file 'INSTALL'
--- INSTALL	2013-11-15 11:21:46 +0000
+++ INSTALL	2021-11-12 09:03:34 +0000
@@ -24,7 +24,7 @@
 
  * MadGraph5_aMC@NLO *
 
-- Python 2.6 (or higher, but not compatible with 3.X)
+- Python 3.7 (or higher, python2.7 is also working for the time being)
 
  * MadEvent *
 Package for the LO cross-section computation and generation of events 
 
=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2021-05-23 14:17:53 +0000
+++ MadSpin/decay.py	2021-11-12 09:03:34 +0000
@@ -178,7 +178,7 @@
             part=self.event2mg[item]
             if part>0:
                 particle_line=self.get_particle_line(self.particle[part])
-                if abs(self.particle[part]["istup"]) == 1:
+                if abs(self.particle[part]["istup"]) == 1 or abs(self.particle[part]["istup"]) == 2:
                     if "pt_scale" in self.particle[part]:
                         scales.append(self.particle[part]["pt_scale"])
                     else:
@@ -187,17 +187,16 @@
                 particle_line=self.get_particle_line(self.resonance[part])
             line+=particle_line        
         
+        if self.diese:
+            line += self.diese
+        if self.rwgt:
+            line += self.rwgt
         if any(scales):
             sqrts = self.particle[1]["pt_scale"]
-            line += "<scales %s></scales>\n" % ' '.join(['pt_clust_%i=\"%s\"' 
-                                                        %(i-1,s if s else sqrts)
+            line += "<scales %s></scales>\n" % ' '.join(['pt_clust_%i=\"%s\"'
+                                                        %(i+1,s if s else self.scale)
                                                        for i,s in enumerate(scales)
                                                        if i>1])
-        
-        if self.diese:
-            line += self.diese
-        if self.rwgt:
-            line += self.rwgt
         line+="</event> \n"
         return line
 
 
=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py	2021-05-23 07:47:20 +0000
+++ MadSpin/interface_madspin.py	2021-11-12 09:03:34 +0000
@@ -85,7 +85,7 @@
     ############################################################################
     ##  Special post-processing of the options                                ## 
     ############################################################################
-    def post_set_ms_dir(self, value, change_userdefine, raiseerror):
+    def post_set_ms_dir(self, value, change_userdefine, raiseerror, *opts):
         """ special handling for set ms_dir """
         
         self.__setitem__('curr_dir', value, change_userdefine=change_userdefine)
@@ -97,7 +97,7 @@
         random.seed(value)
 
     ############################################################################        
-    def post_set_run_card(self, value, change_userdefine, raiseerror):
+    def post_set_run_card(self, value, change_userdefine, raiseerror, *opts):
         """ special handling for set run_card """
         
         if value == 'default':
 
=== modified file 'MadSpin/madspin'
--- MadSpin/madspin	2021-06-29 19:12:21 +0000
+++ MadSpin/madspin	2021-11-12 09:03:34 +0000
@@ -1,4 +1,4 @@
-#! /usr/bin/env python
+#! /usr/bin/env python3
 
 ################################################################################
 #
@@ -20,8 +20,8 @@
 from __future__ import absolute_import
 from __future__ import print_function
 import sys
-if not sys.version_info[0] == 2 or sys.version_info[1] < 6:
-    sys.exit('MadSpin works only with python 2.6 or later (but not python 3.X).\n\
+if  sys.version_info[1] < 6:
+    sys.exit('MadSpin works only with python 2.7 or with python 3.6 (or later).\n\
                Please upgrate your version of python.')
 
 import os
 
=== modified file 'Template/Common/Source/PDF/pdf.inc'
--- Template/Common/Source/PDF/pdf.inc	2021-04-16 08:31:20 +0000
+++ Template/Common/Source/PDF/pdf.inc	2021-11-12 09:03:34 +0000
@@ -6,6 +6,7 @@
 c     pdlabel= string identifying the pdf    
 c***********************************************************************
       character*7 pdlabel,epa_label
+      character*7 pdsublabel(2)
       integer lhaid
-      common/to_pdf/lhaid,pdlabel,epa_label
+      common/to_pdf/lhaid,pdlabel,epa_label,pdsublabel
  	
 
=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2021-01-08 16:12:55 +0000
+++ Template/LO/Cards/run_card.dat	2021-11-12 09:03:34 +0000
@@ -27,8 +27,9 @@
   %(iseed)s   = iseed   ! rnd seed (0=assigned automatically=default))
 #*********************************************************************
 # Collider type and energy                                           *
-# lpp: 0=No PDF, 1=proton, -1=antiproton, 2=photon from proton,      *
-#                3=photon from electron, 4=photon from muon          *
+# lpp: 0=No PDF, 1=proton, -1=antiproton, 2=elastic photon of proton,*
+#             +/-3=PDF of electron/positron beam                     *
+#             +/-4=PDF of muon/antimuon beam                         *
 #*********************************************************************
      %(lpp1)s        = lpp1    ! beam 1 type 
      %(lpp2)s        = lpp2    ! beam 2 type
@@ -37,21 +38,26 @@
 $beam_pol
 $frame
 #*********************************************************************
-# PDF CHOICE: this automatically fixes also alpha_s and its evol.    *
+# 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.       [21yy.zzzzz]            *
+#          none=No PDF, same as lhapdf with lppx=0                   *
 #*********************************************************************
-     %(pdlabel)s    = pdlabel     ! PDF set                                     
+$pdlabel
      %(lhaid)s    = lhaid     ! if pdlabel=lhapdf, this is the lhapdf number
 $ion_pdf
 #*********************************************************************
 # Renormalization and factorization scales                           *
 #*********************************************************************
  %(fixed_ren_scale)s = fixed_ren_scale  ! if .true. use fixed ren scale
- %(fixed_fac_scale)s        = fixed_fac_scale  ! if .true. use fixed fac scale
+$fixed_fact_scale
  %(scale)s  = scale            ! fixed ren scale
  %(dsqrt_q2fact1)s  = dsqrt_q2fact1    ! fixed fact scale for pdf1
  %(dsqrt_q2fact2)s  = dsqrt_q2fact2    ! fixed fact scale for pdf2
  %(dynamical_scale_choice)s = dynamical_scale_choice ! Choose one of the preselected dynamical choices
  %(scalefact)s  = scalefact        ! scale factor for event-by-event scales
+ $eva_scale
 #*********************************************************************
 # Type and output format
 #*********************************************************************
@@ -80,7 +86,8 @@
 #
 #*******************************                                                 
 # Parton level cuts definition *
-#*******************************                                     
+#*******************************
+  %(dsqrt_shat)s  = dsqrt_shat ! minimal shat for full process
 #                                                                    
 #
 #*********************************************************************
 
=== added file 'Template/LO/Source/PDF/ElectroweakFlux.f'
--- Template/LO/Source/PDF/ElectroweakFlux.f	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/PDF/ElectroweakFlux.f	2021-11-12 09:03:34 +0000
@@ -0,0 +1,225 @@
+c     /* ********************************************************* *
+c                 Effective Vector Boson Approximation
+c     /* ********************************************************* *      
+c     File: ElectroweakFlux.f
+c     R. Ruiz (2021 February)
+c     For details, see companion paper by ..., et al [arXiv:]
+c     /* ********************************************************* *
+c     /* ********************************************************* *      
+c     function eva_fX_to_vV(gg2,gL2,gR2,fLpol,mv2,x,mu2,ievo):
+c     call electroweak PDF for for vector boson (hel=V=0,+,-) from fermion
+c     with fractional (0<pol<1) LH (RH) polarization of fLpol (1-fLpol)
+c     - fLpol = 0.5d0 = parent fermion is spin-averaged
+c     /* ********************************************************* *
+c     /* ********************************************************* *      
+c     function eva_fF_to_vV(gg2,gF2,mv2,x,mu2,ievo):
+c     electroweak PDF for vector boson (hel=V=0,+,-) from fermion (hel=F=L/R)
+c     - gg2 = (coupling)^2
+c     - gF2 = (L/R chiral coupling)^2
+c     - mv2 = (mass of boson)^2
+c     - x = momentum fraction = (E_V / E_f)
+c     - mu2 = (max of evolution scale)^2
+c     - ievo = for evolution by virtuality (or !=0 for pT)      
+c     /* ********************************************************* *
+c     /* ********************************************************* *
+c     /* ********************************************************* *            
+c     /* ********************************************************* *
+      double precision function eva_fX_to_vp(gg2,gL2,gR2,fLpol,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo
+      double precision gg2,gL2,gR2,fLpol,mv2,x,mu2
+      double precision eva_fL_to_vp,eva_fR_to_vp
+
+      eva_fX_to_vp =       fLpol*eva_fL_to_vp(gg2,gL2,mv2,x,mu2,ievo)
+     &             + (1d0-fLpol)*eva_fR_to_vp(gg2,gR2,mv2,x,mu2,ievo)
+      return
+      end
+c     /* ********************************************************* *
+      double precision function eva_fX_to_vm(gg2,gL2,gR2,fLpol,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo
+      double precision gg2,gL2,gR2,fLpol,mv2,x,mu2
+      double precision eva_fL_to_vm,eva_fR_to_vm
+
+      eva_fX_to_vm =       fLpol*eva_fL_to_vm(gg2,gL2,mv2,x,mu2,ievo)
+     &             + (1d0-fLpol)*eva_fR_to_vm(gg2,gR2,mv2,x,mu2,ievo)
+      return
+      end
+c     /* ********************************************************* *
+      double precision function eva_fX_to_v0(gg2,gL2,gR2,fLpol,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo
+      double precision gg2,gL2,gR2,fLpol,mv2,x,mu2
+      double precision eva_fL_to_v0,eva_fR_to_v0
+
+      eva_fX_to_v0 =       fLpol*eva_fL_to_v0(gg2,gL2,mv2,x,mu2,ievo)
+     &             + (1d0-fLpol)*eva_fR_to_v0(gg2,gR2,mv2,x,mu2,ievo)
+      return
+      end
+c     /* ********************************************************* *
+      double precision function eva_fX_to_fL(gg2,gL2,gR2,fLpol,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo
+      double precision gg2,gL2,gR2,fLpol,mv2,x,mu2
+      double precision eva_fL_to_fL,eva_fR_to_fL
+
+      eva_fX_to_fL =       fLpol*eva_fL_to_fL(gg2,gL2,mv2,x,mu2,ievo)
+     &             + (1d0-fLpol)*eva_fR_to_fL(gg2,gR2,mv2,x,mu2,ievo)
+      return
+      end
+c     /* ********************************************************* *
+      double precision function eva_fX_to_fR(gg2,gL2,gR2,fLpol,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo
+      double precision gg2,gL2,gR2,fLpol,mv2,x,mu2
+      double precision eva_fL_to_fR,eva_fR_to_fR
+
+      eva_fX_to_fR =       fLpol*eva_fL_to_fR(gg2,gL2,mv2,x,mu2,ievo)
+     &             + (1d0-fLpol)*eva_fR_to_fR(gg2,gR2,mv2,x,mu2,ievo)
+      return
+      end      
+c     /* ********************************************************* *      
+c     EVA (1/6) for f_L > v_+
+      double precision function eva_fL_to_vp(gg2,gL2,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo              ! evolution by q2 or pT2
+      double precision gg2,gL2,mv2,x,mu2
+      double precision coup2,split,xxlog,fourPiSq
+      data fourPiSq/39.47841760435743d0/ ! = 4pi**2
+
+c      print*,'gg2,gL2,mv2,x,mu2,ievo',gg2 !3,gL2,mv2,x,mu2,ievo
+      coup2 = gg2*gL2/fourPiSq
+      split = (1.d0-x)**2 / 2.d0 / x
+      if(ievo.eq.0) then
+         xxlog = dlog(mu2/mv2)
+      else
+         xxlog = dlog(mu2/mv2/(1.d0-x))
+      endif
+      
+      eva_fL_to_vp = coup2*split*xxlog
+      return
+      end
+c     /* ********************************************************* *
+c     EVA (2/6) for f_L > v_-
+      double precision function eva_fL_to_vm(gg2,gL2,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo              ! evolution by q2 or pT2
+      double precision gg2,gL2,mv2,x,mu2
+      double precision coup2,split,xxlog,fourPiSq
+      data fourPiSq/39.47841760435743d0/ ! = 4pi**2
+      
+      coup2 = gg2*gL2/fourPiSq
+      split = 1.d0 / 2.d0 / x
+      if(ievo.eq.0) then
+         xxlog = dlog(mu2/mv2)
+      else
+         xxlog = dlog(mu2/mv2/(1.d0-x))
+      endif
+
+      eva_fL_to_vm = coup2*split*xxlog
+      return
+      end
+c     /* ********************************************************* *
+c     EVA (3/6) for f_L > v_0
+      double precision function eva_fL_to_v0(gg2,gL2,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo              ! evolution by q2 or pT2
+      double precision gg2,gL2,mv2,x,mu2
+      double precision coup2,split,xxlog,fourPiSq
+      data fourPiSq/39.47841760435743d0/ ! = 4pi**2
+c
+      coup2 = gg2*gL2/fourPiSq
+      split = (1.d0-x) / x
+      xxlog = 1.d0
+      
+      eva_fL_to_v0 = coup2*split*xxlog
+      return
+      end
+c     /* ********************************************************* *
+c     EVA (4/6) for f_R > v_+
+      double precision function eva_fR_to_vp(gg2,gR2,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo              ! evolution by q2 or pT2
+      double precision gg2,gR2,mv2,x,mu2
+      double precision eva_fL_to_vm
+
+      eva_fR_to_vp = eva_fL_to_vm(gg2,gR2,mv2,x,mu2,ievo)
+      return
+      end
+c     /* ********************************************************* *
+c     EVA (5/6) for f_R > v_-
+      double precision function eva_fR_to_vm(gg2,gR2,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo              ! evolution by q2 or pT2
+      double precision gg2,gR2,mv2,x,mu2
+      double precision eva_fL_to_vp
+
+      eva_fR_to_vm = eva_fL_to_vp(gg2,gR2,mv2,x,mu2,ievo)
+      return
+      end      
+c     /* ********************************************************* *
+c     EVA (6/6) for f_R > v_0
+      double precision function eva_fR_to_v0(gg2,gR2,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo              ! evolution by q2 or pT2
+      double precision gg2,gR2,mv2,x,mu2
+      double precision eva_fL_to_v0
+
+      eva_fR_to_v0 = eva_fL_to_v0(gg2,gR2,mv2,x,mu2,ievo)
+      return
+      end
+c     /* ********************************************************* *  
+c     EVA () for f_L > f_L
+c     fL_to_fL(z) = fL_to_vp(1-z) + fL_to_vm(1-z) 
+      double precision function eva_fL_to_fL(gg2,gL2,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo              ! evolution by q2 or pT2
+      double precision gg2,gL2,mv2,x,mu2
+      double precision tmpVp,tmpVm,z
+      double precision eva_fL_to_vp,eva_fL_to_vm
+
+      z = 1.d0 - x
+      tmpVp = eva_fL_to_vp(gg2,gL2,mv2,z,mu2,ievo)
+      tmpVm = eva_fL_to_vm(gg2,gL2,mv2,z,mu2,ievo)
+
+      eva_fL_to_fL = tmpVp + tmpVm
+      return
+      end
+c     /* ********************************************************* *  
+c     EVA () for f_R > f_R
+c     fR_to_fR(z) = fR_to_vp(1-z) + fR_to_vm(1-z) 
+      double precision function eva_fR_to_fR(gg2,gR2,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo              ! evolution by q2 or pT2
+      double precision gg2,gR2,mv2,x,mu2
+      double precision tmpVp,tmpVm,z
+      double precision eva_fR_to_vp,eva_fR_to_vm
+
+      z = 1.d0 - x
+      tmpVp = eva_fR_to_vp(gg2,gR2,mv2,z,mu2,ievo)
+      tmpVm = eva_fR_to_vm(gg2,gR2,mv2,z,mu2,ievo)
+
+      eva_fR_to_fR = tmpVp + tmpVm
+      return
+      end      
+c     /* ********************************************************* *  
+c     EVA () for f_L > f_R
+      double precision function eva_fL_to_fR(gg2,gL2,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo              ! evolution by q2 or pT2
+      double precision gg2,gL2,mv2,x,mu2
+
+      eva_fL_to_fR = 0d0
+      return
+      end
+c     /* ********************************************************* *       
+c     EVA () for f_R > f_L
+      double precision function eva_fR_to_fL(gg2,gR2,mv2,x,mu2,ievo)
+      implicit none
+      integer ievo              ! evolution by q2 or pT2
+      double precision gg2,gR2,mv2,x,mu2
+
+      eva_fR_to_fL = 0d0
+      return
+      end
+c     /* ********************************************************* *       
 \ No newline at end of file
 
=== added file 'Template/LO/Source/PDF/ElectroweakFluxDriver.f'
--- Template/LO/Source/PDF/ElectroweakFluxDriver.f	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/PDF/ElectroweakFluxDriver.f	2021-11-12 09:03:34 +0000
@@ -0,0 +1,561 @@
+c     /* ********************************************************* *
+c                 Effective Vector Boson Approximation
+c     /* ********************************************************* *      
+c     File: ElectroweakFluxDriver.f
+c     R. Ruiz (2021 February)
+c     For details, see companion paper by Costantini, et al [arXiv:]
+c     /* ********************************************************* *
+c     function eva_get_pdf_by_PID:
+c     - wrapper for eva_get_pdf_by_PID_evo
+c     function eva_get_pdf_by_PID_evo
+c     - set eva PDF couplings by PIDs
+c     - call V_+,V_-,V_0 PDF by v polarization (vpol)
+c     - call PDF for f_L,f_R by fL polarization (fLpol; fLpol=0.5 = unpolarized)
+c     subroutine eva_get_mv2_by_PID
+c     - assign mass by vPID
+c     subroutine eva_get_mf2_by_PID      
+c     - assign mass by fPID
+c     subroutine eva_get_gg2_by_PID
+c     - assign universal coupling strength by vPID
+c     subroutine eva_get_gR2_by_PID
+c     - assign right couplings of fermion by vPID and fPID      
+c     subroutine eva_get_gL2_by_PID
+c     - assign left couplings of fermion by vPID and fPID
+c     /* ********************************************************* *
+      double precision function eva_get_pdf_by_PID(vPID,fPID,vpol,fLpol,x,mu2,ievo)
+      implicit none
+      integer ievo ! =0 for evolution by q^2 (!=0 for evolution by pT^2)
+      integer vPID,fPID,vpol
+      double precision fLpol,x,mu2
+      double precision eva_get_pdf_by_PID_evo
+      double precision eva_get_pdf_photon_evo
+      double precision eva_get_pdf_neutrino_evo
+   
+      double precision tiny,mu2min
+      double precision QW,Qf
+
+      include 'ElectroweakFlux.inc'
+      
+      tiny  = 1d-8
+      mu2min = 1d2 ! (10 GeV)^2 reset mu2min by vPID  
+
+
+c     do the following checks before calling PDF:
+c     1. momentum fraction, x
+c     2. fermion polarization fraction, fLpol
+c     3. vector boson (or neutrino) polarization by PID, vpol vPID
+c     4. evolution scale, mu2
+c     5. QED conservation check
+c     start checks
+c     1. check momentum fraction
+      if(x.lt.tiny.or.x.gt.(1d0-tiny)) then
+         write(*,*) 'eva: x out of range',x
+         eva_get_pdf_by_PID = 0d0
+         return
+      endif
+c     2. check fermion polarization fraction
+      if(fLpol.lt.0d0.or.fLpol.gt.1d0) then
+         write(*,*) 'eva: fLpol out of range',fLpol
+         stop
+         eva_get_pdf_by_PID = 0d0
+         return
+      endif
+c     3. check vector boson (or neutrino) polarization by PID
+c     also set lower bound on muf2 scale evolution by PID
+      select case (iabs(vPID))
+      case (12,14) ! ve, ve~, vm, vm~
+         mu2min = eva_mw2 ! scale set by W emission
+         if(iabs(vPol).ne.1) then
+            write(*,*) 'vPol out of range for ve/vm',vPol
+            stop 1214
+            eva_get_pdf_by_PID = 0d0
+            return
+         endif
+      case (23) ! z
+         mu2min = eva_mz2
+         if(iabs(vPol).ne.1.and.vPol.ne.0) then
+            write(*,*) 'vPol out of range for Z',vPol
+            stop 23
+            eva_get_pdf_by_PID = 0d0
+            return
+         endif
+      case (24) ! w
+         mu2min = eva_mw2
+         if(iabs(vPol).ne.1.and.vPol.ne.0) then
+            write(*,*) 'vPol out of range for W',vPol
+            stop 24
+            eva_get_pdf_by_PID = 0d0
+            return
+         endif
+      case (7,22) ! photon (special treatment for mu2min)
+         call eva_get_mf2_by_PID(mu2min,fPID) ! set scale to mass of parent fermion
+         if(iabs(vPol).ne.1) then
+            write(*,*) 'vPol out of range for A',vPol
+            stop 25
+            eva_get_pdf_by_PID = 0d0
+            return
+         endif         
+c      case (32) (eva for bsm)
+c         mu2min = eva_mx2
+c         if(iabs(vPol).ne.1.and.vPol.ne.0) then
+c            write(*,*) 'vPol out of range',vPol
+c            stop 26
+c            eva_get_pdf_by_PID = 0d0
+c            return
+c         endif
+      case default         
+         write(*,*) 'vPID out of range',vPID
+         stop 27
+         eva_get_pdf_by_PID = 0d0
+         return
+      end select
+c     4. check evolution scale
+      if(ievo.ne.0) then
+         mu2min = (1.d0-x)*mu2min
+      endif
+      if(mu2.lt.mu2min) then
+         write(*,*) 'muf2 too small. setting muf2 to muf2min:',mu2,mu2min
+         mu2 = mu2min
+      endif
+c     5. QED conservation check
+      if(iabs(vPID).eq.24) then
+         QW = dble(vPID/iabs(vPID))
+         call eva_get_qEM_by_PID(Qf,fPID)
+         if(dabs(Qf-QW).gt.eva_one) then
+            write(*,*) 'Stopping EVA: QED charge violation with emission of vPID=',vPID,' by fPID =',fPID
+            stop 24
+         return
+         endif
+      endif
+      if(iabs(vPID).eq.12.or.iabs(vPID).eq.14) then
+         select case(vPID)
+         case (12)
+            if(fPID.ne.11) then
+               write(*,*) 'Stopping EVA: neutrino mismatch with emission of vPID=',vPID,' by fPID =',fPID
+               stop 1211
+            endif
+         case (-12)
+            if(fPID.ne.-11) then
+               write(*,*) 'Stopping EVA: neutrino mismatch with emission of vPID=',vPID,' by fPID =',fPID
+               stop -1211
+            endif
+         case (14)
+            if(fPID.ne.13) then
+               write(*,*) 'Stopping EVA: neutrino mismatch with emission of vPID=',vPID,' by fPID =',fPID
+               stop 1413
+            endif
+         case (-14)
+            if(fPID.ne.-13) then
+               write(*,*) 'Stopping EVA: neutrino mismatch with emission of vPID=',vPID,' by fPID =',fPID
+               stop -1413
+            endif
+         case default
+            write(*,*) 'Stopping EVA at neutrino check. should not be here with emission of vPID=',vPID,' by fPID =',fPID
+               stop -1412
+         end select
+      endif
+c      if(iabs(vPID).eq.22.and.(
+c     &      iabs(fPID).eq.12.or.
+c     &      iabs(fPID).eq.14.or.
+c     &      iabs(fPID).eq.16)) then
+c            write(*,*) 'QED charge violation with a emission by neutrino'
+c            eva_get_pdf_by_PID = 0d0
+c         return
+c      endif
+c     celebrate by calling the PDF
+c      if(vPID.eq.22.or.vPID.eq.7) then
+c         eva_get_pdf_by_PID = eva_get_pdf_photon_evo(vPID,fPID,vpol,fLpol,x,mu2,ievo)
+c      else
+c         eva_get_pdf_by_PID = eva_get_pdf_by_PID_evo(vPID,fPID,vpol,fLpol,x,mu2,ievo)
+c      endif
+
+      select case (abs(vPID))
+      case (7,22)
+         eva_get_pdf_by_PID = eva_get_pdf_photon_evo(vPID,fPID,vpol,fLpol,x,mu2,ievo)
+      case (12,14)
+         eva_get_pdf_by_PID = eva_get_pdf_neutrino_evo(vPID,fPID,vpol,fLpol,x,mu2,ievo)
+      case default
+         eva_get_pdf_by_PID = eva_get_pdf_by_PID_evo(vPID,fPID,vpol,fLpol,x,mu2,ievo)
+      end select
+      return
+      end      
+c     /* ********************************************************* *
+c     /* ********************************************************* *
+c     /* ********************************************************* *
+c     /* ********************************************************* *
+      double precision function eva_get_pdf_by_PID_evo(vPID,fPID,vpol,fLpol,x,mu2,ievo)
+      implicit none
+      integer vPID,fPID,vpol,ievo
+      double precision fLpol,x,mu2
+      double precision eva_fX_to_vm,eva_fX_to_v0,eva_fX_to_vp
+      
+      double precision gg2,gL2,gR2,mv2,tmpPDF
+      call eva_get_mv2_by_PID(mv2,vPID)
+      call eva_get_gg2_by_PID(gg2,vPID,fPID)
+      if( fPID/iabs(fPID).gt.0 ) then ! particle
+         call eva_get_gR2_by_PID(gR2,vPID,fPID)
+         call eva_get_gL2_by_PID(gL2,vPID,fPID)
+      else  ! antiparticle (invert parity)
+         call eva_get_gR2_by_PID(gL2,vPID,fPID)
+         call eva_get_gL2_by_PID(gR2,vPID,fPID)
+      endif
+      select case (vpol)
+      case (-1)
+         tmpPDF = eva_fX_to_vm(gg2,gL2,gR2,fLpol,mv2,x,mu2,ievo)
+      case (0)
+         tmpPDF = eva_fX_to_v0(gg2,gL2,gR2,fLpol,mv2,x,mu2,ievo)
+      case (+1)
+         tmpPDF = eva_fX_to_vp(gg2,gL2,gR2,fLpol,mv2,x,mu2,ievo)
+      case default
+         write(*,*) 'vPol out of range; should not be here',vPol
+         stop
+         tmpPDF = 0d0
+      end select
+      eva_get_pdf_by_PID_evo = tmpPDF
+      return
+      end
+c     /* ********************************************************* *
+      double precision function eva_get_pdf_photon_evo(vPID,fPID,vpol,fLpol,x,mu2,ievo)
+      implicit none
+      integer vPID,fPID,vpol,ievo
+      double precision fLpol,x,mu2
+      double precision eva_fX_to_vm,eva_fX_to_v0,eva_fX_to_vp
+      
+      double precision gg2,gL2,gR2,mf2,tmpPDF
+      call eva_get_mf2_by_PID(mf2,fPID)
+      call eva_get_gg2_by_PID(gg2,vPID,fPID)
+      if( fPID/iabs(fPID).gt.0 ) then ! particle
+         call eva_get_gR2_by_PID(gR2,vPID,fPID)
+         call eva_get_gL2_by_PID(gL2,vPID,fPID)
+      else  ! antiparticle (invert parity)
+         call eva_get_gR2_by_PID(gL2,vPID,fPID)
+         call eva_get_gL2_by_PID(gR2,vPID,fPID)
+      endif
+      select case (vpol)
+      case (-1)
+         tmpPDF = eva_fX_to_vm(gg2,gL2,gR2,fLpol,mf2,x,mu2,ievo)
+      case (+1)
+         tmpPDF = eva_fX_to_vp(gg2,gL2,gR2,fLpol,mf2,x,mu2,ievo)
+      case default
+         write(*,*) 'vPol out of range; should not be here',vPol
+         stop
+         tmpPDF = 0d0         
+      end select
+      eva_get_pdf_photon_evo = tmpPDF
+      return
+      end      
+c     /* ********************************************************* *   
+c     /* ********************************************************* *
+      double precision function eva_get_pdf_neutrino_evo(vPID,fPID,vpol,fLpol,x,mu2,ievo)
+      implicit none
+      integer vPID,fPID,vpol,ievo
+      logical isAntiNu
+      double precision fLpol,x,mu2
+      double precision eva_fX_to_fR,eva_fX_to_fL
+
+      double precision gg2,gL2,gR2,mv2,tmpPDF
+      call eva_get_mv2_by_PID(mv2,vPID)
+      call eva_get_gg2_by_PID(gg2,vPID,fPID)
+      if( fPID/iabs(fPID).gt.0 ) then ! particle
+         isAntiNu = .false.
+         call eva_get_gR2_by_PID(gR2,vPID,fPID)
+         call eva_get_gL2_by_PID(gL2,vPID,fPID)
+      else  ! antiparticle (invert parity)
+         isAntiNu = .true.
+         call eva_get_gR2_by_PID(gL2,vPID,fPID)
+         call eva_get_gL2_by_PID(gR2,vPID,fPID)
+      endif
+
+      select case (vpol)
+      case (-1)
+         if(isAntiNu) then ! no LH antineutrinos
+            tmpPDF = 0
+         else  
+            tmpPDF = eva_fX_to_fL(gg2,gL2,gR2,fLpol,mv2,x,mu2,ievo)
+         endif
+      case (+1)
+         if(isAntiNu) then ! no RH neutrinos
+            tmpPDF = eva_fX_to_fR(gg2,gL2,gR2,fLpol,mv2,x,mu2,ievo)
+         else
+            tmpPDF = 0
+         endif
+      case default
+         write(*,*) 'vPol out of range; should not be here',vPol
+         stop
+         tmpPDF = 0d0         
+      end select
+      eva_get_pdf_neutrino_evo = tmpPDF
+      return
+      end      
+c     /* ********************************************************* *         
+c     /* ********************************************************* *
+c     /* ********************************************************* *
+c     /* ********************************************************* *
+c     /* ********************************************************* *                  
+      subroutine eva_get_mv2_by_PID(mv2,vPID)
+      implicit none
+      integer vPID
+      double precision mv2
+      include 'ElectroweakFlux.inc'
+      
+      select case (iabs(vPID))
+      case (7,22)
+         mv2 = eva_zero
+      case (23)
+         mv2 = eva_mz2
+      case (24)
+         mv2 = eva_mw2
+      case (12,14,16) ! l > vl splitting
+         mv2 = eva_mw2         
+c      case (25)
+c         mv2 = eva_mh2
+c      case (32)
+c         mv2 = eva_mx2
+      case default
+         write(*,*) 'eva: setting m_v to m_w. unknown vPID:', vPID
+         mv2 = eva_mw2
+      end select
+      return
+      end
+c     /* ********************************************************* *
+c     /* ********************************************************* *                  
+      subroutine eva_get_mf2_by_PID(mf2,fPID)
+      implicit none
+      integer fPID
+      double precision mf2
+      include 'ElectroweakFlux.inc'
+      
+      select case (iabs(fPID))
+      case (1)
+         mf2 = eva_md2
+      case (2)
+         mf2 = eva_mu2
+      case (3)
+         mf2 = eva_ms2
+      case (4)
+         mf2 = eva_mc2
+      case (5)
+         mf2 = eva_mb2
+      case (6)
+         mf2 = eva_mt2
+      case (11)
+         mf2 = eva_me2
+      case (12,14,16)
+         mf2 = eva_zero
+      case (13)
+         mf2 = eva_mm2
+      case (15)
+         mf2 = eva_ml2  
+      case default
+         write(*,*) 'eva: asking for mass of unknown fPID: ', fPID
+         stop 25
+         mf2 = eva_zero
+      end select
+      return
+      end
+c     /* ********************************************************* *            
+c     /* ********************************************************* *      
+      subroutine eva_get_gg2_by_PID(gg2,vPID,fPID)
+      implicit none
+      integer vPID,fPID
+      double precision gg2
+      include 'ElectroweakFlux.inc'
+      
+      select case (iabs(vPID))
+c     ******************************                                       
+      case (12,14) ! ve/vm/ve~/vm~
+         gg2 = eva_gw2/2.d0         
+c     ******************************                                                
+      case (7,22)  ! a
+c     ******************************
+         select case (iabs(fPID)) ! nested select case
+         case (1)               ! down
+            gg2 = eva_ee2*eva_qed2 ! = e^2 * (-1/3)^2
+         case (2)               ! up
+            gg2 = eva_ee2*eva_qeu2
+         case (3)               ! strange
+            gg2 = eva_ee2*eva_qed2
+         case (4)               ! charm
+            gg2 = eva_ee2*eva_qeu2
+         case (5)               ! bottom
+            gg2 = eva_ee2*eva_qed2
+         case (6)               ! top
+            gg2 = eva_ee2*eva_qeu2
+         case (11,13,15)        ! electron/muon/tau
+            gg2 = eva_ee2*eva_qee2
+         case (12,14,16)      ! electron/muon/tau-neutrino
+c            write(*,*) 'eva: nu has zero QED charge.'
+            gg2 = eva_zero
+         case default
+            write(*,*) 'eva: setting QED coup to (e*Q_e). unknown fPID:', fPID
+            gg2 = eva_ee2*eva_qee2
+         end select
+c     ******************************                     
+      case (23) ! z
+         gg2 = eva_gz2
+c     ******************************                              
+      case (24) ! w+/w-
+         gg2 = eva_gw2/2.d0
+         if(vPID.eq.24) then ! w+
+            select case (fPID)
+            case (-1,2,-3,4,-5,6,-11,12,-13,14,-15,16)
+               gg2 = gg2
+            case default
+               write(*,*) 'eva: violation of QED conservation. setting w+ffbar coup to zero'
+               gg2 = eva_zero
+            end select
+         else ! w-
+            select case (fPID)
+            case (1,-2,3,-4,5,-6,11,-12,13,-14,15,-16)
+               gg2 = gg2
+            case default
+               write(*,*) 'eva: violation of QED conservation. setting w-ffbar coup to zero'
+               gg2 = eva_zero
+            end select
+         endif 
+c     ******************************                              
+      case default
+         write(*,*) 'eva: setting coup to zero. unknown vPID:', vPID
+         gg2 = eva_zero
+      end select
+      return
+      end
+c     /* ********************************************************* * 
+c     /* ********************************************************* *      
+      subroutine eva_get_qEM_by_PID(qEM,fPID)
+      implicit none
+      integer fPID
+      double precision qEM
+      include 'ElectroweakFlux.inc'
+
+      select case (iabs(fPID)) ! nested select case
+      case (1)               ! down
+         qEM = eva_qed * fPID/iabs(fPID)
+      case (2)               ! up
+         qEM = eva_qeu * fPID/iabs(fPID)
+      case (3)               ! strange
+         qEM = eva_qed * fPID/iabs(fPID)
+      case (4)               ! charm
+         qEM = eva_qeu * fPID/iabs(fPID)
+      case (5)               ! bottom
+         qEM = eva_qed * fPID/iabs(fPID)
+      case (6)               ! top
+         qEM = eva_qeu * fPID/iabs(fPID)
+      case (11)              ! electron
+         qEM = eva_qee * fPID/iabs(fPID)
+      case (12)              ! electron-neutrino
+         qEM = eva_zero
+      case (13)              ! muon
+         qEM = eva_qee * fPID/iabs(fPID)
+      case (14)              ! muon-neutrino
+         qEM = eva_zero
+      case (15)              ! tau
+         qEM = eva_qee * fPID/iabs(fPID)
+      case (16)              ! tau-neutrino
+         qEM = eva_zero
+      case default
+         write(*,*) 'eva: setting QED charge to zero. unknown fPID:', fPID
+         qEM = eva_zero
+      end select
+c     ******************************                     
+      return
+      end
+c     /* ********************************************************* *      
+c     /* ********************************************************* *      
+      subroutine eva_get_gR2_by_PID(gR2,vPID,fPID)
+      implicit none
+      integer vPID,fPID
+      double precision gR2
+      include 'ElectroweakFlux.inc'
+      
+      select case (iabs(vPID))
+      case (7,22)
+         gR2 = eva_one
+      case (23)
+c     ******************************
+         select case (iabs(fPID)) ! nested select case
+         case (1)               ! down
+            gR2 = eva_zRd**2
+         case (2)               ! up
+            gR2 = eva_zRu**2
+         case (3)               ! strange
+            gR2 = eva_zRd**2
+         case (4)               ! charm
+            gR2 = eva_zRu**2
+         case (5)               ! bottom
+            gR2 = eva_zRd**2
+         case (6)               ! top
+            gR2 = eva_zRu**2            
+         case (11)              ! electron
+            gR2 = eva_zRe**2
+         case (12)              ! electron-neutrino
+            gR2 = eva_zRv**2
+         case (13)              ! muon
+            gR2 = eva_zRe**2
+         case (14)              ! muon-neutrino
+            gR2 = eva_zRv**2
+         case (15)              ! tau
+            gR2 = eva_zRe**2
+         case (16)              ! tau-neutrino
+            gR2 = eva_zRv**2                        
+         case default
+            gR2 = eva_one**2
+         end select
+c     ******************************            
+      case (24)
+         gR2 = eva_zero
+      case default
+         gR2 = eva_one
+      end select
+      return
+      end
+c     /* ********************************************************* *      
+c     /* ********************************************************* *      
+      subroutine eva_get_gL2_by_PID(gL2,vPID,fPID)
+      implicit none
+      integer vPID,fPID
+      double precision gL2
+      include 'ElectroweakFlux.inc'
+      
+      select case (iabs(vPID))
+      case (7,22)
+         gL2 = eva_one
+      case (23)
+c     ******************************
+         select case (iabs(fPID)) ! nested select case
+         case (1)               ! down
+            gL2 = eva_zLd**2
+         case (2)               ! up
+            gL2 = eva_zLu**2
+         case (3)               ! strange
+            gL2 = eva_zLd**2
+         case (4)               ! charm
+            gL2 = eva_zLu**2
+         case (5)               ! bottom
+            gL2 = eva_zLd**2
+         case (6)               ! top
+            gL2 = eva_zLu**2            
+         case (11)              ! electron
+            gL2 = eva_zLe**2
+         case (12)              ! electron-neutrino
+            gL2 = eva_zLv**2
+         case (13)              ! muon
+            gL2 = eva_zLe**2
+         case (14)              ! muon-neutrino
+            gL2 = eva_zLv**2
+         case (15)              ! tau
+            gL2 = eva_zLe**2
+         case (16)              ! tau-neutrino
+            gL2 = eva_zLv**2                        
+         case default
+            gL2 = eva_one**2
+         end select
+c     ******************************            
+      case (24)
+         gL2 = eva_one
+      case default
+         gL2 = eva_one
+      end select
+      return
+      end
+c     /* ********************************************************* *      
 
=== modified file 'Template/LO/Source/PDF/PhotonFlux.f'
--- Template/LO/Source/PDF/PhotonFlux.f	2020-08-07 14:49:25 +0000
+++ Template/LO/Source/PDF/PhotonFlux.f	2021-11-12 09:03:34 +0000
@@ -1,9 +1,14 @@
-c/* ********************************************************* */
+c/* ******************************************************** */
 c/*  Equivalent photon approximation structure function.   * */
-c/*   V.M.Budnev et al., Phys.Rep. 15C (1975) 181          * */
-c/*     Improved Weizsaecker-Williams formula              * */
-c/*     http://inspirehep.net/record/359425                * */  
-c/* ********************************************************* */
+c/*   V.M.Budnev, et al., Phys.Rep. 15C (1975) 181         * */
+c/*   https://inspirehep.net/literature/95445              * */
+c/*   For photon from proton (elastic limit, no DGLAP evo) * */
+c/*                                                        * */
+c/*  Improved Weizsaecker-Williams formula                 * */
+c/*   S. Frixione, et al., Phys.Lett.B 319 (1993) 339-345  * */
+c/*   http://inspirehep.net/record/359425                  * */
+c/*   For photon from electron (WW approx) + higher twist  * */
+c/* ******************************************************** */
 c   provided by Tomasz Pierzchala - UCL
 
       real*8 function epa_lepton(x,q2max, mode)
 
=== modified file 'Template/LO/Source/PDF/makefile'
--- Template/LO/Source/PDF/makefile	2021-04-16 07:24:34 +0000
+++ Template/LO/Source/PDF/makefile	2021-11-12 09:03:34 +0000
@@ -15,14 +15,14 @@
     $(error Bad lhadpfversion version 6 is now required)
   else
     ifeq ($(lhapdfsubversion),1) # 6.1.X
-	PDF         = pdfwrap_lhapdf.o pdf_lhapdf6.o pdg2pdf_lhapdf6.o opendata.o PhotonFlux.o
+	PDF         = pdfwrap_lhapdf.o pdf_lhapdf6.o pdg2pdf_lhapdf6.o opendata.o PhotonFlux.o ElectroweakFluxDriver.o ElectroweakFlux.o
     else # 6.2.X
         CXXFLAGS+=-std=c++11
-	PDF         = pdfwrap_lhapdf.o pdf_lhapdf62.o pdg2pdf_lhapdf6.o opendata.o PhotonFlux.o
+	PDF         = pdfwrap_lhapdf.o pdf_lhapdf62.o pdg2pdf_lhapdf6.o opendata.o PhotonFlux.o ElectroweakFluxDriver.o ElectroweakFlux.o
     endif    
   endif
 else
-  PDF         = Ctq6Pdf.o pdfwrap.o opendata.o pdf.o PhotonFlux.o pdg2pdf.o NNPDFDriver.o eepdf.o gridpdfaux.o dfint.o kerset.o
+  PDF         = Ctq6Pdf.o pdfwrap.o opendata.o pdf.o PhotonFlux.o pdg2pdf.o NNPDFDriver.o eepdf.o gridpdfaux.o dfint.o kerset.o ElectroweakFluxDriver.o ElectroweakFlux.o
 endif
 
 all: $(LIBDIR)$(LIBRARY)
 
=== modified file 'Template/LO/Source/PDF/pdfwrap.f'
--- Template/LO/Source/PDF/pdfwrap.f	2021-04-19 12:59:10 +0000
+++ Template/LO/Source/PDF/pdfwrap.f	2021-11-12 09:03:34 +0000
@@ -261,8 +261,8 @@
 
 c---------------------------------------------------------------
 c---------------------------------------------------------------
-
-
+      elseif (pdlabel.eq.'eva'.or.pdlabel.eq.'iww'.or.pdlabel.eq.'none') then
+            asmz=asmz            
       else
          asmz=0.118d0
 c     do not crash anymore since arbitrary PDF can be added for lepton collision
 
=== modified file 'Template/LO/Source/PDF/pdg2pdf.f'
--- Template/LO/Source/PDF/pdg2pdf.f	2021-04-24 20:53:56 +0000
+++ Template/LO/Source/PDF/pdg2pdf.f	2021-11-12 09:03:34 +0000
@@ -52,12 +52,32 @@
       data pdlabellast/2*'abcdefg'/
       data ihlast/2*-99/
 
-      if (ih.eq.9) then
+c     effective w/z/a approximation (leading log fixed order, not resummed)
+      double precision eva_get_pdf_by_PID
+      external eva_get_pdf_by_PID
+      integer ppid
+      integer ievo,ievo_eva
+      common/to_eva/ievo_eva
+      integer hel,helMulti,hel_picked
+      double precision hel_jacobian
+      common/hel_picked/hel_picked,hel_jacobian
+      integer get_nhel
+      external get_nhel
+      real*8 pol(2),fLPol
+      common/to_polarization/pol
+
+
+c     collider configuration
+      integer lpp(2)
+      double precision ebeam(2),xbk(2),q2fact(2)
+      common/to_collider/ebeam,xbk,q2fact,lpp
+
+      if (iabs(ih).eq.9) then
          pdg2pdf = 1d0
          return
       endif
 
-      nb_hadron = (nb_proton(beamid)+nb_neutron(beamid))
+      nb_hadron = (nb_proton(iabs(beamid))+nb_neutron(iabs(beamid)))
 c     Make sure we have a reasonable Bjorken x. Note that even though
 c     x=0 is not reasonable, we prefer to simply return pdg2pdf=0
 c     instead of stopping the code, as this might accidentally happen.
@@ -101,31 +121,48 @@
            stop 1
         endif
         do i_ee = 1, n_ee
-          ee_components(i_ee) = compute_eepdf(x,omx_ee(beamid),xmu,i_ee,ipart,ih_local)
+          ee_components(i_ee) = compute_eepdf(x,omx_ee(iabs(beamid)),xmu,i_ee,ipart,ih_local)
 	enddo
         pdg2pdf =  ee_components(1) ! temporary to test pdf load
-c        write(*,*), x, beamid ,omx_ee(beamid),xmu,1,ipart,ih_local,pdg2pdf
+c        write(*,*), x, beamid ,omx_ee(iabs(beamid)),xmu,1,ipart,ih_local,pdg2pdf
         return
       endif
       
 
-      if (beamid.gt.0) then
+c     If group_subprocesses is true, then IH=abs(lpp) and ipdg=ipdg*sgn(lpp) in export_v4.
+c     For EVA,  group_subprocesses is false and IH=LPP and ipdg are passed, instead.
+c     If group_subprocesses is false, the following sets ipdg=ipdg*sgn(IH) if not in EVA
+      if(pdlabel.eq.'eva'.or.pdsublabel(iabs(beamid)).eq.'eva') then
+         ipart=ipdg
+      else
          ipart=sign(1,ih)*ipdg
-      else
-         ipart = ipdg
-      endif
+      endif      
 
-      if(iabs(ipart).eq.21) then
+      if(iabs(ipart).eq.21) then ! g      
          ipart=0
-      else if(iabs(ipart).eq.22) then
-         ipart=7
-      else if(iabs(ipart).eq.7) then
-         ipart=7
-c     This will be called for any PDG code, but we only support up to 7
+c      else if(ipart.eq.12) then ! ve
+c         ipart=12
+c      else if(ipart.eq.-12) then ! ve~
+c         ipart=-12
+c      else if(ipart.eq.14) then ! vm
+c         ipart=14
+c      else if(ipart.eq.-14) then ! vm~
+c         ipart=-14
+      else if(ipart.eq.24) then  ! w+
+         ipart=24
+      else if(ipart.eq.-24) then ! w-
+         ipart=-24
+      else if(iabs(ipart).eq.23) then ! z
+         ipart=23
+      else if(iabs(ipart).eq.22) then ! a
+         ipart=7
+      else if(iabs(ipart).eq.7) then  ! a
+         ipart=7
+c     This will be called for any PDG code. We only support (for now) 0-7, and 22-24
       else if(iabs(ipart).gt.7)then
          write(*,*) 'PDF not supported for pdg ',ipdg
          write(*,*) 'For lepton colliders, please set the lpp* '//
-     $    'variables to 0 in the run_card'  
+     $    'variables to 0 in the run_card current is' , ih  
          open(unit=26,file='../../../error',status='unknown')
          write(26,*) 'Error: PDF not supported for pdg ',ipdg
          stop 1
@@ -142,10 +179,10 @@
       enddo
 
 c     Reuse previous result, if possible
-      if (ireuse.gt.0.)then
+      if (ireuse.gt.0.and.iabs(iporg).lt.8)then
          if (pdflast(iporg,ireuse).ne.-99d9) then
-            pdg2pdf = get_ion_pdf(pdflast(-7, ireuse), iporg, nb_proton(beamid),
-     $                         nb_neutron(beamid))
+            pdg2pdf = get_ion_pdf(pdflast(-7, ireuse), iporg, nb_proton(iabs(beamid)),
+     $                         nb_neutron(iabs(beamid)))
             return 
          endif
       endif
@@ -186,31 +223,67 @@
       pdlabellast(ireuse)=pdlabel
       ihlast(ireuse)=ih
 
-      if(iabs(ipart).eq.7.and.ih.gt.1) then
-         q2max=xmu*xmu
-         if(abs(ih).eq.3.or.abs(ih).eq.4) then       !from the electron or muonn
-            pdg2pdf=epa_lepton(x,q2max, ih)
-         elseif(ih .eq. 2) then !from a proton without breaking
-            pdg2pdf=epa_proton(x,q2max,beamid)
-         endif 
-         pdflast(iporg,ireuse)=pdg2pdf
-         return
+      if(pdlabel.eq.'eva'.or.pdsublabel(iabs(beamid)).eq.'eva') then
+         if(iabs(ipart).ne.7.and.
+c     &      iabs(ipart).ne.12.and.
+c     &      iabs(ipart).ne.14.and.     
+     &      iabs(ipart).ne.23.and.
+     &      iabs(ipart).ne.24 ) then
+            write(*,*) 'ERROR: EVA PDF only supported for A/Z/W, not for pdg = ',ipart
+            stop 1
+         else
+c         write(*,*) 'running eva'
+            select case (iabs(ih))
+            case (0:2)
+               write(*,*) 'ERROR: EVA PDF only supported for e+/- and mu+/- beams, not for lpp/ih=',ih
+               stop 24
+            case (3) ! e+/-
+               ppid = 11
+            case (4) ! mu+/-
+               ppid = 13
+            case default
+               write(*,*) 'ERROR: EVA PDF only supported for e+/- and mu+/- beams, not for lpp/ih=',ih
+               stop 24
+            end select
+            ppid  = ppid * ih/iabs(ih) ! get sign of parent
+            fLPol = pol(iabs(beamid))        ! see setrun.f for treatment of polbeam*
+            q2max = xmu*xmu
+            ievo = ievo_eva
+            hel      = GET_NHEL(HEL_PICKED, beamid) ! helicity of v
+            helMulti = GET_NHEL(0, beamid)          ! helicity multiplicity of v to undo spin averaging
+            pdg2pdf  = helMulti*eva_get_pdf_by_PID(ipart,ppid,hel,fLpol,x,q2max,ievo)
+            return
+         endif
+      else ! this ensure backwards compatibility
+         if(iabs(ipart).eq.7.and.iabs(ih).gt.1) then
+            q2max=xmu*xmu
+            if(iabs(ih).eq.3.or.iabs(ih).eq.4) then       !from the electron or muonn
+               pdg2pdf=epa_lepton(x,q2max, iabs(ih))
+            elseif(iabs(ih) .eq. 2) then !from a proton without breaking
+               pdg2pdf=epa_proton(x,q2max,beamid)
+            endif 
+            pdflast(iporg,ireuse)=pdg2pdf
+            return
+         endif         
       endif
+
+
+
       
       if (pdlabel(1:5) .eq. 'cteq6') then
 C        Be carefull u and d are flipped inside cteq6
-         if (nb_proton(beamid).gt.1.or.nb_neutron(beamid).ne.0)then
+         if (nb_proton(iabs(beamid)).gt.1.or.nb_neutron(iabs(beamid)).ne.0)then
             if (ipart.eq.1.or.ipart.eq.2)then
                pdflast(1,ireuse)=Ctq6Pdf(2,x*nb_hadron,xmu) ! remember u/d flipping in cteq
                pdflast(2,ireuse)=Ctq6Pdf(1,x*nb_hadron,xmu)
-               pdg2pdf = get_ion_pdf(pdflast(-7,ireuse), ipart, nb_proton(beamid), nb_neutron(beamid))
+               pdg2pdf = get_ion_pdf(pdflast(-7,ireuse), ipart, nb_proton(iabs(beamid)), nb_neutron(iabs(beamid)))
             else if (ipart.eq.-1.or.ipart.eq.-2)then
                pdflast(-1,ireuse)=Ctq6Pdf(-2,x*nb_hadron,xmu) ! remember u/d flipping in cteq
                pdflast(-2,ireuse)=Ctq6Pdf(-1,x*nb_hadron,xmu)
-               pdg2pdf = get_ion_pdf(pdflast(-7,ireuse), ipart, nb_proton(beamid), nb_neutron(beamid))
+               pdg2pdf = get_ion_pdf(pdflast(-7,ireuse), ipart, nb_proton(iabs(beamid)), nb_neutron(iabs(beamid)))
             else
                pdflast(ipart,ireuse)=Ctq6Pdf(ipart,x*nb_hadron,xmu)
-               pdg2pdf = get_ion_pdf(pdflast(-7,ireuse), ipart, nb_proton(beamid), nb_neutron(beamid))
+               pdg2pdf = get_ion_pdf(pdflast(-7,ireuse), ipart, nb_proton(iabs(beamid)), nb_neutron(iabs(beamid)))
             endif 
          else
             if(iabs(ipart).ge.1.and.iabs(ipart).le.2)
@@ -219,9 +292,9 @@
             pdflast(iporg,ireuse)=pdg2pdf
          endif
       else
-         call pftopdg(ih,x*nb_hadron,xmu,pdflast(-7,ireuse))
-         pdg2pdf = get_ion_pdf(pdflast(-7, ireuse), iporg, nb_proton(beamid),
-     $                         nb_neutron(beamid))
+         call pftopdg(iabs(ih),x*nb_hadron,xmu,pdflast(-7,ireuse))
+         pdg2pdf = get_ion_pdf(pdflast(-7, ireuse), iporg, nb_proton(iabs(beamid)),
+     $                         nb_neutron(iabs(beamid)))
       endif      
 
       return
 
=== modified file 'Template/LO/Source/PDF/pdg2pdf_lhapdf6.f'
--- Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	2021-07-28 14:12:54 +0000
+++ Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	2021-11-12 09:03:34 +0000
@@ -34,7 +34,21 @@
       data imemlast/20*-99/
       data i_replace/20/
 
-      nb_hadron = (nb_proton(beamid)+nb_neutron(beamid))
+c     effective w/z/a approximation (leading log fixed order, not resummed)
+      double precision eva_get_pdf_by_PID
+      external eva_get_pdf_by_PID
+      integer ppid
+      integer ievo,ievo_eva
+      common/to_eva/ievo_eva
+      integer hel,helMulti,hel_picked
+      double precision hel_jacobian
+      common/hel_picked/hel_picked,hel_jacobian
+      integer get_nhel
+      external get_nhel
+      real*8 pol(2),fLPol
+      common/to_polarization/pol
+
+      nb_hadron = (nb_proton(iabs(beamid))+nb_neutron(iabs(beamid)))
 c     Make sure we have a reasonable Bjorken x. Note that even though
 c     x=0 is not reasonable, we prefer to simply return pdg2pdf=0
 c     instead of stopping the code, as this might accidentally happen.
@@ -53,21 +67,85 @@
          endif
       endif
 
-      ipart=ipdg
-      if(iabs(ipart).eq.21) then
+c     If group_subprocesses is true, then IH=abs(lpp) and ipdg=ipdg*sgn(lpp) in export_v4.
+c     For EVA,  group_subprocesses is false and IH=LPP and ipdg are passed, instead.
+c     If group_subprocesses is false, the following sets ipdg=ipdg*sgn(IH) if not in EVA
+      if(pdsublabel(iabs(beamid)).eq.'eva') then
+         ipart=ipdg
+      else 
+         ipart=ipdg*ih/iabs(ih)
+      endif    
+
+      if(iabs(ipart).eq.21) then ! g
          ipart=0
-      else if(iabs(ipart).eq.22) then 
-         ipart=7
-      else if (iabs(ipart).eq.7) then
-         ipart=7
-c     This will be called for any PDG code, but we only support up to 7
-      else if(iabs(ipart).gt.7)then
-         write(*,*) 'PDF not supported for pdg ',ipdg
-         write(*,*) 'For lepton colliders, please set the lpp* '//
-     $    'variables to 0 in the run_card'  
-         open(unit=26,file='../../../error',status='unknown')
-         write(26,*) 'Error: PDF not supported for pdg ',ipdg
-         stop 1
+c      else if(ipart.eq.12) then ! ve
+c         ipart=12
+c      else if(ipart.eq.-12) then ! ve~
+c         ipart=-12
+c      else if(ipart.eq.14) then ! vm
+c         ipart=14
+c      else if(ipart.eq.-14) then ! vm~
+c         ipart=-14   
+      else if(ipart.eq.24) then  ! w+
+         ipart=24
+      else if(ipart.eq.-24) then ! w-
+         ipart=-24
+      else if(iabs(ipart).eq.23) then ! z
+         ipart=23
+      else if(iabs(ipart).eq.22) then ! a
+         ipart=7
+      else if(iabs(ipart).eq.7) then  ! a
+         ipart=7
+c     This will be called for any PDG code. We only support (for now) 0-7, and 22-24
+c      else if(iabs(ipart).gt.7)then
+c         write(*,*) 'PDF not supported for pdg ',ipdg
+c         write(*,*) 'For lepton colliders, please set the lpp* '//
+c     $    'variables to 0 in the run_card'  
+c         open(unit=26,file='../../../error',status='unknown')
+c         write(26,*) 'Error: PDF not supported for pdg ',ipdg
+c         stop 1
+      endif
+      
+      if(pdsublabel(iabs(beamid)).eq.'eva') then
+         if(iabs(ipart).ne.7.and.
+c     &      iabs(ipart).ne.12.and.
+c     &      iabs(ipart).ne.14.and.     
+     &      iabs(ipart).ne.23.and.
+     &      iabs(ipart).ne.24 ) then
+            write(*,*) 'ERROR: EVA PDF only supported for A/Z/W, not for pdg = ',ipart
+            stop 1
+         else
+c         write(*,*) 'running eva'
+            select case (iabs(ih))
+            case (0:2)
+               write(*,*) 'ERROR: EVA PDF only supported for e+/- and mu+/- beams, not for lpp/ih=',ih
+               stop 24
+            case (3) ! e+/-
+               ppid = 11
+            case (4) ! mu+/-
+               ppid = 13
+            case default
+               write(*,*) 'ERROR: EVA PDF only supported for e+/- and mu+/- beams, not for lpp/ih=',ih
+               stop 24
+            end select
+            ppid  = ppid * ih/iabs(ih) ! get sign of parent
+            fLPol = pol(iabs(beamid))        ! see setrun.f for treatment of polbeam*
+c              q2max = xmu*xmu
+            ievo = ievo_eva
+            hel      = GET_NHEL(HEL_PICKED, beamid) ! helicity of v
+            helMulti = GET_NHEL(0, beamid)          ! helicity multiplicity of v to undo spin averaging
+            pdg2pdf  = helMulti*eva_get_pdf_by_PID(ipart,ppid,hel,fLpol,x,xmu*xmu,ievo)
+            return
+         endif
+      else
+         if(iabs(ipart).eq.24.or.iabs(ipart).eq.23) then  ! w/z
+            write(*,*) 'LHAPDF not supported for pdg ',ipdg
+            write(*,*) 'For EVA, check if pdlabel and pdsublabel* '//
+     $    'are set correctly in the run_card'  
+            open(unit=26,file='../../../error',status='unknown')
+            write(26,*) 'Error: PDF not supported for pdg ',ipdg
+            stop 1
+         endif
       endif
 
       iporg=ipart
@@ -86,6 +164,9 @@
       ireuse = 0
       ii=i_replace
       do i=1,20
+         if (abs(ipart).gt.7)then
+            exit
+         endif
 c     Check if result can be reused since any of last twenty
 c     calls. Start checking with the last call and move back in time
          if (ih.eq.ihlast(ii)) then
@@ -105,9 +186,9 @@
       enddo
 
 c     Reuse previous result, if possible
-      if (ireuse.gt.0) then
+      if (ireuse.gt.0.and.abs(ipart).le.7) then
          if (pdflast(ipart,ireuse).ne.-99d9) then
-            pdg2pdf = get_ion_pdf(pdflast(-7,ireuse), ipart, nb_proton(beamid), nb_neutron(beamid))/x
+            pdg2pdf = get_ion_pdf(pdflast(-7,ireuse), ipart, nb_proton(iabs(beamid)), nb_neutron(iabs(beamid)))/x
             return 
          endif
       endif
@@ -117,10 +198,10 @@
 
 c     Call lhapdf and give the current values to the arrays that should
 c     be saved
-      if(ih.eq.1) then
-         if (nb_proton(beamid).eq.1.and.nb_neutron(beamid).eq.0) then
+      if(iabs(ih).eq.1) then
+         if (nb_proton(iabs(beamid)).eq.1.and.nb_neutron(iabs(beamid)).eq.0) then
             call evolvepart(ipart,x,xmu,pdg2pdf)
-            pdflast(ipart, i_replace)=pdg2pdf
+            if (abs(ipart).le.7)   pdflast(ipart, i_replace)=pdg2pdf
          else
             if (ipart.eq.1.or.ipart.eq.2) then
                call evolvepart(1,x*nb_hadron
@@ -136,12 +217,12 @@
                call evolvepart(ipart,x*nb_hadron
      &                         ,xmu,pdflast(ipart, i_replace))
             endif 
-            pdg2pdf = get_ion_pdf(pdflast(-7, i_replace), ipart, nb_proton(beamid), nb_neutron(beamid))
+            pdg2pdf = get_ion_pdf(pdflast(-7, i_replace), ipart, nb_proton(iabs(beamid)), nb_neutron(iabs(beamid)))
          endif
          pdg2pdf=pdg2pdf/x
-      else if(abs(ih).eq.3.or.abs(ih).eq.4) then       !from the electron
-            pdg2pdf=epa_lepton(x,xmu*xmu, ih)
-      else if(ih.eq.2) then ! photon from a proton without breaking
+      else if(iabs(ih).eq.3.or.iabs(ih).eq.4) then       !from the electron
+            pdg2pdf=epa_lepton(x,xmu*xmu, iabs(ih))
+      else if(iabs(ih).eq.2) then ! photon from a proton without breaking
           pdg2pdf = epa_proton(x,xmu*xmu, beamid)
 
       else
 
=== modified file 'Template/LO/Source/cuts.inc'
--- Template/LO/Source/cuts.inc	2016-02-04 19:44:31 +0000
+++ Template/LO/Source/cuts.inc	2021-11-12 09:03:34 +0000
@@ -37,6 +37,7 @@
       REAL*8 misset,missetmax,ptheavy
       REAL*8 ptllmin,ptllmax
       integer maxjetflavor
+      REAl*8 dsqrt_shat
 	
       COMMON /to_min_max_cuts/
      &     PTJmax,PTBmax,PTAmax,PTLmax,
@@ -59,7 +60,7 @@
      &     ht2max,ht3max,ht4max,
      &     htjmin,htjmax,ihtmin,ihtmax,
      &     misset,missetmax,ptheavy,
-     &     ptllmin,ptllmax,
+     &     ptllmin,ptllmax,dsqrt_shat,
      &     maxjetflavor
 
 C
 
=== modified file 'Template/LO/Source/dsample.f'
--- Template/LO/Source/dsample.f	2021-05-04 15:12:08 +0000
+++ Template/LO/Source/dsample.f	2021-11-12 09:03:34 +0000
@@ -1459,19 +1459,20 @@
       COMMON/TO_CUTSDONE/CUTSDONE,CUTSPASSED
  
       CHARACTER*7         PDLABEL,EPA_LABEL
+      character*7 pdsublabel(2)
       INTEGER       LHAID
-      COMMON/TO_PDF/LHAID,PDLABEL,EPA_LABEL
+      COMMON/TO_PDF/LHAID,PDLABEL,EPA_LABEL,pdsublabel
 c     
 c     Begin code
 c
 c       It is important to divide the wgt stored in the grid by the 
 c       corresponding jacobian otherwise it flattens the sampled
 c       distribution.
-C       Also, if HEL_PICKED is equal to -1, it means that MadEvent
+C       Also, if HEL_PICKED is greater than 0, it means that MadEvent
 C       is in the initialization stage where all helicity were probed
 c       and added individually to the grid directly by matrix<i>.f so
 c       that they shouldn't be added here.
-        if(ISUM_HEL.ne.0.and.HEL_PICKED.ne.-1.and.
+        if(ISUM_HEL.ne.0.and.HEL_PICKED.gt.0.and.
      &                            (.NOT.CUTSDONE.or.CUTSPASSED)) then
           call DS_add_entry('Helicity',HEL_PICKED,(wgt/hel_jacobian))
         endif
 
=== removed file 'Template/LO/Source/pdf.inc'
=== modified file 'Template/LO/Source/readgrid.f'
--- Template/LO/Source/readgrid.f	2012-03-22 06:16:35 +0000
+++ Template/LO/Source/readgrid.f	2021-11-12 09:03:34 +0000
@@ -4,7 +4,7 @@
 
 c...global variables
       include 'sudgrid.inc'
-      include 'pdf.inc'
+      include 'PDF/pdf.inc'
       include 'maxparticles.inc'
       include 'run.inc'
 
 
=== modified file 'Template/LO/Source/run.inc'
--- Template/LO/Source/run.inc	2021-07-27 15:27:29 +0000
+++ Template/LO/Source/run.inc	2021-11-12 09:03:34 +0000
@@ -6,10 +6,11 @@
 c
       real*8          scale,scalefact,alpsfact
       logical         fixed_ren_scale,fixed_fac_scale1,fixed_fac_scale2,fixed_couplings,hmult
-      integer         ickkw,nhmult,asrwgtflavor, dynamical_scale_choice
+      integer         ickkw,nhmult,asrwgtflavor,dynamical_scale_choice,ievo_eva
       common/to_scale/scale,scalefact,alpsfact,fixed_ren_scale,fixed_fac_scale1,fixed_fac_scale2,
      $                fixed_couplings,ickkw,nhmult,hmult,asrwgtflavor,
      $                dynamical_scale_choice
+      common/to_eva/ievo_eva
 c
 c     Collider
 c
 
=== modified file 'Template/LO/Source/setrun.f'
--- Template/LO/Source/setrun.f	2021-04-24 20:53:56 +0000
+++ Template/LO/Source/setrun.f	2021-11-12 09:03:34 +0000
@@ -106,6 +106,22 @@
      $     sign((abs(pol(1))-1)*100,pol(1)),
      $     sign((abs(pol(2))-1)*100,pol(2))
 
+
+      if(pdlabel.eq.'eva') then
+            ! pbX=-100 (pure LH beam) => fLpol=1.0 (in eva)
+            ! pbX=0    (RH + LH beam) => fLpol=0.5 (in eva)
+            ! pbX=+100 (pure RH beam) => fLpol=0.0 (in eva)
+            pol(1) = (-1d0/200d0)*pb1 + 0.5d0
+            pol(2) = (-1d0/200d0)*pb2 + 0.5d0
+      else
+            if(pdsublabel(1).eq.'eva') then
+                  pol(1) = (-1d0/200d0)*pb1 + 0.5d0
+            endif
+            if(pdsublabel(2).eq.'eva') then
+                  pol(2) = (-1d0/200d0)*pb2 + 0.5d0
+            endif
+      endif
+
 c !!! Default behavior changed (MH, Aug. 07) !!!
 c If no pdf, read the param_card and use the value from there and
 c order of alfas running = 2
@@ -132,10 +148,10 @@
 c           In principle this should be always the case since the
 c           banner.py is expected to correct such wrong run_card.
       if(use_syst)then
-         if(scalefact.ne.1)then
-            write(*,*) 'Warning: use_syst=T, setting scalefact to 1'
-            scalefact=1
-         endif
+c         if(scalefact.ne.1)then
+c            write(*,*) 'Warning: use_syst=T, setting scalefact to 1'
+c            scalefact=1
+c         endif
          if(alpsfact.ne.1)then
             write(*,*) 'Warning: use_syst=T, setting alpsfact to 1'
             alpsfact=1
@@ -152,6 +168,10 @@
           idbmup(i)=11
         elseif(lpp(i).eq.-3) then
           idbmup(i)=-11
+        elseif(lpp(i).eq.4) then
+            idbmup(i)=13
+        elseif(lpp(i).eq.-4) then
+            idbmup(i)=-13
         elseif(lpp(i).eq.0) then
           idbmup(i)=idup(i,1,1)
         else
@@ -179,10 +199,12 @@
       integer mpdf
       integer npdfs,i,pdfgup(2),pdfsup(2),lhaid
 
-      parameter (npdfs=17)
+      parameter (npdfs=19)
       character*7 pdflabs(npdfs)
       data pdflabs/
-     $     'none',
+     $   'none',
+     $   'eva',
+     $   'iww',
      $     'dressed', 
      $   'mrs02nl',
      $   'mrs02nn',
@@ -201,8 +223,10 @@
      $   'nn23nlo'/
       integer numspdf(npdfs)
       data numspdf/
-     $     00000,
-     $     00000, 
+     $   00000,
+     $   00000,
+     $   00000,
+     $   00000, 
      $   20250,
      $   20270,
      $   19150,
 
=== modified file 'Template/LO/SubProcesses/cuts.f'
--- Template/LO/SubProcesses/cuts.f	2021-07-27 15:27:29 +0000
+++ Template/LO/SubProcesses/cuts.f	2021-11-12 09:03:34 +0000
@@ -299,11 +299,12 @@
 c
 c     Limit S_hat
 c
-c      if (x1*x2*stot .gt. 500**2) then
-c         passcuts=.false.
-c         return
-c      endif
-
+      if (dsqrt_shat.ne.0d0)then
+         if (nincoming.eq.2.and.sumdot(p(0,1),p(0,2),1d0) .lt. dsqrt_shat**2) then
+            passcuts=.false.
+            return
+         endif
+      endif
 C $B$ DESACTIVATE_CUT $E$ !This is a tag for MadWeight
 
       if(debug) write (*,*) '============================='
 
=== modified file 'Template/LO/SubProcesses/genps.f'
--- Template/LO/SubProcesses/genps.f	2021-07-28 14:12:54 +0000
+++ Template/LO/SubProcesses/genps.f	2021-11-12 09:03:34 +0000
@@ -147,7 +147,8 @@
 C     Common blocks
       CHARACTER*7         PDLABEL,EPA_LABEL
       INTEGER       LHAID
-      COMMON/TO_PDF/LHAID,PDLABEL,EPA_LABEL
+      character*7 pdsublabel(2)
+      COMMON/TO_PDF/LHAID,PDLABEL,EPA_LABEL, pdsublabel
 
       
       integer          lwgt(0:maxconfigs,maxinvar)
 
=== modified file 'Template/LO/SubProcesses/setcuts.f'
--- Template/LO/SubProcesses/setcuts.f	2021-07-27 15:27:29 +0000
+++ Template/LO/SubProcesses/setcuts.f	2021-11-12 09:03:34 +0000
@@ -679,7 +679,7 @@
       do i=nincoming+1, nexternal
          smin_p = smin_p + pmass(i)
       enddo
-      smin = max(smin, smin_p**2)
+      smin = max(smin, smin_p**2, dsqrt_shat**2)
       write(*,*) "Define smin to", smin
 
 c Check that results are consistent among all the grouped subprocesses
 
=== modified file 'Template/LO/SubProcesses/unwgt.f'
--- Template/LO/SubProcesses/unwgt.f	2021-03-13 10:03:56 +0000
+++ Template/LO/SubProcesses/unwgt.f	2021-11-12 09:03:34 +0000
@@ -724,7 +724,7 @@
             beam_number =2
          endif
          
-         if(n_pdfrw(1).gt.0.and.abs(lpp(1)).eq.1)then
+         if(n_pdfrw(1).gt.0.and.abs(lpp(1)).ne.2)then
             if(2*n_pdfrw(1).lt.10) then
                write(cfmt,'(a,I1,a,I1,a)') '(a,I1,a,I3,',
      $              n_pdfrw(1),'I9,',2*n_pdfrw(1),'E15.8,a)'
@@ -746,7 +746,7 @@
          if (flip) then
             beam_number	= 1
          endif
-         if(n_pdfrw(2).gt.0.and.abs(lpp(2)).eq.1)then
+         if(n_pdfrw(2).gt.0.and.abs(lpp(2)).ne.2)then
             if(2*n_pdfrw(2).lt.10) then
                write(cfmt,'(a,I1,a,I1,a)') '(a,I1,a,I3,',
      $              n_pdfrw(2),'I9,',2*n_pdfrw(2),'E15.8,a)'
 
=== modified file 'Template/MadWeight/Cards/run_card.dat'
--- Template/MadWeight/Cards/run_card.dat	2015-04-05 23:18:48 +0000
+++ Template/MadWeight/Cards/run_card.dat	2021-11-12 09:03:34 +0000
@@ -23,8 +23,9 @@
   %(iseed)s   = iseed   ! rnd seed (0=assigned automatically=default))
 #*********************************************************************
 # Collider type and energy                                           *
-# lpp: 0=No PDF, 1=proton, -1=antiproton, 2=photon from proton,      *
-#                                         3=photon from electron     *
+# lpp: 0=No PDF, 1=proton, -1=antiproton, 2=elastic photon of proton,*
+#             +/-3=PDF of electron/positron beam                     *
+#             +/-4=PDF of muon/antimuon beam                         *
 #*********************************************************************
      %(lpp1)s        = lpp1    ! beam 1 type 
      %(lpp2)s        = lpp2    ! beam 2 type
 
=== modified file 'Template/MadWeight/src/setrun.f'
--- Template/MadWeight/src/setrun.f	2021-08-20 19:53:32 +0000
+++ Template/MadWeight/src/setrun.f	2021-11-12 09:03:34 +0000
@@ -537,3 +537,11 @@
 
       return
       end
+
+      INTEGER FUNCTION GET_NHEL(HEL, IPART)
+      implicit none
+      integer hel
+      integer ipart
+      get_nhel = 99
+      return 
+      end
 
=== modified file 'Template/NLO/Cards/run_card.dat'
--- Template/NLO/Cards/run_card.dat	2021-04-05 16:07:15 +0000
+++ Template/NLO/Cards/run_card.dat	2021-11-12 09:03:34 +0000
@@ -129,7 +129,10 @@
  %(bwcutoff)s  = bwcutoff
 #***********************************************************************
 # Cuts on the jets. Jet clustering is performed by FastJet.            *
-#  - If gamma_is_j, photons are also clustered                            *
+#  - If gamma_is_j, photons are also clustered with jets.              *
+#    Otherwise, they will be treated as tagged particles and photon    * 
+#    isolation will be applied. Note that photons in the real emission *
+#    will always be clustered with QCD partons.                        *
 #  - 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        *
 
=== modified file 'Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f'
--- Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f	2018-10-26 09:55:25 +0000
+++ Template/NLO/FixedOrderAnalysis/analysis_HwU_general.f	2021-11-12 09:03:34 +0000
@@ -5,8 +5,12 @@
       integer nwgt
       character*(*) weights_info(*)
       integer i,l
-      character*6 cc(2)
-      data cc/'|T@NLO','|T@LO '/
+c      character*6 cc(2)
+c      data cc/'|T@NLO','|T@LO '/
+      character*13 cc(9)
+      data cc/ ' |T@NLO      ',' |T@LO       ',' |T@LO1      '
+     $        ,' |T@LO2      ',' |T@LO3      ',' |T@NLO1     '
+     $        ,' |T@NLO2     ',' |T@NLO3     ',' |T@NLO4     '/
 
 c     Also specific perturbative orders can be directly plotted, adding for examples
 c     the following further entries in the variable data
@@ -18,8 +22,8 @@
 c     
 c     See also line 376 in this file 
       call HwU_inithist(nwgt,weights_info)
-      do i=1,2
-      l=(i-1)*55
+      do i=1,9
+      l=(i-1)*59
 c transverse momenta
       call HwU_book(l+ 1,'total rate                 '//cc(i),1,0.5d0,1.5d0)
       call HwU_book(l+ 2,'1st charged lepton log pT  '//cc(i),25,-0.2d0,3.8d0)
@@ -76,6 +80,14 @@
 c HT
       call HwU_book(l+54,'log HT (partons)           '//cc(i),25,-0.2d0,3.8d0)
       call HwU_book(l+55,'log HT (reconstructed)     '//cc(i),25,-0.2d0,3.8d0)
+
+
+      call HwU_book(l+56,'1st isolated ph log pT     '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+57,'2nd isolated ph log pT     '//cc(i),25,-0.2d0,3.8d0)
+      call HwU_book(l+58,'3rd isolated ph log pT     '//cc(i),25,-0.2d0,3.8d0)
+
+      call HwU_book(l+59,'3 isolated ph log M        '//cc(i),25,-0.2d0,3.8d0)
+
       enddo
       return
       end
@@ -110,7 +122,7 @@
      $     ,ptmpmm,ptepve,ptmmvm,ptepemmpmm,ptepvemmvm,ptt,ptat,pttt
      $     ,etmiss,pth(2),pthh,ptv(3),ptjet(3),Mepem,Mmpmm,Mepve,Mmmvm
      $     ,Mepemmpmm,Mepvemmvm,Mtt,Mhh,Mj1j2,Mj1j3,Mj2j3,Mj1j2j3,Mvvv
-     $     ,HTparton,HTreco,p_reco(0:4,nexternal)
+     $     ,HTparton,HTreco,p_reco(0:4,nexternal),ptphiso(nexternal),Mphphph
       integer nQCD,jet(nexternal),njet,itop,iatop,iem,iep,imp,imm,ive
      $     ,ivm,iv1,iv2,iv3,ih1,ih2,il,ipdg_reco(nexternal)
       double precision getptv4,getptv4_2,getptv4_4,getinvm4_2
@@ -119,25 +131,47 @@
      $     ,getinvm4_4,l10
       integer orders_tag_plot
       common /corderstagplot/ orders_tag_plot
+
+
+c Photon isolation
+      integer nph,nem,nin,nphiso
+      double precision ptg,chi_gamma_iso,iso_getdrv40
+      double precision Etsum(0:nexternal)
+      real drlist(nexternal)
+      double precision pgamma(0:3,nexternal),pgammaiso(0:3,nexternal),pem(0:3,nexternal)
+      logical alliso,isolated
+c Sort array of results: ismode>0 for real, isway=0 for ascending order
+      integer ismode,isway,izero,isorted(nexternal)
+      parameter (ismode=1)
+      parameter (isway=0)
+      parameter (izero=0)
+      integer get_n_tagged_photons      
+
+
+      logical is_a_lp(nexternal),is_a_lm(nexternal),is_a_j(nexternal)
+     $     ,is_a_ph(nexternal)
+      REAL*8 pt,eta
+      external pt,eta,chi_gamma_iso,sortzv
+c      integer iph1,iph2,iph3
 c First, try to recombine photons with leptons      
-      if (.not.quarkphreco) then
-         write (*,*) 'quark-photon recombination is turned off. '/
-     $        /'Do need it'
-         stop
-      endif
-      if (.not. lepphreco) then
-         write (*,*) 'lepton-photon recombination is turned off. '
-     $        //'Do need it.'
-         stop
-      endif
+c      if (.not.quarkphreco) then
+c         write (*,*) 'quark-photon recombination is turned off. '/
+c     $        /'Do need it'
+c         stop
+c      endif
+c      if (.not. lepphreco) then
+c         write (*,*) 'lepton-photon recombination is turned off. '
+c     $        //'Do need it.'
+c         stop
+c      endif
       call recombine_momenta(rphreco, etaphreco, lepphreco, quarkphreco,
      $                       p, iPDG, p_reco, iPDG_reco)
 
 c Put all (light) QCD partons(+photon) in momentum array for jet clustering.
       nQCD=0
       do j=nincoming+1,nexternal
-         if (abs(ipdg_reco(j)).le.5 .or. ipdg_reco(j).eq.21 .or.
-     $        ipdg_reco(j).eq.22)then
+         if (abs(ipdg_reco(j)).le.5 .or. ipdg_reco(j).eq.21 
+     $      .or.  (ipdg_reco(j).eq.22.and.gamma_is_j)) then
             nQCD=nQCD+1
             do i=0,3
                pQCD(i,nQCD)=p_reco(i,j)
@@ -173,6 +207,161 @@
 c
 c******************************************************************************
 
+
+c PHOTON (ISOLATION) CUTS
+c
+c find the photons
+      do i=1,nexternal
+         if (istatus(i).eq.1 .and. ipdg(i).eq.22 .and. .not.gamma_is_j) then
+            is_a_ph(i)=.true.
+         else
+            is_a_ph(i)=.false.
+         endif
+      enddo
+      if (ptgmin.ne.0d0) then
+         nph=0
+         do j=nincoming+1,nexternal
+            if (is_a_ph(j)) then
+               nph=nph+1
+               do i=0,3
+                  pgamma(i,nph)=p(i,j)
+               enddo
+            endif
+         enddo
+         if(nph.eq.0)goto 444
+c         write(*,*) 'ERROR in cuts.f: photon isolation is not working'
+c     $           // ' for mixed QED-QCD corrections'
+c         stop 1
+
+         if(isoEM)then
+            nem=nph
+            do k=1,nem
+               do i=0,3
+                  pem(i,k)=pgamma(i,k)
+               enddo
+            enddo
+            do j=nincoming+1,nexternal
+               if (is_a_lp(j).or.is_a_lm(j)) then
+                  nem=nem+1
+                  do i=0,3
+                     pem(i,nem)=p(i,j)
+                  enddo
+               endif
+            enddo
+         endif
+
+c         alliso=.true.
+         nphiso=0
+
+         j=0
+c         do while(j.lt.nph.and.alliso)
+         do while(j.lt.nph)
+
+c Loop over all photons
+            j=j+1
+
+            ptg=pt(pgamma(0,j))
+            if(ptg.lt.ptgmin)then
+               cycle
+c               return
+            endif
+            if (etagamma.gt.0d0) then
+               if (abs(eta(pgamma(0,j))).gt.etagamma) then
+                  cycle
+c                  return
+               endif
+            endif
+
+c Isolate from hadronic energy
+            do i=1,nQCD
+               drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pQCD(0,i)))
+            enddo
+            call sortzv(drlist,isorted,nQCD,ismode,isway,izero)
+            Etsum(0)=0.d0
+            nin=0
+            do i=1,nQCD
+               if(dble(drlist(isorted(i))).le.R0gamma)then
+                  nin=nin+1
+                  Etsum(nin)=Etsum(nin-1)+pt(pQCD(0,isorted(i)))
+               endif
+            enddo
+            isolated=.True.
+            do i=1,nin
+c               alliso=alliso .and.
+c     $              Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
+c     $              R0gamma,xn,epsgamma,ptg)
+                if(Etsum(i).gt.chi_gamma_iso(dble(drlist(isorted(i))),
+     $              R0gamma,xn,epsgamma,ptg)) then
+                    isolated=isolated.and..False.
+                endif
+            enddo
+            if(.not.isolated)cycle
+
+c Isolate from EM energy
+            if(isoEM.and.nem.gt.1)then
+               do i=1,nem
+                  drlist(i)=sngl(iso_getdrv40(pgamma(0,j),pem(0,i)))
+               enddo
+               call sortzv(drlist,isorted,nem,ismode,isway,izero)
+c First of list must be the photon: check this, and drop it
+               if(isorted(1).ne.j.or.drlist(isorted(1)).gt.1.e-4)then
+                  write(*,*)'Error #1 in photon isolation'
+                  write(*,*)j,isorted(1),drlist(isorted(1))
+                  stop
+               endif
+               Etsum(0)=0.d0
+               nin=0
+               do i=2,nem
+                  if(dble(drlist(isorted(i))).le.R0gamma)then
+                     nin=nin+1
+                     Etsum(nin)=Etsum(nin-1)+pt(pem(0,isorted(i)))
+                  endif
+               enddo
+               isolated=.True.
+               do i=1,nin
+c                  alliso=alliso .and.
+c     $               Etsum(i).le.chi_gamma_iso(dble(drlist(isorted(i))),
+c     $               R0gamma,xn,epsgamma,ptg)
+                if(Etsum(i).gt.chi_gamma_iso(dble(drlist(isorted(i))),
+     $               R0gamma,xn,epsgamma,ptg)) then
+                    isolated=isolated.and..False.
+                endif
+            enddo
+            if(.not.isolated)cycle
+            endif
+c End of loop over photons
+
+         nphiso=nphiso+1
+
+         do i=0,3
+           pgammaiso(i,nphiso)=pgamma(i,j)
+         enddo
+
+
+         enddo
+         if(nphiso.lt.get_n_tagged_photons())then
+            print*,"mismatch with cuts.f"
+            stop
+         endif
+
+
+444     continue
+c End photon isolation
+       endif
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 c Look for the other physics objects
       itop=0
       iatop=0
@@ -189,7 +378,14 @@
       iv3=0
       ih1=0
       ih2=0
+c      iph1=0
+c      iph2=0
+c      iph3=0
+c          print*,"nell'analisi"
       do i=1,nexternal
+c          print*,"idpg di ",i,"=",ipdg(i)
+c          print*,"idpg_reco di ",i,"=",ipdg_reco(i)
+
          if (ipdg_reco(i).eq.6) then
             itop=i
          elseif(ipdg_reco(i).eq.-6) then
@@ -232,8 +428,27 @@
                   stop
                endif
             endif
+c         elseif(abs(ipdg_reco(i)).eq.22.and.i.gt.nincoming) then
+c            if (iph1.eq.0) then
+c               iph1=i
+c            else
+c               if (iph2.eq.0) then
+c                  iph2=i
+c               else
+c                  if (iph3.eq.0) then
+c                      iph3=i
+c                  else  
+c                     write (*,*) 'too many photonss'
+c                     stop
+c                  endif
+c               endif
+c            endif
          endif
+
       enddo
+c      print*,itop,iatop
+c      stop
+
       if (itop.ne.0) ptt=getptv4(p_reco(0,itop))
       if (iatop.ne.0) ptat=getptv4(p_reco(0,iatop))
       if (itop.ne.0 .and. iatop.ne.0) then
@@ -318,11 +533,21 @@
             pth(2)=tmp
          endif
       endif
+
+
       if (iv1.ne.0) ptv(1)=getptv4(p_reco(0,iv1))
       if (iv2.ne.0) ptv(2)=getptv4(p_reco(0,iv2))
       if (iv3.ne.0) ptv(3)=getptv4(p_reco(0,iv3))
       if (iv1.ne.0 .and. iv2.ne.0 .and. iv3.ne.0) then
          Mvvv=getinvm4_3(p_reco(0,iv1),p_reco(0,iv2),p_reco(0,iv3))
+      endif
+
+      do i=1,nphiso
+         ptphiso(i)=getptv4(pgammaiso(0,i))
+      enddo
+
+      if (iv1.ne.0 .and. iv2.ne.0 .and. iv3.ne.0) then
+         Mvvv=getinvm4_3(p_reco(0,iv1),p_reco(0,iv2),p_reco(0,iv3))
 c     order the vector bosons (if there are 3)
          do i=1,2
             do j=1,3-i
@@ -341,6 +566,31 @@
             ptv(2)=tmp
          endif
       endif
+
+
+
+
+      if (nphiso.eq.3) then
+         Mphphph=getinvm4_3(pgammaiso(0,1),pgammaiso(0,2),pgammaiso(0,3))
+c     order the isolated photons (if there are 3)
+         do i=1,2
+            do j=1,3-i
+               if (ptphiso(j).lt.ptphiso(j+1)) then
+                  tmp=ptphiso(j)
+                  ptphiso(j)=ptphiso(j+1)
+                  ptphiso(j+1)=tmp
+               endif
+            enddo
+         enddo
+      elseif (nphiso.eq.2) then
+c     order the isolated photons (if there are 2)
+         if (ptphiso(1).lt.ptphiso(2)) then
+            tmp=ptphiso(1)
+            ptphiso(1)=ptphiso(2)
+            ptphiso(2)=tmp
+         endif
+      endif
+      
       do i=1,njet
          ptjet(i)=getptv4(pjet(0,i))
       enddo
@@ -369,9 +619,18 @@
       enddo
       if (ive.ne.0 .or. ivm.ne.0) HTreco=HTreco+etmiss
          
-      do i=1,2
-         l=(i-1)*55
+      do i=1,9
+         l=(i-1)*59
          if (ibody.ne.3 .and.i.eq.2) cycle
+        if (i.eq. 3.and.orders_tag_plot.ne.204) cycle
+         if (i.eq. 4.and.orders_tag_plot.ne.402) cycle
+         if (i.eq. 5.and.orders_tag_plot.ne.600) cycle
+         if (i.eq. 6.and.orders_tag_plot.ne.206) cycle
+         if (i.eq. 7.and.orders_tag_plot.ne.404) cycle
+         if (i.eq. 8.and.orders_tag_plot.ne.602) cycle
+         if (i.eq. 9.and.orders_tag_plot.ne.800) cycle
+
+
 
 c         How to tag orders (QCD+QED*100)
 c
@@ -493,6 +752,18 @@
 c HT
          call HwU_fill(l+54,l10(HTparton),wgts)
          call HwU_fill(l+55,l10(HTreco),wgts)
+
+         if (nphiso.ge.1) call HwU_fill(l+56,l10(ptphiso(1)),wgts)
+         if (nphiso.ge.2) call HwU_fill(l+57,l10(ptphiso(2)),wgts)
+         if (nphiso.ge.3) call HwU_fill(l+58,l10(ptphiso(3)),wgts)
+
+         if (nphiso.ge.3) call HwU_fill(l+59,l10(Mphphph),wgts)
+
+
+
+
+
+
       enddo
 
  999  return      
 
=== added file 'Template/NLO/MCatNLO/Scripts/JetMatching.h'
--- Template/NLO/MCatNLO/Scripts/JetMatching.h	1970-01-01 00:00:00 +0000
+++ Template/NLO/MCatNLO/Scripts/JetMatching.h	2021-11-12 09:03:34 +0000
@@ -0,0 +1,2105 @@
+// JetMatching.h is a part of the PYTHIA event generator.
+// Copyright (C) 2021 Torbjorn Sjostrand.
+// PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
+// Please respect the MCnet Guidelines, see GUIDELINES for details.
+
+// Authors: Richard Corke (implementation of MLM matching as
+// in Alpgen for Alpgen input)
+// and Stephen Mrenna (implementation of MLM-style matching as
+// in Madgraph for Alpgen or Madgraph 5 input.)
+// and Simon de Visscher, Stefan Prestel (implementation of shower-kT
+// MLM-style matching and flavour treatment for Madgraph input)
+// and Stefan Prestel (FxFx NLO jet matching with aMC@NLO.)
+// This file provides the classes to perform MLM matching of
+// Alpgen or MadGraph 5 input.
+// Example usage is shown in main32.cc, and further details
+// can be found in the 'Jet Matching Style' manual page.
+
+#ifndef Pythia8_JetMatching_H
+#define Pythia8_JetMatching_H
+
+// Includes
+#include "Pythia8/Pythia.h"
+#include "Pythia8Plugins/GeneratorInput.h"
+
+namespace Pythia8 {
+
+//==========================================================================
+
+class HJSlowJet: public SlowJet {
+
+ public:
+
+  HJSlowJet(int powerIn, double Rin, double pTjetMinIn = 0.,
+            double etaMaxIn = 25., int selectIn = 1, int massSetIn = 2,
+            SlowJetHook* sjHookPtrIn = 0, bool useFJcoreIn = false,
+            bool useStandardRin = true) :
+  SlowJet(powerIn, Rin, pTjetMinIn, etaMaxIn, selectIn, massSetIn,
+          sjHookPtrIn, useFJcoreIn, useStandardRin) {}
+
+  // Note: The functions below have been made public to ease the generation
+  // of Python bindings.
+  //protected:
+
+  void findNext();
+
+};
+
+//--------------------------------------------------------------------------
+
+// Find next cluster pair to join.
+
+void HJSlowJet::findNext() {
+
+  // Find smallest of diB, dij.
+  if (clSize > 0) {
+    iMin =  0;
+    jMin = -1;
+    dMin = 1.0/TINY;
+    // Remove the possibility of choosing a beam clustering
+    for (int i = 1; i < clSize; ++i) {
+      for (int j = 0; j < i; ++j) {
+        if (dij[i*(i-1)/2 + j] < dMin) {
+          iMin = i;
+          jMin = j;
+          dMin = dij[i*(i-1)/2 + j];
+        }
+      }
+    }
+
+  // If no clusters left then instead default values.
+  } else {
+    iMin = -1;
+    jMin = -1;
+    dMin = 0.;
+  }
+
+}
+
+//==========================================================================
+
+// Declaration of main JetMatching class to perform MLM matching.
+// Note that it is defined with virtual inheritance, so that it can
+// be combined with other UserHooks classes, see e.g. main33.cc.
+
+class JetMatching : virtual public UserHooks {
+
+public:
+
+  // Constructor and destructor
+ JetMatching() : cellJet(NULL), slowJet(NULL), slowJetHard(NULL),
+    hjSlowJet(NULL) {}
+ ~JetMatching() {
+    if (cellJet) delete cellJet;
+    if (slowJet) delete slowJet;
+    if (slowJetHard) delete slowJetHard;
+    if (hjSlowJet) delete hjSlowJet;
+    // Print error statistics before exiting. Printing code
+    // basically copied from Info class.
+    // Header.
+    cout << "\n *-------  JetMatching Error and Warning Messages Statistics"
+         << "  -----------------------------------------------------* \n"
+         << " |                                                       "
+         << "                                                          | \n"
+         << " |  times   message                                      "
+         << "                                                          | \n"
+         << " |                                                       "
+         << "                                                          | \n";
+
+    // Loop over all messages
+    map<string, int>::iterator messageEntry = messages.begin();
+    if (messageEntry == messages.end())
+      cout << " |      0   no errors or warnings to report              "
+           << "                                                          | \n";
+    while (messageEntry != messages.end()) {
+      // Message printout.
+      string temp = messageEntry->first;
+      int len = temp.length();
+      temp.insert( len, max(0, 102 - len), ' ');
+      cout << " | " << setw(6) << messageEntry->second << "   "
+           << temp << " | \n";
+      ++messageEntry;
+    }
+
+    // Done.
+    cout << " |                                                       "
+         << "                                                          | \n"
+         << " *-------  End JetMatching Error and Warning Messages "
+         << "Statistics  -------------------------------------------------* "
+         << endl;
+  }
+
+  // Initialisation
+  virtual bool initAfterBeams() = 0;
+
+  // Process level vetos
+  bool canVetoProcessLevel() { return doMerge; }
+  bool doVetoProcessLevel(Event& process) {
+    eventProcessOrig = process;
+    return false;
+  }
+
+  // Parton level vetos (before beam remnants and resonance decays)
+  bool canVetoPartonLevelEarly() { return doMerge; }
+  bool doVetoPartonLevelEarly(const Event& event);
+
+  // Shower step vetoes (after the first emission, for Shower-kT scheme)
+  int  numberVetoStep() {return 1;}
+  bool canVetoStep() { return false; }
+  bool doVetoStep(int,  int, int, const Event& ) { return false; }
+
+  // Note: The functions below have been made public to ease the generation
+  // of Python bindings.
+  //protected:
+
+  // Different steps of the matching algorithm.
+  virtual void sortIncomingProcess(const Event &)=0;
+  virtual void jetAlgorithmInput(const Event &, int)=0;
+  virtual void runJetAlgorithm()=0;
+  virtual bool matchPartonsToJets(int)=0;
+  virtual int  matchPartonsToJetsLight()=0;
+  virtual int  matchPartonsToJetsHeavy()=0;
+
+  // Print a message the first few times. Insert in database.
+  void errorMsg(string messageIn) {
+    // Recover number of times message occured. Also inserts new string.
+    int times = messages[messageIn];
+    ++messages[messageIn];
+    // Print message the first few times.
+    if (times < TIMESTOPRINT) cout << " PYTHIA " << messageIn << endl;
+  }
+
+protected:
+
+  // Constants to be changed for debug printout or extra checks.
+  static const bool MATCHINGDEBUG, MATCHINGCHECK;
+
+  enum vetoStatus { NONE, LESS_JETS, MORE_JETS, HARD_JET,
+                    UNMATCHED_PARTON, INCLUSIVE_VETO };
+  enum partonTypes { ID_CHARM=4, ID_BOT=5, ID_TOP=6, ID_LEPMIN=11,
+    ID_LEPMAX=16, ID_GLUON=21, ID_PHOTON=22 };
+
+  // Master switch for merging
+  bool   doMerge;
+  // Switch for merging in the shower-kT scheme. Needed here because
+  // the scheme uses different UserHooks functionality.
+  bool   doShowerKt;
+
+  // Maximum and current number of jets
+  int    nJetMax, nJet;
+
+  // Jet algorithm parameters
+  int    jetAlgorithm;
+  double eTjetMin, coneRadius, etaJetMax, etaJetMaxAlgo;
+
+  // Internal jet algorithms
+  CellJet* cellJet;
+  SlowJet* slowJet;
+  SlowJet* slowJetHard;
+  HJSlowJet* hjSlowJet;
+
+  // SlowJet specific
+  int    slowJetPower;
+
+  // Event records to store original incoming process, final-state of the
+  // incoming process and what will be passed to the jet algorithm.
+  // Not completely necessary to store all steps, but makes tracking the
+  // steps of the algorithm a lot easier.
+  Event eventProcessOrig, eventProcess, workEventJet;
+
+  // Sort final-state of incoming process into light/heavy jets and 'other'
+  vector<int> typeIdx[3];
+  set<int>    typeSet[3];
+
+  // Momenta output of jet algorithm (to provide same output regardless of
+  // the selected jet algorithm)
+  vector<Vec4> jetMomenta;
+
+  // CellJet specific
+  int    nEta, nPhi;
+  double eTseed, eTthreshold;
+
+  // Merging procedure parameters
+  int    jetAllow, jetMatch, exclusiveMode;
+  double coneMatchLight, coneRadiusHeavy, coneMatchHeavy;
+  bool   exclusive;
+
+  // Store the minimum eT/pT of matched light jets
+  double eTpTlightMin;
+
+  // Map for all error messages.
+  map<string, int> messages;
+  // Number of times the same error message is repeated, unless overridden.
+  static const int TIMESTOPRINT = 1;
+
+};
+
+//==========================================================================
+
+// Declaration of main UserHooks class to perform Alpgen matching.
+
+class JetMatchingAlpgen : virtual public JetMatching {
+
+public:
+
+  // Constructor and destructor
+  JetMatchingAlpgen() { }
+  ~JetMatchingAlpgen() { }
+
+  // Initialisation
+  bool initAfterBeams();
+
+private:
+
+  // Different steps of the matching algorithm.
+  void sortIncomingProcess(const Event &);
+  void jetAlgorithmInput(const Event &, int);
+  void runJetAlgorithm();
+  bool matchPartonsToJets(int);
+  int  matchPartonsToJetsLight();
+  int  matchPartonsToJetsHeavy();
+
+  // Sorting utility
+  void sortTypeIdx(vector < int > &vecIn);
+
+  // Constants
+  static const double GHOSTENERGY, ZEROTHRESHOLD;
+
+};
+
+//==========================================================================
+
+// Declaration of main UserHooks class to perform Madgraph matching.
+
+class JetMatchingMadgraph : virtual public JetMatching {
+
+public:
+
+  // Constructor and destructor
+  JetMatchingMadgraph() : slowJetDJR(NULL) { }
+  ~JetMatchingMadgraph() { if (slowJetDJR) delete slowJetDJR; }
+
+  // Initialisation
+  bool initAfterBeams();
+
+  // Process level vetos
+  bool canVetoProcessLevel() { return doMerge; }
+  bool doVetoProcessLevel(Event& process);
+
+  // Shower step vetoes (after the first emission, for Shower-kT scheme)
+  int  numberVetoStep() {return 1;}
+  bool canVetoStep() { return doShowerKt; }
+  bool doVetoStep(int,  int, int, const Event& );
+
+  // Jet algorithm to access the jet separations in the cleaned event
+  // after showering.
+  SlowJet* slowJetDJR;
+  // Functions to return the jet clustering scales and number of ME partons.
+  // These are useful to investigate the matching systematics.
+  vector<double> getDJR() { return DJR;}
+  pair<int,int> nMEpartons() { return nMEpartonsSave;}
+
+  // For systematic variations of the jet matching parameters, it is helpful
+  // to decouple the jet matching veto from the internal book-keeping. The
+  // veto can then be applied in hindsight by an expert user. The functions
+  // below return all the necessary information to do this.
+  Event getWorkEventJet() { return workEventJetSave; }
+  Event getProcessSubset() { return processSubsetSave; }
+  bool  getExclusive() { return exclusive; }
+  double getPTfirst() { return pTfirstSave; }
+
+  // Note: The functions below have been made public to ease the generation
+  // of Python bindings.
+  //protected:
+
+  // Different steps of the matching algorithm.
+  void sortIncomingProcess(const Event &);
+  void jetAlgorithmInput(const Event &, int);
+  void runJetAlgorithm();
+  bool matchPartonsToJets(int);
+  int  matchPartonsToJetsLight();
+  int  matchPartonsToJetsHeavy();
+  int  matchPartonsToJetsOther();
+  bool doShowerKtVeto(double pTfirst);
+
+  // Functions to clear and set the jet clustering scales.
+  void clearDJR() { DJR.resize(0);}
+  void setDJR( const Event& event);
+  // Functions to clear and set the jet clustering scales.
+  void clear_nMEpartons() { nMEpartonsSave.first = nMEpartonsSave.second =-1;}
+  void set_nMEpartons( const int nOrig, const int nMatch) {
+    clear_nMEpartons();
+    nMEpartonsSave.first  = nOrig;
+    nMEpartonsSave.second = nMatch;
+  };
+
+  // Function to get the current number of partons in the Born state, as
+  // read from LHE.
+  int npNLO();
+
+private:
+
+  // Stored values of all inputs necessary to perform the jet matching, as
+  // needed when the veto is applied externally.
+  Event processSubsetSave;
+  Event workEventJetSave;
+  double pTfirstSave;
+
+  // Save if code should apply the veto, or simply store the things necessary
+  // to perform the veto externally.
+  bool performVeto;
+
+  // Variables.
+  vector<int> origTypeIdx[3];
+  int    nQmatch;
+  double qCut, qCutSq, clFact;
+  bool   doFxFx;
+  int    nPartonsNow;
+  double qCutME, qCutMESq;
+
+  // Vector to store the jet clustering scales.
+  vector<double> DJR;
+  // Pair of integers giving the number of ME partons read from LHEF and used
+  // in the matching (can be different if some partons should not be matched)
+  pair<int,int> nMEpartonsSave;
+
+};
+
+//==========================================================================
+
+// Main implementation of JetMatching class.
+// This may be split out to a separate C++ file if desired,
+// but currently included here for ease of use.
+
+//--------------------------------------------------------------------------
+
+// Constants to be changed for debug printout or extra checks.
+const bool JetMatching::MATCHINGDEBUG = false;
+const bool JetMatching::MATCHINGCHECK = false;
+
+//--------------------------------------------------------------------------
+
+// Early parton level veto (before beam remnants and resonance showers)
+
+inline bool JetMatching::doVetoPartonLevelEarly(const Event& event) {
+
+  // 1) Sort the original incoming process. After this step is performed,
+  //    the following assignments have been made:
+  //    eventProcessOrig - the original incoming process
+  //    eventProcess     - the final-state of the incoming process with
+  //                       resonance decays removed (and resonances
+  //                       themselves now with positive status code)
+  //    typeIdx[0/1/2]   - Indices into 'eventProcess' of
+  //                       light jets/heavy jets/other
+  //    typeSet[0/1/2]   - Indices into 'event' of light jets/heavy jets/other
+  //    workEvent        - partons from the hardest subsystem + ISR + FSR only
+  sortIncomingProcess(event);
+
+  // For the shower-kT scheme, do not perform any veto here, as any vetoing
+  // will already have taken place in doVetoStep.
+  if ( doShowerKt ) return false;
+
+  // Debug printout.
+  if (MATCHINGDEBUG) {
+    // Begin
+    cout << endl << "-------- Begin Madgraph Debug --------" << endl;
+    // Original incoming process
+    cout << endl << "Original incoming process:";
+    eventProcessOrig.list();
+    // Final-state of original incoming process
+    cout << endl << "Final-state incoming process:";
+    eventProcess.list();
+    // List categories of sorted particles
+    for (size_t i = 0; i < typeIdx[0].size(); i++)
+      cout << ((i == 0) ? "Light jets: " : ", ")   << setw(3) << typeIdx[0][i];
+    if( typeIdx[0].size()== 0 )
+      cout << "Light jets: None";
+
+    for (size_t i = 0; i < typeIdx[1].size(); i++)
+      cout << ((i == 0) ? "\nHeavy jets: " : ", ") << setw(3) << typeIdx[1][i];
+    for (size_t i = 0; i < typeIdx[2].size(); i++)
+      cout << ((i == 0) ? "\nOther:      " : ", ") << setw(3) << typeIdx[2][i];
+    // Full event at this stage
+    cout << endl << endl << "Event:";
+    event.list();
+    // Work event (partons from hardest subsystem + ISR + FSR)
+    cout << endl << "Work event:";
+    workEvent.list();
+  }
+
+  // 2) Light/heavy jets: iType = 0 (light jets), 1 (heavy jets)
+  int iTypeEnd = (typeIdx[2].empty()) ? 2 : 3;
+  for (int iType = 0; iType < iTypeEnd; iType++) {
+
+    // 2a) Find particles which will be passed from the jet algorithm.
+    //     Input from 'workEvent' and output in 'workEventJet'.
+    jetAlgorithmInput(event, iType);
+
+    // Debug printout.
+    if (MATCHINGDEBUG) {
+      // Jet algorithm event
+      cout << endl << "Jet algorithm event (iType = " << iType << "):";
+      workEventJet.list();
+    }
+
+    // 2b) Run jet algorithm on 'workEventJet'.
+    //     Output is stored in jetMomenta.
+    runJetAlgorithm();
+
+    // 2c) Match partons to jets and decide if veto is necessary
+    if (matchPartonsToJets(iType) == true) {
+      // Debug printout.
+      if (MATCHINGDEBUG) {
+        cout << endl << "Event vetoed" << endl
+             << "----------  End MLM Debug  ----------" << endl;
+      }
+      return true;
+    }
+  }
+
+  // Debug printout.
+  if (MATCHINGDEBUG) {
+    cout << endl << "Event accepted" << endl
+         << "----------  End MLM Debug  ----------" << endl;
+  }
+
+  // If we reached here, then no veto
+  return false;
+
+}
+
+//==========================================================================
+
+// Main implementation of Alpgen UserHooks class.
+// This may be split out to a separate C++ file if desired,
+// but currently included here for ease of use.
+
+//--------------------------------------------------------------------------
+
+// Constants: could be changed here if desired, but normally should not.
+// These are of technical nature, as described for each.
+
+// The energy of ghost particles. For technical reasons, this cannot be
+// set arbitrarily low, see 'Particle::TINY' in 'Event.cc' for details.
+const double JetMatchingAlpgen::GHOSTENERGY   = 1e-15;
+
+// A zero threshold value for double comparisons.
+const double JetMatchingAlpgen::ZEROTHRESHOLD = 1e-10;
+
+//--------------------------------------------------------------------------
+
+// Function to sort typeIdx vectors into descending eT/pT order.
+// Uses a selection sort, as number of partons generally small
+// and so efficiency not a worry.
+
+inline void JetMatchingAlpgen::sortTypeIdx(vector < int > &vecIn) {
+  for (size_t i = 0; i < vecIn.size(); i++) {
+    size_t jMax = i;
+    double vMax = (jetAlgorithm == 1) ?
+      eventProcess[vecIn[i]].eT() :
+      eventProcess[vecIn[i]].pT();
+    for (size_t j = i + 1; j < vecIn.size(); j++) {
+      double vNow = (jetAlgorithm == 1)
+        ? eventProcess[vecIn[j]].eT() : eventProcess[vecIn[j]].pT();
+      if (vNow > vMax) {
+        vMax = vNow;
+        jMax = j;
+      }
+    }
+    if (jMax != i) swap(vecIn[i], vecIn[jMax]);
+  }
+}
+
+//--------------------------------------------------------------------------
+
+// Initialisation routine automatically called from Pythia::init().
+// Setup all parts needed for the merging.
+
+inline bool JetMatchingAlpgen::initAfterBeams() {
+
+  // Read in parameters
+  doMerge         = settingsPtr->flag("JetMatching:merge");
+  jetAlgorithm    = settingsPtr->mode("JetMatching:jetAlgorithm");
+  nJet            = settingsPtr->mode("JetMatching:nJet");
+  nJetMax         = settingsPtr->mode("JetMatching:nJetMax");
+  eTjetMin        = settingsPtr->parm("JetMatching:eTjetMin");
+  coneRadius      = settingsPtr->parm("JetMatching:coneRadius");
+  etaJetMax       = settingsPtr->parm("JetMatching:etaJetMax");
+  doShowerKt      = settingsPtr->flag("JetMatching:doShowerKt");
+
+  // Use etaJetMax + coneRadius in input to jet algorithms
+  etaJetMaxAlgo   = etaJetMax + coneRadius;
+
+  // CellJet specific
+  nEta            = settingsPtr->mode("JetMatching:nEta");
+  nPhi            = settingsPtr->mode("JetMatching:nPhi");
+  eTseed          = settingsPtr->parm("JetMatching:eTseed");
+  eTthreshold     = settingsPtr->parm("JetMatching:eTthreshold");
+
+  // SlowJet specific
+  slowJetPower    = settingsPtr->mode("JetMatching:slowJetPower");
+  coneMatchLight  = settingsPtr->parm("JetMatching:coneMatchLight");
+  coneRadiusHeavy = settingsPtr->parm("JetMatching:coneRadiusHeavy");
+  if (coneRadiusHeavy < 0.) coneRadiusHeavy = coneRadius;
+  coneMatchHeavy  = settingsPtr->parm("JetMatching:coneMatchHeavy");
+
+  // Matching procedure
+  jetAllow        = settingsPtr->mode("JetMatching:jetAllow");
+  jetMatch        = settingsPtr->mode("JetMatching:jetMatch");
+  exclusiveMode   = settingsPtr->mode("JetMatching:exclusive");
+
+  // If not merging, then done
+  if (!doMerge) return true;
+
+  // Exclusive mode; if set to 2, then set based on nJet/nJetMax
+  if (exclusiveMode == 2) {
+
+    // No nJet or nJetMax, so default to exclusive mode
+    if (nJet < 0 || nJetMax < 0) {
+      errorMsg("Warning in JetMatchingAlpgen:init: "
+          "missing jet multiplicity information; running in exclusive mode");
+      exclusive = true;
+
+    // Inclusive if nJet == nJetMax, exclusive otherwise
+    } else {
+      exclusive = (nJet == nJetMax) ? false : true;
+    }
+
+  // Otherwise, just set as given
+  } else {
+    exclusive = (exclusiveMode == 0) ? false : true;
+  }
+
+  // Initialise chosen jet algorithm. CellJet.
+  if (jetAlgorithm == 1) {
+
+    // Extra options for CellJet. nSel = 1 means that all final-state
+    // particles are taken and we retain control of what to select.
+    // smear/resolution/upperCut are not used and are set to default values.
+    int    nSel = 2, smear = 0;
+    double resolution = 0.5, upperCut = 2.;
+    cellJet = new CellJet(etaJetMaxAlgo, nEta, nPhi, nSel,
+                          smear, resolution, upperCut, eTthreshold);
+
+  // SlowJet
+  } else if (jetAlgorithm == 2) {
+    slowJet = new SlowJet(slowJetPower, coneRadius, eTjetMin, etaJetMaxAlgo);
+  }
+
+  // Check the jetMatch parameter; option 2 only works with SlowJet
+  if (jetAlgorithm == 1 && jetMatch == 2) {
+    errorMsg("Warning in JetMatchingAlpgen:init: "
+        "jetMatch = 2 only valid with SlowJet algorithm. "
+        "Reverting to jetMatch = 1");
+    jetMatch = 1;
+  }
+
+  // Setup local event records
+  eventProcessOrig.init("(eventProcessOrig)", particleDataPtr);
+  eventProcess.init("(eventProcess)", particleDataPtr);
+  workEventJet.init("(workEventJet)", particleDataPtr);
+
+  // Print information
+  string jetStr  = (jetAlgorithm ==  1) ? "CellJet" :
+                   (slowJetPower == -1) ? "anti-kT" :
+                   (slowJetPower ==  0) ? "C/A"     :
+                   (slowJetPower ==  1) ? "kT"      : "unknown";
+  string modeStr = (exclusive)         ? "exclusive" : "inclusive";
+  stringstream nJetStr, nJetMaxStr;
+  if (nJet >= 0)    nJetStr    << nJet;    else nJetStr    << "unknown";
+  if (nJetMax >= 0) nJetMaxStr << nJetMax; else nJetMaxStr << "unknown";
+  cout << endl
+       << " *-------  MLM matching parameters  -------*" << endl
+       << " |  nJet                |  " << setw(14)
+       << nJetStr.str() << "  |" << endl
+       << " |  nJetMax             |  " << setw(14)
+       << nJetMaxStr.str() << "  |" << endl
+       << " |  Jet algorithm       |  " << setw(14)
+       << jetStr << "  |" << endl
+       << " |  eTjetMin            |  " << setw(14)
+       << eTjetMin << "  |" << endl
+       << " |  coneRadius          |  " << setw(14)
+       << coneRadius << "  |" << endl
+       << " |  etaJetMax           |  " << setw(14)
+       << etaJetMax << "  |" << endl
+       << " |  jetAllow            |  " << setw(14)
+       << jetAllow << "  |" << endl
+       << " |  jetMatch            |  " << setw(14)
+       << jetMatch << "  |" << endl
+       << " |  coneMatchLight      |  " << setw(14)
+       << coneMatchLight << "  |" << endl
+       << " |  coneRadiusHeavy     |  " << setw(14)
+       << coneRadiusHeavy << "  |" << endl
+       << " |  coneMatchHeavy      |  " << setw(14)
+       << coneMatchHeavy << "  |" << endl
+       << " |  Mode                |  " << setw(14)
+       << modeStr << "  |" << endl
+       << " *-----------------------------------------*" << endl;
+
+  return true;
+}
+
+//--------------------------------------------------------------------------
+
+// Step (1): sort the incoming particles
+
+inline void JetMatchingAlpgen::sortIncomingProcess(const Event &event) {
+
+  // Remove resonance decays from original process and keep only final
+  // state. Resonances will have positive status code after this step.
+  omitResonanceDecays(eventProcessOrig, true);
+  eventProcess = workEvent;
+
+  // Sort original process final state into light/heavy jets and 'other'.
+  // Criteria:
+  //   1 <= ID <= 5 and massless, or ID == 21 --> light jet (typeIdx[0])
+  //   4 <= ID <= 6 and massive               --> heavy jet (typeIdx[1])
+  //   All else                               --> other     (typeIdx[2])
+  // Note that 'typeIdx' stores indices into 'eventProcess' (after resonance
+  // decays are omitted), while 'typeSet' stores indices into the original
+  // process record, 'eventProcessOrig', but these indices are also valid
+  // in 'event'.
+  for (int i = 0; i < 3; i++) {
+    typeIdx[i].clear();
+    typeSet[i].clear();
+  }
+  for (int i = 0; i < eventProcess.size(); i++) {
+    // Ignore nonfinal and default to 'other'
+    if (!eventProcess[i].isFinal()) continue;
+    int idx = 2;
+
+    // Light jets
+    if (eventProcess[i].id() == ID_GLUON
+      || (eventProcess[i].idAbs() <= ID_BOT
+      && abs(eventProcess[i].m()) < ZEROTHRESHOLD)) idx = 0;
+
+    // Heavy jets
+    else if (eventProcess[i].idAbs() >= ID_CHARM
+      && eventProcess[i].idAbs() <= ID_TOP) idx = 1;
+
+    // Store
+    typeIdx[idx].push_back(i);
+    typeSet[idx].insert(eventProcess[i].daughter1());
+  }
+
+  // Extract partons from hardest subsystem + ISR + FSR only into
+  // workEvent. Note no resonance showers or MPIs.
+  subEvent(event);
+}
+
+//--------------------------------------------------------------------------
+
+// Step (2a): pick which particles to pass to the jet algorithm
+
+inline void JetMatchingAlpgen::jetAlgorithmInput(const Event &event,
+  int iType) {
+
+  // Take input from 'workEvent' and put output in 'workEventJet'
+  workEventJet = workEvent;
+
+  // Loop over particles and decide what to pass to the jet algorithm
+  for (int i = 0; i < workEventJet.size(); ++i) {
+    if (!workEventJet[i].isFinal()) continue;
+
+    // jetAllow option to disallow certain particle types
+    if (jetAllow == 1) {
+
+      // Original AG+Py6 algorithm explicitly excludes tops,
+      // leptons and photons.
+      int id = workEventJet[i].idAbs();
+      if ( (id >= ID_LEPMIN && id <= ID_LEPMAX) || id == ID_TOP
+        || id == ID_PHOTON) {
+        workEventJet[i].statusNeg();
+        continue;
+      }
+    }
+
+    // Get the index of this particle in original event
+    int idx = workEventJet[i].daughter1();
+
+    // Start with particle idx, and afterwards track mothers
+    while (true) {
+
+      // Light jets
+      if (iType == 0) {
+
+        // Do not include if originates from heavy jet or 'other'
+        if (typeSet[1].find(idx) != typeSet[1].end() ||
+            typeSet[2].find(idx) != typeSet[2].end()) {
+          workEventJet[i].statusNeg();
+          break;
+        }
+
+        // Made it to start of event record so done
+        if (idx == 0) break;
+        // Otherwise next mother and continue
+        idx = event[idx].mother1();
+
+      // Heavy jets
+      } else if (iType == 1) {
+
+        // Only include if originates from heavy jet
+        if (typeSet[1].find(idx) != typeSet[1].end()) break;
+
+        // Made it to start of event record with no heavy jet mother,
+        // so DO NOT include particle
+        if (idx == 0) {
+          workEventJet[i].statusNeg();
+          break;
+        }
+
+        // Otherwise next mother and continue
+        idx = event[idx].mother1();
+
+      // Other jets
+      } else if (iType == 2) {
+
+        // Only include if originates from other jet
+        if (typeSet[2].find(idx) != typeSet[2].end()) break;
+
+        // Made it to start of event record with no heavy jet mother,
+        // so DO NOT include particle
+        if (idx == 0) {
+          workEventJet[i].statusNeg();
+          break;
+        }
+
+        // Otherwise next mother and continue
+        idx = event[idx].mother1();
+
+      } // if (iType)
+    } // while (true)
+  } // for (i)
+
+  // For jetMatch = 2, insert ghost particles corresponding to
+  // each hard parton in the original process
+  if (jetMatch == 2) {
+    for (int i = 0; i < int(typeIdx[iType].size()); i++) {
+      // Get y/phi of the parton
+      Vec4   pIn = eventProcess[typeIdx[iType][i]].p();
+      double y   = pIn.rap();
+      double phi = pIn.phi();
+
+      // Create a ghost particle and add to the workEventJet
+      double e   = GHOSTENERGY;
+      double e2y = exp(2. * y);
+      double pz  = e * (e2y - 1.) / (e2y + 1.);
+      double pt  = sqrt(e*e - pz*pz);
+      double px  = pt * cos(phi);
+      double py  = pt * sin(phi);
+      workEventJet.append( ID_GLUON, 99, 0, 0, 0, 0, 0, 0, px, py, pz, e);
+
+      // Extra check on reconstructed y/phi values. If many warnings
+      // of this type, GHOSTENERGY may be set too low.
+      if (MATCHINGCHECK) {
+      int lastIdx = workEventJet.size() - 1;
+      if (abs(y   - workEventJet[lastIdx].y())   > ZEROTHRESHOLD ||
+          abs(phi - workEventJet[lastIdx].phi()) > ZEROTHRESHOLD)
+        errorMsg("Warning in JetMatchingAlpgen:jetAlgorithmInput: "
+            "ghost particle y/phi mismatch");
+      }
+
+    } // for (i)
+  } // if (jetMatch == 2)
+}
+
+//--------------------------------------------------------------------------
+
+// Step (2b): run jet algorithm and provide common output
+
+inline void JetMatchingAlpgen::runJetAlgorithm() {
+
+  // Run the jet clustering algorithm
+  if (jetAlgorithm == 1)
+    cellJet->analyze(workEventJet, eTjetMin, coneRadius, eTseed);
+  else
+    slowJet->analyze(workEventJet);
+
+  // Extract four-momenta of jets with |eta| < etaJetMax and
+  // put into jetMomenta. Note that this is done backwards as
+  // jets are removed with SlowJet.
+  jetMomenta.clear();
+  int iJet = (jetAlgorithm == 1) ? cellJet->size() - 1:
+                                   slowJet->sizeJet() - 1;
+  for (int i = iJet; i > -1; i--) {
+    Vec4 jetMom = (jetAlgorithm == 1) ? cellJet->pMassive(i) :
+                                        slowJet->p(i);
+    double eta = jetMom.eta();
+
+    if (abs(eta) > etaJetMax) {
+      if (jetAlgorithm == 2) slowJet->removeJet(i);
+      continue;
+    }
+    jetMomenta.push_back(jetMom);
+  }
+
+  // Reverse jetMomenta to restore eT/pT ordering
+  reverse(jetMomenta.begin(), jetMomenta.end());
+}
+
+//--------------------------------------------------------------------------
+
+// Step (2c): veto decision (returning true vetoes the event)
+
+inline bool JetMatchingAlpgen::matchPartonsToJets(int iType) {
+
+  // Use two different routines for light/heavy jets as
+  // different veto conditions and for clarity
+  if (iType == 0) return (matchPartonsToJetsLight() > 0);
+  else if (iType == 1) return (matchPartonsToJetsHeavy() > 0);
+  else if (iType == 2) return false;
+  return true;
+}
+
+//--------------------------------------------------------------------------
+
+// Step(2c): light jets
+// Return codes are given indicating the reason for a veto.
+// Although not currently used, they are a useful debugging tool:
+//   0 = no veto
+//   1 = veto as number of jets less than number of partons
+//   2 = veto as exclusive mode and number of jets greater than
+//       number of partons
+//   3 = veto as inclusive mode and there would be an extra jet
+//       that is harder than any matched soft jet
+//   4 = veto as there is a parton which does not match a jet
+
+inline int JetMatchingAlpgen::matchPartonsToJetsLight() {
+
+  // Always veto if number of jets is less than original number of jets
+  if (jetMomenta.size() < typeIdx[0].size()) return LESS_JETS;
+  // Veto if in exclusive mode and number of jets bigger than original
+  if (exclusive && jetMomenta.size() > typeIdx[0].size()) return MORE_JETS;
+
+  // Sort partons by eT/pT
+  sortTypeIdx(typeIdx[0]);
+
+  // Number of hard partons
+  int nParton = typeIdx[0].size();
+
+  // Keep track of which jets have been assigned a hard parton
+  vector < bool > jetAssigned;
+  jetAssigned.assign(jetMomenta.size(), false);
+
+  // Jet matching procedure: (1) deltaR between partons and jets
+  if (jetMatch == 1) {
+
+    // Loop over light hard partons and get 4-momentum
+    for (int i = 0; i < nParton; i++) {
+      Vec4 p1 = eventProcess[typeIdx[0][i]].p();
+
+      // Track which jet has the minimal dR measure with this parton
+      int    jMin  = -1;
+      double dRmin = 0.;
+
+      // Loop over all jets (skipping those already assigned).
+      for (int j = 0; j < int(jetMomenta.size()); j++) {
+        if (jetAssigned[j]) continue;
+
+        // DeltaR between parton/jet and store if minimum
+        double dR = (jetAlgorithm == 1)
+          ? REtaPhi(p1, jetMomenta[j]) : RRapPhi(p1, jetMomenta[j]);
+        if (jMin < 0 || dR < dRmin) {
+          dRmin = dR;
+          jMin  = j;
+        }
+      } // for (j)
+
+      // Check for jet-parton match
+      if (jMin >= 0 && dRmin < coneRadius * coneMatchLight) {
+
+        // If the matched jet is not one of the nParton hardest jets,
+        // the extra left over jet would be harder than some of the
+        // matched jets. This is disallowed, so veto.
+        if (jMin >= nParton) return HARD_JET;
+
+        // Mark jet as assigned.
+        jetAssigned[jMin] = true;
+
+      // If no match, then event will be vetoed in all cases
+      } else return UNMATCHED_PARTON;
+
+    } // for (i)
+
+  // Jet matching procedure: (2) ghost particles in SlowJet
+  } else {
+
+    // Loop over added 'ghost' particles and find if assigned to a jet
+    for (int i = workEventJet.size() - nParton;
+        i < workEventJet.size(); i++) {
+      int jMin = slowJet->jetAssignment(i);
+
+      // Veto if:
+      //  1) not one of nParton hardest jets
+      //  2) not assigned to a jet
+      //  3) jet has already been assigned
+      if (jMin >= nParton)               return HARD_JET;
+      if (jMin < 0 || jetAssigned[jMin]) return UNMATCHED_PARTON;
+
+      // Mark jet as assigned
+      jetAssigned[jMin] = true;
+
+    } // for (i)
+  } // if (jetMatch)
+
+  // Minimal eT/pT (CellJet/SlowJet) of matched light jets. Needed
+  // later for heavy jet vetos in inclusive mode.
+  if (nParton > 0)
+    eTpTlightMin = (jetAlgorithm == 1) ? jetMomenta[nParton - 1].eT()
+                                       : jetMomenta[nParton - 1].pT();
+  else
+    eTpTlightMin = -1.;
+
+  // No veto
+  return NONE;
+}
+
+//--------------------------------------------------------------------------
+
+// Step(2c): heavy jets
+// Return codes are given indicating the reason for a veto.
+// Although not currently used, they are a useful debugging tool:
+//   0 = no veto as there are no extra jets present
+//   1 = veto as in exclusive mode and extra jets present
+//   2 = veto as in inclusive mode and extra jets were harder
+//       than any matched light jet
+
+inline int JetMatchingAlpgen::matchPartonsToJetsHeavy() {
+
+  // If there are no extra jets, then accept
+  if (jetMomenta.empty()) return NONE;
+
+  // Number of hard partons
+  int nParton = typeIdx[1].size();
+
+  // Remove jets that are close to heavy quarks
+  set < int > removeJets;
+
+  // Jet matching procedure: (1) deltaR between partons and jets
+  if (jetMatch == 1) {
+
+    // Loop over heavy hard partons and get 4-momentum
+    for (int i = 0; i < nParton; i++) {
+      Vec4 p1 = eventProcess[typeIdx[1][i]].p();
+
+      // Loop over all jets, find dR and mark for removal if match
+      for (int j = 0; j < int(jetMomenta.size()); j++) {
+        double dR = (jetAlgorithm == 1) ?
+            REtaPhi(p1, jetMomenta[j]) : RRapPhi(p1, jetMomenta[j]);
+        if (dR < coneRadiusHeavy * coneMatchHeavy)
+          removeJets.insert(j);
+
+      } // for (j)
+    } // for (i)
+
+  // Jet matching procedure: (2) ghost particles in SlowJet
+  } else {
+
+    // Loop over added 'ghost' particles and if assigned to a jet
+    // then mark this jet for removal
+    for (int i = workEventJet.size() - nParton;
+        i < workEventJet.size(); i++) {
+      int jMin = slowJet->jetAssignment(i);
+      if (jMin >= 0) removeJets.insert(jMin);
+    }
+
+  }
+
+  // Remove jets (backwards order to not disturb indices)
+  for (set < int >::reverse_iterator it  = removeJets.rbegin();
+                                     it != removeJets.rend(); it++)
+    jetMomenta.erase(jetMomenta.begin() + *it);
+
+  // Handle case if there are still extra jets
+  if (!jetMomenta.empty()) {
+
+    // Exclusive mode, so immediate veto
+    if (exclusive) return MORE_JETS;
+
+    // Inclusive mode; extra jets must be softer than any matched light jet
+    else if (eTpTlightMin >= 0.)
+      for (size_t j = 0; j < jetMomenta.size(); j++) {
+        // CellJet uses eT, SlowJet uses pT
+        if ( (jetAlgorithm == 1 && jetMomenta[j].eT() > eTpTlightMin) ||
+             (jetAlgorithm == 2 && jetMomenta[j].pT() > eTpTlightMin) )
+          return HARD_JET;
+      }
+
+  } // if (!jetMomenta.empty())
+
+  // No extra jets were present so no veto
+  return NONE;
+}
+
+//==========================================================================
+
+// Main implementation of Madgraph UserHooks class.
+// This may be split out to a separate C++ file if desired,
+// but currently included here for ease of use.
+
+//--------------------------------------------------------------------------
+
+// Initialisation routine automatically called from Pythia::init().
+// Setup all parts needed for the merging.
+
+inline bool JetMatchingMadgraph::initAfterBeams() {
+
+  // Initialise values for stored jet matching veto inputs.
+  pTfirstSave = -1.;
+  processSubsetSave.init("(eventProcess)", particleDataPtr);
+  workEventJetSave.init("(workEventJet)", particleDataPtr);
+
+  // Read in Madgraph specific configuration variables
+  bool setMad    = settingsPtr->flag("JetMatching:setMad");
+
+  // If Madgraph parameters are present, then parse in MadgraphPar object
+  MadgraphPar par;
+  string parStr = infoPtr->header("MGRunCard");
+  if (!parStr.empty()) {
+    par.parse(parStr);
+    par.printParams();
+  }
+
+  // Set Madgraph merging parameters from the file if requested
+  if (setMad) {
+    if ( par.haveParam("xqcut")    && par.haveParam("maxjetflavor")
+      && par.haveParam("alpsfact") && par.haveParam("ickkw") ) {
+      settingsPtr->flag("JetMatching:merge", par.getParam("ickkw"));
+      settingsPtr->parm("JetMatching:qCut", par.getParam("xqcut"));
+      settingsPtr->mode("JetMatching:nQmatch",
+        par.getParamAsInt("maxjetflavor"));
+      settingsPtr->parm("JetMatching:clFact",
+        clFact = par.getParam("alpsfact"));
+      if (par.getParamAsInt("ickkw") == 0)
+        errorMsg("Error in JetMatchingMadgraph:init: "
+          "Madgraph file parameters are not set for merging");
+
+    // Warn if setMad requested, but one or more parameters not present
+    } else {
+       errorMsg("Warning in JetMatchingMadgraph:init: "
+          "Madgraph merging parameters not found");
+       if (!par.haveParam("xqcut")) errorMsg("Warning in "
+          "JetMatchingMadgraph:init: No xqcut");
+       if (!par.haveParam("ickkw")) errorMsg("Warning in "
+          "JetMatchingMadgraph:init: No ickkw");
+       if (!par.haveParam("maxjetflavor")) errorMsg("Warning in "
+          "JetMatchingMadgraph:init: No maxjetflavor");
+       if (!par.haveParam("alpsfact")) errorMsg("Warning in "
+          "JetMatchingMadgraph:init: No alpsfact");
+    }
+  }
+
+  // Read in FxFx matching parameters
+  doFxFx       = settingsPtr->flag("JetMatching:doFxFx");
+  nPartonsNow  = settingsPtr->mode("JetMatching:nPartonsNow");
+  qCutME       = settingsPtr->parm("JetMatching:qCutME");
+  qCutMESq     = pow(qCutME,2);
+
+  // Read in Madgraph merging parameters
+  doMerge      = settingsPtr->flag("JetMatching:merge");
+  doShowerKt   = settingsPtr->flag("JetMatching:doShowerKt");
+  qCut         = settingsPtr->parm("JetMatching:qCut");
+  nQmatch      = settingsPtr->mode("JetMatching:nQmatch");
+  clFact       = settingsPtr->parm("JetMatching:clFact");
+
+  // Read in jet algorithm parameters
+  jetAlgorithm   = settingsPtr->mode("JetMatching:jetAlgorithm");
+  nJetMax        = settingsPtr->mode("JetMatching:nJetMax");
+  eTjetMin       = settingsPtr->parm("JetMatching:eTjetMin");
+  coneRadius     = settingsPtr->parm("JetMatching:coneRadius");
+  etaJetMax      = settingsPtr->parm("JetMatching:etaJetMax");
+  slowJetPower   = settingsPtr->mode("JetMatching:slowJetPower");
+
+  // Matching procedure
+  jetAllow       = settingsPtr->mode("JetMatching:jetAllow");
+  exclusiveMode  = settingsPtr->mode("JetMatching:exclusive");
+  qCutSq         = pow(qCut,2);
+  etaJetMaxAlgo  = etaJetMax;
+
+  // Read if veto should be performed internally.
+  performVeto    = settingsPtr->flag("JetMatching:doVeto");
+
+  // If not merging, then done
+  if (!doMerge) return true;
+
+  // Exclusive mode; if set to 2, then set based on nJet/nJetMax
+  if (exclusiveMode == 2) {
+
+    // No nJet or nJetMax, so default to exclusive mode
+    if (nJetMax < 0) {
+      errorMsg("Warning in JetMatchingMadgraph:init: "
+        "missing jet multiplicity information; running in exclusive mode");
+      exclusiveMode = 1;
+    }
+  }
+
+  // Initialise chosen jet algorithm.
+  // Currently, this only supports the kT-algorithm in SlowJet.
+  // Use the QCD distance measure by default.
+  jetAlgorithm = 2;
+  slowJetPower = 1;
+  slowJet = new SlowJet(slowJetPower, coneRadius, eTjetMin,
+    etaJetMaxAlgo, 2, 2, NULL, false);
+
+  // For FxFx, also initialise jet algorithm to define matrix element jets.
+  // Currently, this only supports the kT-algorithm in SlowJet.
+  // Use the QCD distance measure by default.
+  slowJetHard = new SlowJet(slowJetPower, coneRadius, qCutME,
+    etaJetMaxAlgo, 2, 2, NULL, false);
+
+  // To access the DJR's
+  slowJetDJR = new SlowJet(slowJetPower, coneRadius, qCutME,
+    etaJetMaxAlgo, 2, 2, NULL, false);
+
+  // A special version of SlowJet to handle heavy and other partons
+  hjSlowJet = new HJSlowJet(slowJetPower, coneRadius, 0.0,
+    100.0, 1, 2, NULL, false, true);
+
+  // Setup local event records
+  eventProcessOrig.init("(eventProcessOrig)", particleDataPtr);
+  eventProcess.init("(eventProcess)", particleDataPtr);
+  workEventJet.init("(workEventJet)", particleDataPtr);
+
+  // Print information
+  string jetStr  = (jetAlgorithm ==  1) ? "CellJet" :
+                   (slowJetPower == -1) ? "anti-kT" :
+                   (slowJetPower ==  0) ? "C/A"     :
+                   (slowJetPower ==  1) ? "kT"      : "unknown";
+  string modeStr = (exclusiveMode)         ? "exclusive" : "inclusive";
+  cout << endl
+       << " *-----  Madgraph matching parameters  -----*" << endl
+       << " |  qCut                |  " << setw(14)
+       << qCut << "  |" << endl
+       << " |  nQmatch             |  " << setw(14)
+       << nQmatch << "  |" << endl
+       << " |  clFact              |  " << setw(14)
+       << clFact << "  |" << endl
+       << " |  Jet algorithm       |  " << setw(14)
+       << jetStr << "  |" << endl
+       << " |  eTjetMin            |  " << setw(14)
+       << eTjetMin << "  |" << endl
+       << " |  etaJetMax           |  " << setw(14)
+       << etaJetMax << "  |" << endl
+       << " |  jetAllow            |  " << setw(14)
+       << jetAllow << "  |" << endl
+       << " |  Mode                |  " << setw(14)
+       << modeStr << "  |" << endl
+       << " *-----------------------------------------*" << endl;
+
+  return true;
+}
+
+//--------------------------------------------------------------------------
+
+// Process level vetos
+
+inline bool JetMatchingMadgraph::doVetoProcessLevel(Event& process) {
+
+  eventProcessOrig = process;
+
+  // Setup for veto if hard ME has too many partons.
+  // This is done to achieve consistency with the Pythia6 implementation.
+
+  // Clear the event of MPI systems and resonace decay products. Store trimmed
+  // event in workEvent.
+  sortIncomingProcess(process);
+
+  // Veto in case the hard input matrix element already has too many partons.
+  if ( !doFxFx && int(typeIdx[0].size()) > nJetMax )
+    return true;
+  if ( doFxFx && npNLO() < nJetMax && int(typeIdx[0].size()) > nJetMax )
+    return true;
+
+  // Done
+  return false;
+
+}
+
+//--------------------------------------------------------------------------
+
+inline bool JetMatchingMadgraph::doVetoStep(int iPos, int nISR, int nFSR,
+  const Event& event)  {
+
+  // Do not perform any veto if not in the Shower-kT scheme.
+  if ( !doShowerKt ) return false;
+
+  // Do nothing for emissions after the first one.
+  if ( nISR + nFSR > 1 ) return false;
+
+  // Do nothing in resonance decay showers.
+  if (iPos == 5) return false;
+
+  // Clear the event of MPI systems and resonace decay products. Store trimmed
+  // event in workEvent.
+  sortIncomingProcess(event);
+
+  // Get (kinematical) pT of first emission
+  double pTfirst = 0.;
+
+  // Get weak bosons, for later checks if the emission is a "QCD emission".
+  vector<int> weakBosons;
+  for (int i = 0; i < event.size(); i++) {
+    if ( event[i].id() == 22
+      && event[i].id() == 23
+      && event[i].idAbs() == 24)
+      weakBosons.push_back(i);
+  }
+
+  for (int i =  workEvent.size()-1; i > 0; --i) {
+    if ( workEvent[i].isFinal() && workEvent[i].colType() != 0
+      && (workEvent[i].statusAbs() == 43 || workEvent[i].statusAbs() == 51)) {
+      // Check if any of the EW bosons are ancestors of this parton. This
+      // should never happen for the first non-resonance shower emission.
+      // Check just to be sure.
+      bool QCDemission = true;
+      // Get position of this parton in the actual event (workEvent does
+      // not contain right mother-daughter relations). Stored in daughters.
+      int iPosOld = workEvent[i].daughter1();
+      for (int j = 0; i < int(weakBosons.size()); ++i)
+        if ( event[iPosOld].isAncestor(j)) {
+          QCDemission = false;
+          break;
+        }
+      // Done for a QCD emission.
+      if (QCDemission){
+        pTfirst = workEvent[i].pT();
+        break;
+      }
+    }
+  }
+
+  // Store things that are necessary to perform the shower-kT veto externally.
+  pTfirstSave   = pTfirst;
+  // Done if only inputs for an external vetoing procedure should be stored.
+  if (!performVeto) return false;
+
+  // Check veto.
+  if ( doShowerKtVeto(pTfirst) ) return true;
+
+  // No veto if come this far.
+  return false;
+
+}
+
+//--------------------------------------------------------------------------
+
+inline bool JetMatchingMadgraph::doShowerKtVeto(double pTfirst) {
+
+  // Only check veto in the shower-kT scheme.
+  if ( !doShowerKt ) return false;
+
+  // Reset veto code
+  bool doVeto = false;
+
+  // Find the (kinematical) pT of the softest (light) parton in the hard
+  // process.
+  int nParton = typeIdx[0].size();
+  double pTminME=1e10;
+  for ( int i = 0; i < nParton; ++i)
+    pTminME = min(pTminME,eventProcess[typeIdx[0][i]].pT());
+
+  // Veto if the softest hard process parton is below Qcut.
+  if ( nParton > 0 && pow(pTminME,2) < qCutSq ) doVeto = true;
+
+  // For non-highest multiplicity, veto if the hardest emission is harder
+  // than Qcut.
+  if ( exclusive && pow(pTfirst,2) > qCutSq ) {
+    doVeto = true;
+  // For highest multiplicity sample, veto if the hardest emission is harder
+  // than the hard process parton.
+  } else if ( !exclusive && nParton > 0 && pTfirst > pTminME ) {
+    doVeto = true;
+  }
+
+  // Return veto
+  return doVeto;
+
+}
+
+//--------------------------------------------------------------------------
+
+// Function to set the jet clustering scales (to be used as output)
+
+inline void JetMatchingMadgraph::setDJR( const Event& event) {
+
+ // Clear members.
+ clearDJR();
+ vector<double> result;
+
+  // Initialize SlowJetDJR jet algorithm with event
+  if (!slowJetDJR->setup(event) ) {
+    errorMsg("Warning in JetMatchingMadgraph:setDJR"
+      ": the SlowJet algorithm failed on setup");
+    return;
+  }
+
+  // Cluster in steps to find all hadronic jets
+  while ( slowJetDJR->sizeAll() - slowJetDJR->sizeJet() > 0 ) {
+    // Save the next clustering scale.
+    result.push_back(sqrt(slowJetDJR->dNext()));
+    // Perform step.
+    slowJetDJR->doStep();
+  }
+
+  // Save clustering scales in reserve order.
+  for (int i=int(result.size())-1; i >= 0; --i)
+    DJR.push_back(result[i]);
+
+}
+
+//--------------------------------------------------------------------------
+
+// Function to get the current number of partons in the Born state, as
+// read from LHE.
+
+inline int JetMatchingMadgraph::npNLO(){
+  string npIn = infoPtr->getEventAttribute("npNLO",true);
+  int np = (npIn != "") ? atoi((char*)npIn.c_str()) : -1;
+  if ( np < 0 ) { ; }
+  else return np;
+  return nPartonsNow;
+}
+
+//--------------------------------------------------------------------------
+
+// Step (1): sort the incoming particles
+
+inline void JetMatchingMadgraph::sortIncomingProcess(const Event &event) {
+
+  // Remove resonance decays from original process and keep only final
+  // state. Resonances will have positive status code after this step.
+  omitResonanceDecays(eventProcessOrig, true);
+  clearDJR();
+  clear_nMEpartons();
+
+  // Step-FxFx-1: remove preclustering from FxFx
+
+ /* 
+
+  // For FxFx, pre-cluster partons in the event into jets.
+  if (doFxFx) {
+
+    // Get final state partons
+    eventProcess.clear();
+    workEventJet.clear();
+    for( int i=0; i < workEvent.size(); ++i) {
+      // Original AG+Py6 algorithm explicitly excludes tops,
+      // leptons and photons.
+      int id = workEvent[i].idAbs();
+      if ((id >= ID_LEPMIN && id <= ID_LEPMAX) || id == ID_TOP
+        || id == ID_PHOTON || id == 23 || id == 24 || id == 25) {
+        eventProcess.append(workEvent[i]);
+      } else {
+        workEventJet.append(workEvent[i]);
+      }
+    }
+
+    // Initialize SlowJetHard jet algorithm with current working event
+    if (!slowJetHard->setup(workEventJet) ) {
+      errorMsg("Warning in JetMatchingMadgraph:sortIncomingProcess"
+        ": the SlowJet algorithm failed on setup");
+      return;
+    }
+
+    // Get matrix element cut scale.
+    double localQcutSq = qCutMESq;
+    // Cluster in steps to find all hadronic jets at the scale qCutME
+    while ( slowJetHard->sizeAll() - slowJetHard->sizeJet() > 0 ) {
+      // Done if next step is above qCut
+      if( slowJetHard->dNext() > localQcutSq ) break;
+      // Done if we're at or below the number of partons in the Born state.
+      if( slowJetHard->sizeAll()-slowJetHard->sizeJet() <= npNLO()) break;
+      slowJetHard->doStep();
+    }
+
+    // Construct a master copy of the event containing only the
+    // hardest nPartonsNow hadronic clusters. While constructing the event,
+    // the parton type (ID_GLUON) and status (98,99) are arbitrary.
+    int nJets = slowJetHard->sizeJet();
+    int nClus = slowJetHard->sizeAll();
+    int nNow = 0;
+    for (int i = nJets; i < nClus; ++i) {
+      vector<int> parts;
+      if (i < nClus-nJets) parts = slowJetHard->clusConstituents(i);
+      else parts = slowJetHard->constituents(nClus-nJets-i);
+      int flavour = ID_GLUON;
+      for(int j=0; j < int(parts.size()); ++j)
+        if (workEventJet[parts[j]].id() == ID_BOT)
+          flavour = ID_BOT;
+      eventProcess.append( flavour, 98,
+        workEventJet[parts.back()].mother1(),
+        workEventJet[parts.back()].mother2(),
+        workEventJet[parts.back()].daughter1(),
+        workEventJet[parts.back()].daughter2(),
+        0, 0, slowJetHard->p(i).px(), slowJetHard->p(i).py(),
+        slowJetHard->p(i).pz(), slowJetHard->p(i).e() );
+      nNow++;
+    }
+
+    // Done. Clean-up
+    workEventJet.clear();
+    
+   */
+
+  // Step-FxFx-2: : set the same eventProcess for both MLM and FxFx. This was only for MLM
+  // Add the type 2 selection also for FxFx
+
+  // For MLM matching, simply take hard process state from workEvent,
+  // without any preclustering.
+  //} else {
+    eventProcess = workEvent;
+  //}
+
+  // Sort original process final state into light/heavy jets and 'other'.
+  // Criteria:
+  //   1 <= ID <= nQmatch, or ID == 21         --> light jet (typeIdx[0])
+  //   nQMatch < ID                            --> heavy jet (typeIdx[1])
+  //   All else that is colored                --> other     (typeIdx[2])
+  // Note that 'typeIdx' stores indices into 'eventProcess' (after resonance
+  // decays are omitted), while 'typeSet' stores indices into the original
+  // process record, 'eventProcessOrig', but these indices are also valid
+  // in 'event'.
+  for (int i = 0; i < 3; i++) {
+    typeIdx[i].clear();
+    typeSet[i].clear();
+    origTypeIdx[i].clear();
+  }
+  for (int i = 0; i < eventProcess.size(); i++) {
+    // Ignore non-final state and default to 'other'
+    if (!eventProcess[i].isFinal()) continue;
+    int idx = -1;
+    int orig_idx = -1;
+
+    // Light jets: all gluons and quarks with id less than or equal to nQmatch
+    if (eventProcess[i].isGluon()
+      || (eventProcess[i].idAbs() <= nQmatch) ) {
+      orig_idx = 0;
+      if (doFxFx) {
+          // Crucial point FxFx: MG5 puts the scale of a not-to-be-matched quark 1 MeV lower than scalup. For
+          // such particles, we should keep the default "2"
+          idx = ( trunc(1000. * eventProcess[i].scale()) == trunc(1000. * infoPtr->scalup()) ) ? 0 : 2;
+
+      }
+      else {
+      // Crucial point: MG puts the scale of a non-QCD particle to eCM. For
+      // such particles, we should keep the default "2"
+      idx = ( eventProcess[i].scale() < 1.999 * sqrt(infoPtr->eA()
+        * infoPtr->eB()) ) ? 0 : 2;
+      }
+    }
+
+    // Heavy jets:  all quarks with id greater than nQmatch
+    else if (eventProcess[i].idAbs() > nQmatch
+      && eventProcess[i].idAbs() <= ID_TOP) {
+      idx = 1;
+      orig_idx = 1;
+    // Update to include non-SM colored particles
+    } else if (eventProcess[i].colType() != 0
+      && eventProcess[i].idAbs() > ID_TOP) {
+      idx = 1;
+      orig_idx = 1;
+    }
+    if( idx < 0 ) continue;
+    // Store
+    typeIdx[idx].push_back(i);
+    typeSet[idx].insert(eventProcess[i].daughter1());
+    origTypeIdx[orig_idx].push_back(i);
+  }
+
+  // Exclusive mode; if set to 2, then set based on nJet/nJetMax
+  if (exclusiveMode == 2) {
+
+    // Inclusive if nJet == nJetMax, exclusive otherwise
+    int nParton = origTypeIdx[0].size();
+    exclusive = (nParton == nJetMax) ? false : true;
+
+  // Otherwise, just set as given
+  } else {
+    exclusive = (exclusiveMode == 0) ? false : true;
+  }
+
+  // Extract partons from hardest subsystem + ISR + FSR only into
+  // workEvent. Note no resonance showers or MPIs.
+  subEvent(event);
+
+  // Store things that are necessary to perform the kT-MLM veto externally.
+  int nParton = typeIdx[0].size();
+  processSubsetSave.clear();
+  for ( int i = 0; i < nParton; ++i)
+    processSubsetSave.append( eventProcess[typeIdx[0][i]] );
+
+}
+
+//--------------------------------------------------------------------------
+
+// Step (2a): pick which particles to pass to the jet algorithm
+
+inline void JetMatchingMadgraph::jetAlgorithmInput(const Event &event,
+  int iType) {
+
+  // Take input from 'workEvent' and put output in 'workEventJet'
+  workEventJet = workEvent;
+
+  // Loop over particles and decide what to pass to the jet algorithm
+  for (int i = 0; i < workEventJet.size(); ++i) {
+    if (!workEventJet[i].isFinal()) continue;
+
+    // jetAllow option to disallow certain particle types
+    if (jetAllow == 1) {
+      // Remove all non-QCD partons from veto list
+      if( workEventJet[i].colType() == 0 ) {
+        workEventJet[i].statusNeg();
+        continue;
+      }
+    }
+
+    // Get the index of this particle in original event
+    int idx = workEventJet[i].daughter1();
+
+    // Start with particle idx, and afterwards track mothers
+    while (true) {
+
+      // Light jets
+      if (iType == 0) {
+
+        // Do not include if originates from heavy jet or 'other'
+        if (typeSet[1].find(idx) != typeSet[1].end() ||
+            typeSet[2].find(idx) != typeSet[2].end()) {
+          workEventJet[i].statusNeg();
+          break;
+        }
+
+        // Made it to start of event record so done
+        if (idx == 0) break;
+        // Otherwise next mother and continue
+        idx = event[idx].mother1();
+
+      // Heavy jets
+      } else if (iType == 1) {
+
+        // Only include if originates from heavy jet
+        if (typeSet[1].find(idx) != typeSet[1].end()) break;
+
+        // Made it to start of event record with no heavy jet mother,
+        // so DO NOT include particle
+        if (idx == 0) {
+          workEventJet[i].statusNeg();
+          break;
+        }
+
+        // Otherwise next mother and continue
+        idx = event[idx].mother1();
+
+      // Other jets
+      } else if (iType == 2) {
+
+        // Only include if originates from other jet
+        if (typeSet[2].find(idx) != typeSet[2].end()) break;
+
+        // Made it to start of event record with no heavy jet mother,
+        // so DO NOT include particle
+        if (idx == 0) {
+          workEventJet[i].statusNeg();
+          break;
+        }
+
+        // Otherwise next mother and continue
+        idx = event[idx].mother1();
+
+      } // if (iType)
+    } // while (true)
+  } // for (i)
+}
+
+//--------------------------------------------------------------------------
+
+// Step (2b): run jet algorithm and provide common output
+// This does nothing, because the jet algorithm is run several times
+//  in the matching algorithm.
+
+inline void JetMatchingMadgraph::runJetAlgorithm() {; }
+
+//--------------------------------------------------------------------------
+
+// Step (2c): veto decision (returning true vetoes the event)
+
+inline bool JetMatchingMadgraph::matchPartonsToJets(int iType) {
+
+  // Use different routines for light/heavy/other jets as
+  // different veto conditions and for clarity
+  if (iType == 0) {
+    // Record the jet separations here, also if matchPartonsToJetsLight
+    // returns preemptively.
+    setDJR(workEventJet);
+    set_nMEpartons(origTypeIdx[0].size(), typeIdx[0].size());
+    // Perform jet matching.
+    return (matchPartonsToJetsLight() > 0);
+  } else if (iType == 1) {
+     return (matchPartonsToJetsHeavy() > 0);
+  } else {
+     return (matchPartonsToJetsOther() > 0);
+  }
+
+}
+
+//--------------------------------------------------------------------------
+
+// Step(2c): light jets
+// Return codes are given indicating the reason for a veto.
+// Although not currently used, they are a useful debugging tool:
+//   0 = no veto
+//   1 = veto as number of jets less than number of partons
+//   2 = veto as exclusive mode and number of jets greater than
+//       number of partons
+//   3 = veto as inclusive mode and there would be an extra jet
+//       that is harder than any matched soft jet
+//   4 = veto as there is a parton which does not match a jet
+
+inline int JetMatchingMadgraph::matchPartonsToJetsLight() {
+
+  // Store things that are necessary to perform the kT-MLM veto externally.
+  workEventJetSave  = workEventJet;
+  // Done if only inputs for an external vetoing procedure should be stored.
+  if (!performVeto) return false;
+
+  // Count the number of hard partons
+  int nParton = typeIdx[0].size();
+
+  // Initialize SlowJet with current working event
+  if (!slowJet->setup(workEventJet) ) {
+    errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
+      "Light: the SlowJet algorithm failed on setup");
+    return NONE;
+  }
+  double localQcutSq = qCutSq;
+  double dOld = 0.0;
+  // Cluster in steps to find all hadronic jets at the scale qCut
+  while ( slowJet->sizeAll() - slowJet->sizeJet() > 0 ) {
+    if( slowJet->dNext() > localQcutSq ) break;
+    dOld = slowJet->dNext();
+    slowJet->doStep();
+  }
+  int nJets = slowJet->sizeJet();
+  int nClus = slowJet->sizeAll();
+
+  // Debug printout.
+  if (MATCHINGDEBUG) slowJet->list(true);
+
+  // Count of the number of hadronic jets in SlowJet accounting
+  int nCLjets = nClus - nJets;
+  // Get number of partons. Different for MLM and FxFx schemes.
+  //int nRequested = (doFxFx) ? npNLO() : nParton;
+  //Step-FxFx-3: Change nRequested subtracting the typeIdx[2] partons
+  //Exclude the highest multiplicity sample in the case of real emissions and all typeIdx[2]
+  //npNLO=multiplicity,nJetMax=njmax(shower_card),typeIdx[2]="Weak" jets
+  int nRequested = (doFxFx && !(npNLO()==nJetMax && npNLO()==(typeIdx[2].size()-1) )) ? npNLO()-typeIdx[2].size() : nParton;
+
+  //Step-FxFx-4:For FxFx veto the real emissions that have only typeIdx=2 partons
+  //From Step-FxFx-3 they already have negative nRequested, so this step may not be necessary
+  //Exclude the highest multiplicity sample
+  if (doFxFx && npNLO()<nJetMax && typeIdx[2].size()>0 && npNLO()==(typeIdx[2].size()-1)) {
+        return MORE_JETS;
+  }
+//----------------
+  //Exclude all Weak Jets for matching
+  //if (doFxFx && typeIdx[2].size()>0) {
+  //      return MORE_JETS;
+  //} // 2A
+  //keep only Weak Jets for matching
+  //if (doFxFx && typeIdx[2].size()==0) {
+  //      return MORE_JETS;
+  //} // 2B
+  //keep only lowest multiplicity sample @0
+  //if (doFxFx && npNLO()==1) {
+  //      return MORE_JETS;
+  //} // 2C
+  //keep only highest multiplicity sample @1
+  //if (doFxFx && npNLO()==0) {
+  //      return MORE_JETS;
+  //} // 2D
+//--------------
+  // Veto event if too few hadronic jets
+  if ( nCLjets < nRequested ) return LESS_JETS;
+
+  // In exclusive mode, do not allow more hadronic jets than partons
+  if ( exclusive && !doFxFx ) {
+    if ( nCLjets > nRequested ) return MORE_JETS;
+  } else {
+
+    // For FxFx, in the non-highest multipicity, all jets need to matched to
+    // partons. For nCLjets > nRequested, this is not possible. Hence, we can
+    // veto here already.
+    // if ( doFxFx && nRequested < nJetMax && nCLjets > nRequested )
+    //Step-FxFx-5:Change the nRequested to npNLO() in the first condition
+    //Before in Step-FxFx-3 it was nRequested=npNLO() for FxFx
+    //This way we correctly select the non-highest multipicity regardless the nRequested
+     if ( doFxFx && npNLO() < nJetMax && nCLjets > nRequested )
+      return MORE_JETS;
+
+    // Now continue in inclusive mode.
+    // In inclusive mode, there can be more hadronic jets than partons,
+    // provided that all partons are properly matched to hadronic jets.
+    // Start by setting up the jet algorithm.
+    if (!slowJet->setup(workEventJet) ) {
+      errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
+        "Light: the SlowJet algorithm failed on setup");
+      return NONE;
+    }
+
+    // For FxFx, continue clustering as long as the jet separation is above
+    // qCut.
+    if (doFxFx) {
+      while ( slowJet->sizeAll() - slowJet->sizeJet() > 0 ) {
+        if( slowJet->dNext() > localQcutSq ) break;
+        slowJet->doStep();
+      }
+    // For MLM, cluster into hadronic jets until there are the same number as
+    // partons.
+    } else {
+      while ( slowJet->sizeAll() - slowJet->sizeJet() > nParton )
+        slowJet->doStep();
+    }
+
+    // Sort partons in pT.  Update local qCut value.
+    //  Hadronic jets are already sorted in pT.
+    localQcutSq = dOld;
+    if ( clFact >= 0. && nParton > 0 ) {
+       vector<double> partonPt;
+       for (int i = 0; i < nParton; ++i)
+         partonPt.push_back( eventProcess[typeIdx[0][i]].pT2() );
+       sort( partonPt.begin(), partonPt.end());
+       localQcutSq = max( qCutSq, partonPt[0]);
+    }
+    nJets = slowJet->sizeJet();
+    nClus = slowJet->sizeAll();
+  }
+  // Update scale if clustering factor is non-zero
+  if ( clFact != 0. ) localQcutSq *= pow2(clFact);
+
+  Event tempEvent;
+  tempEvent.init( "(tempEvent)", particleDataPtr);
+  int nPass = 0;
+  double pTminEstimate = -1.;
+  // Construct a master copy of the event containing only the
+  // hardest nParton hadronic clusters. While constructing the event,
+  // the parton type (ID_GLUON) and status (98,99) are arbitrary.
+  for (int i = nJets; i < nClus; ++i) {
+    tempEvent.append( ID_GLUON, 98, 0, 0, 0, 0, 0, 0, slowJet->p(i).px(),
+      slowJet->p(i).py(), slowJet->p(i).pz(), slowJet->p(i).e() );
+    ++nPass;
+    pTminEstimate = max( pTminEstimate, slowJet->pT(i));
+    if(nPass == nRequested) break;
+  }
+
+  int tempSize = tempEvent.size();
+  // This keeps track of which hadronic jets are matched to parton
+  vector<bool> jetAssigned;
+  jetAssigned.assign( tempSize, false);
+
+  // This keeps track of which partons are matched to which hadronic
+  // jets.
+  vector< vector<bool> > partonMatchesJet;
+  for (int i=0; i < nParton; ++i )
+    partonMatchesJet.push_back( vector<bool>(tempEvent.size(),false) );
+
+  // Begin matching.
+  // Do jet matching for FxFx.
+  // Make sure that the nPartonsNow hardest hadronic jets are matched to any
+  // of the nPartonsNow (+1) partons. This matching is done by attaching a jet
+  // from the list of unmatched hadronic jets, and appending a jet from the
+  // list of partonic jets, one at a time. The partonic jet will be clustered
+  // with the hadronic jet or the beam if the distance measure is below the
+  // cut. The hadronic jet is matched once this happens. Otherwise, another
+  // partonic jet is tried. When a hadronic jet is matched to a partonic jet,
+  // it is removed from the list of unmatched hadronic jets. This process
+  // continues until the nPartonsNow hardest hadronic jets are matched to
+  // partonic jets, or it is not possible to make a match for a hadronic jet.
+  int iNow = 0;
+  int nMatched = 0;
+  while ( doFxFx && iNow < tempSize ) {
+
+    // Check if this shower jet matches any partonic jet.
+    Event tempEventJet;
+    tempEventJet.init("(tempEventJet)", particleDataPtr);
+    for (int i=0; i < nParton; ++i ) {
+
+      //// Only assign a parton once.
+      //for (int j=0; j < tempSize; ++j )
+      //  if ( partonMatchesJet[i][j]) continue;
+
+      // Attach a single hadronic jet.
+      tempEventJet.clear();
+      tempEventJet.append( ID_GLUON, 98, 0, 0, 0, 0, 0, 0,
+        tempEvent[iNow].px(), tempEvent[iNow].py(),
+        tempEvent[iNow].pz(), tempEvent[iNow].e() );
+      // Attach the current parton.
+      Vec4 pIn = eventProcess[typeIdx[0][i]].p();
+      tempEventJet.append( ID_GLUON, 99, 0, 0, 0, 0, 0, 0,
+        pIn.px(), pIn.py(), pIn.pz(), pIn.e() );
+
+      // Setup jet algorithm.
+      if ( !slowJet->setup(tempEventJet) ) {
+        errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
+          "Light: the SlowJet algorithm failed on setup");
+        return NONE;
+      }
+
+      // These are the conditions for the hadronic jet to match the parton
+      //  at the local qCut scale
+      if ( slowJet->iNext() == tempEventJet.size() - 1
+        && slowJet->jNext() > -1 && slowJet->dNext() < localQcutSq ) {
+        jetAssigned[iNow] = true;
+        partonMatchesJet[i][iNow] = true;
+      }
+
+    } // End loop over hard partons.
+
+    // Veto if the jet could not be assigned to any parton.
+    if ( jetAssigned[iNow] ) nMatched++;
+
+    // Continue;
+    ++iNow;
+  }
+
+  // Jet matching veto for FxFx
+  if (doFxFx) {
+//    if ( nRequested <  nJetMax && nMatched != nRequested )
+//      return UNMATCHED_PARTON;
+//    if ( nRequested == nJetMax && nMatched <  nRequested )
+//      return UNMATCHED_PARTON;
+      //Step-FxFx-6:Change the nRequested to npNLO() in the first condition (like in Step-FxFx-5)
+      //Before in Step-FxFx-3 it was nRequested=npNLO() for FxFx
+      //This way we correctly select the non-highest multipicity regardless the nRequested
+      if ( npNLO() <  nJetMax && nMatched != nRequested )
+        return UNMATCHED_PARTON;
+      if ( npNLO() == nJetMax && nMatched <  nRequested )
+        return UNMATCHED_PARTON;
+  }
+
+  // Do jet matching for MLM.
+  // Take the list of unmatched hadronic jets and append a parton, one at
+  // a time. The parton will be clustered with the "closest" hadronic jet
+  // or the beam if the distance measure is below the cut. When a hadronic
+  // jet is matched to a parton, it is removed from the list of unmatched
+  // hadronic jets. This process continues until all hadronic jets are
+  // matched to partons or it is not possible to make a match.
+  iNow = 0;
+  while (!doFxFx && iNow < nParton ) {
+    Event tempEventJet;
+    tempEventJet.init("(tempEventJet)", particleDataPtr);
+    for (int i = 0; i < tempSize; ++i) {
+      if (jetAssigned[i]) continue;
+      Vec4 pIn = tempEvent[i].p();
+      // Append unmatched hadronic jets
+      tempEventJet.append( ID_GLUON, 98, 0, 0, 0, 0, 0, 0,
+        pIn.px(), pIn.py(), pIn.pz(), pIn.e() );
+    }
+
+    Vec4 pIn = eventProcess[typeIdx[0][iNow]].p();
+    // Append the current parton
+    tempEventJet.append( ID_GLUON, 99, 0, 0, 0, 0, 0, 0,
+      pIn.px(), pIn.py(), pIn.pz(), pIn.e() );
+    if ( !slowJet->setup(tempEventJet) ) {
+      errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
+        "Light: the SlowJet algorithm failed on setup");
+      return NONE;
+    }
+    // These are the conditions for the hadronic jet to match the parton
+    //  at the local qCut scale
+    if ( slowJet->iNext() == tempEventJet.size() - 1
+      && slowJet->jNext() > -1 && slowJet->dNext() < localQcutSq ) {
+      int iKnt = -1;
+      for (int i = 0; i != tempSize; ++i) {
+        if (jetAssigned[i]) continue;
+        ++iKnt;
+        // Identify the hadronic jet that matches the parton
+        if (iKnt == slowJet->jNext() ) jetAssigned[i] = true;
+      }
+    } else {
+      return UNMATCHED_PARTON;
+    }
+    ++iNow;
+  }
+
+  // Minimal eT/pT (CellJet/SlowJet) of matched light jets.
+  // Needed later for heavy jet vetos in inclusive mode.
+  // This information is not used currently.
+  if (nParton > 0 && pTminEstimate > 0) eTpTlightMin = pTminEstimate;
+  else eTpTlightMin = -1.;
+
+  // Record the jet separations.
+  setDJR(workEventJet);
+
+  // No veto
+  return NONE;
+}
+
+//--------------------------------------------------------------------------
+
+// Step(2c): heavy jets
+// Return codes are given indicating the reason for a veto.
+// Although not currently used, they are a useful debugging tool:
+//   0 = no veto as there are no extra jets present
+//   1 = veto as in exclusive mode and extra jets present
+//   2 = veto as in inclusive mode and extra jets were harder
+//       than any matched light jet
+
+inline int JetMatchingMadgraph::matchPartonsToJetsHeavy() {
+
+  // Currently, heavy jets are unmatched
+  // If there are no extra jets, then accept
+  // jetMomenta is NEVER used by MadGraph and is always empty.
+  //  This check does nothing.
+  //  Rather, if there is any heavy flavor that is harder than
+  //  what is present at the LHE level, then the event should
+  //  be vetoed.
+
+  // if (jetMomenta.empty()) return NONE;
+  // Count the number of hard partons
+  int nParton = typeIdx[1].size();
+
+  Event tempEventJet(workEventJet);
+
+  double scaleF(1.0);
+  // Rescale the heavy partons that are from the hard process to
+  //  have pT=collider energy.   Soft/collinear gluons will cluster
+  //  onto them, leaving a remnant of hard emissions.
+  for( int i=0; i<nParton; ++i) {
+    scaleF = eventProcessOrig[0].e()/workEventJet[typeIdx[1][i]].pT();
+    tempEventJet[typeIdx[1][i]].rescale5(scaleF);
+  }
+
+  if (!hjSlowJet->setup(tempEventJet) ) {
+    errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
+             "Heavy: the SlowJet algorithm failed on setup");
+    return NONE;
+  }
+
+
+  while ( hjSlowJet->sizeAll() - hjSlowJet->sizeJet() > 0 ) {
+    if( hjSlowJet->dNext() > qCutSq ) break;
+    hjSlowJet->doStep();
+  }
+
+  int nCLjets(0);
+  // Count the number of clusters with pT>qCut.  This includes the
+  //  original hard partons plus any hard emissions.
+  for(int idx=0 ; idx< hjSlowJet->sizeAll(); ++idx) {
+    if( hjSlowJet->pT(idx) > sqrt(qCutSq) ) nCLjets++;
+  }
+
+  // Debug printout.
+  if (MATCHINGDEBUG) hjSlowJet->list(true);
+
+  // Count of the number of hadronic jets in SlowJet accounting
+  //  int nCLjets = nClus - nJets;
+  // Get number of partons. Different for MLM and FxFx schemes.
+  int nRequested = nParton;
+
+  // Veto event if too few hadronic jets
+  if ( nCLjets < nRequested ) {
+    if (MATCHINGDEBUG) cout << "veto : hvy  LESS_JETS " << endl;
+    if (MATCHINGDEBUG) cout << "nCLjets = " << nCLjets << "; nRequest = "
+      << nRequested << endl;
+    return LESS_JETS;
+  }
+
+  // In exclusive mode, do not allow more hadronic jets than partons
+  if ( exclusive ) {
+    if ( nCLjets > nRequested ) {
+      if (MATCHINGDEBUG) cout << "veto : excl hvy  MORE_JETS " << endl;
+      return MORE_JETS;
+    }
+  }
+
+  // No extra jets were present so no veto
+  return NONE;
+}
+
+//--------------------------------------------------------------------------
+
+// Step(2c): other jets
+// Return codes are given indicating the reason for a veto.
+// Although not currently used, they are a useful debugging tool:
+//   0 = no veto as there are no extra jets present
+//   1 = veto as in exclusive mode and extra jets present
+//   2 = veto as in inclusive mode and extra jets were harder
+//       than any matched light jet
+
+inline int JetMatchingMadgraph::matchPartonsToJetsOther() {
+
+  // Currently, heavy jets are unmatched
+  // If there are no extra jets, then accept
+  // jetMomenta is NEVER used by MadGraph and is always empty.
+  //  This check does nothing.
+  //  Rather, if there is any heavy flavor that is harder than
+  //  what is present at the LHE level, then the event should
+  //  be vetoed.
+
+  // if (jetMomenta.empty()) return NONE;
+  // Count the number of hard partons
+  int nParton = typeIdx[2].size();
+
+  Event tempEventJet(workEventJet);
+
+  double scaleF(1.0);
+  // Rescale the heavy partons that are from the hard process to
+  //  have pT=collider energy.   Soft/collinear gluons will cluster
+  //  onto them, leaving a remnant of hard emissions.
+  for( int i=0; i<nParton; ++i) {
+    scaleF = eventProcessOrig[0].e()/workEventJet[typeIdx[2][i]].pT();
+    tempEventJet[typeIdx[2][i]].rescale5(scaleF);
+  }
+
+  if (!hjSlowJet->setup(tempEventJet) ) {
+    errorMsg("Warning in JetMatchingMadgraph:matchPartonsToJets"
+             "Heavy: the SlowJet algorithm failed on setup");
+    return NONE;
+  }
+
+
+  while ( hjSlowJet->sizeAll() - hjSlowJet->sizeJet() > 0 ) {
+    if( hjSlowJet->dNext() > qCutSq ) break;
+    hjSlowJet->doStep();
+  }
+
+  int nCLjets(0);
+  // Count the number of clusters with pT>qCut.  This includes the
+  //  original hard partons plus any hard emissions.
+  for(int idx=0 ; idx< hjSlowJet->sizeAll(); ++idx) {
+    if( hjSlowJet->pT(idx) > sqrt(qCutSq) ) nCLjets++;
+  }
+
+  // Debug printout.
+  if (MATCHINGDEBUG) hjSlowJet->list(true);
+
+  // Count of the number of hadronic jets in SlowJet accounting
+  //  int nCLjets = nClus - nJets;
+  // Get number of partons. Different for MLM and FxFx schemes.
+  int nRequested = nParton;
+
+  // Veto event if too few hadronic jets
+  if ( nCLjets < nRequested ) {
+    if (MATCHINGDEBUG) cout << "veto : other LESS_JETS " << endl;
+    if (MATCHINGDEBUG) cout << "nCLjets = " << nCLjets << "; nRequest = "
+      << nRequested << endl;
+    return LESS_JETS;
+  }
+
+  // In exclusive mode, do not allow more hadronic jets than partons
+  if ( exclusive ) {
+    if ( nCLjets > nRequested ) {
+      if (MATCHINGDEBUG) cout << "veto : excl other MORE_JETS" << endl;
+      return MORE_JETS;
+    }
+  }
+
+  // No extra jets were present so no veto
+  return NONE;
+}
+
+//==========================================================================
+
+} // end namespace Pythia8
+
+#endif // end Pythia8_JetMatching_H
 
=== modified file 'Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA8.Script'
--- Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA8.Script	2021-04-15 15:17:31 +0000
+++ Template/NLO/MCatNLO/Scripts/MCatNLO_MadFKS_PYTHIA8.Script	2021-11-12 09:03:34 +0000
@@ -230,23 +230,33 @@
 
 if [ ! -d $PYTHIA8LOCATION/xmldoc ]
 then
- PY8VER="8.2"
- echo ' USING PYTHIA 8, version >= 8.200'
- chmod +x $PYTHIA8LOCATION/bin/pythia8-config
- read HEPMCINCLIB <<< $($PYTHIA8LOCATION/bin/pythia8-config --hepmc2)
+  if [ -f $PYTHIA8LOCATION/share/Pythia8/xmldoc/Version.xml ]
+  then
+    while read line; do
+    if [[ $line == *"Pythia:versionNumber"*"8.2"* ]]; then PY8VER="8.2" ; echo ' USING PYTHIA 8.2xy' ; fi
+    if [[ $line == *"Pythia:versionNumber"*"8.3"* ]]; then PY8VER="8.3" ; echo ' USING PYTHIA 8.3xy' ; fi
+    done < $PYTHIA8LOCATION/share/Pythia8/xmldoc/Version.xml
+  else
+    echo "Cannot determine Pythia8 version, stopping run"
+    exit -1
+  fi
+  chmod +x $PYTHIA8LOCATION/bin/pythia8-config
+  read HEPMCINCLIB <<< $($PYTHIA8LOCATION/bin/pythia8-config --hepmc2)
 else
- PY8VER="8.1"
- echo ' USING PYTHIA 8, version < 8.200'
+  PY8VER="8.1"
+  echo ' USING PYTHIA 8.1xy'
 fi
 
 if [ "$PY8UTI" != "" ]
 then
  if [ $PY8VER == "8.1" ]; then cp -rf $thisdir/srcPythia8/Pythia8.cc $thisdir/Pythia8.cc; fi
  if [ $PY8VER == "8.2" ]; then cp -rf $thisdir/srcPythia8/Pythia82.cc $thisdir/Pythia82.cc; fi
+ if [ $PY8VER == "8.3" ]; then cp -rf $thisdir/srcPythia8/Pythia83.cc $thisdir/Pythia83.cc; fi
  cp -rf $thisdir/srcPythia8/Makefile $thisdir/Makefile
 else
  if [ $PY8VER == "8.1" ]; then cp -rf $thisdir/srcPythia8/Pythia8_hep.cc $thisdir/Pythia8.cc; fi
  if [ $PY8VER == "8.2" ]; then cp -rf $thisdir/srcPythia8/Pythia82_hep.cc $thisdir/Pythia82.cc; fi
+ if [ $PY8VER == "8.3" ]; then cp -rf $thisdir/srcPythia8/Pythia83_hep.cc $thisdir/Pythia83.cc; fi
  cp -rf $thisdir/srcPythia8/Makefile_hep $thisdir/Makefile
 fi
 if [ -f $thisdir/Pythia8.exe ]
@@ -256,6 +266,7 @@
 
 if [ $PY8VER == "8.1" ]; then cp -rf $PYTHIA8LOCATION/xmldoc $thisdir; fi
 if [ $PY8VER == "8.2" ]; then cp -rf $PYTHIA8LOCATION/share/Pythia8/xmldoc $thisdir; fi
+if [ $PY8VER == "8.3" ]; then cp -rf $PYTHIA8LOCATION/share/Pythia8/xmldoc $thisdir; fi
 PY8UTIMAKE="PY8UTI="$PY8UTI
 FORFILEMAKE="FORFILE="$PY8UTI
 HEPMCMAKE="HEPMCINCLIB="$HEPMCINCLIB
@@ -263,6 +274,7 @@
 
 if [ $PY8VER == "8.1" ]; then make Pythia8 "$PY8UTIMAKE" "$FORFILEMAKE" "$LGFORTRANMAKE" "$LIBSMAKE" "$LIBSMAKEP" "$INCLMAKE"; fi
 if [ $PY8VER == "8.2" ]; then make Pythia82 "$PY8UTIMAKE" "$FORFILEMAKE" "$LGFORTRANMAKE" "$LIBSMAKE" "$LIBSMAKEP" "$INCLMAKE" "$HEPMCMAKE"; fi
+if [ $PY8VER == "8.3" ]; then make Pythia83 "$PY8UTIMAKE" "$FORFILEMAKE" "$LGFORTRANMAKE" "$LIBSMAKE" "$LIBSMAKEP" "$INCLMAKE" "$HEPMCMAKE"; fi
 
 cd $thisdir
 if [ ! -f $thisdir/Pythia8.exe ]
@@ -366,6 +378,7 @@
 ! 3) Beam-parameter settings.
 Beams:frameType = 4                ! LHE initialization
 Beams:LHEF = events.lhe            ! Input LHE file
+Beams:setProductionScalesFromLHEF=on ! Identify the Weak Jets
 
 ! 4) Switch on/off the key event-generation steps.
 EOF
 
=== modified file 'Template/NLO/MCatNLO/srcPythia8/Makefile'
--- Template/NLO/MCatNLO/srcPythia8/Makefile	2015-12-04 17:56:46 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Makefile	2021-11-12 09:03:34 +0000
@@ -51,6 +51,11 @@
 	  -L$(PYTHIA8LOCATION)/$(LIBDIR) -lpythia8 $(LIBGZIP) \
 	  -lHepMCfio $(LIBGFORTRAN) $(HEPMCINCLIB) $(EXTRAPATHS) $(EXTRALIBS) \
 
+Pythia83: $(FORFILE) $(PYTHIA8LOCATION)/$(LIBDIR)/libpythia8.a
+	  $(CXX) $(CXXFLAGS) -I$(PYTHIA8LOCATION)/$(INCDIR) \
+	  -I$(TOPDIR)/$(INCDIR) $(INCLOPTION) $@.cc -o Pythia8.exe $(FORFILE) \
+	  -L$(PYTHIA8LOCATION)/$(LIBDIR) -lpythia8 $(LIBGZIP) \
+	  -lHepMCfio $(LIBGFORTRAN) $(HEPMCINCLIB) $(EXTRAPATHS) $(EXTRALIBS) \
 
 .PHONY: clean
 clean:
 
=== modified file 'Template/NLO/MCatNLO/srcPythia8/Makefile_hep'
--- Template/NLO/MCatNLO/srcPythia8/Makefile_hep	2015-12-04 17:56:46 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Makefile_hep	2021-11-12 09:03:34 +0000
@@ -51,6 +51,12 @@
 	  -L$(PYTHIA8LOCATION)/$(LIBDIR) -lpythia8 $(LIBGZIP) \
 	  $(HEPMCINCLIB) $(EXTRAPATHS) $(EXTRALIBS)  \
 
+Pythia83: $(FORFILE) $(PYTHIA8LOCATION)/$(LIBDIR)/libpythia8.a
+	  $(CXX) $(CXXFLAGS) -I$(PYTHIA8LOCATION)/$(INCDIR) \
+	  -I$(TOPDIR)/$(INCDIR) $(INCLOPTION) $@.cc -o Pythia8.exe $(FORFILE) \
+	  -L$(PYTHIA8LOCATION)/$(LIBDIR) -lpythia8 $(LIBGZIP) \
+	  $(HEPMCINCLIB) $(EXTRAPATHS) $(EXTRALIBS)  \
+
 .PHONY: clean
 clean:
 	rm -rf *.exe
 
=== added file 'Template/NLO/MCatNLO/srcPythia8/Pythia83.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia83.cc	1970-01-01 00:00:00 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia83.cc	2021-11-12 09:03:34 +0000
@@ -0,0 +1,171 @@
+// Driver for Pythia 8. Reads an input file dynamically created on
+// the basis of the inputs specified in MCatNLO_MadFKS_PY8.Script 
+#include "Pythia8/Pythia.h"
+#include "Pythia8Plugins/HepMC2.h"
+#include "Pythia8Plugins/aMCatNLOHooks.h"
+#include "Pythia8Plugins/CombineMatchingInput.h"
+#include "HepMC/GenEvent.h"
+#include "HepMC/IO_GenEvent.h"
+#include "HepMC/IO_BaseClass.h"
+#include "HepMC/IO_HEPEVT.h"
+#include "HepMC/HEPEVT_Wrapper.h"
+#include "fstream"
+#include "LHEFRead.h"
+
+using namespace Pythia8;
+
+extern "C" {
+  extern struct {
+    double EVWGT;
+  } cevwgt_;
+}
+#define cevwgt cevwgt_
+
+extern "C" { 
+  void pyabeg_(int&,char(*)[50]);
+  void pyaend_(double&);
+  void pyanal_(int&,double(*));
+}
+
+int main() {
+  Pythia pythia;
+
+  int cwgtinfo_nn;
+  char cwgtinfo_weights_info[1024][50];
+  double cwgt_ww[1024];
+
+  string inputname="Pythia8.cmd",outputname="Pythia8.hep";
+
+  pythia.readFile(inputname.c_str());
+
+  //Create UserHooks pointer for the FxFX matching. Stop if it failed. Pass pointer to Pythia.
+  CombineMatchingInput combined;
+  //UserHooks* matching            = NULL;
+
+  string filename = pythia.word("Beams:LHEF");
+
+  MyReader read(filename);
+  read.lhef_read_wgtsinfo_(cwgtinfo_nn,cwgtinfo_weights_info);
+  pyabeg_(cwgtinfo_nn,cwgtinfo_weights_info);
+
+  int nAbort=10;
+  int nPrintLHA=1;
+  int iAbort=0;
+  int iPrintLHA=0;
+  int nstep=5000;
+  int iEventtot=pythia.mode("Main:numberOfEvents");
+  int iEventshower=pythia.mode("Main:spareMode1");
+  string evt_norm=pythia.word("Main:spareWord1");
+  int iEventtot_norm=iEventtot;
+  if (evt_norm != "sum"){
+    iEventtot_norm=1;
+  }
+
+  //FxFx merging
+  bool isFxFx=pythia.flag("JetMatching:doFxFx");
+  if (isFxFx) {
+    combined.setHook(pythia);
+    //matching = combined->getHook(pythia);
+    //if (!matching) {
+    //  std::cout << " Failed to initialise jet matching structures.\n"
+    //            << " Program stopped.";
+    //  return 1;
+    //}
+    //pythia.setUserHooksPtr(matching);
+    int nJmax=pythia.mode("JetMatching:nJetMax");
+    double Qcut=pythia.parm("JetMatching:qCut");
+    double PTcut=pythia.parm("JetMatching:qCutME");
+    if (Qcut <= PTcut || Qcut <= 0.) {
+      std::cout << " \n";
+      std::cout << "Merging scale (shower_card.dat) smaller than pTcut (run_card.dat)"
+		<< Qcut << " " << PTcut << "\n";
+      return 0;
+    }
+  }
+
+  // Initialise Pythia.
+  if (!pythia.init()) {
+    cout << "Error: could not initialise Pythia" << endl;
+    return 0;
+  };
+
+  HepMC::IO_BaseClass *_hepevtio;
+  HepMC::Pythia8ToHepMC ToHepMC;
+  HepMC::IO_GenEvent ascii_io(outputname.c_str(), std::ios::out);
+  double nSelected;
+  int nTry;
+  double norm;
+
+  // Cross section
+  double sigmaTotal  = 0.;
+  int iLHEFread=0;
+
+  for (int iEvent = 0; ; ++iEvent) {
+    if (!pythia.next()) {
+      if (++iAbort < nAbort) continue;
+      break;
+    }
+    // the number of events read by Pythia so far
+    nSelected=double(pythia.info.nSelected());
+    // normalisation factor for the default analyses defined in pyanal_
+    norm=iEventtot_norm*iEvent/nSelected;
+
+    if (nSelected >= iEventshower) break;
+    if (pythia.info.isLHA() && iPrintLHA < nPrintLHA) {
+      pythia.LHAeventList();
+      pythia.info.list();
+      pythia.process.list();
+      pythia.event.list();
+      ++iPrintLHA;
+    }
+
+    double evtweight = pythia.info.weight();
+    double normhepmc;
+    // Add the weight of the current event to the cross section.
+    normhepmc = 1. / double(iEventshower);
+    if (evt_norm != "sum") {
+      sigmaTotal  += evtweight*normhepmc;
+    } else {
+      sigmaTotal  += evtweight*normhepmc*iEventtot;
+    }
+
+    HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
+    ToHepMC.fill_next_event( pythia, hepmcevt );
+
+    //define the IO_HEPEVT
+    _hepevtio = new HepMC::IO_HEPEVT;
+    _hepevtio->write_event(hepmcevt);
+    
+    //event weight
+    cevwgt.EVWGT=hepmcevt->weights()[0];
+
+    //call the FORTRAN analysis for this event. First, make sure to
+    //re-synchronize the reading of the weights with the reading of
+    //the event. (They get desynchronised if an event was rejected).
+    nTry=pythia.info.nTried();
+    for (; iLHEFread<nTry ; ++iLHEFread) {
+      read.lhef_read_wgts_(cwgt_ww);
+    }
+    cwgt_ww[0]=cevwgt.EVWGT;
+    pyanal_(cwgtinfo_nn,cwgt_ww);
+
+    if (iEvent % nstep == 0 && iEvent >= 100){
+      pyaend_(norm);
+    }
+    delete hepmcevt;
+  }
+  pyaend_(norm);
+
+  pythia.stat();
+  if (isFxFx){
+    std::cout << " \n";
+    std::cout << "*********************************************************************** \n";
+    std::cout << "*********************************************************************** \n";
+    std::cout << "Cross section, including FxFx merging is: "
+	      << sigmaTotal << "\n";
+    std::cout << "*********************************************************************** \n";
+    std::cout << "*********************************************************************** \n";
+  }
+
+  return 0;
+}
 
=== added file 'Template/NLO/MCatNLO/srcPythia8/Pythia83_hep.cc'
--- Template/NLO/MCatNLO/srcPythia8/Pythia83_hep.cc	1970-01-01 00:00:00 +0000
+++ Template/NLO/MCatNLO/srcPythia8/Pythia83_hep.cc	2021-11-12 09:03:34 +0000
@@ -0,0 +1,121 @@
+// Driver for Pythia 8. Reads an input file dynamically created on
+// the basis of the inputs specified in MCatNLO_MadFKS_PY8.Script 
+#include "Pythia8/Pythia.h"
+#include "Pythia8Plugins/HepMC2.h"
+#include "Pythia8Plugins/aMCatNLOHooks.h"
+#include "Pythia8Plugins/CombineMatchingInput.h"
+#include "HepMC/GenEvent.h"
+#include "HepMC/IO_GenEvent.h"
+
+using namespace Pythia8;
+
+int main() {
+  Pythia pythia;
+
+  string inputname="Pythia8.cmd",outputname="Pythia8.hep";
+
+  pythia.readFile(inputname.c_str());
+
+  //Create UserHooks pointer for the FxFX matching. Stop if it failed. Pass pointer to Pythia.
+  CombineMatchingInput combined;
+  //UserHooks* matching            = NULL;
+
+  int nAbort=10;
+  int nPrintLHA=1;
+  int iAbort=0;
+  int iPrintLHA=0;
+  int iEventtot=pythia.mode("Main:numberOfEvents");
+  int iEventshower=pythia.mode("Main:spareMode1");
+  string evt_norm=pythia.word("Main:spareWord1");
+
+  //FxFx merging
+  bool isFxFx=pythia.flag("JetMatching:doFxFx");
+  if (isFxFx) {
+    combined.setHook(pythia);
+    //matching = combined->getHook(pythia);
+    //if (!matching) {
+    //  std::cout << " Failed to initialise jet matching structures.\n"
+    //            << " Program stopped.";
+    //  return 1;
+    //}
+    //pythia.setUserHooksPtr(matching);
+    int nJmax=pythia.mode("JetMatching:nJetMax");
+    double Qcut=pythia.parm("JetMatching:qCut");
+    double PTcut=pythia.parm("JetMatching:qCutME");
+    if (Qcut <= PTcut || Qcut <= 0.) {
+      std::cout << " \n";
+      std::cout << "Merging scale (shower_card.dat) smaller than pTcut (run_card.dat)"
+		<< Qcut << " " << PTcut << "\n";
+      return 0;
+    }
+  }
+
+  // Initialise Pythia.
+  if (!pythia.init()) {
+    cout << "Error: could not initialise Pythia" << endl;
+    return 0;
+  };
+
+  HepMC::Pythia8ToHepMC ToHepMC;
+  HepMC::IO_GenEvent ascii_io(outputname.c_str(), std::ios::out);
+  // Do not store cross section information, as this will be done manually.
+  ToHepMC.set_store_pdf(false);
+  ToHepMC.set_store_proc(false);
+  ToHepMC.set_store_xsec(false);
+
+  // Cross section an error.
+  double sigmaTotal  = 0.;
+  double errorTotal  = 0.;
+
+  for (int iEvent = 0; ; ++iEvent) {
+    if (!pythia.next()) {
+      if (++iAbort < nAbort) continue;
+      break;
+    }
+    // the number of events read by Pythia so far
+    int nSelected=pythia.info.nSelected();
+
+    if (nSelected > iEventshower) break;
+    if (pythia.info.isLHA() && iPrintLHA < nPrintLHA) {
+      pythia.LHAeventList();
+      pythia.info.list();
+      pythia.process.list();
+      pythia.event.list();
+      ++iPrintLHA;
+    }
+
+    HepMC::GenEvent* hepmcevt = new HepMC::GenEvent();
+    double evtweight = pythia.info.weight();
+    double normhepmc;
+    // ALWAYS NORMALISE HEPMC WEIGHTS TO SUM TO THE CROSS SECTION
+    if (evt_norm != "sum") {
+      normhepmc = 1. / double(iEventshower);
+    } else {
+      normhepmc = double(iEventtot) / double(iEventshower);
+    }
+    sigmaTotal += evtweight*normhepmc;
+    hepmcevt->weights().push_back(evtweight*normhepmc);
+    ToHepMC.fill_next_event( pythia, hepmcevt );
+    // Add the weight of the current event to the cross section.
+    // Report cross section to hepmc
+    HepMC::GenCrossSection xsec;
+    xsec.set_cross_section( sigmaTotal, pythia.info.sigmaErr() );
+    hepmcevt->set_cross_section( xsec );
+    // Write the HepMC event to file. Done with it.
+    ascii_io << hepmcevt;    
+    delete hepmcevt;
+  }
+
+  pythia.stat();
+  if (isFxFx){
+    std::cout << " \n";
+    std::cout << "*********************************************************************** \n";
+    std::cout << "*********************************************************************** \n";
+    std::cout << "Cross section, including FxFx merging is: "
+	      << sigmaTotal << "\n";
+    std::cout << "*********************************************************************** \n";
+    std::cout << "*********************************************************************** \n";
+  }
+
+  return 0;
+}
 
=== removed file 'Template/NLO/Source/genps.inc'
=== modified file 'Template/NLO/Source/make_opts.inc'
--- Template/NLO/Source/make_opts.inc	2020-09-16 12:13:39 +0000
+++ Template/NLO/Source/make_opts.inc	2021-11-12 09:03:34 +0000
@@ -93,23 +93,17 @@
   endif
   llhapdf = $(shell $(lhapdf) --libs)
   reweight_xsec_events_pdf_dummy=
-  # check if we need to activate c++11 (for lhapdf6.2)
-  ifeq ($(origin CXX),default)
-    ifeq ($lhapdfversion$lhapdfsubversion,62)
-      CXX=$(DEFAULT_CPP_COMPILER) -std=c++11
-    else
-      CXX=$(DEFAULT_CPP_COMPILER)
-    endif
-  endif	     
 else
   alfas_functions=alfas_functions
   llhapdf=
   reweight_xsec_events_pdf_dummy=reweight_xsec_events_pdf_dummy.o
-  # Set CXX unless it's defined by an environment variable
-  ifeq ($(origin CXX),default)
-    CXX=$(DEFAULT_CPP_COMPILER)
-  endif
-endif
+endif
+
+# Set CXX unless it's defined by an environment variable
+ifeq ($(origin CXX),default)
+   CXX=$(DEFAULT_CPP_COMPILER) -std=c++11
+endif
+
 
 # Option PineAPPL
 ifneq ($(pineappl),)
 
=== removed file 'Template/NLO/Source/maxconfigs.inc'
=== removed file 'Template/NLO/Source/maxparticles.inc'
=== modified file 'Template/NLO/SubProcesses/BinothLHA.f'
--- Template/NLO/SubProcesses/BinothLHA.f	2020-06-26 09:54:58 +0000
+++ Template/NLO/SubProcesses/BinothLHA.f	2021-11-12 09:03:34 +0000
@@ -8,6 +8,7 @@
       implicit none
       include "nexternal.inc"
       include "coupl.inc"
+      include "../../Source/MODEL/input.inc"
       include 'born_nhel.inc'
       double precision pi, zero,mone
       parameter (pi=3.1415926535897932385d0)
@@ -57,6 +58,7 @@
       integer getordpowfromindex_ml5
       integer orders_to_amp_split_pos
       logical, allocatable, save :: keep_order(:)
+      include 'a0Gmuconv.inc'
       include 'orders.inc'
       integer amp_orders(nsplitorders)
       integer split_amp_orders(nsplitorders), iamp
@@ -74,10 +76,17 @@
 
       logical updateloop
       common /to_updateloop/updateloop
+
+      integer get_n_tagged_photons
+      external get_n_tagged_photons
+      double precision get_virtual_a0Gmu_conv
+      external get_virtual_a0Gmu_conv 
+      integer ntagph, qed_pow_b
 c masses
       include 'pmass.inc'
       data nbad / 0 /
 
+
       IOErrCounter = 0
 c update the ren_scale for MadLoop and the couplings (should be the
 c Ellis-Sexton scale)
@@ -93,8 +102,7 @@
       single  = 0d0
       double  = 0d0
       born_hel_from_virt = 0d0
-C     reset the amp_split array
-      amp_split(1:amp_split_size) = 0d0
+C     reset the various arrays
       amp_split_finite_ML(1:amp_split_size) = 0d0
       amp_split_poles_ML(1:amp_split_size,1) = 0d0
       amp_split_poles_ML(1:amp_split_size,2) = 0d0
@@ -256,6 +264,32 @@
 c      endif
 c      virt_wgt=virt_wgt+conversion*born_wgt*ao2pi
 c======================================================================
+
+c======================================================================
+c If there are tagged photon and other particles in the process, 
+C one must use a mixed Gmu-alpha0 renormalisation.
+      ntagph = get_n_tagged_photons() 
+      if (ntagph.ne.0) then
+        do i = 1, AMP_SPLIT_SIZE_BORN 
+          call amp_split_pos_to_orders(i, amp_orders)
+          born_wgt = amp_split(i)
+          ! this is the number of powers of 'e' in the born
+          qed_pow_b = amp_orders(qed_pos)
+          amp_orders(qed_pos) = amp_orders(qed_pos) + 2
+          if (amp_orders(qed_pos).gt.NLO_ORDERS(qed_pos)) cycle
+
+          amp_split_poles_ML(orders_to_amp_split_pos(amp_orders),1) = 
+     $        amp_split_poles_ML(orders_to_amp_split_pos(amp_orders),1) +
+     $         get_virtual_a0Gmu_conv(qed_pow_b,ntagph,1,born_wgt)
+
+          amp_split_finite_ML(orders_to_amp_split_pos(amp_orders)) = 
+     $        amp_split_finite_ML(orders_to_amp_split_pos(amp_orders)) +
+     $         get_virtual_a0Gmu_conv(qed_pow_b,ntagph,0,born_wgt)
+
+        enddo
+      endif
+c======================================================================
+
 c
 c Check poles for the first PS points when doing MC over helicities, and
 c for all phase-space points when not doing MC over helicities. Skip
@@ -268,7 +302,7 @@
          polecheck_passed = .true.
          ! loop over the full result and each of the amp_split
          ! contribution
-         do iamp=0,amp_split_size
+         do iamp=1,amp_split_size
           ! skip 0 contributions in the amp_split array
             if (iamp.ne.0) then
                if (amp_split_poles_FKS(iamp,1).eq.0d0.and.
 
=== modified file 'Template/NLO/SubProcesses/add_write_info.f'
--- Template/NLO/SubProcesses/add_write_info.f	2018-05-01 13:46:54 +0000
+++ Template/NLO/SubProcesses/add_write_info.f	2021-11-12 09:03:34 +0000
@@ -92,8 +92,8 @@
       data (xdir(i),i=1,3) /0,0,1/
 
 c For (n+1)-body this is the configuration mapping
-      integer            mapconfig(0:lmaxconfigs), this_config
-      common/to_mconfigs/mapconfig, this_config
+      integer            this_config
+      common/to_mconfigs/this_config
 
 c For shifting QCD partons from zero to their mass-shell
       double precision x(99),p(0:3,99)
@@ -119,6 +119,14 @@
      &     ,fksfather_lhe(fks_configs) ,ipartner_lhe(fks_configs)
       common/cto_LHE1/iSorH_lhe,ifks_lhe,jfks_lhe,
      #                fksfather_lhe,ipartner_lhe
+c pt_clust string
+      integer need_matching(nexternal)
+      common /c_need_matching_to_write/ need_matching
+      double precision ptclus
+      CHARACTER temp*600,temp0*7,integ*1,float*18
+      CHARACTER integfour*4      
+      CHARACTER(LEN=1000) ptclusstring
+      common /c_ptclusstring/ ptclusstring
       include 'orders.inc'
       logical is_aorg(nexternal)
       common /c_is_aorg/is_aorg
@@ -581,6 +589,38 @@
             pb(j,ito(i))=pb(j,i)
          enddo
       enddo
+
+
+! write the ptclusstring that knows about which partons should be
+! considered in the MLM-like matching
+      if (ickkw.eq.3) then
+         if (nincoming.ne.2) then
+            write (*,*) 'Need to incoming particles with ickkw=3 '/
+     $           /'(add_write_info.f)'
+            stop
+         endif
+         temp0='<scales '
+         temp=''
+         do i=nincoming+1,nexpart
+            integfour=''
+            float=''
+            if (need_matching(i).ne.-1) then
+               ptclus=shower_scale
+            else
+               if (nincoming.ne.2) then
+                  write (*,*) 'need two incoming particles '/
+     $                 /'in add_write_info.f'
+                  stop 1
+               endif
+               ptclus=shower_scale-0.001d0
+            endif
+            Write(float,'(f16.5)') ptclus
+            write(integfour,'(i4)') ito(i)
+            temp=trim(temp)//' pt_clust_'//trim(adjustl(integfour))/
+     $           /'="'//trim(adjustl(float))//'"'
+         enddo
+         ptclusstring=trim(adjustl(temp0//trim(temp)//'></scales>'))
+      endif
 c
 c     Set the number of particles that needs to be written in event file
 c
@@ -593,108 +633,51 @@
       subroutine set_itree(iconfig,Hevents,itree,sprop_tree,pmass_tree
      &     ,pwidth_tree)
       implicit none
-      integer iconfig
+      integer iconfig,iconf
       logical Hevents
       include "genps.inc"
       include 'nexternal.inc'
       include "coupl.inc"
-      double precision ZERO
-      parameter (ZERO=0d0)
+      include 'nFKSconfigs.inc'
+      include 'real_from_born_configs.inc'
+      double precision pmass(-nexternal:0,lmaxconfigs,0:fks_configs)
+      double precision pwidth(-nexternal:0,lmaxconfigs,0:fks_configs)
+      integer iforest(2,-max_branch:-1,lmaxconfigs,0:fks_configs)
+      integer sprop(-max_branch:-1,lmaxconfigs,0:fks_configs)
+      integer tprid(-max_branch:-1,lmaxconfigs,0:fks_configs)
+      integer mapconfig(0:lmaxconfigs,0:fks_configs)
+      common /c_configurations/pmass,pwidth,iforest,sprop,tprid
+     $     ,mapconfig
       integer itree(2,-max_branch:-1)
       integer sprop_tree(-max_branch:-1)
-      include "born_conf.inc"
-      integer i,j,jj
-      double precision pmass(-nexternal:0,lmaxconfigs)
-      double precision pwidth(-nexternal:0,lmaxconfigs)
-      integer pow(-nexternal:0,lmaxconfigs)
+      integer i,j
       double precision pmass_tree(-nexternal:0)
       double precision pwidth_tree(-nexternal:0)
-      integer i_fks,j_fks
-      common/fks_indices/i_fks,j_fks
-      include "born_props.inc"
-c
+      INTEGER NFKSPROCESS
+      COMMON/C_NFKSPROCESS/NFKSPROCESS
 c Do not really care about t-channels: loop should just go to
-c nexternal-4, even though there is one more when there are t-channels
-c around
-      do j=-(nexternal-4),-1
-         do i=1,2
-            itree(i,j)=iforest(i,j,iconfig)
-         enddo
-         sprop_tree(j)=sprop(j,iconfig)
-         pmass_tree(j)=pmass(j,iconfig)
-         pwidth_tree(j)=pwidth(j,iconfig)
-      enddo
-c
-c When we are doing H-events, we need to add --when j_fks is final
-c state-- the s-channel branching fks_mother -> j_fks + i_fks.  When
-c j_fks is initial state, we need to add a 'bogus' t-channel splitting
-c to make sure that all the loops stop properly
-c
+c nexternal-3 (or 4), even though there is one more when there are
+c t-channels around
       if (Hevents) then
-c must re-label the external particles to get the correct daughters
-         if (i_fks.le.j_fks) then
-            write (*,*) 'ERROR: i_fks should be greater than j_fks'
-            stop
-         endif
+         iconf=real_from_born_conf(iconfig,nFKSprocess)
+         do j=-(nexternal-3),-1
+            do i=1,2
+               itree(i,j)=iforest(i,j,iconf,nFKSprocess)
+            enddo
+            sprop_tree(j)=sprop(j,iconf,nFKSprocess)
+            pmass_tree(j)=pmass(j,iconf,nFKSprocess)
+            pwidth_tree(j)=pwidth(j,iconf,nFKSprocess)
+         enddo
+      else
          do j=-(nexternal-4),-1
             do i=1,2
-               if ( itree(i,j).ge.i_fks ) then
-                  itree(i,j)=itree(i,j)+1
-               endif
+               itree(i,j)=iforest(i,j,iconfig,0)
             enddo
+            sprop_tree(j)=sprop(j,iconfig,0)
+            pmass_tree(j)=pmass(j,iconfig,0)
+            pwidth_tree(j)=pwidth(j,iconfig,0)
          enddo
-c
-         if (j_fks.gt.nincoming) then
-c we must add an extra s-channel. Easiest is to add it all the way at
-c the beginning. Therefore, relabel everything else first. Use the
-c original ones to make sure we are doing it correctly (and not
-c accidentally overwriting something).
-            do j=-(nexternal-4),-1
-               sprop_tree(j-1)=sprop(j,iconfig)
-               pmass_tree(j-1)=pmass(j,iconfig)
-               pwidth_tree(j-1)=pwidth(j,iconfig)
-               do i=1,2
-                  itree(i,j-1)=iforest(i,j,iconfig)
-c Also update the internal references
-                  if ( itree(i,j-1).lt. 0 ) then
-                     itree(i,j-1)=itree(i,j-1)-1
-                  endif
-               enddo
-            enddo
-c
-c Add the new s-channel
-c
-            itree(1,-1)=i_fks
-            itree(2,-1)=j_fks
-c This will never be an on-shell s-channel. Hence, give it some bogus values:
-            sprop_tree(-1)=0
-            pmass_tree(-1)=0d0
-            pwidth_tree(-1)=0d0
-c
-c We have to make sure that the fks_mother (which is equal to the j_fks
-c label of the Born) is replaced by the new s-channel
-c
-            do j=-(nexternal-3),-2 ! do not include the new s-channel
-               do i=1,2
-                  if ( itree(i,j).eq. j_fks ) then
-                     itree(i,j)=-1 ! reference to the new s-channel
-                  endif
-               enddo
-            enddo
-
-         else
-c j_fks is initial state
-            jj=-(nexternal-3)     ! Just add it at the end
-c setting itree to 1 (or 2) makes sure that the loops over s-channel
-c propagators will exit
-            itree(1,jj)=1
-            itree(2,jj)=1
-            sprop_tree(jj)=0
-            pmass_tree(jj)=0d0
-            pwidth_tree(jj)=0d0
-         endif
       endif
-
       return
       end
 
 
=== modified file 'Template/NLO/SubProcesses/analysis_lhe.f'
--- Template/NLO/SubProcesses/analysis_lhe.f	2018-06-04 16:53:37 +0000
+++ Template/NLO/SubProcesses/analysis_lhe.f	2021-11-12 09:03:34 +0000
@@ -87,7 +87,7 @@
 c Auxiliary quantities used when writing events
       integer kwgtinfo
       integer i_wgt, kk, ii, jj, n, nn
-      character*140 buff
+      character*1000 buff
       INTEGER MAXNUP
       PARAMETER (MAXNUP=500)
       INTEGER NUP,IDPRUP,IDUP(MAXNUP),ISTUP(MAXNUP),
 
=== removed file 'Template/NLO/SubProcesses/c_weight.inc.THIS'
=== modified file 'Template/NLO/SubProcesses/chooser_functions.f'
--- Template/NLO/SubProcesses/chooser_functions.f	2018-04-05 17:59:21 +0000
+++ Template/NLO/SubProcesses/chooser_functions.f	2021-11-12 09:03:34 +0000
@@ -8,25 +8,24 @@
       double precision ZERO
       parameter (ZERO=0d0)
       include 'maxparticles.inc'
-      include 'ngraphs.inc'
+      include 'maxconfigs.inc'
       integer i,j,k
       INTEGER NFKSPROCESS
       COMMON/C_NFKSPROCESS/NFKSPROCESS
-      integer iforest(2,-max_branch:-1,n_max_cg)
-      integer sprop(-max_branch:-1,n_max_cg)
-      integer tprid(-max_branch:-1,n_max_cg)
-      integer mapconfig(0:n_max_cg)
+      integer iforest(2,-max_branch:-1,lmaxconfigs)
+      integer sprop(-max_branch:-1,lmaxconfigs)
+      integer tprid(-max_branch:-1,lmaxconfigs)
+      integer mapconfig(0:lmaxconfigs)
       common/c_configs_inc/iforest,sprop,tprid,mapconfig
-      double precision prmass(-max_branch:nexternal,n_max_cg)
-      double precision prwidth(-max_branch:-1,n_max_cg)
-      integer prow(-max_branch:-1,n_max_cg)
+      double precision prmass(-max_branch:nexternal,lmaxconfigs)
+      double precision prwidth(-max_branch:-1,lmaxconfigs)