Master for360

epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/MG5aMC_patches/PROD/patch.P1

@@ -126,34 +126,34 @@ index 1124a9164..27a6e4674 100644
        open(unit=lun,file=tempname,status='old',ERR=20)
        fopened=.true.
 diff --git b/epochX/cudacpp/gg_tt.mad/SubProcesses/P1_gg_ttx/matrix1.f a/epochX/cudacpp/gg_tt.mad/SubProcesses/P1_gg_ttx/matrix1.f
-index 8b73188ef..710fd1297 100644
+index 1acba8200..069c74ef4 100644
 --- b/epochX/cudacpp/gg_tt.mad/SubProcesses/P1_gg_ttx/matrix1.f
 +++ a/epochX/cudacpp/gg_tt.mad/SubProcesses/P1_gg_ttx/matrix1.f
-@@ -72,7 +72,10 @@ C
+@@ -71,7 +71,10 @@ C
        DATA NB_FAIL /0/
        DOUBLE PRECISION GET_CHANNEL_CUT
        EXTERNAL GET_CHANNEL_CUT
 -
 +C
-+      INTEGER NGOODHEL(2) ! -1 if not yet retrieved and printed
++      INTEGER NGOODHEL ! -1 if not yet retrieved and printed
 +      SAVE NGOODHEL
-+      DATA NGOODHEL/-1,-1/
++      DATA NGOODHEL/-1/
  C     
  C     This is just to temporarily store the reference grid for
  C      helicity of the DiscreteSampler so as to obtain its number of
-@@ -230,6 +233,17 @@ C            update.
+@@ -211,6 +214,17 @@ C     ----------
            ENDIF
-           IF(NTRY(IMIRROR).EQ.MAXTRIES)THEN
-             ISHEL(IMIRROR)=MIN(ISUM_HEL,NGOOD(IMIRROR))
+           IF(NTRY(1).EQ.MAXTRIES)THEN
+             ISHEL=MIN(ISUM_HEL,NGOOD)
 +C           Print the number of good helicities
-+            IF (NGOODHEL(IMIRROR).EQ.-1) THEN
-+              NGOODHEL(IMIRROR)=0
++            IF (NGOODHEL.EQ.-1) THEN
++              NGOODHEL=0
 +              DO I=1,NCOMB
-+                IF (GOODHEL(I,IMIRROR)) THEN
-+                  NGOODHEL(IMIRROR)=NGOODHEL(IMIRROR)+1
++                IF (GOODHEL(I,1)) THEN
++                  NGOODHEL=NGOODHEL+1
 +                ENDIF
 +              END DO
-+              WRITE (6,*) 'NGOODHEL =', NGOODHEL(IMIRROR) ! no need to print imirror?
++              WRITE (6,*) 'NGOODHEL =', NGOODHEL
 +              WRITE (6,*) 'NCOMB =', NCOMB
 +            ENDIF
            ENDIF

epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/MG5aMC_patches/PROD/patch.common

@@ -1,17 +1,3 @@
-diff --git b/epochX/cudacpp/gg_tt.mad/Source/genps.inc a/epochX/cudacpp/gg_tt.mad/Source/genps.inc
-index a59181c70..af7e0efbc 100644
---- b/epochX/cudacpp/gg_tt.mad/Source/genps.inc
-+++ a/epochX/cudacpp/gg_tt.mad/Source/genps.inc
-@@ -30,7 +30,8 @@ c*************************************************************************
- c	Parameters for helicity sums in matrixN.f
- c*************************************************************************
-       REAL*8 LIMHEL
--      PARAMETER(LIMHEL=1e-8) ! ME threshold for helicity filtering
-+c     PARAMETER(LIMHEL=1e-8) ! ME threshold for helicity filtering (Fortran default)
-+      PARAMETER(LIMHEL=0) ! ME threshold for helicity filtering (force Fortran to mimic cudacpp, see #419)
-       INTEGER MAXTRIES
-       PARAMETER(MAXTRIES=25)
- C     To pass the helicity configuration chosen by the DiscreteSampler to 
 diff --git b/epochX/cudacpp/gg_tt.mad/SubProcesses/makefile a/epochX/cudacpp/gg_tt.mad/SubProcesses/makefile
 index 348c283be..49e6800ff 100644
 --- b/epochX/cudacpp/gg_tt.mad/SubProcesses/makefile

epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/launch_plugin.py

@@ -98,6 +98,7 @@ def default_setup(self):
         self['vector_size'] = 16 # already setup in default class (just change value)
         self['aloha_flag'] = '--fast-math'
         self['matrix_flag'] = '-O3'
+        self['limhel'] = 0
         self.display_block.append('simd')
         self.display_block.append('psoptim')
 

epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/madgraph/iolibs/template_files/gpu/smatrix_multi.f

@@ -22,9 +22,9 @@
      &      SELECTED_HEL2, SELECTED_COL2, .TRUE.)  ! quit after computing helicities
           FIRST = .FALSE.
 c         ! This is a workaround for https://github.com/oliviermattelaer/mg5amc_test/issues/22 (see PR #486)
-          IF( FBRIDGE_MODE .EQ. 1 ) THEN  ! (CppOnly=1 : SMATRIX1 is not called at all)
-            CALL RESET_CUMULATIVE_VARIABLE()  ! mimic 'avoid bias of the initialization' within SMATRIX1
-          ENDIF
+c          IF( FBRIDGE_MODE .EQ. 1 ) THEN  ! (CppOnly=1 : SMATRIX1 is not called at all)
+c            CALL RESET_CUMULATIVE_VARIABLE()  ! mimic 'avoid bias of the initialization' within SMATRIX1
+c          ENDIF
           CALL FBRIDGEGETNGOODHEL(FBRIDGE_PBRIDGE,NGOODHEL,NTOTHEL)
           IF( NTOTHEL .NE. NCOMB ) THEN
             WRITE(6,*) 'ERROR! Cudacpp/Fortran mismatch',

epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/output.py

@@ -195,7 +195,7 @@ def copy_template(self, model):
                 open(os.path.join('test', 'cudacpp_test.mk'), 'w').write(makefile_test)
 
     # OM - overload export_v4.py version to add additional_clean section (and avoid patchMad.sh for Source/makefile)
-    def write_source_makefile(self, writer, default=None):
+    def write_source_makefile(self, writer, model=None, default=None):
         if default:
             replace_dict = default
         else:

epochX/cudacpp/ee_mumu.mad/Source/MODEL/couplings3.f

@@ -0,0 +1,17 @@
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c      written by the UFO converter
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      SUBROUTINE COUP3( VECID)
+
+      IMPLICIT NONE
+      INTEGER VECID
+      INCLUDE 'model_functions.inc'
+
+      DOUBLE PRECISION PI, ZERO
+      PARAMETER  (PI=3.141592653589793D0)
+      PARAMETER  (ZERO=0D0)
+      INCLUDE 'input.inc'
+      INCLUDE '../vector.inc'
+      INCLUDE 'coupl.inc'
+      END

epochX/cudacpp/ee_mumu.mad/Source/PDF/lep_densities/isrbetll/eepdf.f

@@ -0,0 +1,113 @@
+c     This function calculate the reduced structure function, with energy
+c     fraction given by "x", at scale "Qsquare"
+      real*8 function eepdf_tilde(x,Qsquare,n,partonid,beamid)
+      implicit none
+      real*8 x,Qsquare
+      real*8 me
+      data me /0.511d-3/
+      real*8 PI
+      real*8 alphaem
+c     In Gmu scheme
+      data alphaem/0.007562397d0/
+      real*8 beta
+      integer n,partonid,beamid
+      real*8 isr_tilde_racoon
+      real*8 res
+      data res/0d0/
+
+      PI=4.D0*DATAN(1.D0)
+
+      beta = alphaem/PI * (dlog(Qsquare/me/me)-1)
+
+c     electron beam
+      if (beamid .eq. 11) then
+c     other partons are zero
+          if (partonid .ne. 11) then
+             res = 0d0
+          else
+             if (n .eq. 1) then
+                res = isr_tilde_racoon(x,beta)
+             else
+                res = 0d0
+             endif
+          endif
+      else if (beamid .eq. -11) then
+          if (partonid .ne. -11) then
+              res = 0d0
+          else
+             if (n .eq. 1) then
+                res = isr_tilde_racoon(x,beta)
+             else
+                res = 0d0
+             endif
+          endif
+      endif
+      eepdf_tilde = res
+      end
+
+c     https://arxiv.org/pdf/hep-ph/0302198.pdf, eq.(2.44)
+c     note that beta_e in eq.(2.45) is twice our beta
+c     so eq.(2.44) needs to be corrected by some factor of 2
+      real*8 function isr_tilde_racoon(x,beta)
+      implicit none
+      real*8 x,beta
+      real*8 res
+      real*8 PI
+      real*8 gE
+      real*8 logx, logomx
+      real*8 dlgam,DDILOG
+      external dlgam,DDILOG
+      PI=4.D0*DATAN(1.D0)
+      gE=0.5772156649d0
+      res=0d0
+      if (x .lt. 0.9999999d0) then
+         logx=dlog(x)
+         logomx=dlog(1d0-x)
+c     ----------------------------
+c     order alpha
+         res=-beta*(1d0+x)/2d0
+c     order alpha^2
+         res=res-(beta**2)/8d0*(
+     c        (1d0+3d0*x*x)/(1d0-x)*logx
+     c        +4d0*(1d0+x)*logomx+5d0+x)
+c     order alpha^3
+         res=res-(beta**3)/128d0*(
+     c        (1d0+x)*(6d0*DDILOG(x)+12d0*(logomx**2)-3d0*PI**2)
+     c        +(1d0/(1d0-x)) * (
+     c        (3d0/2d0)*(1d0+8d0*x+3d0*x**2)*logx
+     c        +6d0*(x+5d0)*(1d0-x)*logomx
+     c        +12d0*(1d0+x**2)*logx*logomx
+     c        -1d0/2d0*(1d0+7d0*x**2)*logx**2
+     c        +1d0/4d0*(39d0-24d0*x-15d0*x**2)))
+c     ----------------------------
+         res=res*(1d0-x)**(1-beta)
+      endif
+      res=res+exp(beta*(-gE+3d0/4d0))/exp(dlgam(1d0+beta))*beta
+      isr_tilde_racoon=res
+      end
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
+      function dlgam(xx)
+c real logarithm of gamma function
+      implicit real * 8 (a-h,o-z)
+      real * 8 cof(6),stp,gt,g,cst
+      data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,
+     # -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
+      data gt,g/5.5d0,5.0d0/
+      data cst/4.081061466d0/
+      x = xx - 1
+      xpgt = x + gt
+      xmp5  = x + .5d0
+      s0 = 1
+      do 1 j=1,6
+        x = x + 1
+        tmp = cof(j)/x
+        s0  = s0 + tmp
+  1     continue
+      r10 = log(s0)
+      dlgam = xmp5*(log(xpgt)-1) + r10 - cst
+      return
+      end      
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+

epochX/cudacpp/ee_mumu.mad/Source/PDF/lep_densities/isrbetll/gridpdfaux.f

@@ -0,0 +1,140 @@
+      integer function eepdf_n_components(partonid,beamid)
+      implicit none
+      integer partonid,beamid
+      integer ncom
+c     electron beam
+      if (beamid .eq. 11) then
+c     other partons are zero
+        if (partonid .ne. 11) then
+          ncom=0
+        else
+          ncom=1
+        endif
+      else if (beamid .eq. -11) then
+        if (partonid .ne. -11) then
+          ncom=0
+        else
+          ncom=1
+        endif
+      endif
+      eepdf_n_components=ncom
+      end
+
+c     This function return the power of (1-x)
+      real*8 function eepdf_tilde_power(Q2,n,partonid,beamid)
+      implicit none
+      real*8 me
+      data me /0.511d-3/
+      real*8 PI
+      real*8 alphaem
+c     In Gmu scheme
+      data alphaem/0.007562397d0/
+      real*8 beta,eta,Q2
+      integer n,partonid,beamid
+      real*8 k,b
+
+      PI=4.D0*DATAN(1.D0)
+      beta = alphaem/PI * (dlog(Q2/me/me)-1)
+      eta = alphaem/PI * (dlog(Q2/me/me))
+      b=-2.D0/3.D0
+
+c     electron beam
+      if (beamid .eq. 11) then
+c     other partons are zero
+        if (partonid .ne. 11) then
+          k=0d0
+        else
+          if (n .eq. 1) then
+            k=1d0-beta
+          else
+            k=0d0
+          endif
+        endif
+      else if (beamid .eq. -11) then
+        if (partonid .ne. -11) then
+          k=0d0
+        else
+          if (n .eq. 1) then
+            k=1d0-beta
+          else
+            k=0d0
+          endif
+        endif
+      endif
+      eepdf_tilde_power = k
+      end
+
+c     This function return the type of this component
+      integer function eepdf_tilde_type(n,partonid,beamid)
+      implicit none
+      integer n,partonid,beamid
+      integer res
+
+c     electron beam
+      if (beamid .eq. 11) then
+c     other partons are zero
+        if (partonid .ne. 11) then
+          res=0
+        else
+          if (n .eq. 1) then
+            res=1
+          else
+            res=0
+          endif
+        endif
+      else if (beamid .eq. -11) then
+        if (partonid .ne. -11) then
+          res=0
+        else
+          if (n .eq. 1) then
+            res=1
+          else
+            res=0
+          endif
+        endif
+      endif
+      eepdf_tilde_type = res
+      end
+
+c     This is to calculate the factor for grid implementation
+      real*8 function eepdf_tilde_factor(x,Q2,n,partonid,beamid)
+      implicit none
+      real*8 x,Q2
+      real*8 me
+      data me /0.511d-3/
+      real*8 PI
+      real*8 alphaem
+c     In Gmu scheme
+      data alphaem/0.007562397d0/
+      real*8 beta
+      integer n,partonid,beamid
+      real*8 res
+
+      PI=4.D0*DATAN(1.D0)
+      beta = alphaem/PI * (dlog(Q2/me/me))
+
+c     electron beam
+      if (beamid .eq. 11) then
+c     other partons are zero
+        if (partonid .ne. 11) then
+          res=1d0
+        else
+          if (n .eq. 1) then
+            res = 1d0
+          else
+            res = 1d0
+          endif
+        endif
+      else if (beamid .eq. -11) then
+        if (partonid .ne. -11) then
+          res = 1d0
+        else
+          if (n .eq. 1) then
+            res = 1d0
+          else
+            res = 1d0
+          endif
+        endif
+      endif
+      eepdf_tilde_factor = res
+      end