Merge lp:~maddevelopers/mg5amcnlo/2.7.3_tflipping into lp:mg5amcnlo/lts

[Olivier Mattelaer]

implement helicity recycling at LO and new phase-space integration strategy

This branch implements some significant speed-up to the LO code

with two types of speed-up
1) faster method to evaluation tree-level (for LO) matrix-element --same can be implemented at NLO--
2) improved phase-space integration method.

The Paper is close to be finalized.

The new phase-space improves (even more since 2.8.0 started to improve it already) the speed mainly for VBF processes. This branch has been extensively used for muon collider study where the previous version of the code was not even converging.

The modification of the matrix-element provides a factor of ~2 in speed for complicated process (like g g > t t~ g g) but does not improves the situation much for very simple process.

This will be released as 2.9.0

[marco zaro (marco-zaro)]

Thanks Olivier,
how did you achieve the speed improvement for the matrix elements?
Cheers,

Marco

On Sun, Dec 6, 2020 at 10:08 PM Olivier Mattelaer <
<email address hidden>> wrote:

> Olivier Mattelaer has proposed merging
> lp:~maddevelopers/mg5amcnlo/2.7.3_tflipping into lp:mg5amcnlo.
>
> Commit message:
> implement helicity recycling at LO and new phase-space integration strategy
>
> Requested reviews:
> MadTeam (madteam)
>
> For more details, see:
>
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.7.3_tflipping/+merge/394923
>
> This branch implements some significant speed-up to the LO code
>
> with two types of speed-up
> 1) faster method to evaluation tree-level (for LO) matrix-element --same
> can be implemented at NLO--
> 2) improved phase-space integration method.
>
> The Paper is close to be finalized.
>
>
> The new phase-space improves (even more since 2.8.0 started to improve it
> already) the speed mainly for VBF processes. This branch has been
> extensively used for muon collider study where the previous version of the
> code was not even converging.
>
> The modification of the matrix-element provides a factor of ~2 in speed
> for complicated process (like g g > t t~ g g) but does not improves the
> situation much for very simple process.
>
>
> This will be released as 2.9.0
>
> --
> Your team MadDevelopers is subscribed to branch
> lp:~maddevelopers/mg5amcnlo/2.7.3_tflipping.
>

[Olivier Mattelaer (olivier-mattelaer)]

This is the method that I describe in our last meeting and you have one talk about it here:
https://indico.cern.ch/event/955391/contributions/4075843/attachments/2137914/3601421/beamer.pdf

To sum up the initial idea is to remove the explicit loop over helicty and to implement a recycling avoiding to compute twice the same term (in the same spirit as the current recycling between Feynman Diagram).
This also allows up to split the amplitude routine FFV1_0 in two successive call which the first one can be recycle and the second one being a very simple scalar product of two complex vector. This is crucial since we are dominated by the amplitude function in general.

This can (and should) be implemented for the real in a later stage

[marco zaro (marco-zaro)]

Right, I remember now, thanks a lot!

Cheers,

Marco

> On 9 Dec 2020, at 10:50, Olivier Mattelaer <email address hidden> wrote:
>
> This is the method that I describe in our last meeting and you have one talk about it here:
> https://indico.cern.ch/event/955391/contributions/4075843/attachments/2137914/3601421/beamer.pdf
>
> To sum up the initial idea is to remove the explicit loop over helicty and to implement a recycling avoiding to compute twice the same term (in the same spirit as the current recycling between Feynman Diagram).
> This also allows up to split the amplitude routine FFV1_0 in two successive call which the first one can be recycle and the second one being a very simple scalar product of two complex vector. This is crucial since we are dominated by the amplitude function in general.
>
> This can (and should) be implemented for the real in a later stage
>
> --
> https://code.launchpad.net/~maddevelopers/mg5amcnlo/2.7.3_tflipping/+merge/394923
> Your team MadTeam is requested to review the proposed merge of lp:~maddevelopers/mg5amcnlo/2.7.3_tflipping into lp:mg5amcnlo.

[lp:~maddevelopers/mg5amcnlo/2.7.3_tflipping]

336. By olivier-mattelaer on 2020-12-10

change option name in the run-card

337. By olivier-mattelaer on 2020-12-10

fix logger debugging issue

338. By olivier-mattelaer on 2020-12-16

fix bug reported by Luca

339. By olivier-mattelaer on 2020-12-20

include code for better jamp computation (not default)

340. By olivier-mattelaer on 2020-12-22

pass to exponential suppression for tmin

341. By olivier-mattelaer on 2020-12-23

fix issue with the auto-setting of SDE

342. By olivier-mattelaer on 2020-12-31

fix s-channel suppression factor

343. By olivier-mattelaer on 2021-01-02

set jamp_optim as default + remove denom

344. By olivier-mattelaer on 2021-01-07

better color matrix + fix to t channel ordering

345. By Olivier Mattelaer on 2021-01-07

add t_strategy -1 in the mix

346. By Olivier Mattelaer on 2021-01-08

display sde_strat in the run_card, add option for more ps parameter and better default for sde_strat, and fix issue with update frame

347. By olivier-mattelaer on 2021-01-08

improve smart t-channel ordering

348. By olivier-mattelaer on 2021-01-10

fix ordering related issue in gensym and other

349. By olivier-mattelaer on 2021-01-11

fixing SDE default and fixing MLM/CKKW run_card issue

350. By olivier-mattelaer on 2021-01-14

update the UpdateNote

351. By olivier-mattelaer on 2021-01-17

allow to specify fortran flag via the run_card

352. By olivier-mattelaer on 2021-01-19

merge with latest 2.8.3.py3

353. By olivier-mattelaer on 2021-01-21

fix a lot of tests

354. By olivier-mattelaer on 2021-01-21

merge with 2.8.3

355. By olivier-mattelaer on 2021-01-21

merge with latest 2.8.3

356. By olivier-mattelaer on 2021-01-23

fixing more tests

357. By olivier-mattelaer on 2021-01-23

fix issue with default make_opts

358. By olivier-mattelaer on 2021-01-26

set founds-check for debug mode by default

359. By olivier-mattelaer on 2021-01-26

merge with official 2.8.3

360. By olivier-mattelaer on 2021-01-26

Update Notes update

361. By olivier-mattelaer on 2021-01-27

prevent future issue with renamed directory regexp

362. By olivier-mattelaer on 2021-01-28

fix python2 issue --integer division--, an issue with elastic photon and systematics, and some unittests

363. By olivier-mattelaer on 2021-01-28

merge with some 2.8.3.py2

364. By Kiran Ostrolenk on 2021-01-29

Allow for non-sequential calls to NHEL

365. By olivier-mattelaer on 2021-01-29

update IO test

366. By olivier-mattelaer on 2021-01-29

add an options to force the model conversion to py3

367. By olivier-mattelaer on 2021-01-30

fixing one more tests

368. By olivier-mattelaer on 2021-01-30

add functionality for adding information/message from online

=== modified file 'Template/LO/Cards/run_card.dat'
--- Template/LO/Cards/run_card.dat	2020-08-07 14:49:25 +0000
+++ Template/LO/Cards/run_card.dat	2021-01-30 10:13:29 +0000
@@ -58,16 +58,19 @@
   %(gridpack)s     = gridpack  !True = setting up the grid pack
   %(time_of_flight)s = time_of_flight ! threshold (in mm) below which the invariant livetime is not written (-1 means not written)
   %(event_norm)s =  event_norm       ! average/sum. Normalization of the weight in the LHEF
-$MLM
-$CKKW
+$mlm
+$ckkw
 #*********************************************************************
 #
 #*********************************************************************
-# handling of the helicities:
-#  0: sum over all helicities
-#  1: importance sampling over helicities
+# Phase-Space Optimization strategy (basic options)
 #*********************************************************************
    %(nhel)s  = nhel          ! using helicities importance sampling or not.
+                             ! 0: sum over helicity, 1: importance sampling
+   %(sde_strategy)s  = sde_strategy  ! default integration strategy (hep-ph/2021.xxxxx)
+                             ! 1 is old strategy (using amp square)
+			     ! 2 is new strategy (using only the denominator)
+$psoptim			     
 #*********************************************************************
 # Generation bias, check the wiki page below for more information:   *
 #  'cp3.irmp.ucl.ac.be/projects/madgraph/wiki/LOEventGenerationBias' *

=== modified file 'Template/LO/Source/.make_opts'
--- Template/LO/Source/.make_opts	2019-01-23 16:03:38 +0000
+++ Template/LO/Source/.make_opts	2021-01-30 10:13:29 +0000
@@ -13,10 +13,11 @@
 
 # Rest of the makefile
 ifeq ($(origin FFLAGS),undefined)
-FFLAGS= -O -w -fbounds-check -fPIC
+FFLAGS= -w -fPIC
 #FFLAGS+= -g -fbounds-check -ffpe-trap=invalid,zero,overflow,underflow,denormal -Wall
 endif
 
+FFLAGS += $(GLOBAL_FLAG)
 
 # REMOVE MACFLAG IF NOT ON MAC OR FOR F2PY
 UNAME := $(shell uname -s)

=== modified file 'Template/LO/Source/PDF/pdf.f'
--- Template/LO/Source/PDF/pdf.f	2020-08-07 14:49:25 +0000
+++ Template/LO/Source/PDF/pdf.f	2021-01-30 10:13:29 +0000
@@ -297,9 +297,14 @@
 c  a "diffractive" photon
 c      
       q2max=xmu*xmu
-      if(abs(ih) .eq. 3.or.abs(ih) .eq. 4) then  !from the electron
-          fx(7)=epa_lepton(x,q2max, ih)
-      elseif(ih .eq. 2) then  !from a proton without breaking
+      if(abs(ih) .eq. 3.or.abs(ih) .eq. 4) then !from the electron
+         write(*,*) 'impossible call (or was it) to pdf-> please reporrt'
+         stop 23
+         fx(7)=epa_lepton(x,q2max, ih)
+      elseif(ih .eq. 2) then    !from a proton without breaking
+         write(*,*) 'impossible call (or was it) to pdf-> please reporrt'
+         stop 23
+         ! isssue with next call since this function takes now a third argument (beamid not define here)
           fx(7)=epa_proton(x,q2max)
       endif      
       

=== removed file 'Template/LO/Source/PDF/pdf_lhapdf.f'
=== modified file 'Template/LO/Source/dsample.f'
--- Template/LO/Source/dsample.f	2020-09-09 12:43:23 +0000
+++ Template/LO/Source/dsample.f	2021-01-30 10:13:29 +0000
@@ -257,6 +257,7 @@
 c     Need to start from scratch. This is clunky but I'll just
 c     remove the grid, so we are clean
 c
+      goto 200
       write(*,*) "Trying w/ fresh grid"
       open(unit=25,file='ftn25',status='unknown',err=102)
       write(25,*) ' '
@@ -329,7 +330,7 @@
 c
 c     All done
 c
-      open(unit=66,file='results.dat',status='unknown')
+200   open(unit=66,file='results.dat',status='unknown')
       i=1
       do while(xmean(i) .ne. 0 .and. i .lt. cur_it)
          i=i+1
@@ -2274,7 +2275,11 @@
       integer                                      nsteps
       character*40          result_file,where_file
       common /sample_status/result_file,where_file,nsteps
-
+c
+c     
+c
+       logical init_mode
+       common/to_determine_zero_hel/init_mode
 c----
 c  Begin Code
 c----
@@ -2282,6 +2287,9 @@
          write(*,*) nb_pass_cuts, 
      &    ' points passed the cut but all returned zero'
          write(*,*) 'therefore considering this contribution as zero'
+         if (init_mode) then
+            call print_zero_amp()
+         endif
       else if (nb_pass_cuts.gt.0.and.nb_pass_cuts.lt.1000)then
          write(*,*) 'only', nb_pass_cuts, 
      &    ' points passed the cut and they all returned zero'

=== modified file 'Template/LO/Source/invarients.f'
--- Template/LO/Source/invarients.f	2011-04-16 05:25:14 +0000
+++ Template/LO/Source/invarients.f	2021-01-30 10:13:29 +0000
@@ -185,10 +185,10 @@
       end
 
 
-      subroutine map_invarients(Minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
+      subroutine map_invarients(Minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming,nb_tchannel)
 c****************************************************************************
 c     Determines mappings for each structure of invarients onto integration
-c     variables.  Input: Ninvar, iforest.  Output: Minvar, ninvar
+c     variables.  Input: Ninvar, iforest.  Output: Minvar, ninvar, nb_tchannel
 c****************************************************************************
       implicit none
 c
@@ -201,11 +201,13 @@
 c
       integer Minvar(maxdim,lmaxconfigs),nconfigs,ninvar,nexternal,nincoming
       integer mincfig,maxcfig
+      integer nb_tchannel
 c
 c     Local
 c
       integer iconfig, jgrid,j, nbranch
       logical found,tchannel
+      integer ns_channel
 c
 c     Global
 c
@@ -224,20 +226,30 @@
 c     
 c     Try simple mapping if nconfigs = 1
 c
-
+      nb_tchannel = 0
       if (nconfigs .eq. 1) then
 c         do j=1,3*nbranch-4+2
          do j=1,maxdim
             minvar(j,mincfig)=j
          enddo
          jgrid=j-1
+
+         write(*,*) 'Determine nb_t'
+         ns_channel=1
+         do while((iforest(1,-ns_channel,mincfig) .ne. 1.and.iforest(1,-ns_channel,mincfig) .ne. 2).and.ns_channel.lt.nbranch)
+            ns_channel=ns_channel+1
+         enddo
+         ns_channel=ns_channel - 1
+         nb_tchannel=nbranch-ns_channel-1
+         write(*,*) 'T-channel found: ',nb_tchannel 
+         
       else
 c      if (ep) jgrid=1
 c      if (pp) jgrid=2
       do iconfig=mincfig,maxcfig
          tchannel = .false.         
          do j=1,nbranch-1
-            if (iforest(1,-j,iconfig) .eq. 1) then
+            if (iforest(1,-j,iconfig) .eq. 1.or.(nincoming.eq.2.and.iforest(1,-j,iconfig) .eq. 2)) then
                tchannel=.true.
             endif
             jgrid=jgrid+1

=== modified file 'Template/LO/Source/run.inc'
--- Template/LO/Source/run.inc	2019-03-23 22:04:23 +0000
+++ Template/LO/Source/run.inc	2021-01-30 10:13:29 +0000
@@ -97,4 +97,9 @@
      &etamax4pdg, mxxmin4pdg,mxxpart_antipart
 
        double precision small_width_treatment
-       common/narrow_width/small_width_treatment
\ No newline at end of file
+       common/narrow_width/small_width_treatment
+
+       double precision tmin_for_channel
+       integer sde_strat ! 1 means standard single diagram enhancement strategy,
+c                          2 means approximation by the denominator of the propagator
+       common/TO_CHANNEL_STRAT/tmin_for_channel, sde_strat
\ No newline at end of file

=== modified file 'Template/LO/SubProcesses/cluster.f'
--- Template/LO/SubProcesses/cluster.f	2013-03-21 12:45:09 +0000
+++ Template/LO/SubProcesses/cluster.f	2021-01-30 10:13:29 +0000
@@ -210,7 +210,8 @@
       integer ignum, ipnum, ipids(nexternal,4,2:nexternal)
 C $B$ IFOREST $B$ !this is a tag for MadWeight
       integer i, iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/iforest, tstrategy
       integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
       integer tprid(-max_branch:-1,lmaxconfigs)
       common/to_sprop/sprop,tprid
@@ -404,7 +405,8 @@
       common/to_mconfigs/mapconfig, this_config
 C $B$ IFOREST $B$ !this is a tag for MadWeight
       integer i, iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/iforest,	tstrategy
 C $E$ IFOREST $E$ !this is a tag for MadWeight
 C $B$ DECAYBW $E$ !this is a tag for MadWeight
 

=== modified file 'Template/LO/SubProcesses/genps.f'
--- Template/LO/SubProcesses/genps.f	2020-08-10 19:06:52 +0000
+++ Template/LO/SubProcesses/genps.f	2021-01-30 10:13:29 +0000
@@ -163,7 +163,8 @@
 
 
       integer iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/ iforest, tstrategy
 
       integer            mapconfig(0:lmaxconfigs), this_config
       common/to_mconfigs/mapconfig, this_config
@@ -321,7 +322,7 @@
       endif
       pswgt = 1d0
       jac   = 1d0
-      call one_tree(iforest(1,-max_branch,iconfig),mincfig,
+      call one_tree(iforest(1,-max_branch,iconfig), tstrategy(iconfig),mincfig,
      &     nbranch,P,M,S,X,jac,pswgt)
 c
 c     Add what I think are the essentials
@@ -520,7 +521,7 @@
       integer nbranch,ndim,nconfigs
       integer ninvar
       integer nparticles,nfinal
-
+      integer nb_tchannel
 
 c
 c     Arguments
@@ -595,7 +596,7 @@
          nconfigs = 1
          mincfig=iconfig
          maxcfig=iconfig
-         call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
+         call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming,nb_tchannel)
          maxwgt=0d0
          nparticles   = nexternal
          nfinal       = nparticles-nincoming
@@ -632,7 +633,7 @@
       return
       end
 
-      subroutine one_tree(itree,iconfig,nbranch,P,M,S,X,jac,pswgt)
+      subroutine one_tree(itree,tstrategy,iconfig,nbranch,P,M,S,X,jac,pswgt)
 c************************************************************************
 c     Calculates the momentum for everything below in the tree until
 c     it reaches the end.
@@ -652,8 +653,10 @@
 c     Arguments
 c
       integer itree(2,-max_branch:-1) !Structure of configuration
+      integer tstrategy ! current strategy for t-channel
       integer iconfig                 !Which configuration working on
       double precision P(0:3,-max_branch:max_particles)
+      double precision pother(0:3), ptemp(0:3), pboost(0:3), ptemp2(0:3)
       double precision M(-max_branch:max_particles)
       double precision S(-max_branch:0)
 c      double precision spole(-max_branch:0),swidth(-max_branch:0)
@@ -664,12 +667,13 @@
 c     Local
 c
       logical pass
+      double precision tmass(-max_branch:-1)
       integer ibranch,i,ns_channel,nt_channel,ix  !,nerr
-c      data nerr/0/
+      integer iopposite ! index for t-channel mapping for the part not handle by itree
+c     data nerr/0/
       double precision smin,smax,totmass,totmassin,xa2,xb2,wgt
       double precision costh,phi,tmin,tmax,t
       double precision ma2,mb2,m12,mn2,s1, mi2
-      double precision tmass(-max_branch:-1)
 c
 c     External
 c
@@ -703,7 +707,12 @@
 c     the s channel p1+p2
 c
       ns_channel=1
-      do while(itree(1,-ns_channel) .ne. 1 .and.ns_channel.lt.nbranch)
+      iopposite = 1
+      if (abs(tstrategy).eq.1) then
+         iopposite = 2
+      endif
+      
+      do while(itree(1,-ns_channel) .ne. iopposite .and.ns_channel.lt.nbranch)
          m(-ns_channel)=0d0                 
          ns_channel=ns_channel+1         
       enddo
@@ -713,7 +722,6 @@
       if (nt_channel .eq. 0 .and. nincoming .eq. 2) then
          ns_channel=ns_channel-1
       endif
-
 c
 c     Determine masses for all intermediate states.  Starting
 c     from outer most (real particle) states
@@ -811,7 +819,7 @@
          totmass=totmass-m(itree(2,ibranch))      !for remaining particles
          smin = totmass**2                        !This affects t_min/max
          smax = (m(ibranch) - m(itree(2,ibranch)))**2
-
+         
          if (smin .gt. smax) then
             jac=-3d0
             return
@@ -837,18 +845,239 @@
 c
 c     From here, on we can just pretend this is a 2->2 scattering with
 c     Pa                    + Pb     -> P1          + P2
+
+      if (tstrategy.eq.-2.or.tstrategy.eq.-1) then
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc  
+cc       T-channel ping-pong strategy starting with 2
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+         
+c     No -flipping case:      
+c      p(0,itree(ibranch,1)) + p(0,2) -> p(0,ibranch)+ p(0,itree(ibranch,2))
+c       -  M(ibranch) is the total mass available (Pa+Pb)^2
+c       - M(ibranch-1) is the mass of P2  (all the remaining particles)
+c
+c     With flipping case
+c      p(0,itree(ibranch,1)) + pother -> p(0,ibranch)+ p(0,itree(ibranch,2))
+c       - pother = p(0,itree(ibranch,1)) -p(0,itree(ibranch,2))         
+c       - M(ibranch) is the total mass available (Pa+Pb)^2
+c       - M(ibranch-1) is the mass of P2  (all the remaining particles)      
+c
+c     This assumes that P(0, ibranch) is set to the T-channel propa (likely)
+c      do ibranch = -ns_channel-1,-nbranch,-1
+c         totmass=totmass+m(itree(2,ibranch))
+c      enddo
+      do ibranch=-ns_channel-1,-nbranch+1,-1
+c         totmass=totmass-m(itree(2,ibranch))      !for remaining particles
+c         smin = totmass**2                        !This affects t_min/max
+c         smax = (m(ibranch) - m(itree(2,ibranch)))**2
+
+         if (ibranch.ne.-ns_channel-1)then
+            pother(:) = P(:,ibranch+1)
+            iopposite = ibranch +1
+         else
+            pother(:) = p(:,2)
+            iopposite = abs(tstrategy)
+         endif
+         s1  = m(ibranch)**2                        !Total mass available
+         ma2 = dot(pother, pother)
+         mb2 = dot(P(0,itree(1,ibranch)),P(0,itree(1,ibranch)))
+         m12 = m(itree(2,ibranch))**2
+         mn2 = m(ibranch-1)**2
+c
+c$$$         write(*,*) itree(1, ibranch), '-----------T----------', itree(2, ibranch), 'm=', m(itree(2, ibranch))
+c$$$         write(*,*)        '                       |           '
+c$$$         write(*,*)   	   '                       |           '
+c$$$         write(*,*)   	   '                       | ',ibranch
+c$$$         write(*,*)   	   '                       |           '
+c$$$         write(*,*)        '                       |           '
+c$$$         if (ibranch.ne.-ns_channel-1)then
+c$$$            write(*,*) iopposite, '-----------T---------- m=', m(ibranch-1)
+c$$$         else
+c$$$            write(*,*) 2, '-----------T---------- m=', m(ibranch-1)
+c$$$            write(*,*) m(1), m(2), m(3), m(4), m(5)
+c$$$         endif
+c$$$         write(*,*) 'Pa', P(0,itree(1, ibranch)),P(1,itree(1, ibranch)),P(2,itree(1, ibranch)),P(3,itree(1, ibranch))
+c$$$         if (ibranch.ne.-ns_channel-1) then
+c$$$            write(*,*) 'Pb', P(0,iopposite),P(1,iopposite),P(2,iopposite), P(3,iopposite)
+c$$$            do i=0,3
+c$$$               pother(i) = P(i,itree(1, ibranch)) + P(i,ibranch+1) 
+c$$$            enddo
+c$$$         else
+c$$$            write(*,*) 'Pb', P(0,2),P(1,2),P(2,2),P(3,2)
+c$$$            do i=0,3
+c$$$               pother(i) = P(i,1) + P(i,2)
+c$$$            enddo
+c$$$         endif
+c$$$         do i=0,3
+c$$$            pother(i) = P(i,itree(1, ibranch)) + P(i,iopposite)
+c$$$         enddo
+c$$$         write(*,*) 'DSQRT(s1) = ', m(ibranch), DSQRT(dot(pother, pother))
+c$$$c         if (m(ibranch)**2.ne.dot(pother, pother)) stop 1
+c$$$         write(*,*) 'm12= Pd**2 = ', m12 ,DSQRT(m12)
+c$$$         write(*,*) 'mn2 = Pc**2 =', mn2, DSQRT(mn2)
+         
+C     WRITE(*,*) 'Enertering yminmax',sqrt(s1),sqrt(m12),sqrt(mn2)
+         
+         call yminmax(s1,0d0,m12,ma2,mb2,mn2,tmin,tmax)
+c         call yminmax(s1,0d0,m12,ma2,mb2,smax,tmin_temp,tmax_temp)
+c         if (tmin_temp.lt.tmin) tmin = tmin_temp
+c         if (tmax_temp.gt.tmax) tmax = tmax_temp
+         
+c
+c     Call for 0<x<1
+c
+c         call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+c     &        .5d0*(tmin/stot+1d0),
+c     &        .5d0*(tmax/stot+1d0))
+c         t   = Stot*(x(-ibranch)*2d0-1d0)
+c
+c     call for -1<x<1
+c
+
+c         write(*,*) 'tmin, tmax/ temp',tmin,tmax, tmin_temp, tmax_temp
+
+c         if (nt_channel.ge.2)then
+c            tmin = max(tmin, -stot)
+c         endif
+c      if ((tmax-tmin)/stot.gt.0.1)then
+c            call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+c     $           0d0, 1d0)
+c         t = stot*(-x(-ibranch))
+      
+c      else if (tmax/stot.gt.-0.01.and.tmin/stot.lt.-0.02)then
+c         set tmax to 0. The idea is to be sure to be able to hit zero
+c         and not to be block by numerical inacuracy
+c         tmax = max(tmax,0d0) !This line if want really t freedom
+c         call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+c     $        0d0, -tmin/stot)
+c         t = stot*(-x(-ibranch))
+
+c      else
+         call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+     $        -tmax/stot, -tmin/stot)
+         t = stot*(-x(-ibranch))
+c      endif
+
+c      call yminmax(s1,0d0,m12,ma2,mb2,mn2,tmin_temp,tmax_temp)
+      if (t .lt. tmin .or. t .gt. tmax) then
+         jac=-3d0
+         return
+      endif
+c
+c     tmin and tmax set to -s,+s for jacobian because part of jacobian
+c     was determined from choosing the point x from the grid based on
+c     tmin and tmax.  (ie wgt contains some of the jacobian)
+c
+         tmin=-stot
+         tmax= stot
+         call sample_get_x(wgt,x(nbranch+(-ibranch-1)*2),
+     &        nbranch+(-ibranch-1)*2,iconfig,0d0,1d0)
+         phi = 2d0*pi*x(nbranch+(-ibranch-1)*2)
+         jac = jac*(tmax-tmin)*2d0*pi /2d0 ! I need /2d0 if -1<x<1
+
+         
+c
+c     Finally generate the momentum. The call is of the form
+c     pa+pb -> p1+ p2; t=(pa-p1)**2;   pr = pa-p1
+c     gentcms(pa,pb,t,phi,m1,m2,p1,pr) 
+c
+         if (itree(1,ibranch).gt.-ns_channel-1)then
+            mi2 = m(itree(1,ibranch))**2
+         else
+            mi2 = tmass(itree(1,ibranch))
+         endif
+         tmass(ibranch) = t
+         call gentcms(p(0,itree(1,ibranch)),p(0,iopposite),t,phi,mi2,
+     &        m(itree(2,ibranch)),m(ibranch-1),p(0,itree(2,ibranch)),
+     &        p(0,ibranch),jac)
+c$$$         write(*,*) 'RESULT'
+c$$$         write(*,*) 'pa', p(0,itree(1,ibranch)),p(1,itree(1,ibranch)),p(2,itree(1,ibranch)),p(3,itree(1,ibranch))
+c$$$         write(*,*) 'pb', p(0,iopposite),p(1,iopposite),p(2,iopposite),p(3,iopposite)
+c$$$         write(*,*) '->'
+c$$$         write(*,*) 'pc', p(0,itree(2,ibranch)),p(1,itree(2,ibranch)),p(2,itree(2,ibranch)),p(3,itree(2,ibranch))
+c$$$         do i =0,3
+c$$$            pother(i) = p(i,itree(1,ibranch)) + p(i,iopposite) - p(i,itree(2,ibranch))
+c$$$         enddo
+c$$$         write(*,*) 'pc', p(0,itree(2,ibranch)),p(1,itree(2,ibranch)),p(2,itree(2,ibranch)),p(3,itree(2,ibranch))
+c$$$         write(*,*) 'pd', pother(0), pother(1), pother(2), pother(3) , DSQRT(dot(pother,pother))
+c$$$         write(*,*) 'T channel'
+c$$$         write(*,*) 'pr', p(0,ibranch),p(1,ibranch),p(2,ibranch),p(3,ibranch)
+c$$$         write(*,*) 'pa-pc', p(0,itree(1,ibranch))-p(0,itree(2,ibranch)),p(1,itree(1,ibranch))-p(1,itree(2,ibranch))
+
+
+         
+         if (jac .lt. 0d0) then
+c            nerr=nerr+1
+c            if(nerr.le.5)
+c     $           write(*,*) 'Failed gentcms',iconfig,ibranch
+            return              !Failed, probably due to negative x
+         endif
+
+         pswgt = pswgt/(4d0*dsqrt(lambda(s1,ma2,mb2)))
+      enddo
+
+c     
+c     We need to get the momentum of the last external particle.
+c     This should just be the sum of p(0,2) and the remaining
+c     momentum from our last t channel 2->2
+c
+      if (nt_channel.eq.1) then
+c$$$         write(*,*) 'need to assign last', itree(2,-nbranch)
+c$$$         write(*,*) 'nbranch is at', nbranch
+c$$$         do i=-nbranch,nexternal
+c$$$            write(*,*) 'p',i, p(0,i),p(1,i),p(2,i),p(3,i)
+c$$$         enddo
+         do i=0,3
+            p(i,itree(2,-nbranch)) = p(i,-nbranch+1)+p(i,2)
+         enddo
+      else
+c$$$                  write(*,*) 'need to assign last', itree(2,-nbranch)
+c$$$         write(*,*) 'nbranch is at', nbranch
+c$$$         do i=-nbranch,nexternal
+c$$$            write(*,*) 'p',i, p(0,i),p(1,i),p(2,i),p(3,i)
+c$$$         enddo
+         do i=0,3
+            p(i,itree(2,-nbranch)) = p(i,-nbranch+1)+p(i,-nbranch+2)
+         enddo
+      endif
+
+
+      else if (tstrategy.eq.2.or.tstrategy.eq.1) then
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+cc       T-channel One side eat all strategy ending with 2
+ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c     Now perform the t-channel decay sequence. Most of this comes from: 
+c     Particle Kinematics Chapter 6 section 3 page 166
+c
+c     From here, on we can just pretend this is a 2->2 scattering with
+c     Pa                    + Pb     -> P1          + P2
 c     p(0,itree(ibranch,1)) + p(0,2) -> p(0,ibranch)+ p(0,itree(ibranch,2))
 c     M(ibranch) is the total mass available (Pa+Pb)^2
 c     M(ibranch-1) is the mass of P2  (all the remaining particles)
 c
       do ibranch=-ns_channel-1,-nbranch+1,-1
          s1  = m(ibranch)**2                        !Total mass available
-         ma2 = m(2)**2
+         ma2 = m(tstrategy)**2
          mb2 = dot(P(0,itree(1,ibranch)),P(0,itree(1,ibranch)))
          m12 = m(itree(2,ibranch))**2
          mn2 = m(ibranch-1)**2
 c         write(*,*) 'Enertering yminmax',sqrt(s1),sqrt(m12),sqrt(mn2)
          call yminmax(s1,0d0,m12,ma2,mb2,mn2,tmin,tmax)
+
+         if(.false.) then
+             write(*,*) itree(1, ibranch), 'a----------T----------', itree(2, ibranch), 'm=', m(itree(2, ibranch))
+             write(*,*)        '                       |           '
+             write(*,*)   	   '                       |           '
+             write(*,*)   	   '                       | ',ibranch
+             write(*,*)   	   '                       |           '
+             write(*,*)        '                       |           '
+             write(*,*) tstrategy, '-----------T---------- m=', m(ibranch-1), ibranch-1
+             write(*,*) 'S', dsqrt(s1), 'm_top=', dsqrt(mb2), 'M_bottom', dsqrt(ma2)
+             
+c     write(*,*) m(1), m(2), m(3), m(4), m(5)
+            write(*,*) 'Pa', P(0,itree(1, ibranch)),P(1,itree(1, ibranch)),P(2,itree(1, ibranch)),P(3,itree(1, ibranch))
+         endif
 c
 c     Call for 0<x<1
 c
@@ -861,15 +1090,80 @@
 c
 
 c         write(*,*) 'tmin, tmax',tmin,tmax
-
-      if (tmax.gt.-0.01.and.tmin.lt.-0.02)then
+         if(.false.) then
+            ! NOT VALIDATED METHOD, momenta are ok but not jacobian
+            
+            call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+     $        0d0, 1d0)
+            costh= 2d0*x(-ibranch)-1d0
+            call sample_get_x(wgt,x(nbranch+(-ibranch-1)*2),
+     &        nbranch+(-ibranch-1)*2,iconfig,0d0,1d0)
+            phi = 2d0*pi*x(nbranch+(-ibranch-1)*2)
+            jac = jac * 4d0*pi
+                     m12 = m(itree(2,ibranch))**2
+         mn2 = m(ibranch-1)**2
+         call mom2cx(dsqrt(s1),m(itree(2,ibranch)),m(ibranch-1),costh,phi,
+     &        p(0,itree(2,ibranch)),pother)
+
+         I= itree(2,ibranch)
+         DO I=0,3
+            pboost(I) = P(I,tstrategy) + P(I,itree(1, ibranch))
+         ENDDO
+         
+         call boostm(p(0,itree(2,ibranch)),pboost,m(itree(2,ibranch)),p(0,itree(2,ibranch)))
+         call boostm(pother,pboost,m(ibranch),pother)
+
+         do I=0,3
+            p(I,ibranch) = pother(i) - p(i, tstrategy)
+         enddo
+         if(.false.)then
+         write(*,*) 'input'
+         write(*,*) 'p(tstrategy=',tstrategy,')', p(0,tstrategy), p(1,tstrategy), p(2,tstrategy), p(3,tstrategy)
+         I=itree(1, ibranch)
+         write(*,*) 'p(',I,',)', p(0,i), p(1,i), p(2,i), p(3,i)
+         write(*,*) 'output'
+         I= itree(2,ibranch)
+         write(*,*) 'p(',I,',)', p(0,i), p(1,i), p(2,i), p(3,i)
+         write(*,*) 'pother', pother(0),pother(1),pother(2),pother(3)
+         write(*,*) 'check'
+         do i=0,3
+            ptemp(i) = p(i,tstrategy) + p(i, itree(1, ibranch))
+         enddo
+         write(*,*) 'pa+pb', ptemp(0), ptemp(1), ptemp(2), ptemp(3), dsqrt(dot(ptemp, ptemp))
+         do i=0,3
+            ptemp(i) = pother(i) + p(i,itree(2,ibranch))
+         enddo
+         write(*,*) 'p1+p2', ptemp(0), ptemp(1), ptemp(2), ptemp(3), dsqrt(dot(ptemp, ptemp))
+         
+         write(*,*) 'Tchannel'
+         I = ibranch
+         write(*,*) 'p(',I,',)', p(0,i), p(1,i), p(2,i), p(3,i), dot(p(0,I), p(0,I))
+         endif
+         
+         pswgt = pswgt/(4d0*dsqrt(lambda(s1,ma2,mb2)))
+
+      else
+
+c     test of impact of low t part
+         if (nt_channel.ge.2)then
+            tmin = max(tmin,  -stot)
+         endif
+      if ((tmax-tmin)/stot.gt.0.1)then
+            call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+     $           0d0, 1d0)
+         t = stot*(-x(-ibranch))
+c     if (dabs(tmax - tmin)/stot.gt.0.05d0) then
+c         call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
+c     $        0d0,  1d0)
+c     set tmax to 0 and tmin to -1 The idea is to avoid dimension correlation
+c     the condition ensure a minimum efficiency in the generation of events
+c         t = stot*(-x(-ibranch))
+      else if (tmax.gt.-0.01.and.tmin.lt.-0.02)then
 c         set tmax to 0. The idea is to be sure to be able to hit zero
 c         and not to be block by numerical inacuracy
-c         tmax = max(tmax,0d0) !This line if want really t freedom
          call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
      $        0d0, -tmin/stot)
          t = stot*(-x(-ibranch))
-
       else
          call sample_get_x(wgt,x(-ibranch),-ibranch,iconfig,
      $        -tmax/stot, -tmin/stot)
@@ -903,7 +1197,7 @@
             mi2 = tmass(itree(1,ibranch))
          endif
          tmass(ibranch) = t
-         call gentcms(p(0,itree(1,ibranch)),p(0,2),t,phi, mi2,
+         call gentcms(p(0,itree(1,ibranch)),p(0,tstrategy),t,phi, mi2,
      &        m(itree(2,ibranch)),m(ibranch-1),p(0,itree(2,ibranch)),
      &        p(0,ibranch),jac)
 
@@ -915,6 +1209,7 @@
          endif
 
          pswgt = pswgt/(4d0*dsqrt(lambda(s1,ma2,mb2)))
+      endif
       enddo
 c
 c     We need to get the momentum of the last external particle.
@@ -922,10 +1217,16 @@
 c     momentum from our last t channel 2->2
 c
       do i=0,3
-         p(i,itree(2,-nbranch)) = p(i,-nbranch+1)+p(i,2)
+         p(i,itree(2,-nbranch)) = p(i,-nbranch+1)+p(i,tstrategy)
       enddo
 
 
+      
+      else
+         write(*,*) 'not supported tstrategy'
+         stop 2
+      endif
+
       endif                     !t-channel stuff
 
 
@@ -960,6 +1261,18 @@
          call boostm(p(0,itree(1,i)),p(0,i),m(i),p(0,itree(1,i)))
          call boostm(p(0,itree(2,i)),p(0,i),m(i),p(0,itree(2,i)))
       enddo
+c$$$      write(*,*) '****'
+c$$$      do i=-nbranch,nexternal
+c$$$         write(*,*) 'mass', i, m(i)
+c$$$      enddo
+c$$$      do i=-nbranch,nexternal
+c$$$         write(*,*) 'p',i, p(0,i),p(1,i),p(2,i),p(3,i)
+c$$$      enddo
+c$$$      do i =0,3
+c$$$         pother(i) = p(i,1) + p(i,2) - p(i,3) -p(i,4)
+c$$$      enddo
+c$$$      write(*,*) 'p5 expected', pother(0), pother(1),pother(2),pother(3)
+
       jac = jac*wgt
       if (.not. pass) jac = -99
       end
@@ -1024,7 +1337,7 @@
 c
 c     Arguments
 c
-      double precision t,phi,m1,m2               !inputs
+      double precision t,phi,m1,m2,ma               !inputs
       double precision pa(0:3),pb(0:3),jac
       double precision p1(0:3),pr(0:3)           !outputs
 c
@@ -1314,6 +1627,146 @@
       return
       end
 
-
+      double precision function get_channel_cut(p, config)
+      implicit none
+
+      include 'maxconfigs.inc'
+      include 'nexternal.inc'
+      include 'genps.inc'
+      include 'maxamps.inc'
+      include 'coupl.inc'
+c     include 'run.inc'
+
+      double precision p(0:3, nexternal)
+      integer config
+      
+
+      integer iforest(2,-max_branch:-1,lmaxconfigs)
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/ iforest, tstrategy
+
+      integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
+      integer tprid(-max_branch:-1,lmaxconfigs)
+      common/to_sprop/sprop,tprid
+
+      double precision stot,m1,m2
+      common/to_stot/stot,m1,m2
+
+      double precision tmin_for_channel
+       integer sde_strat ! 1 means standard single diagram enhancement strategy,
+c      	      	      	   2 means approximation by the	denominator of the propagator
+       common/TO_CHANNEL_STRAT/tmin_for_channel,	sde_strat
+      
+      integer            mapconfig(0:lmaxconfigs), this_config
+      common/to_mconfigs/mapconfig, this_config
+
+      double precision      spole(maxinvar),swidth(maxinvar),bwjac
+      common/to_brietwigner/spole          ,swidth          ,bwjac
+
+      double precision ptemp(0:3, -nexternal:nexternal)
+      integer i,j
+      integer d1, d2
+      double precision t
+      double precision dot
+      external dot
+      integer ns_channel
+      integer nb_tchannel
+      integer nbranch
+      double precision tmp, tmp2
+      
+      double precision ZERO
+      parameter (ZERO=0d0)
+      double precision prmass(-nexternal:0,lmaxconfigs)
+      double precision prwidth(-nexternal:0,lmaxconfigs)
+      integer pow(-nexternal:0,lmaxconfigs)
+      logical first_time
+      save prmass,prwidth,pow
+      data first_time /.true./
+
+      double precision Mass, Width
+      
+      include 'configs.inc'
+
+      if(sde_strat.eq.1.and.tmin_for_channel.eq.-1)then
+         get_channel_cut = 1d0
+         return
+      endif
+      
+      if (first_time) then
+         include 'props.inc'
+         first_time=.false.
+      endif
+      
+      do i = 1, nexternal
+         do j =0,3
+            ptemp(j,i) = p(j,i)
+            ptemp(j,-i) = 0d0
+         enddo
+      enddo
+
+      nbranch = nexternal -2
+      ns_channel=1
+      do while((iforest(1,-ns_channel,config) .ne. 1.and.iforest(1,-ns_channel,config) .ne. 2).and.ns_channel.lt.nbranch)
+         ns_channel=ns_channel+1
+      enddo
+      ns_channel=ns_channel - 1
+      nb_tchannel=nbranch-ns_channel-1
+c      write(*,*) 'T-channel found: ',nb_tchannel
+
+
+
+
+
+      
+      get_channel_cut = 1.
+      if (nb_tchannel.lt.2.and.sde_strat.eq.1)then
+         get_channel_cut = 1.
+         return
+      endif
+      
+      do i = 1, nexternal-3
+         d1 = iforest(1, -i, config)
+         d2 = iforest(2, -i, config)
+         do j=0,3
+            if (d1.gt.0.and.d1.le.2) then
+               ptemp(j,-i) = ptemp(j,-i) - ptemp(j, d1)
+            else
+               ptemp(j,-i) = ptemp(j,-i)+ptemp(j, d1)
+            endif
+            if (d2.gt.0.and.d2.le.2) then
+               ptemp(j,-i) = ptemp(j,-i) - ptemp(j, d2)
+            else
+               ptemp(j,-i) = ptemp(j,-i)+ptemp(j, d2)
+            endif
+         enddo
+         if (tprid(-i,config).ne.0)then
+            if(sde_strat.eq.2)then
+               t = dot(ptemp(0,-i), ptemp(0,-i))
+               Mass  = prmass(-i, config)
+               get_channel_cut = get_channel_cut / ((t-Mass)*(t+Mass))**2
+            endif
+c            write(*,*) i, "t, Mass, fact", t, Mass, ((t-Mass)*(t+Mass))**2,get_channel_cut
+            t = t/stot 
+            if (t.lt.tmin_for_channel)then
+                get_channel_cut = get_channel_cut * exp((t-tmin_for_channel)/(t+1))
+c               get_channel_cut = get_channel_cut * (t+1)/(1+tmin_for_channel)
+c            else if(t.gt.2*tmin_for_channel)then
+c               get_channel_cut = get_channel_cut * (2*tmin_for_channel-t)/tmin_for_channel
+            endif
+         else
+            if(sde_strat.eq.2)then
+               t = dot(ptemp(0,-i), ptemp(0,-i))
+               Mass  = prmass(-i, config)
+               Width = prwidth(-i, config)
+               tmp = (t-Mass)*(t+Mass)
+               tmp2 = Mass*Width
+               get_channel_cut = get_channel_cut* (tmp**2 - tmp2**2)/(tmp**2 + tmp2**2)**2 
+            endif
+c            write(*,*) i, "s, Mass, Width, fact", t, Mass, Width, (((t-Mass)*(t+Mass) )**2 + Width**2*Mass**2), get_channel_cut
+         endif
+      enddo
+c      write(*,*) 'final for config', config, get_channel_cut
+      return
+      end
 
 

=== modified file 'Template/LO/SubProcesses/makefile'
--- Template/LO/SubProcesses/makefile	2017-07-13 10:04:07 +0000
+++ Template/LO/SubProcesses/makefile	2021-01-30 10:13:29 +0000
@@ -30,18 +30,27 @@
 
 # Source files
 
+MATRIX_HEL = $(patsubst %.f,%.o,$(wildcard matrix*_orig.f))
+MATRIX = $(patsubst %.f,%.o,$(wildcard matrix*_optim.f))
+ifeq ($(strip $(MATRIX_HEL)),)
+        MATRIX = $(patsubst %.f,%.o,$(wildcard matrix*.f))
+endif
+
+
 PROCESS= driver.o myamp.o genps.o unwgt.o setcuts.o get_color.o \
          cuts.o cluster.o reweight.o initcluster.o addmothers.o setscales.o \
 	 idenparts.o dummy_fct.o \
          $(patsubst %.f,%.o,$(wildcard auto_dsig*.f)) \
-         $(patsubst %.f,%.o,$(wildcard matrix*.f))
 
 SYMMETRY = symmetry.o idenparts.o 
 
 # Binaries
 
-$(PROG): $(PROCESS) auto_dsig.o $(LIBS)
-	$(FC) -o $(PROG) $(PROCESS) $(LINKLIBS) $(LDFLAGS) $(BIASDEPENDENCIES)
+$(PROG): $(PROCESS) auto_dsig.o $(LIBS) $(MATRIX)
+	$(FC) -o $(PROG) $(PROCESS) $(MATRIX) $(LINKLIBS) $(LDFLAGS) $(BIASDEPENDENCIES)
+
+$(PROG)_forhel: $(PROCESS) auto_dsig.o $(LIBS) $(MATRIX_HEL)
+	$(FC) -o $(PROG)_forhel $(PROCESS) $(MATRIX_HEL) $(LINKLIBS) $(LDFLAGS) $(BIASDEPENDENCIES)
 
 gensym: $(SYMMETRY) configs.inc $(LIBDIR)libmodel.$(libext) $(LIBDIR)libgeneric.$(libext)
 	$(FC) -o gensym $(SYMMETRY) -L../../lib/ -lmodel -lgeneric -lpdf $(llhapdf) $(LDFLAGS)
@@ -78,4 +87,4 @@
 initcluster.o: message.inc
 
 clean:
-	$(RM) *.o gensym madevent
+	$(RM) *.o gensym madevent madevent_forhel

=== modified file 'Template/LO/SubProcesses/myamp.f'
--- Template/LO/SubProcesses/myamp.f	2019-12-04 21:34:25 +0000
+++ Template/LO/SubProcesses/myamp.f	2021-01-30 10:13:29 +0000
@@ -38,7 +38,8 @@
 c
       include 'maxamps.inc'
       integer iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/iforest,	tstrategy
       integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
       integer tprid(-max_branch:-1,lmaxconfigs)
       common/to_sprop/sprop,tprid
@@ -77,7 +78,8 @@
          include 'props.inc'
          nbw = 0
          do i=-1,-(nexternal-3),-1
-            if (iforest(1,i,iconfig) .eq. 1 .or. prwidth(i,iconfig).le.0) then
+            if (iforest(1,i,iconfig) .eq. 1 .or. prwidth(i,iconfig).le.0.or.
+     &                       (nincoming.eq.2.and.iforest(1,i,iconfig) .eq. 2)) then
               cycle
             endif
             nbw=nbw+1
@@ -110,7 +112,7 @@
 c     Start loop over propagators
       do i=-1,-(nexternal-3),-1
          onbw(i) = .false.
-         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
+         if (iforest(1,i,iconfig) .eq. 1.or.(nincoming.eq.2.and.iforest(1,i,iconfig).eq.2)) tsgn=-1d0
          do j=0,3
             xp(j,i) = xp(j,iforest(1,i,iconfig))
      $           +tsgn*xp(j,iforest(2,i,iconfig))
@@ -247,7 +249,8 @@
       common/to_smin/ Smin
 
       integer iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/iforest,	tstrategy
 
       integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
       integer tprid(-max_branch:-1,lmaxconfigs)
@@ -354,8 +357,10 @@
 
 c     Start loop over propagators
       do i=-1,-(nexternal-3),-1
-         if (iforest(1,i,iconfig) .eq. 1) tsgn=-1d0
-         if (tsgn .eq. 1d0) then                         !s channel
+         if (iforest(1,i,iconfig) .eq. 1.or.iforest(1,i,iconfig) .eq. 2)then
+              tsgn=-1d0
+         endif
+         if (tsgn .eq. 1d0) then !s channel
             xm(i) = xm(iforest(1,i,iconfig))+xm(iforest(2,i,iconfig))
             xe(i) = xe(iforest(1,i,iconfig))+xe(iforest(2,i,iconfig))
             mtot = mtot - xm(i)

=== modified file 'Template/LO/SubProcesses/reweight.f'
--- Template/LO/SubProcesses/reweight.f	2020-08-07 15:02:30 +0000
+++ Template/LO/SubProcesses/reweight.f	2021-01-30 10:13:29 +0000
@@ -590,6 +590,8 @@
       common /to_specxpt/xptj,xptb,xpta,xptl,xmtc,xetamin,xqcut,deltaeta
       double precision stot,m1,m2
       common/to_stot/stot,m1,m2
+      logical init_mode
+      common/to_determine_zero_hel/init_mode
 
 C   local variables
       integer i, j, idi, idj, k,m
@@ -597,6 +599,7 @@
       real*8 PI
       parameter( PI = 3.14159265358979323846d0 )
       integer iforest(2,-max_branch:-1,lmaxconfigs)
+      integer tstrategy
       double precision asref, pt2prev(n_max_cl),pt2min
       integer n, ibeam(2), iqcd(0:2)
       integer idfl, idmap(-nexternal:nexternal)
@@ -650,11 +653,12 @@
       endif
  100  clustered = cluster(p(0,1))
       if(.not.clustered) then
+         if(init_mode) goto 999
          open(unit=26,file='../../../error',status='unknown',err=999)
          write(26,*) 'Error: Clustering failed in cluster.f.'
          write(*,*) 'Error: Clustering failed in cluster.f.'
          stop
- 999     write(*,*) 'error'
+ 999     write(*,*) 'error for clustering'
          setclscales=.false.
          clustered = .false.
          return
@@ -1349,6 +1353,7 @@
       logical setclscales
       integer mapconfig(0:lmaxconfigs), this_config
       integer iforest(2,-max_branch:-1,lmaxconfigs)
+      integer tstrategy(lmaxconfigs)
       integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
       integer tprid(-max_branch:-1,lmaxconfigs)
       include 'configs.inc'

=== modified file 'Template/LO/SubProcesses/setcuts.f'
--- Template/LO/SubProcesses/setcuts.f	2019-06-17 14:27:09 +0000
+++ Template/LO/SubProcesses/setcuts.f	2021-01-30 10:13:29 +0000
@@ -879,7 +879,8 @@
       double precision pmass(nexternal)
       common/to_mass/  pmass
       integer iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/iforest,	tstrategy
       integer mapconfig(0:lmaxconfigs), this_config
       common/to_mconfigs/mapconfig, this_config
 
@@ -946,7 +947,8 @@
 c     Global
 c
       integer iforest(2,-max_branch:-1,lmaxconfigs)
-      common/to_forest/ iforest
+      integer tstrategy(lmaxconfigs)
+      common/to_forest/iforest,	tstrategy
       integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
       integer tprid(-max_branch:-1,lmaxconfigs)
       common/to_sprop/sprop,tprid

=== modified file 'Template/LO/SubProcesses/unwgt.f'
--- Template/LO/SubProcesses/unwgt.f	2020-07-30 21:04:54 +0000
+++ Template/LO/SubProcesses/unwgt.f	2021-01-30 10:13:29 +0000
@@ -469,6 +469,8 @@
       include 'run.inc'
       include 'run_config.inc'
 
+
+      
 c
 c     Arguments
 c
@@ -505,6 +507,9 @@
 C     
 C     GLOBAL
 C
+      logical init_mode
+      common/to_determine_zero_hel/init_mode
+      
       double precision twgt, maxwgt,swgt(maxevents)
       integer                             lun, nw, itmin
       common/to_unwgt/twgt, maxwgt, swgt, lun, nw, itmin
@@ -552,7 +557,7 @@
 C-----
       
       if ((nw .ge. maxevents).and.do_write_events) return
-
+      if (init_mode) return
 C     if all the necessary inputs to write the events have already been
 C     computed in the bias module, then directly jump to write_events
       if (AlreadySetInBiasModule) then

=== modified file 'UpdateNotes.txt'
--- UpdateNotes.txt	2021-01-26 14:28:08 +0000
+++ UpdateNotes.txt	2021-01-30 10:13:29 +0000
@@ -1,5 +1,30 @@
 Update notes for MadGraph5_aMC@NLO (in reverse time order)
 
+2.9.0 ()    **** Major speed-up update for LO computation ***
+      Kiran+OM: Optimization of the matrix-element at run-time using recycling helicity method
+                - This fasten LO computation by a factor around 2
+		- This can be turned off via the command "output PATH --hel_recycling=False"
+      OM: Various optimization fot T-channel integration
+          - use different ordering for T-channel. Four different ordering have been implemented
+	    at generation time, the code decides which ordering to use channel-per-channel.
+	    This can be turned off via the command "output PATH --t_strategy=2"
+	  - increase number of maximum iteration for double or more T-channel
+	  - implement an alternative to the standard multi-channel.
+	     instead of using the amplitude square it use the product of the denominator
+	     (the default strategy use depends of the process)
+	  - speed-up for VBF type of process  are often around 100 times faster
+	  - This fixes a lot of issue for processes failing to generate the targetted number of events.
+      OM: Better version of the color-computation (but this leads to only modest gain but for very complex processes.)
+          This optimizations leads to a longer generation time of the code (only for complex processes).
+	  This can be turned off  via "output PATH --jamp_optim=False"
+      OM: New parameter in the run_card:
+          - sde_strategy: allows to change the default multi-channel strategy
+	  - a new phase-space optimization parameter are now easily availabe via the command "update ps_optim" 
+      OM: New global parameter "auto_convert_model"
+          if set on True (set auto_convert_model T or via input/mg5_configuration.txt)
+	  all model crashing due to a python3 compatibility issue of the UFO model will be automatically converted
+	  to a python3 compatible model. 
+
 2.8.3 (26/01/21):
       OM: Buch of bunch fixing related to python3 issue (mainly related to unicode encoding)
       OM: Various fix for reweighting with loop (mainly with python3 as well)

=== modified file 'VERSION'
--- VERSION	2021-01-26 14:28:08 +0000
+++ VERSION	2021-01-30 10:13:29 +0000
@@ -1,4 +1,4 @@
-version = 2.8.3
+version = 2.9.0
 date = 2021-01-26
 
 

=== modified file 'aloha/aloha_parsers.py'
--- aloha/aloha_parsers.py	2020-08-20 21:22:22 +0000
+++ aloha/aloha_parsers.py	2021-01-30 10:13:29 +0000
@@ -198,7 +198,7 @@
 
     aloha_object = ['P', 'PBar', 'PVec','Gamma','Gamma5','Sigma','Mass','PSlash',
                     'OverMass2','Width','Scalar','Spinor','Vector',
-                    'Spin2','Spin32','C','Epsilon','Metric','Identity',
+                    'Spin2','Spin32','C','Epsilon','Metric','Identity', 'IdentityL',
                     'ProjM','ProjP','Coup','Norm', 'EPSL', 'EPST1', 'EPST2', 'PT',
                     'UFP', 'UFM', 'UFPC', 'UFMC',
                     'VFP', 'VFM', 'VFPC', 'VFMC',

=== modified file 'aloha/aloha_writers.py'
--- aloha/aloha_writers.py	2020-08-20 15:33:24 +0000
+++ aloha/aloha_writers.py	2021-01-30 10:13:29 +0000
@@ -221,7 +221,9 @@
             self.declaration.add(('complex',coup))
             
         if self.offshell:
-            if aloha.complex_mass:
+            if 'P1N' in self.tag:
+                pass
+            elif aloha.complex_mass:
                 call_arg.append(('complex','M%s' % self.outgoing))              
                 self.declaration.add(('complex','M%s' % self.outgoing))
             else:
@@ -590,7 +592,7 @@
         """Define the Header of the fortran file. This include
             - momentum conservation
             - definition of the impulsion"""
-                    
+                        
         out = StringIO()
         
         # Define all the required momenta
@@ -612,11 +614,15 @@
                 
             if self.declaration.is_used('P%s' % (i+1)):
                 self.get_one_momenta_def(i+1, out)
-                
+        
         # define the resulting momenta
-        if self.offshell:
-            
-            
+        bypass = False
+        if 'P1N' in self.tag:
+            if  not self.declaration.is_used('P%s' % (self.outgoing)):
+                bypass = True
+
+        if self.offshell and not bypass:
+
             energy_pos = out_size -2
             type = self.particles[self.outgoing-1]
             
@@ -816,23 +822,29 @@
                 coeff = 'denom*'    
                 if not aloha.complex_mass:
                     if self.routine.denominator:
-                        out.write('    denom = %(COUP)s/(%(denom)s)\n' % {'COUP': coup_name,\
+                        if 'P1N' not in self.tag:
+                            out.write('    denom = %(COUP)s/(%(denom)s)\n' % {'COUP': coup_name,\
                                 'denom':self.write_obj(self.routine.denominator)}) 
                     else:
                         out.write('    denom = %(COUP)s/(P%(i)s(0)**2-P%(i)s(1)**2-P%(i)s(2)**2-P%(i)s(3)**2 - M%(i)s * (M%(i)s -CI* W%(i)s))\n' % \
                                   {'i': self.outgoing, 'COUP': coup_name})
                 else:
                     if self.routine.denominator:
-                        raise Exception('modify denominator are not compatible with complex mass scheme')                
-
-                    out.write('    denom = %(COUP)s/(P%(i)s(0)**2-P%(i)s(1)**2-P%(i)s(2)**2-P%(i)s(3)**2 - M%(i)s**2)\n' % \
+                        if 'P1N' not in self.tag:
+                            raise Exception('modify denominator are not compatible with complex mass scheme', self.tag)                
+                    if 'P1N' not in self.tag:
+                        out.write('    denom = %(COUP)s/(P%(i)s(0)**2-P%(i)s(1)**2-P%(i)s(2)**2-P%(i)s(3)**2 - M%(i)s**2)\n' % \
                       {'i': self.outgoing, 'COUP': coup_name})
-                self.declaration.add(('complex','denom'))
+                if 'P1N' not in self.tag:
+                    self.declaration.add(('complex','denom'))
                 if aloha.loop_mode:
                     ptype = 'list_complex'
                 else:
                     ptype = 'list_double'
-                self.declaration.add((ptype,'P%s' % self.outgoing))
+                if 'P1N' not in self.tag:
+                    self.declaration.add((ptype,'P%s' % self.outgoing))
+                else:
+                    coeff = 'COUP*'  
             else:
                 if coup_name == 'COUP':
                     coeff = 'COUP*'

=== modified file 'aloha/create_aloha.py'
--- aloha/create_aloha.py	2020-11-25 20:43:03 +0000
+++ aloha/create_aloha.py	2021-01-30 10:13:29 +0000
@@ -312,7 +312,6 @@
                 propa = [t[1:] for t in self.tag if t.startswith('P')]
                 if propa == ['0']: 
                     if spin == 3 and aloha.unitary_gauge == 2:
-                        misc.sprint(spin)
                         lorentz *= complex(0,1) * self.get_custom_propa('1PS', spin, id)
                         continue
                     else:
@@ -482,6 +481,12 @@
         elif propa == "1PS":
             numerator = "(-1*(P(-1,id)*PBar(-1,id)) * Metric(1, 2) + P(1,id)*PBar(2,id) + PBar(1,id)*P(2,id))"
             denominator = "(P(-3,id)*PBar(-3,id))*P(-2,id)**2"
+        elif propa == "1N":
+            if spin == 3:
+                numerator = '-1*IdentityL(1,2)'
+            else:
+                numerator = "-1"
+            denominator = "1"
         else:
             raise Exception
 
@@ -539,6 +544,7 @@
             return eval(numerator) * propaR
         else:
             return eval(numerator)
+        
     
             
 
@@ -920,7 +926,20 @@
         # self.explicit_combine = False
         request = {}
 
+        #handle special outgoing=-1 flag (which means do it for all particle outgoing)
+        # use for P1N
+        for list_l_name, tag, outgoing in data[:]:
+            if outgoing == -1: #means do for all particle.
+                data.remove((list_l_name, tag, outgoing))
+                l_name = list_l_name[0]
+                lorentz = eval('self.model.lorentz.%s' % l_name)
+                for i in range(len(lorentz.spins)):
+                    data.append((list_l_name, tag, i+1))
+                
+                
+
         for list_l_name, tag, outgoing in data:
+            
             #allow tag to have integer for retro-compatibility
             all_tag = tag[:]
             conjugate = [i for i in tag if isinstance(i, int)]

=== modified file 'aloha/template_files/Makefile_F'
--- aloha/template_files/Makefile_F	2016-03-02 03:03:20 +0000
+++ aloha/template_files/Makefile_F	2021-01-30 10:13:29 +0000
@@ -18,6 +18,8 @@
 
 BASIC_OBJS = aloha_functions.o
 
+FFLAGS += $(ALOHA_FLAF)
+
 include ./aloha_file.inc
 
 all: $(LIBDIR)$(LIBRARY)

=== modified file 'aloha/template_files/aloha_functions.f'
--- aloha/template_files/aloha_functions.f	2013-11-29 07:28:53 +0000
+++ aloha/template_files/aloha_functions.f	2021-01-30 10:13:29 +0000
@@ -2036,4 +2036,37 @@
       return
       end
 
-
+      subroutine CombineAmp(nb, ihels, iwfcts, W1, Wall, Amp)
+
+      integer nb ! size of the vectors
+      integer ihels(*), iwfcts(*)
+      double complex W1(6)
+      double complex Wall(6,*)
+      double complex Amp(*)
+c     local variable
+      integer i
+
+      do i = 1, nb
+         Amp(ihels(i)) =  W1(3) * Wall(3,iwfcts(i)) +
+     &                    W1(4) * Wall(4,iwfcts(i))+
+     &                    W1(5) * Wall(5,iwfcts(i))+
+     &                    W1(6) * Wall(6,iwfcts(i))
+      enddo
+      return
+      end
+      
+      subroutine CombineAmpS(nb, ihels, iwfcts, W1, Wall, Amp)
+
+      integer nb ! size of the vectors
+      integer ihels(*), iwfcts(*)
+      double complex W1(3)
+      double complex Wall(6,*)
+      double complex Amp(*)
+c     local variable
+      integer i
+
+      do i = 1, nb
+         Amp(ihels(i)) =  W1(3) * Wall(3,iwfcts(i))
+      enddo
+      return
+      end

=== modified file 'madgraph/core/diagram_generation.py'
--- madgraph/core/diagram_generation.py	2020-11-04 20:30:11 +0000
+++ madgraph/core/diagram_generation.py	2021-01-30 10:13:29 +0000
@@ -574,7 +574,7 @@
                                      process.get('overall_orders')[key])
             except KeyError:
                 process.get('orders')[key] = process.get('overall_orders')[key]
-
+                
         assert model.get('particles'), \
            "particles are missing in model: %s" %  model.get('particles')
 

=== modified file 'madgraph/core/helas_objects.py'
--- madgraph/core/helas_objects.py	2021-01-18 08:53:05 +0000
+++ madgraph/core/helas_objects.py	2021-01-30 10:13:29 +0000
@@ -954,7 +954,6 @@
             return False
         if self['mothers']:
             nb_t_channel= sum(int(wf.is_t_channel()) for wf in self['mothers'])
-            #raise Exception
         else:
             return True
             
@@ -3935,11 +3934,10 @@
                     # Extract all wavefunctions contributing to this amplitude
                     diagram_wfs = HelasWavefunctionList()
                     for amplitude in diagram.get('amplitudes'):
-                        wavefunctions = \
-                          sorted(HelasWavefunctionList.\
-                               extract_wavefunctions(amplitude.get('mothers')),
-                                 key=lambda wf: wf.get('number'))
-                        for wf in wavefunctions:
+                        wavefunctions = HelasWavefunctionList.\
+                               extract_wavefunctions(amplitude.get('mothers'))
+
+                        for wf in reversed(wavefunctions):
                             # Check if wavefunction already used, otherwise add
                             if wf.get('number') not in wf_numbers and \
                                    wf not in diagram_wfs:

=== modified file 'madgraph/interface/common_run_interface.py'
--- madgraph/interface/common_run_interface.py	2021-01-06 14:57:59 +0000
+++ madgraph/interface/common_run_interface.py	2021-01-30 10:13:29 +0000
@@ -3967,7 +3967,7 @@
         return self.deal_multiple_categories(completion, formatting)
         
 
-    def update_make_opts(self):
+    def update_make_opts(self, run_card=None):
         """update the make_opts file writing the environmental variables
         stored in make_opts_var"""
         make_opts = os.path.join(self.me_dir, 'Source', 'make_opts')
@@ -3982,6 +3982,14 @@
 
         self.make_opts_var['MG5AMC_VERSION'] = misc.get_pkg_info()['version']
 
+        if run_card and run_card.LO:
+            if __debug__ and 'global_flag' not in run_card.user_set:
+                self.make_opts_var['GLOBAL_FLAG'] = "-O -fbounds-check"
+            else:
+                self.make_opts_var['ALOHA_FLAG'] = run_card['global_flag']     
+            self.make_opts_var['ALOHA_FLAG'] = run_card['aloha_flag']     
+            self.make_opts_var['MATRIX_FLAG'] = run_card['matrix_flag']
+
         return self.update_make_opts_full(make_opts, self.make_opts_var)
 
     @staticmethod
@@ -5924,6 +5932,16 @@
                         "As your process seems to be impacted by the issue,\n" +\
                       "You can NOT run with MLM matching/merging. Please check if merging outside MG5aMC are suitable or refrain to use merging with this model") 
                 
+        # 
+        if self.run_card and isinstance(self.run_card,banner_mod.RunCardLO):
+            if not 'sde_strategy' in self.run_card.user_set:
+                if proc_charac['single_color']:
+                    self.run_card['SDE_strategy'] = 2
+                else:
+                    self.run_card['SDE_strategy'] = 1
+                logger.debug("set SDE to %s", self.run_card['SDE_strategy'])
+            else:
+                 logger.debug("keep SDE to %s", self.run_card['SDE_strategy'])
 
         ########################################################################
         #       NLO specific check

=== modified file 'madgraph/interface/madevent_interface.py'
--- madgraph/interface/madevent_interface.py	2020-10-09 22:09:19 +0000
+++ madgraph/interface/madevent_interface.py	2021-01-30 10:13:29 +0000
@@ -3359,7 +3359,23 @@
             self.pass_in_difficult_integration_mode()
         elif self.run_card['hard_survey']:
             self.pass_in_difficult_integration_mode(self.run_card['hard_survey'])
-            
+
+        if self.proc_characteristics['hel_recycling'] and self.run_card['hel_recycling']:
+            jobs, P_zero_result = ajobcreator.get_helicity()
+        else:
+            for p in subproc:
+                for f in misc.glob('matrix*_orig.f', pjoin(self.me_dir, 'SubProcesses', p)):
+                    new_file = f.replace('_orig','_optim')
+                    files.cp(f, f.replace('_orig','_optim'))
+                    f = '%s.o' % f[:-2]
+                    if os.path.exists(f):
+                        files.cp(f, f.replace('_orig','_optim'))
+            try:
+                os.remove(pjoin(self.me_dir, 'SubProcesses', p, 'Hel', 'selection'))
+            except Exception as error:
+                logger.debug(error)
+                pass
+                                
         jobs, P_zero_result = ajobcreator.launch()
         # Check if all or only some fails
         if P_zero_result:
@@ -5692,7 +5708,7 @@
             self.treat_ckkw_matching()
 
         # add the make_opts_var to make_opts
-        self.update_make_opts()
+        self.update_make_opts(self.run_card)
         # reset list of Gdirectory
         self.Gdirs = None
             

=== modified file 'madgraph/interface/madgraph_interface.py'
--- madgraph/interface/madgraph_interface.py	2020-11-08 22:23:40 +0000
+++ madgraph/interface/madgraph_interface.py	2021-01-30 10:13:29 +0000
@@ -300,13 +300,14 @@
     """ The Series of help routine for the MadGraphCmd"""
 
     def help_save(self):
-        logger.info("syntax: save %s FILENAME" % "|".join(self._save_opts),'$MG:color:BLUE')
+        logger.info("syntax: save %s FILENAME [OPTIONS]" % "|".join(self._save_opts),'$MG:color:BLUE')
         logger.info("-- save information as file FILENAME",'$MG:BOLD')
         logger.info("   FILENAME is optional for saving 'options'.")
         logger.info('   By default it uses ./input/mg5_configuration.txt')
         logger.info('   If you put "global" for FILENAME it will use ~/.mg5/mg5_configuration.txt')
         logger.info('   If this files exists, it is uses by all MG5 on the system but continues')
         logger.info('   to read the local options files.')
+        logger.info('   if additional argument are defined for save options, only those arguments will be saved to the configuration file.')
 
     def help_load(self):
         logger.info("syntax: load %s FILENAME" % "|".join(self._save_opts),'$MG:color:BLUE')
@@ -466,6 +467,9 @@
         logger.info("      -nojpeg: no jpeg diagrams will be generated.")
         logger.info("      --noeps=True: no jpeg and eps diagrams will be generated.")
         logger.info("      -name: the postfix of the main file in pythia8 mode.")
+        logger.info("      --jamp_optim=[True|False]: [madevent(default:True)|standalone(default:False)] allows a more efficient code computing the color-factor.")
+        logger.info("      --t_strategy: [madevent] allows to change ordering strategy for t-channel.")
+        logger.info("      --hel_recycling=False: [madevent] forbids helicity recycling optimization")
         logger.info("   Examples:",'$MG:color:GREEN')
         logger.info("       output",'$MG:color:GREEN')
         logger.info("       output standalone MYRUN -f",'$MG:color:GREEN')
@@ -801,7 +805,9 @@
         logger.info(" > T channel propagators. Such channel can sometimes be quite slow to integrate")
         logger.info("zerowidth_tchannel <value>",'$MG:color:GREEN')
         logger.info(" > (default: True) [Used ONLY for tree-level output with madevent]")
-        logger.info(" > set the width to zero for all T-channel propagator --no impact on complex-mass scheme mode")        
+        logger.info(" > set the width to zero for all T-channel propagator --no impact on complex-mass scheme mode")
+        logger.info("auto_convert_model <value>",'$MG:color:GREEN')   
+        logger.info(" > (default: False) If set on True any python2 UFO model will be automatically converted to pyton3 format")     
 #===============================================================================
 # CheckValidForCmd
 #===============================================================================
@@ -2461,7 +2467,8 @@
     @cmd.debug()
     def complete_output(self, text, line, begidx, endidx,
                         possible_options = ['f', 'noclean', 'nojpeg'],
-                        possible_options_full = ['-f', '-noclean', '-nojpeg', '--noeps=True']):
+                        possible_options_full = ['-f', '-noclean', '-nojpeg', '--noeps=True','--hel_recycling=False',
+                                                 '--jamp_optim=', '--t_strategy=']):
         "Complete the output command"
 
         possible_format = self._export_formats
@@ -2470,7 +2477,7 @@
                             'Calculators', 'MadAnalysis', 'SimpleAnalysis',
                             'mg5', 'DECAY', 'EventConverter', 'Models',
                             'ExRootAnalysis', 'HELAS', 'Transfer_Fct', 'aloha',
-                            'matchbox', 'matchbox_cpp', 'tests']
+                            'matchbox', 'matchbox_cpp', 'tests', 'launch']
 
         #name of the run =>proposes old run name
         args = self.split_arg(line[0:begidx])
@@ -2897,7 +2904,8 @@
                     'max_npoint_for_channel',
                     'max_t_for_channel',
                     'zerowidth_tchannel',
-                    'default_unset_couplings']
+                    'default_unset_couplings',
+                    ]
     _valid_nlo_modes = ['all','real','virt','sqrvirt','tree','noborn','LOonly']
     _valid_sqso_types = ['==','<=','=','>']
     _valid_amp_so_types = ['=','<=', '==', '>']
@@ -2948,7 +2956,8 @@
                        'cluster_retry_wait':300,
                        'cluster_size':100,
                        'output_dependencies':'external',
-                       'crash_on_error':False
+                       'crash_on_error':False,
+                       'auto_convert_model': False,
                        }
 
     options_madgraph= {'group_subprocesses': 'Auto',
@@ -3328,10 +3337,12 @@
             raise Exception( 'model to convert need to provide a full path')
         model_dir = args[0]
         
-        answer = self.ask('model conversion to support both py2 and py3 are done in place.\n They are NO guarantee of success.\n It can make the model to stop working under PY2 as well.\n Do you want to proceed?',
-                 'y', ['y','n'])
-        if answer != 'y':
-            return 
+        
+        if not ('-f' not in args or self.options['auto_convert_model']): 
+            answer = self.ask('model conversion to support both py2 and py3 are done in place.\n They are NO guarantee of success.\n It can make the model to stop working under PY2 as well.\n Do you want to proceed?',
+                     'y', ['y','n'])
+            if answer != 'y':
+                return 
         
         #Object_library 
         text = open(pjoin(model_dir, 'object_library.py')).read()
@@ -5493,8 +5504,27 @@
                 else:
                     aloha.aloha_prefix=''
                 
-                self._curr_model = import_ufo.import_model(args[1], prefix=prefix,
+                try:
+                    self._curr_model = import_ufo.import_model(args[1], prefix=prefix,
                         complex_mass_scheme=self.options['complex_mass_scheme'])
+                except ufomodels.UFOError as err:
+                    model_path, _,_ = import_ufo.get_path_restrict(args[1])
+                    if six.PY3 and self.options['auto_convert_model']:
+                        logger.info("fail to load model but auto_convert_model is on True. Trying to convert the model")
+                        
+                        self.exec_cmd('convert model %s' % model_path, errorhandling=False, printcmd=True, precmd=False, postcmd=False)
+                        logger.info('retry the load of the model')
+                        tmp_opt = dict(self.options)
+                        tmp_opt['auto_convert_model'] = False
+                        with misc.TMP_variable(self, 'options', tmp_opt):
+                            try:
+                                self.exec_cmd('import %s' % line, errorhandling=False, printcmd=True, precmd=False, postcmd=False)
+                            except Exception:
+                                raise err
+                    elif six.PY3:
+                        raise self.InvalidCmd('UFO model not python3 compatible. You can convert it via the command \nconvert model %s\nYou can also type \"set auto_convert_model T\" to automatically convert all python2 module to be python3 compatible in the future.' % model_path)
+                    else:
+                        raise
                 if os.path.sep in args[1] and "import" in self.history[-1]:
                     self.history[-1] = 'import model %s' % self._curr_model.get('modelpath+restriction')
 
@@ -6013,7 +6043,6 @@
                 self.options['lhapdf_py2'] = pjoin(prefix,'lhapdf6','bin', 'lhapdf-config')
                 self.exec_cmd('save options %s lhapdf_py2' % config_file)
                 self.options['lhapdf'] = self.options['lhapdf_py2']
-            
         elif tool == 'lhapdf5':
             self.options['lhapdf'] = pjoin(prefix,'lhapdf5','bin', 'lhapdf-config')
             self.exec_cmd('save options %s lhapdf' % config_file, printcmd=False, log=False)            
@@ -6606,8 +6635,8 @@
             """function to apply the patch"""
             text = filetext.read().decode()
             
-            pattern = re.compile(r'''=== renamed directory \'(?P<orig>[^\']*)\' => \'(?P<new>[^\']*)\'''')
-            #=== renamed directory 'Template' => 'Template/LO'
+            pattern = re.compile(r'''^=== renamed directory \'(?P<orig>[^\']*)\' => \'(?P<new>[^\']*)\'''')
+            #= = = renamed directory 'Template' => 'Template/LO'
             for orig, new in pattern.findall(text):
                 shutil.copytree(pjoin(MG5DIR, orig), pjoin(MG5DIR, 'UPDATE_TMP'))
                 full_path = os.path.dirname(pjoin(MG5DIR, new)).split('/')
@@ -6790,8 +6819,8 @@
         else:
             raise self.InvalidCmd('Unknown mode for command install update')
 
-        if not os.path.exists(os.path.join(MG5DIR,'input','.autoupdate')) or \
-                os.path.exists(os.path.join(MG5DIR,'.bzr')):
+        if not os.path.exists(os.path.join(MG5DIR,'input','.autoupdate')):# or \
+                #os.path.exists(os.path.join(MG5DIR,'.bzr')):
             error_text = """This version of MG5 doesn\'t support auto-update. Common reasons are:
             1) This version was loaded via bazaar (use bzr pull to update instead).
             2) This version is a beta release of MG5."""
@@ -6807,7 +6836,7 @@
             return
 
         # read the data present in .autoupdate
-        data = {}
+        data = {'last_message':0}
         for line in open(os.path.join(MG5DIR,'input','.autoupdate')):
             if not line.strip():
                 continue
@@ -6837,9 +6866,16 @@
         signal.alarm(timeout)
         to_update = 0
         try:
-            filetext = six.moves.urllib.request.urlopen('http://madgraph.phys.ucl.ac.be/mg5amc_build_nb')
+            filetext = six.moves.urllib.request.urlopen('http://madgraph.physics.illinois.edu/mg5amc_build_nb')
             signal.alarm(0)
-            web_version = int(filetext.read().strip())
+            text = filetext.read().decode().split('\n')
+            web_version = int(text[0].strip())
+            try:
+                msg_version = int(text[1].strip())
+                message = '\n'.join(text[2:])
+            except:
+                msg_version = 0
+                message = ""
         except (TimeOutError, ValueError, IOError):
             signal.alarm(0)
             print('failed to connect server')
@@ -6847,17 +6883,35 @@
                 # wait 24h before next check
                 fsock = open(os.path.join(MG5DIR,'input','.autoupdate'),'w')
                 fsock.write("version_nb   %s\n" % data['version_nb'])
-                fsock.write("last_check   %s\n" % \
-                int(time.time()) - 3600 * 24 * (self.options['auto_update'] -1))
+                fsock.write("last_check   %s\n" % (\
+                int(time.time()) - 3600 * 24 * (self.options['auto_update'] -1)))
+                fsock.write("last_message   %s\n" % data['last_message'])
                 fsock.close()
             return
 
+        if msg_version > data['last_message']:
+            data['last_message'] = msg_version
+            logger.info("************* INFORMATION *************", '$MG:BOLD')
+            logger.info(message.replace('\n','\n    '))
+            logger.info("************* INFORMATION *************", '$MG:BOLD')
+            fsock = open(os.path.join(MG5DIR,'input','.autoupdate'),'w')
+            fsock.write("version_nb   %s\n" % data['version_nb'])
+            fsock.write("last_check   %s\n" % (\
+            int(time.time()) - 3600 * 24 * (int(self.options['auto_update']) -1)))
+            fsock.write("last_message   %s\n" % data['last_message'])
+            fsock.close()
+            
+        if os.path.exists(os.path.join(MG5DIR,'.bzr')):
+            logger.info("bzr version: use bzr pull to update")
+            return 
+        
         if web_version == data['version_nb']:
             logger.info('No new version of MG5 available')
             # update .autoupdate to prevent a too close check
             fsock = open(os.path.join(MG5DIR,'input','.autoupdate'),'w')
             fsock.write("version_nb   %s\n" % data['version_nb'])
             fsock.write("last_check   %s\n" % int(time.time()))
+            fsock.write("last_message   %s\n" % data['last_message'])
             fsock.close()
             return
         elif data['version_nb'] > web_version:
@@ -6865,6 +6919,7 @@
             fsock = open(os.path.join(MG5DIR,'input','.autoupdate'),'w')
             fsock.write("version_nb   %s\n" % data['version_nb'])
             fsock.write("last_check   %s\n" % int(time.time()))
+            fsock.write("last_message   %s\n" % data['last_message'])
             fsock.close()
             return
         else:
@@ -7760,9 +7815,9 @@
             else:
                 raise self.InvalidCmd('expected bool for notification_center')
         # True/False formatting
-        elif args[0] in ['crash_on_error']:
+        elif args[0] in ['crash_on_error', 'auto_convert_model']:
             try:
-                tmp = banner_module.ConfigFile.format_variable(args[1], bool, 'crash_on_error')
+                tmp = banner_module.ConfigFile.format_variable(args[1], bool, args[0])
             except Exception:
                 if args[1].lower() in ['never']:
                     tmp = args[1].lower()
@@ -7833,8 +7888,12 @@
                     
         if '--postpone_model' in args:
             flaglist.append('store_model')
-        
-        line_options = dict(arg[2:].split('=') for arg in args if arg.startswith('--') and '=' in arg)
+        if '--hel_recycling=False' in args:
+            flaglist.append('no_helrecycling')
+                    
+        line_options = dict( (arg[2:].split('=')  if "=" in arg else (arg[2:], True))
+                             for arg in args if arg.startswith('--'))
+#        line_options = dict(arg[2:].split('=') for arg in args if arg.startswith('--') and '=' not in arg)
         main_file_name = ""
         try:
             main_file_name = args[args.index('-name') + 1]
@@ -7960,7 +8019,7 @@
 ||   -> add process <proc_def>
 """)
                     group_processes = False
-
+    
         #Exporter + Template
         if options['exporter'] == 'v4':
             self._curr_exporter = export_v4.ExportV4Factory(self, noclean, 
@@ -8019,7 +8078,7 @@
             if self._curr_amps and self._curr_amps[0].get_ninitial() == 1:
                 options['zerowidth_tchannel'] = False
             
-            self._curr_helas_model = helas_call_writers.FortranUFOHelasCallWriter(self._curr_model)
+            self._curr_helas_model = helas_call_writers.FortranUFOHelasCallWriter(self._curr_model, options=options)
 
         version = [arg[10:] for arg in args if arg.startswith('--version=')]
         if version:
@@ -8272,6 +8331,13 @@
             # these processes
             wanted_lorentz = self._curr_matrix_elements.get_used_lorentz()
             wanted_couplings = self._curr_matrix_elements.get_used_couplings()
+
+            if self._export_format == 'madevent' and not 'no_helrecycling' in flaglist:
+                for (name, flag, out) in wanted_lorentz[:]:
+                    if out == 0:
+                        newflag = list(flag) + ['P1N']
+                        wanted_lorentz.append((name, tuple(newflag), -1))
+                
             # For a unique output of multiple type of exporter need to store this
             # information.             
             if hasattr(self, 'previous_lorentz'):

=== modified file 'madgraph/interface/reweight_interface.py'
--- madgraph/interface/reweight_interface.py	2021-01-21 13:08:04 +0000
+++ madgraph/interface/reweight_interface.py	2021-01-30 10:13:29 +0000
@@ -1293,7 +1293,7 @@
                 scale2 = event.scale**2
             else:
                 scale2 = 0
-                
+
         with misc.chdir(Pdir):
             with misc.stdchannel_redirected(sys.stdout, os.devnull):
                 me_value = module.smatrixhel(pdg, pid, p, event.aqcd, scale2, nhel)

=== modified file 'madgraph/iolibs/export_fks.py'
--- madgraph/iolibs/export_fks.py	2020-09-24 08:52:15 +0000
+++ madgraph/iolibs/export_fks.py	2021-01-30 10:13:29 +0000
@@ -2016,8 +2016,10 @@
         replace_dict['amp2_lines'] = '\n'.join(amp2_lines)
     
         # Extract JAMP lines
-        jamp_lines = self.get_JAMP_lines(matrix_element)
+        jamp_lines, nb_tmp_jamp = self.get_JAMP_lines(matrix_element)
         replace_dict['jamp_lines'] = '\n'.join(jamp_lines)
+        replace_dict['nb_temp_jamp'] = nb_tmp_jamp
+
 
         # Set the size of Wavefunction
         if not self.model or any([p.get('spin') in [4,5] for p in self.model.get('particles') if p]):
@@ -2109,8 +2111,9 @@
         replace_dict['amp2_lines'] = '\n'.join(amp2_lines)
     
         # Extract JAMP lines
-        jamp_lines = self.get_JAMP_lines(matrix_element)
+        jamp_lines, nb_tmp_jamp = self.get_JAMP_lines(matrix_element)
         replace_dict['jamp_lines'] = '\n'.join(jamp_lines)
+        replace_dict['nb_temp_jamp'] = nb_tmp_jamp
 
         # Extract den_factor_lines
         den_factor_lines = self.get_den_factor_lines(fksborn)
@@ -2272,20 +2275,15 @@
         replace_dict['amp2_lines'] = '\n'.join(amp2_lines)
     
         # Extract JAMP lines
-        jamp_lines = self.get_JAMP_lines(matrix_element)
-        new_jamp_lines = []
-        for line in jamp_lines:
-            line = line.replace('JAMP', 'JAMP1')
-            new_jamp_lines.append(line)
-        replace_dict['jamp1_lines'] = '\n'.join(new_jamp_lines)
+        jamp_lines, nb_tmp_jamp = self.get_JAMP_lines(matrix_element, JAMP_format="JAMP1(%s)")
+        replace_dict['jamp1_lines'] = '\n'.join(jamp_lines)
+        replace_dict['nb_temp_jamp'] = nb_tmp_jamp
     
+        
         matrix_element.set('color_basis', link['link_basis'] )
-        jamp_lines = self.get_JAMP_lines(matrix_element)
-        new_jamp_lines = []
-        for line in jamp_lines:
-            line = line.replace( 'JAMP', 'JAMP2')
-            new_jamp_lines.append(line)
-        replace_dict['jamp2_lines'] = '\n'.join(new_jamp_lines)
+        jamp_lines, nb_tmp_jamp = self.get_JAMP_lines(matrix_element, JAMP_format="JAMP2(%s)")
+        replace_dict['jamp2_lines'] = '\n'.join(jamp_lines)
+        replace_dict['nb_temp_jamp'] = max(nb_tmp_jamp, replace_dict['nb_temp_jamp'])
     
     
         # Extract the number of FKS process
@@ -2526,9 +2524,9 @@
             replace_dict['wavefunctionsize'] = 8
     
         # Extract JAMP lines
-        jamp_lines = self.get_JAMP_lines(matrix_element)
-    
+        jamp_lines, nb_tmp_jamp = self.get_JAMP_lines(matrix_element)
         replace_dict['jamp_lines'] = '\n'.join(jamp_lines)
+        replace_dict['nb_temp_jamp'] = nb_tmp_jamp
     
         realfile = open(os.path.join(_file_path, \
                              'iolibs/template_files/realmatrix_fks.inc')).read()
@@ -3389,6 +3387,7 @@
     """Class to take care of exporting a set of matrix elements to
     Fortran (v4) format."""
 
+    jamp_optim = True 
 
     def finalize(self, *args, **opts):
         ProcessExporterFortranFKS.finalize(self, *args, **opts)

=== modified file 'madgraph/iolibs/export_v4.py'
--- madgraph/iolibs/export_v4.py	2021-01-18 22:00:26 +0000
+++ madgraph/iolibs/export_v4.py	2021-01-30 10:13:29 +0000
@@ -12,10 +12,11 @@
 # For more information, visit madgraph.phys.ucl.ac.be and amcatnlo.web.cern.ch
 #
 ################################################################################
-from __future__ import absolute_import
+from __future__ import absolute_import, division
 from madgraph.iolibs.helas_call_writers import HelasCallWriter
 from six.moves import range
 from six.moves import zip
+from fractions import Fraction
 """Methods and classes to export matrix elements to v4 format."""
 
 import copy
@@ -33,6 +34,7 @@
 import sys
 import time
 import traceback
+import  collections
 
 import aloha
 
@@ -164,6 +166,7 @@
                         'output_options':{}
                         }
     grouped_mode = False
+    jamp_optim = False
 
     def __init__(self,  dir_path = "", opt=None):
         """Initiate the ProcessExporterFortran with directory information"""
@@ -174,7 +177,6 @@
         self.opt = dict(self.default_opt)
         if opt:
             self.opt.update(opt)
-        
         self.cmd_options = self.opt['output_options']
         
         #place holder to pass information to the run_interface
@@ -1072,7 +1074,7 @@
                 enumerate(matrix_element.get('color_matrix').\
                                                  get_line_denominators()):
                 # First write the common denominator for this color matrix line
-                ret_list.append("DATA Denom(%i)/%i/" % (index + 1, denominator))
+                #ret_list.append("DATA Denom(%i)/%i/" % (index + 1, denominator))
                 # Then write the numerators for the matrix elements
                 num_list = matrix_element.get('color_matrix').\
                                             get_line_numerators(index, denominator)
@@ -1082,10 +1084,10 @@
                 for k in range(0, len(num_list), n):
                     ret_list.append("DATA (CF(i,%3r),i=%3r,%3r) /%s/" % \
                                     (index + 1, k + 1, min(k + n, len(num_list)),
-                                     ','.join(["%5i" % int(i) for i in num_list[k:k + n]])))
+                                     ','.join([("%.15e" % (int(i)/denominator)).replace('e','d') for i in num_list[k:k + n]])))
+                
                 my_cs.from_immutable(sorted(matrix_element.get('color_basis').keys())[index])
                 ret_list.append("C %s" % repr(my_cs))
-
             return ret_list
 
 
@@ -1113,6 +1115,7 @@
                              ','.join([".true." for i in range(configs)])))
             return ret_list
 
+
         # There is a color basis - create a list showing which JAMPs have
         # contributions to which configs
 
@@ -1138,6 +1141,8 @@
                     # Add this JAMP number to this diag_num
                     diag_jamp[diag_num] = diag_jamp.setdefault(diag_num, []) + \
                                           [ijamp+1]
+                else:
+                    self.proc_characteristic['single_color'] = False
 
         colamps = ijamp + 1
         for iconfig, num_diag in enumerate(mapconfigs):        
@@ -1154,7 +1159,7 @@
 
         return ret_list
 
-    def get_amp2_lines(self, matrix_element, config_map = []):
+    def get_amp2_lines(self, matrix_element, config_map = [], replace_dict=None):
         """Return the amp2(i) = sum(amp for diag(i))^2 lines"""
 
         nexternal, ninitial = matrix_element.get_nexternal_ninitial()
@@ -1195,9 +1200,12 @@
                 # Not using \sum |M|^2 anymore since this creates troubles
                 # when ckm is not diagonal due to the JIM mechanism.
                 if '+' in amp:
-                    line += "(%s)*dconjg(%s)" % (amp, amp)
+                    amp = "(%s)*dconjg(%s)" % (amp, amp)
                 else:
-                    line += "%s*dconjg(%s)" % (amp, amp)
+                    amp = "%s*dconjg(%s)" % (amp, amp)
+                
+                line =  line + "%s" % (amp)
+                #line += " * get_channel_cut(p, %s) " % (config)
                 ret_lines.append(line)
         else:
             for idiag, diag in enumerate(matrix_element.get('diagrams')):
@@ -1321,6 +1329,7 @@
         # construct the list of color_amplitudes for JAMP to be constructed
         # accordingly.
         res_list=[]
+        max_tmp = 0
         for i, amp_order in enumerate(split_order_amps):
             col_amps_order = []
             for jamp in color_amplitudes:
@@ -1332,12 +1341,14 @@
             if self.opt['export_format'] in ['madloop_matchbox']:
                 res_list.extend(self.get_JAMP_lines(col_amps_order,
                                    JAMP_format="JAMP(%s,{0})".format(str(i+1)),
-                                   JAMP_formatLC="LNJAMP(%s,{0})".format(str(i+1))))
+                                   JAMP_formatLC="LNJAMP(%s,{0})".format(str(i+1)))[0])
             else:
-                res_list.extend(self.get_JAMP_lines(col_amps_order,
-                                   JAMP_format="JAMP(%s,{0})".format(str(i+1))))         
+                toadd, nb_tmp = self.get_JAMP_lines(col_amps_order,
+                                   JAMP_format="JAMP(%s,{0})".format(str(i+1)))
+                res_list.extend(toadd)
+                max_tmp = max(max_tmp, nb_tmp)         
 
-        return res_list
+        return res_list, max_tmp
 
 
     def get_JAMP_lines(self, col_amps, JAMP_format="JAMP(%s)", AMP_format="AMP(%s)", 
@@ -1361,7 +1372,7 @@
         else:
             raise MadGraph5Error("Incorrect col_amps argument passed to get_JAMP_lines")
 
-
+        all_element = {}
         res_list = []
         for i, coeff_list in enumerate(color_amplitudes):
             # It might happen that coeff_list is empty if this function was
@@ -1396,6 +1407,11 @@
                 for (coefficient, amp_number) in coefs:
                     if not coefficient:
                         continue
+                    value = (1j if coefficient[2] else 1)* coefficient[0] * coefficient[1] * fractions.Fraction(3)**coefficient[3]
+                    if (i+1, amp_number) not in all_element:
+                        all_element[(i+1, amp_number)] = value
+                    else:
+                        all_element[(i+1, amp_number)] += value
                     if common_factor:
                         res = (res + "%s" + AMP_format) % \
                                                    (self.coeff(coefficient[0],
@@ -1414,8 +1430,151 @@
                     res = res + ')'
     
                 res_list.append(res)
+        
+        if 'jamp_optim' in self.cmd_options:
+            jamp_optim = banner_mod.ConfigFile.format_variable(self.cmd_options['jamp_optim'], bool, 'jamp_optim')
+        else:
+            # class default
+            jamp_optim = self.jamp_optim
+                
+        if not jamp_optim:
+            return res_list, 0
+        else:
+            saved = list(res_list)
+        
+        if len(all_element) > 1000:
+            logger.info("Computing Color-Flow optimization [%s term]", len(all_element))
+            start_time = time.time()
+        else: 
+            start_time = 0
+        
+        res_list = []
+        #misc.sprint(len(all_element))  
+        
+        self.myjamp_count = 0
+        new_mat, defs = self.optimise_jamp(all_element)
+        if start_time:
+            logger.info("Color-Flow passed to %s term in %ss. Introduce %i contraction", len(new_mat), int(time.time()-start_time), len(defs))
+        
+        
+        #misc.sprint("number of iteration", self.myjamp_count)
+        def format(frac):
+            if isinstance(frac, Fraction):
+                if frac.denominator == 1:
+                    return str(frac.numerator)
+                else:
+                    return "%id0/%id0" % (frac.numerator, frac.denominator)
+            elif frac.real == frac:
+                #misc.sprint(frac.real, frac)
+                return str(float(frac.real)).replace('e','d')
+            else:
+                return str(frac).replace('e','d').replace('j','*imag1')
+                
+        
+        
+        for i, amp1, amp2, frac, nb in defs:
+            if amp1 > 0:
+                amp1 = AMP_format % amp1
+            else:
+                amp1 = "TMP_JAMP(%d)" % -amp1
+            if amp2 > 0:
+                amp2 = AMP_format % amp2
+            else:
+                amp2 = "TMP_JAMP(%d)" % -amp2
+                
+            res_list.append(' TMP_JAMP(%d) = %s + (%s) * %s ! used %d times' % (i,amp1, format(frac), amp2, nb))                
+                 
+
+#        misc.sprint(new_mat)
+        jamp_res = collections.defaultdict(list)
+        max_jamp=0
+        for (jamp, var), factor in new_mat.items():
+            if var > 0:
+                name = AMP_format % var
+            else:
+                name = "TMP_JAMP(%d)" % -var
+            jamp_res[jamp].append("(%s)*%s" % (format(factor), name))
+            max_jamp = max(max_jamp, jamp)
+        
+        
+        for i in range(1,max_jamp+1):
+            name = JAMP_format % i
+            res_list.append(" %s = %s" %(name, '+'.join(jamp_res[i])))
+            
+        return res_list, len(defs)
+
+    def optimise_jamp(self, all_element, nb_line=0, nb_col=0, added=0):
+        """ optimise problem of type Y = A X
+                A is a matrix (all_element)
+                X is the fortran name of the input.
+            The code iteratively add sub-expression jtemp[sub_add]
+            and recall itself (this is add to the X size)
+        """
+        self.myjamp_count +=1
+        
+        if not nb_line:
+            for i,j in all_element:
+                if i > nb_line:
+                    nb_line = i+1
+                if j> nb_col:
+                    nb_col = j+1
+        
+        #misc.sprint(nb_line, nb_col)
+        
+
+        max_count = 0
+        all_index = []
+        operation = collections.defaultdict(lambda: collections.defaultdict(int))
+        for i in range(nb_line):
+            for j1 in range(-added, nb_col):
+                v1 = all_element.get((i,j1), 0)
+                if not v1: 
+                    continue                    
+                for j2 in range(j1+1, nb_col):
+                    R = all_element.get((i,j2), 0)/v1
+                    if not R:
+                        continue
                     
-        return res_list
+                    #misc.sprint(j1,j2)
+                    operation[(j1,j2)][R] +=1 
+                    if operation[(j1,j2)][R] > max_count:
+                        max_count = operation[(j1,j2)][R]
+                        all_index = [(j1,j2, R)]
+                    elif operation[(j1,j2)][R] == max_count:
+                        all_index.append((j1,j2, R))
+        if max_count <= 1:
+            return all_element, []
+        #added += 1
+        #misc.sprint(max_count, len(all_index))
+        #misc.sprint(operation)
+        to_add = []
+        for index in all_index:
+            j1,j2,R = index
+            first = True
+            for i in range(nb_line):
+                v1 = all_element.get((i,j1), 0)
+                v2 = all_element.get((i,j2), 0)
+                if not v1 or not v2: 
+                    continue
+                if v2/v1 == R:
+                    if first:
+                        first = False
+                        added +=1
+                        to_add.append((added,j1,j2,R, max_count))
+                        
+                    all_element[(i,-added)] = v1
+                    del all_element[(i,j1)] #= 0
+                    del all_element[(i,j2)] #= 0 
+
+        logger.log(5,"Define %d new shortcut reused %d times", len(to_add), max_count)
+        new_element, new_def =  self.optimise_jamp(all_element, nb_line=nb_line, nb_col=nb_col, added=added)
+        for one_def in to_add:
+            new_def.insert(0, one_def)
+        return new_element, new_def   
+           
+           
+            
+            
 
     def get_pdf_lines(self, matrix_element, ninitial, subproc_group = False):
         """Generate the PDF lines for the auto_dsig.f file"""
@@ -1880,8 +2039,9 @@
             if not version:# not linux 
                 version = 14 # set version to remove MACFLAG
             else:
-                version = int(version.split('.')[1])
-            if version >= 14:
+                majversion, version = [int(x) for x in version.split('.',3)[:2]]
+
+            if majversion >= 11 or (majversion ==10 and version >= 14):
                 for_update['MACFLAG'] = '-mmacosx-version-min=10.8' if is_lc else ''
 
         if not root_dir:
@@ -2536,7 +2696,7 @@
         if len(split_orders)==0:
             replace_dict['nSplitOrders']=''
             # Extract JAMP lines
-            jamp_lines = self.get_JAMP_lines(matrix_element)
+            jamp_lines, nb_tmp_jamp = self.get_JAMP_lines(matrix_element)
             # Consider the output of a dummy order 'ALL_ORDERS' for which we
             # set all amplitude order to weight 1 and only one squared order
             # contribution which is of course ALL_ORDERS=2.
@@ -2545,6 +2705,8 @@
             replace_dict['chosen_so_configs'] = '.TRUE.'
             replace_dict['nSqAmpSplitOrders']=1
             replace_dict['split_order_str_list']=''
+            replace_dict['nb_temp_jamp'] = nb_tmp_jamp
+
         else:
             squared_orders, amp_orders = matrix_element.get_split_orders_mapping()
             replace_dict['nAmpSplitOrders']=len(amp_orders)
@@ -2556,9 +2718,9 @@
             sqamp_so = self.get_split_orders_lines(squared_orders,'SQSPLITORDERS')
             replace_dict['ampsplitorders']='\n'.join(amp_so)
             replace_dict['sqsplitorders']='\n'.join(sqamp_so)           
-            jamp_lines = self.get_JAMP_lines_split_order(\
+            jamp_lines, nb_tmp_jamp = self.get_JAMP_lines_split_order(\
                        matrix_element,amp_orders,split_order_names=split_orders)
-            
+            replace_dict['nb_temp_jamp'] = nb_tmp_jamp
             # Now setup the array specifying what squared split order is chosen
             replace_dict['chosen_so_configs']=self.set_chosen_SO_index(
                               matrix_element.get('processes')[0],squared_orders)
@@ -2733,7 +2895,7 @@
         else:
             raise MadGraph5Error(error_msg % 'col_amps')
 
-        text = super(ProcessExporterFortranMatchBox, self).get_JAMP_lines(col_amps,
+        text, nb = super(ProcessExporterFortranMatchBox, self).get_JAMP_lines(col_amps,
                                             JAMP_format=JAMP_format,
                                             AMP_format=AMP_format,
                                             split=-1)
@@ -2749,12 +2911,13 @@
                     to_add.append( (coefficient, amp_number) )
             LC_col_amps.append(to_add)
            
-        text += super(ProcessExporterFortranMatchBox, self).get_JAMP_lines(LC_col_amps,
+        text2, nb = super(ProcessExporterFortranMatchBox, self).get_JAMP_lines(LC_col_amps,
                                             JAMP_format=JAMP_formatLC,
                                             AMP_format=AMP_format,
                                             split=-1)
+        text += text2 
         
-        return text
+        return text, 0
 
 
 
@@ -2767,6 +2930,7 @@
     MadGraph v4 - MadWeight format."""
 
     matrix_file="matrix_standalone_v4.inc"
+    jamp_optim = False
 
     def copy_template(self, model):
         """Additional actions needed for setup of Template
@@ -3238,7 +3402,7 @@
         replace_dict['helas_calls'] = "\n".join(helas_calls)
 
         # Extract JAMP lines
-        jamp_lines = self.get_JAMP_lines(matrix_element)
+        jamp_lines, nb = self.get_JAMP_lines(matrix_element)
         replace_dict['jamp_lines'] = '\n'.join(jamp_lines)
         
         replace_dict['template_file'] =  os.path.join(_file_path, \
@@ -3543,6 +3707,29 @@
     matrix_file = "matrix_madevent_v4.inc"
     done_warning_tchannel = False
     
+    default_opt = {'clean': False, 'complex_mass':False,
+                        'export_format':'madevent', 'mp': False,
+                        'v5_model': True,
+                        'output_options':{},
+                        'hel_recycling': False
+                        }
+    jamp_optim = True
+    
+    def __init__(self,  dir_path = "", opt=None):
+        
+        super(ProcessExporterFortranME, self).__init__(dir_path, opt)
+        
+        # check and format the hel_recycling options as it should if provided 
+        if opt and isinstance(opt['output_options'], dict) and \
+                                       'hel_recycling' in opt['output_options']:
+            self.opt['hel_recycling'] = banner_mod.ConfigFile.format_variable(
+                  opt['output_options']['hel_recycling'], bool, 'hel_recycling')
+
+        if opt and isinstance(opt['output_options'], dict) and \
+                                       't_strategy' in opt['output_options']:
+            self.opt['t_strategy'] = banner_mod.ConfigFile.format_variable(
+                  opt['output_options']['t_strategy'], int, 't_strategy')
+
     # helper function for customise helas writter
     @staticmethod
     def custom_helas_call(call, arg):
@@ -3708,8 +3895,6 @@
         if not self.model:
             self.model = matrix_element.get('processes')[0].get('model')
 
-
-
         #os.chdir(path)
         # Create the directory PN_xx_xxxxx in the specified path
         subprocdir = "P%s" % matrix_element.get('processes')[0].shell_string()
@@ -3736,8 +3921,12 @@
         self.write_driver(writers.FortranWriter(filename),ncomb,n_grouped_proc=1,
                           v5=self.opt['v5_model'])
 
+
         # Create the matrix.f file, auto_dsig.f file and all inc files
-        filename = pjoin(Ppath, 'matrix.f')
+        if self.opt['hel_recycling']:
+            filename = pjoin(Ppath, 'matrix_orig.f')
+        else:
+            filename = pjoin(Ppath, 'matrix.f')
         calls, ncolor = \
                self.write_matrix_element_v4(writers.FortranWriter(filename),
                       matrix_element, fortran_model, subproc_number = me_number)
@@ -4084,7 +4273,8 @@
         # The proc prefix is not used for MadEvent output so it can safely be set
         # to an empty string.
         replace_dict = {'proc_prefix':''}
-
+ 
+ 
         # Extract helas calls
         helas_calls = fortran_model.get_matrix_element_calls(\
                     matrix_element)
@@ -4189,7 +4379,7 @@
             replace_dict['wavefunctionsize'] = 6
 
         # Extract amp2 lines
-        amp2_lines = self.get_amp2_lines(matrix_element, config_map)
+        amp2_lines = self.get_amp2_lines(matrix_element, config_map, replace_dict)
         replace_dict['amp2_lines'] = '\n'.join(amp2_lines)
 
         # The JAMP definition depends on the splitting order
@@ -4198,6 +4388,8 @@
             squared_orders, amp_orders = matrix_element.get_split_orders_mapping()
             replace_dict['chosen_so_configs']=self.set_chosen_SO_index(
                               matrix_element.get('processes')[0],squared_orders)
+            replace_dict['select_configs_if'] = '          IF (CHOSEN_SO_CONFIGS(SQSOINDEX%(proc_id)s(M,N))) THEN' % replace_dict
+            replace_dict['select_configs_endif'] = ' endif'
         else:
             # Consider the output of a dummy order 'ALL_ORDERS' for which we
             # set all amplitude order to weight 1 and only one squared order
@@ -4205,6 +4397,9 @@
             squared_orders = [(2,),]
             amp_orders = [((1,),tuple(range(1,ngraphs+1)))]
             replace_dict['chosen_so_configs'] = '.TRUE.'
+            # addtionally set the function to NOT be called
+            replace_dict['select_configs_if'] = ''
+            replace_dict['select_configs_endif'] = ''
             
         replace_dict['nAmpSplitOrders']=len(amp_orders)
         replace_dict['nSqAmpSplitOrders']=len(squared_orders)
@@ -4219,10 +4414,11 @@
 
         # Extract JAMP lines
         # If no split_orders then artificiall add one entry called 'ALL_ORDERS'
-        jamp_lines = self.get_JAMP_lines_split_order(\
+        jamp_lines, nb_temp = self.get_JAMP_lines_split_order(\
                              matrix_element,amp_orders,split_order_names=
                         split_orders if len(split_orders)>0 else ['ALL_ORDERS'])
         replace_dict['jamp_lines'] = '\n'.join(jamp_lines)
+        replace_dict['nb_temp_jamp'] = nb_temp
 
         replace_dict['template_file'] = pjoin(_file_path, \
                           'iolibs/template_files/%s' % self.matrix_file)
@@ -4595,7 +4791,7 @@
         For s-channels, we need to output one PDG for each subprocess in
         the subprocess group, in order to be able to pick the right
         one for multiprocesses."""
-
+        
         lines = []
 
         s_and_t_channels = []
@@ -4635,12 +4831,18 @@
                 else:
                     stchannels.append((empty_verts, None))
 
+
             # For t-channels, just need the first non-empty one
             tchannels = [t for s,t in stchannels if t != None][0]
-
+                 
+            # pass to ping-pong strategy for t-channel for 3 ore more T-channel
+            #  this is directly related to change in genps.f
+            tstrat = self.opt.get('t_strategy', 0)
+            tchannels, tchannels_strategy = ProcessExporterFortranME.reorder_tchannels(tchannels, tstrat, self.model)
+            
             # For s_and_t_channels (to be used later) use only first config
             s_and_t_channels.append([[s for s,t in stchannels if t != None][0],
-                                     tchannels])
+                                     tchannels, tchannels_strategy])
 
             # Make sure empty_verts is same length as real vertices
             if any([s for s,t in stchannels]):
@@ -4663,6 +4865,7 @@
             # Correspondance between the config and the diagram = amp2
             lines.append("data mapconfig(%d)/%d/" % (nconfigs,
                                                      mapconfigs[iconfig]))
+            lines.append("data tstrategy(%d)/%d/" % (nconfigs, tchannels_strategy))
             # Number of QCD couplings in this diagram
             nqcd = 0
             for h in helas_diags:
@@ -4715,7 +4918,321 @@
         writer.writelines(lines)
 
         return s_and_t_channels, nqcd_list
-
+    
+
+
+    #===========================================================================
+    # reoder t-channels
+    #===========================================================================
+    
+    #ordering = 0    
+    @staticmethod
+    def reorder_tchannels(tchannels, tstrat, model):
+        # no need to modified anything if 1 or less T-Channel
+        #Note that this counts the number of vertex (one more vertex compare to T)
+        #ProcessExporterFortranME.ordering +=1
+        if len(tchannels) < 3 or tstrat == 2 or not model:
+            return tchannels, 2
+        elif tstrat == 1:
+            return ProcessExporterFortranME.reorder_tchannels_flipside(tchannels), 1
+        elif tstrat == -2:
+            return ProcessExporterFortranME.reorder_tchannels_pingpong(tchannels), -2
+        elif tstrat == -1:
+            return ProcessExporterFortranME.reorder_tchannels_pingpong(tchannels, 1), -1        
+        elif len(tchannels) < 4:
+            #
+            first = tchannels[0]['legs'][1]['number']
+            t1 =  tchannels[0]['legs'][-1]['id']
+            last = tchannels[-1]['legs'][1]['number']
+            t2 = tchannels[-1]['legs'][0]['id']
+            m1  = model.get_particle(t1).get('mass') == 'ZERO'
+            m2  = model.get_particle(t2).get('mass') == 'ZERO'
+            if m2 and not m1:
+                return ProcessExporterFortranME.reorder_tchannels_flipside(tchannels), 1
+            elif m1 and not m2:
+                return tchannels, 2
+            elif first < last:
+                return ProcessExporterFortranME.reorder_tchannels_flipside(tchannels), 1
+            else:
+                return tchannels, 2 
+        else:
+            first = tchannels[0]['legs'][1]['number']
+            t1 =  tchannels[0]['legs'][-1]['id']
+            last = tchannels[-1]['legs'][1]['number']
+            t2 = tchannels[-1]['legs'][0]['id']
+            m1  = model.get_particle(t1).get('mass') == 'ZERO'
+            m2  = model.get_particle(t2).get('mass') == 'ZERO'
+            
+            t12 =  tchannels[1]['legs'][-1]['id']
+            m12 = model.get_particle(t12).get('mass') == 'ZERO'
+            t22 = tchannels[-2]['legs'][0]['id']
+            m22 = model.get_particle(t22).get('mass') == 'ZERO'
+            if m2 and not m1:
+                if m22:
+                    return ProcessExporterFortranME.reorder_tchannels_flipside(tchannels), 1
+                else:
+                    return ProcessExporterFortranME.reorder_tchannels_pingpong(tchannels), -2
+            elif m1 and not m2:
+                if m12:
+                    return tchannels, 2
+                else:
+                    return ProcessExporterFortranME.reorder_tchannels_pingpong(tchannels), -2
+            elif m1 and m2 and  len(tchannels) == 4 and not m12: # 3 T propa
+                return ProcessExporterFortranME.reorder_tchannels_pingpong(tchannels), -2
+                # this case seems quite sensitive we tested method 2 specifically and this was not helping in general 
+            elif not m1 and not m2 and  len(tchannels) == 4 and m12:
+                if first < last:
+                    return ProcessExporterFortranME.reorder_tchannels_flipside(tchannels), 1
+                return tchannels, 2
+            else:
+                return ProcessExporterFortranME.reorder_tchannels_pingpong(tchannels), -2
+
+
+                
+
+    @staticmethod
+    def reorder_tchannels_flipside(tchannels):
+        """change the tchannel ordering to pass to a ping-pong strategy.
+           assume ninitial == 2
+        
+        We assume that we receive something like this
+        
+        1 ----- X ------- -2
+                |
+                | (-X) 
+                |
+                X -------- 4
+                | 
+                | (-X-1)
+                |
+                X --------- -1
+
+                X----------  3
+                | 
+                | (-N+2)
+                |                
+                X --------- L
+                |
+                | (-N+1) 
+                |                
+        -N ----- X ------- P        
+        
+        coded as 
+        (1 -2 > -X) (-X 4 > -X-1) (-X-1 -1 > -X-2) ...
+        ((-N+3) 3 > (-N+2)) ((-n+2) L > (-n+1)) ((-n+1) P > -N)
+        
+        we want to convert this as:
+        -N ----- X ------- -2
+                |
+                | (-N+1) 
+                |
+                X -------- 4
+                | 
+                | (-N+2)
+                |
+                X --------- -1
+
+                X----------  3
+                | 
+                | (-X-1)
+                |                
+                X --------- L
+                |
+                | (-X) 
+                |                
+        2 ----- X ------- P          
+        
+        coded as 
+        ( 2 P > -X) (-X L > -X-1) (-X-1 3 > -X-2)... (-X-L -2 > -N)
+        """
+        
+        # no need to modified anything if 1 or less T-Channel
+        #Note that this counts the number of vertex (one more vertex compare to T)
+        if len(tchannels) < 2:
+            return tchannels
+
+        out = []
+        oldid2new = {}
+        
+        # initialisation
+        # id of the first T-channel (-X)
+        propa_id = tchannels[0]['legs'][-1]['number'] 
+        #
+        # Setup the last vertex to refenence the second id beam
+        # -N (need to setup it to 2.
+        initialid = tchannels[-1]['legs'][-1]['number']       
+        oldid2new[initialid] = 2
+        oldid2new[1] = initialid
+            
+        i = 0 
+        while tchannels:
+            old_vert = tchannels.pop()
+                
+            #copy the vertex /leglist to avoid side effects
+            new_vert = base_objects.Vertex(old_vert)
+            new_vert['legs'] = base_objects.LegList([base_objects.Leg(l) for l in old_vert['legs']])
+            # vertex taken from the bottom we have 
+            # (-N+1 X > -N) we need to flip to pass to 
+            # -N X > -N+1 (and then relabel -N and -N+1  
+            legs = new_vert['legs'] # shorcut
+            id1 = legs[0]['number']
+            id2 = legs[1]['number']
+            id3 = legs[2]['number']
+            # to be secure  we also support (X -N+1 > -N)
+            if id3 == id2 -1 and id1 !=1:
+                legs[0], legs[1] = legs[1], legs[0]
+            #flipping side
+            legs[0], legs[2] = legs[2], legs[0]
+
+            # the only new relabelling is the last element of the list
+            # always thanks to the above flipping
+            old_propa_id = new_vert['legs'][-1]['number'] 
+            oldid2new[old_propa_id] = propa_id
+
+            
+            #pass to new convention for leg numbering:
+            for l in new_vert['legs']:
+                if l['number'] in  oldid2new:
+                    l['number'] = oldid2new[l['number']]  
+                    
+            # new_vert is now ready
+            out.append(new_vert)
+            # prepare next iteration
+            propa_id -=1
+            i +=1
+
+        return out
+    
+    @staticmethod
+    def reorder_tchannels_pingpong(tchannels, id=2):
+        """change the tchannel ordering to pass to a ping-pong strategy.
+           assume ninitial == 2
+        
+        We assume that we receive something like this
+        
+        1 ----- X ------- -2
+                |
+                | (-X) 
+                |
+                X -------- 4
+                | 
+                | (-X-1)
+                |
+                X --------- -1
+
+                X----------  3
+                | 
+                | (-N+2)
+                |                
+                X --------- L
+                |
+                | (-N+1) 
+                |                
+        -N ----- X ------- P        
+        
+        coded as 
+        (1 -2 > -X) (-X 4 > -X-1) (-X-1 -1 > -X-2) ...
+        ((-N+3) 3 > (-N+2)) ((-n+2) L > (-n+1)) ((-n+1) P > -N)
+        
+        we want to convert this as:
+        1 ----- X ------- -2
+                |
+                | (-X) 
+                |
+                X -------- 4
+                | 
+                | (-X-2)
+                |
+                X --------- -1
+
+                X----------  3
+                | 
+                | (-X-3)
+                |                
+                X --------- L
+                |
+                | (-X-1) 
+                |                
+        2 ----- X ------- P          
+        
+        coded as 
+        (1 -2 > -X) (2 P > -X-1) (-X 4 > -X-2) (-X-1 L > -X-3) ...
+        """
+
+        # no need to modified anything if 1 or less T-Channel
+        #Note that this counts the number of vertex (one more vertex compare to T)
+        if len(tchannels) < 2:
+            return tchannels
+
+        out = []
+        oldid2new = {}
+        
+        # initialisation
+        # id of the first T-channel (-X)
+        propa_id = tchannels[0]['legs'][-1]['number'] 
+        #
+        # Setup the last vertex to refenence the second id beam
+        # -N (need to setup it to 2.
+        initialid = tchannels[-1]['legs'][-1]['number']       
+        oldid2new[initialid] = id
+
+
+        
+        i = 0 
+        while tchannels:
+            #ping pong by taking first/last element in aternance
+            if id ==2:
+                if i % 2 == 0:
+                    old_vert = tchannels.pop(0)
+                else:
+                    old_vert = tchannels.pop()
+            else:
+                if i % 2 != 0:
+                    old_vert = tchannels.pop(0)
+                else:
+                    old_vert = tchannels.pop()
+                    
+            #copy the vertex /leglist to avoid side effects
+            new_vert = base_objects.Vertex(old_vert)
+            new_vert['legs'] = base_objects.LegList([base_objects.Leg(l) for l in old_vert['legs']])
+            # if vertex taken from the bottom we have 
+            # (-N+1 X > -N) we need to flip to pass to 
+            # -N X > -N+1 (and then relabel -N and -N+1
+            # to be secure  we also support (X -N+1 > -N)
+            if (i % 2 ==1 and id ==2) or (i %2 == 0 and id ==1): 
+                legs = new_vert['legs'] # shorcut
+                id1 = legs[0]['number']
+                id2 = legs[1]['number'] 
+                if id1 > id2:
+                    legs[0], legs[1] = legs[1], legs[0]
+                else:
+                    legs[0], legs[2] = legs[2], legs[0]
+            
+            # the only new relabelling is the last element of the list
+            # always thanks to the above flipping
+            old_propa_id = new_vert['legs'][-1]['number'] 
+            oldid2new[old_propa_id] = propa_id
+
+            if i==0 and id==1:
+                legs[0]['number'] = 2
+            
+            #pass to new convention for leg numbering:
+            for l in new_vert['legs']:
+                if l['number'] in  oldid2new:
+                    l['number'] = oldid2new[l['number']]    
+            
+            # new_vert is now ready
+            out.append(new_vert)
+            # prepare next iteration
+            propa_id -=1
+            i +=1
+
+        return out
+
+            
+        
+        
+    
     #===========================================================================
     # write_decayBW_file
     #===========================================================================
@@ -4828,8 +5345,9 @@
     def write_symmetry(self, writer, v5=True):
         """Write the SubProcess/driver.f file for ME"""
 
+        
         path = pjoin(_file_path,'iolibs','template_files','madevent_symmetry.f')
-        
+
         if self.model_name == 'mssm' or self.model_name.startswith('mssm-'):
             card = 'Source/MODEL/MG5_param.dat'
         else:
@@ -5036,6 +5554,14 @@
 
     matrix_file = "matrix_madevent_group_v4.inc"
     grouped_mode = 'madevent'
+    default_opt = {'clean': False, 'complex_mass':False,
+                        'export_format':'madevent', 'mp': False,
+                        'v5_model': True,
+                        'output_options':{},
+                        'hel_recycling': True
+                        }
+    
+    
     #===========================================================================
     # generate_subprocess_directory
     #===========================================================================
@@ -5096,16 +5622,54 @@
         self.write_driver(writers.FortranWriter(filename),ncomb,
                                   n_grouped_proc=len(matrix_elements), v5=self.opt['v5_model'])
 
+        self.proc_characteristic['hel_recycling'] = self.opt['hel_recycling']
         for ime, matrix_element in \
                 enumerate(matrix_elements):
-            filename = 'matrix%d.f' % (ime+1)
-            calls, ncolor = \
-               self.write_matrix_element_v4(writers.FortranWriter(filename), 
-                            matrix_element,
-                            fortran_model,
-                            proc_id=str(ime+1),
-                            config_map=subproc_group.get('diagram_maps')[ime],
-                            subproc_number=group_number)
+            if self.opt['hel_recycling']:
+                filename = 'matrix%d_orig.f' % (ime+1)
+                replace_dict = self.write_matrix_element_v4(None, 
+                                matrix_element,
+                                fortran_model,
+                                proc_id=str(ime+1),
+                                config_map=subproc_group.get('diagram_maps')[ime],
+                                subproc_number=group_number)
+                calls,ncolor = replace_dict['return_value']
+                tfile = open(replace_dict['template_file']).read()
+                file = tfile % replace_dict
+                # Add the split orders helper functions.
+                file = file + '\n' + open(replace_dict['template_file2'])\
+                                                            .read()%replace_dict
+                # Write the file
+                writer = writers.FortranWriter(filename)
+                writer.writelines(file)
+                
+                #
+                # write the dedicated template for helicity recycling
+                #
+                tfile = open(replace_dict['template_file'].replace('.inc',"_hel.inc")).read() 
+                file = tfile % replace_dict
+                # Add the split orders helper functions.
+                file = file + '\n' + open(replace_dict['template_file2'])\
+                                                            .read()%replace_dict
+                # Write the file
+                writer = writers.FortranWriter('template_matrix%d.f' % (ime+1))
+                writer.uniformcase = False
+                writer.writelines(file)
+                
+                
+                
+                
+            else:
+                filename = 'matrix%d.f' % (ime+1)
+                calls, ncolor = \
+                   self.write_matrix_element_v4(writers.FortranWriter(filename), 
+                                matrix_element,
+                                fortran_model,
+                                proc_id=str(ime+1),
+                                config_map=subproc_group.get('diagram_maps')[ime],
+                                subproc_number=group_number)
+
+
 
             filename = 'auto_dsig%d.f' % (ime+1)
             self.write_auto_dsig_file(writers.FortranWriter(filename),
@@ -5323,6 +5887,13 @@
         replace_dict['nb_spin_state1'] = s1
         replace_dict['nb_spin_state2'] = s2
         
+        printzeroamp = []
+        for iproc in range(len(matrix_elements)):
+            printzeroamp.append(\
+                "        call print_zero_amp_%i()" % ( iproc + 1))
+        replace_dict['print_zero_amp'] = "\n".join(printzeroamp)
+        
+        
         if writer:
             file = open(pjoin(_file_path, \
                        'iolibs/template_files/super_auto_dsig_group_v4.inc')).read()

=== modified file 'madgraph/iolibs/file_writers.py'
--- madgraph/iolibs/file_writers.py	2020-11-25 20:43:03 +0000
+++ madgraph/iolibs/file_writers.py	2021-01-30 10:13:29 +0000
@@ -202,6 +202,7 @@
     number_re = re.compile('^(?P<num>\d+)\s+(?P<rest>.*)')
     line_cont_char = '$'
     comment_char = 'c'
+    uniformcase = True #force everyting to be lower/upper case 
     downcase = False
     line_length = 71
     max_split = 20
@@ -257,23 +258,24 @@
             # Replace all double quotes by single quotes
             myline = myline.replace('\"', '\'')
             # Downcase or upcase Fortran code, except for quotes
-            splitline = myline.split('\'')
-            myline = ""
-            i = 0
-            while i < len(splitline):
-                if i % 2 == 1:
-                    # This is a quote - check for escaped \'s
-                    while  splitline[i] and splitline[i][-1] == '\\':
-                        splitline[i] = splitline[i] + '\'' + splitline.pop(i + 1)
-                else:
-                    # Otherwise downcase/upcase
-                    if FortranWriter.downcase:
-                        splitline[i] = splitline[i].lower()
+            if self.uniformcase:
+                splitline = myline.split('\'')
+                myline = ""
+                i = 0
+                while i < len(splitline):
+                    if i % 2 == 1:
+                        # This is a quote - check for escaped \'s
+                        while  splitline[i] and splitline[i][-1] == '\\':
+                            splitline[i] = splitline[i] + '\'' + splitline.pop(i + 1)
                     else:
-                        splitline[i] = splitline[i].upper()
-                i = i + 1
-
-            myline = "\'".join(splitline).rstrip()
+                        # Otherwise downcase/upcase
+                        if FortranWriter.downcase:
+                            splitline[i] = splitline[i].lower()
+                        else:
+                            splitline[i] = splitline[i].upper()
+                    i = i + 1
+    
+                myline = "\'".join(splitline).rstrip()
 
             # Check if line starts with dual keyword and adjust indent 
             if self.__keyword_list and re.search(self.keyword_pairs[\

=== modified file 'madgraph/iolibs/gen_infohtml.py'
--- madgraph/iolibs/gen_infohtml.py	2019-04-17 14:39:47 +0000
+++ madgraph/iolibs/gen_infohtml.py	2021-01-30 10:13:29 +0000
@@ -243,7 +243,8 @@
         nb_diag = 0
                 
         pat = re.compile(r'''Amplitude\(s\) for diagram number (\d+)''' )
-       
+        if not os.path.exists(path):
+            path = os.path.join(self.dir, 'SubProcesses', proc, 'matrix%s_orig.f' % id)
         text = open(path).read()
         for match in re.finditer(pat, text):
             pass

=== modified file 'madgraph/iolibs/template_files/auto_dsig_v4.inc'
--- madgraph/iolibs/template_files/auto_dsig_v4.inc	2019-06-10 19:05:58 +0000
+++ madgraph/iolibs/template_files/auto_dsig_v4.inc	2021-01-30 10:13:29 +0000
@@ -156,3 +156,9 @@
 %(read_write_good_hel)s
 
 ## }
+
+
+     SUBROUTINE PRINT_ZERO_AMP%(proc_id)s()
+
+     return
+     end
\ No newline at end of file

=== modified file 'madgraph/iolibs/template_files/b_sf_xxx_fks.inc'
--- madgraph/iolibs/template_files/b_sf_xxx_fks.inc	2016-09-08 23:15:34 +0000
+++ madgraph/iolibs/template_files/b_sf_xxx_fks.inc	2021-01-30 10:13:29 +0000
@@ -99,6 +99,7 @@
       INTEGER I,J
       REAL*8 DENOM(NCOLOR1), CF(NCOLOR2,NCOLOR1)
       COMPLEX*16 ZTEMP, AMP(NGRAPHS), JAMP1(NCOLOR1), JAMP2(NCOLOR2)
+      COMPLEX*16 TMP_JAMP(%(nb_temp_jamp)i)
 C  
 C GLOBAL VARIABLES
 C  

=== modified file 'madgraph/iolibs/template_files/born_fks.inc'
--- madgraph/iolibs/template_files/born_fks.inc	2020-01-16 13:45:20 +0000
+++ madgraph/iolibs/template_files/born_fks.inc	2021-01-30 10:13:29 +0000
@@ -166,8 +166,9 @@
       INTEGER IC(NEXTERNAL-1),nmo
       parameter (nmo=nexternal-1)
       data ic /nmo*1/
-      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
+      REAL*8 CF(NCOLOR,NCOLOR)
       COMPLEX*16 ZTEMP, AMP(NGRAPHS), JAMP(NCOLOR), W(%(wavefunctionsize)d,NWAVEFUNCS), jamph(2, ncolor)
+      COMPLEX*16 TMP_JAMP(%(nb_temp_jamp)i)
 C  
 C GLOBAL VARIABLES
 C  
@@ -249,7 +250,7 @@
           DO J = 1, NCOLOR
             ZTEMP = ZTEMP + CF(J,I)*JAMP(J)
           ENDDO
-          BORNS(2-(1+back_hel*ihel)/2)=BORNS(2-(1+back_hel*ihel)/2)+ZTEMP*DCONJG(JAMP(I))/DENOM(I)   
+          BORNS(2-(1+back_hel*ihel)/2)=BORNS(2-(1+back_hel*ihel)/2)+ZTEMP*DCONJG(JAMP(I))   
         ENDDO
         Do I = 1, NGRAPHS
           amp2(i)=amp2(i)+amp(i)*dconjg(amp(i))
@@ -266,7 +267,7 @@
         DO J = 1, NCOLOR
           ZTEMP = ZTEMP + CF(J,I)*JAMPH(2,J)
         ENDDO
-        BORNTILDE = BORNTILDE + ZTEMP*DCONJG(JAMPH(1,I))/DENOM(I)
+        BORNTILDE = BORNTILDE + ZTEMP*DCONJG(JAMPH(1,I))
       ENDDO
       if (glu_ij.ne.0) nhel(glu_ij) = back_hel
       END

=== modified file 'madgraph/iolibs/template_files/born_fks_hel.inc'
--- madgraph/iolibs/template_files/born_fks_hel.inc	2016-09-08 23:15:34 +0000
+++ madgraph/iolibs/template_files/born_fks_hel.inc	2021-01-30 10:13:29 +0000
@@ -104,6 +104,7 @@
       INTEGER I,J
       REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
       COMPLEX*16 ZTEMP, AMP(NGRAPHS), JAMP(NCOLOR)
+      COMPLEX*16 TMP_JAMP(%(nb_temp_jamp)i)
 C  
 C GLOBAL VARIABLES
 C  

=== modified file 'madgraph/iolibs/template_files/check_sa.f'
--- madgraph/iolibs/template_files/check_sa.f	2015-06-18 08:24:43 +0000
+++ madgraph/iolibs/template_files/check_sa.f	2021-01-30 10:13:29 +0000
@@ -177,7 +177,7 @@
             P(2,2)=0d0
             P(3,2)=-mom
              
-            call rambo(nexternal-2,energy,pmass(3),prambo,WGT)
+            call rambo(nexternal-2,energy,pmass(nincoming+1),prambo,WGT)
             DO I=3, NEXTERNAL
                P(0,I)=PRAMBO(4,I-2)	
                P(1,I)=PRAMBO(1,I-2)

=== modified file 'madgraph/iolibs/template_files/check_sa_splitOrders.f'
--- madgraph/iolibs/template_files/check_sa_splitOrders.f	2015-04-09 22:42:44 +0000
+++ madgraph/iolibs/template_files/check_sa_splitOrders.f	2021-01-30 10:13:29 +0000
@@ -201,7 +201,7 @@
             P(2,2)=0d0
             P(3,2)=-mom
              
-            call rambo(nexternal-2,energy,pmass(3),prambo,WGT)
+            call rambo(nexternal-2,energy,pmass(nincoming+1),prambo,WGT)
             DO I=3, NEXTERNAL
                P(0,I)=PRAMBO(4,I-2)	
                P(1,I)=PRAMBO(1,I-2)

=== modified file 'madgraph/iolibs/template_files/loop/check_sa.inc'
--- madgraph/iolibs/template_files/loop/check_sa.inc	2016-05-22 02:36:25 +0000
+++ madgraph/iolibs/template_files/loop/check_sa.inc	2021-01-30 10:13:29 +0000
@@ -402,7 +402,7 @@
             P(2,2)=0d0
             P(3,2)=-mom
              
-            call rambo(nexternal-2,energy,pmass(3),prambo,WGT)
+            call rambo(nexternal-2,energy,pmass(nexternal+1),prambo,WGT)
             DO I=3, NEXTERNAL
                P(0,I)=PRAMBO(4,I-2)      
                P(1,I)=PRAMBO(1,I-2)

=== modified file 'madgraph/iolibs/template_files/loop/check_sa_loop_induced.inc'
--- madgraph/iolibs/template_files/loop/check_sa_loop_induced.inc	2016-05-22 02:36:25 +0000
+++ madgraph/iolibs/template_files/loop/check_sa_loop_induced.inc	2021-01-30 10:13:29 +0000
@@ -360,7 +360,7 @@
             P(2,2)=0d0
             P(3,2)=-mom
              
-            call rambo(nexternal-2,energy,pmass(3),prambo,WGT)
+            call rambo(nexternal-2,energy,pmass(nexternal+1),prambo,WGT)
             DO I=3, NEXTERNAL
                P(0,I)=PRAMBO(4,I-2)      
                P(1,I)=PRAMBO(1,I-2)

=== modified file 'madgraph/iolibs/template_files/madevent_driver.f'
--- madgraph/iolibs/template_files/madevent_driver.f	2017-06-08 18:43:40 +0000
+++ madgraph/iolibs/template_files/madevent_driver.f	2021-01-30 10:13:29 +0000
@@ -27,6 +27,7 @@
       external NextUnopen
       double precision t_before
       logical fopened
+      integer nb_tchannel
 c
 c     Global
 c
@@ -159,7 +160,7 @@
       maxcfig=mincfig
       minvar(1,1) = 0              !This tells it to map things invarients
       write(*,*) 'Attempting mappinvarients',nconfigs,nexternal
-      call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming)
+      call map_invarients(minvar,nconfigs,ninvar,mincfig,maxcfig,nexternal,nincoming,nb_tchannel)
       write(*,*) "Completed mapping",nexternal
       ndim = 3*(nexternal-nincoming)-4
       if (nincoming.gt.1.and.abs(lpp(1)) .ge. 1) ndim=ndim+1
@@ -173,6 +174,12 @@
             minvar(ndim,j) = ninvar
          endif
       enddo
+c      ncall =  ncall * max(1., min(3., (nb_tchannel+1.)/2.))
+      if (nb_tchannel.gt.1) then
+c         itmin = itmin + 1
+         itmax = itmax + 2
+      endif
+
       write(*,*) "about to integrate ", ndim,ncall,itmax,itmin,ninvar,nconfigs
       call sample_full(ndim,ncall,itmax,itmin,dsig,ninvar,nconfigs)
 
@@ -233,6 +240,9 @@
       common /to_accuracy/accur
       integer           use_cut
       common /to_weight/use_cut
+      logical init_mode
+      common /to_determine_zero_hel/init_mode
+
 
       integer        lbw(0:nexternal)  !Use of B.W.
       common /to_BW/ lbw
@@ -270,11 +280,17 @@
          write(*,*) 'Using full amplitude.'
       endif
 
+      init_mode = .false.
       write(*,10) 'Exact helicity sum (0 yes, n = number/event)? '
       read(*,*) i
       if (i .eq. 0) then
          isum_hel = 0
          write(*,*) 'Explicitly summing over helicities'
+      else if (i.eq.-1) then
+         isum_hel = 0
+         multi_channel = .false.
+         init_mode = .true.
+         write(*,*) 'Determining zero helicities'
       else
          isum_hel= i
          write(*,*) 'Monte-Carlo over helicities'

=== modified file 'madgraph/iolibs/template_files/madevent_symmetry.f'
--- madgraph/iolibs/template_files/madevent_symmetry.f	2020-04-16 19:10:18 +0000
+++ madgraph/iolibs/template_files/madevent_symmetry.f	2021-01-30 10:13:29 +0000
@@ -47,6 +47,7 @@
 c     DATA
 c
       integer iforest(2,-max_branch:-1,lmaxconfigs)
+      integer tstrategy(lmaxconfigs)
       integer sprop(maxsproc,-max_branch:-1,lmaxconfigs)
       integer tprid(-max_branch:-1,lmaxconfigs)
       include 'configs.inc'
@@ -354,7 +355,8 @@
 c     chain and could potentially conflict
 c
       i=1
-      do while (i .lt. nexternal-2 .and. itree(1,-i) .ne. 1)
+      do while (i .lt. nexternal-2 .and. .not. (
+     &    itree(1,-i) .eq. 1.or.(nincoming.eq.2.and.itree(1,-i).eq.2)))
          xmass(-i) = xmass(itree(1,-i))+xmass(itree(2,-i))
          mtot=mtot-xmass(-i)
          if (prwidth(-i,iconfig) .gt. 0d0) then

=== modified file 'madgraph/iolibs/template_files/matrix_madevent_group_v4.inc'
--- madgraph/iolibs/template_files/matrix_madevent_group_v4.inc	2020-10-26 15:18:03 +0000
+++ madgraph/iolibs/template_files/matrix_madevent_group_v4.inc	2021-01-30 10:13:29 +0000
@@ -20,7 +20,7 @@
     Include 'nexternal.inc'
     Include 'maxamps.inc'
     INTEGER                 NCOMB         
-    PARAMETER (             NCOMB=%(ncomb)d)
+    PARAMETER (             NCOMB=%(ncomb)s)
     INTEGER    NGRAPHS
     PARAMETER (NGRAPHS=%(ngraphs)d) 
     INTEGER    NDIAGS
@@ -40,7 +40,9 @@
 
 C  
 C LOCAL VARIABLES 
-C  
+C
+        INTEGER CONFSUB(MAXSPROC,LMAXCONFIGS)
+        INCLUDE 'config_subproc_map.inc'
         INTEGER NHEL(NEXTERNAL,NCOMB)
         INTEGER ISHEL(2)
         REAL*8 T,MATRIX%(proc_id)s
@@ -56,6 +58,8 @@
 	SAVE NB_FAIL
 	DATA THIS_NTRY /0,0/
 	DATA NB_FAIL /0/
+        double precision get_channel_cut
+        external get_channel_cut
 	
 c
 C       This is just to temporarily store the reference grid for helicity of the DiscreteSampler so as to obtain its number of entries with ref_helicity_grid%%n_tot_entries
@@ -63,6 +67,8 @@
 C  
 C GLOBAL VARIABLES
 C  
+    logical init_mode
+    common /to_determine_zero_hel/init_mode
     DOUBLE PRECISION AMP2(MAXAMPS), JAMP2(0:MAXFLOW)
     COMMON/TO_AMPS/  AMP2,       JAMP2
    
@@ -74,7 +80,12 @@
     
     INTEGER IMIRROR
     COMMON/TO_MIRROR/ IMIRROR
- 
+
+       double precision tmin_for_channel	
+       integer sde_strat ! 1 means standard single diagram enhancement strategy,
+c      	      	      	   2 means approximation by the	denominator of the propagator
+       common/TO_CHANNEL_STRAT/tmin_for_channel,	sde_strat
+
     REAL*8 POL(2)
     COMMON/TO_POLARIZATION/ POL
 
@@ -125,8 +136,8 @@
 !   If HEL_PICKED==-1, this means that calls to other matrix<i> where in initialization mode as well for the helicity.
     IF ((ISHEL(IMIRROR).EQ.0.and.ISUM_HEL.eq.0).or.(DS_get_dim_status('Helicity').eq.0).or.(HEL_PICKED.eq.-1)) THEN
       DO I=1,NCOMB
-         IF (GOODHEL(I,IMIRROR) .OR. NTRY(IMIRROR).LE.MAXTRIES.or.(ISUM_HEL.NE.0).or.THIS_NTRY(IMIRROR).le.2) THEN
-            T=MATRIX%(proc_id)s(P ,NHEL(1,I),JC(1))
+         IF (GOODHEL(I,IMIRROR) .OR. NTRY(IMIRROR).LE.MAXTRIES.or.(ISUM_HEL.NE.0).or.THIS_NTRY(IMIRROR).le.10) THEN
+            T=MATRIX%(proc_id)s(P ,NHEL(1,I),JC(1),I)
             DO JJ=1,nincoming
                IF(POL(JJ).NE.1d0.AND.NHEL(JJ,I).EQ.INT(SIGN(1d0,POL(JJ)))) THEN
                   T=T*ABS(POL(JJ))*NB_SPIN_STATE_IN(JJ)/2d0 ! NB_SPIN_STATE(JJ)/2d0 is added for polarised beam
@@ -159,9 +170,13 @@
           endif
 	  ELSE
         JHEL(IMIRROR) = 1
-        IF(NTRY(IMIRROR).LE.MAXTRIES.or.THIS_NTRY(IMIRROR).le.2)THEN
+        IF(NTRY(IMIRROR).LE.MAXTRIES.or.THIS_NTRY(IMIRROR).le.10)THEN
            DO I=1,NCOMB
-              IF (.NOT.GOODHEL(I,IMIRROR) .AND. (DABS(TS(I)).GT.ANS*LIMHEL/NCOMB)) THEN
+	      IF(init_mode) THEN
+                  IF (DABS(TS(I)).GT.ANS*LIMHEL/NCOMB) THEN
+                     PRINT *, 'Matrix Element/Good Helicity: %(proc_id)s ', i, 'IMIRROR', IMIRROR
+	          ENDIF
+	      ELSE IF (.NOT.GOODHEL(I,IMIRROR) .AND. (DABS(TS(I)).GT.ANS*LIMHEL/NCOMB)) THEN
                  GOODHEL(I,IMIRROR)=.TRUE.
                  NGOOD(IMIRROR) = NGOOD(IMIRROR) +1
                  IGOOD(NGOOD(IMIRROR),IMIRROR) = I
@@ -173,11 +188,11 @@
            ISHEL(IMIRROR)=MIN(ISUM_HEL,NGOOD(IMIRROR))
         ENDIF
       ENDIF
-    ELSE             ! random helicity 
+    ELSE IF (.not.init_mode) then            ! random helicity 
 C           The helicity configuration was chosen already by genps and put in a common block defined in genps.inc.
             I = HEL_PICKED
 
-            T=MATRIX%(proc_id)s(P ,NHEL(1,I),JC(1))
+            T=MATRIX%(proc_id)s(P ,NHEL(1,I),JC(1),I)
 
             DO JJ=1,nincoming
               IF(POL(JJ).NE.1d0.AND.NHEL(JJ,I).EQ.INT(SIGN(1d0,POL(JJ)))) THEN
@@ -192,6 +207,9 @@
             ANS = ANS * hel_jacobian
 
           WRITE(HEL_BUFF,'(20i5)')(NHEL(II,I),II=1,NEXTERNAL)
+    ELSE
+	ANS = 1D0
+	RETURN
     ENDIF
     IF (ANS.ne.0d0.and.(ISUM_HEL .NE. 1.or.HEL_PICKED.eq.-1)) THEN
     CALL RANMAR(R)
@@ -209,8 +227,17 @@
     ENDIF
     IF (MULTI_CHANNEL) THEN
         XTOT=0D0
-        DO I=1,NDIAGS
-            XTOT=XTOT+AMP2(I)
+        DO I=1,LMAXCONFIGS
+	    J = CONFSUB(%(proc_id)s, I)
+	    if (J.ne.0) then
+	    if(sde_strat.eq.1) then
+	         AMP2(J) = AMP2(J) * GET_CHANNEL_CUT(P, I)
+                 XTOT=XTOT+AMP2(J)
+	    else
+	    	 AMP2(J) = GET_CHANNEL_CUT(P, I)
+                 XTOT=XTOT+AMP2(J)
+	    endif	
+            endif	
         ENDDO
         IF (XTOT.NE.0D0) THEN
 %(set_amp2_line)s
@@ -228,7 +255,7 @@
     END
  
  
-REAL*8 FUNCTION MATRIX%(proc_id)s(P,NHEL,IC)
+REAL*8 FUNCTION MATRIX%(proc_id)s(P,NHEL,IC, IHEL)
 C  
 %(info_lines)s
 C
@@ -242,7 +269,9 @@
 C CONSTANTS
 C  
     INTEGER    NGRAPHS
-    PARAMETER (NGRAPHS=%(ngraphs)d) 
+    PARAMETER (NGRAPHS=%(ngraphs)d)
+    INTEGER                 NCOMB         
+    PARAMETER (             NCOMB=%(ncomb)s)
     include 'genps.inc'
     include 'nexternal.inc'
     include 'maxamps.inc'
@@ -262,12 +291,13 @@
 C  
     REAL*8 P(0:3,NEXTERNAL)
     INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
+    INTEGER IHEL
 C  
 C LOCAL VARIABLES 
 C  
     INTEGER I,J,M,N
-    COMPLEX*16 ZTEMP
-    REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
+    COMPLEX*16 ZTEMP, TMP_JAMP(%(nb_temp_jamp)i)
+    REAL*8 CF(NCOLOR,NCOLOR)
     COMPLEX*16 AMP(NGRAPHS), JAMP(NCOLOR,NAMPSO)
     COMPLEX*16 W(%(wavefunctionsize)d,NWAVEFUNCS)
 C   Needed for v4 models
@@ -291,6 +321,18 @@
 
 	double precision small_width_treatment
 	common/narrow_width/small_width_treatment
+
+	logical init_mode
+	common/to_determine_zero_hel/init_mode
+
+	logical zeroamp_%(proc_id)s(NCOMB,NGRAPHS)
+	common/to_zeroamp_%(proc_id)s/zeroamp_%(proc_id)s
+
+       double precision tmin_for_channel
+       integer sde_strat ! 1 means standard single diagram enhancement strategy,
+c      	      	      	   2 means approximation by the	denominator of the propagator
+       common/TO_CHANNEL_STRAT/tmin_for_channel,	sde_strat
+       
 C  
 C COLOR DATA
 C  
@@ -301,12 +343,24 @@
 if (first) then
 	first=.false.
 	%(fake_width_definitions)s
+
+	if(init_mode) then
+	   zeroamp_%(proc_id)s(:,:) = .true.
+	endif	
 endif
 
 
 %(helas_calls)s
 %(jamp_lines)s
 
+    if(init_mode)then
+      DO I=1, NGRAPHS	
+        if (AMP(I).ne.0) then
+       	  zeroamp_%(proc_id)s(ihel,i) = .false.
+       endif
+      ENDDO
+    endif	
+
     MATRIX%(proc_id)s = 0.D0 
 	DO M = 1, NAMPSO
       DO I = 1, NCOLOR
@@ -315,22 +369,62 @@
           ZTEMP = ZTEMP + CF(J,I)*JAMP(J,M)
         ENDDO
 		DO N = 1, NAMPSO
-          IF (CHOSEN_SO_CONFIGS(SQSOINDEX%(proc_id)s(M,N))) THEN
-		     MATRIX%(proc_id)s = MATRIX%(proc_id)s + ZTEMP*DCONJG(JAMP(I,N))/DENOM(I)
-		  ENDIF
+		  %(select_configs_if)s
+             MATRIX%(proc_id)s = MATRIX%(proc_id)s + ZTEMP*DCONJG(JAMP(I,N))
+	      %(select_configs_endif)s
 		ENDDO
       ENDDO
 	ENDDO
 
+    if(sde_strat.eq.1)then
 %(amp2_lines)s
+    endif 
+
     Do I = 1, NCOLOR
 	  DO M = 1, NAMPSO	
 		DO N = 1, NAMPSO
-          IF (CHOSEN_SO_CONFIGS(SQSOINDEX%(proc_id)s(M,N))) THEN
+		  %(select_configs_if)s
             Jamp2(i)=Jamp2(i)+DABS(DBLE(Jamp(i,m)*dconjg(Jamp(i,n))))
-		  ENDIF
+		  %(select_configs_endif)s
 		enddo
 	  enddo
     Enddo
 
     END
+
+    SUBROUTINE PRINT_ZERO_AMP_%(proc_id)s()
+
+    IMPLICIT NONE
+    INTEGER    NGRAPHS
+    PARAMETER (NGRAPHS=%(ngraphs)d)
+
+    INTEGER    NCOMB
+    PARAMETER (NCOMB=%(ncomb)d)
+
+    logical zeroamp_%(proc_id)s(NCOMB, NGRAPHS)
+    common/to_zeroamp_%(proc_id)s/zeroamp_%(proc_id)s
+
+    integer i,j
+    logical all_false
+
+    do i=1, NGRAPHS
+       all_false = .true.
+       do j=1,NCOMB
+       	  if (.not.zeroamp_%(proc_id)s(j, i)) then
+	     all_false = .false.
+	     exit
+	  endif
+       enddo	  
+       if (all_false) then
+       	  write(*,*) 'Amplitude/ZEROAMP:', %(proc_id)s, i
+       else
+          do j=1,NCOMB
+               if (zeroamp_%(proc_id)s(j, i)) then
+	         write(*,*) 'HEL/ZEROAMP:', %(proc_id)s, j  , i
+	      endif	    
+           enddo	
+       endif
+    enddo   
+
+    return
+    end
\ No newline at end of file

=== added file 'madgraph/iolibs/template_files/matrix_madevent_group_v4_hel.inc'
--- madgraph/iolibs/template_files/matrix_madevent_group_v4_hel.inc	1970-01-01 00:00:00 +0000
+++ madgraph/iolibs/template_files/matrix_madevent_group_v4_hel.inc	2021-01-30 10:13:29 +0000
@@ -0,0 +1,281 @@
+SUBROUTINE SMATRIX%(proc_id)s(P,ANS)
+C 
+%(info_lines)s
+C 
+C MadGraph5_aMC@NLO for Madevent Version
+C 
+C Returns amplitude squared summed/avg over colors
+c and helicities
+c for the point in phase space P(0:3,NEXTERNAL)
+C  
+%(process_lines)s
+C 
+    IMPLICIT NONE
+C  
+C CONSTANTS
+C  
+    Include 'genps.inc'
+    Include 'maxconfigs.inc'
+    Include 'nexternal.inc'
+    Include 'maxamps.inc'
+    INTEGER                 NCOMB         
+    PARAMETER (             NCOMB=${ncomb})
+    INTEGER    NGRAPHS
+    PARAMETER (NGRAPHS=%(ngraphs)d) 
+    INTEGER    NDIAGS
+    PARAMETER (NDIAGS=%(ndiags)d) 
+    INTEGER    THEL
+    PARAMETER (THEL=2*NCOMB)
+C  
+C ARGUMENTS 
+C  
+    REAL*8 P(0:3,NEXTERNAL),ANS
+C  
+C LOCAL VARIABLES 
+C  
+        INTEGER CONFSUB(MAXSPROC,LMAXCONFIGS)
+        INCLUDE 'config_subproc_map.inc'
+        INTEGER NHEL(NEXTERNAL,NCOMB)
+        REAL*8 T
+        REAL*8 R,SUMHEL,TS(NCOMB)
+        INTEGER I,IDEN
+        INTEGER JC(NEXTERNAL),II
+        REAL*8 XTOT
+        INTEGER  J, JJ
+
+      double precision get_channel_cut
+      external get_channel_cut
+
+C  
+C GLOBAL VARIABLES
+C  
+    DOUBLE PRECISION AMP2(MAXAMPS), JAMP2(0:MAXFLOW)
+    COMMON/TO_AMPS/  AMP2,       JAMP2
+
+
+C
+C   INFORMATION TO WRITE THE HELICITY IN THE EVENT --not memory efficient--
+C   
+    CHARACTER*101         HEL_BUFF
+    COMMON/TO_HELICITY/  HEL_BUFF
+
+    INTEGER NB_SPIN_STATE_in(2)
+    common /nb_hel_state/ nb_spin_state_in
+    
+    REAL*8 POL(2)
+    
+    COMMON/TO_POLARIZATION/ POL
+       double precision tmin_for_channel
+       integer sde_strat ! 1 means standard single diagram enhancement strategy,
+c      	      	      	   2 means approximation by the	denominator of the propagator
+       common/TO_CHANNEL_STRAT/tmin_for_channel,	sde_strat
+       
+    INTEGER          ISUM_HEL
+    LOGICAL                    MULTI_CHANNEL
+    COMMON/TO_MATRIX/ISUM_HEL, MULTI_CHANNEL
+%(define_iconfigs_lines)s
+${helicity_lines}
+%(den_factor_line)s
+
+C ----------
+C BEGIN CODE
+C ----------
+
+    DO I=1,NEXTERNAL
+       JC(I) = +1
+    ENDDO
+     
+    IF (multi_channel) THEN
+        DO I=1,NDIAGS
+            AMP2(I)=0D0
+        ENDDO
+        JAMP2(0)=%(ncolor)d
+        DO I=1,INT(JAMP2(0))
+            JAMP2(I)=0D0
+        ENDDO
+    ENDIF
+    ANS = 0D0
+    WRITE(HEL_BUFF,'(20I5)') (0,I=1,NEXTERNAL)
+c   Kiran please check if you need this:    
+    DO I=1,NCOMB
+       TS(I)=0d0
+    ENDDO
+
+    call MATRIX%(proc_id)s(P ,JC(1), TS)
+    DO I=1,NCOMB     
+       T=TS(I)  
+       DO JJ=1,nincoming
+          IF(POL(JJ).NE.1d0.AND.NHEL(JJ,I).EQ.INT(SIGN(1d0,POL(JJ)))) THEN
+             T=T*ABS(POL(JJ))*NB_SPIN_STATE_IN(JJ)/2d0 ! NB_SPIN_STATE(JJ)/2d0 is added for polarised beam
+          ELSE IF(POL(JJ).NE.1d0)THEN
+             T=T*(2d0-ABS(POL(JJ)))*NB_SPIN_STATE_IN(JJ)/2d0
+          ENDIF
+       ENDDO
+       ANS=ANS+DABS(T)
+       TS(I)=T
+    ENDDO
+
+    IF (ANS.ne.0d0) THEN
+        CALL RANMAR(R)
+        SUMHEL=0d0
+	DO I=1,NCOMB
+            SUMHEL=SUMHEL+DABS(TS(I))/ANS
+            IF(R.LT.SUMHEL)THEN
+                WRITE(HEL_BUFF,'(20i5)')(NHEL(II,I),II=1,NEXTERNAL)
+c         Set right sign for ANS, based on sign of chosen helicity
+                ANS=DSIGN(ANS,TS(I))
+                GOTO 10
+            ENDIF
+        ENDDO
+ 10     CONTINUE   
+    ENDIF
+    IF (MULTI_CHANNEL) THEN
+        XTOT=0D0
+        DO I=1,LMAXCONFIGS
+	    J = CONFSUB(%(proc_id)s, I)
+	    if (J.ne.0)then
+	      if (sde_strat.eq.1)then
+	         AMP2(J) = AMP2(J) * GET_CHANNEL_CUT(P, I)
+               else
+	     	         AMP2(J) = GET_CHANNEL_CUT(P, I)
+	       endif
+             XTOT=XTOT+AMP2(J)
+
+             endif
+        ENDDO
+        IF (XTOT.NE.0D0) THEN
+%(set_amp2_line)s
+        ELSE IF(ANS.ne.0d0) THEN
+			write(*,*) "Problem in the multi-channeling. All amp2 are zero but not the total matrix-element"
+	        stop 1
+        ENDIF
+    ENDIF
+    ANS=ANS/DBLE(IDEN)
+    END
+ 
+ 
+Subroutine  MATRIX%(proc_id)s(P,IC, TS)
+C  
+%(info_lines)s
+C
+C Returns amplitude squared summed/avg over colors
+c for the point with external lines W(0:6,NEXTERNAL)
+C  
+%(process_lines)s
+C  
+    IMPLICIT NONE
+C  
+C CONSTANTS
+C  
+    INTEGER    NGRAPHS
+    PARAMETER (NGRAPHS=%(ngraphs)d) 
+    include 'genps.inc'
+    include 'nexternal.inc'
+    include 'maxamps.inc'
+    INTEGER    NWAVEFUNCS,     NCOLOR
+    PARAMETER (NWAVEFUNCS=${nwavefuncs}, NCOLOR=%(ncolor)d) 
+    REAL*8     ZERO
+    PARAMETER (ZERO=0D0)
+    COMPLEX*16 IMAG1
+    PARAMETER (IMAG1=(0D0,1D0))
+    INTEGER NAMPSO, NSQAMPSO
+    PARAMETER (NAMPSO=%(nAmpSplitOrders)d, NSQAMPSO=%(nSqAmpSplitOrders)d)
+	LOGICAL CHOSEN_SO_CONFIGS(NSQAMPSO)
+	DATA CHOSEN_SO_CONFIGS/%(chosen_so_configs)s/
+	SAVE CHOSEN_SO_CONFIGS
+    INTEGER                 NCOMB         
+    PARAMETER (             NCOMB=${ncomb})
+C  
+C ARGUMENTS 
+C
+    REAL*8 P(0:3,NEXTERNAL)
+    INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
+    REAL*8 TS(NCOMB)
+C  
+C LOCAL VARIABLES 
+C  
+    INTEGER I,J,M,N,K
+    COMPLEX*16 ZTEMP,TMP_JAMP(%(nb_temp_jamp)i)
+    COMPLEX*16 TMP(6)
+    REAL*8 CF(NCOLOR,NCOLOR)
+    COMPLEX*16 AMP(NCOMB,NGRAPHS), JAMP(NCOLOR,NAMPSO)
+    COMPLEX*16 W(%(wavefunctionsize)d,NWAVEFUNCS)
+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
+      INTEGER SQSOINDEX%(proc_id)s
+
+C  
+C GLOBAL VARIABLES
+C  
+    Double Precision amp2(maxamps), jamp2(0:maxflow)
+    common/to_amps/  amp2,       jamp2
+    include 'coupl.inc'
+
+       double precision tmin_for_channel
+       integer sde_strat ! 1 means standard single diagram enhancement strategy,
+c      	      	      	   2 means approximation by the	denominator of the propagator
+       common/TO_CHANNEL_STRAT/tmin_for_channel,	sde_strat
+
+	double precision small_width_treatment
+	common/narrow_width/small_width_treatment
+C  
+C COLOR DATA
+C  
+%(color_data_lines)s
+C ----------
+C BEGIN CODE
+C ----------
+if (first) then
+	first=.false.
+	%(fake_width_definitions)s
+endif
+    AMP(:,:) = (0d0,0d0)
+${helas_calls}
+    DO K = 1, NCOMB
+${jamp_lines}
+        TS(K) = 0.D0 
+    	DO M = 1, NAMPSO
+          DO I = 1, NCOLOR
+            ZTEMP = (0.D0,0.D0)
+            DO J = 1, NCOLOR
+              ZTEMP = ZTEMP + CF(J,I)*JAMP(J,M)
+            ENDDO
+    		DO N = 1, NAMPSO
+			%(select_configs_if)s
+    		     TS(K) = TS(K) + ZTEMP*DCONJG(JAMP(I,N))
+			  %(select_configs_endif)s
+    		ENDDO
+          ENDDO
+    	ENDDO
+	if(sde_strat.eq.1) then
+        ${amp2_lines}
+	endif
+        Do I = 1, NCOLOR
+    	  DO M = 1, NAMPSO
+    		DO N = 1, NAMPSO
+			%(select_configs_if)s
+                Jamp2(i)=Jamp2(i)+DABS(DBLE(Jamp(i,m)*dconjg(Jamp(i,n))))
+			%(select_configs_endif)s
+    		enddo
+    	  enddo
+        Enddo
+    ENDDO
+
+    END
+
+
+    SUBROUTINE PRINT_ZERO_AMP_%(proc_id)s()
+
+    integer i
+    	    i =1
+    return
+    end
\ No newline at end of file

=== modified file 'madgraph/iolibs/template_files/matrix_madevent_v4.inc'
--- madgraph/iolibs/template_files/matrix_madevent_v4.inc	2020-10-26 15:18:03 +0000
+++ madgraph/iolibs/template_files/matrix_madevent_v4.inc	2021-01-30 10:13:29 +0000
@@ -233,8 +233,8 @@
 C LOCAL VARIABLES 
 C  
     INTEGER I,J,M,N
-    COMPLEX*16 ZTEMP
-    REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
+    COMPLEX*16 ZTEMP, TMP_JAMP(%(nb_temp_jamp)i)
+    REAL*8 CF(NCOLOR,NCOLOR)
     COMPLEX*16 AMP(NGRAPHS), JAMP(NCOLOR,NAMPSO)
     COMPLEX*16 W(18,NWAVEFUNCS)
 C   Needed for v4 models
@@ -282,7 +282,7 @@
         ENDDO
 		DO N = 1, NAMPSO
           IF (CHOSEN_SO_CONFIGS(SQSOINDEX%(proc_id)s(M,N))) THEN
-		     MATRIX%(proc_id)s = MATRIX%(proc_id)s + ZTEMP*DCONJG(JAMP(I,N))/DENOM(I)
+		     MATRIX%(proc_id)s = MATRIX%(proc_id)s + ZTEMP*DCONJG(JAMP(I,N))
 		  ENDIF
 		ENDDO
       ENDDO

=== modified file 'madgraph/iolibs/template_files/matrix_madweight_group_v4.inc'
--- madgraph/iolibs/template_files/matrix_madweight_group_v4.inc	2018-10-07 20:06:38 +0000
+++ madgraph/iolibs/template_files/matrix_madweight_group_v4.inc	2021-01-30 10:13:29 +0000
@@ -90,7 +90,7 @@
 C  
       INTEGER I,J
       COMPLEX*16 ZTEMP
-      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
+      REAL*8  CF(NCOLOR,NCOLOR)
       COMPLEX*16 AMP(NGRAPHS), JAMP(NCOLOR)
       COMPLEX*16 W(18,NWAVEFUNCS)
       COMPLEX*16 DUM0,DUM1
@@ -115,6 +115,6 @@
           DO J = 1, NCOLOR
               ZTEMP = ZTEMP + CF(J,I)*JAMP(J)
           ENDDO
-          MATRIX%(proc_id)s = MATRIX%(proc_id)s+ZTEMP*DCONJG(JAMP(I))/DENOM(I)   
+          MATRIX%(proc_id)s = MATRIX%(proc_id)s+ZTEMP*DCONJG(JAMP(I))   
       ENDDO
       END

=== modified file 'madgraph/iolibs/template_files/matrix_standalone_matchbox.inc'
--- madgraph/iolibs/template_files/matrix_standalone_matchbox.inc	2014-11-04 08:17:57 +0000
+++ madgraph/iolibs/template_files/matrix_standalone_matchbox.inc	2021-01-30 10:13:29 +0000
@@ -89,7 +89,7 @@
 C  
       INTEGER I,J
       COMPLEX*16 ZTEMP
-      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
+      REAL*8  CF(NCOLOR,NCOLOR)
       COMPLEX*16 AMP(NGRAPHS), JAMP(NCOLOR), LNJAMP(NCOLOR)
       COMPLEX*16 W(18,NWAVEFUNCS)
       COMPLEX*16 DUM0,DUM1
@@ -115,7 +115,7 @@
           DO J = 1, NCOLOR
               ZTEMP = ZTEMP + CF(J,I)*JAMP(J)
           ENDDO
-          %(proc_prefix)sMATRIX = %(proc_prefix)sMATRIX+ZTEMP*DCONJG(JAMP(I))/DENOM(I)   
+          %(proc_prefix)sMATRIX = %(proc_prefix)sMATRIX+ZTEMP*DCONJG(JAMP(I))   
       ENDDO
 %(amp2_lines)s
       END
@@ -162,7 +162,7 @@
 C  
       INTEGER I,J
       COMPLEX*16 ZTEMP
-      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
+      REAL*8 CF(NCOLOR,NCOLOR)
       COMPLEX*16 AMP(NGRAPHS), JAMP(NCOLOR), LNJAMP(NCOLOR)
       COMMON/%(proc_prefix)sJAMP/JAMP,LNJAMP
 

=== modified file 'madgraph/iolibs/template_files/matrix_standalone_matchbox_splitOrders_v4.inc'
--- madgraph/iolibs/template_files/matrix_standalone_matchbox_splitOrders_v4.inc	2014-11-19 10:57:55 +0000
+++ madgraph/iolibs/template_files/matrix_standalone_matchbox_splitOrders_v4.inc	2021-01-30 10:13:29 +0000
@@ -211,7 +211,7 @@
 C  
       INTEGER I,J,M,N
       COMPLEX*16 ZTEMP
-      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
+      REAL*8 CF(NCOLOR,NCOLOR)
       COMPLEX*16 AMP(NGRAPHS)
       COMPLEX*16 JAMP(NCOLOR,NAMPSO), LNJAMP(NCOLOR,NAMPSO)
       COMMON/%(proc_prefix)sJAMP/JAMP,LNJAMP  
@@ -244,7 +244,7 @@
               ZTEMP = ZTEMP + CF(J,I)*JAMP(J,M)
           ENDDO
 		  DO N = 1, NAMPSO
-            RES(%(proc_prefix)sSQSOINDEX(M,N)) = RES(%(proc_prefix)sSQSOINDEX(M,N)) + ZTEMP*DCONJG(JAMP(I,N))/DENOM(I)
+            RES(%(proc_prefix)sSQSOINDEX(M,N)) = RES(%(proc_prefix)sSQSOINDEX(M,N)) + ZTEMP*DCONJG(JAMP(I,N))
 		  ENDDO
         ENDDO
 	  ENDDO
@@ -319,7 +319,7 @@
 C  
       INTEGER I,J,M,N
       COMPLEX*16 ZTEMP
-      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
+      REAL*8 CF(NCOLOR,NCOLOR)
       COMPLEX*16 AMP(NGRAPHS)
       COMPLEX*16 JAMP(NCOLOR,NAMPSO), LNJAMP(NCOLOR,NAMPSO)
       COMMON/%(proc_prefix)sJAMP/JAMP,LNJAMP

=== modified file 'madgraph/iolibs/template_files/matrix_standalone_msF_v4.inc'
--- madgraph/iolibs/template_files/matrix_standalone_msF_v4.inc	2014-02-18 16:35:27 +0000
+++ madgraph/iolibs/template_files/matrix_standalone_msF_v4.inc	2021-01-30 10:13:29 +0000
@@ -94,7 +94,7 @@
 C  
       INTEGER I,J
       COMPLEX*16 ZTEMP
-      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
+      REAL*8 CF(NCOLOR,NCOLOR)
       COMPLEX*16 AMP(NGRAPHS), JAMP(NCOLOR)
       COMPLEX*16 W(18,NWAVEFUNCS)
       COMPLEX*16 DUM0,DUM1
@@ -119,7 +119,7 @@
           DO J = 1, NCOLOR
               ZTEMP = ZTEMP + CF(J,I)*JAMP(J)
           ENDDO
-          MATRIX = MATRIX+ZTEMP*DCONJG(JAMP(I))/DENOM(I)   
+          MATRIX = MATRIX+ZTEMP*DCONJG(JAMP(I))   
       ENDDO
 
       END

=== modified file 'madgraph/iolibs/template_files/matrix_standalone_msP_v4.inc'
--- madgraph/iolibs/template_files/matrix_standalone_msP_v4.inc	2013-12-13 01:58:39 +0000
+++ madgraph/iolibs/template_files/matrix_standalone_msP_v4.inc	2021-01-30 10:13:29 +0000
@@ -89,7 +89,7 @@
 C  
       INTEGER I,J
       COMPLEX*16 ZTEMP
-      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
+      REAL*8 CF(NCOLOR,NCOLOR)
       COMPLEX*16 AMP(NGRAPHS), JAMP(NCOLOR)
       COMPLEX*16 W(18,NWAVEFUNCS)
       COMPLEX*16 DUM0,DUM1
@@ -116,7 +116,7 @@
           DO J = 1, NCOLOR
               ZTEMP = ZTEMP + CF(J,I)*JAMP(J)
           ENDDO
-          MATRIX_PROD = MATRIX_PROD+ZTEMP*DCONJG(JAMP(I))/DENOM(I)   
+          MATRIX_PROD = MATRIX_PROD+ZTEMP*DCONJG(JAMP(I))   
       ENDDO
 %(amp2_lines)s
       END

=== modified file 'madgraph/iolibs/template_files/matrix_standalone_splitOrders_v4.inc'
--- madgraph/iolibs/template_files/matrix_standalone_splitOrders_v4.inc	2018-05-17 12:03:54 +0000
+++ madgraph/iolibs/template_files/matrix_standalone_splitOrders_v4.inc	2021-01-30 10:13:29 +0000
@@ -254,7 +254,7 @@
 C  
       INTEGER I,J,M,N
       COMPLEX*16 ZTEMP
-      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
+      REAL*8 CF(NCOLOR,NCOLOR)
       COMPLEX*16 AMP(NGRAPHS)
 	  COMPLEX*16 JAMP(NCOLOR,NAMPSO)
       COMPLEX*16 W(20,NWAVEFUNCS)
@@ -286,7 +286,7 @@
               ZTEMP = ZTEMP + CF(J,I)*JAMP(J,M)
           ENDDO
 		  DO N = 1, NAMPSO
-            RES(%(proc_prefix)sSQSOINDEX(M,N)) = RES(%(proc_prefix)sSQSOINDEX(M,N)) + ZTEMP*DCONJG(JAMP(I,N))/DENOM(I)
+            RES(%(proc_prefix)sSQSOINDEX(M,N)) = RES(%(proc_prefix)sSQSOINDEX(M,N)) + ZTEMP*DCONJG(JAMP(I,N))
 		  ENDDO
         ENDDO
 	  ENDDO

=== modified file 'madgraph/iolibs/template_files/matrix_standalone_v4.inc'
--- madgraph/iolibs/template_files/matrix_standalone_v4.inc	2020-11-24 16:20:51 +0000
+++ madgraph/iolibs/template_files/matrix_standalone_v4.inc	2021-01-30 10:13:29 +0000
@@ -194,8 +194,8 @@
 C  
       INTEGER I,J
       COMPLEX*16 ZTEMP
-      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
-      COMPLEX*16 AMP(NGRAPHS), JAMP(NCOLOR)
+      REAL*8 CF(NCOLOR,NCOLOR)
+      COMPLEX*16 AMP(NGRAPHS), JAMP(NCOLOR), TMP_JAMP(%(nb_temp_jamp)i)
       COMPLEX*16 W(20,NWAVEFUNCS)
       COMPLEX*16 DUM0,DUM1
       DATA DUM0, DUM1/(0d0, 0d0), (1d0, 0d0)/
@@ -220,8 +220,8 @@
           DO J = 1, NCOLOR
               ZTEMP = ZTEMP + CF(J,I)*JAMP(J)
           ENDDO
-          %(proc_prefix)sMATRIX = %(proc_prefix)sMATRIX+ZTEMP*DCONJG(JAMP(I))/DENOM(I)   
-      ENDDO
+          %(proc_prefix)sMATRIX = %(proc_prefix)sMATRIX+ZTEMP*DCONJG(JAMP(I))
+	  ENDDO
 %(amp2_lines)s
       END
 

=== modified file 'madgraph/iolibs/template_files/realmatrix_fks.inc'
--- madgraph/iolibs/template_files/realmatrix_fks.inc	2020-01-16 13:45:20 +0000
+++ madgraph/iolibs/template_files/realmatrix_fks.inc	2021-01-30 10:13:29 +0000
@@ -112,8 +112,9 @@
       INTEGER I,J
       INTEGER IC(NEXTERNAL)
       data ic /nexternal*1/
-      REAL*8 DENOM(NCOLOR), CF(NCOLOR,NCOLOR)
+      REAL*8 CF(NCOLOR,NCOLOR)
       COMPLEX*16 ZTEMP, AMP(NGRAPHS), JAMP(NCOLOR), W(%(wavefunctionsize)d,NWAVEFUNCS)
+      COMPLEX*16 TMP_JAMP(%(nb_temp_jamp)i)
 C  
 C COLOR DATA
 C  
@@ -129,6 +130,6 @@
           DO J = 1, NCOLOR
               ZTEMP = ZTEMP + CF(J,I)*JAMP(J)
           ENDDO
-          MATRIX_%(N_me)d = MATRIX_%(N_me)d+ZTEMP*DCONJG(JAMP(I))/DENOM(I)   
+          MATRIX_%(N_me)d = MATRIX_%(N_me)d+ZTEMP*DCONJG(JAMP(I))   
       ENDDO
       END

=== modified file 'madgraph/iolibs/template_files/super_auto_dsig_group_v4.inc'
--- madgraph/iolibs/template_files/super_auto_dsig_group_v4.inc	2020-01-16 07:19:30 +0000
+++ madgraph/iolibs/template_files/super_auto_dsig_group_v4.inc	2021-01-30 10:13:29 +0000
@@ -97,6 +97,8 @@
 	  DATA LAST_ICONF/-1/
 	  COMMON/TO_LAST_ICONF/LAST_ICONF
 
+      logical init_mode
+      common /to_determine_zero_hel/init_mode
 C ----------
 C BEGIN CODE
 C ----------
@@ -180,7 +182,7 @@
 		call DS_set_min_points(10,'grouped_processes')
         DO J=1,SYMCONF(0)
           DO IPROC=1,MAXSPROC
-            IF(CONFSUB(IPROC,SYMCONF(J)).NE.0) THEN
+            IF(init_mode.or.CONFSUB(IPROC,SYMCONF(J)).NE.0) THEN
               DO IMIRROR=1,2
                 IF(IMIRROR.EQ.1.OR.MIRRORPROCS(IPROC))THEN
                   CALL MAP_3_TO_1(J,IPROC,IMIRROR,MAXSPROC,2,LMAPPED)
@@ -201,7 +203,7 @@
       LAST_ICONF=0
       DO J=1,SYMCONF(0)
          DO IPROC=1,MAXSPROC
-            IF(CONFSUB(IPROC,SYMCONF(J)).NE.0) THEN
+            IF(init_mode.or.CONFSUB(IPROC,SYMCONF(J)).NE.0) THEN
                DO IMIRROR=1,2
                   IF(IMIRROR.EQ.1.OR.MIRRORPROCS(IPROC))THEN
 C                 Calculate PDF weight for all subprocesses
@@ -245,7 +247,7 @@
         LAST_ICONF=0
         DO J=1,SYMCONF(0)
           DO I=1,MAXSPROC
-            IF(CONFSUB(I,SYMCONF(J)).NE.0) THEN		  
+            IF(init_mode.or.CONFSUB(I,SYMCONF(J)).NE.0) THEN		  
               DO K=1,2
                 IF(K.EQ.1.OR.MIRRORPROCS(I))THEN			
                   IPROC=I
@@ -264,7 +266,11 @@
                     XBK(2)=XDUM
                     CM_RAP=-CM_RAP
                   ENDIF
-                  SELPROC(K,I,J) = DABS(DSIG*SELPROC(K,I,J))
+		  if(init_mode) then
+		    SELPROC(K,I,J) = 1d0
+		  else
+                    SELPROC(K,I,J) = DABS(DSIG*SELPROC(K,I,J))
+		  endif
                   SUMPROB = SUMPROB + SELPROC(K,I,J) 
 		        ENDIF
               ENDDO
@@ -289,7 +295,7 @@
       TOTWGT=0D0
       DO J=1,SYMCONF(0)
         DO I=1,MAXSPROC
-	  IF(CONFSUB(I,SYMCONF(J)).NE.0) THEN
+	  IF(init_mode.or.CONFSUB(I,SYMCONF(J)).NE.0) THEN
           DO K=1,2
             TOTWGT=TOTWGT+SELPROC(K,I,J)
             IF(R.LT.TOTWGT)THEN
@@ -507,4 +513,8 @@
 %(read_write_good_hel)s
 
 
+     SUBROUTINE PRINT_ZERO_AMP()
 
+     %(print_zero_amp)s
+     return
+     end

=== modified file 'madgraph/madevent/gen_ximprove.py'
--- madgraph/madevent/gen_ximprove.py	2020-09-16 07:01:30 +0000
+++ madgraph/madevent/gen_ximprove.py	2021-01-30 10:13:29 +0000
@@ -29,6 +29,7 @@
 import shutil
 import stat
 import sys
+import six
 from six.moves import range
 from six.moves import zip
 
@@ -44,6 +45,8 @@
     import internal.combine_grid as combine_grid
     import internal.combine_runs as combine_runs
     import internal.lhe_parser as lhe_parser
+    if six.PY3:
+        import internal.hel_recycle as hel_recycle
 else:
     MADEVENT= False
     import madgraph.madevent.sum_html as sum_html
@@ -54,6 +57,8 @@
     import madgraph.madevent.combine_grid as combine_grid
     import madgraph.madevent.combine_runs as combine_runs
     import madgraph.various.lhe_parser as lhe_parser
+    if six.PY3:
+        import madgraph.madevent.hel_recycle as hel_recycle
 
 logger = logging.getLogger('madgraph.madevent.gen_ximprove')
 pjoin = os.path.join
@@ -109,22 +114,193 @@
         #if the user defines it in the run_card:
         if self.run_card['survey_splitting'] != -1:
             self.splitted_grid = self.run_card['survey_splitting']
-        if self.run_card['survey_nchannel_per_job'] != -1:
+        if self.run_card['survey_nchannel_per_job'] != 1 and 'survey_nchannel_per_job' in self.run_card.user_set:
             self.combining_job = self.run_card['survey_nchannel_per_job']        
         elif self.run_card['hard_survey'] > 1:
             self.combining_job = 1
             
-            
         
         self.splitted_Pdir = {}
         self.splitted_for_dir = lambda x,y: self.splitted_grid
         self.combining_job_for_Pdir = lambda x: self.combining_job
         self.lastoffset = {}
     
+    done_warning_zero_coupling = False
+    def get_helicity(self, to_submit=True, clean=True):
+        """launch a single call to madevent to get the list of non zero helicity"""
+    
+        self.subproc = [l.strip() for l in open(pjoin(self.me_dir,'SubProcesses', 
+                                                                 'subproc.mg'))]
+        subproc = self.subproc
+        P_zero_result = []
+        nb_tot_proc = len(subproc)
+        job_list = {}      
+        
+          
+        for nb_proc,subdir in enumerate(subproc):
+            self.cmd.update_status('Compiling for process %s/%s.' % \
+                               (nb_proc+1,nb_tot_proc), level=None)
+
+            subdir = subdir.strip()
+            Pdir = pjoin(self.me_dir, 'SubProcesses',subdir)
+            logger.info('    %s ' % subdir)
+            
+            self.cmd.compile(['madevent_forhel'], cwd=Pdir)
+            if not os.path.exists(pjoin(Pdir, 'madevent_forhel')):
+                raise Exception('Error make madevent_forhel not successful')  
+            
+            if not os.path.exists(pjoin(Pdir, 'Hel')):
+                os.mkdir(pjoin(Pdir, 'Hel'))
+                ff = open(pjoin(Pdir, 'Hel', 'input_app.txt'),'w')
+                ff.write('1000 1 1 \n 0.1 \n 2\n 0\n -1\n 1\n')
+                ff.close()
+            else:
+                try:
+                    os.remove(pjoin(Pdir, 'Hel','results.dat'))
+                except Exception:
+                    pass         
+            # Launch gensym
+            p = misc.Popen(['../madevent_forhel < input_app.txt'], stdout=subprocess.PIPE, 
+                                 stderr=subprocess.STDOUT, cwd=pjoin(Pdir,'Hel'), shell=True)
+            #sym_input = "%(points)d %(iterations)d %(accuracy)f \n" % self.opts
+            (stdout, _) = p.communicate(" ".encode())
+            stdout = stdout.decode('ascii')
+            if os.path.exists(pjoin(self.me_dir,'error')):
+                raise Exception(pjoin(self.me_dir,'error')) 
+                # note a continue is not enough here, we have in top to link
+                # the matrixX_optim.f to matrixX_orig.f to let the code to work
+                # after this error.
+
+            if 'no events passed cuts' in stdout:
+                raise Exception
+
+            all_zamp = set()
+            all_hel = set()
+            zero_gc = list()
+            all_zampperhel = set()
+            all_bad_amps_perhel = set()
+            
+            for line in stdout.splitlines():
+                if 'GC_' in line:
+                    lsplit = line.split()
+                    if float(lsplit[2]) ==0 == float(lsplit[3]):
+                        zero_gc.append(lsplit[0])
+                if 'Matrix Element/Good Helicity:' in line:
+                    all_hel.add(tuple(line.split()[3:5]))
+                if 'Amplitude/ZEROAMP:' in line:
+                    all_zamp.add(tuple(line.split()[1:3]))
+                if 'HEL/ZEROAMP:' in line:
+                    nb_mat, nb_hel, nb_amp = line.split()[1:4]
+                    if (nb_mat, nb_hel) not in all_hel:
+                        continue
+                    if (nb_mat,nb_amp) in all_zamp:
+                        continue
+                    all_zampperhel.add(tuple(line.split()[1:4]))
+
+            if zero_gc and not gensym.done_warning_zero_coupling:
+                gensym.done_warning_zero_coupling = True
+                logger.warning("The optimizer detects that you have coupling evaluated to zero: \n"+\
+                                "%s\n" % (' '.join(zero_gc)) +\
+                               "This will slow down the computation. Please consider using restricted model:\n" +\
+                               "https://answers.launchpad.net/mg5amcnlo/+faq/2312")
+            
+                
+            all_good_hels = collections.defaultdict(list)
+            for me_index, hel in all_hel:
+                all_good_hels[me_index].append(int(hel))                           
+                               
+            #print(all_hel)
+            if self.run_card['hel_zeroamp']:
+                all_bad_amps = collections.defaultdict(list)
+                for me_index, amp in all_zamp:
+                    all_bad_amps[me_index].append(int(amp))
+                    
+                all_bad_amps_perhel = collections.defaultdict(list)
+                for me_index, hel, amp in all_zampperhel:
+                    all_bad_amps_perhel[me_index].append((int(hel),int(amp)))    
+                    
+            elif all_zamp:
+                nb_zero = sum(int(a[1]) for a in all_zamp)
+                if zero_gc:
+                    logger.warning("Those zero couplings lead to %s Feynman diagram evaluated to zero (on 10 PS point),\n" % nb_zero +\
+                                   "This part can optimize if you set the flag  hel_zeroamp to True in the run_card."+\
+                                   "Note that restricted model will be more optimal.")
+                else:
+                    logger.warning("The optimization detected that you have %i zero matrix-element for this SubProcess: %s.\n" % nb_zero +\
+                                   "This part can optimize if you set the flag  hel_zeroamp to True in the run_card.")
+            
+            #check if we need to do something and write associate information"
+            data = [all_hel, all_zamp, all_bad_amps_perhel]
+            if not self.run_card['hel_zeroamp']:
+                data[1] = ''
+            if not self.run_card['hel_filtering']:
+                data[0] = ''
+            data = str(data)
+            if os.path.exists(pjoin(Pdir,'Hel','selection')):
+                old_data = open(pjoin(Pdir,'Hel','selection')).read()
+                if old_data == data:
+                    continue
+                
+            
+            with open(pjoin(Pdir,'Hel','selection'),'w') as fsock:
+                fsock.write(data)        
+                
+        
+            for matrix_file in misc.glob('matrix*orig.f', Pdir):
+    
+                split_file = matrix_file.split('/')
+                me_index = split_file[-1][len('matrix'):-len('_orig.f')]
+
+                basename = split_file[-1].replace('orig', 'optim')
+                split_out = split_file[:-1] + [basename]
+                out_file = pjoin('/', '/'.join(split_out))
+
+                basename = 'template_%s' % split_file[-1].replace("_orig", "")
+                split_templ = split_file[:-1] + [basename]
+                templ_file = pjoin('/', '/'.join(split_templ))
+
+                # Convert to sorted list for reproducibility
+                #good_hels = sorted(list(good_hels))
+                good_hels = [str(x) for x in sorted(all_good_hels[me_index])]
+                if self.run_card['hel_zeroamp']:
+                    
+                    bad_amps = [str(x) for x in sorted(all_bad_amps[me_index])]
+                    bad_amps_perhel = [x for x in sorted(all_bad_amps_perhel[me_index])]
+                else:
+                    bad_amps = [] 
+                    bad_amps_perhel = []
+                if __debug__:
+                    mtext = open(matrix_file).read()
+                    nb_amp = int(re.findall('PARAMETER \(NGRAPHS=(\d+)\)', mtext)[0])
+                    logger.debug('nb_hel: %s zero amp: %s bad_amps_hel: %s/%s', len(good_hels),len(bad_amps),len(bad_amps_perhel), len(good_hels)*nb_amp )
+                if len(good_hels) == 1:
+                    files.cp(matrix_file, matrix_file.replace('orig','optim'))
+                    continue # avoid optimization if onlye one helicity
+                recycler = hel_recycle.HelicityRecycler(good_hels, bad_amps, bad_amps_perhel)
+                # In case of bugs you can play around with these:
+                recycler.hel_filt = self.run_card['hel_filtering']
+                recycler.amp_splt = self.run_card['hel_splitamp']
+                recycler.amp_filt = self.run_card['hel_zeroamp']
+
+                recycler.set_input(matrix_file)
+                recycler.set_output(out_file)
+                recycler.set_template(templ_file)              
+                recycler.generate_output_file()
+                del recycler
+
+            # with misc.chdir():
+            #     pass
+
+            #files.ln(pjoin(Pdir, 'madevent_forhel'), Pdir, name='madevent') ##to be removed
+
+        return {}, P_zero_result
+
+    
     def launch(self, to_submit=True, clean=True):
         """ """
-        
-        self.subproc = [l.strip() for l in open(pjoin(self.me_dir,'SubProcesses', 
+
+        if not hasattr(self, 'subproc'):
+            self.subproc = [l.strip() for l in open(pjoin(self.me_dir,'SubProcesses', 
                                                                  'subproc.mg'))]
         subproc = self.subproc
         
@@ -916,7 +1092,7 @@
             if C.get('axsec') == 0:
                 continue
             if goal_lum/(C.get('luminosity')+1e-99) >= 1 + (self.gen_events_security-1)/2:
-                logger.debug("channel %s is at lum=%s (need to improve by %s) (xsec=%s pb)", C.name,  C.get('luminosity'), goal_lum/(C.get('luminosity')+1e-99), C.get('xsec'))
+                logger.debug("channel %s need to improve by %.2f (xsec=%s pb, iter=%s)", C.name, goal_lum/(C.get('luminosity')+1e-99), C.get('xsec'), int(C.get('maxit')))
                 to_refine.append(C)
             elif C.get('xerr') > max(C.get('axsec'),
               (1/(100*math.sqrt(self.err_goal)))*all_channels[-1].get('axsec')):

=== added file 'madgraph/madevent/hel_recycle.py'
--- madgraph/madevent/hel_recycle.py	1970-01-01 00:00:00 +0000
+++ madgraph/madevent/hel_recycle.py	2021-01-30 10:13:29 +0000
@@ -0,0 +1,924 @@
+#!/usr/bin/env python3
+
+import argparse
+import atexit
+import re
+import collections
+from string import Template
+from copy import copy
+from itertools import product
+from functools import reduce 
+
+try:
+     import madgraph
+except:
+     import internal.misc as misc
+else:
+     import madgraph.various.misc as misc
+import mmap
+try:
+    from tqdm import tqdm
+except ImportError:
+    tqdm = misc.tqdm
+
+    
+def get_num_lines(file_path):
+    fp = open(file_path, 'r+')
+    buf = mmap.mmap(fp.fileno(),0)
+    lines = 0
+    while buf.readline():
+        lines += 1
+    return lines
+
+class DAG:
+
+    def __init__(self):
+        self.graph = {}
+        self.all_wavs = []
+        self.external_wavs = []
+        self.internal_wavs = []
+
+    def store_wav(self, wav, ext_deps=[]):
+        self.all_wavs.append(wav)
+        nature = wav.nature
+        if nature == 'external':
+            self.external_wavs.append(wav)
+        if nature == 'internal':
+            self.internal_wavs.append(wav)
+        for ext in ext_deps:
+            self.add_branch(wav, ext)
+
+    def add_branch(self, node_i, node_f):
+        try:
+            self.graph[node_i].append(node_f)
+        except KeyError:
+            self.graph[node_i] = [node_f]
+
+    def dependencies(self, old_name):
+        deps = [wav for wav in self.all_wavs
+                if wav.old_name == old_name and not wav.dead]
+        return deps
+
+    def kill_old(self, old_name):
+        for wav in self.all_wavs:
+            if wav.old_name == old_name:
+                wav.dead = True
+
+    def old_names(self):
+        return {wav.old_name for wav in self.all_wavs}
+
+    def find_path(self, start, end, path=[]):
+        '''Taken from https://www.python.org/doc/essays/graphs/'''
+
+        path = path + [start]
+        if start == end:
+            return path
+        if start not in self.graph:
+            return None
+        for node in self.graph[start]:
+            if node not in path:
+                newpath = self.find_path(node, end, path)
+                if newpath:
+                    return newpath
+        return None
+
+    def __str__(self):
+        return self.__repr__()
+
+    def __repr__(self):
+        print_str = 'With new names:\n\t'
+        print_str += '\n\t'.join([f'{key} : {item}' for key, item in self.graph.items() ])
+        print_str += '\n\nWith old names:\n\t'
+        print_str += '\n\t'.join([f'{key.old_name} : {[i.old_name for i in item]}' for key, item in self.graph.items() ])
+        return print_str
+
+
+
+class MathsObject:
+    '''Abstract class for wavefunctions and Amplitudes'''
+
+    # Store here which externals the last wav/amp depends on.
+    # This saves us having to call find_path multiple times.
+    ext_deps = None
+
+    def __init__(self, arguments, old_name, nature):
+        self.args = arguments
+        self.old_name = old_name
+        self.nature = nature
+        self.name = None
+        self.dead = False
+        self.nb_used = 0
+        self.linkdag = []
+
+    def set_name(self, *args):
+        self.args[-1] = self.format_name(*args)
+        self.name = self.args[-1]
+
+    def format_name(self, *nums):
+        pass
+
+    @staticmethod
+    def get_deps(line, graph):
+        old_args = get_arguments(line)
+        old_name = old_args[-1].replace(' ','')
+        matches = graph.old_names() & set(old_args)
+        try:
+            matches.remove(old_name)
+        except KeyError:
+            pass
+        old_deps = old_args[0:len(matches)]
+
+        # If we're overwriting a wav clear it from graph
+        graph.kill_old(old_name)
+        return [graph.dependencies(dep) for dep in old_deps]
+
+    @classmethod
+    def good_helicity(cls, wavs, graph, diag_number=None, all_hel=[], bad_hel_amp=[]):
+        exts = graph.external_wavs
+        cls.ext_deps = { i for dep in wavs for i in exts if graph.find_path(dep, i) }
+        this_comb_good = False
+        for comb in External.good_wav_combs:
+            if cls.ext_deps.issubset(set(comb)):
+                this_comb_good = True
+                break
+            
+        if diag_number and this_comb_good and cls.ext_deps:
+            helicity = dict([(a.get_id(), a.hel) for a in cls.ext_deps])
+            this_hel = [helicity[i] for i in range(1, len(helicity)+1)] 
+            hel_number = 1 + all_hel.index(tuple(this_hel))
+            
+            if (hel_number,diag_number) in bad_hel_amp:        
+                this_comb_good = False
+            
+
+            
+        return this_comb_good and cls.ext_deps
+
+    @staticmethod
+    def get_new_args(line, wavs):
+        old_args = get_arguments(line)
+        old_name = old_args[-1].replace(' ','')
+        # Work out if wavs corresponds to an allowed helicity combination
+        this_args = copy(old_args)
+        wav_names = [w.name for w in wavs]
+        this_args[0:len(wavs)] = wav_names
+        # This isnt maximally efficient
+        # Could take the num from wavs that've been deleted in graph
+        return this_args
+
+    @staticmethod
+    def get_number():
+        pass
+
+    @classmethod
+    def get_obj(cls, line, wavs, graph, diag_num = None):
+        old_name = get_arguments(line)[-1].replace(' ','')
+        new_args = cls.get_new_args(line, wavs)
+        num = cls.get_number(wavs, graph)
+        
+        this_obj = cls.call_constructor(new_args, old_name, diag_num)
+        this_obj.set_name(num, diag_num)
+        if this_obj.nature != 'amplitude':
+            graph.store_wav(this_obj, cls.ext_deps)
+        return this_obj
+
+
+    def __str__(self):
+        return self.name
+
+    def __repr__(self):
+        return self.name
+
+class External(MathsObject):
+    '''Class for storing external wavefunctions'''
+
+    good_hel = []
+    nhel_lines = ''
+    num_externals = 0
+    # Could get this from dag but I'm worried about preserving order
+    wavs_same_leg = {}
+    good_wav_combs = []
+
+    def __init__(self, arguments, old_name):
+        super().__init__(arguments, old_name, 'external')
+        self.hel = int(self.args[2])
+        self.hel_ranges = []
+        self.raise_num()
+
+    @classmethod
+    def raise_num(cls):
+        cls.num_externals += 1
+
+    @classmethod
+    def generate_wavfuncs(cls, line, graph):
+        # If graph is passed in Internal it should be done here to so
+        # we can set names
+        old_args = get_arguments(line)
+        old_name = old_args[-1].replace(' ','')
+        graph.kill_old(old_name)
+
+        if 'NHEL' in old_args[2].upper():
+            nhel_index = re.search(r'\(.*?\)', old_args[2]).group()
+            ext_num = int(nhel_index[1:-1]) - 1
+            new_hels = sorted(list(External.hel_ranges[ext_num]), reverse=True)
+            new_hels = [int_to_string(i) for i in new_hels]
+        else:
+            # Spinor must be a scalar so give it hel = 0
+            new_hels = [' 0']
+
+        new_wavfuncs = []
+        for hel in new_hels:
+
+            this_args = copy(old_args)
+            this_args[2] = hel
+
+            this_wavfunc = External(this_args, old_name)
+            this_wavfunc.set_name(len(graph.external_wavs) + len(graph.internal_wavs) +1)
+
+            graph.store_wav(this_wavfunc)
+            new_wavfuncs.append(this_wavfunc)
+
+        cls.wavs_same_leg[ext_num] = new_wavfuncs
+
+        return new_wavfuncs
+
+    @classmethod
+    def get_gwc(cls):
+        num_combs = len(cls.good_hel)
+        gwc = [[] for x in range(num_combs)]
+
+        for n, comb in enumerate(cls.good_hel):
+            for leg, wavs in cls.wavs_same_leg.items():
+                for wav in wavs:
+                    if comb[leg] == wav.hel:
+                        gwc[n].append(wav)
+                        
+        cls.good_wav_combs = gwc
+
+    @staticmethod
+    def format_name(*nums):
+        return f'W(1,{nums[0]})'
+    
+    def get_id(self):
+        """ return the id of the particle under consideration """
+        
+        try:
+           return self.id 
+        except:
+            self.id =  int(re.findall(r'P\(0,(\d+)\)', self.args[0])[0])
+            return self.id
+        
+        
+
+class Internal(MathsObject):
+    '''Class for storing internal wavefunctions'''
+
+    max_wav_num = 0
+    num_internals = 0
+
+    @classmethod
+    def raise_num(cls):
+        cls.num_internals += 1
+
+    @classmethod
+    def generate_wavfuncs(cls, line, graph):
+        deps = cls.get_deps(line, graph)
+        new_wavfuncs = [ cls.get_obj(line, wavs, graph) 
+                         for wavs in product(*deps) 
+                         if cls.good_helicity(wavs, graph) ]
+
+        return new_wavfuncs
+
+
+    # There must be a better way
+    @classmethod
+    def call_constructor(cls, new_args, old_name, diag_num):
+        return Internal(new_args, old_name)
+
+    @classmethod
+    def get_number(cls, *args):
+        num = External.num_externals + Internal.num_internals + 1
+        if cls.max_wav_num < num:
+            cls.max_wav_num = num
+        return num
+
+    def __init__(self, arguments, old_name):
+        super().__init__(arguments, old_name, 'internal')
+        self.raise_num()
+
+
+    @staticmethod
+    def format_name(*nums):
+        return f'W(1,{nums[0]})'
+
+class Amplitude(MathsObject):
+    '''Class for storing Amplitudes'''
+
+    max_amp_num = 0
+
+    def __init__(self, arguments, old_name, diag_num):
+        self.diag_num = diag_num
+        super().__init__(arguments, old_name, 'amplitude')
+
+
+    @staticmethod
+    def format_name(*nums):
+        return f'AMP({nums[0]},{nums[1]})'
+
+    @classmethod
+    def generate_amps(cls, line, graph, all_hel=None, all_bad_hel=[]):
+        old_args = get_arguments(line)
+        old_name = old_args[-1].replace(' ','')
+
+        amp_index = re.search(r'\(.*?\)', old_name).group()
+        diag_num = int(amp_index[1:-1])
+
+        deps = cls.get_deps(line, graph)
+
+        new_amps = [cls.get_obj(line, wavs, graph, diag_num) 
+                        for wavs in product(*deps) 
+                        if cls.good_helicity(wavs, graph, diag_num, all_hel,all_bad_hel)]
+
+        return new_amps
+
+    @classmethod
+    def call_constructor(cls, new_args, old_name, diag_num):
+        return Amplitude(new_args, old_name, diag_num)
+
+    @classmethod
+    def get_number(cls, *args):
+        wavs, graph = args
+        amp_num = -1
+        exts = graph.external_wavs
+        for i, comb in enumerate(External.good_wav_combs):
+            if set(comb) == set(cls.ext_deps):
+                # Offset because Fortran counts from 1
+                amp_num = i + 1
+        if amp_num < 1:
+            set_trace()
+            print('Failed to find amp_num')
+            exit(1)
+        if cls.max_amp_num < amp_num:
+            cls.max_amp_num = amp_num 
+        return amp_num  
+
+class HelicityRecycler():
+    '''Class for recycling helicity'''
+
+    def __init__(self, good_elements, bad_amps=[], bad_amps_perhel=[]):
+
+        External.good_hel = []
+        External.nhel_lines = ''
+        External.num_externals = 0
+        External.wavs_same_leg = {}
+        External.good_wav_combs = []
+
+        Internal.max_wav_num = 0
+        Internal.num_internals = 0
+
+        Amplitude.max_amp_num = 0
+
+        self.good_elements = good_elements
+        self.bad_amps = bad_amps
+        self.bad_amps_perhel = bad_amps_perhel
+
+        # Default file names
+        self.input_file = 'matrix_orig.f'
+        self.output_file = 'matrix_orig.f'
+        self.template_file = 'template_matrix.f'
+        
+        self.template_dict = {}
+        #initialise everything as for zero matrix element
+        self.template_dict['helicity_lines'] = '\n'
+        self.template_dict['helas_calls'] = []
+        self.template_dict['jamp_lines'] = '\n'
+        self.template_dict['amp2_lines'] = '\n'
+        self.template_dict['ncomb'] = '0'  
+        self.template_dict['nwavefuncs'] = '0' 
+
+        self.dag = DAG()
+
+        self.diag_num = 1
+        self.got_gwc = False
+
+        self.procedure_name = self.input_file.split('.')[0].upper()
+        self.procedure_kind = 'FUNCTION'
+
+        self.old_out_name = ''
+        self.loop_var = 'K'
+
+        self.all_hel = []
+        self.hel_filt = True
+
+    def set_input(self, file):
+        if 'born_matrix' in file:
+            print('HelicityRecycler is currently '
+                  f'unable to handle {file}')
+            exit(1)
+        self.procedure_name = file.split('.')[0].upper()
+        self.procedure_kind = 'FUNCTION'
+        self.input_file = file
+
+    def set_output(self, file):
+        self.output_file = file
+
+    def set_template(self, file):
+        self.template_file = file
+
+    def function_call(self, line):
+        # Check a function is called at all
+        if not 'CALL' in line:
+            return None
+
+        # Now check for external spinor
+        ext_calls = ['CALL OXXXXX', 'CALL IXXXXX', 'CALL VXXXXX', 'CALL SXXXXX']
+        if any( call in line for call in ext_calls ):
+            return 'external'
+
+        # Now check for internal
+        # Wont find a internal when no externals have been found...
+        # ... I assume
+        if not self.dag.external_wavs:
+            return None
+
+        # Search for internals by looking for calls to the externals
+        # Maybe I should just get a list of all internals?
+        matches = self.dag.old_names() & set(get_arguments(line))
+        try:
+            matches.remove(get_arguments(line)[-1])
+        except KeyError:
+            pass
+        try:
+            function = (line.split('(', 1)[0]).split()[-1]
+        except IndexError:
+            return None
+        # What if [-1] is garbage? Then I'm relying on needs changing.
+        # Is that OK?
+        if (function.split('_')[-1] != '0'):
+            return 'internal'
+        elif (function.split('_')[-1] == '0'):
+            return 'amplitude'
+        else:
+            print(f'Ahhhh what is going on here?\n{line}')
+            set_trace()
+
+        return None
+
+    # string manipulation
+
+    def add_amp_index(self, matchobj):
+        old_pat = matchobj.group()
+        new_pat = old_pat.replace('AMP(', 'AMP( %s,' % self.loop_var)
+        #misc.sprint(self.loop_var, old_pat)
+        
+        #new_pat = f'{self.loop_var},{old_pat[:-1]}{old_pat[-1]}'
+        return new_pat
+
+    def add_indices(self, line):
+        '''Add loop_var index to amp and output variable. 
+           Also update name of output variable.'''
+        # Doesnt work if the AMP arguments contain brackets
+        new_line = re.sub(r'\WAMP\(.*?\)', self.add_amp_index, line)
+        new_line = re.sub(r'MATRIX\d+', 'TS(K)', new_line)
+        return new_line
+
+    def jamp_finished(self, line):
+        # indent_end = re.compile(fr'{self.jamp_indent}END\W')
+        # m = indent_end.match(line)
+        # if m:
+        #     return True
+        return 'init_mode' in line.lower() 
+        #if f'{self.old_out_name}=0.D0' in line.replace(' ', ''):
+        #    return True
+        #return False
+
+    def get_old_name(self, line):
+        if f'{self.procedure_kind} {self.procedure_name}' in line:
+            if 'SUBROUTINE' == self.procedure_kind:
+                self.old_out_name = get_arguments(line)[-1]
+            if 'FUNCTION' == self.procedure_kind:
+                self.old_out_name = line.split('(')[0].split()[-1]
+
+    def get_amp_stuff(self, line_num, line):
+
+        if 'diagram number' in line:
+            self.amp_calc_started = True
+        # Check if the calculation of this diagram is finished
+        if ('AMP' not in get_arguments(line)[-1]
+                and self.amp_calc_started and list(line)[0] != 'C'):
+            # Check if the calculation of all diagrams is finished
+            if self.function_call(line) not in ['external',
+                                                'internal',
+                                                'amplitude']:
+                self.jamp_started = True
+            self.amp_calc_started = False
+        if self.jamp_started:
+            self.get_jamp_lines(line)
+        if self.in_amp2:
+            self.get_amp2_lines(line)
+        if self.find_amp2 and line.startswith('      ENDDO'):
+            self.in_amp2 = True
+            self.find_amp2 = False
+
+    def get_jamp_lines(self, line):
+        if self.jamp_finished(line):
+            self.jamp_started = False
+            self.find_amp2 = True
+        elif not line.isspace():
+            self.template_dict['jamp_lines'] += f'{line[0:6]}  {self.add_indices(line[6:])}'
+
+    def get_amp2_lines(self, line):
+        if line.startswith('      DO I = 1, NCOLOR'):
+            self.in_amp2 = False
+        elif not line.isspace():
+            self.template_dict['amp2_lines'] += f'{line[0:6]}  {self.add_indices(line[6:])}'
+
+    def prepare_bools(self):
+        self.amp_calc_started = False
+        self.jamp_started = False
+        self.find_amp2 = False
+        self.in_amp2 = False
+        self.nhel_started = False
+
+    def unfold_helicities(self, line, nature):
+
+
+
+        #print('deps',line, deps)
+        if nature not in  ['external', 'internal', 'amplitude']:
+            raise Exception('wrong unfolding')
+        
+        if nature == 'external':
+            new_objs = External.generate_wavfuncs(line, self.dag)
+            for obj in new_objs:
+                obj.line = apply_args(line, [obj.args])
+        else:
+            deps = Amplitude.get_deps(line, self.dag)
+            name2dep = dict([(d.name,d) for d in sum(deps,[])])
+            
+            
+        if nature == 'internal':
+            new_objs = Internal.generate_wavfuncs(line, self.dag)
+            for obj in new_objs:
+                obj.line = apply_args(line, [obj.args])
+                obj.linkdag = []
+                for name in obj.args:
+                    if name in name2dep:
+                        name2dep[name].nb_used +=1
+                        obj.linkdag.append(name2dep[name])
+                
+        if nature == 'amplitude':
+            nb_diag = re.findall(r'AMP\((\d+)\)', line)[0]
+            if nb_diag not in self.bad_amps:
+                new_objs = Amplitude.generate_amps(line, self.dag, self.all_hel, self.bad_amps_perhel)
+                out_line = self.apply_amps(line, new_objs)
+                for i,obj in enumerate(new_objs):
+                    if i == 0: 
+                        obj.line = out_line
+                        obj.nb_used = 1
+                    else:
+                        obj.line = ''
+                        obj.nb_used = 1
+                    obj.linkdag = []
+                    for name in obj.args:
+                        if name in name2dep:
+                            name2dep[name].nb_used +=1
+                            obj.linkdag.append(name2dep[name])
+            else:
+                return ''
+
+          
+        return new_objs
+        #return f'{line}\n' if nature == 'external' else line
+
+    def apply_amps(self, line, new_objs):
+        if self.amp_splt:
+            return split_amps(line, new_objs)  
+        else: 
+
+            return apply_args(line, [i.args for i in new_objs])
+
+    def get_gwc(self, line):
+        if self.got_gwc:
+            return
+        
+        num_found = len({obj.old_name for obj in self.template_dict['helas_calls']
+                         if isinstance(obj, External)})
+            
+        try:
+            num_exts = len(External.good_hel[0])
+        except IndexError:
+            return
+        if num_found <= num_exts:
+            External.get_gwc()
+        if num_found == num_exts :
+            self.got_gwc=True
+
+    def get_good_hel(self, line):
+        if 'DATA (NHEL' in line:
+            self.nhel_started = True
+            this_hel = [int(hel) for hel in line.split('/')[1].split(',')]
+            self.all_hel.append(tuple(this_hel))
+        elif self.nhel_started:
+            self.nhel_started = False
+            
+            if self.hel_filt:
+                External.good_hel = [ self.all_hel[int(i)-1] for i in self.good_elements ]
+            else:
+                External.good_hel = self.all_hel
+
+            External.hel_ranges = [set() for hel in External.good_hel[0]]
+            for comb in External.good_hel:
+                for i, hel in enumerate(comb):
+                    External.hel_ranges[i].add(hel)
+
+            self.counter = 0
+            nhel_array = [self.nhel_string(hel)
+                          for hel in External.good_hel]
+            nhel_lines = '\n'.join(nhel_array)
+            self.template_dict['helicity_lines'] += nhel_lines
+
+            self.template_dict['ncomb'] = len(External.good_hel)
+
+    def nhel_string(self, hel_comb):
+        self.counter += 1
+        formatted_hel = [f'{hel}' if hel < 0 else f' {hel}' for hel in hel_comb]
+        nexternal = len(hel_comb)
+        return (f'      DATA (NHEL(I,{self.counter}),I=1,{nexternal}) /{",".join(formatted_hel)}/')
+
+    def read_orig(self):
+
+        with open(self.input_file, 'r') as input_file:
+
+            self.prepare_bools()
+
+            for line_num, line in tqdm(enumerate(input_file), total=get_num_lines(self.input_file)):
+                if line_num == 0:
+                    line_cache = line
+                    continue
+                
+                if '!SKIP' in line:
+                    continue
+                
+                char_5 = ''
+                try:
+                    char_5 = line[5]
+                except IndexError:
+                    pass
+                if char_5 == '$':
+                    line_cache = undo_multiline(line_cache, line)
+                    continue
+
+                line, line_cache = line_cache, line
+
+                self.get_old_name(line)
+                self.get_good_hel(line)
+                self.get_amp_stuff(line_num, line)
+                self.get_gwc(line)
+
+                call_type = self.function_call(line)
+                if call_type in ['external', 'internal', 'amplitude']:
+                    self.template_dict['helas_calls'] += self.unfold_helicities(
+                        line, call_type)
+
+        self.template_dict['nwavefuncs'] = max(External.num_externals, Internal.max_wav_num)
+        # filter out uselless call
+        for i in range(len(self.template_dict['helas_calls'])-1,-1,-1):
+            obj = self.template_dict['helas_calls'][i]
+            if obj.nb_used == 0:
+                obj.line = ''
+                for dep in obj.linkdag:
+                    dep.nb_used -= 1
+
+        
+        
+        self.template_dict['helas_calls'] = '\n'.join([f'{obj.line.rstrip()} ! count {obj.nb_used}' 
+                                 for obj in self.template_dict['helas_calls']
+                                 if obj.nb_used > 0 and obj.line])
+
+    def read_template(self):
+        out_file = open(self.output_file, 'w+')
+        with open(self.template_file, 'r') as file:
+            for line in file:
+                s = Template(line)
+                line = s.safe_substitute(self.template_dict)
+                line = '\n'.join([do_multiline(sub_lines) for sub_lines in line.split('\n')])
+                out_file.write(line)
+        out_file.close()
+
+    def write_zero_matrix_element(self):
+        out_file = open(self.output_file, 'w+')
+        self.template_dict['ncomb'] = '0'  
+        self.template_dict['nwavefuncs'] = '0'  
+        self.template_dict['helas_calls'] = ''    
+        with open(self.template_file, 'r') as file:
+            for line in file:
+                s = Template(line)
+                line = s.safe_substitute(self.template_dict)
+                line = '\n'.join([do_multiline(sub_lines) for sub_lines in line.split('\n')])
+                out_file.write(line)
+        out_file.close()
+
+
+    def generate_output_file(self):
+        if not self.good_elements:
+            misc.sprint("No helicity", self.input_file)
+            self.write_zero_matrix_element()
+            return
+        
+        atexit.register(self.clean_up)
+        self.read_orig()
+        self.read_template()
+        atexit.unregister(self.clean_up)
+
+    def clean_up(self):
+        pass
+
+
+def get_arguments(line):
+    '''Find the substrings separated by commas between the first
+    closed set of parentheses in 'line'. 
+    '''
+    bracket_depth = 0
+    element = 0
+    arguments = ['']
+    for char in line:
+        if char == '(':
+            bracket_depth += 1
+            if bracket_depth - 1 == 0:
+                # This is the first '('. We don't want to add it to
+                # 'arguments'
+                continue
+        if char == ')':
+            bracket_depth -= 1
+            if bracket_depth == 0:
+                # We've reached the end
+                break
+        if char == ',' and bracket_depth == 1:
+            element += 1
+            arguments.append('')
+            continue
+        if bracket_depth > 0:
+            arguments[element] += char
+    return arguments
+
+
+def apply_args(old_line, all_the_args):
+    function = (old_line.split('(')[0]).split()[-1]
+    old_args = old_line.split(function)[-1]
+    new_lines = [old_line.replace(old_args, f'({",".join(x)})\n')
+                 for x in all_the_args]
+    
+    return ''.join(new_lines)
+
+def split_amps(line, new_amps):
+    fct = line.split('(',1)[0].split('_0')[0]
+    for i,amp in enumerate(new_amps):
+        if i == 0:
+            occur = []
+            for a in amp.args:
+                if "W(1," in a:
+                    tmp = collections.defaultdict(int)
+                    tmp[a] += 1
+                    occur.append(tmp)
+        else:
+            for i in range(len(occur)):
+                a = amp.args[i]
+                occur[i][a] +=1
+    # Each element in occur is the wavs that appear in a column, with
+    # the number of occurences
+    nb_wav =  [len(o) for o in occur]
+    to_remove = nb_wav.index(max(nb_wav)) 
+    # Remove the one that occurs the most
+    occur.pop(to_remove)
+    
+    lines = [] 
+    # Get the wavs per column
+    wav_name = [o.keys() for o in occur]          
+    for wfcts in product(*wav_name):
+        # Select the amplitudes produced by wfcts
+        sub_amps = [amp for amp in new_amps 
+                    if all(w in amp.args for w in wfcts)]
+        if not sub_amps:
+            continue
+        if len(sub_amps) ==1:
+            lines.append(apply_args(line, [i.args for i in sub_amps]).replace('\n',''))
+            
+            continue
+                         
+        # the next line is to make the code nicer 
+        sub_amps.sort(key=lambda a: int(a.args[-1][:-1].split(',',1)[1]))
+        windices = []
+        hel_calculated = []
+        iamp = 0
+        for i,amp in enumerate(sub_amps):
+            args = amp.args[:]   
+            # Remove wav and get its index
+            wcontract = args.pop(to_remove)
+            windex = wcontract.split(',')[1].split(')')[0]
+            windices.append(windex)
+            amp_result,  args[-1]  =  args[-1], 'TMP(1)'
+            
+            if i ==0:
+                # Call the original fct with P1N_...
+                # Final arg is replaced with TMP(1)
+                spin = fct.split(None,1)[1][to_remove]
+                lines.append('%sP1N_%s(%s)' % (fct, to_remove+1, ', '.join(args)))
+
+            hel, iamp = re.findall('AMP\((\d+),(\d+)\)', amp_result)[0]
+            hel_calculated.append(hel)
+            #lines.append(' %(result)s = TMP(3) * W(3,%(w)s) + TMP(4) * W(4,%(w)s)+'
+            #             % {'result': amp_result, 'w':  windex}) 
+            #lines.append('     &             TMP(5) * W(5,%(w)s)+TMP(6) * W(6,%(w)s)'
+            #             % {'result': amp_result, 'w':  windex})
+        if spin in "VF":
+            lines.append("""      call CombineAmp(%(nb)i,
+     & (/%(hel_list)s/), 
+     & (/%(w_list)s/),
+     & TMP, W, AMP(1,%(iamp)s))""" %
+                               {'nb': len(sub_amps),
+                                'hel_list': ','.join(hel_calculated),
+                                'w_list': ','.join(windices),
+                                'iamp': iamp
+                               })
+        elif spin == "S":
+            lines.append("""      call CombineAmpS(%(nb)i, 
+     &(/%(hel_list)s/), 
+     & (/%(w_list)s/), 
+     & TMP, W, AMP(1,%(iamp)s))""" %
+                               {'nb': len(sub_amps),
+                                'hel_list': ','.join(hel_calculated),
+                                'w_list': ','.join(windices),
+                                'iamp': iamp
+                               })            
+        else:
+            raise Exception("split amp are not supported for spin2 and 3/2")
+            
+    #lines.append('')
+    return '\n'.join(lines)
+
+def get_num(wav):
+    name = wav.name
+    between_brackets = re.search(r'\(.*?\)', name).group()
+    num = int(between_brackets[1:-1].split(',')[-1])