Merge lp:~maddevelopers/mg5amcnlo/phase_space_stuff into lp:~maddevelopers/mg5amcnlo/2.5.1

[Rikkert Frederix]

Slight improvement in the phase-space generation for NLO processes. Also, a bit of factorising of that code to make it slightly simpler.

[marco zaro (marco-zaro)]

Hi Rik,
I am looking at this merge, and I have a bunch of questions. Are your changes supposed to improve the performances of the code? Do the changes include the resonance-aware subtraction?
Thanks,
Marco

review: Needs Information

[Rikkert Frederix (frederix)]

Hi Marco,

It's just a phase-space optimisation. It doesn't include the resonance aware stuff. But for processes with competing resonances it is working slightly more efficient than the old way of doing things. It's not a great improvement, though.

Cheers,
Rik

[marco zaro (marco-zaro)]

Ciao Rik,
i have checked g g > h > 4 leptons in the heft model and the results are the same as with the stanadrd branch (within uncertainties). Timings are also similar, i think you can do the merge.

cheers,
Marco

review: Approve

=== modified file 'Template/NLO/SubProcesses/genps_fks.f'
--- Template/NLO/SubProcesses/genps_fks.f	2015-10-07 18:25:55 +0000
+++ Template/NLO/SubProcesses/genps_fks.f	2016-09-30 09:38:53 +0000
@@ -168,8 +168,8 @@
      &     ,ionebody,fksmass,nbranch
 c Conflicting BW stuff
       integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
-      double precision cBW_mass(-nexternal:-1,-1:1),
-     &     cBW_width(-nexternal:-1,-1:1)
+      double precision cBW_mass(-1:1,-nexternal:-1),
+     &     cBW_width(-1:1,-nexternal:-1)
       common/c_conflictingBW/cBW_mass,cBW_width,cBW_level_max,cBW
      $     ,cBW_level
 
@@ -251,13 +251,13 @@
             if(nt_channel.eq.0 .and. qwidth(-ns_channel-1).ne.0.d0 .and.
      $           cBW(-ns_channel-1).ne.2)then
 c Generate tau according to a Breit-Wiger function
-               call generate_tau_BW(stot,x(ndim-4),qmass(-ns_channel-1)
-     $              ,qwidth(-ns_channel-1),cBW(-ns_channel-1),cBW_mass(
-     $              -ns_channel-1,1),cBW_width(-ns_channel-1,1),tau_born
-     $              ,xjac0)
+               call generate_tau_BW(stot,ndim-4,x(ndim-4),qmass(
+     $              -ns_channel-1),qwidth(-ns_channel-1),cBW(-ns_channel
+     $              -1),cBW_mass(-1, -ns_channel-1),cBW_width(-1,
+     $              -ns_channel-1),tau_born,xjac0)
             else
-c not a Breit Wigner
-               call generate_tau(x(ndim-4),tau_born,xjac0)
+c     not a Breit Wigner
+               call generate_tau(stot,ndim-4,x(ndim-4),tau_born,xjac0)
             endif
          endif
 c Generate the rapditity of the Born system
@@ -313,7 +313,8 @@
 c
 c Start by generating all the invariant masses of the s-channels
       call generate_inv_mass_sch(ns_channel,itree,m,sqrtshat_born
-     &     ,totmass,qwidth,qmass,cBW,cBW_mass,cBW_width,s,x,xjac0,pass)
+     $     ,totmass,qwidth,qmass,cBW,cBW_level,cBW_mass,cBW_width,s,x
+     $     ,xjac0,pass)
       if (.not.pass) goto 222
 c If only s-channels, also set the p1+p2 s-channel
       if (nt_channel .eq. 0 .and. nincoming .eq. 2) then
@@ -328,8 +329,8 @@
 c     Next do the T-channel branchings
 c
       if (nt_channel.ne.0) then
-         call generate_t_channel_branchings(ns_channel,nbranch,itree,m,s
-     &        ,x,pb,xjac0,xpswgt0,pass)
+         call generate_t_channel_branchings(ns_channel,nbranch,itree
+     $        ,m,s,x,pb,xjac0,xpswgt0,pass)
          if (.not.pass) goto 222
       endif
 c
@@ -1107,8 +1108,8 @@
          costh_i_fks=y_ij_fks+(1-y_ij_fks**2)*xi_i_fks/sqrt(cffC2)
          if(abs(costh_i_fks).gt.1.d0)costh_i_fks=y_ij_fks
       else
-         costh_i_fks=(x3len_fks_mother**2-x3len_j_fks**2+x3len_i_fks**2)/
-     &               (2*x3len_fks_mother*x3len_i_fks)
+         costh_i_fks=(x3len_fks_mother**2-x3len_j_fks**2+x3len_i_fks**2)
+     $        /(2*x3len_fks_mother*x3len_i_fks)
          if(abs(costh_i_fks).gt.1.d0+qtiny)then
             write(*,*)'Fatal error #8 in one_tree',
      &           costh_i_fks,xi_i_fks,y_ij_fks,xmrec2
@@ -1425,7 +1426,6 @@
 c insert here further importance sampling towards xi_i_hat->0
             xi_i_hat=sstiny+(1-sstiny)*x(1)**2
          endif
-c$$$         xi_i_fks=xi_i_hat*xiimax
          xi_i_fks=xiimin+(xiimax-xiimin)*xi_i_hat
       elseif( (icountevts.eq.-100.or.abs(icountevts).eq.1) .and.
      &        (colltest.and.xi_i_fks_fix.ne.-2.d0) .and.
@@ -1771,8 +1771,8 @@
       end
 
 
-      subroutine generate_tau_BW(stot,x,mass,width,cBW,BWmass,BWwidth
-     $     ,tau,jac)
+      subroutine generate_tau_BW(stot,idim,x,mass,width,cBW,BWmass
+     $     ,BWwidth,tau,jac)
       implicit none
       real*8 pi
       parameter (pi=3.1415926535897932d0)
@@ -1780,9 +1780,10 @@
       parameter (nsamp=1)
       include 'run.inc'
       include 'genps.inc'
-      integer cBW,icount
+      integer cBW,icount,idim
       data icount /0/
-      double precision stot,x,tau,jac,mass,width,BWmass,BWwidth,m,w,a,b
+      double precision stot,x,tau,jac,mass,width,m,w,a,b,BWmass(-1:1)
+     $     ,BWwidth(-1:1),s,s_mass
       double precision smax,smin,xm02,stemp,bwmdpl,bwmdmn,bwfmpl,bwfmmn
       double precision bwdelf,xbwmass3,bwfunc,x0
       external xmwmass3,bwfunc
@@ -1792,154 +1793,54 @@
      &     ,tau_lower_bound
       common/ctau_lower_bound/tau_Born_lower_bound
      &     ,tau_lower_bound_resonance,tau_lower_bound
-      if (cBW.eq.1 .and. width.gt.0d0 .and. BWwidth.gt.0d0) then
-c conflicting Breit-Wigner
-c     use the ratio of the widths to determine which is the most-likely
-c     one to go on-shell
-         a=1d0-width/(width+BWwidth)
-c     use the mass separation to determine which region of phase-space
-c     should be generated by which BW
-         if (BWmass.le.mass) then
-            write (*,*) 'Error in generate_tau_BW #1',mass,BWmass
-            stop
-         endif
-         b=(BWmass-mass)/(width+BWwidth)
-         b=mass+b*width
-         b=b**2
-         if (x.lt.0.5d0) then
-c the resonace BW
-            m=mass
-            w=width
-            x0=x*2d0
-            jac=jac*2d0
-            smax=b
-            smin=tau_Born_lower_bound*stot
-         else
-c the alternative "BW". It's a rather flat distribution peaked above
-c BWmass: we can use the normal generate_tau with the
-c 'tau_lower_bound_resonance' equal to the BW mass, and
-c 'tau_lower_bound' equal to the parameter 'b' (with the correct
-c normalization)
-            x0=2d0*x-1d0
-            jac=jac*2d0
-c     save default tau's
-            tau_lower_bound_resonance_save=tau_lower_bound_resonance
-            tau_Born_lower_bound_save=tau_Born_lower_bound
-            tau_lower_bound_save=tau_lower_bound
-c     overwrite the tau's
-            tau_Born_lower_bound=b/stot
-            tau_lower_bound=tau_Born_lower_bound
-            tau_lower_bound_resonance=BWmass**2/stot
-c     Generate tau
-            call generate_tau(x0,tau,jac)
-c     restore default tau's
-            tau_lower_bound_resonance=tau_lower_bound_resonance_save
-            tau_Born_lower_bound=tau_Born_lower_bound_save
-            tau_lower_bound=tau_lower_bound_save
-            return      ! tau and jac generate: exit this function
-         endif
-      elseif (cBW.eq.0) then
-c Normal Breit-Wigner
-         m=mass
-         w=width
-         x0=x
+      if (cBW.eq.1 .and. width.gt.0d0 .and. BWwidth(1).gt.0d0) then
+         smin=tau_Born_lower_bound*stot
          smax=stot
-         smin=tau_Born_lower_bound*stot
+         s_mass=smin
+         call trans_x(5,idim,x,smin,smax,s_mass,mass,width,BWmass(
+     $        -1),BWwidth(-1),jac,s)
+         tau=s/stot
+         jac=jac/stot
       else
-         write (*,*) 'Error in generate_tau_BW #2',cBW
-         stop
+         smin=tau_Born_lower_bound*stot
+         smax=stot
+         s_mass=smin
+         call trans_x(3,idim,x,smin,smax,s_mass,mass,width,BWmass(
+     $        -1),BWwidth(-1),jac,s)
+         tau=s/stot
+         jac=jac/stot
       endif
-      xm02=m**2
-      bwmdpl=smax-xm02
-      bwmdmn=xm02-smin
-      bwfmpl=atan(bwmdpl/(m*w))
-      bwfmmn=atan(bwmdmn/(m*w))
-      bwdelf=(bwfmpl+bwfmmn)/pi
-      stemp=xbwmass3(x0,xm02,w,bwdelf,bwfmmn)
-      jac=jac*bwdelf/bwfunc(stemp,xm02,w)
-      tau=stemp/stot
-      jac=jac/stot
       return
       end
 
 
-      subroutine generate_tau(x,tau,jac)
+      subroutine generate_tau(stot,idim,x,tau,jac)
       double precision x,tau,jac
       double precision roH,roHs,fract,ximax0,ximin0,tmp,fract1,fract2
-     &     ,roHj
-      integer nsamp
+     &     ,roHj,stot,s,dum,dum3(-1:1),smin,smax,s_mass
+      integer nsamp,idim
       parameter (nsamp=1)
       double precision tau_Born_lower_bound,tau_lower_bound_resonance
-     &     ,tau_lower_bound
+     &     ,tau_lower_bound,tiny
+      parameter (tiny=1d-8)
       common/ctau_lower_bound/tau_Born_lower_bound
      &     ,tau_lower_bound_resonance,tau_lower_bound
       character*4 abrv
       common /to_abrv/ abrv
-c The only possible order for the lower bounds on tau is as follows:
-c tau_Born_lower_bound <= tau_lower_bound <= tau_lower_bound_resonance
-c If the order is different, we should quit.
-      if ( tau_Born_lower_bound.ge.tau_lower_bound+1d8 .or.
-     &     tau_Born_lower_bound.ge.tau_lower_bound_resonance+1d8 .or.
-     &     tau_lower_bound.ge.tau_lower_bound_resonance+1d8 ) then
-         write (*,*) 'Errors on the tau_Born_lower_bound'
-     &        ,tau_Born_lower_bound,tau_lower_bound
-     &        ,tau_lower_bound_resonance
-         stop
-      endif
-      roH=tau_Born_lower_bound
-      roHj=tau_lower_bound
-      roHs=tau_lower_bound_resonance
-c User x below 'fract1' for phase-space region below jet cut-off
-      if (tau_lower_bound.le.tau_Born_lower_bound +1d-8) then
-         fract1=0d0
-      else
-         fract1=0.05d0
-      endif
-c User x between 'fract1' and 'fract2' for phase-space region between
-c jet and soft cut-off
-      if(tau_lower_bound_resonance.le.tau_lower_bound + 1d-8)then
-         fract2=0.0d0
-      else
-         fract2=0.05d0
-      endif
-      if (x.lt.fract1) then
-c Use x^2 importance sampling below jet cut-off
-         if (fract1.lt.1d-5) then
-            write (*,*) 'fract1 is too small'
-            stop
-         endif
-         tau=roHj-(roHj-roH)*(1d0-x/fract1)**2
-         if (tau.le.1d-8) tau=1d-8 ! avoid numerical pathologies
-         jac=jac*(roHj-roH)*2d0*(1d0-x/fract1)/fract1
-         if (abrv.eq.'grid') then
-c If we are setting-up the integration grids, using 'grid' (i.e. only
-c the Born), we have to make sure that we have a viable PS point also
-c when tau_born < tau_lower_bound. Simply set tau equal to the lower
-c bound. Don't worry about putting a correct Jacobian, because the
-c 'grid' cross section is non-physical anyway. We can use the Jacobian
-c to enhance a bit the region close to the cut (and compensate for the
-c Jacobian above, which goes in the other direction).
-            tau=roHj+1d-8
-            jac=jac*(x/fract1)**2
-         endif
-      elseif (x.lt.fract1+fract2) then
-c Flat grid below soft cut-off
-         if (fract2.lt.1d-5) then
-            write (*,*) 'fract2 is too small'
-            stop
-         endif
-         tau=roHj-(roHs-roHj)*(fract1-x)/fract2
-         jac=jac*(roHs-roHj)/fract2
-      else
-c Use 1/x^(nsamp) importance sampling above soft cut-off
-         fract=fract1+fract2
-         ximax0 = roHs**(-nsamp)
-         ximin0 = 1.d0
-         tmp  = ximin0 +(1d0-(x-fract)/(1d0-fract))*(ximax0-ximin0)
-         tau = tmp**(-1/dble(nsamp))
-         jac= jac/nsamp*tau**(nsamp+1)*
-     &        (ximax0-ximin0)/(1d0-fract)
-      endif
+      smin=tau_born_lower_bound*stot
+      smax=stot
+      s_mass=tau_lower_bound_resonance*stot
+      if (s_mass.gt.smin+tiny) then
+         call trans_x(2,idim,x,smin,smax,s_mass,dum,dum
+     $        ,dum3,dum3,jac,s)
+      elseif(abs(s_mass-smin).lt.tiny) then
+         call trans_x(7,idim,x,smin,smax,s_mass,dum,dum
+     $        ,dum3,dum3,jac,s)
+      else
+         write (*,*) 'ERROR #39 in genps_fks.f',s_mass,smin
+      endif
+      tau=s/stot
+      jac=jac/stot
       return
       end
 
@@ -1985,7 +1886,8 @@
 
       
       subroutine generate_inv_mass_sch(ns_channel,itree,m,sqrtshat_born
-     &     ,totmass,qwidth,qmass,cBW,cBW_mass,cBW_width,s,x,xjac0,pass)
+     $     ,totmass,qwidth,qmass,cBW,cBW_level,cBW_mass,cBW_width,s,x
+     $     ,xjac0,pass)
       implicit none
       real*8 pi
       parameter (pi=3.1415926535897932d0)
@@ -1997,24 +1899,27 @@
       double precision s(-max_branch:max_particles)
       double precision sqrtshat_born,totmass,xjac0
       integer itree(2,-max_branch:-1)
-      integer i,j,nsamp
-      parameter (nsamp=1)
+      integer i,j,ii,order(-nexternal:0)
       double precision smin,smax,xm02,bwmdpl,bwmdmn,bwfmpl,bwfmmn,bwdelf
      &     ,totalmass,tmp,ximin0,ximax0
       double precision xbwmass3,bwfunc
       external xbwmass3,bwfunc
       logical pass
       integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
-      double precision cBW_mass(-nexternal:-1,-1:1),
-     &     cBW_width(-nexternal:-1,-1:1)
+      double precision cBW_mass(-1:1,-nexternal:-1),
+     &     cBW_width(-1:1,-nexternal:-1)
       double precision b(-1:1),x0
-      double precision s_mass(-nexternal:-1),xi,fract
+      double precision s_mass(-nexternal:nexternal),xi,fract,vol_new
+     $     ,shat,x_new
       parameter (fract=0.1d0)
       common/to_phase_space_s_channel/s_mass
       pass=.true.
       totalmass=totmass
-      do i = -1,-ns_channel,-1
-c Generate invariant masses for all s-channel branchings of the Born
+      do ii = -1,-ns_channel,-1
+c Randomize the order with which to generate the s-channel masses:
+         call sChan_order(ns_channel,order)
+         i=order(ii)
+c     Generate invariant masses for all s-channel branchings of the Born
          smin = (m(itree(1,i))+m(itree(2,i)))**2
          smax = (sqrtshat_born-totalmass+sqrt(smin))**2
          if(smax.lt.smin.or.smax.lt.0.d0.or.smin.lt.0.d0)then
@@ -2026,139 +1931,67 @@
          if(qwidth(i).ne.0.d0 .and. cBW(i).ne.2)then
 c Breit Wigner
             if (cBW(i).eq.1 .and.
-     &          cBW_width(i,1).gt.0d0 .and. cBW_width(i,-1).gt.0d0) then
+     &          cBW_width(1,i).gt.0d0 .and. cBW_width(-1,i).gt.0d0) then
 c     conflicting BW on both sides
-               do j=-1,1,2
-                  b(j)=(cBW_mass(i,j)-qmass(i))/
-     &                 (qwidth(i)+cBW_width(i,j))
-                  b(j)=qmass(i)+b(j)*qwidth(i)
-                  b(j)=b(j)**2
-               enddo
-               if (x(-i).lt.1d0/3d0) then
-                  x0=3d0*x(-i)
-                  s(i)=(b(-1)-smin)*x0+smin
-                  xjac0=3d0*xjac0*(b(-1)-smin)
-               elseif (x(-i).gt.1d0/3d0 .and. x(-i).lt.2d0/3d0) then
-                  x0=3d0*x(-i)-1d0
-                  xm02=qmass(i)**2
-                  bwmdpl=b(1)-xm02
-                  bwmdmn=xm02-b(-1)
-                  bwfmpl=atan(bwmdpl/(qmass(i)*qwidth(i)))
-                  bwfmmn=atan(bwmdmn/(qmass(i)*qwidth(i)))
-                  bwdelf=(bwfmpl+bwfmmn)/pi
-                  s(i)=xbwmass3(x0,xm02,qwidth(i),bwdelf
-     &                 ,bwfmmn)
-                  xjac0=3d0*xjac0*bwdelf/bwfunc(s(i),xm02,qwidth(i))
-               else
-                  x0=3d0*x(-i)-2d0
-                  s(i)=(smax-b(1))*x0+b(1)
-                  xjac0=3d0*xjac0*(smax-b(1))
-               endif
-            elseif (cBW(i).eq.1.and.cBW_width(i,1).gt.0d0) then
+               call trans_x(6,-i,x(-i),smin,smax,s_mass(i),qmass(i)
+     $              ,qwidth(i),cBW_mass(-1,i),cBW_width(-1,i),xjac0
+     $              ,s(i))
+            elseif (cBW(i).eq.1.and.cBW_width(1,i).gt.0d0) then
 c     conflicting BW with alternative mass larger
-               b(1)=(cBW_mass(i,1)-qmass(i))/
-     &              (qwidth(i)+cBW_width(i,1))
-               b(1)=qmass(i)+b(1)*qwidth(i)
-               b(1)=b(1)**2
-               if (x(-i).lt.0.5d0) then
-                  x0=2d0*x(-i)
-                  xm02=qmass(i)**2
-                  bwmdpl=b(1)-xm02
-                  bwmdmn=xm02-smin
-                  bwfmpl=atan(bwmdpl/(qmass(i)*qwidth(i)))
-                  bwfmmn=atan(bwmdmn/(qmass(i)*qwidth(i)))
-                  bwdelf=(bwfmpl+bwfmmn)/pi
-                  s(i)=xbwmass3(x0,xm02,qwidth(i),bwdelf
-     &                 ,bwfmmn)
-                  xjac0=2d0*xjac0*bwdelf/bwfunc(s(i),xm02,qwidth(i))
-               else
-                  x0=2d0*x(-i)-1d0
-                  s(i)=(smax-b(1))*x0+b(1)
-                  xjac0=2d0*xjac0*(smax-b(1))
-               endif
-            elseif (cBW(i).eq.1.and.cBW_width(i,-1).gt.0d0) then
+               call trans_x(5,-i,x(-i),smin,smax,s_mass(i),qmass(i)
+     $              ,qwidth(i),cBW_mass(-1,i),cBW_width(-1,i),xjac0
+     $              ,s(i))
+            elseif (cBW(i).eq.1.and.cBW_width(-1,i).gt.0d0) then
 c     conflicting BW with alternative mass smaller
-               b(-1)=(cBW_mass(i,-1)-qmass(i))/
-     &              (qwidth(i)+cBW_width(i,-1)) ! b(-1) is negative here
-               b(-1)=qmass(i)+b(-1)*qwidth(i)
-               b(-1)=b(-1)**2
-
-               if (b(-1).gt.smax) then
-                   s(i)=(smax-smin)*x(-i)+smin
-                   xjac0=xjac0*(smax-smin)
-               elseif (x(-i).lt.0.5d0) then
-                  x0=2d0*x(-i)
-                  s(i)=(b(-1)-smin)*x0+smin
-                  xjac0=2d0*xjac0*(b(-1)-smin)
-               else
-                  x0=2d0*x(-i)-1d0
-                  xm02=qmass(i)**2
-                  bwmdpl=smax-xm02
-                  bwmdmn=xm02-b(-1)
-                  bwfmpl=atan(bwmdpl/(qmass(i)*qwidth(i)))
-                  bwfmmn=atan(bwmdmn/(qmass(i)*qwidth(i)))
-                  bwdelf=(bwfmpl+bwfmmn)/pi
-                  s(i)=xbwmass3(x0,xm02,qwidth(i),bwdelf
-     &                 ,bwfmmn)
-                  xjac0=2d0*xjac0*bwdelf/bwfunc(s(i),xm02,qwidth(i))
-               endif
+               call trans_x(4,-i,x(-i),smin,smax,s_mass(i),qmass(i)
+     $              ,qwidth(i),cBW_mass(-1,i),cBW_width(-1,i),xjac0
+     $              ,s(i))
             else
 c     normal BW
-               xm02=qmass(i)**2
-               bwmdpl=smax-xm02
-               bwmdmn=xm02-smin
-               bwfmpl=atan(bwmdpl/(qmass(i)*qwidth(i)))
-               bwfmmn=atan(bwmdmn/(qmass(i)*qwidth(i)))
-               bwdelf=(bwfmpl+bwfmmn)/pi
-               s(i)=xbwmass3(x(-i),xm02,qwidth(i),bwdelf
-     &              ,bwfmmn)
-               xjac0=xjac0*bwdelf/bwfunc(s(i),xm02,qwidth(i))
+               call trans_x(3,-i,x(-i),smin,smax,s_mass(i),qmass(i)
+     $              ,qwidth(i),cBW_mass(-1,i),cBW_width(-1,i),xjac0
+     $              ,s(i))
             endif
          else
 c not a Breit Wigner
             if (smin.eq.0d0 .and. s_mass(i).eq.0d0) then
 c     no lower limit on invariant mass from cuts or final state masses:
 c     use flat distribution
-               s(i) = (smax-smin)*x(-i)+smin
-               xjac0 = xjac0*(smax-smin)
+               call trans_x(1,-i,x(-i),smin,smax,s_mass(i),qmass(i)
+     $              ,qwidth(i),cBW_mass(-1,i),cBW_width(-1,i),xjac0
+     $              ,s(i))
             elseif (smin.ge.s_mass(i) .and. smin.gt.0d0) then
 c     A lower limit on smin, which is larger than lower limit from cuts
-c     or masses. Use 1/x^nsamp importance sampling
-               ximax0 = smin**(-nsamp)
-               ximin0 = smax**(-nsamp)
-               tmp  = ximin0 +(1d0-x(-i))*(ximax0-ximin0)
-               s(i) = tmp**(-1/dble(nsamp))
-               xjac0= xjac0/nsamp*s(i)**(nsamp+1)*(ximax0-ximin0)
+c     or masses. Use 1/x importance sampling
+               call trans_x(7,-i,x(-i),smin,smax,s_mass(i),qmass(i)
+     $              ,qwidth(i),cBW_mass(-1,i),cBW_width(-1,i),xjac0
+     $              ,s(i))
             elseif (smin.lt.s_mass(i) .and. s_mass(i).gt.0d0) then
 c     Use flat grid between smin and s_mass(i), and 1/x^nsamp above
 c     s_mass(i)
-               if (x(-i).lt.fract) then
-                  xi=x(-i)/fract ! between 0 and 1
-                  xjac0=xjac0/fract
-                  s(i) = (s_mass(i)-smin)*xi+smin
-                  xjac0 = xjac0*(s_mass(i)-smin)
-               else
-                  xi=(x(-i)-fract)/(1d0-fract) ! between 0 and 1
-                  xjac0=xjac0/(1d0-fract)
-                  ximax0 = s_mass(i)**(-nsamp)
-                  ximin0 = smax**(-nsamp)
-                  tmp  = ximin0 +(1d0-xi)*(ximax0-ximin0)
-                  s(i) = tmp**(-1/dble(nsamp))
-                  xjac0= xjac0/nsamp*s(i)**(nsamp+1)*(ximax0-ximin0)
-               endif
+               call trans_x(2,-i,x(-i),smin,smax,s_mass(i),qmass(i)
+     $              ,qwidth(i),cBW_mass(-1,i),cBW_width(-1,i),xjac0
+     $              ,s(i))
             else
                write (*,*) "ERROR in genps_fks.f:"/
-     $              /" cannot set s-channel without BW"
+     $              /" cannot set s-channel without BW",i,smin,s_mass(i)
                stop 1
             endif
          endif
 c If numerical inaccuracy, quit loop
-         if (xjac0 .lt. 0d0) then
+         if (xjac0.le.0d0) then
+            if ((xjac0.gt.-400d0 .or. xjac0.le.-500d0) .and.
+     $           xjac0.ne.0d0)then
+               write (*,*) 'WARNING #31 in genps_fks.f',i,s(i),smin,smax
+     $              ,xjac0
+            endif
             xjac0 = -6
             pass=.false.
             return
          endif
-         if (s(i) .lt. smin) then
+         if (s(i).lt.smin .or. s(i).gt.smax) then
+            write (*,*) 'WARNING #32 in genps_fks.f',i,s(i),smin,smax,x(
+     $           -i)
             xjac0=-5
             pass=.false.
             return
@@ -2170,6 +2003,8 @@
          totalmass=totalmass+m(i)-
      &        m(itree(1,i))-m(itree(2,i))
          if ( totalmass.gt.sqrtshat_born )then
+            write (*,*) 'WARNING #33 in genps_fks.f',i,totalmass
+     $           ,sqrtshat_born,s(i)
             xjac0 = -4
             pass=.false.
             return
@@ -2180,7 +2015,7 @@
 
 
       subroutine generate_t_channel_branchings(ns_channel,nbranch,itree
-     &     ,m,s,x,pb,xjac0,xpswgt0,pass)
+     $     ,m,s,x,pb,xjac0,xpswgt0,pass)
 c First we need to determine the energy of the remaining particles this
 c is essentially in place of the cos(theta) degree of freedom we have
 c with the s channel decay sequence
@@ -2198,15 +2033,20 @@
       logical pass
 c
       double precision totalmass,smin,smax,s1,ma2,mbq,m12,mnq,tmin,tmax
-     &     ,t,tmax_temp,phi
-      integer i,ibranch
+     &     ,t,tmax_temp,phi,dum,dum3(-1:1),s_m,tm,tiny
+      parameter (tiny=1d-8)
+      integer i,ibranch,idim
       double precision lambda,dot
       external lambda,dot
+      double precision s_mass(-nexternal:nexternal)
+      common/to_phase_space_s_channel/s_mass
 c 
       pass=.true.
       totalmass=0d0
+      s_m=0d0
       do ibranch = -ns_channel-1,-nbranch,-1
          totalmass=totalmass+m(itree(2,ibranch))
+         s_m=s_m+sqrt(s_mass(itree(2,ibranch)))
       enddo
       m(-ns_channel-1) = dsqrt(S(-nbranch))
 c     
@@ -2215,7 +2055,7 @@
 c t-channel propagator down to the initial-state particle found at the end
 c of the t-channel line.
       do ibranch = -ns_channel-1,-nbranch+2,-1
-         totalmass=totalmass-m(itree(2,ibranch))  
+         totalmass=totalmass-m(itree(2,ibranch))
          smin = totalmass**2                    
          smax = (m(ibranch) - m(itree(2,ibranch)))**2
          if (smin .gt. smax) then
@@ -2223,9 +2063,26 @@
             pass=.false.
             return
          endif
-         m(ibranch-1)=dsqrt((smax-smin)*
-     &        x(nbranch-1+(-ibranch)*2)+smin)
-         xjac0 = xjac0*(smax-smin)
+         idim=(nbranch-1+(-ibranch)*2)
+         s_m=s_m-sqrt(s_mass(itree(2,ibranch)))
+         if (abs(smin-s_m**2).lt.tiny) then
+            call trans_x(1,idim,x(idim),smin,smax,s_m**2,dum
+     $           ,dum,dum3(-1),dum3(-1),xjac0,s1)
+         else
+            call trans_x(1,idim,x(idim),smin,smax,s_m**2,dum
+     $           ,dum,dum3(-1),dum3(-1),xjac0,s1)
+         endif
+         if (xjac0.le.0d0) then
+            if ((xjac0.gt.-400d0 .or. xjac0.le.-500d0) .and.
+     $           xjac0.ne.0d0)then
+               write (*,*) 'WARNING #31a in genps_fks.f',ibranch,s1
+     $              ,smin,smax,s_m**2,xjac0
+            endif
+            xjac0 = -6
+            pass=.false.
+            return
+         endif
+         m(ibranch-1)=sqrt(s1)
          if (m(ibranch-1)**2.lt.smin.or.m(ibranch-1)**2.gt.smax
      &        .or.m(ibranch-1).ne.m(ibranch-1)) then
             xjac0=-1d0
@@ -2256,10 +2113,21 @@
          m12 = m(itree(2,ibranch))**2
          mnq = m(ibranch-1)**2
          call yminmax(s1,t,m12,ma2,mbq,mnq,tmin,tmax)
-         tmax_temp = tmax
-         t = (tmax_temp-tmin)*x(-ibranch)+tmin
-         xjac0=xjac0*(tmax_temp-tmin)
+         call trans_x(1,-ibranch,x(-ibranch),-tmax,-tmin,s_mass(ibranch)
+     $        ,dum,dum,dum3(-1),dum3(-1),xjac0,tm)
+         if (xjac0.le.0d0) then
+            if ((xjac0.gt.-400d0 .or. xjac0.le.-500d0) .and.
+     $           xjac0.ne.0d0)then
+               write (*,*) 'WARNING #31b in genps_fks.f',ibranch,tm
+     $              ,-tmax,-tmin,xjac0
+            endif
+            xjac0 = -6
+            pass=.false.
+            return
+         endif
+         t=-tm
          if (t .lt. tmin .or. t .gt. tmax) then
+            write (*,*) "WARNING #35 in genps_fks.f",t,tmin,tmax
             xjac0=-3d0
             pass=.false.
             return
@@ -2395,3 +2263,218 @@
       endif
       return
       end
+
+
+      
+      subroutine trans_x(itype,idim,x,smin,smax,s_mass,qmass,qwidth
+     $     ,cBW_mass,cBW_width,jac,s)
+c Given the input random number 'x', returns the corresponding value of
+c the invariant mass squared 's'. 
+c
+c     itype=1: flat transformation
+c     itype=2: flat between 0 and s_mass/stot, 1/x above
+c     itype=3: Breit-Wigner
+c     itype=4: Conflicting BW, with alternative mass smaller
+c     itype=5: Conflicting BW, with alternative mass larger
+c     itype=6: Conflicting BW on both sides
+c
+      implicit none
+      include 'nexternal.inc'
+      include 'run.inc'         ! ebeam
+c ARGUMENTS
+      integer itype,idim
+      double precision x,smin,smax,s_mass,qmass,qwidth,cBW_mass(-1:1)
+     $     ,cBW_width(-1:1),jac,s
+c LOCAL      
+      double precision tmp,vol_new,x_min,x_max,x_new,fract,x_mass,xg,A
+     $     ,B,C,D,E,F,G,bs(-1:1),maxi,mini
+      integer j
+c
+      if (itype.eq.1) then
+c     flat transformation:
+         A=smax-smin
+         B=smin
+         s=A*x+B
+         jac=jac*A
+      elseif (itype.eq.2) then
+         fract=0.25d0
+         if (s_mass.eq.0d0) then
+            write (*,*) 's_mass is zero',itype,idim
+         endif
+         if (x.lt.fract) then
+c     flat transformation:
+            if (s_mass.lt.smin) then
+               jac=-421d0
+               return
+            endif
+            maxi=min(s_mass,smax)
+            A=(maxi-smin)/fract
+            B=smin
+            s=A*x+B
+            jac=jac*A
+         else
+c     S=A/(B-x) transformation:
+            if (s_mass.gt.smax) then
+               jac=-422d0
+               return
+            endif
+            mini=max(s_mass,smin)
+            A=mini*smax*(1d0-fract)/(smax-mini)
+            B=(smax-fract*mini)/(smax-mini)
+            s=A/(B-x)
+            jac=jac*s**2/A
+         endif
+      elseif(itype.eq.3) then
+c     Normal Breit-Wigner, i.e.
+c        \int_smin^smax ds g(s)/((s-qmass^2)^2-qmass^2*qwidth^2) =
+c        \int_0^1 dx g(s(x))
+         A=atan((qmass-smin/qmass)/qwidth)
+         B=atan((qmass-smax/qmass)/qwidth)
+         s=qmass*(qmass-qwidth*tan(A-(A-B)*x))
+         jac=jac*qmass*qwidth*(A-B)/(cos(A-(A-B)*x))**2
+      elseif(itype.eq.4) then
+c     Conflicting BW, with alternative mass smaller than current
+c     mass. That is, we need to throw also many events at smaller masses
+c     than the peak of the current BW. Split 'x' at 'bs(-1)', using a
+c     flat distribution below the split, and a BW above the split.
+         fract=0.3d0
+         bs(-1)=(cBW_mass(-1)-qmass)/
+     &        (qwidth+cBW_width(-1)) ! bs(-1) is negative here
+         bs(-1)=qmass+bs(-1)*qwidth
+         bs(-1)=bs(-1)**2
+         if (x.lt.fract) then
+            if(smin.gt.bs(-1)) then
+               jac=-441d0
+               return
+            endif
+            maxi=min(bs(-1),smax)
+            A=(maxi-smin)/fract
+            B=smin
+            s=A*x+B
+            jac=jac*A
+         else
+            if(smax.lt.bs(-1)) then
+               jac=-442d0
+               return
+            endif
+            mini=max(bs(-1),smin)
+            A=atan((qmass-mini/qmass)/qwidth)
+            B=atan((qmass-smax/qmass)/qwidth)
+            C=((1d0-x)*A+(x-fract)*B)/(1d0-fract)
+            s=qmass*(qmass-qwidth*tan(C))
+            jac=jac*qmass*qwidth*(A-B)/((cos(C))**2*(1d0-fract))
+         endif
+      elseif(itype.eq.5) then
+c     Conflicting BW, with alternative mass larger than current
+c     mass. That is, we need to throw also many events at larger masses
+c     than the peak of the current BW. Split 'x' at 'bs(1)' and the
+c     alternative mass. Use a BW below bs(1), a flat distribution
+c     between bs(1) and the alternative mass, and a 1/x above the
+c     alternative mass.
+         fract=0.35d0
+         bs(1)=(cBW_mass(1)-qmass)/
+     &        (qwidth+cBW_width(1))
+         bs(1)=qmass+bs(1)*qwidth
+         bs(1)=bs(1)**2
+         if (x.lt.fract) then
+            if(smin.gt.bs(1)) then
+               jac=-451d0
+               return
+            endif
+            maxi=min(bs(1),smax)
+            A=atan((qmass-smin/qmass)/qwidth)
+            B=atan((qmass-maxi/qmass)/qwidth)
+            C=((B-A)*x+fract*A)/fract
+            s=qmass*(qmass-qwidth*tan(C))
+            jac=jac*qmass*qwidth*(A-B)/((cos(C))**2*fract)
+         elseif (x.lt.1d0-fract) then
+            if(smin.gt.cBW_mass(1)**2 .or. smax.lt.bs(1)) then
+               jac=-452d0
+               return
+            endif
+            maxi=min(cBW_mass(1)**2,smax)
+            mini=max(bs(1),smin)
+            A=(maxi-mini)/(1d0-2d0*fract)
+            B=((1d0-fract)*mini-fract*maxi)/(1d0-2d0*fract)
+            s=A*x+B
+            jac=jac*A
+         else
+            if(smax.lt.cBW_mass(1)**2) then
+               jac=-453d0
+               return
+            endif
+            mini=max(cBW_mass(1)**2,smin)
+            A=mini*smax*fract/(smax-mini)
+            B=(smax-(1d0-fract)*mini)/(smax-mini)
+            s=A/(B-x)
+            jac=jac*s**2/A
+         endif
+      elseif(itype.eq.6) then
+         fract=0.3d0
+c     Conflicting BW on both sides. Use flat below bs(-1); BW between
+c     bs(-1) and bs(1); flat between bs(1) and alternative mass; and 1/x
+c     above alternative mass.
+         do j=-1,1,2
+            bs(j)=(cBW_mass(j)-qmass)/
+     &           (qwidth+cBW_width(j))
+            bs(j)=qmass+bs(j)*qwidth
+            bs(j)=bs(j)**2
+         enddo
+         if (x.lt.fract) then
+            if(smin.gt.bs(-1)) then
+               jac=-461d0
+               return
+            endif
+            maxi=min(bs(-1),smax)
+            A=(maxi-smin)/fract
+            B=smin
+            s=A*x+B
+            jac=jac*A
+         elseif(x.lt.1d0-fract) then
+            if(smin.gt.bs(1) .or. smax.lt.bs(-1)) then
+               jac=-462d0
+               return
+            endif
+            maxi=min(bs(1),smax)
+            mini=max(bs(-1),smin)
+            A=atan((qmass-mini/qmass)/qwidth)
+            B=atan((qmass-maxi/qmass)/qwidth)
+            C=((1d0-fract-x)*A+(x-fract)*B)/(1d0-2d0*fract)
+            s=qmass*(qmass-qwidth*tan(C))
+            jac=-jac*qmass*qwidth*(B-A)/((cos(C))**2*(1d0-2d0*fract))
+         elseif(x.lt.1d0-fract/2d0) then
+            if(smin.gt.cBW_mass(1)**2 .or. smax.lt.bs(1)) then
+               jac=-463d0
+               return
+            endif
+            maxi=min(cBW_mass(1)**2,smax)
+            mini=max(bs(1),smin)
+            A=2d0*(maxi-mini)/fract
+            B=2d0*maxi-mini-2d0*(maxi-mini)/fract
+            s=A*x+B
+            jac=jac*A
+         else
+            if(smax.lt.cBW_mass(1)**2) then
+               jac=-464d0
+               return
+            endif
+            mini=max(cBW_mass(1)**2,smin)
+            A=mini*smax*fract/(2d0*(smax-mini))
+            B=(smax-(1d0-fract/2d0)*mini)/(smax-mini)
+            s=A/(B-x)
+            jac=jac*s**2/A
+         endif
+      elseif (itype.eq.7) then
+c     S=A/(B-x) transformation:
+         if (smin.le.0d0) then
+            jac=-471d0
+            return
+         endif
+         A=smin*smax/(smax-smin)
+         B=smax/(smax-smin)
+         s=A/(B-x)
+         jac=jac*s**2/A
+      endif
+      return
+      end
+

=== modified file 'Template/NLO/SubProcesses/mint-integrator2.f'
--- Template/NLO/SubProcesses/mint-integrator2.f	2015-11-20 15:07:40 +0000
+++ Template/NLO/SubProcesses/mint-integrator2.f	2016-09-30 09:38:53 +0000
@@ -96,6 +96,8 @@
       double precision average_virtual,virtual_fraction
       common/c_avg_virt/average_virtual,virtual_fraction
       character*13 title(nintegrals)
+      logical new_point
+      common /c_new_point/ new_point
       data title(1)/'ABS integral '/
       data title(2)/'Integral     '/
       data title(3)/'Virtual      '/
@@ -237,6 +239,7 @@
 c Loop over PS points
  2    kpoint_iter=kpoint_iter+1
       do kpoint=1,ncalls
+         new_point=.true.
 c find random x, and its random cell
          do kdim=1,ndim
             kfold(kdim)=1
@@ -258,7 +261,7 @@
          ifirst=0
  1       continue
          vol=1
-c compute jacobian ('vol') for the PS point
+c c convert 'flat x' ('rand') to 'vegas x' ('x') and include jacobian ('vol')
          do kdim=1,ndim
             nintcurr=nint_used/ifold(kdim)
             icell(kdim)=ncell(kdim)+(kfold(kdim)-1)*nintcurr
@@ -360,7 +363,8 @@
       endif
 c Goto beginning of loop over PS points until enough points have found
 c that pass cuts.
-      if (non_zero_point(1).lt.ncalls .and. double_events) goto 2
+      if (non_zero_point(1).lt.int(0.99*ncalls)
+     &                        .and. double_events) goto 2
 
 c Iteration done. Update the accumulated results and print them to the
 c screen
@@ -652,36 +656,54 @@
       include "mint.inc"
       integer  ninter,nhits(nintervals)
       real * 8 xacc(0:nintervals),xgrid(0:nintervals)
-      real * 8 xn(nintervals),r,tiny,xl,xu,nl,nu
+      real * 8 xn(nintervals),r,tiny,xl,xu,nl,nu,sum
       parameter ( tiny=1d-8 )
       integer kint,jint
-c Use the same smoothing as in VEGAS uses for the grids (i.e. use the
-c average of the central and the two neighbouring grid points):
-      xl=xacc(1)
-      xu=xacc(2)
-      xacc(1)=(xl+xu)/2d0
-      nl=nhits(1)
-      nu=nhits(2)
-      nhits(1)=nint((nl+nu)/2d0)
-      do kint=2,ninter-1
-         xacc(kint)=xl+xu
-         xl=xu
-         xu=xacc(kint+1)
-         xacc(kint)=(xacc(kint)+xu)/3d0
-         nhits(kint)=nl+nu
-         nl=nu
-         nu=nhits(kint+1)
-         nhits(kint)=nint((nhits(kint)+nu)/3d0)
+c Use the same smoothing as in VEGAS uses for the grids, i.e. use the
+c average of the central and the two neighbouring grid points: (Only do
+c this if we are already at the maximum intervals, because the doubling
+c of the grids also includes a smoothing).
+      if (ninter.eq.nintervals) then
+         xl=xacc(1)
+         xu=xacc(2)
+         xacc(1)=(xl+xu)/2d0
+         nl=nhits(1)
+         nu=nhits(2)
+         nhits(1)=nint((nl+nu)/2d0)
+         do kint=2,ninter-1
+            xacc(kint)=xl+xu
+            xl=xu
+            xu=xacc(kint+1)
+            xacc(kint)=(xacc(kint)+xu)/3d0
+            nhits(kint)=nl+nu
+            nl=nu
+            nu=nhits(kint+1)
+            nhits(kint)=nint((nhits(kint)+nu)/3d0)
+         enddo
+         xacc(ninter)=(xu+xl)/2d0
+         nhits(ninter)=nint((nu+nl)/2d0)
+      endif
+c
+      sum=0d0
+      do kint=1,ninter
+         if (nhits(kint).ne.0) then
+            xacc(kint)=abs(xacc(kint))/nhits(kint)
+            sum=sum+xacc(kint)
+         endif
       enddo
-      xacc(ninter)=(xu+xl)/2d0
-      nhits(ninter)=nint((nu+nl)/2d0)
 c
       do kint=1,ninter
 c xacc (xerr) already contains a factor equal to the interval size
 c Thus the integral of rho is performed by summing up
          if(nhits(kint).ne.0) then
-            xacc(kint)= xacc(kint-1)
-     #           + abs(xacc(kint))/nhits(kint)
+c     take logarithm to help convergence (taken from LO dsample.f)
+            if (xacc(kint).ne.sum) then
+               xacc(kint)=((xacc(kint)/sum-1d0)/
+     &                        log(xacc(kint)/sum))**1.5
+            else
+               xacc(kint)=1d0
+            endif
+            xacc(kint)= xacc(kint-1) + abs(xacc(kint))
          else
             xacc(kint)=xacc(kint-1)
          endif
@@ -791,6 +813,8 @@
       common/c_avg_virt/average_virtual,virtual_fraction
       save icalls,mcalls,icalls_virt,mcalls_virt,xmmm,icalls_nz
      $     ,icalls_virt_nz
+      logical new_point
+      common /c_new_point/ new_point
       if(imode.eq.0) then
          do kdim=1,ndim
             nintcurr=nintervals/ifold(kdim)
@@ -844,6 +868,7 @@
          stop
       endif
  10   continue
+      new_point=.true.
       if (vn.eq.1) then
          icalls_virt=icalls_virt+1
       elseif(vn.eq.2 .or. vn.eq.3) then

=== modified file 'Template/NLO/SubProcesses/mint.inc'
--- Template/NLO/SubProcesses/mint.inc	2013-11-22 15:49:06 +0000
+++ Template/NLO/SubProcesses/mint.inc	2016-09-30 09:38:53 +0000
@@ -4,11 +4,11 @@
 c where 'n' is smaller or equal to 'min_it'.
 
       integer nintervals,ndimmax,nintegrals,nintervals_virt
-      parameter (nintervals=32,ndimmax=60,nintegrals=6,
+      parameter (nintervals=64,ndimmax=60,nintegrals=6,
      &		nintervals_virt=8)
 
       integer min_inter,min_it0,min_it1,max_it,max_points
-      data min_inter /4/
+      data min_inter /8/
       data min_it0 /4/
       data min_it1 /5/
       data max_it /100/

=== modified file 'Template/NLO/SubProcesses/setcuts.f'
--- Template/NLO/SubProcesses/setcuts.f	2015-12-16 16:13:10 +0000
+++ Template/NLO/SubProcesses/setcuts.f	2016-09-30 09:38:53 +0000
@@ -156,9 +156,9 @@
       INTEGER NFKSPROCESS
       COMMON/C_NFKSPROCESS/NFKSPROCESS
       double precision taumin(fks_configs),taumin_s(fks_configs)
-     &     ,taumin_j(fks_configs),stot
+     &     ,taumin_j(fks_configs),stot,xk(-nexternal:0)
       save  taumin,taumin_s,taumin_j,stot
-      integer i,j,k,d1,d2,iFKS
+      integer i,j,k,d1,d2,iFKS,nt
       double precision xm(-nexternal:nexternal),xm1,xm2,xmi
       double precision xw(-nexternal:nexternal),xw1,xw2,xwi
       integer tsign,j_fks
@@ -173,16 +173,16 @@
       integer cBW_FKS_level_max(fks_configs),
      &     cBW_FKS(fks_configs,-nexternal:-1),
      &     cBW_FKS_level(fks_configs,-nexternal:-1)
-      double precision cBW_FKS_mass(fks_configs,-nexternal:-1,-1:1),
-     &     cBW_FKS_width(fks_configs,-nexternal:-1,-1:1)
+      double precision cBW_FKS_mass(fks_configs,-1:1,-nexternal:-1),
+     &     cBW_FKS_width(fks_configs,-1:1,-nexternal:-1)
       save cBW_FKS_level_max,cBW_FKS,cBW_FKS_level,cBW_FKS_mass
      $     ,cBW_FKS_width
       integer cBW_level_max,cBW(-nexternal:-1),cBW_level(-nexternal:-1)
-      double precision cBW_mass(-nexternal:-1,-1:1),
-     &     cBW_width(-nexternal:-1,-1:1)
+      double precision cBW_mass(-1:1,-nexternal:-1),
+     &     cBW_width(-1:1,-nexternal:-1)
       common/c_conflictingBW/cBW_mass,cBW_width,cBW_level_max,cBW
      $     ,cBW_level
-      double precision s_mass(-nexternal:-1)
+      double precision s_mass(-nexternal:nexternal)
      $     ,s_mass_FKS(fks_configs,-nexternal:nexternal)
       save s_mass_FKS
       common/to_phase_space_s_channel/s_mass
@@ -359,9 +359,10 @@
 c Also find the minimum lower bound if all internal s-channel particles
 c were on-shell
             tsign=-1
+            nt=0
             do i=-1,-(nexternal-3),-1 ! All propagators
                if ( itree(1,i) .eq. 1 .or. itree(1,i) .eq. 2 ) tsign=1
-               if (tsign.eq.-1) then ! Only s-channels
+               if (tsign.eq.-1) then ! s-channels
                   d1=itree(1,i)
                   d2=itree(2,i)
 c If daughter is a jet, we should treat the ptj as a mass. Except if
@@ -388,22 +389,31 @@
 c subtract the daughters, because they are already included above or in
 c the previous iteration of the loop
                   taumin_s(iFKS)=taumin_s(iFKS)+xm(i)-xm1-xm2
-               else
-                  xm(i)=0d0
+               else             ! t-channels
+                  if (i.eq.-(nexternal-3)) then
+                     xm(i)=0d0
+                     cycle
+                  endif
+                  nt=nt+1
+                  d1=itree(2,i) ! only use 2nd daughter (which is the outgoing one)
+                  xm1=xm(d1)
+                  if (nt.gt.1) xm1=max(xm1,xk(nt-1))
+                  xk(nt)=xm1
+                  j=i-1         ! this is the closest to p2
+                  d2=itree(2,j)
+                  xm2=xm(d2)
+                  xm(i)=min(xm1,xm2)
                endif
             enddo
-
+      
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c Determine the "minimal" s-channel invariant masses
             do i=nincoming+1,nexternal-1
                s_mass_FKS(iFKS,i)=xm(i)**2
             enddo
-            do i=-1,-(nexternal-3),-1 ! All propagators
-               if ( itree(1,i) .eq. 1 .or. itree(1,i) .eq. 2 ) exit ! only s-channels
-               s_mass_FKS(iFKS,i)=(sqrt(s_mass_FKS(iFKS,itree(1,i)))
-     $              +sqrt(s_mass_FKS(iFKS,itree(2,i))))**2
+            do i=-1,-(nexternal-3),-1 ! All propagators 
+               s_mass_FKS(iFKS,i)=xm(i)**2
             enddo
-
 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
 c Determine the conflicting Breit-Wigner's. Note that xm(i) contains the
 c mass of the BW
@@ -413,10 +423,10 @@
             cBW_FKS_level_max(iFKS)=0
             t_channel=0
             do i=-1,-(nexternal-3),-1 ! All propagators
-               cBW_FKS_mass(iFKS,i,1)=0d0
-               cBW_FKS_width(iFKS,i,1)=0d0
-               cBW_FKS_mass(iFKS,i,-1)=0d0
-               cBW_FKS_width(iFKS,i,-1)=0d0
+               cBW_FKS_mass(iFKS,1,i)=0d0
+               cBW_FKS_width(iFKS,1,i)=0d0
+               cBW_FKS_mass(iFKS,-1,i)=0d0
+               cBW_FKS_width(iFKS,-1,i)=0d0
                masslow(i)=9d99
                widthlow(i)=0d0
                if ( itree(1,i).eq.1 .or. itree(1,i).eq.2 ) t_channel=i
@@ -442,10 +452,10 @@
                      cBW_FKS_level_max(iFKS)=max(cBW_FKS_level_max(iFKS)
      $                    ,cBW_FKS_level(iFKS,i))
 c     Set here the mass (and width) of the alternative mass; it's the
-c     sum of daughter masses. (3rd argument is '1', because this
+c     sum of daughter masses. (2nd argument is '1', because this
 c     alternative mass is LARGER than the resonance mass).
-                     cBW_FKS_mass(iFKS,i,1)=xm(i)
-                     cBW_FKS_width(iFKS,i,1)=xw(i)
+                     cBW_FKS_mass(iFKS,1,i)=xm(i)
+                     cBW_FKS_width(iFKS,1,i)=xw(i)
                   endif
 c     set the daughters also as conflicting (recursively)
                   masslow(i)=pmass(i,iconfig)
@@ -472,9 +482,9 @@
                         if (pwidth(itree(k,j),iconfig).eq.0d0 .or.
      $                       masslow(itree(k,j)).ge.pmass(itree(k,j)
      $                       ,iconfig)) cycle
-                        cBW_FKS_mass(iFKS,itree(k,j),-1)=
+                        cBW_FKS_mass(iFKS,-1,itree(k,j))=
      $                       masslow(itree(k,j))
-                        cBW_FKS_width(iFKS,itree(k,j),-1)=
+                        cBW_FKS_width(iFKS,-1,itree(k,j))=
      $                       widthlow(itree(k,j))
                      enddo
                   enddo
@@ -503,8 +513,8 @@
                      if (itree(2,i).gt.0) cycle
                      if (cBW_FKS(iFKS,itree(2,i)).ne.2) then
                         cBW_FKS(iFKS,itree(2,i))=1
-                        cBW_FKS_mass(iFKS,itree(2,i),-1)=sqrt(stot)/2d0
-                        cBW_FKS_width(iFKS,itree(2,i),-1)=xw(itree(2,i))
+                        cBW_FKS_mass(iFKS,-1,itree(2,i))=sqrt(stot)/2d0
+                        cBW_FKS_width(iFKS,-1,itree(2,i))=xw(itree(2,i))
                      endif
                   enddo
                endif
@@ -543,14 +553,71 @@
          cBW(i)=cBW_FKS(nFKSprocess,i)
          cBW_level(i)=cBW_FKS_level(nFKSprocess,i)
          do j=-1,1,2
-            cBW_mass(i,j)=cBW_FKS_mass(nFKSprocess,i,j)
-            cBW_width(i,j)=cBW_FKS_width(nFKSprocess,i,j)
+            cBW_mass(j,i)=cBW_FKS_mass(nFKSprocess,j,i)
+            cBW_width(j,i)=cBW_FKS_width(nFKSprocess,j,i)
          enddo
+      enddo
+      do i=-nexternal,nexternal
          s_mass(i)=s_mass_FKS(nFKSprocess,i)
       enddo
       cBW_level_max=cBW_FKS_level_max(nFKSprocess)
-         
-      return
-      end
-
-
+      return
+      end
+
+
+      subroutine sChan_order(ns_channel,order)
+      implicit none
+      include 'nexternal.inc'
+      include 'maxparticles.inc'
+      include 'maxconfigs.inc'
+      include 'mint.inc'
+      double precision pmass(-nexternal:0,lmaxconfigs)
+      double precision pwidth(-nexternal:0,lmaxconfigs)
+      integer pow(-nexternal:0,lmaxconfigs)
+      integer itree(2,-max_branch:-1),iconfig
+      common /to_itree/itree,iconfig
+      logical new_point
+      common /c_new_point/new_point
+      double precision ran2,rnd
+      integer i,j,order(-nexternal:0),ipos,ns_channel,npos
+     $     ,pos(nexternal),ord(-nexternal:0)
+      logical done(-nexternal:nexternal)
+      external ran2
+      save ord
+      if (.not. new_point) then
+         do j=-1,-ns_channel,-1
+            order(j)=ord(j)
+         enddo
+         return
+      endif
+      do i=-ns_channel,0
+         done(i)=.false.
+      enddo
+      do i=1,nexternal
+         done(i)=.true.
+      enddo
+      do j=-1,-ns_channel,-1
+         npos=0
+         do i=-1,-ns_channel,-1
+            if((.not. done(i)) .and.
+     &           done(itree(1,i))  .and. done(itree(2,i))) then
+               npos=npos+1
+               pos(npos)=i
+            endif
+         enddo
+         if (npos.gt.1) then
+            rnd=ran2()
+            ipos=min(int(rnd*npos)+1,npos)
+            ord(j)=pos(ipos)
+            done(pos(ipos))=.true.
+         elseif (npos.eq.1) then
+            ord(j)=pos(npos)
+            done(pos(npos))=.true.
+         else
+            write (*,*) 'ERROR in sChan_order',npos
+         endif
+         order(j)=ord(j)
+      enddo
+      new_point=.false.
+      return
+      end

=== modified file 'Template/NLO/SubProcesses/symmetry_fks_test_MC.f'
--- Template/NLO/SubProcesses/symmetry_fks_test_MC.f	2016-06-03 09:49:40 +0000
+++ Template/NLO/SubProcesses/symmetry_fks_test_MC.f	2016-09-30 09:38:53 +0000
@@ -159,7 +159,10 @@
       common /cshowerscale2/shower_S_scale,shower_H_scale,ref_H_scale
      &     ,pt_hardness
 
-c      integer icomp
+      logical new_point
+      common /c_new_point/new_point
+
+c     integer icomp
 c
 c     DATA
 c
@@ -342,6 +345,7 @@
       do jj=1,ndim
          x(jj)=ran2()
       enddo
+      new_point=.true.
       call generate_momenta(ndim,iconfig,wgt,x,p)
       calculatedBorn=.false.
       do while (( wgt.lt.0 .or. p(0,1).le.0d0 .or. p_born(0,1).le.0d0
@@ -349,6 +353,7 @@
          do jj=1,ndim
             x(jj)=ran2()
          enddo
+         new_point=.true.
          wgt=1d0
          call generate_momenta(ndim,iconfig,wgt,x,p)
          calculatedBorn=.false.
@@ -388,12 +393,14 @@
          do jj=1,ndim
             x(jj)=ran2()
          enddo
+         new_point=.true.
          call generate_momenta(ndim,iconfig,wgt,x,p)
          do while (( wgt.lt.0 .or. p(0,1).le.0d0) .and. ntry.lt.1000)
             wgt=1d0
             do jj=1,ndim
                x(jj)=ran2()
             enddo
+            new_point=.true.
             call generate_momenta(ndim,iconfig,wgt,x,p)
             ntry=ntry+1
          enddo
@@ -535,12 +542,14 @@
          do jj=1,ndim
             x(jj)=ran2()
          enddo
+         new_point=.true.
          call generate_momenta(ndim,iconfig,wgt,x,p)
          do while (( wgt.lt.0 .or. p(0,1).le.0d0) .and. ntry.lt.1000)
             wgt=1d0
             do jj=1,ndim
                x(jj)=ran2()
             enddo
+            new_point=.true.
             call generate_momenta(ndim,iconfig,wgt,x,p)
             ntry=ntry+1
          enddo

=== modified file 'Template/NLO/SubProcesses/symmetry_fks_test_ME.f'
--- Template/NLO/SubProcesses/symmetry_fks_test_ME.f	2016-06-03 09:49:40 +0000
+++ Template/NLO/SubProcesses/symmetry_fks_test_ME.f	2016-09-30 09:38:53 +0000
@@ -139,6 +139,8 @@
       LOGICAL  IS_A_PH(NEXTERNAL)
       COMMON /TO_SPECISA/IS_A_J,IS_A_LP,IS_A_LM,IS_A_PH
       
+      logical new_point
+      common /c_new_point/new_point
 c      integer icomp
 c
 c     DATA
@@ -289,6 +291,7 @@
       do jj=1,ndim
          x(jj)=ran2()
       enddo
+      new_point=.true.
       call generate_momenta(ndim,iconfig,wgt,x,p)
       calculatedBorn=.false.
       do while (( wgt.lt.0 .or. p(0,1).le.0d0 .or. p_born(0,1).le.0d0
@@ -296,6 +299,7 @@
          do jj=1,ndim
             x(jj)=ran2()
          enddo
+         new_point=.true.
          wgt=1d0
          call generate_momenta(ndim,iconfig,wgt,x,p)
          calculatedBorn=.false.
@@ -334,12 +338,14 @@
          do jj=1,ndim
             x(jj)=ran2()
          enddo
+         new_point=.true.
          call generate_momenta(ndim,iconfig,wgt,x,p)
          do while (( wgt.lt.0 .or. p(0,1).le.0d0) .and. ntry.lt.1000)
             wgt=1d0
             do jj=1,ndim
                x(jj)=ran2()
             enddo
+            new_point=.true.
             call generate_momenta(ndim,iconfig,wgt,x,p)
             ntry=ntry+1
          enddo
@@ -460,12 +466,14 @@
          do jj=1,ndim
             x(jj)=ran2()
          enddo
+         new_point=.true.
          call generate_momenta(ndim,iconfig,wgt,x,p)
          do while (( wgt.lt.0 .or. p(0,1).le.0d0) .and. ntry.lt.1000)
             wgt=1d0
             do jj=1,ndim
                x(jj)=ran2()
             enddo
+            new_point=.true.
             call generate_momenta(ndim,iconfig,wgt,x,p)
             ntry=ntry+1
          enddo

=== modified file 'Template/NLO/SubProcesses/symmetry_fks_v3.f'
--- Template/NLO/SubProcesses/symmetry_fks_v3.f	2016-03-16 14:21:42 +0000
+++ Template/NLO/SubProcesses/symmetry_fks_v3.f	2016-09-30 09:38:53 +0000
@@ -125,6 +125,8 @@
       logical nbody
       common/cnbody/nbody
 
+      logical new_point
+      common /c_new_point/new_point
 c-----
 c  Begin Code
 c-----
@@ -210,12 +212,14 @@
       do j=1,ndim
          x(j)=ran2()
       enddo
+      new_point=.true.
       call generate_momenta(ndim,iconfig,wgt,x,p)
       do while ((.not.passcuts(p,rwgt) .or. wgt.lt.0 .or. p(0,1).le.0d0
      &           .or. p_born(0,1).le.0d0) .and. ntry.lt.10000)
          do j=1,ndim
             x(j)=ran2()
          enddo
+         new_point=.true.
          wgt=1d0
          call generate_momenta(ndim,iconfig,wgt,x,p)
          ntry=ntry+1