Merge lp:~mg5core1/mg5amcnlo/2.6.4 into lp:mg5amcnlo/lts

[Olivier Mattelaer]

2.6.4 is now ready. I will procede to merging it.

=== modified file 'MadSpin/decay.py'
--- MadSpin/decay.py	2018-03-22 14:00:54 +0000
+++ MadSpin/decay.py	2018-11-09 10:32:12 +0000
@@ -3963,14 +3963,23 @@
                 elif mode=='full':
                     stdin_text="5 0 0 0 \n"  # before closing, write down the seed 
                     external = self.calculator[('full',path)]
-                    external.stdin.write(stdin_text)
+                    try:
+                        external.stdin.write(stdin_text)
+                    except Exception:
+                        continue
                     ranmar_state=external.stdout.readline()
                     ranmar_file=pjoin(path,'ranmar_state.dat')
                     ranmar=open(ranmar_file, 'w')
                     ranmar.write(ranmar_state)
                     ranmar.close()
-                    external.stdin.close()
-                    external.stdout.close()
+                    try:
+                        external.stdin.close()
+                    except Exception:
+                        continue
+                    try:
+                        external.stdout.close()
+                    except Exception:
+                        continue
                     external.terminate()
                     del external
         else:
 
=== modified file 'MadSpin/interface_madspin.py'
--- MadSpin/interface_madspin.py	2018-06-22 06:52:02 +0000
+++ MadSpin/interface_madspin.py	2018-11-09 10:32:12 +0000
@@ -199,10 +199,11 @@
 
         self.inputfile = inputfile
         if self.options['spinmode'] == 'none' and \
-           (self.options['input_format'] != 'lhe' or (self.options['input_format'] == 'auto' and '.lhe' in inputfile[:-5])):  
+           (self.options['input_format'] not in ['lhe','auto'] or 
+             (self.options['input_format'] == 'auto' and '.lhe'  not in inputfile[-7:])):  
             self.banner = banner.Banner()
             self.setup_for_pure_decay()
-            return   
+            return  
         
         if inputfile.endswith('.gz'):
             misc.gunzip(inputfile)
@@ -984,8 +985,11 @@
                 #misc.sprint(i, particle.pdg, particle.pid)
                 #misc.sprint(self.final_state, evt_decayfile)
                 # check if we need to decay the particle 
-                if not (particle.pdg in self.final_state or particle.pdg in evt_decayfile):
-                    continue # nothing to do for this particle
+                if self.final_state and particle.pdg not in self.final_state:
+                    continue # nothing to do for this particle
+                if particle.pdg not in evt_decayfile:
+                    continue # nothing to do for this particle
+                
                 # check how the decay need to be done
                 nb_decay = len(evt_decayfile[particle.pdg])
                 if nb_decay == 0:
@@ -1004,11 +1008,18 @@
                     cumul = 0
                     for j,events in evt_decayfile[particle.pdg].items():
                         cumul += events.cross
-                        if r < cumul:
+                        if r <= cumul:
                             decay_file = events
                             decay_file_nb = j
-                        else:
                             break
+                    else:
+                        # security for numerical accuracy issue... (unlikely but better safe)
+                        if (cumul-tot)/tot < 1e-5:
+                            decay_file = events
+                            decay_file_nb = j
+                        else:
+                            misc.sprint(j,cumul, events.cross, tot, (tot-cumul)/tot)
+                            raise Exception
                 
                 if self.options['new_wgt'] == 'BR':
                     tot_width = float(self.banner.get('param', 'decay', abs(pdg)).value)
@@ -1056,8 +1067,13 @@
                     if len(hepmc_output) == 0:
                         hepmc_output.append(lhe_parser.Particle(event=hepmc_output))
                         hepmc_output[0].color2 = 0
+                        hepmc_output[0].status = -1
+                        hepmc_output.nexternal+=1
                     decayed_particle = lhe_parser.Particle(particle, hepmc_output)
+                    decayed_particle.mother1 = hepmc_output[0]
+                    decayed_particle.mother2 = hepmc_output[0]
                     hepmc_output.append(decayed_particle)
+                    hepmc_output.nexternal+=1
                     decayed_particle.add_decay(decay)
             # change the weight associate to the event
             if self.options['new_wgt'] == 'cross-section':
@@ -1075,7 +1091,6 @@
             else:
                 hepmc_output.wgt = event.wgt
                 hepmc_output.nexternal = len(hepmc_output) # the append does not update nexternal
-                hepmc_output.assign_mother()
                 output_lhe.write(str(hepmc_output))
         else:
             if counter==0:
@@ -1425,7 +1440,7 @@
             if self.options['fixed_order']:
                 production, counterevt= production[0], production[1:]
             if curr_event and curr_event % 1000 == 0 and float(str(curr_event)[1:]) ==0:
-                print "decaying event number %s. Efficiency: %s [%s s]" % (curr_event, 1/self.efficiency, time.time()-start)
+                logger.info("decaying event number %s. Efficiency: %s [%s s]" % (curr_event, 1/self.efficiency, time.time()-start))
             while 1:
                 nb_try +=1
                 decays = self.get_decay_from_file(production, evt_decayfile, nb_event-curr_event)
 
=== modified file 'MadSpin/src/driver.f'
--- MadSpin/src/driver.f	2017-08-03 18:49:25 +0000
+++ MadSpin/src/driver.f	2018-11-09 10:32:12 +0000
@@ -340,10 +340,10 @@
            if (jac.lt.0d0) then
              counter2=counter2+1 
              counter3=counter3+1 
-c             if (counter3.gt.500) then
-c               write(*,*) "500_pts_failed_stop_executation"
-c               stop
-c             endif
+             if (counter3.gt.500*8*100) then
+               write(*,*) "500_pts_failed_stop_executation"
+               stop
+             endif
              if (counter2.ge.8) then ! use another topology to generate PS points
                do k=1,n_max_cg
                  amp2(k)=0d0
@@ -960,7 +960,7 @@
       include 'nexternal_prod.inc'
       !include 'run.inc'
 c arguments
-      double precision jac,x(36),p(0:3,nexternal), p_prod(0:3,nexternal)
+      double precision jac,x(36),p(0:3,nexternal), p_prod(0:3,nexternal_prod)
       integer itree(2,-nexternal:-1)
       double precision qmass(-nexternal:0),qwidth(-nexternal:0)
 c common
 
=== modified file 'Template/LO/Source/run.inc'
--- Template/LO/Source/run.inc	2018-02-01 22:30:56 +0000
+++ Template/LO/Source/run.inc	2018-11-09 10:32:12 +0000
@@ -91,4 +91,7 @@
       double precision mxxmin4pdg(0:25) 
       logical mxxpart_antipart(1:25)
       common/TO_PDG_SPECIFIC_CUT/pdg_cut, ptmin4pdg,ptmax4pdg, Emin4pdg, Emax4pdg, etamin4pdg,
-     &etamax4pdg, mxxmin4pdg,mxxpart_antipart
 \ No newline at end of file
+     &etamax4pdg, mxxmin4pdg,mxxpart_antipart
+
+       double precision small_width_treatment
+       common/narrow_width/small_width_treatment
 \ No newline at end of file
 
=== modified file 'Template/LO/Source/transpole.f'
--- Template/LO/Source/transpole.f	2016-06-15 22:46:42 +0000
+++ Template/LO/Source/transpole.f	2018-11-09 10:32:12 +0000
@@ -28,13 +28,24 @@
       double precision z,zmin,zmax,xmin,xmax,ez
       double precision pole1,width1,x,xc
       double precision a,b
+c
+c     small width treatment
+c
+      double precision small_width_treatment
+      common/narrow_width/small_width_treatment
 c-----
 c  Begin Code
 c-----
       pole=pole1
       width=width1
+
       x = x1
       if (pole .gt. 0d0) then
+         if (width.lt.pole*small_width_treatment)then
+            width = pole * small_width_treatment
+            jac = jac * width/width1
+         endif
+
          zmin = atan((-pole)/width)/width
          zmax = atan((1d0-pole)/width)/width
          if (x .gt. del .and. x .lt. 1d0-del) then
@@ -175,7 +186,11 @@
 c
       double precision pole1,width1,y1,jac
       real*8 x
-
+c
+c     small width treatment
+c
+      double precision small_width_treatment
+      common/narrow_width/small_width_treatment
 c
 c     Local
 c
@@ -191,6 +206,10 @@
       width=width1
       y = y1
       if (pole .gt. 0d0) then                   !BW 
+         if (width.lt.pole*small_width_treatment)then
+            width = pole * small_width_treatment
+            jac = jac * width/width1
+         endif
          zmin = atan((-pole)/width)/width
          zmax = atan((1d0-pole)/width)/width
          z = atan((y-pole)/width)/width
 
=== modified file 'Template/LO/SubProcesses/genps.f'
--- Template/LO/SubProcesses/genps.f	2018-06-20 11:54:42 +0000
+++ Template/LO/SubProcesses/genps.f	2018-11-09 10:32:12 +0000
@@ -41,7 +41,7 @@
       integer mincfig,maxcfig                  !Range of configurations
       integer invar
       double precision wgt                     !(input and output)
-      double precision x(maxdim),p(maxdim)     !x,p (output) [p(0:3,nexternal)]
+      double precision x(*),p(*)     !x,p (output) [p(0:3,nexternal)]
 c
 c     Local
 c
@@ -568,7 +568,7 @@
             if (abs(lpp(1)) .eq. 3) m1 = 0.000511d0
             if (abs(lpp(2)) .eq. 3) m2 = 0.000511d0
             if (mass_ion(1).ge.0d0) m1 = mass_ion(1)
-            if (mass_ion(2).ge.0d0) m1 = mass_ion(2)
+            if (mass_ion(2).ge.0d0) m2 = mass_ion(2)
             if(ebeam(1).lt.m1.and.lpp(1).ne.9) ebeam(1)=m1
             if(ebeam(2).lt.m2.and.lpp(2).ne.9) ebeam(2)=m2
             pi1(0)=ebeam(1)
 
=== modified file 'Template/LO/SubProcesses/myamp.f'
--- Template/LO/SubProcesses/myamp.f	2016-12-17 16:26:56 +0000
+++ Template/LO/SubProcesses/myamp.f	2018-11-09 10:32:12 +0000
@@ -1,91 +1,3 @@
-      double precision function testamp(p)
-c*****************************************************************************
-c     Approximates matrix element by propagators
-c*****************************************************************************
-      implicit none
-c
-c     Constants
-c     
-      include 'genps.inc'
-      include 'maxconfigs.inc'
-      include 'nexternal.inc'
-      double precision   zero
-      parameter (zero = 0d0)
-c
-c     Arguments
-c
-      double precision p(0:3,nexternal)
-c      integer iconfig
-c
-c     Local
-c
-      double precision xp(0:3,-nexternal:nexternal)
-      double precision mpole(-nexternal:0),shat,tsgn
-      integer i,j,iconfig
-
-      double precision prmass(-nexternal:0,lmaxconfigs)
-      double precision prwidth(-nexternal:0,lmaxconfigs)
-      integer pow(-nexternal:0,lmaxconfigs)
-      logical first_time
-c
-c     Global
-c
-      integer iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
-      integer            mapconfig(0:lmaxconfigs), this_config
-      common/to_mconfigs/mapconfig, this_config
-      
-      include 'coupl.inc'
-c
-c     External
-c
-      double precision dot
-
-      save prmass,prwidth,pow
-      data first_time /.true./
-c-----
-c  Begin Code
-c-----      
-      iconfig = this_config
-      if (first_time) then
-c         include 'props.inc'
-         first_time=.false.
-      endif
-
-      do i=1,nexternal
-         mpole(-i)=0d0
-         do j=0,3
-            xp(j,i)=p(j,i)
-         enddo
-      enddo
-c      mpole(-3) = 174**2
-c      shat = dot(p(0,1),p(0,2))/(1800)**2
-      shat = dot(p(0,1),p(0,2))/(10)**2
-c      shat = 1d0
-      testamp = 1d0
-      tsgn    = +1d0
-      do i=-1,-(nexternal-3),-1              !Find all the propagotors
-         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
-         do j=0,3
-            xp(j,i) = xp(j,iforest(1,i,iconfig))
-     $           +tsgn*xp(j,iforest(2,i,iconfig))
-         enddo
-         if (prwidth(i,iconfig) .ne. 0d0 .and. .false.) then
-            testamp=testamp/((dot(xp(0,i),xp(0,i))
-     $                        -prmass(i,iconfig)**2)**2
-     $         -(prmass(i,iconfig)*prwidth(i,iconfig))**2)
-         else
-            testamp = testamp/((dot(xp(0,i),xp(0,i)) -
-     $                          prmass(i,iconfig)**2)
-     $                          **(pow(i,iconfig)))
-         endif
-        testamp=testamp*shat**(pow(i,iconfig))
-c        write(*,*) i,iconfig,pow(i,iconfig),prmass(i,iconfig)
-      enddo
-c      testamp = 1d0/dot(xp(0,-1),xp(0,-1))
-      testamp=abs(testamp)
-c      testamp = 1d0
-      end
 
       logical function cut_bw(p)
 c*****************************************************************************
@@ -114,6 +26,7 @@
 
       double precision prmass(-nexternal:0,lmaxconfigs)
       double precision prwidth(-nexternal:0,lmaxconfigs)
+      double precision prwidth_tmp(-nexternal:0,lmaxconfigs)
       integer pow(-nexternal:0,lmaxconfigs)
       logical first_time, onshell
       double precision xmass
@@ -152,13 +65,14 @@
 c
       double precision dot
 
-      save prmass,prwidth,pow
+      save prmass,prwidth,pow,prwidth_tmp
       data first_time /.true./
 c-----
 c  Begin Code
 c-----      
       cut_bw = .false.    !Default is we passed the cut
       iconfig = this_config
+
       if (first_time) then
          include 'props.inc'
          nbw = 0
@@ -209,10 +123,12 @@
 c            write(*,*) 'xmass',xmass,prmass(i,iconfig)
 c
 c           Here we set if the BW is "on-shell" for LesHouches
-c
+c            
+            prwidth_tmp(i,iconfig) = max(prwidth(i,iconfig), prmass(i,iconfig)*small_width_treatment)
+ 
             onshell = (abs(xmass - prmass(i,iconfig)) .lt.
-     $           bwcutoff*prwidth(i,iconfig).and.
-     $           (prwidth(i,iconfig)/prmass(i,iconfig).lt.0.1d0.or.
+     $           bwcutoff*prwidth_tmp(i,iconfig).and.
+     $           (prwidth_tmp(i,iconfig)/prmass(i,iconfig).lt.0.1d0.or.
      $            gForceBW(i,iconfig).eq.1))
             if(onshell)then
 c     Remove on-shell forbidden s-channels (gForceBW=2) (JA 2/10/11)
@@ -264,10 +180,10 @@
 c     For decay-chain syntax use BWcutoff here too (22/12/14)
             if (gForceBW(i, iconfig).eq.1) then
                onshell = (abs(xmass - prmass(i,iconfig)) .lt.
-     $           bwcutoff*prwidth(i,iconfig))
+     $           bwcutoff*prwidth_tmp(i,iconfig))
             else
                onshell = (abs(xmass - prmass(i,iconfig)) .lt.
-     $           5d0*prwidth(i,iconfig))
+     $           5d0*prwidth_tmp(i,iconfig))
             endif
 
             if (onshell .and. (lbw(nbw).eq. 2) .or.
@@ -313,6 +229,7 @@
 
       double precision prmass(-nexternal:0,lmaxconfigs)
       double precision prwidth(-nexternal:0,lmaxconfigs)
+      double precision prwidth_tmp(-nexternal:0,lmaxconfigs)
       integer pow(-nexternal:0,lmaxconfigs)
 
       integer idup(nexternal,maxproc,maxsproc)
@@ -382,11 +299,21 @@
 
 c-----
 c  Begin Code
-c-----      
+c-----     
+      iconfig = this_config
+c     needs to be initialise to avoid segfault
+      do i = -nexternal,-1
+         prwidth(i,iconfig) = 0
+         prmass(i,iconfig) =0
+      enddo 
       include 'props.inc'
 c      etmin = 10
       nt = 0
-      iconfig = this_config
+      do i = -nexternal,-1
+         prwidth_tmp(i,iconfig) = max(prwidth(i,iconfig), prmass(i,iconfig)*small_width_treatment)
+      enddo
+
+
       mtot = 0d0
       etot = 0d0   !Total energy needed
       spmass = 0d0 !Keep track of BW masses for shat
@@ -423,7 +350,7 @@
 c     Look for identical particles to map radiation processes
       call idenparts(iden_part, iforest(1,-max_branch,iconfig),
      $     sprop(1,-max_branch,iconfig), gForceBW(-max_branch,iconfig),
-     $     prwidth(-nexternal,iconfig))
+     $     prwidth_tmp(-nexternal,iconfig))
 
 c     Start loop over propagators
       do i=-1,-(nexternal-3),-1
@@ -446,15 +373,15 @@
               xm(i)=max(xm(i),max(xqcutij(l1,l2),0d0))
             endif
 c            write(*,*) 'iconfig,i',iconfig,i
-c            write(*,*) prwidth(i,iconfig),prmass(i,iconfig)
-            if (prwidth(i,iconfig) .gt. 0 ) then
+c            write(*,*) prwidth_tmp(i,iconfig),prmass(i,iconfig)
+            if (prwidth_tmp(i,iconfig) .gt. 0 ) then
                nbw=nbw+1
 c              JA 6/8/2011 Set xe(i) for resonances
                if (gforcebw(i,iconfig).eq.1) then
-                  xm(i) = max(xm(i), prmass(i,iconfig)-bwcutoff*prwidth(i,iconfig))
+                  xm(i) = max(xm(i), prmass(i,iconfig)-bwcutoff*prwidth_tmp(i,iconfig))
                   bwcut_for_PS(i) = bwcutoff
                else if (lbw(nbw).eq.1) then
-                  xm(i) = max(xm(i), prmass(i,iconfig)-5d0*prwidth(i,iconfig))
+                  xm(i) = max(xm(i), prmass(i,iconfig)-5d0*prwidth_tmp(i,iconfig))
                   bwcut_for_PS(i) = 5d0
                else
                   bwcut_for_PS(i) = 5d0
@@ -465,19 +392,19 @@
 c     Either: required onshell and daughter masses too large
 c     Or: forced and daughter masses too large
 c     Or: required offshell and forced
-            if(prwidth(i,iconfig) .gt. 0.and.
+            if(prwidth_tmp(i,iconfig) .gt. 0.and.
      $         (lbw(nbw).eq.1.and.
-     $          (prmass(i,iconfig)+bwcut_for_PS(i)*prwidth(i,iconfig).lt.xm(i)
-     $           .or.prmass(i,iconfig)-bwcut_for_PS(i)*prwidth(i,iconfig).gt.dsqrt(stot))
+     $          (prmass(i,iconfig)+bwcut_for_PS(i)*prwidth_tmp(i,iconfig).lt.xm(i)
+     $           .or.prmass(i,iconfig)-bwcut_for_PS(i)*prwidth_tmp(i,iconfig).gt.dsqrt(stot))
      $          .or.gforcebw(i,iconfig).eq.1.and.
-     $              prmass(i,iconfig)+bwcutoff*prwidth(i,iconfig).lt.xm(i)
+     $              prmass(i,iconfig)+bwcutoff*prwidth_tmp(i,iconfig).lt.xm(i)
      $          .or.lbw(nbw).eq.2.and.gforcebw(i,iconfig).eq.1))
      $        then
 c     Write results.dat and quit
                call write_null_results()
                stop
             endif
-            if (prwidth(i,iconfig) .gt. 0 .and. lbw(nbw) .le. 1) then         !B.W.
+            if (prwidth_tmp(i,iconfig) .gt. 0 .and. lbw(nbw) .le. 1) then         !B.W.
                if (i .eq. -(nexternal-(nincoming+1))) then  !This is s-hat
                   j = 3*(nexternal-2)-4+1    !set i to ndim+1
 c-----
@@ -489,14 +416,14 @@
      $                 .or. lbw(nbw).eq.1) then
                      write(*,*) 'Setting PDF BW',j,nbw,prmass(i,iconfig)
                      spole(j)=prmass(i,iconfig)*prmass(i,iconfig)/stot
-                     swidth(j) = prwidth(i,iconfig)*prmass(i,iconfig)/stot
+                     swidth(j) = prwidth(i,iconfig)*prmass(i,iconfig)/stot ! keep the real width here (important for the jacobian)
                   endif
-               else if((prmass(i,iconfig)+bwcut_for_PS(i)*prwidth(i,iconfig)).ge.xm(i)
+               else if((prmass(i,iconfig)+bwcut_for_PS(i)*prwidth_tmp(i,iconfig)).ge.xm(i)
      $                  .and. iden_part(i).eq.0 .or. lbw(nbw).eq.1) then
 c              JA 02/13 Only allow BW if xm below M+5*Gamma
                   write(*,*) 'Setting BW',i,nbw,prmass(i,iconfig)
                   spole(-i)=prmass(i,iconfig)*prmass(i,iconfig)/stot
-                  swidth(-i) = prwidth(i,iconfig)*prmass(i,iconfig)/stot
+                  swidth(-i) = prwidth(i,iconfig)*prmass(i,iconfig)/stot ! keep the real width here (important for the jacobian)
                endif
 c     JA 4/1/2011 Set grid in case there is no BW (radiation process)
                if (swidth(-i) .eq. 0d0 .and.
@@ -509,7 +436,7 @@
 c     Set spmass for BWs
                if (swidth(-i) .ne. 0d0)
      $              spmass=spmass-xm(i) +
-     $              max(xm(i),prmass(i,iconfig)-bwcut_for_PS(i)*prwidth(i,iconfig))
+     $              max(xm(i),prmass(i,iconfig)-bwcut_for_PS(i)*prwidth_tmp(i,iconfig))
             else                                  !1/x^pow
               a=prmass(i,iconfig)**2/stot
 c     JA 4/1/2011 always set grid
@@ -530,7 +457,7 @@
                  endif
               endif
               if (xo.eq.0d0) xo=1d0/stot
-c              if (prwidth(i, iconfig) .eq. 0d0.or.iden_part(i).gt.0) then 
+c              if (prwidth_tmp(i, iconfig) .eq. 0d0.or.iden_part(i).gt.0) then 
               call setgrid(-i,xo,a,1)
 c              else 
 c                 write(*,*) 'Using flat grid for BW',i,nbw,
@@ -610,7 +537,7 @@
             spole(i)= -2.0d0    ! 1/s pole
             write(*,*) "Transforming s_hat 1/s ",i,xo, smin, stot
          else
-            write(*,*) "Transforming s_hat BW ",spole(i),swidth(i)
+            write(*,*) "Transforming s_hat BW ",spole(i), max(swidth(i), spole(i)*small_width_treatment)
          endif
       endif
 
 
=== modified file 'Template/LO/SubProcesses/reweight.f'
--- Template/LO/SubProcesses/reweight.f	2018-06-08 14:42:29 +0000
+++ Template/LO/SubProcesses/reweight.f	2018-11-09 10:32:12 +0000
@@ -6,6 +6,7 @@
 c     ----
 c     p p > t t~ (up to 2jet)
 c     p p > w+ (up to 3 jet)
+c     p p > j j w+ (ordering seems important) 
 c     p p > z t t~ j j  (no MLM needed)
 c
 c
@@ -243,13 +244,15 @@
       integer iddgluon, iddother, idgluon, idother
       logical isqcd
       external isqcd
+      integer get_color
+      external get_color 
 
       idmo=ipdg(imo)
       idda1=ipdg(ida1)
       idda2=ipdg(ida2)
 
       if (btest(mlevel,4)) then
-        write(*,*) ' updating ipart for: ',ida1,ida2,' -> ',imo
+        write(*,*) 'updating ipart for: ',ida1,ida2,' -> ',imo
       endif
 
       if (btest(mlevel,4)) then
@@ -291,7 +294,6 @@
      $        ' (',ipdg(imo),')'
          return
       endif        
-
 c     FS clustering
 c     Transmit parton PDG code for parton vertex
       if(isjet(idmo)) then
@@ -370,14 +372,26 @@
 c     quark -> gluon-quark or Z-quark or h-quark or W-quark
         ipart(1,imo)=ipart(1,ida2)
         ipart(2,imo)=0
-      else
+      else if(iabs(get_color(idmo)).eq.3.and.iabs(get_color(idda1)).eq.3.and.get_color(idda2).eq.1) then
+c       exotic q > q' Scalar         
+        ipart(1,imo)=ipart(1,ida1)
+        ipart(2,imo)=0
+      else if(iabs(get_color(idmo)).eq.3.and.iabs(get_color(idda2)).eq.3.and.get_color(idda1).eq.1) then
+c       exotic q > Scalar q'
+        ipart(1,imo)=ipart(1,ida2)
+        ipart(2,imo)=0
+      else if (get_color(idmo).eq.1) then
 c     Color singlet
          ipart(1,imo)=ipart(1,ida1)
          ipart(2,imo)=ipart(1,ida2)
+      else
+         write(*,*) idmo,'>', idda1, idda2, 'color', get_color(idmo),'>', get_color(idda1), get_color(idda2)
+         write(*,*) "failed for ipartupdate. Please retry without MLM/default dynamical scale"
+         stop 3
       endif
       
       if (btest(mlevel,4)) then
-        write(*,*) ' -> ',(ipart(i,imo),i=1,2),' (',ipdg(imo),')'
+        write(*,*) 'XY -> ',(ipart(i,imo),i=1,2),' (',ipdg(imo),')'
       endif
 
       return
@@ -695,6 +709,7 @@
 c     increasecode gives whether we should increase jcode at next vertex
       increasecode=.false.
       do n=1,nexternal-2
+c         write(*,*) 'QCD jet status (before n= ',n,'):',(iqjets(i),i=3,nexternal)
         do i=1,2 ! index of the child in the interaction
           do j=1,2 ! j index of the beam
             if(idacl(n,i).eq.ibeam(j))then
@@ -769,7 +784,6 @@
                  pdgm   = ipdgcl(imocl(n),igraphs(1),iproc)
                  pdgid1 = ipdgcl(idacl(n,1),igraphs(1),iproc)
                  pdgid2 = ipdgcl(idacl(n,2),igraphs(1),iproc)
-
                  if (.not.isqcd(pdgm).and..not.isqcd(pdgid1).and..not.isqcd(pdgid2)) then
                     ! this is to avoid to do weird stuff for w+ w- z (or h h h) 
                     ! this fix an issue for qq_zttxqq G1594.08
@@ -914,8 +928,8 @@
             endif
          enddo
          njetstore(iconfig)=njets
-         if (btest(mlevel,4))
-     $        write(*,*) 'Storing jets: ',(iqjetstore(i,iconfig),i=1,njets)
+        if (btest(mlevel,4))
+     $       write(*,*) 'Storing jets: ',(iqjetstore(i,iconfig),i=1,njets)
 c     Recluster without requiring chcluster
          goto 100
       else
 
=== modified file 'Template/NLO/MCatNLO/Makefile_MadFKS'
--- Template/NLO/MCatNLO/Makefile_MadFKS	2016-02-24 00:19:50 +0000
+++ Template/NLO/MCatNLO/Makefile_MadFKS	2018-11-09 10:32:12 +0000
@@ -1,7 +1,7 @@
 -include ../../Source/make_opts
 
 DEBUG=
-FF=$(FC) $(FFLAGS) $(DEBUG)
+FF=$(FC) $(FFLAGS) $(DEBUG) -std=legacy
 CPP=$(CXX) $(CXXFLAGS) $(DEBUG) $(INCLOPTION)
 CC=$(CXX) $(CFLAGS) $(DEBUG) $(INCLOPTION)
 
 
=== modified file 'Template/NLO/SubProcesses/check_poles.f'
--- Template/NLO/SubProcesses/check_poles.f	2017-06-18 20:58:19 +0000
+++ Template/NLO/SubProcesses/check_poles.f	2018-11-09 10:32:12 +0000
@@ -19,7 +19,7 @@
       include 'genps.inc'
       include 'nexternal.inc'
       include 'nFKSconfigs.inc'
-      double precision p(0:3, nexternal), prambo(0:3, nexternal)
+      double precision p(0:3, nexternal), prambo(0:3,100)
       double precision p_born(0:3,nexternal-1)
       common/pborn/p_born
       double precision pswgt
@@ -49,7 +49,7 @@
       include 'coupl.inc'
       include 'q_es.inc'
       integer nsqso      
-      double precision pmass(nexternal), pmass_rambo(nexternal)
+      double precision pmass(nexternal), pmass_rambo(100)
       integer nfail
       logical first_time
       data first_time/.TRUE./
 
=== modified file 'Template/NLO/SubProcesses/combine_plots_FO.sh'
--- Template/NLO/SubProcesses/combine_plots_FO.sh	2013-05-25 12:52:41 +0000
+++ Template/NLO/SubProcesses/combine_plots_FO.sh	2018-11-09 10:32:12 +0000
@@ -1,6 +1,4 @@
 #!/bin/bash
-TD=td_mac_intel
-PSOPEN=open
 
 # $? is the value of last executed command. A call to this function
 # after a failure will cause the program to quit the script
@@ -15,126 +13,74 @@
 fi
 }
 
-thisdir=`pwd`
-if [ -f $thisdir/MADatNLO.top ]
-then
-rm -f $thisdir/MADatNLO*.top
-fi
-make read40
-EXENAME=$thisdir/read40
+function combine {
+counter=0
+i=1
 echo -n "" > dir
-counterp=0
-for p in P* ; do
-    cd $thisdir
-    cd $p
-    rm -f MADatNLO*
-    rm -f *read40*.out
-    echo -n "" > dir
-    counter=0
-    i=1
-    for g in $* ; do
-        if [ ! -f $g/MADatNLO*.top ]
-        then
-          echo "No topdrawer file in $p/$g"
-	else
-	counter=$(($counter + 1))
-	echo $g"/MADatNLO.top" >> dir
-	if [[ $(($counter % 40)) == 0 ]]; then
-            echo -n "" > dir2
-	    echo "40" > dir2
-	    cat dir >> dir2
-	    $EXENAME < dir2
-            teststatus Failure in step 1
-	    rm -f dir
-	    rm -f dir2
-	    mv fort.88 MADatNLO$i\.top 
-	    mv read40.out "S1read40_"$i\.out
-	    i=$(($i + 1))
-	    echo -n "" > dir
-        fi
-	fi
-    done
-    if [[ $(($counter % 40)) != 0 ]] ; then
-	echo $(($counter % 40)) > dir2
+for plot_file in "${@:2}" ; do
+    counter=$(($counter + 1))
+    echo $plot_file >> dir
+    if [[ $(($counter % 40)) == 0 ]]; then
+        echo -n "" > dir2
+	echo "40" > dir2
 	cat dir >> dir2
-	while [[ $(($counter % 40)) -ne 0 ]]; do
-	    counter=$[$counter-1]
-	done
-	$EXENAME < dir2
-	teststatus Failure in step 2
+	./read40 < dir2
+        teststatus Failure in step 1
 	rm -f dir
 	rm -f dir2
-	mv fort.88 MADatNLO$i\.top 
-	mv read40.out "S2read40_"$i\.out
-    fi
-	
-    counter=0
-    for m in MADatNLO*.top ; do
-	counter=$[$counter+1]
-	echo $m >> dir
-    done
-    if [[ $counter -gt 1 ]] ; then
-	echo $counter > dir2
-	cat dir >> dir2
-	$EXENAME < dir2
-	teststatus Failure in step 3
-	rm -f dir2
-	mv fort.88 MADatNLO.top 
-	mv read40.out S3read40.out
-    else
-	mv MADatNLO1.top MADatNLO.top
-    fi
-    rm -f dir
-    cd $thisdir
-    counterp=$(($counterp + 1))
-    echo $p"/MADatNLO.top" >> dir
-    echo $counterp $p "done"
-done
-
-cd $thisdir
-make split40
-./split40
-
-i=0
-#pdir* are created by split
-for p in pdir* ; do
-    i=$(($i+1))
-    $EXENAME < $p    
-    teststatus Failure in step 4
-    mv fort.88 MADatNLO$i\.top
-    mv read40.out "S4read40_"$i\.out
-    rm $p
-    echo $i $p "done"
-done
-
-if [[ $i -gt 1 ]] ; then
-    echo $i > dir2
-    for p in MADatNLO*.top ; do
-	echo $p >>dir2
-    done
-    $EXENAME < dir2
-    teststatus Failure in step 5
+	mv fort.88 $1_$i\.top 
+	mv read40.out "S1read40_"$1_$i\.out
+	i=$(($i + 1))
+	echo -n "" > dir
+    fi
+done
+if [[ $(($counter % 40)) != 0 ]] ; then
+    echo $(($counter % 40)) > dir2
+    cat dir >> dir2
+    ./read40 < dir2
+    teststatus Failure in step 2
     rm -f dir
     rm -f dir2
-    mv fort.88 MADatNLO.top 
-    mv read40.out S5read40.out
+    mv fort.88 $1_$i\.top 
+    mv read40.out "S2read40_"$1_$i\.out
+fi
+}
+
+make read40
+rm -f MADatNLO*.top
+
+if [ "$#" -gt 64000 ]; then
+    combine "MADatNLO_A" "$@"
+    combine "MADatNLO_B" MADatNLO_A*.top
+    combine "MADatNLO_C" MADatNLO_B*.top
+    combine "MADatNLO_D" MADatNLO_C*.top
+    mv MADatNLO_D_1.top MADatNLO.top
+elif [ "$#" -gt 1600 ]; then
+    combine "MADatNLO_A" "$@"
+    combine "MADatNLO_B" MADatNLO_A*.top
+    combine "MADatNLO_C" MADatNLO_B*.top
+    mv MADatNLO_C_1.top MADatNLO.top
+elif [ "$#" -gt 40 ]; then
+    combine "MADatNLO_A" "$@"
+    combine "MADatNLO_B" MADatNLO_A*.top
+    mv MADatNLO_B_1.top MADatNLO.top
+elif [ "$#" -gt 1 ]; then
+    combine "MADatNLO_A" "$@"
+    mv MADatNLO_A_1.top MADatNLO.top
 else
-    mv MADatNLO1.top MADatNLO.top
+    cp $1 MADatNLO.top
 fi
 
-if [ -f $thisdir/read40.errors ]
+
+if [ -f read40.errors ]
 then
-rm -f $thisdir/read40.errors
+rm -f read40.errors
 fi
-find . -name S"*"read40"*".out -exec fgrep -il "READ40 ERROR READ40 ERROR" {} \; > $thisdir/read40.errors
-
-if [ ! -s $thisdir/read40.errors ] ; then
+
+find . -name S"*"read40"*".out -exec fgrep -il "READ40 ERROR READ40 ERROR" {} \; > read40.errors
+
+if [ ! -s read40.errors ] ; then
     echo 'no errors found'
 else
     echo 'ERRORS! see read40.errors for details'
 fi
-
-#./NLO_Born3.py
-#$TD MADatNLO_combined.top "POSTSCR,ORIENT=3"
-#$PSOPEN MADatNLO_combined.ps
-##td2pdf MADatNLO_con
 
=== modified file 'Template/NLO/SubProcesses/driver_mintFO.f'
--- Template/NLO/SubProcesses/driver_mintFO.f	2017-09-21 13:47:31 +0000
+++ Template/NLO/SubProcesses/driver_mintFO.f	2018-11-09 10:32:12 +0000
@@ -48,7 +48,7 @@
       include "mint.inc"
       integer nhits_in_grids(maxchannels)
       real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
-     $     ,ndimmax,maxchannels),ymax_virt(maxchannels),ans(nintegrals
+     $     ,ndimmax,maxchannels),ymax_virt(0:maxchannels),ans(nintegrals
      $     ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals
      $     ,0:maxchannels),x(ndimmax),itmax_fl
       integer ixi_i,iphi_i,iy_ij,vn
 
=== modified file 'Template/NLO/SubProcesses/driver_mintMC.f'
--- Template/NLO/SubProcesses/driver_mintMC.f	2017-09-21 13:47:31 +0000
+++ Template/NLO/SubProcesses/driver_mintMC.f	2018-11-09 10:32:12 +0000
@@ -55,7 +55,7 @@
       common /event_normalisation/event_norm
 c For MINT:
       real* 8 xgrid(0:nintervals,ndimmax,maxchannels),ymax(nintervals
-     $     ,ndimmax,maxchannels),ymax_virt(maxchannels),ans(nintegrals
+     $     ,ndimmax,maxchannels),ymax_virt(0:maxchannels),ans(nintegrals
      $     ,0:maxchannels),unc(nintegrals,0:maxchannels),chi2(nintegrals
      $     ,0:maxchannels),x(ndimmax)
       integer ixi_i,iphi_i,iy_ij,vn,nhits_in_grids(maxchannels)
 
=== modified file 'Template/NLO/SubProcesses/makefile'
--- Template/NLO/SubProcesses/makefile	2017-06-15 16:22:23 +0000
+++ Template/NLO/SubProcesses/makefile	2018-11-09 10:32:12 +0000
@@ -7,9 +7,6 @@
 # Files for the read40 combiner of top drawer files
 READ40=read40.o
 
-# Files for the split40 splitter of top drawer files
-SPLIT40=split40.o
-
 # Files to collect all the events in the separate integration channels into a single file
 COLLECT_EVENTS=collect_events.o handling_lhe_events.o fill_MC_mshell.o
 
@@ -21,9 +18,6 @@
 read40: $(READ40)
 	$(FC) $(LDFLAGS) -o read40 $(READ40)
 
-split40: $(SPLIT40)
-	$(FC) $(LDFLAGS) -o split40 $(SPLIT40)
-
 collect_events: $(COLLECT_EVENTS) $(LIBS)
 	$(FC) -o collect_events $(COLLECT_EVENTS) $(LINKLIBS) $(LDFLAGS)
 	rm handling_lhe_events.o
 
=== modified file 'Template/NLO/SubProcesses/mint-integrator2.f'
--- Template/NLO/SubProcesses/mint-integrator2.f	2017-05-09 14:07:18 +0000
+++ Template/NLO/SubProcesses/mint-integrator2.f	2018-11-09 10:32:12 +0000
@@ -77,7 +77,7 @@
      $     ,ans(nintegrals,0:maxchannels),unc(nintegrals,0:maxchannels)
      $     ,ans3(nintegrals,3),unc3(nintegrals,3),ans_l3(nintegrals)
      $     ,unc_l3(nintegrals),chi2_l3(nintegrals),vol_chan
-      real * 8 xint_virt(maxchannels),ymax_virt(maxchannels)
+      real * 8 xint_virt(maxchannels),ymax_virt(0:maxchannels)
      $     ,ans_chan(0:maxchannels)
       real * 8 x(ndimmax),vol
       real * 8 xacc(0:nintervals,ndimmax,maxchannels)
@@ -985,8 +985,8 @@
       integer ndim,imode
       include "mint.inc"
       real * 8 fun,xgrid(0:nintervals,ndimmax,maxchannels)
-     $     ,ymax(nintervals,ndimmax,maxchannels),ymax_virt(maxchannels)
-     $     ,x(ndimmax)
+     $     ,ymax(nintervals,ndimmax,maxchannels)
+     $     ,ymax_virt(0:maxchannels),x(ndimmax)
       real * 8 dx(ndimmax),xx(ndimmax),vol_chan,dummy
       integer icell(ndimmax),ncell(ndimmax),ncell_virt
       integer ifold(ndimmax),kfold(ndimmax)
 
=== modified file 'Template/NLO/SubProcesses/read40.for'
--- Template/NLO/SubProcesses/read40.for	2012-08-28 21:06:34 +0000
+++ Template/NLO/SubProcesses/read40.for	2018-11-09 10:32:12 +0000
@@ -5,7 +5,7 @@
       real * 8 xfacts(1:40),step(1:40)
       real * 8 xv(1:40,1:1000),yv(1:40,1:1000),er(1:40,1:1000)
       real * 8 xa(1:1000),ya(1:1000),ea(1:1000)
-      character * 70 file,file1,filename(1:40)
+      character * 700 file,file1,filename(1:40)
       character * 70 line,vartype,vartype2
       character * 70 strscale,strlimx,strlimy,strtit
       character * 11 chref1,chref7
 
=== removed file 'Template/NLO/SubProcesses/split40.f'
--- Template/NLO/SubProcesses/split40.f	2012-08-28 21:06:34 +0000
+++ Template/NLO/SubProcesses/split40.f	1970-01-01 00:00:00 +0000
@@ -1,51 +0,0 @@
-      PROGRAM SPLIT
-      IMPLICIT NONE
-      character*80 buff
-      integer number,i,j
-
-      open (unit=1, file="dir", status="old")
-
-      number=0
-      do 
-         read(1,*,end=99,err=99) buff
-         number=number+1
-      enddo
- 99   continue
-      rewind(1)
-      i=1
-      do while (number.gt.40)
-         open(unit=2,file='pdir'//char(48+i),status='unknown')
-         write (2,'(i3)') 40
-         do j=1,40
-            read(1,'(a)',end=98,err=98) buff
-            write (2,'(a)') buff
-         enddo
-c         do j=1,9
-c            write (2,'(i1)') 1
-c         enddo
-         close(2)
-         i=i+1
-         number=number-40
-      enddo
-      if (number.gt.0) then
-         open(unit=2,file='pdir'//char(48+i),status='unknown')
-         if (number.gt.99) then
-            write (2,'(i3)') number
-         elseif(number.gt.9)then
-            write (2,'(i2)') number
-         else
-            write (2,'(i1)') number
-         endif
-         do j=1,number
-            read(1,'(a)',end=98,err=98) buff
-            write (2,'(a)') buff
-         enddo
-c         do j=1,number
-c            write (2,'(i1)') 1
-c         enddo
-         close(2)
-      endif
-
-
- 98   return
-      end
 
=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2018-06-22 06:52:02 +0000
+++ UpdateNotes.txt	2018-11-09 10:32:12 +0000
@@ -1,5 +1,22 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
+2.6.4 (//)
+      OM: add specific treatement for small width (at LO only and not for loop-induced)
+          if the width is smaller than 1e-6 times the mass, a fake width (at that value) is used for the
+          numerical evaluation of the matrix-element. S-channel resonances are re-scaled according to 
+          narrow-width approximation to return the correct total cross-section (the distribution of events 
+          will on the other hand follow the new width). 
+          The parameter '1e-6' can be changed by adding to (LO) run_card the parameter: "small_width_treatment"
+      OM: add a new command "install looptools" to trigger the question that is automatically trigger 
+          the first time a loop computation is needed.
+      RF: Fixed a bug when using TopDrawer plots for f(N)LO runs, where the combination of the plots could lead
+          to completely wrong histograms/distributions in case of high-precision runs.
+      OM: Fix some MLM crash for some processes (in particular BSM processes with W'). 
+      OM: Fix a bug in the reweighting due to the new lhe format (the one avoiding some issue with py8)
+      OM: Fix a behavior for negative mass, the width was set to negative in the param_card automatically
+          making the Parton-shower (and other code) to crash since this does not follow the convention.
+      OM: Change compiler flag to support Mojave.
+
 2.6.3.2 (22/06/18)
       OM: Fix a bug in auto-width when mass are below QCD scale.
       OM: Fix a bug for g b initial state where the mass in the lhe file was not always correctly assigned
 
=== modified file 'VERSION'
--- VERSION	2018-06-22 06:52:02 +0000
+++ VERSION	2018-11-09 10:32:12 +0000
@@ -1,5 +1,5 @@
-version = 2.6.3.2
-date = 2018-06-22
+version = 2.6.4
+date = 2018-11-30
 
 
 
 
=== modified file 'aloha/aloha_writers.py'
--- aloha/aloha_writers.py	2018-05-23 14:15:47 +0000
+++ aloha/aloha_writers.py	2018-11-09 10:32:12 +0000
@@ -523,7 +523,7 @@
         if 'MP' in self.tag:
             out.write(' complex*32 CI\n')
             if KERNEL.has_pi:
-                out.write(' double*16 PI\n')
+                out.write(' REAL ( KIND = 16 ) PI\n')
         else:
             out.write(' complex*16 CI\n')
             if KERNEL.has_pi:
@@ -673,6 +673,8 @@
             self.has_model_parameter = True
             if name.lower() in ['pi', 'as', 'mu_r', 'aewm1','g']:
                 return name
+            if name.startswith(aloha.aloha_prefix):
+                return name
             return '%s%s' % (aloha.aloha_prefix, name)
         
         if '_' in name:
 
=== modified file 'madgraph/core/base_objects.py'
--- madgraph/core/base_objects.py	2018-04-23 20:23:22 +0000
+++ madgraph/core/base_objects.py	2018-11-09 10:32:12 +0000
@@ -1487,7 +1487,7 @@
     def change_parameter_name_with_prefix(self, prefix='mdl_'):
         """ Change all model parameter by a given prefix.
         Modify the parameter if some of them are identical up to the case"""
-
+        
         lower_dict={}
         duplicate = set()
         keys = self.get('parameters').keys()
@@ -1568,7 +1568,13 @@
             for key in self['couplings'].keys():
                 for coup in self['couplings'][key]:
                     coup.expr = rep_pattern.sub(replace, coup.expr)
-                    
+
+            # change form-factor
+            ff = [l.formfactors for l in self['lorentz'] if hasattr(l, 'formfactors')]
+            ff = set(sum(ff,[])) # here we have the list of ff used in the model
+            for f in ff:
+                f.value = rep_pattern.sub(replace, f.value)
+
             # change mass/width
             for part in self['particles']:
                 if str(part.get('mass')) in one_change:
 
=== modified file 'madgraph/core/helas_objects.py'
--- madgraph/core/helas_objects.py	2017-08-10 13:09:28 +0000
+++ madgraph/core/helas_objects.py	2018-11-09 10:32:12 +0000
@@ -4569,6 +4569,13 @@
         return sum([d.get('wavefunctions') for d in \
                        self.get('diagrams')], [])
 
+
+    def get_all_mass_widths(self):
+        """Gives a list of all widths used by this ME (from propagator)"""
+
+        return set([(d.get('mass'),d.get('width')) for d in self.get_all_wavefunctions()])
+
+
     def get_all_amplitudes(self):
         """Gives a list of all amplitudes for this ME"""
 
 
=== modified file 'madgraph/interface/amcatnlo_interface.py'
--- madgraph/interface/amcatnlo_interface.py	2018-06-08 13:35:22 +0000
+++ madgraph/interface/amcatnlo_interface.py	2018-11-09 10:32:12 +0000
@@ -445,7 +445,11 @@
         proc_type=self.extract_process_type(line)
         if proc_type[1] not in ['real', 'LOonly']:
             run_interface.check_compiler(self.options, block=False)
+        #validate_model will reset self._generate_info; to avoid
+        #this store it
+        geninfo = self._generate_info
         self.validate_model(proc_type[1], coupling_type=proc_type[2])
+        self._generate_info = geninfo
 
         #now generate the amplitudes as usual
         #self.options['group_subprocesses'] = 'False'
 
=== modified file 'madgraph/interface/amcatnlo_run_interface.py'
--- madgraph/interface/amcatnlo_run_interface.py	2018-06-14 15:24:33 +0000
+++ madgraph/interface/amcatnlo_run_interface.py	2018-11-09 10:32:12 +0000
@@ -1053,7 +1053,10 @@
         
         if self.last_mode in ['LO', 'NLO']:
             self.switch['fixed_order'] = 'ON'
-        self.switch['fixed_order'] = 'OFF'
+        if self.proc_characteristics['ninitial'] == 1:
+            self.switch['fixed_order'] = 'ON'
+        else:
+            self.switch['fixed_order'] = 'OFF'
 
     def color_for_fixed_order(self, switch_value):
          
@@ -1158,6 +1161,11 @@
         if self.last_mode in ['LO', 'NLO', 'noshower', 'noshowerLO']:
             self.switch['shower'] = 'OFF'
             return 
+        
+        if self.proc_characteristics['ninitial'] == 1:
+            self.switch['shower'] = 'OFF'
+            return
+        
          
         if os.path.exists(pjoin(self.me_dir, 'Cards', 'shower_card.dat')):
             self.switch['shower'] = self.run_card['parton_shower']  
@@ -2603,7 +2611,13 @@
         devnull = open(os.devnull, 'w') 
         
         if self.analyse_card['fo_analysis_format'].lower() == 'topdrawer':
-            misc.call(['./combine_plots_FO.sh'] + folder_name, \
+            topfiles = []
+            for job in jobs:
+                if job['dirname'].endswith('.top'):
+                    topfiles.append(job['dirname'])
+                else:
+                    topfiles.append(pjoin(job['dirname'],'MADatNLO.top'))
+            misc.call(['./combine_plots_FO.sh'] + topfiles, \
                       stdout=devnull, 
                       cwd=pjoin(self.me_dir, 'SubProcesses'))
             files.cp(pjoin(self.me_dir, 'SubProcesses', 'MADatNLO.top'),
 
=== modified file 'madgraph/interface/common_run_interface.py'
--- madgraph/interface/common_run_interface.py	2018-06-22 06:52:02 +0000
+++ madgraph/interface/common_run_interface.py	2018-11-09 10:32:12 +0000
@@ -963,13 +963,15 @@
         else:
             return None
 
-    def ask_edit_cards(self, cards, mode='fixed', plot=True, first_cmd=None):
+    def ask_edit_cards(self, cards, mode='fixed', plot=True, first_cmd=None, from_banner=None,
+                       banner=None):
         """ """
         if not self.options['madanalysis_path']:
             plot = False
 
         self.ask_edit_card_static(cards, mode, plot, self.options['timeout'],
-                                  self.ask, first_cmd=first_cmd)
+                                  self.ask, first_cmd=first_cmd, from_banner=from_banner,
+                                  banner=banner)
         
         for c in cards:
             if not os.path.isabs(c):
@@ -3240,7 +3242,19 @@
         - Check that no width are too small (raise a warning if this is the case)
         3) if dependent is on True check for dependent parameter (automatic for scan)"""
         
-        return self.static_check_param_card(path, self, run=run, dependent=dependent)
+        self.static_check_param_card(path, self, run=run, dependent=dependent)
+        
+        card = param_card_mod.ParamCard(path)
+        for param in card['decay']:
+            width = param.value
+            if width == 0:
+                continue
+            try:
+                mass = card['mass'].get(param.lhacode).value
+            except Exception:
+                continue
+        
+        
         
     @staticmethod
     def static_check_param_card(path, interface, run=True, dependent=False, 
@@ -4138,6 +4152,21 @@
             return CommonRunCmd.install_lhapdf_pdfset_static(lhapdf_config, pdfsets_dir, 
                                                               filename.replace('.LHgrid',''), 
                                         lhapdf_version, alternate_path)
+        elif lhapdf_version.startswith('6.'):
+            # try to do a simple wget
+            wwwpath = "http://www.hepforge.org/archive/lhapdf/pdfsets/%s/%s.tar.gz" 
+            wwwpath %= ('.'.join(lhapdf_version.split('.')[:2]), filename)
+            misc.wget(wwwpath, pjoin(pdfsets_dir, '%s.tar.gz' %filename))
+            misc.call(['tar', '-xzpvf', '%s.tar.gz' %filename],
+                      cwd=pdfsets_dir)
+            if os.path.exists(pjoin(pdfsets_dir, filename)) or \
+               os.path.isdir(pjoin(pdfsets_dir, filename)):
+                logger.info('%s successfully downloaded and stored in %s' \
+                        % (filename, pdfsets_dir))  
+            else:
+                raise MadGraph5Error, \
+                'Could not download %s into %s. Please try to install it manually.' \
+                    % (filename, pdfsets_dir)                          
             
         else:
             raise MadGraph5Error, \
@@ -4342,7 +4371,7 @@
 
      
     
-    def __init__(self, question, cards=[], mode='auto', *args, **opt):
+    def __init__(self, question, cards=[], from_banner=None, banner=None, mode='auto', *args, **opt):
 
 
         self.load_default()        
@@ -4365,7 +4394,8 @@
         self.all_vars = set()
         self.modified_card = set() #set of cards not in sync with filesystem
                               # need to sync them before editing/leaving
-
+        self.init_from_banner(from_banner, banner)
+        
         #update default path by custom one if specify in cards
         for card in cards:
             if os.path.exists(card) and os.path.sep in cards:
@@ -4379,11 +4409,31 @@
             new_conflict = self.all_vars.intersection(new_vars)
             self.conflict.union(new_conflict)
             self.all_vars.union(new_vars)
+            
+
+    def init_from_banner(self, from_banner, banner):
+        """ defined card that need to be initialized from the banner file 
+            from_banner should be a list of card to load from the banner object
+        """
+
+        if from_banner is None:
+            self.from_banner = {}
+            return
+        misc.sprint(from_banner)
+        self.from_banner = {}
+        for card in from_banner:
+            self.from_banner[card] = banner.charge_card(card)
+        return self.from_banner
+    
 
     def get_path(self, name, cards):
         """initialise the path if requested"""
 
         defname = '%s_default' % name
+        
+        if name in self.from_banner:
+            return self.from_banner[name]
+        
         if isinstance(cards, list):
             if name in cards:
                 return True
@@ -4417,7 +4467,13 @@
         self.pname2block = {}
         self.restricted_value = {}
         self.param_card = {}
-        if not self.get_path('param', cards):
+        
+        is_valid_path = self.get_path('param', cards)
+        if not is_valid_path:
+            self.param_consistency = False
+            return []
+        if isinstance(is_valid_path, param_card_mod.ParamCard):
+            self.param_card = is_valid_path
             self.param_consistency = False
             return []
 
@@ -4441,10 +4497,15 @@
         return self.pname2block.keys()
         
     def init_run(self, cards):
-        
+      
         self.run_set = []
-        if not self.get_path('run', cards):
-            return []
+        is_valid_path = self.get_path('run', cards)
+        if not is_valid_path:
+            return []
+        if isinstance(is_valid_path, banner_mod.RunCard):
+            self.run_card = is_valid_path
+            return []
+        
         
         try:
             self.run_card = banner_mod.RunCard(self.paths['run'], consistency='warning')
@@ -4455,6 +4516,7 @@
         except IOError:
             run_card_def = {}
 
+
         if run_card_def:
             if self.run_card:
                 self.run_set = run_card_def.keys() + self.run_card.hidden_param
@@ -5129,7 +5191,6 @@
                         logger.warning(str(e))
             return
 
-        
         start = 0
         if len(args) < 2:
             logger.warning('Invalid set command %s (need two arguments)' % line)
@@ -5206,6 +5267,7 @@
             
         if args[0] in ['run_card', 'param_card', 'MadWeight_card', 'shower_card',
                        'delphes_card','madanalysis5_hadron_card','madanalysis5_parton_card']:
+
             if args[1] == 'default':
                 logger.info('replace %s by the default card' % args[0],'$MG:BOLD')
                 files.cp(self.paths['%s_default' %args[0][:-5]], self.paths[args[0][:-5]])
@@ -5262,7 +5324,10 @@
 
         #### RUN CARD
         if args[start] in [l.lower() for l in self.run_card.keys()] and card in ['', 'run_card']:
+
             if args[start] not in self.run_set:
+                if card in self.from_banner or 'run' in self.from_banner:
+                    raise Exception, "change not allowed for this card: event already generated!"
                 args[start] = [l for l in self.run_set if l.lower() == args[start]][0]
 
             if args[start] in self.conflict and card == '':
@@ -5279,9 +5344,10 @@
                     logger.info('remove information %s from the run_card' % args[start],'$MG:BOLD')
                     del self.run_card[args[start]]
             else:
-                if args[0].startswith('sys_') or \
-                   args[0] in self.run_card.list_parameter or \
-                   args[0] in self.run_card.dict_parameter:
+                lower_name = args[0].lower()
+                if lower_name.startswith('sys_') or \
+                   lower_name in self.run_card.list_parameter or \
+                   lower_name in self.run_card.dict_parameter:
                     val = ' '.join(args[start+1:])
                     val = val.split('#')[0]
                 else:
@@ -5626,6 +5692,30 @@
                  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\"."
             
+            # check the status of small width status from LO
+            for param in self.param_card['decay']:
+                width = param.value
+                if width == 0 or isinstance(width,str):
+                    continue
+                try:
+                    mass = self.param_card['mass'].get(param.lhacode).value
+                except Exception:
+                    continue
+                if isinstance(mass,str):
+                    continue
+                
+                if mass:
+                    if abs(width/mass) < self.run_card['small_width_treatment']:
+                        logger.warning("Particle %s will use a fake width  ( %s instead of %s ).\n" +
+                          "Cross-section will be rescaled according to NWA if needed."  +
+                          "To force exact treatment reduce the value of 'small_width_treatment' parameter of the run_card",
+                          param.lhacode[0], mass*self.run_card['small_width_treatment'], width)
+                    elif abs(width/mass) < 1e-12:
+                        logger.error('The width of particle %s is too small for an s-channel resonance (%s). If you have this particle in an s-channel, this is likely to create numerical instabilities .', param.lhacode[0], width)
+                    if CommonRunCmd.sleep_for_error:
+                        time.sleep(5)
+                        CommonRunCmd.sleep_for_error = False
+
 
         ########################################################################
         #       NLO specific check
 
=== modified file 'madgraph/interface/extended_cmd.py'
--- madgraph/interface/extended_cmd.py	2018-05-16 08:43:26 +0000
+++ madgraph/interface/extended_cmd.py	2018-11-09 10:32:12 +0000
@@ -2618,6 +2618,13 @@
 
     def postcmd(self, stop, line):
         
+        # for diamond class arch where both branch defines the postcmd
+        # set it up to be in coop mode
+        try:
+            out = super(ControlSwitch,self).postcmd(stop, line)
+        except AttributeError:
+            pass
+        
         line = line.strip()
         if ';' in line:
             line= [l for l in line.split(';') if l][-1] 
@@ -2676,6 +2683,7 @@
     def check_consistency(self, key, value):
         """check the consistency of the new flag with the old ones"""
         
+        
         if key in self.last_changed:
             self.last_changed.remove(key)
         self.last_changed.append(key)
 
=== modified file 'madgraph/interface/loop_interface.py'
--- madgraph/interface/loop_interface.py	2018-03-22 14:00:54 +0000
+++ madgraph/interface/loop_interface.py	2018-11-09 10:32:12 +0000
@@ -21,6 +21,7 @@
 import time
 import logging
 import re
+import sys
 
 import madgraph
 from madgraph import MG4DIR, MG5DIR, MadGraph5Error
@@ -494,14 +495,19 @@
         aloha.mp_precision = aloha_original_quad_mode
 
 
-    def install_reduction_library(self):
+    def install_reduction_library(self, force=False):
         """Code to install the reduction library if needed"""
         
         opt = self.options
-        
+                
         # Check if first time:
-        if (opt['ninja'] is None) or (os.path.isfile(pjoin(MG5DIR, opt['ninja'],'libninja.a'))): 
+        if not force and ((opt['ninja'] is None) or (os.path.isfile(pjoin(MG5DIR, opt['ninja'],'libninja.a')))): 
             return
+
+        # do not trigger the question for tests
+        if 'test_manager.py' in sys.argv[0]:
+            from unittest.case import SkipTest
+            raise SkipTest
         
         logger.info("First output using loop matrix-elements has been detected. Now asking for loop reduction:", '$MG:BOLD')
         to_install = self.ask('install', '0',  ask_class=AskLoopInstaller, timeout=300, 
@@ -950,9 +956,11 @@
             if os.path.exists(pjoin(install_dir1, 'collier')):
                 self.code['collier'] =  pjoin(install_dir1, 'collier')
             if os.path.exists(pjoin(install_dir2, 'PJFry','bin','qd-config')):
-                self.code['collier'] =  pjoin(install_dir2, 'PJFry')
+                self.code['pjfry'] =  pjoin(install_dir2, 'PJFry')
             if os.path.exists(pjoin(install_dir2, 'golem95')):
-                self.code['collier'] =  pjoin(install_dir2, 'golem95')
+                self.code['glem'] =  pjoin(install_dir2, 'golem95')
+            if os.path.exists(pjoin(install_dir1, 'ninja')):
+                self.code['ninja'] =  pjoin(install_dir2, 'ninja','lib')
         
         # 1. create the question
         question, allowed_answer = self.create_question(first=True)
 
=== modified file 'madgraph/interface/madevent_interface.py'
--- madgraph/interface/madevent_interface.py	2018-06-08 13:35:22 +0000
+++ madgraph/interface/madevent_interface.py	2018-11-09 10:32:12 +0000
@@ -4142,9 +4142,14 @@
             self.configure_directory(html_opening =False)
             self.check_pythia8(args)        
 
+        # Update the banner with the pythia card
+        if not self.banner or len(self.banner) <=1:
+            # Here the level keyword 'pythia' must not be changed to 'pythia8'.
+            self.banner = banner_mod.recover_banner(self.results, 'pythia')
+
         # the args are modify and the last arg is always the mode 
         if not no_default:
-            self.ask_pythia_run_configuration(args[-1], pythia_version=8)
+            self.ask_pythia_run_configuration(args[-1], pythia_version=8, banner=self.banner)
 
         if self.options['automatic_html_opening']:
             misc.open_file(os.path.join(self.me_dir, 'crossx.html'))
@@ -4157,10 +4162,7 @@
              #"The normalisation of the hepmc output file will be wrong (i.e. non-standard).\n"+\
              #"Please use 'event_norm = average' in the run_card to avoid this problem.")
 
-        # Update the banner with the pythia card
-        if not self.banner or len(self.banner) <=1:
-            # Here the level keyword 'pythia' must not be changed to 'pythia8'.
-            self.banner = banner_mod.recover_banner(self.results, 'pythia')
+
         
         if not self.options['mg5amc_py8_interface_path'] or not \
              os.path.exists(pjoin(self.options['mg5amc_py8_interface_path'],
@@ -6126,7 +6128,7 @@
         return switch
     
     ############################################################################
-    def ask_pythia_run_configuration(self, mode=None, pythia_version=6):
+    def ask_pythia_run_configuration(self, mode=None, pythia_version=6, banner=None):
         """Ask the question when launching pythia"""
         
         pythia_suffix = '' if pythia_version==6 else '%d'%pythia_version
@@ -6187,10 +6189,16 @@
         if self.force:
             return mode
         
+        if not banner:
+            banner = self.banner
+        
         if auto:
-            self.ask_edit_cards(cards, mode='auto', plot=(pythia_version==6))
+            self.ask_edit_cards(cards, from_banner=['param', 'run'], 
+                                mode='auto', plot=(pythia_version==6), banner=banner
+                                )
         else:
-            self.ask_edit_cards(cards, plot=(pythia_version==6))
+            self.ask_edit_cards(cards, from_banner=['param', 'run'],
+                                 plot=(pythia_version==6), banner=banner)
 
         return mode
                 
 
=== modified file 'madgraph/interface/madgraph_interface.py'
--- madgraph/interface/madgraph_interface.py	2018-06-22 06:52:02 +0000
+++ madgraph/interface/madgraph_interface.py	2018-11-09 10:32:12 +0000
@@ -2782,7 +2782,8 @@
                    'gauge','lorentz', 'brs', 'cms']
     _import_formats = ['model_v4', 'model', 'proc_v4', 'command', 'banner']
     _install_opts = ['Delphes', 'MadAnalysis4', 'ExRootAnalysis',
-                     'update', 'Golem95', 'PJFry', 'QCDLoop', 'maddm', 'maddump']
+                     'update', 'Golem95', 'PJFry', 'QCDLoop', 'maddm', 'maddump',
+                     'looptools']
     
     # The targets below are installed using the HEPToolsInstaller.py script
     _advanced_install_opts = ['pythia8','zlib','boost','lhapdf6','lhapdf5','collier',
@@ -5840,7 +5841,7 @@
         function will overwrite any existing installation of the tool without 
         warnings.
         """
-        
+
         # Make sure to avoid any border effect on custom_additional_options
         add_options = list(additional_options)
         
@@ -5857,6 +5858,10 @@
         if args[0] == 'update':
             self.install_update(['update']+install_options['update_options'],wget=program)
             return
+        elif args[0] == 'looptools':
+            self.install_reduction_library(force=True)
+            return
+        
 
         plugin = self.install_plugin
         
@@ -5903,12 +5908,20 @@
                     data = urllib.urlopen(cluster_path)
                 except Exception:
                     continue
+                if data.getcode() != 200:
+                    continue
+                
                 break
+                
             else:
                 raise MadGraph5Error, '''Impossible to connect any of us servers.
                 Please check your internet connection or retry later'''
-            for line in data:
-                split = line.split()
+            for wwwline in data:
+                split = wwwline.split()
+                if len(split)!=2:
+                    if '--source' not in line:
+                        source = {0:'uiuc',1:'ucl'}[index]
+                        return self.do_install(line+' --source='+source, paths=paths, additional_options=additional_options)
                 path[split[0]] = split[1]
 
 ################################################################################
@@ -7037,7 +7050,6 @@
             param_card = check_param_card.ParamCard(out_path.getvalue().split('\n'))
             
             for (block, lhacode) in put_to_one:
-                misc.sprint(block, lhacode)
                 try:
                     param_card[block].get(lhacode).value = 1
                 except:
 
=== modified file 'madgraph/interface/reweight_interface.py'
--- madgraph/interface/reweight_interface.py	2018-06-15 15:10:56 +0000
+++ madgraph/interface/reweight_interface.py	2018-11-09 10:32:12 +0000
@@ -233,8 +233,8 @@
             
             if '=' in order:
                 # get the type NLO QCD/QED/...
-                order = order.split('=',1)[1]
-                
+                order = order.split('=',1)[1].strip()
+
             # define the list of particles that are needed for the radiation
             pert = fks_common.find_pert_particles_interactions(model,
                                            pert_order = order)['soft_particles']
@@ -746,7 +746,8 @@
         # Find new tag in the banner and add information if needed
         if 'initrwgt' in self.banner and self.output_type == 'default': 
             if 'name=\'mg_reweighting\'' in self.banner['initrwgt']:
-                blockpat = re.compile(r'''<weightgroup name=\'mg_reweighting\'\s*>(?P<text>.*?)</weightgroup>''', re.I+re.M+re.S)
+                blockpat = re.compile(r'''<weightgroup name=\'mg_reweighting\'\s*weight_name_strategy=\'includeIdInWeightName\'>(?P<text>.*?)</weightgroup>''', re.I+re.M+re.S)
+                misc.sprint(blockpat, self.banner['initrwgt'])
                 before, content, after = blockpat.split(self.banner['initrwgt'])
                 header_rwgt_other = before + after
                 pattern = re.compile('<weight id=\'(?:rwgt_(?P<id>\d+)|(?P<id2>[_\w]+))(?P<rwgttype>\s*|_\w+)\'>(?P<info>.*?)</weight>', re.S+re.I+re.M)
@@ -1787,7 +1788,7 @@
                         for i in range(len(pdg)):
                             if pdg[i] == oldpdg[i]:
                                 continue
-                            if not self.model or not getattr(self.model, 'get_mass'):
+                            if not self.model or not hasattr(self.model, 'get_mass'):
                                 continue
                             if self.model.get_mass(int(pdg[i])) == self.model.get_mass(int(oldpdg[i])):
                                 continue
 
=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2018-06-07 13:31:42 +0000
+++ madgraph/iolibs/export_v4.py	2018-11-09 10:32:12 +0000
@@ -12,6 +12,7 @@
 # For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
 #
 ################################################################################
+from madgraph.iolibs.helas_call_writers import HelasCallWriter
 """Methods and classes to export matrix elements to v4 format."""
 
 import copy
@@ -93,12 +94,28 @@
     #    - None, madgraph do nothing for initialisation
     exporter = 'v4'
     # language of the output 'v4' for Fortran output
-    #                        'cpp' for C++ output 
+    #                        'cpp' for C++ output
     
     
     def __init__(self, dir_path = "", opt=None):
-        # cmd_options is a dictionary with all the optional argurment passed at output time 
-        return
+        # cmd_options is a dictionary with all the optional argurment passed at output time
+        
+        # Activate some monkey patching for the helas call writer.
+        helas_call_writers.HelasCallWriter.customize_argument_for_all_other_helas_object = \
+                self.helas_call_writer_custom
+        
+
+    # helper function for customise helas writter
+    @staticmethod
+    def custom_helas_call(call, arg):
+        """static method to customise the way aloha function call are written
+        call is the default template for the call
+        arg are the dictionary used for the call
+        """
+        return call, arg
+    
+    helas_call_writer_custom = lambda x,y,z: x.custom_helas_call(y,z)
+
 
     def copy_template(self, model):
         return
@@ -160,6 +177,9 @@
         #place holder to pass information to the run_interface
         self.proc_characteristic = banner_mod.ProcCharacteristic()
         
+        # call mother class
+        super(ProcessExporterFortran,self).__init__(dir_path, opt)
+        
         
     #===========================================================================
     # process exporter fortran switch between group and not grouped
@@ -1834,13 +1854,25 @@
         is_clang = misc.detect_if_cpp_compiler_is_clang(compiler)
         is_lc    = misc.detect_cpp_std_lib_dependence(compiler) == '-lc++'
 
+
         # list of the variable to set in the make_opts file
         for_update= {'DEFAULT_CPP_COMPILER':compiler,
                      'MACFLAG':'-mmacosx-version-min=10.7' if is_clang and is_lc else '',
                      'STDLIB': '-lc++' if is_lc else '-lstdc++',
                      'STDLIB_FLAG': '-stdlib=libc++' if is_lc and is_clang else ''
                      }
-        
+
+        # for MOJAVE remove the MACFLAG:
+        if is_clang:
+            import platform
+            version, _, _ = platform.mac_ver()
+            if not version:# not linux 
+                version = 14 # set version to remove MACFLAG
+            else:
+                version = int(version.split('.')[1])
+            if version >= 14:
+                for_update['MACFLAG'] = '-mmacosx-version-min=10.8' if is_lc else ''
+
         if not root_dir:
             root_dir = self.dir_path
         make_opts = pjoin(root_dir, 'Source', 'make_opts')
@@ -3411,6 +3443,7 @@
         return s_and_t_channels
 
 
+
 #===============================================================================
 # ProcessExporterFortranME
 #===============================================================================
@@ -3419,7 +3452,16 @@
     MadEvent format."""
 
     matrix_file = "matrix_madevent_v4.inc"
-
+    
+    # helper function for customise helas writter
+    @staticmethod
+    def custom_helas_call(call, arg):
+        if arg['mass'] == '%(M)s,%(W)s,':
+            arg['mass'] = '%(M)s, fk_%(W)s,'
+        elif '%(W)s' in arg['mass']:
+            raise Exception
+        return call, arg
+    
     def copy_template(self, model):
         """Additional actions needed for setup of Template
         """
@@ -3443,6 +3485,7 @@
         
         
 
+    
 
 
     #===========================================================================
@@ -3933,10 +3976,29 @@
         # Extract helas calls
         helas_calls = fortran_model.get_matrix_element_calls(\
                     matrix_element)
+        
 
         replace_dict['helas_calls'] = "\n".join(helas_calls)
 
 
+        #adding the support for the fake width (forbidding too small width)
+        mass_width = matrix_element.get_all_mass_widths()
+        width_list = set([e[1] for e in mass_width])
+        
+        replace_dict['fake_width_declaration'] = \
+            ('  double precision fk_%s \n' * len(width_list)) % tuple(width_list)
+        replace_dict['fake_width_declaration'] += \
+            ('  save fk_%s \n' * len(width_list)) % tuple(width_list)
+        fk_w_defs = []
+        one_def = ' fk_%(w)s = SIGN(MAX(ABS(%(w)s), ABS(%(m)s*small_width_treatment)), %(w)s)'     
+        for m, w in mass_width:
+            if w == 'zero':
+                if ' fk_zero = 0d0' not in fk_w_defs: 
+                    fk_w_defs.append(' fk_zero = 0d0')
+                continue    
+            fk_w_defs.append(one_def %{'m':m, 'w':w})
+        replace_dict['fake_width_definitions'] = '\n'.join(fk_w_defs)
+
         # Extract version number and date from VERSION file
         info_lines = self.get_mg5_info_lines()
         replace_dict['info_lines'] = info_lines
 
=== modified file 'madgraph/iolibs/helas_call_writers.py'
--- madgraph/iolibs/helas_call_writers.py	2016-03-03 23:17:06 +0000
+++ madgraph/iolibs/helas_call_writers.py	2018-11-09 10:32:12 +0000
@@ -21,11 +21,13 @@
 import aloha.aloha_writers as aloha_writers
 import aloha
 from madgraph import MadGraph5Error
+import madgraph.various.misc as misc
 
 class HelasWriterError(Exception):
     """Class for the error of this module """
     pass
 
+
 #===============================================================================
 # HelasCallWriter
 #===============================================================================
@@ -41,6 +43,19 @@
     # Dictionaries from spin states to letters in Helas call
     mother_dict = {1: 'S', 2: 'O', -2: 'I', 3: 'V', 5: 'T', 4:'OR', -4:'IR',
                    99:'P'}
+    
+    @staticmethod
+    def customize_argument_for_all_other_helas_object(call, arg):
+        """ Place holder for PLUGIN/...
+            (used by madevent output for small width handling) 
+        """
+        return call, arg
+        
+    @staticmethod
+    def default_customize_argument_for_all_other_helas_object(call,arg):
+        return call,arg  
+    #customize_argument_for_all_other_helas_object = fct_customize_argument_for_all_other_helas_object
+    #default_customize_argument_for_all_other_helas_object = fct_customize_argument_for_all_other_helas_object
 
     def default_setup(self):
 
@@ -85,6 +100,11 @@
 
         return ['model', 'wavefunctions', 'amplitudes']
 
+
+
+    
+     
+
     def get_loop_amp_helas_calls(self, matrix_element):
         """Return a list of strings, corresponding to the Helas calls
         for building loop amplitudes (AMPL) only."""
@@ -1197,6 +1217,7 @@
                  arg['second_line'] = ampl+"="+ampl+"*(%(uvct)s)"           
 
         # ALL ARGUMENT FORMATTED ###############################################
+        call, arg = HelasCallWriter.customize_argument_for_all_other_helas_object(call, arg)
         # Store the result.
         call = call % arg
         # Now we have a line correctly formatted
@@ -1209,6 +1230,8 @@
         else:
             self.add_amplitude(argument.get_call_key(), call_function)
 
+
+
     def get_loop_amplitude_helas_calls(self, loop_matrix_element):
         """ Returns a list of strings corresponding to the Helas calls for each 
         loop amplitude of this loop matrix element. This function is placed in
 
=== modified file 'madgraph/iolibs/template_files/matrix_loop_induced_madevent_group.inc'
--- madgraph/iolibs/template_files/matrix_loop_induced_madevent_group.inc	2017-11-01 10:08:52 +0000
+++ madgraph/iolibs/template_files/matrix_loop_induced_madevent_group.inc	2018-11-09 10:32:12 +0000
@@ -39,6 +39,7 @@
 C   When set negative, the security above is removed
 	DOUBLE PRECISION MULTICHANNEL_THRES
 	PARAMETER (MULTICHANNEL_THRES=1.0d-5)
+
 c
 c global (due to reading writting) 
 c
 
=== modified file 'madgraph/iolibs/template_files/matrix_madevent_group_v4.inc'
--- madgraph/iolibs/template_files/matrix_madevent_group_v4.inc	2017-11-01 10:08:52 +0000
+++ madgraph/iolibs/template_files/matrix_madevent_group_v4.inc	2018-11-09 10:32:12 +0000
@@ -68,6 +68,9 @@
  
     REAL*8 POL(2)
     COMMON/TO_POLARIZATION/ POL
+
+	double precision small_width_treatment
+	common/narrow_width/small_width_treatment
     
     INTEGER          ISUM_HEL
     LOGICAL                    MULTI_CHANNEL
@@ -87,6 +90,7 @@
 C ----------
 C BEGIN CODE
 C ----------
+
     NTRY(IMIRROR)=NTRY(IMIRROR)+1
     THIS_NTRY(IMIRROR) = THIS_NTRY(IMIRROR)+1 
     DO I=1,NEXTERNAL
@@ -255,6 +259,11 @@
 C   Needed for v4 models
     COMPLEX*16 DUM0,DUM1
     DATA DUM0, DUM1/(0d0, 0d0), (1d0, 0d0)/
+
+	%(fake_width_declaration)s
+	logical first
+	data first /.true./
+	save first
 C
 C FUNCTION
 C
@@ -265,6 +274,9 @@
     Double Precision amp2(maxamps), jamp2(0:maxflow)
     common/to_amps/  amp2,       jamp2
     include 'coupl.inc'
+
+	double precision small_width_treatment
+	common/narrow_width/small_width_treatment
 C  
 C COLOR DATA
 C  
@@ -272,6 +284,12 @@
 C ----------
 C BEGIN CODE
 C ----------
+if (first) then
+	first=.false.
+	%(fake_width_definitions)s
+endif
+
+
 %(helas_calls)s
 %(jamp_lines)s
 
 
=== modified file 'madgraph/iolibs/template_files/matrix_madevent_v4.inc'
--- madgraph/iolibs/template_files/matrix_madevent_v4.inc	2017-11-01 10:08:52 +0000
+++ madgraph/iolibs/template_files/matrix_madevent_v4.inc	2018-11-09 10:32:12 +0000
@@ -49,6 +49,7 @@
     INTEGER IDUM, NGOOD, IGOOD(NCOMB), JHEL, J, JJ
     REAL     XRAN1
     EXTERNAL XRAN1
+
 C  
 C GLOBAL VARIABLES
 C  
@@ -58,7 +59,7 @@
     CHARACTER*101        HEL_BUFF
     COMMON/TO_HELICITY/  HEL_BUFF
     
-    REAL*8 POL(2)
+	REAL*8 POL(2)
     COMMON/TO_POLARIZATION/ POL
     
     INTEGER          ISUM_HEL
@@ -68,6 +69,7 @@
     DATA IDUM /-1/
     DATA XTRY, XREJ, NGOOD /0,0,0/
     SAVE YFRAC, IGOOD, JHEL
+
 %(helicity_lines)s
 %(den_factor_line)s
 C ----------
@@ -235,6 +237,11 @@
 C   Needed for v4 models
     COMPLEX*16 DUM0,DUM1
     DATA DUM0, DUM1/(0d0, 0d0), (1d0, 0d0)/
+
+	%(fake_width_declaration)s
+	logical first
+	data first /.true./
+	save first
 C
 C FUNCTION
 C
@@ -245,6 +252,9 @@
     Double Precision amp2(maxamps), jamp2(0:maxflow)
     common/to_amps/  amp2,       jamp2
     include 'coupl.inc'
+
+	double precision small_width_treatment
+	common/narrow_width/small_width_treatment
 C  
 C COLOR DATA
 C  
@@ -252,6 +262,11 @@
 C ----------
 C BEGIN CODE
 C ----------
+if (first) then
+	first=.false.
+	%(fake_width_definitions)s
+endif
+
 %(helas_calls)s
 %(jamp_lines)s
 
 
=== modified file 'madgraph/iolibs/template_files/matrix_madweight_group_v4.inc'
--- madgraph/iolibs/template_files/matrix_madweight_group_v4.inc	2012-02-29 12:53:46 +0000
+++ madgraph/iolibs/template_files/matrix_madweight_group_v4.inc	2018-11-09 10:32:12 +0000
@@ -22,6 +22,7 @@
 C ARGUMENTS 
 C  
       REAL*8 P(0:3,NEXTERNAL),ANS
+
 C  
 C LOCAL VARIABLES 
 C  
 
=== modified file 'madgraph/madevent/sum_html.py'
--- madgraph/madevent/sum_html.py	2018-05-25 09:13:24 +0000
+++ madgraph/madevent/sum_html.py	2018-11-09 10:32:12 +0000
@@ -699,7 +699,7 @@
         P_comb = Combine_results(Pdir)
         
         if jobs:
-            for job in filter(lambda j: j['p_dir'] == Pdir, jobs):
+            for job in filter(lambda j: j['p_dir'] in Pdir, jobs):
                     P_comb.add_results(os.path.basename(job['dirname']),\
                                        pjoin(job['dirname'],'results.dat'))
         elif folder_names:
 
=== modified file 'madgraph/various/banner.py'
--- madgraph/various/banner.py	2018-05-17 13:55:48 +0000
+++ madgraph/various/banner.py	2018-11-09 10:32:12 +0000
@@ -483,9 +483,9 @@
     def charge_card(self, tag):
         """Build the python object associated to the card"""
         
-        if tag == 'param_card':
+        if tag in ['param_card', 'param']:
             tag = 'slha'
-        elif tag == 'run_card':
+        elif tag  in ['run_card', 'run']:
             tag = 'mgruncard' 
         elif tag == 'proc_card':
             tag = 'mg5proccard' 
@@ -905,7 +905,10 @@
             out = []
             for line in self:
                 if line.startswith('define'):
-                    name, content = line[7:].split('=',1)
+                    try:
+                        name, content = line[7:].split('=',1)
+                    except ValueError:
+                        name, content = line[7:].split(None,1)
                     out.append((name, content))
             return out 
         else:
@@ -2372,9 +2375,9 @@
                     else:
                         this_group = this_group[0]
                     if block_name in write_block:
-                        text += this_group.on_template % self
+                        text += this_group.template_on % self
                         for name in this_group.fields:
-                            written.add(f)
+                            written.add(name)
                             if name in to_write:
                                 to_write.remove(name)
                     else:
@@ -2949,6 +2952,7 @@
         self.add_param('survey_splitting', -1, hidden=True, include=False, comment="for loop-induced control how many core are used at survey for the computation of a single iteration.")
         self.add_param('survey_nchannel_per_job', 2, hidden=True, include=False, comment="control how many Channel are integrated inside a single job on cluster/multicore")
         self.add_param('refine_evt_by_job', -1, hidden=True, include=False, comment="control the maximal number of events for the first iteration of the refine (larger means less jobs)")
+        self.add_param('small_width_treatment', 1e-6, hidden=True, comment="generation where the width is below VALUE times mass will be replace by VALUE times mass for the computation. The cross-section will be corrected assuming NWA. Not used for loop-induced process")
         
         # parameter allowing to define simple cut via the pdg
         # Special syntax are related to those. (can not be edit directly)
@@ -2970,7 +2974,7 @@
         self.add_param('etamax4pdg',[-1.], system=True)   
         self.add_param('mxxmin4pdg',[-1.], system=True)
         self.add_param('mxxpart_antipart', [False], system=True)
-            
+                     
              
     def check_validity(self):
         """ """
 
=== modified file 'madgraph/various/lhe_parser.py'
--- madgraph/various/lhe_parser.py	2018-06-21 08:57:21 +0000
+++ madgraph/various/lhe_parser.py	2018-11-09 10:32:12 +0000
@@ -167,6 +167,7 @@
         """associate to this particle the decay in the associate event"""
         
         return self.event.add_decay_to_particle(self.event_id, decay_event)
+
             
     def __repr__(self):
         return 'Particle("%s", event=%s)' % (str(self), self.event)
@@ -1263,26 +1264,28 @@
                 
             if 'part' == status:
                 part = Particle(line, event=self)
-                if part.E != 0:
+                if part.E != 0 or part.status==-1:
                     self.append(part)
                 elif self.nexternal:
-                        self.nexternal-=1
+                    self.nexternal-=1
             else:
                 if '</event>' in line:
                     line = line.replace('</event>','',1)
                 self.tag += '%s\n' % line
                 
         self.assign_mother()
-        
+    
+    
     def assign_mother(self):
         """convert the number in actual particle"""
-        #Security if not incoming particle. Define a fake particle and set all particle as 
-        # decaying from that fake particle
+        #Security if not incoming particle. Define a fake particle 
         if all(p.status != -1 for p in self):
+            if not self.nexternal:
+                return
             if self.warning_order:
                 Event.warning_order = False
                 logger.warning("Weird format for lhe format: no incoming particle... adding a fake one")
-            
+                raise Exception
             mother = Particle(event=self)
             mother.status = -1
             mother.pid = 0
@@ -1291,8 +1294,10 @@
             mother.event_id = 0
             self.nexternal += 1
             for p in self[1:]:
-                p.mother1 = 1
-                p.mother2 = 1
+                if isinstance(p.mother1, int) and p.mother1 > 1:
+                    p.mother1 += 1
+                if isinstance(p.mother2, int) and p.mother2 > 1:
+                    p.mother2 += 1
                 p.event_id += 1
             
         
 
=== modified file 'madgraph/various/systematics.py'
--- madgraph/various/systematics.py	2018-03-25 08:30:56 +0000
+++ madgraph/various/systematics.py	2018-11-09 10:32:12 +0000
@@ -885,6 +885,7 @@
             for onewgt in cevent.wgts:
                 if not __debug__ and (dyn== -1 and Dmur==1 and Dmuf==1 and pdf==self.orig_pdf):
                     wgt += onewgt.ref_wgt 
+                    continue
                 
                 if dyn == -1:
                     mur2 = onewgt.scales2[1]
 
=== modified file 'models/check_param_card.py'
--- models/check_param_card.py	2018-04-26 14:24:04 +0000
+++ models/check_param_card.py	2018-11-09 10:32:12 +0000
@@ -178,7 +178,7 @@
     
     def rename_keys(self, change_keys):
         
-        misc.sprint(self.param_dict, change_keys, [p.lhacode for p in self])
+
         for old_key, new_key in change_keys.items():
             
             assert old_key in self.param_dict
@@ -509,15 +509,15 @@
                     if model_value.imag > 1e-5 * model_value.real:
                         raise Exception, "Width should be real number: particle %s (%s) has mass: %s" 
                     model_value = model_value.real
-                if not misc.equal(model_value, param_value, 4):
+                if not misc.equal(abs(model_value), param_value, 4):
                     modify = True
                     if loglevel == 20:
                         logger.info('For consistency, the width of particle %s (%s) is changed to %s.' % (lhacode, particle.get('name'), model_value), '$MG:BOLD')
                     else:
                         logger.log(loglevel,'For consistency, the width of particle %s (%s) is changed to %s.' % (lhacode, particle.get('name'), model_value))
                     #logger.debug('was %s', param_value)
-                if model_value != param_value:   
-                    self.get('decay').get(abs(particle.get_pdg_code())).value = model_value
+                if abs(model_value) != param_value:   
+                    self.get('decay').get(abs(particle.get_pdg_code())).value = abs(model_value)
 
         return modify
 
@@ -640,6 +640,11 @@
                 logger.warning('information about \"%s %s" is missing (full block missing) using default value: %s.' %\
                                    (block, lhaid, value))
             value = str(value).lower()
+            #special handling for negative mass -> set width negative
+            if block == 'decay':
+                if self['mass'].get(tuple(lhaid)).value < 0:
+                    value = '-%s' % value
+
             fout.writelines(' %s = %s' % (variable, ('%e'%float(value)).replace('e','d')))
             if need_mp:
                 fout.writelines(' mp__%s = %s_16' % (variable, value))
@@ -1036,7 +1041,6 @@
                     data.append(info[k])
                 else:
                     data.append(0.)
-            misc.sprint(name, bench, data)
             ff.write(formatting % tuple([name] + bench + data))
                 
         if not path:
 
=== modified file 'models/import_ufo.py'
--- models/import_ufo.py	2018-06-15 13:28:53 +0000
+++ models/import_ufo.py	2018-11-09 10:32:12 +0000
@@ -112,6 +112,8 @@
             data = urllib.urlopen(cluster_path)
         except Exception:
             continue
+        if data.getcode() != 200:
+            continue
         break
     else:
         raise MadGraph5Error, '''Model not found locally and Impossible to connect any of us servers.
@@ -138,7 +140,8 @@
     target = None 
     if 'PYTHONPATH' in os.environ and not local_dir:
         for directory in os.environ['PYTHONPATH'].split(':'):
-            if 'UFO' in os.path.basename(directory) and os.path.exists(directory):
+            if 'UFO' in os.path.basename(directory) and os.path.exists(directory) and\
+                    misc.glob('*/couplings.py', path=directory):
                 target= directory 
     if target is None:
         target = pjoin(MG5DIR, 'models')    
@@ -646,14 +649,15 @@
                 continue
             names = [interaction['lorentz'][i] for i in to_lor[key]]
             names.sort()
-            
+            if self.lorentz_info[names[0]].get('structure') == 'external':
+                continue
             # get name of the new lorentz
             if tuple(names) in self.lorentz_combine:
                 # already created new loretnz
                 new_name = self.lorentz_combine[tuple(names)]
             else:
                 new_name = self.add_merge_lorentz(names)
-                
+
             # remove the old couplings 
             color, coup = key
             to_remove = [(color, lor) for lor in to_lor[key]]  
@@ -696,7 +700,12 @@
         # load the associate lorentz expression
         new_struct = ' + '.join([self.lorentz_info[n].get('structure') for n in names])
         spins = self.lorentz_info[names[0]].get('spins')
-        new_lor = self.add_lorentz(new_name, spins, new_struct)
+        formfactors = sum([ self.lorentz_info[n].get('formfactors') for n in names \
+                            if hasattr(self.lorentz_info[n], 'formfactors') \
+                            and self.lorentz_info[n].get('formfactors') \
+                      ],[])
+                        
+        new_lor = self.add_lorentz(new_name, spins, new_struct, formfactors)
         self.lorentz_info[new_name] = new_lor
         
         return new_name
@@ -1368,15 +1377,19 @@
                     
         return  '' if sign ==1 else '-'
 
-    def add_lorentz(self, name, spins , expr):
+    def add_lorentz(self, name, spins , expr, formfact=None):
         """ Add a Lorentz expression which is not present in the UFO """
 
+        logger.debug('MG5 converter defines %s to %s', name, expr)
         assert name not in [l.name for l in self.model['lorentz']]
         with misc.TMP_variable(self.ufomodel.object_library, 'all_lorentz', 
                                self.model['lorentz']):
             new = self.model['lorentz'][0].__class__(name = name,
                     spins = spins,
                     structure = expr)
+            if formfact:
+                new.formfactors = formfact
+
         assert name in [l.name for l in self.model['lorentz']]
         assert name not in [l.name for l in self.ufomodel.all_lorentz]
         #self.model['lorentz'].append(new) # already done by above command
@@ -1973,7 +1986,7 @@
                 null_parameters.append(name)
             elif value == 1:
                 one_parameters.append(name)
-        
+
         return null_parameters, one_parameters
     
     def apply_conditional_simplifications(self, modified_params,
@@ -2301,7 +2314,9 @@
                 particle['width'] = 'ZERO'
             if particle['width'] in one_parameters:
                 one_parameters.remove(particle['width'])                
-                
+            if particle['mass'] in one_parameters:
+                one_parameters.remove(particle['mass'])                
+
         for pdg, particle in self['particle_dict'].items():
             if particle['mass'] in zero_parameters:
                 particle['mass'] = 'ZERO'
@@ -2340,6 +2355,13 @@
                     for coupling in coupling_list:
                         for use in  re_pat.findall(coupling.expr):
                             used.add(use)
+                
+                # check in form-factor
+                for lor in self['lorentz']:
+                    if hasattr(lor, 'formfactors') and lor.formfactors:
+                        for ff in lor.formfactors:
+                            for use in  re_pat.findall(ff.value):
+                                used.add(use)
         else:
             used = set([i for i in special_parameters if i])
         
@@ -2484,16 +2506,26 @@
         if any( n.startswith('d') for n in names ):
             new_struct += '-' + ' - '.join(['1.*(%s)' %self.lorentz_info[n[1:]].get('structure') for n in names if n.startswith('d')])
         spins = self.lorentz_info[names[0][1:]].get('spins')
-        new_lor = self.add_lorentz(new_name, spins, new_struct)
+        formfact = sum([ self.lorentz_info[n[1:]].get('formfactors') for n in names \
+                            if hasattr(self.lorentz_info[n[1:]], 'formfactors') \
+                              and self.lorentz_info[n[1:]].get('formfactors') \
+                       ],[])
+
+
+
+ 
+        new_lor = self.add_lorentz(new_name, spins, new_struct, formfact)
         self.lorentz_info[new_name] = new_lor
         
         return new_name
     
-    def add_lorentz(self, name, spin, struct):
+    def add_lorentz(self, name, spin, struct, formfact=None):
         """adding lorentz structure to the current model"""
         new = self['lorentz'][0].__class__(name = name,
                                            spins = spin,
                                            structure = struct)
+        if formfact:
+            new.formfactors = formfact
         self['lorentz'].append(new)
         self.create_lorentz_dict()
         
 
=== modified file 'tests/IOTests.py'
--- tests/IOTests.py	2016-09-07 13:39:06 +0000
+++ tests/IOTests.py	2018-11-09 10:32:12 +0000
@@ -296,20 +296,26 @@
         """Clean up the file created. Called at the end of the test run."""
         
         pathsToClean = [self.temporary_folder]
-        
+
         if not self.clean_function is None:
             paths, prevent_cleanUp = self.clean_function(*args, **kwargs)
-            pathsToClean.extend(paths)
+            if isinstance(paths, str):
+                pathsToClean.append(paths)
+            else:
+                pathsToClean.extend(paths)
             if prevent_cleanUp:
                 print colored%(31,
                     "Clean up of the following of temporary folders prevented:")
                 for path in pathsToClean:
                     print colored%(31,"  > %s"%str(path))
-
+       
         try:
             for path in pathsToClean:
-                shutil.rmtree(path)
-        except OSError:
+                if os.path.isdir(path):
+                    shutil.rmtree(path)
+                else:
+                    os.remove(path)
+        except OSError,error:
             pass               
 
 #===============================================================================
 
=== modified file 'tests/acceptance_tests/test_cmd_madloop.py'
--- tests/acceptance_tests/test_cmd_madloop.py	2018-04-24 22:39:50 +0000
+++ tests/acceptance_tests/test_cmd_madloop.py	2018-11-09 10:32:12 +0000
@@ -492,7 +492,7 @@
         # Select the Tensor Integral to include in the test
         misc.deactivate_dependence('pjfry', cmd = self.interface, log='stdout')
         misc.deactivate_dependence('samurai', cmd = self.interface, log='stdout')        
-        misc.activate_dependence('golem', cmd = self.interface, log='stdout')
+        misc.deactivate_dependence('golem', cmd = self.interface, log='stdout')
         misc.activate_dependence('ninja', cmd = self.interface, log='stdout',MG5dir=MG5DIR)
         misc.activate_dependence('collier', cmd = self.interface, log='stdout',MG5dir=MG5DIR)
 
@@ -533,7 +533,7 @@
         # Select the Tensor Integral to include in the test
         misc.deactivate_dependence('pjfry', cmd = interface, log='stdout')
         misc.deactivate_dependence('samurai', cmd = interface, log='stdout')        
-        misc.activate_dependence('golem', cmd = interface, log='stdout')
+        misc.deactivate_dependence('golem', cmd = interface, log='stdout')
         misc.activate_dependence('ninja', cmd = interface, log='stdout',MG5dir=MG5DIR)
 
         run_cmd('generate g g > t t~ [virt=QCD]')
 
=== modified file 'tests/acceptance_tests/test_export_fks.py'
--- tests/acceptance_tests/test_export_fks.py	2017-05-15 10:39:31 +0000
+++ tests/acceptance_tests/test_export_fks.py	2018-11-09 10:32:12 +0000
@@ -75,7 +75,7 @@
         # Select the Tensor Integral to include in the test
         misc.deactivate_dependence('pjfry', cmd = interface, log='stdout')
         misc.deactivate_dependence('samurai', cmd = interface, log='stdout')        
-        misc.activate_dependence('golem', cmd = interface, log='stdout')
+        misc.deactivate_dependence('golem', cmd = interface, log='stdout')
         misc.activate_dependence('ninja', cmd = interface, log='stdout',MG5dir=MG5DIR)        
         
         run_cmd('import model %s' % model)
 
=== modified file 'tests/acceptance_tests/test_madspin.py'
--- tests/acceptance_tests/test_madspin.py	2018-05-25 09:13:24 +0000
+++ tests/acceptance_tests/test_madspin.py	2018-11-09 10:32:12 +0000
@@ -93,7 +93,7 @@
         subprocess.call([pjoin(MG5DIR, 'MadSpin', 'madspin'),
                          pjoin(self.path, 'test_hepmc')],
                         cwd=pjoin(self.path),
-                        stdout=stdout,stderr=stdout)
+                        stdout=stdout,stderr=stderr)
         self.assertTrue(os.path.exists(pjoin(self.path, 'test_decayed.lhe.gz')))
         lhe = lhe_parser.EventFile(pjoin(self.path, 'test_decayed.lhe.gz'))
         self.assertEqual(10, len(lhe))
@@ -111,4 +111,56 @@
                     
         self.assertEqual(nb_dec, 116)
         self.assertEqual(nb_photon, 116)
-        
 \ No newline at end of file
+
+    def test_lhe_none_decay(self):
+        """ """
+        
+        cwd = os.getcwd()
+        
+        files.cp(pjoin(MG5DIR, 'tests', 'input_files', 'test_spinmode_none.lhe.gz'), self.path)
+
+
+        fsock = open(pjoin(self.path, 'test_hepmc'),'w')
+        text = """
+        set spinmode none
+        import ./test_spinmode_none.lhe.gz
+        decay z > mu+ mu-
+        launch
+        """
+        
+        fsock.write(text)
+        fsock.close()
+
+        import subprocess
+        if logging.getLogger('madgraph').level <= 20:
+            stdout=None
+            stderr=None
+        else:
+            devnull =open(os.devnull,'w')
+            stdout=devnull
+            stderr=devnull
+
+        subprocess.call([pjoin(MG5DIR, 'MadSpin', 'madspin'),
+                         pjoin(self.path, 'test_hepmc')],
+                        cwd=pjoin(self.path),
+                        stdout=stdout,stderr=stderr)
+
+        self.assertTrue(os.path.exists(pjoin(self.path, 'test_spinmode_none_decayed.lhe.gz')))
+        lhe = lhe_parser.EventFile(pjoin(self.path, 'test_spinmode_none_decayed.lhe.gz'))
+        self.assertEqual(100, len(lhe))
+        
+        nb_dec = 0
+        nb_muon = 0
+        for event in lhe:
+            muon_in = 0
+            self.assertEqual(event.nexternal, len(event))
+            for particle in event:
+                if particle.pdg == 23:
+                    self.assertEqual(particle.status,2)
+                    nb_dec += 1
+                if particle.pdg == 13:
+                    nb_muon += 1
+                    muon_in +=1
+            self.assertEqual(muon_in, 1)
+        self.assertEqual(nb_dec, 189)
+        self.assertEqual(nb_muon, 100)
 \ No newline at end of file
 
=== modified file 'tests/acceptance_tests/test_model_equivalence.py'
--- tests/acceptance_tests/test_model_equivalence.py	2018-06-15 08:18:31 +0000
+++ tests/acceptance_tests/test_model_equivalence.py	2018-11-09 10:32:12 +0000
@@ -146,15 +146,14 @@
         ufo_model.pass_particles_name_in_mg_default()
         
         # import MG4 model
+        
         model = base_objects.Model()
         if not MG4DIR:
             raise MadGraph5Error, "Please provide a valid MG/ME path with -d"
         v4_path = os.path.join(MG4DIR, 'models', 'mssm_v4')
         if not os.path.isdir(v4_path):
-            v4_path = os.path.join(MG4DIR, 'Models', 'mssm')
-            if not os.path.isdir(v4_path):
-                raise MadGraph5Error, \
-                      "Please provide a valid MG/ME path with -d"
+            import_ufo.import_model_from_db('mssm_v4', local_dir=True)
+                
 
         model.set('particles', files.read_from_file(
                os.path.join(v4_path,'particles.dat'),
 
=== removed file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%GOLEM_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%GOLEM_interface.f	2016-08-07 07:16:07 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%GOLEM_interface.f	1970-01-01 00:00:00 +0000
@@ -1,748 +0,0 @@
-      SUBROUTINE GOLEMLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
-C     
-C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
-C     By the MadGraph5_aMC@NLO Development Team
-C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
-C     
-C     Interface between MG5 and Golem95.
-C     The Golem95 version should be higher than 1.3.0.
-C     It supports RANK = NLOOPLINE + 1 tensor integrals when 1 <
-C      NLOOPLINE < 6.
-C     
-C     Process: d~ u > w+ QED<=1 WEIGHTED<=2 [ all = QCD ]
-C     Process: s~ c > w+ QED<=1 WEIGHTED<=2 [ all = QCD ]
-C     
-C     
-C     MODULES
-C     
-      USE MATRICE_S
-      USE FORM_FACTOR_TYPE, ONLY: FORM_FACTOR
-      USE PRECISION_GOLEM, ONLY: KI
-      USE TENS_COMB
-      USE TENS_REC
-      USE FORM_FACTOR_1P, ONLY: A10
-      USE FORM_FACTOR_2P, ONLY: A20
-      USE FORM_FACTOR_3P, ONLY: A30
-      USE FORM_FACTOR_4P, ONLY: A40
-      USE FORM_FACTOR_5P, ONLY: A50
-      USE FORM_FACTOR_6P, ONLY: A60
-      IMPLICIT NONE
-C     
-C     CONSTANTS 
-C     
-      INTEGER    NEXTERNAL
-      PARAMETER (NEXTERNAL=3)
-      LOGICAL CHECKPCONSERVATION
-      PARAMETER (CHECKPCONSERVATION=.TRUE.)
-      REAL*8 NORMALIZATION
-      PARAMETER (NORMALIZATION = 1.D0/(16.D0*3.14159265358979323846D0*
-     $ *2))
-      REAL(KI),DIMENSION(0:3),PARAMETER::NULL_VEC = (/0.0_KI,0.0_KI
-     $ ,0.0_KI,0.0_KI/)
-C     GOLEM_RUN_MODE = 1: Use directly MadLoop tensorial coefficients
-C     GOLEM_RUN_MODE = 2: Reconstruct the tensorial coefficeints
-C      directly from 
-C     numerator using golem internal reconstruction routine
-C     GOLEM_RUN_MODE = 3: Cross-checked reconstructed coefficients
-C      against
-C     MadLoop internal ones.
-      INTEGER GOLEM_RUN_MODE
-      PARAMETER (GOLEM_RUN_MODE=1)
-C     The following is the acceptance threshold used for
-C      GOLEM_RUN_MODE = 3
-      REAL*8 COEF_CHECK_THRS
-      DATA COEF_CHECK_THRS/1.0D-13/
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-
-      LOGICAL PASS_COEF_CHECK
-C     
-C     ARGUMENTS 
-C     
-      INTEGER NLOOPLINE, RANK
-      REAL*8 PL(0:3,NLOOPLINE)
-      REAL*8 PCT(0:3,0:NLOOPLINE-1), ABSPCT(0:3)
-      REAL*8 REF_P
-      REAL(KI) PGOLEM(NLOOPLINE,0:3)
-      COMPLEX*16 M2L(NLOOPLINE)
-      COMPLEX(KI) M2LGOLEM(NLOOPLINE)
-      COMPLEX*16 RES(3)
-      LOGICAL STABLE
-C     
-C     LOCAL VARIABLES 
-C     
-      INTEGER I, J, K
-      TYPE(FORM_FACTOR)::RES_GOLEM
-
-      COMPLEX(KI)::COEFFS0,COEFFS0_REC
-      TYPE(COEFF_TYPE_1)::COEFFS1,COEFFS1_REC
-      TYPE(COEFF_TYPE_2)::COEFFS2,COEFFS2_REC
-      TYPE(COEFF_TYPE_3)::COEFFS3,COEFFS3_REC
-      TYPE(COEFF_TYPE_4)::COEFFS4,COEFFS4_REC
-      TYPE(COEFF_TYPE_5)::COEFFS5,COEFFS5_REC
-      TYPE(COEFF_TYPE_6)::COEFFS6,COEFFS6_REC
-
-C     The pinch propagator optimization is not used, so for now it is
-C     always 0.
-      INTEGER PINCH
-C     
-C     EXTERNAL FUNCTIONS
-C     
-      COMPLEX(KI) GOLEM_LOOPNUM
-      EXTERNAL GOLEM_LOOPNUM
-      LOGICAL COMPARE_COEFS_0
-      LOGICAL COMPARE_COEFS_1
-      LOGICAL COMPARE_COEFS_2
-      LOGICAL COMPARE_COEFS_3
-      LOGICAL COMPARE_COEFS_4
-      LOGICAL COMPARE_COEFS_5
-      LOGICAL COMPARE_COEFS_6
-C     
-C     GLOBAL VARIABLES
-C     
-      INCLUDE 'coupl.inc'
-      INTEGER CTMODE
-      REAL*8 LSCALE
-      COMMON/CT/LSCALE,CTMODE
-
-      INTEGER ID,SQSOINDEX,R
-      COMMON/LOOP/ID,SQSOINDEX,R
-
-      LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
-     $ ,COLLIERINIT
-      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT,GOLEMINIT,SAMURAIINIT
-     $ ,NINJAINIT,COLLIERINIT
-
-      INTEGER NLOOPGROUPS
-      PARAMETER (NLOOPGROUPS=1)
-      INTEGER NSQUAREDSO
-      PARAMETER (NSQUAREDSO=1)
-      INCLUDE 'loop_max_coefs.inc'
-
-      COMPLEX*16 LOOPCOEFS(0:LOOPMAXCOEFS-1,NSQUAREDSO,NLOOPGROUPS)
-      COMMON/LCOEFS/LOOPCOEFS
-C     ----------
-C     BEGIN CODE
-C     ----------
-
-C     The CT initialization is also performed here if not done already
-C      because it calls MPINIT of OneLOop which is necessary on some
-C      system
-      IF (CTINIT) THEN
-        CTINIT=.FALSE.
-        CALL INITCT()
-      ENDIF
-
-C     INITIALIZE GOLEM IF NEEDED
-      IF (GOLEMINIT) THEN
-        GOLEMINIT=.FALSE.
-        CALL INITGOLEM()
-      ENDIF
-
-C     No stability test intrisic to Golem95 now
-      STABLE=.TRUE.
-
-C     This initialization must be done for each reduction because we
-C     have not setup anyoptimization using pinched propagators yet.
-      CALL INITGOLEM95(NLOOPLINE)
-      PINCH = 0
-
-C     YOU CAN FIND THE DETAILS ABOUT THE DIFFERENT CTMODE AT THE
-C      BEGINNING OF THE FILE CTS_CUTS.F90 IN THE CUTTOOLS DISTRIBUTION
-
-C     CONVERT THE MASSES TO BE COMPLEX
-      DO I=1,NLOOPLINE
-        M2LGOLEM(I)=M2L(I)
-      ENDDO
-
-C     CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO CT CONVENTIONS
-      DO I=0,3
-        ABSPCT(I) = 0.D0
-        DO J=0,(NLOOPLINE-1)
-          PCT(I,J)=0.D0
-        ENDDO
-      ENDDO
-      DO I=0,3
-        DO J=1,NLOOPLINE
-          PCT(I,0)=PCT(I,0)+PL(I,J)
-          ABSPCT(I)=ABSPCT(I)+ABS(PL(I,J))
-        ENDDO
-      ENDDO
-      REF_P = MAX(ABSPCT(0), ABSPCT(1),ABSPCT(2),ABSPCT(3))
-      DO I=0,3
-        ABSPCT(I) = MAX(REF_P*1E-6, ABSPCT(I))
-      ENDDO
-      IF (CHECKPCONSERVATION.AND.REF_P.GT.1D-8) THEN
-        IF ((PCT(0,0)/ABSPCT(0)).GT.1.D-6) THEN
-          WRITE(*,*) 'energy is not conserved ',PCT(0,0)
-          STOP 'energy is not conserved'
-        ELSEIF ((PCT(1,0)/ABSPCT(1)).GT.1.D-6) THEN
-          WRITE(*,*) 'px is not conserved ',PCT(1,0)
-          STOP 'px is not conserved'
-        ELSEIF ((PCT(2,0)/ABSPCT(2)).GT.1.D-6) THEN
-          WRITE(*,*) 'py is not conserved ',PCT(2,0)
-          STOP 'py is not conserved'
-        ELSEIF ((PCT(3,0)/ABSPCT(3)).GT.1.D-6) THEN
-          WRITE(*,*) 'pz is not conserved ',PCT(3,0)
-          STOP 'pz is not conserved'
-        ENDIF
-      ENDIF
-      DO I=0,3
-        DO J=1,(NLOOPLINE-1)
-          DO K=1,J
-            PCT(I,J)=PCT(I,J)+PL(I,K)
-          ENDDO
-        ENDDO
-      ENDDO
-
-C     Now convert the loop momenta to Golem95 conventions
-      DO I=0,3
-        PGOLEM(1,I)=0.0E0_KI
-        DO J=2,NLOOPLINE
-          PGOLEM(J,I)=PCT(I,J-1)
-        ENDDO
-      ENDDO
-
-C     Fill in the kinematic s-matrix while taking care of on-shell
-C      limits.
-      CALL SETUP_KIN_MATRIX(NLOOPLINE,PGOLEM,M2LGOLEM)
-C     Construct the golem internal matrices derived from the kinetic
-C      one.
-      CALL PREPARESMATRIX()
-
-C     Fill in the golem coefficents and compute the loop
-      IF(GOLEM_RUN_MODE.EQ.2)THEN
-        RES_GOLEM = EVALUATE_B(GOLEM_LOOPNUM,PGOLEM,0,RANK)
-      ELSE
-        PASS_COEF_CHECK=.TRUE.
-        SELECT CASE(RANK)
-        CASE(0)
-        CALL FILL_GOLEM_COEFFS_0(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS0)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          COEFFS0_REC = GOLEM_LOOPNUM(NULL_VEC,0.0_KI)
-          PASS_COEF_CHECK=COMPARE_COEFS_0(COEFFS0,COEFFS0_REC)
-        ENDIF
-        CASE(1)
-        CALL FILL_GOLEM_COEFFS_1(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS1)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT1(GOLEM_LOOPNUM,COEFFS1_REC)
-          PASS_COEF_CHECK=COMPARE_COEFS_1(COEFFS1,COEFFS1_REC)
-        ENDIF
-        CASE(2)
-        CALL FILL_GOLEM_COEFFS_2(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS2)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT2(GOLEM_LOOPNUM,COEFFS2_REC)
-          PASS_COEF_CHECK=COMPARE_COEFS_2(COEFFS2,COEFFS2_REC)
-        ENDIF
-        CASE(3)
-        CALL FILL_GOLEM_COEFFS_3(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS3)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT3(GOLEM_LOOPNUM,COEFFS3_REC)
-          PASS_COEF_CHECK=COMPARE_COEFS_3(COEFFS3,COEFFS3_REC)
-        ENDIF
-        CASE(4)
-        CALL FILL_GOLEM_COEFFS_4(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS4)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT4(GOLEM_LOOPNUM,COEFFS4_REC)
-          PASS_COEF_CHECK=COMPARE_COEFS_4(COEFFS4,COEFFS4_REC)
-        ENDIF
-        CASE(5)
-        CALL FILL_GOLEM_COEFFS_5(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS5)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT5(GOLEM_LOOPNUM,COEFFS5_REC)
-          PASS_COEF_CHECK=COMPARE_COEFS_5(COEFFS5,COEFFS5_REC)
-        ENDIF
-        CASE(6)
-        CALL FILL_GOLEM_COEFFS_6(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS6)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT6(GOLEM_LOOPNUM,COEFFS6_REC)
-          PASS_COEF_CHECK=COMPARE_COEFS_6(COEFFS6,COEFFS6_REC)
-        ENDIF
-        CASE DEFAULT
-        WRITE(*,*)'Not yet implemented in Golem95 for rank= ',RANK
-        STOP
-        END SELECT
-
-        IF(.NOT.PASS_COEF_CHECK)THEN
-          WRITE(*,*)'Coefs mismatch for ID ',ID,' and rank ',RANK
-          WRITE(*,*)'Coefs form MadLoop5:'
-          SELECT CASE(RANK)
-          CASE(0)
-          WRITE(*,*)'Constant coef = ',COEFFS0
-          CASE(1)
-          CALL PRINT_COEFFS(COEFFS1)
-          CASE(2)
-          CALL PRINT_COEFFS(COEFFS2)
-          CASE(3)
-          CALL PRINT_COEFFS(COEFFS3)
-          CASE(4)
-          CALL PRINT_COEFFS(COEFFS4)
-          CASE(5)
-          CALL PRINT_COEFFS(COEFFS5)
-          CASE(6)
-          CALL PRINT_COEFFS(COEFFS6)
-          END SELECT
-          WRITE(*,*)'Coefs reconstructed by Golem95:'
-          SELECT CASE(RANK)
-          CASE(0)
-          WRITE(*,*)'Constant coef = ',COEFFS0_REC
-          CASE(1)
-          CALL PRINT_COEFFS(COEFFS1_REC)
-          CASE(2)
-          CALL PRINT_COEFFS(COEFFS2_REC)
-          CASE(3)
-          CALL PRINT_COEFFS(COEFFS3_REC)
-          CASE(4)
-          CALL PRINT_COEFFS(COEFFS4_REC)
-          CASE(5)
-          CALL PRINT_COEFFS(COEFFS5_REC)
-          CASE(6)
-          CALL PRINT_COEFFS(COEFFS6_REC)
-          END SELECT
-          STOP
-        ENDIF
-
-        SELECT CASE(NLOOPLINE)
-        CASE(1)
-        WRITE(*,*)'Golem95 cannot handle with tadpole yet'
-        STOP
-        CASE(2)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A20(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT2_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT2_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT2_3(COEFFS3,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',2,RANK
-        STOP
-        END SELECT
-        CASE(3)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A30(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT3_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT3_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT3_3(COEFFS3,PGOLEM,PINCH)
-        CASE(4)
-        RES_GOLEM = CONTRACT3_4(COEFFS4,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',3,RANK
-        STOP
-        END SELECT
-        CASE(4)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A40(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT4_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT4_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT4_3(COEFFS3,PGOLEM,PINCH)
-        CASE(4)
-        RES_GOLEM = CONTRACT4_4(COEFFS4,PGOLEM,PINCH)
-        CASE(5)
-        RES_GOLEM = CONTRACT4_5(COEFFS5,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',4,RANK
-        STOP
-        END SELECT
-        CASE(5)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A50(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT5_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT5_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT5_3(COEFFS3,PGOLEM,PINCH)
-        CASE(4)
-        RES_GOLEM = CONTRACT5_4(COEFFS4,PGOLEM,PINCH)
-        CASE(5)
-        RES_GOLEM = CONTRACT5_5(COEFFS5,PGOLEM,PINCH)
-        CASE(6)
-        RES_GOLEM = CONTRACT5_6(COEFFS6,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',5,RANK
-        STOP
-        END SELECT
-        CASE(6)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A60(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT6_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT6_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT6_3(COEFFS3,PGOLEM,PINCH)
-        CASE(4)
-        RES_GOLEM = CONTRACT6_4(COEFFS4,PGOLEM,PINCH)
-        CASE(5)
-        RES_GOLEM = CONTRACT6_5(COEFFS5,PGOLEM,PINCH)
-        CASE(6)
-        RES_GOLEM = CONTRACT6_6(COEFFS6,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',6,RANK
-        STOP
-        END SELECT
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N = ',NLOOPLINE
-        STOP
-        END SELECT
-      ENDIF
-
-      RES(1)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%C+2.0*LOG(MU_R)
-     $ *RES_GOLEM%%B+2.0*LOG(MU_R)**2*RES_GOLEM%%A)
-      RES(2)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%B+2.0*LOG(MU_R)
-     $ *RES_GOLEM%%A)
-      RES(3)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%A)
-C     WRITE(*,*) 'Loop ID',ID,' =',RES(1),RES(2),RES(3)
-
-C     Finally free golem memory and cache
-      CALL EXITGOLEM95()
-
-      END
-
-      FUNCTION GOLEM_LOOPNUM(Q,MU2)
-      USE PRECISION_GOLEM, ONLY: KI
-      REAL(KI),DIMENSION(0:3),INTENT(IN)::Q
-      REAL(KI),INTENT(IN)::MU2
-      COMPLEX(KI)::GOLEM_LOOPNUM
-
-      COMPLEX*16 QQ(0:3),NUM
-      INTEGER I
-
-      DO I=0,3
-        QQ(I)=CMPLX(Q(I),0.0D0,KIND=16)
-      ENDDO
-
-      CALL LOOPNUM(QQ,NUM)
-      GOLEM_LOOPNUM=NUM
-      RETURN
-      END FUNCTION
-
-      SUBROUTINE INITGOLEM()
-C     
-C     INITIALISATION OF GOLEM
-C     
-C     
-C     MODULE
-C     
-      USE PARAMETRE
-C     
-C     LOCAL VARIABLES 
-C     
-      REAL*8 THRS
-      LOGICAL EXT_NUM_FOR_R1
-C     
-C     GLOBAL VARIABLES 
-C     
-      INCLUDE 'MadLoopParams.inc'
-C     ----------
-C     BEGIN CODE
-C     ----------
-
-C     DEFAULT PARAMETERS FOR GOLEM
-C     -------------------------------  
-C     One can chose here to have either just the rational R1 piece 
-C     or everything but the R2
-      RAT_OR_TOT_PAR = TOT
-
-      END
-
-      SUBROUTINE SETUP_KIN_MATRIX(NLOOPLINE,PGOLEM,M2L)
-C     
-C     MODULE
-C     
-      USE MATRICE_S
-      USE PRECISION_GOLEM, ONLY: KI
-C     
-C     ARGUMENTS
-C     
-      INTEGER NLOOPLINE
-      REAL(KI) PGOLEM(NLOOPLINE,0:3)
-      COMPLEX(KI) M2L(NLOOPLINE)
-C     
-C     LOCAL VARIABLES
-C     
-      INTEGER I,J
-      COMPLEX*16 S_MAT_FROM_MG(NLOOPLINE,NLOOPLINE)
-C     ----------
-C     BEGIN CODE
-C     ----------
-
-      CALL BUILD_KINEMATIC_MATRIX(NLOOPLINE,PGOLEM,M2L,S_MAT_FROM_MG)
-
-      DO I=1,NLOOPLINE
-        DO J=1,NLOOPLINE
-          S_MAT(I,J)=S_MAT_FROM_MG(I,J)
-        ENDDO
-      ENDDO
-
-      END
-
-      FUNCTION COMPARE_COEFS_0(COEFS_A,COEFS_B)
-
-      USE PRECISION_GOLEM, ONLY: KI
-      COMPLEX(KI) COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_0
-
-      NUM = ABS(COEFS_A-COEFS_B)
-      DENOM = ABS(COEFS_A+COEFS_B)
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_0=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_0=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION COMPARE_COEFS_1(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_1)COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_1
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )
-
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_1=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_1=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION COMPARE_COEFS_2(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_2) COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_2
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_2=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_2=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION COMPARE_COEFS_3(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_3) COEFS_A, COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_3
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3)
-     $ )
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_3=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_3=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION COMPARE_COEFS_4(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_4) COEFS_A, COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_4
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
-     $ +SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3)
-     $ )+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_4=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_4=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION COMPARE_COEFS_5(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_5) COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_5
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
-     $ +SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3)
-     $ )+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_5=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_5=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION COMPARE_COEFS_6(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_6) COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_6
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
-     $ +SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3)
-     $ )+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_6=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_6=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-
-      SUBROUTINE FILL_GOLEM_COEFFS_0(ML_COEFS,GOLEM_COEFS)
-      USE PRECISION_GOLEM, ONLY: KI
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      COMPLEX(KI) GOLEM_COEFS
-      GOLEM_COEFS=ML_COEFS(0)
-      END
-
-      SUBROUTINE FILL_GOLEM_COEFFS_1(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_1
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_1) GOLEM_COEFS
-C     Constant coefficient 
-      GOLEM_COEFS%%C0=ML_COEFS(0)
-C     Coefficient q(0)
-      GOLEM_COEFS%%C1(1,1)=-ML_COEFS(1)
-C     Coefficient q(1)
-      GOLEM_COEFS%%C1(2,1)=-ML_COEFS(2)
-C     Coefficient q(2)
-      GOLEM_COEFS%%C1(3,1)=-ML_COEFS(3)
-C     Coefficient q(3)
-      GOLEM_COEFS%%C1(4,1)=-ML_COEFS(4)
-      END
-
-      SUBROUTINE FILL_GOLEM_COEFFS_2(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_2
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_2) GOLEM_COEFS
-C     Constant coefficient 
-      GOLEM_COEFS%%C0=ML_COEFS(0)
-C     Coefficient q(0)
-      GOLEM_COEFS%%C1(1,1)=-ML_COEFS(1)
-C     Coefficient q(0)^2
-      GOLEM_COEFS%%C1(1,2)= ML_COEFS(5)
-C     Coefficient q(1)
-      GOLEM_COEFS%%C1(2,1)=-ML_COEFS(2)
-C     Coefficient q(1)^2
-      GOLEM_COEFS%%C1(2,2)= ML_COEFS(7)
-C     Coefficient q(2)
-      GOLEM_COEFS%%C1(3,1)=-ML_COEFS(3)
-C     Coefficient q(2)^2
-      GOLEM_COEFS%%C1(3,2)= ML_COEFS(10)
-C     Coefficient q(3)
-      GOLEM_COEFS%%C1(4,1)=-ML_COEFS(4)
-C     Coefficient q(3)^2
-      GOLEM_COEFS%%C1(4,2)= ML_COEFS(14)
-C     Coefficient q(0)*q(1)
-      GOLEM_COEFS%%C2(1,1)= ML_COEFS(6)
-C     Coefficient q(0)*q(2)
-      GOLEM_COEFS%%C2(2,1)= ML_COEFS(8)
-C     Coefficient q(0)*q(3)
-      GOLEM_COEFS%%C2(3,1)= ML_COEFS(11)
-C     Coefficient q(1)*q(2)
-      GOLEM_COEFS%%C2(4,1)= ML_COEFS(9)
-C     Coefficient q(1)*q(3)
-      GOLEM_COEFS%%C2(5,1)= ML_COEFS(12)
-C     Coefficient q(2)*q(3)
-      GOLEM_COEFS%%C2(6,1)= ML_COEFS(13)
-      END
-
-      SUBROUTINE FILL_GOLEM_COEFFS_3(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_3
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_3) GOLEM_COEFS
-C     Dummy routine for FILL_GOLEM_COEFS_3
-      STOP 'ERROR: 3 > 2'
-      END
-
-      SUBROUTINE FILL_GOLEM_COEFFS_4(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_4
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_4) GOLEM_COEFS
-C     Dummy routine for FILL_GOLEM_COEFS_4
-      STOP 'ERROR: 4 > 2'
-      END
-
-      SUBROUTINE FILL_GOLEM_COEFFS_5(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_5
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_5) GOLEM_COEFS
-C     Dummy routine for FILL_GOLEM_COEFS_5
-      STOP 'ERROR: 5 > 2'
-      END
-
-      SUBROUTINE FILL_GOLEM_COEFFS_6(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_6
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_6) GOLEM_COEFS
-C     Dummy routine for FILL_GOLEM_COEFS_6
-      STOP 'ERROR: 6 > 2'
-      END
 
=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f	2016-08-07 07:16:07 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%TIR_interface.f	2018-11-09 10:32:12 +0000
@@ -96,17 +96,9 @@
       ENDIF
 
       IF (MLREDUCTIONLIB(I_LIB).EQ.4) THEN
-C       Using Golem95
-C       PDEN is dummy for Golem95 so we just initialize it to zero
-C        here so as to use it for the function SWITCHORDER
-        DO I=0,3
-          DO J=1,NLOOPLINE-1
-            PDEN(I,J)=0.0D0
-          ENDDO
-        ENDDO
-        CALL SWITCH_ORDER(CTMODE,NLOOPLINE,PL,PDEN,M2L)
-        CALL GOLEMLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
-        RETURN
+C       Golem95 not available
+        WRITE(*,*) 'ERROR:: Golem95 is not interfaced.'
+        STOP
       ENDIF
 
 C     INITIALIZE TIR IF NEEDED
 
=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f	2017-08-14 07:04:15 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_dxu_wp%V0_dxu_wp%loop_matrix.f	2018-11-09 10:32:12 +0000
@@ -244,7 +244,7 @@
 C     A FLAG TO DENOTE WHETHER THE CORRESPONDING LOOPLIBS ARE
 C      AVAILABLE OR NOT
       LOGICAL LOOPLIBS_AVAILABLE(NLOOPLIB)
-      DATA LOOPLIBS_AVAILABLE/.TRUE.,.FALSE.,.TRUE.,.TRUE.,.FALSE.
+      DATA LOOPLIBS_AVAILABLE/.TRUE.,.FALSE.,.TRUE.,.FALSE.,.FALSE.
      $ ,.TRUE.,.TRUE./
       COMMON/LOOPLIBS_AV/ LOOPLIBS_AVAILABLE
 C     A FLAG TO DENOTE WHETHER THE CORRESPONDING DIRECTION TESTS
 
=== removed file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%GOLEM_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%GOLEM_interface.f	2016-08-07 07:16:07 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%GOLEM_interface.f	1970-01-01 00:00:00 +0000
@@ -1,748 +0,0 @@
-      SUBROUTINE GOLEMLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
-C     
-C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
-C     By the MadGraph5_aMC@NLO Development Team
-C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
-C     
-C     Interface between MG5 and Golem95.
-C     The Golem95 version should be higher than 1.3.0.
-C     It supports RANK = NLOOPLINE + 1 tensor integrals when 1 <
-C      NLOOPLINE < 6.
-C     
-C     Process: u d~ > w+ QED<=1 WEIGHTED<=2 [ all = QCD ]
-C     Process: c s~ > w+ QED<=1 WEIGHTED<=2 [ all = QCD ]
-C     
-C     
-C     MODULES
-C     
-      USE MATRICE_S
-      USE FORM_FACTOR_TYPE, ONLY: FORM_FACTOR
-      USE PRECISION_GOLEM, ONLY: KI
-      USE TENS_COMB
-      USE TENS_REC
-      USE FORM_FACTOR_1P, ONLY: A10
-      USE FORM_FACTOR_2P, ONLY: A20
-      USE FORM_FACTOR_3P, ONLY: A30
-      USE FORM_FACTOR_4P, ONLY: A40
-      USE FORM_FACTOR_5P, ONLY: A50
-      USE FORM_FACTOR_6P, ONLY: A60
-      IMPLICIT NONE
-C     
-C     CONSTANTS 
-C     
-      INTEGER    NEXTERNAL
-      PARAMETER (NEXTERNAL=3)
-      LOGICAL CHECKPCONSERVATION
-      PARAMETER (CHECKPCONSERVATION=.TRUE.)
-      REAL*8 NORMALIZATION
-      PARAMETER (NORMALIZATION = 1.D0/(16.D0*3.14159265358979323846D0*
-     $ *2))
-      REAL(KI),DIMENSION(0:3),PARAMETER::NULL_VEC = (/0.0_KI,0.0_KI
-     $ ,0.0_KI,0.0_KI/)
-C     GOLEM_RUN_MODE = 1: Use directly MadLoop tensorial coefficients
-C     GOLEM_RUN_MODE = 2: Reconstruct the tensorial coefficeints
-C      directly from 
-C     numerator using golem internal reconstruction routine
-C     GOLEM_RUN_MODE = 3: Cross-checked reconstructed coefficients
-C      against
-C     MadLoop internal ones.
-      INTEGER GOLEM_RUN_MODE
-      PARAMETER (GOLEM_RUN_MODE=1)
-C     The following is the acceptance threshold used for
-C      GOLEM_RUN_MODE = 3
-      REAL*8 COEF_CHECK_THRS
-      DATA COEF_CHECK_THRS/1.0D-13/
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-
-      LOGICAL PASS_COEF_CHECK
-C     
-C     ARGUMENTS 
-C     
-      INTEGER NLOOPLINE, RANK
-      REAL*8 PL(0:3,NLOOPLINE)
-      REAL*8 PCT(0:3,0:NLOOPLINE-1), ABSPCT(0:3)
-      REAL*8 REF_P
-      REAL(KI) PGOLEM(NLOOPLINE,0:3)
-      COMPLEX*16 M2L(NLOOPLINE)
-      COMPLEX(KI) M2LGOLEM(NLOOPLINE)
-      COMPLEX*16 RES(3)
-      LOGICAL STABLE
-C     
-C     LOCAL VARIABLES 
-C     
-      INTEGER I, J, K
-      TYPE(FORM_FACTOR)::RES_GOLEM
-
-      COMPLEX(KI)::COEFFS0,COEFFS0_REC
-      TYPE(COEFF_TYPE_1)::COEFFS1,COEFFS1_REC
-      TYPE(COEFF_TYPE_2)::COEFFS2,COEFFS2_REC
-      TYPE(COEFF_TYPE_3)::COEFFS3,COEFFS3_REC
-      TYPE(COEFF_TYPE_4)::COEFFS4,COEFFS4_REC
-      TYPE(COEFF_TYPE_5)::COEFFS5,COEFFS5_REC
-      TYPE(COEFF_TYPE_6)::COEFFS6,COEFFS6_REC
-
-C     The pinch propagator optimization is not used, so for now it is
-C     always 0.
-      INTEGER PINCH
-C     
-C     EXTERNAL FUNCTIONS
-C     
-      COMPLEX(KI) GOLEM_LOOPNUM
-      EXTERNAL GOLEM_LOOPNUM
-      LOGICAL COMPARE_COEFS_0
-      LOGICAL COMPARE_COEFS_1
-      LOGICAL COMPARE_COEFS_2
-      LOGICAL COMPARE_COEFS_3
-      LOGICAL COMPARE_COEFS_4
-      LOGICAL COMPARE_COEFS_5
-      LOGICAL COMPARE_COEFS_6
-C     
-C     GLOBAL VARIABLES
-C     
-      INCLUDE 'coupl.inc'
-      INTEGER CTMODE
-      REAL*8 LSCALE
-      COMMON/CT/LSCALE,CTMODE
-
-      INTEGER ID,SQSOINDEX,R
-      COMMON/LOOP/ID,SQSOINDEX,R
-
-      LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
-     $ ,COLLIERINIT
-      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT,GOLEMINIT,SAMURAIINIT
-     $ ,NINJAINIT,COLLIERINIT
-
-      INTEGER NLOOPGROUPS
-      PARAMETER (NLOOPGROUPS=1)
-      INTEGER NSQUAREDSO
-      PARAMETER (NSQUAREDSO=1)
-      INCLUDE 'loop_max_coefs.inc'
-
-      COMPLEX*16 LOOPCOEFS(0:LOOPMAXCOEFS-1,NSQUAREDSO,NLOOPGROUPS)
-      COMMON/LCOEFS/LOOPCOEFS
-C     ----------
-C     BEGIN CODE
-C     ----------
-
-C     The CT initialization is also performed here if not done already
-C      because it calls MPINIT of OneLOop which is necessary on some
-C      system
-      IF (CTINIT) THEN
-        CTINIT=.FALSE.
-        CALL INITCT()
-      ENDIF
-
-C     INITIALIZE GOLEM IF NEEDED
-      IF (GOLEMINIT) THEN
-        GOLEMINIT=.FALSE.
-        CALL INITGOLEM()
-      ENDIF
-
-C     No stability test intrisic to Golem95 now
-      STABLE=.TRUE.
-
-C     This initialization must be done for each reduction because we
-C     have not setup anyoptimization using pinched propagators yet.
-      CALL INITGOLEM95(NLOOPLINE)
-      PINCH = 0
-
-C     YOU CAN FIND THE DETAILS ABOUT THE DIFFERENT CTMODE AT THE
-C      BEGINNING OF THE FILE CTS_CUTS.F90 IN THE CUTTOOLS DISTRIBUTION
-
-C     CONVERT THE MASSES TO BE COMPLEX
-      DO I=1,NLOOPLINE
-        M2LGOLEM(I)=M2L(I)
-      ENDDO
-
-C     CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO CT CONVENTIONS
-      DO I=0,3
-        ABSPCT(I) = 0.D0
-        DO J=0,(NLOOPLINE-1)
-          PCT(I,J)=0.D0
-        ENDDO
-      ENDDO
-      DO I=0,3
-        DO J=1,NLOOPLINE
-          PCT(I,0)=PCT(I,0)+PL(I,J)
-          ABSPCT(I)=ABSPCT(I)+ABS(PL(I,J))
-        ENDDO
-      ENDDO
-      REF_P = MAX(ABSPCT(0), ABSPCT(1),ABSPCT(2),ABSPCT(3))
-      DO I=0,3
-        ABSPCT(I) = MAX(REF_P*1E-6, ABSPCT(I))
-      ENDDO
-      IF (CHECKPCONSERVATION.AND.REF_P.GT.1D-8) THEN
-        IF ((PCT(0,0)/ABSPCT(0)).GT.1.D-6) THEN
-          WRITE(*,*) 'energy is not conserved ',PCT(0,0)
-          STOP 'energy is not conserved'
-        ELSEIF ((PCT(1,0)/ABSPCT(1)).GT.1.D-6) THEN
-          WRITE(*,*) 'px is not conserved ',PCT(1,0)
-          STOP 'px is not conserved'
-        ELSEIF ((PCT(2,0)/ABSPCT(2)).GT.1.D-6) THEN
-          WRITE(*,*) 'py is not conserved ',PCT(2,0)
-          STOP 'py is not conserved'
-        ELSEIF ((PCT(3,0)/ABSPCT(3)).GT.1.D-6) THEN
-          WRITE(*,*) 'pz is not conserved ',PCT(3,0)
-          STOP 'pz is not conserved'
-        ENDIF
-      ENDIF
-      DO I=0,3
-        DO J=1,(NLOOPLINE-1)
-          DO K=1,J
-            PCT(I,J)=PCT(I,J)+PL(I,K)
-          ENDDO
-        ENDDO
-      ENDDO
-
-C     Now convert the loop momenta to Golem95 conventions
-      DO I=0,3
-        PGOLEM(1,I)=0.0E0_KI
-        DO J=2,NLOOPLINE
-          PGOLEM(J,I)=PCT(I,J-1)
-        ENDDO
-      ENDDO
-
-C     Fill in the kinematic s-matrix while taking care of on-shell
-C      limits.
-      CALL SETUP_KIN_MATRIX(NLOOPLINE,PGOLEM,M2LGOLEM)
-C     Construct the golem internal matrices derived from the kinetic
-C      one.
-      CALL PREPARESMATRIX()
-
-C     Fill in the golem coefficents and compute the loop
-      IF(GOLEM_RUN_MODE.EQ.2)THEN
-        RES_GOLEM = EVALUATE_B(GOLEM_LOOPNUM,PGOLEM,0,RANK)
-      ELSE
-        PASS_COEF_CHECK=.TRUE.
-        SELECT CASE(RANK)
-        CASE(0)
-        CALL FILL_GOLEM_COEFFS_0(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS0)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          COEFFS0_REC = GOLEM_LOOPNUM(NULL_VEC,0.0_KI)
-          PASS_COEF_CHECK=COMPARE_COEFS_0(COEFFS0,COEFFS0_REC)
-        ENDIF
-        CASE(1)
-        CALL FILL_GOLEM_COEFFS_1(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS1)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT1(GOLEM_LOOPNUM,COEFFS1_REC)
-          PASS_COEF_CHECK=COMPARE_COEFS_1(COEFFS1,COEFFS1_REC)
-        ENDIF
-        CASE(2)
-        CALL FILL_GOLEM_COEFFS_2(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS2)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT2(GOLEM_LOOPNUM,COEFFS2_REC)
-          PASS_COEF_CHECK=COMPARE_COEFS_2(COEFFS2,COEFFS2_REC)
-        ENDIF
-        CASE(3)
-        CALL FILL_GOLEM_COEFFS_3(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS3)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT3(GOLEM_LOOPNUM,COEFFS3_REC)
-          PASS_COEF_CHECK=COMPARE_COEFS_3(COEFFS3,COEFFS3_REC)
-        ENDIF
-        CASE(4)
-        CALL FILL_GOLEM_COEFFS_4(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS4)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT4(GOLEM_LOOPNUM,COEFFS4_REC)
-          PASS_COEF_CHECK=COMPARE_COEFS_4(COEFFS4,COEFFS4_REC)
-        ENDIF
-        CASE(5)
-        CALL FILL_GOLEM_COEFFS_5(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS5)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT5(GOLEM_LOOPNUM,COEFFS5_REC)
-          PASS_COEF_CHECK=COMPARE_COEFS_5(COEFFS5,COEFFS5_REC)
-        ENDIF
-        CASE(6)
-        CALL FILL_GOLEM_COEFFS_6(LOOPCOEFS(0,SQSOINDEX,ID),COEFFS6)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT6(GOLEM_LOOPNUM,COEFFS6_REC)
-          PASS_COEF_CHECK=COMPARE_COEFS_6(COEFFS6,COEFFS6_REC)
-        ENDIF
-        CASE DEFAULT
-        WRITE(*,*)'Not yet implemented in Golem95 for rank= ',RANK
-        STOP
-        END SELECT
-
-        IF(.NOT.PASS_COEF_CHECK)THEN
-          WRITE(*,*)'Coefs mismatch for ID ',ID,' and rank ',RANK
-          WRITE(*,*)'Coefs form MadLoop5:'
-          SELECT CASE(RANK)
-          CASE(0)
-          WRITE(*,*)'Constant coef = ',COEFFS0
-          CASE(1)
-          CALL PRINT_COEFFS(COEFFS1)
-          CASE(2)
-          CALL PRINT_COEFFS(COEFFS2)
-          CASE(3)
-          CALL PRINT_COEFFS(COEFFS3)
-          CASE(4)
-          CALL PRINT_COEFFS(COEFFS4)
-          CASE(5)
-          CALL PRINT_COEFFS(COEFFS5)
-          CASE(6)
-          CALL PRINT_COEFFS(COEFFS6)
-          END SELECT
-          WRITE(*,*)'Coefs reconstructed by Golem95:'
-          SELECT CASE(RANK)
-          CASE(0)
-          WRITE(*,*)'Constant coef = ',COEFFS0_REC
-          CASE(1)
-          CALL PRINT_COEFFS(COEFFS1_REC)
-          CASE(2)
-          CALL PRINT_COEFFS(COEFFS2_REC)
-          CASE(3)
-          CALL PRINT_COEFFS(COEFFS3_REC)
-          CASE(4)
-          CALL PRINT_COEFFS(COEFFS4_REC)
-          CASE(5)
-          CALL PRINT_COEFFS(COEFFS5_REC)
-          CASE(6)
-          CALL PRINT_COEFFS(COEFFS6_REC)
-          END SELECT
-          STOP
-        ENDIF
-
-        SELECT CASE(NLOOPLINE)
-        CASE(1)
-        WRITE(*,*)'Golem95 cannot handle with tadpole yet'
-        STOP
-        CASE(2)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A20(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT2_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT2_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT2_3(COEFFS3,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',2,RANK
-        STOP
-        END SELECT
-        CASE(3)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A30(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT3_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT3_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT3_3(COEFFS3,PGOLEM,PINCH)
-        CASE(4)
-        RES_GOLEM = CONTRACT3_4(COEFFS4,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',3,RANK
-        STOP
-        END SELECT
-        CASE(4)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A40(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT4_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT4_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT4_3(COEFFS3,PGOLEM,PINCH)
-        CASE(4)
-        RES_GOLEM = CONTRACT4_4(COEFFS4,PGOLEM,PINCH)
-        CASE(5)
-        RES_GOLEM = CONTRACT4_5(COEFFS5,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',4,RANK
-        STOP
-        END SELECT
-        CASE(5)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A50(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT5_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT5_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT5_3(COEFFS3,PGOLEM,PINCH)
-        CASE(4)
-        RES_GOLEM = CONTRACT5_4(COEFFS4,PGOLEM,PINCH)
-        CASE(5)
-        RES_GOLEM = CONTRACT5_5(COEFFS5,PGOLEM,PINCH)
-        CASE(6)
-        RES_GOLEM = CONTRACT5_6(COEFFS6,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',5,RANK
-        STOP
-        END SELECT
-        CASE(6)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A60(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT6_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT6_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT6_3(COEFFS3,PGOLEM,PINCH)
-        CASE(4)
-        RES_GOLEM = CONTRACT6_4(COEFFS4,PGOLEM,PINCH)
-        CASE(5)
-        RES_GOLEM = CONTRACT6_5(COEFFS5,PGOLEM,PINCH)
-        CASE(6)
-        RES_GOLEM = CONTRACT6_6(COEFFS6,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',6,RANK
-        STOP
-        END SELECT
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N = ',NLOOPLINE
-        STOP
-        END SELECT
-      ENDIF
-
-      RES(1)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%C+2.0*LOG(MU_R)
-     $ *RES_GOLEM%%B+2.0*LOG(MU_R)**2*RES_GOLEM%%A)
-      RES(2)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%B+2.0*LOG(MU_R)
-     $ *RES_GOLEM%%A)
-      RES(3)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%A)
-C     WRITE(*,*) 'Loop ID',ID,' =',RES(1),RES(2),RES(3)
-
-C     Finally free golem memory and cache
-      CALL EXITGOLEM95()
-
-      END
-
-      FUNCTION GOLEM_LOOPNUM(Q,MU2)
-      USE PRECISION_GOLEM, ONLY: KI
-      REAL(KI),DIMENSION(0:3),INTENT(IN)::Q
-      REAL(KI),INTENT(IN)::MU2
-      COMPLEX(KI)::GOLEM_LOOPNUM
-
-      COMPLEX*16 QQ(0:3),NUM
-      INTEGER I
-
-      DO I=0,3
-        QQ(I)=CMPLX(Q(I),0.0D0,KIND=16)
-      ENDDO
-
-      CALL LOOPNUM(QQ,NUM)
-      GOLEM_LOOPNUM=NUM
-      RETURN
-      END FUNCTION
-
-      SUBROUTINE INITGOLEM()
-C     
-C     INITIALISATION OF GOLEM
-C     
-C     
-C     MODULE
-C     
-      USE PARAMETRE
-C     
-C     LOCAL VARIABLES 
-C     
-      REAL*8 THRS
-      LOGICAL EXT_NUM_FOR_R1
-C     
-C     GLOBAL VARIABLES 
-C     
-      INCLUDE 'MadLoopParams.inc'
-C     ----------
-C     BEGIN CODE
-C     ----------
-
-C     DEFAULT PARAMETERS FOR GOLEM
-C     -------------------------------  
-C     One can chose here to have either just the rational R1 piece 
-C     or everything but the R2
-      RAT_OR_TOT_PAR = TOT
-
-      END
-
-      SUBROUTINE SETUP_KIN_MATRIX(NLOOPLINE,PGOLEM,M2L)
-C     
-C     MODULE
-C     
-      USE MATRICE_S
-      USE PRECISION_GOLEM, ONLY: KI
-C     
-C     ARGUMENTS
-C     
-      INTEGER NLOOPLINE
-      REAL(KI) PGOLEM(NLOOPLINE,0:3)
-      COMPLEX(KI) M2L(NLOOPLINE)
-C     
-C     LOCAL VARIABLES
-C     
-      INTEGER I,J
-      COMPLEX*16 S_MAT_FROM_MG(NLOOPLINE,NLOOPLINE)
-C     ----------
-C     BEGIN CODE
-C     ----------
-
-      CALL BUILD_KINEMATIC_MATRIX(NLOOPLINE,PGOLEM,M2L,S_MAT_FROM_MG)
-
-      DO I=1,NLOOPLINE
-        DO J=1,NLOOPLINE
-          S_MAT(I,J)=S_MAT_FROM_MG(I,J)
-        ENDDO
-      ENDDO
-
-      END
-
-      FUNCTION COMPARE_COEFS_0(COEFS_A,COEFS_B)
-
-      USE PRECISION_GOLEM, ONLY: KI
-      COMPLEX(KI) COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_0
-
-      NUM = ABS(COEFS_A-COEFS_B)
-      DENOM = ABS(COEFS_A+COEFS_B)
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_0=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_0=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION COMPARE_COEFS_1(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_1)COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_1
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )
-
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_1=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_1=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION COMPARE_COEFS_2(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_2) COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_2
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_2=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_2=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION COMPARE_COEFS_3(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_3) COEFS_A, COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_3
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3)
-     $ )
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_3=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_3=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION COMPARE_COEFS_4(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_4) COEFS_A, COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_4
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
-     $ +SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3)
-     $ )+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_4=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_4=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION COMPARE_COEFS_5(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_5) COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_5
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
-     $ +SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3)
-     $ )+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_5=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_5=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION COMPARE_COEFS_6(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_6) COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL COMPARE_COEFS_6
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
-     $ +SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3)
-     $ )+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
-      IF(DENOM.GT.0D0)THEN
-        COMPARE_COEFS_6=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        COMPARE_COEFS_6=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-
-      SUBROUTINE FILL_GOLEM_COEFFS_0(ML_COEFS,GOLEM_COEFS)
-      USE PRECISION_GOLEM, ONLY: KI
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      COMPLEX(KI) GOLEM_COEFS
-      GOLEM_COEFS=ML_COEFS(0)
-      END
-
-      SUBROUTINE FILL_GOLEM_COEFFS_1(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_1
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_1) GOLEM_COEFS
-C     Constant coefficient 
-      GOLEM_COEFS%%C0=ML_COEFS(0)
-C     Coefficient q(0)
-      GOLEM_COEFS%%C1(1,1)=-ML_COEFS(1)
-C     Coefficient q(1)
-      GOLEM_COEFS%%C1(2,1)=-ML_COEFS(2)
-C     Coefficient q(2)
-      GOLEM_COEFS%%C1(3,1)=-ML_COEFS(3)
-C     Coefficient q(3)
-      GOLEM_COEFS%%C1(4,1)=-ML_COEFS(4)
-      END
-
-      SUBROUTINE FILL_GOLEM_COEFFS_2(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_2
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_2) GOLEM_COEFS
-C     Constant coefficient 
-      GOLEM_COEFS%%C0=ML_COEFS(0)
-C     Coefficient q(0)
-      GOLEM_COEFS%%C1(1,1)=-ML_COEFS(1)
-C     Coefficient q(0)^2
-      GOLEM_COEFS%%C1(1,2)= ML_COEFS(5)
-C     Coefficient q(1)
-      GOLEM_COEFS%%C1(2,1)=-ML_COEFS(2)
-C     Coefficient q(1)^2
-      GOLEM_COEFS%%C1(2,2)= ML_COEFS(7)
-C     Coefficient q(2)
-      GOLEM_COEFS%%C1(3,1)=-ML_COEFS(3)
-C     Coefficient q(2)^2
-      GOLEM_COEFS%%C1(3,2)= ML_COEFS(10)
-C     Coefficient q(3)
-      GOLEM_COEFS%%C1(4,1)=-ML_COEFS(4)
-C     Coefficient q(3)^2
-      GOLEM_COEFS%%C1(4,2)= ML_COEFS(14)
-C     Coefficient q(0)*q(1)
-      GOLEM_COEFS%%C2(1,1)= ML_COEFS(6)
-C     Coefficient q(0)*q(2)
-      GOLEM_COEFS%%C2(2,1)= ML_COEFS(8)
-C     Coefficient q(0)*q(3)
-      GOLEM_COEFS%%C2(3,1)= ML_COEFS(11)
-C     Coefficient q(1)*q(2)
-      GOLEM_COEFS%%C2(4,1)= ML_COEFS(9)
-C     Coefficient q(1)*q(3)
-      GOLEM_COEFS%%C2(5,1)= ML_COEFS(12)
-C     Coefficient q(2)*q(3)
-      GOLEM_COEFS%%C2(6,1)= ML_COEFS(13)
-      END
-
-      SUBROUTINE FILL_GOLEM_COEFFS_3(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_3
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_3) GOLEM_COEFS
-C     Dummy routine for FILL_GOLEM_COEFS_3
-      STOP 'ERROR: 3 > 2'
-      END
-
-      SUBROUTINE FILL_GOLEM_COEFFS_4(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_4
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_4) GOLEM_COEFS
-C     Dummy routine for FILL_GOLEM_COEFS_4
-      STOP 'ERROR: 4 > 2'
-      END
-
-      SUBROUTINE FILL_GOLEM_COEFFS_5(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_5
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_5) GOLEM_COEFS
-C     Dummy routine for FILL_GOLEM_COEFS_5
-      STOP 'ERROR: 5 > 2'
-      END
-
-      SUBROUTINE FILL_GOLEM_COEFFS_6(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_6
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_6) GOLEM_COEFS
-C     Dummy routine for FILL_GOLEM_COEFS_6
-      STOP 'ERROR: 6 > 2'
-      END
 
=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%TIR_interface.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%TIR_interface.f	2016-08-07 07:16:07 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%TIR_interface.f	2018-11-09 10:32:12 +0000
@@ -96,17 +96,9 @@
       ENDIF
 
       IF (MLREDUCTIONLIB(I_LIB).EQ.4) THEN
-C       Using Golem95
-C       PDEN is dummy for Golem95 so we just initialize it to zero
-C        here so as to use it for the function SWITCHORDER
-        DO I=0,3
-          DO J=1,NLOOPLINE-1
-            PDEN(I,J)=0.0D0
-          ENDDO
-        ENDDO
-        CALL SWITCH_ORDER(CTMODE,NLOOPLINE,PL,PDEN,M2L)
-        CALL GOLEMLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
-        RETURN
+C       Golem95 not available
+        WRITE(*,*) 'ERROR:: Golem95 is not interfaced.'
+        STOP
       ENDIF
 
 C     INITIALIZE TIR IF NEEDED
 
=== modified file 'tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%loop_matrix.f'
--- tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%loop_matrix.f	2017-08-14 07:04:15 +0000
+++ tests/input_files/IOTestsComparison/IOExportFKSTest/test_ppw_fksall/%SubProcesses%P0_udx_wp%V0_udx_wp%loop_matrix.f	2018-11-09 10:32:12 +0000
@@ -244,7 +244,7 @@
 C     A FLAG TO DENOTE WHETHER THE CORRESPONDING LOOPLIBS ARE
 C      AVAILABLE OR NOT
       LOGICAL LOOPLIBS_AVAILABLE(NLOOPLIB)
-      DATA LOOPLIBS_AVAILABLE/.TRUE.,.FALSE.,.TRUE.,.TRUE.,.FALSE.
+      DATA LOOPLIBS_AVAILABLE/.TRUE.,.FALSE.,.TRUE.,.FALSE.,.FALSE.
      $ ,.TRUE.,.TRUE./
       COMMON/LOOPLIBS_AV/ LOOPLIBS_AVAILABLE
 C     A FLAG TO DENOTE WHETHER THE CORRESPONDING DIRECTION TESTS
 
=== 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	2017-11-24 12:17:04 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_group/matrix1.f	2018-11-09 10:32:12 +0000
@@ -73,6 +73,9 @@
       REAL*8 POL(2)
       COMMON/TO_POLARIZATION/ POL
 
+      DOUBLE PRECISION SMALL_WIDTH_TREATMENT
+      COMMON/NARROW_WIDTH/SMALL_WIDTH_TREATMENT
+
       INTEGER          ISUM_HEL
       LOGICAL                    MULTI_CHANNEL
       COMMON/TO_MATRIX/ISUM_HEL, MULTI_CHANNEL
@@ -110,6 +113,7 @@
 C     ----------
 C     BEGIN CODE
 C     ----------
+
       NTRY(IMIRROR)=NTRY(IMIRROR)+1
       THIS_NTRY(IMIRROR) = THIS_NTRY(IMIRROR)+1
       DO I=1,NEXTERNAL
@@ -290,6 +294,15 @@
 C     Needed for v4 models
       COMPLEX*16 DUM0,DUM1
       DATA DUM0, DUM1/(0D0, 0D0), (1D0, 0D0)/
+
+      DOUBLE PRECISION FK_ZERO
+      DOUBLE PRECISION FK_WZ
+      SAVE FK_ZERO
+      SAVE FK_WZ
+
+      LOGICAL FIRST
+      DATA FIRST /.TRUE./
+      SAVE FIRST
 C     
 C     FUNCTION
 C     
@@ -300,6 +313,9 @@
       DOUBLE PRECISION AMP2(MAXAMPS), JAMP2(0:MAXFLOW)
       COMMON/TO_AMPS/  AMP2,       JAMP2
       INCLUDE 'coupl.inc'
+
+      DOUBLE PRECISION SMALL_WIDTH_TREATMENT
+      COMMON/NARROW_WIDTH/SMALL_WIDTH_TREATMENT
 C     
 C     COLOR DATA
 C     
@@ -312,26 +328,33 @@
 C     ----------
 C     BEGIN CODE
 C     ----------
+      IF (FIRST) THEN
+        FIRST=.FALSE.
+        FK_WZ = SIGN(MAX(ABS(WZ), ABS(MZ*SMALL_WIDTH_TREATMENT)), WZ)
+        FK_ZERO = 0D0
+      ENDIF
+
+
       CALL IXXXXX(P(0,1),ZERO,NHEL(1),+1*IC(1),W(1,1))
       CALL OXXXXX(P(0,2),ZERO,NHEL(2),-1*IC(2),W(1,2))
       CALL OXXXXX(P(0,3),ZERO,NHEL(3),+1*IC(3),W(1,3))
       CALL IXXXXX(P(0,4),ZERO,NHEL(4),-1*IC(4),W(1,4))
-      CALL FFV1_3(W(1,1),W(1,2),GQQ,ZERO,ZERO,W(1,5))
+      CALL FFV1_3(W(1,1),W(1,2),GQQ,ZERO, FK_ZERO,W(1,5))
 C     Amplitude(s) for diagram number 1
       CALL FFV1_0(W(1,4),W(1,3),W(1,5),GQQ,AMP(1))
-      CALL FFV1_3(W(1,1),W(1,2),GQED,ZERO,ZERO,W(1,5))
+      CALL FFV1_3(W(1,1),W(1,2),GQED,ZERO, FK_ZERO,W(1,5))
 C     Amplitude(s) for diagram number 2
       CALL FFV1_0(W(1,4),W(1,3),W(1,5),GQED,AMP(2))
-      CALL FFV1_2_3(W(1,1),W(1,2),GUZ1,GUZ2,MZ,WZ,W(1,5))
+      CALL FFV1_2_3(W(1,1),W(1,2),GUZ1,GUZ2,MZ, FK_WZ,W(1,5))
 C     Amplitude(s) for diagram number 3
       CALL FFV1_2_0(W(1,4),W(1,3),W(1,5),GUZ1,GUZ2,AMP(3))
-      CALL FFV1_3(W(1,1),W(1,3),GQQ,ZERO,ZERO,W(1,5))
+      CALL FFV1_3(W(1,1),W(1,3),GQQ,ZERO, FK_ZERO,W(1,5))
 C     Amplitude(s) for diagram number 4
       CALL FFV1_0(W(1,4),W(1,2),W(1,5),GQQ,AMP(4))
-      CALL FFV1_3(W(1,1),W(1,3),GQED,ZERO,ZERO,W(1,5))
+      CALL FFV1_3(W(1,1),W(1,3),GQED,ZERO, FK_ZERO,W(1,5))
 C     Amplitude(s) for diagram number 5
       CALL FFV1_0(W(1,4),W(1,2),W(1,5),GQED,AMP(5))
-      CALL FFV1_2_3(W(1,1),W(1,3),GUZ1,GUZ2,MZ,WZ,W(1,5))
+      CALL FFV1_2_3(W(1,1),W(1,3),GUZ1,GUZ2,MZ, FK_WZ,W(1,5))
 C     Amplitude(s) for diagram number 6
       CALL FFV1_2_0(W(1,4),W(1,2),W(1,5),GUZ1,GUZ2,AMP(6))
 C     JAMPs contributing to orders ALL_ORDERS=1
 
=== 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	2017-11-26 00:21:49 +0000
+++ tests/input_files/IOTestsComparison/IOExportV4IOTest/export_matrix_element_v4_madevent_nogroup/matrix.f	2018-11-09 10:32:12 +0000
@@ -51,6 +51,7 @@
       INTEGER IDUM, NGOOD, IGOOD(NCOMB), JHEL, J, JJ
       REAL     XRAN1
       EXTERNAL XRAN1
+
 C     
 C     GLOBAL VARIABLES
 C     
@@ -71,6 +72,7 @@
       DATA IDUM /-1/
       DATA XTRY, XREJ, NGOOD /0,0,0/
       SAVE YFRAC, IGOOD, JHEL
+
       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/
@@ -261,6 +263,13 @@
 C     Needed for v4 models
       COMPLEX*16 DUM0,DUM1
       DATA DUM0, DUM1/(0D0, 0D0), (1D0, 0D0)/
+
+      DOUBLE PRECISION FK_ZERO
+      SAVE FK_ZERO
+
+      LOGICAL FIRST
+      DATA FIRST /.TRUE./
+      SAVE FIRST
 C     
 C     FUNCTION
 C     
@@ -271,6 +280,9 @@
       DOUBLE PRECISION AMP2(MAXAMPS), JAMP2(0:MAXFLOW)
       COMMON/TO_AMPS/  AMP2,       JAMP2
       INCLUDE 'coupl.inc'
+
+      DOUBLE PRECISION SMALL_WIDTH_TREATMENT
+      COMMON/NARROW_WIDTH/SMALL_WIDTH_TREATMENT
 C     
 C     COLOR DATA
 C     
@@ -283,17 +295,22 @@
 C     ----------
 C     BEGIN CODE
 C     ----------
+      IF (FIRST) THEN
+        FIRST=.FALSE.
+        FK_ZERO = 0D0
+      ENDIF
+
       CALL IXXXXX(P(0,1),ZERO,NHEL(1),+1*IC(1),W(1,1))
       CALL OXXXXX(P(0,2),ZERO,NHEL(2),-1*IC(2),W(1,2))
       CALL VXXXXX(P(0,3),ZERO,NHEL(3),+1*IC(3),W(1,3))
       CALL VXXXXX(P(0,4),ZERO,NHEL(4),+1*IC(4),W(1,4))
-      CALL FFV1_3(W(1,1),W(1,2),GQQ,ZERO,ZERO,W(1,5))
+      CALL FFV1_3(W(1,1),W(1,2),GQQ,ZERO, FK_ZERO,W(1,5))
 C     Amplitude(s) for diagram number 1
       CALL VVV1_0(W(1,3),W(1,4),W(1,5),G,AMP(1))
-      CALL FFV1_2(W(1,1),W(1,3),GQQ,ZERO,ZERO,W(1,5))
+      CALL FFV1_2(W(1,1),W(1,3),GQQ,ZERO, FK_ZERO,W(1,5))
 C     Amplitude(s) for diagram number 2
       CALL FFV1_0(W(1,5),W(1,2),W(1,4),GQQ,AMP(2))
-      CALL FFV1_2(W(1,1),W(1,4),GQQ,ZERO,ZERO,W(1,5))
+      CALL FFV1_2(W(1,1),W(1,4),GQQ,ZERO, FK_ZERO,W(1,5))
 C     Amplitude(s) for diagram number 3
       CALL FFV1_0(W(1,5),W(1,2),W(1,3),GQQ,AMP(3))
 C     JAMPs contributing to orders ALL_ORDERS=1
 
=== removed file 'tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%GOLEM_interface.f'
--- tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%GOLEM_interface.f	2016-05-26 03:24:09 +0000
+++ tests/input_files/IOTestsComparison/MadLoop_output_from_the_interface/TIR_output/%ggttx_IOTest%SubProcesses%P0_gg_ttx%GOLEM_interface.f	1970-01-01 00:00:00 +0000
@@ -1,823 +0,0 @@
-      SUBROUTINE ML5_0_GOLEMLOOP(NLOOPLINE,PL,M2L,RANK,RES,STABLE)
-C     
-C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
-C     By the MadGraph5_aMC@NLO Development Team
-C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
-C     
-C     Interface between MG5 and Golem95.
-C     The Golem95 version should be higher than 1.3.0.
-C     It supports RANK = NLOOPLINE + 1 tensor integrals when 1 <
-C      NLOOPLINE < 6.
-C     
-C     Process: g g > t t~ [ virt = QCD ]
-C     
-C     
-C     MODULES
-C     
-      USE MATRICE_S
-      USE FORM_FACTOR_TYPE, ONLY: FORM_FACTOR
-      USE PRECISION_GOLEM, ONLY: KI
-      USE TENS_COMB
-      USE TENS_REC
-      USE FORM_FACTOR_1P, ONLY: A10
-      USE FORM_FACTOR_2P, ONLY: A20
-      USE FORM_FACTOR_3P, ONLY: A30
-      USE FORM_FACTOR_4P, ONLY: A40
-      USE FORM_FACTOR_5P, ONLY: A50
-      USE FORM_FACTOR_6P, ONLY: A60
-      IMPLICIT NONE
-C     
-C     CONSTANTS 
-C     
-      INTEGER    NEXTERNAL
-      PARAMETER (NEXTERNAL=4)
-      LOGICAL CHECKPCONSERVATION
-      PARAMETER (CHECKPCONSERVATION=.TRUE.)
-      REAL*8 NORMALIZATION
-      PARAMETER (NORMALIZATION = 1.D0/(16.D0*3.14159265358979323846D0*
-     $ *2))
-      REAL(KI),DIMENSION(0:3),PARAMETER::NULL_VEC = (/0.0_KI,0.0_KI
-     $ ,0.0_KI,0.0_KI/)
-C     GOLEM_RUN_MODE = 1: Use directly MadLoop tensorial coefficients
-C     GOLEM_RUN_MODE = 2: Reconstruct the tensorial coefficeints
-C      directly from 
-C     numerator using golem internal reconstruction routine
-C     GOLEM_RUN_MODE = 3: Cross-checked reconstructed coefficients
-C      against
-C     MadLoop internal ones.
-      INTEGER GOLEM_RUN_MODE
-      PARAMETER (GOLEM_RUN_MODE=1)
-C     The following is the acceptance threshold used for
-C      GOLEM_RUN_MODE = 3
-      REAL*8 COEF_CHECK_THRS
-      DATA COEF_CHECK_THRS/1.0D-13/
-      COMMON/ML5_0_COEF_CHECK_THRS/COEF_CHECK_THRS
-
-      LOGICAL PASS_COEF_CHECK
-C     
-C     ARGUMENTS 
-C     
-      INTEGER NLOOPLINE, RANK
-      REAL*8 PL(0:3,NLOOPLINE)
-      REAL*8 PCT(0:3,0:NLOOPLINE-1), ABSPCT(0:3)
-      REAL*8 REF_P
-      REAL(KI) PGOLEM(NLOOPLINE,0:3)
-      COMPLEX*16 M2L(NLOOPLINE)
-      COMPLEX(KI) M2LGOLEM(NLOOPLINE)
-      COMPLEX*16 RES(3)
-      LOGICAL STABLE
-C     
-C     LOCAL VARIABLES 
-C     
-      INTEGER I, J, K
-      TYPE(FORM_FACTOR)::RES_GOLEM
-
-      COMPLEX(KI)::COEFFS0,COEFFS0_REC
-      TYPE(COEFF_TYPE_1)::COEFFS1,COEFFS1_REC
-      TYPE(COEFF_TYPE_2)::COEFFS2,COEFFS2_REC
-      TYPE(COEFF_TYPE_3)::COEFFS3,COEFFS3_REC
-      TYPE(COEFF_TYPE_4)::COEFFS4,COEFFS4_REC
-      TYPE(COEFF_TYPE_5)::COEFFS5,COEFFS5_REC
-      TYPE(COEFF_TYPE_6)::COEFFS6,COEFFS6_REC
-
-C     The pinch propagator optimization is not used, so for now it is
-C     always 0.
-      INTEGER PINCH
-C     
-C     EXTERNAL FUNCTIONS
-C     
-      COMPLEX(KI) ML5_0_GOLEM_LOOPNUM
-      EXTERNAL ML5_0_GOLEM_LOOPNUM
-      LOGICAL ML5_0_COMPARE_COEFS_0
-      LOGICAL ML5_0_COMPARE_COEFS_1
-      LOGICAL ML5_0_COMPARE_COEFS_2
-      LOGICAL ML5_0_COMPARE_COEFS_3
-      LOGICAL ML5_0_COMPARE_COEFS_4
-      LOGICAL ML5_0_COMPARE_COEFS_5
-      LOGICAL ML5_0_COMPARE_COEFS_6
-C     
-C     GLOBAL VARIABLES
-C     
-      INCLUDE 'coupl.inc'
-      INTEGER CTMODE
-      REAL*8 LSCALE
-      COMMON/ML5_0_CT/LSCALE,CTMODE
-
-      INTEGER ID,SQSOINDEX,R
-      COMMON/ML5_0_LOOP/ID,SQSOINDEX,R
-
-      LOGICAL CTINIT, TIRINIT, GOLEMINIT, SAMURAIINIT, NINJAINIT
-     $ ,COLLIERINIT
-      COMMON/REDUCTIONCODEINIT/CTINIT, TIRINIT,GOLEMINIT,SAMURAIINIT
-     $ ,NINJAINIT,COLLIERINIT
-
-      INTEGER NLOOPGROUPS
-      PARAMETER (NLOOPGROUPS=26)
-      INTEGER NSQUAREDSO
-      PARAMETER (NSQUAREDSO=1)
-      INCLUDE 'loop_max_coefs.inc'
-
-      COMPLEX*16 LOOPCOEFS(0:LOOPMAXCOEFS-1,NSQUAREDSO,NLOOPGROUPS)
-      COMMON/ML5_0_LCOEFS/LOOPCOEFS
-C     ----------
-C     BEGIN CODE
-C     ----------
-
-C     The CT initialization is also performed here if not done already
-C      because it calls MPINIT of OneLOop which is necessary on some
-C      system
-      IF (CTINIT) THEN
-        CTINIT=.FALSE.
-        CALL ML5_0_INITCT()
-      ENDIF
-
-C     INITIALIZE GOLEM IF NEEDED
-      IF (GOLEMINIT) THEN
-        GOLEMINIT=.FALSE.
-        CALL ML5_0_INITGOLEM()
-      ENDIF
-
-C     No stability test intrisic to Golem95 now
-      STABLE=.TRUE.
-
-C     This initialization must be done for each reduction because we
-C     have not setup anyoptimization using pinched propagators yet.
-      CALL INITGOLEM95(NLOOPLINE)
-      PINCH = 0
-
-C     YOU CAN FIND THE DETAILS ABOUT THE DIFFERENT CTMODE AT THE
-C      BEGINNING OF THE FILE CTS_CUTS.F90 IN THE CUTTOOLS DISTRIBUTION
-
-C     CONVERT THE MASSES TO BE COMPLEX
-      DO I=1,NLOOPLINE
-        M2LGOLEM(I)=M2L(I)
-      ENDDO
-
-C     CONVERT THE MOMENTA FLOWING IN THE LOOP LINES TO CT CONVENTIONS
-      DO I=0,3
-        ABSPCT(I) = 0.D0
-        DO J=0,(NLOOPLINE-1)
-          PCT(I,J)=0.D0
-        ENDDO
-      ENDDO
-      DO I=0,3
-        DO J=1,NLOOPLINE
-          PCT(I,0)=PCT(I,0)+PL(I,J)
-          ABSPCT(I)=ABSPCT(I)+ABS(PL(I,J))
-        ENDDO
-      ENDDO
-      REF_P = MAX(ABSPCT(0), ABSPCT(1),ABSPCT(2),ABSPCT(3))
-      DO I=0,3
-        ABSPCT(I) = MAX(REF_P*1E-6, ABSPCT(I))
-      ENDDO
-      IF (CHECKPCONSERVATION.AND.REF_P.GT.1D-8) THEN
-        IF ((PCT(0,0)/ABSPCT(0)).GT.1.D-6) THEN
-          WRITE(*,*) 'energy is not conserved ',PCT(0,0)
-          STOP 'energy is not conserved'
-        ELSEIF ((PCT(1,0)/ABSPCT(1)).GT.1.D-6) THEN
-          WRITE(*,*) 'px is not conserved ',PCT(1,0)
-          STOP 'px is not conserved'
-        ELSEIF ((PCT(2,0)/ABSPCT(2)).GT.1.D-6) THEN
-          WRITE(*,*) 'py is not conserved ',PCT(2,0)
-          STOP 'py is not conserved'
-        ELSEIF ((PCT(3,0)/ABSPCT(3)).GT.1.D-6) THEN
-          WRITE(*,*) 'pz is not conserved ',PCT(3,0)
-          STOP 'pz is not conserved'
-        ENDIF
-      ENDIF
-      DO I=0,3
-        DO J=1,(NLOOPLINE-1)
-          DO K=1,J
-            PCT(I,J)=PCT(I,J)+PL(I,K)
-          ENDDO
-        ENDDO
-      ENDDO
-
-C     Now convert the loop momenta to Golem95 conventions
-      DO I=0,3
-        PGOLEM(1,I)=0.0E0_KI
-        DO J=2,NLOOPLINE
-          PGOLEM(J,I)=PCT(I,J-1)
-        ENDDO
-      ENDDO
-
-C     Fill in the kinematic s-matrix while taking care of on-shell
-C      limits.
-      CALL ML5_0_SETUP_KIN_MATRIX(NLOOPLINE,PGOLEM,M2LGOLEM)
-C     Construct the golem internal matrices derived from the kinetic
-C      one.
-      CALL PREPARESMATRIX()
-
-C     Fill in the golem coefficents and compute the loop
-      IF(GOLEM_RUN_MODE.EQ.2)THEN
-        RES_GOLEM = EVALUATE_B(ML5_0_GOLEM_LOOPNUM,PGOLEM,0,RANK)
-      ELSE
-        PASS_COEF_CHECK=.TRUE.
-        SELECT CASE(RANK)
-        CASE(0)
-        CALL ML5_0_FILL_GOLEM_COEFFS_0(LOOPCOEFS(0,SQSOINDEX,ID)
-     $   ,COEFFS0)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          COEFFS0_REC = ML5_0_GOLEM_LOOPNUM(NULL_VEC,0.0_KI)
-          PASS_COEF_CHECK=ML5_0_COMPARE_COEFS_0(COEFFS0,COEFFS0_REC)
-        ENDIF
-        CASE(1)
-        CALL ML5_0_FILL_GOLEM_COEFFS_1(LOOPCOEFS(0,SQSOINDEX,ID)
-     $   ,COEFFS1)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT1(ML5_0_GOLEM_LOOPNUM,COEFFS1_REC)
-          PASS_COEF_CHECK=ML5_0_COMPARE_COEFS_1(COEFFS1,COEFFS1_REC)
-        ENDIF
-        CASE(2)
-        CALL ML5_0_FILL_GOLEM_COEFFS_2(LOOPCOEFS(0,SQSOINDEX,ID)
-     $   ,COEFFS2)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT2(ML5_0_GOLEM_LOOPNUM,COEFFS2_REC)
-          PASS_COEF_CHECK=ML5_0_COMPARE_COEFS_2(COEFFS2,COEFFS2_REC)
-        ENDIF
-        CASE(3)
-        CALL ML5_0_FILL_GOLEM_COEFFS_3(LOOPCOEFS(0,SQSOINDEX,ID)
-     $   ,COEFFS3)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT3(ML5_0_GOLEM_LOOPNUM,COEFFS3_REC)
-          PASS_COEF_CHECK=ML5_0_COMPARE_COEFS_3(COEFFS3,COEFFS3_REC)
-        ENDIF
-        CASE(4)
-        CALL ML5_0_FILL_GOLEM_COEFFS_4(LOOPCOEFS(0,SQSOINDEX,ID)
-     $   ,COEFFS4)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT4(ML5_0_GOLEM_LOOPNUM,COEFFS4_REC)
-          PASS_COEF_CHECK=ML5_0_COMPARE_COEFS_4(COEFFS4,COEFFS4_REC)
-        ENDIF
-        CASE(5)
-        CALL ML5_0_FILL_GOLEM_COEFFS_5(LOOPCOEFS(0,SQSOINDEX,ID)
-     $   ,COEFFS5)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT5(ML5_0_GOLEM_LOOPNUM,COEFFS5_REC)
-          PASS_COEF_CHECK=ML5_0_COMPARE_COEFS_5(COEFFS5,COEFFS5_REC)
-        ENDIF
-        CASE(6)
-        CALL ML5_0_FILL_GOLEM_COEFFS_6(LOOPCOEFS(0,SQSOINDEX,ID)
-     $   ,COEFFS6)
-        IF(GOLEM_RUN_MODE.EQ.3)THEN
-          CALL RECONSTRUCT6(ML5_0_GOLEM_LOOPNUM,COEFFS6_REC)
-          PASS_COEF_CHECK=ML5_0_COMPARE_COEFS_6(COEFFS6,COEFFS6_REC)
-        ENDIF
-        CASE DEFAULT
-        WRITE(*,*)'Not yet implemented in Golem95 for rank= ',RANK
-        STOP
-        END SELECT
-
-        IF(.NOT.PASS_COEF_CHECK)THEN
-          WRITE(*,*)'Coefs mismatch for ID ',ID,' and rank ',RANK
-          WRITE(*,*)'Coefs form MadLoop5:'
-          SELECT CASE(RANK)
-          CASE(0)
-          WRITE(*,*)'Constant coef = ',COEFFS0
-          CASE(1)
-          CALL PRINT_COEFFS(COEFFS1)
-          CASE(2)
-          CALL PRINT_COEFFS(COEFFS2)
-          CASE(3)
-          CALL PRINT_COEFFS(COEFFS3)
-          CASE(4)
-          CALL PRINT_COEFFS(COEFFS4)
-          CASE(5)
-          CALL PRINT_COEFFS(COEFFS5)
-          CASE(6)
-          CALL PRINT_COEFFS(COEFFS6)
-          END SELECT
-          WRITE(*,*)'Coefs reconstructed by Golem95:'
-          SELECT CASE(RANK)
-          CASE(0)
-          WRITE(*,*)'Constant coef = ',COEFFS0_REC
-          CASE(1)
-          CALL PRINT_COEFFS(COEFFS1_REC)
-          CASE(2)
-          CALL PRINT_COEFFS(COEFFS2_REC)
-          CASE(3)
-          CALL PRINT_COEFFS(COEFFS3_REC)
-          CASE(4)
-          CALL PRINT_COEFFS(COEFFS4_REC)
-          CASE(5)
-          CALL PRINT_COEFFS(COEFFS5_REC)
-          CASE(6)
-          CALL PRINT_COEFFS(COEFFS6_REC)
-          END SELECT
-          STOP
-        ENDIF
-
-        SELECT CASE(NLOOPLINE)
-        CASE(1)
-        WRITE(*,*)'Golem95 cannot handle with tadpole yet'
-        STOP
-        CASE(2)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A20(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT2_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT2_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT2_3(COEFFS3,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',2,RANK
-        STOP
-        END SELECT
-        CASE(3)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A30(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT3_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT3_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT3_3(COEFFS3,PGOLEM,PINCH)
-        CASE(4)
-        RES_GOLEM = CONTRACT3_4(COEFFS4,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',3,RANK
-        STOP
-        END SELECT
-        CASE(4)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A40(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT4_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT4_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT4_3(COEFFS3,PGOLEM,PINCH)
-        CASE(4)
-        RES_GOLEM = CONTRACT4_4(COEFFS4,PGOLEM,PINCH)
-        CASE(5)
-        RES_GOLEM = CONTRACT4_5(COEFFS5,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',4,RANK
-        STOP
-        END SELECT
-        CASE(5)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A50(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT5_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT5_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT5_3(COEFFS3,PGOLEM,PINCH)
-        CASE(4)
-        RES_GOLEM = CONTRACT5_4(COEFFS4,PGOLEM,PINCH)
-        CASE(5)
-        RES_GOLEM = CONTRACT5_5(COEFFS5,PGOLEM,PINCH)
-        CASE(6)
-        RES_GOLEM = CONTRACT5_6(COEFFS6,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',5,RANK
-        STOP
-        END SELECT
-        CASE(6)
-        SELECT CASE(RANK)
-        CASE(0)
-        RES_GOLEM = COEFFS0*A60(PINCH)
-        CASE(1)
-        RES_GOLEM = CONTRACT6_1(COEFFS1,PGOLEM,PINCH)
-        CASE(2)
-        RES_GOLEM = CONTRACT6_2(COEFFS2,PGOLEM,PINCH)
-        CASE(3)
-        RES_GOLEM = CONTRACT6_3(COEFFS3,PGOLEM,PINCH)
-        CASE(4)
-        RES_GOLEM = CONTRACT6_4(COEFFS4,PGOLEM,PINCH)
-        CASE(5)
-        RES_GOLEM = CONTRACT6_5(COEFFS5,PGOLEM,PINCH)
-        CASE(6)
-        RES_GOLEM = CONTRACT6_6(COEFFS6,PGOLEM,PINCH)
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N,r = ',6,RANK
-        STOP
-        END SELECT
-        CASE DEFAULT
-        WRITE(*,*)'Golem95 cannot handle with: N = ',NLOOPLINE
-        STOP
-        END SELECT
-      ENDIF
-
-      RES(1)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%C+2.0*LOG(MU_R)
-     $ *RES_GOLEM%%B+2.0*LOG(MU_R)**2*RES_GOLEM%%A)
-      RES(2)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%B+2.0*LOG(MU_R)
-     $ *RES_GOLEM%%A)
-      RES(3)=NORMALIZATION*2.0D0*DBLE(RES_GOLEM%%A)
-C     WRITE(*,*) 'Loop ID',ID,' =',RES(1),RES(2),RES(3)
-
-C     Finally free golem memory and cache
-      CALL EXITGOLEM95()
-
-      END
-
-      FUNCTION ML5_0_GOLEM_LOOPNUM(Q,MU2)
-      USE PRECISION_GOLEM, ONLY: KI
-      REAL(KI),DIMENSION(0:3),INTENT(IN)::Q
-      REAL(KI),INTENT(IN)::MU2
-      COMPLEX(KI)::ML5_0_GOLEM_LOOPNUM
-
-      COMPLEX*16 QQ(0:3),NUM
-      INTEGER I
-
-      DO I=0,3
-        QQ(I)=CMPLX(Q(I),0.0D0,KIND=16)
-      ENDDO
-
-      CALL ML5_0_LOOPNUM(QQ,NUM)
-      ML5_0_GOLEM_LOOPNUM=NUM
-      RETURN
-      END FUNCTION
-
-      SUBROUTINE ML5_0_INITGOLEM()
-C     
-C     INITIALISATION OF GOLEM
-C     
-C     
-C     MODULE
-C     
-      USE PARAMETRE
-C     
-C     LOCAL VARIABLES 
-C     
-      REAL*8 THRS
-      LOGICAL EXT_NUM_FOR_R1
-C     
-C     GLOBAL VARIABLES 
-C     
-      INCLUDE 'MadLoopParams.inc'
-C     ----------
-C     BEGIN CODE
-C     ----------
-
-C     DEFAULT PARAMETERS FOR GOLEM
-C     -------------------------------  
-C     One can chose here to have either just the rational R1 piece 
-C     or everything but the R2
-      RAT_OR_TOT_PAR = TOT
-
-      END
-
-      SUBROUTINE ML5_0_SETUP_KIN_MATRIX(NLOOPLINE,PGOLEM,M2L)
-C     
-C     MODULE
-C     
-      USE MATRICE_S
-      USE PRECISION_GOLEM, ONLY: KI
-C     
-C     ARGUMENTS
-C     
-      INTEGER NLOOPLINE
-      REAL(KI) PGOLEM(NLOOPLINE,0:3)
-      COMPLEX(KI) M2L(NLOOPLINE)
-C     
-C     LOCAL VARIABLES
-C     
-      INTEGER I,J
-      COMPLEX*16 S_MAT_FROM_MG(NLOOPLINE,NLOOPLINE)
-C     ----------
-C     BEGIN CODE
-C     ----------
-
-      CALL ML5_0_BUILD_KINEMATIC_MATRIX(NLOOPLINE,PGOLEM,M2L
-     $ ,S_MAT_FROM_MG)
-
-      DO I=1,NLOOPLINE
-        DO J=1,NLOOPLINE
-          S_MAT(I,J)=S_MAT_FROM_MG(I,J)
-        ENDDO
-      ENDDO
-
-      END
-
-      FUNCTION ML5_0_COMPARE_COEFS_0(COEFS_A,COEFS_B)
-
-      USE PRECISION_GOLEM, ONLY: KI
-      COMPLEX(KI) COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/ML5_0_COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL ML5_0_COMPARE_COEFS_0
-
-      NUM = ABS(COEFS_A-COEFS_B)
-      DENOM = ABS(COEFS_A+COEFS_B)
-      IF(DENOM.GT.0D0)THEN
-        ML5_0_COMPARE_COEFS_0=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        ML5_0_COMPARE_COEFS_0=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION ML5_0_COMPARE_COEFS_1(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_1)COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/ML5_0_COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL ML5_0_COMPARE_COEFS_1
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )
-
-      IF(DENOM.GT.0D0)THEN
-        ML5_0_COMPARE_COEFS_1=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        ML5_0_COMPARE_COEFS_1=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION ML5_0_COMPARE_COEFS_2(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_2) COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/ML5_0_COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL ML5_0_COMPARE_COEFS_2
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))
-      IF(DENOM.GT.0D0)THEN
-        ML5_0_COMPARE_COEFS_2=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        ML5_0_COMPARE_COEFS_2=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION ML5_0_COMPARE_COEFS_3(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_3) COEFS_A, COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/ML5_0_COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL ML5_0_COMPARE_COEFS_3
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3)
-     $ )
-      IF(DENOM.GT.0D0)THEN
-        ML5_0_COMPARE_COEFS_3=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        ML5_0_COMPARE_COEFS_3=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION ML5_0_COMPARE_COEFS_4(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_4) COEFS_A, COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/ML5_0_COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL ML5_0_COMPARE_COEFS_4
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
-     $ +SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3)
-     $ )+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
-      IF(DENOM.GT.0D0)THEN
-        ML5_0_COMPARE_COEFS_4=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        ML5_0_COMPARE_COEFS_4=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION ML5_0_COMPARE_COEFS_5(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_5) COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/ML5_0_COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL ML5_0_COMPARE_COEFS_5
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
-     $ +SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3)
-     $ )+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
-      IF(DENOM.GT.0D0)THEN
-        ML5_0_COMPARE_COEFS_5=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        ML5_0_COMPARE_COEFS_5=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-      FUNCTION ML5_0_COMPARE_COEFS_6(COEFS_A,COEFS_B)
-
-      USE TENS_REC
-      TYPE(COEFF_TYPE_6) COEFS_A,COEFS_B
-      REAL*8 COEF_CHECK_THRS
-      COMMON/ML5_0_COEF_CHECK_THRS/COEF_CHECK_THRS
-      REAL*8 DENOM,NUM
-      LOGICAL ML5_0_COMPARE_COEFS_6
-
-      NUM = ABS(COEFS_A%%C0-COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1-COEFS_B%%C1))
-     $ +SUM(ABS(COEFS_A%%C2-COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3-COEFS_B%%C3))
-     $ +SUM(ABS(COEFS_A%%C4-COEFS_B%%C4))
-      DENOM = ABS(COEFS_A%%C0+COEFS_B%%C0)+SUM(ABS(COEFS_A%%C1+COEFS_B%%C1)
-     $ )+SUM(ABS(COEFS_A%%C2+COEFS_B%%C2))+SUM(ABS(COEFS_A%%C3+COEFS_B%%C3)
-     $ )+SUM(ABS(COEFS_A%%C4+COEFS_B%%C4))
-      IF(DENOM.GT.0D0)THEN
-        ML5_0_COMPARE_COEFS_6=((NUM/DENOM).LT.COEF_CHECK_THRS)
-      ELSE
-        ML5_0_COMPARE_COEFS_6=(NUM.LT.COEF_CHECK_THRS)
-      ENDIF
-
-      END
-
-
-      SUBROUTINE ML5_0_FILL_GOLEM_COEFFS_0(ML_COEFS,GOLEM_COEFS)
-      USE PRECISION_GOLEM, ONLY: KI
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      COMPLEX(KI) GOLEM_COEFS
-      GOLEM_COEFS=ML_COEFS(0)
-      END
-
-      SUBROUTINE ML5_0_FILL_GOLEM_COEFFS_1(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_1
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_1) GOLEM_COEFS
-C     Constant coefficient 
-      GOLEM_COEFS%%C0=ML_COEFS(0)
-C     Coefficient q(0)
-      GOLEM_COEFS%%C1(1,1)=-ML_COEFS(1)
-C     Coefficient q(1)
-      GOLEM_COEFS%%C1(2,1)=-ML_COEFS(2)
-C     Coefficient q(2)
-      GOLEM_COEFS%%C1(3,1)=-ML_COEFS(3)
-C     Coefficient q(3)
-      GOLEM_COEFS%%C1(4,1)=-ML_COEFS(4)
-      END
-
-      SUBROUTINE ML5_0_FILL_GOLEM_COEFFS_2(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_2
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_2) GOLEM_COEFS
-C     Constant coefficient 
-      GOLEM_COEFS%%C0=ML_COEFS(0)
-C     Coefficient q(0)
-      GOLEM_COEFS%%C1(1,1)=-ML_COEFS(1)
-C     Coefficient q(0)^2
-      GOLEM_COEFS%%C1(1,2)= ML_COEFS(5)
-C     Coefficient q(1)
-      GOLEM_COEFS%%C1(2,1)=-ML_COEFS(2)
-C     Coefficient q(1)^2
-      GOLEM_COEFS%%C1(2,2)= ML_COEFS(7)
-C     Coefficient q(2)
-      GOLEM_COEFS%%C1(3,1)=-ML_COEFS(3)
-C     Coefficient q(2)^2
-      GOLEM_COEFS%%C1(3,2)= ML_COEFS(10)
-C     Coefficient q(3)
-      GOLEM_COEFS%%C1(4,1)=-ML_COEFS(4)
-C     Coefficient q(3)^2
-      GOLEM_COEFS%%C1(4,2)= ML_COEFS(14)
-C     Coefficient q(0)*q(1)
-      GOLEM_COEFS%%C2(1,1)= ML_COEFS(6)
-C     Coefficient q(0)*q(2)
-      GOLEM_COEFS%%C2(2,1)= ML_COEFS(8)
-C     Coefficient q(0)*q(3)
-      GOLEM_COEFS%%C2(3,1)= ML_COEFS(11)
-C     Coefficient q(1)*q(2)
-      GOLEM_COEFS%%C2(4,1)= ML_COEFS(9)
-C     Coefficient q(1)*q(3)
-      GOLEM_COEFS%%C2(5,1)= ML_COEFS(12)
-C     Coefficient q(2)*q(3)
-      GOLEM_COEFS%%C2(6,1)= ML_COEFS(13)
-      END
-
-      SUBROUTINE ML5_0_FILL_GOLEM_COEFFS_3(ML_COEFS,GOLEM_COEFS)
-      USE TENS_REC, ONLY: COEFF_TYPE_3
-      INCLUDE 'coef_specs.inc'
-      INCLUDE 'loop_max_coefs.inc'
-      COMPLEX*16 ML_COEFS(0:LOOPMAXCOEFS-1)
-      TYPE(COEFF_TYPE_3) GOLEM_COEFS
-C     Constant coefficient 
-      GOLEM_COEFS%%C0=ML_COEFS(0)
-C     Coefficient q(0)
-      GOLEM_COEFS%%C1(1,1)=-ML_COEFS(1)
-C     Coefficient q(0)^2
-      GOLEM_COEFS%%C1(1,2)= ML_COEFS(5)
-C     Coefficient q(0)^3
-      GOLEM_COEFS%%C1(1,3)=-ML_COEFS(15)
-C     Coefficient q(1)
-      GOLEM_COEFS%%C1(2,1)=-ML_COEFS(2)
-C     Coefficient q(1)^2
-      GOLEM_COEFS%%C1(2,2)= ML_COEFS(7)
-C     Coefficient q(1)^3
-      GOLEM_COEFS%%C1(2,3)=-ML_COEFS(18)
-C     Coefficient q(2)
-      GOLEM_COEFS%%C1(3,1)=-ML_COEFS(3)
-C     Coefficient q(2)^2
-      GOLEM_COEFS%%C1(3,2)= ML_COEFS(10)
-C     Coefficient q(2)^3
-      GOLEM_COEFS%%C1(3,3)=-ML_COEFS(24)
-C     Coefficient q(3)
-      GOLEM_COEFS%%C1(4,1)=-ML_COEFS(4)
-C     Coefficient q(3)^2
-      GOLEM_COEFS%%C1(4,2)= ML_COEFS(14)
-C     Coefficient q(3)^3
-      GOLEM_COEFS%%C1(4,3)=-ML_COEFS(34)
-C     Coefficient q(0)*q(1)
-      GOLEM_COEFS%%C2(1,1)= ML_COEFS(6)
-C     Coefficient q(0)*q(1)^2
-      GOLEM_COEFS%%C2(1,2)=-ML_COEFS(17)
-C     Coefficient q(0)^2*q(1)
-      GOLEM_COEFS%%C2(1,3)=-ML_COEFS(16)
-C     Coefficient q(0)*q(2)
-      GOLEM_COEFS%%C2(2,1)= ML_COEFS(8)
-C     Coefficient q(0)*q(2)^2
-      GOLEM_COEFS%%C2(2,2)=-ML_COEFS(22)
-C     Coefficient q(0)^2*q(2)
-      GOLEM_COEFS%%C2(2,3)=-ML_COEFS(19)
-C     Coefficient q(0)*q(3)
-      GOLEM_COEFS%%C2(3,1)= ML_COEFS(11)
-C     Coefficient q(0)*q(3)^2
-      GOLEM_COEFS%%C2(3,2)=-ML_COEFS(31)
-C     Coefficient q(0)^2*q(3)
-      GOLEM_COEFS%%C2(3,3)=-ML_COEFS(25)
-C     Coefficient q(1)*q(2)
-      GOLEM_COEFS%%C2(4,1)= ML_COEFS(9)
-C     Coefficient q(1)*q(2)^2
-      GOLEM_COEFS%%C2(4,2)=-ML_COEFS(23)
-C     Coefficient q(1)^2*q(2)
-      GOLEM_COEFS%%C2(4,3)=-ML_COEFS(21)
-C     Coefficient q(1)*q(3)
-      GOLEM_COEFS%%C2(5,1)= ML_COEFS(12)
-C     Coefficient q(1)*q(3)^2
-      GOLEM_COEFS%%C2(5,2)=-ML_COEFS(32)
-C     Coefficient q(1)^2*q(3)
-      GOLEM_COEFS%%C2(5,3)=-ML_COEFS(27)
-C     Coefficient q(2)*q(3)
-      GOLEM_COEFS%%C2(6,1)= ML_COEFS(13)
-C     Coefficient q(2)*q(3)^2
-      GOLEM_COEFS%%C2(6,2)=-ML_COEFS(33)
-C     Coefficient q(2)^2*q(3)
-      GOLEM_COEFS%%C2(6,3)=-ML_COEFS(30)
-C     Coefficient q(0)*q(1)*q(2)
-      GOLEM_COEFS%%C3(1,1)=-ML_COEFS(20)
-C     Coefficient q(0)*q(1)*q(3)
-      GOLEM_COEFS%%C3(2,1)=-ML_COEFS(26)
-C     Coefficient q(0)*q(2)*q(3)
-      GOLEM_COEFS%%C3(3,1)=-ML_COEFS(28)
-C     Coefficient q(1)*q(2)*q(3)
-      GOLEM_COEFS%%C3(4,1)=-ML_COEFS(29)