Merge lp:~maddevelopers/mg5amcnlo/2.6.5_ewa into lp:mg5amcnlo/lts

[Olivier Mattelaer]

only for the convenient diff for the moment

[lp:~maddevelopers/mg5amcnlo/2.6.5_ewa]

307. By olivier-mattelaer on 2020-08-13

merge with 2.7.4

Unmerged revisions

307. By olivier-mattelaer on 2020-08-13

merge with 2.7.4

306. By olivier-mattelaer on 2020-06-09

merge with 2.7.2

305. By olivier-mattelaer on 2019-08-22

merge with latest 2.6.7

304. By olivier-mattelaer on 2019-08-22

fix the helicity factor for the Z boson in EWA

303. By olivier-mattelaer on 2019-07-03

allow lhapdf + EWA for pW initial state

302. By olivier-mattelaer on 2019-05-24

merge with latest 2.6.6

301. By olivier-mattelaer on 2019-05-23

fixing some tests

300. By olivier-mattelaer on 2019-05-16

merge with latest 2.6.6

299. By olivier-mattelaer on 2019-05-08

implement EWA for proton collision --thanks to xiaoran

298. By olivier-mattelaer on 2019-05-08

allow EWA for proton

=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py	2020-06-16 19:24:01 +0000
+++ MadSpin/interface_madspin.py	2020-08-13 08:01:43 +0000
@@ -150,6 +150,7 @@
         self.mg5cmd = master_interface.MasterCmd()
         self.seed = None
         self.err_branching_ratio = 0
+        self.me_run_name = "" # Events diretory name where to stotre the events (used by madevent) not use internally
         
         
         if event_path:
@@ -580,11 +581,35 @@
         decay go > sq j
         launch
         '''
+        
+        self.parser_launch.print_help()
+
+    def parser_launch(self):
+        usage = """launch [-n RUN_NAME]   
+        """
+        parser = misc.OptionParser(usage=usage)
+        parser.add_option("-n", "--name",
+                  default="",
+                  help="When NOT run in standalone instruct MG5aMC where to store the events file")
+        return parser
+    
+    def parse_launch(self, line):
+        
+        args = self.split_arg(line)
+        return self.parser_launch().parse_args(args)
+        
 
     @misc.mute_logger()
     def do_launch(self, line):
         """end of the configuration launched the code"""
         
+        (options, args) = self.parse_launch(line)
+        
+        if options.name:
+            self.me_run_name = options.name # Only use by MG5aMC
+        else:
+            self.me_run_name = ''
+        
         if self.options["spinmode"] in ["none"]:
             return self.run_bridge(line)
         elif self.options["spinmode"] == "onshell":
 
=== modified file 'Template/LO/Cards/pythia8_card_default.dat'
--- Template/LO/Cards/pythia8_card_default.dat	2017-02-06 18:35:42 +0000
+++ Template/LO/Cards/pythia8_card_default.dat	2020-08-13 08:01:43 +0000
@@ -12,6 +12,8 @@
 ! -------------------------------------------------------------------
 ! Specify the HEPMC output of the Pythia8 shower. You can set it to:
 !   auto      : MG5aMC will automatically place it the run_<i> directory
+!   autoremove: MG5aMC will automatically remove the file at the end of the run.
+!                (usefull when running with Delphes)
 !   /dev/null : to turn off the HEPMC output.
 !   <path>    : to select where the HEPMC file must written. It will 
 !               therefore not be placed in the run_<i> directory. The
 
=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2020-05-27 07:52:34 +0000
+++ Template/LO/Cards/run_card.dat	2020-08-13 08:01:43 +0000
@@ -28,7 +28,7 @@
 #*********************************************************************
 # Collider type and energy                                           *
 # lpp: 0=No PDF, 1=proton, -1=antiproton, 2=photon from proton,      *
-#                                         3=photon from electron     *
+#                3=photon from electron, 4=photon from muon          *
 #*********************************************************************
      %(lpp1)s        = lpp1    ! beam 1 type 
      %(lpp2)s        = lpp2    ! beam 2 type
 
=== 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	2020-08-13 08:01:43 +0000
@@ -0,0 +1,1484 @@
+c/* ********************************************************* */
+c/*      Effective W Approximation Splitting Functions        */
+c/*      Pulled from Barger and Phillips, pg 276              */
+c/* ********************************************************* */
+c
+c List of functions:
+c ewa_wX_ByPID,ewa_wX,ewa_wT,ewa_w0
+c ewa_zX_ByPID,ewa_zT,ewa_z0
+c
+c Input for ewa_vX_ByPID:
+c x: longitudinal momentum fraction carried by boson v
+c q2max: scale^2 of v boson pdf
+c pol: polarization of v boson (pol = +1,-1,0)
+c ppid: pid of v boson's parent fermion
+c/* ********************************************************* */
+
+c/* ********************************************************* */
+c return ewa splitting function for w boson by parent PID
+      double precision function ewa_wX_ByPID(x,q2max,pol,ppid)
+      implicit none
+      integer pol,ppid
+      double precision x,q2max
+      double precision ewa_wX
+      double precision ewa_wp_p, ewa_wm_p
+      external ewa_wX, ewa_wp_p, ewa_wm_p
+
+      include 'ElectroweakFlux.inc'
+
+      if(q2max.lt.ewa_mw2) then
+c         write (*,*) 'ERROR: q2max below MW2',q2max,ewa_mw2 
+         ewa_wX_ByPID = 0d0
+         return
+      endif
+      if(x.lt.eps .or. x.gt.(1.d0-eps)) then
+c         write (*,*) 'ERROR: x out of range',x
+         ewa_wX_ByPID =0d0
+         return
+      endif
+
+c is parent a proton (unbroken)
+      if(ppid.eq.2212) then
+         ewa_wX_ByPID = ewa_wp_p(x,q2max,pol)
+         return
+      else if(ppid.eq.-2212) then
+         ewa_wX_ByPID = ewa_wm_p(x,q2max,pol)
+         return
+c  other particle
+      else
+         ewa_wX_ByPID = ewa_wX(x,q2max,pol)
+      endif
+      return
+      end
+c/* ********************************************************* */
+c return ewa splitting function for w boson
+      double precision function ewa_wX(x,q2max,pol)
+      implicit none
+      integer pol
+      double precision x,q2max
+      double precision ewa_wT,ewa_w0
+      external ewa_wT,ewa_w0
+
+      if(abs(pol).gt.0) then
+         ewa_wX = ewa_wT(x,q2max,pol)
+      else
+         ewa_wX = ewa_w0(x)
+      endif
+      return
+      end
+c/* ********************************************************* */
+c return ewa splitting function for z boson by parent PID
+      double precision function ewa_zX_ByPID(x,q2max,pol,ppid)
+      implicit none
+      integer pol,ppid
+      double precision x,q2max
+      double precision ewa_zT,ewa_z0, ewa_z_p
+      external ewa_zT,ewa_z0,ewa_z_p
+
+      include 'ElectroweakFlux.inc'
+
+      if(q2max.lt.ewa_mz2) then
+c         write (*,*) 'ERROR: q2max below MZ2',q2max,ewa_mz2
+         ewa_zX_ByPID = 0d0
+         return
+      endif
+      if(x.lt.eps .or. x.gt.(1.d0-eps))then
+c         write (*,*) 'ERROR: x out of range',x
+         ewa_zX_ByPID = 0d0
+         return
+      endif
+
+c is parent a charged lepton?
+      if(     abs(ppid).eq.11
+     & .or.   abs(ppid).eq.13
+     & .or.   abs(ppid).eq.15) then
+         if(abs(pol).gt.0) then
+            ewa_zX_ByPID = ewa_zT(x,q2max,pol,ewa_gVl,ewa_gAl)
+         else
+            ewa_zX_ByPID = ewa_z0(x,ewa_gVl,ewa_gAl)
+         endif
+         return
+c is parent a proton (unbroken)
+      elseif(ppid.eq.2212) then
+         ewa_zX_ByPID = ewa_Z_p(x,q2max,pol)
+         return
+c is parent a neutrino?
+      elseif( abs(ppid).eq.12
+     & .or.   abs(ppid).eq.14
+     & .or.   abs(ppid).eq.16) then
+         if(abs(pol).gt.0) then
+            ewa_zX_ByPID = ewa_zT(x,q2max,pol,ewa_gVv,ewa_gAv)
+         else
+            ewa_zX_ByPID = ewa_z0(x,ewa_gVv,ewa_gAv)
+            endif
+         return
+c is parent an up-type quark?
+      elseif( abs(ppid).eq.2
+     & .or.   abs(ppid).eq.4
+     & .or.   abs(ppid).eq.6) then
+         if(abs(pol).gt.0) then
+            ewa_zX_ByPID = ewa_zT(x,q2max,pol,ewa_gVu,ewa_gAu)
+         else
+            ewa_zX_ByPID = ewa_z0(x,ewa_gVu,ewa_gAu)
+         endif
+         return
+c is parent a down-type quark?
+      elseif( abs(ppid).eq.1
+     & .or.   abs(ppid).eq.3
+     & .or.   abs(ppid).eq.5) then
+         if(abs(pol).gt.0) then
+            ewa_zX_ByPID = ewa_zT(x,q2max,pol,ewa_gVd,ewa_gAd)
+         else
+            ewa_zX_ByPID = ewa_z0(x,ewa_gVd,ewa_gAd)
+         endif
+         return
+c or something else?
+      else
+         ewa_zX_ByPID = 0d0
+         return
+      endif
+      end
+c/* ********************************************************* */
+c w boson splitting function: longitudinal polarization
+      double precision function ewa_w0(x)
+      implicit none
+      double precision x
+      double precision coup
+      
+      include 'ElectroweakFlux.inc'
+      
+c     P_W(x,lambda=0) = (gW/4pi)**2 (1-x)/x
+      coup = ewa_gW2/(16d0*pi2)
+      ewa_w0 = coup * (1d0-x)/x
+      return
+      end
+c/* ********************************************************* */
+c z boson splitting function: longitudinal polarization
+      double precision function ewa_z0(x,gV,gA)
+      implicit none
+      double precision x,gV,gA
+      double precision coup
+
+      include 'ElectroweakFlux.inc'
+
+c     P_Z(x,lambda=0) = (gW2/cw2*4pi2) * (gV2 + gA2) (1-x)/x
+      coup = ewa_gW2/(ewa_cw2*4d0*pi2) * (gV**2 + gA**2)
+      ewa_z0 = coup * (1d0-x)/x
+      return
+      end
+c/* ********************************************************* */
+c w boson splitting function: transverse polarization
+      double precision function ewa_wT(x,q2max,pol)
+      implicit none
+      integer pol
+      double precision x,q2max
+      double precision coup,numer
+
+      include 'ElectroweakFlux.inc'
+      
+c     P_W(x,lambda=\pm) = coup * [(gV \mp gA)^2 + (gV \pm gA)^2 (1-x)^2] * log(Q2/MV2)
+      coup  = ewa_gW2/(8d0*16d0*pi2)
+      numer = (ewa_gV-pol*ewa_gA)**2 
+     &      + (ewa_gV+pol*ewa_gA)**2 * (1d0-x)**2
+      ewa_wT = coup * numer * log(q2max/ewa_mw2) / x
+      return
+      end
+c/* ********************************************************* */
+c z boson splitting function: transverse polarization
+      double precision function ewa_zT(x,q2max,pol,gV,gA)
+      implicit none
+      integer pol
+      double precision x,q2max,gV,gA
+      double precision coup,numer
+
+      include 'ElectroweakFlux.inc'
+      
+c     P_Z(x,lambda=\pm) = coup * [(gV \mp gA)^2 + (gV \pm gA)^2 (1-x)^2] * log(Q2/MV2)
+      coup = ewa_gW2/(ewa_cw2*16d0*pi2)
+      numer = (gV-pol*gA)**2
+     &      + (gV+pol*gA)**2 * (1d0-x)**2
+      ewa_zT = coup * numer * log(q2max/ewa_mz2) / x
+      return
+      end
+          FUNCTION dFINT(NARG,ARG,NENT,ENT,TABLE)
+C
+C   INTERPOLATION ROUTINE. AUTHOR C. LETERTRE.
+C   MODIFIED BY B. SCHORR, 1.07.1982.
+C
+          implicit real * 8 (a-h,o-z)
+          INTEGER   NENT(*)
+          REAL*8      ARG(*),   ENT(*),   TABLE(*)
+          INTEGER   INDEX(32)
+          REAL*8      WEIGHT(32)
+          LOGICAL   MFLAG,    RFLAG
+          DFINT  =  0.
+          IF(NARG .LT. 1  .OR.  NARG .GT. 5)  GOTO 300
+          LMAX      =  0
+          ISTEP     =  1
+          KNOTS     =  1
+          INDEX(1)  =  1
+          WEIGHT(1) =  1.
+          DO 100    N  =  1, NARG
+             X     =  ARG(N)
+             NDIM  =  NENT(N)
+             LOCA  =  LMAX
+             LMIN  =  LMAX + 1
+             LMAX  =  LMAX + NDIM
+             IF(NDIM .GT. 2)  GOTO 10
+             IF(NDIM .EQ. 1)  GOTO 100
+             H  =  X - ENT(LMIN)
+             IF(H .EQ. 0.)  GOTO 90
+             ISHIFT  =  ISTEP
+             IF(X-ENT(LMIN+1) .EQ. 0.)  GOTO 21
+             ISHIFT  =  0
+             ETA     =  H / (ENT(LMIN+1) - ENT(LMIN))
+             GOTO 30
+  10         LOCB  =  LMAX + 1
+  11         LOCC  =  (LOCA+LOCB) / 2
+             IF(X-ENT(LOCC))  12, 20, 13
+  12         LOCB  =  LOCC
+             GOTO 14
+  13         LOCA  =  LOCC
+  14         IF(LOCB-LOCA .GT. 1)  GOTO 11
+             LOCA    =  MIN0( MAX0(LOCA,LMIN), LMAX-1 )
+             ISHIFT  =  (LOCA - LMIN) * ISTEP
+             ETA     =  (X - ENT(LOCA)) / (ENT(LOCA+1) - ENT(LOCA))
+             GOTO 30
+  20         ISHIFT  =  (LOCC - LMIN) * ISTEP
+  21         DO 22  K  =  1, KNOTS
+                INDEX(K)  =  INDEX(K) + ISHIFT
+  22            CONTINUE
+             GOTO 90
+  30         DO 31  K  =  1, KNOTS
+                INDEX(K)         =  INDEX(K) + ISHIFT
+                INDEX(K+KNOTS)   =  INDEX(K) + ISTEP
+                WEIGHT(K+KNOTS)  =  WEIGHT(K) * ETA
+                WEIGHT(K)        =  WEIGHT(K) - WEIGHT(K+KNOTS)
+  31            CONTINUE
+             KNOTS  =  2*KNOTS
+  90         ISTEP  =  ISTEP * NDIM
+ 100         CONTINUE
+          DO 200    K  =  1, KNOTS
+             I  =  INDEX(K)
+             DFINT  =  DFINT + WEIGHT(K) * TABLE(I)
+ 200         CONTINUE
+          RETURN
+ 300      CALL KERMTR('E104.1',LGFILE,MFLAG,RFLAG)
+          IF(MFLAG) THEN
+             IF(LGFILE .EQ. 0) THEN
+                WRITE(*,1000) NARG
+             ELSE
+                WRITE(LGFILE,1000) NARG
+             ENDIF
+          ENDIF
+          IF(.NOT. RFLAG) CALL ABEND
+          RETURN
+1000      FORMAT( 7X, 25HFUNCTION DFINT ... NARG =,I6,
+     +              17H NOT WITHIN RANGE)
+          END
+          SUBROUTINE KERSET(ERCODE,LGFILE,LIMITM,LIMITR)
+                    PARAMETER(KOUNTE  =  28)
+          CHARACTER*6         ERCODE,   CODE(KOUNTE)
+          LOGICAL             MFLAG,    RFLAG
+          INTEGER             KNTM(KOUNTE),       KNTR(KOUNTE)
+          DATA      LOGF      /  0  /
+          DATA      CODE(1), KNTM(1), KNTR(1)  / 'C204.1', 100, 100 /
+          DATA      CODE(2), KNTM(2), KNTR(2)  / 'C204.2', 100, 100 /
+          DATA      CODE(3), KNTM(3), KNTR(3)  / 'C204.3', 100, 100 /
+          DATA      CODE(4), KNTM(4), KNTR(4)  / 'C205.1', 100, 100 /
+          DATA      CODE(5), KNTM(5), KNTR(5)  / 'C205.2', 100, 100 /
+          DATA      CODE(6), KNTM(6), KNTR(6)  / 'C205.3', 100, 100 /
+          DATA      CODE(7), KNTM(7), KNTR(7)  / 'C305.1', 100, 100 /
+          DATA      CODE(8), KNTM(8), KNTR(8)  / 'C308.1', 100, 100 /
+          DATA      CODE(9), KNTM(9), KNTR(9)  / 'C312.1', 100, 100 /
+          DATA      CODE(10),KNTM(10),KNTR(10) / 'C313.1', 100, 100 /
+          DATA      CODE(11),KNTM(11),KNTR(11) / 'C336.1', 100, 100 /
+          DATA      CODE(12),KNTM(12),KNTR(12) / 'C337.1', 100, 100 /
+          DATA      CODE(13),KNTM(13),KNTR(13) / 'C341.1', 100, 100 /
+          DATA      CODE(14),KNTM(14),KNTR(14) / 'D103.1', 100, 100 /
+          DATA      CODE(15),KNTM(15),KNTR(15) / 'D106.1', 100, 100 /
+          DATA      CODE(16),KNTM(16),KNTR(16) / 'D209.1', 100, 100 /
+          DATA      CODE(17),KNTM(17),KNTR(17) / 'D509.1', 100, 100 /
+          DATA      CODE(18),KNTM(18),KNTR(18) / 'E100.1', 100, 100 /
+          DATA      CODE(19),KNTM(19),KNTR(19) / 'E104.1', 100, 100 /
+          DATA      CODE(20),KNTM(20),KNTR(20) / 'E105.1', 100, 100 /
+          DATA      CODE(21),KNTM(21),KNTR(21) / 'E208.1', 100, 100 /
+          DATA      CODE(22),KNTM(22),KNTR(22) / 'E208.2', 100, 100 /
+          DATA      CODE(23),KNTM(23),KNTR(23) / 'F010.1', 100,   0 /
+          DATA      CODE(24),KNTM(24),KNTR(24) / 'F011.1', 100,   0 /
+          DATA      CODE(25),KNTM(25),KNTR(25) / 'F012.1', 100,   0 /
+          DATA      CODE(26),KNTM(26),KNTR(26) / 'F406.1', 100,   0 /
+          DATA      CODE(27),KNTM(27),KNTR(27) / 'G100.1', 100, 100 /
+          DATA      CODE(28),KNTM(28),KNTR(28) / 'G100.2', 100, 100 /
+          LOGF  =  LGFILE
+          IF(ERCODE .EQ. ' ')  THEN
+             L  =  0
+          ELSE
+             DO 10  L = 1, 6
+                IF(ERCODE(1:L) .EQ. ERCODE)  GOTO 12
+  10            CONTINUE
+  12         CONTINUE
+          ENDIF
+          DO 14     I  =  1, KOUNTE
+             IF(L .EQ. 0)  GOTO 13
+             IF(CODE(I)(1:L) .NE. ERCODE(1:L))  GOTO 14
+  13         KNTM(I)  =  LIMITM
+             KNTR(I)  =  LIMITR
+  14         CONTINUE
+          RETURN
+          ENTRY KERMTR(ERCODE,LOG,MFLAG,RFLAG)
+          LOG  =  LOGF
+          DO 20     I  =  1, KOUNTE
+             IF(ERCODE .EQ. CODE(I))  GOTO 21
+  20         CONTINUE
+          WRITE(*,1000)  ERCODE
+          CALL ABEND
+          RETURN
+  21      RFLAG  =  KNTR(I) .GE. 1
+          IF(RFLAG  .AND.  (KNTR(I) .LT. 100))  KNTR(I)  =  KNTR(I) - 1
+          MFLAG  =  KNTM(I) .GE. 1
+          IF(MFLAG  .AND.  (KNTM(I) .LT. 100))  KNTM(I)  =  KNTM(I) - 1
+          IF(.NOT. RFLAG)  THEN
+             IF(LOGF .LT. 1)  THEN
+                WRITE(*,1001)  CODE(I)
+             ELSE
+                WRITE(LOGF,1001)  CODE(I)
+             ENDIF
+          ENDIF
+          IF(MFLAG .AND. RFLAG)  THEN
+             IF(LOGF .LT. 1)  THEN
+                WRITE(*,1002)  CODE(I)
+             ELSE
+                WRITE(LOGF,1002)  CODE(I)
+             ENDIF
+          ENDIF
+          RETURN
+1000      FORMAT(' KERNLIB LIBRARY ERROR. ' /
+     +           ' ERROR CODE ',A6,' NOT RECOGNIZED BY KERMTR',
+     +           ' ERROR MONITOR. RUN ABORTED.')
+1001      FORMAT(/' ***** RUN TERMINATED BY CERN LIBRARY ERROR ',
+     +           'CONDITION ',A6)
+1002      FORMAT(/' ***** CERN LIBRARY ERROR CONDITION ',A6)
+          END
+      SUBROUTINE ABEND
+C
+C CERN PROGLIB# Z035    ABEND           .VERSION KERNVAX  1.10  811126
+
+      STOP '*** ABEND ***'
+      END
+      function ewa_wp_p(x,Q2,pol)
+      implicit none
+      real*8 ewa_wp_p
+      real*8 x,Q2
+      integer pol
+      real*8 tmp
+      real*8 y,z
+      real*8 ylow,yupp,zlow,zupp
+      parameter (ylow=-16.1181,yupp=-1e-07)
+      parameter (zlow=4.38683,zupp=13.8155)
+      real*8 ewa_wp_p_1p
+      real*8 ewa_wp_p_1m
+      real*8 ewa_wp_p_0
+      y=log(x)
+      z=0.5*log(Q2)
+      if(pol.eq.1)then
+        tmp=ewa_wp_p_1p(y,z)
+      else if(pol.eq.-1)then
+        tmp=ewa_wp_p_1m(y,z)
+      else
+        tmp=ewa_wp_p_0(y,z)
+      endif
+      ewa_wp_p=tmp/x*y*y
+      end
+      function ewa_wp_p_1p(y,z)
+      implicit none
+      real*8 ewa_wp_p_1p,y,z
+      integer narg,nny,nnz
+      parameter (narg=2)
+      parameter (nny=32)
+      parameter (nnz=8)
+      integer iny,inz,nent(narg)
+      real*8 tmp,dfint,ymap,zmap
+      real*8 arg(narg),ent(nny+nnz)
+      real*8 yv(nny),zv(nnz),gridv(nny,nnz)
+      logical firsttime
+      external dfint,ymap,zmap
+      data yv/
+     # -1.6118096e+01, -1.5094991e+01, -1.4105430e+01, -1.3149414e+01,
+     # -1.2226943e+01, -1.1338015e+01, -1.0482632e+01, -9.6607941e+00,
+     # -8.8725001e+00, -8.1177506e+00, -7.3965455e+00, -6.7088848e+00,
+     # -6.0547686e+00, -5.4341967e+00, -4.8471693e+00, -4.2936863e+00,
+     # -3.7737478e+00, -3.2873536e+00, -2.8345039e+00, -2.4151986e+00,
+     # -2.0294377e+00, -1.6772213e+00, -1.3585493e+00, -1.0734217e+00,
+     # -8.2183848e-01, -6.0379973e-01, -4.1930540e-01, -2.6835549e-01,
+     # -1.5095001e-01, -6.7088947e-02, -1.6772312e-02, -1.0000000e-07/
+      data zv/
+     #  4.3868276e+00,  5.7337823e+00,  7.0807370e+00,  8.4276917e+00,
+     #  9.7746464e+00,  1.1121601e+01,  1.2468556e+01,  1.3815511e+01/
+      data (gridv(iny, 1),iny=1,32)/
+     #  3.5423305e-07,  2.8295336e-07,  2.2903220e-07,  1.8758200e-07,
+     #  1.5559163e-07,  1.3038773e-07,  1.1025515e-07,  9.4027921e-08,
+     #  8.0829168e-08,  7.0061706e-08,  6.1294211e-08,  5.4182560e-08,
+     #  4.8450890e-08,  4.3888687e-08,  4.0319491e-08,  3.7580596e-08,
+     #  3.5520713e-08,  3.3952301e-08,  3.2625426e-08,  3.1187707e-08,
+     #  2.9181849e-08,  2.6154236e-08,  2.1900372e-08,  1.6640020e-08,
+     #  1.1099066e-08,  6.2780040e-09,  2.8842540e-09,  1.0025616e-09,
+     #  2.4219639e-10,  3.2430180e-11,  1.7122412e-12,  4.1237372e-12/
+      data (gridv(iny, 2),iny=1,32)/
+     #  9.6800285e-03,  7.5895745e-03,  6.0144129e-03,  4.8116391e-03,
+     #  3.8906828e-03,  3.1743361e-03,  2.6106310e-03,  2.1637757e-03,
+     #  1.8068388e-03,  1.5208169e-03,  1.2917371e-03,  1.1084546e-03,
+     #  9.6221855e-04,  8.4652524e-04,  7.5595982e-04,  6.8586986e-04,
+     #  6.3206533e-04,  5.8999236e-04,  5.5430903e-04,  5.1815672e-04,
+     #  4.7362998e-04,  4.1384517e-04,  3.3686600e-04,  2.4790947e-04,
+     #  1.5951917e-04,  8.6605399e-05,  3.7927081e-05,  1.2444648e-05,
+     #  2.7972986e-06,  3.3735215e-07,  1.3499999e-08,  3.8170248e-08/
+      data (gridv(iny, 3),iny=1,32)/
+     #  2.6802909e-02,  2.0666858e-02,  1.6084293e-02,  1.2619023e-02,
+     #  9.9957736e-03,  7.9819494e-03,  6.4206742e-03,  5.2026236e-03,
+     #  4.2458765e-03,  3.4918949e-03,  2.8976082e-03,  2.4290655e-03,
+     #  2.0601111e-03,  1.7714145e-03,  1.5472491e-03,  1.3744754e-03,
+     #  1.2417664e-03,  1.1376810e-03,  1.0498817e-03,  9.6384193e-04,
+     #  8.6434504e-04,  7.3958460e-04,  5.8805171e-04,  4.2143198e-04,
+     #  2.6317925e-04,  1.3808501e-04,  5.8100243e-05,  1.8168748e-05,
+     #  3.8451248e-06,  4.2503569e-07,  1.2748665e-08,  4.3387870e-08/
+      data (gridv(iny, 4),iny=1,32)/
+     #  5.2421053e-02,  3.9824965e-02,  3.0512558e-02,  2.3545386e-02,
+     #  1.8330098e-02,  1.4377269e-02,  1.1354532e-02,  9.0302263e-03,
+     #  7.2316381e-03,  5.8351965e-03,  4.7504014e-03,  3.9068753e-03,
+     #  3.2509450e-03,  2.7436667e-03,  2.3536072e-03,  2.0551591e-03,
+     #  1.8269264e-03,  1.6483845e-03,  1.4987780e-03,  1.3553122e-03,
+     #  1.1959351e-03,  1.0052207e-03,  7.8338419e-04,  5.4881055e-04,
+     #  3.3404786e-04,  1.7020660e-04,  6.9197187e-05,  2.0763953e-05,
+     #  4.1724093e-06,  4.2803444e-07,  9.3531305e-09,  3.9817765e-08/
+      data (gridv(iny, 5),iny=1,32)/
+     #  8.7314473e-02,  6.5462234e-02,  4.9470387e-02,  3.7629031e-02,
+     #  2.8860417e-02,  2.2292410e-02,  1.7332180e-02,  1.3567321e-02,
+     #  1.0692431e-02,  8.4897045e-03,  6.8007137e-03,  5.5036549e-03,
+     #  4.5068644e-03,  3.7444546e-03,  3.1638216e-03,  2.7231011e-03,
+     #  2.3880260e-03,  2.1270240e-03,  1.9097372e-03,  1.7046592e-03,
+     #  1.4833082e-03,  1.2275453e-03,  9.4002916e-04,  6.4560015e-04,
+     #  3.8423698e-04,  1.9081605e-04,  7.5273640e-05,  2.1785061e-05,
+     #  4.1829038e-06,  4.0198589e-07,  6.0711799e-09,  3.4349528e-08/
+      data (gridv(iny, 6),iny=1,32)/
+     #  1.3235511e-01,  9.8047579e-02,  7.3187600e-02,  5.4960335e-02,
+     #  4.1600160e-02,  3.1701979e-02,  2.4311750e-02,  1.8768166e-02,
+     #  1.4585527e-02,  1.1418854e-02,  9.0192598e-03,  7.1973844e-03,
+     #  5.8123298e-03,  4.7637364e-03,  3.9726338e-03,  3.3768371e-03,
+     #  2.9266403e-03,  2.5776574e-03,  2.2889023e-03,  2.0198147e-03,
+     #  1.7358250e-03,  1.4167343e-03,  1.0680290e-03,  7.2056798e-04,
+     #  4.2032019e-04,  2.0398263e-04,  7.8308169e-05,  2.1899198e-05,
+     #  4.0241531e-06,  3.6414048e-07,  3.7675413e-09,  2.9031813e-08/
+      data (gridv(iny, 7),iny=1,32)/
+     #  1.9057983e-01,  1.3957407e-01,  1.0297250e-01,  7.6391524e-02,
+     #  5.7105084e-02,  4.2965226e-02,  3.2523923e-02,  2.4779895e-02,
+     #  1.9003972e-02,  1.4680964e-02,  1.1442346e-02,  9.0104755e-03,
+     #  7.1810453e-03,  5.8100245e-03,  4.7852651e-03,  4.0196945e-03,
+     #  3.4450331e-03,  3.0020034e-03,  2.6377816e-03,  2.3022876e-03,
+     #  1.9551214e-03,  1.5745712e-03,  1.1692227e-03,  7.7537535e-04,
+     #  4.4359770e-04,  2.1050710e-04,  7.8557363e-05,  2.0871834e-05,
+     #  3.5743974e-06,  2.9762988e-07,  2.7844229e-09,  2.2782149e-08/
+      data (gridv(iny, 8),iny=1,32)/
+     #  2.6712405e-01,  1.9341146e-01,  1.4102973e-01,  1.0336364e-01,
+     #  7.6305708e-02,  5.6682239e-02,  4.2352056e-02,  3.1846115e-02,
+     #  2.4101065e-02,  1.8371384e-02,  1.4128575e-02,  1.0978416e-02,
+     #  8.6341061e-03,  6.8955860e-03,  5.6086588e-03,  4.6555113e-03,
+     #  3.9452818e-03,  3.4012072e-03,  2.9571084e-03,  2.5527100e-03,
+     #  2.1419401e-03,  1.7020208e-03,  1.2447953e-03,  8.1123628e-04,
+     #  4.5515586e-04,  2.1111600e-04,  7.6224092e-05,  1.8396998e-05,
+     #  2.7018652e-06,  1.8403044e-07,  1.1449565e-09,  1.2394449e-08/
+      data firsttime/.true./
+      save
+      if(firsttime)then
+        firsttime=.false.
+        nent(1)=nny
+        nent(2)=nnz
+        do iny=1,nny
+          ent(iny)=(yv(iny))
+        enddo
+        do inz=1,nnz
+          ent(nny+inz)=(zv(inz))
+        enddo
+      endif
+      arg(1)=(y)
+      arg(2)=(z)
+      tmp=dfint(narg,arg,nent,ent,gridv)
+      ewa_wp_p_1p=tmp
+      return
+      end
+      function ewa_wp_p_1m(y,z)
+      implicit none
+      real*8 ewa_wp_p_1m,y,z
+      integer narg,nny,nnz
+      parameter (narg=2)
+      parameter (nny=32)
+      parameter (nnz=8)
+      integer iny,inz,nent(narg)
+      real*8 tmp,dfint,ymap,zmap
+      real*8 arg(narg),ent(nny+nnz)
+      real*8 yv(nny),zv(nnz),gridv(nny,nnz)
+      logical firsttime
+      external dfint,ymap,zmap
+      data yv/
+     # -1.6118096e+01, -1.5094991e+01, -1.4105430e+01, -1.3149414e+01,
+     # -1.2226943e+01, -1.1338015e+01, -1.0482632e+01, -9.6607941e+00,
+     # -8.8725001e+00, -8.1177506e+00, -7.3965455e+00, -6.7088848e+00,
+     # -6.0547686e+00, -5.4341967e+00, -4.8471693e+00, -4.2936863e+00,
+     # -3.7737478e+00, -3.2873536e+00, -2.8345039e+00, -2.4151986e+00,
+     # -2.0294377e+00, -1.6772213e+00, -1.3585493e+00, -1.0734217e+00,
+     # -8.2183848e-01, -6.0379973e-01, -4.1930540e-01, -2.6835549e-01,
+     # -1.5095001e-01, -6.7088947e-02, -1.6772312e-02, -1.0000000e-07/
+      data zv/
+     #  4.3868276e+00,  5.7337823e+00,  7.0807370e+00,  8.4276917e+00,
+     #  9.7746464e+00,  1.1121601e+01,  1.2468556e+01,  1.3815511e+01/
+      data (gridv(iny, 1),iny=1,32)/
+     #  2.2352243e-07,  1.7777182e-07,  1.4298652e-07,  1.1622085e-07,
+     #  9.5328508e-08,  7.8759840e-08,  6.5499586e-08,  5.4802442e-08,
+     #  4.6136489e-08,  3.9104599e-08,  3.3386967e-08,  2.8725277e-08,
+     #  2.4905274e-08,  2.1731925e-08,  1.9030359e-08,  1.6631364e-08,
+     #  1.4386235e-08,  1.2169804e-08,  9.9111740e-09,  7.6254325e-09,
+     #  5.4230477e-09,  3.4799784e-09,  1.9601523e-09,  9.3903981e-10,
+     #  3.6862928e-10,  1.1297902e-10,  2.5188758e-11,  3.6724956e-12,
+     #  2.9351639e-13,  8.1701272e-15,  7.6912480e-17,  6.8729255e-27/
+      data (gridv(iny, 2),iny=1,32)/
+     #  5.9726888e-03,  4.6455606e-03,  3.6455128e-03,  2.8853282e-03,
+     #  2.3012642e-03,  1.8468005e-03,  1.4906859e-03,  1.2098779e-03,
+     #  9.8764085e-04,  8.1150430e-04,  6.7166904e-04,  5.6040282e-04,
+     #  4.7152850e-04,  3.9983688e-04,  3.4081959e-04,  2.9051542e-04,
+     #  2.4549043e-04,  2.0307075e-04,  1.6177240e-04,  1.2167096e-04,
+     #  8.4473943e-05,  5.2813192e-05,  2.8909650e-05,  1.3417511e-05,
+     #  5.0832560e-06,  1.4959648e-06,  3.1801753e-07,  4.3791104e-08,
+     #  3.2498348e-09,  8.1063886e-11,  6.0039546e-13,  6.3617388e-23/
+      data (gridv(iny, 3),iny=1,32)/
+     #  1.6219214e-02,  1.2379365e-02,  9.5193091e-03,  7.3739795e-03,
+     #  5.7510065e-03,  4.5098861e-03,  3.5553416e-03,  2.8172455e-03,
+     #  2.2446775e-03,  1.7999287e-03,  1.4539007e-03,  1.1841462e-03,
+     #  9.7322588e-04,  8.0691077e-04,  6.7339686e-04,  5.6275965e-04,
+     #  4.6677178e-04,  3.7925827e-04,  2.9674760e-04,  2.1907115e-04,
+     #  1.4910620e-04,  9.1228668e-05,  4.8765008e-05,  2.2043195e-05,
+     #  8.1072382e-06,  2.3064166e-06,  4.7120345e-07,  6.1862016e-08,
+     #  4.3153587e-09,  9.8155908e-11,  5.5917625e-13,  7.2313499e-23/
+      data (gridv(iny, 4),iny=1,32)/
+     #  3.1188696e-02,  2.3420186e-02,  1.7702079e-02,  1.3467734e-02,
+     #  1.0309782e-02,  7.9318460e-03,  6.1325773e-03,  4.7646388e-03,
+     #  3.7215358e-03,  2.9251760e-03,  2.3161962e-03,  1.8496758e-03,
+     #  1.4912341e-03,  1.2139639e-03,  9.9569551e-04,  8.1869337e-04,
+     #  6.6868465e-04,  5.3525022e-04,  4.1252936e-04,  2.9977979e-04,
+     #  2.0061564e-04,  1.2049765e-04,  6.3112103e-05,  2.7890159e-05,
+     #  1.0000143e-05,  2.7633522e-06,  5.4561380e-07,  6.8754820e-08,
+     #  4.5477918e-09,  9.5512862e-11,  4.0175364e-13,  6.6363321e-23/
+      data (gridv(iny, 5),iny=1,32)/
+     #  5.1182299e-02,  3.7889407e-02,  2.8215137e-02,  2.1136435e-02,
+     #  1.5925131e-02,  1.2054633e-02,  9.1677464e-03,  7.0050385e-03,
+     #  5.3803216e-03,  4.1584009e-03,  3.2378903e-03,  2.5431829e-03,
+     #  2.0176932e-03,  1.6173071e-03,  1.3073170e-03,  1.0602661e-03,
+     #  8.5475523e-04,  6.7549247e-04,  5.1388835e-04,  3.6836179e-04,
+     #  2.4290493e-04,  1.4356363e-04,  7.3865824e-05,  3.2002011e-05,
+     #  1.1221536e-05,  3.0227536e-06,  5.7924224e-07,  7.0416470e-08,
+     #  4.4455239e-09,  8.7027993e-11,  2.5174163e-13,  5.7249565e-23/
+      data (gridv(iny, 6),iny=1,32)/
+     #  7.6557895e-02,  5.5955076e-02,  4.1120172e-02,  3.0385490e-02,
+     #  2.2576036e-02,  1.6847439e-02,  1.2629288e-02,  9.5105660e-03,
+     #  7.1985244e-03,  5.4827229e-03,  4.2072308e-03,  3.2573133e-03,
+     #  2.5483598e-03,  2.0155007e-03,  1.6086767e-03,  1.2891809e-03,
+     #  1.0274915e-03,  8.0291084e-04,  6.0382880e-04,  4.2759547e-04,
+     #  2.7827791e-04,  1.6211285e-04,  8.2090171e-05,  3.4938634e-05,
+     #  1.2008139e-05,  3.1607346e-06,  5.8901727e-07,  6.9129330e-08,
+     #  4.1796650e-09,  7.6790396e-11,  1.4801921e-13,  4.8386670e-23/
+      data (gridv(iny, 7),iny=1,32)/
+     #  1.0885701e-01,  7.8595632e-02,  5.7044003e-02,  4.1606207e-02,
+     #  3.0510595e-02,  2.2464186e-02,  1.6611725e-02,  1.2338586e-02,
+     #  9.2105421e-03,  6.9185264e-03,  5.2361582e-03,  3.9989670e-03,
+     #  3.0873318e-03,  2.4109199e-03,  1.9013023e-03,  1.5064199e-03,
+     #  1.1876339e-03,  9.1812852e-04,  6.8290918e-04,  4.7797892e-04,
+     #  3.0715516e-04,  1.7645732e-04,  8.7991353e-05,  3.6807315e-05,
+     #  1.2399813e-05,  3.1850040e-06,  5.7189229e-07,  6.3166456e-08,
+     #  3.5998373e-09,  6.0994275e-11,  1.0798487e-13,  3.7970497e-23/
+      data (gridv(iny, 8),iny=1,32)/
+     #  1.5068353e-01,  1.0749621e-01,  7.7041842e-02,  5.5476031e-02,
+     #  4.0143008e-02,  2.9160969e-02,  2.1271437e-02,  1.5583408e-02,
+     #  1.1472391e-02,  8.4985083e-03,  6.3433700e-03,  4.7785683e-03,
+     #  3.6402566e-03,  2.8065350e-03,  2.1865589e-03,  1.7127763e-03,
+     #  1.3356279e-03,  1.0214490e-03,  7.5140048e-04,  5.1978083e-04,
+     #  3.2977123e-04,  1.8681197e-04,  9.1728471e-05,  3.7697570e-05,
+     #  1.2430821e-05,  3.1016518e-06,  5.2367878e-07,  5.0603141e-08,
+     #  2.5423963e-09,  3.4635376e-11,  4.1186342e-14,  2.0657555e-23/
+      data firsttime/.true./
+      save
+      if(firsttime)then
+        firsttime=.false.
+        nent(1)=nny
+        nent(2)=nnz
+        do iny=1,nny
+          ent(iny)=(yv(iny))
+        enddo
+        do inz=1,nnz
+          ent(nny+inz)=(zv(inz))
+        enddo
+      endif
+      arg(1)=(y)
+      arg(2)=(z)
+      tmp=dfint(narg,arg,nent,ent,gridv)
+      ewa_wp_p_1m=tmp
+      return
+      end
+      function ewa_wp_p_0(y,z)
+      implicit none
+      real*8 ewa_wp_p_0,y,z
+      integer narg,nny,nnz
+      parameter (narg=2)
+      parameter (nny=32)
+      parameter (nnz=8)
+      integer iny,inz,nent(narg)
+      real*8 tmp,dfint,ymap,zmap
+      real*8 arg(narg),ent(nny+nnz)
+      real*8 yv(nny),zv(nnz),gridv(nny,nnz)
+      logical firsttime
+      external dfint,ymap,zmap
+      data yv/
+     # -1.6118096e+01, -1.5094991e+01, -1.4105430e+01, -1.3149414e+01,
+     # -1.2226943e+01, -1.1338015e+01, -1.0482632e+01, -9.6607941e+00,
+     # -8.8725001e+00, -8.1177506e+00, -7.3965455e+00, -6.7088848e+00,
+     # -6.0547686e+00, -5.4341967e+00, -4.8471693e+00, -4.2936863e+00,
+     # -3.7737478e+00, -3.2873536e+00, -2.8345039e+00, -2.4151986e+00,
+     # -2.0294377e+00, -1.6772213e+00, -1.3585493e+00, -1.0734217e+00,
+     # -8.2183848e-01, -6.0379973e-01, -4.1930540e-01, -2.6835549e-01,
+     # -1.5095001e-01, -6.7088947e-02, -1.6772312e-02, -1.0000000e-07/
+      data zv/
+     #  4.3868276e+00,  5.7337823e+00,  7.0807370e+00,  8.4276917e+00,
+     #  9.7746464e+00,  1.1121601e+01,  1.2468556e+01,  1.3815511e+01/
+      data (gridv(iny, 1),iny=1,32)/
+     #  3.5187008e-03,  2.8034859e-03,  2.2605790e-03,  1.8430659e-03,
+     #  1.5186478e-03,  1.2616393e-03,  1.0560624e-03,  8.9012320e-04,
+     #  7.5538508e-04,  6.4582892e-04,  5.5671967e-04,  4.8420795e-04,
+     #  4.2519673e-04,  3.7700036e-04,  3.3717825e-04,  3.0339543e-04,
+     #  2.7329402e-04,  2.4447706e-04,  2.1469168e-04,  1.8222898e-04,
+     #  1.4664818e-04,  1.0943422e-04,  7.3828376e-05,  4.3723959e-05,
+     #  2.1983802e-05,  9.0102753e-06,  2.8448697e-06,  6.3299117e-07,
+     #  8.7487931e-08,  5.2831049e-09,  1.1568847e-10,  1.8415286e-15/
+      data (gridv(iny, 2),iny=1,32)/
+     #  5.2548138e-03,  4.0999760e-03,  3.2295174e-03,  2.5667153e-03,
+     #  2.0584664e-03,  1.6625476e-03,  1.3517402e-03,  1.1060358e-03,
+     #  9.1078460e-04,  7.5538896e-04,  6.3160108e-04,  5.3288902e-04,
+     #  4.5415547e-04,  3.9119341e-04,  3.4038971e-04,  2.9850573e-04,
+     #  2.6252377e-04,  2.2960220e-04,  1.9724634e-04,  1.6373150e-04,
+     #  1.2870167e-04,  9.3622734e-05,  6.1405954e-05,  3.5240336e-05,
+     #  1.7102418e-05,  6.7317976e-06,  2.0267789e-06,  4.2595999e-07,
+     #  5.4725748e-08,  2.9719851e-09,  4.8981810e-11,  9.4455659e-16/
+      data (gridv(iny, 3),iny=1,32)/
+     #  7.1876070e-03,  5.5076672e-03,  4.2546965e-03,  3.3120537e-03,
+     #  2.5990113e-03,  2.0524391e-03,  1.6307512e-03,  1.3034276e-03,
+     #  1.0481689e-03,  8.4882548e-04,  6.9291048e-04,  5.7083971e-04,
+     #  4.7522642e-04,  4.0016654e-04,  3.4076979e-04,  2.9285118e-04,
+     #  2.5271536e-04,  2.1707747e-04,  1.8321030e-04,  1.4933939e-04,
+     #  1.1513883e-04,  8.2005502e-05,  5.2541044e-05,  2.9373325e-05,
+     #  1.3840817e-05,  5.2672006e-06,  1.5241105e-06,  3.0542392e-07,
+     #  3.6915891e-08,  1.8347912e-09,  2.2114402e-11,  5.3685066e-16/
+      data (gridv(iny, 4),iny=1,32)/
+     #  9.2731430e-03,  6.9948236e-03,  5.3144732e-03,  4.0652903e-03,
+     #  3.1328132e-03,  2.4282484e-03,  1.8929862e-03,  1.4841263e-03,
+     #  1.1704474e-03,  9.2943341e-04,  7.4399977e-04,  6.0109301e-04,
+     #  4.9092271e-04,  4.0580385e-04,  3.3953903e-04,  2.8700023e-04,
+     #  2.4383801e-04,  2.0634913e-04,  1.7159359e-04,  1.3773971e-04,
+     #  1.0446314e-04,  7.3071554e-05,  4.5887405e-05,  2.5084061e-05,
+     #  1.1524443e-05,  4.2604502e-06,  1.1914778e-06,  2.2920043e-07,
+     #  2.6286612e-08,  1.2107299e-09,  1.0080946e-11,  3.2845426e-16/
+      data (gridv(iny, 5),iny=1,32)/
+     #  1.1477138e-02,  8.5381747e-03,  6.3939961e-03,  4.8178867e-03,
+     #  3.6556258e-03,  2.7887348e-03,  2.1392407e-03,  1.6500181e-03,
+     #  1.2800074e-03,  9.9975891e-04,  7.8717888e-04,  6.2562654e-04,
+     #  5.0280633e-04,  4.0921770e-04,  3.3736136e-04,  2.8121172e-04,
+     #  2.3579394e-04,  1.9702349e-04,  1.6176823e-04,  1.2814066e-04,
+     #  9.5802895e-05,  6.5966820e-05,  4.0705648e-05,  2.1818543e-05,
+     #  9.8043554e-06,  3.5336308e-06,  9.5911504e-07,  1.7800578e-07,
+     #  1.9496528e-08,  8.4002050e-10,  4.3204403e-12,  2.1251138e-16/
+      data (gridv(iny, 6),iny=1,32)/
+     #  1.3802612e-02,  1.0141026e-02,  7.4971127e-03,  5.5741476e-03,
+     #  4.1720271e-03,  3.1386776e-03,  2.3737894e-03,  1.8049673e-03,
+     #  1.3802137e-03,  1.0625971e-03,  8.2471236e-04,  6.4617450e-04,
+     #  5.1210721e-04,  4.1120934e-04,  3.3467950e-04,  2.7561780e-04,
+     #  2.2846603e-04,  1.8878732e-04,  1.5327653e-04,  1.1999194e-04,
+     #  8.8573452e-05,  6.0136140e-05,  3.6529669e-05,  1.9238247e-05,
+     #  8.4746327e-06,  2.9853000e-06,  7.8837418e-07,  1.4128313e-07,
+     #  1.4820095e-08,  6.0070190e-10,  1.7185758e-12,  1.4369011e-16/
+      data (gridv(iny, 7),iny=1,32)/
+     #  1.6432308e-02,  1.1930946e-02,  8.7133575e-03,  6.3970356e-03,
+     #  4.7264435e-03,  3.5091707e-03,  2.6187361e-03,  1.9645111e-03,
+     #  1.4819077e-03,  1.1254264e-03,  8.6166780e-04,  6.6606423e-04,
+     #  5.2091794e-04,  4.1296488e-04,  3.3203523e-04,  2.7031018e-04,
+     #  2.2163224e-04,  1.8120725e-04,  1.4555330e-04,  1.1266604e-04,
+     #  8.2152609e-05,  5.5024591e-05,  3.2922540e-05,  1.7044970e-05,
+     #  7.3639505e-06,  2.5349841e-06,  6.4817141e-07,  1.0938699e-07,
+     #  1.0789142e-08,  4.0267815e-10,  9.9738412e-13,  9.3965182e-17/
+      data (gridv(iny, 8),iny=1,32)/
+     #  1.9588443e-02,  1.4054534e-02,  1.0140342e-02,  7.3511754e-03,
+     #  5.3617513e-03,  3.9285711e-03,  2.8927241e-03,  2.1409034e-03,
+     #  1.5930991e-03,  1.1934393e-03,  9.0134688e-04,  6.8734040e-04,
+     #  5.3043871e-04,  4.1513075e-04,  3.2970284e-04,  2.6531895e-04,
+     #  2.1516289e-04,  1.7405250e-04,  1.3830950e-04,  1.0585100e-04,
+     #  7.6238215e-05,  5.0368640e-05,  2.9681423e-05,  1.5103253e-05,
+     #  6.3960840e-06,  2.1473293e-06,  5.2444493e-07,  7.7848289e-08,
+     #  6.7448840e-09,  2.0379331e-10,  2.3209743e-13,  4.3818080e-17/
+      data firsttime/.true./
+      save
+      if(firsttime)then
+        firsttime=.false.
+        nent(1)=nny
+        nent(2)=nnz
+        do iny=1,nny
+          ent(iny)=(yv(iny))
+        enddo
+        do inz=1,nnz
+          ent(nny+inz)=(zv(inz))
+        enddo
+      endif
+      arg(1)=(y)
+      arg(2)=(z)
+      tmp=dfint(narg,arg,nent,ent,gridv)
+      ewa_wp_p_0=tmp
+      return
+      end
+      function ewa_wm_p(x,Q2,pol)
+      implicit none
+      real*8 ewa_wm_p
+      real*8 x,Q2
+      integer pol
+      real*8 tmp
+      real*8 y,z
+      real*8 ylow,yupp,zlow,zupp
+      parameter (ylow=-16.1181,yupp=-1e-07)
+      parameter (zlow=4.38683,zupp=13.8155)
+      real*8 ewa_wm_p_1p
+      real*8 ewa_wm_p_1m
+      real*8 ewa_wm_p_0
+      y=log(x)
+      z=0.5*log(Q2)
+      if(pol.eq.1)then
+        tmp=ewa_wm_p_1p(y,z)
+      else if(pol.eq.-1)then
+        tmp=ewa_wm_p_1m(y,z)
+      else
+        tmp=ewa_wm_p_0(y,z)
+      endif
+      ewa_wm_p=tmp/x*y*y
+      end
+      function ewa_wm_p_1p(y,z)
+      implicit none
+      real*8 ewa_wm_p_1p,y,z
+      integer narg,nny,nnz
+      parameter (narg=2)
+      parameter (nny=32)
+      parameter (nnz=8)
+      integer iny,inz,nent(narg)
+      real*8 tmp,dfint,ymap,zmap
+      real*8 arg(narg),ent(nny+nnz)
+      real*8 yv(nny),zv(nnz),gridv(nny,nnz)
+      logical firsttime
+      external dfint,ymap,zmap
+      data yv/
+     # -1.6118096e+01, -1.5094991e+01, -1.4105430e+01, -1.3149414e+01,
+     # -1.2226943e+01, -1.1338015e+01, -1.0482632e+01, -9.6607941e+00,
+     # -8.8725001e+00, -8.1177506e+00, -7.3965455e+00, -6.7088848e+00,
+     # -6.0547686e+00, -5.4341967e+00, -4.8471693e+00, -4.2936863e+00,
+     # -3.7737478e+00, -3.2873536e+00, -2.8345039e+00, -2.4151986e+00,
+     # -2.0294377e+00, -1.6772213e+00, -1.3585493e+00, -1.0734217e+00,
+     # -8.2183848e-01, -6.0379973e-01, -4.1930540e-01, -2.6835549e-01,
+     # -1.5095001e-01, -6.7088947e-02, -1.6772312e-02, -1.0000000e-07/
+      data zv/
+     #  4.3868276e+00,  5.7337823e+00,  7.0807370e+00,  8.4276917e+00,
+     #  9.7746464e+00,  1.1121601e+01,  1.2468556e+01,  1.3815511e+01/
+      data (gridv(iny, 1),iny=1,32)/
+     #  3.5423305e-07,  2.8295336e-07,  2.2903220e-07,  1.8758200e-07,
+     #  1.5559163e-07,  1.3038773e-07,  1.1025515e-07,  9.4027921e-08,
+     #  8.0829168e-08,  7.0061706e-08,  6.1294211e-08,  5.4182560e-08,
+     #  4.8450890e-08,  4.3888687e-08,  4.0319491e-08,  3.7580596e-08,
+     #  3.5520713e-08,  3.3952301e-08,  3.2625426e-08,  3.1187707e-08,
+     #  2.9181849e-08,  2.6154236e-08,  2.1900372e-08,  1.6640020e-08,
+     #  1.1099066e-08,  6.2780040e-09,  2.8842540e-09,  1.0025616e-09,
+     #  2.4219639e-10,  3.2430180e-11,  1.7122412e-12,  4.1237372e-12/
+      data (gridv(iny, 2),iny=1,32)/
+     #  9.6800285e-03,  7.5895745e-03,  6.0144129e-03,  4.8116391e-03,
+     #  3.8906828e-03,  3.1743361e-03,  2.6106310e-03,  2.1637757e-03,
+     #  1.8068388e-03,  1.5208169e-03,  1.2917371e-03,  1.1084546e-03,
+     #  9.6221855e-04,  8.4652524e-04,  7.5595982e-04,  6.8586986e-04,
+     #  6.3206533e-04,  5.8999236e-04,  5.5430903e-04,  5.1815672e-04,
+     #  4.7362998e-04,  4.1384517e-04,  3.3686600e-04,  2.4790947e-04,
+     #  1.5951917e-04,  8.6605399e-05,  3.7927081e-05,  1.2444648e-05,
+     #  2.7972986e-06,  3.3735215e-07,  1.3499999e-08,  3.8170248e-08/
+      data (gridv(iny, 3),iny=1,32)/
+     #  2.6802909e-02,  2.0666858e-02,  1.6084293e-02,  1.2619023e-02,
+     #  9.9957736e-03,  7.9819494e-03,  6.4206742e-03,  5.2026236e-03,
+     #  4.2458765e-03,  3.4918949e-03,  2.8976082e-03,  2.4290655e-03,
+     #  2.0601111e-03,  1.7714145e-03,  1.5472491e-03,  1.3744754e-03,
+     #  1.2417664e-03,  1.1376810e-03,  1.0498817e-03,  9.6384193e-04,
+     #  8.6434504e-04,  7.3958460e-04,  5.8805171e-04,  4.2143198e-04,
+     #  2.6317925e-04,  1.3808501e-04,  5.8100243e-05,  1.8168748e-05,
+     #  3.8451248e-06,  4.2503569e-07,  1.2748665e-08,  4.3387870e-08/
+      data (gridv(iny, 4),iny=1,32)/
+     #  5.2421053e-02,  3.9824965e-02,  3.0512558e-02,  2.3545386e-02,
+     #  1.8330098e-02,  1.4377269e-02,  1.1354532e-02,  9.0302263e-03,
+     #  7.2316381e-03,  5.8351965e-03,  4.7504014e-03,  3.9068753e-03,
+     #  3.2509450e-03,  2.7436667e-03,  2.3536072e-03,  2.0551591e-03,
+     #  1.8269264e-03,  1.6483845e-03,  1.4987780e-03,  1.3553122e-03,
+     #  1.1959351e-03,  1.0052207e-03,  7.8338419e-04,  5.4881055e-04,
+     #  3.3404786e-04,  1.7020660e-04,  6.9197187e-05,  2.0763953e-05,
+     #  4.1724093e-06,  4.2803444e-07,  9.3531305e-09,  3.9817765e-08/
+      data (gridv(iny, 5),iny=1,32)/
+     #  8.7314473e-02,  6.5462234e-02,  4.9470387e-02,  3.7629031e-02,
+     #  2.8860417e-02,  2.2292410e-02,  1.7332180e-02,  1.3567321e-02,
+     #  1.0692431e-02,  8.4897045e-03,  6.8007137e-03,  5.5036549e-03,
+     #  4.5068644e-03,  3.7444546e-03,  3.1638216e-03,  2.7231011e-03,
+     #  2.3880260e-03,  2.1270240e-03,  1.9097372e-03,  1.7046592e-03,
+     #  1.4833082e-03,  1.2275453e-03,  9.4002916e-04,  6.4560015e-04,
+     #  3.8423698e-04,  1.9081605e-04,  7.5273640e-05,  2.1785061e-05,
+     #  4.1829038e-06,  4.0198589e-07,  6.0711799e-09,  3.4349528e-08/
+      data (gridv(iny, 6),iny=1,32)/
+     #  1.3235511e-01,  9.8047579e-02,  7.3187600e-02,  5.4960335e-02,
+     #  4.1600160e-02,  3.1701979e-02,  2.4311750e-02,  1.8768166e-02,
+     #  1.4585527e-02,  1.1418854e-02,  9.0192598e-03,  7.1973844e-03,
+     #  5.8123298e-03,  4.7637364e-03,  3.9726338e-03,  3.3768371e-03,
+     #  2.9266403e-03,  2.5776574e-03,  2.2889023e-03,  2.0198147e-03,
+     #  1.7358250e-03,  1.4167343e-03,  1.0680290e-03,  7.2056798e-04,
+     #  4.2032019e-04,  2.0398263e-04,  7.8308169e-05,  2.1899198e-05,
+     #  4.0241531e-06,  3.6414048e-07,  3.7675413e-09,  2.9031813e-08/
+      data (gridv(iny, 7),iny=1,32)/
+     #  1.9057983e-01,  1.3957407e-01,  1.0297250e-01,  7.6391524e-02,
+     #  5.7105084e-02,  4.2965226e-02,  3.2523923e-02,  2.4779895e-02,
+     #  1.9003972e-02,  1.4680964e-02,  1.1442346e-02,  9.0104755e-03,
+     #  7.1810453e-03,  5.8100245e-03,  4.7852651e-03,  4.0196945e-03,
+     #  3.4450331e-03,  3.0020034e-03,  2.6377816e-03,  2.3022876e-03,
+     #  1.9551214e-03,  1.5745712e-03,  1.1692227e-03,  7.7537535e-04,
+     #  4.4359770e-04,  2.1050710e-04,  7.8557363e-05,  2.0871834e-05,
+     #  3.5743974e-06,  2.9762988e-07,  2.7844229e-09,  2.2782149e-08/
+      data (gridv(iny, 8),iny=1,32)/
+     #  2.6712405e-01,  1.9341146e-01,  1.4102973e-01,  1.0336364e-01,
+     #  7.6305708e-02,  5.6682239e-02,  4.2352056e-02,  3.1846115e-02,
+     #  2.4101065e-02,  1.8371384e-02,  1.4128575e-02,  1.0978416e-02,
+     #  8.6341061e-03,  6.8955860e-03,  5.6086588e-03,  4.6555113e-03,
+     #  3.9452818e-03,  3.4012072e-03,  2.9571084e-03,  2.5527100e-03,
+     #  2.1419401e-03,  1.7020208e-03,  1.2447953e-03,  8.1123628e-04,
+     #  4.5515586e-04,  2.1111600e-04,  7.6224092e-05,  1.8396998e-05,
+     #  2.7018652e-06,  1.8403044e-07,  1.1449565e-09,  1.2394449e-08/
+      data firsttime/.true./
+      save
+      if(firsttime)then
+        firsttime=.false.
+        nent(1)=nny
+        nent(2)=nnz
+        do iny=1,nny
+          ent(iny)=(yv(iny))
+        enddo
+        do inz=1,nnz
+          ent(nny+inz)=(zv(inz))
+        enddo
+      endif
+      arg(1)=(y)
+      arg(2)=(z)
+      tmp=dfint(narg,arg,nent,ent,gridv)
+      ewa_wm_p_1p=tmp
+      return
+      end
+      function ewa_wm_p_1m(y,z)
+      implicit none
+      real*8 ewa_wm_p_1m,y,z
+      integer narg,nny,nnz
+      parameter (narg=2)
+      parameter (nny=32)
+      parameter (nnz=8)
+      integer iny,inz,nent(narg)
+      real*8 tmp,dfint,ymap,zmap
+      real*8 arg(narg),ent(nny+nnz)
+      real*8 yv(nny),zv(nnz),gridv(nny,nnz)
+      logical firsttime
+      external dfint,ymap,zmap
+      data yv/
+     # -1.6118096e+01, -1.5094991e+01, -1.4105430e+01, -1.3149414e+01,
+     # -1.2226943e+01, -1.1338015e+01, -1.0482632e+01, -9.6607941e+00,
+     # -8.8725001e+00, -8.1177506e+00, -7.3965455e+00, -6.7088848e+00,
+     # -6.0547686e+00, -5.4341967e+00, -4.8471693e+00, -4.2936863e+00,
+     # -3.7737478e+00, -3.2873536e+00, -2.8345039e+00, -2.4151986e+00,
+     # -2.0294377e+00, -1.6772213e+00, -1.3585493e+00, -1.0734217e+00,
+     # -8.2183848e-01, -6.0379973e-01, -4.1930540e-01, -2.6835549e-01,
+     # -1.5095001e-01, -6.7088947e-02, -1.6772312e-02, -1.0000000e-07/
+      data zv/
+     #  4.3868276e+00,  5.7337823e+00,  7.0807370e+00,  8.4276917e+00,
+     #  9.7746464e+00,  1.1121601e+01,  1.2468556e+01,  1.3815511e+01/
+      data (gridv(iny, 1),iny=1,32)/
+     #  2.2352243e-07,  1.7777182e-07,  1.4298652e-07,  1.1622085e-07,
+     #  9.5328508e-08,  7.8759840e-08,  6.5499586e-08,  5.4802442e-08,
+     #  4.6136489e-08,  3.9104599e-08,  3.3386967e-08,  2.8725277e-08,
+     #  2.4905274e-08,  2.1731925e-08,  1.9030359e-08,  1.6631364e-08,
+     #  1.4386235e-08,  1.2169804e-08,  9.9111740e-09,  7.6254325e-09,
+     #  5.4230477e-09,  3.4799784e-09,  1.9601523e-09,  9.3903981e-10,
+     #  3.6862928e-10,  1.1297902e-10,  2.5188758e-11,  3.6724956e-12,
+     #  2.9351639e-13,  8.1701272e-15,  7.6912480e-17,  6.8729255e-27/
+      data (gridv(iny, 2),iny=1,32)/
+     #  5.9726888e-03,  4.6455606e-03,  3.6455128e-03,  2.8853282e-03,
+     #  2.3012642e-03,  1.8468005e-03,  1.4906859e-03,  1.2098779e-03,
+     #  9.8764085e-04,  8.1150430e-04,  6.7166904e-04,  5.6040282e-04,
+     #  4.7152850e-04,  3.9983688e-04,  3.4081959e-04,  2.9051542e-04,
+     #  2.4549043e-04,  2.0307075e-04,  1.6177240e-04,  1.2167096e-04,
+     #  8.4473943e-05,  5.2813192e-05,  2.8909650e-05,  1.3417511e-05,
+     #  5.0832560e-06,  1.4959648e-06,  3.1801753e-07,  4.3791104e-08,
+     #  3.2498348e-09,  8.1063886e-11,  6.0039546e-13,  6.3617388e-23/
+      data (gridv(iny, 3),iny=1,32)/
+     #  1.6219214e-02,  1.2379365e-02,  9.5193091e-03,  7.3739795e-03,
+     #  5.7510065e-03,  4.5098861e-03,  3.5553416e-03,  2.8172455e-03,
+     #  2.2446775e-03,  1.7999287e-03,  1.4539007e-03,  1.1841462e-03,
+     #  9.7322588e-04,  8.0691077e-04,  6.7339686e-04,  5.6275965e-04,
+     #  4.6677178e-04,  3.7925827e-04,  2.9674760e-04,  2.1907115e-04,
+     #  1.4910620e-04,  9.1228668e-05,  4.8765008e-05,  2.2043195e-05,
+     #  8.1072382e-06,  2.3064166e-06,  4.7120345e-07,  6.1862016e-08,
+     #  4.3153587e-09,  9.8155908e-11,  5.5917625e-13,  7.2313499e-23/
+      data (gridv(iny, 4),iny=1,32)/
+     #  3.1188696e-02,  2.3420186e-02,  1.7702079e-02,  1.3467734e-02,
+     #  1.0309782e-02,  7.9318460e-03,  6.1325773e-03,  4.7646388e-03,
+     #  3.7215358e-03,  2.9251760e-03,  2.3161962e-03,  1.8496758e-03,
+     #  1.4912341e-03,  1.2139639e-03,  9.9569551e-04,  8.1869337e-04,
+     #  6.6868465e-04,  5.3525022e-04,  4.1252936e-04,  2.9977979e-04,
+     #  2.0061564e-04,  1.2049765e-04,  6.3112103e-05,  2.7890159e-05,
+     #  1.0000143e-05,  2.7633522e-06,  5.4561380e-07,  6.8754820e-08,
+     #  4.5477918e-09,  9.5512862e-11,  4.0175364e-13,  6.6363321e-23/
+      data (gridv(iny, 5),iny=1,32)/
+     #  5.1182299e-02,  3.7889407e-02,  2.8215137e-02,  2.1136435e-02,
+     #  1.5925131e-02,  1.2054633e-02,  9.1677464e-03,  7.0050385e-03,
+     #  5.3803216e-03,  4.1584009e-03,  3.2378903e-03,  2.5431829e-03,
+     #  2.0176932e-03,  1.6173071e-03,  1.3073170e-03,  1.0602661e-03,
+     #  8.5475523e-04,  6.7549247e-04,  5.1388835e-04,  3.6836179e-04,
+     #  2.4290493e-04,  1.4356363e-04,  7.3865824e-05,  3.2002011e-05,
+     #  1.1221536e-05,  3.0227536e-06,  5.7924224e-07,  7.0416470e-08,
+     #  4.4455239e-09,  8.7027993e-11,  2.5174163e-13,  5.7249565e-23/
+      data (gridv(iny, 6),iny=1,32)/
+     #  7.6557895e-02,  5.5955076e-02,  4.1120172e-02,  3.0385490e-02,
+     #  2.2576036e-02,  1.6847439e-02,  1.2629288e-02,  9.5105660e-03,
+     #  7.1985244e-03,  5.4827229e-03,  4.2072308e-03,  3.2573133e-03,
+     #  2.5483598e-03,  2.0155007e-03,  1.6086767e-03,  1.2891809e-03,
+     #  1.0274915e-03,  8.0291084e-04,  6.0382880e-04,  4.2759547e-04,
+     #  2.7827791e-04,  1.6211285e-04,  8.2090171e-05,  3.4938634e-05,
+     #  1.2008139e-05,  3.1607346e-06,  5.8901727e-07,  6.9129330e-08,
+     #  4.1796650e-09,  7.6790396e-11,  1.4801921e-13,  4.8386670e-23/
+      data (gridv(iny, 7),iny=1,32)/
+     #  1.0885701e-01,  7.8595632e-02,  5.7044003e-02,  4.1606207e-02,
+     #  3.0510595e-02,  2.2464186e-02,  1.6611725e-02,  1.2338586e-02,
+     #  9.2105421e-03,  6.9185264e-03,  5.2361582e-03,  3.9989670e-03,
+     #  3.0873318e-03,  2.4109199e-03,  1.9013023e-03,  1.5064199e-03,
+     #  1.1876339e-03,  9.1812852e-04,  6.8290918e-04,  4.7797892e-04,
+     #  3.0715516e-04,  1.7645732e-04,  8.7991353e-05,  3.6807315e-05,
+     #  1.2399813e-05,  3.1850040e-06,  5.7189229e-07,  6.3166456e-08,
+     #  3.5998373e-09,  6.0994275e-11,  1.0798487e-13,  3.7970497e-23/
+      data (gridv(iny, 8),iny=1,32)/
+     #  1.5068353e-01,  1.0749621e-01,  7.7041842e-02,  5.5476031e-02,
+     #  4.0143008e-02,  2.9160969e-02,  2.1271437e-02,  1.5583408e-02,
+     #  1.1472391e-02,  8.4985083e-03,  6.3433700e-03,  4.7785683e-03,
+     #  3.6402566e-03,  2.8065350e-03,  2.1865589e-03,  1.7127763e-03,
+     #  1.3356279e-03,  1.0214490e-03,  7.5140048e-04,  5.1978083e-04,
+     #  3.2977123e-04,  1.8681197e-04,  9.1728471e-05,  3.7697570e-05,
+     #  1.2430821e-05,  3.1016518e-06,  5.2367878e-07,  5.0603141e-08,
+     #  2.5423963e-09,  3.4635376e-11,  4.1186342e-14,  2.0657555e-23/
+      data firsttime/.true./
+      save
+      if(firsttime)then
+        firsttime=.false.
+        nent(1)=nny
+        nent(2)=nnz
+        do iny=1,nny
+          ent(iny)=(yv(iny))
+        enddo
+        do inz=1,nnz
+          ent(nny+inz)=(zv(inz))
+        enddo
+      endif
+      arg(1)=(y)
+      arg(2)=(z)
+      tmp=dfint(narg,arg,nent,ent,gridv)
+      ewa_wm_p_1m=tmp
+      return
+      end
+      function ewa_wm_p_0(y,z)
+      implicit none
+      real*8 ewa_wm_p_0,y,z
+      integer narg,nny,nnz
+      parameter (narg=2)
+      parameter (nny=32)
+      parameter (nnz=8)
+      integer iny,inz,nent(narg)
+      real*8 tmp,dfint,ymap,zmap
+      real*8 arg(narg),ent(nny+nnz)
+      real*8 yv(nny),zv(nnz),gridv(nny,nnz)
+      logical firsttime
+      external dfint,ymap,zmap
+      data yv/
+     # -1.6118096e+01, -1.5094991e+01, -1.4105430e+01, -1.3149414e+01,
+     # -1.2226943e+01, -1.1338015e+01, -1.0482632e+01, -9.6607941e+00,
+     # -8.8725001e+00, -8.1177506e+00, -7.3965455e+00, -6.7088848e+00,
+     # -6.0547686e+00, -5.4341967e+00, -4.8471693e+00, -4.2936863e+00,
+     # -3.7737478e+00, -3.2873536e+00, -2.8345039e+00, -2.4151986e+00,
+     # -2.0294377e+00, -1.6772213e+00, -1.3585493e+00, -1.0734217e+00,
+     # -8.2183848e-01, -6.0379973e-01, -4.1930540e-01, -2.6835549e-01,
+     # -1.5095001e-01, -6.7088947e-02, -1.6772312e-02, -1.0000000e-07/
+      data zv/
+     #  4.3868276e+00,  5.7337823e+00,  7.0807370e+00,  8.4276917e+00,
+     #  9.7746464e+00,  1.1121601e+01,  1.2468556e+01,  1.3815511e+01/
+      data (gridv(iny, 1),iny=1,32)/
+     #  3.5187008e-03,  2.8034859e-03,  2.2605790e-03,  1.8430659e-03,
+     #  1.5186478e-03,  1.2616393e-03,  1.0560624e-03,  8.9012320e-04,
+     #  7.5538508e-04,  6.4582892e-04,  5.5671967e-04,  4.8420795e-04,
+     #  4.2519673e-04,  3.7700036e-04,  3.3717825e-04,  3.0339543e-04,
+     #  2.7329402e-04,  2.4447706e-04,  2.1469168e-04,  1.8222898e-04,
+     #  1.4664818e-04,  1.0943422e-04,  7.3828376e-05,  4.3723959e-05,
+     #  2.1983802e-05,  9.0102753e-06,  2.8448697e-06,  6.3299117e-07,
+     #  8.7487931e-08,  5.2831049e-09,  1.1568847e-10,  1.8415286e-15/
+      data (gridv(iny, 2),iny=1,32)/
+     #  5.2548138e-03,  4.0999760e-03,  3.2295174e-03,  2.5667153e-03,
+     #  2.0584664e-03,  1.6625476e-03,  1.3517402e-03,  1.1060358e-03,
+     #  9.1078460e-04,  7.5538896e-04,  6.3160108e-04,  5.3288902e-04,
+     #  4.5415547e-04,  3.9119341e-04,  3.4038971e-04,  2.9850573e-04,
+     #  2.6252377e-04,  2.2960220e-04,  1.9724634e-04,  1.6373150e-04,
+     #  1.2870167e-04,  9.3622734e-05,  6.1405954e-05,  3.5240336e-05,
+     #  1.7102418e-05,  6.7317976e-06,  2.0267789e-06,  4.2595999e-07,
+     #  5.4725748e-08,  2.9719851e-09,  4.8981810e-11,  9.4455659e-16/
+      data (gridv(iny, 3),iny=1,32)/
+     #  7.1876070e-03,  5.5076672e-03,  4.2546965e-03,  3.3120537e-03,
+     #  2.5990113e-03,  2.0524391e-03,  1.6307512e-03,  1.3034276e-03,
+     #  1.0481689e-03,  8.4882548e-04,  6.9291048e-04,  5.7083971e-04,
+     #  4.7522642e-04,  4.0016654e-04,  3.4076979e-04,  2.9285118e-04,
+     #  2.5271536e-04,  2.1707747e-04,  1.8321030e-04,  1.4933939e-04,
+     #  1.1513883e-04,  8.2005502e-05,  5.2541044e-05,  2.9373325e-05,
+     #  1.3840817e-05,  5.2672006e-06,  1.5241105e-06,  3.0542392e-07,
+     #  3.6915891e-08,  1.8347912e-09,  2.2114402e-11,  5.3685066e-16/
+      data (gridv(iny, 4),iny=1,32)/
+     #  9.2731430e-03,  6.9948236e-03,  5.3144732e-03,  4.0652903e-03,
+     #  3.1328132e-03,  2.4282484e-03,  1.8929862e-03,  1.4841263e-03,
+     #  1.1704474e-03,  9.2943341e-04,  7.4399977e-04,  6.0109301e-04,
+     #  4.9092271e-04,  4.0580385e-04,  3.3953903e-04,  2.8700023e-04,
+     #  2.4383801e-04,  2.0634913e-04,  1.7159359e-04,  1.3773971e-04,
+     #  1.0446314e-04,  7.3071554e-05,  4.5887405e-05,  2.5084061e-05,
+     #  1.1524443e-05,  4.2604502e-06,  1.1914778e-06,  2.2920043e-07,
+     #  2.6286612e-08,  1.2107299e-09,  1.0080946e-11,  3.2845426e-16/
+      data (gridv(iny, 5),iny=1,32)/
+     #  1.1477138e-02,  8.5381747e-03,  6.3939961e-03,  4.8178867e-03,
+     #  3.6556258e-03,  2.7887348e-03,  2.1392407e-03,  1.6500181e-03,
+     #  1.2800074e-03,  9.9975891e-04,  7.8717888e-04,  6.2562654e-04,
+     #  5.0280633e-04,  4.0921770e-04,  3.3736136e-04,  2.8121172e-04,
+     #  2.3579394e-04,  1.9702349e-04,  1.6176823e-04,  1.2814066e-04,
+     #  9.5802895e-05,  6.5966820e-05,  4.0705648e-05,  2.1818543e-05,
+     #  9.8043554e-06,  3.5336308e-06,  9.5911504e-07,  1.7800578e-07,
+     #  1.9496528e-08,  8.4002050e-10,  4.3204403e-12,  2.1251138e-16/
+      data (gridv(iny, 6),iny=1,32)/
+     #  1.3802612e-02,  1.0141026e-02,  7.4971127e-03,  5.5741476e-03,
+     #  4.1720271e-03,  3.1386776e-03,  2.3737894e-03,  1.8049673e-03,
+     #  1.3802137e-03,  1.0625971e-03,  8.2471236e-04,  6.4617450e-04,
+     #  5.1210721e-04,  4.1120934e-04,  3.3467950e-04,  2.7561780e-04,
+     #  2.2846603e-04,  1.8878732e-04,  1.5327653e-04,  1.1999194e-04,
+     #  8.8573452e-05,  6.0136140e-05,  3.6529669e-05,  1.9238247e-05,
+     #  8.4746327e-06,  2.9853000e-06,  7.8837418e-07,  1.4128313e-07,
+     #  1.4820095e-08,  6.0070190e-10,  1.7185758e-12,  1.4369011e-16/
+      data (gridv(iny, 7),iny=1,32)/
+     #  1.6432308e-02,  1.1930946e-02,  8.7133575e-03,  6.3970356e-03,
+     #  4.7264435e-03,  3.5091707e-03,  2.6187361e-03,  1.9645111e-03,
+     #  1.4819077e-03,  1.1254264e-03,  8.6166780e-04,  6.6606423e-04,
+     #  5.2091794e-04,  4.1296488e-04,  3.3203523e-04,  2.7031018e-04,
+     #  2.2163224e-04,  1.8120725e-04,  1.4555330e-04,  1.1266604e-04,
+     #  8.2152609e-05,  5.5024591e-05,  3.2922540e-05,  1.7044970e-05,
+     #  7.3639505e-06,  2.5349841e-06,  6.4817141e-07,  1.0938699e-07,
+     #  1.0789142e-08,  4.0267815e-10,  9.9738412e-13,  9.3965182e-17/
+      data (gridv(iny, 8),iny=1,32)/
+     #  1.9588443e-02,  1.4054534e-02,  1.0140342e-02,  7.3511754e-03,
+     #  5.3617513e-03,  3.9285711e-03,  2.8927241e-03,  2.1409034e-03,
+     #  1.5930991e-03,  1.1934393e-03,  9.0134688e-04,  6.8734040e-04,
+     #  5.3043871e-04,  4.1513075e-04,  3.2970284e-04,  2.6531895e-04,
+     #  2.1516289e-04,  1.7405250e-04,  1.3830950e-04,  1.0585100e-04,
+     #  7.6238215e-05,  5.0368640e-05,  2.9681423e-05,  1.5103253e-05,
+     #  6.3960840e-06,  2.1473293e-06,  5.2444493e-07,  7.7848289e-08,
+     #  6.7448840e-09,  2.0379331e-10,  2.3209743e-13,  4.3818080e-17/
+      data firsttime/.true./
+      save
+      if(firsttime)then
+        firsttime=.false.
+        nent(1)=nny
+        nent(2)=nnz
+        do iny=1,nny
+          ent(iny)=(yv(iny))
+        enddo
+        do inz=1,nnz
+          ent(nny+inz)=(zv(inz))
+        enddo
+      endif
+      arg(1)=(y)
+      arg(2)=(z)
+      tmp=dfint(narg,arg,nent,ent,gridv)
+      ewa_wm_p_0=tmp
+      return
+      end
+      function ewa_z_p(x,Q2,pol)
+      implicit none
+      real*8 ewa_z_p
+      real*8 x,Q2
+      integer pol
+      real*8 tmp
+      real*8 y,z
+      real*8 ylow,yupp,zlow,zupp
+      parameter (ylow=-16.1181,yupp=-1e-07)
+      parameter (zlow=4.38683,zupp=13.8155)
+      real*8 ewa_z_p_1p
+      real*8 ewa_z_p_1m
+      real*8 ewa_z_p_0
+      y=log(x)
+      z=0.5*log(Q2)
+      if(pol.eq.1)then
+        tmp=ewa_z_p_1p(y,z)
+      else if(pol.eq.-1)then
+        tmp=ewa_z_p_1m(y,z)
+      else
+        tmp=ewa_z_p_0(y,z)
+      endif
+      ewa_z_p=tmp/x*y*y
+      end
+      function ewa_z_p_1p(y,z)
+      implicit none
+      real*8 ewa_z_p_1p,y,z
+      integer narg,nny,nnz
+      parameter (narg=2)
+      parameter (nny=32)
+      parameter (nnz=8)
+      integer iny,inz,nent(narg)
+      real*8 tmp,dfint,ymap,zmap
+      real*8 arg(narg),ent(nny+nnz)
+      real*8 yv(nny),zv(nnz),gridv(nny,nnz)
+      logical firsttime
+      external dfint,ymap,zmap
+      data yv/
+     # -1.6118096e+01, -1.5094991e+01, -1.4105430e+01, -1.3149414e+01,
+     # -1.2226943e+01, -1.1338015e+01, -1.0482632e+01, -9.6607941e+00,
+     # -8.8725001e+00, -8.1177506e+00, -7.3965455e+00, -6.7088848e+00,
+     # -6.0547686e+00, -5.4341967e+00, -4.8471693e+00, -4.2936863e+00,
+     # -3.7737478e+00, -3.2873536e+00, -2.8345039e+00, -2.4151986e+00,
+     # -2.0294377e+00, -1.6772213e+00, -1.3585493e+00, -1.0734217e+00,
+     # -8.2183848e-01, -6.0379973e-01, -4.1930540e-01, -2.6835549e-01,
+     # -1.5095001e-01, -6.7088947e-02, -1.6772312e-02, -1.0000000e-07/
+      data zv/
+     #  4.3868276e+00,  5.7337823e+00,  7.0807370e+00,  8.4276917e+00,
+     #  9.7746464e+00,  1.1121601e+01,  1.2468556e+01,  1.3815511e+01/
+      data (gridv(iny, 1),iny=1,32)/
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00/
+      data (gridv(iny, 2),iny=1,32)/
+     #  3.7303045e-03,  2.9227921e-03,  2.3141955e-03,  1.8493581e-03,
+     #  1.4933843e-03,  1.2162782e-03,  9.9810630e-04,  8.2501787e-04,
+     #  6.8659282e-04,  5.7544120e-04,  4.8609674e-04,  4.1418784e-04,
+     #  3.5626496e-04,  3.0973606e-04,  2.7244727e-04,  2.4257706e-04,
+     #  2.1854792e-04,  1.9876081e-04,  1.8151373e-04,  1.6480475e-04,
+     #  1.4652391e-04,  1.2495404e-04,  9.9767621e-05,  7.2413235e-05,
+     #  4.6189116e-05,  2.4986033e-05,  1.0981613e-05,  3.6318686e-06,
+     #  8.2837802e-07,  1.0232602e-07,  4.9850666e-09,  1.2419766e-08/
+      data (gridv(iny, 3),iny=1,32)/
+     #  1.0860874e-02,  8.3686491e-03,  6.5072730e-03,  5.1000356e-03,
+     #  4.0344753e-03,  3.2163074e-03,  2.5819036e-03,  2.0867745e-03,
+     #  1.6975997e-03,  1.3904974e-03,  1.1478197e-03,  9.5566048e-04,
+     #  8.0324203e-04,  6.8259459e-04,  5.8722232e-04,  5.1178459e-04,
+     #  4.5182361e-04,  4.0308935e-04,  3.6137523e-04,  3.2211460e-04,
+     #  2.8093250e-04,  2.3463440e-04,  1.8303741e-04,  1.2940501e-04,
+     #  8.0128195e-05,  4.1902153e-05,  1.7698447e-05,  5.5786779e-06,
+     #  1.1982518e-06,  1.3567906e-07,  5.1862975e-09,  1.4921624e-08/
+      data (gridv(iny, 4),iny=1,32)/
+     #  2.1582709e-02,  1.6384815e-02,  1.2542453e-02,  9.6686760e-03,
+     #  7.5174482e-03,  5.8869725e-03,  4.6402349e-03,  3.6814614e-03,
+     #  2.9392736e-03,  2.3625196e-03,  1.9136197e-03,  1.5633551e-03,
+     #  1.2894368e-03,  1.0755569e-03,  9.0867962e-04,  7.7828223e-04,
+     #  6.7581905e-04,  5.9350250e-04,  5.2402918e-04,  4.5996420e-04,
+     #  3.9470912e-04,  3.2387169e-04,  2.4768512e-04,  1.7121604e-04,
+     #  1.0335642e-04,  5.2501923e-05,  2.1430744e-05,  6.4822592e-06,
+     #  1.3222258e-06,  1.3895999e-07,  4.1233808e-09,  1.3983723e-08/
+      data (gridv(iny, 5),iny=1,32)/
+     #  3.6225205e-02,  2.7139223e-02,  2.0491247e-02,  1.5570792e-02,
+     #  1.1927494e-02,  9.1990020e-03,  7.1389947e-03,  5.5752406e-03,
+     #  4.3810786e-03,  3.4655687e-03,  2.7625306e-03,  2.2211128e-03,
+     #  1.8030325e-03,  1.4806433e-03,  1.2320561e-03,  1.0399620e-03,
+     #  8.9058578e-04,  7.7179147e-04,  6.7267988e-04,  5.8270496e-04,
+     #  4.9308080e-04,  3.9840200e-04,  2.9945131e-04,  2.0297089e-04,
+     #  1.1982930e-04,  5.9340172e-05,  2.3506889e-05,  6.8579981e-06,
+     #  1.3368364e-06,  1.3162902e-07,  2.9711215e-09,  1.2215759e-08/
+      data (gridv(iny, 6),iny=1,32)/
+     #  5.5154342e-02,  4.0827860e-02,  3.0449006e-02,  2.2843096e-02,
+     #  1.7269263e-02,  1.3140804e-02,  1.0059401e-02,  7.7483419e-03,
+     #  6.0046554e-03,  4.6839966e-03,  3.6820525e-03,  2.9195345e-03,
+     #  2.3374414e-03,  1.8936267e-03,  1.5550949e-03,  1.2961420e-03,
+     #  1.0966733e-03,  9.3947224e-04,  8.0960033e-04,  6.9320130e-04,
+     #  5.7933366e-04,  4.6170562e-04,  3.4169800e-04,  2.2756193e-04,
+     #  1.3169713e-04,  6.3745564e-05,  2.4577867e-05,  6.9293186e-06,
+     #  1.2930488e-06,  1.1992460e-07,  2.0948206e-09,  1.0413998e-08/
+      data (gridv(iny, 7),iny=1,32)/
+     #  7.9641332e-02,  5.8283390e-02,  4.2961224e-02,  3.1840880e-02,
+     #  2.3773251e-02,  1.7861019e-02,  1.3496990e-02,  1.0261222e-02,
+     #  7.8482212e-03,  6.0415388e-03,  4.6867971e-03,  3.6676358e-03,
+     #  2.8981380e-03,  2.3178462e-03,  1.8798546e-03,  1.5481330e-03,
+     #  1.2949777e-03,  1.0972312e-03,  9.3539626e-04,  7.9205439e-04,
+     #  6.5410473e-04,  5.1445082e-04,  3.7509234e-04,  2.4557633e-04,
+     #  1.3941666e-04,  6.6000870e-05,  2.4739779e-05,  6.6282306e-06,
+     #  1.1538853e-06,  9.8657074e-08,  1.5984130e-09,  8.2433190e-09/
+      data (gridv(iny, 8),iny=1,32)/
+     #  1.1184163e-01,  8.0917333e-02,  5.8950236e-02,  4.3164542e-02,
+     #  3.1828151e-02,  2.3609789e-02,  1.7611453e-02,  1.3215123e-02,
+     #  9.9748771e-03,  7.5776696e-03,  5.8010863e-03,  4.4800093e-03,
+     #  3.4937519e-03,  2.7582893e-03,  2.2091448e-03,  1.7974773e-03,
+     #  1.4863368e-03,  1.2455380e-03,  1.0503765e-03,  8.7953603e-04,
+     #  7.1769738e-04,  5.5700142e-04,  4.0005545e-04,  2.5742663e-04,
+     #  1.4335613e-04,  6.6352481e-05,  2.4061846e-05,  5.8581041e-06,
+     #  8.7686172e-07,  6.1694293e-08,  7.6078967e-10,  4.5762876e-09/
+      data firsttime/.true./
+      save
+      if(firsttime)then
+        firsttime=.false.
+        nent(1)=nny
+        nent(2)=nnz
+        do iny=1,nny
+          ent(iny)=(yv(iny))
+        enddo
+        do inz=1,nnz
+          ent(nny+inz)=(zv(inz))
+        enddo
+      endif
+      arg(1)=(y)
+      arg(2)=(z)
+      tmp=dfint(narg,arg,nent,ent,gridv)
+      ewa_z_p_1p=tmp
+      return
+      end
+      function ewa_z_p_1m(y,z)
+      implicit none
+      real*8 ewa_z_p_1m,y,z
+      integer narg,nny,nnz
+      parameter (narg=2)
+      parameter (nny=32)
+      parameter (nnz=8)
+      integer iny,inz,nent(narg)
+      real*8 tmp,dfint,ymap,zmap
+      real*8 arg(narg),ent(nny+nnz)
+      real*8 yv(nny),zv(nnz),gridv(nny,nnz)
+      logical firsttime
+      external dfint,ymap,zmap
+      data yv/
+     # -1.6118096e+01, -1.5094991e+01, -1.4105430e+01, -1.3149414e+01,
+     # -1.2226943e+01, -1.1338015e+01, -1.0482632e+01, -9.6607941e+00,
+     # -8.8725001e+00, -8.1177506e+00, -7.3965455e+00, -6.7088848e+00,
+     # -6.0547686e+00, -5.4341967e+00, -4.8471693e+00, -4.2936863e+00,
+     # -3.7737478e+00, -3.2873536e+00, -2.8345039e+00, -2.4151986e+00,
+     # -2.0294377e+00, -1.6772213e+00, -1.3585493e+00, -1.0734217e+00,
+     # -8.2183848e-01, -6.0379973e-01, -4.1930540e-01, -2.6835549e-01,
+     # -1.5095001e-01, -6.7088947e-02, -1.6772312e-02, -1.0000000e-07/
+      data zv/
+     #  4.3868276e+00,  5.7337823e+00,  7.0807370e+00,  8.4276917e+00,
+     #  9.7746464e+00,  1.1121601e+01,  1.2468556e+01,  1.3815511e+01/
+      data (gridv(iny, 1),iny=1,32)/
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00/
+      data (gridv(iny, 2),iny=1,32)/
+     #  2.4722465e-03,  1.9238378e-03,  1.5105106e-03,  1.1960749e-03,
+     #  9.5450636e-04,  7.6644980e-04,  6.1891326e-04,  5.0238454e-04,
+     #  4.0993066e-04,  3.3642652e-04,  2.7788916e-04,  2.3119738e-04,
+     #  1.9389505e-04,  1.6397539e-04,  1.3974473e-04,  1.1973970e-04,
+     #  1.0268876e-04,  8.7531986e-05,  7.3507411e-05,  6.0200203e-05,
+     #  4.7561978e-05,  3.5825675e-05,  2.5363959e-05,  1.6554077e-05,
+     #  9.6919422e-06,  4.9066425e-06,  2.0356874e-06,  6.4136541e-07,
+     #  1.3870147e-07,  1.5934953e-08,  3.8134239e-10,  1.5583767e-09/
+      data (gridv(iny, 3),iny=1,32)/
+     #  7.0798337e-03,  5.4081649e-03,  4.1625029e-03,  3.2271346e-03,
+     #  2.5192925e-03,  1.9775657e-03,  1.5603254e-03,  1.2371174e-03,
+     #  9.8575176e-04,  7.8989545e-04,  6.3703740e-04,  5.1756292e-04,
+     #  4.2407586e-04,  3.5069124e-04,  2.9260099e-04,  2.4581210e-04,
+     #  2.0699146e-04,  1.7346123e-04,  1.4331114e-04,  1.1548538e-04,
+     #  8.9728440e-05,  6.6379120e-05,  4.6055769e-05,  2.9366897e-05,
+     #  1.6731521e-05,  8.2008648e-06,  3.2726090e-06,  9.8342622e-07,
+     #  2.0024596e-07,  2.1085367e-08,  3.1115922e-10,  1.8414779e-09/
+      data (gridv(iny, 4),iny=1,32)/
+     #  1.3868555e-02,  1.0425215e-02,  7.8892735e-03,  6.0089146e-03,
+     #  4.6057219e-03,  3.5480566e-03,  2.7464613e-03,  2.1359653e-03,
+     #  1.6689984e-03,  1.3114613e-03,  1.0371878e-03,  8.2648300e-04,
+     #  6.6446546e-04,  5.3953412e-04,  4.4244464e-04,  3.6574043e-04,
+     #  3.0339116e-04,  2.5068534e-04,  2.0431482e-04,  1.6241791e-04,
+     #  1.2441356e-04,  9.0628951e-05,  6.1799132e-05,  3.8623433e-05,
+     #  2.1495349e-05,  1.0246602e-05,  3.9545153e-06,  1.1410065e-06,
+     #  2.2060684e-07,  2.1556121e-08,  1.5852146e-10,  1.6987942e-09/
+      data (gridv(iny, 5),iny=1,32)/
+     #  2.2986936e-02,  1.7037927e-02,  1.2705430e-02,  9.5306823e-03,
+     #  7.1915567e-03,  5.4521722e-03,  4.1524897e-03,  3.1768102e-03,
+     #  2.4418608e-03,  1.8873270e-03,  1.4682359e-03,  1.1510039e-03,
+     #  9.1077042e-04,  7.2826552e-04,  5.8860114e-04,  4.7999193e-04,
+     #  3.9314928e-04,  3.2098035e-04,  2.5857717e-04,  2.0315066e-04,
+     #  1.5370454e-04,  1.1046511e-04,  7.4187115e-05,  4.5556559e-05,
+     #  2.4836768e-05,  1.1553471e-05,  4.3299295e-06,  1.2055759e-06,
+     #  2.2273199e-07,  2.0385221e-08,  2.5365287e-11,  1.4614457e-09/
+      data (gridv(iny, 6),iny=1,32)/
+     #  3.4607710e-02,  2.5329231e-02,  1.8643112e-02,  1.3797244e-02,
+     #  1.0268406e-02,  7.6763540e-03,  5.7640388e-03,  4.3470540e-03,
+     #  3.2936565e-03,  2.5093254e-03,  1.9243527e-03,  1.4874209e-03,
+     #  1.1607380e-03,  9.1588915e-04,  7.3097517e-04,  5.8910398e-04,
+     #  4.7722044e-04,  3.8555267e-04,  3.0742488e-04,  2.3902283e-04,
+     #  1.7886336e-04,  1.2700213e-04,  8.4135064e-05,  5.0853622e-05,
+     #  2.7216172e-05,  1.2384960e-05,  4.5198910e-06,  1.2164900e-06,
+     #  2.1505007e-07,  1.8520572e-08, -5.3697542e-11,  1.2296759e-09/
+      data (gridv(iny, 7),iny=1,32)/
+     #  4.9445135e-02,  3.5757010e-02,  2.5994493e-02,  1.8993674e-02,
+     #  1.3952679e-02,  1.0293089e-02,  7.6258392e-03,  5.6738387e-03,
+     #  4.2408262e-03,  3.1872474e-03,  2.4113098e-03,  1.8389250e-03,
+     #  1.4163575e-03,  1.1035738e-03,  8.7029523e-04,  6.9356210e-04,
+     #  5.5597102e-04,  4.4470895e-04,  3.5113742e-04,  2.7029621e-04,
+     #  2.0013417e-04,  1.4045648e-04,  9.1828197e-05,  5.4657353e-05,
+     #  2.8729898e-05,  1.2794546e-05,  4.5409281e-06,  1.1613608e-06,
+     #  1.9109911e-07,  1.5092145e-08, -5.2811872e-11,  9.5507973e-10/
+      data (gridv(iny, 8),iny=1,32)/
+     #  6.8709761e-02,  4.9096606e-02,  3.5253307e-02,  2.5432125e-02,
+     #  1.8440272e-02,  1.3424173e-02,  9.8127739e-03,  7.2026791e-03,
+     #  5.3106334e-03,  3.9371633e-03,  2.9384004e-03,  2.2109274e-03,
+     #  1.6805250e-03,  1.2928794e-03,  1.0073789e-03,  7.9380317e-04,
+     #  6.2964750e-04,  4.9861161e-04,  3.8984865e-04,  2.9710009e-04,
+     #  2.1764892e-04,  1.5095775e-04,  9.7391148e-05,  5.7072702e-05,
+     #  2.9454329e-05,  1.2827026e-05,  4.4049387e-06,  1.0235568e-06,
+     #  1.4395353e-07,  9.1928386e-09, -5.1047909e-11,  4.9741292e-10/
+      data firsttime/.true./
+      save
+      if(firsttime)then
+        firsttime=.false.
+        nent(1)=nny
+        nent(2)=nnz
+        do iny=1,nny
+          ent(iny)=(yv(iny))
+        enddo
+        do inz=1,nnz
+          ent(nny+inz)=(zv(inz))
+        enddo
+      endif
+      arg(1)=(y)
+      arg(2)=(z)
+      tmp=dfint(narg,arg,nent,ent,gridv)
+      ewa_z_p_1m=tmp
+      return
+      end
+      function ewa_z_p_0(y,z)
+      implicit none
+      real*8 ewa_z_p_0,y,z
+      integer narg,nny,nnz
+      parameter (narg=2)
+      parameter (nny=32)
+      parameter (nnz=8)
+      integer iny,inz,nent(narg)
+      real*8 tmp,dfint,ymap,zmap
+      real*8 arg(narg),ent(nny+nnz)
+      real*8 yv(nny),zv(nnz),gridv(nny,nnz)
+      logical firsttime
+      external dfint,ymap,zmap
+      data yv/
+     # -1.6118096e+01, -1.5094991e+01, -1.4105430e+01, -1.3149414e+01,
+     # -1.2226943e+01, -1.1338015e+01, -1.0482632e+01, -9.6607941e+00,
+     # -8.8725001e+00, -8.1177506e+00, -7.3965455e+00, -6.7088848e+00,
+     # -6.0547686e+00, -5.4341967e+00, -4.8471693e+00, -4.2936863e+00,
+     # -3.7737478e+00, -3.2873536e+00, -2.8345039e+00, -2.4151986e+00,
+     # -2.0294377e+00, -1.6772213e+00, -1.3585493e+00, -1.0734217e+00,
+     # -8.2183848e-01, -6.0379973e-01, -4.1930540e-01, -2.6835549e-01,
+     # -1.5095001e-01, -6.7088947e-02, -1.6772312e-02, -1.0000000e-07/
+      data zv/
+     #  4.3868276e+00,  5.7337823e+00,  7.0807370e+00,  8.4276917e+00,
+     #  9.7746464e+00,  1.1121601e+01,  1.2468556e+01,  1.3815511e+01/
+      data (gridv(iny, 1),iny=1,32)/
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,
+     #  0.0000000e+00,  0.0000000e+00,  0.0000000e+00,  0.0000000e+00/
+      data (gridv(iny, 2),iny=1,32)/
+     #  2.2972746e-03,  1.7917254e-03,  1.4105851e-03,  1.1202991e-03,
+     #  8.9758753e-04,  7.2401692e-04,  5.8765422e-04,  4.7973220e-04,
+     #  3.9382732e-04,  3.2529169e-04,  2.7051159e-04,  2.2663375e-04,
+     #  1.9144281e-04,  1.6313368e-04,  1.4017673e-04,  1.2123763e-04,
+     #  1.0509369e-04,  9.0611789e-05,  7.6817182e-05,  6.3022584e-05,
+     #  4.9058611e-05,  3.5418441e-05,  2.3106352e-05,  1.3215971e-05,
+     #  6.4045348e-06,  2.5232862e-06,  7.6274114e-07,  1.6120969e-07,
+     #  2.0920899e-08,  1.1554945e-09,  2.1846051e-11,  3.8164747e-16/
+      data (gridv(iny, 3),iny=1,32)/
+     #  3.1443709e-03,  2.4086938e-03,  1.8599149e-03,  1.4469692e-03,
+     #  1.1345306e-03,  8.9494004e-04,  7.0998827e-04,  5.6630210e-04,
+     #  4.5410728e-04,  3.6631536e-04,  2.9748449e-04,  2.4339831e-04,
+     #  2.0085187e-04,  1.6729924e-04,  1.4065172e-04,  1.1916159e-04,
+     #  1.0130455e-04,  8.5745522e-05,  7.1384433e-05,  5.7493734e-05,
+     #  4.3890462e-05,  3.1023010e-05,  1.9770186e-05,  1.1015777e-05,
+     #  5.1834875e-06,  1.9746473e-06,  5.7369065e-07,  1.1562072e-07,
+     #  1.4116885e-08,  7.1378361e-10,  1.0286193e-11,  2.1760568e-16/
+      data (gridv(iny, 4),iny=1,32)/
+     #  4.0583782e-03,  3.0604630e-03,  2.3243967e-03,  1.7771165e-03,
+     #  1.3685115e-03,  1.0596807e-03,  8.2496027e-04,  6.4554632e-04,
+     #  5.0775614e-04,  4.0172285e-04,  3.1996314e-04,  2.5677121e-04,
+     #  2.0788121e-04,  1.6996988e-04,  1.4037525e-04,  1.1693715e-04,
+     #  9.7839570e-05,  8.1555616e-05,  6.6877295e-05,  5.3032557e-05,
+     #  3.9820348e-05,  2.7642036e-05,  1.7265954e-05,  9.4071798e-06,
+     #  4.3162443e-06,  1.5974573e-06,  4.4856401e-07,  8.6784619e-08,
+     #  1.0054849e-08,  4.7125287e-10,  5.0055779e-12,  1.3353325e-16/
+      data (gridv(iny, 5),iny=1,32)/
+     #  5.0242808e-03,  3.7368538e-03,  2.7975112e-03,  2.1069870e-03,
+     #  1.5976618e-03,  1.2178011e-03,  9.3297051e-04,  7.1829351e-04,
+     #  5.5582318e-04,  4.3260555e-04,  3.3896481e-04,  2.6762365e-04,
+     #  2.1321911e-04,  1.7163549e-04,  1.3964927e-04,  1.1469329e-04,
+     #  9.4678118e-05,  7.7901482e-05,  6.3058870e-05,  4.9337909e-05,
+     #  3.6517503e-05,  2.4952978e-05,  1.5315618e-05,  8.1824913e-06,
+     #  3.6722130e-06,  1.3251079e-06,  3.6114078e-07,  6.7413213e-08,
+     #  7.4593247e-09,  3.2711613e-10,  2.4150815e-12,  8.6646158e-17/
+      data (gridv(iny, 6),iny=1,32)/
+     #  6.0434242e-03,  4.4393078e-03,  3.2809706e-03,  2.4384204e-03,
+     #  1.8240018e-03,  1.3711741e-03,  1.0357484e-03,  7.8624158e-04,
+     #  5.9978597e-04,  4.6020082e-04,  3.5548403e-04,  2.7671910e-04,
+     #  2.1741591e-04,  1.7266372e-04,  1.3867688e-04,  1.1250016e-04,
+     #  9.1784840e-05,  7.4666871e-05,  5.9754966e-05,  4.6199819e-05,
+     #  3.3759746e-05,  2.2745966e-05,  1.3743797e-05,  7.2147556e-06,
+     #  3.1743092e-06,  1.1196153e-06,  2.9689509e-07,  5.3518235e-08,
+     #  5.6721821e-09,  2.3409415e-10,  1.1787607e-12,  5.8729404e-17/
+      data (gridv(iny, 7),iny=1,32)/
+     #  7.1959867e-03,  5.2238170e-03,  3.8140515e-03,  2.7991324e-03,
+     #  2.0670251e-03,  1.5335081e-03,  1.1431466e-03,  8.5620909e-04,
+     #  6.4440376e-04,  4.8779317e-04,  3.7174709e-04,  2.8552072e-04,
+     #  2.2138742e-04,  1.7357728e-04,  1.3770573e-04,  1.1041226e-04,
+     #  8.9081969e-05,  7.1687010e-05,  5.6748485e-05,  4.3377928e-05,
+     #  3.1310191e-05,  2.0811054e-05,  1.2385981e-05,  6.3920772e-06,
+     #  2.7583999e-06,  9.5083600e-07,  2.4412955e-07,  4.1454456e-08,
+     #  4.1336315e-09,  1.5732057e-10,  7.2362906e-13,  3.8540532e-17/
+      data (gridv(iny, 8),iny=1,32)/
+     #  8.5794573e-03,  6.1549995e-03,  4.4395854e-03,  3.2173803e-03,
+     #  2.3455029e-03,  1.7173947e-03,  1.2632918e-03,  9.3357781e-04,
+     #  6.9319431e-04,  5.1766046e-04,  3.8920495e-04,  2.9492776e-04,
+     #  2.2566512e-04,  1.7466681e-04,  1.3686160e-04,  1.0845033e-04,
+     #  8.6522845e-05,  6.8873718e-05,  5.3928273e-05,  4.0752757e-05,
+     #  2.9053811e-05,  1.9048548e-05,  1.1165828e-05,  5.6636607e-06,
+     #  2.3959922e-06,  8.0553030e-07,  1.9753531e-07,  2.9524070e-08,
+     #  2.5899490e-09,  8.0207162e-11,  2.5018623e-13,  1.8180313e-17/
+      data firsttime/.true./
+      save
+      if(firsttime)then
+        firsttime=.false.
+        nent(1)=nny
+        nent(2)=nnz
+        do iny=1,nny
+          ent(iny)=(yv(iny))
+        enddo
+        do inz=1,nnz
+          ent(nny+inz)=(zv(inz))
+        enddo
+      endif
+      arg(1)=(y)
+      arg(2)=(z)
+      tmp=dfint(narg,arg,nent,ent,gridv)
+      ewa_z_p_0=tmp
+      return
+      end
 
=== added file 'Template/LO/Source/PDF/ElectroweakFlux.inc'
--- Template/LO/Source/PDF/ElectroweakFlux.inc	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/PDF/ElectroweakFlux.inc	2020-08-13 08:01:43 +0000
@@ -0,0 +1,38 @@
+	double precision pi, pi2, sqrt2, eps
+	double precision ewa_mw2, ewa_mz2, ewa_gf
+	double precision ewa_qu,ewa_qd,ewa_qv,ewa_ql,ewa_t3u,ewa_t3d
+	double precision ewa_cw2,ewa_sw2,ewa_gW2
+	double precision ewa_gV,ewa_gVu,ewa_gVd,ewa_gVv,ewa_gVl
+	double precision ewa_gA,ewa_gAu,ewa_gAd,ewa_gAv,ewa_gAl
+c Math
+	parameter (pi=3.141592653589793d0)
+	parameter (pi2=9.869604401089358d0)
+	parameter (sqrt2=1.4142135623730951d0)
+	parameter (eps=1d-16)
+c Input EW parameters
+	parameter (ewa_mw2 = 80.379d0**2)  ! [GeV]^2 PDG 2018
+	parameter (ewa_mz2 = 91.1876d0**2) ! [GeV]^2 PDG 2018
+	parameter (ewa_gf  = 1.1663787d-5) ! [GeV]^-2 PDG 2018
+	parameter (ewa_qu  =  2d0/3d0)
+	parameter (ewa_qd  = -1d0/3d0)
+	parameter (ewa_qv  =  0d0)
+	parameter (ewa_ql  = -1d0)
+	parameter (ewa_t3u =  0.5d0)
+	parameter (ewa_t3d = -0.5d0)
+c Derived EW parameters
+  	parameter (ewa_cw2 = ewa_mw2/ewa_mz2)
+	parameter (ewa_sw2 = 1d0 - ewa_mw2/ewa_mz2)
+	parameter (ewa_gW2 = sqrt2*4d0*ewa_gf*ewa_mw2)
+c Weak Vector Couplings
+	parameter (ewa_gV  = 1d0)
+	parameter (ewa_gVu = 0.5d0*ewa_t3u - ewa_qu*ewa_sw2)
+	parameter (ewa_gVd = 0.5d0*ewa_t3d - ewa_qd*ewa_sw2)
+	parameter (ewa_gVv = 0.5d0*ewa_t3u)
+	parameter (ewa_gVl = 0.5d0*ewa_t3d - ewa_ql*ewa_sw2)
+c Weak Axialvector Couplings
+	parameter (ewa_gA  = -1d0)
+	parameter (ewa_gAu = -0.5d0*ewa_t3u)
+	parameter (ewa_gAd = -0.5d0*ewa_t3d)
+	parameter (ewa_gAv = -0.5d0*ewa_t3u)
+	parameter (ewa_gAl = -0.5d0*ewa_t3d)
+c
 \ No newline at end of file
 
=== modified file 'Template/LO/Source/PDF/PhotonFlux.f'
--- Template/LO/Source/PDF/PhotonFlux.f	2020-03-27 13:59:33 +0000
+++ Template/LO/Source/PDF/PhotonFlux.f	2020-08-13 08:01:43 +0000
@@ -1,30 +1,35 @@
 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/*  Equivalent photon approximation structure function.    * */
+c/*   V.M.Budnev et al., Phys.Rep. 15C (1975) 181           * */
+c/*                                                         * */
+c/*     Improved Weizsaecker-Williams formula               * */
+c/*   S. Frixione et al.   hep-ph/9310350                   * */
+c/*                                                         * */
 c/* ********************************************************* */
 c   provided by Tomasz Pierzchala - UCL
 
-      real*8 function epa_electron(x,q2max)
-      integer i
+      real*8 function epa_lepton(x,q2max, mode)
+      implicit none
+      integer i, mode, imode
+c     mode is +3/-3  for electron and +4/-4 for muon      
       real*8 x,phi_f
-      real*8 xin
+      real*8 xin(3:4)
       real*8 alpha
       real*8 f, q2min,q2max
       real*8 PI
       data PI/3.14159265358979323846/
 
-      data xin/0.511d-3/ !electron mass in GeV
+      data xin/0.511d-3, 0.105658d0/ !electron mass in GeV
 
       alpha = .0072992701
-
+      imode = abs(mode)
+      
 C     // x = omega/E = (E-E')/E
       if (x.lt.1) then
-         q2min= xin*xin*x*x/(1-x)
+         q2min= xin(imode)*xin(imode)*x*x/(1-x)
          if(q2min.lt.q2max) then 
-            f = alpha/2d0/PI*
-     &           (2d0*xin*xin*x*(-1/q2min+1/q2max)+
+             f = alpha/2d0/PI*
+     &           (2d0*xin(imode)*xin(imode)*x*(-1/q2min+1/q2max)+
      &           (2-2d0*x+x*x)/x*dlog(q2max/q2min))
             
          else
@@ -35,7 +40,7 @@
       endif
 c      write (*,*) x,dsqrt(q2min),dsqrt(q2max),f
       if (f .lt. 0) f = 0
-      epa_electron= f
+      epa_lepton = f
 
       end
 
 
=== modified file 'Template/LO/Source/PDF/makefile'
--- Template/LO/Source/PDF/makefile	2020-04-08 19:22:14 +0000
+++ Template/LO/Source/PDF/makefile	2020-08-13 08:01:43 +0000
@@ -14,15 +14,16 @@
   ifeq ($(lhapdfversion),5)
     $(error Bad lhadpfversion version 6 is now required)
   else
+    PDF         = pdfwrap_lhapdf.o pdf_lhapdf6.o pdg2pdf_lhapdf6.o opendata.o PhotonFlux.o ElectroweakFlux.o
     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 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 ElectroweakFlux.o
     endif    
   endif
 else
-  PDF         = Ctq6Pdf.o pdfwrap.o opendata.o pdf.o PhotonFlux.o pdg2pdf.o NNPDFDriver.o
+  PDF         = Ctq6Pdf.o pdfwrap.o opendata.o pdf.o PhotonFlux.o pdg2pdf.o NNPDFDriver.o ElectroweakFlux.o
 endif
 
 all: $(LIBDIR)$(LIBRARY)
 
=== modified file 'Template/LO/Source/PDF/pdf.f'
--- Template/LO/Source/PDF/pdf.f	2018-02-01 23:17:18 +0000
+++ Template/LO/Source/PDF/pdf.f	2020-08-13 08:01:43 +0000
@@ -29,7 +29,7 @@
       double precision u_val,d_val,u_sea,d_sea,s_sea,c_sea,b_sea,gluon
       double precision Ctq3df,Ctq4Fn,Ctq5Pdf,Ctq6Pdf,Ctq5L
       double precision q2max
-      double precision epa_electron,epa_proton
+      double precision epa_lepton,epa_proton
       include 'pdf.inc'
 
       integer mode,Iprtn,Irt
@@ -297,8 +297,8 @@
 c  a "diffractive" photon
 c      
       q2max=xmu*xmu
-      if(ih .eq. 3) then  !from the electron
-          fx(7)=epa_electron(x,q2max)
+      if(abs(ih) .eq. 3.or.abs(ih) .eq. 4) then  !from the electron
+          fx(7)=epa_lepton(x,q2max, ih)
       elseif(ih .eq. 2) then  !from a proton without breaking
           fx(7)=epa_proton(x,q2max)
       endif      
 
=== modified file 'Template/LO/Source/PDF/pdg2pdf.f'
--- Template/LO/Source/PDF/pdg2pdf.f	2019-12-04 21:34:25 +0000
+++ Template/LO/Source/PDF/pdg2pdf.f	2020-08-13 08:01:43 +0000
@@ -24,7 +24,7 @@
       integer mode,Irt,i,j
       double precision xlast(2),xmulast(2),pdflast(-7:7,2),q2max
       character*7 pdlabellast(2)
-      double precision epa_electron,epa_proton
+      double precision epa_lepton,epa_proton
       integer ipart,ireuse,iporg,ihlast(2)
       save xlast,xmulast,pdflast,pdlabellast,ihlast
       data xlast/2*-99d9/
@@ -32,7 +32,29 @@
       data pdflast/30*-99d9/
       data pdlabellast/2*'abcdefg'/
       data ihlast/2*-99/
-
+c
+c     Special for EWA (need to know current helicity)
+c
+      double precision ewa_wX_ByPID, ewa_zX_ByPID
+      external ewa_wX_ByPID, ewa_zX_ByPID
+      integer hel, ppid
+      double precision hel_jacobian
+      INTEGER HEL_PICKED
+      COMMON/HEL_PICKED/HEL_PICKED,hel_jacobian
+      integer get_nhel
+      external get_nhel
+      REAL*8 POL(2)
+      common/to_polarization/ POL
+c
+c     Collider
+c
+      integer                                        lpp(2)
+      double precision    ebeam(2), xbk(2),q2fact(2)
+      common/to_collider/ ebeam   , xbk   ,q2fact,   lpp
+c
+c  
+c
+      double precision factor
       if (ih.eq.9) then
          pdg2pdf = 1d0
          return
@@ -66,12 +88,14 @@
          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
+         if(iabs(ih).lt.2)then
+            write(*,*) 'PDF not supported for pdg ',ipdg
+            write(*,*) 'For lepton colliders, please set the lpp* '//
+     $           '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
+         endif
       endif
 
       iporg=ipart
@@ -85,7 +109,7 @@
       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))
@@ -129,15 +153,85 @@
       pdlabellast(ireuse)=pdlabel
       ihlast(ireuse)=ih
 
-      if(iabs(ipart).eq.7.and.ih.gt.1) then
+c
+c     Special mode for Improved Weizsacker-Williams,
+c      Equivalent photon, and Effective W/Z Approximations
+c
+      if (ih.gt.1.and.ih.ne.9) then
+c        Effective photon approximation
          q2max=xmu*xmu
-         if(ih.eq.3) then       !from the electron
-            pdg2pdf=epa_electron(x,q2max)
-         elseif(ih .eq. 2) then !from a proton without breaking
-            pdg2pdf=epa_proton(x,q2max,beamid)
-         endif 
-         pdflast(iporg,ireuse)=pdg2pdf
-         return
+         if (iabs(ipart).eq.7) then
+            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
+         elseif(iabs(ipart).eq.23) then ! z boson
+            factor = 3d0        ! 3 to cancel the average over W helicities´
+            if (ABS(ih).eq.3) then
+               ppid = 11
+            elseif(ABS(ih).eq.4) then
+               ppid = 13
+            elseif(ABS(ih).eq.2) then
+               ppid = 2212
+            else
+               write(*,*) "type of beam not yet supported for EWA"
+               stop 4
+            endif
+            hel = GET_NHEL(HEL_PICKED, beamid)
+            pdg2pdf = factor * ewa_zX_ByPID(x,q2max,hel,ppid)
+            if (ABS(POL(beamid)).ne.1d0) then
+               write(*,*) 'beam polarization not supported for EWA with Z'
+               stop 33
+            endif 
+
+            return
+         elseif(ipart.eq.24) then ! w+ boson
+            factor = 3d0        ! 3 to cancel the average over W helicities´
+            if (ih.eq.-3) then    ! from e+ 
+               ppid = -11
+            elseif(ih.eq.-4) then ! from mu+
+               ppid = -13
+            elseif(ih.eq.2) then
+               ppid = 2212
+            else
+               write(*,*) "type of beam not yet supported for EWA"
+               stop 4
+            endif
+            hel = GET_NHEL(HEL_PICKED, beamid)
+            pdg2pdf = factor*ewa_wX_ByPID(x,q2max,hel,ppid) 
+            if(pol(beamid).lt.0d0)then
+               pdg2pdf = pdg2pdf * (2-ABS(POL(beamid)))/2d0 ! couple only to right handed componenta
+            else
+               pdg2pdf = pdg2pdf * ABS(POL(beamid))/2d0 ! couple only to right handed componenta
+            endif
+            return
+         elseif(ipart.eq.-24) then ! w- boson
+            factor = 3d0        ! 3 to cancel the average over W helicities´
+            if (ih.eq.3) then      ! from e-
+               ppid = 11
+            elseif(ih.eq.4) then   ! from mu-
+               ppid = 13
+            elseif(ih.eq.2) then
+               ppid = -2212
+            else
+               write(*,*) "type of beam not yet supported for EWA"
+               stop 4
+            endif
+
+            hel = GET_NHEL(HEL_PICKED, beamid)
+            pdg2pdf = factor*ewa_wX_ByPID(x,q2max,hel,ppid)
+            if(pol(beamid)*lpp(beamid).lt.0d0)then
+               pdg2pdf = pdg2pdf * ABS(POL(beamid))/2d0 ! couple only to left handed componenta
+            else
+               pdg2pdf = pdg2pdf * (2-ABS(POL(beamid)))/2d0 ! couple only to left handed componenta
+            endif
+
+
+            return
+         endif
       endif
       
       if (pdlabel(1:5) .eq. 'cteq6') then
 
=== removed file 'Template/LO/Source/PDF/pdg2pdf_lhapdf.f'
=== modified file 'Template/LO/Source/PDF/pdg2pdf_lhapdf6.f'
--- Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	2020-03-06 20:47:59 +0000
+++ Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	2020-08-13 08:01:43 +0000
@@ -24,7 +24,7 @@
       integer i,j,ihlast(20),ipart,iporg,ireuse,imemlast(20),iset,imem
      &     ,i_replace,ii,ipartlast(20)
       double precision xlast(20),xmulast(20),pdflast(-7:7,20)
-      double precision epa_proton, epa_electron
+      double precision epa_proton, epa_lepton
       save ihlast,xlast,xmulast,pdflast,imemlast,ipartlast
       data ihlast/20*-99/
       data ipartlast/20*-99/
@@ -33,6 +33,37 @@
       data pdflast/300*-99d9/
       data imemlast/20*-99/
       data i_replace/20/
+c
+c     Special for EWA (need to know current helicity)
+c
+      double precision ewa_wX_ByPID, ewa_zX_ByPID
+      external ewa_wX_ByPID, ewa_zX_ByPID
+      integer hel, ppid
+      double precision hel_jacobian
+      INTEGER HEL_PICKED
+      COMMON/HEL_PICKED/HEL_PICKED,hel_jacobian
+      integer get_nhel
+      external get_nhel
+      REAL*8 POL(2)
+      common/to_polarization/ POL
+      double precision q2max
+      double precision epa_electron
+c
+c     Collider
+c
+      integer                                        lpp(2)
+      double precision    ebeam(2), xbk(2),q2fact(2)
+      common/to_collider/ ebeam   , xbk   ,q2fact,   lpp
+c
+c  
+c
+      double precision factor
+
+      if (ih.eq.9) then
+         pdg2pdf = 1d0
+         return
+      endif
+
 
       nb_hadron = (nb_proton(beamid)+nb_neutron(beamid))
 c     Make sure we have a reasonable Bjorken x. Note that even though
@@ -62,12 +93,14 @@
          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
+         if(iabs(ih).lt.2)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
+         endif
       endif
 
       iporg=ipart
@@ -139,10 +172,77 @@
             pdg2pdf = get_ion_pdf(pdflast(-7, i_replace), ipart, nb_proton(beamid), nb_neutron(beamid))
          endif
          pdg2pdf=pdg2pdf/x
-      else if(ih.eq.3) then       !from the electron
-            pdg2pdf=epa_electron(x,xmu*xmu)
-      else if(ih.eq.2) then ! photon from a proton without breaking
-          pdg2pdf = epa_proton(x,xmu*xmu)
+      else if(ih.gt.1.and.ih.ne.9) then ! photon from a proton without breaking
+c        Effective photon approximation
+         q2max=xmu*xmu
+         if (iabs(ipart).eq.7) then
+            if(abs(ih).eq.3.or.abs(ih).eq.4) then       !from the electron
+               pdg2pdf=epa_lepton(x,xmu*xmu, ih)
+            elseif(ih .eq. 2) then !from a proton without breaking
+               pdg2pdf=epa_proton(x,q2max,beamid)
+            endif 
+            pdflast(iporg,ireuse)=pdg2pdf
+            return
+         elseif(iabs(ipart).eq.23) then ! z boson
+            if (ih.eq.3) then
+               ppid = 11
+            elseif(ih.eq.4) then
+               ppid = 13
+            elseif(ih.eq.2) then
+               ppid = 2212
+            else
+               write(*,*) "type of beam not yet supported for EWA"
+               stop 4
+            endif
+            hel = GET_NHEL(HEL_PICKED, beamid)
+            pdg2pdf = ewa_zX_ByPID(x,q2max,hel,ppid)
+            if (ABS(POL(beamid)).ne.1d0) then
+               write(*,*) 
+     &              'beam polarization not supported for EWA with Z'
+               stop 33
+            endif 
+            return
+         elseif(ipart.eq.24) then ! w+ boson
+            factor = 3d0        ! 3 to cancel the average over W helicities´
+            if (ih.eq.-3) then    ! from e+ 
+               ppid = -11
+            elseif(ih.eq.-4) then ! from mu+
+               ppid = -13
+            elseif(ih.eq.2) then
+               ppid = 2212
+            else
+               write(*,*) "type of beam not yet supported for EWA"
+               stop 4
+            endif
+            hel = GET_NHEL(HEL_PICKED, beamid)
+            pdg2pdf = factor*ewa_wX_ByPID(x,q2max,hel,ppid) 
+            if(pol(beamid).lt.0d0)then
+               pdg2pdf = pdg2pdf * (2-ABS(POL(beamid)))/2d0 ! couple only to right handed componenta
+            else
+               pdg2pdf = pdg2pdf * ABS(POL(beamid))/2d0 ! couple only to right handed componenta
+            endif
+            return
+         elseif(ipart.eq.-24) then ! w- boson
+            factor = 3d0        ! 3 to cancel the average over W helicities´
+            if (ih.eq.3) then      ! from e-
+               ppid = 11
+            elseif(ih.eq.4) then   ! from mu-
+               ppid = 13
+            elseif(ih.eq.2) then
+               ppid = -2212
+            else
+               write(*,*) "type of beam not yet supported for EWA"
+               stop 4
+            endif
+
+            hel = GET_NHEL(HEL_PICKED, beamid)
+            pdg2pdf = factor*ewa_wX_ByPID(x,q2max,hel,ppid)
+            if(pol(beamid)*lpp(beamid).lt.0d0)then
+               pdg2pdf = pdg2pdf * ABS(POL(beamid))/2d0 ! couple only to left handed componenta
+            else
+               pdg2pdf = pdg2pdf * (2-ABS(POL(beamid)))/2d0 ! couple only to left handed componenta
+            endif
+         endif 
       else
          write (*,*) 'beam type not supported in lhadpf'
          stop 1
 
=== modified file 'Template/LO/Source/kin_functions.f'
--- Template/LO/Source/kin_functions.f	2016-02-10 14:56:17 +0000
+++ Template/LO/Source/kin_functions.f	2020-08-13 08:01:43 +0000
@@ -590,11 +590,13 @@
 C     4-Vector Dot product
 C****************************************************************************
       implicit none
-      double precision p1(0:3),p2(0:3)
+      double precision p1(0:3),p2(0:3), dot2
       dot=p1(0)*p2(0)-p1(1)*p2(1)-p1(2)*p2(2)-p1(3)*p2(3)
 
-      if(dabs(dot).lt.1d-6)then ! solve numerical problem 
-         dot=0d0
+      if(dabs(dot).lt.1d-6)then ! solve numerical problem
+         dot2 = p1(0)*p2(0)+p1(1)*p2(1)+p1(2)*p2(2)+p1(3)*p2(3)
+         dot2 = max(1e-99,dot2)
+         if (dot/dot2.lt.1e-6) dot=0d0
       endif
 
       end
 
=== modified file 'Template/LO/Source/setrun.f'
--- Template/LO/Source/setrun.f	2016-03-04 16:16:23 +0000
+++ Template/LO/Source/setrun.f	2020-08-13 08:01:43 +0000
@@ -83,12 +83,27 @@
       q2fact(1) = sf1**2      ! fact scale**2 for pdf1
       q2fact(2) = sf2**2      ! fact scale**2 for pdf2     
 
+      if(pb1.ne.0d0)then
+         if (abs(lpp(1)).eq.1.or.abs(lpp(1)).eq.2)then
+            write(*,*) 'proton/anti-proton beam polarization are not allowed'
+            stop 1
+         endif
+         pol(1)=sign(1+abs(pb1)/100d0,pb1)
+      endif
+      if(pb2.ne.0d0)then
+         if (abs(lpp(2)).eq.1.or.abs(lpp(2)).eq.2)then
+            write(*,*) 'proton/anti-proton beam polarization are not allowed'
+            stop 1
+         endif
+         pol(2)=sign(1+abs(pb2)/100d0,pb2)
+      endif
+      
       if(pb1.ne.0d0.and.lpp(1).eq.0) pol(1)=sign(1+abs(pb1)/100d0,pb1)
       if(pb2.ne.0d0.and.lpp(2).eq.0) pol(2)=sign(1+abs(pb2)/100d0,pb2)
 
       if(pb1.ne.0.or.pb2.ne.0) write(*,*) 'Setting beam polarization ',
-     $     sign((abs(pol(1))-1)*100,pol(1)),
-     $     sign((abs(pol(2))-1)*100,pol(2))
+     $     pol(1), pol(2)
+
 
 c !!! Default behavior changed (MH, Aug. 07) !!!
 c If no pdf, read the param_card and use the value from there and
 
=== modified file 'Template/LO/SubProcesses/genps.f'
--- Template/LO/SubProcesses/genps.f	2020-06-15 07:16:51 +0000
+++ Template/LO/SubProcesses/genps.f	2020-08-13 08:01:43 +0000
@@ -572,6 +572,9 @@
             if (abs(lpp(2)) .eq. 1 .or. abs(lpp(2)) .eq. 2) m2 = 0.938d0
             if (abs(lpp(1)) .eq. 3) m1 = 0.000511d0
             if (abs(lpp(2)) .eq. 3) m2 = 0.000511d0
+            if (abs(lpp(1)) .eq. 4) m1 = 0.000511d0
+            if (abs(lpp(1)) .eq. 4) m1 = 0.105658d0
+            if (abs(lpp(2)) .eq. 4) m2 = 0.105658d0
             if (mass_ion(1).ge.0d0) m1 = mass_ion(1)
             if (mass_ion(2).ge.0d0) m2 = mass_ion(2)
             if(ebeam(1).lt.m1.and.lpp(1).ne.9) ebeam(1)=m1
 
=== modified file 'Template/LO/SubProcesses/reweight.f'
--- Template/LO/SubProcesses/reweight.f	2020-05-24 20:16:56 +0000
+++ Template/LO/SubProcesses/reweight.f	2020-08-13 08:01:43 +0000
@@ -1738,9 +1738,9 @@
          if (btest(mlevel,3))
      $        write(*,*)' set fact scales for PS to ',
      $        sqrt(q2fact(1)),sqrt(q2fact(2))
-      else if (abs(lpp(1)).eq.2.or.abs(lpp(1)).eq.3) then
+      else if (abs(lpp(1)).ge.2.and.abs(lpp(1)).le.4) then
          q2fact(1)=q2bck(1)
-      else if (abs(lpp(2)).eq.2.or.abs(lpp(2)).eq.3) then
+      else if (abs(lpp(2)).ge.2.or.abs(lpp(2)).le.4) then
          q2fact(2)=q2bck(2)
       endif
 
 
=== modified file 'Template/LO/SubProcesses/unwgt.f'
--- Template/LO/SubProcesses/unwgt.f	2020-05-20 11:52:45 +0000
+++ Template/LO/SubProcesses/unwgt.f	2020-08-13 08:01:43 +0000
@@ -145,10 +145,19 @@
 C     compute the boost for the requested transformation
       implicit none
       double precision pin(0:3), pout(0:3)
-      double precision denom
+      double precision denom, control
 
       denom = pin(0)*pout(0) - pin(3)*pout(3)
       if (denom.ne.0d0) then
+         control = pin(0)*pout(0) + pin(3)*pout(3)
+         if (denom/control.gt.1d-12)then
+            get_betaz = (pin(3) * pout(0) - pout(3) * pin(0)) / denom
+c        pass in massless case
+         else if (abs(pin(0)+pin(3)).gt.abs(pin(0))) then
+            get_betaz = (pin(0)**2 - pout(0)**2)/(pin(0)**2 + pout(0)**2)
+         else
+            get_betaz = (pout(0)**2 - pin(0)**2)/(pin(0)**2 + pout(0)**2)
+         endif
          get_betaz = (pin(3) * pout(0) - pout(3) * pin(0)) / denom
       else if (pin(0).eq.pin(3)) then
          get_betaz = (pin(0)**2 - pout(0)**2)/(pin(0)**2 + pout(0)**2)
 
=== modified file 'Template/MadWeight/src/setrun.f'
--- Template/MadWeight/src/setrun.f	2018-04-27 21:01:37 +0000
+++ Template/MadWeight/src/setrun.f	2020-08-13 08:01:43 +0000
@@ -536,3 +536,11 @@
 
       return
       end
+
+      INTEGER FUNCTION GET_NHEL(HEL, IPART)
+      implicit none
+      integer hel
+      integer ipart
+      get_nhel = 99
+      return 
+      end
 
=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2020-06-21 12:17:36 +0000
+++ UpdateNotes.txt	2020-08-13 08:01:43 +0000
@@ -1,5 +1,23 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
+2.8.0:
+      OM: In madspin_card you can now replace the line "launch"
+          by "launch -n NAME", this will allow to specify the name of the
+	  directory in EVENTS where that run is stored.
+      OM: Change in the python interface of the standalone output. In top of the pdg of the particles,
+          you can now use the process_id (the one specified with @X) to distinguish process with the same particle
+	  content. This parameter can be set to -1 and the function then ignore that parameter.
+      OM: Adding new option in reweighting to allow the user to use the process_id in presence of ambiguous
+          initial/final state.
+      OM: For LO process, you can now set lpp1 and lpp2 to "4" for process with initial photon in order to get the
+          effective photon approximation. This mode behaves like the "3" one: The cut-off scale of the approximation
+	  is taken from the fixed renormalization factorisation scale: "dsqrt_q2fact1/dsqrt_q2fact2"
+      OM: Update the makefile of standalone interface to python to be able to compile in multicore.
+           (thanks Matthias Komm)
+      OM: Update of the auto-width code to support UFO form-factors
+      OM: Fixing numerical issue at very low invariant mass and another numerical issue with the boost in EPA mode.
+      
+
 2.7.3(21/06/20)
       OM: Fixing some bug for read-only LO gridpacks (wrong cross-section and shape when generating events).
           Thanks to Congqiao Li for this
 
=== modified file 'aloha/aloha_writers.py'
--- aloha/aloha_writers.py	2020-01-09 23:14:56 +0000
+++ aloha/aloha_writers.py	2020-08-13 08:01:43 +0000
@@ -1241,12 +1241,15 @@
     """ build the name for combined aloha function """
 
     def myHash(target_string):
+        suffix = ''
+        if '%(propa)s' in target_string:
+            target_string = target_string.replace('%(propa)s','')
+            suffix = '%(propa)s'
+            
         if len(target_string)<50:
-            return target_string
-        if '%(propa)s' in target_string:
-            return 'ALOHA_'+(str(hash(target_string.lower()))).replace('-','m')+'%(propa)s'
+            return '%s%s' % (target_string, suffix)
         else:
-            return 'ALOHA_'+(str(hash(target_string.lower()))).replace('-','m')
+            return 'ALOHA_%s%s' % (str(hash(target_string.lower())).replace('-','m'), suffix)
 
     if tag and any(t.startswith('P') for t in tag[:-1]):
         # propagator need to be the last entry for the tag
 
=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2020-05-07 13:47:24 +0000
+++ madgraph/core/base_objects.py	2020-08-13 08:01:43 +0000
@@ -1222,7 +1222,8 @@
 
         if name == 'particles':
             # Ensure no doublets in particle list
-            make_unique(value)
+            if value:
+                make_unique(value)
             # Reset dictionaries
             self['particle_dict'] = {}
             self['ref_dict_to0'] = {}
@@ -4053,7 +4054,7 @@
     is needed"""
 
     assert isinstance(doubletlist, list), \
-           "Argument to make_unique must be list"
+           "Argument to make_unique must be list not %s" % type(doubletlist)
     
 
     uniquelist = []
 
=== modified file 'madgraph/interface/common_run_interface.py'
--- madgraph/interface/common_run_interface.py	2020-06-17 07:23:03 +0000
+++ madgraph/interface/common_run_interface.py	2020-08-13 08:01:43 +0000
@@ -3704,11 +3704,21 @@
 
         madspin_cmd.import_command_file(path)
 
-        # create a new run_name directory for this output
-        i = 1
-        while os.path.exists(pjoin(self.me_dir,'Events', '%s_decayed_%i' % (self.run_name,i))):
-            i+=1
-        new_run = '%s_decayed_%i' % (self.run_name,i)
+
+        if not madspin_cmd.me_run_name: 
+            # create a new run_name directory for this output
+            i = 1
+            while os.path.exists(pjoin(self.me_dir,'Events', '%s_decayed_%i' % (self.run_name,i))):
+                i+=1
+            new_run = '%s_decayed_%i' % (self.run_name,i)
+        else:
+            new_run = madspin_cmd.me_run_name
+            if os.path.exists(pjoin(self.me_dir,'Events', new_run)):
+                i = 1
+                while os.path.exists(pjoin(self.me_dir,'Events', '%s_%i' % (new_run,i))):
+                    i+=1
+                new_run = '%s_%i' % (new_run,i)
+
         evt_dir = pjoin(self.me_dir, 'Events')
 
         os.mkdir(pjoin(evt_dir, new_run))
@@ -5828,6 +5838,21 @@
                  self.run_card['mass_ion1'] != self.run_card['mass_ion2']):
                 raise Exception, "Heavy ion profile for both beam are different but the symmetry used forbids it. \n Please generate your process with \"set group_subprocesses False\"."
             
+            # EWA mode should use MC over polarization and no systematics (due to pdf call inreweight)
+            if self.run_card['lpp1'] in [2,3,4,-3,-4] or self.run_card['lpp2'] in [2,3,4,-3,-4]:
+                if any(pdg in proc_charac['pdg_initial1']+ proc_charac['pdg_initial2']
+                       for pdg in [24, -24, 23]):
+                    if self.run_card['nhel'] == 0: 
+                        logger.warning("only MC over polarization supported for EWA mode")
+                        self.do_set('run_card nhel 1') 
+                    if self.run_card['use_syst']:
+                        logger.warning("evaluation of systematics not allowed in EWA mode. set use_syst to False")
+                        self.do_set('run_card use_syst False')                                    
+                
+                if self.run_card['lpp1'] in [4,-4] or self.run_card['lpp2'] in [4,-4]:
+                    if 22 in proc_charac['pdg_initial1']+ proc_charac['pdg_initial2']:
+                        raise Exception, "EPA not implemented from muon beam. (please open a ticket and we will do it)"
+                                   
             # check the status of small width status from LO
             for param in self.param_card['decay']:
                 width = param.value
 
=== modified file 'madgraph/interface/madevent_interface.py'
--- madgraph/interface/madevent_interface.py	2020-06-18 20:23:25 +0000
+++ madgraph/interface/madevent_interface.py	2020-08-13 08:01:43 +0000
@@ -5414,6 +5414,7 @@
         elif mode in [1,2]:
             exename = os.path.basename(exe)
             # For condor cluster, create the input/output files
+            misc.sprint(exename)
             if 'ajob' in exename: 
                 input_files = ['madevent','input_app.txt','symfact.dat','iproc.dat','dname.mg',
                                pjoin(self.me_dir, 'SubProcesses','randinit')]
@@ -5509,7 +5510,7 @@
                 #submitting
                 self.cluster.cluster_submit(exe, stdout=stdout, cwd=cwd, argument=argument,  
                              input_files=input_files, output_files=output_files,
-                             required_output=required_output, **opt)
+                             required_output=required_output, allow_array=True,  **opt)
             elif "refine_splitted.sh" in exename:
                 input_files = ['madevent','symfact.dat','iproc.dat', 'dname.mg',
                                pjoin(self.me_dir, 'SubProcesses','randinit')]                 
@@ -5533,7 +5534,7 @@
                 #submitting
                 self.cluster.cluster_submit(exe, stdout=stdout, cwd=cwd, argument=argument,  
                              input_files=input_files, output_files=output_files,
-                             required_output=required_output, **opt)                
+                             required_output=required_output,allow_array=True, **opt)                
 
                 
             
 
=== modified file 'madgraph/interface/madgraph_interface.py'
--- madgraph/interface/madgraph_interface.py	2020-06-17 06:15:07 +0000
+++ madgraph/interface/madgraph_interface.py	2020-08-13 08:01:43 +0000
@@ -8353,7 +8353,7 @@
             decay_dir = pjoin(path,'temp_decay')
             logger_mg.info('More info in temporary files:\n    %s/index.html' % (decay_dir))
             with misc.MuteLogger(['madgraph','ALOHA','cmdprint','madevent'], [40,40,40,40]):
-                self.exec_cmd('output %s -f' % decay_dir,child=False)
+                self.exec_cmd('output madevent %s -f' % decay_dir,child=False)
                 
                 #modify some parameter of the default run_card
                 run_card = banner_module.RunCard(pjoin(decay_dir,'Cards','run_card.dat'))
 
=== modified file 'madgraph/interface/reweight_interface.py'
--- madgraph/interface/reweight_interface.py	2020-06-16 13:10:57 +0000
+++ madgraph/interface/reweight_interface.py	2020-08-13 08:01:43 +0000
@@ -111,6 +111,7 @@
         self.rwgt_dir = None
         self.exitted = False # Flag to know if do_quit was already called.
         self.keep_ordering = False
+        self.use_eventid = False
         if event_path:
             logger.info("Extracting the banner ...")
             self.do_import(event_path, allow_madspin=allow_madspin)
@@ -356,8 +357,8 @@
             if self.has_standalone_dir:
                 self.terminate_fortran_executables()
                 self.has_standalone_dir = False
-        elif args[0] == "keep_ordering":
-            self.keep_ordering = banner.ConfigFile.format_variable(args[1], bool, "keep_ordering")
+        elif args[0] in ["keep_ordering", "use_eventid"]:
+            setattr(self, args[0], banner.ConfigFile.format_variable(args[1], bool, args[0]))
         elif args[0] == "allow_missing_finalstate":
             self.options["allow_missing_finalstate"] = banner.ConfigFile.format_variable(args[1], bool, "allow_missing_finalstate")
         elif args[0] == "process":
@@ -980,11 +981,12 @@
     def calculate_weight(self, event):
         """space defines where to find the calculator (in multicore)"""
         
-        if not hasattr(self,'pdf'):
-            lhapdf = misc.import_python_lhapdf(self.mg5cmd.options['lhapdf'])
-            self.pdf = lhapdf.mkPDF(self.banner.run_card.get_lhapdf_id())
-            
+
         if self.has_nlo and self.rwgt_mode != "LO":
+            if not hasattr(self,'pdf'):
+                lhapdf = misc.import_python_lhapdf(self.mg5cmd.options['lhapdf'])
+                self.pdf = lhapdf.mkPDF(self.banner.run_card.get_lhapdf_id())
+                
             return self.calculate_nlo_weight(event)
         
         event.parse_reweight()                    
@@ -1206,7 +1208,7 @@
         #    base = "rw_mevirt"
         #else:
         #    base = "rw_me"
-        
+
         if (not self.second_model and not self.second_process and not self.dedicated_path) or hypp_id==0:
             orig_order, Pdir, hel_dict = self.id_to_path[tag]
         else:
@@ -1269,10 +1271,16 @@
         pold = list(p)
         p = self.invert_momenta(p)
         pdg = list(orig_order[0])+list(orig_order[1])
-
+        try:
+            pid = event.ievent
+        except AttributeError:
+            pid = -1
+        if not self.use_eventid:
+            pid = -1
+        
         with misc.chdir(Pdir):
             with misc.stdchannel_redirected(sys.stdout, os.devnull):
-                me_value = module.smatrixhel(pdg,p, event.aqcd, scale2, nhel)
+                me_value = module.smatrixhel(pdg, pid, p, event.aqcd, scale2, nhel)
                                 
         # for loop we have also the stability status code
         if isinstance(me_value, tuple):
@@ -1387,10 +1395,12 @@
                     nlo_order = nlo_order.replace('noborn', 'virt')
                 commandline += "add process %s [%s] %s;" % (base,nlo_order,post)
             commandline = commandline.replace('add process', 'generate',1)
-            logger.info("RETRY with %s", commandline)
-            mgcmd.exec_cmd(commandline, precmd=True)
-            has_nlo = False
+            if commandline:
+                logger.info("RETRY with %s", commandline)
+                mgcmd.exec_cmd(commandline, precmd=True)
+                has_nlo = False
         except Exception, error:
+            misc.sprint(type(error))
             raise
         
         commandline = 'output standalone_rw %s --prefix=int' % pjoin(path_me,data['paths'][0])
@@ -1802,8 +1812,9 @@
                 data = self.id_to_path_second
 
             # get all the information
-            all_pdgs = mymod.get_pdg_order()
-            all_pdgs = [[pdg for pdg in pdgs if pdg!=0] for pdgs in  mymod.get_pdg_order()]
+            allids, all_pids = mymod.get_pdg_order()
+            all_pdgs = [[pdg for pdg in pdgs if pdg!=0] for pdgs in  allids]
+            #all_pids = [pid for (pdgs, pid) in  allids]
             all_prefix = [''.join(j).strip().lower() for j in mymod.get_prefix()]
             prefix_set = set(all_prefix)
 
@@ -1826,7 +1837,7 @@
                     misc.sprint(os.path.exists(pjoin(path_me,onedir,'SubProcesses','MadLoop5_resources', '%sHelConfigs.dat' % prefix.upper())))
                     continue
 
-            for i,pdg in enumerate(all_pdgs):
+            for i,(pdg,pid) in enumerate(zip(all_pdgs,all_pids)):
                 if self.is_decay:
                     incoming = [pdg[0]]
                     outgoing = pdg[1:]
 
=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2020-06-20 21:10:10 +0000
+++ madgraph/iolibs/export_v4.py	2020-08-13 08:01:43 +0000
@@ -166,6 +166,7 @@
         self.mgme_dir = MG5DIR
         self.dir_path = dir_path
         self.model = None
+        self.beam_polarization = True
 
         self.opt = dict(self.default_opt)
         if opt:
@@ -185,12 +186,48 @@
     def export_processes(self, matrix_elements, fortran_model):
         """Make the switch between grouped and not grouped output"""
         
+
         calls = 0
         if isinstance(matrix_elements, group_subprocs.SubProcessGroupList):
+            # check handling for the polarization
+            for m in matrix_elements:
+                for me in m.get('matrix_elements'):
+                    for p in me.get('processes'):
+                        for pid in p.get_initial_ids():
+                            spin = p.get('model').get_particle(pid).get('spin')
+                            if spin != 2:
+                                self.beam_polarization = False
+                                break
+                        else:
+                            continue
+                        break
+                    else:
+                        continue
+                    break
+                else:
+                    continue
+                break
+
             for (group_number, me_group) in enumerate(matrix_elements):
                 calls = calls + self.generate_subprocess_directory(\
                                           me_group, fortran_model, group_number)
         else:
+            # check handling for the polarization
+            self.beam_polarization = True
+            for me in matrix_elements.get_matrix_elements():
+                for p in me.get('processes'):
+                    for pid in p.get_initial_ids():
+                        spin = p.get('model').get_particle(pid).get('spin')
+                        if spin != 2:
+                            self.beam_polarization = False
+                            break
+                    else:
+                        continue
+                    break
+                else:
+                    continue
+                break                
+            
             for me_number, me in enumerate(matrix_elements.get_matrix_elements()):
                 calls = calls + self.generate_subprocess_directory(\
                                                    me, fortran_model, me_number)    
@@ -2066,17 +2103,18 @@
         
         template = """
 %(python_information)s
-  subroutine smatrixhel(pdgs, npdg, p, ALPHAS, SCALE2, nhel, ANS)
+  subroutine smatrixhel(pdgs, procid, npdg, p, ALPHAS, SCALE2, nhel, ANS)
   IMPLICIT NONE
 
 CF2PY double precision, intent(in), dimension(0:3,npdg) :: p
 CF2PY integer, intent(in), dimension(npdg) :: pdgs
+CF2PY integer, intent(in):: procid
 CF2PY integer, intent(in) :: npdg
 CF2PY double precision, intent(out) :: ANS
 CF2PY double precision, intent(in) :: ALPHAS
 CF2PY double precision, intent(in) :: SCALE2
   integer pdgs(*)
-  integer npdg, nhel
+  integer npdg, nhel, procid
   double precision p(*)
   double precision ANS, ALPHAS, PI,SCALE2
   include 'coupl.inc'
@@ -2128,12 +2166,16 @@
     end
       
 
-    subroutine get_pdg_order(PDG)
+    subroutine get_pdg_order(PDG, ALLPROC)
   IMPLICIT NONE
 CF2PY INTEGER, intent(out) :: PDG(%(nb_me)i,%(maxpart)i)  
+CF2PY INTEGER, intent(out) :: ALLPROC(%(nb_me)i)
   INTEGER PDG(%(nb_me)i,%(maxpart)i), PDGS(%(nb_me)i,%(maxpart)i)
+  INTEGER ALLPROC(%(nb_me)i),PIDs(%(nb_me)i)
   DATA PDGS/ %(pdgs)s /
+  DATA PIDS/ %(pids)s /
   PDG = PDGS
+  ALLPROC = PIDS
   RETURN
   END 
 
@@ -2151,31 +2193,33 @@
          
         allids = self.prefix_info.keys()
         allprefix = [self.prefix_info[key][0] for key in allids]
-        min_nexternal = min([len(ids) for ids in allids])
-        max_nexternal = max([len(ids) for ids in allids])
+        min_nexternal = min([len(ids[0]) for ids in allids])
+        max_nexternal = max([len(ids[0]) for ids in allids])
 
         info = []
-        for key, (prefix, tag) in self.prefix_info.items():
-            info.append('#PY %s : %s # %s' % (tag, key, prefix))
+        for (key, pid), (prefix, tag) in self.prefix_info.items():
+            info.append('#PY %s : %s # %s %s' % (tag, key, prefix, pid))
             
 
         text = []
         for n_ext in range(min_nexternal, max_nexternal+1):
-            current = [ids for ids in allids if len(ids)==n_ext]
-            if not current:
+            current_id = [ids[0] for ids in allids if len(ids[0])==n_ext]
+            current_pid = [ids[1] for ids in allids if len(ids[0])==n_ext]
+            if not current_id:
                 continue
             if min_nexternal != max_nexternal:
                 if n_ext == min_nexternal:
                     text.append('       if (npdg.eq.%i)then' % n_ext)
                 else:
                     text.append('       else if (npdg.eq.%i)then' % n_ext)
-            for ii,pdgs in enumerate(current):
+            for ii,pdgs in enumerate(current_id):
+                pid = current_pid[ii]
                 condition = '.and.'.join(['%i.eq.pdgs(%i)' %(pdg, i+1) for i, pdg in enumerate(pdgs)])
                 if ii==0:
-                    text.append( ' if(%s) then ! %i' % (condition, i))
+                    text.append( ' if(%s.and.(procid.le.0.or.procid.eq.%d)) then ! %i' % (condition, pid, i))
                 else:
-                    text.append( ' else if(%s) then ! %i' % (condition,i))
-                text.append(' call %ssmatrixhel(p, nhel, ans)' % self.prefix_info[pdgs][0])
+                    text.append( ' else if(%s.and.(procid.le.0.or.procid.eq.%d)) then ! %i' % (condition,pid,i))
+                text.append(' call %ssmatrixhel(p, nhel, ans)' % self.prefix_info[(pdgs,pid)][0])
             text.append(' endif')
         #close the function
         if min_nexternal != max_nexternal:
@@ -2192,8 +2236,9 @@
                           'maxpart': max_nexternal,
                           'nb_me': len(allids),
                           'pdgs': ','.join(str(pdg[i]) if i<len(pdg) else '0' 
-                                             for i in range(max_nexternal) for pdg in allids),
+                                           for i in range(max_nexternal) for (pdg,pid) in allids),
                           'prefix':'\',\''.join(allprefix),
+                          'pids': ','.join(str(pid) for (pdg,pid) in allids),
                           'parameter_setup': '\n'.join(parameter_setup),
                           }
         formatting['lenprefix'] = len(formatting['prefix'])
@@ -2311,7 +2356,7 @@
                 raise Exception, '--prefix options supports only \'int\' and \'proc\''
             for proc in matrix_element.get('processes'):
                 ids = [l.get('id') for l in proc.get('legs_with_decays')]
-                self.prefix_info[tuple(ids)] = [proc_prefix, proc.get_tag()] 
+                self.prefix_info[(tuple(ids), proc.get('id'))] = [proc_prefix, proc.get_tag()] 
                 
         calls = self.write_matrix_element_v4(
             writers.FortranWriter(filename),
@@ -4186,11 +4231,24 @@
                         split_orders if len(split_orders)>0 else ['ALL_ORDERS'])
         replace_dict['jamp_lines'] = '\n'.join(jamp_lines)
 
+        if self.beam_polarization:
+            replace_dict['beam_polarization'] = """
+                         DO JJ=1,nincoming
+               IF(POL(JJ).NE.1d0.AND.NHEL(JJ,I).EQ.INT(SIGN(1d0,POL(JJ)))) THEN
+                 T=T*ABS(POL(JJ))
+               ELSE IF(POL(JJ).NE.1d0)THEN
+                 T=T*(2d0-ABS(POL(JJ)))
+               ENDIF
+             ENDDO
+            """
+        else:
+            replace_dict['beam_polarization'] = ""
+
+
         replace_dict['template_file'] = pjoin(_file_path, \
                           'iolibs/template_files/%s' % self.matrix_file)
         replace_dict['template_file2'] = pjoin(_file_path, \
-                          'iolibs/template_files/split_orders_helping_functions.inc')      
-        
+                          'iolibs/template_files/split_orders_helping_functions.inc')              
         s1,s2 = matrix_element.get_spin_state_initial()
         replace_dict['nb_spin_state1'] = s1
         replace_dict['nb_spin_state2'] = s2
@@ -4297,14 +4355,21 @@
             replace_dict['define_subdiag_lines'] = ""
             replace_dict['cutsdone'] = "      cutsdone=.false.\n       cutspassed=.false."
 
+        ncomb=matrix_element.get_helicity_combinations()
+
+        replace_dict['ncomb']= ncomb
+        helicity_lines = self.get_helicity_lines(matrix_element)
+        replace_dict['helicity_lines'] = helicity_lines
+        
         if not isinstance(self, ProcessExporterFortranMEGroup):
-            ncomb=matrix_element.get_helicity_combinations()
             replace_dict['read_write_good_hel'] = self.read_write_good_hel(ncomb)
         else:
             replace_dict['read_write_good_hel'] = ""
+            # Extract helicity lines
         
         context = {'read_write_good_hel':True}
         
+        
         if writer:
             file = open(pjoin(_file_path, \
                           'iolibs/template_files/auto_dsig_v4.inc')).read()
 
=== modified file 'madgraph/iolibs/template_files/auto_dsig_v4.inc'
--- madgraph/iolibs/template_files/auto_dsig_v4.inc	2019-06-10 19:05:58 +0000
+++ madgraph/iolibs/template_files/auto_dsig_v4.inc	2020-08-13 08:01:43 +0000
@@ -156,3 +156,18 @@
 %(read_write_good_hel)s
 
 ## }
+
+
+	  integer FUNCTION GET_NHEL%(proc_id)s(hel, ipart)
+	  implicit none
+	  integer hel,i, ipart
+	  Include 'nexternal.inc'
+	  integer one_nhel(nexternal)
+	  INTEGER                 NCOMB
+	  PARAMETER (             NCOMB=%(ncomb)d)
+	  INTEGER NHEL(NEXTERNAL,NCOMB)
+	  %(helicity_lines)s
+
+	  get_nhel%(proc_id)s = nhel(ipart, hel)
+	  return
+	  end
 
=== modified file 'madgraph/iolibs/template_files/loop_optimized/check_sa_all.py.inc'
--- madgraph/iolibs/template_files/loop_optimized/check_sa_all.py.inc	2016-10-20 23:12:29 +0000
+++ madgraph/iolibs/template_files/loop_optimized/check_sa_all.py.inc	2020-08-13 08:01:43 +0000
@@ -49,10 +49,10 @@
 allmatrix2py.set_madloop_path(os.path.abspath(pjoin(root_path,'MadLoop5_resources')))
 # List all processes available in this module
 
-all_procs_available = [[pdg for pdg in proc_pdgs if pdg!=0] for proc_pdgs in allmatrix2py.get_pdg_order()]
+all_procs_available = [[(pdg,pids[i]) for i, pdg in enumerate(proc_pdgs) if pdg!=0] for (proc_pdgs, pids) in allmatrix2py.get_pdg_order()]
 print "Available processes are:"
 for process_pdgs in all_procs_available:
-    print "   --> %s"%str(process_pdgs)
+    print "   --> %s: process id: %i "%str(process_pdgs[0], process_pdgs[1])
 
 def invert_momenta(p):
     """ fortran/C-python do not order table in the same order"""
@@ -62,7 +62,7 @@
 
 # The kinematic configuration in the convention (E, px, py, pz) and with particles ordered as in the process definition.
 # This is the structure of this dictionary.
-p= dict( (tuple(process_pdgs),[[None,]*4]*len(process_pdgs)) for process_pdgs in all_procs_available)
+p= dict( (tuple(process_pdgs[0]),[[None,]*4]*len(process_pdgs[0])) for process_pdgs in all_procs_available)
 
 if (p.values()[0])[0][0] is None:
     if not os.path.isfile(pjoin(root_path,'PS.input')):
@@ -120,7 +120,7 @@
     # Choice of renormalization scale
     renormalization_scale = 91.188
 
-    finite_loop_me, return_code = allmatrix2py.smatrixhel(pdgs=list(process_pdgs),
+    finite_loop_me, return_code = allmatrix2py.smatrixhel(pdgs=list(process_pdgs), procid=-1,
 		                          p=P, alphas=alphas,scales2=renormalization_scale,nhel=nhel)
 
     print '* %-108s *'%' MadLoop evaluation for the process '
 
=== modified file 'madgraph/iolibs/template_files/loop_optimized/loop_matrix_standalone.inc'
--- madgraph/iolibs/template_files/loop_optimized/loop_matrix_standalone.inc	2017-07-31 22:08:40 +0000
+++ madgraph/iolibs/template_files/loop_optimized/loop_matrix_standalone.inc	2020-08-13 08:01:43 +0000
@@ -2737,11 +2737,14 @@
       PARAMETER (NEXTERNAL=%(nexternal)d)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=%(nSquaredSO)d)
+      %(get_nsqso_born)s
 C  
 C ARGUMENTS 
 C  
       %(real_dp_format)s P(0:3,NEXTERNAL)
-      %(real_dp_format)s ANS(0:3,0:NSQUAREDSO)
+	  INTEGER ANS_DIMENSION
+	  PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      %(real_dp_format)s ANS(0:3,0:ANS_DIMENSION)
 	  INTEGER HEL, USERHEL
 	  common/%(proc_prefix)sUSERCHOICE/USERHEL
 C ----------
@@ -2764,7 +2767,10 @@
 C ARGUMENTS 
 C  
       %(real_dp_format)s P(0:3,NEXTERNAL)
-      %(real_dp_format)s ANS(0:3,0:NSQUAREDSO)
+      %(get_nsqso_born)s
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      %(real_dp_format)s ANS(0:3,0:ANS_DIMENSION)
 	  INTEGER HEL, RET_CODE
 	  %(real_dp_format)s PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C
@@ -2875,7 +2881,10 @@
 C ARGUMENTS 
 C  
       %(real_dp_format)s P(0:3,NEXTERNAL)
-      %(real_dp_format)s ANS(0:3,0:NSQUAREDSO)
+      %(get_nsqso_born)s
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      %(real_dp_format)s ANS(0:3,0:ANS_DIMENSION)
 	  %(real_dp_format)s PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 	  INTEGER RET_CODE	  
 C
 
=== modified file 'madgraph/iolibs/template_files/makefile_sa_f2py'
--- madgraph/iolibs/template_files/makefile_sa_f2py	2017-09-25 21:30:36 +0000
+++ madgraph/iolibs/template_files/makefile_sa_f2py	2020-08-13 08:01:43 +0000
@@ -2,7 +2,8 @@
 SHELL = /bin/bash
 LIBDIR = ../lib/
 PROG   = check
-LINKLIBS =  -L../lib/ -ldhelas -lmodel
+LINKLIBS_ME =  -L../lib/ -ldhelas -lmodel
+LINKLIBS_ALL =  -L../lib/ -lmatrix -ldhelas -lmodel
 LIBS = $(LIBDIR)libdhelas.$(libext) $(LIBDIR)libmodel.$(libext)
 PROCESS=  allmatrix.o
 LIBRARY = libmatrix.a
@@ -12,12 +13,19 @@
   MENUM=2
 endif
 
-
-allmatrix$(MENUM)py.so: $(LIBDIR)/$(LIBRARY)  all_matrix.o
-	$(F2PY) --f77exec=$(FC) $(LINKLIBS) -lmatrix  -c all_matrix.f */matrix.f -m allmatrix$(MENUM)py
-
-
-
+%$(MENUM)py.so: %/matrix.f
+	$(F2PY) --f77exec=$(FC) $(LINKLIBS_ME) -c $^ -m $(patsubst %/matrix.f,%$(MENUM)py,$^)
+	
+all_matrix$(MENUM)py.so: $(LIBDIR)/$(LIBRARY) all_matrix.f
+	$(F2PY) --f77exec=$(FC) $(LINKLIBS_ALL) -c all_matrix.f -m all_matrix$(MENUM)py
+
+    
+allmatrix$(MENUM)py.so: $(patsubst %/matrix.f,%$(MENUM)py.so,$(wildcard */matrix.f)) all_matrix$(MENUM)py.so
+	for lib in $(patsubst %.so,%,$^); do \
+		echo "from $$lib import *" >> allmatrix$(MENUM)py.py; \
+	done
+	
+	
 $(LIBDIR)/$(LIBRARY): $(patsubst %.f,%.o,$(wildcard */matrix.f))  all_matrix.o
 	$(call CREATELIB, $@, $^)  
 
 
=== modified file 'madgraph/iolibs/template_files/matrix_madevent_group_v4.inc'
--- madgraph/iolibs/template_files/matrix_madevent_group_v4.inc	2019-06-26 10:35:18 +0000
+++ madgraph/iolibs/template_files/matrix_madevent_group_v4.inc	2020-08-13 08:01:43 +0000
@@ -123,13 +123,7 @@
       DO I=1,NCOMB
          IF (GOODHEL(I,IMIRROR) .OR. NTRY(IMIRROR).LE.MAXTRIES.or.(ISUM_HEL.NE.0).or.THIS_NTRY(IMIRROR).le.2) THEN
             T=MATRIX%(proc_id)s(P ,NHEL(1,I),JC(1))
-            DO JJ=1,nincoming
-               IF(POL(JJ).NE.1d0.AND.NHEL(JJ,I).EQ.INT(SIGN(1d0,POL(JJ)))) THEN
-                  T=T*ABS(POL(JJ))*NB_SPIN_STATE_IN(JJ)/2d0 ! NB_SPIN_STATE(JJ)/2d0 is added for polarised beam
-               ELSE IF(POL(JJ).NE.1d0)THEN
-                  T=T*(2d0-ABS(POL(JJ)))*NB_SPIN_STATE_IN(JJ)/2d0
-               ENDIF
-            ENDDO
+	    %(beam_polarization)s
 			IF (ISUM_HEL.NE.0.and.DS_get_dim_status('Helicity').eq.0.and.ALLOW_HELICITY_GRID_ENTRIES) then
 			  call DS_add_entry('Helicity',I,T)
 			endif
@@ -174,14 +168,7 @@
             I = HEL_PICKED
 
             T=MATRIX%(proc_id)s(P ,NHEL(1,I),JC(1))
-
-            DO JJ=1,nincoming
-              IF(POL(JJ).NE.1d0.AND.NHEL(JJ,I).EQ.INT(SIGN(1d0,POL(JJ)))) THEN
-                T=T*ABS(POL(JJ))
-              ELSE IF(POL(JJ).NE.1d0)THEN
-                T=T*(2d0-ABS(POL(JJ)))
-              ENDIF
-            ENDDO
+	    %(beam_polarization)s
 c           Always one helicity at a time
             ANS = T
 c           Include the Jacobian from helicity sampling
 
=== modified file 'madgraph/iolibs/template_files/matrix_madevent_v4.inc'
--- madgraph/iolibs/template_files/matrix_madevent_v4.inc	2019-09-26 11:09:58 +0000
+++ madgraph/iolibs/template_files/matrix_madevent_v4.inc	2020-08-13 08:01:43 +0000
@@ -103,13 +103,7 @@
         DO I=1,NCOMB
            IF (GOODHEL(I) .OR. NTRY .LE. MAXTRIES.OR.(ISUM_HEL.NE.0)) THEN
                T=MATRIX%(proc_id)s(P ,NHEL(1,I),JC(1))            
-             DO JJ=1,nincoming
-               IF(POL(JJ).NE.1d0.AND.NHEL(JJ,I).EQ.INT(SIGN(1d0,POL(JJ)))) THEN
-                 T=T*ABS(POL(JJ))*NB_SPIN_STATE(JJ)/2d0 ! NB_SPIN_STATE(JJ)/2d0 is added for polarised beam
-               ELSE IF(POL(JJ).NE.1d0)THEN
-                 T=T*(2d0-ABS(POL(JJ)))*NB_SPIN_STATE(JJ)/2d0
-               ENDIF
-             ENDDO
+	       %(beam_polarization)s
 			 IF (ISUM_HEL.NE.0) then
 			   call DS_add_entry('Helicity',I,T)
 			 endif
@@ -153,13 +147,7 @@
             I = HEL_PICKED
 			
 			T=MATRIX%(proc_id)s(P ,NHEL(1,I),JC(1))            
-            DO JJ=1,nincoming
-              IF(POL(JJ).NE.1d0.AND.NHEL(JJ,I).EQ.INT(SIGN(1d0,POL(JJ)))) THEN
-                T=T*ABS(POL(JJ))
-              ELSE IF(POL(JJ).NE.1d0)THEN
-                T=T*(2d0-ABS(POL(JJ)))
-              ENDIF
-            ENDDO
+	      %(beam_polarization)s
 c           Always one helicity at a time
             ANS = T
 c           Include the Jacobian from helicity sampling
 
=== modified file 'madgraph/iolibs/template_files/super_auto_dsig_group_v4.inc'
--- madgraph/iolibs/template_files/super_auto_dsig_group_v4.inc	2020-01-16 07:19:30 +0000
+++ madgraph/iolibs/template_files/super_auto_dsig_group_v4.inc	2020-08-13 08:01:43 +0000
@@ -506,5 +506,10 @@
 C
 %(read_write_good_hel)s
 
-
-
+      integer function get_nhel(hel, partid)
+	implicit none
+	integer hel, partid
+	write(*,*) "do not support this type of pdf with group_subprocess=True"
+	stop 1
+	return
+		end
 
=== modified file 'madgraph/loop/loop_exporters.py'
--- madgraph/loop/loop_exporters.py	2020-06-20 21:10:10 +0000
+++ madgraph/loop/loop_exporters.py	2020-08-13 08:01:43 +0000
@@ -334,11 +334,12 @@
         CALL SETMADLOOPPATH(PATH)
       END
 
-  subroutine smatrixhel(pdgs, npdg, p, ALPHAS, SCALES2, nhel, ANS, RETURNCODE)
+  subroutine smatrixhel(pdgs, procid, npdg, p, ALPHAS, SCALES2, nhel, ANS, RETURNCODE)
   IMPLICIT NONE
 
 CF2PY double precision, intent(in), dimension(0:3,npdg) :: p
 CF2PY integer, intent(in), dimension(npdg) :: pdgs
+CF2PY integer, intent(in):: procid
 CF2PY integer, intent(in) :: npdg
 CF2PY double precision, intent(out) :: ANS
 CF2PY integer, intent(out) :: RETURNCODE
@@ -346,22 +347,25 @@
 CF2PY double precision, intent(in) :: SCALES2
 
   integer pdgs(*)
-  integer npdg, nhel, RETURNCODE
+  integer npdg, nhel, RETURNCODE, procid
   double precision p(*)
   double precision ANS, ALPHAS, PI,SCALES2
-
+ 1 continue
 %(smatrixhel)s
 
       return
       end
   
-  subroutine get_pdg_order(OUT)
+  subroutine get_pdg_order(OUT, ALLPROC)
   IMPLICIT NONE
 CF2PY INTEGER, intent(out) :: OUT(%(nb_me)i,%(maxpart)i)  
-  
+CF2PY INTEGER, intent(out) :: ALLPROC(%(nb_me)i)
   INTEGER OUT(%(nb_me)i,%(maxpart)i), PDGS(%(nb_me)i,%(maxpart)i)
+  INTEGER ALLPROC(%(nb_me)i),PIDs(%(nb_me)i)
   DATA PDGS/ %(pdgs)s /
+  DATA PIDS/ %(pids)s /
   OUT=PDGS
+  ALLPROC = PIDS
   RETURN
   END
   
@@ -378,31 +382,37 @@
          
         allids = self.prefix_info.keys()
         allprefix = [self.prefix_info[key][0] for key in allids]
-        min_nexternal = min([len(ids) for ids in allids])
-        max_nexternal = max([len(ids) for ids in allids])
+        min_nexternal = min([len(ids[0]) for ids in allids])
+        max_nexternal = max([len(ids[0]) for ids in allids])
 
         info = []
-        for key, (prefix, tag) in self.prefix_info.items():
-            info.append('#PY %s : %s # %s' % (tag, key, prefix))
+        for (key,pid), (prefix, tag) in self.prefix_info.items():
+            info.append('#PY %s : %s # %s %s' % (tag, key, prefix, pid))
             
 
         text = []
         for n_ext in range(min_nexternal, max_nexternal+1):
-            current = [ids for ids in allids if len(ids)==n_ext]
-            if not current:
+            current_id = [ids[0] for ids in allids if len(ids[0])==n_ext]
+            current_pid = [ids[1] for ids in allids if len(ids[0])==n_ext]
+            if not current_id:
                 continue
             if min_nexternal != max_nexternal:
                 if n_ext == min_nexternal:
                     text.append('       if (npdg.eq.%i)then' % n_ext)
                 else:
                     text.append('       else if (npdg.eq.%i)then' % n_ext)
-            for ii,pdgs in enumerate(current):
+            for ii,pdgs in enumerate(current_id):
+                pid =  current_pid[ii]
                 condition = '.and.'.join(['%i.eq.pdgs(%i)' %(pdg, i+1) for i, pdg in enumerate(pdgs)])
                 if ii==0:
-                    text.append( ' if(%s) then ! %i' % (condition, i))
+                    text.append( ' if(%s.and.(procid.le.0.or.procid.eq.%d)) then ! %i' % (condition, pid, i))
                 else:
-                    text.append( ' else if(%s) then ! %i' % (condition,i))
-                text.append(' call %sget_me(p, ALPHAS, DSQRT(SCALES2), NHEL, ANS, RETURNCODE)' % self.prefix_info[pdgs][0])
+                    text.append( ' else if(%s.and.(procid.le.0.or.procid.eq.%d)) then ! %i' % (condition,pid,i))
+                text.append(' call %sget_me(p, ALPHAS, DSQRT(SCALES2), NHEL, ANS, RETURNCODE)' % self.prefix_info[(pdgs,pid)][0])
+            text.append( ' else if(procid.gt.0) then !')
+            text.append( ' procid = -1' )
+            text.append( ' goto 1' )
+            
             text.append(' endif')
         #close the function
         if min_nexternal != max_nexternal:
@@ -422,9 +432,10 @@
                           'nb_me': len(allids),
                           'pdgs': ','.join([str(pdg[i]) if i<len(pdg) else '0' 
                                              for i in range(max_nexternal) \
-                                             for pdg in allids]),
+                                             for (pdg,pid) in allids]),
                       'prefix':'\',\''.join(allprefix),
                       'parameter_setup': '\n'.join(parameter_setup),
+                      'pids':  ','.join(str(pid) for (pdg,pid) in allids),
                       }
     
     
@@ -1023,7 +1034,7 @@
         if 'prefix' in self.cmd_options and self.cmd_options['prefix'] in ['int','proc']:
             for proc in matrix_element.get('processes'):
                 ids = [l.get('id') for l in proc.get('legs_with_decays')]
-                self.prefix_info[tuple(ids)] = [dict['proc_prefix'], proc.get_tag()]
+                self.prefix_info[tuple(ids),proc.get('id')] = [dict['proc_prefix'], proc.get_tag()]
 
         # The proc_id is used for MadEvent grouping, so none of our concern here
         # and it is simply set to an empty string.        
@@ -3013,7 +3024,6 @@
         matrix_element.rep_dict['coef_construction']=replace_dict['coef_construction']            
         
         replace_dict['coef_merging']='\n'.join(coef_merging)
-
         file = file % replace_dict
         number_of_calls = len(filter(lambda call: call.find('CALL LOOP') != 0, \
                                                                  loop_CT_calls))   
 
=== modified file 'madgraph/various/banner.py'
--- madgraph/various/banner.py	2020-05-25 21:16:05 +0000
+++ madgraph/various/banner.py	2020-08-13 08:01:43 +0000
@@ -24,6 +24,7 @@
 import math
 import StringIO
 import itertools
+import time
 
 
 pjoin = os.path.join
@@ -2896,10 +2897,10 @@
         self.add_param("nevents", 10000)        
         self.add_param("iseed", 0)
         self.add_param("python_seed", -2, include=False, hidden=True, comment="controlling python seed [handling in particular the final unweighting].\n -1 means use default from random module.\n -2 means set to same value as iseed")
-        self.add_param("lpp1", 1, fortran_name="lpp(1)", allowed=[-1,1,0,2,3,9, -2,-3],
-                        comment='first beam energy distribution:\n 0: fixed energy\n 1: PDF from proton\n -1: PDF from anti-proton\n 2:photon from proton, 3:photon from electron, 9: PLUGIN MODE')
-        self.add_param("lpp2", 1, fortran_name="lpp(2)", allowed=[-1,1,0,2,3,9],
-                       comment='first beam energy distribution:\n 0: fixed energy\n 1: PDF from proton\n -1: PDF from anti-proton\n 2:photon from proton, 3:photon from electron, 9: PLUGIN MODE')
+        self.add_param("lpp1", 1, fortran_name="lpp(1)", allowed=[-1,1,0,2,3,4,9,-2,-3,-4],
+                        comment='first beam energy distribution:\n 0: fixed energy\n 1: PDF from proton\n -1: PDF from anti-proton\n 2:photon from proton, 3:photon/W/Z from electron, 4:photon/W/Z from muon,\n 9: PLUGIN MODE')
+        self.add_param("lpp2", 1, fortran_name="lpp(2)", allowed=[-1,1,0,2,3,4,9,-3,-4],
+                       comment='first beam energy distribution:\n 0: fixed energy\n 1: PDF from proton\n -1: PDF from anti-proton\n 2:photon from proton, 3:photon/W/Z from electron, 4:photon/W/Z from muon,\n 9: PLUGIN MODE')
         self.add_param("ebeam1", 6500.0, fortran_name="ebeam(1)")
         self.add_param("ebeam2", 6500.0, fortran_name="ebeam(2)")
         self.add_param("polbeam1", 0.0, fortran_name="pb1", hidden=True,
@@ -3182,7 +3183,6 @@
         if self['xqcut'] > 0:
             if self['ickkw'] == 0:
                 logger.error('xqcut>0 but ickkw=0. Potentially not fully consistent setup. Be carefull')
-                import time
                 time.sleep(5)
             if self['drjj'] != 0:
                 if 'drjj' in self.user_set:
@@ -3212,14 +3212,16 @@
 
         # check if lpp = 
         for i in [1,2]:
-            if self['lpp%s' % i ] == 3 and self['dsqrt_q2fact%s'%i] > 4:
-                raise InvalidRunCard( "Photon from electron are using fixed scale value of muf [dsqrt_q2fact%s] as the cut off value of the approximation.\n" % i + \
-                                      "For EPA this number should be small (for HERA prediction it should be 2 at most)")
-            if self['lpp%s' % i ] == 2 and self['dsqrt_q2fact%s'%i] == 91.188:
-                raise InvalidRunCard("Since 2.7.1 Photon from proton are using fixed scale value of muf [dsqrt_q2fact%s] as the cut of th Improved Weizsaecker-Williams formula. Please edit it accordingly." % i)
+            if abs(self['lpp%s' % i ]) in [3,4] and self['dsqrt_q2fact%s'%i] == 91.188:
+                logger.warning("Photon from lepton are using fixed scale value of muf [dsqrt_q2fact%s] as the cut of the EPA. Looks like you kept the default value (Mz). Is this really the cut-off of the EPA that you want to use?" % i)
+                time.sleep(5)
         
+            if abs(self['lpp%s' % i ]) == 2 and self['dsqrt_q2fact%s'%i] == 91.188:
+                logger.warning("Since 2.7.1 Photon from proton are using fixed scale value of muf [dsqrt_q2fact%s] as the cut of the Improved Weizsaecker-Williams formula. Please edit it accordingly." % i)
+                time.sleep(5)
+                
         # if both lpp1/2 are on PA mode -> force fixed factorization scale
-        if self['lpp1'] in [2, 3] and self['lpp2'] in [2, 3] and not self['fixed_fac_scale']:
+        if abs(self['lpp1']) in [2, 3,4] and abs(self['lpp2']) in [2, 3,4] and not self['fixed_fac_scale']:
             raise InvalidRunCard("Having both beam in elastic photon mode requires fixec_fac_scale to be on True [since this is use as cutoff]")
 
 
@@ -3293,8 +3295,11 @@
           process 1->N all cut set on off.
           loop_induced -> MC over helicity
           e+ e- beam -> lpp:0 ebeam:500
+          P  X beam (X: neutrino, lepton, EWA)
+          support of EWA beam [ZW] (from proton by default -> 100TeV) 
           p p beam -> set maxjetflavor automatically
           more than one multiplicity: ickkw=1 xqcut=30 use_syst=F
+          interference module -> no systematics
          """
 
         if proc_characteristic['loop_induced']:
@@ -3318,11 +3323,87 @@
                             beam_id.add(leg['id'])
 
             if any(i in beam_id for i in [1,-1,2,-2,3,-3,4,-4,5,-5,21,22]):
-                maxjetflavor = max([4]+[abs(i) for i in beam_id if  -7< i < 7])
+                maxjetflavor = max([3]+[abs(i) for i in beam_id if  -7< i < 7])
                 self['maxjetflavor'] = maxjetflavor
                 self['asrwgtflavor'] = maxjetflavor
             
-            if any(i in beam_id for i in [1,-1,2,-2,3,-3,4,-4,5,-5,21,22]):
+            # check EWA case  -> for pp at 100TeV or EWA-lepton (1+1 TeV)            
+            if any(i in beam_id for i in [23,24,-24]):
+
+                self['use_syst'] = False
+                self['nhel'] = 1
+
+                # check where the W/Z are (on which beam)                
+                EWA_in1 =  any(i in beam_id_split[0] for i in [23,24,-24])
+                EWA_in2 = any(i in beam_id_split[1] for i in [23,24,-24])
+                p_in1 = any(i in beam_id_split[0] for i in [1,-1,2,-2,3,-3,4,-4,5,-5,21])
+                p_in2 = any(i in beam_id_split[1] for i in [1,-1,2,-2,3,-3,4,-4,5,-5,21]) 
+                
+                
+                # case1 W W collision 
+                if EWA_in1 and EWA_in2:
+                    self['lpp1'] = 2    
+                    self['lpp2'] = 2
+                    self['ebeam1'] = '50k'
+                    self['ebeam2'] = '50k'
+                # case2 p W collision
+                elif p_in1 and EWA_in2:
+                    self['lpp1'] = 1    
+                    self['lpp2'] = 2
+                    self['ebeam1'] = '50k'
+                    self['ebeam2'] = '50k' 
+                # case3  W p collision                   
+                elif p_in2 and EWA_in1:
+                    self['lpp1'] = 2    
+                    self['lpp2'] = 1
+                    self['ebeam1'] = '50k'
+                    self['ebeam2'] = '50k'
+                # case3  W X collision   (mainly for X=lepton)                                  
+                elif EWA_in1:
+                    self['ebeam1'] = '1k'  
+                    self['ebeam2'] = '1k' 
+                    if 11 in beam_id_split[1]:
+                        self['lpp1'] = -3  
+                        self['lpp2'] = 0  
+                    elif -11 in beam_id_split[1]:
+                        self['lpp1'] = 3  
+                        self['lpp2'] = 0  
+                    elif 13 in beam_id_split[1]:                                 
+                        self['lpp1'] = -4  
+                        self['lpp2'] = 0  
+                    elif -13 in beam_id_split[1]:
+                        self['lpp1'] = 4  
+                        self['lpp2'] = 0 
+                    else:
+                        self['lpp1'] = 0  
+                        self['lpp2'] = 0  
+                # case3   X W collision   (mainly for X=lepton)                       
+                elif EWA_in2:
+                    self['ebeam1'] = '1k'  
+                    self['ebeam2'] = '1k' 
+                    if 11 in beam_id_split[0]:
+                        self['lpp1'] = 0  
+                        self['lpp2'] = -3  
+                    elif -11 in beam_id_split[0]:
+                        self['lpp1'] = 0 
+                        self['lpp2'] = 3  
+                    elif 13 in beam_id_split[0]:                                 
+                        self['lpp1'] = 0  
+                        self['lpp2'] = -4  
+                    elif -13 in beam_id_split[0]:
+                        self['lpp1'] = 0 
+                        self['lpp2'] = 4
+                    else:
+                        self['lpp1'] = 0  
+                        self['lpp2'] = 0
+                # anything else  should it happen?
+                else:
+                    self['nhel'] = 1
+                    self['lpp1'] = 0  
+                    self['lpp2'] = 0  
+                    self['ebeam1'] = '1k'  
+                    self['ebeam2'] = '1k'               
+            elif any(i in beam_id for i in [1,-1,2,-2,3,-3,4,-4,5,-5,21,22]):
                 # check for e p collision
                 if any(id  in beam_id for id in [11,-11,13,-13]):
                     self.display_block.append('beam_pol')
@@ -3336,20 +3417,25 @@
                         self['lpp2'] = 0  
                         self['ebeam1'] = '6500'  
                         self['ebeam2'] = '1k'  
-            
-            elif 11 in beam_id or -11 in beam_id:
+
+            elif any(id in beam_id for id in [11,-11,13,-13]):
                 self['lpp1'] = 0
                 self['lpp2'] = 0
                 self['ebeam1'] = 500
                 self['ebeam2'] = 500
                 self['use_syst'] = False
+                if beam_id_split[0] == beam_id_split[1]:
+                    self.display_block.append('ecut')
                 self.display_block.append('beam_pol')
-                self.display_block.append('ecut')
+
+                if beam_id_split[0] == beam_id_split[1]:
+                    self.display_block.append('ecut')
             else:
                 self['lpp1'] = 0
                 self['lpp2'] = 0    
                 self['use_syst'] = False   
-                self.display_block.append('beam_pol')         
+                self.display_block.append('beam_pol')  
+                self.display_block.append('ecut')       
             
             # automatic polarisation of the beam if neutrino beam  
             if any(id  in beam_id for id in [12,-12,14,-14,16,-16]):
@@ -3358,25 +3444,25 @@
                     self['lpp1'] = 0   
                     self['ebeam1'] = '1k'  
                     self['polbeam1'] = -100
-                    if not all(id  in beam_id_split[0] for id in [12,14,16]):
-                        logger.warning('Issue with default beam setup of neutrino in the run_card. Please check it up [polbeam1].')
+                    if not all(id  in [12,14,16] for id in beam_id_split[0]):
+                        logger.warning('Issue with default beam setup of neutrino in the run_card. Please check it up [polbeam1]. %s')
                 elif any(id  in beam_id_split[0] for id in [-12,-14,-16]):
                     self['lpp1'] = 0   
                     self['ebeam1'] = '1k'  
                     self['polbeam1'] = 100
-                    if not all(id  in beam_id_split[0] for id in [-12,-14,-16]):
+                    if not all(id  in [-12,-14,-16] for id in beam_id_split[0]):
                         logger.warning('Issue with default beam setup of neutrino in the run_card. Please check it up [polbeam1].')                         
                 if any(id  in beam_id_split[1] for id in [12,14,16]):
                     self['lpp2'] = 0   
                     self['ebeam2'] = '1k'  
                     self['polbeam2'] = -100
-                    if not all(id  in beam_id_split[1] for id in [12,14,16]):
+                    if not all(id  in [12,14,16] for id in beam_id_split[1]):
                         logger.warning('Issue with default beam setup of neutrino in the run_card. Please check it up [polbeam2].')
                 if any(id  in beam_id_split[1] for id in [-12,-14,-16]):
                     self['lpp2'] = 0   
                     self['ebeam2'] = '1k'  
                     self['polbeam2'] = 100
-                    if not all(id  in beam_id_split[1] for id in [-12,-14,-16]):
+                    if not all(id  in [-12,-14,-16] for id in beam_id_split[1]):
                         logger.warning('Issue with default beam setup of neutrino in the run_card. Please check it up [polbeam2].')
             
         # Check if need matching
@@ -4285,7 +4371,7 @@
             maxjetflavor = max([4]+[abs(i) for i in beam_id if  -7< i < 7])
             self['maxjetflavor'] = maxjetflavor
             pass
-        elif 11 in beam_id or -11 in beam_id:
+        elif any(id in beam_id for id in [11,-11,13,-13]):
             self['lpp1'] = 0
             self['lpp2'] = 0
             self['ebeam1'] = 500
 
=== modified file 'madgraph/various/cluster.py'
--- madgraph/various/cluster.py	2017-11-17 14:15:01 +0000
+++ madgraph/various/cluster.py	2020-08-13 08:01:43 +0000
@@ -87,6 +87,7 @@
     """Basic Class for all cluster type submission"""
     name = 'mother class'
     identifier_length = 14
+    allow_arrays = False
 
     def __init__(self,*args, **opts):
         """Init the cluster"""
@@ -97,6 +98,8 @@
         self.submitted_dirs = [] #HTCaaS
         self.submitted_exes = [] #HTCaaS
         self.submitted_args = [] #HTCaaS
+        self.waiting_array_job = [] # for job array submission
+        
 
         if 'cluster_queue' in opts:
             self.cluster_queue = opts['cluster_queue']
@@ -122,7 +125,6 @@
         """How to make one submission. Return status id on the cluster."""
         raise NotImplemented, 'No implementation of how to submit a job to cluster \'%s\'' % self.name
 
-
     @store_input()
     def submit2(self, prog, argument=[], cwd=None, stdout=None, stderr=None, 
                 log=None, input_files=[], output_files=[], required_output=[],
@@ -199,10 +201,15 @@
 
     def cluster_submit(self, prog, argument=[], cwd=None, stdout=None, stderr=None, 
                 log=None, input_files=[], output_files=[], required_output=[],
-                nb_submit=0, packet_member=None):
+                nb_submit=0, packet_member=None, allow_array=False):
             """This function wrap the cluster submition with cluster independant
                method should not be overwritten (but for DAG type submission)"""
-               
+            
+            if allow_array and self.allow_array:
+                self.add_into_array(prog, argument, cwd, stdout, stderr, log, input_files, 
+                             output_files, required_output, nb_submit)
+                return None   
+            
             id = self.submit2(prog, argument, cwd, stdout, stderr, log, input_files, 
                              output_files, required_output, nb_submit)               
                
@@ -458,7 +465,27 @@
                 self.submit2(**args)
             return 'resubmit'
         return 'done'
-            
+
+    def add_into_array(self,*args, **opts):
+        """ """
+        raise NotImplemented, 'No implementation of how to submit  with a job array to cluster \'%s\'' % self.name
+
+    def submit_array(self, list_of_args):
+        """ """
+        raise NotImplemented, 'No implementation of how to submit  with a job array to cluster \'%s\'' % self.name
+
+
+    def release_array(self):
+        
+        if not self.allow_arrays:
+            if self.waiting_array_jobs:
+                raise NotImplemented, 'No implementation of how to submit  with a job array to cluster \'%s\'' % self.name
+            return None
+        elif self.waiting_array_jobs:
+            return self.submit_array(self.waiting_array_jobs)
+        else:
+            return None
+
     @check_interupt()
     def launch_and_wait(self, prog, argument=[], cwd=None, stdout=None, 
                         stderr=None, log=None, required_output=[], nb_submit=0,
@@ -1648,6 +1675,7 @@
     running_tag = ['R', 'CG']
     complete_tag = ['C']
     identifier_length = 8
+    allow_arrays = True
 
     @multiple_try()
     def submit(self, prog, argument=[], cwd=None, stdout=None, stderr=None, log=None,
@@ -1692,6 +1720,42 @@
         self.submitted_ids.append(id)
         return id
 
+    def add_into_array(prog, argument, cwd, stdout, stderr, log, input_files, 
+                             output_files, required_output, nb_submit):
+        
+        self.waiting_array_job.append((prog, argument, cwd, stdout, stderr, log, input_files, 
+                             output_files, required_output, nb_submit))
+        
+
+    def submit_array(self, waiting_array_jobs):
+        """write a temporary file in order to submit an array of jobs"""
+    
+        for arg in waiting_array_jobs):
+            print( arg)
+        misc.sprint(waiting_array_jobs)
+        template= """
+#SBATCH --job-name=%(job_name)s   # Job name
+#SBATCH --mail-type=NONE         # Mail events (NONE, BEGIN, END, FAIL, ALL)
+#SBATCH --mail-user=[e-mail] # Where to send mail    
+#SBATCH --nodes=1                   # Use one node
+#SBATCH --ntasks=1                  # Run a single task
+#SBATCH --output=array_%A-%a.out    # Standard output and error log
+#SBATCH --array=1-5                 # Array range
+        """
+    
+        import os
+        import tempfile
+
+        fd, path = tempfile.mkstemp()
+        try:
+            with os.fdopen(fd, 'w') as tmp:
+                # do stuff with temp file
+                tmp.write('stuff')
+        finally:        
+            os.remove(path)
+
+
+
     @multiple_try()
     def control_one_job(self, id):
         """ control the status of a single job with it's cluster id """
 
=== modified file 'mg5decay/decay_objects.py'
--- mg5decay/decay_objects.py	2018-06-01 08:46:21 +0000
+++ mg5decay/decay_objects.py	2020-08-13 08:01:43 +0000
@@ -4086,15 +4086,28 @@
                 # Evaluate the coupling strength
                 vertex =  model.get('interaction_dict')[abs(vert.get('id'))]
                 lorentz_factor = 0
+
+                for fct in model['functions']:
+                    if ';' in fct.expr:
+                        raise Exception('; is not allowed in function_library')
+                    exec("mdl_%s = lambda %s: %s" %(fct.name, ','.join( fct.arguments), fct.expr))
+                    exec("%s = lambda %s: %s" %(fct.name, ','.join( fct.arguments), fct.expr))                    
+
                 for key, v in vertex['couplings'].items():
                     if not hasattr(model, 'lorentz_dict'):
                         model.lorentz_dict = dict([(l.name, l) for l in model['lorentz']])
                         self.init_regular_expression()
+                    
+                    lorentz =  model.lorentz_dict[vertex['lorentz'][key[1]]]
+                    structure = lorentz.structure
+                    if hasattr(lorentz, 'formfactors') and lorentz.formfactors:
+                        for ff in lorentz.formfactors:
+                            structure = structure.replace(ff.name, '(%s)' % ff.value)
                         
-                    structure = model.lorentz_dict[vertex['lorentz'][key[1]]].structure
                     new_structure = self.lor_pattern.sub(self.simplify_lorentz,
                                                          structure)
-                    lor_value = eval(new_structure % q_dict_lor)
+
+                    lor_value = eval(new_structure % q_dict_lor, locals(), globals())
                     if lor_value == 0:
                         new_structure = new_structure.replace('-','+')
                         lor_value = eval(new_structure % q_dict_lor)
@@ -4131,6 +4144,12 @@
         # A quick estimate of the next-level decay of a off-shell decay
         # Consider all legs are onshell.
         else:
+            for fct in model['functions']:
+                if ';' in fct.expr:
+                    raise Exception('; is not allowed in function_library')
+                exec("mdl_%s = lambda %s: %s" %(fct.name, ','.join( fct.arguments), fct.expr))
+                exec("%s = lambda %s: %s" %(fct.name, ','.join( fct.arguments), fct.expr))                    
+
             M = abs(eval(ini_part.get('mass')))
             # The avg_E is lower by one more particle in the next-level.
             avg_E = (M/(len(self.get_final_legs())+1.))
 
=== modified file 'tests/acceptance_tests/test_cmd_reweight.py'
--- tests/acceptance_tests/test_cmd_reweight.py	2020-06-20 21:10:10 +0000
+++ tests/acceptance_tests/test_cmd_reweight.py	2020-08-13 08:01:43 +0000
@@ -212,8 +212,12 @@
         ff = open(pjoin(self.run_dir, 'Cards', 'reweight_card.dat'),'w')
         ff.write(cmd_lines)
         ff.close()
-        
-        with misc.stdchannel_redirected(sys.stdout, os.devnull):
+        mad_logger = logging.getLogger('madgraph')
+
+        if mad_logger.level > 20: 
+            with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                me_cmd.run_cmd('reweight run_01 --from_cards')
+        else:
             me_cmd.run_cmd('reweight run_01 --from_cards')
         
         lhe = lhe_parser.EventFile(pjoin(self.run_dir,'Events','run_01', 'events.lhe.gz'))
@@ -242,9 +246,12 @@
         ff.write(cmd_lines)
         ff.close()
         
-        if 1: #with misc.stdchannel_redirected(sys.stdout, os.devnull):
+        mad_logger = logging.getLogger('madgraph')
+        if mad_logger.level > 20: 
+            with misc.stdchannel_redirected(sys.stdout, os.devnull):
+                me_cmd.run_cmd('reweight run_01 --from_cards')
+        else:
             me_cmd.run_cmd('reweight run_01 --from_cards')
-        
         #solutions = [41.511565, 41.930505, 41.511565, 41.511565, 41.586169, 41.511565, 41.511565, 41.511565, 41.511565, 42.046806, 41.511565, 44.164503, 41.511565, -41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.778368, 41.511565, 42.05448, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 42.115494, 41.511854, 41.511567, 41.511565, 42.00028, 42.120605, 41.514867, -41.511565, 41.511565, 45.125706, 41.511565, 42.180208, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.997509, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.513926, 41.511565, 41.882499, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 47.532736, 41.511565, 41.809063, 41.511565, 41.511565, 41.511565, 41.927695, 41.511565, 41.511565, 41.555692, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, -41.511565, 41.511565, 42.029675, 41.725129, 41.511565, 41.511565, 41.511778, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 42.050439, 41.511565, 41.511565, 41.511565, -41.511565, 42.071105, 41.511565, 41.511565]
         solutions =  [41.511565, 41.75164, 41.511565, 41.511565, 41.557009, 41.511565, 41.511565, 41.511565, 41.511565, 41.759881, 41.511565, 44.717244, 41.511565, -41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 42.041827, 41.511565, 41.663092, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.920917, 41.511736, 41.511566, 41.511565, 205.29791, 42.538389, 41.513573, -41.511565, 41.511565, 44.650974, 41.511565, 42.549187, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.431553, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.513265, 41.511565, 43.902643, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 8853.0397, 41.511565, 333.38491, 41.511565, 41.511565, 41.511565, 75.701362, 41.511565, 41.511565, 41.54695, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, -41.511565, 41.511565, 41.454923, 41.53984, 41.511565, 41.511565, 41.511743, 41.511565, 41.511565, 41.511565, 41.511565, 41.511565, 41.732627, 41.511565, 41.511565, 41.511565, -41.511565, 41.395865, 41.511565, 41.511565]
         lhe = lhe_parser.EventFile(pjoin(self.run_dir,'Events','run_01', 'rwgt_events_tree_rwgt_1.lhe.gz'))
 
=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f	2020-08-13 08:01:43 +0000
@@ -2553,11 +2553,14 @@
       PARAMETER (NEXTERNAL=3)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=1)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/USERCHOICE/USERHEL
 C     ----------
@@ -2581,7 +2584,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2698,7 +2704,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%loop_matrix.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%loop_matrix.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%loop_matrix.f	2020-08-13 08:01:43 +0000
@@ -2553,11 +2553,14 @@
       PARAMETER (NEXTERNAL=3)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=1)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/USERCHOICE/USERHEL
 C     ----------
@@ -2581,7 +2584,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2698,7 +2704,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/auto_dsig.f'
--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/auto_dsig.f	2019-06-17 21:29:56 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/auto_dsig.f	2020-08-13 08:01:43 +0000
@@ -174,3 +174,33 @@
 
 
 
+
+      INTEGER FUNCTION GET_NHEL1(HEL, IPART)
+      IMPLICIT NONE
+      INTEGER HEL,I, IPART
+      INCLUDE 'nexternal.inc'
+      INTEGER ONE_NHEL(NEXTERNAL)
+      INTEGER                 NCOMB
+      PARAMETER (             NCOMB=16)
+      INTEGER NHEL(NEXTERNAL,NCOMB)
+      DATA (NHEL(I,   1),I=1,4) / 1,-1,-1, 1/
+      DATA (NHEL(I,   2),I=1,4) / 1,-1,-1,-1/
+      DATA (NHEL(I,   3),I=1,4) / 1,-1, 1, 1/
+      DATA (NHEL(I,   4),I=1,4) / 1,-1, 1,-1/
+      DATA (NHEL(I,   5),I=1,4) / 1, 1,-1, 1/
+      DATA (NHEL(I,   6),I=1,4) / 1, 1,-1,-1/
+      DATA (NHEL(I,   7),I=1,4) / 1, 1, 1, 1/
+      DATA (NHEL(I,   8),I=1,4) / 1, 1, 1,-1/
+      DATA (NHEL(I,   9),I=1,4) /-1,-1,-1, 1/
+      DATA (NHEL(I,  10),I=1,4) /-1,-1,-1,-1/
+      DATA (NHEL(I,  11),I=1,4) /-1,-1, 1, 1/
+      DATA (NHEL(I,  12),I=1,4) /-1,-1, 1,-1/
+      DATA (NHEL(I,  13),I=1,4) /-1, 1,-1, 1/
+      DATA (NHEL(I,  14),I=1,4) /-1, 1,-1,-1/
+      DATA (NHEL(I,  15),I=1,4) /-1, 1, 1, 1/
+      DATA (NHEL(I,  16),I=1,4) /-1, 1, 1,-1/
+
+      GET_NHEL1 = NHEL(IPART, HEL)
+      RETURN
+      END
+
 
=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/matrix1.f'
--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/matrix1.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/matrix1.f	2020-08-13 08:01:43 +0000
@@ -150,6 +150,7 @@
      $     .LE.MAXTRIES.OR.(ISUM_HEL.NE.0).OR.THIS_NTRY(IMIRROR).LE.2)
      $      THEN
             T=MATRIX1(P ,NHEL(1,I),JC(1))
+
             DO JJ=1,NINCOMING
               IF(POL(JJ).NE.1D0.AND.NHEL(JJ,I).EQ.INT(SIGN(1D0,POL(JJ))
      $         )) THEN
@@ -158,6 +159,7 @@
                 T=T*(2D0-ABS(POL(JJ)))*NB_SPIN_STATE_IN(JJ)/2D0
               ENDIF
             ENDDO
+
             IF (ISUM_HEL.NE.0.AND.DS_GET_DIM_STATUS('Helicity')
      $       .EQ.0.AND.ALLOW_HELICITY_GRID_ENTRIES) THEN
               CALL DS_ADD_ENTRY('Helicity',I,T)
@@ -216,6 +218,7 @@
             T=T*(2D0-ABS(POL(JJ)))
           ENDIF
         ENDDO
+
 C       Always one helicity at a time
         ANS = T
 C       Include the Jacobian from helicity sampling
 
=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/super_auto_dsig.f'
--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/super_auto_dsig.f	2020-01-19 16:22:42 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/super_auto_dsig.f	2020-08-13 08:01:43 +0000
@@ -595,6 +595,12 @@
 
 
 
-
-
+      INTEGER FUNCTION GET_NHEL(HEL, PARTID)
+      IMPLICIT NONE
+      INTEGER HEL, PARTID
+      WRITE(*,*) 'do not support this type of pdf with'
+     $ //' group_subprocess=True'
+      STOP 1
+      RETURN
+      END
 
 
=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/auto_dsig.f'
--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/auto_dsig.f	2019-06-17 21:29:56 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/auto_dsig.f	2020-08-13 08:01:43 +0000
@@ -221,3 +221,33 @@
 
 
 
+
+      INTEGER FUNCTION GET_NHEL(HEL, IPART)
+      IMPLICIT NONE
+      INTEGER HEL,I, IPART
+      INCLUDE 'nexternal.inc'
+      INTEGER ONE_NHEL(NEXTERNAL)
+      INTEGER                 NCOMB
+      PARAMETER (             NCOMB=16)
+      INTEGER NHEL(NEXTERNAL,NCOMB)
+      DATA (NHEL(I,   1),I=1,4) / 1,-1,-1,-1/
+      DATA (NHEL(I,   2),I=1,4) / 1,-1,-1, 1/
+      DATA (NHEL(I,   3),I=1,4) / 1,-1, 1,-1/
+      DATA (NHEL(I,   4),I=1,4) / 1,-1, 1, 1/
+      DATA (NHEL(I,   5),I=1,4) / 1, 1,-1,-1/
+      DATA (NHEL(I,   6),I=1,4) / 1, 1,-1, 1/
+      DATA (NHEL(I,   7),I=1,4) / 1, 1, 1,-1/
+      DATA (NHEL(I,   8),I=1,4) / 1, 1, 1, 1/
+      DATA (NHEL(I,   9),I=1,4) /-1,-1,-1,-1/
+      DATA (NHEL(I,  10),I=1,4) /-1,-1,-1, 1/
+      DATA (NHEL(I,  11),I=1,4) /-1,-1, 1,-1/
+      DATA (NHEL(I,  12),I=1,4) /-1,-1, 1, 1/
+      DATA (NHEL(I,  13),I=1,4) /-1, 1,-1,-1/
+      DATA (NHEL(I,  14),I=1,4) /-1, 1,-1, 1/
+      DATA (NHEL(I,  15),I=1,4) /-1, 1, 1,-1/
+      DATA (NHEL(I,  16),I=1,4) /-1, 1, 1, 1/
+
+      GET_NHEL = NHEL(IPART, HEL)
+      RETURN
+      END
+
 
=== modified file 'tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f'
--- tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f	2020-08-13 08:01:43 +0000
@@ -122,6 +122,7 @@
           IF (GOODHEL(I) .OR. NTRY .LE. MAXTRIES.OR.(ISUM_HEL.NE.0))
      $      THEN
             T=MATRIX(P ,NHEL(1,I),JC(1))
+
             DO JJ=1,NINCOMING
               IF(POL(JJ).NE.1D0.AND.NHEL(JJ,I).EQ.INT(SIGN(1D0,POL(JJ))
      $         )) THEN
@@ -130,6 +131,7 @@
                 T=T*(2D0-ABS(POL(JJ)))*NB_SPIN_STATE(JJ)/2D0
               ENDIF
             ENDDO
+
             IF (ISUM_HEL.NE.0) THEN
               CALL DS_ADD_ENTRY('Helicity',I,T)
             ENDIF
@@ -175,6 +177,7 @@
         I = HEL_PICKED
 
         T=MATRIX(P ,NHEL(1,I),JC(1))
+
         DO JJ=1,NINCOMING
           IF(POL(JJ).NE.1D0.AND.NHEL(JJ,I).EQ.INT(SIGN(1D0,POL(JJ))))
      $      THEN
@@ -183,6 +186,7 @@
             T=T*(2D0-ABS(POL(JJ)))
           ENDIF
         ENDDO
+
 C       Always one helicity at a time
         ANS = T
 C       Include the Jacobian from helicity sampling
 
=== modified file 'tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QCDsq_eq_4.f'
--- tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QCDsq_eq_4.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QCDsq_eq_4.f	2020-08-13 08:01:43 +0000
@@ -2630,11 +2630,14 @@
       PARAMETER (NEXTERNAL=4)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=3)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/ML5_0_USERCHOICE/USERHEL
 C     ----------
@@ -2658,7 +2661,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2775,7 +2781,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QCDsq_gt_0_QEDAmpAndQEDsq_gt_2.f'
--- tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QCDsq_gt_0_QEDAmpAndQEDsq_gt_2.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QCDsq_gt_0_QEDAmpAndQEDsq_gt_2.f	2020-08-13 08:01:43 +0000
@@ -2801,11 +2801,14 @@
       PARAMETER (NEXTERNAL=4)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=3)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/ML5_0_USERCHOICE/USERHEL
 C     ----------
@@ -2829,7 +2832,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2946,7 +2952,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QCDsq_gt_4.f'
--- tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QCDsq_gt_4.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QCDsq_gt_4.f	2020-08-13 08:01:43 +0000
@@ -2567,11 +2567,14 @@
       PARAMETER (NEXTERNAL=4)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=1)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/ML5_0_USERCHOICE/USERHEL
 C     ----------
@@ -2595,7 +2598,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2712,7 +2718,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QEDsq_le_4.f'
--- tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QEDsq_le_4.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_QEDsq_le_4.f	2020-08-13 08:01:43 +0000
@@ -2817,11 +2817,14 @@
       PARAMETER (NEXTERNAL=4)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=4)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/ML5_0_USERCHOICE/USERHEL
 C     ----------
@@ -2845,7 +2848,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2962,7 +2968,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_WGTsq_le_10_QEDAmpAndQEDsq_gt_2.f'
--- tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_WGTsq_le_10_QEDAmpAndQEDsq_gt_2.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_WGTsq_le_10_QEDAmpAndQEDsq_gt_2.f	2020-08-13 08:01:43 +0000
@@ -2801,11 +2801,14 @@
       PARAMETER (NEXTERNAL=4)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=3)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/ML5_0_USERCHOICE/USERHEL
 C     ----------
@@ -2829,7 +2832,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2946,7 +2952,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_default.f'
--- tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_default.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDQEDpert_default.f	2020-08-13 08:01:43 +0000
@@ -2630,11 +2630,14 @@
       PARAMETER (NEXTERNAL=4)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=3)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/ML5_0_USERCHOICE/USERHEL
 C     ----------
@@ -2658,7 +2661,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2775,7 +2781,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDpert_default.f'
--- tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDpert_default.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QCDpert_default.f	2020-08-13 08:01:43 +0000
@@ -2567,11 +2567,14 @@
       PARAMETER (NEXTERNAL=4)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=1)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/ML5_0_USERCHOICE/USERHEL
 C     ----------
@@ -2595,7 +2598,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2712,7 +2718,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QEDpert_default.f'
--- tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QEDpert_default.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/LoopSquaredOrder_IOTest/Loop_sqso_uux_ddx/loop_matrix_QEDpert_default.f	2020-08-13 08:01:43 +0000
@@ -2565,11 +2565,14 @@
       PARAMETER (NEXTERNAL=4)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=1)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/ML5_0_USERCHOICE/USERHEL
 C     ----------
@@ -2593,7 +2596,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2710,7 +2716,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%loop_matrix.f'
--- tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%loop_matrix.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%loop_matrix.f	2020-08-13 08:01:43 +0000
@@ -2584,11 +2584,14 @@
       PARAMETER (NEXTERNAL=4)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=1)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/ML5_0_USERCHOICE/USERHEL
 C     ----------
@@ -2612,7 +2615,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2729,7 +2735,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/TestCmdMatchBox/MatchBoxOutput/%TEST%SubProcesses%P1_uux_uux%loop_matrix.f'
--- tests/input_files/IOTestsComparison/TestCmdMatchBox/MatchBoxOutput/%TEST%SubProcesses%P1_uux_uux%loop_matrix.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/TestCmdMatchBox/MatchBoxOutput/%TEST%SubProcesses%P1_uux_uux%loop_matrix.f	2020-08-13 08:01:43 +0000
@@ -2574,11 +2574,14 @@
       PARAMETER (NEXTERNAL=4)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=1)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/MG5_1_USERCHOICE/USERHEL
 C     ----------
@@ -2602,7 +2605,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2719,7 +2725,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/dux_mumvmxg/loop_matrix.f'
--- tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/dux_mumvmxg/loop_matrix.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/dux_mumvmxg/loop_matrix.f	2020-08-13 08:01:43 +0000
@@ -2548,11 +2548,14 @@
       PARAMETER (NEXTERNAL=5)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=0)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/ML5_0_USERCHOICE/USERHEL
 C     ----------
@@ -2576,7 +2579,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2693,7 +2699,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/gg_wmtbx/loop_matrix.f'
--- tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/gg_wmtbx/loop_matrix.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/long_ML_SMQCD_optimized/gg_wmtbx/loop_matrix.f	2020-08-13 08:01:43 +0000
@@ -2548,11 +2548,14 @@
       PARAMETER (NEXTERNAL=5)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=0)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/ML5_0_USERCHOICE/USERHEL
 C     ----------
@@ -2576,7 +2579,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2693,7 +2699,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/ddx_ttx/loop_matrix.f'
--- tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/ddx_ttx/loop_matrix.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/ddx_ttx/loop_matrix.f	2020-08-13 08:01:43 +0000
@@ -2548,11 +2548,14 @@
       PARAMETER (NEXTERNAL=4)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=0)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/ML5_0_USERCHOICE/USERHEL
 C     ----------
@@ -2576,7 +2579,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2693,7 +2699,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/gg_ttx/loop_matrix.f'
--- tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/gg_ttx/loop_matrix.f	2020-06-20 21:10:10 +0000
+++ tests/input_files/IOTestsComparison/short_ML_SMQCD_optimized/gg_ttx/loop_matrix.f	2020-08-13 08:01:43 +0000
@@ -2548,11 +2548,14 @@
       PARAMETER (NEXTERNAL=4)
       INTEGER    NSQUAREDSO
       PARAMETER (NSQUAREDSO=0)
+      INCLUDE 'nsqso_born.inc'
 C     
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, USERHEL
       COMMON/ML5_0_USERCHOICE/USERHEL
 C     ----------
@@ -2576,7 +2579,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       INTEGER HEL, RET_CODE
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
 C     
@@ -2693,7 +2699,10 @@
 C     ARGUMENTS 
 C     
       REAL*8 P(0:3,NEXTERNAL)
-      REAL*8 ANS(0:3,0:NSQUAREDSO)
+      INCLUDE 'nsqso_born.inc'
+      INTEGER ANS_DIMENSION
+      PARAMETER(ANS_DIMENSION=MAX(NSQSO_BORN,NSQUAREDSO))
+      REAL*8 ANS(0:3,0:ANS_DIMENSION)
       REAL*8 PREC_ASKED,PREC_FOUND(0:NSQUAREDSO)
       INTEGER RET_CODE
 C     
 
=== modified file 'tests/test_manager.py'
--- tests/test_manager.py	2018-08-10 12:40:47 +0000
+++ tests/test_manager.py	2020-08-13 08:01:43 +0000
@@ -44,6 +44,7 @@
 import time
 import datetime
 import shutil
+import subprocess
 import glob
 from functools import wraps
 
@@ -973,7 +974,8 @@
 
     parser.add_option("-N", "--notification", default=45,
           help="Running time, below which no notification is raised. (-1 for no notification)")        
-
+    parser.add_option("", "--nocaffeinate", action="store_false", default=True, dest='nosleep',
+                  help='For mac user, forbids to use caffeinate when running with a script')
     
     (options, args) = parser.parse_args()
 
@@ -1046,6 +1048,14 @@
         logging.getLogger('tutorial').setLevel('ERROR')
     except:
         pass
+
+    if sys.platform == "darwin" and options.nosleep:
+        logging.getLogger('madgraph').warning("launching caffeinate to prevent idle sleep when MG5aMC is running. Run './bin/mg5_aMC -s' to prevent this.")
+        pid = os.getpid()
+        subprocess.Popen(['caffeinate', '-i', '-w', str(pid)])
+        
+
+
     
     if options.synchronize and IOTestManager._compress_ref_fodler:
         print "The tarball synchronization is not necessary since"+ \