Merge lp:~maddevelopers/mg5amcnlo/2.6.5_ewa into lp:~mg5core1/mg5amcnlo/2.7.4

[Olivier Mattelaer]

[Olivier Mattelaer (olivier-mattelaer)]

put comment for richard

[Olivier Mattelaer (olivier-mattelaer)]

so the important changes are
1) forcing MC over helicity [here simple copy-paste]
2) moving polarization from the parton to the PDF section (when needed--not perfectly done here due to DIS lack of support)
3) include the new PDF [here you should be ok, be carefull with the POL() normalization]
4) will need to worry about the number of helicty allowed for the W (what if only transverse)

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

=== added file 'Template/LO/Source/PDF/ElectroweakFlux.f'
--- Template/LO/Source/PDF/ElectroweakFlux.f	1970-01-01 00:00:00 +0000
+++ Template/LO/Source/PDF/ElectroweakFlux.f	2021-02-05 11:45:49 +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	2021-02-05 11:45:49 +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-08-07 14:49:25 +0000
+++ Template/LO/Source/PDF/PhotonFlux.f	2021-02-05 11:45:49 +0000
@@ -1,8 +1,10 @@
 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
 
 
=== modified file 'Template/LO/Source/PDF/makefile'
--- Template/LO/Source/PDF/makefile	2020-04-08 19:22:14 +0000
+++ Template/LO/Source/PDF/makefile	2021-02-05 11:45:49 +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/pdg2pdf.f'
--- Template/LO/Source/PDF/pdg2pdf.f	2020-08-07 14:49:25 +0000
+++ Template/LO/Source/PDF/pdg2pdf.f	2021-02-05 11:45:49 +0000
@@ -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(abs(ih).eq.3.or.abs(ih).eq.4) then       !from the electron or muonn
-            pdg2pdf=epa_lepton(x,q2max, ih)
-         elseif(ih .eq. 2) then !from a proton without breaking
-            pdg2pdf=epa_proton(x,q2max,beamid)
-         endif 
-         pdflast(iporg,ireuse)=pdg2pdf
-         return
+         if (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
 
=== modified file 'Template/LO/Source/PDF/pdg2pdf_lhapdf6.f'
--- Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	2020-08-07 14:49:25 +0000
+++ Template/LO/Source/PDF/pdg2pdf_lhapdf6.f	2021-02-05 11:45:49 +0000
@@ -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(abs(ih).eq.3.or.abs(ih).eq.4) then       !from the electron
-            pdg2pdf=epa_lepton(x,xmu*xmu, ih)
-      else if(ih.eq.2) then ! photon from a proton without breaking
-          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/setrun.f'
--- Template/LO/Source/setrun.f	2020-08-07 14:49:25 +0000
+++ Template/LO/Source/setrun.f	2021-02-05 11:45:49 +0000
@@ -97,14 +97,13 @@
          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-08-10 19:06:52 +0000
+++ Template/LO/SubProcesses/genps.f	2021-02-05 11:45:49 +0000
@@ -572,6 +572,7 @@
             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)
 
=== modified file 'Template/MadWeight/src/setrun.f'
--- Template/MadWeight/src/setrun.f	2018-04-27 21:01:37 +0000
+++ Template/MadWeight/src/setrun.f	2021-02-05 11:45:49 +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 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2020-08-20 15:12:34 +0000
+++ madgraph/core/base_objects.py	2021-02-05 11:45:49 +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'] = {}
@@ -4068,7 +4069,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-08-17 22:31:42 +0000
+++ madgraph/interface/common_run_interface.py	2021-02-05 11:45:49 +0000
@@ -5839,6 +5839,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	2021-02-05 11:45:49 +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/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2020-08-20 13:17:30 +0000
+++ madgraph/iolibs/export_v4.py	2021-02-05 11:45:49 +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)    
@@ -4197,11 +4234,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
@@ -4308,14 +4358,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	2021-02-05 11:45:49 +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/matrix_madevent_group_v4.inc'
--- madgraph/iolibs/template_files/matrix_madevent_group_v4.inc	2020-05-18 19:26:16 +0000
+++ madgraph/iolibs/template_files/matrix_madevent_group_v4.inc	2021-02-05 11:45:49 +0000
@@ -127,13 +127,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
@@ -178,14 +172,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	2021-02-05 11:45:49 +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	2021-02-05 11:45:49 +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/various/banner.py'
--- madgraph/various/banner.py	2020-08-21 07:00:01 +0000
+++ madgraph/various/banner.py	2021-02-05 11:45:49 +0000
@@ -2897,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,4,-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 from electron, 4: photon from muon, 9: PLUGIN MODE')
-        self.add_param("lpp2", 1, fortran_name="lpp(2)", allowed=[-1,1,0,2,3,9,4,-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 from electron, 4: photon from muon, 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,
@@ -3296,8 +3296,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']:
@@ -3321,11 +3324,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')
@@ -3339,7 +3418,7 @@
                         self['lpp2'] = 0  
                         self['ebeam1'] = '6500'  
                         self['ebeam2'] = '1k'  
-            
+
             elif any(id in beam_id for id in [11,-11,13,-13]):
                 self['lpp1'] = 0
                 self['lpp2'] = 0
@@ -3349,6 +3428,9 @@
                 if set([ abs(i) for i in beam_id_split[0]]) == set([ abs(i) for i in beam_id_split[1]]):
                     self.display_block.append('ecut')
                 self.display_block.append('beam_pol')
+
+                if beam_id_split[0] == beam_id_split[1]:
+                    self.display_block.append('ecut')
             else:
                 self['lpp1'] = 0
                 self['lpp2'] = 0    
 
=== modified file 'madgraph/various/cluster.py'
--- madgraph/various/cluster.py	2017-11-17 14:15:01 +0000
+++ madgraph/various/cluster.py	2021-02-05 11:45:49 +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 '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	2021-02-05 11:45:49 +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-08-20 15:12:34 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/matrix1.f	2021-02-05 11:45:49 +0000
@@ -154,6 +154,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
@@ -162,6 +163,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)
@@ -220,6 +222,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	2021-02-05 11:45:49 +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	2021-02-05 11:45:49 +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	2021-02-05 11:45:49 +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